Got rid of modules by using private functions instead

.gitmodules

-[submodule "+LRTC"]
-	path = +LRTC
-	url = git@gitlab.inria.fr:edib/lrtc.git
\ No newline at end of file

LRTC @ 3f7710a6

-Subproject commit 3f7710a61341949178d3d56acec2c67ec013248e

private/Fold.m

+function [X] = Fold(X, dim, i)
+dim = circshift(dim, [1-i, 1-i]);
+X = shiftdim(reshape(X, dim), length(dim)+1-i);
\ No newline at end of file

private/LRTCADMrho.m

+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% ADM algorithm: tensor completion
+% paper: Tensor completion for estimating missing values in visual data
+% date: 05-22-2011
+% min_X: \sum_i \alpha_i \|X_(i)\|_*
+% s.t.:  X_\Omega = T_\Omega
+%Note MLP : corresponds to algorithm 4 (HaLRTC) in the article
+%(i.e. using augmented lagrangian method and alternating direction).
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%function [X_final, X, M, rank, errList, OmegaPrime] = LRTCADMrho(T, Omega, alpha, betaMult, maxIter, epsilon, TraceNormToRankParam, CorruptCols,AreaMask, maxRank, X)
+function [X_final, X, M, rank, errList,OmegaPrime] = LRTCADMrho(T, Omega, alpha, betaMult, maxIter, epsilon, TraceNormToRankParam, CorruptCols, maxRank, X)
+
+DoCenter = false; %center the data before completion.
+legacyThreshold = false; %true  : threshold value determined for nuclear norm minimization.
+                         %false : threshold value determined for rank minimization.
+                         %This parameter is automatically set to true when TraceNormToRankParam=0 (nuclear norm minimization).
+l0Reg = false;
+
+
+Omega = logical(Omega);
+meanT = sum(T(:).*Omega(:)) / sum(Omega(:));
+
+if(DoCenter)
+    OmegaNan = double(Omega); OmegaNan(~Omega) = nan;
+    meanTCol = nanmean(T.*OmegaNan,2);
+    clear OmegaNan
+    T = bsxfun(@minus,T,meanTCol);
+    
+    normT2 = sum((T(:).*Omega(:)).^2);
+else
+    normT2 = sum(((T(:)-meanT).*Omega(:)).^2);
+end
+
+if(l0Reg)
+    numKnown=sum(Omega(:));
+end
+%numKnownCorrupt=sum(vec(Omega(:,CorruptCols,:)));
+%maxVal2 = 255^2;
+
+if(isscalar(epsilon))
+    epsilonStop = epsilon;
+    epsilonNoise = epsilon;
+else
+    epsilonStop = epsilon(1);
+    epsilonNoise = epsilon(2);
+end
+if(epsilonNoise==0||~exist('CorruptCols','var')||isempty(CorruptCols)||~any(ismember(CorruptCols,1:size(T,2))))
+    NoisyData=false;
+    if(l0Reg)
+        CorruptCols = intersect(CorruptCols,1:size(T,2));
+        NoCorruptCols = setdiff(1:size(T,2),CorruptCols);
+    else
+        CorruptCols = [];
+        NoCorruptCols = 1:size(T,2);
+    end
+else
+    NoisyData=true;
+    CorruptCols = intersect(CorruptCols,1:size(T,2));
+    NoCorruptCols = setdiff(1:size(T,2),CorruptCols);
+    if(isempty(NoCorruptCols))
+        AllCorrupt = true;
+        normTCorrupt2 = normT2;
+    else
+        AllCorrupt = false;
+        if(DoCenter)
+            normTCorrupt2 = sum(vec((T(:,CorruptCols,:).*Omega(:,CorruptCols,:))).^2);
+        else
+            meanTCorrupt = sum(vec((T(:,CorruptCols,:).*Omega(:,CorruptCols,:)))) / sum(vec(Omega(:,CorruptCols,:)));
+            normTCorrupt2 = sum(vec(((T(:,CorruptCols,:)-meanTCorrupt).*Omega(:,CorruptCols,:))).^2);
+        end
+    end
+end
+
+if (~exist('X','var')||isempty(X))
+    X = T;
+    if(DoCenter)
+        X(~Omega) = 0;
+    else
+        X(~Omega) = meanT;%X(logical(1-Omega)) = 256;
+    end
+    %OmegaNan=double(Omega);OmegaNan(~Omega)=nan;
+    %X = T.*Omega + (1-Omega).*repmat(nanmean(T.*OmegaNan),size(T,1),1);
+    %clear OmegaNan
+end
+if (~exist('TraceNormToRankParam','var')|| isempty(TraceNormToRankParam))
+    TraceNormToRankParam=0;
+end
+if(~exist('maxRank','var')||isempty(maxRank))
+    maxRank=inf;
+end
+
+if(TraceNormToRankParam==0)
+    legacyThreshold = true;
+end
+
+
+%Check the valid dimensions for tensor unfolding
+TensorDim = ndims(T);
+if(TensorDim==2)
+    TensorDim=1;
+    alpha=1;
+    ValidDims=1;
+    numValidDims=1;
+else
+    ValidDims = 1:TensorDim; ValidDims(alpha(ValidDims)==0)=[];
+    numValidDims = length(ValidDims);
+end
+
+dim = size(T);
+M = cell(TensorDim, 1);
+Y = M;
+rank = zeros(TensorDim,1);
+beta = zeros(1,TensorDim);
+if(l0Reg)
+    OmegaPrime = Omega;
+else
+    OmegaPrime=[];
+end
+
+trMult=@(x)x'*x;
+for i = ValidDims
+    M{i} = X;
+    Y{i} = zeros(dim, class(T));
+    %beta(i) = 1.25 / lansvd(Unfold(X,dim,i), 1, 'L');
+    SV12 = sqrt(svds(trMult(Unfold(X,dim,i)), 2));
+    if(legacyThreshold)
+        beta(i) = 2 / (SV12(1)+SV12(2));    %initial threshold = 1/beta = mean(SV1+SV2);
+        beta(i) = beta(i)*(alpha(i).^2);
+    else
+        beta(i) = 8 / (SV12(1)+SV12(2)).^2; %initial threshold = sqrt(2/beta) = mean(SV1+SV2);
+        beta(i) = beta(i)*alpha(i);
+    end
+end
+%beta(ValidDims) = max(beta)*ones(1,numValidDims);
+
+
+Ysum = zeros(dim, class(T));
+Msum = zeros(dim, class(T));
+errList=[];
+fprintf('Iteration:                  ');
+for k = 1: maxIter
+    
+    % update M
+    Ysum = 0*Ysum;
+    Msum = 0*Msum;
+    for i = ValidDims
+        if(legacyThreshold)
+            [M{i}, rank(i), ~] = Pro2TraceNorm(Unfold(X-Y{i}/beta(i), dim, i), alpha(i)/beta(i), TraceNormToRankParam);
+        else
+            [M{i}, rank(i), ~] = Pro2TraceNorm(Unfold(X-Y{i}/beta(i), dim, i), alpha(i)*sqrt(2/beta(i)), TraceNormToRankParam);
+        end
+        M{i} = Fold(M{i}, dim, i);
+        Ysum = Ysum + Y{i};
+        Msum = Msum + beta(i)*M{i};
+        fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%5i (rank =%4i)',k,rank(i));
+    end
+    
+    if(any(rank > maxRank))
+        break;
+    end
+    
+    % update X (outside of Omega)
+    X = (Msum + Ysum) / sum(beta);
+    
+    
+    % compute the error and check stop criterion
+    if(l0Reg)
+        [errList(end+1),errCorrupt2] = SQErrors(X,T,OmegaPrime,CorruptCols);
+    else
+        [errList(end+1),errCorrupt2] = SQErrors(X,T,Omega,CorruptCols);
+    end
+    errList(end) = errList(end) / normT2;
+    %errList(end+1) = norm((X-T).*Omega,'fro')^2 / normT2;
+    %errList(end+1) = sum(((X(:)-T(:)).*Omega(:)).^2) / normT2;%/ numKnown / maxVal2;
+    if (errList(end) < epsilonStop || k == maxIter )
+    %if k == maxIter
+        X_final = X;
+        X_final(Omega) = T(Omega);
+        if(DoCenter)
+            X_final = bsxfun(@plus,X_final,meanTCol);
+        end
+        fprintf('\n');
+        break;
+    end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Update X (in Omega)
+if(NoisyData) %l2 tolerance:
+    % Allow a small error epsilon (in l2 norm) on the known elements of the matrix.
+    if(AllCorrupt)
+        lambda = max(sqrt( errList(end) / (epsilonNoise) ) - 1, 0);
+        X = X.*(1-Omega) + ((X + lambda*T)/(1+lambda)).*Omega;
+        %X(Omega) = ( X(Omega) + lambda * T(Omega) ) / (1 + lambda); %Much slower by indexing with Omega.
+    else
+        errCorrupt2 = errCorrupt2 / normTCorrupt2;
+        %errCorrupt2 = norm((X(:,CorruptCols,:)-T(:,CorruptCols,:)).*Omega(:,CorruptCols,:),'fro')^2 / normTCorrupt2;
+        %errCorrupt2 = sum(vec((X(:,CorruptCols,:)-T(:,CorruptCols,:)).*Omega(:,CorruptCols,:)).^2)/ normTCorrupt2;%numKnownCorrupt / maxVal2;
+        lambda = max(sqrt( errCorrupt2 / (epsilonNoise) ) - 1, 0);
+
+        UpdateXOmega(X,T,Omega,lambda,CorruptCols); %faster with mex.
+        %X(:,CorruptCols,:) = X(:,CorruptCols,:).*(1-Omega(:,CorruptCols,:)) + ((X(:,CorruptCols,:) + lambda*T(:,CorruptCols,:))/(1+lambda)).*Omega(:,CorruptCols,:);
+        %X(:,NoCorruptCols,:) = X(:,NoCorruptCols,:).*(1-Omega(:,NoCorruptCols,:)) + T(:,NoCorruptCols,:).*Omega(:,NoCorruptCols,:);
+    end
+else
+    
+    if(l0Reg)
+        
+    OmegaPrime = Omega;
+    
+    %l1 tolerance:
+    %{
+    X0 = Omega(:,CorruptCols,:).*(X(:,CorruptCols,:)-T(:,CorruptCols,:));
+    lambda_l1 = SoftThresholdNum(sort(abs(X0(:)),'descend'),5e5);%3e8
+    
+    X(:,CorruptCols,:) = (1-Omega(:,CorruptCols,:)).*X(:,CorruptCols,:) + Omega(:,CorruptCols,:).*( T(:,CorruptCols,:) + max(abs(X0)-lambda_l1,0).*(2*(X0>0)-1) );
+    X(:,NoCorruptCols,:) = (1-Omega(:,NoCorruptCols,:)).*X(:,NoCorruptCols,:) + Omega(:,NoCorruptCols,:).*T(:,NoCorruptCols,:);
+    %}
+    
+    
+    %l0 tolerance:
+    %%{
+    %if(any(rank>1))
+    OmegaPrimeId=Omega(:,CorruptCols,:).*(abs(X(:,CorruptCols,:)-T(:,CorruptCols,:)));
+    [OmegaPrime0,OmegaPrimeId] = sort(OmegaPrimeId(:),'descend');
+    
+    numTh = HardThresholdNum(OmegaPrime0,3e6);%2e5%1e8
+    
+    OmegaPrime0 = OmegaPrime(:,CorruptCols,:);
+    OmegaPrime0(OmegaPrimeId(1:numTh))=false;
+    %OmegaPrime0(OmegaPrimeId(1:round(numKnown*max(rank)/min(size(T)))))=false;
+    %OmegaPrime0(OmegaPrimeId(1:round(numKnown/180)))=false;
+    
+    OmegaPrime(:,CorruptCols,:) = OmegaPrime0;
+    %end
+    
+    X = X.*(1-OmegaPrime) + T.*OmegaPrime;
+    
+    %}
+    else
+        X = X.*(1-Omega) + T.*Omega;
+        %X(Omega) = T(Omega);%Much slower by indexing with Omega.
+    end
+    
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % update Y
+    for i = ValidDims
+        Y{i} = Y{i} + beta(i)*(M{i} - X);
+    end
+    
+    % update beta
+    beta  = beta  * betaMult;
+    
+end
+

private/Pro2TraceNorm.m

+function [X, n, Sigma2] = Pro2TraceNorm(Z, tau, TraceNormToRankParam)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% min: 1/2*||Z-X||^2 + ||X||_tr
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% [S, V, D, Sigma2] = MySVDtau(Z, tau);
+% V = max(diag(V) - tau, 0);
+% n = sum(V > 0);
+% X = S(:, 1:n) * diag(V(1:n)) * D(:, 1:n)';
+
+%TraceNormToRankParam : - Set to zero (default) for trace norm minimization (soft thresholding of singular values).
+%                       - Set to inf for actual rank minimization (hard thresholding of singular values).
+%                       - Other positive values should give intermediate results (keeping some high frequencies in the final matrix).
+
+%% new
+[m, n] = size(Z);
+
+if(nargin<3)
+    TraceNormToRankParam=0;
+end
+
+if 2*m < n
+    AAT = Z*Z';
+    [S, Sigma2, D] = svd(AAT);
+    Sigma2 = diag(Sigma2);
+    V = sqrt(Sigma2);
+    tol = max(size(Z)) * eps(max(V));
+    n = sum(V > max(tol, tau));
+    
+    if(isinf(TraceNormToRankParam))
+        sub = zeros(1,n);
+    else
+        sub = ([0 [1:n-1]/(n-1)].^TraceNormToRankParam)*tau;
+    end
+    
+    mid = max(V(1:n)-sub', 0) ./ V(1:n) ;
+    X = S(:, 1:n) * diag(mid) * ( S(:, 1:n)' * Z );
+    return;
+end
+if m > 2*n
+    ATA = Z'*Z;
+    [S, Sigma2, D] = svd(ATA);
+    Sigma2 = diag(Sigma2);
+    V = sqrt(Sigma2);
+    tol = max(size(Z)) * eps(max(V));
+    n = sum(V > max(tol, tau));
+    
+    if(isinf(TraceNormToRankParam))
+        sub = zeros(1,n);
+    else
+        sub = ([0 [1:n-1]/(n-1)].^TraceNormToRankParam)*tau;
+    end
+    
+    mid = max(V(1:n)-sub', 0) ./ V(1:n) ;
+    X = (Z *  D(:, 1:n)) * bsxfun(@times, D(:, 1:n)', mid);
+    %X = Z * ( D(:, 1:n) * diag(mid) * D(:, 1:n)' );
+    return;
+    
+    %[X, n, Sigma2] = Pro2TraceNorm(Z', tau, TraceNormToRankParam);
+    %X = X';
+    %return;
+end
+
+%%{
+[S,V,D] = svd(Z);
+Sigma2 = diag(V).^2;
+n = sum(diag(V) > tau);
+
+if(isinf(TraceNormToRankParam))
+    sub = zeros(1,n);
+else
+    sub=([0 [1:n-1]/(n-1)].^TraceNormToRankParam)*tau;
+end
+V(1:n,1:n) = diag(diag(V(1:n,1:n))-sub');
+%V(2:n,2:n) = max(V(2:n,2:n)-subtractval, 0);
+X = S(:, 1:n) * V(1:n,1:n) * D(:, 1:n)';
+%}
+

private/Unfold.m

+function [X] = Unfold( X, dim, i )
+X = reshape(shiftdim(X,i-1), dim(i), []);
\ No newline at end of file