text:
function regressor=task_regressor(TR,T,taskfile)

%taskfile: If using the Custom 3 column format, prepare a plain text file (avoiding editors such as MS Word)
%containing one line for each stimulus associated with the condition to be modelled by this EV.
%Each line must contain three numbers: the onset time (in seconds); the duration (in seconds);
%the relative magnitude of each stimulus (usually set to 1)
% Example:
% taskfile='fear.txt';
% 32.08   18   1
% 74.223  18   1
% 116.365  18  1

%TR: in seconds
%T: number of frames

%Convolve with HRF
Convolve=1; %0=no 1=yes

hrf=spm_hrf(TR);

%Regressor
regressor=zeros(T,1);
t=dlmread(taskfile);
for i=1:size(t,1)
    start=t(i,1);
    finish=start+t(i,2);

    %Convert from seconds to TRs
    start=round(start/TR)+1;
    finish=round(finish/TR)+1;

    regressor(start:finish)=1;
end
if Convolve; regressor=conv(hrf,regressor); end %Convolve with HRF
regressor=regressor(1:T);

author: yetianmed

repo: subcortex

save_path: github-repos/MATLAB/yetianmed-subcortex


text:
function F = cut(varargin)
%CUT Defines a cut constraint
%
% The syntax for CUT is exactly the same as the
% syntax for ordinary constraints.
%
% The difference between a ordinary constraint and
% a cut constraint is that the CUT will not be used
% in the solution of the upper bound problem in a
% global solver, but only in the relaxation for the
% lower problem.

switch nargin
case 0
    F = lmi;
case 1
    F = lmi(varargin{1});
case 2
    F = lmi(varargin{1},varargin{2});
case 3
    F = lmi(varargin{1},varargin{1},varargin{3});
otherwise
end

F = setcutflag(F);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function gauss_seidel_test ( )

%*****************************************************************************80
%
%% GAUSS_SEIDEL_TEST tests the GAUSS_SEIDEL library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2011
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GAUSS_SEIDEL_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the GAUSS_SEIDEL library.\n' );

  gauss_seidel_test01 ( );
  gauss_seidel_test02 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GAUSS_SEIDEL_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Test file for singfun/isfinite.m

function pass = test_isfinite(pref)

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

% The order of the exponents:
a = 0.64;
b = -0.64;
c = 1.28;
d = -1.28;

%%
% Check a few cases.

% fractional root at the left endpoint
data.exponents = [a 0];
data.singType = {'root', 'none'};
f = singfun(@(x) (1 + x).^a.*exp(x), data, pref);
pass(1) = isfinite(f);

% fractional pole at the left endpoint
data.exponents = [d 0];
data.singType = {'sing', 'none'};
f = singfun(@(x) (1 + x).^d.*sin(x), data, pref);
pass(2) = ~isfinite(f);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%%%
%> @brief  computes the homogeneous transformation A of the pose vector v
%> @param  v pose vector
%> @return A homogeneous transformation
%> @author Giorgio Grisetti
%%%
function A = v2t(v)
% V2T vector to homogeneous transformation
c = cos(v(3));
s = sin(v(3));
A = [c, -s, v(1);
     s,  c, v(2);
     0   0  1];
end

repo: graphslam


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Reset:
%   Set the specified variable to zero.
%   Example: robot = Reset(robot, do1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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 [robot] = Reset(robot, variable)

global configuration

robot.variable = 0;

repo: ARTE


text:
function [X, z, mu] = kmeansRnd(d, k, n)
% Generate samples from a Gaussian mixture distribution with common variances (kmeans model).
% Input:
%   d: dimension of data
%   k: number of components
%   n: number of data
% Output:
%   X: d x n data matrix
%   z: 1 x n response variable
%   mu: d x k centers of clusters
% Written by Mo Chen (sth4nth@gmail.com).
alpha = 1;
beta = nthroot(k,d); % k points in volume x^d : x^d=k

X = randn(d,n);
w = dirichletRnd(alpha,ones(1,k)/k);
z = discreteRnd(w,n);
E = full(sparse(z,1:n,1,k,n,n));
mu = randn(d,k)*beta;
X = X+mu*E;

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function model = yalmip2cbc(interfacedata)

F_struc = interfacedata.F_struc;
Q       = interfacedata.Q;
c       = interfacedata.c;
K       = interfacedata.K;
lb      = full(interfacedata.lb);
ub      = full(interfacedata.ub);
x0      = interfacedata.x0;

if isempty(F_struc)
    Aeq = [];
    beq = [];
    A = [];
    b = [];
else
    Aeq = -F_struc(1:1:K.f,2:end);
    beq = F_struc(1:1:K.f,1);
    A =-F_struc(K.f+1:end,2:end);
    b = F_struc(K.f+1:end,1);
end

ivars = repmat('C',length(c),1);
ivars(interfacedata.integer_variables) = 'I';
ivars(interfacedata.binary_variables) = 'B';

% CBC merges equalities and equalities
ru = full([beq;b]);
rl = full([beq;repmat(-inf,length(b),1)]);
A = [Aeq;A];

% SOS currently not supported
sos.type='';
sos.index=[];
sos.weight=[];

model.f = full(c);
model.A = A;
model.rl = rl;
model.ru = ru;
model.lb = lb;
model.ub = ub;
model.xtype = ivars;
model.sos = sos;
model.x0 = [];
model.H = [];


text:
%%                   testinitSpectral.m

% This test file implements the Spectral as well as the Truncated Spectral
% initializer. The code builds a synthetic formulation of the Phase
% Retrieval problem, b = |Ax| and computes an estimate to x. The code
% finally outputs the correlation achieved.

% PAPER TITLE:
%              Phase Retrieval via Wirtinger Flow: Theory and Algorithms.

% ARXIV LINK:
%              https://arxiv.org/abs/1407.1065

% This is a test script for the spectral initializer.

% 1.) Each test script for an initializer starts out by defining the length
% of the unknown signal, n and the number of measurements, m. These
% mesurements can be complex by setting the isComplex flag to be true.

% 2.) We then build the test problem by generating random gaussian
% measurements and using b0 = abs(Ax) where A is the measurement matrix.

% 3.) We run x0 = initSpectral(A,[],b0,n,isTruncated,isScaled) which runs
% the initialier and and recovers the test vector with high correlation.

% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

%% -----------------------------START----------------------------------

clc
clear
close all

% Parameters
n = 500;           % number of unknowns
m = 8*n;          % number of measurements
isComplex = true;  % use complex matrices? or just stick to real?

%  Build the test problem
xt = randn(n,1)+isComplex*randn(n,1)*1i; % true solution
A = randn(m,n)+isComplex*randn(m,n)*1i;  % matrix
b0 = abs(A*xt);                          % data

% Set up Parameters
isTruncated = true;
isScaled = true;

% Invoke the truncated spectral initial method
x0 = initSpectral(A,[],b0,n,isTruncated,isScaled);

% Calculate the correlation between the recovered signal and the true signal
correlation = abs(x0'*xt/norm(x0)/norm(xt));

fprintf('correlation: %f\n', correlation);

author: tomgoldstein

repo: phasepack-matlab


text:
function t_most_w_ds(quiet, solver, verbose)
%T_MOST_W_DS  Test for MOST with dynamical system constraints.

%   MOST
%   Copyright (c) 2015-2020, Power Systems Engineering Research Center (PSERC)
%   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
%   and Ray Zimmerman, PSERC Cornell
%
%   This file is part of MOST.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/most for more info.

if nargin < 3
    verbose = 0;
    if nargin < 2
        solver = '';
        if nargin < 1
            quiet = 0;
        end
    end
end

include_MIPS = 0;   %% set to 1, to attempt even if MIPS is the best solver
                    %% available (takes a LONG time and currently fails)
n_tests = 3;

t_begin(n_tests, quiet);

casefile = 'c118swf';
solnfile = 't_most_w_ds_z';

%% define named indices into data matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
[CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
    CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
    CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
    CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
    CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
    CT_MODCOST_X] = idx_ct;

if have_feature('cplex') || have_feature('gurobi') || ...
        have_feature('mosek') || have_feature('quadprog_ls') || include_MIPS
    mdi = md_init;

    %% choose solver
    if isempty(solver) || strcmp(upper(solver), 'DEFAULT')
        if have_feature('mosek')
            solver = 'MOSEK';
        elseif have_feature('cplex')
            solver = 'CPLEX';
        elseif have_feature('gurobi')
            solver = 'GUROBI';
        elseif have_feature('quadprog_ls')
            solver = 'OT';
        else
            solver = 'MIPS';
        end
    end
    mpopt = mpoption('most.solver', solver, 'verbose', verbose);

    %% set options
    if have_feature('cplex')
        mpopt = mpoption(mpopt, 'cplex.opts.threads', 2);   % set this manually here
    end
    if have_feature('gurobi')
        mpopt = mpoption(mpopt, 'gurobi.method', 2);        %% barrier
        mpopt = mpoption(mpopt, 'gurobi.threads', 4);
        mpopt = mpoption(mpopt, 'gurobi.opts.BarConvTol', 1e-6);        %% 1e-8
        mpopt = mpoption(mpopt, 'gurobi.opts.FeasibilityTol', 1e-4);    %% 1e-6
        mpopt = mpoption(mpopt, 'gurobi.opts.OptimalityTol', 1e-5);     %% 1e-6
    end
    if have_feature('quadprog_ls')
        mpopt = mpoption(mpopt, 'quadprog.TolFun', 1e-13);
    end
    if have_feature('mosek')
        mpopt = mpoption(mpopt, 'mosek.num_threads', 2);
    else
        mpopt = mpoption(mpopt, 'mips.max_it', 500);
        if have_feature('pardiso')
            mpopt = mpoption(mpopt, 'mips.linsolver', 'PARDISO');
        end
    end
    %% use e.g. t_most_w_ds(0, 'MOSEK') instead of uncommenting these lines
    % mpopt = mpoption(mpopt, 'most.solver', 'MOSEK');
    % mpopt = mpoption(mpopt, 'most.solver', 'CPLEX');
    % mpopt = mpoption(mpopt, 'most.solver', 'GUROBI');
    % mpopt = mpoption(mpopt, 'most.solver', 'OT');
    % mpopt = mpoption(mpopt, 'most.solver', 'BPMPD');
    % mpopt = mpoption(mpopt, 'most.solver', 'MIPS');

    mdi.mpc = loadcase(casefile);
    mdi.InitialPg = mdi.mpc.gen(:,PG);
    nt = 24;
    ng = size(mdi.mpc.gen, 1);
    mdi.idx.nt = nt;
    PositiveActiveReservePrice = ones(ng,1);
    PositiveActiveReserveQuantity = 0.25*mdi.mpc.gen(:,PMAX);
    NegativeActiveReservePrice = ones(ng,1);
    NegativeActiveReserveQuantity = PositiveActiveReserveQuantity;
    PositiveActiveDeltaPrice = ones(ng,1);
    NegativeActiveDeltaPrice = ones(ng,1);
    PositiveLoadFollowReservePrice = ones(ng,1);
    PositiveLoadFollowReserveQuantity = 0.5*mdi.mpc.gen(:,PMAX);
    NegativeLoadFollowReservePrice = ones(ng,1);
    NegativeLoadFollowReserveQuantity = PositiveLoadFollowReserveQuantity;
    %mdi.mpc.gen(:,RAMP_10) = 0.20 * mdi.mpc.gen(PMAX);
    %mdi.mpc.gen(:,RAMP_AGC) = 0.20 * mdi.mpc.gen(PMAX);
    %mdi.mpc.gen(:,RAMP_30) = 0.50 * mdi.mpc.gen(PMAX);
    mdi.mpc.gen(:,RAMP_10) = 1.0 * mdi.mpc.gen(:,PMAX);
    mdi.mpc.gen(:,RAMP_AGC) = 1.0 * mdi.mpc.gen(:,PMAX);
    mdi.mpc.gen(:,RAMP_30) = 1.0 * mdi.mpc.gen(:,PMAX);

    mdi.RampWearCostCoeff = 0.05 * ones(ng,1);   % (i, t) note different scheme!
    for t = 2:nt
      mdi.RampWearCostCoeff(:, t) = mdi.RampWearCostCoeff(:, 1);
    end
    mdi.Storage(1).UnitIdx = mdi.mpc.iess;
    ns = length(mdi.Storage.UnitIdx);
    Minstor = zeros(ns,1);
    Maxstor = 200 * ones(ns,1);
    %mdi.Storage.MinStorageLevel    = zeros(ns,1);
    %mdi.Storage.MaxStorageLevel    = 200 * ones(ns,1);
    mdi.Storage.InitialStorage     = 50 * ones(ns,1);
    mdi.Storage.InitialStorageLowerBound = 50*ones(ns,1);
    mdi.Storage.InitialStorageUpperBound = 50*ones(ns,1);
    mdi.Storage.OutEff             = 0.95 * ones(ns,1);
    mdi.Storage.InEff              = 0.9  * ones(ns ,1);
    mdi.Storage.InitialStorageCost         = 35 * ones(ns, 1);
    mdi.Storage.TerminalStoragePrice       = 35 * ones(ns, 1); % applied to psc_tij0, psd_tij0 (non-terminal states)
    mdi.Storage.TerminalChargingPrice0     = 35 * ones(ns, 1); % applied to psc_tijk (contingency terminal states)
    mdi.Storage.TerminalDischargingPrice0  = 35 * ones(ns, 1); % applied to psd_tijk (contingency terminal states)
    mdi.Storage.TerminalChargingPriceK     = 10 * ones(ns, 1); % applied to psc_tij0 (end-of-horizon terminal states)
    mdi.Storage.TerminalDischargingPriceK  = 40 * ones(ns, 1); % applied to psd_tij0 (end-of-horizon terminal states)
    mpopt = mpoption(mpopt, 'most.storage.terminal_target', 0);
    mdi.Storage.ExpectedTerminalStorageAim = mdi.Storage.InitialStorage;  % expected terminal storage if mpopt.most.storage.terminal_target is true
    mdi.Storage.LossFactor         = zeros(ns,1);  % fraction of storage lost in each period
    mdi.Storage.IncludeValueOfTerminalStorage = 1;
    mpopt = mpoption(mpopt, 'most.storage.cyclic', 1);

    for t = 1:nt
      mdi.offer(t).gencost = mdi.mpc.gencost;
      mdi.offer(t).PositiveActiveReservePrice = PositiveActiveReservePrice;
      mdi.offer(t).PositiveActiveReserveQuantity = PositiveActiveReserveQuantity;
      mdi.offer(t).NegativeActiveReservePrice = NegativeActiveReservePrice;
      mdi.offer(t).NegativeActiveReserveQuantity = NegativeActiveReserveQuantity;
      mdi.offer(t).PositiveActiveDeltaPrice = PositiveActiveDeltaPrice;
      mdi.offer(t).NegativeActiveDeltaPrice = NegativeActiveDeltaPrice;
      mdi.offer(t).PositiveLoadFollowReservePrice = PositiveLoadFollowReservePrice;
      mdi.offer(t).PositiveLoadFollowReserveQuantity = PositiveLoadFollowReserveQuantity;
      mdi.offer(t).NegativeLoadFollowReservePrice = NegativeLoadFollowReservePrice;
      mdi.offer(t).NegativeLoadFollowReserveQuantity = NegativeLoadFollowReserveQuantity;
      mdi.Storage.MinStorageLevel(:,t) = Minstor;
      mdi.Storage.MaxStorageLevel(:,t) = Maxstor;
    end
    mdi.Storage.MinStorageLevel(:,nt+1) = Minstor;  % Needed if mpopt.most.storage.cyclic
    mdi.Storage.MaxStorageLevel(:,nt+1) = Maxstor;

    %mdi.Storage.MinStorageLevel(:,4) = [10 ; 10];  % is this enough to create infeasibility?
    %mdi.Storage.MaxStorageLevel(:,4) = [10; 50 ];

    mdi.UC.CommitSched = ones(ng,nt);

    mdi.Delta_T = 1;
    %            0:0  1:00 2:00 3:00 4:00  5:00  6:00  7:00 8:00 9:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00 23:00
    loadprof = [ 0.6  0.6   0.6 0.7  0.75  0.8   0.9   1.1  1.2  1.3   1.4  1.4   1.3   1.4    1.4   1.4   1.4   1.3   1.1   1.1   1.0   0.9   0.8   0.7  ];

    %                label  probability   type      row     column      chg type  newvalue
    partialcontabrow =[ 1       0        CT_TBUS     0        PD         CT_REL ];
    %mdi.tstep(1).OpCondSched(1).tab= [
    %                   1        0        CT_TBUS     0        PD         CT_REL     0.8  ];
    %mdi.tstep(2).OpCondSched(1).tab= [
    %                   1        0        CT_TBUS     0        PD         CT_REL     1.0  ];
    %mdi.tstep(3).OpCondSched(1).tab= [
    %                   1        0        CT_TBUS     0        PD         CT_REL     1.2  ];
    %mdi.tstep(4).OpCondSched(1).tab= [
    %                   1        0        CT_TBUS     0        PD         CT_REL     0.9 ];

    for t = 1:nt
      mdi.tstep(t).OpCondSched(1).tab = [ partialcontabrow   loadprof(t) ];
    end

    for t = 1:nt
       mdi.tstep(t).OpCondSched(2).tab = mdi.tstep(t).OpCondSched(1).tab;
       mdi.tstep(t).OpCondSched(2).tab(1,7) = 1.1*mdi.tstep(t).OpCondSched(2).tab(1,7);
       mdi.tstep(t).OpCondSched(3).tab = mdi.tstep(t).OpCondSched(1).tab;
       mdi.tstep(t).OpCondSched(3).tab(1,7) = 0.9*mdi.tstep(t).OpCondSched(1).tab(1,7);
    end

    contab = [%         1       0.01      CT_TBUS     0        PD         CT_REL     1.05 ;
                       1       0.01      CT_TGEN     2     GEN_STATUS    CT_REP      0    ;
                       2       0.01      CT_TGEN     5     GEN_STATUS    CT_REP      0    ;
                       ];

    for t = 1:nt
      for j = 1:3  % mdi.idx.nj(t)
        mdi.cont(t,j).contab = contab;
      end
    end

    mdi.tstep(1).TransMat = [ 1/3;
                               1/3
                               1/3];
    for t = 2:nt
      mdi.tstep(t).TransMat = 1/3 * ones(3,3);
    end

    ntds = 24;
    mdi.idx.ntds = ntds;
    m1 = 8;
    m2 = 12;
    B = sparse(m1*m2, ng);
    ilist = [ 2 3 4 5   3 4 5 6   3 4 5 7   3 4 5 7   3 4 5 6   4 5 6 7    2 3 4 ];
    jlist = [ 2 2 2 2   3 3 3 3   5 5 5 5   6 6 6 6   7 7 7 7   8 8 8 8    11 11 11 ];
    for i = 1:length(mdi.mpc.icoal)
     B((jlist(i)-1)*m1+ilist(i), mdi.mpc.icoal(i)) = 0.1;
    end
    A = mkdif(m1, m2, 0.5, 0.97, [1.0 0]);
    C = [];
    D = [];
    zmin = zeros(m1*m2, 1);
    zmax = 100*ones(m1*m2, 1);
    ymin = 0;
    ymax = 100;
    for t = 1:ntds
     mdi.dstep(t).A = A;
     mdi.dstep(t).B = B;
     mdi.dstep(t).C = C;
     mdi.dstep(t).D = D;
     mdi.dstep(t).zmin = zmin;
     mdi.dstep(t).zmax = zmax;
     mdi.dstep(t).ymin = ymin;
     mdi.dstep(t).ymax = ymax;
    end
    mdi.z1 = zeros(m1*m2, 1);

    mdo = most(mdi, mpopt);

    s = load(solnfile);

    t = 'success';
    t_is(mdo.results.success, 1, 12, t);

    t = 'objective function value (f)';
    t_is(mdo.QP.f, 1593399.5, -0.5, t);
% 1593399.487 % CPLEX ~25 sec
% 1594112.218 % GUROBI ~183 sec (fails)
% 1593400.863 % MOSEK ~10 sec
% 1593399.481 % OT ~88 sec

    t = 'dynamical system state (Z)';
    t_is(mdo.results.Z, s.Z, 4, t);

    % Z = mdo.results.Z;
    % save t_most_w_ds_z Z
else
    t_skip(3, 'requires MOSEK, CPLEX, Gurobi or quadprog');
end

% YorN = input('Play movie? (y/n) : ', 's');
% if strcmp(upper(YorN(1)), 'Y')
%     domovie;
% end

t_end

function A = mkdif(m1, m2, alpha, r, w)
% A = mkdif(m1, m2, r, w)
%
% computes an A matrix for the difussion equations in an m1 x m2 grid,
% using a difussion speed factor alpha <= 1, a dissipation factor r < 1 and
% a "wind" vector w = [wx wy] that roughly tells where the pollutants go.

% Carlos Murillo.  As naive as can be.

if norm(w) > 1
  w = w / norm(w);
end

n = m1 * m2;
A = sparse(n, n);
north = alpha / 4;
east = north;
south = north;
west = north;
self = (1-alpha);
if w(1) > 0
  west = west + w(1)*self;
  self = (1-w(1)) * self;
elseif w(2) < 0
  east = east + (-w(1))*self;
  self = (1+w(1)) * self;
end
if w(2) > 0
  south = south + w(2)*self;
  self = (1-w(2)) * self;
elseif w(2) < 0
  north = north + (-w(2))*self;
  self = (1+w(2)) * self;
end
tot = self + north + south + east + west;
self = (r/tot) * self;
north = (r/tot) * north;
south = (r/tot) * south;
east = (r/tot) * east;
west = (r/tot) * west;

for i = 1:m1
  for j = 1:m2
    A((j-1)*m1+i, (j-1)*m1+i) = self;
    % North
    if i > 1
      A((j-1)*m1+i, (j-1)*m1+i-1) = north;
    end
    % South
    if i < m1
      A((j-1)*m1+i, (j-1)*m1+i+1) = south;
    end
    % West
    if j > 1
      A((j-1)*m1+i, (j-2)*m1+i) = west;
    end
    % East
    if j < m2
      A((j-1)*m1+i, j*m1+i) = east;
    end
  end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function n = p06_n ( )

%*****************************************************************************80
%
%% P06_N returns the number of variables N for problem 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer N, the number of variables.
%
  n = 10;

  return
end


text:
function prob = RidgeRegressionTest(data, model)

prob = model.beta' * data.feat;
% prob = data.feat' * model.beta;

author: flyers

repo: drone-tracking


text:
function [p,q] = bandwidth(A)
%BANDWIDTH    Upper and lower bandwidth of matrix A
%
%   A   = 0  for  i-j > p  or j-i > q
%    ij
%
%    [p,q] = bandwidth(A)
%

% written  05/21/09     S.M. Rump
%

  [p,q] = bandwidth(reshape(A.t(1,:),A.size));

repo: OpenMAS


text:
function param = bbox_search_cnn(frame, param, tmpl, background_threshold, frameNum, netId)

num_scale = 4;
neg_num = 16;
p = param.est;
ave_len = (p(3) + p(4)) / 2;

% test in batch to get the probability map of both the postive samples and the negative samples
feas = forward_cnn(tmpl, frame, num_scale + neg_num);

% Collect neg samples with high response on the CNN, which are used for finetune the cnn afterwards
neg_sample = [];
for i = num_scale + 1 : num_scale + neg_num
    if (sum(sum(feas(:, :, i))) > background_threshold)
        neg_sample = [neg_sample, tmpl.basis(:, i)];
    end
end

% bestConf = -1e9;
% bestCord = zeros(1,4);

% hierarchical search in the probability map of different scales to find the best bounding box, details of the searching scheme can be found in the paper
for scale = 1 : num_scale
    curfea = feas(:, :, scale);
    i = 1;
    clear x1_ y1_ x2_ y2_;
    % first get a rough estimation of the bounding box by thresholding the probability map
    for threshold = 0.1 : 0.1 : 0.7
        tmp = curfea > threshold;
        if sum(sum(tmp)) < 10
            break;
        end
        [x1_(i),y1_(i),x2_(i),y2_(i)] = find_box_in_map(tmp);
        i = i + 1;
    end
    if ~exist('x1_', 'var')
        continue;
    end
    cur_cord = mean([x1_', y1_', x2_', y2_'], 1);
    fea = feas(:, :, scale);
    x1 = cur_cord(1); x2 = cur_cord(3);
    y1 = cur_cord(2); y2 = cur_cord(4);

    % param.lastUpdate = 1;

    % save the average response of the region, which is then used to determine whether the cnn will be updated
    curconf = mean(mean(fea(x1 : x2, y1 : y2)));

    %% based on the above estimated center location, we then find a more accurate scale of the bounding box
    % calculate the location in the original frame
    scale_ratio = scale * 0.5;
    scale_width = (p(3) + scale_ratio * ave_len ) / 50;
    scale_height = (p(4) + scale_ratio * ave_len ) / 50;
    new_center = [(x1+x2)/2, (y1+y2)/2];
    center_x = (1 + 50)/2;
    center_y = (1 + 50)/2;
    new_p = p;
    new_p(1) = (p(1) + (new_center(1) - center_x ) * scale_width);
    new_p(2) = (p(2) + (new_center(2) - center_y ) * scale_height);

    % search again in the probability map for a more accurate scale.
    h = p(3) + scale_ratio * ave_len;
    w = p(4) + scale_ratio * ave_len;
    oriResizedMap = imresize(fea, [w, h]);

    total_scale = 0;
    total_ratio = 0;
    for th = 1 : 0.05 : 1.3
        resizedMap = max(oriResizedMap, 0) * 2 - th;
        % penalize boundary values, otherwise the estimated scale turns to be bigger
        resizedMap(resizedMap < -th + 0.1) = -2;
        iter_center(1) = new_center(1) * h / 50;
        iter_center(2) = new_center(2) * w / 50;
        forDecision = -1e9;
        curRes = -1e9;
        bestScale = 0;
        bestRatio = 0;
        sumMap = integral(double(resizedMap));
        for i = -0.02 : 0.001 : 0.02
            for j = 0.00 : 0.002 : 0.00
                newH = p(3) * (1 + i);
                newW = p(4) * (1 + i) * (1 + j);
                x1 = round(max(1, iter_center(2) - newW / 2));
                x2 = round(min(w, iter_center(2) + newW / 2));
                y1 = round(max(1, iter_center(1) - newH / 2));
                y2 = round(min(h, iter_center(1) + newH / 2));
                x1 = x1 - 1;
                y1 = y1 - 1;
                inConf = get(sumMap, x2, y2) - get(sumMap,x1, y2) - ...
                    get(sumMap, x2, y1) + get(sumMap, x1, y1);
                fmeasure = inConf * (newH * newW) / p(3) / p(4);

                if fmeasure > curRes
                    curRes = fmeasure;
                    bestScale = i;
                    bestRatio = j;
                    forDecision = inConf / newH / newW;
                end
            end
        end
        if forDecision < 0
            bestScale = 0;
            bestRatio = 0;
        end
        total_scale = total_scale + bestScale;
        total_ratio = total_ratio + bestRatio;
    end
    total_scale = total_scale / 7;
    total_ratio = total_ratio / 7;
    new_p(3) = ( p(3) * (1 + total_scale) );
    new_p(4) = ( p(4) * (1 + total_scale) * (1 + total_ratio));

    best_p = new_p;
    param.conf = curconf;
    param.neg_sample = neg_sample;
    param.est = best_p;
    return;
end

% if the code runs here, it means we cannot find the target in any scale, thus the target is reported missing.
param.neg_sample = [];
param.conf = 0;
end

function res = get(map, x, y)
    if x == 0 || y == 0
        res = 0;
    else
        res = map(x, y);
    end
end

function [x1,y1,x2,y2] = find_box_in_map(tmp)

x1 = find( sum(tmp, 1) > 0, 1 );
y1 = find( sum(tmp, 2) > 0, 1 );
x2 = size(tmp, 2) + 1 - find( fliplr(sum(tmp, 1)) > 0, 1 );
y2 = size(tmp, 1) + 1 - find( flipud(sum(tmp, 2)) > 0, 1 );

end


function figure2dot(filename, varargin)
%FIGURE2DOT    Save figure in Graphviz (.dot) file.
%  FIGURE2DOT(filename) saves the current figure as dot-file.
%
%  FIGURE2DOT(filename, 'object', HGOBJECT) constructs the graph representation
%  of the specified object (default: gcf)
%
%  You can visualize the constructed DOT file using:
%    - [GraphViz](http://www.graphviz.org) on your computer
%    - [WebGraphViz](http://www.webgraphviz.com) online
%    - [Gravizo](http://www.gravizo.com) for your markdown files
%    - and a lot of other software such as OmniGraffle
%
%  See also: matlab2tikz, cleanfigure, uiinspect, inspect

    ipp = m2tInputParser();
    ipp = ipp.addRequired(ipp, 'filename', @ischar);
    ipp = ipp.addParamValue(ipp, 'object', gcf, @ishghandle);
    ipp = ipp.parse(ipp, filename, varargin{:});
    args = ipp.Results;

    filehandle = fopen(args.filename, 'w');
    finally_fclose_filehandle = onCleanup(@() fclose(filehandle));

    % start printing
    fprintf(filehandle, 'digraph simple_hierarchy {\n\n');
    fprintf(filehandle, 'node[shape=box];\n\n');

    % define the root node
    node_number = 0;
    p = get(args.object, 'Parent');
    % define root element
    type = get(p, 'Type');
    fprintf(filehandle, 'N%d [label="%s"]\n\n', node_number, type);

    % start recursion
    plot_children(filehandle, p, node_number);

    % finish off
    fprintf(filehandle, '}');

    % ----------------------------------------------------------------------------
    function plot_children(fh, h, parent_node)

        children = allchild(h);

        for h = children(:)'
            if shouldSkip(h), continue, end;
            node_number = node_number + 1;

            label = {};
            label = addHGProperty(label, h, 'Type', '');
            try
                hClass = class(handle(h));
                label = addProperty(label, 'Class', hClass);
            catch
                % don't do anything
            end
            label = addProperty(label, 'Handle', sprintf('%g', double(h)));
            label = addHGProperty(label, h, 'Title', '');
            label = addHGProperty(label, h, 'Axes', []);
            label = addHGProperty(label, h, 'String', '');
            label = addHGProperty(label, h, 'Tag', '');
            label = addHGProperty(label, h, 'DisplayName', '');
            label = addHGProperty(label, h, 'Visible', 'on');
            label = addHGProperty(label, h, 'HandleVisibility', 'on');

            % print node
            fprintf(fh, 'N%d [label="%s"]\n', ...
                    node_number, m2tstrjoin(label, '\n'));

            % connect to the child
            fprintf(fh, 'N%d -> N%d;\n\n', parent_node, node_number);

            % recurse
            plot_children(fh, h, node_number);
        end
    end
end
% ==============================================================================
function bool = shouldSkip(h)
    %  returns TRUE for objects that can be skipped
    objType = get(h, 'Type');
    bool = ismember(lower(objType), guitypes());
end
% ==============================================================================
function label = addHGProperty(label, h, propName, default)
    % get a HG property and assign it to a GraphViz node label
    if ~exist('default','var') || isempty(default)
        shouldOmit = @isempty;
    elseif isa(default, 'function_handle')
        shouldOmit = default;
    else
        shouldOmit = @(v) isequal(v,default);
    end

    if isprop(h, propName)
        propValue = get(h, propName);
        if numel(propValue) == 1 && ishghandle(propValue) && isprop(propValue, 'String')
            % dereference Titles, labels, ...
            propValue = get(propValue, 'String');
        elseif ishghandle(propValue)
            % dereference other HG objects to their raw handle value (double)
            propValue = double(propValue);
        elseif iscell(propValue)
            propValue = ['{' m2tstrjoin(propValue,',') '}'];
        end

        if ~shouldOmit(propValue)
            label = addProperty(label, propName, propValue);
        end
    end
end
function label = addProperty(label, propName, propValue)
    % add a property to a GraphViz node label
    if isnumeric(propValue)
        propValue = num2str(propValue);
    elseif iscell(propValue)
        propValue = m2tstrjoin(propValue,sprintf('\n'));
    end
    label = [label, sprintf('%s: %s', propName, propValue)];
end
% ==============================================================================


text:
%
%  Univariate factorization that is capable of identifying multiple roots
%  and multiplicities accurately even if the polynomial is perturbed
%
%   Input    p -- coefficient vector of the polynomial to be factored
%          tol -- (optional) backward error tolerance
%    showroots -- (optional, 0 or 1) showing roots or not
%
%  Output    F -- kx3 matrix containing factors with
%                    each row [a,b,m] representing a factor (a*x+b)^m
%
%   Syntax:  >> LHS = RHS
%            where LHS and RHS can be any one of the following:
%
%               LHS options:   |  RHS options:
%               ---------------|------------------
%                 F            |  uvFactor(p)
%                [F,res]       |  uvFactor(p,tol)
%                [F,res,fcnd]  |  uvFactor(p,tol,1)
%

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function Population = EnvironmentalSelection(Population,V,theta)
% The environmental selection of RVEA*

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Select only the non-dominated solutions
    Population = Population(NDSort(Population.objs,1)==1);
    PopObj = Population.objs;
    [N,M]  = size(PopObj);
    NV     = size(V,1);

    %% Translate the population
    PopObj = PopObj - repmat(min(PopObj,[],1),N,1);

    %% Calculate the degree of violation of each solution
    CV = sum(max(0,Population.cons),2);

    %% Calculate the smallest angle value between each vector and others
    cosine = 1 - pdist2(V,V,'cosine');
    cosine(logical(eye(length(cosine)))) = 0;
    gamma  = min(acos(cosine),[],2);

    %% Associate each solution to a reference vector
    Angle = acos(1-pdist2(PopObj,V,'cosine'));
    [~,associate] = min(Angle,[],2);

    %% Select one solution for each reference vector
    Next = zeros(1,NV);
    for i = unique(associate)'
        current1 = find(associate==i & CV==0);
        current2 = find(associate==i & CV~=0);
        if ~isempty(current1)
            % Calculate the APD value of each solution
            APD = (1+M*theta*Angle(current1,i)/gamma(i)).*sqrt(sum(PopObj(current1,:).^2,2));
            % Select the one with the minimum APD value
            [~,best] = min(APD);
            Next(i)  = current1(best);
        elseif ~isempty(current2)
            % Select the one with the minimum CV value
            [~,best] = min(CV(current2));
            Next(i)  = current2(best);
        end
    end
    % Population for next generation
    Population = Population(Next(Next~=0));
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%% Importing Vectors
%
%%
% Large lists of vectors can be imported from a text file by the command

fname = fullfile(mtexDataPath,'vector3d','vectors.txt');
v = vector3d.load(fname,'ColumnNames',{'polar angle','azimuth angle'})

%%
% In order to visualize large lists of specimen directions scatter plots

scatter(v,'upper')

%%
% or contour plots may be helpful

contourf(v,'upper')

author: mtex-toolbox

repo: mtex

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


text:
function [f,fu,fx] = f_fu_fx_Wrapper(u,x,p,parIdx)
% fu and fx are calculated using finite difference by default

% dx/dt = f(u,x,p)
% parIdx: index of the core (for reentrant purpose)

    [xDim,~] = size(x);
    [uDim,~] = size(u);
    f = f_Wrapper(u,x,p,parIdx);
    fu = zeros(xDim,uDim);
    fx = zeros(xDim,xDim);
    h  = 1e-8;
    % fu
    for i=1:uDim
         ei = zeros(uDim,1);
         ei(i,1) = 1;
         fu(:,i) = (f_Wrapper(u+ei*h,x,p,parIdx) - f)/h;
    end
    % fx
    for i=1:xDim
         ei = zeros(xDim,1);
         ei(i,1) = 1;
         fx(:,i) = (f_Wrapper(u,x+ei*h,p,parIdx) - f)/h;
    end
end

author: deng-haoyang

repo: ParNMPC


text:
function [TruePositives, FalseNegatives] = testPutrefactionPathways(model, microbeID, biomassReaction, database)
% Performs an FVA and reports those putrefaction pathway end reactions (exchange reactions)
% that can carry flux in the model and should carry flux according to
% data (true positives) and those putrefaction pathway end reactions that
% cannot carry flux in the model but should be secreted according to in
% vitro data (false negatives).
%
% INPUT
% model             COBRA model structure
% microbeID         Microbe ID in carbon source data file
% biomassReaction   Biomass objective functions (low flux through BOF
%                   required in analysis)
% database          Structure containing rBioNet reaction and metabolite
%                   database
%
% OUTPUT
% TruePositives     Cell array of strings listing all putrefaction reactions
% that can carry flux in the model and in in vitro data.
% FalseNegatives    Cell array of strings listing all putrefaction reactions
% that cannot carry flux in the model but should carry flux according to in
% vitro data.
%
% .. Author:
%      Almut Heinken, Dec 2017
%                     March 2022 - changed code to string-matching to make
%                     it more robust

global CBT_LP_SOLVER
if isempty(CBT_LP_SOLVER)
    initCobraToolbox
end

% read putrefaction product tables
dataTable = readInputTableForPipeline('PutrefactionTable.txt');

corrRxns = {'Histidine degradation (histidine -> glutamate)','GLUFORT';'THF production (histidine -> tetrahydrofolate)','FTCD';'Glutamate production (glutamate -> acetate + pyruvate)','CITMALL';'Putrescine_1','EX_ptrc(e)';'Putrescine_2','EX_ptrc(e)';'Putrescine_3','EX_ptrc(e)';'Spermidine/ Spermine production (methionine -> spermidine)','EX_spmd(e)';'Cadaverine production (lysine -> cadaverine)','EX_15dap(e)';'Cresol production (tyrosine -> cresol)','EX_pcresol(e)';'Indole production (tryptophan -> indole)','EX_indole(e)';'Phenol production (tyrosine -> phenol)','EX_phenol(e)';'H2S_1','EX_h2s(e)';'H2S_2','EX_h2s(e)';'H2S_3','EX_h2s(e)';'H2S_4','EX_h2s(e)';'H2S_5','EX_h2s(e)'};

TruePositives = {};  % true positives (uptake in vitro and in silico)
FalseNegatives = {};  % false negatives (uptake in vitro not in silico)

% find microbe index in putrefaction table
mInd = find(strcmp(dataTable(:,1), microbeID));
if isempty(mInd)
    warning(['Microbe "', microbeID, '" not found in putrefaction product data file.'])
else
    % perform FVA to identify uptake metabolites
    % set BOF
    if ~any(ismember(model.rxns, biomassReaction)) || nargin < 3
        error(['Biomass reaction "', biomassReaction, '" not found in model.'])
    end
    model = changeObjective(model, biomassReaction);
    % set a low lower bound for biomass
    %     model = changeRxnBounds(model, biomassReaction, 1e-3, 'l');
    % list exchange reactions
    exchanges = model.rxns(strncmp('EX_', model.rxns, 3));
    % open all exchanges
    model = changeRxnBounds(model, exchanges, -1000, 'l');
    model = changeRxnBounds(model, exchanges, 1000, 'u');

    % get the reactions to test
    rxns = {};
    for i=2:size(dataTable,2)
        if contains(version,'(R202') % for Matlab R2020a and newer
            if dataTable{mInd,i}==1
                findCorrRxns = find(strcmp(corrRxns(:,1),dataTable{1,i}));
                rxns = union(rxns,corrRxns(findCorrRxns,2:end));
            end
        else
            if strcmp(dataTable{mInd,i},'1')
                findCorrRxns = find(strcmp(corrRxns(:,1),dataTable{1,i}));
                rxns = union(rxns,corrRxns(findCorrRxns,2:end));
            end
        end
    end

    % flux variability analysis on reactions of interest
    rxns = unique(rxns);
    rxns = rxns(~cellfun('isempty', rxns));
    rxnsInModel=intersect(rxns,model.rxns);
    if ~isempty(rxnsInModel)
        currentDir=pwd;
        try
            [minFlux, maxFlux, ~, ~] = fastFVA(model, 0, 'max', 'ibm_cplex', ...
                rxnsInModel, 'S');
        catch
            warning('fastFVA could not run, so fluxVariability is instead used. Consider installing fastFVA for shorter computation times.');
            cd(currentDir)
            [minFlux, maxFlux] = fluxVariability(model, 0, 'max', rxnsInModel);
        end

        % active flux
        flux = rxnsInModel(maxFlux > 1e-6);
    else
        flux = {};
    end

    % which reaction should carry flux according to in vitro data
    for i=2:size(dataTable,2)
        rxn={};
        if contains(version,'(R202') % for Matlab R2020a and newer
            if dataTable{mInd,i}==1
                rxn = corrRxns{find(strcmp(corrRxns(:,1),dataTable{1,i})),2};
            end
        else
            if strcmp(dataTable{mInd,i},'1')
                rxn = corrRxns{find(strcmp(corrRxns(:,1),dataTable{1,i})),2};
            end
        end
        if ~isempty(rxn)
            % add any that are not in model/not carrying flux to the false negatives
            if ~isempty(intersect(rxn,flux))
                TruePositives = union(TruePositives,rxn);
            else
                FalseNegatives=union(FalseNegatives,rxn);
            end
        end
    end
end

% replace reaction IDs with metabolite names
if ~isempty(TruePositives)
    TruePositives = TruePositives(~cellfun(@isempty, TruePositives));
    TruePositives=strrep(TruePositives,'EX_','');
    TruePositives=strrep(TruePositives,'(e)','');

    for i=1:length(TruePositives)
        if ~isempty(find(strcmp(database.metabolites(:,1),TruePositives{i})))
            TruePositives{i}=database.metabolites{find(strcmp(database.metabolites(:,1),TruePositives{i})),2};
        end
    end
end

% warn about false negatives
if ~isempty(FalseNegatives)
    FalseNegatives = FalseNegatives(~cellfun(@isempty, FalseNegatives));
    FalseNegatives=strrep(FalseNegatives,'EX_','');
    FalseNegatives=strrep(FalseNegatives,'(e)','');
    for i = 1:length(FalseNegatives)
        if ~isempty(find(strcmp(database.metabolites(:,1),FalseNegatives{i})))
        FalseNegatives{i}=database.metabolites{find(strcmp(database.metabolites(:,1),FalseNegatives{i})),2};
        end
    end
end

end

author: opencobra

repo: cobratoolbox


text:
function geometry_test062 ( )

%*****************************************************************************80
%
%% TEST062 tests PLANE_NORMAL_TRIANGLE_INT_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;
  ntest = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST062\n' );
  fprintf ( 1, '  PLANE_NORMAL_TRIANGLE_INT_3D finds the\n' );
  fprintf ( 1, '    intersection points of a normal form plane\n' );
  fprintf ( 1, '    and a triangle.\n' );

  normal(1:dim_num) =  [ 1.0, -2.0, -3.0 ];

  pp(1:dim_num) = [ 0.0, 0.0, 2.0 ];

  r8vec_print ( dim_num, pp, '  The point PP:' );

  r8vec_print ( dim_num, normal, '  The normal vector N:' );

  for i = 1 : ntest

    if ( i == 1 )
      t(1:dim_num,1:3) = [ ...
         3.0,  0.0, -7.0; ...
        13.0, -4.0, -1.0; ...
         5.0,  1.0, -2.0 ]';
    elseif ( i == 2 )
      t(1:dim_num,1:3) = [ ...
         3.0,  0.0, -7.0; ...
        13.0, -4.0, -1.0; ...
         9.0,  3.0,  8.0 ]';
    elseif ( i == 3 )
      t(1:dim_num,1:3) = [ ...
        -6.0,  0.0,  0.0; ...
         0.0,  3.0,  0.0; ...
         0.0,  0.0,  2.0 ]';
    elseif ( i == 4 )
      t(1:dim_num,1:3) = [ ...
        -4.0,  1.0,  0.0; ...
         0.0,  6.0, -2.0; ...
         0.0,  0.0,  1.0 ]';
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Case %d\n', i );

    r8mat_transpose_print ( dim_num, 3, t, '  Triangle vertices:' );

    [ num_int, pint ] = plane_normal_triangle_int_3d ( pp, normal, t );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Number of intersection points is %d\n', num_int );

    r8mat_transpose_print ( dim_num, num_int, pint, '  Intersection points:' );

  end

  return
end


text:
function [ flag ed ] = edf(matr,eds,vexs,ed,tem)
flag = 1;
for i = 2:eds
    [dvex f]=flecvexf(matr,i,vexs,eds,ed,tem);
    if f==1;
        flag = 0;
        break;
    end
    if dvex ~= 0
        ed(:,i)=[vexs(1,i) dvex];
        vexs(1,i+1)=dvex;
        matr(vexs(1,i+1),vexs(1,i))=0;
    else
        break;
    end
end

end

author: HuangCongQing


text:
%-------------------------------------------------------------------------------------------------------------
% This is an implementation of the TWSC algorithm for real-world image denoising
%
% Author:  Jun Xu, csjunxu@comp.polyu.edu.hk / nankaimathxujun@gmail.com
%          The Hong Kong Polytechnic University
%
% Please refer to the following paper if you find this code helps:
%
% @article{TWSC_ECCV2018,
% 	author = {Jun Xu and Lei Zhang and David Zhang},
% 	title = {A Trilateral Weighted Sparse Coding Scheme for Real-World Image Denoising},
% 	journal = {ECCV},
% 	year = {2018}
% }
%
% Please see the file License.txt for the license governing this code.
%-------------------------------------------------------------------------------------------------------------
clear;
Original_image_dir = '../dnd_2017/images_srgb/';
fpath = fullfile(Original_image_dir, '*.mat');
im_dir  = dir(fpath);
im_num = length(im_dir);
load '../dnd_2017/info.mat';

method = 'TWSC';
dataset = 'dnd_2017';
% write image directory
write_MAT_dir = ['../' dataset '_Results/'];
write_sRGB_dir = [write_MAT_dir method];
if ~isdir(write_sRGB_dir)
    mkdir(write_sRGB_dir)
end

% set parameters
Par.ps   = 6;       % patch size
Par.step = 3;       % the step of two neighbor patches
Par.win  = 20;      % size of window around the patch
Par.Outerloop = 8;  % iteration number of algorithm
Par.Innerloop = 2;  % iteration number of block matching
Par.nlspini = 70;   % initial number of patches
Par.display = 0;    %
Par.delta   = 0;    %
Par.nlspgap = 0;    % 10
Par.lambda1 = 0;    % set this parameter positive to apply TWSC model
Par.lambda2 = 3;    % parameter for estimating local noise level

alltime  = zeros(im_num, 1, 'double');
for i = 1 :im_num
    Par.image = i;
    load(fullfile(Original_image_dir, im_dir(i).name));
    S = regexp(im_dir(i).name, '\.', 'split');
    [h,w,ch] = size(InoisySRGB);
    % iterate over bounding boxes
    Idenoised_crop_bbs = cell(1,20);
    for j = 1:size(info(1).boundingboxes,1)
        Par.nlsp = Par.nlspini;  % number of non-local patches
        IMinname = [S{1} '_' num2str(j)];
        bb = info(i).boundingboxes(j,:);
        Par.nim = InoisySRGB(bb(1):bb(3), bb(2):bb(4),:);
        Par.I = Par.nim;
        % noise estimation
        for c = 1:ch
            Par.nSig(c) = NoiseEstimation(Par.nim(:, :, c)*255, Par.ps)/255;
        end
        % initial PSNR and SSIM
        fprintf('%s: \n', IMinname);
        % denoising
        t1=clock;
        [IMout, Par]  =  TWSC_Sigma_RW(Par);
        t2=clock;
        % etime(t2,t1)
        alltime(Par.image)  = etime(t2, t1);
        %% output
        IMoutname = sprintf([write_sRGB_dir '/' method '_' dataset '_' IMinname '.png']);
        imwrite(IMout, IMoutname);
        Idenoised_crop_bbs{j} = single(IMout);
    end
    for j = 1:size(info(1).boundingboxes,1)
        Idenoised_crop = Idenoised_crop_bbs{j};
        save(fullfile(write_MAT_dir, sprintf('%04d_%02d.mat', i, j)), 'Idenoised_crop');
    end
    fprintf('Image %d/%d done\n', i,50);
end
% generate submission files
bundle_submission_srgb( write_MAT_dir );

author: csjunxu


text:
function q = ik (K)
% Anthropomorphic arm with 3 DOF
% It calculates the Inverse Kinematic of an Anthropomorphic arm with 3 DOF.
% 'q' is the solutions in radiant and K is the direct Kinematic matrix.
%
%               K = [ n s a p;
%                     0 0 0 1]
% where n, s, a are three vectors fo 3 elements that represents the
% end-effector's orientation, and p is the desired end-effector position.

% Denavit-Hartenberg's Parameters
a1=0;           % [m]
a2=0.2;         % [m]
a3=0.2;         % [m]
alfa1=pi/2;     % [rad]
alfa2=0;        % [rad]
alfa3=0;        % [rad]

dk=K;          % Direct kinematics matrix

% Inverse Kinematic
pw_x=dk(1,4);   % Vector's components that representes the end-effector position
pw_y=dk(2,4);
pw_z=dk(3,4);

c3=(pw_x^2+pw_y^2+pw_z^2-a2^2-a3^2)/(2*a2*a3);  % cos(teta3)
s3=-sqrt(1-c3^2);        % sin(teta3)
teta3=atan2(s3,c3);

c2=(sqrt(pw_x^2+pw_y^2)*(a2+a3*c3)+pw_z*a3*s3)/(a2^2+a3^2+2*a2*a3*c3);      % cos(teta2)
s2=(pw_z*(a2+a3*c3)-sqrt(pw_x^2+pw_y^2)*a3*s3)/(a2^2+a3^2+2*a2*a3*c3);      % sin(teta2)
teta2=atan2((a2+a3*c3)*pw_z-a3*s3*sqrt(pw_x^2+pw_y^2),(a2+a3*c3)*sqrt(pw_x^2+pw_y^2)+a3*s3*pw_z);

teta1=atan2(pw_y,pw_x);

q=[teta1 teta2 teta3]';     % Solutions in radiant

author: Sable

repo: mcbench-benchmarks


text:
function batch_anonymize_CERR_plans(varargin)
%function batch_anonymize_CERR_plans(varargin)
%
%Type "function batch_anonymize_CERR_plans(varargin)" (without quotes) in Command window to run batch conversion. User will be
%prompted to select source and destination directories. This function anonymizes CERR plans under sourceDir and subdirectories
%and places them in destinationDir.
%
%APA, 10/05/2010
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

% Example-run
% sourceDir = 'J:\bioruser\apte\batch code';
% destinationDir = 'J:\bioruser\apte\batch code\OUT';
% batchConvert(sourceDir,destinationDir)

feature accel off

anonymize_file_name_flag = 1; % 1:anonymize file name with research number, 0: retain file name.
researchNumber = 1000; %naming will start with this number + 3.
replacementName = 'LungDoseResponse';

persistent convertedC planNameC
if isempty(varargin)
    %init_ML_DICOM
    convertedC = {'Source Directory:'};
    planNameC = {'Converted Plan Name:'};
    researchNumberC = {'MedPhys Research Number:'};
    sourceDir = uigetdir(cd, 'Select Source directory');
    if isnumeric(sourceDir)
        return;
    end
    destinationDir = uigetdir(cd, 'Select Destination directory');
    if isnumeric(destinationDir)
        return;
    end
else
    convertedC = {'Source Directory:'};
    planNameC = {'Converted Plan Name:'};
    researchNumberC = {'MedPhys Research Number:'};
    sourceDir = varargin{1};
    destinationDir = varargin{2};
end

if ispc
    slashType = '\';
else
    slashType = '/';
end
if strcmpi(sourceDir(end),slashType)
    sourceDir(end) = [];
end
slashIndex = strfind(sourceDir,slashType);
allDirS = dir(sourceDir);
for dirNum = 1:length(allDirS)
    researchNumber = researchNumber + 1;
    [pathStr,nameStr,extStr] = fileparts(allDirS(dirNum).name);
    [pathStr,nameStr,extStr2] = fileparts(nameStr);
    if ~allDirS(dirNum).isdir && (strcmpi(extStr,'.bz2') || strcmpi(extStr,'.mat')) && ~ismember(sourceDir,convertedC)

        if ~strcmpi(sourceDir(end),slashType)
            sourceDir = [sourceDir,slashType];
        end

        fileName = fullfile(sourceDir, allDirS(dirNum).name);

        try

            disp(['Converting ',fileName,' ...'])
            planC = loadPlanC(fileName, tempdir);

            %Anonymize
            PHI = {'studyNumberOfOrigin','PatientName','caseNumber','archive','PatientID','scanDate','DICOMHeaders'};

            if ~isempty(replacementName)

                for i = 1 : length(PHI)
                    str = PHI{i};
                    %planC = anonymize(planC,str,replacementName);
                    planC = anonCERRplanC(planC);
                end

            end

            %Check for duplicate name of sourceDirName
            if anonymize_file_name_flag
                newFileName = [num2str(researchNumber),extStr];
            else
                sourceDirName = sourceDir;
                dirOut = dir(destinationDir);
                allOutNames = {dirOut.name};
                newFileName = allDirS(dirNum).name;
                indexSlash = 1;
                while indexSlash~=length(slashIndex) && any(strcmpi(newFileName, allOutNames))
                    newFileName = [sourceDirName,'_',newFileName];
                    sourceDirName = sourceDir(slashIndex(end-indexSlash)+1:end);
                    sourceDirName(strfind(sourceDirName,slashType)) = deal('_');
                    indexSlash = indexSlash + 1;
                end
                if any(strcmpi(newFileName,allOutNames))
                    newFileName = [newFileName,'duplicate_',num2str(rand(1)),extStr2,extStr];
                end
            end

            %save_planC(planC,[], 'passed', fullfile(destinationDir,[sourceDirName,'.mat.bz2']));
            save_planC(planC,[], 'passed', fullfile(destinationDir,newFileName));

            clear planC
            convertedC{end+1} = fileName;
            planNameC{end+1} = newFileName;
            researchNumberC{end+1} = num2str(researchNumber);
        catch
            convertedC{end+1} = fileName;
            planNameC{end+1} = 'NOT CONVERTED';
            researchNumberC{end+1} = num2str(researchNumber);
            disp(['NOT CONVERTED ',fileName,' ...'])
        end

    elseif allDirS(dirNum).isdir && ~strcmp(allDirS(dirNum).name,'.') && ~strcmp(allDirS(dirNum).name,'..')
        batch_anonymize_CERR_plans(fullfile(sourceDir,allDirS(dirNum).name),destinationDir)
    end
end

% Save keys to Excel file
xlswrite(fullfile(destinationDir,'batch_convert_results.xls'),convertedC','Sheet1','A1')
xlswrite(fullfile(destinationDir,'batch_convert_results.xls'),planNameC','Sheet1','B1')
xlswrite(fullfile(destinationDir,'batch_convert_results.xls'),researchNumberC','Sheet1','C1')

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function q = Qconj(p)

% QCONJ   Quaternion conjugate
%
%   Q = QCONJ(P)returns the quaternion conjugate
%    q = CONJQ(p) returns the quaternion conjugate of the quaternion P. P
%    is a 4-vector representing a quaternion or an array 4*N (column i
%    represents quaternion i) where N is the number of quaternions. If P =
%    p0 + vecp, then P* = p0-vecp. Q has the same size as P.
%
% See also DQCONJ

sp = size(p);
if sp == [1 4], p = p'; sp = size(p); end

% wrong size
if sp(1) ~= 4
    error('DualQuaternion:Qconj:wrongsize',...
        '%d rows in array p. It should be 4.',sp(1));
end

q = p;
q(2:4,:) = - q(2:4,:);


text:
function [di, dj] = gradient_2d(M, mask)
% [di dj] = gradient_2d(M)
% [di dj] = gradient_2d(M, mask)
%
% Approximate the gradient of a two-dimensional function represented by M(i,j), using
% a 3x3 plane fitting algorithm.
%
%                         PARAMETERS
%
% M Matrix of function values, M(i,j) = f(i,j), the spacing is assumed to be uniform.
%
% mask (Optional) Matrix of the same size of M, with zero values in positions of M to be ignored by the
%                 gradient computations.
%
%                        RETURN VALUES
%
% di First component of the gradient: di(i,j) is the approximated [d/di f](i,j)
%
% dj Second component of the gradient: dj(i,j) is the approximated [d/dj f](i,j)
%
%                        LIMITATIONS
%
% To avoid ill-conditioned solutions, the mask is ignored near the "edges" of M,
% so this may result in artifacts in some cases. Leaving a "safety" border
% on the edges should be more than enough to prevent the issue.
%
% Author: Luca Vezzaro (elvezzaro@gmail.com)

	nr = size(M, 1);
	nc = size(M, 2);

	if nargin < 2
		mask = ones(nr, nc);
	end

	v = zeros(3,1);

	% The used algorithm is simple: the evaluation point is set
	% as the origin, and a plane is fitted to the 8 nearest positions,
	% centered on the origin. On the boundary, there are less values,
	% but the method is the same.

	% First the "corners"

	p(1,:) = [1 0]; v(1) = M(2,1) - M(1,1);
	p(2,:) = [0 1]; v(2) = M(1,2) - M(1,1);
	p(3,:) = [1 1]; v(3) = M(2,2) - M(1,1);

	n = p \ v;

	di(1,1) = n(1);
	dj(1,1) = n(2);

	p(1,:) = [0 -1]; v(1) =  M(1,nc-1) - M(1,nc);
	p(2,:) = [1 -1]; v(2) =  M(2,nc-1) - M(1,nc);
	p(3,:) = [1 0];  v(3) =  M(2,nc)   - M(1,nc);

	n = p \ v;

	di(1,nc) = n(1);
	dj(1,nc) = n(2);

	p(1,:) = [-1 0]; v(1) = M(nr-1,1) - M(nr,1);
	p(2,:) = [-1 1]; v(2) = M(nr-1,2) - M(nr,1);
	p(3,:) = [0 1];  v(3) = M(nr,2)   - M(nr,1);

	n = p \ v;

	di(nr,1) = n(1);
	dj(nr,1) = n(2);

	p(1,:) = [-1 -1]; v(1) = M(nr-1,nc-1) - M(nr,nc);
	p(2,:) = [0 -1];  v(2) = M(nr,nc-1)   - M(nr,nc);
	p(3,:) = [-1 0];  v(3) = M(nr-1,nc)   - M(nr,nc);

	n = p \ v;

	di(nr,nc) = n(1);
	dj(nr,nc) = n(2);

	% Then, the boundaries

	for j=2:nc-1
		p(1,:) = [0 -1]; v(1) = M(1,j-1) - M(1,j);
		p(2,:) = [1 -1]; v(2) = M(2,j-1) - M(1,j);
		p(3,:) = [1 0];  v(3) = M(2,j)   - M(1,j);
		p(4,:) = [0 1];  v(4) = M(1,j+1) - M(1,j);
		p(5,:) = [1 1];  v(5) = M(2,j+1) - M(1,j);

		n = p \ v;

		di(1,j) = n(1);
		dj(1,j) = n(2);

		p(1,:) = [-1 -1]; v(1) = M(nr-1,j-1) - M(nr,j);
		p(2,:) = [0 -1];  v(2) = M(nr,j-1)   - M(nr,j);
		p(3,:) = [-1 0];  v(3) = M(nr-1,j)   - M(nr,j);
		p(4,:) = [-1 1];  v(4) = M(nr-1,j+1) - M(nr,j);
		p(5,:) = [0 1];   v(5) = M(nr,j+1)   - M(nr,j);

		n = p \ v;

		di(nr,j) = n(1);
		dj(nr,j) = n(2);
	end

	for i=2:nr-1
		p(1,:) = [-1 0]; v(1) = M(i-1,1) - M(i,1);
		p(2,:) = [1 0];  v(2) = M(i+1,1) - M(i,1);
		p(3,:) = [-1 1]; v(3) = M(i-1,2) - M(i,1);
		p(4,:) = [0 1];  v(4) = M(i,2)   - M(i,1);
		p(5,:) = [1 1];  v(5) = M(i+1,2) - M(i,1);

		n = p \ v;

		di(i,1) = n(1);
		dj(i,1) = n(2);

		p(1,:) = [-1 -1]; v(1) = M(i-1,nc-1) - M(i,nc);
		p(2,:) = [0 -1];  v(2) = M(i,nc-1)   - M(i,nc);
		p(3,:) = [1 -1];  v(3) = M(i+1,nc-1) - M(i,nc);
		p(4,:) = [-1 0];  v(4) = M(i-1,nc)   - M(i,nc);
		p(5,:) = [1 0];   v(5) = M(i+1,nc)   - M(i,nc);

		n = p \ v;

		di(i,nc) = n(1);
		dj(i,nc) = n(2);
	end

	% And finally the remaining function values
	for i=2:nr-1
		for j=2:nc-1
			if mask(i,j) ~= 0
				k = 1;

				if mask(i-1,j-1) ~= 0
					p(k,:) = [-1 -1]; v(k) = M(i-1,j-1) - M(i,j);
					k = k + 1;
				end

				if mask(i-1,j) ~= 0
					p(k,:) = [-1 0];  v(k) = M(i-1,j)   - M(i,j);
					k = k + 1;
				end

				if mask(i-1,j+1) ~= 0
					p(k,:) = [-1 1];  v(k) = M(i-1,j+1) - M(i,j);
					k = k + 1;
				end

				if mask(i,j-1) ~= 0
					p(k,:) = [0 -1];  v(k) = M(i,j-1)   - M(i,j);
					k = k + 1;
				end

				if mask(i,j+1) ~= 0
					p(k,:) = [0 1];   v(k) = M(i,j+1)   - M(i,j);
					k = k + 1;
				end

				if mask(i+1,j-1) ~= 0
					p(k,:) = [1 -1];  v(k) = M(i+1,j-1) - M(i,j);
					k = k + 1;
				end

				if mask(i+1,j) ~= 0
					p(k,:) = [1 0];   v(k) = M(i+1,j)   - M(i,j);
					k = k + 1;
				end

				if mask(i+1,j+1) ~= 0
					p(k,:) = [1 1];   v(k) = M(i+1,j+1) - M(i,j);
					k = k + 1;
				end

				if k > 2
					n = p(1:k-1,:) \ v(1:k-1);

					di(i,j) = n(1);
					dj(i,j) = n(2);
				else
					di(i,j) = 0;
					dj(i,j) = 0;
				end
			else
				di(i,j) = 0;
				dj(i,j) = 0;
			end
		end
	end


text:
function gX = velotransKernGradX(kern, X, X2)

% VELOTRANSKERNGRADX Gradient of VELOTRANS kernel with respect to a point x.
% FORMAT
% DESC computes the gradient of the velocity translate
% kernel with respect to the input positions.
% ARG kern : kernel structure for which gradients are being
% computed.
% ARG x : locations against which gradients are being computed.
% RETURN g : the returned gradients. The gradients are returned in
% a matrix which is numData x numInputs x numData. Where numData is
% the number of data points and numInputs is the number of input
% dimensions in X.
%
% FORMAT
% DESC computes the gradident of the velocity translate
% kernel with respect to the input positions where both the row
% positions and column positions are provided separately.
% ARG kern : kernel structure for which gradients are being
% computed.
% ARG x1 : row locations against which gradients are being computed.
% ARG x2 : column locations against which gradients are being computed.
% RETURN g : the returned gradients. The gradients are returned in
% a matrix which is numData2 x numInputs x numData1. Where numData1 is
% the number of data points in X1, numData2 is the number of data
% points in X2 and numInputs is the number of input
% dimensions in X.
%
% SEEALSO velotransKernParamInit, kernGradX, velotransKernDiagGradX, translateKernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2011

% KERN

  t = X(:, end);
  xPass = X(:, 1:end-1) - t*kern.velocity;
  t2 = X2(:, end);
  x2Pass = X2(:, 1:end-1) - t2*kern.velocity;
  gXRecover = cmpndKernGradX(kern, xPass, x2Pass);
  gX = zeros(size(X2, 1), size(X2, 2), size(X, 1));
  for i = 1:size(X, 1)
    gX(:, 1:end-1, i) = gXRecover(:, :, i);
    gX(:, end, i) = -gXRecover(:, :, i)*kern.velocity'; % Another gradient
                                                        % guess
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [C1, C2, U1, U2, H, K, N] = meshCurvatures(vertices, faces, varargin)
%MESHCURVATURES Compute principal curvatures on mesh vertices.
%
%   [C1, C2] = meshCurvatures(VERTICES, FACES)
%   Computes the principal curvatures C1 and C2 for each vertex of the mesh
%   defined by VERTICES and FACES.
%
%   [C1, C2] = meshCurvatures(..., PNAME, PVALUE)
%   Provides additional input arguments based on a list of name-value pairs
%   of arguments. Parameter names can be:
%   * 'SmoothingSteps'      (integer, default: 3)
%       Specifies the number of steps for smoothing vertex curvature
%       tensors.
%   * 'Verbose'             (boolean, default: true)
%       Displays details about algorithm processing.
%   * 'ShowProgress'        (boolean, default: true)
%       Displays a text-based progress bar.
%
%   Algorithm
%   The function is adapted from the "compute_curvature" function, in the
%   "toolbox_graph" fro Gabriel Peyre.
%   The basic idea is to define a curvature tensor for each edge, by
%   assigning a minimum curvature equal to zero along the edge, and a
%   maximum curvature equal to the dihedral angle across the edge.
%   Averaging around the neighbors of a vertex v yields a summation formula
%   over the neighbor edges to compute the curvature tensor of a vertex:
%           1
%   C(v) = ----     Sum       \beta(e) || e \cap A(v) || ebar ebar^t
%          A(v)  {e \in A(v)}
%   where:
%   * A(v) is the neighborhood region, usually defined as a 'ring' around
%       the vertex v
%   * beta(e) is the dihedral angle between the normals of the two faces
%       incident to edge e
%   * || e \cap A(v) || is the length of e (more exactly, the length of the
%       part of e contained within the neighborhood region
%   * ebar is the normalized edge
%
%   The curvature tensor is then decomposed into C = P D P^-1, with P
%   containing main direction vectors and normal, and D being a diagonal
%   matrix with the two main curvatures and zero along the diagonal.
%
%   References
%   * David Cohen-Steiner and Jean-Marie Morvan (2003).
%       "Restricted Delaunay triangulations and normal cycle".
%       In Proc. 19th Annual ACM Symposium on Computational Geometry,
%       pages 237-246.
%   * Pierre Alliez, David Cohen-Steiner, Olivier Devillers, Bruno Levy,
%       and Mathieu Desbrun (2003). "Anisotropic Polygonal Remeshing".
%       ACM Transactions on Graphics.
%       (SIGGRAPH '2003 Conference Proceedings)
%   * Mario Botsch, Leif Kobbelt, M. Pauly, P. Alliez, B. Levy (2010).
%       "Polygon Mesh Processing", Taylor and Francis Group, New York.
%
%   Example
%     [v, f] = torusMesh;
%     f2 = triangulateFaces(f);
%     [c1, c2] = meshCurvatures(v, f2);
%     figure; hold on; axis equal; view(3);
%     drawMesh(v, f2, 'VertexColor', c1 .* c2);
%
%   See also
%     meshes3d, drawMesh, triangulateFaces
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2021-09-21, using Matlab 9.10.0.1684407 (R2021a) Update 3
% Copyright 2021-2022 INRAE - BIA Research Unit - BIBS Platform (Nantes)

%% Process input arguments

% default values for options
nIters = 3;
verbose = true;
showProgress = true;

while length(varargin) > 1
    name = varargin{1};
    if strcmpi(name, 'SmoothingSteps')
        nIters = varargin{2};
    elseif strcmpi(name, 'Verbose')
        verbose = varargin{2};
    elseif strcmpi(name, 'ShowProgress')
        showProgress = varargin{2};
    else
        error('Unknown option: %s', name);
    end
    varargin(1:2) = [];
end

% validate vertices
if ~isnumeric(vertices) || size(vertices, 2) ~= 3
    error('Requires vertices to be a N-by-3 numeric array');
end

% ensure faces are triangular
if ~isnumeric(faces) || size(faces, 2) > 3
    warning('requires triangle mesh, forces triangulation');
    faces = triangulateFaces(faces);
end

%% Retrieve adjacency relationships

if verbose
    disp('compute adjacencies');
end

% number of elements of each type
nv = size(vertices, 1);
nf = size(faces, 1);

% ev1 and ev2 are indices of source and target vertex of each edge
% (recomputed later)
ev1 = [faces(:,1); faces(:,2); faces(:,3)];
ev2 = [faces(:,2); faces(:,3); faces(:,1)];

% Compute sparse matrix representing edge-to-face adjacency
s = [1:nf 1:nf 1:nf]';
A = sparse(ev1, ev2, s, nv, nv);

% converts sparse matrix to indices of adjacent vertices and faces
[~, ~, ef1] = find(A);         % index of 'right' face
[ev1, ev2, ef2] = find(A');    % index of 'left' face, and of vertices

% edges are consdered twice (one for each vertex)
% keep only the edge with lower source index
inds = find(ev1 < ev2);
ef1 = ef1(inds);
ef2 = ef2(inds);
ev1 = ev1(inds);
ev2 = ev2(inds);

% number of edges
ne = length(ev1);

%% Compute geometry features

% compute edge direction vectors
edgeVectors = vertices(ev2,:) - vertices(ev1,:);

% normalize edge direction vecotrs
d = sqrt(sum(edgeVectors.^2, 2));
edgeVectors = bsxfun(@rdivide, edgeVectors, d);

% avoid too large numerics
d = d ./ mean(d);

% normals to faces
normals = meshFaceNormals(vertices, faces);

% ensure normals point outward the mesh
if meshVolume(vertices, faces) < 0
    normals = -normals;
end

% inner product of normals
dp = sum(normals(ef1, :) .* normals(ef2, :), 2);

% compute the (unsigned) dihedral angle between the normals of the two
% faces incident to each edge
beta = acos(min(max(dp, -1), 1));

% relatice orientation of face normals cross product and edge orientation
cp = crossProduct3d(normals(ef1, :), normals(ef2, :));
si = sign(sum(cp .* edgeVectors, 2));

% compute signed dihedral angle
beta = beta .* si;

%% Compute tensors

if verbose
    disp('compute edge tensors');
end

% curvature tensor of each edge
T = zeros(3, 3, ne);
for i = 1:3
    for j = 1:i
        T(i, j, :) = reshape(edgeVectors(:,i) .* edgeVectors(:,j), 1, 1, ne);
        T(j, i, :) = T(i, j, :);
    end
end
T = bsxfun(@times, T, reshape(d .* beta, [1 1 ne]));

% curvature tensor of each vertex by pooling edge tensors
Tv = zeros(3, 3, nv);
w = zeros(1, 1, nv);
for k = 1:ne
    if showProgress
        displayProgress(k, ne);
    end
    Tv(:,:,ev1(k)) = Tv(:,:,ev1(k)) + T(:,:,k);
    Tv(:,:,ev2(k)) = Tv(:,:,ev2(k)) + T(:,:,k);
    w(:,:,ev1(k)) = w(:,:,ev1(k)) + 1;
    w(:,:,ev2(k)) = w(:,:,ev2(k)) + 1;
end
w(w < eps) = 1;
Tv = Tv ./ repmat(w, [3 3 1]);

if verbose
    disp('average vertex tensors');
end

% apply smoothing on the tensor field
for i = 1:3
    for j = 1:3
        a = Tv(i, j, :);
        a = smoothMeshFunction(vertices, faces, a(:), nIters);
        Tv(i, j, :) = reshape(a, [1 1 nv]);
    end
end

%% Retrieve curvatures and eigen vectors from tensors

if verbose
    disp('retrieve curvatures');
end

% allocate memory
U = zeros(3, 3, nv);
D = zeros(3, nv);

% iterate over vertices
for k = 1:nv
    % display progress
    if showProgress
        displayProgress(k,nv);
    end

    % extract eigenvectors and eigenvalues for current vertex
    [u, d] = eig(Tv(:,:,k));
    d = real(diag(d));

    % sort acording to [normal, min curv, max curv]
    [~, I] = sort(abs(d));
    D(:, k) = d(I);
    U(:, :, k) = real(u(:,I));
end

% retrieve main curvatures and associated directions
C1 = D(2,:)';
C2 = D(3,:)';
U1 = squeeze(U(:,3,:))';
U2 = squeeze(U(:,2,:))';

% enforce C1 < C2
inds = find(C1 > C2);
C1tmp = C1;
U1tmp = U1;
C1(inds) = C2(inds);
C2(inds) = C1tmp(inds);
U1(inds,:) = U2(inds,:);
U2(inds,:) = U1tmp(inds,:);

% compute optional output arguments
if nargout > 4
    % average and gaussian curvatures
    H = (C1 + C2) / 2;
    K = C1 .* C2;

    if nargout > 6
        % normal vector for each vertex
        N = squeeze(U(:,1,:))';
    end
end

function displayProgress(n, N)
% Display the progress of current step using a text-based progress bar.
%
% based on the 'progressbar' function in G. Peyre's Graph Toolbox.

% width of the progress bar
w = 20;

% compute progress ratio as an integer betsween 0 and w
p = min( floor(n/N*(w+1)), w);

global pprev;
if isempty(pprev)
    pprev = -1;
end

if p ~= pprev
    str1 = repmat('*', 1, p);
    str2 = repmat('.', 1, w-p);
    str = sprintf('[%s%s]', str1, str2);
    if n > 1
        % clear previous string
        fprintf(repmat('\b', [1 length(str)]));
    end
    fprintf(str);
end

pprev = p;
if n == N
    fprintf('\n');
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function sdc_coeff = mfcc2sdc(CepCoeff,N,d,P,K)
%--------------------------------------------------------------------------
%Function for Shifted Delta Coefficient Computation.
%
% If you have any query or suggestion please mail to sahidullahmd@gmail.com
% Ussage: sdc_coeff = mfcc2sdc(CepCoeff,N,d,P,k)
%         CepCoeff: MFCC Coefficients stored in row-wise
%         N: NoOfCepstrum i.e. no of column of CepCoeff
%         d: Amount of shift for delta computation
%         P: Amount of shift for next frame whose deltas are to be
%            computed.
%         K: No of frame whose deltas are to be stacked.
%
%         sdc_coeff: Shifted delta coefficient of CepCoeff.
%                    Dimension of the output: NumberOfFrame x N*K
% Example:
% CepCoeff=rand(1000,19); % Randomly generate MFCC of 1000 frames
% sdc_coeff = mfcc2sdc(CepCoeff,7,1,3,7); %Compute SDC Coefficients.
%
%Details of shifted delta computation is available in W.M. Campbell, J.P.
%    Campbell, D.A. Reynolds, E. Singer, P.A. Torres-Carrasquillo, Support
%    vector machines for speaker and language recognition,
%    Computer Speech & Language, Volume 20, Issues 2-3,
%    Odyssey 2004: The speaker and Language Recognition
%    Workshop - Odyssey-04, April-July 2006, Pages 210-229.
%--------------------------------------------------------------------------
ToT=size(CepCoeff,1); %Actual Number of Frames
CepCoeff=(horzcat(CepCoeff', CepCoeff(1:P*(K-1),:)'))'; %Circular padding
[NoOfFrame NoOfCoeff]=size(CepCoeff); %Note the number of frames and
%                                             number of coefficients (N).
% % N=7; d=1; P=3; k=7;      %Default values are hard-coded if necessary.
delt=(mfcc2delta(CepCoeff,d))';   %Delta Feature Computation
sd_temp=cell(1,K);                %Preparation of a cell array for delta's
for i=1:K                         %For k number of shifts
    temp=delt(:,P*(i-1)+1:1:end); %P: Size of shift
    sd_temp{i}=temp(:,1:ToT);     %Take only desired (i.e. ToT) no of deltas.
end
sdc_coeff=cell2mat(sd_temp');     %Stacking the SDCs in a single variable
sdc_coeff=sdc_coeff';
%--------------------------------------------------------------------------


text:
function out = linearremap(x)
% LINEARREMAP Dumb linear remap
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

if ~isreal(x)
    out=abs(double(x));
else
    out=x;
end

end
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function pass = test_sum2( )
% Test spherefun sum2() command.

tol = 1000*chebfunpref().cheb2Prefs.chebfun2eps;

f = @(x,y,z) 1 + x + y + z;
g = spherefun(f);
exact_int = 4*pi;
pass(1) = abs(sum2(g) - exact_int) < tol;

f = @(x,y,z) 0.75*exp(-(9*x-2).^2/4 - (9*y-2).^2/4 - (9*z-2).^2/4)...
   + 0.75*exp(-(9*x+1).^2/49 - (9*y+1)/10 - (9*z+1)/10)...
   + 0.5*exp(-(9*x-7).^2/4 - (9*y-3).^2/4 - (9*z-5).^2/4)...
   - 0.2*exp(-(9*x-4).^2 - (9*y-7).^2 - (9*z-5).^2);
exact_int = 6.6961822200736179523;
g = spherefun(f);
pass(2) = abs(sum2(g) - exact_int) < tol;

end

author: chebfun

repo: chebfun


text:
function [fibercurvature] = dtiComputeFiberCurvature(fg)
%
% [fibercurvature] = dtiComputeFiberCurvature(fg)
%
% Compute the curvature of individual fibers and make a cell-array output
% of curvatures in each node.
%
% % INPUTS:
% fibers   = A fiber group.
%
% % OUTPUTS:
% fibercurvature = A cell array of fiber curvarture.
%
% This code is an updated version of dtiComputeFiberGroupCurvatures.m
%
% (C) Hiromasa Takemura, Stanford VISTA team, 2014

if notDefined('fg'), error('Fiber group required'); end

nfibers=size(fg.fibers, 1);

for i=1:nfibers

    numNodes(i) =size(fg.fibers{i}, 2);

    %Compute curvature representations
    fibercurvature{i}=dtiFiberCurvature(fg.fibers{i});
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function model = lfmExpandParam(model, params)

% LFMEXPANDPARAM Expand the given parameters into a LFM structure.
% FORMAT
% DESC takes the given vector of parameters and places them in the
% model structure, it then updates any stored representations that
% are dependent on those parameters, for example kernel matrices
% etc..
% ARG model : the model structure for which parameters are to be
% updated.
% ARG params : a vector of parameters for placing in the model
% structure.
% RETURN model : a returned model structure containing the updated
% parameters.
%
% SEEALSO : lfmCreate, lfmExtractParam, modelExtractParam, lfmUpdateKernels
%
% COPYRIGHT : Neil D. Lawrence, 2007

% KERN

params = real(params);
if isfield(model, 'fix')
  for i = 1:length(model.fix)
    params(model.fix(i).index) = model.fix(i).value;
  end
end

if length(params) ~= model.numParams
  error('Parameter vector is incorrect length');
end
startVal = 1;
endVal = model.kern.nParams;
model.kern = kernExpandParam(model.kern, params(startVal:endVal));

% The mass, springs and dampers  are actually stored in the kernel.
% We'll put them here as well for convenience.
for i = 1:model.kern.numBlocks
  model.mass(i) = model.kern.comp{i}.mass;
  model.spring(i) = model.kern.comp{i}.spring;
  model.damper(i) = model.kern.comp{i}.damper;
  model.sensitivity(i) = model.kern.comp{i}.sensitivity;
end

model = lfmUpdateKernels(model);
lengthObs = size(model.t, 1);
ind = 1:lengthObs;
for i = 1:model.numDisplacements
  if isfield(model, 'mu')
    model.m(ind) = model.y(ind) - model.mu(i)*ones(lengthObs, 1);
  else
    model.m(ind) = model.y(ind);
  end
  ind = ind + lengthObs;
end


text:
classdef SPX_ConjugateDescent < handle
% Conjugate gradient algorithm implementation

    properties
        % Settable and gettable properties

        % Maximum number of iterations for which the algorithm can run
        MaxIterations
        % Threshold of norm (in terms of percentage)
        NormThreshold
    end

    properties(SetAccess=private)
        % Gettable properties

        % The sparse symmetric positive definite matrix operator
        A
        % The vectors to be solved
        B
        % Problem dimension
        N
        % Number of equations to be solved
        S
        % Number of iterations taken for solving the problem
        Iterations
        % Residual norms at the end
        ResidualNorms
        % The solution vectors
        X
        % Residual vectors
        Residuals
        % Indicates if the problems converged
        % This would be false only if the algorithm crossed max iterations
        Converged
    end

    methods
        % Public methods
        function self = SPX_ConjugateDescent(A, B)
            % Constructor
            if isa(A, 'spx.dict.Operator')
                self.A = A;
            elseif ismatrix(A)
                self.A = spx.dict.MatrixOperator(A);
            else
                error('Unsupported operator.');
            end
            self.B = B;
            [self.N, self.S] = size(B);
            self.MaxIterations = self.N ;
            self.NormThreshold = 1e-6;
        end

        function result = solve(self)
            aa = self.A;
            bb = self.B;
            % Initial estimate vectors are all zeros
            xx = zeros(self.N, self.S);
            result = xx;
            self.Iterations = zeros(self.S, 1);
            self.ResidualNorms = zeros(self.S, 1);
            self.Converged = false(self.S, 1);
            % Initial residual vectors
            rr = bb  - aa * xx;
            self.Residuals  = rr;
            % Norm squared of initial residual vectors
            deltas = spx.norm.inner_product_cw(rr, rr);
            % The factor with which the norm needs to be reduced
            epsilon = self.NormThreshold;
            % Target limits on norm squared of residuals
            limits = epsilon^2 * deltas;
            % Maximum number of iterations for which the algorithm
            % is allowed to run
            imax = self.MaxIterations;
            % Number of problems being solved
            ns = self.S;
            for s=1:ns
                % Initialize iteration counter
                i = 0;
                % The quantities for this problem
                limit = limits(s);
                delta = deltas(s);
                % First residual
                r = rr(:, s);
                % First estimate
                x = xx(:, s);
                % Target b
                b = bb(:, s);
                % First direction
                d = r;
                while i < imax && delta > limit
                    % Compute the intermediate variable
                    q = aa * d;
                    % the line search scale factor in current direction
                    alpha = delta / (d' * q);
                    % Update estimate in current direction
                    x = x + alpha * d;
                    if mod(i , 50) == 0
                        % In order to avoid propagation of floating point
                        % errors, we will recompute the value of residual
                        r = b - aa * x;
                    else
                        % Otherwise we use a shortcut
                        r = r  - alpha * q;
                    end
                    % hold the current residual norm squared
                    delta_old = delta;
                    % Update residual norm squared
                    delta = r' * r;
                    % Compute the ratio
                    beta = delta / delta_old;
                    % choose the new direction
                    d = r + beta * d;
                    % Increase iteration counter
                    i = i + 1;
                end
                % The problem has been solved
                result(:, s) = x;
                % Number of iterations taken to solve this problem
                self.Iterations(s) = i;
                self.ResidualNorms(s) = sqrt(delta);
                self.Residuals(:, s) = r;
                self.Converged(s) = delta <= limit;
            end
            % Maintain the result for reference
            self.X = result;
        end

        function result = hasConverged(self)
            % Returns if all the solutions have converged
            result = all(self.Converged);
        end

        function printResults(self)
            ns = self.S;
            nn = self.N;
            for s = 1:ns
                fprintf('Problem: %d\n', s);
                fprintf('Iterations: %d\n', self.Iterations(s));
                fprintf('Residual norm: %.2f, Converged: %d\n', ...
                    self.ResidualNorms(s), self.Converged(s));
                if nn < 10
                    % We will print the solutions too
                    fprintf('Solution vector: ');
                    fprintf('%.4f ', self.X(:, s));
                    fprintf('\n');
                    fprintf('Residual vector: ');
                    fprintf('%.4f ', self.Residuals(:, s));
                    fprintf('\n');
                end
            end
        end
    end

    methods(Access=private)
        % Private methods

    end

    methods(Static)
        % Public static methods

    end

end

author: indigits

repo: sparse-plex


text:
function jsubs = cfg2jobsubs(item, subs)

% function jsubs = cfg2jobsubs(item, subs)
% Return the subscript into the job tree for a given subscript vector into
% the val part of the cfg tree. This generic function should be called only
% for leafs (cfg_entry, cfg_file, cfg_menu) of the cfg tree. It returns the
% subscripts that remain after item.val{1} has been addressed.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg2jobsubs.m 1862 2008-06-30 14:12:49Z volkmar $

rev = '$Rev: 1862 $'; %#ok

% Only de-reference subscripts into .val{1} field
if isequal(subs(1:2), substruct('.','val','{}',{1}))
    % return possible indices into .val{1} field
    jsubs = subs(3:end);
else
    cfg_message('matlabbatch:cfg2jobsubs:wrongsubs', 'Inappropriate subscript reference in item ''%s''.', item.tag);
    jsubs = struct('type',{},'subs',{});
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function g = polyardKernGradient(kern, x, varargin)

% POLYARDKERNGRADIENT Gradient of POLYARD kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the
% automatic relevance determination polynomial
% kernel's parameters. As well as the kernel structure and the
% input positions, the user provides a matrix PARTIAL which gives
% the partial derivatives of the function with respect to the
% relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x : the input locations for which the gradients are being
% computed.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in X.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x1 : the input locations associated with the rows of the
% kernel matrix.
% ARG x2 : the input locations associated with the columns of the
% kernel matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as X1 and the same number of columns
% as X2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% SEEALSO polyardKernParamInit, kernGradient, polyardKernDiagGradient, kernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006

% KERN

scales = sparse(diag(sqrt(kern.inputScales)));
xScale = x*scales;
if nargin < 4
  innerProd = xScale*xScale';
else
  xScale2 = varargin{1}*scales;
  innerProd = xScale*xScale2';
end
arg = kern.weightVariance*innerProd+kern.biasVariance;
base = kern.variance*kern.degree*arg.^(kern.degree-1);
baseCovGrad = base.*varargin{end};

g(1) = sum(sum(innerProd.*baseCovGrad));
g(2) = sum(sum(baseCovGrad));
g(3) = sum(sum(varargin{end}.*arg.^kern.degree));

if nargin < 4
  for j = 1:kern.inputDimension
    g(3+j) = sum(sum((x(:, j)*x(:, j)').*baseCovGrad))*kern.weightVariance;
  end
else
  for j = 1:kern.inputDimension
    g(3+j) = sum(sum((x(:, j)*varargin{1}(:, j)').*baseCovGrad))*kern.weightVariance;
  end
end


text:
function value = i4_divp ( i, j )

%*****************************************************************************80
%
%% I4_DIVP returns the smallest multiple of J greater than or equal to I.
%
%  Example:
%
%    I  J  VALUE
%
%    0  4    0
%    1  4    1
%    2  4    1
%    3  4    1
%    4  4    1
%    5  4    2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the number to be analyzed.
%
%    Input, integer J, the number, multiples of which will
%    be compared against I.  J may not be zero.
%
%    Output, integer VALUE, the smallest multiple of J that
%    is greater than or equal to I.
%
  if ( j ~= 0 )
    value = 1 + floor ( ( i - 1 ) / j );
  else
    value = 0;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4_DIVP - Fatal error!\n' );
    fprintf ( 1, '  The input value of J was zero!\n' );
    error ( 'I4_DIVP - Fatal error!\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function fm=oz(x)
%OZ FM=OZ(X) is the nonlinear residual for the
%   Ornstein-Zernike equations.
%
% Evaluate the nonlinearity as a sum of the substitution
% and compact parts.
%
fm=ozsub(x)+ozinteg(x);
%
% This is the substitution part.
%
function km=ozsub(x)
global L U rho
n=length(x); n2=n/2; h=x(1:n2); c=x(n2+1:n);
nl=U.*exp(h-c);
mup=h+1- nl;
mdown=c+1-nl;
km=[mup',mdown']';
%
% This is the compact part of the nonlinearity. The upper
% component is zero, the lower is rho*L*c*h.
%
% rho and L are scalars for this problem
%
function km=ozinteg(x)
global L U rho
L3=L*L*L*rho;
n=length(x); n2=n/2; h=x(1:n2); c=x(n2+1:n);
kd=L3*convk(c,h);
kd(1)=2*kd(2)-kd(3);
kd(n2)=2*kd(n2-1)-kd(n2-2);
kup=zeros(n2,1);
km=[kup',kd']';
%
% Discrete convolution with the Hankel transform.
%
function hc=convk(h,c)
Lx=1.0;
n=length(h);
dr=Lx/n;
th=hank2(h); tc=hank2(c); fhc=th.*tc; hc=ihank2(fhc);
%
% Hankel transform using the fast sine transform.
%
function hf=hank2(f)
nf=length(f); n=nf-1; m=nf-2;
h=1/n; beta=2*(m+1)*(h^3); ff=beta;
p=1:m; p=p';
hft=ff*lsint(p.*f(2:n))./p;
hf=[0,hft',0]';
%
% Inverse Hankel transform using the fast sine transform.
%
function ihf=ihank2(f)
nf=length(f); n=nf-1; m=nf-2;
h=1/n; beta=2*(m+1)*(h^3); ff=2/(n*beta);
p=1:m; p=p';
hft=ff*lsint(p.*f(2:n))./p;
p1=1:n; p1=p1'; p1=p1.*p1; p1(n)=p1(n)/2;
ihf=[0,hft',0]';
% LSINT
% Fast sine transform with MATLAB's FFT.
%
function lf=lsint(f)
n=length(f);
ft=-fft([0,f']',2*n+2);
lf=imag(ft(2:n+1));

repo: Aero-matlab


text:
%Trajectory_Examples.m
%
% This script contains a few sample "roads". SetPts are a series of points
% along the road, to which a polynomial is fit. Therer are also a few
% simulation parameters associated with each example.
%

Params.Sim.Example_Number = Example_Number;
switch Example_Number
    case 1
        Params.Traj.Func = false;
        Params.Traj.SetPts = [...
            0,0;
            10,5;
            20,9;
            40,9;
            60,0;
            80,-10];
        Params.Traj.Order = 4;
        Params.Sim.IC_Error = [0;0;0;0];
        Max_Pos_Err = 0.5;
        Params.Ctl.Minimum_Goal_Distance = 1; %(m) how far to "look ahead"
        Simulation_Duration = 60;
        Control_Frequency = 50;   %Hz

    case 2
        %This one appears to work reasonably well
        Params.Traj.Func = false;
        Params.Traj.SetPts = [...
            0,0;
            10,5;
            20,9;
            40,9;
            60,0;
            65,-2.5;
            70, -5;
            75,-8
            80,-12];
        Params.Traj.Order = 4;
        Params.Sim.IC_Error = [0;0;0;0];
        Max_Pos_Err = 0.5;
        Params.Ctl.Minimum_Goal_Distance = 1; %(m) how far to "look ahead"
        Simulation_Duration = 50;
        Control_Frequency = 50;   %Hz

    case 3
        Params.Traj.Func = false;
        Params.Traj.SetPts = [...
            0 0;
            5 0;
            10 1;
            15 3;
            20 2;
            25 0;
            30 -3;
            35 -8;
            40 -13;
            45 -20;
            50 -30;
            55 -40];
        Params.Traj.Order = 4;
        Params.Sim.IC_Error = [0;0;0;0];
        Max_Pos_Err = 1.5;
        Params.Ctl.Minimum_Goal_Distance = 1; %(m) how far to "look ahead"
        Simulation_Duration = 125;
        Control_Frequency = 40;   %Hz

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 4
        %USED IN PRESENTATION - MAIN RESULT

	    Params.Traj.Func = true;
        x = linspace(0,150,Params.Traj.Npts+1)';
        Params.Traj.x = x;
        Params.Traj.y = sin(x/20).*25.*x.^2/(150^2);
        Params.Sim.IC_Error = [0;0;0;0];
        Max_Pos_Err = 1.5;  %used for setting lqr gains
        Params.Ctl.Minimum_Goal_Distance = 1; %(m) how far to "look ahead"
        Simulation_Duration = 90;
        Control_Frequency = 50;   %Hz
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    case 5
        % Works
        Params.Traj.Func = true;
        x = linspace(0,150,Params.Traj.Npts+1)';
        Params.Traj.x = x;
        Params.Traj.y = sin(x/15).*25;
        Params.Sim.IC_Error = [0;0;0;0];
        Max_Pos_Err = 1.5;
        Params.Ctl.Minimum_Goal_Distance = 1; %(m) how far to "look ahead"
        Simulation_Duration = 160;
        Control_Frequency = 10;   %Hz

    case 6
      % USED IN PRESENTATION - Appendix 1
      %Long version of 5 (many turns)  -- good demonstration
        Params.Traj.Func = true;
        Params.Traj.Npts = 250; %Override the number of points
        Params.Traj.Max_Rd_Idx = 245;
        x = linspace(0,350,Params.Traj.Npts+1)';
        Params.Traj.x = x;
        Params.Traj.y = sin(x/20+pi/3).*cos(x/100)*25;
        Params.Sim.IC_Error = [0;0;0;0];
        Max_Pos_Err = 1.5;  %used for setting lqr gains
        Params.Ctl.Minimum_Goal_Distance = 1; %(m) how far to "look ahead"
        Simulation_Duration = 235;
        Control_Frequency = 50;   %Hz

    case 7
        %Cosine Road
        Params.Traj.Npts = 150;
        Freq = 0.35;
        Amp = 25;
        Params.Traj.Func = true;
        x = linspace(0,150,Params.Traj.Npts+1)';
        Params.Traj.x = x;
        Params.Traj.y = Amp*cos((2*pi/Amp)*x*Freq+0.1);
        Params.Sim.IC_Error = [0;0;0;0];
        Max_Pos_Err = 1.5;
        Params.Ctl.Minimum_Goal_Distance = 1; %(m) how far to "look ahead"
        Simulation_Duration = 250;
        Control_Frequency = 50;   %Hz

     case 8
      %REALLY LONG version of 5 (many turns)  -- good demonstration
        Params.Traj.Func = true;
        Params.Traj.Npts = 1000; %Override the number of points
        Params.Traj.Max_Rd_Idx = 995;
        x = linspace(0,350,Params.Traj.Npts+1)';
        Params.Traj.x = x;
        Params.Traj.y = sin(x/20+pi/3).*cos(x/100)*25;
        Params.Sim.IC_Error = [0;0;0;0];
        Max_Pos_Err = 1.5;  %used for setting lqr gains
        Params.Ctl.Minimum_Goal_Distance = 2; %(m) how far to "look ahead"
        Simulation_Duration = 235*4;
        Control_Frequency = 50;   %Hz

    otherwise
        error('Invalid Example Number')

end

author: MatthewPeterKelly

repo: dscTutorials


text:
% MSH_TO_VTK: Export to VTK format for plotting.
%
%  MSH_to_vtk (pts, values, filename, fieldnames)
%
% INPUT:
%
%     pts:        points at which the field was computed
%     values:     cell-array, with values of the fields at the selected points
%     filename:   name of the output file
%     fieldnames: how to name the saved variables in the vtk file
%
% OUTPUT:
%
%    a vtk structured mesh file named <filename> is produced
%
% Copyright (C) 2009, 2010 Carlo de Falco, Rafael Vazquez
% Copyright (C) 2011, 2015 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function msh_to_vtk (pts, values, filename, fieldnames)

  if (iscell (values))
    if (numel (values) ~= numel (fieldnames))
      error ('The number of fields and the number of names should be the same')
    end
  else
    values = {values};
    if (~iscell (fieldnames))
      fieldnames = {fieldnames};
    end
  end

  ndim = numel (size (pts)) - 1;
  rdim = size (pts, 1);

  str1 = cat (2,'<?xml version="1.0"?> \n', ...
'<VTKFile type="StructuredGrid" version="0.1"> \n', ...
'<StructuredGrid WholeExtent="0 %d 0 %d 0 %d"> \n', ...
'<Piece Extent="0 %d 0 %d 0 %d"> \n', ...
'<PointData>\n');

  str1b = cat (2, ...
'<DataArray type="Float32" Name="%s" format="ascii" NumberOfComponents="%d"> \n');

  str1c = cat (2,'</DataArray> \n');

  str2 = cat (2,'</PointData> \n', ...
'<Points> \n', ...
'<DataArray type="Float32" NumberOfComponents="3"> \n');

  str3 = cat (2, '\n', ...
'</DataArray>\n', ...
'</Points> \n', ...
'</Piece> \n', ...
'</StructuredGrid> \n', ...
'</VTKFile> \n');

% Even for 2D data, everything is saved in 3D
% ndims (or size) do not work properly for the 1D case. I remove singleton
% dimensions using this trick

  size_pts = size (pts);
  npts = size_pts (2:end);

  if (ndim < 3)
    npts (ndim+1:3) = 1;
  end
  if (rdim < 3)
    pts(ndim+1:3,:,:) = 0;
  end

  if (length (filename) < 4 || ~strcmp (filename(end-3:end), '.vts'))
    filename = cat (2, filename, '.vts');
  end

  fid = fopen (filename, 'w');
  if (fid < 0)
    error ('msh_to_vtk: could not open file %s', filename);
  end

  fprintf (fid, str1, ...
           npts(1)-1, npts(2)-1, npts(3)-1, ...
           npts(1)-1, npts(2)-1, npts(3)-1);

  for iopt = 1:numel(values)
    if (sum (size (values{iopt}) > 1) == ndim)
      ncomp = 1;
    elseif (sum (size (values{iopt}) > 1) == ndim + 1)
      ncomp = 3;
    else
      ncomp = 9;
    end
    if (ncomp == 3 && rdim < 3)
      values{iopt}(rdim+1:3,:,:) = 0;
    elseif (ncomp == 9 && rdim < 3)
      values{iopt}(rdim+1:3,rdim+1:3,:,:) = 0;
    end
    fprintf (fid, str1b, fieldnames{iopt}, ncomp);
    fprintf (fid, '%g ', values{iopt}(:));
    fprintf (fid, str1c);
  end
  fprintf (fid, str2);
  fprintf (fid, '%g ', pts(:));
  fprintf (fid, str3);

  fclose (fid);

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function k = rbfperiodicKernDiagCompute(kern, x)

% RBFPERIODICKERNDIAGCOMPUTE Compute diagonal of RBFPERIODIC kernel.
% FORMAT
% DESC computes the diagonal of the kernel matrix for the RBF derived periodic kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : a vector containing the diagonal of the kernel matrix
% computed at the given points.
%
% SEEALSO : rbfperiodicKernParamInit, kernDiagCompute, kernCreate, rbfperiodicKernCompute
%
% COPYRIGHT : Neil D. Lawrence, 2007

% KERN

k = repmat(kern.variance, size(x, 1), 1);


text:
function rowMatrix = computeRowsOfToeplitzHankelMatrix(rowNumber,...
        nColumns, crossCorrelationVectors, hankelMatrixIsAdded, dcIsIncluded)
    if rowNumber == 1
        toeplitzRows = crossCorrelationVectors(1:nColumns,:);
    else
        toeplitzRows = ...
            [flip(crossCorrelationVectors(2:rowNumber,:),1);...
            crossCorrelationVectors(1:nColumns-rowNumber+1,:)];
    end

    if dcIsIncluded && hankelMatrixIsAdded
        hankelOffset = 1;
    else
        hankelOffset = 3;
    end

    hankelRows = crossCorrelationVectors((0:nColumns-1)+...
                                         hankelOffset+rowNumber-1,:);

    if hankelMatrixIsAdded
        rowMatrix = toeplitzRows + hankelRows;
    else
        rowMatrix = toeplitzRows - hankelRows;
    end

end

author: LimingShi


text:
function stroud_test23 ( )

%*****************************************************************************80
%
%% TEST23 tests PARALLELIPIPED_VOLUME_ND.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST23\n' );
  fprintf ( 1, '  PARALLELIPIPED_VOLUME_ND computes the volume of a\n' );
  fprintf ( 1, '    parallelipiped in N dimensions.\n' );
  fprintf ( 1, '\n' );

  for n = 2 : 4

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Spatial dimension N = %d\n', n );
%
%  Set the values of the parallelipiped.
%
    v = setsim ( n );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Parallelipiped vertices:\n' );
    fprintf ( 1, '\n' );
    for i = 1 : n+1
      fprintf ( 1, '  ' );
      for j = 1 : n
        fprintf ( 1, '  %6f', v(i,j) );
      end
      fprintf ( 1, '\n' );
    end

    volume = parallelipiped_volume_nd ( n, v );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Volume is %f\n', volume );

  end

  return
end