text:
switch input_sig
    case 1          % 1: synthetic data
        if load_sig
            load matlab
        else
        A_0=randn(n1,data.r)*randn(data.r,n2);  % true low-rank matrix A_0
%         A_max=max(abs(A_0(:)));
        A_max=sqrt(max(n1,n2));
        B_0=((rand(n1,n2)<.5)-.5)*2*A_max;
        B_0(randperm(n1*n2,n1*n2-data.s))=0;    % true sparse matrix B_0
%         X_0=A_0+B_0; Z=X_0;
        Z=A_0+B_0;
%         Z=randn(size(Z));
        end

    case 2          % 2: airport video
        addpath('./video_airport')
        nFrame=200;
        iFrame=1000;
        pathname=['airport',num2str(iFrame),'.bmp'];
        Z1=imread(pathname);
        Z1=double(uint8(round(sum(Z1,3)/3)))/255;
        Z=zeros(size(Z1,1),size(Z1,2),nFrame);
        Z(:,:,1)=Z1;
        for j=1:nFrame-1
            pathname=['airport',num2str(iFrame+j),'.bmp'];
            Z1=imread(pathname);
            Z1=double(uint8(round(sum(Z1,3)/3)))/255;
            Z(:,:,j+1)=Z1;
        end

        [z1,z2,z3]=size(Z);
        n1=z1*z2; n2=z3;
        Z=reshape(Z,n1,n2);

        clear nFrame iFrame pathname Z1

    case 3          % 3: lobby video
        addpath('./video_lobby')
        nFrame=400;
        iFrame=1900;
        pathname=['SwitchLight',num2str(iFrame),'.bmp'];
        Z1=imread(pathname);
        Z1=double(uint8(round(sum(Z1,3)/3)))/255;
        Z=zeros(size(Z1,1),size(Z1,2),nFrame);
        Z(:,:,1)=Z1;
        for j=1:nFrame-1
            pathname=['SwitchLight',num2str(iFrame+j),'.bmp'];
            Z1=imread(pathname);
            Z1=double(uint8(round(sum(Z1,3)/3)))/255;
            Z(:,:,j+1)=Z1;
        end

        [z1,z2,z3]=size(Z);
        n1=z1*z2; n2=z3;
        Z=reshape(Z,n1,n2);

        clear nFrame iFrame pathname Z1
end

% add noise
% noise=.0;
% noise=.05*std(data);
if ~load_sig && noise~=0
    Z=Z+noise*randn(size(Z));                       % Gaussian noise
end

author: andrewssobral

repo: lrslibrary


text:
function [z] = griddata_version_control(gridX,gridY,gridZ,x,y,method,matlab_version)
% This function just computes the gridded data set, but its input is
% modified to work with matlab version newer than 2012.
% By David Bekaert - December 2012

if matlab_version>=2012
    z=griddata(gridX,gridY,gridZ,x,y,method);
else
    z=griddata(gridX,gridY,gridZ,x,y,method,{'QJ'});
end

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function [im,lv,rv,myo] = ventricles_phantom(m,v)
% VENTRICLES_PHANTOM Create phantom image of simplified heart ventricles
%
% [IM,LV,RV,MYO] = ventricles(M,V)
%
%   IM is a 3D image that consists of a heart phantom depicting the left
%   ventricle, right ventricle and myocardium
%
%   LV is a 3D image that consists of the LV cavity phantom
%
%   RV is a 3D image that consists of the RV cavity phantom
%
%   MYO is a 3D image that consists of a myocardium phantom
%
%   M is the mean Gaussian white noise added to the final heart phantom
%   generated (the default value is zero)
%
%   V is the the variance of the Gaussian white noise added to the final
%   heart phantom generated (the default value is zero)

% Author: Christopher Kelly <christopher.kelly28@gmail.com>
% Copyright © 2013 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments

narginchk(0,2);
nargoutchk(0,4);

% default values

if (nargin == 1)
    error('You must either specify all 2 input parameters or none at all')
elseif (nargin < 2)
    m      = 0;
    v      = 0;
end

% TODO: we hard-code these values. In the future, we may want to have these
% as input arguments, but currently changing these values produces odd
% looking results
lvendo = 300;
lvepi  = 400;
rvendo = 300;
rvepi  = 400;

%% create image space

xRange = -(128/2) + 0.5 : (128/2) - 0.5;
yRange = xRange;
zRange = xRange;

[ x, y, z ] = ndgrid( xRange, yRange, zRange );

im = zeros(size(x));

lv = im;
rv = im;
myo = im;
heart = im;

%% define ellipsoids shapes

radius1 = ( (x).^2 + (y).^2 + (z/2).^2 );
radius2 = ( (x).^2 + (y/1.7).^2 + (z/2).^2 );

%% define voxels in left ventricle wall

% define LV ellpsoid edge thickness
ellipsoid1Outer = find(radius1 < lvepi );
ellipsoid1Inner = find(radius1 < lvendo );

% define LV ellipsoid
ellipsoid1 = ismember(ellipsoid1Outer,ellipsoid1Inner);
ellipsoid1 = ellipsoid1Outer(ellipsoid1 == 0);

[ellipsoid1X, ellipsoid1Y,ellipsoid1Z] = ind2sub(size(im),ellipsoid1);

ellipsoid1X(ellipsoid1Z > 70) = [];
ellipsoid1Y(ellipsoid1Z > 70) = [];
ellipsoid1Z(ellipsoid1Z > 70) = [];

% find voxel numbers that will make up LV wall
LV_wall = sub2ind(size(im),ellipsoid1X,ellipsoid1Y,ellipsoid1Z);

%% define voxels in right ventricle wall

% define RV ellpsoid edge thickness
ellipsoid2Outer = find(radius2 < rvepi );
ellipsoid2Inner = find(radius2 < rvendo );

% define RV ellipsoid
ellipsoid2 = ismember(ellipsoid2Outer, ellipsoid2Inner);
ellipsoid2 = ellipsoid2Outer(ellipsoid2 == 0);                                	% find voxel numbers that will make up RV wall

[ellipsoid2X, ellipsoid2Y,ellipsoid2Z] = ind2sub(size(im),ellipsoid2);

ellipsoid2X(ellipsoid2Z > 70) = [];
ellipsoid2Y(ellipsoid2Z > 70) = [];
ellipsoid2Z(ellipsoid2Z > 70) = [];

ellipsoid2X(ellipsoid2Y > (size(im,2)/2)) = [];
ellipsoid2Z(ellipsoid2Y > (size(im,2)/2)) = [];
ellipsoid2Y(ellipsoid2Y > (size(im,2)/2)) = [];

RV_wall = sub2ind(size(im),ellipsoid2X,ellipsoid2Y,ellipsoid2Z);

%% define left ventricle

lv(radius1 < lvendo) = 1;

LVellipsoid = find(lv == 1);
[LVellipsoidX, LVellipsoidY, LVellipsoidZ] = ind2sub(size(im),LVellipsoid);

LVellipsoidX(LVellipsoidZ < 70) = [];
LVellipsoidY(LVellipsoidZ < 70) = [];
LVellipsoidZ(LVellipsoidZ < 70) = [];

LVellipsoid = sub2ind(size(im),LVellipsoidX, LVellipsoidY, LVellipsoidZ);
lv(LVellipsoid) = 0;

%% define right ventricle

rv(radius2 < rvendo) = 1;

RVellipsoid = find(rv == 1);
[RVellipsoidX, RVellipsoidY, RVellipsoidZ] = ind2sub(size(im),RVellipsoid);

RVellipsoidX( (RVellipsoidZ < 70) ) = [];
RVellipsoidY( (RVellipsoidZ < 70) ) = [];
RVellipsoidZ( (RVellipsoidZ < 70) ) = [];

RVellipsoid = sub2ind(size(im),RVellipsoidX, RVellipsoidY, RVellipsoidZ);
rv(RVellipsoid) = 0;

RVellipsoid = find(rv == 1);
[RVellipsoidX, RVellipsoidY, RVellipsoidZ] = ind2sub(size(im),RVellipsoid);

RVellipsoidX(RVellipsoidY < (size(im,2)/2)) = [];
RVellipsoidZ(RVellipsoidY < (size(im,2)/2)) = [];
RVellipsoidY(RVellipsoidY < (size(im,2)/2)) = [];

RVellipsoid = sub2ind(size(im),RVellipsoidX, RVellipsoidY, RVellipsoidZ);
rv(RVellipsoid) = 0;

rv([LV_wall;find(lv == 1)]) = 0;

%% define myocardium

myo([LV_wall;RV_wall]) = 1;

%% define heart

heart( lv == 1 ) = 1;
heart( rv == 1 ) = 2;
heart( myo == 1 ) = 3;

%% add noise

heart = (heart - min(heart(:)))/(max(max(max(heart)))-min(min(min(heart))));

for i = 1:size(im,3)
    heart(:,:,i) = imnoise(heart(:,:,i),'Gaussian',m,v);
    lv(:,:,i) = imnoise(lv(:,:,i),'Gaussian',m,v);
    rv(:,:,i) = imnoise(rv(:,:,i),'Gaussian',m,v);
    myo(:,:,i) = imnoise(myo(:,:,i),'Gaussian',m,v);
end

im = heart;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function bellman_ford_test ( )

%*****************************************************************************80
%
%% BELLMAN_FORD_TEST tests the BELLMAN_FORD library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BELLMAN_FORD_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the BELLMAN_FORD library.\n' );

  bellman_ford_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BELLMAN_FORD_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function  [Z, sigma] =  MCWNNM_ADMM_NL1( Y, NSig, Par )
% This routine solves the following weighted nuclear norm optimization problem with column weights,
%
% min_{X, Z} ||W(Y-X)||_F^2 + ||Z||_w,*  s.t.  X = Z
%
% Inputs:
%        Y        -- 3p^2 x M dimensional noisy matrix, D is the data dimension, and N is the number of image patches.
%        NSig   -- 3p^2 x 1 dimensional vector of weights
%        Par     -- structure of parameters
% Output:
%        Z        -- 3p^2 x M dimensional denoised matrix
%        sigma -- the noise standard deviation

% tol = 1e-8;
if ~isfield(Par, 'maxIter')
    Par.maxIter = 10;
end
if ~isfield(Par, 'rho')
    Par.rho = 1;
end
if ~isfield(Par, 'mu')
    Par.mu = 1;
end
if ~isfield(Par, 'display')
    Par.display = true;
end
%% Initializing optimization variables
% intialize

% Initializing optimization variables
% Intialize the weight matrix W
mNSig = min(NSig);
W = (mNSig+eps) ./ (NSig+eps);

Z = zeros(size(Y));
A = zeros(size(Y));
%% Start main loop
iter = 0;
PatNum       = size(Y,2);
TempC  = Par.Constant * sqrt(PatNum) * mNSig^2;
% TempC  = Par.Constant * sqrt(PatNum);
while iter < Par.maxIter
    iter = iter + 1;

    % update X, fix Z and A
    % min_{X} ||W * Y - W * X||_F^2 + 0.5 * rho * ||X - Z + 1/rho * A||_F^2
    X = diag(1 ./ (W.^2 + 0.5 * Par.rho)) * (diag(W.^2) * Y + 0.5 * Par.rho * Z - 0.5 * A);

    % update Z, fix X and A
    % min_{Z} ||Z||_*,w + 0.5 * rho * ||Z - (X + 1/rho * A)||_F^2
    Temp = X + A/Par.rho;
    [U, SigmaTemp, V] =   svd(full(Temp), 'econ');
    [SigmaZ, svp] = ClosedWNNM(diag(SigmaTemp), 2/Par.rho*TempC, eps);
    Z =  U(:, 1:svp) * diag(SigmaZ) * V(:, 1:svp)';
    %     % check the convergence conditions
    %     stopC = max(max(abs(X - Z)));
    %     if Par.display && (iter==1 || mod(iter,10)==0 || stopC<tol)
    %         disp(['iter ' num2str(iter) ',mu=' num2str(Par.mu,'%2.1e') ...
    %             ',rank=' num2str(rank(Z,1e-4*norm(Z,2))) ',stopALM=' num2str(stopC,'%2.3e')]);
    %     end
    %     if stopC < tol
    %         break;
    %     else
    % update the multiplier A, fix Z and X
    A = A + Par.rho * (X - Z);
    Par.rho = min(1e4, Par.mu * Par.rho);
    %     end
end
sigma = sqrt( mean( reshape(mean((Y - Z).^2, 2)', [Par.ps2 Par.ch])) )';
return;

author: csjunxu


text:
function [fval,optimState,idx] = funlogger_vbmc(fun,x,optimState,state,varargin)
%FUNLOGGER_VBMC Call objective function and do some bookkeeping.
%   [~,OPTIMSTATE] = FUNLOGGER_VBMC(FUN,NVARS,OPTIMSTATE,'init') starts logging
%   function FUN with starting point X and optimization struct OPTIMSTATE.
%
%   [~,OPTIMSTATE] = FUNLOGGER_VBMC(FUN,NVARS,OPTIMSTATE,'init',NMAX) stores up
%   to NMAX function values (default NMAX=1e4).
%
%   [FVAL,OPTIMSTATE] = FUNLOGGER_VBMC(FUN,X,OPTIMSTATE,'iter') evaluates
%   function FUN at X with optimization struct OPTIMSTATE, returns function
%   value FVAL. FUN must take a vector input and return a scalar value and,
%   optionally, the (estimated) SD of the returned value (if heteroskedastic
%   noise handling is on).
%
%   [FVAL,OPTIMSTATE] = FUNLOGGER_VBMC(FUN,X,OPTIMSTATE,'single') as 'iter'
%   but does not store function values in the cache.
%
%   [~,OPTIMSTATE] = FUNLOGGER_VBMC(FUN,X,OPTIMSTATE,'done') finalizes
%   stored function values.

%   Luigi Acerbi 2019

fval = [];

% Switch warnings on
for ff = fields(optimState.DefaultWarnings)'
    warning(optimState.DefaultWarnings.(ff{:}).state,optimState.DefaultWarnings.(ff{:}).identifier);
end

switch lower(state)
    case 'init' % Start new function logging session

        nvars = x;

        % Number of stored function values
        if nargin > 4; nmax = varargin{1}; else; nmax = []; end
        if isempty(nmax); nmax = 500; end

        % Read noise level
        noise_flag = optimState.UncertaintyHandlingLevel > 0;

        optimState.funccount = 0;
        optimState.cachecount = 0;

        % Passing empty or new OPTIMSTATE
        if ~isfield(optimState,'X_orig') || isempty(optimState.X_orig)
            optimState.X_orig = NaN(nmax,nvars);
            optimState.y_orig = NaN(nmax,1);
            optimState.X = NaN(nmax,nvars);
            optimState.y = NaN(nmax,1);
            if noise_flag; optimState.S = NaN(nmax,1); end
            optimState.nevals = zeros(nmax,1);
            optimState.Xn = 0;                    % Last filled entry
            optimState.Xmax = nmax;               % Maximum entry index
            optimState.X_flag = false(nmax,1);
            optimState.ymax = -Inf;

        else % Receiving previous evaluations (e.g. from previous run)

            error('Previous function evaluations not supported yet.');

            % Find last evaluated point
            optimState.Xmax = find(~isnan(optimState.Y),1,'last');
            if ~isfield(optimState,'Xn') || isempty(optimState.Xn)
                optimState.Xn = optimState.Xmax;
            end

            % Current cache is smaller than previous one, keep most recent
            % evaluated points
            if optimState.Xmax > nmax
                optimState.X_orig = circshift(optimState.X_orig,-optimState.Xn);
                optimState.X_orig = optimState.X_orig(end-nmax+1:end,:);
                optimState.y_orig = circshift(optimState.y_orig,-optimState.Xn);
                optimState.y_orig = optimState.y_orig(end-nmax+1:nmax,:);
                if noise_flag
                    optimState.S = circshift(optimState.S,-optimState.Xn);
                    optimState.S = optimState.S(end-nmax+1:nmax,:);
                end
                optimState.Xn = 0;
                optimState.Xmax = nmax;
            else
                offset = nmax - size(optimState.X_orig,1);
                optimState.X_orig = [optimState.X_orig; NaN(offset,nvars)];
                optimState.y_orig = [optimState.y_orig; NaN(offset,1)];
                if noise_flag
                    optimState.S = [optimState.S; NaN(offset,1)];
                end
            end
            optimState.X = warpvars_vbmc(optimState.X_orig,'d',optimState.trinfo);
        end
        optimState.funevaltime = NaN(nmax,1);
        optimState.totalfunevaltime = 0;

    case {'iter','single'} % Evaluate function (and store output for 'iter')

        x_orig = warpvars_vbmc(x,'inv',optimState.trinfo);    % Convert back to original space
        % Heteroscedastic noise?
        noise_flag = optimState.UncertaintyHandlingLevel > 0;

        try
            funtime = tic;
            if noise_flag && optimState.UncertaintyHandlingLevel == 2
                [fval_orig,fsd] = fun(x_orig);
            else
                fval_orig = fun(x_orig);
                if noise_flag
                    fsd = 1;
                else
                    fsd = [];
                end
            end
            t = toc(funtime);

            % Check returned function value
            %if isscalar(fval_orig) && isreal(fval_orig) && fval_orig == -Inf
            %    warning(['funlogger_vbmc:InfiniteFuncValue',...
            %        'The function returned -Inf as function value, which should not be allowed. Trying to continue, but results might be affected.'])
            %    fval_orig = log(realmin);
            if ~isscalar(fval_orig) || ~isfinite(fval_orig) || ~isreal(fval_orig)
                error(['funlogger_vbmc:InvalidFuncValue',...
                    'The returned function value must be a finite real-valued scalar (returned value: ' mat2str(fval_orig) ').']);
            end

            % Check returned function SD
            if noise_flag && ...
                    (~isscalar(fsd) || ~isfinite(fsd) || ~isreal(fsd) || fsd <= 0.0)
                error(['funlogger_vbmc:InvalidNoiseValue',...
                    'The returned estimated SD (second function output) must be a finite, positive real-valued scalar (returned SD: ' mat2str(fsd) ').']);
            end

            % Tempered posterior
            if isfield(optimState,'temperature') && ~isempty(optimState.temperature)
                fval_orig = fval_orig / optimState.temperature;
                fsd = fsd / optimState.temperature;
            end

        catch fun_error
            warning(['funlogger_vbmc:FuncError',...
                'Error in executing the logged function ''' func2str(fun) ''' with input: ' mat2str(x)]);
            fprintf('Error identifier: %s\nMessage: %s\n',fun_error.identifier,fun_error.message);
            fprintf('File: %s\n',fun_error.stack(1).file);
            fprintf('Function: %s, line %d.\n',fun_error.stack(1).name,fun_error.stack(1).line);
            rethrow(fun_error);
        end

        % Update function records
        optimState.funccount = optimState.funccount + 1;
        if strcmpi(state, 'iter')
            [optimState,fval,idx] = record(optimState,x_orig,x,fval_orig,t,fsd);
        end
        optimState.totalfunevaltime = optimState.totalfunevaltime + t;

    case {'add'} % Add previously evaluated function

        fval_orig = varargin{1};
        x_orig = warpvars_vbmc(x,'inv',optimState.trinfo);    % Convert back to original space
        % Heteroscedastic noise?
        noise_flag = isfield(optimState,'S');
        if noise_flag
            fsd = varargin{2};
            if isempty(fsd); fsd = 1; end
        else
            fsd = [];
        end

        % Check function value
        if ~isscalar(fval_orig) || ~isfinite(fval_orig) || ~isreal(fval_orig)
            error(['funlogger_vbmc:InvalidFuncValue',...
                'The provided function value must be a finite real-valued scalar (provided value: ' mat2str(fval_orig) ').']);
        end

        % Check returned function SD
        if noise_flag && (~isscalar(fsd) || ~isfinite(fsd) || ~isreal(fsd) || fsd <= 0.0)
            error(['funlogger_vbmc:InvalidNoiseValue',...
                'The provided estimated SD (second function output) must be a finite, positive real-valued scalar (provided SD: ' mat2str(fsd) ').']);
        end

        % Tempered posterior
        if isfield(optimState,'temperature') && ~isempty(optimState.temperature)
            fval_orig = fval_orig / optimState.temperature;
            fsd = fsd / optimState.temperature;
        end

        % Update function records
        optimState.cachecount = optimState.cachecount + 1;
        [optimState,fval,idx] = record(optimState,x_orig,x,fval_orig,0,fsd);

    case 'done' % Finalize stored table

        noise_flag = isfield(optimState,'S');

        optimState.X_orig = optimState.X_orig(1:optimState.Xn,:);
        optimState.y_orig = optimState.y_orig(1:optimState.Xn);
        optimState.X_flag = optimState.X_flag(1:optimState.Xn);
        if noise_flag; optimState.S = optimState.S(1:optimState.Xn); end
        optimState.funevaltime = optimState.funevaltime(1:optimState.Xn);
        optimState.nevals = optimState.nevals(1:optimState.Xn);

        optimState = rmfield(optimState,'X');
        optimState = rmfield(optimState,'y');

    otherwise
        error('funlogger_vbmc:UnknownAction','Unknown FUNLOGGER action.');
end

% Switch warnings off again
for ff = fields(optimState.DefaultWarnings)'
    warning('off',optimState.DefaultWarnings.(ff{:}).identifier);
end

end

%--------------------------------------------------------------------------
function [optimState,fval,idx] = record(optimState,x_orig,x,fval_orig,t,fsd)
%RECORD Record function evaluation.

if nargin < 6; fsd = []; end

duplicate_flag = all(bsxfun(@eq,x,optimState.X),2);

% Tempered posterior
if isfield(optimState,'temperature') && ~isempty(optimState.temperature)
    T = optimState.temperature;
else
    T = 1;
end

if any(duplicate_flag)
    if sum(duplicate_flag) > 1; error('More than one match.'); end
    idx = find(duplicate_flag);
    N = optimState.nevals(idx);

    if ~isempty(fsd)
        tau_n = 1/optimState.S(idx)^2;
        tau_1 = 1/fsd^2;
        optimState.y_orig(idx) = (tau_n*optimState.y_orig(idx) + tau_1*fval_orig)/(tau_n + tau_1);
        optimState.S(idx) = 1/sqrt(tau_n + tau_1);
    else
        optimState.y_orig(idx) = (N*optimState.y_orig(idx) + fval_orig)/(N+1);
    end

    fval = optimState.y_orig(idx) + warpvars_vbmc(x,'logp',optimState.trinfo)/T;
    optimState.y(optimState.Xn) = fval;
    % if ~isempty(fsd); optimState.S(idx) = sqrt((N^2*optimState.S(idx)^2 + fsd^2)/(N+1)^2); end
    optimState.funevaltime(idx) = (N*optimState.funevaltime(idx) + t)/(N+1);
    optimState.nevals(idx) = optimState.nevals(idx) + 1;

else
    optimState.Xn = optimState.Xn+1;
    idx = optimState.Xn;

    % Expand storage by 50% at a time
    if optimState.Xn > optimState.Xmax
        optimState.Xmax = ceil(optimState.Xmax*1.5);
        optimState.X_orig = expand(optimState.X_orig,optimState.Xmax);
        optimState.X = expand(optimState.X,optimState.Xmax);
        optimState.y_orig = expand(optimState.y_orig,optimState.Xmax);
        optimState.y = expand(optimState.y,optimState.Xmax);
        if ~isempty(fsd); optimState.S = expand(optimState.S,optimState.Xmax); end
        optimState.X_flag = expand(optimState.X_flag,optimState.Xmax,false);
        optimState.funevaltime = expand(optimState.funevaltime,optimState.Xmax);
        optimState.nevals = expand(optimState.nevals,optimState.Xmax);
    end

    optimState.X_orig(optimState.Xn,:) = x_orig;
    optimState.X(optimState.Xn,:) = x;
    optimState.y_orig(optimState.Xn) = fval_orig;
    fval = fval_orig + warpvars_vbmc(x,'logp',optimState.trinfo)/T;
    optimState.y(optimState.Xn) = fval;
    if ~isempty(fsd); optimState.S(optimState.Xn) = fsd; end
    optimState.X_flag(optimState.Xn) = true;
    optimState.funevaltime(optimState.Xn) = t;
    optimState.nevals(optimState.Xn) = max(1,optimState.nevals(optimState.Xn) + 1);
end

optimState.ymax = max(optimState.y(optimState.X_flag));
optimState.N = optimState.Xn;  % Number of training inputs
optimState.Neff = sum(optimState.nevals(optimState.X_flag));

end

%--------------------------------------------------------------------------
function x = expand(x,nmax,fill)
%EXPAND Expand storage for given matrix
if nargin < 3 || isempty(fill); fill = NaN; end

if isnan(fill)
    x = [x; NaN(nmax-size(x,1),size(x,2))];
elseif islogical(fill) && fill
    x = [x; true(nmax-size(x,1),size(x,2))];
elseif islogical(fill) && ~fill
    x = [x; false(nmax-size(x,1),size(x,2))];
else
    x = [x; fill*ones(nmax-size(x,1),size(x,2))];
end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
% Extended Kernel Recursive Least Squares algorithm
%
% W. Liu, I.M. Park. Y. Wang, and J.C. Principe, "Extended Kernel Recursive
% Least Squares Algorithm," IEEE Transactions on Signal Processing, vol.
% 57, no. 10, pp. 3801-3814, Oct. 2009,
% http://dx.doi.org/10.1109/TSP.2009.2022007
%
% Remark: implementation of the tracking model, includes a maximum
% dictionary size M
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

classdef exkrls < kernel_adaptive_filter

    properties (GetAccess = 'public', SetAccess = 'private')
        alphaf = .999; % state forgetting factor, "alpha" in publication
        beta = .995; % data forgetting factor
        lambda = 1E-2; % regularization
        q = 1E-3; % trade-off between modeling variation and measurement disturbance
        M = 500; % maximum dictionary size
        kerneltype = 'gauss'; % kernel type
        kernelpar = 1; % kernel parameter
    end

    properties (GetAccess = 'public', SetAccess = 'private')
        dict = []; % dictionary
        rho = [];
        Q = [];
        i = 0; % iteration number;
        alpha = []; % expansion coefficients, "a" in publication
    end

    methods
        function kaf = exkrls(parameters) % constructor
            allpars = {'alphaf','lambda','beta','q','kerneltype','kernelpar','M'};
            if (nargin > 0)
                for j=1:length(allpars)
                    p = allpars{j};
                    if isfield(parameters,p), kaf.(p) = parameters.(p); end
                end
            end
        end

        function y_est = evaluate(kaf,x) % evaluate the algorithm
            if size(kaf.dict,1)>0
                k = kernel(kaf.dict,x,kaf.kerneltype,kaf.kernelpar);
                y_est = k'*kaf.alpha;
            else
                y_est = zeros(size(x,1),1);
            end
        end

        function train(kaf,x,y) % train the algorithm
            kaf.i = kaf.i + 1;
            k = kernel([kaf.dict; x],x,kaf.kerneltype,kaf.kernelpar);
            kt = k(1:end-1);
            ktt = k(end);
            if numel(kt)==0 % initialize
                kaf.alpha = kaf.alphaf*y/(kaf.lambda*kaf.beta+ktt);
                kaf.rho = kaf.lambda*kaf.beta/(kaf.alphaf^2*kaf.beta + kaf.lambda*kaf.q);
                kaf.Q = kaf.alphaf^2/((kaf.beta*kaf.lambda+ktt)*(kaf.alphaf^2+kaf.beta*kaf.lambda*kaf.q));
                kaf.dict = x;
            else
                if (size(kaf.dict,1)<kaf.M) % avoid infinite growth
                    z = kaf.Q*kt;
                    r = kaf.beta^kaf.i*kaf.rho + ktt - kt'*z;
                    err = y - kt'*kaf.alpha;

                    kaf.alpha = kaf.alphaf*[kaf.alpha - z*err/r; err/r]; % grow
                    kaf.dict = [kaf.dict; x];
                    dummy = kaf.alphaf^2 + kaf.beta^kaf.i*kaf.q*kaf.rho;
                    kaf.rho = kaf.rho/dummy;
                    kaf.Q = kaf.alphaf^2/(r*dummy)*...
                        [kaf.Q*r + z*z', -z; -z', 1];
                end
            end
        end
    end
end

repo: kafbox


text:
function [out1,out2,out3] = ginput_OH(arg1,strpointertype,xx,yy)

% function [X,Y,BUTTON] = ginput_OH(N,pointertype,xx,yy)
%
% ginput_OH.m has the same functionality as Matlab's ginput.m with
% two additional features:
% 1. From the first click onward the pointer is connected to the previous
%    mouse click by a line.
% 2. Selectable cursor pointer (copied from MYGINPUT by F. Moisy)
%
% INPUT:
% N = number of points
% pointertype: 'fullcrosshair','crosshair','arrow','circle'etc.
%           See "Specifying the Figure Pointer" in Matlab's documentation
%           to see the list of available pointers. (default = 'fullcrosshair')
% xx      : reference for line of first point (used by clickfit_OH)
% yy      : reference for line of first point (used by clickfit_OH)
%
% OUTPUT:
% X     : x-coordinates of the clicked series
% Y     : y-coordinates of the clicked series
% BUTTON : see GINPUT
%
% USES:
% ginput.m rev.  5.32.4.13 (M2009a)
% gline.m rev.  2.12.2.5 (M2009a)
% MYGINPUT by F. Moisy (../fileexchange/12770-myginput)
%
% USED by:
% clickfit_OH.m
%
% Author:   Oscar Hartogensis (oscar.hartogensis at wur.nl)
% Date:     Januari 2010
% Rev:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Adjusted by MYGINPUT and ginput_OH
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<2 strpointertype='fullcrosshair'; end
if nargin<4 xx=[]; yy=[]; end
first_run = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

out1 = []; out2 = []; out3 = []; y = [];
c = computer;
if ~strcmp(c(1:2),'PC')
   tp = get(0,'TerminalProtocol');
else
   tp = 'micro';
end

if ~strcmp(tp,'none') && ~strcmp(tp,'x') && ~strcmp(tp,'micro'),
   if nargout == 1,
      if nargin == 1,
         out1 = trmginput(arg1);
      else
         out1 = trmginput;
      end
   elseif nargout == 2 || nargout == 0,
      if nargin == 1,
         [out1,out2] = trmginput(arg1);
      else
         [out1,out2] = trmginput;
      end
      if  nargout == 0
         out1 = [ out1 out2 ];
      end
   elseif nargout == 3,
      if nargin == 1,
         [out1,out2,out3] = trmginput(arg1);
      else
         [out1,out2,out3] = trmginput;
      end
   end
else

   fig = gcf;
   figure(gcf);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Adjusted by ginput_OH
    % Copy-Paste from gline.m rev.  2.12.2.5 (M2009a)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    draw_fig = gcf;
    ax = get(draw_fig,'CurrentAxes');
    if isempty(ax)
       ax = axes('Parent',draw_fig);
    end

    oldtag = get(draw_fig,'Tag');
    figure(draw_fig);
    if any(get(ax,'view')~=[0 90]),
     set(draw_fig, 'Tag', oldtag);
     error('stats:gline:NotTwoDimensional','GLINE works only for 2-D plots.');
    end

    % Can't do this if another mouse mode is active
    hManager = uigetmodemanager(draw_fig);
    if (~isempty(hManager.CurrentMode))
       error('stats:gline:ModeActive',...
             'Cannot use GLINE while the mouse mode ''%s'' is active.',...
             get(hManager.CurrentMode,'name'))
    end

    % Create an invisible line to preallocate the xor line color.
    xlimits = get(ax,'Xlim');
    x_gl = (xlimits + flipud(xlimits))./2;
    ylimits = get(ax,'Ylim');
    y_gl = (ylimits + flipud(ylimits))./2;
    hline = line(x_gl,y_gl,'Parent',ax,'Visible','off');

    bdown = get(draw_fig,'WindowButtonDownFcn');
    bup = get(draw_fig,'WindowButtonUpFcn');
    bmotion = get(draw_fig,'WindowButtonMotionFcn');
    oldud = get(draw_fig,'UserData');
    %    set(draw_fig,'WindowButtonDownFcn','gline(''down'')')
    %    set(draw_fig,'WindowButtonMotionFcn','gline(''motion'')')
    %    set(draw_fig,'WindowButtonupFcn','')

    ud.hline = hline;
    ud.pts = [xx,yy;xx,yy];
    ud.buttonfcn = {bdown; bup; bmotion};
    ud.oldud = oldud;
    ud.oldtag = oldtag;
    set(draw_fig,'UserData',ud);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   if nargin == 0
      how_many = -1;
      b = [];
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Adjusted by ginput_OH
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   elseif isempty(arg1)
      how_many = -1;
      b = [];
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   else
      how_many = arg1;
      b = [];
      if  ischar(how_many) ...
            || size(how_many,1) ~= 1 || size(how_many,2) ~= 1 ...
            || ~(fix(how_many) == how_many) ...
            || how_many < 0
         error('MATLAB:ginput:NeedPositiveInt', 'Requires a positive integer.')
      end
      if how_many == 0
         ptr_fig = 0;
         while(ptr_fig ~= fig)
            ptr_fig = get(0,'PointerWindow');
         end
         scrn_pt = get(0,'PointerLocation');
         loc = get(fig,'Position');
         pt = [scrn_pt(1) - loc(1), scrn_pt(2) - loc(2)];
         out1 = pt(1); y = pt(2);
      elseif how_many < 0
         error('MATLAB:ginput:InvalidArgument', 'Argument must be a positive integer.')
      end
   end

   % Suspend figure functions
   state = uisuspend(fig);

   toolbar = findobj(allchild(fig),'flat','Type','uitoolbar');
   if ~isempty(toolbar)
        ptButtons = [uigettool(toolbar,'Plottools.PlottoolsOff'), ...
                     uigettool(toolbar,'Plottools.PlottoolsOn')];
        ptState = get (ptButtons,'Enable');
        set (ptButtons,'Enable','off');
   end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Adjusted by MYGINPUT
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    set(fig,'pointer',strpointertype);  % modified MYGINPUT
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Adjusted by ginput_OH
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isempty(ud.pts)
        set(draw_fig,'WindowButtonMotionFcn',@motion)
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   fig_units = get(fig,'units');
   char = 0;

   % We need to pump the event queue on unix
   % before calling WAITFORBUTTONPRESS
   drawnow

   while how_many ~= 0

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % Adjusted by ginput_OH
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      if ~first_run
          ud.pts = [pt(1,1),pt(1,2);pt(1,1),pt(1,2)];
          set(draw_fig,'WindowButtonMotionFcn',@motion)
      end
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      % Use no-side effect WAITFORBUTTONPRESS
      waserr = 0;
      try
	keydown = wfbp;
      catch
	waserr = 1;
      end
      if(waserr == 1)
         if(ishghandle(fig))
            set(fig,'units',fig_units);
	    uirestore(state);
            error('MATLAB:ginput:Interrupted', 'Interrupted');
         else
            error('MATLAB:ginput:FigureDeletionPause', 'Interrupted by figure deletion');
         end
      end
        % g467403 - ginput failed to discern clicks/keypresses on the figure it was
        % registered to operate on and any other open figures whose handle
        % visibility were set to off
        figchildren = allchild(0);
        if ~isempty(figchildren)
            ptr_fig = figchildren(1);
        else
            error('MATLAB:ginput:FigureUnavailable','No figure available to process a mouse/key event');
        end
%         old code -> ptr_fig = get(0,'CurrentFigure'); Fails when the
%         clicked figure has handlevisibility set to callback
      if(ptr_fig == fig)
         if keydown
            char = get(fig, 'CurrentCharacter');
            button = abs(get(fig, 'CurrentCharacter'));
            scrn_pt = get(0, 'PointerLocation');
            set(fig,'units','pixels')
            loc = get(fig, 'Position');
            % We need to compensate for an off-by-one error:
            pt = [scrn_pt(1) - loc(1) + 1, scrn_pt(2) - loc(2) + 1];
            set(fig,'CurrentPoint',pt);
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Adjusted by ginput_OH
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            if ~first_run && ~isempty(ud.pts)
                set(ud.hline,'Xdata',ud.pts(:,1),'Ydata',ud.pts(:,2), 'Visible','off');
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         else
            button = get(fig, 'SelectionType');
            if strcmp(button,'open')
               button = 1;
            elseif strcmp(button,'normal')
               button = 1;
            elseif strcmp(button,'extend')
               button = 2;
            elseif strcmp(button,'alt')
               button = 3;
            else
               error('MATLAB:ginput:InvalidSelection', 'Invalid mouse selection.')
            end
         end
         pt = get(gca, 'CurrentPoint');

         how_many = how_many - 1;

         if(char == 13) % & how_many ~= 0)
            % if the return key was pressed, char will == 13,
            % and that's our signal to break out of here whether
            % or not we have collected all the requested data
            % points.
            % If this was an early breakout, don't include
            % the <Return> key info in the return arrays.
            % We will no longer count it if it's the last input.
            break;
         end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Adjusted by ginput_OH
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        first_run = 0;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

         out1 = [out1;pt(1,1)];
         y = [y;pt(1,2)];
         b = [b;button];
      end
   end

   uirestore(state);
   if ~isempty(toolbar) && ~isempty(ptButtons)
        set (ptButtons(1),'Enable',ptState{1});
        set (ptButtons(2),'Enable',ptState{2});
   end
   set(fig,'units',fig_units);

   if nargout > 1
      out2 = y;
      if nargout > 2
         out3 = b;
      end
   else
      out1 = [out1 y];
   end

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Adjusted by ginput_OH
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if ~first_run && ~isempty(ud.pts)
       set(ud.hline,'Xdata',ud.pts(:,1),'Ydata',ud.pts(:,2), 'Visible','off');
   end
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Adjusted by ginput_OH
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function motion(varargin)
   % Copy-Paste from gline.m (case 'motion') rev.  2.12.2.5 (M2009a)
   set(ud.hline,'Xdata',ud.pts(:,1),'Ydata',ud.pts(:,2), ...
        'linestyle','-', 'linewidth', 1.5, 'Color','r','Visible','on');
   Pt2 = get(ax,'CurrentPoint');
   Pt2 = Pt2(1,1:2);
   ud.pts(2,:) = Pt2;
   set(draw_fig,'UserData',ud);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function key = wfbp
%WFBP   Replacement for WAITFORBUTTONPRESS that has no side effects.

fig = gcf;
current_char = [];

% Now wait for that buttonpress, and check for error conditions
waserr = 0;
try
  h=findall(fig,'type','uimenu','accel','C');   % Disabling ^C for edit menu so the only ^C is for
  set(h,'accel','');                            % interrupting the function.
  keydown = waitforbuttonpress;
  current_char = double(get(fig,'CurrentCharacter')); % Capturing the character.
  if~isempty(current_char) && (keydown == 1)           % If the character was generated by the
	  if(current_char == 3)                       % current keypress AND is ^C, set 'waserr'to 1
		  waserr = 1;                             % so that it errors out.
	  end
  end

  set(h,'accel','C');                                 % Set back the accelerator for edit menu.
catch
  waserr = 1;
end
drawnow;
if(waserr == 1)
   set(h,'accel','C');                                % Set back the accelerator if it errored out.
   error('MATLAB:ginput:Interrupted', 'Interrupted');
end

if nargout>0, key = keydown; end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end

end

author: Sable

repo: mcbench-benchmarks


text:
function deg = degree(x,y,z)
%DEGREE Polynomial degree
%
% DEG = DEGREE(p,x,flag,vector)
%
% p      : SDPVAR object.
% x      : Degree w.r.t linear SDPVAR objects.
% flag   : 'max', 'min'. Default 'max'
% vector : If vector = 1, returns degree of each element in p
%
% Examples
% x1 = sdpvar(1,1);x2 = sdpvar(1,1);
% p = [x1;x1*x2+x2^2];
%
% degree(p) returns 2
%
% degree(p,x1) returns 1
%
% degree(p,[x1 x2]) returns [1 2]
%
% degree(p,[x1 x2],'max',1) returns [1 0;1 2]
%
% degree(p,[],1) returns [1;2]

if isa(x,'double')
    deg = 0;
else
    error('DEGREE only defined for SDPVAR and DOUBLE.');
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
filename='Gripping_tetrahedra_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'nonadjoint_compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
constraint_case = 'INEQUALITY';
optimizer = 'IPOPT'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.2;
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:
filename='RVE_Square_Triangle';
ptype = 'MICRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'circleInclusion';
cost={'enforceCh_CCstar_L2','perimeter'};%enforceCh_CCstar_L2
weights=[1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'MMA'; incrementFactor = 1;%mma CON 1E-4
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.7;
Perimeter_target=3;
optimality_final =1e-4;
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;

%Micro
epsilon_isotropy_initial=1e-1;
epsilon_isotropy_final = 1e-3;
selectiveC_Cstar = 'Vfrac07';


text:
function line = extractAreaLine(area_data)

% EXTRACTS THE LINE STYLE USED FOR MATLAB OBJECTS
% OF TYPE "LINE". THESE OBJECTS ARE USED IN LINESERIES,
% STAIRSERIES, STEMSERIES, BASELINESERIES, AND BOXPLOTS

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

%-INITIALIZE OUTPUT-%
line = struct();

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

%-AREA LINE COLOR-%

if(~strcmp(area_data.LineStyle,'none'))

    % marker edge color
    LineColor = area_data.EdgeColor;

    if isnumeric(LineColor)
        col = 255*LineColor;
        linecolor = ['rgb(' num2str(col(1)) ',' num2str(col(2)) ',' num2str(col(3)) ')'];
    else
        linecolor = 'rgba(0,0,0,0)';
    end

    line.color = linecolor;

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

    %-PATCH LINE WIDTH (STYLE)-%
    line.width = area_data.LineWidth;

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

    %-PATCH LINE DASH (STYLE)-%
    switch area_data.LineStyle
        case '-'
            LineStyle = 'solid';
        case '--'
            LineStyle = 'dash';
        case ':'
            LineStyle = 'dot';
        case '-.'
            LineStyle = 'dashdot';
    end

    line.dash = LineStyle;

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

end
end

author: plotly


text:
function test_plotMeanTSeries
%Validate getting plotting the mean TSeries from an INPLANE ROI
%
%  test_plotMeanTSeries()
%
% Tests: meanTSeries, plotMeanTSeries
%
% INPUTS
%  No inputs
%
% RETURNS
%  No returns
%
% Example: test_plotMeanTSeries()
%
% See also MRVTEST
%
% Copyright Stanford team, mrVista, 2011

%% Initialize the key variables and data path
% Data directory (where the mrSession file is located)
dataDir = mrtInstallSampleData('functional','mrBOLD_01');

% This is the validation file
val = mrtGetValididationData('plotMeanTSeriesFromINPLANE');

% These are the items we stored in the validation file
%
% val.detrendDim  = size(d.detrend.frameNumbers);
% val.rawDim      = size(d.raw.frameNumbers);
% val.detrendMn   = mean(d.detrend.tSeries);
% val.detrendMd   = median(d.detrend.tSeries);
% val.detrendMx   = max(d.detrend.tSeries);
% val.detrendMin  = min(d.detrend.tSeries);
% val.rawMn       = mean(d.raw.tSeries);
% val.rawMd       = median(d.raw.tSeries);
% val.rawMx       = max(d.raw.tSeries);
% val.rawMin      = min(d.raw.tSeries);
%
% save(vFile, '-struct', 'val')

%% Retain original directory, change to data directory
curDir = pwd;
cd(dataDir);

% There can be several data types - name the one you want to probe
dataType = 'Original';

% Which scan number from that data type?
scan = 1;

%% Get data structure:
vw = initHiddenInplane(); % Foregoes interface - loads data silently

%% Set dataTYPE:
vw = viewSet(vw, 'Current DataType', dataType); % Data type

%% Load an ROI and coranal
vw = loadROI(vw, 'LV1.mat');

detrend = true;

% open a plot figure
newGraphWin;

% load both raw and detrended tSeries so we can validate both
d.detrend = plotMeanTSeries(vw, scan, [], ~detrend);
d.raw     = plotMeanTSeries(vw, scan, [], detrend);

% close it
closeGraphWin;

%% Go home
cd(curDir)

%% Validate..

% check the number of time points
assertEqual(val.detrendDim, size(d.detrend.frameNumbers));
assertEqual(val.rawDim,     size(d.raw.frameNumbers));

% check detrended t-series
assertElementsAlmostEqual(val.detrendMn,mean(d.detrend.tSeries));
assertElementsAlmostEqual(val.detrendMd,median(d.detrend.tSeries));
assertElementsAlmostEqual(val.detrendMx,max(d.detrend.tSeries));
assertElementsAlmostEqual(val.detrendMin,min(d.detrend.tSeries));

% check raw t-series
assertElementsAlmostEqual(val.rawMn,mean(d.raw.tSeries));
assertElementsAlmostEqual(val.rawMd,median(d.raw.tSeries));
assertElementsAlmostEqual(val.rawMx,max(d.raw.tSeries));
assertElementsAlmostEqual(val.rawMin,min(d.raw.tSeries));

%% Clear workspace

mrvCleanWorkspace;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [pvec] = real_mutate(pvec, pmut_real, eta_m, ...
                                min_realvar, max_realvar)
%   This code applies polynomical mutation over pvec.

nreal = length(pvec);

if(nreal < 10)
    pvec = polymut_looped(pvec, pmut_real, eta_m, ...
                            min_realvar, max_realvar);
else
    pvec = polymut_vectorized(pvec, pmut_real, eta_m, ...
                            min_realvar, max_realvar);
end
end

function [pvec] = polymut_looped(pvec, pmut_real, eta_m, ...
                                    min_realvar, max_realvar)
% This is the original translation from the nsga-2 c-code.
nreal = length(pvec);
for j = 1:nreal
    if (rand(1) <= pmut_real)
    % if (randomperc() <= pmut_real) % SLOW !!!
        y = pvec(j);
        yl = min_realvar(j);
        yu = max_realvar(j);
        delta1 = (y-yl)/(yu-yl);
        delta2 = (yu-y)/(yu-yl);
        mut_pow = 1.0/(eta_m + 1.0);
        r = rand(1) ;
        % r = randomperc() ; % SLOW !!!
        if (r <= 0.5)
            xy = 1.0 - delta1;
            val = 2.0 * r + (1.0 - 2.0 * r) * (xy ^ (eta_m + 1.0));
            deltaq = (val ^ mut_pow) - 1.0;
        else
            xy = 1.0 - delta2;
            val = 2.0 * (1.0 - r) + 2.0 * (r - 0.5) * (xy ^ (eta_m + 1.0));
            deltaq = 1.0 - (val ^ mut_pow);
        end
        y = y + deltaq * (yu - yl);
        if (y < yl)
            y = yl;
        end
        if (y > yu)
            y = yu;
        end
        pvec(j) = y ;
    end
end
end

function [pvec] = polymut_vectorized(pvec, pmut_real, eta_m, ...
                                        min_realvar, max_realvar)
% This is the vectorized version of the above code. This code is generally
% 3 times faster than the above, more gain could be observed if the number
% of variable is bigger and the mutation rate is higher.

nreal = length(pvec);
mut_index = rand(1,nreal) < pmut_real ;
% mut_index = randompercv(1,nreal) < pmut_real ; % SLOW !!!
abs_mut_index = 1:nreal ;
abs_mut_index = abs_mut_index(mut_index);
mlen = length(abs_mut_index);

if(mlen > 0)
    eta_mv = ones(1,mlen) * eta_m ;
    yv = pvec(mut_index) ;
    ylv = min_realvar(mut_index).';
    yuv = max_realvar(mut_index).';
    delta1v = (yv - ylv) ./ (yuv - ylv);
    delta2v = (yuv - yv) ./ (yuv - ylv);
    mut_powv = 1.0 ./ (eta_mv + 1.0);
    rv = rand(1, mlen); % r2s ;
    % rv = randompercv(1, mlen); % r2s ; % SLOW !!!
    rvlthalf = rv < 0.5 ;
    valv1 = 2.0 .* rv + (1.0 - 2.0 .* rv) .* ...
                        ((1.0 - delta1v) .^ (eta_mv + 1.0));
    valv2 = 2.0 .* (1.0 - rv) + 2.0 .* (rv - 0.5) .* ...
                        ((1.0 - delta2v) .^ (eta_mv + 1.0));
    deltaqv = rvlthalf .* ((valv1 .^ mut_powv) - 1.0) + ...
             (~rvlthalf) .* (1.0 - (valv2 .^ mut_powv));
    yv = yv + deltaqv .* (yuv - ylv);
    yltyl = yv < ylv ;
    ygtyu = yv > yuv ;
    yv = (yltyl .* ylv) + (yv .* (~yltyl));
    yv = (ygtyu .* yuv) + (yv .* (~ygtyu));
    pvec = pvec' ;
    pvec(abs_mut_index.',:) = yv' ;
    pvec = pvec' ;
end
end

author: chudur-budur


text:
function [y_cp] = SpanLoading(l_s, l_t, CL, geo, panel)
%Determine center of pressure and spanwise loading
%due to bound vortices and trailing vortices.

BV = [panel.BV];
BV1 = [panel.BV1];
s = [panel.s]';

y_cp.alpha = sum([l_s.alpha].*BV(2,:)' + [l_t.alpha].*BV1(2,:)')/(0.5*CL.alpha*0.5*geo.b);  %Eqn 35 in NASA paper
cl.alpha = ([l_s.alpha] + [l_t.alpha])*geo.S./(CL.alpha*2*s*cos(geo.dih)*geo.c_av);  %Eqn 37 in NASA paper

if CL.c ~= 0  %Prevent divide by zero
y_cp.c = sum([l_s.c].*BV(2,:)' + [l_t.c].*BV1(2,:)')/(0.5*CL.c*0.5*geo.b);  %Eqn 35 in NASA paper
cl.c = ([l_s.c] + [l_t.c])*geo.S./(CL.c*2*s*cos(geo.dih)*geo.c_av);  %Eqn 37 in NASA paper
else
    y_cp.c = 0;
    cl.c = zeros(size(cl.alpha));
end


text:
function c=zak(f,a);
%ZAK  Zak transform
%   Usage:  c=zak(f,a);
%
%   `zak(f,a)` computes the Zak transform of *f* with parameter *a*.  The
%   coefficients are arranged in an $a \times L/a$ matrix, where *L* is the
%   length of *f*.
%
%   If *f* is a matrix then the transformation is applied to each column.
%   This is then indexed by the third dimension of the output.
%
%   Assume that $c=zak(f,a)$, where *f* is a column vector of length *L* and
%   $N=L/a$. Then the following holds for $m=0,\ldots,a-1$ and $n=0,\ldots,N-1$
%
%   ..                     N-1
%     c(m+1,n+1)=1/sqrt(N)*sum f(m-k*a+1)*exp(2*pi*i*n*k/N)
%                          k=0
%
%   .. math:: c(m+1,n+1) = \frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}f(m-ka+1)e^{2\pi ink/M}
%
%   Examples:
%   ---------
%
%   This figure shows the absolute value of the Zak-transform of a Gaussian.
%   Notice that the Zak-transform is 0 in only a single point right in the
%   middle of the plot :::
%
%     a=64;
%     L=a^2;
%     g=pgauss(L);
%     zg=zak(g,a);
%
%     surf(abs(zg));
%
%   This figure shows the absolute value of the Zak-transform of a 4th order
%   Hermite function.  Notice how the Zak transform of the Hermite functions
%   is zero on a circle centered on the corner :::
%
%     a=64;
%     L=a^2;
%     g=pherm(L,4);
%     zg=zak(g,a);
%
%     surf(abs(zg));
%
%   See also:  izak
%
%   References: ja94-4 bohl97-1

%   AUTHOR : Peter L. Søndergaard
%   TESTING: TEST_ZAK
%   REFERENCE: REF_ZAK

complainif_argnonotinrange(nargin,2,2,mfilename);

if (prod(size(a))~=1 || ~isnumeric(a))
  error([callfun,': a must be a scalar']);
end;

if rem(a,1)~=0
  error([callfun,': a must be an integer']);
end;

if size(f,2)>1 && size(f,1)==1
  % f was a row vector.
  f=f(:);
end;

L=size(f,1);
W=size(f,2);
N=L/a;

if rem(N,1)~=0
  error('The parameter for ZAK must divide the length of the signal.');
end;

c=zeros(a,N,W,assert_classname(f));

for ii=1:W
  % Compute it, it can be done in one line!
  % We use a normalized DFT, as this gives the correct normalization
  % of the Zak transform.
  c(:,:,ii)=dft(reshape(f(:,ii),a,N),[],2);
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
filename='Gripping_triangle_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'nonadjoint_compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
constraint_case = 'INEQUALITY';
optimizer = 'IPOPT'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 1;
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;


text:
function test_bug1082

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_prepare_vol_sens ft_prepare_headmodel ft_compute_leadfield ft_plot_topo3d ft_headmodel_openmeeg

fprintf('***  DIAGNOSTICAL INFORMATION ***\n');
fprintf('test script is running on host: %s\n', gethostname());

%%

% generate a unit sphere
[pnt, tri] = mesh_sphere(162);

% create the BEM geometries (in mm)
bnd = [];
bnd.pnt = pnt * 100;
bnd.tri = tri; % normals outwards
% bnd.tri = fliplr(tri); % normals inwards

% create a set of electrodes
sel = find(pnt(:,3)>0);
sens = [];
sens.chanpos = pnt(sel,:) * 100;
sens.elecpos = sens.chanpos;
for i=1:length(sel)
  sens.label{i} = sprintf('chan%03d', i);
end

% this is the position of the dipole
pos = [0 0 50];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% calculate volume conductors

tmpcfg = [];
tmpcfg.conductivity = [1];
tmpcfg.method = 'singlesphere';
vol1 = ft_prepare_headmodel(tmpcfg,bnd);
[vol1, sens1] = ft_prepare_vol_sens(vol1, sens);
lf1 = ft_compute_leadfield(pos, sens1, vol1);

tmpcfg = [];
tmpcfg.conductivity = [1];
tmpcfg.method = 'openmeeg';
vol2 = ft_prepare_headmodel(tmpcfg,bnd);
[vol2, sens2] = ft_prepare_vol_sens(vol2, sens);
lf2 = ft_compute_leadfield(pos, sens2, vol2);

tmpcfg = [];
tmpcfg.conductivity = [1];
tmpcfg.method = 'dipoli';
vol3 = ft_prepare_headmodel(tmpcfg,bnd);
[vol3, sens3] = ft_prepare_vol_sens(vol3, sens);
lf3 = ft_compute_leadfield(pos, sens3, vol3);

%%

figure;
subplot(2,2,1); ft_plot_topo3d(sens1.chanpos, lf1(:,1))
subplot(2,2,2); ft_plot_topo3d(sens1.chanpos, lf1(:,2))
subplot(2,2,3); ft_plot_topo3d(sens1.chanpos, lf1(:,3))
colorbar

figure;
subplot(2,2,1); ft_plot_topo3d(sens2.chanpos, lf2(:,1))
subplot(2,2,2); ft_plot_topo3d(sens2.chanpos, lf2(:,2))
subplot(2,2,3); ft_plot_topo3d(sens2.chanpos, lf2(:,3))
colorbar

figure;
subplot(2,2,1); ft_plot_topo3d(sens3.chanpos, lf3(:,1))
subplot(2,2,2); ft_plot_topo3d(sens3.chanpos, lf3(:,2))
subplot(2,2,3); ft_plot_topo3d(sens3.chanpos, lf3(:,3))
colorbar

% figure;
% subplot(2,2,1); ft_plot_topo3d(senst.chanpos, lft(:,1))
% subplot(2,2,2); ft_plot_topo3d(senst.chanpos, lft(:,2))
% subplot(2,2,3); ft_plot_topo3d(senst.chanpos, lft(:,3))
% colorbar

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% readfid allows one to read in the Varian 4T FID data acquired with:
%
% [K, HDR, BLOCK_HDR] = READFID(FPATH, PAR, DC_CORRECT) reads k-space data
% from fid file if the FID directory FPATH. DC correction based on level
% and tilt can be can be suppressed by setting DC_CORRECT to false. If DC
% correction is applied, the level and tilt are set to zero. PAR is the
% parameters as read by getparam. The main header is returned in HDR and
% the block headers are returned in BLOCK_HDR
%
% see also WRITEFID, READPROCPAR, WRITEPROCPAR

function [k, hdr, block_hdr, par] = readfid(fpath, par, dc_correct, tfisp, ext)

if nargin < 1
   error('cfmm:fid:name','READFID requires fid file name');
end

if nargin < 5
   ext = [];
end

% Read procpar if not already provided
if nargin < 2 || isempty(par)
   par = readprocpar(fpath, ext);
end

% Turn on DC correction by default
if nargin < 3 || isempty(dc_correct)
   dc_correct = true;
end
% Turn off tfisp by default
if nargin < 4 || isempty(tfisp)
   tfisp = false;
end

% Check for validity of fid file
if exist([fpath '/fid' ext], 'file') ~= 2
   if exist([fpath '.fid/fid'], 'file') ~= 2
      error('cfmm:fid:name','Unable to open FID file: %s', fpath);
   else
      fpath = [fpath '.fid'];
   end
end

infid = fopen([fpath '/fid' ext],'r', 'ieee-be');
if infid == -1
   error('cfmm:fid:name','Unable to open FID file: %s', fpath);
end

% Read in all the information from the FID file
try
   hdr.nblocks    = fread(infid, 1, 'int32');
   hdr.ntraces    = fread(infid, 1, 'int32');
   hdr.np         = fread(infid, 1, 'int32');
   hdr.ebytes     = fread(infid, 1, 'int32');
   hdr.tbytes     = fread(infid, 1, 'int32');
   hdr.bbytes     = fread(infid, 1, 'int32');
   hdr.vers_id    = fread(infid, 1, 'int16');
   hdr.status     = fread(infid, 1, 'int16');
   hdr.nbheaders  = fread(infid, 1, 'int32');

   m = hdr.bbytes/hdr.ebytes;
   n = hdr.nblocks;

   % Read in the block headers
   % Imaging data all has only one block header and VERY BAD people use
   % nbheaders as their own spare field. Therefore we ignore the nbheaders
   % and treat everything as having one block header.
   if nargout > 2
      % Create the block_hdr structure
      block_hdr.scale   = 0;
      block_hdr.status  = 0;
      block_hdr.index   = 0;
      block_hdr.mode    = 0;
      block_hdr.ctcount = 0;
      block_hdr.lpval   = 0;
      block_hdr.rpval   = 0;
      block_hdr.lvl     = 0;
      block_hdr.tlt     = 0;
      block_hdr = block_hdr(ones(n, 1));

      % Read in the block header values
      for i = 1:n
         block_hdr(i).scale   = fread(infid, 1, 'short');
         block_hdr(i).status  = fread(infid, 1, 'short');
         block_hdr(i).index   = fread(infid, 1, 'short');
         block_hdr(i).mode    = fread(infid, 1, 'short');
         block_hdr(i).ctcount = fread(infid, 1, 'long');
         block_hdr(i).lpval   = fread(infid, 1, 'float');
         block_hdr(i).rpval   = fread(infid, 1, 'float');
         block_hdr(i).lvl     = fread(infid, 1, 'float');
         block_hdr(i).tlt     = fread(infid, 1, 'float');
         fseek(infid, hdr.bbytes-28*hdr.nbheaders, 'cof');
      end
   end

   % Get the data format from the header
   if bitand(hdr.status, 8)
      % Floating
      format = '*float32';
      if hdr.ebytes ~= 4
         error('cfmm:fid:readError','Inconsistent Format: %d-bytes and %s.', hdr.ebytes, format);
      end
   else
      % Integer
      if bitand(hdr.status, 4)
         % 32-bit integer
         format = '*int32';
         if hdr.ebytes ~= 4
            error('cfmm:fid:readError','Inconsistent Format: %d-bytes and %s.', hdr.ebytes, format);
         end
      else
         % 16-bit integer
         format = '*int16';
         if hdr.ebytes ~= 2
            error('cfmm:fid:readError','Inconsistent Format: %d-bytes and %s.', hdr.ebytes, format);
         end
      end
   end

   % Read in all the data
   fseek(infid, 32, 'bof');
   k = fread(infid, [m,n], format);

   % Strip block header and reorganize into complex values
   k = complex(k(28/hdr.ebytes+1:2:m,:), k(28/hdr.ebytes+2:2:m,:));

   % Apply the lvl and tlt correction
   if dc_correct && rem(par.nt(1),2)
      % Read in the lvl and tlt correction
      fseek(infid, 52, 'bof');
      cor = fread(infid, [2, n], '2*float32', (hdr.bbytes) - 8);

      % Set the block header values to zero after correcting
      if (nargout > 1)
         for i = 1:n
            block_hdr(i).lvl     = 0;
            block_hdr(i).tlt     = 0;
         end
      end

      % nt = 3,7,11,... have DC rotated by 90 degrees
      if rem(par.nt,4) == 3
         cor = complex(-cor(2,:), cor(1,:));
      else
         cor = complex(cor(1,:), cor(2,:));
      end
      if tfisp
         % tfisp switches detection direction on each phase encode
         % acquisition is always compressed phase encodes
         [m, n] = size(k);
         nro = par.np / 2;
         if par.nv == 0
            nv = 1;
         else
            nv = par.nv;
         end
         r = m / (nro * nv);
         k = reshape(k, [nro nv r n]);
         if any(cor ~= cor(1))
            cor = reshape(cor, [1 1 1 n]);
            nv_1 = ceil(nv / 2);
            nv_2 = floor(nv / 2);
            if rem(par.ssc, 2) == 0
               k(:,1:2:end,:,:) = k(:,1:2:end,:,:) + cor(ones(1,nro), ones(1,nv_1), ones(1,r), :);
               k(:,2:2:end,:,:) = k(:,2:2:end,:,:) - cor(ones(1,nro), ones(1,nv_2), ones(1,r), :);
            else
               k(:,2:2:end,:,:) = k(:,2:2:end,:,:) + cor(ones(1,nro), ones(1,nv_2), ones(1,r), :);
               k(:,1:2:end,:,:) = k(:,1:2:end,:,:) - cor(ones(1,nro), ones(1,nv_1), ones(1,r), :);
            end
            % cor = reshape(cor, [1 n]);
         else
            if rem(par.ssc, 2) == 0
               k(:,1:2:end,:) = k(:,1:2:end,:) + cor(1);
               k(:,2:2:end,:) = k(:,2:2:end,:) - cor(1);
            else
               k(:,2:2:end,:) = k(:,2:2:end,:) + cor(1);
               k(:,1:2:end,:) = k(:,1:2:end,:) - cor(1);
            end
         end
         k = reshape(k, [m n]);
      else
         if any(cor ~= cor(1))
            k = k - cor(ones(1,size(k,1)),:);
         else
            k = k - cor(1);
         end
      end
   end

catch err
    fclose(infid);
    rethrow(err);
end
fclose(infid);

author: CIC-methods

repo: FID-A

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


function [prices, stdErrs] = Price_MC_Hindsight_Strikes_func(Spath, call, strikes, M, mult, disc)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% About: Calculates Hindsight option prices for vector of strikes, given the simulatd paths
% Author: Justin Lars Kirkby
%
% -----------------
% Params
% -----------------
% Contract Types
% Fixed Strike (Hindsight) Put: (W - min{S_m: 0<=m<=M})^+
% Fixed Strike (Hindsight) Call: (max{S_m: 0<=m<=M} - W)^+
%
% Spath = paths of underlying, dimension N_sim x M+1, where M = number of time steps (since includes S_0)
% call = 1 for call (else put)
% Kvec = strike vector
% M = number of monitoring points, e.g. 252 for "daily" monitoring
% mult = time partitioning multiplier to reduce bias (e.g. mult = 2 or 5)
% S_0 = initial price
% disc = discount factor (e.g. exp(-r*T))
% T = Time (in years)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N_sim = size(Spath,1);  % number of paths

M_mult = M*mult;  %time partitioning to reduce bias

prices = zeros(length(strikes),1);
stdErrs = zeros(length(strikes),1); %TODO: add code to return std errros

if call == 1
    curr_max = zeros(N_sim, 1);
    for n = 1:N_sim
        curr_max(n) = max(Spath(n, 1:mult:M_mult+1));
    end
else  % find_min
    curr_min = zeros(N_sim, 1);
    for n = 1:N_sim
        curr_min(n) = min(Spath(n, 1:mult:M_mult+1));
    end
end

for k = 1:length(strikes)
    K = strikes(k);
    if call ==1  % Fixed Strike (Hindsight) Call: (max{S_m: 0<=m<=M} - W)^+
       payoffs  = max(0, curr_max - K);

    else  % Fixed Strike (Hindsight) Put: (W - min{S_m: 0<=m<=M})^+
        payoffs  = max(0, K - curr_min);
    end
    prices(k) = disc*mean(payoffs);
    stdErrs(k) = disc*std(payoffs) / sqrt(N_sim);

end

end


text:
% This is the timer function to update the figure.
% Save as timerCallbackFcn.m
function timerCallbackFcn(hTimer, eventData, hfig, I, map)
    figure(hfig);

    % i is the frame index
    persistent i;
    if isempty(i), i=1; end

    % display the i-th frame in the GIF
    imshow(I(:,:,i),map);

    % increment frame index i
    i=i+1;
    numframes=size(I,4);
    if (i>numframes), i=1; end
end

author: kmansouri

repo: OPERA

save_path: github-repos/MATLAB/kmansouri-OPERA


text:
function [tt]=tt_qfromfull(ff,s,d,eps,q)

% Approximates a full-format s-dimensional q^{d} x...x q^{d}-tensor
% full, indexed by
% i^{1}_{1} ,..., i^{d}_{1} ,...,..., i^{1}_{s} ,..., i^{d}_{s},
% by a QTT decomposition _tt_,
% k-th core of which is indexed by
% i^{k}_{1} ,..., i^{k}_{s},
% k=1,...,d
%
% August 11, 2011
% Vladimir Kazeev
% vladimir.kazeev@gmail.com
% INM RAS
% Moscow, Russia
%

if (nargin<4) || isempty (eps)
    eps=1.e-8;
end

if (nargin<5) || isempty (q)
    q=2;
end

ff=reshape(ff,[q*ones(1,s*d),1]);

prm = 1:s*d; prm = reshape(prm,[d,s]); prm = prm'; prm = reshape(prm,[1,d*s]);

ff=permute(ff,[prm,1+s*d]);
ff=reshape(ff,[q^s*ones(1,d),1]);

tt=core(tt_tensor(ff,eps));
if (numel(tt) > d)
	tt{d}=tt{d}*tt{d+1};
	tt=tt(1:d);
end

for k=1:d
	[~,p,r]=size(tt{k});
	tt{k}=reshape(tt{k},[q*ones(1,s),p,r]);
end

return
end

author: oseledets

repo: TT-Toolbox


text:
function value = sphere_area_3d ( r )

%*****************************************************************************80
%
%% SPHERE_AREA_3D computes the area of a sphere in 3D.
%
%  Integration region:
%
%    Points (X,Y,Z) such that
%
%      X**2 + Y**2 + Z**2 = R**2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the sphere.
%
%    Output, real SPHERE_AREA_3D, the area of the sphere.
%
  value = 4.0E+00 * pi * r * r;

  return
end


text:
function out = LegPlot(text,opt)
% =======================================================================
% This function creates a legend outside the chart, centered, & below
% =======================================================================
% out = LegPlot(text,opt)
% -----------------------------------------------------------------------
% INPUT
%   - text: cell, containing the legend text
% -----------------------------------------------------------------------
% OPTIONAL INPUTS
%   - opt: options for the legend, output of LegOption function
% =======================================================================
% Ambrogio Cesa Bianchi, March 2015
% ambrogio.cesabianchi@gmail.com

%% Preliminaries
% =========================================================================
if ~exist('opt','var')
    opt = LegOption;
end

hsize = opt.hsize;
interpreter = opt.interpreter;
direction = opt.direction;
handle = opt.handle;

if hsize<=0 || hsize>1
    disp('Error: hsize must be between zero and one');
    return
end

if isempty(handle)==1
    handle_flag = 0;
else
    handle_flag = 1;
end

%% Legend
% =========================================================================
% Retrieve the position of the figure: [a b c d] = get(gcf,'Position')
% a and b are the (x,y) location of the upper left corner; c and d are the width and height of the window
figpos = get(gcf,'Position');

% Create the legend
if direction==1
    if handle_flag==1
        h1 = legend(handle,text,'Orientation','Horizontal','Fontsize',get(gca,'Fontsize')-1,'Interpreter',interpreter);
    else
        h1 = legend(text,'Orientation','Horizontal','Fontsize',get(gca,'Fontsize')-1,'Interpreter',interpreter);
    end
else
    if handle_flag==1
        h1 = legend(handle,text,'Orientation','Vertical','Fontsize',get(gca,'Fontsize')-1,'Interpreter',interpreter);
    else
        h1 = legend(text,'Orientation','Vertical','Fontsize',get(gca,'Fontsize')-1,'Interpreter',interpreter);
    end
end

% Retrieve the position of the legend:
% out(1) is the location along the x-axis
% out(2) is the location along the y-axis
% out(3) is the width
% out(2) is the height
out = get(h1,'Position');

% ------------- Examples ------------------------------------
% out(1)	out(2)	Vertical alignment	Horizontal alignment
% 0       0           Bottom              Left
% 0.5     0           Bottom              Center
% 0.9     0           Bottom              Right
% 0       0.5         Center              Left
% 0.5     0.5         Center              Center
% 0.9     0.5         Center              Right
% 0       0.9         Top                 Left
% 0.5     0.9         Top                 Center
% 0.9     0.9         Top                 Right
% ------------------------------------------------------------

% Make room for sublegend if font is different from default
axpos = get(gca,'outerPosition');

% % Increase the hsize of the figure to make room for the legend
% FigSize(figpos(3),figpos(4)*(1+out(4)));

% Reduce the y hsize of the plot and move up to make room for legend
axpos(1) = axpos(1);
axpos(2) = axpos(2) + 1.1*out(4);
axpos(3) = axpos(3);
axpos(4) = 1 - 1.1*out(4);
set(gca,'outerPosition',axpos);

% Determine the position for centering the legend
out(3) = hsize; % First set the width of the legend in relative terms (1=whole figure)
out(1) = (figpos(3)-figpos(3)*out(3)) /2 /figpos(3); % Then compute the center
out(2) = 0.005;

% Place the legend
set(h1,'Position',out);

author: ambropo

repo: VAR-Toolbox

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


text:
%ISROT Test if SO(3) rotation matrix
%
% ISROT(R) is true (1) if the argument is of dimension 3x3 or 3x3xN, else false (0).
%
% ISROT(R, 'check') as above, but also checks the validity of the rotation
% matrix.
%
% Notes::
% - A valid rotation matrix has determinant of 1.
%
% See also ISHOMOG, ISROT2, ISVEC.

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math

function h = isrot(R, rtest)

    h = false;
    d = size(R);

    if ndims(R) >= 2
        if ~(all(d(1:2) == [3 3]))
            return %false
        end

        if nargin > 1
            for i = 1:size(R,3)
                RR = R(:,:,i);
                e = RR'*RR - eye(3,3);
                if norm(e) > 10*eps
                    return %false
                end
                e = abs(det(RR) - 1);
                if norm(e) > 10*eps
                    return %false
                end
            end
        end
    end

    h = true;
end

author: petercorke

repo: spatialmath-matlab


text:
% sideways motion of center of rotation
function [data,units] = compute_dv_cor(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  if trx(fly).nframes < 2,
    data{i} = [];
  else
    % center of rotation
    [x_cor_curr,y_cor_curr,x_cor_next,y_cor_next] = rfrac2center(trx,fly,[trx(fly).corfrac_maj;trx(fly).corfrac_min]);
    % change in center of rotation
    dx_cor = x_cor_next - x_cor_curr;
    dy_cor = y_cor_next - y_cor_curr;
    % forward motion of center of rotation
    data{i} = (dx_cor.*cos(trx(fly).theta_mm(1:end-1)+pi/2) + dy_cor.*sin(trx(fly).theta_mm(1:end-1)+pi/2))./trx(fly).dt;

  end

end
units = parseunits('mm/s');

author: kristinbranson

repo: JAABA


text:
%
% LibQPEP: A Library for Globally Optimal Solving Quadratic Pose Estimation Problems (QPEPs),
%          It also gives highly accurate uncertainty description of the solutions.
%
%
% Article:
%      Wu, J., Zheng, Y., Gao, Z., Jiang, Y., Hu, X., Zhu, Y., Jiao, J., Liu, M. (2020)
%           Quadratic Pose Estimation Problems: Unified Solutions,
%           Solvability/Observability Analysis and Uncertainty Description
%           in A Globally Optimal Framework.
%
%
% Authors:      Jin Wu and Ming Liu
% Affiliation:  Hong Kong University of Science and Technology (HKUST)
% Emails:       jin_wu_uestc@hotmail.com; eelium@ust.hk
% Websites:     https://zarathustr.github.io
%               https://ram-lab.com

function R = q2R(q)
q0 = q(1); q1 = q(2); q2 = q(3); q3 = q(4);
R = [
        q0^2 + q1^2 - q2^2 - q3^2,         2*q0*q3 + 2*q1*q2,         2*q1*q3 - 2*q0*q2;
                2*q1*q2 - 2*q0*q3, q0^2 - q1^2 + q2^2 - q3^2,         2*q0*q1 + 2*q2*q3;
                2*q0*q2 + 2*q1*q3,         2*q2*q3 - 2*q0*q1, q0^2 - q1^2 - q2^2 + q3^2];
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
classdef VehicleArticulatedLinear < VehicleDynamicsLateral.VehicleArticulated
    % VehicleArticulatedLinear Linear articulated vehicle model.
    %
    % It inherits properties from VehicleArticulated.

    methods
        % Constructor
        function self = VehicleArticulatedLinear()
            self.mF0 = 5200;
            self.mR0 = 2400;
            self.mF = 6000;
            self.mR = 10000;
            self.mM = 17000;
            self.IT = 46000;
            self.IS = 450000;
            self.lT = 3.5;
            self.lS = 7.7;
            self.c = -0.3;
            self.nF = 2;
            self.nR = 4;
            self.nM = 8;
            self.wT = 2.6;
            self.wS = 2.4;
            self.muy = 0.3;
            self.deltaf = 0;
            self.Fxf = 0;
            self.Fxr = 0;
            self.Fxm = 0;
        end

        function dx = Model(self,t,states,tspan)
            % Vehicle parameters
            mT = self.mT;
            mS = self.mS;
            % IT = self.IT;
            % IS = self.IS;
            a = self.a;
            b = self.b;
            c = self.c;
            d = self.d;
            e = self.e;
            nF = self.nF;
            nR = self.nR;
            nM = self.nM;
            g = 9.81;

            % Vertical forces
            FzF = self.mF * g;
            FzR = self.mR * g;
            FzM = self.mM * g;
            muy = self.muy;

            v0 = 20;

            % State
            X = states(1,1);
            Y = states(2,1);
            PSI     = states(3,1);
            PHI     = states(4,1);
            V       = states(5,1);
            ALPHAT  = states(6,1);
            dPSI    = states(7,1);
            dPHI    = states(8,1);

            if isa(self.deltaf,'function_handle')
                deltaf = self.deltaf([X;Y;PSI;PHI;V;ALPHAT;dPSI;dPHI],t);
            elseif length(self.deltaf)>1
                deltaf = interp1(tspan,self.deltaf,t);
            else
                deltaf = self.deltaf;
            end

            % Slip angles - linear
            ALPHAF = ALPHAT + a/v0*dPSI - deltaf;
            ALPHAR = ALPHAT - b/v0*dPSI;
            ALPHAM = ALPHAT + PHI - (dPSI*(b + c + d + e))/v0 + (dPHI*(d + e))/v0;

            % Longitudinal forces
            if isa(self.Fxf,'function_handle')
                FxF = self.Fxf([X;Y;PSI;PHI;V;ALPHAT;dPSI;dPHI],t);
            elseif length(self.Fxf)>1
                FxF = interp1(tspan,self.Fxf,t);
            else
                FxF = self.Fxf;
            end

            if isa(self.Fxr,'function_handle')
                FxR = self.Fxr([X;Y;PSI;PHI;V;ALPHAT;dPSI;dPHI],t);
            elseif length(self.Fxr)>1
                FxR = interp1(tspan,self.Fxr,t);
            else
                FxR = self.Fxr;
            end

            if isa(self.Fxm,'function_handle')
                FxM = self.Fxm([X;Y;PSI;PHI;V;ALPHAT;dPSI;dPHI],t);
            elseif length(self.Fxm)>1
                FxM = interp1(tspan,self.Fxm,t);
            else
                FxM = self.Fxm;
            end

            % Lateral forces - Characteristic curve
            FyF = nF*self.tire.Characteristic(ALPHAF,FzF/nF,muy);
            FyR = nR*self.tire.Characteristic(ALPHAR,FzR/nR,muy);
            FyM = nM*self.tire.Characteristic(ALPHAM,FzM/nM,muy);

            A = [ 0  0   0  0  1   0                       0  0;...
                  0  0  v0  0  0  v0                       0  0;...
                  0  0   0  0  0   0                       1  0;...
                  0  0   0  0  0   0                       0  1;...
                  0  0   0  0  0   0                       0  0;...
                  0  0   0  0  0   0           (-v0*(mS + mT))  0;...
                  0  0   0  0  0   0  (mS*(d*v0 + v0*(b + c)))  0;...
                  0  0   0  0  0   0                (-d*mS*v0)  0];

            B = [ 0  0  0  0  0   0                0;...
                  0  0  0  0  0   0                0;...
                  0  0  0  0  0   0                0;...
                  0  0  0  0  0   0                0;...
                  0  1  1  1  0   0                0;...
                  0  0  0  0  1   1                1;...
                  0  0  0  0  a  -b  (- b - c - d - e);...
                  0  0  0  0  0   0            (d + e)];

            vetEst = [X ; Y ; PSI ; PHI ; V ; ALPHAT ; dPSI ; dPHI];
            vetEnt = [deltaf ; FxF ; FxR ; FxM ; FyF ; FyR ; FyM];

            % Integrator output
            dx = A*vetEst + B*vetEnt;
        end

        function [value,isterminal,direction] = velocity(~,~,states)
            % If the velocity is less than 0.1m/s the integrator stops.
            % The MassMatrix is singular when the velocity is 0 m/s.
            value = states(5,1) - 0.1;
            isterminal = 1;
            direction = -1;
        end

        function E = MassMatrix(self,~,~)
            % Vehicle parameters
            mT = self.mT;
            mS = self.mS;
            IT = self.IT;
            IS = self.IS;
            % a = self.a;
            b = self.b;
            c = self.c;
            d = self.d;
            % e = self.e;
            % deltaf = self.deltaf;
            % nF = self.nF;
            % nR = self.nR;
            % nM = self.nM;

            % g = 9.81;

            v0 = 20;

            % Mass matrix
            E = [ 1  0  0  0        0                        0                                               0                            0;...
                  0  1  0  0        0                        0                                               0                            0;...
                  0  0  1  0        0                        0                                               0                            0;...
                  0  0  0  1        0                        0                                               0                            0;...
                  0  0  0  0  (mS + mT)                      0                                               0                            0;...
                  0  0  0  0        0             (v0*(mS + mT))                                 (-mS*(b + c + d))                         (d*mS);...
                  0  0  0  0        0  ( -mS*v0*(b + c + d))            (IS + IT + mS*(b + c + d)^2)  (- IS - mS*(d^2 + (b + c)*d));...
                  0  0  0  0        0                  (d*mS*v0)                     (- IS - mS*(d^2 + (b + c)*d))                  (mS*d^2 + IS)];

        end
    end
end

%% See Also
%
% <../../index.html Home>
%

author: andresmendes

repo: Vehicle-Dynamics-Lateral


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

% This is a dummy function

pvec    = [];
pstruct = struct;

return;

author: translationalneuromodeling

repo: tapas


text:
function varargout = squeeze(varargin)
%SQUEEZE   Squeeze a SPHEREFUN to one variable, if possible.
%   G = squeeze(F) returns a SPHEREFUN if F depends on LAMBDA and THETA in
%   spherical coordinates. If F depends only on the lambda-variable a row
%   CHEBFUN is returned and if it depends on just the theta-variable a
%   column CHEBFUN is returned.

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

[varargout{1:nargout}] = squeeze@separableApprox(varargin{:});
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% (Internal) Intersection of two sets with tolerance
%
%     [sft_intersect1 idx1 idx2 ] = soft_set_intersect( val1, val2, win_size )
%
%
% Arguments:
%
%   + val1, val2: data elements to be intersected,
%
%   + win_size: tolerance to consider val1(i) == val2(i). In fact we
%   consider equal two elements if val1(val1 >= (val2(i) - win_size) & val1 <= (val2(i) + win_size) )
%
% Output:
%
%   + sft_intersect1: the elements in the soft intersection
%
%   + idx1, idx2: the indexes of val1(idx1) which are in the soft
%   intersection.
%
% Example:
%
%
% See also soft_set_difference, soft_set_union
%
% Author: Mariano Llamedo Soria (llamedom at frba.utn.edu.ar)
% Version: 0.1 beta
% Birthdate: 17/12/2010
% Last update: 17/12/2010
% Copyright 2008-2015
%
function [sft_intersect1, idx1 ] = soft_set_intersect( val1, val2, win_size )

% Find the intersection between index sequence val1-2 within a windows
% win_size.
idx1 = [];

sft_intersect1 = unique(cell2mat(arrayfun(@(a)( colvec( val1(val1 >= (a - win_size) & val1 <= (a + win_size) ) ) ), colvec(val2), 'UniformOutput', false)));

if( nargout > 1)
    [ ~, idx1] = intersect(val1, sft_intersect1);
end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [mask] = estimate_mask(s, filter_params, iterations)
% Creates binary mask based on SNR estimate of the signal
%   [mask] = estimate_mask(s, filter_params, iterations)
%     s              Source signal
%     filter_params  Filter parameters for smoothing variance estimate
%     iterations     Number of iterations for estimating noise variance

% Copyright (C) 2004, 2005 DSS MATLAB package team (dss@cis.hut.fi).
% Distributed by Laboratory of Computer and Information Science,
% Helsinki University of Technology. http://www.cis.hut.fi/projects/dss/.
% $Id$

filter_h = @denoise_filter;

if nargin<2;
    % default filtering is lowpass DCT
    T = length(s);
    t = 1:T;
    filter_params.filter_dct = exp(-0.5*(t-1).^2 / (T/64).^2 );
end
if nargin<3; iterations = 8; end

% noise variance estimate for gaussian noise
[p, var_smooth] = feval(filter_h, filter_params, randn(1, length(s)).^2, []);
noise_var = est_noise_var(var_smooth, 0);
normalization_c = 1 / noise_var;

% smoothed variance
[p, var_totsm] = feval(filter_h, filter_params, s.^2, []);

% iterate signal noise variance estimate
var_noise = 1;
for i = 1 : iterations
  var_noise = est_noise_var(var_totsm, var_noise)*normalization_c;
end

% create binary mask
mask = var_totsm>var_noise;

%DEBUG
%fprintf('Noise: %d\n', var_noise);
%clf
%subplot(3, 1, 1);
%plot(s);
%subplot(3, 1, 2);
%plot(var_totsm);
%subplot(3, 1, 3);
%plot(mask);
%axis([0 length(s) -0.5 1.5])

% --------
function noise_var = est_noise_var(var_tot, var_noise)
% Estimates noise variance based on total variance estimate and
% previous noise variance estimate.
noise_var = exp(mean(log(var_tot+repmat(var_noise,1,size(var_tot,2))),2))-var_noise;


text:
% Delete cache that had been created for input file
%
function deleteCache(inputFile)
    global gBntData;
    try
        trials = io.parseGeneralFile(inputFile);
    catch
        io.readLinearInputFile(inputFile);
        trials = gBntData;
    end
    for i = 1:length(trials)
        trial = trials{i};

        [baseFolder, firstName] = helpers.fileparts(trial.sessions{1});
        if isempty(baseFolder)
            baseFolder = trial.path;
        end

        trial = io.detectRecordingSystem(trial);
        if strcmpi(trial.system, bntConstants.RecSystem.Neuralynx)
            baseFolder = trial.sessions{1};
        end

        if ~isfield(trial, 'basename') || isempty(trial.basename)
            trial.basename = firstName;
            trial.path = baseFolder;

            if length(trial.sessions) > 1
                if length(trial.sessions) > size(trial.cuts, 2)
                    [~, lastName] = helpers.fileparts(trial.sessions{end});
                    trial.basename = strcat(firstName, '+', lastName(end-1:end));
                end
            end
        end

        %% do the inheritance check
        % units
        if isempty(trial.units)
            if trial.inheritUnits
                if i > 1
                    trial.units = trials{i-1}.units;
                    trials{i}.units = trials{i-1}.units;
                end
            else
                if strcmpi(trial.system, bntConstants.RecSystem.Neuralynx)
                    trial = io.neuralynx.detectTetrodes(trial);
                else
                    trial = io.axona.detectTetrodes(trial);
                end
            end
        end
        if isempty(trial.units)
            error('BNT:noUnits', 'There are no units assigned with trial %u (first session ''%s'').\nCheck your input file!', i, trials{i}.sessions{1});
        end

        % shape info
        if ~isfield(trial.extraInfo, 'shape')
            if i > 1 && isfield(trials{i-1}.extraInfo, 'shape')
                trial.extraInfo.shape = trials{i-1}.extraInfo.shape;
                trials{i}.extraInfo.shape = trials{i-1}.extraInfo.shape;
            end
        end

        trial = io.initCuts(trial);

        tetrodes = unique(trial.units(:, 1), 'stable');
        numTetrodes = length(tetrodes);
        for t = 1:numTetrodes
            if t <= size(trial.cuts, 1)
                cutNames = trial.cuts(t, :);
            else
                if ~isempty(trial.cuts)
                    cutNames = trial.cuts(1, :);
                else
                    cutNames = {[]};
                end
            end
            notEmptyInd = ~cellfun(@isempty, cutNames);

            if ~isempty(trial.cuts) && all(notEmptyInd)
                unitCut = cutNames;
            else
                unitCut = [];
            end
            key = helpers.uniqueKey(trial.basename, unitCut);
            files = fullfile(trial.path, sprintf('*_T%u*_%s.mat', tetrodes(t), key));
            delete(files);
        end

        key = helpers.uniqueKey(trial.basename, trial.cuts);
        files = fullfile(trial.path, sprintf('*_%s.mat', key));
        delete(files);

        posFile = helpers.uniqueFilename(trial, 'posCleanScale');
        if exist(posFile, 'file')
            delete(posFile);
        end

        posFile = helpers.uniqueFilename(trial, 'posClean');
        if exist(posFile, 'file')
            delete(posFile);
        end

        posFile = helpers.uniqueFilename(trial, 'pos');
        if exist(posFile, 'file')
            delete(posFile);
        end
        % here is a small hack because the knowledge of kalmanFile is used
        [filePath, fileName] = helpers.fileparts(posFile);
        id = strfind(fileName, '_');
        fileName(id:end) = [];
        kalmanFile = fullfile(filePath, sprintf('%s_kalman*.mat', fileName));
        if exist(kalmanFile, 'file')
            delete(kalmanFile);
        end
    end
end

author: kavli-ntnu


text:
function fx = p57_fun ( n, x )

%*****************************************************************************80
%
%% P57_FUN evaluates the integrand for problem 57.
%
%  Interval:
%
%    0 <= x <= 1
%
%  Integrand:
%
%    x^(3/2)
%
%  Antiderivative:
%
%    (2/5) * x^(5/2)
%
%  Exact Integral:
%
%    0.4
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 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:n) = sqrt ( x(1:n).^3 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ax, data] = swissRollScatter(Y, shade);

% SWISSROLLSCATTER 3-D scatter plot with colors.
% FORMAT
% DESC produces a 3-D scatter plot with colors of the type peple use for  the 'swiss roll data'.
% ARG Y : scatter points.
% ARG shade : color indicator for each data point so that they may be given different
% colours.
% RETURN ax : the axes handle where the scatter plot was placed.
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006, 2008, 2011
%
% SEEALSO : fgplvmVisualise, lvmTwoDPlot, lvmScatterPlot

% MLTOOLS

  shade = shade - min(shade)+eps;
  shade = shade/max(shade);
  shade = ceil(shade*64);

  ax = gca;
  jt = colormap('jet');
  plot3(Y(:, 1), Y(:, 2), Y(:, 3), '.');

  xLim = get(ax, 'xlim');
  yLim = get(ax, 'ylim');
  zLim = get(ax, 'zlim');
  %cla
  set(ax, 'xLim', xLim);
  set(ax, 'yLim', yLim);
  set(ax, 'zLim', zLim);
  hold on
  for i = 1:size(Y, 1)
    data(i) = plot3(Y(i, 1), Y(i, 2), Y(i, 3), '.');
    set(data(i), 'color', jt(shade(i), :), 'markersize', 10);
  end

  set(ax, 'fontname', 'arial');
  set(ax, 'fontsize', 20);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat