text:
% fig2A: correlation to phototaxis stimulus
% for single regressor, detail
% option set to do correlation on trial-averaged traces
% plot corresponding histogram of corr.coeffs,
% and plot functional as well as anat plot

% compare to fig3A_sensory_motor, which for now plots a pair of regressors
% in a double colormap, and batch saves figures.

hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);
i_fish = 6;
ClusterIDs = [2,1];
[cIX,gIX,M,stim,behavior,M_0] = LoadSingleFishDefault(i_fish,hfig,ClusterIDs);

%% corr with PT regressor, stim-avr

setappdata(hfig,'isStimAvr',0);
[M_0,M,behavior,stim] = UpdateTimeIndex(hfig);

fishset = getappdata(hfig,'fishset');
[~,names,regressors] = GetStimRegressor(stim,fishset,i_fish);
reg_range = 2;
Reg = regressors(reg_range,:);
R = corr(Reg',M_0');

thres_reg = 0.5;
cIX = (find(R>thres_reg))';
wIX = R(cIX); % weight, i.e. corr.coeff
[~,I] = sort(wIX,'descend');
cIX = cIX(I);
gIX = ones(size(cIX));
wIX = wIX(I)';

UpdateIndices_Manual(hfig,cIX,gIX);

%% histogram of corr.coeff, for all cells
figure('Position',[500,200,500,150]);
hold on;
bins = -1:0.05:1;
[N,~] = histcounts(R,bins);
histogram(R,bins,'FaceColor',[0.4 0.4 0.4]);%,'EdgeColor','none'
plot([thres_reg,thres_reg],[0,max(N)],'r--');
xlim([-1,1]);ylim([0,max(N)]);
set(gca,'ycolor','w')

%% histogram with rainbow colors
clrbins = 0.5:0.05:1;
cmap = hot(length(clrbins)-1);
n_blank = length(bins)-length(clrbins);
cmap = vertcat(ones(n_blank,3),cmap);
figure('Position',[500,200,500,150]);
hold on
DrawRainbowHistogram(N,bins,cmap);
plot([thres_reg,thres_reg],[0,max(N)],'r--');
xlim([-1,1]);ylim([0,max(N)]);
set(gca,'ycolor','w')

%% left plot
figure('Position',[50,500,300,200])%[50,100,800,1000]);
% isCentroid,isPlotLines,isPlotBehavior,isPlotRegWithTS
setappdata(hfig,'isPlotBehavior',0);
setappdata(hfig,'isStimAvr',0);
UpdateTimeIndex(hfig);

% set color
clrmap = [0,0,0];
numK = max(gIX);
assert(size(clrmap,1)==numK);

opts = [];
opts.clrmap = clrmap;
opts.isAxTight = 1;
DrawTimeSeries(hfig,cIX,gIX,opts);

%% set up custom colormap
clrbins = 0.5:0.05:1;
clrmap0 = hot(round(1.5*length(clrbins)));
clrmap = clrmap0(end-length(clrbins):end,:);
gIX = interp1(clrbins,1:length(clrbins),wIX,'nearest');

%% right plot
I = LoadCurrentFishForAnatPlot(hfig,cIX,gIX);
I.clrmap = clrmap;
DrawCellsOnAnat(I);

%%
numTicks = 2;
DrawCustomColorbar(clrmap,clrbins,numTicks);

%% draw average plot (even if regression is done on full trace)
gIX = ones(size(cIX));
setappdata(hfig,'isStimAvr',1);
UpdateTimeIndex(hfig);
figure('Position',[50,500,300,200])%[50,100,800,1000]);
DrawTimeSeries(hfig,cIX,gIX,opts);

%%
[~,tot_image] = DrawCellsOnAnat(I);
ax = axes('Position',[0.8,0.8,0.05,0.15],'Units','normalized');
DrawCustomColorbar(clrmap,clrbins,2,ax);

%%

% %% left-right combined plot
% figure('Position',[50,100,1400,800]);
% % isCentroid,isPlotLines,isPlotBehavior,isPlotRegWithTS
% subplot(121)
% setappdata(hfig,'isPlotBehavior',0);
% setappdata(hfig,'isStimAvr',1);
% UpdateTimeIndex(hfig);
% DrawTimeSeries(hfig,cIX,gIX);
%
% % right plot
% subplot(122)
% I = LoadCurrentFishForAnatPlot(hfig,cIX,gIX);
% DrawCellsOnAnat(I);

%%
% dataDir = GetCurrentDataDir;
% saveDir = fullfile(dataDir,'motorsourceplot');
% if ~exist(saveDir, 'dir'), mkdir(saveDir), end;
% filename = fullfile(saveDir, num2str(i_fish));
% saveas(gcf, filename, 'png');
% close(gcf)

%%
%% motor regressor
%             [~,names,regressors] = GetMotorRegressor(behavior,i_fish);

author: xiuyechen

repo: FishExplorer


text:
function y = sumlogs(x,dim)
%SUMLOGS Sum of elements where numbers are represented by their logarithms.
%
%  Description
%    C=SUMLOGS(A) computes C=log(sum(exp(A))) in such a fashion
%    that it works even when elements have large magnitude.
%
%    C=SUMLOGS(A,DIM) sums along the dimension DIM.
%
%    See also SUM
%
% Copyright (c) 2013 Aki Vehtari

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

if nargin<2
  dim=find(size(x)>1,1);
end
maxx=max(x(:));
y=maxx+log(sum(exp(x-maxx),dim));

author: gpstuff-dev

repo: gpstuff

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


text:
function [node,elem,HB,bdFlag,indexMap,tree]= ...
         coarsen3(node,elem,markedElem,HB,bdFlag)
%% COARSEN3 coarse a 3-D triangulation.
%
% [node,elem,HB] = COARSEN3(node,elem,markedElem,HB) removes good-to-coraen
% nodes whose star are marked for coarsening. Unlike the 2-D version
% carsen, one additional data structure HB is needed. The local index of HB
% is 2 -- 1 -- 3, i.e., HB(:,1) is the middle point of the edge formed by
% HB(:,2:3). HB(:,4) is used to store the generation of the node HB(:,1).
%
% [node,elem,HB,bdFlag] = COARSEN3(node,elem,markedElem,HB,bdFlag) updates
% boundary conditions represented by bdFlag.
%
% [node,elem,HB,bdFlag,indexMap,tree] = COARSEN3(node,elem,markedElem,HB,bdFlag)
% outputs two additional information: indexMap and tree.
%
% - indexMap is the map between nodes in the fine mesh (node in the input)
%   to that in the coarse mesh (node in the output). For example,
%   indexMap(10) = 6 means the 10-th node in the fine grid  is now the 6-th
%   node in the coarse one. indexMap is useful for the interpolation of
%   nodewise function; see also nodeinterpolate
%
% - tree(:,1:3) stores the binary tree of the coarsening. tree(:,1) is the
%   index of parent element in coarsened mesh and tree(:,2:3) are two
%   children indices in original mesh.
%
%
% Example
%
%     [node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],2);
%     bdFlag = setboundary3(node,elem,'Dirichlet');
%     [node,elem,bdFlag,HB] = bisect3(node,elem,1,bdFlag,HB);
%     figure(2); subplot(1,2,1);
%     showmesh3(node,elem,[],'FaceAlpha',0.15); view([210 8]);
%     findnode3(node);
%     [node,elem,HB,bdFlag,tree,indexMap] = coarsen3(node,elem,'all',HB,bdFlag);
%     figure(2); subplot(1,2,2);
%     showmesh3(node,elem,[],'FaceAlpha',0.15); view([210 8]);
%
%
% coarsen3doc
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

tree = []; indexMap  = (1:size(node,1))'; % default output
if isempty(markedElem), return; end
if strcmp(markedElem,'all'), markedElem = (1:size(elem,1))'; end

%% Find good-to-coarsen nodes
N = size(node,1);
NT = size(elem,1);
generation = zeros(N,1);
generation(HB(:,1)) = HB(:,4);
valence = accumarray(elem(:),ones(4*NT,1),[N 1]);
valenceNew = accumarray(elem(:,4),ones(NT,1), [N 1]);
markedVal = accumarray(elem(markedElem,4),ones(length(markedElem),1),[N 1]);
isGoodNode=(valence==valenceNew) & (valenceNew==markedVal) & (generation>0);
if ~any(isGoodNode)
    return
end
t = find(isGoodNode(elem(:,4)));

%% Coarsen boundary faces
if (nargin==5) && ~isempty(bdFlag)
    bdElem = t(bdFlag(t,4)~=0);
    if ~isempty(bdElem)
        isBdGoodNode(elem(bdElem,4)) = true;
        bdElem = t(isBdGoodNode(elem(t,4)));
        T = auxstructure3(elem(bdElem,:));  % only the structure of bdElem
        leftNode = HB(elem(bdElem,4),2);
        for k = 1:3
            kidx = (elem(bdElem,k)==leftNode); % idx of bdElem
            tl = bdElem(kidx);                 % left element
            tr = bdElem(T.neighbor(kidx,k));   % matched neighbor
            bdFlag(tl,k) = bdFlag(tr,4);
        end
    end
else
	bdFlag=[];
end

%% Coarsen elements and record tree
leftNode = HB(elem(t,4),2);
idx = (elem(t,1)==leftNode) | (elem(t,2)==leftNode) | (elem(t,3)==leftNode);
tl = t(idx);
tr = t(~idx);
Nr = length(tl); % number of refined elements
tree = zeros(Nr,3);
% sort tl and tr by the common vertices such that tl and tr matches
if nargout == 6
    leftNode = HB(elem(tl,4),2);
    rightNode = HB(elem(tr,4),3);
    temptl = zeros(Nr,2);
    idx = (elem(tl,1) == leftNode);
    temptl(idx,:) = elem(tl(idx),[2 3]);
    idx = (elem(tl,2) == leftNode);
    temptl(idx,:) = elem(tl(idx),[1 3]);
    idx = (elem(tl,3) == leftNode);
    temptl(idx,:) = elem(tl(idx),[1 2]);
    temptl = sort(temptl,2);
    temptr = zeros(Nr,2);
    idx = (elem(tr,1) == rightNode);
    temptr(idx,:) = elem(tr(idx),[2 3]);
    idx = (elem(tr,2) == rightNode);
    temptr(idx,:) = elem(tr(idx),[1 3]);
    idx = (elem(tr,3) == rightNode);
    temptr(idx,:) = elem(tr(idx),[1 2]);
    temptr = sort(temptr,2);
    [temptl, Il] = sortrows([temptl elem(tl,4)]);
    [temptr, Ir] = sortrows([temptr elem(tr,4)]);
    tl = tl(Il);
    tr = tr(Ir);
    tree(:,3) = tr;
    tree(:,1) = tl;
    tree(:,2) = tl;
end
% coarsen tl to t
elem(tl,4) = HB(elem(tl,4),3);

%% Sort element nodes by generations
% the newest node is the node with maxmum generation
if (length(tl)==1)
    idx = max(transpose(generation(elem(t1,:))));
else
    [tempvar,idx] = max(generation(elem(tl,:)),[],2);  %#ok<ASGLU>
end
elem(tl((idx==1)),1:4) = elem(tl((idx==1)),[2 4 3 1]);
elem(tl((idx==2)),1:4) = elem(tl((idx==2)),[3 4 1 2]);
elem(tl((idx==3)),1:4) = elem(tl((idx==3)),[4 2 1 3]);
if (nargin==5) && ~isempty(bdFlag)
    bdFlag(tl((idx==1)),1:4) = bdFlag(tl((idx==1)),[2 4 3 1]);
    bdFlag(tl((idx==2)),1:4) = bdFlag(tl((idx==2)),[3 4 1 2]);
    bdFlag(tl((idx==3)),1:4) = bdFlag(tl((idx==3)),[4 2 1 3]);
end

%% Clean and shift index
elem(tr,:) = [];
inCoarse = true(NT,1);
inCoarse(tr) = false;
elemidxMap = zeros(NT,1);
elemidxMap(inCoarse) = 1:size(elem,1);
tree(:,1) = elemidxMap(tree(:,1));
if (nargin==5) && ~isempty(bdFlag)
    bdFlag(tr,:) = [];
end
node(isGoodNode,:) = [];
HB(isGoodNode,:) = [];
indexMap = zeros(N,1);
indexMap(~isGoodNode) = 1:size(node,1);
elem = indexMap(elem);
HB(:,1:3) = indexMap(HB(:,1:3));

repo: Aero-matlab


text:
function [label, energy] = knKmeansPred(model, Xt)
% Prediction for kernel kmeans clusterng
% Input:
%   model: trained model structure
%   Xt: d x n testing data
% Ouput:
%   label: 1 x n predict label
%   engery: optimization target value
% Written by Mo Chen (sth4nth@gmail.com).
X = model.X;
t = model.label;
kn = model.kn;

n = size(X,2);
k = max(t);
E = sparse(t,1:n,1,k,n,n);
E = bsxfun(@times,E,1./sum(E,2));
Z = bsxfun(@minus,E*kn(X,Xt),diag(E*kn(X,X)*E')/2);
[val, label] = max(Z,[],1);
energy = sum(kn(Xt))-2*sum(val);

author: epfl-lasa


text:
function [U,G,J,BC,SU,L] = slice_isolines(V,F,SV,val,varargin)
  % SLICE_ISOLINES Slice through a triangle mesh (V,F) at isolines val of a given
  % per-vertex scalar function.
  %
  % [U,G] = slice_isolines(V,F,SV,val);
  % [U,G,J,BC,SU] = slice_isolines(V,F,SV,val, ...
  %   'ParameterName',parameter_value, ...)
  %
  % Inputs:
  %   V  #V by dim list of tet mesh vertices
  %   F  #F by 3 list of tet indices into V
  %   S  #V list of scalar values per vertex
  %   val  #val list of isolines values
  %   Optional:
  %     'Manifold' followed by whether to stitch together triangles into a
  %       manifold mesh {false}: results in more compact U but slightly slower.
  % Outputs:
  %   U  #U by 3 list of triangle mesh vertices along slice
  %   G  #G by 3 list of triangles indices into U
  %   J  #G list of indices into F revealing which tet this face came from
  %   BC  #U by #V list of barycentric coordinates (or more generally: linear
  %     interpolation coordinates) so that U = BC*V
  %   SU  #U list of interpolated scalar values at U
  %   L  #G list of indices into val of next value above this triangle
  %

  manifold = false;

  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Manifold'}, ...
    {'manifold'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this
      % workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  % helper assuming val is a scalar
  function [U,G,J,BC,SU,L] = single_val(V,F,SV,val)
    function [U,G,BC,L] = one_below(V,F,SF,val)
      [sSF,sJ] = sort(SF,2);
      sF = F(sub2ind(size(F),repmat(1:size(F,1),size(F,2),1)',sJ));
      lambda = (sSF(:,2:3)-val)./bsxfun(@minus,sSF(:,2:3),sSF(:,1));
      BC = sparse( ...
        repmat((1:size(sF,1)*2)',1,2), ...
        [repmat(sF(:,1),2,1) reshape(sF(:,2:3),size(sF,1)*2,1)], ...
        [lambda(:) 1-lambda(:)], ...
        size(sF,1)*2,size(V,1));
      U = [V;BC * V];
      % Split into three triangles
      G = [ ...
        sF(:,1) size(V,1)+[1:size(F,1);size(F,1)+(1:size(F,1))]'; ...
        fliplr([sF(:,2) size(V,1)+[1:size(F,1);size(F,1)+(1:size(F,1))]']); ...
        sF(:,[2 3]) size(V,1)+[size(F,1)+(1:size(F,1))]'; ...
        ];
      L = [ones(size(sF,1),1);2*ones(2*size(sF,1),1)];
      flip = repmat( ...
        (sJ(:,1)==1 & sJ(:,2)==3) | ...
        (sJ(:,1)==2 & sJ(:,2)==1) | ...
        (sJ(:,1)==3 & sJ(:,2)==2),3,1);
      G(flip,:) = fliplr(G(flip,:));
    end
    SF = reshape(SV(F),size(F));
    I12 = sum(SF<val,2) == 1;
    % U is the running set of vertices
    U = V;
    [U,G12,BC12,L12] = one_below(U,F(I12,:),SF(I12,:),val);
    I21 = sum(SF>=val,2) == 1;
    [U,G21,BC21,L21] = one_below(U,F(I21,:),2*val-SF(I21,:),val);
    BC = [speye(size(V,1));BC12;BC21(:,1:size(V,1))];
    untouched = find(~I12&~I21);
    G = [F(untouched,:);G12;G21];
    J = [untouched;repmat(find(I12),3,1);repmat(find(I21),3,1)];
    SU = BC*SV;
    L = [(sum(SF(untouched,:)>=val,2)==3)+1;L12;3-L21];
  end

  U = V;
  G = F;
  SU = SV;
  J = 1:size(F,1);
  L = ones(size(F,1),1);
  BC = 1;
  for v = 1:numel(val)
    prev_J = J;
    prev_L = L;
    [U,G,J,BCC,SU,L] = single_val(U,G,SU,val(v));
    L = prev_L(J)+(L-1);
    J = prev_J(J);
    BC = BCC*BC;
  end

  if manifold
    % should be able to do this combinatorially
    bbd = normrow(max(V)-min(V));
    flag = [(1:size(V,1))';ones(size(U,1)-size(V,1),1)];
    [U,I,IM] = remove_duplicate_vertices([U flag],1e-14*bbd);
    U = U(:,1:3);
    BC = BC(I,:);
    SU = SU(I,:);
    G = IM(G);
  end

end

author: alecjacobson

repo: gptoolbox


text:
function signal = set_selepos(varargin)
% Select a set of epochs from a data set.
% Signal = set_selepos(Signal, EpochIndices)
%
% In:
%   Signal       : epoched EEGLAB data set from which epochs shall be selected
%
%   EpochIndices : indices of the epochs that should be retained
%
% Out:
%   Signal       : Newly created data set consisting of the retained epochs in the order
%                  specified in EpochIndices.
%
% Examples:
%   % for an epoched data set, select epochs 50-100
%   eeg = set_selepos(eeg,50:100)
%
% See also:
%   set_selinterval, set_partition, set_makepos
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-01
dp;

% set_selepos_version<2.0> -- for the cache

if ~exp_beginfun('editing') return; end

declare_properties('independent_channels',true,'independent_trials',false);

arg_define(varargin, ...
    arg_norep({'signal','Signal'}), ...
    arg({'epoch_range','EpochIndices','epos'},[],uint32([1 1000000]),'Indices of retained epochs.','shape','row'));

% input validation
utl_check_fields(signal,{'epoch','event','data','pnts'},'signal','signal'); %#ok<NODEF>
if isempty(signal.epoch) && size(signal.data,3) > 1
    error('Your signal has epoched data but has an empty .epoch field. This is not permitted. Use set_makepos instead of pop_epoch for epoching in BCILAB.'); end

if ~isempty(signal.epoch)
    if ~isfield(signal.epoch,'event')
        error('The given signal is lacking the required .epoch.event field.'); end
    if isa(epoch_range,'logical') %#ok<NODEF>
        epoch_range = find(epoch_range); end
    if ~isreal(epoch_range)
        error('The given epoch range contains data that is not real-valued.'); end
    if length(signal.epoch) < max(epoch_range)
        error('The epoch range to select exceeds the length of the signal''s .epoch field (%i): %s.',length(signal.epoch),hlp_tostring(epoch_range)); end

    % select epochs within time series fields
    for field = utl_timeseries_fields(signal)
        if ~isempty(signal.(field{1}))
            try
                signal.(field{1}) = signal.(field{1})(:,:,epoch_range,:,:,:,:,:);
            catch e
                error('The given epoch indices could not be extracted from the time-series field .%s with error: %s. Field size was: %i, epoch indices were: %s',field{1},e.message,size(signal.(field{1}),3),hlp_tostring(epoch_range));
            end
        end
    end

    % reindex .epoch field
    signal.epoch = signal.epoch(epoch_range);
    % reindex .event field
    signal.event = signal.event([signal.epoch.event]);
    if ~isempty(epoch_range)
        % get new .event.epoch indices
        epoch_numevents = cellfun('length',{signal.epoch.event});
        for j=length(epoch_numevents):-1:1
            tmp{j} = ones(1,epoch_numevents(j))*j; end;
        epoch_indices = [tmp{:}];
        if ~isempty(signal.event)
            % validate signal.event
            if ~isfield(signal.event,'latency')
                error('The given signal is lacking the required .event.latency field.'); end
            if ~isfield(signal.event,'epoch')
                error('The given signal is lacking the required .event.epoch field.'); end
            latency_numels = cellfun('prodofsize',{signal.event.latency});
            if any(latency_numels == 0)
                error('The given signal has one or more events with empty .latency field. This is not permitted.');
            elseif any(latency_numels ~= 1)
                error('The given signal has one or more events with a .latency value that is not a scalar. This is not permitted.');
            end
            % rewrite pseudo-continuous event latencies based on epoch index change
            [signal.event.latency] = arraydeal([signal.event.latency] + signal.pnts*(epoch_indices - [signal.event.epoch]));
            % relink events and epochs
            [signal.event.epoch] = arraydeal(epoch_indices);
        end
        [signal.epoch.event] = chopdeal(1:length(signal.event),epoch_numevents);
    end
    % update misc meta-data
    signal.trials = size(signal.data,3);
    if ~signal.trials
        disp_once('WARNING: After set_selepos your signal is empty.'); end
elseif ~isempty(epoch_range)
    error('The given epochs could not be selected since the signal is empty.');
end

exp_endfun;

repo: REST


text:
function y2 = y_astronomical_to_common ( y )

%*****************************************************************************80
%
%% Y_ASTRONOMICAL_TO_COMMON converts an Astronomical year to a Common year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the astronomical year.
%
%    Output, integer Y2, the Common year.
%
  if ( y <= 0 )
    y2 = y - 1;
  else
    y2 = y;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [K] = spm_perm_mtx(n)
% Returns a matrix of indices permuted over n
% FORMAT [K] = spm_perm_mtx(n)
%    n   - (scalar) number of indices
%    K   - (2^n x n) permutation matrix
%    n   - (vector) indices
%    K   - (length(n)! x n) permutation matrix
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_perm_mtx.m 5657 2013-09-26 16:53:40Z karl $

% get permutations
%==========================================================================

% permute zeros and ones
%--------------------------------------------------------------------------
if isscalar(n)

    N  = 2^n;
    K  = sparse(N,n);
    x  = sparse(1,1,1,2,1);
    for i = 1:n
        y      = ones(N/length(x),1);
        K(:,i) = kron(x,y);
        x      = [x;x];
    end

% permute indices
%--------------------------------------------------------------------------
elseif isvector(n)

    n  = n(:);
    K  = n;
    while size(K,2) < length(n)
        x     = [];
        for i = 1:size(K,1)
            d = K(i,:);
            r = n;
            for j = 1:length(d)
                r(r == d(j)) = [];
            end
            x = [x; [kron(ones(length(r),1),d) r]];
        end
        K = x;
    end
end

% make logical
%--------------------------------------------------------------------------
K  = logical(K);

author: spm


text:
function b = isPointOnLine3d(point, line, varargin)
%ISPOINTONLINE3D Test if a 3D point belongs to a 3D line.
%
%   B = isPointOnLine3d(POINT, LINE)
%   with POINT being [xp yp zp], and LINE being [x0 y0 z0 dx dy dz].
%   Returns 1 if point lies on the line, 0 otherwise.
%
%   If POINT is an N-by-3 array of points, B is a N-by-1 array of booleans.
%
%   If LINE is a N-by-6 array of lines, B is a N-by-1 array of booleans.
%
%   B = isPointOnLine3d(POINT, LINE, TOL)
%   Specifies the tolerance used for testing location on 3D line.
%
%   See also
%   lines3d, distancePointLine3d, linePosition3d, isPointOnLine
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2003-10-31
% Copyright 2003-2022 INRA - TPV URPOI - BIA IMASTE

% extract computation tolerance
tol = 1e-14;
if ~isempty(varargin)
    tol = varargin{1};
end

% size of inputs
np = size(point,1);
nl = size(line, 1);

if np == 1 || nl == 1 || np == nl
    % test if lines are colinear, using norm of the cross product
    b = bsxfun(@rdivide, vectorNorm3d( ...
        crossProduct3d(bsxfun(@minus, line(:,1:3), point), line(:,4:6))), ...
        vectorNorm3d(line(:,4:6))) < tol;
else
    % same test, but after reshaping arrays to manage difference of
    % dimensionality
    point = reshape(point, [np 1 3]);
    line = reshape(line, [1 nl 6]);
    b = bsxfun(@rdivide, vectorNorm3d( ...
        cross(bsxfun(@minus, line(:,:,1:3), point), line(ones(1,np),:,4:6), 3)), ...
        vectorNorm3d(line(:,:,4:6))) < tol;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function node_num = grid_ql_node_num ( nelemx, nelemy )

%*****************************************************************************80
%
%% GRID_QL_NODE_NUM counts the nodes in a grid of QL elements.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NELEMX, NELEMY, the number of elements along the
%    X and Y directions.  The number of elements generated will be
%    NELEMX * NELEMY.
%
%    Output, integer NODE_NUM, the number of nodes in the grid.
%
  node_num = 2 * nelemx * nelemy + 2 * nelemx + nelemy + 1;

  return
end


text:
function varargout=importsii(varargin)
%
% For calling details please see v_importsii.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_importsii.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: importsii.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_importsii(varargin{:});
else
    v_importsii(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
% forcecol.m
% function to force a vector to be a single column
%
% Brian Birge
% Rev 1.0
% 7/1/98

function[out]=forcecol(in)
len=prod(size(in));
out=reshape(in,[len,1]);

author: HuangCongQing


text:
close all; clear all; clc

addpath('apm')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Configuration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% number of terms
ny = 2; % output coefficients
nu = 1; % input coefficients
% number of inputs
ni = 1;
% number of outputs
no = 1;
% load data and parse into columns
data = load('data_step_test.csv');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% generate time-series model
sysd = apm_id(data,ni,nu,ny);

% plt.figure(1)
% plt.subplot(2,1,1)
% plt.plot(data[:,0],ypred,'r-',LineWidth=2)
% plt.plot(data[:,0],data[:,2],'b--',LineWidth=2)
% plt.legend(['Predicted','Measured'],loc='best')
% plt.ylabel('Temp (K)')
%
% plt.subplot(2,1,2)
% plt.plot(data[:,0],data[:,1],'k-',LineWidth=2)
% plt.legend(['Heater'],loc='best')
% plt.ylabel('Heater')
% plt.show()

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


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

    % load gcTrajData
    [id, gc_x, gc_y, gc_z] = load_gc_traj_data(gcTrajFile);

    subplot(1, 1, 1,'replace');

    % draw cube
    edge1_x = [0 0];
    edge1_y = [0 0];
    edge1_z = [0 36];

    edge2_x = [0 0];
    edge2_y = [0 36];
    edge2_z = [0 0];

    edge3_x = [0 36];
    edge3_y = [0 0];
    edge3_z = [0 0];

    edge4_x = [0 36];
    edge4_y = [36 36];
    edge4_z = [0 0];

    edge5_x = [0 0];
    edge5_y = [36 36];
    edge5_z = [0 36];

    edge6_x = [36 36];
    edge6_y = [0 36];
    edge6_z = [0 0];

    edge7_x = [36 36];
    edge7_y = [0 0];
    edge7_z = [0 36];

    edge8_x = [0 0];
    edge8_y = [0 36];
    edge8_z = [36 36];

    edge9_x = [0 36];
    edge9_y = [0 0];
    edge9_z = [36 36];

    edge10_x = [0 36];
    edge10_y = [36 36];
    edge10_z = [36 36];

    edge11_x = [36 36];
    edge11_y = [0 36];
    edge11_z = [36 36];

    edge12_x = [36 36];
    edge12_y = [36 36];
    edge12_z = [0 36];

    plot3(edge1_x,edge1_y,edge1_z, '-k');
    hold on
    plot3(edge2_x,edge2_y,edge2_z, '-k');
    hold on
    plot3(edge3_x,edge3_y,edge3_z, '-k');
    hold on
    plot3(edge4_x,edge4_y,edge4_z, '-k');
    hold on
    plot3(edge5_x,edge5_y,edge5_z, '-k');
    hold on
    plot3(edge6_x,edge6_y,edge6_z, '-k');
    hold on
    plot3(edge7_x,edge7_y,edge7_z, '-k');
    hold on
    plot3(edge8_x,edge8_y,edge8_z, '-k');
    hold on
    plot3(edge9_x,edge9_y,edge9_z, '-k');
    hold on
    plot3(edge10_x,edge10_y,edge10_z, '-k');
    hold on
    plot3(edge11_x,edge11_y,edge11_z, '-k');
    hold on
    plot3(edge12_x,edge12_y,edge12_z, '-k');

    view(20,20)

    % draw gc trajectory

    hold on;
    plot3(gc_y, gc_x, gc_z, '.r', 'MarkerSize',6);

    hold off;
    xl = xlabel('y', 'FontSize',24);
    yl = ylabel('x', 'FontSize',24);
    z1 = zlabel('z', 'FontSize',24);
    view(33,12)

%     set(xl,'Rotation',15);
    set(yl,'Rotation',30);

    axis([0 36 0 36 0 36]);

    set(gca,'xtick',0:9:36)
    set(gca,'ytick',0:18:36)
    set(gca,'ztick',0:9:36)
    set(gca,'FontSize',24, 'LineWidth',1.5); % axis font

    fig = get(groot,'CurrentFigure');
    set (fig,'Position',[500,300,400,300], 'color','w')
end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function pdflib_test01 ( )

%*****************************************************************************80
%
%% PDFLIB_TEST01 tests R8MAT_POFAC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2013
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PDFLIB_TEST01\n' );
  fprintf ( 1, '  R8MAT_POFAC computes the Cholesky factor R of a\n' );
  fprintf ( 1, '  positive definite matrix A, so that A = R'' * R.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Start with random R1;\n' );
  fprintf ( 1, '  Compute A = R1'' * R1.\n' );
  fprintf ( 1, '  Call R8MAT_POFAC and see if you recover R2 = R1.\n' );
%
%  Generate a random upper triangular matrix with positive diagonal.
%
  r1 = zeros ( n, n );

  for j = 1 : n
    for i = 1 : j
       r1(i,j) = r8_uniform_01_sample ( );
    end
  end

  r8mat_print ( n, n, r1, '  R1:' );
%
%  Compute a positive definite symmetric matrix A.
%
  a = r1' * r1;

  r8mat_print ( n, n, a, '  A:' );

  r2 = r8mat_pofac ( n, a );

  diff = r8mat_norm_fro_affine ( n, n, r1, r2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Frobenius difference between R1 and R2 = %g\n', diff );

  return
end


text:
function [st3, rez] = standalone_detector(rez, spkTh)
% Detects spikes across the entire recording using generic templates.
% Each generic template has rank one (in space-time).
% In time, we use the 1D template prototypes found in wTEMP.
% In space, we use Gaussian weights of several sizes, centered
% on (x,y) positions that are part of a super-resolution grid covering the
% entire probe (pre-specified in the calling function).
% In total, there ~100x more generic templates than channels.

NrankPC = size(rez.wPCA, 2);
ops = rez.ops;

% minimum/base sigma for the Gaussian.
sig = 10;

% grid of centers for the generic tempates
[ycup, xcup] = meshgrid(ops.yup, ops.xup);

% determine prototypical timecourses by clustering of simple threshold crossings.
wTEMP = gpuArray(rez.wTEMP);
wPCA = gpuArray(rez.wPCA);

% Get nearest channels for every template center.
% Template products will only be computed on these channels.
NchanNear = 10;
[iC, dist] = getClosestChannels2(ycup, xcup, rez.yc, rez.xc, NchanNear);

% Templates with centers that are far from an active site are discarded
dNearActiveSite = median(diff(unique(rez.yc)));
igood = dist(1,:)<dNearActiveSite;
iC = iC(:, igood);
dist = dist(:, igood);
ycup = ycup(igood);
xcup = xcup(igood);

% number of nearby templates to compare for local template maximum
NchanNearUp =  min(numel(xcup), 10*NchanNear);
[iC2, dist2] = getClosestChannels2(ycup, xcup, ycup, xcup, NchanNearUp);

% pregenerate the Gaussian weights used for spatial components
nsizes = 5;
v2 = gpuArray.zeros(nsizes, size(dist,2), 'single');
for k = 1:nsizes
    v2(k, :) = sum(exp( - 2 * dist.^2 / (sig * k)^2), 1);
end

% build up Params
NchanUp = size(iC,2);

% preallocate the results
st3 = zeros(1000000, 6);
t0 = 0; %ceil(rez.ops.trange(1) * ops.fs); % I think this should be 0 all the time.
nsp = 0; % counter for total number of spikes
%%
tic
for k = 1:ops.Nbatch
    % get a batch of whitened and filtered data
    dataRAW = get_batch(rez.ops, k);

    Params = [size(dataRAW,1) ops.Nchan ops.nt0 NchanNear NrankPC ops.nt0min spkTh NchanUp NchanNearUp sig];

    % run the CUDA function on this batch
    [dat, kkmax, st, cF] = spikedetector3(Params, dataRAW, wTEMP, iC-1, dist, v2, iC2-1, dist2);

    % upsample the y position using the center of mass of template products
    % coming out of the CUDA function.
    ys = rez.yc(iC);
    cF0 = max(0, cF);
    cF0 = cF0 ./ sum(cF0, 1);
    iChan = st(2, :) + 1;
    yct = sum(cF0 .* ys(:, iChan), 1);

    % build st for the current batch
    st = double(gather(st));
    st(6, :) = st(2,:);
    st(2,:) = gather(yct);

    toff = ops.nt0min + t0 + ops.NT*(k-1);
    st(1,:) = st(1,:) + toff; % these offsets ensure the times are computed correctly

    st(5,:) = k; % batch number

    nsp0 = size(st,2);
    if nsp0 + nsp > size(st3,1)
       st3(nsp + 1e6, 1) = 0; % if we need to preallocate more space
    end

    st3(nsp + [1:nsp0], :) = st';
    nsp = nsp + nsp0;

    if rem(k,100)==1 || k==ops.Nbatch
        fprintf('%2.2f sec, %d batches, %d spikes \n', toc, k, nsp)
    end
end

st3 = st3(1:nsp, :);
%%

rez.iC = gather(iC);
rez.dist =  gather(dist);

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function varargout = transformVector3d(varargin)
%TRANSFORMVECTOR3D Transform a vector with a 3D affine transform.
%
%   V2 = transformVector3d(V1, TRANS);
%   Computes the vector obtained by transforming vector V1 with affine
%   transform TRANS.
%   V1 has the form [x1 y1 z1], and TRANS is a [3x3], [3x4], or [4x4]
%   matrix, with one of the forms:
%   [a b c]   ,   [a b c j] , or [a b c j]
%   [d e f]       [d e f k]      [d e f k]
%   [g h i]       [g h i l]      [g h i l]
%                                [0 0 0 1]
%
%   V2 = transformVector3d(V1, TRANS) also works when V1 is a [Nx3xMxEtc]
%   array of double. In this case, V2 has the same size as V1.
%
%   V2 = transformVector3d(X1, Y1, Z1, TRANS);
%   Specifies vectors coordinates in three arrays with same size.
%
%   [X2 Y2 Z2] = transformVector3d(...);
%   Returns the coordinates of the transformed vector separately.
%
%
%   See also
%   vectors3d, transforms3d, transformPoint3d

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2008-11-25, from transformPoint3d
% Copyright 2008-2022 INRA - TPV URPOI - BIA IMASTE

if nargin~=2 && nargin~=4
    error('Invalid number of input arguments. Type ''help transformVector3d'' for details.');
end

% Extract only the linear part of the affine transform
trans = varargin{end};
trans(1:4,4) = [0; 0; 0; 1];

% Call transformPoint3d using equivalent output arguments
varargout = cell(1, max(1,nargout));
[varargout{:}] = transformPoint3d(varargin{1:end-1}, trans);


text:
function test_laggedcoherence

% MEM 2gb
% WALLTIME 00:20:00
% DEPENDENCY ft_connectivityanalysis ft_laggedcoherence ft_freqanalysis

% this is a test for the low-level function

datain=[];
datain.label={'Ch1','Ch2'};
datain.fsample=1000;
datain.time{1}=0:1/datain.fsample:8.832; % changed by JM, since it does not make sense to have less than 6 s of data if 1 Hz is requested with 3 cycles
datain.trial{1}=randn(2,length(datain.time{1}));

cfg=[];
cfg.foi=[1 1.5 2 3 4 5 6 7 8 17 18 19 20 21];
cfg.loi=[1 2 3];
cfg.numcycles=3;

lcout=ft_laggedcoherence(cfg,datain);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = edgesEval( model, varargin )
% Run and evaluate structured edge detector on BSDS500.
%
% This function first runs the trained structured edge detector on every
% test or validation image in BSDS then call edgesEvalDir.m to perform the
% actual edge evaluation. edgesEval is specific to the structured edge
% detector and BSDS, edgesEvalDir is general purpose edge evaluation code.
% For example usage of edgesEval see edgesDemo.
%
% USAGE
%  varargout = edgesEval( model, prms )
%
% INPUTS
%  model      - structured edge model trained with edgesTrain
%  prms       - parameters (struct or name/value pairs)
%   .dataType   - ['test'] should be either 'test' or 'val'
%   .name       - [''] name to append to evaluation
%   .opts       - {} list of model opts to overwrite
%   .show       - [0] if true plot results using edgesEvalPlot
%   .pDistr     - [{'type','parfor'}] parameters for fevalDistr
%   .cleanup    - [0] if true delete temporary files
%   .thrs       - [99] number or vector of thresholds for evaluation
%   .maxDist    - [.0075] maximum tolerance for edge match
%
% OUTPUTS
%  varargout  - same outputs as edgesEvalDir
%
% EXAMPLE
%
% See also edgesDemo, edgesDetect, edgesEvalDir, edgesEvalPlot
%
% Structured Edge Detection Toolbox      Version 3.01
% Code written by Piotr Dollar, 2014.
% Licensed under the MSR-LA Full Rights License [see license.txt]

% get default parameters
dfs={'dataType','test', 'name','', 'opts',{}, 'show',0, ...
  'pDistr',{{'type','parfor'}}, 'cleanup',0, 'thrs',99, 'maxDist',.0075 };
p=getPrmDflt(varargin,dfs,1);

% load model and update model.opts acoording to opts
if( ischar(model) ), model=load(model); model=model.model; end
for i=1:length(p.opts)/2, model.opts.(p.opts{i*2-1})=p.opts{i*2}; end
rgbd=model.opts.rgbd; model.opts.nms=1;

% get list of relevant directories (image, depth, gt, results)
name = [model.opts.modelFnm,p.name];
imgDir = fullfile(model.opts.bsdsDir,'images',p.dataType);
depDir = fullfile(model.opts.bsdsDir,'depth',p.dataType);
gtDir  = fullfile(model.opts.bsdsDir,'groundTruth',p.dataType);
resDir = fullfile(model.opts.modelDir,p.dataType,name);
assert(exist(imgDir,'dir')==7); assert(exist(gtDir,'dir')==7);

% run edgesDetect() on every image in imgDir and store in resDir
if( ~exist(fullfile([resDir '-eval'],'eval_bdry.txt'),'file') )
  if(~exist(resDir,'dir')), mkdir(resDir); end
  ids=dir(imgDir); ids=ids([ids.bytes]>0); ids={ids.name}; n=length(ids);
  ext=ids{1}(end-2:end); for i=1:n, ids{i}=ids{i}(1:end-4); end
  res=cell(1,n); for i=1:n, res{i}=fullfile(resDir,[ids{i} '.png']); end
  do=false(1,n); for i=1:n, do(i)=~exist(res{i},'file'); end
  ids=ids(do); res=res(do); m=length(ids);
  parfor i=1:m, id=ids{i};
    I = imread(fullfile(imgDir,[id '.' ext])); D=[];
    if(rgbd), D=single(imread(fullfile(depDir,[id '.png'])))/1e4; end
    if(rgbd==1), I=D; elseif(rgbd==2), I=cat(3,single(I)/255,D); end
    E=edgesDetect(I,model); imwrite(uint8(E*255),res{i});
  end
end

% perform actual evaluation using edgesEvalDir
varargout=cell(1,max(1,nargout));
[varargout{:}] = edgesEvalDir('resDir',resDir,'gtDir',gtDir,...
  'pDistr',p.pDistr,'cleanup',p.cleanup,'thrs',p.thrs,'maxDist',p.maxDist);
if( p.show ), figure(p.show); edgesEvalPlot(resDir,name); end

end

author: pdollar

repo: edges

save_path: github-repos/MATLAB/pdollar-edges


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2014 the National Renewable Energy Laboratory and Sandia Corporation
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Run the WEC-Sim RM3 Test Case for 1DOF w/PTO
% This script will run the RM3 in WEC-Sim for irregular waves with Hs=2.5[m]
% and Tp=8[s]. The RM3 WEC models has 2 bodies, restricted to heave motion
% only, and has PTO damping=1200[kN-s/m].

global plotNO
locdir=pwd;
%% Run Simulation
wecSim;                     % Run Simulation

%% Post-Process Data
% Body 1
B1.WEC_Sim_new.time=output.bodies(1).time;
Hshift1=output.bodies(1).position(1,3);
B1.WEC_Sim_new.heave=output.bodies(1).position(:,3)-Hshift1;
%Body 2
B2.WEC_Sim_new.time=output.bodies(2).time;
Hshift2=output.bodies(2).position(1,3);
B2.WEC_Sim_new.heave=output.bodies(2).position(:,3)-Hshift2;
% Relative Motion
Rel.WEC_Sim_new.time=B1.WEC_Sim_new.time;
Rel.WEC_Sim_new.heave=B1.WEC_Sim_new.heave-B2.WEC_Sim_new.heave;
% Spectrum
Sp.WEC_Sim_new.m0 = calcSpectralMoment(waves.omega,waves.spectrum,0);
Sp.WEC_Sim_new.m2 = calcSpectralMoment(waves.omega,waves.spectrum,2);

%% Load Data
irrSSoutput_new=output;        % Keeps new run in workspace
load('irregularSS_org.mat')    % Load Previous WEC-Sim Data

%% Post-Process Data
% Body 1
B1.WEC_Sim_org.time=output.bodies(1).time;
Hshift1=output.bodies(1).position(1,3);
B1.WEC_Sim_org.heave=output.bodies(1).position(:,3)-Hshift1;
%Body 2
B2.WEC_Sim_org.time=output.bodies(2).time;
Hshift2=output.bodies(2).position(1,3);
B2.WEC_Sim_org.heave=output.bodies(2).position(:,3)-Hshift2;
% Relative Motion
Rel.WEC_Sim_org.time=B1.WEC_Sim_org.time;
Rel.WEC_Sim_org.heave=B1.WEC_Sim_org.heave-B2.WEC_Sim_org.heave;
% Spectrum
Sp.WEC_Sim_org.m0 = calcSpectralMoment(waves.omega,waves.spectrum,0);
Sp.WEC_Sim_org.m2 = calcSpectralMoment(waves.omega,waves.spectrum,2);

clear Hshift1 Hshift2

%% Quantify Maximum Difference Between Saved and Current WEC-Sim Runs
[irregularSS.B1_H_max ,irregularSS.B1_H_I]=max(abs(B1.WEC_Sim_org.heave-B1.WEC_Sim_new.heave));
[irregularSS.B2_H_max ,irregularSS.B2_H_I]=max(abs(B2.WEC_Sim_org.heave-B2.WEC_Sim_new.heave));
[irregularSS.Rel_H_max,irregularSS.Rel_H_I]=max(abs(Rel.WEC_Sim_org.heave-Rel.WEC_Sim_new.heave));

irregularSS.B1  = B1;
irregularSS.B2  = B2;
irregularSS.Rel = Rel;
irregularSS.Sp  = Sp;

save('irregularSS','irregularSS')

%% Plot Old vs. New Comparison
if plotNO==1 % global variable
    cd ../..
    plotOldNew(B1,B2,Rel,[irregularSS.B1_H_max ,irregularSS.B1_H_I],...
        [irregularSS.B2_H_max ,irregularSS.B2_H_I],...
        [irregularSS.Rel_H_max,irregularSS.Rel_H_I],'IrregularSS');
    cd(locdir)
end
%% Clear output and .slx directory

try
	rmdir('output','s')
	rmdir('slprj','s')
	delete('RM3.slx.autosave', 'RM3_sfun.mexmaci64')
catch
end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function ehmm = updateAlpha_ehmm(ehmm,rangeK)
if nargin < 2 || isempty(rangeK), rangeK = 1:ehmm.K; end % K+1
ehmm = updateAlpha(ehmm,rangeK);
% setstateoptions;
% for k = 1:ehmm.K
%     for n = 1:ehmm.train.ndim
%         if ~regressed(n), continue; end
%         indr = S(n,:)==1;
%         for i = 1:length(orders)
%             index = (i-1)*ehmm.train.ndim + n + ~train.zeromean;
%             if ehmm.train.ndim > 1
%                 sigmaprior = (ehmm.state(k).sigma.Gam_shape(n,indr) ./ ...
%                     ehmm.state(k).sigma.Gam_rate(n,indr));
%             else
%                 sigmaprior = zeros(1,sum(indr));
%             end
%             ehmm.state(k).alpha.Gam_rate(i) = ehmm.state(k).alpha.Gam_rate(i) + ...
%                 0.5 * ( (ehmm.state(k).W.Mu_W(index,indr) .* ...
%                 sigmaprior ) * ...
%                 ehmm.state(k).W.Mu_W(index,indr)' + sum( ...
%                 sigmaprior .* ...
%                 ehmm.state(k).W.S_W(indr,index,index)'));
%         end
%         ehmm.state(k).alpha.Gam_shape = ehmm.state(k).alpha.Gam_shape + 0.5 * sum(indr);
%     end
% end
end

author: OHBA-analysis

repo: HMM-MAR


text:
function [v1, v2] = lambertBattinVector(r1, r2, dt, numRevs, gmu)
%lambertBattin Summary of this function goes here
%   Detailed explanation goes here
    if(size(r1,1) ~= 3)
        error('r1 not of length 3');
    end

    if(size(r2,1) ~= 3)
        error('r2 not of length 3');
    end

    r1Mag = sqrt(sum(abs(r1).^2,1));
    r2Mag = sqrt(sum(abs(r2).^2,1));

    tm = zeros(size(dt));
    if(any(dt(dt<0)))
        B = dt<0;
        tm(B) = -1.0;
    end
    if(any(dt(dt>0)))
        B = dt>0;
        tm(B) = 1.0;
    end
    if(any(dt(dt==0)))
        error('dt = 0; lambert cannot compute');
    end
    dt = abs(dt);

    cosDeltaTA = dot(r1,r2,1)./(r1Mag.*r2Mag);
    sinDeltaTA = tm .* sqrt(1 - cosDeltaTA.^2);
    deltaTA = atan2(sinDeltaTA,cosDeltaTA);
    deltaTA = AngleZero2Pi(deltaTA);

    c = sqrt(r1Mag.^2 + r2Mag.^2 - 2.*r1Mag.*r2Mag.*cosDeltaTA);
    s = (r1Mag + r2Mag + c)./2;
    epsilon = (r2Mag - r1Mag)./r1Mag;

    intemed1 = r2Mag./r1Mag;
    TanSqr2w = (epsilon.^2/4) ./ (sqrt(intemed1) + intemed1.*(2 + sqrt(intemed1)));

    sinSqrDeltaTAOver4 = (sin(deltaTA/4)).^2;
    cosSqrDeltaTAOver4 = (cos(deltaTA/4)).^2;
    rop = sqrt(r1Mag.*r2Mag) .* (cosSqrDeltaTAOver4 + TanSqr2w);

    l = zeros(size(tm));
    if(any(tm(tm==1.0)))
        bool = tm==1.0;
        l(bool) = (sinSqrDeltaTAOver4(bool) + TanSqr2w(bool))./(sinSqrDeltaTAOver4(bool) + TanSqr2w(bool) + cos(deltaTA(bool)/2));
    end
    if(any(tm(tm==-1.0)))
        bool = tm==-1.0;
        l(bool) = (cosSqrDeltaTAOver4(bool) + TanSqr2w(bool) - cos(deltaTA(bool)/2)) ./ (cosSqrDeltaTAOver4(bool) + TanSqr2w(bool));
    end

    m = (gmu .* dt.^2)./(8*rop.^3);

    x = l;
    x_change = 1;
    loops=0;
    while(any(x_change > 1E-6) && loops <= 30)
            ksi = computeKsi(x, 20);

            h1 = ((l+x).^2.*(1 + 3.*x + ksi)) ./ ((1+2.*x+l).*(4.*x + ksi.*(3+x)));
            h2 = (m.*(x - l + ksi)) ./ ((1+2.*x+l).*(4.*x + ksi.*(3+x)));

            B=27*h2./(4*(1+h1).^3);
            U=B./(2*(sqrt(1+B)+1));
            K_U=Kay(U);
            y=(1+h1)./3.*(2+sqrt(1+B)./(1+2*U.*K_U.^2));

            x_new = sqrt(((1-l)/2).^2 + m./(y.^2)) - (1+l)/2;
            x_change=abs(x-x_new);
            x = x_new;
            loops=loops+1;
    end
    a = (gmu .* dt.^2)./(16*(rop.^2).*x.*(y.^2));

    sinBetaEOver2 = zeros(size(a));
    betaE = zeros(size(a));
    amin = zeros(size(a));
    tmin = zeros(size(a));
    alphaE = zeros(size(a));
    deltaE = zeros(size(a));
    alphaH = zeros(size(a));
    betaH = zeros(size(a));
    deltaH = zeros(size(a));
    f = zeros(size(a));
    g = zeros(size(a));
    g_dot = zeros(size(a));
    if(any(a(a > 0.0)))
        bool = a > 0.0;

        sinBetaEOver2(bool) = sqrt((s(bool)-c(bool))./(2*a(bool)));
        betaE(bool) = real(2*asin(sinBetaEOver2(bool))); %real needed to prevent complex numbers from appearing

        if(any(deltaTA(bool) > pi))
            bool2 = bool & deltaTA > pi;
            betaE(bool2) = -betaE(bool2);
        end

        amin(bool) = s(bool)./2;
        tmin(bool) = sqrt(amin(bool).^3./gmu(bool)) .* (pi - betaE(bool) + sin(betaE(bool)));

        alphaE(bool) = real(2*asin(sqrt(s(bool)./(2*a(bool))))); %real needed to prevent complex numbers from appearing

        if(any(dt > tmin))
            bool2 = bool & dt > tmin;
            alphaE(bool2) = 2*pi - alphaE(bool2);
        end

        deltaE(bool) = alphaE(bool) - betaE(bool);

        f(bool)     = real(1 - (a(bool)./r1Mag(bool)) .* (1 - cos(deltaE(bool)))); %real() necessary to prevent complex doubles from forming
        g(bool)     = real(dt(bool) - sqrt(a(bool).^3./gmu(bool)) .* (deltaE(bool) - sin(deltaE(bool)))); %real() necessary to prevent complex doubles from forming
        g_dot(bool) = real(1 - (a(bool)./r2Mag(bool)) .* (1- cos(deltaE(bool))));  %real() necessary to prevent complex doubles from forming
    end
    if(any(a(a < 0.0)))
        bool = a < 0.0;

        alphaH(bool) = 2*asinh(sqrt(s(bool)./(-2*a(bool))));
        betaH(bool)  = 2*asinh(sqrt((s(bool)-c(bool))./(-2*a(bool))));
        deltaH(bool) = alphaH(bool)-betaH(bool);

        f(bool)     = real(1 - (a(bool)./r1Mag(bool)) .* (1 - cosh(deltaH(bool)))); %real() necessary to prevent complex doubles from forming
        g(bool)     = real(dt(bool) - sqrt(-(a(bool).^3)./gmu(bool)) .* (sinh(deltaH(bool)) - deltaH(bool))); %real() necessary to prevent complex doubles from forming
        g_dot(bool) = real(1 - (a(bool)./r2Mag(bool)) .* (1 - cosh(deltaH(bool)))); %real() necessary to prevent complex doubles from forming
    end
    if(any(a(a == 0.0)))
        error('a = 0.0');
    end

    v1 = bsxfun(@rdivide,r2 - bsxfun(@times,r1,f),g);
    v2 = bsxfun(@rdivide,bsxfun(@times,r2,g_dot) - r1,g);
end

function ksi = computeKsi(x, numLevels)
    eta = x./(sqrt(1+x) + 1).^2;
    num = 8.*(sqrt(1+x) + 1);
    denom = 1;

    if(numLevels>0)
        denom = 3 + 1./(eta + computeKsi(eta, numLevels-1));
    end
    ksi = num ./ denom;
%     disp(ksi);
end

function [K_U]=Kay(U)
    % setup the C variable
    c(1)=4/27;   %@n=0
    c(2)=8/27; %@n=1
    c(3)=208/891; %n=2
    c(4)=340/1287; %@n=3
    c(5)=700/2907;
    c(6)=928/3591;
    c(7)=296/1215;
    c(8)=1804/7047;
    c(9)=2548/10395;
    c(10)=2968/11655;
    c(11)=3904/15867;
    c(12)=884/3483;
    c(13)=5548/22491;
    c(14)=6160/24327;
    c(15)=7480/30267;
    c(16)=8188/32391;
    c(17)=1940/7839;
    c(18)=10504/41607;
    c(19)=12208/49275;
    c(20)=13108/51975;

    % sum up all the intermediate variables
    Z=1+c(20).*U;
    Z=1+c(19).*U./Z;
    Z=1+c(18).*U./Z;
    Z=1+c(17).*U./Z;
    Z=1+c(16).*U./Z;
    Z=1+c(15).*U./Z;
    Z=1+c(14).*U./Z;
    Z=1+c(13).*U./Z;
    Z=1+c(12).*U./Z;
    Z=1+c(11).*U./Z;
    Z=1+c(10).*U./Z;
    Z=1+c(9).*U./Z;
    Z=1+c(8).*U./Z;
    Z=1+c(7).*U./Z;
    Z=1+c(6).*U./Z;
    Z=1+c(5).*U./Z;
    Z=1+c(4).*U./Z;
    Z=1+c(3).*U./Z;
    Z=1+c(2).*U./Z;
    Z=1+c(1).*U./Z;

    K_U=1./(3*Z);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [gp, diagnosis] = svigp(gp, x, y, varargin)
%SVIGP  Stochastic variational inference for GP
%
%  Description
%    GP = SVIGP(GP, X, Y, OPTIONS) optimises the variational, likelihood
%    and covariance function parameters of a sparse SVI GP model
%    given matrix X of training inputs and vector Y of training targets.
%    The model follows the description in Hensman, Fusi and Lawrence
%    (2013). Supported likelihood functions are gaussian (for regression)
%    and probit (for classification).
%
%    [GP, DIAGNOSIS] = SVIGP(GP, X, Y, OPTIONS) also returns values
%    monitored during the iteration. The structure DIAGNOSIS contains
%    the energy in the field e, the likelihood and covariance function
%    parameters in the field w, and the mean log predictive density (if xt
%    and yt is provided) in the field mlpd. The first dimension corresponds
%    to the main iteration in all these arrays. The second dimension
%    corresponds to the minibatch iteration in e and w. The third dimension
%    corresponds to the different parameters in w.
%
%    OPTIONS is optional parameter-value pair
%      xt               - test inputs. Used for monitoring the mean log
%                         predictive density, N.B. the monitoring needs
%                         both xt and yt.
%      yt               - observed yt in test points. Used for monitoring
%                         the mean log predictive density, N.B. the
%                         monitoring needs both xt and yt.
%      z                - optional observed quantity in triplet
%                         (x_i,y_i,z_i) Some likelihoods may use this. For
%                         example, in case of Poisson likelihood we have
%                         z_i=E_i, that is, expected value for ith case.
%      zt               - optional observed quantity in triplet
%                         (xt_i,yt_i,zt_i) Some likelihoods may use this.
%                         For example, in case of Poisson likelihood we
%                         have z_i=E_i, that is, the expected value for the
%                         ith case. N.B. used only in the monitoring of the
%                         mean log predictive density.
%      X_u              - inducing inputs. If omitted, kmeans clustering is
%                         applied and the resulting cluster centroids are
%                         selected. The number of inducing variables is
%                         then controlled by the parameter nu (see below).
%      nu               - the number of inducing variables if X_u is
%                         omitted. Can be given as an absolute value or
%                         relative to the input size. The default is
%                         min(floor(0.1*n),1500), where n is the number
%                         of training inputs.
%      m                - initial mean of the inducing variables. The
%                         default is zero vector.
%      S                - initial covariance of the inducing variables. The
%                         default is 0.1 times identity matrix.
%      n_minibatch      - absolute or relative size of the minibatches
%                         (relative to the number of the training inputs).
%                         Does not have to be divisible with the number of
%                         training inputs. The default is 0.1 (relative).
%      maxiter          - the maximum number of iterations (default 5000).
%                         Providing 0 initialises the gp structure
%                         parameters without optimisation.
%      momentum         - momentum term for the covariance function
%                         parameters (default 0.9)
%      mu1              - initial step size of the variational parameters
%                         (default 0.01).
%      mu2              - initial step size of the likelihood and
%                         covariance function parameters (default 1e-5)
%      tol              - tolerance of energy for the convergence
%                         (default 1e-6).
%      step_size        - function handle for the step size as a function
%                         of the iteration. The default function is f(i) =
%                         1/(1+i/n_minibatch), where n_minibatch is the
%                         size of the minibatches.
%      lik_sigma2       - likelihood variance in the case of non-gaussian
%                         likelihood (default 0.1).
%      lik_sigma2_prior - prior for the likelihood variance in the case of
%                         non-gaussian likelihood. The default is
%                         prior_loggaussian('mu',-2, 's2', 0.5).
%      display          - Control the amount of diagnostic verbosity.
%                         'off' displays nothing, 'final' display the
%                         final output, and 'iter' displays output at
%                         each iteration. Alternatively by providing a
%                         scalar value, fewer iterations can be displayed,
%                         e.g. value 10 displays only every tenth
%                         value (default behaviour).
%
%  See also
%    GP_SET, GPSVI_PRED, GPSVI_PREDGRAD, GPSVI_E, GPSVI_G, DEMO_SVI*
%
%  References:
%    Hensman, J., Fusi, N. and Lawrence, N. D. (2013). Gaussian processes
%    for big data. arXiv preprint arXiv:1309.6835.

% Copyright (c) 2014 Ville Tolvanen
% Copyright (c) 2014 Tuomas Sivula

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

ip=inputParser;
ip.FunctionName = 'SVIGP';
ip.addRequired('gp',@isstruct);
ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addParamValue('xt', [], @(x) isnumeric(x) && isreal(x) && all(isfinite(x(:))))
ip.addParamValue('yt', [], @(x) isnumeric(x) && isreal(x) && all(isfinite(x(:))))
ip.addParamValue('z', [], @(x) isnumeric(x) && isreal(x) && all(isfinite(x(:))))
ip.addParamValue('zt', [], @(x) isnumeric(x) && isreal(x) && all(isfinite(x(:))))
ip.addParamValue('X_u', [], @(x) isnumeric(x) && isreal(x) && all(isfinite(x(:))))
ip.addParamValue('nu', [], @(x) isreal(x) && isscalar(x) && x > 0)
ip.addParamValue('m', [], @(x) isnumeric(x) && isreal(x) && all(isfinite(x(:))))
ip.addParamValue('S', [], @(x) isnumeric(x) && isreal(x) && all(isfinite(x(:))))
ip.addParamValue('n_minibatch', 0.1, @(x) isreal(x) && isscalar(x) && x > 0)
ip.addParamValue('maxiter', 2000, @(x) isreal(x) && isscalar(x) && x >= 0)
ip.addParamValue('momentum', 0.9, @(x) isreal(x) && isscalar(x) && x > 0)
ip.addParamValue('mu1', 0.01, @(x) isreal(x) && isscalar(x) && x > 0)
ip.addParamValue('mu2', 1e-5, @(x) isreal(x) && isscalar(x) && x > 0)
ip.addParamValue('tol', 1e-6, @(x) isreal(x) && isscalar(x))
ip.addParamValue('step_size', [], @(x) isa(x,'function_handle'))
ip.addParamValue('lik_sigma2',0.1, @(x) isscalar(x) && x>=0);
ip.addParamValue('lik_sigma2_prior',prior_loggaussian('mu',-2, 's2', 0.5), ...
  @(x) isstruct(x) || isempty(x));
ip.addParamValue('display', 10, @(x) ismember(x,{'final','iter','off'}) ...
  || (isreal(x) && isscalar(x) && x > 1) )

ip.parse(gp, x,y,varargin{:});
x=ip.Results.x;
y=ip.Results.y;
z=ip.Results.z;
xt=ip.Results.xt;
yt=ip.Results.yt;
zt=ip.Results.zt;
X_u=ip.Results.X_u;
nu=ip.Results.nu;
m=ip.Results.m;
S=ip.Results.S;
n_minibatch=ip.Results.n_minibatch;
momentum=ip.Results.momentum;
mu1=ip.Results.mu1;
mu2=ip.Results.mu2;
step_size=ip.Results.step_size;
maxiter=ip.Results.maxiter;
tol=ip.Results.tol;
lik_sigma2 = ip.Results.lik_sigma2;
lik_sigma2_prior = ip.Results.lik_sigma2_prior;
display = ip.Results.display;

% Initialise the diagnosis output structure
diagnosis = struct();

% Check if latent method SVI has been set
if ~isfield(gp, 'latent_method') || ~isequal(gp.latent_method, 'SVI')
  gp=gp_set(gp, 'latent_method', 'SVI');
end

% Check if latent method is not gaussian or probit
if ~strcmp(gp.lik.type, 'Gaussian') && ~strcmp(gp.lik.type, 'Probit')
  error('Supported likelihoods for SVIGP are gaussian and probit.')
end

% Process parameters
if xor(isempty(xt), isempty(yt))
  warning('Need both xt and yt for monitoring mean log predictive density.');
end
n=size(x,1);

if n_minibatch < 1
  n_minibatch = max(floor(n_minibatch*n),1);
end
if n_minibatch > n
  n_minibatch = max(floor(0.1*n),1);
  warning('Too many minibatches, using floor(0.1*n) = %d instead.', ...
    n_minibatch)
end

if isempty(step_size)
  step_size=@(iter) 1/(1+iter./n_minibatch);
end

% Handle the inducing inputs
if ~ismember('X_u',ip.UsingDefaults)
  gp.X_u = X_u;
end
if isempty(gp.X_u)
  % Assign X_u by clustering
  if isempty(nu)
    nu = min(max(floor(0.1.*n),1),1500);
  elseif nu < 1
    nu = max(floor(nu.*n),1);
  end
  fprintf('Assign inducing inputs by clustering\n')
  Sw=warning('off','stats:kmeans:EmptyCluster');
  [~,X_u] = kmeans(x, nu,'Start','uniform',...
    'EmptyAction','singleton');
  warning(Sw);
  gp.X_u = X_u;
end
gp.nind = size(gp.X_u,1);

% Handle the rest of the parameters
if ~ismember('m',ip.UsingDefaults)
  if length(m) == gp.nind
    gp.m = m;
  else
    error('The size of m does not match with X_u')
  end
elseif ~isfield(gp, 'm') || length(gp.m) ~= gp.nind
  gp.m = zeros(gp.nind,1);
end

if ~ismember('S',ip.UsingDefaults)
  if ismatrix(S) && all(size(S) == gp.nind)
    gp.S = S;
  else
    error('The size of S does not match with X_u')
  end
elseif ~isfield(gp, 'S') || any(size(gp.S) ~= gp.nind)
  gp.S = 0.1*eye(gp.nind);
  % gp.S = gp_trcov(gp,gp.X_u);
end
gp.t1=gp.S\gp.m;
gp.t2=-0.5.*inv(gp.S);

% Handle the likelihood variance
if ~isfield(gp.lik, 'sigma2')
  gp.lik.sigma2 = lik_sigma2;
  gp.lik.p.sigma2 = lik_sigma2_prior;
  gp.lik.fh.lp = @lik_lp;
  gp.lik.fh.lpg = @lik_lpg;
end

% Return prematurely if only initialising
if maxiter == 0
  return
end

% Parameters
w = gp_pak(gp);
w0 =w;
nh1 = numel(gp.t1) + numel(gp.t2);
nh2 = length(w)-nh1;

% Initial step-size vector
mu0 = mu1.*ones(size(w));
mu0(end-nh2+1:end)=mu2;

% Size of minibatches
nb=n_minibatch;
% Number of minibatches
nbb=ceil(n/nb);

% Monitored values
if nargout > 1
  e_all = zeros(maxiter,nbb);
  w_all = zeros(maxiter,nbb,nh2);
  if ~isempty(xt) && ~isempty(yt)
    mlpd_all = zeros(maxiter,1);
    rmse_all = zeros(maxiter,1);
  end
end
g_old=zeros(size(w));
momentum=momentum.*ones(size(w));
etot_old=Inf;

% Preprocess conditions for iteration
disp_iter = strcmp(display, 'iter');
disp_i = 0;
disp_count = display;

% Iterate until convergence or maxiter
converged = 0;
try_fix_mu2 = 0;
for iter = 1:maxiter

  % Divide the data into minibatches
  inds = cell(nbb,1);
  ind = randperm(n);
  for i=1:nbb-1
    inds{i} = ind((i-1)*nb+1:i*nb);
  end
  inds{nbb} = ind((nbb-1)*nb+1:end);

  % Compute the step-size
  mu = mu0.*step_size(iter);

  % Iterate all the minibatches
  etot = 0;
  broken = 0;
  for i=1:nbb
    if isempty(z)
      zi = [];
    else
      zi = z(inds{i},:);
    end
    gp.data_prop=length(inds{i})./n;
    [e,~,~,param] = gpsvi_e(w,gp,x(inds{i},:),y(inds{i},:), 'z', zi);
    etot = etot+e;
    g = gpsvi_g(w,gp,x(inds{i},:),y(inds{i},:), 'z', zi, ...
      'gpsvi_e_param', param);
    g = mu.*g + momentum.*g_old;
    g_old = g;
    w = w+g;
    if nargout > 1
      e_all(iter,i) = e;
      w_all(iter,i,:) = w(end-nh2+1:end);
    end
    if ~isnan(etot) ...
        && all(~isinf(exp(w(end-nh2+1:end)))) ...
        && all(exp(w(end-nh2+1:end))~=0) ...
        && ~any(isnan(g)) ...
        && ~isnan(gpsvi_e(w+g,gp,x(inds{i},:),y(inds{i},:), 'z', zi))
      gp = gp_unpak(gp,w);
    elseif ~try_fix_mu2
      fprintf('Bad parameter values, decreasing mu2.\n');
      try_fix_mu2 = 1;
      mu2 = 0.1*mu2;
      mu0(end-nh2+1:end)=mu2;
      w = w0;
      g_old = zeros(1,nh1+nh2);
      broken = 1;
      break
    else
      fprintf('Bad parameter values, decreasing step-size and momentum.\n');
      mu0 = 0.1.*mu0;
      momentum = 0.5*momentum;
      w = w0;
      g_old = zeros(1,nh1+nh2);
      broken = 1;
      break
    end
    gpsvi_e('clearcache',gp);
  end
  if broken
    continue
  end
  etot=etot/nbb;

  % Analyse and print
  if isscalar(display)
    if disp_count == display
      disp_i = 1;
      disp_count = 1;
    else
      disp_i = 0;
      disp_count = disp_count +1;
    end
  end
  if ~isempty(xt) && ~isempty(yt) ...
      && ( disp_iter || nargout > 1 || disp_i)
    [Eft,~,lpyt] = gpsvi_pred(gp,x,y,xt,'yt',yt, 'z', zi, 'zt', zt);
    lpyt = mean(lpyt);
    if nargout > 1
      mlpd_all(iter) = lpyt;
      rmse = sqrt(mean((yt-Eft).^2));
      rmse_all(iter) = rmse;
    end
    if disp_iter || disp_i
      fprintf('iter=%d/%d, e=%.3f, mlpd=%.4g, rmse=%.4g, de=%.5g\n', ...
        iter, maxiter, etot, lpyt, rmse, abs(etot-etot_old));
    end
  elseif disp_iter || disp_i
    fprintf('iter=%d/%d, e=%.3f, de=%.5g\n', ...
      iter, maxiter, etot, abs(etot-etot_old));
  end

  % Check for convergence
  if abs(etot-etot_old)<tol
    converged = 1;
    if strcmp(display, 'iter')
      fprintf('Energy converged\n')
    end
    break
  end
  etot_old=etot;

end

% Display results
if ~converged
  fprintf('Iteration limit %d reached while optimising the parameters\n', ...
    maxiter)
elseif ~strcmp(display, 'off')
  fprintf('Tolerance reached in %d iterations\n', iter)
end
if strcmp(display, 'final') || isscalar(display)
  if ~isempty(xt) && ~isempty(yt) && ( nargout > 1 || disp_i)
    fprintf(['Final values:\n' ...
      'e=%.3f, mlpd=%.4g, rmse=%.4g\n'], etot, lpyt, rmse);
  elseif ~isempty(xt) && ~isempty(yt)
    [Eft,~,lpyt] = gpsvi_pred(gp,x,y,xt,'yt',yt, 'z', z, 'zt', zt);
    fprintf(['Final values:\n' ...
      'e=%.3f, mlpd=%.4g, rmse=%.4g\n'], ...
      etot, mean(lpyt)), sqrt(mean((yt-Eft).^2));
  else
    fprintf('Final energy: %.3f\n', etot);
  end
end

% Save the monitored values
if nargout > 1
  diagnosis.e = e_all(1:iter,:);
  diagnosis.w = w_all(1:iter,:,:);
  if ~isempty(xt) && ~isempty(yt)
    diagnosis.mlpd = mlpd_all(1:iter);
    diagnosis.rmse = rmse_all(1:iter);
  end
end

end

function lp = lik_lp(lik, varargin)
%LIK_LP  log(prior) of the likelihood parameters
%
%  Description
%    LP = LIK_LP(LIK) takes a likelihood structure LIK and
%    returns log(p(th)), where th collects the parameters. This
%    subfunction is needed if there are likelihood parameters.
%    Added for non-gaussian likelihoods in SVIGP.
%
%  See also
%    LIK_*, SVIGP

% If prior for sigma2sion parameter, add its contribution
lp=0;
if ~isempty(lik.p.sigma2)
  lp = lik.p.sigma2.fh.lp(lik.sigma2, lik.p.sigma2) +log(lik.sigma2);
end
end

function lpg = lik_lpg(lik)
%LIK_LPG  d log(prior)/dth of the likelihood parameters
%
%  Description
%    E = LIK_NEGBIN_LPG(LIK) takes a likelihood structure LIK and
%    returns d log(p(th))/dth, where th collects the parameters.
%    This subfunction is needed if there are likelihood parameters.
%    Added for non-gaussian likelihoods in SVIGP.
%
%  See also
%    LIK_*, SVIGP

lpg=[];
if ~isempty(lik.p.sigma2)
  % Evaluate the gprior with respect to sigma2
  ggs = lik.p.sigma2.fh.lpg(lik.sigma2, lik.p.sigma2);
  lpg = ggs(1).*lik.sigma2 + 1;
  if length(ggs) > 1
    lpg = [lpg ggs(2:end)];
  end
end
end


text:
function [LL, LLS, ht] = tarch_likelihood(parameters, data, fdata, fIdata, p, o, q, error_type, tarch_type, back_cast, T, estim_flag)
% Log likelihood for TARCH(P,O,Q) estimation
%
% USAGE:
%   [LL, LLS, HT] = tarch_likelihood(PARAMETERS, DATA, FDATA, FIDATA, P, O, Q, ERROR_TYPE, TARCH_TYPE, BACK_CAST, T, ESTIM_FLAG)
%
% INPUTS:
%   PARAMETERS    - A vector of GARCH process parameters
%                   [omega alpha gamma beta [nu lambda]]
%   DATA          - Vector of mean zero residuals
%   FDATA         - Either abs(data) or data.^2, depending on tarch_type
%   FIDATA        - fdata times an indicator for negative, e.g. fdata.*(data<0)
%   P             - The lag order length for ARCH
%   O             - The lag order of asymmetric terms
%   Q             - The lag order length for GARCH
%   ERROR_TYPE    - The type of error being assumed, valid types are:
%                     1 if 'NORMAL'
%                     2 if 'STUDENTST'
%                     3 if 'GED'
%                     4 if 'SKEWT'
%   TARCH_TYPE    - 1 for absolute vale of return
%                 - 2 for squared returns (standard case)
%   BACK_CAST     - The value used for variance recursion
%   T             - Length of data
%   ESTIM_FLAG    - [OPTIONAL] Flag (0 or 1) to indicate if the function
%                   is being used in estimation.  If it is 1, then the parameters are
%                   transformed from unconstrained values to constrained by standard
%                   garch model constraints
%
% OUTPUTS:
%   LL             - Minus 1 times the log likelihood
%   LLS            - Time series of log likelihoods (Also multiplied by -1)
%   HT             - Time series of conditional variances
%
% COMMENTS:
%   See also TARCH

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 9/1/2005

if nargin==12 && estim_flag
    %If for estimation, transform the parameters
    [parameters,nu,lambda]=tarch_itransform(parameters,p,o,q,error_type);
else
    %Otherwise the parameters simply must be parsed
    if error_type==2 || error_type==3
        %Seperate nu from the remaning parameters
        nu=parameters(p+o+q+2);
        parameters=parameters(1:1+p+o+q);
    elseif error_type==4
        lambda=parameters(p+o+q+3);
        nu=parameters(p+o+q+2);
        parameters=parameters(1:1+p+o+q);
    end
end

%Backcast length
m  =  max([p o q]);
%Compute the conditional variances
ht=tarch_core(fdata,fIdata,parameters,back_cast,p,o,q,m,T,tarch_type);

%Indices for the relevant opservations
t = (m + 1):T;
ht=ht(t);
data=data(t);
%Compute the log likelihoods
switch error_type
    case 1
        [LL, LLS] = normloglik(data,0,ht);
        LLS = -LLS;
        LL = -LL;
    case 2
        [LL, LLS] = stdtloglik(data,0,ht,nu);
        LLS = -LLS;
        LL = -LL;
    case 3
        [LL, LLS] = gedloglik(data,0,ht,nu);
        LLS = -LLS;
        LL = -LL;
    case 4
        [LL, LLS] = skewtloglik(data,0,ht,nu,lambda);
        LLS = -LLS;
        LL = -LL;
end

author: bashtage

repo: mfe-toolbox


text:
function errplotl(sol,eldata,ev,xy,x,y,fig)
%errplotl   plots solution and error estimate on L-shaped domain
%   errplotl(sol,eldata,ev,xy,x,y,fig)
%   input
%          sol        nodal solution vector
%          eldata     element error vector
%          ev         element mapping matrix
%          xy         vertex coordinate vector
%          x          vector of x-axis interpolation points
%          y          vector of y-axis interpolation points
%          fig        figure number
%
%   IFISS function: DJS; 5 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
fprintf('plotting solution and estimated errors... ')
% interpolate to a cartesian product mesh
[X,Y]=meshgrid(x,y);
xysol = griddata(xy(:,1),xy(:,2),sol,X,Y);
[II,JJ]=find(X<0 & Y<0); xysol(II,JJ)=nan;
figure(fig)
subplot(221),contour(X,Y,xysol,20),axis('square')
axis('off'), ellx, title('Finite Element Solution')
subplot(222),mesh(X,Y,xysol),axis('square')
view(330,30)
%%
xx=xy(:,1); yy=xy(:,2);
nel=length(eldata);
% loop over elements
for ielem = 1:nel
xl = xx(ev(ielem,:));
yl = yy(ev(ielem,:));
xc(ielem,1) = 0.25*sum(xl);
xc(ielem,2) = 0.25*sum(yl);
end
%
% interpolate to a cartesian product mesh
x=0.5*(x(1:end-1)+x(2:end));
y=0.5*(y(1:end-1)+y(2:end));
[X,Y]=meshgrid(x,y);
xysol = griddata(xc(:,1),xc(:,2),eldata,X,Y);
[II,JJ]=find(X<0 & Y<0); xysol(II,JJ)=nan;
subplot(223),contour(X,Y,xysol,15),axis('square')
axis('off'), ellx, title('Estimated Error')
subplot(224),mesh(X,Y,xysol),axis('square')
view(330,30)
subplot(111)
fprintf('done\n')
return


text:
function [fval,optimState,fsd] = funlogger(fun,u,optimState,state,varargin)
%FUNLOGGER Call objective function and do some bookkeeping.
%   [~,OPTIMSTATE] = FUNLOGGER(FUN,U,OPTIMSTATE,'init') starts logging
%   function FUN with starting point U and optimization struct OPTIMSTATE.
%
%   [~,OPTIMSTATE] = FUNLOGGER(FUN,U,OPTIMSTATE,'init',NMAX) stores up to
%   NMAX function values (default NMAX=1e4).
%
%   [~,OPTIMSTATE] = FUNLOGGER(FUN,U,OPTIMSTATE,'init',NMAX,1) also stores
%   heteroskedastic noise (second output argument) from the logged function.
%
%   [FVAL,OPTIMSTATE] = FUNLOGGER(FUN,U,OPTIMSTATE,'iter') evaluates function
%   FUN at U with optimization struct OPTIMSTATE, return function value FVAL.
%   FUN must take a vector input and returns a scalar value and, optionally,
%   the (estimated) SD of the returned value (if heteroskedastic noise
%   handling is on).
%
%   [FVAL,OPTIMSTATE] = FUNLOGGER(FUN,U,OPTIMSTATE,'single') as 'iter' but
%   does not store function values in the cache.
%
%   [~,OPTIMSTATE] = FUNLOGGER(FUN,U,OPTIMSTATE,'done') finalizes stored
%   function values.

%   Luigi Acerbi 2017

fval = [];
fsd = [];

switch lower(state)
    case 'init' % Start new function logging session

        nvars = numel(u);

        % Number of stored function values
        if nargin > 4; nmax = varargin{1}; else; nmax = []; end
        if isempty(nmax); nmax = 1e4; end

        % Heteroscedastic noise provided as 2nd output argument of FUN?
        if nargin > 5; henoise_flag = varargin{2}; else; henoise_flag = []; end
        if isempty(henoise_flag); henoise_flag = false; end

        optimState.funccount = 0;

        % Passing empty or new OPTIMSTATE
        if ~isfield(optimState,'X') || isempty(optimState.X)
            optimState.X = NaN(nmax,nvars);
            optimState.U = NaN(nmax,nvars);
            optimState.Y = NaN(nmax,1);
            if henoise_flag; optimState.S = NaN(nmax,1); end
            optimState.Xn = 0;                    % Last filled entry
            optimState.Xmax = 0;                  % Maximum entry index

        else % Receiving previous evaluations (e.g. from previous run)

            % Find last evaluated point
            optimState.Xmax = find(~isnan(optimState.Y),1,'last');
            if ~isfield(optimState,'Xn') || isempty(optimState.Xn)
                optimState.Xn = optimState.Xmax;
            end

            % Current cache is smaller than previous one, keep most recent
            % evaluated points
            if optimState.Xmax > nmax
                optimState.X = circshift(optimState.X,-optimState.Xn);
                optimState.X = optimState.X(end-nmax+1:end,:);
                optimState.Y = circshift(optimState.Y,-optimState.Xn);
                optimState.Y = optimState.Y(end-nmax+1:nmax,:);
                if henoise_flag
                    optimState.S = circshift(optimState.S,-optimState.Xn);
                    optimState.S = optimState.S(end-nmax+1:nmax,:);
                end
                optimState.Xn = 0;
                optimState.Xmax = nmax;
            else
                offset = nmax - size(optimState.X,1);
                optimState.X = [optimState.X; NaN(offset,nvars)];
                optimState.Y = [optimState.Y; NaN(offset,1)];
                if henoise_flag
                    optimState.S = [optimState.S; NaN(offset,1)];
                end
            end
            optimState.U = gridunits(optimState.X,optimState);
        end
        optimState.funevaltime = NaN(nmax,1);
        optimState.totalfunevaltime = 0;

    case {'iter','single'} % Evaluate function (and store output for 'iter')

        x = origunits(u,optimState);    % Convert back to original space
        % Heteroscedastic noise?
        henoise_flag = isfield(optimState,'S');

        try
            t_fun = tic;
            if henoise_flag
                [fval,fsd] = fun(x);
            else
                fval = fun(x);
            end
            t = toc(t_fun);

            % Check returned function value
            if ~isscalar(fval) || ~isfinite(fval) || ~isreal(fval)
                error(['The returned function value must be a finite real-valued scalar (returned value: ' mat2str(fval) ').']);
            end

            % Check returned function SD
            if henoise_flag && (~isscalar(fsd) || ~isfinite(fsd) || ~isreal(fsd) || fsd <= 0.0)
                error(['The returned estimated SD (second function output) must be a finite, non-negative real-valued scalar (returned SD: ' mat2str(fsd) ').']);
            end

        catch fun_error
            warning('bads:funError', ['Error in executing the logged function ''' func2str(fun) ''' with input: ' mat2str(x)]);
            rethrow(fun_error);
        end

        % Update function records
        optimState.funccount = optimState.funccount + 1;
        if strcmpi(state, 'iter')
            optimState.Xn = max(1,mod(optimState.Xn+1, size(optimState.X,1)));
            optimState.Xmax = min(optimState.Xmax+1, size(optimState.X,1));
            optimState.X(optimState.Xn,:) = x;
            optimState.U(optimState.Xn,:) = u;
            optimState.Y(optimState.Xn) = fval;
            if henoise_flag
                optimState.S(optimState.Xn) = fsd;
            end
            optimState.funevaltime(optimState.Xn) = t;
        end
        optimState.totalfunevaltime = optimState.totalfunevaltime + t;

    case 'done' % Finalize stored table

        henoise_flag = isfield(optimState,'S');

        if optimState.Xmax < size(optimState.X,1)
            optimState.X = optimState.X(1:optimState.Xmax,:);
            optimState.Y = optimState.Y(1:optimState.Xmax);
            if henoise_flag; optimState.S = optimState.S(1:optimState.Xmax); end
            optimState.funevaltime = optimState.funevaltime(1:optimState.Xmax);
        else
            optimState.X = circshift(optimState.X,-optimState.Xn);
            optimState.Y = circshift(optimState.Y,-optimState.Xn);
            if henoise_flag; optimState.S = circshift(optimState.S,-optimState.Xn); end
            optimState.funevaltime = circshift(optimState.funevaltime,-optimState.Xn);
        end
        optimState = rmfield(optimState,'U');

    otherwise
        error('Unknown FUNLOGGER action.');
end

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function [ res ] = pesq3( reference_sig, degraded_sig, Fs, fileNum )
% A wrapper for the objective Perceptual Evaluation of Speech Quality measure
%
% Syntax:	[ res ] = pesq3( reference_sig, degraded_sig, Fs, fileNum )
%
% Inputs:
% 	reference_sig - Reference (clean, talker, sender) speech signal
% 	degraded_sig - Degraded (noisy, listener, receiver) speech signal
% 	Fs - Sampling Frequency
%   fileNum - An ID number to append to the temporary audio files. Useful
%   when several instances are to be run together (in parallel).
%
% Outputs:
% 	res - Raw PESQ result for narrowband and MOS-LQO result for wideband
%
% See also: pesq2mos.m

% Author: Jacob Donley
% University of Wollongong
% Email: jrd089@uowmail.edu.au
% Copyright: Jacob Donley 2017
% Date: 03 October 2015
% Revision: 0.1
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 4
    fileNum = 0;
end
temp_path = [pwd filesep '+Miscellaneous\+Temporary\'];
ref_path = [ 'tmp_ref' num2str(fileNum) '.wav'];
deg_path = [ 'tmp_deg' num2str(fileNum) '.wav'];

if ~exist(temp_path,'dir'); mkdir(temp_path); end

max_val = max(abs([reference_sig(:); degraded_sig(:)]));

audiowrite([temp_path ref_path], reference_sig / max_val, Fs);
audiowrite([temp_path deg_path], degraded_sig / max_val, Fs);

res = pesq2_mtlb(ref_path, ...
                 deg_path, ...
                 Fs, 'wb', [pwd filesep '+Tools\pesq_NoResFile.exe'], ...
                 temp_path);

end

author: jdonley


text:
function phiFaceAverage = tvdMean3D(phi, u, FL)
% This function gets the value of the field variable phi defined
% over the MeshStructure and calculates the upwind average on
% the cell faces, based on the direction of the velocity vector for a uniform mesh.
%
% SYNOPSIS:
%   phiFaceAverage = tvdMean3D(u, phi, FL)
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Written by Ali A. Eftekhari
% See the license file

% extract the velocity data
% note: size(ux) = [1:m+1, 1:n] and size(uy) = [1:m, 1:n+1]
ux = u.xvalue;
uy = u.yvalue;
uz = u.zvalue;

% check the size of the variable and the mesh dimension
Nx = u.domain.dims(1);
Ny = u.domain.dims(2);
Nz = u.domain.dims(3);
dx=repmat(0.5*(u.domain.cellsize.x(1:end-1)+u.domain.cellsize.x(2:end)), 1, Ny, Nz);
dy=repmat(0.5*(u.domain.cellsize.y(1:end-1)+u.domain.cellsize.y(2:end))', Nx, 1, Nz);
dz=zeros(1, 1, Nz+1);
dz(1,1,:)=0.5*(u.domain.cellsize.z(1:end-1)+u.domain.cellsize.z(2:end));
dz=repmat(dz, Nx, Ny, 1);

% define the tvd face vectors
phiX_p = zeros(Nx+1,Ny,Nz);
phiX_m = zeros(Nx+1,Ny,Nz);
phiY_p = zeros(Nx,Ny+1,Nz);
phiY_m = zeros(Nx,Ny+1,Nz);
phiZ_p = zeros(Nx,Ny,Nz+1);
phiZ_m = zeros(Nx,Ny,Nz+1);

% calculate the upstream to downstream gradient ratios for u>0 (+ ratio)
% x direction
dphiX_p = (phi.value(2:Nx+2, 2:Ny+1, 2:Nz+1)-phi.value(1:Nx+1, 2:Ny+1, 2:Nz+1))./dx;
rX_p = dphiX_p(1:end-1,:,:)./fsign(dphiX_p(2:end,:,:));
phiX_p(2:Nx+1,:,:) = phi.value(2:Nx+1, 2:Ny+1, 2:Nz+1)+0.5*FL(rX_p).* ...
    (phi.value(3:Nx+2,2:Ny+1,2:Nz+1)-phi.value(2:Nx+1,2:Ny+1,2:Nz+1));
phiX_p(1,:,:) = (phi.value(1,2:Ny+1,2:Nz+1)+phi.value(2,2:Ny+1,2:Nz+1))/2; % left boundary
% y direction
dphiY_p = (phi.value(2:Nx+1, 2:Ny+2, 2:Nz+1)-phi.value(2:Nx+1, 1:Ny+1, 2:Nz+1))./dy;
rY_p = dphiY_p(:,1:end-1,:)./fsign(dphiY_p(:,2:end,:));
phiY_p(:,2:Ny+1,:) = phi.value(2:Nx+1, 2:Ny+1, 2:Nz+1)+0.5*FL(rY_p).* ...
    (phi.value(2:Nx+1,3:Ny+2,2:Nz+1)-phi.value(2:Nx+1, 2:Ny+1,2:Nz+1));
phiY_p(:,1,:) = (phi.value(2:Nx+1,1,2:Nz+1)+phi.value(2:Nx+1,2,2:Nz+1))/2; % Bottom boundary
% z direction
dphiZ_p = (phi.value(2:Nx+1, 2:Ny+1, 2:Nz+2)-phi.value(2:Nx+1, 2:Ny+1, 1:Nz+1))./dz;
rZ_p = dphiZ_p(:,:,1:end-1)./fsign(dphiZ_p(:,:,2:end));
phiZ_p(:,:,2:Nz+1) = phi.value(2:Nx+1, 2:Ny+1, 2:Nz+1)+0.5*FL(rZ_p).* ...
    (phi.value(2:Nx+1,2:Ny+1,3:Nz+2)-phi.value(2:Nx+1,2:Ny+1,2:Nz+1));
phiZ_p(:,:,1) = (phi.value(2:Nx+1,2:Ny+1,1)+phi.value(2:Nx+1,2:Ny+1,2))/2; % Back boundary

% calculate the upstream to downstream gradient ratios for u<0 (- ratio)
% x direction
rX_m = dphiX_p(2:end,:,:)./fsign(dphiX_p(1:end-1,:,:));
phiX_m(1:Nx,:,:) = phi.value(2:Nx+1, 2:Ny+1, 2:Nz+1)+0.5*FL(rX_m).* ...
    (phi.value(1:Nx, 2:Ny+1, 2:Nz+1)-phi.value(2:Nx+1, 2:Ny+1, 2:Nz+1));
phiX_m(Nx+1,:,:) = (phi.value(end,2:Ny+1,2:Nz+1)+phi.value(end-1,2:Ny+1,2:Nz+1))/2; % right boundary
% y direction
rY_m = dphiY_p(:,2:end,:)./fsign(dphiY_p(:,1:end-1,:));
phiY_m(:,1:Ny,:) = phi.value(2:Nx+1,2:Ny+1,2:Nz+1)+0.5*FL(rY_m).* ...
    (phi.value(2:Nx+1,1:Ny,2:Nz+1)-phi.value(2:Nx+1,2:Ny+1,2:Nz+1));
phiY_m(:,Ny+1,:) = (phi.value(2:Nx+1, end,2:Nz+1)+phi.value(2:Nx+1, end-1,2:Nz+1))/2; % top boundary
% z direction
rZ_m = dphiZ_p(:,:,2:end)./fsign(dphiZ_p(:,:,1:end-1));
phiZ_m(:,:,1:Nz) = phi.value(2:Nx+1,2:Ny+1,2:Nz+1)+0.5*FL(rZ_m).* ...
    (phi.value(2:Nx+1,2:Ny+1,1:Nz)-phi.value(2:Nx+1,2:Ny+1,2:Nz+1));
phiZ_m(:,:,Nz+1) = (phi.value(2:Nx+1,2:Ny+1,end)+phi.value(2:Nx+1,2:Ny+1,end-1))/2; % front boundary

% calculate the average value
xvalue = (ux>0).*phiX_p+ ...
                        (ux<0).*phiX_m+ ...
                        0.5*(ux==0).*(phi.value(1:Nx+1,2:Ny+1,2:Nz+1)+phi.value(2:Nx+2,2:Ny+1,2:Nz+1));
yvalue = (uy>0).*phiY_p+ ...
                        (uy<0).*phiY_m+ ...
                        0.5*(uy==0).*(phi.value(2:Nx+1,1:Ny+1,2:Nz+1)+phi.value(2:Nx+1,2:Ny+2,2:Nz+1));
zvalue = (uz>0).*phiZ_p+ ...
                        (uz<0).*phiZ_m+ ...
                        0.5*(uz==0).*(phi.value(2:Nx+1,2:Ny+1,1:Nz+1)+phi.value(2:Nx+1,2:Ny+1,2:Nz+2));
phiFaceAverage=FaceVariable(phi.domain, xvalue, yvalue, zvalue);
end

function phi_out = fsign(phi_in)
% This function checks the value of phi_in and assigns an eps value to the
% elements that are less than or equal to zero, while keeping the signs of
% the nonzero elements
    phi_out = (abs(phi_in)>=eps).*phi_in+eps*(phi_in==0)+eps*(abs(phi_in)<eps).*sign(phi_in);
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function tagS = patient_module_tags_subset
%"patient_module_tags_subset"
%   Return a subset of the tags specified in C.7.1.1 in PS3.3 of the
%   2006 DICOM specification required for distinguishing patients.
%
% Tags are returned in a struct array with 3 fields:
%   Tag: String containing hex DICOM tag of a field.
%  Type: String describing type of field, with 5 options:
%         '1' Field must exist, data must exist and be valid.
%         '2' Field must exist, data can exist and be valid, or be NULL.
%         '3' Field is optional, if the field exists data can exist and be
%             valid or be NULL.
%         '1C' Field must exist under certain conditions and contain valid
%             data.
%         '2C' Field must exist under certain conditions and can contain
%             valid data or be NULL.
%Children: For sequences, a tagS with the same format as this struct array.
%
% AI 11/20/2020
%
%Usage:
%   tagS = patient_module_tags_subset
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

% Tags: {Patient's Name, PatientID}

tagC = {'00100010','00100020'};

tagDecC = {1048592,1048608};

tagTypeC = {'2','2'};

childC = {[],[]};

tagS = struct('tag', tagC, 'tagdec', tagDecC, 'type', tagTypeC, 'children', childC);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function success=testJabCreationWithScoreFeatureInAlternateOrder()

testDataDirName='/groups/branson/bransonlab/projects/JAABA/test_data';
%testDataDirName='/Users/taylora/jaaba/test_data';
nameOfExpToAdd='GMR_71G01_AE_01_TrpA_Rig2Plate14BowlD_20110707T154929';

% % Copy the fooing scores into the exp dir
% cmd=sprintf('cp "%s/scores_fooing_GMR_71G01_AE_01_TrpA_Rig1Plate15BowlA_20120316T144027.mat" "%s/%s/scores_fooing.mat"', ...
%             testDataDirName, ...
%             testDataDirName, ...
%             nameOfExpToAdd);
% system(cmd);

% Make sure the fooing score are _not_ already in the perframe dir
scorePerframeFileNameAbs= ...
  sprintf('%s/%s/perframe/scores_fooing.mat', ...
          testDataDirName, ...
          nameOfExpToAdd);
if exist(scorePerframeFileNameAbs,'file') ,
  cmd=sprintf('rm -f "%s"',scorePerframeFileNameAbs);
  system(cmd);
end

% Make a Macguffin, from which we'll open a new file in JLabelData
macguffin=Macguffin('flies');
macguffin.setMainBehaviorName('metafooing');
macguffin.setScoreFileName('scores_metafooing.mat');
macguffin.setTrxFileName('registered_trx.mat');
macguffin.setMovieFileName('movie.ufmf');

% Create the JLabelData object
data=JLabelData('setstatusfn',@(str)(fprintf('%s\n',str)), ...
                'clearstatusfn',@()(nop()), ...
                'isInteractive',false);

% Create a new file in the JLabelData object
data.newJabFile(macguffin);

% Add a new exp dir
nameOfExpDirToAdd='/groups/branson/bransonlab/projects/JAABA/test_data/GMR_71G01_AE_01_TrpA_Rig2Plate14BowlD_20110707T154929';
%nameOfExpDirToAdd='/Users/taylora/jaaba/test_data/GMR_71G01_AE_01_TrpA_Rig2Plate14BowlD_20110707T154929';
data.SetGenerateMissingFiles();  % tell JLabelData to generate any missing perframe files
[success,msg]=data.AddExpDir(nameOfExpDirToAdd);
if ~success ,
  error(msg);
end

% Add the fooing score feature
scoreFeatureJabFileNameAbs=fullfile(testDataDirName,'fooing_labels_and_classifier.jab');
fooingJLabelData=JLabelData('setstatusfn',@(~)(nop()), ...
                            'clearstatusfn',@()(nop()));
fooingJLabelData.openJabFile(scoreFeatureJabFileNameAbs,false);
timeStamp=fooingJLabelData.classifierTS;
scoreFeatureName=baseNameFromFileName(fooingJLabelData.scorefilename);
fooingJLabelData.closeJabFile();
fooingJLabelData=[];  %#ok
data.setScoreFeatures({scoreFeatureJabFileNameAbs},timeStamp,{scoreFeatureName});

% Check the status table
perframeDirExists=data.FileExists('perframedir');
if perframeDirExists ,
  fprintf('The JLabelData object says that the perframe directory exists, as it should.\n');
else
  error('The JLabelData object says that the perframe directory does not exist, which is not as it should be.');
end

% Add the score feature to the current vocabulary
featureLexicon = data.featureLexicon;
scoreFeatures = data.scoreFeatures;
toBeCalculatedPFNames=data.allperframefns;
windowFeatureParams = data.GetPerframeParams();
fv= ...
  FeatureVocabularyForSelectFeatures(featureLexicon, ...
                                     scoreFeatures, ...
                                     toBeCalculatedPFNames, ...
                                     windowFeatureParams);
fv.setPFToWFAmount('scores_fooing','normal');
fv.setPFEnablement('scores_fooing',true);
windowFeaturesParams=fv.getInJLabelDataFormat();
data.setWindowFeaturesParams(windowFeaturesParams);

% Change the current target to the first fly of the new experiment
iExp=whichstr(nameOfExpDirToAdd,data.expdirs);
if isempty(iExp)
  error('Added experiment seems to be missing.');
elseif length(iExp)>1
  error('Added experiment seems to be present more than once in the JLabelData object.');
end
iFly=1;
data.setCurrentTarget(iExp,iFly);

% Add a label
ts1=5000+(0:50);  % frame indices to set
iBehavior1=whichstri('metafooing',data.labelnames);
isImportant1=true;
data.SetLabel(iExp,iFly,ts1,iBehavior1,isImportant1);

% Add a label
ts2=5100+(0:50);  % frame indices to set
iBehavior2=whichstri('none',data.labelnames);
isImportant2=true;
data.SetLabel(iExp,iFly,ts2,iBehavior2,isImportant2);

% Train
data.Train();

% The predictions should be correct for the labeled bouts
data.Predict(iExp,iFly,ts1(1),ts1(end));  % Make sure the prediction has been done
pred1=data.GetPredictedIdx(iExp,iFly,ts1(1),ts1(end));
if all(pred1.predictedidx==iBehavior1)
  fprintf('Predictions are correct for first labeled bout!\n');
else
  error('Predictions are wrong for first label');
end
data.Predict(iExp,iFly,ts2(1),ts2(end));  % Make sure the prediction has been done
pred2=data.GetPredictedIdx(iExp,iFly,ts2(1),ts2(end));
if all(pred2.predictedidx==iBehavior2)
  fprintf('Predictions are correct for second labeled bout!\n');
else
  error('Predictions are wrong for second label');
end

% Clean up
data.closeJabFile();
data=[];  %#ok
success=true;

end

author: kristinbranson

repo: JAABA


text:
function [coords] = getCropCoords(thisFrame,varargin)
	% GUI to allow users to select crop coordinates.
	% Biafra Ahanonu
	% started: 2014.01.03 [19:13:01]
	% inputs
		%
	% outputs
		%

	[coords] = ciapkg.image.getCropCoords(thisFrame,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function PowerOverTime(nlp, sys)

    domains = sys.Gamma.Nodes.Domain;
    T  = SymVariable('t',[2, 1]);
    for i=1:numel(domains)
        domain = domains{i};

        u = domain.Inputs.Control.u;
        dq = domain.States.dx;
        B = domain.Gmap.Control.u;

        P = tovector(norm(dq.*(B*u)).^2)./(T(2)-T(1));
        power_fun = SymFunction(['powerOverTime_' sys.Gamma.Nodes.Name{i}],P,{T,u,dq});
        rs_phase = getPhaseIndex(nlp,sys.Gamma.Nodes.Name{i});
        addRunningCost(nlp.Phase(rs_phase),power_fun,{'T','u','dx'});

    end
    nlp.update;
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function [y,idx] = nanmin(a,dim,b)
% FORMAT: [Y,IDX] = NANMIN(A,DIM,[B])
%
%    Minimum ignoring NaNs
%
%    This function enhances the functionality of NANMIN as distributed in
%    the MATLAB Statistics Toolbox and is meant as a replacement (hence the
%    identical name).
%
%    If fact NANMIN simply rearranges the input arguments to MIN because
%    MIN already ignores NaNs.
%
%    NANMIN(A,DIM) calculates the minimum of A along the dimension DIM of
%    the N-D array X. If DIM is omitted NANMIN calculates the minimum along
%    the first non-singleton dimension of X.
%
%    NANMIN(A,[],B) returns the minimum of the N-D arrays A and B.  A and
%    B must be of the same size.
%
%    Comparing two matrices in a particular dimension is not supported,
%    e.g. NANMIN(A,2,B) is invalid.
%
%    [Y,IDX] = NANMIN(X,DIM) returns the index to the minimum in IDX.
%
%    Similar replacements exist for NANMAX, NANMEAN, NANSTD, NANMEDIAN and
%    NANSUM which are all part of the NaN-suite.
%
%    See also MIN

% -------------------------------------------------------------------------
%    author:      Jan Glscher
%    affiliation: Neuroimage Nord, University of Hamburg, Germany
%    email:       glaescher@uke.uni-hamburg.de
%
%    $Revision: 1.1 $ $Date: 2004/07/15 22:42:14 $

if nargin < 1
	error('Requires at least one input argument')
end

if nargin == 1
	if nargout > 1
		[y,idx] = min(a);
	else
		y = min(a);
	end
elseif nargin == 2
	if nargout > 1
		[y,idx] = min(a,[],dim);
	else
		y = min(a,[],dim);
	end
elseif nargin == 3
	if ~isempty(dim)
		error('Comparing two matrices along a particular dimension is not supported')
	else
		if nargout > 1
			[y,idx] = min(a,b);
		else
			y = min(a,b);
		end
	end
elseif nargin > 3
	error('Too many input arguments.')
end

% $Id: nanmin.m,v 1.1 2004/07/15 22:42:14 glaescher Exp glaescher $

author: mtex-toolbox

repo: mtex

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


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

test_suite = functiontests(localfunctions);

function testTriangle(testCase) %#ok<*DEFNU>

pts3d = [3 0 0; 0 6 0;0 0 9];
line1 = [0 0 0 1 2 3];

inter = intersectRayPolygon3d(line1, pts3d);
exp = [1 2 3];
testCase.assertEqual(exp, inter);

line2 = [10 0 0 1 2 3];

[inter, inside] = intersectRayPolygon3d(line2, pts3d); %#ok<ASGLU>
testCase.assertFalse(inside);

function testReverseTriangle(testCase)

pts3d = [3 0 0; 0 6 0;0 0 9];
line1 = [0 0 0 1 2 3];

inter = intersectRayPolygon3d(line1, pts3d);
exp = [1 2 3];
testCase.assertEqual(exp, inter);

function testLineArray(testCase)

pts3d = [3 0 0;0 0 9; 0 6 0];
lines = [0 0 0 3 6 9 ; 10 0 0 3 6 9; 3 6 9 3 6 9];

[inters, inside] = intersectRayPolygon3d(lines, pts3d);
testCase.assertEqual(3, size(inters, 1));
testCase.assertEqual([true;false;false], inside);

author: mattools

repo: matGeom


text:
function [ds, err, dsname, plotok] = dfcreatedataset(varargin)
%DFCREATEDATASET Create dfittool data set
%   DFCREATEDATASET(YEXPR,CEXPR,FEXPR,DSNAME) creates a data set named DSNAME
%   using the data obtained by evaluating the expressions entered in the gui
%   for Y, CENSORING, and FREQUENCY.
%
%   DFCREATEDATASET(YDATA,CDATA,FDATA,DSNAME,YNAME,CNAME,FNAME) creates a data
%   set named DSNAME using the y, censoring, and frequency data passed in via
%   the command line.  The "name" arguments are the names of these data
%   variables, or empty if they are not simple variable names.

%   $Revision: 1.1.6.7 $  $Date: 2004/01/24 09:35:24 $
%   Copyright 2003-2004 The MathWorks, Inc.

% Determine if we had data passed directly in from the command line
if nargin>=1 && isnumeric(varargin{1})
   fromgui = false;
else
   fromgui = true;
end

% Check for a data set name (we used to check for duplicate names
% before calling this function).

if nargin>=4
    dsname = varargin{4};
    if fromgui
        % Get the array of data sets
        dsdb = getdsdb;
        dset = down(dsdb);
        while(~isempty(dset))
            if strcmp(dset.name, dsname);
                err = 'datanamethesame';
                ds = [];
                plotok = false;
                return;
            end
            dset = right(dset);
        end
    end
else
   dsname = '';
end

n = min(length(varargin),3);
nameargs = cell(1,3);
if ~fromgui
   % Input data entered at the command line, maybe input names in args 5-7
   nameargs(:) = {''};
   valargs = cell(1,3);
   for j=1:n
      valargs{j} = varargin{j};
   end
   for j=5:nargin
      nameargs{j-4} = varargin{j};
   end
else
   % Input expressions from gui, get values by evaluating them
   nameargs = varargin(1:n);
   valargs = cell(1,0);
end

% Check input expressions or data vectors
outvecs = cell(1,3);
[err,outvecs{:}] = dfcheckselections(nameargs{:}, valargs{:});
if ~isequal(err, '')
   ds = [];
   return;
end

% Make a new data set
ds = stats.dfdata(nameargs{:},dsname,outvecs{:});
dsname = ds.name;

% Set the function type before plotting
setftype(ds,dfgetset('ftype'));
ds.conflev = dfgetset('conflev');

% Plot the histogram or other empirical curve if possible
dffig = dfgetset('dffig');
if dfcanplotdata(ds,dffig)
   ds.plot = 1;
   ds.plotok = 1;   % this was the default when we created ds
   plotok = 'true'; % use text; boolean causes problems in java caller
else
   % Update the java dialog to show current flag state
   % ds.plot = 0;   % this was the default when we created ds
   ds.plotok = 0;
   plotok = 'false';
end

% Update the plot
dfupdatexlim;
dfupdateylim;
dfupdatelegend(dffig);

author: zouchuhang

repo: LayoutNet


text:
% Copyright 2017 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     https://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

% Returns set of parameters
% set light_mode = true to run the code in a sub optimal but faster mode
% set light_mode = false to obtain the results reported in the RED paper

function params = GetUniformDeblurSDParams(light_mode)

% regularization factor
params.lambda = 0.02;

% number of outer iterations
if light_mode
    params.outer_iters = 500;
else
    params.outer_iters = 1500;
end

% level of noise assumed in the regularization-denoiser
params.effective_sigma = 3.25;

author: google

repo: RED

save_path: github-repos/MATLAB/google-RED


text:
function [indBoundary]=tesBoundary(varargin)

% function [indBoundary]=tesBoundary(F)
% ------------------------------------------------------------------------
%
% Change log:
% 2009
% 2019/04/24 Added varargin support since V can be skipped
% 2019/04/24 Started working on cell (e.g. mixed mesh) support, not
% completed yet
% 2020/07/09 Added mixed mesh support, works for pentahedra, needs work for
% general case
% 2021/10/08 Simplified
% 2021/10/08 No longer needs vertices as input
% ------------------------------------------------------------------------

%% Parse input

switch nargin
    case 1
        F=varargin{1};
    case 2
        F=varargin{1};
        warning('Second input (vertices) no longer required. Update code to avoid future error.');
end

%%

if isa(F,'cell')
    indBoundary=cell(size(F));
    for q=1:1:numel(F)
        indBoundary{q}=getBoundary(F{q});
    end
else
    [indBoundary]=getBoundary(F);
end

end

%%

function [indBoundary]=getBoundary(F)

Fs=sort(F,2); %Sort so faces with same nodes have the same rows
[~,~,~,F_count]=cunique(Fs,'rows'); %get indices for unique faces
indBoundary=find(F_count==1);

% Old method
%     Fbs=sort(F,2);
%     sizVirt=numPoints*ones(1,size(Fbs,2));
%     ind_F = sub2indn(sizVirt,Fbs);
%     [~,indUni1,~]=unique(Fbs,'rows'); %Get indices for unique faces
%     ind_F_uni=ind_F(indUni1,:);
%     ind=1:1:size(Fbs,1);
%     ind=ind(~ismember(ind,indUni1));
%     ind_Fb_cut=ind_F(ind,:);
%     L_uni=~ismember(ind_F_uni,ind_Fb_cut);
%     indBoundary=indUni1(L_uni,:);
end

%%
% _*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