-
Notifications
You must be signed in to change notification settings - Fork 36
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
The v1.0.1 release does not contain any major changes but mainly incremental improvements and bug fixes. ###### LinOp - New LinOpBroadcastMatrix - fix pb in makeComposition of LinOpSum other cases in makeComposition btw LinOpDiag and LinOpConv - LinOpHess is now working for any dimensions, with the choice of the index along which the hessian is computed (Like in LinOpGrad). The method makeHtH return the correct corresponding convolution when the boundary condition is circular - Adapt MixNormSchatten to the new Hessian functionalities - makeHtH in LinOpGrad works now for any number of dimensions. It also returns the correct LinOpConv if the gradient is applied to specific dimensions - LinOpCpx is deprecated as it is not really a LinOp ###### Cost - New CostLinear - CostL2 now accept any linop as weighting (metric) operator - Indicator functions: Move the size parameter at the first place according to what was done for all maps ###### Opti - add methods computeCost and computeSNR in OutputOpti to facillitate deriving classes - Fix copy on write side effect in OptiVMLMB - Fix OptiVMLMB such that the returned value is the last accepted one - Add OutputOptiADMM which evaluates the Lagrangian instead of the cost function - OutputOpti can deal with both logical or scalar flags ###### Bug fixes other changes - Improve doc - Improve examples - Faster operation when memoize is off - Add Conjugate Gradient in Deconv_LS_NoReg01
- Loading branch information
Showing
45 changed files
with
1,562 additions
and
838 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,16 +1,16 @@ | ||
classdef LinOp < Map | ||
% Abstract class for linear operators | ||
% $$ \\mathrm{H}: \\mathrm{X}\\rightarrow \\mathrm{Y}.$$ | ||
% $$ \\mathrm{H}: \\mathrm{X}\\rightarrow \\mathrm{Y}.$$ | ||
% where \\(\\mathrm{X}\\) and \\(\\mathrm{Y}\\) are either | ||
% \\(\\mathbb{R}^N\\) or \\(\\mathbb{C}^N\\). | ||
% | ||
% All attributes of parent class :class:`Map` are inherited | ||
% All attributes of parent class :class:`Map` are inherited | ||
% | ||
% See also :class:`Map` | ||
|
||
%% Copyright (C) 2017 | ||
%% Copyright (C) 2017 | ||
% F. Soulez [email protected] | ||
% M. McCann [email protected] | ||
% M. McCann [email protected] | ||
% E. Soubies [email protected] | ||
% | ||
% This program is free software: you can redistribute it and/or modify | ||
|
@@ -26,23 +26,23 @@ | |
% You should have received a copy of the GNU General Public License | ||
% along with this program. If not, see <http://www.gnu.org/licenses/>. | ||
|
||
|
||
%% Constructor | ||
methods | ||
function this=LinOp() | ||
% Add new fields to memoizeOpts and memoCache | ||
this.memoizeOpts.applyAdjoint=false; | ||
this.memoizeOpts.applyAdjointInverse=false; | ||
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.applyAdjointInverse=struct('in', [], 'out', []); | ||
this.memoCache.applyAdjoint=struct('in', [], 'out', []); | ||
this.memoCache.applyHtH=struct('in', [], 'out', []); | ||
this.memoCache.applyHHt=struct('in', [], 'out', []); | ||
this.isDifferentiable=true; | ||
end | ||
end | ||
|
||
|
||
%% Interface Methods (cannot be overloaded in derived classes: Sealed) | ||
% In addition to inherited methods from Map | ||
|
@@ -62,9 +62,13 @@ | |
if ~checkSize(y, this.sizeout) % check input size | ||
error('Input to applyAdjoint was size [%s], didn''t match stated sizeout: [%s].',... | ||
num2str(size(y)), num2str(this.sizeout)); | ||
end | ||
end | ||
% memoize | ||
x = this.memoize('applyAdjoint', @this.applyAdjoint_, y); | ||
if this.memoizeOpts.applyAdjoint | ||
x = this.memoize('applyAdjoint', @this.applyAdjoint_, y); | ||
else | ||
x= this.applyAdjoint_(y); | ||
end | ||
% check output size | ||
if ~checkSize(x, this.sizein) | ||
warning('Output of applyAdjoint was size [%s], didn''t match stated sizein: [%s].',... | ||
|
@@ -78,9 +82,13 @@ | |
if ~checkSize(x, this.sizein) % check input size | ||
error('Input to applyHtH was size [%s], didn''t match stated sizein: [%s].',... | ||
num2str(size(x)), num2str(this.sizein)); | ||
end | ||
end | ||
% memoize | ||
y = this.memoize('applyHtH', @this.applyHtH_, x); | ||
if this.memoizeOpts.applyHtH | ||
y = this.memoize('applyHtH', @this.applyHtH_, x); | ||
else | ||
y= this.applyHtH_(x); | ||
end | ||
% check output size | ||
if ~checkSize(y, this.sizein) | ||
warning('Output of applyHtH was size [%s], didn''t match stated sizein: [%s].',... | ||
|
@@ -94,9 +102,13 @@ | |
if ~checkSize(y, this.sizeout) % check input size | ||
error('Input to applyHHt was size [%s], didn''t match stated sizeout: [%s].',... | ||
num2str(size(y)), num2str(this.sizeout)); | ||
end | ||
end | ||
% memoize | ||
x = this.memoize('applyHHt', @this.applyHHt_, y); | ||
if this.memoizeOpts.applyHHt | ||
x = this.memoize('applyHHt', @this.applyHHt_, y); | ||
else | ||
x =this.applyHHt_( y); | ||
end | ||
% check output size | ||
if ~checkSize(x, this.sizeout) | ||
warning('Output of applyHHt was size [%s], didn''t match stated sizeout: [%s].',... | ||
|
@@ -105,7 +117,7 @@ | |
end | ||
function L = transpose(this) | ||
% Returns a new :class:`LinOp` which is the Adjoint \\(\\mathrm{H}^{\\star}\\) of the | ||
% current \\(\\mathrm{H}\\). | ||
% current \\(\\mathrm{H}\\). | ||
if this.isComplexOut | ||
warning('Warning: Do you mean adjoint? For LinOp objects transpose() is an alias of adjoint method'); | ||
end | ||
|
@@ -114,7 +126,7 @@ | |
function L = ctranspose(this) | ||
% Do the same as :meth:`transpose` | ||
L = this.makeAdjoint_(); | ||
end | ||
end | ||
function L = makeAdjoint(this) | ||
% Do the same as :meth:`transpose` | ||
L = this.makeAdjoint_(); | ||
|
@@ -129,15 +141,19 @@ | |
num2str(size(x)), num2str(this.sizein)); | ||
end | ||
% memoize | ||
y = this.memoize('applyAdjointInverse', @this.applyAdjointInverse_, x); | ||
if this.memoizeOpts.applyAdjointInverse | ||
y = this.memoize('applyAdjointInverse', @this.applyAdjointInverse_, x); | ||
else | ||
y =this.applyAdjointInverse_(x); | ||
end | ||
% check output size | ||
if ~checkSize(y, this.sizeout) | ||
warning('Output of applyAdjointInverse was size [%s], didn''t match stated sizeout: [%s].',... | ||
num2str(size(y)), num2str(this.sizeout)); | ||
end | ||
end | ||
function M=makeHtH(this) | ||
% Compose the Adjoint Map \\(\\mathrm{H}^{\\star}\\) with | ||
% Compose the Adjoint Map \\(\\mathrm{H}^{\\star}\\) with | ||
% \\(\\mathrm{H}\\). Returns a new map \\(\\mathrm{M=H}^{\\star} \\mathrm{H}\\) | ||
% | ||
% Calls the method :meth:`makeHtH_` | ||
|
@@ -190,7 +206,7 @@ | |
% 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=[]; | ||
if isa(G,'LinOpSummation') | ||
|
@@ -212,7 +228,7 @@ | |
if isempty(M) | ||
M=LinOpSummation({this,G},[1,1]); | ||
end | ||
else | ||
else | ||
M = MapSummation({this,G},[1,1]); | ||
end | ||
end | ||
|
@@ -232,15 +248,15 @@ | |
M=LinOpComposition(this,this'); | ||
end | ||
function M = makeInversion_(this) | ||
% Constructs a :class:`LinOpInversion` corresponding to | ||
% 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` | ||
% object to compose the current LinOp (this) with the given :class:`LinOp`\\(\\mathrm{G}\\). | ||
% object to compose the current LinOp (this) with the given :class:`LinOp`\\(\\mathrm{G}\\). | ||
% Otherwise the composition will be a :class:`MapComposition`. | ||
if isa(G,'LinOp') | ||
if isa(G,'LinOpDiag') && G.isScaledIdentity | ||
|
@@ -249,17 +265,17 @@ | |
else | ||
M=LinOpDiag(this.sizeout,G.diag)*this; % To always have the diag Op in H1 for LinOpComposition | ||
end | ||
% is HHt or HtH | ||
% is HHt or HtH | ||
elseif (isa(G,'LinOpAdjoint') && isequal(G.TLinOp,this)) || (isa(this,'LinOpAdjoint') && isequal(this.TLinOp,G)) | ||
M=this.makeHHt(); | ||
% composition with a sum of LinOps | ||
% composition with a sum of LinOps | ||
elseif isa(G,'LinOpSummation') | ||
M=G.alpha(1)*this*G.mapsCell{1}; | ||
for i=2:G.numMaps | ||
M=M+G.alpha(i)*this*G.mapsCell{i}; | ||
end | ||
% composition with a LinOpComposition | ||
elseif isa(G,'LinOpComposition') | ||
% composition with a LinOpComposition | ||
elseif isa(G,'LinOpComposition') | ||
% to handle properly scalar multiplications | ||
if isa(G.H1,'LinOpDiag') && G.H1.isScaledIdentity | ||
if isa(this,'LinOpDiag') && this.isScaledIdentity | ||
|
@@ -270,7 +286,7 @@ | |
else | ||
M=this*G.H1*G.H2; | ||
end | ||
% composition with other LinOps | ||
% composition with other LinOps | ||
else | ||
M = LinOpComposition(this,G); | ||
end | ||
|
@@ -289,22 +305,6 @@ | |
x = this.applyAdjoint(y); | ||
end | ||
end | ||
|
||
% TODO: ask ferreol: is this needed? when do we have fast HT W H? | ||
% function y = HtWH(this,x,W) | ||
% % Apply \\(\\mathrm{H}^*\\mathrm{WH}\\) | ||
% % | ||
% % :param x: \\(\\in X\\) | ||
% % :param W: a :class:`LinOp` object | ||
% % :returns y: \\(= \\mathrm{H^*WHx}\\) | ||
% % | ||
% | ||
% if (isscalar(W) && isreal(W)) | ||
% y = W.*this.HtH(x); | ||
% else | ||
% assert(isa(W,'LinOp'),'W must be a LinOp'); | ||
% y = this.adjoint(W.apply(this.apply(x))); | ||
% end | ||
% end | ||
|
||
end | ||
|
Oops, something went wrong.