text:
function render = tbx_cfg_render
% Configuration file for toolbox 'Rendering'
%__________________________________________________________________________
% Copyright (C) 2008-2016 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: tbx_cfg_render.m 6960 2016-12-05 17:05:09Z guillaume $

if ~isdeployed, addpath(fullfile(spm('dir'),'toolbox','SRender')); end

%--------------------------------------------------------------------------
% images Input Images
%--------------------------------------------------------------------------
images         = cfg_files;
images.tag     = 'images';
images.name    = 'Input Images';
images.help    = {
    'These are the images that are used by the calculator.'
    'They are referred to as i1, i2, i3, etc in the order that they are specified.'
    }';
images.filter  = 'image';
images.ufilter = '.*';
images.num     = [1 Inf];

%--------------------------------------------------------------------------
% expression Expression
%--------------------------------------------------------------------------
expression         = cfg_entry;
expression.tag     = 'expression';
expression.name    = 'Expression';
expression.help    = {
    'Example expressions (f):'
    '    * Mean of six images (select six images)'
    '       f = ''(i1+i2+i3+i4+i5+i6)/6'''
    '    * Make a binary mask image at threshold of 100'
    '       f = ''i1>100'''
    '    * Make a mask from one image and apply to another'
    '       f = ''i2.*(i1>100)'''
    '             - here the first image is used to make the mask, which is applied to the second image'
    '    * Sum of n images'
    '       f = ''i1 + i2 + i3 + i4 + i5 + ...'''
    }';
expression.strtype = 's';
expression.num     = [2  Inf];
expression.val     = {'i1'};

%--------------------------------------------------------------------------
% thresh Surface isovalue(s)
%--------------------------------------------------------------------------
thresh         = cfg_entry;
thresh.tag     = 'thresh';
thresh.name    = 'Surface isovalue(s)';
thresh.help    = {'Enter the value at which isosurfaces through the resulting image is to be computed.'};
thresh.strtype = 'e';
thresh.num     = [1  1];
thresh.val     = {0.5};

%--------------------------------------------------------------------------
% surface Surface
%--------------------------------------------------------------------------
surface      = cfg_branch;
surface.tag  = 'surface';
surface.name = 'Surface';
surface.val  = {expression thresh };
surface.help = {'An expression and threshold for each of the surfaces to be generated.'};

%--------------------------------------------------------------------------
% Surfaces Surfaces
%--------------------------------------------------------------------------
Surfaces        = cfg_repeat;
Surfaces.tag    = 'Surfaces';
Surfaces.name   = 'Surfaces';
Surfaces.help   = {'Multiple surfaces can be created from the same image data.'};
Surfaces.values = {surface };
Surfaces.num    = [0 Inf];

%--------------------------------------------------------------------------
% SExtract Surface Extraction
%--------------------------------------------------------------------------
SExtract      = cfg_exbranch;
SExtract.tag  = 'SExtract';
SExtract.name = 'Surface Extraction';
SExtract.val  = {images Surfaces };
SExtract.help = {'User-specified algebraic manipulations are performed on a set of images, with the result being used to generate a surface file. The user is prompted to supply images to work on and a number of expressions to evaluate, along with some thresholds. The expression should be a standard matlab expression, within which the images should be referred to as i1, i2, i3,... etc. An isosurface file is created from the results at the user-specified threshold.'};
SExtract.prog = @spm_sextract;
SExtract.vout = @vout_sextract;

%--------------------------------------------------------------------------
% SurfaceFile Surface File
%--------------------------------------------------------------------------
SurfaceFile         = cfg_files;
SurfaceFile.tag     = 'SurfaceFile';
SurfaceFile.name    = 'Surface File';
SurfaceFile.help    = {
    'Filename of the surf_*.gii file containing the rendering information.'
    'This can be generated via the surface extraction routine in SPM.'
    'Normally, a surface is extracted from grey and white matter tissue class images, but it is also possible to threshold e.g. an spmT image so that activations can be displayed.'
    };
SurfaceFile.filter  = 'mesh';
SurfaceFile.ufilter = '.*';
SurfaceFile.num     = [1 1];

%--------------------------------------------------------------------------
% Red Red
%--------------------------------------------------------------------------
Red        = cfg_menu;
Red.tag    = 'Red';
Red.name   = 'Red';
Red.help   = {'The intensity of the red colouring (0 to 1).'};
Red.labels = {'0.0' '0.2' '0.4' '0.6' '0.8' '1.0'};
Red.values = {0 0.2 0.4 0.6 0.8 1};
Red.val    = {1};

%--------------------------------------------------------------------------
% Green Green
%--------------------------------------------------------------------------
Green        = cfg_menu;
Green.tag    = 'Green';
Green.name   = 'Green';
Green.help   = {'The intensity of the green colouring (0 to 1).'};
Green.labels = {'0.0' '0.2' '0.4' '0.6' '0.8' '1.0'}';
Green.values = {0 0.2 0.4 0.6 0.8 1};
Green.val    = {1};

%--------------------------------------------------------------------------
% Blue Blue
%--------------------------------------------------------------------------
Blue        = cfg_menu;
Blue.tag    = 'Blue';
Blue.name   = 'Blue';
Blue.help   = {'The intensity of the blue colouring (0 to 1).'};
Blue.labels = {'0.0' '0.2' '0.4' '0.6' '0.8' '1.0'}';
Blue.values = {0 0.2 0.4 0.6 0.8 1};
Blue.val    = {1};

%--------------------------------------------------------------------------
% Color Color
%--------------------------------------------------------------------------
Color      = cfg_branch;
Color.tag  = 'Color';
Color.name = 'Color';
Color.val  = {Red Green Blue};
Color.help = {
    'Specify the colour using a mixture of red, green and blue.'
    'For example, white is specified by 1,1,1, black is by 0,0,0 and purple by 1,0,1.'
    }';

%--------------------------------------------------------------------------
% DiffuseStrength Diffuse Strength
%--------------------------------------------------------------------------
DiffuseStrength        = cfg_menu;
DiffuseStrength.tag    = 'DiffuseStrength';
DiffuseStrength.name   = 'Diffuse Strength';
DiffuseStrength.help   = {'The strength with which the object diffusely reflects light. Mat surfaces reflect light diffusely, whereas shiny surfaces reflect speculatively.'};
DiffuseStrength.labels = {'0.0' '0.2' '0.4' '0.6' '0.8' '1.0'}';
DiffuseStrength.values = {0 0.2 0.4 0.6 0.8 1};
DiffuseStrength.val    = {0.8};

%--------------------------------------------------------------------------
% AmbientStrength Ambient Strength
%--------------------------------------------------------------------------
AmbientStrength        = cfg_menu;
AmbientStrength.tag    = 'AmbientStrength';
AmbientStrength.name   = 'Ambient Strength';
AmbientStrength.help   = {'The strength with which the object reflects ambient (non-directional) lighting.'};
AmbientStrength.labels = {'0.0' '0.2' '0.4' '0.6' '0.8' '1.0'}';
AmbientStrength.values = {0 0.2 0.4 0.6 0.8 1};
AmbientStrength.val    = {0.2};

%--------------------------------------------------------------------------
% SpecularStrength Specular Strength
%--------------------------------------------------------------------------
SpecularStrength        = cfg_menu;
SpecularStrength.tag    = 'SpecularStrength';
SpecularStrength.name   = 'Specular Strength';
SpecularStrength.help   = {'The strength with which the object specularly reflects light (i.e. how shiny it is). Mat surfaces reflect light diffusely, whereas shiny surfaces reflect speculatively.'};
SpecularStrength.labels = {'0.0' '0.2' '0.4' '0.6' '0.8' '1.0'}';
SpecularStrength.values = {0 0.2 0.4 0.6 0.8 1};
SpecularStrength.val    = {0.2};

%--------------------------------------------------------------------------
% SpecularExponent Specular Exponent
%--------------------------------------------------------------------------
SpecularExponent        = cfg_menu;
SpecularExponent.tag    = 'SpecularExponent';
SpecularExponent.name   = 'Specular Exponent';
SpecularExponent.help   = {'A parameter describing the specular reflectance behaviour. It relates to the size of the high-lights.'};
SpecularExponent.labels = {'0.01' '0.1' '10' '100'};
SpecularExponent.values = {0.01  0.1 10 100};
SpecularExponent.val    = {10};

%--------------------------------------------------------------------------
% SpecularColorReflectance Specular Color Reflectance
%--------------------------------------------------------------------------
SpecularColorReflectance        = cfg_menu;
SpecularColorReflectance.tag    = 'SpecularColorReflectance';
SpecularColorReflectance.name   = 'Specular Color Reflectance';
SpecularColorReflectance.help   = {'Another parameter describing the specular reflectance behaviour.'};
SpecularColorReflectance.labels = {'0.0' '0.2' '0.4' '0.6' '0.8'  '1.0'};
SpecularColorReflectance.values = {0 0.2 0.4 0.6  0.8 1};
SpecularColorReflectance.val    = {0.8};

%--------------------------------------------------------------------------
% FaceAlpha Face Alpha
%--------------------------------------------------------------------------
FaceAlpha        = cfg_menu;
FaceAlpha.tag    = 'FaceAlpha';
FaceAlpha.name   = 'Face Alpha';
FaceAlpha.help   = {
    'The opaqueness of the surface.'
    'A value of 1 means it is opaque, whereas a value of 0 means it is transparent.'
    }';
FaceAlpha.labels = {'0.0' '0.2' '0.4' '0.6' '0.8' '1.0'}';
FaceAlpha.values = {0 0.2 0.4 0.6 0.8 1};
FaceAlpha.val    = {1};

%--------------------------------------------------------------------------
% Object Object
%--------------------------------------------------------------------------
Object      = cfg_branch;
Object.tag  = 'Object';
Object.name = 'Object';
Object.val  = {SurfaceFile Color DiffuseStrength AmbientStrength SpecularStrength SpecularExponent SpecularColorReflectance FaceAlpha };
Object.help = {'Each object is a surface (from a surf_*.gii file), which may have a number of light-reflecting qualities, such as colour and shinyness.'};

%--------------------------------------------------------------------------
% Objects Objects
%--------------------------------------------------------------------------
Objects        = cfg_repeat;
Objects.tag    = 'Objects';
Objects.name   = 'Objects';
Objects.help   = {'Several surface objects can be displayed together in different colours and with different reflective properties.'};
Objects.values = {Object };
Objects.num    = [0 Inf];

%--------------------------------------------------------------------------
% Position Position
%--------------------------------------------------------------------------
Position         = cfg_entry;
Position.tag     = 'Position';
Position.name    = 'Position';
Position.help    = {'The position of the light in 3D.'};
Position.strtype = 'e';
Position.num     = [1  3];
Position.val     = {[100 100 100]};

%--------------------------------------------------------------------------
% Light Light
%--------------------------------------------------------------------------
Light      = cfg_branch;
Light.tag  = 'Light';
Light.name = 'Light';
Light.val  = {Position Color};
Light.help = {'Specification of a light source in terms of position and colour.'};

%--------------------------------------------------------------------------
% Lights Lights
%--------------------------------------------------------------------------
Lights        = cfg_repeat;
Lights.tag    = 'Lights';
Lights.name   = 'Lights';
Lights.help   = {'There should be at least one light specified so that the objects can be clearly seen.'};
Lights.values = {Light};
Lights.val    = {Light};
Lights.num    = [0 Inf];

%--------------------------------------------------------------------------
% SRender Surface Rendering
%--------------------------------------------------------------------------
SRender      = cfg_exbranch;
SRender.tag  = 'SRender';
SRender.name = 'Surface Rendering';
SRender.val  = {Objects Lights};
SRender.help = {'This utility is for visualising surfaces.  Surfaces first need to be extracted and saved in surf_*.gii files using the surface extraction routine.'};
SRender.prog = @spm_srender;

%--------------------------------------------------------------------------
% render Rendering
%--------------------------------------------------------------------------
render        = cfg_choice;
render.tag    = 'render';
render.name   = 'Rendering';
render.help   = {'This is a toolbox that provides a limited range of surface rendering options. The idea is to first extract surfaces from image data, which are saved in rend_*.mat files. These can then be loaded and displayed as surfaces. Note that OpenGL rendering is used, which can be problematic on some computers. The tools are limited - and they do what they do.'};
render.values = {SExtract SRender};

%==========================================================================
function dep = vout_sextract(job)
dep = cfg_dep;
for k=1:numel(job.surface),
    dep(k)            = cfg_dep;
    dep(k).sname      = ['Surface File ' num2str(k)];
    dep(k).src_output = substruct('.','SurfaceFile','()',{k});
    dep(k).tgt_spec   = cfg_findspec({{'filter','mesh'}});
end

author: spm


text:
function [ net,res,opts ] = lstm_bp( net,res,opts )
%LSTM_BP Summary of this function goes here
%   Detailed explanation goes here

    n_frames=opts.parameters.n_frames;
    n_cell_nodes=opts.parameters.n_hidden_nodes;

    %1: calculate the gradients of the data fitting transform

    for f=1:n_frames
        opts.dzdy=res.Fit{f}(numel(net{end}.layers)+1).dzdx;
        [net{4},res.Fit{f},opts] = net_bp(net{4},res.Fit{f},opts);
    end

    %2: BPTT: calculate the gradient wrt memory cell

    dzdc=0;
    for f=n_frames:-1:1
        %gradient from the output gate
        opts.dzdy=res.Gates{f}(end).x(2*n_cell_nodes+1:3*n_cell_nodes,:).*res.Fit{f}(1).dzdx;
        [net{3},res.Cell{f+1},opts] = net_bp(net{3},res.Cell{f+1},opts);
        %
        res.Cell{f+1}(1).dzdx=dzdc+res.Cell{f+1}(1).dzdx;
        dzdc=res.Cell{f+1}(1).dzdx.*res.Gates{f}(end).x(n_cell_nodes+1:2*n_cell_nodes,:);
    end

    res.Cell{1}(end+1).x=0;%just some padding

    %3: calculate the gradients of the input transform
    for f=1:n_frames
        opts.dzdy= res.Gates{f}(end).x(1:n_cell_nodes,:).*res.Cell{f+1}(1).dzdx;
        [net{2},res.Input{f},opts] = net_bp(net{2},res.Input{f},opts);
    end

    %4: calculate the gradients of the input,forget,output gates:
    for f=1:n_frames
        %input gate gradient;forget gate gradient; output gate gradient;
        opts.dzdy=[ res.Input{f}(end).x .*res.Cell{f+1}(1).dzdx;...
            res.Cell{f}(1).x .*res.Cell{f+1}(1).dzdx;...
            res.Fit{f}(1).dzdx.*res.Cell{f+1}(end).x];
        [net{1},res.Gates{f},opts] = net_bp(net{1},res.Gates{f},opts);
    end

    %%%accumulate gradients in all time frames

    res.Fit=average_gradients_in_frames(res.Fit);
    res.Input=average_gradients_in_frames(res.Input);
    res.Gates=average_gradients_in_frames(res.Gates);
    res.Cell=average_gradients_in_frames(res.Cell);

end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function pass = test_plotcoeffs(pref)
% This test just ensures that CHEBMATRIX/PLOTCOEFFS() does not crash.

% Create an invisible figure to plot on:
hfig = figure('Visible', 'off');

%% Setup

d = [-2 2];                   % function domain
x = chebfun(@(x) x, d);

%% Create a couple of CHEBMATRICES and try to call PLOTCOEFFS() on them

% Only CHEBFUNS
cm1 = [sin(x-.32); cos(100*(x+.3)); 3*tanh(50*sin(x-.2))];
cm2 = [sin(x-.32), cos(100*(x+.3)); ...
    3*tanh(50*sin(x-.2)) 4*sin(40*(x+.25))];

% CHEBFUNS and doubles
cm3 = [sin(x-.32); cos(100*(x+.3)); 3];
cm4 = [sin(x-.32), 100; ...
    3 4*sin((x+.25))];

pass(1) = doesNotCrash(@() plotcoeffs(cm1));
pass(2) = doesNotCrash(@() plotcoeffs(cm2));
pass(3) = doesNotCrash(@() plotcoeffs(cm3));
pass(4) = doesNotCrash(@() plotcoeffs(cm4));

%%

close(hfig);

end

function pass = doesNotCrash(fn)

try
    fn();
    pass = true;
catch ME
    rethrow(ME)
    pass = false;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function H = hypot(f, g, pref)
%HYPOT   Robust computation of the square root of the sum of squares.
%   H = HYPOT(F, G) returns SQRT(ABS(F).^2 + ABS(G).^2) for two CHEBFUN objects
%   F and G (or a CHEBFUN and a double) carefully computed to avoid underflow
%   and overflow.
%
% Example:
%       f = chebfun(@(x) 3*[1e300*x 1e-300*x]);
%       g = chebfun(@(x) 4*[1e300*x 1e-300*x]);
%       % h1 = sqrt(f.^2 + g.^2) % This will fail because of overflow
%       h2 = hypot(f, g)
%
% See also ABS, NORM, SQRT.

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

% Grab some preferences:
if ( nargin < 3 )
    pref = chebfunpref();
end

% Insert breaks at the roots:
if ( isa(f, 'chebfun') )
    f = addBreaksAtRoots(f, pref);
end
if ( isa(g, 'chebfun') )
    g = addBreaksAtRoots(g, pref);
end

% Call compose:
H = compose(f, @hypot, g, pref);

end


text:
function [e,cnt] = my_normest(S,St,n,tol, maxiter)
% [norm,iter_count] = my_normest( A, At, n, [tol], [maxiter] )
%   estimates the spectral norm of A using the power method,
%   where A is a function handle to compute A(x) = A*x
%   and At is a function handle to compute At(x) = A'*x
%
% Copied from MATLAB's "normest" function, but allows function handles, not just sparse matrices
    if nargin < 4, tol = 1.e-6; end
    if nargin < 5, maxiter = 20; end
    if isempty(St)
        St = S;  % we assume the matrix is symmetric;
    end
    x = ones(n,1);
    cnt = 0;
    e = norm(x);
    if e == 0, return, end
    x = x/e;
    e0 = 0;
    while abs(e-e0) > tol*e && cnt < maxiter
       e0 = e;
       Sx = S(x);
       if nnz(Sx) == 0
          Sx = rand(size(Sx));
       end
       e = norm(Sx);
       x = St(Sx);
       x = x/norm(x);
       cnt = cnt+1;
    end
end

repo: REST


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         Demonstration of Basic Utilities of JSONlab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rngstate = rand ('state');
randseed=hex2dec('623F9A9E');
clear data2json json2data

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a simple scalar value \n')
fprintf(1,'%%=================================================\n\n')

data2json=pi
saveubjson('',data2json)
json2data=loadubjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a complex number\n')
fprintf(1,'%%=================================================\n\n')

clear i;
data2json=1+2*i
saveubjson('',data2json)
json2data=loadubjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a complex matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=magic(6);
data2json=data2json(:,1:3)+data2json(:,4:6)*i
saveubjson('',data2json)
json2data=loadubjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  MATLAB special constants\n')
fprintf(1,'%%=================================================\n\n')

data2json=[NaN Inf -Inf]
saveubjson('specials',data2json)
json2data=loadubjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a real sparse matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=sprand(10,10,0.1)
saveubjson('sparse',data2json)
json2data=loadubjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a complex sparse matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=data2json-data2json*i
saveubjson('complex_sparse',data2json)
json2data=loadubjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  an all-zero sparse matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=sparse(2,3);
saveubjson('all_zero_sparse',data2json)
json2data=loadubjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  an empty sparse matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=sparse([]);
saveubjson('empty_sparse',data2json)
json2data=loadubjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  an empty 0-by-0 real matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=[];
saveubjson('empty_0by0_real',data2json)
json2data=loadubjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  an empty 0-by-3 real matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=zeros(0,3);
saveubjson('empty_0by3_real',data2json)
json2data=loadubjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a sparse real column vector\n')
fprintf(1,'%%=================================================\n\n')

data2json=sparse([0,3,0,1,4]');
saveubjson('sparse_column_vector',data2json)
json2data=loadubjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a sparse complex column vector\n')
fprintf(1,'%%=================================================\n\n')

data2json=data2json-1i*data2json;
saveubjson('complex_sparse_column_vector',data2json)
json2data=loadubjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a sparse real row vector\n')
fprintf(1,'%%=================================================\n\n')

data2json=sparse([0,3,0,1,4]);
saveubjson('sparse_row_vector',data2json)
json2data=loadubjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a sparse complex row vector\n')
fprintf(1,'%%=================================================\n\n')

data2json=data2json-1i*data2json;
saveubjson('complex_sparse_row_vector',data2json)
json2data=loadubjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a structure\n')
fprintf(1,'%%=================================================\n\n')

data2json=struct('name','Think Different','year',1997,'magic',magic(3),...
                 'misfits',[Inf,NaN],'embedded',struct('left',true,'right',false))
saveubjson('astruct',data2json,struct('ParseLogical',1))
json2data=loadubjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a structure array\n')
fprintf(1,'%%=================================================\n\n')

data2json=struct('name','Nexus Prime','rank',9);
data2json(2)=struct('name','Sentinel Prime','rank',9);
data2json(3)=struct('name','Optimus Prime','rank',9);
saveubjson('Supreme Commander',data2json)
json2data=loadubjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a cell array\n')
fprintf(1,'%%=================================================\n\n')

data2json=cell(3,1);
data2json{1}=struct('buzz',1.1,'rex',1.2,'bo',1.3,'hamm',2.0,'slink',2.1,'potato',2.2,...
              'woody',3.0,'sarge',3.1,'etch',4.0,'lenny',5.0,'squeeze',6.0,'wheezy',7.0);
data2json{2}=struct('Ubuntu',['Kubuntu';'Xubuntu';'Lubuntu']);
data2json{3}=[10.04,10.10,11.04,11.10]
saveubjson('debian',data2json,struct('FloatFormat','%.2f'))
json2data=loadubjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  invalid field-name handling\n')
fprintf(1,'%%=================================================\n\n')

json2data=loadubjson(saveubjson('',loadjson('{"ValidName":1, "_InvalidName":2, ":Field:":3, "项目":"绝密"}')))

rand ('state',rngstate);

author: plotly


text:
%load('BoxTorusMesh.mat')
%meshFit = mesh;
uh = UHfit(:,22);
meshFit = enrichMesh3D(meshFit,0);
bc = [1,1,1,1,1,1];
femFit = genNedFEM3D(meshFit,bc);
ng = 1;
X1 = meshFit.p(meshFit.t(:,1),:); X2 = meshFit.p(meshFit.t(:,2),:);
X3 = meshFit.p(meshFit.t(:,3),:); X4 = meshFit.p(meshFit.t(:,4),:);
gw = gaussWtetra(ng);
[gx,gy,gz] = gaussPtetra(X1,X2,X3,X4,ng);
feEvalBas = @EvalNed1Bas3D;

uhK = uh(femFit.g2ldof);
uK1 = 0; uK2 = 0; uK3 = 0;
for i = 1:6
    uK1 = uK1 + uhK(:,i).*feEvalBas(femFit.bas, ':', gx, gy, gz, i, 0, 1).*femFit.t_e_orit(:,i);
    uK2 = uK2 + uhK(:,i).*feEvalBas(femFit.bas, ':', gx, gy, gz, i, 0, 2).*femFit.t_e_orit(:,i);
    uK3 = uK3 + uhK(:,i).*feEvalBas(femFit.bas, ':', gx, gy, gz, i, 0, 3).*femFit.t_e_orit(:,i);
end

Vhfit = [uK1,uK2,uK3];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% separate evaluation
domain = [-1,1,-1,1,-1,1];
nx = 64;
ny = nx;
nz = nx;
mesh = genMesh3D(domain, nx, ny, nz);
mesh = enrichMesh3D(mesh,2);
tol = 10^(-12);
zCut = -0.3125;
fp1 = mesh.p(mesh.f(:,1),:);
fp2 = mesh.p(mesh.f(:,2),:);
fp3 = mesh.p(mesh.f(:,3),:);
fZid = (abs(fp1(:,3)-zCut)<tol).*(abs(fp2(:,3)-zCut)<tol).*(abs(fp3(:,3)-zCut)<tol);
fZid = find(fZid==1);
fmpt = (fp1(fZid,:)+fp2(fZid,:)+fp3(fZid,:))/3;

fMallfit = (meshFit.p(meshFit.t(:,1),:)+meshFit.p(meshFit.t(:,2),:)+...
    meshFit.p(meshFit.t(:,3),:)+meshFit.p(meshFit.t(:,4),:))/4;

tic
[tr,tj]=findtria(meshFit.p,meshFit.t,fmpt);
toc

VhCut = [uK1(tj),uK2(tj),uK3(tj)];

fmptPlan = fmpt(:,1:2);
zDT = delaunay(fmptPlan);

%%%%% trisulf plot
trisurf(zDT,fmptPlan(:,1),fmptPlan(:,2),VhCut(:,1));
shading interp

%%%%% vector field plot
quiver3(fmpt(:,1),fmpt(:,2),fmpt(:,3),VhCut(:,1),VhCut(:,2),VhCut(:,3))

%%%%%% contour plot
tricontour(fmptPlan,zDT,VhCut(:,1),10)

tid1 = find(mesh.tLoc==1);
tid1 = tid1(1:50:end);
fmpt1 = fMall(tid1,:);
[tr,tj]=findtria(meshFit.p,meshFit.t,fmpt1);
Vh1 = Vhfit(tj,:);
quiver5(fmpt1(:,1),fmpt1(:,2),fmpt1(:,3),...
    Vh1(:,1),Vh1(:,2),Vh1(:,3),5,'filled');

repo: ifem


text:
function Z = ldivide(X,Y)
%LDIVIDE Left array divide for tensor.
%
%   LDIVIDE(A,B) is called for the syntax 'A .\ B' when A or B is a tensor.
%   A and B must have the same size, unless one is a scalar.
%
%   Examples
%   X = tenrand([4 3 2],5);
%   X .\ 3
%   X .\ X
%
%   See also TENSOR, TENSOR/RDIVIDE.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

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

Z = tenfun(@ldivide,X,Y);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%*******************************************************
%
% DESCRIPTION:
% 		This script contains many useful constants for GPS
%       and related work.  It should be kept in only
%       one place so that updates are immediately available
%       to all other scripts/functions.
%
% ARGUMENTS:
% 		None, just call this script to place
% 		the constants in your workspace.
%
% OUTPUT:
% 		Variables in your current workspace.
%
% CALLED BY:
% 		Many other codes.
%
% FUNCTIONS CALLED:
% 		None.
%
% MODIFICATIONS:
%       XX-XX-02  :  Jan Weiss - Original
%       07-25-04  :  Jan Weiss - updated header.
%       10-19-04  :  Jan Weiss - Cleanup and added
%                                some conversion factors.
%                 :  See SVN log for further updates.
%
% Colorado Center for Astrodynamics Research
% Copyright 2005 University of Colorado, Boulder
%*******************************************************

% GENERAL CONSTANTS
% =========================================================================
c = 299792458;          %----> Speed of light (meters/s).
Re = 6378137 ;          %----> Earth Radius (meters)
% =========================================================================

% CONVERSION FACTORS
% =========================================================================
Hz2MHz = 1E-6;
MHz2Hz = 1E6;
s2ns = 1E9;
ns2s = 1E-9;
s2micros = 1E6;
micros2s = 1E-6;
s2ms = 1E3;
ms2s = 1E-3;
dtr = pi / 180;
rtd = 180 / pi;
m2cm = 100;
cm2m = 1 / 100;
m2mm = 1000;
mm2m = 1 / 1000;
ft2m = 0.3048;  % Source: http://www.nodc.noaa.gov/dsdt/ucg/
m2ft = 1 / 0.3048;
ns2m = c * ns2s;  % Converts time in nano-sec to distance,
                  % assuming the speed of light.
% =========================================================================

% GNSS SPECIFIC CONSTANTS
% =========================================================================
L1 = 1575.42e6;         %----> Freqs in Hz.
L2 = 1227.60e6;
L5 = 1176.45e6;

L1MHz = 1575.42;        %----> Freqs in MHz.
L2MHz = 1227.60;
L5MHz = 1176.45;

L1GHz = 1.57542;        %----> Freqs in GHz.
L2GHz = 1.22760;
L5GHz = 1.17645;

LAMBDA_L1 = c / L1;     %----> Wavelengths in meters.
LAMBDA_L2 = c / L2;
LAMBDA_L5 = c / L5;

CA_CODE_RATE = 1.023e6; %----> C/A and P code chipping rate in chips/s.
P_CODE_RATE = 10.23e6;

CA_CHIP_PERIOD = 1 / CA_CODE_RATE;   %----> C/A & P code chip periods in s.
P_CHIP_PERIOD = 1 / P_CODE_RATE;

CA_CHIP_LENGTH = c / CA_CODE_RATE;  %----> C/A & P code chip lengths in meters.
P_CHIP_LENGTH = c / P_CODE_RATE;

CA_CODE_LENGTH = 1023;  % chips
% =========================================================================

author: yandld


text:
function value = s_is_alpha ( s )

%*****************************************************************************80
%
%% S_IS_ALPHA returns TRUE if the string contains only alphabetic characters.
%
%  Discussion:
%
%    Here, alphabetic characters are 'A' through 'Z' and 'a' through 'z'.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be checked.
%
%    Output, logical VALUE, is TRUE if the string contains only
%    alphabetic characters.
%
  for i = 1 : length ( s )

    if ( ~ch_is_alpha ( s(i) ) )
      value = 0;
      return
    end

  end

  value = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [sFile, ChannelMat] = in_fopen_mff(DataFile, ImportOptions, channelsOnly)
% IN_FOPEN_MFF: Open a Philips .MFF file

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

%% ===== INSTALL/LOAD MFF LIBRARY =====
[isInstalled, errMsg] = bst_plugin('Install', 'mff', ImportOptions.DisplayMessages);
if ~isInstalled
    error(errMsg);
end

%% ===== EXTRACT MFF DIRECTORY =====
[parentFolder, file, ext] = bst_fileparts(DataFile);
if ismember(lower(ext), {'.bin', '.xml'})
    DataFile = parentFolder;
end

%% ===== READ MFF FILE WITH EEGLAB PLUGIN =====
hdr = struct();
hdr.filename = DataFile;
if channelsOnly
    hdr.EEG = LoadChanlocsOnly(DataFile);
else
    hdr.EEG = mff_import(DataFile);
end

%% ===== IMPORT FILE USING EEGLAB IMPORTER =====
% Convert electrodes positions from cm to mm, to avoid the interactive question about the spatial scaling
if isfield(hdr.EEG, 'chanlocs') && ~isempty(hdr.EEG.chanlocs) && isfield(hdr.EEG.chanlocs(1), 'X')
    for iChan = 1:length(hdr.EEG.chanlocs)
        hdr.EEG.chanlocs(iChan).X = hdr.EEG.chanlocs(iChan).X * 10;
        hdr.EEG.chanlocs(iChan).Y = hdr.EEG.chanlocs(iChan).Y * 10;
        hdr.EEG.chanlocs(iChan).Z = hdr.EEG.chanlocs(iChan).Z * 10;
    end
end
% Convert EEGLAB structure in Brainstorm structures
if channelsOnly
    sFile = [];
    if isempty(hdr.EEG)
        ChannelMat = [];
    else
        ChannelMat = in_channel_eeglab_set(hdr);
    end
else
    [sFile, ChannelMat] = in_fopen_eeglab(hdr, ImportOptions);
    sFile.format       = 'EEG-EGI-MFF';
    sFile.device       = 'MFF';
end

if ~isempty(ChannelMat)
    ChannelMat.Comment = 'MFF channels';
end

% Rectify sensors positions: mff_import saves the orientation with nose a +Y, brainstorm needs the nose as +X axis
if ~isempty(ChannelMat) && ~isempty(ChannelMat.Channel)
    % Define transformation: Rotation(-90deg/Z)  +  Translation(+40mm/Z)
    angleZ = -pi/2;
    Transf90Z = [cos(angleZ), -sin(angleZ), 0, 0; ...
                 sin(angleZ),  cos(angleZ), 0, 0; ...
                           0,            0, 1, 0.040; ...
                           0,            0, 0, 1];
    % Apply transformation
    ChannelMatFix = channel_apply_transf({ChannelMat}, Transf90Z);
    if ~isempty(ChannelMatFix)
        ChannelMat = ChannelMatFix{1};
    end
end

end

%% ===== LOAD CHANNEL LOC =====
function EEG = LoadChanlocsOnly(mffFile)
    % If we're only importing the channel file, create a EEG structure
    % ourselves, just like calling mff_import() without loading the data.
    EEG = struct();
    EEG.data = [];
    EEG.chaninfo = struct();
    EEG.chaninfo.nosedir = '+Y';
    try
        [EEG.chanlocs, EEG.ref] = mff_importcoordinates(mffFile);
    catch
        EEG = [];
        return;
    end
    try
        EEG = eeg_checkchanlocs(EEG);
    catch
    end
    EEG.nbchan = length(EEG.chanlocs);
end

author: brainstorm-tools


text:
% GP_COV_KRON - Covariance function as a product of two covariance
%                  functions.
%
% Usage:
%
%   COVFUNC = GP_COV_PRODUCT(COVFUNC1, COVFUNC2, ...)
%
% The returned covariance function is called as
%
%   K = COVFUNC(THETA)
%   [K, DK] = COVFUNC(THETA)
%
% where THETA is a vector of parameters. The parameters of COVFUNC1 and
% COVFUNC2 are concatenated into a single parameter vector.
%
% If no arguments are given, dimensionality information is returned:
%
%   [N_THETA, N1, N2] = COVFUNC()
%
% N_THETA : Number of non-fixed parameters
% N1      : Number of rows in the covariance matrix
% N2      : Number of columns in the covariance matrix
%
% See also GP_COV, GP_COV_SUM.

% Last modified 2012-03-14
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@aalto.fi)

function func = gp_cov_kron(varargin)

if nargin < 2
  error('Must give at least two covariance functions')
end

func = @get_covariance;

covfuncs = varargin ;

% Number of covariance functions
n_funcs = nargin ;

% Number of parameters for each covariance function
n_theta = zeros(n_funcs,1) ;
% Dimensionalities of each covariance function
M = zeros(n_funcs,1) ;
N = zeros(n_funcs,1) ;
% Extract the values from the covariance functions
for i = 1:n_funcs
  [n_theta(i), M(i), N(i)] = varargin{i}();
end

% Indices of the hyperparameters for each covariance function
ind_theta = cell(n_funcs,1);
for i = 1:n_funcs
  ind_theta{i} = (1+sum(n_theta(1:(i-1)))):(sum(n_theta(1:i))) ;
end

%[n_theta2, M2, N2] = covfunc2();

%if M1 ~= M2 || N1 ~= N2
%  error('Can''t multiply covariance matrices with different dimensionalities');
%end

  function varargout = get_covariance(theta)

  if nargout == 0
    nout = 1;
  else
    nout = nargout;
  end

  varargout = cell(nout,1);

  out = cell(nout,1);
  %out2 = cell(nout,1);
  varargout = cell(nout,1);

  % Return only dimension information if requested
  if nargin == 0
    if nout >= 1
      varargout{1} = sum(n_theta); % number of parameters
      if nout >= 2
        varargout{2} = prod(M) ; % dimensionalities
        if nout >= 3
          varargout{3} = prod(N); % dimensionalities
        end
      end
    end
    return
  end

  if numel(theta) ~= sum(n_theta)
    error('Wrong number of parameters (%d), should be %d', numel(theta), ...
          sum(n_theta));
  end

  % Initialize covariance matrix and gradients
  varargout{1} = 1 ;
  if nout >= 2
    varargout{2} = cell(sum(n_theta),1) ;
    for i = 1:sum(n_theta)
      varargout{2}{i} = 1 ;
    end
  end

  % Compute the covariance matrix
  for i = 1:n_funcs
    [out{:}] = covfuncs{i}(theta(ind_theta{i}));
    varargout{1} = kron(varargout{1}, out{1});

    % Compute the derivative
    if nout >= 2 && ~isempty(ind_theta{i})
      % No derivative for these hyperparameters
      for n = 1:(ind_theta{i}(1)-1)
        varargout{2}{n} = kron(varargout{2}{n}, out{1}) ;
      end
      % Derivative for these hyperparameters
      for n = 1:n_theta(i)
        ind = ind_theta{i}(1) + n - 1 ;
        varargout{2}{ind} = kron(varargout{2}{ind}, out{2}{n}) ;
      end
      % No derivative for these hyperparameters
      for n = (ind_theta{i}(end)+1):sum(n_theta)
        varargout{2}{n} = kron(varargout{2}{n}, out{1}) ;
      end
    end
  end

  end

end

%function X = kronecker(varargin)

%X = varargin{1} ;
%for i = 2:nargin
%  X = kron(X, varargin{i}) ;
%end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function geometric_cdf_values_test ( )

%*****************************************************************************80
%
%% GEOMETRIC_CDF_VALUES_TEST demonstrates the use of GEOMETRIC_CDF_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GEOMETRIC_CDF_VALUES_TEST:\n' );
  fprintf ( 1, '  GEOMETRIC_CDF_VALUES stores values of \n' );
  fprintf ( 1, '  the Geometric Probability Cumulative Density Function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     X      P       CDF\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, p, cdf ] = geometric_cdf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %4d  %8f  %24.16f\n', x, p, cdf );

  end

  return
end


text:
function R = tapas_physio_scaleorthmean_regressors(cardiac_sess)
% create orthogonalized cardiac regressors (Fourier expansion) to ensure
% numerical stability of SVD in spm_spm, if regressors are highly
% correlated (e.g. after triggered acquisition)
%
% USAGE:
%   R = tapas_physio_scaleorthmean_regressors(cardiac_sess)
%
% -------------------------------------------------------------------------
% INPUT:
%   cardiac_sess    - RETROICOR-generated cardiac regressors for a design
%                     matrix [Nscans X (2*expansion_order)]
%
% -------------------------------------------------------------------------
% OUTPUT:
%   cardiac_sess    - scaled to max, orthogonalized, mean-free cardiac regressors
%
% -------------------------------------------------------------------------
% Lars Kasper, August 2011
% Copyright (C) 2013, Institute for Biomedical Engineering, ETH/Uni Zurich.
%
% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

if isempty(cardiac_sess)
    R = [];
else
    R       = scale_max(spm_orth(mean_centre(cardiac_sess)));
end

function Rout = mean_centre( R )
Rout = R - repmat(mean(R), length(R), 1);

function Rout = scale_max( R )
Rout = R./repmat(max(abs(R)), length(R), 1);

author: translationalneuromodeling

repo: tapas


text:
%kkfraction 'Output = Fractional Part of Input'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kfraction.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'Input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kkfraction(i, {'i','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kfraction - Output = Fractional Part of Input
%
%  DESCRIPTION
%
% The "Fractional Part" operator returns the fractional part of
% each data point in the input data object, \fBInput".
%
% Executing "Fractional Part" runs the program \fIkarith1\fP
% with the -fract flag.
%
%  "\fBData Type"" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_1input
%
%  "\fBMap Data"" 5
% .cI $DATAMANIP/repos/shared/man/sections/map_1input
%
%  "\fBValidity Mask"" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_1input
%
%  "\fBExplicit Location and Time Data"" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%  "\fBFailure Modes""
% .cI $DATAMANIP/repos/shared/man/sections/fail_1input
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::karith1
%
%  RESTRICTIONS
% Operations on complex data are not supported at this time.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkfraction(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkfraction(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output'};
maxval={0,0};
minval={0,0};
istoggle=[0,0];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'karith1"  -fract'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% inhomo_iso.m demonstrates Edge-preserving non-linear
% inhomogeneous linear isotropic filtering for denoising
%
% Copyright (c) Ritwik Kumar, Harvard University 2010
%               www.seas.harvard.edu/~rkkumar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
%read the image
im = image_read('synimgn2');
w = double(im);

% set up finite difference parameters
alpha =0.5;
k = 1;
h = 1;

lambda = (alpha^2)*(k/(h^2));

[m n] = size(w);
% stick image into a vector
w_vec = reshape(w,n*n,1);
w_old = w_vec;
w_new = w_vec;

%smooth the image using Gaussian filtering
im_smth = filter_function(w,3);

%required of g() calculation
[dx_im_smth dy_im_smth] = gradient(im_smth);
gr_im_smth = dx_im_smth.^2 + dy_im_smth.^2;

% element by element g() calculation
[mmm nnn] = size(gr_im_smth);
for i=1:mmm
    for j=1:nnn
        g(i,j) = 1/(1+(gr_im_smth(i,j)/36));
    end
end
g = g';
jj=1;

figure;
for k=1:200   % for each iteration
    for i=1:n*n % solve using JAcobi iteration scheme
        row = ceil(i/n); %compute what row this pixel belongs to in original image
        col = i - (row-1) * n; % compute the column similarly

        %different if conditions handles pixels at different location in
        %the image as depending on their location they may or may not have
        %all their neighbor pixels which will be required for finite
        %differences

        if((col > 1) && (col < n) && (row > 1) && (row < n))
            xx = [(lambda * (-g(row,col) + g(row,col+1) + g(row+1,col) -g(row+1,col+1))) * ...
                    (- w_old(i) + w_old(i+1) - w_old(i+n) + w_old(i+n+1))];

            s = [-(lambda * g(row,col)) * (w_old(i+1) + w_old(i-1) + w_old(i+n) + w_old(i-n))] ...
                + xx;
            w_new(i) = (-s + w_old(i))/(1+4*lambda*g(row,col));
        elseif((row == 1) && (col > 1) && (col < n))
            xx = [(lambda * (-g(row,col) + g(row,col+1) + g(row+1,col) -g(row+1,col+1))) * ...
                    (- w_old(i) + w_old(i+1) - w_old(i+n) + w_old(i+n+1))];

            s = [-(lambda * g(row,col)) * (w_old(i+1) + w_old(i-1) + w_old(i+n))] ...
                + xx;

            w_new(i) = (-s + w_old(i))/(1+4*lambda);
        elseif((row == n) && (col > 1) && (col < n))
            xx = [(lambda * (-g(row,col) + g(row,col+1))) * ...
                    (- w_old(i) + w_old(i+1))];

            s = [-(lambda * g(row,col)) * (w_old(i+1) + w_old(i-1) + w_old(i-n))] ...
                + xx;

            w_new(i) = (-s + w_old(i))/(1+4*lambda);
        elseif((col == 1) && (row > 1) && (row < n))
            xx = [(lambda * (-g(row,col) + g(row,col+1) + g(row+1,col) -g(row+1,col+1))) * ...
                    (- w_old(i) + w_old(i+1) - w_old(i+n) + w_old(i+n+1))];

            s = [-(lambda * g(row,col)) * (w_old(i+1) + w_old(i+n) + w_old(i-n))] ...
                + xx;

            %s = -(lambda) * (w_old(i+1) + w_old(i+n) + w_old(i-n));
            w_new(i) = (-s + w_old(i))/(1+4*lambda);
        elseif((col == n) && (row > 1) && (row < n))
            xx = [(lambda * (-g(row,col) + g(row+1,col) )) * ...
                    (- w_old(i)  - w_old(i+n))];

            s = [-(lambda * g(row,col)) * ( w_old(i-1) + w_old(i+n) + w_old(i-n))] ...
                + xx;

            w_new(i) = (-s + w_old(i))/(1+4*lambda);
        elseif((col==1) && (row==1))
            xx = [(lambda * (-g(row,col) + g(row,col+1) + g(row+1,col) -g(row+1,col+1))) * ...
                    (- w_old(i) + w_old(i+1) - w_old(i+n) + w_old(i+n+1))];

            s = [-(lambda * g(row,col)) * (w_old(i+1) + w_old(i+n))] ...
                + xx;

            w_new(i) = (-s + w_old(i))/(1+4*lambda);
        elseif((col==n) && (row==1))
                        xx = [(lambda * (-g(row,col) + g(row+1,col))) * ...
                (- w_old(i)  - w_old(i+n) )];

            s = [-(lambda * g(row,col)) * (w_old(i-1) + w_old(i+n))] ...
                + xx;

            w_new(i) = (-s + w_old(i))/(1+4*lambda);
        elseif((col==1) && (row==n))
                       xx = [(lambda * (-g(row,col) + g(row,col+1) )) * ...
                (- w_old(i) + w_old(i+1) )];

            s = [-(lambda * g(row,col)) * (w_old(i+1) + w_old(i-n))] ...
                + xx;

            w_new(i) = (-s + w_old(i))/(1+4*lambda);
        elseif((col==n) && (row==n))
                        xx = [(lambda * (-g(row,col) )) * ...
                (- w_old(i) )];

            s = [-(lambda * g(row,col)) * ( w_old(i-1) + w_old(i-n))] ...
                + xx;

            w_new(i) = (-s + w_old(i))/(1+4*lambda);
        end
    end
    w_old = w_new;
    if((k==10)||(k==50)||(k==100)||(k==200)) % outputting the values at different iterations
        subplot(2,2,jj)
        imshow(uint8(reshape(w_new,n,n)));
        imwrite(uint8(reshape(w_new,n,n)),strcat('iso-synimgn2',int2str(jj),'.bmp'));
                title('Linear Inhomogeneous Isotropic Filtering (Implicit Euler Method)');
        xlabel(strcat('Lenna at t=',int2str(k),', Alpha=0.5'));

        jj=jj+1;
    end

end

author: Sable

repo: mcbench-benchmarks


text:
%This Matlab script can be used to reproduce Figures 4.18-19 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%Number of BSs
L = 16;

%Number of UEs per BS
K = 10;

%Define the range of BS antennas
Mrange = 10:10:100;

%Extract maximum number of BS antennas
Mmax = max(Mrange);

%Define the range of pilot reuse factors
fRange = [1 2 4];

%Select the number of setups with random UE locations
nbrOfSetups = 50;

%Select the number of channel realizations per setup
nbrOfRealizations = 500;

%% Propagation parameters

%Communication bandwidth
B = 20e6;

%Total uplink transmit power per UE (mW)
p = 100;

%Total downlink transmit power per UE (mW)
rho = 100;

%Noise figure at the BS (in dB)
noiseFigure = 7;

%Compute noise power
noiseVariancedBm = -174 + 10*log10(B) + noiseFigure;

%Select length of coherence block
tau_c = 200;

%Use the approximation of the Gaussian local scattering model
accuracy = 2;

%Angular standard deviation in the local scattering model (in degrees)
ASDdeg = 10;

%Prepare to save simulation results
sumSE_MR = zeros(length(Mrange),length(fRange),nbrOfSetups);
sumSE_ZF = zeros(length(Mrange),length(fRange),nbrOfSetups);
sumSE_SMMSE = zeros(length(Mrange),length(fRange),nbrOfSetups);
sumSE_RZF = zeros(length(Mrange),length(fRange),nbrOfSetups);
sumSE_MMMSE = zeros(length(Mrange),length(fRange),nbrOfSetups);

%% Go through all setups
for n = 1:nbrOfSetups

    %Output simulation progress
    disp([num2str(n) ' setups out of ' num2str(nbrOfSetups)]);

    %Compute channel statistics for one setup
    [R,channelGaindB] = functionExampleSetup(L,K,Mmax,accuracy,ASDdeg);

    %Compute the normalized average channel gain, where the normalization
    %is based on the noise power
    channelGainOverNoise = channelGaindB - noiseVariancedBm;

    %Go through all number of antennas
    for m = 1:length(Mrange)

        %Output simulation progress
        disp([num2str(m) ' antennas out of ' num2str(length(Mrange))]);

        %Go through all pilot reuse factors
        for s = 1:length(fRange)

            %Extract pilot reuse factor
            f = fRange(s);

            %Generate channel realizations with estimates and estimation
            %error correlation matrices
            [Hhat,C,tau_p,Rscaled,H] = functionChannelEstimates(R(1:Mrange(m),1:Mrange(m),:,:,:),channelGainOverNoise,nbrOfRealizations,Mrange(m),K,L,p,f);

            %Compute SEs with the estimation bound in Theorem 4.6 using
            %Monte Carlo simulations
            [SE_hardening_MR,SE_hardening_RZF,SE_hardening_MMMSE,SE_hardening_ZF,SE_hardening_SMMSE] = functionComputeSE_DL_hardening(H,Hhat,C,Rscaled,tau_c,tau_p,nbrOfRealizations,Mrange(m),K,L,p,rho);

            %Compute SEs with the estimation bound in Theorem 4.9 using
            %Monte Carlo simulations
            [SE_MR,SE_RZF,SE_MMMSE,SE_ZF,SE_SMMSE] = functionComputeSE_DL_estimation(H,Hhat,C,Rscaled,tau_c,tau_p,nbrOfRealizations,Mrange(m),K,L,p,rho);

            %Use the largest of the two bounds
            SE_MR(SE_hardening_MR>SE_MR) = SE_hardening_MR(SE_hardening_MR>SE_MR);
            SE_RZF(SE_hardening_RZF>SE_RZF) = SE_hardening_RZF(SE_hardening_RZF>SE_RZF);
            SE_MMMSE(SE_hardening_MMMSE>SE_MMMSE) = SE_hardening_MMMSE(SE_hardening_MMMSE>SE_MMMSE);
            SE_ZF(SE_hardening_ZF>SE_ZF) = SE_hardening_ZF(SE_hardening_ZF>SE_ZF);
            SE_SMMSE(SE_hardening_SMMSE>SE_SMMSE) = SE_hardening_SMMSE(SE_hardening_SMMSE>SE_SMMSE);

            %Save results
            sumSE_MR(m,s,n) = mean(sum(SE_MR,1));
            sumSE_ZF(m,s,n) = mean(sum(SE_ZF,1));
            sumSE_SMMSE(m,s,n) = mean(sum(SE_SMMSE,1));
            sumSE_RZF(m,s,n) = mean(sum(SE_RZF,1));
            sumSE_MMMSE(m,s,n) = mean(sum(SE_MMMSE,1));

            %Delete large matrices
            clear H Hhat C Rscaled;

        end

    end

    %Delete large matrices
    clear R;

end

%% Plot the simulation results
for s = 1:length(fRange)

    figure(s);
    hold on; box on;

    plot(Mrange,mean(sumSE_MMMSE(:,s,:),3),'rd-','LineWidth',1);
    plot(Mrange,mean(sumSE_SMMSE(:,s,:),3),'b:','LineWidth',1);
    plot(Mrange,mean(sumSE_RZF(:,s,:),3),'k-.','LineWidth',1);
    plot(Mrange,mean(sumSE_ZF(:,s,:),3),'r--','LineWidth',1);
    plot(Mrange,mean(sumSE_MR(:,s,:),3),'bs-','LineWidth',1);

    xlabel('Number of antennas (M)');
    ylabel('Average sum SE [bit/s/Hz/cell]');
    legend('M-MMSE','S-MMSE','RZF','ZF','MR','Location','NorthWest');
    ylim([0 60]);

end

author: emilbjornson

repo: massivemimobook


text:
function [xPred,SPred,exitCode]=sqrtCDEKFPred(xPrev,SPrev,a,D,AJacob,tPrev,tPred,RKOptions)
%%SQRTCDEKFPRED Predict forward a Gaussian state estimate through time when
%               the evolution of the state is described by a
%               continuous-time stochastic diferential equation using a
%               square root continuous-discrete extended Kalman filter.
%               This uses a Gaussian approximation throughout the entire
%               prediction, allowing for a solution based on deterministic
%               differential equations. The differential equations are
%               integrated forward using the RKAdaptiveOverRange function,
%               an explicit Runge-Kutta method, which can take a number of
%               options.
%
%INPUTS: xPrev The xDim X 1 state estimate at time tPrior.
%        SPrev The xDim X xDim square root state covariance matrix at time
%              tPrior.
%            a The drift function in the continuous-time stochastic
%              dynamic model. It takes the state and a time variable as its
%              arguments.
%            D The diffusion function in the continuous-time stochastic
%              dynamic model. It takes the state and a time variable as its
%              arguments.
%       AJacob The derivative wrt x of the drift function in the
%              continuous-time stochastic dynamic model. It takes the
%              state and a time variable as its arguments. If an empty
%              matrix is passed, then AJacob will be found using
%              numerical differentiation via the numDiff function with
%              default parameters.
%        tPrev The time of xPrev and SPrev.
%        tPred The time to which xPrev and SPrev should be predicted.
%    RKOptions An optional structure whose components have the same name
%              and meaning as the corresponding components in the
%              RKAdaptiveOverRange function. The possible compnents of
%              the RKOptions function (with default values if omitted)
%              are
%              RKOptions.initStepSize: (tPred-tPrev)/RKOptions.maxSteps
%              RKOptions.order:           5
%              RKOptions.solutionChoice:  0
%              RKOptions.RelTol:          1e-3.
%              RKOptions.AbsTol:          1e-6
%              RKOptions.maxSteps:        1024
%              If an empty matrix is passed in place of RKOptions, then
%              only the default values will be used.
%
%OUTPUTS: xPred The xDim X 1 predicted state estimate. If the
%               RKAdaptiveOverRange function failed, which might occur,
%               for example, if the maximum number of allowed steps is
%               too small, then an empty matrix is returned.
%         SPred The xDim X xDim predicted square root state covariance
%               matrix, or an empty matrix if the RKAdaptiveOverRange
%               failed.
%      exitCode The exit code from the RKAdaptiveOverRange function. This
%               is zero if the integration was a success. Otherwise, the
%               value identifies what the problem was. See the
%               RKAdaptiveOverRange function for more details.
%
%The algorithm is modified from the mean-covariance Runge-Kutta (MC-RK4)
%method described in [1]. Other methods in this paper are criticized in
%Section III of [2].
%
%Parts of the calculation for the derivative of the square root state
%covariance matrix are taken from Section VI of [3].
%
%REFERENCES:
%[1] P. Frogerais., J. Bellanger, and L. Senhadji. "Various ways to compute
%    the continuous-discrete extended Kalman filter," IEEE
%    Transactions on Automatic Control, vol. 57, no. 4, pp. 1000-1004,
%    2012.
%[2] G. Y. Kulikov and M. V. Kulikova. "Accurate numerical implementation
%    of the continuous-discrete extended Kalman filter," IEEE Transactions
%    on Automatic Control, vol. 59, no. 1, pp. 273-279, Jan. 2014.
%[3] D. F. Crouse, "Basic tracking using nonlinear continuous-time dynamic
%    models," IEEE Aerospace and Electronic Systems Magazine, vol. 30, no.
%    2, Part II, pp. 4-41, Feb. 2015.
%
%March 2015 David Karnick, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    xDim=length(xPrev);

    if(nargin<8)
        RKOptions=[];
    end
    if(isempty(AJacob))
        AJacob=@(x,t)numDiff(x,@(y) a(y,t),xDim);
    end

    %The default options for the RKAdaptiveOverRange function. The
    %initStepSize default option is not set until the other options are
    %read, since the default depends on the value of maxSteps.
    initStepSize=[];
    order=5;
    solutionChoice=0;
    RelTol=1e-3;
    AbsTol=1e-6;
    maxSteps=1024;
    if(~isempty(RKOptions))
        if(isfield(RKOptions,'initStepSize'))
            initStepSize=RKOptions.initStepSize;
        end
        if(isfield(RKOptions,'order'))
            order=RKOptions.order;
        end
        if(isfield(RKOptions,'solutionChoice'))
            solutionChoice=RKOptions.solutionChoice;
        end
        if(isfield(RKOptions,'RelTol'))
            RelTol=RKOptions.RelTol;
        end
        if(isfield(RKOptions,'AbsTol'))
            AbsTol=RKOptions.AbsTol;
        end
        if(isfield(RKOptions,'maxSteps'))
            maxSteps=RKOptions.maxSteps;
        end
    end

    if(isempty(initStepSize))
        initStepSize=(tPred-tPrev)/maxSteps;
    end

    %The zeros in S are not given to the differential equation solver.
    %lowerTriEls holds the indices of the lower-triangular elements.
    lowerTriEls=tril(true(xDim,xDim));
    [xVals,~,~,exitCode]=RKAdaptiveOverRange([xPrev;SPrev(lowerTriEls)],[tPrev tPred],@f,initStepSize,0,order,solutionChoice,RelTol,AbsTol,maxSteps);
    if(isempty(xVals))%If the RKAdaptiveOverRange function failed.
        xPred=[];
        SPred=[];
        return;
    end

    xSVec=xVals(:,end);
    xPred=xSVec(1:xDim,end);
    SEls=xSVec((xDim+1):end,end);
    SPred=zeros(xDim,xDim);
    SPred(lowerTriEls)=SEls;

    function dVal=f(x,t)
        %This function provides the derivatives of the state and the
        %lower-triangular square root covariance matrix at a particular
        %time.
        xCur=x(1:xDim);
        SElsCur=x((xDim+1):end);
        %Only the nonzero elements of S were passed.
        SCur=zeros(xDim,xDim);
        SCur(lowerTriEls)=SElsCur;
        PCur=SCur*SCur';

        dx=a(xCur,t);

        dPoints=D(xCur,t);
        F=AJacob(xCur,t);
        dP=PCur*F'+F*PCur'+dPoints*dPoints';
        SInv=pinv(SCur);
        dS=SCur*triLower(SInv*dP*SInv');

        dVal=[dx;dS(lowerTriEls)];
    end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function kern = linard2KernParamInit(kern)

% LINARD2KERNPARAMINIT LINARD2 kernel parameter initialisation.
% The automatic relevance determination version of the linear
% kernel (LINARD2) is the simple inner product kernel with feature
% selection applied.
%
% k(x_i, x_j) = x_i'*A* x_j
%
% where A is a diagonal matrix of values constrained to positve. These
% parameters are stored in the field 'inputScales'.
%
% SEEALSO : linKernParamInit, rbfardKernParamInit
%
% FORMAT
% DESC initialises the automatic relevance determination linear
%  kernel structure with some default parameters.
% ARG kern : the kernel structure which requires initialisation.
% RETURN kern : the kernel structure with the default parameters placed in.
%
% SEEALSO : kernCreate, kernParamInit
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006
%
% COPYRIGHT : Michalis K. Titsias, 2009

% KERN

% These parameters are restricted to positive
kern.inputScales = 0.999*ones(1, kern.inputDimension);
kern.nParams = kern.inputDimension;

kern.transforms(1).index = [1:kern.nParams];
kern.transforms(1).type = optimiDefaultConstraint('positive');

kern.isStationary = false;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef SettingsMeshUnfitted < AbstractSettings

    properties (Access = protected)
        defaultParamsName = 'paramsMeshUnfitted.json'
    end

    properties (Access = public)
        backgroundMesh
        boundaryMesh
    end

    methods (Access = public)

        function obj = SettingsMeshUnfitted(varargin)
            switch nargin
                case 1
                    obj.loadParams(varargin{1});
                case 2
                    obj.backgroundMesh = varargin{2};
                case 3
                    disp('eis');
                case 4
                    disp('eis');

            end
            obj.init();
        end

    end

    methods (Access = private)

        function init(obj)
            obj.createBackgroundMesh();
        end

        function createBackgroundMesh(obj)
            if ischar(obj.backgroundMesh)
                fileName = obj.backgroundMesh;
                femReader = FemInputReader_GiD();
                s = femReader.read(fileName);
                obj.backgroundMesh = s.mesh;
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function varargout = gradient(varargin)
%GRADIENT  Gradient of a CHEBFUN2.
%   [FX FY] = GRADIENT(F) returns the gradient of the CHEBFUN2 F,
%   where FX is the derivative of F in the x direction and FY is the derivative
%   in the y direction. Both derivatives are returned as CHEBFUN2 objects.
%
%   G = GRADIENT(F) returns a CHEBFUN2V which represents
%
%            G = ( F_x ; F_y )
%
% See also GRAD.

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

[varargout{1:nargout}] = gradient@separableApprox(varargin{:});

end


text:
function [g1, g2, g3] = sdlfmXsdlfmKernGradientBlockILJ(lfmKern1, lfmKern2, ...
    t1, t2, i, j, generalConst, generalConstGrad, covGrad, g1Mean, ...
    g2Mean, gsp1Mean, gsp2Mean, typeMeanParam, typeMeanSwitching, ...
    typeKernParam1, typeKernParam2, typeKernSwitching1, typeKernSwitching2)

% SDLFMXSDLFMKERNGRADIENTBLOCKILJ
% FORMAT
% DESC computes the gradients of the parameters for system 1 and system 2
% when i is lower than j.
% ARG lfmKern1 : structure containing parameters for the system 1
% ARG lfmKern2 : structure containing parameters for the system 2
% ARG t1 : times at which the system 1 is evaluated
% ARG t2 : times at which the system 2 is evaluated
% ARG i : interval to be evaluated for system 1
% ARG j : interval to be evaluated for system 2
% ARG generalConstant : constants evaluated with sdlfmKernComputeConstant.m
% ARG generalConstGrad : derivatives of the constants computed with
% sdlfmKernGradientConstant.m
% ARG covGrad : partial derivatives of the objective function wrt portion
% of the corresponding kernel matrix
% ARG g1Mean : gradients of the parameter of the system 1 obtained from the
% part of the function that uses the funcitons accommpanying the initial
% conditions.
% ARG g2Mean : gradients of the parameter of the system 2 obtained from the
% part of the function that uses the funcitons accommpanying the initial
% conditions.
% ARG gsp1Mean : gradient of the switching point of the system 1 obtained
% from the part of the function that uses the funcitons accommpanying the
% initial conditions.
% ARG gsp2Mean : gradient of the switching point of the system 2 obtained
% from the part of the function that uses the funcitons accommpanying the
% initial conditions.
% ARG typeMeanParam : specify the mean functions used to compute this part
% of the kernel
% ARG typeMeanSwitching : specify the functions used to compute the
% gradients of the swicthing points in part thst uses mean functions
% ARG typeKernParam1 : specify the first kernel function used to compute
% this part of the kernel
% ARG typeKernParam2 : specify the second kernel function used to compute
% this part of the kernel
% ARG typeKernSwitching1 : specify the functions used to compute the
% gradients of the swicthing points in both sides of the kernel function 1
% ARG typeKernSwitching2 : specify the functions used to compute the
% gradients of the swicthing points in both sides of the kernel function 2
% RETURN g1 : gradients of parameters for the system 1
% RETURN g2 : gradients of parameters for the system 2
% RETURN g3 : gradients of switching points
%
% COPYRIGHT : Mauricio A. Alvarez, 2010.

% KERN

if nargin < 14
    typeMeanParam = 'sdlfm';
    typeMeanSwitching = 'sdlfmv';
    typeKernParam1 = 'lfmXlfm';
    typeKernParam2 = 'lfmvXlfm';
    typeKernSwitching1= {'lfmvXlfm', 'lfmvXlfm'};
    typeKernSwitching2 ={'lfmvXlfmv', 'lfmaXlfm'};
end

g3 = [];

fhandleMeanParam = str2func([typeMeanParam 'MeanCompute']);
fhandleMeanGradParam = str2func([typeMeanParam 'MeanGradient']);
fhandleMeanGradSwitching = str2func([typeMeanSwitching 'MeanCompute']);
fhandleKernPosPos = str2func([typeKernParam1 'KernCompute']);
fhandleKernGradPosPos = str2func([typeKernParam1 'KernGradient']);
fhandleKernGradSwitchingPosPos1 = str2func([typeKernSwitching1{1} 'KernCompute']);
fhandleKernGradSwitchingPosPos2 = str2func([typeKernSwitching1{2} 'KernCompute']);
fhandleKernVelPos = str2func([typeKernParam2 'KernCompute']);
fhandleKernGradVelPos = str2func([typeKernParam2 'KernGradient']);
fhandleKernGradSwitchingVelPos1 = str2func([typeKernSwitching2{1} 'KernCompute']);
fhandleKernGradSwitchingVelPos2 = str2func([typeKernSwitching2{2} 'KernCompute']);

c2 = fhandleMeanParam(lfmKern2(1), t2, 'Pos');
e2 = fhandleMeanParam(lfmKern2(1), t2, 'Vel');

if isempty(generalConst{i,j})
    coeffPosPos = c2;
    coeffPosVel = e2;
else
    coeffPosPos = generalConst{i,j}(1,1)*c2 + generalConst{i,j}(2,1)*e2;
    coeffPosVel = generalConst{i,j}(1,2)*c2 + generalConst{i,j}(2,2)*e2;
end
g1Kern = zeros(length(lfmKern1), 5);
g2Kern = zeros(length(lfmKern1), 5);
PosPos = zeros(length(t1),1);
PosVel = zeros(length(t1),1);
for k=1:length(lfmKern1)
    PosPos = PosPos + fhandleKernPosPos(lfmKern1(k), lfmKern2(k), t1,...
        lfmKern1(k).limit);
    [g1KLocal, g2KLocal] = fhandleKernGradPosPos(lfmKern1(k), lfmKern2(k), ...
        t1, lfmKern1(k).limit, covGrad, coeffPosPos.');
    g1Kern(k,:) = g1Kern(k, :) + g1KLocal;
    g2Kern(k,:) = g2Kern(k, :) + g2KLocal;
    PosVel = PosVel + fhandleKernVelPos(lfmKern2(k), lfmKern1(k), ...
        lfmKern1(k).limit, t1).';
    [g2KLocal, g1KLocal] = fhandleKernGradVelPos(lfmKern2(k), lfmKern1(k), ...
        lfmKern1(k).limit, t1, covGrad, coeffPosVel);
    g1Kern(k,:) = g1Kern(k, :) + g1KLocal;
    g2Kern(k,:) = g2Kern(k, :) + g2KLocal;
end
[gcAlpha, gcOmega] = fhandleMeanGradParam(lfmKern2(1), t2, 'Pos');
[geAlpha, geOmega] = fhandleMeanGradParam(lfmKern2(1), t2, 'Vel');
[gradAlpha, gradOmega] = getLocalGradAlphaOmega(lfmKern2);
g2Pos = fhandleMeanGradSwitching(lfmKern2(1), t2, 'Pos');
h2Vel = fhandleMeanGradSwitching(lfmKern2(1), t2, 'Vel');
gsp1Kern = zeros(1, length(lfmKern1));
gsp2Kern = zeros(1, length(lfmKern1));
% This means that are only two switching points involved, t1
% and t0 appear in k, like k_{ff}(t1-t0,t'-t0) and
% k_{mf}(t1-t0,t'-t0). The other points appear in c1(t - t1)
% and e1(t - t1)
for k=1:length(lfmKern1)
    temp = fhandleKernGradSwitchingPosPos1(lfmKern1(k), lfmKern2(k), t1, lfmKern1(k).limit);
    gsp1Kern(k) = gsp1Kern(k) - sum(sum((temp*coeffPosPos.').*covGrad));
    temp = fhandleKernGradSwitchingPosPos2(lfmKern2(k), lfmKern1(k), lfmKern1(k).limit, t1).';
    gsp1Kern(k) = gsp1Kern(k) - sum(sum((temp*coeffPosPos.').*covGrad));
    gsp2Kern(k) = gsp2Kern(k) + sum(sum((temp*coeffPosPos.').*covGrad));
    temp = fhandleKernGradSwitchingVelPos1(lfmKern1(k), lfmKern2(k), t1, lfmKern1(k).limit);
    gsp1Kern(k) = gsp1Kern(k) - sum(sum((temp*coeffPosVel.').*covGrad));
    temp = fhandleKernGradSwitchingVelPos2(lfmKern2(k), lfmKern1(k), lfmKern1(k).limit, t1).';
    gsp1Kern(k) = gsp1Kern(k) - sum(sum((temp*coeffPosVel.').*covGrad));
    gsp2Kern(k) = gsp2Kern(k) + sum(sum((temp*coeffPosVel.').*covGrad));
end
if isempty(generalConst{i,j})
    matGradAlpha = PosPos*gcAlpha.' + PosVel*geAlpha.';
    matGradOmega = PosPos*gcOmega.' + PosVel*geOmega.';
    gCoeff = gradAlpha*sum(sum(matGradAlpha.*covGrad)) + ...
        gradOmega*(sum(sum(matGradOmega.*covGrad)));
    matGradSp2 = PosPos*g2Pos.' + PosVel*h2Vel.';
    gsp2 = - sum(sum(matGradSp2.*covGrad));
    g3(i) = gsp1Mean + sum(gsp1Kern);                % switching points 1
    g3(j) = gsp2Mean + sum(gsp2Kern) + gsp2;         % switching points 2
else
    constGradAlpha = generalConstGrad{1}{i,j};
    constGradOmega = generalConstGrad{2}{i,j};
    constVal = generalConst{i,j};
    % Derivative wrt parameters
    matGradAlpha = PosPos*((constVal(1,1)*gcAlpha.' + constGradAlpha(1,1)*c2.') ...
        + constVal(2,1)*geAlpha.' + constGradAlpha(2,1)*e2.') ...
        + PosVel*((constVal(1,2)*gcAlpha.' + constGradAlpha(1,2)*c2.') ...
        + constVal(2,2)*geAlpha.' + constGradAlpha(2,2)*e2.');

    matGradOmega = PosPos*((constVal(1,1)*gcOmega.' + constGradOmega(1,1)*c2.') ...
        + constVal(2,1)*geOmega.' + constGradOmega(2,1)*e2.') ...
        + PosVel*((constVal(1,2)*gcOmega.' + constGradOmega(1,2)*c2.') ...
        + constVal(2,2)*geOmega.' + constGradOmega(2,2)*e2.');

    gCoeff = gradAlpha*sum(sum(matGradAlpha.*covGrad)) + ...
        gradOmega*(sum(sum(matGradOmega.*covGrad)));
    % Derivative wrt switching points
    % Firts, notice that gsp2Kern doesn't correspond to the switching point
    % of the interval j (or say 2). It's just the derivative of the inner
    % switching point that lead to the actual switching point for the
    % interval j. So we rename it first.
    gspInt = sum(gsp2Kern);
    % Compute the derivative of the swicthing point j, not appearing in the
    % constant
    matGradSp2 = PosPos*(constVal(1,1)*g2Pos.' + constVal(2,1)*h2Vel.') + ...
        PosVel*(constVal(1,2)*g2Pos.' + constVal(2,2)*h2Vel.');
    gsp2 = - sum(sum(matGradSp2.*covGrad));
    % Compute the derivatives for all the switching points apearing in the
    % constant
    constGradSPoint = generalConstGrad{3}{i,j};
    numberSP = size(constGradSPoint,2);
    gspInBetween = zeros(1, numberSP);
    for k=1:numberSP
        temp = PosPos*(constGradSPoint(1,k)*c2.' + constGradSPoint(3,k)*e2.') + ...
            PosVel*(constGradSPoint(2,k)*c2.' + constGradSPoint(4,k)*e2.');
        gspInBetween(k) = sum(sum(temp.*covGrad));
    end
    % Assign derivatives wrt all other switching points, with a correction
    % for the derivative of the innermost switching point in the constant
    gspInBetween(end) =  gspInBetween(end) + gspInt;
    g3(i) = gsp1Mean + sum(gsp1Kern);           % switching point 1
    g3(j) = gsp2Mean + gsp2 + gspInBetween(1);  % switching point 2
    g3(i+1:j-1) = fliplr(gspInBetween(2:end));
end
% Assign derivatives wrt first system
g1{1} = g1Mean + sum(g1Kern(:,1:3), 1);          % mass 1, spring 1, damper 1
g1{2} = g1Kern(:,4).';                            % inverse widths
g1{3} = g1Kern(:,5).';                            % sensitivities 1
% Assign derivatives wrt second system
g2{1} = g2Mean + sum(g2Kern(:,1:3), 1) + gCoeff; % mass 2, spring 2, damper 2
g2{2} = g2Kern(:,4).';                            % inverse widths
g2{3} = g2Kern(:,5).';                            % sensitivities 2


text:
Tsym=1.25e-6;
Fdmax=40;
Ps=1;
SampleTime=3e-6;
SampleTime2=3e-6;


text:
function [rp,nodenew,elemnew] = recoverP02P1(node,elem,p,recoverMethod)
%% recoverP
% Recover the pressure from barycenter to the node.
% The recover methods
% 1) 'LS': least square
% 2) 'LA': solve the Laplacian problem, only this method
%           the new node and new elem will be constructed.
%
% Created by Lin Zhong April, 2013.

if ~exist('recoveryMethod','var')
    recoverMethod = 'LA';
end
NT = size(elem,1);
N = size(node,1);
rp = zeros(N,1);
% baycenter
xnode = 1/3*(node(elem(:,1),:) + node(elem(:,2),:) + node(elem(:,3),:));
NT2N = sparse(repmat(1:NT,1,3),[elem(:,1);elem(:,2);elem(:,3)],true,NT,N);

if strcmp(recoverMethod,'LS')
    % least square fit for every node
    for i = 1:N
        tempx = xnode(NT2N(:,i),:);
        tempp = p(NT2N(:,i));
        tempn = size(tempx,1);
        if tempn == 1
            rp(i) = tempp;
        elseif tempn == 2
            rp(i) = mean(tempp);
        else
            X = ones(tempn,3);
            X(:,1:2) = tempx;
            coefficient = (X'*X)\(X'*tempp);
            rp(i) = node(i,:)*coefficient(1:2) +coefficient(3);
        end
    end
end

if strcmp(recoverMethod,'LA')
    T = auxstructure(elem);
    edge = T.edge;
    edge2elem = T.edge2elem;
    %% find the boundary edge
    isbdEdge = (edge2elem(:,1) == edge2elem(:,2));
    freeEdgeidx = find(~isbdEdge);
    bdEdge = edge(isbdEdge,:);
    freeEdge = edge(freeEdgeidx,:);

    NE = size(edge,1);
    NEbd = size(bdEdge,1);
    NEin = NE - NEbd;
    NTnew = 2*NEin;
    elemnew = zeros(NTnew,3);
    nodenew = [node; xnode];
    %% find the corner point
    acubdNode = accumarray([bdEdge(:,1);bdEdge(:,2)], 1, [N 1]);
    bdNodeIdx = find(acubdNode);
    freeNodeIdx = find(~acubdNode);
    bdNode = node(bdNodeIdx);
    NbdNode = size(bdNode,1);

    % new elems corresponding to the boundary edges
%     elemnew(1:NEbd,:) = [bdEdge(:,1) bdEdge(:,2) N+edge2elem(isbdEdge,1)];

    % new elems corresponding to the interior edges
    elemnew(1:NEin,:) = [freeEdge(:,1) N+edge2elem(freeEdgeidx,1) N+edge2elem(freeEdgeidx,2)];
    elemnew(NEin+1:end,:) = [freeEdge(:,2) N+edge2elem(freeEdgeidx,1) N+edge2elem(freeEdgeidx,2)];

%   A = assemblematrix(nodenew,elemnew,1);
    %------- construc stiffness matrix-------------
    A11 = sparse(N,N);
    A12 = sparse(N,NT);
    [Dlambda,area] = gradbasis(nodenew,elemnew);
    % the original node
    ii = double(elemnew(:,1));
    DiDi =area.*dot(Dlambda(:,:,1),Dlambda(:,:,1),2);
    A11 = A11 +sparse(ii,ii,DiDi,N,N);
    % the original node and the new nodes
    for j = 2:3
        jj = double(elemnew(:,j)-N);
        DiDj =area.*dot(Dlambda(:,:,1),Dlambda(:,:,j),2);
        A12 = A12 +sparse(ii,jj,DiDj,N,NT);
    end
    A = [A11 A12];
    %------------------------------------------
    rp(freeNodeIdx) = (-A(freeNodeIdx,N+1:end)*p)./diag(A(freeNodeIdx,freeNodeIdx));

    % use least square to recover the boundary nodes
    for i = 1:NbdNode
        ii = bdNodeIdx(i);
        tempx = xnode(NT2N(:,ii),:);
        tempp = p(NT2N(:,ii));
        tempn = size(tempx,1);
        if tempn == 1
            rp(ii) = tempp;
        elseif tempn == 2
            rp(ii) = mean(tempp);
        else
            X = ones(tempn,3);
            X(:,1:2) = tempx;
            coefficient = (X'*X)\(X'*tempp);
            rp(ii) = node(ii,:)*coefficient(1:2) +coefficient(3);
        end
    end
end
end


text:
function f = fliplr( f )
%FLIPLR   Flip/reverse a SEPARABLEAPPROX in the x-direction.
%   G = FLIPLR( F ) returns a SEPARABLEAPPROX G with the same domain as F but reversed;
%   that is, G(x,y) = F(a+b-x,y), where the domain is [a, b, c, d].
%
% See also FLIPUD.

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

% Empty check:
if ( isempty( f ) )
    return
end

% Flip the row slices:
f.rows = flipud( f.rows );

end


text:
function spline_test235 ( )

%*****************************************************************************80
%
%% TEST235 tests SPLINE_PCHIP_SET and SPLINE_PCHIP_VAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 August 2005
%
%  Author
%
%    John Burkardt
%
  n = 21;
  ne = 101;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST235\n' );
  fprintf ( 1, '  SPLINE_PCHIP_SET sets up a piecewise cubic \n' );
  fprintf ( 1, '    Hermite interpolant.\n' );
  fprintf ( 1, '  SPLINE_PCHIP_VAL evaluates the interpolant.\n' );
  fprintf ( 1, '\n' );
%
%  Compute Runge's function at N points in [-1,1].
%
  for i = 1 : n
    x(i) = -1.0 + ( i - 1 ) / 10.0;
    f(i) = frunge ( x(i) );
  end
%
%  SPLINE_PCHIP_SET takes the data in X and F, and constructs a table in D
%  that defines the interpolant.
%
  d = spline_pchip_set ( n, x, f );
%
%  Evaluate the interpolant and derivative at NE points from -1 to 0.
%
  for i = 1 : ne
    xe(i) = -1.0 + ( i - 1 ) / ( ne - 1 );
  end

  fe = spline_pchip_val ( n, x, f, d, ne, xe );
%
%  Print the table of X, F(exact) and F(interpolated)
%
  for i = 1 : ne
    diff = fe(i) - frunge ( xe(i) );
    fprintf ( 1, '  %8f  %10f  %10f  %14e\n', ...
      xe(i), frunge ( xe(i) ), fe(i), diff );
  end

  return
end


text:
function pf = rdivide(arg1,arg2)
% implements pf1 ./ b and a ./ pf2
%
% overload the .* operator, i.e. one can now write x .* pf in order to
% scale the @PoleFigure pf by the factor x
%
% See also
% PoleFigure/PoleFigure PoleFigure/plus PoleFigure/minus

if isa(arg2,'double')
  arg2 = 1./arg2;
else
  for i = 1:arg2.numPF, arg2.allI{i} = 1./arg2.allI{i};end
end

pf = arg1 .* arg2;

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
%% demo.m
%
% Shows a couple of sample registrations using
% ICP - Iterative Closest Point
%
% Jakob Wilm and Martin Kjer, Technical University of Denmark, 2012

addpath('../../Utils');

m = 80; % width of grid
n = m^2; % number of points

[X,Y] = meshgrid(linspace(-2,2,m), linspace(-2,2,m));
X = reshape(X,1,[]);
Y = reshape(Y,1,[]);
Z = sin(X).*cos(Y);

% Create the data point-matrix
D = [X; Y; Z];

% Translation values (a.u.):
Tx = 0.5;
Ty = -0.3;
Tz = 0.2;

% Translation vector
T = [Tx; Ty; Tz];

% Rotation values (rad.):
rx = 0.3;
ry = -0.2;
rz = 0.05;

Rx = [1 0 0;
      0 cos(rx) -sin(rx);
      0 sin(rx) cos(rx)];

Ry = [cos(ry) 0 sin(ry);
      0 1 0;
      -sin(ry) 0 cos(ry)];

Rz = [cos(rz) -sin(rz) 0;
      sin(rz) cos(rz) 0;
      0 0 1];

% Rotation matrix
R = Rx*Ry*Rz;

% Transform data-matrix plus noise into model-matrix
M = R * D + repmat(T, 1, n);

% Add noise to model and data
rng(2912673);
M = M + 0.01*randn(3,n);
D = D + 0.01*randn(3,n);

%% Run ICP (standard settings)
[Ricp Ticp ER t] = icp(M, D, 15);

% Transform data-matrix using ICP result
Dicp = Ricp * D + repmat(Ticp, 1, n);

% Plot model points blue and transformed points red
figure;
subplot(2,2,1);
plot3(M(1,:),M(2,:),M(3,:),'bo',D(1,:),D(2,:),D(3,:),'r.');
axis equal;
xlabel('x'); ylabel('y'); zlabel('z');
title('Red: z=sin(x)*cos(y), blue: transformed point cloud');

% Plot the results
subplot(2,2,2);
plot3(M(1,:),M(2,:),M(3,:),'bo',Dicp(1,:),Dicp(2,:),Dicp(3,:),'r.');
axis equal;
xlabel('x'); ylabel('y'); zlabel('z');
title('ICP result');

% Plot RMS curve
subplot(2,2,[3 4]);
plot(0:15,ER,'--x');
xlabel('iteration#');
ylabel('d_{RMS}');
legend('bruteForce matching');
title(['Total elapsed time: ' num2str(t(end),2) ' s']);

%% Run ICP (fast kDtree matching and extrapolation)
[Ricp Ticp ER t] = icp(M, D, 15, 'Matching', 'kDtree', 'Extrapolation', true);

% Transform data-matrix using ICP result
Dicp = Ricp * D + repmat(Ticp, 1, n);

% Plot model points blue and transformed points red
figure;
subplot(2,2,1);
plot3(M(1,:),M(2,:),M(3,:),'bo',D(1,:),D(2,:),D(3,:),'r.');
axis equal;
xlabel('x'); ylabel('y'); zlabel('z');
title('Red: z=sin(x)*cos(y), blue: transformed point cloud');

% Plot the results
subplot(2,2,2);
plot3(M(1,:),M(2,:),M(3,:),'bo',Dicp(1,:),Dicp(2,:),Dicp(3,:),'r.');
axis equal;
xlabel('x'); ylabel('y'); zlabel('z');
title('ICP result');

% Plot RMS curve
subplot(2,2,[3 4]);
plot(0:15,ER,'--x');
xlabel('iteration#');
ylabel('d_{RMS}');
legend('kDtree matching and extrapolation');
title(['Total elapsed time: ' num2str(t(end),2) ' s']);

%% Run ICP (partial data)

% Partial model point cloud
Mp = M(:,Y>=0);

% Boundary of partial model point cloud
b = (abs(X(Y>=0)) == 2) | (Y(Y>=0) == min(Y(Y>=0))) | (Y(Y>=0) == max(Y(Y>=0)));
bound = find(b);

% Partial data point cloud
Dp = D(:,X>=0);

[Ricp Ticp ER t] = icp(Mp, Dp, 50, 'EdgeRejection', true, 'Boundary', bound, 'Matching', 'kDtree');

% Transform data-matrix using ICP result
Dicp = Ricp * Dp + repmat(Ticp, 1, size(Dp,2));

% Plot model points blue and transformed points red
figure;
subplot(2,2,1);
plot3(Mp(1,not(b)),Mp(2,not(b)),Mp(3,not(b)),'bo',...
      Mp(1,b),Mp(2,b),Mp(3,b),'go',...
      Dp(1,:),Dp(2,:),Dp(3,:),'r.')
axis equal;
xlabel('x'); ylabel('y'); zlabel('z');
title('Red: z=sin(x)*cos(y), blue: transformed point cloud');

% Plot the results
subplot(2,2,2);
plot3(Mp(1,not(b)),Mp(2,not(b)),Mp(3,not(b)),'bo',...
      Mp(1,b),Mp(2,b),Mp(3,b),'go',...
      Dicp(1,:),Dicp(2,:),Dicp(3,:),'r.');
axis equal;
xlabel('x'); ylabel('y'); zlabel('z');
title('ICP result');

% Plot RMS curve
subplot(2,2,[3 4]);
plot(0:50,ER,'--x');
xlabel('iteration#');
ylabel('d_{RMS}');
legend('partial overlap');
title(['Total elapsed time: ' num2str(t(end),2) ' s']);

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function losses = loss_func(o)
% Compute the loss associated with the intersection over union
% overlap between a ground-truth bounding box and any other
% windows.
%   losses = loss_func(o)
%
% Return value
%   losses    Loss for each element in the input
%
% Argument
%   o         Vector of overlap values

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

% The PASCAL VOC detection task loss
% Loss is 0 for IoU >= 0.5
% Loss is 1 for IoU < 0.5
losses = zeros(size(o));
I = find(o < 0.5);
losses(I) = 1.0;

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function str = showmydoc(item, indent)

% function str = showmydoc(item, indent)
% Display help text for a cfg_repeat and all of its options, without
% recursive calls to child nodes.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: showmydoc.m 2085 2008-09-12 10:26:59Z volkmar $

rev = '$Rev: 2085 $'; %#ok
str1 = showdoc(item.cfg_item, indent);
citems = subsref(item, substruct('.','values'));
if isfinite(item.num(2))
    str2{1} = sprintf('Between %d and %d options must be selected from:', ...
                         item.num(1), item.num(2));
else
    str2{1} = sprintf('%d or more options must be selected from:', ...
                         item.num(1));
end;
% Display short listing of cfg_repeat value items first
str3 = cell(1,numel(citems));
for k = 1:numel(citems)
    str3{k} = sprintf('* %s', subsref(citems{k}, substruct('.','name')));
end;
valitem = subsref(item, substruct('.','val'));
if ~isempty(valitem)
    str4 = cell(1,numel(valitem)+1);
    str4{1} = sprintf('Currently selected options:');
    for k = 1:numel(valitem)
        str4{k+1} = sprintf('* "%s"', subsref(valitem{k}, substruct('.','name')));
    end;
else
    str4 = {};
end;
str = [str1(:); {''}; str2(:); str3(:); str4(:)];

author: spm


text:
function obj = replace_basis_set(obj, condition_num, xBF_hires)
% Replace a basis set in an fmri_model object with another one of your
% choosing.
%
% This allows one to use a custom basis set, and also to use different
% basis sets for different trial types.
%
% Each condition across all sessions must be modeled with the same basis
% set. That is, there can be only one basis set per condition, e.g., one
% for anticipation (used in each session) and one for pain.
%
% :Usage:
% ::
%
%     obj = replace_basis_set(obj, condition_num, xBF_hires)
%
% :Examples:
% ::
%
%    % generate a custom spline basis set and use that for Condition 1,
%    % and the standard one for Condition 2:
%
%    [xBF_hires, xBF] = fmri_spline_basis(2, 'length', 12, 'nbasis', 3, 'order', 3, 'plot');
%
%    %save this to get info that is not typically in basis set until after
%    %model is built.

oldBF = obj.xBF(1);

% SPM adds these things here, so we will too, for consistency
xBF_hires.T = oldBF.T;
xBF_hires.T0 = oldBF.T0;
xBF_hires.UNITS = oldBF.UNITS;
xBF_hires.Volterra = oldBF.Volterra;

obj.xBF(condition_num) = xBF_hires;

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%% COMPARISON OF CONVERGENCE RATE OF VARIOUS FINITE ELEMENT METHODS
%
% This example is to compare the rate of convergence of various finite
% element approximation of the Poisson equation on the unit square:
%
% $$- \Delta u = f \; \hbox{in } (0,1)^2$$
%
% for the Dirichlet boundary condition. Results holds for other boundary
% conditions.

%% Problem
clear all; close all; clc;
[node,elem] = squaremesh([0,1,0,1],0.25);
pde = sincosdata;
bdFlag = setboundary(node,elem,'Dirichlet');
% pde = sincosNeumanndata;
% bdFlag = setboundary(node,elem,'Neumann');
option.L0 = 2;
option.maxIt = 5;
option.printlevel = 1;
option.plotflag = 0;

%% P1 element
err = femPoisson(node,elem,pde,bdFlag,option);

%%
% The optimal rate of convergence of the H1-norm (1st order), L2-norm (2nd
% order), and discrte maximum norm is observed. The 2nd order convergent
% rate between two discrete functions ||DuI-Duh|| is known as
% superconvergence.

%% CR nonconforming P1 element
option.elemType = 'CR';
err = femPoisson(node,elem,pde,bdFlag,option);

%%
% The optimal rate of convergence of the H1-norm (1st order), L2-norm (2nd
% order), and discrte maximum norm is observed. But no superconvergence of
% ||DuI-Duh||.

%% P2 element
option.elemType = 'P2';
err = femPoisson(node,elem,pde,bdFlag,option);

%% P3 element
option.elemType = 'P3';
err = femPoisson(node,elem,pde,bdFlag,option);

%% Conclusion

repo: Aero-matlab


text:
function test_bug2639
% DEPENDENCY ft_checkdata

% MEM 2gb
% WALLTIME 00:10:00

shufflechan = [1 3 2]';

channelcmb = {
  '1' '2'
  '1' '3'
  '2' '3'
  };

%%
freq1o           = [];
freq1o.freq      = 1;
freq1o.label     = {'1', '2', '3'}';
freq1o.dimord    = 'chan_freq';
freq1o.powspctrm = reshape([1 2 3], [3 1]);

freq1r           = freq1o;
freq1r.label     = freq1o.label(shufflechan);
freq1r.powspctrm = freq1o.powspctrm(shufflechan,:);

freq2o = ft_checkdata(freq1o);
freq2r = ft_checkdata(freq1r);
[selo, selr] = match_str(freq2o.label, freq2r.label);
assert(isequal(freq2o.powspctrm(selo,:), freq2r.powspctrm(selr,:)));

%%
freq1o           = [];
freq1o.freq      = 1;
freq1o.cumtapcnt = 1;
freq1o.label     = {'1', '2', '3'}';
freq1o.dimord    = 'rpt_chan_freq';
freq1o.fourierspctrm = reshape([1 2 3], [1 3 1]);

freq1r           = freq1o;
freq1r.label     = freq1o.label(shufflechan);
freq1r.fourierspctrm = freq1o.fourierspctrm(:,shufflechan,:);

freq2o = ft_checkdata(freq1o);
freq2r = ft_checkdata(freq1r);
[selo, selr] = match_str(freq2o.label, freq2r.label);
assert(isequal(freq2o.fourierspctrm(:,selo,:), freq2r.fourierspctrm(:,selr,:)));

%% full, sparse, fourier, sparsewithpow, fullfast

freq2o = ft_checkdata(freq1o, 'cmbstyle', 'fourier');
freq2r = ft_checkdata(freq1r, 'cmbstyle', 'fourier');
[selo, selr] = match_str(freq2o.label, freq2r.label);
assert(isequal(freq2o.fourierspctrm(:,selo,:), freq2r.fourierspctrm(:,selr,:)));

assert( isequal(freq1r.label,freq2r.label));
assert(~isequal(freq2o.label,freq2r.label));

freq2o = ft_checkdata(freq1o, 'cmbstyle', 'full');
freq2r = ft_checkdata(freq1r, 'cmbstyle', 'full');
[selo, selr] = match_str(freq2o.label, freq2r.label);
assert(isequal(freq2o.crsspctrm(selo,selo), freq2r.crsspctrm(selr,selr)));

assert( isequal(freq1r.label,freq2r.label));
assert(~isequal(freq2o.label,freq2r.label));

freq2o = ft_checkdata(freq1o, 'cmbstyle', 'fullfast');
freq2r = ft_checkdata(freq1r, 'cmbstyle', 'fullfast');
[selo, selr] = match_str(freq2o.label, freq2r.label);
assert(isequal(freq2o.crsspctrm(selo,selo), freq2r.crsspctrm(selr,selr)));

assert( isequal(freq1r.label,freq2r.label));
assert(~isequal(freq2o.label,freq2r.label));

freq2o = ft_checkdata(freq1o, 'cmbstyle', 'sparse', 'channelcmb', channelcmb);
freq2r = ft_checkdata(freq1r, 'cmbstyle', 'sparse', 'channelcmb', channelcmb);
[selo, selr] = match_strcmb(freq2o.labelcmb, freq2r.labelcmb);
assert(isequal(freq2o.crsspctrm(selo,:), freq2r.crsspctrm(selr,:)));

assert( isequal(freq2o.labelcmb,freq2r.labelcmb));

freq2o = ft_checkdata(freq1o, 'cmbstyle', 'sparsewithpow', 'channelcmb', channelcmb);
freq2r = ft_checkdata(freq1r, 'cmbstyle', 'sparsewithpow', 'channelcmb', channelcmb);
[selo, selr] = match_strcmb(freq2o.labelcmb, freq2r.labelcmb);
assert(isequal(freq2o.crsspctrm(selo,:), freq2r.crsspctrm(selr,:)));
[selo, selr] = match_str(freq2o.label, freq2r.label);
assert(isequal(freq2o.powspctrm(selo,:), freq2r.powspctrm(selr,:)));

assert( isequal(freq1r.label,freq2r.label));
assert(~isequal(freq2o.label,freq2r.label));
assert( isequal(freq2o.labelcmb,freq2r.labelcmb));

%%
freq1o           = [];
freq1o.freq      = 1;
freq1o.cumtapcnt = 1;
freq1o.dimord    = 'chancmb_freq';
freq1o.label     = {'1', '2', '3'}';
freq1o.powspctrm = [1 2 3]';
freq1o.labelcmb  = {
  '1' '2'
  '1' '3'
  '2' '3'
  };
freq1o.crsspctrm = [
  2
  3
  6
  ];

freq1r           = freq1o;
freq1r.label     = freq1o.label(shufflechan);
freq1r.labelcmb  = freq1o.labelcmb(shufflechan,:);
freq1r.powspctrm = freq1o.powspctrm(shufflechan,:);
freq1r.crsspctrm = freq1o.crsspctrm(shufflechan,:);

freq2o = ft_checkdata(freq1o);
freq2r = ft_checkdata(freq1r);
[selo, selr] = match_str(freq2o.label, freq2r.label);
assert(isequal(freq2o.powspctrm(selo,:), freq2r.powspctrm(selr,:)));
[selo, selr] = match_strcmb(freq2o.labelcmb, freq2r.labelcmb);
assert(isequal(freq2o.crsspctrm(selo,:), freq2r.crsspctrm(selr,:)));

freq2o = ft_checkdata(freq1o, 'cmbstyle', 'full');
freq2r = ft_checkdata(freq1r, 'cmbstyle', 'full');
[selo, selr] = match_str(freq2o.label, freq2r.label);
assert(isequal(freq2o.crsspctrm(selo,selo), freq2r.crsspctrm(selr,selr)));

% assert( isequal(freq1r.label,freq2r.label)); % this is where the channel reordering becomes clear
% assert(~isequal(freq2o.label,freq2r.label)); % this is where the channel reordering becomes clear

freq2o = ft_checkdata(freq1o, 'cmbstyle', 'fullfast');
freq2r = ft_checkdata(freq1r, 'cmbstyle', 'fullfast');
[selo, selr] = match_str(freq2o.label, freq2r.label);
assert(isequal(freq2o.crsspctrm(selo,selo), freq2r.crsspctrm(selr,selr))); % this is where another problem appears

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [sel1, sel2] = match_strcmb(chancmb1, chancmb2)
for i=1:size(chancmb1,1)
  chan1{i} = sprintf('%s_%s', chancmb1{i,:});
end
for i=1:size(chancmb2,1)
  chan2{i} = sprintf('%s_%s', chancmb2{i,:});
end
[sel1, sel2] = match_str(chan1, chan2);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip