text:
function x = fnlCg(x0,params)
%-----------------------------------------------------------------------
%
% res = fnlCg(x0,params)
%
% implementation of a L1 penalized non linear conjugate gradient reconstruction
%
% The function solves the following problem:
%
% given k-space measurments y, and a fourier operator F the function
% finds the image x that minimizes:
%
% Phi(x) = ||F* W' *x - y||^2 + lambda1*|x|_1 + lambda2*TV(W'*x)
%
%
% the optimization method used is non linear conjugate gradient with fast&cheap backtracking
% line-search.
%
% (c) Michael Lustig 2007
%-------------------------------------------------------------------------
x = x0;

% line search parameters
maxlsiter = params.lineSearchItnlim ;
gradToll = params.gradToll ;
alpha = params.lineSearchAlpha; ,    beta = params.lineSearchBeta;
t0 = params.lineSearchT0;
k = 0;
t = 1;

% copmute g0  = grad(Phi(x))

g0 = wGradient(x,params);

dx = -g0;

% iterations
while(1)

% backtracking line-search

	% pre-calculate values, such that it would be cheap to compute the objective
	% many times for efficient line-search
	[FTXFMtx, FTXFMtdx, DXFMtx, DXFMtdx] = preobjective(x, dx, params);
	f0 = objective(FTXFMtx, FTXFMtdx, DXFMtx, DXFMtdx,x,dx, 0, params);
	t = t0;
        [f1, ERRobj, RMSerr]  =  objective(FTXFMtx, FTXFMtdx, DXFMtx, DXFMtdx,x,dx, t, params);

	lsiter = 0;

	while (f1 > f0 - alpha*t*abs(g0(:)'*dx(:)))^2 & (lsiter<maxlsiter)
		lsiter = lsiter + 1;
		t = t * beta;
		[f1, ERRobj, RMSerr]  =  objective(FTXFMtx, FTXFMtdx, DXFMtx, DXFMtdx,x,dx, t, params);
	end

	if lsiter == maxlsiter
		disp('Reached max line search,.... not so good... might have a bug in operators. exiting... ');
		return;
	end

	% control the number of line searches by adapting the initial step search
	if lsiter > 2
		t0 = t0 * beta;
	end

	if lsiter<1
		t0 = t0 / beta;
	end

	x = (x + t*dx);

	%--------- uncomment for debug purposes ------------------------
	disp(sprintf('%d   , obj: %f, RMS: %f, L-S: %d', k,f1,RMSerr,lsiter));

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

    %conjugate gradient calculation

	g1 = wGradient(x,params);
	bk = g1(:)'*g1(:)/(g0(:)'*g0(:)+eps);
	g0 = g1;
	dx =  - g1 + bk* dx;
	k = k + 1;

	%TODO: need to "think" of a "better" stopping criteria ;-)
	if (k > params.Itnlim) | (norm(dx(:)) < gradToll)
		break;
	end

end

return;

function [FTXFMtx, FTXFMtdx, DXFMtx, DXFMtdx] = preobjective(x, dx, params)

% precalculates transforms to make line search cheap

FTXFMtx = params.FT*(params.XFM'*x);
FTXFMtdx = params.FT*(params.XFM'*dx);

if params.TVWeight
    DXFMtx = params.TV*(params.XFM'*x);
    DXFMtdx = params.TV*(params.XFM'*dx);
else
    DXFMtx = 0;
    DXFMtdx = 0;
end

function [res, obj, RMS] = objective(FTXFMtx, FTXFMtdx, DXFMtx, DXFMtdx, x,dx,t, params);
%calculated the objective function

p = params.pNorm;

obj = FTXFMtx + t*FTXFMtdx - params.data;
obj = obj(:)'*obj(:);

if params.TVWeight
    w = DXFMtx(:) + t*DXFMtdx(:);
    TV = (w.*conj(w)+params.l1Smooth).^(p/2);
else
    TV = 0;
end

if params.xfmWeight
   w = x(:) + t*dx(:);
   XFM = (w.*conj(w)+params.l1Smooth).^(p/2);
else
    XFM=0;
end

TV = sum(TV.*params.TVWeight(:));
XFM = sum(XFM.*params.xfmWeight(:));
RMS = sqrt(obj/sum(abs(params.data(:))>0));

res = obj + (TV) + (XFM) ;

function grad = wGradient(x,params)

gradXFM = 0;
gradTV = 0;

gradObj = gOBJ(x,params);
if params.xfmWeight
gradXFM = gXFM(x,params);
end
if params.TVWeight
gradTV = gTV(x,params);
end

grad = (gradObj +  params.xfmWeight.*gradXFM + params.TVWeight.*gradTV);

function gradObj = gOBJ(x,params)
% computes the gradient of the data consistency

	gradObj = params.XFM*(params.FT'*(params.FT*(params.XFM'*x) - params.data));

gradObj = 2*gradObj ;

function grad = gXFM(x,params)
% compute gradient of the L1 transform operator

p = params.pNorm;

grad = p*x.*(x.*conj(x)+params.l1Smooth).^(p/2-1);

function grad = gTV(x,params)
% compute gradient of TV operator

p = params.pNorm;

Dx = params.TV*(params.XFM'*x);

G = p*Dx.*(Dx.*conj(Dx) + params.l1Smooth).^(p/2-1);
grad = params.XFM*(params.TV'*G);

author: thomaskuestner


text:
function [ center, radii, evecs, v, chi2 ] = ellipsoid_fit_new( X, equals )
%
% Fit an ellispoid/sphere/paraboloid/hyperboloid to a set of xyz data points:
%
%   [center, radii, evecs, pars ] = ellipsoid_fit( X )
%   [center, radii, evecs, pars ] = ellipsoid_fit( [x y z] );
%   [center, radii, evecs, pars ] = ellipsoid_fit( X, 1 );
%   [center, radii, evecs, pars ] = ellipsoid_fit( X, 2, 'xz' );
%   [center, radii, evecs, pars ] = ellipsoid_fit( X, 3 );
%
% Parameters:
% * X, [x y z]   - Cartesian data, n x 3 matrix or three n x 1 vectors
% * flag         - '' or empty fits an arbitrary ellipsoid (default),
%                - 'xy' fits a spheroid with x- and y- radii equal
%                - 'xz' fits a spheroid with x- and z- radii equal
%                - 'xyz' fits a sphere
%                - '0' fits an ellipsoid with its axes aligned along [x y z] axes
%                - '0xy' the same with x- and y- radii equal
%                - '0xz' the same with x- and z- radii equal
%
% Output:
% * center    -  ellispoid or other conic center coordinates [xc; yc; zc]
% * radii     -  ellipsoid or other conic radii [a; b; c]
% * evecs     -  the radii directions as columns of the 3x3 matrix
% * v         -  the 10 parameters describing the ellipsoid / conic algebraically:
%                Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + 2Gx + 2Hy + 2Iz + J = 0
% * chi2      -  residual sum of squared errors (chi^2), this chi2 is in the
%                coordinate frame in which the ellipsoid is a unit sphere.
%
% Author:
% Yury Petrov, Oculus VR
% Date:
% September, 2015
%

%reference
narginchk( 1, 3 ) ;  % check input arguments
if nargin == 1
    equals = ''; % no constraints by default
end

if size( X, 2 ) ~= 3
    error( 'Input data must have three columns!' );
else
    x = X( :, 1 );
    y = X( :, 2 );
    z = X( :, 3 );
end

% need nine or more data points
if length( x ) < 9 && strcmp( equals, '' )
   error( 'Must have at least 9 points to fit a unique ellipsoid' );
end
if length( x ) < 8 && ( strcmp( equals, 'xy' ) || strcmp( equals, 'xz' ) )
   error( 'Must have at least 8 points to fit a unique ellipsoid with two equal radii' );
end
if length( x ) < 6 && strcmp( equals, '0' )
   error( 'Must have at least 6 points to fit a unique oriented ellipsoid' );
end
if length( x ) < 5 && ( strcmp( equals, '0xy' ) || strcmp( equals, '0xz' ) )
   error( 'Must have at least 5 points to fit a unique oriented ellipsoid with two equal radii' );
end
if length( x ) < 4 && strcmp( equals, 'xyz' );
   error( 'Must have at least 4 points to fit a unique sphere' );
end

% fit ellipsoid in the form Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + 2Gx +
% 2Hy + 2Iz + J = 0 and A + B + C = 3 constraint removing one extra
% parameter
if strcmp( equals, '' )
    D = [ x .* x + y .* y - 2 * z .* z, ...
        x .* x + z .* z - 2 * y .* y, ...
        2 * x .* y, ...
        2 * x .* z, ...
        2 * y .* z, ...
        2 * x, ...
        2 * y, ...
        2 * z, ...
        1 + 0 * x ];  % ndatapoints x 9 ellipsoid parameters
elseif strcmp( equals, 'xy' )
    D = [ x .* x + y .* y - 2 * z .* z, ...
        2 * x .* y, ...
        2 * x .* z, ...
        2 * y .* z, ...
        2 * x, ...
        2 * y, ...
        2 * z, ...
        1 + 0 * x ];  % ndatapoints x 8 ellipsoid parameters
elseif strcmp( equals, 'xz' )
    D = [ x .* x + z .* z - 2 * y .* y, ...
        2 * x .* y, ...
        2 * x .* z, ...
        2 * y .* z, ...
        2 * x, ...
        2 * y, ...
        2 * z, ...
        1 + 0 * x ];  % ndatapoints x 8 ellipsoid parameters
    % fit ellipsoid in the form Ax^2 + By^2 + Cz^2 + 2Gx + 2Hy + 2Iz = 1
elseif strcmp( equals, '0' )
    D = [ x .* x + y .* y - 2 * z .* z, ...
          x .* x + z .* z - 2 * y .* y, ...
          2 * x, ...
          2 * y, ...
          2 * z, ...
          1 + 0 * x ];  % ndatapoints x 6 ellipsoid parameters
    % fit ellipsoid in the form Ax^2 + By^2 + Cz^2 + 2Gx + 2Hy + 2Iz = 1,
    % where A = B or B = C or A = C
elseif strcmp( equals, '0xy' )
    D = [ x .* x + y .* y - 2 * z .* z, ...
          2 * x, ...
          2 * y, ...
          2 * z, ...
          1 + 0 * x ];  % ndatapoints x 5 ellipsoid parameters
elseif strcmp( equals, '0xz' )
    D = [ x .* x + z .* z - 2 * y .* y, ...
          2 * x, ...
          2 * y, ...
          2 * z, ...
          1 + 0 * x ];  % ndatapoints x 5 ellipsoid parameters
     % fit sphere in the form A(x^2 + y^2 + z^2) + 2Gx + 2Hy + 2Iz = 1
elseif strcmp( equals, 'xyz' )
    D = [ 2 * x, ...
          2 * y, ...
          2 * z, ...
          1 + 0 * x ];  % ndatapoints x 4 ellipsoid parameters
else
    error( [ 'Unknown parameter value ' equals '!' ] );
end

% solve the normal system of equations
d2 = x .* x + y .* y + z .* z; % the RHS of the llsq problem (y's)
u = ( D' * D ) \ ( D' * d2 );  % solution to the normal equations

% find the residual sum of errors
% chi2 = sum( ( 1 - ( D * u ) ./ d2 ).^2 ); % this chi2 is in the coordinate frame in which the ellipsoid is a unit sphere.

% find the ellipsoid parameters
% convert back to the conventional algebraic form
if strcmp( equals, '' )
    v(1) = u(1) +     u(2) - 1;
    v(2) = u(1) - 2 * u(2) - 1;
    v(3) = u(2) - 2 * u(1) - 1;
    v( 4 : 10 ) = u( 3 : 9 );
elseif strcmp( equals, 'xy' )
    v(1) = u(1) - 1;
    v(2) = u(1) - 1;
    v(3) = -2 * u(1) - 1;
    v( 4 : 10 ) = u( 2 : 8 );
elseif strcmp( equals, 'xz' )
    v(1) = u(1) - 1;
    v(2) = -2 * u(1) - 1;
    v(3) = u(1) - 1;
    v( 4 : 10 ) = u( 2 : 8 );
elseif strcmp( equals, '0' )
    v(1) = u(1) +     u(2) - 1;
    v(2) = u(1) - 2 * u(2) - 1;
    v(3) = u(2) - 2 * u(1) - 1;
    v = [ v(1) v(2) v(3) 0 0 0 u( 3 : 6 )' ];

elseif strcmp( equals, '0xy' )
    v(1) = u(1) - 1;
    v(2) = u(1) - 1;
    v(3) = -2 * u(1) - 1;
    v = [ v(1) v(2) v(3) 0 0 0 u( 2 : 5 )' ];
elseif strcmp( equals, '0xz' )
    v(1) = u(1) - 1;
    v(2) = -2 * u(1) - 1;
    v(3) = u(1) - 1;
    v = [ v(1) v(2) v(3) 0 0 0 u( 2 : 5 )' ];
elseif strcmp( equals, 'xyz' )
    v = [ -1 -1 -1 0 0 0 u( 1 : 4 )' ];
end
v = v';

% form the algebraic form of the ellipsoid
A = [ v(1) v(4) v(5) v(7); ...
      v(4) v(2) v(6) v(8); ...
      v(5) v(6) v(3) v(9); ...
      v(7) v(8) v(9) v(10) ];
% find the center of the ellipsoid
center = -A( 1:3, 1:3 ) \ v( 7:9 );
% form the corresponding translation matrix
T = eye( 4 );
T( 4, 1:3 ) = center';
% translate to the center
R = T * A * T';
% solve the eigenproblem
[ evecs, evals ] = eig( R( 1:3, 1:3 ) / -R( 4, 4 ) );
radii = sqrt( 1 ./ diag( abs( evals ) ) );
sgns = sign( diag( evals ) );
radii = radii .* sgns;

% calculate difference of the fitted points from the actual data normalized by the conic radii
d = [ x - center(1), y - center(2), z - center(3) ]; % shift data to origin
d = d * evecs; % rotate to cardinal axes of the conic;
d = [ d(:,1) / radii(1), d(:,2) / radii(2), d(:,3) / radii(3) ]; % normalize to the conic radii
chi2 = sum( abs( 1+0*x - sum( d.^2, 2 ) ) );

if abs( v(end) ) > 1e-6
    v = -v / v(end); % normalize to the more conventional form with constant term = -1
else
    v = -sign( v(end) ) * v;
end

author: yandld


text:
function pstruct = tapas_hgf_namep(pvec)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

pstruct = struct;

l = length(pvec)/5;

if l ~= floor(l)
    error('tapas:hgf:UndetNumLevels', 'Cannot determine number of levels');
end

pstruct.mu_0      = pvec(1:l);
pstruct.sa_0      = pvec(l+1:2*l);
pstruct.rho       = pvec(2*l+1:3*l);
pstruct.ka        = pvec(3*l+1:4*l-1);
pstruct.om        = pvec(4*l:5*l-1);
pstruct.pi_u      = pvec(5*l);

return;

author: translationalneuromodeling

repo: tapas


text:
function cmap = grayColorCmap(numGrays,numColors)
%
% cmap = grayColorCmap(numGrays,numColors)
%
% Makes colormap array with:
%   gray scale - 1:numGrays
%   jet colors - numGrays+1:numGrays+numColors
%
% This one differs from grayCmap.m as this one has both gray band and color
% band in gray color.

if ~exist('numGrays','var')
    numGrays = 128;
end
if ~exist('numColors','var')
    numColors = 128;
end

cmap = zeros(numGrays+numColors,3);
cmap(1:numGrays+numColors,:) = [gray(numGrays);gray(numColors)];

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function p25_p_set ( p )

%*****************************************************************************80
%
%% P25_P_SET sets the value of a parameter for problem 25.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P, the value of the parameter that multiplies the
%    Gaussian noise.
%
  p25_p_val ( 'SET', p );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% A script for visualizing the orientations in the menpo data, allows to
% see the distributions of the data
clear
load('../menpo_68_pts.mat');
addpath('../../PDM_helpers');

xs = all_pts(1:end/2,:);
ys = all_pts(end/2+1:end,:);
num_imgs = size(xs, 1);

rots = zeros(3, num_imgs);
errs = zeros(1,num_imgs);

pdmLoc = ['../pdm_68_aligned_menpo.mat'];

load(pdmLoc);

pdm = struct;
pdm.M = double(M);
pdm.E = double(E);
pdm.V = double(V);
errs_poss = [];
for i=1:num_imgs

    labels_curr = cat(2, xs(i,:)', ys(i,:)');
    labels_curr(labels_curr==-1) = 0;

    [ a, R, T, ~, l_params, err, shapeOrtho] = fit_PDM_ortho_proj_to_2D(pdm.M, pdm.E, pdm.V, labels_curr);
    errs(i) = err/a;
    rots(:,i) = Rot2Euler(R);

end

hist(rots', 100);

author: TadasBaltrusaitis

repo: OpenFace


text:
function errors = armaxerrors(parameters,p,q,constant,y,x,m,sigma)
% PURPOSE:
%   Compute errors from an ARMAX model for use in a least squares optimizer
%
% USAGE:
%   [ERRORS] = armaxfilter_likelihood(PARAMETERS,P,Q,CONSTANT,Y,X,M,SIGMA)
%
% INPUTS:
%   PARAMETERS - A vector of GARCH process aprams of the form [constant, arch, garch]
%   P          - Vector containing lag indices of the AR component
%   Q          - Vector containing lag indices of the MA component
%   CONSTANT   - Value indicating whether the model contains a constant (1) or not (0)
%   Y          - Data augments with max(max(Q)-max(P),0) zeros
%   X          - Regressors augmented with max(max(Q)-max(P),0) zeros
%   M          - Index to first element to use in the recursive residual calculation
%   SIGMA      - Vector of conditional standard deviations with the same dimension as Y for use in
%                  GLS estimation
%
% OUTPUTS:
%   ERRORS     - Vector of errors with the same size as Y.  First M elements are 0.
%
% COMMENTS:
%
%  See also ARMAXFILTER_LIKELIHOOD

% Author: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 10/19/2009

np = length(p);
nq = length(q);
T = length(y);
errors = zeros(T,1);
k = size(x,2);

for t=m+1:T
    errors(t) = y(t);
    if constant
        errors(t) = errors(t) - parameters(1);
    end
    for i=1:np
        errors(t) = errors(t) - parameters(constant+i)*y(t-p(i));
    end
    for i=1:k
        errors(t) = errors(t) - parameters(constant+np+i)*x(t,i);
    end
    for i=1:nq
        errors(t) = errors(t) - parameters(constant+np+k+i)*errors(t-q(i));
    end
    errors(t) = errors(t);
end
errors = errors./sigma;

author: bashtage

repo: mfe-toolbox


text:
function value = ch_is_alpha ( c )

%*****************************************************************************80
%
%% CH_IS_ALPHA returns TRUE if C is an alphabetic character.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, a character to check.
%
%    Output, logical VALUE is TRUE if C is an alphabetic character.
%
  if ( ( 'a' <= c && c <= 'z' ) || ...
       ( 'A' <= c && c <= 'Z' ) )
    value = 1;
  else
    value = 0;
  end

  return
end


text:
function d = calc(e,d)

if 0
   if isa(d,'cell')
     for i=1:length(d)
       d{i}=feval(e.loss_type,e,d{i});
     end
   else
     d=feval(e.loss_type,e,d);
   end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [gv gh]=sobel_fn(x)

%%% sobel_fn: Computes the vertical & horizontal edges of x using sobel operator.
%%%
%%% [gv gh]=sobel_fn(x)
%%%
%%% Author : B. K. SHREYAMSHA KUMAR
%%% Created on 28-10-2011.
%%% Updated on 28-10-2011.

vtemp=[-1 0 1;-2 0 2;-1 0 1]/8;
htemp=[-1 -2 -1;-0 0 0;1 2 1]/8;

[a b]=size(htemp);
x_ext=per_extn_im_fn(x,a);
[p,q]=size(x_ext);
for ii=2:p-1
   for jj=2:q-1
      gv(ii-1,jj-1)=sum(sum(x_ext(ii-1:ii+1,jj-1:jj+1).*vtemp));
      gh(ii-1,jj-1)=sum(sum(x_ext(ii-1:ii+1,jj-1:jj+1).*htemp));
   end
end

author: Linfeng-Tang

repo: Image-Fusion


text:
function calpak_test51535 ( )

%*****************************************************************************80
%
%% CALPAK_TEST51535 tests YEAR_IS_EMBOLISMIC_EG_LUNAR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST51535\n' );
  fprintf ( 1, '  For the Egyptian Lunar calendar:\n' );
  fprintf ( 1, '  YEAR_IS_EMBOLISMIC_EG_LUNAR determines if a year is\n' );
  fprintf ( 1, '  an embolismic year.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Year  Embolismic?\n' );
  fprintf ( 1, '\n' );

  for y = 1 : 25
    s = y_to_s_eg_lunar ( y );
    fprintf ( 1, '  %s  %d\n', s, year_is_embolismic_eg_lunar ( y ) );
  end

  return
end


text:
function [ h, hmin, hmax ] = p04_stepsize ( option )

%*****************************************************************************80
%
%% P04_STEPSIZE returns step sizes for problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, real H, HMIN, HMAX, suggested values for the
%    initial step, the minimum step, and the maximum step.
%
  h =     0.300;
  hmin =  0.001;
  hmax = 25.000;

  return
end


text:
function [rc] = calc_ratechangeF(mycat,time,timef,bootloops,maepi)
    % function [rc] = calc_ratechangeF(mycat,time,timef,bootloops,maepi);
    % ----------------------------------------------------------------
    % Determines ratechanges within aftershock sequences for defined time window
    %
    % Input parameters:
    %   mycat       earthquake catalog
    %   time_as     delay times (days)
    %   step        number of quakes to determine forecast period
    %   time        learning period
    %   timeF       forecast period
    %   bootloops   Number of bootstraps
    %   ZG.maepi       Mainsock values
    %
    % Output parameters:
    %   rc      Matrix containing: time, absdiff, sigma, numreal, nummod,
    %           pval, pvalstd, cval, cvalstd, kval, kvalstd, t_forecast, fMc
    %
    % Info: Version for ZMAP, original function by ratechange.m by S. Neukomm May 27, 2002
    % J. Woessner
    % last update: 09.07.03

    % if min(time_as) >= 1
    %     rc = [];
    %     return
    % end

    % Initialize
    rc = [];
    % fMc0 = mcwithtime(mycat,time_as,7);
    %
    % %time = mintime;
    %
    % % determination of magnitude of completeness
    % if time < 7
    %     fMc = mcwithtime(mycat,time_as,time);
    % else
    %     fMc = fMc0;
    % end
report_this_filefun(mfilename('fullpath'));
    date_matlab = datenum(mycat.Date);
    date_main = datenum(ZG.maepi.Date);
    time_aftershock = date_matlab-date_main;

    l = time_aftershock(:) > 0;
    tas = time_aftershock(l);
    eqcatalogue = mycat.subset(l);

    % Estimation of Omori parameters from learning period
    l = tas <= time;% & mycat.Magnitude >= fMc;
    time_as=tas(l);
    % [pval, pvalstd, cval, cvalstd, kval, kvalstd, loopout] = brutebootF(time_as, bootloops);
    % Calculate uncertainties and mean values of p,c,and k
    [pval, pvalstd, cval, cvalstd, kval, kvalstd, loopout] = brutebootloglike(time_as,bootloops);
    % Calculate p,c,k for real dataset
    [pval, cval, kval] = bruteforceloglike(sort(time_as));
    pval = round(100*pval)/100;
    cval = round(100*cval)/100;
    kval = round(10*kval)/10;

    [H,P,KSSTAT,fRMS] = calc_llkstest(time_as, pval, cval, kval);
    if H==1
        rc = [NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN];
        return
    end

    if isnan(pval) == 0
        lf = tas <= time+timef ;
        time_asf= [tas(lf) ];
        %     nummod = length(time_asf); % forecasted number of aftershocks in time window
        if pval == 1
            pv = 1-10^(-6);
        else
            pv = pval;
        end

        cumnrf = (1:length(time_asf))';
        cumnr_modelf = [];
        for i=1:length(time_asf)
            if pval ~= 1
                cm = kval/(pval-1)*(cval^(1-pval)-(time_asf(i)+cval)^(1-pval));
            else
                cm = kval*log(time_asf(i)/cval+1);
            end
            cumnr_modelf = [cumnr_modelf; cm];
        end
        %     t_forecast = (-nummod*(pv-1)/kval+(time+cval)^(1-pv))^(1/(1-pv))-cval; % forecast interval
        %     if isreal(t_forecast) == 1 & t_forecast > time
        %         if t_forecast > max(time_as)
        %             return
        %         end
        %l = time_as <= t_forecast & mycat.Magnitude >= fMc & time_as >= time;
        %l = tas <= time+timef;
        % Find amount of events in forecast period for modeled data
        nummod = max(cumnr_modelf)-cumnr_modelf(length(time_as));
        % Find amount of  events in forecast period for observed data
        l = time_asf <=time+timef & time_asf > time;
        numreal = sum(l); % observed number of aftershocks
        absdiff = numreal-nummod;

        % calculate uncertainty sigma in forecasted number of aftershocks by Fehlerfortpflanzungsgesetz
        %time1 = t_forecast;
        time1 = time+timef;
        mpm1 = 1-pv;
        t1c = time1+cval;
        t0c = time+cval;
        sigma = (((-t1c^mpm1+t0c^mpm1)/(pv-1)*kvalstd)^2+...
            (kval/(pv-1)*(-t1c^mpm1*mpm1/t1c+t0c^mpm1*mpm1/t0c)*cvalstd)^2+...
            (kval/(pv-1)*(t1c^mpm1*log(t1c)+t1c^mpm1/(pv-1)-t0c^mpm1*log(t0c)-t0c^mpm1/(pv-1))*pvalstd)^2)^0.5;

        %rc = [rc; time absdiff sigma numreal nummod pval pvalstd cval cvalstd kval kvalstd t_forecast fMc];

        % Compute 2nd moment of forecasted number of events at time
        for j = 1:length(loopout(:,1))
            cumnr = (1:length(time_asf))';
            cumnr_model = [];
            pvalb = loopout(j,1);
            cvalb = loopout(j,2);
            kvalb = loopout(j,3);
            for i=1:length(time_asf)
                if pval ~= 1
                    cm = kvalb/(pvalb-1)*(cvalb^(1-pvalb)-(time_asf(i)+cvalb)^(1-pvalb));
                else
                    cm = kvalb*log(time_asf(i)/cvalb+1);
                end
                cumnr_model = [cumnr_model; cm];
            end
            loopout(j,4) = max(cumnr_model);
        end
        fStdBst = calc_StdDev(loopout(:,4));
        rc = [time absdiff sigma numreal nummod pval pvalstd cval cvalstd kval kvalstd fStdBst];
    else
        rc = [NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN];
    end

author: CelsoReyes


text:
function [swav]=s_wav(x,a_swav,d_swav,t_swav,li)
l=li;
x=x-t_swav;
a=a_swav;
b=(2*l)/d_swav;
n=100;
s1=(a/(2*b))*(2-b);
s2=0;
for i = 1:n
    harm3=(((2*b*a)/(i*i*pi*pi))*(1-cos((i*pi)/b)))*cos((i*pi*x)/l);
    s2=s2+harm3;
end
swav=-1*(s1+s2);

author: Sable

repo: mcbench-benchmarks


text:
function [minFlux, maxFlux, optsol, ret, fbasol, fvamin, fvamax, statussolmin, statussolmax] = fastFVA(model, optPercentage, osenseStr, solverName, rxnsList, matrixAS, cpxControl, strategy, rxnsOptMode, printLevel)
% Flux variablity analysis optimized for the CPLEX solver.
% Solves LPs of the form:
%
% .. math::
%
%    \forall ~ v_j: ~&~ max/min ~&~ v_j\\
%                 ~&~ s.t.    ~&~ Sv = b\\
%                 ~&~         ~&~ l_b \leq v \leq u_b
%
% If the optional fields are supplied, following LPs are solved
%
% .. math::
%
%    \forall ~ v_j: ~&~ max/min ~&~ v_j\\
%                 ~&~ s.t.    ~&~ Av (c_sense) b\\
%                 ~&~         ~&~ l_b \leq v \leq u_b
%
% fastFVA returns vectors for the initial FBA in FBASOL together with matrices FVAMIN and
% FVAMAX containing the flux values for each individual min/max problem.
% Note that for large models the memory requirements may become prohibitive.
% To save large `fvamin` and `fvamax` matrices, toggle v7.3 in Preferences -> General -> MAT-Files
%
% If a `rxnsList` vector is specified then only the corresponding entries in
% minFlux and maxFlux are defined (all remaining entries are zero).
%
% USAGE:
%
%    [minFlux, maxFlux, optsol, ret, fbasol, fvamin, fvamax, statussolmin, statussolmax] = fastFVA(model, optPercentage, osenseStr, solverName, rxnsList, matrixAS, cpxControl, strategy, rxnsOptMode)
%
% INPUTS:
%   model:             COBRA model structure
%
%                        * .S - (required) Stoichiometric matrix
%                        * .b - (required) Right hand side vector
%                        * .c - (required) Objective coefficients
%                        * .lb - (required) Lower bounds
%                        * .ub - (required) Upper bounds
%                        * .A - (optional) Stoichiometric matrix (with constraints)
%                        * .csense - (optional) Type of constraints, `csense` is a vector with elements `E` (equal), `L` (less than) or `G` (greater than).
%   optPercentage:     Only consider solutions that give you at least a certain
%                      percentage of the optimal solution (default = `100`, equivalent to optimal solutions only)
%   osenseStr:         Objective ('min' or 'max') (default 'max')
%   solverName:        name of the solver, default: `ibm_cplex`
%
% OPTIONAL INPUTS:
%   matrixAS:          `A` or `S` - choice of the model matrix, coupled (A) or uncoupled (S)
%   cpxControl:        Parameter set of CPLEX loaded externally
%   rxnsList:          List of reactions to analyze (default all rxns, i.e. `1:length(model.rxns)``)
%   strategy:          Paralell distribution strategy of reactions among workers
%
%                       * 0 = Blind splitting: default random distribution
%                       * 1 = Extremal dense-and-sparse splitting: every worker receives dense and sparse reactions, starting from both extremal indices of the sorted column density vector
%                       * 2 = Central dense-and-sparse splitting: every worker receives dense and sparse reactions, starting from the beginning and center indices of the sorted column density vector
%   rxnsOptMode:       List of min/max optimizations to perform:
%                       * 0 = only minimization;
%                       * 1 = only maximization;
%                       * 2 = minimization & maximization;
%   printLevel:        Verbose level (default: 1)
%                       * 0 = mute
%                       * 1 = default
%
% OUTPUTS:
%   minFlux:           Minimum flux for each reaction
%   maxFlux:           Maximum flux for each reaction
%   optsol:            Optimal solution (of the initial FBA)
%   ret:               Zero if success (global return code from FVA)
%
% OPTIONAL OUTPUTS:
%   fbasol:            Initial FBA in FBASOL
%   fvamin:            matrix with flux values for the minimization problem
%   fvamax:            matrix with flux values for the maximization problem
%   statussolmin:      vector of solution status for each reaction (minimization)
%   statussolmax:      vector of solution status for each reaction (maximization)
%
%
% EXAMPLE:
%
%    load modelRecon1Biomass.mat % Human reconstruction network (Recon1)
%    setWorkerCount(4);
%    [minFlux,maxFlux] = fastFVA(model, 90);
%
% NOTE:
%
%    S. Gudmundsson and I. Thiele, Computationally efficient
%    Flux Variability Analysis. BMC Bioinformatics, 2010, 11:489
%
% NOTE:
%
%    * Matlab R2014a fully tested on UNIX and DOS Systems
%    * Matlab R2015b throws compatibility errors with CPLEX 12.6.3 on DOS Systems
%    * Matlab R2016b and the MinGW64 compiler are not compatible with the CPLEX 12.6.3 library
%
%    The version of fastFVA only supports the CPLEX solver. The code has been tested with
%    CPLEX 12.6.2, 12.6.3, 12.7.0 and 12.7.1. Install
%    CPLEX (64-bit) as explained `here`_.
%    A particular interface, such as TOMLAB, is not needed in order to run fastFVA.
%    Please note that only the 64-bit versions of CPLEX 12.7.1 are supported.
%    In order to run the code on 32-bit systems, the appropriate MEX files need to be generated
%    using generateMexFastFVA().
%
% .. _here: https://opencobra.github.io/cobratoolbox/docs/solvers.html
%
% .. Authors:
%       - Original author: Steinn Gudmundsson.
%       - Contributor: Laurent Heirendt, LCSB
%
% .. Last updated: October 2016

global CBTDIR

% save the userpath
originalUserPath = path;

% save the current directory
currentDir = pwd;

% the root path must be the root directory as the path to the logFiles is hard-coded
cd(CBTDIR);

% determine the printLevel
if (nargin < 10 || isempty(printLevel))
    printLevel = 1;
end

% determine the latest installed CPLEX version
cplexVersion = getCobraSolverVersion('ibm_cplex', printLevel);

% check if the provided fastFVA binaries are compatible with the current system configuration
[newestCplexBin, throwBinGenerationError]= checkFastFVAbin(cplexVersion);

% set a random log filename to avoid overwriting ongoing runs
rng('shuffle');
filenameParfor = ['parfor_progress_', datestr(now, 30), '_', num2str(randi(9)), '.txt'];
filenameParfor = [CBTDIR filesep '.tmp' filesep filenameParfor];

% turn on the load balancing for large problems
loadBalancing = 0;  % 0: off; 1: on

% define if information about the work load distriibution will be shown or not
showSplitting = 1;

% define the input arguments
if (nargin < 8 || isempty(strategy))
    strategy = 0;
end
if (nargin < 7 || isempty(cpxControl))
    cpxControl = struct([]);
end
if (nargin < 6 || isempty(matrixAS)) || ~isfield(model,'C')
    matrixAS = 'S';
end
if (nargin < 5 || isempty(rxnsList))
    rxns = 1:length(model.rxns);
    rxnsList = model.rxns;
else
    % check here if the vector of rxns is sorted or not
    % this needs to be fixed to sort the flux vectors accordingly
    % as the find() function sorts the reactions automatically
    % ->> this is currently an issue on git
    [~, indexRxns] = ismember(model.rxns, rxnsList);
    nonZeroIndices = [];
    for i = 1:length(indexRxns)
        if indexRxns(i) ~= 0
            nonZeroIndices = [nonZeroIndices, indexRxns(i)];
        end
    end
    if issorted(nonZeroIndices) == 0
        error('\n-- ERROR:: Your input reaction vector is not sorted. Please sort your reaction vector first.\n\n')
    end

    rxns = find(ismember(model.rxns, rxnsList))';  % transpose rxns

end
if (nargin < 9 || isempty(rxnsOptMode))
    rxnsOptMode = 2 * ones(length(rxns), 1)';  % status = 2 (min & max) for all reactions
end
if (nargin < 4 || isempty(solverName))
    solverName = 'ibm_cplex';
end
if (nargin < 3 || isempty(osenseStr))
    osenseStr = 'max';
end
if (nargin < 2 || isempty(optPercentage))
    optPercentage = 100;
end

% define extra outputs if required
if nargout > 4 && nargout <= 7
    assert(nargout == 7);
    bExtraOutputs = true;
else
    bExtraOutputs = false;
end

% define extra outputs if required
if nargout > 7
    assert(nargout == 9);
    bExtraOutputs1 = true;
else
    bExtraOutputs1 = false;
end

% print a warning when output arguments are not defined.
if nargout ~= 4 && nargout ~= 7 && nargout ~= 9
    if printLevel > 0
        fprintf('\n-- Warning:: You may only ouput 4, 7 or 9 variables.\n\n')
    end
end

% define the osenseStr
if strcmpi(osenseStr, 'max')
    obj = -1;
elseif strcmpi(osenseStr, 'min')
    obj = 1;
else
    error('Unknown osenseStr');
end

% define the solverName
if strcmp('glpk', solverName)
    error('ERROR : GLPK is not (yet) supported as the binaries are not yet available.')
elseif strcmp('ibm_cplex', solverName)
    if throwBinGenerationError
        %attempt to use the newest cplex binaries available, even though
        %they may not work, it is better than nothing
        FVAc = str2func(['cplexFVA' newestCplexBin]);
    else
        FVAc = str2func(['cplexFVA' cplexVersion]);
    end
else
    error(['Solver ', solverName, ' not supported.']);
end
% define the CPLEX parameter set and the associated values - split the struct
namesCPLEXparams = fieldnames(cpxControl);
nCPLEXparams = length(namesCPLEXparams);
valuesCPLEXparams = zeros(nCPLEXparams, 1);
for i = 1:nCPLEXparams
    valuesCPLEXparams(i) = getfield(cpxControl, namesCPLEXparams{i});
end

% create an LP problem
LPproblem = buildLPproblemFromModel(model);

% define the stoichiometric matrix to be solved
if matrixAS == 'A'
    [A,b,csense,lb,ub,c] = deal(LPproblem.A,LPproblem.b,LPproblem.csense,LPproblem.lb,LPproblem.ub,LPproblem.c);
    if printLevel > 0
        fprintf(' >> Solving Model.A. (coupled) - Generalized\n');
    end
else
    if ~isfield(model,'csense')
        model.csense = repmat('E',size(model.S,1),1);
    end
    [A,b,csense,lb,ub,c] = deal(model.S,model.b,model.csense,model.lb,model.ub,model.c);
    if printLevel > 0
        fprintf(' >> Solving Model.S. (uncoupled) \n');
    end
end

if printLevel > 0
    fprintf(' >> The number of arguments is: input: %d, output %d.\n', nargin, nargout);
end

% define the matrix A as sparse in case it is not as
% C code assumes a sparse stochiometric matrix
if ~issparse(A)
    A = sparse(A);
end

% determine the size of the stoichiometric matrix
[m, n] = size(A);
if printLevel > 0
    fprintf(' >> Size of stoichiometric matrix: (%d,%d)\n', m, n);
end

% determine the number of reactions that are considered
nR = length(rxns);
if nR ~= n
    if printLevel > 0
        fprintf(' >> Only %d reactions of %d are solved (~ %1.2f%%).\n', nR, n, nR * 100 / n);
    end
    n = nR;
else
    if printLevel > 0
        fprintf(' >> All reactions are solved (%d reactions - 100%%).\n', n);
    end
end

% output how many reactions are min, max, or both
totalOptMode = length(find(rxnsOptMode == 0));
if totalOptMode == 1
    if printLevel > 0
        fprintf(' >> %d reaction out of %d is minimized (%1.2f%%).\n', totalOptMode, n, totalOptMode * 100 / n);
    end
else
    if printLevel > 0
        fprintf(' >> %d reactions out of %d are minimized (%1.2f%%).\n', totalOptMode, n, totalOptMode * 100 / n);
    end
end

totalOptMode = length(find(rxnsOptMode == 1));
if totalOptMode == 1
    if printLevel > 0
        fprintf(' >> %d reaction out of %d is maximized (%1.2f%%).\n', totalOptMode, n, totalOptMode * 100 / n);
    end
else
    if printLevel > 0
        fprintf(' >> %d reactions out of %d are maximized (%1.2f%%).\n', totalOptMode, n, totalOptMode * 100 / n);
    end
end

totalOptMode = length(find(rxnsOptMode == 2));
if totalOptMode == 1
    if printLevel > 0
        fprintf(' >> %d reaction out of %d is minimized and maximized (%1.2f%%).\n', totalOptMode, n, totalOptMode * 100 / n);
    end
else
    if printLevel > 0
        fprintf(' >> %d reactions out of %d are minimized and maximized (%1.2f%%).\n', totalOptMode, n, totalOptMode * 100 / n);
    end
end

% count the number of workers
poolobj = gcp('nocreate');  % If no pool, do not create new one.
if isempty(poolobj)
    nworkers = 1;
else
    nworkers = poolobj.NumWorkers;
end

% creates the directory where the log files will be generated
rootDirFastFVA = [CBTDIR filesep 'src' filesep 'analysis' filesep 'FVA' filesep 'fastFVA'];

% define the directory to the logFiles and results directories
logFileDir = [rootDirFastFVA filesep 'logFiles'];

% initialisations
istart(1) = 1;
maxFluxTmp = {};
minFluxTmp = {};

% launch fastFVA on 1 core
if nworkers <= 1

    % define the end of the index vector
    iend(1) = n;

    if ~isempty(rxnsList)
        rxnsKey = find(ismember(model.rxns, rxnsList));
    else
        rxnsKey = (1:n);
    end

    % sequential version
    if printLevel > 0
        fprintf(' \n WARNING: The Sequential Version might take a long time.\n\n');
    end
    if bExtraOutputs1
        [minFlux, maxFlux, optsol, ret, fbasol, fvamin, fvamax, statussolmin, statussolmax] = FVAc(c, A, b, csense, lb, ub, ...
                                                                                optPercentage, obj, rxnsKey, ...
                                                                                1, cpxControl, valuesCPLEXparams, rxnsOptMode, logFileDir, printLevel);
    elseif bExtraOutputs
        [minFlux, maxFlux, optsol, ret, fbasol, fvamin, fvamax] = FVAc(c, A, b, csense, lb, ub, ...
                                                                       optPercentage, obj, rxnsKey, ...
                                                                       1, cpxControl, valuesCPLEXparams, rxnsOptMode, logFileDir, printLevel);
    else
        [minFlux, maxFlux, optsol, ret] = FVAc(c, A, b, csense, lb, ub, ...
                                               optPercentage, obj, rxnsKey, ...
                                               1, cpxControl, valuesCPLEXparams, rxnsOptMode, logFileDir, printLevel);
    end

    if ret ~= 0 && printLevel > 0
        if printLevel > 0
            fprintf('Unable to complete the FVA, return code=%d\n', ret);
        end
    end

    % output the minimum and maximum fluxes
    minFluxTmp{1} = minFlux;
    maxFluxTmp{1} = maxFlux;

    % output the minimum and maximum flux vectors
    if bExtraOutputs || bExtraOutputs1
        fvaminRes{1} = fvamin;
        fvamaxRes{1} = fvamax;
        fbasolRes{1} = fbasol;
    end

    % output the solver status
    if bExtraOutputs1
        statussolminRes{1} = statussolmin;
        statussolmaxRes{1} = statussolmax;
    end
else
    % divide the reactions amongst workers

    % Note:
    %      The load balancing can be improved for certain problems, e.g. in case
    %      of problems involving E-type matrices, some workers will get mostly
    %      well-behaved LPs while others may get many badly scaled LPs.

    if n > 5000 && loadBalancing == 1
        % primitive load-balancing strategy for large problems
        nworkers = 4 * nworkers;
        if printLevel > 0
            fprintf(' >> The load is balanced and the number of virtual workers is %d.\n', nworkers);
        end
    end

    nrxn = repmat(fix(n / nworkers), nworkers, 1);
    i = 1;
    while sum(nrxn) < n
        nrxn(i) = nrxn(i) + 1;
        i = i + 1;
    end

    [Nmets, Nrxns] = size(A);
    assert(sum(nrxn) == n);
    iend(1) = nrxn(1);
    for i = 2:nworkers
        istart(i) = iend(i - 1) + 1;
        iend(i) = istart(i) + nrxn(i) - 1;
    end

    startMarker1 = istart;
    endMarker1 = iend;

    startMarker2 = istart;
    endMarker2 = iend;

    % calculate the column density and row density
    NrxnsList = length(rxnsList);

    % initialize the column density vector
    cdVect = zeros(NrxnsList, 1);

    for i = 1:NrxnsList
        tmpRxnID = findRxnIDs(model, rxnsList(i));

        columnDensity = nnz(A(:, tmpRxnID));
        columnDensity = columnDensity / Nmets * 100;
        cdVect(i) = columnDensity;
    end

    [~, indexcdVect] = sort(cdVect, 'descend');

    rxnsVect = linspace(1, NrxnsList, NrxnsList);

    sortedrxnsVect = rxnsVect(indexcdVect);

    if strategy > 0
        pRxnsHalfWorker = ceil(NrxnsList / (2 * nworkers));

        for i = 1:nworkers

            startMarker1(i) = (i - 1) * pRxnsHalfWorker + 1;
            endMarker1(i) = i * pRxnsHalfWorker;

            if strategy == 1
                startMarker2(i) = startMarker1(i) + ceil(NrxnsList / 2);
                endMarker2(i) = endMarker1(i) + ceil(NrxnsList / 2);
            elseif strategy == 2
                startMarker2(i) = ceil(NrxnsList / 2) + startMarker1(i);
                endMarker2(i) = startMarker2(i) + pRxnsHalfWorker + 1;
            end

            % avoid start indices beyond the total number of reactions
            if startMarker1(i) > NrxnsList
                startMarker1(i) = NrxnsList;
            end
            if startMarker2(i) > NrxnsList
                startMarker2(i) = NrxnsList;
            end

            % avoid end indices beyond the total number of reactions
            if endMarker1(i) > NrxnsList
                endMarker1(i) = NrxnsList;
            end
            if endMarker2(i) > NrxnsList
                endMarker2(i) = NrxnsList;
            end

            % avoid flipped chunks
            if startMarker1(i) > endMarker1(i)
                startMarker1(i) = endMarker1(i);
            end
            if startMarker2(i) > endMarker2(i)
                startMarker2(i) = endMarker2(i);
            end
        end
    end

    minFlux = zeros(length(model.rxns), 1);
    maxFlux = zeros(length(model.rxns), 1);
    iopt = zeros(nworkers, 1);
    iret = zeros(nworkers, 1);

    % initialize extra outputs
    if bExtraOutputs || bExtraOutputs1
        fvaminRes = {};
        fvamaxRes = {};
        fbasolRes = {};
    end

    if bExtraOutputs1
        statussolminRes = {};
        statussolmaxRes = {};
    end

    if printLevel > 0
        fprintf('\n -- Starting to loop through the %d workers. -- \n', nworkers);
        fprintf('\n -- The splitting strategy is %d. -- \n', strategy);
    end

    out = parfor_progress(nworkers, filenameParfor);

    parfor i = 1:nworkers

        rxnsKey = 0;  % silence warning

        % preparation of reactionKey
        if strategy == 1 || strategy == 2
            rxnsKey = [sortedrxnsVect(startMarker1(i):endMarker1(i)), sortedrxnsVect(startMarker2(i):endMarker2(i))];
        else
            rxnsKey = rxns(istart(i):iend(i));
        end

        t = getCurrentTask();
        if printLevel > 0
            fprintf('\n----------------------------------------------------------------------------------\n');
        end
        if strategy == 0
            if printLevel > 0
                fprintf('--  Task Launched // TaskID: %d / %d (LoopID = %d) <> [%d, %d] / [%d, %d].\n', ...
                    t.ID, nworkers, i, istart(i), iend(i), m, n);
            end
        else
            if printLevel > 0
                fprintf('--  Task Launched // TaskID: %d / %d (LoopID = %d) <> [%d:%d] & [%d:%d] / [%d, %d].\n', ...
                    t.ID, nworkers, i, startMarker1(i), endMarker1(i), ...
                    startMarker2(i), endMarker2(i), m, n);
            end
        end

        tstart = tic;

        minf = zeros(length(model.rxns), 1);
        maxf = zeros(length(model.rxns), 1);
        fvamin_single = 0; fvamax_single = 0; fbasol_single = 0; statussolmin_single = 0; statussolmax_single = 0;  % silence warnings

        if bExtraOutputs1
            [minf, maxf, iopt(i), iret(i), fbasol_single, fvamin_single, fvamax_single, ...
             statussolmin_single, statussolmax_single] = FVAc(model.c, A, b, csense, model.lb, model.ub, ...
                                                              optPercentage, obj, rxnsKey', ...
                                                              t.ID, cpxControl, valuesCPLEXparams, ...
                                                              rxnsOptMode(istart(i):iend(i)), logFileDir, printLevel);
        elseif bExtraOutputs
            [minf, maxf, iopt(i), iret(i), fbasol_single, fvamin_single, fvamax_single] = FVAc(model.c, A, b, csense, ...
                                                                                model.lb, model.ub, ...
                                                                                optPercentage, obj, rxnsKey', ...
                                                                                t.ID, cpxControl, valuesCPLEXparams, ...
                                                                                rxnsOptMode(istart(i):iend(i)), ...
                                                                                logFileDir, printLevel);
        else
            if printLevel > 0
                fprintf(' >> Number of reactions given to the worker: %d \n', length(rxnsKey));
            end

            [minf, maxf, iopt(i), iret(i)] = FVAc(model.c, A, b, csense, model.lb, model.ub, ...
                                                  optPercentage, obj, rxnsKey', ...
                                                  t.ID, cpxControl, valuesCPLEXparams, ...
                                                  rxnsOptMode(istart(i):iend(i)), logFileDir, printLevel);
        end
        if printLevel > 0
            fprintf(' >> Time spent in FVAc: %1.1f seconds.', toc(tstart));
        end

        if iret(i) ~= 0 && printLevel > 0
            if printLevel > 0
                fprintf('Problems solving partition %d, return code=%d\n', i, iret(i))
            end
        end

        minFluxTmp{i} = minf;
        maxFluxTmp{i} = maxf;

        if bExtraOutputs || bExtraOutputs1
            fvaminRes{i} = fvamin_single;
            fvamaxRes{i} = fvamax_single;
            fbasolRes{i} = fbasol_single;
        end

        if bExtraOutputs1
            statussolminRes{i} = statussolmin_single;
            statussolmaxRes{i} = statussolmax_single;
        end
        if printLevel > 0
            fprintf('\n----------------------------------------------------------------------------------\n');
        end

        % print out the percentage of the progress
        percout = parfor_progress(-1, filenameParfor);

        if percout < 100
            if printLevel > 0
                fprintf(' ==> %1.1f%% done. Please wait ...\n', percout);
            end
        else
            if printLevel > 0
                fprintf(' ==> 100%% done. Analysis completed.\n', percout);
            end
        end
    end

    % aggregate results
    optsol = iopt(1);
    ret = max(iret);
    out = parfor_progress(0, filenameParfor);
end

% aggregate the results for the maximum and minimum flux vectors
for i = 1:nworkers
    % preparation of reactionKey
    if strategy == 1 || strategy == 2
        indices = [sortedrxnsVect(startMarker1(i):endMarker1(i)), sortedrxnsVect(startMarker2(i):endMarker2(i))];
    else
        indices = rxns(istart(i):iend(i));
    end

    % store the minFlux
    tmp = maxFluxTmp{i};
    maxFlux(indices, 1) = tmp(indices);

    % store the maxFlux
    tmp = minFluxTmp{i};
    minFlux(indices, 1) = tmp(indices);
end

if bExtraOutputs || bExtraOutputs1

    if nworkers > 1
        fbasol = fbasolRes{1};  % initial FBA solutions are identical across workers
    end

    fvamin = zeros(length(model.rxns), length(model.rxns));
    fvamax = zeros(length(model.rxns), length(model.rxns));

    if nworkers > 1
        if bExtraOutputs1
            statussolmin = -1 + zeros(length(model.rxns), 1);
            statussolmax = -1 + zeros(length(model.rxns), 1);
        end
    end

    for i = 1:nworkers
        % preparation of reactionKey
        if strategy == 1 || strategy == 2
            indices = [sortedrxnsVect(startMarker1(i):endMarker1(i)), sortedrxnsVect(startMarker2(i):endMarker2(i))];
        else
            indices = rxns(istart(i):iend(i));
        end

        fvamin(:, indices) = fvaminRes{i};
        fvamax(:, indices) = fvamaxRes{i};

        if bExtraOutputs1
            tmp = statussolminRes{i}';
            statussolmin(indices, 1) = tmp(indices);
            tmp = statussolmaxRes{i}';
            statussolmax(indices, 1) = tmp(indices);
        end
    end
end

if strategy == 0 && ~isempty(rxnsList)
    if bExtraOutputs || bExtraOutputs1
        fvamin = fvamin(:, rxns);  % keep only nonzero columns
        fvamax = fvamax(:, rxns);
    end

    if bExtraOutputs1
        statussolmin = statussolmin(rxns);
        statussolmax = statussolmax(rxns);
    end

    minFlux(find(~ismember(model.rxns, rxnsList))) = [];
    maxFlux(find(~ismember(model.rxns, rxnsList))) = [];
end

% restore the original path
path(originalUserPath);
addpath(originalUserPath);

% change back to the original directory
cd(currentDir);

function [newestCplexBin, throwBinGenerationError] = checkFastFVAbin(cplexVersion)
% determine the version of the CPLEX binaries by
% browsing to the folder with the binaries and the .tmp folder (if it exists)
% and retrieve all versions
%
% USAGE:
%    checkFastFVAbin(cplexVersion)
%
% INPUT:
%    cplexVersion:    CPLEX version (string), obtained using `getCobraSolverVersion`
%

global CBTDIR

% retrieve the contents of the binary directory (architecture dependent)
aDir = [CBTDIR filesep 'binary' filesep computer('arch') filesep 'bin' filesep 'fastFVA'];
if 1 %debug
    cd(aDir);
end

d1 = dir(aDir);

% include CPLEX binaries that already have been generated using generateMexFastFVA
tmpDir = [CBTDIR filesep '.tmp'];
d2 = dir(tmpDir);

% create .tmp if not already present
if ~exist(tmpDir, 'dir')
    mkdir(tmpDir);
end

% concatenate both directories
d = {d1; d2};

for p = 1:length(d)
    tmpD = d{p};
    k = 1;
    for i = 1:numel(tmpD)
        if contains(tmpD(i).name, 'cplexFVA') && ~strcmpi(tmpD(i).name, '.') && ~strcmpi(tmpD(i).name, '..') && isempty(strfind(tmpD(i).name, '.txt'))
            tmpName = tmpD(i).name;
            tmpNameSplit = strsplit(tmpName, '.');
            tmpName = tmpNameSplit{1};
            %look for the binary versio
            binVersion{k} = tmpName(length('cplexFVA')+1:end);
            k = k + 1;
        end
    end
end

%compare the binary version of cplexFVA with the version of cplex
throwBinGenerationError = false;
for k = 1:length(binVersion)
    if ~strcmpi(cplexVersion, binVersion{k})
        throwBinGenerationError = true;
        kp = k;
    end
end

if throwBinGenerationError
    maxLenBinVersion=0;
    for k = 1:length(binVersion)
        if maxLenBinVersion<length(binVersion{k})
            maxLenBinVersion=length(binVersion{k});
        end
    end
    for k = 1:length(binVersion)
        tmp(1:maxLenBinVersion)='0';
        tmp(1:length(binVersion{k}))=binVersion{k};
        equalLenbinVersion(k)=str2num(tmp);
    end
    [equalLenbinVersionMax, equalLenbinVersionMaxInd]=max(equalLenbinVersion);
    newestCplexBin=binVersion{equalLenbinVersionMaxInd};

    if 0
    error(['Official binaries are only available for CPLEX version ', binVersion{kp}, '. ', ...
            'You have installed version ', cplexVersion, '. Please run: ', ...
            '>> generateMexFastFVA() in order to generate a new binary file.']);
    else
     warning(['The most recent cplexFVA binaries are only tested for CPLEX version ', newestCplexBin, '. ', ...
            'However, you have installed version ', cplexVersion, '. If you experience an error using fastFVA, you need to run: ', ...
            '>> generateMexFastFVA() in order to generate a new cplexFVA binary file that matches cplex version' cplexVersion]);
    end
end

author: opencobra

repo: cobratoolbox


function[box,flag] = getBox2(theta,thetaDot,x,xDot)
theta = wrapTo180(rad2deg(theta));
thetaDot = rad2deg(thetaDot);
flag = 0;
if(theta>=-180&&theta<-150)
    thetaBucket = 1;
    flag = -1;
elseif(theta>=-150&&theta<-120)
    thetaBucket = 2;
    flag = -1;
elseif(theta>=-120&&theta<-90)
    thetaBucket = 3;
    flag = -1;
elseif(theta>=-90&&theta<-60)
    thetaBucket = 4;
    flag = -1;
elseif(theta>=-60&&theta<-30)
    thetaBucket = 5;
    flag = -1;
elseif(theta>=-30&&theta<-24)
    thetaBucket = 6;
    flag = -1;
elseif(theta>=-24&&theta<-18)
    thetaBucket = 7;
    flag = -1;
elseif(theta>=-18&&theta<-12)
    thetaBucket = 8;
    flag = -1;
elseif(theta>=-12&&theta<-6)
    thetaBucket = 9;
elseif(theta>=-6&&theta<-1)
    thetaBucket = 10;
elseif(theta>=-1&&theta<0)
    thetaBucket = 11;
elseif(theta>=0&&theta<1)
    thetaBucket  =12;
elseif(theta>=1&&theta<6)
    thetaBucket = 13;
elseif(theta>=6&&theta<12)
    thetaBucket = 14;
elseif(theta>=12&&theta<18)
    thetaBucket = 15;
    flag = -1;
elseif(theta>=18&&theta<24)
    thetaBucket = 16;
    flag = -1;
elseif(theta>=24&&theta<30)
    thetaBucket = 17;
    flag = -1;
elseif(theta>=30&&theta<60)
    thetaBucket = 18;
    flag = -1;
elseif(theta>=60&&theta<90)
    thetaBucket = 19;
    flag = -1;
elseif(theta>=90&&theta<120)
    thetaBucket = 20;
    flag = -1;
elseif(theta>=120&&theta<150)
    thetaBucket = 21;
    flag = -1;
elseif(theta>=150&&theta<=180)
    thetaBucket = 22;
    flag = -1;
end
if (x < -2.4 || x > 2.4)
    flag = -1;
    xBucket = 1;
end
if (x<-0.8&&x>=-2.4)
	xBucket = 1;
elseif (x<=0.8&&x>=-0.8)
	xBucket = 2;
elseif (x<=2.4&&x>0.8)
	xBucket = 3;
end

if (xDot<-0.5)
	xDotBucket = 1;
elseif (xDot>=-0.5&&xDot<=0.5)
	xDotBucket = 2;
else
	xDotBucket = 3;
end

if (thetaDot<-50)
	thetaDotBucket = 1;
elseif (thetaDot>=-50&&thetaDot<=50)
	thetaDotBucket = 2;
else
	thetaDotBucket = 3;
end

box = sub2ind([22,3,3,3],thetaBucket,thetaDotBucket,xBucket,xDotBucket);


text:
% writelocs() - write a file containing channel location, type and gain information
%
% Usage:
%   >> writelocs( chanstruct, filename );
%   >> writelocs( chanstruct, filename, 'key', 'val' );
%
% Inputs:
%   chanstruct - EEG.chanlocs data structure returned by readlocs() containing
%                channel location, type and gain information.
%   filename   - File name for saving channel location, type and gain information
%
% Optional inputs:
%   'filetype'  - ['loc'|'sph'|'sfp'|'xyz'|'polhemus'|'besa'|'chanedit'|'custom']
%                 Type of the file to write. By default the file type is indicated
%                 by the file extension.
%                  'loc' - An EEGLAB 2-D polar coordinates channel locations file
%                          Coordinates are theta and radius (see definitions below).
%                  'sph' - A Matlab spherical coordinates file (Note: spherical
%                          coordinates used by Matlab functions are different
%                          from spherical coordinates used in BESA - see below).
%                  'sfp' - EGI cartesian coordinates (not Matlab cartesian - see below).
%                  'xyz' - MATLAB/EEGLAB cartesian coordinates (Not EGI cartesian;
%                          z is toward nose; y is toward left ear; z is toward vertex).
%                  'polhemus' or 'polhemusx' - Polhemus electrode location file recorded with
%                          'X' on sensor pointing to subject (see below and readelp()).
%                  'polhemusy' - Polhemus electrode location file recorded with
%                          'Y' on sensor pointing to subject (see below and readelp()).
%                  'besa' - BESA'(.elp') spherical coordinate file. (Not MATLAB spherical
%                           - see below).
%                  'chanedit' - EEGLAB channel location files created by pop_chanedit().
%                  'custom' - Ascii files with columns in user-defined 'format' (see below).
%   'format'    - [cell array] Format of a 'custom' channel location file (see above).
%                          Default if no file type is defined. The cell array contains
%                          labels defining the meaning of each column of the input file.
%                           'channum'   [positive integer] channel number
%                           'labels'    [string] channel name (no spaces)
%                           'theta'     [real degrees] 2-D angle in polar coordinates.
%                                       positive => rotating from nose (0) toward left ear
%                           'radius'    [real] radius in 2-D polar coords (0.5 is disk limits)
%                           'X'         [real] Matlab-cartesian X coordinate (to nose)
%                           'Y'         [real] Matlab-cartesian Y coordinate (to left ear)
%                           'Z'         [real] Matlab-cartesian Z coordinate (to vertex)
%                           '-X','-Y','-Z' Matlab-cartesian coordinates pointing away from above
%                           'sph_theta' [real degrees] Matlab spherical horizontal angle.
%                                       positive => rotating from nose (0) toward left ear.
%                           'sph_phi'   [real degrees] Matlab spherical elevation angle;
%                                       positive => rotating from horizontal (0) upwards.
%                           'sph_radius' [real] distance from head center (unused)
%                           'sph_phi_besa' [real degrees] BESA phi angle from vertical.
%                                       positive => rotating from vertex (0) towards right ear.
%                           'sph_theta_besa' [real degrees] BESA theta horiz/azimuthal angle.
%                                       positive => rotating from right ear (0) toward nose.
%     The input file may also contain other channel information fields
%                           'type'      channel type: 'EEG', 'MEG', 'EMG', 'ECG', others ...
%                           'calib'     [real near 1.0] channel calibration value.
%                           'gain'      [real > 1] channel gain.
%                           'custom1'   custom field #1.
%                           'custom2', 'custom3', 'custom4' more custom fields.
%   'unicoord'     - ['on'|'off'] Uniformize all coordinates. Default 'on'.
%   'header'       - ['on'|'off'] Add a header comment line with the name of each column.
%                             Comment lines begin with '%'. Default is 'off'.
%   'customheader' - [string] Add a custom header at the beginning of the file and
%                             preceded by '%'.  If used with 'header' set to 'on',
%                             the column names will be insterted after the custom header.
%   'elecind' - [integer array] Indices of channels to export.
%                             Default is all channels.
%
% Note: for file formats, see readlocs() help  (>> help readlocs)
%
% Author: Arnaud Delorme, Salk Institute, 16 Dec 2002
%
% See also: readlocs(), readelp()

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

function writelocs( chans, filename, varargin );

if nargin < 2
	help writelocs;
	return;
end;

% get infos from readlocs
% -----------------------
%[listtype formatinfo listcolformat formatskip] = readlocs('getinfos');
[chanformat listcolformat] = readlocs('getinfos');
indformat  = [];
for index = 1:length(chanformat),
    if ~isstr(chanformat(index).importformat)
        indformat = [ indformat index ];
    end;
    if isempty(chanformat(index).skipline), chanformat(index).skipline = 0; end;
end;
listtype   = { chanformat(indformat).type };
formatinfo = { chanformat(indformat).importformat };
formatskip = [ chanformat(indformat).skipline ];

g = finputcheck( varargin, ...
                 { 'filetype'	  'string'	 listtype 			'loc';
                   'header'       'string'   { 'on' 'off' } 	'off';
                   'customheader' 'string'   [] 					'';
                   'elecind'      'integer'  [1 Inf]				[];
                   'unicoord'     'string'   { 'on' 'off' } 	'on';
                   'format'		  'cell'	 []					{} }, 'writelocs');
if isstr(g), error(g); end;

if strcmpi(g.unicoord, 'on')
    disp('Uniformizing coordinates');
    chans = convertlocs(chans, 'auto', 'verbose', 'off');
end;

% select channels
% ---------------
if ~isempty(g.elecind)
	chans = chans(g.elecind);
end;

% finding types of input
% ----------------------
if isempty(g.format)
   indexformat = strmatch(lower(g.filetype), listtype, 'exact');
   g.format = formatinfo{indexformat};
   g.skipline = formatskip(indexformat);
else
   g.skipline = 0;
end;

% creating file
% -------------
fid = fopen(filename, 'w');

% exporting header
% ----------------
if ~isempty(g.customheader)
    allstrs = cellstr(g.customheader);
    for index=1:length(allstrs)
        fprintf(fid, '%s\n', allstrs{index});
    end;
end;
if  strcmpi(g.header, 'on') | g.skipline == 2
   for index=1:length(g.format)
      fprintf(fid, '%8s\t', g.format{index});
   end;
   fprintf(fid, '\n');
   for index=1:length(g.format)
      fprintf(fid, '%8s\t', char(ones(1,8)*45));
   end;
   fprintf(fid, '\n');
end;
if g.skipline == 1
   fprintf(fid, '%d\n', length(chans));
end;

% writing infos
% -------------
for indexchan = 1:length(chans)
   for index=1:length(g.format)
      [str, mult] = checkformat(g.format{index});
      if strcmpi(str, 'channum')
         fprintf(fid, '%d', indexchan);
      else
         if ~isfield(chans, str)
            error([ 'Non-existant field: ''' str '''' ]);
         end;
         eval( [ 'chanval = chans(indexchan).' str ';' ] );
         if   isstr(chanval), fprintf(fid, '%8s', chanval);
         else
             if abs(mult*chanval) > 1E-10
                 fprintf(fid, '%8s', num2str(mult*chanval,5));
             else
                 fprintf(fid, '%8s', '0');
             end;
         end;
      end;
      if index ~= length(g.format)
         fprintf(fid, '\t');
      end;
   end;
   fprintf(fid, '\n');
end;
fclose(fid);

return;

% check field format
% ------------------
function [str, mult] = checkformat(str)
	mult = 1;
	if strcmpi(str, 'labels'), str = lower(str); return; end;
	if strcmpi(str, 'channum'), str = lower(str); return; end;
	if strcmpi(str, 'theta'), str = lower(str); return; end;
	if strcmpi(str, 'radius'), str = lower(str); return; end;
	if strcmpi(str, 'sph_theta'), str = lower(str); return; end;
	if strcmpi(str, 'sph_phi'), str = lower(str); return; end;
	if strcmpi(str, 'sph_radius'), str = lower(str); return; end;
	if strcmpi(str, 'sph_theta_besa'), str = lower(str); return; end;
	if strcmpi(str, 'sph_phi_besa'), str = lower(str); return; end;
	if strcmpi(str, 'gain'), str = lower(str); return; end;
	if strcmpi(str, 'calib'), str = lower(str); return; end;
	if strcmpi(str, 'type') , str = lower(str); return; end;
	if strcmpi(str, 'X'), str = upper(str); return; end;
	if strcmpi(str, 'Y'), str = upper(str); return; end;
	if strcmpi(str, 'Z'), str = upper(str); return; end;
	if strcmpi(str, '-X'), str = upper(str(2:end)); mult = -1; return; end;
	if strcmpi(str, '-Y'), str = upper(str(2:end)); mult = -1; return; end;
	if strcmpi(str, '-Z'), str = upper(str(2:end)); mult = -1; return; end;
	if strcmpi(str, 'custum1'), return; end;
	if strcmpi(str, 'custum2'), return; end;
	if strcmpi(str, 'custum3'), return; end;
	if strcmpi(str, 'custum4'), return; end;
   error(['writelocs: undefined field ''' str '''']);

author: PatternRecognition

repo: OpenBMI


text:
function [kip] = ozf2kip(ozf)
% Convert force from ounces force to kip.
% Chad A. Greene 2012
kip = ozf*0.0000625;


text:
function ipFuncCoords = vol2ipXformCoords(gray, inplane, preserveExactValues)
% Return coords from INPLANE functional data that correspond to coords in
% mrVista Gray view
%
%  ipFuncCoords = vol2ipXformCoords(gray, inplane, preserveExactValues)
%
% We first get a 3xn matrix of gray view coordinates, then find the
% corresponding inplane anatomical coordinates, and finally find the
% corresponding inplane functional coordinates. We use this transform when
% we want to convert functional data (e.g., a parameter map, coranal, or
% time series) from the inplane view to the gray view.
%
% INPUTS
%   gray: mrVista view structure (must be a gray view)
%   inplane: mrVista view structure (must be an inplane view)
%   preserveExactValues: boolean. If false, return integer coordinates. If
%                   true, return the calculated (non-integer values). If
%                   non-integer values are returned, then the parent
%                   function will have to deal with these, e.g., via
%                   interpolation.
% OUTPUTS
%   ipFuncCoords: 3xn matrix of coordinates in inplane functional space
%                   corresponding to 3xn matrix of gray view coords
%
% Example:
%   ipFuncCoords = vol2ipXformCoords(gray, inplane)
%
% This code was duplicated in many functions, including ip2volTSeries,
% ip2volParMap, and ip2VolCorAnal. It has now been cropped out of those
% functions and placed in a single routine. This is that routine.
%
% JW 7/2012

% check inputs
if ~exist('preserveExactValues', 'var'), preserveExactValues = false; end

% we need mrSESSION for the alignment matrix
mrGlobals;

% The gray coords are the integer-valued (y,x,z) volume
% coordinates that correspond to the inplanes.  Convert to
% homogeneous form by adding a row of ones.
coords  = viewGet(gray, 'coords');
nVoxels = size(coords, 2);
coords  = double([coords; ones(1,nVoxels)]);

% vol2InplaneXform is the 4x4 homogeneous transform matrix that
% takes volume (y',x',z',1) coordinates into inplane (y,x,z,1)
% coordinates.
vol2InplaneXform = inv(sessionGet(mrSESSION,'alignment'));

% We don't care about the last coordinate in (y,x,z,1), so we
% toss the fourth row of Xform.  Then our outputs will be (y,x,z).
%
vol2InplaneXform = vol2InplaneXform(1:3,:);

% Transform coord positions to the inplanes.  Hence, ipAnatomicalCoords
% contains the inplane position of each of the gray matter voxels.  These
% will generally not fall on integer-valued coordinates, rather they will
% fall between voxels.  Other functions that rely on the output of this
% function will require interpolation to get the data at these
% between-voxel positions.
%
ipAnatomicalCoords = vol2InplaneXform*coords;

% Convert coords from inplane anatomical space to inplane functional space.
% We do this because the anatomical inplane is often higher resolution than
% the functional data.  We preserve the number of coords so that the number
% of output functional voxels is identical to the number of gray or volume
% voxels. If requested, we preserve the exact (non-integer) values, which
% means that the functional coordinates will lie in locations between the
% actual functional data points.
preserveCoords = true;
ipFuncCoords   = ip2functionalCoords(inplane, ipAnatomicalCoords, ...
    [], preserveCoords, preserveExactValues);

% Some confusion about zeros in the output coordinates. It seems best to
% leave the calculations as they are unless there is a problem that someone
% understands. See below.

% HH, 6/16/2012: If we use the code below, we elimintate some voxels
% and then nVoxels is different from the size of ipCoords, leading to a
% mismatch in the size of the expected time series (tSeries) and the
% size of interpolated tSeries, causing an error in  the line,
%  tSeries(frame,:) = interp3(subData, ...
%                                   ipFuncCoords(2,:), ...
%                                   ipFuncCoords(1,:), ...
%                                   ipFuncCoords(3,:), ...
%                                   method);
%
% Hence we comment out the lines below. Leaving in the values with 0
% did not produce an error, at least for the function we tested it on,
% ip2volTSeries.
%
% % ras 12/20/04:
% % occasionally the xformed grayCoords will include a 0
% % coordinate. While it seems this should not happen,
% % I'm applying this band-aid for the time being:
% [badRows badCols] = find(ipFuncCoords==0); %#ok<ASGLU>
% goodCols = setdiff(1:nVoxels, badCols);
% ipFuncCoords = ipFuncCoords(:,goodCols);
% if length(badCols)>1
%     fprintf('%i voxels mapped to slice 0...\n',length(badCols));
% end

return


text:
function [bc] = PoissonIPDGbc3D(ubc)

% function [bc] = PoissonIPDGbc3D(ubc)
% Purpose: Set up the discrete Poisson matrix directly
%          using IP. The operator is set up in the weak form

Globals3D;

% face mass matrix 1
Fm = Fmask(:,1); faceR = r(Fm); faceS = s(Fm);
V2D = Vandermonde2D(N, faceR, faceS);
massEdge(Fm,Fm,1) = inv(V2D*V2D');

% face mass matrix 2
Fm = Fmask(:,2); faceR = r(Fm); faceT = t(Fm);
V2D = Vandermonde2D(N, faceR, faceT);
massEdge(Fm,Fm,2) = inv(V2D*V2D');

% face mass matrix 3
Fm = Fmask(:,3); faceS = s(Fm); faceT = t(Fm);
V2D = Vandermonde2D(N, faceS, faceT);
massEdge(Fm,Fm,3) = inv(V2D*V2D');

% face mass matrix 4
Fm = Fmask(:,4); faceS = s(Fm); faceT = t(Fm);
V2D = Vandermonde2D(N, faceS, faceT);
massEdge(Fm,Fm,4) = inv(V2D*V2D');

% build DG right hand side
bc = zeros(Np, K);

for k1=1:K
  if(~mod(k1,1000)) k1, end;

  for f1=1:Nfaces

    if(EToE(k1,f1)==k1)

      Fm1 = Fmask(:,f1);
      fidM  = (k1-1)*Nfp*Nfaces + (f1-1)*Nfp + (1:Nfp)';

      id = 1+(f1-1)*Nfp + (k1-1)*Nfp*Nfaces;
      lnx = nx(id);  lny = ny(id);  lnz = nz(id);
      lsJ = sJ(id); hinv = Fscale(id);

      Dx = rx(1,k1)*Dr + sx(1,k1)*Ds + tx(1,k1)*Dt;
      Dy = ry(1,k1)*Dr + sy(1,k1)*Ds + ty(1,k1)*Dt;
      Dz = rz(1,k1)*Dr + sz(1,k1)*Ds + tz(1,k1)*Dt;
      Dn1 = lnx*Dx + lny*Dy +lnz*Dz;

      mmE = lsJ*massEdge(:,:,f1);

      gtau = 2*(N+1)*(N+1)*hinv; % set penalty scaling

      bc(:,k1) = bc(:,k1) + (gtau*mmE(:,Fm1) - Dn1'*mmE(:,Fm1))*ubc(fidM);

    end
  end
end
return

repo: Aero-matlab


text:
% Assume a target positioned at x = 1, travelling with speed v = 0.1
state = [1;0.1;2;0;3;0.2;4;0.3];

% Create an instance of a 3D Constant Velocity model
cv = ConstantVelocityX('NumDims',4,'VelocityErrVariance',0.1);

% View the transition matrix and process covariance matrices
F = cv.feval();
Q = cv.covar();

% Predict the target's position and velocity after the interval has passed
newState  = cv.feval(state);

% Do the same as above, but this time add process noise to the prediction
newState2 = cv.feval(state,true);

% Generate 50 random noise samples from the dynamic model
noise = cv.random(50);

% Check how likely the predictions we made are
lik = cv.pdf(newState,state);
lik2 = cv.pdf(newState2,state); % HINT: newState2 should be less likely

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
%DXDEMO Demonstrate distance transform planner using animation
%
% MORPHDEMO(IM, SE, OPTIONS) displays an animation to show the principles
% of the mathematical morphology operations dilation or erosion.  Two
% windows are displayed side by side, input binary image on the left and
% output image on the right.  The structuring element moves over the input
% image and is colored red if the result is zero, else blue.  Pixels in
% the output image are initially all grey but change to black or white
% as the structuring element moves.
%
% OUT = MORPHDEMO(IM, SE, OPTIONS) as above but returns the output image.
%
% Options::
% 'dilate'      Perform morphological dilation
% 'erode'       Perform morphological erosion
% 'delay'       Time between animation frames (default 0.5s)
% 'scale',S     Scale factor for output image (default 64)
% 'movie',M     Write image frames to the folder M
%
% Notes::
% - This is meant for small images, say 10x10 pixels.
%
% See also IMORPH, IDILATE, IERODE.

function out = dxdemo(map, goal, varargin)

    opt.delay = 0;
    opt.movie = [];
    opt.scale = 64;
    opt.metric = {'euclidean', 'manhattan'};
    opt = tb_optparse(opt, varargin);

    set(gcf, 'Position', [90         435        1355         520])
    clf

    goal = [4 8];
    start = [7 2];
    opt.metric = 'euclidean';
    opt.movie = [];%'dxform2.mp4'
    opt.delay = 0;

    % make a simple map
    occgrid = zeros(10,10);
    occgrid(4:6,3:7) = 1;
    %occgrid(7:8,7) = 1;  % extra bit

    cost0 = occgrid;
    cost0(cost0==1) = NaN;

    cost = cost0;
    cost(cost0==0) = Inf;
    cost(goal(2), goal(1)) = 0;

    if ~isempty(opt.movie)
        anim = Animate(opt.movie);
    end

    if ~isempty(opt.movie)
        anim.add();
    end

    switch opt.metric
        case 'cityblock'
            m = [inf 1 inf
                1  0  1
                inf 1 inf];
        case 'euclidean'
            r2 = sqrt(2);
            m = [r2 1 r2
                1 0  1
                r2 1 r2];
        otherwise
            error('unknown distance metric');
    end

    iteration = 0;
    ninf = 0;
    n2 = 1; % half width
    n22 = 1.5;

    newcost = inf(size(cost));

    title('Wavefront path planning simulation');

    while true
        iteration = iteration+1;

        maxval = max(max(cost(isfinite(cost))));

        subplot(121)
        showpixels(cost, 'contrast', maxval*0.7, 'fmt', '%.2g', 'cscale', [0 maxval+2], 'fontsize', 20, 'nancolor', 'nohideinf', 'infsymbol', 'nohidenan', 'infcolor')
        xlabel('x', 'FontSize', 20); ylabel('y', 'FontSize', 20);
            hpatch1 = patch(1, 1,  'y', 'FaceAlpha', 0.5);

        for r=n2+1:numrows(cost)-n2
            for c=n2+1:numcols(cost)-n2

                win = cost(r-n2:r+n2, c-n2:c+n2);

                if isnan(cost(r,c))
                    newcost(r,c) = NaN;
                else
                    newcost(r,c) = min(min(win+m));
                end

                % animate the patch
                hpatch1.XData = [c-n22 c+n22 c+n22 c-n22];
                hpatch1.YData = [r-n22 r-n22 r+n22 r+n22];

                subplot(122)
                cla
                showpixels(newcost, 'contrast', maxval*0.7, 'fmt', '%.2g', 'cscale', [0 maxval+2], 'fontsize', 20, 'nancolor', 'nohideinf', 'infsymbol', 'nohidenan', 'infcolor')
                xlabel('x', 'FontSize', 20); ylabel('y', 'FontSize', 20)
                    hpatch2 = patch([c-0.5 c+0.5 c+0.5 c-0.5], [r-0.5 r-0.5 r+0.5 r+0.5],  'y', 'FaceAlpha', 0.5);

                if ~isempty(opt.movie)
                    anim.add();
                end
                if opt.delay == 0
                    drawnow
                else
                    pause(opt.delay);
                end
            end
        end

        cost = newcost;
        ninfnow = sum(sum( isinf(cost(2:end-1,2:end-1)) )) % current number of Infs
        if ninfnow == 0 || ninfnow == ninf
            % stop if the number of Infs left in the map had stopped reducing
            % it may never get to zero if there are unreachable cells in the map
            break;
        end
        ninf = ninfnow;

    if ~isempty(opt.movie)
        anim.close();
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function x = triangle_cdf_inv ( cdf, a, b, c )

%*****************************************************************************80
%
%% TRIANGLE_CDF_INV inverts the Triangle CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%    0.0 <= CDF <= 1.0.
%
%    Input, real A, B, C, the parameters of the PDF.
%    A <= B <= C and A < C.
%
%    Output, real X, the corresponding argument.
%
  if ( cdf < 0.0 | 1.0 < cdf )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGLE_CDF_INV - Fatal error!\n' );
    fprintf ( 1, '  CDF < 0 or 1 < CDF.\n' );
    error ( 'TRIANGLE_CDF_INV - Fatal error!' );
  end

  d = 2.0 / ( c - a );
  cdf_mid = 0.5 * d * ( b - a );

  if ( cdf <= cdf_mid )
    x = a + sqrt ( cdf * ( b - a ) * ( c - a ) );
  else
    x = c - sqrt ( ( c - b ) * ( ( c - b ) - ( cdf - cdf_mid ) * ( c - a ) ) );
  end

  return
end


text:
function test_bug2685

% MEM 4gb
% WALLTIME 00:10:00
% DEPENDENCY ft_scalpcurrentdensity ft_fetch_sens

%% load data
load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug2685/bug2685.mat'));

%% scalp current density
cfg                 = [];
cfg.method          = 'spline';
cfg.elec            = ERP_standard.elec;
%cfg.elec            = % Do we have electrode positions?

scd   = ft_scalpcurrentdensity(cfg, ERP_standard);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%--------------------------------------------------------------------------
%   [dataout] = read_bin(file_name,file_type,mode_bl)
%--------------------------------------------------------------------------
%   功能:
%   读取bin文件数据
%--------------------------------------------------------------------------
%   输入:
%           file_name                  文件名
%           file_type                  类型,如int8,single(float),double
%           mode_bl                    大小端     'n'系统默认
%                                                 'b' 大端(常用)
%                                                 'l' 小端(常用)
%                                                 's' 大端64bit
%                                                 'a' 小端64bit
%   output:
%           output                     输出读取数据结果
%--------------------------------------------------------------------------
%   例子:
%   read_bin('1.bin','int8','b')
% ans =
%      1
%      2
%      3
%      4
%      5
%--------------------------------------------------------------------------
function [output] = read_bin(file_name,file_type,mode_bl)
if nargin == 2
    mode_bl = 'n';
end
f = fopen(file_name,'r');
output = fread(f,file_type,mode_bl);
fclose(f);

repo: SP


text:
function n = ndims( x )
n = length( x.size_ );

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

author: yu-jiang

repo: radpbook

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