Changes to how the nearest neighbour is picked (min is faster than nnsearch),... (12c3d1c9) · Commits · HEAG-OAR / scsearch

functions/fastcounts.m

+4 −1

Original line number	Original line	Diff line number	Diff line
	function counts = fastcounts(times,tedges)		function counts = fastcounts(times,tedges)
	%FASTCOUNTS Histogram of array of times put in bins with edges tedges		%FASTCOUNTS Histogram of array of times in bins with tedges (use on RAM!)
	% For longish arrays it's much faster than hisctounts and uses only a		% For longish arrays it's much faster than hisctounts and uses only a
	% fraction of its memory; operating on sorted times are where it shines		% fraction of its memory; operating on sorted times are where it shines
	% but it still is faster in unsorted ones.		% but it still is faster in unsorted ones.
	@@ -30,6 +30,9 @@ counts = counts(1:Nbins);

	j = 1;		j = 1;
	i = 1;		i = 1;
			while times(j) < tedges(1)
			j=j+1;
			end
	while i < Nedges		while i < Nedges
	while (j<=Nphot && times(j)<tedges(i+1))		while (j<=Nphot && times(j)<tedges(i+1))
	counts(i) = counts(i)+1;		counts(i) = counts(i)+1;

scsearch.m

+67 −49

Original line number	Original line	Diff line number	Diff line
	@@ -2,8 +2,9 @@ clear;
	close all;		close all;

	import matlab.io.*		import matlab.io.*
			addpath('functions/')

	maxNumCompThreads(48);		maxNumCompThreads(24);

	% pathfi = '/data/Sorgenti/3FGLJ1544.6-1125/HST/';		% pathfi = '/data/Sorgenti/3FGLJ1544.6-1125/HST/';
	pathfi = '/data/Sorgenti/SCORPIUS-X-1/';		pathfi = '/data/Sorgenti/SCORPIUS-X-1/';
	@@ -445,12 +446,12 @@ for m=1:M
	% tedge = (tm(m,1)-1.5*dt-tmid)		% tedge = (tm(m,1)-1.5*dt-tmid)
	Y0 = fft(xtemp);		Y0 = fft(xtemp);

	graph1 = [pathfigu,'PS_',finame,'_m',num2str(m),'.pdf'];		graph1 = [pathfigu,'PS_',finame,'_Tseg',num2str(Tseg),'_m',num2str(m),'.pdf'];
	graph2 = [pathfigu,'CT_',finame,'_m',num2str(m),'.pdf'];		graph2 = [pathfigu,'CT_',finame,'_Tseg',num2str(Tseg),'_m',num2str(m),'.pdf'];
	if (~(isfile(graph1)) \|\| ~(isfile(graph2)))		if (~(isfile(graph1)) \|\| ~(isfile(graph2)))
	[VsuEtot,epsct,hmp] = powdist_f(Nyq,ttemp,false,graph1,graph2);		[VsuEtot,epsct,hmp] = powdist(Nyq,ttemp,false,graph1,graph2);
	end		end
	norman = gather(2./mean((abs(Y0(F1>1000 & F1<Nyq)).^2)));		norman = gather(2.*abs(Y0(1))./mean((abs(Y0(F1>1000 & F1<Nyq)).^2)));
	% norman = 2./(Nphot.*var(xtemp)./mean(xtemp));		% norman = 2./(Nphot.*var(xtemp)./mean(xtemp));
	tau = ones(size(ttemp),"like",ttemp);		tau = ones(size(ttemp),"like",ttemp);

	@@ -477,14 +478,14 @@ for m=1:M

	% Y0 = fft(xtemp);		% Y0 = fft(xtemp);
	% inzo = (interp1(tau(:),(xtemp(:)./dtau),ttemp(:),'linear',0)).*dt;		% inzo = (interp1(tau(:),(xtemp(:)./dtau),ttemp(:),'linear',0)).*dt;
	Y1 = fft(mat).';		Y1 = abs(fft(mat).');
	% sumo = real(Y1(1));		% sumo = real(Y1(1));
	% plot(F1,abs(Y0)./sumx,F1,abs(Y1)./sumo)		% plot(F1,abs(Y0)./sumx,F1,abs(Y1)./sumo)
	% pause		% pause
	% Careful! ni_0 still needs to be defined (it's the spin freq.		% Careful! ni_0 still needs to be defined (it's the spin freq.
	% currently being considered)		% currently being considered)
	% Yenne(:)=Y1(f_ind(:));		% Yenne(:)=Y1(f_ind(:));
	Lambda(:,i) = gather(Y1(f_ind(:)));		Lambda(:,i) = gather(Y1(f_ind(:))./sqrt(Y1(1)));
	% Lambda(:,i) = norman.*(abs(Yenne(:)).^2)./sumo;		% Lambda(:,i) = norman.*(abs(Yenne(:)).^2)./sumo;
	%		%
	% for n = 1:length(f_gr)		% for n = 1:length(f_gr)
	@@ -495,7 +496,7 @@ for m=1:M
	% end		% end
	% disp('1ni')		% disp('1ni')
	end		end
	Lambda = (abs(Lambda).^2).*norman;		Lambda = (Lambda.^2).*norman;
	% disp(max(Lambdacp(:,:)-Lambda(:,:)))		% disp(max(Lambdacp(:,:)-Lambda(:,:)))
	% disp(min(Lambdacp(:,:)-Lambda(:,:)))		% disp(min(Lambdacp(:,:)-Lambda(:,:)))
	% toc		% toc
	@@ -558,12 +559,12 @@ else
	disp(['The length of nibank in greater than ',num2str(intmax('uint64')),', something definitely makes no sense here'])		disp(['The length of nibank in greater than ',num2str(intmax('uint64')),', something definitely makes no sense here'])
	end		end

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

	bles = cospi(1/2);		bles = cospi(1/2);
	@@ -572,7 +573,7 @@ for m = 1:M
	segstart = tic;		segstart = tic;
	lamfiname = sprintf('Lambda_fmax_%d_seg_%d.mat', f_max, m);		lamfiname = sprintf('Lambda_fmax_%d_seg_%d.mat', f_max, m);
	load([pathfi,lamfiname]);		load([pathfi,lamfiname]);
	for i=1:length(parbank)		for i=1:parno

	curpar = [bles,parbank(i,:)];		curpar = [bles,parbank(i,:)];
	% Move from ν,P,a,T_asc to ν,Ω,a,γ		% Move from ν,P,a,T_asc to ν,Ω,a,γ
	@@ -585,7 +586,9 @@ for m = 1:M
	curni = -curni.*curpar(3);		curni = -curni.*curpar(3);
	curni(1) = curni(1) + 1;		curni(1) = curni(1) + 1;

	[Idx,D] = knnsearch(nisearcher,curni);		% [Idx,D] = knnsearch(nisearcher,curni);
			[schi3,fo3] = min(abs(nibank(:,:)-curni),[],1);
			Idx = sum(fo3)-length(fo3)+1;

	% nimask(Idx) = 1;		% nimask(Idx) = 1;

	@@ -607,19 +610,19 @@ totlam = totlam - 1;

	%% Calculate detection threshold for a multi-trial false alarm prob. = pfa		%% Calculate detection threshold for a multi-trial false alarm prob. = pfa
	pfa = 0.01;		pfa = 0.01;
	sigmastar = 2gammaincinv((1-(1-pfa)^(1/(length(parbank)length(f_gr)))),M,'upper');		sigmastar = 2gammaincinv((1-(1-pfa)^(1/(parnolength(f_gr)))),M,'upper');
	%sigmastar = 2gammaincinv((1-(1-pfa)^(1/(sum(nimask)length(f_gr)))),M,'upper');		%sigmastar = 2gammaincinv((1-(1-pfa)^(1/(sum(nimask)length(f_gr)))),M,'upper');
	[bru,tto] = max(totlam,[],2);		[bru,tto] = max(totlam,[],2);
	[maxtotlam,maxinde] = max(bru);		[maxtotlam,maxinde] = max(bru);
	bestpar = [f_gr(maxinde), parbank(tto(maxinde),:), double(maxtotlam)];		bestpar = [f_gr(maxinde), parbank(tto(maxinde),:), double(maxtotlam)];
	bparch = string(bestpar);		bparch = string(bestpar);
	tascmjd = MJDREF+tasc_gr./86400;		tascmjd = MJDREF+tasc_gr./86400;
	pfamax = 1.0-( (1.0-gammainc(0.5double(maxtotlam),M,'upper'))^(length(parbank)length(f_gr)));		pfamax = 1.0-( (1.0-gammainc(0.5double(maxtotlam),M,'upper'))^(parnolength(f_gr)));
	disp('Our most significant peak has parameters:');		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(['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+bestpar(4)/86400)),'.']);
	disp(['Its detection statistics is Σ_max = ',num2str(maxtotlam),',']);		disp(['Its detection statistics is Σ_max = ',num2str(maxtotlam),',']);
	disp(['which considering all parameters, corresponds to a false alarm probability of ',num2str(pfamax),','])		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),'.']);		disp(['while the multi-trial threshold for a false alarm probability = ',num2str(pfa),' is eq. to ',num2str(sigmastar),'.']);
	truwalk = size(unique(niwalk.', 'rows'),1);		truwalk = size(unique(niwalk.', 'rows'),1);
	sigmastar2 = 2gammaincinv((1-(1-pfa)^(1/(truwalklength(f_gr)))),M,'upper');		sigmastar2 = 2gammaincinv((1-(1-pfa)^(1/(truwalklength(f_gr)))),M,'upper');
	@@ -659,8 +662,13 @@ diary off
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%%%%%%%%%% Fai il grafico delle totlam a posteriori %%%%%%%%%%%%%%%%		%%%%%%%%%%%%% Fai il grafico delle totlam a posteriori %%%%%%%%%%%%%%%%

			figu = figure('visible','on');
			% figu = figure('visible','on');
			rina = [pathfigu,'f_maxtotlam.m'];
	plot(f_gr,bru)		plot(f_gr,bru)
	yline(sigmastar,'--','Σ*')		yline(sigmastar,'--','Σ*')
			yline(sigmastar2,'--','Σ*2')
			saveas(figu,rina);

	% for n = 1:length(f_gr)		% for n = 1:length(f_gr)
	% hold on		% hold on
	@@ -675,41 +683,51 @@ yline(sigmastar,'--','Σ*')
	% contour(a_gr,tasc_gr,squeeze(max(provaccia,[],3)),[20:20:120])		% contour(a_gr,tasc_gr,squeeze(max(provaccia,[],3)),[20:20:120])
	% contour(a_gr,porb_gr,squeeze(max(provaccia,[],1)).',[20:20:120])		% contour(a_gr,porb_gr,squeeze(max(provaccia,[],1)).',[20:20:120])

	% provaccia = reshape(totlam.',length(parbank)*(length(f_gr)),1);		if (inGB<350)
	% provaccia = reshape(provaccia,length(tasc_gr),length(a_gr),length(porb_gr),length(f_gr));		provaccia = reshape(totlam.',length(parbank)*(length(f_gr)),1);
	% figstep = (10);		provaccia = reshape(provaccia,length(tasc_gr),length(a_gr),length(porb_gr),length(f_gr));
			minfiggrid = min(maxtotlam0.9,sigmastar20.9);
			maxfiggrid = max(maxtotlam1.1,sigmastar21.1);
			figstep = (maxfiggrid-minfiggrid)/4;
			fig_grid = [minfiggrid:figstep:maxfiggrid];
			figu = figure('visible','on');
			% figu = figure('visible','on');
			rina = [pathfigu,'porb_tasc.m'];
			contourf(porb_gr,tascmjd,squeeze(max(provaccia,[],[2 4])),fig_grid,'ShowText','on')
			hold on
			scatter(bestpar(2),(MJDREF+bestpar(4)/86400))
			if (exist('trupar') == 1)
			scatter(trupar(2),trupar(4),'+')
			end
			hold off
			saveas(figu,rina);
	% figu = figure('visible','off');		% figu = figure('visible','off');
	% rina = "porb_tasc.m";		figu = figure('visible','on');
	% contourf(porb_gr,tascmjd,squeeze(max(provaccia,[],[2 4])),[sigmastar0.9:figstep:1.1maxtotlam],'ShowText','on')		rina = [pathfigu,'porb_a.m'];
	% hold on		contourf(porb_gr,a_gr,squeeze(max(provaccia,[],[1 4])),fig_grid,'ShowText','on')
	% scatter(bestpar(2),(MJDREF+bestpar(4)/86400))		hold on
	% if (exist('trupar') == 1)		scatter(bestpar(2),bestpar(3))
	% scatter(trupar(2),trupar(4),'+')		if (exist('trupar') == 1)
	% end		scatter(trupar(2),trupar(3),'+')
	% hold off		end
	% saveas(figu,rina);		hold off
	% figu = figure('visible','off');		saveas(figu,rina);
	% rina = "porb_a.m";
	% contourf(porb_gr,a_gr,squeeze(max(provaccia,[],[1 4])),[sigmastar:figstep:1.1*maxtotlam],'ShowText','on')
	% hold on
	% scatter(bestpar(2),bestpar(3))
	% if (exist('trupar') == 1)
	% scatter(trupar(2),trupar(3),'+')
	% end
	% hold off
	% saveas(figu,rina);
	% figu = figure('visible','off');		% figu = figure('visible','off');
	% rina = "f_porb.m";		figu = figure('visible','on');
	% contourf(f_gr,porb_gr,squeeze(max(provaccia,[],[1 2])),[sigmastar:figstep:1.1*maxtotlam],'ShowText','on')		rina = [pathfigu,'f_porb.m'];
	% hold on		contourf(f_gr,porb_gr,squeeze(max(provaccia,[],[1 2])),fig_grid,'ShowText','off')
	% scatter(bestpar(1),bestpar(2))		hold on
	% if (exist('trupar') == 1)		scatter(bestpar(1),bestpar(2))
	% scatter(trupar(1),trupar(2),'+')		if (exist('trupar') == 1)
	% end		scatter(trupar(1),trupar(2),'+')
	% hold off		end
	% saveas(figu,rina);		hold off
			saveas(figu,rina);
			end

			% if (ceil(timerapp)<11)
	% exit		% exit
			% end

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

scsearch_1023.m

0 → 100644

+818 −0

File added.

Preview size limit exceeded, changes collapsed.