text:
function gpcf = gpcf_mask(varargin)
%GPCF_MASK  Create a mask covariance function
%
%  Description
%    GPCF = GPCF_MASK('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates a mask covariance function structure in
%    which the named parameters have the specified values. Any
%    unspecified parameters are set to default values.
%
%    GPCF = GPCF_MASK(GPCF,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a covariance function structure with the named
%    parameters altered with the specified values.
%
%    Mask covariance function returns correlation 1 if input
%    values X_i and X_j are both non zero and 0 otherwise.
%
%    Parameters for mask covariance function
%      selectedVariables - vector defining which inputs are used
%
%  See also
%    GP_SET, GPCF_*, PRIOR_*, MEAN_*

% Copyright (c) 2007-2010 Jarno Vanhatalo
% Copyright (c) 2008-2010 Jaakko Riihim�ki
% Copyright (c) 2010 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'GPCF_MASK';
  ip.addOptional('gpcf', [], @isstruct);
  ip.addParamValue('selectedVariables',[], ...
                   @(x) isempty(x) || (isvector(x) && all(x>0)));
  ip.parse(varargin{:});
  gpcf=ip.Results.gpcf;

  if isempty(gpcf)
    % Initialize a covariance function structure
    init=true;
    gpcf.type = 'gpcf_mask';
  else
    % Modify a covariance function structure
    if ~isfield(gpcf,'type') && ~isequal(gpcf.type,'gpcf_mask')
      error('First argument does not seem to be a valid covariance function structure')
    end
    init=false;
  end

  if ~ismember('selectedVariables',ip.UsingDefaults)
    if ~isempty(ip.Results.selectedVariables)
      gpcf.selectedVariables = ip.Results.selectedVariables;
    elseif isfield(gpcf,'selectedVariables')
      gpcf=rmfield(gpcf,'selectedVariables');
    end
  end
  if init
    gpcf.fh.pak = @gpcf_mask_pak;
    gpcf.fh.unpak = @gpcf_mask_unpak;
    gpcf.fh.lp = @gpcf_mask_lp;
    gpcf.fh.lpg = @gpcf_mask_lpg;
    gpcf.fh.cfg = @gpcf_mask_cfg;
    gpcf.fh.ginput = @gpcf_mask_ginput;
    gpcf.fh.cov = @gpcf_mask_cov;
    gpcf.fh.trcov  = @gpcf_mask_trcov;
    gpcf.fh.trvar  = @gpcf_mask_trvar;
    gpcf.fh.recappend = @gpcf_mask_recappend;
  end

end

function [w,s] = gpcf_mask_pak(gpcf, w)
%GPCF_MASK_PAK  Combine GP covariance function parameters into
%              one vector.
%
%  Description
%    W = GPCF_MASK_PAK(GPCF) takes a covariance function
%    structure GPCF and combines the covariance function
%    parameters and their hyperparameters into a single row
%    vector W. This is a mandatory subfunction used for
%    example in energy and gradient computations.
%
%       w = []
%
%  See also
%    GPCF_MASK_UNPAK

  w = []; s = {};
end

function [gpcf, w] = gpcf_mask_unpak(gpcf, w)
%GPCF_MASK_UNPAK  Sets the covariance function parameters into
%                the structure
%
%  Description
%    [GPCF, W] = GPCF_MASK_UNPAK(GPCF, W) takes a covariance
%    function structure GPCF and a parameter vector W, and
%    returns a covariance function structure identical to the
%    input, except that the covariance parameters have been set
%    to the values in W. Deletes the values set to GPCF from W
%    and returns the modified W. This is a mandatory subfunction
%    used for example in energy and gradient computations.
%
%    Assignment is inverse of
%       w = []
%
%  See also
%   GPCF_MASK_PAK

end

function lp = gpcf_mask_lp(gpcf)
%GPCF_MASK_LP  Evaluate the energy of prior of covariance function parameters
%
%  Description
%    LP = GPCF_MASK_LP(GPCF) takes a covariance function
%    structure GPCF and returns log(p(th)), where th collects the
%    parameters. This is a mandatory subfunction used for example
%    in energy computations.
%
%  See also
%    GPCF_MASK_PAK, GPCF_MASK_UNPAK, GPCF_MASK_LPG, GP_E

  lp = 0;

end

function lpg = gpcf_mask_lpg(gpcf)
%GPCF_MASK_LPG  Evaluate gradient of the log prior with respect
%               to the parameters.
%
%  Description
%    LPG = GPCF_MASK_LPG(GPCF) takes a covariance function
%    structure GPCF and returns LPG = d log (p(th))/dth, where th
%    is the vector of parameters. This is a mandatory subfunction
%    used for example in gradient computations.
%
%  See also
%    GPCF_MASK_PAK, GPCF_MASK_UNPAK, GPCF_MASK_LP, GP_G

  lpg = [];

end

function DKff = gpcf_mask_cfg(gpcf, x, x2, mask,i1)
%GPCF_MASK_CFG  Evaluate gradient of covariance function
%               with respect to the parameters.
%
%  Description
%    DKff = GPCF_MASK_CFG(GPCF, X) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X) with
%    respect to th (cell array with matrix elements). This is a
%    mandatory subfunction used in gradient computations.
%
%    DKff = GPCF_MASK_CFG(GPCF, X, X2) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X2) with
%    respect to th (cell array with matrix elements). This subfunction
%    is needed when using sparse approximations (e.g. FIC).
%
%    DKff = GPCF_MASK_CFG(GPCF, X, [], MASK) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the diagonal of gradients of covariance matrix
%    Kff = k(X,X2) with respect to th (cell array with matrix
%    elements). This subfunction is needed when using sparse
%    approximations (e.g. FIC).
%
%    DKff = GPCF_MASK_CFG(GPCF, X, X2, [], i) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X2), or k(X,X)
%    if X2 is empty, with respect to ith hyperparameter. This subfunction
%    is needed when using memory save option in gp_set.
%
%  See also
%    GPCF_MASK_PAK, GPCF_MASK_UNPAK, GPCF_MASK_LP, GP_G

  DKff = {};

end

function [DKff, lpg]  = gpcf_mask_ginput(gpcf, x, x2, i1)
%GPCF_MASK_GINPUT  Evaluate gradient of covariance function with
%                  respect to x.
%
%  Description
%    DKff = GPCF_MASK_GINPUT(GPCF, X) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X) with
%    respect to X (cell array with matrix elements). This subfunction
%    is needed when computing gradients with respect to inducing
%    inputs in sparse approximations.
%
%    DKff = GPCF_MASK_GINPUT(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors
%    and returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2) with respect to X (cell array with matrix elements).
%    This subfunction is needed when computing gradients with
%    respect to inducing inputs in sparse approximations.
%
%    DKff = GPCF_MASK_GINPUT(GPCF, X, X2, i) takes a covariance
%    function structure GPCF, a matrix X of input vectors
%    and returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2), or k(X,X) if X2 is empty, with respect to ith covariate
%    in X. This subfunction is needed when using memory saving option
%    in gp_set.
%
%  See also
%   GPCF_MASK_PAK, GPCF_MASK_UNPAK, GPCF_MASK_LP, GP_G

  [n, m] =size(x);

  if nargin == 2 || isempty(x2)
    ii1 = 0;
    for i=1:m
      for j = 1:n
        ii1 = ii1 + 1;
        DKff{ii1} = zeros(n);
        lpg(ii1) = 0;
      end
    end
  elseif nargin == 3 || nargin == 4
    ii1 = 0;
    for i=1:m
      for j = 1:n
        ii1 = ii1 + 1;
        DKff{ii1} = zeros(n, size(x2,1));
        lpg(ii1) = 0;
      end
    end
  end
  if nargin==4
    DKff=DKff{1};
  end
end

function C = gpcf_mask_cov(gpcf, x1, x2, varargin)
%GP_MASK_COV  Evaluate covariance matrix between two input vectors
%
%  Description
%    C = GP_MASK_COV(GP, TX, X) takes in covariance function of a
%    Gaussian process GP and two matrixes TX and X that contain
%    input vectors to GP. Returns covariance matrix C. Every
%    element ij of C contains covariance between inputs i in TX
%    and j in X. This is a mandatory subfunction used for example in
%    prediction and energy computations.
%
%  See also
%    GPCF_MASK_TRCOV, GPCF_MASK_TRVAR, GP_COV, GP_TRCOV

  if isempty(x2)
    x2=x1;
  end
  [n1,m1]=size(x1);
  [n2,m2]=size(x2);

  if m1~=m2
    error('the number of columns of X1 and X2 has to be same')
  end

  C=repmat(true,n1,n2);
  if isfield(gpcf, 'selectedVariables')
    for j = 1:length(gpcf.selectedVariables)
      jj=gpcf.selectedVariables(j);
      C = C & bsxfun(@and,x1(:,jj)~=0,x2(:,jj)'~=0);
    end
  else
    for j = 1:m1
      C = C & bsxfun(@and,x1(:,j)~=0,x2(:,j)'~=0);
    end
  end
  C=double(C);

end

function C = gpcf_mask_trcov(gpcf, x)
%GP_MASK_TRCOV  Evaluate training covariance matrix of inputs
%
%  Description
%    C = GP_MASK_TRCOV(GP, TX) takes in covariance function of a
%    Gaussian process GP and matrix TX that contains training
%    input vectors. Returns covariance matrix C. Every element ij
%    of C contains covariance between inputs i and j in TX. This
%    is a mandatory subfunction used for example in prediction
%    and energy computations.
%
%  See also
%    GPCF_MASK_COV, GPCF_MASK_TRVAR, GP_COV, GP_TRCOV

  [n,m]=size(x);

  C=repmat(true,n,n);
  if isfield(gpcf, 'selectedVariables')
    for j = 1:length(gpcf.selectedVariables)
      jj=gpcf.selectedVariables(j);
      C = C & bsxfun(@and,x(:,jj)~=0,x(:,jj)'~=0);
    end
  else
    for j = 1:m
      C = C & bsxfun(@and,x(:,j)~=0,x(:,j)'~=0);
    end
  end
  C=double(C);

end

function C = gpcf_mask_trvar(gpcf, x)
%GP_MASK_TRVAR  Evaluate training variance vector
%
%  Description
%    C = GP_MASK_TRVAR(GPCF, TX) takes in covariance function of a
%    Gaussian process GPCF and matrix TX that contains training
%    inputs. Returns variance vector C. Every element i of C
%    contains variance of input i in TX. This is a mandatory
%    subfunction used for example in prediction and energy
%    computations.
%
%  See also
%    GPCF_MASK_COV, GP_COV, GP_TRCOV

  [n,m]=size(x);
  C=true(n,1);
  if isfield(gpcf, 'selectedVariables')
    for j = 1:length(gpcf.selectedVariables)
      jj=gpcf.selectedVariables(j);
      C = C & x(:,jj)~=0;
    end
  else
    for j = 1:m
      C = C & x(:,j)~=0;
    end
  end
  C=double(C);

end

function reccf = gpcf_mask_recappend(reccf, ri, gpcf)
%RECAPPEND  Record append
%
%  Description
%    RECCF = GPCF_MASK_RECAPPEND(RECCF, RI, GPCF) takes a
%    covariance function record structure RECCF, record index RI
%    and covariance function structure GPCF with the current MCMC
%    samples of the parameters. Returns RECCF which contains
%    all the old samples and the current samples from GPCF.
%    This subfunction is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC and GP_MC -> RECAPPEND

  if nargin == 2
    % Initialize the record
    reccf.type = 'gpcf_mask';

    % Initialize parameters
    reccf.coeffSigma2= [];

    % Set the function handles
    reccf.fh.pak = @gpcf_mask_pak;
    reccf.fh.unpak = @gpcf_mask_unpak;
    reccf.fh.lp = @gpcf_mask_lp;
    reccf.fh.lpg = @gpcf_mask_lpg;
    reccf.fh.cfg = @gpcf_mask_cfg;
    reccf.fh.cov = @gpcf_mask_cov;
    reccf.fh.trcov  = @gpcf_mask_trcov;
    reccf.fh.trvar  = @gpcf_mask_trvar;
    reccf.fh.recappend = @gpcf_mask_recappend;
  else
    % Append to the record
    if isfield(gpcf, 'selectedVariables')
      reccf.selectedVariables = gpcf.selectedVariables;
    end
  end
end


text:
% Test file for trigtech/isnan.m

function pass = test_isnan(pref)

% Get preferences.
if ( nargin < 1 )
    pref = trigtech.techPref();
end

testclass = trigtech();

% Test a scalar-valued function.
f = testclass.make(@(x) cos(pi*x), [], pref);
pass(1) = ~isnan(f);

% Test an array-valued function.
f = testclass.make(@(x) [cos(pi*x), cos(pi*x).^2], [], pref);
pass(2) = ~isnan(f);

% Artificially construct and test a NaN-valued function.
f = testclass.make(NaN);
pass(3) = isnan(f);

% Test a NaN scalar-valued function.
try
    f = testclass.make(@(x) cos(pi*x) + NaN);
    pass(4) = isnan(f);
catch ME
    pass(4) = strcmpi(ME.message, 'Cannot handle functions that evaluate to Inf or NaN.');
end

% Test a NaN array-valued function.
try
    f = testclass.make(@(x) [cos(pi*x) + NaN, cos(pi*x)]);
    pass(5) = isnan(f);
catch ME
    pass(5) = strcmpi(ME.message, 'Cannot handle functions that evaluate to Inf or NaN.');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% DEMSILHOUETTEPLOT

% GP

% Show prediction for test data.
yPred = modelOut(model, XTest);
xyzankurAnimCompare(yPred, yTest, 96);

yDiff = (yPred - yTest);
rmsError = sqrt(sum(sum(yDiff.*yDiff))/numel(yDiff));

counter = 0;
if printDiagram
  ind = 1:27:size(yPred, 1);
  for i = ind
    counter = counter + 1;
    figure
    handle = xyzankurVisualise(yPred(i,:), 1);
    printPlot([fileBaseName '_' num2str(counter)], '../tex/diagrams', '../html')
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef IWT2HD < dagnn.ElementWise

  properties (Transient)
      padding = 0
      wavename = 'haart'
      opts = {}
    numInputs
  end

  methods
    function outputs = forward(obj, inputs, params)
      obj.numInputs = numel(inputs) ;
      outputs{1} = vl_nniwt2(inputs{1}, [], ...
          'wavename', obj.wavename, 'padding', obj.padding, obj.opts{:}) ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      derInputs{1} =  vl_nniwt2(inputs{1},  derOutputs{1}, ...
          'wavename', obj.wavename, 'padding', obj.padding, obj.opts{:}) ;
      derParams = {0} ;
    end
    function outputSizes = getOutputSizes(obj, inputSizes)
      outputSizes = {} ;
    end

    function rfs = getReceptiveFields(obj)
        rfs = [] ;
    end

    function obj = IWT2HD(varargin)
      obj.load(varargin) ;
    end
  end
end

repo: MWCNN


text:
function [L] = lap2DDir(N, h);
%
% [L] = lap2DDir(N, h)
%
%  Constructs the 2D Laplacian for a square mesh with zero Dirichlet
%     boundary conditions using the standard 5-point stencil.
%
%  Returns:
%     L = discrete Laplacian
%
%  Input:
%     N = number of mesh points in each direction
%     h = mesh width
%
%
%
%  License: This code is free to use for any purposes, provided
%           any publications resulting from the use of this code
%           reference the original code/author.
%
%  Author:  Samuel Isaacson (isaacson@math.utah.edu)
%  Date:    11/2007
%
%  Please notify the author of any bugs, and contribute any
%  modifications or bug fixes back to the original author.
%
%  Disclaimer:
%   This code is provided as is. The author takes no responsibility
%   for its results or effects.

M = N * N;

% 2D Laplace Operator on Square:
e              = ones(M,1);
ed             = -4*e;
eu1            = ones(M,1);
eu1((N+1):N:M) = 0;
ed1            = ones(M,1);
ed1(N:N:M)     = 0;

L = spdiags([e ed1 ed eu1 e], [-N -1:1 N], M, M);

L = L ./ (h*h);

repo: Aero-matlab


text:
function [taumat, thetamat] ...
         = SimulateControl(thetalist, dthetalist, g, Ftipmat, Mlist, ...
                           Glist, Slist, thetamatd, dthetamatd, ...
                           ddthetamatd, gtilde, Mtildelist, Gtildelist, ...
                           Kp, Ki, Kd, dt, intRes)
% *** CHAPTER 11: ROBOT CONTROL ***
% Takes thetalist: n-vector of initial joint variables,
%       dthetalist: n-vector of initial joint velocities,
%       g: Actual gravity vector g,
%       Ftipmat: An N x 6 matrix of spatial forces applied by the
%                end-effector (If there are no tip forces, the user should
%                input a zero and a zero matrix will be used),
%       Mlist: Actual list of link frames i relative to i? at the home
%              position,
%       Glist: Actual spatial inertia matrices Gi of the links,
%       Slist: Screw axes Si of the joints in a space frame, in the format
%              of a matrix with the screw axes as the columns,
%       thetamatd: An Nxn matrix of desired joint variables from the
%                  reference trajectory,
%       dthetamatd: An Nxn matrix of desired joint velocities,
%       ddthetamatd: An Nxn matrix of desired joint accelerations,
%       gtilde: The gravity vector based on the model of the actual robot
%               (actual values given above),
%       Mtildelist: The link frame locations based on the model of the
%                   actual robot (actual values given above),
%       Gtildelist: The link spatial inertias based on the model of the
%                   actual robot (actual values given above),
%       Kp: The feedback proportional gain (identical for each joint),
%       Ki: The feedback integral gain (identical for each joint),
%       Kd: The feedback derivative gain (identical for each joint),
%       dt: The timestep between points on the reference trajectory.
%       intRes: Integration resolution is the number of times integration
%               (Euler) takes places between each time step. Must be an
%               integer value greater than or equal to 1.
% Returns taumat: An Nxn matrix of the controller commanded joint
%                 forces/torques, where each row of n forces/torques
%                 corresponds to a single time instant,
%         thetamat: An Nxn matrix of actual joint angles.
% The end of this function plots all the actual and desired joint angles.
% Example Usage
%
% clc; clear;
% thetalist = [0.1; 0.1; 0.1];
% dthetalist = [0.1; 0.2; 0.3];
% %Initialize robot description (Example with 3 links)
% g = [0; 0; -9.8];
% M01 = [[1, 0, 0, 0]; [0, 1, 0, 0]; [0, 0, 1, 0.089159]; [0, 0, 0, 1]];
% M12 = [[0, 0, 1, 0.28]; [0, 1, 0, 0.13585]; [-1, 0 ,0, 0]; [0, 0, 0, 1]];
% M23 = [[1, 0, 0, 0]; [0, 1, 0, -0.1197]; [0, 0, 1, 0.395]; [0, 0, 0, 1]];
% M34 = [[1, 0, 0, 0]; [0, 1, 0, 0]; [0, 0, 1, 0.14225]; [0, 0, 0, 1]];
% G1 = diag([0.010267, 0.010267, 0.00666, 3.7, 3.7, 3.7]);
% G2 = diag([0.22689, 0.22689, 0.0151074, 8.393, 8.393, 8.393]);
% G3 = diag([0.0494433, 0.0494433, 0.004095, 2.275, 2.275, 2.275]);
% Glist = cat(3, G1, G2, G3);
% Mlist = cat(3, M01, M12, M23, M34);
% Slist = [[1; 0; 1;      0; 1;     0], ...
%        [0; 1; 0; -0.089; 0;     0], ...
%        [0; 1; 0; -0.089; 0; 0.425]];
% dt = 0.01;
% %Create a trajectory to follow
% thetaend =[pi / 2; pi; 1.5 * pi];
% Tf = 1;
% N = Tf / dt;
% method = 5;
% thetamatd = JointTrajectory(thetalist, thetaend, Tf, N, method);
% dthetamatd = zeros(N, 3);
% ddthetamatd = zeros(N, 3);
% dt = Tf / (N - 1);
% for i = 1: N - 1
%   dthetamatd(i + 1, :) = (thetamatd(i + 1, :) - thetamatd(i, :)) / dt;
%   ddthetamatd(i + 1, :) = (dthetamatd(i + 1, :) ...
%                           - dthetamatd(i, :)) / dt;
% end
% %Possibly wrong robot description (Example with 3 links)
% gtilde = [0.8; 0.2; -8.8];
% Mhat01 = [[1, 0, 0, 0]; [0, 1, 0, 0]; [0, 0, 1, 0.1]; [0, 0, 0, 1]];
% Mhat12 = [[0, 0, 1, 0.3]; [0, 1, 0, 0.2]; [-1, 0 ,0, 0]; [0, 0, 0, 1]];
% Mhat23 = [[1, 0, 0, 0]; [0, 1, 0, -0.2]; [0, 0, 1, 0.4]; [0, 0, 0, 1]];
% Mhat34 = [[1, 0, 0, 0]; [0, 1, 0, 0]; [0, 0, 1, 0.2]; [0, 0, 0, 1]];
% Ghat1 = diag([0.1, 0.1, 0.1, 4, 4, 4]);
% Ghat2 = diag([0.3, 0.3, 0.1, 9, 9, 9]);
% Ghat3 = diag([0.1, 0.1, 0.1, 3, 3, 3]);
% Gtildelist = cat(3, Ghat1, Ghat2, Ghat3);
% Mtildelist = cat(4, Mhat01, Mhat12, Mhat23, Mhat34);
% Ftipmat = ones(N, 6);
% Kp = 20;
% Ki = 10;
% Kd = 18;
% intRes = 8;
% [taumat, thetamat] ...
% = SimulateControl(thetalist, dthetalist, g, Ftipmat, Mlist, Glist, ...
%                 Slist, thetamatd, dthetamatd, ddthetamatd, gtilde, ...
%                 Mtildelist, Gtildelist, Kp, Ki, Kd, dt, intRes);
%

Ftipmat = Ftipmat';
thetamatd = thetamatd';
dthetamatd = dthetamatd';
ddthetamatd = ddthetamatd';
n = size(thetamatd, 2);
taumat = zeros(size(thetamatd));
thetamat = zeros(size(thetamatd));
thetacurrent = thetalist;
dthetacurrent = dthetalist;
eint = zeros(size(thetamatd, 1), 1);
for i=1: n
    taulist ...
    = ComputedTorque(thetacurrent, dthetacurrent, eint, gtilde, ...
                     Mtildelist, Gtildelist, Slist, thetamatd(:, i), ...
                     dthetamatd(:, i), ddthetamatd(:, i), Kp, Ki, Kd);
    for j=1: intRes
        ddthetalist ...
        = ForwardDynamics(thetacurrent, dthetacurrent, taulist, g, ...
                          Ftipmat(:, i), Mlist, Glist, Slist);
        [thetacurrent, dthetacurrent] ...
        = EulerStep(thetacurrent, dthetacurrent, ddthetalist, ...
                    (dt / intRes));
    end
    taumat(:, i) = taulist;
    thetamat(:, i) = thetacurrent;
    eint = eint + (dt*(thetamatd(:, i) - thetacurrent));
end
%Output using matplotlib
links = size(thetamat, 1);
leg = cell(1, 2 * links);
time=0: dt: dt * n - dt;
timed=0: dt: dt * n - dt;
figure
hold on
for i=1: links
    col = rand(1, 3);
    plot(time, (thetamat(i, :)'), '-', 'Color', col)
    plot(timed, (thetamatd(i, :)'), '.', 'Color', col)
    leg{2 * i - 1} = (strcat('ActualTheta', num2str(i)));
    leg{2 * i} = (strcat('DesiredTheta', num2str(i)));
end
title('Plot of Actual and Desired Joint Angles')
xlabel('Time')
ylabel('Joint Angles')
legend(leg, 'Location', 'NorthWest')
taumat = taumat';
thetamat = thetamat';
end

author: ShuoYangRobotics

repo: QuadrupedSim


text:
function [ p, l, u ] = oto_plu ( n )

%*****************************************************************************80
%
%% OTO_PLU returns the PLU factors of the OTO matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 June 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real P(N,N), L(N,N), U(N,N), the PLU factors.
%
  p = zeros ( n, n );
 for j = 1 : n
   for i = 1 : n
      if ( i == j )
        p(i,j) = 1.0;
      else
        p(i,j) = 0.0;
      end
    end
  end

 l = zeros ( n, n );
 for j = 1 : n
   for i = 1 : n
      if ( i == j )
        l(i,j) = 1.0;
      elseif ( i == j + 1 )
        l(i,j) = j / ( j + 1 );
      else
        l(i,j) = 0.0;
      end
    end
  end

  u = zeros ( n, n );
  for j = 1 : n
    for i = 1 : n
      if ( i == j )
        u(i,j) = ( i + 1 ) / i;
      elseif ( i == j - 1 )
        u(i,j) = 1.0;
      else
        u(i,j) = 0.0;
      end
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef CEC2010_F7 < PROBLEM
% <single> <real> <constrained>
% CEC'2010 constrained optimization benchmark problem

%------------------------------- Reference --------------------------------
% R. Mallipeddi and P. N. Suganthan, Problem definitions and evaluation
% criteria for the CEC 2010 competition on constrained real-parameter
% optimization, Nanyang Technological University, Singapore, 2010.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        O;  % Optimal decision vector
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2010.mat'),'Data');
            obj.O = Data{7};
            obj.M = 1;
            if isempty(obj.D); obj.D = 10; end
            obj.D = min(obj.D,length(obj.O));
            obj.lower    = zeros(1,obj.D) - 140;
            obj.upper    = zeros(1,obj.D) + 140;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = 1 + PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopObj = sum(100*(Z(:,1:end-1).^2-Z(:,2:end)).^2+(Z(:,1:end-1)-1).^2,2);
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            Y = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopCon = 0.5 - exp(-0.1*sqrt(mean(Y.^2,2))) - 3*exp(mean(cos(0.1*Y),2)) + exp(1);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%% t_gifti_3
%
% Introduction to the gifti reader/writer.
% GIFTI is from
% http://www.artefact.tk/software/matlab/gifti/
% https://github.com/nno/matlab_GIfTI
%
% The examples here
%
%  1. Run through the gifti team example
%  2. Prepare T1 and itkGray class file for viewing with dtiquery
%  3. Produces a gifti mesh from a gray/white segmentation class file
%  4. Illustrations of mesh computations
%       reducing the number of triangles and computing the normals
%       (isonormals, reducepatch).
%  5. Interoperability of GIFTI and mrMesh meshes
%     * Loads a VISTASOFT mesh, converts it to GIFTI format, colors it,
%   attaches metadata, saves it, and display the results in a Matlab window.
%     * Load a GIFTI mesh, converts it to a VISTASOFT mesh, displays it in a
%  mrMesh window.
%
% This script should be separated into about 3 other simpler scripts.
%
% (c) Stanford VISTA Team

%%  Set the vista data path and change into the gifti data directory

chdir(fullfile(vistaRootPath,'local'));

%% Read T1 anatomical and itkGray class file

% The T1 nifti is called the background image. That data and the two gifti
% meshes created in this cell are used within dtiquery.  The gifti files
% are assigned a transform in the the g.mat slot that aligns the T1 and
% meshes. The package (t1.nii.gz, left_gifti and right_gifti) are all used
% by dtiquery as part of the visualization.
%
% Read the anatomical.  We will use the transform in qto_xyz for
% coregistering.

% Help
%{
 % If youi need to download the nifti files, use this
 rdt.crp('/vistadata/anatomy/T1andMesh');
 niT1File = rdt.readArtifact('t1.nii',...
    'type','gz',...
    'destinationFolder',fullFolderName);
%}

niT1File = fullfile(vistaRootPath,'local','t1.nii.gz');
niT1 = niftiRead(niT1File);

% In this example, we produce the gifti data from an itkGray segmentation.
% The class file with gray identified is written by mrgSaveClassWithGray.
% That routine reads the itkGray class file, grows gray matter separately
% for left and right, and saves the output.
niCFile = fullfile(vistaRootPath,'local','t1_class_5GrayLayers.nii.gz');
niClass = niftiRead(niCFile);

%% Make an isosurface between left gray and white

% These are the ITKGRAY class labels
%    0: unlabeled
%    1: CSF
%    2: Subcortical
%    3: left white matter
%    5: left gray matter
%    4: right white matter
%    6: right gray matter

Ds = uint8(niClass.data);

% To get the boundary between white and gray in left hemisphere set the
% labels for everything that is not left gray or white to unlabeled.
Ds(Ds == 2) = 0;    % Subcortical
Ds(Ds == 4) = 0;    % right white
Ds(Ds == 6) = 0;    % right gray

%% Matlab calculations reducing the patches and computing normals

% Makes the isosurface at level 4 (20 sec)
fv = isosurface(Ds,4);

% Reduce the number of patches (10 sec)
fv = reducepatch(fv,.06);

% A blurry, red brain, no shading.  Ugly.
mrvNewGraphWin;
p1 = patch(fv, 'FaceColor','red','EdgeColor','none');

%%  Make it nicer

% Plots the surface using the normals, I think

mrvNewGraphWin;
p1 = patch(fv, 'FaceColor','red','EdgeColor','none');
fv.normals = isonormals(Ds,fv.vertices);
view(3); daspect([1 1 1]); axis tight
camlight;  camlight(-80,-10); lighting phong;
title('Data Normals')

%% End

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [yCouch, lines] =  getCouchLocationHough(scan3M,minLengthOpt,retryOpt)

% Function: getCouchLocationHough
% Description: Returns anterior coordinate of patient couch surface
%
% Usage
% inputStack is a 3-dim array of CT scan image with array axes [Y,X,Z]. [Note: "Image" notation permutes Y <-> X]
% yCouch, the Y-coordinate of the scanning table/couch
%
% EML 2020-04-13
%

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

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

midptS = floor(size(scan3M,1)/2);

maxM = max(scan3M, [], 3);
histeqM = histeq(maxM);
edgeM1 = edge(histeqM,'sobel',[],'horizontal');
edgeM2 = bwmorph(edgeM1,'thicken');

[H,T,R] = hough(edgeM2);
P = houghpeaks(H,20);

if isempty(minLengthOpt)
    minLength = floor(size(edgeM2,2)/8); % couch covers 1/8th of image
else
    minLength = minLengthOpt;
end

% lines = houghlines(edgeM2,T,R,P,'FillGap',5,'MinLength',minLength);
lines = houghlines(edgeM2,T,R,P);
overlapFraction = zeros(1,numel(lines));
midV = [floor(0.5*midptS):floor(0.5*midptS) + midptS];
% Require couch lines to have same starting & ending point2
yi = zeros(1,numel(lines));
% figure; imagesc(maxM); axis equal; hold on
for i = 1:numel(lines)
    len = norm(lines(i).point1 - lines(i).point2);
    if lines(i).point1(2) == lines(i).point2(2) && len > minLength
%         xy = [lines(i).point1; lines(i).point2];
%         plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');
%         % Plot beginnings and ends of lines
%         plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');
%         plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red');
        lineV = [lines(i).point1(1):lines(i).point2(1)];
        if lines(i).point1(2) > midptS && ~isempty(intersect(lineV,midV))
            yi(i) = lines(i).point2(2);
            overlapFraction(i) = numel(intersect(lineV,midV));
        end
    end
end

if any(overlapFraction)
    [~,I] = max(overlapFraction);
    yCouch = yi(I);
else
    yCouch = min(yi(find(yi > 0)));
end

if retryOpt && isempty(yCouch)
    [yCouch, lines] =  getCouchLocationHough(scan3M,minLength/2);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function som = SOMSimple(nfeatures, ndim, nepochs, ntrainingvectors, eta0, etadecay, sgm0, sgmdecay, showMode)
%SOMSimple Simple demonstration of a Self-Organizing Map that was proposed by Kohonen.
%   sommap = SOMSimple(nfeatures, ndim, nepochs, ntrainingvectors, eta0, neta, sgm0, nsgm, showMode)
%   trains a self-organizing map with the following parameters
%       nfeatures        - dimension size of the training feature vectors
%       ndim             - width of a square SOM map
%       nepochs          - number of epochs used for training
%       ntrainingvectors - number of training vectors that are randomly generated
%       eta0             - initial learning rate
%       etadecay         - exponential decay rate of the learning rate
%       sgm0             - initial variance of a Gaussian function that
%                          is used to determine the neighbours of the best
%                          matching unit (BMU)
%       sgmdecay         - exponential decay rate of the Gaussian variance
%       showMode         - 0: do not show output,
%                          1: show the initially randomly generated SOM map
%                             and the trained SOM map,
%                          2: show the trained SOM map after each update
%
%   For example: A demonstration of an SOM map that is trained by RGB values
%
%       som = SOMSimple(3,60,10,100,0.1,0.05,20,0.05,2);
%       % It uses:
%       %   3    : dimensions for training vectors, such as RGB values
%       %   60x60: neurons
%       %   10   : epochs
%       %   100  : training vectors
%       %   0.1  : initial learning rate
%       %   0.05 : exponential decay rate of the learning rate
%       %   20   : initial Gaussian variance
%       %   0.05 : exponential decay rate of the Gaussian variance
%       %   2    : Display the som map after every update

nrows = ndim;
ncols = ndim;

som = rand(nrows,ncols,nfeatures);

if showMode >= 1
    fig = figure;
    displaySOMmap(fig, 1, 'Randomly initialized SOM', som, nfeatures);
end

% Generate random training data
trainingData = rand(ntrainingvectors,nfeatures);

% Generate coordinate system
[x y] = meshgrid(1:ncols,1:nrows);

for t = 1:nepochs
    % Compute the learning rate for the current epoch
    eta = eta0 * exp(-t*etadecay);

    % Compute the variance of the Gaussian (Neighbourhood) function for the ucrrent epoch
    sgm = sgm0 * exp(-t*sgmdecay);

    % Consider the width of the Gaussian function as 3 sigma
    width = ceil(sgm*3);

    for ntraining = 1:ntrainingvectors
        % Get current training vector
        trainingVector = trainingData(ntraining,:);

        % Compute the Euclidean distance between the training vector and
        % each neuron in the SOM map
        dist = getEuclideanDistance(trainingVector, som, nrows, ncols, nfeatures);

        % Find the best matching unit (bmu)
        [~, bmuindex] = min(dist);

        % transform the bmu index into 2D
        [bmurow bmucol] = ind2sub([nrows ncols],bmuindex);

        % Generate a Gaussian function centered on the location of the bmu
        g = exp(-(((x - bmucol).^2) + ((y - bmurow).^2)) / (2*sgm*sgm));

        % Determine the boundary of the local neighbourhood
        fromrow = max(1,bmurow - width);
        torow   = min(bmurow + width,nrows);
        fromcol = max(1,bmucol - width);
        tocol   = min(bmucol + width,ncols);

        % Get the neighbouring neurons and determine the size of the neighbourhood
        neighbourNeurons = som(fromrow:torow,fromcol:tocol,:);
        sz = size(neighbourNeurons);

        % Transform the training vector and the Gaussian function into
        % multi-dimensional to facilitate the computation of the neuron weights update
        T = reshape(repmat(trainingVector,sz(1)*sz(2),1),sz(1),sz(2),nfeatures);
        G = repmat(g(fromrow:torow,fromcol:tocol),[1 1 nfeatures]);

        % Update the weights of the neurons that are in the neighbourhood of the bmu
        neighbourNeurons = neighbourNeurons + eta .* G .* (T - neighbourNeurons);

        % Put the new weights of the BMU neighbouring neurons back to the
        % entire SOM map
        som(fromrow:torow,fromcol:tocol,:) = neighbourNeurons;

        if showMode == 2
            displaySOMmap(fig, 2, ['Epoch: ',num2str(t),'/',num2str(nepochs),', Training Vector: ',num2str(ntraining),'/',num2str(ntrainingvectors)], som, nfeatures);
        end
    end
end

if showMode == 1
    displaySOMmap(fig, 2, 'Trained SOM', som, nfeatures);
end

function ed = getEuclideanDistance(trainingVector, sommap, nrows, ncols, nfeatures)

% Transform the 3D representation of neurons into 2D
neuronList = reshape(sommap,nrows*ncols,nfeatures);

% Initialize Euclidean Distance
ed = 0;
for n = 1:size(neuronList,2)
    ed = ed + (trainingVector(n)-neuronList(:,n)).^2;
end
ed = sqrt(ed);

function displaySOMmap(fig, nsubplot, description, sommap, nfeatures)
% Display given SOM map

figure(fig);
subplot(1,2,nsubplot);
if nfeatures >= 3
    imagesc(sommap(:,:,1:3));
else
    imagesc(sommap(:,:,1));
end
axis off;axis square;
title(description);

author: Sable

repo: mcbench-benchmarks


text:
function im=mat2im(mat,cmap,limits)
% mat2im - convert to rgb image
%
% function im=mat2im(mat,cmap,maxVal)
%
% PURPOSE
% Uses vectorized code to convert matrix "mat" to an m-by-n-by-3
% image matrix which can be handled by the Mathworks image-processing
% functions. The the image is created using a specified color-map
% and, optionally, a specified maximum value. Note that it discards
% negative values!
%
% INPUTS
% mat     - an m-by-n matrix
% cmap    - an m-by-3 color-map matrix. e.g. hot(100). If the colormap has
%           few rows (e.g. less than 20 or so) then the image will appear
%           contour-like.
% limits  - by default the image is normalised to it's max and min values
%           so as to use the full dynamic range of the
%           colormap. Alternatively, it may be normalised to between
%           limits(1) and limits(2). Nan values in limits are ignored. So
%           to clip the max alone you would do, for example, [nan, 2]
%
%
% OUTPUTS
% im - an m-by-n-by-3 image matrix
%
%
% Example 1 - combine multiple color maps on one figure
% clf, colormap jet, r=rand(40);
% subplot(1,3,1),imagesc(r), axis equal off , title('jet')
% subplot(1,3,2),imshow(mat2im(r,hot(100))) , title('hot')
% subplot(1,3,3),imshow(mat2im(r,summer(100))), title('summer')
% colormap winter %changes colormap in only the first panel
%
% Example 2 - clipping
% p=peaks(128); J=jet(100);
% subplot(2,2,1), imshow(mat2im(p,J)); title('Unclipped')
% subplot(2,2,2), imshow(mat2im(p,J,[0,nan])); title('Remove pixels <0')
% subplot(2,2,3), imshow(mat2im(p,J,[nan,0])); title('Remove pixels >0')
% subplot(2,2,4), imshow(mat2im(p,J,[-1,3])); title('Plot narrow pixel range')
%
% Rob Campbell - April 2009
%
% See Also: ind2rgb, imadjust

%Check input arguments
error(nargchk(2,3,nargin));

if ~isa(mat, 'double')
    mat = double(mat)+1;    % Switch to one based indexing
    limits = limits + 1;
end

if ~isnumeric(cmap)
    error('cmap must be a colormap, such as jet(100)')
end

%Clip if desired
L=length(cmap);
if nargin==3 && length(limits)==1
    warning('limits should be vector of length of 2. Assuming a max value was specified.')
    limits=[nan,limits];
end

if nargin==3
    minVal=limits(1);
    if isnan(minVal), minVal=min(mat(:)); end
    mat(mat<minVal)=minVal;

    maxVal=limits(2);
    if isnan(maxVal), maxVal=max(mat(:)); end
    mat(mat>maxVal)=maxVal;
else
minVal=min(mat(:));
maxVal=max(mat(:));
end

%Normalise
mat=mat-minVal;
mat=(mat/(maxVal-minVal))*(L-1);
mat=mat+1;

%convert to indecies
mat=round(mat);

%Vectorised way of making the image matrix
im=reshape(cmap(mat(:),:),[size(mat),3]);

author: geopavlakos


text:
function y=label_bibtex(bibtex,k)

%LABEL_BIBTEX labels bibliographic entries in bibtex file
%Syntax: y=label_bibtex(bibtex,k)
%Description: adds a label to each entry composed as Author:Year:[k letters of the title]
%             if k argument is not provided 25 will be used by default
%
%for example file file.bib with contents
%
%------------------------------------------------
%@article{de_Almeida:2001-Analysis-of-genomic-seque,
%   Author = {de Almeida, J. S. and Carrico, J. A. and Maretzek, A. and Noble, P. A. and Fletcher, M.},
%   Title = {Analysis of genomic sequences by Chaos Game Representation},
%   Journal = {Bioinformatics},
%   Volume = {17},
%   Number = {5},
%   Pages = {429-37},
%(...)
%    Year = {2001} }
%
%@article{Almeida:2002-Universal-sequence-map--U,
%   Author = {Almeida, J. S. and Vinga, S.},
%   Title = {Universal sequence map (USM) of arbitrary discrete sequences},
%   Journal = {BMC Bioinformatics},
%   Volume = {3},
%   Number = {1},
%
%   etc ...
%------------------------------------------------
%
%will be rewritten by y=endnote2bibtex('file.bib',25) as
%
%------------------------------------------------
%   Author = {Almeida, J. S. and Carrico, J. A. and Maretzek, A. and Noble, P. A. and Fletcher, M.},
%   Title = {Analysis of genomic sequences by Chaos Game Representation},
%   Journal = {Bioinformatics},
%   Volume = {17},
%   Number = {5},
%(...)
%    Year = {2001} }
%
%@article{Almeida:2002:Universal sequence map (U...,
%   Author = {Almeida, J. S. and Vinga, S.},
%   Title = {Universal sequence map (USM) of arbitrary discrete sequences},
%   Journal = {BMC Bioinformatics},
%   Volume = {3},
%   Number = {1},
%
%   etc ...
%------------------------------------------------
%
%Jonas Almeida, almeidaj@musc.edu, December 2004

if nargin<2;k=25;end

%1.Reads all contents and extracts Author, Year and Title
fid=fopen(bibtex,'r');
L=[];i=0;
while ~feof(fid)
    linha=fgetl(fid);i=i+1;%disp(linha);
    [a,b,c,d,e]=regexpi(linha,'@article{.*');
    if length(e)>0%this is the beggining of an entry
        L=[L,i]; %keep track of the lines that have to be changed
        j=length(L);
        %label{j}='';
    else
        [a,b,c,d,e]=regexpi(linha,'\s+Author = {(.+?),.*');
        if length(e)>0
            %disp(e{1});
            %label(j)=e{1};label{j}=[label{j},':'];
            Author{j}=regexprep(e{1}{1},'\W','_');
            %Author{j}=e{1}{1};
        else
            [a,b,c,d,e]=regexpi(linha,'\s+Year = {(\w+)}.*');
            if length(e)>0
                %disp(e{1});
                %label(j)=e{1};label{j}=[label{j},':'];
                Year{j}=e{1}{1};
            else
                [a,b,c,d,e]=regexpi(linha,'\s+Title = {(.*?)}.*');
                if length(e)>0
                    %disp(e{1});
                    %label(j)=e{1};label{j}=[label{j},':'];
                    Title{j}=regexprep(e{1}{1},'\W','-');
                end
            end
        end
    end
end
fclose(fid);
y.L=L;
y.Author=Author;
y.Year=Year;
y.Title=Title;

%2. And now writes it back to the file
%first copy to temporary file
fid1=fopen(bibtex,'r');
fid2=fopen(['temp_',bibtex],'w');
i=0;j=1;
while ~feof(fid1)
    i=i+1;linha=fgetl(fid1);
    if sum(L==i)>0
        if k>0
            fprintf(fid2,'%s\n',[linha(1:9),y.Author{j},':',y.Year{j},'-',y.Title{j}(1:min([k,length(y.Title{j})])),',']);
        else
            fprintf(fid2,'%s\n',[linha(1:9),y.Author{j},':',y.Year{j},',']);
        end
        j=j+1;
    else
        fprintf(fid2,'%s\n',linha);
    end
end
fclose(fid1);fclose(fid2);
%copy temprary file onto original
copyfile(['temp_',bibtex],bibtex);
%delete temporary file
delete(['temp_',bibtex]);


text:
function [reg_min,G,reg_param,minG] = gcv(U,s,b,method)
%GCV Plot the GCV function and find its minimum.
%
% [reg_min,G,reg_param] = gcv(U,s,b,method)
% [reg_min,G,reg_param] = gcv(U,sm,b,method)  ,  sm = [sigma,mu]
%
% Plots the GCV-function
%          || A*x - b ||^2
%    G = -------------------
%        (trace(I - A*A_I)^2
% as a function of the regularization parameter reg_param.
% Here, A_I is a matrix which produces the regularized solution.
%
% The following methods are allowed:
%    method = 'Tikh' : Tikhonov regularization   (solid line )
%    method = 'tsvd' : truncated SVD or GSVD     (o markers  )
%    method = 'dsvd' : damped SVD or GSVD        (dotted line)
% If method is not specified, 'Tikh' is default.
%
% If any output arguments are specified, then the minimum of G is
% identified and the corresponding reg. parameter reg_min is returned.

% Per Christian Hansen, IMM, Dec. 16, 2003.

% Reference: G. Wahba, "Spline Models for Observational Data",
% SIAM, 1990.

% Set defaults.
if (nargin==3), method='Tikh'; end  % Default method.
npoints = 200;                      % Number of points on the curve.
smin_ratio = 16*eps;                % Smallest regularization parameter.

% Initialization.
[m,n] = size(U); [p,ps] = size(s);
beta = U'*b; beta2 = norm(b)^2 - norm(beta)^2;
if (ps==2)
  s = s(p:-1:1,1)./s(p:-1:1,2); beta = beta(p:-1:1);
end
if (nargout > 0), find_min = 1; else find_min = 0; end

if (strncmp(method,'Tikh',4) | strncmp(method,'tikh',4))

  % Vector of regularization parameters.
  reg_param = zeros(npoints,1); G = reg_param; s2 = s.^2;
  reg_param(npoints) = max([s(p),s(1)*smin_ratio]);
  ratio = (s(1)/reg_param(npoints))^(1/(npoints-1));
  for i=npoints-1:-1:1, reg_param(i) = ratio*reg_param(i+1); end

  % Intrinsic residual.
  delta0 = 0;
  if (m > n && beta2 > 0), delta0 = beta2; end

  % Vector of GCV-function values.
  for i=1:npoints
    G(i) = gcvfun(reg_param(i),s2,beta(1:p),delta0,m-n);
  end

  if nargout == 0
      % Plot GCV function.
      loglog(reg_param,G,'-'), xlabel('\lambda'), ylabel('G(\lambda)')
      title('GCV function')
  end

  % Find minimum, if requested.
  if (find_min)
    [minG,minGi] = min(G); % Initial guess.
    reg_min = fminbnd('gcvfun',...
      reg_param(min(minGi+1,npoints)),reg_param(max(minGi-1,1)),...
      optimset('Display','off'),s2,beta(1:p),delta0,m-n); % Minimizer.
    minG = gcvfun(reg_min,s2,beta(1:p),delta0,m-n); % Minimum of GCV function.
    if nargout == 0
        ax = axis;
        HoldState = ishold; hold on;
        loglog(reg_min,minG,'*r',[reg_min,reg_min],[minG/1000,minG],':r')
        title(['GCV function, minimum at \lambda = ',num2str(reg_min)])
        axis(ax)
        if (~HoldState), hold off; end
    end
  end

elseif (strncmp(method,'tsvd',4) | strncmp(method,'tgsv',4))

  % Vector of GCV-function values.
  rho2(p-1) = abs(beta(p))^2;
  if (m > n & beta2 > 0), rho2(p-1) = rho2(p-1) + beta2; end
  for k=p-2:-1:1, rho2(k) = rho2(k+1) + abs(beta(k+1))^2; end
  G = zeros(p-1,1);
  for k=1:p-1
    G(k) = rho2(k)/(m - k + (n - p))^2;
  end
  reg_param = (1:p-1)';

  % Plot GCV function.
  semilogy(reg_param,G,'o'), xlabel('k'), ylabel('G(k)')
  title('GCV function')

  % Find minimum, if requested.
  if (find_min)
    [minG,reg_min] = min(G);
    ax = axis;
    HoldState = ishold; hold on;
    semilogy(reg_min,minG,'*r',[reg_min,reg_min],[minG/1000,minG],':r')
    title(['GCV function, minimum at k = ',num2str(reg_min)])
    axis(ax);
    if (~HoldState), hold off; end
  end

elseif (strncmp(method,'dsvd',4) | strncmp(method,'dgsv',4))

  % Vector of regularization parameters.
  reg_param = zeros(npoints,1); G = reg_param;
  reg_param(npoints) = max([s(p),s(1)*smin_ratio]);
  ratio = (s(1)/reg_param(npoints))^(1/(npoints-1));
  for i=npoints-1:-1:1, reg_param(i) = ratio*reg_param(i+1); end

  % Intrinsic residual.
  delta0 = 0;
  if (m > n & beta2 > 0), delta0 = beta2; end

  % Vector of GCV-function values.
  for i=1:npoints
    G(i) = gcvfun(reg_param(i),s,beta(1:p),delta0,m-n,1);
  end

  % Plot GCV function.
  loglog(reg_param,G,':'), xlabel('\lambda'), ylabel('G(\lambda)')
  title('GCV function')

  % Find minimum, if requested.
  if (find_min)
    [minG,minGi] = min(G); % Initial guess.
    reg_min = fminbnd('gcvfun',...
      reg_param(min(minGi+1,npoints)),reg_param(max(minGi-1,1)),...
      optimset('Display','off'),s,beta(1:p),delta0,m-n,1); % Minimizer.
    minG = gcvfun(reg_min,s,beta(1:p),delta0,m-n,1); % Minimum of GCV function.
    ax = axis;
    HoldState = ishold; hold on;
    loglog(reg_min,minG,'*r',[reg_min,reg_min],[minG/1000,minG],':r')
    title(['GCV function, minimum at \lambda = ',num2str(reg_min)])
    axis(ax)
    if (~HoldState), hold off; end
  end

elseif (strncmp(method,'mtsv',4) | strncmp(method,'ttls',4))

  error('The MTSVD and TTLS methods are not supported')

else
  error('Illegal method')
end

repo: REST


function [vals, c_index_1, c_index_2, y_1, y_2] = price_2d_ctmc( S_0s, T, r, rho, sigmas, qs, params, contractParams, M)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% About: Pricing Function for European/Bermudan/Barrier Options using CTMC approximation method
% Models Supported: 2D Diffusions
% Returns: price of contract
% Author: Justin Lars Kirkby
%
% References:  (1) A General Continuous Time Markov Chain Approximation for
%               Multi-Asset option pricing with systems of correlated diffusions,
%               Applied Math. and Comput., 2020 (JL Kirkby, Duy Nguyen, Dang Nguyen)
%
% ----------------------
% Model Params
% ----------------------
% S_0s   = initial asset prices
% T      = time remaining until maturity (in years, e.g. T=1)
% r      = interest rate (e.g. 0.05)
% rho    = correlation between brownian motions
% sigmas = volatilities per asset
% qs     = div yeilds per asset
% M      = num monitoring points for Barrier/Bermudan (also controls num steps for multi-step European pricing)
%
% ----------------------
% Contract Params  (contractParams)
% ----------------------
%
% contractParams.payoff_type:
%
% 1: Linear, G = S_1  (linear payoff in first underlying)
% 2: Linear, G = S_2  (linear payoff in second underlying)
% 3: Exchange, G = (S_1 - S_2)^+
% 4: Spread,  G = (S_1 - S_2 - K)^+   (NOTE: must set strike, K)
% 5: Geometric Basket Call / Put,  G = (sqrt(S_1) * sqrt(S_2) - K)^+  (for the call)
% 6: Arithmetic Basket Call / Put,  G = (sqrt(S_1) * sqrt(S_2) - K)^+  (for the call)
% 7: Call-on-Max and Put-on-Min, Gcall = (max(S_1,S_2) - K)^+ , Gput = (K - min(S_1,S_2))^+
% 8: Call/put on just S_2, G = (S_2 - K)^+  (for the call)
% 9: Best-of / Worst-of,  G = max(S_1,S_2), G = min(S_1,S_2)
%
% contractParams.contract:
%
% 1: European, Single Step Pricing
% 2: European, Multi Step Pricing (M above controls number of steps)
% 3: Bermudan (M above controls number of monitoring points)
% 4: Barrier (M above controls number of monitoring points)
%
% Note: for barrier option, set:
%       contractParams.barriers_1 = lower/upper barriers on first asset  (e.g. [0 50])
%       contractParams.barriers_2 = lower/upper barriers on second asset  (e.g. [0 50000000000] to disable barrier on S_2)
%
% ----------------------
% Numerical (CTMC) Params
% ----------------------
% params = CTMC parameters
% params.m_0 = num CTCM states
% params.num_devs = num std devs used in the grid
% params.gridMethod = choose the grid method (several RnD versions)
% params.GridMultParam = non-uniformity param, in (0,1)

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

if nargin < 9
    M = 1; % M is only needed for
end
dt = T/M;

contract = contractParams.contract;

if contract == 1  % European
    dt = 1; M = 1;
end

method = 4;
num_devs = params.num_devs;
m_0 = params.m_0;
GridMultParam = params.GridMultParam;
gridMethod = params.gridMethod;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

drifts = r - qs;

R = [1 rho; rho 1];  % Correlation Matrix
[ L, D, C, Cinv ] = get_transform_matrices_2d( R, method );

% Now Define New Uncorrelated System  (the dc underscore)
[drift_dc, sigma_dc ] = decorrelate(sigmas, drifts, C, D );

[Ls_dc, Rs_dc ] = get_CTMC_decorr_boundaries(sigmas, C, T, 0, sigma_dc, num_devs);
Y_0s = [0 0];

% Form CTMC 1
center = Y_0s(1);
mu_func = @(s) drift_dc(1)*[s>-100000];
sig_func = @(s) sigma_dc(1)*[s>-100000];
[Q, y_1, c_index_1] = Q_Matrix(m_0, mu_func,sig_func,Ls_dc(1),Rs_dc(1),gridMethod,center, GridMultParam);
P1 = expm(Q*dt);

% Form CTMC 2
center = Y_0s(2);
mu_func = @(s) drift_dc(2)*[s>-100000];
sig_func = @(s) sigma_dc(2)*[s>-100000];
[Q, y_2, c_index_2] = Q_Matrix(m_0, mu_func,sig_func,Ls_dc(2),Rs_dc(2),gridMethod,center, GridMultParam);
P2 = expm(Q*dt);

G = get_payoff_G_matrix_from_ygrid_2d( y_1, y_2, S_0s, sigmas, rho, contractParams);

if contract == 1  % European, Price by single step
    vals = exp(-r*T)*P1*G*P2.';

elseif contract == 2  % European, Price By Multi Step
    vals = G;
    for m=M-1:-1:0
        vals = exp(-r*dt)*P1*vals*P2.';
    end

elseif contract == 3  % Bermudan
    vals = G;
    for m=M-1:-1:0
        vals = max(exp(-r*dt)*P1*vals*P2.', G);
    end

elseif contract == 4 % Barrier
    b1 = contractParams.barriers_1; L1 = b1(1); U1 = b1(2);
    b2 = contractParams.barriers_2; L2 = b2(1); U2 = b2(2);
    B = ones(m_0, m_0);
    for i = 1:m_0
        y1 = y_1(i);
        S1 = S_0s(1)*exp(sigmas(1)*y1);
        if S1 < L1 || S1 > U1
            B(i,:) = 0;
        else
            for j = 1:m_0
                y2 = y_2(j);
                S2 = S_0s(2)*exp(sigmas(2)*(y2 + rho*y1));
                if S2 < L2 || S2 > U2
                    B(i,j) = 0;
                end
            end
        end
    end
    vals = G.*B;
    for m=M-1:-1:0
        vals = exp(-r*dt)*(B.*P1*vals*P2.');
    end
end

end


text:
%This Progra, develop by Renoald Tang
%University Teknologi Malaysia, Photogrammetry and Laser scanning group
%for academic purpose
%Email:renoald@live.com
function varargout = AMLS(varargin)
% AMLS M-file for AMLS.fig
%      AMLS, by itself, creates a new AMLS or raises the existing
%      singleton*.
%
%      H = AMLS returns the handle to a new AMLS or the handle to
%      the existing singleton*.
%
%      AMLS('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in AMLS.M with the given input arguments.
%
%      AMLS('Property','Value',...) creates a new AMLS or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before AMLS_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to AMLS_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

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

% Last Modified by GUIDE v2.5 01-Jul-2004 00:11:31

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

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

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

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

% Update handles structure
guidata(hObject, handles);
%axes(handles.axes1)
%axis off
%axes(handles.axes3)
%axis off

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

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

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

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

[FileName,PathName] = uigetfile({'*.xyz'},'Select the Point clouds');
%combine the filename and pathname
directory=strcat(PathName ,FileName);
%.xyz check file extension
ext=findstr(FileName,'.xyz');
n=isempty(ext);
if n==0
    point=dlmread(directory);
end
n=length(point);
    dlmwrite('point.xyz', point, 'delimiter', '\t','precision', 6);
 %setappdata(AMLS,'point1',point);
 axes(handles.axes1)
 axis on
plot3(point(:,1),point(:,2),point(:,3),'r.')
grid on
num= num2str(n);
set(handles.text2,'String',num);

function edit1_Callback(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit1 as text
%        str2double(get(hObject,'String')) returns contents of edit1 as a double
input = str2num(get(hObject,'String'));
if (isempty(input))
     set(hObject,'String','2')
end
guidata(hObject, handles);

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

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in Computebutton.
function Computebutton_Callback(hObject, eventdata, handles)
% hObject    handle to Computebutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
para = get(handles.edit1,'String');
parameter=str2num(para);
myname='point.xyz';
out='mypoint';
command=sprintf('NormFet.exe %s  %s -nml -ftr \n',myname,out);
system(command);
lastfile='mypoint_bbr2.50_nn8.xyznw';
outfile='output';
factor=parameter;
command2=sprintf('AMLS.exe %s %s -hr %d\n',lastfile,outfile,factor);
system(command2);
f1='output_adpmls_hr';
f2='.00.xyz';
finalfile=strcat(f1,para,f2);
%disp(finalfile)
surface=dlmread(finalfile);
setappdata(AMLS,'surface',surface);
axes(handles.axes3)
axis on
plot3(surface(:,1),surface(:,2),surface(:,3),'g.');
grid on
% --- Executes on button press in savebutton.
function savebutton_Callback(hObject, eventdata, handles)
% hObject    handle to savebutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[file1,path1] = uiputfile({'*.xyz'},'Save file name');
%combine the filename and pathname

directory=strcat(path1,file1);
mypoint=getappdata(AMLS,'surface');
dlmwrite(directory, mypoint, 'delimiter', '\t', ...
         'precision', 6);

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

% --- Executes on button press in clearpushbutton.
function clearpushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to clearpushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
arrayfun(@cla,findall(0,'type','axes'));


text:
% J. C. Spall, March 1999
% This code evaluates the basic root-finding (R-M) SA algorithm.  Computes sample mean
% and standard deviation of the mean from Navg realizations.
%
n=1000;
randn('seed',71111113)
a=0.081;
A=5.1;
alpha=.501;
s=.1;
theta_0=[1,1]';
Navg=2000;
loss_eq2_5(theta_0)
rslossavg=0;
rslossavgsq=0;
for i=1:Navg
  theta=theta_0;
  for k=1:n
    theta=theta-(a/(k+A)^alpha)*gradloss(theta,s)';
  end
  rslossavg=(i-1)*rslossavg/i+loss_eq2_5(theta)/i;
  rslossavgsq=(i-1)*rslossavgsq/i+(loss_eq2_5(theta)^2)/i;
end
mean=rslossavg
stdevmean=(((Navg/(Navg-1))*(rslossavgsq-mean^2))^.5)/(Navg^.5)


text:
function Constraint = dd(X)

if issymmetric(X)
    W = X-diag(diag(X));
    Z = sdpvar(length(X));
    Constraint = [-Z(:) <= W(:) <= Z(:), diag(X) >= sum(Z,2)];
else
    error('dd requires a symmetric argument.');
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function cvx_clearpath( arg ) %#ok

%CVX_CLEARPATH   Clears the cvx path.
%   CVX_CLEARPATH removes the internal cvx directories from Matlab's path. CVX
%   does this automatically when a model is completed (i.e., after CVX_END), in
%   order to reduce potential naming conflicts with other packages. There is no
%   need to call this function during the normal use of CVX.

global cvx___
cvx_global
if nargin == 0,
    cvx___.path.hold = false;
end
if cvx___.path.hold,
    cvx_setspath;
else
    cvx_clearspath;
    if ~isempty( cvx___.path.string ),
        path(strrep(path,cvx___.path.string,''));
    end
    cvx___.path.active = false;
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function fmri_loc1_GLMcontrasts

% LO Localizer scans processing code
%
% USER IMPORTANT!!!: Before running this script, make sure that all loloc
% scans have the correct annotation fields and scan groups (are grouped
% correctly). Otherwise, fmri_FFA_checkDataTypesForGLM will fail.
%
% 1. Set GLM parameters for all scans
% 2. Check LO localizer experiments for:
%    a. Datatype and scan numbers
%    b. Record parfiles
%    c. Make sure scan groupings are correct
% 3. Run GLM for LO localizer experiments
% 4. Run contrasts for the LO localizer experiments
%
% Calls the functions:
% * fmri_loc1_setGLMparams (subfunction)
% * fmri_loc1_checkDataTypesForGLM
% * fmri_loc1_computeContrastsLocalizer
%
%
% DY 06/05/2008
% Modified 06/18/2008: there are calls to er_setParams and saveSessions all
% over the place, which results in the params for the MotionComp_RefScan1
% being changed for the worse. I'm not sure how to efficiently prevent this
% from happening, so instead, I let it happen, and just fix the ANNOTATION
% and PARFILE fields at the end, since accuracy here is important.
%
% Code taken from fmri_FFA_motionCompFix
% AL 11/08: Adapted for loc1 GLM

clear all

% Set directory and subject list
% Alternatively: subs = {'sess1','sess2','sess3',...};
if ispc
    fmriDir = 'W:\projects\Kids\fmri\localizer\';
else
    %fmriDir = '/biac1/kgs/projects/Kids/fmri/localizer/';
    fmriDir='/biac1/kgs/projects/Kids/fmri/code/Alina/Subject/';
end
cd(fmriDir); s=dir('*0*'); subs={s.name};

% Start a log text file to document successes and failures in preprocessing
dateAndTime=datestr(now); dateAndTime(12)='_'; dateAndTime(15)='h';dateAndTime(18)='m';
logDir=fullfile(fmriDir,'logs');
if ~isdir(logDir), mkdir(logDir), end;
logFile = fullfile(fmriDir,'logs',['LO1_Log_' dateAndTime '.txt']);
fid=fopen(logFile,'w');
[junk,logFileName]=fileparts(logFile)
scriptName = 'fmri_loc1_GLMcontrasts.m';
startTime = clock;

% 1. Set GLM parameters
params=fmri_loc1_setGLMparams(fid,logFileName,scriptName); % Subfunction within this script

for ii=1:length(subs)

    % Get rid of all mrVista globals in workspace. Put this at the top,
    % just in case we break from the rest of the loop.
    clear dataTYPES mrSESSION vANATOMYPATH HOMEDIR hI newDt newScan

    thisDir = fullfile(fmriDir,subs{ii}); cd(thisDir);
    fprintf(fid,'\n ------------------------------------------ \n');
    fprintf(fid,'Processing %s\n\n',thisDir);
    fprintf('Processing %s\n\n',thisDir);

    %% Check if the subject already has been processed by script (and thus
    %% should not be re-entered to this script). If they do, skip to the
    %% next subject by using CONTINUE.
    load mrSESSION.mat
    theglm=find(strcmp({dataTYPES.name},'GLMs'));
    if ~isempty(theglm)
        checkScans=length(dataTYPES(theglm).eventAnalysisParams);
        alreadyDone=[];
        for xx=1:checkScans
            tmp=(strfind(lower(dataTYPES(theglm).eventAnalysisParams(xx).annotation),'lo localizer run 1'));
            alreadyDone=[alreadyDone tmp];
        end
        if(~isempty(alreadyDone))
            fprintf(fid,'Error: Subject already processed (has GLM with ''LO localizer Run 1'' in annotation field), skipping... \n\n');
            fprintf('Error: Subject already processed (has GLM with ''LO localizer Run 1'' in annotation field), skipping... \n\n');
            continue
        end
    end

    %% 2a/b/c. Find datatype number for MotionComp_RefScan1, find scan
    %% numbers for lo localizer scans, record parfiles for these scans and
    %% make sure scan groupings are correct. If everything checks out,
    %% proceed with GLM. If not, SKIP THIS SUBJECT.
    [go,s1]=fmri_loc1_checkDataTypesForGLM(fid,dataTYPES);
    %clear dataTYPES mrSESSION vANATOMYPATH

    % If there was a problem in step 3, skip GLM and proceed to next
    % subject
    if go==0
        continue

    %% 3. Otherwise, run GLM (to debug: fid=1)
    elseif go==1
        try
            % Initialize a hidden view of the 'inplane' data (like invoking
            % mrVista without a GUI, dataTYPES/etc in workspace)
            hI = initHiddenInplane(s1.mc, s1.scan);
            params.parfiles={s1.parfile};
            er_setParams(hI, params, s1.scan, s1.mc); % write params as eventAnalysisParams fields
            tic; [hI, lo_scan] = applyGlm(hI,s1.mc,s1.scangroup,params);
            time=toc;
            fprintf(fid,'GLM COMPLETED -- dt %d, scan %d. Time: %d m %2.2f s \n\n',...
                hI.curDataType,lo_scan,floor(time/60),mod(time,60));
            fprintf('GLM COMPLETED -- dt %d, scan %d. Time: %d m %2.2f s \n\n',...
                hI.curDataType,lo_scan,floor(time/60),mod(time,60));
            glmOK=1;
        catch
            fprintf(fid,'FAILURE TO SUCCESSFULLY RUN GLM \n');
            fprintf('FAILURE TO SUCCESSFULLY RUN GLM \n');
            theerror=lasterror;
            fprintf(fid,'%s \n\n',theerror.message);
            fprintf('%s \n\n',theerror.message);
            glmOK=0;
        end

        %% 4. Compute contrast maps
        if glmOK

            try
                fmri_loc1_computeContrastsLocalizer(hI, lo_scan)
                fprintf(fid,'Computed contrasts using fmri_FFA_computeContrastsLocalizer \n');
                fprintf('Computed contrasts using fmri_FFA_computeContrastsLocalizer \n');
            catch
                fprintf(fid,'FAILURE TO SUCCESSFULLY COMPUTE CONTRASTS \n');
                fprintf('FAILURE TO SUCCESSFULLY COMPUTE CONTRASTS \n');
                theerror=lasterror;
                fprintf(fid,'%s \n\n',theerror.message);
                fprintf('%s \n\n',theerror.message);

            end
        end
    end

    % Fix the annotation and parfile fields for the MotionComp_RefScan
    % dataTYPE, since these are undesirably changed (collateral damage).
    dataTYPES(s1.mc).eventAnalysisParams(s1.scan).annotation='loloc_run1';

    par1=dataTYPES(s1.mc).scanParams(s1.scan).parfile;

    dataTYPES(s1.mc).eventAnalysisParams(s1.scan).parfiles=par1;
    save(fullfile(thisDir,'mrSESSION.mat'),'mrSESSION','dataTYPES', '-append');
    clear par1

end

totalTime=etime(clock,startTime);

fprintf(fid,'\n ------------------------------------------ \n');
fprintf(fid,'Total running time for script: %f minutes \n',totalTime/60);
fprintf('Total running time for script: %f minutes \n',totalTime/60);

fclose(fid); % close out the log file

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function params=fmri_loc1_setGLMparams(fid,logFileName,scriptName)

% Set GLM parameters: We set the PARAMS struct fields to the values on the
% wiki, which are further explained there. Also we print the values to the
% log. http://vpnl.stanford.edu/internal/wiki/index.php/ROIanalysis#GLMs

fprintf(fid,'\n ------------------------------------------ \n');
fprintf(fid,'\n GLM PARAMETERS \n');
fprintf(fid,' ------------------------------------------ \n');

params=er_defaultParams % Initialize param struct with defaults
p.al='params.alpha=0.0500'; eval(p.al); fprintf(fid,'%s \n',p.al);
p.at='params.ampType=''betas'''; eval(p.at); fprintf(fid,'%s \n',p.at);
p.an='params.annotation=''LO localizer Run 1'''; eval(p.an); fprintf(fid,'%s \n',p.an);
p.ap='params.assignParfiles=0'; eval(p.ap); fprintf(fid,'%s \n',p.ap);
p.bp='params.bslPeriod=[-6 -5 -4 -3 -2 -1 0]'; eval(p.bp); fprintf(fid,'%s \n',p.bp);
p.dt='params.detrend=1'; eval(p.dt); fprintf(fid,'%s \n',p.dt);
p.df='params.detrendFrames=20;'; eval(p.df); fprintf(fid,'%s \n',p.df);
p.ea='params.eventAnalysis=1'; eval(p.ea); fprintf(fid,'%s \n',p.ea);
p.eb='params.eventsPerBlock=6'; eval(p.eb); fprintf(fid,'%s \n',p.eb);
p.fp='params.framePeriod=2'; eval(p.fp); fprintf(fid,'%s \n',p.fp);
p.gt='params.glmHRF=3'; eval(p.gt); fprintf(fid,'%s \n',p.gt);
p.gp='params.glmHRF_params=[6 16 1 1 6 0 24]'; eval(p.gp); fprintf(fid,'%s \n',p.gp);
p.gw='params.glmWhiten=0'; eval(p.gw); fprintf(fid,'%s \n',p.gw); % temporal autocorrelation
p.ic='params.inhomoCorrect=1'; eval(p.ic); fprintf(fid,'%s \n',p.ic);
p.lp='params.lowPassFilter=0'; eval(p.lp); fprintf(fid,'%s \n',p.lp);
p.nb='params.normBsl=1'; eval(p.nb); fprintf(fid,'%s \n',p.nb);
p.od='params.onsetDelta=0'; eval(p.od); fprintf(fid,'%s \n',p.od);
p.pf='params.parfiles='''''; eval(p.pf); fprintf(fid,'%s \n',p.pf); % fill in after checkDataTypes...
p.pp='params.peakPeriod=[4 5 6 7 8 9 10 11 12 13 14]'; eval(p.pp); fprintf(fid,'%s \n',p.pp);
p.sh='params.setHRFParams=0'; eval(p.sh); fprintf(fid,'%s \n',p.sh);
p.sc='params.snrConds=1'; eval(p.sc); fprintf(fid,'%s \n',p.sc);
p.tn='params.temporalNormalization=0'; eval(p.tn); fprintf(fid,'%s \n',p.tn);
p.tw='params.timeWindow=[-6:1:24]'; eval(p.tw); fprintf(fid,'%s \n',p.tw);
p.log=['params.logFile=''' logFileName '''']; eval(p.log); fprintf(fid,'%s \n',p.log)
p.script=['params.scriptName=''' scriptName '''']; eval(p.script); fprintf(fid,'%s \n',p.script)

% params.framePeriod = dataTYPES(dt).scanParams(scan).framePeriod;


text:
function [g,a] = wpfbt2filterbank( wt, varargin)
%WPFBT2FILTERBANK  WPFBT equivalent non-iterated filterbank
%   Usage: [g,a] = wpfbt2filterbank(wt)
%
%   Input parameters:
%         wt : Wavelet filter tree definition
%
%   Output parameters:
%         g   : Cell array containing filters
%         a   : Vector of sub-/upsampling factors
%
%   `wpfbt2filterbank(wt)` calculates the impulse responses *g* and the
%   subsampling factors *a* of non-iterated filterbank, which is equivalent
%   to the wavelet packet filterbank tree described by *wt*. The returned
%   parameters can be used directly in |filterbank|, |ufilterbank| or
%   |filterbank|.
%
%   Please see help on |wfbt| for description of *wt*. The function
%   additionally support the following flags:
%
%   `'freq'`(default),`'nat'`
%      The filters are ordered to produce subbands in the same order as
%      |wpfbt| with the same flag.
%
%   `'intsqrt'`(default),`'intnoscale'`, `'intscale'`
%      The filters in the filterbank tree are scaled to reflect the
%      behavior of |wpfbt| and |iwpfbt| with the same flags.
%
%   `'scaling_notset'`(default),`'noscale'`,`'scale'`,`'sqrt'`
%     Support for scaling flags as described in |uwpfbt|. By default,
%     the returned filterbank *g* and *a* is equivalent to |wpfbt|,
%     passing any of the non-default flags results in a filterbank
%     equivalent to |uwpfbt| i.e. scaled and with `a(:)=1`.
%
%   Examples:
%   ---------
%
%   The following two examples create a multirate identity filterbank
%   using a tree of depth 3. In the first example, the filterbank is
%   identical to the DWT tree:::
%
%     [g,a] = wpfbt2filterbank({'db10',3,'dwt'});
%     filterbankfreqz(g,a,1024,'plot','linabs','posfreq');
%
%
%   In the second example, the filterbank is identical to the full
%   wavelet tree:::
%
%     [g,a] = wpfbt2filterbank({'db10',3,'full'});
%     filterbankfreqz(g,a,1024,'plot','linabs','posfreq');
%
%   See also: wfbtinit

% AUTHOR: Zdenek Prusa

complainif_notenoughargs(nargin,1,'WPFBT2FILTERBANK');

definput.import = {'wfbtcommon','uwfbtcommon'};
definput.importdefaults = {'scaling_notset'};
definput.flags.interscaling={'intsqrt','intnoscale','intscale'};
[flags]=ltfatarghelper({},definput,varargin);

% build the tree
wt = wfbtinit({'strict',wt},flags.forder);

wt = comp_wpfbtscale(wt,flags.interscaling);

nIdx = nodesLevelsBForder(wt);
% Now we need to walk the tree by levels
g = {};
a = [];
for ii=1:numel(nIdx)
    rangeLoc = cellfun(@(eEl) 1:numel(eEl.h),wt.nodes(nIdx{ii}),...
                       'UniformOutput',0);
    rangeOut = cellfun(@(eEl) numel(eEl.h),wt.nodes(nIdx{ii}));
    rangeOut = mat2cell(1:sum(rangeOut),1,rangeOut);
    [gtmp,atmp] = nodesMultid(nIdx{ii},rangeLoc,rangeOut,wt);
    g(end+1:end+numel(gtmp)) = gtmp;
    a(end+1:end+numel(atmp)) = atmp;
end
g = g(:);
a = a(:);

if ~flags.do_scaling_notset
   g = comp_filterbankscale(g,a,flags.scaling);
   a = ones(numel(g),1);
end

function nodesIdxs = nodesLevelsBForder(treeStruct)

%find root
nodeNo = find(treeStruct.parents==0);
toGoTrough = [nodeNo];
nodesIdxs = {nodeNo};
inLevel = [1];
counter = 0;
level = 2;
chIdxSum = 0;
while ~isempty(toGoTrough)
   chtmp = find(treeStruct.children{toGoTrough(1)}~=0);
   chIdxtmp = treeStruct.children{toGoTrough(1)}(chtmp);
   counter = counter + 1;

   if(length(nodesIdxs)<level&&~isempty(chIdxtmp))
       nodesIdxs = {nodesIdxs{:},[]};
   end

   chIdxSum = chIdxSum + length(chIdxtmp);
   if(~isempty(chIdxtmp))
       nodesIdxs{level} = [nodesIdxs{level},chIdxtmp];
   end

   toGoTrough = [toGoTrough(2:end),chIdxtmp];

   if(counter==inLevel(level-1))
       counter = 0;
       inLevel(level) = chIdxSum;
       level = level + 1;
       chIdxSum = 0;
   end
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function x=t2x(T,str)

% x=t2x(T,str);
%
% Converts transformation matrix T between B and A coordinate
% frames into a generalized position vector x, which contains
% position and orientation vectors of B with respect to A.
% Orientation can be expressed with quaternions, euler angles
% (xyz or zxz convention), unit vector and rotation angle.
% Also both orientation and position can be expressed with
% Denavitt-Hartemberg parameters.
%
% ---------------------------------------------------------------------------
%
% The transformation matrix T between B and A coordinate
% frames is a 4 by 4 matrix such that:
% T(1:3,1:3) = Orientation matrix between B and A = unit vectors
%              of x,y,z axes of B expressed in the A coordinates.
% T(1:3,4)   = Origin of B expressed in A coordinates.
% T(4,1:3)   = zeros(1,3)
% T(4,4)     = 1
%
% ---------------------------------------------------------------------------
%
% The generalized position vector x contains the origin of B
% expressed in the A coordinates in the first four entries,
% and orientation of B with respect to A in the last four entries.
% In more detail, its shape depends on the value of str as
% specified below :
%
% ---------------------------------------------------------------------------
%
% str='van' : UNIT VECTOR AND ROTATION ANGLE
%
%           [ Ox ]     origin of the B coordinate frame
%  x(1:4) = [ Oy ]     with respect to A.
%           [ Oz ]
%           [  1 ]
%
%           [ Vx ]     Vx,Vy,Vz = unit vector respect to A,
%  x(5:8) = [ Vy ]     which B is rotated about.
%           [ Vz ]
%           [ Th ]     Th = angle which B is rotated (-pi,pi].
% ---------------------------------------------------------------------------
%
% str='qua' : UNIT QUATERNION
%
%           [ Ox ]     origin of the B coordinate frame
%  x(1:4) = [ Oy ]     with respect to A.
%           [ Oz ]
%           [  1 ]
%
%           [ q1 ]     q1,q2,q3 = V*sin(Th/2)
%  x(5:8) = [ q2 ]     q0 = cos(Th/2) where :
%           [ q3 ]     V = unit vector respect to A, which B is
%           [ q0 ]     rotated about, Th = angle which B is rotated (-pi,pi].
% ---------------------------------------------------------------------------
%
% str='erp' : EULER-RODRIGUEZ PARAMETERS
%
%           [ Ox ]     origin of the B coordinate frame
%  x(1:4) = [ Oy ]     with respect to A.
%           [ Oz ]
%           [  1 ]
%
%           [ r1 ]     r1,r2,r3 = V*tan(Th/2), where :
%  x(5:8) = [ r2 ]     V = unit vector with respect to A, which B is
%           [ r3 ]     rotated about.
%           [  0 ]     Th = angle which B is rotated (-pi,pi) (<> pi).
% ---------------------------------------------------------------------------
%
% str='rpy' : ROLL, PITCH, YAW ANGLES (euler x-y-z convention)
%
%           [ Ox ]     origin of the B coordinate frame
%  x(1:4) = [ Oy ]     with respect to A.
%           [ Oz ]
%           [  1 ]
%
%            [ r ]     r = roll angle  ( fi    (-pi,pi], about x,          )
%  x(5:8) =  [ p ]     p = pitch angle ( theta (-pi,pi], about y, <> +-pi/2)
%            [ y ]     y = yaw angle   ( psi   (-pi,pi], about z,          )
%            [ 0 ]
% ---------------------------------------------------------------------------
%
% str='rpm' : ROTATION, PRECESSION, MUTATION ANGLES (euler z-x-z convention)
%
%           [ Ox ]     origin of the B coordinate frame
%  x(1:4) = [ Oy ]     with respect to A.
%           [ Oz ]
%           [  1 ]
%
%            [ r ]     r = rotation angle   ( (-pi,pi] ,about z           )
%  x(5:8) =  [ p ]     p = precession angle ( (-pi,pi] ,about x , <> 0,pi )
%            [ y ]     y = mutation angle   ( (-pi,pi] ,about z           )
%            [ 0 ]
% ---------------------------------------------------------------------------
%
% str='dht' : DENAVITT-HARTEMBERG PARAMETERS
%
%            [ b ]                [ a ]     this four-parameter
%  x(1:4) =  [ d ] ,   x(5:8) =   [ t ] ,   description does not involve
%            [ 0 ]                [ 0 ]     a loss of information if and
%            [ 0 ]                [ 0 ]     only if T has this shape:
%
%           [    ct   -st     0     b ]     where :
%  T =      [ ca*st ca*ct   -sa -d*sa ]
%           [ sa*st sa*ct    ca  d*ca ]     sa = sin(a), ca = cos(a)
%           [     0     0     0     1 ]     st = sin(t), ct = cos(t)
% ---------------------------------------------------------------------------
%
% Example (see also x2t):
% x=[rand(3,1);1;rand(3,1);0];x-t2x(x2t(x,'rpm'),'rpm')
%
%
% Giampiero Campa 1/11/96
%

% ---------------------------------------------------------------------------
% UNIT VECTOR AND ROTATION ANGLE

if [ str=='van' size(T)==[4 4] ],

O=T(1:3,4);
R=T(1:3,1:3);
d=round(.5*(trace(R)-1)*1e12)/1e12;

if d==1,
   v=[0 0 1]';
   th=0;

elseif d==-1,
    v0=sum(R'+eye(3,3))';

    if v0 == 0
	 v0=R(:,2)+R(:,3)+[0 1 1]';
    end;

    if v0 == 0
	 v0=R(:,3)+[0 0 1]';
    end;

    v=v0/norm(v0);
    th=pi;

else

   sg=(vp([1 0 0]',R(:,1))+vp([0 1 0]',R(:,2))+vp([0 0 1]',R(:,3)));

    if norm(sg) < 1e-12
        disp(' ');
        disp('T2x warning: det(R)<>1, unit vector assumed to be [0 0 1]''.');
        disp(' ');
        sg=[0 0 1]';
    end

   v=sg/norm(sg);
   th=atan2(norm(sg)/2,d);

end

x=[O;1;v;th];

% ---------------------------------------------------------------------------
% UNIT QUATERNION

elseif [ str=='qua' size(T)==[4 4] ],

O=T(1:3,4);
R=T(1:3,1:3);
d=round(.5*(trace(R)-1)*1e12)/1e12;

if d==1,
   v=[0 0 1]';
   th=0;

elseif d==-1,
    v0=sum(R'+eye(3,3))';

    if v0 == 0
	v0=R(:,2)+R(:,3)+[0 1 1]';
    end;

    if v0 == 0
	v0=R(:,3)+[0 0 1]';
    end;

    v=v0/norm(v0);
    th=pi;

else

   sg=(vp([1 0 0]',R(:,1))+vp([0 1 0]',R(:,2))+vp([0 0 1]',R(:,3)));

    if norm(sg) < 1e-12
        disp(' ');
        disp('T2x warning: det(R)<>1, unit vector assumed to be [0 0 1]''.');
        disp(' ');
        sg=[0 0 1]';
    end

   v=sg/norm(sg);
   th=atan2(norm(sg)/2,d);

end

q=v*sin(th/2);
q0=cos(th/2);

x=[O;1;q;q0];

% ---------------------------------------------------------------------------
% EULER-RODRIGUEZ PARAMETERS

elseif [ str=='erp' size(T)==[4 4] ],

O=T(1:3,4);
R=T(1:3,1:3);
d=round(.5*(trace(R)-1)*1e12)/1e12;

if d==1,
   v=[0 0 1]';
   th=0;

elseif d==-1,
    v0=sum(R'+eye(3,3))';

    if v0 == 0
	v0=R(:,2)+R(:,3)+[0 1 1]';
    end;

    if v0 == 0
	v0=R(:,3)+[0 0 1]';
    end;

    v=v0/norm(v0);
    th=pi;

else

   sg=(vp([1 0 0]',R(:,1))+vp([0 1 0]',R(:,2))+vp([0 0 1]',R(:,3)));

    if norm(sg) < 1e-12
        disp(' ');
        disp('T2x warning: det(R)<>1, unit vector assumed to be [0 0 1]''.');
        disp(' ');
        sg=[0 0 1]';
    end

   v=sg/norm(sg);
   th=atan2(norm(sg)/2,d);

end

p=v*tan(th/2);
x=[O;1;p;0];

% ---------------------------------------------------------------------------
% ROLL, PITCH, YAW ANGLES (euler x-y-z convention)

elseif [ str=='rpy' size(T)==[4 4] ],

O=T(1:3,4);
R=T(1:3,1:3);
d=round([0 0 1]*R(:,1)*1e12)/1e12;

if d==1,
   y=atan2([0 1 0]*R(:,2),[1 0 0]*R(:,2));
   p=-pi/2;
   r=-pi/2;

elseif d==-1
   y=atan2([0 1 0]*R(:,2),[1 0 0]*R(:,2));
   p=pi/2;
   r=pi/2;

else
   sg=vp([0 0 1]',R(:,1));
   j2=sg/sqrt(sg'*sg);
   k2=vp(R(:,1),j2);

   r=atan2(k2'*R(:,2),j2'*R(:,2));
   p=atan2(-[0 0 1]*R(:,1),[0 0 1]*k2);
   y=atan2(-[1 0 0]*j2,[0 1 0]*j2);
end

y1=y+(1-sign(y)-sign(y)^2)*pi;
p1=p+(1-sign(p)-sign(p)^2)*pi;
r1=r+(1-sign(r)-sign(r)^2)*pi;

% takes smaller values of angles

if norm([y1 p1 r1]) < norm([y p r])
    x=[O;1;r1;-p1;y1;0];
else
    x=[O;1;r;p;y;0];
end

% ---------------------------------------------------------------------------
% ROTATION, PRECESSION, MUTATION ANGLES (euler z-x-z convention)

elseif [ str=='rpm' size(T)==[4 4] ],

O=T(1:3,4);
R=T(1:3,1:3);
d=round([0 0 1]*R(:,3)*1e12)/1e12;

if d==1,
   m=0;
   p=0;
   r=atan2([0 1 0]*R(:,1),[1 0 0]*R(:,1));

elseif d==-1
   m=0;
   p=pi;
   r=atan2([0 1 0]*R(:,1),[1 0 0]*R(:,1));

else
   sg=vp([0 0 1]',R(:,3));
   i2=sg/norm(sg);
   j2=vp(R(:,3),i2);

   m=atan2(j2'*R(:,1),i2'*R(:,1));
   p=atan2([0 0 1]*j2,[0 0 1]*R(:,3));
   r=atan2([0 1 0]*i2,[1 0 0]*i2);

end

r1=r+(1-sign(r)-sign(r)^2)*pi;
p1=p;
m1=m+(1-sign(m)-sign(m)^2)*pi;

% takes the smaller values of angles

if norm([r1 p1 m1]) < norm([r p m])
    x=[O;1;r1;-p1;m1;0];
else
    x=[O;1;r;p;m;0];
end

% ---------------------------------------------------------------------------
% DENAVITT-HARTEMBERG PARAMETERS

elseif [ str=='dht' size(T)==[4 4] ],

% b calculation
b=T(1,4);

% d calculation
if norm(T(3,3)) > norm(T(2,3))
    d=T(3,4)/T(3,3);
else
    d=T(2,4)/T(2,3);
end

% alfa and theta computation by mean of euler-zxz angles
R=T(1:3,1:3);
dl=round([0 0 1]*R(:,3)*1e12)/1e12;

if dl==1,
   th=0;
   alfa=0;
   n=atan2([0 1 0]*R(:,1),[1 0 0]*R(:,1));
elseif dl==-1
   th=0;
   alfa=pi;
   n=atan2([0 1 0]*R(:,1),[1 0 0]*R(:,1));
else
   sg=vp([0 0 1]',R(:,3));
   i2=sg/norm(sg);
   j2=vp(R(:,3),i2);

   th=atan2(j2'*R(:,1),i2'*R(:,1));
   alfa=atan2([0 0 1]*j2,[0 0 1]*R(:,3));
   n=atan2([0 1 0]*i2,[1 0 0]*i2);
end

th1=th+(1-sign(th)-sign(th)^2)*pi;
alfa1=alfa;
n1=n+(1-sign(n)-sign(n)^2)*pi;

% takes the smaller values of angles

if norm([th1 alfa1 n1]) < norm([th alfa n])
   x=[b;d;0;0;-alfa1;th1;0;0];
else
   x=[b;d;0;0;alfa;th;0;0];
end

% ---------------------------------------------------------------------------
% OTHER STRING

else

disp('   ');
disp('   x=T2x(T,str)');
disp('   where T is a 4 by 4 matrix (see help for details)');
disp('   and str can be : ''van'',''qua'',''erp'',''rpy'',''rpm'',''dht''. ');
disp('   ');

end

function z=vp(x,y)

% z=vp(x,y); z = 3d cross product of x and y
% vp(x) is the 3d cross product matrix : vp(x)*y=vp(x,y).
%
% by Giampiero Campa.

z=[  0    -x(3)   x(2);
    x(3)    0    -x(1);
   -x(2)   x(1)    0   ];

if nargin>1, z=z*y; end


text:
function LDR = openldr(FN,PERMISSION,Mode,arg4,arg5,arg6)
% OPENLDR loads neuroscan LDR files
% LDR = OPENLDR(Filename [, PERMISSION [, Mode]]);
% LDR = OPENLDR(LDR [, PERMISSION [, Mode]]);
%
% LDR is a struct with the following fields
%   LDR.FileName 	Name of LDR-file
%   LDR.Label_Out	Labels of output channels
%   LDR.Label_In	Labels of input channels
%   LDR.RR		re-referencing matrix
%   LDR.datatype	'REREF_MATRIX' indicates this datatype
%
% PERMISSION	'r'	reads LDR file
% 		'w'	writes LDR file
% 		'r+w'	reads and writes LDR file
%			(useful in combination with RESCALE-Mode)
%
% Mode [optional] 'RESCALE' performs a rescaling of the weights
%		sum of positive weights becomes +1
%		sum of negative weights becomes -1
%

%	$Id: openldr.m 2473 2010-06-04 09:56:57Z schloegl $
%	Copyright (C) 1997-2003,2008 by Alois Schloegl <a.schloegl@ieee.org>
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/

% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 3
% of the  License, or (at your option) any later version.

if nargin<2, PERMISSION=''; end;
if nargin<3, Mode=''; end;

if ~isstruct(FN),
        LDR.FileName = FN;
        PERMISSION = [PERMISSION,'r'];
else
        LDR = FN;
        PERMISSION = [PERMISSION,'w'];
end;

if any(PERMISSION=='r');
        % load LDR-file
        fid = fopen(LDR.FileName);
        if fid<0, fprintf(2,'File %s not found.\n',LDR.FileName); return; end;

        % reads line 1: size information
        s  = fgetl(fid);
        if ~ischar(s), fprintf(2,'ERROR LOADLDR: file %s corrupted\n',FN); end;
        sz = str2num(s);

        % reads line 2: output labels
        s  = fgetl(fid);
        if ~ischar(s), fprintf(2,'ERROR LOADLDR: file %s corrupted\n',FN); end;
        tmp = reshape(s,12,sz(2)+1)';
        LDR.Label_Out = tmp(2:sz(2)+1,:);

        % read lines 3+: input labels and weights
        for k = 1:sz(1),
                s = fgetl(fid);
                if ~ischar(s), fprintf(2,'ERROR LOADLDR: file %s corrupted\n',FN); end;
                r = reshape(s,12,sz(2)+1)';
                LDR.Label_In(k,1:12) = char(abs(s(1:12)));
                LDR.RR(k,1:sz(2)) = str2num(r(2:size(r,1),:))';
        end;

        fclose(fid);
        LDR.datatype='REREF_MATRIX';

end;

if strcmp(Mode,'RESCALE');
        tmp = LDR.RR;
        tmp(tmp<0) = 0;
        w = sum(tmp);
        RR = zeros(sz);
        if ~all(w==0 | w==1)
                ix = find(w);
                RR(:,ix) = tmp(:,ix)./w(ones(sz(1),1),ix);
                % RR(:,ix) = (tmp(:,ix)>0)./(ones(sz(1),1)*sum(tmp(:,ix)>0));
        end;
        tmp = LDR.RR;
        tmp(tmp>0)=0;
        w = sum(tmp);
        if ~all(w==0 | w==-1)
                ix = find(w);
                RR(:,ix) = tmp(:,ix)./w(ones(sz(1),1),ix);
                % RR(:,ix) = RR(:,ix)-(tmp(:,ix)<0)./(ones(sz(1),1)*sum(tmp(:,ix)<0));
        end;
        LDR.RR=RR;
end;

if any(PERMISSION=='w');
        tmp = isfield(LDR,'datatype') & strcmp(LDR.datatype,'REREF_MATRIX');
        if ~tmp,
                warning('LDR.datatype does not fit');
        end;
        sz = size(LDR.RR);
        if sz(1)~=size(LDR.Label_In,1);
                fprintf(2,'Size of Label_In does not fit RR\n');
                return;
        end;
        if sz(2)~=size(LDR.Label_Out,1);
                fprintf(2,'Size of Label_Out does not fit RR\n');
                return;
        end;

	% open LDR-file
        fid = fopen(LDR.FileName,'w+');
	if fid<0, fprintf(2,'Couldnot open file %s .\n',LDR.FileName); return; end;

        % write line 1: size information
        fprintf(fid,'%i %i\n',sz);

        % condition label information
        %LDR.Label_Out = char(LDR.Label_Out);

        nc = size(LDR.Label_Out,2);
        if nc<12,
                LDR.Label_Out = [LDR.Label_Out,abs(' ')*ones(sz(1),12-nc)];
        elseif nc>12,
                LDR.Label_Out = LDR.Label_Out(:,1:12);
        end;
        %LDR.Label_In = char(LDR.Label_In);
        nc = size(LDR.Label_In,2);
        if nc<12,
                LDR.Label_In = [LDR.Label_In,abs(' ')*ones(sz(1),12-nc)];
        elseif nc>12,
                LDR.Label_In = LDR.Label_In(:,1:12);
        end;

	% write line 2: out-labels
        fwrite(fid,32+zeros(12,1),'uint8');
        %fprintf(fid,'%c',abs(' ')*ones(12,1));
	for k = 1:sz(2),
                fwrite(fid,abs(LDR.Label_Out(k,1:12)),'uint8');
        end;

        % write lines 3+: in-labels and weights
        for k = 1:sz(1),
                fprintf(fid,'\n');
                fwrite(fid,abs(LDR.Label_In(k,1:12)),'uint8');
                fprintf(fid,'%12.5f',LDR.RR(k,:));
        end;
        fclose(fid);
end;

author: PatternRecognition

repo: OpenBMI


text:
function F = display(X)
% DISPLAY Overloaded

try
    X = sdpvar(X);
    display(X);
catch
    disp('Incomplete block variable.');
end


text:
function [Model, Info] = linear_map_sparse_cov(X,Y,Model,parm)
%  Estimate linear weight matrix for input-output mapping
%     Automatic Relevance Prior for each input dimension
%     is imposed to get sparse weight matrix
%
%   [Model, Info] = linear_map_sparse_cov(X,Y,Model,parm)
%
% --- Input
%  X  : Input data  ( M x T )
%  Y  : Output data ( N x T )
%  N  =  # of output
%  M  =  # of input
%  T  =  # of data
%
%  Model : Structure for estimated model
%  Model.SY0 :  Output data variance                 ( 1 x 1 )
%  Model.A0  :  (Output data var)/(Input data var)   ( 1 x 1 )
%
%  parm  : Structure for learning parameter
%  parm.Ntrain :  # of training
%  parm.Nskip  :  skip # for print
%  parm.a_min  :  Min value for pruning small variance component
%  parm.Prune  :  = 1 : Prune small variance & irrelevant input dimension
%
% --- Output
%  Model : Structure for estimated model
%  Model.SY  :  Noise variance         ( 1 x 1 )
%  Model.SW  :  Weight variance        ( M x M )
%  Model.W   :  Weight matrix          ( N x M )
%  Model.A   :  Prior weight variance  ( N x M ) ARD hyper parameter
%
%  Info  : Structure for learning process history
%  Info.FE  = LP + H : Free energy
%  Info.LP  = - (Log error)
%  Info.H   = - (# of effective weight parameters)
%
% 2007/1/26 Made by M. Sato

% Constants
MINVAL  = 1.0e-15;
MinCond = 1.0e-10;

% # of total training iteration
Ntrain = parm.Ntrain;

Nskip  = 100;   % skip steps for display info
a_min  = 1e-10; % Minimum value for weight pruning
Fdiff  = 1e-10; % Threshold for convergence
Ncheck = 100;   % Minimum number of training iteration
Fstep  = 5;     % Free energy convergence check step
Prune  = 1;     % Prune mode

if isfield(parm,'Nskip'), Nskip  = parm.Nskip; end;
if isfield(parm,'Fdiff'), Fdiff   = parm.Fdiff; end;
if isfield(parm,'a_min'), a_min   = parm.a_min ; end;
if isfield(parm,'Prune'), Prune = parm.Prune; end;
if isfield(parm,'Ncheck'), Ncheck = parm.Ncheck; end
if isfield(parm,'Fstep'),Fstep  = parm.Fstep ; end

% # of embedding dimension
if isfield(parm,'Dtau')
	D    = parm.Dtau;
	tau  = parm.Tau;
else
	D    = 1;
	tau  = 1;
end

% Dimension
[M ,Tx ,Nx ]= size(X); % input dim
[N ,Ty ,Ny ]= size(Y); % output dim

if Nx~=Ny, error('Trial number is different for input & output'); end
if Tx~=Ty, error('Time sample is different for input & output'); end

% Reshape into 2D matrix
T = Tx*Nx; % # of data
X = reshape(X, [M T]);
Y = reshape(Y, [N T]);

% # of stable VB-update in initial training
if isfield(parm,'Npre_train')
	Npre_train = parm.Npre_train;
else
	Npre_train = Ntrain;
end
if Npre_train > Ntrain, Npre_train = Ntrain; end;

fprintf('linear map sparse covariance start\n')
if 0 % hacked TH111005
    fprintf('--- Output Dimension  = %d\n',N)
    fprintf('--- Input  Dimension  = %d\n',M)
    fprintf('--- Embedding  Dimension  = %d\n',D)
    fprintf('--- Number of trials  = %d\n',Nx)
    fprintf('--- Number of training sample = %d\n',Tx)
    fprintf('--- Total update iteration    = %d (%d)\n',Ntrain,Npre_train)
else
    fprintf('---O:%d/I:%d/E:%d/Nt:%d/Ns:%d/Ti:%d(%d)/\n',N,M,D,Nx,Tx,Ntrain,Npre_train)
end
%
% --- Initialization
%  A  : Initial variable to use 1st update
%     : 1 x M

% Input/Output variance
sx = mean(repadd(X, - mean(X,2)).^2, 2);
sy = mean(repadd(Y, - mean(Y,2)).^2, 2);

A0  = 1./mean(sx);
SY0 = mean(sy);

if isempty(Model)
	A   = repmat(A0, [1,M]);
	W   = zeros(N,M);
	SY  = SY0;
else
	A   = Model.A ;	 % 1 x M
	W   = Model.W ;  % N x M
	SY  = mean(Model.SY);  % 1 x 1
end

if isfield(parm,'Ta0') && parm.Ta0 > 0,
	Ta0 = parm.Ta0;
	a0  = parm.a0 * A0;
else
	Ta0 = 0;
	a0  = 1;
end

%  --- Initialization by other method ---
% Model.mode = 'scalar': ARD term = alpha * W^2
% Model.mode = 'cov'   : ARD term = alpha * W^2 * SY(^-1)
%
if isfield(Model, 'mode') &&  strcmp(Model.mode,'scalar')==1,
	fprintf('Old result is used as initial value\n')
	fprintf('Old method = %s\n', Model.method)
	A  = sum(A,1)./(sum(SY));
end

A = max(A,MINVAL);

% Original input dimension
M_ALL = M;

if isfield(Model,'ix_act')
	% Active index
	IX_act = Model.ix_act;

	X = X(IX_act,:); 	% M x T
	M = length(IX_act);
else
	IX_act = 1:M;
end

% Initial active index
M_all  = M;
A_all  = A/max(A) ;
ix_act_old = 1:M;

ix_act = find( A_all > a_min );   % effective indices
Mnew   = length(ix_act);  		% # of effective input

if Mnew < M,
    % convert to relative index
    jx_act = trans_index(ix_act,ix_act_old,M_all);

    M   = Mnew;
    A   = A(jx_act) ;  			% 1 x M
    W   = W(:,jx_act) ;  		% N x M
	X 	= X(jx_act,:);			% M x T
end

% Input covariance
%XX  = (X * X')/T;   	% M x M
%YX  = (Y * X')/T;       % N x M
%YY  = sum(Y.^2,2)/T;    % N x 1

% Covariance matrix (not normalised)
YX  = (Y * X');       % N x M
YY  = sum(Y.^2,2);    % N x 1

if 0 % hacked TH111005
fprintf('a_min = %g\n', a_min)
fprintf('SY0   = %g\n', SY0)
fprintf('SY    = %g\n', SY)
end
% Working variable
if T <= M
	XX = [];
	CC = zeros(T,T);
	SW = zeros(T,T);
else
	XX = (X * X');   	  % M x M
	CC = zeros(M,M);
	SW = zeros(M,M);
end

G_A = zeros(1,M);       % 1 x M
log_a = 0;
A_old = A;

% Free energy histry
FE  = zeros(Ntrain,1);
LP  = zeros(Ntrain,1);
H   = zeros(Ntrain,1);
MM  = zeros(Ntrain,1);
Err = zeros(Ntrain,1);

% ARD hyper param. history
if isfield(parm,'Debug') & ~isempty(parm.Debug) & parm.Debug > 0
	Debug = 1;
	A_tmp = zeros(M_all, ceil(Ntrain/Nskip));
else
	Debug = 0;
end

k_save  = 0;

%%%%%% Learning Loop %%%%%%
for k=1:Ntrain
	% ARD hyper variance parameter
	% A = 1/alpha

	if T < M
	    % Weight variance
	    % inv(X*X' + 1./A) = A - A * X * inv(X'*A*X + 1) * X' * A
		%  C = ( X' *A* X + eye(T) );
		XA = repmultiply(X' , A); % T x M
	    CC = XA * X + eye(T);  % T x T

		% Weight update
	    % inv(X*X' + 1./A) = A - A * X * Cinv * X' * A
		% W0 = YX .* A;
		% W  = W0 - (((W0 * X) * Cinv) * X') .* A;
		% W  = W0 - ((W0 * X) / C ) * XA;
		W  = repmultiply(YX , A);
		XC = X / CC;
		W  = W - (W * XC) * XA;

		%  G_A = diag( X * inv(C) * X' *A )
		%      = diag( (X / C) * X') .*A
		G_A  = A .* sum(X .* XC, 2)';

		% Log variance
		log_sw  = - log_det(CC) ;
		if mod(k, Nskip)==0, fprintf('- '); end
	else
		if isempty(XX)
			% covariance matrix in reduced space
			XX = X * X';
			% save original index
			IX_act = IX_act(ix_act);
			% new active index in reduced space
			ix_act = 1:M;
			M_all  = M;
			A_all  = A;
		end

		% Weight covariance
    	SW  = XX + diag(1./A);

		% Weight update
		W  = YX / SW;

		% SW  = XX + diag(1./A)
		% G_A = diag( XX * inv(SW))
		%     = 1 - diag(inv(SW)) ./A
		G_A = diag( XX /SW )';

		log_sw  = - log_det(SW) - sum(log(A));
		if mod(k, Nskip)==0, fprintf('+ '); end
	end

	WW = sum(W.^2, 1);
    % Noise variance update
    SY = (sum(YY) - sum(sum(W.*YX)))/(N*T);

    if (SY/SY0) <= MINVAL,
	    % Error
	    dY  = Y - W * X;        % N x T
	    dYY = sum(dY.^2, 2);  	% N x 1

	    SY  = (sum(dYY) + sum( WW./A ))/(N*T);
	    % Prevent zero variance
	    SY  = max( SY, MINVAL);
	    fprintf('*')
	end

	% Log variance
    log_sy  = N * log(SY) ;

    if Ta0 > 0,
	    log_a   = Ta0 * sum( - log(A./a0) - a0./A + 1);
	end

    % Free energy
    H(k)   =   0.5*N * (log_sw - M);
    LP(k)  = - 0.5*(T * log_sy) ;
    FE(k)  = LP(k) + H(k);
    Err(k) = (SY)./(SY0);
    MM(k)  = M;

    % Estimation Gain
    %jx_act = find(G_A > a_min);
    G_A = max((G_A), MINVAL);

    % Hyper parameter for weight variance (ARD)
	if k <= Npre_train,
		% VB update rule
		% N * A  = (1./SY)' * (W.^2) + N * (A - A.*G_A)  ;
		% A  = (WW./SY + N * (A - A.*G_A) + 2*Ta0*a0)./( N + 2*Ta0 );
		% A^2 = A .* (WW./SY) ./ (G_A * N);
        try % modified by TH130111
            A2=A;
            A  = sqrt(A.*(WW./SY)./(G_A * N));
        catch
            A=A2;
        end
	else
	    % Accelerated update rule
		%	A  = (1./SY)' * (W.^2) ./ (G_A * N);
	    A  = ((WW./SY) + 2*Ta0*a0)./(G_A * N + 2*Ta0);
	end

    % Prune small variance
    if Prune == 1
	    % Find active input dimension
	    ix_act_old = ix_act;

	    % Recover all component
	    switch	Prune
	    case	1
		    A_all(ix_act) = WW/max(WW);    % Prune by Weight
	    case	2
		    A_all(ix_act) = A /max(A);	   % Prune by Alpha
	    case	3
		    A_all(ix_act) = A * (1/SX);    % Prune by Alpha
	    end

	    % Find active input dimension (absolute index)
	    ix_act = find( A_all > a_min ); % effective indices
	    Mnew   = length(ix_act);  		% # of effective input

	    if Mnew < M,
		    % convert to relative index
		    jx_act = trans_index(ix_act,ix_act_old,M_all);

		    M   = Mnew;
		    A   = A(jx_act) ;  			% 1 x M
		    W   = W(:,jx_act) ;  		% N x M
			X 	= X(jx_act,:);			% M x T
			YX  = YX(:,jx_act);  	 	% N x M
			if ~isempty(XX)
				XX	= XX(jx_act,jx_act);	% M x M
			end
		end
	else
		A = max(A,MINVAL);
    end

    if mod(k, Nskip)==0
        % Save history
		if Debug == 1
        	k_save = k_save + 1;
        	A_tmp(:,k_save) = A_all(:);
		end

        fprintf('Iter = %4d, M = %4d, err = %g, F = %g, H = %g\n', ...
               k, M, Err(k), FE(k), - H(k));
    end

	if k > Ncheck && M == MM(k-1)
		Adif = max(abs(A - A_old));
	else
		Adif = 1;
	end
	if Adif < Fdiff,
		fprintf('Converged : Alpha change = %g\n',Adif)
		break;
	end;

	A_old = A;

%		Fdif = (FE(k) - FE(k-Fstep))/abs(FE(k));
%	else
%		Fdif = Fdiff + 1;
%	end
end

% convert to relative index
ix_act = IX_act(ix_act);

% Active index
Model.ix_act = ix_act;
Model.M_all  = M_ALL ;

% Save trained variable
%  W & A is sufficient for cov-method initialization
Model.A  = A ;
Model.W  = W ;
Model.SY = SY;

Model.method = 'linear_map_sparse_cov';
Model.mode   = 'cov';
Model.sparse = 'sparse';

% Save history
Info.FE  = FE(1:k);
Info.LP  = LP(1:k);
Info.H   = H(1:k) ;
Info.Err = Err(1:k);
Info.M   = MM(1:k);

if exist('A_tmp','var')
	Info.A   = A_tmp(:,1:k_save) ;
end

%%% ---- Index transformation from old active_index to current active_index
function	jx = trans_index(ix,ix_old,M)
% ix = ix_old(jx)

N = length(ix_old);
Itrans = zeros(M,1);
Itrans(ix_old) = 1:N;

jx = Itrans(ix);

author: KamitaniLab

repo: GenericObjectDecoding


text:
function cr = pwidentification(simmat, labels)
% pwidentification    Perform pairwise identification from similarity matrix
%
% Inputs:
%
% - simmat : Similarity matrix [pred x correct]
% - labels : Index matrix
%
% Outputs:
%
% - cr : Correct rate
%
%
% Author: Tomoyasu Horikawa <horikawa-t@atr.jp>, Shuntaro C. Aoki <aoki@atr.jp>
% Created: 2016-09-10
% Modified: 2016-10-20
%

% Fix labels to a vertical vector
if size(labels, 1) == 1
    labels = labels';
end

% Remove NaN
if any(any(isnan(simmat), 1))
    simmat(:, any(isnan(simmat), 1)) = [];
end

% Get num of candidate
numCandidate = size(simmat, 2) - 1;

% Sort simmat for each prediction
[sortedSimmat, order] = sort(simmat, 2, 'descend');

% Get num of incorrect
[labelList, numIncorrect] = find(~bsxfun(@minus, order, labels));
numIncorrect = numIncorrect - 1;

% Get index of the original label list.
[sortedLabelList, ind] = sort(labelList, 'ascend');

% Calculate correct rate
cr = (numCandidate - numIncorrect(ind)) ./ numCandidate;
% `cr = 1 - numIncorrect(labelInd) ./ numCandidate` caused slight numerical error
% compared with the original code (~ 10^-16).


text:
function z = svec( x, nrm )

if nargin < 2 || isempty( nrm ) || isequal( nrm, 'fro' ),
    nrm = 2;
elseif ~isnumeric( nrm ) || length( nrm ) ~= 1 || nrm < 1,
    error( 'Second argument must be a number between 1 and Inf, or ''fro''.' );
end

if ~isreal( x ) && nrm ~= 2,
    z = vec( x );
    return
else
    [ xR, y ] = bcompress( x );
    if isempty( y ),
        z = cvx( 0 );
    else
        z = y .* norms( xR, nrm, 2 );
    end
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.