text:
function [x,testdata]=alstpz_solve(A,y,tol,varargin)
%Solution of linear systems in TT-format via ALS(t+z) iteration
%!Note! This method is provided mostly for test purposes.
%!Note! You may prefer to use amen_solve2 instead.
%   [X,somedata]=ALSTPZ_SOLVE(A,Y,TOL,OPTIONS) Attempts to solve the linear
%   system A*X = Y with accuracy/residual TOL using the global (t+z)-enriched ALS.
%   Matrix A has to be given in the TT-format, right-hand side Y should be
%   given in the TT-format also. Options are provided in form
%   'PropertyName1',PropertyValue1,'PropertyName2',PropertyValue2 and so
%   on. The parameters are set to default (in brackets in the following)
%   The list of option names and default values are:
%       o x0 - initial approximation [random rank-2 tensor]
%       o nswp - maximal number of sweeps [10]
%       o max_full_size - maximal size of the local matrix to full solver [50]
%       o local_prec - local preconditioner: '' (no prec.), 'ljacobi',
%         'cjacobi', 'rjacobi' ['']
%       o local_iters - number of local gmres restarts [2]
%       o local_restart - dimension of local gmres [40]
%       o kickrank - compression rank of the residual Z, i.e. expansion
%         size [4]
%       o kicktype - how to approximate the residual: via the SVD rounding ('svd'),
%         or using one approximation ALS sweep ('als'). ['svd']
%       o ismex - shall we use the MEX lib solve3d_2 for local solution [true]
%       o resid_damp - solve local problems with accuracy tol/resid_damp.
%         Larger value may reduce a spurious noise from inexact local
%         solutions, but increase CPU time [2]
%       o symm - shall we consider the symmetrized system (A'A)x=(A'y). [false]
%
%       Example:
%           d=8; f=8;
%           mat=tt_qlaplace_dd(d*ones(1,f)); % Laplace in the QTT-format
%           rhs=tt_ones(2,d*f); % Right-hand side of all ones
%           sol = alstpz_solve(mat, rhs, 1e-5); % solve the Poisson eqn.
%
%********
%   References:
%   S. Dolgov, D. Savostyanov.
%   http://arxiv.org/abs/1301.6068
%   http://arxiv.org/abs/1304.1222
%
%   Please send feedback to: {sergey.v.dolgov,dmitry.savostyanov}@gmail.com
%
%********
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------
% Inner parameters
max_full_size=50;

resid_damp = 2; % Truncation error to true residual treshold

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

local_prec = '';

ismex = true;

% kicktype = 'als';
kicktype = 'svd';

kickrank = 4;
x=[];
symm = false;
verb = 3;

for i=1:2:length(varargin)-1
    switch lower(varargin{i})
        case 'nswp'
            nswp=varargin{i+1};
        case 'x0'
            x=varargin{i+1};
        case 'local_prec'
            local_prec=varargin{i+1};
        case 'local_restart'
            local_restart=varargin{i+1};
        case 'ismex'
            ismex=varargin{i+1};
        case 'local_iters'
            local_iters=varargin{i+1};
        case 'kickrank'
            kickrank=varargin{i+1};
        case 'kicktype'
            kicktype=varargin{i+1};
        case  'max_full_size'
            max_full_size=varargin{i+1};
        case 'resid_damp'
            resid_damp = varargin{i+1};
        case 'symm'
            symm=varargin{i+1};
        case 'verb'
            verb=varargin{i+1};

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

% Disable MEX if it does not exist
if (ismex)&&(exist('solve3d_2', 'file')<2)
    warning('MEX local solver is not found, disabled');
    ismex = false;
end;

% Extract sizes and TT blocks
d = y.d;
n = A.n;
if (isempty(x))
    x = tt_rand(n, A.d, 2);
end;

if (symm)
    % Symmetrize the system if necessary
    A2 = A'*A;
    Ay = A'*y;
end;

if (strcmp(kicktype, 'als'))&&(kickrank>0)
    z = tt_rand(n, d, kickrank);
    rz = z.r;
    phi = cell(d+1,1); phi{1}=1; phi{d+1}=1;
end;

% Test output data -- e.g. for convergence tracking
testdata = cell(3,1);
testdata{2} = cell(d, nswp);
testdata{1} = zeros(d, nswp);
testdata{3} = zeros(1, nswp);
t_corr_amen = tic;

% ALS(t+z) sweeps
for swp=1:nswp
    % Truncate the solution first
    x = round(x,tol);
    if (strcmp(kicktype, 'als'))&&(kickrank>0)
        % Update the residual approximation via one ALS sweep for |z-ZZ|->min
        % Exact residual
        if (symm)
            ZZ = Ay - A2*x;
        else
            ZZ = y - A*x;
        end;
        rZZ = ZZ.r;
        ZZ = core2cell(ZZ);
        z = core2cell(z); % Approximate one
        for i=1:d-1
            % Go backward and prepare projections
            cr = reshape(z{i}, rz(i)*n(i), rz(i+1));
            [cr,rv]=qr(cr,0);
            cr2 = reshape(z{i+1}, rz(i+1), n(i+1)*rz(i+2));
            cr2 = rv*cr2;
            rz(i+1) = size(cr, 2);
            z{i} = reshape(cr, rz(i), n(i), rz(i+1));
            z{i+1} = reshape(cr2, rz(i+1), n(i+1), rz(i+2));

            phi{i+1} = compute_next_Phi(phi{i}, z{i}, [], ZZ{i}, 'lr');
        end;
        for i=d:-1:1
            % Go forward and project exact ZZ
            cr = reshape(ZZ{i}, rZZ(i)*n(i), rZZ(i+1));
            cr = cr*phi{i+1};
            cr = reshape(cr, rZZ(i), n(i)*rz(i+1));
            cr = phi{i}*cr;
            if (i>1)
                [cr,rv]=qr(cr.', 0);
                cr2 = reshape(z{i-1}, rz(i-1)*n(i-1), rz(i));
                cr2 = cr2*rv.';
                rz(i) = size(cr, 2);
                z{i} = reshape(cr.', rz(i), n(i), rz(i+1));
                z{i-1} = reshape(cr2, rz(i-1), n(i-1), rz(i));

                phi{i} = compute_next_Phi(phi{i+1}, z{i}, [], ZZ{i}, 'rl');
            else
                z{i} = reshape(cr, rz(i), n(i), rz(i+1));
            end;
        end;
        z = cell2core(tt_tensor, z);
    else
        % SVD rounding of the exact residual up to kickrank
        if (kickrank>0)
            if (symm)
                z = round(Ay - A2*x, tol, kickrank);
            else
                z = round(y - A*x, tol, kickrank);
            end;
        end;
    end;

    % Augment the solution with Z basis
    if (kickrank>0)
        x = x+0*z;
    end;
    t_old = toc(t_corr_amen);
    % Update via the ALS iteration
    if (symm)
        [x,td_als] = als_solve(A2, Ay, tol,  'max_full_size', max_full_size, 'x0', x, 'verb', verb, 'local_prec', local_prec, 'local_iters', local_iters, 'local_restart', local_restart, 'resid_damp', resid_damp, 'ismex', ismex);
    else
        [x,td_als] = als_solve(A, y, tol,  'max_full_size', max_full_size, 'x0', x, 'verb', verb,  'local_prec', local_prec, 'local_iters', local_iters, 'local_restart', local_restart, 'resid_damp', resid_damp, 'ismex', ismex);
    end;
    testdata{1}(:, swp) = t_old+td_als{1}(:,1);
    testdata{2}(:, swp) = td_als{2};
    testdata{3}(swp) = max(td_als{3});
end;

end

function [x,td]=als_solve(A, y, tol, varargin)
% One sweep of the one-site DMRG solution scheme.
% Uses the same parameters as in the main function

local_prec_char = 0;
verb = 1;
ismex = true;

x=[];

for i=1:2:length(varargin)-1
    switch lower(varargin{i})
        case 'x0'
            x=varargin{i+1};
        case 'verb'
            verb=varargin{i+1};
        case 'local_prec'
            local_prec=varargin{i+1};
        case 'local_restart'
            local_restart=varargin{i+1};
        case 'local_iters'
            local_iters=varargin{i+1};
        case 'ismex'
            ismex=varargin{i+1};
        case  'max_full_size'
            max_full_size=varargin{i+1};
        case 'resid_damp'
            resid_damp = varargin{i+1};

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

if (strcmp(local_prec, 'cjacobi')); local_prec_char = 1;  end;
if (strcmp(local_prec, 'ljacobi')); local_prec_char = 2;  end;
if (strcmp(local_prec, 'rjacobi')); local_prec_char = 3;  end;

% Extract sizes and TT cores
d = y.d;
n = A.n;
if (isempty(x))
    x = tt_rand(n, A.d, 2);
end;

ry = y.r;
ra = A.r;
rx = x.r;

cry = core2cell(y);
crA = core2cell(A);
crx = core2cell(x);

% Interfaces
phia = cell(d+1,1); phia{1}=1; phia{d+1}=1;
phiy = cell(d+1,1); phiy{1}=1; phiy{d+1}=1;

% This is some convergence output for test purposes
td = cell(3,1);
td{1} = zeros(d, 1);
td{2} = cell(d, 1);
td{3} = zeros(d, 1);

t_amr_solve = tic;

% Orthogonalization
for i=d:-1:2
    cr = crx{i};
    cr = reshape(cr, rx(i), n(i)*rx(i+1));
    [cr, rv]=qr(cr.', 0);
    cr2 = crx{i-1};
    cr2 = reshape(cr2, rx(i-1)*n(i-1), rx(i));
    cr2 = cr2*(rv.');
    rx(i) = size(cr, 2);
    cr = reshape(cr.', rx(i), n(i), rx(i+1));
    crx{i-1} = reshape(cr2, rx(i-1), n(i-1), rx(i));
    crx{i} = cr;

    phia{i} = compute_next_Phi(phia{i+1}, cr, crA{i}, cr, 'rl');
    phiy{i} = compute_next_Phi(phiy{i+1}, cr, [], cry{i}, 'rl');
end;

max_res = 0;
max_dx = 0;

% DMRG sweep
for i=1:d
    % Extract elements - matrix
    Phi1 = phia{i}; Phi2 = phia{i+1};
    % Phi1: rx'1, rx1, ra1, or rx'1, ry1
    % Phi2: rx2, ra2, rx'2, or ry'2, rx2
    A1 = crA{i}; y1 = cry{i};
    % RHS - rewrite it in accordance with new index ordering
    rhs = phiy{i}; % rx'1, ry1
    y1 = reshape(y1, ry(i), n(i)*ry(i+1));
    rhs = rhs*y1;
    rhs = reshape(rhs, rx(i)*n(i), ry(i+1));
    rhs = rhs*phiy{i+1};
    rhs = reshape(rhs, rx(i)*n(i)*rx(i+1),1);
    norm_rhs = norm(rhs);
    % sol_prev
    sol_prev = reshape(crx{i}, rx(i)*n(i)*rx(i+1), 1);

    real_tol = (tol/sqrt(d))/resid_damp;

    if (rx(i)*n(i)*rx(i+1)<max_full_size) % Full solution
        %      |     |    |
        % B = Phi1 - A1 - Phi2
        %      |     |    |
        B = reshape(Phi1, rx(i)*rx(i), ra(i));
        B = B*reshape(A1, ra(i), n(i)*n(i)*ra(i+1));
        B = reshape(B, rx(i), rx(i), n(i), n(i)*ra(i+1));
        B = permute(B, [1, 3, 2, 4]);
        B = reshape(B, rx(i)*n(i)*rx(i)*n(i), ra(i+1));
        B = B*reshape(permute(Phi2, [2, 3, 1]), ra(i+1), rx(i+1)*rx(i+1));
        B = reshape(B, rx(i)*n(i), rx(i)*n(i), rx(i+1), rx(i+1));
        B = permute(B, [1, 3, 2, 4]);
        B = reshape(B, rx(i)*n(i)*rx(i+1), rx(i)*n(i)*rx(i+1));

        res_prev = norm(B*sol_prev - rhs)/norm_rhs;

        if (res_prev>real_tol)
            sol = B \ rhs;
            res_new = norm(B*sol-rhs)/norm_rhs;
        else
            sol = sol_prev;
            res_new = res_prev;
        end;

    else % Structured solution.

        res_prev = norm(bfun3(Phi1, A1, Phi2, sol_prev) - rhs)/norm_rhs;

        if (res_prev>real_tol)
            if (~ismex)
                sol = solve3d_2ml(Phi1, A1, Phi2, rhs, real_tol, sol_prev, local_prec_char, local_restart, local_iters);
            else % use MEX
                sol = solve3d_2(Phi1, A1, Phi2, rhs, real_tol, 1, sol_prev, local_prec_char, local_restart, local_iters, 0);
            end;

            res_new = norm(bfun3(Phi1, A1, Phi2, sol) - rhs)/norm_rhs;
        else
            sol = sol_prev;
            res_new = res_prev;
        end;

    end;

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

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

    td{3}(i) = dx;

    % QR
    sol = reshape(sol, rx(i)*n(i), rx(i+1));
    [u,v]=qr(sol, 0);
    r = size(u,2);

    if (verb==2)
        fprintf('=als_solve=   block %d, dx: %3.3e, res: %3.3e, r: %d\n', i, dx, res_prev, r);
    end;

    if (i<d) % left-to-right, kickrank, etc
        cr2 = crx{i+1};
        cr2 = reshape(cr2, rx(i+1), n(i+1)*rx(i+2));
        v = v*cr2; % size r+radd, n2, r3

        u = reshape(u, rx(i), n(i), r);
        v = reshape(v, r, n(i+1), rx(i+2));

        % Recompute phi. Left ones, so permute them appropriately
        phia{i+1} = compute_next_Phi(phia{i}, u, crA{i}, u, 'lr');
        phiy{i+1} = compute_next_Phi(phiy{i}, u, [], cry{i}, 'lr');

        % Stuff back
        rx(i+1) = r;
        crx{i} = u;
        crx{i+1} = v;
    else
        % Just stuff back the last core
        sol = reshape(sol, rx(i), n(i), rx(i+1));
        crx{i} = sol;
    end;

    if (verb>2)
        td{1}(i) = toc(t_amr_solve);
        if (verb>3)||(i==d) % each microstep is returned only if really asked for
            % Otherwise the memory will blow up
            x = cell2core(x, crx); % for test
            td{2}{i} = x;
        end;
    end;

end;

if (verb>0)
     fprintf('=als_solve= max_dx: %3.3e, max_res: %3.3e, erank: %g\n', max_dx, max_res, sqrt(rx(1:d)'*(n.*rx(2:d+1))/sum(n)));
end;

x = cell2core(x, crx);

end

function [Phi] = compute_next_Phi(Phi_prev, x, A, y, direction)
% Performs the recurrent Phi (or Psi) matrix computation
% Phi = Phi_prev * (x'Ay).
% If direction is 'lr', computes Psi
% if direction is 'rl', computes Phi
% A can be empty, then only x'y is computed.

% Phi1: rx1, ry1, ra1, or rx1, ry1
% Phi2: ry2, ra2, rx2, or ry2, rx2

rx1 = size(x,1); n = size(x,2); rx2 = size(x,3);
ry1 = size(y,1); m = size(y,2); ry2 = size(y,3);
if (~isempty(A))
    ra1 = size(A,1); ra2 = size(A,4);
else
    ra1 = 1; ra2 = 1;
end;

if (strcmp(direction, 'lr'))
    %lr: Phi1
    x = reshape(x, rx1, n*rx2);
    Phi = reshape(Phi_prev, rx1, ry1*ra1);
    Phi = x'*Phi;
    if (~isempty(A))
        Phi = reshape(Phi, n*rx2*ry1, ra1);
        Phi = Phi.';
        Phi = reshape(Phi, ra1*n, rx2*ry1);
        A = reshape(A, ra1*n, m*ra2);
        Phi = A.'*Phi;
        Phi = reshape(Phi, m, ra2*rx2*ry1);
    else
        Phi = reshape(Phi, n, rx2*ry1);
    end;
    Phi = Phi.';
    Phi = reshape(Phi, ra2*rx2, ry1*m);

    y = reshape(y, ry1*m, ry2);
    Phi = Phi*y;
    if (~isempty(A))
        Phi = reshape(Phi, ra2, rx2*ry2);
        Phi = Phi.';
    end;
    Phi = reshape(Phi, rx2, ry2, ra2);
else
    %rl: Phi2
    y = reshape(y, ry1*m, ry2);
    Phi = reshape(Phi_prev, ry2, ra2*rx2);
    Phi = y*Phi;
    if (~isempty(A))
        Phi = reshape(Phi, ry1, m*ra2*rx2);
        Phi = Phi.';
        Phi = reshape(Phi, m*ra2, rx2*ry1);
        A = reshape(A, ra1*n, m*ra2);
        Phi = A*Phi;
        Phi = reshape(Phi, ra1*n*rx2, ry1);
        Phi = Phi.';
    end;

    Phi = reshape(Phi, ry1*ra1, n*rx2);
    x = reshape(x, rx1, n*rx2);
    Phi = Phi*x';
    if (~isempty(A))
        Phi = reshape(Phi, ry1, ra1, rx1);
    else
        Phi = reshape(Phi, ry1, rx1);
    end;
end;

end

% new
function [y]=bfun3(Phi1, A, Phi2, x)
% Phi1: ry1, rx1, ra1
ry1 = size(Phi1,1);
rx1 = size(Phi1,2);
ra1 = size(Phi1,3);
% Phi2: rx2, ra2, ry2
ry2 = size(Phi2,3);
rx2 = size(Phi2,1);
ra2 = size(Phi2,2);

n = size(A,2);
m = size(A,3);

y = reshape(x, rx1*m, rx2);
Phi2 = reshape(Phi2, rx2, ra2*ry2);
y = y*Phi2;
y = reshape(y, rx1, m*ra2*ry2);
y = y.';
y = reshape(y, m*ra2, ry2*rx1);
A = reshape(A, ra1*n, m*ra2);
y = A*y;
y = reshape(y, ra1*n*ry2, rx1);
y = y.';
y = reshape(y, rx1*ra1, n*ry2);
Phi1 = reshape(Phi1, ry1, rx1*ra1);
y = Phi1*y;
y = reshape(y, ry1*n*ry2, 1);
end

author: oseledets

repo: TT-Toolbox


text:
function line_integrals_test01 ( )

%*****************************************************************************80
%
%% LINE_INTEGRALS_TEST01 compares exact and estimated monomial integrals.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2014
%
%  Author:
%
%    John Burkardt
%
  m = 1;
  n = 4192;
  test_num = 11;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LINE_INTEGRALS_TEST01\n' );
  fprintf ( 1, '  Use LINE01_SAMPLE to estimate integrals\n' );
  fprintf ( 1, '  along the length of the unit line in 1D.\n' );
%
%  Get sample points.
%
  seed = 123456789;
  [ x, seed ] = line01_sample ( n, seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of sample points used is %d\n', n );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   E     MC-Estimate      Exact           Error\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    e = test - 1;

    value = monomial_value_1d ( n, e, x );

    result = line01_length ( ) * sum ( value(1:n) ) / n;
    exact = line01_monomial_integral ( e );
    error = abs ( result - exact );

    fprintf ( 1, '  %2d  %14.6g  %14.6g  %10.2e\n', ...
      e, result, exact, error );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ raw ] = RawFrom4Channels( fourchan )
%RAWTO4CHANNELS Summary of this function goes here
%   Detailed explanation goes here

[h,w,~] = size(fourchan);

raw = zeros(h*2,w*2);

idx = [1,1; 1,2; 2,1; 2,2];

for c=1:4
    raw(idx(c,1):2:end, idx(c,2):2:end) = fourchan(:,:,c);
end

end

author: AbdoKamel

repo: sidd-ground-truth-image-estimation


text:
function [scan3dd] = elec_load_scan_3ddasc(filename)

% elec_load_scan_3ddasc - read ascii export of Neuroscan 3DD file
%
% Load an ascii electrode file and return
% the electrode labels and coordinates.  This
% function can read and return Cartesian (x,y,z)
% and/or spherical (theta,phi,r) coordinates.
%
% [scan3dd] = elec_load_scan_3ddasc(filename)
%
% where:
%
% file = '<path><filename>' with row format '%s %d %f %f %f'
%
% The file format is that of NeuroScan 3Dspace ascii
% export files.  Each row of the file comprises an electrode label,
% an electrode type code (see below), and the x,y,z coordinates (cm).
% Each field is separated by spaces.  For example,
%
% Fz    69  Xcm Ycm Zcm
%
% Example result:
%
% scan3dd =
%
%      label: {1x128 cell}
%          x: [128x1 double]
%          y: [128x1 double]
%          z: [128x1 double]
%        hsp: [1617x3 double]
%     nasion: [-0.0333 9.0341 0]
%        lpa: [-6.9128 0 0]
%        rpa: [6.9128 0 0]
%     origin: [0 0 0]
%        ref: [0 9.9341 0]
%
% All coordinates are in centimeters.
%
% Notes:
%
% i) Type is defined as follows (from Neuroscan 3Dspace ascii export):
%
%           Electrode               Type
%           ---------               ----
%           Nasion                   110 (or 78)
%           Left                     108 (or 76)
%           Right                    114 (or 82)
%           electrodes                69
%           Centroid (origin)         99 (or 67; eg <0,0,0>)
%           Ref                      120 (or 88)
%           Scalp Points (hsp)        32
%
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:55 $

% Licence:  GNU GPL, no express or implied warranties
% History:  08/1999, Darren.Weber_at_radiology.ucsf.edu
%           08/2003, Darren.Weber_at_radiology.ucsf.edu
%                    complete rewrite, replacing elec_load.m
%                    new version of 3Dspace exports different type numbers ;-)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ver = '$Revision: 1.1 $';
fprintf('\nELEC_LOAD_SCAN3DDASC [v %s]\n',ver(11:15));

tic;

[path,name,ext] = fileparts(filename);
file = fullfile(path,[name ext]);

fprintf('...loading electrodes from:\n\t%s\n', file);

scan3dd = read_3dd(file);

fprintf('...loaded %d electrodes\n', size(scan3dd.x,1));

t = toc; fprintf('...done (%6.2f sec).\n\n',t);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function scan3dd = read_3dd(file),

scan3dd.label = [];
scan3dd.x = [];
scan3dd.y = [];
scan3dd.z = [];
scan3dd.hsp = [];

fid = fopen(file);
if fid < 0,
  msg = sprintf('cannot open file: %s\n',file);
  error(msg);
end

% 3DD ascii files contain position information for
% fiducial, sensor, and head shape fields.

% First get the fiducial points by reading the whole file
% (clumsy, but exhaustive search for field indicators).

% Fiducial points are required for MRI registration. They are
% the nasion, left and right preauricular points, eg:

%      Nasion	78	-0.033340	9.034104	0.000000
%      Left	76	-6.912818	-0.000000	0.000000
%      Right	82	6.912818	0.000000	-0.000000

fprintf('...searching for fiducials...');

n = 0;
while n < 5,
  tmp = fgetl(fid);
  if tmp < 0, break; end

  tmp = lower(tmp);

  if strfind(tmp,'nasion'),
    tmp = sscanf(tmp,'%s %d %f %f %f');
    scan3dd.nasion = [tmp(end-2) tmp(end-1) tmp(end)];
    n = n + 1;
    continue;
  end
  if strfind(tmp,'left'),
    tmp = sscanf(tmp,'%s %d %f %f %f');
    scan3dd.lpa = [tmp(end-2) tmp(end-1) tmp(end)];
    n = n + 1;
    continue;
  end
  if strfind(tmp,'right'),
    tmp = sscanf(tmp,'%s %d %f %f %f');
    scan3dd.rpa = [tmp(end-2) tmp(end-1) tmp(end)];
    n = n + 1;
    continue;
  end
  if strfind(tmp,'centroid'),
    tmp = sscanf(tmp,'%s %d %f %f %f');
    scan3dd.origin = [tmp(end-2) tmp(end-1) tmp(end)];
    n = n + 1;
    continue;
  end
  if strfind(tmp,'ref'),
    tmp = sscanf(tmp,'%s %d %f %f %f');
    scan3dd.ref = [tmp(end-2) tmp(end-1) tmp(end)];
    n = n + 1;
    continue;
  end
end

frewind(fid);
fprintf('done\n');

fprintf('...searching for electrodes...');

ok = 1;
while ok,
  tmp = fgetl(fid);
  if tmp < 0, break; end

  if strfind(lower(tmp),'nasion'),   continue; end
  if strfind(lower(tmp),'left'),     continue; end
  if strfind(lower(tmp),'right'),    continue; end
  if strfind(lower(tmp),'centroid'), continue; end
  if strfind(lower(tmp),'ref'),      continue; end

  tmp = sscanf(tmp,'%s %d %f %f %f');

  if tmp(end-3) == 69,
    scan3dd.label{end+1} = char(tmp(1:end-4))';
    scan3dd.x(end+1,1) = tmp(end-2);
    scan3dd.y(end+1,1) = tmp(end-1);
    scan3dd.z(end+1,1) = tmp(end);
    continue;
  end
  if strfind(char(tmp(1:end-4))','32'),
    break;
    % found a head shape point
    scan3dd.hsp(end+1,:) = [tmp(end-2),tmp(end-1),tmp(end)];
    continue;
  end
end

frewind(fid);
fprintf('done\n');

fprintf('...searching for head shape points...');

ok = 1;
while ok,
  tmp = fgetl(fid);
  if tmp < 0, break; end

  if strfind(lower(tmp),'nasion'),   continue; end
  if strfind(lower(tmp),'left'),     continue; end
  if strfind(lower(tmp),'right'),    continue; end
  if strfind(lower(tmp),'centroid'), continue; end
  if strfind(lower(tmp),'ref'),      continue; end
  if strfind(lower(tmp),'69'),       continue; end

  tmp = sscanf(tmp,'%d %f %f %f');

  if tmp(end-3) == 32,
    % found a head shape point
    scan3dd.hsp(end+1,:) = [tmp(end-2),tmp(end-1),tmp(end)];
    continue;
  end
end

frewind(fid);
fprintf('done\n');

fclose(fid);

return

author: PatternRecognition

repo: OpenBMI


text:
function DCM = VBA_spm_dcm_test(DCM)
% tests whether DCM parameters are zero using Savage-Dickey ratios
% function DCM = spm_dcm_test(DCM)
% IN:
%   - DCM: the DCM structure
% OUT:
%   - DCM: the DCM structure, augmented with field .sdr, which contains the
%   Savage-Dickey ratios in the same format as DCM.Ep

n = size(DCM.Ep.A,1);
nu = size(DCM.Ep.C,2);

E = VBA_spm_vec(DCM.Ep);
V = DCM.Cp;
E0 = VBA_spm_vec(DCM.M.pE);
V0 = DCM.M.pC;

sdr.A = zeros(n,n);
for i=1:n
    for j=1:n
        % 0- identify param index in DCM.Cp
        tmp = DCM.Vp;
        tmp.A(i,j) = 0;
        ind = find(full(VBA_spm_vec(tmp))~=full(VBA_spm_vec(DCM.Vp)));
        % 1- define reduced model
        E0r = E0;
        E0r(ind) = 0;
        V0r = V0;
        V0r(ind,:) = 0;
        V0r(:,ind) = 0;
        % 2- call Savage-Dickey ratios
        if isempty(ind) % this param was fixed
            if isequal(E0,E0r)
                dF = Inf;
            else
                dF = -Inf;
            end
        else
            dF = VBA_spm_log_evidence(E,V,E0,V0,E0r,V0r);
        end
        % 3- store Bayes factor
        sdr.A(i,j) = 1./(1+exp(dF));
    end
end

sdr.B = zeros(n,n,nu);
for i=1:n
    for j=1:n
        for k=1:nu
            tmp = DCM.Vp;
            tmp.B(i,j,k) = 0;
            ind = find(full(VBA_spm_vec(tmp))~=full(VBA_spm_vec(DCM.Vp)));
            E0r = E0;
            E0r(ind) = 0;
            V0r = V0;
            V0r(ind,:) = 0;
            V0r(:,ind) = 0;
            if isempty(ind) % this param was fixed
                if isequal(E0,E0r)
                    dF = Inf;
                else
                    dF = -Inf;
                end
            else
                dF = VBA_spm_log_evidence(E,V,E0,V0,E0r,V0r);
            end
            sdr.B(i,j,k) = 1./(1+exp(dF));
        end
    end
end

sdr.C = zeros(n,nu);
for i=1:n
    for j=1:nu
        tmp = DCM.Vp;
        tmp.C(i,j) = 0;
        ind = find(full(VBA_spm_vec(tmp))~=full(VBA_spm_vec(DCM.Vp)));
        E0r = E0;
        E0r(ind) = 0;
        V0r = V0;
        V0r(ind,:) = 0;
        V0r(:,ind) = 0;
        if isempty(ind) % this param was fixed
            if isequal(E0,E0r)
                dF = Inf;
            else
                dF = -Inf;
            end
        else
            dF = VBA_spm_log_evidence(E,V,E0,V0,E0r,V0r);
        end
        sdr.C(i,j) = 1./(1+exp(dF));
    end
end

nD = size(DCM.Ep.D,3);
if nD > 0
    for i=1:n
        for j=1:n
            for k=1:n
                tmp = DCM.Vp;
                tmp.D(i,j,k) = 0;
                ind = find(full(VBA_spm_vec(tmp))~=full(VBA_spm_vec(DCM.Vp)));
                E0r = E0;
                E0r(ind) = 0;
                V0r = V0;
                V0r(ind,:) = 0;
                V0r(:,ind) = 0;
                if isempty(ind) % this param was fixed
                    if isequal(E0,E0r)
                        dF = Inf;
                    else
                        dF = -Inf;
                    end
                else
                    dF = VBA_spm_log_evidence(E,V,E0,V0,E0r,V0r);
                end
                sdr.D(i,j,k) = 1./(1+exp(dF));
            end
        end
    end
end

DCM.sdr = sdr;

author: MBB-team

repo: VBA-toolbox


text:
function [n,edges,nbedges,xechan] = hist_ic(x,crit)

%HIST_IC  optimal Histogram based on IC information criterion
%
%   [N,EDGES,NBEDGES,XECHAN] = HIST-IC(X,CRIT)
%	bins the elements of X into an optimal number of bins according
%	to a cost function based on Akaike's Criterion.
%
%
%   CRIT = 1 | 2 | 3  (choose one of the 3 possible criterium)  (default 3)
%          4          (returns the initial histogram instead of the optimal one)
%
%
%   N = cell array containing the distribution of each column of X
%   (or a vector if X is a column vector)
%   EDGES = cell array containing the bin edges of each column of X
%   (or a vector if X is a column vector)
%   NBEDGES = vector containing the number of bin edges for each column of X
%   (or a number if X is a column vector)
%   XECHAN = discretized version of X
%
%   Ref : O. Colot et al., Information Criteria and Abrupt Changes in
%         Probability Laws, Signal Processing VII: Theory and Applications
%	  pp.1855-18858, September 1994
%
%   F. El-Matouat, O. Colot 2000 (first version)
%   Revised 01-06-2001 by Ph. Leray - philippe.leray@univ-nantes.fr
%
%
%   Things to do :
%	* Call criteron by a name ('aic','xxx', ...) instead of a number
%

if nargin == 0
    error('Requires one or two arguments.')
end

if nargin == 1
    crit = 3;
end;

if min(size(x))==1, x = x(:); end

if isstr(x)
    error('Input argument must be numeric.')
end

if isempty(x),
	error('No elements to count')
end

[nb_l,nb_c]=size(x);

% Outputs declaration
xechan=zeros(nb_l,nb_c);

edges=cell(nb_c,1);
% Local variables
maxi = max(x);
mini = min(x);

%% Erreur ? ancien code :
%% nb_clas_ini=2*round(sqrt(nb_l)-1);	% article Fatima

nb_clas_ini=round(2*sqrt(nb_l)-1);

pas_ini=(maxi-mini)/nb_clas_ini;	% initial step

for j=1:nb_c,

	% optimal histogram for each column of X
	histo_ini =hist(x(:,j),nb_clas_ini);		% initial histogram

	if (crit~=4)
		[hist_opt,pas_opt]=hist1_ic(histo_ini,nb_l,pas_ini(j),nb_clas_ini,crit);
	else
		fprintf('Histo initial\n');
		hist_opt=histo_ini;
		pas_opt=ones(1,nb_clas_ini)*pas_ini(j);
	end;
	nbedges(j)=size(hist_opt,2);
	edges{j}=mini(j)+cumsum(pas_opt(1:nbedges(j)-1)); %+1e-7;
        [n{j} xechan(:,j)]=histc(x(:,j),[-inf edges{j} inf]);
	n{j}=n{j}(1:end-1);
end

if (nb_c==1)
	n=n{1}; edges=edges{1};
end

% ============================== subfunctions

function [hist_opt,step_opt]=hist1_ic(histo,nb,step_ini,m,critere);

%HIST1_IC  optimal Histogram based on IC information criterion
%
%   [HIST_OPT, STEP_OPT] = HIST1_IC(HISTO, NB, STEP_INI, NBSTEP_INI, CRIT)
%       fusion of an 1D histogramme (HISTO) according to an IC criterion (CRIT)
%
%   This function is mainly an internal function used by HIST_IC
%
%   Ref : O. Colot et al., Information Criteria and Abrupt Changes in
%         Probability Laws, Signal Processing VII: Theory and Applications
%         pp.1855-18858, September 1994
%
%   F. El-Matouat, O. Colot 2000 (first version)
%   Revised 11-06-2001 by Ph. Leray
%
%
%   Things to do :
%       * Call criteron by a name ('aic','xxx', ...) instead of a number
%

aic=[];
aic2=[];

% Initialisation
histt=histo;
teta = histt/nb;
pas = step_ini*ones(1,m);

% Calcul de l'ensemble des histogrammes optimaux

for z=1:m

	% Calcul de AIC pour l'union entre hist(indice,u) et hist(indice,u+1)
	aic2 = [aic2 cal_aic(nb,teta,pas,m+1-z,critere)];

	if (z~=m)
		% Calcul des couples de classes adjacentes
		if critere==1
			penalite=(2*(m-z)-1)/nb;
		elseif critere==2
			penalite=(m-z-1)*(1+log(nb))/nb;
		else
			penalite=(m-z)*(1+log(log(nb)))/nb;
		end
 		aic=cla_adj(nb,teta,pas,step_ini,m-z+1,histt,penalite,aic);
		% Recherche de la valeur min du crit�re pour les classes adjacentes
		[min_aic classe]=min(aic(1:(m-z)));

		% Fusion de hist(classe) et hist(classe+1)
		nb_pas1=pas(classe)/step_ini;
		nb_pas2=pas(classe+1)/step_ini;

		ess=round( nb_pas1*histt(classe)+nb_pas2*histt(classe+1) );
		teta(classe)=ess / nb;
		histt(classe)=ess / (nb_pas1+nb_pas2);
		pas(classe)=pas(classe)+pas(classe+1);

		% Cr�ation du nouvel histogramme
		itemp = setdiff(1:m+1-z,classe+1);
		histt = histt(itemp);
		pas = pas(itemp);
		teta = teta(itemp);
	end
end

% Recherche du crit�re minimun AIC
[min_AIC fusion]=min(aic2(1:m));

% Initialisation de histo
histt=histo;
teta = histt/nb;
pas = step_ini*ones(1,m);

% Calcul de l'histogramme optimal

for z=1:fusion-1

	% Calcul des couples de classes adjacentes
	if critere==1
		penalite=(2*(m-1)-1)/nb;
	elseif critere==2
		penalite=(m-2)*(1+log(nb))/nb;
	else
		penalite=(m-1)*(1+log(log(nb)))/nb;
	end

	aic=cla_adj(nb,teta,pas,step_ini,m,histt,penalite,aic);

	% Recherche de la valeur min du crit�re pour les classes adjacentes
	[min_aic classe]=min(aic(1:m-1));

	% Fusion de hist(indice,classe) et hist(indice,classe+1)

	nb_pas1=pas(classe)/step_ini;
	nb_pas2=pas(classe+1)/step_ini;

	teta(classe)=(round(nb_pas1*histt(classe)+nb_pas2*histt(classe+1)))/nb;
	histt(classe)=(nb_pas1*histt(classe)+nb_pas2*histt(classe+1))/(nb_pas1+nb_pas2);
	pas(classe)=pas(classe)+pas(classe+1);

	% Cr�ation du nouvel histogramme

	itemp=setdiff(1:m,classe+1);
	histt = histt(itemp);
	pas = pas(itemp);
	teta = teta(itemp);
	%aic=zeros(1,m-1);

	m=m-1;
end
hist_opt=histt;
step_opt=pas;

%=====================================================
% Calcul du Critere pour l'ensemble des classes

function akaike=cal_aic(size_ech,teta,pas,m,critere);

if critere==1
	a=(2*m-1)/size_ech;
elseif critere==2
	a=(m-1)*(1+log(size_ech))/size_ech;
else
	a=m*(1+log(log(size_ech)))/size_ech;
end

indu = find(teta);
akaike = a - 2*sum(teta(indu).*log(teta(indu)./pas(indu)));

%=====================================================
% Cla_adj.m
% aic=cla_adj(taille,indice,teta,pas,pas_ini,m,hist,penalite,aic)
% taille=nombre d'�l�ments dans chacune des hypotheses;
% indice=numero de la classe;
% Calcul du critere de Akaike pour l'histogramme totale avec
% fusion de deux classes adjacentes u et (u+1).

function aic=cla_adj(size_ech,teta,pas,pas_ini,m,hist,penalite,aic);

for u=1:m-1

	cumul=0;

	% This loop is faster than a sum of a vectorised computation !
	for x=1:m
		if x~=u & x~=u+1 & teta(x)~=0
			cumul=cumul+teta(x)*log(teta(x)/pas(x));
		end
	end

	nb_pas1=pas(u)/pas_ini;
	nb_pas2=pas(u+1)/pas_ini;

	b=( round(nb_pas1*hist(u)+nb_pas2*hist(u+1) ) ) / size_ech;

	if b~=0
		c=2*b*log( b / ( pas(u) + pas(u+1) ) );
	else
		c=0;
	end

	aic(u)=penalite-2*cumul-c;
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% GetValueOfAssignment Gets the value of a variable assignment in a factor.
%
%   v = GetValueOfAssignment(F, A) returns the value of a variable assignment,
%   A, in factor F. The order of the variables in A are assumed to be the
%   same as the order in F.var.
%
%   v = GetValueOfAssignment(F, A, VO) gets the value of a variable assignment,
%   A, in factor F. The order of the variables in A are given by the vector VO.
%
%   See also SetValueOfAssignment.m and SampleFactors.m
%
% Copyright (C) Daphne Koller, Stanford University, 2012

function v = GetValueOfAssignment(F, A, VO)

if (nargin == 2),
    indx = AssignmentToIndex(A, F.card);
else
    map = zeros(length(F.var), 1);
    for i = 1:length(F.var),
        map(i) = find(VO == F.var(i));
    end;
    indx = AssignmentToIndex(A(map), F.card);
end;

v = F.val(indx);

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function []=addLightButton(varargin)
%% Viewing figures: projecting light
% addLightButton
% addLightButton(hf)

% This script creates a new pushtool in a figure's toolbar to manipulate
% the light projected on a displayed axis
%%
switch nargin
    case 0
        hf=gcf;
    case 1
        hf=varargin{1};
end

% Get icon
filePath=mfilename('fullpath');
toolboxPath=fileparts(fileparts(filePath));
iconPath=fullfile(toolboxPath,'lib_ext','GIBBON','icons');
hb = findall(hf,'Type','uitoolbar');

D=imread(fullfile(iconPath,'lightbulb.jpg'));
S=double(D);
S=S-min(S(:));
S=S./max(S(:));
if size(S,3)==1
    S=repmat(S,[1 1 3]);
end

% Create a uipushtool in the toolbar
uipushtool(hb(1),'TooltipString','light','CData',S,'Tag','light_button','ClickedCallback',{@lightFunc,{hf}});

end

%% Light function  lightFunc

function lightFunc(~,~,inputCell)

hf = inputCell{1};
ha = findobj(hf,'type','axes');

if isempty(ha)
    msgbox('There are no axes in the figure');
else
    prompt = {'Light'};
    title = '';
    format = struct('type','list');
    format.style = 'popupmenu';
    format.items = {' -- Select a lighting option -- ','none','left','right','headlight'};
    default = cell(size(prompt,1),1);
    default{1,1} = 1; % Default selection will always be instruction item
    default = cell2struct(default,prompt,1);
    prompt = repmat(prompt,1,2);
    options.AlignControls = 'on';
    choice = inputsdlg(prompt,title,format,default,options);

    % Unless the user closes the pop up window, clear all existing light
    if choice.Light ~= 1
        existingLight = findobj(ha,'type','light');
        if ~isempty(existingLight)
            delete(existingLight);
        end
        % Apply selected light to all figure axes
        if choice.Light == 3
            for ia = 1:size(ha,1)
                axes(ha(ia)); % Set the current axis
                h(ia) = camlight('left');
            end
        elseif choice.Light == 4
            for ia = 1:size(ha,1)
                axes(ha(ia));
                h(ia) = camlight('right');
            end
        elseif choice.Light == 5
            for ia = 1:size(ha,1)
                axes(ha(ia));
                h(ia) = camlight('headlight');
            end
        end
    end
end
end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% Modified by Rana Odabas 2018
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function p=legendrep(m,x)
%   function which construct Legendre polynomial Pm(x)
%   where M is the degree of polynomial and X is the variable.
%   Result - P is the char string that should be evaluated EVAL(P)
%   Example:
%   P=legendrep(2,'cos(theta)') will produce
%   P='(3*cos(theta).^2 -1)/2' which then can be evaluated as
%   theta=0.3; P=legendrep(2,'cos(theta)'); Lp=eval(P); produce
%   Lp = 0.8690
%   For Matlab R14 the following example can be used:
%   x=-5:.1:5; p=legendrep(5,'x .*cos(x)'); xp = eval(p);
%   figure; plot(x, xp, 'r.-'); grid

%% References:
%   Gradshteyn, Ryzhik "Table of Integrals Series and Products", 6th ed., p.973
%
%__________________________________________________
% 	Sergei Koptenko, Resonant Medical Inc., Toronto
%	sergei.koptenko@resonantmedical.com
%______________March/30/2004_______________________
%%
switch m
case 0
    p='1';     return
case 1
    p=x;     return
case 2
    p=['(3*' x '.^2 -1)/2'];     return
case 3
    p=['(5*' x '.^3 - 3 *' x ')/2'];     return
case 4
    p=['(35 *' x '.^4 - 30 * ' x '.^2 + 3)/8'];     return
case 5
    p=['( 63 * ' x '.^5 - 70 * ' x '.^3 + 15 *' x ' )/8'];     return
case 6
    p=['(231 *' x '.^6 - 315 * ' x '.^4 + 105 * ' x '.^2 -5)/16'];     return
 case 7
   p=['(429 * ' x '.^7 - 693 * ' x '.^5 +315 * ' x '.^3 -35 *' x ' )/16'];     return
case 8
    p=['(6435 *' x '.^8 -12012 *' x '.^6 + 6930 * ' x '.^4 -1260 * ' x '.^2 +35)/128'];     return
case 9
    p=['(12155 * ' x '.^9 -25740 * ' x '.^7 +18018 * ' x '.^5 -4620 * ' x '.^3 +315 *' x ' )/128'];
    return
case 10
       p=['(46189 *' x '.^10 -109395 *' x '.^8 +90090 *' x '.^6 -30030 * ' x '.^4 +3465 * ' x '.^2 -63)/256'];
    return
case 11
   p=['(88179 * ' x '.^11 -230945 * ' x '.^9 +218790 * ' x '.^7 -90090 * ' x '.^5 +15015 * ' x '.^3 -693 *' x ' )/256'];
    return
case 12
    p=['(676039 *' x '.^12 -1939938 *' x '.^10 +2078505 *' x '.^8 -1021020 *' x '.^6 +225225 * ' x '.^4 -18018 * ' x '.^2 +231)/1024'];
    return

otherwise
    iii=m-10;    %shift counter
   pp=strvcat(legendrep(11,x),legendrep(12,x)); % get last two members
    for ii=3:1:iii,   % Begin construct from 13th member
        p_ii=[num2str(1/(ii)) ' * (' num2str(2*ii-1) ' * ' x '.*(' ...
        deblank(pp(ii-1,:)) ')-' num2str(ii-1) '.*(' deblank(pp(ii-2,:)) '))'];
        pp=strvcat(pp, p_ii);
    end
p=deblank(pp(iii,:)); % remove traiing blanks
   return
end

return

author: Sable

repo: mcbench-benchmarks


text:
classdef (Abstract) ProbabilityDistributionX < BaseX
% ProbabilityDistributionX class
%
% Summary of ProbabilityDistributionX:
% This is the base class for all distribution objects.
%
% November 2018 Lyudmil Vladimirov, University of Liverpool.

    properties (Dependent, SetAccess = protected)
        % NumVariables - Number of random variables for this distribution
        NumVariables
    end

    properties (Access = protected)
        % NumVariables_ - Internal storage for number of random variables
        NumVariables_
    end

    methods (Abstract)
        % Random Draw random samples from the distribution
        %
        % Parameters
        % ----------
        % numSamples: scalar
        %   Number of samples to draw from distribution
        %
        samples = random(this, numSamples)

        % Reset Reset the distribution
        %
        % Parameters
        % ----------
        % numVariables: scalar
        %   Number of random variables that are described by the
        %   distribution
        reset(this, numVariables);
    end

    methods
        function this = ProbabilityDistributionX()
        %ProbabilityDistribution Construct a numVariables-variate probability distribution

        end

        function numVariables = get.NumVariables(this)
        %get.NumVariables Custom getter for NumVariables property
            numVariables = getNumVariables(this);
        end

        function set.NumVariables(this, numVariables)
        %set.NumVariables Custom setter for NumVariables property
            this.NumVariables_ = setNumVariables(this, numVariables);
        end
    end

    methods (Access=protected)
        function numVariables = getNumVariables(this)
            numVariables = this.NumVariables_;
        end
        function NumVariables = setNumVariables(this, numVariables)
            NumVariables = numVariables;
        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function [x,state] = struct_ctranspose(z,task)
%STRUCT_CTRANSPOSE Complex conjugate transpose.
%   [x,state] = struct_ctranspose(z) computes x as the complex conjugate
%   transpose of z. The structure state stores information which is reused
%   in computing the right and left Jacobian-vector products.
%
%   struct_ctranspose(z,task) computes the right or left Jacobian-vector
%   product of this transformation, depending on the structure task. Use
%   the structure state and add the field 'r' of the same shape as z or the
%   field 'l' of the same shape as x to obtain the structure task for
%   computing the right and left Jacobian-vector products
%
%      (dF(:)/dz(:).')*task.r(:) and
%      (dF(:)/dz(:).')'*task.l(:) + conj((dF(:)/dconj(z(:)).')'*task.l(:)),
%
%   respectively. Here, F(z) represents this transormation, (:) signifies
%   vectorization and the derivative w.r.t. z (conj(z)) is a partial
%   derivative which treats conj(z) (z) as constant. The output has the
%   same shape as x or z for the right and left Jacobian-vector products,
%   respectively.
%
%   See also struct_conj, struct_transpose.

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

if nargin < 2, task = []; end
state = [];

if isempty(task) || (isempty(task.l) && isempty(task.r))
    x = z';
elseif ~isempty(task.r)
    if ~isreal(z) || ~isreal(task.r)
        error('struct_ctranspose:nonanalytic',['Nonanalytic objective ' ...
            'functions are currently not supported in sdf_nls, please ' ...
            'use sdf_minf instead.']);
    end
    x = task.r.';
elseif ~isempty(task.l)
    x = task.l';
end

end

author: andrewssobral

repo: lrslibrary


text:
function x = sb_solve(sysmat,vecb)

% SB_SOLVE
%
% $Id$

%scalen
disp('Scaling stiffness matrix...')
dkond = 1./(sqrt(diag(sysmat)));
vecb = vecb.*dkond;
[indexi, indexj, s] = find(sysmat);
sys_size = size(sysmat,1);
clear sysmat;
s = (s.*dkond(indexi)).*dkond(indexj);
s(1) = 1;
disp('Preconditioning...')
L = sparse(indexi,indexj,s,sys_size,sys_size,length(s));
%partch
try
    L = ichol(L);
catch
    disp('Could not compute incomplete Cholesky-decompositon. Rescaling stiffness matrix...')
    alpha = 0.5d-6;
    alpha = alpha*8.d0;
    alpha = 1 / (alpha + 1);
    s = alpha*s;
    dia = find(indexi == indexj);
    s(dia) = (1./alpha)*s(dia);
    s(dia) = sqrt(s(dia));
    s(1) = 1;
    L = sparse(indexi,indexj,s,sys_size,sys_size,length(s));
    clear dia;
    L = ichol(L);
end
%startvektor
disp('Finding startvector...')
vecb_ = L \ (-vecb);
vecx = L' \ vecb_;
clear vecb_;
%sonstiges
sysmat = sparse(indexi,indexj,s,sys_size,sys_size,length(s));
sysmat = sysmat + sysmat' - sparse(1:sys_size,1:sys_size,diag(sysmat),sys_size,sys_size,sys_size);
clear indexi indexj s;
%dprod = sysmat * vecx;
%lösen
disp('Solving equation system...')
[x,fl,rr,it,rv]= pcg(sysmat,vecb,10e-9,5000,L,L',vecx);
%fl
rr
%it
%rescal
x = x.*dkond;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function jac = p13_jac2 ( option, m, nvar, lambda, u )

%*****************************************************************************80
%
%% P13_JAC2 computes the jacobian by recasting it on a square grid.
%
%  Discussion:
%
%    Actually, to stave off insanity, we only "recast" the variables into
%    a 2D array that corresponds to the spatial ordering of the grid.
%    We leave the jacobian in its original arrangement, which assumes
%    a linear ordering of variables and equations, and we simply
%    compute the equation and variable indices of the jacobian when
%    we are ready to put entries into it.  This approach seems to produce
%    a smaller amount of cosmic grief than the alternatives.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Input, integer M, the number of grid points on a side of the square.
%
%    Input, integer NVAR, the number of variables.
%
%    Input, real LAMBDA, the value of the parameter.
%
%    Input, real U(M,M), the value of the grid function.
%
%    Output, real JAC(NVAR-1,NVAR), the jacobian matrix evaluated
%    at X.  The NVAR-th row is not set by this routine.
%
  jac(1:nvar-1,1:nvar) = 0.0;

  h = 1.0 / ( m + 1 );

  ieqn = 0;

  for i = 1 : m
    for j = 1 : m

      ieqn = ( j - 1 ) * m + i;

      uc = u(i + (j-1)*m);

      if ( i < m )
        un = u(i+1+(j-1)*m);
      else
        un = 0.0;
      end

      if ( 1 < i )
        us = u(i-1+(j-1)*m);
      else
        us = 0.0;
      end

      if ( j < m )
        ue = u(i+j*m);
      else
        ue = 0.0;
      end

      if ( 1 < j )
        uw = u(i+(j-2)*m);
      else
        uw = 0.0;
      end

      fc = p13_gx ( option, uc );
      fn = p13_gx ( option, un );
      fs = p13_gx ( option, us );
      fe = p13_gx ( option, ue );
      fw = p13_gx ( option, uw );

      del5f = fc + h * h * ( - 4.0 * fc + fn + fs + fe + fw ) / 12.0;

      fcp = p13_gp ( option, uc );
      fnp = p13_gp ( option, un );
      fsp = p13_gp ( option, us );
      fep = p13_gp ( option, ue );
      fwp = p13_gp ( option, uw );

      ivar = ( j - 1 ) * m + i;
      jac(ieqn,ivar) = - 20.0 / ( 6.0 * h * h ) ...
        + lambda * ( fcp - 4.0 * h * h * fcp / 12.0 );

      if ( i < m )
        ivar = ( j - 1 ) * m + i + 1;
        jac(ieqn,ivar) = 4.0 / ( 6.0 * h * h ) ...
          + lambda * h * h * fnp / 12.0;
      end

      if ( 1 < i )
        ivar = ( j - 1 ) * m + i - 1;
        jac(ieqn,ivar) = 4.0 / ( 6.0 * h * h ) ...
          + lambda * h * h * fsp / 12.0;
      end

      if ( j < m )
        ivar = j * m + i;
        jac(ieqn,ivar) = 4.0 / ( 6.0 * h * h ) ...
          + lambda * h * h * fep / 12.0;
      end

      if ( 1 < j )
        ivar = ( j - 2 ) * m + i;
        jac(ieqn,ivar) = 4.0 / ( 6.0 * h * h ) ...
          + lambda * h * h * fwp / 12.0;
      end

      if ( 1 < i & 1 < j )
        ivar = ( j - 2 ) * m + i - 1;
        jac(ieqn,ivar) = 1.0 / ( 6.0 * h * h );
      end

      if ( 1 < i & j < m )
        ivar = j * m + i - 1;
        jac(ieqn,ivar) = 1.0 / ( 6.0 * h * h );
      end

      if ( i < m & 1 < j )
        ivar = ( j - 2 ) * m + i + 1;
        jac(ieqn,ivar) = 1.0 / ( 6.0 * h * h );
      end

      if ( i < m & j < m )
        ivar = j * m + i + 1;
        jac(ieqn,ivar) = 1.0 / ( 6.0 * h * h );
      end

      ivar = nvar;
      jac(ieqn,nvar) = del5f;

    end
  end

  return
end


text:
function y = SmoothRamp(x,alpha)

%y = SmoothRamp(x,alpha)
%
% This function is used to make a smooth version of the ramp function
%
%INPUTS:
% x = a vector of real numbers to be smoothed
% alpha = a positive smoothing parameter. Small values of alpha correspond
%         to little smoothing, large values correspond to heavy smoothing

y = alpha*log(exp(x/alpha)+1);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function writeDoublePendulumDynamics(f,M)

%This function writes the dynamics file for the double pendulum.

filename = 'doublePendulumDynamics';

comments{1} = ['DZ = ' upper(filename) '(T,Z,U,P)'];
comments{2} = ' ';
comments{3} = 'FUNCTION:  This function computes the dynamics of a double';
comments{4} = '    pendulum, and is designed to be called from ode45. The';
comments{5} = '    model allows for arbitrary mass and inertia for each';
comments{6} = '    link, but no friction or actuation';
comments{7} = ' ';
comments{8}  = 'INPUTS: ';
comments{9}  = '    t = time. Dummy input for ode45. Not used.';
comments{10} = '    z = [4xn] matrix of states.';
comments{11} = '    u = [2xn] matrix of inputs';
comments{12} = '    P = struct of parameters';
comments{13} = 'OUTPUTS: ';
comments{14} = '    dz = [4xn] matrix of state derivatives';
comments{15} = ' ';
comments{16} = 'NOTES:';
comments{17} = ['    This file was automatically generated by ' mfilename];

params{1} = {'m1','link one mass'};
params{2} = {'m2','link two mass'};
params{3} = {'g ','gravity'};
params{4} = {'l1','link one length'};
params{5} = {'l2','link two length'};
params{6} = {'I1','link one moment of inertia about its center of mass'};
params{7} = {'I2','link two moment of inertia about its center of mass'};
params{8} = {'d1','distance between link one center of mass and parent joint'};
params{9} = {'d2','distance between link two center of mass and parent joint'};

states{1} = {'th1','link one absolute angle'};
states{2} = {'dth1','link one angular rate'};
states{3} = {'th2','link two absolute angle'};
states{4} = {'dth2','link two angular rate'};

input{1} = {'u1','torque acting on link 1 wrt ground'};
input{2} = {'u2','torque acting on link 2 wrt ground'};

dstates{1} = {'dth1','derivative of link one absolute angle'};
dstates{2} = {'ddth1','derivative of link one angular rate'};
dstates{3} = {'dth2','derivative of link two absolute angle'};
dstates{4} = {'ddth2','derivative of link two angular rate'};

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                               write file                                %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
fid = fopen([filename '.m'],'w');

fprintf(fid, ['function dz = ' filename '(~,z,u,P) \n']);

for i=1:length(comments)
    fprintf(fid,['%%' comments{i} '\n']);
end
fprintf(fid,'\n');

for i=1:length(params)
    fprintf(fid,[params{i}{1} ' = P.' params{i}{1} '; %%' params{i}{2} '\n']);
end
fprintf(fid,'\n');

for i=1:length(states)
    fprintf(fid,[states{i}{1} ' = z(' num2str(i) ',:); %%' states{i}{2} '\n']);
end
fprintf(fid,'\n');

for i=1:length(input)
    fprintf(fid,[input{i}{1} ' = u(' num2str(i) ',:); %%' input{i}{2} '\n']);
end
fprintf(fid,'\n');

for i=1:2;
    fprintf(fid,['f' num2str(i) ' = ' vectorize(char(f(i))) ';\n']);
end
fprintf(fid,'\n');

for i=1:2
    for j=1:2
        fprintf(fid,['M' num2str(i) num2str(j) ' = ' vectorize(char(M(i,j))) ';\n']);
    end
end
fprintf(fid,'\n');

fprintf(fid,'D = M11.*M22 - M12.*M21;\n\n');

fprintf(fid,'ddth1 = (f2.*M12 - f1.*M22)./D;\n');
fprintf(fid,'ddth2 = -(f2.*M11 - f1.*M21)./D;\n');
fprintf(fid,'\n');

fprintf(fid,'dz = [...\n');
for i=1:length(dstates)
    fprintf(fid,['    ' dstates{i}{1} '; %%' dstates{i}{2} '\n']);
end
fprintf(fid,'];\n');

fprintf(fid,'end \n');

fclose(fid);
end


text:
function pred_value = xval_simple_ols_loo(X, Y)
% Check: a very simple leave-one-out cross-validated regression
% pred_value = xval_simple_ols_loo(X, Y)
%
% pred_value: cross-validated predictions of outcome data
% X: n x variables matrix of predictors
% Y: n x 1 vector of outcomes
%
% Tor Wager, June 2010
%
% Go to any LASSO output directory and run this:
%
% maskInfo = iimg_read_img(fullfile(pwd, 'mask.img'), 2);
% dat{1} =  iimg_get_data(maskInfo, anticimages);
% pred_value = xval_simple_ols_loo(X, Y)

[N, k] = size(X);
pred_value = zeros(N, 1);

create_figure('test', 1, 2)

for i = 1:N

    % select training data
    Xi = X; Xi(i, :) = [];
    Yi = Y; Yi(i) = [];

    Xi = [ones(N-1, 1) Xi]; % add intercept

    % Make prediction
    b = pinv(Xi) * Yi;
    pred_value(i, 1) = [1 X(i, :)] * b;
     %pred_value(i, 1) = X(i, :) * b;

    create_figure('test', 1, 2, 1); subplot(1, 2, 1); plot(Xi*b, Yi, 'ko');
    plot(pred_value(i, 1), Y(i), 'ro', 'MarkerFaceColor', 'r');
    subplot(1, 2, 2);
   plot([ones(N, 1) X] * b, 'k');
   title('Black circles = training set; Red = holdout obs');
    %plot(X * b, 'k');
    drawnow;
    %pause(.1);

     %intercept_vals(i, 1) = b(1);
end

cm = colormap(jet(N));
figure; hold on;
for i = 1:N
    plot(pred_value(i), Y(i), 'ko', 'MarkerFaceColor', cm(i, :));
end
xlabel('Predicted outcome (xval)'); ylabel('Outcome');
title('Color = order in data series');

% figure; hold on;
% for i = 1:N
%     plot(intercept_vals(i), pred_value(i), 'ko', 'MarkerFaceColor', cm(i, :));
% end

%
% figure; hold on;
%     plot(pred_value(whorder(1:24)), Y(whorder(1:24)), 'bo', 'MarkerFaceColor', 'b');
% plot(pred_value(whorder(25:end)), Y(whorder(25:end)), 'ro', 'MarkerFaceColor', 'r');

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function images = get_stimulusimages(username, accessKey)
% get_stimulusimages    Get stimulus images from ImageNet
%
% Note:
%
% This script is for downloading all 1200 training and 50 test images of
% experimental stimulus used in Horikawa and Kamitani 2017
% "Generic decoding of seen and imagined objects using hierarchical visual features".
% It will take about one day to complete all processes with single computer,
% but you can accerelate by using multiple computers in parallel.
%
% This function requires `websave`, which is equipped in MATLAB >= 2014b.
%
% Inputs:
%
% - username  [char] : ImageNet user name
% - accessKey [char] : ImageNet access key
%
% Outputs:
%
% - images [struct] : A structure contains image information
%
% Settings:
%
% - imageListFiles    : List of CSV files which contains image IDs and image file names
% - imageType         : List of strings representing image types
% - imageDir          : Directory in which images will be saved
% - imageSize         : The images will be trimmed down with this size (pixel)
% - deleteOriginalImg : If true, remove the downloaded original image files
%

%% Settings %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% List of CSV files which contains image IDs and image file names
imageListFiles = {'./imageID_training.csv', ...
                  './imageID_test.csv'};

% Type of images ('training' or 'test')
% This should correspond to 'imageListFiles'
imageType = {'training', ...
             'test'};

% The directory in which images will be saved
imageDir = './images';

% Image trimming size (500 px)
imageSize = 500;

% Delete the downloaded original images or not
deleteOriginalImg = false;

% Misc options for downloading (websave)
options = weboptions;
options.Timeout = 100;

%% Main %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf('%s started\n', mfilename);

urlImageNet = 'http://image-net.org/download/synset';

%% Check inputs ----------------------------------------------------------------
if ~exist('username','var')
    error('Please give your ImageNet username name.');
end

if ~exist('accessKey','var')
    error('Please give your ImageNet access key.');
end

%% Get image lists -------------------------------------------------------------
for i = 1:length(imageListFiles)
    fid = fopen(imageListFiles{i}, 'r');
    imageList = textscan(fid, '%n%s', 'delimiter', ',');
    fclose(fid);

    % Remove quotations from file names
    for j = 1:length(imageList{2})
        imageList{2}{j} = imageList{2}{j}(2:end-1);
    end

    images(i).type     = imageType{i};
    images(i).listFile = imageListFiles{1};
    images(i).id       = imageList{1};
    images(i).files    = imageList{2};

    %% Get image files for each synset
    images(i).synsets = [];
    synsetList = {};
    for j = 1:length(images(i).files)
        imgFileName = images(i).files{j};
        synsetId = regexp(imgFileName, '^(n\d+)_\d+\.JPEG$', 'tokens');
        synsetId = synsetId{1}{1};

        [ismem, ind] = ismember(synsetId, synsetList);

        if ismem
            images(i).synsets(ind).files{end+1} = imgFileName;
        else
            images(i).synsets(end+1).name = synsetId;
            images(i).synsets(end).files{1} = imgFileName;

            synsetList{end+1} = synsetId;
        end
    end

    images(i).synsetList = unique(synsetList);
end

%% Download and convert images -------------------------------------------------
for i = 1:length(images)
for j = 1:length(images(i).synsets)

    %% Multiprocessing
    imgType = images(i).type;
    synsetId = images(i).synsets(j).name;

    lockFile = sprintf('%s_%s_%s.lock', mfilename, imgType, synsetId);
    if exist(lockFile,'file')
        continue;
    end
    fclose(fopen(lockFile, 'w'));

    fprintf('Image type:\t%s\n', imgType);
    fprintf('Synset:\t\t%s\n', synsetId);

    %% Prepare directories
    synsetDir = fullfile(imageDir, imgType, synsetId);
    if ~exist(synsetDir, 'dir')
        fprintf('Create %s\n', synsetDir);
        mkdir(synsetDir);
    end

    dlDir = fullfile(imageDir, 'download', synsetId);
    dlFile = fullfile(dlDir, [synsetId '.tar']);

    if exist(dlFile,'file');
        fprintf('%s already exists. Download was skipped.\n', dlFile);
    else
        if ~exist(dlDir, 'dir')
            fprintf('Create %s\n', dlDir);
            mkdir(dlDir);
        end

        %% Download synset images
        fprintf('Downloading %s...', synsetId);
        tic;
        res = websave(dlFile, urlImageNet, ...
                      'wnid', synsetId, ...
                      'username', username, ...
                      'accesskey', accessKey, ...
                      'release', 'latest', ...
                      'src', 'stanford', ...
                      options);
        elapsedTime = toc;
        fprintf('done (%f sec)\n', elapsedTime);

        fprintf('Extracting image files...');
        untar(dlFile, dlDir)
        fprintf('done\n');
    end

    %% Copy and trim images
    for k = 1:length(images(i).synsets(j).files)
        imgFile = images(i).synsets(j).files{k};
        sourceFile = fullfile(imageDir, 'download', synsetId, imgFile);
        targetFile = fullfile(synsetDir, imgFile);

        if exist(targetFile,'file')
            continue;
        end

        fprintf('Converting %s\n', imgFile);

        while ~exist(sourceFile,'file')
            % Wait until 'sourceFile' is created by another process
            pause(1);
        end

        img = imread(sourceFile);

        imgTrim = trimImgCenter(img);
        imgResize = imresize(imgTrim, [imageSize, imageSize]);

        imwrite(imgResize, targetFile);

        fprintf('%s saved\n', targetFile);
    end

    delete(lockFile);
end
end

%% Remove the downloaded original images
if deleteOriginalImg
    rmdir(fullfile(imageDir, 'download'), 's');
end

fprintf('Finished\n');

end

%% Subfunctions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [trimmed] = trimImgCenter(img)
% trimimgCenter    Trims an image on the center
%
% Input:
%
% - img: image
%
% Output:
%
% - trimmed: trimmed image

[r, c, d] = size(img);
trimSize = min([r, c]);
lr = floor((c - trimSize) / 2);
ud = floor((r - trimSize) / 2);

trimmed = img((1:min([(trimSize), r - ud])) + ud, ...
              (1:min([(trimSize), c - lr])) + lr, :);

end

author: KamitaniLab

repo: GenericObjectDecoding


text:
classdef chebop2
%CHEBOP2   CHEBOP2 class for representing linear partial differential operators.
%
% Class used to solve linear PDEs defined on rectangular domains that have unique and
% globally smooth solutions.
%
% N = CHEBOP2(@(u) op(u)) constructs an operator N representing the operator
% given in @(u)op(u) acting on functions of two variables on [-1,1] by [-1,1].
%
% N = CHEBOP2(@(u) op(u), [a b c d]) constructs an operator N acting on
% functions of two variables defined on [a,b] by [c,d].
%
% N = CHEBOP2(@(x,y,u) op(x,y,u),...) constructs a variable coefficient PDE
% operator. If a partial differential operator is constant coefficient, we
% recommend the @(u) notation rather than @(x,y,u) as it is more efficient.
%
% Boundary conditions are imposed via the syntax N.lbc, N.rbc, N.ubc, and N.dbc.
%
% Example 1: (Poisson with zero Dirichlet conditions):
%    N = chebop2(@(u) diff(u,2,1) + diff(u,2,2));
%    N.bc = 0;
%    u = N \ 1;
%
% Example 2: (Helmholtz equation with gravity)
%    N = chebop2(@(x,y,u) laplacian(u) - 10*y.^2.*u, [-1 1 -3 0]);
%    N.bc = 1;
%    u = N \ 0;
%
% Example 3: (Klein-Gordon equation)
%    N = chebop2(@(u) diff(u,2,1) - diff(u,2,2) + 5*u,[-1 1 0 3]);
%    N.lbc = 0; N.rbc = 0;
%    N.dbc = @(x,u) [u - exp(-30*x.^2) ; diff(u)];
%    u = N \ 0;
%
% Example 4: (Poisson equation with general Dirichlet conditions and rhs)
%    N = chebop2( @(u) lap(u));
%    N.lbc = @(y) -y.^2;
%    N.rbc = @(y) y.^2;
%    N.ubc = @(x) x;
%    N.dbc = @(x) x;
%    u = N \ chebfun2(@(x,y) cos(x.*y));
%
% For further details about the PDE solver, see:
%
% A. Townsend and S. Olver, The automatic solution of partial differential
% equations using a global spectral method, J. Comput. Phys., 299 (2015),
% pp. 106-123.
%
% Warning: This PDE solver is an experimental new feature. It has not been
% publicly advertised.  Chebop2 cannot do nonlinear problems as more
% algorithmic advances are needed.

% Copyright 2017 by The University of Oxford and The Chebfun2 Developers.
% See http://www.chebfun.org/ for Chebfun information.

    %% PROPERTIES.
    properties ( GetAccess = 'public', SetAccess = 'public' )

        domain = [];  % Domain of the operator.
        op = [];      % The operator.
        ubc = [];     % Up boundary condition(s).
        lbc = [];     % Left boundary condition(s).
        rbc = [];     % Right boundary condition(s).
        dbc = [];     % Down boundary condition(s).
        dim = [];     % Size of the system (number of eqns).
        scale = [];   % Relative solution scale.
        coeffs = [];  % Matrix storing constant coefficients.
        rhs = [];     % Righthand side, if given by user.
        xorder = 0;   % Diff order in the x-variable.
        yorder = 0;   % Diff order in the y-variable.
        U             %
        S             % Low rank form of the partial differential operator.
        V             %

    end

    %% CONSTRUCTOR.
    methods

        function N = chebop2(varargin)
            % CHEBOP2 CONSTRUCTOR.

            % Get CHEBFUN2 preferences.
            pref = chebfunpref();
            tol = pref.cheb2Prefs.chebfun2eps;

            % If empty input arguments then return an empty CHEBOP2 object.
            if ( isempty(varargin) )
                return
            end

            % What domain is the operator defined on?
            if ( numel(varargin) > 1 )
                ends = varargin{2}; % Second argument should be a domain.
                if ( length(ends) == 4 )
                    % Valid domain?
                    if ( diff(ends(1:2)) > 0 && diff(ends(3:4)) > 0 )
                        dom = ends;
                    else
                        error('CHEBFUN:CHEBOP2:chebop2:emptyDomain', ...
                            'Empty domain.');
                    end
                else
                    error('CHEBFUN:CHEBOP2:chebop2:badDomain',...
                        'Argument should be a domain given by four doubles.')
                end
            else
                if ( isa(varargin{1}, 'function_handle') )
                    % Pick the default domain.
                    rect1 = [-1, 1];
                    rect2 = [-1, 1];
                    dom = [rect1, rect2];
                elseif ( isa( varargin{1}, 'double') )
                    % Set up identity operator on the domain.
                    N = chebop2(@(u) u, varargin{1});
                    return
                else
                    error('CHEBFUN:CHEBOP2:chebop2:badArg',...
                        'First argument is not an operator or domain.')
                end
            end

            % First argument in the constructor is the operator. If the
            % operator is univariate then it's a constant coefficient PDE,
            % otherwise assume it is a variable coefficient.
            if ( isa(varargin{1}, 'function_handle') )
                fh = varargin{1};

                if ( nargin(fh) == 1 )  % The PDE has constant coefficients.
                    % Trust that the user has formed the CHEBFUN2 objects
                    % outside of CHEBOP2.

                    % Extract out rhs:
                    x = chebfun2(@(x,y) x, dom);
                    RHS = fh(0*x);
                    % In this case, the RHS must be a constant chebfun2.
                    if ( norm( RHS ) > 0 )
                        fh = @(u) fh(u) - RHS;
                        N.rhs = -RHS;             % store for later.
                    else
                        N.rhs = 0;
                    end

                    u = adchebfun2(chebfun2(@(x,y) x.*y, dom));
                    v = fh(u);
                    % If the PDO has constant coefficients then convert to
                    % double:
                    try
                        A = cell2mat(v.jacobian).';
                    catch
                        % PDO has variable coefficients, keep them in a
                        % cell array:
                        A = v.jacobian;
                    end

                elseif ( nargin(fh) == 2 )
                    error('CHEBFUN:CHEBOP2:chebop2:badOp1', ...
                        'Did you intend to have @(x,y,u)?')
                elseif ( nargin(fh) == 3 )
                    % The coefficients of the PDE are now variable
                    % coefficient.

                    % Setup a chebfun2 on the right domain
                    u = adchebfun2(chebfun2(@(x,y) x.*y, dom));
                    x = chebfun2(@(x,y) x, dom);
                    y = chebfun2(@(x,y) y, dom);

                    % Extract out rhs:
                    RHS = fh(x, y, 0*x);
                    % The RHS is a chebfun2. Check if it is the zero. Do
                    % not change function handle unless we have to because
                    % we want the user to be able to see the PDO
                    % untouched if possible.
                    if ( norm( RHS ) > 0 )
                        fh = @(x, y, u) fh(x, y, u) - RHS;
                        N.rhs = -RHS;             % store for later.
                    else
                        N.rhs = 0;
                    end

                    % Apply it to the operator.
                    v = fh(x, y, u);
                    A = v.jacobian;  % Cell array of variable coefficients.

                    % If we have a variable coefficient PDO, then compute the
                    % separable representation immediately. We need it now.
                    [cellU, matS, cellV] = chebop2.separableFormat( A,...
                                                    size(A,2), size(A,1), dom );
                    N.U = cellU;
                    N.S = matS;
                    N.V = cellV;
                else
                    error('CHEBFUN:CHEBOP2:chebop2:badOp2',...
                        'Operator should be @(u) or @(x,y,u).')
                end

            else
                error('CHEBFUN:CHEBOP2:chebop2:badOp3',...
                    'First argument should be an operator')
            end

            % Often the coefficients are obtained with small rounding errors
            % and it is important to remove the very small non-zero ones to
            % have rank(A) correct.
            if ( iscell(A) )
                for jj = size(A, 1)
                    for kk = size(A, 2)
                        if ( isa(A{jj,kk}, 'double') && abs(A{jj,kk}) < 10*tol )
                            A{jj,kk} = 0;
                        end
                    end
                end
            else
                A(abs(A) < 10*tol) = 0;
            end

            % Construct CHEBOP2 object. The boundary conditions will be
            % given later.
            N.domain = dom;
            N.op = fh;
            N.coeffs = A;

            % Calculate xorder and yorder of PDE.
            % Find the differential order of the PDE operator.
            if ( iscell(A) )
                xdifforder = size(A, 2) - 1;
                ydifforder = size(A, 1) - 1;
            elseif ( min(size(A)) > 1 )
                xdifforder = find(sum(abs(A), 2) > 100*tol, 1, 'last') - 1;
                ydifforder = find(sum(abs(A)) > 100*tol, 1, 'last' ) - 1;
            else
                if ( size(A, 1) == 1 )
                    ydifforder = length(A) - 1;
                    xdifforder = 0;
                else
                    xdifforder = length(A) - 1;
                    ydifforder = 0;
                end
            end
            N.xorder = xdifforder;
            N.yorder = ydifforder;

            % Issue a warning to the user for the first CHEBOP2:
            warning('CHEBFUN:CHEBOP2:chebop2:experimental',...
                ['CHEBOP2 is a new experimental feature.\n'...
                'It has not been tested to the same extent as other'...
                'parts of the software.']);
            % Turn it off:
            warning('off', 'CHEBFUN:CHEBOP2:chebop2:experimental');

        end

    end

    %% STATIC HIDDEN METHODS.
    methods ( Static = true, Hidden = true )

        % Matrix equation solver for AX - XB = F (p and q are ADI shifts):
        X = adi( A, B, F, p, q );

        % Matrix equation solver for A*X-X*B = M*N.' (p and q are ADI shifts):
        [UX, DX, VX] = fadi( A, B, M, N, p, q );

        % ADI shift parameters for ADI method:
        [p, q] = adiShifts(a, b, c, d, tol);

        % Matrix equation solver for AXB^T + CXD^T = E. xsplit, ysplit = 1 if
        % the even and odd modes (coefficients) decouple.
        X = bartelsStewart(A, B, C, D, E, xsplit, ysplit);

        % Use automatic differentiation to pull out the coeffs of the
        % operator:
        deriv = chebfun2deriv(op);

        % This is used to discretize the linear constrains:
        [bcrow, bcvalue] = constructBC(bcArg, bcpos,...
            een, bcn, dom, scl, order);

        % Recover coefficient functions of a linear operator:
        p = recoverCoeffs(L);

        % Checks the type of the boundary conditions:
        [bctype, g] = checkBC(N, m, n);

        % This is used to discretize the PDE:
        [CC, rhs, bb, gg, Px, Py, xsplit, ysplit] =...
            discretize(N, f, m, n, flag);

        % Convert all the different user inputs for bc into uniform format:
        bc = createBC(bcArg, ends);

        % Method for deciding how to solve the matrix equation:
        X = denseSolve(N, f, m, n);

        % Compute the separable representation of a PDO:
        [cellU, S, cellV] = separableFormat(A, xorder, yorder, dom);

        % Setup Laplace operator on domain DOM:
        N = setupLaplace( dom );

        % Remove trailing coefficients.
        a = truncate(a, tol);

    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
#computes the pose vector v from an homogeneous transform A
function v=t2v(A)
  v = [A(1:2,3); atan2(A(2,1),A(1,1))];
end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
% Fig. 2.3  Feedback Control of Dynamic Systems, 6e
%            Franklin, Powell, Emami
%

clear all;
close all;
clf
hold off
num = 1/1000;
den = [1 50/1000];
sys = tf(num,den);
t = 0:100;
y = step(num*500,den,t);
plot(t,y),grid
xlabel('Time (sec)')
ylabel('Amplitude')
title('Fig. 2.3')
nicegrid


text:
function ho=comprose(x,y,n,w,az,varargin)
%COMPROSE Compass rose plot
%
%       COMPROSE(X,Y,N,W,AZ) adds a compass rose on current axis located at
%       position X,Y with N points (N is 1, 4, 8 or 16), width W (radius)
%       and North pointing to azimuth AZ (in degree, AZ = 0 means an arrow
%       pointing to positive Y-axis direction, rotating clockwise).
%
%       COMPROSE(...,FS) adds Cardinal directions with font size FS.
%
%       COMPROSE(...,'param1',value1,'param2',value2,...) specifies any
%       additionnal properties of the Patch using standard parameter/value
%       pairs. Note that 'FaceColor' concerns only the default black-filled
%       parts of the drawing.
%
%	H=COMPROSE(...) returns a Nx3 matrix of object handles: first column
%	addresses solid filled patches, second column for white patches, third
%	column for Cardinal direction text.
%
%       Examples:
%		comprose(0,0,8,.5,10)
%
%		comprose(2,-1,1,2,0,20,'LineWidth',2,'FaceColor',.5*[1,1,1])
%
%		h = comprose(1,2.5,16,1,-10);
%		set(h(:,1),'FaceColor',[0,.5,.5],'EdgeColor',.5*[1,1,1])
%		set(h(:,2),'FaceColor',.9*[1,1,1])
%
%       See also PATCH.
%
%	Author: Francois Beauducel <beauducel@ipgp.fr>
%	Created: 2012-06-24

%	Copyright (c) 2012, François Beauducel, covered by BSD License.
%	All rights reserved.
%
%	Redistribution and use in source and binary forms, with or without
%	modification, are permitted provided that the following conditions are
%	met:
%
%	   * Redistributions of source code must retain the above copyright
%	     notice, this list of conditions and the following disclaimer.
%	   * Redistributions in binary form must reproduce the above copyright
%	     notice, this list of conditions and the following disclaimer in
%	     the documentation and/or other materials provided with the distribution
%
%	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%	AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%	IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%	ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%	LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%	CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%	SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%	INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%	CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%	ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%	POSSIBILITY OF SUCH DAMAGE.

if nargin < 5
	error('Not enough input arguments.')
end

if nargin > 5 & isnumeric(varargin{1})
	fs = varargin{1};
	varargin = varargin(2:end);
else
	fs = 0;
end

if ~isnumeric(x) | ~isnumeric(y) | ~isnumeric(n) | ~isnumeric(w) | ~isnumeric(az)
	error('X,Y,N,W, and AZ must be numeric.')
end

if all(n ~= [1,4,8,16])
	error('N must be equal to 1, 4, 8 or 16.')
end

hh = [];

for k = n:-1:1
	h = branch(x,y,k,w,az,fs,varargin{:});
	hh = [hh;h];
end

if nargout > 0
	ho = hh;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = branch(x,y,k,w,az,fs,varargin)

cpt = {'N','E','S','W','NE','SE','SW','NW', ...
	'NNE','ENE','ESE','SSE','SSW','WSW','WNW','NNW'};
cph = {'center','left','center','right','left','left','right','right', ...
	'left','left','left','left','right','right','right','right'};
cpv = {'bottom','middle','top','middle','bottom','top','top','bottom', ...
	'bottom','bottom','top','top','top','top','bottom','bottom'};

if k <= 4
	a = (k-1)*90 + az;
else
	if k <= 8
		a = (k-8)*90 - 45 + az;
		w = w*.8;
		fs = fs*.8;
	else
		a = (k-16)*45 - 22.5 + az;
		w = w*.6;
		fs = fs*.6;
	end
end

r = .15;
xx = w*[0 0 r];
yy = w*[0 1 r];
x1 =  xx*cosd(a) + yy*sind(a);
y1 = -xx*sind(a) + yy*cosd(a);
x2 = -xx*cosd(a) + yy*sind(a);
y2 =  xx*sind(a) + yy*cosd(a);
h(1) = patch(x1 + x,y1 + y,'k',varargin{:});
h(2) = patch(x2 + x,y2 + y,'k',varargin{:},'FaceColor','w');

% Cardinal points
xp = 1.03*x1(2) + x;
yp = 1.03*y1(2) + y;
h(3) = text(xp,yp,cpt{k}, ...
	'HorizontalAlignment',cph{k},'VerticalAlignment',cpv{k},'Rotation',-az);
if fs > 0
	set(h(3),'FontSize',fs','FontWeight','bold')
else
	set(h(3),'Visible','off')
end


text:
function out=cosh(x)

out=(exp(x)+exp(-x))/2;

author: opencobra

repo: cobratoolbox


text:
%%% ADMM
% A = L + S + Z
% L: low-rank
% S: sparse
% Z: stochastic or deterministic pertubation
function results = ADMM(M)
N = 3;
[mm,nn] = size(M);
g2_max = norm(M(:),inf); % Inf-norm
g3_max = norm(M); % L2-norm
g2 = 0.15*g2_max;
g3 = 0.15*g3_max;
MAX_ITER = 100;
ABSTOL   = 1e-4;
RELTOL   = 1e-2;
lambda = 1;
rho = 1/lambda;
Z = zeros(mm,nn); S = zeros(mm,nn); L = zeros(mm,nn);
z = zeros(mm,N*nn); U = zeros(mm,nn);
fprintf('\n%3s\t%10s\t%10s\t%10s\t%10s\t%10s\n', 'iter', ...
    'r norm', 'eps pri', 's norm', 'eps dual', 'objective');
for k = 1:MAX_ITER
    B = avg(Z, S, L) - M./N + U;

    % x-update
    Z = (1/(1+lambda))*(Z - B);
    S = prox_l1(S - B, lambda*g2);
    L = prox_matrix(L - B, lambda*g3, @prox_l1);

    % (for termination checks only)
    x = [Z S L];
    zold = z;
    z = x + repmat(-avg(Z, S, L) + M./N, 1, N);

    % u-update
    U = B;

    % diagnostics, reporting, termination checks
    results.objval(k)   = objective(Z, g2, S, g3, L);
    results.r_norm(k)   = norm(x - z,'fro');
    results.s_norm(k)   = norm(-rho*(z - zold),'fro');
    results.eps_pri(k)  = sqrt(mm*nn*N)*ABSTOL + RELTOL*max(norm(x,'fro'), norm(-z,'fro'));
    results.eps_dual(k) = sqrt(mm*nn*N)*ABSTOL + RELTOL*sqrt(N)*norm(rho*U,'fro');

    if k == 1 || mod(k,1) == 0
      fprintf('%4d\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.2f\n', k, ...
        results.r_norm(k), results.eps_pri(k), results.s_norm(k), results.eps_dual(k), results.objval(k));
    end
    if results.r_norm(k) < results.eps_pri(k) && results.s_norm(k) < results.eps_dual(k)
       break;
    end
end
results.iter = k;
results.Z = Z;
results.S = S;
results.L = L;
M_hat = L + S + Z;

error = norm(M_hat(:)-M(:))/norm(M(:));
disp(['Error: ' num2str(error)]);
end


text:
function tests = test_spm_ncFpdf
% Unit Tests for spm_ncFpdf
%__________________________________________________________________________
% Copyright (C) 2018 Wellcome Trust Centre for Neuroimaging

% $Id: test_spm_ncFpdf.m 7263 2018-02-21 13:30:16Z guillaume $

tests = functiontests(localfunctions);

function test_spm_ncFpdf_1(testCase)
exp = 0.230361989229139; % ncfpdf(1,1,10,0) = fpdf(1,1,10)
act = spm_ncFpdf(1,[1,10],0);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

exp = 0.339916677089114; % ncfpdf(1,2,12,0)
act = spm_ncFpdf(1,[2,12],0);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

exp = 0.187052298433871; % ncfpdf(2,2,10,2)
act = spm_ncFpdf(2,[2,10],2);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

function test_spm_ncFpdf_2(testCase)
exp = [0 % ncfpdf(0:0.3:5,3,24,pi)
   0.237214867508572
   0.305378558226511
   0.326510412867316
   0.320250126742450
   0.298516612369365
   0.269121488384202
   0.237061968288903
   0.205374511459040
   0.175766182458195
   0.149074704100068
   0.125592498668479
   0.105287389552812
   0.087948186675270
   0.073277611204176
   0.060949393292047
   0.050641638818210]';
act = spm_ncFpdf(0:0.3:5,[3,24],pi);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

author: spm


text:
function timelock = fttimelock(this, chanind, timeind, trialind, freqind)
% Method for converting meeg object to Fieldtrip timelock/freq struct
% FORMAT timelock = fttimelock(this, chanind, timeind, trialind, freqind)
%
% The method support both time and TF data and outputs different variants
% of timelock or freq FT struct depending on the dataset type and requested
% data dimensions.
% _______________________________________________________________________
% Copyright (C) 2008-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: fttimelock.m 5063 2012-11-16 12:10:41Z vladimir $

if nargin < 2 || isempty(chanind)
    chanind = 1:nchannels(this);
end

if nargin < 3 || isempty(timeind)
    timeind = 1:nsamples(this);
end

if nargin < 4 || isempty(trialind)
    trialind = 1:ntrials(this);
end

if strncmpi(transformtype(this),'TF',2) && ...
        (nargin < 5 || isempty(freqind))
     freqind = 1:nfrequencies(this);
end

timelock             = [];
timelock.label       = chanlabels(this, chanind)';

if isequal(transformtype(this), 'time')
   if isequal(type(this), 'continuous')
            error('For continuous data use ftraw method');
   end

   if isequal(type(this), 'single') || length(trialind)>1
       timelock.dimord  = 'rpt_chan_time';
       timelock.trial   =  permute(this.data.y(chanind, timeind, trialind), [3 1 2]);
   else
       timelock.dimord  = 'chan_time';
       timelock.avg     =  spm_squeeze(this.data.y(chanind, timeind, trialind), 3);
   end

   timelock.time       = time(this, timeind);

elseif strncmpi(transformtype(this),'TF',2)
    if length(timeind)>1
        if isequal(type(this), 'single') || length(trialind)>1
            timelock.dimord    = 'rpt_chan_freq_time';
            timelock.powspctrm = permute(this.data.y(chanind, freqind, timeind, trialind), [4 1 2 3]);
        else
            timelock.dimord     = 'chan_freq_time';
            timelock.powspctrm  =  spm_squeeze(this.data.y(chanind, freqind, timeind, trialind), 3);
        end

        timelock.time       = time(this, timeind);
    else
        if isequal(type(this), 'single') || length(trialind)>1
            timelock.dimord    = 'rpt_chan_freq';
            timelock.powspctrm = spm_squeeze(permute(this.data.y(chanind, freqind, timeind, trialind), [4 1 2 3]), 4);
        else
            timelock.dimord     = 'chan_freq';
            timelock.powspctrm  =  spm_squeeze(this.data.y(chanind, freqind, timeind, trialind), [3 4]);
        end
    end

    timelock.freq      = frequencies(this, freqind);
else
    error('Unknown transform type.');
end

if length(trialind)>1

    clist      =  condlist(this);
    condlabels = conditions(this, trialind);
    timelock.trialinfo = 0*trialind;

    for k = 1:numel(clist)
        fprintf('mapping condition label "%s" to condition code %d\n', clist{k}, k);
        sel = strcmp(clist{k}, condlabels);
        timelock.trialinfo(sel) = k;
    end

end

if ~isempty(sensors(this, 'MEG'))
    timelock.grad = sensors(this, 'MEG');
end

if ~isempty(sensors(this, 'EEG'))
    timelock.elec = sensors(this, 'EEG');
end


text:
function [R] = ipsrf(varargin)
%IPSRF Interval Potential Scale Reduction Factor
%
%   [R] = IPSRF(X) or
%   [R] = IPSRF(x1,x2,...,xs)
%   returns "Potential Scale Reduction Factor" (PSRF) for
%   collection of MCMC-simulations. X is a NxDxM matrix
%   which contains M MCMC simulations of length N, each with
%   dimension D. MCMC-simulations can be given as separate
%   arguments x1,x2,... which should have the same length.
%
%   Returns
%     R     PSRF in a row vector of length D
%
%   The idea of the PSRF is that if R is not near 1 (below 1.1 for
%   example) one may conclude that the tested samples were not from
%   the same distribution (chain might not have been converged
%   yet). Instead of normality assumption, 80% empirical intervals
%   are used to compute R.
%
%   If only one simulation is given, the factor is calculated
%   between first and last third of the chain. Note that use of
%   only one chain will produce over-optimistic result.
%
%   Method is from:
%      Brooks, S.P. and Gelman, A. (1998) General methods for
%      monitoring convergence of iterative simulations. Journal of
%      Computational and Graphical Statistics. 7, 434-455.
%
%   See also
%     CIPSRF, PSRF

% Copyright (C) 1999 Simo S�rkk�
% Copyright (C) 2004 Aki Vehtari
%
% This software is distributed under the GNU General Public
% Licence (version 3 or later); please refer to the file
% Licence.txt, included with the software, for details.

% In case of one argument split to two halves (first and last thirds)
onechain=0;
if nargin==1
  X = varargin{1};
  if size(X,3)==1
    n = floor(size(X,1)/3);
    x = zeros([n size(X,2) 2]);
    x(:,:,1) = X(1:n,:);
    x(:,:,2) = X((end-n+1):end,:);
    X = x;
    onechain=1;
  end
elseif nargin==0
  error('Cannot calculate PSRF of scalar');
else
  X = zeros([size(varargin{1}) nargin]);
  for i=1:nargin
    X(:,:,i) = varargin{i};
  end
end

[N,D,M]=size(X);

if N<1
  error('Too few samples');
end

W = zeros(1,D);
V = zeros(1,D);
for d=1:D
  x=X(:,d,:);
  W(1,d)=mean(diff(prctile(x,[10 90])));
  V(1,d)=diff(prctile(x(:),[10 90]));
end
R = V./W;


text:
function res = getset(this, parent, fieldname, ind, values)
% Generic method for getting and setting multiple fields of meeg struct
% FORMAT res = getset(this, parent, fieldname, ind, values)
% _______________________________________________________________________
% Copyright (C) 2008-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: getset.m 6535 2015-08-25 11:45:26Z vladimir $

this = struct(this);

if nargin == 3 || ~isnumeric(ind)
    try
        ind = 1:numel(getfield(this, parent));
    catch
        res = [];
        return;
    end
end

% Get
if nargin <= 4
    res = cell(1, length(ind));
    for i = 1:length(ind)
        res{i} = getfield(this, parent, {ind(i)}, fieldname);
    end

    if isempty(res) || (all(cellfun('isclass', res, 'double') & (cellfun(@numel, res) == 1)))
        res = [res{:}];
    end

    if iscell(res) && (numel(res) == 1) && (numel(getfield(this, parent)) == 1) &&...
            strcmp(parent, 'trials') && strcmp(this.type, 'continuous')
        res = res{1};
    end

    return
end

% Set
if nargin == 5
    % This might fail in some pathological cases, but not in what it's
    % supposed to be used for.
    if (isnumeric(values) || islogical(values)) && (length(values) == length(ind))
        values = num2cell(values);
    end

    if iscell(values) && ~(numel(values)==1 || numel(values) == length(ind))
        error('Illegal assignment: cannot match values and indices.');
    end

    for i = 1:length(ind)
        if iscell(values)
            this = setfield(this, parent, {ind(i)}, fieldname, values{i});
        else
            this = setfield(this, parent, {ind(i)}, fieldname, values);
        end
    end
    % getset is sometimes used on subfields of meeg then checkmeeg should
    % not be used
    if all(isfield(this, {'type', 'Nsamples', 'Fsample', 'timeOnset'}))
        res = meeg(this);
    else
        res = this;
    end
end


text:
classdef Stiffness3DTensor < SymmetricFourthOrder3DTensor ...
                           & StiffnessDescriptor

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%COLLISIONS Conduct collision checking for SerialLink objects
%
% Uses the point data stored in the points property of SerialLink
% objects to conduct point-primitive collision checking with objects of
% the CollisionModel class. The function does not currently check the
% base of the SerialLink object!
%
% Copyright (C) Bryan Moutrie, 2013-2014
% Licensed under the GNU Lesser General Public License
% see full file for full statement
%
% This file requires file(s) from The Robotics Toolbox for MATLAB (RTB)
% by Peter Corke (www.petercorke.com), see file for statement
%
% Syntax:
%  (1) c = robot.collisions(q, cmdl)
%  (2) c = robot.collisions(q, cmdl, dyn_cmdl, dyn_T)
%  (3) c = robot.collisions(q, cmdl, dyn_cmdl)
%  (4) c = robot.collisions(q, [], ...)
%
%  (2) is as per (1) but also checks CollisionModel dyn_cmdl, whose
%       primitives move as specified by the frames in dyn_T
%  (3) is as per (2) but assumes dyn_cmdl is the robot's tool, i.e. the
%       the robot's end-effector frame (including tool transform) is
%       used in lieu of dyn_T
%  (4) is as per (2) or (3) but with no static objects
%
% Outputs:
%  c : mx1 vector, where m is the rows of q. Elements are false if no
%       collision exists or true if a collision exists
%
% Inputs:
%  q        : mxn matrix of m joint configurations of the robot, which
%              has n links
%  cmdl     : CollisionModel object, which is static
%  dyn_cmdl : CollisionModel object which move according to dyn_T
%  dyn_T    : Transforms applied to the primitives before the transform
%              of their own property being applied. The size of dyn_T
%              should be 4x4xmxp, where p is the number of primitives
%              in dyn_cmdl. The ith shape in dyn_cmdl corresponds to the
%              ith value of the 4th dimension. Use cat to create dyn_T.
%
% See also cat CollisionModel SerialLink

% LICENSE STATEMENT:
%
% This file is part of pHRIWARE.
%
% pHRIWARE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as
% published by the Free Software Foundation, either version 3 of
% the License, or (at your option) any later version.
%
% pHRIWARE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public
% License along with pHRIWARE.  If not, see <http://www.gnu.org/licenses/>.
%
% RTB LIBRARY:
%
% Copyright (C) 1993-2014, by Peter I. Corke
% http://www.petercorke.com
% Released under the GNU Lesser General Public license

function c = collisions(robot, q, cmdl, dyn_cmdl, dyn_T)

% VERSION WITH BASE CHECKING
% pts = robot.points;
pts = robot.points(end-robot.n+1:end);
for i = length(pts): -1: 1
    numPts(i) = size(pts{i},1);
    pts{i}(:,4) = 1;
    pts{i} = pts{i}';
end

if isempty(cmdl), checkfuns = []; else checkfuns = cmdl.checkFuns; end
if nargin > 3
    dyn_checkfuns = dyn_cmdl.checkFuns;
    if nargin == 4 || isempty(dyn_T)
        tool = robot.tool;
        dyn_T = [];
    end
else
    dyn_checkfuns = [];
end

% VERSION WITH BASE CHECKING
% base = length(pts) - robot.n;
% switch base
%     case 0
%         % No base
%     case 1
%         T = robot.base;
%         trPts = T * pts{1};
%         points = trPts (1:3,:)';
%         if any(checkfuns{1}(points(:,1),points(:,2),points(:,3)))
%             C(:) = 1;
%             display('Base is colliding');
%             return;
%         end
%     otherwise
%         error('robot has missing or extra points');
% end

poses = size(q,1);
c = false(poses, 1);

nan = any(isnan(q),2);
c(nan) = true;
notnan = find(~nan)';

T0 = robot.base;

for p = notnan
    T = T0;
    prevPts = 0;
    for i = 1: robot.n;
        T = T * robot.links(i).A(q(p,i));
        if numPts(i) % Allows some links to be STLless
            %       VERSION WITH BASE CHECKING
            %       nextPts = prevPts+numPts(i+base);
            %       trPts(:,prevPts+1:nextPts) = T * pts{i+base};
            nextPts = prevPts+numPts(i);
            trPts(:,prevPts+1:nextPts) = T * pts{i};
            prevPts = nextPts;
        end
    end

    % Does same thing as cmdl.collision(trPoints(1:3,:)'), but
    % Does not have to access object every time - quicker
    for i = 1: length(checkfuns)
        if any(checkfuns{i}(trPts(1,:), trPts(2,:), trPts(3,:)))
            c(p) = true;
            break;
        end
    end

    % Then check the dynamic collision models, if any
    for i = 1: length(dyn_checkfuns)
        if isempty(dyn_T)
            dyn_trPts = T*tool \ trPts;
        else
            dyn_trPts = dyn_T(:,:,p,i) \ trPts;
        end
        if any(dyn_checkfuns{i}(dyn_trPts(1,:), dyn_trPts(2,:), ...
                dyn_trPts(3,:)))
            c(p) = true;
            break;
        end
    end

end
end

author: Allopart

repo: rbpf-gmapping


text:
function lo=v_lpcrf2lo(rf)
%V_LPCRF2LO Convert reflection coefficients to log area ratios LO=(RF)
%the output values are limited to about +-14.5

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: v_lpcrf2lo.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

r=max(min(rf,1-1E-6),1E-6-1);
lo=log((1-r)./(1+r));

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [Estimation_X] = EKF_propagate(Estimation_X, OdometryFromThis2Next, odom_sigma )

v = OdometryFromThis2Next(1:3);
w = OdometryFromThis2Next(4:6);

% update position and orientation
Estimation_X.position = Estimation_X.position+Estimation_X.orientation*v;
orientation=Estimation_X.orientation;
Estimation_X.orientation = Estimation_X.orientation*Exp(w);

NumberOfLandmarks = size(Estimation_X.landmarks, 2);
%Jrw = J_r(-w);
%ExpMinusM = Exp(-w);

%G1 = [-orientation*Jrw zeros(3,3);zeros(3,3)  -orientation];
%G1 = [-orientation zeros(3,3);zeros(3,3)  -orientation];
G1=[orientation zeros(3,3); -skew(orientation*v)   orientation];

G= [G1; zeros(3* NumberOfLandmarks ,6)];
odoCov=diag([w.^2;v.^2])*odom_sigma^2;
W = G*odoCov*G';

% compute matrix A_{n}
temp = repmat({ eye(3) }, NumberOfLandmarks+2,1 );
A = blkdiag(temp{:});
A(1:3,1:3) = eye(3);% ExpMinusM;
A(4:6,1:3) = -skew(orientation*v);

% final update the covariance
Estimation_X.cov = A*Estimation_X.cov*A'+W;

end

author: RomaTeng

repo: EKF-SLAM-on-Manifold