text:
function tapas_hgf_categorical_plotTraj(r)
% Plots trajectories estimated by fitModel for the hgf_categorical perceptual model
% Usage:  est = tapas_fitModel(responses, inputs); tapas_hgf_plotTraj(est);
%
% --------------------------------------------------------------------------------------------------
% 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/>.

% Check whether we have a configuration structure
if ~isfield(r,'c_prc')
    error('tapas:hgf:ConfigRequired', 'Configuration required: before calling tapas_hgf_categorical_plotTraj, tapas_hgf_categorical_config has to be called.');
end

% Number of outcomes
no = r.c_prc.n_outcomes;

% Define colors
colors = [1 0 0; 0.67 0 1; 0 0.67 1; 0.67 1 0];

% Set up display
scrsz = get(0,'screenSize');
outerpos = [0.2*scrsz(3),0.2*scrsz(4),0.8*scrsz(3),0.8*scrsz(4)];
figure(...
    'OuterPosition', outerpos,...
    'Name','HGF trajectories');

% Number of trials
t = size(r.u,1);

% Optional plotting of standard deviations (true or false)
plotsd2 = true;
plotsd3 = true;

% Subplots
subplot(3,1,1);

if plotsd3 == true
    upper3prior = r.p_prc.mu3_0 +sqrt(r.p_prc.sa3_0);
    lower3prior = r.p_prc.mu3_0 -sqrt(r.p_prc.sa3_0);
    upper3 = [upper3prior; r.traj.mu(:,3)+sqrt(r.traj.sa(:,3))];
    lower3 = [lower3prior; r.traj.mu(:,3)-sqrt(r.traj.sa(:,3))];

    plot(0, upper3prior, 'ob', 'LineWidth', 1);
    hold all;
    plot(0, lower3prior, 'ob', 'LineWidth', 1);
    fill([0:t, fliplr(0:t)], [(upper3)', fliplr((lower3)')], ...
         'b', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
end
plot(0:t, [r.p_prc.mu3_0; r.traj.mu(:,3)], 'b', 'LineWidth', 2);
hold all;
plot(0, r.p_prc.mu3_0, 'ob', 'LineWidth', 2); % prior
xlim([0 t]);
title('Posterior expectation \mu_3 of log-volatility of tendency x_3', 'FontWeight', 'bold');
xlabel('Trial number');
ylabel('\mu_3');

subplot(3,1,2);
if plotsd2 == true
    for j=1:no
    upper2prior = r.p_prc.mu2_0(j) +sqrt(r.p_prc.sa2_0(j));
    lower2prior = r.p_prc.mu2_0(j) -sqrt(r.p_prc.sa2_0(j));
    upper2 = [upper2prior; r.traj.mu(:,2,j)+sqrt(r.traj.sa(:,2,j))];
    lower2 = [lower2prior; r.traj.mu(:,2,j)-sqrt(r.traj.sa(:,2,j))];

    plot(0, upper2prior, 'o', 'Color', colors(j,:), 'LineWidth', 1);
    hold all;
    plot(0, lower2prior, 'o', 'Color', colors(j,:), 'LineWidth', 1);
    fill([0:t, fliplr(0:t)], [(upper2)', fliplr((lower2)')], ...
         colors(j,:), 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
    end
end
for j=1:no
    plot(0:t, [r.p_prc.mu2_0(j); r.traj.mu(:,2,j)], 'Color', colors(j,:), 'LineWidth', 2);
    hold all;
    plot(0, r.p_prc.mu2_0(j), 'o', 'Color', colors(j,:), 'LineWidth', 2); % prior
end
xlim([0 t]);
title('Posterior expectations \mu_2 of tendencies x_2', 'FontWeight', 'bold');
xlabel({'Trial number', ' '}); % A hack to get the relative subplot sizes right
ylabel('\mu_2');
hold off;

subplot(3,1,3);
for j=1:no
    plot(0:t, [tapas_sgm(r.p_prc.mu2_0(j), 1); tapas_sgm(r.traj.mu(:,2,j), 1)], 'Color', colors(j,:), 'LineWidth', 2);
    hold all;
    plot(0, tapas_sgm(r.p_prc.mu2_0(j), 1), 'o', 'Color', colors(j,:), 'LineWidth', 2); % prior
end
u = r.u(:,1);
for j=1:no
    plot(find(u==j), -0.08*ones([1 length(find(u==j))]), '.', 'Color', colors(j,:)); % inputs
end
if ~isempty(find(strcmp(fieldnames(r),'y'))) && ~isempty(r.y)
    y = r.y(:,1);
    if ~isempty(find(strcmp(fieldnames(r),'irr')))
        y(r.irr) = NaN; % weed out irregular responses
        plot(r.irr,  1.08.*ones([1 length(r.irr)]), 'x', 'Color', [1 0.7 0], 'Markersize', 11, 'LineWidth', 2); % irregular responses
    end
    for j=1:no
        plot(find(y==j), 1.08*ones([1 length(find(y==j))]), '.', 'Color', colors(j,:)); % responses
    end
    title(['Response y (top dot row), input u (bottom dot row), and posterior probability of outcomes s(\mu_2) for \kappa=', ...
           num2str(r.p_prc.ka), ', \omega=', num2str(r.p_prc.om), ', \vartheta=', num2str(r.p_prc.th)], ...
          'FontWeight', 'bold');
    ylabel('y, u, s(\mu_2)');
    axis([0 t -0.1 1.15]);
else
    title(['Input u (bottom dot row) and posterior probability of outcomes s(\mu_2) for \kappa=', ...
           num2str(r.p_prc.ka), ', \omega=', num2str(r.p_prc.om), ', \vartheta=', num2str(r.p_prc.th)], ...
      'FontWeight', 'bold');
    ylabel('u, s(\mu_2)');
    axis([0 t -0.1 1.1]);
end
plot(1:t, 0.5, 'k');
xlabel('Trial number');
hold off;

author: translationalneuromodeling

repo: tapas


text:
function output = callosqp(interfacedata)

options = interfacedata.options;
model = yalmip2quadprog(interfacedata);

if options.savedebug
    save debugfile model
end

if options.showprogress;showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);end

% Define QP
n_var = length(model.c);
P = model.Q;
q = model.c;
eye_n = speye(n_var);
A = [model.Aeq;model.A; eye_n];
l = full([model.beq; -inf(length(model.b),1); model.lb]);
u = full([model.beq; model.b; model.ub]);

% Define verbose option
options.osqp.verbose = options.verbose;

% Solve with OSQP
OSQPSolver = osqp;
OSQPSolver.setup(P, q, A, l, u, options.osqp);
results = OSQPSolver.solve();

switch results.info.status_val
    case 1
        problem = 0;
    case 2
        problem = 0;
    case -2
        problem = 3;
    case -3
        problem = 1;
    case 3
        problem = 1;
    case -4
        problem = 2;
    case 4
        problem = 2;
    case -5
        problem = 16;
    case -10
        problem = 11;
    otherwise
        problem = -10;
end

% Solver time
solvertime = results.info.run_time;

% Standard interface
Primal      = results.x(:);
Dual        = results.y(1:end-n_var);
infostr     = yalmiperror(problem,interfacedata.solver.tag);
if ~options.savesolverinput
    solverinput = [];
else
    solverinput = model;
end
if ~options.savesolveroutput
    solveroutput = [];
else
    solveroutput = results;
end

% Standard interface
output = createOutputStructure(Primal,Dual,[],problem,infostr,solverinput,solveroutput,solvertime);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function vivmRunDataSet(dataSetName, invariance, experimentNo, ...
                             selectionCriterion, dVal, seedVal);

% VIVMRUNDATASET Try the virtual IVM on a data set and save the results.

% IVM

if nargin < 6
  seedVal = [];
  randn('seed', 1e5)
  rand('seed', 1e5)
else
  randn('seed', seedVal);
  rand('seed', seedVal);
end

jobId = getenv('PBS_JOBID');
jobName = getenv('PBS_JOBNAME');
fprintf('Seed %2.0e\n', seedVal);

if isempty(seedVal)
  [X, y, XTest, yTest] = mapLoadData(dataSetName);
else
  [X, y, XTest, yTest] = mapLoadData(dataSetName, seedVal);
end

capitalName = dataSetName;
capitalName(1) = upper(capitalName(1));

options = ivmOptions;
options.extIters = 0; % don't optimise parameters
mu = zeros(size(yTest));
varSigma = zeros(size(yTest));

load(['dem' capitalName '_' ...
      num2str(experimentNo) '_' ...
      num2str(seedVal) '.mat'])

% Create new training set with virtual SVs.
origIndex = ivmInfoStore.I;
[X, y] = ivmVirtual(X(ivmInfoStore.I, :), y(ivmInfoStore.I, ...
                                                  :), invariance);
tic
model = ivmRun(X, y, kernStore, ...
               noiseStore, selectionCriterion, dVal, ...
               options);

runTime = toc;
if ~isempty(XTest) & ~isempty(yTest);
  yPred = ivmOut(model, XTest);
  testError = 1-sum(yPred==yTest)/size(yTest, 1);
  fprintf('Data set %s, test error %2.4f\n', dataSetName, testError);
end
[kern, noise, ivmInfo] = ivmDeconstruct(model);
save(['dem' capitalName '_' invariance '_' num2str(experimentNo) '_d' num2str(dVal) '_seed' num2str(seedVal)], 'testError', ...
     'ivmInfo', 'kern', 'noise', 'runTime', 'jobId', ...
     'jobName', 'origIndex')

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [u_new, V, exitflag, output] = solveOptimalControlProblem (fst, varargin)

% fst.horizon, fst.xmeasure, fst.u0, fst.opt_option, fst.price, fst.net_load, battery
%SOLVEOPTIMALCONTROLPROBLEM Summary of this function goes here
%   solves the optimal control problem of the

    %x = computeOpenloopSolution(fst); %For linear constraints

    % Set control and linear bounds
    A = [];
    b = [];
    Aeq = [];
    beq = [];
    lb = [];
    ub = [];
    for k=1:fst.horizon  %Aggregation
        [Anew, bnew, Aeqnew, beqnew, lbnew, ubnew] = fst.l_constraints( fst, k );

        A = blkdiag(A,Anew);
        b = [b, bnew];
        Aeq = blkdiag(Aeq,Aeqnew);
        beq = [beq, beqnew];
        lb = [lb, lbnew];
        ub = [ub, ubnew];
    end

    % Solve optimization problem
    [u_new, V, exitflag, output] = fmincon( @(u) fst.costfunction( fst, u ), fst.u0 , ...    % Objective
         A, b, Aeq, beq, lb, ub, ...                                            % Linear Constarints
        @(u) fst.nonlinearconstraints(fst, u ), fst.option);                    % Nonlinear Constraints
end

author: juchengquan


text:
classdef GPSO < ALGORITHM
% <single> <real> <large/none> <constrained/none>
% Gradient based particle swarm optimization algorithm

%------------------------------- Reference --------------------------------
% M. M. Noel, A new gradient based particle swarm optimization algorithm
% for accurate computation of global minimum, Applied Soft Computing, 2012,
% 12: 353-359.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB Platform
% for Evolutionary Multi-Objective Optimization [Educational Forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Generate random population
            Population = Problem.Initialization();
            Pbest      = Population;
            [~,best]   = min(FitnessSingle(Pbest));
            Gbest      = LocalSearch(Problem,Pbest(best));

            %% Optimization
            while Algorithm.NotTerminated(Population)
                Population     = OperatorPSO(Problem,Population,Pbest,Gbest);
                replace        = FitnessSingle(Pbest) > FitnessSingle(Population);
                Pbest(replace) = Population(replace);
                [~,best]       = min(FitnessSingle(Pbest));
                Gbest          = LocalSearch(Problem,Pbest(best));
            end
        end
	end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function c = ttv(a,v,dims)
%TTV Tensor times vector for ktensor.
%
%   Y = TTV(X,A,N) computes the product of Kruskal tensor X with a
%   (column) vector A.  The integer N specifies the dimension in X
%   along which A is multiplied.  If size(A) = [I,1], then X must have
%   size(X,N) = I.  Note that ndims(Y) = ndims(X) - 1 because the N-th
%   dimension is removed.
%
%   Y = TTV(X,{A1,A2,...}) computes the product of tensor X with a
%   sequence of vectors in the cell array.  The products are computed
%   sequentially along all dimensions (or modes) of X. The cell array
%   contains ndims(X) vectors.
%
%   Y = TTV(X,{A1,A2,...},DIMS) computes the sequence tensor-vector
%   products along the dimensions specified by DIMS.
%
%   See also TENSOR/TTV, KTENSOR, KTENSOR/TTM.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

%%%%%%%%%%%%%%%%%%%%%%
%%% ERROR CHECKING %%%
%%%%%%%%%%%%%%%%%%%%%%

% Check the number of arguments
if (nargin < 2)
    error('TTV requires at least two arguments.');
end

% Check for 3rd argument
if ~exist('dims','var')
    dims = [];
end

% Check that 2nd argument is cell array. If not, recall with v as a
% cell array with one element.
if ~iscell(v)
    c = ttv(a,{v},dims);
    return;
end

% Get sorted dims and index for multiplicands
[dims,vidx] = tt_dimscheck(dims,ndims(a),numel(v));

% Check that each multiplicand is the right size.
for i = 1:numel(dims)
    if ~isequal(size(v{vidx(i)}),[size(a,dims(i)) 1])
        error('Multiplicand is wrong size');
    end
end

% Figure out which dimensions will be left when we're done
remdims = setdiff(1:ndims(a),dims);

% Collapse dimensions that are being multiplied out
newlambda = a.lambda;
for i = 1:numel(dims)
    newlambda = newlambda .* ( a.u{dims(i)}' * v{vidx(i)} );
end

% Create final result
if isempty(remdims)
    c = sum(newlambda);
else
    c = ktensor(newlambda,a.u{remdims});
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
trx = Trx('trxfilestr','trx.mat',...
                        'moviefilestr','movie.mov',...
                        'perframedir','perframe',...
                        'default_landmark_params',landmark_params,...
                        'perframe_params',perframe_params);
trx.AddExpDir(expdir);

% clean this data to force computation
if forcecompute,
  %deletefns = setdiff(perframefns,Trx.TrajectoryFieldNames());
  trx.CleanPerFrameData();
end

% compute each
for i = 1:nfns,
  fn = perframefns{i};
  fprintf(logfid,'Computing %s...\n',fn);
  trx.(fn);
end

author: kristinbranson

repo: JAABA


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

%*****************************************************************************80
%
%% LAGUERRE_POLYNOMIAL_VALUES returns some values of the Laguerre polynomial.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      LaguerreL[n,x]
%
%  Differential equation:
%
%    X * Y'' + (1-X) * Y' + N * Y = 0
%
%  First terms:
%
%      1
%     -X   +  1
%   (  X^2 -  4 X   +   2 ) / 2
%   ( -X^3 +  9 X^2 -  18 X   +    6 ) / 6
%   (  X^4 - 16 X^3 +  72 X^2 -   96 X +      24 ) / 24
%   ( -X^5 + 25 X^4 - 200 X^3 +  600 X^2 -   600 X   +   120 ) / 120
%   (  X^6 - 36 X^5 + 450 X^4 - 2400 X^3 +  5400 X^2 -  4320 X   +   720 ) / 720
%   ( -X^7 + 49 X^6 - 882 X^5 + 7350 X^4 - 29400 X^3 + 52920 X^2 - 35280 X + 5040 ) / 5040
%
%  Recursion:
%
%    L(0)(X) = 1,
%    L(1)(X) = 1-X,
%    N * L(N)(X) = (2*N-1-X) * L(N-1)(X) - (N-1) * L(N-2)(X)
%
%  Orthogonality:
%
%    Integral ( 0 <= X < oo ) exp ( - X ) * L(N)(X) * L(M)(X) dX
%    = 0 if N /= M
%    = 1 if N == M
%
%  Special values:
%
%    L(N)(0) = 1.
%
%  Relations:
%
%    L(N)(X) = (-1)^N / N! * exp ( x ) * (d/dx)^n ( exp ( - x ) * X^n )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  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, integer N, the order of the polynomial.
%
%    Output, real X, the point where the polynomial is evaluated.
%
%    Output, real FX, the value of the function.
%
  n_max = 17;

  fx_vec = [ ...
      0.1000000000000000E+01, ...
      0.0000000000000000E+00, ...
     -0.5000000000000000E+00, ...
     -0.6666666666666667E+00, ...
     -0.6250000000000000E+00, ...
     -0.4666666666666667E+00, ...
     -0.2569444444444444E+00, ...
     -0.4047619047619048E-01, ...
      0.1539930555555556E+00, ...
      0.3097442680776014E+00, ...
      0.4189459325396825E+00, ...
      0.4801341790925124E+00, ...
      0.4962122235082305E+00, ...
     -0.4455729166666667E+00, ...
      0.8500000000000000E+00, ...
     -0.3166666666666667E+01, ...
      0.3433333333333333E+02  ];

  n_vec = [ ...
     0,  1,  2, ...
     3,  4,  5, ...
     6,  7,  8, ...
     9, 10, 11, ...
    12,  5,  5, ...
     5,  5 ];

  x_vec = [ ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     0.5E+00, ...
     3.0E+00, ...
     5.0E+00, ...
     1.0E+01 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

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

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function preprocess_shape_collection(path_shapes,path_save,global_params, shot_params)
    d = dir([path_shapes,'*.mat']);

    if ~exist(path_save, 'dir')
       mkdir(path_save)
    end

    for i=1:numel(d)
        %Load shapes
        model = load([path_shapes,d(i).name]);
        model.X = model.VERT(:,1); model.Y = model.VERT(:,2); model.Z = model.VERT(:,3);

        %Calculate LBO eigenfunctions and SHOT descriptors
        [model_evecs,~,model_evals,model_S] = extract_eigen_functions_new(model,global_params.num_evecs);
        model_shot = calc_shot(model.VERT', model.TRIV', 1:model.n, shot_params.num_bins, shot_params.radius, 3)';

        %save as single
        model_evecs_trans = single(model_evecs'*model_S);
        model_evecs = single(model_evecs);
        model_S = single(full(diag(model_S)));
        save([path_save,d(i).name],'model_shot', 'model_evecs', 'model_evecs_trans', 'model_S', 'shot_params');

        display(i)
    end

    delete(gcp('nocreate'))
end

author: OshriHalimi


text:
function outList = transform2switches(stimList)
% assumes 2 stimtypes, 1 and 2
% classifies them into no switch (1 and 2) or switch (3 and 4)
outList(1,1) = stimList(1,1);
for i = 2:size(stimList,1)
   switch stimList(i-1,1)		% based on previous trial
   case 1
      if stimList(i,1) == 2, outList(i,1) = 4;,else outList(i,1) = stimList(i,1);,end
   case 2
      if stimList(i,1) == 1, outList(i,1) = 3;,else outList(i,1) = stimList(i,1);,end
   case 5
      if i > 2 & stimList(i,1) == 1 & stimList(i-2,1) == 2, outList(i,1) = 3;
      elseif i > 2 & stimList(i,1) == 2 & stimList(i-2,1) == 1, outList(i,1) = 4;
      else outList(i,1) = stimList(i,1);
      end
   otherwise outList(i,1) = stimList(i,1);
   end
end
return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [TABLE, beta, tstat] = ARDLprint(ARDL,vnames,approx,vnames_ex)
% =======================================================================
% Prints the output of a VAR estimation
% =======================================================================
% [beta, tstat, TABLE] = VARprint(VAR,VARopt,approx)
% -----------------------------------------------------------------------
% INPUT
%   - VAR: structure output of VARmodel function
%   - VARopt: options of the VAR (see VARopt from VARmodel)
% -----------------------------------------------------------------------
% OPTIONAL INPUT
%   approx: number of decimal digits. Default = 4
%------------------------------------------------------------------------
% OUPUT
%   beta: table of estimated coefficients in cell array
%   tstat: table of t-stats in cell array
%   TABLE: table of estimated coefficients and t-stats in cell array
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2012. Updated November 2020
% -----------------------------------------------------------------------

%% Check inputs
%===============================================
% Check they are not empty
if isempty(vnames)
    error('You need to add label for endogenous variables in VARopt');
end
if ~exist('approx','var')
    approx = 4;
end
if ARDL.nvar_ex>0
    if ~exist('vnames_ex','var')
        error('You need to add label for exogenous variables');
    end
end

%% Retrieve and initialize variables
%=============================================================
nlag = ARDL.nlag;
nlag_ex = ARDL.nlag_ex;
const = ARDL.const;
% vnames_ex = vnames(2:end);
% vnames = vnames(1);

%% Additional checks
%=============================================================
% Check size of vnames (and change it if necessary)
htext = vnames;
nvars = 1;
nvars_ex = 1;

%% Labels of deterministic components
%===============================================
switch const
    case 0
        aux = [];
    case 1
        aux = {'c'};
    case 2
        aux = {'c';'trend'};
    case 3
        aux = {'c';'trend';'trend2'};
end
vtext = {' '};
vtext = [vtext; aux];
clear aux

%% Labels of lagged variables
%===============================================
for jj=1:nlag
    for ii=1:nvars
        aux(ii,1) = {[vnames{ii} '(-' num2str(jj) ')' ]};
    end
    vtext = [vtext ; aux];
end
clear aux

%% Labels of exogenous variables
%===============================================
vtext = [vtext ; vnames_ex'];
if nlag_ex > 0
    for jj=1:nlag_ex
        for ii=1:nvars_ex
            aux(ii,1) = {[vnames_ex{ii} '(-' num2str(jj) ')' ]};
        end
        vtext = [vtext ; aux];
    end
    clear aux
end

%% Save
%===============================================
% Save a beta table
beta = roundnum2cell(ARDL.beta,approx);
beta = [htext; beta];
beta = [vtext beta];

% Save a std err table
bstd = roundnum2cell(ARDL.bstd,approx);
bstd = [htext; bstd];
bstd = [vtext bstd];

% Save a tstat table
tstat = ARDL.tstat;
tstat = roundnum2cell(tstat,2);
tstat = [htext; tstat];
tstat = [vtext tstat];

% Save a p-value table
tprob = ARDL.tprob;
tprob = roundnum2cell(tprob,approx);
tprob = [htext; tprob];
tprob = [vtext tprob];

% Save a beta & tstat table
nn = size(beta,1)-1;
TABLE = {''};
index = 1;
for ii=1:nn
    for jj=1:nvars
        TABLE(index,jj) = beta(1+ii,1+jj);
        TABLE(index+1,jj) = bstd(1+ii,1+jj);
        aux1 = cell2mat(tstat(1+ii,1+jj)); % get the numeric value from cell
        aux2 = [ '[' num2str(aux1) ']' ];  % add parenthesis to t-stat value
        TABLE{index+2,jj} = aux2;
        TABLE(index+3,jj) = tprob(1+ii,1+jj);
    end
    index = index+4;
end
clear aux
TABLE = [htext; TABLE];

% Create vertical label
TAB_v = {''};
index = 2;
for ii=1:nn
    TAB_v(index,1) = vtext(1+ii);
    TAB_v(index+1,1) = {['std(' vtext{1+ii} ')']};
    TAB_v(index+2,1) = {['t(' vtext{1+ii} ')']};
    TAB_v(index+3,1) = {['p(' vtext{1+ii} ')']};
    index = index + 4;
end
TABLE = [TAB_v TABLE];
TABLE(index,2) = num2cell(ARDL.rsqr); TABLE(index,1) = {'R2'};    index = index+1;
TABLE(index,2) = num2cell(ARDL.rbar); TABLE(index,1) = {'R2bar'}; index = index+1;
if const>0
    TABLE(index,2) = num2cell(ARDL.F); TABLE(index,1) = {'F'};     index = index+1;
else
    TABLE(index,2) = [];               TABLE(index,1) = {'F'};     index = index+1;
end
TABLE(index,2) = num2cell(ARDL.nobs);  TABLE(index,1) = {'Obs'};   index = index+1;
%
% %% Print the table on screen (only beta)
% %===============================================
% info.cnames = char(htext);
% info.rnames = char(vtext);
% disp(' ')
% %disp('---------------------------------------------------------------------')
% disp(' ')
% disp('Reduced form VAR estimation:')
% disp(' ')
% mprint(ARDL.Ft,info)
% %disp('---------------------------------------------------------------------')

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function varargout = contour3( f, varargin )
%CONTOUR3   3-D contour plot of a SPHEREFUN.
%   CONTOUR3(F) is a contour plot of F treating the values of F as heights
%   above the sphere, analagous to SURF(F, 'projection', 'bumpy'). A
%   contour plot shows the level curves of F for some values V. The values
%   V are chosen automatically.
%
%   CONTOUR3(F, N) draws N contour lines, overriding the automatic number.
%   The values V are still chosen automatically.
%
%   CONTOUR3(F, V) draws LENGTH(V) contour lines at the values specified in
%   the vector V. Use CONTOUR3(F, [V V]) to compute a single contour at the
%   level V.
%
%   CONTOUR3(F, 'NUMPTS', N) plots the contour lines on an N by N uniform
%   grid. If NUMPTS is not given then we plot on a 200 by 200 grid.
%
% See also CONTOUR, CONTOURF.

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

if ( isempty(f) )  % Empty check.
    contour3([]);
    return
end

holdState = ishold;

% Minimum number of plotting points:
minplotnum = 200;

% Extract from the inputs the user defined options:
j = 1;
argin = {};
while ( ~isempty(varargin) )
    if ( strcmpi(varargin{1}, 'numpts') ) % If given numpts then use them.
        minplotnum = varargin{2};
        varargin(1:2) = [];
    else
        argin{j} = varargin{1};
        varargin(1) = [];
        j = j+1;
    end
end

if ( isa(f, 'spherefun') )
    dom = f.domain;
    % Evaluate at equally spaced grid:
    x = linspace(dom(1), dom(2), minplotnum);
    y = linspace(dom(3), dom(4), minplotnum);
    vals = sample(f, minplotnum-1, minplotnum);
    vals = [vals vals(:,1)];
else
    error('CHEBFUN:SPHEREFUN:contour3:inputs', ...
        'Input must be a spherefun.');
end

if ( iscolat(f) )
    xc = @(ll,tt,rr) rr.*cos(ll).*sin(tt);
    yc = @(ll,tt,rr) rr.*sin(ll).*sin(tt);
    zc = @(tt,rr) rr.*cos(tt);
else
    xc = @(ll,tt,rr) rr.*cos(ll).*cos(tt);
    yc = @(ll,tt,rr) rr.*sin(ll).*cos(tt);
    zc = @(tt,rr) rr.*sin(tt);
end

% Use contour rather than contourc so that it can handle parsing the inputs
% correctly.
[c, h] = contour( x', y', vals, argin{:} );

% Extract out the options we need to plot the contours with plot3.
LW = 'LineWidth';
lw = h.LineWidth;
LS = 'LineStyle';
ls = h.LineStyle;
LC = 'Color';
lc = h.LineColor;
levelList = h.LevelList;
clrmap = parula(numel(levelList));

% Remove the contour plot that was generated.
delete(h);

scl = 0.15;           % Match the bumpy parameter in SPHEREFUN/SURF().
lim = [-1-scl 1+scl]; % Pad the axis limits.
m = minandmax2(f);

% If the plot is not being added to another, then set the axis properties.
if ( ~holdState )
    xlim(lim), ylim(lim), zlim(lim)
    daspect([1 1 1])
    grid on, box off
    hold on
end
view(3)

cl = size(c,2);
k = 1;
while ( k < cl )
    kl = c(2,k);
    v = k+1:k+kl;

    % Bump out the radial components according to the function values.
    rr = rescale(c(1,k), 1-scl, 1+scl, 'InputMin', m(1), 'InputMax', m(2));
    xv = xc(c(1,v), c(2,v), rr);
    yv = yc(c(1,v), c(2,v), rr);
    zv = zc(c(2,v), rr);

    % If the line color is a float then we are plotting all contours in a
    % single color.
    if ( isfloat(lc) )
        plot3(xv, yv, zv, LW, lw, LC, lc, LS, ls);
    else
        % We need to plot each contour in a color using the default
        % colormap.
        % Determine the color for the level being plotted.
        clr = clrmap(abs(c(1, k) - levelList) < 10*eps, :);
        plot3(xv, yv, zv, LW, lw, LC, clr, LS, ls);
    end

    k = k+kl+1;
end

if ( ~holdState )
    hold off
end

% Return plot handle if appropriate.
if ( nargout >= 1 )
    warning('CHEBFUN:SPHEREFUN:contour3:outputs', ...
        'Outputs from contour3 are not supported');
    varargout = { [], [] };
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function printCell(fid, cellArray, prefix)
  fprintf(fid, '%s', prefix);
  for ii=1:length(cellArray)
    if ii<length(cellArray)
      fprintf(fid, '%s ', cellArray{ii});
    else
      fprintf(fid, '%s\n', cellArray{ii});
    end
  end
end

author: lmthang


text:
function WG = slnbreconweights(X0, X, G, varargin)
%SLNBRECONWEIGHTS Solve the optimal reconstruction weights on given neighbors
%
% $ Syntax $
%   - WG = slnbreconweights(X0, X, G, ...)
%
% $ Arguments $
%   - X0:       The reference samples to reconstruct the query samples
%   - X:        The query samples
%   - G:        The graph giving the neighborhood relations
%   - WG:       The weighted graph giving the solved weights
%
% $ Description $
%   - WG = slnbreconweights(X0, X, G, ...) solves the optimal weights to
%     reconstruct the samples in X from those in X0. If X is empty, then
%     it would use X0 as X. The graph G indicates the neighborhood
%     relation, having n0 sources and n targets. The WG is a graph with
%     of the same size as G, and the reconstuction weights are placed
%     in the positions corresponding to those in G.
%     You can specify the following properties to control the solving:
%     \*
%     \t    Table. The Properties of Reconstruction Weights Solving
%     \h       name        &     description
%            'constraint'  & The constraint on the solution, it can be
%                            one of the following string to indicate a
%                            single constraint or a cell array of multiple
%                            strings to indicate compound constaints.
%                            - 'nonneg':  non-negative
%                            - 's1':      the weights sum to 1
%                            (default = 's1')
%            'delta'       & The value of regularization. In practice,
%                            regularization is essential to guarantee the
%                            stability of the solution. In implementation,
%                            the diagonal elements of the gram matrix will
%                            be added with a value:
%                               (delta^2) * trace(G) / K
%                            here G is X^T * X, K is the neighbor number.
%                            (default = 0.1)
%            'solver'      & The solver offered by user (function handle).
%                            If the user specify a non-empty solver, then
%                            it will use the user's solver to solve
%                            weights. The solver is like the form:
%                               w = f(X, y)
%                            Here X is d x K neighbor sample matrix, y is
%                            a d x 1 vector representing the target sample.
%                            It should output a K x 1 vector giving the
%                            reconstruction weights.
%                            By default, solver = [], indicating to use
%                            internal solver based on constraint and delta.
%            'thres'       & The thres, if the ratio of a weight value
%                            to the average weight for that reconstruction
%                            is lower than thres, the weight is set
%                            to strictly zeros. This would significantly
%                            reduces the near-zero weights, and thus
%                            reduces the complexity of the graph.
%                            (default = 1e-8)
%     \*
%
% $ Remarks $
%   - When the user specify a non-empty solver, the internal solver will
%     not be used, thus constraint and delta will not take effect.
%
%   - G would be in all acceptable graph form. WG will always be a
%     numeric matrix. If G is a sparse adjmat, then WG would be sparse,
%     otherwise WG is full.
%
%   - With the WG solved, to reconstruct by weighed combination of
%     neighbors, you can simply write it as: Xr = X0 * WG, then Xr
%     is a d x n matrix with the j-th column reconstructed from the
%     referenced samples in Xr using the j-th column's weights in WG.
%
% $ History $
%   - Created by Dahua Lin, on Sep 11st, 2006
%

%% parse and verify input arguments

if nargin < 3
    raise_lackinput('slnbreconweights', 3);
end

if ~isnumeric(X0) || ndims(X0) ~= 2
    error('sltoolbox:invalidarg', ...
        'X should be a 2D numeric matrix');
end
[d, n0] = size(X0);

if isempty(X)
    X = X0;
else
    if ~isnumeric(X) || ndims(X) ~= 2
        error('sltoolbox:invalidarg', ...
            'X0 should be a 2D numeric matrix');
    end
    if size(X, 1) ~= d
        error('sltoolbox:sizmismatch', ...
            'The sample dimension in X is not the same as that in X0');
    end
end
n = size(X, 2);

gi = slgraphinfo(G);
if gi.n ~= n0 || gi.nt ~= n
    error('sltoolbox:sizmismatch', ...
        'The size of the graph is not consisitent with the sample set');
end

opts.constraint = 's1';
opts.delta = 0.1;
opts.solver = [];
opts.thres = 1e-8;
opts = slparseprops(opts, varargin{:});

thres = opts.thres;

%% Prepare parameters

% prepare graph
if ~strcmp(gi.form, 'adjmat')
    G = sladjmat(G, 'sparse', true, 'valtype', 'logical');
end

% prepare solver
if isempty(opts.solver)

    % parse constraints
    cs = opts.constraint;
    if ~isempty(cs)
        if ~iscell(cs)
            cs = {cs};
        end
        constraint = parse_constraints(cs);
    else
        constraint = parse_constraints({});
    end

    % decide solver
    if ~constraint.nonneg
        delta2 = opts.delta^2;
        if ~constraint.s1
            wsolver = @(X, y) internal_wsolver_unc(X, y, delta2);
        else
            wsolver = @(X, y) internal_wsolver_s1(X, y, delta2);
        end
    else
        optimopts = optimset('Display', 'off', 'LargeScale', 'off');
        if ~constraint.s1
            wsolver = @(X, y) internal_wsolver_nonneg(X, y, opts.delta, optimopts);
        else
            wsolver = @(X, y) internal_wsolver_nonneg_s1(X, y, opts.delta, optimopts);
        end
    end

else
    if ~isa(opts.solver, 'function_handle')
        error('The weight solver should be a function handle');
    end
    wsolver = opts.solver;
end

%% main skeleton

% init WG
if issparse(G)
    WG = spalloc(n0, n, nnz(G));
else
    WG = zeros(n0, n);
end

% solve weights
for i = 1 : n
    nbinds = find(G(:,i));
    if ~isempty(nbinds)
        Xnb = X0(:, nbinds);
        y = X(:,i);
        w = wsolver(Xnb, y);
        if thres > 0
            absw = abs(w);
            curthres = thres * sum(absw) / length(w);
            w(absw < curthres) = 0;
        end
        WG(nbinds, i) = w;
    end
end

%% constraint parsing function

function c = parse_constraints(cs)

ncs = length(cs);

c = struct('nonneg', false, 's1', false);

for i = 1 : ncs
    cname = cs{i};
    if ~ischar(cname)
        error('sltoolbox:invalidarg', ...
            'The constraint should be given in char string');
    end
    switch cname
        case 'nonneg'
            c.nonneg = true;
        case 's1'
            c.s1 = true;
        otherwise
            error('sltoolbox:invalidarg', ...
                'Invalid constraint name for weight solving: %s', cname);
    end
end

%% The internal weight solvers

% unconstrained solver
function w = internal_wsolver_unc(X, y, delta2)

[G, Xty] = compute_G_Xty(X, y, delta2);
w = G \ Xty;

% solver with s1 constraint
function w = internal_wsolver_s1(X, y, delta2)

[G, Xty, K] = compute_G_Xty(X, y, delta2);
wu = G \ [Xty, ones(K, 1)];

w = wu(:,1);
u = wu(:,2);
lambda = (1 - sum(w)) / sum(u);
w = w + lambda * u;

% solver with nonnegative constraint
function w = internal_wsolver_nonneg(X, y, delta, optimopts)

[Xa, ya, K] = augformulate(X, y, delta);
if K <= 20
    w = lsqnonneg(Xa, ya, [], optimopts);
else
    lb = zeros(K, 1);
    w = lsqlin(Xa, ya, [], [], [], [], lb, [], [], optimopts);
end

% solver with nonnegative and s1 constraint
function w = internal_wsolver_nonneg_s1(X, y, delta, optimopts)

[Xa, ya, K] = augformulate(X, y, delta);

Aeq = ones(1, K);
beq = 1;
lb = zeros(K, 1);
w = lsqlin(Xa, ya, [], [], Aeq, beq, lb, [], [], optimopts);

% solver preparation function

function [G, Xty, K] = compute_G_Xty(X, y, delta2)

% compute Xt, G, and Xty
K = size(X, 2);
Xt = X';
G = Xt * X;
Xty = Xt * y;

% regularize
if delta2 > 0
    diaginds = (1:K)*(K+1) - K;
    rv = delta2 * sum(G(diaginds)) / K;
    G(diaginds) = G(diaginds) + rv;
end

function [Xa, ya, K] = augformulate(X, y, delta)

K = size(X, 2);
if delta ~= 0
    Xa = [X; delta * eye(K)];
    ya = [y; zeros(K, 1)];
else
    Xa = X;
    ya = y;
end


text:
function ExportAedat2Frames(aedat)

%{
This function exports data to a .aedat file.
The .aedat file format is documented here:

http://inilabs.com/support/software/fileformat/
%}

dbstop if error

if ~exist('aedat', 'var')
	error('Missing input')
end

if ~isfield(aedat, 'data')
	disp('No data to export')
    return
end

% Create the file
if ~isfield(aedat.exportParams, 'filePath')
    error('Missing file path and name')
end

f = fopen(aedat.exportParams.filePath, 'w', 'b');

% Simple - events only - assume DAVIS

% CRLF \r\n is needed to not break header parsing in jAER
fprintf(f,'#!AER-DAT2.0\r\n');
fprintf(f,'# This is a raw AE data file created by an export function in the AedatTools library\r\n');
fprintf(f,'# Data format is int32 address, int32 timestamp (8 bytes total), repeated for each event\r\n');
fprintf(f,'# Timestamps tick is 1 us\r\n');
% Put the source in - use an override if it has been given
if isfield(aedat.exportParams, 'source')
    fprintf(f,['# AEChip: ' aedat.exportParams.source '\r\n']);
else
    fprintf(f,['# AEChip: ' aedat.info.source '\r\n']);
end
fprintf(f,'# End of ASCII Header\r\n');

% DAVIS
% In the 32-bit address (1-based):
% Bit-32   Bit-11   Meaning
% 1        0        APS sample
% 1        1        IMU sample
% 0        1        Special event
% 0        0        Polarity event

yShiftBits = 22;
xShiftBits = 12;
% frameShiftBits = 0;
frameFlagShiftBits = 31;
signalShiftBits = 9;

frameData = aedat.data.frame;

numFrames = frameData.numEvents;
xDim = aedat.info.deviceAddressSpace(1);
yDim = aedat.info.deviceAddressSpace(2);
numPixels = xDim * yDim;

% Allocate horizontal vectors to hold output data.
% Why are the vectors that big? The factor of 2 is because
% we insert dummy 'reset' frames prior to each frame.
samples = uint32(zeros(1, 2 * numFrames * numPixels));
timeStamps = uint32(zeros(1, 2 * numFrames * numPixels));

% The output vector is twice as big again because samples and timeStamps
% will be interspersed in the 'output' vector.
output = uint32(zeros(1, 2 * 2 * numFrames * numPixels));
y = repmat(uint32(yDim - 1 : -1 : 0), 1, xDim * numFrames * 2);
x = repmat(uint32(xDim - 1 : -1 : 0), yDim, numFrames * 2);
x = x(:);
x = x';
% in bit 11 (1-based) 1 means signal read and 0 means reset read.
signalFlag = repmat([zeros(1, numPixels, 'uint32')  ones(1, numPixels, 'uint32') * 2^10], 1, numFrames);
% The last event mask is synonymous with the sample from x=0 y=0; data is
% therefore ordered backwards.
for frameIndex = 1 : numFrames
    samplesTemp = frameData.samples{frameIndex}(:);
    samplesTemp = samplesTemp(end : -1 : 1);
    samples((frameIndex * 2 - 1) * numPixels + 1 : frameIndex * 2 * numPixels) ...
        = samplesTemp ;
    timeStamps((frameIndex - 1) * 2 * numPixels + 1 : frameIndex * 2 * numPixels) ...
        = frameData.timeStampStart(frameIndex);
end
frameFlag = uint32(ones(1, numFrames * 2 * numPixels) * 2 ^ frameFlagShiftBits);
y = y * uint32(2 ^ yShiftBits);
x = x * uint32(2 ^ xShiftBits);
% samples should now be in the range 0-1023 (10-bit).
% subtract samples from 1023. This has the effect of leaving all the reset
% frame samples at 1023 - the highest value, against which the signal frames
% will later be subtracted.
samples = 1023 - samples;

output(1:2:end) = frameFlag + y + x + signalFlag + samples;
output(2:2:end) = timeStamps; % set even elements to timestamps

% write addresses and timestamps
count=fwrite(f, output, 'uint32', 0, 'b')/2; % write 4 byte data
fclose(f);
fprintf('wrote %d events to %s\n', count, aedat.exportParams.filePath);

author: panpanfei


text:
function [sys,x0,str,ts]=Nonlinear_PD_fhan(t,x,u,flag,r,h,h1)

if flag==0

    sys = [0;1;1;2;0;1;1];
     x0 = [0];
    str = [ ];
     ts = [h 0];

elseif flag==2

    sys=-fhan(u(1),u(2),r,h1);

elseif flag==3
    sys=x;

elseif flag==4
    sys=sys+h;

else
    sys=[];
end
function y=fhan(x1,x2,r,h)
d=r*h;
d0=h*d;
y=x1+h*x2;
a0=sqrt(d^2+8*r*abs(y));
if abs(y)>d0
    a=x2+(a0-d)*sign(y)/2;
else
    a=x2+y/h;
end

if abs(a)>d
    y=-r*sign(a);
else
    y=-r*a/d;
end

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
function [ fea, out ] = ex_swirl_flow2( varargin )
%EX_SWIRL_FLOW2 2D Axisymmetric swirl flow in step domain.
%
%   [ FEA, OUT ] = EX_SWIRL_FLOW2( VARARGIN ) Axisymmetric swirl for in tubular step
%   region where the inner cylindrical wall is rotating.
%
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       omega       scalar {100}           Angular rotational velocity (of inner wall)
%       sf_u        string {sflag1}        Shape function for velocity
%       sf_p        string {sflag1}        Shape function for pressure
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { 'omega',    100;
            'sf_u',     'sflag1';
            'sf_p',     'sflag1';
            'iphys',    1;
            'iplot',    1;
            'tol',      [];
            'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Geometry and grid generation.
fea.sdim = {'r' 'z'};

fea.geom.objects = { gobj_rectangle(0.5,1.5,0,3) gobj_rectangle(1.0,1.5,1.5,3,'R2') };
fea = geom_apply_formula( fea, 'R1-R2' );

fea.grid = gridgen( fea, 'hmax', 0.1, 'fid', fid );

% Equation definition.
if ( opt.iphys==1 )

  fea = addphys(fea,@swirlflow);
  fea.phys.sw.eqn.coef{1,end} = { 1 };
  fea.phys.sw.eqn.coef{2,end} = { 1 };
  fea.phys.sw.sfun            = [ repmat( {opt.sf_u}, 1, 3 ) {opt.sf_p} ];
  fea.phys.sw.bdr.sel = [1 1 5 2 1 1];
  fea.phys.sw.bdr.coef{2,end}{2,4} = opt.omega;
  fea.phys.sw.prop.artstab.ps = isequal(opt.sf_u,opt.sf_p);

  fea = parsephys(fea);
  if( isfield(fea,'constr') )
    fea = rmfield(fea,'constr');
  end

else

  opt.sf_u = 'sflag2';
  opt.sf_p = 'sflag1';

  fea.dvar = { 'u', 'v', 'w', 'p' };

  fea.sfun = [ repmat( {opt.sf_u}, 1, 3 ) {opt.sf_p} ];
  c_eqn    = { 'r*rho*u'' - r*miu*(2*ur_r + uz_z  +   wr_z) + r*rho*(u*ur_t + w*uz_t) + r*p_r     = r*Fr - 2*miu/r*u_t + p_t + rho*v*v_t';
               'r*rho*v'' - r*miu*(  vr_r + vz_z) + miu*v_r + r*rho*(u*vr_t + w*vz_t) + rho*u*v_t = r*Fth + miu*(v_r - 1/r*v_t)';
               'r*rho*w'' - r*miu*(  wr_r + uz_r  + 2*wz_z) + r*rho*(u*wr_t + w*wz_t) + r*p_z     = r*Fz';
               'r*ur_t + r*wz_t + u_t = 0' };
  fea.eqn = parseeqn( c_eqn, fea.dvar, fea.sdim );

  fea.coef = { 'rho', 1 ;
               'miu', 1 ;
               'Fr',  0 ;
               'Fth', 0 ;
               'Fz',  0 };

  % Boundary conditions.
  fea.bdr.d = { 0  0  [] 0  0  0 ;
                0  0  [] opt.omega 0 0 ;
                0  0  0  0  0  0 ;
                [] [] [] [] [] [] };
  fea.bdr.n = cell(size(fea.bdr.d));

end

% Fix pressure at p([r,z]=[ro,h/2]) = 0.
[~,ix_p] = min( sqrt( (fea.grid.p(1,:)-1.5).^2 + (fea.grid.p(2,:)-1.5).^2) );
fea.pnt = struct( 'type',  'constr', ...
                  'index', ix_p, ...
                  'dvar',  'p', ...
                  'expr',  '0' );

% Parse and solve problem.
fea = parseprob( fea );
fea.sol.u = solvestat( fea, 'maxnit', 50, 'fid', fid );

% Postprocessing.
if( opt.iplot )
  postplot( fea, 'surfexpr', 'sqrt(u^2+v^2+w^2)', 'isoexpr', 'v' )
end

% Error checking.
out.ref  = [ -6.1 10.5 73 1.25 ];
if( ~got.tol )
  if( opt.sf_u(end) == '2' )
    opt.tol = 0.05;
  else
    opt.tol = 0.3;
  end
end
[u_min,u_max] = minmaxsubd( 'u', fea );
out.val  = [ u_min u_max intsubd('v',fea) intsubd('w',fea) ];
out.pass = mean(abs(out.val-out.ref)./abs(out.ref)) < opt.tol;

if( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
classdef LogarithmicSequence < IncrementalSequence

    properties (Access = private)
       coef
    end

    methods (Access = public)

        function obj = LogarithmicSequence(a,nSteps,initialValue,finalValue)
            obj.coef = a;
            obj.nSteps = nSteps;
            obj.initialValue = initialValue;
            obj.finalValue = finalValue;
            obj.generateAlphaSequence();
        end

    end

    methods (Access = protected)

        function generateAlphaSequence(obj)
            %obj.alpha = logspace(obj.x0,obj.x1,obj.nSteps);
            %              exp = linspace(obj.x0,obj.x1,obj.nSteps);
            %              obj.alpha = 10.^exp;
            %              alpha0 = min(abs(obj.alpha));
            %              alpha1 = max(abs(obj.alpha));
            %              obj.alpha = (obj.alpha - alpha0)/(alpha1-alpha0);
            x = linspace(0,1,obj.nSteps);
            b = 1;
            a = b*(obj.coef+1)/(1-obj.coef);
            obj.alpha = (a*x./((a-1)*x +1)).^b;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function saved=loadEEG(dd, order, band, ival, channel_layout)

[cnt,mrk,mnt]=eegfile_loadMatlab(dd);
cnt = proc_filtButter(cnt,order,band);
cnt = proc_selectChannels(cnt,channel_layout);
% cnt to eo
epo=cntToEpo(cnt,mrk,ival);
%% variables can b changed by the data
classes=size(epo.className,2);
trial=size(epo.x,3)/2/(classes-1);

%% extract the rest class
for ii =1:classes
    if strcmp(epo.className{ii},'Rest')
        epoRest=proc_selectClasses(epo,{epo.className{ii}});
        epoRest.x=datasample(epoRest.x,trial,3,'Replace',false);
        epoRest.y=datasample(epoRest.y,trial,2,'Replace',false);
    else
        epo_check(ii)=proc_selectClasses(epo,{epo.className{ii}});
        % random sampling
        epo_check(ii).x=datasample(epo_check(ii).x,trial,3,'Replace',false);
        epo_check(ii).y=datasample(epo_check(ii).y,trial,2,'Replace',false);
    end
end
if classes<7
    epo_check(size(epo_check,2)+1)=epoRest;
end

%% concatenate the classes
for ii=1:size(epo_check,2)
    if ii==1
        concatEpo=epo_check(ii);
    else
        concatEpo=proc_appendEpochs(concatEpo,epo_check(ii));
    end
end
saved=concatEpo;

author: PatternRecognition

repo: OpenBMI


text:
function dist=mrManDist_new(nodeList,edgeList,startNode,dimDist,noval,radius)
% function [dist,pathIndices]=mrManDist_new(nodeList,edgeList,startNode,dimdist,-1,0)
%
% Obsolete, I think (BW)
%
%% ARW 031201
% This is a wrapper to replace mrManDist with the dijkstra routine from Stanford University
% mrManDist seems to be broken in some way. It certainly doesn't work with small floats
% This is the sort of thing you had to do to call mrManDist....
% [nodeList, edgeList]=generateManDistNodes(mesh);
% mesh.dist = mrManDist(nodeList,edgeList,startNode,dimdist,-1,0);
% dimdist is a scale factor for each linear dimension
% -1 is what's returned for a node with no connection to the start node.
% 0 says we're not sending in a list of node distances
% Edgelist is a 2xn list of node pairs
% Nodelist is a list of 3d node positions.

% To call dijkstra we need to make a weighted sparse connection matrix
% With each entry being the distance between the relevent nodes

disp('mrManDist_new is Obsolete. If you think it is necessary, please rename it and fix mrManDist.')
evalin('caller','mfilename')

return;

mesh.connectionMatrix = buildConnectionMatrix(nodeList, edgeList);

% Get rid of all the extra rubbish
mesh.uniqueVertices = nodeList(1:3,:)';

D = find3DNeighbourDists(mesh,dimDist);
dist = dijkstra(D,startNode);
dist = sqrt(dist);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function mtx = CreateMatrixFromPrimaries(R, G, B, Wp)
%
%       mtx = CreateMatrixFromPrimaries(R, G, B, Wp)
%
%
%        Input:
%           -R: the red primary for the given color space expressed as an
%           XYZ color.
%           -G: the green primary for the given color space expressed as an
%           XYZ color.
%           -B: the blue primary for the given color space expressed as an
%           XYZ color.
%           -Wp: the white-point primary for the given color space expressed as an
%           XYZ color.
%
%        Output:
%           -mtx: a conversion matrix from XYZ color space to the color
%           space defined by the three input primaries (R, G, and B) and
%           white point (Wp).
%
%     Copyright (C) 2018  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

A_R = [R(1) R(2) R(3) 0 0 0 0 0 0; ...
       0 0 0 R(1) R(2) R(3) 0 0 0; ...
       0 0 0 0 0 0 R(1) R(2) R(3)];
b_R = [1; 0; 0];

A_G = [G(1) G(2) G(3) 0 0 0 0 0 0; ...
       0 0 0 G(1) G(2) G(3) 0 0 0; ...
       0 0 0 0 0 0 G(1) G(2) G(3)];
b_G = [0; 1; 0];

A_B = [B(1) B(2) B(3) 0 0 0 0 0 0; ...
       0 0 0 B(1) B(2) B(3) 0 0 0; ...
       0 0 0 0 0 0 B(1) B(2) B(3)];
b_B = [0; 0; 1];

A_Wp = [Wp(1) Wp(2) Wp(3) 0 0 0 0 0 0; ...
       0 0 0 Wp(1) Wp(2) Wp(3) 0 0 0; ...
       0 0 0 0 0 0 Wp(1) Wp(2) Wp(3)];

b_Wp = [1; 1; 1];

A = [A_R; A_G; A_B; A_Wp];
b = [b_R; b_G; b_B;  b_Wp];

mtx = A \ b;

mtx = reshape(mtx, 3, 3)';

end

author: banterle


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

%%
clear; close all; clc;

%% Syntax
% |[M]=rigidTransformationMatrixDirect(V1,V2);|

%% Description
% This function computes the rigid transformation matrix (translation, and
% rotation). For the two point matches input sets V1 and V2.

%% Examples

%%
% Plot settings
fontSize=15;
faceAlpha=1;
edgeColor='k';

%% Example:

%% Determine the rigid transformation to overlay two surfaces
% Below is a demonstration of the features of the |rigidTransformationMatrixDirect| function

% Load example patch data
[F,Vd]=parasaurolophus;

%Translation
OR=[-1 2 -3]; %Translations
T  = [1 0 0 OR(1);...
      0 1 0 OR(2);...
      0 0 1 OR(3);...
      0 0 0 1];
%Rotation
a=[-0.25*pi 0.75*pi 0.1*pi]; %Euler angles
R  = eye(4,4); R(1:3,1:3)=euler2DCM(a);

%Scaling
v=[2 3 1]; %Scaling factors
S  = [v(1) 0    0    0;...
      0    v(2) 0    0;...
      0    0    v(3) 0;...
      0    0    0    1];

M_true  = T * R; %The true transformation matrix

%Point set 1
V1=Vd+0.5;

%Point set 2
V2=tform(M_true,V1);

%%
% Plotting input data

hf=cFigure;
title('The untransformed surfaces','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);

hp=patch('Faces',F,'Vertices',V1,'FaceColor','g','FaceAlpha',faceAlpha);

hp=patch('Faces',F,'Vertices',V2,'FaceColor','r','FaceAlpha',faceAlpha,'edgeColor',edgeColor);

camlight headlight;
set(gca,'FontSize',fontSize);
view(3); axis tight;  axis equal; box on;
drawnow;

%%
% Get the transformation matrix for the point matched data using |rigidTransformationMatrixDirect|

[M_fit]=rigidTransformationMatrixDirect(V1,V2);

V1f=tform((M_fit),V1);

%%
% Plotting results

hf=cFigure;
title('The green surfaces transformed towards the red','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);

hp=patch('Faces',F,'Vertices',V2,'FaceColor','r','FaceAlpha',0.5,'edgeColor',edgeColor);

hp=patch('Faces',F,'Vertices',V1f,'FaceColor','none','FaceAlpha',0.5,'edgeColor','g');
% hp=patch('Faces',F,'Vertices',V2ff,'FaceColor','none','FaceAlpha',0.5*2,'edgeColor','b');

set(gca,'FontSize',fontSize);
view(3); axis tight;  axis equal; box on;
drawnow;

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

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


%  Script file: calc_roots2.m
%
%  Purpose:
%    This program solves for the roots of a quadratic equation
%    of the form a*x**2 + b*x + c = 0.  It calculates the answers
%    regardless of the type of roots that the equation possesses.
%
%  Record of revisions:
%      Date       Programmer          Description of change
%      ====       ==========          =====================
%    02/24/07    S. J. Chapman        Original code
%
% Define variables:
%   a            -- Coefficient of x^2 term of equation
%   b            -- Coefficient of x term of equation
%   c            -- Constant term of equation
%   discriminant -- Discriminant of the equation
%   x1           -- First solution of equation
%   x2           -- Second solution of equation

% Prompt the user for the coefficients of the equation
disp ('This program solves for the roots of a quadratic ');
disp ('equation of the form A*X^2 + B*X + C = 0. ');
a = input ('Enter the coefficient A: ');
b = input ('Enter the coefficient B: ');
c = input ('Enter the coefficient C: ');

% Calculate discriminant
discriminant = b^2 - 4 * a * c;

% Solve for the roots
x1 = ( -b + sqrt(discriminant) ) / ( 2 * a );
x2 = ( -b - sqrt(discriminant) ) / ( 2 * a );

% Display results
disp ('The roots of this equation are:');
fprintf ('x1 = (%f) +i (%f)\n', real(x1), imag(x1));
fprintf ('x2 = (%f) +i (%f)\n', real(x2), imag(x2));


text:
classdef TestMergeMertens
    %TestMergeMertens

    methods (Static)
        function test_1
            fpath = fullfile(mexopencv.root(),'test');
            files = dir(fullfile(fpath,'memorial*.png'));
            imgs = cell(1,numel(files));
            for i=1:numel(files)
                imgs{i} = imread(fullfile(fpath,files(i).name));
            end

            obj = cv.MergeMertens();
            ldr = obj.process(imgs);
            validateattributes(ldr, {'single'}, {'size',size(imgs{1})});
            %TODO: output LDR is not always in [0,1] range
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [K] = kv1v1(x, xp, hyp, i)

logsigmau = hyp(1);
logthetau = hyp(2);
logsigmav = hyp(3);
logthetav = hyp(4);

n_x = size(x,1);
n_xp = size(xp,1);

x = repmat(x,1,n_xp);
xp = repmat(xp',n_x,1);

switch i

case 0

K=exp(1).^(logsigmav+(-1/2).*exp(1).^((-1).*logthetav).*(x+(-1).*xp).^2); ...

case 1 % logsigmau

K=0;

case 2 % logthetau

K=0;

case 3 % logsigmav

K=exp(1).^(logsigmav+(-1/2).*exp(1).^((-1).*logthetav).*(x+(-1).*xp).^2); ...

case 4 % logthetav

K=(1/2).*exp(1).^(logsigmav+(-1).*logthetav+(-1/2).*exp(1).^((-1).* ...
  logthetav).*(x+(-1).*xp).^2).*(x+(-1).*xp).^2;

otherwise

        K = zeros(n_x, n_xp);
end

if K == 0

    K = zeros(n_x, n_xp);

end

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
function [Population,Dec,Mask,FrontNo,CrowdDis] = EnvironmentalSelection(Population,Dec,Mask,N)
% The environmental selection of SLMEA

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Delete duplicated solutions
    objs = Population.objs;
    objs = gather(objs);

    [objs,uni] = unique(objs,'rows');
    Population = Population(uni);
    Dec        = Dec(uni,:);
    Mask       = Mask(uni,:);
    N          = min(N,length(Population));
    %% Non-dominated sorting
    [FrontNo,MaxFNo] = NDSort(objs,gather(Population.cons),N);
    Next = FrontNo < MaxFNo;

    %% Calculate the crowding distance of each solution
    CrowdDis = CrowdingDistance(objs,FrontNo);

    %% Select the solutions in the last front based on their crowding distances
    Last     = find(FrontNo==MaxFNo);
    [~,Rank] = sort(CrowdDis(Last),'descend');
    Next(Last(Rank(1:N-sum(Next)))) = true;

    %% Population for next generation
    Population = Population(Next);
    FrontNo    = FrontNo(Next);
    CrowdDis   = CrowdDis(Next);
    Dec        = Dec(Next,:);
    Mask       = Mask(Next,:);
end

author: BIMK


text:
function [sourceout] = ft_source2sparse(sourcein)

% FT_SOURCE2SPARSE removes the grid locations outside the brain from the source
% reconstruction, thereby saving memory.
%
% This invalidates the fields that describe the grid, and also makes it
% more difficult to make a plot of each of the slices of the source volume.
% The original source structure can be recreated using FT_SOURCE2FULL.
%
% Use as
%   [source] = ft_source2sparse(source)
%
% See also FT_SOURCE2FULL

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% this function does not take a cfg as input
cfg = [];

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble provenance sourcein

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

%%

sourcein = ft_checkdata(sourcein, 'datatype', 'source', 'insidestyle', 'logical');

inside  = find( sourcein.inside(:));
outside = find(~sourcein.inside(:));

fprintf('total number of dipoles        : %d\n', length(inside)+length(outside));
fprintf('number of dipoles inside  brain: %d\n', length(inside));
fprintf('number of dipoles outside brain: %d\n', length(outside));

%% construct the output data

% loop over parameters
fn = fieldnames(sourcein);
fn = setdiff(fn, {'dim', 'transform'});

% keep the descriptive fields that remain valid
sourceout = keepfields(sourcein, {'dim', 'transform', 'time', 'freq'});

for i=1:numel(fn)
  dimord = getdimord(sourcein, fn{i});

  clear tmp
  if startsWith(dimord, 'pos_pos') % this should go first
    tmp = sourcein.(fn{i});
    tmp = tmp(inside,inside,:,:,:,:,:);
    sourceout.(fn{i}) = tmp;
  elseif startsWith(dimord, 'pos')
    tmp = sourcein.(fn{i});
    tmp = tmp(inside,:,:,:,:,:,:);
    sourceout.(fn{i}) = tmp;
  elseif startsWith(dimord, '{pos}')
    tmp = sourcein.(fn{i});
    tmp = tmp(inside);
    sourceout.(fn{i}) = tmp;
  elseif startsWith(dimord, '{pos_pos}')
    tmp = sourcein.(fn{i});
    tmp = tmp(inside,inside);
    sourceout.(fn{i}) = tmp;
  else
    % skipping parameters with unknown dimord
  end

  if exist('tmp', 'var')
    ft_info('making "%s" with dimord "%s" sparse', fn{i}, dimord);
  else
    ft_info('skipping "%s" with dimord "%s"', fn{i}, dimord);
  end

end % for each field

%% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   sourcein
ft_postamble provenance sourceout
ft_postamble history    sourceout

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function so = streamAccumulate(si,start,accumulator)
% STREAMACCUMULATE  Stream of accumulation of values
so = {start,delayEval(@streamAccumulate, ...
        {streamCdr(si), feval(accumulator,start,streamCar(si)),accumulator})};

author: Sable

repo: mcbench-benchmarks


text:
function type = tt_assignment_type(x,subs,rhs)
%TT_ASSIGNMENT_TYPE What type of subsasgn is this?
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

if isequal(class(x),class(rhs))
    type = 'subtensor';
    return;
end

if (numel(subs) >= 2)
    type = 'subtensor';
    return;
end

type = 'subscripts';

author: andrewssobral


load('data/bal2.mat');
[L,Det]=lmnn(xTr,yTr,'quiet',1);
enerr=energyclassify(L,xTr,yTr,xTe,yTe,3);
knnerrL=knnclassify(L,xTr,yTr,xTe,yTe,3);
knnerrI=knnclassify(eye(size(L)),xTr,yTr,xTe,yTe,3);

clc;
fprintf('Bal data set:\n');
fprintf('3-NN Euclidean training error: %2.2f\n',knnerrI(1)*100);
fprintf('3-NN Euclidean testing error: %2.2f\n',knnerrI(2)*100);
fprintf('3-NN Malhalanobis training error: %2.2f\n',knnerrL(1)*100);
fprintf('3-NN Malhalanobis testing error: %2.2f\n',knnerrL(2)*100);
fprintf('Energy classification error: %2.2f\n',enerr*100);
fprintf('Training time: %2.2fs\n (As a reality check: My desktop needs 20.53s)\n\n',Det.time);


text:
function varargout = sample(f, varargin)
%SAMPLE   Samples a DISKFUN object on a tensor product grid.
%   X = SAMPLE(F) returns the matrix of values of F(theta, r) on a
%   Fourier-Chebyshev tensor product grid. (theta, r) are polar coordinates,
%   with -pi <= theta <= pi and 0 <= r <= 1.
%
%   [U, D, V] = SAMPLE(F) returns the low rank representation of the
%   values of F on a tensor product grid where X = U * D * V'.
%
%   [U, D, V] = SAMPLE(F,M,N) returns the values of F on an M-by-N
%   tensor product grid.

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

% Empty check.
if ( isempty(f) )
    varargout = { [] };
    return
end

if ( nargin == 1 )
    % Get degrees:
    [m, n] = length(f);

elseif ( nargin == 2 )
    error('CHEBFUN:DISKFUN:sample:inputs', 'Dimension not specified.');

else
    m = varargin{ 1 };
    n = varargin{ 2 };
    if ( (m <= 0) || (n <= 0) )
        error('CHEBFUN:DISKFUN:sample:inputs', ['Number of sample ' ...
             'points must be positive.']);
    end
end

% Get the low rank representation for f.
[cols, d, rows] = cdr(f);

C = sample(cols, max(2*n-1, 1));
C = C(n:end, :);

R = real( sample(rows, m));

% Evaluate:
if ( nargout <= 1 )
    varargout = {C * d * R.'};
else
    varargout = {C , d, R};
end

end


text:
function pdf = pareto_pdf ( x, a, b )

%*****************************************************************************80
%
%% PARETO_PDF evaluates the Pareto PDF.
%
%  Formula:
%
%    PDF(X)(A,B) = B * A**B / X**(B+1).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%    A <= X
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < A.
%    0.0 < B.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x < a )
    pdf = 0.0;
  else
    pdf = b * a^b / x^( b + 1.0 );
  end

  return
end

author: johannesgerer


text:
function [mean_m1, b1, sig1, av2] =  bmemag(mag);
    %BMEMAG
    maximum_mag = max(mag );
    minimum_mag = min(mag );
    if minimum_mag > 0 ; minimum_mag = 0 ; end

    % calculate the mean magnitude, b(mean) and std
    n = length(mag );
    mean_m1 = mean(mag );
    b1 = (1/(mean_m1-min(mag -0.05)))*log10(exp(1));
    sig1 = (sum((mag -mean_m1).^2))/(n*(n-1));
    sig1 = sqrt(sig1);
    sig1 = 2.30*sig1*b1^2;            % standard deviation
    %disp ([' b-value segment 1 = ' num2str(b1) ]);
    %disp ([' standard dev b_val_1 = ' num2str(sig1) ]);
    av2 = log10(length(mag ))+b1*min(mag );
end

author: geoscience-community-codes

repo: GISMO


text:
function [U] = dualquatlbs(V,DQ,W)
  % DUALQUATLBS Compute dual quaternions linear blend skinning deformation of
  % vertices V, using rigid transformations stored as dual quaternions, DQ, at
  % some control points, propogated to the mesh using weights W.
  %
  % [U] = dualquatlbs(V,DQ,W)
  %
  % Inputs:
  %  V  list of vertex positions
  %  DQ  list of rigid transformations for each control point stored as dual
  %    quaternions
  %    2 by 4 by #controls
  %  W  weights, # vertices by # handles matrix of weights
  % Output:
  %  U  list of new vertex positions
  %
  % Copyright 2011, Alec Jacobson (jacobson@inf.ethz.ch)
  %
  % See also: lbs
  %

  % pad V with zeros if not 3D
  if(size(V,2) ~= 3)
    was_2d = true;
    V = [V zeros(size(V,1),1)];
  else
    was_2d = false;
  end

  % number of control points
  m = size(DQ,3);
  % should be same in W
  assert(m == size(W,2));

  % number of domain vertices
  n = size(V,1);

  % See algorithm 1 in "Geometric skinning with approximate dual quaternion
  % blending" by Kavan et al

  % compute weighted combination of DQs for each domain vertex
  % DQs seen by every vertex in domain, 2 by 4 by m by n
  WDQ = permute(repmat(W',[1,1,2,4]),[3 4 1 2]) .* repmat(DQ,[1,1,1,n]);
  % sum of weighted dual quaternions, 2 by 4 by n
  VDQ = permute(sum(WDQ,3),[1 2 4 3]);
  % regular part, n by 4
  VDQ1 = permute(VDQ(1,:,:),[3 2 1]);
  % dual part, n by 4
  VDQ2 = permute(VDQ(2,:,:),[3 2 1]);
  %VDQ2 = permute(sum(WDQ(2,:,:,:),3),[4 2 1 3]);
  len = repmat(sqrt(sum(VDQ1.^2,2)),1,4);
  VDQ1 = VDQ1./len;
  VDQ2 = VDQ2./len;

  U = ...
    V + ...
    2*cross( ...
      VDQ1(:,2:4), ...
      cross(VDQ1(:,2:4),V,2)+repmat(VDQ1(:,1),1,3).*V,2) + ...
    2*( ...
      repmat(VDQ1(:,1),1,3).*VDQ2(:,2:4) - ...
      repmat(VDQ2(:,1),1,3).*VDQ1(:,2:4) + ...
      cross(VDQ1(:,2:4),VDQ2(:,2:4),2));

  %VDQ = zeros([2 size(VDQ1)]);
  %VDQ(1,:,:) = VDQ1;
  %VDQ(2,:,:) = VDQ2;
  %% convert DQs into transformations, multiply against vertices
  %% convert rotation to quaternion and translation
  %[Q,T] = udq2quattrans(VDQ);
  %% convert quaterion
  %R = quat2mat(Q);

  % unpad U if V was 2D
  if(was_2d)
    U = U(:,1:2);
  end
end

author: alecjacobson

repo: gptoolbox


text:
% Compute the average degree of a node in a graph, defined as
% 2 times the number of edges divided by the number of nodes
%          (every edge is counted towards the degrees twice).
%
% Inputs: adjacency matrix, nxn
% Outputs: float, the average degree, a number between 0 and max(sum(adj))
%
% Note: The average degree is related to the link density, namely:
%       link_density = ave_degree/(n-1), where n is the number of nodes
%
% Other routines used: numNodes.m, numEdges.m
% GB: last update, September 20, 2012

function k=averageDegree(adj)

k=2*numEdges(adj)/numNodes(adj);

author: aeolianine

repo: octave-networks-toolbox


min_cost_0_1000 = min_cost_0;
min_cost_1_1000 = min_cost_1;
min_cost_2_1000 = min_cost_2;
min_cost_3_1000 = min_cost_3;
min_cost_4_1000 = min_cost_4;
min_cost_5_1000 = min_cost_5;
min_cost_6_1000 = min_cost_6;

totalNum_0_1000 = totalNum_0;
totalNum_1_1000 = totalNum_1;
totalNum_2_1000 = totalNum_2;
totalNum_3_1000 = totalNum_3;
totalNum_4_1000 = totalNum_4;
totalNum_5_1000 = totalNum_5;
totalNum_6_1000 = totalNum_6;


text:
function f = spm_Bpdf(x,v,w)
% Probability Density Function (PDF) of Beta distribution
% FORMAT f = spm_Bpdf(x,v,w)
%
% x   - Beta variates (Beta has range [0,1])
% v   - Shape parameter (v>0)
% w   - Shape parameter (w>0)
% F   - PDF of Beta distribution with shape parameters [v,w] at points x
%__________________________________________________________________________
%
% spm_Bpdf implements the Probability Density Function for Beta distributions.
%
% Definition:
%--------------------------------------------------------------------------
% The PDF of the Beta distribution shape parameters v & w, defined
% for positive integer degrees of freedom v>0 & w>0, and for x in
% [0,1] is given by: (See Evans et al., Ch5)
%
%            x^(v-1) * (1-x)^(w-1)
%    f(x) = -----------------------
%                 beta(v,w)
%
% Variate relationships:
%--------------------------------------------------------------------------
% Many: See Evans et al., Ch5
%
% Algorithm:
%--------------------------------------------------------------------------
% Direct computation using logs and MATLAB's implementation of the log
% beta function (betaln).
%
% References:
%--------------------------------------------------------------------------
% Evans M, Hastings N, Peacock B (1993)
%       "Statistical Distributions"
%        2nd Ed. Wiley, New York
%
% Abramowitz M, Stegun IA, (1964)
%       "Handbook of Mathematical Functions"
%        US Government Printing Office
%
% Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
%       "Numerical Recipes in C"
%        Cambridge
%__________________________________________________________________________
% Copyright (C) 1999-2011 Wellcome Trust Centre for Neuroimaging

% Andrew Holmes
% $Id: spm_Bpdf.m 4182 2011-02-01 12:29:09Z guillaume $

%-Format arguments, note & check sizes
%--------------------------------------------------------------------------
if nargin<3, error('Insufficient arguments'), end

ad = [ndims(x);ndims(v);ndims(w)];
rd = max(ad);
as = [[size(x),ones(1,rd-ad(1))];...
      [size(v),ones(1,rd-ad(2))];...
      [size(w),ones(1,rd-ad(3))]];
rs = max(as);
xa = prod(as,2)>1;
if sum(xa)>1 && any(any(diff(as(xa,:)),1))
    error('non-scalar args must match in size');
end

%-Computation
%--------------------------------------------------------------------------
%-Initialise result to zeros
f = zeros(rs);

%-Only defined for x in [0,1] & strictly positive v & w.
% Return NaN if undefined.
md = ( x>=0  &  x<=1  &  v>0  &  w>0 );
if any(~md(:))
    f(~md) = NaN;
    warning('Returning NaN for out of range arguments');
end

%-Special cases: x=0 & x=1
f(md & x==0 & v<1) = Inf;
f(md & x==1 & w<1) = Inf;

%-Non-zero where defined & x in (0,1)
Q  = find( md  &  x>0  &  x<1 );
if isempty(Q), return, end
if xa(1), Qx=Q; else Qx=1; end
if xa(2), Qv=Q; else Qv=1; end
if xa(3), Qw=Q; else Qw=1; end

%-Compute
f(Q) = exp((v(Qv)-1).*log(x(Qx)) + (w(Qw)-1).*log(1-x(Qx)) -betaln(v(Qv),w(Qw)));

author: spm


text:
function prob_test034 ( )

%*****************************************************************************80
%
%% TEST034 tests CHI_SQUARE_MEAN, CHI_SQUARE_SAMPLE, CHI_SQUARE_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST034\n' );
  fprintf ( 1, '  For the central chi square PDF:\n' );
  fprintf ( 1, '  CHI_SQUARE_MEAN computes the mean;\n' );
  fprintf ( 1, '  CHI_SQUARE_SAMPLE samples;\n' );
  fprintf ( 1, '  CHI_SQUARE_VARIANCE computes the variance.\n' );

  a = 10.0;

  check = chi_square_check ( a );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST034 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  mean = chi_square_mean ( a );
  variance = chi_square_variance ( a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A = %14f\n', a );
  fprintf ( 1, '  PDF mean =        %14f\n', mean );
  fprintf ( 1, '  PDF variance =    %14f\n', variance );

  for i = 1 : nsample
    [ x(i), seed ] = chi_square_sample ( a, seed );
  end

  mean = r8vec_mean ( nsample, x );
  variance = r8vec_variance ( nsample, x );
  xmax = max ( x(1:nsample) );
  xmin = min ( x(1:nsample) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %6d\n', nsample );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %14f\n', xmax );
  fprintf ( 1, '  Sample minimum =  %14f\n', xmin );

  return
end

repo: jburkardt-m


text:
function I = make_odd_by_odd(I)
%MAKE_ODD_BY_ODD strips a row and or a column so that it will be
%an odd number of pixels in each dimension.

if ~isodd(size(I,1))
    I = I(1:end -1,:,:);
end

if ~isodd(size(I,2))
    I = I(:, 1:end -1,:);
end
% Copyright 2002 - 2009 The MathWorks, Inc.


text:
function [Ep,Cp,Eh,F,L,dFdp,dFdpp] = spm_nlsi_GN(M,U,Y)
% Bayesian inversion of nonlinear models - Gauss-Newton/Variational Laplace
% FORMAT [Ep,Cp,Eh,F] = spm_nlsi_GN(M,U,Y)
%
% [Dynamic] MIMO models
%__________________________________________________________________________
%
% M.IS - function name f(P,M,U) - generative model
%        This function specifies the nonlinear model:
%          y = Y.y = IS(P,M,U) + X0*P0 + e
%        where e ~ N(0,C). For dynamic systems this would be an integration
%        scheme (e.g. spm_int). spm_int expects the following:
%
%     M.f  - f(x,u,P,M)
%     M.g  - g(x,u,P,M)
%     M.h  - h(x,u,P,M)
%       x  - state variables
%       u  - inputs or causes
%       P  - free parameters
%       M  - fixed functional forms and parameters in M
%
% M.FS - function name f(y,M)   - feature selection
%        This [optional] function performs feature selection assuming the
%        generalized model y = FS(y,M) = FS(IS(P,M,U),M) + X0*P0 + e
%
% M.P  - starting estimates for model parameters [optional]
%
% M.pE - prior expectation      - E{P}   of model parameters
% M.pC - prior covariance       - Cov{P} of model parameters
%
% M.hE - prior expectation      - E{h}   of log-precision parameters
% M.hC - prior covariance       - Cov{h} of log-precision parameters
%
% U.u  - inputs (or just U)
% U.dt - sampling interval
%
% Y.y  - outputs (samples x observations x ...)
% Y.dt - sampling interval for outputs
% Y.X0 - confounds or null space      (over size(y,1) samples or all vec(y))
% Y.Q  - q error precision components (over size(y,1) samples or all vec(y))
%
%
% Parameter estimates
%--------------------------------------------------------------------------
% Ep  - (p x 1)         conditional expectation    E{P|y}
% Cp  - (p x p)         conditional covariance     Cov{P|y}
% Eh  - (q x 1)         conditional log-precisions E{h|y}
%
% log evidence
%--------------------------------------------------------------------------
% F   - [-ve] free energy F = log evidence = p(y|f,g,pE,pC) = p(y|m)
%
%__________________________________________________________________________
% Returns the moments of the posterior p.d.f. of the parameters of a
% nonlinear model specified by IS(P,M,U) under Gaussian assumptions.
% Usually, IS is an integrator of a dynamic MIMO input-state-output model
%
%              dx/dt = f(x,u,P)
%              y     = g(x,u,P)  + X0*P0 + e
%
% A static nonlinear observation model with fixed input or causes u
% obtains when x = []. i.e.
%
%              y     = g([],u,P) + X0*P0e + e
%
% but static nonlinear models are specified more simply using
%
%              y     = IS(P,M,U) + X0*P0 + e
%
% Priors on the free parameters P are specified in terms of expectation pE
% and covariance pC. The E-Step uses a Fisher-Scoring scheme and a Laplace
% approximation to estimate the conditional expectation and covariance of P
% If the free-energy starts to increase,  an abbreviated descent is
% invoked.  The M-Step estimates the precision components of e, in terms
% of log-precisions.  Although these two steps can be thought of in
% terms of E and N steps they are in fact variational steps of a full
% variational Laplace scheme that accommodates conditional uncertainty
% over both parameters and log precisions (c.f. hyperparameters with hyper
% priors)
%
% An optional feature selection can be specified with parameters M.FS.
%
% For generic aspects of the scheme see:
%
% Friston K, Mattout J, Trujillo-Barreto N, Ashburner J, Penny W.
% Variational free energy and the Laplace approximation.
% NeuroImage. 2007 Jan 1;34(1):220-34.
%
% This scheme handels complex data along the lines originally described in:
%
% Sehpard RJ, Lordan BP, and Grant EH.
% Least squares analysis of complex data with applications to permittivity
% measurements.
% J. Phys. D. Appl. Phys 1970 3:1759-1764.
%
%__________________________________________________________________________
% Copyright (C) 2001-2015 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_nlsi_GN.m 7279 2018-03-10 21:22:44Z karl $

% options
%--------------------------------------------------------------------------
try, M.nograph; catch, M.nograph = 0;   end
try, M.noprint; catch, M.noprint = 0;   end
try, M.Nmax;    catch, M.Nmax    = 128; end

% figure (unless disabled)
%--------------------------------------------------------------------------
if ~M.nograph
    Fsi = spm_figure('GetWin','SI');
end

% check integrator
%--------------------------------------------------------------------------
try
    M.IS;
catch
    M.IS = 'spm_int';
end

% composition of feature selection and prediction (usually an integrator)
%--------------------------------------------------------------------------
try
    y  = Y.y;
catch
    y  = Y;
end

try

    % try FS(y,M)
    %----------------------------------------------------------------------
    try
        y  = feval(M.FS,y,M);
        IS = inline([M.FS '(' M.IS '(P,M,U),M)'],'P','M','U');

        % try FS(y)
        %------------------------------------------------------------------
    catch
        y  = feval(M.FS,y);
        IS = inline([M.FS '(' M.IS '(P,M,U))'],'P','M','U');
    end

catch

    % otherwise FS(y) = y
    %----------------------------------------------------------------------
    try
        IS = inline([M.IS '(P,M,U)'],'P','M','U');
    catch
        IS = M.IS;
    end
end

% converted to function handle
%--------------------------------------------------------------------------
IS  = spm_funcheck(IS);

% paramter update eqation
%--------------------------------------------------------------------------
if isfield(M,'f'), M.f = spm_funcheck(M.f);  end
if isfield(M,'g'), M.g = spm_funcheck(M.g);  end
if isfield(M,'h'), M.h = spm_funcheck(M.h);  end

% size of data (samples x response component x response component ...)
%--------------------------------------------------------------------------
if iscell(y)
    ns = size(y{1},1);
else
    ns = size(y,1);
end
ny   = length(spm_vec(y));          % total number of response variables
nr   = ny/ns;                       % number response components
M.ns = ns;                          % number of samples M.ns

% initial states
%--------------------------------------------------------------------------
try
    M.x;
catch
    if ~isfield(M,'n'), M.n = 0;    end
    M.x = sparse(M.n,1);
end

% input
%--------------------------------------------------------------------------
try
    U;
catch
    U = [];
end

% initial parameters
%--------------------------------------------------------------------------
try
    spm_vec(M.P) - spm_vec(M.pE);
    fprintf('\nParameter initialisation successful\n')
catch
    M.P = M.pE;
end

% time-step
%--------------------------------------------------------------------------
try
    dt = Y.dt;
catch
    dt = 1;
end

% precision components Q
%--------------------------------------------------------------------------
try
    Q = Y.Q;
    if isnumeric(Q), Q = {Q}; end
catch
    Q = spm_Ce(ns*ones(1,nr));
end
nh    = length(Q);                  % number of precision components
nq    = ny/length(Q{1});            % for compact Kronecker form of M-step

% prior moments (assume uninformative priors if not specifed)
%--------------------------------------------------------------------------
pE       = M.pE;
try
    pC   = M.pC;
catch
    np   = spm_length(M.pE);
    pC   = speye(np,np)*exp(16);
end

% confounds (if specified)
%--------------------------------------------------------------------------
try
    nb   = size(Y.X0,1);            % number of bins
    nx   = ny/nb;                   % number of blocks
    dfdu = kron(speye(nx,nx),Y.X0);
catch
    dfdu = sparse(ny,0);
end
if isempty(dfdu), dfdu = sparse(ny,0); end

% hyperpriors - expectation (and initialize hyperparameters)
%--------------------------------------------------------------------------
try
    hE  = M.hE;
    if length(hE) ~= nh
        hE = hE + sparse(nh,1);
    end
catch
    hE  = sparse(nh,1) - log(var(spm_vec(y))) + 4;
end
h       = hE;

% hyperpriors - covariance
%--------------------------------------------------------------------------
try
    ihC = spm_inv(M.hC);
    if length(ihC) ~= nh
        ihC = ihC*speye(nh,nh);
    end
catch
    ihC = speye(nh,nh)*exp(4);
end

% unpack covariance
%--------------------------------------------------------------------------
if isstruct(pC);
    pC = spm_diag(spm_vec(pC));
end

% dimension reduction of parameter space
%--------------------------------------------------------------------------
V     = spm_svd(pC,0);
nu    = size(dfdu,2);                 % number of parameters (confounds)
np    = size(V,2);                    % number of parameters (effective)
ip    = (1:np)';
iu    = (1:nu)' + np;

% second-order moments (in reduced space)
%--------------------------------------------------------------------------
pC    = V'*pC*V;
uC    = speye(nu,nu)/1e-8;
ipC   = inv(spm_cat(spm_diag({pC,uC})));

% initialize conditional density
%--------------------------------------------------------------------------
Eu    = spm_pinv(dfdu)*spm_vec(y);
p     = [V'*(spm_vec(M.P) - spm_vec(M.pE)); Eu];
Ep    = spm_unvec(spm_vec(pE) + V*p(ip),pE);

% EM
%==========================================================================
criterion = [0 0 0 0];

C.F   = -Inf;                                   % free energy
v     = -4;                                     % log ascent rate
dFdh  = zeros(nh,1);
dFdhh = zeros(nh,nh);
for k = 1:M.Nmax

    % time
    %----------------------------------------------------------------------
    tStart = tic;

    % E-Step: prediction f, and gradients; dfdp
    %======================================================================
    try

        % gradients
        %------------------------------------------------------------------
        [dfdp,f] = spm_diff(IS,Ep,M,U,1,{V});
        dfdp     = reshape(spm_vec(dfdp),ny,np);

        % check for stability
        %------------------------------------------------------------------
        normdfdp = norm(dfdp,'inf');
        revert   = isnan(normdfdp) || normdfdp > exp(32);

    catch
        revert   = true;
    end

    if revert && k > 1
        for i = 1:4

            % reset expansion point and increase regularization
            %--------------------------------------------------------------
            v        = min(v - 2,-4);

            % E-Step: update
            %--------------------------------------------------------------
            p        = C.p + spm_dx(dFdpp,dFdp,{v});
            Ep       = spm_unvec(spm_vec(pE) + V*p(ip),pE);

            % try again
            %--------------------------------------------------------------
            try
                [dfdp,f] = spm_diff(IS,Ep,M,U,1,{V});
                dfdp     = reshape(spm_vec(dfdp),ny,np);

                % check for stability
                %----------------------------------------------------------
                normdfdp = norm(dfdp,'inf');
                revert   = isnan(normdfdp) || normdfdp > exp(32);

            catch
                revert   = true;
            end

            % break
            %--------------------------------------------------------------
            if ~revert, break, end

        end
    end

    % convergence failure
    %----------------------------------------------------------------------
    if revert
        error('SPM:spm_nlsi_GN','Convergence failure.');
    end

    % prediction error and full gradients
    %----------------------------------------------------------------------
    e     =  spm_vec(y) - spm_vec(f) - dfdu*p(iu);
    J     = -[dfdp dfdu];

    % M-step: Fisher scoring scheme to find h = max{F(p,h)}
    %======================================================================
    for m = 1:8

        % precision and conditional covariance
        %------------------------------------------------------------------
        iS    = sparse(0);
        for i = 1:nh
            iS = iS + Q{i}*(exp(-32) + exp(h(i)));
        end
        S     = spm_inv(iS);
        iS    = kron(speye(nq),iS);
        Pp    = real(J'*iS*J);
        Cp    = spm_inv(Pp + ipC);

        % precision operators for M-Step
        %------------------------------------------------------------------
        for i = 1:nh
            P{i}   = Q{i}*exp(h(i));
            PS{i}  = P{i}*S;
            P{i}   = kron(speye(nq),P{i});
            JPJ{i} = real(J'*P{i}*J);
        end

        % derivatives: dLdh = dL/dh,...
        %------------------------------------------------------------------
        for i = 1:nh
            dFdh(i,1)      =   trace(PS{i})*nq/2 ...
                - real(e'*P{i}*e)/2 ...
                - spm_trace(Cp,JPJ{i})/2;
            for j = i:nh
                dFdhh(i,j) = - spm_trace(PS{i},PS{j})*nq/2;
                dFdhh(j,i) =   dFdhh(i,j);
            end
        end

        % add hyperpriors
        %------------------------------------------------------------------
        d     = h     - hE;
        dFdh  = dFdh  - ihC*d;
        dFdhh = dFdhh - ihC;
        Ch    = spm_inv(real(-dFdhh));

        % update ReML estimate
        %------------------------------------------------------------------
        dh    = spm_dx(dFdhh,dFdh,{4});
        dh    = min(max(dh,-1),1);
        h     = h  + dh;

        % convergence
        %------------------------------------------------------------------
        dF    = dFdh'*dh;
        if dF < 1e-2, break, end

    end

    % E-Step with Levenberg-Marquardt regularization
    %======================================================================

    % objective function: F(p) = log evidence - divergence
    %----------------------------------------------------------------------
    L(1) = spm_logdet(iS)*nq/2  - real(e'*iS*e)/2 - ny*log(8*atan(1))/2;            ...
    L(2) = spm_logdet(ipC*Cp)/2 - p'*ipC*p/2;
    L(3) = spm_logdet(ihC*Ch)/2 - d'*ihC*d/2;
    F    = sum(L);

    % record increases and reference log-evidence for reporting
    %----------------------------------------------------------------------
    try
        F0;
        if ~M.noprint
            fprintf(' actual: %.3e (%.2f sec)\n',full(F - C.F),toc(tStart))
        end
    catch
        F0 = F;
    end

    % if F has increased, update gradients and curvatures for E-Step
    %----------------------------------------------------------------------
    if F > C.F || k < 3

        % accept current estimates
        %------------------------------------------------------------------
        C.p   = p;
        C.h   = h;
        C.F   = F;
        C.L   = L;
        C.Cp  = Cp;

        % E-Step: Conditional update of gradients and curvature
        %------------------------------------------------------------------
        dFdp  = -real(J'*iS*e) - ipC*p;
        dFdpp = -real(J'*iS*J) - ipC;

        % decrease regularization
        %------------------------------------------------------------------
        v     = min(v + 1/2,4);
        str   = 'EM:(+)';

    else

        % reset expansion point
        %------------------------------------------------------------------
        p     = C.p;
        h     = C.h;
        Cp    = C.Cp;

        % and increase regularization
        %------------------------------------------------------------------
        v     = min(v - 2,-4);
        str   = 'EM:(-)';

    end

    % E-Step: update
    %======================================================================
    dp    = spm_dx(dFdpp,dFdp,{v});
    p     = p + dp;
    Ep    = spm_unvec(spm_vec(pE) + V*p(ip),pE);

    % Graphics
    %======================================================================
    if exist('Fsi', 'var')
        spm_figure('Select', Fsi)

        % reshape prediction if necessary
        %------------------------------------------------------------------
        e  = spm_vec(e);
        f  = spm_vec(f);
        try
            e  = reshape(e,ns,nr);
            f  = reshape(f,ns,nr);
        end

        % subplot prediction
        %------------------------------------------------------------------
        x    = (1:size(e,1))*dt;
        xLab = 'time (seconds)';
        try
            if length(M.Hz) == ns
                x    = M.Hz;
                xLab = 'Frequency (Hz)';
            end
        end

        % plot real or complex predictions
        %------------------------------------------------------------------
        tstr = sprintf('%s: %i','prediction and response: E-Step',k);
        if isreal(spm_vec(y))

            subplot(2,1,1)
            plot(x,real(f)), hold on
            plot(x,real(f + e),':'), hold off
            xlabel(xLab)
            title(tstr,'FontSize',16)
            grid on

        else
            subplot(2,2,1)
            plot(x,real(f)), hold on
            plot(x,real(f + e),':'), hold off
            xlabel(xLab)
            ylabel('real')
            title(tstr,'FontSize',16)
            grid on

            subplot(2,2,2)
            plot(x,imag(f)), hold on
            plot(x,imag(f + e),':'), hold off
            xlabel(xLab)
            ylabel('imaginary')
            title(tstr,'FontSize',16)
            grid on
        end

        % subplot parameters
        %--------------------------------------------------------------
        subplot(2,1,2)
        bar(full(V*p(ip)))
        xlabel('parameter')
        tstr = 'conditional [minus prior] expectation';
        title(tstr,'FontSize',16)
        grid on
        drawnow

    end

    % convergence
    %----------------------------------------------------------------------
    dF  = dFdp'*dp;
    if ~M.noprint
        fprintf('%-6s: %i %6s %-6.3e %6s %.3e ',str,k,'F:',full(C.F - F0),'dF predicted:',full(dF))
    end
    criterion = [(dF < 1e-1) criterion(1:end - 1)];
    if all(criterion)
        if ~M.noprint
            fprintf(' convergence\n')
        end
        break
    end

end

if exist('Fsi', 'var')
    spm_figure('Focus', Fsi)
end

% outputs
%--------------------------------------------------------------------------
Ep     = spm_unvec(spm_vec(pE) + V*C.p(ip),pE);
Cp     = V*C.Cp(ip,ip)*V';
Eh     = C.h;
F      = C.F;
L      = C.L;


text:
function [y,deriv] = exp_mv2df(w)
% This is an MV2DF. See MV2DF_API_DEFINITION.readme.
% y = exp(w), vectorized as MATLAB usually does.

if nargin==0
    test_this();
    return;
end

if isempty(w)
    y = @(w)exp_mv2df(w);
    return;
end

if isa(w,'function_handle')
    outer = exp_mv2df([]);
    y = compose_mv(outer,w,[]);
    return;
end

w = w(:);
y = exp(w);
deriv = @(dy) deriv_this(dy,y);

function [g,hess,linear] = deriv_this(dy,y)

linear  = false;
g = dy.*y;
hess = @(d) hess_this(d,dy,y);

function [h,Jv] = hess_this(d,dy,y)

h = dy.*y.*d;
if nargout>1
    Jv = d.*y;
end

function test_this()
f = exp_mv2df([]);
test_MV2DF(f,randn(3,1));

author: nesl


text:
function stop = optiplotlogfval(iter,fval,~)
%OPTIPLOTLOGVAL  Plot log10(objective) vs iterations
%
%   Default plotting callback for scalar objectives
%
%   Based on MATLAB's optimplotfval

persistent liter

%Ensure always create a new plot
if(isempty(liter))
    liter = Inf;
end

%Never stop from this function
stop = false;

%Only plot scalars
if(length(fval) > 1)
    fval = norm(fval);
end
%Take the log
fval = log10(abs(fval));

%Check for new solver run
if(iter < liter)
    h = plot(iter,fval,'b*');
    title(sprintf('Current Function Value: %g',fval));
    xlabel('Iteration /  Function Evaluation');
    set(h,'Tag','optiplotf1');
    ylabel('log10(Function Value)')
else
    h = findobj(get(gca,'Children'),'Tag','optiplotf1');
    newX = [get(h,'Xdata') iter];
    newY = [get(h,'Ydata') fval];
    set(h,'Xdata',newX, 'Ydata',newY);
    set(get(gca,'Title'),'String',sprintf('Current Function Value: %g',fval));
end

%Save for future comparison
liter = iter;

%Call drawnow to refresh the screen
drawnow;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [vertex_coords, faces, magic] = read_surf(fname, flag)
%
% [vertex_coords, faces] = read_surf(fname)
%
% or
%
% [vertex_coords, faces] = read_surf(fname, flag)
%
% reads a the vertex coordinates and face lists from a surface file
% note that reading the faces from a quad file can take a very long
% time due to the goofy format that they are stored in. If the faces
% output variable is not specified, they will not be read so it
% should execute pretty quickly.
%
% if a second optional boolean flag is specified to be true, this will determine whether
% the ASCII formatted coregistration information is read from the bottom of the file.
% This coregistration information is applied to the vertex_coords, to represent the
% surface in the coordinate system of the corresponding volume. The function's default
% is not to use this coregistration information.

% changes:
%  - 20220518: also optionally read the transformation information that is present as
%    ASCII text at the bottom of the binary file. I could not find documentation
%    online about this, but based on some limited trial and error it seems that
%    the translation information (i.e. a shift of the origin), and possibly the
%    voxel scaling can be used to represent the vertex coordinates in the coordinate
%    system of the volume from which the surface was derived. This option makes sense
%    if the surface is to be directly related to a volume, without the resulting
%    headaches of mismatching coordinate systems.
%
% read_surf.m
%
% Original Author: Bruce Fischl
% CVS Revision Info:
%    $Author: fischl $
%    $Date: 2014/04/30 12:59:03 $
%    $Revision: 1.7 $
%
% Copyright (C) 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

if nargin<2
  flag = false;
end

%fid = fopen(fname, 'r') ;
%nvertices = fscanf(fid, '%d', 1);
%all = fscanf(fid, '%d %f %f %f %f\n', [5, nvertices]) ;
%curv = all(5, :)' ;

% open it as a big-endian file

%QUAD_FILE_MAGIC_NUMBER =  (-1 & 0x00ffffff) ;
%NEW_QUAD_FILE_MAGIC_NUMBER =  (-3 & 0x00ffffff) ;

TRIANGLE_FILE_MAGIC_NUMBER =  16777214 ;
QUAD_FILE_MAGIC_NUMBER =  16777215 ;
NEW_QUAD_FILE_MAGIC_NUMBER =  16777213 ;

fid = fopen(fname, 'rb', 'b') ;
if (fid < 0)
  str = sprintf('could not open surface file %s.', fname) ;
  error(str) ;
end
magic = fread3(fid) ;

if((magic == QUAD_FILE_MAGIC_NUMBER) || (magic == NEW_QUAD_FILE_MAGIC_NUMBER))
  vnum = fread3(fid) ;
  fnum = fread3(fid) ;
  vertex_coords = fread(fid, vnum*3, 'int16') ./ 100 ;
  if (nargout > 1)
    for i=1:fnum
      for n=1:4
        faces(i,n) = fread3(fid) ;
      end
    end
  end
elseif (magic == TRIANGLE_FILE_MAGIC_NUMBER)
  fgets(fid) ;
  fgets(fid) ;
  vnum = fread(fid, 1, 'int32') ;
  fnum = fread(fid, 1, 'int32') ;
  % possibly the comment line was not followed by two \n’s, this patch is suggested by J.M.Schoffelen, 20190502
  if round(abs(vnum))~=vnum || round(abs(fnum))~=fnum
    frewind(fid);
    fread3(fid);
    fgets(fid);
    vnum = fread(fid, 1, 'int32') ;
    fnum = fread(fid, 1, 'int32') ;
  end
  vertex_coords = fread(fid, vnum*3, 'float32') ;
  if (nargout > 1)
    faces = fread(fid, fnum*3, 'int32') ;
    faces = reshape(faces, 3, fnum)' ;
  end
else
  fprintf('ERROR: magic number %d unknown\n',magic);
  vertex_coords = [];
  faces = [];
  return;
end

vertex_coords = reshape(vertex_coords, 3, vnum)' ;

if flag
  % JM the below is based a bit on trial and error, but tries to extract
  % coregistration information, that may be present as plain text at the
  % bottom of the file. I could not find any documentation online about this,
  % but anecdotally this info can be used to align the mesh with the
  % specified volume. Note that only a translation seems to be required, I am
  % not sure about the voxelsize scaling
  tline = fgetl(fid);
  if ~isempty(strfind(tline, 'volume info valid'))
    while tline~=-1
      tline = fgetl(fid);

      % the following assumes that the text lines are of a structure:
      % something = something else
      if ischar(tline)
        [v, rest] = strtok(tline);
        [equalsign, rest] = strtok(rest);
        if isequal(v, 'filename')
          [r, rest] = strtok(rest);
          eval(sprintf('%s = ''%s'';', v, rest));
        else
          eval(sprintf('%s = [%s];', v, rest));
        end
      end
    end
    R = diag(voxelsize); % not sure about this
    M = [R cras';0 0 0 1];
    tmp = [vertex_coords ones(size(vertex_coords,1),1)] * M';
    vertex_coords = tmp(:,1:3);
  end
end

fclose(fid);