text:
function computeAllModelChoiceTime_HN(pathSet,pathModels,training,param,maxOrder,nBoot,seed)
% -------------------------------------------------------------------------
% function computeAllModelChoice_HN(pathSet,pathModels,training,param,maxOrder,nBoot,imbalance,batchNum)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes feature set selection for a given feature set type
% and for all experiments with different degrees of freedom. See ref. [1,2]
% for more details.
% -------------------------------------------------------------------------
% REFERENCES:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% [2] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathSet: Full path to the Feature Set folder of the corresponding experiment.
%             --> Ex: '/myProject/WORKSPACE/COHORT-BASED-RESULTS/Experiment1/FSET'
% 2. pathModels: Full path to the Models folder of the corresponding experiment.
%                --> Ex: '/myProject/WORKSPACE/COHORT-BASED-RESULTS/Experiment1/MODELS'
% 3. training: Structure defining all parameters for the given experiment
%              to perform. See masterScript_HN.m for more details.
% 4. param: Cell of two strings, defining 1) The feature set type/name; and
%           2) the outcome to model
% 5. maxOrder: Integer specifying the maximal model order to construct.
%              --> Ex: 10
% 6. nBoot: Number of bootstrap samples to use.
%           --> Ex: 100
% 7. imbalance: String specifying the type of imbalance-adjustement strategy
%               employed. Either 'IABR' for imbalance-adjusted bootstrap
%               resampling (see ref.[1]), or 'IALR' for imbalance-adjusted
%               logistic regression (see ref.[2]).
%               --> Ex: 'IALR'
% 8. batchNum: (optional input). If present, integer that specifies the
%              batch number for parallelization purposes.
%              --> Ex: 6
% -------------------------------------------------------------------------
% OUTPUTS: Mutivariable models are saved in a folder named 'MODELS' in the
% corresponding folder of "pathExperiment"
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: March 2016
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
warning off

cd(pathSet)
fSetName = param{1};
outcomeName = param{2};
outcome = training.(outcomeName).timeToEvent;
censoring = 1 - training.(outcomeName).outcome;
fSet = load(['FSET_',fSetName,'_',outcomeName]); fSet = struct2cell(fSet); fSet = fSet{1};

tic
fprintf(['\n--> SELECTING FEATURES (MODEL ORDERS OF 1 to % u) FOR "',fSetName,'" FEATURE SET AND "',outcomeName,'" OUTCOME ... '],maxOrder)
[models] = featureSelectionTime_HN(fSet.Data,outcome,maxOrder,nBoot,fSet.Info,censoring,seed);
cd(pathModels), save(['MODELS_',fSetName,'_',outcomeName],'models')
fprintf('DONE!\n'), toc

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function model = rmGridFit_oneGaussian(model,prediction,data,params,t)
% rmGridFit_oneGaussian - core of one Gaussian fit
%
% model = rmGridFit_oneGaussian(model,prediction,data,params);
%
% 2008/01 SOD: split of from rmGridFit.

% input check
if nargin < 4,
    error('Not enough arguments');
end

% some variables we need
%rssinf         = inf(size(data(1,:)),'single');
trends         = t.trends;
t_id           = t.dcid+1;

% we compute mean rss from the sum rss (old convention)
% Some explanation of the divisor would be helpful here.
model.rss = single(model.rss./(size(prediction,1)-(size(trends,2)+1)));

%-----------------------------------
%--- fit different receptive fields profiles
%--- another loop --- and a slow one too
%-----------------------------------
tic; progress = 0;

%warning('off', 'MATLAB:lscov:RankDefDesignMat')

for n=1:numel(params.analysis.x0),
    %-----------------------------------
    % progress monitor (10 dots) and time indicator
    %-----------------------------------
    if floor(n./numel(params.analysis.x0).*10)>progress,
        if progress==0,
            % print out estimated time left
            esttime = toc.*10;
            if floor(esttime./3600)>0,
                fprintf(1,'[%s]:Estimated processing time: %d hours.\t(%s)\n',...
                    mfilename, ceil(esttime./3600), datestr(now));
            else
                fprintf(1, '[%s]:Estimated processing time: %d minutes.\t(%s)\n',...
                    mfilename, ceil(esttime./60), datestr(now));
            end;
            fprintf(1,'[%s]:Grid (x,y,sigma) fit:',mfilename);drawnow;
        end;
        % progress monitor
        fprintf(1,'.');drawnow;
        progress = progress + 1;
    end;

    %-----------------------------------
    %--- now apply glm to fit RF
    %-----------------------------------
    % minimum RSS fit

    X    = [prediction(:,n) trends];  % Build the model matrix

    % This line takes up 30% of the time
    % lscov takes as long as the pinv method but provides the rss as
    % well... Let's put in the ~.  All modern Matlab now uses it, and it
    % will save us memory space.
    [b,ci,rss]    = lscov(X,data);

    % Compute RSS only for positive fits. The basic problem is
    % that if you have two complementary locations, you
    % could fit with a postive beta on the one that drives the signal or a
    % negative beta on the portion of the visual field that never sees the
    % stimulus. This would produce the same prediction. We don't like that
    nkeep   = (b(1,:)<0); % Now we only set the negative fits to inf.

    % To save time limit the rss computation to those we care about.
    rss(nkeep) = inf('single');

    %-----------------------------------
    %--- store data with lower rss
    %-----------------------------------
    minRssIndex = (rss < model.rss);  % This is where we use model.rss

    % now update model parameters for those cases in which the computed rss
    % is smaller than the input model.rss
    model.x0(minRssIndex)       = params.analysis.x0(n);
    model.y0(minRssIndex)       = params.analysis.y0(n);
    model.s(minRssIndex)        = params.analysis.sigmaMajor(n);
    model.s_major(minRssIndex)        = params.analysis.sigmaMajor(n);
    model.s_minor(minRssIndex)        = params.analysis.sigmaMajor(n);
    model.s_theta(minRssIndex)        = params.analysis.theta(n);
    model.rss(minRssIndex)      = rss(minRssIndex);
    model.b([1 t_id],minRssIndex) = b(:,minRssIndex);
end

%warning('on', 'MATLAB:lscov:RankDefDesignMat')

% Under some conditions, the grid fit never returns an acceptable fit, For
% example for onegaussian fits with data driven DC component, when the DC
% is artificially high. In this case some of the rss values remain Inf,
% which fails to interpolate and compute correct variance explained values.
% So we check it here and reset any Inf (bad fits) to rawrss, so the
% variance explained will be 0.
model.rss(model.rss==Inf) = model.rawrss(model.rss==Inf);

% Correct lscov. It returns the mean rss. To maintain compatibility with
% the sum rss that this function expects, we have to multiply by the
% divisor. See the lscov docs for details.
model.rss = single(model.rss.*(size(prediction,1)-(size(trends,2)+1)));

% end time monitor
et  = toc;
if floor(esttime/3600)>0,
    fprintf(1,'Done[%d hours].\t(%s)\n', ceil(et/3600), datestr(now));
else
    fprintf(1,'Done[%d minutes].\t(%s)\n', ceil(et/60), datestr(now));
end;
drawnow;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [varargout] = spm_diff(varargin)
% matrix high-order numerical differentiation
% FORMAT [dfdx] = spm_diff(f,x,...,n)
% FORMAT [dfdx] = spm_diff(f,x,...,n,V)
% FORMAT [dfdx] = spm_diff(f,x,...,n,'q')
%
% f      - [inline] function f(x{1},...)
% x      - input argument[s]
% n      - arguments to differentiate w.r.t.
%
% V      - cell array of matrices that allow for differentiation w.r.t.
% to a linear transformation of the parameters: i.e., returns
%
% df/dy{i};    x = V{i}y{i};    V = dx(i)/dy(i)
%
% q      - (char) flag to preclude default concatenation of dfdx
%
% dfdx          - df/dx{i}                     ; n =  i
% dfdx{p}...{q} - df/dx{i}dx{j}(q)...dx{k}(p)  ; n = [i j ... k]
%
%
% This routine has the same functionality as spm_ddiff, however it
% uses one sample point to approximate gradients with numerical (finite)
% differences:
%
% dfdx  = (f(x + dx)- f(x))/dx
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_diff.m 7143 2017-07-29 18:50:38Z karl $

% step size for numerical derivatives
%--------------------------------------------------------------------------
global GLOBAL_DX
if ~isempty(GLOBAL_DX)
    dx = GLOBAL_DX;
else
    dx = exp(-8);
end

% create inline object
%--------------------------------------------------------------------------
f     = spm_funcheck(varargin{1});

% parse input arguments
%--------------------------------------------------------------------------
if iscell(varargin{end})
    x = varargin(2:(end - 2));
    n = varargin{end - 1};
    V = varargin{end};
    q = 1;
elseif isnumeric(varargin{end})
    x = varargin(2:(end - 1));
    n = varargin{end};
    V = cell(1,length(x));
    q = 1;
elseif ischar(varargin{end})
    x = varargin(2:(end - 2));
    n = varargin{end - 1};
    V = cell(1,length(x));
    q = 0;
else
    error('improper call')
end

% check transform matrices V = dxdy
%--------------------------------------------------------------------------
for i = 1:length(x)
    try
        V{i};
    catch
        V{i} = [];
    end
    if isempty(V{i}) && any(n == i);
        V{i} = speye(spm_length(x{i}));
    end
end

% initialise
%--------------------------------------------------------------------------
m     = n(end);
xm    = spm_vec(x{m});
J     = cell(1,size(V{m},2));

% proceed to derivatives
%==========================================================================
if length(n) == 1

    % dfdx
    %----------------------------------------------------------------------
    f0    = f(x{:});
    for i = 1:length(J)
        xi    = x;
        xi{m} = spm_unvec(xm + V{m}(:,i)*dx,x{m});
        J{i}  = spm_dfdx(f(xi{:}),f0,dx);
    end

    % return numeric array for first-order derivatives
    %======================================================================

    % vectorise f
    %----------------------------------------------------------------------
    f  = spm_vec(f0);

    % if there are no arguments to differentiate w.r.t. ...
    %----------------------------------------------------------------------
    if isempty(xm)
        J = sparse(length(f),0);

    % or there are no arguments to differentiate
    %----------------------------------------------------------------------
    elseif isempty(f)
        J = sparse(0,length(xm));
    end

    % differentiation of a scalar or vector
    %----------------------------------------------------------------------
    if isnumeric(f0) && iscell(J) && q
        J = spm_dfdx_cat(J);
    end

    % assign output argument and return
    %----------------------------------------------------------------------
    varargout{1} = J;
    varargout{2} = f0;

else

    % dfdxdxdx....
    %----------------------------------------------------------------------
    f0        = cell(1,length(n));
    [f0{:}]   = spm_diff(f,x{:},n(1:end - 1),V);
    p         = true;

    for i = 1:length(J)
        xi    = x;
        xmi   = xm + V{m}(:,i)*dx;
        xi{m} = spm_unvec(xmi,x{m});
        fi    = spm_diff(f,xi{:},n(1:end - 1),V);
        J{i}  = spm_dfdx(fi,f0{1},dx);
        p     = p & isnumeric(J{i});
    end

    % or differentiation of a scalar or vector
    %----------------------------------------------------------------------
    if p && q
        J = spm_dfdx_cat(J);
    end
    varargout = [{J} f0];
end

function dfdx = spm_dfdx(f,f0,dx)
% cell subtraction
%__________________________________________________________________________
if iscell(f)
    dfdx  = f;
    for i = 1:length(f(:))
        dfdx{i} = spm_dfdx(f{i},f0{i},dx);
    end
elseif isstruct(f)
    dfdx  = (spm_vec(f) - spm_vec(f0))/dx;
else
    dfdx  = (f - f0)/dx;
end

return

function J = spm_dfdx_cat(J)
% concatenate into a matrix
%--------------------------------------------------------------------------
if isvector(J{1})
    if size(J{1},2) == 1
        J = spm_cat(J);
    else
        J = spm_cat(J')';
    end
end

author: spm


text:
function B = ordfilt2(varargin)
%ORDFILT2 Perform 2-D order-statistic filtering.
%   B=ORDFILT2(A,ORDER,DOMAIN) replaces each element in A by the
%   ORDER-th element in the sorted set of neighbors specified by
%   the nonzero elements in DOMAIN.
%
%   B = ORDFILT2(A,ORDER,DOMAIN,S), where S is the same size as
%   DOMAIN, uses the values of S corresponding to the nonzero
%   values of DOMAIN as additive offsets.
%
%   B = ORDFILT2(...,PADOPT) controls how the matrix boundaries
%   are padded.  PADOPT may be 'zeros' (the default) or
%   'symmetric'.  If PADOPT is 'zeros', A is padded with zeros at
%   the boundaries.  If PADOPT is 'symmetric', A is symmetrically
%   extended at the boundaries.
%
%   Class Support
%   -------------
%   The class of A may be numeric or logical.  The class of B is
%   the same as the class of A, unless the additive offset form of
%   ORDFILT2 is used, in which case the class of B is double.
%
%   Example
%   -------
%   Use a maximum filter on snowflakes.png with a [5 5] neighborhood.  This is
%   equivalent to imdilate(image,strel('square',5)).
%
%       A = imread('snowflakes.png');
%       B = ordfilt2(A,25,true(5));
%       imview(A), imview(B)
%
%   Remarks
%   -------
%   DOMAIN is equivalent to the structuring element used for
%   binary image operations. It is a matrix containing only 1's
%   and 0's; the 1's define the neighborhood for the filtering
%   operation.
%
%   For example, B=ORDFILT2(A,5,ONES(3,3)) implements a 3-by-3
%   median filter; B=ORDFILT2(A,1,ONES(3,3)) implements a 3-by-3
%   minimum filter; and B=ORDFILT2(A,9,ONES(3,3)) implements a
%   3-by-3 maximum filter.  B=ORDFILT2(A,4,[0 1 0; 1 0 1; 0 1 0])
%   replaces each element in A by the maximum of its north, east,
%   south, and west neighbors.
%
%   See also MEDFILT2.

%   Copyright 1993-2003 The MathWorks, Inc.
%   $Revision: 5.17.4.5 $  $Date: 2003/08/23 05:53:07 $

[A,order,domain,s,padopt,msg] = ParseInputs(varargin{:});

domainSize = size(domain);
center = floor((domainSize + 1) / 2);
[r,c] = find(domain);
r = r - center(1);
c = c - center(2);
padSize = max(max(abs(r)), max(abs(c)));
originalSize = size(A);
if (strcmp(padopt, 'zeros'))
    A = padarray(A, padSize * [1 1], 0, 'both');
elseif (strcmp(padopt, 'ones'))
    % padopt of 'ones' is for support of medfilt2; it is
    % undocumented
    A = padarray(A, padSize * [1 1], 1, 'both');
else
    A = padarray(A, padSize * [1 1], 'symmetric', 'both');
end
Ma = size(A,1);
offsets = c*Ma + r;

% make sure that offsets are valid
if ~isreal(offsets) || any(floor(offsets) ~= offsets) || any(~isfinite(offsets))
    %should never get here
    eid = sprintf('Images:%s:internalError', mfilename);
    msg = 'Internal error: bad OFFSETS.';
    error(eid,'%s',msg);
end

if isempty(s)
  %ORDFILT2(A,ORDER,DOMAIN)
  B = ordf(A, order, offsets, [padSize padSize] + 1, ...
             originalSize, domainSize);
else
  %ORDFILT2(A,ORDER,DOMAIN,S,PADOPT)
  B = ordf(A, order, offsets, [padSize padSize] + 1, ...
           originalSize, domainSize, s);
end

%%%
%%% ParseInputs
%%%
function [A,order,domain,s,padopt,msg] = ParseInputs(varargin)

A = [];
order = [];
domain = [];
s = [];
padopt = 'zeros';
msg = '';

% checknargin(3,5,nargin,mfilename);

A = varargin{1};
order = varargin{2};
domain = varargin{3};
options = {'zeros', 'ones', 'symmetric'};
% padopt of 'ones' is for supporting medfilt2; it is undocumented.

if (nargin == 4)
  if (ischar(varargin{4}))
    padopt = checkstrs(varargin{4},options,mfilename,'PADOPT',4);
  else
    s = varargin{4};
  end

elseif (nargin == 5)
  s = varargin{4};
  padopt = checkstrs(varargin{5},options,mfilename,'PADOPT',5);
end

% make sure that arguments are valid
% checkinput(order,'double',{'real','scalar','integer'},mfilename, ...
%           'ORDER',2);

if ~isempty(s)
  if (~isa(A, 'double'))
    A = double(A);
  end
%  checkinput(A, 'double', {'2d','real'}, mfilename, 'A', 1);
  s = s(find(domain));
%  checkinput(s, 'double', 'real', mfilename, 'S', 4);
else
%  checkinput(A, {'numeric','logical'}, {'2d','real'}, mfilename, 'A', 1);
end

author: gpeyre

repo: matlab-toolboxes


text:
function determ = dif1cyclic_determinant ( n )

%*****************************************************************************80
%
%% DIF1CYCLIC_DETERMINANT returns the determinant of the DIF1CYCLIC matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real DETERM, the determinant.
%
  determ = 0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [curr_dat_sz, curr_lab_sz] = store2hdf5(filename, data, labels, create, startloc, chunksz)
  % *data* is W*H*C*N matrix of images should be normalized (e.g. to lie between 0 and 1) beforehand
  % *label* is D*N matrix of labels (D labels per sample)
  % *create* [0/1] specifies whether to create file newly or to append to previously created file, useful to store information in batches when a dataset is too big to be held in memory  (default: 1)
  % *startloc* (point at which to start writing data). By default,
  % if create=1 (create mode), startloc.data=[1 1 1 1], and startloc.lab=[1 1];
  % if create=0 (append mode), startloc.data=[1 1 1 K+1], and startloc.lab = [1 K+1]; where K is the current number of samples stored in the HDF
  % chunksz (used only in create mode), specifies number of samples to be stored per chunk (see HDF5 documentation on chunking) for creating HDF5 files with unbounded maximum size - TLDR; higher chunk sizes allow faster read-write operations

  % verify that format is right
  dat_dims=size(data);
  lab_dims=size(labels);
  num_samples=dat_dims(end);

  assert(lab_dims(end)==num_samples, 'Number of samples should be matched between data and labels');

  if ~exist('create','var')
    create=true;
  end

  if create
    %fprintf('Creating dataset with %d samples\n', num_samples);
    if ~exist('chunksz', 'var')
      chunksz=1000;
    end
    if exist(filename, 'file')
      fprintf('Warning: replacing existing file %s \n', filename);
      delete(filename);
    end
    h5create(filename, '/data', [dat_dims(1:end-1) Inf], 'Datatype', 'single', 'ChunkSize', [dat_dims(1:end-1) chunksz]); % width, height, channels, number
    h5create(filename, '/label', [lab_dims(1:end-1) Inf], 'Datatype', 'single', 'ChunkSize', [lab_dims(1:end-1) chunksz]); % width, height, channels, number
    if ~exist('startloc','var')
      startloc.dat=[ones(1,length(dat_dims)-1), 1];
      startloc.lab=[ones(1,length(lab_dims)-1), 1];
    end
  else  % append mode
    if ~exist('startloc','var')
      info=h5info(filename);
      prev_dat_sz=info.Datasets(1).Dataspace.Size;
      prev_lab_sz=info.Datasets(2).Dataspace.Size;
      assert(all(prev_dat_sz(1:end-1)==dat_dims(1:end-1)), 'Data dimensions must match existing dimensions in dataset');
      assert(all(prev_lab_sz(1:end-1)==lab_dims(1:end-1)), 'Label dimensions must match existing dimensions in dataset');
      startloc.dat=[ones(1,length(dat_dims)-1), prev_dat_sz(end)+1];
      startloc.lab=[ones(1,length(lab_dims)-1), prev_lab_sz(end)+1];
    end
  end

  if ~isempty(data)
    h5write(filename, '/data', single(data), startloc.dat, size(data));
    h5write(filename, '/label', single(labels), startloc.lab, size(labels));
  end

  if nargout
    info=h5info(filename);
    curr_dat_sz=info.Datasets(1).Dataspace.Size;
    curr_lab_sz=info.Datasets(2).Dataspace.Size;
  end
end

author: happynear

repo: FaceVerification


text:
%% HILBERT_SCRIPT is a script to call HILBERT_FUN.
%
%  Discussion:
%
%    The BATCH command runs scripts, not functions.  So we have to write
%    this short script if we want to work with BATCH!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2012
%
%  Author:
%
%    John Burkardt
%
  m = 100;
  n = 100;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HILBERT_SCRIPT\n' );
  fprintf ( 1, '  Call HILBERT_FUN to set up the Hilbert matrix.\n', n );

  tic;
  a = hilbert_fun ( m, n );
  toc;


text:
function img_mrtrix=tensor_model_2(gx,gy,gz)
% This script convert vector file to tensor visulization in MRtrix
% vectorFile='Average_Vn2_VectorFile.mat';

%Select eignvector
%Gradient 1, 2 or 3
%1 is eigenvector 2, 2 is eigenvector 3, 3 is eigenvector 4

%Scale factor for all values in the tensor
%Affects all models
ScaleFactor=0.05; %Model 2

%Balance between isotropic and anistropic component in Model 2
%1 fully anisotropic; 0 fully isotropic
%Only affects Model 2
BalanceFactor=0.995;

% Model 2: Outlier detection method with clipping

gx1=gx; gy1=gy; gz1=gz;

%Size of image
N=size(gx);

%Initialize tensors for the three models
for i=1:6
    tnsr{i}=zeros(N); %tensor image
end

mag1=sqrt(gx1.^2+gy1.^2+gz1.^2);

%Mask of non-zero eigenvectors
msk=zeros(N);
ind=find(mag1);
msk(ind)=1;

%Eigenvector magnitudes
mag_vec1=mag1(ind);

%Scale eignvector magnitudes
meth='quartile'; t=5;
mag_vec1x=filloutliers(mag_vec1,'clip',meth,'ThresholdFactor',t);

%V=3;
for i=1:length(ind)
    v1=[gx1(ind(i));gy1(ind(i));gz1(ind(i))];
    v1=v1/mag_vec1(i); %normalize

    vv=v1;
    mag_vecx=mag_vec1x;

    % scale anistropic component by scaled eigenvector magnitude
    tmp2=ScaleFactor*( BalanceFactor*mag_vecx(i)*(vv*vv') + (1-BalanceFactor)*eye(3) );

    tnsr{1}(ind(i))=tmp2(1,1); tnsr{2}(ind(i))=tmp2(1,2);
    tnsr{3}(ind(i))=tmp2(2,2); tnsr{4}(ind(i))=tmp2(1,3);
    tnsr{5}(ind(i))=tmp2(2,3); tnsr{6}(ind(i))=tmp2(3,3);

end

% MRtrix order D11, D22, D33, D12, D13, D23
img_mrtrix=cat(4,tnsr{1},tnsr{3},tnsr{6},tnsr{2},tnsr{4},tnsr{5});

author: yetianmed

repo: subcortex

save_path: github-repos/MATLAB/yetianmed-subcortex


text:
function bc = spm_cfg_eeg_bc
% configuration file for baseline correction
%__________________________________________________________________________
% Copyright (C) 2009-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_cfg_eeg_bc.m 5652 2013-09-25 09:36:22Z volkmar $

%--------------------------------------------------------------------------
% D
%--------------------------------------------------------------------------
D        = cfg_files;
D.tag    = 'D';
D.name   = 'File Name';
D.filter = 'mat';
D.num    = [1 1];
D.help   = {'Select the M/EEG mat file.'};

%--------------------------------------------------------------------------
% timewin
%--------------------------------------------------------------------------
timewin         = cfg_entry;
timewin.tag     = 'timewin';
timewin.name    = 'Baseline';
timewin.help    = {'Start and stop of baseline [ms].'};
timewin.strtype = 'r';
timewin.num     = [1 2];

%--------------------------------------------------------------------------
% prefix
%--------------------------------------------------------------------------
prefix         = cfg_entry;
prefix.tag     = 'prefix';
prefix.name    = 'Filename Prefix';
prefix.help    = {'Specify the string to be prepended to the filenames of the filtered dataset. Default prefix is ''b''.'};
prefix.strtype = 's';
prefix.num     = [1 Inf];
prefix.val     = {'b'};

%--------------------------------------------------------------------------
% bc
%--------------------------------------------------------------------------
bc          = cfg_exbranch;
bc.tag      = 'bc';
bc.name     = 'Baseline correction';
bc.val      = {D, timewin, prefix};
bc.help     = {'Baseline correction of M/EEG time data'}';
bc.prog     = @eeg_bc;
bc.vout     = @vout_eeg_bc;
bc.modality = {'EEG'};

%==========================================================================
function out = eeg_bc(job)
% construct the S struct
S = job;
S.D = S.D{1};

out.D          = spm_eeg_bc(S);
out.Dfname     = {fullfile(out.D.path,out.D.fname)};

%==========================================================================
function dep = vout_eeg_bc(job)
% return dependencies
dep(1)            = cfg_dep;
dep(1).sname      = 'Baseline corrected M/EEG data';
dep(1).src_output = substruct('.','D');
dep(1).tgt_spec   = cfg_findspec({{'strtype','e'}});

dep(2)            = cfg_dep;
dep(2).sname      = 'Baseline corrected M/EEG datafile';
dep(2).src_output = substruct('.','Dfname');
dep(2).tgt_spec   = cfg_findspec({{'filter','mat'}});


text:
function region = p17_region ( )

%*****************************************************************************80
%
%% P17_REGION returns the name of the integration region for problem 17.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string REGION, the name of the integration region.
%
  region = 'box';

  return
end


text:
function [F_S frames_index] = bss_make_frames(S,WINDOW,NOVERLAP)

% decompose some signal(s) into frames
%
% Usage: [F_S frames_index] = bss_make_frames(S,WINDOW,NOVERLAP)
%
% Input:
%   - S: matrix of size n x T (with T>n),
%   - WINDOW: 1 x W window
%   - NOVERLAP: number of samples overlap
%
% Output:
%   - F_S:
%       * if n=1, F_S is a n_frames x W matrix containing the frames (of length W) in
%       rows,
%       * if n>1, F_S is a n_frames x W x n tensor containing the frames
%       decomposition of each row of S.
%   - frames_index: index of the beginning of each frame in the rows of S
%
% Developers:  - Cedric Fevotte (fevotte@tsi.enst.fr) - Emmanuel Vincent
% (emmanuel.vincent@irisa.fr) - Remi Gribonval (remi.gribonval@irisa.fr)

[n,T]=size(S);

if n>T
    disp('Wrong dimensions: must have T>n.')
    return;
end

%%% Default values %%%
W=length(WINDOW); % Length of window

if T < W
    disp('Please choose a window smaller than the signals.')
    return;
end

n_frames = fix((T-NOVERLAP)/(W-NOVERLAP)); % Number of frames
% If needed the very end of the signal is removed.

frames_index = 1 + (0:(n_frames-1))*(W-NOVERLAP); % Index of beginnings of frames

F_S=zeros(n_frames,W,n); % If n=1, F_S is a 2-D tensor (matrix)

for i=1:n
    for k=1:n_frames
        F_S(k,:,i)=(S(i,frames_index(k)+(0:(W-1))).*WINDOW);
    end
end

author: posenhuang

repo: deeplearningsourceseparation


text:
function export_3D_model(model,FileName,psz,should_diffuse)

if(isempty(model))
    disp('Error: export_3D_model: Can not write empty data to a file. Skipping.');
    return;
end

if(~exist('psz','var'))
    psz = [];
end

if(~exist('should_diffuse','var'))
    should_diffuse = 0;
end

data.vertex.x = model.x;
data.vertex.y = model.y;
data.vertex.z = model.z;

data.vertex.nx = model.nx;
data.vertex.ny = model.ny;
data.vertex.nz = model.nz;

if(model_has_colors(model))
    if(should_diffuse)
        data.vertex.diffuse_red = model.r*255;
        data.vertex.diffuse_green = model.g*255;
        data.vertex.diffuse_blue = model.b*255;
    else
        data.vertex.red = model.r*255;
        data.vertex.green = model.g*255;
        data.vertex.blue = model.b*255;
    end
end

if(~isempty(psz))
    data.vertex.psz = 0.1*ones(length(model.x),1);
end

if(isfield(model,'face'))
    data.face = model.face;
end

[pathstr, name, ext] = fileparts(FileName);

if(strcmpi(ext,'.ply'))
    nplyWrite(FileName,data);
elseif(strcmpi(ext,'.mat'))
    save(FileName,'data');
else
    error('export_3D_model: Invalid file extension');
end

author: lmb-freiburg

repo: orion

save_path: github-repos/MATLAB/lmb-freiburg-orion


text:
function v = setLayerParams(obj, layer, params)
%SETLAYEPARAMS  Set or change the parameters of a layer
%   Example: NET.SETLAYERPARAMS('layerName', {'param1', 'param2', ...})

v = [] ;
l = obj.getLayerIndex(layer) ;
for param = params
  v(end+1) = obj.addParam(char(param)) ;
end
obj.layers(l).params = params ;
obj.rebuild() ;

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function results = vl_test_whistc(varargin)
% VL_TEST_WHISTC
vl_test_init ;

function test_acc()
x = ones(1, 10) ;
e = 1 ;
o = 1:10 ;
vl_assert_equal(vl_whistc(x, o, e), 55) ;

function test_basic()
x = 1:10 ;
e = 1:10 ;
o = ones(1, 10) ;
vl_assert_equal(histc(x, e), vl_whistc(x, o, e)) ;

x = linspace(-1,11,100) ;
o = ones(size(x)) ;
vl_assert_equal(histc(x, e), vl_whistc(x, o, e)) ;

function test_multidim()
x = rand(10, 20, 30) ;
e = linspace(0,1,10) ;
o = ones(size(x)) ;

vl_assert_equal(histc(x, e), vl_whistc(x, o, e)) ;
vl_assert_equal(histc(x, e, 1), vl_whistc(x, o, e, 1)) ;
vl_assert_equal(histc(x, e, 2), vl_whistc(x, o, e, 2)) ;
vl_assert_equal(histc(x, e, 3), vl_whistc(x, o, e, 3)) ;

function test_nan()
x = rand(10, 20, 30) ;
e = linspace(0,1,10) ;
o = ones(size(x)) ;
x(1:7:end) = NaN ;

vl_assert_equal(histc(x, e), vl_whistc(x, o, e)) ;
vl_assert_equal(histc(x, e, 1), vl_whistc(x, o, e, 1)) ;
vl_assert_equal(histc(x, e, 2), vl_whistc(x, o, e, 2)) ;
vl_assert_equal(histc(x, e, 3), vl_whistc(x, o, e, 3)) ;

function test_no_edges()
x = rand(10, 20, 30) ;
o = ones(size(x)) ;
vl_assert_equal(histc(1, []), vl_whistc(1, 1, [])) ;
vl_assert_equal(histc(x, []), vl_whistc(x, o, [])) ;
vl_assert_equal(histc(x, [], 1), vl_whistc(x, o, [], 1)) ;
vl_assert_equal(histc(x, [], 2), vl_whistc(x, o, [], 2)) ;
vl_assert_equal(histc(x, [], 3), vl_whistc(x, o, [], 3)) ;

author: Cloud-CV

repo: object-proposals


text:
function exit_reason = checkSimulationStopConditions(n_cells, Phis_t, cs_bar_t, T_t, param)
% checkSimulationStopConditions checks if a particular stop condition for the simulation is met

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

exit_reason = 0;
% Check stop conditions for each cell
for i=1:n_cells
    Phis_pos_cc_t = 1.5*Phis_t{i}(end,1) - 0.5*Phis_t{i}(end,2);
    Phis_neg_cc_t = 1.5*Phis_t{i}(end,end) - 0.5*Phis_t{i}(end,end-1);
    voltage = Phis_pos_cc_t - Phis_neg_cc_t;
    Sout    = internalSOCestimate(cs_bar_t,param,i);
    max_layer_temperature = max(T_t{1}(end,:));
    % Break conditions.
    if(voltage<param{i}.CutoffVoltage)
        if param{i}.suppress_status_prints == 0
            fprintf('\nCell #%d  below its Cutoff voltage. Stopping ...\n',i);
        end
        if exit_reason==0
            exit_reason=[];
        end
        exit_reason = [exit_reason;1];
    end

    if(voltage>param{i}.CutoverVoltage)
        if param{i}.suppress_status_prints == 0
            fprintf('\nCell #%d  above its Cutover voltage. Stopping ...\n',i);
        end
        if exit_reason==0
            exit_reason=[];
        end
        exit_reason = [exit_reason;2];
    end

    if(Sout<param{i}.CutoffSOC)
        if param{i}.suppress_status_prints == 0
            fprintf('\nCell #%d  below its Cutoff SOC. Stopping ...\n',i);
        end
        if exit_reason==0
            exit_reason=[];
        end
        exit_reason = [exit_reason;3];
    end

    if(Sout>param{i}.CutoverSOC)
        if param{i}.suppress_status_prints == 0
            fprintf('\nCell #%d  above its Cutover SOC. Stopping ...\n',i);
        end
        if exit_reason==0
            exit_reason=[];
        end
        exit_reason = [exit_reason;4];
    end

    if(max_layer_temperature > 0.99*param{i}.Tmax)
        if param{i}.suppress_status_prints == 0
            fprintf('\nCell #%d  above its Maximum Permitted Temperature. Stopping ...\n',i);
        end
        if exit_reason==0
            exit_reason=[];
        end
        exit_reason = [exit_reason;5];
    end

    if(param{i}.enable_csneg_Saturation_limit == 1)
        if(param{i}.SolidPhaseDiffusion~=3)
            max_cs_surface_neg = max(cs_bar_t{i}(end,param{i}.Np+1:end));
        else
            max_cs_surface_neg = max(cs_bar_t{i}(end, (param{i}.Np*param{i}.Nr_p) +1:end));
        end

        if (max_cs_surface_neg > param{i}.cs_sat_thresh*param{i}.cs_neg_saturation)
            if param{i}.suppress_status_prints == 0
                fprintf('\n\nCell #%d  above its safe surface concentration (saturation threshold) Stopping ...\n\n',i);
            end
            if exit_reason==0
                exit_reason=[];
            end
            exit_reason = [exit_reason;6];
        end
    end
end
end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
function [lin_restr,nonlin_restr,markov_chains]=create_restrictions_and_markov_chains2(markov_chains)
% create_restrictions_and_markov_chains2 -- creates restrictions and
% markov chains for the SVAR model in which Coefficients and variances have
% different chains, different regimes, and different durations
%
% ::
%
%
%   [lin_restr,nonlin_restr,tpl]=create_restrictions_and_markov_chains2(tpl)
%
% Args:
%
%    - **markov_chains** [empty|struct]: structure of previously defined
%    markov chains
%
% Returns:
%    :
%
%    - **lin_restr** [cell]: cell array of restrictions (see below).
%
%    - **nonlin_restr** [cell]: cell array of inequality restrictions
%
%    - **markov_chains** [struct]: modified markov chains
%
% Note:
%
%    - The syntax to construct a restriction
%      --> ai(eqtn)
%      --> ai(eqtn,vbl)
%      --> ai(eqtn,vbl,chain_name,state)
%      --> a(eqtn)
%      --> a(eqtn,vbl)
%      --> a(eqtn,vbl,chain_name,state)
%      - **eqtn** [integer]: integer
%      - **vbl** [integer|char]: integer or variable name
%      - **i** [integer]: lag
%      - **chain_name** [char]: name of the markov chain
%      - **state** [integer]: state number
%
%    - The lag coefficients are labelled a0, a1, a2,...,ak, for a model with k
%    lags. Obviously, a0 denotes the contemporaneous coefficients.
%
%    - The constant terms labelled c_1_1, c_2_2,...,c_n_n, for a model with n
%    endogenous variables.
%
%    - The standard deviations labelled s_1_1, s_2_2,...,s_n_n, for a
%    model with n endogenous variables.
%
% Example:
%
%    See also:

if nargin==0||isempty(markov_chains)

    markov_chains=struct('name',{},...
    'states_expected_duration',{},...
    'controlled_parameters',{});

end

% We borrow both the restrictions and the markov chains from the model in
% which all coefficients across all equations switch in lockstep.
%--------------------------------------------------------------------------
[lin_restr,nonlin_restr,markov_chains]=create_restrictions_and_markov_chains1(markov_chains);

% Then we add another chain controling all variances across all equations
%-------------------------------------------------------------------------
% We just make sure we do not change the restrictions
last=numel(markov_chains);

markov_chains(last+1)=struct('name','syncvol',...
    'states_expected_duration',[2+1i,2+1i,2+1i],...
    'controlled_parameters',{{'s'}});

end

author: jmaih


text:
function [M, Minit, output] = cp_apr(X, R, varargin)
%CP_APR Compute nonnegative CP with alternating Poisson regression.
%
%   M = CP_APR(X, R) computes an estimate of the best rank-R CP model of a
%   nonnegative tensor X using an alternating Poisson regression. This is
%   most appropriate for sparse count data (i.e., nonnegative integer
%   values) because it uses Kullback-Liebler divergence.  The input X can
%   be a tensor or sptensor. The result M is a ktensor.  Input data must be
%   nonnegative, and the computed ktensor factors are all nonnegative.
%
%   Different algorithm variants are available (selected by the 'alg'
%   parameter):
%     'pqnr' - row subproblems by projected quasi-Newton (default)
%     'pdnr' - row subproblems by projected damped Hessian
%     'mu'   - multiplicative update (default in version 2.5)
%
%   M = CP_APR(X, R, 'param', value, ...) specifies optional parameters and
%   values. Some parameters work in all situations, others apply only for
%   a particular choice of algorithm.
%
%   Valid parameters and their default values are:
%      'alg'           - Algorithm ['mu'|'pdnr'|'pqnr'] {'pqnr'}
%      'stoptol'       - Tolerance on the overall KKT violation {1.0e-4}
%      'stoptime'      - Maximum number of seconds to run {1e6}
%      'maxiters'      - Maximum number of iterations {1000}
%      'init'          - Initial guess [{'random'}|ktensor]
%      'maxinneriters' - Maximum inner iterations per outer iteration {10}
%      'epsDivZero'    - Safeguard against divide by zero {1.0e-10}
%      'printitn'      - Print every n outer iterations; 0 for none {1}
%      'printinneritn' - Print every n inner iterations {0}
%
%   Additional input parameters for algorithm 'mu':
%      'kappa'         - Offset to fix complementary slackness {100}
%      'kappatol'      - Tolerance on complementary slackness {1.0e-10}
%
%   Additional input parameters for algorithm 'pdnr':
%      'epsActive'     - Bertsekas tolerance for active set {1.0e-8}
%      'mu0'           - Initial damping parameter {1.0e-5}
%      'precompinds'   - Precompute sparse tensor indices {true}
%      'inexact'       - Compute inexact Newton steps {true}
%
%   Additional input parameters for algorithm 'pqnr':
%      'epsActive'     - Bertsekas tolerance for active set {1.0e-8}
%      'lbfgsMem'      - Number vector pairs to store for L-BFGS {3}
%      'precompinds'   - Precompute sparse tensor indices {true}
%
%   [M,M0] = CP_APR(...) also returns the initial guess.
%
%   [M,M0,out] = CP_APR(...) also returns additional output:
%      out.kktViolations - maximum KKT violation per iteration
%      out.nInnerIters   - number of inner iterations per outer iteration
%      out.obj           - final negative log-likelihood objective
%      out.ttlTime       - time algorithm took to converge or reach max
%      out.times         - cumulative time through each outer iteration
%    If algorithm is 'mu':
%      out.nViolations   - number of factor matrices needing complementary
%                          slackness adjustment per iteration
%    If algorithm is 'pdnr' or 'pqnr':
%      out.nZeros        - number of zero factor entries per iteration
%
%   REFERENCES:
%   * E. C. Chi and T. G. Kolda. On Tensors, Sparsity, and Nonnegative
%     Factorizations, SIAM J. Matrix Analysis,  33(4):1272-1299, Dec. 2012,
%     http://dx.doi.org/10.1137/110859063
%   * S. Hansen, T. Plantenga and T. G. Kolda, Newton-Based Optimization
%     for Kullback-Leibler Nonnegative Tensor Factorizations,
%     Optimization Methods and Software, 2015,
%     http://dx.doi.org/10.1080/10556788.2015.1009977
%
%   See also CP_ALS, KTENSOR, TENSOR, SPTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Set the algorithm choice and initial guess from input or defaults.
params = inputParser;
params.addParameter('alg', 'pqnr', @(x) (ismember(x,{'mu','pdnr','pqnr'})) );
params.addParameter('init','random', @(x) (isa(x,'ktensor') || ismember(x,{'random'})) );
params.KeepUnmatched = true;
params.parse(varargin{:});

alg   = params.Results.alg;
Minit = params.Results.init;

% Extract the number of modes in tensor X.
N = ndims(X);

if (R <= 0)
    error('Number of components requested must be positive');
end

%% Check that the data is nonnegative.
tmp = find(X < 0.0);
if (size(tmp,1) > 0)
    error('Data tensor must be nonnegative for Poisson-based factorization');
end

%% Set up an initial guess for the factor matrices.
if isa(Minit,'ktensor')
    % User provided an initial ktensor; validate it.

    if (ndims(Minit) ~= N)
        error('Initial guess does not have the right number of modes');
    end
    if (ncomponents(Minit) ~= R)
        error('Initial guess does not have the right number of components');
    end

    for n = 1:N
        if (size(Minit,n) ~= size(X,n))
            error('Mode %d of the initial guess is the wrong size',n);
        end
        if (min(min(Minit.U{n})) < 0.0)
            error('Initial guess has negative element in mode %d',n);
        end
    end
    if (min(Minit.lambda) < 0.0)
        error('Initial guess has a negative ktensor weight');
    end

elseif strcmp(Minit,'random')
    % Choose random values for each element in the range (0,1).
    F = cell(N,1);
    for n = 1:N
        F{n} = rand(size(X,n),R);
    end
    Minit = ktensor(F);
end

%% Call a solver based on the choice of algorithm parameter, passing
%  all the other input parameters.
if strcmp(alg,'mu')
    [M, output] = tt_cp_apr_mu (X, R, Minit, params.Unmatched);
    output.params.alg = 'mu';

elseif strcmp(alg,'pdnr')
    [M, output] = tt_cp_apr_pdnr (X, R, Minit, params.Unmatched);
    output.params.alg = 'pdnr';

elseif strcmp(alg,'pqnr')
    [M, output] = tt_cp_apr_pqnr (X, R, Minit, params.Unmatched);
    output.params.alg = 'pqnr';
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Main algorithm PQNR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [M, out] = tt_cp_apr_pqnr(X, R, Minit, varargin)
%TT_CP_APR_PQNR Compute nonnegative CP with alternating Poisson regression.
%
%   tt_cp_apr_pqnr(X, R, ...) computes an estimate of the best rank-R
%   CP model of a tensor X using an alternating Poisson regression.
%   The algorithm solves "row subproblems" in each alternating subproblem,
%   using a quasi-Newton Hessian approximation.
%   The function is typically called by cp_apr.
%
%   The model is solved by nonlinear optimization, and the code literally
%   minimizes the negative of log-likelihood.  However, printouts to the
%   console reverse the sign to show maximization of log-likelihood.
%
%   The function call can specify optional parameters and values.
%   Valid parameters and their default values are:
%      'stoptol'       - Tolerance on the overall KKT violation {1.0e-4}
%      'stoptime'      - Maximum number of seconds to run {1e6}
%      'maxiters'      - Maximum number of iterations {1000}
%      'maxinneriters' - Maximum inner iterations per outer iteration {10}
%      'epsDivZero'    - Safeguard against divide by zero {1.0e-10}
%      'printitn'      - Print every n outer iterations; 0 for no printing {1}
%      'printinneritn' - Print every n inner iterations {0}
%      'epsActive'     - Bertsekas tolerance for active set {1.0e-8}
%      'lbfgsMem'      - Number vector pairs to store for L-BFGS {3}
%      'precompinds'   - Precompute sparse tensor indices to run faster {true}
%
%   Return values are:
%      M                 - ktensor model with R components
%      out.fnEvals       - number of row obj fn evaluations per outer iteration
%      out.kktViolations - maximum KKT violation per iteration
%      out.nInnerIters   - number of inner iterations per outer iteration
%      out.nZeros        - number of factor elements equal to zero per iteration
%      out.obj           - final log-likelihood objective
%                          (minimization objective is actually -1 times this)
%      out.ttlTime       - time algorithm took to converge or reach max
%      out.times         - cumulative time through each outer iteration
%
%   REFERENCE: Samantha Hansen, Todd Plantenga, Tamara G. Kolda.
%   Newton-Based Optimization for Nonnegative Tensor Factorizations,
%   arXiv:1304.4964 [math.NA], April 2013,
%   URL: http://arxiv.org/abs/1304.4964. Submitted for publication.
%
%   See also CP_APR, KTENSOR, TENSOR, SPTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Set algorithm parameters from input or by using defaults.
params = inputParser;
params.addParamValue('epsActive', 1e-8, @isscalar);
params.addParamValue('epsDivZero',1e-10,@isscalar);
params.addParamValue('lbfgsMem',3,@isscalar);
params.addParamValue('maxinneriters',10,@isscalar);
params.addParamValue('maxiters',1000,@(x) isscalar(x) & x > 0);
params.addParamValue('precompinds',true,@(x) isa(x,'logical'));
params.addParamValue('printinneritn',0,@isscalar);
params.addParamValue('printitn',1,@isscalar);
params.addParamValue('stoptime',1e6,@isscalar);
params.addParamValue('stoptol',1e-4,@isscalar);
params.parse(varargin{:});

%% Copy from params object.
epsActSet               = params.Results.epsActive;
epsDivZero              = params.Results.epsDivZero;
nSizeLBFGS              = params.Results.lbfgsMem;
maxInnerIters           = params.Results.maxinneriters;
maxOuterIters           = params.Results.maxiters;
precomputeSparseIndices = params.Results.precompinds;
printInnerItn           = params.Results.printinneritn;
printOuterItn           = params.Results.printitn;
stoptime                = params.Results.stoptime;
stoptol                 = params.Results.stoptol;

% Extract the number of modes in tensor X.
N = ndims(X);

% If the initial guess has any rows of all zero elements, then modify
% so the row subproblem is not taking log(0).  Values will be restored to
% zero later if the unfolded X for the row has no nonzeros.
for n = 1:N
  rowsum = sum(Minit{n},2);
  tmpIx = find(rowsum == 0);
  if (isempty(tmpIx) == false)
    Minit{n}(tmpIx,1) = 1.0e-8;
  end
end

% Start with the initial guess, normalized using the vector L1 norm.
M = normalize(Minit,[],1);

% Sparse tensor flag affects how Pi and Phi are computed.
if isa(X,'sptensor')
    isSparse = true;
else
    isSparse = false;
end

% Initialize output arrays.
fnEvals = zeros(maxOuterIters,1);
kktViolations = -ones(maxOuterIters,1);
nInnerIters = zeros(maxOuterIters,1);
nzeros = zeros(maxOuterIters,1);
times = zeros(maxOuterIters,1);

if (printOuterItn > 0)
    fprintf('\nCP_PQNR (alternating Poisson regression using quasi-Newton)\n');
end
dispLineWarn = (printInnerItn > 0);

% Start the wall clock timer.
tic;

if (isSparse && precomputeSparseIndices)
    % Precompute sparse index sets for all the row subproblems.
    % Takes more memory but can cut execution time significantly in some cases.
    if (printOuterItn > 0)
        fprintf('  Precomputing sparse index sets...');
    end
    sparseIx = cell(N);
    for n = 1:N
        num_rows = size(M{n},1);
        sparseIx{n} = cell(num_rows,1);
        for jj = 1:num_rows
            sparseIx{n}{jj} = find(X.subs(:,n) == jj);
        end
    end
    if (printOuterItn > 0)
        fprintf('done\n');
    end
end

%% Main Loop: Iterate until convergence or a max threshold is reached.
for iter = 1:maxOuterIters

    isConverged = true;
    kktModeViolations = zeros(N,1);
    countInnerIters = zeros(1,N);

    % Alternate thru each factor matrix, A_1, A_2, ... , A_N.
    for n = 1:N

        % Shift the weight from lambda to mode n.
        M = redistribute(M,n);

        % Calculate Khatri-Rhao product of all matrices but the n-th.
        if (isSparse == false)
            % Data is not a sparse tensor.
            Pi = tt_calcpi_prowsubprob(X, isSparse, M, R, n, N, []);
            X_mat = double(tenmat(X,n));
        end

        num_rows = size(M{n},1);
        isRowNOTconverged = zeros(1,num_rows);

        % Loop over the row subproblems in mode n.
        for jj = 1:num_rows

            % Get data values for row jj of matricized mode n.
            if (isSparse)
                % Data is a sparse tensor.
                if (precomputeSparseIndices == false)
                    sparse_indices = find(X.subs(:,n) == jj);
                else
                    sparse_indices = sparseIx{n}{jj};
                end
                if (isempty(sparse_indices))
                    % The row jj of matricized tensor X in mode n is empty.
                    M{n}(jj,:) = 0;
                    continue
                end
                x_row = X.vals(sparse_indices);

                % Calculate just the columns of Pi needed for this row.
                Pi = tt_calcpi_prowsubprob(X, isSparse, M, ...
                                           R, n, N, sparse_indices);
            else
                x_row = X_mat(jj,:);
            end

            % Get current values of the row subproblem variables.
            m_row = M{n}(jj,:);

            % Initialize L-BFGS storage for the row subproblem.
            delm = zeros(R, nSizeLBFGS);
            delg = zeros(R, nSizeLBFGS);
            rho = zeros(nSizeLBFGS, 1);
            lbfgsPos = 1;
            m_rowOLD = [];
            gradOLD = [];

            % Iteratively solve the row subproblem with projected qNewton steps.
            for i = 1:maxInnerIters
                % Calculate the gradient.
                [gradM, phi_row] = calc_grad(isSparse, Pi, epsDivZero, ...
                                             x_row, m_row);

                if (i == 1)
                    % Original cp_aprPQN_row code (and plb_row) does a gradient
                    % step to prime the L-BFGS approximation.  However, it means
                    % a row subproblem that already converged wastes time
                    % doing a gradient step before checking KKT conditions.
                    % TODO: fix in a future release.
                    m_rowOLD = m_row;
                    gradOLD = gradM;
                    [m_row, f, f_unit, f_new, num_evals] ...
                        = tt_linesearch_prowsubprob(-gradM', gradM', ...
                                                    m_rowOLD, ...
                                                    1, 1/2, 10, 1.0e-4, ...
                                                    isSparse, x_row, Pi, ...
                                                    phi_row, dispLineWarn);
                    fnEvals(iter) = fnEvals(iter) + num_evals;
                    [gradM, phi_row] = calc_grad(isSparse, Pi, epsDivZero, ...
                                                 x_row, m_row);
                end

                % Compute the row subproblem kkt_violation.
                % Experiments in the original paper used this:
                %kkt_violation = norm(abs(min(m_row,gradM')),2);
                % Now we use | KKT |_inf:
                kkt_violation = max(abs(min(m_row,gradM')));

                % Report largest row subproblem initial violation.
                if ((i == 1) && (kkt_violation > kktModeViolations(n)))
                     kktModeViolations(n) = kkt_violation;
                end

                if (mod(i, printInnerItn) == 0)
                    fprintf('    Mode = %1d, Row = %d, InnerIt = %d', ...
                            n, jj, i);
                    if (i == 1)
                        fprintf(', RowKKT = %.2e\n', kkt_violation);
                    else
                        fprintf(', RowKKT = %.2e, RowObj = %.4e\n', ...
                                kkt_violation, -f_new);
                    end
                end

                % Check for row subproblem convergence.
                if (kkt_violation < stoptol)
                    break;
                else
                    % Not converged, so m_row will be modified.
                    isRowNOTconverged(jj) = 1;
                end

                % Update the L-BFGS approximation.
                tmp_delm = m_row - m_rowOLD;
                tmp_delg = gradM - gradOLD;
                tmp_rho = 1 / (tmp_delm * tmp_delg);
                if ((tmp_rho > 0.0) && (isinf(tmp_rho) == false))
                    delm(:,lbfgsPos) = tmp_delm;
                    delg(:,lbfgsPos) = tmp_delg;
                    rho(lbfgsPos) = tmp_rho;
                else
                    % Rho is required to be positive; if not, then skip
                    % the L-BFGS update pair.  The recommended safeguard for
                    % full BFGS is Powell damping, but not clear how to damp
                    % in 2-loop L-BFGS.
                    if (dispLineWarn)
                        fprintf('WARNING: skipping L-BFGS update, rho would be 1 / %.2e\n', ...
                                (tmp_delm * tmp_delg));
                    end
                    % Roll back lbfgsPos since it will increment later.
                    if (lbfgsPos == 1)
                        if (rho(nSizeLBFGS) > 0)
                            lbfgsPos = nSizeLBFGS;
                        else
                            % Fatal error, should not happen.
                            fprintf('ERROR: L-BFGS first iterate is bad\n');
                            return;
                        end
                    else
                        lbfgsPos = lbfgsPos - 1;
                    end
                end

                % Calculate the search direction.
                search_dir = getSearchDirPqnr(m_row, gradM, epsActSet, ...
                                              delm, delg, rho, lbfgsPos, ...
                                              i, dispLineWarn);
                lbfgsPos = mod(lbfgsPos, nSizeLBFGS) + 1;

                m_rowOLD = m_row;
                gradOLD = gradM;

                % Perform a projected linesearch and update variables.
                % Start from a unit step length, decrease by 1/2, stop with
                % sufficient decrease of 1.0e-4 or at most 10 steps.
                [m_row, f, f_unit, f_new, num_evals] ...
                    = tt_linesearch_prowsubprob(search_dir', gradOLD', m_rowOLD, ...
                                                1, 1/2, 10, 1.0e-4, ...
                                                isSparse, x_row, Pi, ...
                                                phi_row, dispLineWarn);
                fnEvals(iter) = fnEvals(iter) + num_evals;
            end

            M{n}(jj,:) = m_row;
            countInnerIters(n) = countInnerIters(n) + i;

        end

        % Test if all row subproblems have converged, which means that
        % no variables in this mode were changed.
        if (sum(isRowNOTconverged) ~= 0)
            isConverged = false;
        end

        % Shift weight from mode n back to lambda.
        M = normalize(M,[],1,n);

        % Total number of inner iterations for a given outer iteration,
        % totalled across all modes and all row subproblems in each mode.
        nInnerIters(iter) = nInnerIters(iter) + countInnerIters(n);
    end

    % Save output items for the outer iteration.
    num_zero = 0;
    for n = 1:N
        num_zero = num_zero + nnz(find(M{n} == 0.0));
    end
    nzeros(iter) = num_zero;
    kktViolations(iter) = max(kktModeViolations);

    % Print outer iteration status.
    if (mod(iter,printOuterItn) == 0)
        fprintf('%4d. Ttl Inner Its: %d, KKT viol = %.2e, obj = %.8e, nz: %d\n', ...
        iter, nInnerIters(iter), kktViolations(iter), tt_loglikelihood(X,M), ...
        num_zero);
    end

    times(iter) = toc;

    % Check for convergence
    if (isConverged)
        break;
    end
    if (times(iter) > stoptime)
        fprintf('Exiting because time limit exceeded\n');
        break;
    end

end

t_stop = toc;

%% Clean up final result and set output items.
M = normalize(M,'sort',1);
loglike = tt_loglikelihood(X,M);

if (printOuterItn > 0)
    % For legacy reasons, compute "fit", the fraction explained by the model.
    % Fit is in the range [0,1], with 1 being the best fit.
    normX = norm(X);
    normresidual = sqrt( normX^2 + norm(M)^2 - 2 * innerprod(X,M) );
    fit = 1 - (normresidual / normX);

    fprintf('===========================================\n');
    fprintf(' Final log-likelihood = %e \n', loglike);
    fprintf(' Final least squares fit = %e \n', fit);
    fprintf(' Final KKT violation = %7.7e\n', kktViolations(iter));
    fprintf(' Total inner iterations = %d\n', sum(nInnerIters));
    fprintf(' Total execution time = %.2f secs\n', t_stop);
end

out = struct;
out.params = params.Results;
out.obj = loglike;
out.kktViolations = kktViolations(1:iter);
out.fnEvals = fnEvals(1:iter);
out.nInnerIters = nInnerIters(1:iter);
out.nZeros = nzeros(1:iter);
out.times = times(1:iter);
out.ttlTime = t_stop;

end

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

function [grad_row, phi_row] = calc_grad(isSparse, Pi, eps_div_zero, x_row, m_row)
%function grad_row = calc_grad(isSparse, Pi, eps_div_zero, x_row, m_row)
% Compute the gradient for a PQNR row subproblem.
%
%   isSparse     - true if x_row is sparse, false if dense
%   Pi           - matrix
%   eps_div_zero - safeguard value to prevent division by zero
%   x_row        - row vector of data values for the row subproblem
%   m_row        - vector of variables for the row subproblem
%
%   Returns the gradient vector for a row subproblem.

    if (isSparse)
        v = m_row * Pi';
        w = x_row' ./ max(v, eps_div_zero);
        phi_row = w * Pi;

    else
        v = m_row * Pi';
        w = x_row ./ max(v, eps_div_zero);
        phi_row = w * Pi;
    end

    grad_row = (ones(size(phi_row)) - phi_row)';
end

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

function [d] = getSearchDirPqnr (m_row, grad, epsActSet, ...
                                 delta_m, delta_g, rho, lbfgs_pos, ...
                                 iters, disp_warn)
% Compute the search direction by projecting with L-BFGS.
%
%   m_row     - current variable values
%   grad      - gradient at m_row
%   epsActSet - Bertsekas tolerance for active set determination
%   delta_m   - L-BFGS array of vector variable deltas
%   delta_g   - L-BFGS array of gradient deltas
%   lbfgs_pos - pointer into L-BFGS arrays
%
%   Returns
%     d       - search direction based on current L-BFGS and grad
%
%   Adapted from MATLAB code of Dongmin Kim and Suvrit Sra written in 2008.
%   Modified extensively to solve row subproblems and use a better linesearch;
%   see the reference at the top of this file for details.

    lbfgsSize = size(delta_m,2);

    % Determine active and free variables.
    % If epsActSet is zero, then the following works:
    %   fixedVars = find((m_row == 0) & (grad' > 0));
    % For the general case this works but is less clear and assumes m_row > 0:
    %   fixedVars = find((grad' > 0) & (m_row <= min(epsActSet,grad')));
    projGradStep = (m_row - grad') .* (m_row - grad' > 0);
    wk = norm(m_row - projGradStep);
    fixedVars = find((grad' > 0) & (m_row <= min(epsActSet,wk)));

    d = -grad;
    d(fixedVars) = 0;

    if ((delta_m(:,lbfgs_pos)' * delta_g(:,lbfgs_pos)) == 0.0)
        % Cannot proceed with this L-BFGS data; most likely the iteration
        % has converged, so this is rarely seen.
        if (disp_warn)
            fprintf('WARNING: L-BFGS update is orthogonal, using gradient\n');
        end
        return;
    end

    alpha = ones(lbfgsSize,1);
    k = lbfgs_pos;

    % Perform an L-BFGS two-loop recursion to compute the search direction.

    for i = 1 : min(iters, lbfgsSize)
        alpha(k) = rho(k) * delta_m(:, k)' * d;
        d = d - alpha(k) * delta_g(:, k);
        k = lbfgsSize - mod(1 - k, lbfgsSize);
    end

    coef = 1 / rho(lbfgs_pos) / (delta_g(:, lbfgs_pos)' * delta_g(:, lbfgs_pos));
    d = coef * d;

    for i = 1 : min(iters, lbfgsSize)
        k = mod(k, lbfgsSize) + 1;
        b = rho(k) * delta_g(:, k)' * d;
        d = d + (alpha(k) - b) * delta_m(:, k);
    end

    d(fixedVars) = 0;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Main algorithm PDNR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [M, out] = tt_cp_apr_pdnr(X, R, Minit, varargin)
%TT_CP_APR_PDNR Compute nonnegative CP with alternating Poisson regression.
%
%   tt_cp_apr_pdnr(X, R, ...) computes an estimate of the best rank-R
%   CP model of a tensor X using an alternating Poisson regression.
%   The algorithm solves "row subproblems" in each alternating subproblem,
%   using a Hessian of size R^2.
%   The function is typically called by cp_apr.
%
%   The model is solved by nonlinear optimization, and the code literally
%   minimizes the negative of log-likelihood.  However, printouts to the
%   console reverse the sign to show maximization of log-likelihood.
%
%   The function call can specify optional parameters and values.
%   Valid parameters and their default values are:
%      'stoptol'       - Tolerance on the overall KKT violation {1.0e-4}
%      'stoptime'      - Maximum number of seconds to run {1e6}
%      'maxiters'      - Maximum number of iterations {1000}
%      'maxinneriters' - Maximum inner iterations per outer iteration {10}
%      'epsDivZero'    - Safeguard against divide by zero {1.0e-10}
%      'printitn'      - Print every n outer iterations; 0 for no printing {1}
%      'printinneritn' - Print every n inner iterations {0}
%      'epsActive'     - Bertsekas tolerance for active set {1.0e-8}
%      'mu0'           - Initial damping parameter {1.0e-5}
%      'precompinds'   - Precompute sparse tensor indices to run faster {true}
%      'inexact'       - Compute inexact Newton steps {true}
%
%   Return values are:
%      M                 - ktensor model with R components
%      out.fnEvals       - number of row obj fn evaluations per outer iteration
%      out.kktViolations - maximum KKT violation per iteration
%      out.nInnerIters   - number of inner iterations per outer iteration
%      out.nZeros        - number of factor elements equal to zero per iteration
%      out.obj           - final log-likelihood objective
%                          (minimization objective is actually -1 times this)
%      out.ttlTime       - time algorithm took to converge or reach max
%      out.times         - cumulative time through each outer iteration
%
%   REFERENCE: Samantha Hansen, Todd Plantenga, Tamara G. Kolda.
%   Newton-Based Optimization for Nonnegative Tensor Factorizations,
%   arXiv:1304.4964 [math.NA], April 2013,
%   URL: http://arxiv.org/abs/1304.4964. Submitted for publication.
%
%   See also CP_APR, KTENSOR, TENSOR, SPTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Set algorithm parameters from input or by using defaults.
params = inputParser;
params.addParamValue('epsActive', 1e-8, @isscalar);
params.addParamValue('epsDivZero',1e-10,@isscalar);
params.addParamValue('maxinneriters',10,@isscalar);
params.addParamValue('maxiters',1000,@(x) isscalar(x) & x > 0);
params.addParamValue('precompinds',true,@(x) isa(x,'logical'));
params.addParamValue('inexact',true,@(x) isa(x,'logical'));
params.addParamValue('mu0',1e-5,@isscalar);
params.addParamValue('printinneritn',0,@isscalar);
params.addParamValue('printitn',1,@isscalar);
params.addParamValue('stoptime',1e6,@isscalar);
params.addParamValue('stoptol',1e-4,@isscalar);
params.parse(varargin{:});

%% Copy from params object.
epsActSet               = params.Results.epsActive;
epsDivZero              = params.Results.epsDivZero;
maxInnerIters           = params.Results.maxinneriters;
maxOuterIters           = params.Results.maxiters;
mu0                     = params.Results.mu0;
precomputeSparseIndices = params.Results.precompinds;
inexactNewton           = params.Results.inexact;
printInnerItn           = params.Results.printinneritn;
printOuterItn           = params.Results.printitn;
stoptime                = params.Results.stoptime;
stoptol                 = params.Results.stoptol;

% Extract the number of modes in tensor X.
N = ndims(X);

% If the initial guess has any rows of all zero elements, then modify
% so the row subproblem is not taking log(0).  Values will be restored to
% zero later if the unfolded X for the row has no nonzeros.
for n = 1:N
  rowsum = sum(Minit{n},2);
  tmpIx = find(rowsum == 0);
  if (isempty(tmpIx) == false)
    Minit{n}(tmpIx,1) = 1.0e-8;
  end
end

% Start with the initial guess, normalized using the vector L1 norm.
M = normalize(Minit,[],1);

% Sparse tensor flag affects how Pi and Phi are computed.
if isa(X,'sptensor')
    isSparse = true;
else
    isSparse = false;
end

% Initialize output arrays.
fnEvals = zeros(maxOuterIters,1);
kktViolations = -ones(maxOuterIters,1);
nInnerIters = zeros(maxOuterIters,1);
nzeros = zeros(maxOuterIters,1);
times = zeros(maxOuterIters,1);

if (printOuterItn > 0)
    fprintf('\nCP_PDNR (alternating Poisson regression using damped Newton)\n');
end
dispLineWarn = (printInnerItn > 0);

% Start the wall clock timer.
tic;

if (isSparse && precomputeSparseIndices)
    % Precompute sparse index sets for all the row subproblems.
    % Takes more memory but can cut execution time significantly in some cases.
    if (printOuterItn > 0)
        fprintf('  Precomputing sparse index sets...');
    end
    sparseIx = cell(N);
    for n = 1:N
        num_rows = size(M{n},1);
        sparseIx{n} = cell(num_rows,1);
        for jj = 1:num_rows
            sparseIx{n}{jj} = find(X.subs(:,n) == jj);
        end
    end
    if (printOuterItn > 0)
        fprintf('done\n');
    end
end

e_vec = ones(1,R);

rowsubprobStopTol = stoptol;

%% Main Loop: Iterate until convergence or a max threshold is reached.
for iter = 1:maxOuterIters

    isConverged = true;
    kktModeViolations = zeros(N,1);
    countInnerIters = zeros(1,N);

    % Alternate thru each factor matrix, A_1, A_2, ... , A_N.
    for n = 1:N

        % Shift the weight from lambda to mode n.
        M = redistribute(M,n);

        % Calculate Khatri-Rhao product of all matrices but the n-th.
        if (isSparse == false)
            % Data is not a sparse tensor.
            Pi = tt_calcpi_prowsubprob(X, isSparse, M, R, n, N, []);
            X_mat = double(tenmat(X,n));
        end

        num_rows = size(M{n},1);
        isRowNOTconverged = zeros(1,num_rows);

        % Loop over the row subproblems in mode n.
        for jj = 1:num_rows
            % Initialize the damped Hessian parameter for the row subproblem.
            mu = mu0;

            % Get data values for row jj of matricized mode n.
            if (isSparse)
                % Data is a sparse tensor.
                if (precomputeSparseIndices == false)
                    sparse_indices = find(X.subs(:,n) == jj);
                else
                    sparse_indices = sparseIx{n}{jj};
                end
                if (isempty(sparse_indices))
                    % The row jj of matricized tensor X in mode n is empty.
                    M{n}(jj,:) = 0;
                    continue
                end
                x_row = X.vals(sparse_indices);

                % Calculate just the columns of Pi needed for this row.
                Pi = tt_calcpi_prowsubprob(X, isSparse, M, ...
                                           R, n, N, sparse_indices);
            else
                x_row = X_mat(jj,:);
            end

            % Get current values of the row subproblem variables.
            m_row = M{n}(jj,:);

            % Iteratively solve the row subproblem with projected Newton steps.
            innerIterMaximum = maxInnerIters;
            if (inexactNewton && (iter == 1))
                innerIterMaximum = 2;
            end
            for i = 1:innerIterMaximum
                % Calculate the gradient.
                [phi_row, ups_row] ...
                    = calc_partials(isSparse, Pi, epsDivZero, x_row, m_row);
                gradM = (e_vec - phi_row)';

                % Compute the row subproblem kkt_violation.
                % Experiments in the original paper used this:
                %kkt_violation = norm(abs(min(m_row,gradM')),2);
                % Now we use | KKT |_inf:
                kkt_violation = max(abs(min(m_row,gradM')));

                % Report largest row subproblem initial violation.
                if ((i == 1) && (kkt_violation > kktModeViolations(n)))
                     kktModeViolations(n) = kkt_violation;
                end

                if (mod(i, printInnerItn) == 0)
                    fprintf('    Mode = %1d, Row = %d, InnerIt = %d', ...
                            n, jj, i);
                    if (i == 1)
                        fprintf(', RowKKT = %.2e\n', kkt_violation);
                    else
                        fprintf(', RowKKT = %.2e, RowObj = %.4e\n', ...
                                kkt_violation, -f_new);
                    end
                end

                % Check for row subproblem convergence.
                if (kkt_violation < rowsubprobStopTol)
                    break;
                else
                    % Not converged, so m_row will be modified.
                    isRowNOTconverged(jj) = 1;
                end

                % Calculate the search direction.
                [search_dir, predicted_red] ...
                    = getSearchDirPdnr(Pi, ups_row, R, gradM, m_row, mu, epsActSet);

                % Perform a projected linesearch and update variables.
                % Start from a unit step length, decrease by 1/2, stop with
                % sufficient decrease of 1.0e-4 or at most 10 steps.
                [m_rowNEW, f_old, f_unit, f_new, num_evals] ...
                    = tt_linesearch_prowsubprob(search_dir', gradM', m_row, ...
                                                1, 1/2, 10, 1.0e-4, ...
                                                isSparse, x_row, Pi, ...
                                                phi_row, dispLineWarn);
                fnEvals(iter) = fnEvals(iter) + num_evals;
                m_row = m_rowNEW;

                % Update damping parameter mu based on the unit step length,
                % which is returned in f_unit.
                actual_red = f_old - f_unit;
                rho = actual_red / (-predicted_red);
                if (predicted_red == 0)
                    mu = 10 * mu;
                elseif (rho < 1/4)
                    mu = (7/2) * mu;
                elseif (rho > 3/4)
                    mu = (2/7) * mu;
                end
            end

            M{n}(jj,:) = m_row;
            countInnerIters(n) = countInnerIters(n) + i;

        end

        % Test if all row subproblems have converged, which means that
        % no variables in this mode were changed.
        if (sum(isRowNOTconverged) ~= 0)
            isConverged = false;
        end

        % Shift weight from mode n back to lambda.
        M = normalize(M,[],1,n);

        % Total number of inner iterations for a given outer iteration,
        % totalled across all modes and all row subproblems in each mode.
        nInnerIters(iter) = nInnerIters(iter) + countInnerIters(n);
    end

    % Save output items for the outer iteration.
    num_zero = 0;
    for n = 1:N
        num_zero = num_zero + nnz(find(M{n} == 0.0));
    end
    nzeros(iter) = num_zero;
    kktViolations(iter) = max(kktModeViolations);
    if (inexactNewton)
        rowsubprobStopTol = max(stoptol, kktViolations(iter) / 100.0);
    end

    % Print outer iteration status.
    if (mod(iter,printOuterItn) == 0)
        fprintf('%4d. Ttl Inner Its: %d, KKT viol = %.2e, obj = %.8e, nz: %d\n', ...
        iter, nInnerIters(iter), kktViolations(iter), tt_loglikelihood(X,M), ...
        num_zero);
    end

    times(iter) = toc;

    % Check for convergence
    if (isConverged && (inexactNewton == false))
        break;
    end
    if (isConverged && (inexactNewton == true) && (rowsubprobStopTol <= stoptol))
        break;
    end
    if (times(iter) > stoptime)
        fprintf('Exiting because time limit exceeded\n');
        break;
    end

end

t_stop = toc;

%% Clean up final result and set output items.
M = normalize(M,'sort',1);
loglike = tt_loglikelihood(X,M);

if (printOuterItn > 0)
    % For legacy reasons, compute "fit", the fraction explained by the model.
    % Fit is in the range [0,1], with 1 being the best fit.
    normX = norm(X);
    normresidual = sqrt( normX^2 + norm(M)^2 - 2 * innerprod(X,M) );
    fit = 1 - (normresidual / normX);

    fprintf('===========================================\n');
    fprintf(' Final log-likelihood = %e \n', loglike);
    fprintf(' Final least squares fit = %e \n', fit);
    fprintf(' Final KKT violation = %7.7e\n', kktViolations(iter));
    fprintf(' Total inner iterations = %d\n', sum(nInnerIters));
    fprintf(' Total execution time = %.2f secs\n', t_stop);
end

out = struct;
out.params = params.Results;
out.obj = loglike;
out.kktViolations = kktViolations(1:iter);
out.fnEvals = fnEvals(1:iter);
out.nInnerIters = nInnerIters(1:iter);
out.nZeros = nzeros(1:iter);
out.times = times(1:iter);
out.ttlTime = t_stop;

end

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

function [phi_row, ups_row] ...
    = calc_partials(isSparse, Pi, eps_div_zero, x_row, m_row)
% Compute derivative quantities for a PDNR row subproblem.
%
%   isSparse     - true if x_row is sparse, false if dense
%   Pi           - matrix
%   eps_div_zero - safeguard value to prevent division by zero
%   x_row        - row vector of data values for the row subproblem
%   m_row        - vector of variables for the row subproblem
%
%   Returns two vectors for a row subproblem:
%     phi_row - gradient of row subproblem, except for a constant
%     ups_row - intermediate quantity (upsilon) used for second derivatives

    if (isSparse)
        v = m_row * Pi';
        w = x_row' ./ max(v, eps_div_zero);
        phi_row = w * Pi;
        u = v .^ 2;
        ups_row = x_row' ./ max(u, eps_div_zero);

    else
        v = m_row * Pi';
        w = x_row ./ max(v, eps_div_zero);
        phi_row = w * Pi;
        u = v .^ 2;
        ups_row = x_row ./ max(u, eps_div_zero);
    end

end

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

function H = getHessian(upsilon, Pi, free_indices)
% Return the Hessian for one PDNR row subproblem of M{n}, for just the rows and
% columns corresponding to the free variables.

    num_free = length(free_indices);
    H = zeros(num_free,num_free);
    for i = 1:num_free
        for j = i:num_free
            c = free_indices(i);
            d = free_indices(j);
            val = sum(upsilon' .* Pi(:,c) .* Pi(:,d));
            H(i,j) = val;
            H(j,i) = val;
        end
    end

end

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

function [search_dir, pred_red] ...
    = getSearchDirPdnr (Pi, ups_row, R, gradM, m_row, mu, epsActSet)
% Compute the search direction for PDNR using a two-metric projection
% with damped Hessian.
%
%   Pi        - matrix
%   ups_row   - intermediate quantity (upsilon) used for second derivatives
%   R         - number of variables for the row subproblem
%   gradM     - gradient vector for the row subproblem
%   m_row     - vector of variables for the row subproblem
%   mu        - damping parameter
%   epsActSet - Bertsekas tolerance for active set determination
%
%   Returns:
%     search_dir - search direction vector
%     pred_red   - predicted reduction in quadratic model

    search_dir = zeros(R,1);
    projGradStep = (m_row - gradM') .* (m_row - gradM' > 0);
    wk = norm(m_row - projGradStep);

    % Determine active and free variables.
    num_free = 0;
    free_indices_tmp = zeros(R,1);
    for r = 1:R
        if ((m_row(r) <= min(epsActSet,wk)) && (gradM(r) > 0) )
            % Variable is not free (belongs to set A or G).
            if (m_row(r) ~= 0)
                % Variable moves according to the gradient (set G).
                search_dir(r) = -gradM(r);
            end
        else
            % Variable is free (set F).
            num_free = num_free + 1;
            free_indices_tmp(num_free) = r;
        end
    end
    free_indices = free_indices_tmp(1:num_free);

    % Compute the Hessian for free variables.
    Hessian_free = getHessian(ups_row, Pi, free_indices);
    grad_free = -gradM(free_indices);

    % Compute the damped Newton search direction over free variables.
    search_dir(free_indices) ...
        = linsolve(Hessian_free + (mu * eye(num_free)), grad_free);

    % If the Hessian is too ill-conditioned, use gradient descent.
    [~, msgid] = lastwarn('MATLAB:noWarning');
    if (strcmp(msgid,'MATLAB:nearlySingularMatrix'))
        fprintf('WARNING: damped Hessian is nearly singular\n');
        search_dir = -gradM;
    end

    % Calculate expected reduction in the quadratic model of the objective.
    q = search_dir(free_indices)' ...
        * (Hessian_free + (mu * eye(num_free))) ...
        * search_dir(free_indices);
    pred_red = (search_dir(free_indices)' * gradM(free_indices)) + (0.5 * q);
    if (pred_red > 0)
        fprintf('ERROR: expected decrease is positive\n');
        search_dir = -gradM;
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Main algorithm MU
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [M, output] = tt_cp_apr_mu(X, R, Minit, varargin)
%TT_CP_APR_MU Compute nonnegative CP with alternating Poisson regression.
%
%   tt_cp_apr_mu(X, R, ...) computes an estimate of the best rank-R
%   CP model of a tensor X using an alternating Poisson regression.
%   The algorithm solves each alternating subproblem using multiplicative
%   updates with adjustments for values near zero.
%   The function is typically called by cp_apr.
%
%   The function call can specify optional parameters and values.
%   Valid parameters and their default values are:
%      'stoptol'       - Tolerance on the overall KKT violation {1.0e-4}
%      'stoptime'      - Maximum number of seconds to run {1e6}
%      'maxiters'      - Maximum number of iterations {1000}
%      'maxinneriters' - Maximum inner iterations per outer iteration {10}
%      'epsDivZero'    - Safeguard against divide by zero {1.0e-10}
%      'printitn'      - Print every n outer iterations; 0 for no printing {1}
%      'printinneritn' - Print every n inner iterations {0}
%      'kappatol'      - Tolerance on complementary slackness {1.0e-10}
%      'kappa'         - Offset to fix complementary slackness {100}
%
%   Return values are:
%      M                 - ktensor model with R components
%      out.kktViolations - maximum KKT violation per iteration
%      out.nInnerIters   - number of inner iterations per outer iteration
%      out.nViolations   - number of factor matrices needing complementary
%                          slackness adjustment per iteration
%      out.obj           - final log-likelihood objective
%      out.ttlTime       - time algorithm took to converge or reach max
%      out.times         - cumulative time through each outer iteration
%
%   REFERENCE: E. C. Chi and T. G. Kolda. On Tensors, Sparsity, and
%   Nonnegative Factorizations, arXiv:1112.2414 [math.NA], December 2011,
%   URL: http://arxiv.org/abs/1112.2414. Submitted for publication.
%
%   See also CP_APR, KTENSOR, TENSOR, SPTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Set algorithm parameters from input or by using defaults.
params = inputParser;
params.addParamValue('epsDivZero',1e-10,@isscalar);
params.addParamValue('kappa',1e-2,@isscalar);
params.addParamValue('kappatol',1e-10,@isscalar);
params.addParamValue('maxinneriters',10,@isscalar);
params.addParamValue('maxiters',1000,@(x) isscalar(x) & x > 0);
params.addParamValue('printinneritn',0,@isscalar);
params.addParamValue('printitn',1,@isscalar);
params.addParamValue('stoptime',1e6,@isscalar);
params.addParamValue('stoptol',1e-4,@isscalar);
params.parse(varargin{:});

%% Extract dimensions of X and number of dimensions of X.
N = ndims(X);

%% Copy from params object.
epsilon       = params.Results.epsDivZero;
tol           = params.Results.stoptol;
stoptime      = params.Results.stoptime;
maxOuterIters = params.Results.maxiters;
kappa         = params.Results.kappa;
kappaTol      = params.Results.kappatol;
maxInnerIters = params.Results.maxinneriters;
printOuterItn = params.Results.printitn;
printInnerItn = params.Results.printinneritn;
kktViolations = -ones(maxOuterIters,1);
nInnerIters   = zeros(maxOuterIters,1);
times         = zeros(maxOuterIters,1);

%% Set up and error checking on initial guess for U.
if isa(Minit,'ktensor')
    if ndims(Minit) ~= N
        error('Initial guess does not have the right number of dimensions');
    end

    if ncomponents(Minit) ~= R
        error('Initial guess does not have the right number of components');
    end

    for n = 1:N
        if size(Minit,n) ~= size(X,n)
            error('Dimension %d of the initial guess is the wrong size',n);
        end
    end
elseif strcmp(Minit,'random')
    F = cell(N,1);
    for n = 1:N
        F{n} = rand(size(X,n),R);
    end
    Minit = ktensor(F);
else
    error('The selected initialization method is not supported');
end

%% Set up for iterations - initializing M and Phi.
M = normalize(Minit,[],1);
Phi = cell(N,1);
kktModeViolations = zeros(N,1);

if printOuterItn > 0
  fprintf('\nCP_APR:\n');
end

nViolations = zeros(maxOuterIters,1);

% Start the wall clock timer.
tic;

% PDN-R and PQN-R benefit from precomputing sparse indices of X for each
% mode subproblem.  However, MU execution time barely changes, so the
% precompute option is not offered.

%% Main Loop: Iterate until convergence.
for iter = 1:maxOuterIters

    isConverged = true;
    for n = 1:N

        % Make adjustments to entries of M{n} that are violating
        % complementary slackness conditions.
        if (iter > 1)
            V = (Phi{n} > 1) & (M{n} < kappaTol);
            if any(V(:))
                nViolations(iter) = nViolations(iter) + 1;
                M{n}(V>0) = M{n}(V>0) + kappa;
            end
        end

        % Shift the weight from lambda to mode n
        M = redistribute(M,n);

        % Calculate product of all matrices but the n-th
        % (Sparse case only calculates entries corresponding to nonzeros in X)
        Pi = calculatePi(X, M, R, n, N);

        % Do the multiplicative updates
        for i = 1:maxInnerIters

            % Count the inner iterations
            nInnerIters(iter) = nInnerIters(iter) + 1;

            % Calculate matrix for multiplicative update
            Phi{n} = calculatePhi(X, M, R, n, Pi, epsilon);

            % Check for convergence
            kktModeViolations(n) = max(abs(vectorizeForMu(min(M.U{n},1-Phi{n}))));
            if (kktModeViolations(n) < tol)
                break;
            else
                isConverged = false;
            end

            % Do the multiplicative update
            M{n} = M{n} .* Phi{n};

            % Print status
             if mod(i, printInnerItn)==0
                 fprintf('    Mode = %1d, Inner Iter = %2d, KKT violation = %.6e\n', n, i, kktModeViolations(n));
             end
        end

        % Shift weight from mode n back to lambda
        M = normalize(M,[],1,n);

    end

    kktViolations(iter) = max(kktModeViolations);

    if (mod(iter,printOuterItn)==0)
        fprintf(' Iter %4d: Inner Its = %2d KKT violation = %.6e, nViolations = %2d\n', ...
        iter, nInnerIters(iter), kktViolations(iter), nViolations(iter));
    end
    times(iter) = toc;

    % Check for convergence
    if (isConverged)
        if printOuterItn>0
            fprintf('Exiting because all subproblems reached KKT tol.\n');
        end
        break;
    end
    if (times(iter) > stoptime)
        if printOuterItn>0
            fprintf('Exiting because time limit exceeded.\n');
        end
        break;
    end
end
t_stop = toc;

%% Clean up final result
M = normalize(M,'sort',1);

obj = tt_loglikelihood(X,M);
if printOuterItn>0
    normX = norm(X);
    normresidual = sqrt( normX^2 + norm(M)^2 - 2 * innerprod(X,M) );
    fit = 1 - (normresidual / normX); %fraction explained by model
    fprintf('===========================================\n');
    fprintf(' Final log-likelihood = %e \n', obj);
    fprintf(' Final least squares fit = %e \n', fit);
    fprintf(' Final KKT violation = %7.7e\n', kktViolations(iter));
    fprintf(' Total inner iterations = %d\n', sum(nInnerIters));
    fprintf(' Total execution time = %.2f secs\n', t_stop);
end

output = struct;
output.params = params.Results;
output.kktViolations = kktViolations(1:iter);
output.nInnerIters = nInnerIters(1:iter);
output.nViolations = nViolations(1:iter);
output.nTotalIters = sum(nInnerIters);
output.times = times(1:iter);
output.ttlTime = t_stop;
output.obj = obj;

end

function Pi = calculatePi(X, M, R, n, N)

if (isa(X,'sptensor'))
    Pi = ones(nnz(X), R);
    for nn = [1:n-1,n+1:N]
        Pi = M{nn}(X.subs(:,nn),:) .* Pi;
    end
else
    U = M.U;
    Pi = khatrirao(U{[1:n-1,n+1:N]},'r');
end

end

function Phi = calculatePhi(X, M, R, n, Pi, epsilon)

if (isa(X,'sptensor'))
    Phi = -ones(size(X,n),R);
    xsubs = X.subs(:,n);
    v = sum(M.U{n}(xsubs,:).*Pi,2);
    wvals = X.vals ./ max(v, epsilon);
    for r = 1:R
        Yr = accumarray(xsubs, wvals .* Pi(:,r), [size(X,n) 1]);
        Phi(:,r) = Yr;
    end
else
    Xn = double(tenmat(X,n));
    V = M.U{n}*Pi';
    W = Xn ./ max(V, epsilon);
    Y = W * Pi;
    Phi = Y;
end

end

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

function y = vectorizeForMu(x)
y = x(:);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Shared Internal Functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Pi = tt_calcpi_prowsubprob(X, isSparse, M, R, n, N, sparse_indices)
% TT_CALCPI_PROWSUBPROB Compute Pi for a row subproblem.
%
%   X              - data tensor
%   isSparse       - true if X is sparse, false if dense
%   M              - current factor matrices
%   R              - number of columns in each factor matrix
%   n              - mode
%   N              - number of modes (equals the number of factor matrices)
%   sparse_indices - indices of row subproblem nonzero elements
%
%   Returns Pi matrix.
%
%   Intended for use by CP_PDN and CP_PQN.
%   Based on calculatePi() in CP_APR, which computes for an entire mode
%   instead of a single row subproblem.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

    if (isSparse)
        % X is a sparse tensor.  Compute Pi for the row subproblem specified
        % by sparse_indices.
        num_row_nnz = length(sparse_indices);

        Pi = ones(num_row_nnz, R);
        for nn = [1:n-1,n+1:N]
            Pi = M{nn}(X.subs(sparse_indices,nn),:) .* Pi;
        end
    else
        % X is a dense tensor.  Compute Pi for all rows in the mode.
        U = M.U;
        Pi = khatrirao(U{[1:n-1,n+1:N]},'r');
    end

end

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

function [m_new, f_old, f_1, f_new, num_evals] ...
    = tt_linesearch_prowsubprob(d, grad, m_old, step_len, step_red, ...
                                max_steps, suff_decr, isSparse, x_row, Pi, ...
                                phi_row, disp_warn)
% TT_LINESEARCH_PROWSUBPROB Perform a line search on a row subproblem.
%
%   d         - search direction
%   grad      - gradient vector at m_old
%   m_old     - current variable values
%   step_len  - initial step length, which is the maximum possible step length
%   step_red  - step reduction factor (suggest 1/2)
%   max_steps - maximum number of steps to try (suggest 10)
%   suff_decr - sufficient decrease for convergence (suggest 1.0e-4)
%   isSparse  - sparsity flag for computing the objective
%   x_row     - row subproblem data, for computing the objective
%   Pi        - Pi matrix, for computing the objective
%   phi_row   - 1-grad, more accurate if failing over to multiplicative update
%   disp_warn - true means warning messages are displayed
%
%   Returns
%     m_new     - new (improved) variable values
%     num_evals - number of times objective was evaluated
%     f_old     - objective value at m_old
%     f_1       - objective value at m_old + step_len * d
%     f_new     - objective value at m_new

    minDescentTol = 1.0e-7;
    smallStepTol = 1.0e-7;

    stepSize = step_len;

    % Evaluate the current objective value.
    f_old = -1 * tt_loglikelihood_row(isSparse, x_row, m_old, Pi);
    num_evals = 1;
    count = 1;

    while (count <= max_steps)
        % Compute a new step and project it onto the positive orthant.
        m_new = m_old + (stepSize .* d);
        m_new = m_new .* (m_new > 0);

        % Check that it is a descent direction.
        gDotd = sum(grad .* (m_new - m_old));
        if (gDotd > 0) || (sum(m_new) < minDescentTol)
            % Don't evaluate the objective if not a descent direction
            % or if all of the elements of m_new are close to zero.
            f_new = Inf;
            if (count == 1)
               f_1 = f_new;
            end

            stepSize = stepSize * step_red;
            count = count + 1;
        else
            % Evaluate objective function at new iterate.
            f_new = -1 * tt_loglikelihood_row(isSparse, x_row, m_new, Pi);
            num_evals = num_evals + 1;
            if (count == 1)
               f_1 = f_new;
            end

            % Check for sufficient decrease.
            if (f_new <= f_old + suff_decr * gDotd)
                break;
            else
                stepSize = stepSize * step_red;
                count = count + 1;
            end
        end
    end

    % Check if the line search failed.
    if (isinf(f_1) == 1)
        % Unit step failed; return a value that yields ared = 0.
        f_1 = f_old;
    end
    if (   ((count >= max_steps) && (f_new > f_old)) ...
        || (sum(m_new) < smallStepTol) )

        % Fall back on a multiplicative update step (scaled steepest descent).
        % Experiments indicate it works better than a unit step in the direction
        % of steepest descent, which would be the following:
        % m_new = m_old - (step_len * grad);     % steepest descent
        % A simple update formula follows, but suffers from round-off error
        % when phi_row is tiny:
        % m_new = m_old - (m_old .* grad);
        % Use this for best accuracy:
        m_new = m_old .* phi_row;                % multiplicative update

        % Project to the constraints and reevaluate the subproblem objective.
        m_new = m_new .* (m_new > 0);
        f_new = -1 * tt_loglikelihood_row(isSparse, x_row, m_new, Pi);
        num_evals = num_evals + 1;

        % Let the caller know the search direction made no progress.
        f_1 = f_old;

        if (disp_warn)
            fprintf('WARNING: line search failed, using multiplicative update step\n');
        end
    end

end

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

function f = tt_loglikelihood_row(isSparse, x, m, Pi)
%TT_LOGLIKELIHOOD_ROW Compute log-likelihood of one row subproblem.
%
%    The row subproblem for a given mode includes one row of matricized tensor
%    data (x) and one row of the model (m) in the same matricized mode.
%    Then
%       (dense case)
%          m:  R-length vector
%          x:  J-length vector
%          Pi: R x J matrix
%       (sparse case)
%          m:  R-length vector
%          x:  p-length vector, where p = nnz in row of matricized data tensor
%          Pi: R x p matrix
%       F = - (sum_r m_r - sum_j x_j * log (m * Pi_j)
%           where Pi_j denotes the j^th column of Pi
%           NOTE: Rows of Pi' must sum to one
%
%   isSparse - true if x is sparse, false if dense
%   x        - vector of data values
%   m        - vector of model values
%   Pi       - matrix
%
%   Returns the log-likelihood probability f.
%
%   Intended for use by CP_PDN and CP_PQN.
%   Similar to tt_loglikelihood() in CP_APR, which computes log likelihood
%   for the entire tensor instead of a single row subproblem.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

    term1 = -sum(m);

    if (isSparse)
        term2 = sum(x' .* log(m * Pi'));
    else
        b_pi = m * Pi';
        term2 = 0;
        for i = 1:length(x)
            if (x(i) == 0)
                % Define zero times log(anything) to be zero.
            else
                term2 = term2 + x(i) .* log(b_pi(i));
            end
        end
    end

    f = term1 + term2;

end

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

function f = tt_loglikelihood(X,M)
%TT_LOGLIKELIHOOD Compute log-likelihood of data X with model M.
%
%   F = TT_LOGLIKELIHOOD(X,M) computes the log-likelihood of model M given
%   data X, where M is a ktensor and X is a tensor or sptensor.
%   Specifically, F = - (sum_i m_i - x_i * log_i) where i is a multiindex
%   across all tensor dimensions.
%
%   See also cp_apr, tensor, sptensor, ktensor.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

N = ndims(X);

if ~isa(M, 'ktensor')
    error('M must be a ktensor');
end

M = normalize(M,1,1);

if isa(X, 'sptensor')
    xsubs = X.subs;
    A = M.U{1}(xsubs(:,1),:);
    for n = 2:N
       A = A .* M.U{n}(xsubs(:,n),:);
    end
    f = sum(X.vals .* log(sum(A,2))) - sum(sum(M.U{1}));
else
%{
% Old code is probably faster, but returns NaN if X and M are both zero
% for some element.
    f = sum(sum(double(tenmat(X,1)) .* log(double(tenmat(M,1))))) - sum(sum(M.U{1}));
%}
    % The check for x==0 is also in tt_loglikelihood_row.
    dX = double(tenmat(X,1));
    dM = double(tenmat(M,1));
    f = 0;
    for i = 1:size(dX,1)
      for j = 1:size(dX,2)
        if (dX(i,j) == 0.0)
          % Define zero times log(anything) to be zero.
        else
          f = f + dX(i,j) * log(dM(i,j));
        end
      end
    end
    f = f - sum(sum(M.U{1}));

end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% Test file for SPINOP2:

function pass = test_spinop2()

% Construction from STRING for GL equation:
S = spinop2('GL');
N = S.nonlin;
pass(1) = strcmpi(func2str(N), '@(u)u-(1+1.5i)*u.*(abs(u).^2)');

% Construction from DOM/TSPAN:
dom = [0 2*pi 0 2*pi];
tspan = [0 1];
S = spinop2(dom, tspan);
pass(2) = isequal(S.domain, dom);
pass(3) = isequal(S.tspan, tspan);

% Test recursive SUBSREF:
pass(4) = isequal(S.domain([2 4]), [2*pi 2*pi]);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun