biomedical-imaging-group / globalbioim Goto Github PK

There is still two classes (OneToMany and StackLinOp) that should be properly integrated in the new framework (v1.0 milestone). We should think to which extend it is possible to have generic broadcasting and concatenation tools for LinOps.

Hello GBI Team,

thanks for creating and maintaining this neat Matlab library!

I was facing some issues regarding the adjoint calculation in applyJacobianT_ of StackMap:

function x = applyJacobianT_(this, y, v)
    % Reimplemented from :class:`Map`   
    x = cell(this.numMaps,1);
    for n = 1:this.numMaps
        x{n} = this.alpha(n) .* this.mapsCell{n}.applyJacobianT(y,v);
    end
    x = cat(length(this.sizeout),x{:});
end

• Size related run-time failure.
• Functional issue: the forward pass at its core is broadcast -> per-channel-op-application -> concat; so the backward pass should be slice -> per-channel-op-transpose-application -> sum-reduction.

function x = applyJacobianT_(this, y, v)
    % Reimplemented from :class:`Map`
    
    % The adjoint op of a StackMap (:= apply & concat) is slice &
    % apply & summation.
    % Therefore we need to slice the upstream gradient along its
    % last dimension. To be working with arbitrary dimension, we
    % could e.g. use approaches like:
    %     https://stackoverflow.com/questions/19955653/matlab-last-dimension-access-on-ndimensions-matrix
    % For now, we take the simple path just flatten it to 2-D,
    % followed by another reshape op back to basic.

    % TODO(cl): quick'n'dirty, to be tested for robustness!
    
    full_dim = size(y);
    plain_dim = full_dim(1:end - 1);
    cat_dim = full_dim(end);
    
    y_flattened = reshape(y, [], cat_dim);
    
    x = cell(this.numMaps,1);

    for n = 1:this.numMaps
        y_n = reshape(y_flattened(:, n), plain_dim);
        x{n} = this.alpha(n) .* this.mapsCell{n}.applyJacobianT(y_n,v);
    end

    x = cat(length(this.sizeout),x{:});
    x = sum(x, length(this.sizeout));
end

What's your thoughts on a preferred solution? Let me know! Of course, this needs testing and cleaning-up...

Somewhat related question. Calculating the adjoint of StackMap * UpstreamOp, the adjoint of UpsteamOp gets executed multiple times, once for each channel of StackMap. Any better solution than wrapping StackMap functionality to first execute and sum all channel contributions before passing the gradient to the upstream node?

Best,
CL

Would be nice to make GlobalBioIm compatible with Octave.

Should we include a default call to the function estimateNorm as soon as an access to the norm field of a LinOp is performed ? Or this should be let to the user ?

Hi,

I got an error when I run one of the examples of Deconv_Ls_NonNeg_NoReg.m, when Installing OptimPackLegacy, the below error came out:

Building with 'Microsoft Visual C++ 2017 (C)'.
Error using mex
opl_vmlmb.c
P:\work\research\toolbox\GlobalBioIm\Opti\OptiUtils\MatlabOptimPack\OptimPackLegacy-master\src\opl_vmlmb.c(100):
error C2059: syntax error: ''
P:\work\research\toolbox\GlobalBioIm\Opti\OptiUtils\MatlabOptimPack\OptimPackLegacy-master\src\opl_vmlmb.c(104):
error C2059: syntax error: ''

Error in installOptimPack (line 37)
eval(['mex ',matDir,'m_opl_vmlmb_get_reason.c ', CFiles,MexOpt]);
Error in OptiVMLMB (line 97)
installOptimPack();
Error in Deconv_LS_NonNeg_NoReg (line 68)
VMLMB=OptiVMLMB(F,0.,[]);

I'm using Win10 and Matlab 2019 a with 'Microsoft Visual C++ 2017 (C)'. Hope this can be resolved soon.

Thanks.

I encountered a bug when trying to sum two "LinOpConv" operators whose dimensions along which the convolution is performed do not match. I'm attaching a code snippet illustrating the problem. I don't think it is possible to rewrite this sum as a single convolution operator in the general case, so an easy fix would be to change the "plus_" method of LinOpConv accordingly (from testing isa(G,'LinOpConv') to something like isa(G,'LinOpConv') && isequal(sort(this.index), sort(G.index))). However I'm not sure I fully understand all the subcases of the "plus_" method so I'm not confortable trying to fix the problem completely myself

f.

Some examples of stuff that have to be improved:

• In Map plus_ and minus_ test if equal to itself (or a scaled version of itself) and return the corresponding scaled version of itself. Do the analogy in LinOp and Cost.
• In Map makeComposition_ test if inverse and return identity (the analogy of makeComposition_ in MapInversion)
• For LinOpConv if we do H^-1*H we get a LinOpConv and not an identity. In makeComposition_ of LinOpConv test if the new mtf is constant and return the appropriate LinOpDiag. This is more general and holds other LinOps that have a new implementation of mpower_
• add the possibility to add a scalar to a Map for adding a scaled identity

After discussing with Michael it seems it would be convenient to acknowledge CIBM by adding a logo somewhere. I think the most natural would be to also add an EPFL one. If you agree I will go ahead and make a PR with respect to that. Thanks ;-)

Hi guys,

Do you have the plan to implement the TwIST algorithm (Two-step Iterative Shrinkage/Thresholding Algorithm for Linear Inverse Problems) ?
Since it has been proved very efficient in compressive sensing.

Thanks.

Python 2.7 will soon no longer be supported on travis.
Our doc is stuck with python 2.7 and sphinx 1.8. and sphinxcontrib-matlabdomain 0.7 (see requirements.txt )
Pull request #5 #6 #7 are only temporary fix