text:
function [trl, event] = ft_trialfun_artinis(cfg)

% FT_TRIALFUN_ARTINIS is adjoining the result of ft_trialfun_general and
% those events found by FT_TRIALFUN_GENERAL.
%
% The trialdef structure can contain the following specifications
%   cfg.trialdef.eventtype  = 'string'
%   cfg.trialdef.eventvalue = number, string or list with numbers or strings
%   cfg.trialdef.oxyproj    = 'string', indicating an oxyproj-file, in
%                             which information about the events for this
%                             oxy3-file are stored
%   cfg.trialdef.prestim    = latency in seconds (optional)
%   cfg.trialdef.poststim   = latency in seconds (optional)
%
% If you want to read all data from a continuous file in segments, you can specify
%    cfg.trialdef.triallength = duration in seconds (can be Inf)
%    cfg.trialdef.ntrials     = number of trials
%
% If you specify
%   cfg.trialdef.eventtype  = '?'
% a list with the events in your datafile will be displayed on screen.
%
% If you specify
%   cfg.trialdef.eventtype = 'gui'
% a graphical user interface will allow you to select events of interest.
%
% See also FT_TRIALFUN_GENERAL, FT_DEFINETRIAL, FT_PREPROCESSING

% You are using the FieldTrip NIRS toolbox developed and maintained by
% Artinis Medical Systems (http://www.artinis.com). For more information
% on FieldTrip, see http://www.fieldtriptoolbox.org
%
% This work is licensed under a Creative Commons Attribution-ShareAlike 4.0
% International License. To view a copy of this license, visit
% http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to
% Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
%
% Creative Commons Attribution-ShareAlike 4.0 International License:
% -----------------------------------
% You are free to:
%
%     Share - copy and redistribute the material in any medium or format
%     Adapt - remix, transform, and build upon the material
%     for any purpose, even commercially.
%
%     The licensor cannot revoke these freedoms as long as you follow the
%     license terms.
%
% Under the following terms:
%
%     Attribution - You must give appropriate credit, provide a link to
%                    the license, and indicate if changes were made. You
%                    may do so in any reasonable manner, but not in any way
%                    that suggests the licensor endorses you or your use.
%
%     ShareAlike - If you remix, transform, or build upon the material,
%                   you must distribute your contributions under the same
%                   license as the original.
%
%     No additional restrictions - You may not apply legal terms or
%                                   technological measures that legally
%                                   restrict others from doing anything the
%                                   license permits.
%
% -----------------------------------
%
% This toolbox is not to be used for medical or clinical purposes.
%
% Copyright (c) 2016 by Artinis Medical Systems.
% Contact: askforinfo@artinis.com
%
% Main programmer:
% Jörn M. Horschig, Artinis Medical Systems BV, http://www.artinis.com
% $Id$

if ~isfield(cfg.trialdef, 'oxyproj') || ~strcmp(cfg.trialdef.oxyproj(end-7:end), '.oxyproj') || ~exist(cfg.trialdef.oxyproj, 'file')
  error('You need to provide an existing oxyproj-file as input to use this function (cfg.trialdef.oxyproj)')
end

cfg.eventformat = ft_getopt(cfg, 'eventformat');
cfg.headerformat = ft_getopt(cfg, 'headerformat');
cfg.dataformat = ft_getopt(cfg, 'dataformat');

fprintf('reading the events from ''%s''\n', cfg.headerfile);
event = ft_read_event(cfg.headerfile, 'headerformat', cfg.headerformat, 'eventformat', cfg.eventformat, 'dataformat', cfg.dataformat);

fprintf('reading the events from ''%s''\n', cfg.trialdef.oxyproj);
measurement = parseoxyproj(cfg.trialdef.oxyproj, cfg.datafile);

if numel(measurement)>1
  warning('More than one measurement found in the oxyproj file - removing all but the first one\n');
  measurement(2:end) = [];
end

% read the header for obtaining the sample rate
hdr = ft_read_header(cfg.headerfile, 'headerformat', cfg.headerformat);

% create the trl and event-structure with the new information
trl = [];
for i=1:numel(measurement.events.names)
  add2trl = strcmp(cfg.trialdef.eventtype, 'oxyproj') && strcmp(cfg.trialdef.eventvalue, measurement.events.names{i});
  for e=1:numel(measurement.events.onsets{i})
    % add event
    event(end+1).type      = 'oxyproj';
    event(end).sample    = measurement.events.onsets{i}(e);
    event(end).value     = measurement.events.names{i};
    event(end).offset    = 0;
    event(end).duration  = measurement.events.durations{i}(e);
    if add2trl
      trloff = round(-cfg.trialdef.prestim * hdr.Fs);
      trlbeg = event(end).sample+trloff;
      % this will not work if prestim was not defined, the code will then crash
      trldur = round((cfg.trialdef.poststim+cfg.trialdef.prestim)*hdr.Fs) - 1;
      trlend = trlbeg + trldur;
      trl(end+1, :) = [trlbeg trlend trloff];
    end
  end
end

% create trials
cfg.event = event;
[trl, event] = ft_trialfun_general(cfg);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% dimension of ambient space
n = 6;
% number of subspaces = number of clusters
ns = 2;
% dimensions of individual subspaces 1 and 2
d1  = 2;
d2 = 2;
% number of signals in individual subspaces
s1 = 10;
s2 = 10;
% total number of signals
S = s1 + s2;
% A random basis for first subspace;
basis1 = randn(n,d1);
% coefficients for s1 vectors chosen randomly in subspace 1
coeffs1 = randn(d1,s1);
% Random signals from first subspace
X1 = basis1 * coeffs1;
% A random basis for second subspace
basis2 = randn(n,d2);
% coefficients for s2 vectors chosen randomly in subspace 2
coeffs2 = randn(d2,s2);
% Random signals from first subspace
X2 = basis2 * coeffs2;
% Prepare the overall set of signals
X = [X1 X2];
% ground through clustering data
true_labels = [1*ones(s1,1) ; 2*ones(s2,1)];
% the largest dimension amongst all subspaces
K = max(d1, d2);
X = spx.norm.normalize_l2(X);
C = spx.fast.gomp_spr(X, K, 2);

solver = spx.pursuit.single.OrthogonalMatchingPursuit(X, K);
C2 = zeros(S, S);
for s=1:S
    solver.IgnoredAtom = s;
    result = solver.solve(X(:, s));
    C2(:, s) = result.z;
end
max(C2 - C)
max(abs(X - X * C))

author: indigits

repo: sparse-plex


text:
function mm_io_test08 ( )

%*****************************************************************************80
%
%% MM_IO_TEST08 tests MM_VALUES_READ.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MM_IO_TEST08:\n' );
  fprintf ( 1, '  MM_VALUE_READ reads the value records \n' );
  fprintf ( 1, '  of a Matrix Market file.\n' );

  input_filename = 'mm_io_test07.mm';

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MM_IO_TEST08 - Fatal error!\n' );
    fprintf ( 1, '  Could not open the input file.\n' );
    error ( 'MM_IO_TEST08 - Fatal error!' );
  end

  [ id, type, rep, field, symm ] = mm_header_read ( input_unit );

  while ( 1 )

    comment = mm_comment_read ( input_unit );

    if ( comment == -1 )
      break;
    end

    if ( comment(1) ~= '%' )
      break;
    end

  end
%
%  The COMMENT line is actually the line containing the size information.
%
  [ nrow, ncol, nnz ] = mm_size_read_string ( comment, rep, symm );

  [ indx, jndx, value ] = mm_values_read ( input_unit, rep, field, nnz );

  fclose ( input_unit );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here is the value information read from the file:\n' );
  fprintf ( 1, '\n' );

  mm_values_print ( rep, field, nnz, indx, jndx, value );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'MM_IO_TEST08:\n' );
  fprintf ( 1, '  The MM file "%s" containing \n', input_filename );
  fprintf ( 1, '  the header, comments, size and values, has been read.\n' );
  fprintf ( 1, '\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function trans = translation3d(varargin)
%TRANSLATION3D return 4x4 matrix of a 3D translation.
%
%   usage :
%   TRANS = translation3d(DX, DY, DZ);
%   return the translation corresponding to DX and DY.
%   The returned matrix has the form :
%   [1 0 0 DX]
%   [0 1 0 DY]
%   [0 0 1 DZ]
%   [0 0 0  1]
%
%   TRANS = translation3d(VECT);
%   return the translation corresponding to the given vector [x y z].
%
%
%   See also:
%   vectors3d, transforms3d, transformPoint, rotation

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2004-04-06
% Copyright 2004 INRA - TPV URPOI - BIA IMASTE

% deprecation warning
warning('geom3d:deprecated', ...
    [mfilename ' is deprecated, use ''createTranslation3d'' instead']);

% call current implementation
trans = createTranslation3d(varargin{:});

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function nu = PoissonRatio(S,x,y)
% Poisson ratio of an elasticity tensor
%
% Syntax
%   nu = PoissonRatio(S,x,y)
%   nu = PoissonRatio(S) % isotropic case
%
% Input
%  S - elastic @complianceTensor
%  x - @vector3d tension direction
%  y - @vector3d perpendicular direction
%
% Output
%  nu - Poisson ratio
%
% Description
%
% $$\nu = \frac{-S_{ijkl} x_i x_j y_k y_l}{S_{mnop} x_m x_n x_o x_p}$$
%

if nargin == 1  % isotrpoic theory

  S_iso_Voigt = mean(S,uniformODF(S.CS));
  nu = S_iso_Voigt.PoissonRatio(xvector,yvector);

  return
end

% generate a function if required
if nargin == 1 || isempty(x)

  nu = S2FunHarmonicSym.quadrature(@(v) PoissonRatio(S,v,y),'bandwidth',4,S.CS);

elseif nargin <= 2 || isempty(y)

  nu = S2FunHarmonicSym.quadrature(@(v) PoissonRatio(S,x,v),'bandwidth',4,S.CS);

else

  % compute tensor product
  nu = -EinsteinSum(S,[-1 -2 -3 -4],x,-1,x,-2,y,-3,y,-4) ./ ...
    EinsteinSum(S,[-1 -2 -3 -4],x,-1,x,-2,x,-3,x,-4);

end

author: mtex-toolbox

repo: mtex

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


text:
function [ y, next ] = p19_equil ( neqn, next )

%*****************************************************************************80
%
%% P19_EQUIL returns equilibrium solutions of problem p19.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, integer NEXT, the index of the previous
%    equilibrium, which should be 0 on first call.
%
%    Output, real Y(NEQN), the "next" equilibrium solution, if any.
%
%    Output, integer NEXT, the index of the current equilibrium,
%    or 0 if there are no more.
%)
  if ( next == 0 )
    next = 1;
    y(1:neqn,1) = 0.0;
  else
    next = 0;
    y = [];
  end

  return
end


text:
classdef MimCapture2DImage < MimGuiPlugin
    % MimCapture2DImage. Gui Plugin for exporting the image currently in the
    % visualisation window to a file
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimCapture2DImage is a Gui Plugin for the MIM Toolkit.
    %     The gui will create a button for the user to run this plugin.
    %     Running this plugin will capture the current viewer image and save
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Capture'
        SelectedText = 'Capture'

        ToolTip = 'Save image and overlay view to files'
        Category = 'Export'
        Visibility = 'Dataset'
        Mode = 'Toolbar'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1

        Icon = 'camera.png'
        Location = 19
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            drawnow;
            gui_app.Capture;
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded;
        end

        function is_selected = IsSelected(gui_app)
            is_selected = false;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [tfr,rtfr,hat] = tfrrpwv(x,t,N,h,trace);
%TFRRPWV Reassigned  pseudo Wigner-Ville distribution.
%	[TFR,RTFR,HAT] = TFRRPWV(X,T,N,H,TRACE)
%	computes the pseudo Wigner-Ville distribution
%	and its reassigned version.
%
%	X     : analysed signal,
%	T     : the time instant(s)      (default : 1:length(X)).
%	N     : number of frequency bins (default : length(X)).
%	H     : frequency smoothing window, H(0) being forced to 1
%	                                 (default : Hamming(N/4)).
%	TRACE : if nonzero, the progression of the algorithm is shown
%	                                 (default : 0).
%	TFR,  : time-frequency representation and its reassigned
%	RTFR    version. When called without output arguments,
%	        TFRRPWV runs TFRQVIEW.
%	HAT   : Complex matrix of the reassignment vectors.
%
%	Example:
%	 sig=fmlin(128,0.1,0.4); t=1:2:128;
%	 h=tftb_window(17,'Kaiser'); tfrrpwv(sig,t,64,h,1);
%
%	See also all the time-frequency representations listed in
%	 the file CONTENTS (TFR*)

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

if (nargin == 0),
 error('At least 1 parameter required');
end;
[xrow,xcol] = size(x);
if (nargin < 1),
 error('At least 1 parameter is required');
elseif (nargin <= 2),
 N=xrow;
end;

hlength=floor(N/4);
if (rem(hlength,2)==0),
 hlength=hlength+1;
end;

if (nargin == 1),
 t=1:xrow; h = tftb_window(hlength); trace=0;
elseif (nargin == 2)|(nargin == 3),
 h = tftb_window(hlength); trace=0;
elseif (nargin == 4),
 trace = 0;
end;

if (N<0),
 error('N must be greater than zero');
end;
[trow,tcol] = size(t);
if (xcol~=1),
 error('X must have only one column');
elseif (trow~=1),
 error('T must only have one row');
elseif (2^nextpow2(N)~=N),
 fprintf('For a faster computation, N should be a power of two\n');
end;

[hrow,hcol]=size(h); Lh=(hrow-1)/2; h=h/h(Lh+1);
if (hcol~=1)|(rem(hrow,2)==0),
 error('H must be a smoothing window with odd length');
end;

if (tcol==1),
 Dt=1;
else
 Deltat=t(2:tcol)-t(1:tcol-1);
 Mini=min(Deltat); Maxi=max(Deltat);
 if (Mini~=Maxi),
  error('The time instants must be regularly sampled.');
 else
  Dt=Mini;
 end;
 clear Deltat Mini Maxi;
end;

tfr= zeros(N,tcol); tf2= zeros(N,tcol);
if trace, disp('Pseudo Wigner-Ville distribution'); end;
Dh=dwindow(h);
for icol=1:tcol,
 ti= t(icol); taumax=min([ti-1,xrow-ti,round(N/2)-1,Lh]);
 tau=-taumax:taumax; indices= rem(N+tau,N)+1;
 if trace, disprog(icol,tcol,10); end;
 tfr(indices,icol)= h(Lh+1+tau).*x(ti+tau).*conj(x(ti-tau));
 tf2(indices,icol)=Dh(Lh+1+tau).*x(ti+tau).*conj(x(ti-tau));
 tau=round(N/2);
 if (ti<=xrow-tau)&(ti>=tau+1)&(tau<=Lh),
  tfr(tau+1,icol) = 0.5 * ( h(Lh+1+tau) * x(ti+tau,1) * conj(x(ti-tau,xcol))  + ...
                            h(Lh+1-tau) * x(ti-tau,1) * conj(x(ti+tau,xcol))) ;
  tf2(tau+1,icol) = 0.5 * (Dh(Lh+1+tau) * x(ti+tau,1) * conj(x(ti-tau,xcol))  + ...
                           Dh(Lh+1-tau) * x(ti-tau,1) * conj(x(ti+tau,xcol))) ;
 end;
end ;
tfr= real(fft(tfr));
tf2=imag(fft(tf2));
tfr=tfr(:);tf2=tf2(:);
avoid_warn=find(tfr~=0);
tf2(avoid_warn)=round(N*tf2(avoid_warn)./tfr(avoid_warn)/(2.0*pi));
%tf2= round(N*imag(fft(tf2))./tfr/(2.0*pi));
if trace, fprintf ('\nreassignment: \n'); end;
tfr=reshape(tfr,N,tcol);
tf2=reshape(tf2,N,tcol);

rtfr= zeros(N,tcol);
Ex=mean(abs(x(min(t):max(t))).^2); Threshold=1.0e-6*Ex;
for icol=1:tcol,
 if trace, disprog(icol,tcol,10); end;
 for jcol=1:N,
  if abs(tfr(jcol,icol))>Threshold,
   jcolhat= jcol - tf2(jcol,icol);
   jcolhat=rem(rem(jcolhat-1,N)+N,N)+1;
   rtfr(jcolhat,icol)=rtfr(jcolhat,icol) + tfr(jcol,icol) ;
   tf2(jcol,icol)=jcolhat;
  else
   tf2(jcol,icol)=inf;
   rtfr(jcol,icol)=rtfr(jcol,icol) + tfr(jcol,icol) ;
  end;
 end;
end;

if trace, fprintf('\n'); end;
if (nargout==0),
 TFTBcontinue=1;
 while (TFTBcontinue==1),
  choice=menu ('Choose the representation:',...
               'stop',...
               'pseudo Wigner-Ville distribution',...
               'reassigned pseudo Wigner-Ville distribution');
  if (choice==1), TFTBcontinue=0;
  elseif (choice==2),
   tfrqview(tfr,x,t,'tfrpwv',h);
  elseif (choice==3),
   tfrqview(rtfr,x,t,'tfrrpwv',h);
  end;
 end;
elseif (nargout>2),
 hat=tf2;
end;

author: HeLiangHIT


text:
clc;
clear all;

load('..\data\nlc_data_5zones.mat')
load('..\data\lc_data_20s_withpoints.mat')

% transfer to related coordinate, in terms of the host
%%% for the lc_data
r_lc_data = {};

Label = [points; 600*ones(length(nlc_data),2)];

for i = 1: length(lc_data)
    trajs = lc_data{i};
    origin_x = trajs.veh_s.x(1);
    origin_y = trajs.veh_s.y(1);

    R_trajs.veh_s_x = trajs.veh_s.x(1:100) - origin_x;
    R_trajs.veh_s_y = trajs.veh_s.y(1:100) - origin_y;
    Label(i,1) = trajs.veh_s.y(points(i,1)) - origin_y;
    Label(i,2) = trajs.veh_s.y(points(i,2)) - origin_y;

    da = [];

    if (trajs.veh_r.x(1) ~= 10000)
        R_trajs.veh_r_x = trajs.veh_r.x(1:100) - origin_x;
        R_trajs.veh_r_y = trajs.veh_r.y(1:100) - origin_y;
    else
        R_trajs.veh_r_x = 0*trajs.veh_r.x(1:100);
        R_trajs.veh_r_y = 0*trajs.veh_r.y(1:100);
    end
    if (trajs.veh_f.x(1) ~= 10000)
        R_trajs.veh_f_x = trajs.veh_f.x(1:100) - origin_x;
        R_trajs.veh_f_y = trajs.veh_f.y(1:100) - origin_y;
    else
        R_trajs.veh_f_x = 0*trajs.veh_f.x(1:100);
        R_trajs.veh_f_y = 0*trajs.veh_f.y(1:100);
    end

    if (trajs.veh_rt.x(1) ~= 10000)
        R_trajs.veh_rt_x = trajs.veh_rt.x(1:100) - origin_x;
        R_trajs.veh_rt_y = trajs.veh_rt.y(1:100) - origin_y;
    else
        R_trajs.veh_rt_x = 0*trajs.veh_rt.x(1:100);
        R_trajs.veh_rt_y = 0*trajs.veh_rt.y(1:100);
    end

    if (trajs.veh_ft.x(1) ~= 10000)
        R_trajs.veh_ft_x = trajs.veh_ft.x(1:100) - origin_x;
        R_trajs.veh_ft_y = trajs.veh_ft.y(1:100) - origin_y;
    else
        R_trajs.veh_ft_x = 0*trajs.veh_ft.x(1:100);
        R_trajs.veh_ft_y = 0*trajs.veh_ft.y(1:100);
    end

    if isfield(trajs, 'veh_st') && (trajs.veh_st.x(1) ~= 10000)
        R_trajs.veh_st_x = trajs.veh_st.x(1:100) - origin_x;
        R_trajs.veh_st_y = trajs.veh_st.y(1:100) - origin_y;
    else
        R_trajs.veh_st_x = 0*trajs.veh_ft.x(1:100);
        R_trajs.veh_st_y = 0*trajs.veh_ft.y(1:100);
    end

    da = [da;R_trajs.veh_s_x(end-19:end-10) ];
    da = [da;R_trajs.veh_f_x(end-19:end-10) ];
    da = [da;R_trajs.veh_rt_x(end-19:end-10) ];
    da = [da;R_trajs.veh_ft_x(end-19:end-10) ];
    da = [da;R_trajs.veh_st_x(end-19:end-10) ];

    r_lc_data{i}=R_trajs;
end

r_nlc_data = {};
for i = 1: length(nlc_data)
    da = [];
    trajs = nlc_data{i};
    origin_x = trajs.veh_s.x(1);
    origin_y = trajs.veh_s.y(1);

    R_trajs.veh_s_x = trajs.veh_s.x - origin_x;
    R_trajs.veh_s_y = trajs.veh_s.y - origin_y;

    if (trajs.veh_r.x(1) ~= 10000)
        R_trajs.veh_r_x = trajs.veh_r.x - origin_x;
        R_trajs.veh_r_y = trajs.veh_r.y - origin_y;
    else
        R_trajs.veh_r_x = 0*trajs.veh_r.x;
        R_trajs.veh_r_y = 0*trajs.veh_r.y;
    end

    if (trajs.veh_f.x(1) ~= 10000)
        R_trajs.veh_f_x = trajs.veh_f.x - origin_x;
        R_trajs.veh_f_y = trajs.veh_f.y - origin_y;
    else
        R_trajs.veh_f_x = 0*trajs.veh_f.x;
        R_trajs.veh_f_y = 0*trajs.veh_f.y;
    end

    if (trajs.veh_rt.x(1) ~= 10000)
        R_trajs.veh_rt_x = trajs.veh_rt.x - origin_x;
        R_trajs.veh_rt_y = trajs.veh_rt.y - origin_y;
    else
        R_trajs.veh_rt_x = 0*trajs.veh_rt.x;
        R_trajs.veh_rt_y = 0*trajs.veh_rt.y;
    end

    if (trajs.veh_ft.x(1) ~= 10000)
        R_trajs.veh_ft_x = trajs.veh_ft.x - origin_x;
        R_trajs.veh_ft_y = trajs.veh_ft.y - origin_y;
    else
        R_trajs.veh_ft_x = 0*trajs.veh_ft.x;
        R_trajs.veh_ft_y = 0*trajs.veh_ft.y;
    end

    if isfield(trajs, 'veh_st') && (trajs.veh_st.x(1) ~= 10000)
        R_trajs.veh_st_x = trajs.veh_st.x - origin_x;
        R_trajs.veh_st_y = trajs.veh_st.y - origin_y;
    else
        R_trajs.veh_st_x = 0*trajs.veh_ft.x;
        R_trajs.veh_st_y = 0*trajs.veh_ft.y;
    end
    da = [da;R_trajs.veh_s_x(end-19:end-10) ];
    da = [da;R_trajs.veh_f_x(end-19:end-10) ];
    da = [da;R_trajs.veh_rt_x(end-19:end-10) ];
    da = [da;R_trajs.veh_ft_x(end-19:end-10) ];
    da = [da;R_trajs.veh_st_x(end-19:end-10) ];

    r_nlc_data{i}=R_trajs;
end

save('../data/processed_relative.mat','Label','r_lc_data','r_nlc_data')

author: donnydcy


text:
function v = cvx_vexity( x )

global cvx___
s = x.size_;
if any( s == 0 ),
    v = cvx_zeros( s );
    return
end
p  = cvx___.vexity;
b  = x.basis_;
n  = length( p );
nb = size( b, 1 );
if nb < n,
    p = p( 1 : nb, 1 );
elseif n < nb,
    p( nb, 1 ) = 0;
end
if ~any( p ),
    v = cvx_zeros( x.size_ );
    if x.slow_,
        v( isnan( x.basis_( 1, : ) ) ) = NaN;
    end
    return
end
b = b( p ~= 0, : );
p = nonzeros(p).';
if cvx___.nan_used,
    b = sparse( b );
end
v = full( p * b );
tt = abs( v ) ~= abs( p ) * abs( b );
if x.slow_,
    v( tt | isnan( x.basis_( 1, : ) ) ) = NaN;
else
    v( tt ) = NaN;
end
if ~isreal( x ),
    v( any( imag( x.basis_ ), 1 ) & v ) = NaN;
end
v = sign( v );
v = reshape( v, x.size_ );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [x, infos] = alternating_onmf(V, rank, in_options)
% Orthogonal two-block coordinate descent (2-BCD) for non-negative matrix factorization (Alt-Orth-NMF).
%
% The problem of interest is defined as
%
%       min || V - WH ||_F^2,
%       where
%       {V, W, H} >= 0, and W or H is orthogonal, i.e., HH^T = I_r.
%
% Given a non-negative matrix V, factorized non-negative matrices {W, H} are calculated.
%
%
% Inputs:
%       matrix      V
%       rank        rank
%
% Output:
%       w           solution of w
%       infos       information
%
% References:
%       F. Pompili, N. Gillis, P.-A. Absil and F. Glineur,
%       "Two Algorithms for Orthogonal Nonnegative Matrix Factorization
%       with Application to Clustering,"
%       Neurocomputing 141, pp. 15-25, 2014.
%
%
% This file is part of NMFLibrary.
%
% This file has been ported from
% alternatingONMF.m at https://gitlab.com/ngillis/nmfbook/-/tree/master/algorithms
% by Nicolas Gillis (nicolas.gillis@umons.ac.be)
%
% Change log:
%
%       June 14, 2022 (Hiroyuki Kasai): Ported initial version
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.orth_h = true;
    local_options.delta = 0;
    local_options.special_stop_condition = @(epoch, infos, options, stop_options) alt_nmf_stop_func(epoch, infos, options, stop_options);

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    W = init_factors.W; % use only W

    % initialize
    method_name = 'ALT-ONMF';
    epoch = 0;
    grad_calc_count = 0;

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    % initialize for this algorithm
    % Xn: normalized version of X, ||Xn(:,j)||_2  1 for all j
    norm2x = sqrt(sum(V.^2,1));
    Vn = V .* repmat(1./(norm2x+1e-16), m, 1);
    normV2 = sum(V(:).^2);
    stop_options.normV2 = normV2;
    H = nnls_orth(V, W, Vn);

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(V, W, H, [], options, [], epoch, grad_calc_count, 0);
    orth_val = norm(H*H' - eye(rank), 'fro');
    [infos.orth] = orth_val;
    infos.prev_error =  Inf;

    if options.verbose > 1
        fprintf('%s: Epoch = 0000, cost = %.16e, optgap = %.4e\n', method_name, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag, infos] = check_stop_condition(epoch, infos, options, stop_options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        % update H
        % H = argmin_H ||X-WH||_F, H >= 0, rows H orthogonal up to a scaling of the rows of H
        H = nnls_orth(V, W, Vn);

        % normalize rows of H
        norm2h = sqrt(sum(H'.^2, 1)) + 1e-16;
        H = repmat(1./norm2h', 1, n) .* H;

        % update W
        W = V * H';

        % measure gradient calc count
        grad_calc_count = grad_calc_count + m*n;

        % measure elapsed time
        elapsed_time = toc(start_time);

        % update epoch
        epoch = epoch + 1;

        % store info
        infos = store_nmf_info(V, W, H, [], options, infos, epoch, grad_calc_count, elapsed_time);
        orth_val = norm(H*H' - eye(rank), 'fro');
        [infos.orth] = [infos.orth orth_val];

        % display info
        display_info(method_name, epoch, infos, options);

        % check convergence (by original code)
        e = sqrt( (normV2-sum(sum(W.^2)))/normV2 );
        if (epoch > 2) && (abs(e_prev-e) < options.delta)
            %break;
        else
            e_prev = e;
        end
    end

    x.W = W;
    x.H = H;

end

function [stop_flag, reason, infos] = alt_nmf_stop_func(epoch, infos, options, stop_options)

    stop_flag = false;
    reason = [];

    normV2 = stop_options.normV2;
    W = infos.final_W;

    error = sqrt( (normV2-sum(sum(W.^2)))/normV2 );
    if (epoch > 2) && (abs(infos.prev_error - error) < options.delta)
        stop_flag = true;
        reason = sprintf('Relative solution change tolerance reached: %.4e < %.4e\n', abs(infos.prev_error - error), options.delta);
    else
        infos.prev_error = error;
    end

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Torr Vision Group (TVG)
%  University of Oxford - UK
%
%  Qizhu Li <liqizhu@robots.ox.ac.uk>
%  August 2018
% ------------------------------------------------------------------------
% This file is part of the weakly-supervised training method presented in:
%    Qizhu Li*, Anurag Arnab*, Philip H.S. Torr,
%    "Weakly- and Semi-Supervised Panoptic Segmentation,"
%    European Conference on Computer Vision (ECCV) 2018.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
% This script demos generation of iterative ground truths for weakly-
% supervised experiments.
%
% It post-processes network predictions to produce iterative GT for
% training in the next iteration.
%  INPUT:
%  - results/pred_flat_feat/*.mat: the softmax scores of the predicted
%  classes.
%  - results/pred_sem_raw/*.png: the prediction made by the current
%  weakly-supervised model.
%  - results/mcg_and_grabcut/<split>/*.png: the combined cues from MCG
%  and Grabcut. Optional. Set opts.run_merge_with_mcg_and_grabcut = false
%  to disable. To reproduce the results in our paper, disable after first
%  5 iterations.
% ------------------------------------------------------------------------

clearvars;
addpath scripts
addpath utils
addpath visualisation

% Extract detections from Cityscapes ground truth file and save as .mat
% this only need to be done once, as they do not change over iterative
% training stages
demo_instanceTrainId_to_dets;

clearvars;
dataset = 'cityscapes';
split = 'train';
opts = get_opts(dataset, split);
opts.list_path = 'lists/demo_id.txt';
opts.visualise_results = true;
[opts, results] = run_sub(opts);

% visualise
visualise_results_iterative_gt(opts, results);

author: qizhuli

repo: Weakly-Supervised-Panoptic-Segmentation


text:
classdef IMMOEA_F2 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP for testing IM-MOEA

%------------------------------- Reference --------------------------------
% R. Cheng, Y. Jin, K. Narukawa, and B. Sendhoff, A multiobjective
% evolutionary algorithm using Gaussian process-based inverse modeling,
% IEEE Transactions on Evolutionary Computation, 2015, 19(6): 838-856.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            t = (1+5*repmat(2:obj.D,size(X,1),1)/obj.D).*X(:,2:obj.D) - repmat(X(:,1),1,obj.D-1);
            g = 1 + 9*mean(t.^2,2);
            PopObj(:,1) = X(:,1);
            PopObj(:,2) = g.*(1-(PopObj(:,1)./g).^2);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - R(:,1).^2;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function flatROI = vol2flatROILevels(volROI,volView,flatView)
%
% flatROI = vol2flatROILevels(volROI,volView,flatView)
%
% Creates a flat ROI from a volume ROI by looking up the
% corresponding coords. Will create coordinates in the
% 'separate gray levels' slices.
%
% flatROI and volROI are ROI structures, like those found in
% view.ROIs
%
% volView must be the VOLUME structure.
% flatView must be the FLAT structure.
%
% ras, 10/04. Version for flat level views, off vol2flatROI.
% Intersect volROI.coords with volView.coords.  When in gray mode,
% this restricts the ROI to the gray matter.
ROIcoords = volROI.coords;
ROIcoords = intersectCols(ROIcoords,volView.coords);

% init coords for flat ROI
coords = [];

% get grayCoords for flat, across all slices
grayCoords = []; flatCoords = [];
for slice = 1:numSlices(flatView)
    grayCoords = [grayCoords flatView.grayCoords{slice}];
    sliceCoords = flatView.coords{slice};
    sliceCoords(3,:) = slice;
    flatCoords = [flatCoords sliceCoords];
end

[bothCoords ia ib] = intersectCols(ROIcoords,grayCoords);
coords = flatCoords(:,ib);

% % error check. ROIcoords should be contained entirely within gray coords.
% if (size(ROIcoords,2) ~= size(coords,2)) & (strcmp(volView.viewType,'Gray'))
%    fprintf(['\nGray nodes loaded from gray classification file are ',...
%          'incompatible with those loaded from the flat.mat file. ',...
%          'Rebuild the Gray and Flat view: rm Gray/*.mat Flat/*.mat']);
% end

% Remove duplicates (there shouldn't be any, but just to make sure).
coords = intersectCols(coords,coords);

% Set the fields
flatROI.coords = coords;
flatROI.name = volROI.name;
flatROI.color = volROI.color;
flatROI.viewType = flatView.viewType;

return

% from the older intersect code:
% flatCoords = round(flatView.coords(:,ib));
% flatCoords = [flatCoords; h*ones([1,size(flatCoords,2)])];
% coords = [coords,flatCoords];

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [point, ind] = util_get_annopoint_by_id(points, id)
point = [];
for i=1:length(points)
    if (points(i).id == id)
        point = points(i);
        ind = i;
        return;
    end
end
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [imgs, sess_images, spersess] = scn_get_image_names(nruns, ndisdaqs)

    if nargin < 2, ndisdaqs = 0;
    else
        fprintf('Removing first %03d image names from each session.\n', ndisdaqs);
    end

 % get image names
    % --------------------------------------------------------------
    filestr = ['r[1-' num2str(nruns) ']' filesep 'vols.V*.img'];
    disp(['Getting images: ' filestr])
    imgs = filenames(filestr, 'char', 'absolute');

    clear sess_images

    for i = 1:nruns
        filestr = ['r[' num2str(i) ']' filesep 'vols.V*.img'];
        sess_images{i} = filenames(filestr, 'char', 'absolute');

        % disdaqs
        if ~isempty(sess_images{i}) && size(sess_images{i}, 1) > ndisdaqs && ndisdaqs
            sess_images{i}(1:ndisdaqs, :) = [];
        end

    end

    wh_empty = false(1, nruns);
    spersess = zeros(1, nruns);
    for i = 1:nruns
        if isempty(sess_images{i}) || size(sess_images{1}, 1) == 0  % do seem to need 2nd expression
            wh_empty(i) = true;
        else
            spersess(i) = size(sess_images{1}, 1);
        end
    end

    if ndisdaqs
        % re-get all image names, w.o disdaqs
        imgs = char(sess_images{:});
    end

    fprintf('Sessions (runs) with images: %03d\tSessions without images:%03d\n', sum(~wh_empty), sum(wh_empty));
    spersess(wh_empty) = [];
    fprintf('Images per Session (spersess):\n');
    fprintf('%3.0f\t', spersess);
    fprintf('\n');

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function p_sign = perm_sign ( n, p )

%*****************************************************************************80
%
%% PERM_SIGN returns the sign of a permutation.
%
%  Discussion:
%
%    A permutation can always be replaced by a sequence of pairwise
%    transpositions.  A given permutation can be represented by
%    many different such transposition sequences, but the number of
%    such transpositions will always be odd or always be even.
%    If the number of transpositions is even or odd, the permutation is
%    said to be even or odd.
%
%  Example:
%
%    Input:
%
%      N = 9
%      P = 2, 3, 9, 6, 7, 8, 5, 4, 1
%
%    Output:
%
%      P_SIGN = +1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2007
%
%  Author:
%
%    Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Albert Nijenhuis, Herbert Wilf,
%    Combinatorial Algorithms,
%    Academic Press, 1978, second edition,
%    ISBN 0-12-519260-6.
%
%  Parameters:
%
%    Input, integer N, the number of objects permuted.
%
%    Input, integer P(N), a permutation, in standard index form.
%
%    Output, integer P_SIGN, the "sign" of the permutation.
%    +1, the permutation is even,
%    -1, the permutation is odd.
%
  ierror = perm_check ( n, p );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PERM_SIGN - Fatal error!\n' );
    fprintf ( 1, '  The input array does not represent\n' );
    fprintf ( 1, '  a proper permutation.  In particular, the\n' );
    fprintf ( 1, '  array is missing the value %d\n', ierror );
    error ( 'PERM_SIGN - Fatal error!' );
  end
%
%  Make a temporary copy of the permutation.
%
  q(1:n) = p(1:n);
%
%  Start with P_SIGN indicating an even permutation.
%  Restore each element of the permutation to its correct position,
%  updating P_SIGN as you go.
%
  p_sign = 1;

  for i = 1 : n - 1

    j = i4vec_index ( n, q, i );

    if ( j ~= i )
      temp = q(i);
      q(i) = q(j);
      q(j) = temp;
      p_sign = - p_sign;
    end

  end

  return
end


text:
function xyz_io_test ( )

%*****************************************************************************80
%
%% XYZ_IO_TEST tests the XYZ_IO library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 January 2009
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'XYZ_IO_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the XYZ_IO library.\n' );

  xyz_io_test01 ( );
  xyz_io_test02 ( );
  xyz_io_test03 ( );
  xyz_io_test04 ( );
  xyz_io_test05 ( );
  xyz_io_test06 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'XYZ_IO_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function [newF, logQaa] = sampleFfromRestrictedSet( aa, anchorObjIDs, LL, propFeatIDs, propPiIN, myFPr, TargetPsi )
% Sample feature assignments to the feature pair indicated by propFeatIDs

switch aa
    case anchorObjIDs(1)
        restrictedFeatureSet = [1 0; 1 1];
    case anchorObjIDs(2)
        restrictedFeatureSet = [0 1; 1 1];
    otherwise
        restrictedFeatureSet = [0 1; 1 0; 1 1];
end

availFeatIDs = propPiIN.availFeatIDs;
kiloc = find( availFeatIDs == propFeatIDs(1) );
kjloc = find( availFeatIDs == propFeatIDs(2) );

logPr = -Inf(1,3);
for featureVals = restrictedFeatureSet'
    featureVals = featureVals'; % make it row vector again
    featureCombo = (featureVals * [10;1] );

    switch featureCombo
        case 01
            % propPi = TS.obj(ii) with featIDs(1) deleted
            %propPi = buildTransStructForProposal( propTS, hyperparams, model, aa, propFeatIDs(1), 0);
            %propPi = propPi.obj(aa);
            %propPi = propPiIN;
            keepIDs = [1:kiloc-1 kiloc+1:length(availFeatIDs)];
        case 10
            keepIDs = [1:kjloc-1 kjloc+1:length(availFeatIDs)];
            % propPi = TS.obj(ii) with featIDs(2) deleted
            %propPi = buildTransStructForProposal( propTS, hyperparams, model, aa, propFeatIDs(2), 0);
            %propPi = propPi.obj(aa);
        case 11
            keepIDs = 1:length(availFeatIDs);
            %propPi = propPiIN; % propTS.obj(aa);
    end
    propPi.availFeatIDs = propPiIN.availFeatIDs( keepIDs );
    propPi.pi_z = propPiIN.pi_z( keepIDs, keepIDs );
    propPi.pi_init = propPiIN.pi_init( keepIDs );

    logPrY = calcLogMargPrObsSeqFAST(  LL( propPi.availFeatIDs, :), propPi.pi_z );
    switch aa
        case anchorObjIDs(1)
            logPrF = featureVals(2) * log( myFPr(2) ) + (1-featureVals(2)).*log( 1-myFPr(2) );
        case anchorObjIDs(2)
            logPrF = featureVals(1) * log( myFPr(1) ) + (1-featureVals(1)).*log( 1-myFPr(1) );
        otherwise
            logPrF = sum( featureVals .* log( myFPr ) + (1-featureVals).*log( 1-myFPr ) );
    end

    switch featureCombo
        case 01
            logPr(1) = logPrY + logPrF;
        case 10
            logPr(2) = logPrY + logPrF;
        case 11
            logPr(3) = logPrY + logPrF;
    end
end
ps = exp( logPr - max(logPr)  );
if isempty( TargetPsi ) || ~exist('TargetPsi', 'var')
    [choice] = multinomial_single_draw( ps );
else
    targetVals = TargetPsi.F( aa,  TargetPsi.activeFeatIDs );
    switch ( targetVals * [10;1] )
        case 01
            choice = 1;
        case 10
            choice = 2;
        case 11
            choice = 3;
    end
end
switch choice
    case 1
        newF = [0 1];
    case 2
        newF = [1 0];
    case 3
        newF = [1 1];
end

logQaa = log( ps(choice)/sum(ps)  );

%assert( ~isinf( logQaa ), 'Bad Calculation of Proposal Prob.' );
%if isinf( logQaa )
%    fprintf( 'WARNING: Forced transition highly improbable (-Inf)\n' );
%end

author: michaelchughes

repo: NPBayesHMM


text:
function varargout = bounds(x,lower,upper,aux)
%BOUNDS Adds implicit bounds on variables.
%
% BOUNDS IS OBSOLETE: Use standard constraints
%
% BOUNDS(x,Lower,Upper)   Adds bound constraints on variables.
%                         These bounds are used when performing
%                         big M formulations and similar things.

variables = getvariables(x);
if nargin == 1
    lb = yalmip('getbounds',variables);
    lower = lb(:,1);
    upper = lb(:,2);
else
    lower = lower(:);
    upper = upper(:);
    if length(lower)==1
        lower = repmat(lower,length(variables),1);
    end
    if length(upper)==1
        upper = repmat(upper,length(variables),1);
    end
    % 0 - No problems
    % 1 - Trying to bound nonlinear variable
    % 2 - Trying to bound nonlinear operator
    fail = yalmip('setbounds',variables,lower,upper);
    switch fail
        case {1,2}
            error('BOUNDS can only be applied to linear unitary SDPVAR variables.')
        otherwise
    end
end

if nargout>0
    varargout{1} = lower;
    varargout{2} = upper;
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function  [mori,hkl,omega,sigma] = CSL(sigma,CS,varargin)
% coincidence site lattice misorientations for cubic symmetry
%
% Syntax
%   q = CSL(sigma,CS)
%
% Input
%  sigma - order of coincidence site lattice misorientation
%  CS - @crystalSymmetry
%
% Options
%  delta    - search radius around angle or axis
%  maxsigma -
%
% Output
%  o - @orientation
%

if nargin < 2 || ~isa(CS,'crystalSymmetry')
  error('Starting with MTEX 4.2 the second argument to CSL should be crystal symmetry.')
end

% we generate first all CSL misorientations up to sigma = 60
csl = generateCubicCSL(varargin{:});
% and then select those we are interested in
ndx = find( ismember([csl.sigma],sigma));

mori = orientation(CS,CS);
hkl = zeros(0,3);
omega = []; sigma = [];

for k = ndx
  hkl(end+1,:) = csl(k).axis;
  omega(end+1) = csl(k).angle;
  sigma(end+1) = csl(k).sigma;
  mori(end+1) = orientation.byAxisAngle(vector3d(hkl(end,:)),omega(end),CS,CS);
end

[mori,id] = unique(mori);
sigma = sigma(id);

end

function csl = generateCubicCSL(varargin)  % only cubic

csl = struct('sigma',{},'axis',{},'angle',{});

maxsigma = get_option(varargin,'maxsigma',60);

% heuristic
for u=0:5
  for v=0:5
    for w=0:5
      maxis = [u v w];
      %if all(maxis==0), continue; end

      % make sure u,v,w have no common divisor
      if gcd(u,gcd(v,w))~=1, continue; end

      for m=1:10

        N = sum(mod([maxis m],2));

        if mod(N,2)
          alpha = 1;
        else
          alpha = N;
        end

        % compute the energie
        sigma = sum([maxis m].^2)/alpha;

        % compute rotational angle
        omega =  2*atan( sqrt( sum(maxis.^2) )/m );

        % store values
        if omega < 63*degree && sigma < maxsigma
          csl(end+1).angle = omega;
          csl(end).axis = maxis;
          csl(end).sigma = sigma;
        end
      end
    end
  end
end

end


text:
function sparse_grid_mixed_weight_test ( dim_num, level_max_min, ...
  level_max_max, rule, alpha, beta, tol )

%*****************************************************************************80
%
%% SPARSE_GRID_MIXED_WEIGHT_TEST checks the sum of the quadrature weights.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 March 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX_MIN, LEVEL_MAX_MAX, the minimum and
%    maximum values of LEVEL_MAX.
%
%    Input, integer RULE(DIM_NUM), the rule in each dimension.
%     1, "CC",  Clenshaw Curtis, Closed Fully Nested rule.
%     2, "F2",  Fejer Type 2, Open Fully Nested rule.
%     3, "GP",  Gauss Patterson, Open Fully Nested rule.
%     4, "GL",  Gauss Legendre, Open Weakly Nested rule.
%     5, "GH",  Gauss Hermite, Open Weakly Nested rule.
%     6, "GGH", Generalized Gauss Hermite, Open Weakly Nested rule.
%     7, "LG",  Gauss Laguerre, Open Non Nested rule.
%     8, "GLG", Generalized Gauss Laguerre, Open Non Nested rule.
%     9, "GJ",  Gauss Jacobi, Open Non Nested rule.
%    10, "GW",  Golub Welsch, (presumed) Open Non Nested rule.
%    11, "CC_SE", Clenshaw Curtis Slow Exponential, Closed Fully Nested rule.
%    12, "F2_SE", Fejer Type 2 Slow Exponential, Closed Fully Nested rule.
%    13, "GP_SE", Gauss Patterson Slow Exponential, Closed Fully Nested rule.
%    14, "CC_ME", Clenshaw Curtis Moderate Exponential, Closed Fully Nested rule.
%    15, "F2_ME", Fejer Type 2 Moderate Exponential, Closed Fully Nested rule.
%    16, "GP_ME", Gauss Patterson Moderate Exponential, Closed Fully Nested rule.
%    17, "CCN", Clenshaw Curtis Nested, Linear, Closed Fully Nested rule.
%
%    Input, real ALPHA(DIM_NUM), BETA(DIM_NUM), parameters used for
%    Generalized Gauss Hermite, Generalized Gauss Laguerre, and Gauss Jacobi rules.
%
%    Input, real TOL, the tolerance for point equality.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_MIXED_WEIGHT_TEST:\n' );
  fprintf ( 1, '  Compute the weights of a sparse grid.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Each sparse grid is of spatial dimension DIM_NUM,\n' );
  fprintf ( 1, '  and is made up of product grids of levels up to LEVEL_MAX.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, ' Dimension      Rule     Alpha          Beta\n' );
  fprintf ( 1, '\n' );

  for dim = 1 : dim_num
    fprintf ( 1, '  %8d  %8d', dim, rule(dim) );
    if ( rule(dim) == 6 | rule(dim) == 8 | rule(dim) == 9 )
      fprintf ( 1, '  %12e', alpha(dim) );
    end
    if ( rule(dim) == 9 )
      fprintf ( 1, '  %12e', beta(dim) );
    end
    fprintf ( 1, '\n' );
  end

  weight_sum_exact = 1.0;

  for dim = 1 : dim_num

    if ( rule(dim) == 1 )
      weight_sum_exact = weight_sum_exact * 2.0;
    elseif ( rule(dim) == 2 )
      weight_sum_exact = weight_sum_exact * 2.0;
    elseif ( rule(dim) == 3 )
      weight_sum_exact = weight_sum_exact * 2.0;
    elseif ( rule(dim) == 4 )
      weight_sum_exact = weight_sum_exact * 2.0;
    elseif ( rule(dim) == 5 )
      weight_sum_exact = weight_sum_exact * sqrt ( pi );
    elseif ( rule(dim) == 6 )
      weight_sum_exact = weight_sum_exact * gamma ( 0.5 * ( alpha(dim) + 1.0 ) );
    elseif ( rule(dim) == 7 )
      weight_sum_exact = weight_sum_exact * 1.0;
    elseif ( rule(dim) == 8 )
      weight_sum_exact = weight_sum_exact * gamma ( alpha(dim) + 1.0 );
    elseif ( rule(dim) == 9 )
      arg1 = - alpha(dim);
      arg2 = 1.0;
      arg3 = beta(dim) + 2.0;
      arg4 = - 1.0;
      value1 = r8_hyper_2f1 ( arg1, arg2, arg3, arg4 );
      arg1 = - beta(dim);
      arg2 = 1.0;
      arg3 = alpha(dim) + 2.0;
      arg4 = - 1.0;
      value2 = r8_hyper_2f1 ( arg1, arg2, arg3, arg4 );
      weight_sum_exact = weight_sum_exact * ( ...
        value1 / ( beta(dim) + 1.0 ) + value2 / ( alpha(dim) + 1.0 ) );
    elseif ( rule(dim) == 10 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SPARSE_GRID_MIXED_WEIGHT_TEST - Fatal error!\n' );
      fprintf ( 1, '  Do not know weight sum for rule 10.\n' );
      error ( 'SPARSE_GRID_MIXED_WEIGHT_TEST - Fatal error!' );
    elseif ( rule(dim) == 11 )
      weight_sum_exact = weight_sum_exact * 2.0;
    elseif ( rule(dim) == 12 )
      weight_sum_exact = weight_sum_exact * 2.0;
    elseif ( rule(dim) == 13 )
      weight_sum_exact = weight_sum_exact * 2.0;
    elseif ( rule(dim) == 14 )
      weight_sum_exact = weight_sum_exact * 2.0;
    elseif ( rule(dim) == 15 )
      weight_sum_exact = weight_sum_exact * 2.0;
    elseif ( rule(dim) == 16 )
      weight_sum_exact = weight_sum_exact * 2.0;
    elseif ( rule(dim) == 17 )
      weight_sum_exact = weight_sum_exact * 2.0;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SPARSE_GRID_MIXED_WEIGHT_TEST - Fatal error!\n' );
      fprintf ( 1, '  Unexpected value of RULE = %d\n', rule(dim) );
      error ( 'SPARSE_GRID_MIXED_WEIGHT_TEST - Fatal error!' );
    end

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  As a simple test, sum these weights.\n' );
  fprintf ( 1, '  They should sum to exactly %f\n', weight_sum_exact );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     Level      Weight sum  Expected sum    Difference\n' );
  fprintf ( 1, '\n' );

  for level_max = level_max_min : level_max_max

    point_total_num = sparse_grid_mixed_size_total ( dim_num, level_max, rule );

    point_num  = sparse_grid_mixed_size ( dim_num, level_max, rule, alpha, beta, tol );

    sparse_unique_index = sparse_grid_mixed_unique_index ( ...
      dim_num, level_max, rule, alpha, beta, tol, point_num, point_total_num );

    sparse_weight = sparse_grid_mixed_weight ( dim_num, level_max, rule, ...
      alpha, beta, point_num, point_total_num, sparse_unique_index );

    weight_sum = sum ( sparse_weight(1:point_num) );

    weight_sum_error = abs ( weight_sum - weight_sum_exact );

    fprintf ( 1, '  %8d  %14e  %14e  %14e\n', ...
      level_max, weight_sum, weight_sum_exact, weight_sum_error );

  end

  return
end


text:
function [varargout]=febioStruct2xml(varargin)

% function [domNode]=febioStruct2xml(parseStruct,fileName,optionStruct)
% ------------------------------------------------------------------------
%
% This function provides the basis for coding FEBio XML input file content
% and helps to generate the XML input file.
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2018/01/15 Created
% 2018/02/09 Found bug in text write mode.
% 2018/02/18 Created several "array parse methods" including based on text
% file exporting providing a significant computational time reduction over
% full XML based parsing.
% 2018/05/15 Create temp directory if it is does not exist
% To do:
% Gracefully handle empty fields e.g. 0x8 element array
% 2020/11/26 Fixed bug in relation to undefined attributeStruct in
% febioStruct2xmlStep
% For febio_spec 4.0 perhaps we can remove 'point' from the arrayLoopKeywords
% Check for for febio_spec 4.0
%------------------------------------------------------------------------

%% Parse input

defaultOptionStruct.attributeKeyword='ATTR';
defaultOptionStruct.valueKeyword='VAL';
defaultOptionStruct.arrayLoopKeywords={'node','elem','face','delem','quad4','quad8','tri3','tri6','tri7','line2','line3','point','pt'};
defaultOptionStruct.arrayRowWrapKeywords={'NodeSet','node_data','element_data'};
defaultOptionStruct.rowWrapLength=8;

%Create number to string conversion option structure (used for mat2strIntDouble)
optionStructConv.formatDouble='%6.7e';
optionStructConv.formatInteger='%d';
optionStructConv.dlmChar=',';
optionStructConv.rowWrapLength=[]; %Is altered in code depending on filedname with arrayRowWrapKeywords
defaultOptionStruct.optionStructConv=optionStructConv; %Add to default option structure

if contains(lower(getFEBioPath),'febio2')
    warning('febio2 and related febio_spec handling are depricated. Update your codes for the latest FEBio version to avoid future errors.')
    defaultOptionStruct.fieldOrder={'Module','Control','Globals','Material',...
        'Geometry','MeshData','Initial','Boundary',...
        'Loads','Contact','Constraints','Rigid Connectors',...
        'Discrete','LoadData','Output','Parameters'};
elseif contains(lower(getFEBioPath),'febio3')
    warning('febio3 and related febio_spec handling are depricated. Update your codes for the latest FEBio version to avoid future errors.')
    defaultOptionStruct.fieldOrder={...
        'Module','Control','Globals','Material',...
        'Mesh','MeshDomains','MeshData','Initial','Boundary',...
        'Rigid','Loads','Contact','Constraints',...
        'Discrete','LoadData','Step','Output'};
else
    defaultOptionStruct.fieldOrder={...
        'Module','Control','Globals','Material',...
        'Mesh','MeshDomains','MeshData','Initial','Boundary',...
        'Rigid','Loads','Contact','Constraints',...
        'Discrete','LoadData','Step','Output'};
end

defaultOptionStruct.arrayParseMethod=1;

switch nargin
    case 1
        parseStruct=varargin{1};
        fileName=[];
        optionStruct=[];
    case 2
        parseStruct=varargin{1};
        fileName=varargin{2};
        optionStruct=[];
    case 3
        parseStruct=varargin{1};
        fileName=varargin{2};
        optionStruct=varargin{3};
    otherwise
        error('Wrong number of intput arguments, between 1 and 3 input arguments supported i.e.: parseStruct,fileName,optionStruct')
end

[optionStruct]=structComplete(optionStruct,defaultOptionStruct,1); %Complement provided with default if missing or empty

defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');

%Create temp folder if it does not exist
if ~exist(savePath,'file')
    mkdir(savePath);
end

%Create output file name if not provided
if isempty(fileName)
    fileName=fullfile(savePath,'temp.xml');
end

%% Check febio_spec version
if isfield(parseStruct,'ATTR')
    if isfield(parseStruct.ATTR,'version')
        if contains(parseStruct.ATTR.version,'2.5')
            warning('Old febio_spec detected. febio_spec 2.5 is depricated. Please upgrade to febio_spec 3.0');
            if contains(lower(getFEBioPath),'febio3')
                warning('febio_spec 2.5 is likely not compatible with FEBio3');
            end
        end
    end
end

%% Initialize XML object

domNode = com.mathworks.xml.XMLUtils.createDocument('febio_spec'); %Create the overall febio_spec field

%% Add comment
commentString = ['Created using GIBBON, ',char(datetime)];
commentNode = domNode.createComment(commentString);
rootNode=domNode.getElementsByTagName('febio_spec').item(0);
rootNode.appendChild(commentNode);

%% Reorder fields
[parseStruct]=reorderStruct(parseStruct,optionStruct.fieldOrder);

%% Add content to XML
[domNode]=febioStruct2xmlStep(domNode,rootNode,parseStruct,optionStruct,fileName);

%% Export XML file
xmlwrite_xerces(fileName,domNode); %Custom XML write function

%Import to cut out bad element names from text write mode
T=txtfile2cell(fileName);

logicKeep=~gcontains(T,'xxxxxxxxxxx_CUT');
cell2txtfile(fileName,T(logicKeep),0,0);

%% Collect output
if nargout>0
    varargout{1}=domNode;
end

end

function [domNode,rootNode]=febioStruct2xmlStep(domNode,rootNode,parseStruct,optionStruct,fileName)

arrayParseMethod=optionStruct.arrayParseMethod;

%%
%Get keywords
attributeKeyword=optionStruct.attributeKeyword;
valueKeyword=optionStruct.valueKeyword;
arrayLoopKeywords=optionStruct.arrayLoopKeywords;

%Check if parseStruct if empty
if isempty(parseStruct) %throw error
    error('Empty structure component encountered. Note this may occur when index is skipped in febio structure')
end

%Get field names occuring in the input structure
fieldNameSet = fieldnames(parseStruct);
attributeStruct=[]; %Initialize as empty
for q_field=1:1:numel(fieldNameSet) %Loop for all field names

    %Get current field names and values
    currentFieldName=fieldNameSet{q_field}; %attribute name or element name
    currentFieldValue=parseStruct.(fieldNameSet{q_field}); %attribute value or element value

    if ~any(strcmp(currentFieldName,{attributeKeyword,valueKeyword}))
        if any(strcmp(currentFieldName,optionStruct.arrayRowWrapKeywords))
            optionStruct.optionStructConv.rowWrapLength=optionStruct.rowWrapLength;
        else
            optionStruct.optionStructConv.rowWrapLength=[];
        end
    end

    %Switch behaviour depending on name, i.e. create attribute, assign
    %value (to current entity), or create an element with a value
    switch currentFieldName
        case attributeKeyword %ATTRIBUTE
            attributeStruct=currentFieldValue; %The attribute value should be a structure e.g. ATTR.attributeName.attributeValue
            if ~isempty(rootNode)
                %Add the current set of attribute for the element at rootNode
                [domNode]=addAttributeSetXML(domNode,rootNode,attributeStruct);
            else
                error('The rootNode is empty, the current attribute appears to lack a parent element');
            end
        case valueKeyword %VALUE
            %Add the current value to rootNode
            [domNode]=addElementValueXML(domNode,rootNode,currentFieldValue,optionStruct.optionStructConv);
        otherwise %ELEMENT

            %Check if the current element tag correspods to one requiring looping
            logicArray=any(strcmp(currentFieldName,arrayLoopKeywords));

            if logicArray %Check for looping across values
                parseStructSub=currentFieldValue; %The current element should be a structure
                if isfield(parseStructSub,valueKeyword) %if the structure contains a fielname which mathes the value keyword
                    arrayDataSet=parseStructSub.(valueKeyword); %Get the array of data
                    parseStructSub = rmfield(parseStructSub,valueKeyword); %Remove the value field from the array
                    numArray=size(arrayDataSet,1); %Measure the number of entries in the array
                else %Some other field
                    arrayDataSet=[];
                    numArray=[];
                end

                if isfield(parseStructSub,attributeKeyword) %if the structure contains a fielname which mathes the attribute keyword
                    attributeStruct=parseStructSub.(attributeKeyword);
                    parseStructSub = rmfield(parseStructSub,attributeKeyword); %Remove the value field from the array
                    fieldNameSetSub=fieldnames(attributeStruct);
                else
                    fieldNameSetSub=[];
                end

                if ~isempty(arrayDataSet) %if array entries exist
                    switch arrayParseMethod
                        case 0
                            for q_array=1:1:numArray %Loop over elements in array
                                %Add the element and its value
                                if iscell(arrayDataSet)
                                    [domNode,elementNode,rootNode]=addElementXML(domNode,rootNode,currentFieldName,arrayDataSet{q_array},optionStruct.optionStructConv);
                                else
                                    [domNode,elementNode,rootNode]=addElementXML(domNode,rootNode,currentFieldName,arrayDataSet(q_array,:),optionStruct.optionStructConv);
                                end
                                if ~isempty(fieldNameSetSub) %If attributes exist
                                    for q_fieldSub=1:1:numel(fieldNameSetSub) %loop over all attributes
                                        currentFieldNameSub=fieldNameSetSub{q_fieldSub}; %Get current attribute name
                                        %Create structure for attribute creation
                                        if iscell(attributeStruct.(currentFieldNameSub))
                                            attributeStructSub.(currentFieldNameSub)=attributeStruct.(currentFieldNameSub){q_array};
                                        else
                                            attributeStructSub.(currentFieldNameSub)=attributeStruct.(currentFieldNameSub)(q_array,:);
                                        end
                                        [domNode]=addAttributeSetXML(domNode,elementNode,attributeStructSub); %Add attribute
                                    end
                                end
                            end
                        otherwise
                            [domNode,rootNode]=textModeElementAttributeParse(domNode,rootNode,currentFieldName,attributeStruct,arrayDataSet,numArray,fieldNameSetSub,fileName,arrayParseMethod);
                    end

                elseif ~isempty(fieldNameSetSub) %If this is not empty there are attributes to add despite empty array

                    %Get length from first
                    currentFieldNameSub=fieldNameSetSub{1};
                    numArray=numel(attributeStruct.(currentFieldNameSub));

                    switch arrayParseMethod
                        case 0
                            for q_array=1:1:numArray %Loop over all entries
                                [domNode,elementNode,rootNode]=addElementXML(domNode,rootNode,currentFieldName,[]); %Create element without value
                                for q_fieldSub=1:1:numel(fieldNameSetSub) %Loop over and add all attributes
                                    currentFieldNameSub=fieldNameSetSub{q_fieldSub};
                                    %Create structure for attribute creation
                                    if iscell(attributeStruct.(currentFieldNameSub))
                                        attributeStructSub.(currentFieldNameSub)=attributeStruct.(currentFieldNameSub){q_array};
                                    else
                                        attributeStructSub.(currentFieldNameSub)=attributeStruct.(currentFieldNameSub)(q_array,:);
                                    end
                                    [domNode]=addAttributeSetXML(domNode,elementNode,attributeStructSub); %Add attribute
                                end
                            end
                        otherwise
                            [domNode,rootNode]=textModeElementAttributeParse(domNode,rootNode,currentFieldName,attributeStruct,arrayDataSet,numArray,fieldNameSetSub,fileName,arrayParseMethod);
                    end
                elseif ~isempty(parseStructSub) %Despite arrayLoopKeywords trigger no values or attributes are to be added, but other fields might exist
                    [domNode,rootNode]=recursiveElementParse(domNode,rootNode,currentFieldName,currentFieldValue,optionStruct,fileName);
                end
            else
                [domNode,rootNode]=recursiveElementParse(domNode,rootNode,currentFieldName,currentFieldValue,optionStruct,fileName);
            end
    end
end

end

function [domNode,rootNode]=recursiveElementParse(domNode,rootNode,currentFieldName,currentFieldValue,optionStruct,fileName)
if any(strcmp(currentFieldName,optionStruct.arrayRowWrapKeywords))
    optionStruct.optionStructConv.rowWrapLength=optionStruct.rowWrapLength;
else
    optionStruct.optionStructConv.rowWrapLength=[];
end

if isnumeric(currentFieldValue) || ischar(currentFieldValue)
    [domNode,~,rootNode]=addElementXML(domNode,rootNode,currentFieldName,currentFieldValue,optionStruct.optionStructConv);
elseif isstruct(currentFieldValue)
    [domNode,elementNode,rootNode]=addElementXML(domNode,rootNode,currentFieldName,[],optionStruct.optionStructConv); %Create the current element
    [domNode,elementNode]=febioStruct2xmlStep(domNode,elementNode,currentFieldValue,optionStruct,fileName); %Add whatever is nested in here
elseif iscell(currentFieldValue)
    for q_cell=1:1:numel(currentFieldValue)
        parseStructSub=currentFieldValue{q_cell};
        [domNode,elementNode,rootNode]=addElementXML(domNode,rootNode,currentFieldName,[],optionStruct.optionStructConv);
        [domNode,elementNode]=febioStruct2xmlStep(domNode,elementNode,parseStructSub,optionStruct,fileName); %Add whatever is nested in here
    end
end

end

%%

function [domNode,rootNode]=textModeElementAttributeParse(domNode,rootNode,currentFieldName,attributeStruct,arrayDataSet,numArray,fieldNameSetSub,fileName,arrayParseMethod)

switch arrayParseMethod
    case 1 %Non-flexible, assumes equal attribute type and size (fastest)
        textFormat=['<',currentFieldName,' '];
        if ~isempty(fieldNameSetSub) %If attributes exist
            D=[];
            for q_fieldSub=1:1:numel(fieldNameSetSub) %loop over all attributes
                attributeName=fieldNameSetSub{q_fieldSub}; %Get current attribute name
                if iscell(attributeStruct.(attributeName))
                    attributeNameValue=attributeStruct.(attributeName);
                else
                    attributeNameValue=attributeStruct.(attributeName);
                end
                D=[D attributeNameValue];
                if all(isrounded(attributeNameValue))
                    t_form=repmat(['%d',', '],1,size(attributeNameValue,2));
                else
                    t_form=repmat(['%6.7e',', '],1,size(attributeNameValue,2));
                end
                t_form=t_form(1:end-2); %Take away last space and comma
                textFormat=[textFormat,' ',attributeName,'="',t_form,'"'];
            end
        else
            D=[];
        end

        if all(isrounded(arrayDataSet))
            t_form=repmat(['%d',', '],1,size(arrayDataSet,2));
        else
            t_form=repmat(['%6.7e',', '],1,size(arrayDataSet,2));
        end
        t_form=t_form(1:end-2); %Take away last space and comma
        textFormat=[textFormat,'>',t_form,'</',currentFieldName,'>','\n'];
        D=[D arrayDataSet];

        %Convert to string
        t=sprintf(textFormat,D');
        t=t(1:end-1); %take away last end of line character \n, would result in empty line in text file
        T_add={t};
    case 2 %Flexible, handles non-uniform attributes, uses cellfun to "loop"
        T_add=cell(numArray,1);
        tquote=repmat({'"'},numArray,1);
        if ~isempty(fieldNameSetSub) %If attributes exist
            for q_fieldSub=1:1:numel(fieldNameSetSub) %loop over all attributes
                attributeName=fieldNameSetSub{q_fieldSub}; %Get current attribute name
                attributeNameValue=attributeStruct.(attributeName);
                T_attributeValue=var2cellstr(attributeNameValue);
                T_add=strcat(T_add,repmat({[' ',attributeName,'=']},numArray,1),tquote,T_attributeValue,tquote);
            end
        end
        T_val=var2cellstr(arrayDataSet);
        T_add=strcat(repmat({['<',currentFieldName,' ']},numArray,1),T_add,repmat({'>'},numArray,1),T_val,repmat({['</',currentFieldName,'>']},numArray,1));
    case 3 %Flexible, handles non-uniform attributes, loop through all entries (slowest)
        T_add=cell(numArray,1);
        for q_array=1:1:numArray %Loop over elements in array
            attributeText=[];
            if ~isempty(fieldNameSetSub) %If attributes exist
                for q_fieldSub=1:1:numel(fieldNameSetSub) %loop over all attributes
                    attributeName=fieldNameSetSub{q_fieldSub}; %Get current attribute name
                    if iscell(attributeStruct.(attributeName))
                        attributeNameValue=attributeStruct.(attributeName){q_array};
                    else
                        attributeNameValue=attributeStruct.(attributeName)(q_array,:);
                    end
                    [attributeNameValue]=vec2strIntDouble(attributeNameValue,'%6.7e');
                    attributeText=[attributeText,' ',attributeName,'="',attributeNameValue,'"'];
                end
            end
            if ~isempty(arrayDataSet)
                [currentValue]=vec2strIntDouble(arrayDataSet(q_array,:),'%6.7e');
            else
                currentValue=[];
            end
            T_add{q_array}=['<',currentFieldName,' ',attributeText,'>',currentValue,'</',currentFieldName,'>'];
        end
end

fieldNameCut=[currentFieldName,'_xxxxxxxxxxx_CUT'];
T_top{1}=['<',fieldNameCut,'>'];
totalTextCell=T_top;
totalTextCell(end+1:end+numel(T_add))=T_add;
totalTextCell(end+1)={['</',fieldNameCut,'>']};

% Export text cell to .feb file
cell2txtfile(fileName,totalTextCell,0,0); %Export to text file

% Reimport XML type
domNode2 = xmlread(fileName); %Reimport docNode

nodeList=domNode2.getElementsByTagName(fieldNameCut);
rootNode2=nodeList.item(0);

firstDocImportedNode = domNode.importNode(rootNode2.getChildNodes,true);
rootNode.appendChild(firstDocImportedNode);

%Delete file
delete(fileName);

end

%%
function [B]=reorderStruct(A,fieldOrder)

%Start with fieldOrder fields
fieldNameSet=fieldnames(A);
for q=1:1:numel(fieldOrder)
    fieldNameNow=fieldOrder{q};
    if any(strcmp(fieldNameNow,fieldNameSet))
        B.(fieldNameNow)=A.(fieldNameNow);
        A=rmfield(A,fieldNameNow);
    end
end

%Add remaining fields
fieldNameSet=fieldnames(A);
for q=1:1:numel(fieldNameSet)
    fieldNameNow=fieldNameSet{q};
    B.(fieldNameNow)=A.(fieldNameNow);
end

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


text:
function [model, lbls] = fgplvmReadFromFID(FID)

% FGPLVMREADFROMFID Load from a FID produced by the C++ implementation.
% FORMAT
% DESC loads in from a file stream the data format produced by the
% C++ implementation of the GP-LVM.
% ARG FID : the file ID from where the data is loaded.
% RETURN model : the model loaded in from the file.
% RETURN lbls : any data labels associated with the data in the
% model.
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006, 2009
%
% SEEALSO : fgplvmReadFromFile

% FGPLVM
try
  version = readVersionFromFID(FID);
catch
  if strcmp(lasterr, 'Incorrect file format')
    version = readDoubleFromFID(FID, 'gplvmVersion');
  else
    error(lasterr)
  end
end

if version > 0.2
  error('Incorrect file version.')
end

if version>0.11
  baseType = readStringFromFID(FID, 'baseType')
  if ~strcmp(baseType, 'datamodel')
    error('Incorrect base type in file.')
  end
  type = readStringFromFID(FID, 'type')
  if ~strcmp(type, 'gplvm')
    error('Incorrect type in file.')
  end
end
numData = readIntFromFID(FID, 'numData')
if version > 0.11
  strv = 'outputDim';
else
  strv = 'numProcesses';
end
dataDim = readIntFromFID(FID, strv)
if version > 0.11
  strv = 'inputDim';
else
  strv = 'latentDim';
end
latentDim = readIntFromFID(FID, strv);
if(version>0.1)
  latentRegularised = readBoolFromFID(FID, 'latentRegularised');
  backConstrained = readBoolFromFID(FID, 'backConstrained');
  dynamicsLearnt = readBoolFromFID(FID,  'dynamicsLearnt');
else
  latentRegularised = 0;
  backConstrained = 0;
  dynamicsLearnt = 0;
end
if version > 0.11
  kern = modelReadFromFID(FID);
else
  kern = kernReadFromFID(FID, version);
end
if dynamicsLearnt
  if version > 0.11
    dynKern = modelReadFromFID(FID);
  else
    dynKern = kernReadFromFID(FID, version);
  end
end
if version > 0.11
  noise = modelReadFromFID(FID);
else
  noise = noiseReadFromFID(FID, version);
end

labelsPresent = 0;
labels=ones(numData, 1);
lineStr = getline(FID);
tokens = tokenise(lineStr,',');
for i=1:length(tokens)
  subTokens = tokenise(tokens{i}, ':');
  if strcmp(subTokens{1}, 'Y')
    if str2num(subTokens{2})~=dataDim
      error('Incorrect file format.');
    end
  elseif strcmp(subTokens{1}, 'X')
    if str2num(subTokens{2})~=latentDim
      error('Incorrect file format.');
    end
  elseif strcmp(subTokens{1}, 'labels')
    labelsPresent = 1;
    if str2num(subTokens{2})~=1
      error('Incorrect file format.');
    end
  else
    error('Incorrect file format.');
  end
end

X = zeros(numData, latentDim);
y = zeros(numData, dataDim);

for i=1:numData
  lineStr = getline(FID);
  tokens = tokenise(lineStr);
  for j=1:dataDim
    y(i, j) = str2num(tokens{j});
  end
  for j=1:latentDim
    X(i, j) = str2num(tokens{j+dataDim});
  end
  if labelsPresent
    labels(i) = str2num(tokens{dataDim+latentDim+1});
  end
end
if(~strcmp(noise.type, 'scale'))
  error('Can only load noise of type scale.')
end
if latentRegularised
  prior = 'gaussian';
else
  prior = [];
end
if backConstrained
  warning('Not loading form of back constraints.')
end
options = fgplvmOptions('ftc');
options.kern = kern;
options.prior = prior;
model = fgplvmCreate(latentDim, size(y, 2), y, options);
model.X = X;
model.scale = noise.scale;
model.bias = noise.bias;
model.m = gpComputeM(model);

% This forces kernel computation.
initParams = fgplvmExtractParam(model);
model = fgplvmExpandParam(model, initParams);
if dynamicsLearnt
  options = gpOptions;
  options.kern = dynKern;
  model = fgplvmAddDynamics(model, 'gp', options);
end
lbls=[];
if labelsPresent
  minLbl = min(labels);
  maxLbl = max(labels);
  counter = 0;
  for i=minLbl:maxLbl
    counter = counter + 1;
    lbl = zeros(size(labels));
    lbl(find(labels==i))=1;
    lbls = [lbls lbl];
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function rr = SFu(data, up)
%SFu fuses RR estimates using an implementation of smart fusion.

bw_est = data.bw.v;
am_est = data.am.v;
fm_est = data.fm.v;

%% Find times at which all three RRs were estimated
mutual_times = intersect(intersect(data.bw.t, data.fm.t), data.am.t);
for mod = {'bw', 'am', 'fm'}
    eval(['temp_data = data.' mod{1,1} ';']);
    [rel_times,rel_els,~] = intersect(temp_data.t,mutual_times);

    eval([mod{1,1} '_est = data.' mod{1,1} '.v(rel_els);']);
end

%% Cycle through times at which RR was estimated
% Find out min length incase one has fewer estimates.
min_length = min([length(bw_est), length(am_est), length(fm_est)]);
rr.v = nan(min_length,1);
rr.t = rel_times;
for s = 1 : min_length

    rrEst_values = [bw_est(s), am_est(s), fm_est(s)];

    rel_els = ~isnan(rrEst_values);

    standev = std(rrEst_values);
    % is - 3/22/2019
    %standev = nanstd(rrEst_values);
    % is - end

    if sum(rel_els) < 3 || standev > 4
        rr.v(s) = nan;
    else
        rr.v(s) = mean(rrEst_values);
    end

end

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
%this examples fits a simple neural network to the function
%z= 0.5*sin(pi*y(1,:).^2).*sin(2*pi*y(2,:));

N_neurons=15;%number of neurons
N_in=2;%number of inputs

N_params=N_neurons*(N_in+2)+1;%number of parameters

N_pop_min=N_params+1;%minimum population

N_pop=round(N_pop_min*1.5);%population used

fcn_name='fitness_function';%name of function to maximize

bounds=[-20*ones(1,N_params); 20*ones(1,N_params)];

gen_max=100000;%number of times to evaluate fcn_name;

x_start=[];%let the complex function initialize the population randomly

fit_start=[];%let the complex funciton initialize the fitness;

fcn_opts.N_neurons=N_neurons;%this parameter is passed to the fitness function

%%%%%%%%%%%%crunch numbers!%%%%%%%%%%
tic;
[x_best, fit_best, x_pop, fit_pop stats]=complexmethod(fcn_name,bounds,gen_max,x_start,fit_start,fcn_opts);
timtoc=toc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf('Total time=%f s. Time per generation=%f s\n',timtoc,timtoc/gen_max);
fprintf('Final RMSE=%f\n',-fit_best)

%%%%the rest is plotting results%%%%%%
figure(1)
loglog(-stats.trace_fitness,'.')
xlabel('Generation')
ylabel('RMSE')

N_p=10;%number of points in each dimension
y=[reshape(linspace(0,1,N_p)'*ones(1,N_p),1,[]);reshape((linspace(0,1,N_p)'*ones(1,N_p))',1,[])];%make mesh

N_in=size(y,1);
z= 0.5*sin(pi*y(1,:).^2).*sin(2*pi*y(2,:));

W1=reshape(x_best(1:N_neurons*N_in),N_neurons,N_in);%extract weights for first layer
B1=reshape(x_best((N_neurons*N_in+1):N_neurons*(N_in+1)),N_neurons,1);
W2=reshape(x_best((N_neurons*(N_in+1)+1):N_neurons*(N_in+2)),1,N_neurons);
B2=x_best(N_neurons*(N_in+2)+1);
z_hat=W2*(tanh(W1*y+B1*ones(1,size(y,2))))+B2; %neural calculation

figure(2);
clf
surf(linspace(0,1,N_p),linspace(0,1,N_p),reshape(z_hat,N_p,N_p));
shading interp;
hold on
plot3(y(2,:),y(1,:),z,'kx')
hold off

author: Sable

repo: mcbench-benchmarks


text:
function ts = loadTimestamps(ts_dir)

fid = fopen([ts_dir '/timestamps.txt']);
col = textscan(fid,'%s\n',-1,'delimiter',',');
ts = col{1};
fclose(fid);

repo: MSCKF


text:
%ARMASA version 1.9
%==================
%
%An AutoRegressive Moving Average Spectral Analysis toolbox for use
%with Matlab.
%
%Release date: May 20, 2009
%
%The core of ARMASA consists of four programs: ARMASEL, ARMA2PSD,
%ARMA2COR and MODERR. ARMASEL is a unique program, that is able to
%estimate a time series model from a stationary stochastic signal with
%unknown characteristics, without any user intervention. The automatic
%estimation and selection of either an AR, a MA, or an ARMA model is
%based on statistical criteria and defines a model with an expected
%minimum model error. This model can be used to generate an optimal
%power spectral density estimate, using ARMA2PSD, or alternatively, an
%optimal autocorrelation function estimate, using ARMA2COR. MODERR is a
%quality measure that can be used in simulation experiments to verify
%the goodness of fit of the estimated model. Moreover it can be applied
%in type selection applications, where estimated models from signal
%segments are confronted with reference processes.
%
%The core programs rely on a set of additional programs, bundled in the
%package, that can also be used individually.
%A group of programs, bundled in directory
%'\ARMASA\ASA' carrying the prefix ASA, are used to support the main
%programs, without implementing any time series algorithm.
%
%Additional information is provided on ARMASA functions. It is
%summarized in the outline below:
%
%Demonstration
%-------------
%
%  demo_armasa        - Demonstration of some ARMASA toolbox features
%  demo_tutor         - Demonstration of some ARMASA toolbox features
%  demo_simple        - Demonstration of core ARMASA toolbox statements
%
%Main functions
%--------------
%
%Core
%  armasel            - ARMAsel model identification
%  arma2cor           - ARMA parameters to autocorrelations
%  arma2psd           - ARMA parameters to power spectral density
%  arma2pred          - ARMA parameters to prediction
%  moderr             - Model error
%
%Estimation
%  sig2ar             - AR model identification
%  sig2ma             - MA model identification
%  sig2arma           - ARMA model identification
%  ASAglob_subtr_mean - Subtraction of signal mean
%  ASAglob_mean_adj   - Subtraction of the signal mean identifier
%  ASAglob_ar_cond    - ARMA/MA to AR order selection conditioning
%
%Estimator Tools
%  Burg               - Burg type AR estimator
%  Burg_s             - Burg type AR estimator for multiple data segments
%  cic                - CIC finite sample order selection criterion
%
%Parameter Conversion
%  ar2arset           - AR parameters to optimal lower-order AR models
%  rc2arset           - AR reflectioncoefficients to AR models
%  cov2arset          - AR autocovariances to AR models
%
%Signal Processing
%  armafilter         - Digital filter with ARMA filter coefficients
%  simuarma           - Generating data with ARMA filter coefficients
%  convol             - Convolution sum
%  convolrev          - Convolution sum with one vector reversed
%  deconvol           - Deconvolution
%
%Support Functions
%-----------------
%
%  ASAarg             - Input argument arrangement
%  ASAcontrol         - Function control variable
%  ASAglob            - Variables with a global scope
%  ASAglobretr        - ASAglob variable retrieval
%  ASAversionchk      - Function version check
%  ASAversion2numstr  - Version identifier to number as string type
%  ASAversionbkup     - Version-tagged function backup
%  ASAwarn            - Warning string from numbered file
%  ASAerr             - Error string from numbered file


text:
function bbox_reg = trans_bbox_reg(bbox_reg, use_gpu)
% translate bbox_reg model for fast compute
%
% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Shaoqing Ren
%
% This file is part of the SPP code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

    for i = 1:length(bbox_reg.models)
        bbox_reg.models{i}.Beta = bbox_reg.models{i}.Beta';
        if use_gpu
            bbox_reg.models{i}.Beta_gpu = gpuArray(bbox_reg.models{i}.Beta);
        end
    end

    bbox_reg.combined_Beta = cell2mat(cellfun(@(x) x.Beta', bbox_reg.models', 'UniformOutput', false))';
    if use_gpu
        bbox_reg.combined_Beta_gpu = gpuArray(bbox_reg.combined_Beta);
    end
    bbox_reg.combined_T_inv = cell2mat(cellfun(@(x) x.T_inv, bbox_reg.models, 'UniformOutput', false));
end

author: ShaoqingRen


text:
function [vert,conn,tria,tnum] = deltri2(varargin)
%DELTRI2 compute a constrained 2-simplex Delaunay triangula-
%tion in the two-dimensional plane.
%   [VERT,CONN,TRIA,TNUM]=DELTRI2(VERT,CONN,NODE,PSLG,PART)
%   computes the Delaunay trianguation {VERT,TRIA}, the con-
%   straints CONN, and the "inside" status vector TNUM. VERT
%   is an V-by-2 array of XY coordinates to be triangulated,
%   TRIA is a T-by-3 array of vertex indexing, where each
%   row defines a triangle, such that VERT(TRIA(II,1),:),
%   VERT(TRIA(II,2),:) and VERT(TRIA(II,3),:) are the coord-
%   inates of the II-TH triangle. CONN is a C-by-2 array of
%   constraining edges, where each row defines an edge, as
%   per TRIA. The additional arguments NODE,PSLG and PART
%   define a (mutliply-connected) polygonal region, where
%   NODE is an N-by-2 array of vertices and PSLG is a P-by-2
%   array of edges (a piecewise-straight-line-graph), where
%   each row defines an edge as a pair of indices into NODE.
%   PART is a cell-array of polygonal "parts", where each
%   element PART{KK} is an array of edge indices defining a
%   polygonal region. PSLG(PART{KK},:) is the set of edges
%   in the KK-TH part. TNUM is a T-by-1 array of part index-
%   ing, such that TNUM(II) is the index of the part in whi-
%   ch the II-TH triangle resides.
%
%   See also DELAUNAYTRIANGULATION, DELAUNAYTRI, DELAUNAYN

%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 08/07/2018

    vert = []; conn = []; node = []; PSLG = [];
    part = {}; kind = 'constrained';

%---------------------------------------------- extract args
    if (nargin>=+1), vert = varargin{1}; end
    if (nargin>=+2), conn = varargin{2}; end
    if (nargin>=+3), node = varargin{3}; end
    if (nargin>=+4), PSLG = varargin{4}; end
    if (nargin>=+5), part = varargin{5}; end
    if (nargin>=+6), kind = varargin{6}; end

%---------------------------------------------- basic checks
    if (~isnumeric(vert) || ~isnumeric(conn) || ...
        ~isnumeric(node) || ~isnumeric(PSLG) || ...
        ~iscell   (part) || ~ischar   (kind) )
        error('deltri2:incorrectInputClass' , ...
            'Incorrect input class.') ;
    end

    nvrt = size(vert,+1) ; nnod = size(node,+1) ;
    nedg = size(PSLG,+1) ;

%---------------------------------------------- basic checks
    if (ndims(vert) ~= +2 || ndims(conn) ~= +2)
        error('deltri2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end
    if (size(vert,2)~= +2 || size(conn,2)~= +2)
        error('deltri2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

    if (min([conn(:)])<+1 || max([conn(:)])>nvrt)
        error('deltri2:invalidInputs', ...
            'Invalid CONN input array.') ;
    end

%---------------------------------------------- basic checks
    if (nargin >= +3)

    if (ndims(node) ~= +2 || ndims(PSLG) ~= +2)
        error('deltri2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end
    if (size(node,2)~= +2 || size(PSLG,2)~= +2)
        error('deltri2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

    if (min([PSLG(:)])<+1 || max([PSLG(:)])>nnod)
        error('deltri2:invalidInputs', ...
            'Invalid EDGE input array.') ;
    end

    pmin = cellfun(@min,part);
    pmax = cellfun(@max,part);

    if (min([pmin(:)])<+1 || max([pmax(:)])>nedg)
        error('deltri2:invalidInputs', ...
            'Invalid PART input array.') ;
    end

    end

%------------------------------------ compute Delaunay tria.
    switch (lower(kind))
    case 'constrained'

        if (exist( ...
        'delaunayTriangulation') == +2 )
    %-------------------------------- use class if available
        dtri = ...
        delaunayTriangulation(vert,conn) ;
        vert = dtri.Points;
        conn = dtri.Constraints;
        tria = dtri.ConnectivityList;
        else
        if (exist('DelaunayTri') == +2 )
    %-------------------------------- use class if available
        dtri = DelaunayTri   (vert,conn) ;
        vert = dtri.X;
        conn = dtri.Constraints;
        tria = dtri.Triangulation;
        else
    %-------------------------------- *fall-back* onto qhull
       [vert,conn,tria] ...
                = cfmtri2(vert,conn) ;
        end
        end

    case 'conforming'

    %-------------------------------- "conforming" delaunay!
       [vert,conn,tria] ...
                = cfmtri2(vert,conn) ;

    otherwise
        error('deltri2:invalidInputs', ...
            'Invalid KIND selection.') ;

    end

%------------------------------------ calc. "inside" status!
    tnum = zeros(size(tria,+1),+1) ;

    if (nargin >= +3)

    tmid = vert(tria(:,1),:) ...
         + vert(tria(:,2),:) ...
         + vert(tria(:,3),:) ;
    tmid = tmid / +3.0;

    for ppos = 1 : length(part)

       [stat] = inpoly2( ...
            tmid,node  , ...
            PSLG(part{ppos},:))  ;

        tnum(stat)  = ppos ;

    end

%------------------------------------ keep "interior" tria's
    tria = tria(tnum>+0,:) ;
    tnum = tnum(tnum>+0,:) ;

    end

%------------------------------------ flip for correct signs
    area = triarea(vert,tria) ;

    tria(area<0.,:) = ...
        tria(area<0.,[1,3,2]) ;

end

author: dengwirda


text:
function out=mpeuler(precision)

if nargin==0
 mp_defaults
 precision=default_precision;
end

out_rval=mpfr_euler(precision);

out=class(struct('rval',out_rval,...
                  'ival','0',...
                  'precision',precision),'mp');

author: opencobra

repo: cobratoolbox


text:
%MAKEFUZZYEDGESYS Script to make MAT-file used by FUZZYFILT.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Input membership functions.
zero = @(z) bellmf(z,-0.3,0);
not_used = @(z) onemf(z);

% Output membership functions.
black = @(z) triangmf(z,0,0,0.75);
white = @(z) triangmf(z,0.25,1,1);

% There are four rules and four inputs, so the inmf matrix is 4x4. Each
% row of the inmf matrix corresponds to one rule.
inmf = {zero, not_used, zero, not_used
   not_used, not_used, zero, zero
   not_used, zero, not_used, zero
   zero, zero, not_used, not_used};

% Specify cell array of output membership functions, OUTMF (see function
% IMPLFCNS). There are four IF-THEN rules in this problem, all resulting
% in WHITE, and one ELSE rule resulting in BLACK (see Fig. 3.40).
outmf = {white,white,white,white,black};

% Inputs to the output membership functions are in the range [0,1].
vrange = [0,1];

F = fuzzysysfcn(inmf,outmf,vrange);

% Compute a lookup-table-based approximation to the fuzzy system
% function. Each of the four inputs is in the range [-1,1].
G = approxfcn(F,[-1 1; -1 1; -1 1; -1 1]);

% Save the fuzzy system approximation function to a MAT-file.
save fuzzyedgesys G

author: dipum

repo: dipum-toolbox

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


text:
%points2contour
%Tristan Ursell
%Sept 2013
%
%[Xout,Yout]=points2contour(Xin,Yin,P,direction)
%[Xout,Yout]=points2contour(Xin,Yin,P,direction,dlim)
%[Xout,Yout,orphans]=points2contour(Xin,Yin,P,direction,dlim)
%[Xout,Yout,orphans,indout]=points2contour(Xin,Yin,P,direction,dlim)
%
%Given any list of 2D points (Xin,Yin), construct a singly connected
%nearest-neighbor path in either the 'cw' or 'ccw' directions.  The code
%has been written to handle square and hexagon grid points, as well as any
%non-grid arrangement of points.
%
%'P' sets the point to begin looking for the contour from the original
%ordering of (Xin,Yin), and 'direction' sets the direction of the contour,
%with options 'cw' and 'ccw', specifying clockwise and counter-clockwise,
%respectively.
%
%The optional input parameter 'dlim' sets a distance limit, if the distance
%between a point and all other points is greater than or equal to 'dlim',
%the point is left out of the contour.
%
%The optional output 'orphans' gives the indices of the original (Xin,Yin)
%points that were not included in the contour.
%
%The optional output 'indout' is the order of indices that produces
%Xin(indout)=Xout and Yin(indout)=Yout.
%
%There are many (Inf) situations where there is no unique mapping of points
%into a connected contour -- e.g. any time there are more than 2 nearest
%neighbor points, or in situations where the nearest neighbor matrix is
%non-symmetric.  Picking a different P will result in a different contour.
%Likewise, in cases where one point is far from its neighbors, it may be
%orphaned, and only connected into the path at the end, giving strange
%results.
%
%The input points can be of any numerical class.
%
%Note that this will *not* necessarily form the shortest path between all
%the points -- that is the NP-Hard Traveling Salesman Problem, for which
%there is no deterministic solution.  This will, however, find the shortest
%path for points with a symmetric nearest neighbor matrix.
%
%see also: bwtraceboundary
%
%Example 1:  continuous points
%N=200;
%P=1;
%theta=linspace(0,2*pi*(1-1/N),N);
%[~,I]=sort(rand(1,N));
%R=2+sin(5*theta(I))/3;
%
%Xin=R.*cos(theta(I));
%Yin=R.*sin(theta(I));
%
%[Xout,Yout]=points2contour(Xin,Yin,P,'cw');
%
%figure;
%hold on
%plot(Xin,Yin,'b-')
%plot(Xout,Yout,'r-','Linewidth',2)
%plot(Xout(2:end-1),Yout(2:end-1),'k.','Markersize',15)
%plot(Xout(1),Yout(1),'g.','Markersize',15)
%plot(Xout(end),Yout(end),'r.','Markersize',15)
%xlabel('X')
%ylabel('Y')
%axis equal tight
%title(['Black = original points, Blue = original ordering, Red = new ordering, Green = starting points'])
%box on
%
%
%Example 2:  square grid
%P=1;
%
%Xin=[1,2,3,4,4,4,4,3,2,1,1,1];
%Yin=[0,0,0,0,1,2,3,3,2,2,1,0];
%
%[Xout,Yout]=points2contour(Xin,Yin,P,'cw');
%
%figure;
%hold on
%plot(Xin,Yin,'b-')
%plot(Xout,Yout,'r-','Linewidth',2)
%plot(Xout(2:end-1),Yout(2:end-1),'k.','Markersize',15)
%plot(Xout(1),Yout(1),'g.','Markersize',15)
%plot(Xout(end),Yout(end),'r.','Markersize',15)
%xlabel('X')
%ylabel('Y')
%axis equal tight
%box on
%
%Example 3:  continuous points, pathological case
%N=200;
%P=1;
%theta=linspace(0,2*pi*(1-1/N),N);
%[~,I]=sort(rand(1,N));
%R=2+sin(5*theta(I))/3;
%
%Xin=(1+rand(1,N)/2).*R.*cos(theta(I));
%Yin=(1+rand(1,N)/2).*R.*sin(theta(I));
%
%[Xout,Yout]=points2contour(Xin,Yin,P,'cw');
%
%figure;
%hold on
%plot(Xin,Yin,'b-')
%plot(Xout,Yout,'r-','Linewidth',2)
%plot(Xout(2:end-1),Yout(2:end-1),'k.','Markersize',15)
%plot(Xout(1),Yout(1),'g.','Markersize',15)
%plot(Xout(end),Yout(end),'r.','Markersize',15)
%xlabel('X')
%ylabel('Y')
%axis equal tight
%title(['Black = original points, Blue = original ordering, Red = new ordering, Green = starting points'])
%box on
%
%Example 4:  continuous points, distance limit applied
%N=200;
%P=1;
%theta=linspace(0,2*pi*(1-1/N),N);
%[~,I]=sort(rand(1,N));
%R=2+sin(5*theta(I))/3;
%R(2)=5; %the outlier
%
%Xin=(1+rand(1,N)/16).*R.*cos(theta(I));
%Yin=(1+rand(1,N)/16).*R.*sin(theta(I));
%
%[Xout,Yout,orphans,indout]=points2contour(Xin,Yin,P,'cw',1);
%
%figure;
%hold on
%plot(Xin,Yin,'b-')
%plot(Xin(orphans),Yin(orphans),'kx')
%plot(Xin(indout),Yin(indout),'r-','Linewidth',2)
%plot(Xout(2:end-1),Yout(2:end-1),'k.','Markersize',15)
%plot(Xout(1),Yout(1),'g.','Markersize',15)
%plot(Xout(end),Yout(end),'r.','Markersize',15)
%xlabel('X')
%ylabel('Y')
%axis equal tight
%title(['Black = original points, Blue = original ordering, Red = new ordering, Green = starting points'])
%box on
%

function [Xout,Yout,varargout]=points2contour(Xin,Yin,P,direction,varargin)

%check to make sure the vectors are the same length
if length(Xin)~=length(Yin)
    error('Input vectors must be the same length.')
end

%check to make sure point list is long enough
if length(Xin)<2
    error('The point list must have more than two elements.')
end

%check distance limit
if ~isempty(varargin)
    dlim=varargin{1};
    if dlim<=0
        error('The distance limit parameter must be greater than zero.')
    end
else
    dlim=-1;
end

%check direction input
if and(~strcmp(direction,'cw'),~strcmp(direction,'ccw'))
    error(['Direction input: ' direction ' is not valid, must be either "cw" or "ccw".'])
end

%check to make sure P is in the right range
P=round(P);
npts=length(Xin);

if or(P<1,P>npts)
    error('The starting point P is out of range.')
end

%adjust input vectors for starting point
if size(Xin,1)==1
    Xin=circshift(Xin,[0,1-P]);
    Yin=circshift(Yin,[0,1-P]);
else
    Xin=circshift(Xin,[1-P,0]);
    Yin=circshift(Yin,[1-P,0]);
end

%find distances between all points
D=zeros(npts,npts);
for q1=1:npts
    D(q1,:)=sqrt((Xin(q1)-Xin).^2+(Yin(q1)-Yin).^2);
end

%max distance
maxD=max(D(:));

%avoid self-connections
D=D+eye(npts)*maxD;

%apply distance contraint by removing bad points and starting over
if dlim>0
    D(D>=dlim)=-1;

    %find bad points
    bad_pts=sum(D,1)==-npts;
    orphans=find(bad_pts);

    %check starting point
    if sum(orphans==P)>0
        error('The starting point index is a distance outlier, choose a new starting point.')
    end

    %get good points
    Xin=Xin(~bad_pts);
    Yin=Yin(~bad_pts);

    %number of good points
    npts=length(Xin);

    %find distances between all points
    D=zeros(npts,npts);
    for q1=1:npts
        D(q1,:)=sqrt((Xin(q1)-Xin).^2+(Yin(q1)-Yin).^2);
    end

    %max distance
    maxD=max(D(:));

    %avoid self-connections
    D=D+eye(npts)*maxD;
else
    orphans=[];
    bad_pts=zeros(size(Xin));
end

%tracking vector (has this original index been put into the ordered list?)
track_vec=zeros(1,npts);

%construct directed graph
Xout=zeros(1,npts);
Yout=zeros(1,npts);
indout0=zeros(1,npts);

Xout(1)=Xin(1);
Yout(1)=Yin(1);
indout0(1)=1;

p_now=1;
track_vec(p_now)=1;
for q1=2:npts
    %get current row of distance matrix
    curr_vec=D(p_now,:);

    %remove used points
    curr_vec(track_vec==1)=maxD;

    %find index of closest non-assigned point
    p_temp=find(curr_vec==min(curr_vec),1,'first');

    %reassign point
    Xout(q1)=Xin(p_temp);
    Yout(q1)=Yin(p_temp);

    %move index
    p_now=p_temp;

    %update tracking
    track_vec(p_now)=1;

    %update index vector
    indout0(q1)=p_now;
end

%undo the circshift
temp1=find(~bad_pts);
indout=circshift(temp1(indout0),[P,0]);

%%%%%%% SET CONTOUR DIRECTION %%%%%%%%%%%%
%contour direction is a *global* feature that cannot be determined until
%all the points have been sequentially ordered.

%calculate tangent vectors
tan_vec=zeros(npts,3);
for q1=1:npts
    if q1==npts
        tan_vec(q1,:)=[Xout(1)-Xout(q1),Yout(1)-Yout(q1),0];
        tan_vec(q1,:)=tan_vec(q1,:)/norm(tan_vec(q1,:));
    else
        tan_vec(q1,:)=[Xout(q1+1)-Xout(q1),Yout(q1+1)-Yout(q1),0];
        tan_vec(q1,:)=tan_vec(q1,:)/norm(tan_vec(q1,:));
    end
end

%determine direction of contour
local_cross=zeros(1,npts);
for q1=1:npts
    if q1==npts
        cross1=cross(tan_vec(q1,:),tan_vec(1,:));
    else
        cross1=cross(tan_vec(q1,:),tan_vec(q1+1,:));
    end
    local_cross(q1)=asin(cross1(3));
end

%figure out current direction
if sum(local_cross)<0
    curr_dir='cw';
else
    curr_dir='ccw';
end

%set direction of the contour
if and(strcmp(curr_dir,'cw'),strcmp(direction,'ccw'))
    Xout=fliplr(Xout);
    Yout=fliplr(Yout);
end

%varargout
if nargout==3
    varargout{1}=orphans;
end

if nargout==4
    varargout{1}=orphans;
    varargout{2}=indout;
end

disp('finished')