Skip to content

Commit

Permalink
## v1.1.2 release notes
Browse files Browse the repository at this point in the history 
The v1.1.2 release does not contain any major changes but mainly incremental improvements and bug fixes.

#### Improvements
- New  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,
-  for estimation of the norm of an operator,
- New  to print SNR during optimization,
- Better summation and multiplication of Map objects with scalars,
- Rename  to,
- Add some kind of unitary tests,
- : Epsilon can be a vector,
- non scalar step size in

#### 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.
  • Loading branch information
@FerreolS
FerreolS committed Apr 3, 2019
1 parent 66e07cb commit 95c7214
Show file tree
Hide file tree
Showing 73 changed files with 1,086 additions and 308 deletions.
1 change: 1 addition & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,3 +4,4 @@
*doctrees/
*Opti/OptiUtils/MatlabOptimPack/
.DS_Store
*Util/UnitTest/Data/
3 changes: 1 addition & 2 deletions Abstract/Cost.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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));
Expand Down Expand Up @@ -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

Expand Down
175 changes: 135 additions & 40 deletions Abstract/Map.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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}\\).
Expand Down
5 changes: 3 additions & 2 deletions Abstract/OperationsOnMaps/CostComposition.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 3 additions & 1 deletion Abstract/OperationsOnMaps/CostPartialSummation.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion Abstract/OperationsOnMaps/MapComposition.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
11 changes: 10 additions & 1 deletion Abstract/OperationsOnMaps/MapSummation.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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
Expand Down
4 changes: 2 additions & 2 deletions Abstract/Opti.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)))
Expand Down
2 changes: 1 addition & 1 deletion Cost/CostGoodRoughness.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
26 changes: 19 additions & 7 deletions Cost/CostHyperBolic.m
View file
Original file line number Diff line number Diff line change
@@ -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.
Expand Down Expand Up @@ -33,6 +33,7 @@
epsilon;
index;
sumOp;
sumEpsilon;
end

%% Constructor
Expand All @@ -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
Expand Down
Loading

0 comments on commit 95c7214

Please sign in to comment.