Skip to content

Commit

Permalink
Reading for v1.0 release
Browse files Browse the repository at this point in the history 
Merging master branch and squashing all commits

Major Changes and improvements
=============
- Remove deprecated `LinOpScaling`,
- Improve tests with `CheckLinOp`, `CheckMap` and `TestsSummationLinOp' functions
- implement methods `makeAdjoint` and `makeInversion` in `LinOp`
- more simplifications in `Summation`, `Composition`, `Ìnversion` for `Map` and hence `LinOp`and `Cost`
- New `LinOpBroadcast` linear operator
- Add `ElementWiseOp` non linear operators: `OpEWSquareRoot`, `OpEWInverse` and `OpEWSSquaredMagnitude`.
- New `CostGoodRoughness` cost function
- Introduce positivity constraint in `CostL1`
- New  `CostTV` cost function
- Rework `Opti` optimization class
- Introduce several `TestCvg` classes to assess convergence of optimization according to different criteria
- Automatic download and compilation of OptimPackLegacy for `OptiVMLMB`
- `OptiConjGrad`no longer takes a `W` parameter
- New ÒptiFGP` (Fast gradient proximal) optimization function

Bug fixes and minor changes
------------------
- Improve doc
- Implement function `cmpSize` to compare size of vectors
- Faster `CostHyperBolic` cost function- Automatic compilation of `C` files needed for `CostMixNormSchatt1` cost function
- Better examples
- Add functions `GenerateData` and `GenerateData` to generate  star phantom- More options in `LinOpConv'
- Implement `make...` methods in LinOpDFT
- Reworked `OutputOpti` printing class

Contributors
-----------------
This release contains contributions from Mike Mc Cann, Thomas Debarre, Thanh-an Pham, Emmanuel Soubies and Ferréol Soulez
  • Loading branch information
ferreol committed Mar 16, 2018
1 parent c798da0 commit 6a5447a
Show file tree
Hide file tree
Showing 143 changed files with 4,663 additions and 4,294 deletions.
1 change: 1 addition & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,4 +2,5 @@
*.mex*
*html/
*doctrees/
*Opti/OptiUtils/MatlabOptimPack/
.DS_Store
16 changes: 12 additions & 4 deletions Abstract/Cost.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,10 @@
y=0; % data y
end

properties (SetAccess = protected,GetAccess = protected)
recProxProxFench=0; % boolean to control infinite recursion between Prox and ProxFench when not implemented in subclasses
end

%% Constructor
methods
function this=Cost(sz,y)
Expand Down Expand Up @@ -135,21 +139,25 @@
% at \\(\\mathrm{z} \\in \\mathrm{X} \\) using the
% Moreau's identity which uses the :meth:`applyProxFench` method
% $$\\mathrm{prox}_{\\alpha C}(\\mathrm{z}) = \\mathrm{z} - \\alpha \\,\\mathrm{prox}_{\\frac{1}{\\alpha}C^*}\\left(\\frac{\\mathrm{z}}{\\alpha}\\right).$$
if this.isConvex
if this.isConvex && ~this.recProxProxFench
this.recProxProxFench=1;
x= z - alpha*(this.applyProxFench(z/alpha,1/alpha));
this.recProxProxFench=0;
else
error('applyProx_ method not implemented');
error('applyProx_ and applyProxFench_ methods not implemented');
end
end
function y=applyProxFench_(this,z,alpha)
% By default, if the cost \\(C\\) :attr:`isConvex`, computes the proximity operator of the Fenchel Transform
% \\(C^*\\) at \\(\\mathrm{z} \\in \\mathrm{Y} \\) using the
% Moreau's identity which uses the :meth:`applyProx` method
% $$\\mathrm{prox}_{\\alpha C^*}(\\mathrm{z}) = \\mathrm{z} - \\alpha \\,\\mathrm{prox}_{\\frac{1}{\\alpha}C}\\left(\\frac{\\mathrm{z}}{\\alpha}\\right).$$
if this.isConvex
if this.isConvex && ~this.recProxProxFench
this.recProxProxFench=1;
y= z - alpha*(this.applyProx(z/alpha,1/alpha));
this.recProxProxFench=0;
else
error('applyProxFench_ method not implemented');
error('applyProx_ and applyProxFench_ methods not implemented');
end
end
function M = plus_(this,G)
Expand Down
81 changes: 42 additions & 39 deletions Abstract/LinOp.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,8 +32,10 @@
function this=LinOp()
% Add new fields to memoizeOpts and memoCache
this.memoizeOpts.applyAdjoint=false;
this.memoizeOpts.applyAdjointInverse=false;
this.memoizeOpts.applyHtH=false;
this.memoizeOpts.applyHHt=false;
this.memoCache.applyAdjointInverse=struct('in', [], 'out', []);
this.memoCache.applyAdjoint=struct('in', [], 'out', []);
this.memoCache.applyHtH=struct('in', [], 'out', []);
this.memoCache.applyHHt=struct('in', [], 'out', []);
Expand Down Expand Up @@ -112,9 +114,13 @@
function L = ctranspose(this)
% Do the same as :meth:`transpose`
L = this.makeAdjoint_();
end
function L = makeAdjoint(this)
% Do the same as :meth:`transpose`
L = this.makeAdjoint_();
end
function y= applyAdjointInverse(this,x)
% Computes \\(\\mathrm{y} = \\mathrm{H}^{-\star} \\mathrm{x}\\) for the given
% Computes \\(\\mathrm{y} = \\mathrm{H}^{-\\star} \\mathrm{x}\\) for the given
% \\(\\mathrm{x} \\in \\mathrm{X}\\). (if applicable)
%
% Calls the method :meth:`applyAdjointInverse_`
Expand All @@ -129,18 +135,6 @@
warning('Output of applyAdjointInverse was size [%s], didn''t match stated sizeout: [%s].',...
num2str(size(y)), num2str(this.sizeout));
end

% **(Abstract method)** Apply \\(\\mathrm{H}^{-*}\\) (if applicable)
%
% :param x: \\(\\in X\\)
% :returns y: \\(= \\mathrm{H^{-*}x}\\)
%

if this.isinvertible
error('adjointInverse not implemented');
else
error('Operator not invertible');
end
end
function M=makeHtH(this)
% Compose the Adjoint Map \\(\\mathrm{H}^{\\star}\\) with
Expand Down Expand Up @@ -168,9 +162,10 @@
% - makeAdjoint_(this)
% - makeHtH_(this)
% - makeHHt_(this)
% - makeInversion_(this)
% - makeComposition_(this, G)
methods (Access = protected) % all the other underscore methods
function x = applyAdjoint_(this, y)
function x = applyAdjoint_(~, ~)
% Not implemented in this Abstract class
error('applyAdjoint_ method not implemented');
end
Expand All @@ -187,45 +182,61 @@
% way to perform computation.
x = this.apply(this.applyAdjoint(y));
end
function y = applyAdjointInverse_(this,x)
function y = applyAdjointInverse_(~,~)
% Not implemented in this Abstract class
error('applyAdjointInverse_ method not implemented');
end
function M = plus_(this,G)
% If \\(\\mathrm{G}\\) is a :class:`LinOp`, constructs a :class:`LinOpSummation` object to sum the
% current :class:`LinOp` \\(\\mathrm{H}\\) with the given \\(\\mathrm{G}\\).
% Otherwise the summation will be a :class:`MapSummation`.

if isa(G,'LinOp')
M = LinOpSummation({this,G},[1,1]);
else
M=[];
if isa(G,'LinOpSummation')
% Find elements of the same type
ind=find(strcmp(this.name,cellfun(@(T) T.name,G.mapsCell,'UniformOutput',false)));
if length(ind)==1
% To avoid infinite loops (normally it should never goes in the else because the sum of
% two LinOp of the same type can always be simplified. If not the sum_ method of the corresponding
% LinOp has to be implemented properly).

M=this + G.alpha(1)*G.mapsCell{ind};
for ii=1:G.numMaps
if ii~=ind
M= M+G.alpha(ii)*G.mapsCell{ii};
end
end
end
end
if isempty(M)
M=LinOpSummation({this,G},[1,1]);
end
else
M = MapSummation({this,G},[1,1]);
end
end
function M = minus_(this,G)
% If \\(\\mathrm{G}\\) is a :class:`LinOp`, constructs a :class:`LinOpSummation` object to subtract to the
% current :class:`LinOp` \\(\\mathrm{H}\\), the given \\(\\mathrm{G}\\).
% Otherwise the summation will be a :class:`MapSummation`.
if isa(G,'LinOp')
M = LinOpSummation({this,G},[1,-1]);
else
M = MapSummation({this,G},[1,-1]);
end
end
function M = makeAdjoint_(this)
% Constructs a :class:`LinOpAdjoint` from the current
% current :class:`LinOp` \\(\\mathrm{H}\\)
% current :class:`LinOp` \\(\\mathrm{H}\\)
M=LinOpAdjoint(this);
end
function M = makeHtH_(this)
% Constructs a :class:`LinOpComposition` corresponding to
% Constructs a :class:`LinOpComposition` corresponding to
% \\(\\mathrm{H}^{\\star}\\mathrm{H}\\)
M=LinOpComposition(this',this);
end
function M = makeHHt_(this)
% Constructs a :class:`LinOpComposition` corresponding to
% Constructs a :class:`LinOpComposition` corresponding to
% \\(\\mathrm{H}\\mathrm{H}^{\\star}\\)
M=LinOpComposition(this,this');
end
function M = makeInversion_(this)
% Constructs a :class:`LinOpInversion` corresponding to
% \\(\\mathrm{H}^{-1}\\)

M=LinOpInversion(this);
end
function M = makeComposition_(this, G)
% Reimplemented from parent class :class:`Map`.
% If \\(\\mathrm{G}\\) is a :class:`LinOp`, constructs a :class:`LinOpComposition`
Expand Down Expand Up @@ -267,20 +278,12 @@
M = makeComposition_@Map(this,G);
end
end
function M = mpower_(this,p)
% Reimplemented from :class:`Map`
if p==-1
M=LinOpInversion(this);
else
M=mpower_@Map(this,p);
end
end
end

%% Methods of superclass Map that do not need to be reimplemented in derived Costs
% - applyJacobianT_(this, y, v)
methods (Access = protected, Sealed)
function x = applyJacobianT_(this, y, v)
function x = applyJacobianT_(this, y, ~)
% Uses the method applyAdjoint (hence do not need to be
% reimplemented in derived classes)
x = this.applyAdjoint(y);
Expand Down
64 changes: 50 additions & 14 deletions Abstract/Map.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,7 @@
% - minus(this,G)
% - mpower(this,p)
% - mtimes(this,G)
% - times(this,G)
% - size(this, [dim])
methods (Sealed)
function y = apply(this, x)
Expand Down Expand Up @@ -135,7 +136,7 @@
% \\(\\mathrm{G}\\). Returns a new map \\(\\mathrm{M=HG}\\)
%
% Calls the method :meth:`makeComposition_`
if ~isequal(this.sizein, G.sizeout)
if ~cmpSize(this.sizein,G.sizeout)
error('Input to makeComposition is a %s of sizeout size [%s], which didn''t match the %s sizein [%s].',...
class(G),num2str(G.sizeout),class(this), num2str(this.sizein));
end
Expand All @@ -146,11 +147,11 @@
% $$ \\mathrm{M}(\\mathrm{x}) := \\mathrm{H}(\\mathrm{x}) + \\mathrm{G}(\\mathrm{x})$$
%
% Calls the method :meth:`plus_`
if ~isequal(this.sizein, G.sizein) % check input size
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 ~isequal(this.sizeout, G.sizeout) % check input size
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
Expand All @@ -161,11 +162,11 @@
% $$ \\mathrm{M}(\\mathrm{x}) := \\mathrm{H}(\\mathrm{x}) - \\mathrm{G}(\\mathrm{x})$$
%
% Calls the method :meth:`minus_`
if ~isequal(this.sizein, G.sizein) % check input size
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 ~isequal(this.sizeout, G.sizeout) % check input size
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
Expand All @@ -174,6 +175,8 @@
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);
end
function M = mtimes(this,G)
Expand All @@ -183,12 +186,16 @@
% - 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 isa(G,'Cost')
M=CostMultiplication(this,G);
if (isnumeric(this) && isscalar(this)) % Left multiplication by scalar
if ~isequal(this,1) % if multiply by 1 do noting
if isa(G,'Cost')
M=CostMultiplication(this,G);
else
H=LinOpDiag(G.sizeout,this);
M=H.makeComposition(G);
end
else
H=LinOpDiag(G.sizeout,this);
M=H.makeComposition(G);
M=G;
end
else
M =this.makeComposition(G);
Expand All @@ -197,7 +204,22 @@
M = this.apply(G);
end
end
% TODO Overload operator .* to multiply term by term two Maps ? (avoid the use of LinOpDiag)
function M = times(this,G)
% Returns a new :class:`Map` which is the element-wise multiplication of the
% current \\(\\mathrm{H}\\) with \\(\\mathrm{G}\\)
% $$ \\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));
end
M=this.times_(G);
end
% TODO Overload operator .^ to do (H(x)).^p ?
function sz = size(this, varargin)
sz = {this.sizeout, this.sizein};
Expand All @@ -216,6 +238,8 @@
% - plus_(this,G)
% - minus_(this,G)
% - mpower_(this,p)
% - times_(this,G)
% - makeInversion_(this)
% - makeComposition_(this, H)
methods (Access = protected)
function y = apply_(this, x)
Expand All @@ -238,26 +262,38 @@
function M = minus_(this,G)
% Constructs a :class:`MapSummation` object to subtract to the
% current :class:`Map` \\(\\mathrm{H}\\), the given \\(\\mathrm{G}\\).
M = MapSummation({this,G},[1,-1]);
M = this + (-1)*G;
end
function M = mpower_(this,p)
% When \\(p=-1\\), constructs a :class:`MapInversion` object which
% is the inverse Map of \\(\\mathrm{H}\\).
% When \\(p\\neq-1\\), this method is not implemented in this Abstract class
if p==-1
M=MapInversion(this);
M=this.makeInversion_();
else
error('mpower_ method not implemented for the given power==%d',p);
end
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}\\).

M=MapMultiplication(this,G);
end
function M = makeInversion_(this)
% Constructs a :class:`MapInversion` corresponding to
% \\(\\mathrm{H}^{-1}\\)

M=MapInversion(this);
end
function M = makeComposition_(this, G)
% Constructs a :class:`MapComposition` object to compose the
% current Map \\(\\mathrm{H}\\) with the given \\(\\mathrm{G}\\).

if isa(G,'MapInversion') && isequal(G.M,this)
M=LinOpScaledIdentity(this.sizein,1);
elseif isa(G,'MapComposition')
M=this*G.H1*G.H2;
M=(this*G.H1)*G.H2;
else
M = MapComposition(this,G);
end
Expand Down
11 changes: 10 additions & 1 deletion Abstract/OperationsOnMaps/CostComposition.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,16 @@
end
function M = makeComposition_(this,G)
% Reimplemented from :class:`Cost`
M=CostComposition(this.H1,this.H2*G);
if isa(G,'LinOp')
T=G*G';
if isa(T,'LinOpDiag') && T.isScaledIdentity
M=CostComposition(this,G);
else
M=CostComposition(this.H1,this.H2*G);
end
else
M=CostComposition(this.H1,this.H2*G);
end
end
end

Expand Down
Loading

0 comments on commit 6a5447a

Please sign in to comment.