text:
function Connect3D(p1,p2,option,pt)

h = plot3([p1(1) p2(1)],[p1(2) p2(2)],[p1(3) p2(3)],option);
set(h,'LineWidth',pt)

author: s-kajita

repo: IntroductionToHumanoidRobotics


text:
function [ t_stop, y_stop ] = p39_stop ( neqn )

%*****************************************************************************80
%
%% P39_STOP returns the stopping point for problem p39.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_STOP, Y_STOP(NEQN), the final data.
%
  y_stop = zeros ( neqn, 1 );

  a = p39_param ( 'GET', 'A', [] );
  b = p39_param ( 'GET', 'B', [] );
  if ( a ~= 0.0 )
    t_stop = ( b + 9.0 ) / a;
    y_stop = -3.0;
  elseif ( b ~= 0.0 )
    t_stop = 0.0;
    y_stop = 0.0;
  else
    t_stop = 10.0;
    y_stop = 0.1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%MAIN_cannon.m
%
% This script is used to find the optimal initial state for the cannon
% example problem, as described below:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Cannon Example:
%
% Find the optimal initial velocity for a projectile such that the
% resulting trajectory passes through a desired target.
%
% States:
%   x = horizontal position
%   y = vertical position
%   dx = horizontal speed
%   dy = vertical speed
%   ddx = horizontal acceleration
%   ddy = vertical acceleration
%
% Parameters:
%   c = quadratic drag coefficient
%   T = final simulation time
%
% Boundary Conditions:
%   x(0) = 0;
%   y(0) = 0;
%   x(T) = xf;
%   y(T) = yf;
%
% Objective Function:
%   J = dx(0)^2 + dy(0)^2;
%
% Dynamics Constraint:
%   ddx = -c*dx*sqrt(dx.*dx + dy.*dy);
%   ddy = -c*dy*sqrt(dx.*dx + dy.*dy) - 1;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc; clear;
%%%% Adjust the following parameters and see what happens.
%
% 1) The problem becomes trivial if the drag coefficient (c) is small, but
% numerically stiff if c is large (>0.4). Which methods can handle large
% values of c the best?
%
%
% 2) What do you think happens if you use more (or fewer) grid points?
%
% 3) What happens if you move the target away?
% Open:  TEST_cannonFeasibility.m   to see the answer.
%
% 4) For either large range or large damping, the problem becomes stiff. As
% a result, the integration methods (either implicit or explicit) start to
% fail near the beginning of the trajectory where the speed is very high.
% If you use a log spaced grid (rather than uniform), then all methods
% perform better.
%
% 5) You can speed up the optimizations by disabling diagnostics or making
% each animation shorter.
%

%Provide an initial guess for the solver
guess.initSpeed = 9.0;
guess.initAngle = (pi/180)*45;

% Set the target (assuming that the trajectory starts at x=0, y=0)
target.x = 6.0;
target.y = 0;

% Set up the parameters for the dynamics function:
param.dynamics.c = 0.4;  %Quadratic drag coefficient

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Set up the grid discretization for each method:
param.singleShooting.nGrid = 100;
param.multipleShooting.nSegment = 10;
param.multipleShooting.nSubStep = 5;
param.collocation.nSegment = 50;
param.gpops.nSegment = 10;

%%% Parameters for diagnostics (visualization only)
param.diagnostics.enable = true;   %Enable plotting and log iterations?
param.diagnostics.animationDuration = 5;  %(seconds) How long is the animation?
param.diagnostics.writeGif = true;   %Save animation to a gif?
param.diagnostics.gifPixelDim = [800,400];  %How big of a gif to make?
param.diagnostics.figNum.singleShooting = 10;
param.diagnostics.figNum.multipleShooting = 11;
param.diagnostics.figNum.collocation = 12;
param.diagnostics.figNum.gpops = 13;

% Use single shooting to find the solution:
soln.singleShooting = cannon_singleShooting(guess,target,param);

% Use multiple shooting to find the solution:
soln.multipleShooting = cannon_multipleShooting(guess,target,param);

% Use collocation to find the solution:
soln.collocation = cannon_collocation(guess,target,param);

% % Use GPOPS to solve the problem   (http://www.gpops2.com/)
% % - adaptive high-order collocation
% % - not included in Matlab. Must be purchased seperately.
% soln.gpops = cannon_gpops(guess,target,param);

author: MatthewPeterKelly

repo: dscTutorials


text:
function [grad]= B_hadamard(input_layers, curr_layer, future_layers)

input = input_layers{1}.a;
input2 = input_layers{2}.a;

future_grad = GetFutureGrad(future_layers, curr_layer);

[n1,n2,n3] = size(future_grad);

grad{1} = future_grad .* input2;
grad{2} = future_grad .* input;

end

author: singaxiong

repo: SignalGraph


text:
% COMPDSP - Display standard info figures for a data decomposition
%             Creates four figure windows showing: Component amplitudes,
%             scalp maps, activations and activation spectra.
% Usage:
%      >> compdsp(data,weights,locfile,[srate],[title],[compnums],[amps],[act]);
%
% Inputs:
%   data     = data matrix used to train the decomposition
%   weights  = the unmixing matrix (e.g., weights*sphere from RUNICA)
%
% Optional:
%   locfile  = 2-d electrode locations file (as in >> topoplot example)
%              {default|[]: default locations file given in icadefs.m
%   srate    = sampling rate in Hz {default|[]: as in icadefs.m}
%   title    = optional figure title text
%              {default|'': none}
%   compnums = optional vector of component numbers to display
%              {default|[] -> all}
%   amps     = all component amplitudes (from RUNICA)
%              {default|[]->recompute}
%   act      = activations matrix (from RUNICA)
%              {default|[]->recompute}
%
% Author: Scott Makeig, SCCN/INC/UCSD, La Jolla, 2000

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

% 02-01-01  replaced STD with RMS -sm
% 02-10-01  made srate optional -sm
% 01-25-02 reformated help & license -ad

function compdsp(data,unmix,locfile,srate,titl,compnums,amps,act)

minHz = 2; % frequency display limits
maxHz = 40;
ACTS_PER_EEGPLOT = 32;

%
%%%%%%%%%%%%%%%%%%%%%% Read and test arguments %%%%%%%%%%%%%%%%%%%%%%%%
%
icadefs % read BACKCOLOR, DEFAULT_SRATE

if nargin<2
  help compdsp
  return
end

chans = size(data,1);
frames = size(data,2);
ncomps = size(unmix,1);

if ncomps < 1
   error('Unmixing matrix must have at least one component');
end
if chans < 2
   error('Data must have at least two channels');
end
if frames < 2
   error('Data must have at least two frames');
end
if size(unmix,2) ~= chans
   error('Sizes of unmix and data do not match');
end
if nargin<3
 locfile=[];
end
if isempty(locfile)
   locfile = DEFAULT_ELOC; % from icsdefs.m
end
if ~exist(locfile)
   error('Cannot find electrode locations file');
end

if nargin<4
  srate = 0; % from icadefs
end
if isempty(srate) || srate==0
  srate = DEFAULT_SRATE; % from icadefs
end
if nargin<5
  titl = '';    % default - no title text
end
if isempty(titl)
  titl = '';
end
if nargin<6
  compnums = 0; % default - all components
end
if isempty(compnums)
  compnums = 0;
end
if nargin<7
  amps = NaN; % default - recompute amps
end
if isempty(amps)
  amps = NaN;
end
if ~isnan(amps) && length(amps) ~= ncomps
   error('Supplied amps does not match the number of unmixed components');
end

if compnums(1) == 0
  compnums = 1:ncomps;
end
if min(compnums) < 1
   error('Compnums must be positive');
end
if min(compnums) > ncomps
   error('Some compnum > number of components in unmix');
end

if nargin<8
   act = NaN; % default - recompute act
end

%
%%%%%%%%%%%%%%%%%%%%%%%%%%% Create plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

if chans == ncomps
  winv = inv(unmix);
elseif chans < ncomps
   error('More components than channels?');
else
  winv = pinv(unmix);
end

if isnan(act)
  fprintf('Computing activations ...')
  act = unmix(compnums,:)*data;
  fprintf('\n');
elseif size(act,2) ~= frames
  error('Supplied activations do not match data length');
elseif size(act,1) ~=  ncomps && size(act,1) ~= length(compnums)
  error('Number of supplied activations matrix does not match data or weights');
elseif size(act,1) ==  ncomps
  act = act(compnums,:); % cut down matrix to desired components
end

%
%%%%%%%%%%%%%%%%%%%%% I. Plot component amps %%%%%%%%%%%%%%%%%%%%%%%%%%%
%
pos = [40,520,550,400]; figure('Position',pos);
if isnan(amps)
  fprintf('Computing component rms amplitudes ');
  amps = zeros(1,length(compnums));
  for j=1:length(compnums)
     amps(j) = rms(winv(:,compnums(j)))*rms(act(j,:)');
     fprintf('.')
  end
  fprintf('\n');
else
  amps = amps(compnums); % truncate passed amps to desired compnums
end
plot(compnums,amps,'k');
hold on
plot(compnums,amps,'r^');
xl=xlabel('Component numbers');
yl=ylabel('RMS Amplitudes');
tl=title([titl ' Amplitudes']);
ax = axis;
axis([min(compnums)-1 max(compnums)+1 0 ax(4)]);

%
%%%%%%%%%%%%%%%%%%%%%%%%% II. Plot component maps %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
pos = [40,40,550,400]; figure('Position',pos);

fprintf('Plotting component scalp maps ...') % compmaps() may make multiple figures
compmap(winv,locfile,compnums,[titl ' Scalp Maps'],0,compnums);
fprintf('\n');

%
%%%%%%%%%%%%%%%%%%%%%% III. EEGPLOT activations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if frames/srate < 10
  dispsecs = ceil(frames/srate);
else
  dispsecs = 10; % defaults - display 10s data per screen
end
range = 0.8*max(max(act')-min(act'));

stact=1;
lact=ACTS_PER_EEGPLOT;
if lact>size(act,1)
   lact = size(act,1);
end

pos = [620,520,550,400]; figure('Position',pos);
while stact <= size(act,1)
  % eegplot(data,srate,spacing,eloc_file,windowlength,title,posn)
  eegplot(act(stact:lact,:),srate,range,compnums(stact:lact),...
               dispsecs,[titl ' Activations'],pos);
  pos = pos + [.02 .02 0 0];
  stact = stact+ACTS_PER_EEGPLOT;
  lact  = lact +ACTS_PER_EEGPLOT;
  if lact>size(act,1)
   lact = size(act,1);
  end
end

%
%%%%%%%%%%%%%%%%%%%%%%% IV. PLOTDATA spectra %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
pos = [620,40,550,400]; figure('Position',pos);

if frames > 2048
  windw = 512;
elseif frames > 1024
  windw = 256;
else
  windw = 128;
end
fprintf('Computing component spectra ')
for j = 1:length(compnums)
  % [Pxx,F] = PSD(X,NFFT,Fs,WINDOW,NOVERLAP)
  [spec,freqs] = psd(act(j,:),1024,srate,windw,ceil(windw*0.5));
  if ~exist('specs')
     specs = zeros(length(compnums),length(freqs));
  end
  specs(j,:) = spec';
  fprintf('.')
end
fprintf('\n');
specs = 10*log10(specs);

tmp = ceil(sqrt(length(compnums)));
tmp2 = ceil(length(compnums)/tmp);
for j=1:length(compnums)
  sbplot(tmp2,tmp,j)
  plot(freqs,specs(j,:))
  set(gca,'box','off')
  set(gca,'color',BACKCOLOR);
  ax=axis;
  axis([minHz,maxHz,ax(3),ax(4)]);
  tl = title(int2str(compnums(j)));
end
xl=xlabel('Frequency (Hz)');
yl=ylabel('Power (dB)');
set(gca,'YAxisLocation','right');
txl=textsc([titl ' Activation Spectra'],'title');
axcopy % pop-up axes on mouse click

% plottopo(specs,[tmp2 tmp],0,[2,70 min(min(specs)) max(max(specs))],...
%                 [titl ' Activation Power Spectra']);

function rmsval = rms(column)
   rmsval = sqrt(mean(column.*column)); % don't remove mean

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function g=comp_pchirp(L,n)
%COMP_PCHIRP  Compute periodic chirp
%   Usage:  g=comp_pchirp(L,n);
%
%   `pchirp(L,n)` returns a periodic, discrete chirp of length *L* that
%   revolves *n* times around the time-frequency plane in frequency. *n* must be
%   an integer number.
%
%   This is a computational routine. Do not call it unless you have
%   verified the input parameters.

%   AUTHOR : Peter L. Søndergaard
%   TESTING: OK
%   REFERENCE: OK

l= (0:L-1).';
X = mod(mod(mod(n*l,2*L).*l,2*L)*(L+1),2*L);
g = exp(pi*1i*X/L);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function pde = Lshapedata3
%% LSHAPEDATA data of Lshape Problem
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

pde = struct('f',0,'exactu',@exactu,'g_D',@exactu,'Du',@Du);

    function u = exactu(p) % exact solution
    x = p(:,1); y = p(:,2);
    r = sqrt(x.^2+y.^2);
    theta = atan2(p(:,2),p(:,1));
    theta = (theta>=0).*theta + (theta<0).*(theta+2*pi);
    u = r.^(2/3).*sin(2*theta/3);
    end

    function uprime = Du(p) % exact solution
    x = p(:,1); y = p(:,2);
    r = sqrt(x.^2+y.^2);
    theta = atan2(y,x);
    theta = (theta>=0).*theta + (theta<0).*(theta+2*pi);
    uprime(:,1) = 2/3*r.^(-1/3).*sin(2*theta/3).*x./r ...
                - 2/3*r.^(2/3).*cos(2*theta/3).*y./r.^2;
    uprime(:,2) = 2/3*r.^(-1/3).*sin(2*theta/3).*y./r ...
                + 2/3*r.^(2/3).*cos(2*theta/3).*x./r.^2;
    uprime(:,3) = 0;
    end
end

repo: ifem


text:
%%*******************************************************************
%% HSDHKMpred: Compute (dX,dy,dZ) for the H..K..M direction.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*******************************************************************

function [par,dX,dy,dZ,coeff,L,hRd] = ...
    HSDHKMpred(blk,At,par,rp,Rd,sigmu,X,Z,invZchol)

global schurfun schurfun_par
%%
%% compute HKM scaling
%%
Zinv = cell(size(blk,1),1);   dd = cell(size(blk,1),1);
gamx = cell(size(blk,1),1);   gamz = cell(size(blk,1),1);
for p = 1:size(blk,1)
    pblk = blk(p,:);
    n = sum(pblk{2});
    numblk = length(pblk{2});
    if strcmp(pblk{1},'l')
        Zinv{p} = 1./Z{p};
        dd{p} = X{p}./Z{p};
    elseif strcmp(pblk{1},'q')
        gaptmp  = qops(pblk,X{p},Z{p},1);
        gamz2   = qops(pblk,Z{p},Z{p},2);
        gamz{p} = sqrt(gamz2);
        Zinv{p} = qops(pblk,-1./gamz2,Z{p},4);
        dd{p}   = qops(pblk,gaptmp./gamz2,ones(n,1),4);
    elseif strcmp(pblk{1},'s')
        if (numblk == 1)
            Zinv{p} = Prod2(pblk,full(invZchol{p}),invZchol{p}',1);
        else
            Zinv{p} = Prod2(pblk,invZchol{p},invZchol{p}',1);
        end
    end
end
par.Zinv = Zinv; par.gamx = gamx; par.gamz = gamz; par.dd = dd;
%%
%% compute schur matrix
%%
m = par.m;
schur = sparse(m+2,m+2);
UU = [];  EE = [];
%%
for p = 1:size(blk,1)
    pblk = blk(p,:);
    if strcmp(pblk{1},'l')
        [schur,UU,EE] = schurmat_lblk(blk,At,par,schur,UU,EE,p,par.dd);
    elseif strcmp(pblk{1},'q');
        [schur,UU,EE] = schurmat_qblk(blk,At,par,schur,UU,EE,p,par.dd,par.Zinv,X);
    elseif strcmp(pblk{1},'s')
        if isempty(schurfun{p})
            schur = schurmat_sblk(blk,At,par,schur,p,X,par.Zinv);
        elseif ischar(schurfun{p})
            if ~isempty(par.permZ{p})
                Zpinv = Zinv{p}(par.permZ{p},par.permZ{p});
                Xp = X{p}(par.permZ{p},par.permZ{p});
            else
                Xp = X{p};
                Zpinv = Zinv{p};
            end
            schurtmp = feval( schurfun{p}, Xp,Zpinv,schurfun_par(p,:));
            schur = schur + schurtmp;
        end
    end
end
%%
%% compute rhs
%%
[rhs,EinvRc,hRd] = HSDHKMrhsfun(blk,At,par,X,Z,rp,Rd,sigmu);
%%
%% solve linear system
%%
par.addschur   = par.kap/par.tau;
schur(m+1,m+1) = schur(m+1,m+1) + par.kap/par.tau;
schur(m+2,m+2) = schur(m+2,m+2) + par.addschur;
[xx,coeff,L] = HSDlinsysolve(par,schur,UU,EE,par.Umat,rhs);
%%
%% compute (dX,dZ)
%%
[par,dX,dy,dZ] = HSDHKMdirfun(blk,At,par,Rd,EinvRc,X,xx);
%%*******************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% Version 1.000
%
% Code provided by Ruslan Salakhutdinov
%
% Permission is granted for anyone to copy, use, modify, or distribute this
% program and accompanying programs and documents for any purpose, provided
% this copyright notice is retained and prominently displayed, along with
% a note saying that the original programs are available from our
% web page.
% The programs and documents are distributed without any warranty, express or
% implied.  As the programs were written for research purposes only, they have
% not been tested to the degree that would be advisable in any important
% application.  All use of these programs is entirely at the user's own risk.

%%%%% Initialize biases of the base-rate model by ML %%%%%%%%%%%%%%%%%%%%%
 [numcases numdims numbatches]=size(batchdata);
 count_int = zeros(numdims,1);
 for batch=1:numbatches
    xx = sum(batchdata(:,:,batch));
    count_int = count_int + xx';
 end

 lp=5;
 p_int = (count_int+lp*numbatches)/(numcases*numbatches+lp*numbatches);
 log_base_rate = log( p_int) - log(1-p_int);

repo: DeepFi


text:
function output = yalmip_default_output;
output.Primal    = [];
output.Dual      = [];
output.Slack     = [];
output.problem = 0;
output.infostr = yalmiperror(0);
output.solverinput  = [];
output.solveroutput = [];
output.solvertime   = 0;


text:
function test_suite = test_intersectPlaneSphere
%TESTINTERSECTPLANESPHERE  One-line description here, please.
%
%   output = testIntersectPlaneSphere(input)
%
%   Example
%   testIntersectPlaneSphere
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-06-21,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testOx(testCase) %#ok<*DEFNU>
% Plane normal to direction [1 0 0]

center = [10 20 30];
radius = 50;
sphere = [center radius];

plane = createPlane(center, [6 0 0]);

circle = intersectPlaneSphere(plane, sphere);

exp = [center radius 90 0 0];
testCase.assertEqual(exp, circle);

function testOxShifted(testCase)
% Plane normal to direction [1 0 0]

center = [10 20 30];
radius = 50;
sphere = [center radius];

center2 = center + [radius/2 0 0];
plane = createPlane(center2, [6 0 0]);

circle = intersectPlaneSphere(plane, sphere);

radius2 = radius * sqrt(3) / 2;
exp = [center2 radius2 90 0 0];
testCase.assertEqual(exp, circle);

function testOy(testCase)
% Plane normal to direction [0 1 0]

center = [10 20 30];
radius = 50;
sphere = [center radius];

plane = createPlane(center, [0 6 0]);

circle = intersectPlaneSphere(plane, sphere);

exp = [center radius 90 90 0];
testCase.assertEqual(exp, circle, 'AbsTol', .01);

function testOyShifted(testCase)
% Plane normal to direction [0 1 0]

center = [10 20 30];
radius = 50;
sphere = [center radius];

center2 = center + [0 radius/2 0];
plane = createPlane(center2, [0 6 0]);

circle = intersectPlaneSphere(plane, sphere);

radius2 = radius * sqrt(3) / 2;
exp = [center2 radius2 90 90 0];
testCase.assertEqual(exp, circle, 'AbsTol', .01);

function testOz(testCase)
% Plane normal to direction [0 1 0]

center = [10 20 30];
radius = 50;
sphere = [center radius];

plane = createPlane(center, [0 0 6]);

circle = intersectPlaneSphere(plane, sphere);

exp = [center radius 0 0 0];
testCase.assertEqual(exp, circle, 'AbsTol', .01);

function testOzShifted(testCase)
% Plane normal to direction [0 0 1]

center = [10 20 30];
radius = 50;
sphere = [center radius];

center2 = center + [0 0 radius/2];
plane = createPlane(center2, [0 0 6]);

circle = intersectPlaneSphere(plane, sphere);

radius2 = radius * sqrt(3) / 2;
exp = [center2 radius2 0 0 0];
testCase.assertEqual(exp, circle, 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function pass = test_sinh( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e4*pref.techPrefs.chebfuneps;

% Example 1
f = sinh(ballfun(@(x,y,z)y));
exact = ballfun(@(x,y,z)sinh(y));
pass(1) = norm( f - exact ) < tol;

if (nargout > 0)
    pass = all(pass(:));
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function linplus_test0196 ( )

%*****************************************************************************80
%
%%  TEST0196 tests R8VEC_TO_R8GE, R8GE_TO_R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  m = 4;
  n = 6;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0196\n' );
  fprintf ( 1, '  For a general matrix,\n' );
  fprintf ( 1, '  R8VEC_TO_R8GE converts a real vector to a R8GE matrix.\n' );
  fprintf ( 1, '  R8GE_TO_R8VEC converts a R8GE matrix to a real vector.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix rows M =    %d\n', m );
  fprintf ( 1, '  Matrix columns N = %d\n', n );

  a = r8ge_indicator ( m, n );

  r8ge_print ( m, n, a, '  The R8GE indicator matrix:' );

  x = r8ge_to_r8vec ( m, n, a );

  k = 0;
  for j = 1 : n
    for i = 1 : m
      k = k + 1;
      fprintf ( 1, '%4d  %4d  %4d  %14f\n', i, j, k, x(k) );
    end
  end

  a = r8vec_to_r8ge ( m, n, x );

  r8ge_print ( m, n, a, '  The recovered R8GE indicator matrix:' );

  return
end


text:
classdef ARMOEA < ALGORITHM
% <multi/many> <real/integer/label/binary/permutation> <constrained/none>
% Adaptive reference points based multi-objective evolutionary algorithm

%------------------------------- Reference --------------------------------
% Y. Tian, R. Cheng, X. Zhang, and Y. Jin, An indicator-based
% multiobjective evolutionary algorithm with reference point adaptation for
% better versatility, IEEE Transactions on Evolutionary Computation, 2018,
% 22(4): 609-622.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Generate the sampling points and random population
            Population = Problem.Initialization();
            W          = UniformPoint(Problem.N,Problem.M);
            [Archive,RefPoint,Range] = UpdateRefPoint(Population(all(Population.cons<=0,2)).objs,W,[]);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = MatingSelection(Population,RefPoint,Range);
                Offspring  = OperatorGA(Problem,Population(MatingPool));
                [Archive,RefPoint,Range] = UpdateRefPoint([Archive;Offspring(all(Offspring.cons<=0,2)).objs],W,Range);
                [Population,Range]       = EnvironmentalSelection([Population,Offspring],RefPoint,Range,Problem.N);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [ml] = yd32ml(yd3)
% Convert volume from cubic yards to milliliters.
% Chad Greene 2012
ml = yd3*764554.85798;

author: Sable

repo: mcbench-benchmarks


text:
function x = log_normal_cdf_inv ( cdf, a, b )

%*****************************************************************************80
%
%% LOG_NORMAL_CDF_INV inverts the Lognormal CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 1999
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%    0.0 <= CDF <= 1.0.
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < B.
%
%    Input, real X, the corresponding argument.
%
  if ( cdf < 0.0 | 1.0 < cdf )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LOG_NORMAL_CDF_INV - Fatal error!\n' );
    fprintf ( 1, '  CDF < 0 or 1 < CDF.\n' );
    error ( 'LOG_NORMAL_CDF_INV - Fatal error!' );
  end

  logx = normal_cdf_inv ( cdf, a, b );

  x = exp ( logx );

  return
end


text:
function [x,M] = spm_x_nmm(P)
% initialises a state structure for a mean field model
% FORMAT [x,M] = spm_x_nmm(P)
%
% P - parameter structure
% M - model structure
%
% x        - array of states
% x(i,j,k) - k-th state of j-th population on i-th source
%
%   population: 1 - excitatory spiny stellate cells (input cells)
%               2 - inhibitory interneurons
%               3 - excitatory pyramidal cells      (output cells)
%
%        state: 1 V  - voltage
%               2 gE - conductance (excitatory)
%               3 gI - conductance (inhibitory)
%
% M - model structure
%
% see also: spm_x_mfm
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_x_nmm.m 2393 2008-10-23 14:58:50Z karl $

% get initialisation from full mean-field model
%==========================================================================
[x M] = spm_x_mfm(P);

% remove dispersion and fix the covariance of the states (Cx)
%--------------------------------------------------------------------------
M.x   = x{1};
M.Cx  = x{2}(:,:,1,1);
x     = x{1};

author: spm


text:
function n = p01_n ( )

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

  return
end


text:
function [D] = spm_eeg_invert_prepro(D,val)
%
%% Preprocessing for inversion stage.
%% includes spatial and temporal dimension reduction
%
%% this version only handles single subject single modality data
%% the removal of many scaling factors makes it easier to compare between forward models
% ReML inversion of multiple forward models for EEG-MEG
% FORMAT [D] = spm_eeg_invert_classic(D)
% ReML estimation of regularisation hyperparameters using the
% spatiotemporal hierarchy implicit in EEG/MEG data
%
% Requires:
% D{i}.inv{val}.inverse:
%
%     inverse.modality - modality to use in case of multimodal datasets
%
%     inverse.trials - D.events.types to invert
%     inverse.type   - 'GS' Greedy search on MSPs
%                      'ARD' ARD search on MSPs
%                      'MSP' GS and ARD multiple sparse priors
%                      'LOR' LORETA-like model
%                      'IID' minimum norm
%                      'EBB' for empirical bayes beamformer

%    inverse.priors{} - a cell array of anatomical and functional priors to be
%    considered
%
%     inverse.woi    - time window of interest ([start stop] in ms)
%     inverse.lpf    - band-pass filter - low frequency cut-off (Hz)
%     inverse.hpf    - band-pass filter - high frequency cut-off (Hz)
%     inverse.Han    - switch for Hanning window
%     inverse.xyz    - (n x 3) locations of spherical VOIs
%     inverse.rad    - radius (mm) of VOIs
%
%     inverse.Nm     - maximum number of channel modes
%     inverse.Nr     - maximum number of temporal modes
%     inverse.Np     - number of sparse priors per hemisphere
%     inverse.smooth - smoothness of source priors (0 to 1)
%     inverse.Na     - number of most energetic dipoles
%     inverse.sdv    - standard deviations of Gaussian temporal correlation
%     inverse.Qe     - any sensor error components (e.g. empty-room data)
%     inverse.Qe0     - minimum amount of sensor noise power relative to
%                        signal eg 0.1 would correspond to power SNR of 10.0
%     inverse.A       - predefined spatial modes (Nchans*Nmodes) to project
%                       sensor data through
%
% Evaluates:
%
%     inverse.M      - MAP projector (reduced)
%     inverse.J{i}   - Conditional expectation (i conditions) J = M*U*Y
%     inverse.L      - Lead field (reduced UL := U*L)
%     inverse.qC     - spatial covariance
%     inverse.qV     - temporal correlations
%     inverse.T      - temporal projector
%     inverse.U(j)   - spatial projector (j modalities) - derived from data
%     inverse.A      - pre-specified spatial projector
%     inverse.Y{i}   - reduced data (i conditions) UY = UL*J + UE
%     inverse.Is     - Indices of active dipoles
%     inverse.It     - Indices of time bins
%     inverse.Ic{j}  - Indices of good channels (j modalities)
%     inverse.Nd     - number of dipoles
%     inverse.pst    - peristimulus time
%     inverse.dct    - frequency range
%     inverse.F      - log-evidence
%     inverse.VE     - variance explained in spatial/temporal subspaces (%)
%     inverse.R2     - variance in subspaces accounted for by model (%)
%     inverse.scale  - scaling of data for each of j modalities
%__________________________________________________________________________
%
% This version is for single subject single modality analysis and therefore
% contains none of the associated scaling factors.
% No symmetric priors are used in this implementation (just single patches)
% There is an option for a Beamforming prior : inversion type 'EBB'
% also added new beamforming method- using GS rather than ARD- from Juan
% David Martinez Vargas 'EBBgs'
%
% The code was used in
% Lopez, J. D., Penny, W. D., Espinosa, J. J., Barnes, G. R. (2012).
% A general Bayesian treatment for MEG source reconstruction incorporating
% lead field uncertainty.
% Neuroimage 60(2), 1194-1204 doi:10.1016/j.neuroimage.2012.01.077.
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Jose David Lopez, Gareth Barnes, Vladimir Litvak
% $Id: spm_eeg_invert_prepro.m 7017 2017-02-15 12:50:58Z karl $

Nl = length(D);

if Nl>1
    error('function only defined for a single subject');
else
    D=D{1};
end

% D - SPM data structure
%==========================================================================
if nargin > 1
    D.val = val;
elseif ~isfield(D, 'val')
    D.val = 1;
end

val=D.val;

inverse   = D.inv{val}.inverse;

% forward model
%--------------------------------------------------------------------------

% defaults
%--------------------------------------------------------------------------

try, type = inverse.type;   catch, type = 'GS';     end
try, s    = inverse.smooth; catch, s    = 0.6;      end
try, Np   = inverse.Np;     catch, Np   = 256;      end
try, Nr   = inverse.Nr;     catch, Nr   = 16;       end % requested number of temporal modes, could be changed depending on svd
try, xyz  = inverse.xyz;    catch, xyz  = [0 0 0];  end
try, rad  = inverse.rad;    catch, rad  = 128;      end
try, hpf  = inverse.hpf;    catch, hpf  = 48;       end % need to one day put these the correct way round
try, lpf  = inverse.lpf;    catch, lpf  = 0;        end
try, sdv  = inverse.sdv;    catch, sdv  = 4;        end
try, Han  = inverse.Han;    catch, Han  = 1;        end
try, woi  = inverse.woi;    catch, woi  = [];       end
try, Nm   = inverse.Nm;     catch, Nm   = [];       end
try, Nt   = inverse.Nt;     catch, Nt   = [];       end % fixed number of temporal modes
try, Ip   = inverse.Ip;     catch, Ip   = [];       end
try, QE   = inverse.QE;     catch, QE   = 1;        end %  empty room noise measurement
try, Qe0  = inverse.Qe0;    catch, Qe0  = exp(-5);  end % set noise floor at 1/100th signal power i.e. assume amplitude SNR of 10
try, inverse.A;             catch, inverse.A   = [];end % orthogonal channel modes

% get specified modalities to invert (default to all)
%--------------------------------------------------------------------------
modalities = D.inv{val}.forward.modality;    % MEG in this case
Nmax       = 16;                             % max number of temporal modes

% check lead fields and get number of dipoles (Nd) and channels (Nc)
%==========================================================================
fprintf('Checking leadfields')
[L,D] = spm_eeg_lgainmat(D);                 % Generate/load lead field
Nd    = size(L,2);

% Check gain or lead-field matrices
%------------------------------------------------------------------
if size(modalities,1)>1,
    error('not defined for multiple modalities');
end;

Ic{1}  = setdiff(D.indchantype(modalities), badchannels(D));
Nd    = size(L,2);      % Number of dipoles

%==========================================================================
% Spatial projectors (adjusting for different Lead-fields)
%==========================================================================

fprintf('Optimising and aligning spatial modes ...\n')

% eliminate low SNR spatial modes
%------------------------------------------------------------------
%disp('FIXING FULL RANK A');
%inverse.A=eye(length(Ic)); Nm=size(inverse.A,1);

if isempty(inverse.A), % no spatial modes prespecified
    if isempty(Nm),    % number of modes not specifiedd
        [U,ss,vv] = spm_svd((L*L'),exp(-16));
        A         = U';                % spatial projector A
        UL        = A*L;

    else % number of modes pre-specified
        [U,ss,vv]    = spm_svd((L*L'),0);
        if length(ss)<Nm,
            disp('number available');
            length(ss)
            error('Not this many spatial modes in lead fields');
        end
        ss    = ss(1:Nm);
        disp('using preselected number spatial modes !');
        A     = U(:,1:Nm)';                 % spatial projector A
        UL    = A*L;
    end;
else %% U was specified in input
    disp('Using pre-specified spatial modes');
    if isempty(Nm),
        error('Need to specify number of spatial modes if U is prespecified');
    end
    A  = inverse.A;
    UL = A*L;
end;

Nm     = size(UL,1);         % Number of spatial projectors

clear ss vv

% Report
%----------------------------------------------------------------------
fprintf('Using %d spatial modes',Nm)

%==========================================================================
% Temporal projector
%==========================================================================
AY    = {};                                      % pooled response for MVB
AYYA  = 0;                                       % pooled response for ReML

% Time-window of interest
%----------------------------------------------------------------------

if isempty(woi)
    w  = 1000*[min(D.time) max(D.time)];
else
    w  = woi; %% in milliseconds
end;

It     = (w/1000 - D.timeonset)*D.fsample + 1;
It     = max(1,It(1)):min(It(end), length(D.time));
It     = fix(It);
disp(sprintf('Number of samples %d',length(It)));

% Peristimulus time
%----------------------------------------------------------------------
pst    = 1000*D.time;                   % peristimulus time (ms)
pst    = pst(It);                       % windowed time (ms)
dur    = (pst(end) - pst(1))/1000;      % duration (s)
dct    = (It - It(1))/2/dur;            % DCT frequencies (Hz)
Nb     = length(It);                    % number of time bins

% Serial correlations
%----------------------------------------------------------------------
K      = exp(-(pst - pst(1)).^2/(2*sdv^2)); %% sdv set to 4 by default
K      = toeplitz(K);
qV     = sparse(K*K'); %% Samples* samples covariance matrix- assumes smooth iid

% Confounds and temporal subspace
%----------------------------------------------------------------------
T      = spm_dctmtx(Nb,Nb);         % use plot(T) here
j      = find( (dct >= lpf) & (dct <= hpf) ); %% THis is the wrong way round but leave for nowfor compatibility with spm_eeg_invert
T      = T(:,j);                    % Apply the filter to discrete cosines
dct    = dct(j);                    % Frequencies accepted

% Hanning window
%----------------------------------------------------------------------
if Han
    W  = sparse(1:Nb,1:Nb,spm_hanning(Nb)); %% use hanning unless specified
else
    W  = 1;
end

% get trials or conditions
%----------------------------------------------------------------------
try
    trial = D.inv{D.val}.inverse.trials;
catch
    trial = D.condlist;
end
Ntrialtypes = length(trial);

% get temporal covariance (Y'*Y) to find temporal modes
%======================================================================
YY    = 0;
N     = 0;
badtrialind=D.badtrials;
for j = 1:Ntrialtypes,              % pool over conditions
    c     = D.indtrial(trial{j});   % and trials
    c     = setxor(c,badtrialind);  % ignore bad trials
    Nk    = length(c);
    for k = 1:Nk
        Y  = A*D(Ic{1},It,c(k));

        YY = YY + Y'*Y;
        N  = N + 1;
    end
end
YY  = YY./N;

% Apply any Hanning and filtering
%------------------------------------------------------------------
YY  = W'*YY*W;     % Hanning
YTY = T'*YY*T;     % Filter

%======================================================================
if isempty(Nt), %% automatically assign appropriate number of temporal modes
    [U E]  = spm_svd(YTY,exp(-8));          % get temporal modes
    if isempty(U), %% fallback
        warning('nothing found using spm svd, using svd');
        [U E]  = svd(YTY);          % get temporal modes
    end;
    E      = diag(E)/trace(YTY);            % normalise variance
    Nr     = min(length(E),Nmax);           % number of temporal modes
    Nr=max(Nr,1); %% use at least one mode
else %% use predefined number of modes
    [U E]  = svd(YTY);          % get temporal modes
    E      = diag(E)/trace(YTY);            % normalise variance
    disp('Fixed number of temporal modes');
    Nr=Nt;
end;

V      = U(:,1:Nr);                     % temporal modes
VE     = sum(E(1:Nr));                  % variance explained

fprintf('Using %i temporal modes, ',Nr)
fprintf('accounting for %0.2f percent average variance\n',full(100*VE))

% projection and whitening
%----------------------------------------------------------------------
S      = T*V;                           % temporal projector
Vq     = S*pinv(S'*qV*S)*S';            % temporal precision

%disp('FIXIN TEMP PROJECTOR TO FULL RANK');
%S=eye(size(YY,1));

% get spatial covariance (Y*Y') for Gaussian process model
%======================================================================
% loop over Ntrialtypes trial types
%----------------------------------------------------------------------
UYYU    = 0;
AYYA    = 0;
Nn      = 0;                             % number of samples
AY      = {};
Ntrials = 0;

for j = 1:Ntrialtypes,
    UY{j} = sparse(0);
    c     = D.indtrial(trial{j});
    c     = setxor(c,badtrialind); %% ignore bad trials
    Nk    = length(c);

    % loop over epochs
    %------------------------------------------------------------------
    for k = 1:Nk

        % stack (scaled aligned data) over modalities
        %--------------------------------------------------------------
        Y = D(Ic{1},It,c(k))*S; %% in temporal subspace
        Y = A*Y; %%  in spatial subspace

        % accumulate first & second-order responses
        %--------------------------------------------------------------
        Nn      = Nn + Nr;             % number of samples

        YY      = Y*Y';                % and covariance
        Ntrials = Ntrials+1;

        % accumulate statistics (subject-specific)
        %--------------------------------------------------------------
        UY{j}   = UY{j} + Y;           % condition-specific ERP
        UYYU    = UYYU + YY;           % subject-specific covariance

        % and pool for optimisation of spatial priors over subjects
        %--------------------------------------------------------------
        AY{end + 1} = Y;               % pooled response for MVB
        AYYA        = AYYA  + YY;      % pooled response for ReML

    end
end
AY  = spm_cat(AY);     %% goes to MVB/GS algorithm
ID  = spm_data_id(AY); %% get a unique ID for these filtered data

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

inverse.AY=AY; %% concatenated data in reduced spatial modes
inverse.AYYA=AYYA;%% sensor level covariance of all trials, conditions and samples (containing Nn temporal modes in total)

inverse.UY     = UY;                    % ERP data (reduced)
inverse.L      = UL;                   % Lead-field (reduced)
inverse.Nn     = Nn; %% number of independent samples (temporal modes* trials* conditions)
inverse.qV     = Vq;                   % temporal correlations
inverse.T      = S;                    % temporal projector
inverse.U      = {A};                    % spatial projector
inverse.Ic     = Ic; %% good channels
inverse.It     = It; %% time indices
inverse.Nd     = Nd;                   % number of dipoles
inverse.pst    = pst;                  % peristimulus time
inverse.dct    = dct;                  % frequency range
inverse.ID     = ID;                   % data ID
inverse.woi    = w;                    % time-window inverted

inverse.modality = modalities;         % modalities inverted

% save in struct
%----------------------------------------------------------------------
D.inv{val}.inverse = inverse;
D.inv{val}.method  = 'Imaging';

disp('Done invert prepro');

return


text:
% VL_BENCH_IMREADJPEG   Evaluates the speed of imreadjpeg

numThreads = 4 ;
base = 'data/bench-imreadjpeg' ;

files = {} ;
files = dir(fullfile(base,'*.jpg')) ;
files = fullfile(base, {files.name}) ;
if numel(files) > 256, files = files(1:256) ; end

for preallocate = [true, false]
  opts={'verbose','verbose', 'preallocate', preallocate} ;
  for t=1:4
    % simple read
    fprintf('direct read single thread\n') ;
    clear ims ;
    tic  ;
    ims = vl_imreadjpeg(files, 'numThreads', 1, opts{:}) ;
    directSingle(t) = toc ;
    fprintf('   done\n') ;
    pause(1) ;

    % simple read
    fprintf('direct read multi thread\n') ;
    clear ims ;
    tic  ;
    ims = vl_imreadjpeg(files, 'numThreads', numThreads, opts{:}) ;
    direct(t) = toc ;
    fprintf('   done\n') ;
    pause(1) ;

    % threaded read
    fprintf('issue prefetch\n') ;
    tic ;
    vl_imreadjpeg(files, 'prefetch', opts{:}) ;
    prefetch(t) = toc ;
    fprintf('   done [pause 6]\n') ;
    pause(6)

    fprintf('prefetched read\n') ;
    clear ims_ ; % do not accoutn for the time requried to delete this
    tic ;
    ims_ = vl_imreadjpeg(files, opts{:}) ;
    indirect(t) = toc ;
    pause(1) ;
  end

  n = numel(ims) ;
  fprintf('** test results preallcoate %d\n', preallocate) ;
  fprintf('\tsingle tread: %.1f pm %.1f\n', mean(n./directSingle), std(n./directSingle)) ;
  fprintf('\t%d threads: %.1f pm %.1f\n', numThreads, mean(n./direct), std(n./direct)) ;
  fprintf('\tissue prefetch: %.1f pm %.1f\n', mean(n./prefetch), std(n./prefetch)) ;
  fprintf('\tretrieve prefetched: %.1f pm %.1f\n', mean(n./indirect), std(n./indirect)) ;
  fprintf('\n\n') ;
end

return

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function writeDoublePendulumDynamics(f,M)

%This function writes the dynamics file for the double pendulum.

filename = 'doublePendulumDynamics';

comments{1} = ['DZ = ' upper(filename) '(T,Z,P)'];
comments{2} = ' ';
comments{3} = 'FUNCTION:  This function computes the dynamics of a double';
comments{4} = '    pendulum, and is designed to be called from ode45. The';
comments{5} = '    model allows for arbitrary mass and inertia for each';
comments{6} = '    link, but no friction or actuation';
comments{7} = ' ';
comments{8}  = 'INPUTS: ';
comments{9}  = '    t = time. Dummy input for ode45. Not used.';
comments{10} = '    z = [4xn] matrix of states.';
comments{12} = '    P = struct of parameters';
comments{13} = 'OUTPUTS: ';
comments{14} = '    dz = [4xn] matrix of state derivatives';
comments{15} = ' ';
comments{16} = 'NOTES:';
comments{17} = ['    This file was automatically generated by ' mfilename];

params{1} = {'m1','link one mass'};
params{2} = {'m2','link two mass'};
params{3} = {'g ','gravity'};
params{4} = {'l1','link one length'};
params{5} = {'l2','link two length'};
params{6} = {'I1','link one moment of inertia about its center of mass'};
params{7} = {'I2','link two moment of inertia about its center of mass'};
params{8} = {'d1','distance between link one center of mass and parent joint'};
params{9} = {'d2','distance between link two center of mass and parent joint'};

states{1} = {'th1','link one absolute angle'};
states{2} = {'dth1','link one angular rate'};
states{3} = {'th2','link two absolute angle'};
states{4} = {'dth2','link two angular rate'};

dstates{1} = {'dth1','derivative of link one absolute angle'};
dstates{2} = {'ddth1','derivative of link one angular rate'};
dstates{3} = {'dth2','derivative of link two absolute angle'};
dstates{4} = {'ddth2','derivative of link two angular rate'};

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                               write file                                %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
fid = fopen([filename '.m'],'w');

fprintf(fid, ['function dz = ' filename '(~,z,P) \n']);

for i=1:length(comments)
    fprintf(fid,['%%' comments{i} '\n']);
end
fprintf(fid,'\n');

for i=1:length(params)
    fprintf(fid,[params{i}{1} ' = P.' params{i}{1} '; %%' params{i}{2} '\n']);
end
fprintf(fid,'\n');

for i=1:length(states)
    fprintf(fid,[states{i}{1} ' = z(' num2str(i) ',:); %%' states{i}{2} '\n']);
end
fprintf(fid,'\n');

for i=1:2;
    fprintf(fid,['f' num2str(i) ' = ' vectorize(char(f(i))) ';\n']);
end
fprintf(fid,'\n');

for i=1:2
    for j=1:2
        fprintf(fid,['M' num2str(i) num2str(j) ' = ' vectorize(char(M(i,j))) ';\n']);
    end
end
fprintf(fid,'\n');

fprintf(fid,'D = M11.*M22 - M12.*M21;\n\n');

fprintf(fid,'ddth1 = (f2.*M12 - f1.*M22)./D;\n');
fprintf(fid,'ddth2 = -(f2.*M11 - f1.*M21)./D;\n');
fprintf(fid,'\n');

fprintf(fid,'dz = [...\n');
for i=1:length(dstates)
    fprintf(fid,['    ' dstates{i}{1} '; %%' dstates{i}{2} '\n']);
end
fprintf(fid,'];\n');

fprintf(fid,'end \n');

fclose(fid);
end


text:
function val=doubleFactorial(x)
%%DOUBLEFACTORIAL Evaluate x!!, which is the double factorial (also known
%                 as semifactorial) of x. If x is odd, then the double
%                 factorial is defined as x*(x-2)*(x-4)*...*5*3*1. If x is
%                 even, then it is x(x-2)*(x-4)*...*6*4*2. The double
%                 factorials of 0 and -1 are defined to be 1.
%
%INPUTS: x A scalar or matrix, where the double factorial of each
%          element should be evaluated. The values must be positive
%          integers, or 0 or -1.
%
%OUTPUTS: val The double factorial of each value of x.
%
%The gamma function/factorial representations of the double factorial for
%even and odd values are taken from [1].
%
%Values <=50 are tabulated. For other values, the functions are implemented
%using logarithmic functions so as to minimize the effects of overflow with
%intermediate results.
%
%REFERENCES:
%[1] Weisstein, Eric W. "Double Factorial." From MathWorld--A Wolfram Web
%    Resource. http://mathworld.wolfram.com/DoubleFactorial.html
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numEl=length(x(:));

val=zeros(size(x));
for curEl=1:numEl
    xCur=x(curEl);

    if(xCur<=50)
        %If we can use a tabulated formula.
        if(xCur<-1)
            error('xCur must be >=-1.')
        end
        table=[1,1,1,2,3,8,15,48,105,384,945,3840,10395,46080,135135,645120,2027025,10321920,34459425,185794560,654729075,3715891200,13749310575,81749606400,316234143225,1961990553600,7905853580625,51011754393600,213458046676875,1428329123020800,6190283353629375,42849873690624000,191898783962510625,1371195958099968000,6332659870762850625,46620662575398912000,221643095476699771875,1678343852714360832000,8200794532637891559375,63777066403145711616000,319830986772877770815625,2551082656125828464640000,13113070457687988603440625,107145471557284795514880000,563862029680583509947946875,4714400748520531002654720000,25373791335626257947657609375,216862434431944426122117120000,1192568192774434123539907640625,10409396852733332453861621760000,58435841445947272053455474390625,520469842636666622693081088000000];

        val(curEl)=table(xCur+2);
    else
        if(xCur==0||xCur==-1)
            val(curEl)=1;
        elseif(mod(xCur,2)==0)%If it is divisible by 2.
            k=xCur/2;
            val(curEl)=exp(k*log(2)+gammaln(k+1));
        else%Assume that it is odd.
            n=(xCur+1)/2;
            val(curEl)=exp(gammaln(xCur+2)-(n*log(2)+gammaln(n+1)));
        end
    end
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function x = emailFeatures(word_indices)
    %% EMAILFEATURES takes in a word_indices vector and produces a feature vector
    %from the word indices
    %   x = EMAILFEATURES(word_indices) takes in a word_indices vector and
    %   produces a feature vector from the word indices.

    % Total number of words in the dictionary
    n = 1899;

    % You need to return the following variables correctly.
    x = zeros(n, 1);

    % ====================== YOUR CODE HERE ======================
    % Instructions: Fill in this function to return a feature vector for the
    %               given email (word_indices). To help make it easier to
    %               process the emails, we have have already pre-processed each
    %               email and converted each word in the email into an index in
    %               a fixed dictionary (of 1899 words). The variable
    %               word_indices contains the list of indices of the words
    %               which occur in one email.
    %
    %               Concretely, if an email has the text:
    %
    %                  The quick brown fox jumped over the lazy dog.
    %
    %               Then, the word_indices vector for this text might look
    %               like:
    %
    %                   60  100   33   44   10     53  60  58   5
    %
    %               where, we have mapped each word onto a number, for example:
    %
    %                   the   -- 60
    %                   quick -- 100
    %                   ...
    %
    %              (note: the above numbers are just an example and are not the
    %               actual mappings).
    %
    %              Your task is take one such word_indices vector and construct
    %              a binary feature vector that indicates whether a particular
    %              word occurs in the email. That is, x(i) = 1 when word i
    %              is present in the email. Concretely, if the word 'the' (say,
    %              index 60) appears in the email, then x(60) = 1. The feature
    %              vector should look like:
    %
    %              x = [ 0 0 0 0 1 0 0 0 ... 0 0 0 0 1 ... 0 0 0 1 0 ..];
    %
    x(word_indices) = 1;
end

author: worldveil

repo: coursera-ml


text:
%IM_SELECT_BLOB Fixed mapping selecting largest blob in binary images (DIP_Image)
%
%       B = IM_SELECT_BLOB(A)
%       B = A*IM_SELECT_BLOB
%
% Just the largest object in the image is returned.
%
% SEE ALSO (PRTools Guide)
% DATASETS, DATAFILES, DIP_IMAGE

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function b = im_select_blob(a)

  if nargin < 1 | isempty(a)
    b = prmapping(mfilename,'fixed');
    b = setname(b,'Select largest blob');
	elseif isa(a,'prdataset') % allows datafiles too
		isobjim(a);
    b = filtim(a,mfilename);
		b = setfeatsize(b,getfeatsize(a));
  elseif isa(a,'double') | isa(a,'dip_image') % here we have a single image
    if ~isa(a,'dip_image')
			a = dip_image(a,'bin');
		end;
		labim = label(a);
		sz = measure(labim,labim,{'size'});
		[cc,ind] = max(double(sz));
		I = double(measure(labim,labim,{'mean'}));
		b = a.*(labim==round(I(ind)));
%		c = measure(labim,labim,{'size','mean'});
%		c = double(c);
%		[cc,ind] = max(c(:,1));
%		b = a.*(labim==round(c(ind,2)));
	end

return

author: marianux

repo: ecg-kit

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


text:
function status = GE_convertADW(inDir,outStem,starttime,nvols)
%
% status = GE_convertADW(inDir,outStem,starttime,nvols)
%
%GE_convertADW
%
% Version 3.1
%
% Converts a series of GE slices acquired on the Advanced
% Development Workstation into Analyze format for SPM.
% inDir is the name of the first directory containing the
%      series, e.g. 003
% outStem is the stem used for naming the Analyze files
% starttime is the time point number to start
% nvols is the number of volumes (time points) to convert
%
% status = 1, error
% status = 0, all is well
%
% eg GE_convert('DATA/00023/003','converted',5,27)
% Will convert 27 time points starting at the 5th time point
% named converted_i0001.img, converted_i0002.img, ...,
% converted_i0027.img and their associated header and mat files.
%
% Assumes the data is stored as 003/I.001 through 003/I.999 and
% then 023/I.001 through 023/I.999, etc.
% Remember to skip the template images ie 1 time point for the
% first run.  You don't need to do this for subsequent runs.
%
% Modified to use spm write functions
% Modified to use spm_write_plane instead of spm_write_vol
%
% Souheil J. Inati
% Dartmouth College
% September 2001
% souheil.inati@dartmouth.edu
%

if (nargin < 4)
        error('Not enough input arguments.')
        return
end

% Create the name of the first file in inDir
firstfile = fullfile(inDir,'I.001');

% Read the Header from the first file
[su_hdr,ex_hdr,se_hdr,im_hdr,pix_hdr,im_offset] = GE_readHeader(firstfile);

% The image Dimensions
Nx =  im_hdr.imatrix_X; % X Voxels
Ny =  im_hdr.imatrix_Y; % Y Voxels
Nz =  im_hdr.slquant;   % Z Voxels
volSize = [Nx Ny Nz];

% Is the first image the first or the last?
if (se_hdr.end_loc - se_hdr.start_loc) > 0
  scandir = 1;
else
  scandir = -1;
end

% Compute the M matrix
M = GE_createSPMmat(im_hdr,scandir);

% Create the template SPM volume structure
V.fname = '';
V.dim = [volSize spm_type('int16')]; % short ints
V.mat = M;
V.pinfo = [1 0 0]';

% Loop over the volumes until there is no more looping to be done.
for i = 1:nvols
  passnum = starttime + i - 1;

  % Read in the volume
  [imageVol, lastfile] = GE_readVolume(inDir, passnum, volSize, ...
				       pix_hdr.img_depth, im_offset);

  % Create analyze file
  vol_str = sprintf('000%d',i);
  vol_str = vol_str(length(vol_str)-3:length(vol_str));
  outName = [outStem sprintf('_i%s.img',vol_str)];
  V.fname = outName;
  V = spm_create_image(V);

  % Write out the SPM Volume
  % Don't use spm_write_vol because that sets a scale factor.
  for i = 1:Nz
    V = spm_write_plane(V,squeeze(imageVol(:,:,i)),i);
  end

  % Update to the screen
  fprintf('Wrote %s\n',outName);

end

% Done
fprintf('\nConversion Finished. \n');

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function len = curveLength(varargin)
%CURVELENGTH return length of a curve (a list of points)
%
%   Compute the length of a curve given as a list of following points.
%
%   L = curveLength(X, Y);
%   L = curveLength(POINTS);
%   POINTS should be a [NxD] array, with N being the numbe of points and D
%   the dimension of the points.
%
%   PT = curveLength(..., TYPE);
%   Specifies if the last point is connected to the first one. TYPE can be
%   either 'closed' or 'open'.
%
%   TODO : specify norm (euclidian, taxi, ...).
%
%   Example:
%   Compute the perimeter of a circle with radius 1
%   curveLength(circleAsPolygon([0 0 1], 500), 'closed')
%   -> return 6.2831
%
%   See also:
%   polygons2d, curveCentroid
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 14/06/2004
%

%   HISTORY
%   22/05/2006 manage any dimension for points, closed and open curves,
%       and update doc accordingly.
%   30/06/2009 deprecate and replace by 'polylineLength'.

% deprecation warning
warning('geom2d:deprecated', ...
    '''curveLength'' is deprecated, use ''polylineLength'' instead');

% check whether the curve is closed
closed = false;
var = varargin{end};
if ischar(var)
    if strcmpi(var, 'closed')
        closed = true;
    end
    varargin = varargin(1:end-1);
end

% extract point coordinates
if length(varargin)==1
    points = varargin{1};
elseif length(varargin)==2
    points = [varargin{1} varargin{2}];
end

% compute lengths of each line segment
if closed
    len = sum(sqrt(sum(diff(points([1:end 1],:)).^2, 2)));
else
    len = sum(sqrt(sum(diff(points).^2, 2)));
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

visual_odo_initial( ...
    'ODO_IMG_TRANS_Y_RANGE', 51:150, ...
    'ODO_IMG_TRANS_X_RANGE', 181:300, ...
    'ODO_IMG_HEIGHT_V_Y_RANGE', 51:150, ...
    'ODO_IMG_HEIGHT_V_X_RANGE', 180:300, ...
    'ODO_IMG_YAW_ROT_Y_RANGE', 51:150, ...
    'ODO_IMG_YAW_ROT_X_RANGE', 180:300, ...
    'ODO_IMG_TRANS_RESIZE_RANGE', [100, 120], ...
    'ODO_IMG_YAW_ROT_RESIZE_RANGE', [100, 120], ...
    'ODO_IMG_HEIGHT_V_RESIZE_RANGE', [100, 120], ...
    'ODO_TRANS_V_SCALE', 1, ...
    'ODO_YAW_ROT_V_SCALE', 1, ...
    'ODO_HEIGHT_V_SCALE', 1, ...
    'MAX_TRANS_V_THRESHOLD', 0.02, ...
    'MAX_YAW_ROT_V_THRESHOLD', 10, ...
    'MAX_HEIGHT_V_THRESHOLD', 0.03, ...
    'ODO_SHIFT_MATCH_HORI', 36, ...
    'ODO_SHIFT_MATCH_VERT', 20, ...
    'FOV_HORI_DEGREE', 90, ...
    'FOV_VERT_DEGREE', 50, ...
    'KEY_POINT_SET', [1644, 1741], ...
    'ODO_STEP', 1);

% global avg;
visual_odo_main('C:\NeuroSLAM_Datasets\01_NeuroSLAM_Datasets\02_SynPanData', ...
    'C:\NeuroSLAM_Datasets\02_NeuroSLAM_Groudtruth\02_SynPanData_GT.txt');

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function [beta_swap]=betaswap(beta_gibbs,n,m,p,k)

% function [beta_swap]=bear.betaswap(beta_gibbs,n,m,p,k)
% reorganizes the matrix of gibbs sampler draws of beta, in order to make it easier to plot with matlab "subplot" function
% used to plot the empirical posterior distributions
% inputs:  - matrix 'beta_gibbs': record of the gibbs sampler draws for the beta vector
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'm': number of exogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - integer 'k': number of coefficients to estimate for each equation in the BVAR model (defined p 7 of technical guide)
% outputs: - matrix 'beta_swap': a reorganised gibbs sampler matrix

for ii=1:n
   for jj=1:n
      for kk=1:p
      beta_swap((ii-1)*k+(jj-1)*p+kk,:)=beta_gibbs((ii-1)*k+n*(kk-1)+jj,:);
      end
   end
end

for ii=1:n
beta_swap(ii*k-m+1:ii*k,:)=beta_gibbs(ii*k-m+1:ii*k,:);
end

author: european-central-bank

repo: BEAR-toolbox


text:
function [ state ] =movePTP_ConditionalTorque_ArcXZ_AC(t_Kuka,theta,c,VEL,joints_indices,max_torque,min_torque)
%% This function is used for moving the endeffector on an arc in the XZ plane, for the KUKA iiwa 7 R 800.
% this motion is interruptible when one or more of the joints-torques exceed the predefined limits

%% Syntax:
% [ state ] =movePTP_ConditionalTorque_ArcXZ_AC(t_Kuka,theta,c,VEL,joints_indices,max_torque,min_torque)

%% About:
% This function is used to move the end-effector on an arc in the XZ plane,

%% Arreguments:
% t_Kuka: is the TCP/IP connection
% theta: is the arc angle, in radians
% c: the XZ coordinates of the center of the circle, it is 1x2 vector.
% VEL : is a double, defines the motion velocity mm/sec.
% joints_indices: is a vector of the indices of the joints where the
% torques limits are to be imposed, the joints are indexed starting from one.
% max_torque: is a vector of the maximum torque limits for the joints
% specified in the (joints_indices) vector.
% min_torque: is a vecotr of the minimum torque limits for the joints
% specified in the (joints_indices) vector.

%% Return value:
% state: a number equals to one if the motion is completed sccessfully or a
% zero if the motion was interrupted due to external contact with the
% robot.
% if there is an error, the return value is minus one

% Copy right, Mohammad SAFEEA, 9th of April 2018

if((size(VEL,1)==1)&&(size(VEL,2)==1))
else
    disp('Error, velocity is a double and shall not be an array');
    state=-1;
    return;
end

k=[0;1;0];
pos=getEEFPos( t_Kuka );
c=colVec(c);
c1=[c(1);pos{2};c(2)];
[ state ] =movePTP_ConditionalTorque_Arc_AC(t_Kuka,theta,c1,k,VEL,joints_indices,max_torque,min_torque);

end

function [ y ] = colVec( v)
% Convert a vector to a column vector:
    if(size(v,2)==1)
        y=v;
    else
        y=v';
    end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function cellarray = remove_str_from_cell(cellarray,str)

    wh = cellfun(@isempty, strfind(cellarray, str));
    cellarray = cellarray(wh);

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function mrk= mrk_appendMarkers(mrk1, mrk2, Tmsec)
%MRK_APPENDMARKERS - Append (timewise) two marker structures
%
%Synopsis:
%  MRK= mrk_appendMarkers(MRK1, MRK2, TMSEC)
%  MRK= mrk_appendMarkers(MRK1, MRK2, CNT)
%
%Arguments:
%  MRK1, MRK2: marker structures (or empty [])
%  TMSEC:      duration of first data set in msec
%  CNT:        struct (like cnt) with fields T (duration of first data set
%              in samples) and fs (sampling rate).
%
%Returns:
%  mrk:        marker structure

misc_checkType(mrk1, 'STRUCT(time)');
misc_checkType(mrk2, 'STRUCT(time)');
misc_checkType(Tmsec, '!DOUBLE[1]|struct(T fs)');

if isempty(mrk1),
  mrk= mrk2;
  return;
end

if isstruct(Tmsec),
  cnt= Tmsec;
  Tmsec= sum(cnt.T)*1000/cnt.fs;
end

mrk2.time= mrk2.time + Tmsec;
if isempty(mrk1),
  mrk= mrk2;
  return;
end

if isfield(mrk1, 'event') && isfield(mrk1.event, 'blkno') && ...
    isfield(mrk2, 'event') && isfield(mrk2.event, 'blkno'),
  offset= max(mrk1.event.blkno) - min(mrk2.event.blkno) + 1;
  mrk2.event.blkno= mrk2.event.blkno + offset;
end

mrk= mrk_mergeMarkers(mrk1, mrk2);

author: bbci


text:
classdef ConstraintValues < matlab.mixin.SetGet
    %ConstraintValues Summary of this class goes here
    %   Detailed explanation goes here

    properties
        c = [];
        ceq = [];
        value = [];
        lb = [];
        ub = [];
        type = {};
        eventNum = [];
        cEventInds = [];
        ceqEventInds = [];
        consts AbstractConstraint = AbstractConstraint.empty(1,0);
        cCInds = [];
        cCeqInds = [];
        valueStateComps = [];
    end

    methods
        function obj = ConstraintValues()

        end

        function updateValues(obj, c, ceq, value, lb, ub, type, eventNum, cEventInds, ceqEventInds, consts, cCInds, cCeqInds, valueStateComps)
            obj.c = c;
            obj.ceq = ceq;
            obj.value = value;
            obj.lb = lb;
            obj.ub = ub;
            obj.type = type;
            obj.eventNum = eventNum;
            obj.cEventInds = cEventInds;
            obj.ceqEventInds = ceqEventInds;
            obj.consts = consts;
            obj.cCInds = cCInds;
            obj.cCeqInds = cCeqInds;
            obj.valueStateComps = valueStateComps;
        end

        function [c, ceq, value, lb, ub, type, eventNum, cEventInds, ceqEventInds, consts, cCInds, cCeqInds, valueStateComps] = getValues(obj)
            c = obj.c;
            ceq = obj.ceq;
            value = obj.value;
            lb = obj.lb;
            ub = obj.ub;
            type = obj.type;
            eventNum = obj.eventNum;
            cEventInds = obj.cEventInds;
            ceqEventInds = obj.ceqEventInds;
            consts = obj.consts;
            cCInds = obj.cCInds;
            cCeqInds = obj.cCeqInds;
            valueStateComps = obj.valueStateComps;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function sqi = ksqi(signal)
%kSQI Kurtosis SQI
%
% Returns the kurtosis of a signal
%
% Reference:
% Li, Q., Mark, R. G., Clifford, G. D., & Li. (2008). Robust heart rate
% estimation from multiple asynchronous noisy sources using signal quality
% indices and a Kalman filter. Physiol. Meas., 29(1), 15–32.
% http://doi.org/10.1088/0967-3334/29/1/002
%
% Input:
%   signal:         single channel (F)ECG [1xN double]
%
% Output:
%   sqi:            resulting kSQI for segment
%
% Fetal Extraction Toolbox, version 1.0, February 2014
% Released under the GNU General Public License
%
% Copyright (C) 2014 Fernando Andreotti
% Dresden University of Technology, Institute of Biomedical Engineering
% fernando.andreotti@mailbox.tu-dresden.de
%
% Last updated : 09-03-2014
%
%
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version.
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details.

sqi=kurtosis(signal);

end

author: fernandoandreotti


text:
function [ lambdaF, lambdaR, alphaF, alphaR ] = TireSlips(DeltaSteer_rad, VehicleStates, tyreradius_front_m, tyreradius_rear_m, l_front_m, l_rear_m)

% calculation of longitudinal and lateral tire slips

% minimum speed for slip calculations
vx_min = 0.5;

% get states
vx = VehicleStates(1);
vy = VehicleStates(2);
dPsi = VehicleStates(3);
omegaF = VehicleStates(4);
omegaR = VehicleStates(5);

% check if velocity is high enough that there are significant slips
if(vx > vx_min)
  % calculate side slip angles based on exact side slip formulas
  alphaF = DeltaSteer_rad - atan((vy + dPsi*l_front_m)/vx);
  alphaR = - atan((vy - dPsi*l_rear_m)/vx);
  lambdaF = (omegaF - vx/tyreradius_front_m)/(vx/tyreradius_front_m);
  lambdaR = (omegaR - vx/tyreradius_rear_m)/(vx/tyreradius_rear_m);
else
  alphaF = 0;
  alphaR = 0;
  lambdaF = 0;
  lambdaR = 0;
end

author: TUMFTM


text:
% File: logsumexp.m
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

function out = logsumexp(A)

% LOGSUMEXP
% Computes log( sum( exp( ) ) ) of each row in A in a way that avoids underflow.
% If A is an N x M matrix, then out is a N x 1 vector.

pi_max = max(A, [], 2);
out = pi_max + log(sum(exp(bsxfun(@minus, A, pi_max)), 2));

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function [panelContainer, bstPanel] = gui_show( bstPanel, contType, contName, contIcon, isModal, isAlwaysOnTop, isMaximized, winPos )
% GUI_SHOW: Display a BstPanel.
%
% USAGE: gui_show(bstPanel, contType, contName, contIcon=[], isModal=0, isAlwaysOnTop=0, isMaximized=0)

% INPUT:
%     - bstPanel : BstPanel object or panel function name (string)
%     - contType : {'JavaWindow', 'BrainstormTab', 'BrainstormPanel'} (Default : 'JavaWindow')
%     - contName : Container name (meaning of this argument depends on the contType value)
%                  => JavaWindow       : name of the Java frame (title of the java window)
%                  => BrainstormFigure : name of the panelContainer in which the panel should be added
%     - contIcon : Icon object (optional)
%     - isModal  : If 1, makes window modal
%                  => If the contType is 'BrainstormTab', the panel will not be added as a tab
%                     in a Brainstorm panel container, but will fill the whole panel container.
%     - isAlwaysOnTop : If 1, makes window always on top (only for type 'JavaWindow')
%     - isMaximized   : If 1, maximizes window to fit all the available figure space (only for type 'JavaWindow')
%     - winPos        : [x,y] Window position, relative with the main Braintorm frame (only for type 'JavaWindow')
%
% SEE ALSO gui_hide gui_show_dialog

% @=============================================================================
% 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, 2008-2017

import org.brainstorm.icon.*;
global GlobalData;

% Initialize returned variables
panelContainer = [];
% Headless mode: exit
if (GlobalData.Program.GuiLevel == -1)
    return
end

%% ===== PARSE INPUTS =====
if (nargin < 2) || isempty(contType)
    contType = 'JavaWindow';
end
if (nargin < 3) || isempty(contName) || ~ischar(contName)
    contName = 'New panel';
end
if (nargin < 4) || isempty(contIcon) || ~isa(contIcon, 'javax.swing.Icon')
    contIcon = [];
end
if (nargin < 5) || isempty(isModal)
    isModal = 0;
end
if (nargin < 6) || isempty(isAlwaysOnTop)
    isAlwaysOnTop = 0;
end
if (nargin < 7) || isempty(isMaximized)
    isMaximized = 0;
end
if (nargin < 8) || isempty(winPos)
    winPos = [];
end

%% ===== GET PANEL =====
if ischar(bstPanel)
    bstPanel = feval(bstPanel, 'CreatePanel');
    if isempty(bstPanel)
        return;
    end
elseif ~isa(bstPanel, 'BstPanel')
    error('First argument must be a ''BstPanel'' object.');
end
panelName    = get(bstPanel, 'name');
jPanelHandle = get(bstPanel, 'jHandle');
% Check that the panel is not already registered
if gui_brainstorm('isTabVisible', panelName)
    % Focus previous panel
    bstPanel = bst_get('Panel', panelName);
    jPanelHandle = get(bstPanel, 'jHandle');
    if ~isempty(jPanelHandle.getRootPane()) && ~isempty(jPanelHandle.getRootPane.getParent())
        jPanelHandle.getTopLevelAncestor().show();
    end
    % Return container
    panelContainer = get(bstPanel, 'container');
    return
end

%% ===== SWITCH BETWEEN CONTAINER TYPES =====
switch (contType)
    % =======================================================================
    % ==== JAVA WINDOW ======================================================
    % =======================================================================
    case 'JavaWindow'
        % Get parent frame
        jParent = bst_get('BstFrame');
        % If modal window : JDIALOG
        if (isModal)
            % Make it dependent on the the main Brainstorm JFrame, if it exists
            jWindow = java_create('javax.swing.JDialog', 'Ljava.awt.Frame;Ljava.lang.String;Z', jParent, contName, 1);
        % If not modal : JFRAME
        else
            jWindow = java_create('javax.swing.JFrame', 'Ljava.lang.String;', contName);
            jWindow.setIconImage(IconLoader.ICON_APP.getImage());
        end
        % Configure window
        jWindow.setDefaultCloseOperation(javax.swing.JFrame.DO_NOTHING_ON_CLOSE);
        java_setcb(jWindow, 'WindowClosingCallback', @ClosePanelCallback);
        % Add the BST Panel
        jWindow.getContentPane.add(jPanelHandle);
        % Always on top
        if isAlwaysOnTop
            jWindow.setAlwaysOnTop(1);
        end
        % Maximized over the figure area (ONLY WHEN TILED LAYOUT)
        if isMaximized && strcmpi(bst_get('Layout', 'WindowManager'), 'TileWindows')
            % Get current layout
            [jBstArea, FigArea, nbScreens, jFigArea, jInsets] = gui_layout('GetScreenBrainstormAreas');
            % Place this figure on the screen to take all the possible figure space
            jWindow.setSize(FigArea(3), FigArea(4));
            jWindow.setLocation(FigArea(1)-1, jBstArea.getY() - jInsets.bottom);
        % Normal display
        else
            % Set window size and location
            jWindow.pack();
            jWindow.setLocationRelativeTo(jParent);
            % Set figure position
            if (length(winPos) >= 2)
                if ~isempty(jParent)
                    winPos(1) = winPos(1) + jParent.getLocation().getX();
                    winPos(2) = winPos(2) + jParent.getLocation().getY();
                end
                jWindow.setLocation(winPos(1),winPos(2));
            end
        end
        % Show window: KEEP IT WITH AWTINVOKE
        awtinvoke(jWindow, 'setVisible(Z)', 1);
        % Returned variable
        panelContainer.type   = 'JavaWindow';
        panelContainer.handle = {jWindow};

    % =======================================================================
    % ==== BRAINSTORM TAB ===================================================
    % =======================================================================
    case 'BrainstormTab'
        % Add a Swing components tree to the main brainstorm window.
        % The component can be added in one of the registered panel containers.
        % The panel can be added either as a new tab in the target area, or
        % as the whole target (in this case, the other target tabbed panels are hidden)

        % Get the 'contName' tabbed panel
        jTabContainer = bst_get('PanelContainer', contName);
        % If container does not exist, return with a warning
        if (isempty(jTabContainer))
            warning('Brainstorm:InvalidContainer', 'Destination "%s" is not a valid Swing tab panel', contName);
            return;
        end
        % If this panel is supposed to tabbed => just add it as a new tab
        if (~isModal)
            % Attach the panel to the panel container
            if strcmpi(contName, 'tools')
                iTab = jTabContainer.getTabCount() - 1;
            else
                iTab = jTabContainer.getTabCount();
            end
            jTabContainer.insertTab(panelName, contIcon, jPanelHandle, '', iTab);
            % Make sure that the tabbed container is visible
            jTabContainer.setVisible(1);
        % Else : 1) add it to the parent of the jTabContainer
        else
            % Add panel in the container's parent
            jTabContainer.getParent().add(jPanelHandle);
            % Hide the tab container
            jTabContainer.setVisible(0);
        end

        % Returned variable
        panelContainer.type   = 'BrainstormTab';
        panelContainer.handle = {contName};

    % =======================================================================
    % ==== BRAINSTORM PANEL ===================================================
    % =======================================================================
    case 'BrainstormPanel'
        % Get the 'contName' tabbed panel
        jContainer = bst_get('PanelContainer', contName);
        % If container does not exist, return with a warning
        if (isempty(jContainer))
            warning('Brainstorm:InvalidContainer', 'Destination "%s" is not a valid Swing panel', contName);
            return;
        end
        % Add panel in the container
        jContainer.add(jPanelHandle);
        % Returned variable
        panelContainer.type   = 'BrainstormPanel';
        panelContainer.handle = {contName};

    otherwise
        error('Invalid container type.');
end

%% ===== REFERENCE NEW PANEL =====
% Update bstPanel structure
bstPanel = set(bstPanel, 'container', {panelContainer.type, panelContainer.handle{:}});
% If GUI context is defined
if isempty(GlobalData.Program.GUI.panels)
    GlobalData.Program.GUI.panels = bstPanel;
else
    % Add the new panel to the panel list
    GlobalData.Program.GUI.panels(end+1) = bstPanel;
end

%% =================================================================================
%  === CALLBACKS  ==================================================================
%  =================================================================================
    function ClosePanelCallback(varargin)
        try
            gui_hide(bstPanel);
        catch
            if strcmpi(contType, 'JavaWindow')
                warning('Brainstorm:PanelNotRegistered', 'Panel ''%s'' was not registered in Brainstorm GUI context, hide forced.', panelName);
                jWindow.dispose();
            end
        end
    end

end

author: brainstorm-tools


text:
function test_subplot_2D_3D
%TEST_SUBPLOT_2D_3D Test mixed subplot 2D/3D with DRAGZOOM

figure;
hax1 = subplot(2,1,1);
x = -pi*2:0.1:pi*2;
y = sin(x);
plot(x,y);

hax2 = subplot(2,1,2);
k = 5;
n = 2^k-1;
[x,y,z] = sphere(n);
surf(x,y,z);

dragzoom();

author: Sable


text:
function jHelp = bst_help(helpfile, isModal, figTitle, figSize, jHelp)
% BST_HELP: Display a help file in HTML format.
%
% USAGE:  jHelp = bst_help(helpfile, isModal=1, figTitle='Help', figSize=[510,590], jHelp=[new])
%         jHelp = bst_help(strHtml, ...)

% @=============================================================================
% 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, 2008-2012

% Initializations
import java.awt.*;
import javax.swing.*;

% ===== PARSE INPUTS =====
if (nargin < 5) || isempty(jHelp)
    jHelp = [];
end
if (nargin < 4) || isempty(figSize)
    figSize = [510,590];
end
if (nargin < 3) || isempty(figTitle)
    figTitle = 'Help';
end
if (nargin < 2) || isempty(isModal)
    isModal = 1;
end
if (helpfile(1) == '<')
    % Get HTML help string
    strHtml = helpfile;
    helpfile = [];
else
    % Figure title
    figTitle = [figTitle ': ' helpfile];
    % Get documents path
    helppath = bst_fullfile(bst_get('BrainstormDocDir'), 'help');
    % Get full filename
    helpfile = bst_fullfile(helppath, helpfile);
end

% ===== DIALOG INITIALIZATION =====
jBstFrame = bst_get('BstFrame');
% Create a new window
if isempty(jHelp)
    jFrame = java_create('javax.swing.JDialog', 'Ljava.awt.Frame;Ljava.lang.String;Z', jBstFrame, figTitle, isModal);
    jFrame.setDefaultCloseOperation(JFrame.DISPOSE_ON_CLOSE);
    jFrame.setPreferredSize(Dimension(figSize(1), figSize(2)));
    jPanelMain = jFrame.getContentPane();
% Re-use an existing window
else
    jFrame = jHelp.jFrame;
    jPanelMain = jFrame.getContentPane();
    jPanelMain.removeAll();
end

% === HELP TEXT AREA ===
% Create HTML viewer component
jTextHtml = JEditorPane();
jTextHtml.setEditable(false);
% Add an html editor kit
kit = javax.swing.text.html.HTMLEditorKit();
jTextHtml.setEditorKit(kit);
% Set content
if ~isempty(helpfile)
    jTextHtml.setPage(java.net.URL(['file:///' helpfile]));
else
    doc = kit.createDefaultDocument();
    jTextHtml.setDocument(doc);
    %jTextHtml.setContentType('text/html');
    jTextHtml.setText(strHtml);
end
% Callback when clicking on hyperlinks
java_setcb(jTextHtml, 'HyperlinkUpdateCallback', @HyperTextUpdate);
% Enclose in scroll panel
jScrollText = java_create('javax.swing.JScrollPane', 'Ljava.awt.Component;', jTextHtml);
jPanelMain.add(jScrollText, BorderLayout.CENTER);

% === FOOTER ===
% Add a "Close" button
jButtonClose = JButton('Close');
java_setcb(jButtonClose, 'ActionPerformedCallback', @CloseDialog);
jPanelMain.add(jButtonClose, BorderLayout.SOUTH);

% Display figure
jFrame.pack();
jFrame.setLocationRelativeTo(jFrame.getParent());
jFrame.setVisible(1);

% Return handle
jHelp.close  = @CloseDialog;
jHelp.jFrame = jFrame;

%% =================================================================================================
%  ====== CALLBACKS ================================================================================
%  =================================================================================================
%% ===== CLOSE DIALOG =====
    function CloseDialog(varargin)
        drawnow
        % If dialog was already closed
        if ~jFrame.isVisible()
            return
        end
        % Release mutex
        bst_mutex('release', 'Help');
        % Close dialog
        jFrame.dispose();
    end

%% ===== HYPERTEXT =====
    function HyperTextUpdate(h, ev)
        % If user clicked link
        if (ev.getEventType() == ev.getEventType().ACTIVATED)
            % Remove the 'Always on top' attribute
            jFrame.setAlwaysOnTop(0);
            % Display web page with Matlab browser
            web(char(ev.getURL()), '-browser');
        end
    end
end


text:
function [out,cl] = imageCluster_block(varargin)
% :Usage:
% ::
%
%    [out,cl] = imageCluster_block(input arguments)
%
% Images a cluster isosurface on an existing 3D head or brain plot
%
% works with tsu (Talairach Space Utility) main window as current figure
% old: clusters = getappdata(gcf, 'clusters');
%
% Inputs (in any order): keyword followed by input argument
%
%   **'cluster':**
%        followed by cluster to image, from SPM or TSU.
%
%   **'getclusters':**
%        no other args necessary - starts gui for cluster selection
%        function returns all clusters.  select with clusters(i)
%
%   **'getfigclusters':**
%        get clusters from TSU main figure.  must be current figure.
%
%   **'figure':**
%        create a new figure to image color on
%
%   **'color':**
%        followed by color value - either text or vector
%
%   **'alpha':**
%        followed by transparency value for cluster surface, 0-1
%        1 is opaque, 0 is completely transparent
%
% :Outputs:
%
%   **out:**
%        Patch handle
%
%   **cl:**
%        cluster struct
%
% ..
%    Works in Matlab 5.3, but with no transparency.
%    By Tor Wager, 10/3/2001, last edit 10/18/01
% ..

% ..
%    Set up arguments and default values
% ..
clusters = [];
mycolor = 'y';
myalpha = 1;

for i = 1:nargin
    if isstr(varargin{i})
        switch varargin{i}
        case 'cluster', cl = varargin{i+1};
        case 'getclusters', clusters = tor_ihb_getClusters;
        case 'getfigclusters', clusters = getappdata(gcf, 'clusters');
        case 'figure', h3dfig = figure; set(h3dfig,'Tag','myFig');
        case 'color',mycolor = varargin{i+1};
        case 'alpha',myalpha = varargin{i+1};
        end
    end
end

% If you chose to get clusters, no imaging done - just returns all clusters.
if ~isempty(clusters), out = clusters;, return, end

% do not do this [X,Y,Z] = meshgrid(x,y,z);

%------------------------------------------------------------------------------
% Prepare volume data and create isosurface patch object
%------------------------------------------------------------------------------
    try
	    V = cl.vTal;
        %V = smooth3(cl.vTal,'gaussian',[5 5 5]);
        V = smooth3(cl.vTal);
    catch
	    % no vTal field, make it here.
	    disp(['No vTal field for cluster, attempting to create it.'])
	    cl = get_cluster_volume(cl);
	    V = cl.vTal;
        %V = smooth3(cl.vTal,'gaussian',[5 5 5]);
        V = smooth3(cl.vTal);
    end

	if isfield(cl,'hThreshold'), height_threshold = cl.hThreshold;
	elseif isfield(cl,'threshold'), height_threshold = cl.threshold;
	else height_threshold = 0;,warning('No height threshold found in cluster.')
	end

    %height_threshold = min(cl.Z);
    % to image the whole cluster.  Smoothing will increase the cluster size if you choose this!
    height_threshold = 0;   % set this to 0 because all above 0 should be imaged.

    % use with smoothing, to smooth out the clusters and keep them a reasonable size
    %height_threshold = .5 * mean(mean(mean(V)));

 	set(0,'CurrentFigure', gcf)
	%set(gcf, 'CurrentAxes', gca);

	FV = isosurface(cl.xTal, cl.yTal, cl.zTal, V, height_threshold);
	hPatch = patch(FV);
    IC = isocaps(cl.xTal, cl.yTal, cl.zTal, V, height_threshold);
    iPatch = patch(IC);
	isonormals(cl.xTal, cl.yTal, cl.zTal, cl.vTal, hPatch);
	try
		set(hPatch, 'Tag', 'ihb_surfaceCluster', 'FaceColor', mycolor, 'EdgeColor', 'none', 'FaceAlpha', myalpha);
        set(iPatch, 'Tag', 'ihb_surfaceCluster', 'FaceColor', mycolor, 'EdgeColor', 'none', 'FaceAlpha', myalpha);
	catch
		set(hPatch, 'Tag', 'ihb_surfaceCluster', 'FaceColor', mycolor, 'EdgeColor', 'none');
        set(iPatch, 'Tag', 'ihb_surfaceCluster', 'FaceColor', mycolor, 'EdgeColor', 'none');
	end

    %delete(hMsg);
    set(gcf, 'Visible', 'on');
    out = hPatch;

return


text:
% pdfbclassify_texture.m
% written by: Duncan Po
% Date: December 3, 2002
% perform texture classification based on contourlets
%
% Usage:    kld = pdfbclassify_texture(qimage, qformat, tdb, tdir, mdir, firsttime)
% Inputs:   qimage      - query texture image file name
%           qformat     - format of the query texture image file
%           tdb         - texture database: vector (using cell structure)
%                           of the names of all the texture images
%                           e.g. {'texture1', 'texture2', 'texture3'}
%           tdir        - texture directory: full path of the directory
%                           that contains the texture images
%           mdir        - model directory: full path of the directory that
%                           either contains the texture models, or that is
%                           for saving of the texture models
%           firsttime   - set to 1 if this is the first time that this
%                           algorithm is run and no model is available
%                           for the database images, set to 0 otherwise
% Output:   kld         - the Kublick Liebler distance between the trees

function kld = pdfbclassify_texture(qimage, qformat, tdb, tdir, mdir, firsttime)

N = 8;
NN = 16;
if tdir(end) ~= '/'
    tdir = strcat(tdir, '/');
end;
if mdir(end) ~= '/'
    mdir = strcat(mdir, '/');
end;

mdir = strcat(mdir, '/');

[kld, kld2] = pdfbtestall_imagekld(qimage, qformat, tdb, tdir, mdir, N, firsttime);

[skld2, ind2] = sort(kld2);

ind2 = ind2(1:NN-1);
original = imread(qimage, qformat);
figure;
subplot(4,4,1);
imshow(original);

for ploti = 2:NN
    dbimagenum = floor((ind2(ploti-1)-1)/N)+1;
    dbimage = tdb{dbimagenum};
    dbimagenum = mod(ind2(ploti-1)-1, N)+1;
    file = sprintf('%s%s%02d', tdir, dbimage, dbimagenum);
    qresult = imread(file, 'bmp');
    subplot(4,4,ploti);
    imshow(qresult);
end;


text:
classdef TestDrawDetectedMarkers
    %TestDrawDetectedMarkers

    methods (Static)
        function test_1
            [img, d] = get_image_markers();
            [corners, ids, rejected] = cv.detectMarkers(img, d.dict);
            out = cv.drawDetectedMarkers(img, corners, 'IDs',ids);
            if ~isempty(rejected)
                out = cv.drawDetectedMarkers(out, rejected, ...
                    'BorderColor',[255 255 0]);
            end
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_error_argnum
            try
                cv.drawDetectedMarkers();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

function [img, d] = get_image_markers()
    % 4 markers in a diamond
    d = struct();
    d.dict = '6x6_50';
    d.ids = randi([0 49], [1 4]);
    d.slen = 80;
    d.mlen = 50;
    img = cv.drawCharucoDiamond(d.dict, d.ids, d.slen, d.mlen);
    img = repmat(img, [1 1 3]);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [Xmat,pobjround,info] = local_search_quasar(Zmat,C,rrPar,opt,roundonly)
%% Round and refine for QUASAR SDP
%% Input:
%% Zmat: SDP primal iterate
%% C: cost matrix of the SDP
%% rrPar: a structure where you can pass any data necessary for local search
%% opt: indices of eigenvectors for rounding
%% roundonly: if true, then only round a solution without NLP
%% Output:
%% Xmat: a rank-one SDP iterate
%% pobjround: primal SDP cost attained by Xmat

if nargin < 5
    roundonly = false;
end
if nargin < 4
    % default round the first two eigenvectors
    opt = [1,2];
end

if roundonly
    [q,theta] = quasar_round(Zmat,[1]);
    xtld      = kron([1;theta],q);
    pobjround = xtld' * C{1} * xtld;
    Xmat      = {xtld * xtld'};
else
    [q,theta] = quasar_round(Zmat,opt);
    Zmatround = {};
    pobjround = zeros(length(opt),1);
    for i = 1:length(opt)
        [Ztmp,~,~,pobjtmp] = quasar_manopt(C,q(:,i),theta(:,i));
        Zmatround{end+1} = Ztmp;
        pobjround(i) = pobjtmp;
    end
    pobjs            = pobjround;
    [pobjround,idx]  = min(pobjround);
    if pobjround == inf
        fprintf('        NLP fails to find a good solution, return a rounded solution only.\n');
        qopt         = q(:,1); thetaopt = theta(:,1);
        xopt         = kron([1;thetaopt],qopt);
        Xmat         = {xopt * xopt'};
        nlpsuccess   = false;
    else
        Xmat         = Zmatround{idx};
        nlpsuccess   = true;
    end
end

if nargout > 2
    info.minidx     = idx;
    info.nlpsuccess = nlpsuccess;
    info.pobjs      = pobjs;
    info.diffpobj   = pobjs(1) - pobjround;
end
end

function [q,theta] = quasar_round(X,opt)
if nargin < 2
    opt = 1;
end
[V,D] = eig(X{1});
[~,I] = sort(diag(D),'descend');
V     = V(:,I);
%% take the opt-th eigenvector
x     = V(:,opt);
q     = [];
theta = [];
for k = 1:size(x,2)
    xk    = x(:,k);
    %% round that eigenvector
    n     = length(xk);
    N     = n/4 - 1;
    qk    = xk(blkIndices(1,4));
    qk    = qk/norm(qk);
    thetak= zeros(N,1);
    for i = 1:N
        inprod = qk' * xk(blkIndices(i+1,4));
        if inprod > 0
            thetak(i) = 1;
        else
            thetak(i) = -1;
        end
    end
    q     = [q,qk];
    theta = [theta,thetak];
end
end

function [Xmat,qnew,thetanew,fopt] = quasar_manopt(C,q,theta)
N           = length(theta);
s           = spherefactory(4,1);
o           = obliquefactory(1,N);
elements.A  = s;
elements.B  = o;
manifold    = productmanifold(elements);

problem.M   = manifold;

warning('off', 'manopt:getHessian:approx')
% Define the problem cost function and its Euclidean gradient.
problem.cost  = @(x) quasar_cost(x,C);
problem.egrad = @(x) quasar_egrad(x,C);

% Numerically check gradient consistency (optional).
% checkgradient(problem);
% Solve.
x0.A                 = q;
x0.B                 = theta';
options.verbosity    = 0;
options.tolgradnorm  = 1e-6;
[xopt, fopt, output, options] = trustregions(problem,x0,options);
qopt = xopt.A;
thetaopt = xopt.B';
constraintviolationq          = abs(norm(qopt) - 1);
constraintviolationtheta      = max(abs(thetaopt.^2 - 1));
firstorderopt                 = output(end).gradnorm;
fprintf('        MANOPT: itr: %3d, constraint violation: %3.2e, %3.2e, gradnorm: %3.2e, cost: %3.8e.\n',...
    length(output),constraintviolationq,constraintviolationtheta,firstorderopt,fopt);
if max(constraintviolationq,constraintviolationtheta) < 1e-8 ...
        && firstorderopt < 1e-6
    % Do nothing
else
    fopt    = inf;
end
qnew = qopt/norm(qopt);
thetanew = sign(thetaopt);
vnew     = kron([1;thetanew],qnew);
Xmat     = {vnew * vnew'};
end

function f = quasar_cost(x,C)
q       = x.A;
theta   = x.B;
v       = kron([1;theta(:)],q);
f       = v'*C{1}*v;
end

function g = quasar_egrad(x,C)
q       = x.A;
theta   = x.B;
N       = length(theta);
v       = kron([1;theta(:)],q);

gv      = 2 * v' * C{1}; % 1 x n

vdq     = kron([1;theta(:)],speye(4)); % n x 4
gq      = (gv * vdq)'; % 4 x 1

vdtheta = [sparse(4,N);kron(speye(N),q)]; % n x N
gtheta  = (gv * vdtheta); % 1 x N

g.A     = gq;
g.B     = gtheta;
end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
%% 3D Orientation Visualizations
%

%% Euler angle space
%
% By default the function <orientation.plot.html plot>
% plots orientations in the three dimensional Bunge Euler angle space

cs = crystalSymmetry('cubic')

ori = orientation.rand(100,cs);
plot(ori)

%%
% Note that the orientations are automatically projected into the
% fundamental region. In the case of cubic symmetry this means that the
% Euler angles $\Phi$ and $\phi_2$ are restricted to 90 degrees. If the
% orientations should be plotted at their specified Euler angles the option
% |'ignoreFundamentalRegion'| has to be used.

plot(ori,'ignoreFundamentalRegion')

%% Axis angle space
%
% Alternatively, orientations can be plotted in the three dimensional axis
% angle space.

plot(ori,'AxisAngle','markerEdgeColor',[0 0 0.8],'markerSize',8)

%%
% The orientations are automatically projected into its fundamental region.
% Again, this can be switched off with the option
% |'ignoreFundamentalRegion'|.

plot(ori,'axisAngle','ignoreFundamentalRegion','markerEdgeColor',[0 0 0.8],'markerSize',8)

% visualize the fundamental region
hold on
oR = fundamentalRegion(ori.CS,ori.SS)
plot(oR,'color',[1 0.5 0.5]),
hold off

author: mtex-toolbox

repo: mtex

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