text:
function S = globMatrixNed3DStiff(fun,mesh,fem1,fem2)

%% USAGE: generate stiffness global matrix on a 3D mesh
%
% INPUTS:
% fun --- coefficient function
% mesh --- a struct data contains very rich mesh information.
% fem1 --- global DoF for test function space
% fem2 --- global DoF for trial function space
%
% OUTPUTS:
% [IN JN XN] --- triplets of the sparse matrix from regular elements.
% [II JI XI] --- triplets of the sparse matrix from interface elements.

% Last Modified: 08/07/2020 by Xu Zhang

%% 0. Initializaiton
feEvalBas1 = @EvalNed1Bas3D;
feEvalBas2 = @EvalNed1Bas3D;

dof1 = fem1.ldof; dof2 = fem2.ldof; nloc = dof1*dof2;
nt = length(mesh.t);
A = fem1.area; gx = fem1.gx; gy = fem1.gy; gz = fem1.gz; gw = fem1.gw;
X = zeros(nloc*nt, 1);

coef = feval(fun,gx,gy,gz);
Ibasx = cell(dof1,1); Ibasy = cell(dof1,1); Ibasz = cell(dof1,1);
Jbasx = cell(dof2,1); Jbasy = cell(dof2,1); Jbasz = cell(dof2,1);
for i = 1:dof1
    Ibasx{i} = feEvalBas1(fem1.bas, ':', gx, gy, gz, i, 1, 1).*fem1.t_e_orit(:,i);
    Ibasy{i} = feEvalBas1(fem1.bas, ':', gx, gy, gz, i, 1, 2).*fem1.t_e_orit(:,i);
    Ibasz{i} = feEvalBas1(fem1.bas, ':', gx, gy, gz, i, 1, 3).*fem1.t_e_orit(:,i);
end
for j = 1:dof2
    Jbasx{j} = feEvalBas2(fem2.bas, ':', gx, gy, gz, j, 1, 1).*fem2.t_e_orit(:,j);
    Jbasy{j} = feEvalBas2(fem2.bas, ':', gx, gy, gz, j, 1, 2).*fem2.t_e_orit(:,j);
    Jbasz{j} = feEvalBas2(fem2.bas, ':', gx, gy, gz, j, 1, 3).*fem2.t_e_orit(:,j);
end

I = reshape(repmat(fem1.g2ldof,6,1),nloc*nt,1);
J = repmat(reshape(fem2.g2ldof,dof2*nt,1),6,1);
ind = 0;
for i = 1:dof1
    for j = 1:dof2
        X(ind+1:ind+nt) = A.*(sum(((Ibasx{i}.*(coef.*Jbasx{j})).*gw'),2) + ...
            sum(((Ibasy{i}.*(coef.*Jbasy{j})).*gw'),2) + ...
            sum(((Ibasz{i}.*(coef.*Jbasz{j})).*gw'),2));
        ind = ind + nt;
    end
end
ID = find(X~=0);
S = sparse(I(ID),J(ID),X(ID),size(fem1.gdof,1),size(fem2.gdof,1));

repo: ifem


text:
function tabMI = motionCompMeanMI(meanMap,ROI)
%
%    gb 04/25/05
%
% Computes a matrix of MI from a cell array of mean maps
% tabMI(i,j) = motionCompMI(meanMap{i},ROI,meanMap{j});
%

if ieNotDefined('ROI')
    ROI = '';
end

n = length(meanMap);
tabMI = zeros(n);

for i = 1:n
    for j = 1:n
        if isempty(meanMap{i}) | isempty(meanMap{j})
            continue;
        end

        tabMI(i,j) = motionCompMI(meanMap{i},ROI,meanMap{j});
    end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%BOXFILTER  Blurs an image using the box filter
%
%     dst = cv.boxFilter(src)
%     dst = cv.boxFilter(src, 'OptionName',optionValue, ...)
%
% ## Input
% * __src__ input image.
%
% ## Output
% * __dst__ output image of the same size and type as `src`.
%
% ## Options
% * __DDepth__ the output image depth (-1 to use `class(src)`). Default -1.
%   See cv.filter2D for details.
% * __KSize__ blurring kernel size. Default [5,5]
% * __Anchor__ anchor point `[x,y]`; default value [-1,-1] means that the
%   anchor is at the kernel center.
% * __Normalize__ flag, specifying whether the kernel is normalized by its
%   area or not. default true
% * __BorderType__ border mode used to extrapolate pixels outside of the
%   image. See cv.copyMakeBorder. Default 'Default'
%
% The function smooths an image using the kernel:
%
%     K = alpha * ones(KSize)
%
% where:
%
%             | 1/prod(KSize)  when Normalize=true
%     alpha = |
%             | 1              otherwise
%
% Unnormalized box filter is useful for computing various integral
% characteristics over each pixel neighborhood, such as covariance matrices
% of image derivatives (used in dense optical flow algorithms, and so on). If
% you need to compute pixel sums over variable-size windows, use cv.integral.
%
% See also: cv.blur, cv.bilateralFilter, cv.GaussianBlur, cv.medianBlur,
%  cv.integral, cv.sqrBoxFilter, imfilter, imboxfilt
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


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

    fileID = fopen(fileName,'w');
    len = length(data(:,1));
    for ind = 1 : len
        fprintf(fileID,'%d %d \r\n',ind,data(ind, :));
    end
    fclose(fileID);
end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
% MultiQuadrics scattered data INTERPOLATION.
%
% SYNTAX: [fi, coefficients] = mq_interpolation( r, f, ri, c )
%
%          where r and f are the scattered argument and function
%          values, fi is the data interpolated at ri, c is the
%          multiquadric parameter and coefficients is a vector
%          containing the coefficients of the MQ expansion.
%
%         Examples:
%         N =  51; c = 0.2;
%         x = 2*pi*sort( rand([1 N]) ); f = sin(x);
%         xi = linspace(0,2*pi,N);
%        [fi,coefficients] = mq_interpolation(x(:),f(:),xi(:),c);
%         figure(1), plot(x,f,xi,fi,'o')
%
%	  N = 51; c = 0;
%	  x = 2*rand([1 N]) - 1; y = 2*rand([1 N]) - 1; r = [x(:) y(:)];
%	  for i = 1:N, f(i) = exp( -x(i)^2 -y(i)^2 ); end
%	  xi = linspace(-1.1,1.1,N); yi = linspace(-1.1,1.1,N);
%	 [Xi,Yi] = meshgrid(xi,yi);
%	  for i = 1:N*N, ri(i,:) = [Xi(i) Yi(i)]; end
%	 [fi,coefficients] = mq_interpolation(r,f(:),ri,c);
%	  S = reshape(fi,N,N);
%	  figure(2)
%	  mesh(xi,yi,S), hold on, plot3(x,y,f,'o'), hold off

% Written by Orlando Camargo Rodriguez

function [fi,coefficients] = mq_interpolation(r,f,ri,c)

fi           = [];
coefficients = [];

N = length( r(:,1) );

for i = 1:N

    rj_minus_ri = ( r(i,:)'*ones([1 N]) )' - r;

    Phi(:,i) = sqrt( sum( rj_minus_ri.^2 , 2 ) + c*c );

end

coefficients = Phi\f(:); % Calculate the coefficients through Gaussian elimination

Ni = length( ri(:,1) );

for i = 1:Ni

    ri_minus_r = ( ri(i,:)'*ones([1 N]) )' - r;

    phi_at_ri = sqrt( sum( ri_minus_r.^2 , 2 ) + c*c );

    fi(i) = coefficients'*phi_at_ri(:);

end

author: Sable

repo: mcbench-benchmarks


text:
function [ sparse_order, sparse_index ] = sparse_grid_mixed_index ( ...
  dim_num, level_max, rule, point_num, point_total_num, sparse_unique_index )

%*****************************************************************************80
%
%% SPARSE_GRID_MIXED_INDEX indexes a sparse grid made from mixed 1D rules.
%
%  Discussion:
%
%    For each "unique" point in the sparse grid, we return its INDEX and ORDER.
%
%    That is, for the I-th unique point P, we determine the product grid which
%    first generated this point, and we return in SPARSE_ORDER the orders of
%    the 1D rules in that grid, and in SPARSE_INDEX the component indexes in
%    those rules that generated this specific point.
%
%    For instance, say P was first generated by a rule which was a 3D product
%    of a 9th order CC rule and a 15th order GL rule, and that to generate P,
%    we used the 7-th point of the CC rule and the 3rh point of the GL rule.
%    Then the SPARSE_ORDER information would be (9,15) and the SPARSE_INDEX
%    information would be (7,3).  This, combined with the information in RULE,
%    is enough to regenerate the value of P.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 March 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Input, integer RULE(DIM_NUM), the rule in each dimension.
%     1, "CC",  Clenshaw Curtis, Closed Fully Nested rule.
%     2, "F2",  Fejer Type 2, Open Fully Nested rule.
%     3, "GP",  Gauss Patterson, Open Fully Nested rule.
%     4, "GL",  Gauss Legendre, Open Weakly Nested rule.
%     5, "GH",  Gauss Hermite, Open Weakly Nested rule.
%     6, "GGH", Generalized Gauss Hermite, Open Weakly Nested rule.
%     7, "LG",  Gauss Laguerre, Open Non Nested rule.
%     8, "GLG", Generalized Gauss Laguerre, Open Non Nested rule.
%     9, "GJ",  Gauss Jacobi, Open Non Nested rule.
%    10, "GW",  Golub Welsch, (presumed) Open Non Nested rule.
%    11, "CC_SE", Clenshaw Curtis Slow Exponential, Closed Fully Nested rule.
%    12, "F2_SE", Fejer Type 2 Slow Exponential, Closed Fully Nested rule.
%    13, "GP_SE", Gauss Patterson Slow Exponential, Closed Fully Nested rule.
%    14, "CC_ME", Clenshaw Curtis Moderate Exponential, Closed Fully Nested rule.
%    15, "F2_ME", Fejer Type 2 Moderate Exponential, Closed Fully Nested rule.
%    16, "GP_ME", Gauss Patterson Moderate Exponential, Closed Fully Nested rule.
%    17, "CCN", Clenshaw Curtis Nested, Linear, Closed Fully Nested rule.
%
%    Input, integer POINT_NUM, the number of unique points in the grid.
%
%    Input, integer POINT_TOTAL_NUM, the total number of points in the grid.
%
%    Input, integer SPARSE_UNIQUE_INDEX(POINT_TOTAL_NUM), associates each
%    point in the grid with its unique representative.
%
%    Output, integer SPARSE_ORDER(DIM_NUM,POINT_NUM), lists, for each point,
%    the order of the 1D rules used in the grid that generated it.
%
%    Output, integer SPARSE_INDEX(DIM_NUM,POINT_NUM), lists, for each point,
%    its index in each of the 1D rules in the grid that generated it.
%    The indices are 1-based.
%

%
%  Special cases.
%
  if ( level_max < 0 )
    sparse_order = [];
    sparse_index = [];
    return
  end

  if ( level_max == 0 )
    sparse_order(1:dim_num,1) = 1;
    sparse_index(1:dim_num,1) = 1;
    return
  end

  sparse_order = zeros ( dim_num, point_num );
  sparse_index = zeros ( dim_num, point_num );

  point_count = 0;
%
%  The outer loop generates values of LEVEL.
%
  level_min = max ( 0, level_max + 1 - dim_num );

  for level = level_min : level_max
%
%  The middle loop generates a GRID,
%  based on the next partition that adds up to LEVEL.
%
    level_1d = [];
    more_grids = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more_grids, h, t ] = comp_next ( level, dim_num, level_1d, ...
        more_grids, h, t );

      order_1d = level_to_order_default ( dim_num, level_1d, rule );
%
%  The inner loop generates a POINT of the GRID of the LEVEL.
%
      point_index = [];
      more_points = 0;

      while ( 1 )

        [ point_index, more_points ] = vec_colex_next3 ( dim_num, order_1d, point_index, ...
          more_points );

        if ( ~more_points )
          break
        end

        point_count = point_count + 1;
        point_unique = sparse_unique_index(point_count);
        sparse_order(1:dim_num,point_unique) = order_1d(1:dim_num);
        sparse_index(1:dim_num,point_unique) = point_index(1:dim_num);

      end

      if ( ~more_grids )
        break
      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
clear all, close all, clc

m = 1;
M = 5;
L = 2;
g = -10;
d = 1;

s = 1; % pendulum up (s=1)
A = [0 1 0 0;
    0 -d/M -m*g/M 0;
    0 0 0 1;
    0 -s*d/(M*L) -s*(m+M)*g/(M*L) 0];

B = [0; 1/M; 0; s*1/(M*L)];

C = [1 0 0 0];

D = zeros(size(C,1),size(B,2));

Q = [1 0 0 0;
    0 1 0 0;
    0 0 1 0;
    0 0 0 1];
R = .000001;
K = lqr(A,B,Q,R);  % design controller u = -K*x

%  Augment system with disturbances and noise
Vdmag = .04;
Vd = Vdmag*eye(4);  % disturbance covariance
Vn = .0002;       % noise covariance

BF = [B Vd 0*B];  % augment inputs to include disturbance and noise

sysC = ss(A,BF,C,[0 0 0 0 0 Vn]);  % build big state space system... with single output

sysFullOutput = ss(A,BF,eye(4),zeros(4,size(BF,2)));  % system with full state output, disturbance, no noise

%  Build Kalman filter
[L,P,E] = lqe(A,eye(4),C,Vd,Vn);  % design Kalman filter
Kf = (lqr(A',C',Vd,Vn))';   % alternatively, possible to design using "LQR" code

sysKF = ss(A-Kf*C,[B Kf],eye(4),0*[B Kf]);  % Kalman filter estimator

author: dynamicslab


text:
function oeprint(mu, oev)

% print six classical orbital elements
% (orbital period in days, sma in kilometers)

% input

%  mu      = gravitational constant (km^3/sec^2)
%  oev(1)  = semimajor axis (km)
%  oev(2)  = orbital eccentricity (non-dimensional)
%            (0 <= eccentricity < 1)
%  oev(3)  = orbital inclination (radians)
%            (0 <= inclination <= pi)
%  oev(4)  = argument of periapsis (radians)
%            (0 <= argument of periapsis <= 2 pi)
%  oev(5)  = right ascension of ascending node (radians)
%            (0 <= raan <= 2 pi)
%  oev(6)  = true anomaly (radians)
%            (0 <= true anomaly <= 2 pi)

% Orbital Mechanics with MATLAB

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

rtd = 180 / pi;

% unload orbital elements array

sma = oev(1);
ecc = oev(2);
inc = oev(3);
argper = oev(4);
raan = oev(5);
tanom = oev(6);

arglat = mod(tanom + argper, 2.0 * pi);

period = 2.0 * pi * sma * sqrt(sma / mu) / 86400.0;

% print orbital elements

fprintf ('\n      sma (km)        eccentricity     inclination (deg)    argper (deg)');

fprintf ('\n %12.10e  %12.10e  %12.10e  %12.10e \n', sma, ecc, inc * rtd, argper * rtd);

fprintf ('\n     raan (deg)     true anomaly (deg)   arglat (deg)       period (days)');

fprintf ('\n %12.10e  %12.10e  %12.10e  %12.10e \n', raan * rtd, tanom * rtd, arglat * rtd, period);


text:
classdef  VideoMaker_TopOpt_density < VideoMaker_TopOpt

    methods (Access = public)

        function obj = VideoMaker_TopOpt_density(cParams)
            obj.init(cParams);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function obj = saveObj(obj)
	% Biafra Ahanonu
	% Started: 2021.03.25 [22:11:25] (branched from ciatah.m)

	if isempty(obj.objSaveLocation)
		[filePath,folderPath,~] = uiputfile('*.*','select folder to save object mat file to','CIAtah_properties.mat');
		% exit if user picks nothing
		% if folderListInfo==0; return; end
		savePath = [folderPath filesep filePath];
		% tmpObj = obj;
		% obj = struct(obj);
		obj.objSaveLocation = savePath;
	else
		savePath = obj.objSaveLocation;
	end
	disp(['saving to: ' savePath])
	try
	  save(savePath,'obj','-v7.3');
	catch
	  disp('Problem saving, choose new location...')
	  obj.objSaveLocation = [];
	  obj.saveObj();
	end
	% obj = tmpObj;
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% runmcs.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% test driver for running MCS
% with advice on how to choose the tuning parameters
%
% applied to the Jones test functions with default boxes
%
% To solve your own problem, copy this file (to ownmcs.m, say)
% and adapt the first part labelled `problem definition'.
% Then run the driver by typing `ownmcs' at the Matlab prompt.
%
% If you are not satisfied with the results, or if you want to play
% with the tuning parameters, you also need to adapt the second part
% labelled `change MCS settings'. In particular, for wide bounds,
% you'll probably get better results if you supply your own
% initialization list.
%
% On typing `runmcs' at the Matlab prompt,
% the unmodified file produces test results for the six-hump camel
% function; by only changing the data assignment you can also get
% results for the other test functions from Jones et al.
% You may also play with the bounds by modifying the default bounds.
%

clear; clear mex; clear global;
format compact;format long

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% problem definition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define objective function
%
% Each function value f=f(x) is obtained by MCS using a call
% f=feval(fcn,data,x)
% where data is an arbitrary data object passed to mcs. If you have
% several data items, put them into a cell array or structure array
% called data.
% If there are no data, use fcn='feval' and the function name as data.
%
fcn = 'feval';
data = 'cam';	% select test function from Jones et al. test set
		% bra  n=2   Branin
		% cam  n=2   six-hump camel
		% gpr  n=2   Goldstein-Price
		% shu  n=2   Shubert
		% hm3  n=3   Hartman 3
		% s10  n=4   Shekel 10
		% sh5  n=4   Shekel 5
		% sh7  n=4   Shekel 7
		% hm6  n=6   Hartman 6
path(path,'jones');	% add path to directory with function

% define bounds on variables (+-inf allowed)
%
% u: column vector of lower bounds
% v: column vector of upper bounds
% u(k)<v(k) is required
%
[u,v,fglob] = defaults(data); 	% returns bounds used by Jones et al.
				% and known global optimum
dimension=length(u)		% show dimension
known_global_opt_value=fglob	% show known global minimum value
% u=[-5,0]';v=[10,15]';    		% bra default bounds
% u=[-3,-2]';v=[3,2]';     		% cam default bounds
% u=[-2,-2]';v=[2,2]';     		% gpr default bounds
% u=[-10,-10]';v=[10,10]'; 		% shu default bounds
% u=[0,0,0]';v=[1,1,1]';   		% hm3 default bounds
% u=[0,0,0,0]';v=[10,10,10,10]';   	% sh5,sh7,s10 default bounds
% u=[0,0,0,0,0,0]';v=[1,1,1,1,1,1]';   	% hm6 default bounds
% modify the problem to be unconstrained by activating the next line
% u = -Inf*ones(size(u));v=-u;

use_defaults=1
% *** If you just want to use the default settings,
% *** you don't need to edit the rest of the file,
% *** If you are not satisfied with the results
% *** (this may happen especially when your box bounds are very wide),
% *** or if you want to play with the tuning parameters,
% *** set use_defaults=0, and modify the rest of the file
% *** according to your curiosity or ingenuity

if use_defaults,
  % easy to use version - all parameters preset
  % defaults are being used, it suffices to call
  %%%%%%%%%%%%%%%%%% simple MCS call %%%%%%%%%%%%%%%%%%
  [xbest,fbest,xmin,fmi,ncall,ncloc]=mcs(fcn,data,u,v);
  % or, with less output,
  % [xbest,fbest]=mcs(fcn,data,u,v);

  xbest	  		% best point found
  fbest     		% best function value
  fglob			% global minimum (known for test functions)
  ncall	  		% number of function values used
  ncloc	  		% number of function values in local searches

  % xmin	  	% columns are points in 'shopping basket'
			% may be good alternative local minima
  % fmi	  		% function values in 'shopping basket'
  nbasket = length(fmi) % number of points in 'shopping basket'

  return;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% change MCS settings %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% flexible use version - all parameters may be modified

%
% define amount of output printed
prt = 1;	% print level
		% prt = 0: no output
		% prt = 1: # sweep, minimal nonempty level, # f-calls,
		%          best point and function value
		% prt > 1: in addition levels and function values of
		%          boxes containing the known global minimizers
		%          of a test function

% define global strategy
%
% smax governs the relative amount of global versus local search. By
% increasing smax, more weight is given to global search.
%
% Increasing or decreasing stop(1) increases or decreases the amount
% of work allowed before concluding that nothing more is gained;
% the default choice is quite conservative, and may try to locate
% a better point long after the global minimum has been found.
% stop(1)=5 works for many easier problems, too, with much fewer
% wasted function values.
%
% Increasing nf allows to spend more time on boxes that have a chance
% on their level to contain better points. This may be important for
% hard problems, or for problems with very wide bounds.
%
% But in each case, it is unclear in advance what change would be most
% beneficial to a particular problem.
% We had very mixed experience; if you have many similar problems to
% solve, the best thing to do is to experiment with a few problems to
% find the best values, and to use these on the others.
%
n = length(u);		% problem dimension
smax = 5*n+10;		% number of levels used
nf = 50*n^2; 		% limit on number of f-calls
stop(1) = 3*n;		% = m, integer defining stopping test
stop(2) = -inf;		% = freach, function value to reach
			% if m>0, run until m sweeps without progress
			% if m=0, run until fbest<=freach
			% (or about nf function calls were used)

if 0, 	% known global optimum, for tests only
	% then the entries of stop have a different meaning
  stop(1) = 1.e-4;	% run until this relative error is achieved
  stop(2) = fglob;	% known global optimum value
  stop(3) = 1.e-10;	% stopping tolerance for tiny fglob
end;

% define initialization strategy
%
% for wide boxes, it is advisable (and for unbounded search regions
% strongly advisable) to define a customized initialization list
% that contains for each coordinate at least three reasonable values.
% Without such an initialization list, mcs makes default assumptions
% that roughly amount to estimating reasonable magnitudes from the
% bounds and in case iinit=1 from assuming order unity if this is
% within the bounds.
%
% for a self-defined initialization list, the user should
% write an m-script file init0.m containing a matrix x0 with n
% rows and at least 3 columns and two n-vectors l and L
% the ith column of x0 contains the initialization list
% values for the ith coordinate, their number is L(i), and
% x0(i,l(i)) is the ith coordinate of the initial point

iinit = 0;	% 0: simple initialization list
		%    (default for finite bounds;
		%     do not use this for very wide bounds)
		% 1: safeguarded initialization list
		%    (default for unbounded search regions)
		% 2: (5*u+v)/6, (u+v)/2, (u+5*v)/6
		% 3: initialization list with line searches
		% else: self-defined initialization list
		%       stored in file init0.m

% parameters for local search
%
% A tiny gamma (default) gives a quite accurate but in higher
% dimensions slow local search. Increasing gamma leads to less work
% per local search but a less accurately localized minimizer
%
% If it is known that the Hessian is sparse, providing the sparsity
% pattern saves many function evaluations since the corresponding
% entries in the Hessian need not be estimated. The default pattern
% is a full matrix.
%
local = 50;		% local = 0: no local search
			% else: maximal number of steps in local search
gamma = eps;		% acceptable relative accuracy for local search
hess = ones(n,n);	% sparsity pattern of Hessian

% defaults are not being used, use the full calling sequence
% (including at least the modified arguments)
%%%%%%%%%%%%%%%%%%%%%%% full MCS call %%%%%%%%%%%%%%%%%%%%%%
[xbest,fbest,xmin,fmi,ncall,ncloc]=...
  mcs(fcn,data,u,v,prt,smax,nf,stop,iinit,local,gamma,hess);

xbest	  		% best point found
fbest     		% best function value
fglob			% global minimum (known for test functions)
ncall	  		% number of function values used
ncloc	  		% number of function values in local searches

% xmin	  		% columns are points in 'shopping basket'
			% may be good alternative local minima
% fmi	  		% function values in 'shopping basket'
nbasket = length(fmi) 	% number of points in 'shopping basket'

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function test_ChallengeDB_SVR(testDatabase, testMethod, Trainingproportion)

    %% set up database
    if ~exist(fullfile('result', testDatabase, testMethod), 'dir')
        mkdir(fullfile('result', testDatabase, testMethod));
    end

    fprintf('processing LIVE ChallengeDB images.....');
    Path = fullfile('databases', testDatabase);
    imageList = load(fullfile(Path,'Data','AllImages_release.mat'));
    mos = load(fullfile(Path,'Data','AllMOS_release.mat'));
    imageList = fullfile(Path,'Images',imageList.AllImages_release(8:end));
    mos = mos.AllMOS_release(8:end);

    %% feature extraction
    if ~exist(fullfile('result', testDatabase, testMethod, 'features.mat'), 'file')
        for i = 1:length(mos)
            fprintf('processing image %d / %d \n',i,length(mos));
            features{i} = CallingTheSourceCodeForFeatureExtraction(imageList{i},testMethod);
        end
        features = cat(2,features{:});
        save(fullfile('result', testDatabase, testMethod, 'features.mat'),'features');
    else
        load(fullfile('result', testDatabase, testMethod, 'features.mat'),'features');
    end

    %% train SVR and test
    tic;
    trainingNum = ceil(length(mos)*Trainingproportion);
    for seed = 1:10
        rng(seed-1);
        randomsplit = randperm(length(mos),length(mos));
        trainingSet = randomsplit(1:trainingNum);
        testingSet = randomsplit(trainingNum+1:end);
        trainX = features(:,trainingSet)';
        testX = features(:,testingSet)';
        [trainX, testX] = featNormalize(trainX, testX);
        switch testMethod
            case 'BRISQUE'
                svm_model = svmtrain(mos(trainingSet)', double(trainX), '-s 4 -t 2 -c 10000 -g 0.05');
            case 'CORNIA'
                svm_model = svmtrain(mos(trainingSet)', double(trainX), '-s 4 -t 0');
            case 'HOSA'
                svm_model = train(mos(trainingSet)', sparse(trainX), '-s 11 -c 128');
            case 'DIIVINE'
                svm_model = svmtrain(mos(trainingSet)', double(trainX), '-s 4 -t 2 -c 4096 -g 0.02');
            case 'FRIQUEE'
                svm_model = svmtrain(mos(trainingSet)', double(trainX), '-s 4 -t 2 -c 256 -g 0.02');
            otherwise
                error('Unsupported Method!');
        end
        if strcmp(testMethod,'HOSA')
            [scores, ~,~] = predict(mos(testingSet)', sparse(testX), svm_model,'-q');
        else
            [scores, ~,~] = svmpredict(mos(testingSet)', double(testX), svm_model,'-q');
        end

        SRCC(seed) = corr(scores, mos(testingSet)', 'type', 'Spearman');
        PLCC(seed) = RegressionLIVE(scores, mos(testingSet)');
    end
    toc;

    fprintf('SRCC = %.4f, std = %.4f;  PLCC = %.4f, std = %.4f \n', abs(median(SRCC)), std(SRCC), abs(median((PLCC))), std((PLCC)));

end

author: HuiZeng


text:
%kkstats 'Compute Statistics of Data Object'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kstats.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'Input object'
% InputFile: igate 'Gating Input', optional: 'Gating input data object'
% OutputFile: f 'ASCII Output ', optional: 'Formatted ASCII output file'
% OutputFile: o 'Binary Output', optional: 'Binary output statistics file'
% Toggle: whole 'Whole Data Set', default: 0: 'compute single set of statistics for entire data set'
% Toggle: w 'Width ', default: 0: 'include width in processing unit'
% Toggle: h 'Height ', default: 0: 'include height in processing unit'
% Toggle: d 'Depth ', default: 0: 'include depth in processing unit'
% Toggle: t 'Time ', default: 0: 'include time in processing unit'
% Toggle: e 'Element', default: 0: 'include elements in processing unit'
% Toggle: all 'Calculate ALL Statistics', default: 0: 'compute and store all statistics information'
% Toggle: mean 'mean', default: 0: 'compute mean'
% Toggle: sum 'total integral', default: 0: 'compute total integral'
% Toggle: var 'variance', default: 0: 'compute variance'
% Toggle: wmin 'width coordinate', default: 0: 'store width coordinate of the minimum value'
% Toggle: psum 'pos part of integral', default: 0: 'compute positive part of integral'
% Toggle: sd 'std deviation', default: 0: 'compute standard deviation'
% Toggle: hmin 'height coordinate', default: 0: 'store height coordinate of the minimum value'
% Toggle: nsum 'neg part of integral', default: 0: 'compute negative part of integral'
% Toggle: rms 'rms', default: 0: 'compute root mean square'
% Toggle: dmin 'depth coordinate', default: 0: 'store depth coordinate of the minimum value'
% Toggle: pts 'total contributing pts', default: 0: 'compute total number of contributing points'
% Toggle: skew 'skewness', default: 0: 'compute skewness'
% Toggle: tmin 'time coordinate', default: 0: 'store time coordinate of the minimum value'
% Toggle: ppts 'positive points', default: 0: 'compute number of positive contributing points'
% Toggle: kur 'kurtosis', default: 0: 'compute kurtosis'
% Toggle: emin 'elements coordinate', default: 0: 'store elements coordinate of the minimum value'
% Toggle: npts 'negative points', default: 0: 'compute number of negative contributing points'
% Toggle: minval 'minimum', default: 0: 'compute minimum value'
% Toggle: zpts 'zero points', default: 0: 'compute number of zero-valued contributing pts'
% Toggle: maxval 'maximum', default: 0: 'compute maximum value'
% Toggle: wmax 'width coordinate', default: 0: 'store width coordinate of the maximum value'
% Toggle: wsize 'data width', default: 0: 'store size of objects width dimension'
% Toggle: hmax 'height coordinate', default: 0: 'store height coordinate of the maximum value'
% Toggle: hsize 'data height', default: 0: 'store size of objects height dimension'
% Toggle: dmax 'depth coordinate', default: 0: 'store depth coordinate of the maximum value'
% Toggle: dsize 'data depth', default: 0: 'store size of objects depth dimension'
% Toggle: tmax 'time coordinate', default: 0: 'store time coordinate of the maximum value'
% Toggle: tsize 'time size', default: 0: 'store size of objects time dimension'
% Toggle: emax 'elements coordinate', default: 0: 'store elements coordinate of the maximum value'
% Toggle: esize 'elements size', default: 0: 'store size of objects elements dimension'
%
% Example: [f, o] = kkstats({i, igate}, {'i','';'igate','';'f','';'o','';'whole',0;'w',0;'h',0;'d',0;'t',0;'e',0;'all',0;'mean',0;'sum',0;'var',0;'wmin',0;'psum',0;'sd',0;'hmin',0;'nsum',0;'rms',0;'dmin',0;'pts',0;'skew',0;'tmin',0;'ppts',0;'kur',0;'emin',0;'npts',0;'minval',0;'zpts',0;'maxval',0;'wmax',0;'wsize',0;'hmax',0;'hsize',0;'dmax',0;'dsize',0;'tmax',0;'tsize',0;'emax',0;'esize',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kstats - Compute Statistics of Data Object
%
%  DESCRIPTION
% .I kstats
% computes the mean, variance, standard deviation, RMS level, skew, kurtosis,
% minimum value, maximum value, and integral (positive, negative, and total
% sums) of the Input data object (i).  Information such as the number of
% unmasked or ungated points (positive, negative, zero-valued, and total counts),
% the coordinates of the minimum and maximum values, and the dimensions of the
% data object are also provided.
% .so $DATAMANIP/repos/shared/man/sections/mask_stats
% .so $DATAMANIP/repos/shared/man/sections/map_force
% All statistics and information will reflect the mapping.  For example,
% if the dimensions of the data are requested (wsize, hsize, etc), kstats
% returns the dimensionality after mapping.
% .cI $DATAMANIP/repos/shared/man/sections/fail_1input
% All statistics are calculated as double, and the Binary Output data type
% is double.
%
% Location and time data do not affect the statistics calculations, and are
% not transferred to the output object.
%
% A flag can be set for each statistic and information option so that a
% subset of all available options can be calculated and stored or printed.
% The flags are mutually exclusive with the "Calculate All Statistics"
% flag (all).  (The command line flag names are: mean, var, sd, rms, skew,
% kur, minval, maxval, wmin, hmin, dmin, tmin, emin, wmax, hmax, dmax, tmax,
% emax, sum, psum, nsum, pts, ppts, npts, zpts, wsize, hsize, dsize, tsize,
% and esize.) If no flag is specified when running kstats from the command
% line, all statistics and information will be stored.
%
% Statistics are computed according to the equations given below.  In the
% equations, N is the number of samples, SUM is the sum from i=0 to i=N-1,
% and x(i) is the sample value of x at i.
%
%  "Mean (-mean)"
% mean = (1/N) * SUM(x(i))  i=0..N-1
%
%  "Variance (-var)"
% variance = (1/(N-1)) * (SUM (x(i) - mean)**2)  i=0..N-1
%
%  "Standard Deviation (-sd)"
% standard deviation = sqrt(variance)
%
%  "RMS (-rms)"
% RMS  = sqrt(1/N * SUM(x(i)**2))  i=0..N-1
%
%  "Skewness (-skew)"
% Skewness is a measure of the tendency of the deviations from the mean to be
% larger in one direction than in the other.  It is a measure of the asymmetry
% of a distribution about its mean.  A positive skew value signifies a
% distribution whose tail extends out towards more positive x, and a negative
% tail signifies a distribution whose tail extends out towards more negative x.
% Population skewness is unitless and is defined as:
% E[((x-mean)/(stddev))**3],
% where stddev is the standard deviation of the data.
% .I kstats
% computes the sample skewness as:
% skew = 1/N * SUM( ((x(i) - mean)/stddev) **3 )  i=0..N-1
% If the variance is equal to zero, skewness will be set to 0.0.
%
%  "Kurtosis"
% Kurtosis is a unitless measure of the tail heaviness of a distribution
% and is defined as:
% E[((x-mean)/(stddev))^4] - 3,
% where stddev is the standard deviation of the data.
% .I kstats
% computes the sample kurtosis as:
% kurtosis= (1/N * SUM( ((x(i) - mean)/stddev) **4)) - 3  i=0..N-1
% If the variance is equal to zero, kurtosis will be set to 0.0.
%
%  "Sums or Integrals (-sum, -psum, -nsum)"
% total integral = SUM x(i)
% positive part of integral = SUM x(i),  for all x(i) > 0
% negative part of integral = SUM x(i),  for all x(i) < 0
%
%  "Define Processing Unit"
% To support analysis of subregions within the data object, such as lines,
% planes, volumes, and vectors, an option for defining processing units, or
% regions, is provided.  These regions are defined by the settings of the
% Processing Unit options, which can be either the Whole Data Set (whole),
% or any combination of Width (w), Height (h), Depth (d), Time (t), and
% Elements (e).  Statistics for each region will be computed and printed
% separately. If none of these flags are supplied, then a single set of
% statistics will be computed for the entire data object.
%
%  "Gating Input"
% .so $DATAMANIP/repos/shared/man/sections/gate_simple
%
%  "ASCII Output"
% The "ASCII Output" (f) option allows the user to specify a device or file
% for printing the specified information in formatted ASCII.  A filename
% of # will send the output to stderr.  If neither an "ASCII Output"
% or a "Binary Output" is supplied, the formatted ASCII will automatically
% go to stdout.
%
%  "Binary Output"
% If the "Binary Output" (-o) is selected, the selected statistics and
% information are stored in as double float data in the given file.
% Each statistic is stored as an element of a N-D vector defined along the
% "elements" dimension of the output object.  N is the number of statistics
% and information specified by the user.  The order in which the information
% is stored in the statistics vector is given in the output object's comment
% attribute.
% When M sets of statistics are calculated for M multiple regions of the
% input data object (see discussion on independent region analysis above),
% the statistics vectors are incremented along the width dimension of the
% output statistics object - the output object would have a resulting
% dimensionality of width=M, elements=N (height = depth = time = 1).  For
% example, if the input data object had the dimensionality width = 256,
% height = 256, elements = 7 (depth = time = 1), and the user defined the
% analysis regions to be width-height, the output statistics object would
% have dimensionality width = 7, elements = N (height = depth = time = 1).
% Where the vector at w=0 would contain information pertaining to the first
% region (starting at w=0, h=0, d=0, t=0, e=0), the vector at w=1 would
% contain information about the second region (starting at w=0, h=0, d=0, t=0,
% e=1), and so forth.
%
%
%
%  EXAMPLES
% kstats -i object.xv -f ascii -all
%
% Will create an ASCII statitics file.  Since the -region option is not specified,
% kstats defaults to calculating statistical information for the object
% as a whole.
% kstats -i object.xv -f ascii_file -region -e
%
% Will create an ASCII statistics file. The statistical information is now
% calculated by slicing the object along the element dimension. For example,
% on a 3-band RGB color image, this would compute the statistics for each
% RGB vector.
%
%  "SEE ALSO"
%
%  RESTRICTIONS
% kstats cannot compute complex statistics at this time.  If you wish to
% process the real and imaginary components of a complex data set
% independently, first separate components using kcmplx2real.
%
% In the current implementation, the statistics are accumulated.  This
% can result in overflow.  kstats will be rewritten in the future
% to better handle large values and large data sets.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkstats(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkstats(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'igate', '__input';'f', '__output';'o', '__output';'whole', 0;'w', 0;'h', 0;'d', 0;'t', 0;'e', 0;'all', 0;'mean', 0;'sum', 0;'var', 0;'wmin', 0;'psum', 0;'sd', 0;'hmin', 0;'nsum', 0;'rms', 0;'dmin', 0;'pts', 0;'skew', 0;'tmin', 0;'ppts', 0;'kur', 0;'emin', 0;'npts', 0;'minval', 0;'zpts', 0;'maxval', 0;'wmax', 0;'wsize', 0;'hmax', 0;'hsize', 0;'dmax', 0;'dsize', 0;'tmax', 0;'tsize', 0;'emax', 0;'esize', 0};
maxval={0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
minval={0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
istoggle=[0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','OutputFile','OutputFile','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kstats"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function a = pascal3_inverse ( n, alpha )

%*****************************************************************************80
%
%% PASCAL3_INVERSE returns the inverse of the PASCAL3 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real ALPHA, the parameter.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n

      if ( i == 1 )

        if ( j == 1 )
          a(i,j) = 1.0;
        else
          a(i,j) = 0.0;
        end

      elseif ( j == 1 )

        a(i,j) = - alpha * a(i-1,j);

      else

        a(i,j) = a(i-1,j-1) - alpha * a(i-1,j);

      end

    end
  end

  return
end


text:
function corrplot(c)

% Called internally by correlation/plot to plot correlation matrix.

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

% PREP PLOT
figure('Color','w','Position',[50 50 600 500]);
set(gcf,'DefaultAxesFontSize',14);
imagesc(c.C);
title('Maximum correlation coefficient');

% ADD DATES TO AXES
n = length(c.trig);
set(gca,'XTick',[1:round(n/25):n]);
set(gca,'YTick',[1:round(n/25):n]);
yt = get(gca,'YTick');
set(gca,'YTickLabel',datestr(c.trig(yt),'yyyy-mm-dd HH:MM'),'FontSize',6);

%xt = get(gca,'XTick');
%set(gca,'XTickLabel',datestr(c.trig(xt),'yyyy-mm-dd HH:MM'),'FontSize',6,'Rotation',90);

% DRESS UP THE FIGURE
caxis([0 1]);
cmap = load('colormap_corr.txt');
colormap(cmap);
colorbar;
xlabel('Event number');
ylabel('Event date');

%PRINT OUT FIGURE
set(gcf, 'paperorientation', 'portrait');
set(gcf, 'paperposition', [1.25 2.5 6 6] );
%print(gcf, '-depsc2', 'FIG_tartan.ps');

author: geoscience-community-codes

repo: GISMO


text:
function HDR=bv2biosig_events(EVENT)
% BV2BIOSIG_EVENTS converts VMRK marker information BioSig Event codes.
%  according to biosig/doc/eventcodes.txt.
%  Currently, the convention of the BerlinBCI is implemented and supported.
%
%  HDR = bv2biosig_events(arg1)
%
%  arg1 can be an HDR-struct containg  HDR.EVENT.Desc
%  or a struct containing EVENT.Desc
%  or a cell-array Desc
%  or a char-array Desc
%
%  Warning: Approximately 32 (out of 510) events are currently not supported.
%  For your data, you can check this with this command:
%  HDR.EVENT.Desc(HDR.EVENT.TYP==0)
%
% see also: doc/eventcodes.txt

%	$Id$
%	Copyright (C) 2006,2007 by Alois Schloegl <a.schloegl@ieee.org>
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/

% This library is free software; you can redistribute it and/or
% modify it under the terms of the GNU Library General Public
% License as published by the Free Software Foundation; either
% Version 2 of the License, or (at your option) any later version.
%
% This library 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
% Library General Public License for more details.
%
% You should have received a copy of the GNU Library General Public
% License along with this library; if not, write to the
% Free Software Foundation, Inc., 59 Temple Place - Suite 330,
% Boston, MA  02111-1307, USA.

HDR.EVENT.Desc = [];
if isfield(EVENT,'EVENT')
	HDR = EVENT;
elseif isfield(EVENT,'Desc')
	HDR.EVENT = EVENT;
elseif iscell(EVENT)
	HDR.EVENT.Desc = EVENT;
elseif ischar(EVENT)
	for k=1:size(EVENT,1),
		HDR.EVENT.Desc{k} = EVENT(k,:);
	end;
else
	error('unknown input argument')
end;

if isfield(HDR.EVENT,'TeegType')
	ix = strmatch('New Segment',HDR.EVENT.TeegType);
	HDR.EVENT.TYP(ix)=hex2dec('7ffe');
end;

if ~isfield(HDR,'NS')
	HDR.NS = NaN;
end;
if ~isfield(HDR,'Label'),
	HDR.Label = repmat({' '},min(HDR.NS,128),1);
end;

FLAG_SEASON2_ARTERAWDATA = 0;
if (HDR.NS==128)
	tmp = strvcat(HDR.Label);
	tmp(55:56,4) = 'f';
	FLAG_SEASON2_ARTERAWDATA = isequal(tmp(1:64,1:end-1),tmp(65:128,2:end)) & all(tmp(65:128,1)=='x');
	FLAG_SEASON2_ARTERAWDATA = FLAG_SEASON2_ARTERAWDATA & strncmp(HDR.FILE.Name,'arte',4);
end;
FLAG_ARTERAWDATA = strncmp(HDR.FILE.Name,'arte',4);

for k1 = 1:length(HDR.EVENT.POS)
	if isfield(HDR.EVENT,'Desc')
		tmp = HDR.EVENT.Desc{k1};
	elseif isfield(HDR.EVENT,'CodeDesc') && (HDR.EVENT.TYP(k1)<=length(HDR.EVENT.CodeDesc))
		tmp = HDR.EVENT.CodeDesc{HDR.EVENT.TYP(k1)};
	else
		continue;
	end;

	if 0,

	elseif strncmp(tmp,'TargetCode',10)
		HDR.EVENT.TYP(k1) = str2double(tmp(11:12))+hex2dec('0300');

	elseif strcmp(tmp,'BeginOfTrial')
		HDR.EVENT.TYP(k1) = hex2dec('0300');

	elseif strcmp(tmp,'hit')
        	HDR.EVENT.TYP(k1) = hex2dec('0381');
        elseif strcmp(tmp,'wrong')
        	HDR.EVENT.TYP(k1) = hex2dec('0382');

% eye movements
        elseif strcmpi(tmp,'augen links')
        	HDR.EVENT.TYP(k1) = hex2dec('0431');
        elseif strcmpi(tmp,'augen rechts')
        	HDR.EVENT.TYP(k1) = hex2dec('0432');
        elseif strcmpi(tmp,'augen hoch') | strcmpi(tmp,'augen oben')
        	HDR.EVENT.TYP(k1) = hex2dec('0433');
        elseif strcmpi(tmp,'augen unten') | strcmpi(tmp,'augen runter')
        	HDR.EVENT.TYP(k1) = hex2dec('0434');
        elseif strcmpi(tmp,'augen offen') | strcmp(tmp,'Augen offen & entspannen')
        	HDR.EVENT.TYP(k1) = hex2dec('0114');
        elseif strcmpi(tmp,'augen zu') | strcmp(tmp,'Augen zu & entspannen')
        	HDR.EVENT.TYP(k1) = hex2dec('0115');
        elseif strcmp(tmp,'blinzeln')
        	HDR.EVENT.TYP(k1) = hex2dec('0439');

% muscle movements
        elseif strcmp(tmp,'EMG links')
        	HDR.EVENT.TYP(k1) = hex2dec('0441');
        elseif strcmp(tmp,'EMG rechts')
        	HDR.EVENT.TYP(k1) = hex2dec('0442');
        elseif strcmpi(tmp,'kopf bewegen')
        	HDR.EVENT.TYP(k1) = hex2dec('0443');
        elseif strcmp(tmp,'zunge an')
        	HDR.EVENT.TYP(k1) = hex2dec('0444');
        elseif strcmp(tmp,'zunge aus')
        	HDR.EVENT.TYP(k1) = hex2dec('8444');
        elseif strcmp(tmp,'Kiefer anspannen')
        	HDR.EVENT.TYP(k1) = hex2dec('0446');
        elseif strcmp(tmp,'beißen') | strcmp(tmp,'beissen') | strcmp(tmp,['bei',223,'en']),
        	HDR.EVENT.TYP(k1) = hex2dec('0446');
        elseif strcmp(tmp,'EMG fuss')
		HDR.EVENT.TYP(k1) = hex2dec('0447');
	elseif strcmp(tmp,'Arme bewegen')
		HDR.EVENT.TYP(k1) = hex2dec('0449');

% encoding der season2-arte* rawdata records
	elseif strncmp(tmp,'S',1) & (FLAG_SEASON2_ARTERAWDATA | FLAG_ARTERAWDATA),
		n = str2double(tmp(2:end));
		if n==11,	% EMG left
			HDR.EVENT.TYP(k1) = hex2dec('0441');
		elseif n==12,	% EMG right
			HDR.EVENT.TYP(k1) = hex2dec('0442');
		elseif n==13,	% EMG (foot)
			HDR.EVENT.TYP(k1) = hex2dec('0447');
		elseif n==1,	% Augen (left)
			HDR.EVENT.TYP(k1) = hex2dec('0431');
		elseif n==2,	% Augen (right)
			HDR.EVENT.TYP(k1) = hex2dec('0432');
		elseif n==3,	% Augen oben
			HDR.EVENT.TYP(k1) = hex2dec('0433');
		elseif n==4,	% Augen unten
			HDR.EVENT.TYP(k1) = hex2dec('0434');
		elseif n==5,	% blinzeln
			HDR.EVENT.TYP(k1) = hex2dec('0439');
		elseif n==6,	% Augen zu & entspannen
			HDR.EVENT.TYP(k1) = hex2dec('0115');
		elseif n==7,	% Augen offen & entspannen
			HDR.EVENT.TYP(k1) = hex2dec('0114');
		elseif n==8,	% beissen
			HDR.EVENT.TYP(k1) = hex2dec('0446');
		elseif n==9,	% kopf bewegen
			HDR.EVENT.TYP(k1) = hex2dec('0443');
		elseif any(n==[10,100])	% ende der aktion
			HDR.EVENT.TYP(k1) = bitxor(hex2dec('8000'),HDR.EVENT.TYP(k1-1));
		else
			HDR.EVENT.TYP(k1) = n;
		end;

% hits and misses, feedback
	elseif strncmp(tmp,'S',1) | strncmp(tmp,'R',1)
		HDR.EVENT.CHN(k1) = (tmp(1)=='R')*64+1;  %%% hack to distinguish Player 1 and 2 in SEASON2 data
		n = str2double(tmp(2:end));
		if n==11,	% hit (left)
		       	HDR.EVENT.TYP(k1) = hex2dec('0381');
		elseif n==12,	% hit (right)
	        	HDR.EVENT.TYP(k1) = hex2dec('0381');
		elseif n==13,	% hit (foot)
	        	HDR.EVENT.TYP(k1) = hex2dec('0381');
        	elseif n==21,	% miss (left)
	        	HDR.EVENT.TYP(k1) = hex2dec('0382');
		elseif n==22,	% miss (right)
	        	HDR.EVENT.TYP(k1) = hex2dec('0382');
		elseif n==23,	% time out
	        	HDR.EVENT.TYP(k1) = hex2dec('830d');
		elseif n==60,	% feedback onset
	        	HDR.EVENT.TYP(k1) = hex2dec('030d');
		elseif any(n==[4,5,7]) %%% ignore these
	        	HDR.EVENT.TYP(k1) = NaN;
        	else
	        	HDR.EVENT.TYP(k1) = n;
		end;

% end of segment
        elseif strcmp(tmp,'s') | strcmp(tmp,'stop') | strcmp(tmp,'stopp'),
        	HDR.EVENT.TYP(k1) = bitxor(hex2dec('8000'),HDR.EVENT.TYP(k1-1));

        elseif ~isempty(tmp)
        	[n,v,s] = str2double(tmp(2:end));
        	if (length(n)==1) & (~v)
        		HDR.EVENT.TYP(k1) = n;
       		end;
        end;
end;
HDR.EVENT.TYP = HDR.EVENT.TYP(:);

if isfield(HDR.EVENT,'POS');
	% remove S4, S5, S7 before computing trial end
	ix = find(HDR.EVENT.TYP(2:end)==hex2dec('7ffe'))+1;
	for k = 1:length(ix),
		if any(HDR.EVENT.TYP(ix(k)-1)==[1:3])
			HDR.EVENT.TYP(ix(k)-1) = NaN;	% if cue before end of segment is 1,2,or 3.
		end;
	end;

       	flag_remove = isnan(HDR.EVENT.TYP);
	HDR.EVENT.TYP = HDR.EVENT.TYP(~flag_remove);
	HDR.EVENT.POS = HDR.EVENT.POS(~flag_remove);
	HDR.EVENT.CHN = HDR.EVENT.CHN(~flag_remove);
	HDR.EVENT.DUR = HDR.EVENT.DUR(~flag_remove);
	if isfield(HDR.EVENT,'TeegType');
		HDR.EVENT.TeegType = HDR.EVENT.TeegType(~flag_remove);
	end;
	if isfield(HDR.EVENT,'Desc')
		HDR.EVENT.Desc = HDR.EVENT.Desc(~flag_remove);
	end;

       	ix1 = find(HDR.EVENT.TYP<10);
       	ix2 = find(HDR.EVENT.TYP==100);
       	if ~isempty(ix2),
		HDR.EVENT.TYP(ix2,1) = HDR.EVENT.TYP(ix2-1) + hex2dec('8000');
	end;
	ix0 = find((HDR.EVENT.TYP>0)&(HDR.EVENT.TYP<10));

	HDR.TRIG = HDR.EVENT.POS(ix0);
	HDR.Classlabel = HDR.EVENT.TYP(ix0);
end;

if any(HDR.EVENT.DUR~=1)
	warning('Duration is not 1')
end;

% convert from Type1 into Type3 table.
if 1, % ~isfield(HDR.EVENT,'CHN') & ~isfield(HDR.EVENT,'DUR'),
	% HDR.EVENT.CHN = zeros(size(HDR.EVENT.POS));
	HDR.EVENT.DUR = zeros(size(HDR.EVENT.POS));

	% convert EVENT.Version 1 to 3, currently used by GDF, BDF and alpha
	flag_remove = zeros(size(HDR.EVENT.TYP));
	types  = unique(HDR.EVENT.TYP);
	for k1 = find(bitand(types(:)',hex2dec('8000')));
	        TYP0 = bitand(types(k1),hex2dec('7fff'));
	        TYP1 = types(k1);
	        ix0  = (HDR.EVENT.TYP==TYP0);
	        ix1  = (HDR.EVENT.TYP==TYP1);

	        if sum(ix0)==sum(ix1),
	                HDR.EVENT.DUR(ix0) = HDR.EVENT.POS(ix1) - HDR.EVENT.POS(ix0);
	                flag_remove = flag_remove | (HDR.EVENT.TYP==TYP1);
                elseif 0,

                else
	                fprintf(2,'Warning BV2BIOSIG_EVENT: number of event onset (TYP=%s) and event offset (TYP=%s) differ (%i-%i) in %s\n',dec2hex(double(TYP0)),dec2hex(double(TYP1)),sum(ix0),sum(ix1),HDR.FileName);
                        %% double(.) operator needed because Matlab6.5 can not fix(uint16(..))
	        end;
	end;
	if any(HDR.EVENT.DUR<0)
	        fprintf(2,'Warning SOPEN: EVENT ONSET later than EVENT OFFSET\n',dec2hex(TYP0),dec2hex(TYP1));
	        %HDR.EVENT.DUR(:) = 0
	end;
	HDR.EVENT.TYP = HDR.EVENT.TYP(~flag_remove);
	HDR.EVENT.POS = HDR.EVENT.POS(~flag_remove);
	HDR.EVENT.CHN = HDR.EVENT.CHN(~flag_remove);
	HDR.EVENT.DUR = HDR.EVENT.DUR(~flag_remove);
	if isfield(HDR.EVENT,'TeegType');
		HDR.EVENT.TeegType = HDR.EVENT.TeegType(~flag_remove);
	end;
	%HDR.EVENT.TeegType = HDR.EVENT.TeegType(~flag_remove);
	if isfield(HDR.EVENT,'Desc')
		HDR.EVENT.Desc = HDR.EVENT.Desc(~flag_remove);
	end;
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function p = cycle_to_perm ( n, ncycle, t, index )

%*****************************************************************************80
%
%% CYCLE_TO_PERM converts a permutation from cycle to array form.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the number of items permuted.
%    N must be positive.
%
%    Input, integer NCYCLE, the number of cycles.
%    1 <= NCYCLE <= N.
%
%    Input, integer T(N), INDEX(NCYCLE), describes the permutation
%    as a collection of NCYCLE cycles.  The first cycle is
%    T(1) -> T(2) -> ... -> T(INDEX(1)) -> T(1).
%
%    Output, integer P(N), describes the permutation using a
%    single array.  For each index I, I -> P(I).
%

%
%  Check.
%
  ierror = cycle_check ( n, ncycle, t, index );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CYCLE_TO_PERM - Fatal error!\n' );
    fprintf ( 1, '  The input array is illegal.\n' );
    fprintf ( 1, '  IERROR = %d\n', ierror );
    error ( 'CYCLE_TO_PERM - Fatal error!' );
  end

  jhi = 0;

  for i = 1 : ncycle

    jlo = jhi + 1;
    jhi = jhi + index(i);

    for j = jlo : jhi

      if ( j < jhi )
        p(t(j)) = t(j+1);
      else
        p(t(j)) = t(jlo);
      end

    end

  end

  return
end


text:
% DEMOIL2 Oil data with fully independent training conditional, and MLP back constraints.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'oil';
experimentNo = 2;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('fitc');
options.optimiser = 'scg';
options.back = 'mlp';
options.backOptions = mlpOptions;
latentDim = 2;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
capName = dataSetName;;
capName(1) = upper(capName(1));
save(['dem' capName num2str(experimentNo) '.mat'], 'model');

if exist('printDiagram') & printDiagram
  fgplvmPrintPlot(model, lbls, capName, experimentNo);
end

% Load the results and display dynamically.
fgplvmResultsDynamic(dataSetName, experimentNo, 'vector')

% compute the nearest neighbours errors in latent space.
errors = fgplvmNearestNeighbour(model, lbls);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [model] = populateModelMetStr(model, metabolite_structure,replaceAllHits)
% This function populates the model structure with information contained in
% the metabolite structure
%
% INPUT
% model                     model structure
% metabolite_structure      metabolite structure
% replaceAllHits            replace existing entries (defaul: false)
%
% OUTPUT
% model                     updated model structure
%
%
% Ines Thiele 10/2021
if ~exist('replaceAllHits','var')
    replaceAllHits = 0;
end

if ~exist('replaceAllHits','var')
    replaceAllHits = 0;
end

F = fieldnames(metabolite_structure);
fields =fieldnames( metabolite_structure.(F{1}));
% load translation table between metabolite structure and COBRA model
% fields
translateMetStr2COBRAmodel;

for i = 1 : length(F)
    if isfield(metabolite_structure.(F{i}),'originalID')
        oriID =  metabolite_structure.(F{i}).originalID;
        % remove [] if present
        [tok,rem] = strtok(oriID,'[');
        % find metabolite in mets
        potHits = model.mets(contains(model.mets,tok));
    else
        % use VMHId
        potHits = model.mets(contains(model.mets,metabolite_structure.(F{i}).VMHId));
        tok = metabolite_structure.(F{i}).VMHId;
    end
    for j = 1 :length(potHits)
        [tokHit,rem] = strtok(potHits{j},'[');
        % check for perfect match now
        if  strcmp(tokHit,tok)
            hit = find(ismember(model.mets,potHits{j}));
            % now assign the new IDs
            for k = 1 : size(translation,1)

                if ~isempty(find(strcmp(fields,translation{k,1})))% check that translation field is part of F
                    % if a field is not present - add field
                    if ~isfield(model,translation{k,2})
                        model.(translation{k,2}) = [];
                    end
                    try % works  for non-numeric entries
                        if replaceAllHits == 1
                            model.(translation{k,2}){hit} = metabolite_structure.(F{i}).(translation{k,1});
                        elseif replaceAllHits ==0
                            if length(model.(translation{k,2}))>=hit
                                if isempty(model.(translation{k,2}){hit}) || length(find(isnan(model.(translation{k,2}){hit})))>0
                                    model.(translation{k,2}){hit} = metabolite_structure.(F{i}).(translation{k,1}) ;
                                end
                            else
                                model.(translation{k,2}){hit} = metabolite_structure.(F{i}).(translation{k,1}) ;

                            end
                        end
                    catch % works  for numeric entries
                        try
                            if replaceAllHits == 1
                                model.(translation{k,2})(hit) = str2num(metabolite_structure.(F{i}).(translation{k,1}));
                            elseif   replaceAllHits == 0
                                if length(model.(translation{k,2}))>=hit
                                    if isempty(model.(translation{k,2})(hit)) %|| length(find(isnan(model.(translation{k,2})(hit))))>0
                                        model.(translation{k,2})(hit) = str2num(metabolite_structure.(F{i}).(translation{k,1})) ;
                                    end
                                else
                                    model.(translation{k,2})(hit) = str2num(metabolite_structure.(F{i}).(translation{k,1})) ;
                                end
                            end
                        catch
                            if replaceAllHits == 1
                                model.(translation{k,2})(hit) = (metabolite_structure.(F{i}).(translation{k,1}));
                            elseif   replaceAllHits == 0
                                if length(model.(translation{k,2}))>=hit
                                    if isempty(model.(translation{k,2})(hit)) %|| length(find(isnan(model.(translation{k,2})(hit))))>0
                                        model.(translation{k,2})(hit) = (metabolite_structure.(F{i}).(translation{k,1})) ;
                                    end
                                else
                                    model.(translation{k,2})(hit) = (metabolite_structure.(F{i}).(translation{k,1})) ;
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end
% add SBO term to all metabolites SBO:0000247 represents the term 'simple chemical'.
for i = 1 : length(model.mets)
    model.metSBOTerms{i} = 'SBO:0000247';
end
for k = 1 : size(translation,1)
    if isfield(model,translation{k,2})
        model.(translation{k,2}) = columnVector(model.(translation{k,2}));
    end
end

% Append CHEBI: in front of CHEBI Ids
for i = 1 : length(model.mets)
    if ~isnan( model.metChEBIID{i})
        model.metChEBIID{i} = strcat('CHEBI:',model.metChEBIID{i});
    end
end

% ensure uniform output
if isfield(model,'metPubChemID')
    model.metPubChemID = cellfun(@num2str, model.metPubChemID,'UniformOutput',false);
end
if isfield(model,'metKEGGID')
    model.metKEGGID = cellfun(@num2str, model.metKEGGID,'UniformOutput',false);
end

if isfield(model,'metChEBIID')
    model.metChEBIID = cellfun(@num2str, model.metChEBIID,'UniformOutput',false);
end

if isfield(model,'metHMDBID')
    model.metHMDBID = cellfun(@num2str, model.metHMDBID,'UniformOutput',false);
end
if isfield(model,'metSEEDID')
    model.metSEEDID = cellfun(@num2str, model.metSEEDID,'UniformOutput',false);
end
if isfield(model,'metInchiKey')
    model.metInchiKey = cellfun(@num2str, model.metInchiKey,'UniformOutput',false);
end
if isfield(model,'metReactomeID')
    model.metReactomeID = cellfun(@num2str, model.metReactomeID,'UniformOutput',false);
end
if isfield(model,'metMetaNetXID')
    model.metMetaNetXID = cellfun(@num2str, model.metMetaNetXID,'UniformOutput',false);
end
if isfield(model,'metBioCycID')
    model.metBioCycID = cellfun(@num2str, model.metBioCycID,'UniformOutput',false);
end
if isfield(model,'metBiGGID')
    model.metBiGGID = cellfun(@num2str, model.metBiGGID,'UniformOutput',false);
end
if isfield(model,'metInchiString')
    model.metInchiString = cellfun(@num2str, model.metInchiString,'UniformOutput',false);
end

if isfield(model,'metLIPIDMAPSID')
    model.metLIPIDMAPSID = cellfun(@num2str, model.metLIPIDMAPSID,'UniformOutput',false);
    model.metLIPIDMAPSID  =  model.metLIPIDMAPSID (cellfun('isclass',  model.metLIPIDMAPSID , 'char'));
end

if isfield(model,'metSmile')
    model.metSmile = cellfun(@num2str, model.metSmile,'UniformOutput',false);
    model.metSmile  =  model.metSmile (cellfun('isclass',  model.metSmile , 'char'));
end

author: opencobra

repo: cobratoolbox


text:
function [ri,ti] = traintestindex( s, l1, l2, l3, l4)
% Copyright (c) 2014-present University of Illinois at Urbana-Champaign
% All rights reserved.
%
% Developed by:     Po-Sen Huang, Paris Smaragdis
%                   Department of Electrical and Computer Engineering
%                   Department of Computer Science
%
% Create a training/testing data set based ont he leave-out info provided

ri = [];
ti = [];
for i = 1:s(1)
	for j = 1:s(2)
		for l = 1:s(3)
			for k = 1:s(4)
				if (ismember( i, l1) || isempty( l1)) && (ismember( j, l2) || isempty( l2)) ...
						&& (ismember( l, l3) || isempty( l3)) && (ismember( k, l4) || isempty( l4))
					ti(end+1) = sub2ind( s, i, j, l, k);
				elseif ~ismember( i, l1) && ~ismember( j, l2) && ~ismember( l, l3) && ~ismember( k, l4)
					ri(end+1) = sub2ind( s, i, j, l, k);
				end
			end
		end
	end
end

author: posenhuang

repo: deeplearningsourceseparation


text:
function pde = SinFlat(am,ap,bm,bp,x0,y0,z0,n)
%% USAGE: polynomial solution for Poisson equation
%  Last Modified: 02/21/2020 by Xu Zhang

%% PDE Structure
pde = struct('intf',@intf,...
    'exactu1',@exactu1,'exactu2',@exactu2,'exactu3',@exactu3,...
    'um1',@um1,'um2',@um2,'um3',@um3,'up1',@up1,'up2',@up2,'up3',@up3,...
    'Dxu',@Dxu,'Dxum',@Dxum,'Dxup',@Dxup,'Dyu',@Dyu,...
    'Dyum',@Dyum,'Dyup',@Dyup,'Dzu',@Dzu,'Dzum',@Dzum,'Dzup',@Dzup,...
    'f1',@f1,'f2',@f2,'f3',@f3,...
    'fm1',@fm1,'fm2',@fm2,'fm3',@fm3,...
    'fp1',@fp1,'fp2',@fp2,'fp3',@fp3,...
    'A',@A,'Am',@Am,'Ap',@Ap,'one',@one,...
    'B',@B,'Bm',@Bm,'Bp',@Bp);

pde.bm = bm;
pde.bp = bp;
pde.am = am;
pde.ap = ap;
%% interface function
    function u = intf(x,y,z)
        u = n(1)*(x-x0)+n(2)*(y-y0)+n(3)*(z-z0);
    end

%% exact solution
    function u = exactu1(x,y,z)
        u = um1(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up1(x(id),y(id),z(id));
    end
    function u = exactu2(x,y,z)
        u = um2(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up2(x(id),y(id),z(id));
    end
    function u = exactu3(x,y,z)
        u = um3(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up3(x(id),y(id),z(id));
    end
    function u = um1(x,y,z)
        u = sin(pi*ker(x,y,z))/am;
    end
    function u = um2(x,y,z)
        u = sin(pi*ker(x,y,z))/am;
    end
    function u = um3(x,y,z)
        u = sin(pi*ker(x,y,z))/am;
    end
    function u = up1(x,y,z)
        u = sin(pi*ker(x,y,z))/ap;
    end
    function u = up2(x,y,z)
        u = sin(pi*ker(x,y,z))/ap;
    end
    function u = up3(x,y,z)
        u = sin(pi*ker(x,y,z))/ap;
    end
    function u = ker(x,y,z)
        u = n(1)*(x-x0)+n(2)*(y-y0)+n(3)*(z-z0);
    end

%% Boundary Function
    function u = gD1(x,y,z)
        u = exactu1(x,y,z);
    end
    function u = gD2(x,y,z)
        u = exactu2(x,y,z);
    end
    function u = gD3(x,y,z)
        u = exactu3(x,y,z);
    end
%% Derivative of the exact solution
    function u = Dxu(x,y,z)
        u = Dxum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dxup(x(id),y(id),z(id));
    end
    function u = Dyu(x,y,z)
        u = Dyum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dyup(x(id),y(id),z(id));
    end
    function u = Dzu(x,y,z)
        u = Dzum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dzup(x(id),y(id),z(id));
    end
    function u = Dxum(x,y,z)
        u = cos(pi*ker(x,y,z))*pi*(n(2)-n(3))/am;
    end
    function u = Dyum(x,y,z)
        u = cos(pi*ker(x,y,z))*pi*(n(3)-n(1))/am;
    end
    function u = Dzum(x,y,z)
        u = cos(pi*ker(x,y,z))*pi*(n(1)-n(2))/am;
    end
    function u = Dxup(x,y,z)
        u = cos(pi*ker(x,y,z))*pi*(n(2)-n(3))/ap;
    end
    function u = Dyup(x,y,z)
        u = cos(pi*ker(x,y,z))*pi*(n(3)-n(1))/ap;
    end
    function u = Dzup(x,y,z)
        u = cos(pi*ker(x,y,z))*pi*(n(1)-n(2))/ap;
    end

%% right hand side function
    function u = f1(x,y,z)
        u = fm1(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp1(x(id),y(id),z(id));
    end
    function u = f2(x,y,z)
        u = fm2(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp2(x(id),y(id),z(id));
    end
    function u = f3(x,y,z)
        u = fm3(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp3(x(id),y(id),z(id));
    end

    function u = fm1(x,y,z)
        u = -sin(pi*ker(x,y,z))*pi^2*(n(1)*(n(2)+n(3))-n(2)^2-n(3)^2)+um1(x,y,z)*bm;
    end
    function u = fm2(x,y,z)
        u = -sin(pi*ker(x,y,z))*pi^2*(n(2)*(n(1)+n(3))-n(1)^2-n(3)^2)+um2(x,y,z)*bm;
    end
    function u = fm3(x,y,z)
        u = -sin(pi*ker(x,y,z))*pi^2*(n(3)*(n(1)+n(2))-n(1)^2-n(2)^2)+um3(x,y,z)*bm;
    end
    function u = fp1(x,y,z)
        u = -sin(pi*ker(x,y,z))*pi^2*(n(1)*(n(2)+n(3))-n(2)^2-n(3)^2)+up1(x,y,z)*bp;
    end
    function u = fp2(x,y,z)
        u = -sin(pi*ker(x,y,z))*pi^2*(n(2)*(n(1)+n(3))-n(1)^2-n(3)^2)+up2(x,y,z)*bp;
    end
    function u = fp3(x,y,z)
        u = -sin(pi*ker(x,y,z))*pi^2*(n(3)*(n(1)+n(2))-n(1)^2-n(2)^2)+up2(x,y,z)*bp;
    end

%% Diffusion coefficient function
    function u = A(x,y,z)
        u = Am(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Ap(x(id),y(id),z(id));
    end
    function u = Am(x,y,z)
        u = am*ones(size(x));
    end
    function u = Ap(x,y,z)
        u = ap*ones(size(x));
    end
%% Mass coefficient function
    function u = B(x,y,z)
        u = Bm(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Bp(x(id),y(id),z(id));
    end
    function u = Bm(x,y,z)
        u = bm*ones(size(x));
    end
    function u = Bp(x,y,z)
        u = bp*ones(size(x));
    end

%% Other function
    function u = one(x,y,z)
        u = ones(size(x));
    end
end


text:
function varargout= bbci_acquire_cognionics(varargin)
%BBCI_ACQUIRE_cognionics - Online data acquisition from Cognionics headset
%
%Synopsis:
%  STATE= bbci_acquire_cognionics('init', <PARAM>)
%  [CNTX, MRKTIME, MRKDESC, STATE]= bbci_acquire_cognionics(STATE)
%  bbci_acquire_cognionics('close')
%  bbci_acquire_cognionics('close', STATE)
%
%Arguments:
%  PARAM - Optional arguments
%    'clab', 'fs', 'filtHd', 'blocksize', 'port', 'timeout', 'verbose'
%
%Output:
%  STATE - Structure characterizing the incoming signals; fields:
%     'fs', 'clab', and intern stuff
%  CNTX - 'acquired' signals [Time x Channels]
%    The following variables hold the markers that have been 'acquired' within
%    the current block (if any).
%  MRKTIME - DOUBLE: [1 nMarkers] position [msec] within data block.
%      A marker occurrence within the first sample would give
%      MRKTIME= 1/STATE.fs.
%  MRKDESC - DOUBLE [1 nMarkers] values of markers

% benjamin.blankertz

if isequal(varargin{1}, 'init'),
  state= opt_proplistToStruct(varargin{2:end});
  default_clab= ...
    {'AF5' 'AF3' 'AF1' 'AFz' 'AF2' 'AF4' 'AF6' ...
     'F5' 'F3', 'F1' 'Fz' 'F2' 'F4' 'F6' ...
     'FC7' 'FC5' 'FC3' 'FC1' 'FCz' 'FC2' 'FC4' 'FC6' 'FC8' ...
     'T7' 'C5' 'C3' 'C1' 'Cz' 'C2' 'C4' 'C6' 'T8' ...
     'CP7' 'CP5' 'CP3' 'CP1' 'CPz' 'CP2' 'CP4' 'CP6' 'CP8' ...
     'P7' 'P5' 'P3' 'P1' 'Pz' 'P2' 'P4' 'P6' 'P8' ...
     'PO5' 'PO3' 'PO1' 'POz' 'PO2' 'PO4' 'PO6' ...
     'O5' 'O3' 'O1' 'Oz' 'O2' 'O4' 'O6' ...
    };
  props= {'fs'            250            '!DOUBLE[1]'
          'clab'          default_clab   'CELL{CHAR}'
          'blocksize'     40             '!DOUBLE[1]'
          'port'          'COM11'        '!CHAR'
          'timeout'       3              '!DOUBLE[1]'
          'filtHd'        []             'STRUCT'
          'verbose'       true           '!BOOL'
         };
  [state, isdefault]= opt_setDefaults(state, props, 1);
  if state.fs~=250,
    error('Only fs=250 allowed');
  end
  if isdefault.filtHd,
    % Fs/4
    filt1.b= [0.85 0 0.85]
    filt1.a= [1 0 0.7]
    %Fs/2
    filt2.b= [0.8 0.8]
    filt2.a= [1 0.6]
    state.filtHd= procutil_catFilters(filt1, filt2);
    state.filtHd.PersistentMemory= true;
  end
  state.nChans= length(state.clab);
  state.nBytesPerPacket= 2+3*state.nChans+4;
  nPacketsPerPoll= ceil(state.blocksize/1000*state.fs);
  state.nBytesPerPoll= nPacketsPerPoll*state.nBytesPerPacket;

  % try to find connected cognionics port
  cog_ports= instrfind('Tag', 'Cognionics');
  if isempty(cog_ports),
    state.sp= serial(state.port, 'BaudRate',115200, ...
                                 'Timeout',state.timeout, ...
                                 'Tag', 'Cognionics');
  elseif length(cog_ports)==1,
    state.sp= cog_ports;
    if strcmp(state.sp.Status, 'open'),
      fclose(state.sp);
    end
  else
    error('multiple ports connected to Cognioncs found');
  end
  state.sp.InputBufferSize= 5*state.nBytesPerPoll;

  fopen(state.sp);
  % Poll data several times to get into a more stable state
  for q= 1:3,
    fread(state.sp, state.sp.InputBufferSize, 'uint8');
  end
  state.packetNo=[];
  state.buffer= [];
  state.lastx= zeros(1, state.nChans);
  state.lastMrkDesc= 256;
  state.scale= 1000000/2^24;
  if isempty(state.filtHd),
    reset(state.filtHd);
  end
  output= {state};

elseif isequal(varargin{1}, 'close'),
  cog_ports= instrfind('Tag', 'Cognionics');
  if ~isempty(cog_ports),
    fclose(cog_ports);
    %delete(cog_ports);  % shutdown port completely?
  end
  output= {};
elseif length(varargin)~=1,
  error('Except for INIT/CLOSE case, only one input argument expected');

else
  if ~isstruct(varargin{1}),
    error('First input argument must be ''init'', ''close'', or a struct');
  end
  state= varargin{1};
  %while state.sp.BytesAvailable<state.nBytesPerPoll,
  %  pause(0.01);
  %end

  % Receive data packet from COM port and check packet IDs
  chunk= fread(state.sp, state.nBytesPerPoll, 'uint8');
  chunk= cat(1, state.buffer, chunk);
  chunk= bbciutil_cognionicsRemoveIncompletePackets(chunk, state);
  idx= find(chunk==255);
  % IDs at the very end are not considered, because packet # is not available
  if ~isempty(idx) && idx(end)==length(chunk),
    idx(end)= [];
  end

  % if last packet is incomplete, put it to the buffer
  % and exclude that packet from current acquisition
  if length(chunk)-idx(end)+1 < state.nBytesPerPacket,
    state.buffer= chunk(idx(end):end);
    idx(end)= [];
  else
    state.buffer= [];
  end
  if isempty(idx),
    if state.verbose,
      fprintf('[ACQ-COG] Did not receive any complete data block.\n');
    end
    state.buffer= chunk;
    output= {[], [], [], state};
    varargout= output(1:nargout);
    return;
  end

  % Extract packet numbers
  if isempty(state.packetNo),   % init case
    state.packetNo= chunk(idx(1)+1)-1;
  end
  packet_counter= [state.packetNo; chunk(idx+1)];
%  fprintf('buffer: #%d, packets: %s\n', state.packetNo, ...
%          str_vec2str(packet_counter(2:end)));
  packet_counter= mod(packet_counter,128);
  unwarp_correction= 128 * (floor(state.packetNo/128) + ...
                            [0; cumsum(diff(packet_counter)<0)]);
  packet_counter= packet_counter + unwarp_correction;
  state.packetNo= packet_counter(end);
%  fprintf('packets unwrapped: %s\n', str_vec2str(packet_counter));

  % Obtain signals from data packets
  np= length(idx);
  idxData= repmat(2+[1:3:3*state.nChans]', [1 np]) + ...
           repmat([0:np-1] * state.nBytesPerPacket, [state.nChans 1]);
  idxData= reshape(idxData, [1 np*state.nChans]);
  msb= chunk(idxData);
  x= bitshift(uint32(msb), 16) + ...
     bitshift(uint32(chunk(idxData+1)), 9) + ...
     bitshift(uint32(chunk(idxData+2)), 2);
  x= x + uint32(bitget(msb, 8)) * bitshift(uint32(255), 24);
  cntx= state.scale * double( typecast(x, 'int32') );
  cntx= reshape(cntx, [state.nChans np])';

  % Obtain markers from data packets
  idxTrg= idx(1) + 2 + 3*state.nChans + 1 + state.nBytesPerPacket*[0:np-1];
  mrkDesc= double( bitshift(uint16(chunk(idxTrg)), 8) + ...
                            uint16(chunk(idxTrg+1)) )';
  iNonvoid= find(mrkDesc & diff([state.lastMrkDesc mrkDesc]));
  state.lastMrkDesc= mrkDesc(end);
  mrkTime= 1/state.fs * iNonvoid;
  if state.verbose && ~isempty(iNonvoid),
    fprintf('[ACQ-COG] Markers: %s\n', str_vec2str(mrkDesc(iNonvoid)));
%     for m= 1:length(iNonvoid),
%       fprintf('[ACQ-COG] %6.3f: %s\n', mrkTime(m), mrkDesc(iNonvoid(m)));
%		 end
  end

  % Interpolate lost packets
  dpc= diff(packet_counter);
  ilost= find(dpc>1);
  if ~isempty(ilost),
    if state.verbose,
      fprintf('[ACQ-COG] interpolating %d lost packets after packet #%d\n', ...
              [dpc(ilost)-1, packet_counter(ilost)]');
    end
    [cntx, mapping]= ...
        bbciutil_cognionicsInterpolate(cntx, packet_counter, state);
    mrkTime= 1/state.fs * mapping(iNonvoid)';
  end
  state.lastx= cntx(end,:);

  % Apply filter if requested (dfilt.filter automatically saves the state)
  if ~isempty(state.filtHd),
    cntx= filter(state.filtHd, cntx, 1);
  end

  output= {cntx, mrkTime, mrkDesc(iNonvoid), state};
end
varargout= output(1:nargout);

author: bbci


text:
function [cc] = mi32cc(mi3)
% Convert volume from cubic miles to cubic centimeters.
% Chad Greene 2012
cc = mi3*4168181825400000 ;


text:
%ISHOMOG2 Test if SE(2) homogeneous transformation
%
% ISHOMOG2(T) is true (1) if the argument T is of dimension 3x3 or 3x3xN, else
% false (0).
%
% ISHOMOG2(T, 'valid') as above, but also checks the validity of the rotation
% sub-matrix.
%
% Notes::
% - The first form is a fast, but incomplete, test for a transform in SE(3).
% - Does not work for the SE(3) case.
%
% See also ISHOMOG, ISROT2, ISVEC.

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

function h = ishomog2(tr, rtest)
    d = size(tr);
    if ndims(tr) >= 2
        h =  all(d(1:2) == [3 3]);

        if h && nargin > 1
            h = abs(det(tr(1:2,1:2)) - 1) < eps;
        end
    else
        h = false;
    end

author: Allopart

repo: rbpf-gmapping


text:
function [Zm]=SVD_filter(Z,P,T)

%Computer SVD
[U,S,V] = svd(Z);

%Normalising singular values
Ss=diag(S)-min(diag(S));
Ss=Ss./max(Ss);

%Creating smoothening parameters
p_max=P(1); %1=No blurring
p_min=P(2); %0=straight line fit
p=(Ss.*(p_max-p_min))+p_min; %Scale towards singular values

Zm=nan(size(Z));
for i=1:1:size(U,2)
    v=V(:,i); u=U(:,i); s=S(i,i); %components
    if Ss(i)<=T %Filter after threshold
        us = csaps(1:numel(u),u,p(i),1:numel(u))'; %Smooth u
        vs = csaps(1:numel(v),v,p(i),1:numel(v))'; %Smooth v
    else
        vs=v; us=u; %Keep unsmoothened
    end
    z=us*s*vs'; %sub-data
    Zm(:,:,i)=z;
end
Zm=sum(Zm,3);

%%
% _*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:
classdef EdgeCutInElemComputer < handle

    properties (Access = private)
       nElem
       nEdgeByElem
    end

    properties (Access = private)
      edgesInElem
      isEdgeCut
    end

    methods (Access = public)

        function obj = EdgeCutInElemComputer(cParams)
            obj.init(cParams)
        end

        function isEdgeCut = compute(obj)
            isEdgeCut = false(obj.nEdgeByElem,obj.nElem);
            for iedge = 1:obj.nEdgeByElem
                edge = obj.edgesInElem(:,iedge);
                isEdgeCut(iedge,:) = obj.isEdgeCut(edge);
            end
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.edgesInElem = cParams.edgesInElem;
            obj.isEdgeCut   = cParams.isEdgeCut;
            obj.nElem       = size(obj.edgesInElem,1);
            obj.nEdgeByElem = size(obj.edgesInElem,2);
        end

    end

end


text:
function varargout = symmetrise(m,varargin)
% symmetrcially equivalent directions and its multiple
%
% Syntax
%   mSym = symmetrise(m)
%
%   % include antipodal symmetry
%   mSym = symmetrise(m,'antipodal')
%
%   % exclude antipodal symmetry
%   mSym = symmetrise(m,'noAntipodal')
%
%   % every symmetrically equivalent direction only once
%   [mSym,l,sym] = symmetrise(m,'unique')
%
%   % every symmetrically equivalent axis only once
%   [mSym,l,sym] = symmetrise(v,S,'unique','noAntipodal')
%
% Input
%  v - @Miller
%
% Output
%  mSym - sym * m  @Miller
%  l    - multiplicity of the crystal directions
%  sym  - @rotation
%
% Flags
%  antipodal   - include <VectorsAxes.html antipodal symmetry>
%  noAntipodal - do not include antipodal symmetry (without option unique)
%  noAntipodal - do not remove antipodal vectors (with option unique)
%  unique      - only return distinct axes or directions (noAntipodal)
%

[varargout{1:nargout}] = symmetrise@vector3d(m,m.CS,varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
%File Description:
%  Unknown parameters in accelerometer error model are calculated using the
%  lm optimization algorithm.
close all
clc
clear

load AccRaw  % load uncalibrated accelerometer data
g = 9.8;
m = length(AccRaw);

y_dat = g*ones(m, 1);  % expected gravitational acceleration
p0 = [1 1 1 0 0 0]';
p_init = [1.0 1.0 1.0 0.1 0.1 0.1]';  % initial value of unknown parameters

y_raw = calFunc(AccRaw, p0);  % gravitational acceleration measured by accelerometer
y_raw = y_raw(:);

r_raw = y_dat - y_raw;
p_fit = lm('calFunc', p_init, AccRaw, y_dat);
y_lm = calFunc(AccRaw, p_fit);  % gravitational acceleration measured by calibrated accelerometer
y_lm = y_lm(:);
r_lm = y_dat - y_lm;

kx = p_fit(1);
ky = p_fit(2);
kz = p_fit(3);
bx = p_fit(4);
by = p_fit(5);
bz = p_fit(6);

Ka9_8 = [kx 0 0; 0 ky 0; 0 0 kz]
ba9_8 = [bx by bz]'
save('calP9_8', 'Ka9_8', 'ba9_8')

figure
bar([r_raw'*r_raw, r_lm'*r_lm])
grid on;
set(gca, 'XTickLabel', {'raw', 'lm'});
ylabel('fa');

t=1:m;
figure
title('Accelerometer Calibration')
plot(t, r_raw, t, r_lm)
legend('Uncalibrated', 'Calibrated-LM')
xlabel('Sampling sequence')
ylabel('Residual')

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
function [vert,face] = MakeZcylinder(pos, radius,len)

a = 10;    % number of side faces
theta = (0:a-1)/a * 2*pi;

x  = radius*cos(theta);
y  = radius*sin(theta);
z1 = len/2 * ones(1,a);
z2 = -z1;

vert    = [x x 0 0;
           y y 0 0;
           z1 z2 len/2 -len/2];
for n = 1:3
   vert(n,:) = vert(n,:) + pos(n);
end

% make index data specifying polygons
face_side = [1:a; a+1:2*a; a+2:2*a a+1; 2:a 1];
face_up   = [1:a; 2:a 1];
face_up(3:4,:) = 2*a+1;  % index of up center
face_down = [a+2:2*a a+1; a+1:2*a];
face_down(3:4,:) = 2*a+2; % index of down center
face = [face_side face_up face_down];

author: s-kajita

repo: IntroductionToHumanoidRobotics


text:
% BOOSTING_EXAMPLE
%
% Constructs a 2-dimensional dataset classifiable by boosting, but not any
% simple linear classifier, because of the thresholding nature of the data.

rand('seed', 0);

% m datapoints in 2-dimensions
mm = 150;
X = rand(mm, 2);

thresh_pos = .6;
y = [X(:, 1) < thresh_pos & X(:, 2) < thresh_pos];
y = 2 * y - 1;

for T = [2, 4, 5, 10]
  figure;
  hpos = plot(X(y == 1, 1), X(y == 1, 2), 'o');
  hold on;
  hneg = plot(X(y == -1, 1), X(y == -1, 2), 'x');
  set(hpos, 'linewidth', 2);
  set(hneg, 'linewidth', 2);

  [theta, feature_inds, thresholds] = stump_booster(X, y, T);

  x1_coords = linspace(0, 1, 100);
  x2_coords = linspace(0, 1, 100);
  Z = zeros(100);
  for ii = 1:100
    for jj = 1:100
      pred = (sign(x1_coords(ii) - thresholds(feature_inds == 1))' * ...
              theta(feature_inds == 1)) + ...
             (sign(x2_coords(jj) - thresholds(feature_inds == 2))' * ...
              theta(feature_inds == 2));
      Z(jj, ii) = sign(pred);
    end
  end

  C = contourc(x1_coords, x2_coords, Z, [0 0]);
  h = plot(C(1, 2:end), C(2, 2:end), 'k-');
  set(h, 'linewidth', 2);
  title(sprintf('Iterations = %d', T));
  set(gca, 'fontsize', 18);
  print('-depsc2', sprintf('boost_plot_%d.eps', T));
end

%% Now solve the logistic regression problem directly

mm = 200;
X = rand(mm, 2);
y = [X(:, 1) < thresh_pos & X(:, 2) < thresh_pos];
y = 2 * y - 1;

theta_log = zeros(3, 1);
X_logit = [ones(mm, 1), X];
for iter = 1:1000
  risk = (1/mm) * sum(log(1 + exp(-y .* (X_logit * theta_log))));
  if (mod(iter, 50) == 0)
    fprintf(1, 'Iter %d, loss %1.4f\n', iter, risk);
  end
  p = 1 ./ (1 + exp(y .* (X_logit * theta_log)));
  g = -(1/mm) * X_logit' * (p .* y);
  theta_log = theta_log - 2 * g;
end

x1_coord = linspace(0, 1, 100);
x2_coord = -(theta_log(1) + theta_log(2) * x1_coord) / theta_log(3);

figure;
hpos = plot(X(y == 1, 1), X(y == 1, 2), 'o');
hold on;
hneg = plot(X(y == -1, 1), X(y == -1, 2), 'x');
set(hpos, 'linewidth', 2);
set(hneg, 'linewidth', 2);
h = plot(x1_coord, x2_coord, 'k-', 'linewidth', 2);
axis([0 1 0 1]);
set(gca, 'fontsize', 18);
print -depsc2 'logistic_plot.eps';

author: cycleuser


text:
function pass = test_inv(pref)
% This test constructs two CHEBFUN objects and uses INV() to invert them. It
% checks the inverse calculated is accurate.

% Taken from Chebfun v4 test, invtest.m, by Nick Hale  07/06/2009.

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

algoList = {'roots', 'newton', 'bisection', 'regulafalsi', 'illinois', 'brent'};

for k = 1:6
    x = chebfun('x');
    f = sin(x);
    g = chebfun(@(x) asin(x), [sin(-1), sin(1)]);
    f_inv = inv(f, pref, 'algorithm', algoList{k});
    tol = 100*eps*vscale(f_inv);
    pass(k,1) = norm(g - f_inv, inf) < tol;

    g = inv(f_inv, 'algorithm', algoList{k});
    [f, g] = tweakDomain(f, g, 2e-15);
    pass(k,2) = norm(f - g, inf) < tol;

    x = chebfun('x');
    f = chebfun(@(x) sausagemap(x));
    f_inv = inv(f, pref, 'algorithm', algoList{k});
    tol = 100*eps*vscale(f_inv);
    pass(k,3) = norm(f(f_inv) - x, inf) + norm(f_inv(f) - x, inf) < tol;

    % Check the 'monocheck' and 'rangecheck' options.
    f = chebfun(@exp);
    f_inv = inv(f, pref, 'algorithm', algoList{k}, 'monocheck', 'on', ...
        'rangecheck', 'on');
    tol = 100*eps*vscale(f_inv);
    xx = linspace(f_inv.domain(1), f_inv.domain(end), 20);
    pass(k,4) = norm(f(f_inv(xx)) - xx, inf) < tol;
    pass(k,5) = all(abs(f_inv.domain([1, end]) - exp([-1 1])) < 10*eps);
end

% Check that 'monocheck' fails for a non-monotonic function.
try
    f = chebfun(@(x) x.^2);
    f_inv = inv(f, pref, 'splitting', 'on', 'monocheck', 'on');
    pass(:,6) = false;
catch ME
    if ( strcmpi(ME.identifier, ...
            'CHEBFUN:CHEBFUN:inv:doMonoCheck:notMonotonic') )
        pass(:,6) = true;
    else
        pass(:,6) = false;
    end
end

% Test the example from the help text:
x = chebfun('x');
f = x + .5*abs(x) + .6*sign(x-.5);
g = inv(f);
xx = linspace(f.domain(1), f.domain(end), 10);
pass(:,7) = norm(g(f(xx)) - xx, inf) < 10*vscale(g)*eps;

% Test the inverse of a decreasing function (see #1098).
x = chebfun('x');
f = -sin(x);
g = inv(f);
xx = linspace(f.domain(1), f.domain(end), 10);
pass(:,8) = norm(g(f(xx)) - xx, inf) < 10*vscale(g)*eps;

end

function g = sausagemap(s,d)
if ( nargin < 2 )
    d = 9; % This can be adjusted
end
c = zeros(1,d+1);
c(d:-2:1) = [1 cumprod(1:2:d-2)./cumprod(2:2:d-1)]./(1:2:d);
c = c/sum(c); g = polyval(c,s);
cp = c(1:d).*(d:-1:1);
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function integt
%INTEGT	Unit test for the function INTEG.

%	O. Lemoine - March 1996.

N=128;

% Integration of a sinusoid on a period
x=(0:N);
y=sin(2*pi*x/N);
som=integ(y,x);
if abs(som)>sqrt(eps),
 error('integ test 1 failed');
end

% Integration of a constant
y=ones(1,N);
som=integ(y);
if abs(som-N+1)>sqrt(eps),
 error('integ test 2 failed');
end

% Integration of a ramp
y=(0:N);
som=integ(y);
if abs(som-N*N/2)>sqrt(eps),
 error('integ test 3 failed');
end

N=113;

% Integration of a sinusoid on a period
x=(0:N);
y=sin(2*pi*x/N);
som=integ(y,x);
if abs(som)>sqrt(eps),
 error('integ test 4 failed');
end

% Integration of a constant
y=ones(1,N);
som=integ(y);
if abs(som-N+1)>sqrt(eps),
 error('integ test 5 failed');
end

% Integration of a ramp
y=(0:N);
som=integ(y);
if abs(som-N*N/2)>sqrt(eps),
 error('integ test 6 failed');
end

author: HeLiangHIT


text:
function s = volumeVisualization(x,y,z,v)
%volumeVisualization Engine for choosing y-z planes to view.
%   s = volumeVisualization(X,Y,Z,V) returns a structure
%   containing information about the visualization of volume
%   data described by X,Y,Z,V.  The fields are:
%          addSlicePlane -- function handle to add a slice
%                           plane at location x
%   deleteLastSlicePLane -- function handle to delete the
%                           last slice plane added
%                   xMin -- minimum x location for a plane
%                   xMax -- maximum x location for a plane
%
%      Example:
%      [x,y,z,v] = flow;
%      s = volumeVisualization(x,y,z,v);
%      s.addSlicePlane(3.7)
%      s.addSlicePlane(7.5)
%      pause
%      s.deleteLastSlicePlane()
%      pause
%      s.deleteLastSlicePlane()

%   Copyright 2007 The MathWorks, Inc.

%% Store handles to the various planes
%initialize handle to axis
%initialize handle to slice plane
hAxis = [];
hSlicePlanes = [];

%% Create data for generic slice through yz-plane
[yd,zd] = meshgrid(linspace(min(y(:)),max(y(:)),100), ...
    linspace(min(z(:)),max(z(:)),100));

%% Plot the volume initially
initDisplay()

%% Nested Functions
    function addSlicePlane(xLoc)
    %addSlicePlane   Add a slice plane xLoc.
        xd            = xLoc*ones(size(yd));
        newSlicePlane = slice(hAxis, x, y, z, v, xd, yd, zd);
        hSlicePlanes   = [ hSlicePlanes, newSlicePlane ];
        set(newSlicePlane,'FaceColor'      ,'interp',...
                          'EdgeColor'      ,'none'  ,...
                          'DiffuseStrength',.8       );
    end

    function deleteLastSlicePlane()
    %deleteLastSlicePlane Delete the last slice plane added.
        if ~isempty(hSlicePlanes)
            delete(hSlicePlanes(end));
            hSlicePlanes = hSlicePlanes(1:end-1);
        end
    end

    function initDisplay()
    %initDisplay  Initialize Display.

        % Draw back and bottom walls
        if isempty(hAxis) || ~ishandle(hAxis)
            hAxis = gca;
            hold on;
        end
        hx = slice(hAxis, x, y, z, v, ...
            max(x(:)),       [],       []) ;
        hy = slice(hAxis, x, y, z, v, ...
            [],       max(y(:)),       []) ;
        hz = slice(hAxis, x, y, z, v, ...
            [],              [],min(z(:))) ;

        % Make everything look nice
        set([hx hy hz],'FaceColor','interp',...
            'EdgeColor','none')
        set(hAxis,'FontSize',18,'FontWeight','Bold');
        xlabel('X');ylabel('Y');zlabel('Z')
        daspect([1,1,1])
        axis tight
        box on
        view(-38.5,16)
        colormap (jet(128))
    end

s.addSlicePlane = @addSlicePlane;
s.deleteLastSlicePlane = @deleteLastSlicePlane;
s.xMin = min(x(:));
s.xMax = max(x(:));

end

author: Sable

repo: mcbench-benchmarks


text:
function dip = fixdipole(dip)

% FIXDIPOLE ensures that the dipole position and moment are
% consistently represented throughout FieldTrip functions.

% Copyright (C) 2009, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: fixdipole.m 2885 2011-02-16 09:41:58Z roboos $

[m, n] = size(dip.pos);

if n==3
  % the input representation is Nx3, which is what we want
elseif m==3
  % it is possible to translate it into a Nx3 unambiguously
  warning('input dipole positions should be specified as Nx3 matrix');
  dip.pos = dip.pos';
elseif m==1
  % it is possible to translate it into a Nx3 unambiguously
  warning('input dipole positions should be specified as Nx3 matrix');
  dip.pos = reshape(dip.pos, 3, n/3)';
else
  % it is not clear how to convert to a Nx3 matrix
  error('input dipole positions should be specified as Nx3 matrix');
end

if isfield(dip, 'mom')
  ndip = size(dip.pos,1);
  if numel(dip.mom)==ndip*3
    ntime = 1;
  else
    ntime = numel(dip.mom)/(ndip*3);
  end
  dip.mom = reshape(dip.mom, ndip*3, ntime);
end

author: PatternRecognition

repo: OpenBMI