text:
  function proj = ellipsoid_proj(cg, params, varargin)
%|function proj = ellipsoid_proj(cg, params, varargin)
%|
%| Compute set of 2d line-integral projection views of ellipsoid(s).
%| Works for both parallel-beam and cone-beam geometry.
%|
%| in
%|	cg			ct_geom()
%|	params [ne 9]		ellipsoid parameters:
%|			[x_center y_center z_center  x_radius y_radius z_radius
%|				xy_angle_degrees z_angle_degrees  amplitude]
%| options
%|	oversample		over-sampling factor for emulating "strips"
%|				(to account for finite detector size)
%|
%| out
%|	proj	[ns nt na]	projection views
%|
%| Copyright 2003-10-22, Patty Laskowsky, Nicole Caparanis, Taka Masuda,
%| and Jeff Fessler, University of Michigan

if nargin == 1 && streq(cg, 'test'), ellipsoid_proj_test, return, end
if nargin < 2, ir_usage, end

arg.oversample = 1;
arg = vararg_pair(arg, varargin);

proj = ellipsoid_proj_do(params, cg.s, cg.t, cg.ar, cg.source_zs, ...
		cg.dso, cg.dod, cg.dfs, arg.oversample);

end % ellipsoid_proj()

% ellipsoid_proj_do()
function proj = ellipsoid_proj_do(params, ss, tt, ...
		beta, ... % [radians]
		source_zs, dso, dod, dfs, oversample)

if size(params, 2) ~= 9, error '9 parameters per ellipsoid', end

if oversample > 1
	ds = ss(2) - ss(1);
	dt = tt(2) - tt(1);
	if any(abs(diff(ss) / ds - 1) > 1e-6) ...
	|| any(abs(diff(tt) / dt - 1) > 1e-6)
		error 'uniform spacing required for oversampling'
	end
	No = oversample;
	% determine new finer sampling positions
	ss = outer_sum([-(No-1):2:(No-1)]'/(2*No)*ds, ss(:)'); % [No ns]
	tt = outer_sum([-(No-1):2:(No-1)]'/(2*No)*dt, tt(:)'); % [No nt]
	proj = ellipsoid_proj_do(params, ss(:), tt(:), beta, source_zs, dso, dod, dfs, 1);
	proj = downsample3(proj, [No No 1]);
return
end

% determine equivalent parallel-beam projection coordinates, at beta=0
ns = length(ss);
nt = length(tt);
[sss ttt] = ndgrid(ss, tt);

if isinf(dso) % parallel beam
	uu = sss;
	vv = ttt;
	azim0 = 0;
	polar = 0;

elseif isinf(dfs) % cone-beam with flat detector
	[uu vv azim0 polar] = ir_coord_cb_flat_to_par(sss, ttt, dso, dod);

elseif dfs == 0 % cone-beam with arc detector
	[uu vv azim0 polar] = ir_coord_cb_arc_to_par(sss, ttt, dso, dod);

else
	fail 'not done'
end

clear sss ttt

cpolar = cos(polar);
spolar = sin(polar);
proj = zeros(ns, nt, numel(beta));

% loop over ellipsoids
for ip = 1:size(params,1)
	par = params(ip,:);

	cx = par(1);	rx = par(4);
	cy = par(2);	ry = par(5);
	cz = par(3);	rz = par(6);
	xang = deg2rad(par(7)); % xy-plane rotation
	zang = deg2rad(par(8)); % z-plane rotation
	if zang, error 'z rotation not done', end
	val = par(9);

	for ib = 1:length(beta)
%		az = beta(ib) + azim0 - xang; % assume source rotate in xy plane
		az = beta(ib) + azim0; % correction due to Lei Zhu of Stanford

		% shift property of 3D transform:
		cz_eff = cz - source_zs(ib); % center relative to source
		ushift = cx * cos(az) + cy * sin(az);
		vshift = (cx * sin(az) - cy * cos(az)) .* spolar + cz_eff * cpolar;

		az = az - xang; % correction due to Lei Zhu of Stanford
		p1 = (uu-ushift) .* cos(az) + (vv-vshift) .* sin(az) .* spolar;
		p2 = (uu-ushift) .* sin(az) - (vv-vshift) .* cos(az) .* spolar;
		p3 = (vv-vshift) .* cpolar;

		e1 = -sin(az) .* cpolar;
		e2 = cos(az) .* cpolar;
		e3 = spolar;

		A = e1.^2 / rx^2 + e2.^2 / ry^2 + e3.^2 / rz^2;
		B = p1 .* e1 / rx^2 + p2 .* e2 / ry^2 + p3 .* e3 / rz^2;
		C = p1.^2 / rx^2 + p2.^2 / ry^2 + p3.^2 / rz^2 - 1;

		proj(:,:,ib) = proj(:,:,ib) + 2 * val * sqrt(B.^2 - A.*C) ./ A;
	end
end

% trick: anywhere proj of a single ellipsoid is imaginary, the real part is 0.
proj = real(proj);
end % ellipsoid_proj_do()

% ellipsoid_proj_test()
function ellipsoid_proj_test

ell = [20 0*50 -40 200 100 50 90 0 10;
        0 50 100 80 80 20 0 0 10];

fun_proj = @(cg) ellipsoid_proj(cg, ell, 'oversample', 2); % analytical
fun_im = @(ig) ellipsoid_im(ig, ell, 'oversample', 2, 'checkfov', true);

ir_proj3_compare1(fun_proj, fun_im, 'chat', 1);

end % ellipsoid_proj_test()

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [nlZ,dnlZ,post,K_mat,Q] = gplite_core(hyp,gp,compute_nlZ,compute_nlZ_grad)
%GPLITE_CORE Core kernel computations for lite GP regression.

%% Initialize GP hyperparameters

[N,D] = size(gp.X);         % Number of training points and dimension

Ncov = gp.Ncov;
Nnoise = gp.Nnoise;
Nmean = gp.Nmean;

% Output warping
outwarp_flag = isfield(gp,'outwarpfun') && ~isempty(gp.outwarpfun);
if outwarp_flag
    Noutwarp = gp.Noutwarp;
    hyp_outwarp = hyp(Ncov+Nnoise+Nmean+1:Ncov+Nnoise+Nmean+Noutwarp);
    if compute_nlZ_grad
        [y,dwarp_dt,dwarp_dtheta,d2warp_dthetadt] = gp.outwarpfun(hyp_outwarp,gp.y);
    else
        [y,dwarp_dt] = gp.outwarpfun(hyp_outwarp,gp.y);
    end
    if ~isempty(gp.s2)
        s2 = gp.s2 .* dwarp_dt.^2;  % Warped noise
    else
        s2 = [];
    end
else
    y = gp.y;
    s2 = gp.s2;
end

% Evaluate observation noise on training inputs
hyp_noise = hyp(Ncov+1:Ncov+Nnoise); % Get noise hyperparameters
if compute_nlZ_grad
    [sn2,dsn2] = gplite_noisefun(hyp_noise,gp.X,gp.noisefun,gp.y,s2);
else
    sn2 = gplite_noisefun(hyp_noise,gp.X,gp.noisefun,gp.y,s2);
end
sn2_mult = 1;  % Effective noise variance multiplier

% Evaluate mean function on training inputs
hyp_mean = hyp(Ncov+Nnoise+1:Ncov+Nnoise+Nmean); % Get mean function hyperparameters
if compute_nlZ_grad
    [m,dm] = gplite_meanfun(hyp_mean,gp.X,gp.meanfun,[],gp.meanfun_extras);
else
    m = gplite_meanfun(hyp_mean,gp.X,gp.meanfun,[],gp.meanfun_extras);
end

%% Observed covariance matrix inversion

% Compute kernel matrix K_mat
if gp.covfun(1) == 1
    ell = exp(hyp(1:D));
    sf2 = exp(2*hyp(D+1));
    K_mat = sq_dist(bsxfun(@rdivide,gp.X',ell));
    K_mat = sf2 * exp(-K_mat/2);
else
    hyp_cov = hyp(1:Ncov); % Get covariance function hyperparameters
    if compute_nlZ_grad
        [K_mat,dK_mat] = gplite_covfun(hyp_cov,gp.X,gp.covfun,[]);
    else
        K_mat = gplite_covfun(hyp_cov,gp.X,gp.covfun,[]);
    end
end

% Use Cholesky representation of posterior for non-small noise
Lchol = min(sn2) >= 1e-6;

if Lchol
    if isscalar(sn2)
        sn2div = sn2;
        sn2_mat = eye(N);
    else
        sn2div = min(sn2);
        sn2_mat = diag(sn2/sn2div);
    end

    for iter = 1:10
        [L,p] = chol(K_mat/(sn2div*sn2_mult)+sn2_mat);
        if p > 0; sn2_mult = sn2_mult*10; else; break; end
    end
    sl = sn2div*sn2_mult;
    if nargout > 2
        pL = L;             % L = chol(eye(n)+sW*sW'.*K)
    end
else
    if isscalar(sn2)
        sn2_mat = sn2*eye(N);
    else
        sn2_mat = diag(sn2);
    end
    for iter = 1:10     % Cholesky decomposition until it works
        [L,p] = chol(K_mat+sn2_mult*sn2_mat);
        if p > 0; sn2_mult = sn2_mult*10; else; break; end
    end
    sl = 1;
    if nargout > 2
        pL = -L\(L'\eye(N));    % L = -inv(K+inv(sW^2))
    end
end

alpha = L\(L'\(y-m)) / sl;     % alpha = inv(K_mat + diag(sn2)) * (y - m)  I

%% Integrated basis functions

if gp.intmeanfun > 0
    bb = gp.intmeanfun_mean(:);
    BB = gp.intmeanfun_var(:);

    H = gplite_intmeanfun(gp.X,gp.intmeanfun);
    plus_idx = (BB > 0);    % Non-delta parameters
    betabar = zeros(1,size(H,1));
    if any(~plus_idx)
        T_plus = diag(1./BB(plus_idx)) + H(plus_idx,:)*(L\(L'\H(plus_idx,:)')/sl);
        T_chol = chol(T_plus);
        betabar(plus_idx) = T_chol \ (T_chol' \ (bb(plus_idx)./BB(plus_idx) + H(plus_idx,:)*alpha));
        betabar(~plus_idx) = bb(~plus_idx);
    else
        T_plus = diag(1./BB) + H*(L\(L'\H')/sl);
        T_chol = chol(T_plus);
%        betabar(:) = T_plus \ (bb./BB + H*alpha);
        betabar(:) = T_chol \ (T_chol' \ (bb./BB + H*alpha));
    end
end

%% Negative log marginal likelihood computation
nlZ = []; dnlZ = []; Q = [];

if compute_nlZ
    Nhyp = size(hyp,1);

    if gp.intmeanfun > 0
        % Negative log marginal likelihood with integrated basis functions
        prec_idx = BB > 0 & isfinite(BB);
        inf_idx = isinf(BB);

        vagueall_flag = all(inf_idx);  % Vague prior on *all* basis functions?
        vagueany_flag = any(inf_idx);  % Some vague priors?
        precany_flag = any(prec_idx);  % Some precise priors?

        % Compute first quadratic term
        nu = y-m - H'*bb;
        if vagueall_flag
            nlZ_1 = nu'*alpha/2;
        else
            if precany_flag
                HBH_prec = H(prec_idx,:)'*bsxfun(@times,BB(prec_idx),H(prec_idx,:));
                N_mat = L'*L*sl + HBH_prec;
                N_chol = chol(N_mat);
                % Ninv = N_mat\eye(N);
                Ninv = N_chol\(N_chol'\eye(N));
                nlZ_1 = nu'*(Ninv*nu)/2;
            else
                nlZ_1 = nu'*(L\(L'\nu))/sl/2;
            end
        end

        % Compute second quadratic term (vague prior contribution)
        if vagueany_flag
            if ~precany_flag
                W = chol(H*(L\(L'\H'))/sl);
                % A = H'*((H*Kinv*H')\H);
                A_mat = H'*(W\(W'\H));
                nlZ_v = -nu'*(L\(L'\(A_mat*(L\(L'\nu)))))/sl^2/2;
            else
                Hinf = H(inf_idx,:);
                W = chol(Hinf*Ninv*Hinf');
                % A = Hinf'*((Hinf*Ninv*Hinf')\Hinf);
                A_mat = Hinf'*(W\(W'\Hinf));
                C = Ninv*A_mat*Ninv;
                nlZ_v = -nu'*C*nu/2;
            end
        else
            nlZ_v = 0;
        end

        % Compute determinants
        nldet = sum(log(diag(L)));                  % First component
        if precany_flag                             % Precise priors
            nldet = nldet + sum(log(BB(prec_idx)))/2;
            Tprec_idx = isfinite(BB(BB > 0));
            nldet = nldet + log(det(T_plus(Tprec_idx,Tprec_idx)))/2;
        end
        if vagueany_flag
            nldet = nldet + sum(log(diag(W)));
        end

        nlZ = nlZ_1 + nlZ_v + nldet + N*log(2*pi*sl)/2 - sum(inf_idx)*log(2*pi)/2;

    else
        % Compute negative log marginal likelihood
        nlZ = (y-m)'*alpha/2 + sum(log(diag(L))) + N*log(2*pi*sl)/2;
    end

    if outwarp_flag     % Jacobian correction for output warping
        nlZ = nlZ - sum(log(abs(dwarp_dt)));
    end

    if compute_nlZ_grad
        % Compute gradient of negative log marginal likelihood

        dnlZ = zeros(Nhyp,1);    % allocate space for derivatives

        if gp.intmeanfun > 0
            % Gradient with integrated basis functions
            Kinv = L\(L'\eye(N))/sl;
            if ~precany_flag; Ninv = Kinv; end
            chi = Ninv*nu;
            Q = Kinv - chi*chi';
            if vagueany_flag
                phi = Ninv*A_mat*chi;
                Q = Q - phi*phi' + 2*chi*phi';
                if ~precany_flag
                    Q = Q - Kinv*A_mat*Kinv;
                else
                    Q = Q - Ninv*A_mat*Ninv;
                end
            else
                phi = 0;
            end
            if precany_flag
                Q = Q - Kinv*H(prec_idx,:)'*(T_plus(Tprec_idx,Tprec_idx)\(H(prec_idx,:)*Kinv));
            end
        else
            Q = L\(L'\eye(N))/sl - alpha*alpha';
        end

        if gp.covfun(1) == 1
            for i = 1:D                             % Grad of cov length scales
                K_temp = K_mat .* sq_dist(gp.X(:,i)'/ell(i));
                dnlZ(i) = sum(sum(Q.*K_temp))/2;
            end
            dnlZ(D+1) = sum(sum(Q.*(2*K_mat)))/2;   % Grad of cov output scale
        else
            for i = 1:Ncov                          % Grad of cov hyperparameters
                dnlZ(i) = sum(sum(Q.*dK_mat(:,:,i)))/2;
            end
        end

        % Gradient of GP likelihood
        if isscalar(sn2)
            trQ = trace(Q);
            for i = 1:Nnoise; dnlZ(Ncov+i) = 0.5*sn2_mult*dsn2(i)*trQ; end
        else
            dgQ = diag(Q);
            for i = 1:Nnoise; dnlZ(Ncov+i) = 0.5*sn2_mult*sum(dsn2(:,i).*dgQ); end
            if outwarp_flag
                error('Input-dependent noise not supported with output warping yet.');
            end
        end

        % Gradient of mean function
        if Nmean > 0
            if gp.intmeanfun > 0
                % Mean function gradient with integrated basis functions
                dnlZ(Ncov+Nnoise+(1:Nmean)) = -dm'*(chi - phi);
            else
                dnlZ(Ncov+Nnoise+(1:Nmean)) = -dm'*alpha;
            end
        end

        % Gradient of output warping function
        if outwarp_flag && Noutwarp > 0
            if gp.intmeanfun > 0
                error('Integrated basis functions are not supported with output warping yet.');
            end
            for i = 1:Noutwarp
                dnlZ(Ncov+Nnoise+Nmean+i) = dwarp_dtheta(:,i)'*alpha ...
                    - sum(d2warp_dthetadt(:,i)./dwarp_dt);
            end
        end

    end
end

%% Output posterior struct if requested
if nargout > 2
    post.hyp = hyp;
    post.alpha = alpha;
    post.sW = ones(N,1)./sqrt(min(sn2)*sn2_mult);   % sqrt of noise precision vector
    post.L = pL;
    post.sn2_mult = sn2_mult;
    post.Lchol = Lchol;
    if gp.intmeanfun > 0
        post.intmean.HKinv = H*(L\(L'\eye(N))/sl);
         % Inverse reduced T (only positive variances)
        post.intmean.Tplusinv = T_chol \ (T_chol' \eye(size(T_chol)));
        post.intmean.betabar = betabar;
    end
end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function l = l1dd_inverse ( n, h )

%*****************************************************************************80
%
%% L1DD_INVERSE stores the inverse of the 1D DD Laplacian.
%
%  Discussion:
%
%    The N grid points are assumed to be evenly spaced by H.
%
%    For N = 5, the discrete Laplacian with Dirichlet boundary conditions
%    at both ends of [0,6] has the matrix form L:
%
%       2 -1  0  0  0
%      -1  2 -1  0  0
%       0 -1  2 -1  0
%       0  0 -1  2 -1
%       0  0  0 -1  2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points.
%    N must be at least 3.
%
%    Input, real H, the spacing between points.
%
%    Output, real L(N,N), the inverse of the Laplacian matrix.
%
  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'L1DD_INVERSE - Fatal error!\n' );
    fprintf ( 1, '  N < 3.\n' );
    error ( 'L1DD_INVERSE - Fatal error!' );
  end

  l = zeros ( n, n );

  for j = 1 : n
    for i = 1 : n
      l(i,j) = min ( i, j ) * ( n + 1 - max ( i, j ) ) * h * h / ( n + 1 );
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function fx = p10_fun ( n, x )

%*****************************************************************************80
%
%% P10_FUN evaluates the integrand for problem 10.
%
%  Interval:
%
%    0 <= x <= 1
%
%  Integrand:
%
%    1 / ( 1 + x )
%
%  Antiderivative:
%
%    ln ( 1 + x )
%
%  Exact Integral:
%
%    ln ( 2 )
%
%  Approximate Integral (25 digits):
%
%    0.6931471805599453094172321...
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    David Kahaner,
%    Comparison of Numerical Quadrature Formulas,
%    in Mathematical Software, edited by John R Rice,
%    Academic Press, 1971.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  fx = 1.0 ./ ( 1.0 + x );

  return
end


text:
filename='CantileverBeam_Triangle_Linear_Fine';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1];
constraint = {'volumeConstraint'};
optimizer = 'SLERP';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.4;
optimality_final =1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
classdef StokesComputer < handle

    properties (Access = public)
        computation
    end

    properties (Access = private)
        testName
    end

    methods (Access = public)
        function obj = StokesComputer(cParams)
            obj.testName = cParams.testName;
        end

        function compute(obj)
            a.fileName = obj.testName;
            s = StokesDataContainer(a);
            femSolver = FEM.create(s);
            femSolver.computeVariables;
            obj.computation = femSolver;
        end
    end

end


text:
% op_B0Correction.m
% Jamie Near, McGill University 2021.
%E dits from Brenden Kadota, 2021.
%
% USAGE:
% [out, phaseMap, freqMap]=op_B0Correction(in, (optional) phaseMap);
%
% DESCRIPTION:
% Corrects for slight B0 shifts by aligning all voxels to the voxel with
% the highest water peak. (Possibly not correct output as highest peak may
% no have the correct ppm shifts)
%
% INPUTS:
% in   = MRSI struct
% phaseMap = phaseMap matrix used to apply corrections
%
% OUTPUTS:
% B0Map        = 2D array of shift intensity at each specific voxel
% phaseMap     = 2D array of phase shift intensity at each specific voxel
% freqMap      = 2D array of frequency shift intensity at each specific voxel

function [MRSIStruct, phaseMap, freqMap] = op_CSIB0Correction(MRSIStruct, phaseMap, ...
        freqMap, plottingArguments)
    arguments
        MRSIStruct (1, 1) struct
        phaseMap double = []
        freqMap double = []
        plottingArguments.isPlot = false;
    end
    %only correct if coils have be combined
    if(getFlags(MRSIStruct, 'addedrcvrs') == 0)
        error("MRSIStruct Error. Add coils first in MRSIStruct");
    end
    if(getFlags(MRSIStruct, 'spatialft') == 0)
        error("Please fourier transform along the spatial dimension");
    end
    if(getFlags(MRSIStruct, 'spectralft') == 0)
        error("Please fourier transform along the spectral dimension");
    end

    %find the coordinate of the highest intensity
    [x, y, a] = findMaxCoord(MRSIStruct);

    %create twix object from coordinate with highest intensity
    referenceMRS = op_CSItoMRS(MRSIStruct, x, y, "averageIndex", a);
    referenceMRS = setFlags(referenceMRS, 'averaged', true);
    referenceMRS.flags.isISIS = 0;

    if(~isempty(phaseMap) && ~isempty(freqMap))
        [MRSIStruct, prevPermute, prevSize] = reshapeDimensions(MRSIStruct, {'t', 'y', 'x'});
        MRSIStruct = applyFreqMap(MRSIStruct, freqMap);
        MRSIStruct = applyPhaseMap(MRSIStruct, phaseMap);
        MRSIStruct = reshapeBack(MRSIStruct, prevPermute, prevSize);

    else

        [MRSIStruct, prevPermute, prevSize] = reshapeDimensions(MRSIStruct, {'t', 'y', 'x'});

        phaseMap = zeros(getSizeFromDimensions(MRSIStruct, {'y', 'x', 'extras'}));
        freqMap = zeros(getSizeFromDimensions(MRSIStruct, {'y', 'x', 'extras'}));

        data = getData(MRSIStruct);
        for e = 1:getSizeFromDimensions(MRSIStruct, {'extras'})
            for y = 1:getSizeFromDimensions(MRSIStruct, {'y'})
                for x = 1:getSizeFromDimensions(MRSIStruct, {'x'})
                    %get twix object at x and y coordinate
                    alignMRS = op_CSItoMRS(MRSIStruct, x, y, 'Extra', e);
                    alignMRS = setFlags(alignMRS, 'averaged', true);
                    alignMRS.flags.isISIS = 0;
                    %aligns twix object to the maximum intensity twix object
                    [alignMRS, phase, freq] = op_alignScans(alignMRS, referenceMRS, referenceMRS.t(end));

                    timeDimension = getDimension(alignMRS, 't');
                    %adds the aligned twix to the specs of MRSI object.
                    data(:, y, x, e) = flip(alignMRS.specs, timeDimension);
                    %add phase and frequency to maps
                    phaseMap(y, x, e) = phase;
                    freqMap(y, x, e) = freq;
                end
            end
        end
        MRSIStruct = setData(MRSIStruct, data);
        MRSIStruct = reshapeBack(MRSIStruct, prevPermute, prevSize);

        %plot maps
        if(plottingArguments.isPlot)
            plotFreqAndPhaseMap(phaseMap, MRSIStruct, freqMap);
        end
        MRSIStruct = setFlags(MRSIStruct, 'phasecorrected', 1);
        MRSIStruct = setFlags(MRSIStruct, 'freqcorrected', 1);
    end
end

%finds the (x,y) coordinate with the highest intensity spectrum.
function [x,y,a] = findMaxCoord(in)
    [~, i] = max(getData(in),[], 'all', 'linear');
    out = cell(size(in.sz));
    [out{:}] = ind2sub(in.sz, i);
    x = out{in.dims.x};
    y = out{in.dims.y};
    if(in.dims.averages)
        a = out{in.dims.averages};
    else
        a = 1;
    end

end

function plotFreqAndPhaseMap(phaseMap, MRSIStruct, freqMap)
    for i = 1:size(phaseMap,3)
        figure;
        subplot(2,1,1);
        imagesc(MRSIStruct.coordinates.x, MRSIStruct.coordinates.y, phaseMap(:,:,i)')
        title("phaseMap")
        subplot(2,1,2);
        imagesc(MRSIStruct.coordinates.x, MRSIStruct.coordinates.y, freqMap(:,:,i)')
        title("freqMap")
    end
end

function MRSIStruct = applyFreqMap(MRSIStruct, freqMap)
    %applying the freqMap map if added as an argument
    for e = 1:getSizeFromDimensions(MRSIStruct, {'extras'})
        for x = 1:getSizeFromDimensions(MRSIStruct, {'x'})
            for y = 1:getSizeFromDimensions(MRSIStruct, {'y'})
                %create twix object from the coordinates
                alignMRS = op_CSItoMRS(MRSIStruct, x, y, 'extraIndex', e);
                %align twix object with op_freqshift to the equivalent
                %freqmap coordinate.
                alignMRS = op_freqshift(alignMRS, freqMap(x, y));
                %add the specs to MRSI object
                MRSIStruct.data(:, x, y, e) = alignMRS.specs;
            end
        end
    end
end

function MRSIStruct = applyPhaseMap(MRSIStruct, phaseMap)
    %applying the freqMap map if added as an argument
    for e = 1:getSizeFromDimensions(MRSIStruct, {'extras'})
        for x = 1:getSizeFromDimensions(MRSIStruct, {'x'})
            for y = 1:getSizeFromDimensions(MRSIStruct, {'y'})
                %create twix object from the coordinates
                alignMRS = op_CSItoMRS(MRSIStruct, x, y, 'extraIndex', e);
                %align twix object with op_freqshift to the equivalent
                %freqmap coordinate.
                alignMRS = op_addphase(alignMRS, phaseMap(x, y));
                %add the specs to MRSI object
                MRSIStruct.data(:, x, y, e) = alignMRS.specs;
            end
        end
    end
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
%% Copyright (C) 2016 Lagu
%% Copyright (C) 2017-2019 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @deftypemethod @@sym {@var{L} =} chol (@var{A})
%% Cholesky factorization of symbolic symmetric matrix.
%%
%% Returns a lower-triangular matrix @var{L}, such that @code{L*L'}
%% is matrix @var{A}.  The matrix @var{A} must be symmetric
%% positive-definite.  Example:
%% @example
%% @group
%% A = sym([1 2 4; 2 13 23; 4 23 43])
%%   @result{} A = (sym 3×3 matrix)
%%
%%       ⎡1  2   4 ⎤
%%       ⎢         ⎥
%%       ⎢2  13  23⎥
%%       ⎢         ⎥
%%       ⎣4  23  43⎦
%%
%% L = chol(A)
%%   @result{} L = (sym 3×3 matrix)
%%
%%       ⎡1  0  0 ⎤
%%       ⎢        ⎥
%%       ⎢2  3  0 ⎥
%%       ⎢        ⎥
%%       ⎣4  5  √2⎦
%%
%% L*L'
%%   @result{} (sym 3×3 matrix)
%%
%%       ⎡1  2   4 ⎤
%%       ⎢         ⎥
%%       ⎢2  13  23⎥
%%       ⎢         ⎥
%%       ⎣4  23  43⎦
%% @end group
%% @end example
%%
%% @seealso{chol, @@sym/qr, @@sym/lu}
%% @end deftypemethod

function y = chol(x)
  if (nargin == 2)
    error('Operation not supported yet.');
  elseif (nargin > 2)
    print_usage ();
  end
  y = pycall_sympy__ ('return _ins[0].cholesky(),', x);
end

%!error <must be> chol (sym ([1 2; 3 4]));
%!error <must be square> chol (sym ([1 2; 3 4; 5 6]));

%!test
%! A = chol(hilb(sym(2)));
%! B = [[1 0]; sym(1)/2 sqrt(sym(3))/6];
%! assert( isequal( A, B ))

repo: octsympy


text:
% std_selectdesign() - select an existing STUDY design.
%                      Use std_makedesign() to add a new STUDY.design.
%
% Usage:
%   >> [STUDY] = std_selectdesign(STUDY, ALLEEG, designind);
%
% Inputs:
%   STUDY     - EEGLAB STUDY structure
%   STUDY     - EEGLAB ALLEEG structure
%   designind - desired (existing) design index
%
% Outputs:
%   STUDY     - EEGLAB STUDY structure with currentdesign set to the input design index.
%
% Author: Arnaud Delorme, Institute for Neural Computation, UCSD, 2010-

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

function STUDY = std_selectdesign(STUDY, ALLEEG, designind);

if nargin < 3
    help std_selectdesign;
    return;
end;

if designind < 1 || designind > length(STUDY.design) || isempty(STUDY.design(designind).name)
    disp('Cannot select an empty STUDY.design');
    return;
end;
STUDY.currentdesign = designind;
STUDY = std_rmalldatafields( STUDY );

% remake setinds and allinds
% --------------------------
STUDY = std_changroup(STUDY, ALLEEG, [], 'on'); % with interpolation HAVE TO FIX THAT

% update the component indices
% ----------------------------
STUDY.cluster(1).setinds = {};
STUDY.cluster(1).allinds = {};
for index = 1:length(STUDY.cluster)
    STUDY.cluster(index) = std_setcomps2cell(STUDY, index);
end;

author: PatternRecognition

repo: OpenBMI


text:
function [a,c] = myhungarian(A)
%MYHUNGARIAN

ndig = 3;

A = round(A*10^ndig);

[a,c] = myhung(A);
c = c *10^(-ndig);

author: kristinbranson

repo: JAABA


text:
function data=extractdatapt(data,t,offset)
% Extract segements of spike times between t(1) and t(2)
% Usage: data=extractdatapt(data,t,offset)
%
% Input:
% data: structural array of spike times for each channel/trial or a single
%       array of spike times
% t   : time as a 2d vector [t(1) t(2)]
% offset: 0/1 - if 1, store the spike times relative to start of window i.e. t(1)
%         if 0, don't reset the times. Default 0.
% Note that all times can be in arbitrary units. But the units have to be
% consistent. So, if E is in secs, win, t have to be in secs, and Fs has to
% be Hz. If E is in samples, so are win and t, and Fs=1. In case of spike
% times, the units have to be consistent with the units of data as well.
%
% Output:
% data: spike times between t(1) and t(2)
if nargin < 2; error('Need data and times'); end;
if t(1) < 0 || t(2)<=t(1);
    error('times cannot be negative and t(2) has to greater than t(1)');
end;
if nargin < 3 || isempty(offset); offset=0; end;
if isstruct(data);
    C=length(data);
elseif min(size(data))~=1;
    error('Can only accept single vector data unless it is a struct array');
else
    C=1;
    data=change_row_to_column(data);
end;
%fnames=fieldnames(data);
d2(1:C)=struct('times',[]);
for c=1:C,
    if isstruct(data)
       fnames=fieldnames(data);
       eval(['dtmp=data(c).' fnames{1} ';'])
    else
       dtmp=data(:);
    end
%     eval(['dtmp=data(c).' fnames{1} ';' ])
    sp=dtmp(dtmp>=t(1) & dtmp<t(2));
    if offset==1; d2(c).times=sp-t(1);
    else d2(c).times=sp;end
end;
data=d2;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%% Calculate mean angular error between source and target images
%
% Copyright (c) 2018-present, Mahmoud Afifi
% York University, Canada
% mafifi@eecs.yorku.ca | m.3afifi@gmail.com
%
% This source code is licensed under the license found in the
% LICENSE file in the root directory of this source tree.
% All rights reserved.
%
% Please cite the following work if this program is used:
% Mahmoud Afifi, Brian Price, Scott Cohen, and Michael S. Brown,
% "When color constancy goes wrong: Correcting improperly white-balanced
% images", CVPR 2019.
%
% Input:
%   -source: image A
%   -target: image B
%   -color_chart_area: If there is a color chart in the image (that is
%   masked out from both images, this variable represents the number of
%   pixels of the color chart.
%
% Output:
%   -f: the mean angular error between image A and image B.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
function f = calc_mae(source, target,color_chart_area)

if size(source,1)<=color_chart_area
    error('Color chart area should be less than the image area');
end
source=double(source);
target=double(target);
target_norm = sqrt(sum(target.^2,2));
source_mapped_norm = sqrt(sum(source.^2,2));
angles=dot(source,target,2)./(source_mapped_norm.*target_norm);
angles(angles>1)=1;
f=acosd(angles);
f(isnan(f))=0;
f=sum(f)/(size(source,1)-color_chart_area);
end

author: mahmoudnafifi


text:
function [D]=patchEdgeLengths(F,V)

% function [D]=patchEdgeLengths(F,V)
% -----------------------------------------------------------------------
% Computes the edge lengths (D) for the patch data specified by the faces
% (F) and vertices (V) arrays. If size(F,2)>2 it is assumed that F indeed
% represents faces. If however size(F,2)==2 it is instead assumed that F is
% an array representing edges. As such it skips the computation of the
% edges array. The edges array used is non-unique by default. See the
% |patchEdges| function for more details if the lengths of a unique set of
% edges is desired.
%
%
% See also: |patchEdges|
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2014/03/17
%------------------------------------------------------------------------

%%

%Derive edge array
if size(F,2)>2 %The input is assumed to represent faces hence an edge array is derived
    E=patchEdges(F);
else %It is assumed that the input array represents an edges array
    E=F;
end

%Derive edge vertex arrays
V_E1=V(E(:,1),:);
V_E2=V(E(:,2),:);

%Derive difference vectors
VD=(V_E1-V_E2);

%Compute the edge lengths
D=sqrt(sum(VD.^2,2));

%%
% _*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) 2017  Kevin Mattheus Moerman
%
% 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: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function listString = getListString(list)
% Produce a string with list values
   listString = sprintf('%g ', list(:)');
   listString = ['[ ' listString ']'];

author: VisLab

repo: EEG-Clean-Tools


text:
function displayBarResults_MetricBased_HN_tempFig4(pathExperiments,fSetNames,nameOutcome,metrics,titlePart)
% - fSetNames: Cell of strings specifying the names of all feature sets
% - nameOutcome: String, specifying the name of a given specific outcome
% - titlePart: optional

startpath = pwd;

nExp = 1;

colors = {'red','blue','green'};

nMetrics = numel(metrics);
if nargin < 5
    prefix = 'testResults';
else
    prefix = ['testResults',titlePart];
end

nFset = numel(fSetNames);
data = cell(1,nMetrics);
for m = 1:nMetrics
    data{m} = zeros(nExp,nFset);
end

% Obtaining data
cd(pathExperiments), pathExperiment = pwd;
for exp = 1:nExp
    for f = 1:nFset
        results = load([prefix,'_',fSetNames{f},'_',nameOutcome]); results = struct2cell(results); results = results{1};
        for m = 1:nMetrics
            data{m}(exp,f) = results.(metrics{m});
        end
    end
end

% Producing bar plot
% scrsz = get(groot,'ScreenSize');
% h = figure('Position',[1 7*scrsz(4)/10 19*scrsz(3)/40 7*scrsz(4)/10]);
ind = strfind(nameOutcome,'Death');
if ~isempty(ind)
    nameOutcome(ind:ind+4) = [];
    nameOutcome = ['Survival',nameOutcome];
end
if nargin < 5
    titleName = nameOutcome;
else
    titleName = [nameOutcome,' -- ',titlePart];
end
if exp > 1
    y = zeros(nMetrics,nFset);
    err = zeros(nMetrics,nFset);
    stdData = data;
    for m = 1:nMetrics
        data{m} = mean(data{m});
        stdData{m} = std(data{m});
        y(m,:) = data{m}(1,:);
        err(m,:) = stdData{m}(1,:);
    end
    model_series = y; model_error = err;
    bar(model_series,1);
    set(gca,'YGrid','on')
    set(gca,'GridLineStyle','-')
    set(gca,'XTicklabel','Modelo1|Modelo2|Modelo3')
    lh = legend('Serie1','Serie2','Serie3');
    set(lh,'Location','BestOutside','Orientation','horizontal')
    hold on;
    numgroups = size(model_series, 1);
    numbars = size(model_series, 2);
    groupwidth = min(0.8, numbars/(numbars+1.5));
    for i = 1:numbars
          % Based on barweb.m by Bolu Ajiboye from MATLAB File Exchange
          x = (1:numgroups) - groupwidth/2 + (2*i-1) * groupwidth / (2*numbars);  % Aligning error bar with individual bar
          errorbar(x, model_series(:,i), model_error(:,i), 'k', 'linestyle', 'none','LineWidth',3);
    end
else
    y = zeros(nMetrics,nFset);
    for m = 1:nMetrics
        y(m,:) = data{m}(1,:);
    end
    if nMetrics > 1
        b = bar(y,1);
        for c = 1:nFset
            b(c).FaceColor = colors{c};
        end
    else
        bar([y;zeros(nFset-1,nFset)],1)
    end
end
set(gca,'XTickLabel',metrics)
set(gca,'FontSize',10)
%set(gca,'XTickLabelRotation',45)
axis([0.5 nMetrics+0.5 0 1])
legend(fSetNames)
grid on
%title(titleName,'FontSize',26,'FontWeight','bold')
% annotation(h,'line',[0.13 0.905],[0.515673981191221 0.515673981191221],...
%     'Color',[0.850980401039124 0.325490206480026 0.0980392172932625],...
%     'LineWidth',5,'LineStyle',':');

hold on
x = 0:0.01:nMetrics+1;
y = ones(1,numel(x))*0.5;
plot(x,y,'.m')
hold off

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [output] = ft_transform_geometry(transform, input)

% FT_TRANSFORM_GEOMETRY applies a homogeneous coordinate transformation to a
% structure with geometric information, for example a volume conduction model for the
% head, gradiometer of electrode structure containing EEG or MEG sensor positions and
% MEG coil orientations, a head shape or a source model.
%
% Use as
%   [output] = ft_transform_geometry(transform, input)
% where the transform should be a 4x4 homogeneous transformation matrix and the input
% data structure can be any of the FieldTrip data structures that describes
% geometrical data.
%
% The units of the transformation matrix must be the same as the units in which the
% geometric object is expressed.
%
% The type of geometric object constrains the type of allowed transformations.
%
% For sensor arrays:
% If the input is an MEG gradiometer array, only a rigid-body translation plus
% rotation are allowed. If the input is an EEG electrode or fNIRS optodes array,
% global rescaling and individual axis rescaling is also allowed.
%
% For volume conduction models:
% If the input is a volume conductor model of the following type:
%   localspheres model
%   singleshell model with the spherical harmonic coefficients already computed
%   BEM model with system matrix already computed
%   FEM model with volumetric elements
% only a rigid-body translation plus rotation are allowed.
%
% If the input is a volume conductor model of the following type:
%   BEM model with the system matrix not yet computed
%   singleshell model with the spherical harmonic coefficients not yet computed
% rotation, translation, global rescaling and individual axis rescaling is allowed.
%
% If the input is a volume conductor model of the following type:
%   single sphere
%   concentric spheres
% rotation, translation and global rescaling is allowed.
%
% For source models, either defined as a 3D regular grid, a 2D mesh or unstructred
% point cloud, rotation, translation, global rescaling and individual axis rescaling
% is allowed.
%
% For anatomical MRIs and functional volumetric data, rotation, translation, global
% rescaling and individual axis rescaling are allowed.
%
% See also FT_WARP_APPLY, FT_HEADCOORDINATES, FT_SCALINGFACTOR

% Copyright (C) 2011-2022, Jan-Mathijs Schoffelen and 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: ft_transform_geometry.m$

% determine the rotation matrix
rotation = eye(4);
rotation(1:3,1:3) = transform(1:3,1:3);

if any(abs(transform(4,:)-[0 0 0 1])>100*eps)
  ft_error('invalid transformation matrix');
end

% check whether the transformation includes a scaling operation
s = svd(rotation);
if abs(abs(det(rotation))-1)>1e-4
  % the transformation increases or decreases the overall volume, allow for
  % some numerical imprecision
  if any(abs(s./s(1)-1)>1e-3)
    % axes scale differently
    globalrescale = false;
    axesrescale   = true;
  else
    % global rescaling
    globalrescale = true;
    axesrescale   = true;
  end
else
  globalrescale = false;
  axesrescale   = false;
end

% check whether the input data combines well with the requested
% transformation
dtype = ft_datatype(input);
switch dtype
  case 'grad'
    if globalrescale || axesrescale, ft_error('only a rigid body transformation without rescaling is allowed'); end
  case 'mesh'
    if isstruct(input) && length(input)>1
      % there are multiple meshes
      for i=1:length(input)
        output(i) = ft_transform_geometry(transform, input(i));
      end
      return
    end
  otherwise
    % could be a volume conductor model with constrained transformation possibilities
    if ft_headmodeltype(input, 'singleshell') && isfield(input, 'forwpar') && (globalrescale || axesrescale)
      ft_error('only a rigid body transformation without rescaling is allowed');
    end
    if ft_headmodeltype(input, 'bem') && isfield(input, 'mat') && (globalrescale || axesrescale)
      ft_error('only a rigid body transformation without rescaling is allowed');
    end
    if ft_headmodeltype(input, 'localspheres') && isfield(input, 'mat') && (globalrescale || axesrescale)
      ft_error('only a rigid body transformation without rescaling is allowed');
    end
    if (isfield(input, 'tetra') || isfield(input, 'hex')) && (globalrescale || axesrescale)
      ft_error('only a rigid body transformation without rescaling is allowed');
    end
end

% tfields must be rotated, translated and scaled
% rfields must only be rotated
% mfields must be simply multipied
% recfields must be recursed into
tfields   = {'pos' 'pnt' 'o' 'coilpos' 'elecpos' 'optopos' 'chanpos' 'chanposold' 'nas' 'lpa' 'rpa' 'zpoint'}; % apply rotation plus translation
rfields   = {'ori' 'nrm'     'coilori' 'elecori' 'optoori' 'chanori' 'chanoriold' 'mom'                     }; % only apply rotation
mfields   = {'transform'};                % plain matrix multiplication
recfields = {'fid' 'bnd' 'orig' 'dip'};   % recurse into these fields
% the field 'r' is not included here, because it applies to a volume
% conductor model, and scaling is not allowed, so r will not change.

fnames = fieldnames(input);
for k = 1:numel(fnames)
  % name = sprintf('%s.%s', inputname(2), fnames{k});
  if ~isempty(input.(fnames{k}))
    if any(strcmp(fnames{k}, tfields))
      % ft_info('applying transformation to %s', name);
      input.(fnames{k}) = apply(transform, input.(fnames{k}));
    elseif any(strcmp(fnames{k}, rfields))
      % ft_info('applying rotation to %s', name);
      input.(fnames{k}) = apply(rotation, input.(fnames{k}));
    elseif any(strcmp(fnames{k}, mfields))
      % ft_info('applying multiplication to %s', name);
      input.(fnames{k}) = transform*input.(fnames{k});
    elseif any(strcmp(fnames{k}, recfields))
      for j = 1:numel(input.(fnames{k}))
        % ft_info('recursing into %s', name);
        input.(fnames{k})(j) = ft_transform_geometry(transform, input.(fnames{k})(j));
      end
    else
      % do nothing
    end
  end
end
output = input;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that applies the homogeneous transformation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [new] = apply(transform, old)

[m, n] = size(old);
if m~=3 && n==3
  % each row is one position
  old(:,4) = 1;
  new = old * transform';
  new = new(:,1:3);
elseif m==3 && n~=3
  % each column is one position
  old(4,:) = 1;
  new = transform * old;
  new = new(1:3,:);
else
  % assume that each row is one position
  old(:,4) = 1;
  new = old * transform';
  new = new(:,1:3);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef PTKCoronalAnalysis < PTKPlugin
    % PTKCoronalAnalysis. Plugin for performing analysis of density using bins
    % along the anterior-posterior axis
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %     PTKAxialAnalysis divides the cranial-caudal axis into bins and
    %     performs analysis of the tissue density, air/tissue fraction and
    %     emphysema percentaein each bin.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Coronal<br>analysis'
        ToolTip = 'Performs density analysis in bins along the anterior-posterior axis'
        Category = 'Analysis'

        Context = PTKContextSet.Any
        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'DoNothing'
        HidePluginInDisplay = true
        FlattenPreviewImage = false
        PTKVersion = '2'
        ButtonWidth = 6
        ButtonHeight = 1
        GeneratePreview = false
    end

    methods (Static)
        function results = RunPlugin(dataset, context, reporting)

            % Get the density image
            roi = dataset.GetResult('PTKLungROI', PTKContext.LungROI);
            if ~roi.IsCT
                reporting.ShowMessage('PTKCoronalAnalysis:NotCTImage', 'Cannot perform density analysis as this is not a CT image');
                return;
            end

            % Get a mask for the current region to analyse
            context_mask = dataset.GetTemplateMask(context);

            % Create a region mask excluding the airways
            context_no_airways = dataset.GetResult('PTKGetMaskForContextExcludingAirways', context);

            % Divide the lung into bins along the cranial-caudal axis
            bins = dataset.GetResult('PTKDivideLungsIntoCoronalBins', PTKContext.Lungs);

            results = PTKMultipleRegionAnalysis(bins, roi, context_mask, context_no_airways, 'Coronal distance from lung edge (mm)', reporting);
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function test96()
% function test96()
% Low-rank Euclidena distance matrix completion
%
% We use the tuned geometry, symfixedrankNewYYquotientfactory to solve the
% EDMCP. This test file is different from test file 'test97'
%

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, Dec. 30, 2012.
% Contributors:
% Change log:

clear all; clc; close all;

m = 500;
r = 5;
Yo = randn(m,r); % True embedding
C = squareform(pdist(Yo)'.^2); % True distances

% Create the problem structure
% quotient YY (tuned for least square problems) geometry
problem.M = symfixedrankNewYYquotientfactory(m,  r);

df = problem.M.dim();
p = 5*df/(m*m);
symm = @(M) .5*(M+M');

% mask = symm(rand(m, m)) <= p;
mask = spones(sprandsym(m, p));

prop_known = sum(sum(mask == 1))/(m*m);
fprintf('Low-rank EDM completion... \n');
fprintf('Fraction of entries given: %f \n', full(prop_known));

J = eye(m) - (1/m)*ones(m,1)*ones(m,1)'; % Centering matrix

problem.cost = @cost;
    function f = cost(X)
        YY = X.Y*X.Y';
        m = size(X.Y, 1);
        KvX = diag(YY)*ones(1,m)  + ones(m,1)*diag(YY)' - 2*YY;
        f = 0.5*(norm(mask.*(KvX - C), 'fro')^2);
    end

problem.grad = @grad;
    function g = grad(X)
        Y = X.Y;
        YY = X.Y*X.Y';
        m = size(X.Y, 1);
        KvX = diag(YY)*ones(1,m)  + ones(m,1)*diag(YY)' - 2*YY;
        J= eye(m) - (1/m)*ones(m,1)*ones(m,1)' ;
        mat = mask.*(KvX - C);
        R = J*(diag(mat*ones(m,1)) - mat)*J;
        R = 4*R;
        r = size(Y, 2);
        YtY = Y'*Y;
        invYtY = eye(r) / YtY;

        g = struct('Y', R*Y*invYtY);
    end

problem.hess = @hess;
    function Hess = hess(X, eta)

        Y = X.Y;
        YY = X.Y*X.Y';
        m = size(X.Y, 1);
        KvX = diag(YY)*ones(1,m)  + ones(m,1)*diag(YY)' - 2*YY;
        J = eye(m) - (1/m)*ones(m,1)*ones(m,1)' ;
        mat = mask.*(KvX - C);
        R = J*(diag(mat*ones(m,1)) - mat)*J;
        R = 4*R;

        Xdot = eta.Y*Y' + Y*eta.Y';
        XdotV = Xdot;
        KvXdot = diag(XdotV)*ones(1, m) + ones(m, 1)*diag(XdotV)' - 2*XdotV;
        mat2 = mask.*KvXdot;
        KvstarMatdot = J*(diag(mat2*ones(m,1)) - mat2)*J;
        KvstarMatdot = 4*KvstarMatdot;

        r = size(Y, 2);
        YtY = Y'*Y;
        invYtY = eye(r) / YtY;

        Hess.Y = R*eta.Y*invYtY  +  KvstarMatdot*Y*invYtY;
        Hess.Y = Hess.Y - 2*R*Y*(invYtY * symm(eta.Y'*X.Y) * invYtY);

        gradY = R*Y*invYtY;

        % I still need a correction factor for the non-constant metric
        Hess.Y = Hess.Y + gradY*symm(eta.Y'*X.Y)*invYtY + eta.Y*symm(gradY'*X.Y)*invYtY - X.Y*symm(eta.Y'*gradY)*invYtY;

        Hess = problem.M.proj(X, Hess);
    end

% % Check numerically whether gradient and Hessian are correct
%     checkgradient(problem);
%     drawnow;
%     pause;
%     checkhessian(problem);
%     drawnow;
%     pause;

% Initialization
[U, S, ~ ] = svds(-0.5*J*(mask.*C)*J, r); % This initialization is based on the strain error minimization of EDM
Y0 = U*(S.^0.5);
X0 = struct('Y', Y0);

% Options (not mandatory)
options.maxiter = inf;
options.maxinner = 30;
options.maxtime = 120;
options.tolgradnorm = 1e-9;
options.Delta_bar = m * r;
options.Delta0 = options.Delta_bar / 8;

% Pick an algorithm to solve the problem
   [Xopt, costopt, info] = trustregions(problem, X0, options);
%  [Xopt costopt info] = steepestdescent(problem, X0, options);
% [Xopt costopt info] = conjugategradient(problem, X0, options);

end

author: NicolasBoumal

repo: manopt


text:
function [r,s] = xytors(x,y)

% function [r,s] = xytors(x, y)
% Purpose : From (x,y) in equilateral triangle to (r,s) coordinates in standard triangle

L1 = (sqrt(3.0)*y+1.0)/3.0;
L2 = (-3.0*x - sqrt(3.0)*y + 2.0)/6.0;
L3 = ( 3.0*x - sqrt(3.0)*y + 2.0)/6.0;

r = -L2 + L3 - L1; s = -L2 - L3 + L1;
return;

repo: Aero-matlab


text:
function [init_sz, final_sz] = ideal_size(net, max_sz)

final_sz = forward(net, max_sz);
init_sz = backward(net, final_sz);
while ~all(init_sz <= max_sz)
    final_sz = final_sz - 1;
    init_sz = backward(net, final_sz);
end

end

function n = forward(net, n)
    for i = 1:numel(net.layers)
        l = net.layers{i};
        switch l.type
        case 'conv'
            m = [size(l.weights{1}, 1), size(l.weights{1}, 2)];
            n = filter(n, l.pad, m, l.stride);
        case 'pool'
            n = filter(n, l.pad, l.pool, l.stride);
        end
    end
end

function n = backward(net, n)
    for i = numel(net.layers):-1:1
        l = net.layers{i};
        switch l.type
        case 'conv'
            m = [size(l.weights{1}, 1), size(l.weights{1}, 2)];
            n = unfilter(n, l.pad, m, l.stride);
        case 'pool'
            n = unfilter(n, l.pad, l.pool, l.stride);
        end
    end
end

function n = filter(n, pad, m, k)
assert(numel(pad) == 1);
n = floor((n + 2*pad - m) / k) + 1;
end

function n = unfilter(n, pad, m, k)
assert(numel(pad) == 1);
n = k*(n - 1) + m - 2*pad;
end

author: bertinetto

repo: cfnet

save_path: github-repos/MATLAB/bertinetto-cfnet


text:
function f = fuison_base_parts(b1, b2)
% ZCA + l1 norm
% featrue1 = whitening_norm(b1, 1, 0);
% featrue2 = whitening_norm(b2, 1, 0);
% figure;imshow(featrue1);
% figure;imshow(featrue2);

% f_sum = featrue1+featrue2;
% f_sum(f_sum==0) = 1;

% w1 = featrue1./f_sum;
% w2 = featrue2./f_sum;
% figure;imshow(w1);
% figure;imshow(w2);

w1 = 0.5;
w2 = 0.5;
f = w1.*b1+w2.*b2;
% figure;imshow(f);
end

author: thfylsty


text:
function [comb] = mne_combine_xyz(vec)
%
% function [comb] = mne_combine_xyz(vec)
%
% Compute the three Cartesian components of a vector together
%
%
% vec         - Input row or column vector [ x1 y1 z1 ... x_n y_n z_n ]
% comb        - Output vector [x1^2+y1^2+z1^2 ... x_n^2+y_n^2+z_n^2 ]
%

%
%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.1  2006/05/05 03:50:40  msh
%   Added routines to compute L2-norm inverse solutions.
%   Added mne_write_inverse_sol_stc to write them in stc files
%   Several bug fixes in other files
%
%

me = 'MNE:mne_combine_xyz';
if nargin ~= 1
    error(me,'Wrong number of arguments');
end
if size(vec,1) > size(vec,2)
    vec = vec';
end
if size(vec,1) ~= 1 || mod(size(vec,2),3) ~= 0
    error(me,'Input must be a row or a column vector with 3N components');
end

s = mne_block_diag(vec,3);
comb = full(diag(s*s'));
if size(vec,1) > size(vec,2)
    comb = comb';
end

return;


text:
function h = mtimes(f, g)
%*   BALLFUN multiplication.
%   a*F and F*a multiplies the BALLFUN F by the scalar a.
%
% See also TIMES.

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

fIsBallfun = isa(f, 'ballfun');
gIsBallfun = isa(g, 'ballfun');
fIsBallfunv = isa(g, 'ballfunv');

if (fIsBallfun && isnumeric(g))
    h = ballfun(g*f.coeffs, 'coeffs');
elseif (isnumeric(f) && gIsBallfun)
    h = ballfun(f*g.coeffs, 'coeffs');
elseif (fIsBallfun && gIsBallfun)
    h = times(f,g);
elseif fIsBallfunv
    h = mtimes(g, f);
else
    error('BALLFUN:mtimes:unknown', ...
          ['Undefined function ''mtimes'' for input arguments of type ' ...
           '%s and %s.'], class(f), class(g));
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function opt = LegOption
% =======================================================================
% Creates structure of inputs for LegPlot or LegSubplot
% =======================================================================
% opt = LegOption
% -----------------------------------------------------------------------
% OUTPUT
%   - opt.hsize : default 0.01, horizontal size of the box of the legend
%   - opt.vsize : default 0.01, vertical position of the box of the legend
%           (0=bottom, 0.5=middle)
%   - opt.handle : default [], no handle specified
%   - opt.interpreter : default 'Tex', otherwise chamge to 'none'
%   - opt.direction : default 1, horizontal legend
% -----------------------------------------------------------------------
% EXAMPLE
%   plot(1:10);
%   opt = LegOption; opt.interpreter = 'Tex'; opt.hsize = 0.4
%   legpos = LegPlot({'Plot of \alpha'},opt)
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2015. Updated November 2020
% -----------------------------------------------------------------------

opt.hsize        = 0.01;
opt.vsize        = 0.01;
opt.handle       = [];
opt.interpreter  = 'Tex';
opt.direction    = 1;

author: ambropo

repo: VAR-Toolbox

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


text:
function r = isempty(a)
%ISEMPTY      Returns 1 if input is empty, i.e. [], in the sense of Matlab
%
%   r = isempty(a)
%
%There are no empty intervals in the mathematical sense in INTLAB. For details,
%  please see Readme.txt
%

% written  10/15/99     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 10/12/05     S.M. Rump  functionality exchanged with @intval\isempty_
% modified 04/22/09     S.M. Rump  isempty_ removed
%

  if a.complex
    r = isempty(a.mid) | isempty(a.rad);
  else
    r = isempty(a.inf) | isempty(a.sup);
  end

repo: OpenMAS


text:
function [p, s, df] = myStatTest(x1, x2, type)
% [p, s, df] = myStatTest(X1, X2, 'test')
% Statistical test of two samples- gives the confidence level for
% different tests of the two samples X1 and X2.
%
%       'test'   Null hypothesis      Assumption on distributions  Type
%       -----------------------------------------------------------------
%         't'    Means are equal      The variances are equal      t-test
%         'u'    Means are equal      The variances are not equal  t-test
%         'p'    Means are equal      The data are paired          t-test
%         'f'    Variances are equal                               F-test
%         'r'    Pearson Correlation                               Z-test
%
%       TEST(X1, X2) performs a t-test (for equal variances) by default.
%
%       [P S] = TEST(X1, X2, ['test']) returns the confidence level in P
%       and the value of the statistic (t or F) in S.
%
%       Ref: Press et al. 1992. Numerical recipes in C. 14.2, Cambridge.
%
if (nargin == 2)
   type = 't' ;
end
nanInds1 = isnan(x1);
nanInds2 = isnan(x2);

if (type == 't')
  [p s df] = ttest(x1(~nanInds1), x2(~nanInds2)) ;
elseif (type == 'u')
  [p s df] = uttest(x1(~nanInds1), x2(~nanInds2)) ;
elseif (type == 'p')
  nanIndsBoth = nanInds1|nanInds2;
  [p s df] = pttest(x1(~nanIndsBoth), x2(~nanIndsBoth)) ;
elseif (type == 'f')
  [p s] = ftest(x1(~nanInds1), x2(~nanInds2)) ;
elseif (type == 'r')
  sz = max([size(x1);size(x2)]);
  for(k=1:sz(1))
    nanIndsBoth(k,:) = nanInds1(k,:)|nanInds2(k,:);
  end
  n = sum(~nanIndsBoth,1);
  df = n-2;
  if(sz(1)>1&&sz(2)>1)
    % vectorized version- assumes that x1 is Nx1 and x2 is NxM. Returns the
    % 1xM correlation coefficients.
    mn = mean(x1,1);
    sd = std(x1,0,1);
    Z1 = (x1-repmat(mn, [sz(1) 1])) ./ repmat(sd, [sz(1) 1]);
    mn = mean(x2,1);
    sd = std(x2,0,1);
    Z2 = (x2-repmat(mn, [sz(1) 1])) ./ repmat(sd, [sz(1) 1]);
    s = 0;
    for(k=1:sz(1))
        s = s + Z1(k,:).*Z2(k,:);
    end
    s = s./n;
    fZ = 0.5*(log((1+s)./(1-s)));
    p = erfc((abs(fZ).*sqrt(df))./sqrt(2));
  else
    % just use matlab's built-in
    [pearson,p] = corrcoef(x1(~nanIndsBoth), x2(~nanIndsBoth));
    s = pearson(1,2);
    p = p(1,2);
    %     t = s/sqrt((1-s^2)/df);
    %     p = betainc( df / (df + t*t), df/2, 0.5);
    % Fischer's z' method (slightly less conservative):
    %z = 0.5*(log((1+s)/(1-s)));
    %df = length(x1)-3;
    %p = erfc((abs(z)*sqrt(df))/sqrt(2));
  end
elseif(type=='k')
  % Kendall's tau
  % This is essentially a vectorized, simplfied version of the
  % implementation in Matlab's 'corr' function.
  sz = max([size(x1);size(x2)]);
  for(k=1:sz(1))
    nanIndsBoth(k,:) = nanInds1(k,:)|nanInds2(k,:);
  end
  n = sum(~nanIndsBoth,1);
  n2const = n.*(n-1)./2;
  [x1Rank,x1Adj] = tiedrank(x1,1);
  [x2Rank,x2Adj] = tiedrank(x2,1);
  K = zeros(1,sz(2));
  for(k=1:sz(1)-1)
    sgn1 = sign(repmat(x1(k,:),sz(1)-k,1)-x1(k+1:sz(1),:));
    % Zero-out Nans to effectively ignore these data points (the
    % missing values are accounted for in the DF)
    sgn1(isnan(sgn1)) = 0;
    sgn2 = sign(repmat(x2(k,:),sz(1)-k,1)-x2(k+1:sz(1),:));
    sgn2(isnan(sgn2)) = 0;
    for(j=1:size(sgn1,1))
      K = K + sgn1(j,:).*sgn2(j,:);
    end
  end
  %for(k=1:sz(1)-1)
  %  tmp = sum(sign(repmat(x1(k,:),sz(1)-k,1)-x1(k+1:sz(1),:)).*sign(repmat(x2(k,:),sz(1)-k,1)-x2(k+1:sz(1),:)),1);
  %  tmp(nanIndsBoth(k,:)) = 0;
  %  K = K + tmp;
  %end
  s = K ./ sqrt((n2const - x1Adj(1,:)).*(n2const - x2Adj(1,:)));
  ties = ((x1Adj(1,:)>0) | (x2Adj(1,:)>0));
  % Following p-value calc is probably only valid for n>=10
  df = n-2;
  nfact = factorial(n);
  stdK = n2const.*(2.*n+5)./9;
  if(any(ties))
    stdK = stdK - (x1Adj(3,:) + x2Adj(3,:))./18 + x1Adj(2,:).*x2Adj(2,:)./(18.*n2const.*(n-2)) + x1Adj(1,:).*x2Adj(1,:)./n2const;
  end
  stdK = sqrt(stdK);
  tail = 'b';
  switch tail
    case 'b' % 'both or 'ne'
        p = normcdf(-(abs(K)-1) ./ stdK);
        p = min(2*p, 1); % Don't count continuity correction at center twice
    case 'r' % 'right' or 'gt'
        p = normcdf(-(K-1) ./ stdK);
    case 'l' % 'left' or 'lt'
        p = normcdf((K+1) ./ stdK);
  end
end
return;

function [p, t, df] = uttest(d1, d2)
%UTTEST Student's t-test for unequal variances.
%       UTTEST(X1, X2) gives the probability that Student's t
%       calculated on data X1 and X2, sampled from distributions
%       with different variances, is higher than observed, i.e.
%       the "significance" level.  This is used to test whether
%       two sample have significantly different means.
%       [P, T] = UTTEST(X1, X2) gives this probability P and the
%       value of Student's t in T. The smaller P is, the more
%       significant the difference between the means.
%       E.g. if P = 0.05 or 0.01, it is very likely that the
%       two sets are sampled from distributions with different
%       means.
%
%       This works if the samples are drawn from distributions with
%       DIFFERENT VARIANCE. Otherwise, use TTEST.
%
%See also: TTEST, PTTEST.
[l1 c1] = size(d1) ;
n1 = l1 * c1 ;
x1 = reshape(d1, l1 * c1, 1) ;
[l2 c2] = size(d2) ;
n2 = l2 * c2 ;
x2 = reshape(d2, l2 * c2, 1) ;
[a1 v1] = avevar(x1) ;
[a2 v2] = avevar(x2) ;
df = (v1 / n1 + v2 / n2) * (v1 / n1 + v2 / n2) / ...
     ( (v1 / n1) * (v1 / n1) / (n1 - 1) + (v2 / n2) * (v2 / n2) / (n2 -1) ) ;
t = (a1 - a2) / sqrt( v1 / n1 + v2 / n2 ) ;
p = betainc( df / (df + t*t), df/2, 0.5) ;
return;

function [p, t, df] = pttest(d1, d2)
%PTTEST Student's paired t-test.
%       PTTEST(X1, X2) gives the probability that Student's t
%       calculated on paired data X1 and X2 is higher than
%       observed, i.e. the "significance" level. This is used
%       to test whether two paired samples have significantly
%       different means.
%       [P, T] = PTTEST(X1, X2) gives this probability P and the
%       value of Student's t in T. The smaller P is, the more
%       significant the difference between the means.
%       E.g. if P = 0.05 or 0.01, it is very likely that the
%       two sets are sampled from distributions with different
%       means.
%
%       This works for PAIRED SAMPLES, i.e. when elements of X1
%       and X2 correspond one-on-one somehow.
%       E.g. residuals of two models on the same data.
%
%See also: TTEST, UTTEST.
[l1 c1] = size(d1) ;
n1 = l1 * c1 ;
x1 = reshape(d1, l1 * c1, 1) ;
[l2 c2] = size(d2) ;
n2 = l2 * c2 ;
if (n1 ~= n2)
   error('PTTEST: paired samples must have the same number of elements !')
end
x2 = reshape(d2, l2 * c2, 1) ;
[a1 v1] = avevar(x1) ;
[a2 v2] = avevar(x2) ;
df  = n1 - 1 ;
cab = (x1 - a1)' * (x2 - a2) / (n1 - 1) ;
if (a1 ~= a2)
  % use abs to avoid numerical errors for very similar data
  % for which v1+v2-2cab may be close to 0.
  t = (a1 - a2) / sqrt(abs(v1 + v2 - 2 * cab) / n1) ;
  p = betainc( df / (df + t*t), df/2, 0.5) ;
else
  t = 0 ;
  p = 1 ;
end
return;

function [p, t, df] = ttest(d1, d2)
%TTEST Student's t-test for equal variances.
%       TTEST(X1, X2) gives the probability that Student's t
%       calculated on data X1 and X2, sampled from distributions
%       with the same variance, is higher than observed, i.e.
%       the "significance" level. This is used to test whether
%       two sample have significantly different means.
%       [P, T] = TTEST(X1, X2) gives this probability P and the
%       value of Student's t in T. The smaller P is, the more
%       significant the difference between the means.
%       E.g. if P = 0.05 or 0.01, it is very likely that the
%       two sets are sampled from distributions with different
%       means.
%
%       This works if the samples are drawn from distributions with
%       the SAME VARIANCE. Otherwise, use UTTEST.
%
%See also: UTTEST, PTTEST.
[l1 c1] = size(d1) ;
n1 = l1 * c1 ;
x1 = reshape(d1, l1 * c1, 1) ;
[l2 c2] = size(d2) ;
n2 = l2 * c2 ;
x2 = reshape(d2, l2 * c2, 1) ;
[a1 v1] = avevar(x1) ;
[a2 v2] = avevar(x2) ;
df = n1 + n2 - 2 ;
pvar = ((n1 - 1) * v1 + (n2 - 1) * v2) / df ;
t = (a1 - a2) / sqrt( pvar * (1/n1 + 1/n2)) ;
p = betainc( df / (df + t*t), df/2, 0.5) ;
return;

function [xbar, varx] = avevar(v)
%AVEVAR average and variance of sample.
%       AVEVAR(X) gives the average of the sample in X.
%       X is a vector of values.
%       [A, V] = AVEVAR(X) returns the average of X in A,
%       and the variance in V. The variance is corrected
%       using the two-pass formula.
%
%       Ref: [1] Chan, Golub and LeVeque. 1983. American
%                Statistician, vol. 37, pp. 242--247.
%            [2] Press et al. 1992. Numerical recipes in C.
%                Cambridge university press.
[n l] = size(v) ;
x = reshape(v, n*l, 1) ;
xbar = sum(x) / n / l ;
d = x - xbar ;
varx = ( sum(d .* d) - sum(d) * sum(d) / n / l) / (n * l - 1) ;
return;

function [p, f, df] = ftest(d1, d2)
%FTEST F-test for two samples.
%       FTEST(X1, X2) gives the probability that the F value
%       calculated as the rati of the variances of the two samples is
%       greater than observed, i.e. the significance level.
%       [P, F] = FTEST(X1, X2) gives the probability P and returns
%       the value of F.
%
%       A small value of P would lead to reject the hypothesis that
%       both data sets are sampled from distributions with the same
%       variances.
%
%See also : TTEST, TEST.
[l1 c1] = size(d1) ;
n1 = l1 * c1 ;
x1 = reshape(d1, l1 * c1, 1) ;
[l2 c2] = size(d2) ;
n2 = l2 * c2 ;
x2 = reshape(d2, l2 * c2, 1) ;
[a1 v1] = avevar(x1) ;
[a2 v2] = avevar(x2) ;
f = v1 / v2 ;
df1 = n1 - 1 ;
df2 = n2 - 1 ;
if (v1 > v2)
   p = 2 * betainc( df2 / (df2 + df1 * f), df2 / 2, df1 / 2) ;
else
   f = 1 / f ;
   p = 2 * betainc( df1 / (df1 + df2 * f), df1 / 2, df2 / 2) ;
end
if (p > 1)
   p = 2 - p ;
end

df = [df1,df2];
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [m,n] = size(a,dim)
%SIZE         Implements  size(a)  for Taylor
%
%   [m,n] = size(a,dim)
%
% functionality as Matlab function size, second argument dim optional
%

% written  05/21/09     S.M. Rump
%

  if nargout==2
    m = a.size(1);
    n = a.size(2);
  else
    if nargin==1
      m = a.size;
    else
      m = a.size(dim);
    end
  end


text:
function y = objfun(x)

y = (1/3)*x(1)^2+(1/2)*x(2)^2;

repo: Optimization-Theory-and-Algorithm


text:
function [q, g, solutionQP] = thermoFlux2QNty(model,solution,param)
% Given a steady state thermodynamically feasible flux vector, v, such that
%
%     S*v = b
%  l <= v <= u
%
% Compute q and g = N'*y such that
%
%     N*diag(q)*g = b - B*w
%
% where the stoichiometric matrix and flux vector are split into internal
% and external components: S = [N B] and v = [z; w],

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

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

if ~isfield(model,'SConsistentRxnBool')  || ~isfield(model,'SConsistentMetBool')
    if ~isfield(param,'SConsistentMethod')
        param.SConsistentMethod = 'findSExRxnInd';
    end
    switch param.SConsistentMethod
        case 'findSExRxnInd'
            %finds the reactions in the model which export/import from the model
            %boundary i.e. mass unbalanced reactions
            %e.g. Exchange reactions
            %     Demand reactions
            %     Sink reactions
            model = findSExRxnInd(model,[],param.printLevel-1);
            model.SConsistentMetBool= model.SIntMetBool;
            model.SConsistentRxnBool= model.SIntRxnBool;
        case 'findStoichConsistentSubset'
            % Finds the subset of `S` that is stoichiometrically consistent using
            % an iterative cardinality optimisation approach
            [SConsistentMetBool, SConsistentRxnBool, ...
                SInConsistentMetBool, SInConsistentRxnBool, ...
                unknownSConsistencyMetBool, unknownSConsistencyRxnBool, model] ...
                = findStoichConsistentSubset(model);
    end
else
    if length(model.SConsistentRxnBool)~=nRxn
        error('Length of model.SConsistentRxnBool must equal the number of cols of model.S')
    end
    if length(model.SConsistentMetBool)~=nMet
        error('Length of model.SConsistentMetBool must equal the number of rows of model.S')
    end
end

if ~isfield(param,'bigNum')
    param.bigNum = 1;
end

[model,rankK,nnzK,timeTaken] = internalNullspace(model);

z = solution.v(model.SConsistentRxnBool);

%                       * .A - LHS matrix
%                       * .b - RHS vector
%                       * .F - positive semidefinite matrix for quadratic part of objective (see above)
%                       * .c - Objective coeff vector
%                       * .lb - Lower bound vector
%                       * .ub - Upper bound vector
%                       * .osense - Objective sense for the linear part (-1 max, +1 min)
%                       * .csense - Constraint senses, a string containing the constraint sense for
%                         each row in A ('E', equality, 'G' greater than, 'L' less than).
[n,m]=size(model.K);
QP.A = [model.K'*diag(z); ones(1,n)];
QP.b = [zeros(m,1);1];
QP.csense(1:m+1) ='E';
QP.lb = zeros(n,1);
QP.ub = param.bigNum*ones(n,1);
QP.c = zeros(n,1);
QP.F = speye(n);
QP.osense = 1;

solutionQP = solveCobraQP(QP);

if solutionQP.stat == 1
    q0 = 1./solutionQP.full;
    q = NaN*ones(nRxn,1);
    q(model.SConsistentRxnBool)=q0;

    g0 = diag(solutionQP.full)*z;
    %N = model.S(:,model.SConsistentRxnBool);
    %y = N'\g0;
    g = NaN*ones(nRxn,1);
    g(model.SConsistentRxnBool)=g0;

    %dg0 = N'*y - g0;
    %dg = NaN*ones(nRxn,1);
    %dg(model.SConsistentRxnBool)=dg0;
else
    solutionQP
    error('thermoFlux2QNty: solveCobraQP did not solve')
end
end

author: opencobra

repo: cobratoolbox


text:
function [uCA,ndx,pos] = uniqueRowsCA(iCA,TREAT_NAN_EQUAL,FIRST_LAST)
%uniqueRowsCA  unique Rows for a Cell Array
%   U_CA = uniqueRowsCA(I_CA) returns Unique rows of the Input Cell Array, I_CA
%
%   [U_CA,I,J] = uniqueRowsCA(...) also returns index vectors I and J such
%   that U_CA = iCA(I,:) and iCA = B(J,:).
%
%   uniqueRowsCA(I_CA,TREAT_NAN_EQUAL), treats NaN comparisions as equal if
%   TREAT_NAN_EQUAL is true, the default is to treat NaN as equal
%
%   uniqueRowsCA(I_CA,[],'first'), returns the first instance
%   of the match of the index where I_CA matches U_CA, the default
%   is 'last'
%
%   NOTE: This function only handles strings and scalars, and N-D matrices
%   that are the same size for every row within the same column
%
%   See also: isequal, isequalwithequalnans

    %INPUT HANDLING
    %===========================================================
    if nargin == 1, TREAT_NAN_EQUAL = true; FIRST_LAST = 'last';
    elseif nargin == 2,
        if isempty(TREAT_NAN_EQUAL)
            TREAT_NAN_EQUAL = true;
        end
    elseif nargin ~= 3
        error('Incorrect number of input arguments to %s',mfilename)
    end

    %SOME BASIC INITIALIZATION
    %=====================================================
    [M,N] = size(iCA);
    COL_ORDER = 1:N; %Could make this an input argument,
    %affects order of unique output => see sortrows

    %HANDLING THE SIMPLE CASE
    %======================================================
    if M == 1, uCA = iCA; ndx = 1; pos = 1; return, end

    %ERROR CHECKING - CHECK CONSISTENCY FOR EACH COLUMN
    %===========================================================
    for iCol = 1:N
        %NOT YET DONE
    end

    %THE SORT_CELL_BACK_TO_FRONT FUNCTION & SORTROWS
    %=================================================================
    %Note, cell2mat is very slow, rewrote sortrows

    ndx = (1:M)';
    for k = N:-1:1
        colUse = abs(COL_ORDER(k));
        %NOTE: right here we make an assumption about the types being consistent in a given column
        if isnumeric(iCA{1, colUse})
            sz = cellfun('prodofsize',iCA(:,colUse)); %sz = cellfun(@numel,iCA(:,k));
            if isempty(find(sz > 1,1))
                tmp = fastC2M(iCA(ndx,colUse));
                ind = sortrowsc(tmp, sign(COL_ORDER(k))*1);
            else %We are dealing with matrices

                %SOME ERROR CHECKING ON THE MATRICES
                %=========================================
                szS = cellfun(@size,iCA(:,colUse),'UniformOutput',false);

                %1) DIMENSION SAME CHECK
                %----------------------------------------------
                nDims = cellfun('ndims',szS);
                if length(unique(nDims)) ~= 1
                    error('%s not setup to handle varying size matrices yet',mfilename)
                end

                %2) EACH DIMENSION THE SAME CHECK
                %----------------------------------------------
                szMat = fastC2M(szS); %NOTE: we had to first check that each input is
                %the same size, i.e. that the dimensions were the same otherwise fastC2M barfs
                %We are checking dimensions here instead of number of elements as perhaps someone
                %cares if the dimensions are different, even though the linearized contents might
                %not be ex. a = [1 3;2 4] vs b = [1 2 3 4], size(a) = [2 2], size(b) = [1 4], but
                %all(a(:) == b(:)) is true, note that the matrix code linearizes the inputs for
                %sorting purposes
                curCol = 1;
                notMatchSz = [];
                while curCol <= nDims(1) && isempty(notMatchSz)
                    notMatchSz = find(diff(szMat(:,curCol)) ~= 0,1);
                    curCol = curCol + 1;
                end

                if ~isempty(notMatchSz)
                    error('%s not setup to handle varying size matrices yet',mfilename)
                end
                %=============================================

                tmp = fastC2M(iCA(ndx,colUse));
                ind = sortrowsc(tmp, sign(COL_ORDER(k))*(1:size(tmp,2)));
            end

        else
            tmp = char(iCA(ndx,k)); %Nice and quick
            ind = sortrowsc(tmp, sign(COL_ORDER(k))*(1:size(tmp,2)));
        end
        ndx = ndx(ind);
    end

    %NOW RESORTING AND OBTAINING UNIQUE
    %=========================================
    %NOTE: ndx is the 2nd output of sortrows
    iCA = iCA(ndx,:);

    sameVector = false(M,N);

    if TREAT_NAN_EQUAL
        for iCol = 1:N
            sameVector(2:end,iCol) = arrayfun(@isequalwithequalnans,iCA(2:end,iCol),iCA(1:end-1,iCol));
        end
    else
        for iCol = 1:N
            sameVector(2:end,iCol) = arrayfun(@isequal,iCA(2:end,iCol),iCA(1:end-1,iCol));
        end
    end

    d = ~all(sameVector,2);
    uCA = iCA(d,:);

    %ADDITIONAL OUTPUT HANDLING
    %===========================================
    if nargout == 3
        pos = cumsum(d);
        pos(ndx) = pos;             % Re-reference POS to indexing of SORT.
    end

    % Create indices if needed.
    if nargout > 1
        if FIRST_LAST(1) == 'l'
            ndx = ndx([d(2:end); true]);  % Final element is always a member of unique list.
        else
            ndx = ndx(d); % First element is always a member of unique list.
        end
    else
        ndx = [];
        pos = [];
    end
end

function m = fastC2M(CA)
%ASSUMPTIONS -> all elements of CA have the same size
%TO REPLACE SLOW CELL2MAT FOR MATRICES

    szCA = size(CA);
    sz = numel(CA{1});

    m = zeros(szCA(1),sz(1));
    for iRow = 1:szCA(1)
        m(iRow,:) = CA{iRow}(:);
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function idx = sortConditions(funIn, domain, maxDiffOrders)
%SORTCONDITIONS   Returns how the results of evaluating BCs should be sorted.
%   Calling sequence:
%      IDX = SORTCONDITIONS(FUNIN, DOMAIN, DIFFORDERS)
%   where the inputs are:
%      FUNIN:         An anonymous function, that usually specifies N.LBC or
%                     N.RBC of a CHEBOP.
%      DOMAIN:        The domain of the problem.
%      MAXDIFFORDERS: The highest order of derivatives which each derivative
%                     appears in a problem.
%   and the output is
%      IDX:    A vector with indices on how to sort the results of evaluating
%              N.LBC/RBC, so that they match the order required by the MATLAB
%              ODE solvers.
%
%   Example:
%      Assume that for a coupled system, we want to specify the ICs
%         u(0) = 1, u'(0) = 2, v(0) = 0, v'(0) = 3.
%      One anonymous function we could create to specify those conditions and
%      append to a CHEBOP is
%         fun = @(u,v) [u-1; diff(v)-3; v; diff(u)-2];
%      However, when we evaluate this function to pick out the initial
%      conditions, we get the vector
%         [1; 3; 0; 2]
%      which is in the incorrect order to be passed to the MATLAB solvers, which
%      require the conditions to be of the order u, u', v, v'. Calling
%      sortConditions(),
%         idx = treeVar.sortConditions(fun, dom)
%         idx =
%             1     4     3     2
%      gives the correct order in which the above vector has to be sorted so
%      that the values are in the correct order for MATLAB.

% Check how many unknowns appear in FUNIN. When BCs are assigned via e.g.
% N.lbc = [1;3], nargin(funIn) == -1 due to the function that gets assigned
% under the hood having varargin as its argument. Hence we look at the length of
% diffOrders as well.
numArgs = max(nargin(funIn), length(maxDiffOrders));
args = cell(numArgs, 1);

% The ID vector to be passed to the TREEVAR constructor.
argsVec = zeros(1, numArgs);

% Populate the args cell with TREEVAR objects.
for argCount = 1:numArgs
    % Set the ID of the current variable to 1:
    argsVec(argCount) = 1;
    % Construct the TREEVAR:
    args{argCount} = treeVar(argsVec, domain);
    % Reset the index vector:
    argsVec = 0*argsVec;
end

% Evaluate FUNIN with the TREEVAR arguments. If we're working with the
% CHEBMATRIX syntax, we should not expand the ARGS cell (but do it otherwise):
if ( ( nargin(funIn) == 1 ) && ( length(args) > 1 ) )
    % CHEBMATRIX syntax in specifying BCs:
    bcResults = funIn(args);
else
    bcResults = funIn(args{:});
end
% Look at the results of evaluating the boundary conditions, find what
% constraint operated on what variable, and what its diffOrder was:
varList = cell(numArgs, 1);
diffOrderList = varList;

% Loop through the resulting syntax trees.
for tCounter = 1:length(bcResults)
    % Current tree we're looking at:
    tempTree = bcResults(tCounter).tree;

    % Check whether more than one variable appear in the condition, which we
    % don't support:
    assert( sum(tempTree.ID) <= 1, ...
        'CHEBFUN:TREEVAR:sortConditions:nonSeparated', ...
        'For initial value problems, only separated conditions are supported.')

    % Ensure that only initial/final conditions on the form we support appear:
    assert(acceptedCondition(tempTree), ...
        'CHEBFUN:TREEVAR:sortConditions:unsupportedCondition', [ ...
        'Initial/final condition not supported. Please ensure that the\n' ...
        'unknown function(s) appear only once in each condition,\n' ...
        'e.g. not as u + u'' - 1, and that the coefficient of the unknown\n'...
        'function is 1, e.g. not 5*u-1. Alternatively, try solving the\n' ...
        'problem in collocation mode (see help cheboppref for details).']);

    % What's the active variable in the current tree (i.e. what variable did the
    % constraint apply to)?
    activeVar = find(tempTree.ID == 1);

    % What's the diffOrder of the current constraint?
    activeDiffOrder = tempTree.diffOrder;
    activeDiffOrder = activeDiffOrder(activeVar);

    % Store in a list what variable the current constraint applies to, and what
    % the current diffOrder is:
    varList{activeVar} = [ varList{activeVar}, tCounter ];
    diffOrderList{activeVar} = [diffOrderList{activeVar}, ...
        activeDiffOrder];
end

% Initalise an index vector to be returned
idx = [];

% Go through the list of what variables appeared in what constraints, and sort
% them based on diffOrders:
for varCounter = 1:numArgs
    [varDiffOrders, diffOrderIndex] = sort(diffOrderList{varCounter});

    % Check that the user is not specifying initial values for too high order
    % derivatives.
    assert( varDiffOrders(end) < maxDiffOrders(varCounter), ...
        'CHEBFUN:TREEVAR:sortConditions:tooHighOrderCondition', ...
        'The value of a derivative of a too high order was specified.')

    % Check that multiple values for the same derivative value was not
    % specified:
    assert( length(varDiffOrders) == 1 || ~(any(diff(varDiffOrders)) == 0), ...
        'CHEBFUN:TREEVAR:sortConditions:multipleConditionsSameVariable', [...
        'Multiple initial conditions on the same variable/derivative '...
        'specified.'])

    % Check that values for all derivatives we expect to appear actually do so.
    % For IVPs/FVPs, the user has to specify the values of the function and all
    % appropriate derivatives.
    assert( all(diff(varDiffOrders) == 1) && ...
        length(varDiffOrders) == maxDiffOrders(varCounter), ...
        'CHEBFUN:TREEVAR:sortConditions:missingConditions', [...
        'Solving an nth order IVP/FVP requires specifying values of the \n' ...
        'solution and all the (n-1)st derivatives at the endpoint.'])
    assert( all( diff(varDiffOrders) == 1) && varDiffOrders(1) == 0, ...
        'myError');
    tempIndex = varList{varCounter}(diffOrderIndex);
    idx = [ idx, tempIndex ];
end

end

function accepted = acceptedCondition(tree)
%ACCEPTEDCONDITION   Check whether we have encountered an unsupported IC
%
% ACCEPTEDCONDITION(TREE) returns FALSE if we are working with functions such as
%   @(u) 5*u - 1
%   @(u) u + diff(u),
% which are not supported by the automatic first order reformulation. Otherwise,
% it returns TRUE.

% If TREE is not a struct, or it has height 0, it can't cause any harm:
if ( ~isstruct(tree) || ( tree.height == 0 ) )
    accepted = true;

% Deal with the case where the operator is + or -
elseif ( any(strcmp(tree.method, {'plus', 'minus'})) )

    % If TREE is of height 1, and it's operator is + or -, we're OK
    if ( tree.height == 1 )
        accepted = true;
        return
    end

    % Otherwise, we have to be more careful. If either argument is not a struct,
    % we only need to check the other. If both arguments are trees, we need
    % ensure that both don't have any IDs, and in the case that happens, check
    % the validity of each subtree.
    if ( ~isstruct(tree.left) )
        accepted = acceptedCondition(tree.right);
    elseif ( ~isstruct(tree.right) )
        accepted = acceptedCondition(tree.left);
    elseif ( any(tree.left.ID) && any(tree.right.ID) )
        accepted = false;
    else
        accepted = acceptedCondition(tree.right) && ...
            acceptedCondition(tree.right);
    end

% We're happy to support a tree with the diff method if it's height is 1:
elseif ( strcmp(tree.method, 'diff') && ( tree.height == 1 ) )
    accepted = true;

% Encountering any other operator implies that we had a condition on the form
% @(u) 5*u-1, which is not supported.
else
    accepted = false;
end

end


text:
function [evtNum, varLocType] = getEventNumberForVar(var, lvdData)
    evtNum = [];
    varLocType = '';

    numEvents = lvdData.script.getTotalNumOfEvents();
    for(i=1:numEvents)
        event = lvdData.script.getEventForInd(i);

        [~, eVars] = event.hasActiveOptVars();

        for(j=1:length(eVars)) %#ok<*NO4LP>
            eVar = eVars(j);

            if(strcmpi(class(var), class(eVar)) && var == eVar)
                evtNum = i;
                varLocType = 'Event';
                break;
            end
        end

        if(not(isempty(evtNum)))
            break;
        end
    end

    if(isempty(evtNum))
        if(lvdData.initStateModel.isVarFromInitialState(var))
            varLocType = 'Initial State';

        elseif(lvdData.launchVehicle.isVarFromLaunchVehicle(var))
            varLocType = 'Launch Vehicle';

        elseif(lvdData.pluginVars.isVarAPluginVar(var))
            varLocType = 'Plugins';

        else
            varLocType = '';

        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function value = r8_gmic ( a, x, alx )

%*****************************************************************************80
%
%% R4_GMIC: complementary incomplete gamma, small X, A near negative integer.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 October 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real A, the parameter.
%
%    Input, real X, the argument.
%
%    Input, real ALX, the logarithm of X.
%
%    Output, real VALUE, the complementary incomplete
%    gamma function.
%
  persistent bot
  persistent eps

  euler = 0.5772156649015329;

  if ( isempty ( eps ) )
    eps = 0.5 * r4_mach ( 3 );
    bot = log ( r4_mach ( 1 ) );
  end

  if ( 0.0 < a )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_GMIC - Fatal error!\n' );
    fprintf ( 1, '  A must be near a negative integer.\n' );
    error ( 'R4_GMIC - Fatal error!' )
  end

  if ( x <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_GMIC - Fatal error!\n' );
    fprintf ( 1, '  X <= 0.\n' );
    error ( 'R4_GMIC - Fatal error!' )
  end

  ma = r4_aint ( a - 0.5 );
  fm = - ma;
  m = - ma;

  te = 1.0;
  t = 1.0;
  s = t;
  converged = 0;
  for k = 1 : 200
    fkp1 = k + 1;
    te = - x * te / ( fm + fkp1 );
    t = te / fkp1;
    s = s + t;
    if ( abs ( t ) < eps * s )
      converged = 1;
      break
    end
  end

  if ( ~ converged )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_GMIC - Fatal error!\n' );
    fprintf ( 1, '  No convergence after 200 iterations.\n' );
    error ( 'R4_GMIC - Fatal error!' )
  end

  value = - alx - euler + x * s / ( fm + 1.0 );

  if ( m == 0 )
    return
  elseif ( m == 1 )
    value = - value - 1.0 + 1.0 / x;
    return
  end

  te = fm;
  t = 1.0;
  s = t;
  mm1 = m - 1;
  for k = 1 : mm1
    fk = k;
    te = - x * te / fk;
    t = te / ( fm - fk );
    s = s + t;
    if ( abs ( t ) < eps * abs ( s ) )
      break
    end
  end

  for k = 1 : m
    value = value + 1.0 / k;
  end

  if ( mod ( m, 2 ) == 1 )
    sgng = - 1.0;
  else
    sgng = + 1.0;
  end

  alng = log ( value ) - r4_lngam ( fm + 1.0 );

  if ( bot < alng )
    value = sgng * exp ( alng );
  else
    value = 0.0;
  end

  if ( s ~= 0.0 )
    value = value ...
      + abs ( exp ( - fm * alx + log ( abs ( s ) / fm ) ) * r8_sign ( s );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function et = set(t,et)
% TIMETIC/SET Set the elapsed time of a timetic object.
%
% Setting a timer object always stops the timer.
%
% Example:
%   t = timetic;
%   set(t,5);
%   start(t);
%   toc(t);

%
% David Gleich
% Stanford University
% 27 September 2006
%

t.set(et);


text:
function [x, w, v] = radaupts(n, alp, bet)
%RADAUPTS   Gauss-Jacobi-Radau quadrature nodes and weights.
%  RADAUPTS(N) returns N Legendre-Radau points X in [-1,1).
%
%  [X, W] = RADAUPTS(N) returns also a row vector W of weights for
%  Gauss-Legendre-Lobatto quadrature.
%
%  [X, W, V] = RADAUPTS(N) returns additionally a column vector V of weights in
%  the barycentric formula corresponding to the points X. The weights are scaled
%  so that max(abs(V)) = 1.
%
%  [...] = RADUAPTS(N, ALP, BET) is similar, but for the Gauss-Jacobi-Radau
%  nodes and weights. Here ALP and BET should be scalars > -1.
%
%  In each case, N should be a positive integer.
%
% See also CHEBPTS, LEGPTS, JACPTS, LEGPOLY, JACPOLY, LOBPTS.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Developer note:
%   The approach used here is to observe that the Gauss-Radau points are
%   precisely the roots of (1+x)P^(0,1)_{n-2}(x), which can be obtained from
%   JACPTS. A similar identity [NIST, (18.9.5) and (18.9.17)] is used for the
%   computation of the quadrature weights from those of JACPTS, and the missing
%   barycentric weights are determined by enforcing the interpolation of f(x) =
%   x at x = 0.
%
%    x_j = roots of (1+x)P^(0,1)_{n-2}(x)
%    w_j = { 2/n^2                                 : x_j = -1
%          { 1/(1-x_j) * 1/[d/dx P_{n-1}(x_j)]^2   : otherwise
%
%   (Note that the weights for n-1 point Gauss-Jacobi with a = 0, b = 1 satisfy
%    u_j = C/(1-x_j^2)/[d/dx P^(0,1)_{n-1}(x_j)]^2)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% TODO: Scaled domains?

if ( nargin < 3 )
    % Default to Gauss-Legendre-Lobatto:
    alp = 0; bet = 0;
end

%% Trivial cases:
if ( n == 1 )
    x = -1;
    w = 2^(1+alp+bet)*beta(1+alp,1+bet);
    v = 1;
    return
end

%% Call JACPTS():
[xi, w, v] = jacpts(n-1, alp, bet+1);

%% Nodes:
x = [-1 ; xi];

%% Quadrature weights:
wi = w./(1 + xi.');
if ( alp == 0 && bet == 0 )
    w = [2/n^2, wi];
else
    % See Walter Gautschi, "Gauss–Radau formulae for Jacobi and Laguerre
    % weight functions", Mathematics and Computers in Simulation, (2000).
    w = [2^(alp+bet+1)*beta(bet+1,n)*beta(alp+n,bet+1)*(bet+1), wi];
end

%% Barycentric weights:
v = v./(1 + xi);
v = v/max(abs(v));
v1 = -sum(v);
v = [v1 ; v];

end


text:
clear all
addpath('../ParNMPC/')
%% Formulate an OCP using Class OptimalControlProblem

% Create an OptimalControlProblem object
OCP = OptimalControlProblem(1,... % dim of inputs
                            2,... % dim of states
                            0,... % dim of parameters
                            12);  % N: num of discritization grids

% Set the prediction horizon T
OCP.setT(1);

% Set the dynamic function f
a = -1;
b = -1;
f = [OCP.x(2); a * OCP.x(1) + b * OCP.x(2) * OCP.u(1)];
OCP.setf(f);
OCP.setDiscretizationMethod('Euler');

% Set the cost function L
Q    = diag([10, 10]);
R    = 0.01;
xRef = [0;0];
uRef = 0.5;
L    =  0.5*(OCP.x-xRef).'*Q*(OCP.x-xRef)...
       +0.5*(OCP.u-uRef).'*R*(OCP.u-uRef);
OCP.setL(L);

% Set the linear constraints G(u,x,p)>=0
G = [OCP.u(1);...
     1 - OCP.u(1)];
OCP.setG(G);

% Generate necessary files
OCP.codeGen();
%% Configrate the solver using Class NMPCSolver

% Create a NMPCSolver object
nmpcSolver = NMPCSolver(OCP);

% Configurate the Hessian approximation method
nmpcSolver.setHessianApproximation('Newton');

% Generate necessary files
nmpcSolver.codeGen();
%% Solve the very first OCP for a given initial state and given parameters using Class OCPSolver

% Set the initial state
x0 =   [1;0];

% Set the parameters
dim = OCP.dim;
N   = OCP.N;
p   = zeros(dim.p,N);

% Solve the very first OCP
solutionInitGuess.lambda = [randn(dim.lambda,1),zeros(dim.lambda,1)];
solutionInitGuess.mu     = randn(dim.mu,1);
solutionInitGuess.u      = uRef;
solutionInitGuess.x      = [x0,xRef];
solutionInitGuess.z      = ones(dim.z,N);
solution = NMPC_SolveOffline(x0,p,solutionInitGuess,1e-4,1000);

plot(solution.x(1,:).');

% Save to file
save GEN_initData.mat dim x0 p N

% Set initial guess
global ParNMPCGlobalVariable
ParNMPCGlobalVariable.solutionInitGuess = solution;
%% Define the controlled plant using Class DynamicSystem

% M(u,x,p) \dot(x) = f(u,x,p)
% Create a DynamicSystem object
plant = DynamicSystem(1,2,0);

g = 9.81;
a = -1;
b = -1;
fPlant = [plant.x(2); a * plant.x(1) + b * plant.x(2) * plant.u(1)];

% Set the dynamic function f
plant.setf(fPlant); % same model

% Generate necessary files
plant.codeGen();

author: deng-haoyang

repo: ParNMPC


text:
avg_odd = 3.93;
max_odd = 4.27;

p_margin = 0.04;
p_real = (1 / avg_odd) - p_margin; % pur estimate of the "true" probability

p_max = 1 / max_odd;

payoff = p_real * max_odd - 1;

if payoff > 0
    msg = sprintf('Payoff = %2.2f. Place bet', payoff);
    disp(msg)
else
    msg = sprintf('Payoff = %2.2f. Do not Place bet', payoff);
    disp(msg)
end

repo: BeatTheBookie


text:
classdef L1VectorNormProximal < handle

    properties (Access = private)
       lambda
       imageSize
       designVariable
    end

    methods (Access = public)

        function obj = L1VectorNormProximal(cParams)
            obj.lambda = cParams.lambda;
            obj.designVariable = cParams.designVariable;
            obj.imageSize = cParams.imageSize;
        end

        function solve(obj)
            lam = obj.lambda;
            normP = obj.computeNormP();
            no  = max(1,normP/lam);
            p = obj.designVariable.value;
            p = p./[no;no];
            obj.designVariable.value = p;
        end

    end

    methods (Access = private)

       function normP = computeNormP(obj)
            p  = obj.designVariable.value;
            mn = obj.imageSize.rowsTimesColumns;
            normP = hypot(p(1:mn),p(mn+1:end));
        end

    end

end

author: SwanLab

repo: Swan


text:
function matreg = Tikhonov_rank_def(mat, rank, lambda)

% Apply Tikhonov regularisation to rank-deficient matrix

% mat: square matrix

% rank: number of singular values to be considered in inversion

% lambda: regularisation parameters

% OH, Sep 2018

% SVD of input matrix

[U,S,V] = svd(mat);

% get singular values

s = diag(S);

% take only relevant values

s2 = s(1:rank);

lambda = (lambda/100) * (sum(s2)/rank);

% regularise eigenvalues with Tikhonov and invert

s2 = s2 ./ (s2.^2 + lambda);

% reconstitute regularised inverse matrix

matreg = V(:,1:rank)*diag(s2)*U(:,1:rank)';

author: spm


text:
function stroud_test35 ( )

%*****************************************************************************80
%
%% TEST35 tests SQUARE_UNIT_SET, RECTANGLE_SUB_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  global FUNC_2D_INDEX;

  num = function_2d_num ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST35\n' );
  fprintf ( 1, '  SQUARE_UNIT_SET sets up a quadrature rule\n' );
  fprintf ( 1, '    on a unit square.\n' );
  fprintf ( 1, '  RECTANGLE_SUB_2D applies it to subrectangles of an\n' );
  fprintf ( 1, '    arbitrary rectangle.\n' );
  fprintf ( 1, '\n' );
%
%  Set the location of the square.
%
  xval(1) = 1.0;
  yval(1) = 2.0;

  xval(2) = 3.0;
  yval(2) = 3.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The corners of the rectangle are:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  %12f  %12f\n', xval(1), yval(1) );
  fprintf ( 1, '  %12f  %12f\n', xval(2), yval(2) );
%
%  Get the quadrature abscissas and weights for a unit square.
%
  rule = 2;

  order = square_unit_set ( rule );

  [ xtab, ytab, weight ] = square_unit_set ( rule, order );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Using unit square integration rule number %d\n', rule );
%
%  Set the function.
%
  for i = 1 : num

    FUNC_2D_INDEX = i;
%
%  Try an increasing number of subdivisions.
%
    fprintf ( 1, '\n' );
    fprintf ( 1, '    Function  Subdivisions  Integral\n' );
    fprintf ( 1, '\n' );

    for j = 1 : 5

      nsub(1) = j;
      nsub(2) = 2 * j;

      result = rectangle_sub_2d ( 'function_2d', xval, yval, nsub, order, xtab, ...
        ytab, weight );

      fname = function_2d_name ( i );

      fprintf ( 1, '  %s  %2d  %2d  %14f\n', fname, nsub(1), nsub(2), result );

    end

  end

  return
end


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

%*****************************************************************************80
%
%% COS_DEGREE_VALUES: the cosine function with argument in degrees.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      Cos[x Degree]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Cambridge University Press, 1999,
%    ISBN: 0-521-64314-7,
%    LC: QA76.95.W65.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0
%    before the first call.  On each call, the routine increments N_DATA by 1,
%    and returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 22;

  fx_vec = [ ...
     0.99619469809174553230, ...
     1.0000000000000000000, ...
     0.99984769515639123916, ...
     0.99939082701909573001, ...
     0.99862953475457387378, ...
     0.99756405025982424761, ...
     0.99619469809174553230, ...
     0.98480775301220805937, ...
     0.96592582628906828675, ...
     0.86602540378443864676, ...
     0.70710678118654752440, ...
     0.50000000000000000000, ...
     0.25881904510252076235, ...
     0.087155742747658173558, ...
     0.069756473744125300776, ...
     0.052335956242943832722, ...
     0.034899496702500971646, ...
     0.017452406437283512819, ...
     0.000000000000000000000, ...
    -0.017452406437283512819, ...
    -0.25881904510252076235, ...
    -1.0000000000000000000 ];
  x_vec = [ ...
     -5.0, ...
      0.0, ...
      1.0, ...
      2.0, ...
      3.0, ...
      4.0, ...
      5.0, ...
     10.0, ...
     15.0, ...
     30.0, ...
     45.0, ...
     60.0, ...
     75.0, ...
     85.0, ...
     86.0, ...
     87.0, ...
     88.0, ...
     89.0, ...
     90.0, ...
     91.0, ...
    105.0, ...
    180.0 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

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

  return
end


text:
function obj = display(obj)
	% Biafra Ahanonu
	% Started: 2021.03.25 [22:11:25] (branched from ciatah.m)

	% Overload display method so can run object by just typing 'obj' in command window.
	obj.runPipeline;
	% display('hello');
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [X, y] = face_gen(num, noise)

n1=floor(0.05*num);
n3=floor(0.15*num);
n4=floor(0.75*num);

interval1=1.25; interval2=1.25;
% left eye
m = [-interval1,interval2];
C = 0.5*noise*eye(2);
x1 = mvnrnd(m,C,n1);
y1 = 1+zeros(n1,1);
% right eye
m = [interval1,interval2];
x2 = mvnrnd(m,C,n1);
y2 = 2+zeros(n1,1);

% nose
x3=zeros(n3,2);
r=1.5;
t=(5/4*pi):pi/(2*n3-1):(7/4*pi); % ÏÂ°ëÔ²
x3(:, 1) = r.*cos(t)'+randn(n3,1)*noise;
x3(:, 2) = r.*sin(t)'+randn(n3,1)*noise;
y3 = 3+zeros(n3,1);

% face
upright = 0.5;
x4=zeros(n4,2);
r = 3;
curve = 2.5;
t = unifrnd(0,0.8,[1,n4]);
x4(:, 1) = r.*sin(curve*pi*t) + noise*randn(1,n4);
x4(:, 2) = r.*cos(curve*pi*t) + noise*randn(1,n4)+upright;;
y4 = 4+zeros(n4,1);

X = [x1;x2;x3;x4];
y = [y1;y2;y3;y4];

author: CHLWR