text:
function [opts] = IPM_SC ( ipts , par)
% ========================================================================
% Spaese Coding of FDDL, Version 1.0
% Copyright(c) 2011  Meng YANG, Lei Zhang, Xiangchu Feng and David Zhang
% All Rights Reserved.
%
% -----------------------------------------------------------------------
% Permission to use, copy, or modify this software and its documentation
% for educational and research purposes only and without fee is here
% granted, provided that this copyright notice and the original authors'
% names appear on all copies and supporting documentation. This program
% shall not be used, rewritten, or adapted as the basis of a commercial
% software or hardware product without first obtaining permission of the
% authors. The authors make no representations about the suitability of
% this software for any purpose. It is provided "as is" without express
% or implied warranty.
%----------------------------------------------------------------------
%
% This is an implementation of the algorithm to do sparse coding
%
% Please refer to the following paper
%
% Meng Yang, Lei Zhang, Xiangchu Feng, and David Zhang,"Fisher Discrimination
% Dictionary Learning for Sparse Representation", In IEEE Int. Conf. on
% Computer Vision, 2011.
% L. Rosasco, A. Verri, M. Santoro, S. Mosci, and S. Villa. Iterative
% Projection Methods for Structured Sparsity Regularization. MIT Technical
% Reports, MIT-CSAIL-TR-2009-050,CBCL-282, 2009.
% J. Bioucas-Dias, M. Figueiredo, ?A new TwIST: two-step iterative shrinkage
% /thresholding  algorithms for image restoration?, IEEE Transactions on
% Image Processing, December 2007.
%----------------------------------------------------------------------
%
%  Inputs :   (1) ipts :    the structre of input data
%                    .D     the dictionary
%                    .tau1    parameter of l1-norm energy of coefficient
%                    .y     the testing sample
%             (2) par :     the struture of input parameters
%                    .eigenv
%
% Outputs:    (1) opts :    the structure of output data
%                    .A     the coefficient matrix
%                    .ert   the total energy sequence
%
%---------------------------------------------------------------------

ipts.cT         =      1e+10;  % stop criterion
ipts.citeT      =      1e-6;   % stop criterion
par.nIter       =      200;    % maximal iteration number
par.cRatio      =      1.05;
ipts.initM      =      'zero'; % coefficiet initialization method

switch lower(ipts.initM)
    case {'zero'}
        x(:,1)  =  zeros(size(ipts.D,2),1);
    case {'transpose'}
        x(:,1)  =  ipts.D'*ipts.y;
    case {'pinv'}
        x(:,1)  =  pinv(ipts.D)*ipts.y;
    otherwise
        error('Nonknown method!');
end

A         =         ipts.D;
tau1      =         ipts.tau1;
y         =         ipts.y;
nIter     =         par.nIter;
if ~isfield(par,{'eigenv'})
c         =         par.cRatio*find_max_eigenv(A'*A);
else
c         =         par.cRatio*par.eigenv ;
end

%%%%%%%%%%%%%%%%%%%%
% TWIST parameter
%%%%%%%%%%%%%%%%%%%%
for_ever           =         1;
IST_iters          =         0;
TwIST_iters        =         0;
sparse_sign        =         1;
verbose            =         1;
enforceMonotone    =         1;
lam1               =         1e-4;   %default eigenvalues
lamN               =         1;      %default eigenvalues
rho0               =         (1-lam1/lamN)/(1+lam1/lamN);
alpha              =         2/(1+sqrt(1-rho0^2));        %default,user can set
beta               =         alpha*2/(lam1+lamN);         %default,user can set

%%%%%%%%%%%%%%%%%%%
%main loop
%%%%%%%%%%%%%%%%%%%
xm2       =      x(:,1);
xm1       =      x(:,1);
prev_f    =      norm(y-A*x(:,1),2)^2+...
                    tau1*norm(x(:,1),1);
for n_it = 2 : nIter;

   ert(n_it-1)      =       norm(y-A*x(:,n_it-1),2)^2+...
        tau1*norm(x(:,n_it-1),1);
   fprintf('Iteration:%f  Total gap:%f\n',n_it,ert(n_it-1));

    while for_ever
       % IST(IPM) estimate
       v        =       1/c*A'*(y-A*xm1)+xm1;
       x_temp   =      soft(v,tau1/c/2);

       if (IST_iters >= 2) | ( TwIST_iters ~= 0)
            % set to zero the past when the present is zero
            % suitable for sparse inducing priors
            if sparse_sign
                mask    =   (x_temp ~= 0);
                xm1     =   xm1.* mask;
                xm2     =   xm2.* mask;
            end
            % two-step iteration
            xm2    =   (alpha-beta)*xm1 + (1-alpha)*xm2 + beta*x_temp;
            % compute residual
            f      =   norm(y-A*xm2,2)^2+...
                       tau1*norm(xm2,1);
            if (f > prev_f) & (enforceMonotone)
                TwIST_iters   =  0;  % do a IST(IPM) iteration if monotonocity fails
            else
                TwIST_iters =   TwIST_iters+1; % TwIST iterations
                IST_iters   =    0;
                x_temp      =   xm2;
                if mod(TwIST_iters,10000) ==0
                   c = 0.9*c;
                end
                break;  % break loop while
            end
        else
            f     =     norm(y-A*x_temp,2)^2+...
                        tau1*norm(x_temp,1);
            if f > prev_f
                % if monotonicity  fails here  is  because
                % max eig (A'A) > 1. Thus, we increase our guess
                % of max_svs
                c         =    2*c;
                if verbose
%                     fprintf('Incrementing c=%2.2e\n',c);
                end
                if  c > ipts.cT
                    break;  % break loop while
                end
                IST_iters = 0;
                TwIST_iters = 0;
            else
                TwIST_iters = TwIST_iters + 1;
                break;  % break loop while
            end
        end
    end

    citerion      =   abs(f-prev_f)/prev_f;
    if citerion < ipts.citeT | c > ipts.cT
%        fprintf('Stop!\n c=%2.2e\n citerion=%2.2e\n',c,citerion);
       break;
    end

    xm2           =   xm1;
    xm1           =   x_temp;
    x(:,n_it)   =   x_temp;
    prev_f        =   f;
end

opts.x     =       x(:,end);
opts.ert   =       ert;

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function visualizeData( X, k, IDX, types, C )
% Project high dim. data unto principal components (PCA) for visualization.
%
% Optionally IDX can be specified to indicate different classes for the
% points; in this case points in different classes are displayed using
% different colors. Up to 12 types are handled (for technical reasons
% involving plot), any cluster with a label>12 is assigned the label 12.
%
% USAGE
%  visualizeData( X, k, [IDX], [types], [C] )
%
% INPUTS
%  X       - column vector of data - N vectors of dimension p (X is Nxp)
%  k       - dimension to which to reduce data (2 or 3)
%  IDX     - [] cluster membership [see kmeans2.m]
%  types   - [] cell array of length ntypes of text labels for each type
%  C       - [] cluster centers (Kxp)
%
% OUTPUTS
%
% EXAMPLE
%  X = [randn(100,5); randn(100,5)+4];
%  C = [mean(X(1:100,:)); mean(X(101:200,:))];
%  IDX = [ones(100,1); 2*ones(100,1)];
%  visualizeData( X, 2, IDX, {'type1','type2' }, C);
%
% See also KMEANS2, DEMOCLUSTER
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<3 || isempty(IDX) ); IDX=[]; end
if( nargin<4 || isempty(types) ); types=[]; end
if( nargin<5 || isempty(C) ); C=[]; end

% apply PCA if necessary
if( size(X,2)~= k )
  [ U, mu ] = pca( X' );
  X = pcaApply( X', U, mu, k )';
  if(~isempty(C)); C = pcaApply( C', U, mu, k )'; end
end

%%% get k
k = size(X,2);
if( k==1 ); X = [X zeros(size(X))]; k = 2; end
if( k>3 ); error( 'k must be <= 3'); end

%%% show points
if( isempty(IDX) )
  if( k==2 )
    plot( X(:,1), X(:,2), '.' );
  elseif( k==3 )
    plot3( X(:,1), X(:,2), X(:,3), '.' );
  end;

else
  IDX(IDX>12)=12;  m=max(IDX);

  if( k==2)
    % plot points
    R = cell(1,3*m+3);
    for i=1:m;
      R((3*i-2):(3*i)) = {X(IDX==i,1), X(IDX==i,2), '.'};
    end;
    R((3*m+1):(3*m+3)) = {X(IDX==-1,1), X(IDX==-1,2), 'k.'};
    plot( R{:} );

    % plot centers
    if( ~isempty(C) )
      R=cell(1,3*m);
      for i=1:m;  R((3*i-2):(3*i)) = {C(i,1), C(i,2), 'x'}; end
      hold('on');  plot( R{:}, 'MarkerSize', 30 );  hold('off');
    end;

  elseif( k==3 )
    % plot points
    R = cell(1,4*m+4);
    for i=1:m;
      R((4*i-3):(4*i)) = {X(IDX==i,1), X(IDX==i,2), X(IDX==i,3), '.'};
    end;
    R((4*m+1):(4*m+4)) = {X(IDX==-1,1), X(IDX==-1,2), X(IDX==-1,3), 'k.'};
    plot3( R{:} );

    % plot centers
    if( ~isempty(C) )
      R=cell(1,4*m);
      for i=1:m;  R((4*i-3):(4*i)) = {C(i,1), C(i,2), C(i,3), 'x'}; end
      hold('on'); plot3( R{:}, 'MarkerSize', 30 );  hold('off');
    end
  end
end
axis('equal');

%%% show legend if types is provided
if(~isempty(types));  legend(types); end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function [params, names] = uniformPriorExtractParam(prior)

% UNIFORMPRIOREXTRACTPARAM Extract params from uniform prior structure.

% PRIOR

params = [];
if nargout > 1
  names = {};
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function Y=runMLP(X,Wx,Wy)
% The matrix implementation of the two-layer Multilayer Perceptron (MLP) neural networks.
%
% Author: Marcelo Augusto Costa Fernandes
% DCA - CT - UFRN
% mfernandes@dca.ufrn.br
%
% Input parameters:
%   X: Input neural network.  X is a (p x K) dimensional matrix, where p is a number of the inputs and K >= 1.
%   Wx: Hidden layer weight matrix. Wx is a (H x p+1) dimensional matrix.
%   Wy: Output layer weight matrix. Wy is a (m x H+1) dimensional matrix.
%
% Output parameters:
%  Y: Outpuy neural network.  Y is a (m x K) dimensional matrix, where m is a number of the output neurons and K >= 1.

[p1 N] = size (X);

bias = -1;

X = [bias*ones(1,N) ; X];

V = Wx*X;
Z = 1./(1+exp(-V));

S = [bias*ones(1,N);Z];
G = Wy*S;

Y = 1./(1+exp(-G));

author: Sable

repo: mcbench-benchmarks


text:
function [ n_data, x, fx ] = shi_values ( n_data )

%*****************************************************************************80
%
%% SHI_VALUES returns some values of the hyperbolic sine integral function.
%
%  Discussion:
%
%    SHI(X) = integral ( 0 <= T <= X ) sinh ( T ) / T dt
%
%    In Mathematica, the function can be evaluated by:
%
%      SinhIntegral[x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Cambridge University Press, 1999,
%    ISBN: 0-521-64314-7,
%    LC: QA76.95.W65.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 16;

  fx_vec = [ ...
    0.5069967498196672, ...
    0.6121303965633808, ...
    0.7193380189288998, ...
    0.8289965633789345, ...
    0.9414978265114335, ...
    1.057250875375729, ...
    1.300250361022057, ...
    1.561713388361002, ...
    1.845814141358504, ...
    2.157290343425901, ...
    2.501567433354976, ...
    3.549340406224435, ...
    4.973440475859807, ...
    6.966162067504942, ...
    9.817326911233034, ...
    13.96788504934715  ];

  x_vec = [ ...
      0.5E+00, ...
      0.6E+00, ...
      0.7E+00, ...
      0.8E+00, ...
      0.9E+00, ...
      1.0E+00, ...
      1.2E+00, ...
      1.4E+00, ...
      1.6E+00, ...
      1.8E+00, ...
      2.0E+00, ...
      2.5E+00, ...
      3.0E+00, ...
      3.5E+00, ...
      4.0E+00, ...
      4.5E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ode,alg] = dae(daefh, ...
  x_struct, z_struct, u_struct, p_struct, x_order, ...
  x, z, u, p, userdata)
% evaluate the system equations for the assigned variables

x = ocl.Variable.create(x_struct,x);
z = ocl.Variable.create(z_struct,z);
u = ocl.Variable.create(u_struct,u);
p = ocl.Variable.create(p_struct,p);

daehandler = ocl.DaeHandler(userdata);
daefh(daehandler,x,z,u,p);

nx = length(x_struct);
nz = length(z_struct);

ode = daehandler.getOde(nx, x_order);
alg = daehandler.getAlg(nz);

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function hole_num = p09_hole_num ( )

%*****************************************************************************80
%
%% P09_HOLE_NUM counts the holes in problem 09.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer HOLE_NUM, the number of holes.
%
  hole_num = 2;

  return
end


text:
function F = metoffice_experiment_gpfa_gpkron2(datanum, anomalies)

% DEBUGGING!!
debug = false;
if debug
  maxiter = 1;
  N_samples = 2;
else
  maxiter = 1000;
  N_samples = 1000;
end

%
% Load the data
%

[data,dataset,folder,maskfile] = metoffice_getdata(datanum, anomalies);
folder = [folder '/gpfa_gpkron2'];
mkdir(folder);
date = datestr(now,'yyyymmdd');

if anomalies
  disp('Model anomalies')
  comps = [0 10 0];
  comps_spatial = [0 10 0];
else
  disp('Don''t remove climatological averages')
  comps = [5 10 0];
  comps_spatial = [2 13 0];
end
% Number of components
D = sum(comps);

% Form the data matrix
Y = data.data;
[M,N] = size(Y);
Obs = ~isnan(Y);
weights = repmat(data.gridsize, [1, N]);

if false
  gpfafile = ['/share/work/jluttine/metoffice/rectest/gpfa_gpkron2/' ...
              'results_rectest_hadsst2d1_gpfa_D=6_anomalies=0_20110318'];
  D = 6;
else
  gpfafile = [];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GPFA inference
%

if isempty(gpfafile)

  %
  % Model for temporal X
  %

  covfunc_x = cell(D,1);
  theta_x = cell(D,1);
  is_pseudos_x = false(D,1);

  % Inputs (assume uniformly spaced time instances which is not exactly correct)
  in_x = linspace(data.time(1), data.time(end), length(data.time));
  pseudo_x = in_x(1:10:end);

  % Squared distances for covariance functions
  D2_pp = sq_dist(pseudo_x);
  D2_px = sq_dist(pseudo_x, in_x);
  d2_x = zeros(size(in_x,2),1); %diag(D2_xx);
  D_pp = sqrt(D2_pp);
  D_px = sqrt(D2_px);
  d_x = sqrt(d2_x);

  % Distance matrices for covariance functions
  d_xx = sqrt(sq_dist(in_x(:,1),in_x));

  ind = 0;

  % Periodic components (1 year period) with decay (rational quadratic)
% $$$ ind = 1:min(D,comps(1));
% $$$ fprintf('%d periodic components for X\n', length(ind));
% $$$ covfunc = @(D,D2) gp_cov_product(gp_cov_periodic(D, 'wavelength', 365), ...
% $$$                                  gp_cov_rq(D2));
% $$$ covfunc_x(ind) = {gp_cov_pseudo(gp_cov_jitter(covfunc(D_pp,D2_pp), 1e-6), ...
% $$$                                 covfunc(D_px,D2_px), ...
% $$$                                 covfunc(d_x,d2_x))};
% $$$ theta_x(ind) = columns_to_cells(...
% $$$     [linspace(1,1,length(ind))      % smoothness of the period
% $$$      365*linspace(20,1,length(ind)) % lengthscale of the decay (RQ)
% $$$      ones(1,length(ind))]);         % alpha (RQ)
% $$$ is_pseudos_x(ind) = true;
% Periodic components (1 year period) WITHOUT decay
  if comps(1) > 0
    ind = 1:min(D,comps(1));
    fprintf('%d periodic components for X\n', length(ind));
    covfunc = @(D) gp_cov_periodic(D, 'wavelength', 365);
    covfunc_x(ind) = {gp_cov_pseudo(gp_cov_jitter(covfunc(D_pp), 1e-6), ...
                                    covfunc(D_px), ...
                                    covfunc(d_x))};
    theta_x(ind) = columns_to_cells(...
        [linspace(1,0.1,length(ind))]);      % smoothness of the period
    is_pseudos_x(ind) = true;
  end

  % Slow components (1-20 years): rational quadratic using pseudo inputs
  if comps(2) > 0
    ind = ind(end) + (1:min(D-ind(end),comps(2)));
    fprintf('%d slow components for X\n', length(ind));
% $$$ covfunc = @(D2) gp_cov_se(D2); % SE
    covfunc = @(D2) gp_cov_rq(D2); % RQ
    covfunc_x(ind) = {gp_cov_pseudo(gp_cov_jitter(covfunc(D2_pp), 1e-3), ...
                                    covfunc(D2_px), ...
                                    covfunc(d2_x))};
% $$$ theta_x(ind) = columns_to_cells(...
% $$$     [365*linspace(20,1,length(ind))]); % lengthscale for SE
    theta_x(ind) = columns_to_cells(...
        [365*linspace(4,1,length(ind)) % lengthscale for RQ
         ones(1,length(ind))]);         % alpha for RQ
    is_pseudos_x(ind) = true;
  end

  % Fast components (4-18 months): piecewise polynomial in 1-D
  % Take advantage of the Toeplitz structure of the covariance matrix
  if comps(3) > 0
    ind = (ind(end)+1):D;
    fprintf('%d fast components for X\n', length(ind));
    covfunc_x(ind) = {gp_cov_jitter(gp_cov_toeplitz(gp_cov_pp(d_xx,1)))};
    theta_x(ind) = columns_to_cells(...
        [30*linspace(12,4,length(ind))]); % lengthscale or cut-off
    is_pseudos_x(ind) = false;
  end

  %
  % Model for spatial W
  %

  covfunc_w = cell(D,1);
  theta_w = cell(D,1);
  is_pseudos_w = false(D,1);

  % Remove land area grid points
  in_w = data.coordinates;

  %% Smooth components (using pseudo inputs)

  % Pseudo inputs (uniformly with respect to area size)
  pseudo_w = points_on_sphere(18); % uniform points by number of latitudes
                                   % Remove pseudo inputs that are on land (the nearest grid point is land)
  ind_pseudo_w = mohsst5_points_to_grid_index(pseudo_w);
  mask = metoffice_get_mask(maskfile);
  pseudo_w(:,~mask(ind_pseudo_w)) = [];
% $$$ % This code shows the pseudo inputs on the map
% $$$ figure
% $$$ map_projection('global-ellipse');
% $$$ map_plot(pseudo_w,'r+');
% $$$ map_coast()
% $$$ map_grid()
% $$$ return

  % Transform inputs to 3-D Euclidean coordinates
  in_w = geographic_to_euclidean(in_w);
  pseudo_w = geographic_to_euclidean(pseudo_w);

  % Squared distance matrices for the covariance functions
  D2_ww = sq_dist(in_w);
  D2_pp = sq_dist(pseudo_w);
  D2_pw = sq_dist(pseudo_w, in_w);
  d2_w = diag(D2_ww);

  ind = 0;

  if comps_spatial(1) > 0
    ind = ind(end) + (1:comps_spatial(1));
    fprintf('%d iid components for W\n', length(ind));

    covfunc_w(ind) = {gp_cov_scale(gp_cov_delta(size(in_w,2)))};
    theta_w(ind) = columns_to_cells(...
        [linspace(1,0.1,length(ind))]);       % magnitudes
  end

  if comps_spatial(2) > 0
    ind = ind(end) + (1:comps_spatial(2));
    fprintf('%d slow components for W (using %d pseudo inputs)\n', length(ind), ...
            size(pseudo_w,2));

    % Covariance function (scaled squared exponential) with pseudo inputs
    covfunc = @(D2) gp_cov_se(D2);
    covfunc_w(ind) = {gp_cov_pseudo(...
        gp_cov_scale(gp_cov_jitter(covfunc(D2_pp), 1e-3)), ...
        gp_cov_scale(covfunc(D2_pw)), ...
        gp_cov_scale(covfunc(d2_w)))};

    % Hyperparameters for the covariance functions
    theta_w(ind) = columns_to_cells(...
        [linspace(1,0.1,length(ind));       % magnitudes
         linspace(4000,1000,length(ind))]); % lengthscales
    is_pseudos_w(ind) = true;
  end

  %
  % Process data
  %

  % Filename for saving the results
  filename = sprintf('%s/results_rectest_%s_gpfa_D=%d_anomalies=%d_%s', ...
                     folder, ...
                     dataset, ...
                     D, ...
                     anomalies, ...
                     date);

  % Component-wise factorization for X
  X_module = factor_module_gp_factorized(N, covfunc_x, theta_x, ...
                                         'update_hyperparameters', [5 10:10:100 100:100:2000], ...
                                         'maxiter_hyperparameters', 5, ...
                                         'is_pseudo', is_pseudos_x, ...
                                         'init', zeros(D,N));

  % Component-wise factorization for W
  W_module = factor_module_gp_factorized(M, covfunc_w, theta_w, ...
                                         'update_hyperparameters', [5 10:10:100 100:100:2000], ...
                                         'maxiter_hyperparameters', 5, ...
                                         'is_pseudo', is_pseudos_w);

  % Isotropic noise (precisions weighted proportionally to grid size)
% $$$ figure
% $$$ mohsst5_mapplot(metoffice_add_land(weights(:,1)));
% $$$ return
  noise_module = noise_module_isotropic(M, N, 1e-3, 1e-3, ...
                                        'init', 10, ...
                                        'weights', weights);

  % Run GPFA
  Q = gpfa(D, Y, W_module, X_module, noise_module, ...
           'maxiter', maxiter, ...
           'rotate', 1:100, ... %[1:50 60:10:2000], ...
           'autosavefile', filename, ...
           'autosave', [10:100:2000]);

  Yh = Q.W'*Q.X;

  % Some performance measures
  fprintf('Weighted training RMSE of the reconstruction: %f\n',  ...
          rmsew(Y(Obs)-Yh(Obs),weights(Obs)));

  % Save the results
  save(filename, '-struct', 'Q');
  fprintf('Saved GPFA results to %s\n', filename);

else

  Q = load(gpfafile);

  fprintf('Loaded GPFA results from %s\n', gpfafile);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Short-scale GP inference for the residuals
%

% Reconstruct
if debug
  Yh = zeros(size(Y));
else
  Yh = Q.W'*Q.X;
end
Yres = Y - Yh;
clear Q;

%
% Temporal covariance function (assume uniformly spaced time instances)
%

d = abs(1-(1:length(data.time)));
covfunc1 = gp_cov_toeplitz(gp_cov_pp(d,1));
theta_temporal = [7];
% $$$ covfunc1 = gp_cov_toeplitz(gp_cov_sum(gp_cov_pp(d,1), ...
% $$$                                       gp_cov_scale(gp_cov_pp(d,1))));
% $$$ theta_temporal = [7;   % length scale 1
% $$$                   1.0; % magnitude 2
% $$$                   3];  % length scale 2

%
% Spatial covariance function
%

[LON,LAT] = meshgrid(data.lon,data.lat);

X = geographic_to_euclidean([LON(:)';LAT(:)']);

% Use block-Toeplitz structure for the covariance function
[lat,lon0] = meshgrid(data.lat,data.lon(1));
X0 = geographic_to_euclidean([lon0(:)';lat(:)']);
D = sqrt(sq_dist(X0,X));
covfunc2 = gp_cov_jitter(gp_cov_toeplitz_block(gp_cov_pp(D,3)), 1e-3);
theta_spatial = nan;
if datanum <= 4
  % 5x5 SST
  theta_spatial = [3000];  % length scale
elseif datanum == 5
  % 1x1 sea ice
  theta_spatial = [50];
else
  error('Unknown datanum');
end

% Select sea areas
sea = metoffice_get_mask(maskfile);
covfunc2 = gp_cov_select(covfunc2, sea);

% $$$ % DEBUG: Test solver stuff
% $$$ K2 = covfunc2(theta_spatial);
% $$$ S = solver_ldlchol();
% $$$ R = S.decompose(K2);
% $$$ L = S.squareroot(R);
% $$$ K2h = L*L';
% $$$ norm_err = sqrt(mean((K2(:)-K2h(:)).^2))
% $$$ L = lchol(K2);
% $$$ x = randn(length(K2),1);
% $$$ y = K2 * x;
% $$$ xh = S.linsolve(R, y);
% $$$ x_err = sqrt(mean((xh-x).^2))
% $$$ return

%
% Inference
%

burnin = floor(N_samples/2);
folder_samples = sprintf('%s/samples_rectest_%s_gpkron2_%s', ...
                         folder, ...
                         dataset, ...
                         date);
mkdir(folder_samples);
filename = sprintf('%s/results_rectest_%s_gpkron2_D=%d_anomalies=%d_%s', ...
                   folder, ...
                   dataset, ...
                   sum(comps), ...
                   anomalies, ...
                   date);
filename_samples = sprintf('%s/samples_rectest_%s_gpkron2_D=%d_anomalies=%d_%s', ...
                           folder_samples, ...
                           dataset, ...
                           sum(comps), ...
                           anomalies, ...
                           date);

a = 1e-3;
b = 1e-3;
logprior_theta = @(theta) sum(gamma_logpdf(theta, a, b));
dlogprior_theta = @(theta) gamma_dlogpdf(theta, a, b);

% Initial guess for covariance parameters
theta_init = [0.5; ...               % total magnitude
              theta_temporal(:); ... % temporal parameters
              0.5; ...               % temporal noise magnitude
              theta_spatial(:); ...  % spatial parameters
              0.5]';                 % spatial noise magnitude

samplefunc = get_sample_function2(numel(theta_init), N_samples, burnin, ...
                                                filename, filename_samples);

% Weights for the noise levels using the respective grid size
w = 1./sqrt(cosd(LAT));
w = w(sea);
%W = repmat(w(sea), [1,size(Y,2)]);

[get_logpdf, get_dlogpdf, get_rand_model, func_theta] = ...
    gp_init_kron1(covfunc1, ...
                  solver_ldlchol(), ...
                  covfunc2, ...
                  solver_ldlchol(), ...
                  logprior_theta, ...
                  dlogprior_theta, ...
                  samplefunc, ...
                  'noise_scale2', w);

% Transform to log-scale
func_theta = @(logtheta, varargin) func_theta_transformed(logtheta, ...
                                                  exp(logtheta), ...
                                                  func_theta, ...
                                                  varargin{:});
get_logpdf = @(f_theta) get_logpdf_transformed(f_theta, ...
                                               get_logpdf, ...
                                               sum(f_theta.theta_transformed));
get_dlogpdf = @(df_theta) get_dlogpdf_transformed(df_theta, ...
                                                  get_dlogpdf, ...
                                                  diag(exp(df_theta.theta_transformed)), ...
                                                  ones(size(df_theta.theta_transformed)));
theta_init = log(theta_init);

[rand_theta, f_theta] = mcmc_init_slicesampling(theta_init, ...
                                                get_logpdf, ...
                                                'fx', func_theta);

rand_model = get_rand_model(rand_theta, f_theta);

% Gibbs sampling
tic
gibbs(Yres,~Obs, N_samples, rand_model);
toc

% Check the results
res = samplefunc();

save(filename, '-struct', 'res');
disp(['Saved short-scale GP results to ', filename]);

% Mean reconstruction
F = Yh + res.F;

filename = sprintf('%s/results_rectest_%s_gpfa_gpkron2_D=%d_anomalies=%d_%s', ...
                   folder, ...
                   dataset, ...
                   sum(comps), ...
                   anomalies, ...
                   date);
save(filename, 'F');
disp(['Saved total reconstruction to ', filename]);

fprintf('Weighted training RMSE of the reconstruction: %f\n',  ...
        rmsew(Y(Obs)-F(Obs),weights(Obs)));

if nargout < 1
  clear F;
end

end

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

function S = solver_ldlchol()

S.decompose = @decompose;
S.linsolve = @linsolve;
S.logdet = @logdet;
S.inv = @inv;
S.squareroot = @squareroot;

  function S = decompose(K)
  [S.LD,p,S.q] = ldlchol(K);
  if p>0
    error('Matrix must be positive definite');
  end
  end

  function x = linsolve(S, y)
  x(S.q,:) = linsolve_ldlchol(S.LD,y(S.q,:));
  end

  function ldet = logdet(S)
  ldet = logdet_ldlchol(S.LD);
  end

  function A = inv(S)
  %A(S.q,S.q) = spinv_ldlchol(S.LD);
  N = length(S.q);
  R = sparse(1:N,S.q,ones(N,1));
  L = R'*spinv_ldlchol(S.LD)*R;
  end

  function L = squareroot(S)
  % WARNING/TODO: This permutation might take A LOT of time..
  % L(S.q,S.q) = ldlchol2lchol(S.LD);
  % Try instead:
  N = length(S.q);
  R = sparse(1:N,S.q,ones(N,1));
  L = R'*ldlchol2lchol(S.LD)*R;
  end

end

function [f_theta, df_theta] = func_theta_transformed(theta_transformed, ...
                                                  theta, func_theta, varargin)
if nargout <= 1
  f_theta = func_theta(theta, varargin{:});
  f_theta.theta_transformed = theta_transformed;
else
  [f_theta, df_theta] = func_theta(theta, varargin{:});
  f_theta.theta_transformed = theta_transformed;
  df_theta.theta_transformed = theta_transformed;
end
end

function logpdf_y = get_logpdf_transformed(fy, get_logpdf, logjacobian)
logpdf = get_logpdf(fy);
logpdf_y = @logpdf_transformed;
  function lpdf = logpdf_transformed(varargin)
  lpdf = logpdf(varargin{:}) + logjacobian;
  end
end

function dlogpdf_y = get_dlogpdf_transformed(dfy, get_dlogpdf, Jacobian, ...
                                                  dlogjacobian)
dlogpdf = get_dlogpdf(dfy);
dlogpdf_y = @dlogpdf_transformed;
  function dlpdf = dlogpdf_transformed(varargin)
  dlpdf = dlogpdf(varargin{:});
  dlpdf = Jacobian*dlpdf + dlogjacobian;
  end
end

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

function samplefunc = get_sample_function2(D_theta, N, burnin, filename, filename_samples)
results.Y = [];
results.F = 0;
results.FF = 0;
results.theta = zeros(D_theta, N);
samplefunc = @process_sample;
n = 1;
  function res = process_sample(Y, F, theta)
  if nargin >= 1
    % Store results
    results.Y = Y;
    if true && n > burnin
      results.F = (F + (n-burnin-1)*results.F) / (n-burnin);
    end
    if true && n > burnin
      results.FF = (F.*F + (n-burnin-1)*results.FF) / (n-burnin);
    end
    results.theta(:,n) = theta(:);
    % Save results
    fprintf('Saving results to %s..', filename)
    save(filename, '-struct', 'results');
    save(sprintf('%s_F%d',filename_samples,n), 'F');
    fprintf(' done.\n')
    n = n + 1;
  end
  if nargout >= 1
    % Return results
    res = results;
  end
  end
end

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

function gibbs(Y, Imv, N, rand_model, samplefunc)

Y(Imv) = 0;

for n=1:N

  t = cputime();
  Y = rand_model(Y,Imv);
  dt = cputime() - t;
  fprintf('Iteration step %d done. (%f seconds)\n', n, dt)

end

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function polinomio=trigpol(N)

%  TRIGPOL generate trigonometric polynomial.
%
%          TRIGPOL(N) generates the following polynomial in e^jw:
%
%          P(e^jw)=sum(k=0,...,N-1){ (N-1+k) ( e^(-jw) - 2 + e^jw )
%                                       k
%
%          The output of the function is a vector of size 2*N-1 holding
%          the coefficients corresponding to the different powers of e^jw.
%
%          See also: NUMCOMB, DAUB, SPLINE

%--------------------------------------------------------
% Copyright (C) 1994, 1995, 1996, by Universidad de Vigo
%
%
% Uvi_Wave is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2, or (at your option) any
% later version.
%
% Uvi_Wave 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 Uvi_Wave; see the file COPYING.  If not, write to the Free
% Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
%
%       Author: Nuria Gonzalez Prelcic
%       e-mail: Uvi_Wave@tsc.uvigo.es
%--------------------------------------------------------

%The polynomial is constructed from a sum.
%coefs holds the coefficients of each term in the sum.

coefs=zeros(N,2*N-1);
coefs(1,N)=1;

for i=1:N-1
	fila=[1 -2 1];
	for j=2:i
		fila=conv(fila,[1 -2 1]);
	end;
	fila=numcomb(N-1+i,i)*(-0.25)^i*fila;
	fila=[ zeros(1,(N-i-1))  fila zeros(1,(N-i-1))];
	coefs(i+1,:)=fila;
end

for i=0:(2*(N-1))
	polinomio(i+1)=0;
	for j=1:N
		polinomio(i+1)=polinomio(i+1)+coefs(j,i+1);
	end
end;

author: alexandrebarachant


text:
function [yi, ypi, yppi] = qhermite(x, y, yp, ypp, xi, v)

% QHERMITE 1-D piecewise quintic Hermite spline
%    QHERMITE(X,Y,YP,YPP,XI,M) interpolates to find YI, the values of
%    the underlying function Y at the points in the array XI, using
%    piecewise quintic Hermite splines.  X and Y must be vectors
%    of length N.
%
%    V specifies how tangents are calculated when derivatives are not
%    specified.  V can be:
%       0 : Finite difference (default)
%       1 : Catmull-Rom spline
%
%    [YI,YPI,YPPI] = QHERMITE() also returns the interpolated
%    quartic derivative and cubic second derivative of the underlying
%    function Y at points XI.

% Joe Henning - Fall 2011

if (nargin < 6)
   v = 0;
end

if (v == 0)
   % precompute finite difference derivatives
   if (isempty(yp) && isempty(ypp))
      m = 0;
      yp = lfindiff (x, y);
      ypp = lfindiff (x, yp);
   elseif (isempty(yp))
      m = 1;
      yp = lfindiff (x, y);
   elseif (isempty(ypp))
      m = 2;
      ypp = lfindiff (x, yp);
   else
      m = -1;
   end
else
   if (isempty(yp) && isempty(ypp))
      m = 0;
   elseif (isempty(yp))
      m = 1;
   elseif (isempty(ypp))
      m = 2;
   else
      m = -1;
   end
end

n = length(x);

for i = 1:length(xi)
   % Find the right place in the table by means of a bisection.
   klo = 1;
   khi = n;
   while (khi-klo > 1)
      k = fix((khi+klo)/2.0);
      if (x(k) > xi(i))
         khi = k;
      else
         klo = k;
      end
   end

   h = x(khi) - x(klo);
   if (h == 0.0)
      fprintf('??? Bad x input to qhermite ==> x values must be distinct\n');
      yi(i) = NaN;
      ypi(i) = NaN;
      yppi(i) = NaN;
      continue;
   end

   if (m == -1)
      a = yp(klo);
      b = yp(khi);
      c = ypp(klo);
      d = ypp(khi);
   elseif (v == 0)
      % Finite difference
      a = yp(klo);
      b = yp(khi);
      c = ypp(klo);
      d = ypp(khi);
   elseif (m == 0)
      if (klo == 1)
         a = (y(khi) - y(klo))/h;
         b = (y(khi+1) - y(klo))/(x(khi+1) - x(klo));
         c = ((y(khi+1) - y(khi))/(x(khi+1) - x(khi)) - (y(khi) - y(klo))/h)/...
             (x(khi+1) - x(klo));
         d = 2*(y(khi+1) - 2*y(khi) + y(klo) - (x(khi+1) - 2*x(khi) + x(klo))*b)/...
             ((x(khi+1) - x(khi))*(x(khi+1) - x(khi)) + h*h);
      elseif (khi == n)
         a = (y(khi) - y(klo-1))/(x(khi) - x(klo-1));
         b = (y(khi) - y(klo))/h;
         c = 2*(y(khi) - 2*y(klo) + y(klo-1) - (x(khi) - 2*x(klo) + x(klo-1))*a)/...
             (h*h + (x(klo) - x(klo-1))*(x(klo) - x(klo-1)));
         d = ((y(khi) - y(klo))/h - (y(klo) - y(klo-1))/(x(klo) - x(klo-1)))/...
             (x(khi) - x(klo-1));
      else
         a = (y(khi) - y(klo-1))/(x(khi) - x(klo-1));
         b = (y(khi+1) - y(klo))/(x(khi+1) - x(klo));
         c = 2*(y(khi) - 2*y(klo) + y(klo-1) - (x(khi) - 2*x(klo) + x(klo-1))*a)/...
             (h*h + (x(klo) - x(klo-1))*(x(klo) - x(klo-1)));
         d = 2*(y(khi+1) - 2*y(khi) + y(klo) - (x(khi+1) - 2*x(khi) + x(klo))*b)/...
             ((x(khi+1) - x(khi))*(x(khi+1) - x(khi)) + h*h);
      end
   elseif (m == 1)
      c = ypp(klo);
      d = ypp(khi);
      if (klo == 1)
         a = (y(khi) - y(klo))/h;
         b = (y(khi+1) - y(klo))/(x(khi+1) - x(klo));
      elseif (khi == n)
         a = (y(khi) - y(klo-1))/(x(khi) - x(klo-1));
         b = (y(khi) - y(klo))/h;
      else
         a = (y(khi) - y(klo-1))/(x(khi) - x(klo-1));
         b = (y(khi+1) - y(klo))/(x(khi+1) - x(klo));
      end
   else
      a = yp(klo);
      b = yp(khi);
      if (klo == 1)
         c = ((y(khi+1) - y(khi))/(x(khi+1) - x(khi)) - (y(khi) - y(klo))/h)/...
             (x(khi+1) - x(klo));
         d = 2*(y(khi+1) - 2*y(khi) + y(klo) - (x(khi+1) - 2*x(khi) + x(klo))*b)/...
             ((x(khi+1) - x(khi))*(x(khi+1) - x(khi)) + h*h);
      elseif (khi == n)
         c = 2*(y(khi) - 2*y(klo) + y(klo-1) - (x(khi) - 2*x(klo) + x(klo-1))*a)/...
             (h*h + (x(klo) - x(klo-1))*(x(klo) - x(klo-1)));
         d = ((y(khi) - y(klo))/h - (y(klo) - y(klo-1))/(x(klo) - x(klo-1)))/...
             (x(khi) - x(klo-1));
      else
         c = 2*(y(khi) - 2*y(klo) + y(klo-1) - (x(khi) - 2*x(klo) + x(klo-1))*a)/...
             (h*h + (x(klo) - x(klo-1))*(x(klo) - x(klo-1)));
         d = 2*(y(khi+1) - 2*y(khi) + y(klo) - (x(khi+1) - 2*x(khi) + x(klo))*b)/...
             ((x(khi+1) - x(khi))*(x(khi+1) - x(khi)) + h*h);
      end
   end

   % Evaluate quintic Hermite polynomial
   t = (xi(i) - x(klo))/h;
   t2 = t*t;
   t3 = t2*t;
   t4 = t3*t;
   t5 = t4*t;
   h2 = h*h;
   z0 = -6*t5 + 15*t4 - 10*t3 + 1;
   z1 = -3*t5 + 8*t4 - 6*t3 + t;
   z2 = 0.5*(-t5 + 3*t4 - 3*t3 + t2);
   z3 = 1 - z0;
   z4 = -3*t5 + 7*t4 - 4*t3;
   z5 = 0.5*(t5 - 2*t4 + t3);

   yi(i) = z0*y(klo) + z1*h*a + z2*h2*c + z3*y(khi) + z4*h*b + z5*h2*d;

   % Differentiate to find the fourth-order interpolant
   z0 = -30*t4 + 60*t3 - 30*t2;
   z1 = -15*t4 + 32*t3 - 18*t2 + 1;
   z2 = 0.5*(-5*t4 + 12*t3 - 9*t2 + 2*t);
   z3 = -z0;
   z4 = -15*t4 + 28*t3 - 12*t2;
   z5 = 0.5*(5*t4 - 8*t3 + 3*t2);

   ypi(i) = z0*y(klo)/h + z1*a + z2*h*c + z3*y(khi)/h + z4*b + z5*h*d;

   % Differentiate to find the third-order interpolant
   z0 = -120*t3 + 180*t2 - 60*t;
   z1 = -60*t3 + 96*t2 - 36*t;
   z2 = 0.5*(-20*t3 + 36*t2 - 18*t + 2);
   z3 = -z0;
   z4 = -60*t3 + 84*t2 - 24*t;
   z5 = 0.5*(20*t3 - 24*t2 + 6*t);

   yppi(i) = z0*y(klo)/h2 + z1*a/h + z2*c + z3*y(khi)/h2 + z4*b/h + z5*d;
end

% 3-pt finite difference
function [yp] = lfindiff(x, y)
% Finite Difference Formulae for Unequal Sub-Intervals Using Lagrange's Interpolation Formula
% Singh, Ashok K. and B. S. Bhadauria
% Int. Journal of Math. Analysis, Vol. 3, 2009, no. 17, 815-827

n = length(x);

if (n ~= length(y))
   fprintf('Error, len(x) != len(y), returning\n');
   yp = NaN;
   return;
end

for k = 1:n
   if (k == 1)
      h1 = x(k+1) - x(k);
      h2 = x(k+2) - x(k+1);
      yp(k) = -(2*h1+h2)/(h1*(h1+h2))*y(k) + (h1+h2)/(h1*h2)*y(k+1) - h1/(h2*(h1+h2))*y(k+2);
   elseif (k == n)
      h1 = x(k-1) - x(k-2);
      h2 = x(k) - x(k-1);
      yp(k) = h2/(h1*(h1+h2))*y(k-2) - (h1+h2)/(h1*h2)*y(k-1) + (h1+2*h2)/(h2*(h1+h2))*y(k);
   else
      h1 = x(k) - x(k-1);
      h2 = x(k+1) - x(k);
      yp(k) = -h2/(h1*(h1+h2))*y(k-1) - (h1-h2)/(h1*h2)*y(k) + h1/(h2*(h1+h2))*y(k+1);
   end
end


text:
function [U,H] = poldec(f)
%POLDEC   Polar decomposition of a CHEBFUN2.
%   [U, H] = POLDEC(F) computes chebfun2 objects U and H such that F = U*H.
%   The domain of U is the same as that of F and it is a partial isometry,
%   which means all its singular values are 1 (a finite number) or 0.
%   The domain of H is the square [a,b]x[a,b] where [a,b] is the
%   x-domain of F, and H is Hermitian positive semidefinite (its
%   eigenvalues are all positive).
%
% See also SVD.

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

[u, S, v] = svd(f); % SVD of chebfun2
U = u*v';
H = v*S*v';

if ( nargout > 1 )
    varargout = { U, H };
else
    varargout = { U };
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [nph,nfh] = sbprogress(varargin)
%SBPROGRESS Creates a progress bar in the "status-bar section" of the figure
%[nph,nfh] = sbprogress(varargin)
%   sbprogress creates or updates a progress bar that is created at the
%   status-bar section of a figure. It can either be created on an existing
%   figure or raised along with a newly created figure. At least 2
%   parameters are
%   required (or zero for demo).
%
%   Syntax:
%       sbprogres
%       sbprogress(ph,x)
%       sbprogress([],x)
%       sbprogress(ph,'Close')
%       sbprogress(ph,'Redraw')
%       sbprogress(ph,x,'PropertyName',PropertyValue,...)
%       sbprogress(ph,'PropertyName',PropertyValue,...)
%       sbprogress([],x,'PropertyName',PropertyValue,...)
%       sbprogress([],'PropertyName',PropertyValue,...)
%       [nph,nfh] = sbprogress(...)
%
%   Description:
%       sbprogress runs a demo
%       sbprogress([],x) creates a NEW figure and a NEW progress bar on it with progress being set to x
%       sbprogress(ph,x) updates the progress bar at handles ph with progress being set/changed to x; this also redraws the progress bar
%       sbprogress(ph,'Close') closes an existing progress bar and also a figure if figure was created together with a progress bar otherwise only the progress bar will be removed from the figure
%       sbprogress(ph,'Redraw') redraws an existing progress bar; e.g., when figure is resized
%       sbprogress(...'PropertyName',PropertyValue,...) sets/resets some properties of the progress bar and the newly cretaed figure
%       [nph,nfh] = sbprogress(...) returns handles to newly created progress bar and figure if figure is created along with the progress bar
%
%   Properties:
%       'Parent':   a figure handle on which the progress is to be created
%                   (if empty or not given a new figure will also be raised)
%       'Msg':      a string message that is displayed on the left of a progress bar
%       'Color':    a color of the progress bar and the percentage text (defaults to blue)
%       'Visible':  'on'/'off'; shows/hides the progress bar and its bevel; if a figure is also created then this affect the figure too
%       'Title':    a string title that is displayed in the figure when the
%                   progress bar is created in a new figure; defaults to an empty string
%       'Size':     height of the progress bar; 1..12 (defaults to 1)
%       'Modal':    1 for newly created figure to be modal and 0 to be
%                   non-modal (default)
%       'DisplayMode': 1 (default) percentage and the remaining time are added
%                   to the title string (only when a new figure is raised);
%                   0 - no time info is displayed; 2 - the remaining time is
%                   added to the progress bar (shown only when x > 0);
%                   creating the sbprogress and/or setting the progress (x)
%                   to zero resets the timer
%
%   Notes:
%       X can be a value between -100..100; given x as a negative value
%       means that the percentage number will not be shown.
%       'Modal' and 'Title' are considered only in cases when a new figure is created alogn with the
%       progress bar.
%
%   EXAMPLES:
%   sbprogress
%           Runs a demo
%   sbprogress([],10,'Msg','Progress:','Title','Test progress')
%           Creates a new figure with a progress bar...
%   ph = sbprogress([],-10,'Msg','Progress:','Title','Test progress')
%           Creates a new figure with a progress bar where the percentage
%           number is not shown; a handle array to the progress-bar axis is also returned
%   sbprogress(ph,'Redraw')
%           Redraws the ph progress bar
%   sbprogress(ph,20)
%           Re-sets the progress-bar's progress
%   ph = sbprogress([],20,'Msg','Test','Parent',fh)
%           Creates a progress bar on an existing figure whose handle is fh
%           and returns its handle array for further manipulation
%           (updating, changing,...)
%   set(fh,'UserData',ph,'ResizeFcn','sbprogress(get(gcf,''UserData''),''redraw'')')
%           Sets figure's (the one with the ph handle) resize function to
%           resize progress bar automatically
%
%   See also: STSBAR, ASCIIPROGRESS
%
%   Primoz Cermelj, 08.12.2003
%   (c) Primoz Cermelj, primoz.cermelj@email.si
%
%   Version: 1.1.0
%   Last revision: 19.02.2005
%--------------------------------------------------------------------------

%----------------
% SBPROGRESS history
%----------------
%
% [v.1.1.0] 19.02.2005
% - NEW: DisplayMode option added to display the remaining time
% - FIX: A bug relating to the Size property fixed
%
% [v.1.0.1] 05.09.2004
% - FIX: Minor update for error-checking
%
% [v.1.0.0] 28.05.2004
% - NEW: First official release

global W border

W = 400;        % Width of a newly created figure
border = 3;

% ---------------------
% First parameter-check and demo run
% ---------------------
if nargin == 0
    % Run a demo
    fh = figure;
    uicontrol('Style','text','String','Try to resize the figure...','Position',[150,200,200,25]);
    [nph,nfh] = sbprogress([],48.5,'Msg','Test:','Title','SBPROGRESS demo','Parent',fh);
    set(fh,'UserData',nph,'ResizeFcn','sbprogress(get(gcf,''UserData''),''redraw'')');
    [nph2,nfh2] = sbprogress([],-0.001,'Msg','Test 2:','Title','SBPROGRESS demo 2','Modal',0,'DisplayMode',1,'Color','r','Size',10);
    for ii=1:10
        sbprogress(nph2,-100*ii/10);
        pause(0.2);
    end
    return
elseif nargin < 2
    error('At least progress length should be given or ''Close'' or ''Redraw'' string');
end

ph = varargin{1};
x = varargin{2};
if ~ishandle(ph)
    error('ph is not a handle');
end

% Default properties
Props.Parent = [];
Props.Msg = '';
Props.Color = 'b';
Props.Title = '';
Props.Visible = 'on';
Props.DisplayMode = 1;
Props.Modal = 0;
Props.Size = 23;     % Default height of a newly created figure; this is aslo related to a progress bar height
propChanged = 0;     % if any property is being changed
phExists = ~isempty(ph);
fhExists = 1;
if phExists
    Data = get(ph,'UserData');
    if ~isfield(Data,'Props')
        error('Props field was not found in progress''s user data. Something was either removed or improper progress-bar handle was given.');
    else
        %%%%%%%%%%%%%%
        Props = Data.Props; % read props from current props
        %%%%%%%%%%%%%%
    end
end

% ---------------------
% Either 'Redraw' or 'Close' is assumed given in the second parameter x
% ---------------------
if ischar(x) & (strcmpi(x,'close') | strcmpi(x,'redraw'))
    if strcmpi(x,'close')
        Data = get(ph,'UserData');
        if ~isempty(Data) && isfield(Data,'fhCreated') && isfield(Data,'fh')
            if Data.fhCreated
                delete(Data.fh);
            else
                delete(ph);
            end
        else
            error('Wrong ph handle given');
        end
    elseif strcmpi(x,'redraw')
        Data = get(ph,'UserData');
        if ~isempty(Data) && isfield(Data,'fhCreated') && isfield(Data,'fh')
            progress(ph,Data.x,Props,propChanged,Data.fh,Data.fhCreated,'update');
        else
            error('Wrong ph handle given');
        end
    else
        error('Use sbprogress(ph,''redraw'') or sbprogress(ph,''close'')');
    end
    return
end

% ---------------------
% x might be ommited, empty, numeric value or a property name
% Sets any default values and any given ones given under properties
% ---------------------
if ischar(x)    % x is assumed being ommited (property is given)
    start = 2;
else
    start = 3;
end
if ~isnumeric(x) | isempty(x)
    if phExists
        x = Data.x;
    else
        x = 0;
    end
end
if nargin > 2       % additional properties are being passed
    nProps = nargin - start + 1;
    propChanged = 1;
    if ~logical( 2*floor(nProps/2)==nProps ), error('Additional properties must be given in pairs'), end;
    for ii=start:2:start+nProps-2
        if strcmpi(varargin{ii},'parent')
            Props.Parent = varargin{ii+1};
        elseif strcmpi(varargin{ii},'msg')
            Props.Msg = varargin{ii+1};
        elseif strcmpi(varargin{ii},'color')
            Props.Color = varargin{ii+1};
        elseif strcmpi(varargin{ii},'title')
            Props.Title = varargin{ii+1};
        elseif strcmpi(varargin{ii},'visible')
            Props.Visible = varargin{ii+1};
        elseif strcmpi(varargin{ii},'modal')
            Props.Modal = logical(varargin{ii+1});
        elseif strcmpi(varargin{ii},'size')
            Props.Size = varargin{ii+1};
        elseif strcmpi(varargin{ii},'displaymode')
            Props.DisplayMode = varargin{ii+1};
        else
            error('Unknown property given');
        end
    end
    if (Props.Size < 1) | (Props.Size > 12)
        Props.Size = 23;
    else
        Props.Size = round(Props.Size) + 22;
    end
end

if isempty(Props.Parent) || ~ishandle(Props.Parent)
    fhExists = 0;
else
    if ~ishandle(ph)
        phExists = 0;
    end
end
if ~phExists    % new progress bar is to be created
    fhCreated = 0;
    if ~fhExists
       nfh = createnewfig(Props);
       fh = nfh;
       Props.Parent = fh;
       fhCreated = 1;
    else
       fh = Props.Parent;
    end
    nph = progress([],x,Props,1,fh,fhCreated,'new');
    ph = nph;
else
    progress(ph,x,Props,propChanged,Data.fh,Data.fhCreated,'update');
end

if nargout > 0
    nph = ph;
    if nargout > 1
        nfh = fh;
    end
end
%-------------------------------------

%-------------------------------------
function fh = createnewfig(Props)
global W
H = Props.Size;
if isempty(Props.Title); Props.Title = ''; end;
set(0,'Units','Pixels');
scr = get(0,'ScreenSize');
fh = figure('Visible','off');
set(fh,'Units','Pixels',...
       'Tag','progressBarFigure',...
       'Position',[scr(3)/2-W/2 scr(4)/2-H/2 W H],...
       'MenuBar','none',...
       'Resize','off',...
       'NumberTitle','off',...
       'Name',Props.Title);
if Props.Modal
    set(fh,'WindowStyle','modal');
end
if strcmpi(Props.Visible,'on')
    set(fh,'Visible','On');
end
%-------------------------------------

%-------------------------------------
function [ph,hMsg] = progress(ph,x,Props,propChanged,fh,fhCreated,state)
% Creates a new progress or updates an existing one
global border

H = Props.Size;
[fPos,pPos] = getcoords(fh,H);
w = pPos(3);
if x < -100; x = 0; end;
if x > 100; x = 100; end;
if propChanged
    if isempty(Props.Msg); Props.Msg = ''; end;
    if isempty(Props.Title); Props.Title = ''; end;
    if isempty(Props.Msg)
        Data.xMsgPos = 0;
    else
        Data.xMsgPos = 2*border;
    end
end
timeStr = '';
if nargin < 5 || strcmpi(state,'new')
    Data.startTime = clock;
    if Props.DisplayMode ~= 0
        timeStr = '?? remaining';
    end
    ph = axes;
    set(ph,'Tag','progressBarAxes','Units','Pixels','Position',pPos,'Box','off','Visible','off');
    Data.hMsg = text(Data.xMsgPos,0.55*H,Props.Msg,'Units','Pixels','Tag','progressBarMsg','Visible','off');
else            % update
    Data = get(ph,'UserData');
    if x == 0
        Data.startTime = clock;
    end
    if Props.DisplayMode ~= 0
        remTime = (100/abs(x+0.00000001)-1)*etime(clock,Data.startTime);
        timeStr = sec2timestr(abs(remTime));
    end
    fh = Data.fh;
    set(ph,'Position',pPos);
    set(Data.hMsg,'String',Props.Msg,'Position',[Data.xMsgPos 0.55*H]);
end
%-----
set(ph,'XLim',[0 w],'YLim',[0 H]);
%-----
tPos = get(Data.hMsg,'Extent');
if (tPos(3)+tPos(1)) > 0.5*w;  pStart = 0.5*w; else;  pStart = (tPos(3)+tPos(1)+2*border); end;
pEnd = pStart + abs(x)/100*abs(w-border-pStart);
%-----
if nargin < 5 || strcmpi(state,'new')
    Data.hRec = rectangle('Position',[pStart border w-border-pStart-1 H-2*border-1],'Tag','progressBarRec','Visible','off');
    Data.hPatch = patch([pStart; pEnd; pEnd; pStart],[border; border; H-border-1; H-border-1],Props.Color,'Tag','progressBarPatch','Visible','off');
    if Props.DisplayMode == 2
        addStr = ['   ' timeStr];
    else
        addStr = '';
    end
    Data.hPercTxt = text( (pStart+w-border-1)/2,0.55*H,[num2str(abs(x),'%2.2f') '%' addStr],'Color',Props.Color,'Units','Pixels','FontWeight','bold','Tag','progressBarPercTxt','Visible','Off','EraseMode','xor','HorizontalAlignment','Center');
    Data.hFrame = panel(fh,ph,'new');
else
    set(Data.hRec,'Position',[pStart border w-border-pStart-1 H-2*border-1]);
    set(Data.hPatch,'XData',[pStart; pEnd; pEnd; pStart],'YData',[border; border; H-border-1; H-border-1],'FaceColor',Props.Color);
    if Props.DisplayMode == 2
        addStr = ['   ' timeStr];
    else
        addStr = '';
    end
    set(Data.hPercTxt,'Position',[(pStart+w-border-1)/2 0.55*H],'String',[num2str(abs(x),'%2.2f') '%' addStr],'Color',Props.Color);
    Data.hFrame = panel(fh,ph,'update');
end
%-----
if fhCreated
    fPos = get(fh,'Position');
    fPos(4) = H;
    if Props.DisplayMode == 1
        addStr = ['[' num2str(abs(x),'%2.2f') '%  ' timeStr ']'];
    else
        addStr = '';
    end
    set(fh,'Name',[Props.Title '  ' addStr],'Position',fPos);
    if propChanged
        if Props.Modal
            set(fh,'WindowStyle','modal');
        else
            set(fh,'WindowStyle','normal');
        end
        if strcmpi(Props.Visible,'on')
            set(fh,'Visible','on');
        else
            set(fh,'Visible','off');
        end
    end
end
%-----
if x >= 0 & strcmpi(Props.Visible,'on')
    set(Data.hPercTxt,'Visible','On');
else
    set(Data.hPercTxt,'Visible','Off');
end
if propChanged
    if strcmpi(Props.Visible,'on')
        set(Data.hRec,'Visible','On');
        set(Data.hMsg,'Visible','On');
        set(Data.hPatch,'Visible','On');
        set(Data.hFrame,'Visible','On');
    else
        set(Data.hRec,'Visible','Off');
        set(Data.hMsg,'Visible','Off');
        set(Data.hPatch,'Visible','Off');
        set(Data.hFrame,'Visible','Off');
    end
end
Data.Props = Props;
Data.pStart = pStart;
Data.fh = fh;
Data.fhCreated = fhCreated;
Data.x = x;
Data.msg = Props.Msg;
set(ph,'UserData',Data);
drawnow;
%-------------------------------------

%-------------------------------------
function frame = panel(fh,ph,state);
% Creates a virtual panel surrounding the progress bar.
% fh is the figure's handdle, ph is the progress' handle (axes).
% It returns a handle array designating the frame.
% See also the self-contained function BEVEL.
pos = get(ph,'Position');
x = pos(1);
y = pos(2);
w = pos(3);
h = pos(4);
col = rgb2hsv(get(fh,'Color'));
downColor = col;   downColor(2) = 0.5*downColor(2);  downColor(3) = 0.9; downColor = hsv2rgb(downColor);
upColor = col;    upColor(3) = 0.4;  upColor = hsv2rgb(upColor);
frame = zeros(4,1);
if nargin < 3 || strcmpi(state,'new')
    frame(1) = line([1 w],[h-1 h-1],[0 0],'Color',upColor,'Visible','off');
    frame(2) = line([1 1],[1 h-1],'Color',upColor,'Visible','off');
    frame(3) = line([2 w-1],[1 1],'Color',downColor,'Visible','off');
    frame(4) = line([w-1 w-1],[1 h-1],'Color',downColor,'Visible','off');
else
    Data = get(ph,'UserData');
    frame = Data.hFrame;
    set(frame(1),'XData',[1 w],'YData',[h-1 h-1]);
    set(frame(2),'XData',[1 1],'YData',[1 h-1]);
    set(frame(3),'XData',[2 w-1],'YData',[1 1]);
    set(frame(4),'XData',[w-1 w-1],'YData',[1 h-1]);
end
%-------------------------------------

%-------------------------------------
function [fPos,pPos] = getcoords(fh,H)
% Returns coordinates for axis, frame, and msg for our sbprogress (for
% creation or updating purposes). fPos is figure's position, pPos is
% progress-bar's position, framPos is frame's position
%
% Progress-bar axis
fUnits = get(fh,'Units');
set(fh,'Units','Pixels');
fPos = get(fh,'Position');
set(fh,'Units',fUnits);
pPos = fPos;
pPos(1) = 1;
pPos(2) = 2;
pPos(3) = pPos(3)-2;
pPos(4) = H-1;
%-------------------------------------

%-------------------------------------
function timeStr = sec2timestr(sec)
% Convert a time measurement from seconds into a human readable string.
dd = floor(sec/86400);
hh = floor((sec-dd*86400)/3600);
mm = floor((sec-dd*86400-hh*3600)/60);
ss = ceil(sec-dd*86400-hh*3600-mm*60);
timeStr = '';
if dd
    timeStr = sprintf('%d days, %d hr remaining',dd,hh);
elseif hh
    timeStr = sprintf('%d hr, %d min remaining',hh,mm);
elseif mm
    timeStr = sprintf('%d min, %d sec remaining',mm,ss);
else
    if sec < 0.01
        timeStr = 'completed';
    else
        timeStr = sprintf('%d seconds remaining',ss);
    end
end
%-------------------------------------


text:
classdef TestImageUtilities < CoreTest
    % TestImageUtilities. Tests for the PTKTextUtilities class.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    methods
        function obj = TestImageUtilities
            obj.TestCT;
        end
    end

    methods (Access = private)
        function TestCT(obj)
            datasets = {};
            datasets{end + 1} = obj.MakeDataset('OT', 12, 'OT0');
            datasets{end + 1} = obj.MakeDataset('CT', 12, 'CT1');
            datasets{end + 1} = obj.MakeDataset('CT', 14, 'CTlarge');
            datasets{end + 1} = obj.MakeDataset('CT', 1, 'CT2');
            datasets{end + 1} = obj.MakeDataset('MR', 13, 'MR1');
            datasets{end + 1} = obj.MakeDataset('OT', 25, 'OT1');
            obj.Assert(strcmp(MimImageUtilities.FindBestSeries(datasets), 'CTlarge'));
        end

        function TestMR(obj)
            datasets = {};
            datasets{end + 1} = obj.MakeDataset('OT', 12, 'OT0');
            datasets{end + 1} = obj.MakeDataset('CT', 12, 'CT1');
            datasets{end + 1} = obj.MakeDataset('CT', 14, 'CTlarge');
            datasets{end + 1} = obj.MakeDataset('CT', 1, 'CT2');
            datasets{end + 1} = obj.MakeDataset('MR', 13, 'MR1');
            datasets{end + 1} = obj.MakeDataset('MR', 100, 'MRlarge');
            datasets{end + 1} = obj.MakeDataset('OT', 250, 'OT1');
            obj.Assert(strcmp(MimImageUtilities.FindBestSeries(datasets), 'MRlarge'));
        end

        function TestOther(obj)
            datasets = {};
            datasets{end + 1} = obj.MakeDataset('AB', 12, 'CT1');
            datasets{end + 1} = obj.MakeDataset('CD', 14, 'CTlarge');
            datasets{end + 1} = obj.MakeDataset('EF', 1, 'CT2');
            datasets{end + 1} = obj.MakeDataset('GH', 1300, 'GHuid');
            datasets{end + 1} = obj.MakeDataset('IJ', 100, 'MRlarge');
            datasets{end + 1} = obj.MakeDataset('KL', 250, 'OT1');
            obj.Assert(strcmp(MimImageUtilities.FindBestSeries(datasets), 'GHuid'));
        end
    end

    methods (Access = private, Static)
        function fake_dataset = MakeDataset(modality, num_images, uid)
            fake_dataset = struct;
            fake_dataset.Modality = modality;
            fake_dataset.NumberOfImages = num_images;
            fake_dataset.SeriesUid = uid;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
domain = [-1,1,-1,1,-1,1];
bc = [1,1,1,1,1,1]; % Dirichelet BC

% epsm = 8.854*10^(-12);
% epsp = 11.68*8.854*10^(-12);
% sigm = 10^5;
% sigp = 10^6;
% mum = 4*pi*10^(-7);
% mup = 4*pi*10^(-7);
epsm = 8.85*10^(-2)*5;
epsp = 8.85*10^(-2);
sigm = 100;
sigp = 1;
mum = (4*pi)^(-1);
mup = (4*pi)^(-1)*2;
x0=0; y0=0; z0=-0.3; r1=0.2; r2=pi/5; omega = 20; a = 1; b= 150; intPt = -1;
pde = TorusTimeInitial2(mum,mup,sigm,sigp,epsm,epsp,omega,x0,y0,z0,r1,r2,a,b,intPt);
TEND = 0.1;
InitCond = 1;
BDCond = 1;

nx = 64;  h=(domain(2) - domain(1))/nx;
ny = nx;
nz = nx;

Ntime = 5*nx;%ceil(sqrt(nx));
deltaT = TEND/Ntime;

mesh = genMesh3D(domain, nx, ny, nz);
mesh = enrichMesh3D(mesh,2); % Mesh detail level = 1 (for IFE).
mesh = genIntfMesh3D(mesh,pde.intf);

Explevel = 3;
levelCount = 0;
BasicEdgeNum = size(mesh.e,1);
% level 1
if Explevel>0
    OldedgeID = unique(reshape(mesh.t_e(mesh.tLoc<0,:),[],1));
    BasicEdge = true(size(mesh.e,1),1); BasicEdge(OldedgeID) = false;
    BasicEdgeNum = sum(BasicEdge);
    TotalOldEdge = zeros(size(mesh.e,1),1);
    TotalOldEdge(BasicEdge) = 1:BasicEdgeNum;
    TotalOldEdge(OldedgeID) = BasicEdgeNum+1:size(mesh.e,1);
    mesh.t_e = TotalOldEdge(mesh.t_e);
    mesh.f_e = TotalOldEdge(mesh.f_e);
    [~,TotalOldEdgeSortID] = sort(TotalOldEdge);
    mesh.eLoc = mesh.eLoc(TotalOldEdgeSortID);
    mesh.e = mesh.e(TotalOldEdgeSortID,:);
    levelCount = levelCount + 1;
end
tId = mesh.tLoc<0;
while (levelCount<Explevel)
    fId = unique(reshape(mesh.t_f(tId,:),[],1));
    tId = unique(reshape(mesh.f_t(fId,:),[],1));
    tId = tId(tId>0);
    OldedgeID = unique(reshape(mesh.t_e(tId,:),[],1));
    BasicEdge = true(size(mesh.e,1),1); BasicEdge(OldedgeID) = false;
    BasicEdgeNum = sum(BasicEdge);
    TotalOldEdge = zeros(size(mesh.e,1),1);
    TotalOldEdge(BasicEdge) = 1:BasicEdgeNum;
    TotalOldEdge(OldedgeID) = BasicEdgeNum+1:size(mesh.e,1);
    mesh.t_e = TotalOldEdge(mesh.t_e);
    mesh.f_e = TotalOldEdge(mesh.f_e);
    [~,TotalOldEdgeSortID] = sort(TotalOldEdge);
    mesh.eLoc = mesh.eLoc(TotalOldEdgeSortID);
    mesh.e = mesh.e(TotalOldEdgeSortID,:);
    levelCount = levelCount + 1;
end

fem = genNedFEM3D(mesh,bc);
bm = epsm/deltaT^2 + sigm/(2*deltaT);
bp = epsp/deltaT^2 + sigp/(2*deltaT);
am = mum^(-1); ap = mup^(-1);
femI = genNed1IFEM3D(mesh,fem,bm,bp,am,ap);

S = globMatrixNedPGIFE3D(pde.Mu,1,mesh,femI,fem,fem);
Me = globMatrixNedPGIFE3D(pde.Epslon,0,mesh,femI,fem,fem);
Ms = globMatrixNedPGIFE3D(pde.Sig,0,mesh,femI,fem,fem);
Atotal = Me/deltaT^2 + Ms/(2*deltaT) + S/2;

Ndof = size(mesh.e,1);
bdidx = zeros(Ndof,1);
isBdEdge = true(Ndof,1);
isBdEdge(fem.mapper) = false;
bdidx(isBdEdge) = 1;
Tbd = spdiags(bdidx,0,Ndof,Ndof);
T = spdiags(1-bdidx,0,Ndof,Ndof);
A = T*Atotal*T + Tbd;

if InitCond == 0
    uppre = zeros(size(mesh.e,1),1);
    upre = zeros(size(mesh.e,1),1);
elseif InitCond ~=0
    tu1 = sum(feval(pde.E1,fem.gex,fem.gey,fem.gez,0).*fem.gew,2);
    tu2 = sum(feval(pde.E2,fem.gex,fem.gey,fem.gez,0).*fem.gew,2);
    tu3 = sum(feval(pde.E3,fem.gex,fem.gey,fem.gez,0).*fem.gew,2);
    tgt = mesh.p(mesh.e(:,2),:) - mesh.p(mesh.e(:,1),:);
    tgt = tgt./sum(tgt.^2,2).^(1/2);
    uppre = tu1.*tgt(:,1) + tu2.*tgt(:,2) + tu3.*tgt(:,3);
    tu1 = sum(feval(pde.Et1,fem.gex,fem.gey,fem.gez,0).*fem.gew,2);
    tu2 = sum(feval(pde.Et2,fem.gex,fem.gey,fem.gez,0).*fem.gew,2);
    tu3 = sum(feval(pde.Et3,fem.gex,fem.gey,fem.gez,0).*fem.gew,2);
    upretmp = tu1.*tgt(:,1) + tu2.*tgt(:,2) + tu3.*tgt(:,3);
    upre = uppre + deltaT*upretmp;
end
UH =zeros(size(mesh.e,1),ceil(Ntime/10)+2);
NumStep = 0;
CurrentT = 0;
SaveCount = 0;
pmt = (mesh.p(mesh.e(:,2),:)+mesh.p(mesh.e(:,1),:))/2;
pf0 = (abs(pmt(:,1)-domain(2))<10^(-10));

while CurrentT < TEND

    NumStep = NumStep + 1;
    CurrentT = CurrentT + deltaT;
    f1Current = @(x,y,z) pde.f1(x,y,z,CurrentT);
    f2Current = @(x,y,z) pde.f2(x,y,z,CurrentT);
    f3Current = @(x,y,z) pde.f3(x,y,z,CurrentT);
    rhsF1 = globNedRHSPGIFE3D(f1Current, mesh, fem, femI, 0, 1);
    rhsF2 = globNedRHSPGIFE3D(f2Current, mesh, fem, femI, 0, 2);
    rhsF3 = globNedRHSPGIFE3D(f3Current, mesh, fem, femI, 0, 3);
    rhsFCurrent = rhsF1 + rhsF2 + rhsF3;

    if BDCond == 0
        tuCurrent = zeros(size(mesh.e,1),1);
    elseif BDCond ~= 0
        exactu1Current = @(x,y,z) pde.E1(x,y,z,CurrentT);
        exactu2Current = @(x,y,z) pde.E2(x,y,z,CurrentT);
        exactu3Current = @(x,y,z) pde.E3(x,y,z,CurrentT);
        tu1 = sum(feval(exactu1Current,fem.gex,fem.gey,fem.gez).*fem.gew,2);
        tu2 = sum(feval(exactu2Current,fem.gex,fem.gey,fem.gez).*fem.gew,2);
        tu3 = sum(feval(exactu3Current,fem.gex,fem.gey,fem.gez).*fem.gew,2);
        tgt = mesh.p(mesh.e(:,2),:) - mesh.p(mesh.e(:,1),:);
        tgt = tgt./sum(tgt.^2,2).^(1/2);
        tuCurrent = tu1.*tgt(:,1) + tu2.*tgt(:,2) + tu3.*tgt(:,3);
        pmt(pf0) = 0;
    end

    ub = tuCurrent;
    ub(fem.mapper) = 0;
    rhsB = Atotal*ub;
    JCurrent = rhsFCurrent - rhsB;
    %JCurrent(isBdEdge) = tuCurrent(isBdEdge);

    fcurrent = Me*(2*upre - uppre)/deltaT^2 + Ms*uppre/(2*deltaT) - S*uppre/2 + JCurrent;
    fcurrent(isBdEdge) = tuCurrent(isBdEdge);

    option.outsolver = 'gmres';
    Eplus = (mesh.eLoc==2);
    Eint = (mesh.eLoc<0);
    alpha = am*ones(size(mesh.e,1),1);
    alpha(Eplus) = ap;
    alpha(Eint) = (am+ap)/2;
    beta = bm*ones(size(mesh.e,1),1);
    beta(Eplus) = bp;
    beta(Eint) = (bm+bp)/2;
    option.alpha = alpha;
    option.beta = beta;
    option.solver = 'amg';
    edge = mesh.e;
    option.smoother = 'BD';
    option.blklevel = 0;
    option.blkId = BasicEdgeNum;
    [x,info] = amgMaxwellinterface2(A,fcurrent,mesh.p,edge,option);
    uppre = upre;
    upre = x;

    if mod(NumStep,10) == 0
       SaveCount = SaveCount+1;
       UH(:,SaveCount) = x;
       disp('save data')
       save('UH5','UH');
    end

end

save('UH5','UH');

%     uh = upre;
%
%     disp(' ');  disp('Start computing error in L2 norm');
%     eNorm = 'L2'; disp(['Start computing error in ',eNorm,'  norm']);
%     [errL2,errL2K1,errL2K2,errL2K3] = getCurlErrIFE3DTime(uh, pde, mesh, fem, femI, eNorm, TEND);
%
%     eNorm = 'Curl'; disp(['Start computing error in ',eNorm,' norm']);
%     [errCurl,errCurl1,errCurl2,errCurl3] = getCurlErrIFE3DTime(uh, pde, mesh, fem, femI, eNorm, TEND);
%     err.l2 = errL2;
%     err.curl = errCurl;
%
%     disp(' ')
%     disp('Errors')
%     disp('L2 norm     H1 norm')
%     formatSpec = '%6.4e  %6.4e\n';
%     fprintf(formatSpec, err.l2, err.curl)
%
%     if i > 1
%         format short
%         rL2 = log(err0.l2/err.l2)./log(h0/h);
%         rcurl = log(err0.curl/err.curl)./log(h0/h);
%         disp(' ')
%         disp('Convergence Rate')
%         disp('L2 norm     Curl norm')
%         formatSpec = '%6.4f      %6.4f\n';
%         fprintf(formatSpec, rL2, rcurl)
%     end
%     err0 = err; h0 = h;

repo: ifem


text:
function ori = calcOrientations(odf,points,varargin)
% draw random orientations from ODF
%
% Syntax
%   ori = calcOrientations(odf,points)
%
% Input
%  odf    - @SO3Fun
%  points - number of orientation to be simualted
%
% Output
%  ori   - @orientation
%
% See also
% ODF_calcPoleFigure, ODF_calcEBSD

warning('The command calcOrientations is depreciated! Please use discreteSample instead.')

ori = discreteSample(odf,points,varargin{:});

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [event] = read_yokogawa_event(filename, varargin)

% READ_YOKOGAWA_EVENT reads event information from continuous,
% epoched or averaged MEG data that has been generated by the Yokogawa
% MEG system and software and allows those events to be used in
% combination with FieldTrip.
%
% Use as
%   [event] = read_yokogawa_event(filename)
%
% See also READ_YOKOGAWA_HEADER, READ_YOKOGAWA_DATA

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

if ~ft_hastoolbox('yokogawa')
    error('cannot determine whether Yokogawa toolbox is present');
end

% get the options
trigindx = keyval('trigindx', varargin); % default is based on chantype helper function

event   = [];
handles = definehandles;

% read the dataset header
hdr = read_yokogawa_header(filename);

if hdr.orig.acq_type==handles.AcqTypeEvokedRaw
  % read the trigger id from all trials
  fid   = fopen(filename, 'r');
  value = GetMeg160TriggerEventM(fid);
  fclose(fid);
  % use the standard FieldTrip header for trial events
  % make an event for each trial as defined in the header
  for i=1:hdr.nTrials
    event(end+1).type     = 'trial';
    event(end  ).sample   = (i-1)*hdr.nSamples + 1;
    event(end  ).offset   = -hdr.nSamplesPre;
    event(end  ).duration =  hdr.nSamples;

    if ~isempty(value)
      event(end  ).value    =  value(i);
    end
  end

elseif hdr.orig.acq_type==handles.AcqTypeEvokedAve
  % make an event for the average
  event(1).type     = 'average';
  event(1).sample   = 1;
  event(1).offset   = -hdr.nSamplesPre;
  event(1).duration =  hdr.nSamples;

elseif hdr.orig.acq_type==handles.AcqTypeContinuousRaw
  % read the trigger channel and detect the flanks
  if isempty(trigindx)
    trigindx = find(hdr.orig.channel_info(:,2)==handles.TriggerChannel);
  end
  event = read_trigger(filename, 'header', hdr, 'chanindx', trigindx, 'detectflank', 'both');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this defines some usefull constants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function handles = definehandles;
handles.output = [];
handles.sqd_load_flag = false;
handles.mri_load_flag = false;
handles.NullChannel         = 0;
handles.MagnetoMeter        = 1;
handles.AxialGradioMeter    = 2;
handles.PlannerGradioMeter  = 3;
handles.RefferenceChannelMark = hex2dec('0100');
handles.RefferenceMagnetoMeter       = bitor( handles.RefferenceChannelMark, handles.MagnetoMeter );
handles.RefferenceAxialGradioMeter   = bitor( handles.RefferenceChannelMark, handles.AxialGradioMeter );
handles.RefferencePlannerGradioMeter = bitor( handles.RefferenceChannelMark, handles.PlannerGradioMeter );
handles.TriggerChannel      = -1;
handles.EegChannel          = -2;
handles.EcgChannel          = -3;
handles.EtcChannel          = -4;
handles.NonMegChannelNameLength = 32;
handles.DefaultMagnetometerSize       = (4.0/1000.0);       % Square of 4.0mm in length
handles.DefaultAxialGradioMeterSize   = (15.5/1000.0);      % Circle of 15.5mm in diameter
handles.DefaultPlannerGradioMeterSize = (12.0/1000.0);      % Square of 12.0mm in length
handles.AcqTypeContinuousRaw = 1;
handles.AcqTypeEvokedAve     = 2;
handles.AcqTypeEvokedRaw     = 3;
handles.sqd = [];
handles.sqd.selected_start  = [];
handles.sqd.selected_end    = [];
handles.sqd.axialgradiometer_ch_no      = [];
handles.sqd.axialgradiometer_ch_info    = [];
handles.sqd.axialgradiometer_data       = [];
handles.sqd.plannergradiometer_ch_no    = [];
handles.sqd.plannergradiometer_ch_info  = [];
handles.sqd.plannergradiometer_data     = [];
handles.sqd.eegchannel_ch_no   = [];
handles.sqd.eegchannel_data    = [];
handles.sqd.nullchannel_ch_no   = [];
handles.sqd.nullchannel_data    = [];
handles.sqd.selected_time       = [];
handles.sqd.sample_rate         = [];
handles.sqd.sample_count        = [];
handles.sqd.pretrigger_length   = [];
handles.sqd.matching_info   = [];
handles.sqd.source_info     = [];
handles.sqd.mri_info        = [];
handles.mri                 = [];

author: PatternRecognition

repo: OpenBMI


text:
function [Estimation_X, FirstPosition] = FEKF_propagate(Estimation_X, OdometryFromThis2Next, Sigma_ODO , FirstPosition )
v = OdometryFromThis2Next(1:3);
w = OdometryFromThis2Next(4:6);

% update position and orientation
Estimation_X.position = Estimation_X.position+Estimation_X.orientation*v;
orientation=Estimation_X.orientation;
Estimation_X.orientation = Estimation_X.orientation*Exp(w);

NumberOfLandmarks = size(Estimation_X.landmarks, 2);
Jrw = J_r(-w);

G1 = [-orientation*Jrw zeros(3,3);zeros(3,3)  -orientation];
G= [G1; zeros(3* NumberOfLandmarks ,6)];
odoCov=diag([w.^2;v.^2])*Sigma_ODO^2;
W = G*odoCov*G';

% compute matrix A_{n}
temp = repmat({ eye(3) }, NumberOfLandmarks+2,1 );
A = blkdiag(temp{:});
A(1:3,1:3) = eye(3);% ExpMinusM;

if isempty(FirstPosition)
   A(4:6,1:3) = -skew(orientation*v);
else
   A(4:6,1:3) = -skew(Estimation_X.position-FirstPosition);
end

% final update the covariance
Estimation_X.cov = A*Estimation_X.cov*A'+W;

FirstPosition=Estimation_X.position;
end

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
function RRest(period)
% RRest runs RR algorithms on ECG and PPG signals using each possible
% combination of options, as specified in "setup_universal_params.m".
%
%               RRest('mimic')
%
%	Inputs:
%		data            data files should be stored in the specified format
%                       in the directory specified in
%                       "setup_universal_params.m". Data can be downloaded
%                       using the "
%       period          this string specifies the dataset to be analysed.
%                       Only the 'mimic' dataset has been used with this
%                       version of the toolbox.
%
%	Outputs:
%       for each subject, N, the following files are made:
%           N_int_respSigs      intermediate respiratory signals,
%           N_respSigs          final respiratory signals
%           N_rrEsts            RR estimates
%           N_rrRef             Reference RR values
%           N_sqi               Signal Quality Index values
%       for the entire dataset, the following files are made:
%           alg_names           Names of RR algorithms tested
%           win_data            Data for every algorithm tested, every
%                               window, and every subject.
%
%   Context:    This is the main file used to run the algorithms. It calls
%               lots of other functions, contained in separate files.
%
%   Further Information:
%       This version of the RRest is provided to facilitate reproduction of
%       the analysis performed in:
%           Charlton P.H. et al., "Waveform Analysis to Estimate
%           Respiratory Rate" [In Press]
%       Further information on this study can be obtained at:
%           http://peterhcharlton.github.io/RRest/waveform_analysis.html
%       In addition, further information on RRest, including future
%       versions, can be obtained at:
%           http://peterhcharlton.github.io/RRest
%
%   Comments, Questions, Criticisms, Feedback, Contributions:
%       See: http://peterhcharlton.github.io/RRest/contributions.html
%
%   Version:
%       v.1 - published on 23rd Feb 2016 by Peter Charlton
%
%   Licence:
%       please see the accompanying file named "LICENSE"
%

%% Setup Universal Parameters
% The universal parameters are used throughout the algorithms
up = setup_universal_params(period);

if up.analysis.run_analysis

    %% Estimate RRs from ECG and PPG
    % Carry out processing for each stage of the algorithms
    for key_comp_no = 1 : length(up.al.key_components)
        feval(up.al.key_components{key_comp_no}, up);
    end

    %% Conduct Signal Quality Assessment of Signals
    % Each window of ECG and PPG is quality assessed
    calculate_sqi(up);

    %% Estimate Reference RRs
    estimate_ref_rr(up);

end

%% Statistical Analysis
calc_stats(up);
create_table_of_algorithms(up);

end

author: peterhcharlton

repo: RRest


text:
function degree = keast_degree ( rule )

%*****************************************************************************80
%
%% KEAST_DEGREE returns the degree of a Keast rule for the tetrahedron.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Patrick Keast,
%    Moderate Degree Tetrahedral Quadrature Formulas,
%    Computer Methods in Applied Mechanics and Engineering,
%    Volume 55, Number 3, May 1986, pages 339-348.
%
%  Parameters:
%
%    Input, integer RULE, the index of the rule.
%
%    Output, integer DEGREE, the degree of the rule.
%
  if ( rule == 1 )
    degree = 0;
  elseif ( rule == 2 )
    degree = 1;
  elseif ( rule == 3 )
    degree = 2;
  elseif ( rule == 4 )
    degree = 3;
  elseif ( rule == 5 )
    degree = 4;
  elseif ( rule == 6 )
    degree = 4;
  elseif ( rule == 7 )
    degree = 5;
  elseif ( rule == 8 )
    degree = 6;
  elseif ( rule == 9 )
    degree = 7;
  elseif ( rule == 10 )
    degree = 8;
  else
    degree = -1;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'KEAST_DEGREE - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of RULE = %d\n', rule );
    error ( 'KEAST_DEGREE - Fatal error!' )
  end

  return
end


text:
function a = triu(a,k)
%TRIU         Implements  triu(a,k)  for gradients
%
%   c = triu(a,k)
%
% functionality as Matlab function triu for matrices
%

% written  10/16/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    improved performance
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if nargin==1
    k = 0;
  end

  a.x = triu(a.x,k);
  index = ( triu( ones(size(a.x)) , k ) == 0 );
  a.dx(index,:) = 0;

repo: OpenMAS


text:
function a = c8mat_identity ( n )

%*****************************************************************************80
%
%% C8MAT_IDENTITY sets the square matrix A to the identity.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, complex A(N,N), the N by N identity matrix.
%
  a(1:n,1:n) = 0.0;

  for i = 1 : n
    a(i,i) = 1.0;
  end

  return
end


text:
function create(obj,wrt)
% Create a NIFTI-1 file
% FORMAT create(obj)
% Write out the header information for the nifti object
%
% FORMAT create(obj,wrt)
% Also write out an empty image volume if wrt==1
%__________________________________________________________________________
% Copyright (C) 2005-2013 Wellcome Trust Centre for Neuroimaging

%
% $Id: create.m 5451 2013-04-26 14:03:05Z guillaume $

for i=1:numel(obj)

    o = obj(i);

    if ~isa(o.dat,'file_array'), error('Data must be a file_array.'); end

    fname = o.dat.fname;
    if isempty(fname), error('No filename to write to.'); end

    %-Write NIFTI header
    sts = write_hdr_raw(fname, o.hdr, o.dat.dtype(end-1)=='B');
    if ~sts, error('Unable to write header for "%s".',fname); end

    %-Write extra information
    write_extras(fname,o.extras);

    %-Create an empty image file if necessary
    if nargin>1 && any(wrt==1)
        [pth,nam] = fileparts(fname);
        if any(strcmp(o.hdr.magic(1:3),{'n+1','n+2'}))
            ext = '.nii';
        else
            ext = '.img';
        end
        o.dat.fname = fullfile(pth,[nam ext]);

        initialise(o.dat);
    end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ x_extend ] = img_extend( x, extension_y, extension_x, extension_z )
% mirrored extension to prevent any edge distortions

% M. Fischer, April 2016

%%
x_extend = zeros(size(x,1)+extension_y,size(x,2)+extension_x,size(x,3)+extension_z);
x_extend(1:end-extension_y,1:end-extension_x,1:end-extension_z) = x;
for j = 1:(extension_z-1)
    x_extend(:,:,end-extension_z+j) = x_extend(:,:,end-extension_z);
end;
x_extend(end-extension_y+1:end,1:end-extension_x,:) = flipud(x(end-extension_y+1:end,:,:));
x_extend(1:end-extension_y,end-extension_x+1:end,:) = fliplr(x(:,end-extension_x+1:end,:));
x_extend(end-extension_y+1:end,end-extension_x+1:end,:) = rot90(x(end-extension_y+1:end,end-extension_x+1:end,:),2);

end

author: thomaskuestner


text:
function isRemoved = db_delete_protocol(isUserConfirm, isRemoveFiles)
% DB_DELETE_PROTOCOL: Remove current protocol from database.
%
% USAGE:  db_delete_protocol(isUserConfirm, isRemoveFiles)  : Remove protocol #iProtocol from protocols list
%
% INPUT:
%     - isUserConfirm : If 0, do not ask user confirmation (default=1)
%     - isRemoveFiles : If 1, delete all the files in this protocol from the hard drive
%                       If 0, keep the files on the hard drive
%                       If empty or not specified: ask user

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

global GlobalData;

%% ===== PARSE INPUTS =====
% Options
if (nargin < 2)
    error('Invalid call to db_delete_protocol().');
end
% Get Protocols list structures (Infos, Subjects, Studies)
sProtocolsListInfo     = GlobalData.DataBase.ProtocolInfo;
sProtocolsListSubjects = GlobalData.DataBase.ProtocolSubjects;
sProtocolsListStudies  = GlobalData.DataBase.ProtocolStudies;
iProtocol              = GlobalData.DataBase.iProtocol;
if isempty(sProtocolsListInfo) || isempty(iProtocol)
    disp('BST> No protocol selected, nothing to delete.');
    return;
end

%% ===== ASK USER CONFIRMATION =====
if isUserConfirm
    % Warning string
    if ~isRemoveFiles
        strWarn = '(Subjects and datasets directories will not be deleted)';
    else
        strWarn = ['
WARNING: All the files will be permanently deleted from your hard drive.
' ...
                   bst_fileparts(sProtocolsListInfo(iProtocol).STUDIES) ''];
    end
    % Display dialog box
    isConfirmed = java_dialog('confirm', ['Remove protocol ''' sProtocolsListInfo(iProtocol).Comment ''' from Brainstorm database ?
' ...
                                          strWarn '

'], sprintf('Remove protocol #%d', iProtocol));
    if ~isConfirmed
        isRemoved = 0;
        return
    end
end

%% ===== REMOVE FILES =====
if isRemoveFiles
    bst_progress('start', 'Delete Protocol', 'Removing all files...');
    % Remove all the contents of STUDIES and SUBJECTS folders
    file_delete( {sProtocolsListInfo(iProtocol).STUDIES, sProtocolsListInfo(iProtocol).SUBJECTS}, 1, 2);
    % If the parent folder (protocol folder) is empty: remove it
    ProtocolDir = bst_fileparts(sProtocolsListInfo(iProtocol).STUDIES);
    listFiles = dir(fullfile(ProtocolDir, '*'));
    listFiles = setdiff({listFiles.name}, {'..','.'});
    if isempty(listFiles)
        rmdir(ProtocolDir);
    else
        warning(['Protocol folder is not empty, cannot be deleted:' ProtocolDir]);
    end
end

%% ===== REMOVE PROTOCOL =====
bst_progress('start', 'Delete Protocol', 'Removing protocol...');
sProtocolsListInfo(iProtocol)     = [];
sProtocolsListSubjects(iProtocol) = [];
sProtocolsListStudies(iProtocol)  = [];
% Update database
GlobalData.DataBase.ProtocolInfo      = sProtocolsListInfo;
GlobalData.DataBase.ProtocolSubjects  = sProtocolsListSubjects;
GlobalData.DataBase.ProtocolStudies   = sProtocolsListStudies;
GlobalData.DataBase.isProtocolLoaded(iProtocol)   = [];
GlobalData.DataBase.isProtocolModified(iProtocol) = [];
% Update protocols ComboBox
gui_brainstorm('UpdateProtocolsList');

% Get new protocol to select
if isempty(sProtocolsListInfo)
    iProtocol = 0;
elseif (iProtocol == 1)
    iProtocol = 1;
else
    iProtocol = iProtocol - 1;
end
% Select current protocol in combo list
gui_brainstorm('SetCurrentProtocol', iProtocol);

% Save database
% db_save();
isRemoved = 1;

bst_progress('stop');

author: brainstorm-tools


text:
function calpak_test016 ( )

%*****************************************************************************80
%
%% CALPAK_TEST016 tests JED_TO_YMDF_COPTIC and YMDF_TO_JED_COPTIC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST016\n' );
  fprintf ( 1, '  For the Coptic calendar:\n' );
  fprintf ( 1, '  JED_TO_YMDF_COPTIC: JED -> YMDF.\n' );
  fprintf ( 1, '  YMDF_TO_JED_COPTIC: YMDF -> JED.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  JED (in)    YMDF               JED (out)\n' );
  fprintf ( 1, '\n' );

  jed_epoch = epoch_to_jed_coptic ( );

  i = 0;

  while ( 1 )

    i = i + 1;
    jed1 = jed_test ( i );

    if ( jed1 < 0.0 )
      break
    end

    if ( jed_epoch <= jed1 )

      [ y2, m2, d2, f2 ] = jed_to_ymdf_coptic ( jed1 );

      s2 = ymdf_to_s_numeric ( y2, m2, d2, f2 );

      jed3 = ymdf_to_jed_coptic ( y2, m2, d2, f2 );

      fprintf ( 1, '  %11.2f  %20s  %11.2f\n', jed1, s2, jed3 );

    end

  end

  return
end


text:
%Terminal organelle assembly test cases
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 7/13/2010
classdef TerminalOrganelleAssembly_Test < edu.stanford.covert.cell.sim.ReactionProcessTestCase
    methods
        function this = TerminalOrganelleAssembly_Test(name)
            this = this@edu.stanford.covert.cell.sim.ReactionProcessTestCase(name);
        end

        function testSimpleFixture(this)
            this.loadSimpleTestFixture();
            m = this.process;

            m.evolveState();
            assertAllEqual(0, m.substrates(:, m.compartmentIndexs_unincorporated));
            assertAllEqual(20, m.substrates(:, m.compartmentIndexs_incorporated));
        end

        function testNoIncorporatedProteins(this)
            m = this.process;
            m.substrates(:, m.compartmentIndexs_incorporated) = 0;
            m.substrates(:, m.compartmentIndexs_unincorporated) = 10;

            m.evolveState();
            assertAllEqual(0, m.substrates(:, m.compartmentIndexs_unincorporated));
            assertAllEqual(10, m.substrates(:, m.compartmentIndexs_incorporated));
        end

        function testNoUnincorporatedProteins(this)
            m = this.process;
            m.substrates(:, m.compartmentIndexs_incorporated) = 10;
            m.substrates(:, m.compartmentIndexs_unincorporated) = 0;

            m.evolveState();
            assertAllEqual(0, m.substrates(:, m.compartmentIndexs_unincorporated));
            assertAllEqual(10, m.substrates(:, m.compartmentIndexs_incorporated));
        end

        function testNoHMW12(this)
            m = this.process;

            m.substrates(:,m.compartmentIndexs_unincorporated) = 1e6;
            m.substrates(:,m.compartmentIndexs_incorporated) = 0;
            m.substrates(m.substrateIndexs_HMW1,:) = 0;
            m.substrates(m.substrateIndexs_HMW2,:) = 0;
            initialState = struct('substrates', m.substrates);
            m.evolveState();
            assertFalse(someUseOfAllSubstrates(m, initialState));

            m.substrates(:) = 100;
            m.substrates(m.substrateIndexs_HMW1, m.compartmentIndexs_incorporated) = 0;
            m.substrates(m.substrateIndexs_HMW2, m.compartmentIndexs_incorporated) = 0;
            initialState = struct('substrates', m.substrates);
            m.evolveState();
            assertTrue(someUseOfAllSubstrates(m, initialState));

            m.substrates(:) = 100;
            m.substrates(m.substrateIndexs_HMW1, m.compartmentIndexs_unincorporated) = 0;
            m.substrates(m.substrateIndexs_HMW2, m.compartmentIndexs_unincorporated) = 0;
            initialState = struct('substrates', m.substrates);
            m.evolveState();
            assertTrue(someUseOfAllSubstrates(m, initialState));
        end

        function testNoHMW1(this)
            m = this.process;

            m.substrates(:) = 100;
            m.substrates(m.substrateIndexs_HMW1,:)=0;
            initialState = struct('substrates', m.substrates);
            m.evolveState();
            assertFalse(someUseOfAllSubstrates(m, initialState));

            m.substrates(:) = 100;
            m.substrates(m.substrateIndexs_HMW1, m.compartmentIndexs_incorporated) = 0;
            initialState = struct('substrates', m.substrates);
            m.evolveState();
            assertTrue(someUseOfAllSubstrates(m, initialState));

            m.substrates(:) = 100;
            m.substrates(m.substrateIndexs_HMW1, m.compartmentIndexs_unincorporated) = 0;
            initialState = struct('substrates', m.substrates);
            m.evolveState();
            assertTrue(someUseOfAllSubstrates(m, initialState));
        end

        function testNoHMW2(this)
            m = this.process;

            m.substrates(:) = 100;
            m.substrates(m.substrateIndexs_HMW2,:) = 0;
            initialState = struct('substrates', m.substrates);
            m.evolveState();
            assertFalse(someUseOfAllSubstrates(m, initialState));

            m.substrates(:) = 100;
            m.substrates(m.substrateIndexs_HMW2, m.compartmentIndexs_incorporated) = 0;
            initialState = struct('substrates', m.substrates);
            m.evolveState();
            assertTrue(someUseOfAllSubstrates(m, initialState));

            m.substrates(:) = 100;
            m.substrates(m.substrateIndexs_HMW2, m.compartmentIndexs_unincorporated) = 0;
            initialState = struct('substrates', m.substrates);
            m.evolveState();
            assertTrue(someUseOfAllSubstrates(m, initialState));
        end

        function testGeneEssentiality(this)
            m = this.process;
            m.substrates(:) = 100;
            this.helpTestGeneEssentiality({
                'MG_218';     %hmw2, high molecular weight cytadherence accessory protein 2
                'MG_312';     %hmw1, high molecular weight cytadherence accessory protein 1
                'MG_317';     %hmw3, high molecular weight cytadherence accessory protein 3
                'MG_318'},... %p32, P32 adhesin
                @someUseOfAllSubstrates);
        end
    end

    methods (Access = private)
        function loadSimpleTestFixture(this)
            m = this.process;

            enzymeSubstrateIndexs = [
                m.enzymeIndexs_HMW1;
                m.enzymeIndexs_HMW2;
                m.enzymeIndexs_HMW3;
                m.enzymeIndexs_P32];

            m.substrateWholeCellModelIDs = sort(m.substrateWholeCellModelIDs__);
            m.enzymeWholeCellModelIDs = m.substrateWholeCellModelIDs(enzymeSubstrateIndexs);

            m.reactionWholeCellModelIDs = {
                'MG_191_TerminalOrganelleAssembly_01';
                'MG_192_TerminalOrganelleAssembly_01';
                'MG_217_TerminalOrganelleAssembly_01';
                'MG_218_TerminalOrganelleAssembly_01';
                'MG_218_TerminalOrganelleAssembly_02';
                'MG_312_TerminalOrganelleAssembly_01';
                'MG_312_TerminalOrganelleAssembly_02';
                'MG_317_TerminalOrganelleAssembly_01';
                'MG_318_TerminalOrganelleAssembly_01';
                'MG_386_TerminalOrganelleAssembly_01'};
            m.reactionNames = m.reactionWholeCellModelIDs;

            nSubstrates = length(m.substrateWholeCellModelIDs);
            nReactions  = length(m.reactionWholeCellModelIDs);

            m.localizationSubstrates = zeros(nReactions, nSubstrates);
            m.localizationSubstrates(1, substrateIdx('MG_191_MONOMER')) = 1;
            m.localizationSubstrates(2, substrateIdx('MG_192_MONOMER')) = 1;
            m.localizationSubstrates(3, substrateIdx('MG_217_MONOMER')) = 1;
            m.localizationSubstrates(4, substrateIdx('MG_218_MONOMER')) = 1;
            m.localizationSubstrates(5, substrateIdx('MG_218_MONOMER')) = 1;
            m.localizationSubstrates(6, substrateIdx('MG_312_MONOMER')) = 1;
            m.localizationSubstrates(7, substrateIdx('MG_312_MONOMER')) = 1;
            m.localizationSubstrates(8, substrateIdx('MG_317_MONOMER')) = 1;
            m.localizationSubstrates(9, substrateIdx('MG_318_MONOMER')) = 1;
            m.localizationSubstrates(10, substrateIdx('MG_386_MONOMER')) = 1;

            i = m.compartmentIndexs_incorporated;
            u = m.compartmentIndexs_unincorporated;

            m.localizationReactions = zeros(nReactions, nSubstrates, 2);
            m.localizationReactions(1, substrateIdx('MG_312_MONOMER'), i) = 1;
            m.localizationReactions(3, substrateIdx('MG_318_MONOMER'), i) = 1;
            m.localizationReactions(4, substrateIdx('MG_312_MONOMER'), u) = 1;
            m.localizationReactions(5, substrateIdx('MG_312_MONOMER'), i) = 1;
            m.localizationReactions(6, substrateIdx('MG_218_MONOMER'), u) = 1;
            m.localizationReactions(7, substrateIdx('MG_218_MONOMER'), i) = 1;
            m.localizationReactions(8, substrateIdx('MG_312_MONOMER'), i) = 1;
            m.localizationReactions(9, substrateIdx('MG_317_MONOMER'), i) = 1;
            m.localizationReactions(10, substrateIdx('MG_312_MONOMER'), i) = 1;

            m.localizationThreshold = sum(sum(m.localizationReactions,3),2);

            m.substrates = repmat(10, length(m.substrateWholeCellModelIDs), 2);
            m.enzymes    = m.substrates(enzymeSubstrateIndexs, :);

            function i = substrateIdx(id)
                [~,i] = ismember(id, m.substrateWholeCellModelIDs);
            end
        end
    end
end

function result = someUseOfAllSubstrates(m, i)
    result = all(...
        i.substrates(:,m.compartmentIndexs_unincorporated)==0 | ...
        i.substrates(:,m.compartmentIndexs_unincorporated) > ...
        m.substrates(:,m.compartmentIndexs_unincorporated));
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
% simple demo to convert a polygon mesh into a voxel representation

% code is 99.99% based on
% [1] http://www.mathworks.com/matlabcentral/fileexchange/24086-polygon2voxel
% [2] http://www.mathworks.com/matlabcentral/fileexchange/21044-3d-voxelizer

%   % Compile the c-coded function
%   mex polygon2voxel_double.c -v

load model;

vertices = vertices - repmat(mean(vertices,1),size(vertices,1),1);

FV.faces = faces;
FV.vertices = vertices;

Volume=polygon2voxel(FV,[20 20 20],'auto');

%% visualization 1
figure
[X,Y,Z]=ind2sub(size(Volume),find(Volume(:)));
plot3(X,Y,Z,'.');
axis equal;
xlabel('x');
ylabel('y');
zlabel('z');

%% visualization 2
% 3d pirnter style visualization to add layer by layer

careMask =  imdilate((Volume),ones(2,2,2));
figure,plot3D(careMask,1,'timed', 0.1)
hold on;%plot3D(Volume,1,'b','*')
 saveas(gcf,sprintf('/Users/shurans/Dropbox/ModelNet/ECCV2014_sub/MaskVideo/%05d.png',21 ));
viewA =33:5:176;
for i =1:length(viewA)
    view(viewA(i),28);
    saveas(gcf,sprintf('/Users/shurans/Dropbox/ModelNet/ECCV2014_sub/MaskVideo/%05d.png',21+i ));
end
%%
%{
%% visualization 3
figure
for i=1:size(Volume,1)
    imagesc(squeeze(Volume(i,:,:)));
    axis equal;
    axis tight;
    axis off
    title(i);
    pause(0.1);
end

for i=1:size(Volume,2)
    imagesc(squeeze(Volume(:,i,:)));
    axis equal;
    axis tight;
    axis off
    title(i);
    pause(0.1);
end

for i=1:size(Volume,3)
    imagesc(squeeze(Volume(:,:,i)));
    axis equal;
    axis tight;
    axis off
    title(i);
    pause(0.1);
end
%}

author: zhirongw


text:
function setT(OCP,T)
    % assert
%     if isa(T,'numeric')
%         if T<0
%             error('T must be a positive number or parameter!');
%         end
%     else
%         if has(T,OCP.u) || has(T,OCP.x)
%             error('T cannot be an input or state!');
%         end
%     end

    OCP.T = T;
    OCP.deltaTau = OCP.T/OCP.N;
end

author: deng-haoyang

repo: ParNMPC


text:
function f_x = ParFor2(in1)
%PARFOR2
%    F_X = PARFOR2(IN1)

%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    20-Sep-2019 09:35:34

u = in1(:,1);
ux = in1(:,4);
uxx = in1(:,5);
uxxx = in1(:,6);
f_x = (u.*1.034234549167718e16+ux.*6.60584588804096e15+uxx.*1.215161942676603e17-uxxx.*4.609471812108288e15-u.*ux.*5.122148099342336e15+u.*uxx.*5.033451929390285e16+u.*uxxx.*5.000975400599552e15-2.603066969792512e15)./(u.*3.126164642771763e16+4.756339916996608e16);

author: dynamicslab

repo: SINDy-PI