text:
classdef AgastFeatureDetector < handle
    %AGASTFEATUREDETECTOR  Wrapping class for feature detection using the AGAST method
    %
    % Detects corners using the AGAST algorithm by [mair2010].
    %
    % ## References
    % [mair2010]:
    % > E. Mair, G. D. Hager, D. Burschka, M. Suppa, and G. Hirzinger.
    % > "Adaptive and generic corner detection based on the accelerated segment test".
    % > In "European Conference on Computer Vision (ECCV'10)", Sept 2010.
    % > [Link](http://www6.in.tum.de/Main/ResearchAgast).
    %
    % See also: cv.AGAST, cv.FAST, cv.FeatureDetector
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % Threshold on difference between intensity of the central pixel
        % and pixels of a circle around this pixel.
        %
        % Default is 10
        Threshold
        % If true, non-maximum suppression is applied to detected corners
        % (keypoints).
        %
        % Default is true.
        NonmaxSuppression
        % One of the four neighborhoods as defined in the paper:
        %
        % * **AGAST_5_8**
        % * **AGAST_7_12d**
        % * **AGAST_7_12s**
        % * **OAST_9_16** (default)
        Type
    end

    methods
        function this = AgastFeatureDetector(varargin)
            %AGASTFEATUREDETECTOR  Constructor
            %
            %     obj = cv.AgastFeatureDetector()
            %     obj = cv.AgastFeatureDetector(..., 'OptionName',optionValue, ...)
            %
            % ## Options
            % * __Threshold__ See cv.AgastFeatureDetector.Threshold, default 10
            % * __NonmaxSuppression__ See
            %   cv.AgastFeatureDetector.NonmaxSuppression, default true
            % * __Type__ See cv.AgastFeatureDetector.Type, default `OAST_9_16`
            %
            % See also: cv.AgastFeatureDetector.detect
            %
            this.id = AgastFeatureDetector_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.AgastFeatureDetector
            %
            if isempty(this.id), return; end
            AgastFeatureDetector_(this.id, 'delete');
        end

        function typename = typeid(this)
            %TYPEID  Name of the C++ type (RTTI)
            %
            %     typename = obj.typeid()
            %
            % ## Output
            % * __typename__ Name of C++ type
            %
            typename = AgastFeatureDetector_(this.id, 'typeid');
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.AgastFeatureDetector.empty,
            %  cv.AgastFeatureDetector.load
            %
            AgastFeatureDetector_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Checks if detector object is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the detector object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.AgastFeatureDetector.clear,
            %  cv.AgastFeatureDetector.load
            %
            b = AgastFeatureDetector_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.AgastFeatureDetector.load
            %
            AgastFeatureDetector_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.AgastFeatureDetector.save
            %
            AgastFeatureDetector_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.AgastFeatureDetector.save, cv.AgastFeatureDetector.load
            %
            name = AgastFeatureDetector_(this.id, 'getDefaultName');
        end
    end

    %% Features2D: FeatureDetector
    methods
        function keypoints = detect(this, img, varargin)
            %DETECT  Detects keypoints in an image or image set
            %
            %     keypoints = obj.detect(img)
            %     keypoints = obj.detect(imgs)
            %     [...] = obj.detect(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __img__ Image (first variant), 8-bit grayscale image where
            %   keypoints (corners) are detected.
            % * __imgs__ Image set (second variant), cell array of images.
            %
            % ## Output
            % * __keypoints__ The detected keypoints. In the first variant, a
            %   1-by-N structure array. In the second variant of the method,
            %   `keypoints{i}` is a set of keypoints detected in `imgs{i}`.
            %
            % ## Options
            % * __Mask__ A mask specifying where to look for keypoints
            %   (optional). It must be a logical or 8-bit integer matrix with
            %   non-zero values in the region of interest. In the second
            %   variant, it is a cell-array of masks for each input image,
            %   `masks{i}` is a mask for `imgs{i}`. Not set by default.
            %
            % See also: cv.AgastFeatureDetector.AgastFeatureDetector
            %
            keypoints = AgastFeatureDetector_(this.id, 'detect', img, varargin{:});
        end
    end

    %% Getters/Setters
    methods
        function value = get.Threshold(this)
            value = AgastFeatureDetector_(this.id, 'get', 'Threshold');
        end
        function set.Threshold(this, value)
            AgastFeatureDetector_(this.id, 'set', 'Threshold', value);
        end

        function value = get.NonmaxSuppression(this)
            value = AgastFeatureDetector_(this.id, 'get', 'NonmaxSuppression');
        end
        function set.NonmaxSuppression(this, value)
            AgastFeatureDetector_(this.id, 'set', 'NonmaxSuppression', value);
        end

        function value = get.Type(this)
            value = AgastFeatureDetector_(this.id, 'get', 'Type');
        end
        function set.Type(this, value)
            AgastFeatureDetector_(this.id, 'set', 'Type', value);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [samples,gammas,Lambda,joints2D,hists] = DGAnyMarginal(pmfs,Sigma,supports,Nsamples)

% [samples,gammas,Lambda,joints2D,hists] = DGAnyMarginal(pmfs,Sigma,supports,Nsamples)
%   Generates samples from a Multivariate Discretized Gaussian with specified marginals
%   and covariance. Also returns parameters of the fitted DG.
%
%   Inputs:
%     pmfs: the probability mass functions of the marginal distribution of the
%       input-random variables. Must be a cell-array with n elements, each of
%       which is a vector which sums to one
%     Sigma: The covariance matrix of the input-random variable. The function
%       does not check for admissability, i.e. results might be wrong if there
%       exists no random variable which has the specified marginals and
%       covariance.
%     supports: The support of each dimension of the input random variable.
%       Must be a cell-array with n elements, each of whcih is a vector with
%       increasing entries giving the possible values of each random variable,
%       e.g. if the first dimension of the rv is 1 with probability .2, 3 with
%       prob .8, then pmfs{1}=[.2,.8], supports{1}=[1,3]; If empty support is
%       specified, then each is taken to be [0:numel(pdfs{k}-1];
%     Nsamples: The number of samples to be generated.
%
%   Outputs:
%     samples: a matrix of size Nsamples by n, where each row is a sample from
%       the DG.
%     gammas: the discretization thresholds, as described in the paper. When
%       sampling. The k-th dimension of the output random variable is f if e.g.
%       supports{k}(1)=f and gammas{k}(f) <= U(k) <= gammas{k}(f+1)
%     Lambda: the covariance matrix of the latent Gaussian random variable U
%     joints2D: An n by n cell array, where each entry contains the 2
%       dimensional joint distribution of  a pair of dimensions of the DG.
%     hists: the empirical marginals of the samples returned in "samples"
%
%   Usage:
%
% Code from the paper: 'Generating spike-trains with specified
% correlations', Macke et al., submitted to Neural Computation
%
% www.kyb.mpg.de/bethgegroup/code/efficientsampling

[gammas,Lambda,joints2D]=FindDGAnyMarginal(pmfs,Sigma,supports);

[samples,hists]=SampleDGAnyMarginal(gammas,Lambda,supports,Nsamples);

samples = samples';

author: Sable

repo: mcbench-benchmarks


text:
function chrpak_test1225 ( )

%*****************************************************************************80
%
%% TEST1225 tests S_S_SUBANAGRAM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2007
%
%  Author:
%
%    John Burkardt
%
  test_num = 4;

  s1_test(1,1:13) = 'Get a priest!';
  s1_test(2,1:13) = 'Get a priest!';
  s1_test(3,1:13) = 'Get a priest!';
  s1_test(4,1:13) = 'Get a priest!';

  s2_test(1,1:6) = 'stripe';
  s2_test(2,1:6) = 'pastor';
  s2_test(3,1:6) = 'a sip ';
  s2_test(4,1:6) = 'tag!  ';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST1225\n' );
  fprintf ( 1, '  S_S_SUBANAGRAM is TRUE if S2 is a "subanagram"\n' );
  fprintf ( 1, '  of S1.\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    s1 = s1_test(test,1:13);
    s2 = s2_test(test,1:6);
    s2_length = s_len_trim ( s2 );

    value = s_s_subanagram ( s1, s2(1:s2_length) );

    fprintf ( 1, '  "%s"  "%s"  %d\n', s1, s2, value );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%*****************************************************************
%% SDPT3data_SEDUMIdata: convert SQLP data in SDPT3 format to
%%                       SeDuMi format
%%
%% [At,b,c,K] = SDPT3data_SEDUMIdata(blk,AAt,CC,bb);
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

  function [At,b,c,K] = SDPT3data_SEDUMIdata(blk,AAt,CC,bb);

  c = []; At = [];
  b = bb;
  mm = length(bb);
%%
  if (~iscell(CC))
     Ctmp = CC; clear CC; CC{1} = Ctmp;
  end
%%
%% extract unrestricted blk
%%
  for p = 1:size(blk,1)
     pblk = blk(p,:);
     if (p==1); K.f = []; end
     if strcmp(pblk{1},'u')
        K.f = [K.f, pblk{2}];
        At  = [At; AAt{p}];
        c   = [c; CC{p}];
     end
  end
  K.f = sum(K.f);
%%
%% extract linear blk
%%
  for p = 1:size(blk,1)
     pblk = blk(p,:);
     if (p==1); K.l = []; end
     if strcmp(pblk{1},'l')
        K.l = [K.l, pblk{2}];
        At  = [At; AAt{p,1}];
        c   = [c; CC{p,1}];
     end
  end
  K.l = sum(K.l);
%%
%% extract second order cone blk
%%
  for p = 1:size(blk,1)
     pblk = blk(p,:);
     if (p==1); K.q = []; end
     if strcmp(pblk{1},'q')
        K.q = [K.q, pblk{2}];
        At  = [At; AAt{p,1}];
        c   = [c; CC{p,1}];
     end
  end
%%
%% extract rotated cone blk
%%
  for p = 1:size(blk,1)
     pblk = blk(p,:);
     if (p==1); K.r = []; end
     if strcmp(pblk{1},'r')
        K.r = [K.r, pblk{2}];
        At  = [At; AAt{p,1}];
        c   = [c; CC{p,1}];
     end
  end
%%
%% extract semidefinite cone blk
%%
  for p = 1:size(blk,1)
     if (p==1); K.s = []; end
     pblk = blk(p,:);
     if strcmp(pblk{1},'s')
        K.s = [K.s, pblk{2}];
        ss = [0,cumsum(pblk{2})];
        idxstart = [0,cumsum(pblk{2}.*pblk{2})];
        numblk = length(pblk{2});
        nnzA = nnz(AAt{p,1});
        II = zeros(2*nnzA,1);
        JJ = zeros(2*nnzA,1);
        VV = zeros(2*nnzA,1);
        m2 = size(AAt{p,1},2);
        if (length(pblk) > 2)
           rr = [0, cumsum(pblk{3})];
	   dd = AAt{p,3};
           idxD = [0; find(diff(dd(:,1))); size(dd,1)];
        end
        count = 0;
        for k = 1:mm
   	   if (k<= m2);
              Ak = smat(pblk,AAt{p,1}(:,k),1);
           else
              idx = [rr(k)+1 : rr(k+1)];
              Vk  = AAt{p,2}(:,idx);
              len = pblk{3}(k);
              if (size(dd,2) == 4)
                 idx2 = [idxD(k)+1:idxD(k+1)];
                 Dk = spconvert([dd(idx2,2:4); len,len,0]);
              elseif (size(dd,2) == 1);
 	         Dk = spdiags(dd(idx),0,len,len);
              end
              Ak = Vk*Dk*Vk';
           end
           for tt = 1:numblk
              if (numblk > 1)
                 idx = [ss(tt)+1: ss(tt+1)];
                 Aksub = full(Ak(idx,idx));
              else
                 Aksub = Ak;
              end
              tmp = Aksub(:);
              nzidx = find(tmp);
              len = length(nzidx);
              II(count+[1:len],1) = idxstart(tt)+nzidx;
              JJ(count+[1:len],1) = k*ones(length(nzidx),1);
              VV(count+[1:len],1) = tmp(nzidx);
              count = count + len;
           end
        end
        II = II(1:count);
        JJ = JJ(1:count);
        VV = VV(1:count);
        At = [At; spconvert([II,JJ,VV; sum(pblk{2}.*pblk{2}), mm, 0])];
	    Cp = CC{p};
        ctmp = [];
        for tt = 1:numblk
   	   if (numblk > 1)
              idx = [ss(tt)+1: ss(tt+1)];
              Csub = full(Cp(idx,idx));
           else
              Csub = Cp;
           end
           ctmp = [ctmp; Csub(:)];
        end
        c = [c; ctmp];
     end
  end
%%**********************************************************

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function yBin = BF_Binarize(y,binarizeHow)
% BF_Binarize    Converts an input vector into a binarized version

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function 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/>.
% ------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
% Check inputs, set defaults:
%-------------------------------------------------------------------------------
if nargin < 2 || isempty(binarizeHow)
    binarizeHow = 'diff';
end

%-------------------------------------------------------------------------------
% Do the binary transformation:
%-------------------------------------------------------------------------------

switch binarizeHow
    case 'diff'
        % Binary signal: 1 for stepwise increases, 0 for stepwise decreases
        yBin = stepBinary(diff(y));

    case 'mean'
        % Binary signal: 1 for above mean, 0 for below mean
        yBin = stepBinary(y - mean(y));

    case 'median'
        % Binary signal: 1 for above median, 0 for below median
        yBin = stepBinary(y - median(y));

    case 'iqr'
        % Binary signal: 1 if inside interquartile range, 0 otherwise
        iqr = quantile(y,[0.25, 0.75]);
        iniqr = (y > iqr(1) & y <= iqr(2));
        yBin = zeros(length(y),1);
        yBin(iniqr) = 1;

    otherwise
        error('Unknown binary transformation setting ''%s''',binarizeHow)
end

%-------------------------------------------------------------------------------
function Y = stepBinary(X)
    % Transform real values to 0 if <=0 and 1 if >0:
    Y = zeros(size(X),'like',X);
    Y(X > 0) = 1;
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [p,infeasible] = detectRedundantInfeasibleSOCPRows(p)
infeasible = 0;
if any(p.K.q)
    newEqualities = [];
    top = startofSOCPCone(p.K);
    for j = 1:length(p.K.q)
        % Looking for norm(x) <= k and k=0 or k<0
        m = p.K.q(j);
        M = p.F_struc(top:top+m-1,:);
        fixed = ~any(M(1,2:end));
        if fixed
            if M(1) <= -p.options.bnb.feastol
                infeasible = 1;
                return
            elseif M(1) <= 0
                % Not numerically infeasible,
                % but at least says x = 0
                equalityRows = M(2:m,:);
                newEqualities = [newEqualities;equalityRows];
                p.F_struc(top:top+m-1,:) = [];
                p.K.q(j) = 0;
            end
            top = top + p.K.q(j);
        end
    end
    p = addEquality(p,newEqualities);
    p.K.q(p.K.q == 0) = [];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function dataset = voc2007_trainval(dataset, usage, use_flip)
% Pascal voc 2007 trainval set
% set opts.imdb_train opts.roidb_train
% or set opts.imdb_test opts.roidb_train

% change to point to your devkit install
devkit                      = voc2007_devkit();

switch usage
    case {'train'}
        dataset.imdb_train    = {  imdb_from_voc(devkit, 'trainval', '2007', use_flip) };
        dataset.roidb_train   = cellfun(@(x) x.roidb_func(x), dataset.imdb_train, 'UniformOutput', false);
    case {'test'}
        dataset.imdb_test     = imdb_from_voc(devkit, 'trainval', '2007', use_flip) ;
        dataset.roidb_test    = dataset.imdb_test.roidb_func(dataset.imdb_test);
    otherwise
        error('usage = ''train'' or ''test''');
end

end

author: ShaoqingRen


text:
function [x]=als_rake_solve(a,y,tol,varargin)
%ALS-MINRES method for the solution of linear systems in QTT-Tucker format
%   [X]=ALS_RAKE_SOLVE(A,Y,TOL,VARARGIN) Attempts to solve the linear
%   system A*X = Y with accuracy EPS using the AMR 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 [tensor of all ones]
%       o nswp - maximal number of AMR sweeps [10]
%       o verb - verbosity level, 0-silent, 1-sweep info, 2-block info [1]
%       o kickrank - rank-increasing parameter [5]
%       o max_full_size - maximal size of the local matrix to full solver
%       [2500]
%       o local_prec: Local preconditioner, '' - none
%         'cjacobi' - block-Jacobi, diagonal over rank indices ['']
%       o local_iters - number of local gmres restarts [2]
%       o local_restart - dimension of local gmres [40]
%       o trunc_norm - truncation and stopping tolerance: 'resid' - the
%       residual norm is used, 'fro' - Frobenius (L2)  ['resid']
%       o resid_damp - gap between the local solver and truncation. Larger
%       values help to filter some noise, but require more local
%       iterations   [2]
%
%
% 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
%---------------------------

Au = [];

% Inner parameters
max_full_size=2500;

resid_damp = 2; % Truncation error to true residual treshold
als_tol_high = 10;
als_tol_low = 4;
als_iters=2;

nswp=10;
local_restart=40;
local_iters=2;

local_prec = '';
% local_prec_char = 0;
% local_prec = 'jacobi';

rmax=Inf;
trunc_norm = 'residual';
% trunc_norm_char = 1;
% trunc_norm = 'fro';

% local_solver = 'gmres';
% local_solver = 'pcg';

verb=1;
kickrank = 5;
x=[];

nswp_f = 1;
nswp_c = 1;

for i=1:2:length(varargin)-1
    switch lower(varargin{i})
        case 'nswp'
            nswp=varargin{i+1};
        case 'rmax'
            rmax=varargin{i+1};
        case 'x0'
            x=varargin{i+1};
        case 'verb'
            verb=varargin{i+1};
        case 'local_prec'
            local_prec=varargin{i+1};
        case 'local_restart'
            local_restart=varargin{i+1};
        case 'local_iters'
            local_iters=varargin{i+1};
%         case 'local_solver'
%             local_solver=varargin{i+1};
        case 'kickrank'
            kickrank=varargin{i+1};
        case 'als_tol_high'
            als_tol_high=varargin{i+1};
        case 'als_tol_low'
            als_tol_low=varargin{i+1};
        case 'als_iters'
            als_iters=varargin{i+1};
        case  'max_full_size'
            max_full_size=varargin{i+1};
%         case 'step_dpow'
%             step_dpow=varargin{i+1};
%         case 'min_dpow'
%             min_dpow=varargin{i+1};
        case 'resid_damp'
            resid_damp = varargin{i+1};
        case 'trunc_norm'
            trunc_norm = varargin{i+1};
%         case 'bot_conv'
%             bot_conv=varargin{i+1};
%         case 'top_conv'
%             top_conv=varargin{i+1};

        otherwise
            error('Unrecognized option: %s\n',varargin{i});
    end
end

tol2 = tol;

d = y.core.d;
yc = core2cell(y.core);
ryc = y.core.r;
n = cell(d,1);
L = zeros(d,1);
yf = cell(d,1);
ryf = cell(d,1);
for i=1:d
    n{i} = y.tuck{i}.n;
    L(i) = y.tuck{i}.d;
    ryf{i} = y.tuck{i}.r;
    ryf{i}(L(i)+2) = ryc(i+1);
    yf{i} = cell(L(i)+1,1);
    yf{i}(1:L(i)) = core2cell(y.tuck{i});
    yf{i}{L(i)+1} = yc{i};
end;

ac = core2cell(a.core);
rac = a.core.r;
af = cell(d,1);
raf = cell(d,1);
for i=1:d
    raf{i} = a.tuck{i}.r;
    raf{i}(L(i)+2) = rac(i+1);
    af{i} = cell(L(i)+1,1);
    af{i}(1:L(i)) = core2cell(a.tuck{i});
    af{i}{L(i)+1} = ac{i};
end;

xf = cell(d,1);
rxf = cell(d,1);
if (isempty(x))
    x = qtt_tucker;
    xc = tt_ones(1, d);
    x.core = xc;
    rxc = xc.r;
    xc = core2cell(xc);
    x.tuck = cell(d,1);
    for i=1:d
        rxf{i} = ones(L(i)+2,1);
        n{i}(L(i)+1) = 1;
        xf{i} = cell(L(i)+1,1);
        x.tuck{i} = tt_ones(n{i}(1:L(i)));
        xf{i}(1:L(i)) = core2cell(x.tuck{i});
        xf{i}{L(i)+1}=xc{i};
    end;
else
    xc = core2cell(x.core);
    rxc = x.core.r;
    for i=1:d
        rxf{i} = x.tuck{i}.r;
        rxf{i}(L(i)+2)=rxc(i+1);
        n{i}(L(i)+1) = rxc(i);
        xf{i} = cell(L(i)+1,1);
        xf{i}(1:L(i)) = core2cell(x.tuck{i});
        xf{i}{L(i)+1}=xc{i};
    end;
end;

max_dx = 0;
max_res = 0;

phiaf = cell(d,1);
phiyf = cell(d,1);
for i=1:d
    phiaf{i} = cell(L(i)+2,1);
    phiaf{i}{1}=1;
    phiyf{i} = cell(L(i)+2,1);
    phiyf{i}{1}=1;
end;
phiac = cell(d+1,1); phiac{1}=1; phiac{d+1}=1;
phiyc = cell(d+1,1); phiyc{1}=1; phiyc{d+1}=1;

acp = cell(d,1);
ycp = cell(d,1);

max_res_prev = Inf;
max_dx_prev = Inf;
max_frank = 0;
regurg_cnt = 0;
last_sweep = false;

for swp=1:nswp

    % Orthog + phi, Factors
    for i=1:d
        n{i}(L(i)+1)=rxc(i);
        rxf{i}(L(i)+2)=rxc(i+1);
        % Permute core blocks: rc1 is mode index now
        xf{i}{L(i)+1} = permute(xf{i}{L(i)+1}, [2,1,3]); % rxt,rxc1,rxc2
        for j=1:L(i)
            cr = xf{i}{j};
            cr = reshape(cr, rxf{i}(j)*n{i}(j), rxf{i}(j+1));
            [cr,rv]=qr(cr,0);
            cr2 = xf{i}{j+1};
            cr2 = reshape(cr2, rxf{i}(j+1), n{i}(j+1)*rxf{i}(j+2));
            cr2 = rv*cr2;
            rxf{i}(j+1) = size(cr, 2);
            cr = reshape(cr, rxf{i}(j), n{i}(j), rxf{i}(j+1));
            xf{i}{j} = cr;
            xf{i}{j+1} = reshape(cr2, rxf{i}(j+1), n{i}(j+1), rxf{i}(j+2));

            phiaf{i}{j+1} = compute_next_Phi(phiaf{i}{j}, cr, af{i}{j}, cr, 'lr');
            phiyf{i}{j+1} = compute_next_Phi(phiyf{i}{j}, cr, [], yf{i}{j}, 'lr');
        end;
        % core-projected matrix block
        acp{i} = permute(phiaf{i}{L(i)+1}, [2, 1, 3]);
        acp{i} = reshape(acp{i}, raf{i}(L(i)+1), rxf{i}(L(i)+1)*rxf{i}(L(i)+1));
        acp{i} = reshape(permute(af{i}{L(i)+1}, [1,3,2]), rac(i)*rac(i+1), raf{i}(L(i)+1)) * acp{i};
        acp{i} = reshape(acp{i}, rac(i), rac(i+1), rxf{i}(L(i)+1), rxf{i}(L(i)+1));
        acp{i} = permute(acp{i}, [1, 3, 4, 2]);
        % core-projected rhs block
        ycp{i} = phiyf{i}{L(i)+1}.';
        ycp{i} = reshape(permute(yf{i}{L(i)+1}, [1,3,2]), ryc(i)*ryc(i+1), ryf{i}(L(i)+1)) * ycp{i};
        ycp{i} = reshape(ycp{i}, ryc(i), ryc(i+1), rxf{i}(L(i)+1));
        ycp{i} = permute(ycp{i}, [1, 3, 2]);
        % return xf{i}{L{i}+1} to the core state
        xf{i}{L(i)+1} = permute(xf{i}{L(i)+1}, [2,1,3]); % rxc1,rxt,rxc2
    end;

    % ALS over core
    nc = zeros(d,1);
    for i=1:d
        nc(i)=rxf{i}(L(i)+1);
    end;
    for k=1:nswp_c
%     max_res_c = 0;
    % Orthog
    for i=1:d-1
        cr = xf{i}{L(i)+1};
        cr = reshape(cr, rxc(i)*nc(i), rxc(i+1));
        [cr,rv]=qr(cr,0);
        cr2 = xf{i+1}{L(i+1)+1};
        cr2 = reshape(cr2, rxc(i+1), nc(i+1)*rxc(i+2));
        cr2 = rv*cr2;
        rxc(i+1) = size(cr,2);
        % Update Factor-sizes too
%         n{i+1}(L(i+1)+1)=rxc(i+1);
%         rxf{i}(L(i)+2)=rxc(i+1);
        cr = reshape(cr, rxc(i), nc(i), rxc(i+1));
        xf{i}{L(i)+1} = cr;
        xf{i+1}{L(i+1)+1} = reshape(cr2, rxc(i+1), nc(i+1), rxc(i+2));

        phiac{i+1} = compute_next_Phi(phiac{i}, cr, acp{i}, cr, 'lr');
        phiyc{i+1} = compute_next_Phi(phiyc{i}, cr, [], ycp{i}, 'lr');
    end;
    % ALS
    for i=d:-1:1
        Phi1 = phiac{i};
        A1 = acp{i};
        Phi2 = phiac{i+1};

        sol_prev = reshape(xf{i}{L(i)+1}, rxc(i)*nc(i)*rxc(i+1), 1);

        dir = -1;
        if (i==1); dir=0; end;
        local_kickrank = kickrank;
        if (last_sweep); local_kickrank=0; end;
        if (verb>1); fprintf('\t core %d\n', i); end;
        [u,v,max_res,max_dx,flg]=local_solve(Phi1,A1,Phi2, phiyc{i},ycp{i},phiyc{i+1}, ...
            tol2/sqrt(sum(L+1))/resid_damp, resid_damp, trunc_norm, sol_prev, ...
            local_prec, local_restart, local_iters, max_full_size, max_res, max_dx, ...
            dir, rmax, local_kickrank, verb);

        if (flg>0); fprintf('-warn- local_solve did not converge at cb {%d}\n', i); end;

        if (i>1)
            cr2 = xf{i-1}{L(i-1)+1};
            cr2 = reshape(cr2, rxc(i-1)*nc(i-1), rxc(i));
            cr2 = cr2*u;
            rxc(i) = size(v,2);
            % Update Factor-sizes
%             n{i}(L(i)+1)=rxc(i);
%             rxf{i-1}(L(i-1)+2)=rxc(i);
            v = reshape(v.', rxc(i), nc(i), rxc(i+1));
            xf{i}{L(i)+1} = v;
            xf{i-1}{L(i-1)+1} = reshape(cr2, rxc(i-1), nc(i-1), rxc(i));

            phiac{i} = compute_next_Phi(phiac{i+1}, v, acp{i}, v, 'rl');
            phiyc{i} = compute_next_Phi(phiyc{i+1}, v, [], ycp{i}, 'rl');
        else
            v = u*(v.');
            xf{i}{L(i)+1} = reshape(v, rxc(i), nc(i), rxc(i+1));
        end;
    end;
    end;
%     max_res = max(max_res, max_res_c);

    % Optimization over factors
    for i=1:d
        % ALS L->1, the orth. holds
        % Copy the extended factor data to the current array
        curx = xf{i};
        % Last mode size <- largest core rank
        if (rxc(i)>rxc(i+1))
            % Reshape the core block to be the factor one
            curx{L(i)+1} = permute(curx{L(i)+1}, [2, 1, 3]); % rtx, rxc1, rxc2

            cura = af{i};
            % last block has to be convolved with phiac-left
            cura{L(i)+1} = permute(phiac{i}, [1,3,2]);
            cura{L(i)+1} = reshape(cura{L(i)+1}, rxc(i)*rxc(i), rac(i));
            cura{L(i)+1} = cura{L(i)+1}*reshape(af{i}{L(i)+1}, rac(i), raf{i}(L(i)+1)*rac(i+1));
            cura{L(i)+1} = reshape(cura{L(i)+1}, rxc(i), rxc(i), raf{i}(L(i)+1), rac(i+1));
            cura{L(i)+1} = permute(cura{L(i)+1}, [3, 1,2, 4]); % tucker rank is now the left

            cury = yf{i};
            cury{L(i)+1} = phiyc{i}*reshape(yf{i}{L(i)+1}, ryc(i), ryf{i}(L(i)+1)*ryc(i+1));
            cury{L(i)+1} = reshape(cury{L(i)+1}, rxc(i), ryf{i}(L(i)+1), ryc(i+1));
            cury{L(i)+1} = permute(cury{L(i)+1}, [2, 1, 3]); % tucker rank is now the left

            % The L+2-th phi comes from the core
            phiaf{i}{L(i)+2} = phiac{i+1}; % sizes rxc2,ra2,rxc2. rxc2 - our last "rank" index
            phiyf{i}{L(i)+2} = phiyc{i+1};

            n{i}(L(i)+1) = rxc(i);
            rxf{i}(L(i)+2) = rxc(i+1);
        else
            % Reshape the core block to be the factor one
            curx{L(i)+1} = permute(curx{L(i)+1}, [2, 3, 1]); %  rtx, rxc2, rxc1

            cura = af{i};
            % last block has to be convolved with phiac-left
            cura{L(i)+1} = permute(phiac{i+1}, [2, 1,3]);
            cura{L(i)+1} = reshape(cura{L(i)+1}, rac(i+1), rxc(i+1)*rxc(i+1));
            cura{L(i)+1} = reshape(af{i}{L(i)+1}, rac(i)*raf{i}(L(i)+1), rac(i+1)) * cura{L(i)+1};
            cura{L(i)+1} = reshape(cura{L(i)+1}, rac(i), raf{i}(L(i)+1), rxc(i+1), rxc(i+1));
            cura{L(i)+1} = permute(cura{L(i)+1}, [2, 3,4, 1]); % tucker rank is now the left

            cury = yf{i};
            cury{L(i)+1} = reshape(yf{i}{L(i)+1}, ryc(i)*ryf{i}(L(i)+1), ryc(i+1)) * (phiyc{i+1}.');
            cury{L(i)+1} = reshape(cury{L(i)+1}, ryc(i), ryf{i}(L(i)+1), rxc(i+1));
            cury{L(i)+1} = permute(cury{L(i)+1}, [2, 3, 1]); % tucker rank is now the left

            % The L+2-th phi comes from the core
            phiaf{i}{L(i)+2} = phiac{i}; % sizes rxc2,ra2,rxc2. rxc2 - our last "rank" index
            phiyf{i}{L(i)+2} = phiyc{i};

            n{i}(L(i)+1) = rxc(i+1);
            rxf{i}(L(i)+2) = rxc(i);
        end;

        for k=1:nswp_f
%         max_res_f = 0;
        % First, orthogonality L->1
        for j=(L(i)+1):-1:2
            cr = curx{j};
            cr = reshape(cr, rxf{i}(j), n{i}(j)*rxf{i}(j+1));
            [cr,rv]=qr(cr.',0);
            cr2 = curx{j-1};
            cr2 = reshape(cr2, rxf{i}(j-1)*n{i}(j-1), rxf{i}(j));
            cr2 = cr2*(rv.');
            rxf{i}(j) = size(cr, 2);
            cr = reshape(cr.', rxf{i}(j), n{i}(j), rxf{i}(j+1));
            curx{j} = cr;
            curx{j-1} = reshape(cr2, rxf{i}(j-1), n{i}(j-1), rxf{i}(j));

            phiaf{i}{j} = compute_next_Phi(phiaf{i}{j+1}, cr, cura{j}, cr, 'rl');
            phiyf{i}{j} = compute_next_Phi(phiyf{i}{j+1}, cr, [], cury{j}, 'rl');
        end;

        % Now, finally, the optimization itself
%         for j=(L(i)+1):-1:1
        for j=1:(L(i)+1)
            Phi1 = phiaf{i}{j};
            A1 = cura{j};
            Phi2 = phiaf{i}{j+1};

            sol_prev = reshape(curx{j}, rxf{i}(j)*n{i}(j)*rxf{i}(j+1), 1);

            dir = 1;
            if (j==(L(i)+1)); dir=0; end;
            local_kickrank = kickrank;
            if (last_sweep); local_kickrank=0; end;
            if (verb>1); fprintf('\t factor {%d,%d}\n', i, j); end;
            [u,v,max_res,max_dx,flg]=local_solve(Phi1,A1,Phi2, phiyf{i}{j},cury{j},phiyf{i}{j+1}, ...
                tol2/sqrt(sum(L+1))/resid_damp, resid_damp, trunc_norm, sol_prev, ...
                local_prec, local_restart, local_iters, max_full_size, max_res, max_dx, ...
                dir, rmax, local_kickrank, verb);

            if (flg>0); fprintf('-warn- local_solve did not converge at fb {%d}{%d}\n', i, j); end;

            if (j<(L(i)+1))
                cr2 = curx{j+1};
                cr2 = reshape(cr2, rxf{i}(j+1), n{i}(j+1)*rxf{i}(j+2));
                cr2 = (v.')*cr2;
                rxf{i}(j+1) = size(u,2);
                max_frank = max(max_frank, rxf{i}(j+1));
                u = reshape(u, rxf{i}(j), n{i}(j), rxf{i}(j+1));
                curx{j} = u;
                curx{j+1} = reshape(cr2, rxf{i}(j+1), n{i}(j+1), rxf{i}(j+2));

                phiaf{i}{j+1} = compute_next_Phi(phiaf{i}{j}, u, cura{j}, u, 'lr');
                phiyf{i}{j+1} = compute_next_Phi(phiyf{i}{j}, u, [], cury{j}, 'lr');
            else
                u = u*(v.');
                curx{j} = reshape(u, rxf{i}(j), n{i}(j), rxf{i}(j+1));
            end;
        end;
        end;
%         max_res = max(max_res, max_res_f);

        if (i<d)
%             % The factor is ready. Now, orthogonalize it l-to-r
%             for j=1:L(i)
%                 cr = curx{j};
%                 cr = reshape(cr, rxf{i}(j)*n{i}(j), rxf{i}(j+1));
%                 [cr,rv]=qr(cr,0);
%                 cr2 = curx{j+1};
%                 cr2 = reshape(cr2, rxf{i}(j+1), n{i}(j+1)*rxf{i}(j+2));
%                 cr2 = rv*cr2;
%                 rxf{i}(j+1) = size(cr, 2);
%                 cr = reshape(cr, rxf{i}(j), n{i}(j), rxf{i}(j+1));
%                 curx{j} = cr;
%                 curx{j+1} = reshape(cr2, rxf{i}(j+1), n{i}(j+1), rxf{i}(j+2));
%
%                 phiaf{i}{j+1} = compute_next_Phi(phiaf{i}{j}, cr, cura{j}, cr, 'lr');
%                 phiyf{i}{j+1} = compute_next_Phi(phiyf{i}{j}, cr, [], cury{j}, 'lr');
%             end;
            % First L cells are all what we need for factor
            xf{i}(1:L(i)) = curx(1:L(i));
            % The last one should go to the core. Moreover, we need xc1->xc2 orth
            if (rxc(i)>rxc(i+1))
                cr = permute(curx{L(i)+1}, [2, 1, 3]); % now rc1,rt,rc2
            else
                cr = permute(curx{L(i)+1}, [3, 1, 2]); % now rc1,rt,rc2
            end;
            cr = reshape(cr, rxc(i)*rxf{i}(L(i)+1), rxc(i+1));
            [cr,rv]=qr(cr,0);
            cr2 = xf{i+1}{L(i+1)+1};
            cr2 = reshape(cr2, rxc(i+1), rxf{i+1}(L(i+1)+1)*rxc(i+2));
            cr2 = rv*cr2;
            rxc(i+1) = size(cr,2);
            % Update Factor-sizes too
%             n{i+1}(L(i+1)+1)=rxc(i+1);
%             rxf{i}(L(i)+2)=rxc(i+1);
            cr = reshape(cr, rxc(i), rxf{i}(L(i)+1), rxc(i+1));
            xf{i}{L(i)+1} = cr;
            xf{i+1}{L(i+1)+1} = reshape(cr2, rxc(i+1), rxf{i+1}(L(i+1)+1), rxc(i+2));

            % Update acp, ycp
            % core-projected matrix block
            acp{i} = permute(phiaf{i}{L(i)+1}, [2, 1, 3]);
            acp{i} = reshape(acp{i}, raf{i}(L(i)+1), rxf{i}(L(i)+1)*rxf{i}(L(i)+1));
            acp{i} = reshape(permute(af{i}{L(i)+1}, [1,3,2]), rac(i)*rac(i+1), raf{i}(L(i)+1)) * acp{i};
            acp{i} = reshape(acp{i}, rac(i), rac(i+1), rxf{i}(L(i)+1), rxf{i}(L(i)+1));
            acp{i} = permute(acp{i}, [1, 3, 4, 2]);
            % core-projected rhs block
            ycp{i} = phiyf{i}{L(i)+1}.';
            ycp{i} = reshape(permute(yf{i}{L(i)+1}, [1,3,2]), ryc(i)*ryc(i+1), ryf{i}(L(i)+1)) * ycp{i};
            ycp{i} = reshape(ycp{i}, ryc(i), ryc(i+1), rxf{i}(L(i)+1));
            ycp{i} = permute(ycp{i}, [1, 3, 2]);

            phiac{i+1} = compute_next_Phi(phiac{i}, cr, acp{i}, cr, 'lr');
            phiyc{i+1} = compute_next_Phi(phiyc{i}, cr, [], ycp{i}, 'lr');
        else
            xf{i}(1:L(i)) = curx(1:L(i));
            if (rxc(i)>rxc(i+1))
                xf{i}{L(i)+1} = permute(curx{L(i)+1}, [2, 1, 3]); % core block
            else
                xf{i}{L(i)+1} = permute(curx{L(i)+1}, [3, 1, 2]); % core block
            end;
        end;
    end;

    % Residual check, etc

%     x_old = x;
%     for i=1:d
%         x.tuck{i} = cell2core(x.tuck{i}, xf{i}(1:L(i)));
%         xc{i} = xf{i}{L(i)+1};
%     end;
%     x.core = cell2core(x.core, xc);

    if (verb>0)
        fprintf('=als_rake_solve= sweep %d, max_dx: %3.3e, max_res: %3.3e, mrank_c: %d, mrank_f: %d\n', swp, max_dx, max_res, max(rxc), max_frank);
    end;

    if (last_sweep)
        break;
    end;

    if (kickrank<0)
        kickrank=kickrank-1;
    end;

    if (strcmp(trunc_norm, 'fro'))
        if (max_dx<tol)&&(kickrank<=-als_iters)
            last_sweep=true;
%             tol2 = tol;
        end;
        if (max_dx_prev<=tol*als_tol_high)
            if (max_dx>max_dx_prev); regurg_cnt=regurg_cnt+1; fprintf('---- Regurgitation %d\n', regurg_cnt); end;
%             kickrank = -1;
%             tol2=tol/4;
        end;
        if ((regurg_cnt>0)||(max_dx<=tol*als_tol_low))&&(kickrank>0); kickrank=-1; end;
    else
%         for i=1:d
%             x.tuck{i} = cell2core(x.tuck{i}, xf{i}(1:L(i)));
%             xc{i} = xf{i}{L(i)+1};
%         end;
%         x.core = cell2core(x.core, xc);
%         x.dphys = d;
%         Au = mvrk2(a, x, tol, 'y0', Au, 'verb', 0);
%         real_res = norm(Au-y)/norm(y);
%         fprintf('=als_rake_solve= sweep %d, \t\t\t real_res: %3.3e\n', swp, real_res);
        if (max_res<tol)
            break;
%             kickrank = 0;
%             last_sweep = true;
        end;
%         if (max_res<tol)&&(kickrank<=-als_iters)
%             last_sweep=true;
%         end;
%         if (max_res<=tol*als_tol_high)
%             if (max_res>max_res_prev); regurg_cnt=regurg_cnt+1; fprintf('---- Regurgitation %d\n', regurg_cnt); end;
%             if ((regurg_cnt>als_iters)||(max_res<tol)); last_sweep=true; end;
% %             if ((regurg_cnt>0)||(max_res<=tol*als_tol_low))&&(kickrank>0); kickrank=-1; end;
%         end;
    end;

    if (swp==nswp-1)||(regurg_cnt>2)
        last_sweep=true;
%         tol2 = tol;
    end;

    max_res_prev = max_res;
    max_dx_prev = max_dx;

    max_res = 0;
    max_dx = 0;
    max_frank=0;
end;

% Stuff back
for i=1:d
    x.tuck{i} = cell2core(x.tuck{i}, xf{i}(1:L(i)));
    xc{i} = xf{i}{L(i)+1};
end;
x.core = cell2core(x.core, xc);
x.dphys = d;
% x.sz = n;
end

function [y]=bfun3(Phi1,B1,Phi2, x)
% Computes (Phi1 * B1 * Phi2)*x
% Phi1 is of sizes ry1, rB1, rx1
% B1 is of sizes rB1, k1, m1, rB2
% Phi2 is of sizes ry2, rB2, rx2
ry1 = size(Phi1,1); ry2 = size(Phi2,1);
rx1 = size(Phi1,3); rx2 = size(Phi2,3);
rb1=size(B1,1); rb2=size(B1,4);
m1 = size(B1,3);
k1 = size(B1,2);

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

function [u,v,max_res,max_dx,flg]=local_solve(Phi1,A1,Phi2, phiy1,y1,phiy2, tol, resid_damp, trunc_norm, sol_prev, local_prec, local_restart, local_iters, max_full_size, max_res, max_dx, dir, rmax, kickrank, verb)
rx1 = size(Phi1,1);
n = size(A1,2);
rx2 = size(Phi2,1);
ra1 = size(Phi1,2);
ra2 = size(Phi2,2);

ry1 = size(phiy1,2);
ry2 = size(phiy2,2);
if (dir>0)
    rhs = reshape(y1, ry1, n*ry2);
    rhs = phiy1*rhs;
    rhs = reshape(rhs, rx1*n, ry2);
    if (kickrank>0); y_save = rhs; end;
    rhs = rhs*(phiy2.');
else
    rhs = phiy2.';
    rhs = reshape(y1, ry1*n, ry2)*rhs;
    rhs = reshape(rhs, ry1, n*rx2);
    if (kickrank>0); y_save = rhs; end;
    rhs = phiy1*rhs;
end;
rhs = rhs(:);

norm_rhs = norm(rhs, 'fro');

if (rx1*n*rx2<max_full_size)
    %      |     |    |
    % B = Phi1 - A1 - Phi2
    %      |     |    |
    B = reshape(permute(Phi1, [1, 3, 2]), rx1*rx1, ra1);
    B = B*reshape(A1, ra1, n*n*ra2);
    B = reshape(B, rx1, rx1, n, n, ra2);
    B = permute(B, [1, 3, 2, 4, 5]);
    B = reshape(B, rx1*n*rx1*n, ra2);
    B = B*reshape(permute(Phi2, [2, 1, 3]), ra2, rx2*rx2);
    B = reshape(B, rx1*n, rx1*n, rx2, rx2);
    B = permute(B, [1, 3, 2, 4]);
    B = reshape(B, rx1*n*rx2, rx1*n*rx2);

    if (norm_rhs==0.0)
        % Ground state
        res_prev = norm(B*sol_prev);
        if (res_prev>tol)
            B2 = B+eye(rx1*n*rx2);
            sol_prev2 = sol_prev;
            for it=1:local_restart
                sol = B2 \ sol_prev2;
                sol = sol/norm(sol);
                res_new = norm(B*sol);
                if (strcmp(trunc_norm, 'fro'))
                    if (norm(sol-sol_prev2)<tol); break; end;
                else
                    if (res_new<tol); break; end;
                end;
                sol_prev2 = sol;
            end;
            flg = 0;
        else
            sol = sol_prev;
            res_new = res_prev;
            flg=0;
        end;
    else
        res_prev = norm(B*sol_prev-rhs)/norm_rhs;
        if (res_prev>tol)
            sol = B \ rhs;
            % If the system was ill-conditioned
            %         [sol,flg] = gmres(B, rhs, local_restart, real_tol, 2, [], [], sol);
            res_new = norm(B*sol-rhs)/norm_rhs;
            flg=0;
            if (res_new>tol); flg=1; end;
        else
            sol = sol_prev;
            res_new = res_prev;
            flg=0;
        end;
    end;
else
    if (norm_rhs==0.0)
        res_prev = norm(bfun3(Phi1, A1, Phi2, sol_prev));
        if (res_prev>tol)
            trunc_norm_char = 1;
            local_prec_char = 0;
            if ((strcmp(local_prec, 'cjacobi'))); local_prec_char = 1; end;
            if (strcmp(local_prec, 'ljacobi')); local_prec_char = 2;  end;
            if (strcmp(local_prec, 'rjacobi')); local_prec_char = 3;  end;

            Psi1 = zeros(rx1, rx1, ra1+1);
            Psi1(:,:,1:ra1) = permute(Phi1, [1,3,2]);
            Psi1(:,:,ra1+1) = eye(rx1);
            B1 = zeros(n, n, ra1+1, ra2+1);
            B1(:,:,1:ra1,1:ra2) = permute(A1, [2,3,1,4]);
            B1(:,:,ra1+1,ra2+1) = eye(n);
            B1 = permute(B1, [3,1,2,4]);
            Psi2 = zeros(rx2,rx2,ra2+1);
            Psi2(:,:,1:ra2) = permute(Phi2, [1,3,2]);
            Psi2(:,:,ra2+1) = eye(rx2);
            Psi2 = permute(Psi2, [2,3,1]);

            sol_prev2 = sol_prev;
            for it=1:local_restart
                sol = solve3d_2(Psi1, B1, Psi2, sol_prev2, tol, trunc_norm_char, sol_prev2, local_prec_char, local_restart, local_iters, 1);
                sol = sol/norm(sol);
                res_new = norm(bfun3(Phi1, A1, Phi2, sol));
                if (strcmp(trunc_norm, 'fro'))
                    if (norm(sol-sol_prev2)<tol); break; end;
                else
                    if (res_new<tol); break; end;
                end;
                sol_prev2 = sol;
            end;

            flg=0;
            if (res_new>tol); flg=1; end;
        else
            sol = sol_prev;
            res_new = res_prev;
            flg = 0;
        end;
    else
        res_prev = norm(bfun3(Phi1, A1, Phi2, sol_prev) - rhs)/norm_rhs;
        if (res_prev>tol)
            trunc_norm_char = 1;
%             if (strcmp(trunc_norm, 'fro')); trunc_norm_char = 0; end;
            local_prec_char = 0;
            if ((strcmp(local_prec, 'cjacobi'))); local_prec_char = 1; end;
            if (strcmp(local_prec, 'ljacobi')); local_prec_char = 2;  end;
            if (strcmp(local_prec, 'rjacobi')); local_prec_char = 3;  end;
            if (res_prev>1)
                sol_prev = zeros(rx1*n*rx2, 1);
                res_prev = 1;
            end;

%             sol = solve3d(permute(Phi1,[1,3,2]), A1, permute(Phi2, [1,3,2]), rhs, tol, trunc_norm_char, sol_prev, local_prec_char, local_restart, local_iters, 1);
            sol = solve3d_2(permute(Phi1,[1,3,2]), A1, permute(Phi2, [3,2,1]), rhs, tol, trunc_norm_char, sol_prev, local_prec_char, local_restart, local_iters, 0);

            res_new = norm(bfun3(Phi1, A1, Phi2, sol) - rhs)/norm_rhs;
            flg=0;
            if (res_new>tol); flg=1; end;
        else
            sol = sol_prev;
            res_new = res_prev;
            flg = 0;
        end;
    end;
end;

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

% if (flg>0)&&(kickrank==0)
%     sol = sol_prev;
%     res_new = res_prev;
% end;

dx = norm(sol-sol_prev)/norm(sol);
max_dx = max(max_dx, dx);

max_res = max(max_res, res_prev);

if (norm_rhs==0.0)
    norm_rhs=1;
end;

% Truncation
if (dir>=0) % left-to-right
    sol = reshape(sol, rx1*n, rx2);
else
    sol = reshape(sol, rx1, n*rx2);
end;

if (kickrank>=0)
[u,s,v]=svd(sol, 'econ');
s = diag(s);

if (strcmp(trunc_norm, 'fro')) % We are happy with L2 truncation (when? but let it be)
    r = my_chop2(s, tol*resid_damp*norm(s));
else
    % Residual trunc; First, bin-search
    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*n*rx2<max_full_size)
            res = norm(B*cursol(:)-rhs)/norm_rhs;
        else
            res = norm(bfun3(Phi1, A1, Phi2, cursol)-rhs)/norm_rhs;
        end;
        if (res<max(tol, res_new)*resid_damp)
            r2 = r;
        else
            r1 = r;
        end;
        r = round((r1+r2)/2);
    end;
    r = max(r-1,1);
    % More accurate Linear search
    while (r<=numel(s))
        cursol = u(:,1:r)*diag(s(1:r))*(v(:,1:r)');
        if (rx1*n*rx2<max_full_size)
            res = norm(B*cursol(:)-rhs)/norm_rhs;
        else
            res = norm(bfun3(Phi1, A1, Phi2, cursol)-rhs)/norm_rhs;
        end;
        if (res<max(tol, res_new)*resid_damp)
            break;
        end;
        r = r+1;
    end;
end;

r = min(r, numel(s));
r = min(r, rmax);

else
    if (dir>=0)
        [u,v]=qr(sol, 0);
        v = v';
        r = size(u,2);
        s = ones(r,1);
    else
        [v,u]=qr(sol.', 0);
        u = u.';
        v = conj(v);
        r = size(v,2);
        s = ones(r,1);
    end;
end;

if (verb>1)
    fprintf('=als_rake_solve= dir %d, dx: %3.3e, res_prev: %3.3e, res_new: %3.3e r: %d\n', dir, dx, res_prev, res_new, r);
end;

    if (dir>0) % left-to-right, kickrank, etc
        u = u(:,1:r);
        v = conj(v(:,1:r))*diag(s(1:r));

        if (kickrank>0)
            % Smarter kick: low-rank PCA in residual
            % Matrix: Phi1-A{i}, rhs: Phi1-y{i}, sizes rx(i)*n - ra(i+1)
            leftresid = reshape(Phi1, rx1*ra1, rx1)*reshape(u*v.', rx1, n*rx2);
            leftresid = reshape(leftresid, rx1, ra1*n, rx2);
            leftresid = reshape(permute(leftresid, [2, 1, 3]), ra1*n, rx1*rx2);
            leftresid = reshape(permute(A1, [2,4,1,3]), n*ra2, ra1*n)*leftresid;
            leftresid = reshape(leftresid, n, ra2, rx1, rx2);
            leftresid = reshape(permute(leftresid, [3,1,2,4]), rx1*n, ra2*rx2);
%             leftA = permute(Phi1, [1, 3, 2]);
%             leftA = reshape(leftA, rx1*rx1, ra1);
%             leftA = leftA*reshape(A1, ra1, n*n*ra2);
%             leftA = reshape(leftA, rx1, rx1, n, n, ra2);
%             leftA = permute(leftA, [1, 3, 5, 2, 4]);
%             leftA = reshape(leftA, rx1*n*ra2, rx1*n);
%             leftresid = leftA*reshape(u*v.', rx1*n, rx2);
%             leftresid = reshape(leftresid, rx1*n, ra2*rx2);
            leftresid = [leftresid, -y_save];
%
%             uk = zeros(rx1*n, min(kickrank,rx1*n));
%             for i=1:min(kickrank, rx1*n)
%                 uk2 = uchol(leftresid.', 2);
%                 uk(:,i) = uk2(:,end);
%                 [uk(:,1:i),rv]=qr(uk(:,1:i), 0);
%                 leftresid = leftA*uk(:,i);
%                 leftresid = reshape(leftresid, rx1*n, ra2);
%             end;

% %             [uk,~,~]=svd(leftresid, 'econ');
% %             uk = uk(:,1:min(kickrank, size(uk,2)));
            uk = uchol(leftresid.', kickrank+1);
            uk = uk(:,size(uk,2):-1:max(size(uk,2)-kickrank+1,1));
%             leftresid = leftA*uk;
%             leftresid = reshape(leftresid, rx1*n, ra2*size(uk,2));
%             uk(:,size(uk,2)+1) = uchol(leftresid.', 1);
% %             uk = randn(rx1*n, kickrank);

            [u,rv]=qr([u,uk], 0);
            radd = size(uk,2);
            v = [v, zeros(rx2, radd)];
            v = v*(rv.');
        end;
    elseif (dir<0) % right-to-left
        u = u(:,1:r)*diag(s(1:r));
        v = conj(v(:,1:r));

        if (kickrank>0)
            % Smarter kick: low-rank PCA in residual
            % Matrix: Phi1-A{i}, rhs: Phi1-y{i}, sizes rx(i)*n - ra(i+1)
            rightresid = reshape(Phi2, rx2*ra2, rx2)*(reshape(u*v.', rx1*n, rx2).');
            rightresid = reshape(rightresid, rx2, ra2, rx1, n);
            rightresid = reshape(permute(rightresid, [4, 2, 3, 1]), n*ra2, rx1*rx2);
            rightresid = reshape(A1, ra1*n, n*ra2)*rightresid;
            rightresid = reshape(rightresid, ra1, n, rx1, rx2);
            rightresid = reshape(permute(rightresid, [2,4,1,3]), n*rx2, ra1*rx1);
%             rightA = permute(Phi2, [2, 1, 3]);
%             rightA = reshape(rightA, ra2, rx2*rx2);
%             rightA = reshape(A1, ra1*n*n, ra2)*rightA;
%             rightA = reshape(rightA, ra1, n, n, rx2, rx2);
%             rightA = permute(rightA, [2, 4, 1, 3, 5]);
%             rightA = reshape(rightA, n*rx2*ra1, n*rx2);
%             rightresid = rightA*(reshape(u*v.', rx1, n*rx2).');
%             rightresid = reshape(rightresid, n*rx2, ra1*rx1);
            rightresid = [rightresid, -(y_save.')];

%             uk = zeros(n*rx2, min(kickrank,n*rx2));
%             for i=1:min(kickrank, n*rx2)
%                 uk2 = uchol(rightresid.', 2);
%                 uk(:,i) = uk2(:,end);
%                 [uk(:,1:i),rv]=qr(uk(:,1:i), 0);
%                 rightresid = rightA*uk(:,i);
%                 rightresid = reshape(rightresid, n*rx2, ra1);
%             end;

% %             [uk,~,~]=svd(rightresid, 'econ');
% %             uk = uk(:,1:min(kickrank, size(uk,2)));

% %             uk = randn(n*rx2, kickrank);
            uk = uchol(rightresid.', kickrank+1);
            uk = uk(:,size(uk,2):-1:max(size(uk,2)-kickrank+1,1));

            [v,rv]=qr([v,uk], 0);
            radd = size(uk,2);
            u = [u, zeros(rx1, radd)];
            u = u*(rv.');
        end;
    else
        % Just stuff back the last core
        u = u(:,1:r);
        v = conj(v(:,1:r))*diag(s(1:r));
    end;

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

author: oseledets

repo: TT-Toolbox


text:
function tagS = patient_module_tags
%"patient_module_tags"
%   Return the tags used to represent a patient as specified by C.7.1.1 in
%   PS3.3 of 2006 DICOM specification. Tags are returned in a struct array
%   with 3 fields:
%
%   Tag: String containing hex DICOM tag of a field.
%  Type: String describing type of field, with 5 options:
%         '1' Field must exist, data must exist and be valid.
%         '2' Field must exist, data can exist and be valid, or be NULL.
%         '3' Field is optional, if the field exists data can exist and be
%             valid or be NULL.
%         '1C' Field must exist under certain conditions and contain valid
%             data.
%         '2C' Field must exist under certain conditions and can contain
%             valid data or be NULL.
%Children: For sequences, a tagS with the same format as this struct array.
%
%JRA 06/06/06
%YWU Modified 03/01/08
%
%Usage:
%   tagS = patient_module_tags
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

% {SudyInstanceUID, Patient's Name, PatientID, Issuer of Patient ID, Patient's Birth Date,
% PatientSex, Patient's Birth Time, Other Patient IDs, Other Patient Names,
% Ethnic Group, Patient Comments}

tagC = {'0020000D','00100010','00100020','00100021','00100030',...
'00100040','00100032','00101000','00101001','00102160','00104000'};

tagdecC = {2097165, 1048592, 1048608, 1048609, 1048624, 1048640, 1048626,...
    1052672, 1052673, 1057120, 1064960};

fieldNameC = {'StudyInstanceUID','PatientName','PatientID',...
    'IssuerOfPatientID','PatientBirthDate','PatientSex',...
    'PatientBirthTime','OtherPatientIDs','OtherPatientNames',...
    'EthnicGroup','PatientComments'};

tagTypeC = {'1','2','2','3','2','2','3','3','3','3','3'};

childC = {[],[],[],[],[],[],[],[],[],[],[]};

tagS = struct('tag', tagC, 'tagdec', tagdecC, 'fieldname', fieldNameC, 'type', tagTypeC, 'children', childC);

%
% %Initialize the tagS structure.
% tagS = struct('tag', {}, 'type', {}, 'children', {});
%
% %Create an empty tagS template for sequence creation.
% template = tagS;
%
% %Add tags based on PS3.3 attribute lists.
%
% %Study Instance UID
% tagS(end+1) = struct('tag', ['0020000D'], 'type', ['1'], 'children', []);
%
% %Patient's Name
% tagS(end+1) = struct('tag', ['00100010'], 'type', ['2'], 'children', []);
%
% %Patient ID
% tagS(end+1) = struct('tag', ['00100020'], 'type', ['2'], 'children', []);
%
% %Issuer of Patient ID
% tagS(end+1) = struct('tag', ['00100021'], 'type', ['3'], 'children', []);
%
% %Patient's Birth Date
% tagS(end+1) = struct('tag', ['00100030'], 'type', ['2'], 'children', []);
%
% %Patient's Sex
% tagS(end+1) = struct('tag', ['00100040'], 'type', ['2'], 'children', []);
%
% % %Referenced Patient Sequence
% % tagS(end+1) = struct('tag', ['00081140'], 'type', ['3'], 'children', []);
% %     child_1        = template;
% %
% %     %Referenced SOP Class UID
% %     child_1(end+1) = struct('tag', ['00081150'], 'type', ['1'], 'children', []);
% %
% %     %Referenced SOP Instance UID
% %     child_1(end+1) = struct('tag', ['00081155'], 'type', ['1'], 'children', []);
% %     tagS(end).children = child_1;
%
% %Patient's Birth Time
% tagS(end+1) = struct('tag', ['00100032'], 'type', ['3'], 'children', []);
%
% %Other Patient IDs
% tagS(end+1) = struct('tag', ['00101000'], 'type', ['3'], 'children', []);
%
% %Other Patient Names
% tagS(end+1) = struct('tag', ['00101001'], 'type', ['3'], 'children', []);
%
% %Ethnic Group
% tagS(end+1) = struct('tag', ['00102160'], 'type', ['3'], 'children', []);
%
% %Patient Comments
% tagS(end+1) = struct('tag', ['00104000'], 'type', ['3'], 'children', []);
%
% %wy %Patient Identity Removed
% % tagS(end+1) = struct('tag', ['00120062'], 'type', ['3'], 'children', []);
% %
% % %De-identification Method
% % tagS(end+1) = struct('tag', ['00120063'], 'type', ['1C'], 'children', []);
% %
% % %De-identification Method Code Sequence
% % tagS(end+1) = struct('tag', ['00120064'], 'type', ['1C'], 'children', []);
% %
% % %Include "Code Sequence Macro"
% % tagS(end).children = code_sequence_macro_tags;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function y=horzcat(varargin)
%% TT-Toolbox 2.2.1, 2009-2016
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel

%for questions on this particular code email Alexey Boyko, Skolkovo Institute of Science and Technology,
% a.boyko@skoltech.ru or alexey.boyko@skolkovotech.ru
%Moscow, Russia, 2016
%%

%now you can assemble rectangular matrices horizontally like
%ttAA=[ttA ttA]

%Attention!
%Due to limitations of MATLAB, for simultaneous horizontal and vertical concatenation
% use
% ttM=tt_blockwise({ttA ttB ttC; ttD ttE ttF ;ttG ttH ttI })
% instead of
% ttM=[ttA ttB ttC; ttD ttE ttF ;ttG ttH ttI ]

y=tt_blockwise(varargin);
end


text:
function model = Learning_WLS_HP(Seqs, model, alg, flag)

D = size(model.A, 1);
L = length(model.landmark);
C = length(Seqs);
XX = [];
YY = [];

for c = 1:C

    Time = Seqs(c).Time;
    if ~isempty(Time)
    Event = Seqs(c).Mark;
    Tstart = Seqs(c).Start;

%     if isempty(alg.Tmax)
%         Tstop = Seqs(c).Stop;
%     else
%         Tstop = alg.Tmax;
%         indt = Time < alg.Tmax;
%         Time = Time(indt);
%         Event = Event(indt);
%     end

    Y = zeros(length(Time), 1);

    if strcmp(flag, 'multisource')
        Xmu = zeros(length(Time), D*C);
    else
        Xmu = zeros(length(Time), D);
    end
    XA = zeros(length(Time), D*L*D);

    for i = 1:length(Time)
        di = Event(i);
        if strcmp(flag, 'multisource')
            Xmu(i, D*(c-1)+di) = 1;
        else
            Xmu(i, di)= 1;
        end

        index = find(Event(1:i)==di);
        Y(i) = length(index);

        if i>1
            dt = Time(i) - Time(1:i-1);
            G = Kernel_Integration(dt, model);

            for j = 1:i-1
                dj = Event(j);
                Gij = G(j,:);
                XA(i, D*L*(di-1)+(dj:L:D*L)) = ...
                    XA(i, D*L*(di-1)+(dj:L:D*L)) + Gij';
            end
        end
    end

    Ytmp = alg.w.*Y./sqrt((Time(:) - Tstart));
    Xtmp = alg.w.*[Xmu, XA]./repmat(sqrt(Time(:) - Tstart), ...
        [1, size(Xmu,2)+size(XA,2)]);
    YY = [YY; Ytmp];
    XX = [XX; Xtmp];
    end
end

res0 = rand(size(XX,2),1);
options = optimoptions('fmincon', 'Algorithm', 'interior-point', ...
    'Display', 'iter', ...
    'MaxFunctionEvaluations', 30000);
res = fmincon(@(x)norm(XX*x - YY)^2, res0, [], [], [], [], ...
    zeros(size(XX,2),1), [], [], options);
%res = XX\YY;
if strcmp(flag, 'multisource')
    model.mu = reshape(res(1:C*D), [D, C]);
    model.A = reshape(res(1+C*D:end), [D, L, D]);
else
    model.mu = reshape(res(1:D), [D, 1]);
    model.A = reshape(res(1+D:end), [D, L, D]);
end

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
% % Search the elements that are lying within a specified interval.
% % Search the indexes of all elements in x (sorted vector of n elements) that
% % lie within the interval.
% % The algorithm uses binary searches, thus it runs in log(n)
% %
% % INPUT:
% % x: vector of numeric values sorted in ascending order
% %    (e.g. 2,7,20,...120)
% % ref: numeric value of the reference point (center of the interval)
% % tol: numeric value corresponding to 1/2 of the width of the interval
% % The fourth input argument: numeric value (optional). Allows to define the maximum
% % number of elements of x that can lie within the specified interval. This is useful
% % in order to speed up the search.
% %
% % OUTPUT:
% % indexes: indexes of elements of x which lie within the interval
% % [ref-tol ref+tol]
% % If ref is not found in x then indexes is empty.
% %
% % Revisions:
% % 14-jun-2010: - Added lines 38-55 in order to consider the case in which lower bound
% %                is less than the minimum value of x or the upper bound is more than
% %                the max value of x
% %              - Reinforced the conditions in line 93 and in line 117 in order to avoid
% %                an error when (to-1)==0 or (from+1)>length(x)
% % 15/apr/2011: - A bug with uint values of input array has been fixed.
% %                Thanks to Igor Varfolomeev

% % % --------------------------------
% % % Author: Dr. Roberto Olmi
% % % Email : robertoolmi at gmail.com
% % % --------------------------------

function indexes = BinaryIntervalSearch(x,ref,tol,varargin)
lbound = ref-tol;
ubound = ref+tol;
from=1;
to=length(x);

if lbound <= x(1)
    if ubound < x(1)
        indexes=[];
        return
    end
    lindex=1;
else
    lindex=0;
end
if ubound >= x(end)
    if lbound > x(end)
        indexes=[];
        return
    end
    uindex=length(x);
else
    uindex=0;
end

go=uindex==0 || lindex==0;
while go
    mid = ceil(from+(to-from)/2);
    %diff = x(mid)-ref;  15/apr/2011
    if x(mid) < ref %diff < 0  15/apr/2011
        if x(mid) >= lbound
            go=false;
        else
            from=mid+1;
            if x(from) >= lbound
                lindex=from;
                go=false;
            end
        end
    else       % x(mid) > ref
        if x(mid) <= ubound
            go=false;
        else
            to=mid-1;
            if x(to) <= ubound
                uindex=to;
                go=false;
            end
        end
    end
end

%remove this if at least one element of x is always in the interval:
if (lindex > 0 && x(lindex) > ubound)...
        || (uindex > 0 && x(uindex) < lbound)
    indexes=[];
    return
end

%search upper index
cfrom=from;
if from==length(x) || x(from+1) > ubound
    uindex=from;
end
if nargin == 4
    to=min([to mid+varargin{1}]);
end
while uindex == 0
    mid = ceil(from+(to-from)/2);
    if x(mid) <= ubound
        from=mid+1;
        if x(from) > ubound
            uindex=mid;
        end
    else
        to=mid-1;
        if x(to) < ubound
            uindex=to;
        end
    end
end

%search lower index
from=cfrom;
to=uindex;
if to==1 || x(to-1) < lbound
    lindex=to;
end
if nargin == 4
    from=max([from uindex-varargin{1}]);
end
while lindex == 0
    mid = ceil(from+(to-from)/2);
    if x(mid) < lbound
        from=mid+1;
        if x(from) > lbound
            lindex=from;
        end
    else
        to=mid-1;
        if x(to) < lbound
            lindex=mid;
        end
    end
end

indexes=lindex:uindex;

% % --------------------------------------------
% % Example code for Testing
% % tol=2;
% % ref=12;
% % maxi=10000;
% % numel=1000;
% % for i=1:maxi
% %     x = sort(randint(1,numel,[0,400]));
% %     indexes = BinaryIntervalSearch(x,ref,tol);
% %     if isempty (indexes)
% %         if any(abs(x-ref) < tol)
% %             disp('Doh!')
% %             break
% %         else
% %             continue
% %         end
% %     end
% %     if indexes(1)>1 && ~(x(indexes(1)-1) < x(indexes(1)))...
% %             || indexes(end)<numel && ~(x(indexes(end)) < x(indexes(end)+1))...
% %             || any(abs(x(indexes)-ref) > tol)
% %         disp('Doh!')
% %         break
% %     end
% %     if i==maxi
% %         disp('OK!!!')
% %     end
% % end
% % --------------------------------------------

% % % --------------------------------
% % % Author: Dr. Roberto Olmi
% % % Email : robertoolmi at gmail.com
% % % --------------------------------


text:
% run cell detection on spatial masks U
function [ops, stat, model] = sourcery(ops)

tic
[ops, U0, model, U2]    = get_svdForROI(ops);

%  U0 = my_conv2(U0, ops.sig, [1 2]);

ops.fig         = getOr(ops, 'fig', 1);
ops.ThScaling   = getOr(ops, 'ThScaling', 1);

% reshape U to be (nMaps x Y x X)
U0 =  reshape(U0, [], size(U0,ndims(U0)))';
Ly = numel(ops.yrange);
Lx = numel(ops.xrange);
[nSVD, Npix] = size(U0);
U0 = reshape(U0, nSVD, Ly, Lx);

% compute neuropil basis functions for cell detection
S = getNeuropilBasis(ops, Ly, Lx, 'Fourier'); % 'raisedcosyne', 'Fourier'
S = normc(S);
nBasis = size(S,2);

StU = S'*U0(:,:)'; % covariance of neuropil with spatial masks
StS = S'*S; % covariance of neuropil basis functions

% make cell mask with ops.diameter
d0   = ceil(ops.diameter); % expected cell diameter in pixels
sig = ceil(d0/4);
dx = repmat([-d0:d0], 2*d0+1, 1);
dy = dx';
rs = dx.^2 + dy.^2 - d0^2;
dx = dx(rs<=0);
dy = dy(rs<=0);

% initialize cell matrices
mPix    = zeros(numel(dx), 1e4);
mLam    = zeros(numel(dx), 1e4);

iter = 0;
icell = 0;
r = rand(1e4,1);

L   = sparse(Ly*Lx, 0);
LtU = zeros(0, nSVD);
LtS = zeros(0, nBasis);

% regress maps onto basis functions and subtract neuropil contribution
% U = Ucell + neu'*S'
% neu = inv(S'*S) * (S'*U')
neu     = StS\StU;
Ucell   =  U0 - reshape(neu' * S', size(U0));

nBasis = size(S,2);

%
while 1
    iter = iter + 1;

    % residual is smoothed at every iteration
    us = my_conv2_circ(Ucell, sig, [2 3]);

    % compute log variance at each location
    V = sq(mean(us.^2,1));
    V = double(V);

    um = sq(mean(Ucell.^2,1));
    um = my_conv2_circ(um, sig, [1 2]);

    V = V./um ;
    %     V = log(V./um);
    V = double(V);
    % do the morphological opening trick
    % take the running max of the running min
    % this normalizes the brightness of the image
    if iter==1
        lbound = -my_min2(-my_min2(V, d0), d0);
    end

    V = V - lbound;

    if iter==1
        % find indices of all maxima  in plus minus 1 range
        % use the median of these peaks to decide stopping criterion
        maxV    = -my_min(-V, 1, [1 2]);
        ix      = (V > maxV-1e-10);

        % threshold is the mean peak, times a potential scaling factor
        pks = V(ix);

        Th  = ops.ThScaling * median(pks(pks>1e-4));

        ops.Vcorr = V;
    end

    % just in case this goes above original value
    V = min(V, ops.Vcorr);

    % find local maxima in a +- d neighborhood
    maxV = -my_min(-V, d0, [1 2]);

    % find indices of these maxima above a threshold
    ix  = (V > maxV-1e-10) & (V > Th);
    ind = find(ix);

    if iter==1
       Nfirst = numel(ind);
    end

    if numel(ind)==0
        break;
    end

    new_codes = normc(us(:, ind));

    ncells = icell;
    LtU(ncells+size(new_codes,2), nSVD) = 0;

    % each source needs to be iteratively subtracted off
    for i = 1:size(new_codes,2)
        icell = icell + 1;
        [ipix, ipos] = getIpix(ind(i), dx, dy, Lx, Ly);

        Usub = Ucell(:, ipix);

        lam = max(0, new_codes(:, i)' * Usub);

        % threshold pixels
        lam(lam<max(lam)/5) = 0;

        mPix(ipos,icell) = ipix;
        mLam(ipos,icell) = lam;

        % extract biggest connected region of lam only
        mLam(:,icell)   = normc(getConnected(mLam(:,icell), rs)); % ADD normc HERE and BELOW!!!
        lam             = mLam(ipos,icell) ;

        L(ipix,icell)   = lam;

        LtU(icell, :)   = U0(:,ipix) * lam;
        LtS(icell, :)   = lam' * S(ipix,:);
    end

    % ADD NEUROPIL INTO REGRESSION HERE
    LtL     = full(L'*L);
    codes   = ([LtL LtS; LtS' StS]+ 1e-3 * eye(icell+nBasis))\[LtU; StU];
    neu     = codes(icell+1:end,:);
    codes   = codes(1:icell,:);
%     codes = (LtL+ 1e-3 * eye(icell))\LtU;

    % subtract off everything
    Ucell = U0 - reshape(neu' * S', size(U0)) - reshape(double(codes') * L', size(U0));

    % re-estimate masks
    L   = sparse(Ly*Lx, icell);
    for j = 1:icell
        ipos = find(mPix(:,j)>0);
        ipix = mPix(ipos,j);

        Usub = Ucell(:, ipix)+ codes(j, :)' * mLam(ipos,j)';

        lam = max(0, codes(j, :) * Usub);
        % threshold pixels
        lam(lam<max(lam)/5) = 0;

        mLam(ipos,j) = lam;

        % extract biggest connected region of lam only
        mLam(:,j) = normc(getConnected(mLam(:,j), rs));
        lam = mLam(ipos,j);

        L(ipix,j) = lam;

        LtU(j, :) = U0(:,ipix) * lam;
        LtS(j, :) = lam' * S(ipix,:);

        Ucell(:, ipix) = Usub - (Usub * lam)* lam';
    end
    %
    err(iter) = mean(Ucell(:).^2);

    Vnew = sq(sum(Ucell.^2,1));

    fprintf('%d total ROIs, err %4.4f, thresh %4.4f \n', icell, err(iter), Th)
    if ops.fig

        figure(1)
        subplot(1,2, 1);
        imagesc(ops.Vcorr, [0 2*Th])
        axis off

        subplot(1,2, 2);
        imagesc(V, [0 2*Th])
        axis off

        figure(2)
        [~, iclust, lam] = drawClusters(ops, r, mPix, mLam, Ly, Lx);

        drawnow
    end

     if (numel(ind)<Nfirst * getOr(ops, 'stopSourcery', 1/10)) || (iter>= getOr(ops, 'maxIterRoiDetection', 100))
        break;
    end
end

% this runs only the mask re-estimation step, on non-smoothed PCs
% (because smoothing is done during clustering to help)
if getOr(ops, 'refine', 1)
	sourceryAddon;
else
	fprintf('no refinement of ROIs done\n')
end

mLam  =  mLam(:, 1:icell);
mPix  =  mPix(:, 1:icell);

mLam = bsxfun(@rdivide, mLam, sum(mLam,1));
%%

% subtract off neuropil only
Ucell = U0 - reshape(neu' * S', size(U0));

% populate stat with cell locations and footprint
stat = getFootprint(ops, codes, Ucell, mPix, mLam);

% compute compactness of ROIs
stat = anatomize(ops, mPix, mLam, stat);

[~, iclust, lam] = drawClusters(ops, r, mPix, mLam, Ly, Lx);

model.L     = L;
model.S     = S;
model.LtS   = LtS;
model.LtL   = LtL;
model.StS   = StS;

% get anatomical projection weights
stat = weightsMeanImage(ops, stat, model);

author: cortex-lab


text:
T_10 = readtable('/Users/lbasek/image-denoising-benchmark/results/final_results_sigma10.csv');
T_10(:,'WNNM') = [];
T_10(:,'NCSR') = [];

T_15 = readtable('/Users/lbasek/image-denoising-benchmark/results/final_results_sigma15.csv');
T_15(:,'WNNM') = [];
T_15(:,'NCSR') = [];

T_20 = readtable('/Users/lbasek/image-denoising-benchmark/results/final_results_sigma20.csv');
T_20(:,'WNNM') = [];
T_20(:,'NCSR') = [];

T_25 = readtable('/Users/lbasek/image-denoising-benchmark/results/final_results_sigma25.csv');
T_25(:,'WNNM') = [];
T_25(:,'NCSR') = [];

T_50 = readtable('/Users/lbasek/image-denoising-benchmark/results/final_results_sigma50.csv');
T_50(:,'WNNM') = [];
T_50(:,'NCSR') = [];

mse_row_10 = T_10(1,:);
mse_row_15 = T_15(1,:);
mse_row_20 = T_20(1,:);
mse_row_25 = T_25(1,:);
mse_row_50 = T_50(1,:);

x = [0,10,15,20,25,50];

bm3d_10 = mse_row_10{1,3};
bm3d_15 = mse_row_15{1,3};
bm3d_20 = mse_row_20{1,3};
bm3d_25 = mse_row_25{1,3};
bm3d_50 = mse_row_50{1,3};

ksvd_10 = mse_row_10{1,4};
ksvd_15 = mse_row_15{1,4};
ksvd_20 = mse_row_20{1,4};
ksvd_25 = mse_row_25{1,4};
ksvd_50 = mse_row_50{1,4};

foe_10 = mse_row_10{1,5};
foe_15 = mse_row_15{1,5};
foe_20 = mse_row_20{1,5};
foe_25 = mse_row_25{1,5};
foe_50 = mse_row_50{1,5};

epll_10 = mse_row_10{1,6};
epll_15 = mse_row_15{1,6};
epll_20 = mse_row_20{1,6};
epll_25 = mse_row_25{1,6};
epll_50 = mse_row_50{1,6};

y1 = [mse_row_10{1,2}, bm3d_10,bm3d_15,bm3d_20,bm3d_25,bm3d_50];
y2 = [mse_row_10{1,2}, ksvd_10,ksvd_15,ksvd_20,ksvd_25,ksvd_50];
y3 = [mse_row_10{1,2}, foe_10,foe_15,foe_20,foe_25,foe_50];
y4 = [mse_row_10{1,2}, epll_10,epll_15,epll_20,epll_25,epll_50];

h = figure;
hold on
h1 = plot(x,y1,'-s','MarkerSize',10);
h2 = plot(x,y2,'-s','MarkerSize',10);
h3 = plot(x,y3,'-s','MarkerSize',10);
h4 = plot(x,y4,'-s','MarkerSize',10);
hold off
set([h1,h2,h3,h4],'LineWidth',2)
ylabel('MSE')
xlabel('Noise Power')
legend('BM3D','KSVD','FOE','EPLL');
grid on
saveas(h,'./graphs/mse','png');

author: lbasek

repo: image-denoising-benchmark


text:
function [varargout]=discQuadMesh(varargin)

% function [F,V,C,indEdge]=discQuadMesh(nElements,r,f)
% ------------------------------------------------------------------------
% This function meshes a circle using quadrilaterial elements.
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 10/05/2016 Updated for GIBBON
%------------------------------------------------------------------------

%% Parse input
switch nargin
    case 1
        nElements=varargin{1};
        r=1;
        f=0.5;
    case 2
        nElements=varargin{1};
        r=varargin{2};
        f=0.5;
    case 3
        nElements=varargin{1};
        r=varargin{2};
        f=varargin{3};
end

%% Creating central regular quad mesh

nElements=nElements+~iseven(nElements);%Force even

[X_centralMesh,Y_centralMesh]=meshgrid(linspace(-1,1,nElements+1));
[F_centralMesh,V_centralMesh] = surf2patch(X_centralMesh,Y_centralMesh,zeros(size(X_centralMesh)));
V_centralMesh=V_centralMesh(:,1:2);

%Edge of central mesh
logicCentralMeshEdge=(X_centralMesh==1)|(Y_centralMesh==1)|(X_centralMesh==-1)|(Y_centralMesh==-1);
nEdge=(nElements*4);

% Scaling radius
[ThetaMesh,RadiusMesh]=cart2pol(V_centralMesh(:,1),V_centralMesh(:,2));
RadiusMesh=f*(1/2)*sqrt(2)*RadiusMesh;
[V_centralMesh(:,1),V_centralMesh(:,2)]=pol2cart(ThetaMesh,RadiusMesh);

%% Creating outer mesh

RadiusOuterEdge=ones(1,nEdge);
ThetaOuterEdge=linspace(0,pi*2,nEdge+1);
ThetaOuterEdge=ThetaOuterEdge(2:end)-pi;

[xOuterEdge,yOuterEdge]=pol2cart(ThetaOuterEdge,RadiusOuterEdge);
V_outerEdge=[xOuterEdge(:) yOuterEdge(:)];

V_innerEdge=V_centralMesh(logicCentralMeshEdge,:);
[ThetaEdge,RadiusEdge]=cart2pol(V_innerEdge(:,1),V_innerEdge(:,2));
[ThetaEdge,sortInd]=sort(ThetaEdge);
RadiusEdge=RadiusEdge(sortInd);
[V_innerEdge(:,1),V_innerEdge(:,2)]=pol2cart(ThetaEdge,RadiusEdge);

[Xr]=linspacen(V_innerEdge(:,1),V_outerEdge(:,1),nElements/2+1); Xr(end+1,:)=Xr(1,:);
[Yr]=linspacen(V_innerEdge(:,2),V_outerEdge(:,2),nElements/2+1); Yr(end+1,:)=Yr(1,:);

[Fs2,Vs2] = surf2patch(Xr,Yr,zeros(size(Xr)));
Vs2=Vs2(:,1:2);

V=[V_centralMesh;Vs2];
F=[F_centralMesh;Fs2+size(V_centralMesh,1)];
C=[ones(size(F_centralMesh,1),1); 2*ones(size(Fs2,1),1); ];

indEdge=((size(V,1)-size(Xr,1))+1):size(V,1);

%% Removing double points

[F,V,~,IND_IND]=mergeVertices(F,V);
indEdge=IND_IND(indEdge(1:end-1));

%Scaling radius
[ThetaMesh,RadiusMesh]=cart2pol(V(:,1),V(:,2));
RadiusMesh=r*RadiusMesh;
[V(:,1),V(:,2)]=pol2cart(ThetaMesh,RadiusMesh);

varargout{1}=F;
varargout{2}=V;
varargout{3}=C;
varargout{4}=indEdge;

%%
% _*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 IncrementalSequence < handle

    properties (GetAccess = public, SetAccess = protected)
        value
    end

    properties (Access = public)
        alpha
        initialValue
        finalValue
    end

    properties (Access = protected)
        x0
        x1
        nSteps
    end

    methods (Access = protected, Abstract)

        generateAlphaSequence(obj)

    end

    methods (Access = public)

        function update(obj,i)
            obj.value = (1-obj.alpha(i))*obj.initialValue + obj.alpha(i)*obj.finalValue;
        end

    end

    methods (Access = protected)

        function init(obj,x0,x1,nSteps,a0,a1)
            obj.x0 = x0;
            obj.x1 = x1;
            obj.nSteps = nSteps;
            obj.initialValue = a0;
            obj.finalValue = a1;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [hog2x2_arr, info] = dense_hog2x2(I, HOG2x2param)
% Adapted from LabelMeToolbox
% URL: http://labelme.csail.mit.edu/Release3.0/browserTools/php/matlab_toolbox.php

grid_spacing = HOG2x2param.grid_spacing;
patch_size = HOG2x2param.patch_size;
half_win = floor(patch_size/2);

I = double(I);
d = pixelwise_hog31(I,half_win);
xmax = size(d,1); ymax = size(d,2);

grid_x = (1:grid_spacing:xmax-half_win);
grid_y = (1:grid_spacing:ymax-half_win);

hog2x2_arr = zeros([length(grid_x),length(grid_y),124]);

hog2x2_arr(:,:, 1: 31) = d(grid_x         ,grid_y         ,:);
hog2x2_arr(:,:,32: 62) = d(grid_x+half_win,grid_y         ,:);
hog2x2_arr(:,:,63: 93) = d(grid_x         ,grid_y+half_win,:);
hog2x2_arr(:,:,94:124) = d(grid_x+half_win,grid_y+half_win,:);

[x y] = meshgrid(grid_x, grid_y);
info.y = x(:);
info.x = y(:);

author: adikhosla

repo: feature-extraction


text:
function [U,output] = btd_rnd(size_tens,size_core,options)
%BTD_RND Pseudorandom initialization for BTD.
%   U = btd_rnd(size_tens,size_core) generates R pseudorandom terms U{r} to
%   initialize algorithms that compute a block term decomposition of an
%   N-th order tensor. Each term U{r} is a cell array of N unitary factor
%   matrices U{r}{n}, followed by a core tensor U{r}{N+1} of size
%   size_core{r}.
%
%   btd_rnd(T,size_core) is shorthand for btd_rnd(size(T),size_core) if
%   T is real. If T is complex, then by default the terms will be generated
%   using pseudorandom complex numbers as well (cf. options).
%
%   btd_rnd(size_tens,size_core,options) and btd_rnd(T,size_core, ...
%   options) may be used to set the following options:
%
%      options.Real =         - The type of random number generator used to
%      [{@randn}|@rand|0]       generate the real part of the core tensor S
%                               and matrices U{n}. If 0, there is no real
%                               part.
%      options.Imag =         - The type of random number generator used to
%      [@randn|@rand|0|...      generate the imaginary part of the core
%       {'auto'}]               tensor S and matrices U{n}. If 0, there is
%                               no imaginary part. On 'auto', options.Imag
%                               is 0 unless the first argument is a complex
%                               tensor T, in which case it is equal to
%                               options.Real.
%      options.Unitary = true - Set equal to false to generate the factor
%                               matrices U{n} without converting them to
%                               unitary matrices afterwards.
%
%   See also cpd_rnd, lmlra_rnd, btdgen.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)

if nargin < 3, options = struct; end
if ~iscell(size_core), size_core = {size_core}; end
function U = capture(size_core)
    [U,S] = lmlra_rnd(size_tens,size_core,options);
    U{end+1} = S;
end
U = cellfun(@capture,size_core,'UniformOutput',false);
output = struct;

end

author: andrewssobral

repo: lrslibrary


text:
function fig()
figure


text:
function [ algs ] = sgd_solver_list(category)
% Return list of solvers.
%
% Inputs:
%       category    category to be returned.
% Output:
%       algs        list of solvers in the category
%
% This file is part of GDLibrary.
%
% Created by H.Kasai on Nov. 1, 2016
% Modified by H.Kasai on Sep. 28, 2017

    % supported algorithms by SGDLibrary
    basic_algs = {'SGD','SGD-CM','SGD-CM-NAG','AdaGrad','RMSProp','AdaDelta','Adam','AdaMax'};
    vr_algs = {'SVRG','SAG','SAGA','SARAH'};
    qnewton_algs = {'SQN','SVRG-SQN','SVRG-LBFGS','SS-SVRG', ...
        'oBFGS-Inf','oBFGS-Lim','Reg-oBFGS-Inf','Reg-oBFGS-Lim','Damp-oBFGS-Inf','Damp-oBFGS-Lim'};
    else_algs = {'SVRG-BB','IQN'};

    switch category
        case 'BASIC'
            algs = basic_algs;
        case 'VR'
            algs = vr_algs;
        case 'QN'
            algs = qnewton_algs;
        case 'ALL'
            algs = [basic_algs, vr_algs, qnewton_algs, else_algs];
        otherwise
    end

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function L = log_prob_node(CPD, self_ev, pev)
% LOG_PROB_NODE Compute prod_m log P(x(i,m)| x(pi_i,m), theta_i) for node i (gaussian)
% L = log_prob_node(CPD, self_ev, pev)
%
% self_ev(m) is the evidence on this node in case m.
% pev(i,m) is the evidence on the i'th parent in case m (if there are any parents).
% (These may also be cell arrays.)

if iscell(self_ev), usecell = 1; else usecell = 0; end

use_log = 1;
ncases = length(self_ev);
nparents = length(CPD.sizes)-1;
assert(ncases == size(pev, 2));

if ncases == 0
  L = 0;
  return;
end

L = 0;
for m=1:ncases
  if isempty(CPD.dps)
    i = 1;
  else
    if usecell
      dpvals = cat(1, pev{CPD.dps, m});
    else
      dpvals = pev(CPD.dps, m);
    end
    i = subv2ind(CPD.sizes(CPD.dps), dpvals(:)');
  end
  if usecell
    y = self_ev{m};
  else
    y = self_ev(m);
  end
  if length(CPD.cps) == 0
    L = L + gaussian_prob(y, CPD.mean(:,i), CPD.cov(:,:,i), use_log);
  else
    if usecell
      x = cat(1, pev{CPD.cps, m});
    else
      x = pev(CPD.cps, m);
    end
    L = L + gaussian_prob(y, CPD.mean(:,i) + CPD.weights(:,:,i)*x, CPD.cov(:,:,i), use_log);
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function grad = getGradient(problem, x, storedb, key)
% Computes the gradient of the cost function at x.
%
% function grad = getGradient(problem, x)
% function grad = getGradient(problem, x, storedb)
% function grad = getGradient(problem, x, storedb, key)
%
% Returns the gradient at x of the cost function described in the problem
% structure.
%
% storedb is a StoreDB object, key is the StoreDB key to point x.
%
% See also: getDirectionalDerivative canGetGradient

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%
%   April 3, 2015 (NB):
%       Works with the new StoreDB class system.
%
%  June 28, 2016 (NB):
%       Works with getPartialGradient.
%
%   Nov. 1, 2016 (NB):
%       Added support for gradient from directional derivatives.
%       Last resort is call to getApproxGradient instead of an exception.
%
%   Sep. 6, 2018 (NB):
%       The gradient is now cached by default. This is made practical by
%       the new storedb 'remove' functionalities that keep the number of
%       cached points down to a minimum. If the gradient is obtained via
%       costgrad, the cost is also cached.
%
%   Feb. 10, 2020 (NB):
%       Allowing M.egrad2rgrad to take (storedb, key) as extra inputs.

    % Allow omission of the key, and even of storedb.
    if ~exist('key', 'var')
        if ~exist('storedb', 'var')
            storedb = StoreDB();
        end
        key = storedb.getNewKey();
    end

    % Contrary to most similar functions, here, we get the store by
    % default. This is for the caching functionality described below.
    store = storedb.getWithShared(key);
    store_is_stale = false;

    % If the gradient has been computed before at this point (and its
    % memory is still in storedb), then we just look up the value.
    force_grad_caching = true;
    if force_grad_caching && isfield(store, 'grad__')
        grad = store.grad__;
        return;
    end

    % We don't normally compute the cost value, but if we get it as a side
    % result, then we may as well take note of it for caching.
    cost_computed = false;

    if isfield(problem, 'grad')
    %% Compute the gradient using grad.

        % Check whether this function wants to deal with storedb or not.
        switch nargin(problem.grad)
            case 1
                grad = problem.grad(x);
            case 2
                [grad, store] = problem.grad(x, store);
            case 3
                % Pass along the whole storedb (by reference), with key.
                grad = problem.grad(x, storedb, key);
                % The store structure in storedb might have been modified
                % (since it is passed by reference), so before caching
                % we'll have to update (see below).
                store_is_stale = true;
            otherwise
                up = MException('manopt:getGradient:badgrad', ...
                    'grad should accept 1, 2 or 3 inputs.');
                throw(up);
        end

    elseif isfield(problem, 'costgrad')
    %% Compute the gradient using costgrad.

        % Check whether this function wants to deal with storedb or not.
        switch nargin(problem.costgrad)
            case 1
                [cost, grad] = problem.costgrad(x);
            case 2
                [cost, grad, store] = problem.costgrad(x, store);
            case 3
                % Pass along the whole storedb (by reference), with key.
                [cost, grad] = problem.costgrad(x, storedb, key);
                store_is_stale = true;
            otherwise
                up = MException('manopt:getGradient:badcostgrad', ...
                    'costgrad should accept 1, 2 or 3 inputs.');
                throw(up);
        end

        cost_computed = true;

    elseif canGetEuclideanGradient(problem)
    %% Compute the Riemannian gradient using the Euclidean gradient.

        egrad = getEuclideanGradient(problem, x, storedb, key);
        % Convert to the Riemannian gradient
        switch nargin(problem.M.egrad2rgrad)
            case 2
                grad = problem.M.egrad2rgrad(x, egrad);
            case 4
                grad = problem.M.egrad2rgrad(x, egrad, storedb, key);
            otherwise
                up = MException('manopt:getGradient:egrad2rgrad', ...
                    'egrad2rgrad should accept 2 or 4 inputs.');
                throw(up);
        end
        store_is_stale = true;

    elseif canGetPartialGradient(problem)
    %% Compute the gradient using a full partial gradient.

        d = problem.ncostterms;
        grad = getPartialGradient(problem, x, 1:d, storedb, key);
        store_is_stale = true;

    elseif canGetDirectionalDerivative(problem)
    %% Compute gradient based on directional derivatives; expensive!

        B = tangentorthobasis(problem.M, x);
        df = zeros(size(B));
        for k = 1 : numel(B)
            df(k) = getDirectionalDerivative(problem, x, B{k}, storedb, key);
        end
        grad = lincomb(problem.M, x, B, df);
        store_is_stale = true;

    else
    %% Attempt the computation of an approximation of the gradient.

        grad = getApproxGradient(problem, x, storedb, key);
        store_is_stale = true;

    end

    % If we are not sure that the store structure is up to date, update.
    if store_is_stale
        store = storedb.getWithShared(key);
    end

    % Cache here.
    if force_grad_caching
        store.grad__ = grad;
    end
    % If we got the gradient via costgrad, then the cost has also been
    % computed and we can cache it.
    if cost_computed
        store.cost__ = cost;
    end

    storedb.setWithShared(store, key);

end

author: NicolasBoumal

repo: manopt


text:
clc; clear all; close all;

meshPath = '../data/woody.obj';
[V,F] = readOBJ(meshPath);
V = V(:,1:2); % the input mesh contains redundant zero third column

% get handles
fig = tsurf(F,V);
axis equal;
fprintf( ...
    ['Point Handle Selection: \n' ...
    '- CLICK the mesh to add point handls \n', ...
    '- BACKSPACE to remvoe the previous selection\n', ...
    '- ENTER to finish selection\n'] ...
    );
try
  [Cx,Cy] = getpts;
catch e
  return  % quit early, stop script
end

C = [Cx,Cy]; % handle locations

% compute pairwise distance
D = zeros(size(V,1), size(C,1));
for ii = 1:size(C,1)
    D(:,ii) = sqrt(sum((V - C(ii,:)).^2,2));
end

% geet handle indices b
[~,b] = min(D);

% TODO: (unbounded) biharmonic weights
W = compute_skinning_weight(V,F,b);

% TODO: linear blended skinning
deform_GUI(V,F,C,W);

author: odedstein

repo: sgi-introduction-course


text:
function [eigvec, eigval] = PTKVectorisedEigenvalues(M, eigenvalues_only)
    % PTKVectorisedEigenvalues. Computes eigenvalues and eigenvectors for many symmetric matrices
    %
    %     PTKVectorisedEigenvalues is similar to Matlab's eigs() function, but can be
    %     used to compute the eigenvalues and eigenvectors for multiple matrices,
    %     which for a large number of points is significantly quicker than using a
    %     for loop. Each input matrix must be symmetric and is represented by a
    %     single row of the input matrix as described below.
    %
    %     The mex function PTKFastEigenvalues is equivalent to function but runs
    %     faster. PTKVectorisedEigenvalues is slower than PTKFastEigenvalues but still
    %     significantly faster than running eigs() in a for loop when a large
    %     number of matrices is involved.
    %
    %
    %         Syntax:
    %             [eigvectors, eigvalues] = PTKFastEigenvalues(M [, eigenvalues_only])
    %
    %         Input:
    %             M is a 6xn matrix. Each column of M represents one 3x3 symmetric matrix as follows
    %
    %                     [V(1) V(2) V(3); V(2) V(4) V(5); V(3) V(5) V(6)]
    %
    %                 where V is a 6x1 column of M
    %
    %             eigenvalues_only is an optional parameter which defaults to
    %                 false. Set to true to only calculate eigenvalues and not
    %                 eigenvectors, which reduces the execution time.
    %
    %          Outputs:
    %              eigenvalues is a 3xn matrix. Each column contains the 3
    %                  eigenvalues of the matrix V described above
    %
    %              eigenvectors is a 3x3xn matrix, where each 3x1 row represents
    %                  an eigenvector (3 for each of the n matrices V described
    %                  above).
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    compute_eigenvectors = true;
    if exist('eigenvalues_only', 'var')
        if (eigenvalues_only)
            compute_eigenvectors = false;
        end
    end

    numvoxels = size(M, 2);
    eigval = zeros(3, numvoxels, 'double');

    m = (M(1,:) + M(4,:) + M(6,:))/3;

    q =  (( M(1,:) - m) .* (M(4,:) - m) .* (M(6,:) - m) + ...
        2 * M(2,:) .* M(5,:) .* M(3,:) - ...
        M(3,:).^2 .* (M(4,:) - m) - ...
        M(5,:).^2 .* (M(1,:) - m) - M(2,:).^2 .* (M(6,:) - m) )/2;

    p = ( ( M(1,:) - m ).^2 + 2 * M(2,:).^2 + 2 * M(3,:).^2 + ...
        ( M(4,:) - m ).^2 + 2 * M(5,:).^2 + ( M(6,:) - m ).^2 )/6;

    p = max(0.01, p);

    phi = 1/3*acos(min(1, q./p.^(3/2)));
    phi(phi<0) = phi(phi<0)+pi/3;

    eigval(1,:) = m + 2*sqrt(p).*cos(phi);
    eigval(2,:) = m - sqrt(p).*(cos(phi) + sqrt(3).*sin(phi));
    eigval(3,:) = m - sqrt(p).*(cos(phi) - sqrt(3).*sin(phi));

    [~, i] = sort(abs(eigval));
    i = i + size(eigval,1)*ones(size(eigval,1), 1)*(0:size(eigval, 2) - 1);
    eigval = eigval(i);

    % Only compute the eigenvectors if requested
    if (compute_eigenvectors)

        eigvec = zeros(3,3, numvoxels, 'double');
        for l = 1 : 2
            Ai = M(1,:) - eigval(l,:);
            Bi = M(4,:) - eigval(l,:);
            Ci = M(6,:) - eigval(l,:);

            eix = ( M(2,:) .* M(5,:) - Bi .* M(3,:) ) .* ( M(3,:) .* M(5,:) - Ci .* M(2,:) );
            eiy = ( M(3,:) .* M(5,:) - Ci .* M(2,:) ) .* ( M(3,:) .* M(2,:) - Ai .* M(5,:) );
            eiz = ( M(2,:) .* M(5,:) - Bi .* M(3,:) ) .* ( M(3,:) .* M(2,:) - Ai .* M(5,:) );

            vec = sqrt(eix.^2+eiy.^2+eiz.^2);
            vec = max(0.01, vec);
            eigvec(:,l,:) = [eix; eiy; eiz] ./ vec([1;1;1], :);
        end

        eigvec(:,3,:) = cross(eigvec(:,1,:), eigvec(:,2,:));
    else
        eigvec = [];
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [ yv, yvp ] = hermite_interpolant_value ( nd, xd, yd, xdp, ydp, nv, ...
  xv )

%*****************************************************************************80
%
%% HERMITE_INTERPOLANT_VALUE evaluates the Hermite interpolant polynomial.
%
%  Discussion:
%
%    In fact, this function will evaluate an arbitrary polynomial that is
%    represented by a difference table.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carl deBoor,
%    A Practical Guide to Splines,
%    Springer, 2001,
%    ISBN: 0387953663,
%    LC: QA1.A647.v27.
%
%  Parameters:
%
%    Input, integer ND, the order of the difference table.
%
%    Input, real XD(ND), YD(ND), the difference table for the
%    interpolant value.
%
%    Input, real XDP(ND-1), YDP(ND-1), the difference table for
%    the interpolant derivative.
%
%    Input, integer NV, the number of evaluation points.
%
%    Input, real XV(NV), the evaluation points.
%
%    Output, real YV(NV), YVP(NV), the value of the interpolant and
%    its derivative at the evaluation points.
%
  ndp = nd - 1;
%
%  Clearly MATLAB is going to be all prissy about combining
%  vectors and scalars and wants me to use repmat or something
%  nonintuitive if I want to vectorize this loop and I say
%  it's spinach and to hell with it.
%
  for j = 1 : nv

    yv(j) = yd(nd);
    for i = nd - 1 : -1 : 1
      yv(j) = yd(i) + ( xv(j) - xd(i) ) * yv(j);
    end

    yvp(j) = ydp(ndp);
    for i = ndp - 1 : -1 : 1
      yvp(j) = ydp(i) + ( xv(j) - xdp(i) ) * yvp(j);
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% This file is part of the project NILM-Eval (https://github.com/beckel/nilm-eval).
% Licence: GPL 2.0 (http://www.gnu.org/licenses/gpl-2.0.html)
% Copyright: ETH Zurich, 2014
% Author: Romano Cicchetti

warning off MATLAB:dispatcher:nameConflict

addpath(genpath('algorithms'));
addpath(genpath('config'));
addpath(genpath('data_access'));
addpath(genpath('framework'));
% some conflict with 'legend' plotting function :(
% addpath(genpath('lib/bnt'));
addpath(genpath('lib/lightspeed'));
addpath(genpath('lib/proxTV-1.0'));
addpath(genpath('lib/spectralClustering'));
addpath(genpath('lib/TVDIP'));
addpath(genpath('lib/YAMLMatlab_0.4.3'));
addpath(genpath('plot'));
addpath('projects/');
addpath('projects/buildsys/');
addpath('projects/buildsys/images/');
addpath(genpath('util'));

set(0,'DefaultTextFontname', 'Times New Roman');
set(0,'DefaultTextFontSize', 9);
set(0,'DefaultAxesFontName', 'Times New Roman');
set(0,'DefaultAxesFontSize', 9);

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval