text:
function [x,state] = struct_cell2mat(z,task)
%STRUCT_CELL2MAT Convert the contents of a cell array into a matrix.
%   [x,state] = struct_cell2mat(z) generates the matrix x as the function
%   cell2mat applied to the cell array z. The structure state stores
%   information which is reused in computing the right and left
%   Jacobian-vector products.
%
%   struct_cell2mat(z,task) computes the right or left Jacobian-vector
%   product of this transformation, depending on the structure task. Use
%   the structure state and add the field 'r' of the same shape as z or the
%   field 'l' of the same shape as x to obtain the structure task for
%   computing the right and left Jacobian-vector products
%
%      (dF(:)/dz(:).')*task.r(:) and
%      (dF(:)/dz(:).')'*task.l(:) + conj((dF(:)/dconj(z(:)).')'*task.l(:)),
%
%   respectively. Here, F(z) represents this transormation, (:) signifies
%   vectorization and the derivative w.r.t. z (conj(z)) is a partial
%   derivative which treats conj(z) (z) as constant. The output has the
%   same shape as x or z for the right and left Jacobian-vector products,
%   respectively.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Structured data fusion,"
%       ESAT-SISTA Internal Report 13-177, KU Leuven, 2013.

if nargin < 2, task = []; end

if isempty(task) || (isempty(task.l) && isempty(task.r))
    tmp = cellfun(@size,z,'UniformOutput',false);
    N = max(cellfun(@length,tmp));
    state.dim = cell(1,N);
    for n = 1:N
        s = size(tmp,n);
        idx = mat2cell([ones(s,n-1) (1:s)' ones(s,N-n)],s,ones(1,N));
        idx = sub2ind(size(tmp),idx{:});
        state.dim{n} = cellfun(@(s)s(n),tmp(idx));
    end
    x = cell2mat(z);
elseif ~isempty(task.r)
    x = cell2mat(task.r);
    state = [];
elseif ~isempty(task.l)
    x = mat2cell(task.l,task.dim{:});
    state = [];
end

end

author: andrewssobral

repo: lrslibrary


text:
function wprob_temp = probwTemperature(BandBT,idwt,h_pt)
%PROBTEMPERATURE calculate temperature probability for water.

    % [temperature test (over water)]
    F_wtemp=BandBT(idwt); % get clear water temperature
    clear idwt;
    t_wtemp=prctile(F_wtemp,100*h_pt);
    clear F_wtemp h_pt;
    wprob_temp=(t_wtemp-BandBT)/400;
    clear t_wtemp BandBT;
    wprob_temp(wprob_temp<0)=0;
    %    wTemp_prob(wTemp_prob > 1) = 1;
end

author: GERSL

repo: Fmask

save_path: github-repos/MATLAB/GERSL-Fmask


% Calculate channel taps based on the 4th or 6th OFDM symbol (both are ZC sequences)
%
% There is almost certainly a better way to do this, but this gets the job done
%
% @param zc_seq Frequency domain ZC sequence from OFDM symbol number 4 or 6 (all FFT bins, no cyclic prefix)
% @param sample_rate Sample rate (in Hz) of `zc_seq`
% @param symbol_idx Which symbol (must be 4 or 6) is in the `zc_seq` vector
% @return taps Result of dividing `zc_seq` into a golden reference copy of the selected ZC sequence
function [taps] = calculate_channel(zc_seq, sample_rate, symbol_idx)
    assert(symbol_idx == 4 || symbol_idx == 6, "Symbol index must be 4 or 6");
    fft_size = get_fft_size(sample_rate);

    % The golden reference needs to be in the frequency domain
    gold_seq = fftshift(fft(reshape(create_zc(fft_size, symbol_idx), size(zc_seq))));

%     figure(400);
%     subplot(1, 2, 1);
%     plot(abs(gold_seq).^2)
%     subplot(1, 2, 2);
%     plot(abs(zc_seq).^2)

    taps = gold_seq ./ zc_seq;
end


text:
% JCR library for cretaceous chalk
data{3}=[0.079 0.1
0.1 0.02
0.3 0.015
0.54 0.0
0.68 -0.4
0.747 -0.7];
% fit model to data
pc1=data{3}(:,2)*1e5; % Pa
sw1=data{3}(:,1);
sw_plot=linspace(0,1,5000);
% pc_plot=interp1(sw, pc, sw_plot, 'linear', 'extrap');
% plot(sw, pc, 'o', sw_plot, pc_plot)
% plot(diff(pc_plot)./diff(sw_plot))
% log(pc)=log(pce)-(1/labda)*log((sw-sw0)/(1-sw0-sor))
% ind0=find(pc==0, 1);
% p1=polyfit(sw(1:ind0-1), log(pc(1:ind0-1)), 1)
% plot(sw(1:ind0-1), log(pc(1:ind0-1)), 'o')

% pc=@(sw, cwi, coi, awi, aoi, swc, sor)(cwi./((sw-swc)/(1-swc)).^awi+coi./((1-sw-sor)/(1-sor)).^aoi);
% swc0=0.0789;
% sor0=0.2529;
% labda=2.4;
% f=@(x, sw)pc(sw, x(1), x(2), 1/labda, 1/labda, swc0, sor0);
% f2=@(x, sw)pc(sw, x(1), x(2), x(3), x(4), swc0, sor0);
% fw=@(x, sw)pc(sw, x(1), 0.0 , x(3), x(4), swc0, sor0);
% fo=@(x, sw)pc(sw, 0.0, x(2), x(3), x(4), swc0, sor0);
% % x=lsqcurvefit(f2, [1e3, -1e3, 1.0, 1.0],sw1, pc1)
% x=lsqnonlin(@(x)(sum((f2(x, sw1)-pc1).^2)), [0.5e3, -1e3, 1.0, 1.0])
% fit a spline
% first, do an extrapolation to assign values to sw=0 and sw=1
sw_end = [0,1];
pc_end = interp1(sw1, pc1, sw_end, 'linear', 'extrap');
pc_pp = makima([sw_end(1); sw1; sw_end(2)], [pc_end(1); pc1; pc_end(2)]);
pc = @(sw)ppval(pc_pp, sw);
pcder = @(sw)ppval(fnder(pc_pp, 1), sw);
plot(sw1, pc1, 'o', sw_plot, pc(sw_plot), '--' )%, ...
% sw_plot, fw(x, sw_plot), '--', ...
% sw_plot, fo(x, sw_plot), '-.');
grid;

% Note:
% Looking at the above data, it is better to go with interpolation of the
% data
% note 2: fnder requires spline toolbox
% see this: https://se.mathworks.com/matlabcentral/answers/95194-how-do-i-find-the-derivative-of-a-spline-curve-in-matlab-7-9-r2009b

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
clear all;
rand('seed',0);
randn('seed',0);
fprintf('test mexProximalPathCoding\n');
p=100;
% generate a DAG
G=sprand(p,p,0.02);
G=mexRemoveCyclesGraph(G);
fprintf('\n');

% generate a data matrix
U=randn(p,10);
U=U-mean(U(:));
U=mexNormalize(U);

% input graph
graph.weights=G;
graph.stop_weights=zeros(1,p);
graph.start_weights=10*ones(1,p);

% FISTA parameters
param.num_threads=-1; % all cores (-1 by default)
param.verbose=true;   % verbosity, false by default
param.lambda=0.05; % regularization parameter

fprintf('Proximal convex path penalty\n');
param.regul='graph-path-conv';
tic
[V1 optim]=mexProximalPathCoding(U,graph,param);
t=toc;
num=mexCountConnexComponents(graph.weights,V1(:,1));
fprintf('Num of connected components: %d\n',num);

fprintf('Proximal non-convex path penalty\n');
param.regul='graph-path-l0';
param.lambda=0.005;
tic
[V2 optim]=mexProximalPathCoding(U,graph,param);
t=toc;
num=mexCountConnexComponents(graph.weights,V2(:,1));
fprintf('Num of connected components: %d\n',num);

graph.start_weights=1*ones(1,p);
param.lambda=0.05;
fprintf('Proximal convex path penalty\n');
param.regul='graph-path-conv';
tic
[V1 optim]=mexProximalPathCoding(U,graph,param);
t=toc;
num=mexCountConnexComponents(graph.weights,V1(:,1));
fprintf('Num of connected components: %d\n',num);

fprintf('Proximal non-convex path penalty\n');
param.regul='graph-path-l0';
param.lambda=0.005;
tic
[V2 optim]=mexProximalPathCoding(U,graph,param);
t=toc;
num=mexCountConnexComponents(graph.weights,V2(:,1));
fprintf('Num of connected components: %d\n',num);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function Pbest = UpdatePbest(Pbest,Population)
% Update the local best position of each particle

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

    replace        = ~all(Population.objs>=Pbest.objs,2);
    Pbest(replace) = Population(replace);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%quad_stokes   set up Stokes problem in quadrilateral domain
%   IFISS scriptfile: DJS; 13 May 2006.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
clear
%% define geometry
pde=3; domain=8; enclosed=0;
quad_domain
load quad_grid
%
%% set up matrices
qmethod=default('Q1-Q1/Q1-P0/Q2-Q1/Q2-P1: 1/2/3/4? (default Q1-P0)',2);
qmethod=qmethod-1;
if qmethod==2,
   [x,y,xy,xyp,mp] = q2q1grid(x,y,xy,mv,bound);
   [A,B,Q,G,Bx,By,f,g] = stokes_q2q1(xy,xyp,mv,mp);
elseif qmethod==3,
   [x,y,xy,xyp] = q2p1grid(x,y,xy,mv,bound);
   [A,B,Q,G,Bx,By,f,g] = stokes_q2p1(xy,xyp,mv);
elseif qmethod==0
   [ev,ee,ebound,xyp] = q1q1grid(x,y,xy,mv,bound,mbound);
   [A,B,Q,C,G,Bx,By,f,g] = stokes_q1q1(xy,ev);
elseif qmethod==1
   [ev,ee,ebound,xyp] = q1p0grid(x,y,xy,mv,bound,mbound);
   [A,B,Q,C,G,Bx,By,f,g] = stokes_q1p0(xy,xyp,mv,ev);
end
gohome
cd datafiles
save square_stokes_nobc pde domain qmethod  A B Q f g xy xyp mbound bound x y
save square_stokes_nobc Bx By  -append
if qmethod==1
   save square_stokes_nobc C G ev ee ebound -append
elseif qmethod==0
   save square_stokes_nobc C G ev ee ebound mv enclosed -append
elseif qmethod==2
   save square_stokes_nobc G mv mp -append
else
   save square_stokes_nobc G mv  -append
end
fprintf('system matrices saved in square_stokes_nobc.mat ...\n')

author: johannesgerer

repo: jburkardt-m


function [ newstr ] = m2tstrjoin( cellstr, delimiter )
%M2TSTRJOIN This function joins a cellstr with a separator
%
% This is an alternative implementation for MATLAB's `strjoin`, since that
% one is not available before R2013a.
%
% See also: strjoin

    %TODO: Unify the private `m2tstrjoin` functions
    %FIXME: differs from src/private/m2tstrjoin in functionality !!!

    nElem = numel(cellstr);
    if nElem == 0
        newstr = '';
        return % m2tstrjoin({}, ...) -> ''
    end

    newstr = cell(2,nElem);
    newstr(1,:)         = reshape(cellstr, 1, nElem);
    newstr(2,1:nElem-1) = {delimiter}; % put delimiters in-between the elements
    newstr(2, end)      = {''}; % for Octave 4 compatibility
    newstr = [newstr{:}];

end


text:
function h = p04_fh ( p, varargin )

%*****************************************************************************80
%
%% P04_FH returns a mesh size function for problem 4.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Modified:
%
%    06 February 2006
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, real P(NP,ND), the point coordinates.
%
%    Input, VARARGIN, room for extra arguments.
%
%    Output, real H(NP,1), the mesh size function.
%
  np = size ( p, 1 );
  h = ones ( np, 1 );

  return
end


text:
function r = eq(p,q)
% this function compares the equality condition of x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

if (isa(p, 'FaceVariable')&&isa(q, 'FaceVariable'))
    r=p;
    r.xvalue = p.xvalue==q.xvalue;
    r.yvalue = p.yvalue==q.yvalue;
    r.zvalue = p.zvalue==q.zvalue;
elseif isa(p, 'FaceVariable')
    r=p;
    r.xvalue = p.xvalue==q;
    r.yvalue = p.yvalue==q;
    r.zvalue = p.zvalue==q;
else
    r=q;
    r.xvalue = p==q.xvalue;
    r.yvalue = p==q.yvalue;
    r.zvalue = p==q.zvalue;
end


text:
function varargout = process_spgranger1( varargin )
% PROCESS_SPGRANGER1: Compute the spectral Granger causality between one signal and all the others, in one file.
%
% USAGE:   OutputFiles = process_spgranger1('Run', sProcess, sInputA)

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Bivariate Granger causality (spectral) 1xN';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Connectivity';
    sProcess.Index       = 665;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Connectivity';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data',     'results',  'matrix'};
    sProcess.OutputTypes = {'timefreq', 'timefreq', 'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;

    % === CONNECT INPUT
    sProcess = process_corr1n('DefineConnectOptions', sProcess, 0);
    % === REMOVE EVOKED REPONSE
    sProcess.options.removeevoked.Comment = 'Remove evoked response from each trial';
    sProcess.options.removeevoked.Type    = 'checkbox';
    sProcess.options.removeevoked.Value   = 0;
    sProcess.options.removeevoked.Group   = 'input';
    % === DIRECTION
    sProcess.options.dirlabel.Comment  = 'Direction of the causality:';
    sProcess.options.dirlabel.Type     = 'label';
    sProcess.options.direction.Comment = {'From the selected node (out)', 'To the selected node (in)', 'Both (generates two files)'};
    sProcess.options.direction.Type    = 'radio';
    sProcess.options.direction.Value   = 3;
    % === GRANGER ORDER
    sProcess.options.grangerorder.Comment = 'Maximum Granger model order (default=10):';
    sProcess.options.grangerorder.Type    = 'value';
    sProcess.options.grangerorder.Value   = {10, '', 0};
    % === MAX FREQUENCY RESOLUTION
    sProcess.options.maxfreqres.Comment = 'Maximum frequency resolution:';
    sProcess.options.maxfreqres.Type    = 'value';
    sProcess.options.maxfreqres.Value   = {2,'Hz',2};
    % === HIGHEST FREQUENCY OF INTEREST
    sProcess.options.maxfreq.Comment = 'Highest frequency of interest:';
    sProcess.options.maxfreq.Type    = 'value';
    sProcess.options.maxfreq.Value   = {100,'Hz',2};
    % === OUTPUT MODE
    sProcess.options.outputmode.Comment = {'Save individual results (one file per input file)', 'Concatenate input files before processing (one file)', 'Save average connectivity matrix (one file)'};
    sProcess.options.outputmode.Type    = 'radio';
    sProcess.options.outputmode.Value   = 1;
    sProcess.options.outputmode.Group   = 'output';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputA) %#ok<DEFNU>
    % Input options
    OPTIONS = process_corr1n('GetConnectOptions', sProcess, sInputA);
    if isempty(OPTIONS)
        OutputFiles = {};
        return
    end

    % Metric options
    OPTIONS.Method = 'spgranger';
    OPTIONS.RemoveEvoked = sProcess.options.removeevoked.Value;
    OPTIONS.GrangerOrder = sProcess.options.grangerorder.Value{1};
    OPTIONS.MaxFreqRes   = sProcess.options.maxfreqres.Value{1};
    OPTIONS.MaxFreq      = sProcess.options.maxfreq.Value{1};

    % Computation depends on the direction
    OutputFiles = {};
    if ismember(sProcess.options.direction.Value, [1 3])
        OPTIONS.GrangerDir = 'out';
        OutputFiles = cat(2, OutputFiles, bst_connectivity({sInputA.FileName}, {sInputA.FileName}, OPTIONS));
    end
    if ismember(sProcess.options.direction.Value, [2 3])
        OPTIONS.GrangerDir = 'in';
        OutputFiles = cat(2, OutputFiles, bst_connectivity({sInputA.FileName}, {sInputA.FileName}, OPTIONS));
    end
end

author: brainstorm-tools


text:
function xdot = dto_rhs (t, x, u)

% first-order equations of motion

% required by dto_trap.m

% input

%  t = current time

%  x = current state vector

%  x(1) = r, x(2) = u, x(3) = v

%  u = current control vector

% output

%  xdot = rhs equations of motion

%  xdot(1) = r dot, xdot(2) = u dot, xdot(3) = v dot

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global acc beta

% current control variable

theta = u(1);

% current thrust acceleration

accm = acc / (1.0 - beta * t);

% evaluate equations of motion at current conditions

xdot(1) = x(2);

xdot(2) = (x(3) * x(3) - 1.0 / x(1)) / x(1) + accm * sin(theta);

xdot(3) = -x(2) * x(3) / x(1) + accm * cos(theta);

author: Sable

repo: mcbench-benchmarks


text:
lambda1= 1;    % [cells/(mm^3 day)] Production source rate of CD4+ T-cells
d      = 0.1;  % Death rate of CD4+ T-cells
alpha1  = 1;   % Infection rate
a      = 0.2;
p1     = 1;
p2     = 1;
c1     = 0.03;
c2     = 0.06;
b1     = 0.1;
b2     = 0.01;
q      = 0.5;
eta    = 0.9799;
h      = .1;

author: eurika-kaiser

repo: SINDY-MPC


text:
function [newton_m2,N_unique,newton_m2_unique] = monomialproducts(N,n);
%MONOMIALPRODUCTS  Internal function used for monomial reduction

% Exponents in squared monomials

N_unique = [];
for i = 1:size(N,1)
    newton_m2{i} = [];
    n = size(N{i},1);
    for j = 1:n
        newton_m2{i} = [newton_m2{i};[(1:n)' repmat(j,n,1) N{i}(1:n,:)+repmat(N{i}(j,:),n,1)]];
    end
    % Whoops, double copies of diagonal (we want double copies of non-diagonals though)
    if isempty(newton_m2{i})
        newton_m2_unique{i} = [];
    else
        [dummy,j,dummy2] = uniquesafe(newton_m2{i}(:,1:2),'rows');
        newton_m2{i} = newton_m2{i}(j,:);
        % Extract unique monomial products
        [dummy,j,dummy2] = uniquesafe(newton_m2{i}(:,3:end),'rows');
        newton_m2_unique{i} = newton_m2{i}(j,:);
    end
    N_unique = [N_unique;newton_m2_unique{i}];
end
if ~isempty(N_unique)
    [dummy,j,dummy2] = uniquesafe(N_unique(:,3:end),'rows');
    N_unique = N_unique(j,:);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% This is a test script which demonstrates the usage of the "MofN_TrackInitiatorX" class.
% =========================================================================>

% Load the ground truth data
load('multiple-robot-tracking.mat');

% Plot settings
ShowPlots = 1;              % Set to 0 to hide plots
numTrueTracks = 3;

% Model parameter shortcuts
lambdaV = 10; % Expected number of clutter measurements over entire surveillance region
V = 10^2;     % Volume of surveillance region (10x10 2D-grid)
V_bounds = [0 10 0 10]; % [x_min x_max y_min y_max]
P_D = 0.9;    % Probability of detection
timestep_duration = duration(0,0,1);

%% Models
transition_model = ConstantVelocityX('VelocityErrVariance', 0.0001,...
                                     'NumDims', 2,...
                                     'TimestepDuration', timestep_duration);
measurement_model = LinearGaussianX('NumMeasDims', 2,...
                                    'NumStateDims', 4,...
                                    'MeasurementErrVariance', 2,...
                                    'Mapping', [1 3]);
clutter_model = PoissonRateUniformPositionX('ClutterRate',lambdaV,...
                                            'Limits',[V_bounds(1:2);...
                                                      V_bounds(3:4)]);
detection_model = ConstantDetectionProbabilityX('DetectionProbability',P_D);

% Compile the State-Space model
model = StateSpaceModelX(transition_model,measurement_model,'Clutter',clutter_model, 'Detection', detection_model);

%% Generate DataList
meas_simulator = MultiTargetMeasurementSimulatorX('Model',model);
% meas_simulator.DetectionProbability = 1;
DataList = meas_simulator.simulate(GroundTruthStateSequence);
N = numel(DataList);

%% Base Filter
obs_covar= measurement_model.covar();
PriorState = GaussianStateX(zeros(4,1), transition_model.covar() + blkdiag(obs_covar(1,1), 0, obs_covar(2,2),0));
base_filter = KalmanFilterX('Model', model, 'StatePrior', PriorState);

%% Data Associator
config.ClutterModel = clutter_model;
config.Clusterer = NaiveClustererX();
config.Gater = EllipsoidalGaterX(2,'GateLevel',10)';
config.DetectionModel = detection_model;
jpdaf = JointProbabilisticDataAssocX(config);

%% Track Initiator

% Initiate Data Associator
config_pdaf.ClutterModel = clutter_model;
config_pdaf.DetectionModel = detection_model;
config_pdaf.Gater = EllipsoidalGaterX(2,'GateLevel',5);
config_pdaf.DetectionProbability = P_D;
pdaf = ProbabilisticDataAssocX(config_pdaf);

% Initiate Tag Generator
tag_gen = RandSampleTagGeneratorX(1:10000);

% Prepare initiator parameters
config_ti.TagGenerator = tag_gen;
config_ti.InitFilter = base_filter;
config_ti.DataAssociator = pdaf;
CovarThreshold = 4*PriorState.Covar;
config_ti.CustomDeleteConditionFcn = ...
    @(x,t) t.Filter.StatePosterior.Covar(1,1)>CovarThreshold(1,1) ...
           || t.Filter.StatePosterior.Covar(3,3)>CovarThreshold(3,3);

% Create the track initiator
myti = LLR_TrackInitiatorX(config_ti);

%% Metric Generators
ospa = OSPAX('CutOffThreshold',1,'Order',2);
gospa = GOSPAX('CutOffThreshold',1,'Order',2);

%% START OF SIMULATION
%  ===================>

% Create figure windows
if(ShowPlots)
    img = imread('maze.png');

    % set the range of the axes
    % The image will be stretched to this.
    min_x = 0;
    max_x = 10;
    min_y = 0;
    max_y = 10;

    % make data to plot - just a line.
    x = min_x:max_x;
    y = (6/8)*x;

    figure('units','normalized','outerposition',[0 0 .5 1])
    ax(1) = gca;
end

TrackList = [];
gospa_vals= zeros(N,4);
for k=2:N
    fprintf('Iteration = %d/%d\n================>\n',k,N);

    %% Extract DataList at time k
    MeasurementList = DataList(k);
    timestamp_km1 = DataList(k-1).Timestamp;
    timestamp_k = MeasurementList.Timestamp;
    dt = timestamp_k - timestamp_km1;
    transition_model.TimestepDuration = dt;
    fprintf('Timestamp = %s\n================>\n',timestamp_k);

    %% Process JPDAF
    jpdaf.MeasurementList = MeasurementList;
    jpdaf.TrackList = TrackList;
    jpdaf.predictTracks();
    jpdaf.associate();
    jpdaf.updateTracks();

    %% Perform Track initiation
    [TrackList, TentativeTrackList] = myti.initiateTracks(jpdaf.TrackList, MeasurementList, jpdaf.AssocWeightsMatrix, jpdaf.LikelihoodMatrix);

    [gospa_vals(k,1), gospa_vals(k,2), gospa_vals(k,3), gospa_vals(k,4)]= ...
        gospa.evaluate(GroundTruthStateSequence{k},jpdaf.TrackList);
    %% Plot update step results
    if(ShowPlots)

        cla(ax(1));
        %imagesc(ax(1),[min_x max_x], [min_y max_y], flipud(img));
        hold on;
        if(exist('data_plot','var'))
            delete(data_plot);
        end
        data_inv = measurement_model.finv(MeasurementList.Vectors);
        data_plot = plot(ax(1), data_inv(1,:), data_inv(3,:),'k*','MarkerSize', 10);

        % Plot confirmed tracks
        for j=1:numel(TrackList)
            means = [TrackList{j}.Trajectory.Mean];
            h2 = plot(ax(1), means(1,:),means(3,:),'-','LineWidth',1);
            h2 = plotgaussellipse(TrackList{j}.Filter.StatePosterior.Mean([1 3]),...
                                  TrackList{j}.Filter.StatePosterior.Covar([1 3],[1 3]),...
                                  'Color','r',...
                                  'Axis',ax(1));
        end

        % Plot tentative tracks
        for j=1:numel(TentativeTrackList)
            h2 = plot(ax(1), TentativeTrackList{j}.Filter.StatePosterior.Mean(1),TentativeTrackList{j}.Filter.StatePosterior.Mean(3),'.','LineWidth',1);
            if j==2
                set(get(get(h2,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
            end
            h2 = plotgaussellipse(TentativeTrackList{j}.Filter.StatePosterior.Mean([1 3]),...
                                  TentativeTrackList{j}.Filter.StatePosterior.Covar([1 3],[1 3]),...
                                  'Color','g',...
                                  'Axis',ax(1));
        end

        % set the y-axis back to normal.
        set(ax(1),'ydir','normal');
        str = sprintf('Robot positions (Update)');
        title(ax(1),str)
        xlabel('X position (m)')
        ylabel('Y position (m)')
        axis(ax(1),V_bounds)
        pause(0.01)
    end
end

figure
title("GOSPA");
subplot(2,3,[1 2 3]), plot(1:k,gospa_vals(1:k,1)), title("GOSPA Metric");
subplot(2,3,4), plot(1:k,gospa_vals(1:k,2)), title("GOSPA Localisation");
subplot(2,3,5), plot(1:k,gospa_vals(1:k,3)), title("GOSPA Missed");
subplot(2,3,6), plot(1:k,gospa_vals(1:k,4)), title("GOSPA False");

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
%%%
%{
Copyright (c) 2010, Jake Hughey
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:

    * Redistributions of source code must retain the above copyright
      notice, this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in
      the documentation and/or other materials provided with the distribution

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
%}
%%%

function [sortedStruct index] = nestedSortStruct(aStruct, fieldNamesCell, directions)
% [sortedStruct index] = nestedSortStruct(aStruct, fieldNamesCell, directions)
% nestedSortStruct returns a nested sort of a (one-dimensional) struct array
% (aStruct), and can also return an index vector. The fields by which to sort are
% specified in a cell array of strings fieldNamesCell. Fields must be single numbers or
% logicals, or chars (usually simple strings).
%
% fieldNamesCell can also be a simple string of one fieldname, in which case
% nestedSortStruct will simply call sortStruct. This will be faster than putting the
% single fieldname in a cell array.
%
% directions is an optional argument to specify whether the struct array should be sorted
% in ascending or descending order for the fields. By default, the struct array will be
% sorted in ascending order for each field. If supplied, directions must be
%       1) a single 1 to sort in ascending order for all fields, or
%       2) a single -1 to sort in descending order for all fields, or
%       3) a vector of 1's and -1's, the same length as fieldNamesCell, where the struct
%          array will be sorted in the order specified by directions(ii) for
%          fieldNamesCell(ii).
%
% nestedSortStruct basically converts the struct array to a cell array, then converts
% relevants parts of the cell array to a matrix, on which sortrows is run.
%
% nestedSortStruct will usually be faster than nestedSortStruct2. For nestedSortStruct,
% the speed of sorting is mostly independent of the order of the fieldnames in fieldNamesCell.
%
% For nestedSortStruct2, the order of the fields in fieldNamesCell affects the speed. The
% sooner a field for which most entries in the struct array have unique values will be
% used to sort the struct array (i.e., the earlier that field's location in
% fieldNamesCell), the faster nestedSortStruct2 will be. If a field with mostly unique
% entries is the first field by which the struct array will be sorted, nestedSortStruct2
% could be faster than nestedSortStruct.

%% check struct
if ~isstruct(aStruct)
    error('first input supplied is not a struct.')
end % if

if sum(size(aStruct)>1)>1 % if more than one non-singleton dimension
    error('I don''t want to sort your multidimensional struct array.')
end % if

%% check fieldnames
if ~iscell(fieldNamesCell)
    if isfield(aStruct, fieldNamesCell) % if fieldNamesCell is a simple string of a valid fieldname
        [sortedStruct index] = sortStruct(aStruct, fieldNamesCell);
        return
    else
        error('second input supplied is not a cell array or simple string of a fieldname.')
    end % if isfield
end % if ~iscell

if ~isfield(aStruct, fieldNamesCell)
    for ii=find(~isfield(aStruct, fieldNamesCell))
        fprintf('%s is not a fieldname in the struct.\n', fieldNamesCell{ii})
    end % for
    error('at least one entry in fieldNamesCell is not a fieldname in the struct.')
end % if

%% check classes of fieldnames
fieldFlag = 0;
for ii=1:length(fieldNamesCell)
    fieldEntry = aStruct(1).(fieldNamesCell{ii});
    if ~( ((isnumeric(fieldEntry) || islogical(fieldEntry)) && numel(fieldEntry)==1) || ischar(fieldEntry) )
        fprintf('%s is not a valid fieldname by which to sort.\n', fieldNamesCell{ii})
        fieldFlag = 1;
    end % if
end % for ii

if fieldFlag
    error('at least one fieldname is not a valid one by which to sort.')
end

%% check directions, create if necessary (1 for ascending, -1 for descending)
if nargin < 3 % if directions doesn't exist
    directions = ones(1, length(fieldNamesCell));
else % check directions if it does exist
    if ~(isnumeric(directions) && all(ismember(directions, [-1 1])))
        error('directions, if given, must be a single number or a vector with 1 (ascending) and -1 (descending).')
    end % if ~(...

    if numel(directions)==1
        directions = directions * ones(1, length(fieldNamesCell)); % create vector from single element
    elseif length(fieldNamesCell)~=length(directions)
        error('fieldNamesCell and directions vector are different lengths.')
    end % if numel...
end % if exist...

%% fieldNamesIdx is a vector of the indices of the fields by which to sort
[dummy fieldNamesIdx] = ismember(fieldNamesCell, fieldnames(aStruct));

%% convert the struct to a cell, squeeze makes sure both row and column arrays are sorted properly, transpose for sortrows
aCell = squeeze(struct2cell(aStruct))';

%% sortrows of aCell, using indices from fieldNamesIdx and directions
[sortedCell index] = sortrows(aCell, fieldNamesIdx .* directions);

sortedStruct = aStruct(index); % apply the index to the struct array

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%% Read shape

[X,T] = readOff('../data/meshes/cat1.off');
M = getMeshData(X,T,10); % compute 10 LB eigenfunctions for fun

%% Design two distributions

blurTime = .0001; % if this gets too small, distances get noisy
blurSteps = 10; % was 3

frontVtx = [18962];% 15966];
backVtx = [22553];%26142

p0 = zeros(M.numVertices,1);
for i=1:length(frontVtx)
    p0(frontVtx(i)) = .5 / M.areaWeights(frontVtx(i));
end
p0 = blurOnMesh(p0,M,blurTime,blurSteps);
showDescriptor(M,p0);

p1 = zeros(M.numVertices,1);
for i=1:length(backVtx)
    p1(backVtx(i)) = .5 / M.areaWeights(backVtx(i));
end
p1 = blurOnMesh(p1,M,blurTime,blurSteps);
showDescriptor(M,p1);

%% Set up Gaussian blur function for barycenter

blurTime = .00015; % if this gets too small, distances get noisy
blurSteps = 10;

blur = @(x) applySymmetricKernel(x,M,blurTime,blurSteps); % faster than pre-factored?
blurTranspose = @(x) blur(x);

%% Take the barycenter

p = [p0 p1];
p(p<1e-10) = 1e-10;
nFunctions = 2;

maxEntropy = max(-sum(p.*log(p).*repmat(M.areaWeights,1,2),1));

entropyLimits = [maxEntropy (maxEntropy+1) (maxEntropy+2) (maxEntropy+3) inf];
nEntropies = length(entropyLimits);

euclideanBarycenter = sum(p,2)/nFunctions;
alpha = [1 1];

options = [];
options.niter = 100; % diminishing returns after that...
%options.unit_area_projection = 1;

barycenter = zeros(M.numVertices,nEntropies);
parfor i=1:length(entropyLimits)
    barycenter(:,i) = convolutionalBarycenter(p,alpha,M.areaWeights,blur,blurTranspose,entropyLimits(i),options);
end

save entropyTest.mat

%%

for i=1:nEntropies
    showDescriptor(M,barycenter(:,i));
end

%% Write meshes

clear full
for i=1:nEntropies
    e = entropyLimits(i);
    e = full(e);
    name = sprintf('entropy_barycenter_%g.obj',e);
    writeTexturedObj(name, M, barycenter(:,i), 'scalar_function.mtl');
end

%% Write boundary

writeTexturedObj('p0.obj', M, p0, 'scalar_function.mtl');
writeTexturedObj('p1.obj', M, p1, 'scalar_function.mtl');

author: gpeyre


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

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

% Last Modified by GUIDE v2.5 30-Oct-2014 21:31:32

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

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

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

% Choose default command line output for ma_AerobrakingCdCalcGUI
handles.output = hObject;
handles.ma_MainGUI = varargin{1};

%GUI setup
populateBodiesCombo(getappdata(handles.ma_MainGUI,'celBodyData'), handles.bodiesCombo);
dragModelCombo_Callback(handles.dragModelCombo, [], handles);

% Update handles structure
guidata(hObject, handles);

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

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

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

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

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

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

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

function preAeroRpText_Callback(hObject, eventdata, handles)
% hObject    handle to preAeroSmaText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of preAeroSmaText as text
%        str2double(get(hObject,'String')) returns contents of preAeroSmaText as a double
newInput = get(hObject,'String');
newInput = attemptStrEval(newInput);
set(hObject,'String', newInput);

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

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

function preAeroSmaText_Callback(hObject, eventdata, handles)
% hObject    handle to preAeroSmaText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of preAeroSmaText as text
%        str2double(get(hObject,'String')) returns contents of preAeroSmaText as a double
newInput = get(hObject,'String');
newInput = attemptStrEval(newInput);
set(hObject,'String', newInput);

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

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

function postAeroSmaText_Callback(hObject, eventdata, handles)
% hObject    handle to postAeroSmaText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of postAeroSmaText as text
%        str2double(get(hObject,'String')) returns contents of postAeroSmaText as a double
newInput = get(hObject,'String');
newInput = attemptStrEval(newInput);
set(hObject,'String', newInput);

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

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

function dragCoeffAnsText_Callback(hObject, eventdata, handles)
% hObject    handle to dragCoeffAnsText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of dragCoeffAnsText as text
%        str2double(get(hObject,'String')) returns contents of dragCoeffAnsText as a double

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

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

% --- Executes on button press in calcCdButton.
function calcCdButton_Callback(hObject, eventdata, handles)
% hObject    handle to calcCdButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    set(handles.dragCoeffAnsText,'String','<Computing>'); drawnow;

    errMsg = validateInputs(handles);

    if(isempty(errMsg))
        computeDragCoefficient(handles);
    else
        set(handles.dragCoeffAnsText,'String','');
        msgbox(errMsg,'Errors were found while computing the drag coefficient.','error');
    end

function errMsg = validateInputs(handles)
    errMsg = {};

    rawDouble = str2double(get(handles.preAeroRpText,'String'));
    enteredStr = get(handles.preAeroSmaText,'String');
    numberName = 'Pre-Aerobrake Periapsis Altitude';
    lb = 0;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(rawDouble, numberName, lb, ub, isInt, errMsg, enteredStr);

    rawDouble = str2double(get(handles.preAeroSmaText,'String'));
    enteredStr = get(handles.preAeroSmaText,'String');
    numberName = 'Pre-Aerobrake Semi-major Axis';
    lb = -Inf;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(rawDouble, numberName, lb, ub, isInt, errMsg, enteredStr);

    rawDouble = str2double(get(handles.postAeroSmaText,'String'));
    enteredStr = get(handles.postAeroSmaText,'String');
    numberName = 'Post-Aerobrake Semi-major Axis';
    lb = -Inf;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(rawDouble, numberName, lb, ub, isInt, errMsg, enteredStr);

    rawDouble = str2double(get(handles.incText,'String'));
    enteredStr = get(handles.incText,'String');
    numberName = 'Orbital Inclination';
    lb = 0;
    ub = 180;
    isInt = false;
    errMsg = validateNumber(rawDouble, numberName, lb, ub, isInt, errMsg, enteredStr);

	contents = cellstr(get(handles.dragModelCombo,'String'));
	dragModel = contents{get(handles.dragModelCombo,'Value')};

    if(~strcmpi(dragModel,'Stock'))
        rawDouble = str2double(get(handles.preAeroMass,'String'));
        enteredStr = get(handles.preAeroMass,'String');
        numberName = 'Pre-Aerobraking Mass';
        lb = eps;
        ub = Inf;
        isInt = false;
        errMsg = validateNumber(rawDouble, numberName, lb, ub, isInt, errMsg, enteredStr);
    end

    celBodyData = getappdata(handles.ma_MainGUI,'celBodyData');

    contents = cellstr(get(handles.bodiesCombo,'String'));
	selBody = contents{get(handles.bodiesCombo,'Value')};
    bodyInfo = celBodyData.(strtrim(lower(selBody)));

    preSMA = str2double(get(handles.preAeroSmaText,'String'));
    postSMA = str2double(get(handles.postAeroSmaText,'String'));

    preEnergy = -bodyInfo.gm/(2*preSMA);
    postEnergy = -bodyInfo.gm/(2*postSMA);
    if(postEnergy >= preEnergy)
        errMsg{end+1} = 'The entered post-aerobrake SMA does not represent a loss in orbital energy w.r.t. the pre-aerobrake SMA.  Your orbital energy must decrease through the aerobrake.';
    end

    preRpAlt = str2double(get(handles.preAeroRpText,'String'));
    if(preRpAlt >= bodyInfo.atmohgt)
        errMsg{end+1} = sprintf('Pre-Aerobrake Periapsis is located outside the atmosphere (%f km).  Lower the  pre-aerobrake periapsis altitude to below %f km and above 0 km.', bodyInfo.atmohgt, bodyInfo.atmohgt);
    end

function computeDragCoefficient(handles)
    celBodyData = getappdata(handles.ma_MainGUI,'celBodyData');

    contents = cellstr(get(handles.dragModelCombo,'String'));
	dragModel = contents{get(handles.dragModelCombo,'Value')};

    contents = cellstr(get(handles.bodiesCombo,'String'));
	selBody = contents{get(handles.bodiesCombo,'Value')};
    bodyInfo = celBodyData.(strtrim(lower(selBody)));

    preRp = str2double(get(handles.preAeroRpText,'String')) + bodyInfo.radius;
    preSMA = str2double(get(handles.preAeroSmaText,'String'));
    postSMA = str2double(get(handles.postAeroSmaText,'String'));
    inc = deg2rad(str2double(get(handles.incText,'String')));
    mass = str2double(get(handles.preAeroMass,'String'));

    preEcc = getEccFromRpAndSma(preSMA, preRp);   %we assume that Rp remains constant

    raan = 0; %arbitary here
    arg = 0;  %arbitary here
    tru = -computeTrueAFromRadiusEcc(bodyInfo.radius+bodyInfo.atmohgt, preSMA, preEcc)-0.001;  %NOT ARBITRARY, SET JUST OUTSIDE ATMOSPHERE

    [rVect,vVectPre]=getStatefromKepler(preSMA, preEcc, inc, raan, arg, tru, bodyInfo.gm);
%     [~,vVectPost]=getStatefromKepler(postSMA, postEcc, inc, raan, arg, tru, bodyInfo.gm);
    %%%
    initialState = buildInitialState(rVect, vVectPre, mass, bodyInfo.id);
    aeroFun = @(dragCoeff) aerobrakeFunc(dragCoeff, initialState, dragModel, postSMA, bodyInfo.gm, celBodyData);

    if(strcmpi(dragModel,'Stock'))
        x0 = 10;
    else
        x0 = 10;
    end

    options = optimset('Display','iter');
    set(handles.dragCoeffAnsText,'Enable','off');
    [Cd, ~, exitflag, ~] = fzero(aeroFun, x0, options);
    Cd = abs(Cd);

    if(exitflag < 0)
        set(handles.dragCoeffAnsText,'String', '<Error>');
        set(handles.dragCoeffAnsText,'Enable','on');
    else
        set(handles.dragCoeffAnsText,'String', num2str(Cd, 10));
        set(handles.dragCoeffAnsText,'Enable','on');
    end

function dEnergy = aerobrakeFunc(dragCoeff, initialState, dragModel, postSMADesired, gmu, celBodyData)
    dragCoeff = abs(dragCoeff);

    aerobrake = ma_createAerobrake('NA', dragCoeff, dragModel, 'r', '-', 1.5);
    eventLog = ma_executeAerobrake(initialState, -1, aerobrake, celBodyData);

    newRVect = eventLog(end,2:4);
    newVVect = eventLog(end,5:7);
    [postSMAActual, ~, ~, ~, ~, ~] = getKeplerFromState(newRVect,newVVect,gmu);

    desiredEnergy = getSpecOrbitEnergyFromSma(postSMADesired, gmu);
    actEnergy = getSpecOrbitEnergyFromSma(postSMAActual, gmu);

    dEnergy = actEnergy - desiredEnergy;

function initialState = buildInitialState(rVect, vVect, mass, bodyId)
    rVect = reshape(rVect,1,3);
    vVect = reshape(vVect,1,3);

    initialState = zeros(1,13);
    initialState(1,:) = [0 rVect vVect bodyId mass 0 0 0 1];

% --- Executes on button press in closeButton.
function closeButton_Callback(hObject, eventdata, handles)
% hObject    handle to closeButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    close(handles.ma_AerobrakingCdCalcGUI);

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

function incText_Callback(hObject, eventdata, handles)
% hObject    handle to incText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of incText as text
%        str2double(get(hObject,'String')) returns contents of incText as a double
newInput = get(hObject,'String');
newInput = attemptStrEval(newInput);
set(hObject,'String', newInput);

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

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

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

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

    if(strcmpi(dragModel,'Stock'))
%         set(handles.preAeroMass,'Enable','off');
%         set(handles.cdLabel,'String','');
        set(handles.preAeroMass,'Enable','on');
        set(handles.cdLabel,'String','m^2');
    elseif(strcmpi(dragModel,'FAR'))
        set(handles.preAeroMass,'Enable','on');
        set(handles.cdLabel,'String','m^2');
    elseif(strcmpi(dragModel,'NEAR'))
        set(handles.preAeroMass,'Enable','on');
        set(handles.cdLabel,'String','m^2');
    else
        error(['Invalid drag model in aerobraking DV calculations: ', dragModel]);
    end

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

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

function preAeroMass_Callback(hObject, eventdata, handles)
% hObject    handle to preAeroMass (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of preAeroMass as text
%        str2double(get(hObject,'String')) returns contents of preAeroMass as a double
newInput = get(hObject,'String');
newInput = attemptStrEval(newInput);
set(hObject,'String', newInput);

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

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

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ G ] = gsp_non_uniform(N, r)
%GSP_NON_UNIFORM Create a random graph from non_uniform sampling of points
%   Usage:  G = gsp_non_uniform(N, r)
%           G = gsp_non_uniform(N)
%           G = gsp_non_uniform()
%
%   Input parameters:
%         N     : Number of nodes (default 200)
%         r     : Factor of non-uniformity (default 3)
%                 the histogram of x axis of samples is proportional to
%                 exp(x*r)
%   Output parameters:
%         G     : Graph structure.
%
%
%   Example:::
%
%          G = gsp_non_uniform(200, 1);
%          figure; gsp_plot_graph(G)
%          G = gsp_non_uniform(200, 5);
%          figure; gsp_plot_graph(G)
%
%
% see also: gsp_sensor, gsp_non_uniform_patch

% Author : Vassilis Kalofolias

% Optional input arguments
% if nargin < 2
%     k = 6;
% end

if nargin < 1
    N = 200;
end
if nargin < 2
    r = 3;
end

% sample x exponentially from [0,5]
x = exp(rand(N, 1) * r);
x = lin_map(sort(x), [0, 5]);
% sample y uniformly from [0,1]
y = rand(N, 1);
y = lin_map(y, [0, 1]);

coords = [x, y];

D = gsp_distanz(coords');

Ds = sort(D, 'ascend');
%max_min_dist = max(min(D + realmax * eye(N)));
s2 = mean(Ds(2,:).^2);
%s2 = mean(Ds(5,:).^2) / 9;

W = exp(- D.^2 / s2 / 4);
W(1:N+1:end) = 0;

% keep a neighbour for the most distant
min_kept_w = min(max(W));

%

% Kill all small connections but always have at least one neighbor
W = W .* (W >= min(min_kept_w, 0.1));

% binary?
%W = double(W>0);

G = gsp_graph(W, coords, [0 5 0 1]);
G.type = 'non_uniform';

end

repo: gspbox


text:
close all;
clearvars;
clc;
rng default;

degree = 5;

n = 10;
sigma = 0.25;
[x_train, t_train] = spx.data.synthetic.func.sinusoid('n', n, 'sigma', sigma);

n2 = 10;
[x_test, t_test] = spx.data.synthetic.func.sinusoid('n', n2, 'sigma', 0);

features_train = spx.ml.features.polynomial(x_train, degree);
model = spx.ml.models.linear.LinearRegression;
model.fit(features_train, t_train);
features_test = spx.ml.features.polynomial(x_test, degree);
y_test = model.predict(features_test);

p = polyfit(x_train, t_train, degree);
py_test = polyval(p,x_test);

fprintf('t_test: ');
spx.io.print.vector(t_test, 3);
fprintf('y_test: ');
spx.io.print.vector(y_test, 3);
fprintf('py_test: ');
spx.io.print.vector(py_test, 3);

author: indigits

repo: sparse-plex


text:
function [divTr, divTe] = sample_leaveOneBlockOut(label, block_idx, strat)
%SAMPLE_LEAVEONEBLOCKOUT - Sampling function: block divisions
%
%Synopsis:
%  [DIVTR, DIVTE] = sample_leaveOneBlockOut(LABEL, BLOCK_IDX, <STRAT>)
%
%Arguments:
% LABEL     - Despensable input argument, required for compability reasons
% BLOCK_IDX - [1 x nSamples] array of indices which defines to which block
%             each sample belongs. This can be e.g. mrk.blkno, when mrk is
%             obtained from mrk_evenlyInBlocks
% STRAT     - [1 x 2] bool indicating to stratify the samples by
%             performing a random sampling on the training and/or test
%             sets, default [0 0]
%
%Returns:
% DIVTR     - Partitions of the training set
% DIVTE     - Partitions of the test set

% 2015-11 Matthias Schultze-Kraft

if nargin<3
    strat = [0 0];
end

divTr = cell(1, 1);
divTe = cell(1, 1);
idx_list = unique(block_idx);
for nn = 1:length(idx_list)
    % test set indices
    idx = find(block_idx==idx_list(nn));
    if strat(2)
        idx = stratify(idx,label);
    end
    divTe{1}(nn) = {idx};
    % training set indices
    idx = find(block_idx~=idx_list(nn));
    if strat(1)
        idx = stratify(idx,label);
    end
    divTr{1}(nn) = {idx};
end

function idx = stratify(idx,label)
idx = {idx(logical(label(1,idx))) idx(logical(label(2,idx)))};
Nmin = min([length(idx{1}) length(idx{2})]);
if Nmin==0
    warning('No stratification performed because the number of samples of one class is zero.')
else
    [Nmax,ci] = max([length(idx{1}) length(idx{2})]);
    ri = randperm(Nmax,Nmin);
    idx{ci} = idx{ci}(ri);
end
idx = [idx{:}];

author: bbci


text:
function [hashVal,hexVal]=genHashVal(inputArray,algorithm)
%%GENHASHVAL  Get a hash value of a character string or a 1D array of
%             bytes. That is, integer values from 0 through 255. The hash
%             value can either be a simple integer value, useful for quick
%             comparison of strings, or one can get a cryptographic hash as
%             an an array of integers from 0 to 255 using various
%             algorithms. Cryptographic hashes are also returned has
%             hexadecimal strings.
%
%INPUTS: inputArray A 1-dimensional string of characters or array of
%                   integer values from 0 to 255 for which one desired a
%                   hash code.
%        algorithm  An optional  string specifying the algorithm to use to
%                   generate the hash. Possible values are
%                   'simple' (the default) Use a simple, non-cryptographic
%                            hash function to get a single integer hash for
%                            the input.
%                   'MD2'    Use the MD2 algorithm to get a weak
%                            cryptographic hash, along with its hexidecimal
%                            string.
%                   'MD5'    Use the MD5 algorithm to get a weak
%                            cryptographic hash, along with its hexidecimal
%                            string.
%                   'SHA-1'  Use the SHA-1 algorithm to get a cryptographic
%                            hash along with its hexidecimal string.
%                  'SHA-384' Use the SHA-384 algorithm to get a
%                            cryptographic hash along with its hexidecimal
%                            string.
%                  'SHA-256' Use the SHA-256 algorithm to get a
%                            cryptographic hash along with its hexidecimal
%                            string.
%                  'SHA-512' Use the SHA-512 algorithm to get a
%                            cryptographic hash along with its hexidecimal
%                            string.
%
%OUTPUTS: hashVal  If algorithmis 'simple', then this is an integer hash
%                  code for the given input. Otherwise, it is an array of
%                  bytes representign the has value.
%         hexVal   A hexidecimal string representing the hash value.
%
%The simple algorithm is executed using the hashCode method in Java's
%java.lang.String class. The cryptographic hashes are executed using the
%java.security.MessageDigest class. In Matlab, one can directly access Java
%functions, which is what was done to implement this function.
%
%April 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(~ischar(inputArray)&&any(inputArray~=fix(inputArray))||any(inputArray>255)||any(inputArray<0))
    error('Hash values can only be obtained from character arrays or arrays of positive integer values from 0 through 255');
end

if(nargin<2)
    algorithm='simple';
end

%If the simple aglorithm is chosen. otherwise, a cryptographic has is being
%used.
if(strcmp(algorithm,'simple'))
    %Convert the array to a character string.
    if(~ischar(inputArray))
        inputArray=native2unicode(inputArray);
    end

    myString=java.lang.String(inputArray);
    hashVal=myString.hashCode;
    %Java's integers are always 32-bits.
    hexVal=dec2hex(typecast(int32(hashVal),'uint32'),8);
    return;
end

%The cryptographic has functions work on bytes, not strings.
if(ischar(inputArray))
    inputArray=unicode2native(inputArray);
end

myDigest=java.security.MessageDigest.getInstance(algorithm);

myDigest.update(inputArray);
hashVal=myDigest.digest;
hexVal=javax.xml.bind.DatatypeConverter.printHexBinary(hashVal);
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [ lo, hi ] = p06_box ( m )

%*****************************************************************************80
%
%% P06_BOX returns a bounding box for problem 06.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Output, real LO(M), HI(M), the low and high corners of the box.
%
  r1 = 1.0;
  r2 = 0.5;

  lo(1:m) = [ -r1, -r1 ];
  hi(1:m) = [ +r1, +r1 ];

  return
end


text:
function convertBPEData(xlsFile,outPath)
% convertBPEData.m
% Function to import T1, T1 FS-pre, T1 FS_post sequences for BPE analysis from paths in spreadsheet
%
% iyera@mskcc.org  5/4/18
%--------------------------------------------------------------------------

%% Read spreadsheet
[~,~,X] = xlsread(xlsFile,1);
X = X(2:end,:);

convertedC = {};
planNameC = {};

%% Import images to CERR
init_ML_DICOM;
for n = 1:size(X,1)
    outname = num2str(X{n,1});
    selDirC = {X{n,2},X{n,3},X{n,4}};
    try
        for k = 1:length(selDirC)
            hWaitbar = waitbar(0,'Scanning Directory Please wait...');
            sourceDir = selDirC{k};
            temp(k) = scandir_mldcm(sourceDir, hWaitbar, 1);
            patient.PATIENT(k) = temp(k).PATIENT;
            close(hWaitbar);
        end
        for j = 1:length(patient.PATIENT)
            dcmdirS.(['patient_' num2str(j)]) = patient.PATIENT(j);
        end
        patNameC = fieldnames(dcmdirS);
        selected = 'all';
        if strcmpi(selected,'all')
            combinedDcmdirS = struct('STUDY',dcmdirS.(patNameC{1}).STUDY,'info',dcmdirS.(patNameC{1}).info);
            for i = 2:length(patNameC)
                for j = 1:length(dcmdirS.(patNameC{i}).STUDY.SERIES)
                    combinedDcmdirS.STUDY.SERIES(end+1) = dcmdirS.(patNameC{i}).STUDY.SERIES(j);
                end
            end
            % Pass the java dicom structures to function to create CERR plan
            planC = dcmdir2planC(combinedDcmdirS,'no');
        else
            planC = dcmdir2planC(patient.PATIENT);
        end

        %% Save to CERR files
        saved_fullFileName = fullfile(outPath,[outname,'.mat']);
        save_planC(planC,[], 'passed', saved_fullFileName);
        convertedC{end+1} = sprintf('%0.8s',outname);
        planNameC{end+1} = saved_fullFileName;

    catch e

        convertedC{end+1} = sprintf('%0.8d',outname);
        planNameC{end+1} = ['Failed with error ', e.message];

    end

    clear planC

end

%% Book-keeping
for i=1:length(convertedC)
    xlswrite(fullfile(outPath,'batch_convert_results.xls'),{convertedC{i}},'Sheet1',['A',num2str(i)])
    xlswrite(fullfile(outPath,'batch_convert_results.xls'),{planNameC{i}},'Sheet1',['B',num2str(i)])
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%% CUBESTOKES Stokes equations on the unit cube
%
%   SQUARESTOKE computes CR-P0 approximations of the Stokes equations in
%   the unit cube.
%
% Added by Shuhao Cao. Apr, 2020.
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

close all; clear;

%% Set up
maxIt = 4;
N = zeros(maxIt,1);
h = zeros(maxIt,1);
erruH1 = zeros(maxIt,1);
erruL2 = zeros(maxIt,1);
errp = zeros(maxIt,1);

%% Generate initial mesh
[node,elem] = cubemesh([-0.5,0.5,-0.5,0.5,-0.5,0.5],0.25);
bdFlag = setboundary3(node,elem,'Dirichlet');

%% PDE and options
pde.exactu = @(p)[(1-cos(2*pi*p(:,1))).*sin(2*pi*p(:,2)).*sin(pi*p(:,3)),...
    -(1-cos(2*pi*p(:,2))).*sin(2*pi*p(:,1)).*sin(pi*p(:,3)), 0*p(:,3)];
pde.f = @(p) [-pi^2*(9*cos(2*pi*p(:,1))-5).*sin(2*pi*p(:,2)).*sin(pi*p(:,3))+ p(:,1).^2,...
    pi^2*(9*cos(2*pi*p(:,2))-5).*sin(2*pi*p(:,1)).*sin(pi*p(:,3)), 0*p(:,3)];
pde.g = @(p) zeros(size(p,1),1);
pde.exactp = @(p) p(:,1).^3/3;
pde.g_D = @(p) pde.exactu(p);
% solver
option.solver = 'diag'; % diagonal preconditioner

%% Finite Element Method
for k = 1:maxIt

    [soln,eqn] = Stokes3CRP0(node,elem,bdFlag,pde,option);
    uh = soln.u;
    uhvec = reshape(uh,size(uh,1)/3,3);
    ph = soln.p;
    N(k) = length(uh)+length(ph);
    h(k) = 1./((size(node,1)).^(1/3)-1);
    % compute error;
    uI = pde.exactu((node(eqn.face(:,1),:)+node(eqn.face(:,2),:)+node(eqn.face(:,3),:))/3);
    erruH1(k) = sqrt((uh-uI(:))'*eqn.A*(uh-uI(:)));
    erruL2(k) = getL2error3(node,elem,pde.exactu,uhvec);
    errp(k) = getL2error3(node,elem,pde.exactp,ph);
    fprintf('Number of Dof %d \n', N(k));

    if k < maxIt
        [node,elem,bdFlag] = uniformrefine3(node,elem,bdFlag);
    end
end

%% Plot convergence rates and display error
figure(2);
showrateh3(h,erruH1,1,'-*','| u_I-u_h |_1',...
           h,erruL2,1,'-*','|| u-u_h ||',...
           h,errp,1,'-+','|| p-p_h||');

fprintf('\n');
disp('Table: Error')
colname = {'#Dof','h','|u_I-u_h|_1','||u-u_h||','||p-p_h||'};
disptable(colname,N,[],h,'%0.3e',erruH1,'%0.5e',erruL2,'%0.5e',errp,'%0.5e');

repo: ifem


text:
function cdmfile(file)
    [pathstr,fname,ext]=fileparts(which(file));
    cd (pathstr);
end


text:
function [dtdx,dtrp,vart,stat]=model_trop(time,pos,azel,rtk,x,nav) %#ok

global glc;
stat=0;
trp=zeros(3,1); dtdx=zeros(3,1); ERR_SAAS=0.3;

if rtk.opt.tropopt==glc.TROPOPT_SAAS
    dtrp=trop_saa(pos,azel,0.7); %humi=0.7
    vart=ERR_SAAS^2;
    stat=1; return;
end

if rtk.opt.tropopt==glc.TROPOPT_EST||rtk.opt.tropopt==glc.TROPOPT_ESTG
    if rtk.opt.tropopt==glc.TROPOPT_EST,trp(1)=x(rtk.it+1);
    else,trp(1:3)=x(rtk.it+1:rtk.it+3);
    end
    [dtdx,dtrp,vart]=trop_model_prec(time,pos,azel,trp);
    stat=1;return;
end

return

repo: GINav


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create JSON using MATLAB for MEG BIDS:
%
% This template is for MEG data of any kind,
% including but not limited to task-based, resting-state, and noise recordings.
%
% If multiple Tasks were performed within a single Run,
% the task description can be set to "task-multitask".
%
% The _meg.json SHOULD contain details on the Tasks.
%
% Some manufacturers data storage
% conventions use folders which contain data files of various nature:
% for example: CTF's .ds format, or 4D/BTi.
%
% Please refer to Appendix VI for examples from a selection of MEG manufacturers
%
% By @Cofficer, Created 9/03/2018
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Writing json files relies on the JSONio library
% https://github.com/bids-standard/bids-matlab
%
% Make sure it is in the matab/octave path
try
    bids.bids_matlab_version;
catch
    warning('%s\n%s\n%s\n%s', ...
            'Writing the JSON file seems to have failed.', ...
            'Make sure that the following library is in the matlab/octave path:', ...
            'https://github.com/bids-standard/bids-matlab');
end

clear;

this_dir = fileparts(mfilename('fullpath'));
root_dir = fullfile(this_dir, '..', filesep, '..');

project = 'templates';

sub_label = '01';
ses_label = '01';
task_label = 'FullExample';
acq_label = 'CTF';
run_label = '1';

% A "proc" (processed) label has been added, especially useful for files coming from Maxfilter
% (for example: sss, tsss, trans, quat, mc, etc.)
proc_label = 'sss';

name_spec.modality = 'meg';
name_spec.suffix = 'meg';
name_spec.ext = '.json';
name_spec.entities = struct('sub', sub_label, ...
                            'acq', acq_label, ...
                            'ses', ses_label, ...
                            'task', task_label, ...
                            'run', run_label, ...
                            'proc', proc_label);

% using the 'use_schema', true
% ensures that the entities will be in the correct order
bids_file = bids.File(name_spec, 'use_schema', true);

% Contrust the fullpath version of the filename
json_name = fullfile(root_dir, project, bids_file.bids_path, bids_file.filename);

%% Adding metadata
% Assign the fields in the Matlab structure that can be saved as a json.
% The following fields must be defined:

% to get the definition of each metadata,
% you can use the bids.Schema class from bids matlab
% For example
schema = bids.Schema;
def = schema.get_definition('TaskName');
fprintf(def.description);

json.TaskName = task_label;
json.InstitutionName = '';
json.InstitutionAddress = '';
json.Manufacturer = '';
json.ManufacturersModelName = '';
json.SoftwareVersions = '';
json.TaskDescription = '';
json.Instructions = '';
json.CogPOID = 'http://www.cogpo.org/ontologies/CogPOver1.owl#FIXME';
json.DeviceSerialNumber = '';

%%
% The following MEG specific fields must also be defined:
json.SamplingFrequency = '';
json.PowerLineFrequency = '';
json.DewarPosition = '';
json.SoftwareFilters = '';
json.DigitizedLandmarks = '';

json.DigitizedHeadPoints = '';

%%
% The following fields should be present for MEG:
json.MEGChannelCount = [];
json.MEGREFChannelCount = [];
json.EEGChannelCount = [];
json.ECOGChannelCount = [];
json.SEEGChannelCount = [];
json.EOGChannelCount = [];
json.ECGChannelCount = [];
json.EMGChannelCount = [];
json.MiscChannelCount = [];
json.TriggerChannelCount = [];
json.RecordingDuration = [];
json.RecordingType = '';
json.EpochLength = '';
json.ContinuousHeadLocalization = '';
json.HeadCoilFrequency = [];
json.MaxMovement = '';
json.SubjectArtefactDescription = '';
json.AssociatedEmptyRoom = '';

%%
% Specific EEG fields if recorded with MEG
% SHOULD be present:
json.EEGSamplingFrequency = [];
json.EEGPlacementScheme = '';
json.ManufacturersAmplifierModelName = '';
json.EEGReference = '';

%% Write JSON
% Make sure the directory exists
bids.util.mkdir(fileparts(json_name));
bids.util.jsonencode(json_name, json);

author: bids-standard

repo: bids-starter-kit


text:
function [incorrectRxns,incorrectMets] = validateCompartments(model)
% This function checks if there are any reactions with incorrect
% compartments in the refined reconstructions. Appropriate compartments
% include [c], [e], and [p].
%
% USAGE:
%
%    [incorrectRxns,incorrectMets] = validateCompartments(model)
%
% INPUTS
% model             COBRA model structure
%
% OUTPUTS
%
% incorrectRxns     Reactions including metabolites in inappropriate compartments
% incorrectMets     Metabolites located in inappropriate compartments
%
% AUTHOR:
%   Almut Heinken, 10/2020

findIncorrectMets = find(~contains(model.mets,{'[c]','[e]','[p]'}));
if ~isempty(findIncorrectMets)
    incorrectMets = model.mets(findIncorrectMets);
    incorrectRxns = {};

    for i=1:length(incorrectMets)
        [Rxns]=findRxnsFromMets(model,incorrectMets{i});
        incorrectRxns=union(incorrectRxns,Rxns);
    end

else
    incorrectMets = {};
    incorrectRxns = {};
end

end

author: opencobra

repo: cobratoolbox


text:
function varargout = StrainMeasurement(varargin)

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

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

function StrainMeasurement_OpeningFcn(hObject, eventdata, handles, varargin)

handles.output = hObject;
guidata(hObject, handles);

global Status Parameters Files

%ENTER DEFAULT DIRECTORY:
Files.DefaultDir = 'C:\';

Files.OldDir = cd;

Status.Loop = 2;
while Status.Loop ~= 0
    if Status.Loop == 2
        Files.AFile = 0;
        Files.APath = 0;
        Files.ZFile = 0;
        Files.ZPath = 0;
        Status.InPar(1:5)=[0 0 0 0 0];
        Status.Analysis = 0;
        cla(handles.Axes_Image)
        cla(handles.Axes_Intensity)
        xlabel(handles.Axes_Intensity, 'pixel  [-]', 'Color', 'w')
        ylabel(handles.Axes_Intensity, 'intensity  [-]', 'Color', 'w')
        set(handles.Button_FileChange, 'Enable', 'off')
        set(handles.Button_TBSet, 'Enable', 'off')
        set(handles.Button_LRSet, 'Enable', 'off')
        set(handles.Button_Initialize, 'Enable', 'off')
        set(handles.Text_FileA, 'String', '- no data selected -', 'ForegroundColor', [1 0 0])
        set(handles.Text_FileZ, 'String', '- no data selected -', 'ForegroundColor', [1 0 0])
        set(handles.Text_TBound, 'String', '---', 'ForegroundColor', [1 0 0])
        set(handles.Text_BBound, 'String', '---', 'ForegroundColor', [1 0 0])
        set(handles.Text_LBound, 'String', '---', 'ForegroundColor', [1 0 0])
        set(handles.Text_RBound, 'String', '---', 'ForegroundColor', [1 0 0])
        set(handles.Text_CurrentFile, 'String', '- - -', 'ForegroundColor', [1 1 0])
        Files.CurrentDir = Files.DefaultDir;
        cd (Files.CurrentDir);

    end
    if (Status.InPar - 1) == zeros
        set(handles.Button_Start, 'Enable', 'on')
    else
        set(handles.Button_Start, 'Enable', 'off')
    end
    if Status.InPar(5) == 0
        set(handles.Text_Status, 'String', '-  N O T   I N I T I A L I Z E D  -', 'ForegroundColor', [1 0 0])
    else
        set(handles.Text_Status, 'String', '-  R E A D Y  -', 'ForegroundColor', [0 1 0])
    end
    Status.Loop = 0;
    figure(handles.Main)
    uiwait(handles.Main)
end

close all
cd (Files.OldDir);

function varargout = StrainMeasurement_OutputFcn(hObject, eventdata, handles)
%varargout{1} = handles.output;

function Button_FileALoad_Callback(hObject, eventdata, handles)
    global Status Parameters Files
    Status.Loop = 1;
    [Files.AFile, Files.APath] = uigetfile({'*.bmp;*.tif;*.jpg;*.TIF;*.BMP;*.JPG','image files (*.bmp,*.tif,*.jpg)';'*.*',  'all files (*.*)'}, 'SELECT FIRST IMAGE', Files.CurrentDir);
    if Files.APath ~= 0
        set(handles.Text_FileA, 'String', Files.AFile, 'ForegroundColor', [0 1 0])
        Files.CurrentDir = Files.APath;
        Status.InPar(1) = 1;
        [Files.ZFile, Files.ZPath] = uigetfile({'*.bmp;*.tif;*.jpg;*.TIF;*.BMP;*.JPG','image files (*.bmp,*.tif,*.jpg)';'*.*',  'all files (*.*)'}, 'SELECT LAST FILE', Files.CurrentDir);
        if Files.ZPath ~= 0
            set(handles.Text_FileZ, 'String', Files.ZFile, 'ForegroundColor', [0 1 0])
            Files.CurrentDir = Files.ZPath;
            Status.InPar(2) = 1;
        end
    end
    if Status.InPar(1:2)-1 == 0
        CreateFilelist(handles)
    end
    uiresume

function Button_FileChange_Callback(hObject, eventdata, handles)
    global Status Parameters Files
    Status.Loop = 1;
    set(handles.Panel_ChangeFile, 'Visible', 'on')
    set(handles.Listbox_File, 'String', Files.Filelist)

function Button_TBSet_Callback(hObject, eventdata, handles)
    global Status Parameters Files
    Status.Loop = 1;
    set(handles.Text_Status, 'String', sprintf('Click TOP and BOTTOM of the integration range \nin the UPPER IMAGE.'), 'ForegroundColor', [1 1 0])
    axes(handles.Axes_Image)
    Input = ginput(2);
    Input = sort(Input, 1, 'ascend');
    Parameters.TB = round(Input(1, 2));
    Parameters.BB = round(Input(2, 2));
    Status.InPar(3) = 1;
    DisplayImage(handles, Files.ActualFile)
    uiresume

function Button_LRSet_Callback(hObject, eventdata, handles)
    global Status Parameters Files
    Status.Loop = 1;
    set(handles.Text_Status, 'String', sprintf('Click LEFT BORDER and RIGHT BORDER \nof the integration range.'), 'ForegroundColor', [1 1 0])
    axes(handles.Axes_Intensity)
    Input = ginput(2);
    Input = sort(Input,1, 'ascend');
    Parameters.LB = round(Input(1, 1));
    Parameters.RB = round(Input(2, 1));
    Status.InPar(4) = 1;
    DisplayImage(handles, Files.ActualFile)
    uiresume

function Checkbox_IIP_Callback(hObject, eventdata, handles)

function Button_Start_Callback(hObject, eventdata, handles)
    global Status Parameters Files
    set(handles.Button_FileALoad, 'Enable', 'off')
    set(handles.Button_FileChange, 'Enable', 'off')
    set(handles.Button_TBSet, 'Enable', 'off')
    set(handles.Button_LRSet, 'Enable', 'off')
    set(handles.Button_Start, 'Enable', 'off')
    set(handles.Button_Initialize, 'Enable', 'off')
    Status.Loop = 1;
    FileCount = 0;
    AnalysisLoop = 1;
    Status.Analysis = 1;
    Parameters.VAR = get(handles.Checkbox_VAR, 'Value');
    Parameters.ICB = get(handles.Checkbox_ICB, 'Value');
    xtest = [1 : Parameters.ImageSize(1)];

    while AnalysisLoop == 1
        FileCount = FileCount + 1;
        DisplayImage(handles, Files.Filelist(FileCount, :))
        xdata = Parameters.IntegratedIntensity(:,1);
        ydata = Parameters.IntegratedIntensity(:,2);
        axes(handles.Axes_Intensity)
        hold on
        % data fit with two error functions, fixed at y(x0) = y(x1) = 0
        if FileCount == 1
            FitGuess = [max(ydata)/2, 0.1, Parameters.InitialGuess(1), 0.1, Parameters.InitialGuess(2)];
        end
        if Parameters.ICB == 1
            %data correction
            MaxLeft = 2/FitGuess(2)+FitGuess(3);
            MaxRight = FitGuess(5) - 2/FitGuess(4);
            Range = [0,0];
            i = 0;
            while Range(1) == 0
                i = i + 1;
                if xdata(i) >= MaxLeft
                	Range(1) = i;
                end
            end
            while Range(2) == 0
            	i = i + 1;
                if xdata(i) >= MaxRight
                	Range(2) = i-1;
                end
            end
            ydata(Range(1):Range(2)) = Parameters.UpperValue;
            FitGuess(1) = Parameters.UpperValue/2;
            plot(xdata(Range(1):Range(2)), ydata(Range(1):Range(2)), 'g-')
        end
        Exitflag = 0;
        LineHandle = 0;
        j = 0;
        while Exitflag == 0
            j = j + 1;
            [FitData, Resnom, Residual, Exitflag] = lsqcurvefit(@(FitData, xdata)  FitData(1)*(erf(FitData(2)*(xdata-FitData(3)))-erf(FitData(4)*(xdata-FitData(5)))), FitGuess, xdata, ydata);
            ytest = FitData(1)*(erf(FitData(2)*(xtest-FitData(3)))-erf(FitData(4)*(xtest-FitData(5))));
            if LineHandle ~= 0
                delete(LineHandle)
            end
            LineHandle = plot (xtest, ytest, 'y-');
            FitGuess = FitData;
            if j == 10
                Exitflag = 1;
            end
        end
        clear xdata ydata ytest
        if FileCount == 1
            RelAnalysisRange = [Parameters.LB - FitData(3), Parameters.RB - FitData(5)];
        end
        if Parameters.VAR == 1
            Parameters.LB = round(FitData(3) + RelAnalysisRange(1));
            Parameters.RB = round(FitData(5) + RelAnalysisRange(2));
            if Parameters.LB < 1
                Parameters.LB = 1;
            end
            if Parameters.RB > Parameters.ImageSize(1)
                Parameters.RB = Parameters.ImageSize(1);
            end
        end
        MarkerPosition(FileCount, 1:2) = [FitData(3), FitData(5)];
        if FileCount == (size(Files.Filelist, 1))
            AnalysisLoop = 0;
        end
        Completed = round(FileCount/size(Files.Filelist, 1) * 100);
        set(handles.Text_Status, 'String', ['A N A L Y S I S   R U N N I N G   -   ', num2str(Completed), '%   C O M P L E T E D'], 'ForegroundColor', [1 0 0])

    end
    Status.Analysis = 0;
    ImgNo(:,1) = [1:FileCount];
    MarkerDistance(:,1) = MarkerPosition(:,2) - MarkerPosition(:,1);
    Strain(:,1) = MarkerDistance(:,1)/MarkerDistance(1,1)-1;
    OutData(:, 1) = ImgNo(:,1);
    OutData(:, 2) = MarkerPosition(:,1);
    OutData(:, 3) = MarkerPosition(:,2);
    OutData(:, 4) = MarkerDistance(:,1);
    OutData(:, 5) = Strain(:,1);
    save([Files.Path, 'StrainMeasurement.dat'], 'OutData', '-ascii', '-tabs')
    save([Files.Path, 'StrainX.dat'], 'Strain', '-ascii', '-tabs')
    PlotHandle = figure;
    plot(ImgNo(:,1), Strain(:,1), 'bo')
    ylabel('strain  [-]')
    xlabel('image no.  [-]')
    set(gca, 'Box', 'on')
    saveas(PlotHandle, [Files.Path, '\StrainX.fig'])
    figure(handles.Main)
    set(handles.Button_FileALoad, 'Enable', 'on')
    set(handles.Button_FileChange, 'Enable', 'on')
    set(handles.Button_TBSet, 'Enable', 'on')
    set(handles.Button_LRSet, 'Enable', 'on')
    set(handles.Button_Start, 'Enable', 'on')
    set(handles.Button_Initialize, 'Enable', 'on')
    set(handles.Text_Status, 'String', '-  A N A L Y S I S   C O M P L E T E D  -', 'ForegroundColor', [0 1 0])
    clear xtest MarkerPosition ImgNo MarkerDistance Strain OutData
    DisplayImage(handles, Files.Filelist(1,:))

    uiresume

function Button_Initialize_Callback(hObject, eventdata, handles)
    global Status Parameters Files
    Status.Loop = 1;
    set(handles.Text_Status, 'String', sprintf('Click CENTER of the LEFT SLOPE \nand CENTER of the RIGHT SLOPE \nin the LOWER GRAPH.'), 'ForegroundColor', [1 1 0])
    axes(handles.Axes_Intensity)
    Input = ginput(2);
    Input = sort(Input,1, 'ascend');
    Parameters.InitialGuess(1:2) = Input(1:2, 1);
    Status.InPar(5) = 1;
    uiresume

function CreateFilelist(handles)
    global Status Parameters Files
    String = get(handles.Text_Status, 'String');
    set(handles.Text_Status, 'String', 'WAIT.....')
    Files.Filelist = [];
    FileLength(1) = size(Files.AFile,2);
    FileLength(2) = size(Files.ZFile,2);
    FileType(1,:) = Files.AFile((FileLength(1)-3):FileLength(1));
    FileType(2,:) = Files.ZFile((FileLength(2)-3):FileLength(2));
    FileStart = str2num(Files.AFile((FileLength(1)-6):(FileLength(1)-4)));
    FileEnd = str2num(Files.ZFile((FileLength(2)-6):(FileLength(2)-4)));
    FileLabel(1,:) = Files.AFile(1:(FileLength(1)-7));
    FileLabel(2,:) = Files.ZFile(1:(FileLength(2)-7));

    if FileType(1,:) ~= FileType(2,:)
        fprintf('CRITICAL FILENAME ERROR - SCRIPT TERMINATED \n');
        uiresume
    end
    if FileLabel(1,:) ~= FileLabel(2,:)
        fprintf('CRITICAL FILENAME ERROR - SCRIPT TERMINATED \n');
        uiresume
    end
    if Files.APath ~= Files.ZPath
        fprintf('CRITICAL FILENAME ERROR - SCRIPT TERMINATED \n');
        uiresume
    end

    j = 0;
    for i = FileStart:FileEnd
        j = j + 1;
        FileNo = num2str(i + 1000);
        Filename(j,:) = [FileLabel(1,:), FileNo(:,2:4), FileType(1,:)];
    end
    Files.Filelist = Filename;
    Parameters.FileType = FileType(1,:);
    Parameters.FileLabel = FileLabel(1,:);
    Files.Path = Files.APath;
    clear Filename FileLength FileType FileLabel

    Files.ActualFile = Files.Filelist(1,:);
    DisplayImage(handles, Files.ActualFile)

    set(handles.Button_FileChange, 'Enable', 'on')
    set(handles.Button_TBSet, 'Enable', 'on')
    set(handles.Button_LRSet, 'Enable', 'on')
    Status.InPar(3:5) = [0 0 0];

    set(handles.Button_Initialize, 'Enable', 'on')

function DisplayImage(handles, Filename)
    global Status Parameters Files

    set(handles.Text_CurrentFile, 'String', Filename)
    TempData = imread([Files.Path, Filename]);
    ImageData = mean(double(TempData), 3);
    Parameters.ImageSize = [size(ImageData, 2) size(ImageData, 1)];
    clear TempData

    axes(handles.Axes_Image)
    cla(handles.Axes_Image)
    axis ij
    hold on
    imagesc(ImageData);
    if Status.InPar(3) == 0
        Parameters.TB = 1;
        Parameters.BB = size(ImageData, 1);
    else
        plot([1 Parameters.ImageSize(1)], [Parameters.TB Parameters.TB], 'w-')
        plot([1 Parameters.ImageSize(1)], [Parameters.BB Parameters.BB], 'w-')
    end
    if Status.InPar(4) == 0
        Parameters.LB = 1;
        Parameters.RB = size(ImageData, 2);
    else
        plot([Parameters.LB Parameters.LB], [1 Parameters.ImageSize(2)], 'w-')
        plot([Parameters.RB Parameters.RB], [1 Parameters.ImageSize(2)], 'w-')
    end
    set(handles.Axes_Image, 'Box', 'on', 'XColor', 'w', 'YColor', 'w', 'XTickLabelMode', 'manual', 'XTickLabel', [], 'XTickMode', 'manual', 'XTick', [],'YTickLabelMode', 'manual', 'YTickLabel', [], 'YTickMode', 'manual', 'YTick', [], 'XLim', [1 Parameters.ImageSize(1)], 'YLim', [1 Parameters.ImageSize(2)])
    set(handles.Text_TBound, 'String', Parameters.TB, 'ForegroundColor', [0 1 0])
    set(handles.Text_BBound, 'String', Parameters.BB, 'ForegroundColor', [0 1 0])
    set(handles.Text_LBound, 'String', Parameters.LB, 'ForegroundColor', [0 1 0])
    set(handles.Text_RBound, 'String', Parameters.RB, 'ForegroundColor', [0 1 0])
    Status.InPar(3:4) = [1 1];
    Status.InPar(5) = 0;

    IntegratedIntensity = [];
    Parameters.IntegratedIntensity = [];
    Parameters.UpperValue = 0;
    IntegratedIntensity(1:(Parameters.RB-Parameters.LB+1), 1:2) = zeros;
    IntegratedIntensity(1:(Parameters.RB-Parameters.LB+1), 1) = [Parameters.LB:Parameters.RB];
    for i = Parameters.TB:Parameters.BB
        TempData(:,1) = ImageData(i, Parameters.LB:Parameters.RB);
        IntegratedIntensity(:,2) = IntegratedIntensity(:,2) + TempData;
        if Status.Analysis == 1 & Parameters.ICB == 1
        	TempData = sort(TempData, 'descend');
            Parameters.UpperValue = Parameters.UpperValue + mean(TempData(1:round(0.1*size(TempData, 1)), 1));
        end
    end
    axes(handles.Axes_Intensity)
    hold off
    plot(IntegratedIntensity(:,1), IntegratedIntensity(:,2), 'g.')
    set(handles.Axes_Intensity, 'Color', [0 0 0], 'Box', 'on', 'XColor', 'w', 'YColor', 'w', 'XLim', [1, Parameters.ImageSize(1)])
    xlabel(handles.Axes_Intensity, 'pixel  [-]', 'Color', 'w')
    ylabel(handles.Axes_Intensity, 'intensity  [-]', 'Color', 'w')

    Parameters.IntegratedIntensity = IntegratedIntensity;
    clear IntegratedIntensity TempData UpperValueTemp

function Listbox_File_Callback(hObject, eventdata, handles)
    global Status Parameters Files
    FileNo = get(handles.Listbox_File, 'Value');
    Files.ActualFile = Files.Filelist(FileNo, :);
    set(handles.Panel_ChangeFile, 'Visible', 'off')
    DisplayImage(handles, Files.ActualFile)

    uiresume

function Button_Reset_Callback(hObject, eventdata, handles)
    global Status Parameters Files
    Status.Loop = 2;
    Files.Filelist = [];

    uiresume


text:
obj    = videoinput('winvideo', 1, 'RGB24_320x240');

vidRes = get(obj, 'VideoResolution');
nBands = get(obj, 'NumberOfBands');
hImage = image( zeros(vidRes(2), vidRes(1), nBands) );
preview(obj, hImage);

repo: OpenTLD


text:
classdef PluginOptimizationVariable < AbstractOptimizationVariable
    %LvdPluginOptimizationVariable Summary of this class goes here
    %   Detailed explanation goes here

    properties
        varObj(1,1) LvdPluginOptimVarWrapper = LvdPluginOptimVarWrapper();

        lb(1,1) double = -1;
        ub(1,1) double = 1;

        useTf(1,1) = false;
    end

    methods
        function obj = PluginOptimizationVariable(varObj)
            arguments
                varObj(1,1) LvdPluginOptimVarWrapper
            end

            obj.varObj = varObj;
            obj.varObj.optVar = obj;

            obj.id = rand();
        end

        function x = getXsForVariable(obj)
            x = [];

            if(obj.useTf)
                x = obj.varObj.value;
            end
        end

        function [lb, ub] = getBndsForVariable(obj)
            lb = obj.lb(obj.useTf);
            ub = obj.ub(obj.useTf);
        end

        function [lb, ub] = getAllBndsForVariable(obj)
            lb = obj.lb;
            ub = obj.ub;
        end

        function setBndsForVariable(obj, lb, ub)
            obj.lb = lb;
            obj.ub = ub;
        end

        function useTf = getUseTfForVariable(obj)
            useTf = obj.useTf;
        end

        function setUseTfForVariable(obj, useTf)
            obj.useTf = useTf;
        end

        function updateObjWithVarValue(obj, x)
            obj.varObj.value = x;
        end

        function nameStrs = getStrNamesOfVars(obj, evtNum, varLocType)
            nameStrs = {sprintf('Plugin Variable "%s"', obj.varObj.name)};
        end
    end
end


text:
function F = full_sparse(I,J,V,S1,S2)
  % FULL_SPARSE Faster than calling:
  %
  % F = full_sparse(I,J,V,S1,S2)
  %
  % F = full(sparse(I,J,V,S1,S2))
  %
  % See also: sparse, accumarray
  if nargin>3
    F = accumarray([I(:) J(:)],V(:),[S1 S2]);
  else
    F = accumarray([I(:) J(:)],V(:));
  end
end

author: alecjacobson

repo: gptoolbox


text:
function nearest = find_closest3 ( m, nr, r, ns, s )

%*****************************************************************************80
%
%% FIND_CLOSEST3 finds the nearest R point to each S point.
%
%  Discussion:
%
%    We are given R, a set of NR points in M dimensions.
%
%    We are given S, a set of NS points in M dimensions.
%
%    For each S(I) in S, we seek the index J of the point R(J)
%    which is nearest to S(I) over all points in R.
%
%    Modified in accordance with suggestions by Gene Cliff, 08 July 2010.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 July 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer NR, the number of cell generators.
%
%    Input, real R(M,NR), the cell generators.
%
%    Input, integer NS, the number of sample points.
%
%    Input, real S(M,NS), the points to be checked.
%
%    Output, integer NEAREST(NS), the index of the cell generator nearest
%    to the sample point.
%
  ones_k = ones ( 1, nr );
  nearest = NaN ( 1, ns );

  for i = 1 : ns
    d1 = r - s(:,i) * ones_k;
    d2 = sum ( d1 .* d1 );
    [ min_val, min_loc ] = min ( d2 );
    nearest(i) = min_loc;
  end

  return
end

repo: jburkardt-m


text:
function plot_mix_gaussian( u,sig,prob,X )
%
% plot_mix_gaussian - plot the samples and the estimation.
%                     for 1D distribution, plot the normalized histogram and the 1D distribution function
%                     for 2D distribution, plot the samples and the contour of FWHM
%                     for mD distribution (m>2) - do not plot nothing.
%
% format:   plot_mix_gaussian( u,sig,prob,X )
%
% input:    u       - mean of each gaussian in the distribution. 1xM vector or 2xM matrix.
%                     each gaussian mean is stored in a separate column.
%           sig     - for 1D distribution -> standard deviation of each gaussian in the distribution
%                     each gaussian mean is stored in a separate column -> a 1xM vector
%                     for 2D distribution -> covariance matrix for each gaussian in the distribution
%                     2x2xM matrix, the 3rd dimension is the gaussians index,
%                     the 1st and 2nd dimensions are the covariance matrix
%           prob    - probability of each gaussian in the distribution to create the current sample.
%                     this is a 1xM vector
%           X       - the samples, 1xN vector or 2xN matrix, depends on the dimension of the
%                     distribution (i.e. 1D or 2D)
%
% output:   to the graphic current axis.
%
%
%

% check input
if (nargin<3)
    error( 'plot_mix_gaussian - insufficient input parameters' );
end
if ~exist( 'X' )
    X = [];
end

% constants
nbins = 200;

% check the size of the input to determin if it's 1D or 2D distribution
if (size(u,1)==2) & (size(u,2)==size(prob,2))
    plot_mix_gaussian_2d( u,sig,prob,X );
else
    plot_mix_gaussian_1d( u,sig,prob,X,nbins );
end
drawnow;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                              Inner function implementation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function plot_mix_gaussian_1d( u,sig,prob,X,nbins )
% plot normalized histogram and normalized distribution on top of it

% constant
points = 1000;

if ~isempty( X )
    [n,x]   = hist( X,nbins );              % calc the histogram
    dx      = x(2)-x(1);                    % calc a single bin width
    n       = n / sum( n*dx );              % normalize histogram to have area of 1
    bar( x,n,'hist' );                      % plot normalized histogram
    xlim( [x(1)-dx/2,x(end)+dx/2] );        % make sure that the axis is squeezed to it's limits
    g       = gca;
    c       = [1 0 0];                      % choose the red color
    x_c     = linspace( x(1)-dx/2,x(end)+dx/2,points );
else
    g       = gca;
    X       = xlim( g );                    % get the axis limits
    c       = [0 1 0];                      % choose the green color
    x_c     = linspace( X(1),X(2),points );
end

% plot the distribution
for m = 1:length(prob)
    y   = prob(m) / sqrt( 2*pi*sig(m)^2 ) * exp( -((x_c-u(m)).^2)/(2*sig(m)^2) );
    line( x_c,y,'color',c,'parent',g,'linewidth',2 );
end
shg;
drawnow;

% ----------------------------------------------------------------------------------------

function plot_mix_gaussian_2d( u,covar,prob,X )
% plot 2D samples and distribution data on top of it

if ~isempty( X )
    if size(X,2)<size(X,1),
        X = X.';
    end
    plot( X(1,:),X(2,:),'.k' );             % plot samples
    xlim( [min(X(1,:)) max(X(1,:))] );      % squeeze axis limits
    ylim( [min(X(2,:)) max(X(2,:))] );       % squeeze axis limits
    g       = gca;
    c       = [1 0 0];                      % choose the red color
else
    g       = gca;
    c       = [0 1 0];                      % choose the green color
end

% plot each gaussian's FWHM and mean (center)
a           = linspace(0,2*pi,361); % degree vector
unit_circle = [cos(a);sin(a)];      % a unit circle for the case of sig_x=sig_y=1
for m = 1:length(prob)
    U       = u(:,m);               % the mean of the mth gaussian
    [V,D]   = eig(covar(:,:,m));    % eigen value and vectors for the covariance of the mth gaussian
    A       = V*sqrt(D);            % factorization: covar = A*A' = positive matrix
    outline = A * unit_circle * sqrt(2*log(2)) + U*ones(1,length(unit_circle)); % the contour
    line( U(1),U(2),'color',c,'parent',g,'marker','o' );
    line( outline(1,:),outline(2,:),'color',c,'parent',g,'linewidth',2 );
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function kern = pathKernExpandParam(kern, params)

% PATHKERNEXPANDPARAM Create kernel structure from PATH kernel's parameters.
% FORMAT
% DESC returns a path kernel structure filled with the
% parameters in the given vector. This is used as a helper function to
% enable parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG kern : the kernel structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% kernel structure.
% RETURN kern : kernel structure with the given parameters in the
% relevant locations.
%
% SEEALSO : pathKernParamInit, pathKernExtractParam, kernExpandParam
%
% COPYRIGHT : Andrea Baisero, Carl Henrik Ek, 2013

% SHEFFIELDML

kern.cd=(params(1))^2;
kern.chv=params(2);
kern.gkern=kernExpandParam(kern.gkern,params(3:end));

kern = pathKernUpdateWMat(kern,[],true);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% The COBRAToolbox: testcomputeFluxSplits.m
%
% Purpose:
%     - tests the basic functionality of computeFluxSplits
%       Tests 1 optional input with 3 possibilities: coeffSign=[] or 0 or 1
%       returns 1 if all tests were completed succesfully, 0 if not
%
% Authors:
%     - Original file: Diana El Assal 03/08/2017
%

% % save the current path
currentDir = pwd;

% % initialize the test
fileDir = fileparts(which('testComputeFluxSplits'));
cd(fileDir);
global CBTDIR

% define the solver packages to be used to run this test
solverPkgs = {'gurobi', 'tomlab_cplex', 'glpk'};

% load the model
model = getDistributedModel('Recon2.0model.mat');

% define the metabolites of interest
mets = {'atp[c]','atp[m]', 'atp[e]'};

% obtain a flux vector using e.g. FBA or uniform sampling
model = changeObjective(model, 'DM_atp_c_');
FBAsolution = optimizeCbModel(model, 'max');
V = FBAsolution.x;

for k = 1:length(solverPkgs);

    % change the COBRA solver (LP)
    solverOK = changeCobraSolver(solverPkgs{k}, 'LP', 0);

    if solverOK == 1
        fprintf('   Testing compute flux splits using %s ... ', solverPkgs{k});

        % check the function without optional input
        fprintf('\n>> without optional input\n');
        [P,C,vP,vC,s] = computeFluxSplits(model,mets,V,[]);
        assert(isequal(sum([model.S(ismember(model.mets,mets),:);...
            sparse(1,size(model.S,2))],1)',s));

        % check the optional input of coeffSign = 0
        fprintf('\n>> with coeffSign = false (stoichiometric coefficient)\n');
        [P,C,vP,vC,s] = computeFluxSplits(model,mets,V,0);
        assert(isequal(sum([model.S(ismember(model.mets,mets),:);...
            sparse(1,size(model.S,2))],1)',s));

        % check the optional input of coeffSign = 1
        fprintf('\n>> with coeffSign = true (sign of stoichiometric coefficient)\n');
        [P,C,vP,vC,s] = computeFluxSplits(model,mets,V,1);
        sInd = find(s);
        for i = length(sInd);
            j = sInd(i);
            assert(s(j) == 1 || s(j) == -1);
        end

        % output a success message
        fprintf('Done.\n');
    end
end

clear i j

% change the directory
cd(currentDir)