text:
function figure_num = trim_test ( figure_num )

%*****************************************************************************80
%
%% TRIM_TEST tries to make a plot without all that wasted margin.
%
%  Discussion:
%
%    Surprisingly, a typical MATLAB plot can involve a lot of empty
%    space.  In the simple case where you are plotting a circle, you
%    will get white space in the grid region itself if you enforce
%    equal axis units, since MATLAB wants to display the grid within
%    a rectangle.  But you will also get a substantial margin of white
%    space around the plot area itself.  This can be annoying when you
%    are trying to prepare graphics for use in a publication, or
%    if you want to put several images side by side - the white space,
%    which contains no information, eats up lots of space.
%
%    Can you cut back that white space?  In particular, if you want to
%    plot a circle (or anything else whose shape doesn't correspond to
%    MATLAB's "golden rectangle"), can you make an image that is
%    mathematically correct, and economic in terms of wasted white space?
%
%    You may try various combinations of "axis equal", "axis square",
%    "axis tight", "axis ( [0,1,0,1] )" and so on, all the time coming
%    halfway to your goal, without ever actually achieving an image
%    that is square (or whatever aspect ratio you want), doesn't have
%    wasted graph space, and doesn't have that excessive margin around
%    the plot.
%
%    After wasting a lot of time, I stumbled across a few commands that,
%    used together, seem to achieve the correct result.  However, in my
%    case, the image still looks rectangular in the MATLAB interactive
%    viewer - but comes out nice and square in the PNG image rendered
%    by the PRINT command.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 May 2013
%
%  Author:
%
%    John Burkardt
%
  if ( nargin < 1 )
    figure_num = 0;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIM_TEST\n' );
  fprintf ( 1, '  Try to suppress the wasted margin.\n' );
%
%  Set up the circle data once.
%
  r = 1.0;
  theta = linspace ( 0.0, 2.0 * pi, 121 );
  x = r * cos ( theta );
  y = r * sin ( theta );
%
%  Plot the circle with the default axis.
%
  figure_num = figure_num + 1;
  gcf = figure ( figure_num );
  plot ( x, y, 'b-', 'LineWidth', 3 );
  grid on
  xlabel ( '<--- X --->' );
  ylabel ( '<--- Y --->' );
  title ( 'A circle drawn with the default aspect ratio and margin' )

  filename = 'trim_test01.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saved graphics file "%s".\n', filename );
%
%  Plot the circle again, using commands I DON'T REALLY UNDERSTAND.
%  But they do seem to result in a PNG image that has lost the margin junk.
%
  figure_num = figure_num + 1;
  gcf = figure ( figure_num );

  set ( gcf, 'PaperUnits', 'inches','PaperPosition', [0, 0, 10, 10] );
  pbaspect ( [ 1.0, 1.0, 1.0 ] );
  set ( gcf, 'Units', 'normal' );
  set ( gca, 'Position', [ 0.05, 0.05, 0.90, 0.90 ] )

  plot ( x, y, 'r-', 'LineWidth', 3 );
  grid on
  xlabel ( '<--- X --->' );
  ylabel ( '<--- Y --->' );
  title ( 'A circle drawn with white space trimmed' )

  filename = 'trim_test02.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Saved graphics file "%s".\n', filename );

  return
end


text:
function score = Task1_Min_value(Population,~)
% <min> <single> <real/integer/label/binary/permutation> <large/none> <constrained/none> <expensive/none> <sparse/none> <multitask>
% The minimum objective value of the first task (for multitask optimization)

%------------------------------- 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".
%--------------------------------------------------------------------------

    PopDec = Population.decs;
    score  = Population(PopDec(:,end)==1).best.objs;
    if isempty(score)
        score = nan;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [ctf] = ctf_channel_sets(ctf);

% ctf_channel_sets - Define CTF MEG sensor regions
%
% ctf = ctf_channel_sets(ctf)
%
% This function parses the resource information in ctf.sensor,
% which is returned from ctf_read_res4
%
% INPUTS
%
%   ctf - the struct created by ctf_read_res4
%
% OUTPUTS
%
%   ctf.sensor.index - various regional sensor indices
%
%      <>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%      <                                                      >
%      <                    DISCLAIMER:                       >
%      <                                                      >
%      < THIS PROGRAM IS INTENDED FOR RESEARCH PURPOSES ONLY. >
%      < THIS PROGRAM IS IN NO WAY INTENDED FOR CLINICAL OR   >
%      <                    OFFICIAL USE.                     >
%      <                                                      >
%      <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>
%

% $Revision: 1.1 $ $Date: 2009-01-30 03:49:26 $

% Copyright (C) 2004  Darren L. Weber
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

% Modified: 02/2004, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%ver = '$Revision: 1.1 $';
%fprintf('\nCTF_CHANNEL_SETS [v %s]\n',ver(11:15)); tic;

if ~exist('ctf','var'),
  ctf = ctf_folder;
end

if ~isfield(ctf,'setup'),
  ctf = ctf_read_res4(ctf.folder);
end

%-------------------------------------------------------------
% Find channel types and define channel sets, see the
% System Administrators .pdf, 'Channel Sets Configuration'

%types = unique([ctf.sensor.info.index])
% 0     1     5     9    11    17
%for typeValue = types,
%  index  = find([ctf.sensor.info.index] == typeValue);
%  labels = ctf.sensor.info(index).label;
%end

% ctf.SensorNames = {...
%     'Ref Magnetometer',... % Sensor Type Index of 0
%     'Ref Gradiometer' ,... % Index of 1
%     ''  ,...               % 2
%     '' ,...                % 3
%     '' ,...                % 4
%     'MEG Sensor',...       % 5
%     '' ,...                % 6
%     '',...                 % 7
%     '',...                 % 8
%     'EEG Sensor',...       % 9
%     'ADC Input Current',...% 10 ADC Input Current (Amps)
%     'Stimulation input',...% 11
%     'Video Time',...       % 12
%     '',...                 % 13
%     '',...                 % 14
%     'SAM Sensor',...       % 15
%     'Virtual Channel',...  % 16
%     'System Clock',...     % 17 System Time Ref
%     'ADC Input Voltage',...% 18 ADC Input Voltage (Volts)
%   };

% find the indices of the MEG reference sensors

ctf.sensor.type.meg_ref = [0 1];
ctf.sensor.type.refMagnetometers = 0; % include B*
ctf.sensor.type.refGradiometers  = 1; % include G*,P*,Q*,R*

magref  = find([ctf.sensor.info.index] == 0);
gradref = find([ctf.sensor.info.index] == 1);

ctf.sensor.index.meg_ref_mag  = magref;
ctf.sensor.index.meg_ref_grad = gradref;
ctf.sensor.index.meg_ref  = [magref,gradref];

% find the indices of the MEG head sensors

ctf.sensor.type.meg_sens = 5; % include B*
ctf.sensor.index.meg_sens = find([ctf.sensor.info.index] == 5);

% find the indices of the EEG head sensors

ctf.sensor.type.eeg_sens = 9; % include EEG*
ctf.sensor.index.eeg_sens = find([ctf.sensor.info.index] == 9);

% find the indices of the ADC channels

ctf.sensor.type.adc = 10;
ctf.sensor.index.adc = find([ctf.sensor.info.index] == 10);

% find the indices of the STIM channels

ctf.sensor.type.stim_ref = 11;
ctf.sensor.index.stim_ref = find([ctf.sensor.info.index] == 11);

% find the indices of the STIM channels

ctf.sensor.type.video_time = 12;
ctf.sensor.index.video_time = find([ctf.sensor.info.index] == 12);

% find the indices of the SAM channels (Synthetic Aperture Magnetometry)

ctf.sensor.type.sam = 15;
ctf.sensor.index.sam = find([ctf.sensor.info.index] == 15);

% find the indices of the virtual channels

ctf.sensor.type.virtual_channels = 16;
ctf.sensor.index.virtual_channels = find([ctf.sensor.info.index] == 16);

% find the indices of the system clock

ctf.sensor.type.sclk_ref = 17; % include SCLK*
ctf.sensor.index.sclk_ref = find([ctf.sensor.info.index] == 17);

% combine eeg and meg channels

ctf.sensor.index.all_sens = [ ctf.sensor.index.meg_sens, ctf.sensor.index.eeg_sens ];

all = [ ctf.sensor.index.all_sens, ctf.sensor.index.meg_ref, ...
        ctf.sensor.index.stim_ref, ctf.sensor.index.sclk_ref ];

ctf.sensor.index.other = setdiff([1:ctf.setup.number_channels],all);

% define regions of MEG sensors
labels = {ctf.sensor.info.label};
ctf.sensor.index.meg_left_central    = strmatch('MLC',labels)';
ctf.sensor.index.meg_left_frontal    = strmatch('MLF',labels)';
ctf.sensor.index.meg_left_parietal   = strmatch('MLP',labels)';
ctf.sensor.index.meg_left_occipital  = strmatch('MLO',labels)';
ctf.sensor.index.meg_left_parietal   = strmatch('MLP',labels)';
ctf.sensor.index.meg_left_temporal   = strmatch('MLT',labels)';
ctf.sensor.index.meg_right_central   = strmatch('MRC',labels)';
ctf.sensor.index.meg_right_frontal   = strmatch('MRF',labels)';
ctf.sensor.index.meg_right_occipital = strmatch('MRO',labels)';
ctf.sensor.index.meg_right_parietal  = strmatch('MRP',labels)';
ctf.sensor.index.meg_right_temporal  = strmatch('MRT',labels)';
ctf.sensor.index.meg_mid_central     = strmatch('MZC',labels)';
ctf.sensor.index.meg_mid_frontal     = strmatch('MZF',labels)';
ctf.sensor.index.meg_mid_occipital   = strmatch('MZO',labels)';
ctf.sensor.index.meg_mid_parietal    = strmatch('MZP',labels)';

ctf.sensor.index.meg_left  = [ ...
    ctf.sensor.index.meg_left_central, ...
    ctf.sensor.index.meg_left_frontal, ...
    ctf.sensor.index.meg_left_occipital, ...
    ctf.sensor.index.meg_left_parietal, ...
    ctf.sensor.index.meg_left_temporal ];

ctf.sensor.index.meg_right  = [ ...
    ctf.sensor.index.meg_right_central, ...
    ctf.sensor.index.meg_right_frontal, ...
    ctf.sensor.index.meg_right_occipital, ...
    ctf.sensor.index.meg_right_parietal, ...
    ctf.sensor.index.meg_right_temporal ];

ctf.sensor.index.meg_central = [ ...
    ctf.sensor.index.meg_left_central, ...
    ctf.sensor.index.meg_right_central, ...
    ctf.sensor.index.meg_mid_central ];

ctf.sensor.index.meg_frontal = [ ...
    ctf.sensor.index.meg_left_frontal, ...
    ctf.sensor.index.meg_right_frontal, ...
    ctf.sensor.index.meg_mid_frontal ];

ctf.sensor.index.meg_occipital = [ ...
    ctf.sensor.index.meg_left_occipital, ...
    ctf.sensor.index.meg_right_occipital, ...
    ctf.sensor.index.meg_mid_occipital ];

ctf.sensor.index.meg_parietal = [ ...
    ctf.sensor.index.meg_left_parietal, ...
    ctf.sensor.index.meg_right_parietal, ...
    ctf.sensor.index.meg_mid_parietal ];

ctf.sensor.index.meg_temporal = [ ...
    ctf.sensor.index.meg_left_temporal, ...
    ctf.sensor.index.meg_right_temporal ];

%t = toc; fprintf('...done (%6.2f sec)\n\n',t);

return

author: PatternRecognition

repo: OpenBMI


text:
%DEMGPARD Demonstrate ARD using a Gaussian Process.
%
%	Description
%	The data consists of three input variables X1, X2 and X3, and one
%	target variable  T. The  target data is generated by computing
%	SIN(2*PI*X1) and adding Gaussian  noise, x2 is a copy of x1 with a
%	higher level of added noise, and x3 is sampled randomly from a
%	Gaussian distribution. A Gaussian Process, is trained by optimising
%	the hyperparameters  using the scaled conjugate gradient algorithm.
%	The final values of the hyperparameters show that the model
%	successfully identifies the importance of each input.
%
%	See also
%	DEMGP, GP, GPERR, GPFWD, GPGRAD, GPINIT, SCG
%

%	Copyright (c) Ian T Nabney (1996-2001)

clc;
randn('state', 1729);
rand('state', 1729);
disp('This demonstration illustrates the technique of automatic relevance')
disp('determination (ARD) using a Gaussian Process.')
disp(' ');
disp('First, we set up a synthetic data set involving three input variables:')
disp('x1 is sampled uniformly from the range (0,1) and has a low level of')
disp('added Gaussian noise, x2 is a copy of x1 with a higher level of added')
disp('noise, and x3 is sampled randomly from a Gaussian distribution. The')
disp('single target variable is given by t = sin(2*pi*x1) with additive')
disp('Gaussian noise. Thus x1 is very relevant for determining the target')
disp('value, x2 is of some relevance, while x3 should in principle be')
disp('irrelevant.')
disp(' ');
disp('Press any key to see a plot of t against x1.')
pause;

ndata = 100;
x1 = rand(ndata, 1);
x2 = x1 + 0.05*randn(ndata, 1);
x3 = 0.5 + 0.5*randn(ndata, 1);
x = [x1, x2, x3];
t = sin(2*pi*x1) + 0.1*randn(ndata, 1);

% Plot the data and the original function.
h = figure;
plotvals = linspace(0, 1, 200)';
plot(x1, t, 'ob')
hold on
xlabel('Input x1')
ylabel('Target')
axis([0 1 -1.5 1.5])
[fx, fy] = fplot('sin(2*pi*x)', [0 1]);
plot(fx, fy, '-g', 'LineWidth', 2);
legend('data', 'function');

disp(' ');
disp('Press any key to continue')
pause; clc;

disp('The Gaussian Process has a separate hyperparameter for each input.')
disp('The hyperparameters are trained by error minimisation using the scaled.')
disp('conjugate gradient optimiser.')
disp(' ');
disp('Press any key to create and train the model.')
disp(' ');
pause;

net = gp(3, 'sqexp');
% Initialise the parameters.
prior.pr_mean = 0;
prior.pr_var = 0.1;
net = gpinit(net, x, t, prior);

% Now train to find the hyperparameters.
options = foptions;
options(1) = 1;
options(14) = 30;

[net, options] = netopt(net, options, x, t, 'scg');

rel = exp(net.inweights);

fprintf(1, ...
  '\nFinal hyperparameters:\n\n  bias:\t\t%10.6f\n  noise:\t%10.6f\n', ...
  exp(net.bias), exp(net.noise));
fprintf(1, '  Vertical scale: %8.6f\n', exp(net.fpar(1)));
fprintf(1, '  Input 1:\t%10.6f\n  Input 2:\t%10.6f\n', ...
  rel(1), rel(2));
fprintf(1, '  Input 3:\t%10.6f\n\n', rel(3));
disp(' ');
disp('We see that the inverse lengthscale associated with')
disp('input x1 is large, that of x2 has an intermediate value and the variance')
disp('of weights associated with x3 is small.')
disp(' ');
disp('This implies that the Gaussian Process is giving greatest emphasis')
disp('to x1 and least emphasis to x3, with intermediate emphasis on')
disp('x2 in the covariance function.')
disp(' ')
disp('Since the target t is statistically independent of x3 we might')
disp('expect the weights associated with this input would go to')
disp('zero. However, for any finite data set there may be some chance')
disp('correlation between x3 and t, and so the corresponding hyperparameter remains')
disp('finite.')
disp('Press any key to continue.')
pause

disp('Finally, we plot the output of the Gaussian Process along the line')
disp('x1 = x2 = x3, together with the true underlying function.')
xt = linspace(0, 1, 50);
xtest = [xt', xt', xt'];

cn = gpcovar(net, x);
cninv = inv(cn);
[ytest, sigsq] = gpfwd(net, xtest, cninv);
sig = sqrt(sigsq);

figure(h); hold on;
plot(xt, ytest, '-k');
plot(xt, ytest+(2*sig), '-b', xt, ytest-(2*sig), '-b');
axis([0 1 -1.5 1.5]);
fplot('sin(2*pi*x)', [0 1], '--m');

disp(' ');
disp('Press any key to end.')
pause; clc; close(h); clear all

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%compute length of tailinflvector

function [data,units]=compute_tailinflmag(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
tailinflmag=cell(1,numlarvae);

for i=1:numlarvae
    larva=larvae(i);
    tailinflmag{1,i}=bsxfun(@hypot,trx(larva).xinflection_mm-trx(larva).xtail_mm,trx(larva).yinflection_mm-trx(larva).ytail_mm);
end

units=parseunits('mm');
data=tailinflmag;

author: kristinbranson

repo: JAABA


text:
function dx = susp_sys(x,u)

% Dynamics of the syspension system
x1 = x(1);
x2 = x(2);
x3 = x(3);
x4 = x(4);

% Coefficients
mb = 300;    % kg
mw = 60;     % kg
bs = 1000;   % N/m/s
ks = 16000 ; % N/m
kt = 190000; % N/m
kn = ks/10;  % N/m

% System Dynamics
dx1 = x2;
dx2 = -(ks*(x1-x3)+kn*(x1-x3)^3)/mb - (bs*(x2-x4)-10000*u)/mb;
dx3 = x4;
dx4 =  (ks*(x1-x3)+kn*(x1-x3)^3)/mw + (bs*(x2-x4)-kt*x3-10000*u)/mw;

% Combine the output
dx = [dx1;
	  dx2;
	  dx3;
	  dx4];
end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function sequences = config_sequence(type)
% config sequences for evaluation
% the configuration files are placed in ./sequence_evaluation_config/;
switch type
    case 'test_set'
        dataset_name = 'testing_set.txt';
    case 'all'
        dataset_name = 'all_dataset.txt';
    otherwise
        error('Error in evaluation dataset type! Either ''testing_set'' or ''all''.')
end

% check if the file exists
if ~exist(dataset_name, 'file')
    error('%s is not found!', dataset_name);
end

% load evaluation sequences
fid = fopen(dataset_name, 'r');
i = 0;
sequences = cell(100000, 1);
while ~feof(fid)
    i = i + 1;
    sequences{i, 1} = fgetl(fid);
end
sequences(i+1:end) = [];
fclose(fid);
end

author: HengLan


text:
function disp(t,name)
%DISP Command window display of a matricized tensor (tenmat).
%
%   DISP(T) displays a tensor as matrix with no name.
%
%   DISP(T,NAME) display a tensor as matrix with the given name.
%
%   See also TENMAT, TENMAT/DISPLAY.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if ~exist('name','var')
    name = 'ans';
end

fprintf('%s is a matrix corresponding to a tensor of size %s\n',...
        name,tt_size2str(t.tsize));
fprintf('\t%s.rindices = %s (modes of tensor corresponding to rows)\n',...
        name,['[ ' num2str(t.rindices) ' ]']);
fprintf('\t%s.cindices = %s (modes of tensor corresponding to columns)\n',...
        name,['[ ' num2str(t.cindices) ' ]']);

if isempty(t.data)
    fprintf('\t%s.data = []\n',name);
else
    fprintf('\t%s.data = \n',name);
    output = tt_matrix2cellstr(t.data);
    fprintf('\t\t%s\n',output{:});
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
classdef PTKAirwayGrowing < PTKPlugin
    % PTKAirwayGrowing. Plugin for generating an artificial airway tree
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Volume filling<br> airways'
        ToolTip = 'Grow the airways into the lobes'
        Category = 'Airways'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = true
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = false
        Visibility = 'Developer'
    end

    methods (Static)
        function results = RunPlugin(dataset, reporting)

            left_and_right_lungs = dataset.GetResult('PTKLeftAndRightLungs');
            lobes = dataset.GetResult('PTKLobes');

            airways_by_lobe = dataset.GetResult('PTKReallocateAirwaysLabelledByLobe');

            upper_right_start_segment = airways_by_lobe.StartBranches.RightUpper;
            middle_right_start_segment = airways_by_lobe.StartBranches.RightMid;
            lower_right_start_segment = airways_by_lobe.StartBranches.RightLower;
            upper_left_start_segment = airways_by_lobe.StartBranches.LeftUpper;
            lower_left_start_segment = airways_by_lobe.StartBranches.LeftLower;

            right_lung = left_and_right_lungs.Copy;
            right_lung.ChangeRawImage(left_and_right_lungs.RawImage == 1);

            left_lung = left_and_right_lungs.Copy;
            left_lung.ChangeRawImage(left_and_right_lungs.RawImage == 2);

            template = left_and_right_lungs.BlankCopy;
            approx_number_points = 31000;
            airway_generator = PTKAirwayGenerator(left_and_right_lungs, airways_by_lobe.StartBranches.Trachea, approx_number_points, reporting);

            reporting.ShowProgress('RIGHT - upper');
            reporting.UpdateProgressStage(0, 5);
            lobes_fill = lobes.Copy;
            lobes_fill.ChangeRawImage(lobes.RawImage == PTKColormapLabels.RightUpperLobe);
            lobes_fill.CropToFit;
            lobes_fill.AddBorder(2);
            airway_generator.GrowTree(lobes_fill, upper_right_start_segment, reporting)

            reporting.ShowProgress('RIGHT - middle');
            reporting.UpdateProgressStage(1, 5);
            lobes_fill = lobes.Copy;
            lobes_fill.ChangeRawImage(lobes.RawImage == PTKColormapLabels.RightMiddleLobe);
            lobes_fill.CropToFit;
            lobes_fill.AddBorder(2);
            airway_generator.GrowTree(lobes_fill, middle_right_start_segment, reporting)

            reporting.ShowProgress('RIGHT - lower');
            reporting.UpdateProgressStage(2, 5);
            lobes_fill = lobes.Copy;
            lobes_fill.ChangeRawImage(lobes.RawImage == PTKColormapLabels.RightLowerLobe);
            lobes_fill.CropToFit;
            lobes_fill.AddBorder(2);
            airway_generator.GrowTree(lobes_fill, lower_right_start_segment, reporting)

            reporting.ShowProgress('LEFT - upper');
            reporting.UpdateProgressStage(3, 5);
            lobes_fill = lobes.Copy;
            lobes_fill.ChangeRawImage(lobes.RawImage == PTKColormapLabels.LeftUpperLobe);
            lobes_fill.CropToFit;
            lobes_fill.AddBorder(2);
            airway_generator.GrowTree(lobes_fill, upper_left_start_segment, reporting)

            reporting.ShowProgress('LEFT - lower');
            reporting.UpdateProgressStage(4, 5);
            lobes_fill = lobes.Copy;
            lobes_fill.ChangeRawImage(lobes.RawImage == PTKColormapLabels.LeftLowerLobe);
            lobes_fill.CropToFit;
            lobes_fill.AddBorder(2);
            airway_generator.GrowTree(lobes_fill, lower_left_start_segment, reporting)

            % Compute radius values based on Strahler orders
            airway_generator.AirwayTree.ComputeStrahlerOrders;

%             airway_generator.AirwayTree.GenerateBranchParameters;

            % Add values of tissue density
            density = dataset.GetResult('PTKDensityAverage');
            airway_generator.AirwayTree.AddDensityValues(density.DensityAverage);

            results = [];
            results.Airways = airway_generator.AirwayTree;

            count_terminal_points = results.Airways.CountTerminalBranches;
            disp(['Number of terminal points: ' int2str(count_terminal_points)]);

            count_branches = results.Airways.CountBranches;
            disp(['Number of branches: ' int2str(count_branches)]);

            results.InitialImage = airway_generator.InitialApexImage;
            if isempty(results.InitialImage)
                results.InitialImage = zeros(template.ImageSize, 'uint8');
            end

        end

        function results = GenerateImageFromResults(airway_results, image_templates, reporting)
            template_image = image_templates.GetTemplateImage(PTKContext.LungROI);

            % Visualising the entire airway tree is Matlab is slow if the number
            % of branches is greater than a few thousand
            % PTKVisualiseAirwayGrowingTree(airway_results.Airways, reporting);

            results = template_image;
            results.ChangeRawImage(zeros(results.ImageSize, 'uint8'));
            results = PTKDrawAirwayGrowingBranchesAsSegmentation(airway_results.Airways, template_image, reporting);
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function checkopt(inopt, dfopt)

% checkopt -- Check options structure for unrecognized fields
%
% Usage:
%       checkopt(inopt, dfopt)
%
% Input:
%       inopt         Input options structure
%       dfopt         Default options structure
%
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2015-06-17
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'License' file distributed with
% the library.

if ~isempty(inopt) && (~isfield(inopt,'NoOptionCheck') || ~inopt.NoOptionCheck),
  ifnc = fieldnames(inopt);
  for k = 1:length(ifnc),
    ifn = ifnc{k};
    if ~strcmp(ifn, 'NoOptionCheck') && ~isfield(dfopt, ifn),
      warning('SPORCO:UnknownOption', 'Unknown option field %s', ifn);
    end
  end
end

return

author: thfylsty


text:
function clean_dir(dir_arg)
% CLEAN_DIR Clean an entire directory of models
%
%
% Inputs:
%   dir_arg  argument for `dir` like '*.OBJ'
%

%addpath ../
files = dir(dir_arg);
for f = 1:numel(files);
  file = files(f).name;
  [path,name,ext] = fileparts(file);
  fprintf('%s...',file);
  [V,F] = load_mesh(file,'Quiet',true);
  mesh_fix_success = true;
  try
    [CV,CF] = meshfix(V,F);
  catch
    fprintf('meshfix crashed...');
    mesh_fix_success = false;
  end
  if size(CF,1) < 0.5*size(F,1)
    fprintf('meshfix removed more than half...');
    tsurf(CF,CV);
    axis equal;
    title(file,'FontSize',15);
    drawnow;
    mesh_fix_success = false;
  end
  if mesh_fix_success
    output_name = fullfile(path,[name '_fixed.obj']);
  else
    try
      [CV,CF] = winding_number_clean(V,F,'SurfaceOnly',true);
      output_name = fullfile(path,[name '_solid.obj']);
    catch(e)
      fprintf('failure.\n');
      continue;
    end
  end
  try
    [CCV,CCT,CCF] = tetgen(CV,CF);
    fprintf('success!\n');
  catch
    fprintf('tetgen failed on CV,CF...');
    fprintf('semi-success?\n');
  end
  tsurf(CF,CV);
  axis equal;
  title(file,'FontSize',15);
  drawnow;
  writeOBJ(output_name,CV,CF);
end

end

author: alecjacobson

repo: gptoolbox


text:
function [metList, noMolMetList] = mol2sdf(mets, molfileDir, sdfFileName, includeRs)
% Concatenates molfiles in `molfileDir` into one SDF file.
%
% USAGE:
%
%    [metList, noMolMetList] = mol2sdf(mets, molfileDir, sdfFileName, includeRs)
%
% INPUTS:
%    mets:            `m x 1` cell array of metabolite identifiers (e.g., BiGG
%                     abbreviations)
%    molfileDir:      Path to directory containing molfiles for metabolites
%                     in mets. Molfile names should match the metabolite
%                     identifiers in mets without compartment assignments.
%
% OPTIONAL INPUTS:
%    sdfFileName:     Name of SDF file. Default is `MetStructures.sdf`.
%    includeRs:       {0, (1)}. If 0, variable structures such as R groups and
%                     repeat units will not be included in SDF.
%
% OUTPUTS:
%    metList:         Cell array listing metabolites in SDF.
%    noMolMetList:    Cell array listing metabolites without mol file.
%
% Written output is `sdfFileName.sdf` - SDF with metabolite structures in same order as in
% metList. Metabolite identifiers in the SDF are the same as in `MetList`.
%
% .. Authors:
%       - Ronan M.T. Fleming
%       - Hulda SH, Nov. 2012   Changed output format from CDF to SDF. Renamed function. Simplified code.

if ~strcmp(molfileDir(end),filesep) % Format inputs
    molfileDir = [molfileDir filesep];
end
if ~exist('sdfFileName','var') ||  isempty(sdfFileName)
    sdfFileName = 'MetStructures.sdf';
end
if ~strcmp(sdfFileName(end-3:end),'.sdf')
    sdfFileName = [sdfFileName '.sdf'];
end
sdfFileName = [molfileDir sdfFileName];

mets = reshape(mets,length(mets),1);
omets = mets; % store original metabolite identifiers
if ischar(mets)
    mets = strtrim(cellstr(mets));
end
if isnumeric(mets)
    mets = strtrim(cellstr(num2str(mets)));
end
mets = regexprep(mets,'(\[\w\])$',''); % Remove compartment assignment

if ~exist('includeRs','var')
    includeRs = 1;
end

% List all molfiles in molfile directory
d = dir(molfileDir);
dirbool = cat(1,d.isdir);
d = d(~dirbool);
molfileNames = {d.name};
molfileNames = regexprep(molfileNames,'(\.mol)$','')';

molfileNames = molfileNames(ismember(molfileNames,mets)); % Only include molfiles for metabolites in mets

noMolMetList = unique(mets(~ismember(mets,molfileNames)));
noMolMetList = reshape(noMolMetList,length(noMolMetList),1);

metList = molfileNames;

% List all elements
elements = {'H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', 'K', 'Ca',...
            'Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', 'Y',...
            'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', 'I', 'Xe', 'Cs', 'Ba', 'La', 'Ce',...
            'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir',...
            'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th', 'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm',...
            'Bk', 'Cf', 'Es', 'Fm', 'Md', 'No', 'Lr', 'Rf', 'Db', 'Sg', 'Bh', 'Hs', 'Mt', 'Ds', 'Rg', 'Cn', 'Uut', 'Fl', 'Uup', 'Lv', 'Uus', 'Uuo'};

% Concatenate molfiles into single SDF file
fid = fopen(sdfFileName,'w+'); % Create SDF
first = true;
for m = 1:length(molfileNames)
    molfile = fileread([molfileDir molfileNames{m} '.mol']); % Read in molfile
    molfile = reshape(molfile, 1, length(molfile)); % for compatibility with Octave (where fileread returns a column vector of chars)
    molfile = regexprep(molfile,'\r',''); % Remove carriage returns
    molfile = regexprep(molfile,'[^\n]*\n',sprintf('%s\n',molfileNames{m}),'once'); % Replace top line with metabolite identifier
    molfile = regexprep(molfile,'\n+$',''); % Ensure there is only one newline caracter at end of each molfile in the SDF

    if includeRs == 0
        % Get atom list
        lines = regexp(molfile,'\n','split');
        try
            atomCount = str2double(lines{4}(1:3));
        catch
            disp('Failure to read mol file:')
            disp(lines)
            continue;
        end
        if atomCount > 0
            atomBlock = lines(5:5 + atomCount - 1);
            pat = '[^a-z_A-z\*]+(?<atom>[a-z_A-Z\*]+)[^a-z_A-z\*]+';
            atoms = regexp(atomBlock,pat,'names');
            atoms = [atoms{:}]';
            atoms = {atoms.atom}';

            % If structure contains a non-elemental atom or repeat units it
            % will not be appended to SDF
            if ~all(ismember(atoms,elements)) || ~isempty(strmatch('M  STY',lines))
                metList = setdiff(metList,molfileNames{m});
                continue;
            end
        end
    end

    % Add molfile to SDF
    if first
        fprintf(fid,'%s\n',molfile);
        first = false;
    else
        fprintf(fid,'%s\n','$$$$');
        fprintf(fid,'%s\n',molfile);
    end
end
fclose(fid);

noMolFileCount = sum(~ismember(mets,molfileNames'));
fprintf('Percentage of metabolites without mol files: %.1f%%\n', 100*noMolFileCount/length(mets));

metList = reshape(metList,length(metList),1);
if isnumeric(omets)
    metList = str2double(metList);
end

author: opencobra

repo: cobratoolbox


text:
function [At,b,c,K,prep,origcoeff] = pretransfo(At,b,c,K,pars)

% [At,b,c,K,prep] = pretransfo(At,b,c,K)
%
% PRETRANSFO  Checks data and then transforms into internal SeDuMi format.
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi

% Nearly complete rewrite for SeDumI 1.4
% Copyright (c) 2013 Michael C. Grant
%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

% -------------------------------------------------------------------------
% Make sure that all fields exist in K, and verify that they are valid
% -------------------------------------------------------------------------

if ~isfield(K,'f') || isempty(K.f)
    K.f = 0;
elseif numel(K.f) ~= 1 || K.f ~= floor(K.f) || K.f < 0 || ~isreal(K.f)
    error('K.f should be nonnegative integer')
end
if ~isfield(K,'l') || isempty(K.l)
    K.l = 0;
elseif K.l ~= floor(K.l) || K.l < 0 || ~isreal(K.l)
    error('K.l should be nonnegative integer')
end
if ~isfield(K,'q') || ~nnz(K.q)
    K.q = zeros(1,0);
else
    K.q = K.q(:)';
    if any(K.q ~= floor(K.q)) || any(K.q<2) || ~isreal(K.q)
        error('K.q should contain only integers bigger than 1')
    end
end
if ~isfield(K,'r') || ~nnz(K.r)
    K.r = zeros(1,0);
else
    K.r = K.r(:)';
    if any(K.r ~= floor(K.r)) || any(K.r<3) || ~isreal(K.r)
        error('K.r should contain only integers bigger than 2')
    end
end
if ~isfield(K,'s') || ~nnz(K.s)
    K.s = zeros(1,0);
else
    K.s = K.s(:)';
    if any(K.s ~= floor(K.s)) || any(K.s<1) || ~isreal(K.s)
        error('K.s should contain only positive integers')
    end
end
% As an alternative to the 'scomplex' flag, I've added a 'z' parameter
% containing a list of Hermitian semidefinite cone sizes. For now, this
% is just translated to 'scomplex' for you. In the future, we may use
% this internally *instead* of scomplex or rsdpN.
if ~isfield(K,'z') || ~nnz(K.z)
    K.z = zeros(1,0);
else
    K.z = K.z(:)';
    if any(K.z ~= floor(K.z)) || any(K.z<1) || ~isreal(K.z)
        error('K.z should contain only positive integers')
    end
end

N_f    = K.f;
N_l    = K.l;
N_fl   = N_f + N_l;
L_q    = length(K.q);
N_q    = sum(K.q);
L_r    = length(K.r);
N_r    = sum(K.r);
N_qr   = N_q + N_r;
L_qr   = L_q + L_r;
L_s    = length(K.s);
L_z    = length(K.z);
L_sz   = L_s + L_z;
N_s    = sum((K.s).^2);
N_z    = sum((K.z).^2);
N_sz   = N_s + N_z;
L_qrsz = L_qr + L_sz;
N_flqr = N_fl + N_qr;
N      = N_flqr + N_sz;

if ~isfield(K,'ycomplex') || isempty(K.ycomplex)
    K.ycomplex = zeros(1,0);
else
    K.ycomplex = sort(K.ycomplex(:))';
    K.ycomplex(find(~diff(K.ycomplex))+1) = [];
    if any(K.ycomplex ~= floor(K.ycomplex)) || any(K.ycomplex<1)
        error('K.ycomplex should contain only positive integers');
    elseif any(K.ycomplex > numel(b))
        error('Elements of K.ycomplex are out of range');
    end
end
if ~isfield(K,'xcomplex') || isempty(K.xcomplex)
    K.xcomplex = zeros(1,0);
else
    K.xcomplex = sort(K.xcomplex(:))';
    K.xcomplex(find(~diff(K.xcomplex))+1) = [];
    if any(K.xcomplex ~= floor(K.xcomplex)) || any(K.xcomplex<1)
        error('K.xcomplex should contain only positive integers');
    elseif any(K.xcomplex>N_flqr)
        error('Elements of K.xcomplex are out of range');
    end
end
if ~isfield(K,'scomplex') || isempty(K.scomplex)
    K.scomplex = zeros(1,0);
else
    K.scomplex = sort(K.scomplex(:))';
    K.scomplex(find(~diff(K.scomplex))+1) = [];
    if any(K.scomplex~=floor(K.scomplex)) || any(K.scomplex<1)
        error('K.scomplex should contain only positive integers');
    elseif any(K.scomplex>L_s)
        error('Elements of K.xcomplex are out of range');
    end
end
if L_z,
    K.s = [ K.s, K.z ];
    K.scomplex = [ K.scomplex, L_s + 1 : L_sz ];
    K.z = zeros(1,0);
    L_s = L_s + L_z;
    N_s = N_s + N_z;
    L_z = 0; %#ok
    N_z = 0; %#ok
end

% -------------------------------------------------------------------------
% Verify the size and validity of At, b, and c
% N = # variables
% -------------------------------------------------------------------------

% SeDuMi assumes that if At is not consistent with K but At' is, that the
% user supplied the transpose of the coefficient matrix. This introduces a
% rare ambiguity in the case where At happens to be square. Past versions
% of SeDuMi would not have allowed this, instead rejecting the case where
% m >= N; however, with complex variables, the situation is not quite as
% clear, and there may technically be cases where m >= N is acceptable.

if ndims(At) > 2 %#ok
    error('A must be a matrix');
elseif nnz(isnan(At)) || nnz(isinf(At)),
    error('A contains NaN or Inf');
elseif size(At,1) == N,
    % nothing
elseif size(At,2) == N,
    At = At';
else
    error('(At,K) size mismatch');
end
if all(size(b)>1)
    error('Parameter b must be a vector');
elseif any(isnan(b)) || any(isinf(b)),
    error('b contains NaN or Inf');
elseif length(b) ~= size(At,2),
    error('(At,b) size mismatch');
else
    b = b(:);
end
if all(size(c)>1)
    error('Parameter c must be a vector');
elseif any(isnan(c)) || any(isinf(c))
    error('c contains NaN or Inf');
elseif length(c) ~= N
    error('(c,K) size mismatch');
else
    c = c(:);
end

% -------------------------------------------------------------------------
% Save the standardized data for further use if needed
% -------------------------------------------------------------------------

if isfield(pars,'errors') && pars.errors==1
    origcoeff.At=At;
    origcoeff.c=c;
    origcoeff.b=b;
    origcoeff.K=K;
else
    origcoeff=[];
end

% -------------------------------------------------------------------------
% Flag diagonal SDP blocks for removal
% -------------------------------------------------------------------------

% There is some serious MATLAB trickery here (if I do say so myself) that
% merits explanation. "spattern" contains the indices of symbolically
% nonzero elements of the dual variable z = c - At * y. This is a single
% vector across all LMI constraints, so "sblk" tells us which indices
% belong to which constraints. The "rem" statement is what determines if
% a particular element is off-diagonal. If there is even one off-diagonal
% element in an SDP, then its corresponding element of "sdiag" is false.
%
% This new version of the analysis replaces the old preprocessSDP(), and
% seems inexpensive enough to apply to all LMIs regardless of size/count.
% The previous version was applied more sparingly.
%
% In theory one could do a more complex analysis of the block structure of
% an LMI, potentially breaking a larger into smaller ones. But this would
% likely be significantly more expensive.

if L_s && ( ~isfield(pars,'sdp') || pars.sdp ),
    ssiz        = (K.s).^2;
    strt        = cumsum([1,ssiz(1:end-1)]);
    sblk        = zeros(1,N_s);
    sblk(strt)  = 1;
    sblk        = cumsum(sblk);
    spattern    = find(c(N_flqr+1:N)~=0|any(At(N_flqr+1:N,:),2))';
    sblk        = sblk(spattern);
    sblk        = sblk(rem(spattern-strt(sblk),K.s(sblk)+1)~=0);
    sdiag       = true(1,L_s);
    sdiag(sblk) = false;
else
    % Even if we disable SDP processing, we're going to move 1x1 SDPs into
    % the nonnegative variable block. It just doesn't make sense to deploy
    % all of that SDP machinery for nonnegative variables.
    sdiag = K.s == 1;
end

% -------------------------------------------------------------------------
% Handle K.ycomplex by splitting apart the complex constraints into pairs
% of real constraints
% -------------------------------------------------------------------------

if ~isempty(K.ycomplex)
    b  = [ real(b) ; imag(b(K.ycomplex)) ];
    At = [ At, 1j * At(:,K.ycomplex) ];
else
    b  = real(b);
end

% -------------------------------------------------------------------------
% Find the locations of the the complex data, so we can convert into
% SeDuMi's internal format, which uses only MATLAB's real representation.
% -------------------------------------------------------------------------
% Strictly speaking, nonnegative variables, the first variable in a Lorentz
% cone, and the first two variables in a rotated Lorentz cone are real. But
% But SeDuMi has allowed them all to be specified as complex; the imaginary
% portions are interpreted as free variables. We have kept that behavior.

% This code replaces whichcpx.c in its entirety. It actually fixes a bug in
% rotated Lorentz cone handling that was probably never exercised.
if isempty(K.xcomplex) && isempty(K.scomplex),
    K.fcplx = zeros(1,0);
    K.qcplx = zeros(1,0);
    K.rcplx = zeros(1,0);
    scplx   = false(1,L_s);
    sreal   = ~sdiag;
    K.rsdpN = L_s;
    N_fc = 0;
    K.cdim = 0;
else
    xc = K.xcomplex;
    tt = xc <= N_fl;
    K.fcplx = xc(tt);
    xc = xc(~tt) -  N_fl;
    tt = xc <= N_q;
    K.qcplx = xc(tt);
    xc = xc(~tt) - N_q;
    tt = xc <= N_r;
    K.rcplx = xc(tt);
    if ~isempty(K.qcplx),
        ndxs = cumsum([1,K.q(1:end-1)]);
        t2 = any(bsxfun(@eq,K.qcplx,ndxs'),1);
        K.fcplx = [ K.fcplx, K.qcplx(t2) + N_fl ];
        K.qcplx(t2) = [];
        t2 = sum(bsxfun(@gt,K.qcplx,ndxs'),1);
        K.q = K.q + full(sparse(1,t2,1,1,L_q));
        K.qcplx = K.qcplx + (1:length(K.qcplx));
        N_q = N_q + length(K.qcplx);
    end
    if ~isempty(K.rcplx),
        ndxs = cumsum([1,K.r(1:end-1)]);
        t2 = any(bsxfun(@eq,K.rcplx,[ndxs,ndxs+1]'),1);
        K.fcplx = [ K.fcplx, K.rcplx(t2) + N_fl + N_q ];
        K.rcplx(t2) = [];
        t2 = sum(bsxfun(@gt,K.rcplx,ndxs'),1);
        K.r = K.r + full(sparse(1,t2,1,1,L_r));
        % This 2*t2 offset is required because QR makes two accesses
        % each of the first two elements of a rotated Lorentz cone.
        K.rcplx = K.rcplx + (1:length(K.rcplx)) + 2 * t2;
        N_r = N_r + length(K.rcplx);
    end
    N_fc = length(K.fcplx);
    N_f  = N_f + N_fc;
    N_fl = N_f + N_l; %#ok
    N_qr = N_q + N_r;
    scplx = false(1,L_s);
    scplx(K.scomplex&~sdiag) = true;
    sreal = ~scplx & ~sdiag;
    K.rsdpN = nnz(sreal);
    K.cdim = length(K.xcomplex) + sum(K.s(scplx).^2);
end

% -------------------------------------------------------------------------
% We have significantly rewritten this section of the code. This section
% constructs a sparse matrix that represents the following transformations:
%   --- Free variables split into differences of nonnegative variables; OR
%   --- Free variables placed in a Lorentz-cone
%   --- Rotated lorentz cones translated to standard Lorentz cones
%   --- Lorentz cones rearranged to trace block + norm-bound blocks
%   --- Conversion of diagonal SDPs to nonnegative variables
%   --- SDP coefficients moved to the lower triangle for increased sparsity
% This code replaces the rotlorenz, qreshape, and vectril MEX files.
% -------------------------------------------------------------------------

newL = 0;
newQ = zeros(1,0);
ii = {}; jj = {}; vv = {};

% Split free variables into the difference of nonnegatives
if ~isfield( pars, 'free' ) || pars.free == 2 && L_qrsz,
    pars.free = 1;
end
if N_f && ~pars.free,
    jt = [ 1 : K.f, K.fcplx ; 1 : K.f, K.fcplx ];
    vt = [ ones(1,K.f), -1j*ones(1,N_fc) ; -ones(1,K.f), 1j*ones(1,N_fc) ];
    ii{end+1} = 1 : 2 * N_f;
    jj{end+1} = jt(:)';
    vv{end+1} = vt(:)';
    newL = 2 * N_f;
    prep.freeL = N_f;
end

% Copy nonnegative variables without change
if K.l,
    ii{end+1} = newL + 1 : newL + K.l;
    jj{end+1} = K.f + 1 : K.f + K.l;
    vv{end+1} = ones(1,K.l);
    newL = newL + K.l;
end

% Convert diagonal SDPs to nonnegative variables
if any(sdiag),
    dsize = K.s(sdiag);
    sdpL = sum(dsize);
    prep.sdiag = dsize;
    jstrt = cumsum([N_flqr+1,K.s(1:end-1).^2]);
    jstrt = jstrt(sdiag);
    istrt = cumsum([1,dsize(1:end-1)]);
    dsize = dsize + 1;
    dblks = cumsum(full(sparse(1,istrt,1,1,sdpL)));
    ii{end+1} = newL + 1 : newL + sdpL;
    jj{end+1} = jstrt(dblks) + dsize(dblks) .* ( ( 1 : sdpL ) - istrt(dblks) );
    vv{end+1} = ones(1,sdpL);
    newL = newL + sdpL;
end

% Stuff free variables into a Lorentz cone
tr_off = newL;
nb_off = newL + L_qr;
if N_f && pars.free,
    tr_off = tr_off + 1;
    nb_off = nb_off + 1;
    ii{end+1} = nb_off + 1 : nb_off + N_f;
    jj{end+1} = [ 1 : K.f, K.fcplx ];
    vv{end+1} = [ ones(1,K.f), -1j*ones(1,N_fc) ];
    nb_off = nb_off + K.f;
    newQ = N_f + 1;
end

% Rearrange Lorentz cones to trace block + norm-bound blocks
if N_q,
    ndxs      = cumsum([1,K.q(1:end-1)]);
    it        = zeros(1,N_q);
    it(ndxs)  = tr_off + 1 : tr_off + L_q;
    it(it==0) = nb_off + 1 : nb_off + ( N_q - L_q );
    jt = K.f + K.l + 1 : K.f + K.l + N_q;
    vt = ones(1,N_q);
    if ~isempty(K.qcplx),
        jt = jt - cumsum(full(sparse(1,K.qcplx,1,1,N_q)));
        vt(K.qcplx) = -1j;
    end
    ii{end+1} = it(:)';
    jj{end+1} = jt;
    vv{end+1} = vt;
    tr_off    = tr_off + L_q;
    nb_off    = nb_off + N_q - L_q;
end

% Transform rotated Lorentz cones to standard Lorentz cones, and rearrange
% to trace block + norm-bound blocks.
if N_r,
    ndxr       = cumsum([1,K.r(1:end-1)]);
    ndxp       = ndxr + 2*(0:L_r-1);
    it         = zeros(1,N_r+2*L_r);
    it(ndxp)   = tr_off + 1 : tr_off + L_r;
    it(ndxp+1) = -1;
    it(ndxp+2) = it(ndxp);
    it(it==0)  = nb_off + 1 : nb_off + ( N_r - L_r );
    it(ndxp+1) = it(ndxp+3);
    jt = [ K.f + K.l + N_q + 1, ones(1,N_r+2*L_r-1) ];
    jt([ndxp+1,ndxp+3]) = 0;
    vt = ones(1,N_r+2*L_r);
    vt([ndxp,ndxp+1,ndxp+2]) = sqrt(0.5);
    vt(ndxp+3) = -sqrt(0.5);
    if ~isempty(K.rcplx),
        jt(K.rcplx) = 0;
        vt(K.rcplx) = -1j;
    end
    ii{end+1} = it;
    jj{end+1} = cumsum(jt);
    vv{end+1} = vt;
    nb_off = nb_off + N_r - L_r;
end

% Replace non-diagonal real SDP coefficients with tril(X) + tril(X',-1).
% This cuts the number of nonzeros approximately in half.
if K.rsdpN,
    dsize = K.s(sreal);
    sdpL  = sum(dsize.^2);
    jstrt = cumsum([N_flqr+1,K.s(1:end-1).^2]);
    jstrt = jstrt(sreal);
    istrt = cumsum([1,dsize(1:end-1).^2]);
    dblks = cumsum(full(sparse(1,istrt,1,1,sdpL)));
    istrt = istrt + nb_off;
    dsize = dsize(dblks);
    istrt = istrt(dblks);
    jndxs = ( nb_off + 1 : nb_off + sdpL ) - istrt;
    cols  = floor(jndxs ./ dsize);
    rows  = jndxs - dsize .* cols;
    ii{end+1} = max(rows,cols) + min(rows,cols) .* dsize + istrt;
    jj{end+1} = jndxs + jstrt(dblks);
    vv{end+1} = ones(1,sdpL);
    nb_off = nb_off + sdpL;
    clear dsize jstrt istrt dblks jndxs rows cols
end

% Replace Hermitian SDP coefficients with tril(X) + tril(X',-1). This one's
% a bit trickier because we have the real and complex values interleaved,
% and the imaginary values along the diagonal are zero.
if K.rsdpN < length(K.s),
    dsize = K.s(scplx);
    jsize = dsize .^ 2;
    sdpL  = 2 * sum(jsize);
    jstrt = cumsum([N_flqr+1,K.s(1:end-1).^2]);
    jstrt = jstrt(scplx);
    bstrt = cumsum([1,2*jsize(1:end-1)]);
    dblks = cumsum(full(sparse(1,bstrt,1,1,sdpL)));
    istrt = bstrt + nb_off;
    dsize = dsize(dblks);
    istrt = istrt(dblks);
    bndxs = ( nb_off + 1 : nb_off + sdpL ) - istrt;
    cols  = floor( bndxs ./ dsize );
    rows  = bndxs - dsize .* cols;
    imgv  = cols >= dsize;
    cols  = cols - imgv .* dsize;
    indxs = max(rows,cols) + min(rows,cols) .* dsize + imgv .* jsize(dblks) + istrt;
    vals  = ( 1 - 2 * ( cols > rows ) ) .* ( 1 - ( 1 + 1j ) .* imgv );
    keep  = ~imgv | ( rows ~= cols );
    jndxs = rows + cols .* dsize + jstrt(dblks);
    ii{end+1} = indxs(keep);
    jj{end+1} = jndxs(keep);
    vv{end+1} = vals(keep);
    clear dsize jsize jstrt bstrt istrt bndxs rows cols vals imgv keep
end

% Update free, nonnegative, and Lorentz variable counts
K.f = 0;
K.l = newL;
K.q = [ newQ, K.q, K.r ];
K.r = zeros(0,1);
K.s = [ K.s(:,~scplx&~sdiag), K.s(scplx&~sdiag) ];
K.rsdpN = nnz(~scplx&~sdiag);
K.N = K.l + sum(K.q) + sum(K.s(1:K.rsdpN).^2)+2*sum(K.s(K.rsdpN+1:end).^2);

% Create the artificial (x0,z0) variable for the self-dual model by
% appending a zero row to At and c. This is accomplished in QR by adding
% 1 to all of the row indices created above.
K.N  = K.N + 1;
K.l  = K.l + 1;
K.m  = length(b);

% Transform At, c
% The transformation matrix QR does not satisfy QR'*QR=I. But it does, in
% fact, serve as a reverse transformation:
%    --- For Lorentz cones, QR applies a permutation; QR' reverses it.
%    --- For rotated Lorentz cones, QR applies a unitary, self-adjoint
%        rotation on the first two variables, so QR' reverses it.
%    --- For split free variables, QR creates the positive and negative
%        parts; QR' combines them back together.
%    --- For free variables placed in a Lorentz cone, QR moves them to the
%        Lorentz cone block, adding an extra epigraph variable. QR' moves
%        them back and drops the extra variable.
%    --- For semidefinte cones, QR adds the strict upper triangle to the
%        strict lower triangle; QR' copies the strict lower triangle to the
%        strict upper triangle, ensuring symmetry.
%    --- QR adds a row for the self-dual variable; QR' removes it.
[dummy,ndxs] = sort(cellfun(@(x) x(1),jj)); %#ok
QR = sparse( horzcat(ii{ndxs})+1, horzcat(jj{ndxs}), horzcat(vv{ndxs}), K.N, length(c) );
At = real( sparse( QR * At ) );
c  = real( sparse( QR * c  ) );
b  = sparse( b );
prep.QR = QR;
clear ii jj vv

% -------------------------------------------------------------------------
% Now K has field K.{l,q,s}
% Generate a more detailed description of cone K:
% Let K.blkstart(i):K.blkstart(i+1)-1 give the index range of block i.
% Compute maxN=max(order), len=sum(order) for LORENTZ, real PSD, herm PSD
% yields: K.{l,q,s,rsdpN,blkstart,rLen,hLen,qMaxn,rMaxn,hMaxn}
% -------------------------------------------------------------------------
Ksr = K.s(1:K.rsdpN);
Ksc = K.s(K.rsdpN+1:end);
K.blkstart = cumsum([K.l+1,length(K.q)+length(K.r),K.q-1,Ksr.^2,2*Ksc.^2]);
K.rLen = sum(Ksr);
K.hLen = sum(Ksc);
K.qMaxn = max([0,K.q]);
K.rMaxn = max([0,Ksr]);
K.hMaxn = max([0,Ksc]);
K.mainblks = K.blkstart(cumsum([1 1 length(K.q)]));
K.qblkstart = K.blkstart(2:2+length(K.q));  % Also include blkend
K.sblkstart = K.blkstart(2+length(K.q):end);
K.lq = K.mainblks(end)-1;


text:
function fdallfig()
    % Plot the correlation integral versus the distance calculated at each grid point.
    % Francesco Pacchiani 1/2000
    %
    %
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    HCIfig=findobj('Type','Figure','-and','Name','Correlation Integral');

    if ~isempty(HCIfig)
        fig = 'addfig';
    else
        fig = 'orifig';
    end

    D = coef(1,1);
    col = 0.6:0.45:3.3;
    colin = jet(64);
    Db = D-0.6;
    Dc = round((Db*64)/2.7);

    switch(fig);

        case 'orifig'

            HCIfig = figure_w_normalized_uicontrolunits('Numbertitle','off','Name','Correlation Integral', 'Visible','on');
            Haxes = gca;

            if D < 0.6 | 3.3 < D

                str5 = 'The fractal dimension calculated is not comprised in the interval 0.6-3.3. The scaling range might be the problem. ';
                msg3 = msgbox(str5, 'Input Error');
                waitforbuttonpress;
                close(msg3);
                close(HCIfig);
            end

            Hline = loglog(r, corint,'color', [colin(Dc,1) colin(Dc,2) colin(Dc,3)]);
            set(Hline,'Linewidth',1);
            xlabel('Interevent Distance R [km]', 'fontsize',12);
            ylabel('Correlation Integral C(R)', 'fontsize',12);
            title('Correlation Integral of All Subsets', 'fontsize',14);
            set(Haxes, 'fontsize', 11);
            set(gca,'NextPlot','add');

        case 'addfig'

            figure(HCIfig);
            set(gcf,'visible','on');
            axes(Haxes);
            set(gca,'NextPlot','add');
            if Dc < 1; Dc = 1; end
            if Dc > 63; Dc = 63; end
            Hline = loglog(r, corint,'color', [colin(Dc,1) colin(Dc,2) colin(Dc,3)]);
            axis([0.001 100 0.000001 5]);
            set(Haxes, 'fontsize', 11);

    end %switch

end

author: CelsoReyes


text:
function [fx,dfdx,dfdp] = f_gen(Xt,Theta,ut,inF)

% Generic evolution function (up to quadratic terms)

% [fx,dfdx,dfdp] = f_gen(Xt,Theta,ut,inF)

deltat = inF.deltat;
n = size(Xt,1);
nc = factorial(n)./(factorial(2)*factorial(n-2));
A = reshape(Theta(1:n^2),n,n);
B = reshape(Theta(n^2+1:n*(2*n+nc)),n,n+nc);

xij = zeros(n+nc,1);
ind = cell(n,1);
ind2 = zeros(n,1);
k = 0;
for i=1:n
    for j=1:n
        if j >= i
            k = k+1;
            xij(k) = Xt(i).*Xt(j);
            if i == j
                ind2(i) = k;
            else
                ind{i} = [ind{i},k];
                ind{j} = [ind{j},k];
            end
        end
    end
end
dbxdx = zeros(n,n);
for i=1:n
    for j=1:n
        dbxdx(i,j) = 2*B(i,ind2(j)).*Xt(j) ...
            + B(i,ind{j})*Xt(setdiff(1:n,j));
    end
end

f = A*Xt + B*xij;
fx = Xt + deltat.*f;
dfdx = eye(n) + deltat*(A + dbxdx)';

dfdp = zeros(n,n*(2*n+nc));
dfdp(:,1:n^2) = kron(Xt',eye(n));
dfdp(:,n^2+1:n*(2*n+nc)) = kron(xij',eye(n));
dfdp = deltat*dfdp';

author: MBB-team

repo: VBA-toolbox


text:
function h = fspecial3(type, sz, param)
% FSPECIAL3  Create predefined 3-dimensional filters
%
% H = FSPECIAL3(TYPE, SZ, PARAM)
%
%   H is a 3-dimensional (3D) filter of TYPE:
%
%     'gaussian'  Rotationally symmetric Gaussian low-pass filter (default)
%                 PARAM: standard deviation in each dimension (default
%                 PARAM = [1.0 1.0 1.0])
%
%   SZ is a vector with the size of the output filter in each dimension, in
%   order [rows, columns, slices]. By default, SZ = [3 3 3]. Sizes have to
%   be odd numbers so that the filter can be centered around 0.
%
% See also: fspecial.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2011 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
error(nargchk(0, 3, nargin, 'struct'));
error(nargoutchk(0, 1, nargout, 'struct'));

% defaults
if (nargin < 1 || isempty(type))
    type = 'gaussian';
end
if (nargin < 2 || isempty(sz))
    sz = [3 3 3];
end
if (nargin < 3 || isempty(param))
    switch type
        case 'gaussian'
            param = [1.0 1.0 1.0];
        otherwise
            error('Filter type not implemented')
    end
end

% check that sizes are odd numbers
if (any(~rem(sz,2)))
    error('Sizes have to be odd numbers')
end

% filter size to each side of 0
l = floor(sz/2);

% compute filter
switch type

    case 'gaussian'

        % filter domain
        [gr, gc, gs] = ndgrid(-l(1):l(1), -l(2):l(2), -l(3):l(3));

        % compute gaussian function
        gr = gr / param(1);
        gc = gc / param(2);
        gs = gs / param(3);
        h = exp(-(gr.*gr + gc.*gc + gs.*gs)*.5);

        % normalize filtered intensity
        h = h / sum(h(:));

    otherwise

        error('Filter type not implemented')
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% This m file contains the list of allowable fields in the
% metabolite_structure
%
% Ines Thiele, 2020/2021

field2Add={
    'VMHId'
    'metNames'
    'charge'
    'neutralFormula'
    'chargedFormula'
    'phenolExplorer'
    'neutralFormula'
    'description'
    'chembl'
    'miriam'
    'hepatonetId'
    'iupac'
    'echa_id' % use https://echa.europa.eu/information-on-chemicals/cl-inventory-database/-/discli/substance/external/echa_id AND https://echa.europa.eu/substance-information/-/substanceinfo/echa_id
    'fda_id'
    'iuphar_id'
    'mesh_id'
    'chodb_id'
    'gtopdb'
    'inchiKey'
    'avgmolweight'
    'monoisotopicweight'
    'keggId'
    'inchiString'
    'pubChemId'
    'cheBIId'
    'hmdb'
    'pdmapName'
    'reconMap'
    'reconMap3'
    'food_db'
    'chemspider'
    'biocyc'
    'biggId'
    'wikipedia'
    'drugbank'
    'seed'
    'metanetx'
    'knapsack'
    'metlin'
    'casRegistry'
    'epa_id'
    'inchiKey'
    'smile'
    'lipidmaps'
    'reactome'
    'GNPS'
    'Recon3D'
    'Agora2'
    'massbank'
    'MoNa'
    'bindingdb'
    'metabolights'
    'rhea'
    'swisslipids'
    'actor'%https://actor.epa.gov/actor/chemical.xhtml?casrn=32981-86-5
    'unii'% https://fdasis.nlm.nih.gov/srs/unii/9G2MP84A8W AND https://druginfo.nlm.nih.gov/drugportal/unii/WDT5SLP0HQ AND https://chem.nlm.nih.gov/chemidplus/unii/WDT5SLP0HQ
    'rxnav' %https://mor.nlm.nih.gov/RxNav/search?searchBy=RXCUI&searchTerm=7052
    'zinc'  %http://zinc.docking.org/substances/ZINC000003812983/
    'pdbeLigand' %https://www.ebi.ac.uk/pdbe-srv/pdbechem/chemicalCompound/show/MOI
    'drugs_com' %https://www.drugs.com/morphine.html
    'chemidplus' %https://chem.nlm.nih.gov/chemidplus/rn/32981-86-5
    'clinicaltrials'    %https://clinicaltrials.gov/search/intervention=Oxypurinol
    'ctd'   %http://ctdbase.org/detail.go?type=chem&acc=D010117 % are they using keggId's? AND http://ctdbase.org/detail.go?type=chem&acc=D010117&view=disease
    'hasmolfile'
    'wikidata'
    'lipidbank'
    'classyFire_Kingdom'
    'classyFire_Superclass'
    'classyFire_Class'
    'classyFire_Subclass'
    'classyFire_Level5'
    'metabolon' % as provided in their input file, currently metabolon_crossmatch_IT_withUpdatedInchiKey
    'exposomeExplorer'%http://exposome-explorer.iarc.fr/search?utf8=%E2%9C%93&query=2-aminophenol+sulfate&button=
    };


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inverse dynamics for the 1dof planar robot
%
%   tau = inversedynamics_1dofplanar(robot, q, qd, qdd, fext)
%
%   Where robot stores the kinematic and dynamic parameters for this robot.
%   q: joint positions.
%   qd: joint velocities.
%   qdd: joint accelerations.
%   fext: vector of external forces. Defined in the last reference system.
%
%   This function just executes the inverse dynamic model for this robot.
%   The equations to compute this dynamic model can be found in:
%   "ROBOT ANALYSIS. The mechanics of Serial and Parallel
%        manipulators". Lung Weng Tsai. John Wiley and Sons, inc. ISBN:
%        0-471-32593-7. page 405.
%
%
%   Author: Arturo Gil Aparicio arturo.gil@umh.es
%   Date: 23/11/2016
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function tau = exercise_inv_dynamics_1dofplanar(robot, q, qd, qdd, g, fext)
a = eval(robot.DH.a);
a1=a(1);

g=abs(g);
m=robot.dynamics.masses(1);

%In this case we must define the Inertia with respect to the rotating axis.
J = (1/3)*m*a1^2;
tau = J*qdd + m*g*a1*cos(q(1))/2;

repo: ARTE


text:
% DEMOILFGPLVM1 Oil data with fully independent training conditional.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'oil';
experimentNo = 1;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('fitc');
options.optimiser = 'scg';
latentDim = 2;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
modelWriteResult(model, dataSetName, experimentNo);

if exist('printDiagram') & printDiagram
  lvmPrintPlot(model, lbls, dataSetName, experimentNo);
end

% Load the results and display them.
lvmScatterPlot(model, lbls);

% compute the nearest neighbours errors in latent space.
errors = lvmNearestNeighbour(model, lbls);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef AddDeltaVAction < AbstractEventAction
    %AddDeltaVAction Summary of this class goes here
    %   Detailed explanation goes here

    properties
        deltaVVect(3,1) double = [0;0;0]; %store as km/s
        frame(1,1) DeltaVFrameEnum = DeltaVFrameEnum.Inertial;
        useDeltaMass(1,1) logical = false;

        optVar AbstractOptimizationVariable
    end

    methods
        function obj = AddDeltaVAction(deltaVVect, frame, useDeltaMass)
            if(nargin > 0)
                obj.deltaVVect = deltaVVect;
                obj.frame = frame;
                obj.useDeltaMass = useDeltaMass;
            end

            obj.id = rand();
        end

        function newStateLogEntry = executeAction(obj, stateLogEntry)
            newStateLogEntry = stateLogEntry;

            if(obj.frame == DeltaVFrameEnum.Inertial)
                dvKmsVect = obj.deltaVVect;

            elseif(obj.frame == DeltaVFrameEnum.OrbitNtw)
                dVVectECI = getNTW2ECIdvVect(obj.deltaVVect, newStateLogEntry.position, newStateLogEntry.velocity);
                dvKmsVect = dVVectECI;

            else
                error('Unknown reference frame found while executing action AddDeltaVAction.');
            end

            newStateLogEntry.velocity = newStateLogEntry.velocity + dvKmsVect;

            if(obj.useDeltaMass)
                [tankMDots, totalThrust, tankStates] = AddDeltaVAction.getTankMDotsAndTotalThrustForStateLogEntry(newStateLogEntry);

                if(abs(sum(tankMDots)) > 0)
                    tankMDotsKgS = tankMDots * 1000;
                    totalMDotKgS = sum(tankMDotsKgS);
                    totalThrustN = totalThrust * 1000;
                    effIsp = totalThrustN / (getG0() * abs(totalMDotKgS)); %sec

                    dvVectMag = norm(dvKmsVect);
                    m0 = newStateLogEntry.getTotalVehicleMass();
                    m1 = revRocketEqn(m0, effIsp, dvVectMag);
                    deltaMassMT = m0 - m1;

                    deltaMassPerTankMT = deltaMassMT * (abs(tankMDots) / abs(sum(tankMDots)));

                    for(i=1:length(tankStates))
                        tankStates(i).setTankMass(tankStates(i).getTankMass() - deltaMassPerTankMT(i));
                    end
                end
            end
        end

        function initAction(obj, initialStateLogEntry)
            %none
        end

        function name = getName(obj)
            name = sprintf('Add Delta-V ([%0.3f %0.3f %0.3f] m/s %s)', obj.deltaVVect(1)*1000, obj.deltaVVect(2)*1000, obj.deltaVVect(3)*1000, obj.frame.nameStr);
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end

        function tf = usesExtremum(obj, extremum)
            tf = false;
        end

        function tf = usesTankToTankConn(obj, tankToTank)
            tf = false;
        end

        function [tf, vars] = hasActiveOptimVar(obj)
            tf = false;
            vars = AbstractOptimizationVariable.empty(0,1);

            if(not(isempty(obj.optVar)))
                tf = any(obj.optVar.getUseTfForVariable());
                vars(end+1) = obj.optVar;
            end
        end

        function data = getUploadDvToKspData(obj, stateLogEntry)
            time = stateLogEntry.time;
            rVect = stateLogEntry.position;
            vVect = stateLogEntry.velocity;

            if(obj.frame == DeltaVFrameEnum.Inertial)
                deltaVNTW = 1000*getNTWdvVect(obj.deltaVVect, rVect(:), vVect(:));

            elseif(obj.frame == DeltaVFrameEnum.OrbitNtw)
                deltaVNTW = 1000*obj.deltaVVect;

            else
                error('Unknown reference frame found while executing action AddDeltaVAction.');
            end

            data(1) = 0;
            data(2) = time;
            data(3) = deltaVNTW(1);
            data(4) = deltaVNTW(2);
            data(5) = deltaVNTW(3);
        end
    end

    methods(Static)
        function addActionTf = openEditActionUI(action, lv)
%             addActionTf = lvd_AddDeltaVActionGUI(action, lv.lvdData);

            output = AppDesignerGUIOutput({false});
            lvd_AddDeltaVActionGUI_App(action, lv.lvdData, output);
            addActionTf = output.output{1};
        end

        function [tankMDots, totalThrust, tankStates] = getTankMDotsAndTotalThrustForStateLogEntry(newStateLogEntry)
            tankStates = newStateLogEntry.getAllActiveTankStates();
            tankStatesMasses = [tankStates.tankMass];
            stageStates = newStateLogEntry.stageStates;
            lvState = newStateLogEntry.lvState;
            ut = newStateLogEntry.time;
            rVect = newStateLogEntry.position;
            vVect = newStateLogEntry.velocity;
            bodyInfo = newStateLogEntry.centralBody;
            steeringModel = newStateLogEntry.steeringModel;

            altitude = newStateLogEntry.altitude;
            pressure = getPressureAtAltitude(bodyInfo, altitude);
            throttle = 1.0;

            powerStorageStates = newStateLogEntry.getAllActivePwrStorageStates();
            storageSoCs = NaN(size(powerStorageStates));
            for(j=1:length(powerStorageStates)) %#ok<*NO4LP>
                storageSoCs(j) = powerStorageStates(j).getStateOfCharge();
            end

            attState = LaunchVehicleAttitudeState();
            attState.dcm = steeringModel.getBody2InertialDcmAtTime(ut, rVect, vVect, bodyInfo);

            [tankMDots, totalThrust, ~, ~] = newStateLogEntry.getTankMassFlowRatesDueToEngines(tankStates, tankStatesMasses, stageStates, throttle, lvState, pressure, ut, rVect, vVect, bodyInfo, steeringModel, storageSoCs, powerStorageStates, attState);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ vX, paramLambda ] = SolveLsNormConst( mA, vB, normConst )
% ----------------------------------------------------------------------------------------------- %
%[ vX ] = SolveLsNormConst( mA, vB, normConst )
% Solves norm constrained Least Squares problem by finding the optimal Dual
% Variable (paramLambda) of the KKT Conditions by succesively solving
% Tikhonov Regularized Least Squares problems.
% The objective Funciton is given by:
% \arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2}
% subject to {\left\| x \right\|}_{2} \leq normSquaredConst
% Input:
%   - mA                -   Input Matrix.
%                           The given matrix of the problem 0.5 || A x - b|| ^ 2.
%                           Structure: Vector (Column Vector).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vB                -   Input Vector.
%                           The given vector of the problem 0.5 || A x - b|| ^ 2.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, ...}.
%   - normConst         -   Norm Constraint.
%                           The solution must obey || x || <= normConst.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: (0, inf).
% Output:
%   - vX                -   Solution Vector.
%                           The optimal solution of the problem 0.5 || A x
%                           - b || ^ 2 subject to || x || <= normConst.
%                           Structure: Vector (Column Vector).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References
%   1.  See - https://stats.stackexchange.com/questions/401212.
% Remarks:
%   1.  T
% TODO:
%   1.  U
% Release Notes:
%   -   1.0.000     14/04/2019
%       *   First realease version.
% ----------------------------------------------------------------------------------------------- %

paramLambda = 0; %<! Initialization

vX = pinv(mA) * vB;
mI = eye(size(vX, 1));

mAA = mA.' * mA;
mAb = mA.' * vB;

if(norm(vX, 2) <= normConst)
    return;
end

hObjFun = @(paramLambda) norm((mAA + (2 * paramLambda * mI)) \ mAb, 2) - normConst;

paramLambda = fzero(hObjFun, paramLambda);

vX = (mAA + (2 * paramLambda * mI)) \ mAb;

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function [ahp, AHP_idp] = idp2ahp(idp)

% IDP2AHP IDP to AHP point transformation.
%   IDP2AHP(IDP) transforms the IDP point into its AHP representation.
%
%   [ahp, AHP_idp] = IDP2AHP(...) returns the Jacobian.

%   Copyright 2010 Joan Sola

x0 = idp(1:3,:);
py = idp(4:5,:);
r = idp(6,:);

if nargout == 1

    ahp = [x0; py2vec(py); r];

else

    if size(idp,2) == 1

        [v,V_py] = py2vec(py);
        ahp = [x0;v;r];
        AHP_idp = [...
            eye(3) zeros(3)
            zeros(3) V_py zeros(3,1)
            zeros(1,5) 1];

    else

        error ('Jacobians not defined for multiple idps.');

    end

end

return

%%
syms x y z e a r real
idp = [x;y;z;e;a;r];
[ahp,AHP_idp] = idp2ahp(idp);

simplify(AHP_idp - jacobian(ahp,idp))

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [ug] = amu2ug(amu)
% Convert mass from atomic mass units to micrograms.
% Chad Greene 2012
ug = amu*1.6605402e-18;

author: Sable

repo: mcbench-benchmarks


text:
function [DT, tE] = compute_CostDataTerm_SENSE(Data, SP3, x, params)
%% Compute the data part of the cost function
% 1/2 * (y-Bx)^H * W * (y-Bx); B = S * F

smap = params.smap;
ncoils = params.ncoils;

%% Compute the Data Term
tS = tic;
Sx = smap .* repmat(x, [1, 1, ncoils]);
FSx = fft2(Sx) .* SP3;
dif = Data - FSx;

DT = sum(abs(dif(:)).^2) / 2;
tE = toc(tS);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [x,g,j,gg] = v_kmeans(d,k,x0,l)
%V_KMEANS Vector quantisation using K-means algorithm [X,ESQ,J]=(D,K,X0,L)
%
%  Inputs:
%
%    D(N,P)  contains N data vectors of dimension P
%    K       is number of centres required
%    X0(K,P) are the initial centres (optional)
%
%      or alternatively
%
%    X0      gives the initialization method
%            'f'   pick K random elements of D as the initial centres [default]
%            'p'   randomly divide D into K sets and choose the centroids
%    L       gives max number of iterations (use 0 if you just want to calculate G and J)
%
%  Outputs:
%
%    X(K,P)  is output row vectors (omitted if L=0)
%    G       is mean square error
%    J(N)    indicates which centre each data vector belongs to
%    GG(L)   gives the mean square error at the start of each iteration (omitted if L=0)
%
% It is often a good idea to scale the input data so that it has equal variance in each
% dimension before calling V_KMEANS.

%  Originally based on a routine by Chuck Anderson, anderson@cs.colostate.edu, 1996

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: v_kmeans.m 4497 2014-04-23 10:28:55Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

memsize=voicebox('memsize');
[n,p] = size(d);
nb=min(n,max(1,floor(memsize/(8*p*k))));    % block size for testing data points
nl=ceil(n/nb);                  % number of blocks
if nargin<4
    l=300;                  % very large max iteration count
    if nargin<3
        x0='f';             % use 'f' initialization mode
    end
end
if ischar(x0)
    if k<n
        if any(x0)=='p'                  % Initialize using a random partition
            ix=ceil(rand(1,n)*k);       % allocate to random clusters
            ix(rnsubset(k,n))=1:k;      % but force at least one point per cluster
            x=zeros(k,p);
            for i=1:k
                x(i,:)=mean(d(ix==i,:),1);
            end
        else                                % Forgy initialization: choose k random points [default]
            x=d(rnsubset(k,n),:);         % sample k centres without replacement
        end
    else
        x=d(mod((1:k)-1,n)+1,:);    % just include all points several times
    end
else
    x=x0;
end
m=zeros(n,1);           % minimum distance to a centre
j=zeros(n,1);           % index of closest centre
gg=zeros(l,1);
wp=ones(1,p);
kk=1:p;
kk=kk(ones(n,1),:);
kk=kk(:);

if l>0
    for ll=1:l                 % loop until x==y causes a break

        % find closest centre to each data point [m(:),j(:)] = distance, index

        ix=1;
        jx=n-nl*nb;
        for il=1:nl
            jx=jx+nb;        % increment upper limit
            ii=ix:jx;
            z = disteusq(d(ii,:),x,'x');
            [m(ii),j(ii)] = min(z,[],2);
            ix=jx+1;
        end
        y = x;              % save old centre list

        % calculate new centres as the mean of their assigned data values (or zero for unused centres)

        nd=full(sparse(j,1,1,k,1));         % number of points allocated to each centre
        md=max(nd,1);                       % remove zeros
        jj=j(:,wp);
        x=full(sparse(jj(:),kk,d(:),k,p))./md(:,wp);    % calculate the new means
        fx=find(nd==0);

        % if any centres are unused, assign them to data values that are not exactly on centres
        % choose randomly if there are more such points than needed

        if ~isempty(fx)
            q=find(m~=0);
            if length(q)<=length(fx)
                x(fx(1:length(q)),:)=d(q,:);
            else
                if length(fx)>1
                    [rr,ri]=sort(rand(length(q),1));
                    x(fx,:)=d(q(ri(1:length(fx))),:);
                else
                    x(fx,:) = d(q(ceil(rand(1)*length(q))),:);
                end
            end
        end

        % quit if the centres are unchanged

        gg(ll)=sum(m,1);
        if x==y
            break
        end
    end
    gg=gg(1:ll)/n;
%     ll % *** DEBUG ***
%     gg' % *** DEBUG ***
    g=gg(end);
else            % if l==0 then just calculate G and J (but rename as X and G)
    ix=1;
    jx=n-nl*nb;
    for il=1:nl
        jx=jx+nb;        % increment upper limit
        ii=ix:jx;
        z = disteusq(d(ii,:),x,'x');
        [m(ii),j(ii)] = min(z,[],2);
        ix=jx+1;
    end
    x=sum(m,1)/n;
    g=j;
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function [prob,sol,fmin] = bilp_prob(varargin)
%BILP_PROB  Return an OPTI BILP
%
%   prob = bilp_prob(no) return a pre-built optiprob of a saved BILP.
%
%   [prob,sol,fmin] = bilp_prob(no) returns the optimum solution and
%   function eval at the optimum
%
%   no = bilp_prob() returns the number of problems available for testing.

%   (C) 2011 Jonathan Currie (I2C2)

%Check if just returning no problems
if(nargin < 1)
    prob = 2; sol = []; fmin = [];
    return;
else
    no = varargin{1};
end

%Big switch yard
switch(no)
    case 1
        f = -[6 5]';
        A = [-3,5; 6,4; 3, -5; -6, -4];
        b = [6;9;1;3];
        xint = 'BB';
        prob = optiprob('f',f,'ineq',A,b,'int',xint);
        sol = [0;1];
        fmin = -5;

    case 2
        f = -[9 5 6 4]';
        A = [6 3 5 2; 0 0 1 1; -1 0 1 0; 0 -1 0 1];
        b = [9; 1; 0; 0];
        xint = 'BBBB';
        prob = optiprob('f',f,'ineq',A,b,'int',xint);
        sol = [1;1;0;0];
        fmin = -14;

    otherwise
        error('Problem not available or not implemented yet');
end


text:
function check = uniform_check ( a, b )

%*****************************************************************************80
%
%% UNIFORM_CHECK checks the parameters of the Uniform CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    A < B.
%
%    Output, logical CHECK, is true if the parameters are legal.
%
  if ( b <= a )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'UNIFORM_CHECK - Fatal error!\n' );
    fprintf ( 1, '  B <= A.\n' );
    check = 0;
    return
  end

  check = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
g=-50;
Q=1;
f=10e3;

[x, Fs] = wavread('tea_16bit_48kHz_short.wav');

[b a]=get_high_shelving_filter(g,Q,f,Fs);

y=filter(b,a,x);
sound(y,Fs);


text:
function gridlines_test02 ( )

%*****************************************************************************80
%
%% GRIDLINES_TEST02 tests GRID_POLAR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GRIDLINES_TEST02\n' );
  fprintf ( 1, '  GRID_POLAR draws a polar grid.\n' );
  fprintf ( 1, '  Compare two methods of sampling the circle.\n' );

  r = sqrt ( rand ( 400, 1 ) );
  t = pi * rand ( 400, 1 );
  plot ( r .* cos ( t ), r .* sin ( t ), 'b*' );

  grid_polar ( 0.0, 0.0, 0.1, 1.0, 10, 0, 180, 13 );

  hold on

  r = - rand ( 400, 1 );
  t = pi * rand ( 400, 1 );
  plot ( r .* cos ( t ), r .* sin ( t ), 'r*' );

  grid_polar ( 0.0, 0.0, 0.1, 1.0, 10, 180, 360, 13 );

  xlabel ( 'X' )
  ylabel ( 'Y' )
  title ( 'Blue: R = sqrt(rand), Red: R = - rand' )
  axis equal

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Press return to continue.\n' )
  hold off
  pause
  clf

  return
end


text:
function createInstallFile(options)
global installfilename
global platform

platform = [];

% Start with a clean slate
cleanup('','','start');

installfilename = sprintf('%s_install', lower(getAppname()));
[~, exename] = getAppname();

setNamespace(exename)

% Generate a LastCheckForUpdates.dat in case that we  haven't run Homer
% before creating a build
checkForHomerUpdates()

if ~exist('options','var') || isempty(options)
    options = 'all';
end

% Find installation path and add it to matlab search paths
dirnameApp = getAppDir();
if isempty(dirnameApp)
    MessageBox('Cannot create installation package. Could not find root application folder.');
    deleteNamespace(exename)
    return;
end
dirnameInstall = filesepStandard(fileparts(which('createInstallFile.m')));
if isempty(dirnameInstall)
    MessageBox('Cannot create installation package. Could not find root installation folder.');
    deleteNamespace(exename)
    return;
end
addpath(dirnameInstall, '-end')
cd(dirnameInstall);

% Set the executable names based on the platform type
platform = setplatformparams();

if exist([dirnameInstall, installfilename],'dir')
    rmdir_safe([dirnameInstall, installfilename]);
end
if exist([dirnameInstall, installfilename, '.zip'],'file')
    delete([dirnameInstall, installfilename, '.zip']);
end
mkdir([dirnameInstall, installfilename]);
mkdir([dirnameInstall, installfilename, '/FuncRegistry']);
mkdir([dirnameInstall, installfilename, '/FuncRegistry/UserFunctions']);
mkdir([dirnameInstall, installfilename, '/SubjDataSample']);

% Generate executables
if ~strcmp(options, 'nobuild')
    Buildme_Setup();
    Buildme();
    if ~ispc()
        c = str2cell(version(),'.');
        mcrver = sprintf('v%s%s', c{1}, c{2});
        if islinux()
            perl('./makesetup.pl','./run_setup.sh','./setup.sh', mcrver);
        elseif ismac()
            perl('./makesetup.pl','./run_setup.sh','./setup.command', mcrver);
        end
    end
end

dirnameDb2DotMat = findWaveletDb2([dirnameInstall, installfilename]);

% Copy files to installation package folder
for ii = 1:length(platform.exename)
    if exist([dirnameInstall, platform.exename{ii}],'file')
        copyfile([dirnameInstall, platform.exename{ii}], [dirnameInstall, installfilename, '/', platform.exename{ii}]);
    end
end
if exist([dirnameInstall, platform.setup_script],'file')==2
    copyfile([dirnameInstall, platform.setup_script], [dirnameInstall, installfilename]);
end
for ii = 1:length(platform.setup_exe)
    if exist([dirnameInstall, platform.setup_exe{ii}],'file')
        if ispc()
            copyfile([dirnameInstall, platform.setup_exe{1}], [dirnameInstall, installfilename, '/installtemp']);
        else
            copyfile([dirnameInstall, platform.setup_exe{ii}], [dirnameInstall, installfilename, '/', platform.setup_exe{ii}]);
        end
    end
end

if exist([dirnameApp, 'SubjDataSample'],'dir')
    copyfile([dirnameApp, 'SubjDataSample'], [dirnameInstall, installfilename, '/SubjDataSample']);
end

for ii=1:length(platform.createshort_script)
    if exist([dirnameInstall, platform.createshort_script{ii}],'file')
        copyfile([dirnameInstall, platform.createshort_script{ii}], [dirnameInstall, installfilename]);
    end
end

dirnameSrc = filesepStandard(fileparts(which([exename, '.m'])));
cfg = ConfigFileClass();
for ii = 1:length(cfg.filenames)
    p = filesepStandard(fileparts(cfg.filenames{ii}));
    k = strfind(p, dirnameSrc);
    pathRelative = p(k+length(dirnameSrc):end);
    fprintf('Copying  %s  to  %s\n', cfg.filenames{ii}, [dirnameInstall, installfilename, '/', pathRelative]);
    if ~ispathvalid([dirnameInstall, installfilename, '/', pathRelative])
        mkdir([dirnameInstall, installfilename, '/', pathRelative])
    end
    copyfile(cfg.filenames{ii}, [dirnameInstall, installfilename, '/', pathRelative]);
end

if exist([dirnameApp, 'AppSettings.cfg'],'file')
    copyfile([dirnameApp, 'AppSettings.cfg'], [dirnameInstall, installfilename]);
end

if exist([dirnameDb2DotMat, 'db2.mat'],'file')
    copyfile([dirnameDb2DotMat, 'db2.mat'], [dirnameInstall, installfilename]);
end

if exist([dirnameInstall, 'makefinalapp.pl'],'file')
    copyfile([dirnameInstall, 'makefinalapp.pl'], [dirnameInstall, installfilename]);
end

if exist([dirnameInstall, 'generateDesktopPath.bat'],'file')
    copyfile([dirnameInstall, 'generateDesktopPath.bat'], [dirnameInstall, installfilename]);
end

if exist([dirnameApp, 'README.md'],'file')
    copyfile([dirnameApp, 'README.md'], [dirnameInstall, installfilename]);
end

if exist([dirnameApp, 'FuncRegistry/UserFunctions'],'dir')
    copyfile([dirnameApp, 'FuncRegistry/UserFunctions'], [dirnameInstall, installfilename, '/FuncRegistry/UserFunctions']);
end

if ispathvalid([dirnameInstall, 'uninstall.bat'])
    copyfile([dirnameInstall, 'uninstall.bat'], [dirnameInstall, installfilename, '/uninstall.bat']);
end

if exist([dirnameApp, 'SDGcolors.csv'],'file')
    copyfile([dirnameApp, 'SDGcolors.csv'], [dirnameInstall, installfilename]);
end

if exist([dirnameApp, 'Version.txt'],'file')
    copyfile([dirnameApp, 'Version.txt'], [dirnameInstall, installfilename]);
end

if exist([dirnameApp, 'LastCheckForUpdates.dat'],'file')
    copyfile([dirnameApp, 'LastCheckForUpdates.dat'], [dirnameInstall, installfilename]);
end

% Zip it all up into a single installation file
zip([dirnameInstall, installfilename, '.zip'], [dirnameInstall, installfilename]);

% Clean up
deleteNamespace(exename)
fclose all;
cleanup(dirnameInstall, dirnameApp);

author: BUNPC


text:
function isobars

% eeg_topo_murks_contour - illustrate isocontours of a surface

clear; close all;

s = 40;
w = 1/150;
cvec = 10:60:300;

p = -s:s;
[x,y] = meshgrid(p);
f = exp(-w*(x.^2 + y.^2));
f = max(f,0);
f = sqrt(f);

surf(x,y,f); shading interp; hold on; rotate3D

V = (x + 10).^2 + (y + 10).^2 + f.^2;

for k = 1:length(cvec)
  c = cvec(k);
  Vc = (V > c);

%  Vcdil = dilate(Vc); Vcdil = double(Vcdil);

  se = [0 1 0 ; 1 1 1 ; 0 1 0];
%  se = ones(1);
  Vcdil = conv2(Vc,se,'same');
  Vcdil = Vcdil > 0;

  Vc = Vcdil - Vc;
  Vcvec = find(Vc);

  plot3(x(Vcvec),y(Vcvec),f(Vcvec),'w.')

end

hold off

author: PatternRecognition

repo: OpenBMI


text:
% Mostly bottom-up method for generating a set
% of good candidate parses to optimize further.
%
% Input
%  I: [n x n bool] binary image
%  lib: library
%  ninit: number of parses to choose
%  verbose: true/false
%
% Output
%  bestMP: [ninit x 1 cell] list of best candidate motor programs
%  score_sorted: [n x 1] score of all parses considered in sorted order
%
function [bestMP,score_sorted] = generate_random_parses(I,lib,ninit,verbose)
    if ~exist('verbose','var')
       verbose = false;
    end

    % Check that image is in the right format
    assert(UtilImage.check_black_is_true(I));

    % If we have a blank image
    if sum(sum(I))==0
       bestMP = [];
       return
    end

    % Get character skeleton from the fast bottom-up method
    G = extract_skeleton(I,verbose);

    % Create a set of random parses through random walks
    ps = defaultps_bottomup;
    if verbose, fprintf(1,'\ngenerating random walks...\n'); end
    RW = RandomWalker(G);
    PP = ProcessParses(I,lib,verbose);

    % Add deterministic minimum angle walks
    for i=1:ps.nwalk_det
        PP.add(RW.det_walk);
    end

    % Sample random walks until we reach capacity.
    walk_count = PP.nwalks;
    while (PP.nl < ps.max_nstroke) && (walk_count < ps.max_nwalk)
        list_walks = RW.sample(1);
        PP.add(list_walks{1});
        walk_count = walk_count + 1;
        if verbose && mod(walk_count,5)==0
            fprintf(1,'%d,',walk_count);
            if mod(walk_count,20)==0
               fprintf(1,'\n');
            end
        end
    end
    if verbose, fprintf(1,'done.\n'); end

    % Turn parses into motor programs with additional search
    % and processing
    [bestMP,score_sorted] = parses_to_MPs(I,PP,ninit,lib,verbose);

end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
function [Y, spectrum] = sllle_wg(G, d)
%SLLLE_WG Solves the Locally Linear Embedding from weight graph
%
% $ Syntax $
%   - Y = sllle_wg(G, d)
%   - [Y, spectrum] = sllle_wg(G, d)
%
% $ Arguments $
%   - G:        The weights graph
%   - d:        The dimension of the embedding
%   - Y:        The sample coordinates in the embedded space
%   - spectrum: The spectrum of the embeded space dimensions
%
% $ Description $
%   - Y = sllle_wg(G, d) solves the locally linear embedding from
%     a given weight graph. G should be a graph of n nodes, where
%     the edge value from i-th node to j-th node, means the weights
%     on the i-th sample in constructing the j-th sample, (or the
%     construction of i-th sample to the j-th sample).
%     Y is solved by taking the eigenvectors of (I - W)(I - W)^T,
%     corresponding to the (d+1) smallest eigenvalues, and discarding
%     the smallest one.  In our output, the dimension is sorted in the
%     ascending order of eigenvalues.
%
%   - [Y, spectrum] = sllle_wg(G, d) also returns the spectrum of the
%     corresponding dimensions, a column vector of the corresponding
%     eigenvalues.
%
% $ Remarks $
%   - The dimension d should be strictly less than n.
%
% $ History $
%   - Created by Dahua Lin, on Sep 11st, 2006
%

%% parse and verify input

if nargin < 2
    raise_lackinput('sllle_wg', 2);
end

gi = slgraphinfo(G, {'square'});
if ~gi.valued
    error('sltoolbox:invalidarg', 'The graph G should be a valued graph');
end
if isnumeric(G)
    W = G;
else
    W = sladjmat(G, 'sparse', true);
end

n = gi.n;
if d >= n
    error('sltoolbox:invalidarg', 'd should be strictly less than n');
end

%% compute

if issparse(W)
    M = speye(n) - W;
else
    M = eye(n) - W;
end
clear W;
M = M * M';

[spectrum, Y] = slsymeig(M, d+1, 'ascend');
clear M;
spectrum = spectrum(2:d+1);
Y = Y(:, 2:d+1);
Y = Y';

author: lmthang


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

% Script chap::2::script
% Density NIG Model with Gamma Ornstein Uhlenbeck stochastic clock
%
%
%
T = 5;                                      % maturity
f0 = 100;                                  % spot value
r=0;
d=0;
ad = 600;                            % spot value
N = 1024;                            % number of grid points
x = ( (0:N-1) - N/2 ) / ad;          % range

%f = 90:.01:110;                            % range

alpha = 2;
beta = 0;                            % CEV exponent base scenario
delta = .3;
lambda = 3;
a = 1;
b = 1;

legend = 'Base';
title_plot = 'NIG-OU Density';

func = @(x)  cf_nigou(x,T,0,r,d,alpha, beta, delta, lambda,a,b);
y = fftdensity(func,ad,N);
%% Changing a
a_low = .5;
a_high = 2;

func_low = @(x) cf_nigou(x,T,0,r,d,alpha, beta, delta, lambda, a_low,b);
y_low = fftdensity(func_low,ad,N);
func_high = @(x) cf_nigou(x,T,0,r,d,alpha, beta, delta, lambda, a_high,b);
y_high = fftdensity(func_high,ad,N);
legend_low = 'Changing a low';
legend_high = 'Changing a high';

createfigure_density(x,y,y_low,y_high,title_plot,legend,legend_low,legend_high);

%% Changing b
b_low = .5;
b_high = 2;

func_low = @(x) cf_nigou(x,T,0,r,d,alpha, beta, delta, lambda, a,b_low);
y_low = fftdensity(func_low,ad,N);
func_high = @(x) cf_nigou(x,T,0,r,d,alpha, beta, delta, lambda, a,b_high);
y_high = fftdensity(func_high,ad,N);
legend_low = 'Changing b low';
legend_high = 'Changing b high';

createfigure_density(x,y,y_low,y_high,title_plot,legend,legend_low,legend_high);

%% Changing lambda
lambda_low = 1;
lambda_high = 20;

func_low = @(x) cf_nigou(x,T,0,r,d,alpha, beta, delta, lambda_low, a,b);
y_low = fftdensity(func_low,ad,N);
func_high = @(x) cf_nigou(x,T,0,r,d,alpha, beta, delta, lambda_high, a,b);
y_high = fftdensity(func_high,ad,N);
legend_low = 'Changing \lambda low';
legend_high = 'Changing \lambda high';

createfigure_density(x,y,y_low,y_high,title_plot,legend,legend_low,legend_high);


text:
function [denoised, noise] = denoise(noisy_signal,prior,std_noise)
  % Given a noisy image and a similar prior image and the standard deviation of
  % the noise (assumes gaussian distribution) denoise the iamge by building a
  % steerable freq pyramid and using probabilistic algorithms to build a coring
  % function (highest 5 bands only) to remap the image to the original. Also
  % prints and shows plots of the estimators.
  %
  % [denoised, noise] = denoise(noisy_signal,prior,std_noise)
  %
  % Inputs:
  %   noisy_signal  image with noise
  %   prior  similar image without noise
  %   std_noise  the standard deviation of the noise distribution
  % Output:
  %   denoised  the denoised result image
  %   noise  the generated noise image
  %

  % create noise image
  noise = std_noise.*randn(size(noisy_signal));

  % build steerable pyramids for each image
  [prior_pyr,prior_ind]=buildSFpyr(prior,3,3);
  [noisy_signal_pyr,noisy_signal_ind]=buildSFpyr(noisy_signal,3,3);
  [noise_pyr,noise_ind]=buildSFpyr(noise,3,3);

  % reconstruct right away (for debugging)
  %before = reconSFpyr(noisy_signal_pyr,noisy_signal_ind);

  % open a new figure for plotting the estimators
  %screen_size = get(0, 'ScreenSize');
  estimator_figure = figure;
  %set(estimator_figure, 'Position', [0 0 screen_size(3) screen_size(4) ] );
  set(estimator_figure, 'Position', [0 0 1000 200 ] );

  for band_index = 1:5
    % extract top five bands of the prior and noise pyramids
    noise_band = pyrBand(noise_pyr,noise_ind,band_index);
    prior_band = pyrBand(prior_pyr,prior_ind,band_index);
    noisy_signal_band = pyrBand(noisy_signal_pyr,noisy_signal_ind,band_index);

    % small value to add to bins to avoid mayhem
    %EPSILON = 0.1;
    EPSILON = 0.000001;
    % histogram the coefficients
    histogram_bins = -0.3:0.01:0.3;
    histogram_prior = EPSILON+ ...
      hist(reshape(prior_band,prod(size(prior_band)),1),histogram_bins);
    histogram_noise = EPSILON+ ...
      hist(reshape(noise_band,prod(size(noise_band)),1),histogram_bins);
    % compute estimate of bayesian estimator using formula from simoncelli96c.pdf:
    %        ∫ Pn(y-z) z Px(z) dz
    % b(y) = ---------------------
    %         ∫ Pn(y-z) Px(z) dz
    %
    % where Pn is the pdf of noise and Px of image. We estimate Pn with
    % histogram_noise and Px with histogram_prior. Then both top and bottom are
    % convolutions ( z is histogram_bins)
    %
    bayesian_estimator = ...
      conv2(histogram_noise,histogram_bins.*histogram_prior,'same') ./ ...
      conv2(histogram_noise,histogram_prior,'same');

    % plot  the estimator with a straight line for comparison
    subplot(1,5,band_index);
    plot(histogram_bins,bayesian_estimator,histogram_bins,histogram_bins);
    grid on
    %legend('bayesian estimator','straight line for comparison');

    if(band_index==1)
      title('hifi residual');
    else
      title(['hifi orientation ' num2str(band_index-1)])
    end

    % (re)map signal coefficients using bayesian_estimator
    denoised_band = interp1(histogram_bins,bayesian_estimator,noisy_signal_band);

    % add back to pyramid
    noisy_signal_pyr = ...
      setPyrBand(noisy_signal_pyr,noisy_signal_ind,denoised_band,band_index);

  end

  % print the plots
  estimator_figure;
  print('-depsc','bayesian_estimator_plots')

  % reconstruct image from modified pyramid
  denoised = reconSFpyr(noisy_signal_pyr,noisy_signal_ind);

end

author: alecjacobson


text:
clear all; close all;
I=imread('pout.tif');
J=log(im2double(I)+1);
K=fft2(J);
n=5;
D0=0.1*pi;
rh=0.7;
rl=0.4;
[row, column]=size(J);
for i=1:row
    for j=1:column
        D1(i,j)=sqrt(i^2+j^2);
        H(i,j)=rl+(rh/(1+(D0/D1(i,j))^(2*n)));
    end
end
L=K.*H;
M=ifft2(L);
N=exp(M)-1;
figure;
subplot(121);
imshow(I);
subplot(122);
imshow(real(N));

author: HuangCongQing


text:
%DIPIMAGECHECK  Check whether the DIP_IMAGE toolbox is in the path

function dipimagecheck

	if exist('dip_image','file') ~= 3
		error([newline 'The DIPIMAGE package is needed and not found.' ...
      newline 'Please download and install it from http://www.diplib.org/'])
	end

author: marianux

repo: ecg-kit

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