text:
function result = TS_GetFromData(dataSource,dataField)
% TS_GetFromData   Load a given field from either a data file, or data structure
%                   (loaded from file).
%
%---INPUTS:
% dataSource: either a .mat filename or a structure generated from loading an
%               HCTSA file, e.g.,
%               dataSource = load('HCTSA.mat');
% dataField: the variable to extract
%
%---OUTPUTS:
% result: the variable extracted from either the file or structure (depending
%           on the type of dataSource)
%
%---EXAMPLE USAGE:
% fromDatabase = TS_GetFromData('HCTSA_EEGDataset.mat','fromDatabase');
% MasterOperations = TS_GetFromData('HCTSA_EEGDataset.mat','MasterOperations');

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

%-------------------------------------------------------------------------------
% In some cases, you provide a structure with the pre-loaded data already in it
% e.g., as a whatDataFile = load('HCTSA.mat');
%-------------------------------------------------------------------------------
if isstruct(dataSource)
    if isfield(dataSource,dataField)
        result = dataSource.(dataField);
    else
        result = [];
    end

%-------------------------------------------------------------------------------
% Often you provide a .mat file name: load in from file
%-------------------------------------------------------------------------------
elseif ischar(dataSource)

    % Decode shorthands for default filenames 'raw', 'loc', 'norm', and 'cl':
    switch dataSource
    case {'raw','loc'} % the raw, un-normalized data:
        dataSource = 'HCTSA.mat';
    case 'norm' % the normalized data:
        dataSource = 'HCTSA_N.mat';
    case 'cl' % the clustered data:
        dataSource = 'HCTSA_N.mat';
    end

    try
        loadAgain = load(dataSource,dataField);
        result = loadAgain.(dataField);
    catch
        result = [];
    end
    % fileVarsStruct = whos('-file',dataSource);
    % fileVars = {fileVarsStruct.name};
    % if ismember(dataField,fileVars)
    %     loadAgain = load(dataSource,dataField);
    %     result = loadAgain.(dataField);
    % else
    %     result = [];
    % end
else
    error('Unknown data source type: %s',class(dataSource));
end

% Check for legacy structure array format and convert to table:
if ismember(dataField,{'TimeSeries','Operations','MasterOperations'}) && isstruct(result)
    warning('Loaded metadata, %s, is still in the old structure-array format; converted to table',dataField)
    result = struct2table(result);
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function y = log(x)
%LOG          Implements  log(x)  for intervals
%
%   y = log(x)
%
%interval standard function implementation
%

% written  12/30/98     S.M. Rump
% modified 08/31/99     S.M. Rump  complex allowed, following N.C. Boersken:
%                                  Komplexe Kreis-Standardfunktionen,
%                                  Freiburger Intervallberichte 78/2,
%                                  NaN input, sparse input, log(0),
%                                  major revision, improved accuracy
% modified 12/06/99                branch cut with warning
% modified 09/02/00     S.M. Rump  rounding unchanged after use
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    accelaration for sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 12/04/05     S.M. Rump  'realstdfctsexcptnignore' added and some
%                                     improvements, tocmplx replaced by cintval,
%                                     extreme values for approximate part
% modified 09/06/07     S.M. Rump  approximate std fcts removed
% modified 10/23/07     S.M. Rump  complex numbers
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 10/20/08     S.M. Rump  check for zero omitted, improved performance
% modified 10/18/08     S.M. Rump  StdFctsException ignore/NaN
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/13/12     S.M. Rump  INTLAB_INTVAL_STDFCTS
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if x.complex
    if issparse(x.mid)
      x.mid = full(x.mid);
      x.rad = full(x.rad);
    end
    y = x;

%   y.mid = log(x.mid);
%   y.rad = - log( 1 - x.rad./abs(x.mid) );

    INTLAB_STDFCTS_PI = getappdata(0,'INTLAB_STDFCTS_PI');

    xmidre = real(x.mid);
    xmidim = imag(x.mid);
    Xmidim = intval(xmidim);
    Mim = atan(Xmidim./xmidre);      % -pi/2 <= Mim <= pi/2

    % special treatment of imaginary axis
    index = ( xmidre==0 );
    if any(index(:))
      Mim.inf(index) = INTLAB_STDFCTS_PI.PI2INF;          % pi/2
      Mim.sup(index) = INTLAB_STDFCTS_PI.PI2SUP;
      indexneg = index & ( xmidim<0 );
      if any(indexneg(:))
        Mim.inf(indexneg) = -INTLAB_STDFCTS_PI.PI2SUP;    % -pi/2
        Mim.sup(indexneg) = -INTLAB_STDFCTS_PI.PI2INF;
      end
    end

    index = ( xmidre < 0 );
    if any(index(:))
      % correct to atan2:  -pi <= Mim <= pi
      Pi = infsup( INTLAB_STDFCTS_PI.PIINF, ...
                   INTLAB_STDFCTS_PI.PISUP );
      signxmidim = sign(xmidim(index));
      corr = Pi.*(signxmidim+(signxmidim==0));
      Mim.inf(index) = Mim.inf(index) + corr.inf;
      Mim.sup(index) = Mim.sup(index) + corr.sup;
    end

    % x.mid = r*exp(j*phi)  ==>  log(x.mid) = log(r) + j*phi
    R = abs(x.mid);
    Mreinf = log_rnd(R,-1);
    Mresup = log_rnd(R,1);

    setround(1)
    mre = Mreinf + 0.5*(Mresup-Mreinf);
    mim = Mim.inf + 0.5*(Mim.sup-Mim.inf);
    y.mid = mre + j*mim;
    mrad = abs( mre-Mreinf + j*(mim-Mim.inf) );
    % log(R)  in  mre + j*mim +/- mrad

    wng = warning;
    warning off
    % x.rad < R ,  otherwise zero interval
    y.rad = log_rnd( -( x.rad./R - 1 ) , -1);
    warning(wng);
    y.rad = -y.rad + mrad;
    setround(0)

    % special treatment of branch cut
    index0 = ( R <= x.rad );               % zero interval
    index = ( ~index0 ) & ...
            (   ( ( xmidre<0 ) & ( xmidim>=0 ) & ( x.rad>xmidim ) ) ...
              | ( ( xmidre<0 ) & ( xmidim<0 ) & ( x.rad>=-xmidim ) ) ...
            );
    if any(index(:))
      warning('Complex Log: Input interval intersects with branch cut')
    end

    if any(index0(:))
      y.mid(index0) = complex(NaN,NaN);
      y.rad(index0) = NaN;
    end

    if rndold
      setround(rndold)
    end

    return
  end

  if issparse(x.inf)
    x.inf = full(x.inf);
    x.sup = full(x.sup);
  end
  % input x real and full
  % real range of definition:  [0,inf]
  INTLAB_STDFCTS_EXCPTN = getappdata(0,'INTLAB_STDFCTS_EXCPTN');
  index = ( x.inf<0 );                  % (partially) exceptional indices
  if any(index(:))                      % handle input out-of-range
    if INTLAB_STDFCTS_EXCPTN<=1  % out-of-range input handled as complex
      if INTLAB_STDFCTS_EXCPTN==1
        warning('LOG: Real interval input out of range changed to be complex')
      end
      y = x;
      %VVVV  y(index) = log(cintval(x(index)));
      s.type = '()'; s.subs = {index}; y = subsasgn(y,s,log(cintval(subsref(x,s))));
      %AAAA  Matlab bug fix
      index = ~index;
      if any(index(:))
        %VVVV  y(index) = log(x(index));
        s.type = '()'; s.subs = {index}; y = subsasgn(y,s,log(subsref(x,s)));
        %AAAA  Matlab bug fix
      end
      if rndold
        setround(rndold)
      end
      return
    end
    setappdata(0,'INTLAB_STDFCTS_EXCPTN_',1);
    if INTLAB_STDFCTS_EXCPTN==3    % ignore input out of range (ignore-mode)
      x.inf(index) = realmin*eps;;        % completely exceptional indices treated below
      indexneg = index & ( x.sup<0);      % completely exceptional indices
    end
  else
    indexneg = [];                        % make sure indexneg is not undefined
  end

  % input x real and full
  y = x;
  wng = warning;
  warning off

  % treat non-exceptional cases
  y.inf = log_rnd(x.inf,-1);
  y.sup = log_rnd(x.sup,1);

  if INTLAB_STDFCTS_EXCPTN==3      % ignore input out of range (ignore-mode)
    if ~isempty(find(indexneg))           % completely exceptional arguments to NaN
      y.inf(indexneg) = NaN;
      y.sup(indexneg) = NaN;
    end
  else                                    % any input out of range to NaN (NaN-mode)
    if ~isempty(find(index))              % exceptional arguments to NaN
      y.inf(index) = NaN;
      y.sup(index) = NaN;
    end
  end

  warning(wng)
  setround(rndold)

repo: OpenMAS


text:
function results = startSimulation(t0,tf,initialState,input_density,startParameters)
%	startSimulation  Starts the simulation of the Li-ion battery.
%
%   results = startSimulation(t0,tf,initialStates,input_density,startParameters)
%   Starts the simulation of the Li-ion Cell with LIONSIMBA.
%
%   Input:
%       - t0 : initial integration time
%       - tf : final integration time
%       - initialStates : structure containing data for initializing the
%       states of the battery.
%       - input_density  : Applied current/power density. If negative, the battery gets
%       discharged. If positive, the battery gets charged.
%       - startParameters : if provided, it defines the cell array containing the parameters
%       structures to be used in the simulation. Every single structure has to be obtained through the Parameters_init script.
%		If a cell array of N parameters structures is used, the simulator will perform a simulation with N cells in series.
%		If a cell of 1 parameters structure is used, then a single cell will be simulated.
%
%   Output:
%       - results : structure containing the solution of the dependent
%       variables. The stored results have as many rows as time instants
%       and columns as number of discrete volumes. If multiple cells are
%       simulated, the index i is used to access to the data of the i-th
%       cell.
%
%     results.Phis{i}:                      Solid Phase potential
%     results.Phie{i}:                      Liquid Phase potential
%     results.ce{i}:                        Electrolyte concentration
%     results.cs_surface{i}:                Electrode surface concentration
%     results.cs_average{i}:                Electrode average concentration
%     results.time{i}:                      Interpolated simulation time
%     results.int_internal_time{i}:         Integrator time steps
%     results.ionic_flux{i}:                Ionic flux
%     results.side_reaction_flux{i}:        Side reaction flux
%     results.SOC{i}:                       State of Charge
%     results.SOC_estimated{i}:             State of Charge estimate according to the
%                                           user-defined function
%     results.Voltage{i}:                   Cell voltage
%     results.Temperature{i}:               Cell Temperature
%     results.Qrev{i}:                      Reversible heat generation rate
%     results.Qrxn{i}:                      Reaction heat generation rate
%     results.Qohm{i}:                      Ohmic heat generation rate
%     results.film{i}:                      Side reaction film resistance
%     results.R_int{i}:                     Internal resistance
%     results.Up{i}:                        Cathode open circuit potential
%     results.Un{i}:                        Anode open circuit potential
%     results.etap{i}:                      Cathode overpotential
%     results.etan{i}:                      Anode overpotential
%     results.parameters{i}:                Parameters used for the simulation
%     results.JacobianFun:                  If evaluated, contains the
%                                           Jacobian matrix

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

try
    test_lionsimba_folder
catch
    error('It seems that you did not add to the Matlab path the battery_model_files directory and the folders therein. Please fix this problem and restart the simulation.')
end

% Version of LIONSIMBA
version       = '2.1';

if(isempty(startParameters))
    % Load battery's parameters if not provided by the user
    param{1} = Parameters_init;
else
    % Use user provided parameters
    param = startParameters;
end

% Validate the input current density/power_density values
if(param{1}.OperatingMode==1 || param{1}.OperatingMode==2)
    if(~isreal(input_density) || isnan(input_density) || isinf(input_density) || isempty(input_density))
        error('The input current/power densities provided by user is complex valued or is NaN. Please check the values and restart simulation.')
    end
end

% Check for environmental tool availability
checkEnvironment(param,nargin);

% If enabled, print the header information
if(param{1}.PrintHeaderInfo==1)
    headerInfo(version)
end

% If everything is ok, let's start to simulate.
results = mainCore(t0,tf,initialState,input_density,param);
end

function results = mainCore(t0,tf,initialState,input_density,param)

% Store the original parameters structure in order to return it at the
% end of simulations.
param_original  = param;

% Get the total number of cells that have to be simulated.
n_cells         = length(param);

% Check if more cells are simulated when potentiostatic conditions are required.
if(param{1}.OperatingMode==3 && n_cells~=1)
    clc;
    error('!!!ERROR!!! -- Potentiostatic simulations are only possible with single cell packs -- !!!ERROR!!!')
end

% Check if the initial state structure is given
[Y0_existence,YP0_existence,Y0,YP0] = checkInitialStates(initialState);

% Switch among the selected operating modes for defining a suitable getCurrentDensity/getPowerDensity
% function. If multiple cells are being simulated, the variable current/power profile
% or its constant value are retreived from the first element of the
% parameters structures. This is valid because in a series string, all cells carry the same current.
switch(param{1}.OperatingMode)
    case 1
        param{1}.getCurrentDensity = @(t,t0,tf,x,param,extra)input_density;
    case 2
        param{1}.getPowerDensity = @(t,t0,tf,x,param,extra)input_density;
    case 3
        param{1}.getCurrentDensity  = @(t,t0,tf,x,param,extra)0; % dummy values, will be over-written by solver
        param{1}.getPowerDensity    = @(t,t0,tf,x,param,extra)0; % unsure of this?
    case 4
        param{1}.getCurrentDensity = param{1}.CurrentDensityFunction; % external, i.e. user-supplied function values from external files loaded in parameterisation file
    case 5
        param{1}.getPowerDensity = param{1}.PowerDensityFunction; % external, i.e. user-supplied function values from external files loaded in parameterisation file
    otherwise
        error('Operating mode not supported');
end

% Define absolute and relative tolerances. If more cells are required,
% these values are taken from the first parameters structure.
opt.AbsTol      = param{1}.AbsTol;
opt.RelTol      = param{1}.RelTol;

n_diff          = zeros(n_cells,1);
n_alg           = zeros(n_cells,1);
start_x_index   = 1;
start_xp_index  = 1;

% For each cell, allcoate memory to store external functions used to
% estimate the SOC.
SOC_estimate    = cell(n_cells,1);

% Perform several checks over the cells
for i=1:n_cells
    % Check the daeFormulation flag in case startSimulation is called
    if(param{i}.daeFormulation~=1)
        error(['Make sure that the daeFormulation flag is set to 1 for each cell parameters structure. Cell ', num2str(i),' does not respect this constraint.'])
    end

    % When Fick's law of diffusion is used, at least 10 discretization
    % points are required. Raise an error if this condition is not met.
    if((param{i}.Nr_p<10 || param{i}.Nr_n<10) && param{i}.SolidPhaseDiffusion==3)
        error('The number of discrete points for the particles must be at least 10 in both cathode and anode.')
    end
    % Check if the SOC estimation function handle have been set. In case that
    % the function handle has not been defined or it does not have the right
    % number of input arguments, then return empty values.
    if(isempty(param{i}.SOC_estimation_function) || nargin(param{i}.SOC_estimation_function)~=6)
        SOC_estimate{i} = @(a,b,c,d,e,f,g,h,i,j,k)[];
    else
        SOC_estimate{i} = @SOCestimation;
    end
    param{i}.Nsum      = param{i}.Np + param{i}.Ns + param{i}.Nn;
    param{i}.Nsum_nos  = param{i}.Np + param{i}.Nn;

    % Define the discretization steps.
    param{i}.deltax_al     = 1 / param{i}.Nal;
    param{i}.deltax_p      = 1 / param{i}.Np;
    param{i}.deltax_s      = 1 / param{i}.Ns;
    param{i}.deltax_n      = 1 / param{i}.Nn;
    param{i}.deltax_cu     = 1 / param{i}.Ncu;

    % Compute the indices used to store the positions of the differential
    % and algebraic variables.
    param{i} = computeVariablesIndices(param{i});

    % Preallocate the differentiation matrices used for the solid phase
    % potential. This can be done here because such matrices are considered
    % to be time invariant.
    param{i} = solidPhaseDifferentiationMatrices(param{i});

    % Initialize Param.I_density or Param.P_density, using the value of the current density/power density (set in lines of code a few lines above
    if param{1}.OperatingMode==1 || param{1}.OperatingMode==4
        param{i}.I_density = 0;%param{1}.getCurrentDensity(0,t0,tf,x,param,param{1}.extraData);
    elseif param{1}.OperatingMode==2 || param{1}.OperatingMode==5
        param{i}.P_density = 0;%param{1}.getPowerDensity(0,t0,tf,x,param,param{1}.extraData);
    else
        param{i}.I_density = 0; % this is a dummy value
        param{i}.P_density = 0;
    end

    % Preallocate the differentiation matrices used for the solid phase
    % diffusion in case of Fick's law.
    param{i} = solidPhaseDiffusionDifferentiationMatrices(param{i});

    % Get the number of unknowns for the differential states
    [~, ~, ~, ~, ~, n_diff(i)] =   differentialInitialConditions(param{i});
    % Get the number of unknowns for the algebraic states
    [~, n_alg(i), ~] = initialise_model(param{i});

    % Store the number of differential and algebraic variables for each cell.
    param{i}.ndiff = n_diff(i);
    param{i}.nalg  = n_alg(i);

    % The x_index variable will be used in the battery model file
    % for indexing purposes.
    param{i}.x_index    = (start_x_index:n_diff(i)+n_alg(i)+start_x_index-1);

    param{i}.xp_index   = (start_xp_index:n_diff(i)+start_xp_index-1);

    % Update the starting x_index value for the (possible) next cell
    start_x_index       = n_diff(i)+n_alg(i)+start_x_index;

    start_xp_index      = n_diff(i)+n_alg(i)+start_xp_index;
end

for i=1:n_cells
    if((Y0_existence==0) && (YP0_existence==0))
        % Get the initial conditions for the differential states of the i-th cell
        [cs_average_init, ce_init, T_init, film_init, Q_init] =   differentialInitialConditions(param{i});
        % Solve the algebraic equations to find a set of semi-consistent initial
        % conditions for the algebraic equations. This will help the DAE solver as a warm startup.
        init_point = initialise_model(param{i});

        % Build the initial values array for the integrator
        Yt0 = [ce_init;cs_average_init;T_init;film_init;Q_init;init_point];
        Y0  = [Y0;Yt0];
        YP0 = [YP0;zeros(size(Yt0))];
    end
end

if(n_cells==1)
    nc = ' cell';
else
    nc = ' cells';
end

% Empty the used arrays
ce_t            = cell(n_cells,1);
cs_bar_t        = cell(n_cells,1);
T_t             = cell(n_cells,1);
jflux_t         = cell(n_cells,1);
Phis_t          = cell(n_cells,1);
Phie_t          = cell(n_cells,1);
cs_star_t       = cell(n_cells,1);
t_tot           = cell(n_cells,1);
Qrev_t          = cell(n_cells,1);
Qrxn_t          = cell(n_cells,1);
Qohm_t          = cell(n_cells,1);
SOC_t           = cell(n_cells,1);
Voltage_t       = cell(n_cells,1);
SOC_estimated_t = cell(n_cells,1);
film_t          = cell(n_cells,1);
js_t            = cell(n_cells,1);
R_int_t         = cell(n_cells,1);
curr_density_t  = cell(n_cells,1);
Up_t            = cell(n_cells,1);
Un_t            = cell(n_cells,1);
etap_t          = cell(n_cells,1);
etan_t          = cell(n_cells,1);
dudtp_t         = cell(n_cells,1);
dudtn_t         = cell(n_cells,1);
Q_t             = cell(n_cells,1);
yp_original     = YP0';

% This flag is used to notify the reason of the simulation stop. If 0
% everything went well.
exit_reason     = 0;

% Define the structure to be passed to the residual function
ida_user_data.param  = param;
ida_user_data.t0     = t0;
ida_user_data.tf     = tf;

% Define algebraic and differential variables.
% id:1-> differential variables,
% id:0-> algebraic variables.
id = [];
constraint_vector=[];
for i=1:n_cells
    id                                              = [id;ones(n_diff(i),1);zeros(n_alg(i),1)];
    temp_constraint_vector                          = zeros(n_diff(i)+n_alg(i),1);
    temp_constraint_vector(param{i}.Phis_indices)   = 1; % enforce positivity on solid phase potential in all nodes
    constraint_vector                               = [constraint_vector;temp_constraint_vector];
end

JacFun = [];

% This statement checks if the user wants to make use of the Jacobian
% matrix and (if yes) it has been already provided or not as part of the parameters structure.
if(param{1}.UseJacobian==1 && isempty(param{1}.JacobianFunction))
    % If the user wants to make use of the Jacobian, but it was not
    % provided in the parameters structure, then evaluate a new Jacobian.
    disp('Evaluating the analytical form of the Jacobian matrix. Please wait...')
    % Import casadi framework
    import casadi.*
    % Define the symbolic variables.
    xsym    = SX.sym('x',[sum(n_diff)+sum(n_alg),1]);
    xpsym   = SX.sym('xp',[sum(n_diff)+sum(n_alg),1]);
    cj      = SX.sym('cj',1);

    % Get the model equations written in an implicit form in a symbolic way.
    [dx_tot, ~, ~] = batteryModel(0,xsym,xpsym,ida_user_data);

    % Evaluate the Jacobian matrix. (Please refer to the Sundials guide for
    % further information about the Jacobian structure).
    J = jacobian(dx_tot,xsym) + cj*jacobian(dx_tot,xpsym);

    % Define a function for the Jacobian evaluation for a given set of
    % differential and algebraic variables.
    JacFun = Function('fJ',{xsym,cj},{J});

    % Store the function into a structure such that IDA will use it for the
    % evaluation of the Jacobian matrix (see jacobianFunction.m in simulator_tools).
    ida_user_data.fJ = JacFun;

    % Define the options for Sundials
    options = IDASetOptions('RelTol', opt.RelTol,...
        'AbsTol'        , opt.AbsTol,...
        'MaxNumSteps'   , 1500,...
        'VariableTypes' , id,...
        'UserData'      , ida_user_data,...
        'JacobianFn'    , @jacobianFunction,...
        'ConstraintTypes', constraint_vector);

elseif(param{1}.UseJacobian==1 && ~isempty(param{1}.JacobianFunction))
    % If the Jacobian needs to be used and it has also been provided in the
    % parameters structure, use it directly.

    disp('Analytical function of the Jacobian matrix provided by the user.')
    % If the Jacobian has been provided from the user, use it.
    JacFun = param{1}.JacobianFunction;

    % Pass this function pointer to the routine that IDA will call for the
    % evaluation of the Jacobian values.
    ida_user_data.fJ = JacFun;

    % Define the options for Sundials
    options = IDASetOptions('RelTol', opt.RelTol,...
        'AbsTol'        , opt.AbsTol,...
        'MaxNumSteps'   , 1500,...
        'VariableTypes' , id,...
        'UserData'      , ida_user_data,...
        'JacobianFn'    , @jacobianFunction,...
        'ConstraintTypes', constraint_vector);
else
    % In this case the user does not want to make use of the Jacobian
    % matrix. A numerical approximation will be calculated instead.

    % Define the options for Sundials
    options = IDASetOptions('RelTol', opt.RelTol,...
        'AbsTol'        , opt.AbsTol,...
        'MaxNumSteps'   , 1500,...
        'VariableTypes' , id,...
        'UserData'      , ida_user_data,...
        'ConstraintTypes', constraint_vector);
end

% Enable the capability to deal with piecewise inputs only when custom
% current profiles are used
if(param{1}.OperatingMode == 4)
    % Add to IDA options structure the pointer to the function used to
    % identify the presence of events (in this particular case
    % discontinuities in the applied input currents)
    options.RootsFn     = @rootFinder;
    options.NumRoots    = 1;
end

% Initialise solver
IDAInit(@batteryModel,t0,Y0,YP0,options);

if param{i}.Scope == 1
    disp(['Finding a set of consistent ICs for ',num2str(n_cells),nc,' battery pack. Please wait..'])
end
% Find consistent initial conditions
[~, yy, ~] = IDACalcIC(t0+10,'FindAlgebraic');

% Init the starting integration time
t = t0;

% Store in the results the initial states values.
y = yy';

% Store the total number of complete set of states
y_tot = y;

[ce_t,cs_bar_t,T_t,jflux_t,Phis_t, Phie_t, cs_star_t, SOC_t, film_t, js_t,Up_t,Un_t,R_int_t,curr_density_t,Voltage_t,SOC_estimated_t,Qrev_t,Qrxn_t,Qohm_t,Q_t,~,dudtp_t, dudtn_t,t_tot] =...
    storeSimulationResults(n_cells,ce_t,cs_bar_t,T_t,jflux_t,Phis_t, Phie_t, cs_star_t, SOC_t, film_t, js_t,curr_density_t,Voltage_t,SOC_estimated_t,Up_t,Un_t,R_int_t,Qrev_t,Qrxn_t,Qohm_t,Q_t,dudtp_t, dudtn_t, t_tot, y, t,SOC_estimate,t0,tf, param);

sim_time = 0;  % Simulation time (i.e. wall) for reporting purposes only. Not used in controlling solver or time-loop
% Loop until the integration time reaches tf.
while(t<tf)

    % Check if simulation stop conditions are met or not
    exit_reason = checkSimulationStopConditions(n_cells, Phis_t, cs_bar_t, T_t, param);
    % If a stop condition is reached, stop the simulation
    if(exit_reason~=0)
        break;
    end

    %% Solve the set of DAEs
    % The solver IDA is used to solve the resulting set of DAEs. Please
    % refer to IDA manual for more information about syntax and its usage.
    tic
    [status, t, y]   = IDASolve(tf,'OneStep');

    % If status > 0 it means that roots have been found during the
    % resolution of the equations
    if(status>0)
        % In case of custom current profiles, the roots are determined by
        % the presence of a discontinuity in the applied current density
        if(param{1}.OperatingMode==4)
            % Define a new time instant at which re-initialize the solver
            % using the last known values of the states
            t = t*(1+1e-5);
            IDAReInit(t,y,0*y,options);

            % Find consistent initial conditions starting from the new
            % point and keep on integrating
            [~, y, yp] = IDACalcIC(t+10,'FindAlgebraic');
        else
            error(e.message);
        end
    end
    sim_time    = sim_time+toc;
    y           = y';
        % Store the matrix containing all the states at all the integration
    % time steps.
    y_tot       = [y_tot;y];
    if(status==0)
        % Store derivative info at each time step
        yp_original = [yp_original;IDAGet('DerivSolution',t,1)'];
    else
        yp_original = [yp_original;yp'];
    end

    [ce_t,cs_bar_t,T_t,jflux_t,Phis_t, Phie_t, cs_star_t, SOC_t, film_t, js_t,Up_t,Un_t,R_int_t,curr_density_t,Voltage_t,SOC_estimated_t,Qrev_t,Qrxn_t,Qohm_t,Q_t,tot_voltage,dudtp_t, dudtn_t,t_tot] =...
        storeSimulationResults(n_cells,ce_t,cs_bar_t,T_t,jflux_t,Phis_t, Phie_t, cs_star_t, SOC_t, film_t, js_t,curr_density_t,Voltage_t,SOC_estimated_t,Up_t,Un_t,R_int_t,Qrev_t,Qrxn_t,Qohm_t,Q_t,dudtp_t, dudtn_t, t_tot, y, t,SOC_estimate,t0,tf, param);

    % If the output scope is active, show additional information to the user
    if(param{1}.Scope==1)
        if(n_cells==1)
            temperature     = T_t{1}(end,end);
            % If Fick's law of diffusion is used, before to evaluate the
            % SOC, it is necessary to compute the average solid
            % concentration in each particle.
            Sout = internalSOCestimate(cs_bar_t,param,1);
            clc
            fprintf(['No. of cells in the pack \t',num2str(n_cells),'\n']);
            fprintf(['Time \t\t\t\t\t',num2str(t),' s\n']);
            % If potentiostatic mode is running, applied current comes as
            % solution of DAEs. Otherwise it is provided by the user.
            if(param{1}.OperatingMode==1) || (param{1}.OperatingMode==4)
                fprintf(['Applied current density \t\t',num2str(y(end)),' A/m^2\n']);
            elseif (param{1}.OperatingMode==2) || (param{1}.OperatingMode==5)
                fprintf(['Applied power density \t\t',num2str(param{1}.getPowerDensity(t,t0,tf,param{1}.extraData)),' W/m^2\n']);
            end

			switch param{1}.edge_values
				% Linear interpolation
				case 2
					output_voltage = (1.5*Phis_t{1}(end,1)-0.5*Phis_t{1}(end,2)) - (1.5*Phis_t{1}(end,end) - 0.5*Phis_t{1}(end,end-1));
				% Value at the centroid
				otherwise
					output_voltage = Phis_t{1}(end,1) - Phis_t{1}(end,end);
			end
            fprintf(['Voltage \t\t\t\t',          num2str(output_voltage),   ' V\n']);
            fprintf(['Temperature \t\t\t',        num2str(temperature),                           ' K\n']);
            fprintf(['SOC \t\t\t\t\t',            num2str(Sout),                                  ' %% \n']);
            fprintf(['Cutoff Voltage \t\t\t',     num2str(param{1}.CutoffVoltage),                ' V\n']);
            fprintf(['Cutover Voltage \t\t',      num2str(param{1}.CutoverVoltage),               ' V\n']);
            fprintf(['Internal Resistance \t',    num2str(R_int_t{1}(end)),                       ' Ohm m^2\n']);
            fprintf(['Absolute tolerance \t\t',   num2str(param{1}.AbsTol),                       '\n']);
            fprintf(['Relative tolerance \t\t',   num2str(param{1}.RelTol),                       '\n']);
            fprintf(['Initial int. time \t\t',    num2str(t0),                                    ' s\n']);
            fprintf(['Final int. time \t\t',      num2str(tf),                                    ' s\n']);
            fprintf(['N. of unknowns \t\t\t',     num2str(length(y)),                             ' \n']);
        else
            clc
            fprintf(['No. of cells in the pack \t',num2str(n_cells),'\n']);
            fprintf(['Time \t\t\t\t\t',num2str(t),' s\n']);
            % If potentiostatic mode is running, applied current comes as
            % solution of DAEs. Otherwise it is provided by the user.
            if(param{1}.OperatingMode==1) || (param{1}.OperatingMode==4)
                fprintf(['Applied current density \t\t',num2str(y(end)),' A/m^2\n']);
            elseif (param{1}.OperatingMode==2) || (param{1}.OperatingMode==5)
                fprintf(['Applied power density \t\t',num2str(param{1}.getPowerDensity(t,t0,tf,param{1}.extraData)),' W/m^2\n']);
            end

            fprintf(['Voltage \t\t\t\t',          num2str(tot_voltage),       ' V\n']);
            fprintf(['Absolute tolerance \t\t',   num2str(param{1}.AbsTol),   ' \n']);
            fprintf(['Relative tolerance \t\t',   num2str(param{1}.RelTol),   ' \n']);
            fprintf(['Initial int. time \t\t',    num2str(t0),                ' s\n']);
            fprintf(['Final int. time \t\t',      num2str(tf),                ' s\n']);
            fprintf(['N. of unknowns \t\t\t',     num2str(length(y)),         ' \n']);
        end
    end
end

disp(['Elasped time: ',num2str(sim_time),' s']);

% Interpolate for fixed time step values
t_tot_original = t_tot;

% Build the time vector used for interpolation
time_vector = (t0:param{i}.sim_datalog_interval:tf);

% In case of the simulation has stopped before the final time set by the
% user, change the tf variable in order to interpolate only available values.

if(t<tf)
    % Set final time equal to the last integration step
    tf          = t;
    % Redefine the time vector used for interpolation
    time_vector = (t0:param{i}.sim_datalog_interval:tf);
end

% If at least one integration step has been done, retreive the first order
% time derivative information. Otherwise use the initial data.
if(time_vector(end)>t0)
    % Retreive derivative information at the last time step
    yp          = interp1(t_tot{i},yp_original,time_vector(end))';
else
    % If the integration step carried out by SUNDIALS is less than the
    % parametrized step size, then return the initial data as set of initial states.
    yp          = YP0;
    yp_original = YP0;
end

% Free memory allocated by IDA solver
IDAFree

% These variables will be used to store the original results of the integration process.
Phis_t_o            = cell(n_cells,1);
Phie_t_o            = cell(n_cells,1);
ce_t_o              = cell(n_cells,1);
cs_star_t_o         = cell(n_cells,1);
cs_average_t_o      = cell(n_cells,1);
jflux_t_o           = cell(n_cells,1);
SOC_t_o             = cell(n_cells,1);
T_t_o               = cell(n_cells,1);
Voltage_t_o         = cell(n_cells,1);
SOC_estimated_t_o   = cell(n_cells,1);
film_t_o            = cell(n_cells,1);
js_t_o              = cell(n_cells,1);
R_int_t_o           = cell(n_cells,1);
Up_t_o              = cell(n_cells,1);
Un_t_o              = cell(n_cells,1);
dudtp_t_o           = cell(n_cells,1);
dudtn_t_o           = cell(n_cells,1);
Qrev_t_o            = cell(n_cells,1);
Qrxn_t_o            = cell(n_cells,1);
Qohm_t_o            = cell(n_cells,1);
etap_t_o            = cell(n_cells,1);
etan_t_o            = cell(n_cells,1);
app_current_t_o     = cell(n_cells,1);
Q_t_o               = cell(n_cells,1);
y                   = [];
y_original          = [];
for i=1:n_cells
    % Save the overpotentials
    etap_t{i} = Phis_t{i}(:,1:param{i}.Np)-Phie_t{i}(:,1:param{i}.Np)-Up_t{i};
    if(param{i}.EnableAgeing==1)
        etan_t{i} = Phis_t{i}(:,param{i}.Np+1:end)-Phie_t{i}(:,param{i}.Np+param{i}.Ns+1:end)-Un_t{i} - param{1}.F*jflux_t{i}(:,param{i}.Np+1:end).*(param{i}.R_SEI+film_t{i}./param{i}.k_n_aging);
    else
        etan_t{i} = Phis_t{i}(:,param{i}.Np+1:end)-Phie_t{i}(:,param{i}.Np+param{i}.Ns+1:end)-Un_t{i};
    end
    if(param{i}.sim_datalog_interval>0)
        % Store original results
        Phis_t_o{i}            = Phis_t{i};
        Phie_t_o{i}            = Phie_t{i};
        ce_t_o{i}              = ce_t{i};
        cs_star_t_o{i}         = cs_star_t{i};
        cs_average_t_o{i}      = cs_bar_t{i};
        jflux_t_o{i}           = jflux_t{i};
        SOC_t_o{i}             = SOC_t{i};
        T_t_o{i}               = T_t{i};
        Voltage_t_o{i}         = Voltage_t{i};
        SOC_estimated_t_o{i}   = SOC_estimated_t{i};
        film_t_o{i}            = film_t{i};
        js_t_o{i}              = js_t{i};
        R_int_t_o{i}           = R_int_t{i};
        app_current_t_o{i}     = curr_density_t{i};
        Up_t_o{i}              = Up_t{i};
        Un_t_o{i}              = Un_t{i};
        dudtp_t_o{i}           = dudtp_t{i};
        dudtn_t_o{i}           = dudtn_t{i};
        etap_t_o{i}            = etap_t{i};
        etan_t_o{i}            = etan_t{i};
        Qrev_t_o{i}            = Qrev_t{i};
        Qrxn_t_o{i}            = Qrxn_t{i};
        Qohm_t_o{i}            = Qohm_t{i};
        Q_t_o{i}               = Q_t{i};

        if(time_vector(end)>t0)
            % Interpolate the results
            Phis_t{i}          = interp1(t_tot{i},Phis_t{i},time_vector');
            Phie_t{i}          = interp1(t_tot{i},Phie_t{i},time_vector');
            ce_t{i}            = interp1(t_tot{i},ce_t{i},time_vector');
            cs_star_t{i}       = interp1(t_tot{i},cs_star_t{i},time_vector');
            cs_bar_t{i}        = interp1(t_tot{i},cs_bar_t{i},time_vector');
            jflux_t{i}         = interp1(t_tot{i},jflux_t{i},time_vector');
            SOC_t{i}           = interp1(t_tot{i},SOC_t{i},time_vector');
            SOC_estimated_t{i} = interp1(t_tot{i},SOC_estimated_t{i},time_vector');
            Voltage_t{i}       = interp1(t_tot{i},Voltage_t{i},time_vector');
            film_t{i}          = interp1(t_tot{i},film_t{i},time_vector');
            js_t{i}            = interp1(t_tot{i},js_t{i},time_vector');
            R_int_t{i}         = interp1(t_tot{i},R_int_t{i},time_vector');
            T_t{i}             = interp1(t_tot{i},T_t{i},time_vector');
            curr_density_t{i}  = interp1(t_tot{i},curr_density_t{i},time_vector');
            Up_t{i}            = interp1(t_tot{i},Up_t{i},time_vector');
            Un_t{i}            = interp1(t_tot{i},Un_t{i},time_vector');
            Qrev_t{i}          = interp1(t_tot{i},Qrev_t{i},time_vector');
            Qrxn_t{i}          = interp1(t_tot{i},Qrxn_t{i},time_vector');
            Qohm_t{i}          = interp1(t_tot{i},Qohm_t{i},time_vector');
            etap_t{i}          = interp1(t_tot{i},etap_t{i},time_vector');
            etan_t{i}          = interp1(t_tot{i},etan_t{i},time_vector');
            dudtp_t{i}         = interp1(t_tot{i},dudtp_t{i},time_vector');
            dudtn_t{i}         = interp1(t_tot{i},dudtn_t{i},time_vector');
            Q_t{i}             = interp1(t_tot{i},Q_t{i},time_vector');
            t_tot{i}           = time_vector';
        end
    end
    % Store results. If integration steps are enabled, store the interpolated data.
    results.Phis{i}                        = Phis_t{i};
    results.Phie{i}                        = Phie_t{i};
    results.ce{i}                          = ce_t{i};
    results.cs_surface{i}                  = cs_star_t{i};
    results.cs_average{i}                  = cs_bar_t{i};
    results.time{i}                        = t_tot{i};
    results.int_internal_time{i}           = t_tot_original{i};
    results.ionic_flux{i}                  = jflux_t{i};
    results.side_reaction_flux{i}          = js_t{i};
    results.SOC{i}                         = SOC_t{i};
    results.SOC_estimated{i}               = SOC_estimated_t{i};
    results.Voltage{i}                     = Voltage_t{i};
    results.Temperature{i}                 = T_t{i};
    results.Qrev{i}                        = Qrev_t{i};
    results.Qrxn{i}                        = Qrxn_t{i};
    results.Qohm{i}                        = Qohm_t{i};
    results.film{i}                        = film_t{i};
    results.R_int{i}                       = R_int_t{i};
    results.Up{i}                          = Up_t{i};
    results.Un{i}                          = Un_t{i};
    results.etap{i}                        = etap_t{i};
    results.etan{i}                        = etan_t{i};
    results.dudtp{i}                       = dudtp_t{i};
    results.dudtn{i}                       = dudtn_t{i};
    results.Q{i}                           = Q_t{i};
    results.parameters{i}                  = param{i};
    results.JacobianFun                    = JacFun;

    % Store original data.
    results.original.Phis{i}               = Phis_t_o{i};
    results.original.Phie{i}               = Phie_t_o{i};
    results.original.ce{i}                 = ce_t_o{i};
    results.original.cs_surface{i}         = cs_star_t_o{i};
    results.original.cs_average{i}         = cs_average_t_o{i};
    results.original.ionic_flux{i}         = jflux_t_o{i};
    results.original.side_reaction_flux{i} = js_t_o{i};
    results.original.SOC{i}                = SOC_t_o{i};
    results.original.SOC_estimated{i}      = SOC_estimated_t_o{i};
    results.original.Voltage{i}            = Voltage_t_o{i};
    results.original.Temperature{i}        = T_t_o{i};
    results.original.film{i}               = film_t_o{i};
    results.original.R_int{i}              = R_int_t_o{i};
    results.original.Up{i}                 = Up_t_o{i};
    results.original.Un{i}                 = Un_t_o{i};
    results.original.etap{i}               = etap_t_o{i};
    results.original.etan{i}               = etan_t_o{i};
    results.original.Q{i}                  = Q_t_o{i};
    results.original.parameters{i}         = param_original{i};

    % Store initial states data
    y           = [y;ce_t{i}(end,:)';cs_bar_t{i}(end,:)';T_t{i}(end,:)';film_t{i}(end,:)';Q_t{i}(end,:)';jflux_t{i}(end,:)';Phis_t{i}(end,:)';Phie_t{i}(end,:)';js_t{i}(end,:)';curr_density_t{i}(end)];
    % Store initial states original data
    y_original  = [y_original;ce_t_o{i}(end,:)';cs_average_t_o{i}(end,:)';T_t_o{i}(end,:)';film_t_o{i}(end,:)';Q_t_o{i}(end,:)';jflux_t_o{i}(end,:)';Phis_t_o{i}(end,:)';Phie_t_o{i}(end,:)';js_t_o{i}(end,:)';app_current_t_o{i}(end)];
end

% Store the array of last results
results.Y                           = y;
results.YP                          = yp;

results.original.Y                  = y_tot;
results.original.YP                 = yp_original;

results.original.initialState.Y     = y_original;
results.original.initialState.YP    = yp_original(end,:);

results.initialState.Y              = y;
results.initialState.YP             = yp;

% Store simulation time
results.simulation_time             = sim_time;

% Exit reason
results.exit_reason                 = exit_reason;

% Check the operating mode and store the results accordingly.
if(param{1}.OperatingMode==2)
    results.power_density          = param{1}.P_density * ones(size(t_tot{1},1),1);
    results.original.power_density = param{1}.P_density * ones(size(t_tot_original{1},1),1);
    % Variable current profile
elseif(param{1}.OperatingMode==5)
    results.power_density          = param{1}.getPowerDensity(t_tot{1},t0,tf,param{1}.extraData);
    results.original.power_density = param{1}.getPowerDensity(t_tot_original{1},t0,tf,param{1}.extraData);
end

results.curr_density          = curr_density_t{1};
results.original.curr_density = app_current_t_o{1};
end

function estimate = SOCestimation(t,t0,tf,param,ce_t,cs_bar_t,cs_star_t,Phie_t,Phis_t,jflux_t,T_t)
% Build the states struct which will be passed to the function
states.ce           = ce_t(end,:);
states.cs_average   = cs_bar_t(end,:);
states.cs_surface   = cs_star_t(end,:);
states.Phie         = Phie_t(end,:);
states.Phis         = Phis_t(end,:);
states.ionic_flux   = jflux_t(end,:);
states.Temperature  = T_t(end,:);
% Call the estimation procedure
estimate = param.SOC_estimation_function(t,t0,tf,states,param.extraData,param);
end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
function column_num = file_column_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_COLUMN_COUNT counts the columns in the first line of a file.
%
%  Discussion:
%
%    The file is assumed to be a simple text file.
%
%    Most lines of the file are presumed to consist of COLUMN_NUM words,
%    separated by spaces.  There may also be some blank lines, and some
%    comment lines, which have a "#" in column 1.
%
%    The routine tries to find the first non-comment non-blank line and
%    counts the number of words in that line.
%
%    If all lines are blanks or comments, it goes back and tries to analyze
%    a comment line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the file.
%
%    Output, integer COLUMN_NUM, the number of columns in the file.
%
  FALSE = 0;
  TRUE = 1;
%
%  Open the file.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_COLUMN_COUNT - Error!' );
    column_num = -1;
    return;
  end
%
%  Read one line, but skip blank lines and comment lines.
%  Use FGETL so we drop the newline character!
%
  got_one = FALSE;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

    elseif ( line(1) == '#' )

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

  if ( got_one == FALSE )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Warning!\n' );
    fprintf ( 1, '  The file does not seem to contain any data.\n' );
    column_num = -1;
    return;
  end

  column_num = s_word_count ( line );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function value = i4vec_multinomial_pdf ( n, p, m, x )

%*****************************************************************************80
%
%% I4VEC_MULTINOMIAL_PDF evaluates the multinomial PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2013
%
%  Author:
%
%    John Burkardt.
%
%  Parameters:
%
%    Input, integer N, the number of trials.
%
%    Input, real P(M), the probability of each outcome
%    on any single trial.
%
%    Input, integer M, the number of possible outcomes
%    of a single trial.
%
%    Input, integer X(M), the results of N trials,
%    with X(I) the number of times outcome I occurred.
%
%    Output, real VALUE, the probability
%    density function evaluated at X.
%

%
%  The combinatorial coefficient is an integer.
%
  c = 1;
  top = n;
  for i = 1 : m
    bot = 1;
    do j = 1 : x(i)
      c = ( c * top ) / bot;
      top = top - 1;
      bot = bot + 1;
    end
  end

  value = c;
  for i = 1 : m
    value = value * p(i) ^ x(i);
  end

  return
end


text:
%ReactionProcessTestCase
% Base class for whole cell reaction process test classes. Currently does
% not provide anything beyond the ProcessTestCase class.
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 7/13/2010
classdef ReactionProcessTestCase < edu.stanford.covert.cell.sim.ProcessTestCase
    %constructor
    methods
        function this = ReactionProcessTestCase(name)
            this = this@edu.stanford.covert.cell.sim.ProcessTestCase(name);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
classdef CatalogExplorationPlot < handle
    % CATALOGEXPLORATIONPLOT create a plot where x,y,z,color, and size are modifiable
    properties
        x_by (1,:)      char
        y_by (1,:)      char
        z_by (1,:)      char
        color_by (1,:)  char            = ZmapGlobal.Data.mainmap_plotby%'Date'
        size_by (1,:)   char            ='Magnitude'
        colorFcn        function_handle = @datenum
        sizeFcn         function_handle = @mag2dotsize
        catalogFcn      function_handle;
        axes_choices    cell            = {};
        myscatter;
        ax              matlab.graphics.axis.Axes;
        conversions;
        curview;
        marker                          = ZmapGlobal.Data.event_marker;
    end
    methods

        function obj=CatalogExplorationPlot(ax, catalogFcn)
            % creates a plot with arbitrarily modifiable axes
            % obj=CATALOGEXPLORATIONPLOT(ax, catalogFcn)
            mycat = catalogFcn();
            obj.x_by = mycat.XLabel;
            obj.y_by = mycat.YLabel;
            obj.z_by = mycat.ZLabel;
            obj.catalogFcn = catalogFcn;
            obj.set_valid_axes_choices();
            obj.set_conversions();
            obj.ax=ax;
            c=ax.UIContextMenu;
            if isempty(c)
                c=uicontextmenu(ancestor(ax,'figure'),'Tag','catexplot');
                addLegendToggleContextMenuItem(c,'top','below');
                ax.UIContextMenu=c;
            else
                addLegendToggleContextMenuItem(c,'top','below');
            end
        end

        function scatter(obj, tag, varargin)
            % scatter plot with interactive axes
            c=obj.catalogFcn();
            x=c.(obj.x_by);
            y=c.(obj.y_by);
            z=c.(obj.z_by);
            s=c.(obj.size_by);
            if ~iscategorical(s)
                s=obj.sizeFcn(s);
            end
            if obj.color_by == "-none-"
                cl=[0 0 0];
            else
                cl=c.(obj.color_by);
                if ~iscategorical(cl)
                    cl=obj.colorFcn(cl);
                end
            end
            % delete(findobj(obj.ax,'Tag',tag));
            if isempty(obj.myscatter)
                obj.myscatter=scatter3(obj.ax,x, y, z, s, cl,'Marker',obj.marker,'Tag',tag);
                obj.myscatter.DisplayName=sprintf('size:%s\ncolor:%s',obj.size_by,obj.color_by);
                grid(obj.ax,'on');
                box(obj.ax,'on');
                fig=ancestor(obj.ax,'figure');
                xl = xlabel(obj.x_by,'interpreter','none');
                yl = ylabel(obj.y_by,'interpreter','none');
                zl = zlabel(obj.z_by,'interpreter','none');
                obj.xContextMenu(xl, tag, fig);
                obj.yContextMenu(yl, tag, fig);
                obj.zContextMenu(zl, tag, fig);
                obj.scatterContextMenu(obj.myscatter, tag);
                obj.ax.UserData.cep = obj;
            else
                obj.ax.NextPlot='replace';
                obj.myscatter.XData=x;
                obj.myscatter.YData=y;
                obj.myscatter.ZData=z;
                obj.myscatter.SizeData=s;
                obj.myscatter.CData=cl;
                obj.myscatter.DisplayName=sprintf('size:%s\ncolor:%s',obj.size_by,obj.color_by);
            end
            if isempty(obj.curview)
                view(obj.ax,2);
            else
                view(obj.ax,obj.curview);
            end

        end

        function update(obj, specific)
            % UPDATE updates the scatter plot, optionally changing only one axis (or color or size)
            %
            % obj.update() update all aspects of the scatter plot
            %
            % obj.update( SPECIFIC ) updates only thes specific part of the plot, where SPECIFIC can
            % be 'x_by', 'y_by', 'z_by', 'size_by', 'sin
            %
            c=obj.catalogFcn();
            [obj.curview(1), obj.curview(2)] = view(obj.ax);
            if ~exist('specific','var')
                switch obj.color_by
                    case '-none-'
                        cdata=obj.colorFcn(1);
                        if ~isequal(size(cdata),[1,3])
                            cdata = [0 0 0];
                        end
                    case 'MagnitudeType'
                        cdata= c.(obj.color_by);
                        set(obj.myscatter,'CData',cdata);
                        cc = categories(cdata);
                        ncats = numel(cc);
                        colormap(obj.ax, lines(ncats) );
                        centers = linspace(1.5, ncats - 0.5, ncats);
                        colorbar('Limits',[1,ncats],...
                            'Ytick', centers,...
                            'YtickLabel', cc)
                    otherwise
                        cdata= obj.colorFcn(c.(obj.color_by));

                end
                switch obj.size_by
                    case 'Single Size'
                        sdata=obj.sizeFcn(1);
                    case 'MagnitudeType'
                        sdata=c.(obj.size_by);
                    otherwise
                        sdata=obj.sizeFcn(c.(obj.size_by));
                end
                set( obj.myscatter,...
                    'XData', c.(obj.x_by),...
                    'YData', c.(obj.y_by),...
                    'ZData', c.(obj.z_by),...
                    'SizeData',sdata,...
                    'CData',cdata...
                    );
            else
                switch specific
                    case 'x_by'
                        doit('XAxis','XData', obj.x_by);
                    case 'y_by'
                        doit('YAxis','YData', obj.y_by);
                    case 'z_by'
                        doit('ZAxis','ZData', obj.z_by);
                    case 'size_by'
                        switch obj.size_by
                            case 'Single Size'
                                set(obj.myscatter,'SizeData',obj.sizeFcn(1));
                    case 'MagnitudeType'
                                set(obj.myscatter,'SizeData',c.(obj.size_by));
                            otherwise
                                set(obj.myscatter,'SizeData',obj.sizeFcn(c.(obj.size_by)));
                        end
                    case 'color_by'
                        switch obj.color_by
                            case '-none-'
                                set(obj.myscatter,'CData',obj.colorFcn(1));
                            case 'MagnitudeType'
                                cdata= c.(obj.color_by);
                                set(obj.myscatter,'CData',cdata);
                                cc = categories(cdata);
                                ncats = numel(cc);
                                colormap(obj.ax, lines(ncats) );
                                centers = linspace(1.5, ncats - 0.5, ncats);
                                colorbar('Limits',[1,ncats],...
                                    'Ytick', centers,...
                                    'YtickLabel', cc)
                            otherwise
                                set(obj.myscatter,'CData',obj.colorFcn(c.(obj.color_by)));
                        end

                end
            end
            %{
            function val = use_correct_ruler(obj, axisName, val)
                switch class(val)
                    case 'datetime'
                        if ~isa(obj.ax.(axisName),'matlab.graphics.axis.decorator.DatetimeRuler')
                            set(obj.ax, axisName, matlab.graphics.axis.decorator.DatetimeRuler);
                        end
                    case 'duration'
                        if ~isa(obj.ax.(axisName),'matlab.graphics.axis.decorator.DurationRuler')
                            set(obj.ax, axisName, matlab.graphics.axis.decorator.DurationRuler);
                        end
                    case 'categorical'
                        if ~isa(obj.ax.(axisName),'matlab.graphics.axis.decorator.CategoricalRuler')
                            set(obj.ax, axisName, matlab.graphics.axis.decorator.CategoricalRuler);
                        end
                    otherwise % numeric
                        if ~isa(obj.ax.(axisName),'matlab.graphics.axis.decorator.NumericRuler')
                            set(obj.ax, axisName, matlab.graphics.axis.decorator.NumericRuler);
                        end
                end
            end
            %}
            function doit(axAx, where, fld)
                % doit(axAx, where, fld) poor name.
                % axAx: 'XAxis', etc...
                % where: 'XData', etc...
                % fld: 'Longitude', etc... which is the result of obj.x_by, etc...
                % axisH = obj.ax.(axAx);
                    enforce_linear_scale_if_necessary();
                    try
                        obj.myscatter.(where) = c.(fld);
                        axruler=obj.ax.(axAx);
                        assert(...
                            (isnumeric(c.fld) &&  isa(axruler,'matlab.graphics.axis.decorator.NumericRuler')) || ...
                            (isdatetime(c.fld) && isa(axruler,'matlab.graphics.axis.decorator.DatetimeRuler')) || ...
                            (iscategorical(c.fld) && isa(axruler,'matlab.graphics.axis.decorator.CategoricalRuler')) ||...
                            (isduration(c.fld) && isa(axruler,'matlab.graphics.axis.decorator.DurationRuler'))...
                            );

                    catch
                        % unable to reuse existing axes rulers. Since they are read-only, we'll have
                        % to recreate.  However, first stash some information that will be lost, so
                        % that the figure appears to only change instead of being recreated.

                        % what is this scatter plot called again?
                        t=obj.myscatter.Tag;
                        ttl=get(obj.ax.Title);
                        % what is the state of the legend?
                        reshowLegend= ~isempty(obj.ax.Legend) && obj.ax.Legend.Visible =="on";
                        if reshowLegend
                            legendLocation = obj.ax.Legend.Location;
                        end

                        % recreate the axes and scatter plot
                        cla(obj.ax);
                        obj.myscatter=[];
                        obj.scatter(t);

                        obj.ax.Title.String=ttl.String;
                        obj.ax.Title.Color=ttl.Color;
                        obj.ax.Title.FontSize=ttl.FontSize;
                        obj.ax.Title.FontName=ttl.FontName;
                        obj.ax.Title.FontWeight=ttl.FontWeight;

                        % restore the legend
                        if reshowLegend
                            legend(obj.ax,'show','Location',legendLocation);
                        end
                    end

                %end

                function enforce_linear_scale_if_necessary()
                    xyzscale=[where(1) 'Scale'];
                    if obj.ax.(xyzscale)=="linear"
                        return
                    end

                    if iscell(c.(fld)) || iscategorical(c.(fld))
                        beep;
                        disp(['enforcing linear ' xyzscale ' because of data type']);
                        obj.ax.(xyzscale)='linear';
                    elseif isnumeric(c.(fld))&&any(c.(fld)<=0) && obj.ax.(xyzscale) == "log"
                        beep;
                        disp(['enforcing linear ' xyzscale ' because of negative values']);
                        obj.ax.(xyzscale)='linear';
                    end
                end

                function [n, scalename] = duration2numbers(d)
                    % put the duration on a reasonable scale
                    persistent logic
                    if isempty(logicparts)
                        % do change how durations are displayed, change this table
                        logictable={... scalename , converterFcn , minValue
                            'years', @years, years(4);
                            'months',@months, months(6);
                            'days',  @days, days(5);
                            'hours', @hours, hours(3);
                            'minutes', @minutes, minutes(2);
                            'seconds', @seconds, seconds(-inf)};

                        logic = cell2struct(logictable,{'name','fn','minval'},2);
                    end
                    idx = 1;
                    mymax = max(d);

                    % min_dur gives us the minimum duration value based on the index

                    while mymax < logic(idx).minval
                        idx = idx + 1;
                    end

                    scalename = logic(idx).name;
                    n = logic(idx).fn(d);
                end

            end % DOIT
        end

    end
    methods(Hidden)

        function set_valid_axes_choices(obj)
            c=obj.catalogFcn();
            p = properties(c);
            usable=true(size(p));
            nEvents=  c.Count;
            for i=1:numel(p)
                try
                    usable(i)=length(c.(p{i}))==nEvents;
                catch
                    usable(i)=false;
                end
                try
                    if all(isnan(c.(p{i})))
                        usable(i)=false;
                    end
                end
            end
            obj.axes_choices=p(usable);
        end

        function set_conversions(obj)
            for j=1:numel(obj.axes_choices)
                ch=obj.axes_choices{j};
            end

        end

        function updateCheckedStatus(obj,h, checkmask)
            labels = {h.Children.Label};
            choices=string(obj.axes_choices);
            for j=1:numel(h.Children)
                whichChoice = labels(j) == choices;
                set(h.Children(whichChoice), 'Checked', tf2onoff(checkmask(whichChoice)) );
            end
        end

        function xContextMenu(obj,xl, tag, fig)
            %delete(findobj(ancestor(ax),['xsel_ctxt' tag]));
            h=xl.UIContextMenu;
            checkmask = strcmp(obj.axes_choices, obj.x_by);
            if isempty(h)
                mytag = ['xsel_ctxt ' tag];
                delete(findobj(fig,'Type','uicontextmenu','-and','Tag',mytag));
                h=uicontextmenu('Tag',mytag);
                for i=1:numel(obj.axes_choices)
                    uimenu(h,'Label',obj.axes_choices{i}, 'MenuSelectedFcn', @(s,~)obj.change(s,'x_by'));
                end
                obj.add_axes_toggles(h,'X');
                xl.UIContextMenu=h;
            end
            obj.updateCheckedStatus(h, checkmask);
        end

        function yContextMenu(obj,yl,tag, fig)
            h=yl.UIContextMenu;
            checkmask = strcmp(obj.axes_choices, obj.y_by);
            if isempty(h)
                mytag = ['ysel_ctxt ' tag];
                delete(findobj(fig,'Type','uicontextmenu','-and','Tag',mytag));
                h=uicontextmenu('Tag',mytag);
                for i=1:numel(obj.axes_choices)
                    uimenu(h,'Label',obj.axes_choices{i},'MenuSelectedFcn',@(s,~)obj.change(s,'y_by'));
                end
                obj.add_axes_toggles(h,'Y');
                yl.UIContextMenu=h;
            end
            obj.updateCheckedStatus(h, checkmask);
        end

        function zContextMenu(obj,zl,tag, fig)
            h=zl.UIContextMenu;
            checkmask = strcmp(obj.axes_choices, obj.z_by);
            if isempty(h)
                mytag = ['zsel_ctxt ' tag];
                delete(findobj(fig,'Type','uicontextmenu','-and','Tag',mytag));
                h=uicontextmenu('Tag',mytag);
                for i=1:numel(obj.axes_choices)
                    uimenu(h,'Label',obj.axes_choices{i}, 'MenuSelectedFcn', @(s,v)obj.change(s,v,'z_by'));
                end
                obj.add_axes_toggles(h,'Z');
                zl.UIContextMenu=h;
            end
            obj.updateCheckedStatus(h, checkmask);
        end

        function scatterContextMenu(obj,sc,tag)
            tag=['ssel_ctxt ' tag];
            f = ancestor(obj.ax,'figure');
            delete(findobj(f,'Tag',tag));
            h=uicontextmenu(f,'Tag',tag);
            szm = uimenu(h,'Label','Size by...',...
                'MenuSelectedFcn', @(src,~)obj.cleanChildren_cb(src, 'size_by'));
            clm = uimenu(h,'Label','Color by...',...
                'MenuSelectedFcn', @(src,~)obj.cleanChildren_cb(src, 'color_by'));
            obj.sizeContextMenu(szm);
            obj.colorContextMenu(clm);
            sc.UIContextMenu=h;
        end
        function cleanChildren_cb(obj,src,bywhat)
            m=findobj(src.Children,'Type','uimenu');
            labels = get(m,'Label');
            ison = get(m,'Checked') == "on";
            isoff = ~ison;
            checkmask = strcmp(labels, obj.(bywhat));
            disp(labels{checkmask})
            set(m(~checkmask & ison),'Checked','off');
            set(m(checkmask & isoff),'Checked','on');
        end

        function sizeContextMenu(obj,h)
            for i=1:numel(obj.axes_choices)
                uimenu(h,'Label',obj.axes_choices{i},'MenuSelectedFcn',@obj.changeSize);
            end
            uimenu(h,'Separator','on','Label','Single Size','MenuSelectedFcn',@obj.changeSize);

        end

        function colorContextMenu(obj,h)
            for i=1:numel(obj.axes_choices)
                uimenu(h,'Label',obj.axes_choices{i},'MenuSelectedFcn',@obj.changeColor);
            end
            uimenu(h,'Separator','on','Label','-none-','MenuSelectedFcn',@obj.changeColor);
        end

        function add_axes_toggles(obj,h,letter)
            uimenu(h,'Label','Flip axes direction','Separator','on',...
                'MenuSelectedFcn',@(src,~)cb_axisdir(letter));
            uimenu(h,'Label','Toggle Log/Linear scale','Separator','on',...
                'MenuSelectedFcn',@(src,~)cb_axisscale(letter));

            function cb_axisdir(letter)
                dirs={'normal','reverse'};
                prop=[letter 'Dir'];
                dirs(strcmp(obj.ax.(prop),dirs))=[];
                obj.ax.(prop) = dirs{1};
            end

            function cb_axisscale(letter)
                scales={'linear','log'};
                prop=[letter 'Scale'];

                scales(strcmp(obj.ax.(prop),scales))=[];

                % if any data is non-positive, axes must remain linear
                if any(obj.myscatter.([letter 'Data']) <=0) && scales{1} == "log"
                    beep
                    disp(['enforcing linear ' prop ' because of negative values'])
                    obj.ax.(prop) = 'linear';
                else
                    obj.ax.(prop) = scales{1};
                end
            end
        end
        function change(obj,src,whatby)
            % whatby is x_by, y_by, etc...

            % remove checkmarks
            set(src.Parent.Children,'Checked','off');

            % change the plotting value
            obj.(whatby) = src.Label;

            % add new checkmark
            src.Checked='on';

            % relabel
            obj.ax.([upper(whatby(1)), 'Label']).String=src.Label;
            %h.String=src.Label;

            %replot
            obj.update(whatby);
        end

        function changeSize(obj,src,~)
            % whatby is x_by, y_by, etc...

            % change the plotting value
            obj.size_by= src.Label;

            % relabel
            set(gco,'DisplayName',sprintf('size:%s\ncolor:%s', obj.size_by, obj.color_by));
            switch(obj.size_by)
                case 'Date'
                    obj.sizeFcn=@(x) normalize(x,2,8,@(x)x.^2 );
                case 'Magnitude'
                    obj.sizeFcn=@mag2dotsize;
                case 'Single Size'
                    sz=num2str(round(mode(obj.myscatter.SizeData)));
                    val = str2double(inputdlg('Choose Marker Size','',1,{sz}));
                    if ~isempty(val) && ~isnan(val)
                        obj.sizeFcn=@(x)val;
                    end
                otherwise
                    obj.sizeFcn=@(x) normalize(x,2,8,@(x)x.^2);
            end
            %replot
            obj.update('size_by');

            function x=normalize(x, minval, scalar, modifier)
                if isa(x,'datetime')
                    x=x - min(x); % now duration
                end
                if isa(x,'duration')
                    x=days(x);
                end
                x = (x-min(x)) ./ range(x);
                x = modifier(x .* scalar + minval);
            end

        end
        function changeColor(obj,src,~)
            % whatby is x_by, y_by, etc...

            % remove checkmarks
            set(src.Parent.Children,'Checked','off');

            % change the plotting value
            obj.color_by = src.Label;

            % add new checkmark
            src.Checked='on';

            % relabel
            set(gco,'DisplayName',sprintf('size:%s\ncolor:%s',obj.size_by,obj.color_by));
            %h.String=src.Label;

            switch(obj.color_by)
                case 'Date'
                    obj.colorFcn=@datenum;
                case '-none-'
                    val=uisetcolor();
                    obj.colorFcn=@(x)val;
                otherwise
                    obj.colorFcn=@(x) normalize(x, 0, 1,@(x)x);
            end

            function x=normalize(x, minval, scalar, modifier)
                if iscategorical(x)
                    return;
                end
                if isdatetime(x)
                    x=x - min(x); % now duration
                end
                if isduration(x)
                    x=days(x);
                end
                x = (x-min(x)) / range(x);
                x = modifier(x.* scalar + minval);
            end
            %replot
            obj.update('color_by');
        end
    end % HIDDEN methods

    methods (Static)
        function s=instructions(varargin)
            % Explore your data in 4-5 dimensions by choosing a parameter for each axis.
            % Right-click on the X, Y, or Z axes for a list of available variables.
            % Right-click on data points to choose how they will be sized or colored.
            helpwin('CatalogExplorationPlot.instructions');
        end
    end
end

author: CelsoReyes


text:
function A = import_data(fname)
%IMPORT_DATA Import tensor-related data to a file.
%
%   A = IMPORT_DATA(FNAME) imports an object A from the file named FNAME.
%   The supported data types and formatting of the file are explained in
%   EXPORT_DATA.
%
%   See also TENSOR, EXPORT_DATA
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

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

%% Open file
fid = fopen(fname,'r');
if (fid == -1)
    error('Cannot open file %s',fname);
end

%% Get the type of object
type = import_type(fid);

%% Import the object

if strcmpi(type,'tensor')

    sz = import_size(fid);
    data = import_array(fid, prod(sz));
    A = tensor(data, sz);

elseif strcmpi(type,'sptensor')

    sz = import_size(fid);
    nz = import_nnz(fid);
    [subs, vals] = import_sparse_array(fid, length(sz), nz);
    A = sptensor(subs, vals, sz);

elseif strcmpi(type,'matrix') || strcmpi(type,'matrix')

    sz = import_size(fid);
    data = import_array(fid, prod(sz));
    A = reshape(data, sz);

elseif strcmpi(type,'ktensor')

    sz = import_size(fid);
    r = import_rank(fid);
    lambda = import_array(fid, r);
    U = {};
    for n = 1:length(sz)
        line = fgets(fid);
        fac_type = import_type(fid);
        fac_sz = import_size(fid);
        fac_data = import_array(fid, prod(fac_sz));
        % row wise reshape
        fac = reshape(fac_data, fliplr(fac_sz))';
        U{n} = fac;
    end
    A = ktensor(lambda,U);

else

    error('Invalid data type for export');

end

%% Close file
fclose(fid);

function type = import_type(fid)
% Import IO data type
line = fgets(fid);
typelist = regexp(line, '\s+', 'split');
type = typelist(1);

function sz = import_size(fid)
% Import the size of something from a file
line = fgets(fid);
n = sscanf(line, '%d');
line = fgets(fid);
sz = sscanf(line, '%d');
sz = sz';
if (size(sz,2) ~= n)
    error('Imported dimensions are not of expected size');
end

function nz = import_nnz(fid)
% Import the size of something from a file
line = fgets(fid);
nz = sscanf(line, '%d');

function r = import_rank(fid)
% Import the rank of something from a file
line = fgets(fid);
r = sscanf(line, '%d');

function data = import_array(fid, n)
% Import dense data that supports numel and linear indexing
data = fscanf(fid, '%e', n);

function [subs, vals] = import_sparse_array(fid, n, nz)
% Import sparse data subs and vals from coordinate format data
data = textscan(fid,[repmat('%d',1,n) '%n']);
subs = cell2mat(data(1:n));
vals = data{n+1};
if (size(subs,1) ~= nz)
    error('Imported nonzeros are not of expected size');
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% MP_INTERFACE_HCURL: create a global numbering of basis functions in multipatch domains, assigning the correct orientation for tangential traces.
%
%   [glob_num, glob_ndof, dofs_ornt] = mp_interface_hcurl (interfaces, sp);
%
% INPUT:
%
%  interfaces: structure with the information of the interfaces between patches (see mp_geo_read_file)
%  sp:         object representing the space of discrete functions, with the curl-preserving transform (see sp_vector)
%
% OUTPUT:
%
%  glob_num:  global numbering of the discrete basis functions
%  glob_ndof: total number of degrees of freedom
%  dofs_ornt: a cell-array with the orientation of the basis functions
%
% Copyright (C) 2010, 2011, 2015 Rafael Vazquez
%
%    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/>.

function [glob_num, glob_ndof, dofs_ornt] = mp_interface_hcurl (interfaces, sp)

  npatch = numel (sp);

  ndim = size (sp{1}.ndof_dir, 2);

  if (~isempty (interfaces))
    glob_num   = cell (npatch, 1);
    dofs_ornt  = cell (npatch, 1);
% patch_intrfc is only used for ndim == 3
    patch_intrfc = cell (npatch, 1);

    ttform = cell (ndim-1, npatch, numel (interfaces));
    ppnum  = cell (ndim-1, numel (interfaces));

    for iptc = 1:npatch
      glob_num{iptc} = zeros (1, sp{iptc}.ndof);
      dofs_ornt{iptc} = ones (1, sp{iptc}.ndof);
      patch_intrfc{iptc} = union (find([interfaces.patch1] == iptc), ...
                                  find([interfaces.patch2] == iptc));
    end

    for intrfc = 1:numel(interfaces)
      iptc1  = interfaces(intrfc).patch1;
      iptc2  = interfaces(intrfc).patch2;
      iside1 = interfaces(intrfc).side1;
      iside2 = interfaces(intrfc).side2;
      for idim = 1:ndim-1
        ttform{idim,iptc1, intrfc} = sp{iptc1}.boundary(iside1).dofs(sp{iptc1}.boundary(iside1).comp_dofs{idim});
        nghbr_dofs{idim} = reshape (sp{iptc2}.boundary(iside2).dofs(sp{iptc2}.boundary(iside2).comp_dofs{idim}), ...
				    [sp{iptc2}.boundary(iside2).ndof_dir(idim,:), 1]);
      end
      if (ndim == 2)
        if (interfaces(intrfc).ornt == -1)
          nghbr_dofs{1} = flipud (nghbr_dofs{1});
        end
      elseif (ndim == 3)
        if (interfaces(intrfc).flag == -1)
          aux = nghbr_dofs{1};
          nghbr_dofs{1} = nghbr_dofs{2}';
          nghbr_dofs{2} = aux';
        end
        if (interfaces(intrfc).ornt1 == -1)
          nghbr_dofs{1} = flipud (nghbr_dofs{1});
          nghbr_dofs{2} = flipud (nghbr_dofs{2});
        end
        if (interfaces(intrfc).ornt2 == -1)
          nghbr_dofs{1} = fliplr (nghbr_dofs{1});
          nghbr_dofs{2} = fliplr (nghbr_dofs{2});
        end
      end
      for idim = 1:ndim-1
        ttform{idim,iptc2, intrfc} = nghbr_dofs{idim}(:)';
        ppnum{idim,intrfc} = zeros (1, numel(sp{iptc1}.boundary(iside1).comp_dofs{idim}));
      end
    end

    glob_ndof = 0;
% We start with the dofs that do not belong to any interface
    for iptc = 1:npatch
      non_intrfc_dofs = setdiff(1:sp{iptc}.ndof, [ttform{:,iptc,:}]);
      glob_num{iptc}(non_intrfc_dofs) = glob_ndof + (1:numel(non_intrfc_dofs));
      glob_ndof = glob_ndof + numel (non_intrfc_dofs);
    end

% Then we set the interfaces
    for intrfc = 1:numel(interfaces)
      iptc     = interfaces(intrfc).patch1;
      iptc2    = interfaces(intrfc).patch2;

      for idim = 1:ndim-1
        new_dofs{idim} = find (glob_num{iptc}(ttform{idim, iptc, intrfc}) == 0);
        new_dofs_number = glob_ndof + (1:numel (new_dofs{idim}));
        glob_num{iptc}(ttform{idim, iptc, intrfc}(new_dofs{idim})) = new_dofs_number;
        dofs_ornt{iptc}(ttform{idim, iptc, intrfc}(new_dofs{idim})) = 1;
        ppnum{idim,intrfc}(new_dofs{idim}) = new_dofs_number;
        glob_ndof = glob_ndof + numel (new_dofs{idim});
      end

      if (ndim == 2)
        glob_num{iptc2}(ttform{1, iptc2, intrfc}) = ...
             glob_num{iptc}(ttform{1, iptc, intrfc});
        dofs_ornt{iptc2}(ttform{1, iptc2, intrfc}) = ...
             interfaces(intrfc).ornt * dofs_ornt{iptc}(ttform{1, iptc, intrfc});
      elseif (ndim == 3)
        [glob_num, dofs_ornt, ppnum] = set_same_interface (iptc, intrfc, ...
            ttform, new_dofs, interfaces, glob_num, dofs_ornt, ppnum, patch_intrfc);
        [glob_num, dofs_ornt, ppnum] = set_same_patch (iptc, intrfc, ...
            ttform, new_dofs, interfaces, glob_num, dofs_ornt, ppnum, patch_intrfc);

      end
    end

  else
    glob_ndof = 0;
    glob_num = cell (npatch, 1);
    dofs_ornt  = cell (npatch, 1);
    for iptc = 1:npatch
      glob_num{iptc} = glob_ndof + (1:sp{iptc}.ndof);
      glob_ndof = glob_ndof + sp{iptc}.ndof;
      dofs_ornt{iptc} = ones (1, sp{iptc}.ndof);
    end
  end

end

function [glob_num, dofs_ornt, ppnum] = set_same_patch (iptc, intrfc, ...
            ttform, new_dofs, interfaces, glob_num, ...
            dofs_ornt, ppnum, patch_intrfc);

  for idim = 1:size(ttform, 1)
    intrfc_dofs{idim} = ttform{idim, iptc, intrfc}(new_dofs{idim});
  end

  for ii = setdiff (patch_intrfc{iptc}, intrfc)
    for idim = 1:size(ttform, 1)
      jdim = mod (idim, 2) + 1;
      [~, pos1{idim}, pos2{idim}] = intersect (ttform{idim, iptc, ii}, intrfc_dofs{idim});
      not_set_p1{idim} = find (ppnum{idim,ii}(pos1{idim}) == 0);

      [~, pos1b{idim}, pos2b{jdim}] = ...
          intersect (ttform{idim, iptc, ii}, intrfc_dofs{jdim});
      not_set_p1b{idim} = find (ppnum{idim, ii}(pos1b{idim}) == 0);
    end

    if (~(all (cellfun (@isempty, not_set_p1))))
      for idim = 1:size(ttform, 1)
        ppnum{idim,ii}(pos1{idim}(not_set_p1{idim})) = ...
          ppnum{idim,intrfc}(new_dofs{idim}(pos2{idim}(not_set_p1{idim})));
        renew_dofs{idim} = pos1{idim}(not_set_p1{idim});
      end
      [glob_num, dofs_ornt, ppnum] = set_same_interface (iptc, ii, ...
          ttform, renew_dofs, interfaces, glob_num, dofs_ornt, ppnum, patch_intrfc);
    elseif (~(all (cellfun (@isempty, not_set_p1b))))
      for idim = 1:size(ttform, 1)
        jdim = mod (idim, 2) + 1;
        ppnum{idim,ii}(pos1b{idim}(not_set_p1b{idim})) = ...
          ppnum{jdim,intrfc}(new_dofs{jdim}(pos2b{jdim}(not_set_p1b{idim})));
        renew_dofs{idim} = pos1b{idim}(not_set_p1b{idim});
      end
      [glob_num, dofs_ornt, ppnum] = set_same_interface (iptc, ii, ...
          ttform, renew_dofs, interfaces, glob_num, dofs_ornt, ppnum, patch_intrfc);
    end
  end
end

function [glob_num, dofs_ornt, ppnum] = set_same_interface (iptc, intrfc, ...
            ttform, new_dofs, interfaces, glob_num, ...
            dofs_ornt, ppnum, patch_intrfc)

  ornt(1) = interfaces(intrfc).ornt1;
  ornt(2) = interfaces(intrfc).ornt2;

  iptc2 = setdiff ([interfaces(intrfc).patch1 interfaces(intrfc).patch2], iptc);
  for idim = 1:size(ttform, 1)
    intrfc_dofs{idim} = ttform{idim, iptc, intrfc}(new_dofs{idim});
    glob_num{iptc2}(ttform{idim, iptc2, intrfc}(new_dofs{idim})) = ...
                              glob_num{iptc}(intrfc_dofs{idim});
    dofs_ornt{iptc2}(ttform{idim, iptc2, intrfc}(new_dofs{idim})) = ...
             ornt(idim) * dofs_ornt{iptc}(intrfc_dofs{idim});
%             interfaces(intrfc).ornt(idim) * dofs_ornt{iptc}(intrfc_dofs{idim});
  end

  [glob_num, dofs_ornt, ppnum] = set_same_patch (iptc2, intrfc, ...
           ttform, new_dofs, interfaces, glob_num, ...
           dofs_ornt, ppnum, patch_intrfc);

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function f = simplify(f, varargin)
%SIMPLIFY  Simplifys a DELTAFUN object F.
%   F = SIMPLIFY(F) simplifies the FUNPART and removes any trivial delta
%   functions from the delta fun F. Preferences can be passed via VARARGIN.
%
% See also SIMPLIFY, CLEANROWS, CLEANCOLUMNS

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

f.funPart = simplify(f.funPart, varargin{:});
f = simplifyDeltas(f, varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function i4vec_decrement_test ( )

%*****************************************************************************80
%
%% I4VEC_DECREMENT_TEST tests I4VEC_DECREMENT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 January 2015
%
%  Author:
%
%    John Burkardt
%
  n = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_DECREMENT_TEST\n' );
  fprintf ( 1, '  I4VEC_DECREMENT decrements an I4VEC.\n' );

  v_lo = -5;
  v_hi = 10;
  seed = 123456789;
  [ v, seed ] = i4vec_uniform_ab ( n, v_lo, v_hi, seed );
  i4vec_print ( n, v, '  The I4VEC:' );
  v = i4vec_decrement ( n, v );
  i4vec_print ( n, v, '  The I4VEC after decrementing:' );

  return
end


text:
function [reg_c,rho_c,eta_c] = l_corner(rho,eta,reg_param,U,s,b,method,M)
%L_CORNER Locate the "corner" of the L-curve.
%
% [reg_c,rho_c,eta_c] =
%        l_corner(rho,eta,reg_param)
%        l_corner(rho,eta,reg_param,U,s,b,method,M)
%        l_corner(rho,eta,reg_param,U,sm,b,method,M) ,  sm = [sigma,mu]
%
% Locates the "corner" of the L-curve in log-log scale.
%
% It is assumed that corresponding values of || A x - b ||, || L x ||,
% and the regularization parameter are stored in the arrays rho, eta,
% and reg_param, respectively (such as the output from routine l_curve).
%
% If nargin = 3, then no particular method is assumed, and if
% nargin = 2 then it is issumed that reg_param = 1:length(rho).
%
% If nargin >= 6, then the following methods are allowed:
%    method = 'Tikh'  : Tikhonov regularization
%    method = 'tsvd'  : truncated SVD or GSVD
%    method = 'dsvd'  : damped SVD or GSVD
%    method = 'mtsvd' : modified TSVD,
% and if no method is specified, 'Tikh' is default.  If the Spline Toolbox
% is not available, then only 'Tikh' and 'dsvd' can be used.
%
% An eighth argument M specifies an upper bound for eta, below which
% the corner should be found.

% Per Christian Hansen, IMM, July 26, 2007.

% Set default regularization method.
if (nargin <= 3)
  method = 'none';
  if (nargin==2), reg_param = (1:length(rho))'; end
else
  if (nargin==6), method = 'Tikh'; end
end

% Set this logical variable to 1 (true) if the corner algorithm
% should always be used, even if the Spline Toolbox is available.
alwayscorner = 0;

% Set threshold for skipping very small singular values in the
% analysis of a discrete L-curve.
s_thr = eps;  % Neglect singular values less than s_thr.

% Set default parameters for treatment of discrete L-curve.
deg   = 2;  % Degree of local smooting polynomial.
q     = 2;  % Half-width of local smoothing interval.
order = 4;  % Order of fitting 2-D spline curve.

% Initialization.
if (length(rho) < order)
  error('Too few data points for L-curve analysis')
end
if (nargin > 3)
  [p,ps] = size(s); [m,n] = size(U);
  beta = U'*b;
  if (m>n), b0 = b - U*beta; end
  if (ps==2)
    s = s(p:-1:1,1)./s(p:-1:1,2);
    beta = beta(p:-1:1);
  end
  xi = beta./s;
end

% Restrict the analysis of the L-curve according to M (if specified).
if (nargin==8)
  index = find(eta < M);
  rho = rho(index); eta = eta(index); reg_param = reg_param(index);
end

if (strncmp(method,'Tikh',4) | strncmp(method,'tikh',4))

  % The L-curve is differentiable; computation of curvature in
  % log-log scale is easy.

  % Compute g = - curvature of L-curve.
  g = lcfun(reg_param,s,beta,xi);

  % Locate the corner.  If the curvature is negative everywhere,
  % then define the leftmost point of the L-curve as the corner.
  [gmin,gi] = min(g);
  reg_c = fminbnd('lcfun',...
    reg_param(min(gi+1,length(g))),reg_param(max(gi-1,1)),...
    optimset('Display','off'),s,beta,xi); % Minimizer.
  kappa_max = - lcfun(reg_c,s,beta,xi); % Maximum curvature.

  if (kappa_max < 0)
    lr = length(rho);
    reg_c = reg_param(lr); rho_c = rho(lr); eta_c = eta(lr);
  else
    f = (s.^2)./(s.^2 + reg_c^2);
    eta_c = norm(f.*xi);
    rho_c = norm((1-f).*beta);
    if (m>n), rho_c = sqrt(rho_c^2 + norm(b0)^2); end
  end

elseif (strncmp(method,'tsvd',4) | strncmp(method,'tgsv',4) | ...
        strncmp(method,'mtsv',4) | strncmp(method,'none',4))

  % Use the adaptive pruning algorithm to find the corner, if the
  % Spline Toolbox is not available.
  if ~exist('splines','dir') | alwayscorner
    %error('The Spline Toolbox in not available so l_corner cannot be used')
    reg_c = corner(rho,eta);
    rho_c = rho(reg_c);
    eta_c = eta(reg_c);
    return
  end

  % Othersise use local smoothing followed by fitting a 2-D spline curve
  % to the smoothed discrete L-curve. Restrict the analysis of the L-curve
  % according to s_thr.
  if (nargin > 3)
    if (nargin==8)       % In case the bound M is in action.
      s = s(index,:);
    end
    index = find(s > s_thr);
    rho = rho(index); eta = eta(index); reg_param = reg_param(index);
  end

  % Convert to logarithms.
  lr = length(rho);
  lrho = log(rho); leta = log(eta); slrho = lrho; sleta = leta;

  % For all interior points k = q+1:length(rho)-q-1 on the discrete
  % L-curve, perform local smoothing with a polynomial of degree deg
  % to the points k-q:k+q.
  v = (-q:q)'; A = zeros(2*q+1,deg+1); A(:,1) = ones(length(v),1);
  for j = 2:deg+1, A(:,j) = A(:,j-1).*v; end
  for k = q+1:lr-q-1
    cr = A\lrho(k+v); slrho(k) = cr(1);
    ce = A\leta(k+v); sleta(k) = ce(1);
  end

  % Fit a 2-D spline curve to the smoothed discrete L-curve.
  sp = spmak((1:lr+order),[slrho';sleta']);
  pp = ppbrk(sp2pp(sp),[4,lr+1]);

  % Extract abscissa and ordinate splines and differentiate them.
  % Compute as many function values as default in spleval.
  P     = spleval(pp);  dpp   = fnder(pp);
  D     = spleval(dpp); ddpp  = fnder(pp,2);
  DD    = spleval(ddpp);
  ppx   = P(1,:);       ppy   = P(2,:);
  dppx  = D(1,:);       dppy  = D(2,:);
  ddppx = DD(1,:);      ddppy = DD(2,:);

  % Compute the corner of the discretized .spline curve via max. curvature.
  % No need to refine this corner, since the final regularization
  % parameter is discrete anyway.
  % Define curvature = 0 where both dppx and dppy are zero.
  k1    = dppx.*ddppy - ddppx.*dppy;
  k2    = (dppx.^2 + dppy.^2).^(1.5);
  I_nz  = find(k2 ~= 0);
  kappa = zeros(1,length(dppx));
  kappa(I_nz) = -k1(I_nz)./k2(I_nz);
  [kmax,ikmax] = max(kappa);
  x_corner = ppx(ikmax); y_corner = ppy(ikmax);

  % Locate the point on the discrete L-curve which is closest to the
  % corner of the spline curve.  Prefer a point below and to the
  % left of the corner.  If the curvature is negative everywhere,
  % then define the leftmost point of the L-curve as the corner.
  if (kmax < 0)
    reg_c = reg_param(lr); rho_c = rho(lr); eta_c = eta(lr);
  else
    index = find(lrho < x_corner & leta < y_corner);
    if ~isempty(index)
      [dummy,rpi] = min((lrho(index)-x_corner).^2 + (leta(index)-y_corner).^2);
      rpi = index(rpi);
    else
      [dummy,rpi] = min((lrho-x_corner).^2 + (leta-y_corner).^2);
    end
    reg_c = reg_param(rpi); rho_c = rho(rpi); eta_c = eta(rpi);
  end

elseif (strncmp(method,'dsvd',4) | strncmp(method,'dgsv',4))

  % The L-curve is differentiable; computation of curvature in
  % log-log scale is easy.

  % Compute g = - curvature of L-curve.
  g = lcfun(reg_param,s,beta,xi,1);

  % Locate the corner.  If the curvature is negative everywhere,
  % then define the leftmost point of the L-curve as the corner.
  [gmin,gi] = min(g);
  reg_c = fminbnd('lcfun',...
    reg_param(min(gi+1,length(g))),reg_param(max(gi-1,1)),...
    optimset('Display','off'),s,beta,xi,1); % Minimizer.
  kappa_max = - lcfun(reg_c,s,beta,xi,1); % Maximum curvature.

  if (kappa_max < 0)
    lr = length(rho);
    reg_c = reg_param(lr); rho_c = rho(lr); eta_c = eta(lr);
  else
    f = s./(s + reg_c);
    eta_c = norm(f.*xi);
    rho_c = norm((1-f).*beta);
    if (m>n), rho_c = sqrt(rho_c^2 + norm(b0)^2); end
  end

else
  error('Illegal method')
end

repo: REST


text:
% DEMUSPSVARGPLVM1 Demonstrate variational GPLVM on USPS data.
% Copyright: Michalis Titsias, Neil Lawrence, Andreas Damianou, 2010 - 2014
% VARGPLVM

% Fix seeds
rng(1e5, 'v4');

dataSetName = 'usps';
experimentNo = 1;
printDiagram = 1;

% load data
[YTrain, lblsTrain, YTest, lblsTest] = lvmLoadData(dataSetName);

nClasses = size(lblsTrain,2);
% Set up model
options = vargplvmOptions('dtcvar');
options.kern = {'rbfard2', 'white'};
options.numActive = 50;
%options.tieParam = 'tied';

options.optimiser = 'scg';
latentDim = 10;
d = size(YTrain, 2);

capName = dataSetName;
capName(1) = upper(capName(1));
modelType = 'vargplvm'; % varmodel{1}.type;
modelType(1) = upper(modelType(1));
modelFile = ['dem' capName modelType num2str(experimentNo) '.mat'];

if exist(modelFile, 'file') % Training data exists, just load it.
    load(modelFile);
else % Do the training:
    iters = 1000;
    display = 1;

    varmodel = cell(nClasses,1);
    % create a separate vargplvm for each digit
    for i=1:nClasses
        %
        fprintf('\n\n# Training for Class # %d\n\n',i)
        Y = YTrain(lblsTrain(:,i)==1,:);

        model = vargplvmCreate(latentDim, d, Y, options);
        %
        model = vargplvmParamInit(model, model.m, model.X);

        model = vargplvmOptimise(model, display, iters);

        varmodel{i} = model;
        %
    end
    save(modelFile, 'varmodel');
end

iters = 100;
display = 0;

% measure performance on test data
prob = zeros(size(YTest,1),nClasses);
TestError = 0;

% New and faster way: give the whole test matrix at once, rather than using
% a for loop through each test element. This is almost 10 times faster then
% the old way (using 2 Matlab workers on a 4-core machine).

% Compute the approximate class conditional density for each digit
tic
for i=1:nClasses
    sprintf('--------\nComputing test probability for class=%d\n', i);
    prob(:,i) = vargplvmProbabilityCompute(varmodel{i}, YTest, 0, iters);
end
[maxP C] = max(prob,[],2);
for n=1:size(YTest,1)
    TestError = TestError + ~lblsTest(n,C(n));
end
fprintf('TestError=%d,\tTime=%fseconds.\n', TestError, toc);
%
save(['dem' capName modelType num2str(experimentNo) '.mat'], 'varmodel', 'prob', 'TestError');

%% Nearest Neighbour baseline
dists = dist2(YTrain,YTest);
[~,positions] = min(dists);
Ypred_NN = lblsTrain(positions,:);
NNError = 0;
for n=1:size(YTest,1)
    pp = find(Ypred_NN(n,:));
    NNError = NNError + ~lblsTest(n, pp);
end

%% Training of logistic regression classifier
labelsTrain = transformLabels(lblsTrain)';
labelsTest = transformLabels(lblsTest)';
Nstar = size(YTest,1);

for i=1:nClasses
    fprintf('\n # LogReg training for class # %d\n', i)
    lb = zeros(size(YTrain,1),1);
    lb(labelsTrain == i) = 1;
    B{i} = glmfit(YTrain, lb,'binomial','logit'); % Logistic regression
end

% Prediction of each binary classifier
Ypred_logReg = zeros(size(lblsTest));
for i=1:nClasses
    Ypred_logReg(:,i) = glmval(B{i},YTest,'logit')';
end

% Replace predictions with maximum probability (ie, make a decision)
[~,ind]=max(Ypred_logReg');
LogRegError = 0;
for i=1:size(YTest,1)
    LogRegError = LogRegError + (ind(i) ~= labelsTest(i));
end

% Print all results
fprintf('Var-GPLVM = %d (%f %%) \n', TestError, TestError*100/Nstar);
fprintf('NN        = %d mistakes (%f %% misclassified) \n', NNError, NNError*100/Nstar);
fprintf('Log. Reg. = %d mistakes (%f %% misclassified) \n', LogRegError, LogRegError*100/Nstar);

author: SheffieldML

repo: vargplvm


text:
function optimo = LocalSearch(Problem,pos,w)

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

    MaxIter = 5;
    Tol     = 1e-3;
    step    = 1;
    k       = 1;
    error   = 10;
    while error>Tol && k<MaxIter
        grad(1,:)    = FiniteDifference(pos,w,Problem);
        offspringdec = pos.dec - step*grad(1,:);
        offspringdec = min(max(offspringdec,Problem.lower),Problem.upper);
        offspring    = Problem.Evaluation(offspringdec);
        grad(2,:)    = FiniteDifference(offspring,w,Problem);
        step         = abs((offspring.dec-pos.dec)*(grad(2,:)-grad(1,:))')/norm(grad(2,:)-grad(1,:))^2;
        error = norm(offspring.dec-pos.dec);
        pos   = offspring;
        k     = k + 1;
    end
    optimo = pos;
end

function df = FiniteDifference(X,W,Problem)
    if any(X.con>0)
        df = Problem.CalConGrad(X.dec)';
        df = sum(df,2);
    else
        df = Problem.CalObjGrad(X.dec)';
        df = df*W';
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function tf = strlexcmp(a, b)
%STRLEXCMP Lexicographic comparison of two strings.
%
%   STRLEXCMP(A, B) returns -1, 0, or 1 depending on whether the left argument
%   is stringwise less than, equal to, or greater than the right argument.
%
%   This is a MATLAB version of the Perl `cmp' operator.
%
%   See also EQ, ISEQUAL.

%   Author:      Peter J. Acklam
%   Time-stamp:  2003-10-13 11:16:31 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   % check arguments
   error(nargchk(2, 2, nargin));
   if ~ischar(a) | ~ischar(b)
      error('Both arguments must be char arrays (strings).');
   end

   % get lengths of strings
   na = length(a);
   nb = length(b);
   n = min(na, nb);

   % find characters that differ
   k = find(a(1:n) ~= b(1:n));
   if isempty(k)
      % all characters are identical -- compare lengths
      tf = sign(na - nb);
   else
      % compare first character that is different
      k = k(1);
      tf = sign(a(k) - b(k));
   end


text:
%% Modelo CNC
%% cargar modelo
    model_cnc = vrworld('cnc5.wrl');
    open(model_cnc)
    cnc3D = view(model_cnc);

%% Definicion de Nodos
    H.puntero = vrnode(model_cnc, 'puntero');
    H.riel_horiz = vrnode(model_cnc, 'riel_horizontal');

% %% Inicializacion de objectos
%     H.riel_horiz.translation = [-5 0.3 0];
%     H.puntero.translation = [ 0 0 -0.5];

 %% rutina
    H.riel_horiz.translation = [-5 0.3 0];
    H.puntero.translation = [ 0 2.5 -0.5];

    %% desplazamiento del puntero atravez del riel que lo traslada

    pause(3)
    despl_p = 2.5:-0.01:-2.5;
    tan_desp = length(despl_p);
    global Puntero
        Puntero = H.puntero;
    % Avanzar de izquierda a derecha
%     fin = Despl_Puntero(1,100);
%     pause(2)
%     Bajar_Puntero()
%     pause(2)
%     fin1 = Despl_Puntero(100, 250);
%     pause(2)
%     Subir_Puntero()
%     fin2 = Despl_Puntero(250, 349);
%     Bajar_Puntero()
%     pause(2)
%     fin3 = Despl_Puntero(350, 400);
%     pause(2)
% %     Subir_Puntero()
% %
% %     msgbox('Retornando', 'Retornando')
% %     Retorno_Puntero();

%% Avanzar un paso
    global Riel
        Riel = H.riel_horiz;
%     Avanzar()
%     pause(4)
%     Avanzar()
%     Subir_Puntero()
%     msgbox('Retornando', 'Retornando')
%     Retorno_Puntero();

author: Sable

repo: mcbench-benchmarks


text:
function rtk=udstate_pppins(rtk,obs,nav)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% update the state parameters (PPP/INS mode)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global glc

% update positon
rtk=udpos_pppins(rtk);

% update clock(include GLONASS icb)
rtk=udclk_pppins(rtk);

% update tropospheric parameters
if rtk.opt.tropopt==glc.TROPOPT_EST||rtk.opt.tropopt==glc.TROPOPT_ESTG
    rtk=udtrop_pppins(rtk);
end

% update ionospheric parameters
if rtk.opt.ionoopt==glc.IONOOPT_EST
    rtk=udiono_pppins(rtk,obs,nav);
end

% update L5-receiver-dcb parameters
if rtk.opt.nf>=3
    rtk=uddcb_pppins(rtk);
end

% update ambiguity
rtk=udamb_pppins(rtk,obs,nav);

return

repo: GINav