Commit 96379c7b authored by Riccardo La Placa's avatar Riccardo La Placa
Browse files

changesss

parent e550d421

scsearch_1023_multiobs.m

+16 −11
Original line number Diff line number Diff line
@@ -6,7 +6,7 @@ addpath('functions/') 


maxNumCompThreads(48);



pathfi = '/data/Sorgenti/J1023/30-31_01_2020/';

pathfi = '/data/Sorgenti/J1023/03_2016/';

folname = [pathfi,'SP_search_',char(datetime('now','Format','dd_MM'))];

figfolname = [folname,'/figures'];

if (exist(folname,"dir")==7)
@@ -77,8 +77,8 @@ filoc = [pathfi,'../',finame]; 
file=fits.openFile(filoc);

mjdref=fits.readKey(file,'MJDREF');

MJDREF(1)=str2double(mjdref);

tasc_min = (58878.9911009 - MJDREF(1)).*86400 - porb_max/2

tasc_max = (58878.9911009 - MJDREF(1)).*86400 + porb_max/2

tasc_min = (57449.72579 - MJDREF(1)).*86400 - porb_max/2

tasc_max = (57449.72579 - MJDREF(1)).*86400 + porb_max/2

tasc_mid = (tasc_max + tasc_min)/2

fits.movAbsHDU(file,2);

tstart = str2double(fits.readKey(file,'TSTART'));
@@ -112,7 +112,7 @@ countsname = ['counts_',finame,'_Tseg_',num2str(Tseg),'.mat']; 
save([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw","tstart","nbin","-v7.3","-nocompression");



toff = tstart - tasc_mid;

delporb(a_max,porb_min,f_max,toff,Tseg,M,tol);

% delporb(a_max,porb_min,f_max,toff,Tseg,M,tol);

% (porb_min^2)*sqrt(0.1*(1-tol))/(a_max*2*pi*Tseg*f_max);

% % stepparino = max(delporb(a_max,porb_min,f_max,toff,Tseg,M,tol),(porb_min^2)*sqrt(0.1*(1-tol))/(a_max*2*pi*Tseg*f_max))

% % porb_step = (porb_max-porb_min)/ceil((porb_max-porb_min)/stepparino)
@@ -267,6 +267,7 @@ nifax = gpuArray(nibank(:,1:s_s).*infax(1:s_s)); 
nino = length(nibank);

disp(['length(nibank) = ',num2str(nino),'']);

s_s

Mtot

% Lambda = zeros(length(f_gr),length(nibank),'single','gpuArray'); % MATLAB

% stores array elements in a column-major layout

Lambda = ones(length(f_gr),length(nibank),'single');
@@ -281,6 +282,7 @@ for o = 1:numfiles 
        finame = char(mario(o))

        countsname = ['counts_',finame,'_Tseg_',num2str(Tseg),'.mat'];

        load([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw","tstart","nbin");

        Lambda = ones(length(f_gr),length(nibank),'single');

    end

    tm=gpuArray(tm);

    tedge = ones(1,N+1,'gpuArray');
@@ -340,7 +342,7 @@ for o = 1:numfiles 
        end

        Lambda = (Lambda.^2).*norman;



        lamfiname = ['Lambda_',finame,'_fmax_',num2str(f),'_seg_',num2str(m),'.mat'];

        lamfiname = ['Lambda_',finame,'_fmax_',num2str(f_max),'_seg_',num2str(m),'.mat'];

        save([pathfi,lamfiname],"Lambda","-v7.3","-nocompression");

        Lambda = ones(size(Lambda),"like",Lambda);

        % Lambda = zeros(length(f_gr),length(nibank),'single','gpuArray');
@@ -354,11 +356,10 @@ for o = 1:numfiles 


    %% The first time in this loop we need to set some general variables such as niwalk and totlam 

    if (o == 1)

        nodesnum = length(porb_gr)*length(a_gr)*length(tasc_gr)*length(f_gr);

        nodesnum = length(porb_gr)*length(a_gr)*length(tasc_gr)*length(f_gr)

        nodesmem = (nodesnum*8.0)/1024^3; %->to bytes -> to GBs



        memthres = 600.0 - (length(f_gr)*parno*8.0)/1024^3;



        memthres = 600.0 - (length(f_gr)*nino*8.0)/1024^3;

        if (nodesmem > memthres)

            % THIS DISCLAIMER IS STILL HERE AFTER MONTHS

            % FIX IIIIT
@@ -405,7 +406,7 @@ for o = 1:numfiles 
    avetime = 0.0;

    for m = 1:M

        segstart = tic;

        lamfiname = ['Lambda_',finame,'_fmax_',num2str(f),'_seg_',num2str(m),'.mat'];

        lamfiname = ['Lambda_',finame,'_fmax_',num2str(f_max),'_seg_',num2str(m),'.mat'];

        load([pathfi,lamfiname]);

        % toc

        % disp(m);
@@ -455,11 +456,11 @@ sigmastar = 2*gammaincinv((1-(1-pfa)^(1/(parno*length(f_gr)))),Mtot,'upper'); 
[maxtotlam,maxinde] = max(bru);

bestpar = [f_gr(maxinde), parbank(tto(maxinde),:), double(maxtotlam)]; 

bparch = string(bestpar);

tascmjd = MJDREF+tasc_gr./86400;

tascmjd = MJDREF(1)+tasc_gr./86400;

pfamax = 1.0-( (gammainc(0.5*double(maxtotlam),Mtot,'lower'))^(parno*length(f_gr)));

disp('Our most significant peak has parameters:');

disp(['f_spin = '+bparch(1)+' Hz, p_orb = '+bparch(2)+' s, a*sin(i)/c = '+bparch(3)+' s, t_asc = '+bparch(4)+' s.']);

disp(['t_asc in days is = ',num2str((MJDREF+bestpar(4)/86400)),'.']);

disp(['t_asc in days is = ',num2str((MJDREF(1)+bestpar(4)/86400)),'.']);

disp(['Its detection statistics is Σ_max = ',num2str(maxtotlam),',']);

disp(['which considering all parameters (',num2str(parno),') corresponds to a false alarm probability of ',num2str(pfamax),','])

disp(['while the multi-trial threshold for a false alarm probability = ',num2str(pfa),' is eq. to ',num2str(sigmastar),'.']);
@@ -471,6 +472,10 @@ disp(['The number of combinations of ni combinations actually walked is ',num2st 
% pfamax2 = 1.0-( (1.0-gammainc(0.5*double(maxtotlam),M,'upper'))^(sum(nimask)*length(f_gr)));

disp(['therefore, when considering just them instead, Σ_max corresponds to a false alarm probability of ',num2str(pfamax2),','])

disp(['while the multi-trial threshold for a false alarm probability = ',num2str(pfa),' is eq. to ',num2str(sigmastar2),'.']);

% To accurately (I hope) calculate the false alarm prob. use Mathematica

% with the following notation (here 73 is Mtot, 1860 is maxtotlam and 282

% is the number of digits to be used to obtain the numerical value):

% N[1- (1-GammaRegularized[73, 1860/2])^(30*12801*11886628), 282]



memused = whos;

summina = 0;