text:
function [mu,isig] = spm_affine_priors(typ)
% Distribution of the priors used in affine registration
%
% The parameters for this distribution were derived empirically from 227
% scans, that were matched to the ICBM space.
%__________________________________________________________________________
% Copyright (C) 2003-2018 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_affine_priors.m 7377 2018-07-23 13:56:24Z guillaume $

% % Values can be derived by...
% sn = spm_select(Inf,'.*seg_inv_sn.mat$');
% X  = zeros(size(sn,1),12);
% for i=1:size(sn,1),
%     p  = load(deblank(sn(i,:)));
%     M  = p.VF(1).mat*p.Affine/p.VG(1).mat;
%     J  = M(1:3,1:3);
%     V  = sqrtm(J*J');
%     R  = V\J;
%     lV      =  logm(V);
%     lR      = -logm(R);
%     P       = zeros(12,1);
%     P(1:3)  = M(1:3,4);
%     P(4:6)  = lR([2 3 6]);
%     P(7:12) = lV([1 2 3 5 6 9]);
%     X(i,:)  = P';
% end
% mu   = mean(X(:,7:12));
% XR   = X(:,7:12) - repmat(mu,[size(X,1),1]);
% isig = inv(XR'*XR/(size(X,1)-1))

switch deblank(lower(typ))

case 'mni' % For registering with MNI templates...
    mu   = [0.0667 0.0039 0.0008 0.0333 0.0071 0.1071]';
    isig = 1e4 * [
        0.0902   -0.0345   -0.0106   -0.0025   -0.0005   -0.0163
       -0.0345    0.7901    0.3883    0.0041   -0.0103   -0.0116
       -0.0106    0.3883    2.2599    0.0113    0.0396   -0.0060
       -0.0025    0.0041    0.0113    0.0925    0.0471   -0.0440
       -0.0005   -0.0103    0.0396    0.0471    0.2964   -0.0062
       -0.0163   -0.0116   -0.0060   -0.0440   -0.0062    0.1144];

case 'imni' % For registering with MNI templates...
    mu   = -[0.0667 0.0039 0.0008 0.0333 0.0071 0.1071]';
    isig = 1e4 * [
        0.0902   -0.0345   -0.0106   -0.0025   -0.0005   -0.0163
       -0.0345    0.7901    0.3883    0.0041   -0.0103   -0.0116
       -0.0106    0.3883    2.2599    0.0113    0.0396   -0.0060
       -0.0025    0.0041    0.0113    0.0925    0.0471   -0.0440
       -0.0005   -0.0103    0.0396    0.0471    0.2964   -0.0062
       -0.0163   -0.0116   -0.0060   -0.0440   -0.0062    0.1144];

case 'rigid' % Constrained to be almost rigid...
    mu   = zeros(6,1);
    isig = eye(6)*1e8; % spm_affreg used 1e9

case 'subj' % For inter-subject registration...
    mu   = zeros(6,1);
    isig = 1e3 * [
        0.8876    0.0784    0.0784   -0.1749    0.0784   -0.1749
        0.0784    5.3894    0.2655    0.0784    0.2655    0.0784
        0.0784    0.2655    5.3894    0.0784    0.2655    0.0784
       -0.1749    0.0784    0.0784    0.8876    0.0784   -0.1749
        0.0784    0.2655    0.2655    0.0784    5.3894    0.0784
       -0.1749    0.0784    0.0784   -0.1749    0.0784    0.8876];

case 'eastern' % For East Asian brains to MNI...
    mu   = [0.0719   -0.0040   -0.0032    0.1416    0.0601    0.2578]';
    isig = 1e4 * [
        0.0757    0.0220   -0.0224   -0.0049    0.0304   -0.0327
        0.0220    0.3125   -0.1555    0.0280   -0.0012   -0.0284
       -0.0224   -0.1555    1.9727    0.0196   -0.0019    0.0122
       -0.0049    0.0280    0.0196    0.0576   -0.0282   -0.0200
        0.0304   -0.0012   -0.0019   -0.0282    0.2128   -0.0275
       -0.0327   -0.0284    0.0122   -0.0200   -0.0275    0.0511];

case 'none' % No regularisation...
    mu   = zeros(6,1);
    isig = zeros(6);

otherwise
    error('"%s" not recognised as type of regularisation.',typ);
end

author: spm


text:
function y = coshminus1(x)
%COSHMINUS1 Hyperbolic cosine minus one.
%
%   COSHMINUS1(X) is COSH(X)-1 calculated in a way that is numerically better
%   when X is close zero.
%
%   This illustrates the difference
%
%      x = 5e-8*(-1:0.01:1);
%      plot(x, cosh(x)-1, 'b-', x, coshminus1(x), 'r-');
%
%   See also ACOSHPLUS1.

%   Author:      Peter J. Acklam
%   Time-stamp:  2003-10-13 14:57:06 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   % check number of input arguments
   error(nargchk(1, 1, nargin));

   y = 2 * sinh(x / 2).^2;

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function seed = latin_center_test01 ( seed )

%*****************************************************************************80
%
%% LATIN_CENTER_TEST01 tests LATIN_CENTER.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 April 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer SEED, an initial seed for the random number generator.
%
%    Output, integer SEED, the updated random number seed.
%
  dim_num = 2;
  point_num = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  LATIN_CENTER chooses a Latin cell arrangement,\n' );
  fprintf ( 1, '  and returns the centers of those cells.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension = %d\n', dim_num );
  fprintf ( 1, '  Number of points =  %d\n', point_num );
  fprintf ( 1, '  Using seed = %d\n', seed );

  [ x, seed ] = latin_center ( dim_num, point_num, seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The Latin center points:\n' );
  fprintf ( 1, '\n' );

  for j = 1: point_num
    for i = 1: dim_num
      fprintf ( 1, '%10f  ', x(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  JOINT_REFERENCES
%  Obtain q when the robot makes a linear trajectory of the end effector
%  in cartesian space
%
%  See also INVERSEKINEMATIC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function q=joint_references(time)

global robot

time;
if time == 0
    robot.q=[0 0 0 0 0 0]';
end

simulation_time = 0.5;%s

%this is important, as this will be the initial joint
%values in the first integrator
%initial position
T1 = directkinematic(robot, [0 0 0 0 0.1 0]);
%final position
T2 = directkinematic(robot, [0.3 0.3 0.3 0.1 0.1 0.1]);

Q1 = T2quaternion(T1);
Q2 = T2quaternion(T2);
%interpolate between the two orientations
[Qm] = slerp(Q1, Q2, time/simulation_time, 0.01);

initial_point=T1(1:3,4);
end_point=T2(1:3,4);

% %NOA matrix end point
% T2=[1 0 0 0.5;
%     0 1 0 0.4;
%     0 0 1 0.3;
%     0 0 0  1];

%speed of the movement
speed = norm(initial_point-end_point)/simulation_time;%m/s

v=(end_point-initial_point);
v=v/norm(v); %unit vector in the direction of the line

point = initial_point + speed*time*v;

%T1(1:3,4)=point(1:3);

Ttotal = quaternion2T(Qm, point(1:3));
qinv = inversekinematic(robot, Ttotal);
q = select_closest_joint_coordinates(qinv, robot.q);

%update current robot joints
robot.q=q;

repo: ARTE


text:
function [indTheta, Xi, numterms] = ADMinitvary(nT,lambda,MaxIter,tol, plotflag)

% normalize the collumns of the null space of Theta to use in the
% initial conditoin for ADM search
for ii = 1:size(nT,2)
    nTn(:,ii) = nT(:,ii)/mean(nT(:,ii));
end

% run ADM algorithm on each row of nTn
for jj = 1:size(nTn,1)
    jj;
    qinit = nTn(jj,:)'; % intial conditions
    q(:,jj) = ADM(nT, qinit, lambda, MaxIter, tol); % algrorithm for
    % finding coefficients resutling in  sparsest vector in null space
    out(:,jj) = nT*q(:,jj); % compose sparsets vectors
    nzeros(jj)= length(find(abs(out(:,jj))<lambda)); % chech how many zeros each
    % of the found sparse vectors have
end

indsparse = find(nzeros==max(nzeros)); % find the vector with the largest number of zeros
indTheta = find(abs(out(:,indsparse(1)))>=lambda); % save the indices of non zero coefficients
% Xi = out(indTheta, indsparse(1)); % get those coefficients
Xi = out(:, indsparse(1)); % get sparsest vector
smallinds =(abs(out(:,indsparse(1)))<lambda);
Xi(smallinds)=0; % set thresholded coefficients to zero

% check that the solution found by ADM is unique.
if length(indsparse)>1
    Xidiff= (out(indTheta, indsparse(1))-out(indTheta,indsparse(2)));
    if Xidiff>tol
        warning('ADM has discovered two different sparsest vectors')
    end
end
% calculate how many terms are in the sparsest vector
numterms = length(indTheta);

if plotflag == 2
    figure(121)
    semilogy(size(nT,1)-nzeros, 'o');
    drawnow
end

author: dynamicslab

repo: SINDy-PI


text:
function title = p17_title ( option )

%*****************************************************************************80
%
%% P17_TITLE sets the title for problem 17.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, string TITLE, the title of the problem.
%
  if ( option == 1 )
    title = 'Bremermann Propane Combustion System, fixed pressure.';
  else
    title = 'Bremermann Propane Combustion System, fixed concentration.';
  end

  return
end


text:
function v = normc(v)
% normalize along first axis (works for GPU arrays)
v = v./repmat(sum(v.^2, 1), size(v,1),1).^.5;

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function [gx] = g_LinDecomp(x,P,u,in)
X = NaN(in.size.X(1),in.size.X(2));
Y = NaN(in.size.Y(1),in.size.Y(2));
for i=1:in.n
    X(:,i) = P(in.ind(i).X);
    Y(i,:) = P(in.ind(i).Y)';
end
gx = VBA_vec(X*Y) + P(in.ind0);

author: MBB-team

repo: VBA-toolbox


text:
classdef TestBlendLinear
    %TestBlendLinear

    methods (Static)
        function test_1
            src1 = imread(fullfile(mexopencv.root(),'test','left01.jpg'));
            src2 = imread(fullfile(mexopencv.root(),'test','right01.jpg'));
            assert(isequal(size(src1), size(src2)));

            [rows,cols,~] = size(src1);
            weights1 = 0.5*ones(rows, cols);
            weights2 = 0.5*ones(rows, cols);
            dst = cv.blendLinear(src1, src2, weights1, weights2);
            validateattributes(dst, {class(src1)}, {'size',size(src1)});
        end

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

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function rcnn_make_window_file(imdb, out_dir, mode)
% rcnn_make_window_file(imdb, out_dir)
%   Makes a window file that can be used by the caffe WindowDataLayer
%   for finetuning.
%
%   The window file format contains repeated blocks of:
%
%     # image_index
%     img_path
%     channels
%     height
%     width
%     num_windows
%     class_index overlap x1 y1 x2 y2
%     <... num_windows-1 more windows follow ...>

% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------
if nargin<3
    mode = 'cnn';
end

roidb = imdb.roidb_func(imdb);

window_file = sprintf('%s/window_file_%s.txt', ...
    out_dir, imdb.name);
fid = fopen(window_file, 'wt');

channels = 3; % three channel images

for i = 1:length(imdb.image_ids)
  tic_toc_print('make window file: %d/%d\n', i, length(imdb.image_ids));
  img_path = imdb.image_at(i);

  roi = roidb.rois(i);

  switch mode
    case 'cnn'
        % do nothing
    case 'ctx'
        roi.boxes = roi.ctx;
    otherwise
        disp('Unspecified mode!');
        keyboard;
  end

  num_boxes = size(roi.boxes, 1);
  fprintf(fid, '# %d\n', i-1);
  fprintf(fid, '%s\n', img_path);
  fprintf(fid, '%d\n%d\n%d\n', ...
      channels, ...
      imdb.sizes(i, 1), ...
      imdb.sizes(i, 2));
  fprintf(fid, '%d\n', num_boxes);
  for j = 1:num_boxes
    [ov, label] = max(roi.overlap(j,:));
    % zero overlap => label = 0 (background)
    if ov < 1e-5
      label = 0;
      ov = 0;
    end
    bbox = roi.boxes(j,:)-1;
    fprintf(fid, '%d %.3f %d %d %d %d\n', ...
        label, ov, bbox(1), bbox(2), bbox(3), bbox(4));
  end
end

fclose(fid);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% Author:  Javier Lopez-Calderon
% 2009
%
function m = modeone(A)
if iscell(A)
        while any(cellfun('isclass',A,'cell'))
                A = cat(2,A{:});
        end

        n = length(A);
        h = cellfun(@isnumeric,A);
        g = nnz(h);

        if g>0 && g~=n % mixed
                error('Elements have to be of same class')
        elseif g>0 && g==n  % only numbers
                A = cell2mat(A);
        end
end
u = unique_bc2(A);
[tf, ind] = ismember_bc2(A, u);
if ind>0
        k = mode(ind);
        m = u(k);
else
        m = NaN;
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
% The COBRAToolbox: testShowprogress.m
%
% Purpose:
%     - testShowprogress tests the showprogress function (0: silent, 1: textprogressbar, 2: waitbar)
%
% Author:
%     - Laurent Heirendt - March 2017

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testShowprogress'));
cd(fileDir);

global WAITBAR_TYPE

for m = 0:2
    WAITBAR_TYPE = m;
    fprintf('Testing show progress, mode = %i:\n', m)

    % run multiple tests with the textprogressbar (WAITBAR_TYPE = 1)
    if WAITBAR_TYPE == 1
        kend = 3;
    else
        kend = 1;
    end

    for k = 1:kend
        showprogress(0,'Testing showprogress ...');
        for i = 1:k * 10
            showprogress(i/(k * 10));
            pause(0.02)
        end
    end
    if WAITBAR_TYPE == 1
        fprintf('\n')
    end
end

% update with 1 input argument
WAITBAR_TYPE = 1;
showprogress(0);

% update with WAITBAR_TYPE=[]
WAITBAR_TYPE = [];
showprogress(0,'Testing showprogress ...');

% reset the default WAITBAR_TYPE
WAITBAR_TYPE = 1;
if ~isempty(strfind(getenv('HOME'), 'jenkins'))
    WAITBAR_TYPE = 0;
end

% return to original directory
cd(currentDir);

author: opencobra

repo: cobratoolbox


text:
function z = snowflake(n,a)
%SNOWFLAKE Koch Snowflake Curve
%   Z = SNOWFLAKE(N,A) is a closed curve in the complex plane
%   with 3*2^N+1 points. N is a nonnegative integer and A is a
%   complex number with |A| < 1 and |1-A| < 1.
%   Default is A = 1/2 + i*sqrt(3)/6.
%
%   % Examples
%   plot(snowflake(10)), axis equal
%   plot(snowflake(10,0.45+0.35i)), axis equal

%   Author: Jonas Lundgren <splinefit@gmail.com> 2010

if nargin < 1, n = 0; end
if nargin < 2, a = 1/2 + sqrt(-3)/6; end

% Constants
b = 1 - a;
c = 1/2 + sqrt(-3)/2;
d = 1 - c;

% Generate point sequence
z = 1;
for k = 1:n
    z = conj(z);
    z = [a*z; b*z+a];
end

% Close snowflake
z = [0; z; 1-c*z; 1-c-d*z];

author: Sable

repo: mcbench-benchmarks


text:
function root = newton(f,g,xo,epsilon)
% Program : newton.m
%
% Purpose : find the root of a given single variable function within the
%           given interval
%
% Author :  Aamir Alaud Din
% Date :    20.09.2013
%
% Inputs :  At least three input arguments are required viz., function,
%           derivative of the function, and initial guess
%
% Syntax :  root = newton(f,g,xo)
%           where f is an inline function and is continuous in the given
%           interval
%           e.g., f = inline('sin(x)','x');
%           g is the first derivative of f in the form of inline function
%           and is continuous in the given interval
%           e.g., g = inline('cos(x)','x');
%           xo is the initial guess which lies in the given interval
%
% Example:  root = newton(f,g,5)
%           The fourth input argument is the stopping criteria
%           If the function value deceeds epsilon, the loop terminates
%           giving the root. Default is 1e-6.
%
% Syntax :  root = newton(f,g,xo,epsilon)
%           where f, g, xo, and epsilon are explained above.
%
% Example:  root = bisection(f,g,-2,1e-9)

if (nargin < 3)
    error('Less input arguments provided.');

elseif (nargin == 3)
    epsilon = 1e-6;
    Iter = 0;
    fprintf('Iteration\t\t   xo\t\t   x1\t\t\tf(x)\t\t\tg(x)\n');
    fprintf('=========\t\t ======\t\t ======\t\t===========\t\t===========\t\t\n');
    while (abs(f(xo)) >= epsilon)
        Iter = Iter + 1;
        if (g(xo) == 0)
            fprintf('g(xo) is zero on iteration number %d',Iter);
            fprintf('\n');
            disp('We can''t compute the root');
            return;
        end
        x1 = xo - (f(xo)/g(xo));
        fprintf('%3d',Iter);
        fprintf('%20.4f',xo);
        fprintf('%12.4f',x1);
        fprintf('%16.4f',f(xo));
        fprintf('%16.4f',g(xo));
        fprintf('\n');
        xo = x1;
    end
    root = xo;
elseif (nargin == 4)
    Iter = 0;
    fprintf('Iteration\t\t   xo\t\t   x1\t\t\tf(x)\t\t\tg(x)\n');
    fprintf('=========\t\t ======\t\t ======\t\t===========\t\t===========\t\t\n');
    while (abs(f(xo)) >= epsilon)
        Iter = Iter + 1;
        if (g(xo) == 0)
            fprintf('g(xo) is zero on iteration number %d',Iter);
            fprintf('\n');
            disp('We can''t compute the root');
            return;
        end
        x1 = xo - (f(xo)/g(xo));
        fprintf('%3d',Iter);
        fprintf('%20.4f',xo);
        fprintf('%12.4f',x1);
        fprintf('%16.4f',f(xo));
        fprintf('%16.4f',g(xo));
        fprintf('\n');
        xo = x1;
    end
    root = xo;
end


text:
% Make an HMM with Gaussian observations
%   X1 -> X2
%   |     |
%   v     v
%   Y1    Y2

intra = zeros(2);
intra(1,2) = 1;
inter = zeros(2);
inter(1,1) = 1;
n = 2;

Q = 2; % num hidden states
O = 2; % size of observed vector
ns = [Q O];
bnet = mk_dbn(intra, inter, ns, 'discrete', 1, 'observed', 2);

prior0 = normalise(rand(Q,1));
transmat0 = mk_stochastic(rand(Q,Q));
mu0 = rand(O,Q);
Sigma0 = repmat(eye(O), [1 1 Q]);
bnet.CPD{1} = tabular_CPD(bnet, 1, prior0);
%% we set the cov prior to 0 to give same results as HMM toolbox
%bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', mu0, 'cov', Sigma0, 'cov_prior_weight', 0);
bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', mu0, 'cov', Sigma0);
bnet.CPD{3} = tabular_CPD(bnet, 3, transmat0);

T = 5; % fixed length sequences

engine = {};
engine{end+1} = smoother_engine(jtree_2TBN_inf_engine(bnet));
engine{end+1} = smoother_engine(hmm_2TBN_inf_engine(bnet));
engine{end+1} = hmm_inf_engine(bnet);
engine{end+1} = jtree_unrolled_dbn_inf_engine(bnet, T);
%engine{end+1} = frontier_inf_engine(bnet);
engine{end+1} = bk_inf_engine(bnet, 'clusters', {[1]});
engine{end+1} = jtree_dbn_inf_engine(bnet);

inf_time = cmp_inference_dbn(bnet, engine, T);

ncases = 2;
max_iter = 2;
[learning_time, CPD, LL, cases] = cmp_learning_dbn(bnet, engine, T, 'ncases', ncases, 'max_iter', max_iter);

% Compare to HMM toolbox

data = zeros(O, T, ncases);
for i=1:ncases
  data(:,:,i) = cell2num(cases{i}(bnet.observed, :));
end

tic
[LL2, prior2, transmat2, mu2, Sigma2] = mhmm_em(data, prior0, transmat0, mu0, Sigma0, [],  'max_iter', max_iter);
t=toc;
disp(['HMM toolbox took ' num2str(t) ' seconds '])

e = 1;
assert(approxeq(prior2, CPD{e,1}.CPT))
assert(approxeq(mu2, CPD{e,2}.mean))
assert(approxeq(Sigma2, CPD{e,2}.cov))
assert(approxeq(transmat2, CPD{e,3}.CPT))
assert(approxeq(LL2, LL{e}))

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function y=check_states(x,s,t_start,t)
x=x+s.*(t-t_start);
dx=min(diff(x));
if (dx<-1e-6)
	y=0;
	dx
else
	y=1;
end;

repo: Aero-matlab


text:
function [ ALL_HYPS, CAN_POOL ] = groundTruthHypothesisNewFunc( gnd, omap, gc, cn, rotImg, typenum, bufName)
%GROUNDTRUTHHYPOTHESISNEWFUNC Summary of this function goes here
%   Detailed explanation goes here
global config;
%% compute the checking direction on omap
% [candiSetXYZ, ~] = icosahedron2sphere(6);
[candiSetXYZ, ~] = getUniformVector( 6 );
[ohei, owid, ~] = size(omap);
candiSetUV = uv2coords(xyz2uvN(candiSetXYZ), owid, ohei);
candiInd = sub2ind([ohei owid], candiSetUV(:,2), candiSetUV(:,1));
% convert gc to omap
gndGC = zeros(ohei, owid, 3);
gndGC(:,:,1) = max(gc(:,:,5),gc(:,:,6));
gndGC(:,:,2) = max(gc(:,:,1),gc(:,:,3));
gndGC(:,:,3) = max(gc(:,:,2),gc(:,:,4));
normGndGC = sum(gndGC(candiInd))+sum(gndGC(candiInd+1*ohei*owid))+sum(gndGC(candiInd+2*ohei*owid));
% omap
gndOmap = omap;
gndOmap(gndOmap>0) = 1;
normGndOmap = sum(gndOmap, 3);
gndOmap = gndOmap./(repmat(normGndOmap+0.0001, [1 1 3]));
normGndOmap = sum(gndOmap(candiInd))+sum(gndOmap(candiInd+1*ohei*owid))+sum(gndOmap(candiInd+2*ohei*owid));
numPoint = size(candiSetXYZ,1);

% merge
gndMerge = zeros(ohei, owid, 3);
gndMerge(1:ohei/2,:,:) = gndOmap(1:ohei/2,:,:);
gndMerge(ohei/2+1:ohei, :, :) = gndGC(ohei/2+1:ohei, :, :);
normMegMap = sum(gndMerge(candiInd))+sum(gndMerge(candiInd+1*ohei*owid))+sum(gndMerge(candiInd+2*ohei*owid));

% optimal merge
gndOptim = zeros(ohei, owid, 3);
gndOptim(1:685,:,:) = gndOmap(1:685,:,:);
gndOptim(686:ohei, :, :) = gndGC(686:ohei, :, :);
normOptMap = sum(gndOptim(candiInd))+sum(gndOptim(candiInd+1*ohei*owid))+sum(gndOptim(candiInd+2*ohei*owid));

%%
room = gnd(1);
% compute room hyp omap feature
hypVp = false(numPoint,3);
align = room.align;
valid1 = insideCone(align([6;5;8;7],:), candiSetXYZ, 0);
valid2 = insideCone(align([7;8;4;3],:), candiSetXYZ, 0);
valid3 = insideCone(align([6;7;3;2],:), candiSetXYZ, 0);
valid4 = insideCone(align([3;4;1;2],:), candiSetXYZ, 0);
valid5 = insideCone(align([5;6;2;1],:), candiSetXYZ, 0);
valid6 = insideCone(align([8;5;1;4],:), candiSetXYZ, 0);
hypVp(valid1,1) = true; hypVp(valid4,1) = true;
hypVp(valid2,2) = true; hypVp(valid5,2) = true;
hypVp(valid3,3) = true; hypVp(valid6,3) = true;
response = sum(gndOmap(candiInd(hypVp(:,1))+0*ohei*owid)) ...
         + sum(gndOmap(candiInd(hypVp(:,2))+1*ohei*owid)) ...
         + sum(gndOmap(candiInd(hypVp(:,3))+2*ohei*owid));
room.omapScr = response/normGndOmap;
response = sum(gndGC(candiInd(hypVp(:,1))+0*ohei*owid)) ...
         + sum(gndGC(candiInd(hypVp(:,2))+1*ohei*owid)) ...
         + sum(gndGC(candiInd(hypVp(:,3))+2*ohei*owid));
room.gcScr = response/normGndGC;
response = sum(gndMerge(candiInd(hypVp(:,1))+0*ohei*owid)) ...
         + sum(gndMerge(candiInd(hypVp(:,2))+1*ohei*owid)) ...
         + sum(gndMerge(candiInd(hypVp(:,3))+2*ohei*owid));
room.mgScr = response/normMegMap;
response = sum(gndOptim(candiInd(hypVp(:,1))+0*ohei*owid)) ...
         + sum(gndOptim(candiInd(hypVp(:,2))+1*ohei*owid)) ...
         + sum(gndOptim(candiInd(hypVp(:,3))+2*ohei*owid));
room.opScr = response/normOptMap;

%%
% load('./object_classifier/objectclassifier.mat');
load(config.classifierName);
objects = gnd(2:end);
[ scores, valid ] = objectClassifier(room, objects, B);
obj_xyz = zeros(length(objects),6);
for i = 1:length(objects)
    obj_xyz(i,:) = [objects(i).align(1,:) objects(i).align(7,:)];
end
sel_hyps.objects = objects;
sel_hyps.scores = scores;
sel_hyps.valid = valid;
sel_hyps.obj_xyz = obj_xyz;

CAN_POOL{1}.sel_hyps = sel_hyps;
CAN_POOL{1}.room = room;

sel_hyps = repmat(struct('fixed',false,'selID',[]), typenum, 1);
for i = 1:length(objects)
    if objects(i).objtype>12
        continue;
    end
    sel_hyps(objects(i).objtype).fixed = true;
    sel_hyps(objects(i).objtype).selID = [sel_hyps(objects(i).objtype).selID i];
end
ALL_HYPS(1).sel_hyps = sel_hyps;

[objfea, roomfea, anglesid] = compObjHypsFeature( rotImg, omap, gc, cn, CAN_POOL{1}, typenum, bufName );
CAN_POOL{1}.sel_hyps.objfea = objfea;
CAN_POOL{1}.room.roomfea = roomfea;
CAN_POOL{1}.sel_hyps.anglesid = anglesid;

[ sceneImgFea ] = compSceneHypsFeatureA( ALL_HYPS, CAN_POOL{1}.sel_hyps, CAN_POOL{1}.room, omap, gc, cn );
ALL_HYPS(1).sceneImgFea = sceneImgFea;

end

author: zouchuhang

repo: LayoutNet


text:
function kd=kronDelta(i,j)

kd=double(i==j);

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function r = mrdivide(p,q)
%RDIVIDE this function divides the x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

if isa(p, 'CellVector') && (~isa(q, 'CellVector'))
    r = (1/q)*p;
else
    error('FVMtool: Wrong use of mrdivide for a face variable. Try using ./ instead.');
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
%
%
% Copyright (c) 2012 University of Crete - Computer Science Department (UOC-CSD)
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function p = wrappednormcdf(x, m, s, N)

    % N=10 gives a log error smaller than 1e-12 for s=8
    % with s=8 distribution is "almost" uniform (std < 1e-12)
    if nargin<4; N=10; end

    p = 0;
    for k=-N:N
        p = p + normcdf(x+2.*pi.*k, m, s);
    end

    p = p - N;

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [out] = naturallabel(n)
% Returns a natural label `idv` of `n` carbons.
% Assumes 1.1% C13
%
% USAGE:
%
%    [out] = naturallabel(n)
%
% INPUT:
%    n:      size of label
%
% OUTPUT:
%    out:    natural label idv of n carbons

if n <= 0
    out = 1;
    return;
end

out = zeros(2^n,1);
for i = 0:(2^n-1)
    t = dec2bin(i,n);
    c13 = sum(t-48); % subtract 48 for the '0' offset.
    c12 = n-c13;
    out(i+1) = .989^c12*.011^c13;
end


text:
classdef WhitespaceTokenizer < bert.tokenizer.internal.Tokenizer
    % WhitespaceTokenizer   The simplest type of tokenization, split on
    % whitespace characters.

    % Copyright 2020 The MathWorks, Inc.

    methods
        function text = tokenize(~,text)
            % tokens = tokenize(tok,str)   Returns an array of tokens formed
            %                              by splitting str on whitespace.
            arguments
                ~
                text
            end
            text = strip(text);
            text = split(text).';
        end
    end
end

author: matlab-deep-learning

repo: transformer-models


text:
%#eml
function [qdd,Fzf,Fzr,Flr,Ftr,Flf,Ftf,vlr,vtr,vlf,vtf,r0,r7,IMU,fx,hx]         = mdl_motorcycle(u,x,p)
% from nonlinear state space to generalized coordinates
qd=x(1:13);
q =x(14:26);
% coordinates
qx = q(1)  ; qxd = qd(1) ; qy = q(2)  ; qyd = qd(2) ; qz = q(3)  ; qzd = qd(3);
q0 = q(4)  ; q0d = qd(4) ; q1 = q(5)  ; q1d = qd(5) ; q2 = q(6)  ; q2d = qd(6);
q4 = q(7)  ; q4d = qd(7) ; q5 = q(8)  ; q5d = qd(8) ; q6 = q(9)  ; q6d = qd(9);
q7 = q(10) ; q7d = qd(10); qf = q(11) ; qfd = qd(11);
q8 = q(12) ; q8d = qd(12); q9 = q(13) ; q9d = qd(13);

% parameters
a1 = p(1)  ; l2 = p(11)  ; l3 = p(23) ; m5 = p(35)  ; a6 = p(47) ;
b1 = p(2)  ; m2 = p(12)  ; m3 = p(24) ; x5 = p(36)  ; b6 = p(48) ;
m1 = p(3)  ; x2 = p(13)  ; x3 = p(25) ; z5 = p(37)  ; m6 = p(49) ;
i1 = p(4)  ; i2 = p(14)  ; z3 = p(26) ; i5 = p(38)  ; i6 = p(50) ;
j1 = p(5)  ; t2 = p(15)  ; i3 = p(27) ; t5 = p(39)  ; j6 = p(51) ;
f1 = p(6)  ; b2 = p(16)  ; k3 = p(28) ; b5 = p(40)  ; f6 = p(52) ;
e1 = p(7)  ; p2 = p(17)  ; d3 = p(29) ; p5 = p(41)  ; e6 = p(53) ;
k1 = p(8)  ; d2 = p(18)  ; l4 = p(30) ; d5 = p(42)  ; k6 = p(54) ;
d1 = p(9)  ; e2 = p(19)  ; m4 = p(31) ; e5 = p(43)  ; d6 = p(55) ;
t1 = p(10) ; k2 = p(20)  ; x4 = p(32) ; k5 = p(44)  ; t6 = p(56) ;
             n2 = p(21)  ; z4 = p(33) ; n5 = p(45)  ;
             f2 = p(22)  ; i4 = p(34) ; f5 = p(46)  ;

u7=p(57); v7=p(58); w7=p(59); x7=p(60); y7=p(61); z7=p(62);
I1x=j1;I1y=i1;I1z=j1;I2x=i2;I2y=i2;I2z=i2;I3x=i3;I3y=i3;I3z=i3;
I4x=i4;I4y=i4;I4z=i4;I5x=i5;I5y=i5;I5z=i5;I6x=j6;I6y=i6;I6z=j6;

% trigonometry
s0 = sin(q0); c0 = cos(q0); s1 = sin(q1); c1 = cos(q1); s2 = sin(q2);c2 = cos(q2);
s4 = sin(q4); c4 = cos(q4); s5 = sin(q5); c5 = cos(q5); s6 = sin(q6);c6 = cos(q6);
s7 = sin(q7); c7 = cos(q7); s8 = sin(q8); c8 = cos(q8); s9 = sin(q9);c9 = cos(q9);

% basic rotations
R0 = [c0 -s0 0;s0 c0 0;0 0 1]; R0q = [-s0 -c0 0;c0 -s0 0;0 0 0]; %rear yaw
R1 = [1 0 0;0 c1 -s1;0 s1 c1]; R1q = [0 0 0;0 -s1 -c1;0 c1 -s1]; %rear roll
R2 = [c2 0 s2;0 1 0;-s2 0 c2]; R2q = [-s2 0 c2;0 0 0;-c2 0 -s2]; %rear pitch
R4 = [c4 -s4 0;s4 c4 0;0 0 1]; R4q = [-s4 -c4 0;c4 -s4 0;0 0 0]; %front yaw
R5 = [1 0 0;0 c5 -s5;0 s5 c5]; R5q = [0 0 0;0 -s5 -c5;0 c5 -s5]; %front roll
R6 = [c6 0 s6;0 1 0;-s6 0 c6]; R6q = [-s6 0 c6;0 0 0;-c6 0 -s6]; %front pitch
R7 = [c7 0 s7;0 1 0;-s7 0 c7]; R7q = [-s7 0 c7;0 0 0;-c7 0 -s7]; %swingarm angle
R8 = [c8 0 s8;0 1 0;-s8 0 c8]; R8q = [-s8 0 c8;0 0 0;-c8 0 -s8]; %rear wheel angle
R9 = [c9 0 s9;0 1 0;-s9 0 c9]; R9q = [-s9 0 c9;0 0 0;-c9 0 -s9]; %front wheel angle

% inertia matrix;
sa = m3*(x3-l3)-(m1+m2)*l3;
sb = m3*z3;
sc = m2*(x2-l2)-m1*l2;
sd = m4*x4+m5*(l4+x5)+m6*l4;
se = m4*z4+m5*z5;
sf = -m5-m6;
sg = m1*l2^2+m2*(x2-l2)^2+I2y;
sh = m6*l4+m5*(l4+x5);
si = 2*m5*z5;
sj = m4*(x4^2+z4^2)+m5*((l4+x5)^2+z5^2)+I4y+I5y+m6*l4^2;
sk = m5*(l4+x5)+m6*l4;
sl = (m1+m2)*l3^2+m3*((-l3+x3)^2+z3^2)+I3y;
sm = l3*(m1*l2+m2*(l2-x2));
sn = m3*(z3^2+x3^2)+l3*((m1+m2+m3)*l3-2*m3*x3)+I3y;
so = l3*sc;
sp = I5z-I5x+I4z-I4x;

M          = (m1+m2+m3+m4+m5+m6)*eye(13);
M(1:3,4:11)= [R0q*(R1*(R2*[sa;0;sb]+R7*[sc;0;0])) R0*(R1q*(R2*[sa;0;sb]+R7*[sc;0;0])) R0*(R1*(R2q*[sa;0;sb])) R4q*(R5*(R6*[sd;0;se+sf*qf])) R4*(R5q*(R6*[sd;0;se+sf*qf])) R4*(R5*(R6q*[sd;0;se+sf*qf])) R0*(R1*(R7q*[sc;0;0])) R4*(R5*(R6*[0;0;sf]))];
M(4,4)     = [1 1 0]*R1.^2*((m1*(R2*[-l3;0;0]+R7*[-l2;0;0]).^2+m2*(R2*[-l3;0;0]+R7*[-l2+x2;0;0]).^2+m3*(R2*[-l3+x3;0;z3]).^2))+I3y+I2y+I1y+c1^2*(I1x-I1y+I2x-I2y+I3x-I3y+c7^2*(I2z-I2x)+c2^2*(I3z-I3x));
M(4,5)     = -c1*(m1*(-c2*l3-c7*l2)*(s2*l3+s7*l2)+m2*(-c2*l3-c7*l2+c7*x2)*(s2*l3+s7*l2-s7*x2)+m3*(-c2*l3+c2*x3+s2*z3)*(s2*l3-s2*x3+c2*z3))+c1*s7*c7*(I2z-I2x)+c1*s2*c2*(I3z-I3x);
M(4,6)     = s1*(sn-cos(q2-q7)*so);
M(4,10)    = s1*(sg+cos(q2-q7)*sm);
M(4,12)    = I1y*s1;
M(5,5)     = m1*(s2*l3+s7*l2)^2+m2*(s2*l3-s7*(-l2+x2))^2+m3*(-s2*(-l3+x3)+c2*z3)^2+c7^2*(I2x-I2z)+c2^2*(I3x-I3z)+I3z+I2z+I1x;
M(6,6)     = sl;
M(6,10)    = cos(q2-q7)*sm;
M(7,7)     = [1 1 0]*(R5.^2*(m4*(R6*[x4;0;z4]).^2+m5*(R6*[l4+x5;0;-qf+z5]).^2+m6*(R6*[l4;0;-qf]).^2))+I6y+I5y+I4y+c5^2*(I6x-I6y+I5x-I5y+I4x-I4y+c6^2*sp);
M(7,8)     = -c5*(c6*s6*(m4*(z4^2-x4^2)+m5*((z5-qf)^2-(l4+x5)^2)+m6*(qf^2-l4^2))+cos(2*q6)*(m4*x4*z4+m5*(l4+x5)*(z5-qf)-m6*l4*qf))+c5*s6*c6*sp;
M(7,9)     = s5*(sj-sf*qf^2-si*qf);
M(7,11)    = s5*sk;
M(7,13)    = I6y*s5;
M(8,8)     = m4*(c6*z4-s6*x4)^2+m5*(-s6*(l4+x5)+c6*(z5-qf))^2+m6*(-s6*l4-c6*qf)^2-c6^2*sp+I4z+I5z+I6x;
M(9,9)     = sj-sf*qf^2-si*qf;
M(9,11)    = sh;
M(10,10)   = sg;
M(11,11)   = -sf;
M(12,12)   = I1y;
M(13,13)   = I6y;
M          = M.'+triu(M,1);

% Christoffel matrix
s=[m6*(-s6*l4-c6*qf)*c6;-s6*(l4+x5)+c6*(-qf+z5);c5*(c6*s6*(m5*(2*qf-2*z5)+2*m6*qf)+cos(2*q6)*(-m5*(l4+x5)-m6*l4));c6*(l4+x5)+s6*(-qf+z5);s1*sin(q2-q7)*sm;s1*sin(q2-q7)*so;s2*l3-s7*(-l2+x2);m1*(s2*l3+s7*l2)*c7*l2;c1*(sg+cos(q2-q7)*sm);-c2*l3-c7*l2+c7*x2;s2*l3+s7*l2-s7*x2;c1*(m1*s7*l2*(s2*l3+s7*l2)+m1*...
   (-c2*l3-c7*l2)*c7*l2+m2*(s7*l2-s7*x2)*(s2*l3+s7*l2-s7*x2)+m2*(-c2*l3-c7*l2+c7*x2)*(c7*l2-c7*x2));-c2*l3+c7*(-l2+x2);m6*(-s6*l4-c6*qf)*(-c6*l4+s6*qf);m5*(-s6*(l4+x5)+c6*(-qf+z5))*(-c6*(l4+x5)-s6*(-qf+z5));m4*(-s6*x4+c6*z4)*(-c6*x4-s6*z4);c5*(sj-sf*qf^2-si*qf);m4*x4*z4+m5*(l4+x5)*(-qf+z5)-m6*l4*qf;...
   m4*(z4^2-x4^2)+m5*((-qf+z5)^2-(l4+x5)^2)+m6*(qf^2-l4^2);c5*(-s6^2*(m4*(z4^2-x4^2)+m5*((-qf+z5)^2-(l4+x5)^2)+m6*(qf^2-l4^2))+c6^2*(m4*(z4^2-x4^2)+m5*((-qf+z5)^2-(l4+x5)^2)+m6*(qf^2-l4^2))-2*sin(2*q6)*(m4*x4*z4+m5*(l4+x5)*(-qf+z5)-m6*l4*qf));-c6*sd-s6*(se+sf*qf);-s6*sd+c6*(se+sf*qf);c5*(I6x-I6y+I5x...
   -I5y+I4x-I4y+c6^2*sp)*s5;s5*(m4*(-s6*x4+c6*z4)^2+m5*(-s6*(l4+x5)+c6*(-qf+z5))^2+m6*(-s6*l4-c6*qf)^2)*c5;-s2*(-l3+x3)+c2*z3;m3*(-s2*(-l3+x3)+c2*z3)*(-c2*(-l3+x3)-s2*z3);m1*(s2*l3+s7*l2)*c2*l3;c1*(sn-cos(q2-q7)*so);-c2*l3+c2*x3+s2*z3;s2*l3-s2*x3+c2*z3;c1*(m1*s2*l3*(s2*l3+s7*l2)+m1*(-c2*l3-c7*l2)*c2*...
   l3+m2*s2*l3*(s2*l3+s7*l2-s7*x2)+m2*(-c2*l3-c7*l2+c7*x2)*c2*l3+m3*(s2*l3-s2*x3+c2*z3)^2+m3*(-c2*l3+c2*x3+s2*z3)*(c2*l3-c2*x3-s2*z3));c1*(I1x-I1y+I2x-I2y+I3x-I3y+c7^2*(I2z-I2x)+c2^2*(I3z-I3x))*s1;s1*(m1*(s2*l3+s7*l2)^2+m2*(s2*l3-s7*(-l2+x2))^2+m3*(-s2*(-l3+x3)+c2*z3)^2)*c1;-s2*sa+c2*sb-s7*sc];
C=[0 0 0 (-c0*(c2*sa+s2*sb+c7*sc)-s0*s1*s(34))*q0d+c0*c1*s(34)*q1d+(-s0*(-s2*sa+c2*sb)+c0*s1*(-c2*sa-s2*sb))*q2d+(s0*s7*sc-c0*s1*c7*sc)*q7d c0*c1*s(34)*q0d-s0*s1*s(34)*q1d+s0*c1*(-c2*sa-s2*sb)*q2d-s0*c1*c7*sc*q7d (-s0*(-s2*sa+c2*sb)+c0*s1*(-c2*sa-s2*sb))*q0d+s0*c1*(-c2*sa-s2*sb)*q1d+(c0*(-c2*sa-s2*...
   sb)+s0*s1*(s2*sa-c2*sb))*q2d (-c4*(c6*sd+s6*(se+sf*qf))-s4*s5*s(22))*q4d+c4*c5*s(22)*q5d+(-s4*s(22)+c4*s5*s(21))*q6d+(-s4*s6*sf+c4*s5*c6*sf)*qfd c4*c5*s(22)*q4d-s4*s5*s(22)*q5d+s4*c5*s(21)*q6d+s4*c5*c6*sf*qfd (-s4*s(22)+c4*s5*s(21))*q4d+s4*c5*s(21)*q5d+(c4*s(21)+s4*s5*(s6*sd-c6*(se+sf*qf)))*q6d+...
   (c4*c6*sf-s4*s5*s6*sf)*qfd (s0*s7*sc-c0*s1*c7*sc)*q0d-s0*c1*c7*sc*q1d+(-c0*c7*sc+s0*s1*s7*sc)*q7d (-s4*s6*sf+c4*s5*c6*sf)*q4d+s4*c5*c6*sf*q5d+(c4*c6*sf-s4*s5*s6*sf)*q6d 0 0 ;0 0 0 (-s0*(c2*sa+s2*sb+c7*sc)+c0*s1*s(34))*q0d+s0*c1*s(34)*q1d+(c0*(-s2*sa+c2*sb)+s0*s1*(-c2*sa-s2*sb))*q2d+(-c0*s7*sc-s0*...
   s1*c7*sc)*q7d s0*c1*s(34)*q0d+c0*s1*s(34)*q1d-c0*c1*(-c2*sa-s2*sb)*q2d+c0*c1*c7*sc*q7d (c0*(-s2*sa+c2*sb)+s0*s1*(-c2*sa-s2*sb))*q0d-c0*c1*(-c2*sa-s2*sb)*q1d+(s0*(-c2*sa-s2*sb)-c0*s1*(s2*sa-c2*sb))*q2d (-s4*(c6*sd+s6*(se+sf*qf))+c4*s5*s(22))*q4d+s4*c5*s(22)*q5d+(c4*s(22)+s4*s5*s(21))*q6d+(c4*s6*sf+...
   s4*s5*c6*sf)*qfd s4*c5*s(22)*q4d+c4*s5*s(22)*q5d-c4*c5*s(21)*q6d-c4*c5*c6*sf*qfd (c4*s(22)+s4*s5*s(21))*q4d-c4*c5*s(21)*q5d+(s4*s(21)-c4*s5*(s6*sd-c6*(se+sf*qf)))*q6d+(s4*c6*sf+c4*s5*s6*sf)*qfd (-c0*s7*sc-s0*s1*c7*sc)*q0d+c0*c1*c7*sc*q1d+(-s0*c7*sc-c0*s1*s7*sc)*q7d (c4*s6*sf+s4*s5*c6*sf)*q4d-c4*...
   c5*c6*sf*q5d+(s4*c6*sf+c4*s5*s6*sf)*q6d 0 0 ;0 0 0 0 -c1*s(34)*q1d-s1*(-c2*sa-s2*sb)*q2d+s1*c7*sc*q7d -s1*(-c2*sa-s2*sb)*q1d+c1*(s2*sa-c2*sb)*q2d 0 -c5*s(22)*q5d-s5*s(21)*q6d-s5*c6*sf*qfd -s5*s(21)*q5d+c5*(s6*sd-c6*(se+sf*qf))*q6d-c5*s6*sf*qfd s1*c7*sc*q1d+c1*s7*sc*q7d -s5*c6*sf*q5d-c5*s6*sf*q6d ...
   0 0 ;0 0 0 (s(33)-s(32))*q1d+(m1*(-c2*l3-c7*l2)*s2*l3+m2*s(13)*s2*l3+m3*(c2*(-l3+x3)+s2*z3)*s(25)+1/2*s1^2*(2*s(27)+2*m2*s(7)*c2*l3+2*s(26))-c1^2*c2*(I3z-I3x)*s2)*q2d+(m1*(-c2*l3-c7*l2)*s7*l2-m2*s(13)*s7*(-l2+x2)+1/2*s1^2*(2*s(8)-2*m2*s(7)*c7*(-l2+x2))-c1^2*c7*(I2z-I2x)*s7)*q7d (s(33)-s(32))*q0d+...
   (s1*(m1*(-c2*l3-c7*l2)*(s2*l3+s7*l2)+m2*s(10)*s(11)+m3*s(29)*s(30))-s1*s7*c7*(I2z-I2x)-s1*s2*c2*(I3z-I3x))*q1d+(1/2*s(28)-1/2*s(31)+1/2*c1*c2^2*(I3z-I3x)-1/2*c1*s2^2*(I3z-I3x))*q2d+(1/2*s(9)-1/2*s(12)+1/2*c1*c7^2*(I2z-I2x)-1/2*c1*s7^2*(I2z-I2x))*q7d+1/2*I1y*c1*q8d (m1*(-c2*l3-c7*l2)*s2*l3+m2*s(13)...
   *s2*l3+m3*(c2*(-l3+x3)+s2*z3)*s(25)+1/2*s1^2*(2*s(27)+2*m2*s(7)*c2*l3+2*s(26))-c1^2*c2*(I3z-I3x)*s2)*q0d+(1/2*s(28)-1/2*s(31)+1/2*c1*c2^2*(I3z-I3x)-1/2*c1*s2^2*(I3z-I3x))*q1d+s1*sin(q2-q7)*so*q2d+(-1/2*s(5)-1/2*s(6))*q7d 0 0 0 (m1*(-c2*l3-c7*l2)*s7*l2-m2*s(13)*s7*(-l2+x2)+1/2*s1^2*(2*s(8)-2*m2*...
   s(7)*c7*(-l2+x2))-c1^2*c7*(I2z-I2x)*s7)*q0d+(1/2*s(9)-1/2*s(12)+1/2*c1*c7^2*(I2z-I2x)-1/2*c1*s7^2*(I2z-I2x))*q1d+(-1/2*s(5)-1/2*s(6))*q2d+s1*sin(q2-q7)*sm*q7d 0 1/2*I1y*c1*q1d 0 ;0 0 0 (-s(33)+s(32))*q0d+(-1/2*s(31)+1/2*c1*c2^2*(I3z-I3x)-1/2*c1*s2^2*(I3z-I3x)-1/2*s(28))*q2d+(-1/2*s(12)+1/2*c1*...
   c7^2*(I2z-I2x)-1/2*c1*s7^2*(I2z-I2x)-1/2*s(9))*q7d-1/2*I1y*c1*q8d (s(27)+m2*s(7)*c2*l3+s(26)-c2*(I3x-I3z)*s2)*q2d+(s(8)-m2*s(7)*c7*(-l2+x2)-c7*(I2x-I2z)*s7)*q7d (-1/2*s(31)+1/2*c1*c2^2*(I3z-I3x)-1/2*c1*s2^2*(I3z-I3x)-1/2*s(28))*q0d+(s(27)+m2*s(7)*c2*l3+s(26)-c2*(I3x-I3z)*s2)*q1d 0 0 0 (-1/2*s(12)...
   +1/2*c1*c7^2*(I2z-I2x)-1/2*c1*s7^2*(I2z-I2x)-1/2*s(9))*q0d+(s(8)-m2*s(7)*c7*(-l2+x2)-c7*(I2x-I2z)*s7)*q1d 0 -1/2*I1y*c1*q0d 0 ;0 0 0 (-m1*(-c2*l3-c7*l2)*s2*l3-m2*s(13)*s2*l3-m3*(c2*(-l3+x3)+s2*z3)*s(25)-1/2*s1^2*(2*s(27)+2*m2*s(7)*c2*l3+2*s(26))+c1^2*c2*(I3z-I3x)*s2)*q0d+(1/2*s(28)+1/2*s(31)-1/2*...
   c1*c2^2*(I3z-I3x)+1/2*c1*s2^2*(I3z-I3x))*q1d+(-1/2*s(6)+1/2*s(5))*q7d (1/2*s(28)+1/2*s(31)-1/2*c1*c2^2*(I3z-I3x)+1/2*c1*s2^2*(I3z-I3x))*q0d+(-s(27)-m2*s(7)*c2*l3-s(26)+c2*(I3x-I3z)*s2)*q1d 0 0 0 0 (-1/2*s(6)+1/2*s(5))*q0d+sin(q2-q7)*sm*q7d 0 0 0 ;0 0 0 0 0 0 (s(24)-s(23))*q5d+(m4*(c6*x4+s6*z4)*...
   (-s6*x4+c6*z4)+m5*s(4)*s(2)+m6*(c6*l4-s6*qf)*(-s6*l4-c6*qf)+1/2*s5^2*(2*s(16)+2*s(15)+2*s(14))-c5^2*c6*sp*s6)*q6d+(-m5*s(4)*s6-m6*(c6*l4-s6*qf)*s6+1/2*s5^2*(-2*m5*s(2)*c6-2*s(1)))*qfd (s(24)-s(23))*q4d+(s5*(c6*s6*s(19)+cos(2*q6)*s(18))-s5*s6*c6*sp)*q5d+(1/2*s(17)-1/2*s(20)+1/2*c5*c6^2*sp-1/2*c5*...
   s6^2*sp)*q6d+(1/2*c5*sk-1/2*s(3))*qfd+1/2*I6y*c5*q9d (m4*(c6*x4+s6*z4)*(-s6*x4+c6*z4)+m5*s(4)*s(2)+m6*(c6*l4-s6*qf)*(-s6*l4-c6*qf)+1/2*s5^2*(2*s(16)+2*s(15)+2*s(14))-c5^2*c6*sp*s6)*q4d+(1/2*s(17)-1/2*s(20)+1/2*c5*c6^2*sp-1/2*c5*s6^2*sp)*q5d+1/2*s5*(-2*sf*qf-si)*qfd 0 (-m5*s(4)*s6-m6*(c6*l4-s6*qf)...
   *s6+1/2*s5^2*(-2*m5*s(2)*c6-2*s(1)))*q4d+(1/2*c5*sk-1/2*s(3))*q5d+1/2*s5*(-2*sf*qf-si)*q6d 0 1/2*I6y*c5*q5d ;0 0 0 0 0 0 (-s(24)+s(23))*q4d+(-1/2*s(20)+1/2*c5*c6^2*sp-1/2*c5*s6^2*sp-1/2*s(17))*q6d+(-1/2*s(3)-1/2*c5*sk)*qfd-1/2*I6y*c5*q9d (s(16)+s(15)+s(14)+c6*sp*s6)*q6d+(-m5*s(2)*c6-s(1))*qfd ...
   (-1/2*s(20)+1/2*c5*c6^2*sp-1/2*c5*s6^2*sp-1/2*s(17))*q4d+(s(16)+s(15)+s(14)+c6*sp*s6)*q5d 0 (-1/2*s(3)-1/2*c5*sk)*q4d+(-m5*s(2)*c6-s(1))*q5d 0 -1/2*I6y*c5*q4d ;0 0 0 0 0 0 (-m4*(c6*x4+s6*z4)*(-s6*x4+c6*z4)-m5*s(4)*s(2)-m6*(c6*l4-s6*qf)*(-s6*l4-c6*qf)-1/2*s5^2*(2*s(16)+2*s(15)+2*s(14))+c5^2*c6*sp*...
   s6)*q4d+(1/2*s(17)+1/2*s(20)-1/2*c5*c6^2*sp+1/2*c5*s6^2*sp)*q5d+1/2*s5*(-2*sf*qf-si)*qfd (1/2*s(17)+1/2*s(20)-1/2*c5*c6^2*sp+1/2*c5*s6^2*sp)*q4d+(-s(16)-s(15)-s(14)-c6*sp*s6)*q5d (-sf*qf-1/2*si)*qfd 0 1/2*s5*(-2*sf*qf-si)*q4d+(-sf*qf-1/2*si)*q6d 0 0 ;0 0 0 (-m1*(-c2*l3-c7*l2)*s7*l2+m2*s(13)*s7*...
   (-l2+x2)-1/2*s1^2*(2*s(8)-2*m2*s(7)*c7*(-l2+x2))+c1^2*c7*(I2z-I2x)*s7)*q0d+(1/2*s(9)+1/2*s(12)-1/2*c1*c7^2*(I2z-I2x)+1/2*c1*s7^2*(I2z-I2x))*q1d+(-1/2*s(5)+1/2*s(6))*q2d (1/2*s(9)+1/2*s(12)-1/2*c1*c7^2*(I2z-I2x)+1/2*c1*s7^2*(I2z-I2x))*q0d+(-s(8)+m2*s(7)*c7*(-l2+x2)+c7*(I2x-I2z)*s7)*q1d (-1/2*s(5)+...
   1/2*s(6))*q0d-sin(q2-q7)*sm*q2d 0 0 0 0 0 0 0 ;0 0 0 0 0 0 (m5*s(4)*s6+m6*(c6*l4-s6*qf)*s6-1/2*s5^2*(-2*m5*s(2)*c6-2*s(1)))*q4d+(1/2*c5*sk+1/2*s(3))*q5d-1/2*s5*(-2*sf*qf-si)*q6d (1/2*c5*sk+1/2*s(3))*q4d+(m5*s(2)*c6+s(1))*q5d -1/2*s5*(-2*sf*qf-si)*q4d+(sf*qf+1/2*si)*q6d 0 0 0 0 ;0 0 0 1/2*I1y*c1*...
   q1d 1/2*I1y*c1*q0d 0 0 0 0 0 0 0 0 ;0 0 0 0 0 0 1/2*I6y*c5*q5d 1/2*I6y*c5*q4d 0 0 0 0 0 ];

% preparation on force calculation
g=9.81;
q70=-.2;
qf0=0.3;

r0d =[qxd+(-s0*(-c2*l3-c7*l2)-c0*(-c1*s1*a1-s1*(s2*l3+s7*l2-b1-c1*a1)))*q0d+s0*c1*(s2*l3+s7*l2-b1-c1*a1)*q1d+(c0*s2*l3+s0*s1*c2*l3)*q2d+(c0*s7*l2+s0*s1*c7*l2)*q7d+c0*(-b1-c1*a1)*q8d
      qyd+(c0*(-c2*l3-c7*l2)-s0*(-c1*s1*a1-s1*(s2*l3+s7*l2-b1-c1*a1)))*q0d-c0*c1*(s2*l3+s7*l2-b1-c1*a1)*q1d+(s0*s2*l3-c0*s1*c2*l3)*q2d+(s0*s7*l2-c0*s1*c7*l2)*q7d+s0*(-b1-c1*a1)*q8d
      qzd-s1*(s2*l3+s7*l2-b1-c1*a1)*q1d+c1*c2*l3*q2d+c1*c7*l2*q7d];
r7d =[qxd+(-s4*(c6*l4-s6*qf)-c4*(-c5*s5*a6-s5*(-s6*l4-c6*qf-b6-c5*a6)))*q4d+s4*c5*(-s6*l4-c6*qf-b6-c5*a6)*q5d+(c4*(-s6*l4-c6*qf)+s4*s5*(-c6*l4+s6*qf))*q6d+(-c4*s6-sin(q4)*s5*c6)*qfd+c4*(-b6-c5*a6)*q9d
      qyd+(c4*(c6*l4-s6*qf)-s4*(-c5*s5*a6-s5*(-s6*l4-c6*qf-b6-c5*a6)))*q4d-c4*c5*(-s6*l4-c6*qf-b6-c5*a6)*q5d+(s4*(-s6*l4-c6*qf)-c4*s5*(-c6*l4+s6*qf))*q6d+(-s4*s6+c4*s5*c6)*qfd+s4*(-b6-c5*a6)*q9d
      qzd-s5*(-s6*l4-c6*qf-b6-c5*a6)*q5d+c5*(-c6*l4+s6*qf)*q6d-c5*c6*qfd];

r0  = [0;0;qz+c1*(s2*l3+s7*l2-b1)-a1];
r7  = [0;0;qz-c5*(s6*l4+c6*qf+b6)-a6];
vtr =-s0*r0d(1)+c0*r0d(2);
vlr = c0*r0d(1)+s0*r0d(2)+.01;
vtf =-s4*r7d(1)+c4*r7d(2);
vlf = c4*r7d(1)+s4*r7d(2)+.01;
v1  = [c0*s2+s0*s1*c2;s0*s2-c0*s1*c2;c1*c2];
v2  = [c4*s6+s4*s5*c6;s4*s6-c4*s5*c6;c5*c6];

Rm3 = R0*(R1*R2); % frame
Rm4 = R4*R5*R6; % steering head

    % 1. find the axis about to rotate and the angle
    v1 = Rm3*[0;0;1];
    v2 = Rm4*[0;0;1];
    w1= cross(v1,v2);
    % 2. use Rodrigues rotation formula to get the new x and y vectors
    x_old = Rm3*[1;0;0];
    y_old = Rm3*[0;1;0];
    x_new = x_old*v1.'*v2+cross(w1,x_old)+w1*w1.'*x_old/(1+v1.'*v2);
    y_new = y_old*v1.'*v2+cross(w1,y_old)+w1*w1.'*y_old/(1+v1.'*v2);
    % 3. calculate angle between the vectors
    qxc = real(x_new.'*Rm4(:,1));
    qyc = real(y_new.'*Rm4(:,2));
    qxs = real(-x_new.'*Rm4(:,2));
    qys = real(y_new.'*Rm4(:,1));
% q3=atan2(qxs,qxc);
% q3=atan2(qxs,qyc);
q3=atan2(qys,qxc);
% q3=atan2(qys,qyc);

T3  = -u(2)*(v1+v2)/norm(v1+v2);

% forces/torques
Frs = k2*(q2-q7-q70)+b2*(q2d-q7d);
Ffs = k5*(qf-p5)+d5*qfd*(qfd<0)+e5*qfd*(qfd>0)-k3*(qf<b5)+k3*(qf>t5);

Frb = u(1)*(u(1)>0);
Ffb = u(1)*(u(1)<0);

Fzr = max((-k1*r0(3)-d1*r0d(3))*(r0(3)<0),0) ;
Ftr = t1*Fzr*atan(-vtr/vlr);
Flr = Fzr*min(max(f1*((b1+a1*cos(q1))*q8d/vlr-1),-1),1);

Fzf = max((-k1*r7(3)-d1*r7d(3))*(r7(3)<0),0) ;
Ftf = t6*Fzf*atan(-vtf/vlf);
Flf = Fzf*min(max(f6*((b6+a6*cos(q5))*q9d/vlf-1),-1),1);

T1  =  k3*(pi/2-dot(v1,v2))*cross(v1,v2);

%% experiment
sr    =  r0d(1:2);
Fr    =  -Fzr*(sr./((sr.'*sr)^4+1)^(1/8)+sr./((sr.'*sr)^2+1));
Flr =  c0*Fr(1)+s0*Fr(2);
Ftr = -s0*Fr(1)+c0*Fr(2);

ss    =  r7d(1:2);
Ff    =  -Fzf*(ss./((ss.'*ss)^4+1)^(1/8)+ss./((ss.'*ss)^2+1));
Flf =  c4*Ff(1)+s4*Ff(2);
Ftf = -s4*Ff(1)+c4*Ff(2);
%% end of experiment

% work done by forces/torques
Qfs = [0 0 0 0 0 0 0 0 0 0 -Ffs 0 0].';
Qrs = [0 0 0 0 0 -Frs 0 0 0 Frs 0 0 0].';
Qfb = [0 0 0 0 0 0 0 0 -Ffb 0 0 0 Ffb].';
Qrb = [0 0 0 0 0 0 0 0 0 -Frb 0 Frb 0].';
QT1 = [0 0 0 T1(3) c0*T1(1)+s0*T1(2) -s0*c1*T1(1)+c0*c1*T1(2)+s1*T1(3) -T1(3) -c4*T1(1)-s4*T1(2) s4*c5*T1(1)-c4*c5*T1(2)-s5*T1(3) 0 0 0 0].';
QT3 = [0 0 0 T3(3) c0*T3(1)+s0*T3(2) -s0*c1*T3(1)+c0*c1*T3(2)+s1*T3(3) -T3(3) -c4*T3(1)-s4*T3(2) s4*c5*T3(1)-c4*c5*T3(2)-s5*T3(3) 0 0 0 0].';
Qrz = [0 0 1 0 -s1*(s2*l3+s7*l2-b1) c1*c2*l3 0 0 0 c1*c7*l2 0 0 0                    ].'*Fzr;
Qfz = [0 0 1 0 0 0 0 s5*(s6*l4+c6*qf+b6) c5*(s6*qf-c6*l4) 0 -c5*c6 0 0               ].'*Fzf;
Qlf = [c4 s4 0 0 0 0 -s5*(s6*l4+c6*qf+b6) 0 -s6*l4-c6*qf 0 -s6 0 -c5*a6-b6           ].'*Flf;
Qlr = [c0 s0 0 -s1*(-s2*l3-s7*l2+b1) 0 s2*l3 0 0 0 s7*l2 0 -c1*a1-b1 0               ].'*Flr;
Qtf = [-s4 c4 0 0 0 0 c6*l4-s6*qf c5*(s6*l4+c6*qf+b6)+a6 s5*(c6*l4-s6*qf) 0 s5*c6 0 0].'*Ftf;
Qtr = [-s0 c0 0 -c2*l3-c7*l2 c1*(b1-s2*l3-s7*l2)+a1 -s1*c2*l3 0 0 0 -s1*c7*l2 0 0 0  ].'*Ftr;
Vq  = -g*[0 0 m1+m2+m3+m4+m5+m6 0 -m1*s1*(s2*l3+s7*l2)-m2*s1*(s2*l3-s7*(-l2+x2))-m3*s1*(-s2*(-l3+x3)+c2*z3) (m1+m2)*c1*c2*l3+m3*c1*(-c2*(-l3+x3)-s2*z3) 0 (-m4*s5*(-s6*x4+c6*z4)-m5*s5*(-s6*(l4+x5)+c6*(-qf+z5))-m6*s5*(-s6*l4-c6*qf)) m4*c5*(-c6*x4-s6*z4)+m5*c5*(-c6*(l4+x5)-s6*(-qf+z5))+m6*c5*(-c6*l4+s6*qf) m1*c1*c7*l2-m2*c1*c7*(-l2+x2) c5*c6*sf 0 0];
% Q   = Vq.'+Qfs+Qrs+Qrz+Qfz+QT1+QT3+Qlf+Qlr+Qtf;
Q   = Vq.'+Qfs+Qrs+Qrz+Qfz+QT1+QT3+Qtf+Qlf+Qlr+Qtr+Qfb+Qrb;
% Q   = Vq.'+Qfs+Qrs+Qrz+Qfz+QT1+QT3+Qflon+Qrlon+Qflat+Qrlat;

qdd = M\(Q-C*qd);

IMU = [0;0;0;0;0;0];

qxdd=qdd(1);
qydd=qdd(2);
qzdd=qdd(3);
q0dd=qdd(4);
q1dd=qdd(5);
q2dd=qdd(6);

su = sin(u7); cu = cos(u7); sv = sin(v7); cv = cos(v7); sw = sin(w7);cw = cos(w7);

% basic rotations
R0 = [c0 -s0 0;s0 c0 0;0 0 1]; R0q = [-s0 -c0 0;c0 -s0 0;0 0 0]; R0qq = [-c0 s0 0;-s0 -c0 0;0 0 0];
R1 = [1 0 0;0 c1 -s1;0 s1 c1]; R1q = [0 0 0;0 -s1 -c1;0 c1 -s1]; R1qq = [0 0 0;0 -c1 s1;0 -s1 -c1];
R2 = [c2 0 s2;0 1 0;-s2 0 c2]; R2q = [-s2 0 c2;0 0 0;-c2 0 -s2]; R2qq = [-c2 0 -s2;0 0 0;s2 0 -c2];

Ru = [cu -su 0;su cu 0;0 0 1];
Rv = [1 0 0;0 cv -sv;0 sv cv];
Rw = [cw 0 sw;0 1 0;-sw 0 cw];

% position vector
rm7 = [qx;qy;qz]+R0*(R1*(R2*[x7;y7;z7]));

% % velocity vector
% rm7d=jacobian(rm7,q)*qd
% rm7d1= [qxd;qyd;qzd]+q0d*R0q*(R1*(R2*[x7;y7;z7]))+R0*(q1d*R1q*(R2*[x7;y7;z7])+R1*(q2d*R2q*[x7;y7;z7]))
% rm7d2= [qxd;qyd;qzd]+R0q*(R1*(R2*q0d*[x7;y7;z7]))+R0*(R1q*(R2*q1d*[x7;y7;z7])+R1*(R2q*q2d*[x7;y7;z7]))
% temp=simple(rm7d1-rm7d2)
% acceleration vector
% rm7dd=(jacobian(rm7d,[q;qd])*[qd;qdd])
% rm7dd1=(jacobian(rm7d1,[q;qd])*[qd;qdd])% this vector is in global coordinates, but the sensor measures in local coordinates.
% rm7dd2=[qxdd;qydd;qzdd] ...
%     +R0qq*(R1*(R2*q0d^2*[x7;y7;z7]))   +R0q*(R1q*(R2*q0d*q1d*[x7;y7;z7])) +R0q*(R1*(R2q*q0d*q2d*[x7;y7;z7])) +R0q*(R1*(R2*q0dd*[x7;y7;z7])) ...
%     +R0q*(R1q*(R2*q0d*q1d*[x7;y7;z7])) +R0*(R1qq*(R2*q1d^2*[x7;y7;z7]))   +R0*(R1q*(R2q*q1d*q2d*[x7;y7;z7])) +R0*(R1q*(R2*q1dd*[x7;y7;z7])) ...
%     +R0q*(R1*(R2q*q0d*q2d*[x7;y7;z7])) +R0*(R1q*(R2q*q1d*q2d*[x7;y7;z7])) +R0*(R1*(R2qq*q2d^2*[x7;y7;z7]))   +R0*(R1*(R2q*q2dd*[x7;y7;z7]))
% temp=simple(rm7dd-rm7dd2)

% absolute position measured in local body coordinates
% - rm7loc = R2'*R1'*R0'*([qx;qy;qz]+R0*R1*R2*[x7;y7;z7])
% - rm7loc = R2'*R1'*R0'*[qx;qy;qz]+ R2'*R1'*R0'*R0*R1*R2*[x7;y7;z7]
% - rm7loc = R2'*R1'*R0'*[qx;qy;qz]+ [x7;y7;z7]
% rm7loc = R2.'*R1.'*R0.'*[qx;qy;qz]+ [x7;y7;z7]
% how can this be possible? this implies (or not?) that the velocity of the sensor is
% independent of its local position [x7;y7;z7].

% absolute velocity measured in local body coordinates
% rm7dloc8 = R2.'*(R1.'*(((R0.'*[qxd;qyd;qzd])) + (R0tR0q *(R1*(R2*q0d*[x7;y7;z7]))))+(R1tR1q*(R2*q1d*[x7;y7;z7])))+(R2tR2q*q2d*[x7;y7;z7])

% syms s0 s1 s2 c0 c1 c2
% rm7ddloc6 =  R2.'*R1.'*R0.'*[qxdd;qydd;qzdd]+ ...
%   [z7*(q2dd + (q1d^2*sin(2*q2))/2 + q0dd*s1 - (q0d^2*sin(2*q2)*c1^2)/2 + 2*q0d*q1d*cos(q1)*cos(q2)^2) - x7*(q0d^2*(c2^2 + s1^2*s2^2) + q1d^2*s2^2 + q2d^2 + 2*q0d*q2d*s1 + q0d*q1d*sin(2*q2)*c1) - y7*(c1*s1*s2*q0d^2 - 2*q1d*c2*s1*q0d + q1dd*s2 + q0dd*c1*c2);
%    x7*(-c1*s1*s2*q0d^2 - 2*q2d*c1*s2*q0d + q1dd*s2 + 2*q1d*q2d*c2 + q0dd*c1*c2) - y7*(q0d^2*c1^2 + q1d^2) + z7*((sin(2*q1)*c2*q0d^2)/2 + 2*q2d*c1*c2*q0d - q1dd*c2 + 2*q1d*q2d*s2 + q0dd*c1*s2);
%    y7*((sin(2*q1)*c2*q0d^2)/2 + 2*q1d*s1*s2*q0d + q1dd*c2 - q0dd*c1*s2) - z7*(q1d^2*c2^2 - q0d^2*(c1^2*c2^2 - 1) + q2d^2 + 2*q0d*q2d*s1 - q0d*q1d*sin(2*q2)*c1) - x7*(q2dd - (q1d^2*sin(2*q2))/2 + q0dd*s1 + (q0d^2*sin(2*q2)*c1^2)/2 + 2*q0d*q1d*c1*s2^2)]

rm7ddloc6 =  R2.'*(R1.'*(R0.'*[qxdd;qydd;qzdd+9.81]))+ ...
[ z7*(q2dd+q0dd*s1+c2*s2*(q1d^2-q0d^2*c1^2)+2*c1*c2^2*q0d*q1d) ...
- y7*(s2*(c1*s1*q0d^2+q1dd)+c2*(c1*q0dd-2*q1d*s1*q0d)) ...
- x7*(q1d^2*s2^2+q0d^2*s1^2*s2^2+q0d^2*c2^2+2*q0d*q2d*s1+2*c1*c2*q0d*q1d*s2+q2d^2); ...
  z7*((c1*(s1*q0d^2+2*q2d*q0d)-q1dd)*c2+(c1*q0dd+2*q1d*q2d)*s2) ...
+ x7*(-c1*s1*s2*q0d^2-2*c1*q2d*s2*q0d+q1dd*s2+c1*c2*q0dd+2*c2*q1d*q2d) ...
- y7*(c1^2*q0d^2+q1d^2); ...
  y7*(c2*(c1*s1*q0d^2+q1dd)+s2*(2*s1*q0d*q1d-c1*q0dd)) ...
- x7*(q2dd+q0dd*s1+c2*s2*(q0d^2*c1^2-q1d^2)+2*c1*q0d*q1d*s2^2) ...
- z7*(q0d^2*(1-c2^2*c1^2)+c2^2*q1d^2+2*q0d*(q2d*s1-c1*c2*q1d*s2)+q2d^2)]...
;

rm7ddloc6 = Rw.'*(Rv.'*(Ru.'*rm7ddloc6));

wrel=R2.'*R1.'*R0.'*[0;0;q0d]+R2.'*R1.'*[q1d;0;0]+R2.'*[0;q2d;0];
wrel = Rw.'*(Rv.'*(Ru.'*wrel));

IMU=[wrel;rm7ddloc6];

% from generalized coordinates to nonlinear state space
fx=[qdd;qd];
hx=[IMU];


text:
 function dt = distance(S,p1,p2)
        dt = norm(S(:,p1) - S(:,p2));
 end

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function ndbi = NDBI( nir, swir )
%NDBI Calculate Normalized Difference Build-up Index (NDBI) using NIR and
%     SWIR bands.
% Syntax
%
%     ndbi = NDBI(swir,nir)
%
% Description
%
%     This function calculates Normalized Difference Build-up Index (NDBI)
%     using SWIR and NIR bands (as following equation).
%
% Input arguments
%
%     swir        Short-wave infrared band
%     nir         Near infrared band
%
% Output arguments
%
%     ndbi        Normalized Difference Build-up Index
%
%
% Author: Shi Qiu (shi.qiu@ttu.edu)
% Date: 19. Dec., 2017

    % calculate NDBI
    ndbi=(swir-nir)./(swir+nir);

    % fix unnormal pixels
    % not 0.01 any more because we will identify urban pixel using ndbi more than 0.
%     ndbi((nir+swir)==0)=0.0;

end

author: GERSL

repo: Fmask

save_path: github-repos/MATLAB/GERSL-Fmask


text:
filename='test2d_micro';
%filename='RVE_Square_Triangle_FineFine';
%filename = 'MicroQuad';
ptype = 'MICRO';
method = 'SIMPALL';
%method = 'SIMP_P3';
materialType = 'ISOTROPIC';
initial_case = 'circleInclusion';
cost={'chomog_alphabeta'};
weights=[1];
constraint = {'volumeConstraint'};
constraint_case = 'EQUALITY';
%incrementFactor = 1;
designVariable = 'Density';
%designVariable = 'LevelSet';
filterType = 'P1';
fracRadius = 0.4;
%optimizer = 'IPOPT';
optimizer = 'DualNestedInPrimal';
%optimizer = 'AlternatingPrimalDual';

optimizerUnconstrained = 'PROJECTED GRADIENT';
line_search_initiator = 'INCREASING LAST STEP';
incrementFactor = 2.95;

%optimizerUnconstrained = 'SLERP';

nsteps = 1;
Vfrac_final = 0.25;
Perimeter_target=1;
optimality_final = 0.2*1e-3;
constr_final =1e-12;

Vfrac_initial = 0.8;
optimality_initial = 0.2*1e-3;
constr_initial = 1e-12;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

%Micro
epsilon_isotropy_initial=1e-1;
epsilon_isotropy_final = 1e-3;
micro.alpha =[1 1 0]';
micro.beta =[1 1 0]';

%nsteps = 10;

% For all tests
plotting = true;
printing = true;
monitoring = true;
monitoring_interval = 1;
maxiter = 2000;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [hd_estimates]=HDestimatesols(hd_record,n,T,HDband,strctident)

% function [hd_estimates]=bear.hdestimates(hd_record,n,T,HDband)
% calculates the point estimate (median), lower bound and upper bound of the historical decomposition from the posterior distribution
% inputs:  - cell 'hd_record': record of the gibbs sampler draws for the historical decomposition
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'T': number of sample time periods (defined p 7 of technical guide)
%          - scalar 'HDband': confidence level for forecasts
% outputs: - cell 'hd_estimates': lower bound, point estimates, and upper bound for the historical decomposition

% create first the cell that will contain the estimates
hd_estimates=cell(length(hd_record),n);

if strctident.MM==0
% deal with shocks in turn
for ii=1:n
   % loop over variables
   for jj=1:length(hd_record)
      % loop over time periods
      for kk=1:T
      % consider the higher and lower confidence band for the hd
      % lower bound
      hd_estimates{jj,ii}(1,kk)=quantile(hd_record{jj,ii}(2,kk),(1-HDband)/2);
      %mean value
      hd_estimates{jj,ii}(2,kk)=quantile(hd_record{jj,ii}(2,kk),0.5);
      % upper bound
      hd_estimates{jj,ii}(3,kk)=quantile(hd_record{jj,ii}(2,kk),HDband+(1-HDband)/2);
      end
   end
end

elseif strctident.MM==1 %Median Model
for ii=1:n
   % loop over variables
   for jj=1:length(hd_record)
      % loop over time periods
      for kk=1:T
      % consider the higher and lower confidence band for the hd
      % lower bound
      hd_estimates{jj,ii}(1,kk)=quantile(hd_record{jj,ii}(:,kk),(1-HDband)/2);
      %medianmodel
      hd_estimates{jj,ii}(2,kk)= hd_record{jj,ii}(medianmodel,kk); %get the best performing model in terms of IRFs
      % upper bound
      hd_estimates{jj,ii}(3,kk)=quantile(hd_record{jj,ii}(:,kk),HDband+(1-HDband)/2);
      % upper bound
      end
   end
end
end

author: european-central-bank

repo: BEAR-toolbox


text:
%% Transient diffusion equation
%% PDE and boundary conditions
% The transient diffusion equation reads
%
% $$\alpha\frac{\partial c}{\partial t}+\nabla.\left(-D\nabla c\right)=0,$$
%
% where $c$ is the independent variable (concentration, temperature, etc)
% , $D$ is the diffusion coefficient, and $\alpha$ is a constant.
% Written by Ali A. Eftekhari
% Last checked: June 2021
clc

%% Define the domain and create a mesh structure
L = 50;  % domain length
Nx = 10; % number of cells
m = createMesh3D(Nx,Nx+3,Nx+6, L/2,2*pi,L);
%% Create the boundary condition structure
BC = createBC(m); % all Neumann boundary condition structure
% BC.left.a(:) = 0; BC.left.b(:)=1; BC.left.c(:)=0; % left boundary
BC.right.a(:) = 0; BC.right.b(:)=1; BC.right.c(:)=0; % right boundary
% BC.top.a(:) = 0; BC.top.b(:)=1; BC.top.c(:)=3; % top boundary
% BC.bottom.a(:) = 0; BC.bottom.b(:)=1; BC.bottom.c(:)=0; % bottom boundary
BC.top.periodic=1;
BC.back.a(:) = 0; BC.back.b(:)=1; BC.back.c(:)=3; % back boundary
BC.front.a(:) = 0; BC.front.b(:)=1; BC.front.c(:)=0; % front boundary
%% define the transfer coeffs
D_val = 1;
D = createCellVariable(m, D_val);
alfa = createCellVariable(m, 1);
u = createFaceVariable(m, [0,0,0.5]);
%% define initial values
c_init = 1;
c_old = createCellVariable(m, c_init,BC); % initial values
c = c_old; % assign the old value of the cells to the current values
%% loop
Dave = harmonicMean(D);
Mdiff = diffusionTerm(Dave);
[Mbc, RHSbc] = boundaryCondition(BC);
FL = fluxLimiter('Superbee');
Mconv = convectionTvdTerm(u, c, FL);
dt = 1; % time step
final_t = 50;
for t=dt:dt:final_t
    [M_trans, RHS_trans] = transientTerm(c_old, dt, alfa);
    M = M_trans-Mdiff+Mbc+Mconv;
    RHS = RHS_trans+RHSbc;
    c = solvePDE(m,M, RHS);
    c_old = c;
    figure(1);visualizeCells(c);drawnow;
end
%% visualization
 figure(1);visualizeCells(c);


text:
function [ Tt ] = updateTransform2(T,w,v,dt)
% This function takes the input of the
normw=norm(w);
if normw>0
    R1=T(1:3,1:3);
    wCrossR1=zeros(3,3);
    for i=1:3
        wCrossR1(:,i)=cross(w,R1(:,i));
    end
    Rt=R1+wCrossR1*dt;
else
    Rt=T(1:3,1:3);
end

dx=v*dt;
xt=T(1:3,4)+dx;

Tt=[Rt,xt];
Tt=[Tt;
    0 0 0 1];
Tt=normalizeTransformMatrix(Tt);
end

function R=quat2rot(q)
R=[(1-2*q(3)*q(3)-2*q(4)*q(4)), 2*(q(1)*q(3)-q(1)*q(4)), 2*(q(2)*q(4)+q(1)*q(3));
2*(q(2)*q(3)+q(1)*q(4)), (1-2*q(2)*q(2)-2*q(4)*q(4)), 2*(q(3)*q(4)-q(1)*q(2));
2*(q(2)*q(4)-q(1)*q(3)), 2*(q(3)*q(4)+q(1)*q(2)), (1-2*q(2)*q(2)-2*q(3)*q(3))];
end

function T=normalizeTransformMatrix(t)
T=zeros(4,4);
T(:,4)=t(:,4);
for i=1:3
    T(1:3,i)=t(1:3,i)/norm(t(1:3,i));
end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function stepfwd(H)
%STEPFWD Step forward one frame of movie in movie player

feval(H.fcns.step_fwd,[],[],H.hfig);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [x]=als_solve_rx(mat, rhs, tol, drx, nswp, addswp)
%Computes an approximate low-rank solution for 2D case
%   [x]=ALS_SOLVE_RX(MAT, RHS, [TOL], [RX], [NSWP])
%   Finds a solution to 2D TTM matrix MAT using the ALS to a 2D TT tensor
%   with rank rx, but the RHS and X are represented as full vectors.
%   TOL is the tolerance for ||x_{i+1}-x_i||/||x_i||,
%   DRX is the random kick rank,
%   NSWP - number of ALS sweeps.
%   default values:
%   tol: 1e-12
%   rx: 1
%   nswp: 10
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

a1 = mat{1}; a2 = mat{2};
n1 = size(a1,1); m1 = size(a1,2);
n2 = size(a2,1); m2 = size(a2,2);
ra = size(a1,3);

% tol2 = 1e-3;

rhs = reshape(rhs, n1, n2);

if (nargin<3)||(isempty(tol))
    tol=1e-12;
end;
if (nargin<4)||(isempty(drx))
    drx=1;
end;
if (nargin<5)||(isempty(nswp))
    nswp=10;
end;
if (nargin<6)||(isempty(addswp))
    addswp=2;
end;

if (drx>m1)||(drx>m2)
    drx = min(m1,m2);
end;

rx=1;

curx = cell(2,1);
curx{1}=rand(m1,rx);
curx{2}=rand(m2,rx);
x = curx{1}*(curx{2}.');
derr = 2;
sp = 0;
resid_old = 1;
for swp=1:nswp

    % QR 2-1
    [q,rv]=qr(curx{2},0); % m2,rx' - rx',rx
    rx = size(q,2);
    curx{2}=q;
%     curx{1}=curx{1}*(rv.');

    % compute phi
    a2 = permute(mat{2}, [1 3 2]);
    a2 = reshape(a2, n2*ra, m2);
    phi = a2*curx{2}; % size n2*ra, rx
    phi = reshape(phi, n2, ra*rx);
    phi = (phi')*phi; % size ra*rx, ra*rx <-- for cplx should also work
    phi = reshape(phi, ra, rx, ra, rx);
    phi = reshape(permute(phi, [1 3 2 4]), ra*ra, rx*rx);
%     phi = reshape(permute(phi, [3 1 2 4]), ra, ra*rx*rx);
%     a2 = reshape(mat{2}, n2, m2*ra);
%     phi = (a2.')*phi; % size m2*ra, ra*rx
%     phi = reshape(phi, m2, ra*ra*rx);
%     phi = (curx{2}.')*phi; % size rx, ra*ra*rx
%     phi = reshape(permute(reshape(phi, rx, ra, ra, rx), [2 3 1 4]), ra*ra, rx*rx);
%     % And the projection of the matrix
    a1 = reshape(permute(mat{1}, [3 2 1]), ra*m1, n1);
    a1 = conj(a1)*reshape(mat{1}, n1, m1*ra); % size ra*m1, m1*ra <-- conjugate!
    a1 = reshape(a1, ra, m1, m1, ra);
    a1 = reshape(permute(a1, [1 4 2 3]), ra*ra, m1*m1);
%     a1 = reshape(permute(mat{1}, [3 2 1]), ra, m1*n1);
    a1 = (phi.')*a1; % size rx*rx, m1*n1
    a1 = reshape(a1, rx, rx, m1, m1);
    a1 = reshape(permute(a1, [2 4 1 3]), rx*m1, rx*m1);
%     a1 = reshape(mat{1}, n1*m1, ra)*phi; % size n1*m1, ra*rx*rx
%     a1 = reshape(a1, n1, m1, ra, rx, rx);
%     a1 = reshape(permute(a1, [1 3 2 4 5]), n1*ra, m1*rx*rx);
%     a1 = conj(reshape(permute(mat{1}, [2 1 3]), m1, n1*ra))*a1; % size m1, m1*rx*rx
%     a1 = reshape(a1, m1, m1, rx, rx);
%     a1 = reshape(permute(a1, [3 1 4 2]), rx*m1, rx*m1);

    %rhs:

    rhs1 = rhs*conj(reshape(mat{2}, n2, m2*ra)); % size n1, m2*ra <-- conjugate
    rhs1 = reshape(rhs1, n1, m2, ra);
    rhs1 = reshape(permute(rhs1, [1 3 2]), n1*ra, m2);

    rhs1 = conj(reshape(permute(mat{1}, [2 1 3]), m1, n1*ra))*rhs1; % size m1, m2
%     rhs1 = rhs;
    rhs1 = rhs1*conj(curx{2}); % size m1,rx
    rhs1 = reshape(rhs1.', rx*m1, 1);

    curx{1}=a1 \ rhs1; % new first block
%     cond_a1 = cond(a1)
    curx{1}=reshape(curx{1}, rx, m1).';

    % Now, let's try the kickass by rank drx:
    if (mod(swp,addswp)==0)
%     if (sp>5)
        curx{1}=[curx{1}, randn(m1,drx)];
%         sp=0;
    end;
%     rx=rx+1;

    % Now, let's compute the second block
    [q,rv]=qr(curx{1},0); % m1,rx' - rx',rx
    rx = size(q,2);
    curx{1}=q;

    % compute phi
    a1 = permute(mat{1}, [1 3 2]);
    a1 = reshape(a1, n1*ra, m1);
    phi = a1*q; % size n1*ra, rx
    phi = reshape(phi, n1, ra*rx);
    phi = (phi')*phi; % size ra*rx, ra*rx
    phi = reshape(phi, ra, rx, ra, rx);
    phi = reshape(permute(phi, [1 3 2 4]), ra*ra, rx*rx);
%     a1 = reshape(mat{1}, n1, m1*ra);
%     phi = (a1.')*phi; % size m1*ra, ra*rx
%     phi = reshape(phi, m1, ra*ra*rx);
%     phi = (curx{1}.')*phi; % size rx, ra*ra*rx
%     phi = reshape(permute(reshape(phi, rx, ra, ra, rx), [2 3 1 4]), ra*ra, rx*rx);
    % And the projection of the matrix
    a2 = reshape(permute(mat{2}, [3 2 1]), ra*m2, n2);
    a2 = conj(a2)*reshape(mat{2}, n2, m2*ra); % size ra*m2, m2*ra
    a2 = reshape(a2, ra, m2, m2, ra);
    a2 = reshape(permute(a2, [1 4 2 3]), ra*ra, m2*m2);
%     a2 = reshape(permute(mat{2}, [3 2 1]), ra, m2*n2);
    a2 = (phi.')*a2; % size rx*rx, m2*n2
    a2 = reshape(a2, rx, rx, m2, m2);
    a2 = reshape(permute(a2, [2 4 1 3]), rx*m2, rx*m2);

    %rhs:
    rhs2 = rhs*conj(reshape(mat{2}, n2, m2*ra)); % size n1, m2*ra
    rhs2 = reshape(rhs2, n1, m2, ra);
    rhs2 = reshape(permute(rhs2, [1 3 2]), n1*ra, m2);

    rhs2 = conj(reshape(permute(mat{1}, [2 1 3]), m1, n1*ra))*rhs2; % size m1, m2
%     rhs2 = rhs;
    rhs2 = (curx{1}')*rhs2; % size rx,m2
    rhs2 = reshape(rhs2, rx*m2, 1);

    curx{2}=a2 \ rhs2; % new first block
    curx{2}=reshape(curx{2}, rx, m2).';

    x_new = curx{1}*(curx{2}.'); % size m1,m2
    derr = norm(x_new(:)-x(:))/norm(x(:));

    x = x_new;

    resid = norm(tt_mat_full_vec(mat, x(:))-rhs(:))/norm(rhs(:));
    conv_fact = resid_old/resid;
    if (conv_fact-1<1e-4)
        sp=sp+1;
    end;

    fprintf('als_solve: swp=%d, derr=%3.3e, rx=%d, resid=%3.3e, conv-1=%3.5e\n', swp, derr, rx, resid, conv_fact-1);
    if (derr<tol)
        break;
    end;
    resid_old = resid;
end;

x = x(:);

end

author: oseledets

repo: TT-Toolbox


text:
function value = ch_roman_to_i4 ( c )

%*****************************************************************************80
%
%% CH_ROMAN_TO_I4 returns the integer value of a single digit of a Roman numeral.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, a Roman numeral.
%
%    Output, integer VALUE, the value of the Roman numeral.
%    If the Roman numeral was not recognized, 0 is returned.
%
  if ( s_eqi ( c, 'M' ) )
    value = 1000;
  elseif ( s_eqi ( c, 'D' ) )
    value = 500;
  elseif ( s_eqi ( c, 'C' ) )
    value = 100;
  elseif ( s_eqi ( c, 'L' ) )
    value = 50;
  elseif ( s_eqi ( c, 'X' ) )
    value = 10;
  elseif ( s_eqi ( c, 'V' ) )
    value = 5;
  elseif ( s_eqi ( c, 'I' ) || s_eqi ( c, 'J' ) )
    value = 1;
  else
    value = 0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%SF_SIMP_P1BUB Linear Lagrange shape function for simplices with bubble (P1+).
%
%   [ VBASE, NLDOF, XLDOF, SFUN ] = SF_SIMP_P1BUB( I_EVAL, N_SDIM, N_VERT, I_DOF, XI, AINVJAC, VBASE )
%   Evaluates conforming linear P1 Lagrange shape functions on simplices
%   an additional with bubble function. XI Barycentric coordinates.
%
%       Input       Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       i_eval      scalar:  1             Evaluate function values
%                           >1             Evaluate values of derivatives
%       n_sdim      scalar: 1-3            Number of space dimensions
%       n_vert      scalar: 2-4            Number of vertices per cell
%       i_dof       scalar: 1-n_ldof       Local basis function to evaluate
%       xi          [n_sdim+1]             Local coordinates of evaluation point
%       aInvJac     [n,n_sdim+1*n_sdim]    Inverse of transformation Jacobian
%       vBase       [n]                    Preallocated output vector
%                                                                                .
%       Output      Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       vBase       [n]                    Evaluated function values
%       nLDof       [4]                    Number of local degrees of freedom on
%                                          vertices, edges, faces, and cell interiors
%       xLDof       [n_sdim,n_ldof]        Local coordinates of local dofs
%       sfun        string                 Function name of called shape function
%
%   See also SF_SIMP_P1

% Copyright 2013-2022 Precise Simulation, Ltd.

% Evaluation type flag.

author: precise-simulation

repo: featool-multiphysics


text:
%DEMO_MODELASSESMENT1 Demonstration for model assessment with WAIC,
% DIC, number of effective parameters and ten-fold cross validation
%
%  Description
%    We will consider the regression problem in demo_regression1.
%    The analysis is conducted with full Gaussian process, and FIC
%    and PIC sparse approximations. The performance of these models
%    are compared by evaluating ten-fold cross validation,
%    leave-one-out cross-validation, WAIC, DIC and the effective
%    number of parameters. The inference will be conducted using
%    maximum a posterior (MAP) estimate for the parameters, via
%    full Markov chain Monte Carlo (MCMC) and with an integration
%    approximation (IA) for the parameters.
%
%    This demo is organised in three parts:
%     1) data analysis with full GP model
%     2) data analysis with FIC approximation
%     3) data analysis with PIC approximation
%
%  See also DEMO_REGRESSION1, DEMO_REGRESSION_SPARSE1
%
% Copyright (c) 2009-2010 Jarno Vanhatalo
% Copyright (c) 2010-2012 Aki Vehtari

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

%========================================================
% PART 1 data analysis with full GP model
%========================================================
disp('Full GP model with Gaussian noise model')

% Load the data
S = which('demo_regression1');
L = strrep(S,'demo_regression1.m','demodata/dat.1');
data=load(L);
x = [data(:,1) data(:,2)];
y = data(:,3);
[n, nin] = size(x);

DIC=repmat(NaN,1,9);DIC2=repmat(NaN,1,9);DIC_latent=repmat(NaN,1,9);
p_eff=repmat(NaN,1,9);p_eff2=repmat(NaN,1,9);p_eff_latent=repmat(NaN,1,9);p_eff_latent2=repmat(NaN,1,9);

% ---------------------------
% --- Construct the model ---
gpcf = gpcf_sexp('lengthScale', [1 1], 'magnSigma2', 0.2^2);
lik = lik_gaussian('sigma2', 0.2^2);
gp = gp_set('lik', lik, 'cf', gpcf, 'jitterSigma2', 1e-9);

% -----------------------------
% --- Conduct the inference ---
%
% We will make the inference first by finding a maximum a posterior estimate
% for the parameters via gradient based optimization. After this we will
% perform an extensive Markov chain Monte Carlo sampling for the parameters.
%
disp('MAP estimate for the parameters')

% --- MAP estimate ---
%     (see gp_optim for more details)

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the BFGS quasi-Newton method
gp=gp_optim(gp,x,y,'opt',opt,'optimf',@fminlbfgs);

% Evaluate the effective number of parameters, DIC and WAIC with focus on
% latent variables.
models{1} = 'full_MAP';
p_eff_latent = gp_peff(gp, x, y);
% For easier comparison to other methods, compute mean log
% predictive density (mlpd) instead of deviance (-2n*mlpd)
[DIC_latent, p_eff_latent2] = gp_dic(gp, x, y);
WAICV(1) = gp_waic(gp,x,y);
WAICG(1) = gp_waic(gp,x,y, 'method', 'G');

% Evaluate the 10-fold cross-validation results.
disp('MAP estimate for the parameters - k-fold-CV')
cvres =  gp_kfcv(gp, x, y, 'display', 'fold');
mlpd_cv(1) = cvres.mlpd_cv;
rmse_cv(1) = cvres.mrmse_cv;

disp('MAP estimate for the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(gp, x, y);
mlpd_loo(1) = mean(lpy);
rmse_loo(1) = sqrt(mean((y-Ef).^2));

% --- MCMC approach ---
disp('MCMC integration over the parameters')

% Do the sampling (this takes about 1 minute)
[rfull,g,opt] = gp_mc(gp, x, y, 'nsamples', 220, 'display', 20);

% After sampling delete the burn-in and thin the sample chain
rfull = thin(rfull, 21, 2);

% Evaluate the effective number of parameters, DIC and WAIC.
models{2} = 'full_MCMC';
% For easier comparison to other methods, compute mean log
% predictive density (mlpd) instead of deviance (-2n*mlpd)
[DIC(2), p_eff(2)] =  gp_dic(rfull, x, y, 'focus', 'hyper');
[DIC2(2), p_eff2(2)] =  gp_dic(rfull, x, y);
WAICV(2) = gp_waic(rfull,x,y);
WAICG(2) = gp_waic(rfull,x,y, 'method', 'G');

% Evaluate the 10-fold cross validation results.
%
% We reduce the number of samples so that the sampling takes less time.
% 50 is too small sample size, though, and for reliable results the 10-CV
% should be run with larger sample size.
disp('MCMC integration over the parameters - k-fold-CV')
opt.nsamples= 50;
cvres =  gp_kfcv(gp, x, y, 'inf_method', 'MCMC', 'opt', opt, 'rstream', 1, 'display', 'fold');
mlpd_cv(2) = cvres.mlpd_cv;
rmse_cv(2) = cvres.mrmse_cv;

disp('MCMC integration over the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(rfull, x, y);
mlpd_loo(2) = mean(lpy);
rmse_loo(2) = sqrt(mean((y-Ef).^2));

% --- Integration approximation approach ---
disp('Grid integration over the parameters')

gp_array = gp_ia(gp, x, y, 'int_method', 'grid');

models{3} = 'full_IA';
% For easier comparison to other methods, compute mean log
% predictive density (mlpd) instead of deviance (-2n*mlpd)
[DIC(3), p_eff(3)] =  gp_dic(gp_array, x, y, 'focus', 'hyper');
[DIC2(3), p_eff2(3)] =  gp_dic(gp_array, x, y);
WAICV(3) = gp_waic(gp_array,x,y);
WAICG(3) = gp_waic(gp_array,x,y, 'method', 'G');

% Then the 10 fold cross-validation.
disp('Grid integration over the parameters - k-fold-CV')
clear opt
opt.int_method = 'grid';
cvres = gp_kfcv(gp, x, y, 'inf_method', 'IA', 'opt', opt, 'display', 'fold');
mlpd_cv(3) = cvres.mlpd_cv;
rmse_cv(3) = cvres.mrmse_cv;

disp('Grid integration over the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(gp_array, x, y);
mlpd_loo(3) = mean(lpy);
rmse_loo(3) = sqrt(mean((y-Ef).^2));

%========================================================
% PART 2 data analysis with FIC GP
%========================================================
disp('GP with FIC sparse approximation')

% ---------------------------
% --- Construct the model ---

% Here we conduct the same analysis as in part 1, but this time we
% use FIC approximation

% Initialize the inducing inputs in a regular grid over the input space
[u1,u2]=meshgrid(linspace(-1.8,1.8,6),linspace(-1.8,1.8,6));
X_u = [u1(:) u2(:)];

% Create the FIC GP structure
gp_fic = gp_set('type', 'FIC', 'lik', lik, 'cf', gpcf, 'jitterSigma2', 1e-6, 'X_u', X_u)

% -----------------------------
% --- Conduct the inference ---

% --- MAP estimate using scaled conjugate gradient algorithm ---
disp('MAP estimate for the parameters')

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the BFGS quasi-Newton method
gp_fic=gp_optim(gp_fic,x,y,'opt',opt,'optimf',@fminlbfgs);

% Evaluate the effective number of parameters and DIC with focus on
% latent variables.
models{4} = 'FIC_MAP';
p_eff_latent(4) = gp_peff(gp_fic, x, y);
[DIC_latent(4), p_eff_latent2(4)] = gp_dic(gp_fic, x, y);
WAICV(4) = gp_waic(gp_fic,x,y);
WAICG(4) = gp_waic(gp_fic,x,y, 'method', 'G');

% Evaluate the 10-fold cross validation results.
disp('MAP estimate for the parameters - k-fold-CV')
cvres = gp_kfcv(gp_fic, x, y, 'display', 'fold');
mlpd_cv(4) = cvres.mlpd_cv;
rmse_cv(4) = cvres.mrmse_cv;

disp('MAP estimate for the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(gp_fic, x, y);
mlpd_loo(4) = mean(lpy);
rmse_loo(4) = sqrt(mean((y-Ef).^2));

% --- MCMC approach ---
% (the inducing inputs are fixed)
disp('MCMC integration over the parameters')

% Do the sampling (this takes about 1 minute)
rfic = gp_mc(gp_fic, x, y, 'nsamples', 220, 'display', 20);

% After sampling we delete the burn-in and thin the sample chain
rfic = thin(rfic, 21, 2);

% Evaluate the effective number of parameters, DIC and WAIC. Note that
% the effective number of parameters as a second output, but here
% we use explicitly the gp_peff function
models{5} = 'FIC_MCMC';
[DIC(5), p_eff(5)] =  gp_dic(rfic, x, y, 'focus', 'hyper');
[DIC2(5), p_eff2(5)] =  gp_dic(rfic, x, y);
WAICV(5) = gp_waic(rfic,x,y);
WAICG(5) = gp_waic(rfic,x,y, 'method', 'G');

% We reduce the number of samples so that the sampling takes less time.
% 50 is too small sample size, though, and for reliable results the 10-CV
% should be run with larger sample size. We also set the save option to 0.
clear opt
opt.nsamples= 50; opt.display=20;
disp('MCMC integration over the parameters - k-fold-CV')
cvres = gp_kfcv(gp_fic, x, y, 'inf_method', 'MCMC', 'opt', opt, 'display', 'fold');
mlpd_cv(5) = cvres.mlpd_cv;
rmse_cv(5) = cvres.mrmse_cv;

disp('MCMC integration over the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(rfic, x, y);
mlpd_loo(5) = mean(lpy);
rmse_loo(5) = sqrt(mean((y-Ef).^2));

% --- Integration approximation approach ---
disp('Grid integration over the parameters')
gpfic_array = gp_ia(gp_fic, x, y, 'int_method', 'grid');

models{6} = 'FIC_IA';
[DIC(6), p_eff(6)] =  gp_dic(gpfic_array, x, y, 'focus', 'hyper');
[DIC2(6), p_eff2(6)] =  gp_dic(gpfic_array, x, y);
WAICV(6) = gp_waic(gpfic_array,x,y);
WAICG(6) = gp_waic(gpfic_array,x,y, 'method', 'G');

% Then the 10 fold cross-validation.
disp('Grid integration over the parameters - k-fold-CV')
clear opt
opt.int_method = 'grid';
cvres = gp_kfcv(gp_fic, x, y, 'inf_method', 'IA', 'opt', opt, 'display', 'fold');
mlpd_cv(6) = cvres.mlpd_cv;
rmse_cv(6) = cvres.mrmse_cv;

disp('Grid integration over the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(gpfic_array, x, y);
mlpd_loo(6) = mean(lpy);
rmse_loo(6) = sqrt(mean((y-Ef).^2));

%========================================================
% PART 3 data analysis with PIC approximation
%========================================================
disp('GP with PIC sparse approximation')

[u1,u2]=meshgrid(linspace(-1.8,1.8,6),linspace(-1.8,1.8,6));
X_u = [u1(:) u2(:)];

% Initialize test points
[p1,p2]=meshgrid(-1.8:0.1:1.8,-1.8:0.1:1.8);
p=[p1(:) p2(:)];

% set the data points into clusters. Here we construct two cell arrays.
%  trindex  contains the block index vectors for training data. That is
%           x(trindex{i},:) and y(trindex{i},:) belong to the i'th block.
b1 = [-1.7 -0.8 0.1 1 1.9];
mask = zeros(size(x,1),size(x,1));
trindex={};
for i1=1:4
  for i2=1:4
    ind = 1:size(x,1);
    ind = ind(: , b1(i1)<=x(ind',1) & x(ind',1) < b1(i1+1));
    ind = ind(: , b1(i2)<=x(ind',2) & x(ind',2) < b1(i2+1));
    trindex{4*(i1-1)+i2} = ind';
  end
end

% Create the PIC GP structure and set the inducing inputs and block indexes
gpcf = gpcf_sexp('lengthScale', [1 1], 'magnSigma2', 0.2^2);
lik = lik_gaussian('sigma2', 0.2^2);

gp_pic = gp_set('type', 'PIC', 'lik', lik, 'cf', gpcf, 'jitterSigma2', 1e-6, 'X_u', X_u);
gp_pic = gp_set(gp_pic, 'tr_index', trindex);

% -----------------------------
% --- Conduct the inference ---

% --- MAP estimate using scaled conjugate gradient algorithm ---
disp('MAP estimate for the parameters')

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the BFGS quasi-Newton method
gp_pic=gp_optim(gp_pic,x,y,'opt',opt,'optimf',@fminlbfgs);

models{7} = 'PIC_MAP';
p_eff_latent(7) = gp_peff(gp_pic, x, y);
[DIC_latent(7), p_eff_latent2(7)] = gp_dic(gp_pic, x, y);
WAICV(7) = gp_waic(gp_pic, x, y);
WAICG(7) = gp_waic(gp_pic, x, y, 'method', 'G');

% Evaluate the 10-fold cross validation results.
disp('MAP estimate for the parameters - k-fold-CV')
cvres = gp_kfcv(gp_pic, x, y, 'display', 'fold');
mlpd_cv(7) = cvres.mlpd_cv;
rmse_cv(7) = cvres.mrmse_cv;

disp('MAP estimate for the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(gp_pic, x, y);
mlpd_loo(7) = mean(lpy);
rmse_loo(7) = sqrt(mean((y-Ef).^2));

% --- MCMC approach ---
disp('MCMC integration over the parameters')

% Do the sampling (this takes about 1 minute)
rpic = gp_mc(gp_pic, x, y, 'nsamples', 220, 'display', 20);

% After sampling we delete the burn-in and thin the sample chain
rpic = rmfield(rpic, 'tr_index');
rpic = thin(rpic, 21, 2);
rpic.tr_index = trindex;

% Evaluate the effective number of parameters and DIC. Note that
% the effective number of parameters as a second output, but here
% we use explicitly the gp_peff function
models{8} = 'PIC_MCMC';
[DIC(8), p_eff(8)] =  gp_dic(rpic, x, y, 'focus', 'hyper');
[DIC2(8), p_eff2(8)] =  gp_dic(rpic, x, y);
WAICV(8) = gp_waic(rpic, x, y);
WAICG(8) = gp_waic(rpic, x, y, 'method', 'G');

% We reduce the number of samples so that the sampling takes less time.
% 50 is too small sample size, though, and for reliable results the 10-CV
% should be run with larger sample size. We also set the save option to 0.
clear opt
opt.nsamples= 50; opt.display=20;
disp('MCMC integration over the parameters - k-fold-CV')
cvres = gp_kfcv(gp_pic, x, y, 'inf_method', 'MCMC', 'opt', opt, 'display', 'fold');
mlpd_cv(8) = cvres.mlpd_cv;
rmse_cv(8) = cvres.mrmse_cv;

disp('MCMC integration over the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(rpic, x, y);
mlpd_loo(8) = mean(lpy);
rmse_loo(8) = sqrt(mean((y-Ef).^2));

% --- Integration approximation approach ---
disp('Grid integration over the parameters')

gppic_array = gp_ia(gp_pic, x, y, 'int_method', 'grid');

models{9} = 'PIC_IA';
[DIC(9), p_eff(9)] =  gp_dic(gppic_array, x, y, 'focus', 'hyper');
[DIC2(9), p_eff2(9)] =  gp_dic(gppic_array, x, y);
WAICV(9) = gp_waic(gppic_array, x, y);
WAICG(9) = gp_waic(gppic_array, x, y, 'method', 'G');

% Then the 10 fold cross-validation.
disp('Grid integration over the parameters - k-fold-CV')
clear opt
opt.int_method = 'grid';
cvres = gp_kfcv(gp_pic, x, y, 'inf_method', 'IA', 'opt', opt, 'display', 'fold');
mlpd_cv(9) = cvres.mlpd_cv;
rmse_cv(9) = cvres.mrmse_cv;

disp('Grid integration over the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(gppic_array, x, y);
mlpd_loo(9) = mean(lpy);
rmse_loo(9) = sqrt(mean((y-Ef).^2));

%========================================================
% PART 4 Print the results
%========================================================
disp('Summary of the results')
S = '      ';
for i = 1:length(models)
    S = [S '  ' models{i}];
end

S = sprintf([S '\n CV-mlpd  %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], mlpd_cv);
S = sprintf([S '\n CV-rmse  %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], rmse_cv);
S = sprintf([S '\n LOO-mlpd %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], mlpd_loo);
S = sprintf([S '\n LOO-rmse %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], rmse_loo);
S = sprintf([S '\n ']);
S = sprintf([S '\n WAIC_V   %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], WAICV);
S = sprintf([S '\n WAIC_G   %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], WAICG);
S = sprintf([S '\n ']);
S = sprintf([S '\n DIC_h    %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], DIC);
S = sprintf([S '\n DIC_a    %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], DIC2);
S = sprintf([S '\n DIC_l    %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], DIC_latent);
S = sprintf([S '\n peff_h   %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], p_eff);
S = sprintf([S '\n peff_a   %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], p_eff2);
S = sprintf([S '\n peff_l   %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], p_eff_latent);
S = sprintf([S '\n peff_l2  %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], p_eff_latent2);
S = sprintf([S '\n ']);
S = sprintf([S '\n ']);
S = sprintf([S '\n The notation is as follows:']);
S = sprintf([S '\n CV-mlpd  = mean log predictive density from the 10-fold CV. ']);
S = sprintf([S '\n CV-rmse  = root mean squared error from the 10-fold LOO-CV. ']);
S = sprintf([S '\n LOO-mlpd = mean log predictive density from the LOO-CV. ']);
S = sprintf([S '\n LOO-rmse = root mean squared error from the 10-fold CV. ']);
S = sprintf([S '\n WAIC_V   = WAIC via variance method ']);
S = sprintf([S '\n WAIC_G   = WAIC via Gibbs training utility method ']);
S = sprintf([S '\n DIC_h    = DIC with focus on parameters. ']);
S = sprintf([S '\n DIC_a    = DIC with focus on parameters and latent variables (all). ']);
S = sprintf([S '\n DIC_l    = DIC with focus on latent variables. ']);
S = sprintf([S '\n peff_h   = effective number of parameters (latent variables marginalized). ']);
S = sprintf([S '\n peff_a   = effective number of parameters and latent variables. ']);
S = sprintf([S '\n peff_l   = effective number of latent variables evaluated with gp_peff. ']);
S = sprintf([S '\n peff_l2  = effective number of latent variables evaluated with gp_dic. ']);
S = sprintf([S '\n '])

author: gpstuff-dev

repo: gpstuff

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


text:
% fieldtripchan2eeglab() - convert Fieldtrip channel location structure
%                          to EEGLAB channel location structure
%
% Usage:
%   >> chanlocs = fieldtripchan2eeglab( fieldlocs );
%
% Inputs:
%   fieldlocs - Fieldtrip channel structure. See help readlocs()
%
% Outputs:
%   chanlocs  - EEGLAB channel location structure.
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, 2006-
%
% See also: readlocs()

% Copyright (C) 2003 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function chanlocs = fieldtripchan2eeglab( loc );

    if nargin < 1
        help fieldtripchan2eeglab;
        return;
    end;

    chanlocs = struct('labels', loc.label(:)', 'X', mattocell(loc.pnt(:,1)'), ...
                                            'Y', mattocell(loc.pnt(:,2)'), ...
                                            'Z', mattocell(loc.pnt(:,3)'));
    chanlocs = convertlocs(chanlocs, 'cart2all');

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ferns,hsPr] = fernsClfChangeNFerns( data, hs, ferns, Mnew, varargin )
% Train random fern classifier.
%
% See "Fast Keypoint Recognition in Ten Lines of Code" by Mustafa Ozuysal,
% Pascal Fua and Vincent Lepetit, CVPR07.
%
% Dimensions:
%  M - number ferns
%  S - fern depth
%  F - number features
%  N - number input vectors
%  H - number classes
%
% USAGE
%  [ferns,hsPr] = fernsClfTrain( data, hs, [varargin] )
%
% INPUTS
%  data     - [NxF] N length F feature vectors
%  hs       - [Nx1] target output labels in [1,H]
%  varargin - additional params (struct or name/value pairs)
%   .S        - [10] fern depth (ferns are exponential in S)
%   .M        - [50] number of ferns to train
%   .thrr     - [0 1] range for randomly generated thresholds
%   .bayes    - [1] if true combine probs using bayes assumption
%   .ferns    - [] if given reuse previous ferns (recompute pFern)
%
% OUTPUTS
%  ferns    - learned fern model w the following fields
%   .fids     - [MxS] feature ids for each fern for each depth
%   .thrs     - [MxS] threshold corresponding to each fid
%   .pFern    - [2^SxHxM] learned log probs at fern leaves
%   .bayes    - if true combine probs using bayes assumption
%   .inds     - [NxM] cached indices for original training data
%   .H        - number classes
%  hsPr     - [Nx1] predicted output labels
%
% EXAMPLE
%  N=5000; H=5; d=2; [xs0,hs0,xs1,hs1]=demoGenData(N,N,H,d,1,1);
%  fernPrm=struct('S',4,'M',50,'thrr',[-1 1],'bayes',1);
%  tic, [ferns,hsPr0]=fernsClfTrain(xs0,hs0,fernPrm); toc
%  tic, hsPr1 = fernsClfApply( xs1, ferns ); toc
%  e0=mean(hsPr0~=hs0); e1=mean(hsPr1~=hs1);
%  fprintf('errors trn=%f tst=%f\n',e0,e1); figure(1);
%  subplot(2,2,1); visualizeData(xs0,2,hs0);
%  subplot(2,2,2); visualizeData(xs0,2,hsPr0);
%  subplot(2,2,3); visualizeData(xs1,2,hs1);
%  subplot(2,2,4); visualizeData(xs1,2,hsPr1);
%
% See also fernsClfApply, fernsInds
%
% Piotr's Image&Video Toolbox      Version 2.50
% Copyright 2010 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Lesser GPL [see external/lgpl.txt]

dfs={'thrr',[0 1],'tmp',''};
[thrr,~]=getPrmDflt(varargin,dfs,1);

[Mold,S] = size(ferns.fids);
[N,F]=size(data); assert(length(hs)==N);
H=max(hs); assert(all(hs>0)); assert(S<=20);

if Mnew == Mold,
elseif Mnew < Mold,
  Mremove = Mold - Mnew;
  % randomly choose some ferns to remove
  idx_remove = randsample(Mold,Mremove);
  ferns.fids(idx_remove,:) = [];
  ferns.thrs(idx_remove,:) = [];
  ferns.counts(:,:,idx_remove) = [];
  ferns.pFern(:,:,idx_remove) = [];
  ferns.inds(:,idx_remove) = [];
else

  Madd = Mnew - Mold;
  % create some new ferns
  thrs_add=rand(Madd,S)*(thrr(2)-thrr(1))+thrr(1);
  fids_add=uint32(floor(rand(Madd,S)*F+1));
  inds_add=fernsInds(data,fids_add,thrs_add);

  % store new ferns
  ferns.fids(Mold+1:Mnew,:) = fids_add;
  ferns.thrs(Mold+1:Mnew,:) = thrs_add;
  ferns.inds(:,Mold+1:Mnew) = inds_add;

  % get counts for each leaf for each class for each new fern
  pFern_add = nan(2^S,H,Madd);
  edges = 1:2^S;
  for m = 1:Madd,
    for h = 1:H,
      pFern_add(:,h,m) = histc(inds_add(hs==h,m),edges);
    end
  end
  pFern_add = pFern_add + ferns.bayes;
  ferns.counts(:,:,Mold+1:Mnew) = pFern_add;

  % convert fern leaf class counts into probabilities
  if( ferns.bayes<=0 )
    norm = 1./sum(pFern_add,2);
    pFern_add = bsxfun(@times,pFern_add,norm);
  else
    norm = 1./sum(pFern_add,1);
    pFern_add = bsxfun(@times,pFern_add,norm);
    pFern_add=log(pFern_add);
  end
  ferns.pFern(:,:,Mold+1:Mnew) = pFern_add;

  clear pFern;
end

if(nargout==2),
  hsPr=fernsClfApply([],ferns,ferns.inds);
end

end

author: kristinbranson

repo: JAABA


text:
function excelrecord6fcn(identified, IRFperiods, endo, fevd_estimates, n, pref)
% script excelrecord6
% records the information contained in the worksheet 'FEVD' of the excel spreadsheet 'results.xls'

% create the cell that will be saved on excel
fevdcell={};

% build preliminary elements: space between the tables
horzspace=repmat({''},2,5*identified);
vertspace=repmat({''},IRFperiods+3,1);

% loop over variables (vertical dimension)
for ii=1:n
    tempcell={};
    % loop over shocks (horizontal dimension)
    for jj=1:identified
        % create cell of fevd record for variable ii in response to shock jj
        temp=['part of ' endo{ii,1} ' fluctuation due to ' endo{jj,1} ' shocks'];
        fevd_ij=[temp {''} {''} {''};{''} {''} {''} {''};{''} {'lw. bound'} {'median'} {'up. bound'};num2cell((1:IRFperiods)') num2cell((fevd_estimates{ii,jj})')];
        tempcell=[tempcell fevd_ij vertspace];
    end
    fevdcell=[fevdcell;horzspace;tempcell];
end

% trim
fevdcell=fevdcell(3:end,1:end-1);

% write in excel
bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),fevdcell,'FEVD','B2');


text:
function [list,dif]=locDepFun(Fun)
%LOCDEPFUN determines the local m-files that a function depends upon
%    [LIST,DIF]=LOCDEPFUN(FUN)
%    FUN is a string with the function name
%    LIST is a cell array of strings with the local m-files that FUN
%    depends upon
%    DIF is a cell array of strings with the local m-files that FUN does
%    not depend upon
%
%    Example
%        [list,dif]=locDepFun('foo');
%    Determines which local m-files foo.m ooes and does not depend upon
%
%    See also DEPFUN.

%R.F. Tap, May 2006

list = depfun(Fun,'-quiet');
d = cd;
lcd = length(d);
I = [];
for i=1:length(list)
    if length(list{i})>=lcd
        if strcmp(d,list{i}(1:lcd))
            I=[I i];
        end
    end
end

dm=dir('*.m');
names = {dm.name};

ieq = [];
for i=1:length(names)
    LN = length(names{i});
    for j=1:length(I)
        if strcmp(names{i},list{I(j)}(end-LN+1:end))
            ieq = [ieq i];
            break
        end
    end
end

idif = setdiff(1:length(names),ieq);
dif = names(idif)';
list = names(ieq)';

author: Sable

repo: mcbench-benchmarks


text:
function intVal=monomialIntCrossPolytope(alpha)
%%MONOMIALINTCROSSPOLYTOPE Evaluate the integral of
%          prod_{i=1}^N x(i)^alpha(i) over the space defined by
%          sum(abs(x))=1. That is, this function gives the value of the
%          desired monomial integral taken over an N-dimensional cross
%          polytope. The region is designated as G_n in [1]. In two
%          dimensions, this shape is a square, in three an octohedron, in
%          N, a cross polytope (or N-dimensional hyperoctahedron).
%
%INPUTS: alpha An NX1 or 1XN vector of the integer exponents of the
%              monomial term. All elements must be >=0.
%
%OUTPUTS: intVal The value of the specified integral.
%
%The formula is taken from Chapter 7.7 of [1]. These types of explicit
%moment formulae are useful for testing cubature integration points.
%
%EXAMPLE:
%Here we verify that the moment value produced here equals that obtained
%using cubature points of an appropriate order.
% [xi,w]=fifthOrderCrossPolyCubPoints(3);
% alpha=[0;2;2];
% theMoment=findMomentFromSamp(alpha,xi,w)
% intVal=monomialIntCrossPolytope(alpha)
%One will find that theMoment and intVal are both near 0.0063.
%
%REFERENCES:
%[1] A.H. Stroud, Approximate Calculation of Multiple Integrals. Cliffs,
%    NJ: Prentice-Hall, Inc., 1971.
%
%February 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(all(mod(alpha,2)==0))
    n=length(alpha);

    intVal=exp(n*log(2)+sum(gammaln(alpha+1))-gammaln(n+sum(alpha)+1));
else
    intVal=0;
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 h = circle_segment_height_from_area ( r, area )

%*****************************************************************************80
%
%% CIRCLE_SEGMENT_HEIGHT_FROM_AREA: height of a circle segment from area.
%
%  Discussion:
%
%    Begin with a circle of radius R.  Choose two points P1 and P2 on the
%    circle, and draw the chord P1:P2.  This chord divides the circle
%    into two pieces, each of which is called a circle segment.
%    Consider one of the pieces.  The "angle" of this segment is the angle
%    P1:C:P2, where C is the center of the circle.  Let Q be the point on
%    the chord P1:P2 which is closest to C.  The "height" of the segment
%    is the distance from Q to the perimeter of the circle.
%
%    This function is given the radius R and area of the segment, and
%    determines the corresponding height.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%    0 < R.
%
%    Input, real AREA, the area of the circle segment.
%    0 <= AREA <= 2.0 * PI * R^2.
%
%    Output, real H, the height of the circle segment.
%
  if ( area < 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CIRCLE_SEGMENT_HEIGHT_FROM_AREA - Fatal error!\n' );
    fprintf ( 1, '  AREA < 0.0.\n' );
    error ( 'CIRCLE_SEGMENT_HEIGHT_FROM_AREA - Fatal error!' );
  end

  area_circle = 2.0 * pi * r^2;

  if ( area == 0.0 )
    h = 0.0;
    return
  end

  if ( area == area_circle )
    h = 2.0 * r;
    return
  end

  if ( area_circle < area )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CIRCLE_SEGMENT_HEIGHT_FROM_AREA - Fatal error!\n' );
    fprintf ( 1, '  2.0 * pi * r^2 < AREA.\n' );
    error ( 'CIRCLE_SEGMENT_HEIGHT_FROM_AREA - Fatal error!' );
  end

  if ( r <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CIRCLE_SEGMENT_HEIGHT_FROM_AREA - Fatal error!\n' );
    fprintf ( 1, '  R <= 0.0.\n' );
    error ( 'CIRCLE_SEGMENT_HEIGHT_FROM_AREA - Fatal error!' );
  end

  h1 = 0.0;
  a1 = circle_segment_area_from_height ( r, h1 );
  h2 = 2.0 * r;
  a2 = circle_segment_area_from_height ( r, h2 );

  it = 0;

  while ( it < 30 )

    h = 0.5 * ( h1 + h2 );
    a = circle_segment_area_from_height ( r, h );
    it = it + 1;

    if ( abs ( a - area ) < sqrt ( eps ) * area )
      break
    end

    if ( a < area )
      h1 = h;
      a1 = a;
    else
      h2 = h;
      a2 = a;
    end

  end

  return
end


text:
function x = cellrowassign(x, y, indx)

for k = 1:numel(x)
  x{k}(indx,:) = y{k};
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [MPa] = ftH2O2MPa(ftH2O)
% Convert pressure from feet of water column at 4 degrees to megapascals
% Chad Greene 2012
MPa = ftH2O*0.00298907;


text:
function coeffs = alias(coeffs, m)
%ALIAS   Alias Fourier coefficients on equally spaced grid.
%   ALIAS(C, M) aliases the Fourier coefficients stored in the column vector C
%   to have length M. If M > LENGTH(C), the coefficients are padded with zeros.
%   If C is a matrix of coefficients, each of the columns is aliased to length
%   M.
%
% See also PROLONG.

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

% Get the number of coefficients.
n = size(coeffs, 1);

% Pad with zeros:
if ( m > n )

    k = ceil((m-n)/2);
    z = zeros(k, size(coeffs, 2));

    % Need to handle the odd vs. even case separately.
    if ( ~mod(n, 2) ) % n even.

        % First, account for the asymmetry in the coefficients when n is even.
        % This will account for the cos(N/2) coefficient, which is stored
        % in the coeffs(n,:) entry, using properties of the complex
        % exponential.
        coeffs = [ coeffs(1,:)/2; coeffs(2:n,:); coeffs(1,:)/2 ];
        coeffs = [ z; coeffs; z(1:end-1,:) ];

        % Next, check if m is odd. If it is, then coeffs is too long and we
        % need to remove the first row (the lowest degree
        % coefficients).
        if ( mod(m, 2) ) % m odd.
            coeffs = coeffs(2:end,:);
        end

    else % n odd.

        % There is no asymmetry in the coefficients, just pad them.
        coeffs = [ z; coeffs; z ];

        % Only need to check if m is even, in which case coeffs is too
        % long and we need to remove the last row (the highest degree
        % coefficients).
        if ( ~mod(m, 2) ) % m even.
            coeffs = coeffs(1:end-1,:);
        end
    end

    return

end

% If the number of coefficients is even then we extend them by one so they
% are odd by exploiting the symmetry property.  This makes the code below a
% little cleaner since fewer cases need to be handled.
if ( ~mod(n, 2) ) % n even.
    coeffs(1,:) = 0.5*coeffs(1,:);
    coeffs = [ coeffs; coeffs(1,:) ];
    n = n + 1;
end

% Need to handle the odd and even case (for m) differently.
if ( mod(m, 2) ) % m is odd.

    if ( m == 1 )

        % Reduce to a single point:
        constCoeffs = coeffs((n-1)/2+1,:);
        posCoeffs = coeffs((n-1)/2:-1:1,:);
        negCoeffs = coeffs((n-1)/2+2:n,:);
        e = ones(1, ceil((n-1)/2));
        e(1:2:end) = -1;
        coeffs = constCoeffs + (e*posCoeffs + e*negCoeffs);

    else

        m2 = (m-1)/2;
        n2 = (n-1)/2;
        % Extract coefficients:
        aliasedCoeffs = coeffs(n2-m2+1:n2+m2+1,:);

        % The code below aliases the coefficients from the higher modes
        % onto the lower modes. The principle behind the formula is figure
        % out which of the higher Fourier modes are indistinguishable from
        % the lower Fouirer modes on the grid consisting of m equally
        % spaced points from [-1,1). In general, when m and n are odd the
        % following will be equal for j=-(n-1)/2 to -(m+1)/2
        % exp(1i*pi*j*x) = sgn*exp(1i*pi*k*x) where
        % k = mod(j+(m+1)/2, -m) + (m-1)/2 and sgn = (-1)^mod(j+k, 2);
        for j = -n2:-m2-1
            k = mod(j+m2+1, -m) + m2;
            coeffIndexK = k + m2 + 1;  % Index into aliasedCoeffs for mode k.
            coeffIndexJ = j + n2 + 1;  % Index into coeffs for mode j.
            sgn = (-1)^mod(j+k, 2);
            aliasedCoeffs(coeffIndexK,:) = aliasedCoeffs(coeffIndexK,:) + sgn*coeffs(coeffIndexJ,:);
            coeffIndexK = -k + m2 + 1;
            coeffIndexJ = -j + n2 + 1;
            aliasedCoeffs(coeffIndexK,:) = aliasedCoeffs(coeffIndexK,:) + sgn*coeffs(coeffIndexJ,:);
        end
        coeffs = aliasedCoeffs;

    end

else % m is even.

    m2 = m/2;
    n2 = (n-1)/2;

    % Put the coefficient for the cos(m/2*pi*x) in the first entry of the
    % coefficient vector.  This corresopnds to the exp(-1i*pi*m/2*x).
    aliasedCoeffs = coeffs(n2+1-m2:n2+m2,:);
    % Extend the aliased cofficient vector so it is symmetric.  This allows
    % us to easily account for aliasing on the exp(-1i*pi*m/2*x) and
    % exp(1i*pi*m/2*x) term in the code below without the need for a
    % special if check.  The negative allows contribution for sin(m/2*pi*x)
    % to be removed.
    aliasedCoeffs = [ aliasedCoeffs; -aliasedCoeffs(1,:) ];

    % Follow a similar structure to the odd m case above.  The main
    % difference is that the higher order modes do not change sign when
    % aliased onto an even point grid.
    for j = -n2:-m2
        k = mod(j+m2,-m) + m2;
        coeffIndexK = k + m2 + 1;  % Index into aliasedCoeffs for mode k.
        coeffIndexJ = j + n2 + 1;  % Index into coeffs for mode j.
        aliasedCoeffs(coeffIndexK,:) = aliasedCoeffs(coeffIndexK,:) + coeffs(coeffIndexJ,:);
        coeffIndexK = -k + m2 + 1;
        coeffIndexJ = -j + n2 + 1;
        aliasedCoeffs(coeffIndexK,:) = aliasedCoeffs(coeffIndexK,:) + coeffs(coeffIndexJ,:);
    end
    % Collapse the aliased coefficient vector down to an even number of
    % terms by adding in the aliasing of the exp(1i*pi*m/2*x) terms to the
    % exp(-1i*pi*m/2*x).
    aliasedCoeffs(1,:) = aliasedCoeffs(1,:) + aliasedCoeffs(end,:);
    aliasedCoeffs(end,:) = [];

    coeffs = aliasedCoeffs;

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun