text:
function [dS,t,f]=mtdspecgrampb(data,movingwin,phi,params)
% Multi-taper derivatives of time-frequency spectrum - binned point process
%
% Usage:
%
% [dS,t,f]=mtdspecgrampb(data,movingwin,phi,params)
% Input:
%   Note that all times can be in arbitrary units. But the units have to be
%   consistent. So, if E is in secs, win, t have to be in secs, and Fs has to
%   be Hz. If E is in samples, so are win and t, and Fs=1. In case of spike
%   times, the units have to be consistent with the units of data as well.
%       data        (in form samples x channels/trials or a single vector) -- required
%       movingwin         (in the form [window winstep] i.e length of moving
%                                                 window and step size.
%                                                 Note that units here have
%                                                 to be consistent with
%                                                 units of Fs
%       phi         (angle for evaluation of derivative) -- required.
%                       e.g. phi=[0,pi/2] giving the time and frequency
%                       derivatives
%       params: structure with fields tapers, pad, Fs, fpass,trialave
%       -optional
%           tapers : precalculated tapers from dpss or in the one of the following
%                    forms:
%                    (1) A numeric vector [TW K] where TW is the
%                        time-bandwidth product and K is the number of
%                        tapers to be used (less than or equal to
%                        2TW-1).
%                    (2) A numeric vector [W T p] where W is the
%                        bandwidth, T is the duration of the data and p
%                        is an integer such that 2TW-p tapers are used. In
%                        this form there is no default i.e. to specify
%                        the bandwidth, you have to specify T and p as
%                        well. Note that the units of W and T have to be
%                        consistent: if W is in Hz, T must be in seconds
%                        and vice versa. Note that these units must also
%                        be consistent with the units of params.Fs: W can
%                        be in Hz if and only if params.Fs is in Hz.
%                        The default is to use form 1 with TW=3 and K=5
%                    Note that T has to be equal to movingwin(1).
%
%	        pad		    (padding factor for the FFT) - optional (can take values -1,0,1,2...).
%                    -1 corresponds to no padding, 0 corresponds to padding
%                    to the next highest power of 2 etc.
%			      	 e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
%			      	 to 512 points, if pad=1, we pad to 1024 points etc.
%			      	 Defaults to 0.
%           Fs   (sampling frequency) - optional. Default 1.
%           fpass    (frequency band to be used in the calculation in the form
%                                   [fmin fmax])- optional.
%                                   Default all frequencies between 0 and
%                                   Fs/2
%           trialave (average over trials when 1, don't average when 0) -
%           optional. Default 0
% Output:
%       dS      (spectral derivative in form phi x time x frequency x channels/trials if trialave=0;
%               phi x time x frequency if trialave=1)
%       t       (times)
%       f       (frequencies)

if nargin < 3; error('Need data, window parameters and angle'); end;
if nargin < 4; params=[]; end;

if length(params.tapers)==3 & movingwin(1)~=params.tapers(2);
    error('Duration of data in params.tapers is inconsistent with movingwin(1), modify params.tapers(2) to proceed')
end

[tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
clear err
data=change_row_to_column(data);
[N,Ch]=size(data);
Nwin=round(Fs*movingwin(1)); % number of samples in window
Nstep=round(movingwin(2)*Fs); % number of samples to step through
nfft=max(2^(nextpow2(Nwin)+pad),Nwin);
f=getfgrid(Fs,nfft,fpass); Nf=length(f);
params.tapers=dpsschk(tapers,Nwin,Fs); % check tapers

winstart=1:Nstep:N-Nwin+1;
nw=length(winstart);
if trialave==0; dS=zeros(length(phi),nw,Nf,Ch); else dS=zeros(length(phi),nw,Nf); end;
for n=1:nw;
   indx=winstart(n):winstart(n)+Nwin-1;
   datawin=data(indx,:);
   [ds,f]=mtdspectrumpb(datawin,phi,params);
   dS(:,n,:,:)=ds;
end;
sz=size(ds);
dS=squeeze(dS);
% if length(sz)==3;
%    dS=permute(dS,[2 1 3 4]);
% elseif length(phi)>1
%    dS=permute(dS,[2 1 3]);
% end;
winmid=winstart+round(Nwin/2);
t=winmid/Fs;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function init_axes3d(axes_handle, grid3d, withinterp, withpml)
chkarg(~isempty(axes_handle) && ishandle(axes_handle), '"axes_handle" should be handle.');
chkarg(istypesizeof(grid3d, 'Grid3d'), '"grid3d" should be instance of Grid3d.');
chkarg(istypesizeof(withinterp, 'logical'), '"withinterp" should be logical.');
chkarg(istypesizeof(withpml, 'logical'), '"withpml" should be logical.');

hold(axes_handle, 'on');

if withinterp
	lplot = grid3d.lplot(GT.prim, withinterp, withpml);
else
	lplot = grid3d.lvoxelbound(GT.prim, withpml);
end
axis(axes_handle, [lplot{Axis.x}([1 end]), lplot{Axis.y}([1 end]), lplot{Axis.z}([1 end])]);
daspect(axes_handle, [1 1 1]);  % axis(ha, 'image');
box(axes_handle, 'on');
view(axes_handle, -38.5, 16)

hr = rotate3d;
set(hr, 'Enable', 'on');
setAllowAxesRotate(hr, axes_handle, true);

axes_str = {char(Axis.x), char(Axis.y), char(Axis.z)};
[like2d, n] = is2dlike(grid3d.N);
[h, v] = cycle(n);
h_str = [char(h), '  (', char(hex2dec('00D7'))];
if grid3d.unitvalue < 1e5 && grid3d.unitvalue > 1e-3
	h_str = [h_str, num2str(grid3d.unitvalue), ')'];
else
	h_str = [h_str, num2str(grid3d.unitvalue, '%.2e'), ')'];
end
axes_str{h} = h_str;

rotations = [0 0 0];
if like2d
	rotations(v) = -eps;  % small perturbation to show label correctly
end

% set(axes_handle, 'FontSize', 18, 'FontWeight', 'Bold');
xlabel(axes_handle, axes_str(Axis.x), 'Rotation', rotations(Axis.x));
ylabel(axes_handle, axes_str(Axis.y), 'Rotation', rotations(Axis.y));
zlabel(axes_handle, axes_str(Axis.z), 'Rotation', rotations(Axis.z));

set(axes_handle, 'TickDir', 'out');

% L0 = scalar3d.gi.L0;
% dlu = scalar3d.gi.display_length_unit;
% dlu_scale = scalar3d.gi.dlu_scale;
% xlabel(strcat(AxisName(Xx), ' (', char(hex2dec('00D7')), ...  % 00D7 is a unicode character for the multiplication sign.
% 	num2str(L0*dlu_scale), dlu, ')'));
% ylabel(strcat(AxisName(Yy), ' (', char(hex2dec('00D7')), ...  % 00D7 is a unicode character for the multiplication sign.
% 	num2str(L0*dlu_scale), dlu, ')'));
% zlabel(strcat(AxisName(Zz), ' (', char(hex2dec('00D7')), ...  % 00D7 is a unicode character for the multiplication sign.
% 	num2str(L0*dlu_scale), dlu, ')'));
% title(plot_title);

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function M = dimindex(A,dim,idx)

% DIMINDEX makes a selection from a multi-dimensional array where the dimension is
% selected by a scalar, not by the place between the brackets.
%
% Use as
%   M = dimindex(A,dim,idx)
%
% The purpose of the function is shown by the following example:
%
% A(:,:,:,23,:,:,...) is the same as dimindex(A,4,23)
% A(2,4,3)            is the same as dimindex(A,[1,2,3],[2,4,3])
% A(4,:,[5:10])       is the same as dimindex(A,[1,3],{4,[5:10]})
%
% See also the function DIMASSIGN

% Copyright (C) 2005, Geerten Kramer
%
% 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$

if (~iscell(idx))
  if (~any(size(dim)==1) || ~any(size(idx)==1) || ndims(dim)>2 || ndims(idx)>2 || length(dim)~=length(idx))
    ft_error('dim and idx must be both scalars or both vectors of the same size');
  end
  dummi = [];
  for i=1:length(idx)
    dummi{i} = idx(i);
  end
  idx = dummi;
  clear dummi;
end
if (~any(size(dim)==1) || ~any(size(idx)==1) || ndims(dim)>2 || ndims(idx)>2 || length(dim)~=length(idx))
  ft_error('dim and idx must be both scalars or both must have the same length');
end

if (length(dim)>ndims(A)||any(dim>ndims(A)))
  ft_error('dim, or one of its contents are larger than the number of dimentions in A');
end
if (~isequal(unique(dim),sort(dim)))
  ft_error('dim must be unique, every dimention can be addressed only once');
end

N = ndims(A);
for i=1:N
  ref = find(dim==i);
  if (isempty(ref))
    C{i} = ':';
  else
    C{i} = idx{ref};
  end
end
M = A(C{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef InstallDriverBase < matlab.hwmgr.internal.hwsetup.TemplateBase

%   Copyright 2017 The MathWorks, Inc.

  properties
    InfoText1
    InfoText2
    InfoText3
  end

  methods
    function obj = InstallDriverBase(varargin)
      obj@matlab.hwmgr.internal.hwsetup.TemplateBase(varargin{:});

      obj.Workflow.InstallDriverWorkflow = true;

      parentContainer = obj.ContentPanel;

      obj.InfoText1 = matlab.hwmgr.internal.hwsetup.Label.getInstance(parentContainer);
      obj.InfoText1.Position = [20 290 400 80];
      obj.InfoText1.Text = '';
      obj.InfoText1.Color = matlab.hwmgr.internal.hwsetup.util.Color.WHITE;

      obj.InfoText2 = matlab.hwmgr.internal.hwsetup.Label.getInstance(parentContainer);
      obj.InfoText2.Position = [20 230 400 50];
      obj.InfoText2.Text = '';
      obj.InfoText2.Color = matlab.hwmgr.internal.hwsetup.util.Color.WHITE;

      obj.InfoText3 = matlab.hwmgr.internal.hwsetup.Label.getInstance(parentContainer);
      obj.InfoText3.Position = [20 170 400 50];
      obj.InfoText3.Text = '';
      obj.InfoText3.Color = matlab.hwmgr.internal.hwsetup.util.Color.WHITE;

      obj.HelpText.AboutSelection = '';
      obj.HelpText.WhatToConsider = '';

      if strcmp(obj.Workflow.FirstScreenID, class(obj))
        obj.BackButton.Visible = 'off';
        obj.BackButton.Enable = 'off';
      end
    end

    function restoreScreen(obj)
      obj.enableScreen();
    end
  end

  methods (Access = protected, Abstract)
    status = installDriverImpl(obj)
  end

  methods (Access = protected)
    function installDriver(obj)
      obj.disableScreen();
      restoreOnCleanup = onCleanup(@obj.restoreScreen);

      p = matlab.hwmgr.internal.hwsetup.Panel.getInstance(obj.ContentPanel);
      p.Title = '';
      b = matlab.hwmgr.internal.hwsetup.BusyOverlay.getInstance(p);
      b.Text = 'Installing USB driver...';
      b.show();
      deleteOverlayOnCleanup = onCleanup(@()delete(p));

      [status,msg] = installDriverImpl(obj);

      if status
        error(message('plutoradio:hwsetup:DriverInstallFailed', msg));
      end

      refocus(obj.Workflow);
    end

    function out = getNextScreenIDImpl(obj)
      out = 'adi.ConnectHardware';
      installDriver(obj);
    end
  end
end

author: analogdevicesinc


text:
function X = MultiSolver(M, RHS)
% function X = MultiSolver(M, RHS)
%
% PURPOSE: Multiple linear (least-square) solver (usually small) for
%          systems having the same size
%
% INPUT
%   M  : 3D array (m x n x p)
%   RHS: 3D array (m x p x q)
% OUTPUT
%   X  : 3D array (n x p x q) if p>1,
%      : 2D array (n x q) if p==1 (squeezed, and equal to standard
%        backslash M\RHS)
%
% Solve p systems of linear equations:
%   M(:,:,k) * squeeze(X(:,k,:)) = squeeze(RHS(:,k,:)) for all k=1,2,...,p
%
% NOTE 1 -- This function is likely used with q == 1: single RHS for each
%           system). In this case
%       M  : 3D array (m x n x p)
%       RHS: 2D array (m x p)
%       X  : 2D array (n x p)
%       M(:,:,k) * X(:,k) = RHS(:,k) for all k=1,2,...,p
%
% NOTE 2 -- Special call with squeezed RHS: common RHS for all systems.
%   Input argument RHS can be squeezed as (m x q), or leave as 3D with
%   singleton in the second dimension: (m x 1 x q).
%       M   is 3D array (m x n x p)
%       RHS is 2D array (m x q) or (m x 1 x q)
%       X   is 2D array (n x p x q)
%       M(:,:,k) * squeeze(X(:,k,:)) = squeeze(RHS) for all k=1,2,...,p
%
% NOTE  3 -- Underdetermined system (more unknowns than equation)
%   The solution is basic solution obtained with sparse mldivide
%   which is not the same as basic solution when calling for full matrix.
%
% See also: SliceMultiSolver, MultiProd
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% History: original 26-May-2009
%          11-Auguts-2009, add help
% Acknowlegement: Tim Davis for the tip of using sparse to speedup
%                 from FOR loop method

if ndims(M)>3 || ndims(RHS)>3
    error('MultiSolver: M or RHS cannot have more than 3 dimensions');
end

[m n p] = size(M);

if size(RHS,1)~=m
    error('MultiSolver: M and RHS dimensions are not compatible');
end

if (size(RHS,2)~=p)
    if ndims(RHS)==2 % (m x q)
        RHS = reshape(RHS, m, 1, []);
    elseif size(RHS,2)~=1
        error('MultiSolver: M and RHS dimensions are not compatible');
    end
    % common RHS, expanded for all systems
    RHS = repmat(RHS, [1 p 1]); % (m x p x q)
end

RHS = reshape(RHS, m*p, []); % (m*p) x q
q = size(RHS,2);

% Build sparse matrix and solve
I = repmat(reshape(1:m*p,m,1,p),[1 n 1]); % m x n x p
J = repmat(reshape(1:n*p,1,n,p),[m 1 1]); % m x n x p
A = sparse(I(:),J(:),M(:));
X = reshape(A \ RHS, [n p q]);

% Squeeze for single system
if p==1
    X = reshape(X, [n q]);
end

end % MultiSolver

author: Sable

repo: mcbench-benchmarks


text:
function p = cvx_create_problem( varargin )
global cvx___

error( nargchk( 0, 2, nargin ) );
cvx_global
cvx_problem = evalin( 'caller', 'cvx_problem', '[]' );
if isa( cvx_problem, 'cvxprob' ) && cvx___.problems( index( cvx_problem ) ).depth == length( dbstack ) - 2,
    error( 'A cvx problem already exists in this scope.\n(To clear it and start a new one, use the command ''cvx_clear''.)', 1 ); %#ok
end
cvx_problem = cvxprob;
if nargin > 0,
    p = index( cvx_problem );
    for k = 1 : nargin,
        mode = varargin{k};
        if isempty( mode ),
            continue;
        elseif ~ischar( mode ) || size( mode, 1 ) ~= 1,
            cvx_pop( cvx_problem, 'clear' );
            error( 'Arguments must be strings.' );
        end
        switch lower( mode ),
            case 'quiet',
                cvx___.problems( p ).quiet = true;
            case 'set',
                cvx___.problems( p ).complete  = false;
                cvx___.problems( p ).direction = 'find';
            case 'sdp',
                if cvx___.problems( p ).gp,
                    cvx_pop( cvx_problem, 'clear' );
                    error( 'The GP and SDP modifiers cannot be used together.' );
                end
                cvx___.problems( p ).sdp = true;
            case 'gp',
                if cvx___.problems( p ).sdp,
                    cvx_pop( cvx_problem, 'clear' );
                    error( 'The GP and SDP modifiers cannot be used together.' );
                end
                cvx___.problems( p ).gp = true;
            case 'separable',
                cvx___.problems( p ).separable = true;
        end
    end
end
assignin( 'caller', 'cvx_problem', cvx_problem );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function res = vect_isPointInMesh(point, v, f, varargin)
%ISPOINTINMESH Check if a point is inside a 3D mesh.
%
%   B = isPointInMesh(PT, V, F)
%   Check if the point PT (given as a 1-by-3 array) is inside the mesh
%   defined by the vertices V and the face array F. The result is a
%   boolean.
%
%   If PT is a N-by-3 point array, the result is a N-by-1 array of logical.
%
%   Example
%     [v, f] = torusMesh([50 50 50 30 10 30 45]);
%     [x, y, z] = meshgrid(5:5:100, 5:5:100, 5:5:100);
%     res = false(size(x));
%     res(:) = isPointInMesh([x(:) y(:) z(:)], v, f);
%     figure; plot3(x(res), y(res), z(res), 'b.'); axis equal;
%
%   Algorithm:
%   The method computes the intersection with a ray starting from the
%   point(s) and with a random orientation. Some errors are possible if
%   rays crosses the mesh between two or three faces.
%
%   See also
%     meshes3d, intersectLineMesh3d

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2018-01-26,    using Matlab 9.3.0.713579 (R2017b)
% Copyright 2018 INRA - Cepia Software Platform.

% choose a random vector
vect = rand(1, 3);

% initialize array for result
np = size(point, 1);
res = false(np, 1);

% iterate over the various points
lines = createLine3d(point, vect);
[~, pos] = vect_intersectLineMesh3d(lines, v, f);
for i = 1:np
    res(i) = mod(sum(pos(:,:,i) > 0), 2) > 0;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function launch = ma_createLaunch(name, launchEpoch, launchLat, launchLong, launchAlt, ...
                                        launchToF, burnoutLat, burnoutLong, burnoutAlt, ...
                                        lineColor, lineStyle, lineWidth, bodyId)

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % vars = [bool bool bool bool bool bool bool bool;
    %         lower lower lower lower lower lower lower lower;
    %         upper upper upper upper upper upper upper upper]
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    launch = struct();
    launch.name = name;
    launch.type = 'Launch';
    launch.bodyId = bodyId;
    launch.launchValue = [launchEpoch, launchLat, launchLong, launchAlt, launchToF, burnoutLat, burnoutLong, burnoutAlt];
    launch.lineColor = lineColor;
    launch.lineStyle = lineStyle;
    launch.lineWidth = lineWidth;
    launch.id = rand(1);
end


text:
%                           solveWirtFlow.m
%
%  Implementation of the Wirtinger Flow (WF) algorithm proposed in the
%  paper.
%
%% I/O
%  Inputs:
%     A:    m x n matrix or a function handle to a method that
%           returns A*x.
%     At:   The adjoint (transpose) of 'A'. If 'A' is a function handle,
%           'At' must be provided.
%     b0:   m x 1 real,non-negative vector consists of all the measurements.
%     x0:   n x 1 vector. It is the initial guess of the unknown signal x.
%     opts: A struct consists of the options for the algorithm. For details,
%           see header in solvePhaseRetrieval.m or the User Guide.
%
%     Note: When a function handle is used, the
%     value of 'At' (a function handle for the adjoint of 'A') must be
%     supplied.
%
%  Outputs:
%     sol:  n x 1 vector. It is the estimated signal.
%     outs: A struct consists of the convergence info. For details,
%           see header in solvePhaseRetrieval.m or the User Guide.
%
%
%  See the script 'testWirtFlow.m' for an example of proper usage of this
%  function.
%
%% Notations
%  x is the estimation of the signal. y is the vector of measurements such
%  that yi = |<ai,x>|^2 for i = 1,...,m
%
%% Algorithm Description
%  WF successively refines the estimate via an update rule that bears a
%  strong resemblance to a gradient descent scheme. Specifically, at each
%  iteration, x = x + mu/m * gradient log-likelihood of x given y For the
%  detailed formulation of "gradient log-likelihood of x given y" and a
%  detailed explanation of the theory, see the WF paper referenced below.
%
%% References
%  Paper Title:   Phase Retrieval via Wirtinger Flow: Theory and Algorithms
%  Place:         Chapter 2.3
%  Authors:       Emmanuel Candes, Xiaodong Li, Mahdi Soltanolkotabi
%  arXiv Address: https://arxiv.org/abs/1407.1065
%
% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

%% -----------------------------START-----------------------------------

function [sol, outs] = solveWirtFlow(A, At, b0, x0, opts)
    innerOpts = struct;
    innerOpts.maxIters = opts.maxIters;
    innerOpts.maxTime = opts.maxTime;
    innerOpts.tol = opts.tol;
    innerOpts.verbose = opts.verbose;
    innerOpts.recordTimes = opts.recordTimes;
    innerOpts.recordResiduals = opts.recordResiduals;
    innerOpts.recordMeasurementErrors = opts.recordMeasurementErrors;
    innerOpts.recordReconErrors = opts.recordReconErrors;
    innerOpts.xt = opts.xt;

    innerOpts.searchMethod = opts.searchMethod;
    innerOpts.betaChoice = opts.betaChoice;

    [sol, outs] = gradientDescentSolver(A, At, x0, b0, @updateObjective, innerOpts);

    function [f, gradf] = updateObjective(~, ~)
        f = @(z) 0.5 * norm(abs(z).^2 - b0.^2)^2;
        gradf = @(z) (abs(z).^2 - b0.^2) .* z;
    end
end

author: tomgoldstein

repo: phasepack-matlab


text:
% compile_mex   Compiles mex files required for hctsa package
%
% This script must be run in the Toolboxes directory.
%
% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
% Check we're in the correct directory
% ------------------------------------------------------------------------------

currentDir = pwd;

% Path split using platform-dependent separator
weHere = regexp(currentDir,filesep,'split');

if ~strcmp(weHere{end},'Toolboxes')
    error('This code must be run in the ''Toolboxes'' directory of the HCTSA package...')
end

% Sweet. Toolbox path is:
toolDir = currentDir;

% ------------------------------------------------------------------------------
% Max Little's fastdfa code
% ------------------------------------------------------------------------------
fprintf(1,'fastdfa...');
try
    cd(fullfile(toolDir,'Max_Little','fastdfa'));
	mex ML_fastdfa_core.c
    fprintf(1,' done.\n');
catch emsg
    fprintf(1,'%s\n\n',emsg.message);
    cd(toolDir);
    errMsg = sprintf(['An error occurred while compiling ML_Fastdfa_core C code.\n' ...
        'It appears that mex is not set up to work on this system (cf. ''doc mex'' and ''mex -setup'').\n' ...
        'Get ''mex ML_fastdfa_core.c'' to work, and then re-run compile_mex.m']);
    error(errMsg)
end

% ------------------------------------------------------------------------------
% Max Little's Steps Bumps Toolkit
% ------------------------------------------------------------------------------
fprintf(1,'Max Little''s ''Steps and bumps'' toolkit...');
cd(fullfile(toolDir,'Max_Little','steps_bumps_toolkit'))
anyErrors = false;
try
    mex ML_kvsteps_core.cpp
catch
    fprintf(1,'ERROR: Max Little''s ''Steps and bumps'' C++ code failed to compile correctly.\n');
end
if ~anyErrors
    fprintf(1,' done.\n');
end

% ------------------------------------------------------------------------------
% Max Little's RPDE toolkit
% ------------------------------------------------------------------------------
fprintf(1,'Max Little''s ''RPDE'' code...');
cd(fullfile(toolDir,'Max_Little','rpde'))
anyErrors = false;
try
    mex ML_close_ret.c
catch
    fprintf(1,'ERROR: Max Little''s ''RPDE'' C code failed to compile correctly\n');
end
if ~anyErrors
    fprintf(1,' done.\n');
end

% ------------------------------------------------------------------------------
% Michael Small's code
% ------------------------------------------------------------------------------
fprintf(1,'Michael Small''s code...');
cd(fullfile(toolDir,'Michael_Small'))
anyErrors = false;
try
    mex MS_complexitybs.c % compile Michael Small's complexitybs C code
catch
    fprintf(1,'ERROR: Michael Small''s ''complexitybs'' C code failed to compile correctly\n');
end
try
    mex MS_nearest.c      % compile Michael Small's nearest C code
catch
    fprintf(1,'ERROR: Michael Small''s ''nearest'' C code failed to compile correctly\n');
end
try
    mex MS_shannon.c      % compile Michael Small's shannon C code
catch
    fprintf(1,'ERROR: Michael Small''s ''shannon'' C code failed to compile correctly\n');
end
if ~anyErrors
    fprintf(1,' done.\n');
end

% ------------------------------------------------------------------------------
% Gaussian Process code, gpml
% ------------------------------------------------------------------------------
fprintf(1,'Gaussian Process Toolbox, Carl Edward Rasmussen and Hannes Nickisch...');
cd(fullfile(toolDir,'gpml','util'))
anyErrors = false;
try
    make
catch
    fprintf(1,'ERROR: Gaussian Process Toolbox failed to compile correctly.\n');
end
if ~anyErrors
    fprintf(1,' done.\n');
end

%-------------------------------------------------------------------------------
% Physionet sample entropy code (turned to mex)
%-------------------------------------------------------------------------------
fprintf(1,'Sample entropy...');
cd(fullfile(toolDir,'Physionet'))
anyErrors = false;
try
    mex sampen_mex.c
catch
    fprintf(1,'ERROR: Physionet implementation of sample entropy failed to compile.\n');
end
if ~anyErrors
    fprintf(1,' done.\n');
end

% ------------------------------------------------------------------------------
% TSTOOL routines (such a mess)
% ------------------------------------------------------------------------------
fprintf(1,'TSTOOL...');
cd(fullfile(toolDir,'OpenTSTOOL','mex-dev'))
anyErrors = false;
try
    makemex
    settspath(fullfile(toolDir,'OpenTSTOOL'));
catch emsg
    fprintf(1,'---ERROR: %s\n',emsg.message);
    fprintf(1,'The TSTOOL package failed to compile correctly.\n');
end
if ~anyErrors
    fprintf(1,'Done!\n');
end

% ------------------------------------------------------------------------------
% TISEAN
% ------------------------------------------------------------------------------
fprintf(1,['NB: To use TISEAN routines, you need to compile and install them on your system\n' ...
    'In the commandline, navigate to the Toolboxes%sTisean_3.0.1 directory\n  ''./configure'',' ...
    ' ''make'' and ''make install'' commands\n(cf. Documentation for instructions,' ...
    ' including implementations for Windows)...\n'],filesep);

%-------------------------------------------------------------------------------
% CATCH22
%-------------------------------------------------------------------------------
cd(fullfile(toolDir,'catch22','wrap_Matlab'))
mexAll

% Return to base directory
cd(toolDir);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
classdef TP10 < PROBLEM
% <multi> <real> <large/none> <constrained> <robust>
% Test problem for robust multi-objective optimization
% delta --- 0.05 --- Maximum disturbance degree
% H     ---   50 --- Number of disturbances

%------------------------------- Reference --------------------------------
% A. Gaspar-Cunha, J. Ferreira, and G. Recio, Evolutionary robustness
% analysis for multi-objective optimization: benchmark problems, Structural
% and Multidisciplinary Optimization, 2014, 49: 771-793.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        delta;      % Maximum disturbance degree
        H;          % Number of disturbances
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            [obj.delta,obj.H] = obj.ParameterSet(0.05,50);
            obj.M = 2;
            obj.D = 3;
            obj.lower    = [0,0,1];
            obj.upper    = [10,10,3];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(~,PopDec)
            PopObj(:,1) = PopDec(:,1).*sqrt(16+PopDec(:,3).^2) + PopDec(:,2).*sqrt(1+PopDec(:,3).^2);
            PopObj(:,2) = 20*sqrt(16+PopDec(:,3).^2)./PopDec(:,1)./PopDec(:,3);
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            PopCon(:,1) = 20*sqrt(16+PopDec(:,3).^2) - 100*PopDec(:,1).*PopDec(:,3);
            PopCon(:,2) = 80*sqrt(1+PopDec(:,3).^2) - 100*PopDec(:,2).*PopDec(:,3);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = [100,100];
        end
        %% Calculate the metric value
        function score = CalMetric(obj,metName,Population)
            switch metName
                case {'Mean_IGD','Mean_HV','Worst_IGD','Worst_HV'}
                    score = feval(metName,Population,obj);
                otherwise
                    score = feval(metName,Population,obj.optimum);
            end
        end
        %% Perturb solutions multiple times
        function PopX = Perturb(obj,PopDec,N)
            if nargin < 3; N = obj.H; end
            Delta = repmat(obj.delta.*(obj.upper-obj.lower),N*size(PopDec,1),1);
            w     = UniformPoint(N,obj.D,'Latin');
            Dec   = 2*Delta.*w(reshape(repmat(1:end,size(PopDec,1),1),1,[]),:) + repmat(PopDec,N,1) - Delta;
            Dec   = obj.CalDec(Dec);
            PopX  = SOLUTION(Dec,obj.CalObj(Dec),obj.CalCon(Dec));
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [ s, seed ] = polygon_sample ( nv, v, n, seed )

%*****************************************************************************80
%
%% POLYGON_SAMPLE uniformly samples a polygon.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NV, the number of vertices.
%
%    Input, real V(2,NV), the vertices of the polygon, listed in
%    counterclockwise order.
%
%    Input, integer N, the number of points to create.
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Output, real S(2,N), the points.
%

%
%  Triangulate the polygon.
%
  triangles = polygon_triangulate ( nv, v(1,1:nv), v(2,1:nv) );
%
%  Determine the areas of each triangle.
%
  for i = 1 : nv - 2
    area_triangle(i) = triangle_area ( ...
      v(1,triangles(1,i)), v(2,triangles(1,i)), ...
      v(1,triangles(2,i)), v(2,triangles(2,i)), ...
      v(1,triangles(3,i)), v(2,triangles(3,i)) );
  end
%
%  Normalize the areas.
%
  area_polygon = sum ( area_triangle(1:nv-2) );
  area_relative(1:nv-2) = area_triangle(1:nv-2) / area_polygon;
%
%  Replace each area by the sum of itself and all previous ones.
%
  area_cumulative(1) = area_relative(1);
  for i = 2 : nv - 2
    area_cumulative(i) = area_relative(i) + area_cumulative(i-1);
  end

  for j = 1 : n
%
%  Choose triangle I at random, based on areas.
%
    [ area_percent, seed ] = r8_uniform_01 ( seed );

    for k = 1 : nv - 2

      i = k;

      if ( area_percent <= area_cumulative(k) )
        break
      end

    end
%
%  Now choose a point at random in triangle I.
%
    [ r, seed ] = r8vec_uniform_01 ( 2, seed );

    if ( 1.0 < sum ( r(1:2) ) )
      r(1:2) = 1.0 - r(1:2);
    end

    s(1:2,j) = ( 1.0 - r(1) - r(2) ) * v(1:2,triangles(1,i)) ...
                     + r(1)          * v(1:2,triangles(2,i)) ...
                            + r(2)   * v(1:2,triangles(3,i));
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function b = symm_random_inverse ( n, d, key )

%*****************************************************************************80
%
%% SYMM_RANDOM_INVERSE returns the inverse of the SYMM_RANDOM matrix.
%
%  Discussion:
%
%    The input value of SEED must be the same as the input value used
%    when defining the original matrix.
%
%    For the inverse to exist, no entry of D may be zero.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real D(N), the desired eigenvalues for the matrix.
%
%    Input, integer KEY, a positive integer that selects the data.
%
%    Output, real B(N,N), the matrix.
%

%
%  Get a random orthogonal matrix Q.
%
  q = orth_random ( n, key );
%
%  Set B = Q * 1/Lambda * Q'.
%
  l_inv = diag ( 1.0 ./ d );

  b = q * l_inv * q';

  return
end


text:
function v = conj(v)
%CONJ Complex conjugate of a BALLFUNV.
%   CONJ(F) returns the complex conjugate of F. For a complex F, CONJ(F) =
%   REAL(F) - i*IMAG(F).
%
%   See also IMAG, REAL.

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

if isempty( v )
    return
end

% Take conjugate part of each component:
v = ballfunv( conj(v.comp{1}), conj(v.comp{2}), conj(v.comp{3}) );
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function gpcf = gpcf_mask(varargin)
%GPCF_MASK  Create a mask covariance function
%
%  Description
%    GPCF = GPCF_MASK('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates a mask covariance function structure in
%    which the named parameters have the specified values. Any
%    unspecified parameters are set to default values.
%
%    GPCF = GPCF_MASK(GPCF,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a covariance function structure with the named
%    parameters altered with the specified values.
%
%    Mask covariance function returns correlation 1 if input
%    values X_i and X_j are both non zero and 0 otherwise.
%
%    Parameters for mask covariance function
%      selectedVariables - vector defining which inputs are used
%
%  See also
%    GP_SET, GPCF_*, PRIOR_*, MEAN_*
%
% Copyright (c) 2007-2010 Jarno Vanhatalo
% Copyright (c) 2008-2010 Jaakko Riihimäki
% Copyright (c) 2010 Aki Vehtari

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

  ip=inputParser;
  ip.FunctionName = 'GPCF_MASK';
  ip.addOptional('gpcf', [], @isstruct);
  ip.addParamValue('selectedVariables',[], ...
                   @(x) isempty(x) || (isvector(x) && all(x>0)));
  ip.parse(varargin{:});
  gpcf=ip.Results.gpcf;

  if isempty(gpcf)
    % Initialize a covariance function structure
    init=true;
    gpcf.type = 'gpcf_mask';
  else
    % Modify a covariance function structure
    if ~isfield(gpcf,'type') && ~isequal(gpcf.type,'gpcf_mask')
      error('First argument does not seem to be a valid covariance function structure')
    end
    init=false;
  end

  if ~ismember('selectedVariables',ip.UsingDefaults)
    if ~isempty(ip.Results.selectedVariables)
      gpcf.selectedVariables = ip.Results.selectedVariables;
    elseif isfield(gpcf,'selectedVariables')
      gpcf=rmfield(gpcf,'selectedVariables');
    end
  end
  if init
    gpcf.fh.pak = @gpcf_mask_pak;
    gpcf.fh.unpak = @gpcf_mask_unpak;
    gpcf.fh.lp = @gpcf_mask_lp;
    gpcf.fh.lpg = @gpcf_mask_lpg;
    gpcf.fh.cfg = @gpcf_mask_cfg;
    gpcf.fh.ginput = @gpcf_mask_ginput;
    gpcf.fh.cov = @gpcf_mask_cov;
    gpcf.fh.trcov  = @gpcf_mask_trcov;
    gpcf.fh.trvar  = @gpcf_mask_trvar;
    gpcf.fh.recappend = @gpcf_mask_recappend;
  end

end

function [w,s,h] = gpcf_mask_pak(gpcf, w)
%GPCF_MASK_PAK  Combine GP covariance function parameters into
%              one vector.
%
%  Description
%    W = GPCF_MASK_PAK(GPCF) takes a covariance function
%    structure GPCF and combines the covariance function
%    parameters and their hyperparameters into a single row
%    vector W. This is a mandatory subfunction used for
%    example in energy and gradient computations.
%
%       w = []
%
%  See also
%    GPCF_MASK_UNPAK

  w = []; s = {}; h = [];
end

function [gpcf, w] = gpcf_mask_unpak(gpcf, w)
%GPCF_MASK_UNPAK  Sets the covariance function parameters into
%                the structure
%
%  Description
%    [GPCF, W] = GPCF_MASK_UNPAK(GPCF, W) takes a covariance
%    function structure GPCF and a parameter vector W, and
%    returns a covariance function structure identical to the
%    input, except that the covariance parameters have been set
%    to the values in W. Deletes the values set to GPCF from W
%    and returns the modified W. This is a mandatory subfunction
%    used for example in energy and gradient computations.
%
%    Assignment is inverse of
%       w = []
%
%  See also
%   GPCF_MASK_PAK

end

function lp = gpcf_mask_lp(gpcf)
%GPCF_MASK_LP  Evaluate the energy of prior of covariance function parameters
%
%  Description
%    LP = GPCF_MASK_LP(GPCF) takes a covariance function
%    structure GPCF and returns log(p(th)), where th collects the
%    parameters. This is a mandatory subfunction used for example
%    in energy computations.
%
%  See also
%    GPCF_MASK_PAK, GPCF_MASK_UNPAK, GPCF_MASK_LPG, GP_E

  lp = 0;

end

function lpg = gpcf_mask_lpg(gpcf)
%GPCF_MASK_LPG  Evaluate gradient of the log prior with respect
%               to the parameters.
%
%  Description
%    LPG = GPCF_MASK_LPG(GPCF) takes a covariance function
%    structure GPCF and returns LPG = d log (p(th))/dth, where th
%    is the vector of parameters. This is a mandatory subfunction
%    used for example in gradient computations.
%
%  See also
%    GPCF_MASK_PAK, GPCF_MASK_UNPAK, GPCF_MASK_LP, GP_G

  lpg = [];

end

function DKff = gpcf_mask_cfg(gpcf, x, x2, mask,i1)
%GPCF_MASK_CFG  Evaluate gradient of covariance function
%               with respect to the parameters.
%
%  Description
%    DKff = GPCF_MASK_CFG(GPCF, X) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X) with
%    respect to th (cell array with matrix elements). This is a
%    mandatory subfunction used in gradient computations.
%
%    DKff = GPCF_MASK_CFG(GPCF, X, X2) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X2) with
%    respect to th (cell array with matrix elements). This subfunction
%    is needed when using sparse approximations (e.g. FIC).
%
%    DKff = GPCF_MASK_CFG(GPCF, X, [], MASK) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the diagonal of gradients of covariance matrix
%    Kff = k(X,X2) with respect to th (cell array with matrix
%    elements). This subfunction is needed when using sparse
%    approximations (e.g. FIC).
%
%    DKff = GPCF_MASK_CFG(GPCF, X, X2, [], i) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X2), or k(X,X)
%    if X2 is empty, with respect to ith hyperparameter. This subfunction
%    is needed when using memory save option in gp_set.
%
%  See also
%    GPCF_MASK_PAK, GPCF_MASK_UNPAK, GPCF_MASK_LP, GP_G

  DKff = {};

end

function [DKff, lpg]  = gpcf_mask_ginput(gpcf, x, x2, i1)
%GPCF_MASK_GINPUT  Evaluate gradient of covariance function with
%                  respect to x.
%
%  Description
%    DKff = GPCF_MASK_GINPUT(GPCF, X) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X) with
%    respect to X (cell array with matrix elements). This subfunction
%    is needed when computing gradients with respect to inducing
%    inputs in sparse approximations.
%
%    DKff = GPCF_MASK_GINPUT(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors
%    and returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2) with respect to X (cell array with matrix elements).
%    This subfunction is needed when computing gradients with
%    respect to inducing inputs in sparse approximations.
%
%    DKff = GPCF_MASK_GINPUT(GPCF, X, X2, i) takes a covariance
%    function structure GPCF, a matrix X of input vectors
%    and returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2), or k(X,X) if X2 is empty, with respect to ith covariate
%    in X. This subfunction is needed when using memory saving option
%    in gp_set.
%
%  See also
%   GPCF_MASK_PAK, GPCF_MASK_UNPAK, GPCF_MASK_LP, GP_G

  [n, m] =size(x);

  if nargin == 2 || isempty(x2)
    ii1 = 0;
    for i=1:m
      for j = 1:n
        ii1 = ii1 + 1;
        DKff{ii1} = zeros(n);
        lpg(ii1) = 0;
      end
    end
  elseif nargin == 3 || nargin == 4
    ii1 = 0;
    for i=1:m
      for j = 1:n
        ii1 = ii1 + 1;
        DKff{ii1} = zeros(n, size(x2,1));
        lpg(ii1) = 0;
      end
    end
  end
  if nargin==4
    DKff=DKff{1};
  end
end

function C = gpcf_mask_cov(gpcf, x1, x2, varargin)
%GP_MASK_COV  Evaluate covariance matrix between two input vectors
%
%  Description
%    C = GP_MASK_COV(GP, TX, X) takes in covariance function of a
%    Gaussian process GP and two matrixes TX and X that contain
%    input vectors to GP. Returns covariance matrix C. Every
%    element ij of C contains covariance between inputs i in TX
%    and j in X. This is a mandatory subfunction used for example in
%    prediction and energy computations.
%
%  See also
%    GPCF_MASK_TRCOV, GPCF_MASK_TRVAR, GP_COV, GP_TRCOV

  if isempty(x2)
    x2=x1;
  end
  [n1,m1]=size(x1);
  [n2,m2]=size(x2);

  if m1~=m2
    error('the number of columns of X1 and X2 has to be same')
  end

  C=repmat(true,n1,n2);
  if isfield(gpcf, 'selectedVariables')
    for j = 1:length(gpcf.selectedVariables)
      jj=gpcf.selectedVariables(j);
      C = C & bsxfun(@and,x1(:,jj)~=0,x2(:,jj)'~=0);
    end
  else
    for j = 1:m1
      C = C & bsxfun(@and,x1(:,j)~=0,x2(:,j)'~=0);
    end
  end
  C=double(C);

end

function C = gpcf_mask_trcov(gpcf, x)
%GP_MASK_TRCOV  Evaluate training covariance matrix of inputs
%
%  Description
%    C = GP_MASK_TRCOV(GP, TX) takes in covariance function of a
%    Gaussian process GP and matrix TX that contains training
%    input vectors. Returns covariance matrix C. Every element ij
%    of C contains covariance between inputs i and j in TX. This
%    is a mandatory subfunction used for example in prediction
%    and energy computations.
%
%  See also
%    GPCF_MASK_COV, GPCF_MASK_TRVAR, GP_COV, GP_TRCOV

  [n,m]=size(x);

  C=repmat(true,n,n);
  if isfield(gpcf, 'selectedVariables')
    for j = 1:length(gpcf.selectedVariables)
      jj=gpcf.selectedVariables(j);
      C = C & bsxfun(@and,x(:,jj)~=0,x(:,jj)'~=0);
    end
  else
    for j = 1:m
      C = C & bsxfun(@and,x(:,j)~=0,x(:,j)'~=0);
    end
  end
  C=double(C);

end

function C = gpcf_mask_trvar(gpcf, x)
%GP_MASK_TRVAR  Evaluate training variance vector
%
%  Description
%    C = GP_MASK_TRVAR(GPCF, TX) takes in covariance function of a
%    Gaussian process GPCF and matrix TX that contains training
%    inputs. Returns variance vector C. Every element i of C
%    contains variance of input i in TX. This is a mandatory
%    subfunction used for example in prediction and energy
%    computations.
%
%  See also
%    GPCF_MASK_COV, GP_COV, GP_TRCOV

  [n,m]=size(x);
  C=true(n,1);
  if isfield(gpcf, 'selectedVariables')
    for j = 1:length(gpcf.selectedVariables)
      jj=gpcf.selectedVariables(j);
      C = C & x(:,jj)~=0;
    end
  else
    for j = 1:m
      C = C & x(:,j)~=0;
    end
  end
  C=double(C);

end

function reccf = gpcf_mask_recappend(reccf, ri, gpcf)
%RECAPPEND  Record append
%
%  Description
%    RECCF = GPCF_MASK_RECAPPEND(RECCF, RI, GPCF) takes a
%    covariance function record structure RECCF, record index RI
%    and covariance function structure GPCF with the current MCMC
%    samples of the parameters. Returns RECCF which contains
%    all the old samples and the current samples from GPCF.
%    This subfunction is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC and GP_MC -> RECAPPEND

  if nargin == 2
    % Initialize the record
    reccf.type = 'gpcf_mask';

    % Initialize parameters
    reccf.coeffSigma2= [];

    % Set the function handles
    reccf.fh.pak = @gpcf_mask_pak;
    reccf.fh.unpak = @gpcf_mask_unpak;
    reccf.fh.lp = @gpcf_mask_lp;
    reccf.fh.lpg = @gpcf_mask_lpg;
    reccf.fh.cfg = @gpcf_mask_cfg;
    reccf.fh.cov = @gpcf_mask_cov;
    reccf.fh.trcov  = @gpcf_mask_trcov;
    reccf.fh.trvar  = @gpcf_mask_trvar;
    reccf.fh.recappend = @gpcf_mask_recappend;
  else
    % Append to the record
    if isfield(gpcf, 'selectedVariables')
      reccf.selectedVariables = gpcf.selectedVariables;
    end
  end
end

author: gpstuff-dev

repo: gpstuff

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


text:
%CODINGCONVENTIONS specify coding conventions for MatGeom files
%
%   General
%
%   This file specifies coding convention for future files. Current files
%   do not necessarily follows the requirements, as they were created
%   before conventions.
%
%   There are three levels of conventions:
%   MUST: an obligation for new files, and existing files not respecting
%     this convention should be changed to respect it.
%   SHOULD: it is better to do it, but not necessary.
%   MAY: it is optional, but commonly used by some developpers.
%
%   Name of functions
%
%   Java convention for file names must be used:
%   * Each function name must be written in lower case, using english
%       denomination.
%   * For compound function names, the first letter of each component after
%       the first one is capitalized.
%       Ex: 'distancePointPlane'.
%
%   Other convention should be used to avoid name conflicts:
%   * For functions which work on a complex data type, the name of the data
%       type should be used as a prefix of the function name.
%       Ex: functions working on polygons have a prototype on the form
%       'polygonLength', 'polygonArea'...
%   * Functions working with geometrical primitives can use a suffix for
%       indicating the dimension of application.
%       Ex: 'distancePointLine3d'.
%       This convention was not used for the beginning of the library
%       development, but new functions will try to conform to it.
%   * Conversion functions use 'To' between radicals. Ex:
%       ellipseToPolygon, edgeToLine, circleArcToPolyline...
%
%   Function header
%
%   * First line of documentation must contain the name of the function,
%   followed by a short description. Name and function may be separated by
%   colon. Name must be capitalized, without punctuation before
%   description (to remain valid with matlab documentation tools).
%   Ex:   'SLICER interactive visualization of 3D images'
%
%   * Main part of header describes usage of function/script with various
%   parameters. Function name should be the real name. Variable names
%   should be capitalized. This does not follow matlab's standard, but
%   allows to copy-paste pieces of code.
%   Ex: 'RES = myBasicFunction(PARAM1, PARAM2, OPTION);'
%
%   * The description of each calling syntax should be written using
%   conjugated verbs.
%
%   * End of header should contain name of author, with date of creation.
%   It may also contains author affiliation (company ...), and date of last
%   revision.
%
%
%   Function code
%
%   * Please comment the code ! For algorithmic parts, as much comments as
%   code is a good start ...
%   * Indentation level should be marked by four space characters.
%
%
%   Tips & Tricks
%
%   * The use of '...' at the end of a line can cause compatibility
%   problems from Linux to Windows.
%   * matlab help uses the file 'Contents.m', not 'contents.m' !
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - MIA MathCell
%   created the 23/03/2005
%

%   HISTORY
%   2007-03-05 add convention for datatype prefix for function names, and
%       library organization
%   2011-01-11 update for matGeom project

help codingConventions;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function featureNamesNew=featureNamesNewFromOld(featureNamesOld, ...
                                                scoreFeaturesOld, ...
                                                scoreFeaturesNew)
% Figures out the feature names in the new feature lexicon, given the old
% feature names, the old scoreFeatures, and the new
% scoreFeatures.

% figure out which of the old were retained, which of the new are not
% present in the old
[kept,added]= ...
  setDifferencesScoreFeatures(scoreFeaturesOld,scoreFeaturesNew);

% Get a list of the deleted feature names
scoreFeatureNamesOld={scoreFeaturesOld(:).scorefilename};
scoreFeatureNamesDeleted=scoreFeatureNamesOld(~kept);

% Make a function that computes whether a feature name was kept
featureNameKept=@(featureName)(~ismember(featureName,scoreFeatureNamesDeleted));

% Get the feature names common to both old a new by filtering the old
% list
featureNamesCommon= ...
  cellFilter(featureNameKept,featureNamesOld);

% Get the score feature names that are truly novel
scoreBaseNameListNew={scoreFeaturesNew(:).scorefilename};
featureNamesAdded=scoreBaseNameListNew(added);

% Merge the common and added lists
featureNamesNew = ...
  [featureNamesCommon;featureNamesAdded];

end

author: kristinbranson

repo: JAABA


text:
%PDF2EPS  Convert a pdf file to eps format using pdftops
%
% Examples:
%   pdf2eps source dest
%
% This function converts a pdf file to eps format.
%
% This function requires that you have pdftops, from the Xpdf suite of
% functions, installed on your system. This can be downloaded from:
% http://xpdfreader.com
%
% Inputs:
%   source - filename of the source pdf file to convert. The filename is
%            assumed to already have the extension ".pdf".
%   dest - filename of the destination eps file. The filename is assumed to
%          already have the extension ".eps".

% Copyright (C) Oliver Woodford 2009-2010, Yair Altman 2015-

% Thanks to Aldebaro Klautau for reporting a bug when saving to
% non-existant directories.

% 22/09/2018 - Xpdf website changed to xpdfreader.com

function pdf2eps(source, dest)
    % Construct the options string for pdftops
    options = ['-q -paper match -eps -level2 "' source '" "' dest '"'];

    % Convert to eps using pdftops
    [status, message] = pdftops(options);

    % Check for error
    if status
        % Report error
        if isempty(message)
            error('Unable to generate eps. Check destination directory is writable.');
        else
            error(message);
        end
    end

    % Fix the DSC error created by pdftops
    fid = fopen(dest, 'r+');
    if fid == -1
        % Cannot open the file
        return
    end
    fgetl(fid); % Get the first line
    str = fgetl(fid); % Get the second line
    if strcmp(str(1:min(13, end)), '% Produced by')
        fseek(fid, -numel(str)-1, 'cof');
        fwrite(fid, '%'); % Turn ' ' into '%'
    end
    fclose(fid);
end

author: altmany


text:
function r = normal_dataset ( m, n, seed, mu, a )

%*****************************************************************************80
%
%% NORMAL_DATASET generates a multivariate normal dataset and writes it to a file.
%
%  Usage:
%
%    x = normal_dataset ( m, n, seed, mu, a )
%
%    where
%
%    * M the spatial dimension;
%    * N the number of points to generate;
%    * SEED the seed, a positive integer;
%    * MU is the mean vector of dimension M;
%    * A is the MxM positive definite symmetric variance-covariance matrix;
%    * R is the M by N array created.
%
%    The command creates an M by N multivariate normal dataset and writes
%    it to the file "normal_M_N.txt".
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2009
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'NORMAL_DATASET\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Generate a multivariate normal random dataset.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The program requests input values from the user:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  * M, the spatial dimension,\n' );
  fprintf ( 1, '  * N, the number of points to generate,\n' );
  fprintf ( 1, '  * SEED, a seed for the random number generator.\n' );
  fprintf ( 1, '  * MU, the mean vector of length M.\n' );
  fprintf ( 1, '  * A, the MxM positive definite symmetric variance-covariance matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The program generates the data and writes it to the file\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    normal_M_N.txt\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  where "M" and "N" are the numeric values.\n' );
%
%  Get the spatial dimension.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    m = input ( '  Enter the spatial dimension M: ' );
  else
    m = str2num ( m );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension M = %d\n', m );
%
%  Get the number of points.
%
  if ( nargin < 2 )
    fprintf ( 1, '\n' );
    n = input ( '  Enter the number of points N: ' );
  else
    n = str2num ( n );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of points N = %d\n', n );
%
%  Get the seed.
%
  if ( nargin < 3 )
    fprintf ( 1, '\n' );
    seed = input ( '  Enter the seed: ' );
  else
    seed = str2num ( seed );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The seed = %d\n', seed );
%
%  Get the mean vector.
%
  if ( nargin < 4 )
    fprintf ( 1, '\n' );
    mu(1:m) = input ( '  Enter the mean vector MU' );
  else
    mu(1:m) = str2num ( mu(1:m) );
  end

  r8vec_print ( m, mu, '  The mean vector MU:' )
%
%  Get the variance-covariance matrix A.
%
  if ( nargin < 5 )
    fprintf ( 1, '\n' );
    a(1:m,1:m) = input ( '  Enter the variance-covariance matrix A' );
  else
    a(1:m,1:m) = str2num ( a(1:m,1:m) );
  end

  r8mat_print ( m, m, a, '  The variance-covariance matrix A:' )
%
%  Compute the data.
%
  r = multinormal_sample ( m, n, a, mu, seed );
%
%  Write it to a file.
%
  output_filename = ...
    strcat ( 'normal_', num2str ( m ), '_', num2str ( n ), '.txt' );

  r8mat_write ( output_filename, m, n, r );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The data was written to the file "%s".\n', ...
    output_filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NORMAL_DATASET:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function [ x, seed ] = multinormal_sample ( m, n, a, mu, seed )

%*****************************************************************************80
%
%% MULTINORMAL_SAMPLE samples a multivariate normal distribution.
%
%  Discussion:
%
%    The multivariate normal distribution for the M dimensional vector X
%    has the form:
%
%      pdf(X) = (2*pi*det(A))**(-M/2) * exp(-0.5*(X-MU)'*inverse(A)*(X-MU))
%
%    where MU is the mean vector, and A is a positive definite symmetric
%    matrix called the variance-covariance matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the dimension of the space.
%
%    Input, integer N, the number of points.
%
%    Input, real A(M,M), the variance-covariance
%    matrix.  A must be positive definite symmetric.
%
%    Input, real MU(M), the mean vector.
%
%    Input, integer SEED, the random number seed.
%
%    Output, real X(M), the points.
%
%    Output, integer SEED, the random number seed.
%

%
%  Compute the upper triangular Cholesky factor R of the variance-covariance
%  matrix.
%
  [ r, info ] = r8po_fa ( m, a );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MULTINORMAL_SAMPLE - Fatal error!\n' );
    fprintf ( 1, ...
      '  The variance-covariance matrix is not positive definite symmetric.\n' );
    error ( 'MULTINORMAL_SAMPLE - Fatal error!' );
  end
%
%  Get an MxN vector of samples of the 1D normal distribution with mean 0
%  and variance 1.
%
  [ x(1:m*n), seed ] = r8vec_normal_01 ( m*n, seed );
%
%  Reshape the vector to an MxN matrix.
%
  x = reshape ( x, m, n );
%
%  Compute MU + R' * X.
%
  for j = 1 : n
    x(1:m,j) = mu(1:m)' + r(1:m,1:m)' * x(1:m,j);
  end

  return
end
function r8mat_print ( m, n, a, title )

%*****************************************************************************80
%
%% R8MAT_PRINT prints an R8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in A.
%
%    Input, integer N, the number of columns in A.
%
%    Input, real A(M,N), the matrix.
%
%    Input, string TITLE, a title to be printed.
%
  r8mat_print_some ( m, n, a, 1, 1, m, n, title );

  return
end
function r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% R8MAT_PRINT_SOME prints out a portion of an R8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 September 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, real A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, a title.
%
  incx = 5;

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

  for j2lo = max ( jlo, 1 ): incx : min ( jhi, n )

    j2hi = j2lo + incx - 1;
    j2hi = min ( j2hi, n );
    j2hi = min ( j2hi, jhi );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col: ' );

    for j = j2lo : j2hi
      fprintf ( 1, '%7d       ', j );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row\n' );

    i2lo = max ( ilo, 1 );
    i2hi = min ( ihi, m );

    for i = i2lo : i2hi

      fprintf ( 1, '%7d  ', i );

      for j = j2lo : j2hi
        fprintf ( 1, '%12g  ', a(i,j) );
      end

      fprintf ( 1, '\n' );

    end

  end

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Discussion:
%
%    An R8MAT is an array of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
%  For smaller data files, and less precision, try:
%
%     fprintf ( output_unit, '  %14.6f', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %24.16f', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function [ r, info ] = r8po_fa ( n, a )

%*****************************************************************************80
%
%% R8PO_FA factors a R8PO matrix.
%
%  Discussion:
%
%    The R8PO storage format is appropriate for a symmetric positive definite
%    matrix and its inverse.  (The Cholesky factor of a R8PO matrix is an
%    upper triangular matrix, so it will be in R8GE storage format.)
%
%    Only the diagonal and upper triangle of the square array are used.
%    This same storage scheme is used when the matrix is factored by
%    R8PO_FA, or inverted by R8PO_INVERSE.  For clarity, the lower triangle
%    is set to zero.
%
%    The positive definite symmetric matrix A has a Cholesky factorization
%    of the form:
%
%      A = R' * R
%
%    where R is an upper triangular matrix with positive elements on
%    its diagonal.  This routine overwrites the matrix A with its
%    factor R.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 February 2004
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Bunch, Moler, Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real A(N,N), the matrix in R8PO storage.
%
%    Output, real R(N,N), the Cholesky factor R in R8GE storage.
%
%    Output, integer INFO, error flag.
%    0, normal return.
%    K, error condition.  The principal minor of order K is not
%    positive definite, and the factorization was not completed.
%
  r(1:n,1:n) = a(1:n,1:n);

  for j = 1 : n

    for k = 1 : j - 1
      t = 0.0;
      for i = 1 : k-1
        t = t + r(i,k) * r(i,j);
      end
      r(k,j) = ( r(k,j) - t ) / r(k,k);
    end

    t = 0.0;
    for i = 1 : j - 1
      t = t + r(i,j)^2;
    end

    s = r(j,j) - t;

    if ( s <= 0.0 )
      info = j;
      return;
    end

    r(j,j) = sqrt ( s );

  end

  info = 0;
%
%  Since the Cholesky factor is stored in R8GE format, be sure to
%  zero out the lower triangle.
%
  for i = 1 : n
    for j = 1 : i-1
      r(i,j) = 0.0;
    end
  end

  return
end
function [ x, seed ] = r8vec_normal_01 ( n, seed )

%*****************************************************************************80
%
%% R8VEC_NORMAL_01 returns a unit pseudonormal R8VEC.
%
%  Discussion:
%
%    The standard normal probability distribution function (PDF) has
%    mean 0 and standard deviation 1.
%
%    This routine can generate a vector of values on one call.  It
%    has the feature that it should provide the same results
%    in the same order no matter how we break up the task.
%
%    Before calling this routine, the user may call RANDOM_SEED
%    in order to set the seed of the random number generator.
%
%    The Box-Muller method is used, which is efficient, but
%    generates an even number of values each time.  On any call
%    to this routine, an even number of new values are generated.
%    Depending on the situation, one value may be left over.
%    In that case, it is saved for the next call.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 July 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of values desired.  If N is negative,
%    then the code will flush its internal memory; in particular,
%    if there is a saved value to be used on the next call, it is
%    instead discarded.  This is useful if the user has reset the
%    random number seed, for instance.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X(N), a sample of the standard normal PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
%  Local parameters:
%
%    Local, integer MADE, records the number of values that have
%    been computed.  On input with negative N, this value overwrites
%    the return value of N, so the user can get an accounting of
%    how much work has been done.
%
%    Local, real R(N+1), is used to store some uniform random values.
%    Its dimension is N+1, but really it is only needed to be the
%    smallest even number greater than or equal to N.
%
%    Local, integer SAVED, is 0 or 1 depending on whether there is a
%    single saved value left over from the previous call.
%
%    Local, integer X_LO_INDEX, X_HI_INDEX, records the range of entries of
%    X that we need to compute.  This starts off as 1:N, but is adjusted
%    if we have a saved value that can be immediately stored in X(1),
%    and so on.
%
%    Local, real Y, the value saved from the previous call, if
%    SAVED is 1.
%
  persistent made;
  persistent saved;
  persistent y;
%
%  I'd like to allow the user to reset the internal data.
%  But this won't work properly if we have a saved value Y.
%  I'm making a crock option that allows the user to signal
%  explicitly that any internal memory should be flushed,
%  by passing in a negative value for N.
%
  if ( n < 0 )
    made = 0;
    saved = 0;
    y = 0.0;
    x = [];
    return
  elseif ( n == 0 )
    x = [];
    return
  end
%
%  Record the range of X we need to fill in.
%
  x_lo_index = 1;
  x_hi_index = n;
%
%  Use up the old value, if we have it.
%
  if ( saved == 1 )
    x(1) = y;
    saved = 0;
    x_lo_index = 2;
  end
%
%  Maybe we don't need any more values.
%
  if ( x_hi_index - x_lo_index + 1 == 0 )
%
%  If we need just one new value, do that here to avoid null arrays.
%
  elseif ( x_hi_index - x_lo_index + 1 == 1 )

    [ r, seed ] = r8vec_uniform_01 ( 2, seed );

    x(x_hi_index) = ...
      sqrt ( -2.0 * log ( r(1) ) ) * cos ( 2.0 * pi * r(2) );
    y =      sqrt ( -2.0 * log ( r(1) ) ) * sin ( 2.0 * pi * r(2) );

    saved = 1;

    made = made + 2;
%
%  If we require an even number of values, that's easy.
%
  elseif ( mod ( x_hi_index - x_lo_index + 1, 2 ) == 0 )

    m = floor ( ( x_hi_index - x_lo_index + 1 ) / 2 );

    [ r, seed ] = r8vec_uniform_01 ( 2*m, seed );

    x(x_lo_index:2:x_hi_index-1) = ...
      sqrt ( -2.0 * log ( r(1:2:2*m-1) ) ) ...
      .* cos ( 2.0 * pi * r(2:2:2*m) );

    x(x_lo_index+1:2:x_hi_index) = ...
      sqrt ( -2.0 * log ( r(1:2:2*m-1) ) ) ...
      .* sin ( 2.0 * pi * r(2:2:2*m) );

    made = made + x_hi_index - x_lo_index + 1;
%
%  If we require an odd number of values, we generate an even number,
%  and handle the last pair specially, storing one in X(N), and
%  saving the other for later.
%
  else

    x_hi_index = x_hi_index - 1;

    m = floor ( ( x_hi_index - x_lo_index + 1 ) / 2 ) + 1;

    [ r, seed ] = r8vec_uniform_01 ( 2*m, seed );

    x(x_lo_index:2:x_hi_index-1) = ...
      sqrt ( -2.0 * log ( r(1:2:2*m-3) ) ) ...
      .* cos ( 2.0 * pi * r(2:2:2*m-2) );

    x(x_lo_index+1:2:x_hi_index) = ...
      sqrt ( -2.0 * log ( r(1:2:2*m-3) ) ) ...
      .* sin ( 2.0 * pi * r(2:2:2*m-2) );

    x(n) = sqrt ( -2.0 * log ( r(2*m-1) ) ) ...
      * cos ( 2.0 * pi * r(2*m) );

    y = sqrt ( -2.0 * log ( r(2*m-1) ) ) ...
      * sin ( 2.0 * pi * r(2*m) );

    saved = 1;

    made = made + x_hi_index - x_lo_index + 2;

  end

  return
end
function r8vec_print ( n, a, title )

%*****************************************************************************80
%
%% R8VEC_PRINT prints a real vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the vector.
%
%    Input, real A(N), the vector to be printed.
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '%6d  %12f\n', i, a(i) );
  end

  return
end
function [ r, seed ] = r8vec_uniform_01 ( n, seed )

%*****************************************************************************80
%
%% R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Second Edition,
%    Springer, 1987,
%    ISBN: 0387964673,
%    LC: QA76.9.C65.B73.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, December 1986, pages 362-376.
%
%    Pierre L'Ecuyer,
%    Random Number Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley, 1998,
%    ISBN: 0471134031,
%    LC: T57.62.H37.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, Number 2, 1969, pages 136-143.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real R(N), the vector of pseudorandom values.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  i4_huge = 2147483647;

  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_UNIFORM_01 - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'R8VEC_UNIFORM_01 - Fatal error!' );
  end

  r = zeros ( n, 1 );

  for i = 1 : n

    k = floor ( seed / 127773 );

    seed = 16807 * ( seed - k * 127773 ) - k * 2836;

    if ( seed < 0 )
      seed = seed + i4_huge;
    end

    r(i) = seed * 4.656612875E-10;

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function [ehmm,Gamma,crithist] = ehmmtrain(data,T,ehmm,Gamma,residuals)
%
% Train ehmm using using Variational Framework
%
% INPUTS:
%
% data          observations - a struct with X (time series) and C (classes)
% T             Number of time points for each time series
% ehmm          ehmm structure with options specified in ehmm.train
% Gamma         Initial state courses
% residuals     in case we train on residuals, the value of those.
%
% OUTPUTS
% ehmm           estimated ehmm
% Gamma         estimated p(state | data)
% crithist     historic of Gamma amount of change across iterations
%
% Author: Diego Vidaurre,
%         CFIN, Aarhus University / OHBA, University of Oxford (2021)

setxx_ehmm;
crithist = [];

for cycle = 1:ehmm.train.cyc

%     tt =  799676-10000:799875+10000;
%     figure(2);imagesc(Gamma(tt,:)');colorbar;%xlim([tt(1) tt(end)])
%     W = zeros(length(ehmm.state(1).W.Mu_W),ehmm.train.K);
%     for k = 1:ehmm.train.K+1, W(:,k) = ehmm.state(k).W.Mu_W; end;
%     fit = hmmspectramar(residuals,size(residuals,1),ehmm);
%     plot_hmmspectra (fit,[],[],5);legend
    %figure(5);imagesc(W);colorbar

    %%% E step - state inference
    if ehmm.train.updateGamma
        if cycle > 1 && strcmpi(ehmm.train.stopcriterion,'ChGamma')
            Gamma0 = Gamma;
        end
        if cycle == 1
            [Gamma,~,Xi] = hsinference(data,T,ehmm,residuals,[],XX,Gamma);
        else
            [Gamma,~,Xi] = hsinference(data,T,ehmm,residuals,[],XX,Gamma,Xi);
        end
    else
        Xi = approximateXi_ehmm(Gamma,size(Gamma,1)+ehmm.train.order,ehmm.train.order);
    end

    %%% M STEP

    % Observation model
    if ehmm.train.updateObs
        ehmm = obsupdate_ehmm(Gamma,ehmm,residuals,XX);
    end

    % Transition matrices and initial state
    if ehmm.train.updateP
        ehmm = hsupdate_ehmm(Xi,Gamma,T,ehmm);
    end

    % Stopping conditions and reporting
    if strcmpi(ehmm.train.stopcriterion,'FreeEnergy')
        % computation of free energy is not exact
        crithist(end+1) = sum(evalfreeenergy_ehmm(T,Gamma,Xi,ehmm,residuals,XX));
        if ehmm.train.verbose
            fprintf('cycle %i Approx free energy = %.10g \n',cycle,crithist(end));
        end
        if cycle > 1
            chgFrEn = (crithist(end) - crithist(end-1)) ...
                / abs(crithist(1) - crithist(end));
            if (abs(chgFrEn) < ehmm.train.tol), break; end
        end

    elseif strcmpi(ehmm.train.stopcriterion,'ChGamma')
        if cycle > 1
            crithist(end+1) = 100 * mean(sum(abs(Gamma0 - Gamma),2)/2 );
            if ehmm.train.verbose
                fprintf('cycle %i Gamma change = %.3g %% \n',...
                    cycle,crithist(end));
            end
            if (crithist(end) < ehmm.train.tol), break; end
        else
            crithist(end+1) = 0;
            if ehmm.train.verbose
                fprintf('cycle 1  \n')
            end
        end
    else % log likelihood
        crithist(end+1) = -sum(evalfreeenergy_ehmm(T,Gamma,Xi,ehmm,residuals,...
            XX,[0 1 1 0 0]));
        if ehmm.train.verbose
            fprintf('cycle %i Approx log likelihood = %.10g \n',cycle,crithist(end));
        end
        if cycle > 1
            chL = (crithist(end) - crithist(end-1)) ...
                / abs(crithist(1) - crithist(end));
            if  (abs(chL) < ehmm.train.tol), break; end  % (chL < 0) ||
         end
    end

    % plot state time courses if requested
    if ehmm.train.plotGamma > 0
        figure(100);clf(100);
        if ehmm.train.plotGamma == 1 % continuous data
            plot_Gamma (Gamma,T,1);
        elseif ehmm.train.plotGamma == 2 % full plot
            plot_Gamma (Gamma,T,0);
        end
        drawnow
    end

    if ~ehmm.train.updateGamma
        break % one iteration is enough
    end

end

for k = 1:K
    if isfield(ehmm.state(k),'cache')
        ehmm.state(k) = rmfield(ehmm.state(k),'cache');
    end
end

if ehmm.train.verbose
    str = 'ehmm '; str2 = 'chains';
    if ~isfield(ehmm.train,'distribution') || strcmp(ehmm.train.distribution,'Gaussian')
        fprintf('%s Model: %d %s, %d data samples, order %d \n', ...
            str,K,str2,sum(T),ehmm.train.order);
    elseif strcmp(ehmm.train.distribution,'logistic')
        fprintf('%s Model: %d %s, %d data samples, logistic regression model. \n', ...
            str,K,str2,sum(T));
    end
    if ehmm.train.useMEX==0
        fprintf('MEX file was not used \n')
    else
        fprintf('MEX file was used for acceleration \n')
    end
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function [llh] = tapas_mdlinear_hier_llh(data, theta, ptheta)
%% Likelihood of the nodes of a linear multivariate model with diagonal
% precisions on each parameter.
%
% It is the likelihood of
%
% p(y|x, beta, pi) = N(y - x*b, 1/pi)
%
% Note that this is a multivariate regression problem, so beta is a matrix
% in which each column correspond to one parameter of the model.
%

%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

y = data.y;
u = data.u;

x = data.u.x;
theta = theta.y;

[np, nc] = size(y);

llh = zeros(np, nc);

ln2pi = log(2 * pi);

for j = 1:nc
    % Mean
    beta = theta{j}.mu;
    % Precision
    pe = theta{j}.pe;
    % Regressors

    lpe = log(pe);

    % number of subjects x number of parameters
    py = u.x * beta;
    [ins, inp] = size(py);
    % TODO this is a problem because it is not consistent
    ty = cell2mat(y(:, j));

    % Residuals
    r = (reshape(ty, inp, ins)' - py);

    % Prediction error
    e = sum(bsxfun(@times, pe, r .* r), 2);

    % This is a vector. For each subject needs to include the precision.
    llh(:, j) = sum(-0.5 * ln2pi +  0.5 * lpe) - 0.5 * e;
end

end

author: translationalneuromodeling

repo: tapas


text:
function s = y_to_s_eg_civil ( y )

%*****************************************************************************80
%
%% Y_TO_S_EG_CIVIL writes an Egyptian Civil year into a string.
%
%  Format:
%
%    EN YearNumber
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Discussion:
%
%    "EN" stands for the Era of Nabonassar, a Babylonian king who
%    acceded in 747 BC, used by the astronomer Ptolemy to assign
%    an artificial starting year for the Egyptian calendar.
%
%  Modified:
%
%    28 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year.
%
%    Output, string S, a representation of the year.
%
  s = sprintf ( 'EN %d', y );

  return
end


text:
function [hrf, fit, e, param] = Fit_Spline(tc, TR, Run, T)
% function [hrf, fit, e, param] = Fit_Spline(tc, TR, Run, T)
%
% Fits Spline model
%
% INPUTS:
%
% tc    - time course
% TR    - time resolution
% Runs  - expermental design
% T     - length of estimated HRF
%
% OUTPUTS:
%
% hrf   - estimated hemodynamic response function
% fit   - estimated time course
% e     - residual time course
% param - estimated amplitude, height and width
%
% Created by Martin Lindquist on 03/06/23

numstim = length(Run);  % Number conditions
len = length(Run{1});   % length of run
t=1:TR:T;
tlen = length(t);       % Number of time points in HRF

K = 8;                         % Number of b-spline basis (This cxan be set to specification)
norder = 4;                    % Order of b-spline basis (This cxan be set to specification)

% Create design matrix
basis = create_bspline_basis([0,tlen], K+3, norder);
B = eval_basis((1:tlen),basis);
B = B(:,3:end-1);

Wi = zeros(len, numstim*K);
for j=1:numstim
    Wji = tor_make_deconv_mtx3(Run{j},tlen,1);
    Wi(:,(j-1)*K+1:j*K) = Wji(:,1:tlen)*B;
end

X = [ones(len,1) Wi];

% Fit model
b = pinv(X)*tc;
e = tc-X*b;
fit = X*b;

b2 = reshape(b(2:end),K,numstim);

% Get parameters

hrf = B*b2;

for i=1:numstim
    param(:,i) = get_parameters2(hrf(:,i),1:length(t));
end

end

% END MAIN FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function A = spm_mesh_area(M,PF)
% Compute the surface area of a triangle mesh
% M        - patch structure: vertices and faces must be mx3 and nx3 arrays
%            or 3xm array of edge distances
% PF       - logical indicating whether to return surface area as a whole
%            or per face [default: false]
%
% A        - surface area
%__________________________________________________________________________
%
% Computed using numerically stable version of Heron's formula:
% See https://www.wikipedia.org/wiki/Heron%27s_formula
%__________________________________________________________________________
% Copyright (C) 2010-2018 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_mesh_area.m 7367 2018-07-06 16:04:27Z guillaume $

if isnumeric(M)
    A = M;
else
    A = M.vertices(M.faces',:);
    A = permute(reshape(A',3,3,[]),[2 1 3]);
    A = squeeze(sqrt(sum((A([1 2 3],:,:) - A([2 3 1],:,:)).^2,2)));
end
A = sort(A,1,'descend');
A = ( A(1,:) + ( A(2,:) + A(3,:) ) ) .* ...
    ( A(3,:) - ( A(1,:) - A(2,:) ) ) .* ...
    ( A(3,:) + ( A(1,:) - A(2,:) ) ) .* ...
    ( A(1,:) + ( A(2,:) - A(3,:) ) );
A(A<0) = 0;
A = 1/4 * sqrt(A);

if nargin < 2 || ~PF
    A = sum(A);
end

author: spm


text:
function [x,fval,exitflag,output]=cplexmiqcp(H, f, Aineq, bineq, Aeq, beq, l, Q, r, sostype, sosind, soswt, lb, ub, ctype, x0, options)
%%
% Purpose
% Solve quadratically constrained linear or quadratic integer programming
% problems.
%
% Syntax
%    x = cplexmiqcp(H,f,Aineq,bineq)
%    x = cplexmiqcp(H,f,Aineq,bineq,Aeq,beq)
%    x = cplexmiqcp(H,f,Aineq,bineq,Aeq,beq,l,Q,r)
%    x = cplexmiqcp(H,f,Aineq,bineq,Aeq,beq,l,Q,r,sostype,sosind,soswt)
%    x = cplexmiqcp(H,f,Aineq,bineq,Aeq,beq,l,Q,r,sostype,sosind,soswt,lb,
%        ub)
%    x = cplexmiqcp(H,f,Aineq,bineq,Aeq,beq,l,Q,r,sostype,sosind,soswt,lb,
%        ub,ctype,x0)
%    x = cplexmiqcp(H,f,Aineq,bineq,Aeq,beq,l,Q,r,sostype,sosind,soswt,lb,
%        ub,ctype,x0,options)
%    x = cplexmiqcp(problem)
%    [x,fval] = cplexmiqcp(...)
%    [x,fval,exitflag] = cplexmiqcp(...)
%    [x,fval,exitflag,output] = cplexmiqcp(...)
%    [x,fval,exitflag,output,lambda] = cplexmiqcp(...)
%
% Description
% Finds the minimum of a problem specified by
%    min      0.5*x'*H*x+f*x or f*x
%    st.      Aineq*x      <= bineq
%             Aeq*x         = beq
%             l*x + x'*Q*x <= r
%             lb <= x <= ub
%             x belongs to BICSN
%
% f, bineq, beq, l, r, lb and ub are column vectors.
% H, Aineq, Aeq, and Q are matrices.
% x is a BICSN vector -- that is, its individual entries are each required
% to be binary, general integer, continuous, semi-continuous or
% semi-integer.
%
% x = cplexmiqcp(H,f,Aineq,bineq) solves the mixed integer programming
% problem min 1/2*x'*H*x + f*x subject to Aineq*x <= bineq. If no
% quadratic objective term exists, set H=[].
%
% x = cplexmiqcp(H,f,Aineq,bineq,Aeq,beq) solves the preceding problem
% while additionally satisfying the equality constraints Aeq*x = beq. If no
% inequalities exist, set Aineq=[] and bineq=[].
%
% x = cplexmiqcp(H,f,Aineq,bineq,Aeq,beq,l,Q,r) solves the preceding
% problem while additionally satisfying the quadratic inequality
% constraints l*x + x'*Q*x <= r. If no equalities exist, set Aeq=[] and beq=[].
%
% x = cplexmiqcp(f,Aineq,bineq,Aeq,beq,l,Q,r,sostype,sosind,soswt) solves
% the preceding problem with the additional requirement that the SOS
% constraints are satisfied. If no quadratic inequalities exist, set l=[],
% Q=[] and r=[].
%
% x = cplexmiqcp(H,f,Aineq,bineq,Aeq,beq,l,Q,r,sostype,sosind,soswt,lb,ub)
% defines a set of lower and upper bounds on the design variables, x, so
% that the solution is in the range lb <= x <= ub. If no SOS constraints
% exist, set sostype=[],sosind=[] and soswt=[].
%
% x =
% cplexmiqcp(H,f,Aineq,bineq,Aeq,beq,l,Q,r,sostype,sosind,soswt,lb,ub,ctype
% ) defines the types for each of the design variables. If no bounds exist,
% set lb=[] and ub=[].
%
% x =
% cplexmiqcp(H,f,Aineq,bineq,Aeq,beq,l,Q,r,sostype,sosind,soswt,lb,ub,x0)
% sets the starting point to x0. If all design variables are continuous,
% set ctype=[].
%
% x = cplexmiqcp(H,f,Aineq,bineq,Aeq,beq,l,Q,r,lb,ub,ctype,x0,options) minimizes
% with the optimization options specified in the structure options, which
% can be created using the function cplexoptimset If you do not wish to
% give an initial point, set x0=[].
%
% x = cplexmiqcp(problem) where problem is a structure.
%
% [x,fval] = cplexmiqcp(...) returns the value of the objective function at
% the solution x: fval = 0.5*x'*H*x + f*x.
%
% [x,fval,exitflag] = cplexmiqcp(...) returns a value exitflag that
% describes the exit condition of cplexmiqcp.
%
% [x,fval,exitflag,output] = cplexmiqcp(...) returns a structure output
% that contains information about the optimization.
%
% Input Arguments
% H         Double matrix for objective function
% f         Double column vector for objective function
% Aineq     Double matrix for linear inequality constraints
% bineq     Double column vector for linear inequality constraints
% Aeq       Double matrix for linear equality constraints
% beq       Double column vector for linear equality constraints
% l         Double column vector or matrix
%           Linear part of quadratic constraints
% Q         Double matrix or double matrix cell for quadratic
%           constraints
% r         Double or double row vector
%           Righthand side of quadratic inequality constraints
% sostype   String with possible char values  '1', '2'
% sosind    Double column vector or column vector cell of indices for
%           the SOSs to be added
% soswt     Double column vector or column vector cell of weights for
%           the SOSs to be added
% lb        Double column vector of lower bounds
% ub        Double column vector of upper bounds
% ctype     String with possible char values 'B','I','C','S','N'
%           ctype(j) to 'B', 'I','C', 'S', or 'N' to indicate
%           that x(j) should be binary, general integer,
%           continuous, semi-continuous or semi-integer
%           (respectively).
% x0        Double column vector for initial point for x
% options   Options structure created with cplexoptimset
%
% problem   Structure containing the following fields:
%           H         Double matrix for objective function
%           f         Double column vector for objective function
%           Aineq     Double matrix for linear inequality constraints
%           bineq     Double column vector for linear inequality constraints
%           Aeq       Double matrix for linear equality constraints
%           beq       Double column vector for linear equality constraints
%      	    qc 	      Struct vector
%  	        qc(i).a 	 Double column vector for linear part of the
%                        quadratic constraint
%  	        qc(i).rhs    Double for righthand side for quadratic constraint
%  	        qc(i).Q 	 Double matrix for quadratic part of the
%                        quadratic constraint
%        	sos 	  Struct vector representing the SOSs
%  	        sos(i).type  String with possible char values  '1', '2'
%       	sos(i).ind   Double column vector of indices for the SOSs
%                        to be added
%  	        sos(i).wt    Double column vector of weights for the SOSs
%                        to be added
%           lb        Double column vector of lower bounds
%           ub        Double column vector of upper bounds
%           ctype     String with possible char values 'B','I','C','S','N'
%                     ctype(j) to 'B', 'I','C', 'S', or 'N' to indicate
%                     that x(j) should be binary, general integer,
%                     continuous, semi-continuous or semi-integer
%                     (respectively).
%           x0        Double column vector for initial point for x
%           options   Options structure created with cplexoptimset
%
% Output Arguments
% x         Solution found by the optimization function. If exitflag > 0,
%           then x is a solution; otherwise, x is the value of the
%           optimization routine when it terminated prematurely.
% fval      Value of the objective function at the solution x
% exitflag  Integer identifying the reason the optimization algorithm
%           terminated
% output    Structure containing information about the optimization. The
%           fields of the structure are:
%           iterations         Number of iterations
%           algorithm          Optimization algorithm used
%           message            Exit message
%           time               Execution time of the algorithm
%           cplexstatus        Status code of the solution
%           cplexstatusstring  Status string of the solution
%
%
%  See also cplexoptimset
%

% ---------------------------------------------------------------------------
% File: cplexmiqcp.m
% ---------------------------------------------------------------------------
% Licensed Materials - Property of IBM
% 5725-A06 5725-A29 5724-Y48 5724-Y49 5724-Y54 5724-Y55
% Copyright IBM Corporation 2008, 2010. All Rights Reserved.
%
% US Government Users Restricted Rights - Use, duplication or
% disclosure restricted by GSA ADP Schedule Contract with
% IBM Corp.
% ---------------------------------------------------------------------------

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function h = hypot3(dx, dy, dz)
%HYPOT3 Diagonal length of a cuboidal 3D box
%
%   h = hypot3(a, b, c)
%   computes the quantity sqrt(a^2 + b^2 + c^2), by avoiding roundoff
%   errors.
%
%   Example
%     % Compute diagonal of unit cube
%     hypot3(1, 1, 1)
%     ans =
%          1.7321
%
%     % Compute more complicated diagonal
%     hypot3(3, 4, 5)
%     ans =
%         7.0711
%
%   See also
%   hypot, vectorNorm3d
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2012-04-29,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2012 INRA - Cepia Software Platform.

h = hypot(hypot(dx, dy), dz);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function poly_product_test ( )

%*****************************************************************************80
%
%% POLY_PRODUCT_TEST tests POLY_PRODUCT.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/triangle_integrals/poly_product_test.m
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLY_PRODUCT_TEST:\n' );
  fprintf ( 1, '  POLY_PRODUCT computes the product of two X,Y polynomials.\n' );
%
%  P1 = ( 1 + 2 x + 3 y )
%  P2 = ( 4 + 5 x )
%  P3 = 4 + 13x + 12y + 10x^2 + 15xy + 0y^2
%
  d1 = 1;
  p1 = [ 1.0, 2.0, 3.0 ];
  fprintf ( 1, '\n' );
  poly_print ( d1, p1, '  p1(x,y)' );

  d2 = 1;
  p2 = [ 4.0, 5.0, 0.0 ];
  fprintf ( 1, '\n' );
  poly_print ( d2, p2, '  p2(x,y)' );

  [ d3, p3 ] = poly_product ( d1, p1, d2, p2 );
  fprintf ( 1, '\n' );
  poly_print ( d3, p3, '  p3(x,y) = p1(x,y) * p2(x,y)' );

  dc = 2;
  pc = [ 4.0, 13.0, 12.0, 10.0, 15.0, 0.0 ];
  fprintf ( 1, '\n' );
  poly_print ( dc, pc, '  Correct answer: p3(x,y)=p1(x,y)*p2(x,y)' );
%
%  P1 = ( 1 - 2 x + 3 y - 4x^2 + 5xy - 6y^2)
%  P2 = ( 7 + 3x^2 )
%  P3 = 7
%    - 14x + 21y
%    - 25 x^2 + 35 xy -  42y^2
%    -  6x^3  + 9x^2y +  0xy^2   + 0y^3
%    - 12x^4 + 15x^3y - 18x^2y^2 + 0 xy^3 + 0y^4
%
  d1 = 2;
  p1 = [ 1.0, -2.0, 3.0 -4.0, +5.0, -6.0 ];
  fprintf ( 1, '\n' );
  poly_print ( d1, p1, '  p1(x,y)' );

  d2 = 2;
  p2 = [ 7.0, 0.0, 0.0, 3.0, 0.0, 0.0 ];
  fprintf ( 1, '\n' );
  poly_print ( d2, p2, '  p2(x,y)' );

  [ d3, p3 ] = poly_product ( d1, p1, d2, p2 );
  fprintf ( 1, '\n' );
  poly_print ( d3, p3, '  p3(x,y) = p1(x,y) * p2(x,y)' );

  dc = 4;
  pc = [ ...
      7.0, ...
    -14.0,  21.0, ...
    -25.0, +35.0, -42.0, ...
     -6.0,   9.0,   0.0, 0.0, ...
    -12.0, +15.0, -18.0, 0.0, 0.0 ];
  fprintf ( 1, '\n' );
  poly_print ( dc, pc, '  Correct answer: p3(x,y)=p1(x,y)*p2(x,y)' );

  return
end


text:
function d = dcircle ( p, xc, yc, r )

%*****************************************************************************80
%
%% DCIRCLE returns the signed distance of one or more points to a circle.
%
%  Discussion:
%
%    The corresponding routine in 3D is DSPHERE.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, real P(NP,2), the point coordinates.
%
%    Input, real XC, YC, the coordinates of the center of the circle.
%
%    Input, real R, the radius of the circle.
%
%    Output, real D(NP), the signed distance of each point to the
%    circle.  The point is inside, on, or outside the circle depending
%    on whether D is negative, zero, or positive.
%
  d = sqrt ( ( p(:,1) - xc ).^2 + ( p(:,2) - yc ).^2 ) - r;

  return
end


text:
%% patchPathAngles
% Below is a demonstration of the features of the |patchPathAngles| function

%%
clear; close all; clc;

%% Syntax
% |[indAngles]=patchPathAngles(F,V,ind,isClosedLoop);|

%% Description
% The |patchPathAngles| function computes the angles between adjacent edges
% on a curve on a patch.

%% Examples

%% Example 1: Get angles on a closed path defined on the the boundary of a patch
% Create test data set
w=1;
[X,Y]=ndgrid(linspace(0,w,15));
Z=ones(size(X));
C=tril(Z);
[F,V,C]=surf2patch(X,Y,Z,C);
C=vertexToFaceMeasure(F,C)>0;

logicKeep=C>0;
F=F(logicKeep,:);
C=C(logicKeep,:);
[F,V]=patchCleanUnused(F,V);

%%
% Get boundary curve

Eb=patchBoundary(F);
indBoundaryCurve=edgeListToCurve(Eb);
indBoundaryCurve=indBoundaryCurve(1:end-1)'; %Start=End for closed curve so remove double entry

%%
% Calculate mesh path angles

isClosedPath=1;
[A]=patchPathAngles(F,V,indBoundaryCurve,isClosedPath);

%%

A=180*(A./pi); % Conver to degrees

% Display unique angles in set
unique(A)

%%

cFigure; hold on;
gpatch(F,V,'kw');
plotV(V(indBoundaryCurve,:),'k-','LineWidth',3);
scatterV(V(indBoundaryCurve,:),75,A,'filled');
colormap gjet; colorbar; caxis([0 360]);
axisGeom; view(2);
drawnow;

%% Example 2: Study angles for altered shape

V(:,1)=V(:,1)-V(:,2);

%%
% Calculate mesh path angles

isClosedPath=1;
[A]=patchPathAngles(F,V,indBoundaryCurve,isClosedPath);

%%

A=180*(A./pi); % Conver to degrees

% Display unique angles in set
unique(A)

%%

cFigure; hold on;
gpatch(F,V,'kw');
plotV(V(indBoundaryCurve,:),'k-','LineWidth',3);
scatterV(V(indBoundaryCurve,:),75,A,'filled');
colormap gjet; colorbar; caxis([0 360]);
axisGeom; view(2);
drawnow;

%% Example 3: Get angles on a non-closed path defined on the the boundary of a patch

%%
% Create path segment
indBoundaryCurve=indBoundaryCurve(1:6);

%%
% Calculate mesh path angles
isClosedPath=0;
[A]=patchPathAngles(F,V,indBoundaryCurve,isClosedPath);

%%

A=180*(A./pi); % Conver to degrees

% Display unique angles in set
unique(A(~isnan(A)))

%%

cFigure; hold on;
gpatch(F,V,'kw');
plotV(V(indBoundaryCurve,:),'k-','LineWidth',3);
scatterV(V(indBoundaryCurve,:),75,A,'filled');
colormap gjet; colorbar; caxis([0 360]);
axisGeom; view(2);
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function trnslate(H,X)
    %
    %    trnslate(Handle,offset)
    %    translate Handle by amount offset=[x_offset, y_offset, z_offest]
    %
    %    See also ROTATE, SCALHAND
    %
    %    Richard G. Cobb    3/96
    %
    dad=get(H,'Parent');
    xlim=get(dad,'Xlim');
    xm=max(get(H,'Xdata'));
    xn=min(get(H,'Xdata'));
    if xm + X(1) > max(xlim)
        X(1) = max(xlim) - xm ;
    elseif xn + X(1) < min(xlim)
        X(1) = min(xlim) - xn ;
    end
    xlim=get(dad,'Ylim');
    xm=max(get(H,'Ydata'));
    xn=min(get(H,'Ydata'));
    if xm + X(2) > max(xlim)
        X(2) = max(xlim) - xm ;
    elseif xn + X(2) < min(xlim)
        X(2) = min(xlim) - xn ;
    end
    xlim=get(dad,'Zlim');
    xm=max(get(H,'Zdata'));
    xn=min(get(H,'Zdata'));
    if xm + X(3) > max(xlim)
        X(3) = max(xlim) - xm ;
    elseif xn + X(3) < min(xlim)
        X(3) = min(xlim) - xn ;
    end
    nx = (get(H,'Xdata')+X(1));
    ny = (get(H,'Ydata')+X(2));
    nz = (get(H,'Zdata')+X(3));
    set(H,'Xdata',nx)
    set(H,'Ydata',ny)
    set(H,'Zdata',nz)
    ud=get(H,'Userdata');
    ud(1)=nx(6);
    ud(2)=ny(6);
    ud(3)=nz(6);
    ud(4)=nx(1);
    ud(5)=ny(1);
    ud(6)=nz(1);

    set(H,'Userdata',ud);

    %eof

author: CelsoReyes


text:
function outputWaveforms = load_antelope(obj, specificDatabase)
   % load a waveform from antelope
   %  w = obj.load_antelope(specificDatabase)
   %   database is the antelope database

   % AUTHOR: Celso Reyes
   % MODIFICATIONS: Glenn Thompson, Carl Tape

   %TODO: maybe order of operations can be changed to avoid unpacking
   %request when unnecessary.

   %call this routine for each database, then return the waveforms.
   if exist('specificDatabase','var')
      obj.dbname = specificDatabase;
      outputWaveforms = emptyWaveform();

      %convert each channeltag into
      obj = obj.buildSearchCriteria();
      for i = 1:numel(obj.searchcriteria)
         %if multiple traces will result, then there may be multiple records for tr
         [tr, dbname, fdb] = get_antelope_traces(obj, obj.searchcriteria{i});
         w = cycleThroughTraces(tr, obj.combinetraces);
         outputWaveforms = [outputWaveforms; w]; %#ok<AGROW>
      end
      dbclose(fdb);
   else
      outputWaveforms = RecursivelyLoadFromEachDatabase(obj);
   end
end

author: geoscience-community-codes

repo: GISMO


text:
function disp( f )
%DISP Display a BALLFUN object to the command line.
%
% See also DISPLAY.

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

Nfunctions = numel(f);
% Vectorize the array
f = f(:);

% Display the object
if ( Nfunctions == 1 )
    % Get information that we want to display:
    display_ballfun(f);
else
    disp('   ballfun array containing')
    fprintf('\n');
    for i = 1:Nfunctions
        display_ballfun(f(i));
    end
end

end

% Display one ballfun function
function display_ballfun(f)

loose = strcmp(get(0, 'FormatSpacing'), 'loose');

if isempty(f)
    fprintf('    empty ballfun\n')
else
S = size(f);
vscl = vscale(f);
% Display the information:
disp('   ballfun object:')
fprintf('      domain           r    lambda    theta    vertical scale\n');
fprintf('     unit ball    %6i  %6i   %6i           %3.2g\n', S(1), S(2), S(3),vscl);
end

if ( loose )
    fprintf('\n');
end

end


text:
function edges = knnGraph(nodes, varargin)
%KNNGRAPH Create the k-nearest neighbors graph of a set of points
%
%   EDGES = knnGraph(NODES)
%
%   Example
%   nodes = rand(10, 2);
%   edges = knnGraph(nodes);
%   drawGraph(nodes, edges);
%
%   See also
%     graphs

% ------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% Created: 2008-07-28,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2008 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

% get number of neighbors for each node
k = 2;
if ~isempty(varargin)
    k = varargin{1};
end

% init size of arrays
n       = size(nodes, 1);
edges   = zeros(k*n, 2);

% iterate on nodes
for i = 1:n
    dists = distancePoints(nodes(i,:), nodes);
    [dists, inds]    = sort(dists); %#ok<ASGLU>
    for j = 1:k
        edges(k*(i-1)+j, :) = [i inds(j+1)];
    end
end

% remove double edges
edges = unique(sort(edges, 2), 'rows');