text:
% POP_LIMO - prepare and convert EEGLAB data and structure to be
%              processed by LIMO.
%
% Usage:
%   >> STUDY = pop_limo( STUDY, ALLEEG, 'key', val)
%   >> STUDY = pop_limo( STUDY, ALLEEG, 'dat', 'key', val)
%   >> STUDY = pop_limo( STUDY, ALLEEG, 'ica', 'key', val)
%
% Inputs:
%   STUDY        - an EEGLAB STUDY set of loaded EEG structures
%   ALLEEG       - ALLEEG vector of one or more loaded EEG dataset structures
%   'dat'|'ica'  - show the interface for data channels or for ICA. The
%                  default is to use data.
%
% Optional inputs: same as STD_LIMO
%
% Graphical interface:
%   "See GLM factors" - [push button] See all the GLM factors or columns in
%                   the design matrix. This list depends on your study
%                   design. Adding contrast will change this list. In
%                   practice, it is important to understand that every
%                   single factor will be fitted by the GLM, and the choice
%                   of these factors is critical. It is recommended to
%                   include as many factors as possible, especially because
%                   this allows you to define new contrast later on without
%                   having to recompute the GLM on each subject (as long as
%                   these new contrast do not involve new factors).
%
%  "Interaction model for categorical indep. var." - When using more than
%                   one categorical variable, clicking this option forces
%                   to have factors which are the conjunction of the
%                   different independent var. values. This is useful only
%                   if you want to calculate interactions at the subject.
%                   The 'best' option is typically to have a design with
%                   all conditions (no factors) and create factors at the
%                   group level.
%
%  "Split regressions (continuous indep. var.)" - This options allows to
%                   split continuous variables for the different
%                   categorical variables. This is useful to compute
%                   interaction between continuous and categorical
%                   variables.
%
%  "Input data to use for the GLM" - [pop up menu] measure to use as input
%                   for the GLM. Currently, only "ERP" and "spectrum" are
%                   supported.
%
%  "Optimization method" - [pop up menu] may be Ordinary Least Squares (OLS),
%                   Weighted Least Squares (WLS), or Iterative Reweighted Least
%                   Squares. WTS should be used as it is more robust to trials
%                   with different time course. OLS is a standard solution and
%                   WLS/IRLS are solution that automatically weight down outliers
%                   (trials/data point).
%
%  "Options"        - [edit box] additional options. These are given
%                   directly as input to the std_limo.m function. They may
%                   be ['freqlim', value] frequency trimming in Hz or
%                   ['timelim', value] time trimming in millisecond. These
%                   allow to speed up computation by not including all the
%                   data as input. Example: 'timelim', [-50 650] only
%                   include data within the -50 ms and 650 ms for ERPs.
%
% "Erase previous model" - [checkbox] erase previous model for this
%                   measure.
%
% Outputs:
%   STUDY       - an EEGLAB STUDY set of loaded EEG structures
%
% Author: Arnaud Delorme, SCCN, UCSD, 2015-
%         Cyril Pernet, LIMO Team - edit info and defaults
%
% See also: STD_LIMO

% Copyright (C) Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [STUDY,com,limofiles] = pop_limo(STUDY, ALLEEG, measureflag, varargin)

if nargin < 2
    help pop_limo;
    return;
end

com = '';
if nargin > 2
    if strcmpi(measureflag,'dat')
        measureflagindx = 1;
    elseif strcmpi(measureflag,'ica')
        measureflagindx = 2;
    else
        varargin = { measureflag varargin{:} };
        measureflagindx = 1;
    end
else
    measureflagindx = 1;
end

% check that channel location are present
ALLEEG = eeg_checkset(ALLEEG, 'chanloc');

if nargin < 4
    dataMeasures = { 'ERP' 'Spectrum' 'ERSP'};
    fileMeasures = { 'daterp' 'datspec' 'dattimef'; 'icaerp' 'icaspec' 'icatimef'};
    methods      = { 'WLS' 'OLS' 'IRLS'};
    cb_measure   = [ 'if get(gcbo, ''value'') == 1,' ...
                     '   set(findobj(gcbf, ''tag'', ''options''), ''string'', '''');' ...
                     'else,' ...
                     '   set(findobj(gcbf, ''tag'', ''options''), ''string'', ''''''freqlim'''', [1 25]'');' ...
                     'end;' ];
    cb_listfactors = [ 'pop_listfactors(STUDY, ''gui'', ''on'', ' ...
                               '''level'', ''both'',' ...
                               '''splitreg''   , fastif(get(findobj(gcbf, ''tag'', ''splitreg''   ), ''value''), ''on'', ''off''),' ...
                               '''interaction'', fastif(get(findobj(gcbf, ''tag'', ''interaction''), ''value''), ''on'', ''off''));' ];
    uilist = { ...
        {'style' 'text'       'string' 'LInear MOdeling of EEG data' 'fontweight' 'bold' 'fontsize', 12} ...
        {'style' 'pushbutton' 'string' 'See GLM variables' 'callback' cb_listfactors } ...
        {'style' 'checkbox'   'string' 'Interaction model for categorical indep. var.' 'value' 0 'tag' 'interaction' } ...
        {'style' 'checkbox'   'string' 'Split regressions (continuous indep. var.)' 'tag' 'splitreg' } {} ...
        {'style' 'text'       'string' 'Input data to use for the GLM' } ...
        {'style' 'popupmenu'  'string' dataMeasures 'tag' 'measure' 'callback' cb_measure} ...
        {'style' 'text'       'string' 'Estimation method' } ...
        {'style' 'popupmenu'  'string' methods 'tag' 'method' } ...
        {'style' 'text'       'string' 'Options' } ...
        {'style' 'edit'       'string' '' 'tag' 'options' } ...
        {'style' 'checkbox'   'string' 'Erase previous model' 'tag' 'erase' 'value' true } ...
        };

    cline = [1.1 0.8];
    geometry = { [1.6 1] 1 1 1 cline cline cline 1 };
    geomvert = [ 1 1 1     1 1 1     1     1 ];

    [out_param userdat tmp res] = inputgui( 'geometry' , geometry, 'uilist', uilist, 'geomvert', geomvert, ...
                                            'title', 'LInear MOdeling of EEG data -- pop_limo()', 'helpcom', 'pophelp(''pop_limo'');');
    if isempty(res), return; end
    opttmp  = eval( [ '{ ' res.options ' }' ]);
    if length(opttmp) > 0 && isnumeric(opttmp{1})
        error([ 'Wrong options. Options must be in the format' 10 '''key'', val. For example ''timelim'', [-100 600]' ]);
    end
    options = { 'method' methods{res.method} 'measure' fileMeasures{measureflagindx,res.measure} opttmp{:} ...
                'erase'       fastif(res.erase, 'on', 'off') ...
                'splitreg'    fastif(res.splitreg, 'on', 'off') ...
                'interaction' fastif(res.interaction, 'on', 'off') };
else
    options = varargin;
end

[STUDY,limofiles] = std_limo(STUDY, ALLEEG, options{:});
com = sprintf('pop_limo(STUDY, ALLEEG, %s);', vararg2str(options));

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
%   This source code is part of the graph optimization package
%   deveoped for the lectures of robotics2 at the University of Freiburg.
%
%     Copyright (c) 2007 Giorgio Grisetti, Gian Diego Tipaldi
%
%   It is licences under the Common Creative License,
%   Attribution-NonCommercial-ShareAlike 3.0
%
%   You are free:
%     - to Share - to copy, distribute and transmit the work
%     - to Remix - to adapt the work
%
%   Under the following conditions:
%
%     - Attribution. You must attribute the work in the manner specified
%       by the author or licensor (but not in any way that suggests that
%       they endorse you or your use of the work).
%
%     - Noncommercial. You may not use this work for commercial purposes.
%
%     - Share Alike. If you alter, transform, or build upon this work,
%       you may distribute the resulting work only under the same or
%       similar license to this one.
%
%   Any of the above conditions can be waived if you get permission
%   from the copyright holder.  Nothing in this license impairs or
%   restricts the author's moral rights.
%
%   This software is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied
%   warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
%   PURPOSE.

%ls-slam.m
%this file is released under the creative common license

%solves a graph-based slam problem via least squares
%vmeans: matrix containing the column vectors of the poses of the vertices
%	 the vertices are odrered such that vmeans[i] corresponds to the ith id
%eids:	 matrix containing the column vectors [idFrom, idTo]' of the ids of the vertices
%	 eids[k] corresponds to emeans[k] and einfs[k].
%emeans: matrix containing the column vectors of the poses of the edges
%einfs:  3d matrix containing the information matrices of the edges
%	 einfs(:,:,k) refers to the information matrix of the k-th edge.
%n:	 number of iterations
%newmeans: matrix containing the column vectors of the updated vertices positions

function newmeans = ls_slam(vmeans, eids, emeans, einfs, n)

for i = 1:n
    vmeans = linearize_and_solve(vmeans, eids, emeans, einfs);
end

newmeans = vmeans;

end

%computes the taylor expansion of the error function of the k_th edge
%vmeans: vertices positions
%eids:   edge ids
%emeans: edge means
%k:	 edge number
%e:	 e_k(x)
%A:	 d e_k(x) / d(x_i)
%B:	 d e_k(x) / d(x_j)
function [e, A, B] = linear_factors(vmeans, eids, emeans, k)
%extract the ids of the vertices connected by the kth edge
id_i = eids(1,k);
id_j = eids(2,k);
%extract the poses of the vertices and the mean of the edge
v_i = vmeans(:,id_i);
v_j = vmeans(:,id_j);
z_ij = emeans(:,k);

%compute the homoeneous transforms of the previous solutions
zt_ij = v2t(z_ij);
vt_i = v2t(v_i);
vt_j = v2t(v_j);

%compute the displacement between x_i and x_j
f_ij=(inv(vt_i) * vt_j);

%this below is too long to explain, to understand it derive it by hand
theta_i = v_i(3);
ti = v_i(1:2,1);
tj = v_j(1:2,1);
dt_ij = tj-ti;

si = sin(theta_i);
ci = cos(theta_i);

A= [-ci, -si, [-si, ci]*dt_ij; si, -ci, [-ci, -si]*dt_ij; 0, 0, -1 ];
B =[  ci, si, 0           ; -si, ci, 0            ; 0, 0, 1 ];

ztinv = inv(zt_ij);
e = t2v(ztinv * f_ij);
ztinv(1:2,3) = 0;
A = ztinv*A;
B = ztinv*B;
end

%linearizes and solves one time the ls-slam problem specified by the input
%vmeans:   vertices positions at the linearization point
%eids:     edge ids
%emeans:   edge means
%einfs:    edge information matrices
%newmeans: new solution computed from the initial guess in vmeans
function newmeans = linearize_and_solve(vmeans, eids, emeans, einfs)
disp('allocating workspace...');
% H and b are respectively the system matrix and the system vector
H = zeros(size(vmeans,2)*3);
b = zeros(size(vmeans,2)*3,1);

disp('linearizing');
% this loop constructs the global system by accumulating in H and b the contributions
% of all edges (see lecture)
for k = 1:size(eids,2),
    id_i = eids(1,k);
    id_j = eids(2,k);
    [e, A, B] = linear_factors(vmeans, eids, emeans,  k);
    omega = einfs(:,:,k);
    %compute the blocks of H^k
    b_i = -A' * omega * e;
    b_j = -B' * omega * e;
    H_ii=  A' * omega * A;
    H_ij=  A' * omega * B;
    H_jj=  B' * omega * B;

    %accumulate the blocks in H and b
    H((id_i-1)*3+1:id_i*3,(id_i-1)*3+1:id_i*3) = ...
        H((id_i-1)*3+1:id_i*3,(id_i-1)*3+1:id_i*3)+ H_ii;
    H((id_j-1)*3+1:id_j*3,(id_j-1)*3+1:id_j*3) = ...
        H((id_j-1)*3+1:id_j*3,(id_j-1)*3+1:id_j*3) + H_jj;
    H((id_i-1)*3+1:id_i*3,(id_j-1)*3+1:id_j*3) = ...
        H((id_i-1)*3+1:id_i*3,(id_j-1)*3+1:id_j*3) + H_ij;
    H((id_j-1)*3+1:id_j*3,(id_i-1)*3+1:id_i*3) = ...
        H((id_j-1)*3+1:id_j*3,(id_i-1)*3+1:id_i*3) + H_ij';
    b((id_i-1)*3+1:id_i*3,1) = ...
        b((id_i-1)*3+1:id_i*3,1) + b_i;
    b((id_j-1)*3+1:id_j*3,1) = ...
        b((id_j-1)*3+1:id_j*3,1) + b_j;

    %NOTE on Matlab compatibility: note that we use the += operator which is octave specific
    %using H=H+.... results in a tremendous overhead since the matrix would be entirely copied every time
    %and the matrix is huge
end;
disp('Done');
%note that the system (H b) is obtained only from
%relative constraints. H is not full rank.
%we solve the problem by anchoring the position of
%the the first vertex.
%this can be expressed by adding the equation
%  deltax(1:3,1)=0;
%which is equivalent to the following
H(1:3,1:3) = H(1:3,1:3) + eye(3);

SH = sparse(H);
disp('System size: '),disp(size(H));
disp('solving (may take some time) ...');
deltax = SH\b;
disp('Done! ');

%split the increments in nice 3x1 vectors and sum them up to the original matrix
newmeans = vmeans + reshape(deltax, 3, size(vmeans,2));

disp('Normalizing the angles');
%normalize the angles between -PI and PI
for i = 1:size(newmeans,2)
    s = sin(newmeans(3,i));
    c = cos(newmeans(3,i));
    newmeans(3,i) = atan2(s,c);
end
disp('Done');
end

repo: graphslam


text:
function rtk = detslip_LLI(rtk,obs,i,rcv)

global glc;
sat=obs(i).sat;

for f=1:rtk.opt.nf

    if obs(i).L(f)==0||abs(timediff(obs(i).time,rtk.sat(sat).pt(rcv,f)))<glc.DTTOL
        continue;
    end

    % restore previous LLI
    if rcv==1
        LLI=getbitu(rtk.sat(sat).slip(f),0,2);
    else
        LLI=getbitu(rtk.sat(sat).slip(f),2,2);
    end

    % detect slip by cycle slip flag in LLI
    if rtk.tt>=0
        slip=obs(i).LLI(f);
    else
        slip=LLI;
    end

    if (bitand(LLI,2)&&~bitand(obs(i).LLI(f),2))||(~bitand(LLI,2)&&bitand(obs(i).LLI(f),2))
        slip=bitor(slip,1);
    end

    if rcv==1
        rtk.sat(sat).slip(f)=setbitu(rtk.sat(sat).slip(f),0,2,obs(i).LLI(f));
    else
        rtk.sat(sat).slip(f)=setbitu(rtk.sat(sat).slip(f),2,2,obs(i).LLI(f));
    end

    rtk.sat(sat).slip(f)=bitor(rtk.sat(sat).slip(f),slip);
    if bitand(obs(i).LLI(f),2)
        rtk.sat(sat).half(f)=0;
    else
        rtk.sat(sat).half(f)=1;
    end

end

return

repo: GINav


text:
function H = displ_field_gradient_Nd(D, pix_resolution)
    Nd = size(D, 4);
    H = zeros([Nd,Nd, size(D, 1),size(D, 2), size(D, 3)]);
    for i = 1 : Nd
%         GS = cell(2,1);
        GS = {};
        if Nd == 2
            [GS{1}, GS{2}] = my_gradient( squeeze(D(:,:, 1, i)), pix_resolution, 0, 2);
        elseif Nd == 3
            [GS{1}, GS{2}, GS{3}] = my_gradient( squeeze(D(:,:,:, i)), pix_resolution, 0, 2);
        end

        for j = 1 : Nd
%             tmpD = imgaussfilt(D(:,:,i), 1.5);
%             tmpD = imgaussfilt(D(:,:,i), 1.9);
%             tmpD = D(:,:,i);
%             G = DGradient(tmpD, [], j) / pix_resolution(j);
              G = GS{j};
%             G = medfilt2(G, [3,3]);
%             G = imgaussfilt(G, 1.5);
            H(i, j, :, :, :) = G;
        end
    end
end

repo: pTVreg


text:
% test for TV projection on the TV ball

path(path, 'toolbox/');

clear options;
options.bound = 'sym';

rep = 'results/tv-projection/';
if not(exist(rep))
    mkdir(rep);
end
repiter = [rep 'iter/'];
if not(exist(repiter))
    mkdir(repiter);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% geometrical image
n = 256;
J = 4;
sigma = 100;
options.bound = 'per';
A = ones(n);
for j=0:J-1
    B = A;
    for k=1:2^j
        I = find(B==k);
        U = perform_blurring(randn(n),sigma,options);
        s = median(U(I));
        I1 = find( (B==k) & (U>s) );
        I2 = find( (B==k) & (U<=s) );
        A(I1) = 2*k-1;
        A(I2) = 2*k;
    end
    clf; imageplot(A);
    imwrite(rescale(A), [rep 'geometrical-' num2str(j) '.png'], 'png');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1D test
name = 'random';
n = 1024*4;
f = randn(n,1);
tv = compute_total_variation(f, options);

tv_divide = [2 4  8 16 32];
m = length(tv_divide);
clf;
imageplot(f, 'Original', m+1,1,1);
for i=1:m
    tau = tv/tv_divide(i);
    [f1,err_tv,err_l2] = perform_tv_projection(f,tau,options);
    tv1 = compute_total_variation(f1, options);
    imageplot(f1, ['TV/' num2str(tv_divide(i))], m+1,1,i+1);
end
saveas(gcf, [rep name '-1d-tv-projection.png'], 'png');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load image
name = 'boat';
name = 'noise';
n = 200;
M = rescale( crop( load_image(name), n ) );

sigma = 0; % smoothing
if sigma>0
    M = perform_blurring(M,sigma);
end

% reduce TV norm by projection
tv = compute_total_variation(M, options);
options.niter = 2000;

tv_divide = [2 4 8 16 32 64 128];
m = length(tv_divide);
a = 2; b = ceil((m+1)/2);
options.x = [];
clf;
if strcmp(name, 'noise')
    MM = clamp( 0.45*(M-mean(M(:)))/std(M(:)), -1,1);
else
    MM = M;
end
warning off;
imwrite( rescale(MM), [repiter name '-2d-tv-projection-sig' num2str(sigma) '-01.png'] );
warning on;
imageplot(M, 'Original', a,b,1);
for i=1:m
    tau = tv/tv_divide(i);
    [M1,err_tv,err_l2] = perform_tv_projection(M,tau,options);
    tv1 = compute_total_variation(M1, options);
    disp( ['Final TV error: ' num2str( (tv1-tau)/tv1 ) '.'] );
    options.x = M1;
    if strcmp(name, 'noise')
        MM = clamp( 0.45*(M1-mean(M1(:)))/std(M1(:)), -1,1);
    else
        MM = M1;
    end
    imageplot(MM, ['TV/' num2str(tv_divide(i))], a,b,i+1);
    warning off;
    imwrite( rescale(MM), [repiter name '-2d-tv-projection-sig' num2str(sigma) '-' num2string_fixeddigit(tv_divide(i),2) '.png'] );
    warning on;
end
saveas(gcf, [rep name '-2d-tv-projection-sig' num2str(sigma) '.png'], 'png');

author: gpeyre

repo: matlab-toolboxes


text:
%% IPOPT Install for OPTI Toolbox
% Copyright (C) 2012 Jonathan Currie (I2C2)

% This file will help you compile Interior Point OPTimizer (IPOPT) for use
% with MATLAB.

% My build platform:
% - Windows 8 x64
% - Visual Studio 2012
% - Intel Compiler XE (FORTRAN)
% - Intel Math Kernel Library

% NOTE - From OPTI v1.71 IPOPT is now dynamically linked against the
% MathWorks supplied libmwma57.dll (HSL MA57). This step is optional
% as we are also compiling MUMPS. Alternatively you can skip MUMPS,
% and just use MA57! Also be aware libmwma57.dll does not play well
% on unconstrained problems due in part to missing MeTiS, thus the
% ma57 pivot order option is overidden automatically.

% To recompile you will need to get / do the following:

% 1) Get IPOPT
% IPOPT is available from http://www.coin-or.org/download/source/Ipopt/.
% Download the latest version of the source.

% 2) Compile IPOPT
% The easiest way to compile IPOPT is to use the Visual Studio Project
% Builder included with OPTI. Use the following commands, substituting the
% required path on your computer:

% %% Visual Studio Builder Commands
% path = 'full path to IPOPT here'; %e.g. 'C:\Solvers\IPOPT'
% mumpspath = 'full path to MUMPS here'; %e.g. 'C:\Solvers\MUMPS'
% sdir = [path '\src'];
% hdrs = {[mumpspath '\include'], [mumpspath '\libseq'], [path '\..\BuildTools\headers']};
% name = 'libipopt';
% opts = [];
% opts.exPP = {'IPOPT_BUILD','FUNNY_MA57_FINT','COINHSL_HAS_MA57','_CRT_SECURE_NO_WARNINGS'};
% opts.exFolder = {'Apps','Algorithm\Inexact','contrib\LinearSolverLoader'};
% opts.exclude = 'AmplTNLP.cpp';
% VS_WriteProj(sdir,name,hdrs,opts)
% %%

% Once complete, you will have a directory called IPOPT\libipopt. Open the
% Visual Studio 2012 project file, then complete the following steps:
%   a) Open config_default.h and make the following changes:
%       - Add #define COIN_HAS_METIS 1 under the METIS library line
%       - Add #define COIN_HAS_MUMPS 1 under the MUMPS library line
%       - Comment the line #define COIN_HAS_HSL 1
%       - Comment the line #define COIN_HAS_ASL 1
%   b) Build a Win32 or x64 Release to compile the code.
%   c) Copy the generated .lib file to the following folder:
%
%   OPTI/Solvers/ipopt/Source/lib/win32 or win64
%
%   You will also need to copy all header files from IPOPT/src to
%   the following folder (inc sub-directories):
%
%   OPTI/Solvers/ipopt/Source/Include

% 3) Complete MUMPS Compilation
% Supplied with this toolbox is the MUMPS linear system solver and METIS
% libraries. These are required when linking IPOPT so make sure you have
% built these first! See opti_MUMPS_Install.m in Solvers/mumps for more
% information.

% 4) libmwma57 Compilation
% As MathWorks only supply a dll, I have manually created an import library
% for libmwma57.dll, so we can link against it. The project and source for
% this is included in libmwma57.zip. It is optional to recompile this, as I
% have included compiled libraries for you!

% 5) IPOPT MEX Interface
% Supplied with the IPOPT source is a MEX interface by Dr. Peter Carbonetto.
% This is in Ipopt/contrib/MatlabInterface/src. However I have modified
% parts of this interface, so you must use the version supplied with OPTI.

% 6) Compile the MEX File
% The code below will automatically include all required libraries and
% directories to build the IPOPT MEX file. Once you have completed all the
% above steps, simply run this file to compile IPOPT! You MUST BE in the
% base directory of OPTI!

clear ipopt

% Modify below function if it cannot find Intel MKL on your system.
[mkl_link,mkl_for_link] = opti_FindMKL();
% Get Arch Dependent Library Path
libdir = opti_GetLibPath();
% Dependency Paths
mumpsdir = '..\..\mumps\Source\';

fprintf('\n------------------------------------------------\n');
fprintf('IPOPT MEX FILE INSTALL\n\n');

%Get IPOPT Libraries
post = ' -IInclude -IInclude/Common -IInclude/BuildTools -IInclude/Interfaces -IInclude/LinAlg -IInclude/LinAlg/TMatrices -IInclude/Algorithm ';
post = [post ' -L' libdir ' -llibIpopt -llibmwma57 -DhaveMA57'];
%Get MUMPS Libraries
post = [post ' -I' mumpsdir '\Include -L' mumpsdir libdir '-llibdmumps_c -llibdmumps_f -llibseq_c -llibseq_f -llibmetis -llibpord'];
%Get Intel Fortran Libraries (for MUMPS build) & MKL Libraries (for BLAS)
post = [post mkl_link mkl_for_link];

%Common Args
post = [post ' -DIPOPT_BUILD -output ipopt'];

%CD to Source Directory
cdir = cd;
cd 'Solvers/ipopt/Source';

%Compile & Move
pre = 'mex -v -largeArrayDims matlabexception.cpp matlabfunctionhandle.cpp matlabjournal.cpp iterate.cpp ipoptoptions.cpp options.cpp sparsematrix.cpp callbackfunctions.cpp matlabinfo.cpp matlabprogram.cpp ipopt.cpp';
try
    eval([pre post])
    movefile(['ipopt.' mexext],'../','f')
    fprintf('Done!\n');
catch ME
    cd(cdir);
    error('opti:ipopt','Error Compiling IPOPT!\n%s',ME.message);
end
cd(cdir);
fprintf('------------------------------------------------\n');

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function this = set_from_affine_geometry(this, affineTransformation, nVoxels, TR_s)
% Sets dimInfo from affine affineTransformation, assuming 4D nifti data
%
%   Y = MrDimInfo()
%   Y.set_from_affine_geometry(affineTransformation, nVoxels, TR_s)
%
% This is a method of class MrDimInfo.
%
% IN
%
% OUT
%
% EXAMPLE
%   set_from_affine_geometry
%
%   See also MrDimInfo

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

if nargin < 3
    nVoxels = [1 1 1 1];
end

if nargin < 4
    TR_s = 1;
end

% or add them...
dimLabelsGeom = {'x','y','z', 't'};
units = {'mm', 'mm', 'mm', 's'};
iDimGeom = 1:4;
% update existing geom dimensions, add new ones for
% non-existing
iValidDimLabels = this.get_dim_index(dimLabelsGeom);
iDimGeomExisting = find(iValidDimLabels);
iDimGeomAdd = setdiff(iDimGeom, iDimGeomExisting);

% need nifti to reference first sampling point as offcenter
resolutions = [affineTransformation.scaling TR_s];

% voxel position by voxel center, time starts at 0
firstSamplingPoint = [affineTransformation.scaling 0]/2;

% if dimension labels exist, just update values
this.set_dims(dimLabelsGeom(iDimGeomExisting), ...
    'resolutions', resolutions(iDimGeomExisting), ...
    'nSamples', nVoxels(iDimGeomExisting), ...
    'firstSamplingPoint', firstSamplingPoint(iDimGeomExisting), ...
    'units', units(iDimGeomExisting));

% if they do not exist, create dims
this.add_dims(dimLabelsGeom(iDimGeomAdd), ...
    'resolutions', resolutions(iDimGeomAdd), ...
    'nSamples', nVoxels(iDimGeomAdd), ...
    'firstSamplingPoint', firstSamplingPoint(iDimGeomAdd), ...
    'units', units(iDimGeomAdd));

author: translationalneuromodeling

repo: tapas


text:
function [easting,northing]=latlon2eastingsNorthings(sourcelat, sourcelon, lat, lon)
% Convert lat,lon to eastings,northings with sourcelat, sourcelon as origin
% inputs are in degrees, output is in metres relative to origin
% [easting,northing]=latlon2eastingsNorthings(sourcelat, sourcelon, lat, lon)
deg2m = deg2km(1) * 1000;
for c=1:length(lat)
    e = distance(lat(c), lon(c), lat(c), sourcelon) * deg2m;
    easting(c) = e * sign(lon(c)-sourcelon);
    n = distance(lat(c), lon(c), sourcelat, lon(c)) * deg2m;
    northing(c) = n * sign(lat(c)-sourcelat);
end

author: geoscience-community-codes

repo: GISMO


text:
function Mn = powermanifold(M, n)
% Returns a structure describing a power manifold M^n = M x M x ... x M.
%
% function Mn = powermanifold(M, n)
%
% Input: a manifold structure M and an integer n >= 1.
%
% Output: a manifold structure Mn representing M x ... x M (n copies of M)
% with the metric of M extended element-wise. Points and vectors are stored
% as cells of size nx1.
%
% This code is for prototyping uses. The structures returned are often
% inefficient representations of power manifolds owing to their use of
% for-loops, but they should allow to rapidly try out an idea.
%
% Example (an inefficient representation of the oblique manifold (3, 10)):
% Mn = powermanifold(spherefactory(3), 10)
% disp(Mn.name());
% x = Mn.rand()
%
% See also: productmanifold

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%   NB, July 4, 2013: Added support for vec, mat, tangent.
%                     Added support for egrad2rgrad and ehess2rhess.

    assert(n >= 1, 'n must be an integer larger than or equal to 1.');

    Mn.name = @() sprintf('[%s]^%d', M.name(), n);

    Mn.dim = @() n*M.dim();

    Mn.inner = @inner;
    function val = inner(x, u, v)
        val = 0;
        for i = 1 : n
            val = val + M.inner(x{i}, u{i}, v{i});
        end
    end

    Mn.norm = @(x, d) sqrt(Mn.inner(x, d, d));

    Mn.dist = @dist;
    function d = dist(x, y)
        sqd = 0;
        for i = 1 : n
            sqd = sqd + M.dist(x{i}, y{i})^2;
        end
        d = sqrt(sqd);
    end

    Mn.typicaldist = @typicaldist;
    function d = typicaldist()
        sqd = 0;
        for i = 1 : n
            sqd = sqd + M.typicaldist()^2;
        end
        d = sqrt(sqd);
    end

    Mn.proj = @proj;
    function u = proj(x, u)
        for i = 1 : n
            u{i} = M.proj(x{i}, u{i});
        end
    end

    Mn.tangent = @tangent;
    function u = tangent(x, u)
        for i = 1 : n
            u{i} = M.tangent(x{i}, u{i});
        end
    end

    if isfield(M, 'tangent2ambient')
        Mn.tangent2ambient = @tangent2ambient;
    else
        Mn.tangent2ambient = @(x, u) u;
    end
    function u = tangent2ambient(x, u)
        for i = 1 : n
            u{i} = M.tangent2ambient(x{i}, u{i});
        end
    end

    Mn.egrad2rgrad = @egrad2rgrad;
    function g = egrad2rgrad(x, g)
        for i = 1 : n
            g{i} = M.egrad2rgrad(x{i}, g{i});
        end
    end

    Mn.ehess2rhess = @ehess2rhess;
    function h = ehess2rhess(x, eg, eh, h)
        for i = 1 : n
            h{i} = M.ehess2rhess(x{i}, eg{i}, eh{i}, h{i});
        end
    end

    Mn.exp = @expo;
    function x = expo(x, u, t)
        if nargin < 3
            t = 1.0;
        end
        for i = 1 : n
            x{i} = M.exp(x{i}, u{i}, t);
        end
    end

    Mn.retr = @retr;
    function x = retr(x, u, t)
        if nargin < 3
            t = 1.0;
        end
        for i = 1 : n
            x{i} = M.retr(x{i}, u{i}, t);
        end
    end

    if isfield(M, 'log')
        Mn.log = @loga;
    end
    function u = loga(x, y)
        u = cell(n, 1);
        for i = 1 : n
            u{i} = M.log(x{i}, y{i});
        end
    end

    Mn.hash = @hash;
    function str = hash(x)
        str = '';
        for i = 1 : n
            str = [str M.hash(x{i})]; %#ok<AGROW>
        end
        str = ['z' hashmd5(str)];
    end

    Mn.lincomb = @lincomb;
    function x = lincomb(x, a1, u1, a2, u2)
        if nargin == 3
            for i = 1 : n
                x{i} = M.lincomb(x{i}, a1, u1{i});
            end
        elseif nargin == 5
            for i = 1 : n
                x{i} = M.lincomb(x{i}, a1, u1{i}, a2, u2{i});
            end
        else
            error('Bad usage of powermanifold.lincomb');
        end
    end

    Mn.rand = @rand;
    function x = rand()
        x = cell(n, 1);
        for i = 1 : n
            x{i} = M.rand();
        end
    end

    Mn.randvec = @randvec;
    function u = randvec(x)
        u = cell(n, 1);
        for i = 1 : n
            u{i} = M.randvec(x{i});
        end
        u = Mn.lincomb(x, 1/sqrt(n), u);
    end

    Mn.zerovec = @zerovec;
    function u = zerovec(x)
        u = cell(n, 1);
        for i = 1 : n
            u{i} = M.zerovec(x{i});
        end
    end

    if isfield(M, 'transp')
        Mn.transp = @transp;
    end
    function u = transp(x1, x2, u)
        for i = 1 : n
            u{i} = M.transp(x1{i}, x2{i}, u{i});
        end
    end

    if isfield(M, 'pairmean')
        Mn.pairmean = @pairmean;
    end
    function y = pairmean(x1, x2)
        y = cell(n, 1);
        for i = 1 : n
            y{i} = M.pairmean(x1{i}, x2{i});
        end
    end

    % Compute the length of a vectorized tangent vector of M at x, assuming
    % this length is independent of the point x (that should be fine).
    if isfield(M, 'vec')
        rand_x = M.rand();
        zero_u = M.zerovec(rand_x);
        len_vec = length(M.vec(rand_x, zero_u));

        Mn.vec = @vec;

        if isfield(M, 'mat')
            Mn.mat = @mat;
        end

    end

    function u_vec = vec(x, u_mat)
        u_vec = zeros(len_vec, n);
        for i = 1 : n
            u_vec(:, i) = M.vec(x{i}, u_mat{i});
        end
        u_vec = u_vec(:);
    end

    function u_mat = mat(x, u_vec)
        u_mat = cell(n, 1);
        u_vec = reshape(u_vec, len_vec, n);
        for i = 1 : n
            u_mat{i} = M.mat(x{i}, u_vec(:, i));
        end
    end

    if isfield(M, 'vecmatareisometries')
        Mn.vecmatareisometries = M.vecmatareisometries;
    else
        Mn.vecmatareisometries = @() false;
    end

end

author: andrewssobral

repo: lrslibrary


text:
function error_handler(ME, optional)
    % unhandled_error provides a one-stop place to track certain errors
    %
    % it is designed to be a drop-in replacement for simple try-catch display
    %
    % usage:
    %   try
    %      do something that might fail
    %   catch ME
    %      unhandled_error(ME) % prints the error message
    %      unhandled_error(ME, "oops") % prints "oops"
    %      unhandled_error(ME, @do_nothing) % does nothing
    %      unhandled_error(ME, @function_name) % executes function function_name
    %   end
    %
    %  Celso G Reyes, 2017

    if exist('optional','var') && ~isempty(optional)
        if isa('optional','function_handle')
            optional();
        else
            disp(optional)
        end
    elseif exist('ME', 'var')
        disp('ignored exception:');
        disp(ME)
    else
        % do nothing
    end

end

author: CelsoReyes


text:
function logPr = calcJointLogPr_BPHMMState( Psi, data )
%  OUTPUT
%     logPr : struct with fields
%                .F     : log p( F | hypers)
%                .obs   : log p( obs_ii | F_ii, z_ii )  sum over all ii
%                .z     : log p( z_ii | F_ii, hypers )  sum over all ii
%                .all   : log p( obs, F, z | hypers )

F = Psi.F > 0;
gamma = Psi.bpM.gamma;
c     = Psi.bpM.c;
TransM = Psi.TransM;
ThetaM = Psi.ThetaM;
stateSeq = Psi.stateSeq;

% -----------------------------   Compute prob. of binary feature mat F
logPr.F = calcLogPrFeatureMatrix( F, gamma, c );

logPr.z = TransM.calcMargPrStateSeq( F, stateSeq );

% Remember, no arguments means use stored suff stats
if ~exist( 'data', 'var' )
    logPr.obs = ThetaM.calcMargPrData( );
else
    logPr.obs = ThetaM.calcMargPrData( data, stateSeq );
end

% ============== combine all logPr into joint prob of chain state
logPrFields = fieldnames( logPr );
logPr.all = 0;
for ff = 1:length( logPrFields )
    logPr.all = logPr.all + logPr.( logPrFields{ff} );
end

author: michaelchughes

repo: NPBayesHMM


text:
% @Author: amishkin
% @Date:   2018-07-10T13:36:34-07:00
% @Email:  amishkin@cs.ubc.ca
% @Last modified by:   aaronmishkin
% @Last modified time: 2018-07-25T13:18:48-07:00

% #################################################
% ############# Reproduce Figure 2(c) #############
% #################################################

% Load the dataset with a fixed seed.
addpath(genpath('.'))

output_dir = './breast_cancer_experiment_data'
mkdir(output_dir)

dataset_name = 'breast_cancer_scale';
[y, X, y_te, X_te] = get_data_log_reg(dataset_name, 1);
[N,D] = size(X);

% Set initial conditions. These will be the same for all methods and restarts.
mu_start = zeros(D,1);
s = 1.*ones(D,1);
sigma_start = diag(1./s);

% Use the same split of the dataset for each random restart.
random_split = 0;

% Run MF_Exact.
method = 'mf_exact'
num_restarts = 1;

% List batch sizes and epochs for each batch size.
% In the case of mf-exact (which is a batch method), the number of epochs controls the number of function
% evaluations L-BFGS is allowed.
epoch_list = [500];
% M_list does not matter for mf-exact because it is a batch method.
M_list = [1]

run_experiment(dataset_name, method, M_list, epoch_list, 0, 0, 0, 0, num_restarts, mu_start, sigma_start, random_split, output_dir)

% Run Vadam with different minibatch sizes:
method = 'Vadam'

% Define the experiment parameters:
num_samples = 1
beta = 0.01;
alpha = 0.01;
decay_rate = 0.55   % controls the rate at which alpha and beta decay during training.

num_restarts = 20;

%  List batch sizes and epochs for each batch size.
M_list = [1,8,16,32,64];
epoch_list = [10000, 10000, 10000, 10000, 10000];

run_experiment(dataset_name, method, M_list, epoch_list, alpha, beta, decay_rate, num_samples, num_restarts, mu_start, sigma_start, random_split, output_dir);

% Run VOGN with a minibatch size of one:
method = 'VOGN'
num_samples = 1
beta = 0.0005;
alpha = 0.0005;
decay_rate = 0;

num_restarts = 20;

%  List batch sizes and epochs for each batch size.
M_list = [1];
epoch_list = [2000];

run_experiment(dataset_name, method, M_list, epoch_list, alpha, beta, decay_rate, num_samples, num_restarts, mu_start, sigma_start, random_split, output_dir);

% generate figure 2(c) using the data from this experiment:
make_fig_two_c

author: emtiyaz

repo: vadam

save_path: github-repos/MATLAB/emtiyaz-vadam


text:
function gnuplot_taylor_test ( )

%*****************************************************************************80
%
%% GNUPLOT_TAYLOR_TEST generates a field on a regular grid and plots it.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/navier_stokes_2d_exact/gnuplot_taylor_test.m
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GNUPLOT_TAYLOR_TEST:\n' );
  fprintf ( 1, '  Taylor Vortex Flow:\n' );
  fprintf ( 1, '  Generate a velocity field on a regular grid.\n' );
  fprintf ( 1, '  Store in GNUPLOT data and command files.\n' );

  x_lo = 0.5;
  x_hi = 2.5;
  x_num = 21;

  y_lo = 0.5;
  y_hi = 2.5;
  y_num = 21;

  [ x, y ] = grid_2d ( x_num, x_lo, x_hi, y_num, y_lo, y_hi );

  nu = 1.0;
  rho = 1.0;
  n = x_num * y_num;
  t = 0.0;

  [ u, v ] = uvp_taylor ( nu, rho, n, x, y, t );

  header = 'taylor';
  s = 0.10;
  ns2de_gnuplot ( header, n, x, y, u, v, s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function a = propa_no_random ( prob, k, n, key )

%*****************************************************************************80
%
%% PROPA_NO_RANDOM returns a random matrix that does not have property A.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real PROB, the probability that a link between
%    two eligible nodes will be made.
%
%    Input, integer K, the number of illegal links between nodes
%    to make.  The routine will TRY to make this many illegal links.  However,
%    it is obviously possible to make K too big.
%
%    Input, integer N, the order of A.
%
%    Input, integer KEY, a positive value that selects the data.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );
%
%  Assign each index randomly to one of two sets.
%  SET(I) is 0 if I is in set 0, and 1 if it is in set 1.
%
  seed = key;
  [ set, seed ] = sub_random ( n, seed );

  for i = 1 : n
    for j = 1 : n
      if ( set(i) ~= set(j) )
        [ chance, seed ] = r8_uniform_01 ( seed );
        if ( chance <= prob )
          a(i,j) = 1.0;
        end
      end
    end
  end
%
%  Now repeatedly pick a pair of indices, and consider setting the
%  corresponding entry of A to 1.
%
  bad = 0;
  tries = 0;

  while ( 1 )

    tries = tries + 1;

    if ( 1000 < tries )
      break
    end

    if ( k <= bad )
      break
    end

    [ i, seed ] = i4_uniform_ab ( 1, n, seed );
    [ j, seed ] = i4_uniform_ab ( 1, n, seed );

    if ( i == j )
      continue
    end

    if ( set(i) ~= set(j) )
      continue
    end

    if ( a(i,j) ~= 0.0 & a(j,i) ~= 0.0 )
      continue
    end

    if ( a(i,j) == 0.0 )
      a(i,j) = 1.0;
    else
      a(j,i) = 1.0;
    end

    bad = bad + 1;

  end

  return
end


text:
function [G,F0] = spm_voice_filter(Y,FS,F1,F2)
% Time frequency decomposition to characterise acoustic spectral envelope
% FORMAT [G,F0] = spm_voice_filter(Y,FS)
%
% Y    - timeseries
% FS   - sampling frequency
% F1   - lower frequency bound [default: 1024  Hz]
% F2   - upper frequency bound [default: 16096 Hz]
%
% G    - power at acoutic frequencies
% F0   - fundamental frequency
%
% This auxiliary routine uses a wavelet decomposition (complex Gaussian wavelets) to
% assess the power frequency range (F1 - F2 Hz). This can be used to
% identify the onset of a word or fast modulations of spectral energy at a
% fundamental frequency F0 of 256 Hz.
%
% This routine is not used for voice recognition but can be useful for
% diagnostics and plotting spectral envelope.
%
% see also: spm_voice_check.m
%__________________________________________________________________________
% Copyright (C) 2019 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_voice_filter.m 7750 2019-12-05 17:54:29Z spm $

% defaults
%--------------------------------------------------------------------------
if nargin < 3; F1 = 1024;  end
if nargin < 4; F2 = 16096; end

% find acoutic energy using spm_wft
%==========================================================================
F0    = 256;
k     = 1:round(F2/F0);                      % cycles per window       % Hz
k     = k(k > F1/F0 & k < F2/F0);            % acoustic range
n     = round(FS/(F0/2))*2;                  % window length (F0 Hz)
g     = abs(spm_wft(Y,k,n));                 % wavlet transform

% instantaneous power (250 - 5000 Hz)
%--------------------------------------------------------------------------
G     = sum(g)';

% find fundamental frequencies
%==========================================================================

% get maxium in the range of F0 (100 - 300Hz)
%--------------------------------------------------------------------------
fG    = abs(fft(G(:)));
nf    = length(fG);
w     = (1:nf)/(nf/FS);
i     = find(w > 64 & w < 300);
[~,j] = max(fG(i));
F0    = w(i(1) + j - 1);

% graphics
%==========================================================================

% time-frequency analysis
%--------------------------------------------------------------------------
subplot(3,1,1)
i    = 1:round(min(size(g,2),FS));
pst  = i*1000/FS;
imagesc(pst,k*F0,g(:,i))
title('time-frequency analysis','FontSize',16), xlabel('time (ms)')

% Energy
%--------------------------------------------------------------------------
subplot(3,1,2);
i    = 1:size(g,2);
pst  = i/FS;
plot(pst,G), title('Energy','FontSize',16), xlabel('time (s)')

% Fundamental frequency
%--------------------------------------------------------------------------
subplot(3,1,3)
i     = find(w > 64 & w < 300);
plot(w(i),fG(i)), hold on
plot([F0 F0],[0 max(fG(i))],'r'), hold off
title('Fundamental frequency','FontSize',16)
xlabel('time (seconds)')

author: spm


text:
function [sys,x0,str,ts]=ADRC_3(t,x,u,flag,h,TD,ESO,NLSEF,b0)

switch flag
    case 0
        [sys,x0,str,ts]=mdlInitializeSizes(h);
    case 2
        sys=mdlUpdate(x,u,h,TD,ESO,b0);
    case 3
        sys=mdlOutputs(x,NLSEF,b0);
    case 4,
        sys=mdlGetTimeOfNextVarHit(t,h);
    case {1,9}
        sys=[];
    otherwise
        error(['Unhandled flag=',num2str(flag)]);
end
function [sys,x0,str,ts]=mdlInitializeSizes(h)
    sizes=simsizes;
    sizes.NumContStates=0;
    sizes.NumDiscStates=7;
    sizes.NumOutputs=2;
    sizes.NumInputs=3;
    sizes.DirFeedthrough=1;
    sizes.NumSampleTimes=1;
    sys=simsizes(sizes);
    x0=[0;0;0;0;0;0;0];
    str=[];
    ts=[h 0];
function sys=mdlUpdate(x,u,h,TD,ESO,b0)
    e1=x(1)-u(1);
    fh=fhan(e1,x(2),TD(1),TD(2));
    sys(1)=x(1)+h*x(2);
    sys(2)=x(2)+h*x(3);
    sys(3)=x(3)+h*fh;
    e2=x(4)-u(2);
    sys(4)=x(4)+h*(x(5)-ESO(1)*e2);
    sys(5)=x(5)+h*(x(6)-ESO(2)*fal(e2,0.5,ESO(5)));
    sys(6)=x(6)+h*(x(7)-ESO(3)*fal(e2,0.25,ESO(5))+b0*u(3));
    sys(7)=x(7)+h*(-ESO(4)*fal(e2,0.125,ESO(5)));

%     fe=fal(e,0.5,Delta);
%     fe1=fal(e,0.25,Delta);
%     sys(1,1)=x(1)+h2*(x(2)-BB(1)*e);
%     sys(2,1)=x(2)+h2*(x(3)-BB(2)*fe+u(1));
%     sys(3,1)=x(3)+h2*(-BB(3)*fe1);
function sys=mdlOutputs(x,NLSEF,b0)
    e3=x(1)-x(4);
    e4=x(2)-x(5);
    e5=x(3)-x(6);
    %sys(1)=NLSEF(1)*e3+NLSEF(2)*e4+NLSEF(3)*e5-x(7)/b0;
    sys(1)=NLSEF(1)*fal(e3,0.5,NLSEF(4))+NLSEF(2)*fal(e4,0.75,NLSEF(4))+NLSEF(3)*fal(e5,1.5,NLSEF(4))-x(7)/b0;
    sys(2)=x(1);
function sys=mdlGetTimeOfNextVarHit(t,h)
    sys=t+h;

function y=fhan(x1,x2,r,h)
d=r*h;
d0=h*d;
y=x1+h*x2;
a0=sqrt(d^2+8*r*abs(y));
if abs(y)>d0
    a=x2+(a0-d)*sign(y)/2;
else
    a=x2+y/h;
end

if abs(a)>d
    y=-r*sign(a);
else
    y=-r*a/d;
end

function f=fal(e,a,d)
    if abs(e)<d
        f=e*d^(a-1);
    else f=(abs(e))^a*sign(e);
    end

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
x=1:20; y=1:10; z=-10:10;
[x,y,z]=meshgrid(x,y,z);
v=x.^2.*y.*(z+1);
isosurface(x,y,z,v)

author: Eurus-Holmes


text:
function K = covSM(Q, hyp, x, z, i)

% Gaussian Spectral Mixture covariance function. The covariance function
% parametrization depends on the sign of Q.
%
% Let t(Dx1) be an offset vector in dataspace e.g. t = x_i - z_j. Then w(DxP)
% are the weights and m(Dx|Q|) = 1/p, v(Dx|Q|) = (2*pi*ell)^-2 are spectral
% means (frequencies) and variances, where p is the period and ell the length
% scale of the Gabor function h(t2v,tm) given by the expression
%   h(t2v,tm) = exp(-2*pi^2*t2v).*cos(2*pi*tm)
%
% Then, the two covariances are obtained as follows:
%
% SM, spectral mixture:  Q>0 => P = 1
%   k(x_i,z_j) = w'*h(v'*(t.*t),m'*t)
%
% SMP, spectral mixture product: Q<0 => P = D
%   k(x_i,z_j) = prod(w'*h(T*T*v,T*m)), T = diag(t)
%
% The hyperparameters are:
%
% hyp = [ log(w(:))
%         log(m(:))
%         log(sqrt(v(:))) ]
%
% For more help on design of covariance functions, try "help covFunctions".
%
% Note that the spectral density H(s) = F[ h(t) ] of covGaboriso is given by
% H(s) = N(s|m,v)/2 + N(s|-m,v)/2 where m=1/p is the mean and v=(2*pi*ell)^-2
% is the variance of a symmetric Gaussian mixture. Hence the covGaboriso
% covariance forms a basis for the class of stationary covariances since a
% weighted sum of covGaboriso covariances corresponds to an isotropic
% location-scale mixture of a symmetric Gaussian mixture in the spectral domain.
%
% Internally, covSM constructs a weighted sum of products of 1d covGaboriso
% covariances using covMask, covProd, covScale and covSum.
%
% For more details, see
% 1) Gaussian Process Kernels for Pattern Discovery and Extrapolation,
% ICML, 2013, by Andrew Gordon Wilson and Ryan Prescott Adams.
% 2) GPatt: Fast Multidimensional Pattern Extrapolation with Gaussian
% Processes, arXiv 1310.5288, 2013, by Andrew Gordon Wilson, Elad Gilboa,
% Arye Nehorai and John P. Cunningham, and
% http://mlg.eng.cam.ac.uk/andrew/pattern
%
% For Q>0, covSM corresponds to Eq. 12 in Ref (1)
% For Q<0, covSM corresponds to Eq. 14 in Ref (2) (but w here = w^2 in (14))
%
% Copyright (c) by Andrew Gordon Wilson and Hannes Nickisch, 2014-09-24.
%
% See also COVFUNCTIONS.M, COVGABORISO.M, COVGABORARD.M.

smp = Q<0; Q = abs(Q);                    % switch between covSM and covSMP mode
if nargin<3                                            % report no of parameters
  if smp, K = '3*D*'; else K = '(1+2*D)*'; end, K = [K,sprintf('%d',Q)]; return
end
if nargin<4, z = []; end                                   % make sure, z exists

D = size(x,2); P = smp*D+(1-smp);                   % dimensionality, P=D or P=1
lw = reshape(hyp(         1:P*Q) ,P,Q);                    % log mixture weights
lm = reshape(hyp(P*Q+    (1:D*Q)),D,Q);                     % log spectral means
ls = reshape(hyp(P*Q+D*Q+(1:D*Q)),D,Q);       % log spectral standard deviations

% In the following, we construct nested cell arrays to finally obtain either
if smp % 1) the product of weighted sums of 1d covGabor covariance functions or
  fac = cell(1,D);
  for d=1:D
    add = cell(1,Q); % a) addends for weighted sum of univariate Gabor functions
    for q=1:Q, add{q} = {'covScale',{'covMask',{d,{'covGaboriso'}}}}; end
    fac{d} = {'covSum',add};                       % b) combine addends into sum
  end
  cov = {'covProd',fac};                       % c) combine factors into product
else   % 2) the weighted sum of multivariate covGaborard covariance functions.
                                  % weighted sum of multivariate Gabor functions
  add = cell(1,Q); for q=1:Q, add{q} = {'covScale',{'covGaborard'}};  end
  cov = {'covSum',add};                                       % combine into sum
end
if smp      % assemble hyp; covGabor is parametrised using -ls-log(2*pi) and -lm
  hyp = [lw(:)'/2; -ls(:)'-log(2*pi); -lm(:)'];
else
  hyp = [lw/2;     -ls-log(2*pi);     -lm    ];
end

if nargin<5                                       % evaluation of the covariance
  K = feval(cov{:},hyp(:),x,z);
else
  % We compute the indices j in the new hyperparameter vector hyp. The
  % correction constants c are needed because some hyperparameters in hyp are
  % powers of hyperparameters accepted by covSM.
  if i<=P*Q                                               % derivatives w.r.t. w
    c =  0.5;
    if smp, j = 1+3*(i-1); else j = 1+(i-1)*(2*D+1); end
  elseif i<=(P+  D)*Q                                     % derivatives w.r.t. m
    c = -1.0; j = i- P   *Q; [j1,j2] = ind2sub([D,Q],j);
    if smp, j = 3+3*(j-1); else j = 1+j1+D+(j2-1)*(2*D+1); end
  elseif i<=(P+2*D)*Q                                     % derivatives w.r.t. v
    c = -1.0; j = i-(P+D)*Q; [j1,j2] = ind2sub([D,Q],j);
    if smp, j = 2+3*(j-1); else j = 1+j1+  (j2-1)*(2*D+1); end
  else
    error('Unknown hyperparameter')
  end
  K = c*feval(cov{:},hyp(:),x,z,j);
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function Cbar = computeCbar(A, B, ny, nu)
    Cbar_tmp = cell(ny,nu);
    tmp = B;

    for i=1:nu
       for k=1:(i-1)
           Cbar_tmp{k,i} = zeros(size(B));
       end

       for j=i:ny
           Cbar_tmp{j,i} = tmp;
           tmp = A*tmp;
       end

       tmp = B;
    end

    Cbar = cell2mat(Cbar_tmp);
end

author: ccalas

repo: mpc

save_path: github-repos/MATLAB/ccalas-mpc


text:
function obj=stat_ellipse(obj,varargin)
%stat_ellipse() Create confidence ellipses around 2D groups of
% points
%
% Parameters:
% 'type': The default '95percentile' displays an ellipse that
% contains 95% of the points (assuming a bivariate normal
% distribution). The option 'ci' will first compute boostrapped
% 2D means and plot the 95% ellipse around these means
% 'geom': Sets how to display the result 'area' for a shaded
% area or 'line' for a simple contour line.
% 'patch_opts': Provide additional patch properties as name-value
% pairs in a cell array (as if those were options for Matlab's
% built in patch() function)

p=inputParser;
my_addParameter(p,'type','95percentile'); %ci
my_addParameter(p,'geom','area'); %line
my_addParameter(p,'patch_opts',{});
parse(p,varargin{:});

obj.geom=vertcat(obj.geom,{@(dobj,dd)my_ellipse(dobj,dd,p.Results)});
obj.results.stat_ellipse={};
end

function hndl=my_ellipse(obj,draw_data,params)

persistent elpoints;
persistent sphpoints;

%Cache unity ellipse points
if isempty(elpoints)
    res=30;
    ang=0:pi/(0.5*res):2*pi;
    elpoints=[cos(ang); sin(ang)];
    [x,y,z]=sphere(10);
    sphpoints=surf2patch(x,y,z);
end

combx=shiftdim(comb(draw_data.x));
comby=shiftdim(comb(draw_data.y));
combz=shiftdim(comb(draw_data.z));

%If we have "enough" points
if sum(~isnan(combx))>2 && sum(~isnan(comby))>2

    if isempty(draw_data.z)

        r=[combx comby];
        %Using a chi square with 2 degrees of freedom is proper
        %here (tested: generated ellipse do contain 1-alpha of the
        %points)
        k=@(alpha) sqrt(chi2inv(1-alpha,2));

        %If a CI on the mean is requested, we replace the original points
        %with bootstrapped mean samples
        if strcmp(params.type,'ci')
            r=bootstrp(obj.stat_options.nboot,@nanmean,r);
        end

        %Extract mean and covariance
        m=nanmean(r);
        cv=nancov(r);

        %Compute ellipse points
        conf_elpoints=sqrtm(cv)*elpoints*k(obj.stat_options.alpha);

        %Compute ellipse axes
        [evec,eval]=eig(cv);
        if eval(2,2)>eval(1,1) %Reorder
            evec=fliplr(evec);
            eval=fliplr(flipud(eval));
        end
        elaxes=sqrtm(cv)*evec*k(obj.stat_options.alpha);

        obj.results.stat_ellipse{obj.result_ind,1}.mean=m;
        obj.results.stat_ellipse{obj.result_ind,1}.cv=cv;
        obj.results.stat_ellipse{obj.result_ind,1}.major_axis=elaxes(:,1)';
        obj.results.stat_ellipse{obj.result_ind,1}.minor_axis=elaxes(:,2)';

        %plot([0 elaxes(1,1)]+m(1),[0 elaxes(2,1)]+m(2),'k')
        %plot([0 elaxes(1,2)]+m(1),[0 elaxes(2,2)]+m(2),'k')

        switch params.geom
            case 'area'
                hndl=patch(conf_elpoints(1,:)+m(1),conf_elpoints(2,:)+m(2),draw_data.color,'FaceColor',draw_data.color,'EdgeColor',draw_data.color,'LineWidth',2,'FaceAlpha',0.2);

            case 'line'
                hndl=patch(conf_elpoints(1,:)+m(1),conf_elpoints(2,:)+m(2),draw_data.color,'FaceColor','none','EdgeColor',draw_data.color,'LineWidth',2);
        end
        set(hndl,params.patch_opts{:});
        %One matlab version displayed stuff if no output value was set (but
        %crashes 2014a and earlier versions)
        %tmp = set(hndl,params.patch_opts{:});

        center_hndl=plot(m(1),m(2),'+','MarkerFaceColor',draw_data.color,'MarkerEdgeColor',draw_data.color,'MarkerSize',10);
    else

        r=[combx comby combz];
        k=@(alpha) sqrt(chi2inv(1-alpha,3));

        %If a CI on the mean is requested, we replace the original points
        %with bootstrapped mean samples
        if strcmp(params.type,'ci')
            r=bootstrp(obj.stat_options.nboot,@nanmean,r);
        end

        %Extract mean and covariance
        m=nanmean(r);
        cv=nancov(r);

        obj.results.stat_ellipse{obj.result_ind,1}.mean=m;
        obj.results.stat_ellipse{obj.result_ind,1}.cv=cv;
        obj.results.stat_ellipse{obj.result_ind,1}.major_axis=[];
        obj.results.stat_ellipse{obj.result_ind,1}.minor_axis=[];

        conf_sphpoints=sphpoints;
        conf_sphpoints.vertices=bsxfun(@plus,sqrtm(cv)*conf_sphpoints.vertices'*k(obj.stat_options.alpha),m')';
        hndl=patch(conf_sphpoints,'FaceColor',draw_data.color,'EdgeColor','none','LineWidth',2,'FaceAlpha',0.2);

        center_hndl=plot3(m(1),m(2),m(3),'+','MarkerFaceColor',draw_data.color,'MarkerEdgeColor',draw_data.color,'MarkerSize',10);
    end
    obj.results.stat_ellipse{obj.result_ind,1}.ellipse_handle=hndl;
    obj.results.stat_ellipse{obj.result_ind,1}.center_handle=center_hndl;
else
    warning('Not enough points for ellipse')

    obj.results.stat_ellipse{obj.result_ind,1}.mean=NaN;
    obj.results.stat_ellipse{obj.result_ind,1}.cv=NaN;
    obj.results.stat_ellipse{obj.result_ind,1}.major_axis=[];
    obj.results.stat_ellipse{obj.result_ind,1}.minor_axis=[];
    obj.results.stat_ellipse{obj.result_ind,1}.ellipse_handle=[];
    obj.results.stat_ellipse{obj.result_ind,1}.center_handle=[];
end
end

author: piermorel

repo: gramm

save_path: github-repos/MATLAB/piermorel-gramm


text:
function ir = compensate_headphone(ir,conf)
%COMPENSATE_HEADPHONE applies a headphone compensation to an impulse response
%
%   Usage: ir = compensate_headphone(ir,conf)
%
%   Input parameters:
%       ir      - Impulse response to which the compensation should be applied
%       conf    - configuration struct (see SFS_config)
%
%   Output:
%       ir      - Impulse response which is compensated for the given headphone
%
%   COMPENSATE_HEADPHONE(ir,conf) applies a headphone compensation to the
%   given impulse response. Which headphone compensation it should use is
%   mentioned in the conf struct. The compensation filter can be a one-channel
%   (same filter for left and right) or two-channel signal (1st signal: left,
%   2nd signal: right) with signals stored as columns in a matrix.
%   The compensation is only applied, if the conf.ir.usehcomp value is not false.
%
%   See also: ir_wfs, ir_point_source, ir_generic

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 2;
nargmax = 2;
narginchk(nargmin,nargmax);
isargmatrix(ir);
isargstruct(conf);

%% ===== Configuration ===================================================
usehcomp = conf.ir.usehcomp;

%% ===== Computation =====================================================
if(usehcomp)
    lenir = size(ir,1);
    % Read headphone compensation filter
    hcomp = audioread(conf.ir.hcompfile);
    % Check if the IR has the right length for the filter
    if lenir<length(hcomp)
        warning(['The length of the used IR is shorter than the headphone ', ...
            'compensation filter.']);
    end
    % Apply filter
    ir = convolution(hcomp,ir);
    ir = fix_length(ir,lenir);
end

author: sfstoolbox

repo: sfs-matlab


text:
function [ySqrtEst, PInvSqrtEst]=sqrtInfoBatchSmoother(ySqrtPred,PInvSqrtPred,z,u,H,F,SR,SQ,Gamma,kD)
%%SQRTINFOBATCHSMOOTHER Run the square root information smoother for linear
%                  dynamic and measurement models on a batch of
%                  measurements. The smoothed result at one time step or
%                  along the entire batch are available. The initial
%                  predicted states can not be uninformative.
%
%INPUTS: ySqrtPred The xDimX1 predicted square root information state at
%                 the time of the initial measurement in z. The predicted
%                 information state is always PInvSqrtPred times the
%                 predicted target state estimate.
%    PInvSqrtPred The inverse square root information matrix associated
%                 with the predicted information state at the time of the
%                 initial measurement in z. If P is the covariance matrix
%                 of a Gaussian state x, then P=PSqrt*PSqrt' and
%                 PInvSqrtPred=inv(PSqrt). This can be either upper
%                 triangular or lower triangular.
%               z The zDim X N matrix of measurements for the whole batch.
%               u The xDim X(N-1) matrix of control inputs for the whole
%                 batch. If there are no control inputs, then set u=[];
%               H The zDim X xDim X N hypermatrix of measurement matrices
%                 such that H(:,:,k)*x+w is the measurement at time k,
%                 where x is the state and w is zero-mean Gaussian noise
%                 with covariance matrix R (:,:,k). Alternatively, if all
%                 of the measurement matrices are the same, one can just
%                 pass a single zDim X xDim matrix.
%               F The xDim X xDim X (N-1) hypermatrix of state transition
%                 matrices. The state at discrete-time k+1 is modeled as
%                 F(:,:,k) times the state at time k plus zero-mean
%                 Gaussian process noise with covariance matrix Q(:,:,k).
%                 Alternatively, if all of the state transition matrices
%                 are the same, one can just pass a single xDim X xDim
%                 matrix.
%              SR The zDim X zDim X N hypermatrix of lower-triangular
%                 square-root of the measurement covariance matrices. The
%                 matrices must be invertible. Alternatively, if all of the
%                 measurement covariance matrices are the same, one can
%                 just pass a single zDim X zDim matrix.
%              SQ The xDim X xDim X (N-1) hypermatrix of lower-triangular
%                 square-root of the process noise covariance matrices.
%                 The matrices must be invertible. Alternatively, if all of
%                 the measurement covariance matrices are the same, one can
%                 just pass a single xDim X xDim matrix.
%           Gamma An optional xDim X xDim X (N-1) hypermatrix of matrices
%                 that transform the process noise to the state domain if
%                 the process noise covariance matrix is singular.
%                 Alternatively, if all of the transform matrices are the
%                 same, one can just pass a single xDim X xDim matrix. If
%                 this is omitted entirely, an identity matrix is used
%                 (i.e. there is no Gamma).
%              kD The discrete time-step at which the smoothed state
%                 estimate is desired, where z(:,1) is at discrete
%                 time-step 1 (not 0). If kD is omitted ot an empty matrix
%                 is passed, then results along the entire batch are
%                 obtained.
%
%OUTPUTS: ySqrtEst The xDimXN smoothed square root information state
%                  estimates at all steps if kD is not provided or the
%                  xDimX1 smoothed information state estimate at step kD if
%                  kD is provided.
%      PSqrtInvEst The inverse square root covariance matrices associated
%                  with the smoothed information state estimates. This is
%                  xDimXxDimXN for the whole batch if kD is not provided
%                  and is xDimXxDim if kD is provided.
%
%The algorithm is that of the linear square root information filter and
%smoother that is described in the book [1]. The forward pass steps can be
%found in Chapter IV, and the backwards smoother in Chapter X.
%
%The algorithm works by running a square root information filter forward
%in time and then smoothing backwards. The smoothing step depends on
%information stored during the forward pass. The smoothing step here is
%performed on a state residual rather than the state itself and applies the
%smoothed output to the previous updated measurement. This was based off of
%NASA's Orbit Determination Toolbox.
%
%z*(j+1)=PInvSqrt*(j+1)(x*(j+1)-xPred(j+1)
%A=[Rw(j)+Rwx(j),   Rwx(j)F(j+1,j),         u(j)
%   PInvSqrt*(j+1), PInvSqrt*(j+1)F(j+1,j), z*(j+1)]
%qr(A)=[Rw*(j),   Rwx*(j),      u*(j)
%       0,        PInvSqrt*(j), z*(j)]
%x*(j)=xUpd(j)+PInvSqrt*(j)\z*(j)
%
%where xPred is the forward predicted state, xUpd is the forward updated
%state, Rw and Rwx are outputs from the forward prediction step, u is the
%the control input (often zero), and x* and PInvSqrt* are the smoothed
%states. x*(N)=xUpd(N) and PInvSqrt*(N)=PInvSqrtUpd(N)
%
%REFERENCES:
%[1] G. J. Bierman, "Factorization Methods for Discrete Sequential
%    Estimation. Academic Press, New York, 1977.
%
%March 2015, David Karnick, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

xDim=size(H,2);
N=size(z,2);

if(nargin<9||isempty(Gamma))
    Gamma=eye(xDim);
end
if(nargin<10||isempty(kD))
    kD=[];
end
if(isempty(ySqrtPred))
    ySqrtPred=zeros(xDim,1);
end
if(isempty(PInvSqrtPred))
    PInvSqrtPred=zeros(xDim,xDim);
end
if(isempty(u))
    u=zeros(xDim,N-1);
end
if(size(H,3)==1)
    H=repmat(H,[1,1,N]);
end
if(size(F,3)==1)
    F=repmat(F,[1,1,N-1]);
end
if(size(SR,3)==1)
    SR=repmat(SR,[1,1,N]);
end
if(size(SQ,3)==1)
    SQ=repmat(SQ,[1,1,N-1]);
end
if(size(Gamma,3)==1)
    Gamma=repmat(Gamma,[1,1,N-1]);
end

%Run the SRIF forward until we have the predictions of step kD|kD-1 (going forwards).
ySqrtFwdPred=zeros(xDim,N);
PInvSqrtFwdPred=zeros(xDim,xDim,N);
ySqrtFwdUpd=zeros(xDim,N);
PInvSqrtFwdUpd=zeros(xDim,xDim,N);
Rw=zeros(xDim,xDim,N);
Rwx=zeros(xDim,xDim,N);

%The first step uses the priors
ySqrtFwdPred(:,1)=ySqrtPred;
PInvSqrtFwdPred(:,:,1)=PInvSqrtPred;

%The rest of the steps
for curStep=1:(N-1)
    [ySqrtFwdUpd(:,curStep),PInvSqrtFwdUpd(:,:,curStep)]=sqrtInfoFilterUpdate(ySqrtFwdPred(:,curStep),PInvSqrtFwdPred(:,:,curStep),z(:,curStep),SR(:,:,curStep),H(:,:,curStep));
    [ySqrtFwdPred(:,curStep+1),PInvSqrtFwdPred(:,:,curStep+1),Rw(:,:,curStep+1),Rwx(:,:,curStep+1)]=sqrtInfoFilterDiscPred(ySqrtFwdUpd(:,curStep),PInvSqrtFwdUpd(:,:,curStep),F(:,:,curStep),SQ(:,:,curStep),u(:,curStep),Gamma(:,:,curStep));
end

%Run the SRIS backwards until we have the updated information state of
%step kD|kD (going backwards).
ySqrtEst=zeros(xDim,N);
PInvSqrtEst=zeros(xDim,xDim,N);

%ODTBX method
[ySqrtEst(:,end),PInvSqrtEst(:,:,end)]=sqrtInfoFilterUpdate(ySqrtFwdPred(:,N),PInvSqrtFwdPred(:,:,N),z(:,N),SR(:,:,N),H(:,:,N));
for curStep=(N-1):-1:1
    xPred=PInvSqrtFwdPred(:,:,curStep+1)\ySqrtFwdPred(:,curStep+1); %Forward predicted state_j+1
    xStar=PInvSqrtEst(:,:,curStep+1)\ySqrtEst(:,curStep+1); %Backward smoothed state_j+1
    zStar=PInvSqrtEst(:,:,curStep+1)*(xStar-xPred);

    [zSmooth,PInvSqrtEst(:,:,curStep)]=sqrtInfoSmoothStep(zStar,PInvSqrtEst(:,:,curStep+1),F(:,:,curStep),Rw(:,:,curStep),Rwx(:,:,curStep),u(:,curStep),Gamma(:,:,curStep));

    xUpd=PInvSqrtFwdUpd(:,:,curStep)\ySqrtFwdUpd(:,curStep); %Forward updated state_j
    ySqrtEst(:,curStep)=PInvSqrtEst(:,:,curStep)*(xUpd+PInvSqrtEst(:,:,curStep)\zSmooth);
end

if(~isempty(kD))
    ySqrtEst=ySqrtEst(:,kD);
    PInvSqrtEst=PInvSqrtEst(:,:,kD);
end
end

function [ySqrtPred, PInvSqrtPred]=sqrtInfoSmoothStep(ySqrtPrev,PInvSqrtPrev,F,Ru,Rux,u,Gamma)
%This is the basic Smoothing step shown in Bierman's book
xDim=size(ySqrtPrev,1);

A=[Ru+Rux*Gamma, Rux*F, u;
    PInvSqrtPrev*Gamma,   PInvSqrtPrev*F,   ySqrtPrev];
[~,T]=qr(A);

ySqrtPred=T((xDim+1):end,end);
PInvSqrtPred=T((xDim+1):end,(end-xDim):(end-1));
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:
% User-Specified Parameters
% ===============================================================

% Conditions and stimuli
% ---------------------------------------------------------------
GA.conditions = [1];
GA.freqConditions = [.5];
GA.scanLength = 480;
GA.ISI = 2;
GA.TR = 2;

% Genetic algorithm parameters
% ---------------------------------------------------------------
nmodels = 1;
GA.cbalColinPowerWeights = [1 1 1 0];	% 1 = cbal, 2 = eff, 3 = hrf shape, 4 = freq
GA.numGenerations = 2;
GA.sizeGenerations = 20;
GA.maxTime = 600;						% max time to run in s, or Inf
GA.alph = 2.1;
GA.plotFlag = 0;

% Filtering, counterbalancing, and design tolerance
% ---------------------------------------------------------------
GA.lowerLimit = [];
GA.HPlength = [];
GA.LPsmooth = [];
GA.maxOrder = 1;
GA.NumStimthresh = [];
GA.maxCbalDevthresh = []; 	% counterbalancing deviation
GA.maxFreqDevthresh = [];

% Contrast setup
% ---------------------------------------------------------------
GA.contrasts = [1];
GA.contrastweights = [1];	% or predictor weights, if no contrasts

% Autocorrelation and special options
% ---------------------------------------------------------------
AutocorrelationFileName = 'myscannerxc';
GA.restlength = [];
GA.restevery = [];
GA.trans2switch = 0;
GA.trans2block = 0;
GA.dofirst = 0;
GA.nonlinthreshold = [];

eval(['load ' AutocorrelationFileName]);
GA.xc = myscannerxc;

% =============================================
% * vary by parameter - set this in this script
% * if running multiple models, nmodels > 1
% ---------------------------------------------

varyparam = 0;
fieldname = 'alph';		% or 'freqConditions(end), etc.
incrementby = 0.3;
incrementevery = 5;		% every n models
contrastsofinterest = 1;% contrasts or predictors, if no cons specified

% =============================================

if varyparam
	eval(['paramvalue = GA.' fieldname ';'])
	disp(' ');disp('* *********************************************************************')
	disp('*  You have selected to vary a parameter across models in ga_gui_script')
	disp(['*  GA.' fieldname ' starting at ' num2str(paramvalue)])
	disp('* *********************************************************************')
end

for nm = 1:nmodels

   eval(['diary model' num2str(nm) '.log'])
   disp(['Starting Model ' num2str(nm)])
    % freqConditions

   M = optimize_rand_search(GA);

   Models(:,nm) = M.stimlist;
   MM{nm} = M;
   save GAworkspace

   if varyparam
   	if isfield(M,'consebeta'),
		Fit(1,nm) = mean(M.consebeta(contrastsofinterest));
	else
		Fit(1,nm) = mean(M.sebeta(contrastsofinterest));
	end
   	eval(['Fit(2,nm) = GA.' fieldname ';'])
   	FC(nm,:) = GA.freqConditions;
   	ALPH(nm,:) = GA.alph;
	str = ['if mod(nm,' num2str(incrementevery) ') == 0,GA.' fieldname ' = GA.' fieldname '+' num2str(incrementby) ';,end']
   	eval(str);
	eval(['paramvalue = GA.' fieldname ';'])
	disp('* *********************************************************************')
	disp('*  Varying parameter option on.')
	disp(['*  GA.' fieldname ' is now ' num2str(paramvalue)])
	disp('* *********************************************************************')
   end

   eval(['save model' num2str(nm) '_' M.date ' M'])
   save GAworkspace
   diary off
   fprintf('\n')
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function modelConstrained = constrainRxnListAboveBound(model, rxnList, c, d, ineqSense)
% Constrains one (weighted) sum of fluxes to be above a lower bound.
% Appends to existing inequality constraints if they are present
%
% USAGE:
%
%    modelConstrained = constrainRxnListAboveBound(model, rxnList, c, d, ineqSense)
%
% INPUTS:
%    model:               model structure
%    rxnList:             cell array of reaction names
%
% OPTIONAL INPUTS:
%    c:                   `k x 1` vector :math:`c*v \geq d`
%    d:                   `n x 1` vector :math:`c*v \geq d`
%    ineqSense:           `k x 1` inequality sense {'L','G'}
%
% OUTPUT:
%    modelConstrained:    constrained model:
%
%                           * S - Stoichiometric matrix
%                           * b - Right hand side = dx/dt
%                           * C - Inequality constraint matrix
%                           * d - Inequality constraint right hand side
%                             :math:`[S; C] * v {=, \leq, \geq } [dxdt, d]`

% EXAMPLE:
%
%    rxnList = {'PCHOLP_hs_f', 'PLA2_2_f', 'SMS_f','PCHOLP_hs_b', 'PLA2_2_b', 'SMS_b'};
%    c = [1, 1, 1, 1, 1, 1];
%    d = 10;
%    ineqSense = 'G';
%    modelConstrained = constrainRxnListAboveBound(modelIrrev, rxnList, C, d, ineqSense);

modelConstrained = addCOBRAConstraints(model, rxnList, d, 'c', c, 'dsense', ineqSense);

author: opencobra

repo: cobratoolbox


text:
function lineStr = getline(FID, comment)

% GETLINE Get a line from a file.
% FORMAT
% DESC gets a line from a file, but ignore it if it starts with a
% comment character.
% ARG fid : the identity of the file to load from.
% ARG comment : the character that indicates a line is a comment
% (default #).
% RETURN lineStr : the next line in the file.
%
% COPYRIGHT : Neil D. Lawrence, 2005
%
% SEEALSO : fgetl, fopen

% NDLUTIL

if nargin < 2
  comment = '#';
end
if length(comment)~=1
  error('Comment can only be one character in length');
end
lineStr = fgetl(FID);
if length(lineStr)==0
  return
end
while lineStr(1)==comment
  lineStr = fgetl(FID);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%% housekeeping
clear
close all
clc
%% RISE the model and assign the steady state file
linear=~true;
if linear
    m=rise('emosp_linear','rise_flags',{'original',false},...
        'solve_linear',true);
    % N.B: original =false removes the variables that are both
    % predetermined and forward-looking, making it easier for Peter Ireland
    % to set up the matrices of his system.
else
    m=rise('emosp','steady_state_file','sstate_file',...
        'steady_state_imposed',true);
end
%% Select the model

sticky_price_model=true;
%% get the parameters
[p,priors]=create_parameters(sticky_price_model);
%% push the baseline calibration
m=set(m,'parameters',p);
%% do various things (solving, irfs, vardec, simulation, etc...)
clc
[ms,retcode]=solve(m);
ms.print_solution
%% Select the sample to use for estimation/filtering
choice=2;
samples={
    [],'1979Q2'
    '1979Q3',[]
    [],[]
    };
start_date=samples{choice,1};
end_date=samples{choice,2};
%% load the data
data=create_data(start_date,end_date);

%% estimate the model
log_vars={};
if ~linear
    % The nonlinear model creates problem for the computation of the
    % covariance matrix. Here we take a log-expansion of the state
    % variables
    log_vars=get(m,'endo_list(state)'); % log_vars={'A','E','K','M','V','X','Z'};
end
clc
profile off
profile on
ms=estimate(m,'data',data,...
    'estim_priors',priors,...
    'estim_start_date',start_date,...
    'solve_log_approx_vars',log_vars);%,'kf_tol',0
profile off
profile viewer
%% Redo filtering?
% clc
% [mfilt,LogLik,Incr,retcode]=filter(m,'data',data,'estim_start_date',start_date);%,'kf_tol',0
%% Take a log-linear approximation from a linear approximation
ml=rise('emosp_linear','rise_flags',{'original',true});
mnl=rise('emosp','steady_state_file','sstate_file');
%% Solve and print the solution
clc
M=set([ml,...
    set(mnl,'solve_log_approx_vars',['C',get(mnl,'endo_list(state)')])...
    ],...
    'parameters',p);

M=solve(M);
M.print_solution(['C',M(1).observables.name])%

author: jmaih


text:
function result = triangle_sub ( func, xval, yval, nsub, norder, xtab, ytab, ...
  weight )

%*****************************************************************************80
%
%% TRIANGLE_SUB carries out quadrature over subdivisions of a triangular region.
%
%  Integration region:
%
%    Points (X,Y) such that:
%
%      (X,Y) =
%          ALPHA                * ( XVAL(1), YVAL(1) )
%        + BETA                 * ( XVAL(2), YVAL(2) )
%        + ( 1 - ALPHA - BETA ) * ( XVAL(3), YVAL(3) )
%      0 <= ALPHA <= 1 - BETA
%      0 <= BETA <= 1 - ALPHA
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied function of
%    two variables which is to be integrated, of the form:
%      function value = func ( x, y )
%
%    Input, real XVAL(3), YVAL(3), the coordinates of the triangle vertices.
%
%    Input, integer NSUB, the number of subdivisions of each side of the
%    input triangle to be made.  NSUB = 1 means no subdivisions are made.
%    NSUB = 3 means that each side of the triangle is subdivided into
%    three portions, and that the original triangle is subdivided into
%    NSUB**2 triangles.  NSUB must be at least 1.
%
%    Input, integer NORDER, the order of the rule.
%
%    Input, real XTAB(NORDER), YTAB(NORDER), the abscissas.
%
%    Input, real WEIGHT(NORDER), the weights of the rule.
%
%    Output, real RESULT, the approximate integral of the function.
%

%
%  Initialize RESULT, the approximate integral.
%
  result = 0.0E+00;
%
%  NSUB must be positive.
%
  if ( nsub <= 0 )
    return;
  end
%
%  Initialize QUAD, the quadrature sum.
%
  quad = 0.0E+00;
%
%  The sub-triangles can be grouped into NSUB strips.
%
  for i = 1 : nsub

    temp1 = 0.0;
    temp2 = i / nsub;

    x2 = xval(2) + temp1 * ( xval(3) - xval(2) ) ...
                 + temp2 * ( xval(1) - xval(2) );

    y2 = yval(2) + temp1 * ( yval(3) - yval(2) ) ...
                 + temp2 * ( yval(1) - yval(2) );

    temp1 = 0.0;
    temp2 = ( i - 1 ) / nsub;

    x3 = xval(2) + temp1 * ( xval(3) - xval(2) ) ...
                 + temp2 * ( xval(1) - xval(2) );

    y3 = yval(2) + temp1 * ( yval(3) - yval(2) ) ...
                 + temp2 * ( yval(1) - yval(2) );
%
%  There are 2*I-1 triangles in strip number I.
%  The next triangle in the strip shares two nodes with the previous one.
%  Compute its corners, (X1,Y1), (X2,Y2), (X3,Y3).
%
    for j = 1 : 2*i-1

      x1 = x2;
      y1 = y2;
      x2 = x3;
      y2 = y3;
      temp1 = floor ( ( j + 1 ) / 2 ) / nsub;
      temp2 = floor ( i - 1 - ( j / 2 ) ) / nsub;

      x3 = xval(2) + temp1 * ( xval(3) - xval(2) ) ...
                   + temp2 * ( xval(1) - xval(2) );

      y3 = yval(2) + temp1 * ( yval(3) - yval(2) ) ...
                   + temp2 * ( yval(1) - yval(2) );
%
%  Now integrate over the triangle, mapping the points ( XTAB(K), YTAB(K) )
%  into the triangle.
%
      for k = 1 : norder

        x = x2 + xtab(k) * ( x3 - x2 ) + ytab(k) * ( x1 - x2 );
        y = y2 + xtab(k) * ( y3 - y2 ) + ytab(k) * ( y1 - y2 );
        quad = quad + weight(k) * feval ( func, x, y );

       end

    end

  end

  volume = triangle_volume ( xval, yval ) / ( nsub * nsub );
  result = quad * volume;

  return
end


text:
function value = p16_exact ( dim_num )

%*****************************************************************************80
%
%% P16_EXACT returns the exact integral for problem 16.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Output, real VALUE, the exact value of the integral.
%

%
%  Get the limits of integration.
%
  [ a, b ] = p16_lim ( dim_num );
%
%  Get the location of Z.
%
  z = [];
  z = p16_r8vec ( 'G', 'Z', dim_num, z );
%
%  The value of the DIM_NUM dimensional integral can be broken down
%  into the weighted sum of 1 dimensional integrals.
%
  value = 0.0;

  volume = prod ( b(1:dim_num) - a(1:dim_num) );

  for i = 1 : dim_num
%
%  Z < A < B
%
    if ( z(i) < a(i) )

      integral = 0.5 * ( b(i) - a(i) ) * ( b(i) + a(i) - 2.0 * z(i) );
%
%  A < Z < B
%
    elseif ( z(i) < b(i) )

      integral = 0.5 * ( ( b(i) - z(i) )^2 + ( z(i) - a(i) )^2 );
%
%  A < B < Z
%
    else

      integral = 0.5 * ( b(i) - a(i) ) * ( 2.0 * z(i) - a(i) - b(i) );

    end

    value = value + volume * integral / ( b(i) - a(i) );

  end

  return
end


text:
function write_neuralynx_nts(filename, nts)

% WRITE_NEURALYNX_NTS writes spike timestamps to a NTS file
%
% Use as
%   write_neuralynx_nts(filename, nts)
%
% See also READ_NEURALYNX_NTS

% Copyright (C) 2007, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if ~isa(nts.TimeStamp, 'uint64')
  ft_error('timestamps should be uint64');
end

fid  = fopen_or_error(filename, 'wb', 'ieee-le');

% construct the header
buf  = [];
buf  = [buf sprintf('######## Neuralynx Data File Header\r\n')];
buf  = [buf sprintf('## File Name: %s\r\n', filename)];
buf  = [buf sprintf('## Time Opened: (m/d/y): %s  At Time: %s\r\n', datestr(clock, 'mm/dd/yy'), datestr(clock, 'HH:MM:SS'))];
f = fieldnames(nts.hdr);
for i=1:length(f)
  v = getfield(nts.hdr, f{i});
  switch class(v)
    case 'char'
      buf = [buf sprintf('-%s\t%s\r\n', f{i}, v)];
    case 'double'
      buf = [buf sprintf('-%s\t%g\r\n', f{i}, v)];
    otherwise
      ft_error('unknown class in writing header');
  end
end

% pad the rest of the header with zeros and write it to file
buf((end+1):16384) = 0;
fwrite(fid, buf);

% The format of a clustered electrode record is
%   int64 TimeStamp

fwrite(fid, nts.TimeStamp, 'uint64');

% close the file
fclose(fid);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function YESNO = even(x,order,d)
%EVEN Check if all numbers are even

if (nargin==2) | (nargin==3)
    switch order
        case 'rows'
           YESNO = ~any(rem(x,d),2);
        case 'cols'
            YESNO = ~any(rem(x,d),1);
        otherwise
            error
    end
else
YESNO = nnz(rem(x,2))==0;
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
classdef IpoptApplication < SolverApplication
    % IpoptApplication defines an interface application class for IPOPT solvers
    %
    %
    % @author Ayonga Hereid @date 2016-10-22
    %
    % Copyright (c) 2016, AMBER Lab
    % All right reserved.
    %
    % Redistribution and use in source and binary forms, with or without
    % modification, are permitted only in compliance with the BSD 3-Clause
    % license, see
    % http://www.opensource.org/licenses/bsd-license.php

    %% Private properties
    properties (Access=public)

        % A structure contains the information of objective functions
        %
        % @type struct
        Objective

        % A sturcture contains the information of constraints
        %
        % @type struct
        Constraint

        % The nonlinear programming object
        %
        % @type NonlinearProgram
        Nlp
    end

    %% Public methods
    methods (Access = public)

        function obj = IpoptApplication(nlp, new_opts)
            % The default constructor function
            %
            % Parameters:
            %  nlp: the NLP problem object to be solved
            %  @type NonlinearProgram
            %  new_opts: the new solver options @type struct

            obj = obj@SolverApplication();

            % ipopt options
            options.initialguess = 'typical';

            % default IPOPT options
            options.ipopt.mu_strategy      = 'adaptive';
            options.ipopt.max_iter         = 1000;
            options.ipopt.tol              = 1e-3;
            options.ipopt.linear_solver    = 'ma57';
            options.ipopt.ma57_automatic_scaling = 'yes';
            options.ipopt.linear_scaling_on_demand = 'no';
            options.ipopt.ma57_pre_alloc = 2;
            %             options.ipopt.alpha_for_y = 'bound-mult';
            options.ipopt.recalc_y = 'yes';
            options.ipopt.recalc_y_feas_tol = 1e-3;
            %             options.ipopt.bound_relax_factor = 1e-3;
            %             options.ipopt.fixed_variable_treatment = 'RELAX_BOUNDS';
            %             options.ipopt.derivative_test = 'first-order';
            %             options.ipopt.point_perturbation_radius = 0;
            %             options.ipopt.derivative_test_perturbation = 1e-8;
            % options.ipopt.derivative_test_print_all = 'yes';
            options.ipopt.print_level = 5;
            options.ipopt.print_timing_statistics = 'yes';
            %             options.ipopt.print_options_documentation = 'yes';

            if nlp.Options.DerivativeLevel == 2
                % user-defined Hessian function is provide
                options.ipopt.hessian_approximation = 'exact';
            else
                options.ipopt.hessian_approximation = 'limited-memory';
                options.ipopt.limited_memory_update_type = 'bfgs';  % {bfgs}, sr1
                options.ipopt.limited_memory_max_history = 10;  % {6}
            end

            if nargin > 1
                options.ipopt = struct_overlay(options.ipopt,new_opts,{'AllowNew',true});
            end

            if exist(['sparse2.',mexext],'file') == 3
                options.UseMexSparse = true;
            else
                options.UseMexSparse = false;
            end

            obj.Options = options;
            obj.Nlp = nlp;
            obj = initialize(obj);

        end

    end

    % function definitions
    methods

        [obj] = initialize(obj);

        [sol, info] = optimize(obj, x0);

    end

end

author: ayonga

repo: frost-dev

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


text:
function F = fit_metric_field(DiffDir, arange,brange, Q)

% fit_metric_field - fit a field of tensor
%
%    F = fit_metric_field(DiffDir, arange,brange, Q);
%
%   F is a (Q,Q,2,2) field of 2x2 SDP matrices.
%   DiffDir(x,h) is the directional value of the metric at x in direction h.
%
%   Copyright (c) 2014 Gabriel Peyre

F = zeros(Q,Q,2,2);
alist = linspace(arange(1), arange(2), Q);
blist = linspace(brange(1), brange(2), Q);
for i=1:Q
    % progressbar(i,Q);
    for j=1:Q
        ab = [alist(i) blist(j)];
        M = fit_metric(@(h)DiffDir(ab,h));
        F(i,j,:,:) = reshape(M, [1 1 2 2]);
    end
end

end

author: gpeyre


text:
function hyperball_monte_carlo_test02 ( )

%*****************************************************************************80
%
%% HYPERBALL_MONTE_CARLO_TEST02 uses HYPERBALL01_SAMPLE in 6D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 January 2014
%
%  Author:
%
%    John Burkardt
%
  m = 6;

  e_test = [ ...
    0, 0, 0, 0, 0, 0; ...
    1, 0, 0, 0, 0, 0; ...
    0, 2, 0, 0, 0, 0; ...
    0, 2, 2, 0, 0, 0; ...
    0, 0, 0, 4, 0, 0; ...
    2, 0, 0, 0, 2, 2; ...
    0, 0, 0, 0, 0, 6 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  Use the Monte Carlo method to estimate integrals \n' );
  fprintf ( 1, '  over the interior of the unit hyperball in M dimensions.\n' );
  fprintf ( 1, 'n' );
  fprintf ( 1, '  Spatial dimension M = %d\n', m );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '         N' );
  fprintf ( 1, '        1      ' );
  fprintf ( 1, '        U      ' );
  fprintf ( 1, '         V^2   ' );
  fprintf ( 1, '         V^2W^2' );
  fprintf ( 1, '         X^4   ' );
  fprintf ( 1, '         Y^2Z^2' );
  fprintf ( 1, '         Z^6\n' );
  fprintf ( 1, '\n' );

  n = 1;

  while ( n <= 65536 )

    [ x, seed ] = hyperball01_sample ( m, n, seed );

    fprintf ( 1, '  %8d', n );
    for j = 1 : 7
      e(1:m) = e_test(1:m,j);
      value = monomial_value ( m, n, e, x );
      result(j) = hyperball01_volume ( m ) * sum ( value(1:n) ) / n;
      fprintf ( 1, '  %14.6g', result(j) );
    end
    fprintf ( 1, '\n' );

    n = 2 * n;

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '     Exact' );
  for j = 1 : 7
    e(1:m) = e_test(1:m,j);
    result(j) = hyperball01_monomial_integral ( m, e );
    fprintf ( 1, '  %14.6g', result(j) );
  end
  fprintf ( 1, '\n' );

  return
end

author: johannesgerer


text:
%% patchIND
% Below is a demonstration of the features of the |patchIND| function

%%
clear; close all; clc;

%% Syntax
% |[IND_F,IND_V]=patchIND(F,V,formOpt);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function r8mat_is_solution_test ( )

%*****************************************************************************80
%
%% R8MAT_IS_SOLUTION_TEST tests R8MAT_IS_SOLUTION.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 March 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8MAT_IS_SOLUTION_TEST:\n' );
  fprintf ( 1, '  R8MAT_IS_SOLUTION tests whether X is the solution of\n' );
  fprintf ( 1, '  A*X=B by computing the Frobenius norm of the residual.\n' );
%
%  Get random shapes.
%
  i4_lo = 1;
  i4_hi = 10;
  seed = 123456789;
  [ m, seed ] = i4_uniform_ab ( i4_lo, i4_hi, seed );
  [ n, seed ] = i4_uniform_ab ( i4_lo, i4_hi, seed );
  [ k, seed ] = i4_uniform_ab ( i4_lo, i4_hi, seed );
%
%  Get a random A.
%
  r8_lo = -5.0;
  r8_hi = +5.0;
  [ a, seed ] = r8mat_uniform_ab ( m, n, r8_lo, r8_hi, seed );
%
%  Get a random X.
%
  r8_lo = -5.0;
  r8_hi = +5.0;
  [ x, seed ] = r8mat_uniform_ab ( n, k, r8_lo, r8_hi, seed );
%
%  Compute B = A * X.
%
  b = a * x;
%
%  Compute || A*X-B||
%
  enorm = r8mat_is_solution ( m, n, k, a, x, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  A is %d by %d\n', m, n );
  fprintf ( 1, '  X is %d by %d\n', n, k );
  fprintf ( 1, '  B is %d by %d\n', m, k );
  fprintf ( 1, '  Frobenius error in A*X-B is %g\n', enorm );

  return
end


text:
function p = extractBounds(p)
if ~isempty(p.F_struc)
    [lb,ub,used_rows_eq,used_rows_lp] = find_lp_bounds(p.F_struc,p.K);
    if ~isempty(used_rows_lp)
        used_rows_lp = used_rows_lp(~any((p.F_struc(p.K.f + used_rows_lp,1+p.nonlinear)),2));
        if ~isempty(used_rows_lp)
            lower_defined = find(~isinf(lb));
            if ~isempty(lower_defined)
                p.lb(lower_defined) = max(p.lb(lower_defined),lb(lower_defined));
            end
            upper_defined = find(~isinf(ub));
            if ~isempty(upper_defined)
                p.ub(upper_defined) = min(p.ub(upper_defined),ub(upper_defined));
            end
            p.F_struc(p.K.f + used_rows_lp,:)=[];
            p.K.l = p.K.l - length(used_rows_lp);
        end
    end
    if ~isempty(used_rows_eq)
        used_rows_eq = used_rows_eq(~any(full(p.F_struc(used_rows_eq,1+p.nonlinear)),2));
        if ~isempty(used_rows_eq)
            lower_defined = find(~isinf(lb));
            if ~isempty(lower_defined)
                p.lb(lower_defined) = max(p.lb(lower_defined),lb(lower_defined));
            end
            upper_defined = find(~isinf(ub));
            if ~isempty(upper_defined)
                p.ub(upper_defined) = min(p.ub(upper_defined),ub(upper_defined));
            end
            p.F_struc(used_rows_eq,:)=[];
            p.K.f = p.K.f - length(used_rows_eq);
        end
    end
end


text:
function RtOut = transformCameraRt(RtIn)

RtOut = [RtIn(1:3,1:3)', - RtIn(1:3,1:3)'* RtIn(1:3,4)];

author: jianxiongxiao

repo: ProfXkit


text:
function obj_val = evaluate_objective(objective_function,scores,classf, ...
				  prior)
% Returns the result of the objective function evaluated on the
% scores.
%
% Inputs:
%   objective_function: a function handle to the objective function
%                       to feed the scores into
%   scores: length T vector of scores to be evaluated where T is
%           the number of trials
%   classf: length T vector with entries +1 for target scores; -1
%           for non-target scores
%   prior: the prior (given to the system that produced the scores)
%
% Outputs
%   obj_val: the value returned by the objective function

if nargin==0
    test_this();
    return;
end

if ~exist('objective_function','var') || isempty(objective_function)
    objective_function = @(w,T,weights,logit_prior) cllr_obj(w,T,weights,logit_prior);
end

logit_prior = logit(prior);
prior_entropy = objective_function([0;0],[1,-1],[prior,1-prior],logit_prior);

ntar = length(find(classf>0));
nnon = length(find(classf<0));
N = nnon+ntar;

weights = zeros(1,N);

weights(classf>0) = prior/(ntar*prior_entropy);
weights(classf<0) = (1-prior)/(nnon*prior_entropy);

obj_val = objective_function(scores,classf,weights,logit_prior);

end

function test_this()
num_trials = 20;
scores = randn(1,num_trials);
classf = [ones(1,num_trials/2),-ones(1,num_trials/2)];
prior = 0.5;
res = evaluate_objective(@cllr_obj,scores,classf,prior)
end

author: nesl


text:
function b = r83p_vxm ( n, a, x )

%*****************************************************************************80
%
%% R83P_VXM multiplies a vector by a R83P matrix.
%
%  Discussion:
%
%    The R83P storage format stores a periodic tridiagonal matrix as
%    a 3 by N array, in which each row corresponds to a diagonal, and
%    column locations are preserved.  The matrix value
%    A(1,N) is stored as the array entry A(3,N), and the matrix value
%    A(N,1) is stored as the array entry A(1,1).
%
%  Example:
%
%    Here is how a R83P matrix of order 5 would be stored:
%
%      A51 A12 A23 A34 A45
%      A11 A22 A33 A44 A55
%      A21 A32 A43 A54 A15
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be at least 3.
%
%    Input, real A(3,N), the R83P matrix.
%
%    Input, real X(N), the vector to be multiplied by A.
%
%    Output, real B(N), the product X * A.
%
  b(1) =   a(1,1)   * x(n)   + a(2,1) * x(1) + a(3,1)   * x(2);

  for i = 2 : n-1
    b(i) = a(1,i)   * x(i-1) + a(2,i) * x(i) + a(3,i)   * x(i+1);
  end

  b(n) =   a(1,n)   * x(n-1) + a(2,n) * x(n) + a(3,n)   * x(1);

  return
end


text:
function [Xhat,XXstar] = computeStateResponses(XX,ehmm,Gamma,states)

regressed = sum(ehmm.train.Sind==1,1)>0;
K = size(Gamma,2); np = size(XX,2); ndim = size(ehmm.train.S,1);

if nargin < 4 || isempty(states), states = 1:K+1; end

Gamma = [Gamma prod(1-Gamma,2) ];
Gamma = rdiv(Gamma,sum(Gamma,2));

XXstar = zeros(size(XX,1),np * length(states));
ind = false((K+1)*np,1);

for ik = 1:length(states)
    k = states(ik);
    XXstar(:,(1:np) + (ik-1)*np) = bsxfun(@times, XX, Gamma(:,k));
    ind( (1:np) + (k-1)*np) = true;
end

Xhat = zeros(size(XX,1),ndim);
for n = find(regressed)
    Xhat(:,n) = XXstar * ehmm.state_shared(n).Mu_W(ind);
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function xcovar=covarf(x,cov_wsize)

%%% covarf: computes covariance of a signal.
%%% covarf(X), if X is a vector, returns the variance.
%%% For matrices, where each row is an observation, and each column a variable,
%%% covarf(X) is the covariance matrix.
%%% diag(covarf(X)) is a vector of variances for each column.
%%% sqrt(diag(covarf(X))) is a vector of standard deviations.
%%% wsize should be odd.

%%% Author : B. K. SHREYAMSHA KUMAR
%%% Created on 15-02-2012.
%%% Updated on 15-02-2012.

tr=x-repmat(mean(x),cov_wsize,1);
xcovar=tr'*tr/(cov_wsize-1);

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
%% Author: epokh
%% Website: www.epokh.org/drupy
%% This software is under GPL

%%This example show the use of DH parameters, omogeneous transformations
%%jacobian,static force analysis and plot functions.

%%An example of a 6 degree of freedom robot: 6 revolute joints
%% Solve a DH problem for forward kinematics
%%Joint variables
clf
theta1=-30;
theta2=-25;
theta3=-50;
theta4=-45;
theta5=-20;
theta6=95;

%%link length
a3=2.5;
a4=2.1;
a5=1.8;

%%First way to calculate the forward kinematics
T01=DHmatrix(theta1,0,0,0);
T12=DHmatrix(0,0,0,-90);
T23=DHmatrix(theta2,0,a3,0);
T34=DHmatrix(theta3,0,a4,0);
T45=DHmatrix(theta4,0,a5,0);
T56=RotZ(theta5)*RotX(theta6);

Tuh1=T01*T12*T23*T34*T45*T56;

%%A second fast way to calculate it
T01=RotZ(theta1);
T12=RotX(-90);
T23=RotZ(theta2)*Tras(a3,0,0);
T34=RotZ(theta3)*Tras(a4,0,0);
T45=RotZ(theta4)*Tras(a5,0,0);
T56=RotZ(theta5)*RotX(theta6);

% %%This shows how to use the plot system
% plotT(T01);
% pause(2);
% plotT(T01*T12);
% pause(2);
% plotT2(T01*T12,T01*T12*T23);
% pause(2);
% plotT2(T01*T12*T23,T01*T12*T23*T34);
% pause(2);
% plotT2(T01*T12*T23*T34,T01*T12*T23*T34*T45);
% pause(2);
% plotT2(T01*T12*T23*T34*T45,T01*T12*T23*T34*T45*T56);

Tuh2=T01*T12*T23*T34*T45*T56;

%%now calculate the jacobian
J6=jacobianT6([T01,T12,T23,T34,T45,T56],['R','R','R','R','R','R']);

F=[2.1;3.2;4.5;10.1;0.5;0.07];
T=staticForce(J6,F);

author: Sable

repo: mcbench-benchmarks


text:
% Rolls NN params vector into the matrix.
function nn_params_rolled = nn_params_roll(nn_params_unrolled, layers)
    % Calculate the total number of layers.
    L = length(layers);

    nn_params_rolled = {};
    unrolled_shift = 0;

    for layer_number=1:(L-1)
        L_in = layers(layer_number);
        L_out = layers(layer_number + 1);

        params_width = L_in + 1;  % We need to remember about bias unit.
        params_height = L_out;
        params_volume = params_height * params_width;

        % We need to remember about bias units when rolling up params.
        layer_params_unrolled = nn_params_unrolled((unrolled_shift + 1):(unrolled_shift + params_volume), 1);
        nn_params_rolled{layer_number} = reshape(layer_params_unrolled, params_height, params_width);

        unrolled_shift = unrolled_shift + params_volume;
    end
end

author: trekhleb

repo: machine-learning-octave