text:
% Converts a MatLab array into a java.util.ArrayList.
%
% Input arguments:
% array:
%    a MatLab row or column vector (with elements of any type)
%
% Output arguments:
% list:
%    a java.util.ArrayList instance
%
% See also: javaArray

% Copyright 2010 Levente Hunyadi
function list = javaArrayList(array)

assert(isvector(array), 'javaArrayList:DimensionMismatch', ...
	'Row or column vector expected.');
list = java.util.ArrayList;
if iscell(array)  % convert cell array into ArrayList
	for k = 1 : numel(array)
		list.add(array{k});
	end
else  % convert (numeric) array into ArrayList
	for k = 1 : numel(array)
		list.add(array(k));
	end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function hermite_polynomial_test14 ( )

%*****************************************************************************80
%
%% HERMITE_POLYNOMIAL_TEST14 tests H_POLYNOMIAL_COEFFICIENTS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2012
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_POLYNOMIAL_TEST14\n' );
  fprintf ( 1, '  H_POLYNOMIAL_COEFFICIENTS determines the physicist''s Hermite \n' );
  fprintf ( 1, '  polynomial coefficients.\n' );

  c = h_polynomial_coefficients ( n );

  for i = 0 : n
    fprintf ( 1, '\n' );
    fprintf ( 1, '  H(%d) = \n', i );
    fprintf ( 1, '\n' );
    for j = i : -1 : 0
      if ( c(i+1,j+1) == 0.0 )

      elseif ( j == 0 )
        fprintf ( 1, '  %f\n', c(i+1,j+1) );
      elseif ( j == 1 )
        fprintf ( 1, '  %f * x\n', c(i+1,j+1) );
      else
        fprintf ( 1, '  %f * x^%d\n', c(i+1,j+1), j );
      end
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [c] = mstep_update_c(P, S_bar, V, E_z, E_zz, RO, Tr)
%MSTEP_UPDATE_C Summary of this function goes here
%   Detailed explanation goes here

%equation 48 of PAMI paper

[K, T] = size(E_z);
J = size(S_bar, 2);
c = zeros(T,1);
Pc = P - Tr(:)*ones(1,J);

parfor t=1:T

    zz_hat_t = [1 E_z(:,t)'; E_z(:,t) E_zz((t-1)*K+1:t*K,:)];
    z_hat_t = [1;E_z(:,t)];
    num = 0;
    den = 0;
    for j=1:J
        M_jt = [1 0 0 ; 0 1 0]*RO{t}*[S_bar(:,j) reshape(V(:,j), 3, K)];
        num = num + z_hat_t'*M_jt'*[Pc(t,j); Pc(t+T,j)];
        den = den+trace(M_jt*zz_hat_t*M_jt');
    end
    c(t,1) = num/den;
end

end

repo: CategoryShapes


text:
function string_maker(block)
% Level-2 MATLAB file S-Function for times two demo.
%   Copyright 1990-2009 The MathWorks, Inc.
%   $Revision: 1.1.6.2 $

  setup(block);

%endfunction

function setup(block)

  %% Register number of input and output ports
  block.NumDialogPrms  = 3;
  block.NumInputPorts  = 3;
  block.NumOutputPorts = 1;

  %% Setup functional port properties to dynamically
  %% inherited.
  block.SetPreCompInpPortInfoToDynamic;
  block.SetPreCompOutPortInfoToDynamic;

  block.InputPort(1).DirectFeedthrough = true;
  block.InputPort(1).DimensionsMode   = 'Fixed';
  block.InputPort(1).Dimensions       = 1;

  block.InputPort(2).DirectFeedthrough = true;
  block.InputPort(2).DimensionsMode   = 'Fixed';
  block.InputPort(2).Dimensions       = 1;

  block.InputPort(3).DirectFeedthrough = true;
  block.InputPort(3).DimensionsMode   = 'Fixed';
  block.InputPort(3).Dimensions       = 1;

  block.OutputPort(1).DimensionsMode   = 'Fixed';
  block.OutputPort(1).Dimensions       = 15;
  block.OutputPort(1).DatatypeID       = -1;
  %% Set block sample time to inherited
  block.SampleTimes = [-1 0];

  %% Set the block simStateCompliance to default (i.e., same as a built-in block)
  block.SimStateCompliance = 'DefaultSimState';

  %% Run accelerator on TLC
  block.SetAccelRunOnTLC(true);

  %% Register methods
  block.RegBlockMethod('Outputs',                 @SetOutPortDims);
  block.RegBlockMethod('Outputs',                 @Output);

%endfunction

function SetOutPortDims(block)   %% SetOutputPortDims - Check and set output port dimensions
    block.InputPort(1).Dimensions  = 1;
    block.InputPort(2).Dimensions  = 1;
    block.InputPort(3).Dimensions  = 1;

    block.OutputPort(1).Dimensions = 15;

%endfunction

  %% Output
function Output(block)
%out_temp = 'A0';
temp_str_1 = num2str(block.InputPort(1).Data);
temp_str_2 = num2str(block.InputPort(2).Data);
temp_str_3 = num2str(block.InputPort(3).Data);
%string1 = strcat(out_temp,temp_str);
  block.OutputPort(1).CurrentDimensions = 15;
  block.OutputPort(1).Data(1) = block.DialogPrm(1).Data;
  block.OutputPort(1).Data(2) = 48;
  block.OutputPort(1).Data(3) = temp_str_1(1);
  block.OutputPort(1).Data(4) = temp_str_1(2);
  block.OutputPort(1).Data(5) = temp_str_1(3);

  block.OutputPort(1).Data(6) = block.DialogPrm(2).Data;
  block.OutputPort(1).Data(7) = 48;
  block.OutputPort(1).Data(8) = temp_str_2(1);
  block.OutputPort(1).Data(9) = temp_str_2(2);
  block.OutputPort(1).Data(10) = temp_str_2(3);

  block.OutputPort(1).Data(11) = block.DialogPrm(3).Data;
  block.OutputPort(1).Data(12) = 48;
  block.OutputPort(1).Data(13) = temp_str_3(1);
  block.OutputPort(1).Data(14) = temp_str_3(2);
  block.OutputPort(1).Data(15) = temp_str_3(3);

%endfunction

author: Sable

repo: mcbench-benchmarks


text:
function res = hbg(im,thresh)

% inputs
% im is the fourier transform of the image
% thresh is the cutoff circle radius

%outputs
% res is the boosted image

[r,c]=size(im);
d0=thresh;
d=zeros(r,c);
h=zeros(r,c);

for i=1:r
    for j=1:c
     d(i,j)=  sqrt( (i-(r/2))^2 + (j-(c/2))^2);
    end
end

A=1.75; % boost factor or coefficient

for i=1:r
    for j=1:c

         h(i,j)= 1-exp ( -( (d(i,j)^2)/(2*(d0^2)) ) );
         h(i,j)=(A-1)+h(i,j);
    end

end

for i=1:r
    for j=1:c
    res(i,j)=(h(i,j))*im(i,j);
    end
end


text:
% Generate sampling points on HEALPix
%
% OUT = HealpixGenerateSampling(n, mode)
% Parameters
% n : resolution of the grid (N_side)
% mode = 'rindex' : two-dimensional ring index
%      = 'scoord' : spherical coordinates

function OUT = HealpixGenerateSampling(n, mode)

n_total = 12 * n^2;

OUT = zeros(n_total, 2);
for pn = 1:n_total
    OUT(pn, :) = HealpixNestedIndex(n, pn - 1);
end

if mode == 'rindex'
    OUT = uint16(OUT);
    return
end

for pn = 1:n_total
    i = OUT(pn, 1);
    j = OUT(pn, 2);
    [theta, phi] = HealpixGetSphCoord(n, i, j);
    OUT(pn, :) = [theta phi];
end

author: mtex-toolbox

repo: mtex

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


text:
function [x,pnum,pgauss] = spm_mci_postslices (post,M,U,Y,Nbins)
% Univariate slices through posterior density
% FORMAT [x,pnum,pgauss] = spm_mci_postslices (post,M,U,Y,Nbins)
%
% post      posterior data structure
% M,U,Y     as usual
% Nbins     Number of bins per dimension
%
% x         [Np x Nbins] matrix where x(p,:) is domain for pth variable
% pnum      [Np x Nbins] where pnum(p,j) = p(x(p)=xj|x(\p),Y) ie. the posterior
%           density of variable p conditioned on the posterior mean of the other
%           variables. This is estimated numerically from evaluation of log joint
% pgauss    As pnum but under assumption that posterior is multivariate Gaussian
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_mci_postslices.m 6697 2016-01-27 14:57:28Z spm $

try, Nbins=Nbins; catch, Nbins=50; end

k=5; % Defines width of domain in number of SDs

% For computing log prior term
M = spm_mci_priors (M);
M = spm_mci_minit (M);

% Need high tolerances to avoid drop-outs in log-joint plots
M.reltol=1e-6;
M.abstol=1e-6;

Np=length(post.Ep);
Ep=post.Ep;
s=sqrt(diag(post.Cp));
pE=spm_vec(M.pE);
for p=1:Np,
    xmin=Ep(p)-k*s(p);
    xmax=Ep(p)+k*s(p);
    x(p,:)=linspace(xmin,xmax,Nbins);
    P = Ep;
    for j=1:Nbins,
        % Get parameters in eigenspace of prior
        P(p) = x(p,j);
        par = M.V'*(P-pE);
        eq(j) = spm_mci_joint (par,M,U,Y);
        pg(j) = spm_mvNpdf(P,Ep,post.Cp);
    end
    eq=exp(eq);
    pnum(p,:)=eq/sum(eq);
    pgauss(p,:)=pg/sum(pg);
end

figure
for i=1:Np,
    subplot(Np,1,i);
    plot(x(i,:),pgauss(i,:),'r');
    hold on
    plot(x(i,:),pnum(i,:),'k');
    legend('Gaussian','Numeric');
    xlabel(sprintf('P(%d)',i));
end

author: spm


text:
function cx = laguerre_poly ( n, x )

%*****************************************************************************80
%
%% LAGUERRE_POLY evaluates the Laguerre polynomials at X.
%
%  Differential equation:
%
%    X * Y'' + (1-X) * Y' + N * Y = 0
%
%  First terms:
%
%      1
%     -X    +  1
%   (  X^2 -  4 X     +  2 ) / 2
%   ( -X^3 +  9 X^2 -  18 X    +    6 ) / 6
%   (  X^4 - 16 X^3 +  72 X^2 -   96 X +      24 ) / 24
%   ( -X^5 + 25 X^4 - 200 X^3 +  600 X^2 -  600 x    +  120 ) / 120
%   (  X^6 - 36 X^5 + 450 X^4 - 2400 X^3 + 5400 X^2 - 4320 X + 720 ) / 720
%   ( -X^7 + 49 X^6 - 882 X^5 + 7350 X^4 - 29400 X^3
%      + 52920 X^2 - 35280 X + 5040 ) / 5040
%
%  Recursion:
%
%    L(0)(X) = 1,
%    L(1)(X) = 1-X,
%    N * L(N)(X) = (2*N-1-X) * L(N-1)(X) - (N-1) * L(N-2)(X)
%
%  Orthogonality:
%
%    Integral ( 0 <= X < Infinity ) exp ( - X ) * L(N)(X) * L(M)(X) dX
%    = 0 if N /= M
%    = 1 if N == M
%
%  Special values:
%
%    L(N)(0) = 1.
%
%  Relations:
%
%    L(N)(X) = (-1)^N / N! * exp ( x ) * (d/dx)^n ( exp ( - x ) * X^n )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%  Parameters:
%
%    Input, integer N, the highest order polynomial to compute.
%    Note that polynomials 0 through N will be computed.
%
%    Input, real X, the point at which the polynomials are to be evaluated.
%
%    Output, real CX(1:N+1), the Laguerre polynomials of degree 0 through
%    N evaluated at the point X.
%
  if ( n < 0 )
    cx = [];
    return
  end

  cx(1) = 1.0;

  if ( n == 0 )
    return
  end

  cx(2) = 1.0 - x;

  for i = 2 : n

    cx(i+1) = ( ( ( 2 * i - 1 ) - x ) * cx(i  )   ...
                - (     i - 1 )       * cx(i-1) ) ...
                / (     i );

  end

  return
end


text:
classdef deLaValleePoussinKernel
% This class is obsolet since MTEX 5.9. Use the class @SO3DeLaValleePoussinKernel
% instead. Anyway the class is preserved, so that the kernels of saved
% unimodal @ODFs can be loaded.

methods (Hidden = true)

  function psi = deLaValleePoussinKernel(varargin)

    warning(['The syntax "deLaValleePoussinKernel" is obsolete. ' ...
      'Please use "SO3DeLaValleePoussinKernel" instead.'])

    psi = SO3DeLaValleePoussinKernel(varargin{:});

  end
end

methods (Static = true, Hidden=true)
  function psi = loadobj(s)
    psi = SO3DeLaValleePoussinKernel(s.kappa);
  end
end

end


text:
clear all, close all, clc

A = [-.75 1; -.3 -.75];
B = [2; 1];
C = [1 2];
D = 0;

sys = ss(A,B,C,D);

Wc = gram(sys,'c'); % Controllability Gramian
Wo = gram(sys,'o'); % Observability Gramian

[sysb,g,Ti,T] = balreal(sys); % Balance the system

BWc = gram(sysb,'c') % Balanced Gramians
BWo = gram(sysb,'o')

%% Plot Gramians
theta = 0:.01:2*pi;
xc = cos(theta);
yc = sin(theta);
CIRC = [xc; yc];

ELLIPb = Ti*sqrt(BWc)*T*CIRC;
ELLIPc = sqrt(Wc)*CIRC;
ELLIPo = sqrt(Wo)*CIRC;

plot(xc,yc,'k--','LineWidth',2)
hold on

% Draw controllability Gramian (unbalanced)
plot(ELLIPc(1,:),ELLIPc(2,:),'r','LineWidth',2)
patch(ELLIPc(1,:),ELLIPc(2,:),'r','FaceAlpha',.75)

% Draw observability Gramian (unbalanced)
plot(ELLIPo(1,:),ELLIPo(2,:),'b','LineWidth',2)
patch(ELLIPo(1,:),ELLIPo(2,:),'b','FaceAlpha',.75)

% Draw balanced Gramians
patch(ELLIPb(1,:),ELLIPb(2,:),'k','FaceColor',[.5 0 .5],'FaceAlpha',.25)
plot(ELLIPb(1,:),ELLIPb(2,:),'Color',[.35 0 .35],'LineWidth',2)

%% Formatting
axis equal, grid on
set(gcf,'Position',[1 1 550 400])
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../figures/FIG_BT_GRAM');

%% Manually compute scaled balancing transformation
[Tu,D] = eig(Wc*Wo); % Tu are unscaled e-vecs

Atu = inv(Tu)*A*Tu;
Btu = inv(Tu)*B;
Ctu = C*Tu;
Dtu = 0;
syst = ss(Atu,Btu,Ctu,Dtu);

Sigmac = gram(syst,'c') % Diagonal Gramians
Sigmao = gram(syst,'o') % (but not equal)
Sigmas = diag(Sigmac)./diag(Sigmao);

% Scaled balancing transformation
T = Tu*diag(Sigmas.^(1/4));

% Permute columns of T to order Sigma
Sigma = diag(Sigmac).^(1/2).*diag(Sigmao).^(1/2);
[sigsort,permind] = sort(Sigma,'descend');
T = T(:,permind); % Hierarchical

% Compute balanced system
At = inv(T)*A*T;
Bt = inv(T)*B;
Ct = C*T;
Dt = 0;
sysBal = ss(At,Bt,Ct,Dt);

BWc = gram(sysBal,'c') % Balanced Gramians
BWo = gram(sysBal,'o')

%%
% ELLIPb = sqrt(BWc)*CIRC;

% A = [-1 2; 0 -1];
% B = [1; 1];
% C = [1 1];
% D = 0;

% A = [-1 2; -.5 -1];
% B = [1; 1];
% C = [1 1];
% D = 0;

author: dynamicslab


text:
function grid_index = multigrid_scale_closed ( dim_num, order_nd, level_max, ...
  level_1d, grid_index )

%*****************************************************************************80
%
%% MULTIGRID_SCALE_CLOSED renumbers a grid as a subgrid on a higher level.
%
%  Discussion:
%
%    This routine takes a grid associated with a given value of
%    LEVEL, and multiplies all the indices by a power of 2, so that
%    the indices reflect the position of the same points, but in
%    a grid of level LEVEL_MAX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 November 2007
%
%  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 ORDER_ND, the number of points in the grid.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Input, integer LEVEL_1D(DIM_NUM), the level in each dimension.
%
%    Input, integer GRID_INDEX(DIM_NUM,POINT_NUM), the index
%    values for each grid point, based in the level for which the grid
%    was generated.
%
%    Output, integer GRID_INDEX(DIM_NUM,POINT_NUM), the index
%    values for each grid point, appropriate for the grid as a subgrid
%    of a grid of level LEVEL_MAX.
%
  for dim = 1 : dim_num

    if ( level_1d(dim) == 0 )

      if ( 0 == level_max )
        order_max = 1;
      else
        order_max = 2^level_max + 1;
      end

      grid_index(dim,1:order_nd) = floor ( ( order_max - 1 ) / 2 );

    else

      factor = 2^( level_max - level_1d(dim) );

      grid_index(dim,1:order_nd) = grid_index(dim,1:order_nd) * factor;

    end

  end

  return
end


text:
function [WolfPops, Tracks] = BoundAdjust(WolfPops, ProbPoints, UAV)
%BOUNDADJUST 约束处理
Tracks = Pops2Tracks(WolfPops, UAV);

% ①  不满足约束："删除"航迹点，即取前后平均值
dim = UAV.PointDim;                                    % 维度
for agent = 1 : size(WolfPops.Pos, 1)
    Track = Tracks{agent};                               % 轨迹 (struct结构)
    ProbPoint = ProbPoints{agent};               % 问题点
    Position = [];                                               % 新狼群编码
    for i =1:UAV.num
          PointNum = UAV.PointNum(i);
          % 删除一条航迹上的问题点
          for k = 1 : PointNum
                flag = ProbPoint{i}(k);
                if flag == 1
                    if k == 1
                        P1 = UAV.S(i, :)' ;
                    else
                        P1 = Track.P{i}(:, k-1);
                    end
                    if k == PointNum
                        P2 = UAV.G(i, :)' ;
                    else
                        P2 = Track.P{i}(:, k+1);
                    end
                    Track.P{i}(:, k) = (P1+P2) / 2;  % 删除航迹点
                end
          end

          % 转换为狼群编码形式
          p = Track.P{i} ;
          p = reshape(p, 1, dim*PointNum);
          Position = [Position, p];
    end
    V = Track.V';
    Position = [Position, V];

    % 新的位置信息
    WolfPops.Pos(agent, :) = Position;
end

% ②  边界处理：越界取边界值
WolfPops.Pos = BoundClamp(WolfPops.Pos, WolfPops.lb, WolfPops.ub);

% 生成新航迹
Tracks = Pops2Tracks(WolfPops, UAV);

end

% 边界裁剪
function x = BoundClamp(x, lb, ub)
    Flag4ub = x > ub;
    Flag4lb = x < lb;
    x = x .* ( ~(Flag4ub + Flag4lb) ) + ub .* Flag4ub + lb .* Flag4lb;
end

repo: Grey-Wolf-Optimizer-for-Path-Planning


text:
function stats = optMyopic(returns, gamma, rf, accuracy)
%OPTMYOPIC: Calculate the optimal portfolio weight through the Myopic
%strategy or alternatively the Buy-and-Hold strategy
%
%   STATS = OPTMYOPIC(returns, ...) retrieves the optimal portfolio weight
%   using the myopic strategy. A grid search is used, based on the input
%   variable 'returns'. At any grid point, the conditional expectation is
%   calculated. In the end, the portfolio weight with the highest
%   conditional expectation is selected.
%   Alternatively, one can use this function to calculate the optimal
%   portfolio with the buy-and-hold strategy. For an investor with
%   investment horizon of 'K' periods, define the cumulative returns over
%   the past 'K' periods and use K*rf instead of rf.
%
%   Note that it is optimal to invest myopically when there are not state
%   variables and when the future states and returns are independent
%   (conditional on the state variables), in continuous time. Moreover,
%   long-term investors should invest myopically when returns are identical
%   and independent distributed.
%
%   'returns'	:   is a 'n-by-1' vector. One could use realized returns or
%                   simulated path returns (see also example).
%
%   The following assumptions are relevant:
%   -   A tradeoff is made between the input 'returns' and the risk free rate 'rf'
% 	-   The input consists only of the (simulated) returns.
%   -   Weights are restricted between 0 and 1.
%
%   STATS = OPTMYOPIC(returns, 'gamma', ...) allows you to specify an own
%   risk aversion parameter 'gamma'.
%       'gamma'   :     risk-aversion constant, initially set to 5.
%
%   STATS = OPTMYOPIC(returns, gamma, 'rf', ...) allows you to specify the
%   risk free rate, relevant for your optimization problem.
%       'rf'      :     the risk free rate constant, initial set to 0.
%
%   STATS = OPTMYOPIC(returns, gamma, rf, 'accuracy') allows you to extend
%   the accuracy of the grid, initially set to 0.01.
%       'accuracy':     the single step length in the grid
%
%   STATS = OPTMYOPIC(...) returns the structure 'stats' consisting of:
%       'stats.optWeight'   :  the optimal portfolio weight (in risky asset)
%       'stats.optCU'       :  optimal conditional expectation
%
%   NB:     -   For three assets we need a two dimensional grid, which
%               leads appearantly to the curse of dimensionality. There are
%               methods that solve this problem, however such solutions are
%               beyond the scope of this code.
%
%  Semin Ibisevic (2012)
%  $Date: 02/05/2012$
%
% -------------------------------------------------------------------------
% References
%
% B.F. Diris. Portfolio Management. Econometric Institute, 2012. Lecture
% FEM21010.
% -------------------------------------------------------------------------

[n, m] = size(returns);
if m >= n,  error('optMyopic:InvalidInput','The number of rows should be sufficiently larger than the number of columns'); end
if m > 1,   error('optMyopic:InvalidInput','The input should consist only of one column. Note: to calculate the portfolio weights using the Buy-and-Hold strategy, use cumulative returns over the past k periods instead.'); end

if nargin < 1, error('optMyopic:InsufficientParameters','Input variable is necessary'); end
if nargin < 2, gamma = 5; end   % default risk aversion
if nargin < 3, rf = 0; end      % default risk free rate
if nargin < 4, accuracy = 0.01; end

% One-dimensional grid
grid = 0.00:accuracy:1.00;

% init
U = zeros(n,size(grid,2));
CU = zeros(1,size(grid,2));
CUopt = -10^(9);
gridopt = 0;

% -------------------------------------------------------------------------

% pick portfolio weight from the grid
wIter = 1;
for w = grid

   % define the realized utility value for all i = 1,...,n
   U(:,wIter) = (1/(1-gamma)) * ( w*exp( returns ) + (1-w)*exp( rf ) ).^(1-gamma);

   % approximate the conditional expected utlity for grid point w
   CU(wIter) = (1/n)*sum( U(:,wIter) );

   % update
   if CU(wIter) > CUopt
       CUopt = CU(wIter);
       gridopt = w;
   end

   wIter = wIter + 1;

end

stats.optWeight = gridopt;
stats.optCU = CUopt;


text:
function line_fekete_rule_bos_levenberg_test ( m )

%*****************************************************************************80
%
%% LINE_FEKETE_RULE_BOS_LEVENBERG_TESTS: Bos Levenberg code for Fekete points.
%
%  Discussion:
%
%    I wanted to run the code in the Bos Levenberg paper for
%    comparison with mine.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 March 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Len Bos, Norm Levenberg,
%    On the calculation of approximate Fekete points: the univariate case,
%    Electronic Transactions on Numerical Analysis,
%    Volume 30, pages 377-397, 2008.
%
%  Parameters:
%
%    Input, integer M, the dimension of the polynomial space.
%    In the paper, M is 21.
%
  if ( nargin < 1 )
    m = 5;
  end
  n = 1000;
  a = -1.0;
  b = +1.0;
  x = linspace ( a, b, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LINE_FEKETE_RULE_BOS_LEVENBERG_TEST:\n' );
  fprintf ( 1, '  Seek Fekete points in [%g,%g]\n', a, b );
  fprintf ( 1, '  using %d equally spaced sample points\n', n );
  fprintf ( 1, '  for polynomial space of dimension M = %d\n', m );
  fprintf ( 1, '  with the Chebyshev basis\n' );
  fprintf ( 1, '  and weight 1/sqrt(1-x^2).\n' );

  A = chebvand ( m, x );
  b = rand ( m, 1 );
  y = A \ b;
  pp = ( y ~= 0.0 );
  xf = x(pp);
  nf = length ( xf );
  r8vec_print ( nf, xf, '  Estimated Fekete points XF:' );

  yf = ones ( nf, 1 );
  plot ( xf, yf, '*' );
  title ( 'Estimated Fekete point locations, Bos-Levenberg' );
  grid on

  filename = 'line_fekete_rule_bos_levenberg_test.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saved plot in file "%s"\n', filename );

  return
end


text:
% isPalindromeNumber: returns true if the given number are a palindrome
%                     otherwise false.
% assumes: the given number is positive
function ans = isPalindromeNumber(number)
  % precondition
  assert(number >= 0,'number must be positive')
  ans = false;
  strNumber = num2str(number); % convert the given number to character array
  if (strNumber == flip(strNumber))
    ans = true;
  endif
endfunction

author: TheAlgorithms

repo: MATLAB-Octave


text:
function patches = load_hpatches(impath)
%LOAD_PATCHES Load and cut HPatches patch file

patches = imread(impath);
npatches = size(patches, 1) ./ size(patches, 2);
assert(mod(npatches, 1) == 0, 'Invalid patch image height.');
patches = mat2cell(patches, size(patches, 2)*ones(1, npatches), size(patches, 2));
patches = permute(patches, [2, 3, 4, 1]);
patches = cell2mat(patches);
end

author: hpatches

repo: hpatches-benchmark


text:
function drawBezier(p,color,weight)
%drawBezier Draw a bezier curve in a figure
%
% drawBezier(p,color,weight)
%
%INPUT
% p         Coordinates
% color     Color of bezier curve
% weight    Weight of bezier curve
%
global CB_MAP_OUTPUT
global mapHandle
if nargin < 2
    color = [0 0 255];
    display('No color specified');
end
numpoints = 50;
%a = zeros(2,numpoints);
%b = zeros(2,numpoints);
%c = zeros(2,numpoints);
i = 0:numpoints;
a = (p(:,1)* i + p(:,2)*(numpoints-i))/numpoints;
b = (p(:,2)* i + p(:,3)*(numpoints-i))/numpoints;
c = (a .* (ones(2,1)* i) + b .*(ones(2,1)  *(numpoints-i)))  /numpoints;
if strcmp(CB_MAP_OUTPUT, 'matlab')
    if find(color>1)
        color = color/255;
    end
    line(c(1,:), -c(2,:),'Color',color,'LineWidth',weight);
elseif strcmp(CB_MAP_OUTPUT, 'java')
    %line(c(1,:), c(2,:),'Color',color,'LineWidth',weight);
    % fill in code
     setDataBezier(mapHandle,c(1,:),c(2,:));
elseif strcmp(CB_MAP_OUTPUT, 'svg')
    %determine type of color input
    if ischar(color)
        colorStroke = color;
    else if isvector(color)
            colorStroke = strcat('rgb(',num2str(color(1)),',',num2str(color(2)),',',num2str(color(3)),')');
        end
    end
    fprintf(mapHandle,'<g id="" stroke="%s" stroke-width="%d" stroke-linecap="round">\n',colorStroke,ceil(weight));
%     fprintf(mapHandle,'<g id="" stroke="deepskyblue" stroke-width="6" stroke-linecap="round">\n');
    %fprintf(mapHandle,'<path style="fill: none;" d="M%8.2f %8.2f C%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f"/>\n',p2(1),-p2(2),p2(1),-p2(2),ptemp(1),-ptemp(2),p1(1),-p1(2));
    fprintf(mapHandle,'<path style="fill: none;" d="M%8.2f %8.2f C%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f"/>\n',p(1,3),p(2,3) ,p(1,3),p(2,3),p(1,2),p(2,2),p(1,1),p(2,1));
    fprintf(mapHandle,'</g>\n');
else
    display('error CB_MAP_OUTPUT in bezier');
end

author: opencobra

repo: cobratoolbox


text:
function [beta_gibbs,sigma_gibbs]=ndgibbstotal(It,Bu,X,Y,y,Bhat,n,T,q)

% modified ndgibbs sampler: insensitive to hyperparameters, total diffuse, "flat" prior in spirit of Uhlig (2005)

% function [beta_gibbs sigma_gibbs]=bear.ndgibbs(It,Bu,beta0,omega0,X,Y,y,Bhat,n,T,q)
% performs the Gibbs algorithm 1.5.2 for the normal-diffuse prior, and returns draws from posterior distribution
% inputs:  - integer 'It': total number of iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - integer 'Bu': number of burn-in iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - vector 'beta0': vector of prior values for beta (defined in 1.3.4)
%          - matrix 'omega0': prior covariance matrix for the VAR coefficients (defined in 1.3.8)
%          - matrix 'X': matrix of regressors for the VAR model (defined in 1.1.8)
%          - matrix 'Y': matrix of regressands for the VAR model (defined in 1.1.8)
%          - vector 'y': vectorised regressands for the VAR model (defined in 1.1.12)
%          - matrix 'Bhat': OLS VAR coefficients, in non vectorised form (defined in 1.1.9)
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'T': number of sample time periods (defined p 7 of technical guide)
%          - integer 'q': total number of coefficients to estimate for the BVAR model (defined p 7 of technical guide)
% outputs: - matrix 'beta_gibbs': record of the gibbs sampler draws for the beta vector
%          - matrix'sigma_gibbs': record of the gibbs sampler draws for the sigma matrix (vectorised)

% preliminary tasks

% % invert omega0, as it will be used repeatedly during step 4
% invomega0=diag(1./diag(omega0));

% set initial values for B (step 2); use OLS estimates
B=Bhat;

% create the progress bar
hbar = bear.parfor_progressbar(It,'Progress of the Gibbs sampler.');

% start iterations
for ii=1:It

% Step 3: at iteration ii, first draw sigma from IW, conditional on beta from previous iteration
% obtain first Shat, defined in (1.6.10)
Shat=(Y-X*B)'*(Y-X*B);
% Correct potential asymmetries due to rounding errors from Matlab
C=chol(bear.nspd(Shat));
Shat=C'*C;

% next draw from IW(Shat,T)
sigma=bear.iwdraw(Shat,T);

% step 4: with sigma drawn, continue iteration ii by drawing beta from a multivariate Normal, conditional on sigma obtained in current iteration
% first invert sigma
C=chol(bear.nspd(sigma));
invC=C\speye(n);
invsigma=invC*invC';

% then obtain the omegabar matrix
invomegabar=kron(invsigma,X'*X);
C=chol(bear.nspd(invomegabar));
invC=C\speye(q);
omegabar=invC*invC';

% following, obtain betabar
betabar=omegabar*(kron(invsigma,X')*y);

% draw from N(betabar,omegabar);
beta=betabar+chol(bear.nspd(omegabar),'lower')*mvnrnd(zeros(q,1),eye(q))';

% % update matrix B with each draw
% B=reshape(beta,size(B));

% record the values if the number of burn-in iterations is exceeded
if ii>Bu
% values of vector beta
beta_gibbs(:,ii-Bu)=beta;
% values of sigma (in vectorized form)
sigma_gibbs(:,ii-Bu)=sigma(:);
% if current iteration is still a burn iteration, do not record the result
else
end

% update progress by one iteration
hbar.iterate(1);

% go for next iteration
end
% close progress bar
close(hbar);

author: european-central-bank

repo: BEAR-toolbox


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% Problem 7- Frequency Response and Impulse Response of an ideal bandpass
% filter

td=.05;
B1=10 ;
B2=20;
w1=-30:-B2;
H1=zeros(size(w1));
w2=-B2:-B1;
H2=exp(-j*w2*td);
w3=-B1:B1;
H3=zeros(size(w3));
w4=B1:B2;
H4= exp(-j*w4*td);
w5=B2:30;
H5=zeros(size(w5));
w=[w1 w2 w3 w4 w5];
H=[H1 H2 H3 H4 H5];
plot(w,abs(H))
ylim([-.2 1.2]);
legend(' |H(\Omega) | ')

figure
plot(w,angle(H));
ylim([-1.1 1.1]);
legend('\angle H(\Omega)')

figure
syms t w
H1=exp(-j*w*td)*(heaviside(w+B2)-heaviside(w+B1));
H2=exp(-j*w*td)*(heaviside(w-B1)-heaviside(w-B2));
h1=ifourier(H1,t);
h2=ifourier(H2,t);
h=h1+h2;
ezplot(h, [-3 3]);
ylim([-4 4]);
legend('Impulse response h(t) ')


text:
function  [value, isterminal, direction] = events_flight(~,Z,P)

LL = P.FullLength;
th = Z(1,:);
y0 = Z(3,:);
y2 = y0 + LL*cos(th);

%%%% CONSTRAINTS %%%%
% 1) y0 > 0
% 2) y2 > 0

% As a rule, each constraint will be satisfied if it's event function value
% is positive. This makes things easier at the FSM level.
n = length(th);
value = zeros(2,n);
isterminal = true(size(value));
direction = -ones(size(value));

value(1,:) = y0;
value(2,:) = y2;

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function model = linearParamInit(model)

% LINEARPARAMINIT Initialise the parameters of an LINEAR model.
% FORMAT
% DESC sets the initial weight vectors and biases to small random
% values.
% ARG model : the input model to initialise.
% RETURN model : the initialised model.
%
% SEEALSO : modelParamInit, linearCreate
%
% COPYRIGHT : Neil D. Lawrence, 2006

% MLTOOLS

model.W = randn(model.inputDim, model.outputDim)/sqrt(model.inputDim + 1);
model.b = randn(1, model.outputDim)/sqrt(model.inputDim + 1);
model.beta = 1;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function value = e_constant ( )

%*****************************************************************************80
%
%% E_CONSTANT returns the value of E.
%
%  Discussion:
%
%   "E" was named in honor of Euler.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real E_CONSTANT, the base of the natural
%    logarithm system.
%
  value = 2.718281828459045;

  return
end


text:
classdef SettingsNumericalHomogenizer < AbstractSettings

    properties (Access = public)
        interpParams
        levelSetCreatorParams
        volumeShFuncParams
        elementDensityCreatorType
        outFileName
        testName
        print
        iter
        pdim
    end

    properties (Access = protected)
        defaultParamsName = 'paramsNumericalHomogenizer'
    end

    methods (Access = public)

        function obj = SettingsNumericalHomogenizer(varargin)
            if nargin == 1
                    obj.loadParams(varargin{1});
            end
            obj.setParams();
        end

    end

    methods (Access = private)

        function setParams(obj)
            obj.volumeShFuncParams.filename = obj.testName;
            obj.volumeShFuncParams.scale = 'MICRO';
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
classdef SeqObsModel_Gaussian < SeqObsModel
    % SeqObsModel : defines Gaussian emission model for observed data in HMM
    %   Allows for multivariate observations of any dimension.
    %   Assumes conjugate prior on mean and covariance (Normal-Inverse Wishart).
    %   Extends the base object "SeqObsModel", filling in Gaussian details
    methods
        % ====================================================== CONSTRUCTOR
        function obj = SeqObsModel_Gaussian( K, D)
            obj.K = K;
            obj.D = D;
            obj.theta = repmat( struct('mu',zeros(D,1),'invSigma', zeros(D,D) ), 1, K );
        end

        % setPrior() sets prior on emission parameters to specified values.
        % USAGE:
        %   To set to default values given data object and model params
        %     ObsM = ObsM.setPrior( data, modelParamStruct )
        %     where modelParamStruct has attributes:
        %        degFree : # deg of freedom
        %        precMu  : precision parameter for mean
        %        Scoef   : scalar multiplier for Scale Matrix
        %        doEmpCovScalePrior: boolean.
        %               true  = set ScaleMat to Scoef*emp covariance
        %               false = set ScaleMat to Scoef*eye(D)
        %  To initialize directly to given parameters
        %      ObsM = ObsM.setPrior( mu, precMu, degFree, ScaleMat )
        function obj = setPrior( obj, varargin )
            obj.priorDef = struct();
            obj.priorDef.doEmpCovScalePrior =0;
            obj.priorDef.Scoef = 0;
            if strfind( class(varargin{1} ), 'SeqData' )
                data = varargin{1};
                obsM = varargin{2};
                obj.prior.mu = zeros( data.D, 1 );
                obj.prior.precMu = obsM.precMu;
                obj.prior.degFree = max( obsM.degFree, data.D+2 );
                obj.priorDef.doEmpCovScalePrior = obsM.doEmpCovScalePrior;
                obj.priorDef.Scoef = obsM.Scoef;
                if obsM.doEmpCovScalePrior
                    obj.prior.ScaleMat = obsM.Scoef * cov( data.Xdata', 1);
                else
                    obj.prior.ScaleMat = obsM.Scoef * eye( data.D );
                end
            elseif length( varargin ) == 1
                PStruct = varargin{1};
                obj.prior.mu = PStruct.Mu(:);
                obj.prior.precMu = PStruct.precMu;
                obj.prior.degFree = PStruct.degFree;
                obj.prior.ScaleMat = PStruct.ScaleMat;
            else
                obj.prior.mu = varargin{1}(:);
                obj.prior.precMu = varargin{2};
                obj.prior.degFree = varargin{3};
                obj.prior.ScaleMat = varargin{4};
            end
        end

        % ==================================================== GET Xstats
        function Xstats = getXSuffStats( obj, Xkk )
            Xstats.nObs = size(Xkk,2);
            if Xstats.nObs > 0
                Xstats.Xsum = sum( Xkk, 2);
                Xstats.XXsum = Xkk*Xkk';
            else
                Xstats.Xsum = [];
                Xstats.XXsum = [];
            end
        end

        function obj = incXStats( obj, ii, data, stateSeq, featIDs )
            for kk = featIDs
                Xkk = data.seq(ii);
                Xkk = Xkk(:, stateSeq(ii).z == kk );
                nNew = size(Xkk,2);
                if length(obj.Xstats)<kk || obj.Xstats(kk).nObs == 0
                    obj.Xstats(kk).nObs = nNew;
                    obj.Xstats(kk).Xsum = sum(Xkk,2);
                    obj.Xstats(kk).XXsum = Xkk*Xkk';
                else
                    obj.Xstats(kk).nObs = obj.Xstats(kk).nObs + nNew;
                    obj.Xstats(kk).Xsum = obj.Xstats(kk).Xsum + sum(Xkk,2);
                    obj.Xstats(kk).XXsum = obj.Xstats(kk).XXsum + Xkk*Xkk';
                end
            end
        end

        function obj = decXStats( obj, ii, data, stateSeq, featIDs )
            for kk = featIDs
                Xkk = data.seq(ii);
                Xkk = Xkk(:, stateSeq(ii).z == kk );
                nNew = size(Xkk,2);
                if obj.Xstats(kk).nObs-nNew <= 0
                    obj.Xstats(kk).nObs=0;
                    obj.Xstats(kk).Xsum=[];
                    obj.Xstats(kk).XXsum=[];
                else
                    obj.Xstats(kk).nObs = obj.Xstats(kk).nObs - nNew;
                    obj.Xstats(kk).Xsum = obj.Xstats(kk).Xsum - sum(Xkk,2);
                    obj.Xstats(kk).XXsum = obj.Xstats(kk).XXsum - Xkk*Xkk';
                end
            end
        end

        % ==================================================== GET params
        function retStr = getParamDescr( obj, PP )
            if ~exist('PP','var')
                PP = obj.prior;
            end
            retStr = 'Norm-InvWish. Mu0=0, dFree=%d,';
            if obj.priorDef.doEmpCovScalePrior
                retStr = [retStr ' S0=%.2f*EmpCov.'];
            else
                retStr = [retStr ' S0=%.2f*eye.'];
            end
            retStr = sprintf( retStr, obj.prior.degFree, obj.priorDef.Scoef );
        end

        function PP = getPosteriorParams( obj, Xstats )
            if Xstats.nObs > 0
                N = Xstats.nObs;
                degFreeN = N + obj.prior.degFree;
                precMuN  = N + obj.prior.precMu;
                Xdiff = Xstats.Xsum/N - obj.prior.mu;
                XdiffMat = ( obj.prior.precMu * N  )/( obj.prior.precMu + N ) .* ( Xdiff*Xdiff' );
                CovMat = (Xstats.XXsum - Xstats.Xsum*Xstats.Xsum'/N);
                ScaleMatN = obj.prior.ScaleMat + CovMat + XdiffMat;
                MuN = (Xstats.Xsum + obj.prior.mu*obj.prior.precMu)/ ( precMuN );

                PP.mu = MuN;
                PP.precMu = precMuN;
                PP.degFree = degFreeN;
                PP.ScaleMat = ScaleMatN;
            else
                PP = obj.prior;
            end
        end

        % ==================================================== GET theta mean

        % Retrieve the expected value of theta under the given params
        % In Gaussian case,
        %   mu = PP.mu
        %   Sigma = PP.ScaleMat / (degFree-D-1 )
        function [theta] = getTheta_Mean( obj, PP )
            if ~exist('PP','var')
                PP = obj.prior;
            end
            theta.mu  = PP.mu;
            theta.invSigma = (PP.degFree-obj.D-1)*( PP.ScaleMat \ eye(obj.D) );
        end

        % ==================================================== GET theta samps
        function [theta, PP] = sampleTheta_FromParams( obj, PP )
            if ~exist('PP','var')
                PP = obj.prior;
            end
            [~, sqrtInvSigma] = randiwishart( PP.ScaleMat, PP.degFree );
            theta.mu = ( sqrt( PP.precMu ) .* sqrtInvSigma ) \randn(obj.D,1)  + PP.mu;
            theta.invSigma = sqrtInvSigma'*sqrtInvSigma;
        end

        function [theta, PN] = sampleTheta( obj, Xkk )
            if ~exist( 'Xkk', 'var' )
                Xkk = [];
            end
            Xstats = obj.getXSuffStats( Xkk );
            PN = obj.getPosteriorParams( Xstats );
            [theta] = obj.sampleTheta_FromParams( PN );
        end

        % ==================================================== GET log probs

        % Calc probability of a realization of a emission parameter "theta"
        %   given the provided prior parameter struct "PP"
        function logPr = calcLogPrTheta( obj, theta, PP )
            if ~exist('PP','var')
                PP = obj.prior;
            end
            logPr = calcLogPrNormalInvWishart( theta.mu, theta.invSigma, PP );
        end

        % Calculate soft evidence for particular sequence "ii"
        %   given all available emission parameters (indicated by "kIDs")
        %OUTPUT
        %  logSoftEv : K x T matrix, where entry k,t gives log probability
        %                of the t-th observation under theta(kk)
        function logSoftEv = calcLogSoftEv(obj, ii, data, kIDs )
            if ~exist( 'kIDs', 'var' )
                kIDs = 1:obj.K;
            end
            Xseq = data.seq(ii);
            T = size( Xseq,2);
            logSoftEv = -inf( obj.K, T );
            for kk = kIDs
                cholInvSigma = chol( obj.theta(kk).invSigma );
                logDetInvSigma = 2*sum( log( diag( cholInvSigma) ) );
                XdiffMu = bsxfun(@minus, Xseq, obj.theta(kk).mu );
                U = XdiffMu'*cholInvSigma';
                logSoftEv(kk,:) = 0.5*logDetInvSigma - 0.5*sum( U.^2,2);
            end
            logSoftEv = logSoftEv - 0.5*obj.D*log(2*pi);
        end

        % Calculate marginal probability of observed data
        %  given the state sequence assignments of every observation
        % USAGE:
        %  logPr = ThetaM.calcMargPrData()
        %    calculate using *stored* sufficient statistics in this object
        %  logPr = ThetaM.calcMargPrData( data, stateSeq )
        %    calculate using freshly computed suff stats
        function logPr = calcMargPrData(obj, data, stateSeq, ks )
            if ~exist( 'ks', 'var' )
                ks = 1:obj.K;
            end
            if exist( 'data', 'var' ) && ~isempty( data )
                obj = obj.updateAllXSuffStats(  horzcat(stateSeq(:).z), data, ks);
            end
            LOG_PI =  1.144729885849400;
            P0 = obj.prior;
            D = obj.D;
            logPr = zeros(1, obj.K);
            for kk = ks
                Nkk = obj.Xstats(kk).nObs;
                if Nkk > 0
                    PN = obj.getPosteriorParams( obj.Xstats(kk) );

                    logPr(kk) = ...
                        logMvGamma(0.5*PN.degFree,D) - logMvGamma(0.5*P0.degFree,D) ...
                        + 0.5*P0.degFree*log( det( P0.ScaleMat ) )   ...
                        - 0.5*PN.degFree*log( det( PN.ScaleMat ) ) ...
                        + 0.5*D*log( P0.precMu / PN.precMu );
                end
            end
            N = sum( [obj.Xstats(ks).nObs] );
            logPr = sum(logPr) - 0.5*N*D*LOG_PI;
        end

    end
end

author: michaelchughes

repo: NPBayesHMM


text:
function [r,ir,iu] = unique(r,varargin)
% disjoint list of rotations
%
% Syntax
%   u = unique(r)
%   u = unique(r,'tolerance',0.01)
%   [u,ir,iu] = unique(r)
%
% Input
%  r   - @rotation
%  tol - double (default 1e-3)
%
% Output
%  u - @rotation
%  ir - index such that u = r(ir)
%  iu - index such that r = u(iu)
%
% Flags
%  stable - prevent sorting
%
% See also
% unique
%

a = r.a(:); b = r.b(:); c = r.c(:); d = r.d(:); i = r.i(:);

if length(r) < 1000 && ~check_option(varargin,'tolerance')

  abcd = [a.^2,b.^2,c.^2,d.^2,a.*b,a.*c,a.*d,b.*c,b.*d,c.*d, i(:)];

  tol = get_option(varargin,'tolerance',1e-3);

  % in case it should not be sorted
  if check_option(varargin,'stable')
    [~,ir,iu] = unique(round(abcd./tol),'rows','stable');
  else
    [~,ir,iu] = uniquetol(abcd,tol,'ByRows',true,'DataScale',1);
  end

else % faster but less accurate

  tol = get_option(varargin,'tolerance',5*degree);

  % ensure upper hemisphere
  ind = d < 0;
  a(ind) = -a(ind); b(ind) = -b(ind);
  c(ind) = -c(ind); d(ind) = -d(ind);

  [~,ir,iu] = uniquetol([a,b,c,d,i],tol / 100 / degree, ...
    'ByRows',true,'DataScale',1);

end

% remove duplicated points
r.a = a(ir);
r.b = b(ir);
r.c = c(ir);
r.d = d(ir);
r.i = i(ir);


text:
function [ ] = DrawQuadrotor( x,y )
%DRAWQUADROTOR Summary of this function goes here
%   Detailed explanation goes here
rotor_R=0.15;
wheelbase=0.2;
offset=(wheelbase/2+rotor_R)/sqrt(2);
rectangle('Position',[x-offset-rotor_R,y-offset-rotor_R,2*rotor_R,2*rotor_R],'Curvature',[1,1],'edgecolor','k'),axis equal
rectangle('Position',[x-offset-rotor_R,y+offset-rotor_R,2*rotor_R,2*rotor_R],'Curvature',[1,1],'edgecolor','r'),axis equal
rectangle('Position',[x+offset-rotor_R,y-offset-rotor_R,2*rotor_R,2*rotor_R],'Curvature',[1,1],'edgecolor','k'),axis equal
rectangle('Position',[x+offset-rotor_R,y+offset-rotor_R,2*rotor_R,2*rotor_R],'Curvature',[1,1],'edgecolor','r'),axis equal
end

author: heartxuxuxu


text:
function ll = robOneDynamicsLogLikelihood(model)

% ROBONEDYNAMICSLOGLIKELIHOOD Give the log likelihood of the robot one dynamics part.

% FGPLVM

thetaDiff = model.theta(2:end)-model.theta(1:end-1);
while any(thetaDiff>pi)
  ind = find(thetaDiff>pi);
  thetaDiff(ind) = thetaDiff(ind) - 2*pi;
end
while any(thetaDiff<-pi)
  ind = find(thetaDiff<-pi);
  thetaDiff(ind) = thetaDiff(ind) + 2*pi;
end

logLikeTheta = log(2*pi*model.sigma2)...
    +(thetaDiff.*thetaDiff)/model.sigma2;
likeTheta = exp(-0.5*logLikeTheta)*model.mixTheta;
likeTheta = likeTheta + (1-model.mixTheta)/(2*pi);
ll = sum(log(likeTheta));

logLikeR1 = model.a*log(model.b) - gammaln(model.a) ...
    +(model.a-1)*log(model.r) - model.b*model.r + log(model.mixR);
logLikeR2 = log(model.b) - model.b*model.r + log(1-model.mixR);

logLikeR1(find(logLikeR1<-316))=-316;
logLikeR2(find(logLikeR2<-316))=-316;

ind = find(logLikeR1>logLikeR2);
logLikeR(ind) = log(1+exp(logLikeR2(ind))./exp(logLikeR1(ind)))+logLikeR1(ind);

ind2 = find(logLikeR1<=logLikeR2);
logLikeR(ind2) = log(1+exp(logLikeR1(ind2))./exp(logLikeR2(ind2)))+logLikeR2(ind2);

ll = ll + sum(logLikeR);


text:
function v = cheby_t_poly ( m, n, x )

%*****************************************************************************80
%
%% CHEBY_T_POLY evaluates Chebyshev polynomials T(n,x).
%
%  Discussion:
%
%    Chebyshev polynomials are useful as a basis for representing the
%    approximation of functions since they are well conditioned, in the sense
%    that in the interval [-1,1] they each have maximum absolute value 1.
%    Hence an error in the value of a coefficient of the approximation, of
%    size epsilon, is exactly reflected in an error of size epsilon between
%    the computed approximation and the theoretical approximation.
%
%    Typical usage is as follows, where we assume for the moment
%    that the interval of approximation is [-1,1].  The value
%    of N is chosen, the highest polynomial to be used in the
%    approximation.  Then the function to be approximated is
%    evaluated at the N+1 points XJ which are the zeroes of the N+1-th
%    Chebyshev polynomial.  Let these values be denoted by F(XJ).
%
%    The coefficients of the approximation are now defined by
%
%      C(I) = 2/(N+1) * sum ( 1 <= J <= N+1 ) F(XJ) T(I,XJ)
%
%    except that C(0) is given a value which is half that assigned
%    to it by the above formula,
%
%    and the representation is
%
%    F(X) approximated by sum ( 0 <= J <= N ) C(J) T(J,X)
%
%    Now note that, again because of the fact that the Chebyshev polynomials
%    have maximum absolute value 1, if the higher order terms of the
%    coefficients C are small, then we have the option of truncating
%    the approximation by dropping these terms, and we will have an
%    exact value for maximum perturbation to the approximation that
%    this will cause.
%
%    It should be noted that typically the error in approximation
%    is dominated by the first neglected basis function (some multiple of
%    T(N+1,X) in the example above).  If this term were the exact error,
%    then we would have found the minimax polynomial, the approximating
%    polynomial of smallest maximum deviation from the original function.
%    The minimax polynomial is hard to compute, and another important
%    feature of the Chebyshev approximation is that it tends to behave
%    like the minimax polynomial while being easy to compute.
%
%    To evaluate a sum like
%
%      sum ( 0 <= J <= N ) C(J) T(J,X),
%
%    Clenshaw's recurrence formula is recommended instead of computing the
%    polynomial values, forming the products and summing.
%
%    Assuming that the coefficients C(J) have been computed
%    for J = 0 to N, then the coefficients of the representation of the
%    indefinite integral of the function may be computed by
%
%      B(I) = ( C(I-1) - C(I+1))/2*(I-1) for I=1 to N+1,
%
%    with
%
%      C(N+1)=0
%      B(0) arbitrary.
%
%    Also, the coefficients of the representation of the derivative of the
%    function may be computed by:
%
%      D(I) = D(I+2)+2*I*C(I) for I=N-1, N-2, ..., 0,
%
%    with
%
%      D(N+1) = D(N)=0.
%
%    Some of the above may have to adjusted because of the irregularity of C(0).
%
%  Differential equation:
%
%    (1-X*X) Y'' - X Y' + N N Y = 0
%
%  Formula:
%
%    T(N,X) = COS(N*ARCCOS(X))
%
%  First terms:
%
%    T(0,X) =  1
%    T(1,X) =  1 X
%    T(2,X) =  2 X^2 -   1
%    T(3,X) =  4 X^3 -   3 X
%    T(4,X) =  8 X^4 -   8 X^2 +  1
%    T(5,X) = 16 X^5 -  20 X^3 +  5 X
%    T(6,X) = 32 X^6 -  48 X^4 + 18 X^2 - 1
%    T(7,X) = 64 X^7 - 112 X^5 + 56 X^3 - 7 X
%
%  Inequality:
%
%    abs ( T(N,X) ) <= 1 for -1 <= X <= 1
%
%  Orthogonality:
%
%    For integration over [-1,1] with weight
%
%      W(X) = 1 / sqrt(1-X*X),
%
%    if we write the inner product of T(I,X) and T(J,X) as
%
%      < T(I,X), T(J,X) > = integral ( -1 <= X <= 1 ) W(X) T(I,X) T(J,X) dX
%
%    then the result is:
%
%      0 if I /= J
%      PI/2 if I == J /= 0
%      PI if I == J == 0
%
%    A discrete orthogonality relation is also satisfied at each of
%    the N zeroes of T(N,X):  sum ( 1 <= K <= N ) T(I,X) * T(J,X)
%                              = 0 if I /= J
%                              = N/2 if I == J /= 0
%                              = N if I == J == 0
%
%  Recursion:
%
%    T(0,X) = 1,
%    T(1,X) = X,
%    T(N,X) = 2 * X * T(N-1,X) - T(N-2,X)
%
%    T'(N,X) = N * ( -X * T(N,X) + T(N-1,X) ) / ( 1 - X^2 )
%
%  Special values:
%
%    T(N,1) = 1
%    T(N,-1) = (-1)^N
%    T(2N,0) = (-1)^N
%    T(2N+1,0) = 0
%    T(N,X) = (-1)^N * T(N,-X)
%
%  Zeroes:
%
%    M-th zero of T(N,X) is cos((2*M-1)*PI/(2*N)), M = 1 to N
%
%  Extrema:
%
%    M-th extremum of T(N,X) is cos(PI*M/N), M = 0 to N
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of evaluation points.
%
%    Input, integer N, the highest polynomial to compute.
%
%    Input, real X(M,1), the evaluation points.
%
%    Output, real V(M,N+1), the values of the Chebyshev polynomials
%    0 through N at X(1:M).
%
  if ( n < 0 )
    c = [];
    return
  end

  v = zeros ( m, n + 1 );

  v(1:m,1) = 1.0;

  if ( n < 1 )
    return
  end

  x = x(:);

  v(1:m,2) = x(1:m,1);

  for i = 2 : n
    v(1:m,i+1) = 2.0 * x(1:m,1) .* v(1:m,i) - v(1:m,i-1);
  end

  return
end


text:
function CMap = cmap_jetinv(nbVals)
% CMAP_JETINV: Inversed JET colormap

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

CMap = bst_flip(jet(nbVals), 1);

author: brainstorm-tools


text:
function [hess, storedb] = getHessian(problem, x, d, storedb)
% Computes the Hessian of the cost function at x along d.
%
% function [hess, storedb] = getHessian(problem, x, d, storedb)
%
% Returns the Hessian at x along d of the cost function described in the
% problem structure. The cache database storedb is passed along, possibly
% modified and returned in the process.
%
% If an exact Hessian is not provided, an approximate Hessian is returned
% if possible, without warning. If not possible, an exception will be
% thrown. To check whether an exact Hessian is available or not (typically
% to issue a warning if not), use canGetHessian.
%
% See also: getPrecon getApproxHessian canGetHessian

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:

    if isfield(problem, 'hess')
    %% Compute the Hessian using hess.

		is_octave = exist('OCTAVE_VERSION', 'builtin');
		if ~is_octave
			narg = nargin(problem.hess);
		else
			narg = 3;
		end

        % Check whether the hess function wants to deal with the store
        % structure or not.
        switch narg
            case 2
                hess = problem.hess(x, d);
            case 3
                % Obtain, pass along, and save the store structure
                % associated to this point.
                store = getStore(problem, x, storedb);
                [hess store] = problem.hess(x, d, store);
                storedb = setStore(problem, x, storedb, store);
            otherwise
                up = MException('manopt:getHessian:badhess', ...
                    'hess should accept 2 or 3 inputs.');
                throw(up);
        end

    elseif isfield(problem, 'ehess') && canGetEuclideanGradient(problem)
    %% Compute the Hessian using ehess.

        % We will need the Euclidean gradient for the conversion from the
        % Euclidean Hessian to the Riemannian Hessian.
        [egrad, storedb] = getEuclideanGradient(problem, x, storedb);

		is_octave = exist('OCTAVE_VERSION', 'builtin');
		if ~is_octave
			narg = nargin(problem.ehess);
		else
			narg = 3;
		end

        % Check whether the ehess function wants to deal with the store
        % structure or not.
        switch narg
            case 2
                ehess = problem.ehess(x, d);
            case 3
                % Obtain, pass along, and save the store structure
                % associated to this point.
                store = getStore(problem, x, storedb);
                [ehess store] = problem.ehess(x, d, store);
                storedb = setStore(problem, x, storedb, store);
            otherwise
                up = MException('manopt:getHessian:badehess', ...
                    'ehess should accept 2 or 3 inputs.');
                throw(up);
        end

        % Convert to the Riemannian Hessian
        hess = problem.M.ehess2rhess(x, egrad, ehess, d);

    else
    %% Attempt the computation of an approximation of the Hessian.

        [hess, storedb] = getApproxHessian(problem, x, d, storedb);

    end

end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function Obj = SOFAgetConventions(sofaconventions,flags,version)
%SOFAgetConventions - Get an empty SOFA object for a given convention
%   Usage: List = SOFAgetConventions()
%          Obj = SOFAgetConventions(conv)
%
%    List = SOFAgetConventions() returns a list with supported conventions.
%
%    Obj = SOFAgetConventions(conv) returns a SOFA object Obj
%    with the predefined mandatory metadata and data for the SOFA convention conv.
%    Obj is empty if conv is not supported. The latest version of the convention
%    conv is provided.
%
%    Obj = SOFAgetConventions(conv, flags) returns only selected
%    metadata with flags encoding the following:
%      'm' : mandatory metadata only (default)
%      'r' : readonly metadata only
%      'a' : all metadata
%
%    SOFAgetConventions(conv, ver) returns the version ver of the convention conv.
%
%    SOFAgetConventions(conv, flag, ver) returns the version ver and selected
%    metadata of the convention conv.

% #Author: Piotr Majdak
% #Author: Michael Mihocic: doc & header documentation updated (28.10.2021)
%
% SOFA Toolbox - function SOFAgetConventions
% Copyright (C) Acoustics Research Institute - Austrian Academy of Sciences
% Licensed under the EUPL, Version 1.2 or - as soon they will be approved by the European Commission - subsequent versions of the EUPL (the "License")
% You may not use this work except in compliance with the License.
% You may obtain a copy of the License at: https://joinup.ec.europa.eu/software/page/eupl
% Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing  permissions and limitations under the License.

%% Cache
persistent AllObj;

%% If flags not provided, return the conventions with all metadata
if ~exist('version','var')
  if ~exist('flags','var')
    flags='m'; % no flags, no version --> flags = a, version = any.
  else
    switch flags % no version, check if flags is flags or version
      case {'a', 'm', 'r'}
        % flags is flags --> do nothing
      otherwise
        version=flags;  % flags is version --> copy to version
        flags='m';
    end
  end
end

%% If sofaconventions not provided, return the list with supported conventions
if ~exist('sofaconventions','var')
  p=mfilename('fullpath');
  d=dir([p(1:length(p)-length(mfilename)) 'conventions' filesep '*_m_*.mat']);
  Obj={};
  for ii=1:length(d)
    dn=d(ii).name;
    idx=strfind(dn,'_');
    Obj{ii,1}=dn(1:idx(1)-1);
  end
  Obj=unique(Obj);
  AllObj=[];
	return;
end

%% Load cached object

if isfield(AllObj,flags)
  if isfield(AllObj.(flags).Obj,'GLOBAL_SOFAConventions')
    if strcmp(AllObj.(flags).Obj.GLOBAL_SOFAConventions,sofaconventions)
      if ~exist('version','var')
        Obj=AllObj.(flags).Obj; % return cached convention object
      else
        if strcmp(AllObj.(flags).Obj.GLOBAL_SOFAConventionsVersion,version)
          Obj=AllObj.(flags).Obj; % return cached convention object
        end
      end
    end
  end
end

if ~exist('Obj','var')
        % cached object not loaded yet
  p=mfilename('fullpath');
  if exist('version','var')
      % load a specific version but do not store in the cache
    if ~isempty(dir([p(1:length(p)-length(mfilename)) 'conventions' filesep sofaconventions '_' flags '_' version '.mat']))
      load([p(1:length(p)-length(mfilename)) 'conventions' filesep sofaconventions '_' flags '_' version '.mat']);
    else
      error(['Convention ' sofaconventions ' with the version ' version ' not found.']);
    end
  else
    allver=dir(fullfile(p(1:length(p)-length(mfilename)), 'conventions' , [sofaconventions '_' flags '_*.mat']));
    if isempty(allver)
      warning(['Convention ' sofaconventions ' not found.']');
      Obj=[];
    else
      vermax='0.0';
      idxmax=0;
      for ii=1:length(allver)
        x=load([p(1:length(p)-length(mfilename)) 'conventions' filesep allver(ii).name]);
        if compareversions(x.Obj.GLOBAL_SOFAConventionsVersion,vermax)>0
          vermax=x.Obj.GLOBAL_SOFAConventionsVersion;
          idxmax=ii;
        end
      end
        % store in cache
      AllObj.(flags)=load([p(1:length(p)-length(mfilename)) 'conventions' filesep allver(idxmax).name]);
      Obj=AllObj.(flags).Obj; % return the cached version
    end
  end
end

%% Overwrite some special fields
if isfield(Obj,'GLOBAL_DateCreated'), Obj.GLOBAL_DateCreated=datestr(now,SOFAdefinitions('dateFormat')); end
if isfield(Obj,'GLOBAL_APIVersion'), Obj.GLOBAL_APIVersion=SOFAgetVersion; end
if isfield(Obj,'GLOBAL_APIName'), Obj.GLOBAL_APIName=SOFAdefinitions('APIName'); end

end

function res=compareversions(ver1, ver2)

  [~,maj1,min1]=fileparts(ver1);
  [~,maj2,min2]=fileparts(ver2);
  maj1=str2num(maj1);
  min1=str2num(min1(2:end));
  maj2=str2num(maj2);
  min2=str2num(min2(2:end));
  if maj1>maj2
     res=1;   % ver1 > ver2
  elseif maj1==maj2 && min1>min2
     res=1;   % ver1 > ver2
  elseif maj2==maj2 && min1==min2
     res=0;   % ver1 == ver2
  else
     res=-1;  % ver < ver2
  end
end

author: sofacoustics

repo: SOFAtoolbox


text:
function msg_out = CreateAppend16BitCRC(msg_no_zeros)

valueCRC = 65535;
genPoly = 4129;

msg_in = [msg_no_zeros 0 0];
for i1 = 1:length(msg_in)
    for i2 = 1:8
        b = mod(floor(msg_in(i1)/(2^(8-i2))),2);
        valueCRCsh1 = bitsll(valueCRC,1);
        valueCRCadd1 = bitor(valueCRCsh1,b);
        if floor(valueCRCadd1/2^16) == 1
            valueCRC = bitxor(valueCRCadd1,genPoly);
        else
            valueCRC = valueCRCadd1;
        end
        valueCRC = mod(valueCRC,2^16);
        2;
    end
end

%msg_out = [msg_in mod(valueCRC,2^8) mod(floor(valueCRC/2^8),2^8)];
msg_out = [msg_no_zeros mod(floor(valueCRC/2^8),2^8) mod(valueCRC,2^8)];