text:
function lineIdx = ...
  mrGeodesic(grayNodes,grayEdges,dimdist,radius,spacing,crit,seed1,seed2)
%
% lineIdx =
%    mrGeodesic(grayNodes,grayEdges,dimdist, ...
%        radius,spacing,crit,seed1,[seed2])
%
% AUTHOR:  Heeger, Wandell
% DATE:
% PURPOSE:
%   Find the gray node indices along a geodesic between seed1 and
% seed2.  The gray node indices are separated by spacing.
%
% grayNodes, grayEdges, dimdist:  These are the usual parameters
%      from the unfolding analysis.
%
% radius:  If one seed, choose the second near this radius
% spacing: Spacing of sample points along the geodesics
% crit:    Criterion level that determines how close does the
%           point have to be to satisfying the "on the line"
% seed1:   Starting point for one end of the geodesic
% [seed2:] Second point on the geodesic.
%
% TODO:
%   Order of arguments could be better.  I am not sure we handle
% to two seed case properly.  BW
%

errordlg('This function is obsolete.  Use the method described in mrManDist.  Then implement that method here.');

noVol = -1;

% Find the distances from seed1 and a second seed that is as far
% away as possible
%
dist1 = mrManDist(grayNodes, grayEdges, seed1, dimdist, noVol, radius + 1);

if (~exist('seed2'))
  [dist12 idx] = max(dist1);
  seed2 = idx;
else
  dist12 = dist1(seed2);
  if (radius < dist12)
    disp('Increasing radius')
    radius = dist12 + 1;
  end
end

% The points along a geodesic between 1 and 2 should have a
% distance from 1 and 2 that sums to the distance between 1 and 2
%
dist2 = ...
    mrManDist(grayNodes, grayEdges, seed2, dimdist, noVol,radius + 1);
between12 = find( abs(dist1 + dist2 - dist12) < crit);

% This finds the points spaced by an amount sep along the
% geodesic
%

sep = [spacing:spacing:dist12];
nSep = length(sep);
lineIdx = ones(nSep,1)*(-1);

for count = 1:nSep;
  [v idx] = min(abs(dist1(between12) - sep(count)));
  lineIdx(count) = between12(idx);
end

return;

% DEBUGGING:
% cd /usr/local/matlab/toolbox/stanford/mri/unfold/Example
% [grayNodes grayEdges] = readGrayGraph('Graph.gray');
% seed1 = 100
% radius = 30
% dimdist = [1.0667 1.0667 1.0000];
% spacing = 1;
% lineIdx = mrGeodesic(grayNodes,grayEdges,dimdist,radius,spacing,crit,seed1,seed2);
%
% grayNodes(1:3,lineIdx(1:nSep,1))
% To view this, we need to plot these points on top of the volume
% anatomy somehow.

author: vistalab

repo: vistasoft


text:
function setcustompwd(option_pwd)

% SETCUSTOMPWD is a helper function that updates the present working directory.

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

% these are for remembering the path on subsequent calls with the same input arguments
persistent previous_argin previous_pwd

if isequal(previous_argin, option_pwd) && isequal(previous_pwd, pwd)
  % no reason to change the pwd again
  return
end

if ~isempty(option_pwd)
  try
    cd(option_pwd);
  catch
    % the "catch me" syntax is broken on MATLAB74, this fixes it
    cd_error = lasterror;
    % don't throw an error, just give a warning (and hope for the best...)
    warning(cd_error.message);
  end
end

% remember the current settings for the next call
previous_argin = option_pwd;
previous_pwd   = pwd;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef testSymmetryForIAniTensorInVoigt < handle

    properties (Access = protected)
        ChVoigtSym
        Ch
        ChVoigt
    end

    properties (Access = public)
        tol = 1e-12;
    end

    methods (Access = public)

        function obj = testSymmetryForIAniTensorInVoigt()
            obj.computeFourthOrderTensor();
            obj.computeVoigtRepresentation();
            obj.computeSymetricVoigthTensor();
        end

        function error = computeError(obj)
            c    = obj.ChVoigt.getValue();
            cSym = obj.ChVoigtSym;
            error = norm(c(:) - cSym(:));
        end

    end

    methods (Access = private)
        function computeFourthOrderTensor(obj)
            obj.Ch = Stiffness3DTensor();
            obj.Ch.createRandomTensor();
        end

        function computeVoigtRepresentation(obj)
            obj.ChVoigt = Tensor2VoigtConverter.convert(obj.Ch);
        end

        function computeSymetricVoigthTensor(obj)
            t = obj.ChVoigt.getValue();
            obj.ChVoigtSym = 0.5*(t + t');
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [ point, seed ] = p05_sample ( m, n, seed )

%*****************************************************************************80
%
%% P05_SAMPLE samples points from the region in problem 05.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real POINT(M,N), the coordinates
%    of the points.
%
%    Output, integer SEED, a seed for the random number generator.
%
  center1 = [  0.0, 0.0 ];
  center2 = [ -0.4, 0.0 ];
  r1 =   1.00;
  r2 =   0.55;

  x1 = center1(1) - r1;
  x2 = center1(1) + r1;

  y1 = center1(2);
  y2 = center1(2) + r1;

  have = 0;
%
%  We are going to generate batches of sample points.
%
  sample_num = min ( 1000, n );

  while ( 1 )
%
%  Generate a batch of points in [0,1]x[0,1].
%
    [ sample, seed ] = r8mat_uniform_01 ( m, sample_num, seed );
%
%  Remap the points to the box [X1,X2] x [Y1,Y2].
%
    sample(1,1:sample_num) = x1 + sample(1,1:sample_num) * ( x2 - x1 );
    sample(2,1:sample_num) = y1 + sample(2,1:sample_num) * ( y2 - y1 );
%
%  Accept those points which are in the big circle and not in the
%  small circle.
%
    for j = 1 : sample_num

      if (                                                          ...
                   ( sample(1,j) - center1(1) ).^2                  ...
                 + ( sample(2,j) - center1(2) ).^2 <= r1 * r1 &     ...
        r2 * r2 <= ( sample(1,j) - center2(1) ).^2                  ...
                 + ( sample(2,j) - center2(2) ).^2 )

        have = have + 1;
        point(1:m,have) = sample(1:m,j);

        if ( have == n )
          return
        end

      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
% UPDATE_ESS Update the Expected Sufficient Statistics of a hhmmF node.
% function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)

% Figure out the node numbers associated with each parent
% so we extract evidence from the right place
dom = fmarginal.domain; % Q(1) .. Q(d) F(d+1) F(d)
Qps = fmarginal.domain(1:end-2);
Q = Qps(end);
Qps = Qps(1:end-1);

Qsz = CPD.Qsizes(CPD.Q);
Qpsz = prod(CPD.Qsizes(CPD.Qps)); % may be 1

% We assume the F node are always hidden, but allow some of the Q nodes
% to be observed. We do case analysis for speed.
%We only extract prob from fmarginal.T when F(d+1)=2 i.e., model below has finished.
% wrong -> % We sum over the possibilities that F(d+1) = 1 or 2

obs_self = ~hidden_bitv(Q);
if obs_self
  self_val = evidence{Q};
end

if isempty(Qps) % independent of parent context
  counts = zeros(Qsz, 2);
  %fmarginal.T(Q(d), F(d+1), F(d))
  if obs_self
    marg = myreshape(fmarginal.T, [1 2 2]);
    counts(self_val,:) = marg(1,2,:);
    %counts(self_val,:) = marg(1,1,:) + marg(1,2,:);
  else
    marg = myreshape(fmarginal.T, [Qsz 2 2]);
    counts = squeeze(marg(:,2,:));
    %counts = squeeze(marg(:,2,:)) + squeeze(marg(:,1,:));
  end
else
  counts = zeros(Qpsz, Qsz, 2);
  %fmarginal.T(Q(1:d-1), Q(d), F(d+1), F(d))
  obs_Qps = ~any(hidden_bitv(Qps));  % we assume that all or none of the Q  parents are observed
  if obs_Qps
    Qps_val = subv2ind(Qpsz, cat(1, evidence{Qps}));
  end
  if obs_self & obs_Qps
    marg = myreshape(fmarginal.T, [1 1 2 2]);
    counts(Qps_val, self_val, :) = squeeze(marg(1,1,2,:));
    %counts(Qps_val, self_val, :) = squeeze(marg(1,1,2,:)) + squeeze(marg(1,1,1,:));
  elseif ~obs_self & obs_Qps
    marg = myreshape(fmarginal.T, [1 Qsz 2 2]);
    counts(Qps_val, :, :) = squeeze(marg(1,:,2,:));
    %counts(Qps_val, :, :) = squeeze(marg(1,:,2,:)) + squeeze(marg(1,:,1,:));
  elseif obs_self & ~obs_Qps
    error('not yet implemented')
  else
    marg = myreshape(fmarginal.T, [Qpsz Qsz 2 2]);
    counts(:, :, :) = squeeze(marg(:,:,2,:));
    %counts(:, :, :) = squeeze(marg(:,:,2,:)) + squeeze(marg(:,:,1,:));
  end
end

CPD.sub_CPD_term = update_ess_simple(CPD.sub_CPD_term, counts);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function barycentric_interp_1d_test02 ( prob, nd )

%*****************************************************************************80
%
%% BARYCENTRIC_INTERP_1D_TEST02 tests LAGCHEBY2_INTERP_1D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROB, the problem index.
%
%    Input, integer ND, the number of data points to use.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BARYCENTRIC_INTERP_1D_TEST02:\n' );
  fprintf ( 1, '  Interpolate data from TEST_INTERP_1D problem #%d.\n', prob );
  fprintf ( 1, '  Use Chebyshev Type 2 spacing for data points.\n' );
  fprintf ( 1, '  Number of data points = %d\n', nd );
%
%  Define the data.
%
  a =  0.0;
  b = +1.0;
  xd = r8vec_cheby2space ( nd, a, b );
  yd = p00_f ( prob, nd, xd );

  if ( nd < 10 )
    r8vec2_print ( nd, xd, yd, '  Data array:' );
  end
%
%  #1:  Does the interpolant match the function at the interpolation points?
%
  ni = nd;
  xi = xd;
  yi = lagcheby2_interp_1d ( nd, xd, yd, ni, xi );

  int_error = norm ( yi - yd ) / ni;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  L2 interpolation error averaged per interpolant node = %g\n', int_error );
%
%  #2: Plot the data.
%
  ni = 500;
  xi = r8vec_linspace ( ni, a, b );
  yi = lagcheby2_interp_1d ( nd, xd, yd, ni, xi );

  clf
  hold on
  plot ( xd, yd, 'b.', 'MarkerSize', 30 );
  plot ( xi, yi, 'r-', 'LineWidth', 3 );
  grid on
  xlabel ( '<---X--->' );
  ylabel ( '<---Y--->' );
  title ( sprintf ( 'Barycentric Lagrange Chebyshev2 Interpolant for %d nodes', nd ) )

  filename = sprintf ( 'p%02d_lagcheby2_%02d.png', prob, nd );
  print ( '-dpng', filename );
  fprintf ( 1, '  Created plot file "%s".\n', filename );

  return
end


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: plots a B-spline
%
%==============================================================================

FAIRfigure(1,'color','w'); clf
tt = linspace(-3,11,1001);
p1 = plot(tt,spline1D(0,tt)); hold on; axis([-3,12,0,1])
set(p1,'linewidth',3,'color','k')
p2 = plot(tt,spline1D(2,tt),'--');
p3 = plot(tt,spline1D(7,tt),'--');
set([p2;p3],'linestyle','--','linewidth',1.5,'color','k'); axis off;
%==============================================================================


text:
function ex4bvp
%EX4BVP  Example 4 of the BVP tutorial.
%   This is a nerve impulse model considered in Example 7.1 of
%   R. Seydel, From Equilibrium to Chaos, Elsevier, 1988.  The
%   differential equations
%
%      y1' = 3*(y1 + y2 - 1/3*y1^3 + lambda)
%      y2' = -(y1 - 0.7 + 0.8*y2)/3
%
%   are to be solved subject to periodic boundary conditions
%
%      y1(0) = y1(T)
%      y2(0) = y2(T)
%
%   This is an example of non-separated boundary conditions, meaning
%   that some conditions involve values of the solution at both ends
%   of the interval.  Periodic boundary conditions are a common source
%   of such problems.  BVP4C is unusual in that it accepts problems
%   with non-separated boundary conditions.
%
%   The parameter lambda has the value -1.3 in the text cited.  The
%   period T is unknown, which is to say that the length of the
%   interval is not known. Such problems require some preparation.
%   By scaling the independent variable to tau = t/T the problem is
%   posed on the fixed interval [0,1]. When this is done, T becomes
%   an unknown parameter, but BVP4C provides for unknown parameters.

% Copyright 1999, The MathWorks, Inc.

%  Periodic functions evaluated in EX4INIT are used as a guess for the
%  solution on a crude mesh of 5 equally spaced points.  A guess of 2*pi
%  is provided for the unknown parameter T.
solinit = bvpinit(linspace(0,1,5),@ex4init,2*pi);

options = bvpset('stats','on');

sol = bvp4c(@ex4ode,@ex4bc,solinit,options);
T = sol.parameters;

%  The independent variable needs to be rescaled to its original value
%  before plotting the solution.  The initial guess is also plotted.
clf reset
plot(T*sol.x,sol.y(1,:),T*solinit.x,solinit.y(1,:),'o')
legend('solution found','initial guess');
axis([0 T -2.2 2]);
title('Nerve impulse model');
ylabel('solution y_1(t)');
xlabel(['Periodic solution with period ',num2str(sol.parameters,4)]);
shg

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

function v = ex4init(x)
%EX4INIT  Guess function for Example 4 of the BVP tutorial.
%   V = EX4INIT(X) returns a column vector V that is a guess for y(x).
v = [ sin(2*pi*x)
      cos(2*pi*x)];

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

function dydt = ex4ode(t,y,T);
%EX4ODE  ODE funcion for Example 4 of the BVP tutorial.
dydt = [ 3*T*(y(1) + y(2) - 1/3*(y(1)^3) - 1.3);
         (-1/3)*T*(y(1) - 0.7 + 0.8*y(2)) ];

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

function res = ex4bc(ya,yb,T)
%EX4BC  Boundary conditions for Example 4 of the BVP tutorial.
res = [ya(1) - yb(1)
       ya(2) - yb(2)
       T*(-1/3)*(ya(1) - 0.7 + 0.8*ya(2)) - 1];

author: Sable

repo: mcbench-benchmarks


text:
% Check whether a graph is regular, i.e. whether every node has the same degree.
%
% INPUTS: adjacency matrix, nxn
% OUTPUTS: Boolean, 0 or 1
%
% Note: Defined for unweighted graphs only.
% GB: last updated, Sep 23, 2012

function S=isRegular(adj)

S=false;

degs=sum(adj>0); % remove weights and sum columns

if degs == degs(1)*ones(size(degs)); S = true; end

author: aeolianine

repo: octave-networks-toolbox


text:
function model=convertToCobraV2(model)
% Generate a model structure with metabolite and reaction subfields instead of vectors.
%
% create a COBRA v2 model structure with model.met(m) for metabolites and
% model.rxn(n) for reactions. This keeps all the information about a
% metabolite or reaction in one place. The primary keys model.mets and
% model.rxns are kept as is, and these can be used to access groups of
% model.met or model.rxn quickly.
%
% Replaces (*) with [*] to denote compartment in metabolite abbreviations
%
%
%INPUT
% model.S                   Stoichiometric matrix
% model.rxns                Abbreviation for each reaction (primary key)
% model.rxnNames            Full name for each reaction
% model.mets                Abbreviation for each metabolite (primary key)
% model.metNames            Full name for each metabolite
% model.biomassRxnAbbr      abbreviation of biomass reaction
%
%OPTIONAL INPUT
% model.metFormulas         m x 1 cell array of strings with chemical
%                           formula for metabolite
% model.metCharges          m x 1 numeric charge for each metabolite
%                           species
%
%OUTPUT
% model.S                       Sparse Stoichiometric matrix +/- with
%                               adjustment of reactions involving C02
% model.SIntRxnBool             Boolean of internal reactions, i.e.
%                               non-mass balanced reactions
% model.met(m).abbreviation     metabolite abbreviation (primary key)
% model.met(m).officialName     metabolite name
% model.rxn(n).abbreviation     reaction abbreviation (primary key)
% model.rxn(n).officialName     reaction name
% model.rxn(n).equation         reaction equation
% model.rxn(n).directionality   qualitative directionality
% model.rxn(n).regulationStatus {('On'),'Off'} Off if lb=ub=0
%
%OPTIONAL OUTPUT
% model.met(m).formula          metabolite elemental formula
%
% Ronan M. T. Fleming

[nMet,nRxn]=size(model.S);

%make it sparse
model.S=sparse(model.S);

numChar=1;
[allMetCompartments,uniqueCompartments]=getCompartment(model.mets,numChar);

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

for m=1:nMet
    %find the rows with columns that are all zeros
    if nnz(model.S(m,:))==0
        fprintf('%s\n',['metabolite ' model.mets{m} ' without a reaction???']);
%         error(['metabolite ' model.mets{m} ' without a reaction.']);
    end

    %replace (*) with [*] to denote compartment
    model.mets{m}(end-2)='[';
    model.mets{m}(end)=']';
end

% start the conversion to a cobra v2 style model
fprintf('%s\n','Converting to cobra toolbox v2 style model.')
[nMet,nRxn]=size(model.S);

for n=1:nRxn
    model.rxn(n).abbreviation=model.rxns{n};
    model.rxn(n).officialName=model.rxnNames{n};
    %equation
    equation=printRxnFormula(model,model.rxns(n),0);
    model.rxn(n).equation=equation{1};
    %directionality
    if model.lb(n)<0 && model.ub(n)>0
        model.rxn(n).directionality='reversible';
    end
    if model.lb(n)<0 && model.ub(n)<=0
        model.rxn(n).directionality='reverse';
    end
    if model.lb(n)>=0 && model.ub(n)>0
        model.rxn(n).directionality='forward';
    end
    if model.lb(n)==0 && model.ub(n)==0
        model.rxn(n).directionality='off';
    end
end

%new metabolite structure
for m=1:nMet
    model.met(m).abbreviation=model.mets{m};
    model.met(m).officialName=model.metNames{m};
%    model=rmfield(model,'metNames');
end

if isfield(model,'metFormulas')
    for m=1:nMet
         model.met(m).formula=model.metFormulas{m};
    end
    model=rmfield(model,'metFormulas');
end

if isfield(model,'metCharges')
    for m=1:nMet
         model.met(m).charge=model.metCharges(m);
    end
    model=rmfield(model,'metCharges');
end

% %replace any dashes in metabolite names with underscores
% for m=1:nMet
%     x = strfind(model.met(m).abbreviation,'-');
%     if x~=0
%         %cobra v2 mets have underscores
%         model.met(m).abbreviation(x)='_';
%         %make cobra v1 mets have all underscores also
%         abbr=model.mets{m};
%         abbr(x)='_';
%         model.mets{m}=abbr;
%     end
% end

%finds the reactions in the model which export/import from the model
%boundary i.e. mass unbalanced reactions
%e.g. Exchange reactions
%     Demand reactions
%     Sink reactions
model=findSExRxnInd(model);

fprintf('%s\n\n','...finished converting.')

author: opencobra

repo: cobratoolbox


text:
function T = myones(sizes)
% MYONES Like the built-in ones, except myones(k) produces a k*1 vector instead of a k*k matrix,
% T = myones(sizes)

if length(sizes)==0
  T = 1;
elseif length(sizes)==1
  T = ones(sizes, 1);
else
  T = ones(sizes(:)');
end


text:
function scanManagementGui(command, varargin)
%"scanManagementGui" GUI
%   Create a GUI to manage doses.
%
%   APA 11/21/05
%
%Usage:
%   scanManagementGui()
%  based on doseManagementGui.m by JRA
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

%The cursed globals. Need em.
global planC stateS
indexS = planC{end};

%Use a static window size, by pixels.  Do not allow resizing.
screenSize = get(0,'ScreenSize');
y = 380;
x = 640;
units = 'normalized';

%Height of a single row of text.
rowHeight = .06;

%If no command given, default to init.
if ~exist('command','var') || (exist('command','var') && isempty(command))
    command = 'init';
end

%Find handle of the gui figure.
% h = findobj('tag', 'scanManagementGui');
h = [];
if isfield(stateS.handle,'scanManagementFig') && ishandle(stateS.handle.scanManagementFig)
    h = stateS.handle.scanManagementFig;
end

%Set framecolor for uicontrols and pseudoframes.
frameColor = [0.8314 0.8157 0.7843];

switch upper(command)
    case 'INIT'
        %If gui doesnt exist, create it, else refresh it.
        if isempty(h)
            %Set up a new GUI window.
            h = figure('units', 'pixels', 'position',[(screenSize(3)-x)/2 (screenSize(4)-y)/2 x y], 'MenuBar', 'none', 'NumberTitle', 'off', 'resize', 'off', 'Tag', 'scanManagementGui', 'Color', [.75 .75 .75], 'WindowButtonUpFcn', 'scanManagementGui(''FIGUREBUTTONUP'')');
            stateS.handle.scanManagementFig = h;
            set(h, 'Name','Scan Management');

            %Create pseudo frames.
            axes('Position', [.02 .05 + (rowHeight + .02) .96 .87 - (rowHeight + .02)/2], 'color', frameColor, 'ytick',[],'xtick',[], 'box', 'on');
            line([.5 .5], [0 1], 'color', 'black');
            axes('Position', [.02 .03 .96 rowHeight+.02], 'color', frameColor, 'ytick',[],'xtick',[], 'box', 'on');

            %Create controls for displaying the GUI status.
            ud.handles.status = uicontrol(h, 'units', units, 'string', 'Status:', 'Position', [.04 .03 .05 rowHeight], 'Style', 'text', 'BackgroundColor', frameColor);
            ud.handles.status = uicontrol(h, 'units', units, 'string', '', 'Position', [.1 .03 .8 rowHeight], 'Style', 'text', 'HorizontalAlignment', 'left', 'foregroundcolor', 'red', 'BackgroundColor', frameColor);
            ud.currentScan = 1;
            set(h, 'userdata', ud);
        end
        scanManagementGui('refresh', h);
        figure(h);

    case 'REFRESHFIELDS'
        %Refresh the field list for the current scan.
        ud = get(h, 'userdata');
        scanNum = ud.currentScan;
        set(ud.handles.name, 'string', planC{indexS.scan}(scanNum).scanType);
        %         set(ud.handles.units, 'string', planC{indexS.scan}(scanNum).scanUnits);
        [dA, isCompress, isRemote] = getScanArray(scanNum, planC);
        if isCompress
            set(ud.handles.compbutton, 'string', 'Decompress');
            set(ud.handles.storageMethod, 'string', 'Memory, Compressed');
        else
            set(ud.handles.compbutton, 'string', 'Compress');
        end

        if isRemote
            set(ud.handles.remotebutton, 'string', 'Use Memory');
            set(ud.handles.storageMethod, 'string', 'On Disk');
        else
            set(ud.handles.remotebutton, 'string', 'Use Disk');
        end

        if ~isCompress & ~isRemote
            set(ud.handles.storageMethod, 'string', 'Memory');
        end
        %         set(ud.handles.maxscan, 'string', ud.maxDoses{scanNum});
        scanSize = getByteSize(planC{indexS.scan}(scanNum));
        set(ud.handles.scansize, 'string', [num2str(scanSize/(1024*1024), '%6.2f') 'MB']);
        % ESpezi MAY 2013
        % refresh scan date and time
        if isfield(planC{indexS.scan}(scanNum).scanInfo(1),'DICOMHeaders')
            scanDate = '';
            if isfield(planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders,'AcquisitionDate')
                scanDate = planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders.AcquisitionDate;
            end
            if ~isempty(scanDate)
                set(ud.handles.scanDate, 'string', datestr(datenum(scanDate,'yyyymmdd'),2));
            else
                set(ud.handles.scanDate, 'string', scanDate);
            end
            scanTime = '';
            if isfield(planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders,'AcquisitionTime')
                scanTime = planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders.AcquisitionTime;
            end
            if ~isempty(scanTime)
                [token, remain] = strtok(num2str(scanTime),'.');
                if ~isempty(strfind(token,':'))
                    [~,aqTime] = strtok(datestr(datenum(token,'HH:MM:SS')));
                else
                    [~,aqTime] = strtok(datestr(datenum(token,'HHMMSS')));
                end
                set(ud.handles.scanTime, 'string', aqTime);
            else
                set(ud.handles.scanTime, 'string', scanTime);
            end
        else
            % RTOG scan
            scanDate = planC{indexS.scan}(scanNum).scanInfo(1).scanDate;
            set(ud.handles.scanDate, 'string', scanDate);
        end

    case 'REFRESH'
        %Recreate and redraw the entire scanManagementGui.
        if isempty(h)
            return;
        end

        %Save the current figure so focus can be returned.
        hFig = gcf;

        %Focus on scanManagementGui for the moment.
        set(0, 'CurrentFigure', h);
        ud = get(h, 'userdata');
        try
            %Scrap all old buttons and axes.
            delete(ud.handles.thumbaxis);
            delete(ud.handles.name);
            delete(ud.handles.scansize);
            delete(ud.handles.storageMethod);
            % Espezi MAY 2013
            delete(ud.handles.scanDate);
            delete(ud.handles.scanTime);
            %delete(ud.handles.maxscan);
            %delete(ud.handles.units);
            delete(ud.handles.compbutton);
            delete(ud.handles.copybutton);
            delete(ud.handles.delbutton);
            ud.handles.thumbaxis = [];
            ud.handles.name = [];
            ud.handles.scansize = [];
            ud.handles.storageMethod = [];
            % Espezi MAY 2013
            ud.handles.scanDate = [];
            ud.handles.scanTime = [];
            %ud.handles.maxscan = [];
            %ud.handles.units = [];
            ud.handles.compbutton = [];
            ud.handles.copybutton = [];
            ud.handles.delbutton = [];
        end

        nScans      = length(planC{indexS.scan});
        if ud.currentScan > nScans
            ud.currentScan = nScans;
        end

        %Downsample colormap, redraws much faster.
        % cM = CERRColorMap(stateS.optS.doseColormap);
        cM = CERRColorMap('gray');
        ud.cM = cM;

        %         n  = size(cM, 1);
        %         newSize = 32;
        %         interval = (n-1) / 32;
        %         b = interp1(1:n, cM(:,1), 1:interval:n);
        %         c = interp1(1:n, cM(:,2), 1:interval:n);
        %         d = interp1(1:n, cM(:,3), 1:interval:n);
        %         ud.cM = [b' c' d'];

        %Setup thumbnail pane, with NxN axes.
        dx = 1/x; %pixel width in x,
        dy = 1/y; %pixel width in y, for margins.
        thumbRegion = [.52 .17 .44 .75];
        subPlotSize = ceil(sqrt(nScans));
        dh = thumbRegion(4)/subPlotSize;
        dw = thumbRegion(3)/subPlotSize;
        for i=1:subPlotSize^2
            row = subPlotSize - ceil(i/subPlotSize) + 1;
            col = mod(i-1,subPlotSize)+1;
            ud.handles.thumbaxis(i) = axes('position', [thumbRegion(1) + dw*(col-1) thumbRegion(2) + dh*(row-1) dw-dx dh-dy]);
            set(ud.handles.thumbaxis(i), 'ytick',[],'xtick',[], 'color', 'black', 'box', 'on', 'xcolor', 'white', 'ycolor', 'white');
            colormap(ud.handles.thumbaxis(i), ud.cM);
        end

        for i=1:nScans
            set(ud.handles.thumbaxis(i), 'ButtonDownFcn', ['scanManagementGui(''CHANGESCAN'', ' num2str(i) ');']);
            maxScan{i} = num2str(drawThumb(ud.handles.thumbaxis(i), planC, i, h));
            ud.previewSlice(i) = 1; %%%%%%%%%%%%%%%%%
        end

        txtLeft = .05;
        textWidth = .13;
        fieldLeft = .35;
        fieldWidth = .13;

        %Make text to describe uicontrols.
        uicontrol(h, 'units',units,'Position',[txtLeft 1-.15 textWidth rowHeight],'String', 'Scan name:', 'Style', 'text', 'horizontalAlignment', 'left', 'BackgroundColor', frameColor);
        uicontrol(h, 'units',units,'Position',[txtLeft 1-.21 textWidth rowHeight],'String', 'RAM Used:', 'Style', 'text', 'horizontalAlignment', 'left', 'BackgroundColor', frameColor);
        uicontrol(h, 'units',units,'Position',[txtLeft 1-.27 textWidth rowHeight],'String', 'Storage Method:', 'Style', 'text', 'horizontalAlignment', 'left', 'BackgroundColor', frameColor);
        % ESpezi MAY 2013
        uicontrol(h, 'units',units,'Position',[txtLeft 1-.33 textWidth rowHeight],'String', 'Acquisition date:', 'Style', 'text', 'horizontalAlignment', 'left', 'BackgroundColor', frameColor);
        uicontrol(h, 'units',units,'Position',[txtLeft 1-.39 textWidth rowHeight],'String', 'Acquisition time:', 'Style', 'text', 'horizontalAlignment', 'left', 'BackgroundColor', frameColor);
        % No max Scan needed like max Dose
        %         uicontrol(h, 'units',units,'Position',[txtLeft 1-.33 textWidth rowHeight],'String', 'Max Dose:', 'Style', 'text', 'horizontalAlignment', 'left', 'BackgroundColor', frameColor);
        %         uicontrol(h, 'units',units,'Position',[txtLeft 1-.39 textWidth rowHeight],'String', 'Units:', 'Style', 'text', 'horizontalAlignment', 'left', 'BackgroundColor', frameColor);

        %Make uicontrols for managing the scans, and displaying info.
        ud.handles.name          = uicontrol(h, 'units',units,'Position',[fieldLeft-fieldWidth 1-.15+.02 fieldWidth*2 rowHeight-.01],'String','', 'Style', 'edit', 'callback', 'scanManagementGui(''NAMEFIELD'');', 'userdata', i, 'horizontalAlignment', 'right');
        ud.handles.scansize      = uicontrol(h, 'units',units,'Position',[fieldLeft 1-.21 fieldWidth rowHeight],'String', '',  'Style', 'text', 'horizontalAlignment', 'right', 'BackgroundColor', frameColor);
        ud.handles.storageMethod = uicontrol(h, 'units',units,'Position',[fieldLeft-.05 1-.27 fieldWidth+.05 rowHeight],'String', '',  'Style', 'text', 'horizontalAlignment', 'right', 'BackgroundColor', frameColor);
        ud.handles.scanDate = uicontrol(h, 'units',units,'Position',[fieldLeft-.05 1-.33 fieldWidth+.05 rowHeight],'String', '',  'Style', 'text', 'horizontalAlignment', 'right', 'BackgroundColor', frameColor);
        ud.handles.scanTime = uicontrol(h, 'units',units,'Position',[fieldLeft-.05 1-.39 fieldWidth+.05 rowHeight],'String', '',  'Style', 'text', 'horizontalAlignment', 'right', 'BackgroundColor', frameColor);
        % No max Scan needed like max Dose
        %         ud.handles.maxdose       = uicontrol(h, 'units',units,'Position',[fieldLeft 1-.33+.02 fieldWidth rowHeight-.01],'String', '', 'Style', 'edit', 'callback', 'scanManagementGui(''DOSESCALE'');', 'userdata', i, 'horizontalAlignment', 'right');
        %         ud.handles.units         = uicontrol(h, 'units',units,'Position',[fieldLeft 1-.39 fieldWidth rowHeight],'String', '',  'Style', 'text', 'horizontalAlignment', 'right', 'BackgroundColor', frameColor);

        ud.handles.remotebutton  = uicontrol(h, 'units',units,'Position',[fieldLeft 1-.59 fieldWidth rowHeight],'String', 'Make Remote', 'Style', 'pushbutton', 'callback', 'scanManagementGui(''REMOTE'');', 'userdata', i);
        ud.handles.compbutton    = uicontrol(h, 'units',units,'Position',[fieldLeft 1-.67 fieldWidth rowHeight],'String', 'Compress', 'Style', 'pushbutton', 'callback', 'scanManagementGui(''COMPRESS'');', 'userdata', i);
        ud.handles.copybutton    = uicontrol(h, 'units',units,'Position',[fieldLeft 1-.75 fieldWidth rowHeight],'String','Save', 'Style', 'pushbutton', 'callback', 'scanManagementGui(''SAVE'');', 'userdata', i);
        % Scan cannot be deleted
        ud.handles.delbutton     = uicontrol(h, 'units',units,'Position',[fieldLeft 1-.83 fieldWidth rowHeight],'String','Delete', 'Style', 'pushbutton', 'callback', 'scanManagementGui(''DELETE'');', 'userdata', i);

        ud.handles.previewAxis   = axes('position', [txtLeft 1-.83 .25 .25*x/y], 'parent', h, 'box', 'on', 'ytick', [], 'xtick', [], 'buttondownfcn', 'scanManagementGui(''PREVIEWBUTTONDOWN'')');
        ud.handles.previewSliceNum = uicontrol(h, 'units',units,'Position',[fieldLeft-0.03 1-.47 fieldWidth rowHeight],'String', '',  'Style', 'text', 'horizontalAlignment', 'left', 'BackgroundColor', frameColor);

        text(.5, .6, 'After selecting a scan, ', 'fontsize', 8, 'hittest', 'off', 'parent', ud.handles.previewAxis, 'horizontalAlignment', 'center');
        text(.5, .5, 'click and drag mouse up', 'fontsize', 8, 'hittest', 'off', 'parent', ud.handles.previewAxis, 'horizontalAlignment', 'center');
        text(.5, .4, 'and down here to preview.', 'fontsize', 8, 'hittest', 'off', 'parent', ud.handles.previewAxis, 'horizontalAlignment', 'center');

        ud.maxScans = maxScan;
        set(h, 'userdata', ud);
        set(0, 'CurrentFigure', hFig);

        if ~isempty(maxScan)
            scanManagementGui('REFRESHFIELDS');
        end

    case 'PREVIEWBUTTONDOWN'
        %Button clicked in the preview window.
        ud = get(h, 'userdata');
        ud.previewDown = 1;
        set(h, 'WindowButtonMotionFcn', 'scanManagementGui(''PREVIEWMOTION'')');
        set(h, 'userdata', ud)

    case 'FIGUREBUTTONUP'
        %Mouse up, if in preview window disable motion fcn.
        ud = get(h, 'userdata');
        if ~isfield(ud, 'previewDown') | ud.previewDown == 1;
            ud.previewDown = 0;
            set(h, 'WindowButtonMotionFcn', '');
            set(h, 'userdata', ud);
        end

    case 'PREVIEWMOTION'
        %Motion in the preview, with mouse down. Change preview slice.
        ud = get(h, 'userdata');
        cp = get(h, 'currentpoint');
        if isfield(ud, 'previewY')
            if ud.previewY > cp(2)
                ud.previewSlice(ud.currentScan) = ud.previewSlice(ud.currentScan)+1;%min(ud.previewSlice(ud.currentDose)+1, size(getDoseArray(ud.currentDose), 3));
                set(h, 'userdata', ud);
                scanManagementGui('refreshpreviewandfields');
            elseif ud.previewY < cp(2)
                ud.previewSlice(ud.currentScan) = ud.previewSlice(ud.currentScan)-1;%max(ud.previewSlice(ud.currentDose)-1,1);
                set(h, 'userdata', ud);
                scanManagementGui('refreshpreviewandfields');
            end
            ud = get(h, 'userdata');
            ud.previewY = cp(2);
        else
            ud.previewY = cp(2);
        end
        set(h, 'userdata', ud);

    case 'REFRESHPREVIEWANDFIELDS'
        %Refresh both the preview and the fields.  Not modular with RefreshFields
        %to ensure that only one call to getDoseArray is made.
        ud = get(h, 'userdata');
        scanNum = ud.currentScan;
        bdf = get(ud.handles.previewAxis, 'buttondownfcn');
        [dA, isCompress, isRemote] = getScanArray(scanNum, planC);
        minSlice = 1;
        maxSlice = size(dA, 3);
        previewSlice = clip(ud.previewSlice(scanNum), minSlice, maxSlice, 'limits');
        ud.previewSlice(scanNum) = previewSlice;
        set(ud.handles.previewSliceNum, 'string', num2str(previewSlice));
        maxScan = str2num(ud.maxScans{scanNum});
        cLim =  [0 maxScan];
        if maxScan == 0;
            cLim = [0 1];
        end
        imagesc(dA(:,:,previewSlice), 'parent', ud.handles.previewAxis, 'hittest', 'off');
        set(ud.handles.previewAxis, 'buttondownfcn', bdf, 'box', 'on', 'ytick', [], 'xtick', [], 'CLim', cLim);
        set(h, 'userdata', ud);

        %Refresh fields as well.
        set(ud.handles.name, 'string', planC{indexS.scan}(scanNum).scanType);
        %         set(ud.handles.units, 'string', planC{indexS.scan}(scanNum).scanUnits);
        if isCompress
            set(ud.handles.compbutton, 'string', 'Decompress');
            set(ud.handles.storageMethod, 'string', 'Memory, Compressed');
        else
            set(ud.handles.compbutton, 'string', 'Compress');
        end

        if isRemote
            set(ud.handles.remotebutton, 'string', 'Use Memory');
            set(ud.handles.storageMethod, 'string', 'On Disk');
        else
            set(ud.handles.remotebutton, 'string', 'Use Disk');
        end

        if ~isCompress & ~isRemote
            set(ud.handles.storageMethod, 'string', 'Memory');
        end

        %Refresh borders.
        set(ud.handles.thumbaxis, 'xcolor', 'white', 'ycolor', 'white');
        set(ud.handles.thumbaxis(scanNum), 'xcolor', 'yellow', 'ycolor', 'yellow');

        %         set(ud.handles.maxscan, 'string', ud.maxDoses{doseNum});
        scanSize = getByteSize(planC{indexS.scan}(scanNum));
        set(ud.handles.scansize, 'string', [num2str(scanSize/(1024*1024), '%6.2f') 'MB']);

        % ESpezi MAY 2013
        % refresh scan date and time
        if isfield(planC{indexS.scan}(scanNum).scanInfo(1),'DICOMHeaders')
            scanDate = '';
            if isfield(planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders,'AcquisitionDate')
                scanDate = planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders.AcquisitionDate;
            end
            if ~isempty(scanDate)
                set(ud.handles.scanDate, 'string', datestr(datenum(scanDate,'yyyymmdd'),2));
            else
                set(ud.handles.scanDate, 'string', scanDate);
            end
            scanTime = '';
            if isfield(planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders,'AcquisitionTime')
                scanTime = planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders.AcquisitionTime;
            end
            if ~isempty(scanTime)
                [token, remain] = strtok(num2str(scanTime),'.');
                if ~isempty(strfind(token,':'))
                    [~,aqTime] = strtok(datestr(datenum(token,'HH:MM:SS')));
                else
                    [~,aqTime] = strtok(datestr(datenum(token,'HHMMSS')));
                end
                set(ud.handles.scanTime, 'string', aqTime);
            else
                set(ud.handles.scanTime, 'string', scanTime);
            end
        else
            % RTOG scan
            scanDate = planC{indexS.scan}(scanNum).scanInfo(1).scanDate;
            set(ud.handles.scanDate, 'string', scanDate);
        end

    case 'CHANGESCAN'
        %New scan has been clicked on.
        ud = get(h, 'userdata');
        newScan = varargin{1};
        ud.currentScan = newScan;
        set(h, 'userdata', ud);
        scanManagementGui('refreshpreviewandfields');

    case 'NAMEFIELD'
        %Dose name has changed, update in planC.
        ud = get(h, 'userdata');
        scanNum = ud.currentScan;
        oldString = planC{indexS.scan}(scanNum).scanType;
        string = get(gcbo, 'string');
        planC{indexS.scan}(scanNum).scanType = string;
        statusString = ['Renamed scan number ' num2str(scanNum) ' from ''' oldString ''' to ''' string '''.'];
        scanManagementGui('status', statusString);

    case 'STATUS'
        %Display passed string in status bar.
        statusString = varargin{1};
        ud = get(gcbf, 'userdata');
        h = ud.handles.status;
        set(h, 'string', statusString);

    case 'COMPRESS'
        %Compress/decompress selected scan.
        ud = get(h, 'userdata');
        scanNum = ud.currentScan;
        scanName = planC{indexS.scan}(scanNum).scanType;

        if ~isCompressed(planC{indexS.scan}(scanNum).scanArray)
            statusString = ['Compressing scan number ' num2str(scanNum) ', ''' scanName ''' please wait...'];
            scanManagementGui('status', statusString);
            planC{indexS.scan}(scanNum).scanArray = compress(getScanArray(scanNum, planC));
            drawThumb(ud.handles.thumbaxis(scanNum), planC, scanNum, h);
            statusString = ['Compressed scan number ' num2str(scanNum)  ', ''' scanName '''.'];
            scanManagementGui('status', statusString);
            set(ud.handles.compbutton, 'string', 'Decompress');
        else
            statusString = ['Decompressing scan number ' num2str(scanNum) ', ''' scanName ''', please wait...'];
            scanManagementGui('status', statusString);
            %Use getScanArray and not decompress to use the cached value.
            planC{indexS.scan}(scanNum).scanArray = getScanArray(scanNum, planC);
            maxScan = drawThumb(ud.handles.thumbaxis(scanNum), planC, scanNum, h);
            %Update the maxScans value, but be sure to get a fresh ud since
            %a user could have clicked during compression.
            ud = get(h, 'userdata');
            ud.maxScans{scanNum} = num2str(maxScan);
            set(h, 'userdata', ud);
            statusString = ['Decompressed scan number ' num2str(scanNum)  ', ''' scanName '''.'];
            scanManagementGui('status', statusString);
            set(ud.handles.compbutton, 'string', 'Compress');
        end
        scanManagementGui('refreshpreviewandfields');

    case 'REMOTE'
        %Make/unmake selected dose remote.
        ud = get(h, 'userdata');
        scanNum = ud.currentScan;
        scanName = planC{indexS.scan}(scanNum).scanType;

        scanUID = planC{indexS.scan}(scanNum).scanUID;

        if isLocal(planC{indexS.scan}(scanNum).scanArray)
            statusString = ['Writing to disk scan number ' num2str(scanNum) ', ''' scanName ''' please wait...'];
            scanManagementGui('status', statusString);
            [fpath,fname] = fileparts(stateS.CERRFile);
            planC{indexS.scan}(scanNum).scanArray = setRemoteVariable(getScanArray(scanNum, planC), 'LOCAL',fullfile(fpath,[fname,'_store']),['scanArray_',scanUID,'.mat']);
            % Also make remote the scanArraySuperior and scanArrayInferior matrices
            planC{indexS.scan}(scanNum).scanArraySuperior = setRemoteVariable(getScanArraySuperior(scanNum, planC), 'LOCAL',fullfile(fpath,[fname,'_store']),['scanArraySuperior_',scanUID,'.mat']);
            planC{indexS.scan}(scanNum).scanArrayInferior = setRemoteVariable(getScanArrayInferior(scanNum, planC), 'LOCAL',fullfile(fpath,[fname,'_store']),['scanArrayInferior_',scanUID,'.mat']);
            drawThumb(ud.handles.thumbaxis(scanNum), planC, scanNum, h);
            statusString = ['Wrote to disk scan number ' num2str(scanNum)  ', ''' scanName '''.'];
            scanManagementGui('status', statusString);
            set(ud.handles.remotebutton, 'string', 'Use Memory');
            uiwait(msgbox(['scanArray stored in folder ',fullfile(fpath,[fname,'_store']),'. Note the Location'],'modal'));
        else
            statusString = ['Reading from disk scan number ' num2str(scanNum) ', ''' scanName ''', please wait...'];
            scanManagementGui('status', statusString);

            remotePath = planC{indexS.scan}(scanNum).scanArray.remotePath;
            filenam = planC{indexS.scan}(scanNum).scanArray.filename;
            stateS.reqdRemoteFiles(strcmp(fullfile(remotePath,filenam),stateS.reqdRemoteFiles)) = [];

            remotePathSup = planC{indexS.scan}(scanNum).scanArraySuperior.remotePath;
            filenameSup = planC{indexS.scan}(scanNum).scanArraySuperior.filename;
            stateS.reqdRemoteFiles(strcmp(fullfile(remotePathSup,filenameSup),stateS.reqdRemoteFiles)) = [];

            remotePathInf = planC{indexS.scan}(scanNum).scanArrayInferior.remotePath;
            filenameInf = planC{indexS.scan}(scanNum).scanArrayInferior.filename;
            stateS.reqdRemoteFiles(strcmp(fullfile(remotePathInf,filenameInf),stateS.reqdRemoteFiles)) = [];

            planC{indexS.scan}(scanNum).scanArray = getScanArray(scanNum,planC);
            planC{indexS.scan}(scanNum).scanArraySuperior = getScanArraySuperior(scanNum, planC);
            planC{indexS.scan}(scanNum).scanArrayInferior = getScanArrayInferior(scanNum, planC);

            if ~ismember(fullfile(remotePath,filenam),stateS.reqdRemoteFiles)
                delete(fullfile(remotePath,filenam))
            end
            if ~ismember(fullfile(remotePathSup,filenameSup),stateS.reqdRemoteFiles)
                delete(fullfile(remotePathSup,filenameSup))
            end
            if ~ismember(fullfile(remotePathInf,filenameInf),stateS.reqdRemoteFiles)
                delete(fullfile(remotePathInf,filenameInf))
            end

            %remove remote storage directory if it is empty
            dirRemoteS = dir(remotePath);
            if ~any(~cellfun('isempty',strfind({dirRemoteS.name},'.mat')))
                rmdir(remotePath)
            end

            maxScan = drawThumb(ud.handles.thumbaxis(scanNum), planC, scanNum, h);
            statusString = ['Read from disk scan number ' num2str(scanNum)  ', ''' scanName '''.'];
            scanManagementGui('status', statusString);
            set(ud.handles.remotebutton, 'string', 'Use Disk');

            %Update the maxDoses value, but be sure to get a fresh ud since
            %a user could have clicked during remote writing to disk.
            ud = get(h, 'userdata');
            ud.maxScans{scanNum} = num2str(maxScan);
            set(h, 'userdata', ud);
        end
        scanManagementGui('refreshpreviewandfields');

    case 'SAVE'
        %Open dialog to save dose array as .mat file.
        ud = get(h, 'userdata');
        scanNum = ud.currentScan;
        scanName = planC{indexS.scan}(scanNum).scanType;
        [filename, pathname] = uiputfile('*.mat', ['Save (uncompressed) scan array number ' num2str(scanNum) ' as:']);
        if filename==0
            return;
        end
        scan3D = getScanArray(scanNum, planC);
        save(fullfile(pathname, filename), 'scan3D');
        statusString = ['Saved scan number ' num2str(scanNum) ', ''' scanName ''' to ' [filename '.mat'] '.'];
        scanManagementGui('status', statusString);

    case 'DELETE'
        %Delete selected scan.  If being displayed, verify deletion with user.
        ud = get(h, 'userdata');
        scanNum = ud.currentScan;
        scanName = planC{indexS.scan}(scanNum).scanType;

        refreshViewer = 0;
        axesV = checkDisplayedScans(scanNum);
        if ~isempty(axesV)
            choice = questdlg('One or more CERR axes are currently displaying this scan.  If you delete it, these axes will be set to display no scan.  Proceed?', 'Continue?', 'Continue', 'Abort', 'Continue');
            if strcmpi(choice, 'Abort')
                statusString = ['Delete aborted.'];
                scanManagementGui('status', statusString);
                return;
            else
                %Set the axes scan value to null.
                %setAxisScanToNull(axesV);
            end
        end

        del = questdlg(['Are you sure you want to delete scan number' num2str(scanNum) ', ''' scanName ''' and the associated structures?'], 'Continue?', 'Continue', 'Abort', 'Continue');
        if strcmpi(del, 'Continue')
            statusString = ['Deleted scan number ' num2str(scanNum) ', ''' scanName '''.'];
            scanManagementGui('status', statusString);
            %Delete the structures associated with this scan
            assocScanV = getStructureAssociatedScan(1:length(planC{indexS.structures}), planC);
            structToDelete = find(assocScanV == scanNum);
            planC{indexS.structures}(structToDelete) = [];
            %Delete structureArray
            %indAssoc = find(strcmpi({planC{indexS.structureArray}.assocScanUID},planC{indexS.scan}(scanNum).scanUID));
            %planC{indexS.structureArray}(indAssoc) = [];
            %if ~isempty(structToDelete)
            if length(planC{indexS.structureArray}) >= scanNum
                planC{indexS.structureArray}(scanNum) = [];
                planC{indexS.structureArrayMore}(scanNum) = [];
            end
            %end
            stateS.structsChanged = 1;
            %Update doses associated with this scan
            while ~isempty(find(strcmpi({planC{indexS.dose}.assocScanUID},planC{indexS.scan}(scanNum).scanUID)))
                indAssoc = find(strcmpi({planC{indexS.dose}.assocScanUID},planC{indexS.scan}(scanNum).scanUID));
                n = indAssoc(1);
                transM = getTransM(planC{indexS.dose}(n),planC);
                planC{indexS.dose}(n).assocScanUID = [];
                planC{indexS.dose}(n).transM = transM;
            end
            %Delete the scan
            planC{indexS.scan}(scanNum) = [];

            %If scan below displayed scan deleted, its number has changed.
            if scanNum < stateS.scanSet
                stateS.scanSet = stateS.scanSet - 1;
                stateS.structSet = stateS.scanSet;
            end
            if stateS.scanSet == 0
                stateS.scanSet = 1;
                stateS.structSet = 1;
            end
            refreshViewer = 1;
            updateAxesForDeletedScan(scanNum);
            scanManagementGui;
        else
        end

        %Refresh CERR axes in the case of changed dose sets.
        if refreshViewer
            stateS.CTDisplayChanged = 1;
            sliceCallBack('refresh');
        end

    case 'DOSESCALE'
        ud          = get(h, 'userdata');
        doseNum     = ud.currentDose;
        doseName    = planC{indexS.dose}(doseNum).fractionGroupID;
        newMaxDose  = str2num(get(gcbo, 'string'));

        %Get the dose array and its compression state.
        [dA, isCompress, isRemote] = getDoseArray(doseNum, planC);

        maxScan = max(dA(:));
        if maxScan == 0
            set(gcbo, 'string', '0')
            statusString = ['Unable to rescale all zero dose distribution.'];
            scanManagementGui('status', statusString);
            return;
        end

        if isCompress | isRemote
            statusString = ['Rescaling compressed or remote dose can take a moment, please wait...'];
            scanManagementGui('status', statusString);
        end

        %Perform the rescale.
        dA = dA *  (1/maxDose * newMaxDose);

        planC = setDoseArray(doseNum, dA, planC);
        statusString = ['Rescaled dose number ' num2str(doseNum) ', ''' doseName ''' from [0 ' num2str(maxDose) '] Gy to [0 ' num2str(newMaxDose) '] Gy.'];
        ud.maxDoses{doseNum} = num2str(newMaxDose);
        set(h, 'userdata', ud);
        scanManagementGui('status', statusString);

    case 'QUIT'
        close;
end

function nBytes = getByteSize(data)
%"getByteSize"
%Returns the number of bytes in the passed data
infoStruct = whos('data');
nBytes = infoStruct.bytes;

function maxScan = drawThumb(hAxis, planC, index, hFigure)
%"drawThumb"
%In passed dose array, find slice with highest dose and draw in hAxis.
%Also denote the index in the corner.  If compressed show compressed.
set(hFigure, 'CurrentAxes', hAxis);
toDelete = get(hAxis, 'children');
delete(toDelete);

%Get the dose array and its compression state.
[dA, isCompress, isRemote] = getScanArray(index, planC);

bdf = get(hAxis, 'buttondownfcn');

maxScan = arrayMax(dA);
% 	maxLoc = find(dA == maxScan);
% 	[r,c,s] = ind2sub(size(dA), maxLoc(1));
% set the scan to median of z-values
indexS = planC{end};
[xV, yV, zV] = getScanXYZVals(planC{indexS.scan}(index));
s = ceil(median(1:length(zV)));
thumbImage = dA(:,:,s(1));
imagesc(thumbImage, 'hittest', 'off', 'parent', hAxis);
set(hAxis, 'ytick',[],'xtick',[]);

if isCompress & isRemote
    text(.1, .1, 'Compressed', 'units', 'normalized', 'fontsize', 8, 'color', 'white', 'horizontalAlignment', 'left', 'hittest', 'off', 'parent', hAxis);
    text(.1, .2, 'Remote', 'units', 'normalized', 'fontsize', 8, 'color', 'white', 'horizontalAlignment', 'left', 'hittest', 'off', 'parent', hAxis);
elseif isRemote
    text(.1, .1, 'Remote', 'units', 'normalized', 'fontsize', 8, 'color', 'white', 'horizontalAlignment', 'left', 'hittest', 'off', 'parent', hAxis);
elseif isCompress
    text(.1, .1, 'Compressed', 'units', 'normalized', 'fontsize', 8, 'color', 'white', 'horizontalAlignment', 'left', 'hittest', 'off', 'parent', hAxis);
end

xLim = get(hAxis, 'xlim');
yLim = get(hAxis, 'ylim');
x = (xLim(2) - xLim(1)) * .87 + xLim(1);
y = (yLim(2) - yLim(1)) * .15 + yLim(1);
text(x, y, num2str(index), 'fontsize', 8, 'color', 'white', 'hittest', 'off', 'parent', hAxis);
set(hAxis, 'buttondownfcn', bdf);
axis ij;
drawnow;
set(hAxis,'xcolor', 'white', 'ycolor', 'white')

function [nAxesV] = checkDisplayedScans(scanNum)
%"checkDisplayedDoses"
%   Check to see if any CERR axes are displaying the requested dose num.
%   Returns the numbers of any axes that are.

global stateS
nAxesV = [];

%Iterate over axes.
for i=1:length(stateS.handle.CERRAxis)
    %Get axis info for this axis.
    %aI = get(stateS.handle.CERRAxis(i), 'userdata');
    aI = stateS.handle.aI(i);
    if ismember(scanNum, aI.scanSets)
        nAxesV = union(nAxesV, i);
    end
end

function setAxisScanToNull(nAxesV)
%"setAxisDoseToNull"
%   Sets the passed axis number's userdata doseSet fields to null if they
%   were using 'manual'.

global stateS

%Iterate over axes.
for i=1:length(nAxesV)
    %Get axis info for this axis.
    % aI = get(stateS.handle.CERRAxis(nAxesV(i)), 'userdata');
    aI = stateS.handle.aI(nAxesV(i));
    if strcmpi(aI.scanSelectMode, 'manual')
        aI.scanSets = [];
        set(stateS.handle.CERRAxis(nAxesV(i)), 'userdata', aI);
    end
end

function updateAxesForDeletedScan(delIndex)
%"updateAxesForDeletedDose"
%   Shifts all doseNums being displayed in CERR axes to account for deleted
%   dose distributions, given the number of the deleted dose.

global stateS

%Iterate over axes.
for i=1:length(stateS.handle.CERRAxis)

    %Get axis info for this axis.
    %aI = get(stateS.handle.CERRAxis(i), 'userdata');
    aI = stateS.handle.aI(i);
    if ~isempty(aI.scanObj)
        %scanSets = aI.scanObj.scanSet;
        if strcmpi(aI.scanSelectMode,'auto')
            aI.scanObj.scanSet = stateS.scanSet; %max(1,scanSets(scanSets >= delIndex)- 1);
        else
            if delIndex <= aI.scanObj.scanSet
                aI.scanObj.scanSet = stateS.scanSet;
            end
        end
        %set(stateS.handle.CERRAxis(i), 'userdata', aI);
        stateS.handle.aI(i) = aI;
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [X,S] = calculator2(D)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

global K_S S_0 mu_max Y_XS

X = Y_XS ./ 3 .* (S_0 - D .* K_S ./ (mu_max - D));

S = D .* K_S ./ (mu_max - D);

end


text:
function ar=lpcls2ar(ls)
%LPCLS2AR convert line spectrum pair frequencies to ar polynomial AR=(LS)
% input vector elements should be in the range 0 to 0.5

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

[nf,p]=size(ls);
p1=p+1;
p2 = p1*2;
ar=zeros(nf,p1);
for k=1:nf
  le=exp(ls(k,:)*pi*2i);
  lf=[1 le -1 conj(fliplr(le))];
  y=real(poly(lf(1:2:p2)));
  x=real(poly(lf(2:2:p2)));
  ar(k,:)=(x(1:p1)+y(1:p1))/2;
end

author: decouples


text:
function root_num = p01_root_num ( )

%*****************************************************************************80
%
%% P01_ROOT_NUM returns the number of known roots for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer ROOT_NUM, the number of known roots.
%
  root_num = 3;

  return
end


text:
% V = est_calcInvCovMatFourierPDC(Rinv,E,foi,fs,N,p, verb)
%
% Obtain the frequency domain transform of the inverse covariance matrix of
% an M-variate VAR[p] process
% This is needed for calculating the analytic statistics for the PDC [1]
%
% Inputs:
%
%     Rinv: inverse process covariance matrix obtained from est_calcInvCovMat
%     E:    noise covariance matrix
%     foi:  frequencies of interest (Hz)
%     fs:   sampling rate
%     N:    # chans
%     p:    model order
%     verb: verbosity level. 0 = no output, 1=text.
%
% Outputs:
%
%     V:    Frequency-domain transform of the inverse covariance matrix
%
% References:
%
% [1] Schelter et al, (2009). Testing for directed influences among neural
% signals using partial directed coherence. J. Neuroscience Methods. 152:210-9.
%
% see also: est_calcInvCovMat(), est_calcInvCovMatFourier()
%
% Author: Tim Mullen 2010, SCCN/INC, UCSD
% Email:  tim@sccn.ucsd.edu

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

function V = est_calcInvCovMatFourierPDC(Rinv,E,foi,fs,N,p,verb,DEBUG)

if nargin<8
    DEBUG = false;
end

%% extract the diagonal elements of H=Rinv
% structure of Hd is (e.g., for p=2):
% [diag(H(1,1)), diag(H(2,1)), diag(H(1,2)), diag(H(2,2))]
% where diag(H(u,v)) is the column vector formed by the diagonal of
% submatrix H(u,v) of H
Hd = zeros(N,p^2);
cnt=1;
for v=1:p
    for u=1:p
        Hd(:,cnt)=diag(Rinv((u-1)*N+1:u*N,(v-1)*N+1:v*N));
        cnt=cnt+1;
    end
end

if DEBUG
    for j=1:N
        try chol(reshape(Hd(j,:),[p p])); catch; fprintf('Hj<0: j=%d - ',j); keyboard; end;
    end
end

%% create u,v index vectors
% us = [1 2 ... p 1 2 ... p ... ]'   (p^2 length)
% vs = [1 1 ... 1 2 2 ... 2 ... ]'   (p^2 length)
us = repmat(1:p,1,p)';
vs = zeros(p^2,1);
for ii=1:p
    vs((ii-1)*p+1:(ii-1)*p+p)=ones(p,1)*ii;
end

%% construct the V matrix for all freqs
fi=0;
COS=zeros(p^2,1);
freqs=(2*pi*foi)/fs;
V = zeros(length(freqs),N,N); % NOTE: V will end up (N,N,freqs)
if verb, h=waitbar(0,'calculating V^-1...'); end
for f=freqs
    fi=fi+1;

    %% construct cosine matrix
    COS(:,1) = cos(us*f).*cos(vs*f)+sin(us*f).*sin(vs*f);
%     COS(:,2) = sin(us*f).*sin(vs*f);
%     COS(:,3) = cos(us*f).*sin(vs*f);
%     COS(:,4) = sin(us*f).*sin(vs*f);

    %% multiply Hd and COS matrices to get matrix where row j is
    % [sum_{u,v=1 : p} Hjj(u,v)COS_11(u,v), ...
    %  sum_{u,v=1 : p} Hjj(u,v)COS_21(u,v), ...
    %  sum_{u,v=1 : p} Hjj(u,v)COS_12(u,v), ...
    %  sum_{u,v=1 : p} Hjj(u,v)COS_22(u,v)]
    %
    % where COSab(u,v) is the a,bth element of the sine-cos transform
    % matrix evaluated at u,v.
    %
    % NOTE: reshaping row j to 2x2 yeilds a matrix proportional to V_ij(f)
    % for some specific i
    Vm = Hd*COS;

    % now we multiply in the variances of the i's (E(i,i)) to generate the
    % full V matrix
    Vm = kron(diag(E)',Vm);

    % next, reshape to desired structure (NxNx2x2)
    V(fi,:,:) = Vm';

    if verb, waitbar(fi/length(freqs),h); end
end % for freqs

% permute V: (chs,chs,freqs,2,2)
V = permute(V,[2 3 1]);
if verb, close(h); end

repo: REST


text:
function pde = elasticity3datapoly(para)
%% ELASTICITYDATA3 data for the elasticity problem in three dimensions
%
% Modified from elasticitydata by Huayi Wei.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Lame constants
if nargin == 0
    lambda = 1;
    mu = 1;
else
    if ~isstruct(para)
        exit('we need a struct data');
    end
    if ~isfield(para,'lambda') || isempty(para.lambda)
        lambda = 1;
    else
        lambda = para.lambda;
    end
    if ~isfield(para,'mu') || isempty(para.mu)
        mu = 1;
    else
        mu = para.mu;
    end
end

pde = struct('lambda',lambda,'mu',mu, 'f', @f, 'exactu',@exactu,'g_D',@g_D);
%%%%%% subfunctions %%%%%%
    function s = f(p)
% f = - mu \Delta u - (mu + lambda) grad(div u)
       x = p(:,1); y = p(:,2); z = p(:,2);
       f1 = lambda*(32*y.*z.*(y - 1).*(z - 1) + 64*x.*y.*(2*z - 1).*(y - 1) ...
           + 32*x.*z.*(2*y - 1).*(z - 1) + 64*y.*(2*z - 1).*(x - 1).*(y - 1) ...
           + 32*z.*(2*y - 1).*(x - 1).*(z - 1)) + 32*mu*z.*(z - 1).*(4*x.*y - 2*y - 3*x + x.^2 + 1) ...
           + 32*mu*y.*(y - 1).*(8*x.*z - 4*z - 5*x + x.^2 + 2) + 64*mu*y.*z.*(y - 1).*(z - 1);
       f2 = lambda*(64*x.*z.*(x - 1).*(z - 1) + 64*x.*y.*(2*z - 1).*(x - 1) ...
           + 16*y.*z.*(2*x - 1).*(z - 1) + 64*x.*(2*z - 1).*(x - 1).*(y - 1) ...
           + 16*z.*(2*x - 1).*(y - 1).*(z - 1)) + 64*mu*x.*(x - 1).*(4*y.*z - 2*z - 3*y + y.^2 + 1) ...
           + 16*mu*z.*(z - 1).*(4*x.*y - 6*y - 2*x + 4*y.^2 + 1) + 128*mu*x.*z.*(x - 1).*(z - 1);
       f3 = lambda*(128*x.*y.*(x - 1).*(y - 1) + 32*x.*z.*(2*y - 1).*(x - 1) ...
           + 16*y.*z.*(2*x - 1).*(y - 1) + 32*x.*(2*y - 1).*(x - 1).*(z - 1) ...
           + 16*y.*(2*x - 1).*(y - 1).*(z - 1)) + 32*mu*x.*(x - 1).*(4*y.*z - 6*z - 2*y + 4*z.^2 + 1) ...
           + 16*mu*y.*(y - 1).*(4*x.*z - 10*z - 2*x + 8*z.^2 + 1) + 256*mu*x.*y.*(x - 1).*(y - 1);
       s = [f1 f2 f3];
    end
    function u = exactu(p)
       x = p(:,1); y = p(:,2); z = p(:,2);
       u1 = 16*x.*(1-x).*y.*(1-y).*z.*(1-z);
       u2 = 32*x.*(1-x).*y.*(1-y).*z.*(1-z);
       u3 = 64*x.*(1-x).*y.*(1-y).*z.*(1-z);
       u = [u1 u2 u3];
    end
    function s = g_D(p)
       s = exactu(p);
    end
end

repo: ifem


text:
function [eq,dict] = exp_match(exp,form,dict)
% Check whether an expression is matched by a pattern.
% [Equals,Assignment] = exp_match(Expression, Pattern)
%
% The pattern may contain blanks, named sub-patterns and conditional expressions. Assignment is an
% optional output which contains sub-expressions matched by named parts of the pattern.
%
% See also:
%   exp_same
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-19

if ~exp_beginfun('symbolic') return; end

if ~exist('dict','var')
    dict = struct(); end

[eq,dict] = utl_match(exp,form,dict);

exp_endfun;


text:
function [TWARes] = TWA_MMA(Param, freq, ecg, q,r,s, stlen)
%OVERVIEW, This function checks for arrhythmia in the current analysis window across all
% leads. If no arrhythmia is detected it then calls sub functions which first align the beats in the analysis
% window and find correlations of each beat with a template, to determine if the ecg data is clean enough for analysis.
% If alignment is successful, TWAs are evaluated using the MMA method
% coupled with a surrogate non parametric reshuffling method used for
% reducing false positive TWA detections.
%
% INPUTS        MANDATORY           DESCRIPTION
%               Param               Structure contains parameters used for
%                                   computing TWAs.
%
%               freq                Sampling frequency of the ecg data,
%                                   needs to be 1000 Hz.
%
%               ecg                 N by M array of ecg data. Contains M
%                                   channels of ECG with N datapoints each.
%
%               q                   array of q point locations
%
%               r                   array of r point locations
%
%               s                   array of s point locations
%
%               stlen               scalar, estimate of s-t segment length
%
% OUTPUT
%
%               TWARes              Structure contains the result for TWA
%                                   computation. Contains the following fields,
%
%               .VAlt               TWA amplitude estimates for the analysis
%                                   window.
%
%               .VAlt_Sig           TWA amplitude estimates which are statistically
%                                   significant compared to the noise threshold
%                                   estimate.
%
%               .noise_median       median of the gamma distribution used to model the noise.
%
%               .noise_95           95th percentile of the gamma distribution
%                                   used to model the noise.
%
%               .VAltPt             The location of the point with maximum
%                                   difference between the average even and odd beats.
%
%   REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Written by Shamim Nemati,
%       editted by Ismail Sadiq on 10/26/2019.
%	COPYRIGHT (C) 2019
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information. The license may be found in
%       the Documents folder of the Physionet-Cardiovascular-Signal-Toolbox.

TWARes.successfull = false;

WIN_SIZE = Param.Interval; %e.g., 60beats
OVERLAP_SIZE = floor(WIN_SIZE/2);
TWARes=[];
Num_Ch = size(ecg, 2);

Align.fidBase = []; Align.fid = []; Align.valid = [];

if (length(r) < WIN_SIZE) % need at least one window of beats
    disp('TWAbyMMAOnAFile: Not enough beats...')
    TWARes.HR=NaN;
    TWARes.VAlt=NaN*ones(1,Num_Ch);TWARes.VAlt_Sig=NaN*ones(1,Num_Ch);
    TWARes.VAltPt=NaN*ones(1,Num_Ch);
    TWARes.Noise_Median= NaN*ones(1,Num_Ch);TWARes.Noise_95 = NaN*ones(1,Num_Ch);
    return;
end

i = 0;
for j = 1:(WIN_SIZE-OVERLAP_SIZE):length(r)-WIN_SIZE+1
    i=i+1;
    disp(['Analysis window index = ' num2str(i)])
    beats = j:j+WIN_SIZE-1;

    % Check for arrhythmia

    cur_r = r(beats); % Check current window of beats for arrhythmia
    features = AF_features(diff(cur_r),freq);
    Index_test = SVM_AFdetection_withoutTrainingModel(features,1);
    if (Index_test)
        disp('Arrhythmia detected...')
        TWARes.HR(i)=NaN;
        TWARes.VAlt(i, :)=NaN*ones(1,Num_Ch);TWARes.VAlt_Sig(i, :)=NaN*ones(1,Num_Ch);
        TWARes.VAltPt(i,:)=NaN*ones(1,Num_Ch);
        TWARes.Noise_Median(i, :) = NaN*ones(1,Num_Ch);TWARes.Noise_95(i, :) = NaN*ones(1,Num_Ch);
        % If arrhythmia detected continue to next window
        continue;
    end

    try
        % Align beats and check validity for analysis before computing TWAs
        [Align hrt_rate] = Align_Beats(beats, q(beats), r(beats), s(beats), stlen, ecg, 1000, Param);
    catch
        % If failed to align due to catch continue to next window.
        disp('TWAbyMMAOnAFile: Catch says alignment failed...')
        TWARes.HR(i)=NaN;
        TWARes.VAlt(i, :)=NaN*ones(1,Num_Ch);TWARes.VAlt_Sig(i, :)=NaN*ones(1,Num_Ch);
        TWARes.VAltPt(i,:)=NaN*ones(1,Num_Ch);
        TWARes.Noise_Median(i, :) = NaN*ones(1,Num_Ch);TWARes.Noise_95(i, :) = NaN*ones(1,Num_Ch);
        continue;
    end
    if (isempty(Align.fid))
        disp('TWAbyMMAOnAFile: alignment failed ...');
        TWARes.HR(i)=NaN;
        TWARes.VAlt(i, :)=NaN*ones(1,Num_Ch); TWARes.VAlt_Sig(i, :)=NaN*ones(1,Num_Ch);
        TWARes.VAltPt(i,:)=NaN*ones(1,Num_Ch);
        TWARes.Noise_Median(i, :)=NaN*ones(1,Num_Ch);TWARes.Noise_95(i, :) = NaN*ones(1,Num_Ch);
        continue;
    else
        disp('TWAbyMMAOnAFile: alignment succeeded ...');
    end

    disp('TWAbyMMAOnAFile: looking for alternans ...');
    TWARes.HR(i) = hrt_rate;
    try
        [TWARes.VAlt(i, :), TWARes.VAlt_Sig(i, :), TWARes.Noise_Median(i,:),TWARes.Noise_95(i,:),...
            TWARes.VAltPt(i,:)]= TWA_by_MMA(Param, TWARes, ecg, Align); % *********************** Estimates T-wave alternans from the beat Matrix
    catch
        disp('TWAbyMMAOnAFile: search for alternans failed ')
        TWARes.VAlt(i, :)=NaN*ones(1,Num_Ch); TWARes.VAlt_Sig(i, :)=NaN*ones(1,Num_Ch);
        TWARes.VAltPt(i,:)=NaN*ones(1,Num_Ch);
        TWARes.Noise_Median(i, :)=NaN*ones(1,Num_Ch);TWARes.Noise_95(i, :) = NaN*ones(1,Num_Ch);
        continue;
    end
end

TWARes.successfull = true;

return;

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
close all; clear all;

%% Setting of the problem
global s
pde = fonedata; % f = 1;
pde.L = 1;
option.theta = 0.3;
option.estType = 'star';
option.maxIt = 20;
option.maxN = 2e5;
option.solver = 'mg';
option.yrefinement = 10; % ratio of h_y/h_x
option.gNquadorder = 5;
option.tol = 1e-6;
[node,elem] = squaremesh([-1,1,-1,1],0.5);
[node,elem] = delmesh(node,elem,'x>0 & y<0');
bdFlag = setboundary(node,elem,'Dirichlet');

%% s = 0.2
s = 0.2;
afemfracLap(node,elem,pde,bdFlag,option);

%% s = 0.4
s = 0.4;
afemfracLap(node,elem,pde,bdFlag,option);

%% s = 0.6
s = 0.6;
afemfracLap(node,elem,pde,bdFlag,option);

%% s = 0.8
s = 0.8;
afemfracLap(node,elem,pde,bdFlag,option);


text:
function [s,Ori_s,Unmix_s,SetupStruc] = readData(Angle,SetupStruc,channel)
Num = length(Angle);
SetupStruc.Angle = Angle;
T60 = SetupStruc.T60;
Sign_compared = SetupStruc.Sign_compared;
Channel_Num = SetupStruc.Channel_Num;
Len = zeros(Num,1);
if(~exist('channel','var'))
    channel = 1;
end
if(SetupStruc.Sign_Mid==1)
    channel_line = 0:Channel_Num-1;
else
    channel_line = 1:Channel_Num;
end
for n = 1:Num
    if(SetupStruc.Height==1)
        dataPath = 'simu_data/DifferentHeight/T60_';
    else
        dataPath = 'simu_data/SameHeight/T60_';
    end
    for i = channel_line%0:Channel_Num-1
        wavName =strcat(dataPath,num2str(T60),'/',num2str(Angle(n)),'/Mic',num2str(i),'.wav');
        if(SetupStruc.Sign_Mid==1)
            eval(strcat('[s_temp', num2str(n), '(:,', num2str(i), '+1),fs] = audioread(''', wavName, ''');'));
        else
            eval(strcat('[s_temp', num2str(n), '(:,', num2str(i), '),fs] = audioread(''', wavName, ''');'));
        end
    end
    eval(strcat('Len(n) = size(s_temp', num2str(n),',1);'));
end
SetupStruc.fs = fs;
len = min(Len);
s = s_temp1(1:len,:);
Unmix_s(:,:,1) = s;
if(Num>1)
    Ori_s = zeros(len,Num);
    Ori_s(:,1) = s_temp1(1:len,channel);
else
    Ori_s = zeros(len,2);
    Ori_s_temp = audioread(strcat(dataPath,'0/',num2str(Angle(1)),'/Mic0.wav'));
    Ori_s(:,1) = Ori_s_temp(1:len);
    Ori_s(:,2) = s_temp1(1:len,channel);
end
for n = 2:Num
    eval(strcat('s_temp = s_temp', num2str(n), ';'));
    s = s+s_temp(1:len,:);
    Unmix_s(:,:,n) = s_temp(1:len,:);
    Ori_s(:,n) = s_temp(1:len,channel);
end
if(Sign_compared<0 && Num>1)
    for i = 1:Num
        %%%% In the last version, if Sign_compared>-1, choose the signal in
        %%%% the T60=0s, else in the select T60
%         Ori_s_temp = audioread(strcat(dataPath,'0/',num2str(Angle(i)),'/Mic',num2str(Sign_compared),'.wav'));
        %%%% In this version, if Sign_compared<0, chosse the original
        %%%% signal, else choose the select T60
        Ori_s_temp = audioread(strcat('simu_data/oral/',num2str(Angle(i)),'.wav'));
        Ori_s(:,i) = Ori_s_temp(1:len);
    end
end
if(SetupStruc.AddNoiseFlag==1)
    av_pow = mean( sum(s.^2,1)/len );       % Average mic power across all received signals.
    sigma_noise = sqrt( av_pow/(10^(SetupStruc.NoiseSNR/10)) );		% st. dev. of white noise component to achieve desired SNR.
    s = s + sigma_noise*randn(size(s));      % Add some random noise
end
% autoPlot(Ori_s,fs,-1);   %%%the 3rd parameter is used to control whether adjusting the range of y axis when equal -1
% autoPlot(s(:,1),fs);
return;

repo: Voice-Separation-and-Enhancement


text:
function title = p03_title ( )

%*****************************************************************************80
%
%% P03_TITLE returns a title for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'The gaussian function';

  return
end


%==============================================================================
% This code is part of the VAMPIRE app for the Matlab-based toolbox
%  FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR/FAIR.m
%==============================================================================
%
% VAMPIRE - Variational Algorithm for Mass-Preserving Image REgistration
%
%
%         xxxxxxxx                                           xxxx
%      xxxxxxxxxxxxxx                                     xxxxxxxxxxx
%    xxxxxxxxxxxxxxxxx                                 xxxxxxxxxxxxxxxx
%   xxxxxxxxxxxxxxxxxxxx                             xxxxxxxxxxxxxxxxxxxx
%  xxxxxxxxxxxxxxxxxxxxxxx                          xxxxxxxxxxxxxxxxxxxxxx
%       xxxxxxxxxxxxxxxxxxx                       xxxxxxxxxxxxxxxxxxxxx
%         xxxxxxxxxxxxxxxxxxx     xxx   xxx      xxxxxxxxxxxxxxxxxxx
%         xxxxxxxxxxxxxxxxxxxx    xxxxxxxxx    xxxxxxxxxxxxxxxxxxxx
%          xxxxxxxxxxxxxxxxxxxxx  xxxxxxxxx  xxxxxxxxxxxxxxxxxxxxx
%          xxxx  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
%                   xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx       xx
%                     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
%                      xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
%                      x       xxxxxxxxxxxxxxxx    xxx
%                               xxxxxxxxxxxxxx
%                               xxx  xxxxxxxx
%                               x          xx
%
%
% For a quick example, run: EV_3Dmouse_VAMPIRE.m
%
% Contents of VAMPIRE toolbox:
%
% VAMPIRENPIRobjFctn.m   - objective function for non-parametric registration
% VAMPIREPIRobjFctn.m    - objective function for non-parametric registration
% VAMPIREsolveGN_PCG.m   - matrix-free PCG solver for Gauss-Newton system
%
% Examples:
%
% EV_2DGaussian_VAMPIRE.m        - 2D example, non-parametric mass-preserving
% EV_2DGaussian_VAMPIRE_MLPIR.m  - 2D example, parametric mass-preserving
% EV_2DGaussianNoise_VAMPIRE.m   - 2D example, non-parametric mass-preserving
% EV_3Dmouse_VAMPIRE.m           - 3D example, non-parametric mass-preserving
% EV_3Dmouse_VAMPIRE_MLPIR.m     - 3D example, parametric mass-preserving
%
% ==================================================================================

function debit = contents
if nargout == 0, help(mfilename); return; end;

debit = {
          'VAMPIRENPIRobjFctn.m'
          'VAMPIREPIRobjFctn.m'
          'VAMPIREsolveGN_PCG.m'
          'contents.m'
          'README.md'
          'examples'
        };


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

%% Syntax
% |[Fq,Vq,Fc,Vc]=dualClad(F,V,shrinkFactor,cladMethod);|

%% Description
% This function creates patch data defining a cladding on a surface. The
% cladding is obtained by shrinking the input faces around their mean and
% by connecting the face sides to form new faces. The output faces cover
% the edges of the dual of the input surface (hence the name). The output
% is the set of shrunk faces (same type as input) and a set of new
% quadrilateral faces connecting the shrunk faces.
% The shrink factor, which can be a constant or a spatially varying metric
% on the nodes or faces, defines the face and edge shrink used. The
% clad method determines whether the output mesh is connected:
% 1: From shrunk face to shrunk face
% 2: From shrunk face to shrunk edge to shrunk face
% 3: From shrunk face to shrunk face through edge (computes intersection at
% edge which may not be the centre of edge. This method avoids potential
% "kinks" seen for method 2.

%% Examples
clear; close all; clc;

%%
% Plot settings

figStruct.ColorDef='black';
figStruct.Color='k';

%%
% Specify test surface. Alter settings to test for different geometries and
% surface types.

%Testing settings
testCase=2; %1= sphere, 2=bunny, 3=dino
meshType=1; %1=triangles, 2=quadrilaterials
cutMesh=0; %0=not cut, 1=cut in half

switch testCase
    case 1
        [F,V,~]=geoSphere(2,1); % Building a geodesic dome surface model
        shrinkFactor=0.25;
    case 2
        [F,V]=stanford_bunny('g'); %Bunny
        V_mean=mean(V,1);
        V=V-V_mean(ones(size(V,1),1),:);
        shrinkFactor=0.25;
    case 3
        [F,V]=parasaurolophus; %dino
        V_mean=mean(V,1);
        V=V-V_mean(ones(size(V,1),1),:);
        shrinkFactor=0.25;
    case 4
        defaultFolder = fileparts(fileparts(mfilename('fullpath')));
        pathName=fullfile(defaultFolder,'data','libSurf');

        dataStruct=load(fullfile(pathName,'sprocket.mat'));
        F=dataStruct.F;
        V=dataStruct.V;
        shrinkFactor=0.5;
    case 5
        defaultFolder = fileparts(fileparts(mfilename('fullpath')));
        pathName=fullfile(defaultFolder,'data','libSurf');

        dataStruct=load(fullfile(pathName,'enginePart_p1.mat'));
        F=dataStruct.F;
        V=dataStruct.V;
        shrinkFactor=0.5;
end

if meshType==2
    [F,V]=tri2quad(F,V);
end

if cutMesh==1
    logicKeep=V(:,1)<mean(V(:,1));
    logicKeep=all(logicKeep(F),2);
    F=F(logicKeep,:);
    [F,V]=patchCleanUnused(F,V);
end

%% Example 1: Explaining the clad method

%%
% Define |dualClad| settings, i.e. the cladMethod and the shrink factor.
cladMethods=[1 2 3]; %1= fact-to-face connections, 2=

%%
% Visualize results

cFigure;

for q=1:1:numel(cladMethods)

    [Fq,Vq,Fc,Vc]=dualClad(F,V,shrinkFactor,cladMethods(q));

    subplot(1,3,q); hold on;
    title(['Clad method: ',num2str(cladMethods(q))]);
    h(1)=gpatch(F,V,'kw','k',0.25);
    h(2)=gpatch(Fc,Vc,'gw','g',1);
    h(3)=gpatch(Fq,Vq,'rw','r',1);
    h(4)=plotV([Vc;Vq],'b.','MarkerSize',15);
    axisGeom;
    camlight headlight;
    view(2);
    hl=legend(h,{'Input','Shrunk faces','Quads','Output vertices'});
    hl.Location='SouthOutside';
end

drawnow;

%% Example 2: Explaining the shrink factor

%%
% Define |dualClad| settings, i.e. the cladMethod and the shrink factor.
cladMethod=3; %1= fact-to-face connections, 2=

shrinkFactors=linspace(0.1,0.75,3); %A range of shrink factors

%%
% Visualize results

cFigure;
subplot(2,2,1); hold on;
title('Input surface');
gpatch(F,V,'gw','g',1);

axisGeom;
camlight headlight;
view(2);

for q=1:1:numel(shrinkFactors)

    [Fq,Vq,Fc,Vc]=dualClad(F,V,shrinkFactors(q),cladMethod);

    subplot(2,2,q+1); hold on;
    title(['Shrink factor: ',num2str(shrinkFactors(q))]);
    gpatch(Fc,Vc,'rw','r',1);
    gpatch(Fq,Vq,'rw','r',1);
    gpatch(F,V,'kw','none',0.25);
    axisGeom;
    camlight headlight;
    view(2);
end

drawnow;

%% Example 3: Spatially varying shrink factors
% An animation will be created to show effect of a spatially varying shrink
% factor

cladMethod=3;

% Define spatially varying shrink factor
shrinkFactor=-V(:,1);
shrinkFactor=mean(shrinkFactor(F),2);
shrinkFactor=shrinkFactor-min(shrinkFactor(:));
shrinkFactor=shrinkFactor./max(shrinkFactor(:));
shrinkFactor=shrinkFactor*1.2;
shrinkFactor=shrinkFactor+0.05;
shrinkFactor(shrinkFactor>1)=1;

[Fq,Vq,Fc,Vc]=dualClad(F,V,shrinkFactor,cladMethod);

%%
% Animating the effect of the shrink factor

%%
% Initialize scene
hf=cFigure(figStruct);
hold on;

hp1=gpatch(Fc,Vc,shrinkFactor,'k',1);
hp2=gpatch(Fq,Vq,'r','none',1);

axisGeom;
camlight headlight;
colormap(viridis(250));
clim([0 1]);
drawnow;
axis off
axis manual;

%%
% Animate scene

nSteps=40; %Number of animation steps
animStruct.Time=linspace(0,1,nSteps); %Create the time vector

t=linspace(0,2,nSteps);
shrinkFactor=-V(:,1);
shrinkFactor=mean(shrinkFactor(F),2);
shrinkFactor=shrinkFactor-min(shrinkFactor(:));
shrinkFactor=shrinkFactor./max(shrinkFactor(:));
shrinkFactor=shrinkFactor+1;
minLevel=0.05;

for q=1:1:nSteps

      shrinkFactorNow=shrinkFactor;
      shrinkFactorNow=shrinkFactorNow-t(q);
      shrinkFactorNow(shrinkFactorNow<minLevel)=minLevel;
      shrinkFactorNow(shrinkFactorNow>1)=1;

    [Fq,Vq,Fc,Vc]=dualClad(F,V,shrinkFactorNow,cladMethod);

    %Set entries in animation structure
    animStruct.Handles{q}=[hp1 hp2 hp1]; %Handles of objects to animate
    animStruct.Props{q}={'Vertices','Vertices','CData'}; %Properties of objects to animate
    animStruct.Set{q}={Vc,Vq,shrinkFactorNow}; %Property values for to set in order to animate
end

anim8(hf,animStruct); %Initiate animation

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

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
classdef mme_xfmr3p_opf < mp.mme_xfmr3p

%   MATPOWER
%   Copyright (c) 2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end

    methods
        function x0 = interior_x0(obj, mm, nm, dm, x0)
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%% Global Robust Adaptive Dynamic Programming
global W W1 F Q noise_on rho
rho = 0.5 ; % Robust Redesign Gain

F = [0 -3/2 -1/2];
Q = [5 0 0;0 0 0;0 0 0];
W = [2 -1.4 -.45];
wsave = W;
W1 = W;
P = eye(2)*10;
Pold = -100*eye(2);
noise_on=1;

Psave=[];
Qsave=[];
ua=[];

Trjsave=[];
Tsave=[];

T=0.005;
xinit=-10;
rinit=5;
x=[xinit;zeros(9,1);rinit]';
for i=0:9
    %CXX=[];
    C1=[];
    C2=[];
    C3=[];
    CQ=[];
    %CZZ=[];
    %x=[1;zeros(9,1)]';
    for j=0:49
        [t,x] = ode45(@polysys,[j*T,j*T+T]+50*i*T,[x(end,1) zeros(1,9) x(end,11)]');
        %[t,x] = ode_yuri(j*T,j*T+T,[x(end,1) zeros(1,9)]',0.001);
        C1 = [C1; 1/2*(x(end,1)^2-x(1,1)^2) 1/3*(x(end,1)^3-x(1,1)^3) 1/4*(x(end,1)^4-x(1,1)^4)];
        C2 = [C2; x(end,2:6)];
        C3 = [C3; x(end,7:9)];
        CQ = [CQ; x(end,10)];
        Tsave=[Tsave;t];
        Trjsave=[Trjsave;x(:,[1 11])];
        for k=1:length(t)
            ua=[ua; W(:)'*x(k,1).^[1 2 3]'+(0.01*sin(10*t(k))+0.01*sin(3*t(k))+0.01*sin(100*t(k)))*noise_on];
        end
    end

    if norm(P(:)-Pold(:))>0.1
        cvx_begin sdp
        variable Wn(3,1)
        variable dQs(3,3) symmetric
        Qv=-([C2 -C3]'*[C2 -C3])\[C2 -C3]'*(CQ+C1*Wn(:));
        dQs(1,1)==Qv(1);
        dQs(1,2)+dQs(2,1)==Qv(2);
        dQs(1,3)+dQs(3,1)+dQs(2,2)==Qv(3);
        dQs(3,2)+dQs(2,3)==Qv(4);
        dQs(3,3)==Qv(5);
        dQs>=0;
        Pn = [1/2*(Wn(1)) 1/6*(Wn(2));  1/6*(Wn(2)) 1/4*(Wn(3))];
        Pn<=P;

        %minimize([-1 100]*Pn*[-1 ;100]+[1 100]*Pn*[1 ;100])
        minimize(Pn(1,1)+Pn(2,2))
        W=Qv(6:8);
        wsave=[wsave;W(:)' ];
        cvx_end
        noise_on=1;
        Psave=[Psave;P(:)'];
        Qsave=[Qsave;dQs(:)'];
        Pold=P;
        P=Pn;
    else
        noise_on=0;
        disp(num2str(i))
    end

end
Psave
Qsave;
%%
figure(1)
[t0,y0] = ode45(@polysys0,[0 Tsave(end)],[xinit, rinit]);
for i=1:length(t0)
 u0(i) = W1*y0(i,1).^[1 2 3]';
end
plot(Tsave,Trjsave(:,1), 'b-', t0,y0(:,1), 'r-.', 'linewidth', 2)
legend('With GRADP-based controller', 'With initial controller')
xlabel('time (sec)')
ylabel('\phi')
ylim([-12 5])
% Create textarrow
annotation(figure(1),'textarrow',[0.267857142857143 0.228200808625336],...
	[0.245238095238095 0.433832086450543],'String',{'1st iteration'},...
	'FontSize',12);

% Create textarrow
annotation(figure(1),'textarrow',[0.351785714285714 0.291071428571429],...
	[0.511904761904762 0.666666666666667],'String',{'2nd iteration'},...
	'FontSize',12);

% Create textarrow
annotation(figure(1),'textarrow',[0.38409703504043 0.386738544474392],...
	[0.763163519772001 0.700619878874244],'String',{'3rd iteration'},...
	'FontSize',12);

% Create textarrow
annotation(figure(1),'textarrow',[0.603571428571428 0.535714285714285],...
	[0.533333333333334 0.673809523809524],'String',{'5th (final) iteration'},...
	'FontSize',12);

% Create textarrow
annotation(figure(1),'textarrow',[0.480357142857143 0.427966101694915],...
	[0.528571428571429 0.669064748201439],'String',{'4th iteration'},...
	'FontSize',12);

figure(2)
plot(Tsave,Trjsave(:,2), 'b-', t0,y0(:,2), 'r-.', 'linewidth', 2)
legend('With GRADP-based controller', 'With initial controller')
xlabel('time (sec)')
ylabel('r')
ylim([-1 5])
% Create textarrow
annotation(figure(2),'textarrow',[0.243684992570579 0.1996691805209],...
	[0.676982591876209 0.278255039384132],'String',{'1st iteration'},...
	'FontSize',12);

% Create textarrow
annotation(figure(2),'textarrow',[0.325408618127786 0.287523619950097],...
	[0.560928433268859 0.281012459180867],'String',{'2nd iteration'},...
	'FontSize',12);

% Create textarrow
annotation(figure(2),'textarrow',[0.407132243684993 0.373938714289719],...
	[0.502901353965184 0.289622365749069],'String',{'3rd iteration'},...
	'FontSize',12);

% Create textarrow
annotation(figure(2),'textarrow',[0.473997028231798 0.450802097059071],...
	[0.433268858800774 0.291719058253785],'String',{'4th iteration'},...
	'FontSize',12);

% Create textarrow
annotation(figure(2),'textarrow',[0.592867756315007 0.537890044576523],...
	[0.586073500967118 0.286266924564797],'String',{'5th (final) iteration'},...
	'FontSize',12);

% %
syms v(y)
%vsx=dsolve(diff(v)*(F(1)*y+F(2)*y^2+F(3)*y^3)+0.01*(y^2+y^4)-1/4*(diff(v))^2==0, v(0)==0)

%x=-1.5:0.05:1.5;
x=-2.5:.01:2.5;
vn=[];
v1=[];
vs=[];
us=[];
u1=[];
un=[];
P1=[Psave(2,1) Psave(2,2);Psave(2,3) Psave(2,4)] ;

vsxt=0;
for i=1:length(x)-1
y=x(i);
vn=[vn [y y^2]*P*[y ;y^2]];
v1=[v1 [y y^2]*P1*[y ;y^2]];
%vsx=(y^2 + 2)^(3/2)/15 - (2*2^(1/2))/15 - y^2/10;

vsx = 2.0*y*(0.25*y^4 + 1.5*y^3 + 2.25*y^2 + 5.0)^(1/2) - 3.0*y^2 - 1.0*y^3;
vsxt = vsxt+vsx*(x(i+1)-x(i));
% vsx=y^3/150 + (101*y^2 + 100)^(3/2)/15150 - 20/303;
vs=[vs vsxt];
end

figure(3)
plot(x(2:end),v1,'g:',x(2:end),vn,'r-.',x(2:end),vs+21.1238,'b','linewidth',2)
legend('V_1 : Initial cost', 'V^5: Improved cost', 'V^o: Optimal cost')
xlabel('\phi')

% figure(4)
% plot(x,u1,'g:',x,un,'r-.',x,us,'b','linewidth',2)
% legend('u_1: Initial control policy', 'u^5: Improved control policy', 'u^o: Optimal control policy')

author: yu-jiang

repo: radpbook

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


text:
function seq = sample_dbn(bnet, varargin)
% SAMPLE_DBN Generate a random sequence from a DBN.
% seq = sample_dbn(bnet, ...)
%
% seq{i,t} contains the values of the i'th node in the t'th slice.
%
% Optional arguments:
%
% length - length of sequence  to be generated (can also just use sample_dbn(bnet,T))
% stop_test - name of a function which is used to decide when to stop;
%   This will be called as feval(stop_test, seq(:,t))
%   i.e., stop_test is passed a cell array containing all the nodes in the current slice.
% evidence - initial evidence; if evidence{i,t} is non-empty, this node won't be sampled.

args = varargin;
nargs = length(args);

if (nargs == 1) & ~isstr(args{1})
  % Old syntax: sample_dbn(bnet, T)
  T = args{1};
else
  % get length
  T = 1;
  for i=1:2:nargs
    switch args{i},
     case 'length',      T = args{i+1};
     case 'evidence',    T = size(args{i+1}, 2);
    end
  end
end

ss = length(bnet.intra);
% set default arguments
seq = cell(ss, T);
stop_test = [];
for i=1:2:nargs
  switch args{i},
   case 'evidence',    seq = args{i+1}; % initialise observed nodes
   case 'stop_test',  stop_test = args{i+1};
  end
end

t = 1;
for i=1:ss
  if ~isempty(stop_test) | isempty(seq{i,t})
    ps = parents(bnet.dag, i);
    e = bnet.equiv_class(i,1);
    pvals = seq(ps);
    seq{i,t} = sample_node(bnet.CPD{e}, pvals);
    %fprintf('sample i=%d,t=%d,val=%d,ps\n', i, t, seq(i,t)); pvals(:)'
  end
end
t = 2;
done = 0;
while ~done
  for i=1:ss
    if ~isempty(stop_test) | isempty(seq{i,t})
      ps = parents(bnet.dag, i+ss) + (t-2)*ss;
      e = bnet.equiv_class(i,2);
      pvals = seq(ps);
      seq{i,t} = sample_node(bnet.CPD{e}, pvals);
      %fprintf('sample i=%d,t=%d,val=%d,ps\n', i, t, seq(i,t)); pvals(:)'
    end
  end
  if ~isempty(stop_test)
    done = feval(stop_test, seq(:,t));
  else
    if t==T
      done = 1;
    end
  end
  t = t + 1;
end


text:
function y = cvxdual( prob, name )
error( nargchk( 2, 2, nargin ) );

%
% Check problem
%

if ~isa( prob, 'cvxprob' ),
    error( 'First argument must be a cvxprob object.' );
end

%
% Create the object
%

y = class( struct( 'problem_', prob, 'name_', name ), 'cvxdual', cvxobj );

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


text:
%%%
%
% For attention-based models, given:
%   grad_contextVecs: lstmSize * batchSize
%   alignWeights: numPositions * batchSize.
%   srcHidVecs: lstmSize * batchSize * numPositions.
% compute the following grads:
%   grad_srcHidVecs: lstmSize * batchSize * numPositions
%   grad_alignWeights: numPositions * batchSize
%
% Thang Luong @ 2015, <lmthang@stanford.edu>
%
%%%
function [grad_alignWeights, grad_srcHidVecs] = contextLayerBackprop(grad_contextVecs, alignWeights, srcHidVecs, unmaskedIds, params)
  % change from grad_contextVecs lstmSize*batchSize -> lstmSize*batchSize*1
  grad_contextVecs = permute(grad_contextVecs, [1, 2, 3]);

  %% Grad formulae:
  %   contextVecs = H_src* a_t
  %   grad_srcHidVecs: outGrad * alignWeights'
  %   grad_alignWeights = H_src' * contexGrad (per example, to scale over multiple examples, i.e., batchSize, need to use bsxfun)

  %% grad_srcHidVecs
  grad_srcHidVecs = zeroMatrix(size(srcHidVecs), params.isGPU, params.dataType);
  grad_srcHidVecs(:, unmaskedIds, :) = bsxfun(@times, grad_contextVecs(:, unmaskedIds), permute(alignWeights(:, unmaskedIds), [3, 2, 1])); % change alignWeights -> 1 * batchSize * numPositions

  %% grad_alignWeights
  grad_alignWeights = zeroMatrix(size(alignWeights), params.isGPU, params.dataType);
  if size(srcHidVecs, 3)==1
    grad_alignWeights(:, unmaskedIds) = sum(bsxfun(@times, srcHidVecs(:, unmaskedIds, :), grad_contextVecs(:, unmaskedIds)), 1); % sum across lstmSize: numPositions * batchSize
  else
    grad_alignWeights(:, unmaskedIds) = squeeze(sum(bsxfun(@times, srcHidVecs(:, unmaskedIds, :), grad_contextVecs(:, unmaskedIds)), 1))'; % sum across lstmSize: numPositions * batchSize
  end

  %% assert
  if params.assert
    assert(isequal(size(grad_alignWeights), size(alignWeights)));

    % compute grad_srcHidVec in a different way
    grad_srcHidVecs1 = zeroMatrix(size(grad_srcHidVecs), params.isGPU, params.dataType);
    for jj=1:length(unmaskedIds)
      ii = unmaskedIds(jj);
      grad_srcHidVecs1(:, ii, :) = grad_contextVecs(:, ii, 1) * alignWeights(:, ii)';
    end
    assert(computeSum(grad_srcHidVecs1-grad_srcHidVecs, params.isGPU)==0);
  end
end

author: lmthang


text:
function x = tril( x, k )

%   Disciplined convex/geometric programming information for TRIL:
%       TRIL imposes no convexity restrictions on its arguments.

%
% Check inputs
%

if nargin < 2, k = 0; end
s = x.size_;
if length( s ) > 2,
    cvx_throw( 'The first argument must be 2-D.' );
elseif ~isnumeric( k ) || length( k ) ~= 1,
    cvx_throw( 'The second argument must be an integer scalar.' );
end

%
% Zero out the elements outside of the desired triangle
%

b = x.basis_;
b( :, ~tril(ones(s),k) ) = 0;
x = cvx( s, b );

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


text:
function varargout = unique(varargin)
% - Adds support for 'stable' flag, with all three outputs. If 'stable' or 'sorted'
% are used they must be the last input
% - The first ocurrence is returned by default, as in new Matlab versions. Octave,
% like old Matlab version, returns the last y default.
% - Second and third outputs are column vectors, as in new Matlab versions.
if nargin>=2 && strcmp(varargin{end},'stable')
    if strcmp(varargin{2},'rows')
        [a, ~, x] = builtin('unique', varargin{1}, 'rows', 'first');
    else
        [a, ~, x] = builtin('unique', varargin{1}, 'first');
    end
    if ~isempty(x)
        [~, ind] = max(bsxfun(@eq, x(:), (1:max(x)))); [y, ind2] = sort(ind(:));
        [~, ind3] = sort(ind2); z =  ind3(x);
        for n = 1:numel(x)
            x(n) = x(n) * ~any(x(n)==x(1:n-1));
        end
        x = nonzeros(x).';
    end
    if strcmp(varargin{2},'rows')
        x = a(x,:);
    else
        x = a(x);
    end
    varargout{1} = x;
    if nargout>=2, varargout{2} = y; end
    if nargout>=3, varargout{3} = z; end
else
    if strcmp(varargin{end},'sorted'), varargin(end) = []; end
    varargout = cell(1,nargout);
    if ~any(cellfun(@(x) any(strcmp(x, {'first', 'last'})), varargin(2:end))), varargin{end+1} = 'first'; end
    [varargout{:}] = builtin('unique', varargin{:});
end
if nargout>=2, varargout{2} = varargout{2}(:); end
if nargout>=3, varargout{3} = varargout{3}(:); end
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
function [mask]= mrAnatClassifyCleanMask(mask,minClusterSize)
% mrAnatClassifyCleanMask - remove small clusters
% [mask]= mrAnatClassifyCleanMask(mask,minClusterSize);

% 16-Jun-2005 SOD wrote it

if nargin < 1, help(mfilename); return; end;

if ieNotDefined('minClusterSize'), minClusterSize = 9; end;

% keep input
orgmask = mask;

% clean mask forground (1)
[LabMask remove]=myClipCluster(orgmask,minClusterSize);
fprintf('[%s]:Removing clusters (1):',mfilename);
for n=1:length(remove),
    mask(LabMask==remove(n))=0;
    fprintf('.');drawnow;
end
fprintf('Done.\n');drawnow;

% clean mask background (0)
[LabMask remove]=myClipCluster(-1.*orgmask+1,minClusterSize);
fprintf('[%s]:Removing clusters (0):',mfilename);
for n=1:length(remove),
    mask(LabMask==remove(n))=1;
    fprintf('.');drawnow;
end
fprintf('Done.\n');drawnow;

%-------------------------
function [imgLabel, remove]=myClipCluster(m,t)

% connectivity
%        6     three-dimensional six-connected neighborhood
%        18    three-dimensional 18-connected neighborhood
%        26    three-dimensional 26-connected neighborhood
conn = 6;

% get clusters
[imgLabel,numObjects] = bwlabeln(m, 6);

% get sizes
[imgHist,labelNum] = hist(imgLabel(:),0:numObjects);

% clean mask
remove = labelNum(find(imgHist(2:end)<t))+1;

return;

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [Population,W] = updateWeight(Population,W,Z,EP,nus)
% Delete overcrowded subproblems and add new subproblems

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

% This function is written by Lucas Farias

	%% Parameter setting
    [N,M] = size(Population.objs);

    %% Delete the overcrowded subproblems
    Dis = pdist2(Population.objs,Population.objs);
    Dis(logical(eye(length(Dis)))) = inf;
    Del = false(1,length(Population));
    while sum(Del) < min(nus,length(EP))
        Remain = find(~Del);
        subDis = sort(Dis(Remain,Remain),2);
        [~,worst] = min(prod(subDis(:,1:min(M,length(Remain))),2));
        Del(Remain(worst)) = true;
    end
    Population = Population(~Del);
    W = W(~Del,:);

    %% Add new subproblems
    % Determine the new solutions be added
    Combine  = [Population,EP];
    Selected = false(1,length(Combine));
    Selected(1:length(Population)) = true;
    Dis = pdist2(Combine.objs,Combine.objs);
    Dis(logical(eye(length(Dis)))) = inf;
    while sum(Selected) < min(N,length(Selected))
        subDis = sort(Dis(~Selected,Selected),2);
        [~,best] = max(prod(subDis(:,1:min(M,size(subDis,2))),2));
        Remain = find(~Selected);
        Selected(Remain(best)) = true;
    end
    % Add new subproblems
    newObjs = EP(Selected(length(Population)+1:end)).objs;
    temp    = 1./(newObjs-repmat(Z,size(newObjs,1),1));
    temp(temp==inf) = 0.999999; % when (temp == Z) then 0
    W = [W;temp./repmat(sum(temp,2),1,size(temp,2))];
    % Add new solutions
    Population = Combine(Selected);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO