Commit 1c52a5a9 authored by Riccardo La Placa's avatar Riccardo La Placa
Browse files

Optimized tau, changed metric

parent 5dc01cda

power_dist_test.m

+11 −2
Original line number Diff line number Diff line
@@ -6,6 +6,7 @@ diary(['log_pow_dist_',char(datetime('now','Format','dd_MM')),'.log']) 
mario = split(ls('*.fits'));

mario = mario(~cellfun('isempty',mario));

mario = sort(mario);

% mario = 'SiFAP2_20220429_3FGLJ1544_N1_Bary.fits'

import matlab.io.*

Nyq = 2000

Tseg = 512
@@ -33,6 +34,7 @@ epsvalues = zeros(length(mario),100); 
varvalues = zeros(length(mario),100,3);

for i = 1:length(mario)

    filoc = [pwd,'/',char(mario(i))]

    % filoc = [pwd,'/',mario]

    % file = fits.openFile(filoc)

    % fits.closeFile(file);

    t_raw = fitsread(filoc,"binarytable");
@@ -55,7 +57,7 @@ for i = 1:length(mario) 
        end

        Yqlbin(:) = Yqlbin(:)./(Fqlbin(:).*Tseg);



        if (i == 4 || i == 6 || i == 7)

        if (i == 7)

            tiledlayout(2,1)

            nexttile
@@ -161,7 +163,8 @@ for i = 1:length(mario) 
            extrtot1 = ptotk(p1_1,extrpoiss);

            



            if (i == 4 || i == 6 || i == 7)



            if (i == 7)

            tiledlayout(2,1)

            nexttile

            semilogy([1:maxcount],extrtot1,[1:maxcount],extrtot4,[1:maxcount],extrtot8,[1:maxcount],datadist,'+')
@@ -188,6 +191,12 @@ epsnans(epsnans == 0) = NaN; 
epsnans(epsnans > 0.1) = NaN;

minans = mivalues;

minans(minans == 0) = NaN;



plot(mean(minans,2,"omitmissing"),mean(epsnans,2,"omitmissing"),'+')

xlim([0 max(minans,[],"all")])

ylim([0 max(epsnans,[],"all")])

text(mean(minans,2,"omitmissing"),mean(epsnans,2,"omitmissing"),mario,"Interpreter",'none')

ylim([min(epsnans,[],"all") max(epsnans,[],"all")])

%% 

% disp('Finished with the fits files, move to dark obs.?')

% pause

scsearch.m

+24 −14
Original line number Diff line number Diff line
@@ -195,8 +195,8 @@ end 


% Parameters for Sco X-1

 % s 

f_max = 1000

f_min = 950

f_max = 1050

f_min = 948

f_step = 1/Tseg

f_gr = (f_min:f_step:f_max).';

Porb = 68023.919 % s
@@ -369,6 +369,7 @@ for s = 1:s_s 
    franco=nismin(s):((nismax(s)-nismin(s))/(Nni(s))):nismax(s);

    nis{s}=franco;

end

infax = gpuArray(1./factorial(1:s_s)); % MATLAB's factorial is slow



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Changing nis to stop at s*(= s_s): this way nibank is length(f_gr)
@@ -387,6 +388,7 @@ tic 
%%%[kung,fu,fight] = ndgrid(nis{1},nino{2},nino{3});

tm=gpuArray(tm);

nibank = gpuArray(combinations(nis{:}).Variables);

nifax = nibank(:,1:s_s).*infax(1:s_s);

% 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','gpuArray');
@@ -411,8 +413,8 @@ nino = length(nibank); 
s_s

M

avetime = 0.0;

guess = nino*M*3071.5583/(3100.0*28);

disp(['Rough estimate of the time needed for the first loop = ',num2str(guess),' s']);  disp(' ');

guess = nino*M*Tseg*450/(192.0*31*256);

disp(['Rough estimate of the time needed for the first loop = ',num2str(ceil(guess)),' s']);  disp(' ');

%%

for m=1:M
@@ -435,9 +437,11 @@ for m=1:M 
    % ttemp = tm(m,:);

    ttemp = gpuArray(t_raw((t_raw>=gather(tm(m,1))-0.5*dt) & t_raw<=gather(tm(m,N))+0.5*dt));

    Nphot = length(ttemp);

    ttempdif = (ttemp(:)-tmid(m)).';

    ttempdifl(1:N) = (tm(m,1:N)-tmid(m)-0.5*dt).';

    ttempdifl(N+1) = (tm(m,N)-tmid(m)+0.5*dt);

    ttempdifs = zeros(Nphot,s_s,'gpuArray');

    ttempdifs(:,1) = (ttemp(:)-tmid(m)).';

    for s = 2:s_s

        ttempdifs(:,s) = ttempdifs(:,1).^s;

    end 

    % xtemp = (x(:,m).');

    xtemp = gpuArray(x(:,m).');

    % ttemp = tmrado(m,:);
@@ -445,9 +449,10 @@ for m=1:M 
    % tedge(N+1) = ttemp(end) + 0.5*dt;

    tedge(1:N) = tm(m,1:N) - 0.5*dt -tmid(m);  

    tedge(N+1) = tm(m,N) + 0.5*dt -tmid(m);

    % tedge = (tm(m,1)-1.5*dt-tmid)

    Y0 = fft(xtemp);

    % norman = 2./mean((abs(Y0(F1>1000 & F1<Nyq)).^2)./sumx);

    norman = 2./(Nphot.*var(xtemp)./mean(xtemp));

    norman = 2./mean((abs(Y0(F1>1000 & F1<Nyq)).^2));

    % norman = 2./(Nphot.*var(xtemp)./mean(xtemp));

    tau = ones(size(ttemp),"like",ttemp);
@@ -463,7 +468,7 @@ for m=1:M 
        % Per prima cosa, generiamo la nuova coordinata temporale per 

        % questa templ combini 

        % tau(1:Nphot) = tmid(m) + sum((nibank(i,1:s_s)./(factorial(1:s_s))).*((ttempdif(1:Nphot)).^((1:s_s).').'),2);

        tau(1:Nphot) = sum((nibank(i,1:s_s)./(factorial(1:s_s))).*((ttempdif(1:Nphot)).^((1:s_s).').'),2);

        tau(1:Nphot) = sum((nifax(i,1:s_s)).*(ttempdifs(1:Nphot,1:s_s)),2);

        % taul(1:N+1) = tmid(m) + sum((nibank(i,1:s_s)./(factorial(1:s_s))).*((ttempdifl(1:N+1)).^((1:s_s).').'),2);

        % dtau(1:N) = (taul(2:N+1)-taul(1:N));
@@ -515,12 +520,17 @@ end 
%%

%tic

% nisearcher = KDTreeSearcher(nibank,'BucketSize',100);

nisearcher = KDTreeSearcher(gather(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)));

sqg = zeros(s_s,1);

for s = 1:s_s

   sqg(s) = ((Tseg^(s-1))/factorial(s-1));

end

sqdist = @(x,Z)sqrt(((x-Z).^2)*sqg);

% nisearcher = ExhaustiveSearcher(gather(nibank),'Distance',gdistance);

% nisearcher = ExhaustiveSearcher(nibank,'Distance',gdistance);

nisearcher = ExhaustiveSearcher(nibank,'Distance',sqdist);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -554,7 +564,7 @@ nimask = false(length(nibank),1); 
% totlam = zeros(length(f_gr),length(parbank),'single');

totlam = ones(length(f_gr),length(parbank),'single');

guess = length(parbank)*M*16192.4539/(4.0*7546032.0);

disp(['(Probably terrible) Estimate of the time needed for the second loop = ',num2str(guess),' s']);  disp(' ');

disp(['(Probably terrible) Estimate of the time needed for the second loop = ',num2str(ceil(guess)),' s']);  disp(' ');



bles = cospi(1/2);

avetime = 0.0;
@@ -731,7 +741,7 @@ yline(sigmastar,'--','Σ*') 
% hold off

% saveas(figu,rina);



exit

% exit



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%