text:
function bms = spm_cfg_bms_map
% Configuration file for BMS interface
%__________________________________________________________________________
% Copyright (C) 2008-2016 Wellcome Trust Centre for Neuroimaging

% Maria Joao Rosa
% $Id: spm_cfg_bms_map.m 6952 2016-11-25 16:03:13Z guillaume $

%--------------------------------------------------------------------------
% dir Directory
%--------------------------------------------------------------------------
dir         = cfg_files;
dir.tag     = 'dir';
dir.name    = 'Directory';
dir.help    = {['Select the directory where the files containing the '...
               'results from BMS (BMS.mat) will be written.']};
dir.filter  = 'dir';
dir.ufilter = '.*';
dir.num     = [1 1];

%--------------------------------------------------------------------------
% mod_map Log evidence maps
%--------------------------------------------------------------------------
mod_map         = cfg_files;
mod_map.tag     = 'mod_map';
mod_map.name    = 'Models';
mod_map.help    = {['Specify the log. evidence map for each model. '...
                    'Log-evidence maps should be specified '...
                    'in the same order for each subject and session.']};
mod_map.filter  = 'image';
mod_map.ufilter = '.*';
mod_map.num     = [1 Inf];

%--------------------------------------------------------------------------
% sess_map Sessions (Maps)
%--------------------------------------------------------------------------
sess_map      = cfg_branch;
sess_map.tag  = 'sess_map';
sess_map.name = 'Session';
sess_map.val  = {mod_map };

%--------------------------------------------------------------------------
% subj_dcm Subject (Maps)
%--------------------------------------------------------------------------
subj_map         = cfg_repeat;
subj_map.tag     = 'subj_map';
subj_map.name    = 'Subject';
subj_map.values  = {sess_map };

%--------------------------------------------------------------------------
% map Data
%--------------------------------------------------------------------------
map         = cfg_repeat;
map.tag     = 'map';
map.name    = 'Data';
map.help    = {['Select the log. evidence maps for each '...
               'model, session and subject.']}';
map.values  = {subj_map };
map.num     = [1 Inf];

%--------------------------------------------------------------------------
% mod_name Name
%--------------------------------------------------------------------------
mod_name         = cfg_entry;
mod_name.tag     = 'mod_name';
mod_name.name    = 'Name';
mod_name.help    = {'Specify name for each model (optional).'};
mod_name.strtype = 's';
mod_name.num     = [0 Inf];
mod_name.val     = {''};

%--------------------------------------------------------------------------
% name_mod Name models
%--------------------------------------------------------------------------
name_mod         = cfg_repeat;
name_mod.tag     = 'name_mod';
name_mod.name    = 'Name models';
name_mod.help    = {'Specify name for each model (optional).'}';
name_mod.values  = {mod_name };
name_mod.num     = [0 Inf];

%--------------------------------------------------------------------------
% method_maps Inference Method (maps)
%--------------------------------------------------------------------------
method_maps         = cfg_menu;
method_maps.tag     = 'method_maps';
method_maps.name    = 'Inference method';
method_maps.help    = {['Specify inference method: random effects '...
                   '(2nd-level, RFX) or fixed effects (1st-level, FFX) analysis. '...
                   'RFX uses a Variational Bayes approach.']};
method_maps.labels  = {
                  'Fixed effects (FFX)'
                  'Random effects (RFX)'
}';
method_maps.values  = {
                  'FFX'
                  'RFX'
}';

% %--------------------------------------------------------------------------
% % priors Priors
% %--------------------------------------------------------------------------
% priors         = cfg_menu;
% priors.tag     = 'priors';
% priors.name    = 'Priors';
% priors.help    = {['Specify priors for family-level inference (RFX only).
% '...
%                    'Options: ''Family'' sets alpha0=1 for each family '...
%                    'while ''Model'' sets alpha0=1 for each model (not '...
%                    'advised).']};
% priors.labels  = {
%                   'Model'
%                   'Family'
% }';
% priors.values  = {
%                   'M-unity'
%                   'F-unity'
% }';
% priors.val      = {'F-unity'};

%--------------------------------------------------------------------------
% out_file Output files
%--------------------------------------------------------------------------
out_file         = cfg_menu;
out_file.tag     = 'out_file';
out_file.name    = 'Output files (RFX)';
out_file.help    = {['Specify which output files to save (only valid for'...
                     'RFX analyses). ']...
                     ''...
                    ['Default option (and faster option): '...
                     'PPM = xppm.<ext> (Expected Posterior Probability Maps) '...
                     'for each model ie. posterior mean.']...
                     ''...
                    ['Second option: PPM + EPM = xppm.<ext> + '...
                     'epm.<ext> (Expected Posterior Probability '...
                     'Maps + Exceedance Probability Maps) for each model.']...
                     ''...
                    ['Third option: PPM + EPM + Alpha = xppm.<ext> + '...
                     'epm.<ext> + alpha.<ext> (PPM, EPM and Map of Dirichlet '...
                     'Parameters) for each model.']};
out_file.labels  = {
                   'PPM'
                   'PPM + EPM'
                   'PPM + EPM + Alpha'

}';
out_file.values  = {
                  0
                  1
                  2
}';
out_file.val     = {0};

%--------------------------------------------------------------------------
% mask Mask Image
%--------------------------------------------------------------------------
mask         = cfg_files;
mask.tag     = 'mask';
mask.name    = 'Mask Image';
mask.help    = {['Specify an image for explicitly masking the analysis. '...
                '(optional). '...
                'A sensible option here is to use a segmention of '...
                'structural images to specify a within-brain mask. '...
                'If you select that image as an explicit mask then only '...
                'those voxels in the brain will be analysed. This both '...
                'speeds the inference process and restricts BMS to '...
                'within-brain voxels. Alternatively, if such structural '...
                'images are unavailble or no masking is required, then '...
                'leave this field empty.']};
mask.filter  = 'image';
mask.ufilter = '.*';
mask.val     = {{''}};
mask.num     = [0 1];

%--------------------------------------------------------------------------
% nsamp Number of samples
%--------------------------------------------------------------------------
nsamp         = cfg_entry;
nsamp.tag     = 'nsamp';
nsamp.name    = 'Number of samples';
nsamp.help    = {['Number of samples used to compute exceedance '...
                  'probabilities (default: 1e6). '...
                  'To make computations faster reduce the number of '...
                  'samples when number of models is bigger than 3.']};
nsamp.strtype = 's';
nsamp.num     = [1 Inf];
nsamp.val     = {'1e6'};

%--------------------------------------------------------------------------
% file BMS.mat
%--------------------------------------------------------------------------
file         = cfg_files;
file.tag     = 'file';
file.name    = 'BMS.mat';
file.help    = {['Specify the BMS (.mat) file obtained from previous BMS '...
               'analysis (optional). Leave field empty to work on '...
               'serial mode.']};
file.filter  = 'mat';
file.ufilter = '.*';
file.val     = {{''}};
file.num     = [0 1];

%--------------------------------------------------------------------------
% img Map to display
%--------------------------------------------------------------------------
img         = cfg_files;
img.tag     = 'img';
img.name    = 'Map to display';
img.help    = {['Specify map obtained from BMS Maps '...
               '(optional). Leave field empty to work on serial mode.']};
img.filter  = 'image';
img.ufilter = '.*';
img.val     = {{''}};
img.num     = [0 1];

%--------------------------------------------------------------------------
% thres Probability Threshold
%--------------------------------------------------------------------------
thres         = cfg_entry;
thres.tag     = 'thres';
thres.name    = 'Probability threshold';
thres.help    = {['Specify the probability threshold to apply to the '...
                 'image (optional). Leave field empty to work on '...
                 'serial mode.']};
thres.strtype = 'r';
thres.num     = [0 Inf];
thres.val     = {[]};

%--------------------------------------------------------------------------
% k Extent threshold
%--------------------------------------------------------------------------
k         = cfg_entry;
k.tag     = 'k';
k.name    = 'Extent threshold';
k.help    = {['Specify extent threshold (minimum number of voxels '...
                 'per cluster).']};
k.strtype = 'w';
k.num     = [0 Inf];
k.val     = {[]};

%--------------------------------------------------------------------------
% scale Map Scale
%--------------------------------------------------------------------------
scale         = cfg_menu;
scale.tag     = 'scale';
scale.name    = 'Map scale';
scale.help    = {['Specify scale to display maps (optional). Default: '...
                 'empty field to work on serial mode. Other options: '...
                 '''None'' will display image with original scale and '...
                 '''Log-odds'' will display image in a log-odds '...
                 ' scale (in this case image should be a '...
                 'probability map).']};
scale.labels  = {
                  'Empty'
                  'None'
                  'Log-odds'
}';
scale.values  = {
                  []
                  0
                  1
}';
scale.val     = {[]};

%--------------------------------------------------------------------------
% bms_map_inf BMS: Maps (Inference), output is BMS map
%--------------------------------------------------------------------------
bms_map_inf      = cfg_exbranch;
bms_map_inf.tag  = 'inference';
bms_map_inf.name = 'BMS: Maps (Inference)';
bms_map_inf.val  = {dir map name_mod method_maps out_file mask nsamp };
bms_map_inf.help = {'Bayesian Model Selection for Log-Evidence Maps.'...
    ''...
    ['Input: log-evidence maps for each model, session and '...
    'subject. Note that there must be identical numbers of models for '...
    'all sessions, and identical numbers of sessions for all '...
    'subjects.']...
    ''...
    ['Output: For the fixed effects analysis, posterior probability maps '...
    'are created for each model. '...
    'For the random effects analysis, expected posterior probability '...
    'and exceedance probability (i.e. the probability that this model '...
    'is more likely than any other model) maps are created for each '...
    'model. If there are multiple sessions per subject, the random '...
    'effects analysis operates on the subject-specific sums of log '...
    'evidences across sessions. In addition, a BMS.mat file will be save '...
    'in the specified directory for both methods']};
bms_map_inf.prog = @spm_run_bms_map;
bms_map_inf.vout = @vout;

%--------------------------------------------------------------------------
% bms_map_vis BMS: Maps (Results), visualisation of BMS Maps results
%--------------------------------------------------------------------------
bms_map_vis      = cfg_exbranch;
bms_map_vis.tag  = 'results';
bms_map_vis.name = 'BMS: Maps (Results)';
bms_map_vis.val  = {file img thres k scale};
bms_map_vis.help = {['Bayesian Model Selection Maps (Results).'...
                    'Show results from BMS Maps (Inference).']};
bms_map_vis.prog = @spm_run_bms_vis;

%--------------------------------------------------------------------------
% bms Bayesian Model Selection
%--------------------------------------------------------------------------
bms         = cfg_choice;
bms.tag     = 'bms_map';
bms.name    = 'Bayesian Model Selection';
bms.help    = {['Bayesian Model Selection for group studies (fixed '...
               'effects and random effects analysis).']};
bms.values  = { bms_map_inf bms_map_vis };

%==========================================================================
function dep = vout(varargin)
% Output file names will be saved in a struct with field .files
dep(1)            = cfg_dep;
dep(1).sname      = 'BMS.mat File';
dep(1).src_output = substruct('.','files');
dep(1).tgt_spec   = cfg_findspec({{'filter','mat','strtype','e'}});


text:
function y = det_inv( varargin )

%DET_INV determinant of the inverse of an SPD matrix.
%   For a square matrix X, DET_INV(X) returns 1.0./DET(X) if X is symmetric
%   (real) or Hermitian (complex) and positive defininte, and +Inf otherwise.
%
%   This function can be used in many convex optimization problems that call
%   for LOG(DET(X)) instead. For example, if the objective function is
%      maximize(logdet(X))
%   then it can be replaced with
%      maximize(-det_inv(X))
%   and the same optimal point will be produced.
%
%   DET_INV(X,p) computes DET_INV(X)^p. p must be a positive real scalar.
%
%   Disciplined convex programming information:
%       DET_INV(X) is convex and nonmonotonic in X; therefore, when used in
%       CVX specifications, its argument must be affine.

persistent P
if isempty( P ),
    P.nargs     = 2;
    P.args      = @det_inv_args;
    P.empty     = 1;
    P.constant  = @det_inv_diag;
    P.diag      = @det_inv_diag;
    P.affine    = @det_inv_aff;
    P.structure = 'psdeig';
end
y = cvx_matrix_op( P, varargin );

function [ X, p ] = det_inv_args( X, p )
if isempty( p ),
    p = 1;
elseif ~( isnumeric(p) && isreal(p) && numel(p)==1 && p>=0 ),
    cvx_throw( 'Second argument must be a positive scalar.' );
end

function y = det_inv_diag( D, p )
y = prod_inv( D, p );

function y = det_inv_aff( X, p )
cvx_begin sdp
    epigraph variable z nonnegative_
    variable Y(n,n) lower_triangular complex_if(X)
    prod_inv( real(diag(Y)), 2*p ) <= z;
    [ diag( D ), Y' ; Y, X ] >= 0;
cvx_end

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

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [u, V, exitflag, output] = snd_solveOptimalControlProblem (snd,  varargin)
%UNTITLED2 Summary of this function goes here
%   solves the optimal control problem of the

    % Set control and linear bounds
    A = [];
    b = [];
    Aeq = [];
    beq = [];
    lb = [];
    ub = [];
    for k=1:snd.horizon %Aggregation
        [Anew, bnew, Aeqnew, beqnew, lbnew, ubnew] = ...
               snd.l_constraints( k, snd.net_load, snd.battery, snd.u0_ref);
        A = blkdiag(A,Anew);
        b = [b, bnew];
        Aeq = blkdiag(Aeq,Aeqnew);
        beq = [beq, beqnew];
        lb = [lb, lbnew];
        ub = [ub, ubnew];
    end

    % Solve optimization problem
    [u, V, exitflag, output] = fmincon( @(u) snd.costfunction( snd, u ), ...
        snd.u0, A, b, Aeq, beq, lb, ub, ...
        @(u) snd.nonlinearconstraints(snd, u ), snd.option);
end

author: juchengquan


text:
function [altitude,gasTable,t,d]=NRLMSISE00Alt4Pres(dayOfYear,secondOfTheDay,pressure,latLon,AP,F107,F107A)
%%NRLMSISE00ALT4PRES Given the pressure at a particular latitude and
%                   longitude and time, obtain the approximate ellipsoidal
%                   height using the NRLMSISE-00 atmospheric model. This
%                   function also returns a standard temperature and the
%                   constitutents of the atmosphere, assuming dry air. The
%                   pressure returned is consistent with the function
%                   NRLMSISE00GasTemp, but will not be perfectly
%                   consistent with the function standardAtmosParam.
%
%INPUTS: dayOfYear The integer day of the year in the Gregorian calendar
%                  in universal coordinated time (UTC). Counting starts at
%                  1. the resolution of the model is not sufficient for it
%                  to matter whether 365 or 366 is given at the day if it
%                  isn't/is a leap year.
%   secondOfTheDay The second of the day. This starts at zero. The
%                  resolution of the model is not high enough for leap
%                  seconds to matter, so values above 86400.0 are just
%                  clipped to 86400.0.
%         pressure The measured atmospheirc pressure in Pascals.
%           latLon WGS-84 latitude and longitude of the point where the
%                  pressure is measured.
%               Ap An optional parameter specifying the average daily
%                  geomagnetic index. If omitted, the value 4.0, which is
%                  suitable for models below 90km is used. The index at a
%                  particular time can be obtained from a service of the
%                  International Service of Geomagnetic Indices at
%                  http://www-app3.gfz-potsdam.de/kp_index/
%                  A scalar value of AP indicates standard mode. If a 7X1
%                  vector is passed, then it is assumed that storm mode is
%                  desired. In storm mode, Ap(1) is the daily Ap index.
%                  Ap(2) through Ap(5) are 3 hours Ap indices respectively
%                  for the current time and 3, 6, and 9 hours before the
%                  current time. Ap(6) is the average of 8 3-hour Ap
%                  indices from 12 to 33 hours prior to the current time
%                  and Ap(7) is the average of 8 3-hour Ap indices from 36
%                  to 59 hours prior to the current time.
%             F107 An optional parameter specifying the 10.7cm solar radio
%                  flux for the previous day. If omitted, a value of 150
%                  is used, which should be suitable for models below
%                  90km. Values of the index can be obtained from
%                  http://www.swpc.noaa.gov/ftpdir/indices/quar_DSD.txt
%            F107A An optional parameter specifying the 81 day average of
%                  the 10.7cm radio solar flux. If omitted, a value of 150
%                  is used, which should be suitable for models below 90km.
%
%OUTPUTS: altitude The approximate ellipsoidal altitude associated with
%                  the pressure given on the input in meters.
%         gasTable A cell array of constituent elements of the atmosphere
%                  and their number densities. gasTable{i,1} is a string
%                  listing the name of the ith constituent element of the
%                  atmosphere. This can take the values
%                  'He' Helium
%                  'O'  Elemental Oxygen
%                  'N2' Nitrogen
%                  'O2' Oxygen
%                  'Ar' Argon
%                  'H'  Elemental Hydrogen
%                  'N'  Elemental Nitrogen
%                  'O*' Anomalous Oxygen
%                  gasTable{i,2} is the corresponding number density of
%                  the ith constituent element in units of atoms per cubic
%                  meter.
%                t Temperature variables. t(1) is the exospheric
%                  temperature and t(2) is the temperature at altitude.
%                  The temperature is in units of Kelvin.
%                d A 9X1 vector containing the same information as
%                  gasTable but without labels. The entries in d are
%                  d(1) - HE NUMBER DENSITY
%                  d(2) - O NUMBER DENSITY
%                  d(3) - N2 NUMBER DENSITY
%                  d(4) - O2 NUMBER DENSITY
%                  d(5) - AR NUMBER DENSITY
%                  d(6) - TOTAL MASS DENSITY without d(9)
%                  d(7) - H NUMBER DENSITY
%                  d(8) - N NUMBER DENSITY
%                  d(9) - Anomalous oxygen NUMBER DENSITY(M-3)
%                  where all number densities are in units of atoms per
%                  cubic meter and the units of d(6) are kilograms per
%                  cubic meter. Note that d(6) is computed using a
%                  definition of the atomic mass unit and definitions of
%                  the atomic masses that are not necessarily the most up
%                  to date.
%
%This function is essentially a wrapper for the C-code implementation of
%the  NRLMSISE-00 atmospheric model using appropriate functions to put the
%time and location in the correct format. The model is described in [1].
%
%The algorithm can be compiled for use in Matlab  using the
%CompileCLibraries function.
%
%The algorithm is run in Matlab using the command format
%[altitude,gasTable,t,d]=NRLMSISE00Alt4Pres(dayOfYear,secondOfTheDay,pressure,latLon);
%or
%[altitude,gasTable,t,d]=NRLMSISE00Alt4Pres(dayOfYear,secondOfTheDay,pressure,latLon,AP);
%or
%[altitude,gasTable,t,d]=NRLMSISE00Alt4Pres(dayOfYear,secondOfTheDay,pressure,latLon,AP,F107);
%or
%[altitude,gasTable,t,d]=NRLMSISE00Alt4Pres(dayOfYear,secondOfTheDay,pressure,latLon,AP,F107,F107A);
%
%EXAMPLE:
% dayOfYear=172;
% secondOfTheDay=29000;
% latLon=[60;-70]*(pi/180);
% F107=150;
% F107A=150;
% AP=4;
% pressure=1;%1 Pascal pressure.
% [altitude,gasTable,t,d]=NRLMSISE00Alt4Pres(dayOfYear,secondOfTheDay,pressure,latLon,AP,F107,F107A)
% %One can verify that the altitude is consistent with the standard model:
% latLonAlt=[latLon;altitude];
% [gasTable1,t1,d1]=NRLMSISE00GasTemp(dayOfYear,secondOfTheDay,latLonAlt,AP,F107,F107A)
% %One will see that t1 and d1 match t and d.
%
%REFERENCES:
%[1] J. M. Picone, A. E. Hedin, D. P. Drob, and A. C. Aikin, "NRLMSISE-00
%    empirical model of the atmosphere: Statistical comparisons and
%    scientific issues," Journal of Geophysical Research: Space Physics,
%    vol. 107, no. A12, Dec. 2002.
%
%June 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

error('This function is only implemented as a mexed C or C++ function. Please run CompileCLibraries.m to compile the function for use.')

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [Ystar,Xstar,Tstar]=doprior(Y,X,n,m,p,T,ar,arvar,lambda1,lambda3,lambda4,priorexo)

% generate Yd, using (XXX)
%Yd=[diag(ar(1:n,1).*arvar/lambda1);zeros(n*(p-1),n);zeros(m,n);diag(arvar)];

Yd=[diag(ar(1:n,1).*arvar/lambda1);zeros(n*(p-1),n);(priorexo./(lambda1.*lambda4))';diag(arvar)];

% generate Xd, using (XXX)
Jp=diag([1:p].^lambda3);
Xd=[kron(Jp,diag(arvar/lambda1)) zeros(n*p,m);zeros(m,n*p) diag(1./(lambda1*lambda4(1,:)));zeros(n,n*p) zeros(n,m)]; % error if m is equal to zero
% Compute Td, using (XXX)
Td=n*(p+1)+m;

% finally generate Ystar, Xstar and Tstar
Ystar=[Y;Yd];
Xstar=[X;Xd];
Tstar=T+Td;

author: european-central-bank

repo: BEAR-toolbox


text:
function [inc]=look2inc(la,height,lat)
%LOOK2INC look angle to incidence angle
%
%    [INCIDENCE]=LOOK2INC(LOOK_ANGLE,HEIGHT,LATITUDE)
%
%    LOOK_ANGLE = Look angle (radians) (can be vector or matrix)
%    HEIGHT     = Height of satellite (m)
%    LATITUDE   = mean latitude of ground (degrees)
%
%     Copyright (C) 2015  Bekaert David - University of Leeds
%     Email: eedpsb@leeds.ac.uk or davidbekaert.com
%     With permission from Andy Hooper
%
%     This program is free software; you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation; either version 2 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License along
%     with this program; if not, write to the Free Software Foundation, Inc.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
% by Andrew Hooper, 2010

if nargin<3
   lat=40
end

lat1=lat*pi/180;

WGS84_A=6378137.0; % semimajor axis wgs84
WGS84_B=6356752.314; % semiminor axis wgs84
Re=WGS84_A*WGS84_B./sqrt(WGS84_A^2*sin(lat1).^2+WGS84_B^2*cos(lat1).^2);

a=Re+height; % Earth centre to satellite
inc=la;

[R] = fminsearch(@(p) (a^2+p^2-2*a*p*cos(mean(la(:)))-Re^2)^2,[600])
for i=1:length(la(:))
    R = fminsearch(@(p) (a^2+p^2-2*a*p*cos(la(i))-Re^2)^2,R);
    inc(i)=pi-acos((Re^2+R^2-a^2)/2/Re/R);
end

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
function [R h] = rmCompareModelsGUI_paramXCorr(M, plotFlag);
% Compute the cross-correlation across parameter estimates for the RM compare
% models GUI.
%
%   [R h] = rmCompareModelsGUI_paramXCorr([M=get from cur figure], [plotFlag=1]);
%
% This correlation matrices, in which the basic pRF parameters
% (x0, y0, sigma, polar angle, eccentricity, and variance
% explained) are correlated across the loaded models. The correlation
% matrices are returned in the matrix R.
%
%
% OUTPUTS:
% R is a 3D matrix of size nModels by nModels by 6. The six slices are
% cross-correlation matrices for, respectively:
%	1) x0
%	2) y0
%	3) sigma major
%	4) polar angle
%	5) eccentricity
%	6) variance explained.
%
% For each cross-correlation matrix, the entry R(m,n) is the Pearson's
% correlation coefficient (R) between the pattern across all voxels of the
% relevant parameter, between model m and model n. The diagonal values
% R(n,n) are all set to 1, since each set of parameters is correlated to
% itself.
%
% INPUTS:
%
% M: compare models GUI data structure. [Default is to get it from the
% current figure.]
%
% plotFlag: indicates whether to show these matrices in a new figure. If 1 or 2,
% will return the handle of the new figure(s) in h.
%	0: don't plot
%	1: plot color-coded images of the correlation matrices
%	2: plot scatter plots for each comparison. (Will create a separate
%	figure for each pRF parameter, with the scatter plots arranged as per
%	the lower quadrant of the cross-correlation matrices.)
%
%
% ras, 04/2009.
if notDefined('M'),		M = get(gcf, 'UserData');			end
if notDefined('plotFlag'),	plotFlag = 1;					end

R = [];
h = [];

%% compute the correlation coefficients
fields = {'x0' 'y0' 'sigma' 'pol' 'ecc' 'varexp'};
for z = 1:6
	f = fields{z};  % field name

	data =  reshape( [M.(f){:}], [M.nVoxels M.nModels] );

	R(:,:,z) = corrcoef(data);
end

%% visualize if selected
if plotFlag==1
	%% plot color-coded images of the X-corr matrices
	nm = sprintf('Cross-Model Correlations %s', M.roi.name);
	h = figure('Color', 'w', 'Name', nm);

	for z = 1:6
		subplot(2, 3, z);

		drawXCorrMatrix( R(:,:,z) );

		% label each entry in the lower-left-hand plot
		if z==4
			axis on
			set(gca, 'Box', 'off', 'XTick', 1:M.nModels, ...
				'YTick', 1:1:M.nModels, 'FontSize', 9)
			xlabel('Model #', 'FontSize', 12);
			ylabel('Model #', 'FontSize', 12);
		end

		title( fields{z}, 'FontSize', 14 );
	end

elseif plotFlag==2
	%% plot multiple figures, with scatter plots for each correlation
	N = M.nModels - 1;  % # rows/columns of subplots in this figure

	for z = 1:6
		f = fields{z};  % field name

		nm = sprintf('Cross-Model Correlations %s', M.roi.name);
		h(z) = figure('Color', 'w', 'Name', nm);

		for x = 1:N
			for y = x+1:N+1
				row = y-1;
				subplot(N, N, [(row-1)*N + x]);
				regressPlot(M.(f){x}, M.(f){y}, 'x=y', 'title');

				if z < 3  % x0 or y0: add x=0 / y=0 line
					AX = axis;
					line(AX(1:2), [0 0], 'Color', [.5 .5 .5]);
					line([0 0], AX(3:4), 'Color', [.5 .5 .5]);
				end

% 				xlabel([M.dtList{x} ' ' f], 'FontSize', 12);
% 				ylabel([M.dtList{y} ' ' f], 'FontSize', 12);
				xlabel([M.params{x}.analysis.pRFmodel ' ' f], 'FontSize', 12);
				ylabel([M.params{y}.analysis.pRFmodel ' ' f], 'FontSize', 12);
            end
		end
	end

end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function matlab_time = unixtime2mat(unix_time);
% unixtime2mat  Converts unix time stamps (seconds since Jan 1, 1970) to
%               Matlab serial date number (decimal days since Jan 1 0000).
%
%               USAGE:
%                      unixtime2mat(unix_time)
%
%
%               The function may not handle leap years or leap seconds
%               appropriately.
%
%               Val Schmidt
%               Center for Coastal and Ocean Mapping
%               2007

unix_epoch = datenum(1970,1,1,0,0,0);
matlab_time = unix_time./86400 + unix_epoch;

author: Sable

repo: mcbench-benchmarks


text:
% SCRIPT TO TEST THE KINEMATICS OF THE 5R robot
%
% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.

close all

%load the robot
%the robot is represented (interiorly) by two planar 2dof arms
robot=load_robot('example','Delta');

adjust_view(robot)

%Represent a final point in space that will be achieved by the end effector
%P=[Px, Py, Pz]=[0.6 0.4 -0.2]

T=eye(4);
T(1,4)=-0.2;%change, for example, to 0.3
T(2,4)=0.3; %change, for example to 0.6
T(3,4)=-0.5;

%find 8 solutions for the INVERSE KINEMATIC in position
q=inversekinematic(robot, T);

drawrobot3d(robot, q(:,8)), pause(2);
drawrobot3d(robot, q(:,7)), pause(2);
drawrobot3d(robot, q(:,6)), pause(2);
drawrobot3d(robot, q(:,5)), pause(2);
drawrobot3d(robot, q(:,4)), pause(2);
drawrobot3d(robot, q(:,3)), pause(2);
drawrobot3d(robot, q(:,2)), pause(2);
drawrobot3d(robot, q(:,1)), pause(2);

% now solve the DIRECTKINEMATIC problem for each of the above computed q
% In this case, only the "down" solution is returned
T=directkinematic(robot,q(:,1))
T=directkinematic(robot,q(:,2))
T=directkinematic(robot,q(:,3))
T=directkinematic(robot,q(:,4))
T=directkinematic(robot,q(:,5))
T=directkinematic(robot,q(:,6))
T=directkinematic(robot,q(:,7))
T=directkinematic(robot,q(:,8))

repo: ARTE


text:
function t = isFactored(disc)
%ISFACTORED   Check if a useful factorization of a VALSDISCRETIZATION is stored.
%   ISFACTORED(DISC) returns true if the discretization has stored LU
%   factors whose dimensions are compatible with the current discretization
%   size.

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

% Assume false:
t = false;

% Check whether any factorisation data is stored:
if ( ~isempty(disc.mldivideData) )
    % Total discretization size is the sum over the pieces, times the
    % number of function variables, plus the number of scalar variables.
    isFun = isFunVariable(disc.source);
    len = sum(isFun)*sum(disc.dimension) + sum(~isFun);

    % Compare length to that of stored data.
    L = disc.mldivideData{1};
    if ( len == size(L, 1) )
        t = true;
    end
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Function that outputs a struct "settings" with the settings used in the
% GNSS-aided INS
%
% Edit: Isaac Skog (skog@kth.se), 2016-09-01
% Revised: Bo Bernhardsson 2018-01-01
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function settings = gnss_imu_local_tan_example_settings()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%              GENERAL PARAMETERS         %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
settings.gnss_outage = 'on';
settings.outagestart = 120;
settings.outagestop = 150;
settings.non_holonomic = 'off';
settings.speed_aiding = 'off';

settings.init_heading = 0*pi/180;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%             FILTER PARAMETERS           %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Process noise covariance (Q)
% Standard deviations, need to be squared
settings.sigma_acc = 0.2;
settings.sigma_gyro = deg2rad(0.5);
settings.sigma_acc_bias = 0.00;
settings.sigma_gyro_bias = deg2rad(0.00);

% GNSS position measurement noise covariance (R)
% Standard deviations, need to be squared
settings.sigma_gps =2/sqrt(3); %[m]
settings.sigma_speed = 1; %[m/s]  Trim here
%settings.sigma_speed = 24; %[m/s]  Trim here
settings.sigma_non_holonomic = 20; %[m/s] Trim here
%settings.sigma_non_holonomic = 3; %[m/s] Trim here

% Initial Kalman filter uncertainties (standard deviations)
settings.factp(1) = 100;                                 % Position [m]
settings.factp(2) = 100;                                  % Velocity [m/s]
settings.factp(3:5) = (pi/180*[50 50 50]');     % Attitude (roll,pitch,yaw) [rad]
settings.factp(6) = 0.02;                               % Accelerometer biases [m/s^2]
settings.factp(7) = deg2rad(2);                     % Gyro biases [rad/s]

settings.gravity  = [0, 0, 9.8184]';        % NED frame 重力为正

end

author: yandld


text:
function v = cvx_vexity( x )

global cvx___
s = x.size_;
if any( s == 0 ),
    v = cvx_zeros( s );
    return
end
p  = cvx___.vexity;
b  = x.basis_;
n  = length( p );
nb = size( b, 1 );
if nb < n,
    p = p( 1 : nb, : );
elseif n < nb,
    p( n+1:nb, : ) = 0;
end
b = b( p ~= 0, : );
if isempty( b ),
    v = cvx_zeros( x.size_ );
    if x.slow_,
        v( isnan( x.basis_( 1, : ) ) ) = NaN;
    end
    return
end
p = nonzeros(p).';
if cvx___.nan_used,
    b = sparse( b );
end
v = full( p * b );
tt = abs( v ) ~= abs( p ) * abs( b );
if x.slow_,
    v( tt | isnan( x.basis_( 1, : ) ) ) = NaN;
else
    v( tt ) = NaN;
end
if ~isreal( x ),
    v( any( imag( x.basis_ ), 1 ) & v ) = NaN;
end
v = sign( v );
v = reshape( v, x.size_ );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
% David Krause, david.krause@ece.queensu.ca
% Queen's University
% July 4, 2006
% Create a two axes plot

% Clean up!
close all
clear all
clc

% Some plot data
x_bottom = [0, 1, 2, 3, 4] * 1e-10;
x_top = [-4, -2, -1, 0, 1] * 1e9;
y_left = [0.5, 0.6, 0.3, -0.1, -0.2] * 1e9;
y_right = [0, 4, -5, 5, 6] * 1e4;

% Create the figure
figure(1);
axes;
plot(x_bottom, y_left, 'k.-');
xlimits = get(gca,'XLim');
ylimits = get(gca,'YLim');
xinc = (xlimits(2)-xlimits(1))/5;
yinc = (ylimits(2)-ylimits(1))/5;
set(gca,'XTick',[xlimits(1):xinc:xlimits(2)],...
        'YTick',[ylimits(1):yinc:ylimits(2)])
xlabel('Time (s)');
ylabel('Intensity (W/m)');
text(0.7e-10, 0e9, '$\bullet$ Black Line');
text(0.7e-10, -0.1e9, '$\circ$ Blue Line');

axes;
plot(x_top, y_right, 'bo-');
set(gca, 'color', 'none', 'YAxisLocation', 'right', ...
                          'XAxisLocation', 'top');
xlimits = get(gca,'XLim');
ylimits = get(gca,'YLim');
xinc = (xlimits(2)-xlimits(1))/5;
yinc = (ylimits(2)-ylimits(1))/5;
set(gca,'XTick',[xlimits(1):xinc:xlimits(2)],...
        'YTick',[ylimits(1):yinc:ylimits(2)]);
xlabel('Time (s)');
ylabel('Widgets ($\mu$ W)');

% Create a structure with the output parameters
st.filename = 'the_figure.tex';
st.comments = 'Hello world';
st.figure_box = [25, 15, 100, 90]; % In mm, the size of the plot area (not including the labels)
st.x_tick_weight = 0.25; % The weight (line thickness) of the ticks
st.x_tick_length = 2; % The length of the ticks
st.x_label_y_offset = 8; % How far (in mm), to offset the xlabel
st.x_ticklabel_y_offset = 3;

st.y_tick_weight = 0.25;
st.y_tick_length = 2;
st.y_label_x_offset = 18;
st.y_ticklabel_x_offset = 2;
st.y_ticklabel_pow10 = 6;

st.plot_line_thickness = 0.35;

% Convert to .tex
figure2latex(1, st);


text:
classdef LineSearch_Backtracking < LineSearch

    properties (Access = private)
       rate
    end

    methods (Access = public)

        function obj = LineSearch_Backtracking(cParams)
            obj.init(cParams);
            obj.rate = cParams.rate;
        end

        function update(obj)
            obj.value = obj.value*obj.rate;
            obj.nTrials = obj.nTrials + 1;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%SerialLink.friction Friction force
%
% TAU = R.friction(QD) is the vector of joint friction forces/torques for the
% robot moving with joint velocities QD.
%
% The friction model includes:
% - Viscous friction which is a linear function of velocity.
% - Coulomb friction which is proportional to sign(QD).
%
% Notes::
% - The friction value should be added to the motor output torque, it has a
%   negative value when QD>0.
% - The returned friction value is referred to the output of the gearbox.
% - The friction parameters in the Link object are referred to the motor.
% - Motor viscous friction is scaled up by G^2.
% - Motor Coulomb friction is scaled up by G.
% - The appropriate Coulomb friction value to use in the non-symmetric case
%   depends on the sign of the joint velocity, not the motor velocity.
% - The absolute value of the gear ratio is used.  Negative gear ratios are
%   tricky: the Puma560 has negative gear ratio for joints 1 and 3.
%
% See also Link.friction.

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function  tau = friction(robot, qd)

	L = robot.links;

    tau = zeros(1,robot.n);
    if robot.issym
        tau = sym(tau);
    end

	for j=1:robot.n
		tau(j) = L(j).friction(qd(j));
	end

author: petercorke

repo: robotics-toolbox-matlab


text:
function marginal = marginal_nodes(engine, b, nodes, t, add_ev)
% MARGINAL_NODES Compute the marginal on the specified nodes (hmm_2TBN)
% marginal = marginal_nodes(engine, b, nodes, t, add_ev)
%
% nodes must be a singleton set

assert(length(nodes)==1)
ss = engine.slice_size;

i = nodes(1);
bigT = b.gamma;
dom = i + (t-1)*ss;

%id = engine.marg_singleton_ndx_id(i);
%global SD_NDX
%ndx = SD_NDX{id};
%marginal.T = marg_table_ndxSD(bigT, engine.maximize, ndx);

ns = engine.eff_node_sizes(:);
bigdom = 1:ss;
marginal.T = marg_table(bigT, bigdom + (t-1)*ss, ns(bigdom), dom, engine.maximize);

marginal.domain = dom;
assert(~add_ev);
%if add_ev
%  marginal = add_ev_to_dmarginal(marginal, engine.evidence, engine.node_sizes);
%end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function out=mpExpForm(xval,xexp)

out=[xval,'e',sprintf('%i',xexp)];

author: opencobra

repo: cobratoolbox


text:
function pass = test_linearInit(pref)
% A linear CHEBOP test. This test tests that if we pass an initial guess to a
% linear CHEBOP N that the solution is actually correct.

%% Setup
dom = [0 pi];
if ( nargin == 0 )
    pref = cheboppref;
end

%% Simple scalar problem
N = chebop(@(x,u) diff(u,2) + x.*u, dom);
N.lbc = 2;
N.rbc = 3;

x = chebfun(@(x) x, dom);
rhs = sin(x);

%% Simple
% Start with chebcolloc2
pref.discretization = @chebcolloc2;
N.init = sin(20*x);
u1 = solvebvp(N, rhs, pref);
err(1) = norm(N(u1) - rhs);
err(2) = abs(feval(N.lbc(u1), dom(1))) + abs(feval(N.rbc(u1), dom(end)));
%% Add a breakpoint, change discretization
pref.discretization = @ultraS;
N.domain = [0 1 pi];
u2 = solvebvp(N, rhs, pref);
err(3) = norm(N(u2) - rhs);
err(4) = abs(feval(N.lbc(u1), dom(1))) + abs(feval(N.rbc(u1), dom(end)));
%% System:
d = [-pi pi];
A = chebop(@(x,u,v) [u - diff(v) ; diff(u) + v],d);
A.lbc = @(u,v) u + 1;
A.rbc = @(u,v) v;
x = chebfun('x',d);
N.init = [cos(20*x); sin(20*x)];
rhs = [exp(-x.^2); 1+exp(-x.^2)];
uv = mldivide(A, rhs, pref);
err(5) = norm(A(uv)-rhs);
err(6) = abs(feval(A.lbc(uv{1}, uv{2}), d(1))) + ...
    abs(feval(A.rbc(uv{1}, uv{2}), d(end)));
%% Happy?
tol = 1e3*pref.bvpTol;
pass = err < tol;


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

%-----------------------------------------------------------------------
% This script demonstrates the use of function QUADPROG on the basis of
% the Markowitz mean-variance portfolio optimization problen

clear
close all
clc
warning 'off'

% load matrix data
% of annualized linear stock returns
load 'LinRetMat'

% M: number of samples
% N: number of assets
[~,N] = size(data);
% sample mean
mu = mean(data)';
% sample covariance
Sigma = cov(data);

% Objective function 0.5*x'Hx + c'x
%---------------------------
H = 2.0*Sigma;
c = zeros(N,1);

% Constraints
%---------------------------
% weights sum to one
Aeq = ones(1,N);
beq = 1.0;
% no short-sales
A = -eye(N);
b = zeros(N,1);

% Minimum Variance Portfolio
%---------------------------
[weightsP,sigmaP2] = quadprog(H,c,A,b,Aeq,beq);
muP = weightsP'*mu;

% maximum return
muMax = max(mu);

% number portfolios
M = 50;
% stepsize of returns
delta = (muMax - muP)/(M-1);

% preallocate cache
weightsMat = zeros(N,M);    % matrix of portfolio weights
sigmaVec = zeros(M,1);      % expected volatility of portfolios
muVec = zeros(M,1);         % expected return of portfolios

% add minimum variance portfolio
weightsMat(:,1) = weightsP;
sigmaVec(1) = sqrt(sigmaP2);
muVec(1) = muP;

% compute efficient portfolios
for i = 2:M
    % set equality constraints w1 + w2 +...+ wN = 1 and w'*mu = R
    Aeq = [ones(1,N); mu'];
    beq = [1.0; muP + (i-1)*delta];
    % start optimization
    [weightsMat(:,i),sig2] = quadprog(H,c,A,b,Aeq,beq);

    sigmaVec(i) = sqrt(sig2);
    muVec(i) = weightsMat(:,i)'*mu;
end

weightsMat(weightsMat < 0.001) = 0.0;

% Create figure
figure1 = figure('Color',[1 1 1]);
colormap('gray');

% plot efficient frontier
subplot(1,2,1)
plot(sigmaVec,muVec,'k','LineWidth',1.5)
legend('Efficient Frontier','Location','NorthWest')
set(gca,'xlim',[0 sigmaVec(end)],'ylim',[0 muVec(end)])
title('Mean-Variance Efficient Frontier','FontSize',17)
xlabel('$\sigma$','Interpreter','latex','FontName','Helvetia','FontSize',16)
ylabel('$\mu$','Interpreter','latex','FontName','Helvetia','FontSize',16,'Rotation',0)
ticks = get(gca,'YTick');
set(gca,'YTickLabel',[num2str(ticks'*100),repmat('%',length(ticks),1)])
ticks = get(gca,'XTick');
set(gca,'XTickLabel',[num2str(ticks'*100),repmat('%',length(ticks),1)])
% plot portfolio weights
subplot(1,2,2)
Data = cumsum(weightsMat,1);
step = (0.8-0.3)/(N-1);
colMat = repmat((0.3:step:0.8)',1,3);
for k = 1:N
    hold on
    fill([sigmaVec(1);sigmaVec;sigmaVec(end)],[0,Data(N-k+1,:),0]',colMat(k,:));
end
set(gca,'xlim',[sigmaVec(1) sigmaVec(end)],'ylim',[0 max(max(Data))])
title('Weights of the Efficient Portfolios','FontSize',17)
xlabel('$\sigma$','Interpreter','latex','FontName','Helvetia','FontSize',16)
ylabel('$w$','Interpreter','latex','FontName','Helvetia','FontSize',16,'Rotation',0)
ticks = get(gca,'YTick');
set(gca,'YTickLabel',[num2str(ticks'*100),repmat('%',length(ticks),1)])
ticks = get(gca,'XTick');
set(gca,'XTickLabel',[num2str(ticks'*100),repmat('%',length(ticks),1)])
legend([repmat('w_{',N,1),num2str((1:N)'),repmat('}',N,1)])


text:
%%%  AttitudeControl_Sim
clear
path('./icon/',path);
Init;

%Constant value
RAD2DEG = 57.2957795;
DEG2RAD = 0.0174533;
%throttle when UAV is hovering
THR_HOVER = 0.609;

%% Initial condition
ModelInit_PosE = [0, 0, 0];
ModelInit_VelB = [0, 0, 0];
ModelInit_AngEuler = [0, 0, 0];
ModelInit_RateB = [0, 0, 0];
ModelInit_Rads = 0;
%% control parameter
%Attitude PID parameters
Kp_PITCH_ANGLE = 6.5;
Kp_PITCH_AngleRate = 0.1;
Ki_PITCH_AngleRate = 0.02;
Kd_PITCH_AngleRate = 0.001;
Kp_ROLL_ANGLE = 6.5;
Kp_ROLL_AngleRate = 0.1;
Ki_ROLL_AngleRate = 0.02;
Kd_ROLL_AngleRate = 0.001;

Kp_YAW_AngleRate = 0.5;
Ki_YAW_AngleRate = 0.01;
Kd_YAW_AngleRate = 0.00;
%integral saturation
Saturation_I_RP_Max = 0.3;
Saturation_I_RP_Min = -0.3;
Saturation_I_Y_Max = 0.2;
Saturation_I_Y_Min = -0.2;
%max control angle,default 35deg
MAX_CONTROL_ANGLE_ROLL = 35;
MAX_CONTROL_ANGLE_PITCH  = 35;
%max control angle rate,rad/s
MAX_CONTROL_ANGLE_RATE_PITCH = 220;
MAX_CONTROL_ANGLE_RATE_ROLL = 220;
MAX_CONTROL_ANGLE_RATE_Y = 200;
%% run simulink model
AttitudeControl_Sim

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
function g = p25_g ( n, x )

%*****************************************************************************80
%
%% P25_G evaluates the gradient for problem 25.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real G(N), the gradient of the objective function.
%
  g = zeros ( n, 1 );

  for i = 1 : n
    g(i) = i * 4.0 * x(i)^3;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function h = rgb2hsv(m)

% RGB2HSV converts red-green-blue colors to hue-saturation-value.
%
% this code is based on the comments in
% http://stackoverflow.com/questions/3018313/algorithm-to-convert-rgb-to-hsv-and-hsv-to-rgb-in-range-0-255-for-both

if isvector(m) && length(m)==3
  r = m(1);
  g = m(2);
  b = m(3);
elseif ismatrix(m) && size(m,2)==3
  r = m(:,1);
  g = m(:,2);
  b = m(:,3);
elseif ndims(m)==3
  r = m(:,:,1);
  g = m(:,:,2);
  b = m(:,:,3);
end

% ensure these to be double precision column vectors
r = double(r(:));
g = double(g(:));
b = double(b(:));

minrgb = min([r g b], [], 2);
maxrgb = max([r g b], [], 2);
delta  = maxrgb - minrgb;

h = nan(size(delta)); % will be determined further down
v = maxrgb;
s = delta./maxrgb;

% set for each element in the array the proper hue, depending on whether red, green or blue is the largest value
% the 0, 2 and 4 are there to end up in the right corner of the color circle
sel = (r==maxrgb);
hue = 0.0 + (g-b)./delta; % between yellow & magenta
h(sel) = hue(sel);
sel = (g==maxrgb);
hue = 2.0 + (b-r)./delta; % between cyan & yellow
h(sel) = hue(sel);
sel = (b==maxrgb);
hue = 4.0 + (r-g)./delta; % between magenta & cyan
h(sel) = hue(sel);

% these should not be NaN
h(delta==0) = 0;
s(delta==0) = 0;

% convert to degrees between 0 and 360
h = h * 60;
h = mod(h, 360);

% convert to a fraction between 0 and 1
h = h/360;

if nargout==3
  % return the three output arguments separately
else
  % combine the hue-saturation-value in a single array
  h = reshape([h(:) s(:) v(:)], size(m));
  clear s v
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%  internal function

%  'xform_nii.m' is an internal function called by "load_nii.m", so
%  you do not need run this program by yourself. It does simplified
%  NIfTI sform/qform affine transform, and supports some of the
%  affine transforms, including translation, reflection, and
%  orthogonal rotation (N*90 degree).
%
%  For other affine transforms, e.g. any degree rotation, shearing
%  etc. you will have to use the included 'reslice_nii.m' program
%  to reslice the image volume. 'reslice_nii.m' is not called by
%  any other program, and you have to run 'reslice_nii.m' explicitly
%  for those NIfTI files that you want to reslice them.
%
%  Since 'xform_nii.m' does not involve any interpolation or any
%  slice change, the original image volume is supposed to be
%  untouched, although it is translated, reflected, or even
%  orthogonally rotated, based on the affine matrix in the
%  NIfTI header.
%
%  However, the affine matrix in the header of a lot NIfTI files
%  contain slightly non-orthogonal rotation. Therefore, optional
%  input parameter 'tolerance' is used to allow some distortion
%  in the loaded image for any non-orthogonal rotation or shearing
%  of NIfTI affine matrix. If you set 'tolerance' to 0, it means
%  that you do not allow any distortion. If you set 'tolerance' to
%  1, it means that you do not care any distortion. The image will
%  fail to be loaded if it can not be tolerated. The tolerance will
%  be set to 0.1 (10%), if it is default or empty.
%
%  Because 'reslice_nii.m' has to perform 3D interpolation, it can
%  be slow depending on image size and affine matrix in the header.
%
%  After you perform the affine transform, the 'nii' structure
%  generated from 'xform_nii.m' or new NIfTI file created from
%  'reslice_nii.m' will be in RAS orientation, i.e. X axis from
%  Left to Right, Y axis from Posterior to Anterior, and Z axis
%  from Inferior to Superior.
%
%  NOTE: This function should be called immediately after load_nii.
%
%  Usage: [ nii, idem? ] = xform_nii(nii, [tolerance], [preferredForm])
%
%  nii	- NIFTI structure (returned from load_nii)
%
%  tolerance (optional) - distortion allowed for non-orthogonal rotation
%	or shearing in NIfTI affine matrix. It will be set to 0.1 (10%),
%	if it is default or empty.
%
%  preferredForm (optional)  -  selects which transformation from voxels
%	to RAS coordinates; values are s,q,S,Q.  Lower case s,q indicate
%	"prefer sform or qform, but use others if preferred not present".
%	Upper case indicate the program is forced to use the specificied
%	transform or fail loading.  'preferredForm' will be 's', if it is
%	default or empty.	- Jeff Gunter
%
%  NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function [nii, idem] = xform_nii(nii, tolerance, preferredForm)

%  save a copy of the header as it was loaded.  This is the
%  header before any sform, qform manipulation is done.
%
nii.original.hdr = nii.hdr;

if ~exist('tolerance','var') | isempty(tolerance)
    tolerance = 0.1;
elseif(tolerance<=0)
    tolerance = eps;
end

if ~exist('preferredForm','var') | isempty(preferredForm)
    preferredForm= 's';				% Jeff
end

%  if scl_slope field is nonzero, then each voxel value in the
%  dataset should be scaled as: y = scl_slope * x + scl_inter
%  I bring it here because hdr will be modified by change_hdr.
%
if nii.hdr.dime.scl_slope ~= 0 && ...
        ismember(nii.hdr.dime.datatype, [2,4,8,16,64,256,512,768]) && ...
        (nii.hdr.dime.scl_slope ~= 1 | nii.hdr.dime.scl_inter ~= 0)

    nii.img = ...
        nii.hdr.dime.scl_slope * double(nii.img) + nii.hdr.dime.scl_inter;

    if nii.hdr.dime.datatype == 64

        nii.hdr.dime.datatype = 64;
        nii.hdr.dime.bitpix = 64;
    else
        nii.img = single(nii.img);

        nii.hdr.dime.datatype = 16;
        nii.hdr.dime.bitpix = 32;
    end

    nii.hdr.dime.glmax = max(double(nii.img(:)));
    nii.hdr.dime.glmin = min(double(nii.img(:)));

    %  set scale to non-use, because it is applied in xform_nii
    %
    nii.hdr.dime.scl_slope = 0;

end

%  However, the scaling is to be ignored if datatype is DT_RGB24.

%  If datatype is a complex type, then the scaling is to be applied
%  to both the real and imaginary parts.
%
if nii.hdr.dime.scl_slope ~= 0 && ...
        ismember(nii.hdr.dime.datatype, [32,1792])

    nii.img = ...
        nii.hdr.dime.scl_slope * double(nii.img) + nii.hdr.dime.scl_inter;

    if nii.hdr.dime.datatype == 32
        nii.img = single(nii.img);
    end

    nii.hdr.dime.glmax = max(double(nii.img(:)));
    nii.hdr.dime.glmin = min(double(nii.img(:)));

    %  set scale to non-use, because it is applied in xform_nii
    %
    nii.hdr.dime.scl_slope = 0;

end

%  There is no need for this program to transform Analyze data
%
if nii.filetype == 0 && exist([nii.fileprefix '.mat'],'file')
    load([nii.fileprefix '.mat']);	% old SPM affine matrix
    R=M(1:3,1:3);
    T=M(1:3,4);
    T=R*ones(3,1)+T;
    M(1:3,4)=T;
    nii.hdr.hist.qform_code=0;
    nii.hdr.hist.sform_code=1;
    nii.hdr.hist.srow_x=M(1,:);
    nii.hdr.hist.srow_y=M(2,:);
    nii.hdr.hist.srow_z=M(3,:);
elseif nii.filetype == 0
    nii.hdr.hist.rot_orient = [];
    nii.hdr.hist.flip_orient = [];
    return;				% no sform/qform for Analyze format
end

hdr = nii.hdr;

[hdr,orient]=change_hdr(hdr,tolerance,preferredForm);

%  flip and/or rotate image data
%
if ~isequal(orient, [1 2 3])
    idem=false;
    old_dim = hdr.dime.dim([2:4]);

    %  More than 1 time frame
    %
    if ndims(nii.img) > 3
        pattern = 1:prod(old_dim);
    else
        pattern = [];
    end

    if ~isempty(pattern)
        pattern = reshape(pattern, old_dim);
    end

    %  calculate for rotation after flip
    %
    rot_orient = mod(orient + 2, 3) + 1;

    %  do flip:
    %
    flip_orient = orient - rot_orient;

    for i = 1:3
        if flip_orient(i)
            if ~isempty(pattern)
                pattern = flipdim(pattern, i);
            else
                nii.img = flipdim(nii.img, i);
            end
        end
    end

    %  get index of orient (rotate inversely)
    %
    [tmp rot_orient] = sort(rot_orient);

    new_dim = old_dim;
    new_dim = new_dim(rot_orient);
    hdr.dime.dim([2:4]) = new_dim;

    new_pixdim = hdr.dime.pixdim([2:4]);
    new_pixdim = new_pixdim(rot_orient);
    hdr.dime.pixdim([2:4]) = new_pixdim;

    %  re-calculate originator
    %
    tmp = hdr.hist.originator([1:3]);
    tmp = tmp(rot_orient);
    flip_orient = flip_orient(rot_orient);

    for i = 1:3
        if flip_orient(i) && ~isequal(tmp(i), 0)
            tmp(i) = new_dim(i) - tmp(i) + 1;
        end
    end

    hdr.hist.originator([1:3]) = tmp;
    hdr.hist.rot_orient = rot_orient;
    hdr.hist.flip_orient = flip_orient;

    %  do rotation:
    %
    if ~isempty(pattern)
        pattern = permute(pattern, rot_orient);
        pattern = pattern(:);

        if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792 | ...
                hdr.dime.datatype == 128 | hdr.dime.datatype == 511

            tmp = reshape(nii.img(:,:,:,1), [prod(new_dim) hdr.dime.dim(5:8)]);
            tmp = tmp(pattern, :);
            nii.img(:,:,:,1) = reshape(tmp, [new_dim       hdr.dime.dim(5:8)]);

            tmp = reshape(nii.img(:,:,:,2), [prod(new_dim) hdr.dime.dim(5:8)]);
            tmp = tmp(pattern, :);
            nii.img(:,:,:,2) = reshape(tmp, [new_dim       hdr.dime.dim(5:8)]);

            if hdr.dime.datatype == 128 | hdr.dime.datatype == 511
                tmp = reshape(nii.img(:,:,:,3), [prod(new_dim) hdr.dime.dim(5:8)]);
                tmp = tmp(pattern, :);
                nii.img(:,:,:,3) = reshape(tmp, [new_dim       hdr.dime.dim(5:8)]);
            end

        else
            nii.img = reshape(nii.img, [prod(new_dim) hdr.dime.dim(5:8)]);
            nii.img = nii.img(pattern, :);
            nii.img = reshape(nii.img, [new_dim       hdr.dime.dim(5:8)]);
        end
    else
        if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792 | ...
                hdr.dime.datatype == 128 | hdr.dime.datatype == 511

            nii.img(:,:,:,1) = permute(nii.img(:,:,:,1), rot_orient);
            nii.img(:,:,:,2) = permute(nii.img(:,:,:,2), rot_orient);

            if hdr.dime.datatype == 128 | hdr.dime.datatype == 511
                nii.img(:,:,:,3) = permute(nii.img(:,:,:,3), rot_orient);
            end
        else
            nii.img = permute(nii.img, rot_orient);
        end
    end
else
    idem=true;
    hdr.hist.rot_orient = [];
    hdr.hist.flip_orient = [];

end

nii.hdr = hdr;

return;					% xform_nii

%-----------------------------------------------------------------------
function [hdr, orient] = change_hdr(hdr, tolerance, preferredForm)

orient = [1 2 3];
affine_transform = 1;

%  NIFTI can have both sform and qform transform. This program
%  will check sform_code prior to qform_code by default.
%
%  If user specifys "preferredForm", user can then choose the
%  priority.					- Jeff
%
useForm=[];					% Jeff

if isequal(preferredForm,'S')
    if isequal(hdr.hist.sform_code,0)
        error('User requires sform, sform not set in header');
    else
        useForm='s';
    end
end						% Jeff

if isequal(preferredForm,'Q')
    if isequal(hdr.hist.qform_code,0)
        error('User requires qform, qform not set in header');
    else
        useForm='q';
    end
end						% Jeff

if isequal(preferredForm,'s')
    if hdr.hist.sform_code > 0
        useForm='s';
    elseif hdr.hist.qform_code > 0
        useForm='q';
    end
end						% Jeff

if isequal(preferredForm,'q')
    if hdr.hist.qform_code > 0
        useForm='q';
    elseif hdr.hist.sform_code > 0
        useForm='s';
    end
end						% Jeff

if isequal(useForm,'s')
    R = [hdr.hist.srow_x(1:3)
        hdr.hist.srow_y(1:3)
        hdr.hist.srow_z(1:3)];

    T = [hdr.hist.srow_x(4)
        hdr.hist.srow_y(4)
        hdr.hist.srow_z(4)];

    if det(R) == 0 || ~isequal(R(find(R)), sum(R)')
        hdr.hist.old_affine = [ [R;[0 0 0]] [T;1] ];
        R_sort = sort(abs(R(:)));
        if tolerance ==1
            R(2:3,1) = 0; R([1,3],2) = 0; R(1:2,3) = 0;
        else
            R( find( abs(R) < tolerance*min(R_sort(end-2:end)) ) ) = 0;
        end
        hdr.hist.new_affine = [ [R;[0 0 0]] [T;1] ];

        if det(R) == 0 || ~isequal(R(find(R)), sum(R)')
            msg = [char(10) char(10) '   Non-orthogonal rotation or shearing '];
            msg = [msg 'found inside the affine matrix' char(10)];
            msg = [msg '   in this NIfTI file. You have 3 options:' char(10) char(10)];
            msg = [msg '   1. Using included ''reslice_nii.m'' program to reslice the NIfTI' char(10)];
            msg = [msg '      file. I strongly recommand this, because it will not cause' char(10)];
            msg = [msg '      negative effect, as long as you remember not to do slice' char(10)];
            msg = [msg '      time correction after using ''reslice_nii.m''.' char(10) char(10)];
            msg = [msg '   2. Using included ''load_untouch_nii.m'' program to load image' char(10)];
            msg = [msg '      without applying any affine geometric transformation or' char(10)];
            msg = [msg '      voxel intensity scaling. This is only for people who want' char(10)];
            msg = [msg '      to do some image processing regardless of image orientation' char(10)];
            msg = [msg '      and to save data back with the same NIfTI header.' char(10) char(10)];
            msg = [msg '   3. Increasing the tolerance to allow more distortion in loaded' char(10)];
            msg = [msg '      image, but I don''t suggest this.' char(10) char(10)];
            msg = [msg '   To get help, please type:' char(10) char(10) '   help reslice_nii.m' char(10)];
            msg = [msg '   help load_untouch_nii.m' char(10) '   help load_nii.m'];
            error(msg);
        end
    end

elseif isequal(useForm,'q')
    b = hdr.hist.quatern_b;
    c = hdr.hist.quatern_c;
    d = hdr.hist.quatern_d;

    if 1.0-(b*b+c*c+d*d) < 0
        if abs(1.0-(b*b+c*c+d*d)) < 1e-5
            a = 0;
        else
            error('Incorrect quaternion values in this NIFTI data.');
        end
    else
        a = sqrt(1.0-(b*b+c*c+d*d));
    end

    qfac = hdr.dime.pixdim(1);
    if qfac==0, qfac = 1; end
    i = hdr.dime.pixdim(2);
    j = hdr.dime.pixdim(3);
    k = qfac * hdr.dime.pixdim(4);

    R = [a*a+b*b-c*c-d*d     2*b*c-2*a*d        2*b*d+2*a*c
        2*b*c+2*a*d         a*a+c*c-b*b-d*d    2*c*d-2*a*b
        2*b*d-2*a*c         2*c*d+2*a*b        a*a+d*d-c*c-b*b];

    T = [hdr.hist.qoffset_x
        hdr.hist.qoffset_y
        hdr.hist.qoffset_z];

    %  qforms are expected to generate rotation matrices R which are
    %  det(R) = 1; we'll make sure that happens.
    %
    %  now we make the same checks as were done above for sform data
    %  BUT we do it on a transform that is in terms of voxels not mm;
    %  after we figure out the angles and squash them to closest
    %  rectilinear direction. After that, the voxel sizes are then
    %  added.
    %
    %  This part is modified by Jeff Gunter.
    %
    if det(R) == 0 | ~isequal(R(find(R)), sum(R)')

        %  det(R) == 0 is not a common trigger for this ---
        %  R(find(R)) is a list of non-zero elements in R; if that
        %  is straight (not oblique) then it should be the same as
        %  columnwise summation. Could just as well have checked the
        %  lengths of R(find(R)) and sum(R)' (which should be 3)
        %
        hdr.hist.old_affine = [ [R * diag([i j k]);[0 0 0]] [T;1] ];
        R_sort = sort(abs(R(:)));
        if tolerance ==1
            R(2:3,1) = 0; R([1,3],2) = 0; R(1:2,3) = 0;
        else
            R( find( abs(R) < tolerance*min(R_sort(end-2:end)) ) ) = 0;
        end
        R = R * diag([i j k]);
        hdr.hist.new_affine = [ [R;[0 0 0]] [T;1] ];

        if det(R) == 0 | ~isequal(R(find(R)), sum(R)')
            msg = [char(10) char(10) '   Non-orthogonal rotation or shearing '];
            msg = [msg 'found inside the affine matrix' char(10)];
            msg = [msg '   in this NIfTI file. You have 3 options:' char(10) char(10)];
            msg = [msg '   1. Using included ''reslice_nii.m'' program to reslice the NIfTI' char(10)];
            msg = [msg '      file. I strongly recommand this, because it will not cause' char(10)];
            msg = [msg '      negative effect, as long as you remember not to do slice' char(10)];
            msg = [msg '      time correction after using ''reslice_nii.m''.' char(10) char(10)];
            msg = [msg '   2. Using included ''load_untouch_nii.m'' program to load image' char(10)];
            msg = [msg '      without applying any affine geometric transformation or' char(10)];
            msg = [msg '      voxel intensity scaling. This is only for people who want' char(10)];
            msg = [msg '      to do some image processing regardless of image orientation' char(10)];
            msg = [msg '      and to save data back with the same NIfTI header.' char(10) char(10)];
            msg = [msg '   3. Increasing the tolerance to allow more distortion in loaded' char(10)];
            msg = [msg '      image, but I don''t suggest this.' char(10) char(10)];
            msg = [msg '   To get help, please type:' char(10) char(10) '   help reslice_nii.m' char(10)];
            msg = [msg '   help load_untouch_nii.m' char(10) '   help load_nii.m'];
            error(msg);
        end

    else
        R = R * diag([i j k]);
    end					% 1st det(R)

else
    affine_transform = 0;	% no sform or qform transform
end

if affine_transform == 1
    voxel_size = abs(sum(R,1));
    inv_R = inv(R);
    originator = inv_R*(-T)+1;
    orient = get_orient(inv_R);

    %  modify pixdim and originator
    %
    %hdr.dime.pixdim(2:4) = voxel_size;
    hdr.hist.originator(1:3) = originator;

    %  set sform or qform to non-use, because they have been
    %  applied in xform_nii
    %
    hdr.hist.qform_code = 0;
    hdr.hist.sform_code = 0;
end

%  apply space_unit to pixdim if not 1 (mm)
%
space_unit = get_units(hdr);

if space_unit ~= 1
    hdr.dime.pixdim(2:4) = hdr.dime.pixdim(2:4) * space_unit;

    %  set space_unit of xyzt_units to millimeter, because
    %  voxel_size has been re-scaled
    %
    hdr.dime.xyzt_units = char(bitset(hdr.dime.xyzt_units,1,0));
    hdr.dime.xyzt_units = char(bitset(hdr.dime.xyzt_units,2,1));
    hdr.dime.xyzt_units = char(bitset(hdr.dime.xyzt_units,3,0));
end

hdr.dime.pixdim = abs(hdr.dime.pixdim);

return;					% change_hdr

%-----------------------------------------------------------------------
function orient = get_orient(R)

orient = [];

for i = 1:3
    switch find(R(i,:)) * sign(sum(R(i,:)))
        case 1
            orient = [orient 1];		% Left to Right
        case 2
            orient = [orient 2];		% Posterior to Anterior
        case 3
            orient = [orient 3];		% Inferior to Superior
        case -1
            orient = [orient 4];		% Right to Left
        case -2
            orient = [orient 5];		% Anterior to Posterior
        case -3
            orient = [orient 6];		% Superior to Inferior
    end
end

return;					% get_orient

%-----------------------------------------------------------------------
function [space_unit, time_unit] = get_units(hdr)

switch bitand(hdr.dime.xyzt_units, 7)	% mask with 0x07
    case 1
        space_unit = 1e+3;		% meter, m
    case 3
        space_unit = 1e-3;		% micrometer, um
    otherwise
        space_unit = 1;			% millimeter, mm
end

switch bitand(hdr.dime.xyzt_units, 56)	% mask with 0x38
    case 16
        time_unit = 1e-3;			% millisecond, ms
    case 24
        time_unit = 1e-6;			% microsecond, us
    otherwise
        time_unit = 1;			% second, s
end

return;					% get_units

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function surfPoints = getSurfacePoints(mask3M)
%"getSurfacePoints"
%   Return the [row,col,slice] of surface points in mask3M, where a surface
%   point is defined as any voxel in the structure that is adjacent to a
%   voxel not in the structure.  Two voxels are considered adjacent when
%   their faces touch, and NOT adjacent if only their corners or edges
%   touch.
%
%   To create the 3D surface analog of mask3M use:
%
%   surfPoints = getSurfacePoints(mask3M);
%   surf3M = repmat(logical(0), size(mask3M));
%   for i=1:size(surfPoints,1)
%        surf3M(surfPoints(i,1),surfPoints(i,2), surfPoints(i,3)) = 1;
%   end
%
%   Constructing a 2nd 3D mask the size of mask3M may pose a memory
%   problem, so substituting mask3M for surf3M is a good idea if
%   mask3M is disposable.
%
% JRA 11/20/03
%     03/23/05 - New algorithm for speed, ~10x faster
%
% Usage: surfPoints = getSurfacePoints(mask3M)

surfPoints = [];

[r,c,s] = find3d(mask3M);

%Find minimum rows, cols, slices.
minR = min(r); maxR = max(r);
minC = min(c); maxC = max(c);
minS = min(s); maxS = max(s);

%Restrict surface calculation to region with structure.
mask3M = mask3M(minR:maxR, minC:maxC, minS:maxS);

%Construct the "allNeighborsOn" matrix, which for any voxel not on the edge
%of maskM, tells whether or not ALL of its neighbors are on.
plusRowShift  = mask3M(3:end, 2:end-1, 2:end-1);
allNeighborsOn = plusRowShift;
clear plusRowShift

minusRowShift = mask3M(1:end-2, 2:end-1, 2:end-1);
allNeighborsOn = allNeighborsOn & minusRowShift;
clear minusRowShift

plusColShift  = mask3M(2:end-1, 3:end, 2:end-1);
allNeighborsOn = allNeighborsOn & plusColShift;
clear plusColShift

minusColShift = mask3M(2:end-1, 1:end-2, 2:end-1);
allNeighborsOn = allNeighborsOn & minusColShift;
clear minusColShift

plusSlcShift  = mask3M(2:end-1, 2:end-1, 3:end);
allNeighborsOn = allNeighborsOn & plusSlcShift;
clear plusSlcShift

minusSlcShift = mask3M(2:end-1, 2:end-1, 1:end-2);
allNeighborsOn = allNeighborsOn & minusSlcShift;
clear minusSlcShift

%Now find all surface points (except those on the edge of maskM), defined
%as those points that are ON in mask3M and don't have ALL their neighbors
%on.
kernal = mask3M(2:end-1, 2:end-1, 2:end-1) & ~allNeighborsOn;

%Finally drop the kernal back into the middle of mask3M.  All points on the
%first/last row, column and slice of mask3M are by definion surface points.
mask3M(2:end-1, 2:end-1, 2:end-1) = kernal;

%Find the location of the surface points.
[r,c,s] = find3d(mask3M);

%Correct for taking a subset of mask3M when we began.
r = r + (minR - 1);
c = c + (minC - 1);
s = s + (minS - 1);

surfPoints = [r;c;s]';

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% OtTableNew Connect to an Omnitrans table.
%
% Usage:
%    htable = OtTableNew()
% or:
%    htable = OtTableNew(<path to .db file>)
%    This form includes a call to OtTableSetName()
%
% To release a table handle, call OtTableFree().

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 2006-02-28 Created (Johan Meijdam, VORtech)

function htable = OtTableNew(dbpath)
   if nargin == 1
      htable = Ot('table', dbpath);
   else
      htable = Ot('table');
   end
end


text:
function [inputs, labels, set] = patches_generation(sigma,size_input,size_label,stride,folder,mode,max_numPatches,batchSize)

inputs  = zeros(size_input, size_input, 1, 1,'single');
labels  = zeros(size_label, size_label, 1, 1,'single');
count   = 0;
padding = abs(size_input - size_label)/2;

ext               =  {'*.jpg','*.png','*.bmp'};
filepaths           =  [];

for i = 1 : length(ext)
    filepaths = cat(1,filepaths, dir(fullfile(folder, ext{i})));
end

for i = 1 : length(filepaths)
    image = imread(fullfile(folder,filepaths(i).name)); % uint8
    %[~, name, exte] = fileparts(filepaths(i).name);
    if size(image,3) == 3
        image = rgb2gray(image); % uint8
    end

    for j = 1:8
        image_aug = data_augmentation(image, j);  % augment data
        im_label  = im2single(image_aug); % single
        [hei,wid] = size(im_label);
        im_input  = im_label; % single
        for x = 1 : stride : (hei-size_input+1)
            for y = 1 :stride : (wid-size_input+1)
                subim_input = im_input(x : x+size_input-1, y : y+size_input-1);
                subim_label = im_label(x+padding : x+padding+size_label-1, y+padding : y+padding+size_label-1);
                count       = count+1;
                inputs(:, :, 1, count)   = subim_input + single(sigma/255*randn(size(subim_input)));
                labels(:, :, 1, count) = subim_label;
            end
        end
    end
end

inputs = inputs(:,:,:,1:(size(inputs,4)-mod(size(inputs,4),batchSize)));
labels = labels(:,:,:,1:(size(labels ,4)-mod(size(labels ,4),batchSize)));
labels = shave(inputs,[padding,padding])-labels; %%% residual image patches; pay attention to this!!!

order  = randperm(size(inputs,4));
inputs = inputs(:, :, 1, order);
labels = labels(:, :, 1, order);

set    = uint8(ones(1,size(inputs,4)));
if mode == 1
    set = uint8(2*ones(1,size(inputs,4)));
end

disp('-------Original Datasize-------')
disp(size(inputs,4));

subNum = min(size(inputs,4),max_numPatches);
inputs = inputs(:,:,:,1:subNum);
labels = labels(:,:,:,1:subNum);
set    = set(1:subNum);

disp('-------Now Datasize-------')
disp(size(inputs,4));

author: cszn

repo: DnCNN

save_path: github-repos/MATLAB/cszn-DnCNN


text:
% Library for Quadratic Programming.
%
% Copyright (C) 2006-2008 Vojtech Franc, xfrancv@cmp.felk.cvut.cz
% Center for Machine Perception, CTU FEL Prague
%
%
% libqp_gsmo         Generalized Sequential Minimal Optimizer (SMO) algorithm.
% libqp_splx         Solves Quadratic Programming task with simplex constraints.
%
% libqp_compile      Compiles MEX interface for LIBQP solvers.
%
% libqp_gsmo_test    Test script for LIBQP_GSMO solver.
% libqp_splx_test    Test script for LIBQP_SPLX solver.
%

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function [cost, x, A, info] = test01(n, usestore)
% function [cost, x, A] = test01(n, usestore)
% All intputs are optional.
%
% Typical call:
%
% profile clear; profile on;
% test1(10000, true);
% profile off; profile report;
%
% If activated (search for 'work!' in the code):
% 'work!' is printed each time a matrix-vector product with A is computed.
% Observe how setting 'usestore' to true of false affects the number of
% products.
%

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:

    clc;
    reset(RandStream.getDefaultStream);
    randnfoo = randn(123456, 1); %#ok<NASGU>

    if ~exist('n', 'var') || isempty(n)
        n = 1042;
    end

    if ~exist('usestore', 'var') || isempty(usestore)
        usestore = true;
    end

    % Define the problem data
%     A = randn(n);
    A = magic(n)/n^3;
    A = (A+A')/2;

    % Create the problem structure
    problem.M = spherefactory(n);

    if usestore
        % These functions use the store capability
        problem.cost = @(x, store)    objective(A, x, store);
        problem.grad = @(x, store)    gradient (A, x, store);
        problem.hess = @(x, h, store) hessian  (A, x, h, store);
    else
        % These functions do not use the store capability
%         problem.cost = @(x)    objective(A, x, struct());
%         problem.grad = @(x)    gradient (A, x, struct());
        problem.costgrad = @(x) costgrad(A, x);
        problem.hess = @(x, h) hessian  (A, x, h, struct());
    end

    % Check consistency of cost, grad and hess.
    debug = 0;
    if debug
        checkgradient(problem);
        pause;
        checkhessian(problem);
        pause;
    end

    % Define a few (optional) options
    options.maxtime = 20; % [seconds]
    options.storedepth = 25;
    options.tolgradnorm = 1e-8;
    options.maxinner = 200;

%     options.stopfun = @stopfun;
    function stop = stopfun(problem, x, info, last)
        if mod(last, 150) == 0
            plot([info.time], [info.cost]);
            xlim([0 options.maxtime]);
            drawnow;
        end
        stop = false;
    end

    options.statsfun = @statsfun;
    function stats = statsfun(problem, x, stats)
        stats.pt = problem.M.hash(x);
    end

    % Solve
    x0 = (1:n)'; x0 = x0/norm(x0);
%     options.linesearch = @linesearch;
    options.linesearch = @linesearch_adaptive;
%     options.ls_max_steps = 10;
%     [x cost info] = steepestdescent(problem, x0, options);
    [x cost info] = conjugategradient(problem, x0, options);
    [x cost info] = trustregions(problem, [], options);
%     [x cost info] = pso(problem, [], options);
%     [x cost info] = neldermead(problem, [], options);

%     figure;
%     subplot(1, 2, 1);
%     semilogy([info.iter], [info.gradnorm], '.-');
%     subplot(1, 2, 2);
%     semilogy([info.time], [info.gradnorm], '.-');

    if isfield(info, 'linesearch')
        figure;
        lsstats = [info.linesearch];
        lscostevals = [lsstats.costevals];
        hist(lscostevals, min(lscostevals):max(lscostevals));
        title('Histogram of cost evaluations per line search');
%         keyboard;
    end

%     keyboard;

end

function [val store] = objective(A, x, store)

    if ~isfield(store, 'Ax')
        store.Ax = A*x; % disp('work!');
    end
    Ax = store.Ax;

    if ~isfield(store, 'val')
        store.val = -.5*(x'*Ax);
    end

    val = store.val;

end

function [grad store] = gradient(A, x, store)

    if ~isfield(store, 'Ax') || ~isfield(store, 'val')
        [~, store] = objective(A, x, store);
    end
    Ax = store.Ax;
    val = store.val;

    grad = -(2*val*x + Ax);

end

function [hess store] = hessian(A, x, h, store)

    if ~isfield(store, 'val')
        [~, store] = objective(A, x, store);
    end
    val = store.val;

    Ah = A*h; % disp('work!');
    hess = -(2*val*h + Ah);
    hess = hess - (x'*hess)*x;         % projection

end

function [cost grad] = costgrad(A, x)
    Ax = A*x;
    cost = -.5*(x'*Ax);
    if nargout == 2
        grad = -(2*cost*x + Ax);
    end
end

author: NicolasBoumal

repo: manopt


text:
% THE DEMO PRESENTS AN ERROR PROPAGATION USING A MONTE-CARLO METHOD
% The error in the positioning of the joints is propagated to an error
% in the position of the end effector

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function draw_errors_monte_carlo
close all

M=3000; %number of particles

%load arm parameters
robot=load_robot('ABB', 'IRB140');
%robot=load_robot('example', '2dofplanar');
robot.graphical.draw_transparent=0;

%standard deviation AT EACH JOINT
sigmaq=0.05;%rad

%find errors around this pose
q=[pi/4 pi/2 -pi/8 0 0 0]';
%q=[pi/2 -pi/2]';

puntos=[];
for i=1:M
    %normrnd is included in the Statistics and Machine Learning Toolbox
    %change to the next line if this package is not included in your
    %qi = q + [normrnd(0, sigmaq, robot.DOF, 1)];
    qi = q + [mygaussian(0, sigmaq, robot.DOF, 1)];
    T=directkinematic(robot, qi);
    puntos=[puntos; T(1,4) T(2,4) T(3,4)];
end

adjust_view(robot)
drawrobot3d(robot,q), hold on
plot3(puntos(:,1),puntos(:,2), puntos(:,3),'r.')

function R=mygaussian(mu, sigma, n, m)

global z1
global generate
%close all

%init randomizer
z1 = rand();
generate = 0;

R=[];
for i=1:n
    for j=1:m
        gauss = generate_gaussian(mu, sigma);
        R(i,j)=gauss;
    end
end


text:
%% DEMO of vectorization of elementwise matrix vector multiplication in MATLAB
% D(i,:,:) is a 3x3 or 2x2 matrix on the i-th element
% D can be the tensorial diffusion coefficient of a PDE
% v(i,:) is a 3x1 or 2x1 vector on the i-th element
% v can be the elementwise gradient of a certain nodal basis, e.g.,
% v(i,:) = Dphi(i,:,k) which is the k-th nodal basis's gradient
% We want to compute b(i,:) = D(i,:,:)*v(i,:)
% Even though MATLAB greatly improves the efficiency of for loop after 2016b
% using built-in vectorization operation is still faster.
% Depending on your computer the speed up is about 100x to 150x.
%
% see also Maxwell

%% set up mesh and gradient matrix
[node,elem] = cubemesh([0,1,0,1,0,1],1/32);
center = (node(elem(:,1),:) + node(elem(:,2),:) + ...
    node(elem(:,3),:) + node(elem(:,4),:))/4;
Dlambda = gradbasis3(node,elem);
fprintf('\nSize of gradient of lambda (%d, %d, %d)\n', size(Dlambda));
NT = size(elem,1);

%% set up coefficient matrix (unfortunately arrayfun is not vectorization)
tic;
c = @(p) [1+p(:,1).^2/10 p(:,1).*p(:,2) 0.*p(:,1); ...
    p(:,1).*p(:,2) 1+p(:,2).^2/10 0*p(:,2); ...
    0*p(:,2) 0*p(:,2) 1+p(:,3).^2/10];
tmp = arrayfun(@(rowidx) c(center(rowidx,:)), 1:size(center,1), 'UniformOutput',0);
c2elem = cat(3,tmp{:});
c2elem = permute(c2elem,[3,1,2]);
% permute switch the element idx to the 1st dim
% if the Dlambda tensor has D(:,:,j) corresponds to the j-th element
% then this operation is not neccessary
fprintf('Time to generate elementwise coefficient matrix %5.4g s\n',toc);

%% benchmark the elementwise product
tic;
b = zeros(NT,3,4);
for j = 1:4
    for i = 1:size(elem,1)
        b(i,:,j) = Dlambda(i,:,j)*squeeze(c2elem(i,:,:));
    end
end
t1 = toc;
fprintf('Time to perform elementwise operation %5.4g s\n',t1);

%% benchmark the vectorized routine
tic;
b = zeros(NT,3,4);
for j = 1:4
    b(:,:,j) = sum(bsxfun(@times, c2elem, Dlambda(:,:,j)), 2);
end
t2= toc;
fprintf('Time to perform vectorized operation %5.4g s\n',t2);
fprintf('Speed up factor is %4d.\n', floor(t1/t2));

%% benchmark the vectorized routine
% label for K   1 4 6
%               4 2 5
%               5 6 3

tic;
K(:,4) = center(:,1).*center(:,2);
K(:,3) = 1+center(:,3).^2/10;
K(:,1) = 1+center(:,1).^2/10;
K(:,2) = 1+center(:,2).^2/10;
toc;

tic;
% b1 = K1*D1 + K4*D2 + K6*D3
% b2 = K4*D1 + K2*D2 + K5*D3
% b3 = K5*D1 + K6*D2 + K3*D3
b = zeros(NT,3,4);
for j = 1:4
    b(:,1,j) = K(:,1).*Dlambda(:,1,j) + K(:,4).*Dlambda(:,3,j);
    b(:,2,j) = K(:,4).*Dlambda(:,1,j) + K(:,2).*Dlambda(:,2,j);
    b(:,3,j) = K(:,3).*Dlambda(:,3,j);
end
t3= toc;
fprintf('Time to perform vectorized operation %5.4g s\n',t3);
fprintf('Speed up factor is %4d.\n', floor(t1/t3));

repo: ifem


text:
function [t,f_new,g_new,funEvals,H] = WolfeLineSearch(...
    x,t,d,f,g,gtd,c1,c2,LS,maxLS,tolX,debug,doPlot,saveHessianComp,funObj,varargin)
%
% Bracketing Line Search to Satisfy Wolfe Conditions
%
% Inputs:
%   x: starting location
%   t: initial step size
%   d: descent direction
%   f: function value at starting location
%   g: gradient at starting location
%   gtd: directional derivative at starting location
%   c1: sufficient decrease parameter
%   c2: curvature parameter
%   debug: display debugging information
%   LS: type of interpolation
%   maxLS: maximum number of iterations
%   tolX: minimum allowable step length
%   doPlot: do a graphical display of interpolation
%   funObj: objective function
%   varargin: parameters of objective function
%
% Outputs:
%   t: step length
%   f_new: function value at x+t*d
%   g_new: gradient value at x+t*d
%   funEvals: number function evaluations performed by line search
%   H: Hessian at initial guess (only computed if requested

% Evaluate the Objective and Gradient at the Initial Step
if nargout == 5
    [f_new,g_new,H] = funObj(x + t*d,varargin{:});
else
    [f_new,g_new] = funObj(x+t*d,varargin{:});
end
funEvals = 1;
gtd_new = g_new'*d;

% Bracket an Interval containing a point satisfying the
% Wolfe criteria

LSiter = 0;
t_prev = 0;
f_prev = f;
g_prev = g;
gtd_prev = gtd;
done = 0;

while LSiter < maxLS

    %% Bracketing Phase
    if ~isLegal(f_new) || ~isLegal(g_new)
        if 0
            if debug
                fprintf('Extrapolated into illegal region, Bisecting\n');
            end
            t = (t + t_prev)/2;
            if ~saveHessianComp && nargout == 5
                [f_new,g_new,H] = funObj(x + t*d,varargin{:});
            else
                [f_new,g_new] = funObj(x + t*d,varargin{:});
            end
            funEvals = funEvals + 1;
            gtd_new = g_new'*d;
            LSiter = LSiter+1;
            continue;
        else
            if debug
                fprintf('Extrapolated into illegal region, switching to Armijo line-search\n');
            end
            t = (t + t_prev)/2;
            % Do Armijo
            if nargout == 5
                [t,x_new,f_new,g_new,armijoFunEvals,H] = ArmijoBacktrack(...
                  x,t,d,f,f,g,gtd,c1,max(0,min(LS-2,2)),tolX,debug,doPlot,saveHessianComp,...
                  funObj,varargin{:});
            else
                [t,x_new,f_new,g_new,armijoFunEvals] = ArmijoBacktrack(...
                  x,t,d,f,f,g,gtd,c1,max(0,min(LS-2,2)),tolX,debug,doPlot,saveHessianComp,...
                  funObj,varargin{:});
            end
            funEvals = funEvals + armijoFunEvals;
            return;
        end
    end

    if f_new > f + c1*t*gtd || (LSiter > 1 && f_new >= f_prev)
        bracket = [t_prev t];
        bracketFval = [f_prev f_new];
        bracketGval = [g_prev g_new];
        break;
    elseif abs(gtd_new) <= -c2*gtd
        bracket = t;
        bracketFval = f_new;
        bracketGval = g_new;
        done = 1;
        break;
    elseif gtd_new >= 0
        bracket = [t_prev t];
        bracketFval = [f_prev f_new];
        bracketGval = [g_prev g_new];
        break;
    end
    temp = t_prev;
    t_prev = t;
    minStep = t + 0.01*(t-temp);
    maxStep = t*10;
    if LS == 3
        if debug
            fprintf('Extending Braket\n');
        end
        t = maxStep;
    elseif LS ==4
        if debug
            fprintf('Cubic Extrapolation\n');
        end
        t = polyinterp([temp f_prev gtd_prev; t f_new gtd_new],doPlot,minStep,maxStep);
    else
        t = mixedExtrap(temp,f_prev,gtd_prev,t,f_new,gtd_new,minStep,maxStep,debug,doPlot);
    end

    f_prev = f_new;
    g_prev = g_new;
    gtd_prev = gtd_new;
    if ~saveHessianComp && nargout == 5
        [f_new,g_new,H] = funObj(x + t*d,varargin{:});
    else
        [f_new,g_new] = funObj(x + t*d,varargin{:});
    end
    funEvals = funEvals + 1;
    gtd_new = g_new'*d;
    LSiter = LSiter+1;
end

if LSiter == maxLS
    bracket = [0 t];
    bracketFval = [f f_new];
    bracketGval = [g g_new];
end

%% Zoom Phase

% We now either have a point satisfying the criteria, or a bracket
% surrounding a point satisfying the criteria
% Refine the bracket until we find a point satisfying the criteria
insufProgress = 0;
Tpos = 2;
LOposRemoved = 0;
while ~done && LSiter < maxLS

    % Find High and Low Points in bracket
    [f_LO LOpos] = min(bracketFval);
    HIpos = -LOpos + 3;

    % Compute new trial value
    if LS == 3 || ~isLegal(bracketFval) || ~isLegal(bracketGval)
        if debug
            fprintf('Bisecting\n');
        end
        t = mean(bracket);
    elseif LS == 4
        if debug
            fprintf('Grad-Cubic Interpolation\n');
        end
        t = polyinterp([bracket(1) bracketFval(1) bracketGval(:,1)'*d
            bracket(2) bracketFval(2) bracketGval(:,2)'*d],doPlot);
    else
        % Mixed Case %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        nonTpos = -Tpos+3;
        if LOposRemoved == 0
            oldLOval = bracket(nonTpos);
            oldLOFval = bracketFval(nonTpos);
            oldLOGval = bracketGval(:,nonTpos);
        end
        t = mixedInterp(bracket,bracketFval,bracketGval,d,Tpos,oldLOval,oldLOFval,oldLOGval,debug,doPlot);
    end

    % Test that we are making sufficient progress
    if min(max(bracket)-t,t-min(bracket))/(max(bracket)-min(bracket)) < 0.1
        if debug
            fprintf('Interpolation close to boundary');
        end
        if insufProgress || t>=max(bracket) || t <= min(bracket)
            if debug
                fprintf(', Evaluating at 0.1 away from boundary\n');
            end
            if abs(t-max(bracket)) < abs(t-min(bracket))
                t = max(bracket)-0.1*(max(bracket)-min(bracket));
            else
                t = min(bracket)+0.1*(max(bracket)-min(bracket));
            end
            insufProgress = 0;
        else
            if debug
                fprintf('\n');
            end
            insufProgress = 1;
        end
    else
        insufProgress = 0;
    end

    % Evaluate new point
    if ~saveHessianComp && nargout == 5
        [f_new,g_new,H] = funObj(x + t*d,varargin{:});
    else
        [f_new,g_new] = funObj(x + t*d,varargin{:});
    end
    funEvals = funEvals + 1;
    gtd_new = g_new'*d;
    LSiter = LSiter+1;

    if f_new > f + c1*t*gtd || f_new >= f_LO
        % Armijo condition not satisfied or not lower than lowest
        % point
        bracket(HIpos) = t;
        bracketFval(HIpos) = f_new;
        bracketGval(:,HIpos) = g_new;
        Tpos = HIpos;
    else
        if abs(gtd_new) <= - c2*gtd
            % Wolfe conditions satisfied
            done = 1;
        elseif gtd_new*(bracket(HIpos)-bracket(LOpos)) >= 0
            % Old HI becomes new LO
            bracket(HIpos) = bracket(LOpos);
            bracketFval(HIpos) = bracketFval(LOpos);
            bracketGval(:,HIpos) = bracketGval(:,LOpos);
            if LS == 5
                if debug
                    fprintf('LO Pos is being removed!\n');
                end
                LOposRemoved = 1;
                oldLOval = bracket(LOpos);
                oldLOFval = bracketFval(LOpos);
                oldLOGval = bracketGval(:,LOpos);
            end
        end
        % New point becomes new LO
        bracket(LOpos) = t;
        bracketFval(LOpos) = f_new;
        bracketGval(:,LOpos) = g_new;
        Tpos = LOpos;
    end

    if sum(abs(bracket(1)-bracket(2))*gtd_new) < tolX
        if debug
            fprintf('Line Search can not make further progress\n');
        end
        break;
    end

end

%%
if LSiter == maxLS
    if debug
        fprintf('Line Search Exceeded Maximum Line Search Iterations\n');
    end
end

[f_LO LOpos] = min(bracketFval);
t = bracket(LOpos);
f_new = bracketFval(LOpos);
g_new = bracketGval(:,LOpos);

% Evaluate Hessian at new point
if nargout == 5 && funEvals > 1 && saveHessianComp
    [f_new,g_new,H] = funObj(x + t*d,varargin{:});
    funEvals = funEvals + 1;
end

end

%%
function [t] = mixedExtrap(x0,f0,g0,x1,f1,g1,minStep,maxStep,debug,doPlot);
alpha_c = polyinterp([x0 f0 g0; x1 f1 g1],doPlot,minStep,maxStep);
alpha_s = polyinterp([x0 f0 g0; x1 sqrt(-1) g1],doPlot,minStep,maxStep);
if alpha_c > minStep && abs(alpha_c - x1) < abs(alpha_s - x1)
    if debug
        fprintf('Cubic Extrapolation\n');
    end
    t = alpha_c;
else
    if debug
        fprintf('Secant Extrapolation\n');
    end
    t = alpha_s;
end
end

%%
function [t] = mixedInterp(bracket,bracketFval,bracketGval,d,Tpos,oldLOval,oldLOFval,oldLOGval,debug,doPlot);

% Mixed Case %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nonTpos = -Tpos+3;

gtdT = bracketGval(:,Tpos)'*d;
gtdNonT = bracketGval(:,nonTpos)'*d;
oldLOgtd = oldLOGval'*d;
if bracketFval(Tpos) > oldLOFval
    alpha_c = polyinterp([oldLOval oldLOFval oldLOgtd
        bracket(Tpos) bracketFval(Tpos) gtdT],doPlot);
    alpha_q = polyinterp([oldLOval oldLOFval oldLOgtd
        bracket(Tpos) bracketFval(Tpos) sqrt(-1)],doPlot);
    if abs(alpha_c - oldLOval) < abs(alpha_q - oldLOval)
        if debug
            fprintf('Cubic Interpolation\n');
        end
        t = alpha_c;
    else
        if debug
            fprintf('Mixed Quad/Cubic Interpolation\n');
        end
        t = (alpha_q + alpha_c)/2;
    end
elseif gtdT'*oldLOgtd < 0
    alpha_c = polyinterp([oldLOval oldLOFval oldLOgtd
        bracket(Tpos) bracketFval(Tpos) gtdT],doPlot);
    alpha_s = polyinterp([oldLOval oldLOFval oldLOgtd
        bracket(Tpos) sqrt(-1) gtdT],doPlot);
    if abs(alpha_c - bracket(Tpos)) >= abs(alpha_s - bracket(Tpos))
        if debug
            fprintf('Cubic Interpolation\n');
        end
        t = alpha_c;
    else
        if debug
            fprintf('Quad Interpolation\n');
        end
        t = alpha_s;
    end
elseif abs(gtdT) <= abs(oldLOgtd)
    alpha_c = polyinterp([oldLOval oldLOFval oldLOgtd
        bracket(Tpos) bracketFval(Tpos) gtdT],...
        doPlot,min(bracket),max(bracket));
    alpha_s = polyinterp([oldLOval sqrt(-1) oldLOgtd
        bracket(Tpos) bracketFval(Tpos) gtdT],...
        doPlot,min(bracket),max(bracket));
    if alpha_c > min(bracket) && alpha_c < max(bracket)
        if abs(alpha_c - bracket(Tpos)) < abs(alpha_s - bracket(Tpos))
            if debug
                fprintf('Bounded Cubic Extrapolation\n');
            end
            t = alpha_c;
        else
            if debug
                fprintf('Bounded Secant Extrapolation\n');
            end
            t = alpha_s;
        end
    else
        if debug
            fprintf('Bounded Secant Extrapolation\n');
        end
        t = alpha_s;
    end

    if bracket(Tpos) > oldLOval
        t = min(bracket(Tpos) + 0.66*(bracket(nonTpos) - bracket(Tpos)),t);
    else
        t = max(bracket(Tpos) + 0.66*(bracket(nonTpos) - bracket(Tpos)),t);
    end
else
    t = polyinterp([bracket(nonTpos) bracketFval(nonTpos) gtdNonT
        bracket(Tpos) bracketFval(Tpos) gtdT],doPlot);
end
end

%%
function [legal] = isLegal(v)
legal = sum(any(imag(v(:))))==0 & sum(isnan(v(:)))==0 & sum(isinf(v(:)))==0;
end