text:
function [p,f] = spm_powell(p,xi,tolsc,func,varargin)
% Powell optimisation method
% FORMAT [p,f] = spm_powell(p,xi,tolsc,func,varargin)
%   p        - Starting parameter values
%   xi       - columns containing directions in which to begin searching
%   tolsc    - stopping criteria, optimisation stops when
%                sqrt(sum(((p-p_prev)./tolsc).^2))<1
%   func     - name of evaluated function
%   varargin - remaining arguments to func (after p)
%
%   p        - final parameter estimates
%   f        - function value at minimum
%__________________________________________________________________________
%
% Method is based on Powell's optimisation method described in
% Numerical Recipes (Press, Flannery, Teukolsky & Vetterling).
%__________________________________________________________________________
% Copyright (C) 2001-2011 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_powell.m 4423 2011-08-04 16:28:51Z guillaume $

p = p(:);
f = feval(func,p,varargin{:});
for iter=1:512
    %if numel(p)>1, fprintf('iteration %d...\n', iter); end;            %-#
    ibig = numel(p);
    pp   = p;
    fp   = f;
    del  = 0;
    for i=1:length(p)
        ft = f;
        [p,junk,f] = min1d(p,xi(:,i),func,f,tolsc,varargin{:});
        if abs(ft-f) > del,
            del  = abs(ft-f);
            ibig = i;
        end
    end
    if numel(p)==1 || sqrt(sum(((p(:)-pp(:))./tolsc(:)).^2))<1, return; end
    ft = feval(func,2.0*p-pp,varargin{:});
    if ft < f
        [p,xi(:,ibig),f] = min1d(p,p-pp,func,f,tolsc,varargin{:});
    end
end
warning('Too many optimisation iterations');

%==========================================================================
% function [p,pi,f] = min1d(p,pi,func,f,tolsc,varargin)
%==========================================================================
function [p,pi,f] = min1d(p,pi,func,f,tolsc,varargin)
% Line search for minimum.

global lnm % used in funeval
lnm      = struct('p',p,'pi',pi,'func',func,'args',[]);
lnm.args = varargin;

min1d_plot('Init', 'Line Minimisation','Function','Parameter Value');
min1d_plot('Set', 0, f);

tol      = 1/sqrt(sum((pi(:)./tolsc(:)).^2));
t        = bracket(f);
[f,pmin] = search(t,tol);
pi       = pi*pmin;
p        = p + pi;

%if length(p)<12,
%    for i=1:length(p), fprintf('%-8.4g ', p(i)); end;                  %-#
%    fprintf('| %.5g\n', f);                                            %-#
%else
%    fprintf('%.5g\n', f);                                              %-#
%end
min1d_plot('Clear');

%==========================================================================
% function f = funeval(p)
%==========================================================================
function f = funeval(p)
% Reconstruct parameters and evaluate.

global lnm % defined in min1d
pt = lnm.p+p.*lnm.pi;
f  = feval(lnm.func,pt,lnm.args{:});
min1d_plot('Set',p,f);

%==========================================================================
% function t = bracket(f)
%==========================================================================
function t = bracket(f)
% Bracket the minimum (t(2)) between t(1) and t(3)

gold   = (1+sqrt(5))/2; % Golden ratio

t(1)   = struct('p',0,'f',f);
t(2).p = 1;
t(2).f = funeval(t(2).p);

% if t(2) not better than t(1) then swap
if t(2).f > t(1).f
    t(3) = t(1);
    t(1) = t(2);
    t(2) = t(3);
end

t(3).p = t(2).p + gold*(t(2).p-t(1).p);
t(3).f = funeval(t(3).p);

while t(2).f > t(3).f

    % fit a polynomial to t
    tmp = cat(1,t.p)-t(2).p;
    pol = pinv([ones(3,1) tmp tmp.^2])*cat(1,t.f);

    % minimum is when gradient of polynomial is zero
    % sign of pol(3) (the 2nd deriv) should be +ve
    if pol(3)>0
        % minimum is when gradient of polynomial is zero
        d    = -pol(2)/(2*pol(3)+eps);

        % A very conservative constraint on the displacement
        if d > (1+gold)*(t(3).p-t(2).p),
            d = (1+gold)*(t(3).p-t(2).p);
        end
        u.p  = t(2).p+d;
    else
        % sign of pol(3) (the 2nd deriv) is not +ve
        % so extend out by golden ratio instead
        u.p  = t(3).p+gold*(t(3).p-t(2).p);
    end

    % FUNCTION EVALUATION
    u.f  = funeval(u.p);

    if (t(2).p < u.p) == (u.p < t(3).p)

        % u is between t(2) and t(3)
        if u.f < t(3).f
            % minimum between t(2) and t(3) - done
            t(1) = t(2);
            t(2) = u;
            return
        elseif u.f > t(2).f
            % minimum between t(1) and u - done
            t(3) = u;
            return;
        end
    end

    % Move all 3 points along
    t(1) = t(2);
    t(2) = t(3);
    t(3) = u;
end

%==========================================================================
% function [f,p] = search(t, tol)
%==========================================================================
function [f,p] = search(t, tol)
% Brent's method for line searching - given that minimum is bracketed

gold1 = 1-(sqrt(5)-1)/2;

% Current and previous displacements
d     = Inf;
pd    = Inf;

% sort t into best first order
[junk,ind] = sort(cat(1,t.f));
t   = t(ind);
brk = [min(cat(1,t.p)) max(cat(1,t.p))];

for iter=1:128
    % check stopping criterion
    if abs(t(1).p - 0.5*(brk(1)+brk(2)))+0.5*(brk(2)-brk(1)) <= 2*tol
        p = t(1).p;
        f = t(1).f;
        return;
    end

    % keep last two displacents
    ppd = pd;
    pd  = d;

    % fit a polynomial to t
    tmp = cat(1,t.p)-t(1).p;
    pol = pinv([ones(3,1) tmp tmp.^2])*cat(1,t.f);

    % minimum is when gradient of polynomial is zero
    d   = -pol(2)/(2*pol(3)+eps);
    u.p = t(1).p+d;

    % check so that displacement is less than the last but two,
    % that the displaced point is between the brackets
    % and that the solution is a minimum rather than a maximum
    eps2 = 2*eps*abs(t(1).p)+eps;
    if abs(d) > abs(ppd)/2 || u.p < brk(1)+eps2 || u.p > brk(2)-eps2 || pol(3)<=0
        % if criteria are not met, then golden search into the larger part
        if t(1).p >= 0.5*(brk(1)+brk(2)),
            d = gold1*(brk(1)-t(1).p);
        else
            d = gold1*(brk(2)-t(1).p);
        end
        u.p = t(1).p+d;
    end

    % FUNCTION EVALUATION
    u.f = funeval(u.p);

    % Insert the new point into the appropriate position and update
    % the brackets if necessary
    if u.f <= t(1).f
        if u.p >= t(1).p, brk(1)=t(1).p; else brk(2)=t(1).p; end
        t(3) = t(2);
        t(2) = t(1);
        t(1) = u;
    else
        if u.p < t(1).p, brk(1)=u.p; else brk(2)=u.p; end
        if u.f <= t(2).f
            t(3) = t(2);
            t(2) = u;
        elseif u.f <= t(3).f
            t(3) = u;
        end
    end
end

%==========================================================================
% function min1d_plot(action,arg1,arg2,arg3)
%==========================================================================
function min1d_plot(action,arg1,arg2,arg3)
% Visual output for line minimisation
persistent min1dplot

if ~nargin, action = 'Init'; end

% Find the Interactive window and exit if not
%--------------------------------------------------------------------------
fg = spm_figure('FindWin','Interactive');
if isempty(fg), return; end

%-Initialize
%--------------------------------------------------------------------------
if strcmpi(action,'init')
    if nargin<4, arg3 = 'Function';          end
    if nargin<3, arg2 = 'Value';             end
    if nargin<2, arg1 = 'Line minimisation'; end

    min1dplot = struct('pointer',get(fg,'Pointer'),...
                       'name',   get(fg,'Name'),...
                       'ax',     [],...
                       'buffer', get(fg,'DoubleBuffer'));
    min1d_plot('Clear');
    set(fg,'Pointer','Watch');
    set(fg,'DoubleBuffer','on');
    min1dplot.ax = axes('Position', [0.15 0.1 0.8 0.75],...
                        'Box',      'on',...
                        'Parent',   fg);
    lab = get(min1dplot.ax,'Xlabel');
    set(lab,'string',arg3,'FontSize',10);
    lab = get(min1dplot.ax,'Ylabel');
    set(lab,'string',arg2,'FontSize',10);
    lab = get(min1dplot.ax,'Title');
    set(lab,'string',arg1);
    line('Xdata',[], 'Ydata',[],...
        'LineWidth',2,'Tag','LinMinPlot',...
        'LineStyle','-','Marker','o',...
        'Parent',min1dplot.ax);
    drawnow;

%-Reset
%--------------------------------------------------------------------------
elseif strcmpi(action,'set')
    br = findobj(fg,'Tag','LinMinPlot');
    if ~isempty(br)
        [xd,indx] = sort([get(br,'Xdata') arg1]);
        yd = [get(br,'Ydata') arg2];
        yd = yd(indx);
        set(br,'Ydata',yd,'Xdata',xd);
        drawnow;
    end

%-Clear
%--------------------------------------------------------------------------
elseif strcmpi(action,'clear')
    fg = spm_figure('FindWin','Interactive');
    if isstruct(min1dplot)
        if ishandle(min1dplot.ax), delete(min1dplot.ax); end
        set(fg,'Pointer',min1dplot.pointer);
        set(fg,'Name',min1dplot.name);
        set(fg,'DoubleBuffer',min1dplot.buffer);
    end
    spm_figure('Clear',fg);
    drawnow;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function DOMnode = xml_write(filename, tree, RootName, Pref)
%XML_WRITE  Writes Matlab data structures to XML file
%
% DESCRIPTION
% xml_write( filename, tree) Converts Matlab data structure 'tree' containing
% cells, structs, numbers and strings to Document Object Model (DOM) node
% tree, then saves it to XML file 'filename' using Matlab's xmlwrite
% function. Optionally one can also use alternative version of xmlwrite
% function which directly calls JAVA functions for XML writing without
% MATLAB middleware. This function is provided as a patch to existing
% bugs in xmlwrite (in R2006b).
%
% xml_write(filename, tree, RootName, Pref) allows you to specify
% additional preferences about file format
%
% DOMnode = xml_write([], tree) same as above except that DOM node is
% not saved to the file but returned.
%
% INPUT
%   filename     file name
%   tree         Matlab structure tree to store in xml file.
%   RootName     String with XML tag name used for root (top level) node
%                Optionally it can be a string cell array storing: Name of
%                root node, document "Processing Instructions" data and
%                document "comment" string
%   Pref         Other preferences:
%     Pref.ItemName - default 'item' -  name of a special tag used to
%                     itemize cell or struct arrays
%     Pref.XmlEngine - let you choose the XML engine. Currently default is
%       'Xerces', which is using directly the apache xerces java file.
%       Other option is 'Matlab' which uses MATLAB's xmlwrite and its
%       XMLUtils java file. Both options create identical results except in
%       case of CDATA sections where xmlwrite fails.
%     Pref.CellItem - default 'true' - allow cell arrays to use 'item'
%       notation. See below.
%    Pref.RootOnly - default true - output variable 'tree' corresponds to
%       xml file root element, otherwise it correspond to the whole file.
%     Pref.StructItem - default 'true' - allow arrays of structs to use
%       'item' notation. For example "Pref.StructItem = true" gives:
%
%
%             <item> ... <\item>
%             <item> ... <\item>
%           <\b>
%         <\a>
%       while "Pref.StructItem = false" gives:
%
%            ... <\b>
%            ... <\b>
%         <\a>
%
%
% Several special xml node types can be created if special tags are used
% for field names of 'tree' nodes:
%  - node.CONTENT - stores data section of the node if other fields
%    (usually ATTRIBUTE are present. Usually data section is stored
%    directly in 'node'.
%  - node.ATTRIBUTE.name - stores node's attribute called 'name'.
%  - node.COMMENT - create comment child node from the string. For global
%    comments see "RootName" input variable.
%  - node.PROCESSING_INSTRUCTIONS - create "processing instruction" child
%    node from the string. For global "processing instructions" see
%    "RootName" input variable.
%  - node.CDATA_SECTION - stores node's CDATA section (string). Only works
%    if Pref.XmlEngine='Xerces'. For more info, see comments of F_xmlwrite.
%  - other special node types like: document fragment nodes, document type
%    nodes, entity nodes and notation nodes are not being handled by
%    'xml_write' at the moment.
%
% OUTPUT
%   DOMnode      Document Object Model (DOM) node tree in the format
%                required as input to xmlwrite. (optional)
%
% EXAMPLES:
%   MyTree=[];
%   MyTree.MyNumber = 13;
%   MyTree.MyString = 'Hello World';
%   xml_write('test.xml', MyTree);
%   type('test.xml')
%   %See also xml_tutorial.m
%
% See also
%   xml_read, xmlread, xmlwrite
%
% Written by Jarek Tuszynski, SAIC, jaroslaw.w.tuszynski_at_saic.com

%% Check Matlab Version
v = ver('MATLAB');
v = str2double(regexp(v.Version, '\d.\d','match','once'));
if (v<7)
  error('Your MATLAB version is too old. You need version 7.0 or newer.');
end

%% default preferences
DPref.TableName  = {'tr','td'}; % name of a special tags used to itemize 2D cell arrays
DPref.ItemName   = 'item'; % name of a special tag used to itemize 1D cell arrays
DPref.StructItem = true;  % allow arrays of structs to use 'item' notation
DPref.CellItem   = true;  % allow cell arrays to use 'item' notation
DPref.StructTable= 'Html';
DPref.CellTable  = 'Html';
DPref.XmlEngine  = 'Matlab';  % use matlab provided XMLUtils
%DPref.XmlEngine  = 'Xerces';  % use Xerces xml generator directly
DPref.PreserveSpace = false; % Preserve or delete spaces at the beggining and the end of stings?
RootOnly         = true;  % Input is root node only
GlobalProcInst = [];
GlobalComment  = [];
GlobalDocType  = [];

%% read user preferences
if (nargin>3)
  if (isfield(Pref, 'TableName' )),  DPref.TableName  = Pref.TableName; end
  if (isfield(Pref, 'ItemName'  )), DPref.ItemName   = Pref.ItemName;   end
  if (isfield(Pref, 'StructItem')), DPref.StructItem = Pref.StructItem; end
  if (isfield(Pref, 'CellItem'  )), DPref.CellItem   = Pref.CellItem;   end
  if (isfield(Pref, 'CellTable')),   DPref.CellTable  = Pref.CellTable; end
  if (isfield(Pref, 'StructTable')), DPref.StructTable= Pref.StructTable; end
  if (isfield(Pref, 'XmlEngine' )), DPref.XmlEngine  = Pref.XmlEngine;  end
  if (isfield(Pref, 'RootOnly'  )), RootOnly         = Pref.RootOnly;   end
  if (isfield(Pref, 'PreserveSpace')), DPref.PreserveSpace = Pref.PreserveSpace; end
end
if (nargin<3 || isempty(RootName)), RootName=inputname(2); end
if (isempty(RootName)), RootName='ROOT'; end
if (iscell(RootName)) % RootName also stores global text node data
  rName = RootName;
  RootName = char(rName{1});
  if (length(rName)>1), GlobalProcInst = char(rName{2}); end
  if (length(rName)>2), GlobalComment  = char(rName{3}); end
  if (length(rName)>3), GlobalDocType  = char(rName{4}); end
end
if(~RootOnly && isstruct(tree))  % if struct than deal with each field separatly
  fields = fieldnames(tree);
  for i=1:length(fields)
    field = fields{i};
    x = tree(1).(field);
    if (strcmp(field, 'COMMENT'))
      GlobalComment = x;
    elseif (strcmp(field, 'PROCESSING_INSTRUCTION'))
      GlobalProcInst = x;
    elseif (strcmp(field, 'DOCUMENT_TYPE'))
      GlobalDocType = x;
    else
      RootName = field;
      t = x;
    end
  end
  tree = t;
end

%% Initialize jave object that will store xml data structure
RootName = varName2str(RootName);
if (~isempty(GlobalDocType))
  %   n = strfind(GlobalDocType, ' ');
  %   if (~isempty(n))
  %     dtype = com.mathworks.xml.XMLUtils.createDocumentType(GlobalDocType);
  %   end
  %   DOMnode = com.mathworks.xml.XMLUtils.createDocument(RootName, dtype);
  warning('xml_io_tools:write:docType', ...
    'DOCUMENT_TYPE node was encountered which is not supported yet. Ignoring.');
end
DOMnode = com.mathworks.xml.XMLUtils.createDocument(RootName);

%% Use recursive function to convert matlab data structure to XML
root = DOMnode.getDocumentElement;
struct2DOMnode(DOMnode, root, tree, DPref.ItemName, DPref);

%% Remove the only child of the root node
root   = DOMnode.getDocumentElement;
Child  = root.getChildNodes; % create array of children nodes
nChild = Child.getLength;    % number of children
if (nChild==1)
  node = root.removeChild(root.getFirstChild);
  while(node.hasChildNodes)
    root.appendChild(node.removeChild(node.getFirstChild));
  end
  while(node.hasAttributes)            % copy all attributes
    root.setAttributeNode(node.removeAttributeNode(node.getAttributes.item(0)));
  end
end

%% Save exotic Global nodes
if (~isempty(GlobalComment))
  DOMnode.insertBefore(DOMnode.createComment(GlobalComment), DOMnode.getFirstChild());
end
if (~isempty(GlobalProcInst))
  n = strfind(GlobalProcInst, ' ');
  if (~isempty(n))
    proc = DOMnode.createProcessingInstruction(GlobalProcInst(1:(n(1)-1)),...
      GlobalProcInst((n(1)+1):end));
    DOMnode.insertBefore(proc, DOMnode.getFirstChild());
  end
end
% Not supported yet as the code below does not work
% if (~isempty(GlobalDocType))
%   n = strfind(GlobalDocType, ' ');
%   if (~isempty(n))
%     dtype = DOMnode.createDocumentType(GlobalDocType);
%     DOMnode.insertBefore(dtype, DOMnode.getFirstChild());
%   end
% end

%% save java DOM tree to XML file
if (~isempty(filename))
  if (strcmpi(DPref.XmlEngine, 'Xerces'))
    xmlwrite_xerces(filename, DOMnode);
  else
    xmlwrite(filename, DOMnode);
  end
end

%% =======================================================================
%  === struct2DOMnode Function ===========================================
%  =======================================================================
function [] = struct2DOMnode(xml, parent, s, TagName, Pref)
% struct2DOMnode is a recursive function that converts matlab's structs to
% DOM nodes.
% INPUTS:
%  xml - jave object that will store xml data structure
%  parent - parent DOM Element
%  s - Matlab data structure to save
%  TagName - name to be used in xml tags describing 's'
%  Pref - preferenced
% OUTPUT:
%  parent - modified 'parent'

% perform some conversions
if (ischar(s) && min(size(s))>1) % if 2D array of characters
  s=cellstr(s);                  % than convert to cell array
end
% if (strcmp(TagName, 'CONTENT'))
%   while (iscell(s) && length(s)==1), s = s{1}; end % unwrap cell arrays of length 1
% end
TagName  = varName2str(TagName);

%% == node is a 2D cell array ==
% convert to some other format prior to further processing
nDim = nnz(size(s)>1);  % is it a scalar, vector, 2D array, 3D cube, etc?
if (iscell(s) && nDim==2 && strcmpi(Pref.CellTable, 'Matlab'))
  s = var2str(s, Pref.PreserveSpace);
end
if (nDim==2 && (iscell  (s) && strcmpi(Pref.CellTable,   'Vector')) || ...
               (isstruct(s) && strcmpi(Pref.StructTable, 'Vector')))
  s = s(:);
end
if (nDim>2), s = s(:); end % can not handle this case well
nItem = numel(s);
nDim  = nnz(size(s)>1);  % is it a scalar, vector, 2D array, 3D cube, etc?

%% == node is a cell ==
if (iscell(s)) % if this is a cell or cell array
  if ((nDim==2 && strcmpi(Pref.CellTable,'Html')) || (nDim< 2 && Pref.CellItem))
    % if 2D array of cells than can use HTML-like notation or if 1D array
    % than can use item notation
    if (strcmp(TagName, 'CONTENT')) % CONTENT nodes already have <TagName> ... </TagName>
      array2DOMnode(xml, parent, s, Pref.ItemName, Pref ); % recursive call
    else
      node = xml.createElement(TagName);   % <TagName> ... </TagName>
      array2DOMnode(xml, node, s, Pref.ItemName, Pref ); % recursive call
      parent.appendChild(node);
    end
  else % use  <TagName>...<\TagName> <TagName>...<\TagName> notation
    array2DOMnode(xml, parent, s, TagName, Pref ); % recursive call
  end
%% == node is a struct ==
elseif (isstruct(s))  % if struct than deal with each field separatly
  if ((nDim==2 && strcmpi(Pref.StructTable,'Html')) || (nItem>1 && Pref.StructItem))
    % if 2D array of structs than can use HTML-like notation or
    % if 1D array of structs than can use 'items' notation
    node = xml.createElement(TagName);
    array2DOMnode(xml, node, s, Pref.ItemName, Pref ); % recursive call
    parent.appendChild(node);
  elseif (nItem>1) % use  <TagName>...<\TagName> <TagName>...<\TagName> notation
    array2DOMnode(xml, parent, s, TagName, Pref ); % recursive call
  else % otherwise save each struct separatelly
    fields = fieldnames(s);
    node = xml.createElement(TagName);
    for i=1:length(fields) % add field by field to the node
      field = fields{i};
      x = s.(field);
      switch field
        case {'COMMENT', 'CDATA_SECTION', 'PROCESSING_INSTRUCTION'}
          if iscellstr(x)  % cell array of strings -> add them one by one
            array2DOMnode(xml, node, x(:), field, Pref ); % recursive call will modify 'node'
          elseif ischar(x) % single string -> add it
            struct2DOMnode(xml, node, x, field, Pref ); % recursive call will modify 'node'
          else % not a string - Ignore
            warning('xml_io_tools:write:badSpecialNode', ...
             ['Struct field named ',field,' encountered which was not a string. Ignoring.']);
          end
        case 'ATTRIBUTE' % set attributes of the node
          if (isempty(x)), continue; end
          if (isstruct(x))
            attName = fieldnames(x);       % get names of all the attributes
            for k=1:length(attName)        % attach them to the node
              att = xml.createAttribute(varName2str(attName(k)));
              att.setValue(var2str(x.(attName{k}),Pref.PreserveSpace));
              node.setAttributeNode(att);
            end
          else
            warning('xml_io_tools:write:badAttribute', ...
              'Struct field named ATTRIBUTE encountered which was not a struct. Ignoring.');
          end
        otherwise                            % set children of the node
          struct2DOMnode(xml, node, x, field, Pref ); % recursive call will modify 'node'
      end
    end  % end for i=1:nFields
    parent.appendChild(node);
  end
%% == node is a leaf node ==
else  % if not a struct and not a cell than it is a leaf node
  switch TagName % different processing depending on desired type of the node
    case 'COMMENT'   % create comment node
      com = xml.createComment(s);
      parent.appendChild(com);
    case 'CDATA_SECTION' % create CDATA Section
      cdt = xml.createCDATASection(s);
      parent.appendChild(cdt);
    case 'PROCESSING_INSTRUCTION' % set attributes of the node
      OK = false;
      if (ischar(s))
        n = strfind(s, ' ');
        if (~isempty(n))
          proc = xml.createProcessingInstruction(s(1:(n(1)-1)),s((n(1)+1):end));
          parent.insertBefore(proc, parent.getFirstChild());
          OK = true;
        end
      end
      if (~OK)
        warning('xml_io_tools:write:badProcInst', ...
          ['Struct field named PROCESSING_INSTRUCTION need to be',...
          ' a string, for example: xml-stylesheet type="text/css" ', ...
          'href="myStyleSheet.css". Ignoring.']);
      end
    case 'CONTENT' % this is text part of already existing node
      txt  = xml.createTextNode(var2str(s, Pref.PreserveSpace)); % convert to text
      parent.appendChild(txt);
    otherwise      % I guess it is a regular text leaf node
      txt  = xml.createTextNode(var2str(s, Pref.PreserveSpace));
      node = xml.createElement(TagName);
      node.appendChild(txt);
      parent.appendChild(node);
  end
end % of struct2DOMnode function

%% =======================================================================
%  === array2DOMnode Function ============================================
%  =======================================================================
function [] = array2DOMnode(xml, parent, s, TagName, Pref)
% Deal with 1D and 2D arrays of cell or struct. Will modify 'parent'.
nDim = nnz(size(s)>1);  % is it a scalar, vector, 2D array, 3D cube, etc?
switch nDim
  case 2 % 2D array
    for r=1:size(s,1)
      subnode = xml.createElement(Pref.TableName{1});
      for c=1:size(s,2)
        v = s(r,c);
        if iscell(v), v = v{1}; end
        struct2DOMnode(xml, subnode, v, Pref.TableName{2}, Pref ); % recursive call
      end
      parent.appendChild(subnode);
    end
  case 1 %1D array
    for iItem=1:numel(s)
      v = s(iItem);
      if iscell(v), v = v{1}; end
      struct2DOMnode(xml, parent, v, TagName, Pref ); % recursive call
    end
  case 0 % scalar -> this case should never be called
    if ~isempty(s)
      if iscell(s), s = s{1}; end
      struct2DOMnode(xml, parent, s, TagName, Pref );
    end
end

%% =======================================================================
%  === var2str Function ==================================================
%  =======================================================================
function str = var2str(object, PreserveSpace)
% convert matlab variables to a string
switch (1)
  case isempty(object)
    str = '';
  case (isnumeric(object) || islogical(object))
    if ndims(object)>2, object=object(:); end  % can't handle arrays with dimention > 2
    str=mat2str(object);           % convert matrix to a string
    % mark logical scalars with [] (logical arrays already have them) so the xml_read
    % recognizes them as MATLAB objects instead of strings. Same with sparse
    % matrices
    if ((islogical(object) && isscalar(object)) || issparse(object)),
      str = ['[' str ']'];
    end
    if (isinteger(object)),
      str = ['[', class(object), '(', str ')]'];
    end
  case iscell(object)
    if ndims(object)>2, object=object(:); end  % can't handle cell arrays with dimention > 2
    [nr nc] = size(object);
    obj2 = object;
    for i=1:length(object(:))
      str = var2str(object{i}, PreserveSpace);
      if (ischar(object{i})), object{i} = ['''' object{i} '''']; else object{i}=str; end
      obj2{i} = [object{i} ','];
    end
    for r = 1:nr, obj2{r,nc} = [object{r,nc} ';']; end
    obj2 = obj2.';
    str = ['{' obj2{:} '}'];
  case isstruct(object)
    str='';
    warning('xml_io_tools:write:var2str', ...
      'Struct was encountered where string was expected. Ignoring.');
  case isa(object, 'function_handle')
    str = ['[@' char(object) ']'];
  case ischar(object)
    str = object;
  otherwise
    str = char(object);
end

%% string clean-up
str=str(:); str=str.';            % make sure this is a row vector of char's
if (~isempty(str))
  str(str<32|str==127)=' ';       % convert no-printable characters to spaces
  if (~PreserveSpace)
    str = strtrim(str);             % remove spaces from begining and the end
    str = regexprep(str,'\s+',' '); % remove multiple spaces
  end
end

%% =======================================================================
%  === var2Namestr Function ==============================================
%  =======================================================================
function str = varName2str(str)
% convert matlab variable names to a sting
str = char(str);
p   = strfind(str,'0x');
if (~isempty(p))
  for i=1:length(p)
    before = str( p(i)+(0:3) );          % string to replace
    after  = char(hex2dec(before(3:4))); % string to replace with
    str = regexprep(str,before,after, 'once', 'ignorecase');
    p=p-3; % since 4 characters were replaced with one - compensate
  end
end
str = regexprep(str,'_COLON_',':', 'once', 'ignorecase');
str = regexprep(str,'_DASH_' ,'-', 'once', 'ignorecase');

author: TadasBaltrusaitis

repo: OpenFace


text:
function [mfile,idx]=get_mfile(obj, fn)
% get mfile either by name or by index
%
% [mfile,idx]=get_mfile(obj, fn)
%
% Inputs:
%   obj                 @MOcovMFileCollection instance
%   fn                  filename or integer
%
% Outputs:
%   mfile               @MOcovMFile that either has name fn (if fn is a
%                       string) or is the fn-th file (if fn is numeric)
%   idx                 index of mfile within obj
%

    if isnumeric(fn)
        lookup_func=@get_mfile_numeric;
    else
        lookup_func=@get_mfile_by_name;
    end

    [mfile,idx]=lookup_func(obj,fn);

function [mfile,idx]=get_mfile_numeric(obj,idx)
    mfile=obj.mfiles{idx};

function [mfile,idx]=get_mfile_by_name(obj,fn)
    root_dir=obj.root_dir;
    abs_fn=mocov_get_absolute_path(fullfile(root_dir,fn));

    mfiles=obj.mfiles;
    n=numel(mfiles);
    for k=1:n
        mfile=mfiles{k};
        mfile_fn=get_filename(mfile);
        if strcmp(mfile_fn, abs_fn)
            [mfile,idx]=get_mfile_numeric(obj,k);
            return;
        end
    end

    error('Not found: %s', fn);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function obj = updateSlice(obj, dataIndex)

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

    %-INITIALIZATIONS-%

    axIndex = obj.getAxisIndex(obj.State.Plot(dataIndex).AssociatedAxis);
    plotData = get(obj.State.Plot(dataIndex).Handle);
    axisData = get(plotData.Parent);
    figureData = get(obj.State.Figure.Handle);
    [xSource, ~] = findSourceAxis(obj,axIndex);

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

    %-get trace data-%
    xData = plotData.XData;
    yData = plotData.YData;
    zData = plotData.ZData;
    cData = plotData.CData;

    xDataSurf = zeros(2*(size(xData)-1));
    yDataSurf = zeros(2*(size(xData)-1));
    zDataSurf = zeros(2*(size(xData)-1));
    cDataSurf = zeros(2*(size(xData)-1));

    for n = 1:size(xData,2)-1
        n2 = 2*(n-1) + 1;

        for m = 1:size(xData,1)-1
            m2 = 2*(m-1) + 1;

            xDataSurf(m2:m2+1,n2:n2+1) = xData(m:m+1,n:n+1);
            yDataSurf(m2:m2+1,n2:n2+1) = yData(m:m+1,n:n+1);
            zDataSurf(m2:m2+1,n2:n2+1) = zData(m:m+1,n:n+1);

            if strcmp(plotData.FaceColor, 'flat')
                cDataSurf(m2:m2+1,n2:n2+1) = ones(2,2)*cData(m,n);
            elseif strcmp(plotData.FaceColor, 'interp')
                cDataSurf(m2:m2+1,n2:n2+1) = cData(m:m+1,n:n+1);
            end

        end
    end

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

    %-set trace-%
    obj.data{dataIndex}.type = 'surface';
    obj.data{dataIndex}.name = plotData.DisplayName;
    obj.data{dataIndex}.visible = strcmp(plotData.Visible,'on');
    obj.data{dataIndex}.scene = sprintf('scene%d', xSource);
    obj.data{dataIndex}.showscale = false;
    obj.data{dataIndex}.surfacecolor = cDataSurf;

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

    %-set trace data-%
    obj.data{dataIndex}.x = xDataSurf;
    obj.data{dataIndex}.y = yDataSurf;
    obj.data{dataIndex}.z = zDataSurf;

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

    %-update face color-%
    updateSurfaceFaceColor(obj, dataIndex, cDataSurf);

    %-update edge color-%
    if isnumeric(plotData.EdgeColor)
        updateSurfaceEdgeColor(obj, dataIndex);
    end

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

function updateScene(obj, dataIndex)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

function updateSurfaceEdgeColor(obj, dataIndex)

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

    %-INITIALIZATIONS-%

    axIndex = obj.getAxisIndex(obj.State.Plot(dataIndex).AssociatedAxis);
    plotData = get(obj.State.Plot(dataIndex).Handle);
    axisData = get(plotData.Parent);

    xData = plotData.XData;
    yData = plotData.YData;
    zData = plotData.ZData;
    cData = plotData.CData;
    edgeColor = plotData.EdgeColor;
    cData = plotData.CData;
    cLim = axisData.CLim;
    colorMap = axisData.Colormap;

    xConst = ( xData(:) - min(xData(:)) ) <= 1e-6;
    yConst = ( yData(:) - min(yData(:)) ) <= 1e-6;

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

    %-edge lines in x direction-%
    xContourSize = mean(diff(xData(1,:)));
    xContourStard = min(xData(1,:));
    xContourEnd = max(xData(1,:));

    obj.data{dataIndex}.contours.x.show = true;
    obj.data{dataIndex}.contours.x.start = xContourStard;
    obj.data{dataIndex}.contours.x.end = xContourEnd;
    obj.data{dataIndex}.contours.x.size = xContourSize;

    %-edge lines in y direction-%
    yContourSize = mean(diff(yData(:,1)));
    yContourStard = min(yData(:,1));
    yContourEnd = max(yData(:,1));

    obj.data{dataIndex}.contours.y.show = true;
    obj.data{dataIndex}.contours.y.start = yContourStard;
    obj.data{dataIndex}.contours.y.end = yContourEnd;
    obj.data{dataIndex}.contours.y.size = yContourSize;

    %-edge lines in z direction-%

    if all(xConst) || all(yConst)
        zContourSize = mean(diff(zData(1,:)));
        zContourStard = min(zData(1,:));
        zContourEnd = max(zData(1,:));

        obj.data{dataIndex}.contours.z.show = true;
        obj.data{dataIndex}.contours.z.start = zContourStard;
        obj.data{dataIndex}.contours.z.end = zContourEnd;
        obj.data{dataIndex}.contours.z.size = zContourSize;
    end

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

    %-coloring-%
    numColor = 255 * edgeColor;
    stringColor = getStringColor(numColor);

    obj.data{dataIndex}.contours.x.color = stringColor;
    obj.data{dataIndex}.contours.y.color = stringColor;
    obj.data{dataIndex}.contours.z.color = stringColor;

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

function updateSurfaceFaceColor(obj, dataIndex, surfaceColor)

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

    %-INITIALIZATIONS-%

    axIndex = obj.getAxisIndex(obj.State.Plot(dataIndex).AssociatedAxis);
    plotData = get(obj.State.Plot(dataIndex).Handle);
    axisData = get(plotData.Parent);

    faceColor = plotData.FaceColor;
    cLim = axisData.CLim;
    colorMap = axisData.Colormap;

    obj.data{dataIndex}.cauto = false;
    obj.data{dataIndex}.autocolorscale = false;

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

    if isnumeric(faceColor)
        numColor = 255 * faceColor;
        stringColor = getStringColor(numColor);

        colorScale{1} = {0, stringColor};
        colorScale{2} = {1, stringColor};
        obj.data{dataIndex}.colorscale = colorScale;

    elseif ismember(faceColor, {'flat', 'interp'})

        nColors = size(colorMap, 1);

        for c = 1:nColors
            stringColor = getStringColor(255*colorMap(c,:));
            colorScale{c} = {(c-1)/(nColors-1), stringColor};
        end

        obj.data{dataIndex}.cmin = cLim(1);
        obj.data{dataIndex}.cmax = cLim(2);

    end

    obj.data{dataIndex}.surfacecolor = surfaceColor;
    obj.data{dataIndex}.colorscale = colorScale;

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

function stringColor = getStringColor(numColor)

    stringColor = sprintf('rgb(%f,%f,%f)', numColor);
end

author: plotly


text:
function [pval m] = circ_otest(alpha, sz, w)
%
% [pval, m] = circ_otest(alpha,sz,w)
%   Computes Omnibus or Hodges-Ajne test for non-uniformity of circular data.
%   H0: the population is uniformly distributed around the circle
%   HA: the population is not distributed uniformly around the circle
%
%   Alternative to the Rayleigh and Rao's test. Works well for unimodal,
%   bimodal or multimodal data. If requirements of the Rayleigh test are
%   met, the latter is more powerful.
%
%   Input:
%     alpha	sample of angles in radians
%     [sz   step size for evaluating distribution, default 1 degree
%     [w		number of incidences in case of binned angle data]

%   Output:
%     pval  p-value
%     m     minimum number of samples falling in one half of the circle
%
% PHB 3/16/2009
%
% References:
%   Biostatistical Analysis, J. H. Zar
%   A bivariate sign test, J. L. Hodges et al., 1955
%   A simple test for uniformity of a circular distribution, B. Ajne, 1968
%
% Circular Statistics Toolbox for Matlab

% By Philipp Berens, 2009
% berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html

if size(alpha,2) > size(alpha,1)
	alpha = alpha';
end

if nargin < 2 || isempty(sz)
  sz = circ_ang2rad(1);
end

if nargin < 3
  w = ones(size(alpha));
else
  if length(alpha)~=length(w)
    error('Input length does not match.')
  end
  w =w(:);
end

alpha = mod(alpha,2*pi);
n = sum(w);
dg = 0:sz:pi;

m1 = zeros(size(dg));
m2 = zeros(size(dg));
for i=1:length(dg)
  m1(i) = sum((alpha > dg(i) & alpha < pi + dg(i)).*w);
  m2(i) = n - m1(i);
end
m = min(min([m1;m2]));

if n > 50
  % approximation by Ajne (1968)
  A = pi*sqrt(n) / 2 / (n-2*m);
  pval = sqrt(2*pi) / A * exp(-pi^2/8/A^2);
else
  % exact formula by Hodges (1955)
  pval = 2^(1-n) * (n-2*m) * nchoosek(n,m);
end

author: brainstorm-tools


text:
%  Figure 7.34      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
%% fig7_34.m
%% SRL for estimator design
clf;
num=[1];
den=conv([1 0 1],[1 0 1]);
rlocus(num,den);
grid;
title('Fig. 7.34: SRL for Estimator Design');
text(0,0.9,'q --> 0');
text(0,-0.9,'q --> 0');
text(-5,4,'q-->\infty');
text(-5,-4,'q-->\infty');

author: Sable

repo: mcbench-benchmarks


text:
function varargout = process_decoding_maxcorr( varargin )
% PROCESS_DECODING_MAXCORR: Decoding of multiple conditions using max-correlation

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Martin Cousineau, 2019; Dimitrios Pantazis, 2019

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Max-correlation decoding';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Decoding';
    sProcess.Index       = 702;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Decoding';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data'};
    sProcess.OutputTypes = {'matrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 2;

    % Definition of the options
    sProcess.options.description.Comment = ['Apply max-correlation classifier on MEG trials.
' ...
                                            'Uses trial subaverages and permutations.

'];
    sProcess.options.description.Type    = 'label';
    % === Sensor types
    sProcess.options.sensortypes.Comment = 'Sensor types or names (empty=all): ';
    sProcess.options.sensortypes.Type    = 'text';
    sProcess.options.sensortypes.Value   = 'MEG';
    % === lowpass filtering
    sProcess.options.lowpass.Comment = 'Low-pass cutoff frequency (0=disabled): ';
    sProcess.options.lowpass.Type    = 'value';
    sProcess.options.lowpass.Value   = {30,'Hz',2};
    % === number of permutations
    sProcess.options.num_permutations.Comment = 'Number of permutations: ';
    sProcess.options.num_permutations.Type    = 'value';
    sProcess.options.num_permutations.Value   = {50,'',0};
    % === trial bin size for sub-averaging
    sProcess.options.kfold.Comment = 'Number of folds: ';
    sProcess.options.kfold.Type    = 'value';
    sProcess.options.kfold.Value   = {5,'',0};
    % === decoding method
    sProcess.options.method.Comment = {'Pairwise', 'Temporal generalization', 'Multiclass', 'Decoding method:'};
    sProcess.options.method.Type    = 'radio_line';
    sProcess.options.method.Value   = 1;
    % === decoding model
    sProcess.options.model.Comment = 'Decoding model: ';
    sProcess.options.model.Type    = 'text';
    sProcess.options.model.Value   = 'maxcorr';
    sProcess.options.model.Hidden  = 1;
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    % Call mother process with actual implementation
    OutputFiles = process_decoding_svm('Run', sProcess, sInputs);
end


text:
function [x,fval,exitflag,info,Opt] = opti_lsqcurvefit(fun,x0,xdata,ydata,lb,ub,opts)
%OPTI_LSQCURVEFIT Solve a NLS using an OPTI NLS Solver (Matlab Overload)
%
%   [x,fval,exitflag,info] = opti_lsqcurvefit(fun,x0,xdata,ydata,lb,ub)
%   solves the nonlinear least squares problem sum[(f(x,xdata)-ydata)^2]
%   where fun is the nonlinear function to be fitted [fun(x,xdata)],
%   starting at x0, to the sample data ydata. Optional bounds lb and ub can
%   be placed on the decision variables x.
%
%   [x,fval,exitflag,info] = opti_lsqcurvefit(fun,...,ub,opts) allows the
%   user to specify optiset options. This includes specifying a solver via
%   the 'solver' field of optiset.
%
%   [x,...,info,Opt] = opti_lsqcurvefit(fun,...) returns the internally
%   built OPTI object.

%   Copyright (C) 2011 Jonathan Currie (I2C2)

% Handle missing arguments
if nargin < 7, opts = optiset; end
if nargin < 6, ub = []; end
if nargin < 5, lb = []; end
if nargin < 4, error('You must supply at least 4 arguments to opti_lsqcurvefit'); end

%Build OPTI Object
Opt = opti('fun',fun,'data',xdata,ydata,'bounds',lb,ub,'x0',x0,'options',opts);

%Solve
[x,fval,exitflag,info] = solve(Opt);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [Q_k,T_k,r,anorm,ierr,work] = lanpro(A,nin,kmax,r,options,...
    Q_k,T_k,anorm)

%LANPRO   Lanczos tridiagonalization with partial reorthogonalization
%   LANPRO computes the Lanczos tridiagonalization of a real symmetric
%   matrix using the symmetric Lanczos algorithm with partial
%   reorthogonalization.
%
%   [Q_K,T_K,R,ANORM,IERR,WORK] = LANPRO(A,K,R0,OPTIONS,Q_old,T_old)
%   [Q_K,T_K,R,ANORM,IERR,WORK] = LANPRO('Afun',N,K,R0,OPTIONS,Q_old,T_old)
%
%   Computes K steps of the Lanczos algorithm with starting vector R0,
%   and returns the K x K tridiagonal T_K, the N x K matrix Q_K
%   with semiorthonormal columns and the residual vector R such that
%
%        A*Q_K = Q_K*T_K + R .
%
%   Partial reorthogonalization is used to keep the columns of Q_K
%   semiorthogonal:
%        MAX(DIAG((eye(k) - Q_K'*Q_K))) <= OPTIONS.delta.
%
%
%   The first input argument is either a real symmetric matrix, a struct with
%   components A.L and A.U or a string containing the name of an M-file which
%   applies a linear operator to the columns of a given matrix.  In the latter
%   case, the second input argument must be N, the order of the problem.
%
%   If A is a struct with components A.L and A.U, such that
%   L*U = (A - sigma*I), a shift-and-invert Lanczos iteration is performed
%
%   The OPTIONS structure is used to control the reorthogonalization:
%     OPTIONS.delta:  Desired level of orthogonality
%                     (default = sqrt(eps/K)).
%     OPTIONS.eta  :  Level of orthogonality after reorthogonalization
%                     (default = eps^(3/4)/sqrt(K)).
%     OPTIONS.cgs  :  Flag for switching between different reorthogonalization
%                     algorithms:
%                      0 = iterated modified Gram-Schmidt  (default)
%                      1 = iterated classical Gram-Schmidt
%     OPTIONS.elr  :  If OPTIONS.elr = 1 (default) then extended local
%                     reorthogonalization is enforced.
%     OPTIONS.Y    :  The lanczos vectors are reorthogonalized against
%                     the columns of the matrix OPTIONS.Y.
%
%   If both R0, Q_old and T_old are provided, they must contain
%   a partial Lanczos tridiagonalization of A on the form
%
%        A Q_old = Q_old T_old + R0 .
%
%   In this case the factorization is extended to dimension K x K by
%   continuing the Lanczos algorithm with R0 as starting vector.
%
%   On exit ANORM contains an approximation to ||A||_2.
%     IERR = 0  :  K steps were performed succesfully.
%     IERR > 0  :  K steps were performed succesfully, but the algorithm
%                  switched to full reorthogonalization after IERR steps.
%     IERR < 0  :  Iteration was terminated after -IERR steps because an
%                  invariant subspace was found, and 3 deflation attempts
%                  were unsuccessful.
%   On exit WORK(1) contains the number of reorthogonalizations performed, and
%   WORK(2) contains the number of inner products performed in the
%   reorthogonalizations.
%
%   See also LANEIG, REORTH, COMPUTE_INT

% References:
% R.M. Larsen, Ph.D. Thesis, Aarhus University, 1998.
%
% G. H. Golub & C. F. Van Loan, "Matrix Computations",
% 3. Ed., Johns Hopkins, 1996.  Chapter 9.
%
% B. N. Parlett, ``The Symmetric Eigenvalue Problem'',
% Prentice-Hall, Englewood Cliffs, NJ, 1980.
%
% H. D. Simon, ``The Lanczos algorithm with partial reorthogonalization'',
% Math. Comp. 42 (1984), no. 165, 115--142.

% Rasmus Munk Larsen, DAIMI, 1998

% Check input arguments.
if nargin<1, error('Not enough input arguments.');  end
if isnumeric(A) | isstruct(A)
  if isnumeric(A)
    [m n] = size(A);
    if m~=n | ~isequal(A,A') | ~isreal(A)
      error('A must be real symmetric')
    end
  elseif isstruct(A)
    [m n] = size(A.L);
  end

  if nargin<7 | isempty(T_k),
    anorm = []; est_anorm=1;
  else
    anorm = T_k; est_anorm=0;
  end
  if nargin<6,  Q_k=[]; T_k=[]; else,  T_k = Q_k; Q_k = options; end
  if nargin<4 | isempty(r),  options = []; else,  options = r;  end
  if nargin<3 | isempty(kmax),
    r = rand(n,1)-0.5;
  else
    r = kmax;
  end
  if nargin<2 | isempty(nin);  kmax = max(10,n/10); else,  kmax = nin;  end
else
  if nargin<2
    error('Not enough input arguments.');
  end
  % Check input functions and parse to create an internal object
  % if an explicit expression is given.
  [A, msg] = fcnchk(A);
  if ~isempty(msg)
    error(msg);
  end
  n = nin;
  if nargin<8 | isempty(anorm), anorm = []; est_anorm=1; else est_anorm=0; end
  if nargin<7,  Q_k=[]; T_k=[]; end
  if nargin<5 | isempty(options),  options = [];          end
  if nargin<4 | isempty(r),  r = rand(n,1)-0.5;   end
  if nargin<3 | isempty(kmax);  kmax = max(10,n/10); end
end

% Set options.
delta = sqrt(eps/kmax); % Desired level of orthogonality.
eta = eps^(3/4)/sqrt(kmax);     % Level of orth. after reorthogonalization.
cgs = 0;                % Flag for switching between iterated CGS and MGS.
elr = 1;                % Flag for switching extended local
                        % reorthogonalization on and off.
deflate = 0;              % Flag for deflation against OPTIONS.Y

% Parse options struct
if ~isempty(options) & isstruct(options)
  c = fieldnames(options);
  for i=1:length(c)
    if strmatch(c(i),'delta'), delta = getfield(options,'delta');  end
    if strmatch(c(i),'eta'), eta = getfield(options,'eta'); end
    if strmatch(c(i),'cgs'), cgs = getfield(options,'cgs'); end
    if strmatch(c(i),'elr'), elr = getfield(options,'elr'); end
    if strmatch(c(i),'Y'), deflate = ~isempty(options.Y);  end
  end
end

np = 0;  nr = 0; ierr=0;

% Rule-of-thumb estimate on the size of round-off terms:
eps1 = sqrt(n)*eps/2; % Notice that {\bf u} == eps/2.
gamma = 1/sqrt(2);

% Prepare Lanczos iteration
if isempty(Q_k) % New Lanczos tridiagonalization.
  % Allocate space
  alpha = zeros(kmax+1,1);  beta = zeros(kmax+1,1);
  Q_k = zeros(n,kmax);
  q = zeros(n,1); beta(1)=norm(r);
  omega = zeros(kmax,1); omega_max = omega;  omega_old = omega;
  omega(1) = 0;   force_reorth= 0;
  j0 = 1;
else            % Extending existing Lanczos tridiagonalization.
  j = size(Q_k,2); % Size of existing factorization
  % Allocate space
  Q_k = [Q_k zeros(n,kmax-j)];
  alpha = zeros(kmax+1,1);  beta = zeros(kmax+1,1);
  alpha(1:j) = diag(T_k);
  if j>1
    beta(2:j) = diag(T_k,-1);
  end
  q = Q_k(:,j);
  % Reorthogonalize r.
  beta(j+1) = norm(r);
  if j<kmax & beta(j+1)*delta < anorm*eps1,
    fro = 1;
  end
  if isfinite(delta)
    int = 1:j;
    [r,beta(j+1),rr] = reorth(Q_k,r,beta(j+1),int,gamma,cgs);
    np = rr*j;    nr = 1;   force_reorth = 1;
  else
     force_reorth = 0;
  end
  % Compute Gerscgorin bound on ||T_k||_2 as SQRT(||T_k'*T_k||_1)
  if est_anorm
    anorm = sqrt(norm(T_k'*T_k,1));
  end
  omega = eps1*ones(kmax,1); omega_max = omega;  omega_old = omega;
  j0 = j+1;
end

if delta==0
  fro = 1; % The user has requested full reorthogonalization.
else
  fro = 0;
end

for j=j0:kmax,
  % Lanczos Step:
  q_old = q;
  if beta(j)==0
    q = r;
  else
    q = r / beta(j);
  end
  Q_k(:,j) = q;
  if isnumeric(A)
    u = A*q;
  elseif isstruct(A)
    u = A.U \ ( A.L \ q);
  else
    u = feval(A,q);
  end
  r = u - beta(j)*q_old;
  alpha(j) = q'*r;
  r = r - alpha(j)*q;

  % Extended local reorthogonalization:
  beta(j+1) = sqrt(r'*r); % Quick and dirty estimate.
  if beta(j+1)<gamma*beta(j) & elr
    if  j==1
      t1=0;
      for i=1:2
	t = q'*r;
	r = r-q*t;
	t1 = t1+t;
      end
      alpha(j) = alpha(j) + t1;
    elseif j>1
      t1 = q_old'*r;
      t2 = q'*r;
      r = r  - (q_old*t1 + q*t2); % Add small terms together first to
      if beta(j)~=0               % reduce risk of cancellation.
	beta(j) = beta(j) + t1;
      end
      alpha(j) = alpha(j) + t2;
    end
    beta(j+1) = sqrt(r'*r); % Quick and dirty estimate.
  end

  % Update Gersgorin estimate of ||T_k|| if required
%  if est_anorm & beta(j+1)~=0
%    T_k = spdiags([[beta(2:j);0] alpha(1:j) beta(1:j)],-1:1,j,j);
%    anorm = sqrt(norm(T_k'*T_k,1))
%  end
  if  est_anorm & beta(j+1)~=0
    anorm = update_gbound(anorm,alpha,beta,j);
  end

  % Update omega-recurrence
  if j>1 & ~fro & beta(j+1)~=0
    [omega,omega_old] = update_omega(omega,omega_old,j,alpha,beta,...
	eps1,anorm);
    omega_max(j) = max(abs(omega));
  end

  % Reorthogonalize if required
  if j>1 & (fro  | force_reorth | omega_max(j)>delta) & beta(j+1)~=0
    if fro
      int = 1:j;
    else
      if force_reorth == 0
	force_reorth= 1; % Do forced reorth to avoid spill-over from q_{j-1}.
	int = compute_int(omega,j,delta,eta,0,0,0);
      else
	force_reorth= 0;
      end
    end
    [r,beta(j+1),rr] = reorth(Q_k,r,beta(j+1),int,gamma,cgs);
    omega(int) = eps1;
    np = np + rr*length(int(:));    nr = nr + 1;
  else
    beta(j+1) = norm(r); % compute norm accurately.
  end

  if deflate
    [r,beta(j+1),rr] = reorth(options.Y,r,beta(j+1),1:size(options.Y,2), ...
			      gamma,cgs);
  end

  if  j<kmax & beta(j+1) < n*anorm*eps  ,
    % If beta is "small" we deflate by setting the off-diagonals of T_k
    % to 0 and attempt to restart with a basis for a new
    % invariant subspace by replacing r with a random starting vector:
    beta(j+1) = 0;
    bailout = 1;
    for attempt=1:3
      r = rand(n,1)-0.5;
      if isnumeric(A)
	r = A*r;
      elseif isstruct(A)
	r = A.U \ ( A.L \ r);
      else
	r = feval(A,r);
      end
      nrm=sqrt(r'*r); % not necessary to compute the norm accurately here.
      int = 1:j;
      [r,nrmnew,rr] = reorth(Q_k,r,nrm,int,gamma,cgs);
      omega(int) = eps1;
      np = np + rr*length(int(:));    nr = nr + 1;
      if nrmnew > 0
	% A vector numerically orthogonal to span(Q_k(:,1:j)) was found.
	% Continue iteration.
	bailout=0;
	break;
      end
    end
    if bailout
      ierr = -j;
      break;
    else
      r=r/nrmnew; % Continue with new normalized r as starting vector.
      force_reorth = 1;
      if delta>0
	fro = 0;    % Turn off full reorthogonalization.
      end
    end
  elseif j<kmax & ~fro & beta(j+1)*delta < anorm*eps1,
    % If anorm*eps1/beta(j+1) > delta then  omega(j+1) will
    % immediately exceed delta, and thus forcing a reorth. to occur at the
    % next step. The components of omega will mainly be determined
    % by the initial value and not the recurrence, and therefore we
    % cannot tell reliably which components exceed eta => we might
    % as well switch to full reorthogonalization to avoid trouble.
    % The user is probably trying to determine pathologically
    % small ( < sqrt(eps)*||A||_2 ) eigenvalues.
    %    warning(['Semiorthogonality cannot be maintained at iteration ', ...
    %	  num2str(j),'. The matrix is probably ill-conditioned.', ...
    %	  ' Switching to full reorthogonalization.'])
    fro = 1;
    ierr = j;
  end
end

% Set up tridiagonal T_k in sparse matrix data structure.
T_k = spdiags([[beta(2:j);0] alpha(1:j) beta(1:j)],-1:1,j,j);
if nargout<2
  Q_k = T_k;
elseif j~=size(Q_k,2)
  Q_k = Q_k(:,1:j);
end
work = [nr np];

function [omega,omega_old] = update_omega(omega, omega_old, j, ...
    alpha,beta,eps1,anorm)
% UPDATE_OMEGA:  Update Simon's omega_recurrence for the Lanczos vectors.
%
% [omega,omega_old] = update_omega(omega, omega_old,j,eps1,alpha,beta,anorm)
%

% Rasmus Munk Larsen, DAIMI, 1998.

% Estimate of contribution to roundoff errors from A*v
%   fl(A*v) = A*v + f,
% where ||f|| \approx eps1*||A||.
% For a full matrix A, a rule-of-thumb estimate is eps1 = sqrt(n)*eps.
T = eps1*anorm;
binv = 1/beta(j+1);

omega_old = omega;
% Update omega(1) using omega(0)==0.
omega_old(1)= beta(2)*omega(2)+ (alpha(1)-alpha(j))*omega(1) -  ...
    beta(j)*omega_old(1);
omega_old(1) = binv*(omega_old(1) + sign(omega_old(1))*T);
% Update remaining components.
k=2:j-2;
omega_old(k) = beta(k+1).*omega(k+1) + (alpha(k)-alpha(j)).*omega(k) ...
     + beta(k).*omega(k-1) - beta(j)*omega_old(k);
omega_old(k) = binv*(omega_old(k) + sign(omega_old(k))*T);
omega_old(j-1) = binv*T;
% Swap omega and omega_old.
temp = omega;
omega = omega_old;
omega_old = omega;
omega(j) =  eps1;

function anorm = update_gbound(anorm,alpha,beta,j)
%UPDATE_GBOUND   Update Gerscgorin estimate of 2-norm
%  ANORM = UPDATE_GBOUND(ANORM,ALPHA,BETA,J) updates the Gerscgorin bound
%  for the tridiagonal in the Lanczos process after the J'th step.
%  Applies Gerscgorins circles to T_K'*T_k instead of T_k itself
%  since this gives a tighter bound.

if j==1 % Apply Gerscgorin circles to T_k'*T_k to estimate || A ||_2
  i=j;
  % scale to avoid overflow
  scale = max(abs(alpha(i)),abs(beta(i+1)));
  alpha(i) = alpha(i)/scale;
  beta(i+1) = beta(i+1)/scale;
  anorm = 1.01*scale*sqrt(alpha(i)^2+beta(i+1)^2 + abs(alpha(i)*beta(i+1)));
elseif j==2
  i=1;
  % scale to avoid overflow
  scale = max(max(abs(alpha(1:2)),max(abs(beta(2:3)))));
  alpha(1:2) = alpha(1:2)/scale;
  beta(2:3) = beta(2:3)/scale;

  anorm = max(anorm, scale*sqrt(alpha(i)^2+beta(i+1)^2 + ...
      abs(alpha(i)*beta(i+1) + alpha(i+1)*beta(i+1)) + ...
      abs(beta(i+1)*beta(i+2))));
  i=2;
  anorm = max(anorm,scale*sqrt(abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
      beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
      abs(alpha(i)*beta(i+1))) );
elseif j==3
  % scale to avoid overflow
  scale = max(max(abs(alpha(1:3)),max(abs(beta(2:4)))));
  alpha(1:3) = alpha(1:3)/scale;
  beta(2:4) = beta(2:4)/scale;
  i=2;
  anorm = max(anorm,scale*sqrt(abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
      beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
      abs(alpha(i)*beta(i+1) + alpha(i+1)*beta(i+1)) + ...
      abs(beta(i+1)*beta(i+2))) );
  i=3;
  anorm = max(anorm,scale*sqrt(abs(beta(i)*beta(i-1)) + ...
      abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
      beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
      abs(alpha(i)*beta(i+1))) );
else
  % scale to avoid overflow
  %  scale = max(max(abs(alpha(j-2:j)),max(abs(beta(j-2:j+1)))));
  %  alpha(j-2:j) = alpha(j-2:j)/scale;
  %  beta(j-2:j+1) = beta(j-2:j+1)/scale;

  % Avoid scaling, which is slow. At j>3 the estimate is usually quite good
  % so just make sure that anorm is not made infinite by overflow.
  i = j-1;
  anorm1 = sqrt(abs(beta(i)*beta(i-1)) + ...
      abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
      beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
      abs(alpha(i)*beta(i+1) + alpha(i+1)*beta(i+1)) + ...
      abs(beta(i+1)*beta(i+2)));
  if isfinite(anorm1)
    anorm = max(anorm,anorm1);
  end
  i = j;
  anorm1 = sqrt(abs(beta(i)*beta(i-1)) + ...
      abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
      beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
      abs(alpha(i)*beta(i+1)));
  if isfinite(anorm1)
    anorm = max(anorm,anorm1);
  end
end

author: andrewssobral

repo: lrslibrary


text:
function dz = dynamics_6_link(~,z,P)
%DZ = DYNAMICS_6_LINK(T,Z,P)
%
%FUNCTION:  This function computes the dynamics of a 6
%    link pendulum, and is designed to be called from ode45.
%    The model allows for arbitrary mass and inertia for each
%    link, but no friction or actuation
%
%INPUTS:
%    t = time. Dummy input for ode45. Not used.
%    z = [12 X nTime]  matrix of states
%    P = struct of parameters
%OUTPUTS:
%    dz = [12 X nTime]  matrix of state derivatives
%
%NOTES:
%    This file was automatically generated by writeDynamics.m

g  = P.g ; %gravity
m1 = P.m(1); % Link 1 mass
m2 = P.m(2); % Link 2 mass
m3 = P.m(3); % Link 3 mass
m4 = P.m(4); % Link 4 mass
m5 = P.m(5); % Link 5 mass
m6 = P.m(6); % Link 6 mass
l1 = P.l(1); % Link 1 length
l2 = P.l(2); % Link 2 length
l3 = P.l(3); % Link 3 length
l4 = P.l(4); % Link 4 length
l5 = P.l(5); % Link 5 length
l6 = P.l(6); % Link 6 length
I1 = P.I(1); % Link 1 moment of inertia about its center of mass
I2 = P.I(2); % Link 2 moment of inertia about its center of mass
I3 = P.I(3); % Link 3 moment of inertia about its center of mass
I4 = P.I(4); % Link 4 moment of inertia about its center of mass
I5 = P.I(5); % Link 5 moment of inertia about its center of mass
I6 = P.I(6); % Link 6 moment of inertia about its center of mass
d1 = P.d(1); % Link 1 distance between center of mass and parent joint
d2 = P.d(2); % Link 2 distance between center of mass and parent joint
d3 = P.d(3); % Link 3 distance between center of mass and parent joint
d4 = P.d(4); % Link 4 distance between center of mass and parent joint
d5 = P.d(5); % Link 5 distance between center of mass and parent joint
d6 = P.d(6); % Link 6 distance between center of mass and parent joint

nTime = size(z,2);
dz = zeros(size(z));
M = zeros(6,6,nTime);
f = zeros(6,nTime);

th1 = z(1,:);
th2 = z(2,:);
th3 = z(3,:);
th4 = z(4,:);
th5 = z(5,:);
th6 = z(6,:);

dth1 = z(7,:);
dth2 = z(8,:);
dth3 = z(9,:);
dth4 = z(10,:);
dth5 = z(11,:);
dth6 = z(12,:);

dz(1,:) = dth1;
dz(2,:) = dth2;
dz(3,:) = dth3;
dz(4,:) = dth4;
dz(5,:) = dth5;
dz(6,:) = dth6;

M(1,1,:) = - I1 - d1.^2.*m1 - l1.^2.*m2 - l1.^2.*m3 - l1.^2.*m4 - l1.^2.*m5 - l1.^2.*m6;
M(1,2,:) = -l1.*cos(th1 - th2).*(d2.*m2 + l2.*m3 + l2.*m4 + l2.*m5 + l2.*m6);
M(1,3,:) = -l1.*cos(th1 - th3).*(d3.*m3 + l3.*m4 + l3.*m5 + l3.*m6);
M(1,4,:) = -l1.*cos(th1 - th4).*(d4.*m4 + l4.*m5 + l4.*m6);
M(1,5,:) = -l1.*cos(th1 - th5).*(d5.*m5 + l5.*m6);
M(1,6,:) = -d6.*l1.*m6.*cos(th1 - th6);
M(2,1,:) = -l1.*cos(th1 - th2).*(d2.*m2 + l2.*m3 + l2.*m4 + l2.*m5 + l2.*m6);
M(2,2,:) = - I2 - d2.^2.*m2 - l2.^2.*m3 - l2.^2.*m4 - l2.^2.*m5 - l2.^2.*m6;
M(2,3,:) = -l2.*cos(th2 - th3).*(d3.*m3 + l3.*m4 + l3.*m5 + l3.*m6);
M(2,4,:) = -l2.*cos(th2 - th4).*(d4.*m4 + l4.*m5 + l4.*m6);
M(2,5,:) = -l2.*cos(th2 - th5).*(d5.*m5 + l5.*m6);
M(2,6,:) = -d6.*l2.*m6.*cos(th2 - th6);
M(3,1,:) = -l1.*cos(th1 - th3).*(d3.*m3 + l3.*m4 + l3.*m5 + l3.*m6);
M(3,2,:) = -l2.*cos(th2 - th3).*(d3.*m3 + l3.*m4 + l3.*m5 + l3.*m6);
M(3,3,:) = - I3 - d3.^2.*m3 - l3.^2.*m4 - l3.^2.*m5 - l3.^2.*m6;
M(3,4,:) = -l3.*cos(th3 - th4).*(d4.*m4 + l4.*m5 + l4.*m6);
M(3,5,:) = -l3.*cos(th3 - th5).*(d5.*m5 + l5.*m6);
M(3,6,:) = -d6.*l3.*m6.*cos(th3 - th6);
M(4,1,:) = -l1.*cos(th1 - th4).*(d4.*m4 + l4.*m5 + l4.*m6);
M(4,2,:) = -l2.*cos(th2 - th4).*(d4.*m4 + l4.*m5 + l4.*m6);
M(4,3,:) = -l3.*cos(th3 - th4).*(d4.*m4 + l4.*m5 + l4.*m6);
M(4,4,:) = - I4 - d4.^2.*m4 - l4.^2.*m5 - l4.^2.*m6;
M(4,5,:) = -l4.*cos(th4 - th5).*(d5.*m5 + l5.*m6);
M(4,6,:) = -d6.*l4.*m6.*cos(th4 - th6);
M(5,1,:) = -l1.*cos(th1 - th5).*(d5.*m5 + l5.*m6);
M(5,2,:) = -l2.*cos(th2 - th5).*(d5.*m5 + l5.*m6);
M(5,3,:) = -l3.*cos(th3 - th5).*(d5.*m5 + l5.*m6);
M(5,4,:) = -l4.*cos(th4 - th5).*(d5.*m5 + l5.*m6);
M(5,5,:) = - I5 - d5.^2.*m5 - l5.^2.*m6;
M(5,6,:) = -d6.*l5.*m6.*cos(th5 - th6);
M(6,1,:) = -d6.*l1.*m6.*cos(th1 - th6);
M(6,2,:) = -d6.*l2.*m6.*cos(th2 - th6);
M(6,3,:) = -d6.*l3.*m6.*cos(th3 - th6);
M(6,4,:) = -d6.*l4.*m6.*cos(th4 - th6);
M(6,5,:) = -d6.*l5.*m6.*cos(th5 - th6);
M(6,6,:) = - I6 - d6.^2.*m6;

f(1,:) = - d1.*g.*m1.*cos(th1) - g.*l1.*m2.*cos(th1) - g.*l1.*m3.*cos(th1) - g.*l1.*m4.*cos(th1) - g.*l1.*m5.*cos(th1) - g.*l1.*m6.*cos(th1) - d2.*dth2.^2.*l1.*m2.*sin(th1 - th2) - d3.*dth3.^2.*l1.*m3.*sin(th1 - th3) - d4.*dth4.^2.*l1.*m4.*sin(th1 - th4) - d5.*dth5.^2.*l1.*m5.*sin(th1 - th5) - d6.*dth6.^2.*l1.*m6.*sin(th1 - th6) - dth2.^2.*l1.*l2.*m3.*sin(th1 - th2) - dth2.^2.*l1.*l2.*m4.*sin(th1 - th2) - dth2.^2.*l1.*l2.*m5.*sin(th1 - th2) - dth2.^2.*l1.*l2.*m6.*sin(th1 - th2) - dth3.^2.*l1.*l3.*m4.*sin(th1 - th3) - dth3.^2.*l1.*l3.*m5.*sin(th1 - th3) - dth3.^2.*l1.*l3.*m6.*sin(th1 - th3) - dth4.^2.*l1.*l4.*m5.*sin(th1 - th4) - dth4.^2.*l1.*l4.*m6.*sin(th1 - th4) - dth5.^2.*l1.*l5.*m6.*sin(th1 - th5);
f(2,:) = d2.*dth1.^2.*l1.*m2.*sin(th1 - th2) - g.*l2.*m3.*cos(th2) - g.*l2.*m4.*cos(th2) - g.*l2.*m5.*cos(th2) - g.*l2.*m6.*cos(th2) - d2.*g.*m2.*cos(th2) - d3.*dth3.^2.*l2.*m3.*sin(th2 - th3) - d4.*dth4.^2.*l2.*m4.*sin(th2 - th4) - d5.*dth5.^2.*l2.*m5.*sin(th2 - th5) - d6.*dth6.^2.*l2.*m6.*sin(th2 - th6) + dth1.^2.*l1.*l2.*m3.*sin(th1 - th2) + dth1.^2.*l1.*l2.*m4.*sin(th1 - th2) + dth1.^2.*l1.*l2.*m5.*sin(th1 - th2) + dth1.^2.*l1.*l2.*m6.*sin(th1 - th2) - dth3.^2.*l2.*l3.*m4.*sin(th2 - th3) - dth3.^2.*l2.*l3.*m5.*sin(th2 - th3) - dth3.^2.*l2.*l3.*m6.*sin(th2 - th3) - dth4.^2.*l2.*l4.*m5.*sin(th2 - th4) - dth4.^2.*l2.*l4.*m6.*sin(th2 - th4) - dth5.^2.*l2.*l5.*m6.*sin(th2 - th5);
f(3,:) = d3.*dth1.^2.*l1.*m3.*sin(th1 - th3) - g.*l3.*m4.*cos(th3) - g.*l3.*m5.*cos(th3) - g.*l3.*m6.*cos(th3) - d3.*g.*m3.*cos(th3) + d3.*dth2.^2.*l2.*m3.*sin(th2 - th3) - d4.*dth4.^2.*l3.*m4.*sin(th3 - th4) - d5.*dth5.^2.*l3.*m5.*sin(th3 - th5) - d6.*dth6.^2.*l3.*m6.*sin(th3 - th6) + dth1.^2.*l1.*l3.*m4.*sin(th1 - th3) + dth1.^2.*l1.*l3.*m5.*sin(th1 - th3) + dth1.^2.*l1.*l3.*m6.*sin(th1 - th3) + dth2.^2.*l2.*l3.*m4.*sin(th2 - th3) + dth2.^2.*l2.*l3.*m5.*sin(th2 - th3) + dth2.^2.*l2.*l3.*m6.*sin(th2 - th3) - dth4.^2.*l3.*l4.*m5.*sin(th3 - th4) - dth4.^2.*l3.*l4.*m6.*sin(th3 - th4) - dth5.^2.*l3.*l5.*m6.*sin(th3 - th5);
f(4,:) = d4.*dth1.^2.*l1.*m4.*sin(th1 - th4) - g.*l4.*m5.*cos(th4) - g.*l4.*m6.*cos(th4) - d4.*g.*m4.*cos(th4) + d4.*dth2.^2.*l2.*m4.*sin(th2 - th4) + d4.*dth3.^2.*l3.*m4.*sin(th3 - th4) - d5.*dth5.^2.*l4.*m5.*sin(th4 - th5) - d6.*dth6.^2.*l4.*m6.*sin(th4 - th6) + dth1.^2.*l1.*l4.*m5.*sin(th1 - th4) + dth1.^2.*l1.*l4.*m6.*sin(th1 - th4) + dth2.^2.*l2.*l4.*m5.*sin(th2 - th4) + dth2.^2.*l2.*l4.*m6.*sin(th2 - th4) + dth3.^2.*l3.*l4.*m5.*sin(th3 - th4) + dth3.^2.*l3.*l4.*m6.*sin(th3 - th4) - dth5.^2.*l4.*l5.*m6.*sin(th4 - th5);
f(5,:) = d5.*dth1.^2.*l1.*m5.*sin(th1 - th5) - g.*l5.*m6.*cos(th5) - d5.*g.*m5.*cos(th5) + d5.*dth2.^2.*l2.*m5.*sin(th2 - th5) + d5.*dth3.^2.*l3.*m5.*sin(th3 - th5) + d5.*dth4.^2.*l4.*m5.*sin(th4 - th5) - d6.*dth6.^2.*l5.*m6.*sin(th5 - th6) + dth1.^2.*l1.*l5.*m6.*sin(th1 - th5) + dth2.^2.*l2.*l5.*m6.*sin(th2 - th5) + dth3.^2.*l3.*l5.*m6.*sin(th3 - th5) + dth4.^2.*l4.*l5.*m6.*sin(th4 - th5);
f(6,:) = d6.*m6.*(dth1.^2.*l1.*sin(th1 - th6) - g.*cos(th6) + dth2.^2.*l2.*sin(th2 - th6) + dth3.^2.*l3.*sin(th3 - th6) + dth4.^2.*l4.*sin(th4 - th6) + dth5.^2.*l5.*sin(th5 - th6));

for i=1:nTime
    MM = M(:,:,i);  ff = f(:,i);
    dz(7:12,i) = -MM \ ff;
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function L = evaluate_log_posterior(this, uv)
%EVALUATE_LOG_POSTERIOR computes the log-posterior (negative energy) of the
%   flow fields UV
%   Actually only proportional to the log posterior since the variance of neither the
%   spatial nor the data terms is considered
%
%   This is a member function of the class 'alt_ba_optical_flow'.
%
% Authors: Deqing Sun, Department of Computer Science, Brown University
% Contact: dqsun@cs.brown.edu
% $Date: 2009 $
% $Revision: $
%
% Copyright 2009-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

% spatial term
S = this.spatial_filters;
p = 0;

for i = 1:length(S)

    u_ = conv2(uv(:,:,1), S{i}, 'valid');
    v_ = conv2(uv(:,:,2), S{i}, 'valid');

    if isa(this.rho_spatial_u{i}, 'robust_function')

        p = p - sum(evaluate(this.rho_spatial_u{i}, u_(:)))...
            - sum(evaluate(this.rho_spatial_v{i}, v_(:)));

    elseif isa(this.rho_spatial_u{i}, 'gsm_density')

        p   = p + sum(evaluate_log(this.rho_spatial_u{i}, u_(:)'))...
                    + sum(evaluate_log(this.rho_spatial_v{i}, v_(:)'));

    else
        error('evaluate_log_posterior: unknown rho function!');
    end;
end;

% likelihood
It  = partial_deriv(this.images, uv, this.interpolation_method);

if isa(this.rho_data, 'robust_function')

    l   = -sum(evaluate(this.rho_data, It(:)));

elseif isa(this.rho_data, 'gsm_density')

    l   = sum(evaluate_log(this.rho_data, It(:)'));

else
    error('evaluate_log_posterior: unknown rho function!');
end;

L = this.lambda*p + l;

if this.display
    fprintf('spatial\t%3.2e\tdata\t%3.2e\n', this.lambda*p, l);
end;

author: kristinbranson

repo: JAABA


text:
function [dimensionsout, datatypesout]=byte2any(dimensions, datatypes)

% BYTE2ANY - Mask Initialization function for UDP Byte Un-Pack Block
% Derived from The MathWorks, Inc. MXPCREADME
% $Revision: 1 $ $Date: 2008/03/07 $

if ~isa(dimensions,'cell')
  error('dimensions argument must be a cell array');
end

if ~isa(datatypes,'cell')
  error('datatypes argument must be a cell array');
end

if length(dimensions)~=length(datatypes)
  error('dimension cell array and datatypes cell array must have the same length');
end

dimensionsout=[];
datatypesout=[];
for i=1:length(dimensions)
  dimension=dimensions{i};
  if ~isa(dimension,'double')
    error('dimensions elements must be double arrays');
  end
  if size(dimension,1)>1 | size(dimension,2)>2
    error('dimensions elements must be either scalar or a row vector with two elements');
  end
  if length(dimension)==1
    dimensionsout=[dimensionsout,dimension,0];
  else
    dimensionsout=[dimensionsout,dimension];
  end
  datatype=datatypes{i};
  if ~isa(datatype,'char')
    error('datatypes elements must be char arrays');
  end
  if strcmp(datatype,'double')
     dtypeout=0;
  elseif strcmp(datatype,'single')
    dtypeout=1;
  elseif strcmp(datatype,'int8')
    dtypeout=2;
  elseif strcmp(datatype,'uint8')
    dtypeout=3;
  elseif strcmp(datatype,'int16')
    dtypeout=4;
  elseif strcmp(datatype,'uint16')
    dtypeout=5;
  elseif strcmp(datatype,'int32')
    dtypeout=6;
  elseif strcmp(datatype,'uint32')
    dtypeout=7;
  elseif strcmp(datatype,'boolean')
    dtypeout=8;
  else
    error('the supported data types are: double, single, int, uint8, int16, uint16, int32, uint32, boolean');
  end
  datatypesout=[datatypesout,dtypeout];
end


text:
function [ y, m, d, ierror ] = ymd_check_gregorian ( y, m, d )

%*****************************************************************************80
%
%% YMD_CHECK_GREGORIAN checks a Gregorian YMD date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, D, the YMD date, which may
%    be corrected if necessary and possible.
%
%    Output, integer IERROR, is 0 if the date is legal.
%

%
%  Check the month.
%
  [ y, m, ierror ] = ym_check_gregorian ( y, m );

  if ( ierror ~= 0 )
    return
  end
%
%  Check the day.
%
  [ y, m, d ] = day_borrow_gregorian ( y, m, d );

  [ y, m, d ] = day_carry_gregorian ( y, m, d );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function linplus_test0192 ( )

%*****************************************************************************80
%
%% TEST0192 tests R8CC_WRITE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = 5;
  nz_num = 12;

  a_file = 'r8cc_a.txt';
  col = [ 1, 4, 6, 8, 10, 13 ];
  col_file = 'r8cc_col.txt';
  row = [ 1, 2, 4, 1, 2, 3, 5, 4, 5, 1, 2, 5 ];
  row_file = 'r8cc_row.txt';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0192\n' );
  fprintf ( 1, '  For a matrix in the R8CC format,\n' );
  fprintf ( 1, '  (double precision compressed column sparse)\n' );
  fprintf ( 1, '  R8CC_WRITE writes the matrix to 3 files.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix rows M     = %d\n', m );
  fprintf ( 1, '  Matrix columns N  = %d\n', n );
  fprintf ( 1, '  Nonzeros NZ_NUM   = %d\n', nz_num );

  i4vec_print ( n+1, col, '  The COL vector:' );

  i4vec_print ( nz_num, row, '  The ROW vector:' );

  a = r8cc_indicator ( m, n, nz_num, col, row );

  r8cc_print ( m, n, nz_num, col, row, a, '  The R8CC matrix:' );

  r8cc_write ( col_file, row_file, a_file, m, n, nz_num, col, row, a );

  return
end


text:
function [VOLLOCS, LOCS] = tapas_physio_create_scan_timing_nominal(t, ...
    sqpar, align_scan, durationPhyslogAfterEndOfLastScan)
% Creates locations of scan volume and slice events in time vector of SCANPHYSLOG-files
%
%   [VOLLOCS, LOCS] = tapas_physio_create_scan_timing_nominal(t, sqpar);
%
% In cases where the SCANPHYSLOG-file has no gradient entries (column
% 7-9), the actual time-course of the sequence has to be inferred from the
% nominal sequence parameters in sqpar. Here, the corresponding slice scan
% events are generated for the resampling of regressors in the GLM under
% the assumption that the SCANPHYSLOG-file ended exactly when the scan ended
% ...which is usually the case if a scan is not stopped manually
%
% Additionally, one can set a buffer end time, i.e., the duration of the
% phys logging lasting longer after the end of the last scan
%
% IN
%   t       - timing vector of SCANPHYSLOG-file, usually sampled with 500
%             Hz (Philips)
%   sqpar                   - sequence timing parameters
%           .Nslices        - number of slices per volume in fMRI scan
%           .NslicesPerBeat - usually equals Nslices, unless you trigger with the heart beat
%           .TR             - repetition time in seconds
%           .Ndummies       - number of dummy volumes
%           .Nscans         - number of full volumes saved (volumes in nifti file,
%                             usually rows in your design matrix)
%           .Nprep          - number of non-dummy, volume like preparation pulses
%                             before 1st dummy scan. If set, logfile is read from beginning,
%                             otherwise volumes are counted from last detected volume in the logfile
%           .time_slice_to_slice - time between the acquisition of 2 subsequent
%                             slices; typically TR/Nslices or
%                             minTR/Nslices, if minimal temporal slice
%                             spacing was chosen
%           .onset_slice    - slice whose scan onset determines the adjustment of the
%                             regressor timing to a particular slice for the whole volume
%   align_scan              'first' or 'last' (default)
%                           'first' t == 0 will be aligned to first scan
%                                   volume, first slice
%                           'last'  t(end) will be aligned to last scan
%                                   volume, last slice
%   durationPhyslogAfterEndOfLastScan
%                            duration (in seconds) of physiological logfile
%                            after end of last scan volume in the run
%                           default: 0
% OUT
%           VOLLOCS         - locations in time vector, when volume scan
%                             events started
%           LOCS            - locations in time vector, when slice or volume scan
%                             events started
% EXAMPLE
%   [VOLLOCS, LOCS] = tapas_physio_create_scan_timing_nominal(t, sqpar);
%
%   See also

% Author: Lars Kasper
% Created: 2013-02-07
% Copyright (C) 2013 Institute for Biomedical Engineering, ETH/Uni Zurich.
%
% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

if nargin < 3
    align_scan = 'last';
end

if nargin < 4
    durationPhyslogAfterEndOfLastScan = 0;
end

Nscans          = sqpar.Nscans;
Ndummies        = sqpar.Ndummies;

NallVols = (Ndummies+Nscans);
VOLLOCS = NaN(NallVols,1);
TR = sqpar.TR;

%% First, find volume starts either forward or backward through time series
do_count_from_start = strcmpi(align_scan, 'first');
if do_count_from_start % t = 0 is assumed to be the start of the scan
    for n = 1:NallVols
        [tmp, VOLLOCS(n)] = min(abs(t - TR*(n-1)));
    end
else
    tEndLastScan = t(end)-durationPhyslogAfterEndOfLastScan;

    tStartPhys = t(1);
    for n = 1:NallVols

        tStartVol = (tEndLastScan-TR*(NallVols-n+1));

        if tStartPhys > tStartVol
            VOLLOCS(n) = NaN;
        else
            [tmp, VOLLOCS(n)] = min(abs(t - tStartVol));
        end
    end
end

%% Then, find slice starts between determined volume starts

LOCS = tapas_physio_create_LOCS_from_VOLLOCS(VOLLOCS, t, sqpar);

end

author: translationalneuromodeling

repo: tapas


text:
function [x,state] = struct_rational(z,task,t)
%STRUCT_RATIONAL Matrix with columns as rational functions.
%   [x,state] = struct_rational(z,[],t) computes a matrix x in which the
%   jth column is equal to the rational function
%
%      polyval(z{1}(j,:),s)./polyval([z{2}(j,:) 1],s)
%
%   evaluated at the points s, defined as
%
%      (t-0.5*(min(t)+max(t)))/(0.5*(max(t)-min(t))).
%
%   The degree of the numerator and denomator are equal to size(z{1},2)-1
%   and size(z{2},2), respectively. The structure state stores information
%   which is reused in computing the right and left Jacobian-vector
%   products.
%
%   struct_rational(z,task,t) computes the right or left Jacobian-vector
%   product of this transformation, depending on the structure task. Use
%   the structure state and add the field 'r' of the same shape as z or the
%   field 'l' of the same shape as x to obtain the structure task for
%   computing the right and left Jacobian-vector products
%
%      (dF(:)/dz(:).')*task.r(:) and
%      (dF(:)/dz(:).')'*task.l(:) + conj((dF(:)/dconj(z(:)).')'*task.l(:)),
%
%   respectively. Here, F(z) represents this transormation, (:) signifies
%   vectorization and the derivative w.r.t. z (conj(z)) is a partial
%   derivative which treats conj(z) (z) as constant. The output has the
%   same shape as x or z for the right and left Jacobian-vector products,
%   respectively.
%
%   See also struct_poly, struct_rbf.

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

if nargin < 2, task = []; end
if nargin < 3
    error('struct_rational:t','Please supply evaluation points.');
end

% Center and scale evaluations points to fit in [-1,1].
mx = max(t);
mn = min(t);
t = (t-0.5*(mn+mx))/(0.5*(mx-mn));

if isempty(task) || (isempty(task.l) && isempty(task.r))
    x = polyval(z{1},t(:));
    if ~isempty(z{2})
        state.denom = polyval([z{2} ones(size(z{2},1),1)],t(:));
        x = x./state.denom;
        state.deriv = -x./state.denom;
    else
        state = [];
    end
elseif ~isempty(task.r)
    x = polyval(task.r{1},t(:));
    if ~isempty(z{2})
        x = x./task.denom+ ...
            task.deriv.*polyval([task.r{2} zeros(size(z{2},1),1)],t(:));
    end
    state = [];
elseif ~isempty(task.l)
    x = {zeros(size(z{1})),zeros(size(z{2}))};
    tmp = ones(length(t),1);
    if isempty(z{2})
        for i = size(z{1},2):-1:1
            x{1}(:,i) = sum(bsxfun(@times,conj(tmp),task.l),1);
            tmp = tmp.*t(:);
        end
    else
        for i = size(z{1},2):-1:1
            x{1}(:,i) = sum(conj(bsxfun(@rdivide,tmp,task.denom)).* ...
                task.l,1);
            tmp = tmp.*t(:);
        end
        tmp = t(:);
        for i = size(z{2},2):-1:1
            x{2}(:,i) = sum(conj(bsxfun(@times,tmp,task.deriv)).*task.l,1);
            tmp = tmp.*t(:);
        end
    end
    state = [];
end

end

function y = polyval(p,t)
    y = p(:,ones(length(t),1)).';
    for i = 2:size(p,2)
        y = bsxfun(@plus,bsxfun(@times,y,t),p(:,i).');
    end
end


text:
function dummy = makeDummyModel(numMets, numRxns)
% Makes an empty model with numMets rows for metabolites and numRxns columns
% for reactions. Includes all fields that are necessary to join models.
%
% USAGE:
%
%   dummy = makeDummyModel(numMets, numRxns)
%
% INPUTS:
%   numMets:       Number of metabolites
%   numRxns:       Number of reactions
%
% OUTPUT:
%   dummy:         Empty COBRA model structure
%
% .. Authors: SM June 2016

dummy.mets = cell(numMets, 1);
dummy.rxns = cell(numRxns, 1);
dummy.c = zeros(numRxns, 1);
dummy.S = zeros(numMets, numRxns);
dummy.lb = zeros(numRxns, 1);
dummy.ub = zeros(numRxns, 1);
dummy.b = zeros(numMets, 1);
dummy.metNames = cell(numMets, 1);
dummy.metFormulas = cell(numMets, 1);
dummy.rxnNames = cell(numRxns, 1);
dummy.subSystems = cell(numRxns, 1);
dummy.genes = {};
dummy.rxnGeneMat = sparse(numRxns, size(dummy.genes, 1));
dummy.rules = cell(numRxns, 1);
dummy.grRules = cell(numRxns, 1);
dummy.comments = cell(numRxns, 1);
dummy.citations = cell(numRxns, 1);
dummy.confidenceScores = cell(numRxns, 1);
dummy.ecNumbers = cell(numRxns, 1);
dummy.rxnKeggID = cell(numRxns, 1);
dummy.metCHEBIID = cell(numMets, 1);
dummy.metCharge = zeros(numMets, 1);
dummy.metFormulas = cell(numMets, 1);
dummy.metHMDB = cell(numMets, 1);
dummy.metInchiString = cell(numMets, 1);
dummy.metKeggID = cell(numMets, 1);
dummy.metSmile = cell(numMets, 1);
dummy.metPubChemID = cell(numMets, 1);
dummy.description.organism = {}; dummy.description.author = {}; dummy.description.name = {};
dummy.description.geneindex = {}; dummy.description.genedate = {}; dummy.description.genesource = {}; dummy.description.notes = {};
dummy.disabled = {};

author: opencobra

repo: cobratoolbox


text:
%findSubfunctionTests Utility script used for subfunction-based tests
%   This file is a script that is called at the top of M-files containing
%   subfunction-based tests.
%
%   The top of a typical M-file using this script looks like this:
%
%       function test_suite = testFeatureA
%
%       findSubfunctionTests;
%
%   IMPORTANT NOTE
%   --------------
%   The output variable name for an M-file using this script must be test_suite.

%   Steven L. Eddins
%   Copyright 2008-2009 The MathWorks, Inc.

[ST,I] = dbstack;
caller_name = ST(I + 1).name;
subFcns = which('-subfun', caller_name);

setup_fcn_name = subFcns(xunit.utils.isSetUpString(subFcns));
if numel(setup_fcn_name) > 1
    error('findSubfunctionTests:tooManySetupFcns', ...
        'Found more than one setup subfunction.')
elseif isempty(setup_fcn_name)
    setup_fcn = [];
else
    setup_fcn = str2func(setup_fcn_name{1});
end

teardown_fcn_name = subFcns(xunit.utils.isTearDownString(subFcns));
if numel(teardown_fcn_name) > 1
    error('findSubfunctionTests:tooManyTeardownFcns', ...
        'Found more than one teardown subfunction.')
elseif isempty(teardown_fcn_name)
    teardown_fcn = [];
else
    teardown_fcn = str2func(teardown_fcn_name{1});
end

test_fcns = cellfun(@str2func, subFcns(xunit.utils.isTestString(subFcns)), ...
    'UniformOutput', false);

suite = TestSuite;
suite.Name = caller_name;
suite.Location = which(caller_name);
for k = 1:numel(test_fcns)
    suite.add(FunctionHandleTestCase(test_fcns{k}, setup_fcn, teardown_fcn));
end

if nargout > 0
    test_suite = suite;
else
    suite.run();
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
%{
load('dataset/trafficdb/traffic_patches.mat');
[M,m,n,p] = convert_video3d_to_2d(im2double(imgdb{100}));
out = run_algorithm('MC', 'SVP', M, [])
show_results(M.*out.Omega,out.L,out.S,out.O,p,m,n);
%}

warning('off','all');
MIdx = M(Idx);
rank = 2;
[U,S,V] = svp(Idx,MIdx,size(M,1),size(M,2),rank);
L = U*diag(S)*V'; % low-rank
S = (M - L); % sparse
warning('on','all');


text:
function nlp  = LoadProblem(system, bounds, gait, load_path)

    if nargin < 4
        load_path = [];
    end

    num_grid.RightStance1 = 10;
    num_grid.LeftStance1 = 10;
    num_grid.RightStance2 = 10;
    num_grid.LeftStance2 = 10;

    options = {'EqualityConstraintBoundary', 1e-4,...
        'DistributeTimeVariable', false,...
        'DistributeParameters',false};

    nlp = HybridTrajectoryOptimization('atrias_opt', system, num_grid, [], options{:});

    if isempty(load_path)
        nlp.configure(bounds);
    else
        nlp.configure(bounds, 'LoadPath',load_path);
    end
    % cost function
    %     opt.cost.Power(nlp, system);
    trans_opt.cost.StateDeviation(nlp, system, gait);
    nlp.update;

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of the
%   FANUC LR MATE 10iA, FANUC Robotics Europe.
%   Author: David Garcia Munyoz, Nestor Gomez Lopez, Teresa Pomares
%   Palorames.
%   Universidad Miguel Hernandez de Elche.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function q = inversekinematic_fanuc_mate(robot, T)
%initialize q,
%eight possible solutions are generally feasible
q=zeros(6,8);

%Evaluate the parameters
d = eval(robot.DH.d);

%See geometry at the reference for this robot
L6=d(6);

%T= [ nx ox ax Px;
%     ny oy ay Py;
%     nz oz az Pz];
Px=T(1,4);
Py=T(2,4);
Pz=T(3,4);

%Compute the position of the wrist, being W the Z component of the end effector's system
W = T(1:3,3);

% Pm: wrist position
Pm = [Px Py Pz]' - L6*W;

%first joint, two possible solutions admited:
% if q(1) is a solution, then q(1) + pi is also a solution
q1=atan2(Pm(2), Pm(1));

%solve for q2
q2_1=solve_for_theta2(robot, [q1 0 0 0 0 0 0], Pm);

q2_2=solve_for_theta2(robot, [q1+pi 0 0 0 0 0 0], Pm);

%solve for q3
q3_1=solve_for_theta3(robot, [q1 0 0 0 0 0 0], Pm);

q3_2=solve_for_theta3(robot, [q1+pi 0 0 0 0 0 0], Pm);

%Arrange solutions, there are 8 possible solutions so far.
% if q1 is a solution, q1* = q1 + pi is also a solution.
% For each (q1, q1*) there are two possible solutions
% for q2 and q3 (namely, elbow up and elbow up solutions)
% So far, we have 4 possible solutions. Howefer, for each triplet (theta1, theta2, theta3),
% there exist two more possible solutions for the last three joints, generally
% called wrist up and wrist down solutions. For this reason,
%the next matrix doubles each column. For each two columns, two different
%configurations for theta4, theta5 and theta6 will be computed. These
%configurations are generally referred as wrist up and wrist down solution
q = [q1         q1         q1        q1       q1+pi   q1+pi   q1+pi   q1+pi;
     q2_1(1)    q2_1(1)    q2_1(2)   q2_1(2)  q2_2(1) q2_2(1) q2_2(2) q2_2(2);
     q3_1(1)    q3_1(1)    q3_1(2)   q3_1(2)  q3_2(1) q3_2(1) q3_2(2) q3_2(2);
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0];

%leave only the real part of the solutions
q=real(q);

%Note that in this robot, the joint q3 has a non-simmetrical range. In this
%case, the joint ranges from 60 deg to -219 deg, thus, the typical normalizing
%step is avoided in this angle (the next line is commented). When solving
%for the orientation, the solutions are normalized to the [-pi, pi] range
%only for the theta4, theta5 and theta6 joints.

%normalize q to [-pi, pi]
q(1,:) = normalize(q(1,:));
q(2,:) = normalize(q(2,:));
% solve for the last three joints
% for any of the possible combinations (theta1, theta2, theta3)
for i=1:2:size(q,2),
    qtemp = solve_spherical_wrist(robot, q(:,i), T, 1,'geometric'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i)=qtemp;

    qtemp = solve_spherical_wrist(robot, q(:,i), T, -1, 'geometric'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i+1)=qtemp;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for second joint theta2, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q2 = solve_for_theta2(robot, q, Pm)

%Evaluate the parameters
d = eval(robot.DH.d);
a = eval(robot.DH.a);

%See geometry
L2=a(2);
L3=d(4);

A3=a(3); % desfase

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Eslab�n equivalente para simplificar el desfase
l3 = sqrt (A3^2 + L3^2);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

beta = atan2(-p1(2), p1(1));
gamma = (acos((L2^2+r^2-l3^2)/(2*r*L2)));

if ~isreal(gamma)
    disp('WARNING:inversekinematic_fanuc_mate: the point is not reachable for this configuration, imaginary solutions');
    %gamma = real(gamma);
end

%return two possible solutions
%elbow up and elbow down
%the order here is important and is coordinated with the function
%solve_for_theta3
q2(1) = pi/2 - beta - gamma; %elbow up
q2(2) = pi/2 - beta + gamma; %elbow down

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for third joint theta3, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q3 = solve_for_theta3(robot, q, Pm)

%Evaluate the parameters
d = eval(robot.DH.d);
a = eval(robot.DH.a);

%See geometry
L2=a(2);
L3=d(4);

A3= a(3); %desfase

%See geometry of the robot
%compute L4
l3 = sqrt(A3^2 + L3^2);

%the angle phi is fixed
phi=acos((A3^2+l3^2-L3^2)/(2*A3*l3));

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

eta = (acos((L2^2 + l3^2 - r^2)/(2*L2*l3)));

if ~isreal(eta)
   disp('WARNING:inversekinematic_fanuc_mate: the point is not reachable for this configuration, imaginary solutions');
   %eta = real(eta);
end

%return two possible solutions
%elbow up and elbow down solutions
%the order here is important
q3(1) = pi - phi - eta;
q3(2) = pi - phi + eta;

%remove complex solutions for q for the q1+pi solutions
function  qreal = arrange_solutions(q)
qreal=q(:,1:4);

%sum along rows if any angle is complex, for any possible solutions, then v(i) is complex
v = sum(q, 1);

for i=5:8,
    if isreal(v(i))
        qreal=[qreal q(:,i)]; %store the real solutions
    end
end

qreal = real(qreal);

repo: ARTE


text:
function M_hat = run_mc(params)
  M = params.M;
  Idx = params.Idx;
  maxrank = 1;
  maxCycles = 100;
  step_size = 0.1;
  [numr,numc] = size(M);
  [Usg,Vsg,err_reg] = grouse(M,Idx,numr,numc,maxrank,step_size,maxCycles);
  M_hat = Usg*Vsg';
end


text:
%% Run a PREP report on an individual dataset -- good for producing figures.

%% Specific examples
% dataFile = 'N:\ARLAnalysis\VEP\VEPRobust_1Hz_Post_Median_Unfiltered\vep_03.set';
% summaryFolder = 'D:\\Temp';
% basename = 'vep';
% publishOn = false;

% dataFile = 'N:\ARLAnalysis\VEP\VEPRobust_1Hz_Post_Median_Unfiltered\vep_15.set';
% summaryFolder = 'D:\\Temp';
% basename = 'vep';
% publishOn = false;

% dataFile = 'N:\ARLAnalysis\VEP\VEPRobust_1Hz_Post_Median_Unfiltered\vep_02.set';
% summaryFolder = 'D:\\Temp';
% basename = 'vep';
% publishOn = false;

% dataFile = 'N:\\ARLAnalysis\\NCTU\NCTURobust_1Hz_New\\session\\67\\eeg_studyLevel2_NCTU_Lane-Keeping_Task_session_67_subject_1_task_motionless_s54_081209n_recording_1.set';
% summaryFolder = 'D:\\Temp';
% basename = 'nctu-lk';
% publishOn = false;

% dataFile = 'N:\ARLAnalysis\NCTU\NCTURobust_1Hz_New\session\79\eeg_studyLevel2_NCTU_Lane-Keeping_Task_session_79_subject_1_task_motion_s77_120504m_recording_1.set';
% summaryFolder = 'D:\\Temp';
% basename = 'nctu-lk';
% publishOn = false;

% dataFile = 'N:\\ARLAnalysis\\NCTU\\NCTURobust_1Hz_New\\session\\80\\eeg_studyLevel2_NCTU_Lane-Keeping_Task_session_80_subject_1_task_motion_s80_120522m_recording_1.set';
% summaryFolder = 'D:\\Temp';
% basename = 'nctu-lk';
% publishOn = false;

% dataFile = 'N:\ARLAnalysis\Shooter\Shooter_Robust_1Hz_Unfiltered\\ARIT_0131_CNT_events.set';
% summaryFolder = 'D:\\Temp';
% basename = 'shooter-1';
% publishOn = false;

dataFile = 'N:\\ARLAnalysis\\RSVP_HeadIT\\RSVP_Robust_1Hz\\rsvp_01.set';
summaryFolder = 'D:\\Temp';
basename = 'rsvp-1';
publishOn = true;

%% Setup up the summary folder
summaryReportName = [basename '_summary.html'];
sessionFolder = '.';
reportSummary = [summaryFolder filesep summaryReportName];
if exist(reportSummary, 'file')
   delete(reportSummary);
end

%% Run the pipeline
[pathstr, fname, ext] = fileparts(dataFile);
sessionReportName = [fname '.pdf'];
load(dataFile, '-mat');
tempReportLocation = [summaryFolder filesep sessionFolder ...
    filesep 'prepPipelineReport.pdf'];
actualReportLocation = [summaryFolder filesep sessionFolder ...
    filesep sessionReportName];
summaryReportLocation = [summaryFolder filesep summaryReportName];
summaryFile = fopen(summaryReportLocation, 'a+', 'n', 'UTF-8');
relativeReportLocation = [sessionFolder filesep sessionReportName];
consoleFID = 1;
publishPrepReport(EEG, reportSummary, actualReportLocation, 1, publishOn);

%% Calculate the mean correlations
meanBefore = EEG.etc.noiseDetection.reference.noisyStatisticsOriginal.maximumCorrelations;
meanBeforeInterp = EEG.etc.noiseDetection.reference.noisyStatisticsBeforeInterpolation.maximumCorrelations;
meanAfter = EEG.etc.noiseDetection.reference.noisyStatistics.maximumCorrelations;
fprintf('Mean before: %g\n', mean(meanBefore(:)));
fprintf('mean before Interpolation: %g\n', mean(meanBeforeInterp(:)));
fprintf('Mean after: %g\n', mean(meanAfter(:)));

author: VisLab

repo: EEG-Clean-Tools


text:
%--- help for find_posterior_mode ---
%
%  INTERNAL FUNCTION: Finds the mode of the posterior
%
%  Note:
%     - Though one can start MCMC from any point in theory, a good starting point with a good estimate of the covariance/Hessian is requirement to create a good proposal distribution for Bayesian analysis of DSGE models. This function finds the mode of the posterior for a candidate starting point of the Bayesian analysis.
%     - Though available for modularity of the RISE toolbox, one should use :func:`estimate <dsge.estimate>` for most use cases.
%
%

author: jmaih


text:
function uI = faceinterpolate(u,node,elem,elemType)
%% FACEINTERPOLATE interpolate to face elements (RT0 or BDM1).
%
% uI = faceinterpolate(u,node,elem,elemType) interpolates a given function u
% into the lowesr order RT0 or BDM1 finite element spaces. The coefficient
% is given by the line integral int_e u*n ds. The input elemType can be 'RT0'
% or 'BDM1'.
%
% Example
%
%    maxIt = 5;
%    node = [1,0; 1,1; 0,1; -1,1; -1,0; -1,-1; 0,-1; 0,0]; % nodes
%    elem = [1,2,8; 3,8,2; 8,3,5; 4,5,3; 7,8,6; 5,6,8];    % elements
%    bdEdge = setboundary(node,elem,'Dirichlet');
%    pde = mixBCdata;
%    err = zeros(maxIt,2); N = zeros(maxIt,1);
%    for i =1:maxIt
%      [node,elem,bdEdge] = uniformrefine(node,elem,bdEdge);
%      [u,sigma,M] = PoissonRT0(node,elem,pde,bdEdge);
%      sigmaI = faceinterpolate(pde.Du,node,elem,'RT0');
%      err(i,1) = getL2errorRT0(node,elem,pde.Du,sigmaI,pde.d);
%      err(i,2)=sqrt((sigma-sigmaI)'*M*(sigma-sigmaI));
%      N(i) = size(u,1);
%    end
%     figure;
%     showrate2(N,err(:,1),2,'r-+','||Du - \sigma_I||',...
%               N,err(:,2),2,'b-+','||\sigma^{RT_0} - \sigma_I||');
%
% See also edgeinterpolate, edgeinterpolate1, edgeinterpolate2
%
% Created by Ming Wang at Mar 29, 2011, M-lint modified at May 14, 2011.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('elemType','var'), elemType = 'RT0'; end
%% edge
if size(elem,2) == 2 % the input elem is edge
    edge = elem;
else
    [tempvar,edge] = dofedge(elem);
end
NE = size(edge,1);
edgeVec = node(edge(:,2),:) - node(edge(:,1),:);
nVec = zeros(NE,2);
nVec(:,1) = edgeVec(:,2);
nVec(:,2) = -edgeVec(:,1);

%% dof for RT0
[lambda,weight] = quadpts1(4);
nQuad = size(lambda,1);
uI = zeros(NE,1);
for i = 1:nQuad
    pxy = lambda(i,1)*node(edge(:,1),:)+lambda(i,2)*node(edge(:,2),:);
    flux = u(pxy);
    uI = uI + weight(i)*dot(flux,nVec,2);
end

%% dof for BDM1
if strcmp(elemType,'BDM1')
   uI(NE+1:2*NE) = zeros(NE,1);
   for i = 1:nQuad
        pxy = lambda(i,1)*node(edge(:,1),:)+lambda(i,2)*node(edge(:,2),:);
        flux = u(pxy);
        uI(NE+1:2*NE) = uI(NE+1:2*NE)+ ...
                     weight(i)*3*(lambda(i,1)-lambda(i,2))*dot(flux,nVec,2);
   end
end

repo: Aero-matlab