text:
function [volInfo, dat] = iimg_read_img(inputimgs, extended_output_flag, reading_data, first_vol)
% :Usage:
% ::
%
%     [volInfo, dat] = iimg_read_img(inputimgs, [extended_output_flag (1|2)], [reading_data (0|1)], [first vol only (0|1)]);
%
% - reads 3-D or 4-D img/nii files, with extended volInfo output
% - flag to read first volume only
% - can read iimg_data vectorized form as well as file names (just passes dat to output)
% - uses SPM
%
% volInfo fields:
%   - all fields returned by spm_vol, plus:
%   - .nvox - number of voxels in 3d image
%  - .image_indx - logical index (e.g., [0 1 1 0 1 ...]) of all nonzero voxels
%  - .wh_inmask (extended_output_flag > 0) - linear index (e.g., [2 3 5 ...]) - same as image_indx, but as a result of find
%  - .n_inmask (extended_output_flag > 0) - length of .wh_inmask
%  - .xyzlist (extended_output_flag > 0) - voxel coords of .wh_inmask voxels
%  - .cluster (extended_output_flag > 1) - cluster structure of in-mask voxels
%
% Extended output flag values:
%   1: Add xyz coord list and linear index of which voxels nonzero & non-nan in mask
%   2: Add clusters from spm_clusters
%
% Enforces that dat is a vector of data values, regardless of input.
%
% :Examples:
% ::
%
%    imname = 'rob_tmap_0001_filt_t_3-05_k5_pos.img';
%    volInfo = iimg_read_img(imname);
%
% WARNING:
% The behavior of this function is SPM version dependent.
% volInfo only contains the spm_vol structure of the FIRST image!!!
% Data in SPM5/8 should be 4-D, but in SPM2 should be 3-D
%
% Example of image reading/writing
% -------------------------------------------------------------------------
%   1) Get volume info from first volume of a 3-D or 4-D image
%       volInfo structure has necessary info for converting to/from
%       "vectorized" format
%       dat returns 4-D data from entire image (all volumes)
%       ::
%
%           iimg_name = 'test_run1_pca.img';
% [maskInfo, dat] = iimg_read_img(img_name, 2);
%
%   2) Get data in "vectorized" image format for each volume in the
%      image. Works for a list of images too. data is in-mask voxels x volumes
%      this can be useful when you want to return whole-brain data for many
%      images, but in a search volume only (i.e., no extra-brain voxels)
%      ::
%
%           idata = iimg_get_data(maskInfo, img_name);
%           data = data';  % make sure columns are volumes
%
% % 3) Write out a 4-D image with the same data, called test_run1_pca2.img
%       ::
%
%           ivoldat3D = iimg_reconstruct_vols(data, maskInfo, 'outname',
%                                                 'test_run1_pca2.img');

% ..
%    Set up arguments
% ..
if ~exist('inputimgs', 'var') || isempty(inputimgs)

    error('Input image not found. Check file and path names.');

end

if ~exist('extended_output_flag', 'var') || isempty(extended_output_flag)
    extended_output_flag = 0;
end
if ~exist('reading_data', 'var') || isempty(reading_data)
    reading_data = 1;
end

if ~exist('first_vol', 'var') || isempty(first_vol)
    first_vol = 0;
end

if(~reading_data && extended_output_flag)
    warning('If not reading in file data, the extended_output_flag has no meaning');
end

imgtype = get_img_type(inputimgs);

volInfo = [];
num_imgs = size(inputimgs, 1);
if first_vol, num_imgs = 1; end

switch imgtype
    case 'name'
        % ---------------------------------------------------------------
        % * String matrix of filenames: Load first volume, then all data
        % ---------------------------------------------------------------
        % use iimg_read_from_file

        space_defining_image = inputimgs(1, :);

%         if first_vol
%             whcomma = find(space_defining_image == ',');
%             if ~isempty(whcomma)
%                 space_defining_image = space_defining_image(1:whcomma-1);
%             end
%
%             space_defining_image = [space_defining_image ',1'];
%
%         end

        % Read first image, or first set of volumes in 4-d file:
        % If firstvol, the ",1" will be added to space defining image after
        % (possible) unzipping and checks on whether it exists.

        [volInfo, dat] = iimg_read_from_file(space_defining_image, extended_output_flag, reading_data, first_vol);

        if(reading_data && num_imgs > 1)
            % load other images; check to make sure they're in same dims as first

            dat = [dat zeros(volInfo.nvox, num_imgs-1)];

            for i = 2:num_imgs
                [indx data] = iimg_read_from_file(inputimgs(i, :), 0, reading_data, 0); % read all, not just first vol

                if length(data) ~= volInfo.nvox
                    fprintf('Warning: img dims for %s do not match 1st image. Resampling.\n', inputimgs(i, :));

                    data = scn_map_image(inputimgs(i, :), inputimgs(1, :));
                    data = data(:);

                end

                dat(:, i) = data;
            end
        end

    case '3dvols'
        % --------------------------------------
        % * matrix of 3-D image data volumes
        % --------------------------------------

        [x, y, z, n] = size(inputimgs);
        nvox = prod(x, y, z);
        dat = zeros(nvox, n);

        for i = 1:num_imgs
            tmp = inputimgs(:, :, :, i);
            dat(:, i) = tmp(:);
        end

        dat(isnan(dat)) = 0;

    case 'indx'
        % --------------------------------------
        % * Already in index -- target format
        % --------------------------------------
        dat = inputimgs;
end
end

function imgtype = get_img_type(inputimgs)
if ischar(inputimgs)
    imgtype = 'name';

    return
end

[x, y, z] = size(inputimgs);

if isa(inputimgs, 'statistic_image')
    imgtype = 'statistic_image';
elseif z > 1
    % 3-D image(s)
    imgtype = '3dvols';
elseif x > 1
    % voxels x images 2-D index matrix
    imgtype = 'indx';
else
    error('I don''t recognize the data format of inputimgs.');
end
end

% iimg_read_from_file
% -------------------------------------------------------------------------

function [volInfo, dat] = iimg_read_from_file(imname, extended_output_flag, reading_data, dofirstvolonly)
% Read image data from file using spm_vol and spm_read_vols.
% Call vol_file_check, which checks zipped and unzipped versions and
% gunzips if needed

% Returns filename for unzipped image with formatting ,# to read first vol
% only if requested
imname = vol_file_check(imname, dofirstvolonly);

% %  tor: oct 2010: speed up by adding ,1 to get only first vol in spm
% [pth, nam, ext] = fileparts(imname);
% wh_comma = find(ext == ',');
% if isempty(wh_comma)
%     imname = [imname ',1'];
% end

volInfo = spm_vol(imname);

% Changed Feb 2008 to deal with 4-D image reading in SPM5 and above
nimgs = length(volInfo);

for i = 1:nimgs
    volInfo(i).nvox = prod(volInfo(i).dim(1:3));
end

dat = [];

if(reading_data)

    switch spm('Ver')
        case 'SPM2'
            % spm_defaults is a script
            disp('WARNING: spm defaults not set for spm2. Make sure your defaults are set correctly');

        case {'SPM5', 'SPM8', 'SPM12'}
            % spm_defaults is a function
            spm_defaults()

        otherwise
            % unknown SPM
            disp('Unknown version of SPM!');
            spm_defaults()
    end

    % pre-allocate array
    dat = zeros(volInfo(1).nvox, nimgs);

    for i = 1:nimgs

        dat_vol = spm_read_vols(volInfo(i));
        dat(:, i) = dat_vol(:);

        volInfo(i).image_indx = dat(:) ~= 0; % locate all voxels in-mask

    end

    if extended_output_flag
        % Return this for first image only right now.
        volInfo(1).wh_inmask = find(volInfo(1).image_indx);    % indices of all voxels in mask

        volInfo(1).n_inmask = size(volInfo(1).wh_inmask, 1);

        [i, j, k] = ind2sub(volInfo(1).dim(1:3), volInfo(1).wh_inmask);
        volInfo(1).xyzlist = [i j k];

        if extended_output_flag > 1
            if volInfo(1).n_inmask < 50000
                volInfo(1).cluster = spm_clusters(volInfo(1).xyzlist')';
            else
                volInfo(1).cluster = ones(volInfo(1).n_inmask, 1);
            end
        end
    end

    dat(isnan(dat)) = 0;
end

% Return only the first volInfo struture because it contains generic info
% about all volumes.  Don't need to return every volume for our
% purposes.
volInfo = volInfo(1);

end

function [imname, waszipped] = vol_file_check(imname, varargin)
% Check that filename exists, in either gz (zipped) or unzipped format, and
% return file name of file with full path name.
% If both zipped and unzipped exist, prefer unzipped.
% Return flag waszipped so file can be re-zipped after loading.

waszipped = false;

if length(varargin) > 0
    dofirstvolonly = varargin{1};
end

if(~ischar(imname) || isempty(imname))
    error('Image name is not string or is empty.');
end

imname = deblank(imname);

% strip ,# volume numbers sometimes added by SPM
imname = regexprep(imname, ',\w*', '');

% strip .gz version for now and create .gz and non-.gz version. check for both.
imname = regexprep(imname, '.gz', '');

gzipexists = exist([imname '.gz'], 'file');
nozipexists = exist(imname, 'file');

if nozipexists
    % use imname

elseif gzipexists
    % use .gz - unzip and use previous value of imname
    % flag for re-gzipping later.
    waszipped = true;
    system(['gunzip ' imname '.gz']);

else
    error('Image "%s" does not exist or cannot be found on path (unzipped or .gz versions)!', imname);
end

full_imname = which(imname); % sometimes ok, sometimes empty if full path already included

if exist(full_imname, 'file') % if not empty and ok, else default to previous value of imname
    imname = full_imname;
end

% Add volume number if first vol only
if dofirstvolonly

    imname = [imname ',1'];

end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [model,b] = robustlssvm(model,ab,X,Y)
% Robust training in the case of non-Gaussian noise or outliers
%(only possible with the object oriented interface)
%
% >> model     = robustlssvm(model)
%
% Robustness towards outliers can be achieved by reducing the
% influence of support values corresponding to large errors.
%
%
% Full syntax
%
%     1. Using the object oriented interface:
%
% >> model = robustlssvm(model)
%
%       Outputs
%         model : Robustly trained object oriented representation of the LS-SVM model
%       Inputs
%         model : Object oriented representation of the LS-SVM model
%
% See also:
%   trainlssvm, tunelssvm, crossvalidate

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

if iscell(model),
    func = 1;
    model = initlssvm(model{:});
else
    func = 0;
end

if model.type(1)~='f',
    error('Robustly weighted least squares only implemented for regression case...');
end

if nargin>1,
    if iscell(ab) && ~isempty(ab),
        model.alpha = ab{1};
        model.b = ab{2};
        model.status = 'trained';
        if nargin>=4,
            model = trainlssvm(model,X,Y);
        end
    else
        model = trainlssvm(model,ab,X);
    end
else
    model = trainlssvm(model);
end

% model errors
ek = model.alpha./model.gam';
g = model.gam;
%
% robust estimation of the variance
%
eval('delta=model.delta;','delta=[];')
for j=1:500
    vare = 1.483*median(abs((ek)-median(ek)));
    alphaold = model.alpha;
    %
    % robust re-estimation of the alpha's and the b
    %
    cases = reshape((ek./vare),1,model.nb_data);
    W = weightingscheme(cases,model.weights,delta);
    W = g*W;

    model = changelssvm(model,'gam',W);
    %     model = changelssvm(model,'implementation','MATLAB');
    model = trainlssvm(model);
    ek = model.alpha./model.gam';

    if norm(abs(alphaold-model.alpha),'fro')<=1e-4,
        fprintf('Converged after %.0f iteration(s)', j);
        if func && nargout~=1,
            b = model.b;
            model = model.alpha;
        end
        return
    end
    model.status = 'changed';
end

author: peterhcharlton

repo: RRest


text:
function trbvp
%TRBVP  Exercise for Example 3 of the BVP tutorial.
%   This problem is studied in section 5.4 of B.A. Finlayson, The Method of
%   Weighted Residuals and Variational Principles, Academic, New York, 1972.
%   It arises when modelling a tubular reactor with axial dispersion.  An
%   isothermal situation with n-th order irreversible reaction leads to
%
%      y'' = Pe*(y' - R*y^n)
%
%   Here Pe is the axial Peclet number and R is the reaction rate group.
%   The boundary conditions are
%
%      y'(0) = Pe*(y(0) - 1),    y'(1) = 0.
%
%   Finlayson compares results he obtains with an orthogonal collocation
%   method to results obtained by others with finite differences when
%   Pe = 1, R = 2, and n = 2.  These results are compared here to results
%   obtained with BVP4C. BVPVAL is used to get a smoother graph of y(x).

% Copyright 1999, The MathWorks, Inc.

% Known parameter
Pe = 1;

options = bvpset('stats','on');
solinit = bvpinit(linspace(0,1,5),[0.5 0]);
sol = bvp4c(@trode,@trbc,solinit,options,Pe);

fprintf('\n');
fprintf('Other authors report y(0) = 0.63678, y(1) = 0.45759.\n');
fprintf('Values computed are  y(0) = %7.5f, y(1) = %7.5f\n',sol.y(1,1),sol.y(1,end));

clf reset
xint = linspace(0,1);
Sxint = bvpval(sol,xint);
plot(xint,Sxint(1,:))
title('Mass transfer in a tubular reactor.')
xlabel('x')
ylabel('y')
shg

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

function dydx = trode(x,y,Pe)
%TRODE  ODE function for the exercise of Example 3 of the BVP tutorial.
dydx = [ y(2)
         Pe*(y(2) + 2*y(1)^2)];

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

function res = trbc(ya,yb,Pe)
%TRBC  Boundary conditions for the exercise of Example 3 of the BVP tutorial.
res = [ ya(2) - Pe*(ya(1) - 1)
        yb(2) ];

author: Sable

repo: mcbench-benchmarks


text:
function varargout = sift(varargin)
% VL_SIFT  Scale-Invariant Feature Transform
%   F = VL_SIFT(I) computes the SIFT frames [1] (keypoints) F of the
%   image I. I is a gray-scale image in single precision. Each column
%   of F is a feature frame and has the format [X;Y;S;TH], where X,Y
%   is the (fractional) center of the frame, S is the scale and TH is
%   the orientation (in radians).
%
%   [F,D] = VL_SIFT(I) computes the SIFT descriptors [1] as well. Each
%   column of D is the descriptor of the corresponding frame in F. A
%   descriptor is a 128-dimensional vector of class UINT8.
%
%   VL_SIFT() accepts the following options:
%
%   Octaves:: [maximum possible]
%     Set the number of octave of the DoG scale space.
%
%   Levels:: [3]
%     Set the number of levels per octave of the DoG scale space.
%
%   FirstOctave:: [0]
%     Set the index of the first octave of the DoG scale space.
%
%   PeakThresh:: [0]
%     Set the peak selection threshold.
%
%   EdgeThresh:: [10]
%     Set the non-edge selection threshold.
%
%   NormThresh:: [-inf]
%     Set the minimum l2-norm of the descriptors before
%     normalization. Descriptors below the threshold are set to zero.
%
%   Magnif:: [3]
%     Set the descriptor magnification factor. The scale of the
%     keypoint is multiplied by this factor to obtain the width (in
%     pixels) of the spatial bins. For instance, if there are there
%     are 4 spatial bins along each spatial direction, the
%     ``side'' of the descriptor is approximatively 4 * MAGNIF.
%
%   WindowSize:: [2]
%     Set the variance of the Gaussian window that determines the
%     descriptor support. It is expressend in units of spatial
%     bins.
%
%   Frames:: [not specified]
%     If specified, set the frames to use (bypass the detector). If
%     frames are not passed in order of increasing scale, they are
%     re-orderded.
%
%   Orientations::
%     If specified, compute the orietantions of the frames overriding
%     the orientation specified by the 'Frames' option.
%
%   Verbose::
%     If specfified, be verbose (may be repeated to increase the
%     verbosity level).
%
%   REFERENCES::
%     [1] D. G. Lowe, Distinctive image features from scale-invariant
%     keypoints. IJCV, vol. 2, no. 60, pp. 91-110, 2004.
%
%   See also: VL_UBCMATCH(), VL_DSIFT(), VL_HELP().
[varargout{1:nargout}] = vl_sift(varargin{:});

author: jianxiongxiao

repo: ProfXkit


text:
function save_avw(img,fname,vtype,vsize)
% SAVE_AVW(img,fname,vtype,vsize)
%
%  Create and save an analyse header (.hdr) and image (.img) file
%   for either a 2D or 3D or 4D array (automatically determined).
%  fname is the filename (must be inside single quotes)
%
%  vtype is 1 character: 'b'=unsigned byte, 's'=short, 'i'=int, 'f'=float
%                        'd'=double or 'c'=complex
%  vsize is a vector [x y z tr] containing the voxel sizes in mm and
%  the tr in seconds  (defaults: [1 1 1 3])
%
%  See also: READ_AVW
%

   if ((~isreal(img)) & (vtype~='c')),
     disp('WARNING:: Overwriting type - saving as complex');
     save_avw_complex(img,fname,vsize);
   else
     if (vtype=='c'),
       save_avw_complex(img,fname,vsize);
     else
       save_avw_hdr(img,fname,vtype,vsize);
       save_avw_img(img,fname,vtype);
     end
   end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function sigm = nonlinearity_nil(x)
    sigm = x;
end

author: jimmy-ren


text:
function dUpsilon = lfmvpGradientUpsilonMatrix(gamma, sigma2, t1, ...
    t2, mode, upsilon)

% LFMVPGRADIENTUPSILONMATRIX Gradient upsilon matrix vel. pos.
% FORMAT
% DESC computes the gradient of a portion of the LFM kernel.
% ARG gamma : Gamma value for system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG upsilon : precomputation of the upsilon matrix.
% ARG mode : operation mode, according to the derivative (mode 0,
% derivative wrt t1, mode 1 derivative wrt t2)
% RETURN upsilon : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : Mauricio A. Alvarez, 2010
%
% SEEALSO : lfmvpComputeUpsilonMatrix.m

% KERN

sigma = sqrt(sigma2);

if nargin<6
    upsilon = lfmComputeUpsilonMatrix(gamma, sigma2, t1, t2);
    if nargin <5
        mode =0;
    end
end

if mode ==0
    dUpsilon = -upsilon - gamma*lfmGradientUpsilonMatrix(gamma, sigma2, t1, t2);
else
    dUpsilon = upsilon + gamma*lfmGradientUpsilonMatrix(gamma, sigma2, t1, t2) ...
        - (2/(sqrt(pi)*sigma))*(t1.*exp(-gamma*t1))*(exp(-(t2.^2)/sigma2)).';
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function z = distSqr_fast(x,y,x2,y2) %x2 = sum(x.^2,1)'; %y2 = sum(y.^2,1);
% function z = distSqr_fast(x,y,x2,y2)
%
% Return matrix of all-pairs squared distances between the vectors
% in the columns of x and y.
%
% INPUTS
% 	x 	dxn matrix of vectors
% 	y 	dxm matrix of vectors
%
% OUTPUTS
% 	z 	nxm matrix of squared distances
%
% This routine is faster when m<n than when m>n. In other words y should be
% smaller than x.
%
% David Martin <dmartin@eecs.berkeley.edu>
% March 2003

% Based on dist2.m code,
% Copyright (c) Christopher M Bishop, Ian T Nabney (1996, 1997)

if ~exist('y','var')
  y = x;
end

if ~exist('x2','var')
  x2 = sum(x.^2,1)';
end

if ~exist('y2','var')
  y2 = sum(y.^2,1);
end

if size(x,1) ~= size(y,1),
  error('size(x,1)~=size(y,1)');
end

[d,n] = size(x);
[d,m] = size(y);

%z = x'*y
%z = repmat(x2,1,m) ...
%    + repmat(y2,n,1) ...
%    - 2*x'*y;
%return

z = x'*y;

for i = 1:m,
  z(:,i) = x2 + y2(i) - 2*z(:,i);
end

%z = zeros(n,m);
%for i = 1:m,
%  z(:,i) = x2 + y2(i) - 2*x'*y(:,i);
%end

author: quantombone

repo: exemplarsvm


text:
%GETSIZE Return size of a mapping
%
%   SIZE = GETSIZE(W)
%
% This function is identical to SIZE(W).

author: marianux

repo: ecg-kit

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


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

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

% Last Modified by GUIDE v2.5 16-Jul-2010 19:34:10

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

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

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

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

% Update handles structure
guidata(hObject, handles);

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

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

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

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
programma
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
dtmf_story
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes during object creation, after setting all properties.
function axes1_CreateFcn(hObject, eventdata, handles)
imshow(imread('cellulare.jpg'))
% hObject    handle to axes1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate axes1

% --- Executes during object creation, after setting all properties.
function axes2_CreateFcn(hObject, eventdata, handles)
imshow(imread('logo-sun.jpg'))
% hObject    handle to axes2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate axes2


text:
function FB = FbMake( dim, flag, show )
% Various 1D/2D/3D filterbanks (hardcoded).
%
% USAGE
%  FB = FbMake( dim, flag, [show] )
%
% INPUTS
%  dim     - dimension
%  flag    - controls type of filterbank to create
%          - if d==1
%            1: gabor filter bank for spatiotemporal stuff
%          - if d==2
%            1: filter bank from Serge Belongie
%            2: 1st/2nd order DooG filters.  Similar to Gabor filterbank.
%            3: similar to Laptev&Lindberg ICPR04
%            4: decent seperable steerable? filterbank
%            5: berkeley filterbank for textons papers
%            6: symmetric DOOG filters
%          - if d==3
%            1: decent seperable steerable filterbank
%  show    - [0] figure to use for optional display
%
% OUTPUTS
%
% EXAMPLE
%  FB = FbMake( 2, 1, 1 );
%
% See also FBAPPLY2D
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<3 || isempty(show) ); show=0; end

% create FB
switch dim
  case 1
    FB = FbMake1D( flag );
  case 2
    FB = FbMake2D( flag );
  case 3
    FB = FbMake3d( flag );
  otherwise
    error( 'dim must be 1 2 or 3');
end

% display
FbVisualize( FB, show );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function FB = FbMake1D( flag )
switch flag
  case 1  %%% gabor filter bank for spatiotemporal stuff
    omegas = 1 ./ [3 4 5 7.5 11];
    sigmas =      [3 4 5 7.5 11];
    FB = FbMakegabor1D( 15, sigmas, omegas );

  otherwise
    error('none created.');
end

function FB = FbMakegabor1D( r, sigmas, omegas )
for i=1:length(omegas)
  [feven,fodd]=filterGabor1d(r,sigmas(i),omegas(i));
  if( i==1 ); FB=repmat(feven,[2*length(omegas) 1]); end
  FB(i*2-1,:)=feven; FB(i*2,:)=fodd;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function FB = FbMake2D( flag )

switch flag
  case 1  %%% filter bank from Berkeley / Serge Belongie
    r=15;
    FB = FbMakegabor( r, 6, 3, 3, sqrt(2) );
    FB2 = FbMakeDOG( r, .6, 2.8, 4);
    FB = cat(3, FB, FB2);
    %FB = FB(:,:,1:2:36); %include only even symmetric filters
    %FB = FB(:,:,2:2:36); %include only odd symmetric filters

  case 2 %%% 1st/2nd order DooG filters.  Similar to Gabor filterbank.
    FB = FbMakeDooG( 15, 6, 3, 5, .5) ;

  case 3 %%% similar to Laptev&Lindberg ICPR04
    % Wierd filterbank of Gaussian derivatives at various scales
    % Higher order filters probably not useful.
    r = 9;   dims=[2*r+1 2*r+1];
    sigs = [.5 1 1.5 3]; % sigs = [1,1.5,2];

    derivs = [];
    %derivs = [ derivs; 0 0 ]; % 0th order
    %derivs = [ derivs; 1 0; 0 1 ]; % first order
    %derivs = [ derivs; 2 0; 0 2; 1 1 ]; % 2nd order
    %derivs = [ derivs; 3 0; 0 3; 1 2; 2 1 ]; % 3rd order
    %derivs = [ derivs; 4 0; 0 4; 1 3; 3 1; 2 2 ]; % 4th order
    derivs = [ derivs; 0 1; 0 2; 0 3; 0 4; 0 5]; % 0n order
    derivs = [ derivs; 1 0; 2 0; 3 0; 4 0; 5 0]; % n0 order
    cnt=1;  nderivs = size(derivs,1);
    for s=1:length(sigs)
      for i=1:nderivs
        dG = filterDoog( dims, [sigs(s) sigs(s)], derivs(i,:), 0 );
        if(s==1 && i==1); FB=repmat(dG,[1 1 length(sigs)*nderivs]); end
        FB(:,:,cnt) = dG; cnt=cnt+1;
        %dG = filterDoog( dims, [sigs(s)*3 sigs(s)], derivs(i,:), 0 );
        %FB(:,:,cnt) = dG; cnt=cnt+1;
        %dG = filterDoog( dims, [sigs(s) sigs(s)*3], derivs(i,:), 0 );
        %FB(:,:,cnt) = dG; cnt=cnt+1;
      end
    end

  case 4 % decent seperable steerable? filterbank
    r = 9;   dims=[2*r+1 2*r+1];
    sigs = [.5 1.5 3];
    derivs = [1 0; 0 1; 2 0; 0 2];
    cnt=1;  nderivs = size(derivs,1);
    for s=1:length(sigs)
      for i=1:nderivs
        dG = filterDoog( dims, [sigs(s) sigs(s)], derivs(i,:), 0 );
        if(s==1 && i==1); FB=repmat(dG,[1 1 length(sigs)*nderivs]); end
        FB(:,:,cnt) = dG; cnt=cnt+1;
      end
    end
    FB2 = FbMakeDOG( r, .6, 2.8, 4);
    FB = cat(3, FB, FB2);

  case 5  %%% berkeley filterbank for textons papers
    FB = FbMakegabor( 7, 6, 1, 2, 2 );

  case 6  %%% symmetric DOOG filters
    FB = FbMakeDooGSym( 4, 2, [.5 1] );

  otherwise
    error('none created.');
end

function FB = FbMakegabor( r, nOrient, nScales, lambda, sigma )
% multi-scale even/odd gabor filters. Adapted from code by Serge Belongie.
cnt=1;
for m=1:nScales
  for n=1:nOrient
    [F1,F2]=filterGabor2d(r,sigma^m,lambda,180*(n-1)/nOrient);
    if(m==1 && n==1); FB=repmat(F1,[1 1 nScales*nOrient*2]); end
    FB(:,:,cnt)=F1;  cnt=cnt+1;   FB(:,:,cnt)=F2;  cnt=cnt+1;
  end
end

function FB = FbMakeDooGSym( r, nOrient, sigs )
% Adds symmetric DooG filters.  These are similar to gabor filters.
cnt=1; dims=[2*r+1 2*r+1];
for s=1:length(sigs)
  Fodd = -filterDoog( dims, [sigs(s) sigs(s)], [1 0], 0 );
  Feven = filterDoog( dims, [sigs(s) sigs(s)], [2 0], 0 );
  if(s==1); FB=repmat(Fodd,[1 1 length(sigs)*nOrient*2]); end
  for n=1:nOrient
    theta = 180*(n-1)/nOrient;
    FB(:,:,cnt) = imrotate( Feven, theta, 'bil', 'crop' );  cnt=cnt+1;
    FB(:,:,cnt) = imrotate( Fodd,  theta, 'bil', 'crop' );  cnt=cnt+1;
  end
end

function FB = FbMakeDooG( r, nOrient, nScales, lambda, sigma )
% 1st/2nd order DooG filters.  Similar to Gabor filterbank.
% Defaults: nOrient=6, nScales=3, lambda=5, sigma=.5,
cnt=1; dims=[2*r+1 2*r+1];
for m=1:nScales
  sigma = sigma * m^.7;
  Fodd = -filterDoog( dims, [sigma lambda*sigma^.6], [1,0], 0 );
  Feven = filterDoog( dims, [sigma lambda*sigma^.6], [2,0], 0 );
  if(m==1); FB=repmat(Fodd,[1 1 nScales*nOrient*2]); end
  for n=1:nOrient
    theta = 180*(n-1)/nOrient;
    FB(:,:,cnt) = imrotate( Feven, theta, 'bil', 'crop' );  cnt=cnt+1;
    FB(:,:,cnt) = imrotate( Fodd,  theta, 'bil', 'crop' );  cnt=cnt+1;
  end
end

function FB = FbMakeDOG( r, sigmaStr, sigmaEnd, n )
% adds a serires of difference of Gaussian filters.
sigs = sigmaStr:(sigmaEnd-sigmaStr)/(n-1):sigmaEnd;
for s=1:length(sigs)
  FB(:,:,s) = filterDog2d(r,sigs(s),2); %#ok<AGROW>
  if( s==1 ); FB=repmat(FB,[1 1 length(sigs)]); end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function FB = FbMake3d( flag )

switch flag
  case 1 % decent seperable steerable filterbank
    r = 25; dims=[2*r+1 2*r+1 2*r+1];
    sigs = [.5 1.5 3];
    derivs = [0 0 1; 0 1 0; 1 0 0; 0 0 2; 0 2 0; 2 0 0];
    cnt=1; nderivs = size(derivs,1);
    for s=1:length(sigs)
      for i=1:nderivs
        dG = filterDoog( dims, repmat(sigs(s),[1 3]), derivs(i,:), 0 );
        if(s==1 && i==1); FB=repmat(dG,[1 1 1 nderivs*length(sigs)]); end
        FB(:,:,:,cnt) = dG; cnt=cnt+1;
      end
    end

  otherwise
    error('none created.');
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function scrollplot(s)

    % Created by Steven Lord, slord@mathworks.com
    % Uploaded to MATLAB Central
    % http://www.mathworks.com/matlabcentral
    % 7 May 2002
    %
    % Permission is granted to adapt this code for your own use.
    % However, if it is reposted this message must be intact.

    % Generate and plot data
    x=s.dnum();
    y=s.data();
    dx=1;
    %% dx is the width of the axis 'window'
    a=gca;
    p=plot(x,y);

    % Set appropriate axis limits and settings
    set(gcf,'doublebuffer','on');
    %% This avoids flickering when updating the axis
    set(a,'xlim',[min(x) min(x)+dx]);
    set(a,'ylim',[min(y) max(y)]);

    % Generate constants for use in uicontrol initialization
    pos=get(a,'position');
    Newpos=[pos(1) pos(2)-0.1 pos(3) 0.05];
    %% This will create a slider which is just underneath the axis
    %% but still leaves room for the axis labels above the slider
    xmax=max(x);
    xmin=min(x);
    xmin=0;
    %gs = get(gcbo,'value')+[min(x) min(x)+dx]
    S=sprintf('set(gca,''xlim'',get(gcbo,''value'')+[%f %f])',[xmin xmin+dx])
    %% Setting up callback string to modify XLim of axis (gca)
    %% based on the position of the slider (gcbo)
    % Creating Uicontrol
    h=uicontrol('style','slider',...
        'units','normalized','position',Newpos,...
        'callback',S,'min',xmin,'max',xmax-dx);
        %'callback',S,'min',0,'max',xmax-dx);
end

author: geoscience-community-codes

repo: GISMO


text:
function asa063_test01 ( )

%*****************************************************************************80
%
%% ASA063_TEST01 demonstrates the use of BETAIN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 January 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA063_TEST01:\n' );
  fprintf ( 1, '  BETAIN computes the incomplete beta function.\n' );
  fprintf ( 1, '  Compare to tabulated values.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      A       B       X       ' );
  fprintf ( 1, 'FX                       FX2\n' );
  fprintf ( 1, '                              ' );
  fprintf ( 1, '(Tabulated)              (BETAIN)                DIFF\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, a, b, x, fx ] = beta_inc_values ( n_data );

    if ( n_data == 0 )
      break
    end

    beta_log = gammaln ( a ) ...
             + gammaln ( b ) ...
             - gammaln ( a + b );

    [ fx2, ifault ] = betain ( x, a, b, beta_log );

    fprintf ( 1, '  %6.2f  %6.2f  %6.2f  %24.16e  %24.16e  %10.4e\n', ...
    a, b, x, fx, fx2, abs ( fx - fx2 ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [sol,output] = sdf_nls(model,options)
%SDF_NLS Structured data fusion by nonlinear least squares.
%   [sol,output] = sdf_nls(model) solves the data fusion problem described
%   by the structure model and returns the solution as the structure sol.
%   The model describes a data fusion problem with three fields:
%
%      model.variables
%
%         A structure or cell array of initializations for the variables in
%         the data fusion problem. Each field or cell in model.variables
%         represents a variable. A variable may be an array (such as
%         scalars, vectors, matrices and tensors) or a (nested) cell array
%         of arrays.
%
%      model.factors
%
%         A structure or cell array of factors. Each field or cell in
%         model.factors represents a factor. A factor is described by a
%         cell array of subfactors. After each subfactor has been
%         generated, the factor is constructed as the cell2mat of its
%         subfactors. A subfactor is a cell array in which the first
%         element is a reference to a variable, and the following elements
%         represent a sequence of transformations of that variable. If
%         model.variables is a structure, then a reference to a variable is
%         a string corresponding to the field name of the desired variable.
%         If model.variables is a cell array, then a reference to a
%         variable is the index of the cell containing the desired
%         variable. Transformations are supplied by functions which
%         describe their (linearized) behaviour. For example, the
%         transformation @struct_inv computes the factor as the matrix
%         inverse of the selected variable. All transformations must have
%         the function signature struct_mytrans(z,task). Instead of
%         supplying a reference to a variable, it is also possible to
%         supply a constant factor in the form of an array.
%
%      model.factorizations
%
%         A structure of data sets to jointly factorize. Each field in
%         model.factorizations represents a factorization. A factorization
%         is described by a structure containing two fields. The first
%         field has the field name 'data' and contains the (dense, sparse
%         or incomplete) array which is to be factorized. If the array
%         contains many zeros or a NaN, it is internally converted to a
%         sparse or incomplete tensor with fmt. The second field's
%         field name designates the type of factorization to compute of the
%         data set. Currently, two types are supported:
%
%         1. Canonical polyadic decomposition: a field 'cpd' should contain
%            a cell array of references to factors in model.factors. The
%            nth reference in the cell array corresponds to the nth factor
%            matrix of the CPD.
%
%         2. Block term decomposition: a field 'btd' should contain a cell
%            array of terms. Each term is itself a cell array containing
%            references to factors in model.factors. The nth reference in a
%            term corresponds to the nth factor matrix of that term. The
%            (N+1)th reference, where N is the number of dimensions of the
%            data set, in a term corresponds to the core tensor of that
%            term.
%
%         Additionally, two types of regularization are available:
%
%         3. L2 regularization: a field 'regL2' should contain a cell array
%            of references to factors in model.factors. This appends a term
%            to the objective function of the form 0.5*norm(F(:)-D(:),2)^2,
%            where F(:) and D(:) are the serialized factors in regL2 and
%            the data field respectively. The data field may be omitted, in
%            which case it is an all-zero vector.
%
%         4. L1 regularization: a field 'regL1' should contain a cell array
%            of references to factors in model.factors. This appends a term
%            to the objective function which is a smooth approximation of
%            0.5*norm(F(:)-D(:),1), where F(:) and D(:) are the serialized
%            factors in regL1 and the data field respectively. The data
%            field may be omitted, in which case it is an all-zero vector.
%
%   The algorithm proceeds to compute the joint factorization of the data
%   sets provided by models.factorizations by minimizing the sum of square
%   magnitude residuals of these factorizations. The output sol contains
%   the optimized variables in sol.variables and the corresponding factors
%   in sol.factors.
%
%   The structure output returns additional information:
%
%      output.Name  - The name of the selected algorithm.
%      output.<...> - The output of the selected algorithm.
%
%   sdf_nls(model,options) may be used to set the following options:
%
%      options.Algorithm =   - The desired optimization method.
%      [@nls_gncgs| ...
%       {@nls_gndl}|@nls_lm]
%      options.PC = true     - Whether or not to use a preconditioner, if
%                              available, for computing the Gauss-Newton
%                              step.
%      options.RelWeights =  - By supplying relative weights, the weights
%      ones(1,F)               options.Weights are computed as follows:
%                              options.Weights(f) = options.RelWeights(f)/
%                              (sum(options.RelWeights)*numel(data_f)) for
%                              each of the F factorizations in the model.
%      options.Weights       - The weight of the F factorizations in the
%                              data fusion model. The SDF objective
%                              function is \sum_f 0.5*options.Weights(f)*
%                              frob(model_f-data_f)^2, where model_f and
%                              data_f are the fth factorization and
%                              corresponding tensor. By default, weights
%                              are provided by options.RelWeights, but
%                              options.Weights has precedence if supplied.
%      options.<...>         - Parameters passed to the selected method,
%                              e.g., options.TolFun, options.TolX.
%                              See also help [options.Algorithm].
%
%   See also sdf_minf.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Structured data fusion,"
%       ESAT-SISTA Internal Report 13-177, KU Leuven, 2013.

% Check the options structure.
if nargin < 2, options = struct; end
isfunc = @(f)isa(f,'function_handle');
xsfunc = @(f)isfunc(f)&&exist(func2str(f),'file');
if ~isfield(options,'Algorithm')
    funcs = {@nls_gndl,@nls_gncgs,@nls_lm};
    options.Algorithm = funcs{find(cellfun(xsfunc,funcs),1)};
end
if ~isfield(options,'MaxIter'), options.MaxIter = 5000; end
if ~isfield(options,'CGMaxIter'), options.CGMaxIter = 15; end
if ~isfield(options,'Display'), options.Display = 0; end
if ~isfield(options,'JHasFullRank'), options.JHasFullRank = false; end
if ~isfield(options,'PC'), options.PC = true; end
if ~isfield(options,'TolLargeScale'), options.TolLargeScale = 0.02; end

% Convert model to internal format.
if ~iscell(model.variables)
    cache.variables.names = fieldnames(model.variables);
    model.variables = reshape(struct2cell(model.variables),1,[]);
else
    model.variables = model.variables(:).';
end
if ~iscell(model.factors)
    cache.factors.names = fieldnames(model.factors);
    model.factors = reshape(struct2cell(model.factors),1,[]);
else
    model.factors = model.factors(:).';
end
if ~iscell(model.factorizations)
    cache.factorizations.names = fieldnames(model.factorizations);
    model.factorizations = reshape(struct2cell(model.factorizations),1,[]);
else
    model.factorizations = model.factorizations(:).';
end
for I = 1:length(model.factors)
    if ~iscell(model.factors{I})
        model.factors{I} = model.factors(I);
    end
    if any(cellfun(@(f)isa(f,'function_handle'),model.factors{I}))
        model.factors{I} = model.factors(I);
    else
        for J = find(~cellfun(@iscell,model.factors{I}(:).'))
            model.factors{I}{J} = model.factors{I}(J);
        end
    end
end

% Set functions for saving state, computing the objective function,
% gradient and fast matrix-vector products. New models need only implement
% these four functions.
for I = 1:length(model.factorizations)
    fn = fieldnames(model.factorizations{I});
    model.factorizations{I}.type = find(~strcmp('data',fn));
    model.factorizations{I}.type = fn{model.factorizations{I}.type};
    model.factorizations{I}.factors = ...
        model.factorizations{I}.(model.factorizations{I}.type);
    model.factorizations{I} = ...
        rmfield(model.factorizations{I},model.factorizations{I}.type);
    switch model.factorizations{I}.type
        case 'cpd'
            model.factorizations{I}.state = @state_cpd;
            model.factorizations{I}.objfun = @objfun_cpd;
            model.factorizations{I}.grad = @grad_cpd;
            model.factorizations{I}.JHJx = @JHJx_cpd;
        case 'btd'
            model.factorizations{I}.state = @state_btd;
            model.factorizations{I}.objfun = @objfun_btd;
            model.factorizations{I}.grad = @grad_btd;
            model.factorizations{I}.JHJx = @JHJx_btd;
        case 'regL2'
            model.factorizations{I}.state = @state_regL2;
            model.factorizations{I}.objfun = @objfun_regL2;
            model.factorizations{I}.grad = @grad_regL2;
            model.factorizations{I}.JHJx = @JHJx_regL2;
        case 'regL1'
            model.factorizations{I}.state = @state_regL1;
            model.factorizations{I}.objfun = @objfun_regL1;
            model.factorizations{I}.grad = @grad_regL1;
            model.factorizations{I}.JHJx = @JHJx_regL1;
        otherwise
            error('sdf_nls:model','Model %s is not supported', ...
                model.factorizations{I}.type);
    end
end

% Fill in constants and dereference pointers from factors to variables.
cache.factors.isconst = cell(1,length(model.factors));
for I = 1:length(model.factors)
    cache.factors.isconst{I} = ...
        cellfun(@(f)isnumeric(f{1})&&~isscalar(f{1}),model.factors{I});
    for J = 1:numel(model.factors{I})
        if cache.factors.isconst{I}(J)
            const = model.factors{I}{J}{1};
            for K = 2:length(model.factors{I}{J})
                const = model.factors{I}{J}{K}(const);
            end
            model.factors{I}{J} = {const};
        elseif ischar(model.factors{I}{J}{1})
            model.factors{I}{J}{1} = ...
                find(strcmp(model.factors{I}{J}{1},cache.variables.names));
        end
    end
end

% Dereference and flatten pointers from factorizations to factors.
for I = 1:length(model.factorizations)
	model.factorizations{I}.factors = ...
        deref(model.factorizations{I}.factors);
end
function array = deref(array)
    for i = 1:length(array)
        if iscell(array{i})
            array{i} = deref(array{i});
        elseif ischar(array{i})
            array{i} = find(strcmp(array{i},cache.factors.names));
        end
    end
end

% Convert full data to incomplete/sparse format where appropriate.
cache.factorizations.isincomplete = false(1,length(model.factorizations));
cache.factorizations.issparse = false(1,length(model.factorizations));
for I = 1:length(model.factorizations)
    if ~isfield(model.factorizations{I},'data') || ...
       iscell(model.factorizations{I}.data)
        continue;
    end
    model.factorizations{I}.data = ...
        fmt(model.factorizations{I}.data,true);
    if isstruct(model.factorizations{I}.data)
        cache.factorizations.isincomplete(I) = ...
            model.factorizations{I}.data.incomplete;
        cache.factorizations.issparse(I) = ...
            model.factorizations{I}.data.sparse;
    end
end

% Cache expansion of variables into factors.
cache.factors.expanded = cell(size(model.factors));
cache.factors.sequence = cell(size(model.factors));
cache.factors.state = cell(size(model.factors));
for I = 1:length(model.factors)
    cache.factors.sequence{I} = cell(size(model.factors{I}));
    cache.factors.state{I} = cell(size(model.factors{I}));
    for J = 1:numel(model.factors{I})
        cache.factors.sequence{I}{J} = cell(1,length(model.factors{I}{J}));
        cache.factors.state{I}{J} = cell(1,length(model.factors{I}{J})-1);
    end
end
expand(model.variables,'cache');

% Cache factors' structure.
cache.factors.structure = ...
    cellfun(@structure,cache.factors.sequence,'UniformOutput',false);

% Cache offsets for variables and expanded factors.
cache.variables.offset = ...
    cumsum([1 cellfun(@(v)numel(serialize(v)),model.variables)]);
cache.factorizations.serialized = cell(size(model.factorizations));
cache.factorizations.offset = cell(size(model.factorizations));
cache.factorizations.suboffset = cell(size(model.factorizations));
for I = 1:length(model.factorizations)
    cache.factorizations.serialized{I} = ...
        serialize(model.factorizations{I}.factors).';
    cache.factorizations.offset{I} = ...
        cumsum([1 cellfun(@(v)numel(serialize(v)), ...
        cache.factors.expanded(cache.factorizations.serialized{I}))]);
    cum = cumsum([0 cellfun(@(f)sum(serialize( ...
        cellfun(@(s)numel(s{end}),f))),cache.factors.sequence( ...
        cache.factorizations.serialized{I}))]);
    fct = cellfun(@(f)cellfun(@(s)false(size(s{end})),f,'UniformOutput',...
        false),cache.factors.sequence( ...
        cache.factorizations.serialized{I}),'UniformOutput',false);
    cache.factorizations.suboffset{I} = cell(1,length(fct));
    for J = 1:length(fct)
        cache.factorizations.suboffset{I}{J} = cell(size(fct{J}));
        for K = 1:numel(fct{J})
            subfct = fct{J};
            subfct{K} = true(size(subfct{K}));
            subfct = find(cell2mat(subfct));
            if all(subfct(2:end) == subfct(1:end-1)+1)
                subfct = [subfct(1) subfct(end)];
            end
            cache.factorizations.suboffset{I}{J}{K} = subfct+cum(J);
        end
    end
end

% Set (relative) weights.
if isfield(options,'RelWeights') && isfield(options,'Weights')
    warning('sdf_nls:weights',['Both relative and absolute weights ' ...
        'are supplied, proceeding with absolute weights.']);
end
if ~isfield(options,'RelWeights')
    options.RelWeights = ones(1,length(model.factorizations));
end
options.RelWeights = options.RelWeights./sum(options.RelWeights);
if ~isfield(options,'Weights')
    for I = 1:length(model.factorizations)
        if isfield(model.factorizations{I},'data') && ...
           ~iscell(model.factorizations{I}.data)
            if isstruct(model.factorizations{I}.data)
                NUM = numel(model.factorizations{I}.data.val);
            else
                NUM = numel(model.factorizations{I}.data);
            end
        else
            NUM = cache.factorizations.offset{I}(end);
        end
        options.Weights(I) = 2*options.RelWeights(I)/NUM;
    end
end
if length(options.Weights) ~= length(model.factorizations)
    error('sdf_nls:weights',['The number of weights must equal the ' ...
        'number of factorizations.']);
end

% Initialize each factorization's state.
cache.factorizations.state = cell(size(model.factorizations));
for I = 1:length(model.factorizations)
    model.factorizations{I}.state(I,true);
end

% Run optimization algorithm.
dF.JHF = @grad;
dF.JHJx = @JHJx;
if options.PC && isfield(options,'M') && isa(options.M,'function_handle')
    dF.M = options.M;
end
[z,output] = options.Algorithm(@objfun,dF,model.variables,options);
output.Name = func2str(options.Algorithm);

% Return output.
if isfield(cache.variables,'names')
    sol.variables = cell2struct(z(:),cache.variables.names);
else
    sol.variables = z;
end
if isfield(cache.factors,'names')
    sol.factors = cell2struct(reshape(expand(z),[],1),cache.factors.names);
else
    sol.factors = expand(z);
end

function x = expand(z,where)
% Expands the variables into factors stored in the field factors.expanded,
% and stores their computational state in the field factors.state.

    % Save some references for speed.
    factors = model.factors;
    isconst = cache.factors.isconst;

    % Expand into cache or into the output variable x.
    if nargin == 2 && ischar(where)

        % For each ith factor...
        for i = 1:length(factors)
            % For each jth subfactor...
            for j = 1:numel(factors{i})
                if isconst{i}(j)
                    % Constant subfactor.
                    cache.factors.sequence{i}{j}{1} = factors{i}{j}{1};
                else
                    % For each kth transformation...
                    cache.factors.sequence{i}{j}{1} = z{factors{i}{j}{1}};
                    for k = 2:length(factors{i}{j})
                        [cache.factors.sequence{i}{j}{k}, ...
                            cache.factors.state{i}{j}{k-1}] = ...
                            factors{i}{j}{k}( ...
                            cache.factors.sequence{i}{j}{k-1},[]);
                    end
                end
            end
            if numel(cache.factors.sequence{i}) == 1
                cache.factors.expanded{i} = ...
                    cache.factors.sequence{i}{j}{end};
            else
                cache.factors.expanded{i} = ...
                    cell2mat(cellfun(@(f)f{end}, ...
                    cache.factors.sequence{i},'UniformOutput',false));
            end
        end

    else

        % For each ith factor...
        x = cell(size(factors));
        for i = 1:length(factors)
            % For each jth subfactor...
            sub = cell(size(factors{i}));
            for j = 1:numel(factors{i})
                if isconst{i}(j)
                    % Constant subfactor.
                    sub{j} = factors{i}{j}{1};
                else
                    % For each kth transformation...
                    sub{j} = z{factors{i}{j}{1}};
                    for k = 2:length(factors{i}{j})
                        sub{j} = factors{i}{j}{k}(sub{j},[]);
                    end
                end
            end
            if numel(sub) == 1, x{i} = sub{1};
            else x{i} = cell2mat(sub);
            end
        end

    end

end

function x = derivexpand(r)
% Linearly expand variables using their transformations' Jacobians.

    % Save some references for speed.
    factors = model.factors;
    voffset = cache.variables.offset;
    isconst = cache.factors.isconst;
    sequence = cache.factors.sequence;
    state = cache.factors.state;
    structure = cache.factors.structure;

    % For each ith factor...
    x = cell(size(factors));
    for i = 1:length(factors)
        % For each jth subfactor...
        sub = cell(size(factors{i}));
        for j = 1:numel(factors{i})
            if isconst{i}(j)
                % Constant subfactor.
                sub{j} = zeros(size(sequence{i}{j}{1}));
            else
                % For each kth transformation...
                seq = factors{i}{j};
                ftr = sequence{i}{j};
                stt = state{i}{j};
                dim = structure{i}{j}{1};
                sub{j} = r(voffset(seq{1}):voffset(seq{1}+1)-1);
                if isnumeric(dim)
                    if ~isempty(dim), sub{j} = reshape(sub{j},dim); end
                else sub{j} = deserialize(sub{j},dim);
                end
                for k = 2:length(seq)
                    task = stt{k-1};
                    task.l = [];
                    task.r = sub{j};
                    sub{j} = seq{k}(ftr{k-1},task);
                end
            end

        end
        if numel(sub) == 1, x{i} = sub{1};
        else x{i} = cell2mat(sub);
        end
    end

end

function y = derivcontract(f,y,l)
% Linearly contract factors using their transformations' Jacobians.

    % Save some references for speed.
    factors = model.factors;
    voffset = cache.variables.offset;
    isconst = cache.factors.isconst;
    sequence = cache.factors.sequence;
    state = cache.factors.state;
    structure = cache.factors.structure;
    serialized = cache.factorizations.serialized{f};
    soffset = cache.factorizations.suboffset{f};

    % Update y with Jacobian-contracted variables.
    for i = 1:length(serialized)
        % For each jth subfactor...
        idx = serialized(i);
        for j = 1:numel(factors{idx})

            % Skip this subfactor if it's constant.
            if isconst{idx}(j), continue; end

            % Apply sequence of Jacobian-vector products.
            seq = factors{idx}{j};
            ftr = sequence{idx}{j};
            stt = state{idx}{j};
            if size(soffset{i}{j},2) == 2
                sub = l(soffset{i}{j}(1):soffset{i}{j}(2));
            else
                sub = l(soffset{i}{j});
            end
            if length(seq) > 1
                dim = structure{idx}{j}{end};
                if isnumeric(dim)
                    if ~isempty(dim), sub = reshape(sub,dim); end
                else sub = deserialize(sub,dim);
                end
                for k = length(seq):-1:2
                    task = stt{k-1};
                    task.l = sub;
                    task.r = [];
                    sub = seq{k}(ftr{k-1},task);
                end
                sub = serialize(sub);
            end

            % Update y.
            jdx = voffset(seq{1}):voffset(seq{1}+1)-1;
            y(jdx) = y(jdx)+sub;

        end
    end

end

function z = getfactors(f,x)
% Retrieves factors of the fth factorization, given the factors x.

    % Deserialize the factors.
    z = deserialize_local(x,model.factorizations{f}.factors);
    function z = deserialize_local(x,dim)
        z = cell(size(dim));
        for i = 1:numel(dim)
            if iscell(dim{i})
                z{i} = deserialize_local(x,dim{i});
            else
                z{i} = x{dim{i}};
            end
        end
    end

end

function fval = objfun(z)

    % Expand the variables into factors.
    x = expand(z);

    % Compute objective function value.
    fval = 0;
    for i = find(options.Weights(:).' ~= 0)
        fval = fval+options.Weights(i)* ...
            model.factorizations{i}.objfun(i,getfactors(i,x));
    end

end

function grad = grad(z)

    % Expand the variables into factors.
    expand(z,'cache');

    % Let each factorization save intermediate computations in the cache.
    for i = 1:length(model.factorizations)
        model.factorizations{i}.state(i);
    end

    % Compute the gradient.
    grad = zeros(cache.variables.offset(end)-1,1);
    for i = find(options.Weights(:).' ~= 0)

        % Compute model's gradient.
        tmp = getfactors(i,cache.factors.expanded);
        JHF = model.factorizations{i}.grad(i,tmp);

        % Contract the factor matrices into variables.
        if i == 1
            grad = options.Weights(i)*derivcontract(i,grad,JHF);
        else
            grad = grad+options.Weights(i)*derivcontract(i,grad,JHF);
        end

    end

end

function y = JHJx(~,x)

    % Expand the variables into factors.
    x = derivexpand(x);

    % Compute (J(z)'*J(z))*x after expansion.
    y = zeros(cache.variables.offset(end)-1,1);
    for i = find(options.Weights(:).' ~= 0)

        % Apply fast matrix-vector product.
        tmpa = getfactors(i,cache.factors.expanded);
        tmpb = getfactors(i,x);
        JHJx = model.factorizations{i}.JHJx(i,tmpa,tmpb);

        % Contract the factor matrices into variables.
        if i == 1
            y = options.Weights(i)*derivcontract(i,y,JHJx);
        else
            y = y+options.Weights(i)*derivcontract(i,y,JHJx);
        end

    end

end

function [z,offset] = deserialize(z,dim,offset)
    if iscell(dim)
        v = z;
        z = cell(size(dim));
        if nargin < 3, offset = 0; end
        for i = 1:numel(z)
            if iscell(dim{i})
                [z{i},offset] = deserialize(v,dim{i},offset);
            else
                n = prod(dim{i}(:));
                z{i} = reshape(v(offset+(1:n)),dim{i});
                offset = offset+n;
            end
        end
    elseif ~isempty(dim)
        z = reshape(z,dim);
    end
end

function z = serialize(z)
    if iscell(z)
        for i = find(cellfun(@iscell,z(:).'))
            z{i} = serialize(z{i});
        end
        s = cellfun(@numel,z(:)); o = [0; cumsum(s)];
        c = z; z = zeros(o(end),1);
        for i = 1:length(s), z(o(i)+(1:s(i))) = c{i}(:); end
    else
        z = z(:);
    end
end

function dim = structure(z)
    if iscell(z)
        dim = cellfun(@size,z,'UniformOutput',false);
        for i = find(cellfun(@iscell,z(:).'))
            dim{i} = structure(z{i});
        end
    else
        dim = size(z);
        if numel(z) == dim(1), dim = []; end
    end
end

% Model: canonical polyadic decomposition ---------------------------------

function state_cpd(f,firstrun)
% Can read from model, cache and options, can save state by writing to the
% structure cache.factorizations.state{f}.

    if nargin == 2 && firstrun

        % Store the fraction of known elements.
        if cache.factorizations.isincomplete(f)
            cache.factorizations.scale{f} = ...
                length(model.factorizations{f}.data.val)./...
                prod(model.factorizations{f}.data.size);
        end

        % Set Block-Jacobi preconditioner if
        % - the model is a single CPD and
        % - each factor consists of exactly one subfactor and
        % - no subfactor is transformed and
        % - every subfactor is a reference to a variable and
        % - all variable references are unique and
        % - all factor references are unique and
        % - an NLS algorithm is used.
        options.BJ = false;
        var = cellfun(@(v)v{1}{1},model.factors,'UniformOutput',false);
        ftr = cell2mat(model.factorizations{f}.factors);
        if options.PC && numel(model.factorizations) == 1 && ...
                all(cellfun(@(v)length(v) == 1,model.factors)) && ...
                all(cellfun(@(v)length(v{1}) == 1,model.factors)) && ...
                all(cellfun(@isscalar,var)) && ...
                length(unique(cell2mat(var))) == length(cell2mat(var))&&...
                length(unique(ftr)) == length(ftr) && ...
                ~isempty(strfind(func2str(options.Algorithm),'nls'))
            options.M = @pc_cpd;
            options.BJ = true;
        end

    end

    % Cache the factor matrices' Gramians.
    idx = cache.factorizations.serialized{f};
    N = length(idx);
    R = size(cache.factors.expanded{idx(1)},2);
    cache.factorizations.state{f}.UHU = zeros(N,R*R);
    for n = 1:N
        tmp = cache.factors.expanded{idx(n)};
        tmp = conj(tmp'*tmp);
        cache.factorizations.state{f}.UHU(n,:) = tmp(:);
    end

    % Optionally cache the inverses of the Gramians for the preconditioner.
    % In a faster language, this should be the Cholesky factor instead.
    if options.BJ
        cache.factorizations.state{f}.invW = cell(1,N);
        for n = 1:N
            tmp = cache.factorizations.state{f}.UHU([1:n-1 n+1:N],:);
            if N > 2, tmp = prod(tmp,1); end
            cache.factorizations.state{f}.invW{n} = ...
                inv(reshape(tmp,[R R]));
        end
    end

end

function fval = objfun_cpd(f,z)

    % CPD objective function.
    T = model.factorizations{f}.data;
    isincomplete = cache.factorizations.isincomplete(f);
    issparse = cache.factorizations.issparse(f);
    if ~isincomplete || length(T.ind)/prod(T.size) > options.TolLargeScale
        fval = z{1}*kr(z(end:-1:2)).';
        if isincomplete, fval = fval(T.ind)-T.val;
        elseif issparse
            if ~isempty(T.ind), fval(T.ind) = fval(T.ind)-T.val; end
        else fval = fval-reshape(T,size(fval));
        end
    else
        fval = -T.val;
        for r = 1:size(z{1},2)
            tmp = z{1}(T.sub{1},r);
            for n = 2:length(z), tmp = tmp.*z{n}(T.sub{n},r); end
            fval = fval+tmp;
        end
    end
    if isincomplete
        T.val = fval;
        if ~isempty(T.matrix)
            T.matrix = sparse(double(T.sub{1}), ...
                double(1+idivide(T.ind-1,int64(size(T.matrix,1)))), ...
                double(fval),size(T.matrix,1),size(T.matrix,2));
        end
        cache.factorizations.residual{f} = T;
    else
        if issparse, size_tens = T.size;
        else size_tens = size(T); end
        cache.factorizations.residual{f} = reshape(fval,size_tens);
    end
    fval = 0.5*(fval(:)'*fval(:));

end

function grad = grad_cpd(f,z)

    % CPD scaled conjugate cogradient.
    E = cache.factorizations.residual{f};
    offset = cache.factorizations.offset{f};
    grad = zeros(offset(end)-1,1);
    for n = 1:length(z)
        tmp = full(mtkrprod(E,z,n));
        grad(offset(n):offset(n+1)-1) = tmp(:);
    end

end

function y = JHJx_cpd(f,z,x)

    % CPD fast Jacobian's Gramian vector product.
    % Ignores the fact that the tensor might be incomplete.
    R = size(z{1},2);
    N = length(z);
    offset = cache.factorizations.offset{f};
    UHU = cache.factorizations.state{f}.UHU;
    XHU = zeros(R,R,N);
    y = zeros(offset(end)-1,1);
    for n = 1:N
        Wn = UHU([1:n-1 n+1:N],:);
        if N > 2, Wn = prod(Wn,1); end
        XHU(:,:,n) = conj(x{n}'*z{n});
        y(offset(n):offset(n+1)-1) = x{n}*reshape(Wn,[R R]);
    end
    for n = 1:N-1
        idxn = offset(n):offset(n+1)-1;
        Wn = zeros(R);
        for m = n+1:N
            idxm = offset(m):offset(m+1)-1;
            if N == 2
                Wn = Wn+XHU(:,:,m);
                JHJmnx = z{m}*XHU(:,:,n);
            else
                Wnm = UHU([1:n-1 n+1:m-1 m+1:N],:);
                if N > 3, Wnm = prod(Wnm,1); end
                Wnm = reshape(Wnm,[R R]);
                Wn = Wn+Wnm.*XHU(:,:,m);
                JHJmnx = z{m}*(Wnm.*XHU(:,:,n));
            end
            y(idxm) = y(idxm)+JHJmnx(:);
        end
        JHJnx = z{n}*Wn;
        y(idxn) = y(idxn)+JHJnx(:);
    end

    % If incomplete, approximate the effect of missing entries.
    if cache.factorizations.isincomplete(f)
        y = y*cache.factorizations.scale{f};
    end

end

function x = pc_cpd(~,b)

    % Solve M*x = b, where M is a block-diagonal approximation for JHJ.
    % Equivalent to simultaneous ALS updates for each of the factors.
    x = zeros(size(b));
    offset = cache.factorizations.offset{1};
    invW = cache.factorizations.state{1}.invW;
    for n = 1:length(offset)-1
        idx = offset(n):offset(n+1)-1;
        tmp = reshape(b(idx),[],size(invW{1},1))*invW{n};
        x(idx) = tmp(:);
    end
    x = x/options.Weights;

    % If incomplete, approximate the effect of missing entries.
    if cache.factorizations.isincomplete(1)
        x = x/cache.factorizations.scale{1};
    end

end

% Model: block term decomposition -----------------------------------------

function state_btd(f,firstrun)
% Can read from model, cache and options, can save state by writing to the
% structure cache.factorizations.state{f}.

    if nargin == 2 && firstrun

        % Store the fraction of known elements.
        if cache.factorizations.isincomplete(f)
            cache.factorizations.scale{f} = ...
                length(model.factorizations{f}.data.val)./...
                prod(model.factorizations{f}.data.size);
        end

        % Set Block-Jacobi preconditioner if
        % - the model is a single BTD and
        % - each factor consists of exactly one subfactor and
        % - no subfactor is transformed and
        % - every subfactor is a reference to a variable and
        % - all variable references are unique and
        % - all factor references are unique and
        % - an NLS algorithm is used.
        options.BJ = false;
        var = cellfun(@(v)v{1}{1},model.factors,'UniformOutput',false);
        ftr = cache.factorizations.serialized{1};
        if options.PC && numel(model.factorizations) == 1 && ...
                all(cellfun(@(v)length(v) == 1,model.factors)) && ...
                all(cellfun(@(v)length(v{1}) == 1,model.factors)) && ...
                all(cellfun(@isscalar,var)) && ...
                length(unique(cell2mat(var))) == length(cell2mat(var))&&...
                length(unique(ftr)) == length(ftr) && ...
                ~isempty(strfind(func2str(options.Algorithm),'nls'))
            options.M = @pc_btd;
            options.BJ = true;
        end

    end

    % Cache the factor matrices' Gramians.
    U = getfactors(f,cache.factors.expanded);
    R = length(U);
    N = length(U{1})-1;
    [idx,jdx,kdx] = ndgrid(1:R,1:N,1:R);
    cache.factorizations.state{f}.UHU = ...
        arrayfun(@(i,n,j)U{i}{n}'*U{j}{n}, ...
        idx,jdx,kdx,'UniformOutput',false);

    % Optionally cache some results for the block-Jacobi preconditioner.
    if options.BJ
        [idx,jdx] = ndgrid(1:R,1:N);
        UHU = cache.factorizations.state{f}.UHU;
        cache.factorizations.state{f}.invSKS = arrayfun( ...
            @(r,n)inv(mtkronprod(U{r}{end},UHU(r,:,r),n)* ...
            conj(reshape(permute(U{r}{end},[1:n-1 n+1:N n]), ...
            [],size(U{r}{end},n)))),idx,jdx,'UniformOutput',false);
        cache.factorizations.state{f}.invUHU = arrayfun( ...
            @(r,n)inv(UHU{r,n,r}),idx,jdx,'UniformOutput',false);
    end

end

function fval = objfun_btd(f,z)

    % BTD objective function.
    T = model.factorizations{f}.data;
    isincomplete = cache.factorizations.isincomplete(f);
    issparse = cache.factorizations.issparse(f);
    if ~isincomplete || length(T.ind)/prod(T.size) > options.TolLargeScale
        fval = z{1}{1}*mtkronprod(z{1}{end},z{1}(1:end-1),1,'H');
        for r = 2:length(z)
            fval = fval+z{r}{1}*mtkronprod(z{r}{end},z{r}(1:end-1),1,'H');
        end
        if isincomplete, fval = fval(T.ind)-T.val;
        elseif issparse
            if ~isempty(T.ind), fval(T.ind) = fval(T.ind)-T.val; end
        else fval = fval-reshape(T,size(fval));
        end
    else
        fval = -T.val;
        for r = 1:length(z)
            size_core = cellfun('size',z{r}(1:end-1),2);
            idx = cell(1,length(size_core));
            S = z{r}{end};
            for i = 1:numel(S)
                [idx{:}] = ind2sub(size_core,i);
                tmp = S(idx{:})*z{r}{1}(T.sub{1},idx{1});
                for n = 2:length(size_core)
                    tmp = tmp.*z{r}{n}(T.sub{n},idx{n});
                end
                fval = fval+tmp;
            end
        end
    end
    if isincomplete
        T.val = fval;
        if ~isempty(T.matrix)
            T.matrix = sparse(double(T.sub{1}), ...
                double(1+idivide(T.ind-1,int64(size(T.matrix,1)))), ...
                double(fval),size(T.matrix,1),size(T.matrix,2));
        end
        cache.factorizations.residual{f} = T;
    else
        if issparse, size_tens = T.size;
        else size_tens = size(T); end
        cache.factorizations.residual{f} = reshape(fval,size_tens);
    end
    fval = 0.5*(fval(:)'*fval(:));

end

function grad = grad_btd(f,z)

    % BTD scaled conjugate cogradient.
    N = length(z{1})-1;
    E = cache.factorizations.residual{f};
    offset = cache.factorizations.offset{f};
    grad = zeros(offset(end)-1,1);
    cnt = 1;
    for r = 1:length(z)
        U = z{r}(1:N);
        S = conj(z{r}{end});
        for n = 1:N
            tmp = full(mtkronprod(E,U,n))* ...
                reshape(permute(S,[1:n-1 n+1:N n]),[],size(S,n));
            grad(offset(cnt):offset(cnt+1)-1) = tmp(:);
            cnt = cnt+1;
        end
        tmp = full(mtkronprod(E,U,0));
        grad(offset(cnt):offset(cnt+1)-1) = tmp;
        cnt = cnt+1;
    end

end

function y = JHJx_btd(f,z,x)

    % BTD fast Jacobian's Gramian vector product.
    % Ignores the fact that the tensor might be incomplete.
    R = length(z);
    N = length(z{1})-1;
    offset = cache.factorizations.offset{f};
    UHU = cache.factorizations.state{f}.UHU;
    [idx,jdx,kdx] = ndgrid(1:R,1:N,1:R);
    XHU = arrayfun(@(i,n,j)x{i}{n}'*z{j}{n}, ...
        idx,jdx,kdx,'UniformOutput',false);
    y = zeros(offset(end)-1,1);
    cnt = 1;
    for ri = 1:R

        % Factor matrices.
        for ni = 1:N
            idx = offset(cnt):offset(cnt+1)-1;
            Sri = permute(z{ri}{end},[1:ni-1 ni+1:N ni]);
            Sri = conj(reshape(Sri,[],size(Sri,N)));
            for rj = 1:R
                Srj = z{rj}{end};
                tmp = mtkronprod(x{rj}{end},UHU(rj,:,ri),ni);
                for nj = [1:ni-1 ni+1:N]
                    proj = UHU(rj,:,ri);
                    proj{nj} = XHU{rj,nj,ri};
                    tmp = tmp+mtkronprod(Srj,proj,ni);
                end
                tmp = z{rj}{ni}*(tmp*Sri);
                tmp = tmp+x{rj}{ni}* ...
                    (mtkronprod(z{rj}{end},UHU(rj,:,ri),ni)*Sri);
                if rj == 1, y(idx) = tmp(:);
                else y(idx) = y(idx)+tmp(:); end
            end
            cnt = cnt+1;
        end

        % Core tensor.
        idx = offset(cnt):offset(cnt+1)-1;
        for rj = 1:R
            Srj = z{rj}{end};
            tmp = mtkronprod(x{rj}{end},UHU(rj,:,ri),0);
            for nj = 1:N
                proj = UHU(rj,:,ri);
                proj{nj} = XHU{rj,nj,ri};
                tmp = tmp+mtkronprod(Srj,proj,0);
            end
            if rj == 1, y(idx) = tmp(:);
            else y(idx) = y(idx)+tmp(:); end
        end
        cnt = cnt+1;

    end

    % If incomplete, approximate the effect of missing entries.
    if cache.factorizations.isincomplete(f)
        y = y*cache.factorizations.scale{f};
    end

end

function x = pc_btd(~,b)

    % Solve M*x = b, where M is a block-diagonal approximation for JHJ.
    x = zeros(size(b));
    R = length(model.factorizations{1}.factors);
    N = length(model.factorizations{1}.factors{1})-1;
    offset = cache.factorizations.offset{1};
    invSKS = cache.factorizations.state{1}.invSKS;
    invUHU = cache.factorizations.state{1}.invUHU;
    cnt = 1;
    for r = 1:R
        for n = 1:N
            idx = offset(cnt):offset(cnt+1)-1;
            tmp = reshape(b(idx),[],size(invSKS{r,n},1))*invSKS{r,n};
            x(idx) = tmp(:);
            cnt = cnt+1;
        end
        idx = offset(cnt):offset(cnt+1)-1;
        size_core = cellfun('size',invUHU(r,:),1);
        x(idx) = mtkronprod(reshape(b(idx),size_core),invUHU(r,:),0);
        cnt = cnt+1;
    end
    x = x/options.Weights;

    % If incomplete, approximate the effect of missing entries.
    if cache.factorizations.isincomplete(1)
        x = x/cache.factorizations.scale{1};
    end

end

% Model: L2 regularization ------------------------------------------------

function state_regL2(f,firstrun)

    % Format right hand side, if available.
    if nargin == 2 && firstrun
         cache.factorizations.hasdata(f) = ...
             isfield(model.factorizations{f},'data');
         if cache.factorizations.hasdata(f) && ...
            ~iscell(model.factorizations{f}.data)
             model.factorizations{f}.data = ...
                 {model.factorizations{f}.data};
         end
         if cache.factorizations.hasdata(f)
             model.factorizations{f}.data = cellfun(@(d)full(d), ...
                 model.factorizations{f}.data,'UniformOutput',false);
         end
    end

end

function fval = objfun_regL2(f,z)
    % L2 regularization objective function.
    fval = 0;
    for i = 1:length(z)
        e = z{i}(:);
        if cache.factorizations.hasdata(f)
            e = e-model.factorizations{f}.data{i}(:);
        end
        fval = fval+(e'*e);
    end
    fval = 0.5*fval;
end

function grad = grad_regL2(f,z)
    % L2 regularization scaled conjugate cogradient.
    if cache.factorizations.hasdata(f)
        grad = cell2mat(cellfun(@(f,d)f(:)-d(:), ...
            z(:),model.factorizations{f}.data(:),'UniformOutput',false));
    else
        grad = cell2mat(cellfun(@(f)f(:),z(:),'UniformOutput',false));
    end
end

function y = JHJx_regL2(~,~,x)
    % L2 regularization Hessian vector product.
    y = cell2mat(cellfun(@(f)f(:),x(:),'UniformOutput',false));
end

% Model: L1 regularization ------------------------------------------------

function state_regL1(f,firstrun)

    % Format right hand side, if available.
    if nargin == 2 && firstrun
         cache.factorizations.hasdata(f) = ...
             isfield(model.factorizations{f},'data');
         if cache.factorizations.hasdata(f) && ...
            ~iscell(model.factorizations{f}.data)
             model.factorizations{f}.data = ...
                 {model.factorizations{f}.data};
         end
         if cache.factorizations.hasdata(f)
             model.factorizations{f}.data = cellfun(@(d)full(d), ...
                 model.factorizations{f}.data,'UniformOutput',false);
         end
         if ~isfield(options,'mu')
             if cache.factorizations.hasdata(f)
                m = max(abs(serialize(model.factorizations{f}.data)));
             else
                m = 0;
             end
             options.Mu = max(m,1)/100;
         end
    end

end

function fval = objfun_regL1(f,z)
    % Approximate L1 regularization objective function.
    fval = 0;
    for i = 1:length(z)
        e = z{i}(:);
        if cache.factorizations.hasdata(f)
            e = e-model.factorizations{f}.data{i}(:);
        end
        far = abs(e) > options.Mu;
        e(far) = abs(e(far));
        e2 = e(~far);
        e2 = e2.*conj(e2);
        e(~far) = 2*options.Mu*e2./(e2+options.Mu^2);
        fval = fval+sum(e);
    end
    fval = 0.5*fval;
end

function grad = grad_regL1(f,z)
    % Approximate L1 regularization scaled conjugate cogradient.
    if cache.factorizations.hasdata(f)
        grad = cell2mat(cellfun(@(f,d)f(:)-d(:), ...
            z(:),model.factorizations{f}.data(:),'UniformOutput',false));
    else
        grad = cell2mat(cellfun(@(f)f(:),z(:),'UniformOutput',false));
    end
    far = abs(grad) > options.Mu;
    grad(far) = grad(far)./(2*abs(grad(far)));
    grad(~far) = 2*options.Mu^3*grad(~far)./ ...
        (grad(~far).*conj(grad(~far))+options.Mu^2).^2;
end

function y = JHJx_regL1(f,z,x)
    % Approximate L1 regularization Hessian vector product.
    if cache.factorizations.hasdata(f)
        hess = cell2mat(cellfun(@(f,d)f(:)-d(:), ...
            z(:),model.factorizations{f}.data(:),'UniformOutput',false));
    else
        hess = cell2mat(cellfun(@(f)f(:),z(:),'UniformOutput',false));
    end
    if ~isreal(hess)
        error('sdf_nls:regL1',['Please use sdf_minf when applying L1 ' ...
            'regularization on complex factors.']);
    end
    far = abs(hess) > options.Mu;
    hess(far) = 0;
    if any(hess)
        hess(~far) = 2*options.Mu^3*(options.Mu^2-3*hess(~far).^2)./ ...
            (options.Mu^2+hess(~far).^2).^3;
        y = hess.*cell2mat(cellfun(@(f)f(:),x(:),'UniformOutput',false));
    else
        y = hess;
    end
end

end

author: andrewssobral

repo: lrslibrary


text:
function [p] = dot(qt1,qt2, do_qr)
%Dot product of two QTT-Tuckers
%   [PR]=DOT(TT1,TT2) -- dot product of two TT-tensors
%
%   [PR]=DOT(TT1,TT2, DO_QR) if DO_QR==true is specified, perform the
%   left-to-right QRs of TT1,TT2
%   before the scalar product. It increases the  accuracy in some cases.
%
% In general, returns a 4D tensor of sizes
% r0(tt1), r0(tt2), rd(tt1), rd(tt2)
% If r0(tt1) = r0(tt2) = 1 it returns a matrix of size rd(tt1) x rd(tt2)
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

if (nargin<3)||(isempty(do_qr))
    do_qr = false;
end;

d=qt1.dphys;
for i=1:d
    % dot product of factors gives rt1 times rt2 matrix, to be convolved
    % between cores
    % Make it more precise by QRs? <-- already inside @tt_tensor/dot
%     [qt1.tuck{i}, rv1]=qr(qt1.tuck{i}, 'lr');
%     [qt2.tuck{i}, rv2]=qr(qt2.tuck{i}, 'lr');
%     rv1 = eye(qt1.tuck{i}.r(end));
%     rv2 = eye(qt2.tuck{i}.r(end));
    Pfac = squeeze(dot(qt1.tuck{i}, qt2.tuck{i}, do_qr)); % size rt1,rt2
%     rt1 = size(rv1, 2); rt2 = size(rv2, 2);
%     rt1new = size(rv1,1); rt2new = size(rv2,1);
    % Now, merge Pfac to the core of qt1
    curcr = qt2.core{i};
    rc1 = size(curcr, 1); rt2 = size(curcr, 2); rc2 = size(curcr, 3);
    curcr = permute(curcr, [1, 3, 2]);
    curcr = reshape(curcr, rc1*rc2, rt2);
%     curcr = curcr*(rv1.');
    curcr = curcr*(Pfac.'); % Now, core2 has the tucker ranks of qt1
    curcr = reshape(curcr, rc1, rc2, qt1.core.n(i));
    qt2.core{i} = permute(curcr, [1, 3, 2]);

%     curcr = qt2.core{i};
%     rc1 = size(curcr, 1); rc2 = size(curcr, 3);
%     curcr = permute(curcr, [1, 3, 2]);
%     curcr = reshape(curcr, rc1*rc2, rt2);
%     curcr = curcr*(rv2.');
%     curcr = reshape(curcr, rc1, rc2, rt2new);
%     qt2.core{i} = permute(curcr, [1, 3, 2]);
end;
% Finaly, dot product of cores. It is consistent, since we merged Pfacs
p = dot(qt1.core, qt2.core, do_qr);
end

author: oseledets

repo: TT-Toolbox


text:
function [init_point, n_alg, initial_terminal_voltage] = initialise_model(param)
% initialise_model performs analytical initialisation of all variables of interest

%   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.

sign_input_density = evaluate_sign_input_density(param); % evaluates sign of input current/power density as per operating mode

T = param.T_init * ones(param.Nsum+param.Nal+param.Ncu,1);

cs_star_p   = param.cs_p_init*ones(param.Np, 1);
cs_star_n   = param.cs_n_init*ones(param.Nn, 1);
cs_star     = [cs_star_p; cs_star_n];

[Phis_p_init,~,Phis_n_init,~] = param.OpenCircuitPotentialFunction(cs_star,T,param,sign_input_density);

Phis_init = [Phis_p_init; Phis_n_init];

if param.OperatingMode==1 || param.OperatingMode==4
    I_density           = param.I_density;
elseif param.OperatingMode==2 || param.OperatingMode==5
    I_density           = param.P_density/(Phis_init(1)-Phis_init(end)); % No need for linear interpolation since starting from equilibrium
else
    I_density = 0; % dummy value for CV mode
end

ce_init = param.ce_init*[ones(param.Np,1);ones(param.Ns,1);ones(param.Nn,1)];
Phie_init = zeros(param.Np + param.Ns + param.Nn, 1); % designated as the ground potential for this system

solverFlux  = zeros(param.Np+param.Nn, 1);
film = zeros(param.Nn, 1);

jflux_init = ionicFlux(ce_init, cs_star, Phis_init, Phie_init, T, solverFlux, film, param,sign_input_density,I_density);

if param.EnableAgeing == 1
    js_init = 0.483e-5*ones(param.Nn, 1); % totally arbitrary/random value for side-reaction flux
else
    js_init = zeros(param.Nn, 1);
end

init_point = [
    jflux_init;...
    Phis_init;...
    Phie_init;...
    js_init;...
    I_density
    ];

n_alg = length(init_point);
initial_terminal_voltage = diff([Phis_n_init(end);Phis_p_init(1)]); % No need for linear interpolation since starting from equilibrium
end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
function y = subsasgn( x, varargin ) %#ok
error( 'cvx tuple objects are read-only.' );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function r = zerofh()
r = @(varargin)0;

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function[teps]=maternedge(varargin)
%MATERNEDGE  Long-time cutoff edge for the Matern impulse response function.
%
%   TE=MATERNEDGE(ALPHA,H,EPSILON) determines a long-time cutoff for the
%   Green's function of a Matern process characterized by slope parameter
%   ALPHA and range parameter H.
%
%   TE is the (approximately) the first time at which the ratio of the
%   time-integrated magnitude of the Green's function, to the value it
%   obtains at infinity, has risen to within EPSILON of unity.
%
%   TE is approximate because of discretization of time.  It will always be
%   greater than the actual time at which the EPSILON level is reached.
%
%   See Lilly et al. (2017) for details.
%
%   MATERNEDGE is a low-level function that is called by MATERNOISE.
%
%   See also MATERNSPEC, MATERNCOV, MATERNIMP, MATERNOISE.
%
%   'maternedge --t' runs some tests.
%
%   Usage: te=maternedge(alpha,h,epsilon);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2013--2017 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--t')
    maternedge_test,return
end

alpha=varargin{1};
h=varargin{2};
epsilon=varargin{3};

arrayify(alpha,h,epsilon);
teps=zeros(size(alpha));

tnorm=[0:0.001:20]';
%tnorm=[0:0.01:20]';
for i=1:length(alpha)

    gint=gammainc(tnorm,alpha(i));  %This is integral divided by its total
    %Note Matlab defines the incomplete gamma function as a ratio to gamma
    ii=find(gint>1-epsilon(i),1,'first');
    %plot(tnorm,gint),hlines(1),hold on
    if ~isempty(ii)
        teps(i)=tnorm(ii)./h(i);
    else
        %disp('The Green''s function does not cross this threshold within 20 e-folding times.')
        teps(i)=inf;
    end
end

if aresame(size(varargin{1}),size(varargin{2}))
    teps=reshape(teps,size(varargin{1},1),size(varargin{1},2));
end

function[]=maternedge_test

alpha=[1 1.5 2 3 4];
%alpha=[3/4 1 1.5 2 3 4];
h=[.01 .02 .05 .2 1];
[alpha,h]=meshgrid(alpha,h);
h=h.*alpha;
epsilon=1e-2;

vcolon(alpha,h);
te=maternedge(alpha,h,epsilon);

t=[1e-6:0.0001:20]';
tnorm=oprod(t,1./h);
g=zeros(size(tnorm));
gint=zeros(size(tnorm));
for i=1:length(alpha)
    d=sqrt(frac(h(i).^(2*alpha(i)-1),materncfun(alpha(i))));
    g(:,i)=frac(1,d)*maternimp(tnorm(:,i),alpha(i),h(i));
    gint(:,i)=cumsum(g(:,i)).*(tnorm(2,i)-tnorm(1,i));
    gint(:,i)=gint(:,i).*h(i).^alpha(i);
end

te2=zeros(size(te));
for i=1:length(alpha)
     ii=find(gint(:,i)>1-epsilon,1,'first');
     if ~isempty(ii)
         te2(i)=tnorm(ii,i);
     else
         te2(i)=nan;
     end
end

%abs(te2-te)./te

reporttest('MATERNEDGE matches explicit calculation to within one percent',allall(abs(te2-te)./te<.01))

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function fg = dtiCreateQuenchStats(fg, statFiberName, statPointName, perPoint, img, fiberSummary, pointCompute,mpRange)
%
% fg = dtiCreateQuenchStats(fg, [statFiberName='avg'], [statPointName='fa'],...
%     [perPoint=1], img, [fiberSummary='avg'],[pointCompute='none']);
%
% Add satistics to a fiber group that are used for Quench visualization and
% queries. Every path gets a number (its statistic) for that path. Some
% paths also have statistics per point along the path (e.g. in the array
% condition). They are stored in the returned fiber group struct.
%
% INPUTS:
%  fg            - Fiber Group. The statistics will be attached to this
%                  structure.
%  statFiberName - Name for a fiber summary statistic. Typically this
%                  will be  something like length or "average FA". THis
%                  statistic is used to limit the range of fibers we see
%                  in the Quench window.
%  statPointName - Name for the per point statistic (e.g., TF, FA,
%                  eccentricity)
%  perPoint      - 1 if computed per point, 0 if only per fiber
%  img           - Data for statistical overlay. The data form is either a
%                  nifti image structure or dti image structure.
%                  If no data are supplied the routine computes a default
%                  statistic determined by the statFiberName. At present,
%                  the only statFiberName that generates a statistic is
%                  'length'.
%  fiberSummary  - If data are supplied, this variable specifies
%                  which statFiber statistic to compute. We currently
%                  support can be 'avg', 'min','max'.  This fiber statistic
%
%  pointCompute  - Method for taking per point color and remapping to new
%                  per point colors.  Either 'none' or 'matchEnd' are in
%                  place.  More methods will arise.  The matchEnd is a
%                  hack.
%
% WEB:
%  mrvBrowseSVN('dtiCreateQuenchStats');
%
%  This function computes output parameters pathStat and point_stat_array.
%  The max and min of pathStat are used to set the range for colormaps etc.
%  in Quench.
%
%Delete this comment soon ---
% Statistics can be one of three types:
%   1) geometry - path length, curvature
%   2) image    - FA value, other pointwise statistic (e.g., eccentricity)
%   3) array    - An array of vectors. each the same length as a
%      corresponding path in the fiber group, and values of the vectors are
%      assigned to each position along the paths.
%
% EXAMPLES:
%
%   Image Stat:
%       fg = mtrImportFibers('fibers.pdb');
%       fa_img = niftiRead('faMap.nii.gz');
%       fg = dtiCreateQuenchStats(fg,'FA_avg','FA', 1, fa_img, 'avg');
%     OR...
%       h = guidata(gcf); fa_img = h.bg(1,2,3 or 5);
%       fg = dtiCreateQuenchStats(fg,'FA_avg','FA', 1, fa_img, 'avg');
%
%   Given Array Stat
%       fg = dtiCreateQuenchStats(fg,'Group',stat_array);
%
%   Geometric Stat
%       fg = dtiCreateQuenchStats(fg,'Length','Length', 1);
%
%
% HISTORY:
% 2009.06.17 : SA wrote it
% 2009.08.04 : AJS fixed the mixing between explicit declaration of
% perPoint and the implicit assumption in most of the code.
% Now we truly have the variable report explicit info.  Also added the
% ability to have a statistic be added by passing the array of per pathway
% stats that you would like to have.
% 2009.10.19 : AJS significantly sped up the code by using cellfun.
%

%% Check if input params are ok
if notDefined('fg') || ~isstruct(fg)
    disp('fg should be a proper fibergroup struct');
    return;
end

if notDefined('statFiberName'), statFiberName = 'avg'; end
if notDefined('statPointName'), statPointName = 'fa'; end
if notDefined('perPoint'), perPoint = 1; end
if notDefined('img')
    if ~strcmpi(statFiberName,'length')
        error('Data for overlay are required');
    end
end
if notDefined('fiberSummary'), fiberSummary = 'avg'; end
if notDefined('pointCompute'), pointCompute = 'none'; end

if ~isfield(fg,'pathwayInfo'), fg.pathwayInfo=[];  end

% Indices into the number of statistical images. Quench counts from 0.  We
% start with Matlab counting from 1, and then we subtract 1 later when we
% do the Quench parameter field.
if ~isfield(fg,'params') || isempty(fg.params),  idx = 1;
else                                             idx = length(fg.params)+1;
end

% If there are only three input variables, the third must be a data set,
% not a name.  So we copy the data set to the statistics array, we set the
% name variable to not applicable, and we do not compute per point.
% Instead, we will be computing per fiber.
% if nargin == 3
%     stat_array = statPointName;
%     statPointName = 'NA';
%     perPoint = 0;
% end

%%
% This statistics parameter field gets attached to the fg.  It contains
% information that Quench uses for visualization and for dynamic queries.
% We should probably write a function for initializing this parameter
% field.
param = struct;              % Initialize structure
param.name = statFiberName;  % Struct name?
param.uid = idx-1;           % Subtract 1 so Quench starts from 0 (C-counting)
param.ile = 1;
param.icpp= perPoint;        % Visualization per point, not per fiber
param.ivs = 1;               %
param.agg = statFiberName;   % A stats name
param.lname = statPointName; % Local Name
param.stat = [];

%dummy.  Same to you.
% c        = struct;
% c.fibers = cell(1);

if notDefined('img')
    % No overlay data, so we use the statFiberName to determine how we color
    % the fibers.

    switch lower(statFiberName)
        case 'length'
            % For each fiber, figure out its length.
            for ii=1:length(fg.fibers)
                fiber = fg.fibers{ii};      % Get the ith fiber.
                %if(size(fiber,2)<size(fiber,1)); fiber=fiber'; end

                % Find the step size - could be be moved out of the loop.  Could
                % use cellfun to find all these.  Mostly these are the same.
                mmPerStep  = sqrt(sum((fiber(:,1)-fiber(:,2)).^2));

                % Find the length for the whole fiber
                param.stat(ii) = mmPerStep*(size(fiber,2)-1);

                % Store the length for this stat image (idx) and the ith pathway.
                % Each fiber group has a pathway info that summarizes all of the
                % stats in stat images.
                fg.pathwayInfo(ii).pathStat(idx) = param.stat(ii);

                % If we are computing per point, then we compute the length to each
                % point on the fiber, so that we don't assign a single length.
                % Rather we assign a length from the end point to each point along
                % the path.  And we put that into the point_stat_array.  Up to now,
                % we had a single length assigned to all of the fibers.
                if perPoint
                    fg.pathwayInfo(ii).point_stat_array(idx,:) = mmPerStep * (0:1:size(fiber,2)-1);
                end
            end
        otherwise
            error('Unknown statFiberName %s\n',statFiberName);
    end

else

    %isNiftiImg = NaN; %mat =[]; data=[];
    minVal=0; maxVal = 1;

    % Test whether the img variable is a nifti file structure or a dti file
    % structure.
    if isNiftiStruct(img);
        % NIFTI land.  Get what you need.
        isNiftiImg = true;
        mat  = img.qto_ijk;
        data = img.data;
    else
        isNiftiImg = false;
        if isfield(img,'img') && isfield(img,'mat')
            % It has the dti fields. So we copy the
            % fields we need and set isNifti to false.
            mat  = inv(img.mat);
            data = img.img;
            if isfield(img,'minVal'),  minVal = img.minVal; end
            if isfield(img,'maxVal'),  maxVal = img.maxVal; end
        else
            error('img must be either nifti or dti');
        end
    end

    % Interpolate from the values in the statistics to the position of the
    % fiber nodes.  For speed, we use nearest neighbor interpolation
    vals = dtiGetValFromFibers(data, fg, mat,[],'nearest');

    % Apply the scaling function function to each entry of vals. This will
    % take a 0,1 range and put it with the scale specified by max and min
    % vals.
    % 'UniformOutput' insists that the return must be a cell array.
    if ~isNiftiImg
        vals = cellfun(@(x) x*(maxVal-minVal)+minVal, vals, 'UniformOutput', false);
    end

    % For each of the cell arrays, we can determine a single value that
    % will be applied to the whole fiber.  This might be the avg, min or
    % max. We compute this value and store it as a single scalar here in
    % the stat field.  The stat field should have one scalar for every
    % fiber in the fiber group.  This is for the non per-point calculation.
    %
    switch(lower(fiberSummary))
        case 'avg'
            param.stat = cellfun(@mean,vals);
        case 'min'
            param.stat = cellfun(@min,vals);
        case 'max'
            param.stat = cellfun(@max,vals);
        otherwise
            error('Unknown fiberSummary: %s\n',fiberSummary);

    end

    % In principle, we could write point_stat_array into param.statPoint
    % and the other could be param.statFiber.  Leave it for now, think
    % later.

    % Copy the statistics to the fg pathwayInfo location.  The stat field
    % is for plotting for a whole fiber, and the per point is for the per
    % point case (not always set).
    for ii=1:length(fg.fibers)
        fg.pathwayInfo(ii).pathStat(idx)=param.stat(ii);

        if perPoint
            % We can write algorithms that take the vals, smooth them, use
            % the max, whatever, to create the point coloring.  This is how
            % we will get perFiber coloring.
            switch lower(pointCompute)
                case 'none'
                    fg.pathwayInfo(ii).point_stat_array(idx,:) = vals{ii}(:);
                case 'matchend'
                    tmp = vals{ii}(:);
                    last = min(length(tmp),10);
                    foo = tmp(1:last);
                    foo = foo(foo ~= 0);
                    if ~isempty(foo)
                        tmp(:) = max(foo(:));
                    else
                        tmp(:) = 0;
                    end
                    fg.pathwayInfo(ii).point_stat_array(idx,:) = tmp;
                otherwise
                    error('Unknown pointCompute method %s\n',pointCompute);
            end
        end
    end

end
%This is a hack to edit the range of the colorbar display in Quench. This
%will only be a problem if you clean fibers with the refine selection tool
%in quench
if ~notDefined('mpRange')
    fg.pathwayInfo(1).pathStat(idx)=mpRange(1);
    fg.pathwayInfo(2).pathStat(idx)=mpRange(2);
    param.stat(1)=mpRange(1);
    param.stat(2)=mpRange(2);
end

fg.params{idx} = param;

end

%% Make this better over time
function isNifti = isNiftiStruct(img)

isNifti = true;
if ~isfield(img,'qto_ijk') || ~isfield(img,'data')
    isNifti = false;
end

end


text:
function i4vec_frac_test ( )

%*****************************************************************************80
%
%% I4VEC_FRAC_TEST tests I4VEC_FRAC;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  b = 1;
  c = 2 * n;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_FRAC_TEST\n' );
  fprintf ( 1, '  I4VEC_FRAC: K-th smallest integer vector entry.\n' );
  fprintf ( 1, '  Using initial random number seed = %d\n', seed );

  [ a, seed ] = i4vec_uniform_ab ( n, b, c, seed );

  i4vec_print ( n, a, '  The array to search:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Fractile    Value\n' );
  fprintf ( 1, '\n' );

  for k = 1 : floor ( n / 2 ) : n

    afrac = i4vec_frac ( n, a, k );

    fprintf ( 1, '  %6d  %6d\n', k, afrac );

  end

  return
end


text:
function [Fi,Gi,details] = mpt_project_back_equality(Matrices,Fi,Gi,details,OriginalMatrices);
if isempty(Matrices.getback)
    return
else
    %Original solution given by S1*z+S2*x+S3
    S1 = Matrices.getback.S1;
    S2 = Matrices.getback.S2;
    S3 = Matrices.getback.S3;
    for i=1:length(Fi)
        if isempty(Fi{i})
            Fi{i} = S2;
            Gi{i} = S3;
        else
        Fi{i} = S1*Fi{i} + S2;
        Gi{i} = S1*Gi{i} + S3;
        end
    end
    if nargin>2
        if OriginalMatrices.qp
            for i=1:length(Fi)
                % FIX : Check this...
                details.Ai{i} = 0.5*Fi{i}'*OriginalMatrices.H*Fi{i} + 0.5*(OriginalMatrices.F*Fi{i}+Fi{i}'*OriginalMatrices.F') + OriginalMatrices.Y;
                details.Bi{i} = OriginalMatrices.Cf*Fi{i}+Gi{i}'*OriginalMatrices.F' + Gi{i}'*OriginalMatrices.H*Fi{i} + OriginalMatrices.Cx;
                details.Ci{i} = OriginalMatrices.Cf*Gi{i}+0.5*Gi{i}'*OriginalMatrices.H*Gi{i} + OriginalMatrices.Cc;
            end
        else
            for i=1:length(Fi)
                details.Bi{i} = OriginalMatrices.H*Fi{i} + OriginalMatrices.F;
                details.Ci{i} = OriginalMatrices.H*Gi{i};
            end
        end
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function value = quat_norm ( q )

%*****************************************************************************80
%
%% QUAT_NORM computes the norm of a quaternion.
%
%  Discussion:
%
%    A quaternion is a quadruplet (A,B,C,D) of real numbers, which
%    may be written as
%
%      Q = A + Bi + Cj + Dk.
%
%    The norm of Q is
%
%      norm(Q) = sqrt ( A**2 + B**2 + C**2 + D**2 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real Q(4), the quaternion.
%
%    Output, real VALUE, the norm of the quaternion.
%
  value = sqrt ( sum ( q(1:4).^2 ) );

  return
end