text:
function xyz = bvh2xyz(skel, channels, noOffset)

% BVH2XYZ Compute XYZ values given structure and channels.
% FORMAT
% DESC Computes X, Y, Z coordinates given a BVH skeleton structure and
% an associated set of channels.
% ARG skel : a skeleton for the bvh file.
% ARG channels : the channels for the bvh file.
% ARG noOffset : don't add the offset in.
% RETURN xyz : the point cloud positions for the skeleton.
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2008, 2012
%
% SEEALSO : acclaim2xyz, skel2xyz

% MOCAP

if nargin< 3
  noOffset = false;
end
for i = 1:length(skel.tree)
  if ~isempty(skel.tree(i).posInd)
    xpos = channels(skel.tree(i).posInd(1));
    ypos = channels(skel.tree(i).posInd(2));
    zpos = channels(skel.tree(i).posInd(3));
  else
    xpos = 0;
    ypos = 0;
    zpos = 0;
  end
  xyzStruct(i) = struct('rotation', [], 'xyz', []);
  if nargin < 2 | isempty(skel.tree(i).rotInd)
    xangle = 0;
    yangle = 0;
    zangle = 0;
  else
    xangle = deg2rad(channels(skel.tree(i).rotInd(1)));
    yangle = deg2rad(channels(skel.tree(i).rotInd(2)));
    zangle = deg2rad(channels(skel.tree(i).rotInd(3)));
  end
  thisRotation = rotationMatrix(xangle, yangle, zangle, skel.tree(i).order);
  thisPosition = [xpos ypos zpos];
  if ~skel.tree(i).parent
      xyzStruct(i).rotation = thisRotation;
      xyzStruct(i).xyz = thisPosition + skel.tree(i).offset;
  else
      if ~noOffset
          thisPosition = skel.tree(i).offset + thisPosition;
      end
      xyzStruct(i).xyz = ...
          thisPosition*xyzStruct(skel.tree(i).parent).rotation ...
          + xyzStruct(skel.tree(i).parent).xyz;
      xyzStruct(i).rotation = thisRotation*xyzStruct(skel.tree(i).parent).rotation;

  end
end
xyz = reshape([xyzStruct(:).xyz], 3, length(skel.tree))';

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [rank,WT]=rankSetPartition(q,WT)
%%RANKSETPARTITION Return the set partition with the given rank in the
%           lexicographic ordering of all set parititions having length(n).
%           The partition is a length-n vector q that specifies which
%           partition each item belongs to. The order of the items in the
%           partitions is not important.
%
%INPUTS: q The nX1 set partition vector that should be ranked. Paritions
%          are numbered in increasing order starting from 1. That is, the
%          first occurence of partition numbered i in q is always before
%          the first occurences of partitions numbered j>1.
%       WT The is, optionally, the Williamson table. It is the table T in
%          [1]. The table is a function of length(q). If omitted, then the
%          table will be computed and returned so that futue calls to the
%          function can use
%
%OUTPUTS: rank The rank of the set partition in lexicographic order.
%              0<=rank<BellNumber(n).
%           WT The Williamson table for a length-n set partition. Passing
%              this to rankSetPartition on subsequent calls will speed up
%              the algorithm.
%
%This function implements the ranking algorithm of [1]. See
%getNextSetPartition for more information on set partitions.
%
%REFERENCES:
%[1] S. G. Williamson, "Ranking algorithms for lists of partitions," SIAM
%    Journal on Computing, vol. 5, no. 4, pp. 602-617, 1976.
%
%January 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

n=length(q);

if(nargin<3||isempty(WT))
    %Compute the Williamson triangle
    WT=zeros(n,n);
    WT(1,1:n)=1;
    for v=2:n
        for mu=1:(n-v+1)
            WT(v,mu)=mu*WT(v-1,mu)+WT(v-1,mu+1);
        end
    end
end

q=q-1;

rank=0;

m=cummax(q);

for t=0:(n-2)
    idx=n-t;

    rank=rank+q(idx)*WT(t+1,m(idx-1)+1);
end

rank=rank+1;

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%% Runge-Kutta 2
%
% for solving
%
% $\frac{dv}{dt}=-\alpha (t)v+\beta (t)$
%
% where $\alpha (t) = \frac{3t}{1-t}$ and $\beta (t)=2(1-t)^3 e^{-t}$
%
% Assuming v(0)=1.0 for the period 0<t<15

clc
clear
%close all

%% Discretization of t
ti=0;    % Initial Time
tf=150;   % Final Time

% Time step
% h = 0.2;
% h = 0.8;
% h = 1.1;
 h = 0.68;

% Discretization of t
i=1;
t(1)=ti;
while t(i)<tf
        i=i+1;
        t(i)=t(i-1)+h;
end
t(i)=tf; clear i;
n=length(t);

%% Discretization of v

v = zeros(n,1);
alpha = zeros(n,1);
beta  = zeros(n,1);
t2 = zeros(n,1);
alpha2 = zeros(n,1);
beta2  = zeros(n,1);

for i=1:n-1
    t2(i)     = t(i)+1/2*h;
    % compute variables alpha(t) and beta(t)
    alpha(i)  = 3*t(i)/(1+t(i));
    alpha2(i) = 3*t2(i)/(1+t2(i)); % alpha(t + 1/2h)
    beta(i)   = 2*(1+t(i))^3*exp(-t(i));
    beta2(i)  = 2*(1+t2(i))^3*exp(-t2(i)); % beta(t + 1/2h)
end

% initial condition (v(0))
v(1)=1;

for i=1:n-1
    % main discrete equation
    k1 = h*(-alpha(i)*v(i)+beta(i));
    k2 = h*(-alpha2(i)*(v(i)+1/2*k1)+beta2(i));
    v(i+1) = v(i)+k2;
end

%% Exact solution
v2 = zeros(n,1);
for i=1:n
     v2(i) = exp(-t(i))*(1+t(i))^3;
end
error=v2-v;
figure
plot(t,v,'.',t,v2,'-')
%figure
%plot(t,error,'.')

repo: Aero-matlab


text:
function Fx = computeFx(x)
    Fx = [0 0 -x(4)*sin(x(3))-x(5)*cos(x(3)) cos(x(3)) -sin(x(3)) 0;...
          0 0 x(4)*cos(x(3))-x(5)*sin(x(3))  sin(x(3))  cos(x(3)) 0;...
          0 0               0                   0           0     1;...
          0 0               0                   0           0     0;...
          0 0               0                   0           0     0;...
          0 0               0                   0           0     0];
end

author: ccalas

repo: mpc

save_path: github-repos/MATLAB/ccalas-mpc


text:
function [varargout]=regionTriMesh3D(varargin)

% function [F,V,boundaryInd]=regionTriMesh3D(regionCell,pointSpacing,resampleCurveOpt,interpMethod)
% ------------------------------------------------------------------------
% This function creates a 3D triangulation for the region specified in the
% variable regionCell. The mesh aims to obtain a point spacing as defined
% by the input pointSpacing.
% The function output contains the triangular faces in F, the vertices in V
% and the per triangle region indices in regionInd. By setting plotOn to 0
% or 1 plotting can be switched on or off.
%
% More on the specification of the region:
% The input variable regionCell is a cell array containing all the boundary
% curves, e.g. for a two curve region 1 we would have something like
% regionSpec{1}={V1,V2} where V1 and V2 are the boundary curves. Multiple
% curves may be given here. The first curve should form the outer boundary
% of the entire region, the curves that follow should define holes inside
% this boundary and the space inside them is therefore not meshed.
%
% See also: regionTriMesh2D
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2013/11/21: Created
% 2013/11/21: Added varargin input parsing
%------------------------------------------------------------------------

%% Parse input

switch nargin
    case 1
        regionCell=varargin{1};
        pointSpacing=[];
        resampleCurveOpt=[];
        interpMethod=[];
    case 2
        regionCell=varargin{1};
        pointSpacing=varargin{2};
        resampleCurveOpt=[];
        interpMethod=[];
    case 3
        regionCell=varargin{1};
        pointSpacing=varargin{2};
        resampleCurveOpt=varargin{3};
        interpMethod=[];
    case 4
        regionCell=varargin{1};
        pointSpacing=varargin{2};
        resampleCurveOpt=varargin{3};
        interpMethod=varargin{4};
end

if isempty(pointSpacing)
    D=zeros(1,numel(regionCell));
    for q=1:1:numel(regionCell)
        D(q)=mean(diff(polyCurveLength(regionCell{1})));
    end
    pointSpacing=mean(D);
end

if isempty(resampleCurveOpt)
    resampleCurveOpt=0;
end

if isempty(interpMethod)
    interpMethod='natural';
end

%% Parse 3D regions and convert to planar 2D

% TO DO: get sizes and do it properly, hint cellfun(@(x) size(x,1),regionCell)
npSets=cellfun(@(x) size(x,1),regionCell);
npTotal=sum(npSets);

%Collect all points
Vt=zeros(npTotal,3);
startInd=1;
endInd=npSets(1);
for qCurve=1:1:numel(regionCell)
    Vs=regionCell{qCurve};
    Vt(startInd:endInd,:)=Vs;
    startInd=startInd+npSets(qCurve);
    if qCurve<numel(regionCell)
        endInd=startInd+npSets(qCurve+1)-1;
    end
end

Vt_mean=mean(Vt,1);
Vtc=Vt-Vt_mean(ones(size(Vt,1),1),:);

%Get orientation
[R_fit]=pointSetPrincipalDir(Vtc);

%Create planar 2D sets
Vtcr=Vt;
startInd=1;
endInd=npSets(1);
regionCell2D=regionCell;
for qCurve=1:1:numel(regionCell)
    Vs=Vt(startInd:endInd,:);
    Vs=Vs-Vt_mean(ones(size(Vs,1),1),:);
    Vs=(R_fit'*Vs')';
    regionCell2D{qCurve}=Vs(:,[1 2]);
    Vtcr(startInd:endInd,:)=Vs;
    startInd=startInd+npSets(qCurve);
    if qCurve<numel(regionCell)
        endInd=startInd+npSets(qCurve+1)-1;
    end
end

%% Use regionTriMesh2D to mesh region in 2D

[F,V,boundaryInd]=regionTriMesh2D(regionCell2D,pointSpacing,resampleCurveOpt,0);

%% Interpolate to obtain Z coordinates

interpFunc = scatteredInterpolant(Vtcr(:,[1 2]),Vtcr(:,3),interpMethod);
V(:,3)=interpFunc(V(:,[1 2]));
V=(R_fit*V')';
V=V+Vt_mean(ones(size(V,1),1),:);

%% Check normal directions against input curve(s)

% V_test=regionCell{1};
% V_test=V_test(1:2,:);

% [~,e]=minDist(V_test,V);
% E=patchEdges(F,1);
% if ~any(E(:,1)==e(1) & E(:,2)==e(2))
%    F=fliplr(F); %Invert if edge does not seem member of faces
% end

%% Collect ouput
varargout{1}=F;
varargout{2}=V;
varargout{3}=boundaryInd;

%%
% _*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:
% rigid registration of frames with offsets ds
function dreg = register_movie(data, ops, ds)

orig_class = class(data);

% compute maximum batch size for GPU
[h, w, nFrames] = size(data);
if getOr(ops, 'useGPU', 0)
    nFramesPerBatch = getBatchSize(w*h);
else
    nFramesPerBatch = 1000; % if not on GPU, should have plenty of RAM available
end
%%
nFramesPerBatch = 8;

nBatches = ceil(nFrames/nFramesPerBatch);
startFrame = 1:nFramesPerBatch:nFrames;
endFrame = min(startFrame+nFramesPerBatch-1, nFrames);
dreg = zeros(size(data), orig_class);

for iBatch = 1:nBatches
    idx = startFrame(iBatch):endFrame(iBatch);
    if ops.useGPU
        dataBatch = gpuArray(single(data(:,:,idx)));
    else
        dataBatch = data(:,:,idx);
    end
    [Ly, Lx, NT] = size(dataBatch);

    Ny = ifftshift([-fix(Ly/2):ceil(Ly/2)-1]);
    Nx = ifftshift([-fix(Lx/2):ceil(Lx/2)-1]);
    [Nx,Ny] = meshgrid(Nx,Ny);
    Nx = Nx / Lx;
    Ny = Ny / Ly;

    if ops.useGPU
        dsBatch = gpuArray(permute(ds(idx, :), [3, 2, 1]));
        Nx = gpuArray(single(Nx));
        Ny = gpuArray(single(Ny));
    else
        dsBatch = ds(idx, :);
    end

    if ops.useGPU % do it batch-by-batch
        dph         = 2*pi*(bsxfun(@times, dsBatch(1,1,:), Ny) + ...
            bsxfun(@times, dsBatch(:,2,:), Nx));
        fdata       = fft2(dataBatch);
        dregBatch   = gather_try(real(ifft2(fdata .* exp(1i * dph))));
    else % do it frame-by-frame
        dregBatch = zeros(size(dataBatch), orig_class);
        for i = 1:NT
            dph         = 2*pi*(dsBatch(i,1)*Ny + dsBatch(i,2)*Nx);
            fdata       = fft2(single(dataBatch(:,:,i)));
            dregBatch(:,:,i) = real(ifft2(fdata .* exp(1i * dph)));
        end
    end

    dreg(:,:,idx) = dregBatch;
end

author: cortex-lab


text:
function xyz = skel2xyz(skel, channels)

% SKEL2XYZ Compute XYZ values given skeleton structure and channels.
% FORMAT
% DESC Computes X, Y, Z coordinates given a BVH or acclaim skeleton
% structure and an associated set of channels.
% ARG skel : a skeleton for the bvh file.
% ARG channels : the channels for the bvh file.
% RETURN xyz : the point cloud positions for the skeleton.
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% SEEALSO : acclaim2xyz, bvh2xyz

% MOCAP

fname = str2func([skel.type '2xyz']);
xyz = fname(skel, channels);


text:
classdef ObjectiveFunctionEnum < matlab.mixin.SetGet
    %ObjectiveFunctionEnum Summary of this class goes here
    %   Detailed explanation goes here
    %%%% See ConstraintEnum!
    enumeration
        NoObjectiveFunction('Satisfy Constraints Only','NoOptimizationObjectiveFcn')
        MinDistToBody('Minimize Distance to Body','MinDistanceToBodyObjectiveFcn')
        MaximizeVehicleMass('Maximize Vehicle Mass','MaximizeLaunchVehicleMassObjectiveFcn')
    end

    properties
        name char = '';
        class char = '';
    end

    methods
        function obj = ObjectiveFunctionEnum(name, class)
            obj.name = name;
            obj.class = class;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('ObjectiveFunctionEnum');
            listBoxStr = {m.name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('ObjectiveFunctionEnum');
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [ind, enum] = getIndForClass(objFcnObject)
            m = enumeration('ObjectiveFunctionEnum');
            ind = find(ismember({m.class},class(objFcnObject)),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function varargout = ImageGUI(varargin)
% IMAGEGUI MATLAB code for ImageGUI.fig
%      IMAGEGUI, by itself, creates a new IMAGEGUI or raises the existing
%      singleton*.
%
%      H = IMAGEGUI returns the handle to a new IMAGEGUI or the handle to
%      the existing singleton*.
%
%      IMAGEGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in IMAGEGUI.M with the given input arguments.
%
%      IMAGEGUI('Property','Value',...) creates a new IMAGEGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before ImageGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to ImageGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help ImageGUI

% Copyright 2011-2013 The MathWorks, Inc.

% Last Modified by GUIDE v2.5 31-Mar-2011 16:55:25

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @ImageGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @ImageGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before ImageGUI is made visible.
function ImageGUI_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to ImageGUI (see VARARGIN)

% Choose default command line output for ImageGUI
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes ImageGUI wait for user response (see UIRESUME)
% uiwait(handles.figure1);

movegui(handles.figure1, 'center');

% --- Outputs from this function are returned to the command line.
function varargout = ImageGUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

findBallFcn(handles.im);

% --------------------------------------------------------------------
function uipushtool1_ClickedCallback(hObject, eventdata, handles)
% hObject    handle to uipushtool1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

[filename, pathname] = uigetfile('*.jpg', 'Select image file');

im = imread(fullfile(pathname, filename));

imshow(im);

handles.im = im;

guidata(hObject, handles);

author: Sable

repo: mcbench-benchmarks


text:
classdef LIRCMOP5 < PROBLEM
% <multi> <real> <large/none> <constrained>
% Constrained benchmark MOP with large infeasible regions

%------------------------------- Reference --------------------------------
% Z. Fan, W. Li, X. Cai, H. Huang, Y. Fang, Y. You, J. Mo, C. Wei, and E.
% Goodman, An improved epsilon constraint-handling method in MOEA/D for
% CMOPs with large infeasible regions, Soft Computing, 2019, 23:
% 12491-12510.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Wenji Li

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values and constraint violations
        function Population = Evaluation(obj,varargin)
            X = varargin{1};
            X = max(min(X,repmat(obj.upper,size(X,1),1)),repmat(obj.lower,size(X,1),1));
            [popsize,variable_length] = size(X);
            sum1 = zeros(popsize,1);
            sum2 = zeros(popsize,1);
            for j = 2 : variable_length
                if mod(j,2) == 1
                    sum1 = sum1+(X(:,j)-sin((0.5*j/variable_length*pi)*X(:,1))).^2;
                else
                    sum2 = sum2+(X(:,j)-cos((0.5*j/variable_length*pi)*X(:,1))).^2;
                end
            end
            gx          = 0.7057;
            PopObj(:,1) = X(:,1)+10*sum1+gx;
            PopObj(:,2) = 1-X(:,1).^0.5+10.*sum2+gx;
            Population  = SOLUTION(X,PopObj,Constraint(PopObj),varargin{2:end});
            obj.FE      = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - sqrt(R(:,1));
            R      = R + 0.7057;
            R(any(Constraint(R)>0,2),:) = [];
        end
        %% Generate the feasible region
        function R = GetPF(obj)
            [x,y] = meshgrid(linspace(0.7057,5,400));
            fes   = all(Constraint([x(:),y(:)])<=0,2);
            z     = nan(size(x));
            z(reshape(fes,size(z)) & sqrt(x-0.7057)+y>=1.7057) = 0;
            R = {x,y,z};
        end
    end
end

function PopCon = Constraint(PopObj)
    p     = [1.6,2.5];
    q     = [1.6,2.5];
    a     = [2,2];
    b     = [4,8];
    r     = 0.1;
    theta = -0.25 * pi;
    for k = 1 : 2
        PopCon(:,k) = r - ((PopObj(:,1)-p(k))*cos(theta)-(PopObj(:,2)-q(k))*sin(theta)).^2/(a(k)^2) - ...
                      ((PopObj(:,1)-p(k))*sin(theta)+(PopObj(:,2)-q(k))*cos(theta)).^2/(b(k)^2);
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function	[xdata,parm] = normalize_data(xdata,norm_mode,parm)
% Normalize data
% --- Usage
%  1. Normalize data and set mean & norm to 'parm'
%  [xdata, parm] = normalize_data(xdata)
%  [xdata, parm] = normalize_data(xdata, norm_mode)
%
%  2.  Normalize data using 'parm'
%  [xdata] = normalize_data(xdata, norm_mode, parm)
%
% --- Input
% xdata : [Xdim x Nsample x Ntrial]
% --- Optional input
% norm_mode : normalization mode
%           = 0: zero mean & std is calculated
%           = 1: zero mean & normalized [Default]
%           = 2: normalized (mean is unchanged)
%
% parm  : if the following fields exist, these values are used for mean & std
%    parm.xmean = mean of xdata
%    parm.xnorm = std (root of variance) of xdata
%
% --- Output
%    parm.xmean = mean of xdata                    [Xdim x 1]
%    parm.xnorm = std (root of variance) of xdata  [Xdim x 1]
%  if norm_mode == 0,
%     xdata : zero mean
%  if norm_mode = 1
%     xdata : zero mean & normalized
%  if norm_mode = 2
%     xdata : normalized  (mean is unchanged)
%
%  if 'parm.xmean' & 'parm.xnorm' are given,
%     these values are used as mean & std
%
% 2007/1/27 M. Sato
% 2007-4-8 Modified by M. Sato
% 2009-11-2 Modified by M. Sato, norm_mode = 2 is added

if ~exist('norm_mode','var'),
	norm_mode = 1;
end
if ~exist('parm','var'),
	parm = [];
end

[Xdim,Ndata,Ntrial] = size(xdata);

xdata = reshape(xdata, [Xdim, Ndata*Ntrial]);

% Mean
if ~isfield(parm,'xmean') | isempty(parm.xmean)
	parm.xmean = sum(xdata,2)/(Ndata*Ntrial);
end

xdim = length(parm.xmean);

if xdim > Xdim,
	parm.xmean = parm.xmean(1:Xdim);
elseif xdim < Xdim,
	parm.xmean = repmat(parm.xmean,[Xdim/xdim 1]);
end

% Zero mean
if norm_mode < 2
% 	xdata = xdata - repmat(parm.xmean, [1 Ndata]);
	xdata = repadd( double(xdata) , - double(parm.xmean) );
end

% Variance
if ~isfield(parm,'xnorm') | isempty(parm.xnorm)
	parm.xnorm = sqrt(sum(xdata.^2,2)/((Ndata-1)*Ntrial)); % modified by TH: use unbiased variance
end

xdim = length(parm.xnorm);

if xdim > Xdim,
	parm.xnorm = parm.xnorm(1:Xdim);
elseif xdim < Xdim,
	parm.xnorm = repmat(parm.xnorm,[Xdim/xdim 1]);
end

% Normalize amplitude
if norm_mode > 0
% 	xdata = xdata./repmat(parm.xnorm, [1 Ndata]);
	xdata = repmultiply( double(xdata) , 1./double(parm.xnorm) );
end

xdata = reshape(xdata, [Xdim, Ndata, Ntrial]);

author: KamitaniLab

repo: GenericObjectDecoding


text:
function [g,a,fc,L,info]=audfilters(fs,Ls,varargin)
%AUDFILTERS Generates filters equidistantly spaced on auditory frequency scales
%   Usage:  [g,a,fc,L]=audfilters(fs,Ls);
%           [g,a,fc,L]=audfilters(fs,Ls,...);
%
%   Input parameters:
%      fs    : Sampling rate (in Hz).
%      Ls    : Signal length.
%   Output parameters:
%      g     : Cell array of filters.
%      a     : Downsampling rate for each channel.
%      fc    : Center frequency of each channel.
%      L     : Next admissible length suitable for the generated filters.
%
%   `[g,a,fc,L]=audfilters(fs,Ls)` constructs a set of filters *g* that are
%   equidistantly spaced on a perceptual frequency scale (see |freqtoaud|) between
%   0 and the Nyquist frequency. The filter bandwidths are proportional to the
%   critical bandwidth of the auditory filters |audfiltbw|. The filters are intended
%   to work with signals with a sampling rate of *fs*. The signal length *Ls* is
%   mandatory, since we need to avoid too narrow frequency windows.
%
%   By default the ERB scale is chosen but other frequency scales are
%   possible. Currently supported scales are 'erb', 'erb83', 'bark', 'mel'
%   and 'mel1000', and can be changed by passing the associated string as
%   an optional parameter. See |freqtoaud| for more information on the
%   supported frequency scales.
%
%   By default, a Hann window shape is chosen as prototype frequency
%   response for all filters. The prototype frequency response can be
%   changed by passing any of the window types from |firwin| or |freqwin|
%   as an optional parameter.
%
%   `[g,a,fc,L]=audfilters(fs,Ls,fmin,fmax)` constructs a set of filters
%   between *fmin* and *fmax*. The filters are equidistantly spaced on the
%   selected frequency scale. One additional filter will be positioned at
%   the 0 and Nyquist frequencies each, so as to cover the full range of
%   positive frequencies.
%   The values of *fmin* and *fmax* can be instead specified using a
%   key/value pair as::
%
%       [g,a,fc,L]=audfilters(fs,Ls,...,'fmin',fmin,'fmax',fmax)
%
%   Default values are *fmin=0* and *fmax=fs/2*.
%
%   For more details on the construction of the filters, please see the
%   given references.
%
%   Downsampling factors
%   --------------------
%
%   The integer downsampling rates of the channels must all divide the
%   signal length, |filterbank| will only work for input signal lengths
%   being multiples of the least common multiple of the downsampling rates.
%   See the help of |filterbanklength|.
%   The fractional downsampling rates restrict the filterbank to a single
%   length *L=Ls*.
%
%   `[g,a]=audfilters(...,'regsampling')` constructs a non-uniform
%   filterbank with integer subsampling factors.
%
%   `[g,a]=audfilters(...,'uniform')` constructs a uniform filterbank
%   where the integer downsampling rate is the same for all the channels. This
%   results in most redundant representation which produces nice plots.
%
%   `[g,a]=audfilters(...,'fractional')` constructs a filterbank with
%   fractional downsampling rates *a*.
%   This results in the least redundant system.
%
%   `[g,a]=audfilters(...,'fractionaluniform')` constructs a filterbank with
%   fractional downsampling rates *a*, which are uniform for all filters
%   except the "filling" low-pass and high-pass filters which can have different
%   fractional downsampling rates. This is useful when uniform subsampling
%   and low redundancy at the same time are desirable.
%
%   Additional parameters
%   ---------------------
%
%   `audfilters` accepts the following optional parameters:
%
%     'spacing',b        Specify the spacing between the filters, measured in
%                        scale units. Default value is *b=1* for the scales
%                        'erb', 'erb83' and 'bark'; the default is *b=100* for
%                        'mel' and 'mel1000'.
%
%     'bwmul',bwmul      Bandwidth of the filters relative to the bandwidth
%                        returned by |audfiltbw|. Default value is *bwmul=1* for
%                        the scales 'erb', 'erb83' and 'bark'; the default is
%                        *b=100* for 'mel' and 'mel1000'.
%
%     'redmul',redmul    Redundancy multiplier. Increasing the value of this
%                        will make the system more redundant by lowering the
%                        channel downsampling rates. It is only used if the
%                        filterbank is a non-uniform filterbank. Default
%                        value is *1*. If the value is less than one, the
%                        system may no longer be painless.
%
%     'redtar',redtar    Target redundancy. The downsampling factors will be
%                        adjusted to achieve a redundancy as close as possible
%                        to 'redtar'.
%
%     'M',M              Specify the total number of filters between *fmin* and
%                        *fmax*. If this parameter is specified, it overwrites the
%                        `'spacing'` parameter.
%
%     'symmetric'        Create filters that are symmetric around their centre
%                        frequency. This is the default.
%
%     'warped'           Create asymmetric filters that are symmetric on the
%                        auditory scale.
%
%     'complex'          Construct a filterbank that covers the entire
%                        frequency range instead of just the positive
%                        frequencies this allows the analysis of complex
%                        valued signals.
%
%     'nosubprec'        Disable subsample window positions.
%
%     'trunc_at'         When using a prototype defined in |freqwin|, a hard
%                        thresholding of the filters at the specified threshold
%                        value is performed to reduce their support size.
%                        The default value is *trunc_at=10e-5*. When no
%                        truncation is desired, *trunc_at=0* should be chosen.
%                        This value is ignored when a prototype shape from
%                        |firwin| was chosen.
%
%     'min_win',min_win  Minimum admissible window length (in samples).
%                        Default is *4*. This restrict the windows not
%                        to become too narrow when *L* is low.
%
%   Examples:
%   ---------
%
%   In the first example, we construct a highly redudant uniform
%   filterbank on the ERB scale and visualize the result:::
%
%     [f,fs]=greasy;  % Get the test signal
%     [g,a,fc,L]=audfilters(fs,length(f),'uniform','M',100);
%     c=filterbank(f,g,a);
%     plotfilterbank(c,a,fc,fs,90,'audtick');
%
%   In the second example, we construct a non-uniform filterbank with
%   fractional sampling that works for this particular signal length, and
%   test the reconstruction. The plot displays the response of the
%   filterbank to verify that the filters are well-behaved both on a
%   normal and an ERB-scale. The second plot shows frequency responses of
%   filters used for analysis (top) and synthesis (bottom). :::
%
%     [f,fs]=greasy;  % Get the test signal
%     L=length(f);
%     [g,a,fc]=audfilters(fs,L,'fractional');
%     c=filterbank(f,{'realdual',g},a);
%     r=2*real(ifilterbank(c,g,a));
%     norm(f-r)
%
%     % Plot the response
%     figure(1);
%     subplot(2,1,1);
%     R=filterbankresponse(g,a,L,fs,'real','plot');
%
%     subplot(2,1,2);
%     semiaudplot(linspace(0,fs/2,L/2+1),R(1:L/2+1));
%     ylabel('Magnitude');
%
%     % Plot frequency responses of individual filters
%     gd=filterbankrealdual(g,a,L);
%     figure(2);
%     subplot(2,1,1);
%     filterbankfreqz(gd,a,L,fs,'plot','linabs','posfreq');
%
%     subplot(2,1,2);
%     filterbankfreqz(g,a,L,fs,'plot','linabs','posfreq');
%
%
%   See also: filterbank, ufilterbank, ifilterbank, ceil23
%
%   References: ltfatnote027 nehobaprpide18

% Authors: Peter L. Søndergaard (original 'erbfilters' function)
% Modified by: Thibaud Necciari, Nicki Holighaus
% Comments updated by: Nicki Holighaus (09.05.22)

% Date: 16.12.16

complainif_notenoughargs(nargin,2,upper(mfilename));
complainif_notposscalar(fs,'fs',upper(mfilename));
complainif_notposint(Ls,'Ls',upper(mfilename));

firwinflags=getfield(arg_firwin,'flags','wintype');
freqwinflags=getfield(arg_freqwin,'flags','wintype');

definput.flags.wintype = [ firwinflags, freqwinflags];
definput.keyvals.M=[];
definput.keyvals.redmul=1;
definput.keyvals.min_win = 4;
definput.keyvals.bwmul=[];
definput.keyvals.spacing=[];
definput.keyvals.trunc_at=10^(-5);
definput.keyvals.fmin=0;
definput.keyvals.fmax=fs/2;
definput.keyvals.redtar=[];
definput.flags.subprec={'subprec','nosubprec'};
definput.flags.audscale={'erb','erb83','bark','mel','mel1000'};
definput.flags.warp     = {'symmetric','warped'};
definput.flags.real     = {'real','complex'};
definput.flags.sampling = {'regsampling','uniform','fractional',...
                           'fractionaluniform'};

[varargin,winCell] = arghelper_filterswinparser(definput.flags.wintype,varargin);

[flags,kv]=ltfatarghelper({'fmin','fmax'},definput,varargin);
if isempty(winCell), winCell = {flags.wintype}; end

switch flags.audscale
    case {'mel','mel1000'} % The mel scales are very fine, therefore default spacing is adjusted
        definput.keyvals.bwmul=100;
        definput.keyvals.spacing=100;
    otherwise
        definput.keyvals.bwmul=1;
        definput.keyvals.spacing=1;
end
[flags,kv]=ltfatarghelper({'fmin','fmax','redtar'},definput,varargin);

if flags.do_bark && (fs > 44100)
    error(['%s: Bark scale is not suitable for sampling rates higher than 44.1 kHz. ',...
    'Please choose another scale.'],upper(mfilename));
end

if ~isscalar(kv.bwmul) || kv.bwmul <= 0
    error('%s: bwmul must be a positive scalar.',upper(mfilename));
end

if ~isscalar(kv.redmul) || kv.redmul <= 0
    error('%s: redmul must be a positive scalar.',upper(mfilename));
end

if ~isempty(kv.redtar)
    if ~isscalar(kv.redtar) || kv.redtar <= 0
        error('%s: redtar must be a positive scalar.',upper(mfilename));
    end
end

if kv.redtar <= 1
    warning('%s: redtar is very low; the resulting system might be unstable.',upper(mfilename));
end

if kv.fmax <= kv.fmin || kv.fmin < 0 || kv.fmax > fs/2
    error('%s: fmax must be bigger than fmin and in the range [0,fs/2].',upper(mfilename));
end

if kv.trunc_at > 1 || kv.trunc_at < 0
    error('%s: trunc_at must be in range [0,1].',upper(mfilename));
end

if ~isscalar(kv.min_win) || rem(kv.min_win,1) ~= 0 || kv.min_win < 1
    error('%s: min_win must be an integer bigger or equal to 1.',upper(mfilename));
end

if ~isempty(kv.M)
    complainif_notposint(kv.M,'M',upper(mfilename));
    kv.spacing = (freqtoaud(kv.fmax,flags.audscale) - freqtoaud(kv.fmin,flags.audscale))/(kv.M-1);
end

% Construct function handle for filter prototype and determine its ERB-type
% bandwidth
[filterfunc,winbw] = helper_filtergeneratorfunc(...
                          flags.wintype,winCell,fs,kv.bwmul,kv.min_win,kv.trunc_at,...
                          flags.audscale,flags.do_subprec,flags.do_symmetric,flags.do_warped);

% Construct the AUD filterbank
fmin = max(kv.fmin,audtofreq(kv.spacing,flags.audscale));
fmax = min(kv.fmax,fs/2);

innerChanNum = floor((freqtoaud(fmax,flags.audscale)-freqtoaud(fmin,flags.audscale))/kv.spacing)+1;

fmax = audtofreq(freqtoaud(fmin,flags.audscale)+(innerChanNum-1)*kv.spacing,flags.audscale);

% Make sure that fmax < fs/2, and F_ERB(fmax) = F_ERB(fmin)+k/spacing, for
% some k.
count = 0;
while fmax >= fs/2
    count = count+1;
    fmax = audtofreq(freqtoaud(fmin,flags.audscale)+(innerChanNum-count-1)*kv.spacing,flags.audscale);
end
innerChanNum = innerChanNum-count;

if fmax <= fmin || fmin > fs/4 || fmax < fs/4
    error(['%s: Bad combination of fs, fmax and fmin.'],upper(mfilename));
end

% Center frequencies are given as equidistantly spaced points on auditory
% scale
fc=audspace(fmin,fmax,innerChanNum,flags.audscale).';
fc = [0;fc;fs/2];
M2 = innerChanNum+2;

ind = (2:M2-1)';
%% Compute the frequency support
% fsupp is measured in Hz

fsupp=zeros(M2,1);
aprecise=zeros(M2,1);

if flags.do_symmetric
    fsupp(ind)=audfiltbw(fc(ind),flags.audscale)/winbw*kv.bwmul;

    % Generate lowpass filter parameters
    fsupp(1) = 0; % Placeholder value
    % Determine border of passband
    fps0 = audtofreq(freqtoaud(fc(2),flags.audscale)+3*kv.spacing,flags.audscale);% f_{p,s}^{-}
    % Determine lowpass width
    fsupp_temp1 = audfiltbw(fps0,flags.audscale)/winbw*kv.bwmul;
    % Determine bandwidth-adapted decimation factor
    aprecise(1) = max(fs./(2*max(fps0,0)+fsupp_temp1*kv.redmul),1);

    % Generate highpass filter parameters
    fsupp(end) = 0; % Placeholder value
    % Determine border of passband
    fps0 = audtofreq(freqtoaud(fc(end-1),flags.audscale)-3*kv.spacing,flags.audscale);% f_{p,s}^{+}
    % Determine highpass width
    fsupp_temp1 = audfiltbw(fps0,flags.audscale)/winbw*kv.bwmul;
    % Determine bandwidth-adapted decimation factor
    aprecise(end) = max(fs./(2*(fc(end)-min(fps0,fs/2))+fsupp_temp1*kv.redmul),1);
else
    % fsupp_scale is measured on the selected auditory scale
    % The scaling is incorrect, it does not account for the warping (NH:
    % I do think it is correct.)
    fsupp_scale=1/winbw*kv.bwmul;

    % Convert fsupp into the correct widths in Hz, necessary to compute
    % "a" in the next if-statement
    fsupp(ind)=audtofreq(freqtoaud(fc(ind),flags.audscale)+fsupp_scale/2,flags.audscale)-...
               audtofreq(freqtoaud(fc(ind),flags.audscale)-fsupp_scale/2,flags.audscale);

    % Generate lowpass filter parameters
    fsupp(1) = 0; % Placeholder value
    % Determine border of passband
    fps0 = audtofreq(freqtoaud(fc(2),flags.audscale)+3*kv.spacing,flags.audscale);% f_{p,s}^{-}
    % Determine lowpass width
    fsupp_temp1 = audfiltbw(fps0,flags.audscale)/winbw*kv.bwmul;
    % Determine bandwidth-adapted decimation factor
    aprecise(1) = max(fs./(2*max(fps0,0)+fsupp_temp1*kv.redmul),1);

    % Generate highpass filter parameters
    fsupp(end) = 0; % Placeholder value
    % Determine border of passband
    fps0 = audtofreq(freqtoaud(fc(end-1),flags.audscale)-3*kv.spacing,flags.audscale);% f_{p,s}^{+}
    % Determine highpass width
    fsupp_temp1 = audfiltbw(fps0,flags.audscale)/winbw*kv.bwmul;
    % Determine bandwidth-adapted decimation factor
    aprecise(end) = max(fs./(2*(fc(end)-min(fps0,fs/2))+fsupp_temp1*kv.redmul),1);
end;

% Do not allow lower bandwidth than keyvals.min_win
fsuppmin = kv.min_win/Ls*fs;
for ii = 1:numel(fsupp)
    if fsupp(ii) < fsuppmin;
        fsupp(ii) = fsuppmin;
    end
end

% Find suitable channel subsampling rates
aprecise(ind)=fs./fsupp(ind)/kv.redmul;
aprecise=aprecise(:);

if any(aprecise<1)
    error('%s: Invalid subsampling rates. Decrease redmul.',upper(mfilename))
end

%% Compute the downsampling rate
if flags.do_regsampling
    % Shrink "a" to the next composite number
    a=floor23(aprecise);

    % Determine the minimal transform length
    L=filterbanklength(Ls,a);

    % Heuristic trying to reduce lcm(a)
    while L>2*Ls && ~(all(a)==a(1))
        maxa = max(a);
        a(a==maxa) = 0;
        a(a==0) = max(a);
        L = filterbanklength(Ls,a);
    end

% Determine true decimation factors from "aprecise" according to chosen
% subsampling scheme (see help above)
elseif flags.do_fractional
    L = Ls;
    N=ceil(Ls./aprecise);
    a=[repmat(Ls,M2,1),N];
elseif flags.do_fractionaluniform
    L = Ls;
    N=ceil(Ls./min(aprecise));
    a= repmat([Ls,N],M2,1);
elseif flags.do_uniform
    a=floor(min(aprecise));
    L=filterbanklength(Ls,a);
    a = repmat(a,M2,1);
end;

% Get an expanded "a" / Convert "a" to LTFAT 2-column fractional format
afull=comp_filterbank_a(a,M2,struct());

%% Compute the scaling of the filters
% Filters are scaled such that the energy of the subband coefficients
% remains approximately constant independent of the decimation factor
scal=sqrt(afull(:,1)./afull(:,2));

%% Construct the real or complex filterbank

if flags.do_real
    % Scale the first and last channels
    scal(1)=scal(1)/sqrt(2);
    scal(M2)=scal(M2)/sqrt(2);
else
    % Replicate the centre frequencies and sampling rates, except the first and
    % last
    a=[a;flipud(a(2:M2-1,:))];
    scal=[scal;flipud(scal(2:M2-1))];
    fc  =[fc; -flipud(fc(2:M2-1))];
    fsupp=[fsupp;flipud(fsupp(2:M2-1))];
    ind = [ind;numel(fc)+2-(M2-1:-1:2)'];
end;

%% Compute the filters
% This is actually much faster than the vectorized call.
g = cell(1,numel(fc));
for m=ind.'
    g{m}=filterfunc(fsupp(m),fc(m),scal(m));
end

% Generate lowpass filter
g{1} = audlowpassfilter(g(1:M2),a(1:M2,:),fc(1:M2),fs,scal(1),kv,flags);

% Generate highpass filter
g{M2} = audhighpassfilter(g(1:M2),a(1:M2,:),fc(1:M2),fs,scal(M2),kv,flags);

% Adjust the downsampling rates in order to achieve 'redtar' (see help
% above)
if ~isempty(kv.redtar)
    if flags.do_uniform
        % Compute and display redundancy for verification
        org_red = (M2-2)/a(1);
        a_new = floor(a*org_red/kv.redtar);
        scal_new = org_red/kv.redtar*ones(numel(g),1);
        % new_red = (M2-2)/a_new(1);
    else
        % The decimation factors of all filters except for lowpass and
        % highpass are adjusted proportional to their bandwidth. For
        % lowpass and highpass, only the tapered part is considered for
        % adjustment to better preserve stability. Please consult the
        % references for details.
        dk_old = a(:,1)./a(:,2);
        org_red = sum(2./dk_old(2:end-1));
        a_new = [a(1,:);[a(2:end-1,1),ceil(a(2:end-1,2)*kv.redtar/org_red)];a(end,:)];
        % Adjust d0 and dK to the new redundancy
        cbw = 2*sum(audfiltbw(fc(2:M2-1),flags.audscale)/winbw*kv.bwmul)/(kv.redtar*fs);
        % Low-pass
        fps0 = audtofreq(freqtoaud(fc(2),flags.audscale)+3*kv.spacing,flags.audscale);% f_{p,s}^{-}
        fsupp_temp0 = audfiltbw(fps0,flags.audscale)/winbw*kv.bwmul;
        a_new(1,2) = ceil(Ls/max(fs./(2*fps0+fsupp_temp0/cbw),1));
        % High-pass
        fps1 = audtofreq(freqtoaud(fc(end-1),flags.audscale)-3*kv.spacing,flags.audscale);% f_{p,s}^{+}
        fsupp_temp1 = audfiltbw(fps1,flags.audscale)/winbw*kv.bwmul;
        a_new(end,2) = ceil(Ls/max(fs./(2*(fc(end)-fps1)+fsupp_temp1/cbw),1));
        % Finally re-scale all filters
        dk_new = a_new(:,1)./a_new(:,2);
        scal_new = dk_new./dk_old;
        % new_red = sum(2./dk_new)-sum(1./dk_new([1,end]));
    end
    g_new = filterbankscale(g,sqrt(scal_new)'); % Perform rescaling
    a = a_new;
    g = g_new;
    if 0
        % Compute and display redundancy for verification
        fprintf('Original redundancy: %g \n', org_red);
        fprintf('Target redundancy: %g \n', kv.redtar);
        fprintf('Actual redundancy: %g \n', new_red);
    end
end

winbwrat = winbw/0.754;
basebw = 1.875657;

info.fc  = 2*fc/fs;
info.tfr = @(L)(1/L)*1./((2*fsupp*winbwrat/fs)./basebw).^2;

function glow = audlowpassfilter(g,a,fc,fs,scal,kv,flags)

% Make a probe, compute the restricted filter bank response to check if low-pass filter is needed
Lprobe = 10000;
FBresp0 = filterbankresponse(g(2:end-1),a(2:end-1,:),Lprobe,'real');
eps_thr = 1e-3;
ind_f1 = floor(fc(2)*Lprobe/fs);
ind_fK = floor(fc(end-1)*Lprobe/fs);
if ind_f1 == 0 || ...
    min(FBresp0(1:ind_f1)) >= (1-eps_thr)*min(FBresp0(ind_f1:ind_fK))
%       Not required
    glow.H = @(L) 0;
    glow.foff = @(L) 0;
    glow.realonly = 0;
    glow.delay = 0;
    glow.fs = g{2}.fs;
else

%       Required
    % Compute the transition frequencies f_{p,s}^{-} and f_{p,e}^{-}
    % Determines the width of the plateau
    fps = audtofreq(freqtoaud(fc(2),flags.audscale)+3*kv.spacing,flags.audscale);
    % Determines the cosine transition frequency
    fpe = audtofreq(freqtoaud(fc(2),flags.audscale)+4*kv.spacing,flags.audscale);
    fsupp_LP = 2*fpe;
    ratio = 2*(fpe-fps)/fsupp_LP;
    Lw = @(L) min(ceil(fsupp_LP*L/fs),L);

    P0 = blfilter({'hann','taper',ratio},fsupp_LP,'fs',fs,'inf','min_win',kv.min_win);
    temp_fbresp = @(L) filterbankresponse(g(2:end-1),a(2:end-1,:),L,'real');
    Hinv = @(L) sqrt(max(temp_fbresp(L))-temp_fbresp(L));

%     Compute the final low-pass filter
    glow.H = @(L) fftshift(long2fir(...
        filterbankfreqz(P0,a(1,:),L).*Hinv(L),Lw(L)))*scal;
    glow.foff = @(L) -floor(Lw(L)/2);
    glow.realonly = 0;
    glow.delay = 0;
    glow.fs = g{2}.fs;
end

function ghigh = audhighpassfilter(g,a,fc,fs,scal,kv,flags)

% Make a probe, compute the restricted filter bank response to check if hi-pass filter is needed
Lprobe = 10000;
FBresp0 = filterbankresponse(g(2:end-1),a(2:end-1,:),Lprobe,'real');
eps_thr = 1e-3;
ind_f1 = floor(fc(2)*Lprobe/fs);
ind_fK = floor(fc(end-1)*Lprobe/fs);
if ind_f1 == 0 ||...
    min(FBresp0(ind_fK:floor(Lprobe/2))) >= (1-eps_thr)*min(FBresp0(ind_f1:ind_fK))
%       Not required
    ghigh.H = @(L) 0;
    ghigh.foff = @(L) 0;
    ghigh.realonly = 0;
    ghigh.delay = 0;
    ghigh.fs = g{2}.fs;
else

    %     Compute the transition frequencies f_{p,s}^{+} and f_{p,e}^{+}
    % Determines the width of the plateau
    fps = audtofreq(freqtoaud(fc(end-1),flags.audscale)-3*kv.spacing,flags.audscale);
    % Determines the cosine transition frequency
    fpe = audtofreq(freqtoaud(fc(end-1),flags.audscale)-4*kv.spacing,flags.audscale);

    %     plateauWidth = 2*(fs/2-fps);
    fsupp_HP = 2*(fs/2-fpe);
    ratio = 2*(fps-fpe)/fsupp_HP;
    Lw = @(L) min(ceil(fsupp_HP*L/fs),L);

    PK = blfilter({'hann','taper',ratio},fsupp_HP,'fc',fs/2,'fs',fs,'inf','min_win',kv.min_win);
    temp_fbresp = @(L) filterbankresponse(g(2:end-1),a(2:end-1,:),L,'real');
    Hinv = @(L) sqrt(max(temp_fbresp(L))-temp_fbresp(L));

    %     Compute the final high-pass filter
    ghigh.H = @(L) fftshift(long2fir(fftshift(...
        filterbankfreqz(PK,a(1,:),L).*Hinv(L)),Lw(L)))*scal;

    ghigh.foff = @(L) ceil(L/2)-floor(Lw(L)/2)-1;
    ghigh.realonly = 0;
    ghigh.delay = 0;
    ghigh.fs = g{2}.fs;
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function test_min_test02 ( )

%*****************************************************************************80
%
%% TEST_MIN_TEST02 evaluates the objective function at each starting point.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_MIN_TEST02\n' );
  fprintf ( 1, '  For each problem, evaluate the function\n' );
  fprintf ( 1, '  at the starting point and the solution.\n' );
%
%  Get the number of problems.
%
  problem_num = p00_problem_num ( );

  for problem = 1 : problem_num

    title = p00_title ( problem );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Problem %d\n', problem );
    fprintf ( 1, '  %s\n', title );
    fprintf ( 1, '\n' );

    x = p00_start ( problem );

    f_start = p00_f ( problem, x );

    fprintf ( 1, '    F(X_START) = %e\n', f_start );

    [ know, x ] = p00_sol ( problem );

    if ( 0 < know )
      f_sol = p00_f ( problem, x );
      fprintf ( 1, '    F(X_SOL) = %e\n', f_sol );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [status,results] = mrtrix_cmd(cmd_str, bkgrnd, verbose)
% function [status,results] = mrtrix_cmd(cmd_str, [bkgrnd=true], [verbose=true])
%
% Send a command to an mrtrix function and get back status and results
%
% INPUTS:
%   cmd_str - A string containing the mrtrix command to run
%   bkgrnd  - [true/false] Whether to execute the command in the background (only possible on unix)
%   verbose - Whether to display stdout to the command window (when it's done).
%
% OUTPUTS:
%   status  - Whether the operation succeeded (1) or not (0)
%   results - The results of the operation in stdout.
%
% Notes:
%  When bkgrnd is set to true, the command will be executed in another
%  terminal.
%
%
% Franco Pestilli, Ariel Rokem & Bod Dougherty Stanford Univesity

if notDefined('bkgrnd'), bkgrnd = false;end
if notDefined('verbose'),verbose = true;end

% Need to bypass the matlab libraries at the top of the path, by screwing
% with it:
orig_ld_path = mrtrix_set_ld_path;

% This opens another xterm, runs it there and gives you back control of the
% matlab session:
if bkgrnd
    if ~isunix
        bkgrnd = false;
        warning('Cannot run in background on MS Windows. Defaulting to run mrtrix command in the Matlab session');
    else
        cmd_str = ['xterm -e ' cmd_str ' &'];
    end
end
fprintf('\n[%s] Running the following command: \n%s\n',mfilename,cmd_str);

% Run the command and get back status and results:
[status,results] = system(cmd_str);

% If there was a failure, throw a warning:
if (status ~=0)
    warning('[%s] There was a failure running the command. \nCommand: %s\nError %s\n',mfilename,results);
end

% If running the command within the matlab session and verbose, display
% results, when it returns:
if (~bkgrnd && verbose), disp(results); end

% Reset the ld_path:
setenv('LD_LIBRARY_PATH', orig_ld_path);

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [fg, fgRemains, fgLeft, fgRight]=dtiSplitInterhemisphericFibers(fg, dt, maxZ)
%Split fibers that cross between hemispheres below a certain z-coordinate.
%
%   fg=dtiSplitInterhemisphericFibers(fg, dt, [Z=-10])
%
% Due to tractography artifacts, whole brain tractography often produces
% fibers which connect homologous motor cortices, crossing between
% hemispheres at the pons level. This is anatomically implausible: these
% fibers mostl likely cross at the pons level, connecting spinal tract and
% contralateral moter cortex. This function splits interhemisperic fibers,
% sparing (dep. on option chosen) either CC, or ACPC & above. Note:
% cerebellar fibers WILL be split as well. Too bad.
%
% Input parameters:
% fg      - fiber group structure
% dt      - dt6 (tensor) data structure
% maxZ    - (a) maxZ=-10 (default) will split fibers, cutting them by a
%           saggital plane 10 mm below acpc line.  We assume that no
%           interhemispheric fiber should connect below ACPC line. Fornix &
%           CC are comissures located above. Provide your own maxZ (the
%           most superior z coordinate of the cutting midsaggital plane)
%           otherwise. (b) maxZ='AllButCC' allows splitting every
%           interhemispheric fiber except callosal.
%
% Example:
% Split intehemispheric fibers below ACPC
%   dt=dtiLoadDt6(dtFileName);
%   fg = dtiLoadFibers(fgFileName);
%   fg=dtiSplitInterhemisphericFibers(fg, dt, 0]);
%
% (c) Vistalab

%HISTORY:
% 08/2009: ER & LMP   wrote it
% 08/24/2009 ER: drop fibers of 1 node long.

if notDefined('maxZ')
    maxZ=-10;
end

if isnumeric(maxZ) && (maxZ<dt.bb(1, 3) ||maxZ>dt.bb(2, 3))|| (~isnumeric(maxZ) && ~strcmp(maxZ, 'AllButCC'))
    error(['The most superior point Z for cutting should be between ' num2str(dt.bb(1, 3)) ' and ' num2str(dt.bb(2, 3)) ' or "AllButCC"']);
end

if notDefined('dt')
    [f,p] = uigetfile({'*.nii.gz';'*.*'},'Select dt6 file...');
    if(isnumeric(f)), disp('User canceled.'); return; end
    dt = fullfile(p,f);
end

midSagThresh = 5;
fg=dtiCleanFibers(fg);
fgname=fg.name;

if isnumeric(maxZ)

   fprintf(1, 'dtiSplitInterhemisphericFibers: Splitting every fiber below Z=%s \n', num2str(maxZ));
   % Bounding Box = dt.bb  x = [-80 0]; y = [-120 90]; z = [-60 90];
    x = [0];
    y = [dt.bb(1, 2):dt.bb(2, 2)];
    z = [dt.bb(1, 3): maxZ];

    [X,Y,Z] = meshgrid(x, y, z);
    roiCoords = [X(:), Y(:), Z(:)];
    roi = dtiNewRoi('MidSaggitalBelowACPC', 'b', roiCoords);

    [fgToChop,contentiousFibers, keep] = dtiIntersectFibersWithRoi([], {'and'}, [], roi, fg);
    fgRemains = dtiNewFiberGroup([fg.name 'Remains']); fgRemains.fibers=fg.fibers(~keep);
    if isfield(fg, 'subgroup') && ~isempty(fg.subgroup)
        fgRemains.subgroup=fg.subgroup(~keep);
    end

    fgRight = dtiNewFiberGroup([fg.name '_ChoppedR']);
    fgLeft = dtiNewFiberGroup([fg.name '_ChoppedL']);
        clear fg;

%emptyleft=0;
%emptyright=0;

    keepRightID=[];keepLeftID=[];
    for i=1:size(fgToChop.fibers)
        pointsRight=(fgToChop.fibers{i}(1, :)>midSagThresh);
        pointsLeft=(fgToChop.fibers{i}(1, :)<-midSagThresh);
        RightChunk{1}=fgToChop.fibers{i}(:, pointsRight);
        if ~isempty(RightChunk{1})
            fgRight.fibers= [fgRight.fibers(:);RightChunk] ; %Right chunk
            keepRightID=[keepRightID; i];
 %       else emptyright=emptyright+1;
        end
        LeftChunk{1}=fgToChop.fibers{i}(:,  pointsLeft);

        if ~isempty(LeftChunk{1})
            fgLeft.fibers = [fgLeft.fibers(:); LeftChunk]; %Left chunk
            keepLeftID=[keepLeftID; i];
 %       else emptyleft=emptyleft+1;
        end
    end
      if isfield(fgToChop, 'subgroupNames') && ~isempty(fgToChop.subgroupNames)
    fgLeft.subgroupNames=fgToChop.subgroupNames;
    fgRight.subgroupNames=fgToChop.subgroupNames;
      end
    if isfield(fgToChop, 'subgroup') && ~isempty(fgToChop.subgroup)
    fgLeft.subgroup=fgToChop.subgroup(keepLeftID);
    fgRight.subgroup=fgToChop.subgroup(keepRightID);
    end

%    emptyleft
%    emptyright
else

    % Code that LMP wrote to take a fiber group, remove callosal fibers, split
    % the remaining groups, then merge the two groups to create one set that
    % is cut down the mid-line but retains the callosal fibers.
    %Limitation: assumes CC is the only interhemispheric connection
    fprintf(1, 'Splitting every fiber except in the CC \n');

    ccCoords = dtiFindCallosum(dt.dt6,dt.b0,dt.xformToAcpc);
    ccRoi = dtiNewRoi('CC','c',ccCoords);
    ccRoi = dtiRoiClean(ccRoi, 3, {'dilate'});

    fgRemains = dtiIntersectFibersWithRoi([], {'and'}, [], ccRoi, fg);
    fgToChop= dtiIntersectFibersWithRoi([], {'not'}, [], ccRoi, fg);
    clear fg;
    fgRight = dtiClipFiberGroup(fgToChop, [-80 midSagThresh],[],[]);
    fgLeft = dtiClipFiberGroup(fgToChop, [-midSagThresh 80],[],[]);

end

fg = dtiMergeFiberGroups(fgRight,fgLeft,fgname);
fg = dtiMergeFiberGroups(fg,fgRemains,fgname);
fgRemains.name  = [fgname '_NoChop'];
fgLeft.name  = [fgname '_L'];
fgRight.name  = [fgname '_R'];
fg.name=[fgname '_BilaterallySplit'];

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%dtiWriteFibersPdb(newFg,[],newFg.name);

%fgToChop.name=['allConnectingGM_withCST_Mori_ToChop'];
%dtiWriteFiberGroup(fgRemains,fgRemains.name);
%dtiWriteFiberGroup(newFgLeft,newFgLeft.name);
%dtiWriteFiberGroup(newFgRight,newFgRight.name);
%dtiWriteFiberGroup(fgToChop,fgToChop.name);


text:
function opto = snirf2opto(probe, measurementList)

% SNIRF2OPTO converts the SNIRF probe and measurementList structures to a FieldTrip
% optode structure.
%
% See https://github.com/fNIRS/snirf/blob/master/snirf_specification.md
%
% The FieldTrip optode structure is defined in FT_DATATYPE_SENS
%
% See also OPTO2HOMER, BTI2GRAD, CTF2GRAD, FIF2GRAD, ITAB2GRAD, MNE2GRAD, NETMEG2GRAD, YOKOGAWA2GRAD, FT_DATATYPE_SENS

% Copyright (C) 2020, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% 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$

nSrc = length(probe.sourceLabels);
nDet = length(probe.detectorLabels);

opto = [];
opto.optolabel = cat(1, probe.sourceLabels(:), probe.detectorLabels(:)); % ensure that they are columns

if ~isempty(probe.sourcePos3D)
  % use the 3D positions if available
  opto.optopos = cat(1, probe.sourcePos3D, probe.detectorPos3D);
else
  % use the 2D positions
  opto.optopos = cat(1, probe.sourcePos2D, probe.detectorPos2D);
  opto.optopos(:,3) = 0; % convert the positions to 3D
end

opto.optotype = cat(1, repmat({'transmitter'}, nSrc, 1), repmat({'receiver'}, nDet, 1));
opto.wavelength = probe.wavelengths(:)';

M = length(measurementList);    % number of channels
N = nSrc+nDet;                  % number of optodes
K = numel(probe.wavelengths);   % number of wavelengths

opto.label = cell(M,1);
for i=1:M
  tx = measurementList(i).sourceIndex;                        % transmitter
  rx = measurementList(i).detectorIndex;                      % receiver
  wl = probe.wavelengths(measurementList(i).wavelengthIndex); % wavelength in nm
  opto.label{i} = sprintf('%s-%s [%dnm]', probe.sourceLabels{tx}, probe.detectorLabels{rx}, round(wl));
end

% the following specifies for each of the M channels at which wavelength each of the
% N optodes transmits (positive integer from 1 to K), or receives (negative ingeger
% from 1 to K), or does not contribute at all (zeros)

opto.tra = zeros(M, N);
for i=1:M
  tx = measurementList(i).sourceIndex;
  rx = measurementList(i).detectorIndex + nSrc; % the transmitters are first in the list of optodes
  opto.tra(i, tx) = +measurementList(i).wavelengthIndex;
  opto.tra(i, rx) = -measurementList(i).wavelengthIndex;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function matrix_output=padded_cell2mat(cell_input,pad)
    % Convert a ragged cell array of arrays in a padded matrix
    % Pierre Morel 2015.

    %With cellfun
%     ticID=tic;
%     maxLength=max(cellfun(@numel,cell_input));
%     tempcell=cellfun(@(x)horzcat(2,x,NaN*zeros(1,maxLength-length(x))),cell_input,'UniformOutput',false);
%     matrix_output=cell2mat(tempcell);
%     disp(['padded matrix (cellfun) created in ' num2str(toc(ticID)) ' s'])

    if nargin<2
        pad=NaN;
    end

    %With for loop (faster than cellfun)
    %ticID=tic;
    lengths=cellfun(@numel,cell_input);
    maxLength=max(lengths);
    matrix_output=zeros(length(cell_input),maxLength)+pad;
    for k=1:length(cell_input)
        matrix_output(k,1:lengths(k))=cell_input{k};
    end
    %disp(['padded matrix (for loop)  created in ' num2str(toc(ticID)) ' s'])
end

author: piermorel

repo: gramm

save_path: github-repos/MATLAB/piermorel-gramm


text:
function allBb2d = crop2DBB(allBb2d,imh,imw)
        allBb2d_comp1 = allBb2d(:,[1,2]);
        allBb2d_comp1(allBb2d_comp1<1) =1;
        allBb2d_comp3 = allBb2d(:,1)+allBb2d(:,3);
        allBb2d_comp3(allBb2d_comp3>imw)=imw;
        allBb2d_comp4 = allBb2d(:,2)+allBb2d(:,4);
        allBb2d_comp4(allBb2d_comp4>imh)=imh;
        if size(allBb2d,2)>4
            allBb2d = [allBb2d_comp1,allBb2d_comp3,allBb2d_comp4,allBb2d(:,5:end)];
        else
            allBb2d = [allBb2d_comp1,allBb2d_comp3,allBb2d_comp4];
        end
        allBb2d(:,[3,4])= allBb2d(:,[3,4])-allBb2d(:,[1,2]);
        allBb2d(allBb2d(:,3)<0| allBb2d(:,4)<0,:)=[];
end

author: thusiyuan


text:
classdef Rx < adi.AD9361.Rx

    methods
        %% Constructor
        function obj = Rx(varargin)
            % Returns the matlabshared.libiio.base object
            coder.allowpcode('plain');
            obj = obj@adi.AD9361.Rx(varargin{:});
        end
    end

end

author: analogdevicesinc


text:
function [celBodyData] = processINIBodyInfo(celBodyDataFromINI, varargin)
%processINIBodyInfo Summary of this function goes here
%   Detailed explanation goes here
    clear findRow

    if(length(varargin)>=1)
        showWaitbar = varargin{1};

        if(length(varargin)>=2)
            dataType = varargin{2};

            if(length(varargin) >= 3)
                hFig = varargin{3};
            else
                hFig = [];
            end
        else
            dataType = 'bodyInfo';
            hFig = [];
        end
    else
        showWaitbar = true;
        dataType = 'bodyInfo';
        hFig = [];
    end

    if(showWaitbar)
        if(not(isempty(hFig)))
            hWaitbar = uiprogressdlg(hFig,'Message','Loading bodies file, please wait...', 'Title','KSP Trajectory Optimizatoon Tool');
        else
            hWaitbar = waitbar(0,'Loading bodies file, please wait...');
        end
    end

    celBodyData = struct();
    for(i=1:size(celBodyDataFromINI,1)) %#ok<*NO4LP>
        row = celBodyDataFromINI(i,:);
        row{1} = lower(matlab.lang.makeValidName(row{1}));
        row{3} = lower(matlab.lang.makeValidName(row{3}));

        if(~isfield(celBodyData,row{1}))
%             celBodyData.(row{1}) = struct();
            celBodyData.(row{1}) = getNewObjForDataType(dataType);
        end

        if(strcmpi(row{3}, 'atmoPressAlts'))
            celBodyData = createCurveFitFromRows(celBodyDataFromINI, row{1}, 'atmoPressAlts', 'atmoPressPresses', 'atmoPressCurve', celBodyData);

        elseif(strcmpi(row{3}, 'atmoTempAlts'))
            celBodyData = createCurveFitFromRows(celBodyDataFromINI, row{1}, 'atmoTempAlts', 'atmoTempTemps', 'atmoTempCurve', celBodyData);

        elseif(strcmpi(row{3}, 'atmoTempSunMultAlts'))
            celBodyData = createCurveFitFromRows(celBodyDataFromINI, row{1}, 'atmoTempSunMultAlts', 'atmoTempSunMults', 'atmoTempSunMultCurve', celBodyData);

        elseif(strcmpi(row{3}, 'latTempBiasLats'))
            celBodyData = createCurveFitFromRows(celBodyDataFromINI, row{1}, 'latTempBiasLats', 'latTempBiases', 'latTempBiasCurve', celBodyData);

        elseif(strcmpi(row{3}, 'latTempSunMultLats'))
            celBodyData = createCurveFitFromRows(celBodyDataFromINI, row{1}, 'latTempSunMultLats', 'latTempSunMults', 'latTempSunMultCurve', celBodyData);

        elseif(strcmpi(row{3}, 'axialtemperaturesunbiasdeg'))
           celBodyData = createCurveFitFromRows(celBodyDataFromINI, row{1}, 'axialtemperaturesunbiasdeg', 'axialtemperaturesunbiases', 'axialTempSunBiasCurve', celBodyData);

        elseif(strcmpi(row{3}, 'axialtemperaturesunmultlats'))
            celBodyData = createCurveFitFromRows(celBodyDataFromINI, row{1}, 'axialtemperaturesunmultlats', 'axialtemperaturesunmults', 'axialTempSunMultCurve', celBodyData);

        elseif(strcmpi(row{3}, 'eccentricitytemperaturebiaspts'))
            celBodyData = createCurveFitFromRows(celBodyDataFromINI, row{1}, 'eccentricitytemperaturebiaspts', 'eccentricitytemperaturebiases', 'eccTempBiasCurve', celBodyData);

        elseif(any(strcmpi(row{3}, {'atmoPressPresses', 'atmoTempTemps', 'atmoTempSunMults', 'latTempBiases', 'latTempSunMults', 'axialtemperaturesunbiases', 'axialtemperaturesunmults', 'eccentricitytemperaturebiases'})))
            continue;
        else
            if(~checkStrIsNumeric(row{4}))
                entry = row{4};
            else
                entry = str2double(row{4});
            end

            if((strcmpi(row{3},'name') || strcmpi(row{3},'parent')) && ischar(entry))
                entry = matlab.lang.makeValidName(entry);
            end

            if(isfield(celBodyData,row{1}))
                celBodyData.(row{1}).(row{3}) = entry;
            else
                celBodyData.(row{1}) = struct();
                celBodyData.(row{1}) = getNewObjForDataType(dataType);
                celBodyData.(row{1}).(row{3}) = entry;
            end
        end

        if(showWaitbar)
            waitbarValue = i/size(celBodyDataFromINI,1);
            waitbarMsg = sprintf('Loading bodies file, please wait... [%u/%u]', i, size(celBodyDataFromINI,1));

            if(not(isempty(hFig)))
                hWaitbar.Value = waitbarValue;
                hWaitbar.Message = waitbarMsg;
            else
                waitbar(waitbarValue, hWaitbar, waitbarMsg);
            end
        end
    end

    if(strcmpi(dataType,'bodyInfo'))
        names = fieldnames(celBodyData);
        for(i=1:length(names))
            name = names{i};
            celBodyData.(name).celBodyData = celBodyData;
            [~] = celBodyData.(name).getParBodyInfo(celBodyData); %set that parent info now so that we don't have to handle it later
            [~] = celBodyData.(name).getChildrenBodyInfo(celBodyData);

            celBodyData.(name).setPropTypeEnum();
            celBodyData.(name).setBaseBodySurfaceTexture();
        end
    end

    if(showWaitbar)
        close(hWaitbar);
    end

    return;
end

function celBodyData = createCurveFitFromRows(celBodyDataFromINI, bodyName, indVarRow, depVarRow, celBodyFieldName, celBodyData)
    indRow = findRow(celBodyDataFromINI, bodyName, indVarRow);
    depRow = findRow(celBodyDataFromINI, bodyName, depVarRow);

    indVar = str2double(strsplit(indRow{4},','));
    indVar = reshape(indVar, length(indVar), 1);

    depVar = str2double(strsplit(depRow{4},','));
    depVar = reshape(depVar, length(depVar), 1);

    if(length(indVar) > 1 && length(depVar) > 1)
        fitobject = griddedInterpolant(indVar,depVar,'linear', 'nearest');
%         fitobject = fit(indVar,depVar,'smoothingspline');
%         fitobject = @(xi) interp1qr(indVar,depVar,xi);
    elseif(length(indVar) == 1 && length(depVar) == 1)
        fitobject = griddedInterpolant([indVar(1) indVar(1)+1]',[depVar(1) depVar(1)]','linear', 'nearest');
%         fitobject = fit([indVar(1) indVar(1)+1]',[depVar(1) depVar(1)]','smoothingspline');
%         fitobject = @(xi) interp1qr([indVar(1) indVar(1)+1]',[depVar(1) depVar(1)]',xi);
    else
        fitobject = griddedInterpolant([0 1]',[0 0]','linear', 'nearest');
%         fitobject = fit([0 1]',[0 0]','smoothingspline');
%         fitobject = @(xi) interp1qr([0 1]',[0 0]',xi);
    end

    celBodyData.(bodyName).(lower(celBodyFieldName)) = fitobject;
end

function row = findRow(celBodyDataFromINI, col1, col3)
    persistent celBodyDataFromINICol1
    if(isempty(celBodyDataFromINICol1) || size(celBodyDataFromINI,1) ~= length(celBodyDataFromINICol1))
        celBodyDataFromINICol1 = matlab.lang.makeValidName(celBodyDataFromINI(:,1));
    end

    row = celBodyDataFromINI(strcmpi(col1,celBodyDataFromINICol1) & strcmpi(col3,celBodyDataFromINI(:,3)), :);
end

function obj = getNewObjForDataType(dataType)
    switch dataType
        case 'bodyInfo'
            obj = KSPTOT_BodyInfo();
        case 'appOptions'
            obj = KSPTOT_AppOptions();
        otherwise
            error('Unknown data type %s in processINIBodyInfo().', dataType);
    end
end


text:
%{
B_LoopDeLoopCam.m

Animate several camera motions around a surface.

Matt Sheen, mws262@cornell.edu
%}

close all; clear all;
fig = figure;

%% Compound Patch - same from folder 2, example b
pX = [-1 1 0;
    0 0 0]';
pY = [-1/3 -1/3 2/3;
    -1/3 -1/3 2/3]';
pZ = [0 0 0;
    0 0.5 0]';

p1 = patch(pX,pY,pZ,'red');
p1.FaceAlpha = 0.5;

OrigVerts = p1.Vertices;

%% Set initial camera position.
axis(2*[-10 10 -10 10 -10 10])
campos([8,8,8]);
camtarget([0,0,0]);
camva(40);

%Turning the axis grid on to help get perspective of movement
fig.Children.ZGrid = 'on';
fig.Children.YGrid = 'on';
fig.Children.XGrid = 'on';

vel = 0.1*[0 1 0];
rot = angle2dcm(0,0,-0.02);
dt = 0.5;
center = [0 0 0];

totalRot = eye(3,3);

%% Cam mode 1: fixed everything. Camera stays exactly the same throughout.
input('Press enter to see default camera orientation.');
for i = 1:1000
    %Moving the plane is the same as previous example
    totalRot = totalRot*rot;
    center = dt*(totalRot*vel')' + center; %transform the velocity, keep track of the center position.
    p1.Vertices = (totalRot*OrigVerts')' + repmat(center,[size(OrigVerts,1),1]); %Rotate the patch object about its center and add new center point.

    pause(0.01);
end

%% Cam mode 2: let's follow the plane. Camera keeps the plane in the center of view, but does not change orientation.
input('Press enter to see plane-centering camera.');
for i = 1:1000
    %Moving the plane is the same as previous example
    totalRot = totalRot*rot;
    center = dt*(totalRot*vel')' + center; %transform the velocity, keep track of the center position.
    p1.Vertices = (totalRot*OrigVerts')' + repmat(center,[size(OrigVerts,1),1]); %Rotate the patch object about its center and add new center point.

    %CAMERA CHANGES:
    camtarget(center); %make the target be the center of the plane.
    campos(center+[5,5,5]); %make the camera position be a fixed offset from the center of the plane.
    pause(0.01);
end

%% Cam mode 3: let's fly the plane. Camera also changes orientation with the plane.
input('Press enter to see plane-following camera.');
for i = 1:1000
    %Moving the plane is the same as previous example
    totalRot = totalRot*rot;
    center = dt*(totalRot*vel')' + center; %transform the velocity, keep track of the center position.
    p1.Vertices = (totalRot*OrigVerts')' + repmat(center,[size(OrigVerts,1),1]); %Rotate the patch object about its center and add new center point.

    %CAMERA CHANGES:
    camtarget(center); %Camera target is the center of the plane.
    campos(center+5*(totalRot*[0.5 -1 0.5]')'); %Camera position is a little behind, up, and to the right of the plane (helps see what's going on). Note that we have to transform the camera offset from plane frame to world frame.
    camup(totalRot*[0 0 1]'); % The camera should think "up" is the up vector of the plane. So we take the plane's up and transform it into the world frame.

    pause(0.01);
end

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
classdef CPD
    % StaticCPD
    %  Defines the conditional probability distributions that make up the
    %  BPL model

    methods (Static)

        % --
        % Model of global statistics
        % --
        function [logwt,x_ink,x_cov] = score_stat(libclass,M)
            % Compute the log-score of a motor program M, in terms
            % of the simple ink and coverage statistics.
            x_ink = get_ink_overlap(M,'stroke');
            x_cov = get_ratio_coverage_type(M);
            log_ink = libclass.stat.q_ink_emp.pdf(x_ink) - libclass.stat.q_ink_prior.pdf(x_ink);
            log_canvas = libclass.stat.q_canvas_emp.pdf(x_cov) - libclass.stat.q_canvas_prior.pdf(x_cov);
            logwt = log_ink + log_canvas;
        end

        % ---
        % Number of strokes model (Kappa)
        % ---

        function samps = sample_number(libclass,nsamp)
            if ~exist('nsamp','var')
               nsamp = 1;
            end
            samps = zeros(nsamp,1);
            for i = 1:nsamp
                samps(i) = find(mnrnd(1,libclass.pkappa));
            end
        end

        function ll = score_number(libclass,data)
        %   data: vector [n x 1] of number of strokes
        % returns
        %   ll : vector [n x 1] of log-likelihood scores
            assert(isvector(data));
            if any(data>length(libclass.pkappa))
               ll = -inf;
               return
            end
            ll = log(libclass.pkappa(data));
        end

        % ---
        % Sequence of sub-strokes model (z)
        % ---

        function nsub = sample_nsub(libclass,ns)
        % ns    : number of strokes in the character
            pvec = libclass.pmat_nsub(ns,:);
            nsub = find(mnrnd(1,pvec));
            assert(isscalar(nsub));
            assert(~isnan(nsub));
        end

        function samps = sample_sequence(libclass,ns,nsub,nsamp)
        % ns    : number of strokes in the character
        % nsub  : number of sub-strokes in the character (optional)
        % nsamp : number of samples
            if ~exist('nsamp','var')
               nsamp = 1;
            end
            samps = cell(nsamp,1);

            % Produce each of the samples
            for i=1:nsamp

                %% number of sub-strokes
                if (~exist('nsub','var') || isempty(nsub))
                    nsub = CPD.sample_nsub(libclass,ns);
                end

                %% sample the sequence
                sq = zeros(nsub,1);
                pStart = exp(libclass.logStart);
                sq(1) = find(mnrnd(1,pStart));
                for bid=2:nsub
                    prev = sq(bid-1);
                    sq(bid) = find(mnrnd(1,libclass.pT(prev)));
                end
                samps{i} = sq;
            end

            if nsamp==1
               samps = samps{1};
            end

        end

        function [ll,lv] = score_sequence(libclass,ns,ids)
        %  ids: [k x 1] data sequence
        %  ns: [scalar] how many strokes are in the character. If ns<0,
        %    this means to leave out scoring the number of sub-strokes
        %  output is a scalar log-likelihood value
            assert(isvector(ids));
            nsub = length(ids);
            nstate = libclass.N;
            assert(all(ids<=nstate));

            %% score the number of sub-strokes
            sz = size(libclass.pmat_nsub);
            if ns < 0 % do not include this part of the score...
               lp = 0;
            elseif ns > sz(1) || nsub > sz(2)
               lp = -inf;
            else
               lp = log(libclass.pmat_nsub(ns,nsub));
            end

            %% score the sub-stroke ids
            lv = zeros(nsub,1);
            lv(1) = libclass.logStart(ids(1));
            for i=2:nsub
               lv(i) = libclass.logT(ids(i-1),ids(i));
            end
            ll = lp + sum(lv);
            assert(~any(isnan(ll)));

        end

        % ---
        % Shape model (x)
        % ---
        function bspline_stack = sample_shape_type(libclass,subid)
        % Given a vector of cluster ids (id), sample
        % from the vanilla gaussian model
        % associated with it
        %
        % Input
        %  subid [k x 1]: the id (index) of the sub-stroke
        %
        % Output
        %  bspline_stack: [ncpt x 2 x k] sampled spline
             assert(isvector(subid));
             k = length(subid);

             if isunif(libclass)
               bspline_stack = CPDUnif.sample_shape_type(libclass,subid);
               return
             end

             Cov = libclass.shape.Sigma(:,:,subid);
             rows_bspline = mvnrnd(libclass.shape.mu(subid,:),Cov);
             ncpt = size(rows_bspline,2)./2;
             bspline_stack = zeros(ncpt,2,k);
             for i=1:k
                 bspline_stack(:,:,i) = reshape(rows_bspline(i,:),[ncpt 2]);
             end
        end

        function ll = score_shape_type(libclass,bspline_stack,subid)
        %
        % Input
        %  subid [k x 1]: the id (index) of the sub-stroke
        %  bspline_stack: (ncpt x 2 x k) shapes of bsplines
        %
        % Output
        %  ll : [k x 1]: vector of scores
            assert(isvector(subid));
            k = length(subid);
            assert(size(bspline_stack,3)==k);

            if isunif(libclass)
               ll = CPDUnif.score_shape_type(libclass,bspline_stack,subid);
               return
            end

            ll = score_shape_type_helper(libclass,bspline_stack,subid);
        end

        function bspline_stack = sample_shape_token(libclass,bspline_stack)
        %  bspline_stack: (ncpt x 2 x k) shapes of bsplines
             sz = size(bspline_stack);
             bspline_stack = bspline_stack + libclass.tokenvar.sigma_shape .* randn(sz);
        end

        function ll = score_shape_token(libclass,bspline_stack_token,bspline_stack_type)
        %
        % Input
        %  bspline_stack_type: (ncpt x 2 x k) shapes of bsplines
        %  bspline_stack_token: (ncpt x 2 x k) shapes of bsplines
        %
        % Output
        %  ll : [k x 1]: vector of scores
            k = size(bspline_stack_token,3);
            assert(k==size(bspline_stack_type,3));
            sz = size(bspline_stack_token);
            X = bspline_stack_token(:)';
            mu = bspline_stack_type(:)';
            llvec = mvnormpdfln(X,mu,libclass.tokenvar.sigma_shape*ones(size(X)));
            llvec = reshape(llvec,sz);
            ll = sum(llvec,1);
            ll = sum(ll,2);
            ll = squeeze(ll);
        end

        function ll = score_shape_marginalize(libclass,bspline_stack_token,subid)
        %
        % Input
        %  subid [k x 1]: the id (index) of the sub-stroke
        %  bspline_stack: (ncpt x 2 x k) shapes of bsplines
        %
        % Output
        %  ll : [k x 1]: vector of scores
            assert(isvector(subid));
            k = length(subid);
            assert(size(bspline_stack_token,3)==k);

            if isunif(libclass)
               ll = CPDUnif.score_shape_marginalize(libclass,bspline_stack_token,subid);
               return
            end

            ll = score_shape_type_helper(libclass,bspline_stack_token,subid,libclass.tokenvar.sigma_shape^2);
        end

        % --
        % Scale model (y)
        % --

        function invscales = sample_invscale_type(libclass,subid)
            % subid: [k x 1] vector of sub-stroke ids
            %
            assert(isvector(subid));

            if isunif(libclass)
               invscales = CPDUnif.sample_invscale_type(libclass,subid);
               return
            end

            theta = libclass.scale.theta(subid,:);
            invscales = gamrnd(theta(:,1),theta(:,2));
        end

        function lprob = score_invscale_type(libclass,invscales_type,subid)
            % subid: [k x 1] vector of sub-stroke ids
            assert(isvector(invscales_type));
            assert(isvector(subid));
            k = length(subid);
            assert(numel(invscales_type)==k);

            if isunif(libclass)
                lprob = CPDUnif.score_invscale_type(libclass,invscales_type,subid);
                return
            end

            theta = libclass.scale.theta(subid,:);
            prob = gampdf(invscales_type(:),theta(:,1),theta(:,2));
            lprob = log(prob);
            assert(numel(lprob)==k);
        end

        function invscales_token = sample_invscale_token(libclass,invscales_type)
        % Gaussian noise, but don't allow negative scales.
        % Sampling is done by rejection sampling
            sz = size(invscales_type);
            invscales_token = invscales_type + libclass.tokenvar.sigma_invscale.*randn(sz);
            while any(isinf(CPD.score_invscale_token(libclass,invscales_token,invscales_type)))
                invscales_token = invscales_type + libclass.tokenvar.sigma_invscale.*randn(sz);
            end
        end

        function ll = score_invscale_token(libclass,invscales_token,invscales_type)
            % Score the token-level inverse scales.
            % Computes the right normalization constant,
            % since negative values are not allowed.
            %
            % Input
            %  invscales_token: [n x 1] vector
            %  invscales_type: [n x 1] vector
            %
            % Output
            %  ll: [n x 1] vector of log-likelihood
            ll = mvnormpdfln(invscales_token(:)',invscales_type(:)',libclass.tokenvar.sigma_invscale);
            ll = ll(:);
            ll(invscales_token <= 0) = -inf; % don't allow invscales that are negative

            % correction for positive only invscales
            p_below = normcdf(0,invscales_type,libclass.tokenvar.sigma_invscale);
            p_above = 1-p_below;
            ll = ll - log(p_above);
        end

        % --
        % Relation model (R)
        % ---
        function R = sample_relation_type(libclass,previous_strokes)
            nprev = length(previous_strokes);
            stroke_num = nprev+1;
            types = {'unihist','start','end','mid'};
            ncpt = libclass.ncpt;
            if nprev == 0
                indx = 1;
            else
                out = mnrnd(1,libclass.rel.mixprob);
                if any(isnan(out))
                   error('improper relation mixing prob.');
                end
                indx = find(out);
            end
            type = types{indx};
            switch (type)
                case 'unihist'
                    gpos = libclass.Spatial.sample(stroke_num);
                    R = RelationIndependent(type,nprev,gpos);
                case {'start','end'}
                    attach_spot = randint(1,1,[1 nprev]);
                    R = RelationAttach(type,nprev,attach_spot);
                case 'mid'
                    attach_spot = randint(1,1,[1 nprev]);
                    nsub = previous_strokes{attach_spot}.nsub;
                    subid_spot = randint(1,1,[1 nsub]);
                    R = RelationAttachAlong(type,nprev,attach_spot,nsub,subid_spot,ncpt);
                    [~,lb,ub] = bspline_gen_s(ncpt,1);
                    R.eval_spot_type = lb + rand*(ub-lb);
                otherwise
                    error('invalid relation');
            end
        end

        function ll = score_relation_type(libclass,R)
            stroke_num = R.nprev+1;
            types = {'unihist','start','end','mid'};
            logp = log(libclass.rel.mixprob);
            if stroke_num > 1 % prior probability on a particular relation
                indx = strcmp(R.type,types);
                ll = logp(indx);
            else % there was only one option
                ll = 0;
            end
            switch (R.type)
                case 'unihist'
                    ll = ll + libclass.Spatial.score(R.gpos,stroke_num);
                case {'start','end'}
                    ll = ll - log(R.nprev);
                case 'mid'
                    ll = ll - log(R.nprev);
                    ll = ll - log(R.nsub);
                    eval_missing = isempty(R.eval_spot_type);
                    if ~eval_missing
                        ncpt = libclass.ncpt;
                        [~,lb,ub] = bspline_gen_s(ncpt,1);
                        ll = ll - log(ub-lb);
                        if (R.eval_spot_type < lb || R.eval_spot_type > ub)
                            ll = -inf;
                        end
                    end
                otherwise
                    error('invalid relation');
            end
        end

        function eval_spot_token = sample_relation_token(libclass,eval_spot_type)
        % sample an attachment, but within the bounds as defined by the model
            eval_spot_token = eval_spot_type + libclass.tokenvar.sigma_attach.*randn;
            while isinf(CPD.score_relation_token(libclass,eval_spot_token,eval_spot_type))
                eval_spot_token = eval_spot_type + libclass.tokenvar.sigma_attach.*randn;
            end
        end

        function ll = score_relation_token_approx_marginalize(libclass,eval_spot_token)
        % score the token evaluation spot, for attached relations, assuming
        % we are missing the type-level spot.
            assert(isscalar(eval_spot_token));
            ncpt = libclass.ncpt;
            [~,lb,ub] = bspline_gen_s(ncpt,1);
            if (eval_spot_token < lb || eval_spot_token > ub)
                ll = -inf;
                return;
            end

            % cached structure that evaluates this complex integral
            ll = libclass.score_eval_marg(eval_spot_token);
            % ll = -log(ub-lb);
        end

        function ll = score_relation_token(libclass,eval_spot_token,eval_spot_type)
        % score the token evaluation spot, for attached relations, assuming
        % the type-level spot is given
            assert(isscalar(eval_spot_token));
            assert(~isempty(eval_spot_type));
            ncpt = libclass.ncpt;
            [~,lb,ub] = bspline_gen_s(ncpt,1);
            if (eval_spot_token < lb || eval_spot_token > ub)
                ll = -inf;
                return
            end
            ll = mvnormpdfln(eval_spot_token,eval_spot_type,libclass.tokenvar.sigma_attach);

            % correction for bounds
            p_within = normcdf(ub,eval_spot_type,libclass.tokenvar.sigma_attach) - normcdf(lb,eval_spot_type,libclass.tokenvar.sigma_attach);
            ll = ll - log(p_within);
        end

        % ---
        % Local position model (L)
        % ---
        function pos = sample_position(libclass,R,previous_strokes)
            % Given a relation R and the previous strokes,
            %  sample where the position of this stroke should be
            base = getAttachPoint(R,previous_strokes);
            pos = [normrnd(base(1),libclass.rel.sigma_x) normrnd(base(2),libclass.rel.sigma_y)];
        end

        function ll = score_position(libclass,pos,R,previous_strokes)
            base = getAttachPoint(R,previous_strokes);
            ll = mvnormpdfln(pos(1),base(1),libclass.rel.sigma_x) + mvnormpdfln(pos(2),base(2),libclass.rel.sigma_y);
        end

        % ---
        % Image blur model (\sigma_b)
        % ---
        function blur_sigma = sample_image_blur(PM)
            % sample uniformly in [min_blur_sigma, max_blur_sigma]
            rg = PM.max_blur_sigma - PM.min_blur_sigma;
            blur_sigma = rg*rand + PM.min_blur_sigma;
        end

        function ll = score_image_blur(blur_sigma,PM)
            if (blur_sigma > PM.max_blur_sigma || blur_sigma < PM.min_blur_sigma)
                ll = -inf;
                return
            end
            rg = PM.max_blur_sigma - PM.min_blur_sigma;
            ll = -log(rg);
        end

        % ---
        % Pixel noise model (\epsilon)
        % ---
        function epsilon = sample_image_noise(PM)
            % sample uniformly in [min_epsilon, max_epsilon]
            rg = PM.max_epsilon - PM.min_epsilon;
            epsilon = rg*rand + PM.min_epsilon;
        end

        function ll = score_image_noise(epsilon,PM)
            if (epsilon > PM.max_epsilon || epsilon < PM.min_epsilon)
               ll = -inf;
               return
            end
            rg = PM.max_epsilon - PM.min_epsilon;
            ll = -log(rg);
        end

        % ---
        % Image model (I)
        % ---
        function I = sample_image(pimg)
        %  sample a binary image
            I = binornd(1,pimg);
            I = logical(I);
        end

        function ll = score_image(I,pimg)
        %  score the image model
            assert(~isempty(I));
            if islogical(I)
                on = I;
                off = ~I;
            elseif isnumeric(I) % if we having missing data
                missing = isinf(I);
                on = I > 0.5;
                off = I < 0.5;
                on(missing) = false;
                off(missing) = false;
            end
            prob_on = pimg;
            ll_on = sum(log(prob_on(on)));
            ll_off = sum(log(1-prob_on(off)));
            ll = ll_on + ll_off;
        end

%         function ll = score_image(I,pimg)
%         %  score the image model
%             assert(~isempty(I));
%             prob_on = pimg;
%             on = I;
%             off = ~I;
%             ll_on = sum(log(prob_on(on)));
%             ll_off = sum(log(1-prob_on(off)));
%             ll = ll_on + ll_off;
%         end

        % ---
        % Affine model (A)
        % ---
        function sample_A = sample_affine(libclass,nsamp)
            % affine transformation [x-scale,y-scale,x-translate,y-translate]
            % the translation is relative to the center of mass
            % (x-scale and y-scale cannot be 0 or negative)
            %
            % sample_A: [nsamp x 4] each row is a sampled affine warp
            if ~exist('nsamp','var')
                nsamp = 1;
            end

            sample_A = zeros(nsamp,4);

            % sample the image scale
            m_scale = libclass.affine.mu_scale;
            S_scale = libclass.affine.Sigma_scale;
            sample_A(:,1:2) = mvnrnd(m_scale,S_scale,nsamp);

            % sample the translation
            m_x = libclass.affine.mu_xtranslate;
            m_y = libclass.affine.mu_ytranslate;
            s_x = libclass.affine.sigma_xtranslate;
            s_y = libclass.affine.sigma_ytranslate;
            sample_A(:,3) = normrnd(m_x,s_x,nsamp,1);
            sample_A(:,4) = normrnd(m_y,s_y,nsamp,1);

            % scales should not be zero or negative
            if ~(all(vec(sample_A(:,1:2)) > 0))
               warning('sampled scale variable is less than zero');
            end

        end

        % NOTE: doesn't bother to truncate Gaussian properly,
        %       since the variance is so tight it won't matter
        function ll = score_affine(libclass,A)
            if isempty(A)
               A = [1; 1; 0; 0];
            end
            % score the affine warp
            % (the number of elements is 4)
            assert(isvector(A) && numel(A)==4);
            A = A(:);
            if (A(1) <= 0 || A(2) <= 0)
               ll = -inf;
               return
            end

            ll = 0;

            % sample the image scale
            m_scale = libclass.affine.mu_scale(:);
            S_scale = libclass.affine.Sigma_scale;
            ll = ll + mvnormpdfln(A(1:2),m_scale,[],S_scale);

            % sample the translation
            m_x = libclass.affine.mu_xtranslate;
            m_y = libclass.affine.mu_ytranslate;
            s_x = libclass.affine.sigma_xtranslate;
            s_y = libclass.affine.sigma_ytranslate;
            ll = ll + mvnormpdfln(A(3),m_x,s_x);
            ll = ll + mvnormpdfln(A(4),m_y,s_y);
        end

    end
end

function y = isunif(libclass)
    y = any(isnan(libclass.shape.mu));
end

% fast vectorized helper MVN likelihood
function ll = score_shape_type_helper(libclass,bspline_stack,subid,regcov)
    mymu = libclass.shape.mu(subid,:);
    mysd = libclass.shape.vsd(subid,:);
    if exist('regcov','var')
       mysd = sqrt( (mysd.^2) + regcov );
    end
    sz = size(mymu);
    vbspline = reshape(bspline_stack,sz(2),sz(1));
    vbspline = vbspline';
    llv = mvnormpdfln(vbspline(:)',mymu(:)',mysd(:)');
    llv = reshape(llv,sz);
    ll = sum(llv,2);
end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
function ta = acf(y,p)
% ACF - Compute Autocorrelations Through p Lags
% >> myacf = acf(y,p)
%
% Inputs:
% y - series to compute acf for, nx1 column vector
% p - total number of lags, 1x1 integer
%
% Output:
% myacf - px1 vector containing autocorrelations
%        (First lag computed is lag 1. Lag 0 not computed)
%
%
% A bar graph of the autocorrelations is also produced, with
% rejection region bands for testing individual autocorrelations = 0.
%
% Note that lag 0 autocorelation is not computed,
% and is not shown on this graph.
%
% Example:
% >> acf(randn(100,1), 10)
%

% --------------------------
% USER INPUT CHECKS
% --------------------------

[n1, n2] = size(y) ;
if n2 ~=1
    error('Input series y must be an nx1 column vector')
end

[a1, a2] = size(p) ;
if ~((a1==1 & a2==1) & (p<n1))
    error('Input number of lags p must be a 1x1 scalar, and must be less than length of series y')
end

% -------------
% BEGIN CODE
% -------------

ta = zeros(p,1) ;
global N
N = max(size(y)) ;
global ybar
ybar = mean(y);

% Collect ACFs at each lag i
for i = 1:p
   ta(i) = acf_k(y,i) ;
end

% Plot ACF
% Plot rejection region lines for test of individual autocorrelations
% H_0: rho(tau) = 0 at alpha=.05
bar(ta)
line([0 p+.5], (1.96)*(1/sqrt(N))*ones(1,2))
line([0 p+.5], (-1.96)*(1/sqrt(N))*ones(1,2))

% Some figure properties
line_hi = (1.96)*(1/sqrt(N))+.05;
line_lo = -(1.96)*(1/sqrt(N))-.05;
bar_hi = max(ta)+.05 ;
bar_lo = -max(ta)-.05 ;

if (abs(line_hi) > abs(bar_hi)) % if rejection lines might not appear on graph
    axis([0 p+.60 line_lo line_hi])
else
    axis([0 p+.60 bar_lo bar_hi])
end
title({' ','Sample Autocorrelations',' '})
xlabel('Lag Length')
set(gca,'YTick',[-1:.20:1])
% set number of lag labels shown
if (p<28 & p>4)
    set(gca,'XTick',floor(linspace(1,p,4)))
elseif (p>=28)
    set(gca,'XTick',floor(linspace(1,p,8)))
end
set(gca,'TickLength',[0 0])

% ---------------
% SUB FUNCTION
% ---------------
function ta2 = acf_k(y,k)
% ACF_K - Autocorrelation at Lag k
% acf(y,k)
%
% Inputs:
% y - series to compute acf for
% k - which lag to compute acf
%
global ybar
global N
cross_sum = zeros(N-k,1) ;

% Numerator, unscaled covariance
for i = (k+1):N
    cross_sum(i) = (y(i)-ybar)*(y(i-k)-ybar) ;
end

% Denominator, unscaled variance
yvar = (y-ybar)'*(y-ybar) ;

ta2 = sum(cross_sum) / yvar ;


text:
%% Load data
clear all; close all; clc;

nst_on = load('msckf_NST_on.mat');
nst_off = load('msckf_NST_off.mat');

kStart = 500; kEnd = 1000;

%% Compute RMSE
nst_off_trans_rmse = sqrt(mean(nst_off.msckf_trans_err.^2,1));
nst_off_rot_rmse = sqrt(mean(nst_off.msckf_rot_err.^2,1));

nst_on_trans_rmse = sqrt(mean(nst_on.msckf_trans_err.^2,1));
nst_on_rot_rmse = sqrt(mean(nst_on.msckf_rot_err.^2,1));

%% Plot stuff
figure(1); clf;
fontSize = 14;
lineWidth = 2;
pos = [200,200,640,400];
xLim = [kStart, kEnd];
k = kStart:kEnd;
xticks = linspace(kStart,kEnd,11);

transLim = [0 0.6];
rotLim = [0 0.25];

set(gcf,'Position',pos);

% Translational RMSE
subplot(2,1,1);
plot(k,nst_off_trans_rmse, '-r','LineWidth',lineWidth); hold on;
plot(k,nst_on_trans_rmse, '-b','LineWidth',lineWidth);
xlim(xLim);
set(gca,'XTick',xticks);
% ylim(transLim)
legend('Noise Correlated to State', 'Noise De-correlated from State', ...
    'Location','northwest');
ylabel('Trans. RMSE (m)')
title('Effect of Nullspace Projection');
set(gca,'FontSize',fontSize);
grid minor;
box on;

% Rotational RMSE
subplot(2,1,2);
plot(k,nst_off_rot_rmse,'-r','LineWidth',lineWidth); hold on;
plot(k,nst_on_rot_rmse,'-b','LineWidth',lineWidth);
xlim(xLim);
set(gca,'XTick',xticks);
% ylim(rotLim)
legend('Noise Correlated to State', 'Noise De-correlated from State', ...
    'Location','northwest');
ylabel('Rot. RMSE (Axis-Angle)')
xlabel('Timestep')
set(gca,'FontSize',fontSize)
grid minor;
box on;

%% Export figure
export_fig(gcf, 'NST_RMSE.pdf', '-transparent');

repo: MSCKF


text:
function K = sdlfmaXsdlfmaKernCompute(sdlfmaKern1, sdlfmaKern2, t1, t2, covIC)

% SDLFMAXSDLFMAKERNCOMPUTE Cross kernel between a SDLFMA and a SDLFMA kernels.
% FORMAT
% DESC computes cross kernel terms between the accel. and the accel. of
% two switching dynamical LFM kernels for the multiple output kernel.
% ARG sdlfmaKern1 : the kernel structure associated with the accel. of the
% first SDLFM kernel.
% ARG sdlfmaKern2 : the kernel structure associated with the accel. of the
% second SDLFM kernel.
% ARG t : inputs for which kernel is to be computed.
% ARG covIC : covariance for the initial conditions
% RETURN K : block of values from kernel matrix.
%
% FORMAT
% DESC computes cross kernel terms between the accel. and the accel. of
% two SDLFM kernels for the multiple output kernel.
% ARG sdlfmaKern1 : the kernel structure associated with the accel. of
% the first SDLFM kernel.
% ARG sdlfmaKern2 : the kernel structure associated with the accel. of
% the second SDLFM kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG covIC : covariance for the initial conditions
% RETURN K : block of values from kernel matrix.
%
% SEEALSO : sdlfmaKernParamInit, sdlfmaKernCompute, sdlfmKernParamInit
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if nargin < 4
    t2 = t1;
end

K = sdlfmXsdlfmKernCompute(sdlfmaKern1, sdlfmaKern2, t1, t2, covIC, 'AccelAccel');

author: SheffieldML

repo: GPmat


text:
function cfArray = ERBSpace(lowFreq, highFreq, N)
% function cfArray = ERBSpace(lowFreq, highFreq, N)
% This function computes an array of N frequencies uniformly spaced between
% highFreq and lowFreq on an ERB scale.  N is set to 100 if not specified.
%
% See also linspace, logspace, MakeERBCoeffs, MakeERBFilters.
%
% For a definition of ERB, see Moore, B. C. J., and Glasberg, B. R. (1983).
% "Suggested formulae for calculating auditory-filter bandwidths and
% excitation patterns," J. Acoust. Soc. Am. 74, 750-753.

if nargin < 1
	lowFreq = 100;
end

if nargin < 2
	highFreq = 44100/4;
end

if nargin < 3
	N = 100;
end

% Change the following three parameters if you wish to use a different
% ERB scale.  Must change in MakeERBCoeffs too.
EarQ = 9.26449;				%  Glasberg and Moore Parameters
minBW = 24.7;
order = 1;

% All of the followFreqing expressions are derived in Apple TR #35, "An
% Efficient Implementation of the Patterson-Holdsworth Cochlear
% Filter Bank."  See pages 33-34.
cfArray = -(EarQ*minBW) + exp((1:N)'*(-log(highFreq + EarQ*minBW) + ...
		log(lowFreq + EarQ*minBW))/N) * (highFreq + EarQ*minBW);

author: detly

repo: gammatone

save_path: github-repos/MATLAB/detly-gammatone


function seed(new_seed)
%SEED Set new seed for function RANDOM0
% Function SEED sets a new seed for function
% RANDOM0.  The new seed should be a positive
% integer.

% Define variables:
%   iseed    -- Random number seed (global)
%   new_seed -- New seed
%
%  Record of revisions:
%      Date       Programmer          Description of change
%      ====       ==========          =====================
%    02/04/07    S. J. Chapman        Original code

% Declare global values
global iseed            % Seed for random number generator

% Check for a legal number of input arguments.
msg = nargchk(1,1,nargin);
error(msg);

% Save seed
new_seed = round(new_seed);
iseed = abs(new_seed);


text:
% solve the problem:
%   min_{x,e} 0.5*|z-Dx-e|_2^2 + 0.5*lambda1*|x|_2^2 + lambda2*|e|_1
%
% solve the projection by APG
% input:
%   z - data point
%   D - basis matrix
%   lambda1, lambda2 - tradeoff parameters
% output:
%   r - projection coefficient
%   e - sparse noise
% copyright Jiashi Feng (jshfeng@gmail.com)
%
function [x,e] = solve_proj2(z,D,lambda1,lambda2)
  % initialization
  [ndim,ncol] = size(D);
  e = zeros(ndim,1);
  x = zeros(ncol,1);
  I = eye(ncol);
  converged = false;
  maxIter = inf;
  iter = 0;
  % alternatively update
  DDt = inv(D'*D+lambda1*I)*D';
  while ~converged
    iter = iter + 1;
    xtemp = x;
    x = DDt*(z-e);
    %     x = (D'*D + lambda1*I)\(D'*(z-e));
    etemp = e;
    e = thres(z-D*x,lambda2);
    stopc = max(norm(e-etemp), norm(x-xtemp))/ndim;
    if stopc < 1e-6 || iter > maxIter
      converged = true;
    end
    %     fval = func_proj(z,D,lambda1,lambda2,x,e);
    %     fprintf('fval = %f\n', fval);
  end
end

function x = thres(y,mu)
  x = max(y-mu, 0);
  x = x + min(y + mu, 0);
end

function fval = func_proj(z,D,lambda1,lambda2,x,e)
  fval = 0;
  fval = fval + 0.5*norm(z-D*x-e)^2;
  fval = fval + 0.5*lambda1*norm(x)^2;
  fval = fval + lambda2*sum(abs(e));
end

author: andrewssobral

repo: lrslibrary