Commit 01513018 authored by riccardo's avatar riccardo
Browse files

Added first version of lambdasum, trivial changes to scsearch

parent 4ff8802f
Loading
Loading
Loading
Loading

functions/lambdasum.m

0 → 100644
+64 −0
Original line number Diff line number Diff line
function [outputArg1,outputArg2] = lambdasum(inputArg1,inputArg2,ParBank,NiSearcher)
%lambdasum Sums the Lambdas for all parameter combinations over obs.s/seg.s
%   Detailed explanation goes here


% 
% ParBank = combinations(porb_gr,a_gr,tasc_gr).Variables;

parno = length(ParBank);
freqno = length(f_gr);
nino = length(NiSearcher.X);
if (nino<intmax('uint16'))
    niwalk = ones(M,parno,'uint16');
elseif (nino<intmax('uint32'))
    niwalk = ones(M,parno,'uint32');
elseif (nino<intmax('uint64'))
    niwalk = ones(M,parno,'uint64');
else
    disp(['The length of nibank in greater than ',num2str(intmax('uint64')),', something definitely makes no sense here'])
end

disp(['length(parbank) = ',num2str(parno),'']);
% Writing on matrix of ones is faster in MATLAB than on zeroes, we subtract
% one after all calculations are done
totlam = ones(freqno,parno,'single');
guess = parno*M*16192.4539/(4.0*7546032.0);
disp(['(Probably terrible) Estimate of the time needed for the second loop = ',num2str(ceil(guess)),' s']);  disp(' ');

avetime = 0.0;
for m = 1:M
    segstart = tic;
    lamfiname = sprintf('Lambda_fmax_%d_seg_%d.mat', f_max, m);
    load([pathfi,lamfiname],"Lambda","");
    for i=1:parno
        
        curpar = ParBank(i,:);
        % Move from ν,P,a,T_asc to ν,Ω,a,γ
        % Move from P,a,T_asc to Ω,a,γ
        curpar(1) = 2*pi/curpar(1);
        curpar(3) = curpar(1)*(tmid(m) - curpar(3));
        curni=zeros(1,s_s);
        for s=1:s_s
            curni(s) = (curpar(1)^s)*sin(curpar(3)+0.5*s*pi);
        end
        curni = -curni.*curpar(2);
        curni(1) = curni(1) + 1;

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

        niwalk(m,i) = Idx;

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

    end
    segend = toc(segstart);
    avetime = avetime + segend;
    disp(['Just finished segment num.: ',num2str(m),'.']);
    disp(['Elapsed time since loop began = ',num2str(avetime),' s, estimated time left = ',num2str(((M-m)*avetime/m)),' s.']); disp(' ');
    clear Lambda
end

outputArg1 = inputArg1;
outputArg2 = inputArg2;
end
 No newline at end of file
+12 −66
Original line number Diff line number Diff line
@@ -302,36 +302,13 @@ nodesnum = length(porb_gr)*length(a_gr)*length(tasc_gr)*length(f_gr)
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);

% 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
% We take 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;
nismin(s) = (1 - a_max*Omega_max);
nismax(s) = (1 + a_max*Omega_max);
@@ -403,8 +380,8 @@ tic
% vari ni1,ni2,...,nis_s
%%%[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);
nibank = combinations(nis{:}).Variables;
nifax = gpuArray(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');
@@ -470,7 +447,7 @@ for m=1:M

    graph1 = [pathfigu,'PS_',finame,'_m',num2str(m),'.pdf'];
    graph2 = [pathfigu,'CT_',finame,'_m',num2str(m),'.pdf'];
    if (~(isfile(graph1)) | ~(isfile(graph2)))
    if (~(isfile(graph1)) || ~(isfile(graph2)))
        [VsuEtot,epsct,hmp] = powdist_f(Nyq,ttemp,false,graph1,graph2);
    end
    norman = gather(2./mean((abs(Y0(F1>1000 & F1<Nyq)).^2)));
@@ -582,8 +559,9 @@ else
end

disp(['length(parbank) = ',num2str(length(parbank)),'']);
nimask = false(length(nibank),1);
% totlam = zeros(length(f_gr),length(parbank),'single');
% nimask = false(length(nibank),1);
% Writing on matrix of ones is faster in MATLAB than on zeroes, we subtract
% one after all calculations are done
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(ceil(guess)),' s']);  disp(' ');
@@ -594,9 +572,6 @@ for m = 1:M
    segstart = tic;
    lamfiname = sprintf('Lambda_fmax_%d_seg_%d.mat', f_max, m);
    load([pathfi,lamfiname]);
    % toc
    % disp(m);
    % tic
    for i=1:length(parbank)
        
        curpar = [bles,parbank(i,:)];
@@ -611,54 +586,25 @@ for m = 1:M
        curni(1) = curni(1) + 1;

        [Idx,D] = knnsearch(nisearcher,curni);
        % if (~nimask(Idx))
        nimask(Idx) = 1;
        % end

        % nimask(Idx) = 1;

        niwalk(m,i) = Idx;

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

    end
    % toc 
    segend = toc(segstart);
    avetime = avetime + segend;
    disp(['Just finished segment num.: ',num2str(m),'.']);
    disp(['Elapsed time since loop began = ',num2str(avetime),' s, estimated time left = ',num2str(((M-m)*avetime/m)),' s.']); disp(' ');
    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
%% 

% totlam was initialized as ones
totlam = totlam - 1;

%% Calculate detection threshold for a multi-trial false alarm prob. = pfa
pfa = 0.01;
sigmastar = 2*gammaincinv((1-(1-pfa)^(1/(length(parbank)*length(f_gr)))),M,'upper');