text:
function [x]=dmrg_rake_solve2(A, y, tol, varargin)
%DMRG-type method for the solution of linear systems in QTT-Tucker format
%   [X]=DMRG_RAKE_SOLVE2(A,Y,TOL,VARARGIN) Attempts to solve the linear
%   system A*X = Y with accuracy EPS using the two-sided DMRG iteration.
%   Matrix A has to be given in the QTT-Tucker, right-hand side Y should be
%   given in the QTT-Tucker format also. Options are provided in form
%   'PropertyName1',PropertyValue1,'PropertyName2',PropertyValue2 and so
%   on. The parameters are set to default (in brackets in the following)
%   The list of option names and default values are:
%       o x0 - initial approximation [random rank-2 tensor]
%       o nswp - maximal number of DMRG sweeps [10]
%       o rmax - maximal TT-rank of the solution [1000]
%       o verb - verbosity level, 0-silent, 1-sweep info, 2-block info [1]
%       o kick_rank - stabilization parameter [2]
%       o max_full_size - maximal size of the local matrix to full solver
%       [2500]
%       o local_prec: Local preconditioner, 'als' - ALS-Richardson
%       iteration, 'selfprec' (Saad selfpreconditioner) ['als']
%       o gmres_iters - number of local gmres restarts [2]
%       o nrestart - dimension of local gmres [25]
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

nswp = 20;
local_format = 'full';
% local_format = 'tt';
max_full_size = 2500;
max_full_size2 = Inf;
nrestart = 40;
gmres_iters = 2;
verb = 1;
kickrank = 2;
% checkrank = 1;
resid_damp_loc = 2;
rmax = Inf;
trunc_norm = 'matrix';

tol2 = tol;

x = [];

for i=1:2:length(varargin)-1
    switch lower(varargin{i})
        case 'nswp'
            nswp=varargin{i+1};
        case 'rmax'
           rmax=lower(varargin{i+1});
        case 'x0'
            x=varargin{i+1};
        case 'verb'
            verb=varargin{i+1};
%         case 'local_prec'
%             local_prec=varargin{i+1};
        case 'nrestart'
            nrestart=varargin{i+1};
        case 'gmres_iters'
            gmres_iters=varargin{i+1};
        case 'kickrank'
            kickrank=varargin{i+1};
        case  'max_full_size'
            max_full_size=varargin{i+1};
        case  'resid_damp'
            resid_damp_loc=varargin{i+1};
        case  'trunc_norm'
            trunc_norm=varargin{i+1};
%         case 'prec_compr'
%             prec_compr=varargin{i+1};
%         case 'prec_tol'
%             prec_tol=varargin{i+1};
%         case 'prec_iters'
%             prec_iters=varargin{i+1};
%         case 'use_self_prec'
%             use_self_prec=varargin{i+1};
%         case 'ddpow'
%             ddpow=varargin{i+1};
%         case 'ddrank'
%             ddrank=varargin{i+1};
%         case 'd_pow_check'
%             d_pow_check=varargin{i+1};
%         case 'bot_conv'
%             bot_conv=varargin{i+1};
%         case 'top_conv'
%             top_conv=varargin{i+1};
%         case 'min_dpow'
%             min_dpow=varargin{i+1};
%         case 'min_drank'
%             min_drank=varargin{i+1};
        otherwise
            error('Unrecognized option: %s\n',varargin{i});
    end
end

d = y.dphys; % Physical dim.
yc = y.core;
yf = y.tuck;
Af = A.tuck;
Ac = A.core;

L = zeros(1,d); % Quantics dims
n = zeros(max(L), d); % Physical mode sizes
for i=1:d
    L(i) = yf{i}.d;
    n(1:L(i), i) = yf{i}.n;
end;

if (isempty(x))
    xc = tt_rand(2,d,2);
    xf = cell(d,1);
    for i=1:d
        xf{i} = tt_rand(n(1:L(i),i), L(i), [1;2*ones(L(i),1)]);
    end;
else
    xc = x.core;
    xf = x.tuck;
end;

% Extract ranks; Note that rf(L(i)+1,i) = r_tuck(i)
rcy = yc.r;
rfy = zeros(max(L)+1, d);
for i=1:d
    rfy(1:L(i)+1, i) = yf{i}.r;
end;
rcA = Ac.r;
rfA = zeros(max(L)+1, d);
for i=1:d
    rfA(1:L(i)+1, i) = Af{i}.r;
end;
rcx = xc.r;
rfx = zeros(max(L)+1, d);
for i=1:d
    rfx(1:L(i)+1, i) = xf{i}.r;
end;

% Init phis. Thousands of them... (c)
phcA = cell(d+1,1); phcA{1} = 1; phcA{d+1}=1; % Core, matrix
phcy = cell(d+1,1); phcy{1} = 1; phcy{d+1}=1; % Core, rhs
phfA = cell(d,1); % factors, matrix
phfy = cell(d,1); % factors, rhs
phAfc = cell(d,1); % factor-core, matrix
phyfc = cell(d,1); % factor-core, rhs
for i=1:d
    phfA{i} = cell(L(i)+1,1);
    phfA{i}{1} = 1; phfA{i}{L(i)+1} = 1;
    phfy{i} = cell(L(i)+1,1);
    phfy{i}{1} = 1; phfy{i}{L(i)+1} = 1;
end;

% For random check
% cphcA = cell(d+1,1); cphcA{1} = 1; cphcA{d+1}=1; % Core, matrix
% cphcy = cell(d+1,1); cphcy{1} = 1; cphcy{d+1}=1; % Core, rhs
% cphfA = cell(d,1); % factors, matrix
% cphfy = cell(d,1); % factors, rhs
% cphAfc = cell(d,1); % factor-core, matrix
% cphyfc = cell(d,1); % factor-core, rhs
% for i=1:d
%     cphfA{i} = cell(L(i)+1,1);
%     cphfA{i}{1} = 1; cphfA{i}{L(i)+1} = 1;
%     cphfy{i} = cell(L(i)+1,1);
%     cphfy{i}{1} = 1; cphfy{i}{L(i)+1} = 1;
% end;

last_sweep = false;

for swp=1:nswp
    % init check vector
%     rcchk = [1; checkrank*ones(d-1,1); 1];
%     rfchk = zeros(max(L)+1, d);
%     for i=1:d
%         rfchk(1:L(i)+1, i) = [1; checkrank*ones(L(i),1)];
%     end;

    dx_max = 0;
    res_max = 0;
    r_max = 0;
%     chk_res_max = 0;
    % bottom-to-top QR and phis
    for i=d:-1:1 % physical dims/core
        for j=1:L(i) % quantics dims
            cr = xf{i}{j};
            cr = reshape(cr, rfx(j,i)*n(j,i), rfx(j+1,i));
            [cr, rv] = qr(cr, 0);
            % What is our next core?
            if (j<L(i))
                % we are still on a "tooth"
                cr2 = xf{i}{j+1};
                cr2 = reshape(cr2, rfx(j+1,i), n(j+1,i)*rfx(j+2,i));
                cr2 = rv*cr2;
                rfx(j+1,i) = size(cr, 2);
                xf{i}{j} = reshape(cr, rfx(j,i), n(j,i), rfx(j+1,i));
                xf{i}{j+1} = reshape(cr2, rfx(j+1,i), n(j+1,i), rfx(j+2,i));
            else
                % We have to convlove rv to the tucker core
                cr2 = xc{i};
                cr2 = permute(cr2, [2, 1, 3]);
                cr2 = reshape(cr2, rfx(j+1,i), rcx(i)*rcx(i+1));
                cr2 = rv*cr2;
                rfx(j+1,i) = size(cr, 2);
                cr2 = reshape(cr2, rfx(j+1,i), rcx(i), rcx(i+1));
                cr2 = permute(cr2, [2, 1, 3]);
                xf{i}{j} = reshape(cr, rfx(j,i), n(j,i), rfx(j+1,i));
                xc{i} = cr2;
            end;
            % Update bottom phis
            cr = reshape(cr, rfx(j,i), n(j,i), rfx(j+1,i));
            if (j<L(i))
                phfA{i}{j+1} = compute_next_Phi(phfA{i}{j}, cr, Af{i}{j}, cr, 'lr');
                phfy{i}{j+1} = compute_next_Phi(phfy{i}{j}, cr, [], yf{i}{j}, 'lr');
            else
                phAfc{i} = compute_next_Phi(phfA{i}{j}, cr, Af{i}{j}, cr, 'lr');
                phyfc{i} = compute_next_Phi(phfy{i}{j}, cr, [], yf{i}{j}, 'lr');
            end;

            % check vector
%             ccr = ones(1, n(j,i), 1);
%             if (j<L(i))
%                 cphfA{i}{j+1} = compute_next_Phi(cphfA{i}{j}, ccr, Af{i}{j}, cr, 'lr');
%                 cphfy{i}{j+1} = compute_next_Phi(cphfy{i}{j}, ccr, [], yf{i}{j}, 'lr');
%             else
%                 cphAfc{i} = compute_next_Phi(cphfA{i}{j}, ccr, Af{i}{j}, cr, 'lr');
%                 cphyfc{i} = compute_next_Phi(cphfy{i}{j}, ccr, [], yf{i}{j}, 'lr');
%             end;
        end;
    end;

    % QRs and phis over the core
    % Project the system on the factors
    Acr = tt_matrix(Ac, Ac.n, ones(d,1));
%     cAcr = tt_matrix(Ac, Ac.n, ones(d,1));
    ycr = yc;
%     cycr = yc;
    for i=1:d
        Acr{i} = core_matrix(Ac{i}, phAfc{i});
        ycr{i} = core_vector(yc{i}, phyfc{i});
%         cAcr{i} = core_matrix(Ac{i}, cphAfc{i});
%         cycr{i} = core_vector(yc{i}, cphyfc{i});
    end;
    for i=d:-1:2
        rtx = rfx(L(i)+1, i);
        cr = xc{i};
        cr = reshape(cr, rcx(i), rtx*rcx(i+1));
        [cr, rv] = qr(cr.', 0);
        cr2 = xc{i-1};
        rtx2 = rfx(L(i-1)+1, i-1);
        cr2 =  reshape(cr2, rcx(i-1)*rtx2, rcx(i));
        cr2 = cr2*(rv.');
        rcx(i) = size(cr, 2);
        cr = reshape(cr.', rcx(i), rtx, rcx(i+1));
        xc{i-1} = reshape(cr2, rcx(i-1), rtx2, rcx(i));
        xc{i} = cr;
        % Update right phi
        phcA{i} = compute_next_Phi(phcA{i+1}, cr, Acr{i}, cr, 'rl');
        phcy{i} = compute_next_Phi(phcy{i+1}, cr, [], ycr{i}, 'rl');

        % check vector
%         ccr = 1;
%         % Update right phi
%         cphcA{i} = compute_next_Phi(cphcA{i+1}, ccr, cAcr{i}, cr, 'rl');
%         cphcy{i} = compute_next_Phi(cphcy{i+1}, ccr, [], cycr{i}, 'rl');
    end;

%     % DMRG over the core
%     % Compute the reduced matrix and rhs;
%     rta = zeros(d,1);
%     for i=1:d
%         rta(i) = rfA(L(i)+1,i);
%     end;
%     Acr = tt_matrix(Ac, rta, ones(d,1));
%     ycr = yc;
%     for i=1:d
%         a1 = Ac{i};
%         a1 = permute(a1, [2, 1, 3]);
%         a1 = reshape(a1, rfA(L(i)+1,i), rcA(i)*rcA(i+1));
%         ph2 = phfAb{i};
%         ph2 = permute(ph2, [1,3,2]);
%         ph2 = reshape(ph2, rfx(L(i)+1,i)*rfx(L(i)+1,i), rfA(L(i)+1,i));
%         a1 = ph2*a1;
%         a1 = reshape(a1, rfx(L(i)+1,i), rfx(L(i)+1,i), rcA(i), rcA(i+1));
%         a1 = permute(a1, [3, 1, 2, 4]);
%         Acr{i} = a1;
%
%         y1 = yc{i};
%         y1 = permute(y1, [2, 1, 3]);
%         y1 = reshape(y1, rfy(L(i)+1,i), rcy(i)*rcy(i+1));
%         ph2 = phfyb{i};
%         y1 = ph2*y1;
%         y1 = reshape(y1, rfx(L(i)+1,i), rcy(i), rcy(i+1));
%         y1 = permute(y1, [2, 1, 3]);
%         ycr{i} = y1;
%     end;
%     fprintf('=rake_solve========= core processing, sweep %d =======\n', swp);
%     xcold = xc;
%     xc = dmrg_solve2(Acr, ycr, tol, 'x0', xcold, 'max_full_size', max_full_size, 'nrestart', nrestart, 'gmres_iters', gmres_iters, 'nswp', 1);
%     res = norm(Acr*xcold-ycr)/norm(ycr);
%     dx = norm(xcold-xc)/norm(xc);
%     dx_max = max(dx_max, dx);
%     r_max = max([r_max; rank(xc)]);
%     fprintf('=rake_solve========= core res_prev: %3.3e, dx: %3.3e, rmax: %d\n', res, dx, max(rank(xc)));
%     res_max = max(res_max, res);

%     % Horisontal phis
%     rcx = xc.r;
%     for i=d:-1:2
%         % n here is in fact a tucker rank
%         cr2 = xc{i};
%         cr2 = reshape(cr2, rcx(i), rfx(L(i)+1,i)*rcx(i+1));
%         [cr2, rv]=qr(cr2.', 0);
%         cr3 = xc{i-1};
%         cr3 = reshape(cr3, rcx(i-1)*rfx(L(i-1)+1,i-1), rcx(i));
%         cr3 = cr3*(rv.');
%         rcx(i) = size(cr2, 2);
%         cr2 = cr2.'; % sizes rcx1, rtuck*rcx2
%         xc{i} = reshape(cr2, rcx(i), rfx(L(i)+1,i), rcx(i+1));
%         xc{i-1} = reshape(cr3, rcx(i-1), rfx(L(i-1)+1,i-1), rcx(i));
%
%         % Update right phis
%         cr2 = reshape(cr2, rcx(i), rfx(L(i)+1,i), rcx(i+1));
%         phcAr{i} = compute_next_Phi(phcAr{i+1}, cr2, Acr{i}, cr2, 'rl');
%         % New size: rcx1, rcA1, rcx1
%         phcyr{i} = compute_next_Phi(phcyr{i+1}, cr2, [], ycr{i}, 'rl');
%         % New size: rcx1, rcy1
%     end;

    % DMRG over the factors
    for i=1:d
        % Convolve the tucker block to the last physical
        % Preparing the matrix and rhs in this case smells like
        % manure.
        % phcAl - Ac{i} - phcAr     <- this has to be convolved
        %          |                <- this has to be convolved
        %         Af{i}{j}
        %          |
        %         Af{i}{j-1}
        %          |                <- this has to be convolved
        %         phfAb{j-1}
        rta = rfA(L(i)+1,i); rtx = rfx(L(i)+1,i); rty = rfy(L(i)+1,i); % rtchk = rfchk(L(i)+1,i);
        a1left = phcA{i};
        a1left = permute(a1left, [1,3,2]);
        a1left = reshape(a1left, rcx(i)*rcx(i), rcA(i));
        a1left = a1left*reshape(Ac{i}, rcA(i), rta*rcA(i+1));

%         % check
%         ca1left = cphcA{i};
%         ca1left = permute(ca1left, [1,3,2]);
%         ca1left = reshape(ca1left, rcchk(i)*rcx(i), rcA(i));
%         ca1left = ca1left*reshape(Ac{i}, rcA(i), rta*rcA(i+1));

        % old
%         % a1left is to use later to compute phcAl
%         ph2 = phcAr{i+1};
%         ph2 = permute(ph2, [2, 1, 3]);
%         ph2 = reshape(ph2, rcA(i+1), rcx(i+1)*rcx(i+1));
%         a1top = reshape(a1left, rcx(i)*rcx(i)*rta, rcA(i+1))*ph2;
%         a1top = reshape(a1top, rcx(i), rcx(i), rta, rcx(i+1), rcx(i+1));
%         a1top = permute(a1top, [1, 4, 2, 5, 3]);
%         % we need a1 as well, to check the  residual in tucker block
%         % splitting
%         a1top = reshape(a1top, rcx(i)*rcx(i+1)*rcx(i)*rcx(i+1), rta);
        % new
        a1left = reshape(a1left, rcx(i)*rcx(i), rta, rcA(i+1));
        a1left = permute(a1left, [2, 1, 3]);
        a1left = reshape(a1left, rta, rcx(i)*rcx(i)*rcA(i+1));

%         ca1left = reshape(ca1left, rcchk(i)*rcx(i), rta, rcA(i+1));
%         ca1left = permute(ca1left, [2, 1, 3]);
%         ca1left = reshape(ca1left, rta, rcchk(i)*rcx(i)*rcA(i+1));

        a2 = Af{i}{L(i)};
        a2 = reshape(a2, rfA(L(i),i)*n(L(i),i)*n(L(i),i), rta);
        % old
%         a2 = a2*a1top.';
%         a2 = reshape(a2, rfA(L(i),i), n(L(i),i), n(L(i),i), rcx(i)*rcx(i+1), rcx(i)*rcx(i+1));
%         a2 = permute(a2, [1, 2, 4, 3, 5]);
%         a2 = reshape(a2, rfA(L(i),i), n(L(i),i)*rcx(i)*rcx(i+1), n(L(i),i)*rcx(i)*rcx(i+1), 1);
        % new
        a2 = a2*a1left;
        a2 = reshape(a2, rfA(L(i),i), n(L(i),i), n(L(i),i), rcx(i), rcx(i), rcA(i+1));
        a2 = permute(a2, [1, 2, 4, 3, 5, 6]);
        a2 = reshape(a2, rfA(L(i),i), n(L(i),i)*rcx(i), n(L(i),i)*rcx(i), rcA(i+1));

%         % check
%         ca2 = Af{i}{L(i)};
%         ca2 = reshape(ca2, rfA(L(i),i)*n(L(i),i)*n(L(i),i), rta);
%         ca2 = ca2*ca1left;
%         ca2 = reshape(ca2, rfA(L(i),i), n(L(i),i), n(L(i),i), rcchk(i), rcx(i), rcA(i+1));
%         ca2 = permute(ca2, [1, 2, 4, 3, 5, 6]);
%         ca2 = reshape(ca2, rfA(L(i),i), n(L(i),i)*rcchk(i), n(L(i),i)*rcx(i), rcA(i+1));

        Afr = Af{i};
        Afr{L(i)} = a2;

%         cAfr = Af{i};
%         cAfr{L(i)} = ca2;

        y1left = phcy{i};
        y1left = y1left*reshape(yc{i}, rcy(i), rty*rcy(i+1));
        % y1left is to use later to compute phcyl
        ph2 = phcy{i+1};
        ph2 = ph2.';
        y1top = reshape(y1left, rcx(i)*rty, rcy(i+1))*ph2;
        y1top = reshape(y1top, rcx(i), rty, rcx(i+1));
        y1top = permute(y1top, [1, 3, 2]);
        y1top = reshape(y1top, rcx(i)*rcx(i+1), rty);
        y2 = yf{i}{L(i)};
        y2 = reshape(y2, rfy(L(i),i)*n(L(i),i), rty);
        y2 = y2*y1top.';
        y2 = reshape(y2, rfy(L(i),i), n(L(i),i)*rcx(i)*rcx(i+1), 1);
        yfr = yf{i};
        yfr{L(i)} = y2;

%         cy1left = cphcy{i};
%         cy1left = cy1left*reshape(yc{i}, rcy(i), rty*rcy(i+1));
%         % y1left is to use later to compute phcyl
%         ph2 = cphcy{i+1};
%         ph2 = ph2.';
%         cy1top = reshape(cy1left, rcchk(i)*rty, rcy(i+1))*ph2;
%         cy1top = reshape(cy1top, rcchk(i), rty, rcchk(i+1));
%         cy1top = permute(cy1top, [1, 3, 2]);
%         cy1top = reshape(cy1top, rcchk(i)*rcchk(i+1), rty);
%         cy2 = yf{i}{L(i)};
%         cy2 = reshape(cy2, rfy(L(i),i)*n(L(i),i), rty);
%         cy2 = cy2*cy1top.';
%         cy2 = reshape(cy2, rfy(L(i),i), n(L(i),i)*rcchk(i)*rcchk(i+1), 1);
%         cyfr = yf{i};
%         cyfr{L(i)} = cy2;

        x1 = xc{i};
        x1 = permute(x1, [2, 1, 3]);
        x1 = reshape(x1, rtx, rcx(i)*rcx(i+1));
        x2 = xf{i}{L(i)};
        x2 = reshape(x2, rfx(L(i),i)*n(L(i),i), rtx);
        x2 = x2*x1;
        x2 = reshape(x2, rfx(L(i),i), n(L(i),i)*rcx(i)*rcx(i+1), 1);
        xfr = xf{i};
        xfr{L(i)} = x2;

%          chk1 = chkc{i};
%          chk1 = permute(chk1, [2, 1, 3]);
%          chk1 = reshape(chk1, rtchk, rcchk(i)*rcchk(i+1));
%          chk2 = chkf{i}{L(i)};
%          chk2 = reshape(chk2, rfchk(L(i),i)*n(L(i),i), rtchk);
%          chk2 = chk2*chk1;
%          chk2 = reshape(chk2, rfchk(L(i),i), n(L(i),i)*rcchk(i)*rcchk(i+1), 1);
%          chkfr = chkf{i};
%          chkfr{L(i)} = chk2;

        curn = xfr.n;
%         curm = xfr.n;
%         curm(L(i)) = n(L(i),i); % *rcchk(i)*rcchk(i+1);
        curra = Afr.r;
        curry = yfr.r;
        currx = xfr.r;
%         currchk = [rfchk(1:L(i),i); 1];
%          currchk = ones(L(i)+1,1);
%          currchk = chkfr.r;
        % DMRG over the factor from L(i) to 1
        for j=L(i):-1:2
            % new
            if (j<L(i))
                Phi2 = phfA{i}{j+1};
%                 cPhi2 = cphfA{i}{j+1};
            else
                Phi2 = phcA{i+1};
%                 cPhi2 = cphcA{i+1};
            end;
            a2 = Afr{j};
            a1 = Afr{j-1};
            Phi1 = phfA{i}{j-1};
%             % check
%             ca2 = cAfr{j};
%             ca1 = cAfr{j-1};
%             cPhi1 = cphfA{i}{j-1};
            % old
%             a2 = phfAt{i}{j+1};
%             a2 = permute(a2, [2, 1, 3]);
%             a2 = reshape(a2, curra(j+1), currx(j+1)*currx(j+1));
%             a2 = reshape(Afr{j}, curra(j)*curn(j)*curn(j), curra(j+1))*a2;
%             a2 = reshape(a2, curra(j), curn(j), curn(j), currx(j+1), currx(j+1));
%             a2 = permute(a2, [2, 4, 3, 5, 1]);
%             a2 = reshape(a2, curn(j)*currx(j+1), curn(j)*currx(j+1), curra(j));
%             a1 = phfAb{i}{j-1};
%             a1 = permute(a1, [1, 3, 2]);
%             a1 = reshape(a1, currx(j-1)*currx(j-1), curra(j-1));
%             a1 = a1*reshape(Afr{j-1}, curra(j-1), curn(j-1)*curn(j-1)*curra(j));
%             a1 = reshape(a1, currx(j-1), currx(j-1), curn(j-1), curn(j-1), curra(j));
%             a1 = permute(a1, [1, 3, 2, 4, 5]);
%             a1 = reshape(a1, currx(j-1)*curn(j-1), currx(j-1)*curn(j-1), curra(j));

            y2 = phfy{i}{j+1}.';
            y2 = reshape(yfr{j}, curry(j)*curn(j), curry(j+1))*y2;
            y2 = reshape(y2, curry(j), curn(j)*currx(j+1));
            y2 = y2.';

            y1 = phfy{i}{j-1};
            y1 = y1*reshape(yfr{j-1}, curry(j-1), curn(j-1)*curry(j));
            y1 = reshape(y1, currx(j-1)*curn(j-1), curry(j));

            x2 = reshape(xfr{j}, currx(j), curn(j)*currx(j+1));
            x2 = x2.';
            x1 = reshape(xfr{j-1}, currx(j-1)*curn(j-1), currx(j));

%             % check
%             cy2 = cphfy{i}{j+1}.';
%             cy2 = reshape(cyfr{j}, curry(j)*curm(j), curry(j+1))*cy2;
%             cy2 = reshape(cy2, curry(j), curm(j)*currchk(j+1));
%             cy1 = cphfy{i}{j-1};
%             cy1 = cy1*reshape(cyfr{j-1}, curry(j-1), curm(j-1)*curry(j));
%             cy1 = reshape(cy1, currchk(j-1)*curm(j-1), curry(j));
%             cy = cy1*cy2;

            % new
            if (j==L(i))
                currx(j+1)=rcx(i+1);
%                 currchk(j+1)=rcchk(i+1);
                curn(j) = n(L(i),i)*rcx(i);
%                 curm(j) = n(L(i),i)*rcchk(i);
            end;

            if (verb>1)
                fprintf('=rake_solve2= swp %d, factor {%d}{%d}, ', swp, i, j);
            end;
            local_format = 'full';
            if (currx(j-1)*curn(j-1)*curn(j)*currx(j+1)>max_full_size2)
                local_format = 'tt';
            end;
            % old
%             [u,s,v,r,dx_max,res_max]=local_solve(a1, a2, y1, y2, x1, x2, ...
%                 currx(j-1), curn(j-1), curn(j), currx(j+1), curra(j), ...
%                 tol/sqrt(L(i))/sqrt(d)/2, res_max, dx_max, ...
%                 local_format, max_full_size, nrestart, gmres_iters, verb);
            % new
            [u,s,v,r,dx_max,res_max]=local_solve(Phi1,a1, a2, Phi2, y1, y2, x1, x2, ...
                currx(j-1), curn(j-1), curn(j), currx(j+1), curra(j), ...
                tol2/sqrt(sum(L)), res_max, dx_max, resid_damp_loc, trunc_norm, ...
                local_format, max_full_size, nrestart, gmres_iters, verb);
            % old rounding tol: tol2/sqrt(L(i))/sqrt(d)/2
	    r = min(rmax, r);
	    u = u(:,1:r); s = s(1:r,1:r); v = v(:,1:r);
            u = u*s;

%             % check
%             Asol = bfun3(cPhi1,ca1,ca2,cPhi2, u*(v.'));
%             chk_res_max = max(chk_res_max, norm(Asol-cy(:))/norm(cy(:)));

            % kick
            if (~last_sweep)
%                 Axprev = bfun3(Phi1,a1, a2, Phi2, x1*x2.');
%                 Axprev = reshape(Axprev, currx(j-1)*curn(j-1), curn(j)*currx(j+1));
%                 [unew,snew,vnew]=svd(Axprev, 'econ');
%                 v = reort(v, vnew(:,1:min(kickrank, size(vnew,2))));
%                 v = reort(v, rand(curn(j)*currx(j+1), kickrank));
                [v,rv]=qr([v, rand(curn(j)*currx(j+1), kickrank)], 0);
                radd = kickrank;
                u = [u, zeros(currx(j-1)*curn(j-1), radd)];
                u = u*(rv.');
            end;
            r = size(v,2);
            xfr{j} = reshape(v.', r, curn(j), currx(j+1));
            xfr{j-1} = reshape(u, currx(j-1), curn(j-1), r);
            currx(j) = r;
            rfx(j,i)=r;
            r_max = max(r_max, r);
            % old
            % new phis
%             a2 = permute(a2, [1, 3, 2]);
%             a2 = reshape(a2, curn(j)*currx(j+1)*curra(j), curn(j)*currx(j+1));
%             a2 = a2*v;
%             a2 = reshape(a2, curn(j)*currx(j+1), curra(j)*r);
%             a2 = (v')*a2;
%             phfAt{i}{j} = reshape(a2, r, curra(j), r);
            phfy{i}{j} = (v')*y2;
            %new
            phfA{i}{j} = compute_next_Phi(Phi2, xfr{j}, a2, xfr{j}, 'rl');

%             % check vector
%             ccr = ones(n(j,i), 1);
%             cphfy{i}{j} = (ccr')*(cy2.');
%             cphfA{i}{j} = compute_next_Phi(cPhi2, ccr.', ca2, xfr{j}, 'rl');
        end;

%         fprintf('=rake_solve========= factor {%d} processing, Sweep %d =======\n', i, swp);
%         xfrold = xfr;
%         xfr = dmrg_solve2(Afr, yfr, tol, 'x0', xfrold, 'max_full_size', max_full_size, 'nrestart', nrestart, 'gmres_iters', gmres_iters, 'nswp', 1);
%         res = norm(Afr*xfrold-yfr)/norm(yfr);
%         dx = norm(xfrold-xfr)/norm(xfr);
%         dx_max = max(dx_max, dx);
%         r_max = max([r_max; rank(xfr)]);
%         fprintf('=rake_solve========= factor {%d} res_prev: %3.3e, dx: %3.3e, rmax: %d\n', i, res, dx, max(rank(xfr)));
%         res_max = max(res_max, res);
%         rfx(1:L(i)+1,i) = xfr.r;

        % We have to split the tucker block, and compute new phf*b
        for j=1:L(i)
            cr = xfr{j};
            % What is our next core?
            if (j<L(i))
                % we are still on a "tooth"
                cr = reshape(cr, rfx(j,i)*n(j,i), rfx(j+1,i));
                [cr, rv] = qr(cr, 0);
                cr2 = xfr{j+1};
                ncur = size(cr2, 2);
                r3cur = size(cr2, 3);
                cr2 = reshape(cr2, rfx(j+1,i), ncur*r3cur);
                cr2 = rv*cr2;
                rfx(j+1,i) = size(cr, 2);
                xfr{j} = reshape(cr, rfx(j,i), n(j,i), rfx(j+1,i));
                xfr{j+1} = reshape(cr2, rfx(j+1,i), ncur, r3cur);
            else
                % We have to split the tucker core and update the tucker
                % rank
                cr = reshape(cr, rfx(j,i)*n(j,i), rcx(i)*rcx(i+1));
                [u,s,v]=svd(cr, 'econ');
                if (strcmp(trunc_norm, 'fro'))
                    r = my_chop2(diag(s), tol/sqrt(L(i))/sqrt(d)/2*norm(diag(s)));
                else
                % Prepare the local matrix and rhs for residue check
                % new
                Phi2 = phcA{i+1};
                curA2 = reshape(a1left, rta, rcx(i), rcx(i), rcA(i+1));
                curA1 = Af{i}{L(i)};
                Phi1 = phfA{i}{j};
                % old
%                 curA = cell(2,1);
%                 curA{2} = reshape(a1top, rcx(i)*rcx(i+1), rcx(i)*rcx(i+1), rta);
%                 curA{1} = Af{i}{L(i)};
%                 curA{1} = reshape(curA{1}, rfA(j,i), n(j,i)*n(j,i)*rta);
%                 ph2 = phfAb{i}{j};
%                 ph2 = permute(ph2, [1,3,2]);
%                 ph2 = reshape(ph2, rfx(j,i)*rfx(j,i), rfA(j,i));
%                 curA{1} = ph2*curA{1};
%                 curA{1} = reshape(curA{1}, rfx(j,i), rfx(j,i), n(j,i), n(j,i), rta);
%                 curA{1} = permute(curA{1}, [1, 3, 2, 4, 5]);
%                 curA{1} = reshape(curA{1}, rfx(j,i)*n(j,i), rfx(j,i)*n(j,i), rta);

                rhs = reshape(yf{i}{j}, rfy(j,i), n(j,i)*rty);
                rhs = phfy{i}{j}*rhs;
                rhs = reshape(rhs, rfx(j,i)*n(j,i), rty);
                rhs = rhs*(y1top.');
                rhs = reshape(rhs, rfx(j,i)*n(j,i)*rcx(i)*rcx(i+1),1);
                r = 1;
                normy = norm(rhs);
                % old
%                 res_true = norm(bfun2(curA, cr, rfx(j,i), n(j,i), rcx(i), rcx(i+1), rfx(j,i), n(j,i), rcx(i), rcx(i+1))-rhs)/normy;
                % new
                res_true = norm(bfun3(Phi1, curA1, curA2, Phi2, cr)-rhs)/normy;
                while (r<=size(s,1))
                    cursol = u(:,1:r)*s(1:r,1:r)*(v(:,1:r)');
                    % old
%                     res = norm(bfun2(curA, cursol, rfx(j,i), n(j,i), rcx(i), rcx(i+1), rfx(j,i), n(j,i), rcx(i), rcx(i+1))-rhs)/normy;
                    % new
                    res = norm(bfun3(Phi1, curA1, curA2, Phi2, cursol)-rhs)/normy;
                    if (res<max(tol/sqrt(sum(L)), res_true*resid_damp_loc))
                        % old: tol/sqrt(L(i))/sqrt(d)/2
                        break;
                    end;
                    r = r+1;
                end;
                end;

                if (verb>1)
                    fprintf('=rake_solve2= swp %d, tuckerrank {%d}, res: %3.3e, r: %d\n', swp, i, res, r);
                end;
                r = min(rmax, r);
                u = u(:,1:r);
                v = conj(v(:,1:r));
                s = s(1:r,1:r);
                v = v*s;
                if (~last_sweep)
%                     Axprev = bfun3(Phi1, curA1, curA2, Phi2, cr);
%                     Axprev = reshape(Axprev, rfx(j,i)*n(j,i), rcx(i)*rcx(i+1));
%                     [unew,snew,vnew]=svd(Axprev, 'econ');
%                     u = reort(u, unew(:,1:min(kickrank, size(unew,2))));
%                     u = reort(u, rand(rfx(j,i)*n(j,i), kickrank));
                    [u,rv]=qr([u, rand(rfx(j,i)*n(j,i), kickrank)], 0);
                    radd = kickrank;
                    v = [v, zeros(rcx(i)*rcx(i+1), radd)];
                    v = v*(rv.');
                end;
                r = size(u,2);
                rfx(j+1,i) = r;
                cr = u;
                xfr{j} = reshape(cr, rfx(j,i), n(j,i), r);
                v = reshape(v.', r, rcx(i), rcx(i+1));
                v = permute(v, [2, 1, 3]);
                xc{i} = v;
                r_max = max(r_max, r);
            end;
            % Update bottom phis
            cr = reshape(cr, rfx(j,i), n(j,i), rfx(j+1,i));
            if (j<L(i))
                phfA{i}{j+1} = compute_next_Phi(phfA{i}{j}, cr, Af{i}{j}, cr, 'lr');
                phfy{i}{j+1} = compute_next_Phi(phfy{i}{j}, cr, [], yf{i}{j}, 'lr');
            else
                phAfc{i} = compute_next_Phi(phfA{i}{j}, cr, Af{i}{j}, cr, 'lr');
                phyfc{i} = compute_next_Phi(phfy{i}{j}, cr, [], yf{i}{j}, 'lr');
            end;

%             % check vector
%             ccr = ones(1, n(j,i), 1);
%             if (j<L(i))
%                 cphfA{i}{j+1} = compute_next_Phi(cphfA{i}{j}, ccr, Af{i}{j}, cr, 'lr');
%                 cphfy{i}{j+1} = compute_next_Phi(cphfy{i}{j}, ccr, [], yf{i}{j}, 'lr');
%             else
%                 cphAfc{i} = compute_next_Phi(cphfA{i}{j}, ccr, Af{i}{j}, cr, 'lr');
%                 cphyfc{i} = compute_next_Phi(cphfy{i}{j}, ccr, [], yf{i}{j}, 'lr');
%             end;
        end;
        xf{i} = xfr;

        % Now, perform the DMRG over the core to the next factor, update phc*l as well
        if (i<d)
            rtx = rfx(L(i)+1,i); rtx2 = rfx(L(i+1)+1,i+1);
            rx1 = rcx(i); rx2 = rcx(i+1); rx3 = rcx(i+2);
            ra2 = rcA(i+1); ra3 = rcA(i+2);
            ry2 = rcy(i+1); ry3 = rcy(i+2);

%             rtchk = rfchk(L(i)+1,i); rtchk2 = rfchk(L(i+1)+1,i+1);
%             rchk1 = rcchk(i); rchk3 = rcchk(i+2);

            % Acr and ycr are okey, except the i-th core
            % new
            Phi1 = phcA{i};
            a1 = core_matrix(Ac{i}, phAfc{i});
            a2 = Acr{i+1};
            Phi2 = phcA{i+2};

%             cPhi1 = cphcA{i};
%             ca1 = core_matrix(Ac{i}, cphAfc{i});
%             ca2 = cAcr{i+1};
%             cPhi2 = cphcA{i+2};
            % old
%             a1 = reshape(a1left, rx1*rx1, rta, ra2);
%             a1 = core_matrix(a1, phAfc{i});
%             a1 = reshape(a1, rx1, rx1, rtx, rtx, ra2);
%             a1 = permute(a1, [1, 3, 2, 4, 5]);
%             a1 = reshape(a1, rx1*rtx, rx1*rtx, ra2);
%             a2 = phcAr{i+2};
%             a2 = permute(a2, [2, 1, 3]);
%             a2 = reshape(a2, ra3, rx3*rx3);
%             a2 = reshape(Acr{i+1}, ra2*rtx2*rtx2, ra3)*a2;
%             a2 = reshape(a2, ra2, rtx2, rtx2, rx3, rx3);
%             a2 = permute(a2, [2, 4, 3, 5, 1]);
%             a2 = reshape(a2, rtx2*rx3, rtx2*rx3, ra2);

            y1 = reshape(y1left, rx1, rty, ry2);
            y1 = core_vector(y1, phyfc{i});
            y1 = reshape(y1, rx1*rtx, ry2);

            y2 = phcy{i+2}.';
            y2 = reshape(ycr{i+1}, ry2*rtx2, ry3)*y2;
            y2 = reshape(y2, ry2, rtx2*rx3);
            y2 = y2.';

            x1 = reshape(xc{i}, rx1*rtx, rx2);
            x2 = reshape(xc{i+1}, rx2, rtx2*rx3);
            x2 = x2.';

            if (verb>1)
                fprintf('=rake_solve2= swp %d, core {%d}, ', swp, i);
            end;
            local_format = 'full';
            if (rx1*rtx*rtx2*rx3>max_full_size2)
                local_format = 'tt';
            end;
            % new
            [u,s,v,r,dx_max,res_max]=local_solve(Phi1, a1, a2, Phi2, y1, y2, x1, x2, ...
                rx1, rtx, rtx2, rx3, ra2, ...
                tol2/sqrt(sum(L)), res_max, dx_max, resid_damp_loc, trunc_norm, ...
                local_format, max_full_size, nrestart, gmres_iters, verb);
            % old tol: tol2/sqrt(d)/2
            % old
%             [u,s,v,r,dx_max,res_max]=local_solve(a1, a2, y1, y2, x1, x2, ...
%                 rx1, rtx, rtx2, rx3, ra2, ...
%                 tol/sqrt(d)/2, res_max, dx_max, ...
%                 local_format, max_full_size, nrestart, gmres_iters, verb);
            r = min(rmax, r);
            u = u(:,1:r); s = s(1:r,1:r); v = v(:,1:r);
            v = v*s;

%             % check
%             Asol = bfun3(cPhi1, ca1, ca2, cPhi2, u*(v.'));
%             cy1 = reshape(cy1left, rchk1, rty, ry2);
%             cy1 = core_vector(cy1, cphyfc{i});
%             cy1 = reshape(cy1, rchk1*rtchk, ry2);
%             cy2 = cphcy{i+2}.';
%             cy2 = reshape(cycr{i+1}, ry2*rtchk2, ry3)*cy2;
%             cy2 = reshape(cy2, ry2, rtchk2*rchk3);
%             cy = cy1*cy2;
%             chk_res_max = max(chk_res_max, norm(Asol-cy(:))/norm(cy(:)));

            % kick
            if (~last_sweep)
%                 Axprev = bfun3(Phi1, a1, a2, Phi2, x1*x2.');
%                 Axprev = reshape(Axprev, rx1*rtx, rtx2*rx3);
%                 [unew,snew,vnew]=svd(Axprev, 'econ');
%                 u = reort(u, unew(:,1:min(kickrank, size(unew,2))));
%                 u = reort(u, rand(rx1*rtx, kickrank));
                [u,rv]=qr([u, rand(rx1*rtx, kickrank)], 0);
                radd = kickrank;
                v = [v, zeros(rtx2*rx3, radd)];
                v = v*(rv.');
            end;
            r = size(u,2);
            xc{i} = reshape(u, rx1, rtx, r);
            xc{i+1} = reshape(v.', r, rtx2, rx3);
            rcx(i+1)=r;
            r_max = max(r_max, r);
            % new phis
            % old
%             a1 = permute(a1, [1, 3, 2]);
%             a1 = reshape(a1, rx1*rtx*ra2, rx1*rtx);
%             a1 = a1*u;
%             a1 = reshape(a1, rx1*rtx, ra2*r);
%             a1 = (u')*a1;
%             phcAl{i+1} = reshape(a1, r, ra2, r);
            phcy{i+1} = (u')*y1;
            % new
            phcA{i+1} = compute_next_Phi(phcA{i}, xc{i}, a1, xc{i}, 'lr');

%             ccr = 1;
%             cphcy{i+1} = (ccr')*cy1;
%             ccr = reshape(ccr, rchk1, rtchk, rcchk(i+1));
%             cphcA{i+1} = compute_next_Phi(cphcA{i}, ccr, ca1, xc{i}, 'lr');

%             cr = xc{i};
%             cr = reshape(cr, rcx(i)*rtx, rcx(i+1));
%             [cr, rv]=qr(cr, 0);
%             cr2 = xc{i+1};
%             cr2 = reshape(cr2, rcx(i+1), rfx(L(i+1)+1,i+1)*rcx(i+2));
%             cr2 = rv*cr2;
%             rcx(i+1) = size(cr, 2);
%             xc{i+1} = reshape(cr2, rcx(i+1), rfx(L(i+1)+1,i+1), rcx(i+2));
%             xc{i} = reshape(cr, rcx(i), rtx, rcx(i+1));
%
%             % We have a1left, y1left of sizes rcx(i)(^2), rtuck, rc*(i+1)
%             curph = reshape(a1left, rcx(i)*rcx(i), rta, rcA(i+1));
%             curph = permute(curph, [2, 1, 3]);
%             curph = reshape(curph, rta, rcx(i)*rcx(i)*rcA(i+1));
%             ph2 = phfAb{i};
%             ph2 = permute(ph2, [1, 3, 2]);
%             ph2 = reshape(ph2, rtx*rtx, rta);
%             curph = ph2*curph;
%             curph = reshape(curph, rtx, rtx, rcx(i), rcx(i), rcA(i+1));
%             curph = permute(curph, [3, 1, 5, 4, 2]);
%             curph = reshape(curph, rcx(i)*rtx*rcA(i+1), rcx(i)*rtx);
%             curph = curph*cr;
%             curph = reshape(curph, rcx(i)*rtx, rcA(i+1)*rcx(i+1));
%             curph = (cr')*curph;
%             phcAl{i+1} = reshape(curph, rcx(i+1), rcA(i+1), rcx(i+1));
%
%             curph = reshape(y1left, rcx(i), rty, rcy(i+1));
%             curph = permute(curph, [2, 1, 3]);
%             curph = reshape(curph, rty, rcx(i)*rcy(i+1));
%             ph2 = phfyb{i};
%             curph = ph2*curph;
%             curph = reshape(curph, rtx, rcx(i), rcy(i+1));
%             curph = permute(curph, [2, 1, 3]);
%             curph = reshape(curph, rcx(i)*rtx, rcy(i+1));
%             curph = (cr')*curph;
%             phcyl{i+1} = curph;
        end;
    end;

    if (verb>0)
        real_res = NaN;
%  	x = qtt_tucker;
%  	x.dphys = d;
%  	x.tuck = xf;
%  	x.core = xc;
%
%  	real_res = mvrk(A, x, tol, 'verb', 0);
%  %          real_res = A*x;
%  	real_res = norm(real_res-y)/norm(y);
        fprintf('\n=rake_solve2= swp %d, dx_max: %3.3e, res_max: %3.3e, r_max: %d, real_res: %3.3e\n\n', swp, dx_max, res_max, r_max, real_res);
    end;
    if (last_sweep)
        break;
    end;
    if (strcmp(trunc_norm, 'fro'))
        if (dx_max<tol)||(swp==nswp-1)
            last_sweep = true;
        end;
    else
        if (res_max<tol)||(swp==nswp-1)
            last_sweep=true;
        end;
%         if (chk_res_max<tol)||(swp==nswp-1)
%             last_sweep = true;
%         end;
    end;
end;

x = qtt_tucker;
x.dphys = d;
x.tuck = xf;
x.core = xc;
end

function [Phi] = compute_next_Phi(Phi_prev, x, A, y, direction)
% Performs the recurrent Phi (or Psi) matrix computation
% Phi = Phi_prev * (x'Ay).
% If direction is 'lr', computes Psi
% if direction is 'rl', computes Phi
% A can be empty, then only x'y is computed.

if (strcmp(direction, 'rl'))
  % Revert ranks to perform the right-to-left recursion
  x = permute(x, [3, 2, 1]);
  y = permute(y, [3, 2, 1]);
  if (~isempty(A))
    A = permute(A, [4, 2, 3, 1]);
  end
end

rx1 = size(x,1); n = size(x,2); rx2 = size(x,3);
ry1 = size(y,1); m = size(y,2); ry2 = size(y,3);
if (~isempty(A))
  ra1 = size(A,1); ra2 = size(A,4);
else
  ra1 = 1; ra2 = 1;
end

Phi = reshape(Phi_prev, [rx1*ra1, ry1]);
y = reshape(y, [ry1, m*ry2]);
Phi = Phi*y;	% complexity §\mcommentfont$\mathcal{O}(n  r_x r_A r_y^2)$§
Phi = reshape(Phi, [rx1, ra1, m, ry2]);
Phi = permute(Phi, [2, 3, 1, 4]);
if (~isempty(A))
  Phi = reshape(Phi, [ra1*m, rx1*ry2]);
  A = permute(A, [4, 2, 1, 3]);
  A = reshape(A, [ra2*n, ra1*m]);
  Phi = A*Phi;	% complexity §\mcommentfont$\mathcal{O}(n^2  r_x r_A^2 r_y)$§
  Phi = reshape(Phi, [ra2, n, rx1, ry2]);
end
Phi = permute(Phi, [3, 2, 1, 4]);
Phi = reshape(Phi, [rx1*n, ra2*ry2]);
x = reshape(x, [rx1*n, rx2]);
Phi = (x')*Phi;	% complexity §\mcommentfont$\mathcal{O}(n  r_x^2 r_A r_y)$§
if (~isempty(A))
  Phi = reshape(Phi, [rx2, ra2, ry2]);
end
end

function [A] = core_matrix(core_block, Phi_factor)
% Computes the matrix to use in DMRG by the core, by convolving the core of
% the Tucker block and the last factor Phi matrix over the Tucker rank
% core_block is rc1 x rtuck x rc2, and
% Phi_factor is n1 x rtuck x n2
% A is then rc1 x n1 x n2 x rc2

r1 = size(core_block, 1); rtuck = size(core_block, 2); r2 = size(core_block, 3);
n1 = size(Phi_factor, 1); n2 = size(Phi_factor, 3);

Phi_factor = permute(Phi_factor, [1, 3, 2]);
Phi_factor = reshape(Phi_factor, n1*n2, rtuck);
A = permute(core_block, [2, 1, 3]);
A = reshape(A, rtuck, r1*r2);
A = Phi_factor*A;
A = reshape(A, n1, n2, r1, r2);
A = permute(A, [3, 1, 2, 4]);

end

function [A] = core_vector(core_block, Phi_factor)
% Computes the vector to use in DMRG by the core, by convolving the core of
% the Tucker block and the last factor Phi matrix over the Tucker rank
% core_block is rc1 x rtuck x rc2, and
% Phi_factor is n1 x rtuck
% A is then rc1 x n1 x rc2

r1 = size(core_block, 1); rtuck = size(core_block, 2); r2 = size(core_block, 3);
n1 = size(Phi_factor, 1);

A = permute(core_block, [2, 1, 3]);
A = reshape(A, rtuck, r1*r2);
A = Phi_factor*A;
A = reshape(A, n1, r1, r2);
A = permute(A, [2, 1, 3]);

end

function [y]=bfun2(B, x, rxm1, m1, m2, rxm3, rxn1, k1, k2, rxn3)
% Computes (B{1} \otimes B{2})x
% B{1} is of sizes rxn1*k1, rxm1*m1, rB
% B{2} is of sizes k2*rxn3, m2*rxm3, rB
rB=size(B{1},3);
x = reshape(x, rxm1*m1, m2*rxm3);
B1 = permute(B{1}, [3 1 2]);
B1 = reshape(B1, rB*rxn1*k1, rxm1*m1);
y = B1*x; % size rB*rxn1*k1,m2*rxm3  % cplx rB*rx^3*n^3
y = reshape(y, rB, rxn1*k1, m2*rxm3);
y = permute(y, [3 1 2]);
y = reshape(y, m2*rxm3*rB, rxn1*k1);
B2 = reshape(B{2}, k2*rxn3, m2*rxm3*rB);
y = B2*y; % size k2*rxn3,rxn1*k1 % cplx rB*rx^3*n^3
y = reshape(y.', rxn1*k1*k2*rxn3, 1);
end

function [y]=bfun3(Phi1,B1,B2,Phi2, x)
% Computes (Phi1 * B1 * B2 * Phi2)*x
% Phi1 is of sizes ry1, rB1, rx1
% B1 is of sizes rB1, k1, m1, rB2
% B2 is of sizes rB2, k2, m2, rB3
% Phi2 is of sizes ry3, rB3, rx3
ry1 = size(Phi1,1); ry3 = size(Phi2,1);
rx1 = size(Phi1,3); rx3 = size(Phi2,3);
rb1=size(B1,1); rb2=size(B1,4); rb3 = size(B2, 4);
m1 = size(B1,3); m2 = size(B2,3);
k1 = size(B1,2); k2 = size(B2,2);

y = reshape(x, rx1, m1*m2*rx3);
Phi1 = reshape(Phi1, ry1*rb1, rx1);
y = Phi1*y; % size ry1*rb1,m1*m2*rx3 % cplx rb*rx^3*m^2
y = reshape(y, ry1, rb1*m1, m2, rx3);
y = permute(y, [2, 1, 3, 4]);
y = reshape(y, rb1*m1, ry1*m2*rx3);
B1 = permute(B1, [2, 4, 1, 3]);
B1 = reshape(B1, k1*rb2, rb1*m1);
y = B1*y; % size k1*rb2, ry1*m2*rx3 % cplx rb^2*rx^2*n^3
y = reshape(y, k1, rb2, ry1, m2, rx3);
y = permute(y, [2, 4, 3, 1, 5]);
y = reshape(y, rb2*m2, ry1*k1*rx3);
B2 = permute(B2, [2, 4, 1, 3]);
B2 = reshape(B2, k2*rb3, rb2*m2);
y = B2*y; % size k2*rb3, ry1*k1*rx3 % cplx rb^2*rx^2*n^3
y = reshape(y, k2, rb3, ry1*k1, rx3);
y = permute(y, [2, 4, 3, 1]);
y = reshape(y, rb3*rx3, ry1*k1*k2);
Phi2 = reshape(Phi2, ry3, rb3*rx3);
y = Phi2*y; % size ry3, ry1*k1*k2 % cplx rb*rx^3*n^2
y = y.';
y = reshape(y, ry1*k1*k2*ry3, 1);
end

% old
% function [u,s,v,r,dx_max,res_max]=local_solve(a1, a2, y1, y2, x1, x2, ...
% new
function [u,s,v,r,dx_max,res_max]=local_solve(Phi1,a1, a2, Phi2, y1, y2, x1, x2, ...
    rx1, n1, n2, rx3, ra2, ...
    real_tol, res_max, dx_max,  resid_damp, trunc_norm, ...
    local_format, max_full_size, nrestart, gmres_iters, verb)

if (strcmp(local_format, 'full'))
    sol_prev = x1*(x2.');
    sol_prev = sol_prev(:);
    rhs = y1*(y2.');
    rhs = rhs(:);
    normy = norm(rhs);
    if (rx1*n1*n2*rx3<max_full_size)
        % new
        B = permute(Phi1, [1,3,2]);
        B = reshape(B, rx1*rx1, size(a1,1));
        a1 = reshape(a1, size(a1,1), n1*n1*ra2);
        B = B*a1;
        B = reshape(B, rx1, rx1, n1, n1, ra2);
        B = permute(B, [1, 3, 2, 4, 5]);
        B = reshape(B, rx1*n1*rx1*n1, ra2);
        ra3 = size(a2,4);
        a2 = reshape(a2, ra2, n2*n2*ra3);
        B = B*a2;
        B = reshape(B, rx1*n1*rx1*n1*n2*n2, ra3);
        Phi2 = permute(Phi2, [2, 1, 3]);
        Phi2 = reshape(Phi2, ra3, rx3*rx3);
        B = B*Phi2;
        B = reshape(B, rx1*n1, rx1*n1, n2, n2, rx3, rx3);
        B = permute(B, [1, 3, 5, 2, 4, 6]);
        % old
%         B = reshape(a1, rx1*n1*rx1*n1, ra2);
%         B = B*(reshape(a2, n2*rx3*n2*rx3, ra2).');
%         B = reshape(B, rx1*n1,  rx1*n1, n2*rx3, n2*rx3);
%         B = permute(B, [1, 3, 2, 4]);

        B = reshape(B, rx1*n1*n2*rx3, rx1*n1*n2*rx3);
%         nB = norm(B, 'fro');
%         B = B/nB;
%         rhs = rhs/nB;
%          sol = (B'*B+1e-15*eye(size(B))) \ (B'*rhs);
%
%          B = B*nB;
%          rhs = rhs*nB;
         sol = B \ rhs;

%          sol = (B'*B + 1e-16*norm(B'*B, 'fro')) \ (B'*rhs);
%         sol = pinv(B)*rhs;
        flg = 0;
%         [sol,flg]=gmres(B, rhs, ...
%             nrestart, real_tol/resid_damp, gmres_iters, [], [], sol);

        res_true = norm(B*sol-rhs)/normy;
        res_prev = norm(B*sol_prev-rhs)/normy;

%         if (res_true/res_prev>1)&&(res_true>real_tol)
% %             If the direct solution sucked
%             [sol,flg]=gmres(B, rhs, ...
%                 nrestart, real_tol/resid_damp, gmres_iters, [], [], sol_prev);
% %             sol = sol_prev;
%             res_true = norm(B*sol-rhs)/normy;
%         else
%             flg = 0;
%         end;

        if (flg>0)
            fprintf('--warn-- gmres did not converge\n');
        end;
    else
        B = cell(2,1);
        B{1} = a1;
        B{2} = a2;
        % old
%         res_prev = norm(bfun2(B, sol_prev, rx1, n1, n2, rx3, rx1, n1, n2, rx3)-rhs)/normy;
        % new
        drhs = bfun3(Phi1, a1, a2, Phi2, sol_prev)-rhs;
        res_prev = norm(drhs)/normy;

%         sol = als_solve_rx_2(B, rhs, real_tol, 10, sol_prev, [], 3);
%         [sol,flg]=bicgstab(@(v)bfun2(B, v, rx1, n1, n2, rx3, rx1, n1, n2, rx3), rhs, ...
%             max(real_tol,res_prev*0.1), nrestart*gmres_iters, [], [], sol_prev);
        % old
%         [sol,flg]=gmres(@(v)bfun2(B, v, rx1, n1, n2, rx3, rx1, n1, n2, rx3), rhs, ...
%             nrestart, max(real_tol,res_prev*0.05), gmres_iters, [], [], sol_prev);
        % new
        [dsol,flg]=gmres(@(v)bfun3(Phi1, a1, a2, Phi2, v), drhs, ...
            nrestart, min(real_tol/resid_damp/res_prev,1), gmres_iters);
        sol = sol_prev-dsol;
        if (flg>0)
            fprintf('--warn-- gmres did not converge\n');
        end;
        % old
%         res_true = norm(bfun2(B, sol, rx1, n1, n2, rx3, rx1, n1, n2, rx3)-rhs)/normy;
        % new
        res_true = norm(bfun3(Phi1, a1, a2, Phi2, sol)-rhs)/normy;
    end;

    if ((res_prev/res_true)<resid_damp)&&(res_true>real_tol/resid_damp)
        fprintf('--warn-- the residual damp by gmres was smaller than in the truncation\n');
    end;

    dx = norm(sol-sol_prev)/norm(sol);
    dx_max = max(dx_max, dx);
    if (rx1*n1*n2*rx3<max_full_size)
        Bx = B*sol; Bx_prev = B*sol_prev;
    else
        Bx = bfun3(Phi1, a1, a2, Phi2, sol);
        Bx_prev = bfun3(Phi1, a1, a2, Phi2, sol_prev);
    end;
    res_max = max(res_max, norm(Bx-Bx_prev)/norm(Bx));
%      res_max = max(res_max, res_prev);

    sol = reshape(sol, rx1*n1, n2*rx3);
    [u,s,v]=svd(sol, 'econ');
    s = diag(s);
    if (strcmp(trunc_norm, 'fro'))
        r = my_chop2(s, norm(s)*max(real_tol, res_true*resid_damp));
    else
        r1 = 1; r2 = numel(s); r = round((r1+r2)/2);
        while (r2-r1>1)
            cursol = u(:,1:r)*diag(s(1:r))*(v(:,1:r)');
            if (rx1*n1*n2*rx3<max_full_size)
                res = norm(B*cursol(:)-rhs)/normy;
            else
                % old
                %             res = norm(bfun2(B, cursol, rx1, n1, n2, rx3, rx1, n1, n2, rx3)-rhs)/normy;
                % new
                res = norm(bfun3(Phi1, a1, a2, Phi2, cursol)-rhs)/normy;
            end;
            if (res<max(real_tol, res_true*resid_damp))
                r2 = r;
            else
                r1 = r;
            end;
            r = round((r1+r2)/2);
        end;
        %     r = 1;
        while (r<=numel(s))
            cursol = u(:,1:r)*diag(s(1:r))*(v(:,1:r)');
            if (rx1*n1*n2*rx3<max_full_size)
                res = norm(B*cursol(:)-rhs)/normy;
            else
                % old
                %             res = norm(bfun2(B, cursol, rx1, n1, n2, rx3, rx1, n1, n2, rx3)-rhs)/normy;
                % new
                res = norm(bfun3(Phi1, a1, a2, Phi2, cursol)-rhs)/normy;
            end;
            if (res<max(real_tol, res_true*resid_damp))
                break;
            end;
            r = r+1;
        end;
    end;
    r = min(r, numel(s));
    if (verb>1)
        fprintf('dx: %3.3e, res: %3.3e, res_prev: %3.3e, r: %d\n', dx, res, res_prev, r);
    end;

    s = diag(s(1:r));
    u = u(:,1:r);
    v = conj(v(:,1:r));
else
    % implement tt-gmres here
    B = cell(2,1);
    B{1} = a1;
    B{2} = a2;
%     iB = tt_minres_selfprec(B, 1e-1, 1e-2, 10, 'right');
    iB = [];
    sol_prev = cell(2,1);
    sol_prev{1} = x1;
    sol_prev{2} = x2;
    rhs = cell(2,1);
    rhs{1} = y1;
    rhs{2} = y2;
    normy = tt_dist3(rhs, tt_scal(rhs,0));
    drhs = tt_mv(B, sol_prev);
    res_prev = tt_dist3(drhs, rhs)/normy;
    drhs = tt_add(rhs, tt_scal(drhs, -1));
    drhs = tt_compr2(drhs, real_tol);
    dsol = tt_gmres(B, drhs, real_tol*resid_damp_loc/res_prev, gmres_iters, nrestart, real_tol, real_tol, iB, [], [], [], 1);
%     sol = tt_gmres(B, rhs, real_tol*2, gmres_iters*5, nrestart/5, real_tol, real_tol, iB, [], [], sol_prev, 1);
    sol = tt_add(sol_prev, dsol);
    sol = tt_compr2(sol, real_tol);
    normsol = tt_dist3(sol, tt_scal(sol,0));
    dx = tt_dist3(sol, sol_prev)/normsol;
    res = tt_dist3(tt_mv(B, sol), rhs)/normy;

    dx_max = max(dx_max, dx);
    res_max = max(res_max, res_prev);
    [v, s]=qr(sol{2}, 0);
    sol{1} = sol{1}*(s.');
    [u, s]=qr(sol{1}, 0);
    r = size(sol{1},2);
    if (verb>1)
        fprintf('dx: %3.3e, res: %3.3e, res_prev: %3.3e, r: %d\n', dx, res, res_prev, r);
    end;
end;
end

author: oseledets

repo: TT-Toolbox


text:
function [T, score_mat] = learn_struct_mwst(data, discrete, node_sizes, node_type, scoring_fn, root)
% LEARN_STRUCT_MWST Learn an oriented tree using the MSWT algorithm
% T = learn_struct_mwst(data, discrete, node_sizes, node_type, scoring_fn, root)
%
% Input :
%   data(i,m) is the node i in the case m,
%   discrete = [ 1 if discret-node 0 if not ],
%   node_sizes = 1 if gaussian node,
%   node_type = {'tabular','gaussian',...},
%   score = 'bic' (for complete data and any node types) or 'mutual_info' (tabular nodes),
%   root is the futur root-node of the tree T.
%
% Output :
%	T = adjacency matrix of the tree
%
% V1.2 : 17 feb 2003 (O. Francois - francois.olivier.c.h@gmail.com, Ph. Leray - philippe.leray@univ-nantes.fr)
%
%
% See Chow&Liu 1968 for the original algorithm using Mutual Information scoring.
% Or Heckerman 1994.

if nargin <4
    error('Requires at least 4 arguments.')
end

if nargin == 4
    scoring_fn='bic'; root=1;
end;

if nargin == 5
    root=1;
end;

N=size(data,1);
score_mat=zeros(N,N);

switch scoring_fn
case 'bic',
    for i=1:(N-1)
            score2 = score_family(i, [], node_type{i}, scoring_fn, node_sizes, discrete, data,[]);
        for j=(i+1):N
            score1 = score_family(i, [j], node_type{i}, scoring_fn, node_sizes, discrete, data,[]);
            score = score2-score1;
            score_mat(i,j)=score;
            score_mat(j,i)=score;
        end
    end
case 'mutual_info',
    for i=1:(N-1)
        for j=(i+1):N
            score_mat(i,j)= -mutual_info_score(i,node_sizes(i),j,node_sizes(j),data);
            score_mat(j,i)=score_mat(i,j);
        end
    end
otherwise,
    error(['unrecognized scoring fn ' scoring_fn]);
end

G = minimum_spanning_tree(score_mat);
T = mk_rooted_tree(G, root);
T=full(T);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [varargout]=faces2elementSideSets(varargin)
% function [elementId,elementSides,elemSets,elemSetSides]=faces2elementSideSets(f,E)
% ------------------------------------------------------------------------
%
% Kevin Mattheus Moerman
% 2020/05/30 Created
% ------------------------------------------------------------------------

%%

switch nargin
    case 2
        f=varargin{1};
        E=varargin{2};
        reorderOpt=0;
    case 3
        f=varargin{1};
        E=varargin{2};
        reorderOpt=varargin{3};
end
%%
logicTouch=sum(ismember(E,f),2)==size(f,2);
ind=(1:1:size(E,1))';
[Fs,C,CF]=element2patch(E(logicTouch,:),ind(logicTouch));
logicMember=isrowmember(Fs,f,reorderOpt);
elementId=C(logicMember);
elementSides=CF(logicMember);
varargout{1}=elementId;
varargout{2}=elementSides;
if nargout>2
    setSideUni=unique(elementSides);
    numSets=numel(setSideUni);
    elemSets=cell(1,numSets);
    elemSetSides=zeros(1,numSets);
    for q=1:1:numSets
        sideType=setSideUni(q);
        logicNow=elementSides==sideType;
        indElem=unique(elementId(logicNow));
        elemSets{q}=indElem(:)';
        elemSetSides(q)=sideType;
    end
    varargout{3}=elemSets;
    varargout{4}=elemSetSides;
end

end
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
classdef SettingsLevelSetSmoothRectangleInclusion < SettingsLevelSetCreator

    properties (Access = public)
        widthH
        widthV
        pnorm
    end

    methods (Access = public)

        function obj = SettingsLevelSetSmoothRectangleInclusion(varargin)
            obj.loadParams('paramsLevelSetCreator_SmoothRectangle.json')
            if nargin == 1
                obj.loadParams(varargin{1})
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [ structout ] = bz_CollapseStruct( structin,dim,combine,NEST )
%structout = CollapseStruct( structin,dim,combine,NEST ) Combines elements in a
%structure array
%
%INPUT
%   structin    struct(N).fields structure array with N elements where each
%               of the N elements has the same fields and field structure.
%               structin can have nested cell arrays or stuctures, but be
%               careful with this... may not work for >1-dimensional
%               structures/cell arrays
%   dim         dimension along which to combine each element in the
%               structure array. (default: 2)
%               option: 'match'. CollapseStruct will find the dimension
%               along which the elements are concatenatable. If there are
%               multiple, will pick the first one.
%   (optional)
%       combine     can take 'mean' or 'median' instead of concatenating
%                   or 'std'. default: 'justcat'
%       NEST        true/false - go into nested structures?
%
%
%% A place for input options
if ~exist('NEST','var')
    NEST = false;
end

if ~exist('combine','var')
    combine = 'justcat';
end

if ~exist('dim','var')
    dim = 2;
end

%%
if ~isstruct(structin)
    structout=structin;
    warning('Not a structure...')
    return
end
fields = fieldnames(structin);

%%

for ff = 1:length(fields)
    currentfield = fields{ff};

	if isstruct(structin(1).(currentfield)) & NEST %For Nested Structures
        %Here: need to match fields, either 'remove' or 'addempty'
        %bz_Matchfields
       %structout.(currentfield) = cat(1,structin(:).(currentfield));
       structout.(currentfield) = bz_Matchfields({structin(:).(currentfield)},[],'remove');
       structout.(currentfield) = bz_CollapseStruct(structout.(currentfield),dim,combine,true);
       continue
    elseif iscell(structin(1).(currentfield)) & NEST %For cell array in field
        if strcmp(dim,'match')
            catdim = bz_FindCatableDims({structin(:).(currentfield)});
            if length(catdim)>1
                catdim = catdim(1); %Change this to be the largest dimension
            end
        else
            catdim = dim;
        end
        try
            structout.(currentfield) = cat(catdim,structin(:).(currentfield));
        catch
            display(['Failed to concatenate field ',currentfield])
        end
    elseif (isstring(structin(1).(currentfield))||ischar(structin(1).(currentfield))) & NEST %For string in field
        structout.(currentfield) = {structin(:).(currentfield)};
    else %For simple array in field
        try
            %Concatenate Array
            if strcmp(dim,'match')
                catdim = bz_FindCatableDims({structin(:).(currentfield)});
            else
                catdim = dim;
            end
        structout.(currentfield) = cat(catdim(1),structin(:).(currentfield));
        catch
             display(['Failed to concatenate field ',currentfield])
%            keyboard
%             continue
         end
	end

        if ~isfield(structout,currentfield)
            continue
        end
    switch combine
        case 'mean'
            structout.(currentfield) = nanmean(structout.(currentfield),catdim);
        case 'median'
            structout.(currentfield) = nanmedian(structout.(currentfield),catdim);
        case 'std'
            structout.(currentfield) = nanstd(structout.(currentfield),[],catdim);
        case 'sem'
            structout.(currentfield) =  nanstd(structout.(currentfield),[],catdim)./sqrt(sum(~isnan(structout.(currentfield)),catdim));
    end
end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function imsegs = processSuperpixelImage(fn)
% imsegs = processSuperpixelImage(fn)
% Creates the imsegs structure from a segmentation image
%
% INPUT:
% fn - 1 or 3 channel pixel maps of segmentation images.
% Segments are denoted by different RGB/Grayscale colors.
%
% OUTPUT:
% imsegs - image segmentation data
%
% Copyright(C) Derek Hoiem, Carnegie Mellon University, 2006

if ~iscell(fn)
    fn = {fn};
end

imsegs(length(fn)) = struct('imname', '', 'imsize', [0 0]);
for f = 1:length(fn)
    im = double(fn{f});

    imsegs(f).imname = [];
    imsegs(f).imsize = size(im);
    imsegs(f).imsize = imsegs(f).imsize(1:2);
    if(size(im,3)==3) % i.e. the raw output of FH segmentation
       im = im(:, :, 1) + im(:, :, 2)*256 + im(:, :, 3)*256^2;
    end

    [gid, gn] = grp2idx(im(:));
    imsegs(f).segimage = uint16(reshape(gid, imsegs(f).imsize));
    imsegs(f).nseg = length(gn);
end
imsegs = APPgetSpStats(imsegs);

author: Cloud-CV

repo: object-proposals


text:
function [ G ] = gsp_update_weights( G, W )
%GSP_UPDATE_WEIGHTS update weight matrix of the graph G
%   Usage: G = gsp_update_weight( G, W );
%
%   Input parameters
%       G   : Graph
%       W   : W new weight matrix
%   Output parameters
%       G   : Graph
%
%   This function will update the weight of graph and recompute the
%   Laplacian the new Laplacian. It will also recompute the value *G.lmax*
%   if it is present in the graph.
%
%   For now this function does not recompute the Fourier basis!
%
%   Example::
%
%          G.W = Wnew;
%          G = gsp_graph_default_parameters( G );
%
%
%   See also: gsp_graph_default_parameters gsp_graph
%

% Author: Nathanael Perraudin
% Date  : 19.11.2014

G.W = W;

if isfield(G,'U')
    G = rmfield(G,'U');
end

if isfield(G,'e')
    G =rmfield(G,'e');
end

if isfield(G,'mu')
    G =rmfield(G,'mu');
end

G = gsp_graph_default_parameters(G);

if isfield(G,'lmax')
    G = rmfield(G,'lmax');
    G = gsp_estimate_lmax(G);
end

end

repo: gspbox


text:
function c = r8mat_add ( m, n, alpha, a, beta, b )

%*****************************************************************************80
%
%% R8MAT_ADD computes C = alpha * A + beta * B for R8MAT's.
%
%  Discussion:
%
%    An R8MAT is an array of R8 values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real ALPHA, the multiplier for A.
%
%    Input, real A(M,N), the first matrix.
%
%    Input, real BETA, the multiplier for A.
%
%    Input, real B(M,N), the second matrix.
%
%    Output, real C(M,N), the sum of alpha*A+beta*B.
%
  c(1:m,1:n) = alpha * a(1:m,1:n) + beta * b(1:m,1:n);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [xi,w]=thirdOrderExpCubPoints(numDim,algorithm)
%%THIRDORDEREXPCUBPOINTS Generate third-order cubature points for
%               integration over real space involving the weighting
%               function w(x)=exp(-sqrt(sum(x.*x))).
%
%INPUTS: numDim An integer specifying the dimensionality of the points to
%               be generated.
%     algorithm A value indicating which algorithm should be used. Possible
%               values are:
%               0 (The default if omitted or an empty matrix is passed)
%                 Formula E_n^r 3-1 in [1], pg. 329, 2*numDim points.
%               1 Formula E_n^r 3-2 in [1], pg. 329, 2^numDim points.
%
%OUTPUTS: xi A numDim X numCubaturePoints matrix containing the cubature
%            points. (Each "point" is a vector)
%          w A numCubaturePoints X 1 vector of the weights associated with
%            the cubature points.
%
%REFERENCES:
%[1] A.H. Stroud, Approximate Calculation of Multiple Integrals. Cliffs,
%    NJ: Prentice-Hall, Inc., 1971.
%
%October 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release

if(nargin<2||isempty(algorithm))
   algorithm=0;
end

switch(algorithm)
    case 0%E_n^r 3-1 in [1], pg. 329, 2*numDim points.
        V=2*pi^(numDim/2)*exp(gammaln(numDim)-gammaln(numDim/2));
        r=sqrt(numDim*(numDim+1));

        w=V/(2*numDim)*ones(2*numDim,1);
        xi=fullSymPerms([r;zeros(numDim-1,1)]);
    case 1%E_n^r 3-2 in [1], pg. 329, 2^numDim points.
        V=2*pi^(numDim/2)*exp(gammaln(numDim)-gammaln(numDim/2));
        r=sqrt(numDim+1);

        xi=PMCombos(r*ones(numDim,1));
        w=2^(-numDim)*V*ones(2^numDim,1);
    otherwise
        error('Unknown algorithm specified');
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%LOAD_ACQ  load BIOPAC's AcqKnowledge file format for Windows PC
%
%  Usage: acq = load_acq(filename, [force_chan_by_chan_load])
%
%  acq - AcqKnowledge file structure containing ACQ header field,
%	and data matrix.
%
%  filename - BIOPAC's AcqKnowledge file
%
%  force_chan_by_chan_load - Optional. By default, this optional flag will
%	be set to 1, which means that when you use acq = load_acq(fname),
%	the data will be loaded one channel after the other. This can
%	avoid memory crash when you load very large ACQ data. If your
%	ACQ data is not huge, I suggest that you set this optional flag to
%	0, i.e. acq = load_acq(fname, 0). In this case, the program will
%	read data depending on the data type. If the program detects that
%	the data type in ACQ file are different from channel to channel,
%	it will still read data channel by channel. Otherwise, it will
%	read whole data in one block (a lot faster than using traditional
%	way from channel to channel with the same result).
%
%  Notes:
%
%  This program is based on Application Note #156 from BIOPAC web site:
%  http://www.biopac.com/ResearchNotes.asp?Aid=&ANid=82&Level=4
%  ( 156 - ACQKNOWLEDGE FILE FORMATS FOR PC WITH WINDOWS )
%
%  The note mentioned that: "This document describes file formatting for
%  all Windows versions of AcqKnowledge 3.9.x or below". Thanks to the
%  open Python source code provided by Nathan Vack, this program can also
%  read AcqKnowledge 4.0 & 4.1 data (with no documentation from BIOPAC).
%  Compressed data is not supported by this program.
%
%  Created on 5-APR-2007 by Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
%  Modify on 31-JUL-2007 by Julio Cruz (julio.cruz@juliocruz.info) to
%	make the program also work on MAC OS.
%
%  Modify on 01-MAY-2008 by Antonio Molins (amolins@mit.edu) to
%	make the program work with files recorded in Biopac student data
%	acquisition system. Changes include:
%	- added reading of markers, with "markers" field included in
%	returned structure..
%	- proper handling of foreign data (was not getting quite there with
%	the code (as it was downloaded from MathWorks).
%	- proper handling of var_sampling_rate = 0, as produced by the
%	Biopac software in occasions.
%
function acq = load_acq(filename, chan_by_chan)

   if ~exist('chan_by_chan','var')
      chan_by_chan = 1;
   end

   if ~exist('filename','var')
      error('Usage: acq = load_acq(filename)');
   end

   if ~exist(filename,'file')
      error([filename, ': Can''t open file']);
   end

   fid = fopen(filename,'r', 'l');

   if fid < 0,
      msg = sprintf('Cannot open file %s.',filename);
      error(msg);
   end

   fread(fid, 1, 'int16')';
   file_version = fread(fid, 1, 'int32')';

   %  try different endian
   %
   if file_version < 0 | file_version > 200
      fclose(fid);
      fid = fopen(filename,'r', 'b');

      fread(fid, 1, 'int16')';
      file_version = fread(fid, 1, 'int32')';
      fseek(fid, 0, 'bof');

      if file_version < 0 | file_version > 200
         error('This ACQ file is not supported');
      end
   end

   fprintf('Loading %s ', filename);

   %  read header
   %
   acq.hdr = read_acq_hdr(fid);

   %  read data
   %
   acq.data = read_acq(fid, acq.hdr, chan_by_chan);

   %  read markers written by AM (Antonio Molins) for ACQ 3
   %
   if file_version < 68
      acq.markers = read_markers(fid, acq.hdr, size(acq.data));
   end

   fclose(fid);

   fprintf(' Done!\n');

   return;					% load_acq

%---------------------------------------------------------------------
function hdr = read_acq_hdr(fid)

   fseek(fid, 0, 'bof');
   hdr.graph = graph(fid);
   acc_chan_header_len = 0;

   for i = 1:hdr.graph.num_channels
      fseek(fid, hdr.graph.ext_item_header_len+acc_chan_header_len, 'bof');
      hdr.per_chan_data(i) = per_chan_data(fid, hdr.graph.file_version);
      acc_chan_header_len = acc_chan_header_len + hdr.per_chan_data(i).chan_header_len;
   end

   % added by AM: foreign data section was being started to read one short,
   % this fseek takes care of that
   fseek(fid, hdr.graph.ext_item_header_len+acc_chan_header_len, 'bof');

   hdr.foreign = foreign(fid, hdr.graph.file_version);

   if hdr.graph.file_version >= 68 & hdr.graph.file_version < 80
      unused = fread(fid, hdr.foreign.length2-12, 'uint8')';
   end

   for i = 1:hdr.graph.num_channels
      hdr.per_chan_type(i) = per_chan_type(fid);
   end

   return;					% read_acq_hdr

%---------------------------------------------------------------------
function hdr = graph(fid, file_version)

   %  Struct						% off + size
   unused = fread(fid, 1, 'int16')';			% 0 + 2
   hdr.file_version = fread(fid, 1, 'int32');		% 2 + 4
   hdr.ext_item_header_len = fread(fid, 1, 'int32');	% 6 + 4
   hdr.num_channels = fread(fid, 1, 'int16');		% 10 + 2
   hdr.horiz_axis_type = fread(fid, 1, 'int16');	% 12 + 2
   hdr.curr_channel = fread(fid, 1, 'int16');		% 14 + 2
   hdr.sample_time = fread(fid, 1, 'double');		% 16 + 8
   hdr.time_offset = fread(fid, 1, 'double');		% 24 + 8
   hdr.time_scale = fread(fid, 1, 'double');		% 32 + 8
   hdr.time_cursor1 = fread(fid, 1, 'double');		% 40 + 8
   hdr.time_cursor2 = fread(fid, 1, 'double');		% 48 + 8
   hdr.chart_window = fread(fid, 4, 'int16');		% 56 + 8
   hdr.measurement = fread(fid, 6, 'int16');		% 64 + 12
   hdr.hilite = fread(fid, 1, 'int16');			% 76 + 2
   hdr.first_time_offset = fread(fid, 1, 'double');	% 78 + 8
   hdr.rescale = fread(fid, 1, 'int16');		% 86 + 2
   hdr.horiz_units1 = deblank(fread(fid, 40, '*char')'); % 88 + 40
   hdr.horiz_units2 = deblank(fread(fid, 10, '*char')'); % 128 + 10
   hdr.in_memory = fread(fid, 1, 'int16');		% 138 + 2
   hdr.grid = fread(fid, 1, 'int16');			% 140 + 2
   hdr.markers = fread(fid, 1, 'int16');		% 142 + 2
   hdr.plot_draft = fread(fid, 1, 'int16');		% 144 + 2
   hdr.display_mode = fread(fid, 1, 'int16');		% 146 + 2
   hdr.reserved = fread(fid, 1, 'int16');		% 148 + 2

   if hdr.file_version > 33 & hdr.file_version < 68
      hdr.show_toolbar = fread(fid, 1, 'int16');	% 150 + 2
      hdr.show_chan_butt = fread(fid, 1, 'int16');	% 152 + 2
      hdr.show_measurement = fread(fid, 1, 'int16');	% 154 + 2
      hdr.show_marker = fread(fid, 1, 'int16');		% 156 + 2
      hdr.show_journal = fread(fid, 1, 'int16');	% 158 + 2
      hdr.cur_x_channel = fread(fid, 1, 'int16');	% 160 + 2
      hdr.mmt_precision = fread(fid, 1, 'int16');	% 162 + 2
   end

   if hdr.file_version > 34 & hdr.file_version < 68
      hdr.measurement_row = fread(fid, 1, 'int16');	% 164 + 2
      hdr.mmt = fread(fid, 40, 'int16');		% 166 + 80
      hdr.mmt_chan = fread(fid, 40, 'int16');		% 246 + 80
   end

   if hdr.file_version > 35 & hdr.file_version < 68
      hdr.mmt_calc_opnd1 = fread(fid, 40, 'int16');	% 326 + 80
      hdr.mmt_calc_opnd2 = fread(fid, 40, 'int16');	% 406 + 80
      hdr.mmt_calc_op = fread(fid, 40, 'int16');	% 486 + 80
      hdr.mmt_calc_constant = fread(fid, 40, 'double');	% 566 + 320
   end

   if hdr.file_version > 37 & hdr.file_version < 68
      hdr.new_grid_minor = fread(fid, 1, 'int32');	% 886 + 4
      hdr.color_major_grid = fread(fid, 1, 'int32');	% 890 + 4
      hdr.color_minor_grid = fread(fid, 1, 'int32');	% 894 + 4
      hdr.major_grid_style = fread(fid, 1, 'int16');	% 898 + 2
      hdr.minor_grid_style = fread(fid, 1, 'int16');	% 900 + 2
      hdr.major_grid_width = fread(fid, 1, 'int16');	% 902 + 2
      hdr.minor_grid_width = fread(fid, 1, 'int16');	% 904 + 2
      hdr.fixed_units_div = fread(fid, 1, 'int32');	% 906 + 4
      hdr.mid_range_show = fread(fid, 1, 'int32');	% 910 + 4
      hdr.start_middle_point = fread(fid, 1, 'double');	% 914 + 8
      hdr.offset_point = fread(fid, 60, 'double');	% 922 + 480
      hdr.h_grid = fread(fid, 1, 'double');		% 1402 + 8
      hdr.v_grid = fread(fid, 60, 'double');		% 1410 + 480
      hdr.enable_wave_tool = fread(fid, 1, 'int32');	% 1890 + 4

      %  interpret color_major_grid
      %
      color_str = sprintf('%06s', dec2hex(hdr.color_major_grid));
      hdr.color_major_grid = ...
         [hex2dec(color_str(5:6)) hex2dec(color_str(3:4)) hex2dec(color_str(1:2))]/255;

      %  interpret color_minor_grid
      %
      color_str = sprintf('%06s', dec2hex(hdr.color_minor_grid));
      hdr.color_minor_grid = ...
         [hex2dec(color_str(5:6)) hex2dec(color_str(3:4)) hex2dec(color_str(1:2))]/255;

   end

   if hdr.file_version > 38 & hdr.file_version < 68
      hdr.horiz_precision = fread(fid, 1, 'int16');	% 1894 + 2
   end

   if hdr.file_version > 40 & hdr.file_version < 68
      hdr.reserved2 = fread(fid, 20, 'int8');		% 1896 + 20
      hdr.overlap_mode = fread(fid, 1, 'int32');	% 1916 + 4
      hdr.show_hardware = fread(fid, 1, 'int32');	% 1920 + 4
      hdr.x_auto_plot = fread(fid, 1, 'int32');	% 1924 + 4
      hdr.x_auto_scroll = fread(fid, 1, 'int32');	% 1928 + 4
      hdr.start_butt_visible = fread(fid, 1, 'int32');	% 1932 + 4
      hdr.compressed = fread(fid, 1, 'int32');		% 1936 + 4
      hdr.always_start_butt_visible = fread(fid, 1, 'int32');	% 1940 + 4
   end

   if hdr.file_version > 42 & hdr.file_version < 68
      hdr.path_video = deblank(fread(fid, 260, '*char')'); % 1944 + 260
      hdr.opt_sync_delay = fread(fid, 1, 'int32');	% 2204 + 4
      hdr.sync_delay = fread(fid, 1, 'double');		% 2208 + 8
      hdr.hrp_paste_measurement = fread(fid, 1, 'int32'); % 2216 + 4
   end

   if hdr.file_version > 44 & hdr.file_version < 68
      hdr.graph_type = fread(fid, 1, 'int32');		% 2220 + 4
      hdr.mmt_calc_expr = fread(fid, [40 256], '*char'); % 2224 + 10240
      hdr.mmt_moment_order = fread(fid, 40, 'int32');	% 12464 + 160
      hdr.mmt_time_delay = fread(fid, 40, 'int32');	% 12624 + 160
      hdr.mmt_embed_dim = fread(fid, 40, 'int32');	% 12784 + 160
      hdr.mmt_mi_delay = fread(fid, 40, 'int32');	% 12944 + 160
   end

   if hdr.file_version >= 68
      unused = fread(fid, 411, 'int16')';		% 150 + 822
      hdr.compressed = fread(fid, 1, 'int32');		% 972 + 4
   end

   return;						% graph

%---------------------------------------------------------------------
function hdr = per_chan_data(fid, file_version)

   %  Struct						% off + size
   hdr.chan_header_len = fread(fid, 1, 'int32')';	% 0 + 4
   hdr.num = fread(fid, 1, 'int16')';			% 4 + 2
   hdr.comment_text = deblank(fread(fid, 40, '*char')'); % 6 + 40
   hdr.rgb_color = fread(fid, 1, 'int32')';		% 46 + 4
   hdr.disp_chan = fread(fid, 1, 'int16')';		% 50 + 2
   hdr.volt_offset = fread(fid, 1, 'double')';		% 52 + 8
   hdr.volt_scale = fread(fid, 1, 'double')';		% 60 + 8
   hdr.units_text = deblank(fread(fid, 20, '*char')');	% 68 + 20
   hdr.buf_length = fread(fid, 1, 'int32')';		% 88 + 4
   hdr.ampl_scale = fread(fid, 1, 'double')';		% 92 + 8
   hdr.ampl_offset = fread(fid, 1, 'double')';		% 100 + 8
   hdr.chan_order = fread(fid, 1, 'int16')';		% 108 + 2
   hdr.disp_size = fread(fid, 1, 'int16')';		% 110 + 2

   if file_version >= 68

      unused = fread(fid, 5, 'double')';		% 112 + 40
      hdr.var_sample_divider = fread(fid, 1, 'int16')'; % 152 + 2

      % FIX according to comments on fileexchange
      hdr.var_sample_divider = 1;

   else

   hdr.plot_mode = fread(fid, 1, 'int16')';		% 112 + 2
   hdr.mid = fread(fid, 1, 'double')';			% 114 + 8

   %  interpret rbg_color
   %
   color_str = sprintf('%06s', dec2hex(hdr.rgb_color));
   hdr.rgb_color = ...
      [hex2dec(color_str(5:6)) hex2dec(color_str(3:4)) hex2dec(color_str(1:2))]/255;

   if file_version > 37
      hdr.description = deblank(fread(fid, 128, '*char')'); % 122 + 128
      hdr.var_sample_divider = fread(fid, 1, 'int16')'; % 250 + 2

      % AM does not make sense to have a zero divider, set those to 1. No
      % warranties this is standard, but .ACQ files generated with biopac
      % programs produce files with var_sample_divider = 0 although this is
      % not documented. This fix works for the files seen so far.
      if hdr.var_sample_divider==0
          hdr.var_sample_divider = 1;
      end
   else
      hdr.var_sample_divider = 1;
   end

   if file_version > 38
      hdr.vert_precision = fread(fid, 1, 'int16')';	% 252 + 2
   end

   if file_version > 42
      hdr.active_seg_color = fread(fid, 1, 'int32')';	% 254 + 4
      hdr.active_seg_style = fread(fid, 1, 'int32')';	% 258 + 4

      %  interpret active_seg_color
      %
      color_str = sprintf('%06s', dec2hex(hdr.active_seg_color));
      hdr.active_seg_color = ...
         [hex2dec(color_str(5:6)) hex2dec(color_str(3:4)) hex2dec(color_str(1:2))]/255;
   end

   end						% if file_version < 68

   return;						% per_chan_data

%---------------------------------------------------------------------
function hdr = foreign(fid, file_version)

   %  Struct						% off + size
   hdr.length = fread(fid, 1, 'short')';		% 0 + 2
   hdr.id = fread(fid, 1, 'short')';			% 2 + 2

   if file_version < 68
      hdr.by_foreign_data = fread(fid, hdr.length-4, 'int8')'; % 4 + x
      hdr.length2 = hdr.length;
   elseif file_version < 80
      unused = fread(fid, 1, 'int32')';			% 4 + 4
      hdr.length2 = fread(fid, 1, 'int32')' + 8;	% 8 + 4
   else
      unused = fread(fid, 1, 'int32')';			% 4 + 4
      hdr.length2 = hdr.length + 8;
   end

   return;						% foreign

%---------------------------------------------------------------------
function hdr = per_chan_type(fid)

   %  Struct						% off + size
   hdr.size = fread(fid, 1, 'short')';			% 0 + 2
   hdr.type = fread(fid, 1, 'short')';			% 2 + 2

   return;						% per_chan_type

%---------------------------------------------------------------------
function data = read_acq(fid, hdr, chan_by_chan)

   start_real_chan = hdr.graph.ext_item_header_len;

   for i = 1:hdr.graph.num_channels
      start_real_chan = start_real_chan + hdr.per_chan_data(i).chan_header_len;

      if hdr.per_chan_type(i).type ~= 1			% if integer
         hdr.per_chan_type(i).size = 2;			% only use int16
      end
   end

   start_real_chan = start_real_chan + hdr.foreign.length2 + 4*hdr.graph.num_channels;

   sample_divider = [hdr.per_chan_data.var_sample_divider];

   if length(unique(sample_divider))==1 & unique(sample_divider)==1
      min_len = min([hdr.per_chan_data.buf_length]);
   else
      min_len = min([hdr.per_chan_data.buf_length].*sample_divider);
   end

   if length(unique([hdr.per_chan_type.type])) & ~chan_by_chan & ...
	length(unique(sample_divider))==1 & unique(sample_divider)==1

      half_chan = round(hdr.graph.num_channels/2);

      for i = 1:half_chan
         fprintf('.');
      end

      if hdr.per_chan_type(i).type == 1			% double

         data=fread(fid,[hdr.graph.num_channels min_len],'double');
         data=data';
      else						% int

         data=fread(fid,[hdr.graph.num_channels min_len],'int16');
         data=data'.*(ones(min_len,1)*[hdr.per_chan_data.ampl_scale]) ...
		    + ones(min_len,1)*[hdr.per_chan_data.ampl_offset];
      end

      for i = 1:(hdr.graph.num_channels-half_chan)
         fprintf('.');
      end

   else				% if we have to do it chan_by_chan

      if length(unique(sample_divider))==1 & unique(sample_divider)==1

         %  Since data are arranged like: "channel in sample"
         %  { s1 of {ch1 ch2 ...}, s2 of {ch1 ch2 ...} ... }
         %  We can either read data sample by sample (all chan at once),
         %  or, channel by channel, but need to skip the rest of chan
         %
         size_all_chan_per_sample = 0;

         for i = 1:hdr.graph.num_channels
            size_all_chan_per_sample = size_all_chan_per_sample + ...
					hdr.per_chan_type(i).size;
         end

         for i = 1:hdr.graph.num_channels

            fprintf('.');
            fseek(fid, start_real_chan, 'bof');

            %  First jump to the start point of the right channel
            %
            start_chan = 0;

            for j = 1 : (i-1)
               start_chan = start_chan + hdr.per_chan_type(j).size;
            end

            fseek(fid, start_chan, 'cof');

            %  Need to skip the rest of chan, in order to read each sample point
            %
            skip_chan = size_all_chan_per_sample - hdr.per_chan_type(i).size;

            if hdr.per_chan_type(i).type == 1		% double

               tmp=fread(fid,hdr.per_chan_data(i).buf_length,'double',skip_chan);

            else					% int

               %  Need to be scaled & shifted by ampl_scale & ampl_offset
               %  for integer data
               %
               tmp=fread(fid,hdr.per_chan_data(i).buf_length,'int16',skip_chan) ...
			* hdr.per_chan_data(i).ampl_scale ...
			+ hdr.per_chan_data(i).ampl_offset;
            end

            data(:,i) = tmp(1:min_len);
         end

      else				% sample_divider>1 or different

         %  If there is a sample_divider>1 in any channel, we have to
         %  read data sample by sample (very slow!)
         %
         data = zeros(min_len, hdr.graph.num_channels);
         mask=zeros(min_len, hdr.graph.num_channels);
         sample_divider = abs(sample_divider);

         for j = 1:hdr.graph.num_channels
            mask(1:sample_divider(j):min_len, j)=1;
         end

         for i = 1:min_len

            if mod(i-1, ceil(min_len/hdr.graph.num_channels))==0
               fprintf('.');
            end

            for j = 1:hdr.graph.num_channels
               if mask(i,j) | sample_divider(j)==1
                  if hdr.per_chan_type(j).type == 1		% double
                     data(i,j) = fread(fid,1,'double');
                  else
                     data(i,j) = fread(fid,1,'int16');		% int

                     %  Need to be scaled & shifted by ampl_scale & ampl_offset
                     %  for integer data
                     %
                     data(i,j) = data(i,j) * hdr.per_chan_data(j).ampl_scale ...
					+ hdr.per_chan_data(j).ampl_offset;
                  end
               else
                  data(i,j) = data(i-1,j);
               end		% if mask ... read sample_divider

            end		% for j
         end		% for i

      end	% if length(unique(sample_divider
   end	% length(unique([hdr.per_chan_type.type

   return;					% read_acq

 %---------------------------------------------------------------------
 function info = read_markers(fid, hdr, data_size)

     % AM foolproof way of getting to the markers: see where data starts and
     % add the size of the data..

     start_real_chan = hdr.graph.ext_item_header_len;

     for i = 1:hdr.graph.num_channels
         start_real_chan = start_real_chan + hdr.per_chan_data(i).chan_header_len;

         if hdr.per_chan_type(i).type ~= 1			% if integer
             hdr.per_chan_type(i).size = 2;			% only use int16
         end
     end

     start_real_chan = start_real_chan + hdr.foreign.length2 + 4*hdr.graph.num_channels;

     size_channel_data = sum([hdr.per_chan_type.size])*data_size(1);
     start_markers = start_real_chan + size_channel_data;

     % AM place the file position in the position where data ends
     fseek(fid,start_markers,'bof');

     % AM then follow specifications; borrowed from code of the
     % non-functioning (at least for me) ACQREAD,
     %    ACQREAD, version 2.0 (2007-08-21)
     %    Copyright (C) 2006-2007  Sebastien Authier and Vincent Finnerty

     info.lLength = fread(fid,1,'*int32');
     info.lMarkers = fread(fid,1,'*int32');		% Number of markers
     if (info.lLength > 0) & (info.lMarkers > 0)
         for n = 1:double(info.lMarkers)
             info.lSample(n) = fread(fid,1,'*int32');	% Location of marker
             info.fSelected(n) = fread(fid,1,'*int16');
             info.fTextLocked(n) = fread(fid,1,'*int16');
             info.fPositionLocked(n) = fread(fid,1,'*int16');
             info.nTextLength(n) = fread(fid,1,'*int16');  % Length of marker text string
             info.szText{n} = deblank(fread(fid,double(info.nTextLength(n))+1,'*char')');  % Marker text string
         end
     else
         info.lSample = [];
         info.fSelected = [];
         info.fTextLocked = [];
         info.fPositionLocked = [];
         info.nTextLength = [];
         info.szText = [];
     end

     return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function favar = initializeFavarResults(opts)
% INITIALIZESTRCIDENT this function initalizes a strcident struct to store
% the results based on the already existing fields on the settings objet
arguments
    opts (1,1) bear.settings.BASEsettings
end

favar.FAVAR     = false;
favar.HDplot    = false;
favar.IRFplot   = false;
favar.FEVDplot  = false;

if isprop(opts, 'favar') && isa(opts.favar, 'bear.settings.favar.FAVARsettings') && opts.favar.FAVAR
    favar = initalizeFavar(opts.favar);
end
end

function favar = initalizeFavar(favarOpt)

p = properties(favarOpt);
for i = 1 : length(p)

    favar.(p{i}) = favarOpt.(p{i});

end

end

author: european-central-bank

repo: BEAR-toolbox


text:
function [Rmean, Pmean, frq, Pmean2] = NormedSingleCurveLengthWindowed(x,winlen,shiftlen,lag,fs,nrmdegree)

beVerbose = false;

% shiftlen is winlen - overlaplen

M = floor((length(x) - winlen) /  shiftlen ) + 1;

Psum = zeros(1,2*lag+1);
Rsum = zeros(1,lag+1);

for ii = 1:M

    xseg = x((ii-1)*shiftlen+1:(ii-1)*shiftlen+winlen);

    [rx,px] = NormedSingleCurveLength(xseg,lag,fs,0,nrmdegree);

    Psum = Psum + px;

    Rsum = Rsum + rx;

    if beVerbose
        display(ii);
    end
end

Pmean = Psum ./ M;
Rmean = Rsum ./ M;

Pmean2 = abs(fft([fliplr(Rmean) Rmean(2:end)]));

frq = linspace(0,0.5*fs,lag+1);

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [out] = evap_2(p1,S,Smax,Ep,dt)
%evap_2

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Evaporation at a scaled, plant-controlled rate
% Constraints:  f <= Ep
%               f <= S/dt
% @(Inputs):    p1   - plant-controlled base evaporation rate [mm/d]
%               S    - current storage [mm]
%               Smax - maximum storage [mm]
%               Ep   - potential evapotranspiration rate [mm/d]
%               dt   - time step size [d]

out = min([p1*S/Smax,Ep,S/dt]);

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text: %UPLUS Dataset overload

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
classdef Norm_DCF_Plot < handle
% This class is for plotting normalized DCF curves from the
% scores of a given detector, as a function of the effective target
% prior. The DCF is parametrized by the prior and has unity cost
% for both misses and false alarms. The algorithms to compute DCF
% have been optimized to allow the use of large score sets, in
% order to be able to obtain meaningful results at extreme
% operating points.
%
% There are six curves that can be plotted:
% (i) minimum DCF, where the evaluator chooses the optimal decision
%     threshold at every operating point. (The operating point is defined
%     by the effective target prior).
% (ii) actual DCF, where the decision threshold is the
%      minimum-probability-of-error Bayes decision threshold, assuming that
%      the given detector scores are log-likelihood-ratios. This decision
%      threshold is just -logit(Ptar).
% (iii) actual misses, the component of actual DCF due to miss errors.
% (iv) actual false alarms, the component of actual DCF due to false-alarm errors.
% (v) minimum misses, the component of minimum DCF due to miss errors.
% (vi) miniumu false alarms, the component of minimum DCF due to false-alarm errors.
%
% All curves are normalized by dividing by the DCF of the default system
% that makes decisions based on the prior alone. The normalization is
% min(Ptar,1-Ptar).
%
% The x-axis is logit(Ptar) = log(Ptar) - log(1-Ptar). The logit
% transformation maps the effective target prior to the whole
% extended real line from -inf to inf, but this plot is confined to
% the given limits xmin and xmax.
%
% For more information on DCF, type:
% > help fast_actDCF
% > help fast_minDCF
%
% Use set_system with the target and non-target scores for a system
% before calling the plotting functions.  The curves plotted will
% all be for the current set system.  Curves for more than one
% system can be plotted by calling set_system for the first system,
% plotting all of its curves and points of interest, calling
% set_system for the second system, plotting, etc.
%
% Inputs:
%    plot_axes: [xmin,xmax,ymin,ymax]: the range for the x-axis and
%      y-axis of this plot, where x = logit(Ptar).
%
%    plot_title: An optional title string for the plot.
%

properties (Access = private)
  fh
  plot_axes
  sys_name
  plo
  actDCF
  actPmiss
  actPfa
  minDCF
  minPmiss
  minPfa
  dr30Miss
  dr30FA
  Ptar_norm
  Pnon_norm
  handles_vec
  legend_strings = {};
end

methods (Access = public)
  set_system(plot_obj,tar,non,sys_name)
  set_system_from_scores(plot_obj,scores,key,sys_name)
  plot_fa_rate_min(plot_obj,plot_args,legend_string)
  plot_fa_rate_act(plot_obj,plot_args,legend_string)
  plot_miss_rate_min(plot_obj,plot_args,legend_string)
  plot_miss_rate_act(plot_obj,plot_args,legend_string)
  plot_dcf_curve_min(plot_obj,plot_args,legend_string)
  plot_dcf_curve_act(plot_obj,plot_args,legend_string)
  plot_DR30_fa(plot_obj,plot_args,legend_string)
  plot_DR30_miss(plot_obj,plot_args,legend_string)
  plot_DR30_both(plot_obj,plot_args_fa,plot_args_miss)
  plot_curves(plot_obj,mask,line_info)
  plot_operating_point(plot_obj,value,plot_args,legend_string)
  display_legend(plot_obj)
  save_as_pdf(plot_obj,outfilename)
end

methods (Access = public)
  % plot_obj = Norm_DCF_Plot()
  % plot_obj = Norm_DCF_Plot(plot_axes)
  % plot_obj = Norm_DCF_Plot(plot_axes,plot_title)
  function plot_obj = Norm_DCF_Plot(plot_axes,plot_title)
  % constructor
  if exist('plot_axes','var') && ~isempty(plot_axes)
    assert(length(plot_axes)==4)
    plot_obj.plot_axes = plot_axes;
  else
    plot_obj.plot_axes = [-10,0,0,1.2];
  end
  xmin = plot_obj.plot_axes(1);
  xmax = plot_obj.plot_axes(2);
  assert(xmin<xmax,'Illegal parameters xmin and xmax.')
  step = (xmax-xmin)/1000;
  plot_obj.plo = xmin:step:xmax;
  Ptar = sigmoid(plot_obj.plo);
  Pnon = sigmoid(-plot_obj.plo);
  refPe = min(Ptar,Pnon);
  plot_obj.Ptar_norm = Ptar ./ refPe;
  plot_obj.Pnon_norm = Pnon ./ refPe;
  plot_obj.actDCF = [];
  plot_obj.handles_vec = [];
  plot_obj.legend_strings = {};

  plot_obj.fh = figure();
  hold on
  ylabel('normalized DCF');
  xlabel('logit P_{tar}');
  grid
  axis(plot_obj.plot_axes);
  if exist('plot_title','var') && ~isempty(plot_title)
    title(plot_title);
  end
  end
end

methods (Access = private)
  add_legend_entry(plot_obj,lh,legend_string,append_name)
  plot_miss_rate(plot_obj,Pmiss,plot_args,legend_string)
  plot_fa_rate(plot_obj,Pfa,plot_args,legend_string)
  plot_dcf_curve(plot_obj,dcf,plot_args,legend_string)
end

end

author: nesl


text:
function varargout = JAABASplashScreen(varargin)
% JAABASPLASHSCREEN MATLAB code for JAABASplashScreen.fig
%      JAABASPLASHSCREEN, by itself, creates a new JAABASPLASHSCREEN or raises the existing
%      singleton*.
%
%      H = JAABASPLASHSCREEN returns the handle to a new JAABASPLASHSCREEN or the handle to
%      the existing singleton*.
%
%      JAABASPLASHSCREEN('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in JAABASPLASHSCREEN.M with the given input arguments.
%
%      JAABASPLASHSCREEN('Property','Value',...) creates a new JAABASPLASHSCREEN or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before JAABASplashScreen_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to JAABASplashScreen_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help JAABASplashScreen

% Last Modified by GUIDE v2.5 13-May-2012 20:22:42

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @JAABASplashScreen_OpeningFcn, ...
                   'gui_OutputFcn',  @JAABASplashScreen_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before JAABASplashScreen is made visible.
function JAABASplashScreen_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to JAABASplashScreen (see VARARGIN)

% Choose default command line output for JAABASplashScreen
handles.output = hObject;

try
  imname = 'janelia_logo.png';
  versionfile = 'version.txt';
  if exist(imname,'file'),
    im = imread(imname);
  else
    im = 0;
  end
  image(im,'Parent',handles.axes1);
  axis(handles.axes1,'image','off');
  if exist(versionfile,'file'),
    fid = fopen(versionfile,'r');
    while true,
      version = fgetl(fid);
      if ~ischar(version),
        break;
      end
      version = strtrim(version);
      if isempty(version),
        continue;
      end
      break;
    end
    fclose(fid);
  else
    version = '???';
  end

  set(handles.text_info,'String',sprintf('JAABA: The Janelia Automatic Animal Behavior Annotator\nversion %s',version));
end

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes JAABASplashScreen wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = JAABASplashScreen_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;
varargout{2} = handles.text_status;

author: kristinbranson

repo: JAABA


text:
function obj = updateIsosurface(obj, isoIndex)

    %-------------------------------------------------------------------------%

	%-INITIALIZATIONS-%

	axIndex = obj.getAxisIndex(obj.State.Plot(isoIndex).AssociatedAxis);
	plotData = get(obj.State.Plot(isoIndex).Handle);
	axisData = get(plotData.Parent);
	[xSource, ySource] = findSourceAxis(obj, axIndex);

	%-update scene-%
	updateScene(obj, isoIndex)

	%-get mesh data-%
	xData = plotData.Vertices(:, 1);
	yData = plotData.Vertices(:, 2);
	zData = plotData.Vertices(:, 3);

	iData = plotData.Faces(:, 1) - 1;
	jData = plotData.Faces(:, 2) - 1;
	kData = plotData.Faces(:, 3) - 1;

	%-------------------------------------------------------------------------%

	%-get trace-%
	obj.data{isoIndex}.type = 'mesh3d';
	obj.data{isoIndex}.name = plotData.DisplayName;
	obj.data{isoIndex}.showscale = false;

	%-------------------------------------------------------------------------%

	%-set mesh data-%
	obj.data{isoIndex}.x = xData;
	obj.data{isoIndex}.y = yData;
	obj.data{isoIndex}.z = zData;

	obj.data{isoIndex}.i = iData;
	obj.data{isoIndex}.j = jData;
	obj.data{isoIndex}.k = kData;

	%-------------------------------------------------------------------------%

	%-mesh coloring-%
	faceColor = getFaceColor(plotData, axisData);

	if iscell(faceColor)
		obj.data{isoIndex}.facecolor = faceColor;
	else
		obj.data{isoIndex}.color = faceColor;
	end

	%-------------------------------------------------------------------------%

	%-lighting settings-%
	if ~strcmp(plotData.FaceLighting, 'flat')
		obj.data{isoIndex}.lighting.diffuse = plotData.DiffuseStrength;
		obj.data{isoIndex}.lighting.ambient = plotData.AmbientStrength;
		obj.data{isoIndex}.lighting.specular = plotData.SpecularStrength;
		obj.data{isoIndex}.lighting.roughness = 0.2;
		obj.data{isoIndex}.lighting.fresnel = 0.5;
		obj.data{isoIndex}.lighting.vertexnormalsepsilon = 1e-12;
		obj.data{isoIndex}.lighting.facenormalsepsilon = 1e-6;
	end

	%-------------------------------------------------------------------------%

	%-associate scene to trace-%
	obj.data{isoIndex}.scene = sprintf('scene%d', xSource);

	%-------------------------------------------------------------------------%
end

function updateScene(obj, isoIndex)

	%-INITIALIZATIONS-%
	axIndex = obj.getAxisIndex(obj.State.Plot(isoIndex).AssociatedAxis);
	plotData = get(obj.State.Plot(isoIndex).Handle);
	axisData = get(plotData.Parent);
	[xSource, ySource] = findSourceAxis(obj, axIndex);
	scene = eval( sprintf('obj.layout.scene%d', xSource) );

	aspectRatio = axisData.PlotBoxAspectRatio;
	cameraPosition = axisData.CameraPosition;
	dataAspectRatio = axisData.DataAspectRatio;
	cameraUpVector = axisData.CameraUpVector;
	cameraEye = cameraPosition./dataAspectRatio;
	normFac = 0.5*abs(min(cameraEye));

	%-aspect ratio-%
	scene.aspectratio.x = aspectRatio(1);
	scene.aspectratio.y = aspectRatio(2);
	scene.aspectratio.z = aspectRatio(3);

	%-camera eye-%
	scene.camera.eye.x = cameraEye(1) / normFac;
    scene.camera.eye.y = cameraEye(2) / normFac;
    scene.camera.eye.z = cameraEye(3) / normFac;

    %-camera up-%
    scene.camera.up.x = cameraUpVector(1);
    scene.camera.up.y = cameraUpVector(2);
    scene.camera.up.z = cameraUpVector(3);

    %-camera projection-%
    % scene.camera.projection.type = axisData.Projection;

	%-scene axis configuration-%
	scene.xaxis.range = axisData.XLim;
	scene.yaxis.range = axisData.YLim;
	scene.zaxis.range = axisData.ZLim;

	scene.xaxis.zeroline = false;
	scene.yaxis.zeroline = false;
	scene.zaxis.zeroline = false;

	scene.xaxis.showline = true;
	scene.yaxis.showline = true;
	scene.zaxis.showline = true;

	scene.xaxis.ticklabelposition = 'outside';
	scene.yaxis.ticklabelposition = 'outside';
	scene.zaxis.ticklabelposition = 'outside';

	scene.xaxis.title = axisData.XLabel.String;
	scene.yaxis.title = axisData.YLabel.String;
	scene.zaxis.title = axisData.ZLabel.String;

	%-tick labels-%
	scene.xaxis.tickvals = axisData.XTick;
	scene.xaxis.ticktext = axisData.XTickLabel;
	scene.yaxis.tickvals = axisData.YTick;
	scene.yaxis.ticktext = axisData.YTickLabel;
	scene.zaxis.tickvals = axisData.ZTick;
	scene.zaxis.ticktext = axisData.ZTickLabel;

	scene.xaxis.tickcolor = 'rgba(0,0,0,1)';
	scene.yaxis.tickcolor = 'rgba(0,0,0,1)';
	scene.zaxis.tickcolor = 'rgba(0,0,0,1)';
	scene.xaxis.tickfont.size = axisData.FontSize;
	scene.yaxis.tickfont.size = axisData.FontSize;
	scene.zaxis.tickfont.size = axisData.FontSize;
	scene.xaxis.tickfont.family = matlab2plotlyfont(axisData.FontName);
	scene.yaxis.tickfont.family = matlab2plotlyfont(axisData.FontName);
	scene.zaxis.tickfont.family = matlab2plotlyfont(axisData.FontName);

	%-grid-%
	if strcmp(axisData.XGrid, 'off'), scene.xaxis.showgrid = false; end
	if strcmp(axisData.YGrid, 'off'), scene.yaxis.showgrid = false; end
	if strcmp(axisData.ZGrid, 'off'), scene.zaxis.showgrid = false; end

	%-SET SCENE TO LAYOUT-%
	obj.layout = setfield(obj.layout, sprintf('scene%d', xSource), scene);
end

function fillColor = getFaceColor(plotData, axisData)

	%-initializations-%
	faceColor = plotData.FaceColor;
	cData = plotData.CData;
	cLim = axisData.CLim;
	colorMap = axisData.Colormap;

	%-get face color depending of faceColor attribute
	if isnumeric(faceColor)
		numColor = 255 * faceColor;
		fillColor = sprintf('rgb(%f,%f,%f)', numColor);

	elseif strcmpi(faceColor, 'flat')
		fillColor = getStringColor(cData, colorMap, cLim);

    elseif strcmpi(faceColor, 'interp')
    	if size(cData, 1) ~= 1
    		for n = 1:size(cData, 2)
				fillColor{n} = getStringColor(mean(cData(:, n)), colorMap, cLim);
    		end
    	else
    		% TODO
    	end
	end
end

function stringColor = getStringColor(cData, colorMap, cLim)
	nColors = size(colorMap, 1);
	cIndex = max( min( cData, cLim(2) ), cLim(1) );
    scaleColor = (cIndex - cLim(1)) / diff(cLim);
    cIndex = 1 + floor(scaleColor*(nColors-1));
    numColor =  255 * colorMap(cIndex, :);
    stringColor = sprintf('rgb(%f,%f,%f)', numColor);
end

author: plotly


text:
function Fcell = get_signals(ops, iplane)

try
   load(sprintf('%s/F_%s_%s_plane%d_Nk%d.mat',...
       ops.ResultsSavePath,ops.mouse_name, ops.date, iplane, ops.Nk))
catch
   error('Could not find cell detection file \n')
end

Nk = numel(stat);
Nkpar = ops.Nk;

%% get signals
[Ly Lx] = size(ops.mimg);

nimgbatch = 2000;

ix = 0;
fclose all;
fid = fopen(ops.RegFile, 'r');

tic
F = zeros(Nk, sum(ops.Nframes), 'single');
while 1
    data = fread(fid,  Ly*Lx*nimgbatch, '*int16');
    if isempty(data)
       break;
    end
    data = reshape(data, Ly, Lx, []);
    data = data(ops.yrange, ops.xrange, :);
    data = single(data);
    NT= size(data,3);

    data = single(reshape(data, [], NT));

    for k = 1:Nk
       ipix = stat(k).ipix;
       if ~isempty(ipix)
%            F(k,ix + (1:NT)) = stat(k).lambda' * data(ipix,:);
            F(k,ix + (1:NT)) = mean(data(ipix,:), 1);
       end
    end

    ix = ix + NT;
    if rem(ix, 3*NT)==0
        fprintf('Frame %d done in time %2.2f \n', ix, toc)
    end
end
fclose(fid);
% F = F(:, 1:ix);

csumNframes = [0 cumsum(ops.Nframes)];
Fcell = cell(1, length(ops.Nframes));
for i = 1:length(ops.Nframes)
    Fcell{i} = F(:, csumNframes(i) + (1:ops.Nframes(i)));
end

save(sprintf('%s/F_%s_%s_plane%d_Nk%d.mat', ops.ResultsSavePath, ...
    ops.mouse_name, ops.date, iplane, ops.Nk),  'ops', 'res', 'stat', 'stat0', 'res0', 'Fcell', 'clustrules')

author: cortex-lab


text:
function disp(f)
%DISP  Display a BALLFUNV to the command line.

% Copyright 2019 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

loose = strcmp(get(0,'FormatSpacing'),'loose');

% Compact version:
if ( isempty(f) )
    fprintf('empty ballfunv\n\n')
    return
end

disp('   ballfunv object containing')
if ( loose )
    fprintf('\n');
end

F = f.comp;
for j = 1:3
    disp(F{j});
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun