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

Pushing last version (now included and made better in scsearch.m)

parent 835a2298

test_efficiency.m

+131 −105
Original line number Diff line number Diff line
@@ -150,61 +150,27 @@ Omega_min = 2*pi/porb_max; 
% 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

% Is (1-Ome*a) more than fmin-fmax? to ask myself

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

nismin(s) = f_max*(1 - a_max*Omega_max);

nismax(s) = f_max*(1 + a_max*Omega_max);

for s = 2:4

    nismin(s) = -f_max*a_max*(Omega_max^s);

    nismax(s) =  f_max*a_max*(Omega_max^s);
@@ -253,7 +219,7 @@ for s = 1:s_s 
        Nni(s) = Nni(s)+1;

    end

    franco=nismin(s):((nismax(s)-nismin(s))/(Nni(s))):nismax(s);

    nis{s}=franco;

    nis{s}=franco/f_max;

end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -274,6 +240,12 @@ tic 
tm=gpuArray(tm);

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

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

F1=gpuArray((0:N-1)./(N*dt_psd)).';

f_ind = zeros(1,length(f_gr),'uint32');

for n = 1:length(f_gr)

    [bles,f_ind(n)] = min(abs(F1-f_gr(n)));

end

clear bles

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

% Lambda = zeros(length(nibank),length(f_gr));

toc
@@ -304,7 +276,9 @@ for m=1:M 
    ttemp = tm(m,:);

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

    ttempdifl = (ttemp(:)-tmid(m)-0.5*dt).';

    xtemp = x(:,m).';

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

    sumx = sum(xtemp);

    



    % toc
@@ -312,54 +286,45 @@ for m=1:M 
    % 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

    for i = 1:length(nibank)

    % for i = 1:1



        % Tento ricampionamento (Eq. 17 MP2015)

        % 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);

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

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

        dtau = ((taul(2)-taul(1)));



        % inzo = (interp1(tau(:),(xtemp(:)./dtau),ttemp(:),'linear',0)).*dt;

        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

        % Careful! ni_0 still needs to be defined (it's the spin freq. 

        % currently being considered)

        for n = 1:length(f_gr)



            Yenne=Y1(f_ind(n)); 

            Lambda(i,n) = 2*(abs(Yenne)^2)/sumx; %CREDO (oppure prendono la potenza massima?)

            

        end

        % disp('1ni')

    end

    toc

    disp(length(f_gr));

    disp(length(f_gr)*i);

    disp(m);

    lamfiname = sprintf('Lambda_seg_%d.mat', m);

    % save(lamfiname,"Lambda");

    save([pathfi,lamfiname],"Lambda");

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

    Lambda = zeros(length(nibank),length(f_gr));

    disp(m);



end





%% 





tic

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

nisearcher = KDTreeSearcher(gather(nibank),'BucketSize',100);
@@ -376,40 +341,101 @@ toc 


totlam = zeros(length(parbank),length(f_gr));

tic

for n=1:length(f_gr)

% 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(3);

%             curni(1) = curni(1) + 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

bles = cospi(1/2);

for m = 1:M

    tic

    lamfiname = sprintf('Lambda_seg_%d.mat', m);

    % load(lamfiname);

    load([pathfi,lamfiname]);

    toc

    disp(m);

    tic

    for i=1:length(parbank)

        curpar = [f_gr(n),parbank(i,:)];

        

        curpar = [bles,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]);



        curni = -curni.*curpar(3);

        curni(1) = curni(1) + 1;

        [Idx,D] = knnsearch(nisearcher,curni);

            totlam(i,n) = totlam(i,n)+Lambda(Idx,n);



            clear Lambda

        totlam(i,:) = totlam(i,:)+Lambda(Idx,:);



    end

        % % if (totlam > bestpar(1))

        % %     bestpar = [totlam,f_gr(n),parbank(i,:)];

        % % end

    end

    toc 

    disp(m);

    clear Lambda

end

% % 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']);

save([pathfi,'risultelli_07_03_',num2str(Tseg),'s.mat'],"-v7.3");

%% 





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