Loading test_efficiency.m 0 → 100644 +437 −0 Original line number Diff line number Diff line clear; close all; % Define variables, M is segm. number pathfi = '/data/Sorgenti/3FGLJ1544.6-1125/'; reset(gpuDevice(2)); gpuDevice(2); %%%%%%%%%%%%%%%% tic finame = 'SiFAP2_20220429_3FGLJ1544_N1_Bary_N3_NOTCORRECTED.fits'; % finame = 'J1023_B_2017_Bary.fits'; % finame = 'EPN_0744840201_bary.fits'; t_raw = fitsread([pathfi,finame],"binarytable"); %t_raw = fitsread('C:\Users\Filippo\Desktop\J1544_2022\SiFAP2_20220429_3FGLJ1544_N1_Bary.fits','binarytable'); t_raw = t_raw{1}; MJDREF=0.943055555602768; % t_raw=t_raw./86400+50814; % MJDREF=fix(t_raw(1)); % t_raw=(t_raw-MJDREF).*86400; % t_raw=t_raw(t_raw>=t_raw(1) & t_raw<(t_raw(1)+ 6*256 + 1)); % info = fitsinfo(finame).BinaryTable.Keywords; % for i = 1:length(info) % if (isequal(info{i},'TIMEDEL')) % res = info{i,2}; % end % end %%%%%%%%%%%%%%%% toc %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Lines added just to test what was done until now % Parameters for J1023_B_2017_Bary.fits % f_tru = 592.42146827248556; %Hz % porb_tru = 17115.5216592; %s % a_tru = 0.343356; %lt-s % tasc_tru = (.009477-.2326620370367891)*86400; %MJD? % Parameters for EPN_0744840201_bary.fits % f_tru = 598.8921309; %Hz % porb_tru = 8844.08; %s % a_tru = 0.0649905; %lt-s % % tasc_tru = 57231.437581; %MJD? % tasc_tru = (57231.437581-MJDREF)*86400; %in secondi % f_gr = zeros(5,1); % porb_gr = zeros(5,1); % a_gr = zeros(5,1); % tasc_gr = zeros(5,1); Tseg=490; %segments' length in seconds %f_gr=f_tru+(-2:2).'; f_gr = (768:1/Tseg:780).'; %f_gr=(772:1/Tseg:774).'; % porb_gr = (18000:10:21600).'; porb_gr=((20868.72-3*0.31):0.1:(20868.72+3*0.31)); a_gr = (0.01:8.0e-4:0.26).'; tasc_gr = (((0.724:8e-4:0.922)-MJDREF).*86400).'; % for j = 1:5 % %f_gr(j) = f_tru*((j^2)/9); % porb_gr(j) = porb_tru*((j^2)/9); % a_gr(j) = a_tru*((j^2)/9); % tasc_gr(j) = tasc_tru*((j^2)/9); % end f_min = min(f_gr); f_max = max(f_gr); porb_min = min(porb_gr); porb_max = max(porb_gr); a_min = min(a_gr); a_max = max(a_gr); tasc_min = min(tasc_gr); tasc_max = max(tasc_gr); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %NOTE: t0 used to define the orbital phase γ is equal to the mid point %of the observation span for dataset; M2015, Sec.7 FAAAAAAAAAAAAAAAAALSE % FAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAALSE t0 = (t_raw(end)+t_raw(1))/2; tic %Rebin ------------------------------------------------------------------- Nyq=2000; %cambiare in caso dt_psd=1/(2*Nyq); %risoluzione spettro di potenza, diverso dal dt dei segmenti nbin=round((t_raw(end)-t_raw(1))/dt_psd); [C,t]=(histcounts(t_raw,nbin)); C=C.'; %conteggi t=(t(1:end-1)+(t(2)-t(1))/2).'; % vettore tempi rebinnato, prendo il centro del bin % t=(t(2:end)).'; %vettore tempi rebinnato, assegno al tempo l'edge destro di ogni bin %VEDERE SE E' NECESSARIO FARE TUTTO IL CICLO SOTTO PER I tm %Come fatto qui sopra, tutti i conteggi sono assegnati al centro di ogni %bin che dura dt toc tic M=fix((t(end)-t(1))/Tseg); %number of segments %con fix prendo la parte intera, scarto l'ultimo segmento che tanto non %sarà mai di lunghezza Tseg (molto improbabile) % Time matrix with bin midpoints for each segment tm(m,j) --------------- % tmid(m) is the midpoint in time for the m-th segment % dt = Tseg/N; % dt = dt_psd; N = fix(Tseg/dt); tm = zeros(M,N); tmid = zeros(M,1); for m = 1:M tm(m,1) = t(1)+(m-1)*Tseg; %t(j) è già centrato for j = 2:N tm(m,j)= tm(m,1)+(j-1)*dt; end tmid(m) = (tm(m,N)+tm(m,1))/2; end %CONTROLLARE SE CON LA FUNZIONE RESHAPE PUOI EVITARE IL CICLO FOR % Time series for each segment x = zeros(N,M); toc tic aux1 = M*N; aux2 = C(1:aux1); clear aux1 x = reshape(aux2,[N,M]); clear aux2 toc clear C % % Add step to generate par. extrema param = [f_min,f_max,porb_min,porb_max,a_min,a_max,tasc_min,tasc_max]; % Aux par Omega_max = 2*pi/porb_min; Omega_min = 2*pi/porb_max; % tdiff = tmid - [tasc_max,tasc_min]; % gam = tdiff(:).*[Omega_min,Omega_max]; % % % gamma_min = Omega_min*(t0-tasc_max); % % gamma_max = Omega_max*(t0-tasc_min); % % % Find largest singam % % M2015 eq. 15 % singal = zeros(4,1); % singah = zeros(4,1); nismin = zeros(4,1); nismax = zeros(4,1); % for s = 1:4 % singah(s) = max(sin(gam + s*pi/2),[],'all'); % singal(s) = min(sin(gam + s*pi/2),[],'all'); % end % % %% Tutto sto macello e poi ora mi viene da pensare che valori intermedi di % %% tasc possono portare a sin(gam) anche più alti o più bassi... % %% Guarda come ora avrei fatto molto prima a mettere semplicemente +1 e -1 % % Adding special case for s = 1 % s = 1; % if singal(s)>0 % nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s) + f_max; % nismax(s) = -f_min*a_min*(Omega_min^s)*singal(s) + f_min; % elseif singah(s)>0 % nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s) + f_max; % nismax(s) = -f_max*a_max*(Omega_max^s)*singal(s) + f_max; % else % nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s) + f_max; % nismax(s) = -f_min*a_min*(Omega_min^s)*singal(s) + f_min; % end % for s = 2:4 % % This range is computed by finding the maximum span of Equation (15) after varying the search % % parameters over their respective ranges (given in Table 2). This % % is done with the exception of ν which is held fixed at its % % maximum value within sub-bands over the frequency search space. % %% RECHECK EVERY COMBINATION % if singal(s)>0 % nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s); % nismax(s) = -f_min*a_min*(Omega_min^s)*singal(s); % elseif singah(s)>0 % nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s); % nismax(s) = -f_max*a_max*(Omega_max^s)*singal(s); % else % nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s); % nismax(s) = -f_min*a_min*(Omega_min^s)*singal(s); % end % end % Find the range in ν^s (i.e. the maximum span of Eq. (15) covered by % combining the search parameters over their respective ranges) % Sin(gamma) going from -1 to 1 % Adding special case for s = 1 s = 1; if a_max*Omega_max>1 nismin(s) = -f_max*a_max*(Omega_max) + f_max; if a_min*Omega_min>1 nismax(s) = -f_min*a_min*(Omega_min) + f_min; else nismax(s) = -f_max*a_min*(Omega_min) + f_max; end else nismin(s) = -f_min*a_max*(Omega_max) + f_min; nismax(s) = -f_max*a_min*(Omega_min) + f_max; end % The other s's are easier for s = 2:4 nismin(s) = -f_max*a_max*(Omega_max^s); nismax(s) = f_max*a_max*(Omega_max^s); end toc tic %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %for m = 1:M % for j = 1:N % pippo = t>=tm(m,j)-dt/2 & t<tm(m,j)+dt/2; % Cmj=C(pippo); clear pippo % x(m,j)=sum(Cmj); % end %end % Try s* and check \nu_s range g_jj=(((pi*Tseg)^2)/3).*[1; (Tseg^2)/60; (Tseg^4)/1344; (Tseg^6)/172800]; %eq. 22 M2015 + calcoli da eq. 21 M2015 mu_s=0.05; %massimo mismatch sulla griglia coerente da scegliere %s_s=uint8(4); s_s = 4; while(1) if((nismax(s_s)-nismin(s_s))<0.5*delta_ni(mu_s,s_s,g_jj(s_s))) s_s=s_s-1; else break; end end %% Create a coherent template bank in the ν_s space %% Save this template bank, we'll need it later % Nni+1 is the number of ni_s points in the grid (counting borders, % Nni is the number of steps) for each s couple. % We make sure to oversample with respect to eq. 23 M2015 to ensure the % max mismatch is not exceeded. ------------- % nis is a structure which holds the various vectors of ni_s^(m) (our % coherent template bank) % each nis{s} is the array giving the respective grid, and its i-th % element is accessed as nis{s}(i) Nni = zeros(s_s,1); for s = 1:s_s Nni(s) = ceil((nismax(s)-nismin(s))/delta_ni(mu_s,s,g_jj(s))); if (mod(Nni(s),2)==0) Nni(s) = Nni(s)+1; end franco=nismin(s):((nismax(s)-nismin(s))/(Nni(s))):nismax(s); nis{s}=franco; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Changing nis to stop at s*(= s_s): this way nibank is length(f_gr) % times smaller, which also helps in the NNsearch %%%%%%%%%%%% -> commenting the next line and adding a cycle on f_gr % nis{s_s+1}=f_gr.'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% toc tic % Resampling the times over nism{m,s} templates % ---- Andrà fatto per ogni possibile combinazione di nism a un dato m % mi sa % Cercheremo un modo per combinare in tutti i possibili modi vettori di % vari ni1,ni2,...,nis_s %%%[kung,fu,fight] = ndgrid(nis{1},nino{2},nino{3}); tm=gpuArray(tm); nibank = gpuArray(combinations(nis{:}).Variables); Lambda = zeros(length(nibank),length(f_gr),'gpuArray'); % nibank = combinations(nis{:}).Variables; % Lambda = zeros(length(nibank),length(f_gr)); toc %Fourier transform on original time-series -------------------------------- %per mantenere l'informazione di fase, non faccio il valore assoluto al quadrato della fft %for each segment (lavoro su tm) % dtau = zeros(1,N); % taul = zeros(N+1,1); taul = zeros(2,1); disp(M); for m=1:M % tic % [Cm,edges]=(histcounts(x(m,:),round((tm(m,end)-tm(m,1))/dt_psd))); % R - convincitene % edges=edges(end)-edges(1); %mi dà il tempo preciso di tutta la TdF, che sarà leggermente diversa da length(C)*dt per come è definito histcounts % Y=fft(x(m,:)).'; % F=((0:length(Y)-1)./(tm(m,end)-tm(m,1))).'; %freq. da 0 a 2Nyq. % L=length(F); %lunghezza iniziale, servirà per lo zero-padding % % %Choose freq. region of interest (RoI) and inverse-fourier transf. over RoI % cond = F>=f_min & F<=f_max; % F=F(cond); % Y=Y(cond); % X=ifft(Y); %inverse-fourier transf. % ttemp = gpuArray(tm(m,:)); % xtemp = gpuArray(x(:,m).'); ttemp = tm(m,:); ttempdif = (ttemp(:)-tmid(m)).'; ttempdifl = (ttemp(:)-tmid(m)-0.5*dt).'; xtemp = x(:,m).'; % toc %FARE IL RICAMPIONAMENTO!!!! % Dev'essere fatto per ciascun template, ergo mettiamo un bel ciclo % sulle possibili combinazioni di ni_s tic % for i = 1:length(nibank) for i = 1:1 % Tento ricampionamento (Eq. 17 MP2015) % Per prima cosa, generiamo la nuova coordinata temporale per % questa templ combini % Careful! ni_0 still needs to be defined (it's the spin freq. % currently being considered) for n = 1:length(f_gr) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nizero = f_gr(n); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % tau = zeros(N,1); % tic % for j = 1:N % % ta(j) = ; No Tau_zero: indifferente per il calcolo della FFT % for s = 1:s_s % tau(j) = tau(j) + (nibank(i,s)/(nizero*factorial(s)))*(tm(m,j)-tmid(m))^s; % end % end % tic tau = tmid(m) + sum((nibank(i,1:s_s)./(nizero*factorial(1:s_s))).*((ttempdif(1:N)).^((1:s_s).').'),2); taul(1:2) = tmid(m) + sum((nibank(i,1:s_s)./(nizero*factorial(1:s_s))).*((ttempdifl(1:2)).^((1:s_s).').'),2); dtau = ((taul(2)-taul(1))); Y1 = fft((interp1(tau(:),(xtemp(:)./dtau),ttemp(:),'linear',0)).*dt).'; F1=((0:length(Y1)-1)./(length(Y1)*dt_psd)).'; [nnfreq,nnind] = min(abs(F1-nizero)); Y1=Y1(nnind); clear F1 Lambda(i,n) = 2*(abs(max(Y1)).^2)/sum(xtemp(:)); %CREDO (oppure prendono la potenza massima?) % toc end % disp('1ni') end toc disp(length(f_gr)); disp(m); lamfiname = sprintf('Lambda_seg_%d.mat', m); % save(lamfiname,"Lambda"); save([pathfi,lamfiname],"Lambda"); Lambda = zeros(length(nibank),length(f_gr)); disp(m); end tic % nisearcher = KDTreeSearcher(nibank,'BucketSize',100); nisearcher = KDTreeSearcher(gather(nibank),'BucketSize',100); %% The line above will have to be substituted with something along the lines %% of the ones below to account for the metric g_jj in the phase derivatives % gdistance = @(a,b)sqrt(((a-b).^2)*(g_jj(1:s_s))); % nisearcher = ExhaustiveSearcher(gather(nibank),'Distance',gdistance); % nisearcher = ExhaustiveSearcher(nibank,'Distance',gdistance); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% parbank = combinations(porb_gr,a_gr,tasc_gr).Variables; bestpar = zeros(1,5); toc totlam = zeros(length(parbank),length(f_gr)); tic for n=1:length(f_gr) for i=1:length(parbank) curpar = [f_gr(n),parbank(i,:)]; % Move from ν,P,a,T_asc to ν,Ω,a,γ curpar(2) = 2*pi/curpar(2); for m = 1:M curpar(4) = curpar(2)*(tmid(m) - parbank(i,3)); curni=zeros(1,s_s); for s=1:s_s curni(s) = (curpar(2)^s)*sin(curpar(4)+0.5*s*pi); end curni = -curni.*curpar(1).*curpar(3); curni(1) = curni(1) + curpar(1); lamfiname = sprintf('Lambda_seg_%d.mat', m); % load(lamfiname); load([pathfi,lamfiname]); [Idx,D] = knnsearch(nisearcher,curni); totlam(i,n) = totlam(i,n)+Lambda(Idx,n); clear Lambda end % % if (totlam > bestpar(1)) % % bestpar = [totlam,f_gr(n),parbank(i,:)]; % % end end end toc % disp(bestpar); % save('C:\Users\Filippo\Desktop\XMM_Jxxx\risultelli.mat'); save([pathfi,'risultelli_',num2str(Tseg),'s.mat']); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%% Fai il grafico delle totlam a posteriori %%%%%%%%%%%%%%%% % for n = 1:length(f_gr) % hold on % txt = ['f_{gr} = ',num2str(f_gr(n))]; % plot(totlam(:,n),'DisplayName',txt) % end % legend %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% COMPLETARE DOPO AVER STUDIATO I RANDOM TEMP. BANKS %% E VEDERE CODICE https://pyfstat.readthedocs.io/en/latest/pyfstat.html#pyfstat.core.SemiCoherentSearch %-------------------------------------------------------------------------- %FUNCTIONS %-------------------------------------------------------------------------- function res = delta_ni(mu_s,s,g) %Compute eq. 23 M2015 res= 2.*sqrt(mu_s./(s.*g)); end No newline at end of file Loading
test_efficiency.m 0 → 100644 +437 −0 Original line number Diff line number Diff line clear; close all; % Define variables, M is segm. number pathfi = '/data/Sorgenti/3FGLJ1544.6-1125/'; reset(gpuDevice(2)); gpuDevice(2); %%%%%%%%%%%%%%%% tic finame = 'SiFAP2_20220429_3FGLJ1544_N1_Bary_N3_NOTCORRECTED.fits'; % finame = 'J1023_B_2017_Bary.fits'; % finame = 'EPN_0744840201_bary.fits'; t_raw = fitsread([pathfi,finame],"binarytable"); %t_raw = fitsread('C:\Users\Filippo\Desktop\J1544_2022\SiFAP2_20220429_3FGLJ1544_N1_Bary.fits','binarytable'); t_raw = t_raw{1}; MJDREF=0.943055555602768; % t_raw=t_raw./86400+50814; % MJDREF=fix(t_raw(1)); % t_raw=(t_raw-MJDREF).*86400; % t_raw=t_raw(t_raw>=t_raw(1) & t_raw<(t_raw(1)+ 6*256 + 1)); % info = fitsinfo(finame).BinaryTable.Keywords; % for i = 1:length(info) % if (isequal(info{i},'TIMEDEL')) % res = info{i,2}; % end % end %%%%%%%%%%%%%%%% toc %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Lines added just to test what was done until now % Parameters for J1023_B_2017_Bary.fits % f_tru = 592.42146827248556; %Hz % porb_tru = 17115.5216592; %s % a_tru = 0.343356; %lt-s % tasc_tru = (.009477-.2326620370367891)*86400; %MJD? % Parameters for EPN_0744840201_bary.fits % f_tru = 598.8921309; %Hz % porb_tru = 8844.08; %s % a_tru = 0.0649905; %lt-s % % tasc_tru = 57231.437581; %MJD? % tasc_tru = (57231.437581-MJDREF)*86400; %in secondi % f_gr = zeros(5,1); % porb_gr = zeros(5,1); % a_gr = zeros(5,1); % tasc_gr = zeros(5,1); Tseg=490; %segments' length in seconds %f_gr=f_tru+(-2:2).'; f_gr = (768:1/Tseg:780).'; %f_gr=(772:1/Tseg:774).'; % porb_gr = (18000:10:21600).'; porb_gr=((20868.72-3*0.31):0.1:(20868.72+3*0.31)); a_gr = (0.01:8.0e-4:0.26).'; tasc_gr = (((0.724:8e-4:0.922)-MJDREF).*86400).'; % for j = 1:5 % %f_gr(j) = f_tru*((j^2)/9); % porb_gr(j) = porb_tru*((j^2)/9); % a_gr(j) = a_tru*((j^2)/9); % tasc_gr(j) = tasc_tru*((j^2)/9); % end f_min = min(f_gr); f_max = max(f_gr); porb_min = min(porb_gr); porb_max = max(porb_gr); a_min = min(a_gr); a_max = max(a_gr); tasc_min = min(tasc_gr); tasc_max = max(tasc_gr); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %NOTE: t0 used to define the orbital phase γ is equal to the mid point %of the observation span for dataset; M2015, Sec.7 FAAAAAAAAAAAAAAAAALSE % FAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAALSE t0 = (t_raw(end)+t_raw(1))/2; tic %Rebin ------------------------------------------------------------------- Nyq=2000; %cambiare in caso dt_psd=1/(2*Nyq); %risoluzione spettro di potenza, diverso dal dt dei segmenti nbin=round((t_raw(end)-t_raw(1))/dt_psd); [C,t]=(histcounts(t_raw,nbin)); C=C.'; %conteggi t=(t(1:end-1)+(t(2)-t(1))/2).'; % vettore tempi rebinnato, prendo il centro del bin % t=(t(2:end)).'; %vettore tempi rebinnato, assegno al tempo l'edge destro di ogni bin %VEDERE SE E' NECESSARIO FARE TUTTO IL CICLO SOTTO PER I tm %Come fatto qui sopra, tutti i conteggi sono assegnati al centro di ogni %bin che dura dt toc tic M=fix((t(end)-t(1))/Tseg); %number of segments %con fix prendo la parte intera, scarto l'ultimo segmento che tanto non %sarà mai di lunghezza Tseg (molto improbabile) % Time matrix with bin midpoints for each segment tm(m,j) --------------- % tmid(m) is the midpoint in time for the m-th segment % dt = Tseg/N; % dt = dt_psd; N = fix(Tseg/dt); tm = zeros(M,N); tmid = zeros(M,1); for m = 1:M tm(m,1) = t(1)+(m-1)*Tseg; %t(j) è già centrato for j = 2:N tm(m,j)= tm(m,1)+(j-1)*dt; end tmid(m) = (tm(m,N)+tm(m,1))/2; end %CONTROLLARE SE CON LA FUNZIONE RESHAPE PUOI EVITARE IL CICLO FOR % Time series for each segment x = zeros(N,M); toc tic aux1 = M*N; aux2 = C(1:aux1); clear aux1 x = reshape(aux2,[N,M]); clear aux2 toc clear C % % Add step to generate par. extrema param = [f_min,f_max,porb_min,porb_max,a_min,a_max,tasc_min,tasc_max]; % Aux par Omega_max = 2*pi/porb_min; Omega_min = 2*pi/porb_max; % tdiff = tmid - [tasc_max,tasc_min]; % gam = tdiff(:).*[Omega_min,Omega_max]; % % % gamma_min = Omega_min*(t0-tasc_max); % % gamma_max = Omega_max*(t0-tasc_min); % % % Find largest singam % % M2015 eq. 15 % singal = zeros(4,1); % singah = zeros(4,1); nismin = zeros(4,1); nismax = zeros(4,1); % for s = 1:4 % singah(s) = max(sin(gam + s*pi/2),[],'all'); % singal(s) = min(sin(gam + s*pi/2),[],'all'); % end % % %% Tutto sto macello e poi ora mi viene da pensare che valori intermedi di % %% tasc possono portare a sin(gam) anche più alti o più bassi... % %% Guarda come ora avrei fatto molto prima a mettere semplicemente +1 e -1 % % Adding special case for s = 1 % s = 1; % if singal(s)>0 % nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s) + f_max; % nismax(s) = -f_min*a_min*(Omega_min^s)*singal(s) + f_min; % elseif singah(s)>0 % nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s) + f_max; % nismax(s) = -f_max*a_max*(Omega_max^s)*singal(s) + f_max; % else % nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s) + f_max; % nismax(s) = -f_min*a_min*(Omega_min^s)*singal(s) + f_min; % end % for s = 2:4 % % This range is computed by finding the maximum span of Equation (15) after varying the search % % parameters over their respective ranges (given in Table 2). This % % is done with the exception of ν which is held fixed at its % % maximum value within sub-bands over the frequency search space. % %% RECHECK EVERY COMBINATION % if singal(s)>0 % nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s); % nismax(s) = -f_min*a_min*(Omega_min^s)*singal(s); % elseif singah(s)>0 % nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s); % nismax(s) = -f_max*a_max*(Omega_max^s)*singal(s); % else % nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s); % nismax(s) = -f_min*a_min*(Omega_min^s)*singal(s); % end % end % Find the range in ν^s (i.e. the maximum span of Eq. (15) covered by % combining the search parameters over their respective ranges) % Sin(gamma) going from -1 to 1 % Adding special case for s = 1 s = 1; if a_max*Omega_max>1 nismin(s) = -f_max*a_max*(Omega_max) + f_max; if a_min*Omega_min>1 nismax(s) = -f_min*a_min*(Omega_min) + f_min; else nismax(s) = -f_max*a_min*(Omega_min) + f_max; end else nismin(s) = -f_min*a_max*(Omega_max) + f_min; nismax(s) = -f_max*a_min*(Omega_min) + f_max; end % The other s's are easier for s = 2:4 nismin(s) = -f_max*a_max*(Omega_max^s); nismax(s) = f_max*a_max*(Omega_max^s); end toc tic %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %for m = 1:M % for j = 1:N % pippo = t>=tm(m,j)-dt/2 & t<tm(m,j)+dt/2; % Cmj=C(pippo); clear pippo % x(m,j)=sum(Cmj); % end %end % Try s* and check \nu_s range g_jj=(((pi*Tseg)^2)/3).*[1; (Tseg^2)/60; (Tseg^4)/1344; (Tseg^6)/172800]; %eq. 22 M2015 + calcoli da eq. 21 M2015 mu_s=0.05; %massimo mismatch sulla griglia coerente da scegliere %s_s=uint8(4); s_s = 4; while(1) if((nismax(s_s)-nismin(s_s))<0.5*delta_ni(mu_s,s_s,g_jj(s_s))) s_s=s_s-1; else break; end end %% Create a coherent template bank in the ν_s space %% Save this template bank, we'll need it later % Nni+1 is the number of ni_s points in the grid (counting borders, % Nni is the number of steps) for each s couple. % We make sure to oversample with respect to eq. 23 M2015 to ensure the % max mismatch is not exceeded. ------------- % nis is a structure which holds the various vectors of ni_s^(m) (our % coherent template bank) % each nis{s} is the array giving the respective grid, and its i-th % element is accessed as nis{s}(i) Nni = zeros(s_s,1); for s = 1:s_s Nni(s) = ceil((nismax(s)-nismin(s))/delta_ni(mu_s,s,g_jj(s))); if (mod(Nni(s),2)==0) Nni(s) = Nni(s)+1; end franco=nismin(s):((nismax(s)-nismin(s))/(Nni(s))):nismax(s); nis{s}=franco; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Changing nis to stop at s*(= s_s): this way nibank is length(f_gr) % times smaller, which also helps in the NNsearch %%%%%%%%%%%% -> commenting the next line and adding a cycle on f_gr % nis{s_s+1}=f_gr.'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% toc tic % Resampling the times over nism{m,s} templates % ---- Andrà fatto per ogni possibile combinazione di nism a un dato m % mi sa % Cercheremo un modo per combinare in tutti i possibili modi vettori di % vari ni1,ni2,...,nis_s %%%[kung,fu,fight] = ndgrid(nis{1},nino{2},nino{3}); tm=gpuArray(tm); nibank = gpuArray(combinations(nis{:}).Variables); Lambda = zeros(length(nibank),length(f_gr),'gpuArray'); % nibank = combinations(nis{:}).Variables; % Lambda = zeros(length(nibank),length(f_gr)); toc %Fourier transform on original time-series -------------------------------- %per mantenere l'informazione di fase, non faccio il valore assoluto al quadrato della fft %for each segment (lavoro su tm) % dtau = zeros(1,N); % taul = zeros(N+1,1); taul = zeros(2,1); disp(M); for m=1:M % tic % [Cm,edges]=(histcounts(x(m,:),round((tm(m,end)-tm(m,1))/dt_psd))); % R - convincitene % edges=edges(end)-edges(1); %mi dà il tempo preciso di tutta la TdF, che sarà leggermente diversa da length(C)*dt per come è definito histcounts % Y=fft(x(m,:)).'; % F=((0:length(Y)-1)./(tm(m,end)-tm(m,1))).'; %freq. da 0 a 2Nyq. % L=length(F); %lunghezza iniziale, servirà per lo zero-padding % % %Choose freq. region of interest (RoI) and inverse-fourier transf. over RoI % cond = F>=f_min & F<=f_max; % F=F(cond); % Y=Y(cond); % X=ifft(Y); %inverse-fourier transf. % ttemp = gpuArray(tm(m,:)); % xtemp = gpuArray(x(:,m).'); ttemp = tm(m,:); ttempdif = (ttemp(:)-tmid(m)).'; ttempdifl = (ttemp(:)-tmid(m)-0.5*dt).'; xtemp = x(:,m).'; % toc %FARE IL RICAMPIONAMENTO!!!! % Dev'essere fatto per ciascun template, ergo mettiamo un bel ciclo % sulle possibili combinazioni di ni_s tic % for i = 1:length(nibank) for i = 1:1 % Tento ricampionamento (Eq. 17 MP2015) % Per prima cosa, generiamo la nuova coordinata temporale per % questa templ combini % Careful! ni_0 still needs to be defined (it's the spin freq. % currently being considered) for n = 1:length(f_gr) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nizero = f_gr(n); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % tau = zeros(N,1); % tic % for j = 1:N % % ta(j) = ; No Tau_zero: indifferente per il calcolo della FFT % for s = 1:s_s % tau(j) = tau(j) + (nibank(i,s)/(nizero*factorial(s)))*(tm(m,j)-tmid(m))^s; % end % end % tic tau = tmid(m) + sum((nibank(i,1:s_s)./(nizero*factorial(1:s_s))).*((ttempdif(1:N)).^((1:s_s).').'),2); taul(1:2) = tmid(m) + sum((nibank(i,1:s_s)./(nizero*factorial(1:s_s))).*((ttempdifl(1:2)).^((1:s_s).').'),2); dtau = ((taul(2)-taul(1))); Y1 = fft((interp1(tau(:),(xtemp(:)./dtau),ttemp(:),'linear',0)).*dt).'; F1=((0:length(Y1)-1)./(length(Y1)*dt_psd)).'; [nnfreq,nnind] = min(abs(F1-nizero)); Y1=Y1(nnind); clear F1 Lambda(i,n) = 2*(abs(max(Y1)).^2)/sum(xtemp(:)); %CREDO (oppure prendono la potenza massima?) % toc end % disp('1ni') end toc disp(length(f_gr)); disp(m); lamfiname = sprintf('Lambda_seg_%d.mat', m); % save(lamfiname,"Lambda"); save([pathfi,lamfiname],"Lambda"); Lambda = zeros(length(nibank),length(f_gr)); disp(m); end tic % nisearcher = KDTreeSearcher(nibank,'BucketSize',100); nisearcher = KDTreeSearcher(gather(nibank),'BucketSize',100); %% The line above will have to be substituted with something along the lines %% of the ones below to account for the metric g_jj in the phase derivatives % gdistance = @(a,b)sqrt(((a-b).^2)*(g_jj(1:s_s))); % nisearcher = ExhaustiveSearcher(gather(nibank),'Distance',gdistance); % nisearcher = ExhaustiveSearcher(nibank,'Distance',gdistance); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% parbank = combinations(porb_gr,a_gr,tasc_gr).Variables; bestpar = zeros(1,5); toc totlam = zeros(length(parbank),length(f_gr)); tic for n=1:length(f_gr) for i=1:length(parbank) curpar = [f_gr(n),parbank(i,:)]; % Move from ν,P,a,T_asc to ν,Ω,a,γ curpar(2) = 2*pi/curpar(2); for m = 1:M curpar(4) = curpar(2)*(tmid(m) - parbank(i,3)); curni=zeros(1,s_s); for s=1:s_s curni(s) = (curpar(2)^s)*sin(curpar(4)+0.5*s*pi); end curni = -curni.*curpar(1).*curpar(3); curni(1) = curni(1) + curpar(1); lamfiname = sprintf('Lambda_seg_%d.mat', m); % load(lamfiname); load([pathfi,lamfiname]); [Idx,D] = knnsearch(nisearcher,curni); totlam(i,n) = totlam(i,n)+Lambda(Idx,n); clear Lambda end % % if (totlam > bestpar(1)) % % bestpar = [totlam,f_gr(n),parbank(i,:)]; % % end end end toc % disp(bestpar); % save('C:\Users\Filippo\Desktop\XMM_Jxxx\risultelli.mat'); save([pathfi,'risultelli_',num2str(Tseg),'s.mat']); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%% Fai il grafico delle totlam a posteriori %%%%%%%%%%%%%%%% % for n = 1:length(f_gr) % hold on % txt = ['f_{gr} = ',num2str(f_gr(n))]; % plot(totlam(:,n),'DisplayName',txt) % end % legend %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% COMPLETARE DOPO AVER STUDIATO I RANDOM TEMP. BANKS %% E VEDERE CODICE https://pyfstat.readthedocs.io/en/latest/pyfstat.html#pyfstat.core.SemiCoherentSearch %-------------------------------------------------------------------------- %FUNCTIONS %-------------------------------------------------------------------------- function res = delta_ni(mu_s,s,g) %Compute eq. 23 M2015 res= 2.*sqrt(mu_s./(s.*g)); end No newline at end of file
