diff --git a/.gitignore b/.gitignore index 2ebc232..543ebcc 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,4 @@ *doctrees/ *Opti/OptiUtils/MatlabOptimPack/ .DS_Store +*Util/UnitTest/Data/ \ No newline at end of file diff --git a/Abstract/Cost.m b/Abstract/Cost.m index 8fa6e23..fbc5de6 100644 --- a/Abstract/Cost.m +++ b/Abstract/Cost.m @@ -98,7 +98,7 @@ % $$ \\mathrm{prox}_{\\alpha C}(\\mathrm{z}) = \\mathrm{arg} \\, \\mathrm{min}_{\\mathrm{u} \\in \\mathrm{X}} \\; \\frac{1}{2\\alpha} \\| \\mathrm{u} - \\mathrm{z} \\|_2^2 + C(\\mathrm{u}). $$ % % Calls the method :meth:`applyProx_` - assert(isscalar(alpha),'alpha must be a scalar'); + assert(isnumeric(alpha),'alpha must be a numeric'); if ~checkSize(z, this.sizein) % check input size error('Input to applyProx was size [%s], didn''t match stated sizein [%s].',... num2str(size(z)), num2str(this.sizein)); @@ -209,7 +209,6 @@ % Uses the method applyGrad (hence do not need to be % reimplemented in derived classes) x=y.*this.applyGrad(v); - x=sum(x(:)); end end diff --git a/Abstract/Map.m b/Abstract/Map.m index ea4db25..8e47ef1 100644 --- a/Abstract/Map.m +++ b/Abstract/Map.m @@ -159,59 +159,128 @@ % $$ \\mathrm{M}(\\mathrm{x}) := \\mathrm{H}(\\mathrm{x}) + \\mathrm{G}(\\mathrm{x})$$ % % Calls the method :meth:`plus_` - if ~cmpSize(this.sizein, G.sizein) % check input size - error('Input to plus is a %s of sizein [%s], which didn''t match the added %s sizein [%s].',... - class(G),num2str(G.sizein),class(this), num2str(this.sizein)); - end - if ~cmpSize(this.sizeout, G.sizeout) % check input size - error('Input to plus is a %s of sizeout [%s], which didn''t match the added %s sizeout [%s].',... - class(G),num2str(G.sizeout),class(this), num2str(this.sizeout)); + + if isnumeric(this) && isscalar(this) + if isequal(this,0) + M = G; + elseif cmpSize(G.sizein, G.sizeout) + M = LinOpDiag(G.sizein,this) +G; + else + error('Cannot sum a scalar and a non-square matrix'); + end + elseif isnumeric(G) && isscalar(G) + if isequal(G,0) + M = this; + elseif cmpSize(this.sizein, this.sizeout) + M = this + LinOpDiag(this.sizein,G); + else + error('Cannot sum a scalar and a non-square matrix'); + end + else + if ~cmpSize(this.sizein, G.sizein) % check input size + error('Input to plus is a %s of sizein [%s], which didn''t match the added %s sizein [%s].',... + class(G),num2str(G.sizein),class(this), num2str(this.sizein)); + end + if ~cmpSize(this.sizeout, G.sizeout) % check input size + error('Input to plus is a %s of sizeout [%s], which didn''t match the added %s sizeout [%s].',... + class(G),num2str(G.sizeout),class(this), num2str(this.sizeout)); + end + M = this.plus_(G); end - M = this.plus_(G); end function M = minus(this,G) % Overload operator (-) for :class:`Map` objects % $$ \\mathrm{M}(\\mathrm{x}) := \\mathrm{H}(\\mathrm{x}) - \\mathrm{G}(\\mathrm{x})$$ % % Calls the method :meth:`minus_` - if ~cmpSize(this.sizein, G.sizein) % check input size - error('Input to plus is a %s of sizein [%s], which didn''t match the substracted %s sizein [%s].',... - class(G),num2str(G.sizein),class(this), num2str(this.sizein)); - end - if ~cmpSize(this.sizeout, G.sizeout) % check input size - error('Input to plus is a %s of sizeout [%s], which didn''t match the substracted %s sizeout [%s].',... - class(G),num2str(G.sizeout),class(this), num2str(this.sizeout)); + if isnumeric(this) && isscalar(this) && isequal(this,0) % if multiply by 0 return 0 + M = -1*G; + elseif isnumeric(G) && isscalar(G) + if isequal(G,0) + M = this; + else + M = this - LinOpDiag(this.sizeout,G); + end + else + + + if ~cmpSize(this.sizein, G.sizein) % check input size + error('Input to plus is a %s of sizein [%s], which didn''t match the substracted %s sizein [%s].',... + class(G),num2str(G.sizein),class(this), num2str(this.sizein)); + end + if ~cmpSize(this.sizeout, G.sizeout) % check input size + error('Input to plus is a %s of sizeout [%s], which didn''t match the substracted %s sizeout [%s].',... + class(G),num2str(G.sizeout),class(this), num2str(this.sizeout)); + end + M = this.minus_(G); end - M = this.minus_(G); end function M = mpower(this,p) % Returns a new :class:`Map` which is the power p \\(\\mathrm{H}^{p}\\) of the % current \\(\\mathrm{H}\\). % % Calls the method :meth:`mpower_` - M=this.mpower_(p); + if ~cmpSize(this.sizeout, this.sizein) + error('Input to mpower must be a square map. Here %s x %s',... + num2str(this.sizeout),num2str(this.sizein)); + end + if isnumeric(p) && isscalar(p) + if isequal(p,0) + M = LinOpIdentity(G.sizein); + elseif isequal(p,1) + M = this; + else + M=this.mpower_(p); + end + else + error('p must be a scalar'); + end + end + + function M = power(this,p) + % Returns a new :class:`Map` which is the power p \\(\\mathrm{H}.^{p}\\) of the + % current \\(\\mathrm{H}\\). + % + % Calls the method :meth:`mpower_` + if isnumeric(p) && isscalar(p) + if isequal(p,1) + M = this; + else + M=this.power_(p); + end + else + error('p must be a scalar'); + end end function G = mtimes(this,G) % Overload operator (*) for :class:`Map` objects % $$ \\mathrm{M}(\\mathrm{x}) := \\mathrm{H}(\\mathrm{G}(\\mathrm{x}))$$ - % - If \\(\\mathrm{G}\\) is numeric of size sizein, then :meth:`apply` is called - % - If \\(\\mathrm{G}\\) is a :class:`Map`, then a - % :class:`MapComposition`is intanciated - if isa(G,'Map') - if (isnumeric(this) && isscalar(this)) % Left multiplication by scalar - if ~isequal(this,1) % if multiply by 1 do noting - if isa(G,'Cost') - G=CostMultiplication(this,G); - else - H=LinOpDiag(G.sizeout,this); - G=H.makeComposition(G); + % + % - If \\(\\mathrm{G}\\) is numeric of size sizein, then :meth:`apply` is called + % + % - If \\(\\mathrm{G}\\) is a :class:`Map`, then a :class:`MapComposition` is intanciated + + if isa(this,'LinOp') && isnumeric(G) && isscalar(G)&& isequal(G,0) + G = 0; + else + if isa(G,'Map') + if (isnumeric(this) && isscalar(this)) % Left multiplication by scalar + if isequal(this,0) % if multiply by 0 return 0 + G = 0; + elseif ~isequal(this,1) % if multiply by 1 do noting + if isa(G,'Cost') + G=CostMultiplication(this,G); + else + H=LinOpDiag(G.sizeout,this); + G=H.makeComposition(G); + end end + else + G =this.makeComposition(G); end else - G =this.makeComposition(G); + G = this.apply(G); end - else - G = this.apply(G); end end function M = times(this,G) @@ -220,17 +289,23 @@ % $$ \\mathrm{M}(\\mathrm{x}) := \\mathrm{H}(\\mathrm{x}) \\times \\mathrm{G}(\\mathrm{x})$$ % % Calls the method :meth:`times_` - if ~cmpSize(this.sizein, G.sizein) % check input size - error('Input to times is a %s of sizein [%s], which didn''t match the multiplied %s sizein [%s].',... - class(G),num2str(G.sizein),class(this), num2str(this.sizein)); - end - if ~cmpSize(this.sizeout, G.sizeout) % check input size - error('Input to times is a %s of sizeout [%s], which didn''t match the multiplied %s sizeout [%s].',... - class(G),num2str(G.sizeout),class(this), num2str(this.sizeout)); + + if isnumeric(this) && isscalar(this) && isequal(this,0) % if multiply by 0 return 0 + M = 0; + elseif isa(this,'LinOp')&& isnumeric(G)&& isscalar(G) && isequal(G,0) + M = 0; + else + if ~cmpSize(this.sizein, G.sizein) % check input size + error('Input to times is a %s of sizein [%s], which didn''t match the multiplied %s sizein [%s].',... + class(G),num2str(G.sizein),class(this), num2str(this.sizein)); + end + if ~cmpSize(this.sizeout, G.sizeout) % check input size + error('Input to times is a %s of sizeout [%s], which didn''t match the multiplied %s sizeout [%s].',... + class(G),num2str(G.sizeout),class(this), num2str(this.sizeout)); + end + M=this.times_(G); end - M=this.times_(G); end - % TODO Overload operator .^ to do (H(x)).^p ? function sz = size(this, varargin) sz = {this.sizeout, this.sizein}; if length(varargin) == 1 @@ -280,10 +355,30 @@ % When \\(p\\neq-1\\), this method is not implemented in this Abstract class if p==-1 M=this.makeInversion_(); + elseif mod(p,1)==0 + if p>0 + M = this; + while p>1 + M = M *this; + p = p-1; + end + else + M = this; + p = abs(p); + while p>1 + M = M *this; + p = p-1; + end + M= M.makeInversion_(); + end else error('mpower_ method not implemented for the given power==%d',p); end end + function M = power_(~,p) + % this method is not implemented in this Abstract class + error('power_ method not implemented for the given power==%d',p); + end function M = times_(this,G) % Constructs a :class:`MapMultiplication` object to element-wise multiply the % current :class:`Map` \\(\\mathrm{H}\\) with the given \\(\\mathrm{G}\\). diff --git a/Abstract/OperationsOnMaps/CostComposition.m b/Abstract/OperationsOnMaps/CostComposition.m index 795b28a..c55838f 100644 --- a/Abstract/OperationsOnMaps/CostComposition.m +++ b/Abstract/OperationsOnMaps/CostComposition.m @@ -77,8 +77,9 @@ % If this.H2 is a :class:`LinOp` and \\(\\mathrm{H} % \\mathrm{H}^{\\star}\\) is a :class:`LinOpScaledIdentity` % - if this.isConvex && this.isH2LinOp && this.isH2SemiOrtho - x = z + 1/this.nu*this.H2.applyAdjoint(this.H1.applyProx(this.H2*z,alpha*this.nu)-this.H2*z); + if this.isConvex && this.isH2LinOp && this.isH2SemiOrtho + H2z=this.H2*z; + x = z + 1/this.nu*this.H2.applyAdjoint(this.H1.applyProx(H2z,alpha*this.nu)-H2z); else x = applyProx_@Cost(this,z,alpha); end diff --git a/Abstract/OperationsOnMaps/CostPartialSummation.m b/Abstract/OperationsOnMaps/CostPartialSummation.m index dba0e4d..d814908 100644 --- a/Abstract/OperationsOnMaps/CostPartialSummation.m +++ b/Abstract/OperationsOnMaps/CostPartialSummation.m @@ -56,6 +56,8 @@ function setLsub(this,Lsub) methods (Access = protected) function updateSubset(this) switch this.partialGrad + case 0 + this.subset = 1:this.numMaps; case 1 this.subset = sort(randi(this.numMaps,this.Lsub,1)); case 2 @@ -86,7 +88,7 @@ function updateSubset(this) function g=applyGrad_(this,x) % Reimplemented from :class:`Cost` if this.partialGrad > 0 - g = zeros(size(x)); + g = zeros_(size(x)); for kk = 1:this.Lsub ind = this.subset(kk); g = g + this.alpha(ind)*this.mapsCell{ind}.applyGrad(x); diff --git a/Abstract/OperationsOnMaps/MapComposition.m b/Abstract/OperationsOnMaps/MapComposition.m index 83aacbe..30f89fe 100644 --- a/Abstract/OperationsOnMaps/MapComposition.m +++ b/Abstract/OperationsOnMaps/MapComposition.m @@ -82,7 +82,7 @@ function M = makeComposition_(this,G) % Reimplemented from :class:`Map` H2G = this.H2*G; % Since H1*H2 is not simplified, first try to compose H2 with G - if ~isa(H2G, 'MapComposition'); + if ~isa(H2G, 'MapComposition') M=this.H1*H2G; else M = MapComposition(this, G); diff --git a/Abstract/OperationsOnMaps/MapSummation.m b/Abstract/OperationsOnMaps/MapSummation.m index 5f75a94..306c602 100644 --- a/Abstract/OperationsOnMaps/MapSummation.m +++ b/Abstract/OperationsOnMaps/MapSummation.m @@ -69,6 +69,15 @@ this.isInvertible=false; this.sizein = this.mapsCell{1}.sizein; this.sizeout = this.mapsCell{1}.sizeout; + this.norm=0; + for n=1:this.numMaps + if this.mapsCell{n}.norm~=-1 + this.norm=this.norm+abs(this.alpha(n))*this.mapsCell{n}.norm; + else + this.norm=-1; + break + end + end for n =2:this.numMaps assert(cmpSize(this.sizein,this.mapsCell{n}.sizein),'%d-th input does not have consistent sizein', n) ; assert(cmpSize(this.sizeout,this.mapsCell{n}.sizeout),'%d-th input does not have the consistent sizeout ', n); @@ -91,7 +100,7 @@ end function x = applyJacobianT_(this, y, v) % Reimplemented from :class:`Map` - x = zeros_(this.sizeout); + x = zeros_(this.sizein); for n = 1:this.numMaps x = x + this.alpha(n) .* this.mapsCell{n}.applyJacobianT(y,v); end diff --git a/Abstract/Opti.m b/Abstract/Opti.m index 6016733..ccd458c 100644 --- a/Abstract/Opti.m +++ b/Abstract/Opti.m @@ -80,11 +80,11 @@ function run(this,x0) % - Algorithm iteration flag=this.doIteration(); - if(flag == this.OPTI_NEXT_IT) + if(flag == this.OPTI_NEXT_IT ) this.niter=this.niter+1; % - Convergence test - if this.CvOp.testConvergence(this), break; end + if this.niter>1 && this.CvOp.testConvergence(this), break; end % - Call OutputOpti object if this.ItUpOut>0 && (((mod(this.niter,this.ItUpOut)==0)|| (this.niter==1))) diff --git a/Cost/CostGoodRoughness.m b/Cost/CostGoodRoughness.m index 222cdb4..a1a1d59 100644 --- a/Cost/CostGoodRoughness.m +++ b/Cost/CostGoodRoughness.m @@ -4,7 +4,7 @@ % with \\( \\vert (\\nabla \\mathrm{x})_{.,k} \\vert^2 \\) the gradient % magnitude at pixel k. % - % :param G: :class:`LinOpGrad`object + % :param G: :class:`LinOpGrad` object % :param bet: smoothing parameter (default 1e-1) % % All attributes of parent class :class:`Cost` are inherited. diff --git a/Cost/CostHyperBolic.m b/Cost/CostHyperBolic.m index c3080a5..3c25f40 100644 --- a/Cost/CostHyperBolic.m +++ b/Cost/CostHyperBolic.m @@ -1,9 +1,9 @@ classdef CostHyperBolic < Cost % CostHyperBolic: Hyperbolic cost function - % $$C(\\mathrm{x}) := \\sum_{k=1}^K \\sqrt{\\sum_{l=1}^L (\\mathrm{x}-y)_{k,l}^2 + \\varepsilon^2}$$ + % $$C(\\mathrm{x}) := \\sum_{k=1}^K \\sqrt{\\sum_{l=1}^L (\\mathrm{x}-y)_{k,l}^2 + \\varepsilon_k^2} - \\sum_{k=1}^K\\varepsilon_k $$ % % :param index: dimensions along which the l2-norm will be applied (inner sum over l) - % :param epsilon: \\(\\in \\mathbb{R}_+\\) smoothing parameter (default + % :param epsilon: \\(\\in \\mathbb{R}^K_+\\) smoothing parameter (default % \\(10^{-3}\\)) % % All attributes of parent class :class:`Cost` are inherited. @@ -33,6 +33,7 @@ epsilon; index; sumOp; + sumEpsilon; end %% Constructor @@ -47,18 +48,29 @@ if nargin<3|| isempty(index) %no sum index=0; end - if nargin<2|| isempty(epsilon) - epsilon=1e-3; - end - this.epsilon = epsilon; this.index = index; - if index~=0 this.sumOp = LinOpSum(sz,index); else this.sumOp = LinOpDiag(sz); end + + if nargin<2|| isempty(epsilon) + epsilon=1e-3; + end + if ~(isscalar(epsilon)||checkSize(epsilon, this.sumOp.sizeout)) + error('epsilon is of size [%s] and didn''t match [%s].',... + num2str(size(epsilon)), num2str(this.sumOp.sizeout)); + end + + this.epsilon = epsilon; + if isscalar(this.epsilon) + this.sumEpsilon = prod(this.sumOp.sizeout(:)).*this.epsilon; + else + this.sumEpsilon = sum(this.epsilon(:)); + end + this.memoizeOpts.computeF=true; this.memoCache.computeF= struct('in',[],'out', []); end diff --git a/Cost/CostL2.m b/Cost/CostL2.m index 44a792f..c9aa905 100644 --- a/Cost/CostL2.m +++ b/Cost/CostL2.m @@ -18,6 +18,7 @@ % This program is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation, either version 3 of the License, or + % % (at your option) any later version. % % This program is distributed in the hope that it will be useful, @@ -65,10 +66,10 @@ % Reimplemented from parent class :class:`Cost`. if (isscalar(this.y)&&(this.y==0)) wr=this.W*(x); - y=0.5*dot(x(:),wr(:)); + y=0.5*real(dot(x(:),wr(:))); else wr=this.W*(x-this.y); - y=0.5*dot(x(:)-this.y(:),wr(:)); + y=0.5*real(dot(x(:)-this.y(:),wr(:))); end end function g=applyGrad_(this,x) diff --git a/Cost/CostL2Composition.m b/Cost/CostL2Composition.m index c489ba2..f6b0118 100644 --- a/Cost/CostL2Composition.m +++ b/Cost/CostL2Composition.m @@ -95,7 +95,7 @@ if ~isfield(this.precomputeCache,'ytWy') this.precomputeCache.ytWy= this.H1.y(:)' * reshape(this.H1.W * this.H1.y,numel(this.H1.y), 1); end - y=0.5 * x(:)' * reshape(this.H1.W*this.H2.applyHtH(x),numel(this.H1.y),1) - x(:)' * this.precomputeCache.WHty(:) + 0.5 * this.precomputeCache.ytWy; + y=0.5 * x(:)' * reshape(this.H1.W*this.H2.applyHtH(x),prod(this.H2.sizein),1) - x(:)' * this.precomputeCache.WHty(:) + 0.5 * this.precomputeCache.ytWy; y = real(y); % due to numerical error else y=apply_@CostComposition(this,x); @@ -121,6 +121,7 @@ end function y=applyProx_(this,x,alpha) % Reimplemented from parent class :class:`CostComposition`. + % % Implemented % - if the operator \\(\\alpha\\mathrm{H^{\\star}WH + I} \\) is invertible: % $$ \\mathrm{y} = (\\alpha\\mathrm{H^{\\star}WH + I} )^{-1} (\\alpha \\mathrm{H^TWy +x})$$ diff --git a/Cost/CostMixNormSchatt1.m b/Cost/CostMixNormSchatt1.m index 05624ae..c497e32 100644 --- a/Cost/CostMixNormSchatt1.m +++ b/Cost/CostMixNormSchatt1.m @@ -11,18 +11,17 @@ % % All attributes of parent class :class:`Cost` are inherited. % - % **Note** The actual implementation works for size (sz) having one of the two following forms: + % **Note** The actual implementation works for size (sz) having one of the two following forms: % - % - (NxMx3) such that the Sp norm will be applied on each symetric 2x2 - % $$ \\begin{bmatrix} \\mathrm{x}_{n m 1} & \\mathrm{x}_{n m 2} \\newline - % \\mathrm{x}_{n m 2} & \\mathrm{x}_{n m 3} \\end{bmatrix}$$ - % and then the \\(\\ell_1\\) norm on the two other dimensions. - % - % - (NxMxKx6) such that the Sp norm will be applied on each symetric 3x3 - % $$ \\begin{bmatrix} \\mathrm{x}_{n m k 1} & \\mathrm{x}_{n m k 2} & \\mathrm{x}_{n m k 3} \\newline - % \\mathrm{x}_{n m k 2} & \\mathrm{x}_{n m k 4} & \\mathrm{x}_{n m k 5} \\newline - % \\mathrm{x}_{n m k 3} & \\mathrm{x}_{n m k 5} & \\mathrm{x}_{n m k 6} \\newline \\end{bmatrix}$$ - % and then the \\(\\ell_1\\) norm on the three other dimensions. + % * (NxMx3) such that the Sp norm will be applied on each symetric 2x2 + % $$ \\begin{bmatrix} \\mathrm{x}_{n m 1} & \\mathrm{x}_{n m 2} \\newline + % \\mathrm{x}_{n m 2} & \\mathrm{x}_{n m 3} \\end{bmatrix}$$ + % and then the \\(\\ell_1\\) norm on the two other dimensions. + % * (NxMxKx6) such that the Sp norm will be applied on each symetric 3x3 + % $$ \\begin{bmatrix} \\mathrm{x}_{n m k 1} & \\mathrm{x}_{n m k 2} & \\mathrm{x}_{n m k 3} \\newline + % \\mathrm{x}_{n m k 2} & \\mathrm{x}_{n m k 4} & \\mathrm{x}_{n m k 5} \\newline + % \\mathrm{x}_{n m k 3} & \\mathrm{x}_{n m k 5} & \\mathrm{x}_{n m k 6} \\newline \\end{bmatrix}$$ + % and then the \\(\\ell_1\\) norm on the three other dimensions. % % **References** % [1] Lefkimmiatis, S., Ward, J. P., & Unser, M. (2013). Hessian Schatten-norm regularization @@ -97,25 +96,14 @@ function y=applyProx_(this,x,alpha) % Reimplemented from parent class :class:`Cost`. - dim=size(x); if this.p==1 [E,V]=this.svdDecomp(x); E=max(abs(E)-alpha,0).*sign(E); y=reshape(this.svdRecomp(E,V),this.sizein); - % y=this.svdRecomp(E,V); elseif this.p==2 - % Emmanuel (05/07/2018): this part should also be updated - % in order to use the last dimension (not :,:,k or :,:,:,k, - % but more general - if dim(end)==3 % 2D - N=sqrt(x(:,:,1).^2+2*x(:,:,2).^2+x(:,:,3).^2); - y=repmat((N-1)./N,[ones(1,length(this.sizein)-1),3]).*x; - elseif dim(end)==6 % 3D - N=sqrt(x(:,:,:,1).^2+2*x(:,:,:,2).^2+2*x(:,:,:,3).^2+x(:,:,:,4)+2*x(:,:,:,5)+x(:,:,:,6)); - y=repmat((N-1)./N,[ones(1,length(this.sizein)-1),6]).*x; - else - error('last dimension of x should be 3 or 6'); - end + [E,V]=this.svdDecomp(x); + E=E./(1+alpha); + y=reshape(this.svdRecomp(E,V),this.sizein); else y=applyProx_@Cost(this,x,alpha); end diff --git a/Cost/CostUtils/HessianSchatten/buildHessianSchatten.m b/Cost/CostUtils/HessianSchatten/buildHessianSchatten.m index f501b46..a7e026d 100644 --- a/Cost/CostUtils/HessianSchatten/buildHessianSchatten.m +++ b/Cost/CostUtils/HessianSchatten/buildHessianSchatten.m @@ -1,6 +1,9 @@ -function buildHessianSchatten() +function buildHessianSchatten(options) %% buildHessianSchatten function % build the HessianSchatten norm mexgl files for CostMixNormSchatt1 +% +% You can give as a parameter of this function the path to your GCC +% compiler. Ex: buildHessianSchatten('GCC=/usr/bin/gcc-6') % Copyright (C) 2018 F. Soulez ferreol.soulez@epfl.ch % @@ -16,16 +19,26 @@ function buildHessianSchatten() % % You should have received a copy of the GNU General Public License % along with this program. If not, see . +if nargin==0 + options=[]; +end -[mpath,~,~] = fileparts(which('buildHessianSchatten')); disp('Installing HessianSchatten'); -disp('WARNING: depending on your system and compiler, OPENMP may not be activated leading to slow computation.'); +get_architecture; +if linux + options = [ options, ' CXXFLAGS='' -fopenmp ''',' LDFLAGS=''$LDFLAGS -fopenmp ''']; +else + disp('On your system and compiler, OPENMP is desactivated leading to slow computation. This can be tuned using the options parameter:'); + disp('Example: options = CXXFLAGS= -fopenmp '); +end + +[mpath,~,~] = fileparts(which('buildHessianSchatten')); pth = cd; cd(mpath); -MexOpt= ['-DUSE_BLAS_LIB ' '-DNEW_MATLAB_BLAS ' '-DINT_64BITS ' '-largeArrayDims ' 'COMPFLAGS=''$COMPFLAGS -Wall -mtune=native -fomit-frame-pointer -O2 -fopenmp ''']; -eval(['mex ','svd2D_recomp.cpp ',MexOpt]); -eval(['mex ','svd2D_decomp.cpp ',MexOpt]); -eval(['mex ','svd3D_recomp.cpp ',MexOpt]); -eval(['mex ','svd3D_decomp.cpp ',MexOpt]); +MexOpt= ['-DUSE_BLAS_LIB ' '-DNEW_MATLAB_BLAS ' '-DINT_64BITS ' '-largeArrayDims ' ,options, ' CXXFLAGS=''$CXXFLAGS -fPIC -Wall -mtune=native -fomit-frame-pointer -O2 ''' ' LDFLAGS=''$LDFLAGS ''']; +eval(['mex ',' svd2D_recomp.cpp ',MexOpt]); +eval(['mex ',' svd2D_decomp.cpp ',MexOpt]); +eval(['mex ',' svd3D_recomp.cpp ',MexOpt]); +eval(['mex ',' svd3D_decomp.cpp ',MexOpt]); cd(pth); end \ No newline at end of file diff --git a/Cost/CostUtils/HessianSchatten/svd2D_decomp.cpp b/Cost/CostUtils/HessianSchatten/svd2D_decomp.cpp index bcafd76..45c138b 100644 --- a/Cost/CostUtils/HessianSchatten/svd2D_decomp.cpp +++ b/Cost/CostUtils/HessianSchatten/svd2D_decomp.cpp @@ -20,6 +20,8 @@ -linux: mex -v svd2D_decomp.cpp CFLAGS="\$CFLAGS -openmp" LDFLAGS="\$LDFLAGS -openmp" -largeArrayDims -mac : mex svd2D_decomp.cpp -DUSE_BLAS_LIB -DNEW_MATLAB_BLAS -DINT_64BITS -largeArrayDims CXX=/usr/local/Cellar/gcc/6.3.0_1/bin/g++-6 CXXOPTIMFLAGS="-O3 -mtune=native -fomit-frame-pointer -fopenmp" LDOPTIMFLAGS=" -O " LINKLIBS="$LINKLIBS -lmwblas -lmwlapack -L"/usr/local/Cellar/gcc/6.3.0_1/lib/gcc/6" -L/ -fopenmp" + -mac (another option): "brew install llvm" on terminal, then use the command + mex svd2D_decomp.cpp -DUSE_BLAS_LIB -DNEW_MATLAB_BLAS -DINT_64BITS -largeArrayDims CXX="/usr/local/opt/llvm/bin/clang++" Copyright (C) 2017 E. Soubies emmanuel.soubies@epfl.ch diff --git a/Cost/IndicatorFunctions/CostComplexCircle.m b/Cost/IndicatorFunctions/CostComplexCircle.m index ee5e964..a7177ca 100644 --- a/Cost/IndicatorFunctions/CostComplexCircle.m +++ b/Cost/IndicatorFunctions/CostComplexCircle.m @@ -7,7 +7,7 @@ % All attributes of parent class :class:`CostComplexRing` are inherited % % - % :param radius radius of the circle (default 1) + % :param radius: radius of the circle (default 1) % % **Example** C=CostComplexCircle(radius,sz,y) % diff --git a/Cost/IndicatorFunctions/CostComplexDisk.m b/Cost/IndicatorFunctions/CostComplexDisk.m index ede6f7b..0f39ba5 100644 --- a/Cost/IndicatorFunctions/CostComplexDisk.m +++ b/Cost/IndicatorFunctions/CostComplexDisk.m @@ -6,7 +6,7 @@ % % All attributes of parent class :class:`CostComplexRing` are inherited % - % :param radius radius of the disk (default 1) + % :param radius: radius of the disk (default 1) % % **Example** C=CostComplexDisk(sz,radius, y) % diff --git a/Doc/source/conditionsuse.rst b/Doc/source/conditionsuse.rst index b9b89ee..8c701e3 100644 --- a/Doc/source/conditionsuse.rst +++ b/Doc/source/conditionsuse.rst @@ -10,7 +10,10 @@ If not, see . Whenever you present or publish results that are based on this library, please cite: -M. Unser, E. Soubies, F. Soulez, M. McCann, L. Donati, -`GlobalBioIm: A Unifying Computational Framework for Solving Inverse Problems `_ -Proceedings of the OSA Imaging and Applied Optics Congress on Computational Optical Sensing and Imaging (COSI'17), San Francisco CA, USA, June 26-29, 2017, paper no. CTu1B. +| `Pocket Guide to Solve Inverse Problems with GlobalBioIm `_, +| Submitted, 2019. +| E. Soubies, F. Soulez, M. T. McCann, T-A. Pham, L. Donati, T. Debarre, D. Sage, and M. Unser. +| `GlobalBioIm: A Unifying Computational Framework for Solving Inverse Problems `_, +| Proceedings of the OSA Computational Optical Sensing and Imaging congres, 2017. +| M. Unser, E. Soubies, F. Soulez, M. McCann, and L. Donati. \ No newline at end of file diff --git a/Doc/source/conf.py b/Doc/source/conf.py index 6987e04..66cdd09 100644 --- a/Doc/source/conf.py +++ b/Doc/source/conf.py @@ -62,9 +62,9 @@ # built documents. # # The short X.Y version. -version = u'1.1.1' +version = u'1.1.2' # The full version, including alpha/beta/rc tags. -release = u'1.1.1' +release = u'1.1.2' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/Doc/source/contact.rst b/Doc/source/contact.rst index fb35e60..e4294ae 100644 --- a/Doc/source/contact.rst +++ b/Doc/source/contact.rst @@ -2,7 +2,7 @@ Contact ******* -For any question concerning the GlobalBioIm library, :email:`contact us `. +For any question concerning the GlobalBioIm library, :email:`contact us `. diff --git a/Doc/source/cost.rst b/Doc/source/cost.rst index 8b19248..6c7f9e6 100644 --- a/Doc/source/cost.rst +++ b/Doc/source/cost.rst @@ -16,7 +16,7 @@ CostGoodRoughness applyGrad_, applyProx_, applyProxFench_, set_y CostHyperBolic ---------------------- +-------------- .. autoclass:: CostHyperBolic :show-inheritance: @@ -68,7 +68,7 @@ CostMixNorm21 applyGrad_, applyProx_, applyProxFench_, set_y CostMixNorm21NonNeg -------------- +------------------- .. autoclass:: CostMixNorm21NonNeg :show-inheritance: diff --git a/Doc/source/examples.rst b/Doc/source/examples.rst index 53fc456..0d3ec12 100644 --- a/Doc/source/examples.rst +++ b/Doc/source/examples.rst @@ -76,11 +76,11 @@ constraint, respectively. We are now ready to instanciate and run our algorithm Hn={H,Opreg,Id}; % Associated operators H_n rho_n=[1e-3,1e-3,1e-3]; % Multipliers rho_n ADMM=OptiADMM([],Fn,Hn,rho_n); % Declare optimizer - ADMM.OutOp=OutputOpti(1,im,round(maxIt/10),[1 2]); % build the output object + ADMM.OutOp=OutputOpti(1,[],round(maxIt/10),[1 2]); % build the output object ADMM.CvOp=TestCvgCombine(TestCvgCostRelative(1e-4), 'StepRelative',1e-4); % Set the convergence tests ADMM.ItUpOut=dispIt; ADMM.maxiter=maxIt; - ADMM.run(xopt); + ADMM.run(zeros(sznew)); Here, three splittings have been done: \\(\\mathrm{u_1=Hx}, \\; \\mathrm{u_2=\\nabla x}\\), and \\(\\mathrm{u_3=x}\\). We do not need to provide a solver to the ADMM algorithm (5th argument) since the operator algebra ensures that the operator @@ -96,3 +96,4 @@ The deconvolved image is shown in Figure 2. :align: center Fig 2. Deconvolved C. Elegans embryo. + diff --git a/Doc/source/index.rst b/Doc/source/index.rst index 1cff590..7a7cc26 100644 --- a/Doc/source/index.rst +++ b/Doc/source/index.rst @@ -25,19 +25,24 @@ automatically from comments within M-files. Releases -------- - - `v 1.1.1 `_ (September 2018). + - `v 1.1.2 `_ (April 2018). + - `v 1.1.1 `_ (September 2018). - `v 1.1 `_ (July 2018). *Speed up your codes using the library with GPU* (:ref:`read more `). - `v 1.0.1 `_ (May 2018). - `v 1.0 `_ milestone (March 2018). - `v 0.2 `_ (November 2017). *New tools, more flexibility, improved composition*. - - v 0.1 (June 2017). First public release of the library. + - v 0.1 (June 2017). *First public release of the library.* -Reference ---------- +References +---------- - - M. Unser, E. Soubies, F. Soulez, M. McCann, L. Donati, - `GlobalBioIm: A Unifying Computational Framework for Solving Inverse Problems `_, - Proceedings of the OSA Imaging and Applied Optics Congress on Computational Optical Sensing and Imaging (COSI'17), San Francisco CA, USA, June 26-29, 2017, paper no. CTu1B. +| `Pocket Guide to Solve Inverse Problems with GlobalBioIm `_, +| Submitted, 2019. +| E. Soubies, F. Soulez, M. T. McCann, T-A. Pham, L. Donati, T. Debarre, D. Sage, and M. Unser. + +| `GlobalBioIm: A Unifying Computational Framework for Solving Inverse Problems `_, +| Proceedings of the OSA Computational Optical Sensing and Imaging congres, 2017. +| M. Unser, E. Soubies, F. Soulez, M. McCann, and L. Donati. Contents ******** @@ -46,9 +51,10 @@ Contents :maxdepth: 1 :caption: General - Download or Clone (v 1.1.1) + Download or Clone (v 1.1.2) infos - examples + examples + relatedPapers conditionsuse .. toctree:: diff --git a/Doc/source/linop.rst b/Doc/source/linop.rst index 401b1e1..57c396a 100644 --- a/Doc/source/linop.rst +++ b/Doc/source/linop.rst @@ -88,6 +88,15 @@ LinOpSum :members: apply_, applyJacobianT_, applyInverse_, plus_, minus_, mpower_, makeComposition_, applyAdjoint_, applyHtH_, applyHHt_, applyAdjointInverse_, makeAdjoint_, makeHtH_, makeHHt_, makeInversion_ + +LinOpXRay +--------- + +.. autoclass:: LinOpXRay + :show-inheritance: + :members: apply_, applyJacobianT_, applyInverse_, plus_, minus_, mpower_, makeComposition_, + applyAdjoint_, applyHtH_, applyHHt_, applyAdjointInverse_, makeAdjoint_, makeHtH_, makeHHt_, makeInversion_ + .. automodule:: LinOp.SelectorLinOps diff --git a/Doc/source/methodssummary.rst b/Doc/source/methodssummary.rst index eeacc46..764d37f 100644 --- a/Doc/source/methodssummary.rst +++ b/Doc/source/methodssummary.rst @@ -129,9 +129,7 @@ Opti +-------------------------+------------------------------------------------------+ | updateParams() | | Update algorithm parameters (e.g. descent step). | +-------------------------+------------------------------------------------------+ -| test_convergence() | | Test if the algorithm has converged. | -+-------------------------+------------------------------------------------------+ | starting_verb() | | Display starting message. | +-------------------------+------------------------------------------------------+ | ending_verb() | | Display ending message. | -+-------------------------+------------------------------------------------------+ \ No newline at end of file ++-------------------------+------------------------------------------------------+ diff --git a/Doc/source/nonlinop.rst b/Doc/source/nonlinop.rst index 1a21a4b..481c853 100644 --- a/Doc/source/nonlinop.rst +++ b/Doc/source/nonlinop.rst @@ -3,17 +3,17 @@ Non-Linear Operators This section contains non-linear operator classes which all derive directly from the abstract class :class:`Map`. -.. automodule:: NonLinop +.. automodule:: NonLinOp -.. automodule:: NonLinop.ElementWiseOp +.. automodule:: NonLinOp.ElementWiseOp Element-wise Operators ---------------------- -OpEWSquareRoot -.............. +OpEWSqrt +........ -.. autoclass:: OpEWSquareRoot +.. autoclass:: OpEWSqrt :show-inheritance: :members: apply_, applyJacobianT_, applyInverse_, plus_, minus_, mpower_, makeComposition_, diff --git a/Doc/source/relatedPapers.rst b/Doc/source/relatedPapers.rst new file mode 100644 index 0000000..447ca00 --- /dev/null +++ b/Doc/source/relatedPapers.rst @@ -0,0 +1,86 @@ +.. _ref-relatedPapers: + +Related Papers +************** + +This list gathers works that use the GlobalBioIm library. +If you want to add your work in this list, do not hesitate to +:email:`contact us`. + +Papers with open-source code +---------------------------- + +| **Improving 3D MA-TIRF Reconstruction with Deconvolution and Background Estimation** +| Proceedings of ISBI, 2019. + `[Paper] `_ + `[Code] `_ +| E. Soubies, L. Blanc-Feraud, S. Schaub, and E. Van Obberghen-Schilling. + +| **Nanometric Axial Resolution of Fibronectin Assembly Units Achieved with an Efficient Reconstruction Approach for Multi-Angle-TIRF Microscopy** +| Scientific Reports, 2019. + `[Paper] `_ + `[Code] `_ +| E. Soubies, A. Radwanska, D. Grall, L. Blanc-Feraud, E. Van Obberghen-Schilling, and S. Schaub. + +| **Region of interest X-ray computed tomography via corrected back projection** +| Proceedings of ISBI, 2018. + `[Paper] `_ + `[Code] `_ +| M. McCann, L. Vilaclara, and M. Unser. + +| **Efficient Inversion of Multiple-Scattering Model for Optical Diffraction Tomography** +| Optics Express 25-18, pp. 21786-21800, 2017. + `[Paper] `_ + `[Code] `_ +| E. Soubies, T.-A. Pham, and M. Unser. + + +Other papers +------------ + +| **Inner-Loop-Free ADMM for Cryo-EM** +| Proceedings of ISBI, 2019. + `[Paper] `_ +| L. Donati, E. Soubies, and M. Unser. + +| **Computational Super-Sectioning for Single-Slice Structured Illumination Microscopy** +| IEEE Transactions on Computational Imaging, 2018. + `[Paper] `_ +| E. Soubies and M. Unser. + +| **Phaseless Diffraction Tomography with Regularized Beam Propagation** +| Proceedings of ISBI, 2018. + `[Paper] `_ +| T.-A. Pham, E. Soubies, J. Lim, A. Goy, F. Soulez, D. Psaltis and M. Unser. + +| **Imaging neural activity in the ventral nerve cord of behaving adult Drosophila** +| Nature communications, 9-1, pp. 4390, 2018. + `[Paper] `_ +| C.-L. Chen, L. Hermans, M. Viswanathan, D. Fortun, F. Aymanns, M. Unser, A. Cammarato, M. Dickinson, and P. Ramdya. + +| **Versatile Reconstruction Framework for Diffraction Tomography with Intensity Measurements and Multiple Scattering** +| Optics Express, 26-3, pp. 2749-2763, 2018. + `[Paper] `_ +| T.-A. Pham, E. Soubies, J. Lim, A. Goy, F. Soulez, D. Psaltis and M. Unser. + +| **Deep Convolutional Neural Network for Inverse Problems in Imaging** +| IEEE Transactions on Image Processing, 2017. + `[Paper] `_ +| K.H. Jin, M.T. McCann, E. Froustey, and M. Unser. + +| **High-Quality Parallel-Ray X-Ray CT Back Projection Using Optimized Interpolation** +| IEEE Transactions on Image Processing, 2017. + `[Paper] `_ +| M.T. McCann and M. Unser. + +| **Compact lensless phase imager** +| Optics Express, Optical Society of America, 2017, 25 (4), pp.889 - 895. + `[Paper] `_ +| M. Rostykus, F. Soulez, M. Unser, C. Moser. + +| **Compact in-line lensfree digital holographic microscope** +| Methods, Elsevier, 2017, + `[Paper] `_ +| M. Rostykus, F. Soulez, M. Unser, C. Moser. + + diff --git a/Example/DeconvEx/2D/Deconv_KL_HessSchatt_NonNeg.m b/Example/DeconvEx/2D/Deconv_KL_HessSchatt_NonNeg.m index 24219d2..2b7d30e 100644 --- a/Example/DeconvEx/2D/Deconv_KL_HessSchatt_NonNeg.m +++ b/Example/DeconvEx/2D/Deconv_KL_HessSchatt_NonNeg.m @@ -53,7 +53,7 @@ Hn={H,Hess,LinOpDiag(sz)}; rho_n=[1e-2,1e-2,1e-2]; ADMM=OptiADMM([],Fn,Hn,rho_n,[]); -ADMM.OutOp=OutputOpti(1,im,40,[1 2]); +ADMM.OutOp=OutputOptiSNR(1,im,40,[1 2]); ADMM.ItUpOut=5; % call OutputOpti update every ItUpOut iterations ADMM.maxiter=200; % max number of iterations ADMM.run(y); % run the algorithm @@ -63,7 +63,7 @@ Fn={lamb*R_1sch,F}; Hn={Hess,H}; PDC=OptiPrimalDualCondat([],R_POS,Fn,Hn); -PDC.OutOp=OutputOpti(1,im,40,[2 3]); +PDC.OutOp=OutputOptiSNR(1,im,40,[2 3]); PDC.tau=100; % set algorithm parameters PDC.sig=1e-2; % PDC.rho=1.2; % diff --git a/Example/DeconvEx/2D/Deconv_KL_NonNeg_NoReg.m b/Example/DeconvEx/2D/Deconv_KL_NonNeg_NoReg.m index ef3a689..dea8ad2 100644 --- a/Example/DeconvEx/2D/Deconv_KL_NonNeg_NoReg.m +++ b/Example/DeconvEx/2D/Deconv_KL_NonNeg_NoReg.m @@ -47,7 +47,7 @@ % -- FISTA KL + NonNeg FBS=OptiFBS(KL*H,R_POS); -FBS.OutOp=OutputOpti(1,im,40); +FBS.OutOp=OutputOptiSNR(1,im,40); FBS.ItUpOut=5; % call OutputOpti update every ItUpOut iterations FBS.fista=true; % activate fista FBS.maxiter=200; % max number of iterations @@ -56,7 +56,7 @@ % -- Richardson-Lucy KL + NonNeg (implicit) RL=OptiRichLucy(KL*H); -RL.OutOp=OutputOpti(1,im,40); +RL.OutOp=OutputOptiSNR(1,im,40); RL.ItUpOut=5; % call OutputOpti update every ItUpOut iterations RL.maxiter=200; % max number of iterations RL.run(y); % run the algorithm diff --git a/Example/DeconvEx/2D/Deconv_KL_TV_NonNeg.m b/Example/DeconvEx/2D/Deconv_KL_TV_NonNeg.m index 2be91e1..fbd32da 100644 --- a/Example/DeconvEx/2D/Deconv_KL_TV_NonNeg.m +++ b/Example/DeconvEx/2D/Deconv_KL_TV_NonNeg.m @@ -57,7 +57,7 @@ Hn={H,G,LinOpDiag(sz)}; rho_n=[1e-3,1e-3,1e-3]; ADMM=OptiADMM([],Fn,Hn,rho_n); -ADMM.OutOp=OutputOpti(1,im,40,[1 2]); +ADMM.OutOp=OutputOptiSNR(1,im,40,[1 2]); ADMM.ItUpOut=2; % call OutputOpti update every ItUpOut iterations ADMM.maxiter=200; % max number of iterations ADMM.run(y); % run the algorithm @@ -66,7 +66,7 @@ Fn={lamb*R_N12,F}; Hn={G,H}; PDC=OptiPrimalDualCondat([],R_POS,Fn,Hn); -PDC.OutOp=OutputOpti(1,im,40,[2 3]); +PDC.OutOp=OutputOptiSNR(1,im,40,[2 3]); PDC.tau=100; % set algorithm parameters PDC.sig=1e-2; % PDC.rho=1.2; % @@ -76,7 +76,7 @@ % -- Richardson-Lucy-TV KL + TV + NonNeg (implicit) RLTV=OptiRichLucy(F*H,1,lamb); -RLTV.OutOp=OutputOpti(1,im,40); +RLTV.OutOp=OutputOptiSNR(1,im,40); RLTV.ItUpOut=2; % call OutputOpti update every ItUpOut iterations RLTV.maxiter=200; % max number of iterations RLTV.run(y); % run the algorithm @@ -86,7 +86,7 @@ C = F*H+ lamb*hyperB; C.memoizeOpts.apply=true; VMLMB=OptiVMLMB(C,0.,[]); -VMLMB.OutOp=OutputOpti(1,im,40); +VMLMB.OutOp=OutputOptiSNR(1,im,40); VMLMB.ItUpOut=2; VMLMB.maxiter=200; % max number of iterations VMLMB.m=3; % number of memorized step in hessian approximation diff --git a/Example/DeconvEx/2D/Deconv_LS_HessSchatt.m b/Example/DeconvEx/2D/Deconv_LS_HessSchatt.m index 3fdd397..d7cb6d8 100644 --- a/Example/DeconvEx/2D/Deconv_LS_HessSchatt.m +++ b/Example/DeconvEx/2D/Deconv_LS_HessSchatt.m @@ -51,7 +51,7 @@ % -- Chambolle-Pock LS + ShattenHess CP=OptiChambPock(lamb*R_1sch,Hess,F); -CP.OutOp=OutputOpti(1,im,20); +CP.OutOp=OutputOptiSNR(1,im,20); CP.CvOp=TestCvgCombine(TestCvgCostRelative(1e-4), 'StepRelative',1e-4); CP.tau=1; % algorithm parameters CP.sig=0.02; % @@ -65,7 +65,7 @@ rho_n=[1e-1]; ADMM=OptiADMM(F,Fn,Hn,rho_n); ADMM.CvOp=TestCvgCombine(TestCvgCostRelative(1e-4), 'StepRelative',1e-4); -ADMM.OutOp=OutputOpti(1,im,10); +ADMM.OutOp=OutputOptiSNR(1,im,10); ADMM.ItUpOut=1; % call OutputOpti update every ItUpOut iterations ADMM.maxiter=200; % max number of iterations ADMM.run(zeros(size(y))); % run the algorithm diff --git a/Example/DeconvEx/2D/Deconv_LS_HessSchatt_NonNeg.m b/Example/DeconvEx/2D/Deconv_LS_HessSchatt_NonNeg.m index 0bdbb08..7221a61 100644 --- a/Example/DeconvEx/2D/Deconv_LS_HessSchatt_NonNeg.m +++ b/Example/DeconvEx/2D/Deconv_LS_HessSchatt_NonNeg.m @@ -57,7 +57,7 @@ rho_n=[1e-1,1e-1]; ADMM=OptiADMM(F,Fn,Hn,rho_n); ADMM.CvOp=TestCvgCombine(TestCvgCostRelative(1e-4,[1 2]), 'StepRelative',1e-4); -ADMM.OutOp=OutputOpti(1,im,10,[1 2]); +ADMM.OutOp=OutputOptiSNR(1,im,10,[1 2]); ADMM.ItUpOut=1; % call OutputOpti update every ItUpOut iterations ADMM.maxiter=200; % max number of iterations ADMM.run(zeros(size(y))); % run the algorithm @@ -67,7 +67,7 @@ Hn={Hess}; PDC=OptiPrimalDualCondat(F,R_POS,Fn,Hn); PDC.CvOp=TestCvgCombine(TestCvgCostRelative(1e-4,[1 3]), 'StepRelative',1e-4); -PDC.OutOp=OutputOpti(1,im,40,[1 3]); +PDC.OutOp=OutputOptiSNR(1,im,40,[1 3]); PDC.tau=1; % set algorithm parameters PDC.sig=1e-2; % PDC.rho=1.7; % diff --git a/Example/DeconvEx/2D/Deconv_LS_NoReg.m b/Example/DeconvEx/2D/Deconv_LS_NoReg.m index aa9ea4b..2b362e8 100644 --- a/Example/DeconvEx/2D/Deconv_LS_NoReg.m +++ b/Example/DeconvEx/2D/Deconv_LS_NoReg.m @@ -57,14 +57,14 @@ % -- Gradient Descent LS GD=OptiGradDsct(F); -GD.OutOp=OutputOpti(1,im,40); +GD.OutOp=OutputOptiSNR(1,im,40); GD.ItUpOut=2; % call OutputOpti update every ItUpOut iterations GD.maxiter=200; % max number of iterations GD.run(y); % run the algorithm (Note that gam is fixed automatically to 1/F.lip here since F.lip is defined and since we do not have setted gam) %% -- VMLMB LS VMLMB=OptiVMLMB(F,[],[]); -VMLMB.OutOp=OutputOpti(1,im,40); +VMLMB.OutOp=OutputOptiSNR(1,im,40); VMLMB.ItUpOut=2; VMLMB.maxiter=200; % max number of iterations VMLMB.m=2; % number of memorized step in hessian approximation (one step is enough for quadratic function) @@ -112,11 +112,18 @@ valVMLMB=VMLMB.OutOp.evolcost;save('Util/UnitTest/Data/Deconv_LS_NoReg_VMLMB','valVMLMB'); valCG=CG.OutOp.evolcost;save('Util/UnitTest/Data/Deconv_LS_NoReg_CG','valCG'); else + if (exist('Util/UnitTest/Data/Deconv_LS_NoReg_GD.mat','file')==2) && ... + (exist('Util/UnitTest/Data/Deconv_LS_NoReg_VMLMB.mat','file')==2) && ... + (exist('Util/UnitTest/Data/Deconv_LS_NoReg_CG.mat','file')==2) + load('Util/UnitTest/Data/Deconv_LS_NoReg_GD'); load('Util/UnitTest/Data/Deconv_LS_NoReg_VMLMB'); load('Util/UnitTest/Data/Deconv_LS_NoReg_CG'); stateTest=isequal(valGD,GD.OutOp.evolcost) && isequal(valVMLMB,VMLMB.OutOp.evolcost) && isequal(valCG,CG.OutOp.evolcost); message=['Max error between costs evolution : ',num2str(max([norm(valGD-GD.OutOp.evolcost),... norm(valVMLMB-VMLMB.OutOp.evolcost),norm(valCG-CG.OutOp.evolcost)]))]; + else + error('you must generate data before setting : set generateDataUnitTests=1 in UnitScript.m'); + end end end \ No newline at end of file diff --git a/Example/DeconvEx/2D/Deconv_LS_NonNeg_NoReg.m b/Example/DeconvEx/2D/Deconv_LS_NonNeg_NoReg.m index f33324e..e6a11a4 100644 --- a/Example/DeconvEx/2D/Deconv_LS_NonNeg_NoReg.m +++ b/Example/DeconvEx/2D/Deconv_LS_NonNeg_NoReg.m @@ -3,10 +3,10 @@ % the Least-Squares function plus the NonNegativity constraint % without any regularizer: % 0.5 ||Hx - y||^2 + i_{>0}(x) -% using FISTA and VMLMB +% using FISTA, Douglas-Rachford and VMLMB % % See LinOp, LinOpConv, Cost, CostL2, CostNonNeg, Opti, -% OptiFBS, OptiVMLMB, OutpuOpti +% OptiFBS, OptiVMLMB, OptiDouglasRachford, OutpuOpti %------------------------------------------------------------ close all; help Deconv_Ls_NonNeg_NoReg @@ -48,19 +48,25 @@ % -- FISTA LS + NonNeg FBS=OptiFBS(F,R_POS); -FBS.OutOp=OutputOpti(1,im,10); +FBS.OutOp=OutputOptiSNR(1,im,10); FBS.CvOp=TestCvgCombine(TestCvgCostRelative(1e-4), 'StepRelative',1e-4); FBS.ItUpOut=2; % call OutputOpti update every ItUpOut iterations FBS.fista=true; % activate fista FBS.maxiter=200; % max number of iterations FBS.run(zeros(size(y)));% run the algorithm (Note that gam is fixed automatically to 1/F.lip here since F.lip is defined and since we do not have setted gam) +% -- Douglas-Rachford LS + NonNeg +DR=OptiDouglasRachford(F,R_POS,[],10,1.5); +DR.OutOp=OutputOptiSNR(1,im,10); +DR.CvOp=TestCvgCombine(TestCvgCostRelative(1e-4), 'StepRelative',1e-4); +DR.ItUpOut=2; % call OutputOpti update every ItUpOut iterations +DR.maxiter=200; % max number of iterations +DR.run(y); -%% -- VMLMB LS + NonNeg - +% - VMLMB LS + NonNeg H.memoizeOpts.applyHtH=true; VMLMB=OptiVMLMB(F,0.,[]); -VMLMB.OutOp=OutputOpti(1,im,10); +VMLMB.OutOp=OutputOptiSNR(1,im,10); VMLMB.CvOp=TestCvgCombine('CostRelative',1e-4, 'StepRelative',1e-4); % identical to VMLMB.CvOp=TestCvgCombine(TestCvgCostRelative(1e-4),TestCvgStepRelative(1e-5)); VMLMB.ItUpOut=2; VMLMB.maxiter=200; % max number of iterations @@ -70,18 +76,22 @@ %% -- Display imdisp(FBS.xopt,'LS + NonNeg (FISTA)',1); imdisp(VMLMB.xopt,'LS+POS (VMLMB)',1); +imdisp(DR.xopt,'LS+POS (DR)',1); figure;plot(FBS.OutOp.iternum,FBS.OutOp.evolcost,'LineWidth',1.5);grid; set(gca,'FontSize',12); -hold all;loglog(VMLMB.OutOp.iternum,VMLMB.OutOp.evolcost,'LineWidth',1.5);grid; set(gca,'FontSize',12); -xlabel('Iterations');ylabel('Cost');legend('LS+POS (FISTA)','LS+POS (VMLMB)');title('Cost evolution'); +hold all;loglog(VMLMB.OutOp.iternum,VMLMB.OutOp.evolcost,'LineWidth',1.5); +loglog(DR.OutOp.iternum,DR.OutOp.evolcost,'LineWidth',1.5); +xlabel('Iterations');ylabel('Cost');legend('LS+POS (FISTA)','LS+POS (VMLMB)','LS+POS (DR)');title('Cost evolution'); figure;subplot(1,2,1); grid; hold all; title('Evolution SNR');set(gca,'FontSize',12); semilogy(FBS.OutOp.iternum,FBS.OutOp.evolsnr,'LineWidth',1.5); semilogy(VMLMB.OutOp.iternum,VMLMB.OutOp.evolsnr,'LineWidth',1.5); +semilogy(DR.OutOp.iternum,DR.OutOp.evolsnr,'LineWidth',1.5); xlabel('Iterations');ylabel('SNR (dB)');legend('LS+POS (FISTA)','LS+POS (VMLMB)','Location','southeast'); subplot(1,2,2);hold on; grid; title('Runing Time (200 iterations)');set(gca,'FontSize',12); orderCol=get(gca,'ColorOrder'); bar(1,[FBS.time],'FaceColor',orderCol(1,:),'EdgeColor','k'); -bar(2,[VMLMB.time],'FaceColor',orderCol(1,:),'EdgeColor','k'); +bar(2,[VMLMB.time],'FaceColor',orderCol(2,:),'EdgeColor','k'); +bar(3,[DR.time],'FaceColor',orderCol(3,:),'EdgeColor','k'); set(gca,'xtick',[1 2]);ylabel('Time (s)'); set(gca,'xticklabels',{'LS+POS (FISTA)','LS+POS (VMLMB)'});set(gca,'XTickLabelRotation',45) @@ -91,11 +101,13 @@ if generateDataUnitTests valFBS=FBS.OutOp.evolcost;save('Util/UnitTest/Data/Deconv_LS_NonNeg_NoReg_FBS','valFBS'); valVMLMB=VMLMB.OutOp.evolcost;save('Util/UnitTest/Data/Deconv_LS_NonNeg_NoReg_VMLMB','valVMLMB'); + valDR=DR.OutOp.evolcost;save('Util/UnitTest/Data/Deconv_LS_NonNeg_NoReg_DR','valDR'); else load('Util/UnitTest/Data/Deconv_LS_NonNeg_NoReg_FBS'); load('Util/UnitTest/Data/Deconv_LS_NonNeg_NoReg_VMLMB'); - stateTest=isequal(valFBS,FBS.OutOp.evolcost) && isequal(valVMLMB,VMLMB.OutOp.evolcost); - message=['Max error between costs evolution : ',num2str(norm(valFBS-FBS.OutOp.evolcost),... - norm(valVMLMB-VMLMB.OutOp.evolcost))]; + load('Util/UnitTest/Data/Deconv_LS_NonNeg_NoReg_DR'); + stateTest=isequal(valFBS,FBS.OutOp.evolcost) && isequal(valVMLMB,VMLMB.OutOp.evolcost) && isequal(valDR,DR.OutOp.evolcost); + message=['Max error between costs evolution : ',num2str(max([norm(valFBS-FBS.OutOp.evolcost),... + norm(valVMLMB-VMLMB.OutOp.evolcost),norm(valDR-DR.OutOp.evolcost)]))]; end end \ No newline at end of file diff --git a/Example/DeconvEx/2D/Deconv_LS_TV.m b/Example/DeconvEx/2D/Deconv_LS_TV.m index a3cc53a..176278b 100644 --- a/Example/DeconvEx/2D/Deconv_LS_TV.m +++ b/Example/DeconvEx/2D/Deconv_LS_TV.m @@ -51,7 +51,7 @@ % -- Chambolle-Pock LS + TV CP=OptiChambPock(lamb*R_N12,G,F); -CP.OutOp=OutputOpti(1,im,10); +CP.OutOp=OutputOptiSNR(1,im,10); CP.CvOp=TestCvgCombine(TestCvgCostRelative(1e-4), 'StepRelative',1e-4); CP.tau=15; % algorithm parameters CP.sig=1/(CP.tau*G.norm^2)*0.99; % @@ -65,7 +65,7 @@ % Here no solver needed in ADMM since the operator H'*H + alpha*G'*G is invertible ADMM=OptiADMM(F,Fn,Hn,rho_n); ADMM.CvOp=TestCvgCombine(TestCvgCostRelative(1e-4), 'StepRelative',1e-4); -ADMM.OutOp=OutputOpti(1,im,10); +ADMM.OutOp=OutputOptiSNR(1,im,10); ADMM.ItUpOut=2; % call OutputOpti update every ItUpOut iterations ADMM.maxiter=200; % max number of iterations ADMM.run(zeros(size(y))); % run the algorithm diff --git a/Example/DeconvEx/2D/Deconv_LS_TV_NonNeg.m b/Example/DeconvEx/2D/Deconv_LS_TV_NonNeg.m index 2787733..3d026ab 100644 --- a/Example/DeconvEx/2D/Deconv_LS_TV_NonNeg.m +++ b/Example/DeconvEx/2D/Deconv_LS_TV_NonNeg.m @@ -62,7 +62,7 @@ Hn={G,LinOpIdentity(sz)}; rho_n=[1e-1,1e-1]; ADMM=OptiADMM(F,Fn,Hn,rho_n); -ADMM.OutOp=OutputOpti(1,im,10,[1 2]); +ADMM.OutOp=OutputOptiSNR(1,im,10,[1 2]); % STOP when the sum successives C = F*x + Fn{1}*Hn{1}*x is lower than 1e-4 or when the distance between two successive step is lower than 1e-5 ADMM.CvOp=TestCvgCombine(TestCvgCostRelative(1e-4,[1 2]), 'StepRelative',1e-4); ADMM.ItUpOut=2; % call OutputOpti update every ItUpOut iterations @@ -73,7 +73,7 @@ Fn={lamb*R_N12}; Hn={G}; PDC=OptiPrimalDualCondat(F,R_POS,Fn,Hn); -PDC.OutOp=OutputOpti(1,im,40,[1 3]); +PDC.OutOp=OutputOptiSNR(1,im,40,[1 3]); % STOP when the sum successives C = F*x + Fn*Hn*x is lower than 1e-4 or when the distance between two successive step is lower than 1e-5 PDC.CvOp=TestCvgCombine(TestCvgCostRelative(1e-4,[1 3]) , 'StepRelative',1e-4); PDC.tau=1; % set algorithm parameters @@ -89,7 +89,7 @@ C = F+ lamb*hyperB; C.memoizeOpts.apply=true; VMLMB=OptiVMLMB(C,0.,[]); -VMLMB.OutOp=OutputOpti(1,im,10); +VMLMB.OutOp=OutputOptiSNR(1,im,10); VMLMB.CvOp=TestCvgCombine('CostRelative',1e-4, 'StepRelative',1e-4); % identical to VMLMB.CvOp=TestCvgCombine(TestCvgCostRelative(1e-4),TestCvgStepRelative(1e-5)); VMLMB.ItUpOut=2; VMLMB.maxiter=200; % max number of iterations diff --git a/Example/DeconvEx/3D/DeconvScript.m b/Example/DeconvEx/3D/DeconvScript.m index c357920..ade59c1 100644 --- a/Example/DeconvEx/3D/DeconvScript.m +++ b/Example/DeconvEx/3D/DeconvScript.m @@ -93,7 +93,7 @@ rho_n=[1e-3,1e-3,1e-3]; % Multipliers rho_n ADMM=OptiADMM([],Fn,Hn,rho_n); end -ADMM.OutOp=OutputOpti(1,im,round(maxIt/10),[1 2]); +ADMM.OutOp=OutputOptiSNR(1,im,round(maxIt/10),[1 2]); ADMM.CvOp=TestCvgCombine(TestCvgCostRelative(1e-4), 'StepRelative',1e-4); ADMM.ItUpOut=round(maxIt/10); ADMM.maxiter=maxIt; diff --git a/LinOp/Deprecated/LinOpCpx.m b/LinOp/Deprecated/LinOpCpx.m index f7c3b7d..e331600 100644 --- a/LinOp/Deprecated/LinOpCpx.m +++ b/LinOp/Deprecated/LinOpCpx.m @@ -35,6 +35,7 @@ %% Constructor methods function this = LinOpCpx(sz) + warning('LinOpCpx Deprecated'); this.name ='LinOpCpx'; assert(issize(sz),'The input size sz should be a conformable to a size '); this.isInvertible = true; diff --git a/LinOp/LinOpBroadcast.m b/LinOp/LinOpBroadcast.m index 28ac004..d36e20c 100644 --- a/LinOp/LinOpBroadcast.m +++ b/LinOp/LinOpBroadcast.m @@ -16,6 +16,7 @@ %% Copyright (C) 2018 % F. Soulez ferreol.soulez@epfl.ch + % Anthony Berdeu (Laboratoire Hubert Curien) % % This program is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by @@ -96,7 +97,7 @@ for n=this.index x = sum(x,n); end - y = squeeze(x); + y = reshape(x, this.sizein); end function y = applyHtH_(this,x) % Reimplemented from parent class :class:`LinOp`. diff --git a/LinOp/LinOpConv.m b/LinOp/LinOpConv.m index d30bf83..051604f 100644 --- a/LinOp/LinOpConv.m +++ b/LinOp/LinOpConv.m @@ -162,7 +162,7 @@ if centering for n=1:ndims(psf) if any(index==n) - psf = fftshift(psf,n); + psf = ifftshift(psf,n); end end diff --git a/LinOp/LinOpGrad.m b/LinOp/LinOpGrad.m index 1325594..6c6c715 100644 --- a/LinOp/LinOpGrad.m +++ b/LinOp/LinOpGrad.m @@ -10,8 +10,8 @@ % All attributes of parent class :class:`LinOp` are inherited. % % **Note** When circular boundary conditions are selected, the method - % makeHtH (or equivalently the composition H'*H) returns a convolution - % linear operator :class:`LinOp` + % makeHtH (or equivalently the composition ``H'*H``) returns a convolution + % linear operator :class:`LinOp` % % **Example** G = LinOpGrad(sz,index,bc,res) % @@ -210,7 +210,7 @@ % Reimplemented from parent class :class:`LinOp`. if strcmp(this.bc,'circular') fHtH=zeros_(this.sizein) ; - idxAll = repmat({1}, 1, this.ndms); + idxAll = repmat({1}, 1, max(this.ndms,2)); % The max(this.ndms,2) is to deal with the vectorial case (i.e. this.ndms==1) rep=this.sizein; sel=idxAll; for ii = this.index diff --git a/LinOp/LinOpHess.m b/LinOp/LinOpHess.m index aea9c22..96195a0 100644 --- a/LinOp/LinOpHess.m +++ b/LinOp/LinOpHess.m @@ -10,13 +10,15 @@ % the index vector. The last dimension of length lidx(lidx+1)/2 % represents the coefficients of upper-right part of the symmetric hessian % matrix ordered from top to bottom and left to right. For instance: - % - in 3D with index=[1 2 3], [d^2F/dxx;d^2F/dxy;d^2F/dxz;d^2F/dyy;d^2F/dyz;d^2F/dzz] - % - in 3D with index=[2 3], [d^2F/dyy;d^2F/dyz;d^2F/dzz] + % + % - in 3D with index=[1 2 3], [d^2F/dxx;d^2F/dxy;d^2F/dxz;d^2F/dyy;d^2F/dyz;d^2F/dzz] + % + % - in 3D with index=[2 3], [d^2F/dyy;d^2F/dyz;d^2F/dzz] % % All attributes of parent class :class:`LinOp` are inherited. % % **Note** When circular boundary conditions are selected, the method - % makeHtH (or equivalently the composition H'*H) returns a convolution + % makeHtH (or equivalently the composition ``H'*H``) returns a convolution % linear operator :class:`LinOp` % % **Example** H = LinOpHess(sz,bc,index) diff --git a/LinOp/LinOpMatrix.m b/LinOp/LinOpMatrix.m index 2904d1b..85d67c8 100644 --- a/LinOp/LinOpMatrix.m +++ b/LinOp/LinOpMatrix.m @@ -44,7 +44,11 @@ this.M = M; this.sz = size(M); this.ndms = ndims(M); - this.norm = norm(M); + if issparse(M) + this.norm = normest(M); + else + this.norm = norm(M); + end if nargin == 1 index = []; end diff --git a/LinOp/LinOpSum.m b/LinOp/LinOpSum.m index cd1432f..361a9ff 100644 --- a/LinOp/LinOpSum.m +++ b/LinOp/LinOpSum.m @@ -78,6 +78,8 @@ this.imdims = this.sizein; this.imdims(~T)=1; + this.norm=sqrt(this.sizein(this.index)); + end end diff --git a/LinOp/LinOpXRay.m b/LinOp/LinOpXRay.m new file mode 100644 index 0000000..47fd1c3 --- /dev/null +++ b/LinOp/LinOpXRay.m @@ -0,0 +1,117 @@ +classdef LinOpXRay < LinOp + % LinOpXray: 2D, parallel x-ray operator based on MATLAB's ``radon`` function. + % + % + % :param x: structure defining the reconstruction domain (always 2D) + % :param x.size: 1x2 vector giving dimensions of the domain in pixels, must be of the form ``N*ones(1,2)``. + % :param x.step: 1x2 vector giving the spacing of the pixels + % (in an apropirate unit of length, e.g. milimeteres). + % :param thetas: 1x(number of projections) vector giving the projection directions in radians. +% Specifically, ``[cos(theta) sin(theta)]`` points in the direction of the positive axis in the projection coordinate system. + % + % **Note** We take rows as $$x$$ and columns as $$y$$ in the reconstruction domain. + % **Note** The center of the reconstruction (and rotation axis) is assumed to be at ``floor((x.size+1)/2)``. + % + % **Note** The center of each projection is at ``ceil(sizeout/2)``. + % + % **Note** Be warned that the adjoint is only accurate to around 20 dB for + % for images of size 500x500 px. The error is smaller for larger + % images. + + %% Copyright (C) 2019 + % M. McCann michael.thompson.mccann@gmail.com +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% This program is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. If not, see . + + + properties + x + y + thetas + end + + methods + + function this = LinOpXRay(x, thetas) + this.name ='XRay'; + + % input checking for x + if length(unique(x.size))>1 + error('only square x domains are supported'); + end + + if length(x.size) ~= 2 + error('only 2d supported'); + end + + if length(unique(x.step))>1 + error('only isotropic x domains are supported'); + end + + this.x = x; + this.sizein = x.size; + + % setting output size + y.size = 2*ceil(norm(x.size-floor((x.size-1)/2)-1))+3; + this.sizeout= [y.size length(thetas)]; + + this.y = y; + this.thetas = thetas; + + % make copies of private image processing toolbox functions so + % that we can call them directly + if isempty(which('iradonmex')) + filename = which('iradonmex', '-all'); + privateDir = fullfile(fileparts(mfilename('fullpath')), 'private'); + if ~exist(privateDir, 'dir') + mkdir(privateDir) + end + copyfile(filename{1}, privateDir) + end + + % estimate norm, assuming shift invariance + d = zeros(x.size); + d(ceil(end/2), ceil(end/2)) = 1; + + HTHd = fft2(this.applyAdjoint_(this.apply_(d))); + + this.norm = sqrt(max(abs(HTHd(:)))); + + + end + + + end + + methods (Access = protected) + + function y=apply_(this,x) % apply the operator + y=radon(x,this.thetas/pi*180 - 90, this.y.size) * this.x.step(1); + end + + function g = applyAdjoint_(this,c) % apply the adjoint + xVects = (1:this.x.size(1)) - floor((this.x.size(1)+1)/2); + [X0, X1] = ndgrid(xVects); + g = iradonmex(size(X0,1), this.thetas - pi/2, X1./this.x.step(2), -X0./this.x.step(1), c, 1); % 1 - linear interp + g = g * this.x.step(1); + end + + function y = applyHtH_(this,x) % apply the HtH matrix + y = this.applyAdjoint(this.apply(x)); + end + + end + + +end + diff --git a/LinOp/SelectorLinOps/LinOpSelector.m b/LinOp/SelectorLinOps/LinOpSelector.m index 619e597..5089c4c 100644 --- a/LinOp/SelectorLinOps/LinOpSelector.m +++ b/LinOp/SelectorLinOps/LinOpSelector.m @@ -59,6 +59,7 @@ function y = apply_(this,x) % Reimplemented from parent class :class:`LinOpSelector`. y =x(this.sel{:}); + y = reshape(y, this.sizeout); end function y = applyAdjoint_(this,x) % Reimplemented from parent class :class:`LinOpSelector`. diff --git a/LinOp/SelectorLinOps/LinOpSelectorPatch.m b/LinOp/SelectorLinOps/LinOpSelectorPatch.m index 67b1854..f5588da 100644 --- a/LinOp/SelectorLinOps/LinOpSelectorPatch.m +++ b/LinOp/SelectorLinOps/LinOpSelectorPatch.m @@ -8,6 +8,7 @@ % :param sz: size of \\(\\mathrm{x}\\) on which the :class:`LinOpDownsample` applies. % :param idxmin: array containing the first kept index in each direction % :param idxmax: array containing the last kept index in each direction + % :param squeezeflag: true to squeeze singleton dimensions (default: false) % % All attributes of parent class :class:`LinOpSelector` are inherited. % @@ -18,6 +19,8 @@ %% Copyright (C) 2017 % E. Soubies emmanuel.soubies@epfl.ch + % F. Soulez + % A. Berdeu % % This program is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by @@ -39,7 +42,7 @@ %% Constructor methods - function this = LinOpSelectorPatch(sz,idxmin,idxmax) + function this = LinOpSelectorPatch(sz,idxmin,idxmax,squeezeflag) this.name ='LinOp SelectorPatch'; assert(cmpSize(size(sz),size(idxmin)) && cmpSize(size(sz),size(idxmax)),'Parameters sz idxmin and idxmax must have the same size'); assert(~any(idxmin<=0) && ~any(idxmin>sz),'idxmin out of sz range'); @@ -48,7 +51,16 @@ this.idxmin=idxmin; this.idxmax=idxmax; this.sizeout=idxmax-idxmin+1; - if this.sizeout(end)==1, this.sizeout=this.sizeout(1:end-1);end; + % flag_squeeze + if (nargin>3)&& ~isscalar(squeezeflag ) + error('squeezeflag must be a single boolean'); + end + if (nargin>3) && squeezeflag + this.sizeout = this.sizeout(this.sizeout~=1); + end + if numel(this.sizeout)==1 + this.sizeout = [this.sizeout, 1]; + end this.sizein=sz; this.sel=cell(length(sz),1); for ii=1:length(sz) diff --git a/LinOp/Tests/testLinOpXray.m b/LinOp/Tests/testLinOpXray.m new file mode 100644 index 0000000..31ede72 --- /dev/null +++ b/LinOp/Tests/testLinOpXray.m @@ -0,0 +1,35 @@ +x.size = [125 125]; +x.step = [1 1]; + +numThetas = 100; +thetas = linspace(0, pi, numThetas+1); +thetas(end) = []; % we don't want to have both 0 and pi + +H = LinOpXRay(x, thetas); + +% make some arbitrary phantom +[X, Y] = ndgrid(1:x.size(1), 1:x.size(2)); +c = (X-50).^2 + (Y-50).^2 < 40^2; +c = c & (((X-40)/1.5).^2 + (Y-24).^2 > 10^2); +c = c & (((X-40)/1.5).^2 + (Y-60).^2 > 20^2); +c = c & ((X-60).^2 + (Y-35).^2 > 5^2); +c = c & ((X-77).^2 + ((Y-35)./1.5).^2 > 8^2); +c = c .* (sqrt( (X - x.size(1)).^2 + (Y- x.size(2)).^2 ) ); +c = reshape( double(c), x.size); + +%% + +g = H*c; +HTg = H'*g; + +%% +d = zeros(x.size); +d(ceil(end/2), ceil(end/2)) = 1; + +HTHd = H'*H*d; + + +%% +checkLinOp(H) + + diff --git a/NEWS.md b/NEWS.md index fa41a63..71f649c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,7 +1,29 @@ +## v1.1.2 release notes +The v1.1.2 release does not contain any major changes but mainly incremental improvements and bug fixes. + +#### Improvements +- New `LinOpXRay` for 2D radon transform (Note that the adjoint is not really the adjoint, but only an approximation (up to about 20 dB), +- LinOp compatibility with sparse matrices, +- `Util/estimateNorm.m` for estimation of the norm of an operator, +- New `OutputOptiSNR` to print SNR during optimization, +- Better summation and multiplication of Map objects with scalars, +- Rename `OpEWSquareRoot` to` OpEWSqrt`, +- Add some kind of unitary tests, +- `CostHyperBolic`: Epsilon can be a vector, +- non scalar step size in `OptiFBS` + + +#### Bug fixes +Many bug fixes, among them: +- Generalize Compiler option for C/C++ files and reactivate openmp for the Schatten norm +- Better handling of singletons dimensions. + + + ## v1.1.1 release notes The v1.1.1 release does not contain any major changes but mainly incremental improvements and bug fixes. -### Improvements +#### Improvements - Sfft functions: faster implementation, work with any dimension, new padding option argument. - Merge LinOpDFT into LinOpSDFT. LinOpDFT remains now as an alias. - Remove code duplication in SelectorLinOps. diff --git a/NonLinOp/ElementWiseOp/OpEWSquareRoot.m b/NonLinOp/ElementWiseOp/OpEWSqrt.m similarity index 93% rename from NonLinOp/ElementWiseOp/OpEWSquareRoot.m rename to NonLinOp/ElementWiseOp/OpEWSqrt.m index a9d1591..67a1efd 100644 --- a/NonLinOp/ElementWiseOp/OpEWSquareRoot.m +++ b/NonLinOp/ElementWiseOp/OpEWSqrt.m @@ -1,4 +1,4 @@ -classdef OpEWSquareRoot < Map +classdef OpEWSqrt < Map % Element-wise square root operator: % $$[\\mathrm{H}(x)]_k = \\sqrt{(x_k)} $$ % @@ -28,8 +28,8 @@ %% Constructor methods - function this = OpEWSquareRoot(sz) - this.name ='OpEWSquareRoot'; + function this = OpEWSqrt(sz) + this.name ='OpEWSqrt'; this.sizein=sz; this.sizeout=sz; this.isDifferentiable=true; diff --git a/NonLinOp/Tests/TestCompositionEWop.m b/NonLinOp/Tests/TestCompositionEWop.m index d814ba2..ea90eb9 100644 --- a/NonLinOp/Tests/TestCompositionEWop.m +++ b/NonLinOp/Tests/TestCompositionEWop.m @@ -9,7 +9,7 @@ sq=OpEWSquaredMagnitude(G.sizeout); % Elem-wise square magnitude for numerator S=LinOpSum(G.sizeout,3); % To sum the two squared gradient components iv=OpEWInverse(sz); % Elem-wise inverse -sqr=OpEWSquareRoot(sz); % Elem-wise square root +sqr=OpEWSqrt(sz); % Elem-wise square root sq2=OpEWSquaredMagnitude(sz); % Elem-wise square magnitude for denominator (different size) SS=LinOpSum(sz); % Final sum of all pixels @@ -28,7 +28,7 @@ G=LinOpGrad(sz); % Gradient operator sq=OpEWSquaredMagnitude(G.sizeout); % Elem-wise square magnitude for numerator S=LinOpSum(G.sizeout,3); % To sum the two squared gradient components -sqr=OpEWSquareRoot(sz); % Elem-wise square root +sqr=OpEWSqrt(sz); % Elem-wise square root SS=LinOpSum(sz); % Final sum of all pixels op=SS*sqr*S*sq*G; % Combinaison of Maps diff --git a/Opti/OptiADMM.m b/Opti/OptiADMM.m index 00322c5..5fc2717 100644 --- a/Opti/OptiADMM.m +++ b/Opti/OptiADMM.m @@ -90,8 +90,10 @@ this.Hn=Hn; this.rho_n=rho_n; if ~isempty(F0) - assert(~isempty(solver) || isa(F0,'CostL2') || isa(F0,'CostL2Composition') || ... - ( isa(F0,'CostSummation') && all(cellfun(@(x) (isa(x,'CostL2Composition') || isa(x,'CostL2')), F0.mapsCell))), ... + assert(~isempty(solver) || isa(F0,'CostL2') || isa(F0,'CostL2Composition') || ... + (isa(F0,'CostMultiplication') && F0.isnum && (isa(F0.cost2,'CostL2') || isa(F0.cost2,'CostL2Composition'))) || ... + ( isa(F0,'CostSummation') && all(cellfun(@(x) (isa(x,'CostL2Composition') || isa(x,'CostL2') || ... + (isa(x,'CostMultiplication') && x.isnum && (isa(x.cost2,'CostL2') || isa(x.cost2,'CostL2Composition'))) ), F0.mapsCell))), ... 'when F0 is nonempty and is not CostL2 or CostL2Composition (or a sum of them), a solver must be given (see help)'); this.cost=F0 +Fn{1}*Hn{1}; else @@ -113,14 +115,30 @@ elseif isa(this.F0,'CostL2') this.A=this.A + LinOpDiag(this.A.sizein); this.b0=this.F0.y; + elseif isa(this.F0,'CostMultiplication') && F0.isnum + if isa(F0.cost2,'CostL2') + this.A=this.A + this.F0.cost1 * LinOpDiag(this.A.sizein); + this.b0=this.F0.cost1 * this.F0.cost2.y; + elseif isa(F0.cost2,'CostL2Composition') + this.A=this.A + this.F0.cost1 * this.F0.cost2.H2'*(this.F0.cost2.H1.W*this.F0.cost2.H2); + this.b0=this.F0.cost1 * this.F0.cost2.H2'*this.F0.cost2.H1.y; + else + error('If F0 is not a CostL2 / CostL2Composition / CostSummation of them, a solver is required in ADMM'); + end elseif isa(this.F0,'CostSummation') for n=1:length(this.F0.mapsCell) if isa(this.F0.mapsCell{n},'CostL2Composition') this.A=this.A + this.F0.alpha(n)*this.F0.mapsCell{n}.H2'*this.F0.mapsCell{n}.H2; this.b0=this.b0 + this.F0.alpha(n)*this.F0.mapsCell{n}.H2'*this.F0.mapsCell{n}.H1.y; - else + elseif isa(this.F0.mapsCell{n},'CostL2') this.A=this.A + this.F0.alpha(n)*LinOpDiag(this.A.sizein); this.b0=this.b0 + this.F0.alpha(n)*this.F0.mapsCell{n}.y; + elseif isa(this.F0.mapsCell{n}.cost2,'CostL2Composition') % If CostMultiplication with a scalar + this.A=this.A + this.F0.alpha(n)*this.F0.mapsCell{n}.cost1*this.F0.mapsCell{n}.cost2.H2'*this.F0.mapsCell{n}.cost2.H2; + this.b0=this.b0 + this.F0.alpha(n)*this.F0.mapsCell{n}.cost1*this.F0.mapsCell{n}.cost2.H2'*this.F0.mapsCell{n}.cost2.H1.y; + elseif isa(this.F0.mapsCell{n}.cost2,'CostL2') + this.A=this.A + this.F0.alpha(n)*this.F0.mapsCell{n}.cost1*LinOpDiag(this.A.sizein); + this.b0=this.b0 + this.F0.alpha(n)*this.F0.mapsCell{n}.cost1*this.F0.mapsCell{n}.cost2.y; end end else diff --git a/Opti/OptiDouglasRachford.m b/Opti/OptiDouglasRachford.m index 46665b7..9110034 100644 --- a/Opti/OptiDouglasRachford.m +++ b/Opti/OptiDouglasRachford.m @@ -31,37 +31,22 @@ % along with this program. If not, see . properties - nu = 1; - useL = 0; F1; F2; lambda gamma - L; end properties (SetAccess = protected,GetAccess = protected) - y; - Ly; + y; end methods function this = OptiDouglasRachford(F1, F2, L, gamma, lambda) + if isempty(L), L=LinOpIdentity(F2.sizein); end; % F1, F2 this.F1 = F1; - this.F2 = F2; - - % L - if exist('L', 'var') && ~isempty(L) - this.useL = 1; - r = randn(L.sizeout); - nu = r ./ L.HHt(r); - assert(std(nu(:)) <1e-6, 'LLt != nu I'); - this.nu = real(mean(nu(:))); - if this.nu==1 - this.useL = 2; - end - this.L=L; - end + this.F2 = F2*L; + this.cost=this.F1+this.F2; % gamma if numel(gamma)==1 @@ -84,18 +69,9 @@ function initialize(this,x0) function flag=doIteration(this) % Reimplementation from :class:`Opti`. - if this.useL - this.Ly = this.L*this.y; - if this.useL==2 - this.xopt = this.L.adjoint( this.F2.prox(this.Ly, this.gamma(2))); - else - this.xopt = this.y + (1./this.nu).* this.L.adjoint( this.F2.prox(this.Ly, this.nu.*this.gamma(2)) - this.Ly); - end - else - this.xopt = this.F2.prox(this.y, this.gamma(2)); - end - this.y = this.y + this.lambda .* ( this.F1.prox(2.*this.xopt- this.y,this.gamma(1)) - this.xopt); - flag=OPTI_NEXT_IT; + this.xopt = this.F2.applyProx(this.y, this.gamma(2)); + this.y = this.y + this.lambda .* ( this.F1.applyProx(2.*this.xopt- this.y,this.gamma(1)) - this.xopt); + flag=this.OPTI_NEXT_IT; end end end diff --git a/Opti/OptiFBS.m b/Opti/OptiFBS.m index d563d68..dbf968d 100644 --- a/Opti/OptiFBS.m +++ b/Opti/OptiFBS.m @@ -6,6 +6,7 @@ % :param G: a :class:`Cost` with an implementation of the :meth:`applyProx`. % :param gam: descent step % :param fista: boolean true if the accelerated version FISTA [3] is used (default false) + % :param momRestart: boolean true if the moment restart strategy is used [4] is used (default false) % % All attributes of parent class :class:`Opti` are inherited. % @@ -26,6 +27,8 @@ % [3] Amir Beck and Marc Teboulle, "A Fast Iterative Shrinkage-Thresholding Algorithm for Linear inverse Problems", % SIAM Journal on Imaging Science, vol 2, no. 1, pp 182-202 (2009) % + % [4] Brendan O'donoghue and Emmanuel Candès. 2015. Adaptive Restart for Accelerated Gradient Schemes. Found. Comput. Math. 15, 3 (June 2015), 715-732. + % % **Example** FBS=OptiFBS(F,G) % % See also :class:`Opti` :class:`OutputOpti` :class:`Cost` @@ -57,7 +60,7 @@ properties F; % Cost F G; % Cost G - + momRestart = false; % Donoho put grad gam=[]; % descent step fista=false; % FISTA option [3] @@ -93,15 +96,23 @@ function initialize(this,x0) % Reimplementation from :class:`Opti`. For details see [1-3]. if this.fista % if fista - this.xopt=this.G.applyProx(this.y - this.gam*this.F.applyGrad(this.y),this.gam); - told=this.tk; - this.tk=0.5*(1+sqrt(1+4*this.tk^2)); - this.y=this.xopt + this.alpha*(told-1)/this.tk*(this.xopt - this.xold); + this.xopt=this.G.applyProx(this.y - this.gam.*this.F.applyGrad(this.y),this.gam); + if this.momRestart && dot(this.y(:) - this.xopt(:),this.xopt(:) - this.xold(:)) > 0 + fprint('Restarting\n'); + this.y = this.xold; + this.tk = 1; + this.xopt=this.G.applyProx(this.y - this.gam.*this.F.applyGrad(this.y),this.gam); + else + told=this.tk; + this.tk=0.5*(1+sqrt(1+4*this.tk^2)); + this.y=this.xopt + this.alpha*(told-1)/this.tk*(this.xopt - this.xold); + end else - this.xopt=this.G.applyProx(this.xopt - this.gam*this.F.applyGrad(this.xopt),this.gam); + this.xopt=this.G.applyProx(this.xopt - this.gam.*this.F.applyGrad(this.xopt),this.gam); end - + flag=this.OPTI_NEXT_IT; - end + + end end end diff --git a/Opti/OptiFGP.m b/Opti/OptiFGP.m index a4f6c13..0033071 100644 --- a/Opti/OptiFGP.m +++ b/Opti/OptiFGP.m @@ -72,7 +72,7 @@ this.TV = TV; this.cost = this.F0 + this.TV; - this.gam = 1/8; + this.gam = 1/this.D.norm^2; end function setLambda(this,new_l) diff --git a/Opti/OptiPrimalDualCondat.m b/Opti/OptiPrimalDualCondat.m index 6c4d399..cd71008 100644 --- a/Opti/OptiPrimalDualCondat.m +++ b/Opti/OptiPrimalDualCondat.m @@ -63,9 +63,9 @@ end % Full public properties properties - tau; % parameter of the algorithm - sig; % parameter of the algorithm - rho; % parameter of the algorithm + tau; % parameter of the algorithm + sig; % parameter of the algorithm + rho=1.95; % parameter of the algorithm end methods @@ -125,11 +125,10 @@ function initialize(this,x0) % Update xopt this.xopt=this.rho*xtilde+(1-this.rho)*this.xopt; % Update ytilde and y - ytilde = cell(length(this.Hn),1); for n=1:length(this.Fn) - ytilde{n}=this.Fn{n}.applyProxFench(this.y{n}+this.sig*this.Hn{n}.apply(2*xtilde-this.xold),this.sig); - this.y{n}=this.rho*ytilde{n}+(1-this.rho)*this.y{n}; - end + ytilde=this.Fn{n}.applyProxFench(this.y{n}+this.sig*this.Hn{n}.apply(2*xtilde-this.xold),this.sig); + this.y{n}=this.rho*ytilde +(1-this.rho)*this.y{n}; + end flag=this.OPTI_NEXT_IT; end end diff --git a/Opti/OptiUtils/MatlabOptimPack/installOptimPack.m b/Opti/OptiUtils/MatlabOptimPack/installOptimPack.m index 764d43a..e6d12b5 100644 --- a/Opti/OptiUtils/MatlabOptimPack/installOptimPack.m +++ b/Opti/OptiUtils/MatlabOptimPack/installOptimPack.m @@ -1,6 +1,9 @@ -function installOptimPack() +function installOptimPack(options) %% installOptimPack function % Install OptimPackLegacy and its matlab API from https://github.com/emmt/OptimPackLegacy +% +% You can give as a parameter of this function the path to your GCC +% compiler. Ex: installOptimPack('/usr/bin/gcc-6') % Copyright (C) 2018 F. Soulez ferreol.soulez@epfl.ch % @@ -16,6 +19,10 @@ function installOptimPack() % % You should have received a copy of the GNU General Public License % along with this program. If not, see . +if nargin==0 + options=[]; +end + [mpath,~,~] = fileparts(which('installOptimPack')); disp('Installing OptimPackLegacy'); @@ -25,12 +32,12 @@ function installOptimPack() unzip('OptimPackLegacy'); matDir = 'OptimPackLegacy-master/matlab/'; srcDir = 'OptimPackLegacy-master/src/'; -MexOpt= ['-largeArrayDims ' '-DUSE_BLAS_LIB ' '-DNEW_MATLAB_BLAS ' '-DINT_64BITS ' '-largeArrayDims ' 'COMPFLAGS=''$COMPFLAGS -Wall -mtune=native -fomit-frame-pointer -O2 ''']; +MexOpt= [ '-largeArrayDims ' '-DUSE_BLAS_LIB ' '-DNEW_MATLAB_BLAS ' '-DINT_64BITS ' '-largeArrayDims ' ,options, ' CFLAGS=''-fPIC -Wall -mtune=native -fomit-frame-pointer -O2 ''' ]; CFiles = [srcDir,'opl_vmlmb.c ',srcDir,'opl_algebra.c ',srcDir,'opl_lnsrch.c ',srcDir,'opl_utils.c ','-IOptimPackLegacy-master/src/ ']; eval(['mex ',matDir,'m_opl_vmlmb_get_reason.c ', CFiles,MexOpt]); eval(['mex ',matDir,'m_opl_vmlmb_create.c ', CFiles,MexOpt]); eval(['mex ',matDir,'m_opl_vmlmb_iterate.c ', CFiles,MexOpt]); -eval(['mex ',matDir,'m_opl_vmlmb_restore.c ', CFiles,MexOpt]); +eval(['mex ',matDir,'m_opl_vmlmb_restore.c ', CFiles,MexOpt]); movefile('OptimPackLegacy-master/matlab/*.m'); delete('makeOptimPack.m'); delete('OptimPackLegacy.zip'); diff --git a/Opti/OptiVMLMB.m b/Opti/OptiVMLMB.m index c70c313..56d525b 100644 --- a/Opti/OptiVMLMB.m +++ b/Opti/OptiVMLMB.m @@ -62,19 +62,20 @@ sxtol=0.1; epsilon=0.01; % a small value, in the range [0,1), equals to the cosine of the maximum angle between the search direction and the anti-gradient. The BFGS recursion is restarted, whenever the search direction is not sufficiently "descending". delta=0.1; % DELTA is a small nonegative value used to compute a small initial step. + xmin=[]; + xmax=[]; + active; end properties (SetAccess = protected,GetAccess = public) nparam; task; - xmin=[]; - xmax=[]; neval; bounds =0; ws; end properties (SetAccess = protected,GetAccess = protected) nbeval; - active; + %active; grad; cc; end @@ -127,8 +128,8 @@ function initialize(this,x0) end - this.cc = this.cost.apply(this.xopt); - this.grad = this.cost.applyGrad(this.xopt); + this.cc = gather(this.cost.apply(this.xopt)); + this.grad = gather(real(this.cost.applyGrad(this.xopt))); this.nbeval=this.nbeval+1; @@ -139,7 +140,7 @@ function initialize(this,x0) % Computes next step: this.xopt(1) = this.xopt(1); %Warning : side effect on x0 if x=x0 (due to the fact that x is passed-by-reference in the mexfiles) - this.task = m_opl_vmlmb_iterate(this.ws,this.xopt,this.cc,this.grad,this.active); + this.task = m_opl_vmlmb_iterate(this.ws,gather(this.xopt),this.cc,this.grad,this.active); flag=this.OPTI_REDO_IT; if (this.task == this.OPL_TASK_FG) @@ -155,8 +156,8 @@ function initialize(this,x0) end - this.cc = this.cost.apply(this.xopt); - this.grad = this.cost.applyGrad(this.xopt); + this.cc = gather(this.cost.apply(this.xopt)); + this.grad = gather(real(this.cost.applyGrad(this.xopt))); this.nbeval=this.nbeval+1; @@ -185,12 +186,13 @@ function initialize(this,x0) case 3 this.active = int32( ( (this.xopt>this.xmin) + (this.grad<0) ).*( (this.xopt0) ) ); end + this.active = gather(this.active); else % Convergence, or error, or warning this.endingMessage = ['Convergence, or error, or warning : ',this.task,m_opl_vmlmb_get_reason(this.ws)]; flag=this.OPTI_STOP; - this.task = m_opl_vmlmb_restore(this.ws,this.xopt,this.cc,this.grad); + this.task = m_opl_vmlmb_restore(this.ws,gather(this.xopt),this.cc,this.grad); return; end diff --git a/Opti/OutputOpti/OutputOpti.m b/Opti/OutputOpti/OutputOpti.m index 09ba718..318ae97 100644 --- a/Opti/OutputOpti/OutputOpti.m +++ b/Opti/OutputOpti/OutputOpti.m @@ -17,15 +17,13 @@ % % :param name: name of the :class:`OutputOpti` % :param computecost: boolean, if true the cost function will be computed - % :param xtrue: ground truth to compute the error with the solution (if provided) % :param evolcost: array to save the evolution of the cost function - % :param evolsnr: array to save the evolution of the SNR % :param saveXopt: boolean (defaul true) to save the evolution of the optimized variable xopt. % :param evolxopt: cell saving the optimization variable xopt % :param iterVerb: message will be displayed every iterVerb iterations (must be a multiple of the :attr:`ItUpOut` parameter of classes :class:`Opti`) % :param costIndex: select a specific cost function among a sum in the case where the optimized cost function is a sum of cost functions % - % **Example** OutOpti=OutputOpti(computecost,xtrue,iterVerb) + % **Example** OutOpti=OutputOpti(computecost,iterVerb,costIndex) % % **Important** The update method should have an unique imput that is the :class:`Opti` object in order to % be generic for all Optimization routines. Hence the update method has access (in reading mode) @@ -52,79 +50,93 @@ properties (SetAccess = protected,GetAccess = public) name = 'OutputOpti'% name of the optimization algorithm evolcost; % array saving the evolution of the cost function - evolsnr; % array saving the evolution of the error with the groud truth evolxopt; % cell saving the optimization variable xopt iternum; % array saving the iteration number corresponding to evolcost, evolxopt and evolerr entries + count; % internal counter + + %% To be removed (deprecated) normXtrue; % Norm of the true signal (to compute snr) isgt=false; % Boolean true if Ground Truth is provided - count; % internal counter - xtrue; % Ground Thruth + evolsnr; % array saving the evolution of the error with the groud truth- + xtrue; % Ground Truth + snrOp=[]; % Map to which the coefficients must be applied to compute reconstructed signal + OptiSNR_ = false; + end properties computecost=false; % Boolean, if true the cost function will be computed iterVerb=0; % message will be displayed every iterVerb iterations costIndex=0; % index of the cost function - saveXopt=false; % save evolution of the optimized variable + saveXopt=false; % save evolution of the optimized variable end methods %% Constructor - function this=OutputOpti(computecost,xtrue,iterVerb,costIndex) + function this=OutputOpti(varargin)%(computecost,xtrue,iterVerb,costIndex,snrOp) if nargin>=1 + computecost = varargin{1}; if isscalar(computecost) computecost = (computecost ~= 0); - end - + end assert(islogical(computecost),'Parameter computecost must be logical'); this.computecost=computecost; end - if nargin>=2, this.xtrue=xtrue;end - if nargin>=3 - assert(isscalar(iterVerb) && iterVerb>=0,'Parameter iterVerb must be a positive integer'); + + + if nargin>=2 + if ~isscalar( varargin{2}) + warning('xtrue parameter OutputOpti is deprecated use OutputOptiSNR'); + this=this.OutputOptiSNR_(varargin{:}); + return + end + iterVerb = varargin{2}; + assert(isscalar(iterVerb) && iterVerb>=0 && mod(iterVerb,1)==0,'Parameter iterVerb must be a positive integer'); this.iterVerb=iterVerb; end - if nargin==4, this.costIndex=costIndex;end - if ~isempty(this.xtrue) - this.isgt=true; - this.xtrue=xtrue; - this.normXtrue=norm(this.xtrue(:)); + + if nargin>=3 + costIndex = varargin{3}; + assert(isnumeric(iterVerb) ,'CostIndex must be numeric'); + this.costIndex=costIndex; end + end %% Initialization function init(this) % Initialize the arrays and counters. this.count=1; this.evolcost=zeros_(1); - this.evolsnr=zeros_(1); this.iternum = []; - this.evolxopt = {}; + this.evolxopt = {}; end %% Update method function update(this,opti) % Computes SNR, cost and display evolution. - str=sprintf('Iter: %5i',opti.niter); - if this.computecost + str=sprintf('Iter: %5i',opti.niter); + if this.computecost cc=this.computeCost(opti); - str=sprintf('%s | Cost: %4.4e',str,cc); - this.evolcost(this.count)=cc; - end - if this.isgt + str=sprintf('%s | Cost: %4.4e',str,cc); + this.evolcost(this.count)=cc; + end + + if this.OptiSNR_ && this.isgt snr=this.computeSNR(opti); - str=sprintf('%s | SNR: %4.4e dB',str,snr); - this.evolsnr(this.count)=snr; + str=sprintf('%s | SNR: %4.4e dB',str,snr); + this.evolsnr(this.count)=snr; end + if this.saveXopt this.evolxopt{this.count}=opti.xopt; end - this.iternum(this.count)=opti.niter; - this.count=this.count+1; - if opti.verbose && (opti.niter~=0 && (mod(opti.niter,this.iterVerb)==0) || (opti.niter==1 && this.iterVerb~=0)), - disp(str); - end + this.iternum(this.count)=opti.niter; + this.count=this.count+1; + if opti.verbose && (opti.niter~=0 && (mod(opti.niter,this.iterVerb)==0) || (opti.niter==1 && this.iterVerb~=0)) + disp(str); + end end function cc=computeCost(this,opti) % Evaluate the cost function at the current iterate xopt of - % the given :class:`Opti` opti object + % the given :class:`Opti` opti object if (any(this.costIndex>0) && isa(opti.cost,'CostSummation')) cc = 0; for n=1:numel(this.costIndex) @@ -134,10 +146,45 @@ function update(this,opti) cc = opti.cost*opti.xopt; end end + %% Deprecated functions + function this= OutputOptiSNR_(this,computecost,xtrue,iterVerb,costIndex,snrOp) + this.OptiSNR_=true; + this.name = 'OutputOptiSNR';% name + + if isscalar(computecost) + computecost = (computecost ~= 0); + end + + assert(islogical(computecost),'Parameter computecost must be logical'); + this.computecost=computecost; + + this.xtrue=xtrue; + + if nargin>=4 + assert(isscalar(iterVerb) && iterVerb>=0,'Parameter iterVerb must be a positive integer'); + this.iterVerb=iterVerb; + end + if nargin==5, this.costIndex=costIndex;end + if ~isempty(this.xtrue) + this.isgt=true; + this.xtrue=xtrue; + this.normXtrue=norm(this.xtrue(:)); + end + if nargin==6, this.snrOp = snrOp;end + if ~isempty(this.snrOp) + assert(~isempty(this.xtrue), 'Ground truth must be provided in order to compute SNR'); + assert(isa(this.snrOp, 'Map') && isequal(this.snrOp.sizeout, size(this.xtrue)), ... + 'The SNR operator must have the same output size as the ground truth'); + elseif ~isempty(this.xtrue) + this.snrOp = LinOpDiag(size(this.xtrue)); % Identity operator by default + end + end + function snr=computeSNR(this,opti) % Evaluate the snr for the current iterate xopt of - % the given :class:`Opti` opti object - snr=20*log10(this.normXtrue/norm(this.xtrue(:)-opti.xopt(:))); + % the given :class:`Opti` opti object + reconstruction = this.snrOp.apply(opti.xopt); + snr=20*log10(this.normXtrue/norm(this.xtrue(:)-reconstruction(:))); end end end diff --git a/Opti/OutputOpti/OutputOptiConjGrad.m b/Opti/OutputOpti/OutputOptiConjGrad.m index be061bf..1056bc2 100644 --- a/Opti/OutputOpti/OutputOptiConjGrad.m +++ b/Opti/OutputOpti/OutputOptiConjGrad.m @@ -13,7 +13,9 @@ % :param xtrue: ground truth to compute the error with the solution (if provided) % :param iterVerb: message will be displayed every iterVerb iterations (must be a multiple of the :attr:`ItUpOut` parameter of classes :class:`Opti`) % :param ytWy: weighted norm of $y$ : $ \\mathrm{ytWy} = \\mathrm{ y^T\\,W\\,y}$ + % % See also :class:`OptiConjGrad` :class:`OutputOpti` + % %% Copyright (C) 2018 % F. Soulez ferreol.soulez@epfl.ch @@ -37,25 +39,10 @@ methods function this=OutputOptiConjGrad(computecost,yty,xtrue,iterVerb) + this@OutputOpti(computecost,xtrue,iterVerb); if nargin>=1 - if isscalar(computecost) - computecost = (computecost ~= 0); - end - - assert(islogical(computecost),'Parameter computecost must be logical'); - this.computecost=computecost; this.ytWy = 0.5.*yty; - end - if nargin>=3, this.xtrue=xtrue;end - if nargin>=4 - assert(isscalar(iterVerb) && iterVerb>=0,'Parameter iterVerb must be a positive integer'); - this.iterVerb=iterVerb; - end - if ~isempty(this.xtrue) - this.isgt=true; - this.xtrue=xtrue; - this.normXtrue=norm(this.xtrue(:)); - end + end end function cc=computeCost(this,opti) % Evaluate the cost function at the current iterate xopt of diff --git a/Opti/OutputOpti/OutputOptiSNR.m b/Opti/OutputOpti/OutputOptiSNR.m new file mode 100644 index 0000000..9bb428a --- /dev/null +++ b/Opti/OutputOpti/OutputOptiSNR.m @@ -0,0 +1,145 @@ +classdef OutputOptiSNR < OutputOpti + % OutputOptiSNR class for algorithms displayings and savings + % + % At each :attr:`ItUpOut` iterations of an optimization algorithm (see :class:`Opti` generic class), + % the update method of an :class:`OutputOptiSNR` object will be executed in order to acheive user + % defined computations, e.g., + % + % - compute cost / SNR + % - store current iterate / cost value + % - plot/display stuffs + % + % The present generic class implements a basic update method that: + % + % - display the iteration number + % - computes & display thhandlee cost (if activated) + % - computes & display the SNR if ground truth is provided + % + % :param name: name of the :class:`OutputOptiSNR` + % :param computecost: boolean, if true the cost function will be computed + % :param xtrue: ground truth to compute the error with the solution (if provided) + % :param evolcost: array to save the evolution of the cost function + % :param evolsnr: array to save the evolution of the SNR + % :param saveXopt: boolean (defaul true) to save the evolution of the optimized variable xopt. + % :param evolxopt: cell saving the optimization variable xopt + % :param iterVerb: message will be displayed every iterVerb iterations (must be a multiple of the :attr:`ItUpOut` parameter of classes :class:`Opti`) + % :param costIndex: select a specific cost function among a sum in the case where the optimized cost function is a sum of cost functions + % + % **Example** OutOpti=OutputOptiSNR(computecost,xtrue,iterVerb) + % + % **Important** The update method should have an unique imput that is the :class:`Opti` object in order to + % be generic for all Optimization routines. Hence the update method has access (in reading mode) + % to all the properties of :class:`Opti` objects. + % + % See also :class:`Opti` + + %% Copyright (C) 2017 + % E. Soubies emmanuel.soubies@epfl.ch + % + % This program is free software: you can redistribute it and/or modify + % it under the terms of the GNU General Public License as published by + % the Free Software Foundation, either version 3 of the License, or + % (at your option) any later version. + % + % This program is distributed in the hope that it will be useful, + % but WITHOUT ANY WARRANTY; without even the implied warranty of + % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + % GNU General Public License for more details. + % + % You should have received a copy of the GNU General Public License + % along with this program. If not, see . + + properties (SetAccess = protected,GetAccess = public) + % evolsnr; % array saving the evolution of the error with the groud truth + % xtrue; % Ground Truth + % snrOp=[]; % Map to which the coefficients must be applied to compute reconstructed signal + % normXtrue; % Norm of the true signal (to compute snr) + end + + methods + %% Constructor + function this=OutputOptiSNR (computecost,xtrue,iterVerb,costIndex,snrOp) + + this.name = 'OutputOptiSNR';% name + + if isscalar(computecost) + computecost = (computecost ~= 0); + end + + assert(islogical(computecost),'Parameter computecost must be logical'); + this.computecost=computecost; + + this.xtrue=xtrue; + + if nargin>=3 + assert(isscalar(iterVerb) && iterVerb>=0,'Parameter iterVerb must be a positive integer'); + this.iterVerb=iterVerb; + end + if nargin==4, this.costIndex=costIndex;end + if ~isempty(this.xtrue) + this.xtrue=xtrue; + this.normXtrue=norm(this.xtrue(:)); + end + if nargin==5, this.snrOp = snrOp;end + if ~isempty(this.snrOp) + assert(~isempty(this.xtrue), 'Ground truth must be provided in order to compute SNR'); + assert(isa(this.snrOp, 'Map') && isequal(this.snrOp.sizeout, size(this.xtrue)), ... + 'The SNR operator must have the same output size as the ground truth'); + elseif ~isempty(this.xtrue) + this.snrOp = LinOpDiag(size(this.xtrue)); % Identity operator by default + end + + + end + %% Initialization + function init(this) + % Initialize the arrays and counters. + this.count=1; + this.evolcost=zeros_(1); + this.evolsnr=zeros_(1); + this.iternum = []; + this.evolxopt = {}; + end + %% Update method + function update(this,opti) + % Computes SNR, cost and display evolution. + str=sprintf('Iter: %5i',opti.niter); + if this.computecost + cc=this.computeCost(opti); + str=sprintf('%s | Cost: %4.4e',str,cc); + this.evolcost(this.count)=cc; + end + + snr=this.computeSNR(opti); + str=sprintf('%s | SNR: %4.4e dB',str,snr); + this.evolsnr(this.count)=snr; + + if this.saveXopt + this.evolxopt{this.count}=opti.xopt; + end + this.iternum(this.count)=opti.niter; + this.count=this.count+1; + if opti.verbose && (opti.niter~=0 && (mod(opti.niter,this.iterVerb)==0) || (opti.niter==1 && this.iterVerb~=0)) + disp(str); + end + end + function cc=computeCost(this,opti) + % Evaluate the cost function at the current iterate xopt of + % the given :class:`Opti` opti object + if (any(this.costIndex>0) && isa(opti.cost,'CostSummation')) + cc = 0; + for n=1:numel(this.costIndex) + cc = cc+opti.cost.mapsCell{this.costIndex(n)}*opti.xopt; + end + else + cc = opti.cost*opti.xopt; + end + end + function snr=computeSNR(this,opti) + % Evaluate the snr for the current iterate xopt of + % the given :class:`Opti` opti object + reconstruction = this.snrOp.apply(opti.xopt); + snr=20*log10(this.normXtrue/norm(this.xtrue(:)-reconstruction(:))); + end + end +end diff --git a/Opti/TestCvg/TestCvgADMM.m b/Opti/TestCvg/TestCvgADMM.m index 91929c8..f222274 100644 --- a/Opti/TestCvg/TestCvgADMM.m +++ b/Opti/TestCvg/TestCvgADMM.m @@ -37,6 +37,11 @@ eps_abs; % Termination criterion tolerances. eps_rel; end + properties + evolResPrim + evolResDual + count=0; + end methods %% Constructor @@ -52,7 +57,7 @@ %% Update method function stop = testConvergence(this,opti) % Reimplemented from parent class :class:`TestCvg`. - + this.count=this.count+1; stop = false; p = 0; % Number of constraints @@ -63,15 +68,19 @@ adjHnwn_norm = 0; % ||H'w|| for n = 1:length(opti.wn) p = p + length(opti.yn{n}); - r_norm = r_norm + norm(opti.Hnx{n}-opti.yn{n})^2; - s_norm = s_norm + norm(opti.rho_n(n)*opti.Hn{n}.applyAdjoint(opti.yn{n} - opti.yold{n}))^2; - Hnx_norm = Hnx_norm + norm(opti.Hnx{n})^2; - y_norm = y_norm + norm(opti.yn{n})^2; - adjHnwn_norm = adjHnwn_norm + norm(opti.Hn{n}.applyAdjoint(opti.wn{n}))^2; + r_norm = r_norm + norm(opti.Hnx{n}(:)-opti.yn{n}(:))^2; + tmp=opti.rho_n(n)*opti.Hn{n}.applyAdjoint(opti.yn{n} - opti.yold{n}); + s_norm = s_norm + norm(tmp(:))^2; + Hnx_norm = Hnx_norm + norm(opti.Hnx{n}(:))^2; + y_norm = y_norm + norm(opti.yn{n}(:))^2; + tmp=opti.Hn{n}.applyAdjoint(opti.wn{n}); + adjHnwn_norm = adjHnwn_norm + norm(tmp(:))^2; end eps_primal = sqrt(p)*this.eps_abs + this.eps_rel*sqrt(max(Hnx_norm, y_norm)); eps_dual = sqrt(length(opti.xopt))*this.eps_abs + this.eps_rel*sqrt(adjHnwn_norm); + this.evolResPrim(this.count)=r_norm; + this.evolResDual(this.count)=s_norm; if (sqrt(r_norm) <= eps_primal) && (sqrt(s_norm) <= eps_dual) stop = true; endingMessage = [this.name,' convergence reached']; diff --git a/Opti/TestCvg/TestCvgCostRelative.m b/Opti/TestCvg/TestCvgCostRelative.m index ef986a9..fbc6254 100644 --- a/Opti/TestCvg/TestCvgCostRelative.m +++ b/Opti/TestCvg/TestCvgCostRelative.m @@ -47,7 +47,7 @@ % Tests algorithm convergence from the relative difference between two successive value of the cost function % % :return: boolean true if - % $$ \\frac{\\left| C(\\mathrm{x}^{k}) - C(\\mathrm{x}^{k-1})\\right|}{\\left|C(\\mathrm{x}^{k-1})\\right|} < \\mathrm{costRelativeTol}$$ + % $$ \\frac{\\left| C(\\mathrm{x}^{k}) - C(\\mathrm{x}^{k-1})\\right|}{\\left|C(\\mathrm{x}^{k-1})\\right|} < \\mathrm{costRelativeTol}$$ stop = false; diff --git a/Opti/TestCvg/TestCvgMaxSnr.m b/Opti/TestCvg/TestCvgMaxSnr.m index 63d6e12..1e910e7 100644 --- a/Opti/TestCvg/TestCvgMaxSnr.m +++ b/Opti/TestCvg/TestCvgMaxSnr.m @@ -43,7 +43,7 @@ % Tests algorithm convergence using the SNR with % % :return: boolean true if - % $$ 20\\log\\left(\\frac{\\| \\mathrm{ref}\\|}{\\| \\mathrm{x}^{k} - \\mathrm{ref}\\|}\\right)< 20\\log\\left(\\frac{\\| \\mathrm{ref}\\|}{\\| \\mathrm{x}^{k-1} - \\mathrm{ref}\\|}\\right)$$ + % $$ 20\\log\\left(\\frac{\\| \\mathrm{ref}\\|}{\\| \\mathrm{x}^{k} - \\mathrm{ref}\\|}\\right)< 20\\log\\left(\\frac{\\| \\mathrm{ref}\\|}{\\| \\mathrm{x}^{k-1} - \\mathrm{ref}\\|}\\right)$$ stop = false; assert(checkSize(this.sizeRef, size(opti.xopt)), ' Reference signal and x should be conformable'); diff --git a/Opti/TestCvg/TestCvgStepRelative.m b/Opti/TestCvg/TestCvgStepRelative.m index 4c81450..e45c363 100644 --- a/Opti/TestCvg/TestCvgStepRelative.m +++ b/Opti/TestCvg/TestCvgStepRelative.m @@ -39,7 +39,7 @@ % Tests algorithm convergence from the relative difference between two successive value of the step function % % :return: boolean true if - % $$ \\frac{\\| \\mathrm{x}^{k} - \\mathrm{x}^{k-1}\\|}{\\|\\mathrm{x}^{k-1}\\|} < \\text{stepRelativeTol}.$$ + % $$ \\frac{\\| \\mathrm{x}^{k} - \\mathrm{x}^{k-1}\\|}{\\|\\mathrm{x}^{k-1}\\|} < \\text{stepRelativeTol}.$$ stop = false; if ~isempty(opti.xold) diff --git a/README.md b/README.md index 41672cf..f89ea1a 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ This Matlab toolbox provides a set of tools (operators, cost functions, optimiza ## Getting started -Download or clone this repository and run the script *setMatlabPath.m* which will update your Matlab path with the library. +Download or clone this repository and run the script *setGlobalBioImPath.m* which will update your Matlab path with the library. ## Documentation @@ -15,12 +15,11 @@ A detailled documentation can be found