minor updates and improvments

README.md

@@ -3,18 +3,27 @@ SpESO
 
 Spectrum Estimation by Sparse Optimization (SpESO). MATLAB code for energy spectrum estimation of turbulence by compressive sampling (compressed sensing). The core of the method is Quasi-Oracle Multilevel Orthogonal Matching Pursuit (QOMOMP), which is a multilevel variation of OMP using a priori information about the signal.
 
-Requires a wavelet package (MATLAB's Wavelet Toolbox, WaveLab) but can also be used with the supplied CDF 9/7 transform. Package type is set with the DWTTYPE parameter.
-
 
 About
 ----------------------
 
-Written by Gudmundur F. Adalsteinsson (2013-2014). This code is intended as a supplement to a journal article and a thesis to promote reproducible research in computational science, the third branch of science, the other two being deductive and empirical science.
+Written by Gudmundur F. Adalsteinsson (2013-2014). This code is intended as a supplement to a journal article and a thesis to promote reproducible research in computational science.
+
+Copyright (c) 2014 Gudmundur Adalsteinsson (MIT License). See file LICENSE for license and warranty details.
+
+Installation and requirements
+----------------------
+
+To use, it is simplest to change the current directory to the folder SpESO. The main files can then be executed directly (a path command temporarily adds the functions folder to the current path). The execution time might be about an hour for SpESO. Once complete, the plotting scripts can be used to visualize the saved results (which are in .mat files).
+
+Requires a wavelet package (MATLAB's Wavelet Toolbox, WaveLab) but can also be used with the supplied CDF 9/7 transform. Package type is set with the DWTTYPE parameter.
 
 
 Reproducible Research
 ----------------------
 
+The classical branches of science are deductive and empirical science. In the last few decades, a third branch has emerged: computational science.
+
 "Scientific Computation is emerging as absolutely central to the scientific method. Unfortunately,
 it is error-prone and currently immature: traditional scientific publication is
 incapable of finding and rooting out errors in scientific computation; this must be recognized

functions/cs/QOMOMPt_2d_gfa.m

@@ -21,7 +21,6 @@
 x=zeros(Ns,Ns);
 ind=1:pow2(J0);
 ind=xyindex2x(ind,ind,Ns,Ns);
-%if (nargin < 6), L=1; end
 if (nargin < 7), tol=1e-4; end  % default is 1e-6 for lsqr
 if (nargin < 8), tolf=1e-9; end % final iteration
 L(L<0)=0; % negative Lj = 0
@@ -34,18 +33,14 @@
 x0=x(ind)';
 b=psyt(ind); b=b(:);
 [xt,flag] = symmlq(@(xx) itfunc(xx,Psi,Psit,ind,Ns),b,tol,200,[],[],x0);
-%[xt,~]  = tfqmr(@(xx) itfunc(xx,Psi,Psit,ind,N),b,tol,200,[],[],x0);
 
 % solution of oracle
 x(ind)=xt;
 r=y-Psi(x);
 
 
 for j=J0:log2(Ns)-1
-	%fprintf('QOMOMP j=%d \n',j);
-	%for stage=1:nstage
-
-
+
 	itdiff = xyindex2x(1:2^(j-1),1:2^(j-1),Ns,Ns); %to remove
 	omeg = setdiff(ind,itdiff,'stable'); % Omega_(j-1)
 	omegs=find(abs(x)>alph(x(omeg))); 
@@ -62,12 +57,7 @@
 
 
 	%============
-
-
-% 	% index set
-% 	temp=abs(Psit(r)); %figure;plot(x),figure;plot(temp)
-% 	temp2=temp(1:2^(j+1),1:2^(j+1));
-
+
 
 	[~,i] = sort(temp2(:));
 	itdiff = xyindex2x(1:2^j,1:2^j,pow2(j+1),pow2(j+1)); %to remove
@@ -109,22 +99,17 @@
 	b=psyt(ind); b=b(:);
 
 	[xt,~]  = symmlq(@(xx) itfunc(xx,Psi,Psit,ind,Ns),b,tol,200,[],[],x0);
-	%[xt,~]  = tfqmr(@(xx) itfunc(xx,Psi,Psit,ind,N),b,tol,200,[],[],x0);
-	%[xt,~]
-
+
 	% update
 	x(:)=0;
 	x(ind)=xt;
 	r=y-Psi(x);
 
-	%plot(x),xlim([-5,100])
-%end
+
 end
-%fprintf('QOMOMP final \n',j);
 
 % one final low tolerance iteration
 [xt,flag,relres,iter] = symmlq(@(xx) itfunc(xx,Psi,Psit,ind,Ns),b,tolf,200,[],[],x0);
-%[xt,flag,relres,iter] = tfqmr(@(xx) itfunc(xx,Psi,Psit,ind,N),b,tolf,200,[],[],x0);
 x(:)=0;
 x(ind)=xt;
 

functions/cs/QOMOMPt_gfa.m

@@ -18,7 +18,6 @@
 psyt=Psit(y);
 x=zeros(N,1);
 ind=1:2^J0;
-%if (nargin < 6), L=1; end
 if (nargin < 8+1), tol=1e-4; end  % default is 1e-6 for lsqr
 if (nargin < 9+1), tolf=1e-9; end % final iteration
 L(L<0)=0; % negative Lj = 0
@@ -32,17 +31,12 @@
 x0=x(ind);
 b=psyt(ind);
 
-
-%b=sign(b).*(abs(b)).^(0.8);
 [xt,~] = tfqmr(@(xx) itfunc(xx,Psi,Psit,ind,N),b,tol,200,[],[],x0);
-%xt = nsoli(x0,@(xx) nlfunc(xx,Psi,Psit,b,ind,N),[tol,0]);
-
 
 % solution of oracle
 x(ind)=xt;
 r=y-Psi(x);
 
-
 for j=J0:log2(N)-1
 	% index set
 	omeg=intersect(ind,2^(j-1)+1:2^(j)); % Omega_(j-1)
@@ -51,39 +45,30 @@
 	omeg=[omegs*2,omegs*2-1]; % Omega_(j)
 
 
-	temp=abs(Psit(r)); %figure;plot(x),figure;plot(temp)
+	temp=abs(Psit(r)); 
 	temp(omeg)=temp(omeg)*beta(j);
 	[~,i] = sort(temp(2^j+1:2^(j+1)));
 	i=i(end-L(j)+1:end) + 2^j;
 	ind=union(ind,i);
 
-	% figure,plot_coef_im(abs(x)>0,1:length(x))
-
 	% iterative method (with pseudo-inverse LS) no transpose
 	% (bicgstab, bicgstabl, cgs, gmres, minres, pcg, symmlq, tfqmr)
 	x0=x(ind);
 	b=psyt(ind);
 
 
-	%b=sign(b).*(abs(b)).^(0.8);
 	[xt,~]  = tfqmr(@(xx) itfunc(xx,Psi,Psit,ind,N),b,tol,200,[],[],x0);
-	%xt = nsoli(x0,@(xx) nlfunc(xx,Psi,Psit,b,ind,N),[tol,0]);
 
 	% update
 	x(1:end)=0;
 	x(ind)=xt;
 	r=y-Psi(x);
 
-	%plot(x),xlim([-5,100])
+
 end
 
 % one final low tolerance iteration
-%[xt,~] = tfqmr(@(xx) itfunc(xx,Psi,Psit,ind,N),b,tolf,200,[],[],x0);
 [xt,flag,relres,iter] = tfqmr(@(xx) itfunc(xx,Psi,Psit,ind,N),b,tolf,200,[],[],x0);
-%xt = nsoli(x0,@(xx) nlfunc(xx,Psi,Psit,b,ind,N),[tolf,0]);
-
-%fprintf('flag=%d, iter=%d\n',flag,iter);
-
 
 x(1:end)=0;
 x(ind)=xt;
@@ -108,25 +93,6 @@
 val=Psit(Psi(x2));
 val=val(ind);
 
-%val=sign(val).*(abs(val)).^(0.8);
-
-end
-
-% dd
-function val=nlfunc(x,Psi,Psit,b,ind,N)
-
-x2=zeros(N,1);
-x2(ind)=x;
-
-val=Psit(Psi(x2));
-val=val(ind)-b;
-
-%val=val(ind);
-%val=sign(val).*log(abs(val))-sign(b).*log(abs(b));
-val=sign(val).*(abs(val)).^(0.5);
-
-%val(abs(val)>1e9)=-1e9;
-
 end
 
 

functions/cs/spectrum_optim.m

@@ -31,22 +31,6 @@
 Ldist=exp(linspace(-d,d,Res));
 
 
-
-
-% for i=1:param.resf
-% 	Li=L;
-% 	Li(jvec)=Li(jvec).*exp(randn(length(jvec),1)*param.df);
-% 	Li=round(Li);
-% 	% QOMOMP
-% 	u=D(Li,tolf);
-% 	% spectrum
-% 	[~,Ea]=E(u);
-%
-% 	Eopt=Eopt+log(Ea);
-%
-% end
-
-
 tolf=tolf*10;
 
 Ljend=L(jvec(end));
@@ -85,17 +69,12 @@
 			[~,Ea]=E(At(A(u)));
 			% error in log space
 			flag = outp.flag ~= 0;  % 0 if 0, 1 if nonzero
-			%confact=((1-flag)+flag*outp.relres/tolf); % convergence factor
 			confact=1;  % do nothing
 			e(i)=norm(log(Ea(wj))-logspec)*confact;
-			%e(i)=norm(Ea(wj)-logspec)*confact;  % linear
-
+
 			flag(i)=outp.flag;
 
-			fprintf('i');
-
 		end
-		%plot(Lcoef,e,'-*'),hold all
 
 		if L0(j)<L0(j-1)
 			eind=Lcoef<=L(j-1)*dimsc;  % restrict L(j)<=L(j-1)*dimsc
@@ -105,36 +84,18 @@
 		imin=find(e(eind)==min(e(eind)),1,'first');
 		L(j)=Lcoef(imin);
 
-		fprintf(' j %d ,',L(j));
 		disp(flag(:)')
 	end
-	fprintf('\n');
-end
-%figure,plot(Lcoef,e,'-*')
 
-disp('final optim')
-Eopt=Eg*0;
-% average spectra
-for i=1:param.resf
-	Li=L;
-	%Li(jvec)=Li(jvec).*exp(randn(length(jvec),1)*param.df);  % randomize?
-	Li=round(Li);
-	% QOMOMP
-	%u=D(Li,tolf);
-	[u,outp]=D(Li,tolf);
-	disp(outp)
-	% spectrum
-	[~,Ea]=E(u);
-
-	Eopt=Eopt+log(Ea);
-
 end
 
-Eopt=exp(Eopt/param.resf);
-Lopt=round(L);
-
-
-
+Li=round(L);
+% QOMOMP
+[u,outp]=D(Li,tolf);
+disp(outp)
+% spectrum
+[~,Eopt]=E(u);
+Lopt=Li;
 
 
 

main_1d_other_run.m

@@ -48,7 +48,7 @@
 % =====================================
 
 % undersampling ratio
-do1vec=4;		% for SpESO (DECODER=82) use half:  82:(32,16,8), 1:(16,8,4)
+do1vec=8;		% for SpESO (DECODER=82) use half:  82:(32,16,8), 1:(16,8,4)
 SIMN=length(do1vec);
 SIMN0=4;		% number of experiments for each case
 

main_1d_speso_run.m

@@ -48,7 +48,7 @@
 % =====================================
 
 % undersampling ratio
-do1vec=8;		% for SpESO (DECODER=82) use half:  82:(32,16,8), 1:(16,8,4)
+do1vec=16;		% for SpESO (DECODER=82) use half:  82:(32,16,8), 1:(16,8,4)
 SIMN=length(do1vec);
 SIMN0=4;		% number of experiments for each case
 

plot_1d_results.m

@@ -4,7 +4,7 @@
 % set parameters according to those set in main files
 SAVEDLOC='./';
 N=1024*32;
-SIMN0=64;  %4
+SIMN0=4;  %4
 SIMN0t=num2str(SIMN0);
 mstr='synth';
 K=num2str(284);
@@ -14,7 +14,7 @@
 M2=num2str(2066); %speso
 
 % which simulation cases to plot?
-PLOTCASES=1:4;  % all: 1:SIMN0
+PLOTCASES=1:SIMN0;  % all: 1:SIMN0
 
 
 var1=load([SAVEDLOC,'saved_1d_',num2str(N),'_',M1,'_',SIMN0t,'_',mstr,'_',K,'_1_31_',EXTRA,'.mat']);

plot_2d_results.m

@@ -5,7 +5,7 @@
 SAVEDLOC='./';
 Ns=1024;
 N=Ns^2;
-SIMN0=8; %2
+SIMN0=2; %2
 SIMN0t=num2str(SIMN0);
 seedstr='1_';  % SEED
 mstr='synth';
@@ -15,7 +15,7 @@
 M2=num2str(32761); %speso
 
 % which simulation cases to plot?
-PLOTCASES=1:2;  % all: 1:SIMN0
+PLOTCASES=1:SIMN0;  % all: 1:SIMN0
 
 
 var1=load([SAVEDLOC,'saved_2d_',num2str(N),'_',M1,'_',SIMN0t,'_',mstr,'_1_2_31_',seedstr,EXTRA,'.mat']);