text:
function [varargout] = webwrite(varargin)

% WEBWRITE is a drop-in replacement for the function with the same
% name that was introduced in MATLAB R2014b. This function is only
% partially compatible with the original.
%
% Use as
%   status = webwrite(url, [data], [options])
%
% This function posts data (which can be a struct or a string) to url.
% The options argument is currently ignored.
%
% Unlike builtin MATLAB webwrite, this function
% - sets headers always to 'Content-Type: application/json; charset=UTF-8'
%   irrespective of the contents of the options argument.
% - the status returned is zero if succesful, and nonzero otherwise
%
% This requires that curl is available on the command-line. If curl
% is not available, or an other error occurs, then this function
% raises an error.

% Copyright (C) 2017-2019, Robert Oostenveld, Nikolaas N. Oosterhof
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% try to automatically remove incorrect compat folders from the path, see https://github.com/fieldtrip/fieldtrip/issues/899

if isempty(strfind(mfilename('fullpath'), matlabroot))
  % automatic cleanup does not work when FieldTrip is installed inside the MATLAB path, see https://github.com/fieldtrip/fieldtrip/issues/1527

  alternatives = which(mfilename, '-all');
  if ~iscell(alternatives)
    % this is needed for octave, see https://github.com/fieldtrip/fieldtrip/pull/1171
    alternatives = {alternatives};
  end

  keep = true(size(alternatives));
  for i=1:numel(alternatives)
    keep(i) = keep(i) && ~any(alternatives{i}=='@');  % exclude methods from classes
    keep(i) = keep(i) && alternatives{i}(end)~='p';   % exclude precompiled files
  end
  alternatives = alternatives(keep);

  if exist(mfilename, 'builtin') || any(strncmp(alternatives, matlabroot, length(matlabroot)) & cellfun(@isempty, strfind(alternatives, fullfile('private', mfilename))))
    % remove this directory from the path
    p = fileparts(mfilename('fullpath'));
    warning('removing "%s" from your path, see http://bit.ly/2SPPjUS', p);
    rmpath(p);
    % call the original MATLAB function
    if exist(mfilename, 'builtin')
      [varargout{1:nargout}] = builtin(mfilename, varargin{:});
    else
      [varargout{1:nargout}] = feval(mfilename, varargin{:});
    end
    return
  end

end % automatic cleanup of compat directories

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is where the actual replacement code starts
% function status = webwrite(url, data, options)

% deal with the input arguments
if nargin==1
  [url               ] = deal(varargin{1:1});
elseif nargin==2
  [url, data         ] = deal(varargin{1:2});
elseif nargin==3
  [url, data, options] = deal(varargin{1:3});
else
  error('incorrect number of input arguments')
end

if nargin<2
  data = '';
end

if isstruct(data)
  % convert it to JSON string
  data = struct2json(data);
  assert(ischar(data));

  to_remove=sprintf('[\n ]'); % newlines and spaces
  data = regexprep(data, to_remove, '');
end

if ~ischar(data)
  error('unsupported data: input must be string or struct');
end

cmd = sprintf(['curl '...
                '-s '...                                % quiet (no output)
                '-H "Content-Type: application/json;'...  % JSON data
                    'charset=UTF-8" '...                % charset
                '-X POST '...                           % method
                '-d ''%s'' %s'], ...                    % data and URL
                data, url);

[status, str] = system(cmd);

if status
  error('An error occured when running this command: ''%s'':\n%s',cmd,str);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% deal with the output arguments

varargout = {status};

author: fieldtrip

repo: fieldtrip


text:
classdef ConstraintEvalTypeEnum < matlab.mixin.SetGet
    %ConstraintEvalTypeEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        FixedBounds('Fixed Bounds')
        StateComparison('State Comparison');
    end

    properties
        name char = '';
    end

    methods
        function obj = ConstraintEvalTypeEnum(name)
            obj.name = name;
        end
    end

    methods(Static)
        function [listBoxStr, enums] = getListBoxStr()
            m = enumeration('ConstraintEvalTypeEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
            enums = m(I);
        end

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

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

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function bounds = computeBounds(config)
% bounds = computeBounds(config)
%
% This function returns the bounds on the decision variables, to be passed
% to the trajectory optimization solver.
%
% INPUTS:
%   config.dyn = physical parameters of the model
%       .m1 = cart mass
%       .m2 = pole mass
%       .g = gravity
%       .l = pendulum length
%
% OUTPUTS:
%   bounds = struct of bounds for decision variables:
%       .initialTime.lower = scalar
%       .initialTime.upper = scalar
%       .finalTime.lower = scalar
%       .finalTime.upper = scalar
%       .initialState.lower = [ns, 1]
%       .initialState.upper = [ns, 1]
%       .finalState.lower = [ns, 1]
%       .finalState.upper = [ns, 1]
%       .state.lower = [ns, 1]
%       .state.upper = [ns, 1]
%       .control.lower = [nc, 1]
%       .control.upper = [nc, 1]
%

% Unpack physical parameters
m1 = config.dyn.m1;
% m2 = config.dyn.m2;
g = config.dyn.g;
l = config.dyn.l;

% Boundry values
bounds.initialState.lower = [0;0;0;0];
bounds.initialState.upper = bounds.initialState.lower;
bounds.finalState.lower = [config.guess.state(1:2,end);0;0];
bounds.finalState.upper = bounds.finalState.lower;

% Bounds on duration (total trajectory time)
wn = sqrt(l/g);
bounds.initialTime.lower = 0;
bounds.initialTime.upper = 0;
bounds.finalTime.lower = 0.5*(2*pi*wn);
bounds.finalTime.upper = 4*(2*pi*wn);

% Bounds on state:
bounds.state.lower = [-1.0*l; -2*pi; -10*l/wn; -10*pi/wn];
bounds.state.upper = [1.0*l; 2*pi; 10*l/wn; 10*pi/wn];

% Bounds on control:
bounds.control.lower = -5*m1*g;
bounds.control.upper = 5*m1*g;

end

author: MatthewPeterKelly

repo: dscTutorials


text:
% Simulation of Simple Pendulum
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Warning : On running this the workspace memory will be deleted. Save if
% any data present before running the code !!
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%-------------------------------------------------------------------------
% Code written by : Siva Srinivas Kolukula                                |
%                   Senior Research Fellow                                |
%                   Structural Mechanics Laboratory                       |
%                   Indira Gandhi Center for Atomic Research              |
%                   India                                                 |
% E-mail : allwayzitzme@gmail.com                                         |
%          http://sites.google.com/site/kolukulasivasrinivas/             |
%-------------------------------------------------------------------------
clear ;clc ;
% Properties of Pendulum (Can be altered)
g = 9.81 ;           % Acceleration due to gravity
M = 10 ;              % Mass of the pendulum
L = 1 ;              % Length of the Pendulum
C = 0.7 ;             % Damping

% Initial Boundary Conditions (Can be altered)
Phi = 0.5;            % Position
dtPhi = 1.0;        % Velocity

duration = 60;      % Duration of the Simulation
fps = 10;           % Frames per second
%movie = true;         % true if wanted to save animation as avi file
movie = false ;         % false if only want to view animation
arrow = 'ShowArrow' ; % Shows the direction of phase plane plot
%arrow = 'NoNoArrow' ;   % Will not show the direction of phase plane plot
interval = [0, duration];               % Time span
ivp=[Phi; dtPhi; g; M; L ; C];     % Initial value's for the problem
% Simulation of Simple Pendulum
Animation(ivp, duration, fps, movie,arrow);

author: Sable

repo: mcbench-benchmarks


text:
function i4row_sort_a_test ( )

%*****************************************************************************80
%
%% I4ROW_SORT_A_TEST tests I4ROW_SORT_A;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  m = 10;
  n = 4;
  b = 0;
  c = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4ROW_SORT_A_TEST\n' );
  fprintf ( 1, '  For a rectangular integer matrix:\n' );
  fprintf ( 1, '  I4ROW_SORT_A sorts the rows;\n' );

  seed = 123456789;

  [ a, seed ] = i4mat_uniform_ab ( m, n, b, c, seed );

  i4mat_print ( m, n, a, '  The original matrix:' );

  a = i4row_sort_a ( m, n, a );

  i4mat_print ( m, n, a, '  The row-sorted matrix:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [logp, yhat, res] = tapas_cdfgaussian_obs(r, infStates, ptrans)
% Calculates the log-probability of response y under a cumulative Gaussian distribution. This
% model has no free parameters.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2015 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Initialize returned log-probabilities as NaNs so that NaN is
% returned for all irregualar trials
n = size(infStates,1);
logp = NaN(n,1);
yhat = NaN(n,1);
res  = NaN(n,1);

% Weed irregular trials out from inferred states and responses
mu2 = infStates(:,2,3);
mu2(r.irr) = [];
sa2 = infStates(:,2,4);
sa2(r.irr) = [];
y = r.y(:,1);
y(r.irr) = [];

% Probability mass for x2 < 0
x2lt0 = 0.5*(1 +erf((0 -mu2)./(sa2.*sqrt(2))));

% Probability of observed choice
probc = y.*(1 -x2lt0) +(1 -y).*x2lt0;

% Calculate log-probabilities for non-irregular trials
% Note: 8*atan(1) == 2*pi (this is used to guard against
% errors resulting from having used pi as a variable).
reg = ~ismember(1:n,r.irr);
logp(reg) = log(probc);
yh = 1 -x2lt0;
yhat(reg) = yh;
res(reg) = (y -yh)./sqrt(yh.*(1 -yh));

return;

author: translationalneuromodeling

repo: tapas


text:
function [ new_images ] = nxel( n, old_image )
%Nxel - takes image and pixelizes into n * n array
    [a,b,~] = size(old_image);
    new_images = cell(n,n);
    for i = 1:n
        for j = 1:n
            new_images{i,j} = old_image((1+(i-1)*floor(a/n)):(floor(i*a/n)-1),...
                                        (1+(j-1)*floor(b/n)):(floor(j*b/n)-1),...
                                   :);
        end
    end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [posR, triR, keeppos] = remove_double_vertices(pos, tri)

% REMOVE_DOUBLE_VERTICES removes double vertices from a triangular mesh
% renumbering the vertex-indices for the triangles and removing all
% triangles with one of the specified vertices.
%
% Use as
%   [pos, tri] = remove_double_vertices(pos, tri)

[pos1,i1,i2] = unique(pos, 'rows');
%keeppos   = find(ismember(pos,pos1,'rows'));
%removepos = setdiff([1:size(pos,1)],keeppos);
keeppos   = i1;
removepos = setdiff(1:size(pos,1), keeppos);

npos = size(pos,1);
ntri = size(tri,1);

if all(removepos==0 | removepos==1)
  removepos = find(removepos);
end

% remove the vertices and determine the new numbering (indices) in numb
keeppos = setdiff(1:npos, removepos);
numb    = zeros(1,npos);
numb(keeppos) = 1:length(keeppos);

% look for triangles referring to removed vertices
removetri = false(ntri,1);
removetri(ismember(tri(:,1), removepos)) = true;
removetri(ismember(tri(:,2), removepos)) = true;
removetri(ismember(tri(:,3), removepos)) = true;

% remove the vertices and triangles
posR = pos(keeppos, :);
triR = tri(~removetri,:);

% renumber the vertex indices for the triangles
triR = numb(triR);

author: spm


text:
function outFiles = out_demultiplex(DataFile, ChannelFile, OutputDir, UseSsp, ram, parallel)
% OUT_DEMULTIPLEX: Load a raw data file and creates separate electrode files.

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

% If the output directory doesn't exist: create it
if ~exist(OutputDir, 'dir')
    mkdir(OutputDir);
end

% Load channel file
ChannelMat = in_bst_channel(ChannelFile);
numChannels = length(ChannelMat.Channel);
% Channel names: Remove any special characters
cleanNames = str_remove_spec_chars({ChannelMat.Channel.Name});

% Assemble output filenames
outFiles = cellfun(@(c)bst_fullfile(OutputDir, ['raw_elec_', c]), cleanNames, 'UniformOutput', 0);
% If all files already exist: nothing else to do in this function
isFileOk = cellfun(@(c)exist([c, '.mat'], 'file'), outFiles);
if all(isFileOk)
    % Add the .mat extension to the file names
    disp(['BST> Channels already demultiplexed in: ' OutputDir]);
    outFiles = cellfun(@(x) [x '.mat'], outFiles, 'UniformOutput', 0);
    return;
% If some files already exist: delete all intermediate existing file, before generating them again
elseif any(isFileOk)
    delete(outFiles{isFileOk});
end

% Load input data file
DataMat = in_bst_data(DataFile, 'F');
sFile = DataMat.F;
sr = sFile.prop.sfreq;

% Apply SSP/ICA when reading from data files
ImportOptions = db_template('ImportOptions');
ImportOptions.UseCtfComp = 0;
ImportOptions.UseSsp     = UseSsp;
% Special case for supported acquisition systems: Save temporary files
% using single precision instead of double to save disk space
if ismember(sFile.format, {'EEG-AXION', 'EEG-BLACKROCK', 'EEG-INTAN', 'EEG-PLEXON'})
    precision = 'single';
    nBytes = 4;
else
    precision = 'double';
    nBytes = 8;
end
ImportOptions.Precision = precision;

% Separate the file to max length based on RAM
max_samples = ram / nBytes / numChannels;
total_samples = round((sFile.prop.times(2) - sFile.prop.times(1)) .* sFile.prop.sfreq); % (Blackrock/Ripple complained). Removed +1
num_segments = ceil(total_samples / max_samples);
num_samples_per_segment = ceil(total_samples / num_segments);

% Loop on segments
for iSegment = 1:num_segments
    sampleBounds(1) = (iSegment - 1) * num_samples_per_segment + round(sFile.prop.times(1)* sFile.prop.sfreq);
    if iSegment < num_segments
        sampleBounds(2) = sampleBounds(1) + num_samples_per_segment - 1;
    else
        sampleBounds(2) = total_samples + round(sFile.prop.times(1)* sFile.prop.sfreq);
    end
    % Read recordings
    F = in_fread(sFile, ChannelMat, [], sampleBounds, [], ImportOptions);
    % Append segment to individual channel file
    if parallel
        bst_progress('start', 'Spike-sorting', 'Demultiplexing raw file...');
        parfor iChannel = 1:numChannels
            electrode_data = F(iChannel,:);
            fid = fopen([outFiles{iChannel} '.bin'], 'a');
            fwrite(fid, electrode_data, precision);
            fclose(fid);
        end
    else
        bst_progress('start', 'Spike-sorting', 'Demultiplexing raw file...', 0, num_segments * numChannels);
        for iChannel = 1:numChannels
            electrode_data = F(iChannel,:);
            fid = fopen([outFiles{iChannel} '.bin'], 'a');
            fwrite(fid, electrode_data, precision);
            fclose(fid);
            bst_progress('inc', 1);
        end
    end
end

% Convert binary files per channel to Matlab files
if parallel
    bst_progress('start', 'Spike-sorting', 'Converting demultiplexed files...');
    parfor iChannel = 1:numChannels
        convert2mat(outFiles{iChannel}, sr, precision);
    end
else
    bst_progress('start', 'Spike-sorting', 'Converting demultiplexed files...', 0, numChannels);
    for iChannel = 1:numChannels
        convert2mat(outFiles{iChannel}, sr, precision);
        bst_progress('inc', 1);
    end
end

% Add the .mat extension to the file names
outFiles = cellfun(@(x) [x '.mat'], outFiles, 'UniformOutput', 0);

end

%% ===== CONVERT BIN TO MAT =====
function convert2mat(chanFile, sr, precision)
    % Read .bin file
    fid = fopen([chanFile '.bin'], 'rb');
    data = fread(fid, precision);
    fclose(fid);
    % Save .mat file
    save([chanFile '.mat'], 'data', 'sr');
    % Delete .bin file
    delete([chanFile '.bin']);
end

author: brainstorm-tools


text:
function calpak_test012 ( )

%*****************************************************************************80
%
%% CALPAK_TEST012 tests JED_TO_YJF_ROMAN and YJF_TO_JED_ROMAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 March 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST012\n' );
  fprintf ( 1, '  For the Roman calendar:\n' );
  fprintf ( 1, '  JED_TO_YJF_ROMAN: JED -> YJF.\n' );
  fprintf ( 1, '  YJF_TO_JED_ROMAN: YJF -> JED.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  JED (in)    YJF                JED (out)\n' );
  fprintf ( 1, '\n' );

  jed_epoch = epoch_to_jed_roman ( );

  i = 0;

  while ( 1 )

    i = i + 1;
    jed1 = jed_test ( i );

    if ( jed1 < 0.0 )
      break
    end

    if ( jed_epoch <= jed1 )

      [ y2, j2, f2 ] = jed_to_yjf_roman ( jed1 );

      s2 = yjf_to_s_roman ( y2, j2, f2 );

      jed3 = yjf_to_jed_roman ( y2, j2, f2 );

      fprintf ( 1, '  %11.2f  %20s  %11.2f\n', jed1, s2, jed3 );

    end

  end

  return
end


text:
classdef chebfun2 < separableApprox
%CHEBFUN2   CHEBFUN2 class for representing functions on [a,b]x[c,d].
%
%   Class for approximating functions defined on finite rectangles. The
%   functions should be smooth.
%
% CHEBFUN2(F) constructs a CHEBFUN2 object representing the function F on
% [-1,1]x[-1 1]. F can be a string, e.g., 'sin(x.*y)', a function handle, e.g.,
% @(x,y) x.*y + cos(x), or a matrix of numbers. For the first two, F should in
% most cases be "vectorized" in the sense that it may be evaluated at a matrix
% of points and returns a matrix output.
%
% CHEBFUN2(F, [A B C D]) specifies a rectangle [A B]x[C D] where the
% function is defined. A, B, C, D must all be finite.
%
% If F is a matrix, F = (f_ij), the numbers f_ij are used as function
% values at tensor Chebyshev points of the 2nd kind. CHEBFUN2(F, 'equi')
% assumes the function values f_ij come from an equispaced tensor grid.
% CHEBFUN2(F, 'equix') assumes the values come from a tensor grid that is
% equispaced in x and 2nd-kind Chebyshev in y. CHEBFUN2(F, 'equiy') assumes
% the values come from a tensor grid that is 2nd-kind Chebyshev in x and
% equispaced in y.
%
% CHEBFUN2(F, k) returns a rank k approximation to F.
%
% CHEBFUN2(F, [m n]) returns a representation of a bivariate polynomial
% with m coefficients in x and n coefficients in y. The polynomial is
% compressed in low rank form and the rank k is still determined adaptively
% (satisfying k<=min(m,n)+1).
%
% CHEBFUN2(F, k, [A B C D]) or CHEBFUN2(F, [m,n], [A B C D]) is nonadaptive in
% rank or degrees, as above, returning a chebfun2 on the domain [A B]x[C D].
%
% CHEBFUN2(F, 'coeffs') where F is a matrix uses the matrix as coefficients
% in a Chebyshev tensor expansion. CHEBFUN2(F, 'coeffsx') uses an expansion
% of coefficients in x and values in y. CHEBFUN2(F, 'coeffsy') uses an
% expansion of values in x and coefficients in y.
%
% CHEBFUN2(F, 'trig') constructs a CHEBFUN2 object representing a smooth
% and periodic function F on [-1,1]x[-1 1]. The resulting CHEBFUN2 is
% represented using a bivariate Fourier expansion. CHEBFUN2(F, 'trigx')
% constructs a CHEBFUN2 that is periodic in x only, represented as a
% bivariate Fourier-Chebyshev expansion. CHEBFUN2(F, 'trigy') constructs a
% CHEBFUN2 that is periodic in y only, represented as a bivariate
% Chebyshev-Fourier expansion.
%
% CHEBFUN2(F, 'periodic')  is the same as CHEBFUN2(F, 'trig').
% CHEBFUN2(F, 'periodicx') is the same as CHEBFUN2(F, 'trigx').
% CHEBFUN2(F, 'periodicy') is the same as CHEBFUN2(F, 'trigy').
%
% CHEBFUN2(F, pref) constructs a CHEBFUN2 object from F with the options
% determined by the CHEBFUNPREF object pref. CHEBFUN2(F, {prefx, prefy})
% uses the CHEBFUNPREF objects prefx in x and prefy in y, where one of
% prefx and prefy may be [].
%
% The Chebfun2 software system is based on:
%
% A. Townsend and L. N. Trefethen, An extension of Chebfun to two dimensions,
% SIAM J. Sci. Comput., 35 (2013), C495-C518.
%
% See also CHEBFUN, CHEBFUN2V.

% Copyright 2017 by The University of Oxford and The Chebfun2 Developers.
% See http://www.chebfun.org/ for Chebfun2 information.

% TODO: Improve documentation of input options.

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS CONSTRUCTOR:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        function f = chebfun2(varargin)
            % The main CHEBFUN2 constructor!

            % Return an empty CHEBFUN:
            if ( (nargin == 0) || isempty(varargin{1}) )
                return
            end

            % Call the constructor, all the work is done here:
            f = constructor(f, varargin{:});

        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% STATIC METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = true )

        % Convert Chebyshev coefficients to values:
        varargout = coeffs2vals(U, varargin);

        % Convert values to Chebyshev coefficients:
        varargout = vals2coeffs(U, varargin);

        % Padua points to tensor grid:
        [C, V, X, Y] = paduaVals2coeffs( F, dom );

        % Tensor product of Chebyshev points:
        [xx, yy] = chebpts2(nx, ny, domain, kind);

        % Outer-product of two chebfuns:
        F = outerProduct(f, g);

        % Fast spectrally-accurate Poisson solver:
        u = poisson(f, varargin);

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Private Static methods implemented by CHEBFUN2 class.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = private, Static = true )

    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function p = kspaceLineRecon(p, dy, dt, c, varargin)
%KSPACELINERECON    2D linear FFT reconstruction.
%
% DESCRIPTION:
%       kspaceLineRecon takes an acoustic pressure time-series p_ty
%       recorded over an evenly spaced array of sensor points on a line,
%       and constructs an estimate of the initial acoustic pressure
%       distribution that gave rise to those measurements using an
%       algorithm based on the FFT. The input p_ty must be indexed
%       p_ty(time step, sensor y position), where the sensor spacing is given
%       by dy, the temporal spacing given by dt, and the sound speed in the
%       propagation medium (which is assumed to be acoustically
%       homogeneous) is given by c. The output p_xy is indexed as
%       p_xy(x position, y position).
%
%       The code uses a k-space algorithm which performs (1) a Fourier
%       transform on the data p_ty along both t and y dimensions (into
%       wavenumber-frequency space, where the wavenumber component is along
%       the detector line), (2) a mapping, based on the dispersion relation
%       for a plane wave in an acoustically homogeneous medium, from
%       wavenumber-frequency space to wavenumber-wavenumber space (where
%       the second component is orthogonal to the detector line), and
%       finally (3) an inverse Fourier transform back from the wavenumber
%       domain to the spatial domain. The result is an estimate of the
%       initial acoustic pressure distribution from which the acoustic
%       waves originated.
%
%       Steps (1) and (3) can be performed efficiently using the fast
%       Fourier transform (FFT); they are therefore fastest when the number
%       of samples and number of detector points are both powers of 2. The
%       mapping in step (2) requires an interpolation of the data from an
%       evenly spaced grid of points in the wavenumber-frequency domain to
%       an evenly-spaced grid of points in the wavenumber-wavenumber
%       domain. The option 'Interp' may be used to choose the interpolation
%       method.
%
%       The physics of photoacoustics requires that the acoustic pressure
%       is initially non-negative everywhere. The estimate of the initial
%       pressure distribution generated by this code may have negative
%       regions due to artefacts arising from differences between the
%       assumed model and the real situation, e.g., homogeneous medium vs.
%       real, somewhat heterogeneous, medium; infinite measurement surface
%       vs. finite-sized region-of-detection, etc. A positivity (or
%       non-negativity) condition can be enforced by setting the optional
%       'PosCond' to true which simply sets any negative parts of the final
%       image to zero.
%
% USAGE:
%       p_xy = kspaceLineRecon(p_ty, dy, dt, c)
%       p_xy = kspaceLineRecon(p_ty, dy, dt, c, ...)
%
% INPUTS:
%       p_ty        - pressure time-series recorded over an evenly spaced
%                     array of sensor points on a line (indexed as t, y)
%       dy          - spatial step [m]
%       dt          - time step [s]
%       c           - acoustically-homogeneous sound speed [m/s]
%
% OPTIONAL INPUTS:
%       Optional 'string', value pairs that may be used to modify the
%       default computational settings.
%
%       'DataOrder' - String input which sets the data order (default =
%                     'ty'). Valid inputs are 'ty' and 'yt'
%       'Interp'    - string input controlling the interpolation method
%                     used by interp2 in the reconstruction (default =
%                     '*nearest')
%       'Plot'      - Boolean controlling whether a plot of the
%                     reconstructed estimate of the initial acoustic
%                     pressure distribution is produced (default = false).
%       'PosCond'   - Boolean controlling whether a positivity condition is
%                     enforced on the reconstructed estimate of the initial
%                     acoustic pressure distribution (default = false).
%
% OUTPUTS:
%       p_xy        - estimate of the initial acoustic pressure
%                     distribution (indexed as x, y)
%
% ABOUT:
%       author      - Bradley Treeby and Ben Cox
%       date        - 11th January 2009
%       last update - 13th August 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also interp2, kspacePlaneRecon, makeGrid

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

% Error: interp2 error in Matlab R2012b
matlab_ver = ver('matlab');
if strcmp(matlab_ver.Release, '(R2012b)')
    error('This function will not run in Matlab R2012b because of a bug in Matlab''s function interp2. It does work in prior and subsequent versions of Matlab.')
end

% start timer
tic;

% define defaults
num_req_inputs = 4;
data_order = 'ty';
interp_method = '*nearest';
plot_recon = false;
positivity_cond = false;

% replace with user defined values if provided
if nargin < num_req_inputs
    error('Incorrect number of inputs');
elseif ~isempty(varargin)
    for input_index = 1:2:length(varargin)
        switch varargin{input_index}
            case 'DataOrder'
                data_order = varargin{input_index + 1};
                if ~strcmp(data_order, 'ty') && ~strcmp(data_order, 'yt')
                    error('Unknown setting for optional input DataOrder');
                end
            case 'Interp'
                interp_method = varargin{input_index + 1};
            case 'Plot'
                plot_recon = varargin{input_index + 1};
            case 'PlotRecon'
                plot_recon = varargin{input_index + 1};
            case 'PosCond'
                positivity_cond = varargin{input_index + 1};
            otherwise
                error('Unknown optional input');
        end
    end
end

% reorder the data if needed (p_ty)
if strcmp(data_order, 'yt')
    p = p.';
end

% mirror the time domain data about t = 0 to allow the cosine transform to
% be computed using an FFT (p_ty)
p = [flipdim(p, 1); p(2:end, :)];

% extract the size of mirrored input data
[Nt, Ny] = size(p);

% update command line status
disp('Running k-Wave line reconstruction...');
disp(['  grid size: ' num2str(Ny) ' by ' num2str((Nt+1)/2) ' grid points']);
disp(['  interpolation mode: ' interp_method]);

% create a computational grid that is evenly spaced in w and ky, where
% Nx = Nt and dx = dt*c
kgrid = makeGrid(Nt, dt*c, Ny, dy);

% from the grid for kx, create a computational grid for w using the
% relation dx = dt*c; this represents the initial sampling of p(w, ky)
w = c*kgrid.kx;

% remap the computational grid for kx onto w using the dispersion
% relation w/c = (kx^2 + ky^2)^1/2. This gives an w grid that is
% evenly spaced in kx. This is used for the interpolation from p(w, ky)
% to p(kx, ky). Only real w is taken to force kx (and thus x) to be
% symmetrical about 0 after the interpolation.
w_new = (c*kgrid.k);

% calculate the scaling factor using the value of kx, where
% kx = sqrt( (w/c).^2 - kgrid.ky.^2 ) and then manually
% replacing the DC value with its limit (otherwise NaN results)
sf = c^2*sqrt( (w/c).^2 - kgrid.ky.^2)./(2*w);
sf(w == 0 & kgrid.ky == 0) = c/2;

% compute the FFT of the input data p(t, y) to yield p(w, ky) and scale
p = sf.*fftshift(fftn(fftshift(p)));

% remove unused variables
clear sf;

% exclude the inhomogeneous part of the wave
p(abs(w) < abs(c*kgrid.ky)) = 0;

% compute the interpolation from p(w, ky) to p(kx, ky)and then force to be
% symmetrical
p = interp2(kgrid.ky, w, p, kgrid.ky, w_new, interp_method);

% remove unused variables
clear kgrid w;

% set values outside the interpolation range to zero
p(isnan(p)) = 0;

% compute the inverse FFT of p(kx, ky) to yield p(x, y)
p = real(ifftshift(ifftn(ifftshift(p))));

% remove the left part of the mirrored data which corresponds to the
% negative part of the mirrored time data
p = p( (Nt + 1)/2:Nt, :);

% correct the scaling - the forward FFT is computed with a spacing of dt
% and the reverse requires a spacing of dy = dt*c, the reconstruction
% assumes that p0 is symmetrical about y, and only half the plane collects
% data (first approximation to correcting the limited view problem)
p = 2*2*p./c;

% enfore positivity condition
if positivity_cond
    disp('  applying positivity condition...');
    p(p < 0) = 0;
end

% update command line status
disp(['  computation completed in ' scaleTime(toc)]);

% plot the reconstruction
if plot_recon

    % allocate axis dimensions
    x_axis = [0 (Nt/2)*dt*c];
    y_axis = [0 Ny*dy];

    % select suitable axis scaling factor
    [x_sc, scale, prefix] = scaleSI(max([Ny*dy, (Nt/2)*dt*c]));

    % select suitable plot scaling factor
    plot_scale = max(p(:));

    % create the figure
    figure;
    imagesc(y_axis*scale, x_axis*scale, p, [-plot_scale, plot_scale]);
    axis image;
    colormap(getColorMap);
    xlabel(['Sensor Position [' prefix 'm]']);
    ylabel(['Depth [' prefix 'm]']);
    colorbar;
end

repo: Aero-matlab


text:
function [xr,Yr]=noisergeo(x,dim,tau,r,q,p,theta)
%Syntax: [xr,Yr]=noisergeo(x,dim,tau,r,q,p,theta)
%________________________________________________
%
% Noise reduction by Local Geometric Projection.
%
% xr is the vector/matrix with the cleaned time series.
% Yr is the phase space of the last cleaned xr.
% x is the time series.
% dim is the embedding dimension.
% tau is the time delay.
% r can be either
%   real defining the neighborhood range.
%   integer defining the number of nearest neighbors.
% q can take one of the following values
%  'wAV' for the weighted average.
%  [an integer from 0 to dim-1] for the local geometric projection.
%  'mod' for the adaptive selection of the local neighborhood dimensions.
% p defines the norm.
% theta is the correction paprameter [0,1]
%   1: full correction.
%   0: no correction.
%
%
% References:
%
% Kantz H, Schreiber T, Hoffmann I, Buzug T, Pfister G, Flepp L G, Simonet
% J, Badii R, Brun E (1993): Nonlinear noise reduction: A case study on
% experimental data. Physical Review E 48: 1529-1538
%
% Leontitsis A., Bountis T., Pange J. (2004): An adaptive way for improving
% noise reduction using local geometric projection. CHAOS 14(6): 106-110
%
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 14 Jul 2001

if nargin<1 | isempty(x)==1
   error('You should provide a time series.');
else
   % x must be a vector
   if min(size(x))>1
      error('Invalid time series.');
   end
   x=x(:);
   % n is the time series length
   n=length(x);
end

if nargin<2 | isempty(dim)==1
   dim=2;
else
   % dim must be either a scalar or a vector
   if min(size(dim))>1
      error('dim must be a scalar or a vector.');
   end
   % dim must be an integer
   if round(dim)-dim~=0
      error('dim must be an integer.');
   end
   % dim values must be above 1
   if any(dim<1)==1
      error('dim values must be above 1');
   end
end

if nargin<3 | isempty(tau)==1
   tau=1;
else
   % tau must be either a scalar or a vector
   if min(size(tau))>1
      error('tau must be a scalar or a vector.');
   end
   % tau must be an integer
   if round(tau)-tau~=0
      error('tau must be an integer.');
   end
   % tau values must be above 1
   if any(tau<1)==1
      error('tau values must be above 1');
   end
end

if nargin<4 | isempty(r)==1
   r=dim+1;
else
   % r must be either a scalar or a vector
   if min(size(r))>1
      error('r must be a scalar or a vector.');
   end
   % r values must be above 0
   if any(r<=0)==1
      error('r values must be greater than 0');
   end
end

if nargin<5 | isempty(q)==1
   q=1;
end

if nargin<6 | isempty(p)==1
   p=2;
else
   % p must be either a scalar or a vector
   if min(size(p))>1
      error('p must be a scalar or a vector.');
   end
end

if nargin<7 | isempty(theta)==1
   theta=1;
else
   % theta must be either a scalar or a vector
   if min(size(theta))>1
      error('theta must be a scalar or a vector.');
   end
   % theta must be above 0 and bellow 1
   if theta<0 | theta>1
      error('theta must be above 0 and bellow 1.')
   end
end

% Only one of dim, tau, r, p or theta should be vector
l=[length(dim),length(tau),length(r),length(p),length(theta)];
if length(find(l>1))>1
   error('Only one of dim, tau, r, p, or theta should be vector.');
end

% Make the phase-space
[Y,T]=phasespace(x,dim,tau);

m=max(l);
dim=ones(1,m).*dim;
tau=ones(1,m).*tau;
r=ones(1,m).*r;
p=ones(1,m).*p;
theta=ones(1,m).*theta;

for i=1:m

    % Initialize Yr
    Yr=zeros(T,dim(i));

    % For every phase-space point
    for j=1:T

        % Locate the j-th point
        y=Y(j,:);

        % Check neighborhood or neighbors
        if mod(r(i),floor(r(i)))==0
            lock=Knearest(y,Y,r(end)+1,p(i));
        else
            lock=radnearest(y,Y,T,r(i),p(i));
            lock(find(j==lock))=[];
        end

        if isempty(lock)==1
            Yr(j,:)=y;
        elseif q=='wAV'
            % The calculations for the weighted average
            Ynearest=Y(lock,:);
            w=[];
            for k=1:length(lock)
                w(k)=norm(y-Ynearest(k,:))/2/r(i);
                w(k)=exp(-w(k)^2);
                Ynearest(k,:)=Ynearest(k,:)*w(k);
            end
            Yr(j,:)=sum(Ynearest)/sum(w);
        else
            % All the neighboring points have equal weight
            if length(lock)<dim(i)
                Yr(j,:)=mean(Y(lock,:));
            else
                % Make the matrix of the local neighborhood
                Ynearest=Y(lock,:);
                %Ynearest=quality3(Ynearest,2);
                %Ynearest=quality2(Ynearest,2,.99,r,y);
                Ynearest=quality(Ynearest,2,.7);
                % Neigborhood average
                ybar=mean(Ynearest);
                % Subtract the average from the neighborhood poionts
                Ynearest=Ynearest-ones(size(Ynearest,1),1)*ybar;
                % SVD on the local neighborhood
                [u,s,v]=svd(Ynearest,0);
                % Either use the Adaptive selection ...
                if q=='mod'
                    s=diag(s);
                    d=s(1:dim(i)-1)./(s(2:dim(i))+eps);
                    q=find(d==max(d));
                    Yr(j,:)=y-(y-ybar)*(v(:,q+1:dim(i))*v(:,q+1:dim(i))');
                % ... or use a predifined q
                else
                    Yr(j,:)=y-(y-ybar)*(v(:,dim(i)-q:dim(i))*v(:,dim(i)-q:dim(i))');
                end
            end
        end
        Yr(j,:)=theta(i)*Yr(j,:)+(1-theta(i))*Y(j,:);
    end

    % Initialize xr
    xr(:,i)=zeros(n,1);
    times(:,i)=zeros(n,1);

    % Calculate the reconstructed time series
    for j=1:dim
        xr(1+(j-1)*tau:T+(j-1)*tau,i)=xr(1+(j-1)*tau:T+(j-1)*tau,i)+Yr(:,j);
        times(1+(j-1)*tau:T+(j-1)*tau,i)=times(1+(j-1)*tau:T+(j-1)*tau,i)+1;
    end
    xr(:,i)=xr(:,i)./times(:,i);
    xr(:,i)=xr(:,i)-mean(xr(:,i)-x);

end


text:
function edges = crackPattern2(box, points, alpha, varargin)
%CRACKPATTERN2 Create a (bounded) crack pattern tessellation.
%
%   E = crackPattern2(BOX, POINTS, ALPHA)
%   create a crack propagation pattern wit following parameters :
%   - pattern is bounded by area BOX which is a polygon.
%   - each crack originates from points given in POINTS
%   - directions of each crack is given by a [NxM] array ALPHA, where M is
%   the number of rays emanating from each seed/
%   - a crack stop when it reaches another already created crack.
%   - all cracks stop when they reach the border of the frame, given by box
%   (a serie of 4 points).
%   The result is a collection of edges, in the form [x1 y1 x2 y2].
%
%   E = crackPattern2(BOX, POINTS, ALPHA, SPEED)
%   Also specify speed of propagation of each crack.
%
%
%   See the result with :
%     figure;
%     drawEdge(E);
%
%   See also drawEdge

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2004-05-25
% Copyright 2004-2022 INRA - TPV URPOI - BIA IMASTE

if ~isempty(varargin)
    speed = varargin{1};
else
    speed = ones(size(points, 1), 1);
end

% Compute line equations for each initial crack.
% The 'Inf' at the end correspond to the position of the limit.
% If an intersection point is found with another line, but whose position
% is after this value, this means that another crack stopped it before it
% reach the intersection point.
NP = size(points, 1);
lines = zeros(0, 5);
for i = 1:size(alpha, 2)
    lines = [lines; points speed.*cos(alpha(:,i)) speed.*sin(alpha(:,i)) Inf*ones(NP, 1)]; %#ok<AGROW>
end
NL = size(lines, 1);

% initialize lines for borders, but assign a very high speed, to be sure
% borders will stop all cracks.
dx = (box([2 3 4 1],1)-box([1 2 3 4],1))*max(speed)*5;
dy = (box([2 3 4 1],2)-box([1 2 3 4],2))*max(speed)*5;

% add borders to the lines set
lines = [lines ; createLine(box, dx, dy) Inf*ones(4,1)];

edges = zeros(0, 4);

while true
    modif = 0;

    % try to update each line
	for i=1:NL

        % initialize first point of edge
        edges(i, 1:2) = lines(i, 1:2);

        % compute intersections with all other lines
        pi = intersectLines(lines(i,:), lines);

        % compute position of all intersection points on the current line
        pos = linePosition(pi, lines(i,:));

        % consider points to the right (positive position), and sort them
        indr = find(pos>1e-12 & pos~=Inf);
        [posr, indr2] = sort(pos(indr));

        % look for the closest intersection to the right
        for i2=1:length(indr2)

            % index of intersected line
            il = indr(indr2(i2));

            % position of point relative to intersected line
            pos2 = linePosition(pi(il, :), lines(il, :));

            % depending on the sign of position, tests if the line2 can
            % stop the current line, or if it was stopped before
            if pos2>0
                if pos2<abs(posr(i2)) && pos2<lines(il, 5)
                    if lines(i, 5) ~= posr(i2)
                        edges(i, 3:4) = pi(il,:);
                        lines(i, 5) = posr(i2);
                        modif = 1;
                    end
                    break;
                end
            end
        end   % end processing of right points of the line

	end % end processing of all lines

    % break the infinite loop if no more modification was made
    if ~modif
        break;
    end
end

% add edges of the surronding box.
edges = [edges ; box(1,:) box(2,:) ; box(2,:) box(3,:); ...
                 box(3,:) box(4,:) ; box(4,:) box(1,:)  ];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
%SCALEM Compute scaling map
%
%   W = SCALEM(A,TYPE)
%   W = A*SCALEM([],TYPE)
%   W = A*SCALEM(TYPE)
%
% INPUT
%   A     Dataset
%   TYPE  Type of scaling (optional; default: the class priors weighted
%         mean of A is shifted to the origin)
%
% OUTPUT
%   W  Scaling mapping
%
% DESCRIPTION
% Computes a scaling map W, whose type depends on the parameter T:
%
% [], 'c-mean'  - The mean of A is shifted to the origin. Class priors are taken
%                 into account.
% 'mean'        - The mean of A is shifted to the origin. This is computed for
%                 the data as given in A, neglecting class priors.
% 'c-variance'  - The mean of A  is shifted to the origin and the average class
%                 variances (within-scatter) are normalized.  Class priors are
%                 taken into account.
% 'variance'    - The mean of A is shifted to the origin and the total variances
%                 of all features are scaled to 1. This is computed for the
%                 data as given in A, neglecting class priors.
% '2-sigma'     - For each feature f, the 2-sigma interval, [f-2*std(f),f+2*std(f)]
%                 is rescaled to [0,1]. Values outside this domain are clipped.
%                 Class priors are taken into account.
% 'domain'      - The domain for all features in the dataset A is set to [0,1].
%                 This is computed for the data as given in A, neglecting class priors.
%
% The map W may be applied to a new dataset B by B*W.
%
% SEE ALSO (PRTools Guide)
% PRMAPPING, DATASET

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: scalem.m,v 1.9 2009/11/27 08:52:02 duin Exp $

function W = scalem(varargin)

  argin = shiftargin(varargin,'char');
  argin = setdefaults(argin,[],'');
  if mapping_task(argin,'definition')
  	W = define_mapping(argin,'untrained');
    return
  end

  [a,t] = deal(argin{:});
  %DXD Make the naming even better:
  dname = 'Scaling mapping';
  switch t
  case 'mean'
	  dname = 'zero-mean';
  case {'','c-mean'}
	  dname = 'zero-mean (+class priors)';
  case 'variance'
	  dname = 'unit-var';
  case 'c-variance'
	  dname = 'unit-var (+class priors)';
  case '2-sigma'
	  dname = '0-1 scaling (+-2 standard dev.)';
  case 'domain'
	  dname = '0-1 scaling';
  end
  % No data, return an untrained mapping.
	if (nargin < 1) | isempty(a)
		W = prmapping('scalem',{t});
		W = setname(W,dname);
		return;
	end

	[a,c,lablist] = cdats(a,1);
	[m,k] = size(a);

	switch t
	case {'','c-mean'}
    u = meancov(a);
		u =getprior(a)*(+u);
		s = ones(1,k);
		clip = 0;

	case 'mean'
		a = remclass(a);
		u = mean(a,1);
		s = ones(1,k);
		clip = 0;

	case 'c-variance'
		p = getprior(a);
		U = meancov(a);
		u = p*(+U);
		s = zeros(1,k);
		for j=1:c
			s = s + p(j)*var(seldat(a,j)*cmapm(+U(j,:),'shift'));
		end
		s = sqrt(s);
		clip = 0;

	case 'variance'
		%[u,G] = meancov(+a);
		%s = sqrt(diag(G))';
  	u = mean(a,1);
		s = std(a,0,1);
		clip = 0;

	case '2-sigma'
		p = getprior(a);
		U = meancov(a);
		u = p*(+U);
		s = zeros(1,k);
		for j=1:c
			s = s + p(j)*var(seldat(a,j));
		end
		s = sqrt(s);
		s = 4*s;
		u = u-0.5*s;
		clip = 1;

	case 'domain'
		mx = max(a,[],1)+eps;
		mn = min(a,[],1)-eps;
		u = mn;
		s = (mx - mn)*(1+eps);
		clip = 0;

	otherwise
		error('Unknown option')
	end

	J = find(s==0);
	s(J) = ones(size(J)); 				% No scaling for zero variances.

	W = cmapm({u,s,clip},'scale');
	W.size_out = a.featsize;
	W = setlabels(W,getfeatlab(a));

return;

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
%% (Internal) Get the time samples of a (xx,yy) sequence below two restrictions (x_thr, y_thr)
%
% This function calculates the time samples (xx) of a two valued sequence
% (xx, yy) below two restrictions (x_thr, y_thr)
%
%       selected_segments = get_segments_from_sequence(yy, xx, y_thr, x_thr )
%
% Arguments:
%
%       xx, yy: xx and yy values.
%
%       x_thr, y_thr: the thresholds to search for within the sequence.
%
% Output:
%
%       selected_segments: a (n x 2) matrix with the segments start-end
%                          where n is the amount of segments found in the
%                          sequence
%
% Example
%
% See also RR_calculation, MedianFiltSequence
%
% Author: Mariano Llamedo Soria (llamedom at frba.utn.edu.ar)
% Version: 0.1 beta
% Birthdate  : 9/5/2017
% Last update: 9/5/2017
% Copyright 2008-2017
%
function selected_segments = get_segments_from_sequence(xx, yy, x_thr, y_thr )

    selected_segments = [];

    % sequence values below the thr
    bAux = yy <= y_thr;
    % indexes above the thr
    aux_idx = find(~bAux);

    aux_val = diff(xx(aux_idx));
    % starts of segments at least x_thr long that are below y_thr
    aux_idx3 = find(aux_val >= x_thr);
    below_idx = xx(bAux);
    above_idx = xx(~bAux);

% figure(3); plot(xx, yy ); xlims = xlim(); ylims = ylim(); ylims = ylims + [0.1 -0.1] * diff(ylims); hold on; plot(xlims, [y_thr y_thr] , '--r' ); hold off;

    for jj = 1:length(aux_idx3)
        start_idx = xx(aux_idx(aux_idx3(jj))+1);
        end_idx = find( above_idx > start_idx, 1,'first');
        if isempty(end_idx)
            continue
        end
        end_idx = below_idx( find( below_idx < above_idx(end_idx) ,1, 'last') );

        if( ~isempty(start_idx) && ~isempty(end_idx) )
            selected_segments = [selected_segments; [start_idx end_idx]];
        end
    end

end


text:
function test_bug2563
% DEPENDENCY ft_selectdata getdimord

% WALLTIME 00:10:00
% MEM 1gb

nsubj = 3;
nchan = 4;
nfreq = 5;

freq = [];
freq.dimord = 'subj_chan_freq';
freq.MIspctrm = randn(nsubj,nchan,nfreq);
freq.freq = 1:nfreq;
for i=1:nchan
  freq.label{i} = num2str(i);
end

cfg = [];
cfg.avgoverfreq = 'yes';
output = ft_selectdata(cfg, freq);

assert(isequal(size(output.MIspctrm), [3 4]));

disp(freq)
disp(output)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the following was reported as http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2561
% but seems to be related, i.e. both have subjects in the data

timelock = [];
timelock.individual = randn(7,204,1500);
timelock.time = 1:1500;
for i=1:204
  timelock.label{i} = num2str(i);
end
timelock.dimord = 'subj_chan_time';

cfg = [];
output = ft_selectdata(cfg, timelock);

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%DECOMPOSEHOMOGRAPHYMAT  Decompose a homography matrix to rotation(s), translation(s) and plane normal(s)
%
%     [motions, nsols] = cv.decomposeHomographyMat(H, K)
%
% ## Input
% * __H__ The input homography matrix between two images, 3x3.
% * __K__ The input intrinsic camera calibration matrix, 3x3.
%
% ## Output
% * __motions__ Decomposed `H`. A scalar struct with the following fields:
%   * __R__ Array of rotation matrices. Cell array of 3x3 rotations.
%   * __t__ Array of translation matrices. Cell array of 3x1 translations.
%   * __n__ Array of plane normal matrices. Cell array of 3x1 normals.
% * __nsols__ number of solutions.
%
% This function extracts relative camera motion between two views observing a
% planar object from the homography `H` induced by the plane. The intrinsic
% camera matrix `K` must also be provided. The function may return up to four
% mathematical solution sets. At least two of the solutions may further be
% invalidated if point correspondences are available by applying positive
% depth constraint (all points must be in front of the camera). The
% decomposition method is described in detail in [Malis]:
%
%     H = s * K * (R{i} + t{i} * n{i}') * inv(K),  for i=1:nsols
%
% `s` is a scaling factor approximated as:
%
%     s = median(svd(inv(K) * H * K))
%
% ## References
% [Malis]:
% > Ezio Malis, Manuel Vargas, and others. "Deeper understanding of the
% > homography decomposition for vision-based control". 2007.
%
% See also: cv.findHomography, cv.decomposeEssentialMat
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 	TEST THE DYNAMICS OF A 2DOF robot under different circumstances.
%
%   Do exercises 1, 2 and 3
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2016, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function exercise_inv_dynamics_2dofplanar()
close all

%Just load the robot once
robot=load_robot('example', '2dofplanar');

%TODO: uncommment as you solve the exercises
exercise1(robot)
%exercise2(robot)
%exercise3(robot)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Given the following specifications, find the maximum torques if:
%
%   g = [0 -9.81 0]'
%   amax = [2 2] rad/s^2, max angular acceleration
%   wmax = [3 3], max angular speed
%
%   OBTAIN tau for different cases and find the maximum value.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function exercise1(robot)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%TODO: CHANGE q to find the "worst" pose
<<<<<<< HEAD
q=[pi/2 0]';
=======
q=[0 pi/2]';
>>>>>>> develop_new_simulations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
qd=[0 0]'; %speed --> does it make any difference?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TODO: CHANGE qdd [-3 3] to find whether higher torques are achieved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
<<<<<<< HEAD
qdd = [0 0]'; %--> max acceleration
=======
qdd = [1 1]'; %--> max acceleration
>>>>>>> develop_new_simulations

g=[0 -9.81 0]'; % acting on the Y0 axis.
fext = [0 0 0 0 0 0]'; %[fx fy fz nx ny nz]
figure, drawrobot3d(robot, q)

%Compute the torques
%TODO: store the results and find the maximum torque
tau = inversedynamic(robot, q, qd, qdd, g, fext)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Given the following specifications, find the maximum torques if:
%
%   g = [0 -9.81 0]'
%   amax = [2 2] rad/s^2, max angular acceleration
%   wmax = [3 3], max angular speed
%
%   The robot carries a payload with m=2kg at the end effector.
%
%   OBTAIN tau for different cases and find the maximum value.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function exercise2(robot)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%TODO: CHANGE q to find the "worst" pose
q=[0 0]';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
qd=[0 0]'; %speed --> does it make any difference?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TODO: CHANGE qdd [-3 3] to find whether higher torques are achieved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
qdd = [0 0]'; %--> max acceleration

m = 2;          %kg
g=[0 -9.81 0]'; % acting on the Y0 axis.

%forces and moments are expressed in the base reference frame
% TODO express the external forces and moments
fext = [  ]';

figure, drawrobot3d(robot, q)

%Compute the torques
tau = inversedynamic(robot, q, qd, qdd, g, fext)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Given the following specifications, find the maximum torques if the robot
%   follows a trajectory from points in cartesian coordinates:
%
%   initial = (1.5, 0.2)
%   final = (0.2,1.3)
%
%   assume that acceleration stays constant at:
%
%       amax = [2 2] rad/s^2, max angular acceleration
%       wmax = [3 3], max angular speed
%   The robot carries a piece at the end effector with mass=2kg.
%   g = [0 -9.81 0]'
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function exercise3(robot)
robot.q = [0 0]
p0 = [1.5, 0.2]; %initial point
p1 = [0.2  1.3]
speed = 0.5; %m/s, linear speed

distance = norm(p1-p0);
%unit vector in the direction of the trajectory normalized
u = (p1-p0)/norm(p1-p0);
dist= linspace(0,distance,100);

%initialize T
T=eye(3);
T(1:2,4) = p0';
q_line = []
%Generate a number of
for j = 1:length(dist),
    %find all the possible solutions to reach T
    qinv = inversekinematic(robot, T);

    %Select the closest coordinates from the two possible solutions of the
    %inverse kinematics.
    q=select_closest_joint_coordinates(qinv, robot.q);

    %store the joint coordinates in a vector for ulterior animation
    q_line = [q_line q];
    robot.q=q;

    pj = p0 + u*dist(j);
    T(1:2,4) = pj';
end

%
% At this point, the joint coordinates needed are stored at q_line
%

figure, plot(q_line(1,:)), hold on
plot(q_line(2,:)), xlabel('time step'), title('Joint coordinates for a line in space')
legend('q_1', 'q_2')

q=[0 0]';
qd=[3 3]'; %speed --> does it make any difference?
qdd=[2 2]'; %--> max acceleration
m = 2;          %kg
g=[0 -9.81 0]'; % acting on the Y0 axis.

%forces and moments are expressed in the last reference frame
% TODO: express the external force.
fext = [  ]';

figure, drawrobot3d(robot, q)

%
%   TODO:   Compute tau for each q_line().
%           Store tau.
%

figure, plot(taus(1,:)), hold on
plot(taus(2,:)), xlabel('time step'), title('Torques when the robot follows a line in space')
legend('\tau_1', '\tau_2')

repo: ARTE


text:
function [Cov]=vbfa(y,nl,nem,fig)

% Output
% Regularized noise covariance from pre-stimulus data

% Input
% y(nk,nt) = Baseline or Prestimulus data
% b(nk,nl) = mixing matrix
% lam(nk,1) = sensor noise precision
% bet(nl,1) = hyperparamaeters
% sig(nk,nk) = data covariance matrix = b*b'+1/lam (approx)
% nk = data dimensionality
% nl  = factor dimensionality (try 5)
% nt = number of time points
% nem = number of EM iterations  (try 50)
% plot = flag for plotting

% experiment with nl,nem using the plots generated by the code:
% set nem such that the likelihood converges (top plot)
% set nl such that some hyperparameters bet approach infinity, i.e. some
% 1/bet vanish (middle plot)

%===============================

nk=size(y,1);nt=size(y,2);

ryy=y*y';
[p d]=svd(ryy/nt);
d=diag(d);
b=p*diag(sqrt(d));
b=b(:,1:nl);
lam=1./diag(ryy/nt);

bet=ones(nl,1); %diff from VBFA

likeli=zeros(nem,1);
rbb=eye(nl)/nt;

for iem=1:nem
   dbet=diag(bet);ldbet=sum(log(bet));
   dlam=diag(lam);ldlam=sum(log(lam));

   gam=b'*dlam*b+eye(nl)+nk*rbb*0; %diff from VBFA
   igam=inv(gam);
   ubar=igam*b'*dlam*y;
   ryu=y*ubar';
   ruu=ubar*ubar'+nt*igam;

   [p d q]=svd(gam);ldgam=sum(log(diag(d)));
   temp1=-.5*ldgam*ones(1,nt)+.5*sum(ubar.*(gam*ubar),1);
   temp2=.5*ldlam*ones(1,nt)-.5*lam'*(y.^2);
   f=temp1+temp2;
   f3=.5*nl*ldlam+.5*nk*ldbet-.5*trace(b'*dlam*b*dbet); %diff
   likeli(iem)=mean(f)+f3/nt;

   betbar=ruu+dbet;
   ibetbar=inv(betbar);
   b=ryu*ibetbar;

   ilam=diag(ryy-b*ryu')/(nt+nl); % diff
   lam=1./ilam;
   dlam=diag(lam);

   bet=1./(diag(b'*dlam*b)/nk+diag(ibetbar));
    if nargin>3
       figure(fig)
       hsub=subplot(3,3,1);plot((1:iem)',likeli(1:iem));title('likelihood')
       hsub=subplot(3,3,4);plot((1:nl)',sqrt([mean(b.^2,1)' 1./bet]));title('1/bet');
       hsub=subplot(3,3,7);plot(1./lam);title('1/lam');
       drawnow;
    end

   rbb=ibetbar;
end

weight=b*igam*b'*dlam;
sig=b*b'+diag(1./lam);
yc=b*ubar;
cy=b*ruu*b'+diag(ilam*trace(ruu*ibetbar));
mlike=likeli(iem);

if nargin>3
    figure(fig)
    subplot(3,3,2);imagesc(ryy/nt);title('ryy/nt');colorbar;
    subplot(3,3,3);imagesc(cy/nt);title('cy/nt');colorbar;
    subplot(3,3,5);imagesc((ryy-cy)/nt);title('(ryy-cy)/nt');colorbar;
    subplot(3,3,6);imagesc(b*b');title('b*bT');colorbar;
    subplot(3,3,8);imagesc(sig);title('sig');colorbar;
end

Cov=cy;

return

%===============================


text:
function vol = imVolumeEstimate(img, varargin)
% Estimate volume of a 3D binary structure with edge correction.
%
%   Vest = imVolumeEstimate(IMG);
%   Vest = imVolumeEstimate(IMG, DELTA);
%
%   Example
%   imVolumeEstimate
%
%   See also
%     imVolume, imSurfaceEstimate

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2010-01-21,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

%% Input arguments processing

% check image dimension
if ndims(img) ~= 3
    error('first argument should be a 3D binary or label image');
end

% in case of a label image, return a vector with a set of results
if ~islogical(img)
    labels = unique(img);
    labels(labels==0) = [];
    vol = zeros(length(labels), 1);
    for i = 1:length(labels)
        vol(i) = imVolumeEstimate(img==labels(i), varargin{:});
    end
    return;
end

% check image resolution
delta = [1 1 1];
if ~isempty(varargin)
    delta = varargin{1};
end

%% main processing

% compute volume in whole image
vol = sum(img(:));

% compute volume on border faces
f1 = sum(sum(sum(img([1 end], :, :))));
f2 = sum(sum(sum(img(:, [1 end], :))));
f3 = sum(sum(sum(img(:, :, [1 end]))));

% compute volume on border edges
e1 = sum(sum(sum(img(:, [1 end], [1 end]))));
e2 = sum(sum(sum(img([1 end], :, [1 end]))));
e3 = sum(sum(sum(img([1 end], [1 end], :))));

% compute volume on corners
v = sum(sum(sum(img([1 end], [1 end], [1 end]))));

% estimate area using edge weighting according to multiplicity
vol = vol -(f1+f2+f3)/2 + (e1+e2+e3)/4 - v/8;

% multiply by volume of a single voxel
vol = vol * prod(delta);

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


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

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

    % Last Modified by GUIDE v2.5 28-Jun-2015 14:21:02

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

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

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

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

    celBodyData = varargin{1};
    setappdata(hObject,'celBodyData',celBodyData);

    pos = get(handles.dispAxes,'Position');
    set(handles.displayAtmoRadio,'UserData',pos);

    setupCelBodyListbox(handles, celBodyData);
    catalogListbox_Callback(handles.catalogListbox, [], handles);

    % Update handles structure
    guidata(hObject, handles);

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

function setupCelBodyListbox(handles, celBodyData)
    bodiesStr = ma_getSortedBodyNames(celBodyData);

    set(handles.catalogListbox, 'String', bodiesStr);
    if(length(bodiesStr) >= 1)
        set(handles.catalogListbox, 'value', 1);
    end

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

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

    % --- Executes on selection change in catalogListbox.
function catalogListbox_Callback(hObject, eventdata, handles)
    % hObject    handle to catalogListbox (see GCBO)
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)

    % Hints: contents = cellstr(get(hObject,'String')) returns catalogListbox contents as cell array
    %        contents{get(hObject,'Value')} returns selected item from catalogListbox
    contents = cellstr(get(hObject,'String'));
    selected = contents{get(hObject,'Value')};

    celBodyData = getappdata(handles.ma_CelBodyCatalogGUI,'celBodyData');
    bodyInfo = celBodyData.(strtrim(lower(selected)));
    parentBodyInfo = bodyInfo.getParBodyInfo(celBodyData);
    bodyInfoStr = generateBodyInfoStr(bodyInfo, celBodyData);
    set(handles.infoLabel,'String',bodyInfoStr);

    dispBodyValue = get(handles.displayBodyRadio, 'Value');
    dispOrbValue = get(handles.displayOrbitRadio, 'Value');
    dispAtmoValue = get(handles.displayAtmoRadio, 'Value');

    if(bodyInfo.atmohgt == 0)
        set(handles.displayAtmoRadio, 'Enable', 'off');
        if(dispAtmoValue==1)
            set(handles.displayAtmoRadio, 'Value', 0);
            set(handles.displayBodyRadio, 'Value', 1);

            dispBodyValue = get(handles.displayBodyRadio, 'Value');
            dispAtmoValue = get(handles.displayAtmoRadio, 'Value');
        end
    else
        set(handles.displayAtmoRadio, 'Enable', 'on');
    end

    cla(handles.dispAxes);
    xlabel('');
    ylabel('');
    if(dispBodyValue==1)
        plotBody(handles.ma_CelBodyCatalogGUI,handles.dispAxes,bodyInfo);
        set(handles.dispAxes,'Position',get(handles.displayAtmoRadio,'UserData'));
    elseif(dispOrbValue==1)
        if(~isempty(parentBodyInfo))
            plotBodyOrbit(bodyInfo, 'r', getParentGM(bodyInfo, celBodyData));
            plotBody(handles.ma_CelBodyCatalogGUI,handles.dispAxes,parentBodyInfo);
        else
            plotBody(handles.ma_CelBodyCatalogGUI,handles.dispAxes,bodyInfo);
        end
        set(handles.dispAxes,'Position',get(handles.displayAtmoRadio,'UserData'));
    elseif(dispAtmoValue==1)
        plotAtmoCurve(handles.dispAxes, bodyInfo, celBodyData);
        set(handles.dispAxes,'OuterPosition',get(handles.displayAtmoRadio,'UserData'));
    end

function plotAtmoCurve(orbitDispAxes, bodyInfo, celBodyData)
    atmoHgt = bodyInfo.atmohgt;

    if(atmoHgt > 0)
        hArr = linspace(0,atmoHgt,100);
        density = zeros(size(hArr));
        for(i=1:length(hArr))
            density(i) = getAtmoDensityAtAltitude(bodyInfo, hArr(i), deg2rad(45), 0, 0);
        end

        plot(orbitDispAxes,hArr,density,'r-','LineWidth',2);
        v = ylim(orbitDispAxes);
        grid on;
        xlim(orbitDispAxes, [min(hArr), max(hArr)]);
        ylim(orbitDispAxes, [0, v(2)]);
        xlabel('Altitude [km]');
        ylabel('Density [kg/m^3]');
    end

function plotBody(hFig,orbitDispAxes, bodyInfo)
    ma_initOrbPlot(hFig, orbitDispAxes, bodyInfo);

function bodyInfoStr = generateBodyInfoStr(bodyInfo, celBodyData)
    hrule = '-----------------------------------';

    [year, day, hour, minute, sec] = convertSec2YearDayHrMnSec(bodyInfo.epoch);
    [dateStr] = formDateStr(year, day, hour, minute, sec);

    parentBodyInfo = bodyInfo.getParBodyInfo(celBodyData);
    if(bodyInfo.sma > 0 && ~isempty(parentBodyInfo))
        [~, ~, secInDay, ~] = getSecondsInVariousTimeUnits();
        period = computePeriod(bodyInfo.sma, getParentGM(bodyInfo, celBodyData));
        periodStr = [num2str(period/secInDay), ' days'];

        [rAp, rPe] = computeApogeePerigee(bodyInfo.sma, bodyInfo.ecc);
        rApStr = [num2str(rAp), ' km'];
        rPeStr = [num2str(rPe), ' km'];
    else
        periodStr = 'N/A';

        [~, rPe] = computeApogeePerigee(bodyInfo.sma, bodyInfo.ecc);
        rApStr = 'N/A';
        rPeStr = [num2str(rPe), ' km'];
    end

    if(bodyInfo.sma == 0)
        rPeStr = 'N/A';
    end

    if(isempty(bodyInfo.parent))
        parentStr = 'N/A';
        rSOIStr = 'N/A';
        rSOI = Inf;
    else
        parentStr = bodyInfo.parent;

        rSOI = getSOIRadius(bodyInfo, parentBodyInfo);
        rSOIStr = [num2str(rSOI), ' km'];
    end

    sPhr = floor(bodyInfo.rotperiod/3600);
    sPmin = floor((bodyInfo.rotperiod - 3600*sPhr)/60);
    sPsec = bodyInfo.rotperiod - 3600*sPhr - 60*sPmin;
    sPStr = sprintf('%02.0f:%02.0f:%07.4f',sPhr,sPmin,sPsec);

    syncSMA = computeSMAFromPeriod(bodyInfo.rotperiod, bodyInfo.gm);
    if(syncSMA > rSOI)
        syncNotPoss = '*';
    else
        syncNotPoss = '';
    end

    bodyInfoStr{1} = ['Name: ', bodyInfo.name];
    bodyInfoStr{end+1} = ['ID: ', num2str(bodyInfo.id)];
    bodyInfoStr{end+1} = ['Parent: ', parentStr];
    bodyInfoStr{end+1} = hrule;

    bodyInfoStr{end+1} = 'Orbit';
    bodyInfoStr{end+1} = [' Epoch = ', dateStr];
    bodyInfoStr{end+1} = [' SMA = ', num2str(bodyInfo.sma), ' km'];
    bodyInfoStr{end+1} = [' Ecc = ', num2str(bodyInfo.ecc)];
    bodyInfoStr{end+1} = [' Inc = ', num2str(bodyInfo.inc), ' deg'];
    bodyInfoStr{end+1} = [' RAAN = ', num2str(bodyInfo.raan), ' deg'];
    bodyInfoStr{end+1} = [' Arg = ', num2str(bodyInfo.arg), ' deg'];
    bodyInfoStr{end+1} = [' Mean = ', num2str(bodyInfo.mean), ' deg'];
    bodyInfoStr{end+1} = hrule;

    bodyInfoStr{end+1} = 'Physical Properties';
    bodyInfoStr{end+1} = [' GM = ', num2str(bodyInfo.gm), ' km^3/s^2'];
    bodyInfoStr{end+1} = [' Radius = ', num2str(bodyInfo.radius), ' km'];
    bodyInfoStr{end+1} = [' Atmo Height = ', num2str(bodyInfo.atmohgt), ' km'];
    bodyInfoStr{end+1} = [' Atmo SL Press = ', num2str(getPressureAtAltitude(bodyInfo, 0)), ' kPa'];
    bodyInfoStr{end+1} = [' Sidereal Day = ', num2str(bodyInfo.rotperiod), ' sec'];
    bodyInfoStr{end+1} = ['              = ', sPStr];
    bodyInfoStr{end+1} = hrule;

    bodyInfoStr{end+1} = 'Derived';
    bodyInfoStr{end+1} = [' Orbit Period = ', periodStr];
    bodyInfoStr{end+1} = [' Orbit RPe = ', rPeStr];
    bodyInfoStr{end+1} = [' Orbit RAp = ', rApStr];
    bodyInfoStr{end+1} = [' SoI Radius = ', rSOIStr];
    bodyInfoStr{end+1} = [' Sync. SMA = ', num2str(syncSMA), ' km',syncNotPoss];

    % --- Executes during object creation, after setting all properties.
function catalogListbox_CreateFcn(hObject, eventdata, handles)
    % hObject    handle to catalogListbox (see GCBO)
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    empty - handles not created until after all CreateFcns called

    % Hint: listbox controls usually have a white background on Windows.
    %       See ISPC and COMPUTER.
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end

    % --- Executes when selected object is changed in dispAxesButtonGrp.
function dispAxesButtonGrp_SelectionChangeFcn(hObject, eventdata, handles)
    % hObject    handle to the selected object in dispAxesButtonGrp
    % eventdata  structure with the following fields (see UIBUTTONGROUP)
    %	EventName: string 'SelectionChanged' (read only)
    %	OldValue: handle of the previously selected object or empty if none was selected
    %	NewValue: handle of the currently selected object
    % handles    structure with handles and user data (see GUIDATA)
    catalogListbox_Callback(handles.catalogListbox, [], handles);

    % --- Executes on key press with focus on ma_CelBodyCatalogGUI or any of its controls.
function ma_CelBodyCatalogGUI_WindowKeyPressFcn(hObject, eventdata, handles)
    % hObject    handle to ma_CelBodyCatalogGUI (see GCBO)
    % eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
    %	Key: name of the key that was pressed, in lower case
    %	Character: character interpretation of the key(s) that was pressed
    %	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
    % handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'escape'
            close(handles.ma_CelBodyCatalogGUI);
    end


text:
function [rec,prec,ap] = VOCevaldet_modified(DATAopts, cls, loadName, draw, flipBoxes)

% load test set
tic;
cp=sprintf(DATAopts.annocachepath,DATAopts.testset);
if exist(cp,'file')
    fprintf('%s: pr: loading ground truth\n',cls);
    load(cp,'gtids','recs');
else
%     [gtids,t]=textread(sprintf(DATAopts.imgsetpath,DATAopts.testset),'%s %d');
    gtids = GetImagesPlusLabels(DATAopts.testset);
    for i=1:length(gtids)
        % display progress
        if toc>1
            fprintf('%s: pr: load: %d/%d\n',cls,i,length(gtids));
            drawnow;
            tic;
        end

        % read annotation
        recs(i)=PASreadrecord(sprintf(DATAopts.annopath,gtids{i}));
    end
    save(cp,'gtids','recs');
end

fprintf('%s: pr: evaluating detections\n',cls);

% hash image ids
hash=VOChash_init_modified(gtids);

% extract ground truth objects

npos=0;
gt(length(gtids))=struct('BB',[],'diff',[],'det',[]);
for i=1:length(gtids)
    % extract objects of class
    clsinds=strmatch(cls,{recs(i).objects(:).class},'exact');
    gt(i).BB=cat(1,recs(i).objects(clsinds).bbox)';
    gt(i).diff=[recs(i).objects(clsinds).difficult];
    gt(i).det=false(length(clsinds),1);
    npos=npos+sum(~gt(i).diff);
end

% load results
% [ids,confidence,b1,b2,b3,b4]=textread(sprintf(DATAopts.detrespath,id,cls),'%s %f %f %f %f %f');
[ids,confidence,b1,b2,b3,b4]=textread(loadName,'%s %f %f %f %f %f');

if exist('flipBoxes', 'var') && flipBoxes == true
    BB=[b2 b1 b4 b3]';
else
    BB=[b1 b2 b3 b4]';
end

% sort detections by decreasing confidence
[sc,si]=sort(-confidence);
ids=ids(si);
BB=BB(:,si);

% assign detections to ground truth objects
nd=length(confidence);
tp=zeros(nd,1);
fp=zeros(nd,1);
tic;
for d=1:nd
    % display progress
    if toc>1
        fprintf('%s: pr: compute: %d/%d\n',cls,d,nd);
        drawnow;
        tic;
    end

    % find ground truth image
    i=VOChash_lookup_modified(hash,ids{d});
    if isempty(i)
        error('unrecognized image "%s"',ids{d});
    elseif length(i)>1
        error('multiple image "%s"',ids{d});
    end

    % assign detection to ground truth object if any
    bb=BB(:,d);
    ovmax=-inf;
    for j=1:size(gt(i).BB,2)
        bbgt=gt(i).BB(:,j);
        bi=[max(bb(1),bbgt(1)) ; max(bb(2),bbgt(2)) ; min(bb(3),bbgt(3)) ; min(bb(4),bbgt(4))];
        iw=bi(3)-bi(1)+1;
        ih=bi(4)-bi(2)+1;
        if iw>0 & ih>0
            % compute overlap as area of intersection / area of union
            ua=(bb(3)-bb(1)+1)*(bb(4)-bb(2)+1)+...
               (bbgt(3)-bbgt(1)+1)*(bbgt(4)-bbgt(2)+1)-...
               iw*ih;
            ov=iw*ih/ua;
            if ov>ovmax
                ovmax=ov;
                jmax=j;
            end
        end
    end
    % assign detection as true positive/don't care/false positive
    if ovmax>=DATAopts.minoverlap
        if ~gt(i).diff(jmax)
            if ~gt(i).det(jmax)
                tp(d)=1;            % true positive
		gt(i).det(jmax)=true;
            else
                fp(d)=1;            % false positive (multiple detection)
            end
        end
    else
        fp(d)=1;                    % false positive
    end
end

% compute precision/recall
fp=cumsum(fp);
tp=cumsum(tp);
rec=tp/npos;
prec=tp./(fp+tp);

ap=VOCap(rec,prec);

if draw
    % plot precision/recall
    plot(rec,prec,'-');
    grid;
    xlabel 'recall'
    ylabel 'precision'
    title(sprintf('class: %s, subset: %s, AP = %.3f',cls,DATAopts.testset,ap));
end

author: nightrome

repo: matconvnet-calvin


text:
function point_num = sphere_icos_point_num ( factor )

%*****************************************************************************80
%
%% SPHERE_ICOS_POINT_NUM sizes an icosahedral grid on a sphere.
%
%  Discussion:
%
%    With FACTOR = 1, the grid has 20 triangular faces, 30 edges, and 12 nodes.
%
%    With FACTOR = 2, each triangle of the icosahedron is subdivided into
%    2x2 subtriangles, resulting in 80 faces, 120 edges, and
%    42 = 12 + 20 * 3 * (1)/2 + 20 * 0 ) nodes.
%
%    With FACTOR = 3, each triangle of the icosahedron is subdivided into
%    3x3 subtriangles, resulting in 180 faces, 270 edges, and
%    92 ( = 12 + 20 * 3 * (2)/2 + 20 * 1 ) nodes.
%
%    In general, each triangle is subdivided into FACTOR*FACTOR subtriangles,
%    resulting in 20 * FACTOR * FACTOR faces, 30 * FACTOR * FACTOR edges, and
%      12
%    + 20 * 3          * (FACTOR-1) / 2
%    + 20 * (FACTOR-2) * (FACTOR-1) / 2 nodes.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer FACTOR, the subdivision factor, which must
%    be at least 1.
%
%    Output, integer POINT_NUM, the number of nodes.
%
  point_num =   12                                   ...
              + 10 * 3              * ( factor - 1 ) ...
              + 10 * ( factor - 2 ) * ( factor - 1 );

  return
end


text:
function mrAnatConvertVAnatToT1Nifti(vAnat,outFileName)
%
%  mrAnatConvertVAnatToT1Nifti([vAnat],[outFileName])
%
% This function takes a vAnatomy.dat file converts it to a NIFTI and saves
% that NIFTI struct as vAnatomy.nii.gz by default.
%
% vAnat is a file name to a vAnatomy.dat
% outFilenmae defaults to 'vAnatomy.nii.gz'.
%
% Example:
%  vAnatPath = mrtInstallSampleData('anatomy/anatomyV','vAnatomy', ...
%   [], [], 'filetype', 'dat');
%  outFileName = [];
%  vAnatFileName = sprintf('%s.dat', vAnatPath);
%  mrAnatConvertVAnatToT1Nifti(vAnatFileName,outFileName);
%
% HISTORY:
% 4/20/2009 LMP wrote the thing
%
% LMP (c) Stanford VISTALAB

if(~exist('vAnat','var') || isempty(vAnat))
    [f,p] = uigetfile({'*.dat';'*.*'},'Select a vAnatomy.dat file for input...');
    if(isnumeric(f)), disp('User canceled.'); return; end
    vAnat = fullfile(p,f);
end

if(~exist('outFileName','var') || isempty(outFileName))
    p = fileparts(vAnat);
    outFileName = fullfile(p,'vAnatomy.nii.gz');
end

[vData,mmPerVox] = readVAnat(vAnat);

% Convert mrgray sagittal format to our preferred axial format for NIFTI
vData = flip(flip(permute(vData,[3 2 1]),2),3);
mmPerVox = mmPerVox([3 2 1]);
xform = [diag(1./mmPerVox), size(vData)'/2; 0 0 0 1];

% Modern version - BW, October 2012
ni = niftiCreate('fname',outFileName,'qto_xyz',inv(xform), 'data', vData);

% ni.qto_xyz = inv(xform);
% ni.qto_ijk = xform;
% ni.sto_xyz = inv(xform);
% ni.sto_ijk = xform;

niftiWrite(ni);

return

end

function [vData,mmPerVox] = readVAnat(vFile)

% set this to nan's in case the vAnatomy is old-style, in which
% case it won't have mmPerPix in the header.
mmPerVox = [nan,nan,nan];
volSize = [nan,nan,nan];

fid = fopen(vFile);
% read header: volSize and mmPerPix (if available)
tmp = fscanf(fid,'rows %f (%f mm/pixel)\n');
volSize(1) = tmp(1);
if length(tmp)>1, mmPerVox(1) = tmp(2); end
tmp = fscanf(fid,'cols %f (%f mm/pixel)\n');
volSize(2) = tmp(1);
if length(tmp)>1, mmPerVox(2) = tmp(2); end
tmp = fscanf(fid,'planes %f (%f mm/pixel)\n');
volSize(3) = tmp(1);
if length(tmp)>1, mmPerVox(3) = tmp(2); end

% Check that this is a valid header for a vAnatomy file. The next line should be '* \n'
endOfHeader = '*';
nextLine = fgets(fid);
if ~(length(nextLine)>=2 && nextLine(1)==endOfHeader)
    myErrorDlg(['vAnatomy file: ',vFile,' has invalid header']);
end

% If mmPerPix was not in the vAnatomy header, try to get it from UnfoldParams
if(isnan(mmPerVox(1)) || isnan(mmPerVox(2)) || isnan(mmPerVox(3)))
    % Try to get the voxel size from the old UnfoldParams.mat file.
    % warning('volume_pix_size not available. Try to get it from the old UnfoldParams file.');
    eval('ufp = load(fullfile(fpath,''UnfoldParams.mat''));', 'ufp = [];');
    if(isfield(ufp, 'volume_pix_size'))
        mmPerVox = 1./ufp.volume_pix_size;
    else
        % Error: can't find mmPerPix
        %myErrorDlg('Can not determine volume voxel size. You need to create an UnfoldParams file.');
        % Display a warning and carry on:
        disp('Warning : Can not determine volume voxel size. Setting to 1x1x1mm.');
        mmPerVox=[1 1 1];
    end
end

% read volume
vData = fread(fid,prod(volSize),'uint8');
fclose(fid);

% *** HACK!  Sometimes the vAnatomy is missing the last byte (???)
if length(vData) == prod(volSize)-1
    vData(end+1) = 0;
end

% Return vData permuted to maintain correct orientations. The old way was
% very inefficient.
vData=reshape(vData,[volSize(2),volSize(1),volSize(3)]);
vData=double(permute(vData,[2,1,3])); % This double cast is required for routines that expect readVolAnat to return a double.

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function headmodel = ft_headmodel_infinite(varargin)

% FT_HEADMODEL_INFINITE returns an infinitely large homogenous
% volume conduction model. For EEG the volume conductor can be used
% to compute the leadfield of electric current dipoles, for MEG it
% can be used for computing the leadfield of magnetic dipoles.
%
% Use as
%   headmodel = ft_headmodel_infinite;
%
% See also FT_PREPARE_VOL_SENS, FT_COMPUTE_LEADFIELD

% Copyright (C) 2012, Donders Centre for Cognitive Neuroimaging, Nijmegen, NL
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

sourcemodel = ft_getopt(varargin, 'sourcemodel', 'default');

% this is an easy one
headmodel = [];

switch sourcemodel
  case 'monopole'
    headmodel.type = 'infinite_monopole';
  case 'magneticdipole'
    headmodel.type = 'infinite_magneticdipole';
  case 'currentdipole'
    headmodel.type = 'infinite_currentdipole';
  case 'default'
    % let the EEG or MEG leadfield code decide
    headmodel.type = 'infinite';
end % switch


text:
function [ vertex_coordinate, vertex_label, edge_vertex, edge_label, ...
  triangle_vertex, triangle_label, quadrilateral_vertex, ...
  quadrilateral_label, tetrahedron_vertex, tetrahedron_label, ...
  hexahedron_vertex, hexahedron_label ] = mesh_data_read ( filename, dim, ...
  vertices, edges, triangles, quadrilaterals, tetrahedrons, hexahedrons )

%*****************************************************************************80
%
%% MESH_READ reads a MESH file defining a mesh and returns the data.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 November 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Pascal Frey,
%    MEDIT: An interactive mesh visualization software,
%    Technical Report RT-0253,
%    Institut National de Recherche en Informatique et en Automatique,
%    03 December 2001.
%
%  Parameters:
%
%    Input, string FILENAME, the name of the MESH file.
%
%    Input, integer DIM, the spatial dimension, which should be 2 or 3.
%
%    Input, integer VERTICES, the number of vertices.
%
%    Input, integer EDGES, the number of edges (may be 0).
%
%    Input, integer TRIANGLES, the number of triangles (may be 0).
%
%    Input, integer QUADRILATERALS, the number of quadrilaterals (may be 0).
%
%    Input, integer TETRAHEDRONS, the number of tetrahedrons (may be 0).
%
%    Input, integer HEXAHEDRONS, the number of hexahedrons (may be 0).
%
%    Output, real VERTEX_COORDINATE(DIM,VERTICES), the coordinates
%    of each vertex.
%
%    Output, integer VERTEX_LABEL(VERTICES), a label for each vertex.
%
%    Output, integer EDGE_VERTEX(2,EDGES), the vertices that form each edge.
%
%    Output, integer EDGE_LABEL(EDGES), a label for each edge.
%
%    Output, integer TRIANGLE_VERTEX(3,TRIANGLES), the vertices that form
%    each triangle.
%
%    Output, integer TRIANGLE_LABEL(TRIANGLES), a label for each triangle.
%
%    Output, integer QUADRILATERAL_VERTEX(4,QUADRILATERALS), the vertices that
%    form each quadrilateral.
%
%    Output, integer QUADRILATERAL_LABEL(QUADRILATERALS), a label for
%    each quadrilateral.
%
%    Output, integer TETRAHEDRON_VERTEX(4,TETRAHEDRONS), the vertices that
%    form each tetrahedron.
%
%    Output, integer TETRAHEDRON_LABEL(TETRAHEDRONS), a label for
%    each tetrahedron.
%
%    Output, integer HEXAHEDRON_VERTEX(8,HEXAHEDRONS), the vertices that form
%    each hexahedron.
%
%    Output, integer HEXAHEDRON_LABEL(HEXAHEDRONS), a label for each hexahedron.
%

%
%  Initialize everything to nothing.
%
  vertex_coordinate = [];
  vertex_label = [];
  edge_vertex = [];
  edge_label = [];
  triangle_vertex = [];
  triangle_label = [];
  quadrilateral_vertex = [];
  quadrilateral_label = [];
  tetrahedron_vertex = [];
  tetrahedron_label = [];
  hexahedron_vertex = [];
  hexahedron_label = [];
%
%  Open the file.
%
  unit = fopen ( filename, 'rt' );

  if ( unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MESH_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the input file "%s".\n', filename );
    error ( 'MESH_DATA_READ - Error!' );
    return
  end
%
%  Read lines til you get alphanumerics and determine a "mode"
%
  line_num = 0;
  keyword = 'NONE';

  while ( 1 )

    text = fgetl ( unit );

    if ( text == -1 )
      break
    end

    line_num = line_num + 1;

    if ( s_len_trim ( text ) == 0 )
      keyword = 'NONE';
      continue
    end

    if ( text(1) == '#' )
      continue
    end
%
%  Remove initial blanks.
%

%
%  Expecting a keyword.
%
        if ( s_eqi ( text, 'CORNERS' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'DIMENSION' ) )

      keyword = 'DIMENSION';

    elseif ( s_eqi ( text, 'EDGES' ) )

      keyword = 'EDGES';

    elseif ( s_eqi ( text, 'END' ) )

      fprintf ( 1, '\n' );
      fprintf ( 1, '  END statement encountered.\n' );
      break

    elseif ( s_eqi ( text, 'HEXAHEDRA' ) || s_eqi ( text, 'HEXAHEDRONS' ) )

      keyword = 'HEXAHEDRONS';

    elseif ( s_begin ( text, 'MESHVERSIONFORMATTED' ) )

    elseif ( s_eqi ( text, 'NORMALATQUADRILATERALVERTICES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'NORMALATTRIANGLEVERTICES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'NORMALATVERTICES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'NORMALS' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'QUADRILATERALS' ) )

      keyword = 'QUADRILATERALS';

    elseif ( s_eqi ( text, 'REQUIREDEDGES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'REQUIREDVERTICES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'RIDGES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'TANGENTATEDGES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'TANGENTS' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'TETRAHEDRA' ) || s_eqi ( text, 'TETRAHEDRONS' ) )

      keyword = 'TETRAHEDRONS';

    elseif ( s_eqi ( text, 'TRIANGLES' ) )

      keyword = 'TRIANGLES';

    elseif ( s_eqi ( text, 'VERTICES' ) )

      keyword = 'VERTICES';
%
%  Presumably, numeric data to be processed by keyword.
%
    elseif ( s_eqi ( keyword, 'DIMENSION' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
%       dim = value;
      end

      keyword = 'NONE';

    elseif ( s_eqi ( keyword, 'EDGES' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
%       edges = value;
      end

      keyword = 'EDGE_VERTEX';
      edge = 0;
      edge_vertex = zeros ( 2, edges );
      edge_label = zeros ( 1, edges );

    elseif ( s_eqi ( keyword, 'EDGE_VERTEX' ) )

      [ value, count ] = sscanf ( text, '%d  %d  %d' );

      edge = edge + 1;
      edge_vertex(1:2,edge) = value(1:2);
      edge_label(edge) = value(3);

    elseif ( s_eqi ( keyword, 'HEXAHEDRONS' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
%       hexahedrons = value;
      end

      keyword = 'HEXAHEDRON_VERTEX';
      hexahedron = 0;
      hexahedron_vertex = zeros ( 8, hexahedrons );
      hexahedron_label = zeros ( 1, hexahedrons );

    elseif ( s_eqi ( keyword, 'HEXAHEDRON_VERTEX' ) )

      [ value, count ] = sscanf ( text, '%d  %d  %d  %d  %d  %d  %d  %d  %d' );

      hexahedron = hexahedron + 1;
      hexahedron_vertex(1:8,hexahedron) = value(1:8);
      hexahedron_label(hexahedron) = value(9);

    elseif ( s_eqi ( keyword, 'QUADRILATERALS' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
%       quadrilaterals = value;
      end

      keyword = 'QUADRILATERAL_VERTEX';
      quadrilateral = 0;
      quadrilateral_vertex = zeros ( 4, quadrilateral );
      quadrilateral_label = zeros ( 1, quadrilateral );

    elseif ( s_eqi ( keyword, 'QUADRILATERAL_VERTEX' ) )

      [ value, count ] = sscanf ( text, '%d  %d  %d  %d  %d' );

      quadrilateral = quadrilateral + 1;
      quadrilateral_vertex(1:4,quadrilateral) = value(1:4);
      quadrilateral_label(quadrilateral) = value(5);

    elseif ( s_eqi ( keyword, 'TETRAHEDRONS' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
%       tetrahedrons = value;
      end

      keyword = 'TETRAHEDRON_VERTEX';
      tetrahedron = 0;
      tetrahedron_vertex = zeros ( 4, tetrahedron );
      tetrahedron_label = zeros ( 1, tetrahedron );

    elseif ( s_eqi ( keyword, 'TETRAHEDRON_VERTEX' ) )

      [ value, count ] = sscanf ( text, '%d  %d  %d  %d  %d' );

      tetrahedron = tetrahedron + 1;
      tetrahedron_vertex(1:4,tetrahedron) = value(1:4);
      tetrahedron_label(tetrahedron) = value(5);

    elseif ( s_eqi ( keyword, 'TRIANGLES' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
%       triangles = value;
      end

      keyword = 'TRIANGLE_VERTEX';
      triangle = 0;
      triangle_vertex = zeros ( 3, triangle );
      triangle_label = zeros ( 1, triangle );

    elseif ( s_eqi ( keyword, 'TRIANGLE_VERTEX' ) )

      [ value, count ] = sscanf ( text, '%d  %d  %d  %d' );

      triangle = triangle + 1;
      triangle_vertex(1:3,triangle) = value(1:3);
      triangle_label(triangle) = value(4);

    elseif ( s_eqi ( keyword, 'VERTICES' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
%       vertices = value;
      end

      keyword = 'VERTEX_COORDINATE';
      vertex = 0;
      vertex_coordinate = zeros ( dim, vertices );
      vertex_label = zeros ( 1, vertices );

    elseif ( s_eqi ( keyword, 'VERTEX_COORDINATE' ) )

      if ( dim == 2 )
        [ value, count ] = sscanf ( text, '%f  %f  %d' );
      elseif ( dim == 3 )
        [ value, count ] = sscanf ( text, '%f  %f  %f  %d' );
      end

      vertex = vertex + 1;
      vertex_coordinate(1:dim,vertex) = value(1:dim);
      vertex_label(vertex) = value(dim+1);

    elseif ( s_eqi ( keyword, 'SKIP' ) )

    else

      fprintf ( 1, '\n' );
      fprintf ( 1, 'MESH_DATA_READ - Fatal error!\n' );
      fprintf ( 1, '  Could not find keyword while reading line %d:\n', line_num );
      fprintf ( 1, '"%s"\n', text );
      error ( 'MESH_DATA_READ - Fatal error!\n' );

    end
  end
%
%  Close the file.
%
  fclose ( unit );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read %d lines from "%s".\n', line_num, filename );

  return
end


text:
function Ht = riskmetrics(data,lambda,backCast)
% Computes the Riskmetrics or other EWMA covariance
%
% USAGE:
%  [HT] = riskmetrics(DATA,LAMBDA,BACKCAST)
%
% INPUTS:
%   DATA     - A T by K matrix of zero mean residuals -OR-
%                K by K by T array of covariance estimators (e.g. realized covariance)
%   LAMBDA   - EWMA smoothing parameter 0<LAMBDA<1
%   BACKCAST - [OPTIONAL] Covariance matrix to use as initial value.  If
%                omitted, a backward EWMA is used
%
% OUTPUTS:
%   HT       - A [K K T] dimension matrix of conditional covariances
%
% COMMENTS:
%   The conditional variance, H(t), of an EWMA covariance model
%      H(t) = (1-lambda)*r(t-1,:)'*r(t-1,:) + lambda*H(t-1)
%
% EXAMPLES:
%   % The standard RiskMetrics EWMA covariance is computed using
%   Ht = riskmetrics(data,.94)
%   % The standard RiskMetrics EWMA covariance using the unconditional
%   % covariance as the backcast
%   Ht = riskmetrics(data,.94,cov(data))

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 03/10/2011

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch nargin
    case 2
        backCast = [];
    case 3
        % nothing
    otherwise
        error('2 or 3 inputs required.')
end

if ndims(data)==2
    [T,k] = size(data);
    temp = zeros(k,k,T);
    for t=1:T
        temp(:,:,t) = data(t,:)'*data(t,:);
    end
    data = temp;
end

if lambda<=0 || lambda>=1
    error('LAMBDA must be between 0 and 1.')
end

if isempty(backCast)
    endPoint = max(min(floor(log(.01)/log(lambda)),T),k);
    weights = (1-lambda).*lambda.^(0:endPoint-1);
    weights = weights/sum(weights);
    backCast = zeros(k);
    for i=1:endPoint
        backCast = backCast + weights(i)*data(:,:,i);
    end

end

backCast = (backCast+backCast)/2;
if min(eig(backCast))<0
    error('BACKCAST must be positive semidefinite if provided.')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ht = zeros(k,k,T);
Ht(:,:,1) = backCast;
for i=2:T
    Ht(:,:,i) = (1-lambda)*data(:,:,i-1) + lambda * Ht(:,:,i-1);
end

author: bashtage

repo: mfe-toolbox