text:
function [seg, t] = retroProjPlkEndPnts(Rob,Sen,Lmk,Obs)

% RETROPROJPLKENDPNTS  Retro project Plucker endpoints.
%   [SEG,T] = RETROPROJPLKENDPNTS(Rob,Sen,Lmk,Obs) retroprojects the
%   segment endpoints in Bos.meas.y onto the Plucker line Lmk.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

global Map

% Plucker line
l = Map.x(Lmk.state.r);

% rays in sensor frame
r1s = pix2PluckerRay(Sen.par.k,Obs.meas.y(1:2));
r2s = pix2PluckerRay(Sen.par.k,Obs.meas.y(3:4));

% rays in world frame
r1 = fromFramePlucker(Rob.frame,fromFramePlucker(Sen.frame,r1s));
r2 = fromFramePlucker(Rob.frame,fromFramePlucker(Sen.frame,r2s));

% endpoints and abscissas
[e1,t1] = intersectPlucker(l,r1);
[e2,t2] = intersectPlucker(l,r2);

% build segment and abscissas vector
if t2 > t1
    seg = [e1;e2];
    t   = [t1;t2];
else
    seg = [e2;e1];
    t   = [t2;t1];
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE


text:
function [mh,mS] = pred_coxphhs(gp, x, y, xt, varargin)
%PRED_COXPHHS  Return hazard and survival functions
%
%  Description
%    [H,S] = PRED_COXPHHS(GP,X,Y,XT)
%    Returns hazard and survival functions for the inputs XT. GP is the
%    Gaussian process structure, X and Y are the training inputs and
%    outputs.
%
% Copyright (c) 2012-2013 Ville Tolvanen, Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

ntime=size(gp.lik.xtime,1);
if iscell(gp)
  % prediction for GP_IA cell array
  nGP = numel(gp);
  mh=zeros(size(xt,1),ntime,nGP);
  mS=zeros(size(xt,1),ntime,nGP);
  P_TH=zeros(1,nGP);
  for i1=1:nGP
    % make prediction for each gp in cell array
    Gp=gp{i1};
    P_TH(:,i1) = Gp.ia_weight;
    [mh(:,:,i1), mS(:,:,i1)]=pred_coxphhs(Gp, x, y, xt, varargin{:});
    mh(:,:,i1)=P_TH(:,i1)*mh(:,:,i1);
    mS(:,:,i1)=P_TH(:,i1)*mS(:,:,i1);
  end
  % combine predictions
  mh=sum(mh,3);
  mS=sum(mS,3);
  return
elseif numel(gp.jitterSigma2)>1
  nmc=size(gp.jitterSigma2,1);
  mh=zeros(size(xt,1),ntime,nmc);
  mS=zeros(size(xt,1),ntime,nmc);
  for i1=1:nmc
    Gp = take_nth(gp,i1);
    [mh(:,:,i1), mS(:,:,i1)]=pred_coxphhs(Gp, x, y, xt, varargin{:});
  end
  % combine predictions
  mh=mean(mh,3);
  mS=mean(mS,3);
  return
end

[Ef1, Ef2, Covf] = pred_coxph(gp,x,y,xt, varargin{:});
nsamps = 10000;
if isfield(gp.lik, 'stratificationVariables')
  ind_str=gp.lik.stratificationVariables;
  ux=unique([x(:,ind_str); xt(:,ind_str)],'rows');
  nu=size(ux,1);
  for i1=1:size(ux,1)
    uind{i1}=find(xt(:,ind_str)==ux(i1,:));
  end
  nf1=ntime*nu;
else
  nf1=ntime;
end
sd=gp.lik.stime(2)-gp.lik.stime(1);

Sigm_tmp=Covf;
Sigm_tmp=(Sigm_tmp+Sigm_tmp')./2;
% f_star=mvnrnd(Ef1, Sigm_tmp(1:ntime,1:ntime), nsamps);
f_star=mvnrnd([Ef1;Ef2], Sigm_tmp, nsamps);

f1=f_star(:,1:nf1);
f2=f_star(:,(nf1+1):end);

la1=exp(f1);
eta2=exp(f2);

if ~isfield(gp.lik, 'stratificationVariables')
  hb=(la1'*sd);
  cumsumtmp=cumsum(hb)';
  for i1=1:size(eta2,2)
    mh(i1,:)=mean(bsxfun(@times,hb',eta2(:,i1)));
    Stime=exp(-bsxfun(@times,cumsumtmp,eta2(:,i1)));
    mS(i1,:)=mean(Stime);
  end
else
  for i2=1:length(uind)
    hb=(la1(:,(i2-1)*ntime+1:i2*ntime)'*sd);
    cumsumtmp=cumsum(hb)';
    for i1=1:size(uind{i2},1)
      ind=uind{i2}(i1);
      mh(ind,:)=mean(bsxfun(@times,hb',eta2(:,ind)));
      Stime=exp(-bsxfun(@times,cumsumtmp,eta2(:,ind)));
      mS(ind,:)=mean(Stime);
    end
  end
end

end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function out = cfg_run_dir_move(job)

% Move, copy or delete directory
%
% 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: cfg_run_dir_move.m 6920 2016-11-02 14:48:01Z guillaume $

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

action = fieldnames(job.action);
action = action{1};
% delete
if strcmp(action, 'delete')
    rmdir(char(job.dir),'s');
    out = [];
else
    [dummy, dirname] = fileparts(char(job.dir));
    if strncmp(action, 'copy', 4)
        % copy
        cmd = @copyfile;
        if strcmp(action,'copyto')
            tgt = job.action.copyto{1};
        else
            tgt = job.action.(action).copyto{1};
        end
        tgt = fullfile(tgt, dirname);
        out.dir = {tgt};
    else
        % move
        cmd = @movefile;
        if strcmp(action,'moveto')
            tgt = job.action.moveto{1};
        else
            tgt = job.action.(action).moveto{1};
        end
        if exist(tgt, 'dir')
            out.dir = {fullfile(tgt, dirname)};
        else
            out.dir = tgt;
        end
    end

    feval(cmd, char(job.dir), tgt);
end

author: spm


text:
% By Roger Aarenstrup, roger.aarenstrup@mathworks.com
% 2006-08-18
%
% This is a state-space DC motor model of a
% Maxon RE25 10 Watt, precious metal brushes, 118743
%
% This model also have a weak connection to a load.
%

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%            The full dc motor model                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Vin is the input voltage to the motor
% i is the motor current
% th_m is the rotor angle, theta
% dth_m is the rotor angular velocity sometimes called omega
% th_l is the load angle
% dth_l is the load angular velocity

% Controller Sample Time
Ts = 100e-6;

% PARAMETERS DC MOTOR
Rm  = 2.06;             % Motor resistance (ohm)
Lm  = 0.000238;         % motor inductance (Henrys)
Kb = 1/((406*2*pi)/60); % Back EMF constant (Volt-sec/Rad)
Kt = 0.0235;            % Torque constand (Nm/A)
Jm = 1.07e-6;           % Rotor inertia (Kg m^2)
bm = 12e-7;             % MEchanical damping (linear model of
                        % friction: bm * dth)
% PARAMETERS LOAD
Jl = 10.07e-6;   % Load inertia (10 times the rotor)
bl = 12e-6;      % Load damping (friction)
Ks = 100;        % Spring constant for connection rotor/load
b = 0.0001;      % Spring damping for connection rotor/load

% SYSTEM MATRICES
%
% States:  [i dth_m th_m dth_l th_l]'
% Input:   Vin the motor voltage
% Outputs: same as states
%
Afull = [-Rm/Lm    -Kb/Lm       0       0       0;
         Kt/Jm   -(bm+b)/Jm   -Ks/Jm    b/Jm     Ks/Jm;
           0          1         0       0        0;
           0         b/Jl    Ks/Jl   -(b+bl)/Jl    -Ks/Jl;
           0          0         0       1        0];

Bfull = [1/Lm 0 0 0 0]';

Cfull = [0 1 0 0 0;
         0 0 1 0 0;
         0 0 0 1 0;
         0 0 0 0 1];

Dfull = [0 0 0 0]';

sys_full = ss(Afull, Bfull, Cfull, Dfull);

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     The reduced dc motor model for current control      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% SYSTEM MATRICES
%
% States:  [dth_m th_m dth_l th_l]'
% Input:   I The current to the dc motor
% Outputs: same as states
%
Ared = [ -(bm+b)/Jm   -Ks/Jm    b/Jm     Ks/Jm;
              1         0       0        0;
             b/Jl    Ks/Jl   -(b+bl)/Jl    -Ks/Jl;
              0         0       1        0];

Bred = [Kt/Jm 0 0 0]';

Cred = eye(4);

Dred = [0 0 0 0]';

sys_red = ss(Ared, Bred, Cred, Dred);

% Discrete version of the model (as seen from the controller)
sys_red_d = c2d(sys_red, Ts, 'zoh');

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         Scaling of reduced state-space model            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% u = Nu * un
% x = Nx * xn
% y = Ny * yn
% u, x, y are real inputs, states and outputs
% un, xn, yn are normalized

max_current = 0.01; % 10 mA as input
max_pos = 4.1;
max_v = 4.1;

Nu = max_current;
Nx = [max_v 0 0 0; 0 max_pos 0 0; 0 0 max_v 0; 0 0 0 max_pos];
Ny = Nx;

An = inv(Nx)*Ared*Nx;
Bn = inv(Nx)*Bred*Nu;
Cn = inv(Ny)*Cred*Nx;
Dn = 0;

sys_red_n = ss(An, Bn, Cn, Dn);

sys_red_d_n = c2d(sys_red_n, Ts, 'zoh');

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        State-space controller design attempt 1          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The discrete poles of the closed loop system
dpole1 = exp(-2*pi*Ts*[20 22 24 26]);

% Calculate the control parameters
L1 = place(sys_red_d_n.a, sys_red_d_n.b, dpole1);

% Keep the static gain of the closed loop system to 1
F=feedback(sys_red_d_n, L1);   % The closed loop system, F.
Kstatic = freqresp(F, 0);    % Get the static gain of F.
Kstat = 1/Kstatic(4);        % It is the fourth output (load position)

% The above commands requires toolboxes, incase you don't have them
% you can simulate the system with a unit step is see the output static
% gain. Kstat will be the inverse of that.

% to verify the pole placement
%pzmap(F);

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  State-space controller design attempt 2 - Integrator   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The discrete poles of the closed loop system
dpole2 = exp(-2*pi*Ts*[20 22 24 26 5]);

% We only consider one output here
Cred_one = [0 1 0 0];

% Calculate the control parameters
Aint = [sys_red_d_n.a [0 0 0 0]';
        -Cred_one 1];

Bint= [sys_red_d_n.b' 0]';

Cint = [Cred_one 0];

Dint = [0]';

sys_int_d2 = ss(Aint, Bint, Cint, Dint, Ts);

% Controllability VS Observability analysis
%ob = obsv(sys_int_d2.a, sys_int_d2.c);
%cr = ctrb(sys_int_d2.a, sys_int_d2.b);
%rank(ob)
%rank(cr)

[L2, a, m] = place(sys_int_d2.a, sys_int_d2.b, dpole2);

% Feed Forward gain
Kff = L2(5)/(dpole2(5) - 1);

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  State-space controller design attempt 3 - Observer     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% These should be about twice as fast as the state feedback
% for the controller.
dpole3 = exp(-2*pi*Ts*[1400 1450 1500 1550]);
%dpole3 = [-0.3 -0.33 -0.35 -0.37];

% Recalculate the reduced model with one ouput

Ared2 = inv(Nx)*Ared*Nx; % Scaling
Bred2 = inv(Nx)*Bred*Nu; % Scaling
Cred2 = [0 1 0 0];
Dred2 = 0;

sys_red2 = ss(Ared2, Bred2, Cred2, Dred2);

% Discrete version of the model (as seen from the controller)
sys_red_d2 = c2d(sys_red2, Ts, 'zoh');

set(sys_red_d2, 'inputname', {'Um'}, ...
                'outputname', {'Enc_out'});

% Calculate the observer state feedback
K = place(sys_red_d2.a', sys_red_d2.c', dpole3);

% Put the observer together
Aobs = [sys_red_d2.a-K'*sys_red_d2.c];
Bobs = [sys_red_d2.b K'];
Cobs = eye(4);
Dobs = [0 0 0 0; 0 0 0 0]';

observer = ss(Aobs, Bobs, Cobs, Dobs, Ts, 'inputname', {'Ue', 'Enc_in'}, ...
              'outputname', {'dth_m', 'th_m', 'dth_l', 'th_l'});

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  State-space controller design attempt 4 - The Servo Case  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Kinv = (Jl+Jm)/Kt;

dpole4 = exp(-2*pi*Ts*[200 204 220 224 300]);
[L3, a, m] = place(sys_int_d2.a, sys_int_d2.b, dpole4);

author: Sable

repo: mcbench-benchmarks


text:
function fnames = subs_fields(item)

% function fnames = subs_fields(item)
% This function works as a "class-based switch" to return the value of
% the private mysubs_fields function for the appropriate class.
% This function is identical for all classes derived from cfg_item, but
% it has to be in the class directory to access the proper private
% function mysubs_fields.
%
% 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: subs_fields.m 4867 2012-08-30 13:04:51Z volkmar $

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

fnames = mysubs_fields;


text:
function A = countTransitions(data, nstates)
%% Count the discrete transitions that occur in data
% data is a cell array of cases of potentially different lengths.
%%

% This file is from pmtk3.googlecode.com

if ~iscell(data)
   if isvector(data)
       data = {data};
   else
       data = mat2cellRows(data);
   end
end

A = zeros(nstates, nstates);
for i=1:numel(data)
   obs = colvec(data{i});
   A = A + accumarray([obs(1:end-1), obs(2:end)], 1, [nstates, nstates]);
end

end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function test_suite = test_randomize_targets()
% tests for cosmo_randomize_targets
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_randomize_targets_basics()
    ds=cosmo_synthetic_dataset('ntargets',4,'nchunks',10);
    [x1,perm1]=cosmo_randomize_targets(ds);
    assertEqual(size(ds.sa.targets),[40 1])
    assertEqual(x1,ds.sa.targets(perm1))
    assertEqual(sort(perm1),(1:size(ds.samples,1))');

    x2=cosmo_randomize_targets(ds);
    assert(any(x1~=x2)); % probablity of failing less than 1e-13
    assert(any(x1~=ds.sa.targets));

    ds_small=cosmo_slice(ds,1:8);
    x3=cosmo_randomize_targets(ds_small,'seed',1);
    assertEqual(x3,[ 3 4 1 2 2 1 3 4]');

    x4=cosmo_randomize_targets(ds_small,'seed',314);
    assertEqual(x4,[ 3 2 4 1 4 3 2 1]');

    ds_single_target=cosmo_slice(ds,1:4:20);
    x5=cosmo_randomize_targets(ds_single_target,'seed',314);
    assertEqual(x5,ones(5,1));

    ds_between=cosmo_slice(ds,1:5:40);
    x6=cosmo_randomize_targets(ds_between,'seed',1);
    assertEqual(x6,[3 2 1 3 4 4 1 2]');
    x7=cosmo_randomize_targets(ds_between);
    assertEqual(histc(x7,1:4),[2 2 2 2]');

    % test exceptions
    aet=@(varargin)assertExceptionThrown(@()cosmo_randomize_targets(...
    varargin{:}),'');

    ds_missing=cosmo_slice(ds,1:3:20);
    aet(ds_missing);

    ds.sa=rmfield(ds.sa,'targets');
    aet(ds);

    ds=rmfield(ds,'sa');
    aet(ds);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function y = logn_lpdf(x,mu,sigma)
%LOGNPDF Lognormal log-probability density function (lpdf).
%
%   Y = LOGN_LPDF(X,MU,SIGMA) Returns the log of the lognormal pdf at
%   the values in X. The mean and standard deviation of log(Y) are
%   MU and SIGMA.
%
%   The size of Y is the common size of the input arguments. A scalar input
%   functions as a constant matrix of the same size as the other inputs.
%
%   Default values for MU and SIGMA are 0 and 1 respectively.

% Copyright (c) 2000 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

y = -0.5 * ((log(x) - mu)./sigma).^2 - log(x .* sqrt(2*pi) .* sigma);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function p45=p45(h)
sldata;
limitpoints;
slopes;
T45=2.78e-3*h+139.1;
p2=p1*(T2/T1)^(-g/(m12*R));
p3=p2*exp(-g*(h3-h2)/(R*T2));
p4=p3*(T4/T3)^(-g/(m34*R));
p45=p4*(T45/T4)^(-g/(m45*R));


text:
 function [xs, Ps] = qpwls_qn(x, G, W, yy, C, P, niter)
%function [xs, Ps] = qpwls_qn(x, G, W, yy, C, P, niter)
%
% quadratic penalized weighted least squares (QPWLS) via
% (preconditioned) quasi-Newton (QN) algorithm
% cost(x) = (y-Gx)'W(y-Gx) / 2  + x'C'Cx / 2
% in
%	x	[np,1]		initial estimate
%	G	[nn,np]		system matrix
%	W	[nn,nn]		data weighting matrix
%	yy	[nn,1]		noisy data
%	C	[nc,np]		penalty 'derivatives' (R = \Half C'*C)
%	P	[np,np]		preconditioner (matrix or object)
%	niter			# total iterations
% out
%	xs	[np,niter]	estimates each iteration
%	ws	[np,niter]	P rank-one update terms each iteration
%
% Copyright July 2000, Jeff Fessler, The University of Michigan

if nargin < 3 || nargin > 7, ir_usage, end
np = ncol(G);

if ~isvar('C')		| isempty(C),		C = 0;			end
if ~isvar('P')		| isempty(P),		P = speye(np);		end
if ~isvar('niter')	| isempty(niter),	niter = 2;		end
if ~isvar('x')		| isempty(x),		x = zeros(np,1);	end

yy = yy(:);
xs = zeros(np, niter);
xs(:,1) = x;

ws = zeros(np, niter);
bs = zeros(1, niter);

%
% initialize projections
%
Gx = G * x;
Cx = C * x;

grad_next = G' * (W * (yy - Gx)) - C' * Cx;

%
% iterate
%
for iter = 2:niter

	grad = grad_next;

	%
	% compute ddir = P * grad, which is -s
	%
	ddir = Pmul(grad, P, bs(1:(iter-2)), ws(:,1:(iter-2)));
	x = x + ddir;

	% check if descent direction
	if ddir' * grad < 0
		warning('wrong direction')
%		keyboard
	end

	%
	% next (negative) gradient
	%
	Gx = G * x;
	Cx = C * x;
	grad_next = G' * (W * (yy - Gx)) - C' * Cx;

	q = grad_next - grad;
	w = ddir + Pmul(q, P, bs(1:(iter-2)), ws(:,1:(iter-2)));

	tmp = w' * q;
	if tmp == 0, error zero, end
	bs(iter-1) = 1 / tmp;
	ws(:,iter-1) = w;

	%
	% step size in search direction
	%
	if 0
		Gdir = G * ddir;
		Cdir = C * ddir;

		step = (ddir' * grad) / (Gdir'*(W*Gdir) + Cdir'*Cdir);
		steps(iter,1) = step;
		if step < 0
			warning('downhill?')
			keyboard
		end

		%
		% update
		%
		x	= x + step * ddir;
		Gx	= Gx  + step * Gdir;
		Cx	= Cx  + step * Cdir;
	end

	xs(:,iter) = x;
end

%
% multiply v = P_n * u where P_n = P_0 + \sum_i b_i w_i w_i'
%
function v = Pmul(u, P, bs, ws)
v = P * u;
bs = bs .* (u' * ws);
v = v - ws * bs';

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%
% Computes manipulability of poses given in qs
%
function manips = compute_manip(robot, qs)
manips = [];
for i=1:size(qs,2)
    %compute current manip
    J = manipulator_jacobian(robot, qs(:,i));
    manip = sqrt(det(J*J'));
    manips = [manips manip];
end

repo: ARTE


text:
function mono_between_random_test ( )

%*****************************************************************************80
%
%% MONO_BETWEEN_RANDOM_TEST tests MONO_BETWEEN_RANDOM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 November 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MONO_BETWEEN_RANDOM_TEST\n' );
  fprintf ( 1, '  MONO_BETWEEN_RANDOM selects at random a monomial\n' );
  fprintf ( 1, '  in M dimensions of total degree between N1 and N2.\n' );

  m = 3;
  n1 = 2;
  n2 = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Let M =  %d\n', m );
  fprintf ( 1, '      N1 = %d\n', n1 );
  fprintf ( 1, '      N2 = %d\n', n2 );
  fprintf ( 1, '\n' );

  seed = 123456789;
  test_num = 5;

  for test = 1 : test_num

    [ x, rank, seed ] = mono_between_random ( m, n1, n2, seed );

    fprintf ( 1, '  %3d:', rank );
    for j = 1 : m
      fprintf ( 1, '  %1d', x(j) );
    end
    fprintf ( 1, '\n' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Z = mldivide(X,Y)
%MLDIVIDE Slash left division for tensors.
%
%   MLDIVIDE(A,B) is called for the syntax 'A \ B' when A is a scalar and B
%   is a tensor.
%
%   Example
%   X = tenrand([4 3 2],5);
%   3 \ X
%
%   See also TENSOR, TENSOR/LDIVIDE.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

if isscalar(X)
    Z = tenfun(@ldivide,X,Y);
    return;
end

error('MLDIVIDE only supports the scalar case for tensors');

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function pass = test_chebfun_lu( pref )
% Test Chebfun LU command.

if ( nargin == 0)
    pref = chebfunpref;
end

tol = pref.chebfuneps;
j = 1;

% Check accuracy on [-1,1]
x = chebfun(@(x) x);
A = [1 x x.^2 x.^3 x.^4 x.^5];
[L, U, p] = lu( A );

pass(j) = norm( triu(U) - U ) < tol; j = j + 1;
pass(j) = norm( A - L * U ) < tol; j = j + 1;
pass(j) = norm( diag( L(p,:) ) - ones(size(L,2),1) ) < 10*tol; j = j + 1;
pass(j) = norm( tril(L(p,:)) - L(p,:)) < 10*tol; j = j + 1;

% Check accuracy on [-2,3]
x = chebfun(@(x) x, [-2 3]);
A = [1 x x.^2 x.^3 x.^4 x.^5];
[L, U, p] = lu( A );

pass(j) = norm( triu(U) - U ) < tol; j = j + 1;
pass(j) = norm( A - L * U ) < 1e3*tol; j = j + 1;

pass(j) = norm( diag( L(p,:) ) - ones(size(L,2),1) ) < 10*tol; j = j + 1;

pass(j) = norm( tril(L(p,:)) - L(p,:)) < 20*tol; j = j + 1;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [C_trialAvr,C_trialRes,C_score,C_d2var_perstim,C_p] = GetTrialAvrLongTrace(hfig,C)
%%

isFulllength = 0; % isFulllength = true for analyses before 12/18/17 - XC

fishset = getappdata(hfig,'fishset');
periods = getappdata(hfig,'periods');

if fishset == 1
    period = periods;
    C_3D_0 = reshape(C,size(C,1),period,[]);

%     C_3D = zscore(C_3D_0,0,2);
%     C_d2var_perstim = nanmean(nanstd(C_3D,0,3),2);
%     C_score = C_d2var_perstim;

    C_period = mean(C_3D_0,3);%prctile(C_3D_0,20,3);%mean(C_3D_0,3);
    nPeriods = round(size(C,2)/period);
    C_trialAvr = repmat(C_period,1,nPeriods);

    C_trialRes = C-C_trialAvr;
    C_3D_tRes = reshape(C_trialRes,size(C,1),period,[]);
%     C_3D_tRes = zscore(C_3D_tRes_0,0,2);
    d2var = (nanstd(C_3D_tRes,0,3)).^2;
    C_score = nanmean(d2var,2);

    C_p = C;
    C_d2var_perstim = [];
else
    stimset = getappdata(hfig,'stimset');
    stimrange = getappdata(hfig,'stimrange');
    timelists = getappdata(hfig,'timelists');

    C_trialAvr = [];
%     C_d2std = [];
    C_d2var_perstim = [];
    C_p = []; % 12/18/17 - XC
    for i = 1:length(stimrange)
        i_stim = stimrange(i);
        if stimset(i_stim).nReps(1)>1
            offset = length(horzcat(timelists{stimrange(1:i-1)}));% works for i=0 too
            tIX_ = 1+offset:length(timelists{stimrange(i)})+offset;

            period = periods(i_stim);
            C_this = C(:,tIX_);

            C_this = zscore(C_this,0,2); % 05/11/18

            C_3D_0 = reshape(C_this,size(C,1),period,[]);
            C_period = mean(C_3D_0,3);%median(C_3D_0,3);
            nPeriods = length(tIX_)/period;
            C_trialAvr_this = repmat(C_period,1,nPeriods);
            C_trialAvr = horzcat(C_trialAvr,C_trialAvr_this); %#ok<AGROW>
            C_trialRes_this = C_this-C_trialAvr_this;

            C_3D_tRes = reshape(C_trialRes_this,size(C,1),period,[]);
            d2var = (nanstd(C_3D_tRes,0,3)).^2; % d2var = (mad(C_3D_tRes,0,3)).^2;
            C_d2var_perstim = horzcat(C_d2var_perstim,nanmean(d2var,2)); %#ok<AGROW>

            C_p = horzcat(C_p,C_this); % 12/18/17 - XC

%             C_3D = zscore(C_3D_0,0,2);
% %             C_d2std = horzcat(C_d2std,nanstd(C_3D,0,3));
%             d2var = (nanstd(C_3D,0,3)).^2;
%             C_d2var_perstim = horzcat(C_d2var_perstim,nanmean(d2var,2));
%             %             C_3D = zscore(C_3D_0,0,2);
%             %             H_raw = horzcat(H_raw,nanmean(nanstd(C_3D,0,3),2));
        else % taking out 12/18/17 - XC
            if isFulllength
                offset = length(horzcat(timelists{stimrange(1:i-1)}));% works for i=0 too
                tIX_ = 1+offset:length(timelists{stimrange(i)})+offset;
                C_trialAvr = horzcat(C_trialAvr,zeros(size(C,1),length(tIX_))); %#ok<AGROW>
            end
        end
    end

    % condense C_d2var_perstim for multiple stimsets down to 1 score
    if isempty(C_d2var_perstim)
        errordlg('chosen stimulus range not suitable for stim-lock analysis');
%         C_score = 1:numU;
        return;
    end
    C_score = zeros(size(C_d2var_perstim,1),1);
    for i = 1:size(C_d2var_perstim,1)
        C_score(i) = min(C_d2var_perstim(i,:));
    end

    if isFulllength
        C_trialRes = C-C_trialAvr;
    else
        C_trialRes = C_p-C_trialAvr; % 12/18/17 - XC
    end
end

end

author: xiuyechen

repo: FishExplorer


text:
function fx = p04_fun ( option, nvar, x )

%*****************************************************************************80
%
%% P04_FUN evaluates the function for problem 4.
%
%  Title:
%
%    The Broyden function
%
%  Description:
%
%    The function F is derived via homotopy from a simpler function G:
%
%      F(X(1),X(2),X(3)) = g(X(1),X(2)) + (X(3)-1) * G(Y1,Y2).
%
%    with
%
%      (Y1,Y2) some starting point,
%
%    and
%
%      G(1) = 0.5*sin(X(1)*X(2)) - X(2)/PI - X(1)
%      G(2) = (1-1/(4*PI))*(exp(2*X(1))-E) + E*X(2)/PI- 2*E*X(1)
%
%    where "E" = exp(1).
%
%  Options:
%
%    The only option starts with (0.4, 3, 0), and seeks the target
%    solution whose third component is 1.  The correct value of the
%    target solution is
%
%      ( -0.2207014, 0.8207467, 1.0D+00 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Charles Broyden,
%    A New Method of Solving Nonlinear Simultaneous Equations,
%    The Computer Journal,
%    Volume 12, 1969, pages 94-99.
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Input, real X(NVAR), the argument of the function.
%
%    Output, real FX(NVAR-1), the value of the function at X.
%

%
%  Get the starting point.
%
  y = p04_start ( option, nvar );
%
%  Get the function value at the starting point and at the
%  current point.
%
  gy = p04_gx ( y );
  gx = p04_gx ( x );
%
%  Use X3 to compute a homotopy.
%
  for i = 1 : nvar-1
    fx(i) = gx(i) + ( x(3) - 1.0 ) * gy(i);
  end

  return
end


text:
function convertLeksellFileToHex(fileName)
%"convertLeksellFileToHex"
%   A simple file reader that reads a binary file and writes it out in
%   hexidecimal for easy decription.
%
%KRK 05/31/07
%
%Usage:
%   convertLeksellFileToHex(fileName)
%
% 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/>.

% Open the given file and read the binary data in bytes
fid = fopen(fileName, 'r', 'b');
binfile = fread(fid, inf, 'uint8');
fclose(fid);

hexfile = dec2hex(binfile, 2); % Convert the bin array to 2 digit hex

wfid = fopen(strcat(fileName,'_HEX.txt'), 'w');
num = size(hexfile); % Gets the size of the hex array
% Write it to a file
for i = 1:1:num(1)
    fwrite(wfid, hexfile(i,:));

    % Space the file out
    if mod(i,10)
        fprintf(wfid, ' ');
    else
        fprintf(wfid, '\n');
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ centroids ] = Kmeans_update( data, centroids, assignments )
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here

for i = 1:size(centroids,3);
    pos = find(assignments == i);
    newCenter = zeros(size(centroids,1), size(centroids,2));
    for j = pos
        newCenter(:,:) = newCenter(:,:) + data(:,:,j);
    end
    newCenter = newCenter ./ length(pos);
    centroids(:,:,i) = newCenter;
end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [ indx, i, j ] = sort_heap_external ( n, indx, isgn )

%*****************************************************************************80
%
%% SORT_HEAP_EXTERNAL externally sorts a list of items into ascending order.
%
%  Discussion:
%
%    The actual list of data is not passed to the routine.  Hence this
%    routine may be used to sort integers, reals, numbers, names,
%    dates, shoe sizes, and so on.  After each call, the routine asks
%    the user to compare or interchange two items, until a special
%    return value signals that the sorting is completed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Albert Nijenhuis, Herbert Wilf.
%    Combinatorial Algorithms,
%    Academic Press, 1978, second edition,
%    ISBN 0-12-519260-6.
%
%  Parameters:
%
%    Input, integer N, the number of items to be sorted.
%
%    Input, integer INDX, the main communication signal.
%    The user must set INDX to 0 before the first call.
%    Thereafter, the user should set the input value of INDX
%    to the output value from the previous call.
%
%    Input, integer ISGN, results of comparison of elements I and J.
%    (Used only when the previous call returned INDX less than 0).
%    ISGN <= 0 means I is less than or equal to J;
%    0 <= ISGN means I is greater than or equal to J.
%
%    Output, integer INDX, the main communication signal.
%    If INDX is
%
%      greater than 0, the user should:
%      * interchange items I and J;
%      * call again.
%
%      less than 0, the user should:
%      * compare items I and J;
%      * set ISGN = -1 if I < J, ISGN = +1 if J < I;
%      * call again.
%
%      equal to 0, the sorting is done.
%
%    Output, integer I, J, the indices of two items.
%    On return with INDX positive, elements I and J should be interchanged.
%    On return with INDX negative, elements I and J should be compared, and
%    the result reported in ISGN on the next call.
%
  persistent i_save;
  persistent j_save;
  persistent k;
  persistent k1;
  persistent n1;

  if ( isempty ( i_save ) )
    i_save = -1;
  end

  if ( isempty ( j_save ) )
    j_save = -1;
  end
%
%  INDX = 0: This is the first call.
%
  if ( indx == 0 )

    k = floor ( n / 2 );
    k1 = k;
    n1 = n;
%
%  INDX < 0: The user is returning the results of a comparison.
%
  elseif ( indx < 0 )

    if ( indx == -2 )

      if ( isgn < 0 )
        i_save = i_save + 1;
      end

      j_save = k1;
      k1 = i_save;
      indx = -1;
      i = i_save;
      j = j_save;
      return;
    end

    if ( 0 < isgn )
      indx = 2;
      i = i_save;
      j = j_save;
      return;
    end

    if ( k <= 1 )

      if ( n1 == 1 )
        i_save = 0;
        j_save = 0;
        indx = 0;
      else
        i_save = n1;
        n1 = n1 - 1;
        j_save = 1;
        indx = 1;
      end

      i = i_save;
      j = j_save;
      return;

    end

    k = k - 1;
    k1 = k;
%
%  0 < INDX, the user was asked to make an interchange.
%
  elseif ( indx == 1 )

    k1 = k;

  end

  while ( 1 )

    i_save = 2 * k1;

    if ( i_save == n1 )
      j_save = k1;
      k1 = i_save;
      indx = -1;
      i = i_save;
      j = j_save;
      return;
    elseif ( i_save <= n1 )
      j_save = i_save + 1;
      indx = -2;
      i = i_save;
      j = j_save;
      return;
    end

    if ( k <= 1 )
      break;
    end

    k = k - 1;
    k1 = k;

  end

  if ( n1 == 1 )
    i_save = 0;
    j_save = 0;
    indx = 0;
    i = i_save;
    j = j_save;
  else
    i_save = n1;
    n1 = n1 - 1;
    j_save = 1;
    indx = 1;
    i = i_save;
    j = j_save;
  end

  return
end


text:
function [err, indx, SubLabel] = WhichSubBricks (InfoStat, Type, Label)
%
%   [err, indx, SubLabel] = WhichSubBricks (InfoStat, Type, Label)
%
%Purpose:
%
%
%
%Input Parameters:
%
%
%
%Output Parameters:
%   err : 0 No Problem
%       : 1  Problems
%   indx: a vector containing sub-brick indices (+1 from AFNI indexing)
%         with labels containing both of Type and Label strings
%   SubLabel: a vector of structures containing the label strings.
%         For example if indx = [1 4] then the labels you were searching
%         for are SubLabel(1).str and SubLabel(4).str which correspond to
%         AFNI sub-bricks 0 and 3 in AFNI's indexing world.
%
%
%Key Terms:
%
%More Info :
%
%
%
%
%     Author : Ziad Saad
%     Date : Fri Jul 18 16:56:58 EDT 2003
%     SSCC/NIMH/ National Institutes of Health, Bethesda Maryland

%Define the function name for easy referencing
FuncName = 'WhichSubBricks';

%Debug Flag
DBG = 1;

%initailize return variables
err = 1;

indx = zeros(1,InfoStat.DATASET_RANK(2));

%change the format of BRICK_LABS
cend = 0;
for (i=1:1:InfoStat.DATASET_RANK(2)),
   [err, SubLabel(i).str, cend] = NextString(InfoStat.BRICK_LABS,'~',cend+1);
end

cnt = 0;
for (i=1:1:InfoStat.DATASET_RANK(2)),
   if (~isempty(findstr(SubLabel(i).str, Label)) & ~isempty(findstr(SubLabel(i).str, Type))),
      cnt = cnt + 1;
      indx(cnt) = i;
   end
end

if (cnt == 0), indx = [];
else indx = indx(1:cnt);
end

err = 0;
return;


text:
function [results] = hyperGlrt(M, t)
% HYPERGLRT Performs the generalized liklihood test ratio algorithm
%   Performs the generalized liklihood test ratio algorithm for target
% detection.
%
% Usage
%   [results] = hyperGlrt(M, U, target)
% Inputs
%   M - 2d matrix of HSI data (p x N)
%   t - target of interest (p x 1)
% Outputs
%   results - vector of detector output (N x 1)
%
% References
%   T F AyouB, "Modified GLRT Signal Detection Algorithm," IEEE
% Transactions on Aerospace and Electronic Systems, Vol 36, No 3, July
% 2000.

[p, N] = size(M);

% Remove mean from data
u = mean(M.').';
M = M - repmat(u, 1, N);
t = t - u;

R = inv(hyperCov(M));

results = zeros(1, N);
for k=1:N
    x = M(:,k);
    results(k) = ((t'*R*x)^2) / ((t'*R*t)*(1 + x'*R*x));
end

author: davidkun

repo: HyperSpectralToolbox


text:
function values = fast_intersection(rows, cols, values, target, unknown_dist, far_dist)
%FAST_INTERSECTION Under the assumption of categorical distance for the
% intersecting simplicial set perform a fast intersection.
%
% values = FAST_INTERSECTION(rows, cols, values, target, unknown_dist, far_dist)
%
% Parameters
% ----------
% rows: array
%     An array of the row of each non-zero in the sparse matrix
%     representation.
%
% cols: array
%     An array of the column of each non-zero in the sparse matrix
%     representation.
%
% values: array
%     An array of the value of each non-zero in the sparse matrix
%     representation.
%
% target: array of shape (n_samples, 1)
%     The categorical labels to use in the intersection.
%
% unknown_dist: double (optional, default 1)
%     The distance an unknown label (-1) is assumed to be from any point.
%
% far_dist: double (optional, default 5)
%     The distance between unmatched labels.
%
% Returns
% -------
% values: array
%     The non-zero entries resulting from the intersection.
%
%   AUTHORSHIP
%   Math Lead & Primary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Secondary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause

    if nargin < 6
        far_dist = 5;
        if nargin < 5
            unknown_dist = 1;
        end
    end

    row_labels = target(rows);
    col_labels = target(cols);

    unknown = row_labels == -1 | col_labels == -1;
    far = (row_labels ~= col_labels) & ~unknown;

    values(unknown) = values(unknown)*exp(-unknown_dist);
    values(far) = values(far)*exp(-far_dist);

    end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function y = cvx_s_diagonal( m, n )
%CVX_S_DIAGONAL Diagonal matrices.
y = cvx_s_banded( m, n, 0, 0 );

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

repo: REST


text:
function v = zvector(varargin)
% vector (0,0,1)

v = vector3d(0,0,1,varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
function urdf = urdfparse(filename)

    root = xmlread(filename);

    robot = root.getElementsByTagName('robot').item(0);

    %% process all the materials elements
    % - these are optional but we do them first because link elements
    %   reference them
    materialNodes = getChildrenByTagName(robot, 'material');
    materials = [];
    materialNames = [];
    if ~isempty(materialNodes)

        for i = 1:length(materialNodes)
            node = materialNodes(i);
            materials(i).name = string(node.getAttribute('name'));
            materials(i).rgba = str2num(node.getElementsByTagName('color').item(0).getAttribute('rgba'));
            materials(i).node = node;
        end

        urdf.materials = materials;
        materialNames = [materials.name];
    end

    %% process all the link elements

    % create a vector of joint structures, fields are:
    % name       name of the joint
    % type       type of joint (fixed, continuous, prismatic)
    % T          SE3 object for pose of joint wrt link frame
    % xyz        position of joint wrt link frame
    % rpy        RPY angles of joint frame wrt link frame
    % node       reference to XML DOM node for this <joint>
    % child      link index for downstream link
    % parent     link index for downstream link

    linkNodes = getChildrenByTagName(robot, 'link');

    for i = 1:length(linkNodes)
        node = linkNodes(i);
        links(i).name = string(node.getAttribute('name'));
        try
            links(i).geometry = node.getElementsByTagName('geometry').item(0);
            links(i).mesh = links(i).geometry.getElementsByTagName('mesh').item(0).getAttribute('filename');
        end
        try
            links(i).material = find(strcmp((node.getElementsByTagName('material').item(0).getAttribute('name')), materialNames));
        end
        try
            inertial = node.getElementsByTagName('inertial').item(0);
            try
                links(i).m = str2num(inertial.getElementsByTagName('mass').item(0).getAttribute('value'));
            end
            try
                I = inertial.getElementsByTagName('inertia').item(0);
                links(i).I.xx = str2num(I.getAttribute('ixx'));
                links(i).I.xy = str2num(I.getAttribute('ixy'));
                links(i).I.xz = str2num(I.getAttribute('ixz'));
                links(i).I.yy = str2num(I.getAttribute('iyy'));
                links(i).I.yz = str2num(I.getAttribute('iyz'));
                links(i).I.zz = str2num(I.getAttribute('izz'));
            end
            try
                t = str2num( node.getElementsByTagName('origin').item(0).getAttribute('xyz'));
                rpy = str2num( node.getElementsByTagName('origin').item(0).getAttribute('rpy'));
                links(i).com = SE3(t) * SE3.rpy(rpy);
            end
        end
        t = [0 0 0]; rpy = [0 0 0]; % in case not provided
        try
            t = str2num( node.getElementsByTagName('origin').item(0).getAttribute('xyz'));
            rpy = str2num( node.getElementsByTagName('origin').item(0).getAttribute('rpy'));
        end
        if isempty(t)
            t = [0 0 0]; % in case the tag has no attribute
        end
        if isempty(rpy)
            rpy = [0 0 0];  % in case the tag has no attribute
        end
        links(i).T = SE3(t) * SE3.rpy(rpy);
        links(i).node = node;
        links(i).children = [];
        links(i).parent = [];
    end

    linkNames = [links.name];  % a list of link names

    %% process all the joint elements

    % create a vector of link structures, fields are are:
    % name       name of the link
    % geometry   reference to XML DOM node
    % T          SE3 object for origin of visual wrt link frame
    % node       reference to XML DOM node for this <link>
    % children   list of joint indexes attaching this link to its children
    % parent     index of joint attaching this link to its parent

    jointNodes = getChildrenByTagName(robot, 'joint');
    joints = [];
    for i = 1:length(jointNodes)
        node = jointNodes(i);
        joints(i).name = string(node.getAttribute('name'));
        joints(i).type = string(node.getAttribute('type'));
        try
            joints(i).axis = str2num(node.getElementsByTagName('axis').item(0).getAttribute('xyz'));
        end
        p = find(strcmp(node.getElementsByTagName('parent').item(0).getAttribute('link'), linkNames));
        c = find(strcmp(node.getElementsByTagName('child').item(0).getAttribute('link'), linkNames));
        joints(i).parent = p;
        joints(i).child = c;
        links(p).children = [links(p).children i];
        links(c).parent = i;

        t = [0 0 0]; rpy = [0 0 0];  % in case not provided
        try
            t = str2num( node.getElementsByTagName('origin').item(0).getAttribute('xyz'));
            rpy = str2num( node.getElementsByTagName('origin').item(0).getAttribute('rpy'));
        end
        if isempty(t)
            t = [0 0 0]; % in case the tag has no attribute
        end
        if isempty(rpy)
            rpy = [0 0 0]; % in case the tag has no attribute
        end
        joints(i).T = SE3(t) * SE3.rpy(rpy);
        joints(i).xyz = t;
        joints(i).rpy = rpy;
        joints(i).node = node;
    end
    urdf.joints = joints;
    urdf.links = links;

    urdf.endlinks = find( cellfun('isempty', {links.children}) );
    urdf.ets = [];

    %% create a transform list

    % the mechanism might have multiple endpoints (if branched), for each
    % endpoint compute an ETS chain and return a string array

    for endlink = urdf.endlinks
        transforms = "";
        joint = joints(links(endlink).parent);
        while true
            transforms = xform2s(joint) + axis2s(joint) + transforms;
            joint = joints(links(joint.parent).parent);
            if isempty(joint)
                break;
            end
        end
        transforms = strtrim(transforms);
        urdf.ets = [urdf.ets transforms];
    end

    %% display a summary
    fprintf('\nLinks::\n');
    for i=1:length(links)
        name = links(i).name;
        if ismember(i, urdf.endlinks)
            name = "*" + name;
        end
        geomlist = "";
        if ~isempty(links(i).geometry)
            g = links(i).geometry.getChildNodes;
            for k=0:g.getLength-1
                c = g.item(k);
                if c.getNodeType == 1  % skip text nodes
                    if geomlist == ""
                        geomlist = "shape=";
                    else
                        geomlist = geomlist + ", ";
                    end
                    geomlist = geomlist + string(c.getNodeName);
                end
            end
        end
        fprintf('%2d: %10s; parentjoint=%s, childjoints=%s; %s\n', i, name, ...
            formatj(links(i).parent), formatj(links(i).children), geomlist);
    end
    if ~isempty(joints)
        fprintf('\nJoints::\n');
        for i=1:length(joints)
            fprintf('%2d: %10s; type=%s, parentlink=%d, childlink=%d\n', i, joints(i).name, joints(i).type, joints(i).parent, joints(i).child);
        end
    end
    if ~isempty(materials)
        fprintf('\nMaterials::\n');
        for i=1:length(materials)
            fprintf('%10s: %f %f %f %f\n', materials(i).name, materials(i).rgba);
        end
    end
    if strlength(transforms) > 0
        fprintf('\nForward kinematics::\n  %s\n', transforms);
    end

end

function s = formatj(j)
    if isempty(j)
        s = '-';
    else
        s = strip(num2str(j, '%d,'), 'right', ',');
    end
end

function t = xform2s(joint)

    t = "";
    if joint.xyz(1) ~= 0
        t = t + sprintf('Tx(%.12g) ', joint.xyz(1));
    end
    if joint.xyz(2) ~= 0
        t = t + sprintf('Ty(%.12g) ', joint.xyz(2));
    end
    if joint.xyz(3) ~= 0
        t = t + sprintf('Tz(%.12g) ', joint.xyz(3));
    end
    if joint.rpy(3) ~= 0
        t = t + sprintf('Rz(%.12g) ', round(joint.rpy(3)*180/pi, 3));
    end
    if joint.rpy(2) ~= 0
        t = t + sprintf('Ry(%.12g) ', round(joint.rpy(2)*180/pi,3));
    end
    if joint.rpy(1) ~= 0
        t = t + sprintf('Rx(%.12g) ', round(joint.rpy(1)*180/pi, 3));
    end
end

function t = axis2s(joint)

    t = "";
    switch joint.type
        case {'revolute', 'continuous'}
            assert(length(find(joint.axis)) == 1, 'Joint axis must be about x, y or z only');
            if joint.axis(1) ~= 0
                t = t + sprintf('Rx(%s) ', joint.name);
            end
            if joint.axis(2) ~= 0
                t = t + sprintf('Ry(%s) ', joint.name);
            end
            if joint.axis(3) ~= 0
                t = t + sprintf('Rz(%s) ', joint.name);
            end
        case 'prismatic'
            assert(length(find(joint.axis)) == 1, 'Joint axis must be about x, y or z only');
            if joint.axis(1) ~= 0
                t = t + sprintf('Tx(%s) ', joint.name);
            end
            if joint.axis(2) ~= 0
                t = t + sprintf('Ty(%s) ', joint.name);
            end
            if joint.axis(3) ~= 0
                t = t + sprintf('Tz(%s) ', joint.name);
            end
        case 'fixed'
    end

end

function children = getChildrenByTagName(parent, name)
    % return array of immediate child nodes with node name equal to name
    % REF https://stackoverflow.com/questions/18776408/get-all-the-children-for-a-given-xml-in-java
    children = [];
    childnodes = parent.getChildNodes;
    for i=1:childnodes.getLength
        node = childnodes.item(i-1);
        if node.getNodeType == node.ELEMENT_NODE && strcmp(node.getNodeName, name)
            children = [children node];
        end
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function model = expvarMeanCreate(inputDim, outputDim, options)

% EXPVARMEANCREATE Creates an the mean function for the EXP kernel.
% FORMAT
% DESC creates a model for returning the first moment of an
% 'exponentiated Gaussian process'. If the output of a Gaussian
% process is exponentiated the resulting process is no longer
% Gaussian. However, it can be approximated by a Gaussian process
% by matching its first and second moments. This function returns
% the mean of that approximating process for a given kernel
% (specified in the options vector). It should be used in tandem
% with the EXP kernel for approximating these Gaussian processes.
% ARG inputDimension : dimension of input to function.
% ARG outputDim : dimension of output from mean function data.
% ARG options : options structure. The structure contains the type
% of kernel that the function is based on. A set of default options
% are given by the file expvarMeanOptions.
% RETURN model : model structure containing the mapping.
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006
%
% SEEALSO : expvarMeanOptions

% SHEFFIELDML

model.type = 'expvarMean';
if isstruct(options.kern)
  model.kern = options.kern;
else
  model.kern = kernCreate(inputDim, options.kern);
end
model.q = inputDim;
model.d = outputDim;
model.numParams = model.kern.nParams;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function dy = FalknerSkanSys(t, y)
    global ETA_INF;
    global BETA0;
    global BETA;
    dy = [ETA_INF*y(2) ETA_INF*y(3) -ETA_INF*(BETA0*y(1)*y(3) + BETA*(1-y(2)^2))]';

repo: mcbench-benchmarks


text:
function MD5 = CalcMD5(Data, InClass, OutClass)  %#ok<STOUT,INUSD>
% 128 bit MD5 checksum: file, string, byte stream [MEX]
% This function calculates a 128 bit checksum for arrays and files.
% Digest = CalcMD5(Data, [InClass], [OutClass])
% INPUT:
%   Data:   Data array or file name. Either numerical or CHAR array.
%           Currently only files and arrays with up to 2^32 bytes (2.1GB) are
%           accepted.
%   InClass: String to declare the type of the 1st input.
%           Optional. Default: 'Char'.
%           'File': [Data] is a file name as string. The digest is calculated
%                   for this file.
%           'Char': [Data] is a char array to calculate the digest for. Only the
%                   ASCII part of the Matlab CHARs is used, such that the digest
%                   is the same as if the array is written to a file as UCHAR,
%                   e.g. with FWRITE.
%           'Unicode': All bytes of the input [Data] are used to calculate the
%                   digest. This is the standard for numerical input.
%   OutClass: String, format of the output. Just the first character matters.
%           Optional, default: 'hex'.
%           'hex': [1 x 32] string as lowercase hexadecimal number.
%           'HEX': [1 x 32] string as uppercase hexadecimal number.
%           'Dec': [1 x 16] double vector with UINT8 values.
%           'Base64': [1 x 22] string, encoded to base 64 (A:Z,a:z,0:9,+,/).
%
% OUTPUT:
%   Digest: A 128 bit number is replied in a format depending on [OutClass].
%           The chance, that different data sets have the same MD5 sum is about
%           2^128 (> 3.4 * 10^38). Therefore MD5 can be used as "finger-print"
%           of a file rather than e.g. CRC32.
%
% EXAMPLES:
%   Three methods to get the MD5 of a file:
%   1. Direct file access (recommended):
%     MD5 = CalcMD5(which('CalcMD5.m'), 'File')
%   2. Import the file to a CHAR array (binary mode for exact line breaks!):
%     FID = fopen(which('CalcMD5.m'), 'rb');
%     S   = fread(FID, inf, 'uchar=>char');
%     fclose(FID);
%     MD5 = CalcMD5(S, 'char')
%   3. Import file as a byte stream:
%     FID = fopen(which('CalcMD5.m'), 'rb');
%     S   = fread(FID, inf, 'uint8=>uint8');
%     fclose(FID);
%     MD5 = CalcMD5(S, 'unicode');  // 'unicode' can be omitted here
%
%   Test string:
%     CalcMD5(char(0:511), 'char', 'HEX')
%       => F5C8E3C31C044BAE0E65569560B54332
%     CalcMD5(char(0:511), 'unicode', 'HEX')
%       => 3484769D4F7EBB88BBE942BB924834CD
%
% Tested: Matlab 6.5, 7.7, 7.8, WinXP, [UnitTest]
% Author: Jan Simon, Heidelberg, (C) 2009-2010 J@n-Simon.De
% License: This program is derived from the RSA Data Security, Inc.
%          MD5 Message Digest Algorithm, RFC 1321, R. Rivest, April 1992
%
% See also CalcCRC32.
% Michael Kleder has published a Java call to compute the MD5 (and further
% check sums): http://www.mathworks.com/matlabcentral/fileexchange/8944

% $JRev: R5.00j V:015 Sum:zh2gTrvHwbd7 Date:17-Dec-2009 02:46:53 $
% $File: CalcMD5\CalcMD5.m $
% History:
% 015: 15-Dec-2009 16:53, BUGFIX: UINT32 has 32 bits on 64 bit systems now.
%      Thanks to Sebastiaan Breedveld!

% If the current Matlab path is the parent folder of this script, the
% MEX function is not found - change the current directory!
error(['JSim:', mfilename, ':NoMex'], 'Cannot find MEX script.');


text:
function [ xa, xb, fxa, fxb ] = secant ( fatol, step_max, prob, xatol, xmin, ...
  xmax, xa, xb, fxa, fxb )

%*****************************************************************************80
%
%% SECANT carries out the secant method to seek a root of F(X) = 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real FATOL, an absolute error tolerance for the
%    function value of the root.  If an approximate root X satisfies
%      ABS ( F ( X ) ) <= FATOL, then X will be accepted as the
%    root and the iteration will be terminated.
%
%    Input, integer STEP_MAX, the maximum number of steps allowed
%    for an iteration.
%
%    Input, integer PROB, the index of the function whose root is
%    to be sought.
%
%    Input, real XATOL, an absolute error tolerance for the root.
%
%    Input, real XMAX, XMIN, the interval in which the root should
%    be sought.
%
%    Input/output, real XA, XB, two points at which the
%    function differs in sign.  On output, these values have been adjusted
%    to a smaller interval.
%
%    Input/output, real FXA, FXB, the value of the function
%    at XA and XB.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SECANT\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Step             X             F(X)\n' );
  fprintf ( 1, '\n' );

  step_num = -1;
  fprintf ( 1, '  %4d  %14g  %14g\n', step_num, xa, fxa );

  if ( abs ( fxa ) <= fatol )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Function small enough for convergence.\n' );
    return
  end

  step_num = 0;
  fprintf ( 1, '  %4d  %14g  %14g\n', step_num, xb, fxb );

  for step_num = 1 : step_max

    if ( abs ( fxb ) <= fatol )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Function small enough for convergence.\n' );
      return
    end

    if ( abs ( xa - xb ) < xatol )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Interval small enough for convergence.\n' );
      return
    end

    if ( xb < xmin || xmax < xb )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Iterate has left the region [XMIN,XMAX].\n' );
      return
    end

    if ( fxa == fxb )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  F(A) = F(B), algorithm fails.\n' );
      return
    end

    xc = ( fxa * xb - fxb * xa ) / ( fxa - fxb );

    fxc = p00_fx ( prob, xc );

    xa = xb;
    fxa = fxb;
    xb = xc;
    fxb = fxc;
    fprintf ( 1, '  %4d  %14g  %14g\n', step_num, xb, fxb );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Took maximum number of steps.\n' );

  return
end


text:
function varargout = kullbackleibler(varargin)
% KULLBACKLEIBLER
%
% y = KULLBACKLEIBLER(x,y)
%
% Computes/declares the convex Kullback-Leibler divergence sum(x.*log(x./y))
% Alternatively -sum(x.*log(y/x)), i.e., negated perspectives of log(y)
%
% See also ENTROPY, CROSSENTROPY

switch class(varargin{1})

    case 'double'
        if nargin == 1
            z = varargin{1};
            z = reshape(z,[],2);
            x = z(:,1);
            y = z(:,2);
        else
            x = varargin{1}(:);
            y = varargin{2}(:);
        end
        l = log(x./y);
        l(x<=0) = 0;
        l = real(l);
        varargout{1} = sum(x.*l);

    case {'sdpvar','ndsdpvar'}

        varargin{1} = reshape(varargin{1},[],1);
        varargin{2} = reshape(varargin{2},[],1);

        if length(varargin{1})~=length(varargin{2})
            if length(varargin{1})==1
                varargin{1} = repmat(varargin{1},length(varargin{2}),1);
            elseif  length(varargin{2})==1
                varargin{2} = repmat(varargin{2},length(varargin{1}),1);
            else
                error('Dimension mismatch in kullbackleibler')
            end
        end

        varargout{1} = yalmip('define',mfilename,[varargin{1};varargin{2}]);

    case 'char'

        X = varargin{3};
        F = [X >= 0];

        operator = struct('convexity','convex','monotonicity','none','definiteness','none','model','callback');
        operator.range = [-inf inf];
        operator.domain = [0 inf];
        operator.derivative = @derivative;

        varargout{1} = F;
        varargout{2} = operator;
        varargout{3} = X;

    otherwise
        error('SDPVAR/KULLBACKLEIBLER called with CHAR argument?');
end

function df = derivative(x)
z = abs(reshape(x,[],2));
x = z(:,1);
y = z(:,2);
% Use KL = -Entropy + Cross Entropy
df = [1+log(x);0*y] + [-log(y);-x./y];

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function i4_is_odd_test ( )

%*****************************************************************************80
%
%% I4_IS_ODD_TEST tests I4_IS_ODD.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 December 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_IS_ODD_TEST\n' );
  fprintf ( 1, '  I4_IS_ODD reports whether an I4 is odd.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I     I4_IS_ODD(I)\n' );
  fprintf ( 1, '\n' );

  for i = -2 : 25
    fprintf ( 1, '  %6d  %1d\n', i, i4_is_odd ( i ) );
  end

  return
end


text:
% weighting
function [gen, weight_a, weight_b] = fusion_strategy(features_a, features_b, source_a, source_b)

[m1,n1] = size(source_a);
% resize
resize_temp1 = imresize(features_a, [m1, n1]);
resize_temp2 = imresize(features_b, [m1, n1]);
% soft-max
weight_ave_temp1 = resize_temp1./(resize_temp1+resize_temp2);
weight_ave_temp2 = resize_temp2./(resize_temp1+resize_temp2);
% reconstruction
gen = source_a.*weight_ave_temp1 + source_b.*weight_ave_temp2;
% figure;imshow(gen);

weight_a = weight_ave_temp1;
weight_b = weight_ave_temp2;
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
hydro = struct();

hydro = readWAMIT(hydro,'coer_comp.out',[]);
hydro = radiationIRF(hydro,10,[],[],[],[]);
hydro = radiationIRFSS(hydro,[],[]);
hydro = excitationIRF(hydro,10,[],[],[],[]);
writeBEMIOH5(hydro)
plotBEMIO(hydro)

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% Template for setup data
%     (here: use brain MRI's; see also setup3DbrainsData.m for references)
%
% This file initializes the following data (which can then be saved):
%   dataT       template  image, a d-array of size mD,
%   dataR       reference image, a d-array of size nD,
%   omega       domain specification
%             omega = (omega(1),omega(2)) x  (omega(3),omega(24))
%   m           initial discretization
%   ML          multi-level representation of the data
%   LM        landmarks, if available
%
% For representation and visualization
%   viewPara      options for image viewer
%   imgPara     options for image interpolation
% see also setup3DbrainsData and E2_setupHandsData
%==============================================================================

% load 3D data
load brain3D;
whos
viewPara  = {'viewImage','imgmontage'};
viewImage('reset',viewPara{:});

imgPara = {'imgModel','linearInter'};
imgModel('reset',imgPara{:});
FAIRfigure(2); colormap(gray(256));
ML    = getMultilevel({dataT,dataR},omega,m);
save(fullfile(FAIRpath,'temp','brain3DML.mat'),...
	'dataT','dataR','omega','m','ML','viewPara','imgPara');
%==============================================================================


text:
function if_smooth=filt_if_law(if_law,window_length)
%---------------------------------------------------------------------
% filter IF law by convolving with Hamming window
%---------------------------------------------------------------------
if(nargin<2 || isempty(window_length)) window_length=15; end

w=hamming( floor(window_length) ); w=w./sum(w);

if_mean=mean(if_law);
if_smooth=conv(if_law-if_mean,w,'same');
if_smooth=if_smooth+if_mean;

author: HeLiangHIT


text:
% 计算欧式距离
function dist = distEclud(vecA,vecB)
    dist  = sum(power((vecA-vecB),2));
end

function dist = softDist(vecA,vecB)
    dist = sum(abs(vecB-vecA));
end

repo: MachineLearning


text:
function varargout = bst_meanvar(varargin)
%BST_MEANVAR: Mex-file to compute mean and variance along first dimension (with imprecise algorithm)
%
% USAGE: [mean,var,nAvg] = bst_meanvar(x, isZeroBad)
%
% INPUTS:
%    - x         : [NxM] double matrix with values to process
%    - isZeroBad : If 1, excludes all the zero values from the computation
%
% OUTPUTS:
%    - mean : [1xM] averages values
%    - var  : [1xM] unbiased estimator of the variance (computed with an algorithm prone to rounding errors)
%    - nAvg : [1xM] number of non-zero values that were averaged
%
% COMPILE:
%    mex -v bst_meanvar.c

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

error('Mex-function bst_meanvar.c not compiled.');

author: brainstorm-tools


text:
function [y,u,f_fname,g_fname,dim,options,in]    = VBA_getDefaults()
% outputs default input to VBA_NLStateSpaceModel.m
% function [y,u,f_fname,g_fname,dim,options,in]    = VBA_getDefaults()
%
% This function provides the defaults inputs to the main
% VBA_NLStateSpaceModel.m routine. See VBA_check.m routine.

%------- Main inputs ------------%
y                   = 'pxn_t measured data matrix';
u                   = 'n_uxn_t driving input matrix (can be left empty)';
f_fname             = 'name/handle of the evolution function';
g_fname             = 'name/handle of the observation function';
dim.n               = 'dimension of the state space';
dim.n_theta         = 'dimension of the evolution parameter vector';
dim.n_phi           = 'dimension of the observation parameter vector';
in                  = 'only for inversion refinement (can be unspecified)';

%------- Options structure -------%
options.priors = 'structure containing the first two moments of all prior pdfs [see VBA_defaultPriors()]';
options.decim = 'for micro-time resolution';
options.microU = 'flag for micro-time input u';
options.inF = '(internal) parameters of the evolution function';
options.inG = '(internal) parameters of the observation function';
options.checkGrads = 'Check analytical gradients against numerical gradients';
options.updateX0 = 'Flag for initial conditions update';
options.updateHP = 'Flag for precision hyperparameters update';
options.MaxIter = 'Maximum number of iterations';
options.MinIter = 'Minimum number of iterations';
options.TolFun = 'Minimum change in the free energy ';
options.DisplayWin = 'VB display window';
options.GnMaxIter = 'Gauss-Newton inner-loops maximum number of iterations';
options.GnTolFun = 'Gauss-Newton inner-loops threshold on relative variational energy';
options.GnFigs = 'Gauss-Newton inner-loops figures';
options.verbose = 'summary messages in MATLAB main window';
options.finalEval = 'Externally-specified function evaluation: end of VB algo';
options.figName = 'Name of the display figure';
options.delays = 'dim.nX1 vector for delay embedding';
options.sources = 'Type and dimension of the data distribution(s)';

author: MBB-team

repo: VBA-toolbox


text:
function detections = nmsProcess(det, isNMS, nmsThre)
if(isNMS)
    detections = [];
    numFr = max(det(:,1));
    for i = 1:numFr
        idx = det(:,1) == i;
        curdet = det(idx,:);
        pick = nms(curdet(:,3:7), nmsThre);
        detections = cat(1, detections, curdet(pick,:));
    end
else
    detections = det;
end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
% Cdiff_test.m
% test Cdiff object

if 1 % look at 2nd order case
	nx = 8; ny = 6;
	C = Cdiff(ones(nx,ny), 'order', 2);
	Cf = C(:,:);
	t = reshape(Cf', nx*ny, nx*ny, []);
	im(t)
	t = zeros(nx,ny); t(nx/2+1,ny/2+1) = 1;
	t = C' * (C * t);
	im(t)
end

if 1
	nx = 16; ny = 14;
	nx = 512; ny = 500; % for large images, the mex file is much faster!
	ig = image_geom('nx', nx, 'ny', ny, 'dx', 1);
	ig.mask = [0 0 [nx ny]/2-5 0 1];
	ig.mask = conv2(double(ellipse_im(ig, []) > 0), ones(2), 'same') > 0;
%	mask = ones(nx,ny); % all

	if 1, printm 'test penalty_mex'
		x = single(ig.mask);
		offsets = [nx-1];
		offsets = int32(offsets);
		d1 = penalty_mex('diff2,forw1', x, offsets);
		d2 = penalty_mex('diff2,forw1', x, offsets, int32(ndims(x)));
		if any(d1(:) ~= d2(:)), error 'bug', end

		x1 = penalty_mex('diff2,back1', d1, offsets);
		x2 = penalty_mex('diff2,back1', d1, offsets, int32(ndims(x)));
		if any(x1(:) ~= x2(:)), error 'bug', end
	end

	ctype = 'leak';
	ctype = 'tight';
	order = 1;
	tic
	C1 = Cdiff(ig.mask, 'edge_type', ctype, 'offsets', '2d,hvd', ...
		'distance_power', 1., 'order', 1);
	printm('make C1 time %g', toc)

	tic
	[C2 wjk] = C2sparse(ctype, ig.mask, 8);
	printm('make C2 time %g', toc)
	C2 = spdiag(sqrt(wjk)) * C2;
	C2 = C2(:,ig.mask(:));

	if 0 % test old Cmask
		cpu etic
		b1 = Cmask('tight,2d,hvd', ig.mask);
		cpu etoc 'make scale time:'
		b1(:,:,[3 4]) = b1(:,:,[3 4]) / sqrt(sqrt(2));
		if 0
			b1 = reshape(sqrt(wjk), [nx ny 4]);
		end

		cpu etic
		b2 = penalty_mex('scales,tight', single(ig.mask), C1.arg.offsets, 1.);
		b2 = double(b2);
		cpu doc 'make scale time:'

		im plc 1 3, im(1, b1), im(2, b2), im(3, b1-b2)
		printm('old vs new: %g%%', max_percent_diff(b1, b2))
		equivs(b1, b2)
	return
	end

	rng(0)
	x = rand(nx, ny);
	x = dsingle(x);
	x = x .* ig.mask;

%prompt
end

if 1
	xm = double(x(ig.mask(:)));
	cpu etic
	d1 = C1 * x;
	cpu etoc 'C1 forw time:'

	cpu etic
	d2 = C2 * xm;
	cpu etoc 'C2 forw time:'
	d2 = reshape(d2, [nx ny 4]);
%	printm('Cx vs penalty_mex: %g%%', max_percent_diff(d1, d2))
	equivs(d1, d2)

	if im
		im plc 1 3
		im(1, d1, 'C * x'), cbar
		im(2, d2, 'penalty--mex'), cbar
		im(3, d2-d1, 'err'), cbar
	prompt
	end
end

if 1 % [x x]
	xx = double([xm xm]);
	cpu etic
	d11 = C1 * xx;
	cpu etoc 'C1 forw time:'

	cpu etic
	d22 = C2 * xx;
	cpu etoc 'C2 forw time:'
%	printm('Cxx vs penalty_mex: %g%%', max_percent_diff(d11, d22))
	equivs(d11, d22)
end

if 1
	d = double(d1(:));
	%d = zeros(nx,ny,4);
	%d(end/2,end/2,1) = 1;

	cpu etic
	x1 = C1' * d;
	cpu etoc 'C1''d time:'
	x1 = embed(x1, ig.mask);

	cpu etic
	x2 = C2' * d;
	cpu etoc 'C2''d time:'
	x2 = embed(x2, ig.mask);

%	printm('C1 vs C2: %g%%', max_percent_diff(x1, x2))
	equivs(x1, x2)

	if im
		im plc 1 3
		im(1, x1, 'C''d'), cbar
		im(2, x2, 'penalty--mex'), cbar
		im(3, x2-x1, 'err'), cbar
	end
end

if 1 % C' [d d]
	dd = double([d(:) d(:)]);

	cpu etic
	x11 = C1' * dd;
	cpu etoc 'C1''d time:'

	cpu etic
	x22 = C2' * dd;
	cpu etoc 'C2''d time:'

%	printm('C1 vs C2 for dd: %g%%', max_percent_diff(x11, x22))
	equivs(x11, x22)
end

author: JeffFessler

repo: mirt


text:
function calpak_test26 ( )

%*****************************************************************************80
%
%% CALPAK_TEST26 tests JED_TO_YMDF_SAKA and YMDF_TO_JED_SAKA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST26\n' );
  fprintf ( 1, '  For the Saka calendar:\n' );
  fprintf ( 1, '  JED_TO_YMDF_SAKA: JED -> YMDF.\n' );
  fprintf ( 1, '  YMDF_TO_JED_SAKA: YMDF -> JED.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  JED (in)    YMDF                JED (out)\n' );
  fprintf ( 1, '\n' );

  jed_epoch = epoch_to_jed_saka ( );

  i = 0;

  while ( 1 )

    i = i + 1;
    jed1 = jed_test ( i );

    if ( jed1 < 0.0 )
      break
    end

    if ( jed_epoch <= jed1 )

      [ y2, m2, d2, f2 ] = jed_to_ymdf_saka ( jed1 );

      s2 = ymdf_to_s_numeric ( y2, m2, d2, f2 );

      jed3 = ymdf_to_jed_saka ( y2, m2, d2, f2 );

      fprintf ( 1, '  %11.2f  %20s  %11.2f\n', jed1, s2, jed3 );

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function asa266_test06 ( )

%*****************************************************************************80
%
%% TEST06 tests GAMAIN, GAMMDS, GAMMAD.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 August 2010
%
%  Author:
%
%    John Burkardt
%
  ntest = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST06\n' );
  fprintf ( 1, '  GAMAIN,\n' );
  fprintf ( 1, '  GAMMDS and \n' );
  fprintf ( 1, '  GAMMAD compute the incomplete Gamma integral.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  X  P  GAMMDS  GAMMAD  GAMAIN\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ntest

    x = i / ntest;

    fprintf ( 1, '\n' );

    for j = 1 : ntest

      p = j / ntest;
      [ g1, ifault ] = gammds ( x, p );
      if ( ifault ~= 0 )
        g1 = -99.0;
      end

      [ g2, ifault ] = gammad ( x, p );
      if ( ifault ~= 0 )
        g2 = -99.0;
      end

      [ g3, ifault ] = gamain ( x, p );
      if ( ifault ~= 0 )
        g3 = - 99.0;
      end

      fprintf ( 1, '  %12f  %12f  %12f  %12f  %12f\n', x, p, g1, g2, g3 );

    end

  end

  return
end


text:
%% This file is the main file of using the SINDy-PI method to
% infer the Single Pendulum on Cart Model.
%
% Date: 2019/07/26
% Coded By: K

%% Close all, clear all, clc
close all;clear all; clc;
dir
set(0,'defaulttextInterpreter','latex')
[status,message,messageid] = mkdir('Results');
addpath('Function');
%% Simulate the budworm population growth and gather the simulation data

%Define the model parameters
M=1;m=1;L=1;g=9.81;

% Define noise level and add gaussian noise to the data
noise=0.02;

% Define whehter you have control, if you have it, please define it
Control=1;

%Define whether you want to shuffel the final data
Shuffle=0;

% Peform simulation
state0=[0.3;0;1;0];state0_test=[0.1;0;0.3;0];
dt=0.001;T=16;T_test=2;
tspan=0:dt:T;tspan_test=0:dt:T_test;
u=-0.2+0.5*sin(6*(tspan'));
u_test=-1+1*sin((tspan_test'))+3*sin(2*(tspan_test'));
[dData,Data]=Get_Sim_Data(@(t,y,u)SinglePendulum_ODE(t,y,u,M,m,L,g),state0,u,tspan,noise,Control,Shuffle);
dData(:,1)=Data(:,3);dData(:,2)=Data(:,4);
[dData_test,Data_test]=Get_Sim_Data(@(t,y,u)SinglePendulum_ODE(t,y,u,M,m,L,g),state0_test,u_test,tspan_test,noise,Control,Shuffle);
dData_test(:,1)=Data_test(:,3);dData_test(:,2)=Data_test(:,4);
%% Plot the simulation data
% close all
% HorizontalRange=5*L;
% VerticalRange=2.5*L;
% displayrate=1;
% cartoonOffLineAcce(Data(:,2),Data(:,1),L,HorizontalRange,VerticalRange,displayrate)
%%
% figure(1)
% plot(tspan,Data(:,1))
% title('Theta')
% grid on
% %
% figure(2)
% plot(tspan,Data(:,2))
% title('Position')
% grid on
% %
% figure(3)
% plot(tspan,Data(:,3))
% title('dTheta')
% grid on
% %
% figure(4)
% plot(tspan,Data(:,4))
% title('Velocity')
% grid on
% %
% figure(5)
% plot(tspan,u,'linewidth',3,'color','black')
% box('off')
% axis('on')

%% Now perform sparse regression of non-linear dynamics

% Get the number of states we have
[dtat_length,n_state]=size(Data);

% Define the control input(Should be zero in our example)
n_control=1;

% Choose whether you want to display actual ODE or not
disp_actual_ode=1;

% If the ODEs you want to display is the actual underlyting dynamics of the
% system, please set actual as 1
actual=1;

% Print the actual ODE we try to discover
digits(4)
z_vars=sym('z',[1,n_state]);
u_vars=sym('u',[1,n_control]);
d_vars=sym('dz',[1,n_state]);
ODEsP=SinglePendulum_ODE(0,z_vars,u_vars,M,m,L,g);
fprintf('The actual ODE of the system is/are :\n')
for i=1:n_state
    fprintf(strcat(char(d_vars(1,i)),'=',char(ODEsP(i,1)),'\n'));
end

% The implicit ODE has the following form:
% dz1=z3
% dz2=z4
% dz3=((981*sin(z1))/50 + u1*cos(z1) + z3^2*cos(z1)*sin(z1))/(cos(z1)^2 - 2)
% dz4=-(u1 + z3^2*sin(z1) + (981*cos(z1)*sin(z1))/100)/(cos(z1)^2 - 2)

% Create symbolic states
dz=sym('dz',[n_state,1]);
z=sym('z',[n_state,1]);
uc=sym('u',[n_control,1]);

% Now we first create the parameters of the function right hand side
Highest_Poly_Order_Guess=1;
Highest_Trig_Order_Guess=2;
Highest_U_Order_Guess=0;

% Then create the right hand side library parameters
Highest_Poly_Order=2;
Highest_Trig_Order=2;
Highest_U_Order=0;
Highest_dPoly_Order=1;

%% Define parameters for the sparese regression
lam=[1e-4;5e-4;1e-3;2e-3;3e-3;4e-3;5e-3;6e-3;7e-3;8e-3;9e-3;1e-2;2e-2;3e-2;4e-2;5e-2;...
    6e-2;7e-2;8e-2;9e-2;1e-1;2e-1;3e-1;4e-1;5e-1;6e-1;7e-1;8e-1;9e-1;1;1.5;2;2.5;3;3.5;4;4.5;5;...
    6;7;8;9;10;20;30;40;50;100;200];

N_iter=20;
disp=0;
LibPoly=[1;1;0;0];
LibTrig=[0;0;3;4];
%
GuessPoly=[1;1;0;0];
GuessTrig=[0;0;1;1];

% Normalize the library?
NormalizeLib=0;

% Define a cell matrix to store the variable
Xi_Final=cell(n_state,1);
Score_Final=cell(n_state,1);

for iter=1:n_state
    % Change the library base on what equation you want to work on
    Highest_Poly_Order=LibPoly(iter);
    Highest_Poly_Order_Guess=GuessPoly(iter);
    Highest_Trig_Order=LibTrig(iter);
    Highest_Trig_Order_Guess=GuessTrig(iter);

    fprintf('Calculating the %i expression \n',iter)

    % According to the previous parameter generate the left hand side guess
    [LHS_Data,LHS_Sym]=GuessLib(Data,dData(:,iter),iter,u,Highest_Poly_Order_Guess,Highest_Trig_Order_Guess,Highest_U_Order_Guess);

    %Generate the corresponding data
    [SINDy_Data,SINDy_Struct]=SINDyLib(Data,dData(:,iter),iter,u,Highest_Poly_Order,Highest_Trig_Order,Highest_U_Order,Highest_dPoly_Order);

    % Run the for loop and try all the left hand guess
    for i=1:length(LHS_Sym)
        fprintf('\n\t\t Testing the left hand side as %s... \n',cell2sym(LHS_Sym(1,i)))

        % Exclude the guess from SINDy-PI library
        [RHS_Data,RHS_Struct]=ExcludeGuess(SINDy_Data,SINDy_Struct,LHS_Sym{i});

        if iter==1
            Xi=cell(length(LHS_Sym),length(lam));
            Score=zeros(length(LHS_Sym),length(lam));
            ODE=cell(length(LHS_Sym),length(lam));
        end

        % Define dummy variables for parfor
        LHS_Sym_Dum=LHS_Sym{i};
        LHS_Data_Dum=LHS_Data(:,i);
        dz_Dum=dz(iter);
        dData_test_dum=dData_test(:,iter);
        for j=1:length(lam)
            fprintf('\n\t\t Testing the lambda as %d... \n',lam(j,1))
            % Perform the sparse regression problem
            [Xi{i,j},ODE{i,j}]=sparsifyDynamics(RHS_Data,LHS_Data_Dum,LHS_Sym_Dum,lam(j,1),N_iter,RHS_Struct,disp,NormalizeLib);

            % Perform sybolic calculation and solve for dX
            Eqn=LHS_Sym_Dum==ODE{i,j};
            digits(4)
            ODE_Guess=simplify(vpa(solve(Eqn,dz_Dum)));
            % Store the result of each valuation
            try
                ODEs(i,j)=ODE_Guess;
                func=matlabFunction(ODEs(i,j),'Vars',{z,uc});
                funcVal=func(Data_test',u_test')';
                Score(i,j)=norm(dData_test_dum-funcVal)/norm(dData_test(:,iter));
            catch
                ODEs(i,j)=0;
                Score(i,j)=NaN;
            end
        end
    end

    % Calculate the minimum score
    [minVal,minIndex]=min(Score,[],2);
    [minVal2,minIndex2]=min(minVal);

    % Store the best ODE approximation
    ODE_Best(iter,1)=ODEs(minIndex2,minIndex(minIndex2));

    % Print the result
    fprintf('\n \tThe SINDy-PI discovered Best ODE for %i equation is:\n',iter)
    digits(4)
    fprintf(strcat(char(dz(iter)),'=',char(simplify(ODE_Best(iter,1))),'\n'));

    % Save this result
    Xi_Final{iter,1}=Xi;
    Score_Final{iter,1}=Score;

    % Plot the score
    figure(iter)
    hold on
    for k=1:length(LHS_Sym)
        plot(lam,Score(k,:),'color',[1 iter*0.2 iter*0.25],'linewidth',2.5)
        scatter(lam,Score(k,:),100,'MarkerFaceColor',[1 iter*0.2 iter*0.25],'MarkerEdgeColor',[0 0 0],'linewidth',2.5)
        set(gca,'XScale','log')
    end

end

%% Now generate the ODE function file and test the accuracy of the
% identified system
Noise=0;

% Now generate this best guess ODE
Generate_ODE_RHS(ODE_Best(:,1),n_state,n_control);

% Simulate the system with the SINDy-PI identified ODE
%state0_test=[0.5;0.1;0.1;-0.1];
state0_test=[pi;0;0;0];
u_test=-0.5+0.2*sin((tspan_test'))+0.3*sin(2*(tspan_test'));
LHS_Sym=0;
[dData_test,Data_test]=Get_Sim_Data(@(t,y,u)SinglePendulum_ODE(t,y,u,M,m,L,g),state0_test,u_test,tspan_test,Noise,Control,Shuffle);
[d_Data_Es,Data_Es]=Get_Sim_Data(@(t,y,u)Sindy_ODE_RHS(t,y,u),state0_test,u_test,tspan_test,Noise,Control,Shuffle);

%% Print the Result
disp_best=1;
if disp_best==1
    fprintf('The SINDy-PI discovered Best ODE is:\n')
    digits(4)
    for i=1:n_state
        fprintf(strcat(char(dz(i)),'=',char(simplify(ODE_Best(i,1))),'\n'));
    end
end
%% Save the result
File_Name=strcat('Results/Noise_Level_',num2str(noise),'.mat');
save(File_Name,'Xi_Final','Score_Final','ODE_Best','state0_test',...
    'u_test','dData_test','Data_test','d_Data_Es','Data_Es','tspan_test')

%% Plot the simulation data
close all
% Create the new directory to save the plot
[fld_status, fld_msg, fld_msgID]=mkdir('Figures');
%
figure(1)
plot(tspan_test,Data_test(:,1),'linewidth',4.5,'Color','black')
hold on
plot(tspan_test,Data_Es(:,1),'linewidth',4.5,'linestyle','--','color','blue')
% legend('Actual Dynamics','Best Model')
% title('Validation','FontSize',24)
% xlabel('Time $(t)$','FontSize', 24)
% ylabel('$\theta(t)$','FontSize', 24)
set(gca,'FontSize',34);
grid on
h = gca;
set(gca,'xticklabel',[])
set(gcf,'Position',[100 100 600 400]);
set(gcf,'PaperPositionMode','auto');
box('on')
%print('-depsc2', '-loose', 'Figures/SinglePendulum_Thetat1.eps');

%
figure(2)
plot(tspan_test,Data_test(:,2),'linewidth',4.5,'Color','black')
hold on
plot(tspan_test,Data_Es(:,2),'linewidth',4.5,'linestyle','--','color','blue')
% legend('Actual Dynamics','Best Model')
% title('Validation','FontSize',18)
% xlabel('Time $(t)$','FontSize', 18)
% ylabel('$x(t)$','FontSize', 18)
set(gca,'FontSize',34);
grid on
%h = gca; h.XAxis.Visible = 'off';
set(gcf,'Position',[100 100 600 400]);
set(gcf,'PaperPositionMode','auto');
box('on')
print('-depsc2', '-loose', 'Figures/SinglePendulum_x.eps');

%%
figure(3)
plot(tspan_test,Data_test(:,3),'linewidth',2,'Color','green')
hold on
plot(tspan_test,Data_Es(:,3),'linewidth',2,'linestyle','--','color','blue')
legend('Actual Dynamics','Best Model')
title('Validation','FontSize',18)
xlabel('Time $(t)$','FontSize', 18)
ylabel('$\dot{\theta}(t)$','FontSize', 18)
set(gca,'FontSize',18);
grid on

set(gcf,'Position',[100 100 600 400]);
set(gcf,'PaperPositionMode','auto');
print('-depsc2', '-loose', 'Figures/SinglePendulum_dThetat.eps');

%%
figure(4)
plot(tspan_test,Data_test(:,4),'linewidth',2,'Color','green')
hold on
plot(tspan_test,Data_Es(:,4),'linewidth',2,'linestyle','--','color','blue')
legend('Actual Dynamics','Best Model')
title('Validation','FontSize',18)
xlabel('Time $(t)$','FontSize', 18)
ylabel('$\dot{x}(t)$','FontSize', 18)
set(gca,'FontSize',18);
grid on

set(gcf,'Position',[100 100 600 400]);
set(gcf,'PaperPositionMode','auto');
print('-depsc2', '-loose', 'Figures/SinglePendulum_dx.eps');

%
figure(5)
plot(Data_test(:,1),Data_test(:,3),'linewidth',2,'Color','green')
hold on
plot(Data_Es(:,1),Data_Es(:,3),'linewidth',2,'linestyle','--','color','blue')
legend('Actual Dynamics','Best Model')
title('Validation','FontSize',18)
xlabel('$\theta(t)$','FontSize', 18)
ylabel('$\dot{\theta}(t)$','FontSize', 18)
set(gca,'FontSize',18);
grid on

set(gcf,'Position',[100 100 600 400]);
set(gcf,'PaperPositionMode','auto');
print('-depsc2', '-loose', 'Figures/SinglePendulum_theta_vs_dtheta.eps');

%
figure(6)
plot(Data_test(:,2),Data_test(:,4),'linewidth',2,'Color','green')
hold on
plot(Data_Es(:,2),Data_Es(:,4),'linewidth',2,'linestyle','--','color','blue')
legend('Actual Dynamics','Best Model')
title('Validation','FontSize',18)
xlabel('$x(t)$','FontSize', 18)
ylabel('$\dot{x}(t)$','FontSize', 18)
set(gca,'FontSize',18);
grid on

set(gcf,'Position',[100 100 600 400]);
set(gcf,'PaperPositionMode','auto');
print('-depsc2', '-loose', 'Figures/SinglePendulum_x_vs_dx.eps');

author: dynamicslab

repo: SINDy-PI


text:
function [surfVel] = getSurfaceVelocity(bodyInfo, ut, rVectECI, vVectECI)
%getSurfaceVelocity Summary of this function goes here
%   Detailed explanation goes here
    [~, ~, ~, ~, horzVel, ~] = getLatLongAltFromInertialVect(ut, rVectECI, bodyInfo, vVectECI);
%     bodyCircum = 2*pi*bodyInfo.radius;
%     equatSpeed = bodyCircum / bodyInfo.rotperiod; %km/s
%
%     speedAtLat = equatSpeed * cos(lat);
%     bodySurfVelECI = -normVector(cross(rVectECI, [0;0;1])) * speedAtLat;
    surfVel = horzVel;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
clear
%Process definition
ar = rc2arset([1 -.4 .2]);
ma = rc2arset([1 .8]);
%End of Process definition
n_obs = 10;

rch = [1 -.6 -.3 -.2];
Lar = length(rch)-1;
arh_set = rc2arset(rch,0:Lar);

%function pe_set = prederrAR(rc,cov)
%PREDERRAR Prediction error of AR models of increasing order
vareps = 1/pgain(ar,ma);
cor = arma2cor(ar,ma,Lar);

r = [cor(end:-1:1) cor(2:end)];
pe(1) = 1;
eo = r; %cor;
fo = r; %cor;
pe_old(1) = 1;
for p = 1:Lar
    carh = [arh_set{p+1} zeros(1,Lar-p)];
    crc = rch(p+1);
    ahahp = convol(fliplr(carh),carh)
    ahah  = convol(carh,carh);
    e = xcorr(ahah,r);  e = e(Lar+1:3*Lar+1);
    f = xcorr(ahahp,r); f = f(2*Lar+1:3*Lar+1);
    %Recursieve berekening van e en f.
    er = NaN*ones(1,Lar+1);
    for q = p:Lar
        er(q+1) = eo(q+1) + 2*crc*fo(q-p+1) + crc^2*eo(2*p-q+1);
    end
    fr = NaN*ones(1,Lar+1);
%     for q = 0:Lar-p
%         fr(q+1) = (1+crc^2)*fo(q+1) + crc*eo(p+q+1); %Klopt alleen voor witte ruis!
%     end
%     e,er
%     f,fr;
    disp(' ')
    pe(p+1) = f(Lar+1) %werkt

    %Preparations for next step
    eo = e;
    fo = f;

    %Old PE
    pe_old(p+1) = (moderr(carh,1,ar,ma,1)+1)*vareps;
end

author: Sable


text:
% Conduct a binary search
%
% binSearch: Conducts a binary search and continues until the score is
% close enough to the target. An upper limit should be set on the number of
% iterations
%
%     [SearchPoint, Step, ReLoop] = iosr.auditory.binSearch(Score, ...
%       TargetScore, ...
%       CloseEnough, ...
%       NextSearchPointa, ...
%       Step, ...
%       LoopCount, ...
%       MaxLoops);
% end
%
% inputs:
% - Score: the test value
% - Target: the 'finished' test value
% - CloseEnough: the distance from the Target at which it is acceptable to
% discontinue the binary search
% - SearchPoint: the next work input value to try
% - Step: the distance which the SearchPoint can move by
% - LoopCount: a counter for the number of iterations completed so far
% - MaxLoops: The upper limit on the number of iterations allowed
%
% outputs:
% - SearchPoint: the next work input value you should try
% - Step: the distance which can be stepped on the NEXT iteration. You
%   should store this value.
% - ReLoop: If this flag is set to 0, the upper level while loop will be
%   terminated
%
%
% example: you wish to use an audibility model to find the level at which
% you can be 50% confident of detecting the signal (-+ 1%). You choose to
% start with an input level of 20dB and take no more than 10 steps of
% 40,20,10 dB etc.
%
% You would implement this in the following way:
%
% TargetPercent = 0.5;
% CloseEnough = 0.01;
% Signal Level = 40;
% Step = 40;
% MaxLoops = 10;
% LoopCount = 0;
%
% while ReLoop = 1
%
%   LoopCount = LoopCount + 1;
%
%   Percent = RunAudibilityModel(SignalLevel)
%
%   [SignalLevel Step ReLoop] = iosr.auditory.binSearch(Percent, ...
%       TargetPercent, ...
%       CloseEnough, ...
%       SignalLevel, ...
%       Step, ...
%       LoopCount, ...
%       MaxLoops);
%
% end
%

%   Copyright 2016 University of Surrey.

function [SearchPoint, Step, ReLoop] = binSearch(Score, Target, CloseEnough, SearchPoint, Step, LoopCount, MaxLoops)

% input tests
% test input types
assert(isnumeric(Score) ...
    &isnumeric(Target) ...
    &isnumeric(CloseEnough) ...
    &isnumeric(SearchPoint) ...
    &isnumeric(Step) ...
    &isnumeric(LoopCount) ...
    &isnumeric(MaxLoops),'input arguments must be numeric!');

% Define next search point
if Score > Target
    SearchPoint = SearchPoint - Step;
else
    SearchPoint = SearchPoint + Step;
end

% half the distance of the next step
Step = Step / 2;

% if we got close enough to our taget, or it was the last iteration
if (abs((Score-Target))<CloseEnough)||(LoopCount>=MaxLoops)
    ReLoop = 0;
else
    ReLoop = 1;
end

end

author: IoSR-Surrey

repo: MatlabToolbox


text:
%{
 * Copyright (C) 2020-2030, The Regents of The University of Michigan.
 * All rights reserved.
 * This software was developed in the Biped Lab (https://www.biped.solutions/)
 * under the direction of Jessy Grizzle, grizzle@umich.edu. This software may
 * be available under alternative licensing terms; contact the address above.
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 * The views and conclusions contained in the software and documentation are those
 * of the authors and should not be interpreted as representing official policies,
 * either expressed or implied, of the Regents of The University of Michigan.
 *
 * AUTHOR: Bruce JK Huang (bjhuang[at]umich.edu)
 * WEBSITE: https://www.brucerobot.com/
%}

function data = splitPointsBasedOnRing(opts, point_cloud)
    data(opts.num_beams) = struct();
    data(opts.num_beams).points = [];
    data(opts.num_beams).point_with_I = [];
    for i = opts.pc_iter : opts.pc_iter+opts.num_scan
        total = 0;
        points = point_cloud(i, :, :);
        points = (reshape(points, size(points, 2),[]))';
        points = removeZeros(points);
        original_size = size(points,2);
        for n = 0:opts.num_beams-1
            if ~any(points(5,:)==n)
    %             disp("n")
    %             disp(n)
    %             disp("points")
    %             disp(points(:, points(5,:)==n))
    %             pause;
                continue
            end
            ring_point = points(:, points(5,:)==n);
            if size(ring_point,2) < opts.min_points
                original_size = original_size - size(ring_point,2);
                continue;
            end
            m = n+1;
            ring_point_homogeneous = [ring_point(1:3,:); ones(1,size(ring_point,2))];
            data(m).points = [data(m).points ring_point_homogeneous];
            data(m).point_with_I = [data(m).point_with_I ring_point];
            total= total+size(ring_point,2);
        end
        if total ~= original_size
            warning("inconsistent number of poins in splitPointsBasedOnRing.m");
            warning("split total: %i", total)
            warning("original total %i", size(points,2))
        end
    end
end

author: UMich-BipedLab


text:
function [result_features]=back_find_low_node_Nopara(result_dir,nSubj,k_times,nROI,w,cross_val,meth_FEX,meth_FS,varargin)
% This function aims to find features used in the classification, suitable
% for the PC, aHOFC, tHOFC brain network construction method. And it is also
% suitable for the LOOCV or the 10-fold cross validation. The found
% features can be applied on the visualization software to present the
% important brain regions or links.
% Input:
%         result_dir: the directory you want to store all the result files;
%         nSubj: number of subjects;
%         k_times: times of 10-fold cross validation;
%         nROI: number of ROIs
%         w: weight of each selected features;
%         cross_val: loocv or 10-fold;
%         meth_FEX: feature extraction method(connection coefficients or local clustering coefficients);
%         meth_FS: feature selection method (ttest,lasso,ttest+lasso);
%         varargin: the feature selection indexes;
%
% Output:
%        result_features: cell array, consisting of the averaged weight and occurrence of the features;
% Written by Zhen Zhou, zzstefan@email.unc.edu
% IDEA lab, https://www.med.unc.edu/bric/ideagroup
% Department of Radiology and BRIC, University of North Carolina at Chapel Hill

%clear all
% load SR_10_fold_middle.mat;
% load SR_loocv_middle.mat;
% cross_val='loocv';
% ttest_p=feature_index_ttest;
% midw_lasso=feature_index_lasso;
%load tHOFC_loocv_middle_clus.mat;
%load tHOFC_coef_LASSO_kfold.mat;
% load tHOFC_clus_LASSO_kfold.mat;
%load PC_loocv_coef_middle.mat;
% cross_val='10-fold';
%ttest_p=feature_index_ttest;
% midw_lasso=feature_index_lasso;
fprintf('Begin contributing feature identification\n');

switch meth_FS
    case 't-test + LASSO'
        ttest_p=varargin{1};
        midw_lasso=varargin{2};
        if strcmpi(cross_val,'10-fold')
            for i=1:size(ttest_p,1)
                for j=1:size(ttest_p,2)
                    tmp=find(ttest_p{i,j}<0.05);
                    tmp_b{i,j}=find(midw_lasso{i,j});
                    index=tmp(tmp_b{i,j});
                    if strcmpi(meth_FEX,'coef')
                        for k=1:length(index)
                            [first{i,j,k}(1),first{i,j,k}(2)]=find_elements(nROI,index(k));
                        end
                    elseif strcmpi(meth_FEX,'clus')
                        for k=1:length(index)
                            first{i,j,k}=index(k);
                        end
                    end
                end
            end
        elseif strcmpi(cross_val,'loocv')
            for i=1:length(ttest_p)
                tmp=find(ttest_p{i}<0.05);
                tmp_b{i}=find(midw_lasso{i});
                index=tmp(tmp_b{i});
                if strcmpi(meth_FEX,'coef')
                    for k=1:length(index)
                        [first{i,k}(1),first{i,k}(2)]=find_elements(nROI,index(k));
                    end
                elseif strcmpi(meth_FEX,'clus')
                    for k=1:length(index)
                        first{i,k}=index(k);
                    end
                end
            end
        end
    case 't-test'
        ttest_p=varargin{1};
        if strcmpi(cross_val,'10-fold')
            for i=1:size(ttest_p,1)
                for j=1:size(ttest_p,2)
                    tmp=find(ttest_p{i,j}<0.05);
                    if strcmpi(meth_FEX,'coef')
                        for k=1:length(tmp)
                            [first{i,j,k}(1),first{i,j,k}(2)]=find_elements(nROI,tmp(k));
                        end
                    elseif strcmpi(meth_FEX,'clus')
                        for k=1:length(tmp)
                            first{i,j,k}=tmp(k);
                        end
                    end
                end
            end
        elseif strcmpi(cross_val,'loocv')
            for i=1:length(ttest_p)
                tmp=find(ttest_p{i}<0.05);
                if strcmpi(meth_FEX,'coef')
                    for k=1:length(tmp)
                        [first{i,k}(1),first{i,k}(2)]=find_elements(nROI,tmp(k));
                    end
                elseif strcmpi(meth_FEX,'clus')
                    for k=1:length(tmp)
                        first{i,k}=tmp(k);
                    end
                end
            end
        end
    case 'LASSO'
        midw_lasso=varargin{1};
        if strcmpi(cross_val,'10-fold')
            for i=1:size(midw_lasso,1)
                for j=1:size(midw_lasso,2)
                    tmp=find(midw_lasso{i,j});
                    if strcmpi(meth_FEX,'coef')
                        for k=1:length(tmp)
                            [first{i,j,k}(1),first{i,j,k}(2)]=find_elements(nROI,tmp(k));
                        end
                    elseif strcmpi(meth_FEX,'clus')
                        for k=1:length(tmp)
                            first{i,j,k}=tmp(k);
                        end
                    end
                end
            end
        elseif strcmpi(cross_val,'loocv')
            for i=1:length(midw_lasso)
                tmp=find(midw_lasso{i});
                if strcmpi(meth_FEX,'coef')
                    for k=1:length(tmp)
                        [first{i,k}(1),first{i,k}(2)]=find_elements(nROI,tmp(k));
                    end
                elseif strcmpi(meth_FEX,'clus')
                    for k=1:length(tmp)
                        first{i,k}=tmp(k);
                    end
                end
            end
        end
end

%first=cellfun(@(x) x',first,'UniformOutput',false);
temp=first(:);

[a1,b,c]=unique(cellfun(@char,temp,'un',0));
lo=histc(c,1:max(c));
loo=lo(:)>1;
out=[temp(b(loo)),num2cell(lo(loo))];
out(1,:)=[];

All_link=out(find(cell2mat(out(:,2))>0),1);%find all links which exists in any one  LOOCV fold
if strcmpi(cross_val,'10-fold')
    for j=1:length(All_link)
        weight=[];
        temp_index=cell2mat(cellfun(@(x)isequal(x,All_link{j}),first,'un',0));
        ind=find(temp_index);
        [A,B,C]=ind2sub(size(temp_index),ind);
        for k=1:length(A)
            weight=[weight,w{A(k),B(k)}(C(k))];
        end
        all_weight(j)=mean(weight);
    end
elseif strcmpi(cross_val,'loocv')
    for j=1:length(All_link)
        weight=[];
        temp_index=cell2mat(cellfun(@(x)isequal(x,All_link{j}),first,'un',0));
        [A,B]=find(temp_index);
        for k=1:length(A)
            weight=[weight,w{A(k)}(B(k))];
        end
        all_weight(j)=mean(weight);
    end
end

if strcmpi(meth_FEX,'coef')
    matrix_1=zeros(nROI,nROI);
    matrix_2=zeros(nROI,nROI);
    for m=1:size(out,1)
        matrix_1(out{m,1}(1),out{m,1}(2))=all_weight(m);
        matrix_1(out{m,1}(2),out{m,1}(1))=all_weight(m);
        matrix_2(out{m,1}(1),out{m,1}(2))=out{m,2};
        matrix_2(out{m,1}(2),out{m,1}(1))=out{m,2};
    end
    result_features{1}=matrix_1;
    %result_features{2}=matrix_2;

elseif strcmpi(meth_FEX,'clus')
    matrix_1=zeros(1,nROI);
    matrix_2=zeros(1,nROI);
    for m=1:size(out,1)
        matrix_1(out{m,1})=all_weight(m);
        matrix_2(out{m,1})=out{m,2};
    end
    result_features{1}=matrix_1;
    %result_features{2}=matrix_2;
end
if strcmpi(cross_val,'loocv')
    matrix_2=matrix_2/nSubj;
elseif strcmpi(cross_val,'10-fold')
    matrix_2=matrix_2/(10*k_times);
end
result_features{2}=matrix_2;

if strcmpi(meth_FEX,'coef')
    figure('visible','off');
    subplot(1,2,1);
    imagesc(matrix_1);
    colormap jet
    colorbar
    axis square
    xlabel('ROI');
    ylabel('ROI');
    title('Averaged weight');

    subplot(1,2,2);
    imagesc(matrix_2);
    colormap jet
    colorbar
    axis square
    xlabel('ROI');
    ylabel('ROI');
    title('Normalized occurence');
    print(gcf,'-r1000','-dtiff',char(strcat(result_dir,'/result_features_weigthAndOccurence.tiff')));
end

fprintf('End contributing feature identification\n');
% node_matrix=eye(116,116);
% for i=1:length(NEW)
%     for j=1:length(NEW{i,3})
%         node_matrix(NEW{i,3}(j,1),NEW{i,3}(j,2))=1;
%         node_matrix(NEW{i,3}(j,2),NEW{i,3}(j,1))=1;
%     end
% end

author: zzstefan

repo: BrainNetClass


text:
% ********************************************************************
% * mex File compiling code for Random Forest (for linux)
% * mex interface to Andy Liaw et al.'s C code (used in R package randomForest)
% * Added by Abhishek Jaiantilal ( abhishek.jaiantilal@colorado.edu )
% * License: GPLv2
% * Version: 0.02
% ********************************************************************/
function compile_linux

    system('rm *.mexglx *.mexa64;');

    system('make clean;make mex;');

    %the fortran code makes it hard to NOT use the Makefile

repo: sr-metric


text:
%% Simulation Data
simu = simulationClass();                       % Initialize simulationClass
simu.simMechanicsFile = 'OSWEC.slx';            % Simulink Model File
simu.startTime = 0;                             % Simulation Start Time [s]
simu.rampTime = 100;                            % Wave Ramp Time [s]
simu.endTime=400;                               % Simulation End Time [s]
simu.dt = 0.1;                                  % Simulation Time-Step [s]

%% Wave Information
% Regular Waves
waves = waveClass('regular');                   % Initialize waveClass
waves.height = 2.5;                                  % Wave Height [m]
waves.period = 8;                                    % Wave Period [s]

%% Body Data
% Flap
body(1) = bodyClass('hydroData/oswec.h5');      % Initialize bodyClass for Flap
body(1).geometryFile = 'geometry/flap.stl';     % Geometry File
body(1).mass = 127000;                          % Mass [kg]
body(1).inertia = [1.85e6 1.85e6 1.85e6];  % Moment of Inertia [kg-m^2]

% Base
body(2) = bodyClass('hydroData/oswec.h5');      % Initialize bodyClass for Base
body(2).geometryFile = 'geometry/base.stl';     % Geometry File
body(2).mass = 'fixed';                         % Creates Fixed Body

%% PTO and Constraint Parameters
% Fixed
constraint(1)= constraintClass('Constraint1');  % Initialize constraintClass for Constraint1
constraint(1).location = [0 0 -10];                  % Constraint Location [m]

% Rotational PTO
pto(1) = ptoClass('PTO1');                      % Initialize ptoClass for PTO1
pto(1).stiffness = 0;                           % PTO Stiffness [Nm/rad]
pto(1).damping = 0;                                   % PTO Damping [Nsm/rad]
pto(1).location = [0 0 -8.9];                        % PTO Location [m]