text:
function determ = circulant2_determinant ( n )

%*****************************************************************************80
%
%% CIRCULANT2_DETERMINANT returns the determinant of the CIRCULANT2 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real DETERM, the determinant.
%
  w = c8vec_unity ( n );

  lambda = zeros ( n, 1 );
  lambda(1:n) = n;
  for i = n-1 : -1 : 1
    lambda(1:n) = lambda(1:n) .* w(1:n) + i;
  end
%
%  First eigenvalue is "special".
%
  determ = real ( lambda(1) );
%
%  Eigenvalues 2, 3     through ( N + 1 ) / 2 are paired with complex conjugates.
%
  for i = 2 : floor ( ( n + 1 ) / 2 )
    determ = determ * ( abs ( lambda(i) ) )^2;
  end
%
%  If N is even, there is another unpaired eigenvalue.
%
  if ( mod ( n, 2 ) == 0 )
    determ = determ * real ( lambda((n/2)+1) );
  end

  return
end


text:
function out = floor(in)

out = cellfun(@floor, in, 'uniformoutput', 0);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Sest = cosamp_cgls(Phi,u,K,tol,maxiterations, err_tol)

% Cosamp algorithm
%   Input
%       K : sparsity of Sest
%       Phi : measurement matrix
%       u: measured vector
%       tol : tolerance for approximation between successive solutions.
%   Output
%       Sest: Solution found by the algorithm
%
% Algorithm as described in "CoSaMP: Iterative signal recovery from
% incomplete and inaccurate samples" by Deanna Needell and Joel Tropp.
%

% This implementation was written by David Mary,
% but modified 20110707 by Bob L. Sturm to make it much clearer,
% and corrected multiple times again and again.
% To begin with, see: http://media.aau.dk/null_space_pursuits/2011/07/ ...
% algorithm-power-hour-compressive-sampling-matching-pursuit-cosamp.html
%
% This script/program is released under the Commons Creative Licence
% with Attribution Non-commercial Share Alike (by-nc-sa)
% http://creativecommons.org/licenses/by-nc-sa/3.0/
% Short Disclaimer: this script is for educational purpose only.
% Longer Disclaimer see  http://igorcarron.googlepages.com/disclaimer

%%Code modified by Praneeth 14March2017 to include case when there is no
%%signal to estimate;replaced pinv() by A\b syntax for speed.
% Initialization
Sest = zeros(size(Phi,2),1);
v = u;
t = 1;
numericalprecision = 1e-12;
T = [];
while ((t <= maxiterations) && (norm(v)/norm(u) > tol) ...
        && (norm(Phi'*v) > err_tol))

  y = abs(Phi'*v);
  [vals,z] = sort(y,'descend');
  Omega = find(y >= vals(2*K) & y > numericalprecision);

  T = union(Omega,T);
  x = v(T);
  %b = pinv(Phi(:,T))*u;
%   b = Phi(:,T) \ u;     %% edit PKN -- 3/3/17
  b = cgls(Phi(:, T), u, 0, err_tol, 3, 0, x);
  [vals,z] = sort(abs(b),'descend');
  Kgoodindices = (abs(b) >= vals(K) & abs(b) > numericalprecision);
  T = T(Kgoodindices);
  Sest = zeros(size(Phi,2),1);
  b = b(Kgoodindices);
  Sest(T) = b;
  v = u - Phi(:,T)*b;
  t = t+1;
end

author: andrewssobral

repo: lrslibrary


text:
classdef prtRvMvn < prtRv & prtRvMemebershipModel
    % prtRvMvn  Multivariate normal random variable
    %
    %   RV = prtRvMvn creates a prtRvMvn object with empty mean and
    %   covariance matrices. The mean and covariance matrices must be set
    %   either directly, or by calling the MLE method.
    %
    %   RV = prtRvMvn('covarianceStructure', VALUE) enforces a covariance
    %   structure, which may be either 'full', 'spherical', or 'diagonal'.
    %   Setting this property to 'spherical' or 'diagonal' will enforce
    %   this structure onto the existing covariance matrix, or one
    %   estimated by calling the MLE method.
    %
    %   RV = prtRvMvn(PROPERTY1, VALUE1,...) creates a prtRvMv object RV
    %   with properties as specified by PROPERTY/VALUE pairs.
    %
    %   A prtRvMvn object inherits all properties from the prtRv class. In
    %   addition, it has the following properties:
    %
    %   covarianceStructure - A string specifying the structure of the
    %                         covariance matrix to estimate or enforce.
    %                         Valid values are 'full','spherical', or
    %                         'diagonal'
    %   mu                  - The mean of the distribution, which is
    %                         a 1 x nDimensions vector.
    %   sigma               - The covariance matrix of the distribution,
    %                         which is a nDimensions x nDimensions
    %                         matrix.
    %
    %  A prtRvMvn object inherits all methods from the prtRv class. The MLE
    %  method can be used to estimate the distribution parameters from
    %  data.
    %
    %  Example:
    %
    %  dataSet    = prtDataGenUnimodal;   % Load a dataset consisting of 2
    %                                     % classes
    %  % Extract one of the classes from the dataSet
    %  dataSetOneClass = prtDataSetClass(dataSet.getObservationsByClass(1));
    %
    %  RV = prtRvMvn;                       % Create a prtRvMvn object
    %  RV = RV.mle(dataSetOneClass.getX);   % Compute the maximum
    %                                       % likelihood estimate from the
    %                                       % data
    %  RV.plotPdf                           % Plot the pdf
    %
    %  RVspec = prtRvMvn;                   % Create another prtRvMvn
    %                                       % object
    %  RVspec.mu = [1 2];                   % Specify the mean
    %  RVspec.sigma = [2 -1; -1 2]          % Specify the covariance
    %  figure;
    %  RVspec.plotPdf                       % Plot the pdf
    %  sample = RVspec.draw(1)              % Draw 1 random sample from the
    %                                       % Distribution
    %
    %   See also: prtRv, prtRvGmm, prtRvMultinomial, prtRvUniform,
    %   prtRvUniformImproper, prtRvVq, prtRvDiscrete

    properties (SetAccess = 'private')
        name = 'Multi-Variate Normal';
        nameAbbreviation = 'RVMVN';
    end

    properties (SetAccess = 'protected')
        isSupervised = false;
        isCrossValidateValid = true;
    end

    properties (Dependent)
        covarianceStructure  % The covariance structure
        mu                   % The mean vector
        sigma                % The covariance matrix
    end
    properties (SetAccess = 'private', GetAccess = 'private', Hidden = true)
        muDepHelper
        sigmaDepHelper
        covarianceStructureDepHelper = 'full';
    end

    properties (Hidden = true)
        covarianceBias = [];
    end

    properties (Hidden = true, Dependent = true)
        nDimensions
    end

    properties (SetAccess = 'private', Hidden = true)
        trueCovariance
        badCovariance = false;
    end

    methods
        function R = prtRvMvn(varargin)
            R = constructorInputParse(R,varargin{:});
        end

        function R = mle(R,X)
           % MLE Compute the maximum likelihood estimate
            %
            % RV = RV.mle(X) computes the maximum likelihood estimate based
            % the data X. X should be nObservations x nDimensions.

            X = R.dataInputParse(X); % Basic error checking etc

            R.mu = mean(X,1);
            if size(X,1) == 1
                % A single observation
                % This is bad..
                % You can't call cov for this case
                % Since we always want to output sum( (x_i - u)'*(x_i - u))
                % We will output a matrix of zeros of the correct size.
                % Error checking for a propert covariance will happen later
                % after the bias has been applied.
                R.sigma = zeros(size(X,2));
            else
                R.sigma = cov(X);
            end

        end

        function vals = pdf(R,X)
            % PDF Output the pdf of the random variable evaluated at the points specified
            %
            % pdf = RV.pdf(X) returns  the pdf of the prtRv
            % object evaluated at X. X must be an N x nDims matrix, where
            % N is the number of locations to evaluate the pdf, and nDims
            % is the same as the number of dimensions, nDimensions, of the
            % prtRv object RV.

            [isValid, reasonStr] = R.isValid;
            assert(isValid,'PDF cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            assert(~R.badCovariance,'Covariance matrix is not positive definite. Consider modifying "covarianceStructure".');
            X = R.dataInputParse(X); % Basic error checking etc

            assert(size(X,2) == R.nDimensions,'Data, RV dimensionality missmatch. Input data, X, has dimensionality %d and this RV has dimensionality %d.', size(X,2), R.nDimensions)
            assert(isnumeric(X) && ndims(X)==2,'X must be a 2D numeric array.');

            vals = exp(prtRvUtilMvnLogPdf(X,R.mu,R.sigma));
        end

        function vals = logPdf(R,X)
            % LOGPDF Output the log pdf of the random variable evaluated at the points specified
            %
            % logpdf = RV.logpdf(X) returns the logarithm of value of the
            % pdf of the prtRv object evaluated at X. X must be an N x
            % nDims matrix, where N is the number of locations to evaluate
            % the pdf, and nDims is the same as the number of dimensions,
            % nDimensions, of the prtRv object RV.

            X = R.dataInputParse(X); % Basic error checking etc
            try
               vals = prtRvUtilMvnLogPdf(X,R.mu,R.sigma);
            catch ME
                [isValid, reasonStr] = R.isValid;
                assert(isValid,'LOGPDF cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);
                assert(size(X,2) == R.nDimensions,'Data, RV dimensionality missmatch. Input data, X, has dimensionality %d and this RV has dimensionality %d.', size(X,2), R.nDimensions)
                throw(ME);
            end
        end

        function varargout = plotCdf(R,varargin)
            % PLOTCDF Plots the CDF of the prtRv
            assert(R.nDimensions == 1,'prtRvMvn.plotCdf can only be used for 1D RV objects.');

            varargout = cell(nargout,1);

            [varargout{:}] = plotCdf@prtRv(R,varargin{:});

        end

        function vals = cdf(R,X)
            % CDF Output the cdf of the random variable evaluated at the points specified
            %
            % cdf = RV.cdf(X) returns the value of the cdf of the prtRv
            % object evaluated at X. X must be an N x nDims matrix, where
            % N is the number of locations to evaluate the pdf, and nDims
            % is the same as the number of dimensions, nDimensions, of the
            % prtRv object RV.

            [isValid, reasonStr] = R.isValid;
            assert(isValid,'CDF cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            assert(~R.badCovariance,'Covariance matrix is not positive definite. Consider modifying "covarianceStructure"');
            X = R.dataInputParse(X); % Basic error checking etc

            assert(size(X,2) == R.nDimensions,'Data, RV dimensionality missmatch. Input data, X, has dimensionality %d and this RV has dimensionality %d.', size(X,2), R.nDimensions)
            vals = prtRvUtilMvnCdf(X,R.mu,R.sigma);
        end

        function vals = draw(R,N)
            % DRAW  Draw random samples from the distribution described by the prtRv object
            %
            % VAL = RV.draw(N) generates N random samples drawn from the
            % distribution described by the prtRv object RV. VAL will be a
            % N x nDimensions vector, where nDimensions is the number of
            % dimensions of RV.

            assert(numel(N)==1 && N==floor(N) && N > 0,'N must be a positive integer scalar.')
            [isValid, reasonStr] = R.isValid;
            assert(isValid,'DRAW cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            assert(~R.badCovariance,'Covariance matrix is not positive definite. Consider modifying "covarianceStructure"');
            vals = prtRvUtilMvnDraw(R.mu,R.sigma,N);
        end
    end

    methods (Hidden=true)
        function [val, reasonStr] = isValid(R)

            if numel(R) > 1
                val = false(size(R));
                for iR = 1:numel(R)
                    [val(iR), reasonStr] = isValid(R(iR));
                end
                return
            end

            val = ~isempty(R.sigma) && ~isempty(R.mu);

            if val
                reasonStr = '';
            else
                badCov = isempty(R.sigma);
                badMean = isempty(R.mu);

                if badCov && ~badMean
                    reasonStr = 'because sigma has not been set';
                elseif ~badCov && badMean
                    reasonStr = 'because mu has not been set';
                elseif badCov && badMean
                    reasonStr = 'because mu and sigma have not been set';
                else
                    reasonStr = 'because of an unknown reason';
                end
            end
        end

        function val = plotLimits(R)
            [isValid, reasonStr] = R.isValid;
            if isValid
                minX = min(R.mu, [], 1)' - 2*sqrt(diag(R.sigma));
                maxX = max(R.mu, [], 1)' + 2*sqrt(diag(R.sigma));

                val = zeros(1,2*R.nDimensions);
                val(1:2:R.nDimensions*2-1) = minX;
                val(2:2:R.nDimensions*2) = maxX;
            else
                error('prtRvMvn:plotLimits','Plotting limits can not be determined for this RV. It is not yet valid %s',reasonStr)
            end
        end

        function R = weightedMle(R,X,weights)
            assert(numel(weights)==size(X,1),'The number of weights must mach the number of observations.');

            weights = weights(:);

            Nbar = sum(weights);
            R.mu = 1/Nbar*sum(bsxfun(@times,X,weights),1);
            X = bsxfun(@times,bsxfun(@minus,X,R.mu),sqrt(weights));
            R.sigma = 1/Nbar*(X'*X);
        end
    end

    % Get methods
    methods
        function val = get.nDimensions(R)
            val = getNumDimensions(R);
        end

        function val = get.mu(R)
            val = R.muDepHelper;
        end
        function val = get.sigma(R)
            val = R.sigmaDepHelper;
        end
        function val = get.covarianceStructure(R)
            val = R.covarianceStructureDepHelper;
        end
    end
    methods (Access = 'protected')
        function val = getNumDimensions(R)
            if ~isempty(R.mu)
                val = length(R.mu);
            elseif ~isempty(R.sigma)
                val = size(R.sigma,2);
            else
                val = [];
            end
        end
    end
    % Set Methods
    methods
        function R = set.covarianceStructure(R,covarianceStructure)
            % Find and fix known abbreviations

            assert(ischar(covarianceStructure),'covarianceStructure must be a string that is either, full, diagonal, or spherical');

            covarianceStructure = lower(covarianceStructure);

            if strcmpi(covarianceStructure,'diag')
                covarianceStructure = 'diagonal';
            end

            % Limit the options for the covariance structure
            if ~(strcmpi(covarianceStructure,'full') || ...
                    strcmpi(covarianceStructure,'diagonal') || ...
                    strcmpi(covarianceStructure,'spherical'))
                error('%s is not a valid covariance structure. Possible types are, full, diagonal, and spherical',covarianceStructure);
            end
            R.covarianceStructureDepHelper = covarianceStructure;

            % Redo the covariance to reflect the updated covarianceStructure
            if ~isempty(R.sigma)
                R.sigma = R.trueCovariance; % This will call set.sigma()
            end
        end

        function R = set.mu(R,meanVal)
            if ~isempty(R.sigma) && size(meanVal,2) ~= size(R.sigma,2)
                error('prtRvMvn:dimensions','Dimensions mismatch between supplied mu and prtRvMvn dimensionality');
            end
            R.muDepHelper = meanVal;
        end

        function R = set.sigma(R,sigma)
            if size(sigma,1) ~= size(sigma,2)
                error('Covariance matrix must be square.')
            end

            if ~isempty(R.mu) && size(sigma,1) ~= R.nDimensions
                error('prtRvMvn:dimensions','Dimensions mismatch between sigma and prtRvMvn dimensionality')
            end

            % Save this input as a true hidden sigma
            R.trueCovariance = sigma;

            % Enforce the covariance structure
            switch lower(R.covarianceStructure)
                case 'full'
                    R.sigmaDepHelper = sigma;
                case 'diagonal'
                    R.sigmaDepHelper = eye(size(sigma)).*sigma;
                case 'spherical'
                    R.sigmaDepHelper = eye(size(sigma))*mean(diag(sigma));
            end

            if ~isempty(R.covarianceBias)
                R.sigmaDepHelper = R.sigmaDepHelper + R.covarianceBias.*eye(size(sigma));
            end

            [dontNeed, cholErr] = chol(R.sigmaDepHelper); %#ok<ASGLU>
            if cholErr ~=0
                R.badCovariance = true;
                warning('prt:prtRvMvn','Covariance matrix is not positive definite. This may cause errors. Consider modifying "covarianceStructure".');
            end
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
%% This file is the ODE file used to simulate the Michaelis-Menten kinetics.
% Date: 04/24/2019
% Coded By: K.Kahirman

% Here Vmax is rhe maximum rate of reaction time
% Km is the concentration of half-maximal reaction rate
%

%% ODE functions
function dx=MMK_ODE(t,x,jx,Vmax,Km)
dx=jx-(Vmax*x(1,:)./(Km+x(1,:)));

author: dynamicslab

repo: SINDy-PI


text:
clc;
clearvars;
close all;
M = 512;
N = 1024;
K = 5;
S = 1000;
L = 2;
Phi = spx.dict.simple.gaussian_mtx(M, N);
G = Phi' * Phi;
X = spx.data.synthetic.SparseSignalGenerator(N, K, S).biGaussian;
Y = Phi * X;
options.verbose = 1;
tstart = tic;
X1 = spx.fast.omp(Phi, Y, K, 1e-3, options);
t1 = toc(tstart);
options.verbose = 1;
options.ls_method = 'chol';
tstart = tic;
X2 =  spx.fast.gomp(Phi, Y, K, L, 1e-3, options);
t2 = toc(tstart);
cmpare = spx.commons.SparseSignalsComparison(X1, X2, K);
cmpare.summarize();

failures = sum(cmpare.support_similarity_ratios() ~= 1);
fprintf('Failed detections %d (%.2f %%)\n', failures, failures* 100 / S);
gain_x = t1/t2;
if t1 > t2
    gain_perc = (t1 - t2) * 100/ t1;
else
    gain_perc = (t1 - t2) * 100/ t2;
end
fprintf('Time: OMP, %.2f s, GOMP: %.2f s, gain: %.2f (%.1f %%)\n', t1, t2, gain_x, gain_perc);

author: indigits

repo: sparse-plex


text:
%  Figure 10.26      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
%  fig10_26.m is a script to generate Fig. 10.26,
%  root locus of the colocated
%  design for the satellite with PD compensation

% parameter values
m=[1, .1]; k0=[0, .091] ; d0=[0, .0036]; k1=[0, .4];

% call function
[f,g,h,j]=twomass(m,k0,d0); [f1,g,h,j] = twomass(m,k1,d0);
h1=[0, 0, 1, 0];

nc1=0.25*[2, 1];
dc1=[1/40, 1];

[ac,bc,cc,dc]=tf2ss(nc1, dc1);
sysc=ss(ac, bc,cc,dc);
sysp1=ss(f,g,h1,j);
sysol= series(sysc,sysp1);
[aol,bol,col,dol]=ssdata(sysol);
[acl]=aol-bol*col;
ccl2=[0*cc h];
sysp11=ss(f1,g,h1,j);
sysol1= series(sysc,sysp11);
[aol1,bol1,col1,dol1] =ssdata(sysol1);
acl1=aol1-bol1*col1;
hold off ; clf
rlocus(aol,bol,col,dol);
v =[-3,2,-1.5,1.5];
axis(v);
title('Fig. 10.26 Root locus for D_5(s)G_{co}(s).')
grid

author: Sable

repo: mcbench-benchmarks


text:
cleardataquestion   % clearing the data
savethis = 1;       % saves the results
params              % parameters for the computitation

delta_vec = [.1 .3 .5 .8 1]; %all...
ncomp_vec = [2];
for nc_ix = 1:length(ncomp_vec)
    peval.ncomp = ncomp_vec(nc_ix); %number of components without the background component....
    for delta = delta_vec
        % data location for each computaution:
        peval.data_path = [getenv('HOME') '/project/data/qdots/D2ptBMrandOblique' '/'];
        peval.data_dir= ['D_N2_delta' num2str(delta*100) '_NoNoise'];
        peval.data_name=peval.data_dir;

        % results name:
        peval.res_nameappendix = ['_nc' num2str(peval.ncomp)];
        peval.res_dir = [peval.data_dir peval.res_nameappendix];
        peval.res_name = peval.res_dir;

        %reads data:
        [dpixc, dpixc_ind, blinkmat, peval, p] = readSimulSimple(peval);
        peval.bg = peval.offset; %cheating here

        for iteration=1:10 % several evaluations of the same data
            for nn = 1 : peval.nrestartH
                % initialization:
                winit_pix = init_w_meandata1res;
                hinit = rand(peval.ncomp, size(dpixc,3));
                %nmf separation:
                [res, peval] = separcomp(dpixc, peval, winit_pix, hinit);

                %saving results
                if savethis
                    saveresultscompleteSimple
                end
            end

            % MAP fiting with nmf initialization
            peval.filenamebase = ['iter' num2str(iteration)];
            % peval.sigmaPSF = 1.3; %sigma of the PSF gauss approx (in pixels)
            if savethis; cd (peval.res_dir); end
            [ihs, x_mu, y_mu, sig] = plotreswh2(res, peval, dpixc, p, savethis, 1, 0, 1, 1,0);
            if savethis; cd ..; end
            % initialization form nmf ressults:
            x1in=reshape(res.h(1:peval.ncomp,:), 1,peval.ncomp*peval.nt);
            x2in=-1+[x_mu, y_mu];
            [x1in,x2in]=initx1x2(Hinit,Winit,x1,x2,x1rand,x2rand);
            Wres=x2in;
            Hres=x1in;
            pointlogWall=[]; flog=[]; pointlogHall=[];
            for kk=1:peval.nIterAlter %alternating optimization

                % W fixed, optimizing H (intensities)
                [Hres, optionsH, flogH, pointlogH] = conjgrad(peval.Hupdate, Hres, optionsH, ['grad' peval.Hupdate],reshape(dpixc, peval.nx*peval.ny, peval.nt), peval.sigmaPSF, peval.alpha, peval.beta, peval, Wres(1:peval.ncomp), Wres(peval.ncomp+1:end));
                flog=[flog; flogH];
                pointlogHall=[pointlogHall; pointlogH];
                % H fixed, optimizing W (positions)
                [Wres, optionsW, flogW, pointlogW] = conjgrad(peval.Wupdate, Wres, optionsW, ['grad' peval.Wupdate],reshape(dpixc, peval.nx*peval.ny, peval.nt), peval.sigmaPSF, peval.alpha, peval.beta, peval, Hres);
                flog=[flog; flogW];
                pointlogWall=[pointlogWall; pointlogW];

                peval.Wres=Wres;
                %% plotting & saves results
                if savethis; cd (peval.res_dir); end
                showimage=0;
                res.dpixc=dpixc;
                plotloglikGaP;
                if savethis; cd ..; end

            end
        end
    end
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function tpl = set(tpl,action,varargin)
%TEMPLATE/SET Edit data stored in a Template object
%  TPL = SET(TPL,ACTION,VARARGIN)
%     ACTION 'root'
%     ACTION 'unknowns'
%     ACTION 'file'
%     ACTION 'block'
%     ACTION 'var'

%  Copyright (C) 2003 Guillaume Flandin <Guillaume@artefact.tk>
%  $Revision: 1.0 $Date: 2003/05/05 22:19:51 $

error(nargchk(3,5,nargin));

switch lower(action)
	case 'root'
		error(nargchk(3,3,nargin));
		if exist(varargin{1},'dir')
			tpl.root = varargin{1};
		else
			error('[Template] No such directory.');
		end
	case 'unknowns'
		error(nargchk(3,3,nargin));
		if ismember(varargin{1},{'remove' 'comment' 'keep'})
			tpl.unknowns = varargin{1};
		else
			error('[Template] Unknowns: ''remove'', ''comment'' or ''keep''.');
		end
	case 'file'
		error(nargchk(4,4,nargin));
		if iscellstr(varargin{1})
			for i=1:length(varargin{1})
				ind = find(ismember(tpl.handles,varargin{1}{i}));
				if isempty(ind)
					tpl.handles{end+1} = varargin{1}{i};
					if strcmp(varargin{2}{i}(1),filesep) %- absolute path (Unix)
						tpl.file{end+1} = varargin{2}{i};
					else %- relative path
						tpl.file{end+1} = fullfile(tpl.root,varargin{2}{i});
					end
				else
					if strcmp(varargin{2}{i}(1),filesep) %- absolute path (Unix)
						tpl.file{ind} = varargin{2}{i};
					else %- relative path
						tpl.file{ind} = fullfile(tpl.root,varargin{2}{i});
					end
				end
			end
		elseif ischar(varargin{1})
			tpl = set(tpl,'file',cellstr(varargin{1}),cellstr(varargin{2}));
		else
			error('[Template] Badly formed handles.');
		end
	case 'block'
		error(nargchk(4,5,nargin));
		tpl = loadtpl(tpl,varargin{1});
		if nargin == 4
			name = varargin{2};
		else
			name = varargin{3};
		end
		str = get(tpl,'var',varargin{1});
		blk = '';
		%- look like this (keep the same number (1) of spaces between characters!)
		%  <!-- BEGIN ??? -->
		%  <!-- END ??? -->

		%%%%%%%%%%%%%%%%%%%%%%%%% WIH REGEXP %%%%%%%%%%%%%%%%%%%%%%%%
		% reg = ['<!--\s+BEGIN ' varargin{2} '\s+-->(.*)\n\s*<!--\s+END ' varargin{2} '\s+-->'];
		% [b, e] = regexp(str,reg,'once');
		% if ~isempty(b), blk = str(b:e); end %- should also remove BEGIN and END comments
		% str = regexprep(str,reg,['{' name '}']);
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

		%%%%%%%%%%%%%%%%%%%%%% WIHTOUT REGEXP %%%%%%%%%%%%%%%%%%%%%%%
		indbegin = findstr(str,['<!-- BEGIN ' varargin{2} ' -->']);
		indend   = findstr(str,['<!-- END ' varargin{2} ' -->']);
		if ~isempty(indbegin) & ~isempty(indend)
		   blk = str(indbegin+length(['<!-- BEGIN ' varargin{2} ' -->'])+1:indend-1);
		   str = [str(1:indbegin-1) '{' name '}' str(indend+length(['<!-- END ' varargin{2} ' -->'])+1:end)];
		end
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

		tpl = set(tpl,'var',varargin{2},blk);
		tpl = set(tpl,'var',varargin{1},str);
	case 'var'
		error(nargchk(3,4,nargin));
		if iscellstr(varargin{1})
			for i=1:length(varargin{1})
				ind = find(ismember(tpl.varkeys,varargin{1}{i}));
				if isempty(ind)
					tpl.varkeys{end+1} = varargin{1}{i};
					if nargin == 4
						tpl.varvals{end+1} = varargin{2}{i};
					else
						tpl.varvals{end+1} = '';
					end
				else
					tpl.varvals{ind} = varargin{2}{i};
				end
			end
		elseif ischar(varargin{1})
			tpl = set(tpl,'var',cellstr(varargin{1}),cellstr(varargin{2}));
		else
			error('[Template] Badly formed variable names.');
		end
	otherwise
		error('[Template] Unknown action to perform.');
end

author: zouchuhang

repo: LayoutNet


text:
function plotOptStates(model, nlp, gait, indices)

    if nargin < 5
        indices = 1:length(model.Joints);
    else
        if isempty(indices), return; end
    end

    joint_names = {model.Joints.Name};

    if length(gait)==1 && isa(nlp,'TrajectoryOptimization')

        t = gait.tspan;

        x_opt = gait.states.x;
        x_lb  = [nlp.OptVarTable.x.LowerBound];
        x_ub  = [nlp.OptVarTable.x.UpperBound];

        dx_opt = gait.states.dx;
        dx_lb  = [nlp.OptVarTable.dx.LowerBound];
        dx_ub  = [nlp.OptVarTable.dx.UpperBound];

        ddx_opt = gait.states.ddx;
        ddx_lb  = [nlp.OptVarTable.ddx.LowerBound];
        ddx_ub  = [nlp.OptVarTable.ddx.UpperBound];
    else
        cont_domain_idx = find(cellfun(@(x)isa(x,'ContinuousDynamics'),{nlp.Phase.Plant}));

        t = [];
        x_opt = []; x_lb = [];  x_ub = [];
        dx_opt = []; dx_lb = [];  dx_ub = [];
        ddx_opt = []; ddx_lb = [];  ddx_ub = [];

        for j=cont_domain_idx
            t = [t,gait(j).tspan];

            x_opt = [x_opt,gait(j).states.x];
            x_lb  = [x_lb,[nlp.Phase(j).OptVarTable.x.LowerBound]];
            x_ub  = [x_ub,[nlp.Phase(j).OptVarTable.x.UpperBound]];

            dx_opt = [dx_opt,gait(j).states.dx];
            dx_lb  = [dx_lb,[nlp.Phase(j).OptVarTable.dx.LowerBound]];
            dx_ub  = [dx_ub,[nlp.Phase(j).OptVarTable.dx.UpperBound]];

            ddx_opt = [ddx_opt,gait(j).states.ddx];
            ddx_lb  = [ddx_lb,[nlp.Phase(j).OptVarTable.ddx.LowerBound]];
            ddx_ub  = [ddx_ub,[nlp.Phase(j).OptVarTable.ddx.UpperBound]];
        end
    end

    ax = [];
    for i=indices
        f = figure;clf;
        set(f, 'WindowStyle', 'docked');
        ax = [ax, subplot(3, 1, 1)]; %#ok<*AGROW>
        hold on;
        plot(t, x_opt(i,:), 'b');
        plot(t, x_lb(i,:), 'r--');
        plot(t, x_ub(i,:), 'g--');

        title('Joint Displacement');
        legend('q', 'lb', 'ub');

        ax = [ax, subplot(3, 1, 2)];
        hold on;
        plot(t, dx_opt(i,:), 'b');
        plot(t, dx_lb(i,:), 'r--');
        plot(t, dx_ub(i,:), 'g--');
        title('Joint Velocity');
        legend('dq', 'lb', 'ub');

        ax = [ax, subplot(3, 1, 3)];
        hold on;
        plot(t, ddx_opt(i,:), 'b');
        plot(t, ddx_lb(i,:), 'r--');
        plot(t, ddx_ub(i,:), 'g--');

        title('Joint Acceleration');
        legend('ddq', 'lb', 'ub');

        f.Name = [joint_names{i},'_state'];
    end

    linkaxes(ax, 'x');

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function [E1,IESM_NewSpin,E2,DelE] = EnergyChange_3D_QPOTTS(IESM,Q);

E1 = Energy1_3D_QPOTTS(IESM);
IESM_NewSpin = IndexElementStateMatrix_3D_QPOTTS_NewSpin(IESM,Q);
E2 = Energy2_3D_QPOTTS(IESM_NewSpin);
DelE = E1-E2;

path: github-repos/MATLAB/Sable-mcbench-benchmarks/mcbench-benchmarks-ba13b2f0296ef49491b95e3f984c7c41fccdb6d8/34985-monte-carlo-simulation-of-three-dimensional-grain-growth-code-version-no-1-basic/Monte Carlo Simulation Q-state Potts model 3D square-lattice - microstructure/EnergyChange_3D_QPOTTS.m


text:
%
% This code calculates the interevent distances and the correlation integral
% of a given earthquake distribution.
% Francesco Pacchiani 3/2000
%
%
% Calculation of the 3D distances between all possible pairs
% (combination of n epicenters taken 2 at a time) of earthquakes of
% the given dataset.
%
%
% Variables
%
N = size(E,1);				% N= # of events in the catalogue; E= Earthquake catalogue
pairdist = []; 			% pairdist= Vector of interevent distances
j = nchoosek(N,2);			% j= # of interevent distances calculated
pairdist = zeros(j,1);
k = 0;
%E.Latitude= (max(Da(:,2))+min(Da(:,2)))/2;
%
%
% Calculation of the interevent distances in 2D plus the depths differences.
%
%
for i = 1:(N-1)

    lon1 = repmat(E(i,1), [(N-i),1]);
    lat1 = repmat(E(i,2), [(N-i),1]);

    lon2 = E((i+1):end, 1);
    lat2 = E((i+1):end, 2);

    pairdist(k+1:k + size(lon1, 1)) = distance(lat1,lon1,lat2,lon2);

    k = k + size(lon1,1);

end

clear i j k;
%
%
% Conversion of the interevent distances from degrees to kilometers and calculates
% the interevent distances in three dimensions.
%
%
if dtokm == 1
    pairdist = pairdist.*111;
end
%
%
% Calculation of the correlation integral using as input the
% pair distances computed above.
%
%
% Variables
%
d = 2;						%d = the dimension of the embedding volume.
rmax = max(pairdist);
rmin = min(pairdist);

if rmin == 0
    rmin = 0.01;
end

lrmin = log10(rmin);
lrmax = log10(max(pairdist));

%u = (log10(rmin):0.15:log10(rmax))';
%
% Defining the distance vector r in order that on the
% log-log graph all the points plot at equal distances from one another.
%
r = (logspace(lrmin, lrmax, 50))';
%r = zeros(size(u,1),1);
%r = 10.^u;
%
%
corint = [];						% corint= Vector of ?cumulative? correlation integral values for increasing interevent radius
corint = zeros(size(r,1),1);
k = 1;

for i = 1:size(r,1)

    j = [];
    j = pairdist < r(i);
    corint (k,1) = (2/(N*(N-1)))*sum(j);
    k = k + 1;

end

clear i j k;
%
%
% Plotting of the correlation integral in function of the interevent
% distance r.
%
%
dofdnofig;

author: CelsoReyes


text:
function [resliced] = ft_volumereslice(cfg, mri)

% FT_VOLUMERESLICE flips, permutes, interpolates and/or reslices a volume along the
% principal axes of the coordinate system according to a specified resolution.
%
% Use as
%   mri = ft_volumereslice(cfg, mri)
% where the input MRI should be a single anatomical or functional MRI volume that
% results from FT_READ_MRI or FT_VOLUMEREALIGN. You can visualize the the input and
% output using FT_SOURCEPLOT.
%
% The configuration structure can contain
%   cfg.method     = string, 'flip', 'nearest', 'linear', 'cubic' or 'spline' (default = 'linear')
%   cfg.downsample = integer number (default = 1, i.e. no downsampling)
%
% If you specify the method as 'flip', it will only permute and flip the volume, but
% not perform any interpolation. For the other methods the input volumetric data will
% also be interpolated on a regular voxel grid.
%
% For the interpolation methods you should specify
%   cfg.resolution = number, in units of distance (e.g. mm)
%   cfg.xrange     = [min max], in units of distance (e.g. mm)
%   cfg.yrange     = [min max], in units of distance (e.g. mm)
%   cfg.zrange     = [min max], in units of distance (e.g. mm)
% or alternatively with
%   cfg.dim        = [nx ny nz], size of the volume in each direction
%
% If the input MRI has a coordsys-field and you don't specify explicit the
% xrange/yrange/zrange, the centre of the volume will be shifted (with respect to the
% origin of the coordinate system), for the brain to fit nicely in the box.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_VOLUMEREALIGN, FT_VOLUMEDOWNSAMPLE, FT_SOURCEINTERPOLATE, FT_SOURCEPLOT

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

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

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

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

% check if the input data is valid for this function and ensure that the structures correctly describes a volume
if isfield(mri, 'inside')
  mri = ft_checkdata(mri, 'datatype', 'volume', 'feedback', 'yes', 'hasunit', 'yes', 'insidestyle', 'logical');
else
  mri = ft_checkdata(mri, 'datatype', 'volume', 'feedback', 'yes', 'hasunit', 'yes');
end

% set the defaults
cfg.method     = ft_getopt(cfg, 'method', 'linear');
cfg.downsample = ft_getopt(cfg, 'downsample', 1);

if isequal(cfg.method, 'flip')
  % these do not apply when flipping
  cfg = ft_checkconfig(cfg, 'forbidden', {'resolution', 'xrange', 'yrange', 'zrange', 'dim'});
else
  % these only applies when interpolating
  cfg.resolution = ft_getopt(cfg, 'resolution', 1 * ft_scalingfactor('mm', mri.unit)); % default is 1 mm, but the actual number depends on the units. See bug2906
  cfg.xrange     = ft_getopt(cfg, 'xrange', []);
  cfg.yrange     = ft_getopt(cfg, 'yrange', []);
  cfg.zrange     = ft_getopt(cfg, 'zrange', []);
  cfg.dim        = ft_getopt(cfg, 'dim', []);

  if isfield(mri, 'coordsys')
    % use some prior knowledge to optimize the location of the bounding box
    % with respect to the origin of the coordinate system
    switch mri.coordsys
      case {'ctf', '4d', 'bti', 'eeglab'}
        xshift = 30./cfg.resolution;
        yshift = 0;
        zshift = 40./cfg.resolution;
      case {'neuromag', 'itab'}
        xshift = 0;
        yshift = 30./cfg.resolution;
        zshift = 40./cfg.resolution;
      case {'acpc', 'spm', 'mni', 'tal'}
        ft_warning('FIXME, the bounding box needs a better default');
        xshift = 0;
        yshift = 0;
        zshift = 0;
      otherwise
        xshift = 0;
        yshift = 0;
        zshift = 0;
    end
  else % if no coordsys is present
    xshift = 0;
    yshift = 0;
    zshift = 0;
  end

  if ~isempty(cfg.dim)
    xrange = [-cfg.dim(1)/2+0.5 cfg.dim(1)/2-0.5] * cfg.resolution + xshift;
    yrange = [-cfg.dim(2)/2+0.5 cfg.dim(2)/2-0.5] * cfg.resolution + yshift;
    zrange = [-cfg.dim(3)/2+0.5 cfg.dim(3)/2-0.5] * cfg.resolution + zshift;
  else % if no cfg.dim is specified, use defaults
    range = [-127.5 127.5] * cfg.resolution; % 255 mm^3 bounding box, assuming human brain
    xrange = range + xshift;
    yrange = range + yshift;
    zrange = range + zshift;
  end

  % if ranges have not been specified by the user
  if isempty(cfg.xrange)
    cfg.xrange = xrange;
  end
  if isempty(cfg.yrange)
    cfg.yrange = yrange;
  end
  if isempty(cfg.zrange)
    cfg.zrange = zrange;
  end

end % if method~=fip

if cfg.downsample~=1
  % optionally downsample the anatomical and/or functional volumes
  tmpcfg = keepfields(cfg, {'downsample', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  mri = ft_volumedownsample(tmpcfg, mri);
  % restore the provenance information
  [cfg, mri] = rollback_provenance(cfg, mri);
end

if strcmp(cfg.method, 'flip')
  % this uses some private functions that change the volumes and the transform
  resliced = volumepermute(mri); % this makes the transform approximately diagonal
  flipvec = false(1,3);
  flipvec(1) = resliced.transform(1,1)<0;
  flipvec(2) = resliced.transform(2,2)<0;
  flipvec(3) = resliced.transform(3,3)<0;
  resliced = volumeflip(resliced, flipvec); % this flips along each of the dimensions
  if ~isequal(mri.transform, resliced.transform)
    ft_info('flipped the volume to make it consistent with the axes of the coordinate system');
  end

else
  % compute the desired grid positions
  xgrid = cfg.xrange(1):cfg.resolution:cfg.xrange(2);
  ygrid = cfg.yrange(1):cfg.resolution:cfg.yrange(2);
  zgrid = cfg.zrange(1):cfg.resolution:cfg.zrange(2);

  resliced           = [];
  resliced.dim       = [length(xgrid) length(ygrid) length(zgrid)];
  resliced.transform = translate([cfg.xrange(1) cfg.yrange(1) cfg.zrange(1)]) * scale([cfg.resolution cfg.resolution cfg.resolution]) * translate([-1 -1 -1]);
  resliced.anatomy   = zeros(resliced.dim, 'int8');
  resliced.unit      = mri.unit;

  % these take a lot of memory
  clear xgrid ygrid zgrid

  % these are the same in the resliced as in the input anatomical MRI
  if isfield(mri, 'coordsys')
    resliced.coordsys = mri.coordsys;
  end

  fprintf('reslicing from [%d %d %d] to [%d %d %d]\n', mri.dim(1), mri.dim(2), mri.dim(3), resliced.dim(1), resliced.dim(2), resliced.dim(3));

  % determine the fields to reslice
  fn = parameterselection('all', mri);

  % the actual work is being done by ft_sourceinterpolate
  % this interpolates the real volume on the resolution that is defined for the resliced volume
  tmpcfg                = [];
  tmpcfg.parameter      = fn;
  tmpcfg.interpmethod   = cfg.method;
  resliced              = ft_sourceinterpolate(tmpcfg, mri, resliced);
  resliced.cfg.previous = resliced.cfg.previous{1}; % the 2nd input is a dummy variable
  cfg.method = resliced.cfg.interpmethod;           % remember the method that was used
  % restore the provenance information
  [cfg, resliced] = rollback_provenance(cfg, resliced);

  % remove fields that were not present in the input
  % this applies specifically to the 'inside' field that may have been added by ft_sourceinterpolate
  resliced = keepfields(resliced, fieldnames(mri));

  % convert any non-finite values to 0 to avoid problems later on
  for i=1:numel(fn)
    resliced.(fn{i})(~isfinite((fn{i}))) = 0;
  end

end % if method=flip or interpolate

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   mri
ft_postamble provenance resliced
ft_postamble history    resliced
ft_postamble savevar    resliced


text:
function [H W] = createHashTable(K,b,t)
%% function [H W] = createHashTable(K,b,t)
%  Create a hash table for kernelized LSH.
%
%  Inputs: K (kernel matrix)
%          b (number of bits)
%          t (number of Gaussian approximation elements)
%
%  Outputs: H (hash table over elements of K, size p x b)
%           W (weight matrix for computing hash keys over queries, size p x b)
%%

[p,p] = size(K);
%center the kernel matrix
K = K - (1/p)*(K*ones(p,1))*ones(1,p) - (1/p)*ones(p,1)*(ones(1,p)*K) + (1/p^2)*sum(sum(K));
if nargin < 3
    t = min(floor(p/4), 30);
end

[V_K,D_K] = eig(K);
d_k = diag(D_K);
ind_k = find(d_k > 1e-8);
d_k(ind_k) = d_k(ind_k).^(-1/2);
K_half = V_K*diag(d_k)*V_K';

%create indices for the t random points for each hash bit
%then form weight matrix
for i = 1:b
    rp = randperm(p);
    I_s(i,:) = rp(1:t);
    e_s = zeros(p,1);
    e_s(I_s(i,:)) = 1;
    W(:,i) = sqrt((p-1)/t)*K_half*e_s;
end
H = (K*W)>0;
W = real(W);

author: ZJULearning

repo: MatlabFunc


text:
function mrSESSION = UpdateSessionFunctionals(mrSESSION,scanParams,doCrop);
% mrSESSION = UpdateSessionFunctionals(mrSESSION,scanParams,doCrop);
%
% Set the fields of mrSESSION.functionals
%
% djh & dbr, 9/2001
% ras, 4/04: fixed (in a somewhat ugly way) a bug that
% occurs when you don't want to crop -- it sometimes
% ends up setting the cropSize off by 0.5, for some
% really weird reason I don't get.
if ~exist('doCrop','var')   doCrop = 1;     end

if isfield(mrSESSION, 'functionals')
    nOldScans = length(mrSESSION.functionals);
    nNewScans = length(scanParams);
    if nOldScans > 0
        oldNames = {mrSESSION.functionals(:).PfileName};
    else
        oldNames = {};
    end
    for iNew=1:nNewScans
        newName = scanParams(iNew).PfileName;
        if ~any(strcmp(newName, oldNames))
            % If the Pfile name is new, append it
            mrSESSION.functionals = [mrSESSION.functionals scanParams(iNew)];
        end
    end
else
    mrSESSION.functionals = scanParams;
end

% Update functionals.crop & cropSize
% (This is very confusing, and doesn't
% seem to work if you don't want to crop:
% e.g. if your crop is [1 1; 128 128], it
% makes the cropSize [128.5 128.5], and
% breaks everything.
%   Therefore, if doCrop is not selected,
% don't do this). -ras, 04/04
if doCrop
    ipCrop = mrSESSION.inplanes.crop;
    ipCrop(1,:) = ipCrop(1,:)-1;
    ipFullSize = mrSESSION.inplanes.fullSize;
    for iScan = 1:length(mrSESSION.functionals)
        funcFullSize = mrSESSION.functionals(iScan).fullSize;
        sizeRatio = ipFullSize ./ funcFullSize;
        crop(1,:) = ipCrop(1,:)./sizeRatio;
        crop(2,:) = ipCrop(2,:)./sizeRatio;
        if any(sizeRatio - floor(sizeRatio) <0)
            FatalInitError('Crop is busted, inplane size not a multiple of functional size. Start over and redo the crop.');
        end
        if any(crop - floor(crop) <0)
            FatalInitError('Crop is busted, inplane crop not a multiple of sizeRatio. Start over and redo the crop.');
        end
        crop(1,:) = crop(1,:) + 1;
        mrSESSION.functionals(iScan).crop = crop;
        cropSize = diff(fliplr(crop)) + 1;
        mrSESSION.functionals(iScan).cropSize = cropSize;
    end
else
    for iScan = 1:length(mrSESSION.functionals)
        cropSize = mrSESSION.functionals(iScan).crop(2,:);
        mrSESSION.functionals(iScan).cropSize = cropSize;
    end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [u,edge,eqn,info] = Maxwellsaddle(node,elem,HB,pde,bdFlag,option)
%% MAXWELL Maxwell equation: lowest order edge element.
%
% u = Maxwellsaddle(node,elem,HB,pde,bdFlag) produces the lowest order edge
%   element approximation of the electric field of the time harmonic
%   Maxwell equation.
%
%                     curl(mu^(-1)curl u)  = J    in \Omega,
%                                  div  u  = 0
%                                   n � u = n � g_D  on \Gamma_D,
%                     n � (mu^(-1)curl u) = n � g_N  on \Gamma_N.
%
% based on the weak formulation
%
% (mu^{-1}curl u, curl v) +(grad p ,v)= (J,v) - <n � g_N,v>_{\Gamma_N}.
% (u,grad q)                          = 0
% Assume P\in H_0^1.
% The data of the equation is enclosed in the pde structure:
%   - pde.mu      : permeability, i.e., magnetic constant/tensor
%   - pde.omega   : wave number
%   - pde.J       : current density
%   - pde.g_D     : Dirichlet boundary condition
%   - pde.g_N     : Neumann boundary condition
%
% The mesh is given by (node,elem) and HB is needed for fast solvers. The
% boundary faces is specified by bdFlag; see bddoc.
%
% The function Maxwell assembes the matrix equation (A-M)*u = b and solves
% it by the direct solver (small size dof <= 2e3) or the HX preconditioned
% Krylov iterative methods (large size dof > 2e3).
%
% u = Maxwellsaddle(node,elem,HB,pde,bdFlag,option) specifies the solver options.
%   - option.solver == 'direct': the built in direct solver \ (mldivide)
%   - option.solver == 'mg':     multigrid-type solvers mg is used.
%   - option.solver == 'notsolve': the solution u = u_D.
% The default setting is to use the direct solver for small size problems
% and multigrid solvers for large size problems. For more options on the
% multigrid solver mg, type help mg.
%
% [u,edge] = Maxwellsaddle(node,elem,pde,bdEdge) returns also the edge array
% which is essential for edge elements.
%
% [u,edge,eqn] = Maxwellsaddle(node,elem,pde,bdEdge) returns also the equation
% structure eqn, which includes:
% - eqn.A: matrix for differential operator;
% - eqn.M: mass matrix;
% - eqn.f: right hand side
% - eqn.g: vector enclosed the Neumann boundary condition
%
% [u,edge,eqn,info] = Maxwellsaddle(node,elem,pde,bdEdge) returns also the
% information on the assembeling and solver, which includes:
% - info.assembleTime: time to assemble the matrix equation
% - info.solverTime:   time to solve the matrix equation
% - info.itStep:       number of iteration steps for the mg solver
% - info.error:        l2 norm of the residual b - A*u
% - info.flag:         flag for the mg solver.
%   flag = 0: converge within max iteration
%   flag = 1: iterated maxIt times but did not converge
%   flag = 2: direct solver
%   flag = 3: no solve
%
% Example
%   cubeMaxwell
%
% See also Maxwell1, Maxwell2, cubeMaxwell, mgMaxwell
%
% Reference page in Help browser
%       Maxwelldoc
%
% Created by Jie Zhou based on Maxwell(node,elem,pde,bdEdge) on
% 09,Sep,2013.
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Set up optional input arguments
if ~exist('bdFlag','var'), bdFlag = []; end
if ~exist('option','var'), option = []; end

%% Sort elem to ascend ordering
[elem,bdFlag] = sortelem3(elem,bdFlag);

%% Construct Data Structure
[elem2dof,edge] = dof3edge(elem);
locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];
N = size(node,1);   NT = size(elem,1);  Ne = size(edge,1);

%% Compute coefficients
if ~isfield(pde,'mu'), pde.mu = 1; end
if ~isempty(pde.mu) && isnumeric(pde.mu)
    mu = pde.mu;                % mu is an array
else                            % mu is a function
    center = (node(elem(:,1),:) + node(elem(:,2),:) + ...
              node(elem(:,3),:) + node(elem(:,4),:))/4;
    mu = pde.mu(center);
end
if ~isfield(pde,'epsilon'), pde.epsilon = 0; end
if ~isempty(pde.epsilon) && isnumeric(pde.epsilon)
    epsilon = pde.epsilon;      % epsilon is an array
else                            % epsilon is a function
    center = (node(elem(:,1),:) + node(elem(:,2),:) + ...
              node(elem(:,3),:) + node(elem(:,4),:))/4;
    epsilon = pde.epsilon(center);
end
if isfield(pde,'omega')
    omega = pde.omega;
else
    omega = 1;
end
epsilon = omega^2*epsilon;

tstart = tic;
%% Element-wise basis
% edge indices of 6 local bases:
% [1 2], [1 3], [1 4], [2 3], [2 4], [3 4]
% phi = lambda_iDlambda_j - lambda_jDlambda_i;
% curl phi = 2*Dlambda_i � Dlambda_j;
[Dlambda,volume] = gradbasis3(node,elem);
curlPhi(:,:,6) = 2*mycross(Dlambda(:,:,3),Dlambda(:,:,4),2);
curlPhi(:,:,1) = 2*mycross(Dlambda(:,:,1),Dlambda(:,:,2),2);
curlPhi(:,:,2) = 2*mycross(Dlambda(:,:,1),Dlambda(:,:,3),2);
curlPhi(:,:,3) = 2*mycross(Dlambda(:,:,1),Dlambda(:,:,4),2);
curlPhi(:,:,4) = 2*mycross(Dlambda(:,:,2),Dlambda(:,:,3),2);
curlPhi(:,:,5) = 2*mycross(Dlambda(:,:,2),Dlambda(:,:,4),2);
DiDj = zeros(NT,4,4);
for i = 1:4
    for j = i:4
        DiDj(:,i,j) = dot(Dlambda(:,:,i),Dlambda(:,:,j),2);
        DiDj(:,j,i) = DiDj(:,i,j);
    end
end

%% Assemble matrices
ii = zeros(21*NT,1); jj = zeros(21*NT,1);
sA = zeros(21*NT,1); sM = zeros(21*NT,1);
index = 0;
for i = 1:6
    for j = i:6
        % local to global index map
        % curl-curl matrix
        Aij = dot(curlPhi(:,:,i),curlPhi(:,:,j),2).*volume./mu;
        ii(index+1:index+NT) = double(elem2dof(:,i));
        jj(index+1:index+NT) = double(elem2dof(:,j));
        sA(index+1:index+NT) = Aij;
        % mass matrix
        % locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];
        i1 = locEdge(i,1); i2 = locEdge(i,2);
        j1 = locEdge(j,1); j2 = locEdge(j,2);
        Mij = 1/20*volume.*( (1+(i1==j1))*DiDj(:,i2,j2) ...
                           - (1+(i1==j2))*DiDj(:,i2,j1) ...
                           - (1+(i2==j1))*DiDj(:,i1,j2) ...
                           + (1+(i2==j2))*DiDj(:,i1,j1));
        %Mij = Mij.*epsilon;
        sM(index+1:index+NT) = Mij;
        index = index + NT;
    end
end
clear curlPhi % clear large size data
diagIdx = (ii == jj);   upperIdx = ~diagIdx;
A = sparse(ii(diagIdx),jj(diagIdx),sA(diagIdx),Ne,Ne);
AU = sparse(ii(upperIdx),jj(upperIdx),sA(upperIdx),Ne,Ne);
A = A + AU + AU';
M = sparse(ii(diagIdx),jj(diagIdx),sM(diagIdx),Ne,Ne);
MU = sparse(ii(upperIdx),jj(upperIdx),sM(upperIdx),Ne,Ne);
M = M + MU + MU';
grad =icdmat(double(edge),[-1,1]);
G    = grad'*M;

% bigA = A + M;

%% Assemble right hand side
f = zeros(Ne,1);
if ~isfield(pde,'J') || (isfield(pde,'J') && isreal(pde.J) && all(pde.J==0))
    pde.J = [];
end
if ~isfield(option,'fquadorder')
    option.fquadorder = 2;   % default order is 3
end
if isfield(pde,'J') && ~isempty(pde.J)
    [lambda,w] = quadpts3(option.fquadorder);
    nQuad = size(lambda,1);
    bt = zeros(NT,6);
    for p = 1:nQuad
        % quadrature points in the x-y-z coordinate
        pxyz = lambda(p,1)*node(elem(:,1),:) ...
             + lambda(p,2)*node(elem(:,2),:) ...
             + lambda(p,3)*node(elem(:,3),:) ...
             + lambda(p,4)*node(elem(:,4),:);
        Jp = pde.J(pxyz);
    %   locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];
        for k = 1:6
            i = locEdge(k,1); j = locEdge(k,2);
            % phi_k = lambda_iDlambda_j - lambda_jDlambda_i;
            phi_k = lambda(p,i)*Dlambda(:,:,j)-lambda(p,j)*Dlambda(:,:,i);
            rhs = dot(phi_k,Jp,2);
            bt(:,k) = bt(:,k) + w(p)*rhs;
        end
    end
    bt = bt.*repmat(volume,1,6);
    f = accumarray(elem2dof(:),bt(:),[Ne 1]);
end
clear pxyz Jp bt rhs phi_k

%% Set up solver
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if Ne <= 1e4  % Direct solver for small size systems
        option.solver = 'direct';
    else            % Multigrid-type  solver for large size systems
        option.solver = 'cg';
    end
end
solver = option.solver;

%% Assembeling corresponding matrices for HX preconditioner
if ~strcmp(solver,'direct') && ~strcmp(solver,'nosolve')
    AP = sparse(N,N);  % AP = - div(mu^{-1}grad) + |Re(epsilon)| I
    BP = sparse(N,N);  % BP = - div(|Re(epsilon)|grad)
    for i = 1:4
        for j = i:4
            temp = DiDj(:,i,j).*volume;
            Aij = 1./mu.*temp;
            Bij = abs(real(epsilon)).*temp;
            Mij = 1/20*abs(real(epsilon)).*volume;
            if (j==i)
                AP = AP + sparse(elem(:,i),elem(:,j),Aij+2*Mij,N,N);
                BP = BP + sparse(elem(:,i),elem(:,j),Bij,N,N);
            else
                AP = AP + sparse([elem(:,i);elem(:,j)],[elem(:,j);elem(:,i)],...
                                 [Aij+Mij; Aij+Mij],N,N);
                BP = BP + sparse([elem(:,i);elem(:,j)],[elem(:,j);elem(:,i)],...
                                 [Bij; Bij],N,N);
            end
        end
    end
end
clear Aij Bij Mij

%% Boundary conditions
if ~isfield(pde,'g_D'), pde.g_D = []; end
if ~isfield(pde,'g_N'), pde.g_N = []; end
if ~isfield(pde,'g_R'), pde.g_R = []; end
if (isempty(pde.g_D) && isempty(pde.g_N) && isempty(pde.g_R))
    % no boundary data is given = homogenous Neumann boundary condition
    bdFlag = [];
end

%% Part 1: Find Dirichlet dof and modify the matrix
% Find Dirichlet boundary dof: fixedDof
isBdEdge = [];
if isempty(bdFlag) && ~isempty(pde.g_D) && isempty(pde.g_N)
    % Dirichlet boundary condition only
    bdFlag = setboundary3(node,elem,'Dirichlet');
end
if ~isempty(bdFlag)
    % Find boundary edges and nodes
    isBdEdge = false(Ne,1);
    isBdNode = false(N,1);
    isBdEdge(elem2dof(bdFlag(:,1) == 1,[4,5,6])) = true;
    isBdEdge(elem2dof(bdFlag(:,2) == 1,[2,3,6])) = true;
    isBdEdge(elem2dof(bdFlag(:,3) == 1,[1,3,5])) = true;
    isBdEdge(elem2dof(bdFlag(:,4) == 1,[1,2,4])) = true;
    bdEdge = edge(isBdEdge,:);
    isBdNode(bdEdge) = true;
end
% modify the matrix to include the Dirichlet boundary condition
if any(isBdEdge)  % contains Dirichlet boundary condition
    bdidx = zeros(Ne,1);
    bdidx(isBdEdge) = 1;
    Tbd = spdiags(bdidx,0,Ne,Ne);
    Te = spdiags(1-bdidx,0,Ne,Ne);
    AD = Te*(A - epsilon*M)*Te + Tbd;
    if ~strcmp(solver,'direct') && ~strcmp(solver,'nosolve')
        % modify the corresponding Poisson matrix
        bdidx = zeros(N,1);
        bdidx(isBdNode) = 1;
        Tbd = spdiags(bdidx,0,N,N);
        Tv = spdiags(1-bdidx,0,N,N);
        AP = Tv*AP*Tv + Tbd;
        BP = Tv*BP*Tv + Tbd;
       % G  = Te*G*Tv  ;
    end
else      % pure Neumann boundary condition
    AD = A - epsilon*M  ;
    if ~strcmp(solver,'direct') && ~strcmp(solver,'nosolve')
       BP = BP + 1e-8*speye(N);  % make B non-singular
    end
end

%% Part 2: Find boundary edges and modify the load b
g = zeros(Ne,1);
% Find Neumann boundary faces
if isempty(bdFlag) && (~isempty(pde.g_N) || ~isempty(pde.g_R))
    bdFlag = setboundary3(node,elem,'Neumann');
end
% non-zero Neumann boundary condition
if ~isempty(bdFlag) && ~isempty(pde.g_N)
    % face 1
    isBdElem = find(bdFlag(:,1) == 2); %#ok<*NASGU>
    face = [2 3 4]; face2locdof = [6 5 4];
    if ~isempty(isBdElem)
        bdb = bdfaceintegral(isBdElem,face,face2locdof);
        g = bdb;
    end
    % face 2
    isBdElem = find(bdFlag(:,2) == 2);
    face = [1 4 3]; face2locdof = [6 2 3];
    if ~isempty(isBdElem)
        bdb = bdfaceintegral(isBdElem,face,face2locdof);
        g = g + bdb;
    end
    % face 3
    isBdElem = find(bdFlag(:,3) == 2);
    face = [1 2 4]; face2locdof = [5 3 1];
    if ~isempty(isBdElem)
        bdb = bdfaceintegral(isBdElem,face,face2locdof);
        g = g + bdb;
    end
    % face 4
    isBdElem = find(bdFlag(:,4) == 2);
    face = [1 3 2]; face2locdof = [4 1 2];
    if ~isempty(isBdElem)
        bdb = bdfaceintegral(isBdElem,face,face2locdof);
        g = g + bdb;
    end
    f = f - g;
end
% nonzero Dirichlet boundary condition
u = zeros(Ne,1);
if ~isempty(bdEdge) && ~isempty(pde.g_D) && ...
   ~(isnumeric(pde.g_D) && all(pde.g_D == 0))
    % else no bddof or g_D = 0 (no modification needed)
    if (isnumeric(pde.g_D) && length(pde.g_D) == Ne)
        u(isBdEdge) = pde.g_D(isBdEdge);
    else
        u(isBdEdge) = edgeinterpolate(pde.g_D,node,bdEdge);
    end
    f = f - (A - epsilon*M)*u;
    f(isBdEdge) = u(isBdEdge);
end
%% We always assume the Langrange multer is in zeros boundary condition
p     = zeros(N,1);
g0    = zeros(N,1);
g0    = -G*u;

%% Remark
% The order of assign Neumann and Dirichlet boundary condition is
% important to get the right setting of the intersection of Dirichlet and
% Neumann faces.

%% Record assembling time
info.assembleTime = toc(tstart);
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 1
    fprintf('Time to assemble matrix equation %4.2g s\n',info.assembleTime);
end

%% Solve the system of linear equations
% if strcmp(solver,'direct')
%     % exact solver
%     tstart = tic;
%     freeDof = find(~isBdEdge);
%     u(freeDof) = AD(freeDof,freeDof)\f(freeDof);
%     time = toc(tstart); itStep = 0; flag = 2; err = norm(f - AD*u);
% elseif strcmp(solver,'nosolve')
%     eqn = struct('A',A,'M',M,'f',f,'g',g,'bigA',AD,'isBdEdge',isBdEdge);
%     info = [];
%     return;
% else
% %     u0 = edgeinterpolate(pde.g_D,node,edge);
%     u0 = u;
%     option.x0 = u0;
%     [u,flag,itStep,err,time] = mgMaxwell(AD,f,AP,BP,node,elem,edge,HB,isBdEdge,option);
% end
temp  = zeros(N+Ne,1);
bigAD = [AD G';G sparse(N,N)];
freeEdge = find(~isBdEdge);
freeNode = find(~isBdNode);
freeDof  = [freeEdge;Ne+freeNode];
tstart = tic;
temp(freeDof) = bigAD(freeDof,freeDof)\[f(freeEdge);g0(freeNode)];
time = toc(tstart);itStep = 0; flag = 2; err = norm([f(freeEdge);g0(freeNode)] - bigAD(freeDof,freeDof)*temp(freeDof));
u(freeEdge) = temp(freeEdge);
p(freeNode) = temp(Ne + freeNode);

%check the Langrange is correct or not.
normp  = norm(p);
if(normp>1.0/N)
disp('the Langrange multer is wrong')
end

%% Output
eqn = struct('A',A,'M',M,'f',f,'g',g,'bigA',AD,'isBdEdge',isBdEdge);
info = struct('solverTime',time,'itStep',itStep,'error',err,'flag',flag);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions bdfaceintegral
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function bdb = bdfaceintegral(isBdElem,face,face2locdof)
    %% Compute boundary surface integral of lowest order edge element.
    %  bdb(k) = \int_{face} (n�g_N, phi_k) dS

    %% Compute scaled normal
    faceIdx = true(4,1);
    faceIdx(face) = false;
    normal = -3*repmat(volume(isBdElem),1,3).*Dlambda(isBdElem,:,faceIdx);

    %% Data structure
    tetLocEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4]; % edge of a tetrahedral [1 2 3 4]
    face2locEdge = [2 3; 3 1; 1 2]; % edge of the face [1 2 3]

    %% Compute surface integral
    Nbd = length(isBdElem);
    bt = zeros(Nbd,3);
    idx = zeros(Nbd,3,'int32');
    [lambda,w] = quadpts(3); % quadrature order is 3
    nQuad = size(lambda,1);
    for pp = 1:nQuad
        % quadrature points in the x-y-z coordinate
        pxyz = lambda(pp,1)*node(elem(isBdElem,face(1)),:) ...
             + lambda(pp,2)*node(elem(isBdElem,face(2)),:) ...
             + lambda(pp,3)*node(elem(isBdElem,face(3)),:);
        gNp = pde.g_N(pxyz,normal);
        for s = 1:3
            kk = face2locdof(s);
            pidx = face(face2locEdge(s,1))< face(face2locEdge(s,2));
            % phi_k = lambda_iDlambda_j - lambda_jDlambda_i;
            % lambda_i is associated to the local index of the face [1 2 3]
            % Dlambda_j is associtated to the index of tetrahedron
            % - when the direction of the local edge s is consistent with the
            %   global oritentation given in the triangulation,
            %           s(1) -- k(1),  s(2) -- k(2)
            % - otherwise
            %           s(2) -- k(1),  s(1) -- k(2)
            if pidx
                phi_k = lambda(pp,face2locEdge(s,1))*Dlambda(isBdElem,:,tetLocEdge(kk,2)) ...
                      - lambda(pp,face2locEdge(s,2))*Dlambda(isBdElem,:,tetLocEdge(kk,1));
            else
                phi_k = lambda(pp,face2locEdge(s,2))*Dlambda(isBdElem,:,tetLocEdge(kk,2)) ...
                      - lambda(pp,face2locEdge(s,1))*Dlambda(isBdElem,:,tetLocEdge(kk,1));
            end
            rhs = dot(phi_k,gNp,2);
            bt(:,s) = bt(:,s) + w(pp)*rhs; % area is included in normal; see line 28
            idx(:,s) = elem2dof(isBdElem,kk);
        end
    end
    %% Distribute to DOF
    bdb = accumarray(idx(:),bt(:),[Ne 1]);
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

repo: Aero-matlab


text:
function charlier_test ( )

%*****************************************************************************80
%
%% CHARLIER_TEST tests CHARLIER.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 March 2009
%
%  Author:
%
%    John Burkardt
%
  test_num = 5;
  a_test = [ 0.25, 0.5, 1.0, 2.0, 10.0 ];
  n = 5;
  OFFSET = 1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHARLIER_TEST:\n' );
  fprintf ( 1, '  CHARLIER evaluates Charlier polynomials.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       N      A         X        P(N,A,X)\n' );

  for test = 1 : test_num

    n = 5;
    a = a_test(test);

    for j = 0 : 5

      x = j / 2.0;

      value = charlier ( n, a, x );

      fprintf ( 1, '\n' );

      for i = 0 : n

        fprintf ( 1, '  %8d  %8f  %8f  %14f\n', i, a, x, value(i+OFFSET) );

      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function cells2file(inputC,fnameS)
%function cells2file.  Write out cell contents to a file.
%JOD.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

fid=fopen(fnameS,'wt+');
i=1;
for i=1:length(inputC)
  lineS=inputC{i};
  fwrite(fid,lineS,'char');
  fprintf(fid,'\n',[]);
end

fclose(fid);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function lf=openmeeg_helper(pos,vol,splitdim)
%
% Splits the calculation of leadfields in smaller chunks
% pos: positions of dipoles
% vol: volume conductor
% splitdim: splitting factor

lf   = [];
ndip = size(pos,1);
numsol   = fix(ndip./splitdim);
nleft    = rem(ndip,splitdim);
lf = zeros(size(vol.mat,1),ndip*3);
if isfield(vol,'mat')
    for i = 1:splitdim
      dsm = openmeeg_dsm(pos((numsol*(i - 1) + 1):numsol*i,:),vol);
      tmp = vol.mat*dsm;
      lf(:,(numsol*(i - 1)*3 + 1):numsol*i*3) = tmp;
    end
    dsm = openmeeg_dsm(pos(ndip-nleft+1:ndip,:),vol);
    tmp = vol.mat*dsm;
    lf(:,(ndip-nleft)*3+1:ndip*3) = tmp;
else
  error('No system matrix is present, BEM head model not calculated yet')
end


text:
function g = p10_g ( n, x )

%*****************************************************************************80
%
%% P10_G evaluates the gradient for problem 10.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real G(N), the gradient of the objective function.
%
  g = zeros ( n, 1 );

  g(1) = 2.0 * x(1) - 2000000.0 + 2.0 * x(1) * x(2) * x(2) - 4.0 * x(2);

  g(2) = 2.0 * x(2) - 0.000004 + 2.0 * x(1)^2 * x(2) - 4.0 * x(1);

  return
end


text:
function [A,B,Ypredict] = bls(X,Y,nhidden,algorithm)
% BLS  Bottleneck least squares aka sparse orthogonalized least squares
%
% X: ninput x nsamples input data matrix
% Y: noutput x nsamples output data matrix
% nhidden: dimension of hidden
%
% A: noutput x nhidden weight matrix
% B: ninput x nhidden weight matrix

if nargin < 3,
  nhidden = min(size(Y,1),2);
end

if nargin < 4,
  algorithm = 'direct';
end

[ninput,nsamples] = size(X);

Sxx = symmetric(X*X')/nsamples;

switch algorithm

  case 'direct'
    Syx = Y*X'/nsamples;
    Sxy = Syx';
    Sxxinv = symmetric(pinv(Sxx));
    S = symmetric(Syx * Sxxinv * Sxy);
    opts.disp = 0;
    opts.issym = 'true';
    [A,D] = eigs(S,nhidden,'LM',opts);   % first nhidden eigenvectors of S
    B = Sxxinv*Sxy*A;

  case 'indirect'          % same result and subspace as 'direct', but then different linear combination

    Syx = Y*X'/nsamples;
    Sxy = Syx';
    Sxxinv = symmetric(pinv(Sxx));

    % initialize randomly

    B = randn(ninput,nhidden);

    iter = 0;
    maxiter = 1000;
    tol = nhidden*ninput*(1e-10);
    Bold = B;
    while iter < maxiter,
      Z = B'*X;
      %		Szz = symmetric(Z*Z')/nsamples;
      Szz = symmetric(B'*Sxx*B);
      %		Syz = symmetric(Y*Z')/nsamples;
      Syz = Syx*B;
      A = Syz*pinv(Szz);

      B = Sxxinv*Sxy*A*pinv(A'*A);
      if sumsqr(B - Bold) < tol,
        fprintf('done after %d iterations!!\n',iter+1);
        iter = maxiter;
      else
        iter = iter + 1;
        Bold = B;
      end
    end

  case 'sequential'      % same result as 'direct', up to minus sign per hidden unit

    R = Y;
    noutput = size(Y,1);
    A = zeros(noutput,nhidden);
    B = zeros(ninput,nhidden);
    for i=1:nhidden,
      [A(:,i),B(:,i)] = bls(X,R,1,'direct');
      if i < nhidden,
        R = R - A(:,i)*B(:,i)'*X;
      end
    end

end

% transform to make latent variable orthonormal

Szz = symmetric(B'*Sxx*B);
Szzinv = symmetric(pinv(Szz));
sqrtSzzinv = chol(Szzinv(end:-1:1,end:-1:1));    % such that the rescaled B(i,:) is a linear combi of the unscaled B(1:i,:); this only makes sense for the direct method
sqrtSzzinv = sqrtSzzinv(end:-1:1,end:-1:1);
B = B*sqrtSzzinv;
A = A*pinv(sqrtSzzinv);

if nargout > 2,
  Ypredict = A*B'*X;
end

function A = symmetric(A)

A = (A+A')/2;


text:
function [XN,IOTA,EPS,U]=quasi_idempotence(X,K)
%QUASI_IDEMPOTENCE     Connection matrix quasi-idempotence
%
%   [XN,IOTA,EPS,U]=quasi_idempotence(X)
%   [XN,IOTA,EPS,U]=quasi_idempotence(X,K)
%
%   The degree of quasi-idempotence of a matrix represents how close it
%   is to being idempotent, i.e. invariant to squaring. In turn, this
%   reflects how closely related the edges in the graph it corresponds to
%   are to the sums of all triangles between the corresponding nodes,
%   spanning the entirety of the network. This probes a form of
%   self-similarity, intended not as between nodes and modules, but
%   between as edges and triangles (or more generally paths).
%   Networks wherein the edge strengths are weakly related to the
%   nodal strengths have low quasi-idempotence, and networks wherein the
%   edge strengths are significantly influenced by the strengths of the
%   corresponding nodes have high quasi-idempotence. In other words the
%   degree of quasi-idempotence may be viewed as a measure of
%   "collectivity", meaning how strongly individual nodes participate
%   in collective dynamics.
%
%   Inputs:
%       X,
%           undirected weighted/binary connection matrix with
%           non-negative weights,
%       K,
%           number of iterations (optional),
%               K<inf,        iterate a predetermined number of times,
%               K=inf,        iterate until convergence is attained (default).
%
%       Note: see Minati et al. (2017) for a discussion of the issues
%       associated with applying this measure to binary graphs.
%
%   Outputs:
%       XN,
%           the final matrix, from the last iteration,
%       IOTA,
%           the vector of correlation coefficients, one per iteration,
%       EPS,
%           the vector of errors, one per iteration,
%       U,
%           the number of iterations performed.
%
%       Note: see Minati et al. (2017) for a discussion of the
%       significance of IOTA(1) and IOTA(end).
%
%   Example:
%       % Create a random undirected weighted connection matrix
%       N=100;
%       X=rand(N);
%       X=X+X';
%       [XN,IOTA,EPS,U]=quasi_idempotence(X);
%       % observe that IOTA~0
%       disp(IOTA);
%       %
%       % make the connection matrix idempotent by squaring it
%       X=X^2;
%       [XN,IOTA,EPS,U]=quasi_idempotence(X);
%       % observe that IOTA~1
%       disp(IOTA);
%
%   Reference:
%       Minati et al. (2017) CHAOS, 043115
%
%   Ludovico Minati, University of Trento, Trento, Italy
%   and Institute of Nuclear Physics, Polish Academy of Sciences, Krakow,
%   Poland, 2017-2018
%   lminati@ieee.org
%
%   Modification history
%   2018: Original

X=cast(X,'double'); % enforce double-precision format
if any(any(isnan(X))) % bail out if any nan found
    error('found nan, giving up!');
end
if any(any(X<0)) % bail out if any negative elements found
    error('found a negative element, giving up!');
end
if ~exist('K','var')||isempty(K)
    K=inf;
end
N=length(X); % get matrix size
X(eye(N)>0)=0; % null the diagonal in the initial matrix
X=X/norm(X); % normalize to unit norm
XN=X;
mask=triu(ones(N),1)>0; % create mask for superdiagonal elements
U=0; % initialize iterations counter
IOTA=[]; % this vector will contain the correlation coefficients
EPS=inf; % and this will contain the errors
if isinf(K)
    while EPS(end)>eps('double') % iterate until error below precision
        U=U+1; % increase iteration counter
        XN_hat=XN; % save the initial matrix
        XN=XN^2; % square the matrix
        XN=XN/norm(XN); % normalize it again (for numerical reasons)
        IOTA(end+1)=corr(X(mask),XN(mask)); % calculate correlation coefficient
        EPS(end+1)=norm(XN_hat-XN); % calculate error
    end
else
    while U<K % iterate a prescribed number of times
        U=U+1; % increase iteration counter
        XN_hat=XN; % save the initial matrix
        XN=XN^2; % square the matrix
        XN=XN/norm(XN); % normalize it again (for numerical reasons)
        IOTA(end+1)=corr(X(mask),XN(mask)); % calculate correlation coefficient
        EPS(end+1)=norm(XN_hat-XN); % calculate error
    end
end
EPS(1)=[];

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [sts, val] = subsasgn_check(item,subs,val)

% function [sts, val] = subsasgn_check(item,subs,val)
% Do a check for proper assignments of values to fields. This routine
% will be called for derived objects from @cfg_.../subsasgn with the
% original object as first argument and the proposed subs and val fields
% before an assignment is made. It is up to each derived class to
% implement assignment checks for both its own fields and fields
% inherited from cfg_item. If used for assignment checks for inherited
% fields, these must be dealt with in special cases in @cfg_.../subsasgn
%
% This routine is a both a check for cfg_item fields and a fallback
% placeholder in cfg_item if a derived class does not implement its own
% checks. In this case it always returns true. A derived class may also
% check assignments to cfg_item fields (e.g. to enforce specific validity
% checks for .val fields). subsasgn_check of the derived class is called
% before this generic subsasgn_check is called and both checks need to be
% passed.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: subsasgn_check.m 6133 2014-08-07 10:35:08Z volkmar $

rev = '$Rev: 6133 $'; %#ok

sts = true;
switch class(item)
    case {'cfg_item'}
        % do subscript checking for base class
        switch subs(1).subs
            case {'tag', 'name'},
                if ~ischar(val)
                    cfg_message('matlabbatch:check:tagname', '%s: Value must be a string.', subsasgn_checkstr(item,subs));
                    sts = false;
                end
                %case {'val'},
                % No checks. All checks must be run in derived classes.
            case {'help'},
                if isempty(val)
                    val = {};
                elseif ~iscellstr(val)
                    cfg_message('matlabbatch:check:help', '%s: Value must be a cell string.', subsasgn_checkstr(item,subs));
                    sts = false;
                end
            case {'check','def','rewrite_job','preview'},
                if isempty(val)
                    val = [];
                elseif ~subsasgn_check_funhandle(val)
                    cfg_message('matlabbatch:check:funhandle', '%s: Value must be a function or function handle.', subsasgn_checkstr(item,subs));
                    sts = false;
                end
            case {'hidden', 'expanded'},
                if ~islogical(val)
                    cfg_message('matlabbatch:check:hiddenexpanded', '%s: Value must be ''true'' or ''false''.', subsasgn_checkstr(item,subs));
                    sts = false;
                end
        end
    otherwise
        % fall back for derived classes
        sts = true;
end


text:
%TestSuiteTest Unit tests for TestSuite class

classdef TestRunLoggerTest < TestCaseInDir

   methods
      function self = TestRunLoggerTest(name)
         self = self@TestCaseInDir(name, ...
             fullfile(fileparts(which(mfilename)), 'helper_classes'));
      end

      function testTwoPassingTests(self)
         logger = TestRunLogger;
         suite = TestSuite('TwoPassingTests');
         suite.run(logger);

         assertEqual(logger.Log, ...
             {'TestRunStarted', ...
             'TestComponentStarted', ...
             'TestComponentStarted', 'TestComponentFinished', ...
             'TestComponentStarted', 'TestComponentFinished', ...
             'TestComponentFinished', ...
             'TestRunFinished'});

         assertEqual(logger.NumTestCases, 2);
         assertEqual(logger.NumFailures, 0);
         assertEqual(logger.NumErrors, 0);
         assertTrue(isempty(logger.Faults));
      end

      function testFailingTestCase(self)
         logger = TestRunLogger;
         suite = TestSuite('FailingTestCase');
         suite.run(logger);

         assertEqual(logger.Log, ...
             {'TestRunStarted', ...
             'TestComponentStarted', ...
             'TestComponentStarted', 'TestCaseFailure', 'TestComponentFinished', ...
             'TestComponentFinished', ...
             'TestRunFinished'});

         assertEqual(logger.NumTestCases, 1);
         assertEqual(logger.NumFailures, 1);
         assertEqual(logger.NumErrors, 0);
         assertEqual(numel(logger.Faults), 1);
         assertEqual(logger.Faults(1).Type, 'failure');
      end

   end

end

author: ricedsp


text:
function vol=elemvolume(node,elem,option)
%
% vol=elemvolume(node,elem,option)
%
% calculate the volume for a list of simplexes
%
% author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
% date: 2007/11/21
%
% input:
%    node:  node coordinates
%    elem:  element table of a mesh
%    option: if option='signed', the volume is the raw determinant,
%            else, the results will be the absolute values
%
% output:
%    vol:   volume values for all elements
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(size(elem,2)==size(node,2))
	enum=size(elem,1);
	vol=zeros(enum,1);
	acol=ones(3,1);
	for i=1:enum
		e1=det([node(elem(i,:),2),node(elem(i,:),3),acol]);
		e2=det([node(elem(i,:),3),node(elem(i,:),1),acol]);
		e3=det([node(elem(i,:),1),node(elem(i,:),2),acol]);
		vol(i)=sqrt(e1*e1+e2*e2+e3*e3)/2;
	end
	return;
end
dim=size(elem,2);
enum=size(elem,1);
vol=zeros(enum,1);
for i=1:enum
    detmat=[node(elem(i,:),:)';ones(1,dim)];
    vol(i)=det(detmat);
end
if(nargin==3 && strcmp(option,'signed'))
  vol=vol/prod(1:size(node,2));
else
  vol=abs(vol)/prod(1:size(node,2));
end

repo: REST


text:
%  Figure 7.40      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% script to generate Fig. 7.40
clf;
clear all;
np=1;
dp=[1, 0, 0];
nc=8.07*[1, 0.619];
dc=[1, 6.41];
num=conv(np,nc);
den=conv(dp,dc);
w=logspace(-1,2);
sys=tf(num,den);
[mag, ph]=bode(sys,w);
mag1=[mag(:,:); ones(size(mag(:,:)))];
ph1=[ph(:,:); -180*ones(size(ph(:,:)))];
subplot(211), loglog(w,mag1), grid;
xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Fig.7.40 Bode plot for reduced-order controller for 1/s^2')
subplot(212) , semilogx(w,ph1) , grid
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');


text:
function [dat] = set_name(dat,nam)

%   data = set_name(data,name) sets the name of the data_global object to
%                              name
dat.algorithm.name=nam;


text:
function taking_file_input()

end

author: anishLearnsToCode

repo: introduction-to-programming-with-matlab


text:
function rgb_3x3 = image_denoise_rgb_3x3 ( rgb )

%*****************************************************************************80
%
%% IMAGE_DENOISE_RGB_3X3 filters out noise from an RGB image.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, uint8 RGB(:,:,3), the noisy RGB data.
%
%    Output, uint8 RGB_3X3(:,:,3), the RGB data for the filtered image.
%
  r = image_denoise_gray_3x3 ( rgb(:,:,1) );
  g = image_denoise_gray_3x3 ( rgb(:,:,2) );
  b = image_denoise_gray_3x3 ( rgb(:,:,3) );

  rgb_3x3 = zeros ( size ( rgb ) );

  rgb_3x3 = uint8 ( rgb_3x3 );

  rgb_3x3(:,:,1) = r(:,:);
  rgb_3x3(:,:,2) = g(:,:);
  rgb_3x3(:,:,3) = b(:,:);

  return
end


text:
function param = quenchInitParams(varargin)
%
% param = quenchInitParams(varargin)
%
% Initializes the params structure which gets attached to a quench fg. It
% contains information that Quench uses for visualization and for dynamic
% queries. Pulled out of dtiCreateQuenchStats.m
%
% PARAMS STRUCTURE:
%       name    -  Name for a fiber summary statistic. Typically this
%                  will be  something like length or "average FA". THis
%                  statistic is used to limit the range of fibers we see in
%                  the Quench window.
%       uid     -  Unique ID [=1] Subtract 1 so Quench starts from 0 (C-counting)
%       ile     -  Is luminance encoding (used by Quench)
%       icpp    -  Is computed per point [default=1]
%       ivs     -  Is viewable stat (used by Quench)
%       agg     -  Name for fiber summary aggregate statistic (defaults to
%                  param.name
%       lname   -  Local name name for the per point statistic (e.g., TF,
%                  FA, eccentricity)
%       stat    -  A (1 x nFibers) array of statistical values
%
% WEB:
%       mrvBrowseSVN('quenchInitParams');
%
% EXAMPLE USAGE:
%       param = quenchInitParams('name','averageFA');
%
% See also:
%       dtiCreateQuenchStats.m , addStatisticHeader.m
%

%% This statistics parameter field gets attached to the fg.
% It contains information that Quench uses for visualization and for
% dynamic queries. We should probably write a function for initializing
% this parameter field. (see also addStatisticHeader.m)

param       = struct;       % Initialize structure
param.name  = 'nameStat';   % Stat name
param.uid   = 1;            % Subtract 1 so Quench starts from 0 (C-counting)
param.ile   = 1;            % Is luminace encoding
param.icpp  = 1;            % Visualization per point, not per fiber
param.ivs   = 1;            % is_viewable_stat
param.agg   = param.name;   % Stat name
param.lname = 'localStat';  % Local Name
param.stat  = [];

param = mrVarargin( param, varargin );

% If the users only passes in 'name' then set param.agg = param.name since
% they are essentially the same thing. If the user passes in both, or just
% param.agg then leave it be.
if ~isempty(strmatch('name',varargin)) && isempty(strmatch('agg',varargin))
    param.agg = param.name;
end

return


text:
function mean = hypergeometric_mean ( n, m, l )

%*****************************************************************************80
%
%% HYPERGEOMETRIC_MEAN returns the mean of the Hypergeometric PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of balls selected.
%    0 <= N <= L.
%
%    Input, integer M, the number of white balls in the population.
%    0 <= M <= L.
%
%    Input, integer L, the number of balls to select from.
%    0 <= L.
%
%    Output, real MEAN, the mean of the PDF.
%
  mean = n * m / l;

  return
end


text:
function triangle_ncc_rule_test03 ( )

%*****************************************************************************80
%
%% TEST03 tests TRIANGLE_NCC_RULE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST03\n' );
  fprintf ( 1, '  TRIANGLE_NCC_RULE returns the points and weights\n' );
  fprintf ( 1, '  of an NCC rule for the triangle.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  In this test, we simply check that, for each\n' );
  fprintf ( 1, '  quadrature point, the barycentric coordinates\n' );
  fprintf ( 1, '  sum to 1.\n' );

  rule_num = triangle_ncc_rule_num ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '      Rule   Suborder    Sum of coordinates\n' );
  fprintf ( 1, '\n' );

  for rule = 1 : rule_num

    suborder_num = triangle_ncc_suborder_num ( rule );

    [ suborder_xyz, suborder_w ] = triangle_ncc_subrule ( rule, suborder_num );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %8d  %8d\n', rule, suborder_num );
    for suborder = 1 : suborder_num
      xyz_sum = sum ( suborder_xyz(1:3,suborder) );
      fprintf ( 1, '                      %25e\n', xyz_sum );
    end

  end

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   demo script for mesh generation from binary volumetric image
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% preparation
% user must add the path of iso2mesh to matlab path list
% addpath('../');

% user need to add the full path to .../iso2mesh/bin directory
% to windows/Linux/Unix PATH environment variable

%% load the sample data
load rat_head.mat

% volimage is a volumetric image such as an X-ray or MRI image
% A,b are registration matrix and vector, respectively
%% perform mesh generation

%% use the alternative 'cgalmesh' method. This will call
% cgalmesher to process labled volume to produce surfaces
% and tetrahedral mesh in a single run.
clear opt
opt.radbound=2;
[node,elem,face]=v2m(uint8(volimage),0.5,opt,100,'cgalmesh');

%% visualize the resulting mesh

plotmesh(node,face(:,1:3));
axis equal;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [y,L] = mdwt(x,h,L)
%    [y,L] = mdwt(x,h,L);
%
%    Function computes the discrete wavelet transform y for a 1D or 2D input
%    signal x using the scaling filter h.
%
%    Input:
%	x : finite length 1D or 2D signal (implicitly periodized)
%       h : scaling filter
%       L : number of levels. In the case of a 1D signal, length(x) must be
%           divisible by 2^L; in the case of a 2D signal, the row and the
%           column dimension must be divisible by 2^L. If no argument is
%           specified, a full DWT is returned for maximal possible L.
%
%    Output:
%       y : the wavelet transform of the signal
%           (see example to understand the coefficients)
%       L : number of decomposition levels
%
%    1D Example:
%       x = makesig('LinChirp',8);
%       h = daubcqf(4,'min');
%       L = 2;
%       [y,L] = mdwt(x,h,L)
%
%    1D Example's  output and explanation:
%
%       y = [1.1097 0.8767 0.8204 -0.5201 -0.0339 0.1001 0.2201 -0.1401]
%       L = 2
%
%    The coefficients in output y are arranged as follows
%
%       y(1) and y(2) : Scaling coefficients (lowest frequency)
%       y(3) and y(4) : Band pass wavelet coefficients
%       y(5) to y(8)  : Finest scale wavelet coefficients (highest frequency)
%
%    2D Example:
%
%       load test_image
%       h = daubcqf(4,'min');
%       L = 1;
%       [y,L] = mdwt(test_image,h,L);
%
%    2D Example's  output and explanation:
%
%       The coefficients in y are arranged as follows.
%
%              .------------------.
%              |         |        |
%              |    4    |   2    |
%              |         |        |
%              |   L,L   |   H,L  |
%              |         |        |
%              --------------------
%              |         |        |
%              |    3    |   1    |
%              |         |        |
%              |   L,H   |  H,H   |
%              |         |        |
%              `------------------'
%
%       where
%            1 : High pass vertically and high pass horizontally
%            2 : Low pass vertically and high pass horizontally
%            3 : High pass vertically and low  pass horizontally
%            4 : Low pass vertically and Low pass horizontally
%                (scaling coefficients)
%
%
%
%
%    See also: midwt, mrdwt, mirdwt
%
%Author: Markus Lang  <lang@jazz.rice.edu>
if exist('OCTAVE_VERSION', 'builtin')
  x = x * 1.0;
  if (exist('L'))
    [y,L] = omdwt(x,h,L);
  else
    [y,L] = omdwt(x,h);
  end
else
  error('You must compile wavelet toolbox before use')
end


text:
function r = eq(p,q)
% MEAS/EQ  Implement p == q for meas.

% find the difference between the two
diff = meas();
diff = p - q;

% see if p == q at the 1-sigma level
if abs(diff.value) < diff.error
    r = 1;
else
    r = 0;
end


text:
%% subplotPPT _ Demo script
% Demonstrates different ways in which the function subplotPPT can work
%
% Copyright 2008, The MathWorks, Inc. MATLAB and Simulink are registered
% trademarks of The MathWorks, Inc. See www.mathworks.com/trademarks for a
% list of additional trademarks. Other product or brand names may be
% trademarks or registered trademarks of their respective holders.

%% Individual usage
% Here we plot successive windows to the third  slide in a  presentation

close all
clear
h(1) = figure;
peaks;
h(2) = figure;
membrane;
h(3) = figure;
spy;
h(4) = figure;
image;

for ii = 1:numel(h)
    subplotPPT(2,2,ii,h(ii),[pwd,filesep,'test.ppt'],1);
end

%% Call together
% Instead of calling the function in a for loop we can group the figures
% and call the function once
close all
clear
h(1) = figure;
peaks;
h(2) = figure;
membrane;
h(3) = figure;
spy;
h(4) = figure;
image;

subplotPPT(2,2,[4:-1:1],h,[pwd,filesep,'test.ppt'],[]);

%% Try merging plots (Similar to subplot)
close all
clear
h(1) = figure;
peaks;
h(2) = figure;
membrane;
h(3) = figure;
spy;
h(4) = figure;
image;
h(5) = figure;
tori4;

subplotPPT(3,3,[1 4; 2 3; 5 5; 6 9; 7 8],h,[pwd,filesep,'test.ppt'],3,...
    'region',[.1 .1 .8 .8],'hgap',10,'vgap',10);

%% Try different way of specifying region
close all
clear
h(1) = figure;
peaks;
h(2) = figure;
membrane;
h(3) = figure;
spy;
h(4) = figure;
image;
h(5) = figure;
tori4;

subplotPPT(3,3,[1 4; 2 3; 5 5; 6 9; 7 8],h,[pwd,filesep,'TestImage.gif'],4,...
    'region',[200 100 500 340],'hgap',10,'vgap',10,'units','pixels');

%% Create image files
% rather than save to clipboard now create image files in a separate
% directory akin to the puiblish command

subplotPPT(3,3,[1 4; 2 3; 5 5; 6 9; 7 8],h,[pwd,filesep,'test.ppt'],7,...
    'region',[200 100 500 340],'hgap',.1,'vgap',.1,'gapunits','norm','units','pixels',...
    'ImageFormat','emf');

%% Create image files
% rather than save to clipboard now create image files in a separate
% directory akin to the puiblish command
close all
clear
h(1) = figure;
peaks;
h(2) = figure;
membrane;
h(3) = figure;
spy;
h(4) = figure;
image;
h(5) = figure;
tori4;
%%
subplotPPT(3,3,[1 4; 2 3; 5 5; 6 9; 7 8],h,[],5,...
    'region',[200 100 500 340],'hgap',.1,'vgap',.1,'gapunits','norm','units','pixels',...
    'ImageFormat','emf');

%% Create image files
% rather than save to clipboard now create image files in a separate
% directory akin to the puiblish command
close all
clear
h(1) = figure;
peaks;
h(2) = figure;
membrane;
h(3) = figure;
spy;
h(4) = figure;
image;
h(5) = figure;
tori4;

ext = {'.ppt','.bmp','.emf','.gif','.jpg','.png','.tif'};
Idx = 3;

subplotPPT(3,3,[1 4; 2 3; 5 5; 6 9; 7 8],h,['test',ext{Idx}],3,...
    'region',[200 100 500 340],'hgap',.1,'vgap',.1,'gapunits','norm','units','pixels',...
    'ImageFormat','emf');


text:
function spline_test04 ( )

%*****************************************************************************80
%
%% TEST04 tests BASIS_MATRIX_BETA_UNI, BASIS_MATRIX_TMP.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2009
%
%  Author
%
%    John Burkardt
%
  n = 4;
  ndata = 4;
  nsample = 4;
  tdata = [ -1.0E+00, 0.0E+00, 1.0E+00, 2.0E+00 ]';
  ydata = [  4.0E+00, 7.0E+00, 12.0E+00, 19.0E+00 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST04\n' );
  fprintf ( 1, '  BASIS_MATRIX_BETA_UNI sets up the basis matrix\n' );
  fprintf ( 1, '    for the uniform beta spline.\n' );
%
%  First test
%
  beta1 = 1.0E+00;
  beta2 = 0.0E+00;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  BETA1 = %14f\n', beta1 );
  fprintf ( 1, '  BETA2 = %14f\n', beta2 );

  mbasis = basis_matrix_beta_uni ( beta1, beta2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '    TDATA, YDATA\n' );
  fprintf ( 1, '\n' );
  for i = 1 : ndata
    fprintf ( 1, '%12f  %12f\n', tdata(i), ydata(i) );
  end

  left = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, '    T, Spline(T)\n' );
  fprintf ( 1, '\n' );

  for i = 0 : ndata

    if ( i == 0 )
      tlo = tdata(1) - 0.5E+00 * ( tdata(2) - tdata(1) );
      thi = tdata(1);
    elseif ( i < ndata )
      tlo = tdata(i);
      thi = tdata(i+1);
    elseif ( ndata <= i )
      tlo = tdata(ndata);
      thi = tdata(ndata) + 0.5E+00 * ( tdata(ndata) - tdata(ndata-1) );
    end

    if ( i < ndata )
      jhi = nsample - 1;
    else
      jhi = nsample;
    end

    for j = 0 : jhi

      tval = ( ( nsample - j ) * tlo   ...
             + (           j ) * thi ) ...
             /   nsample;

      yval = basis_matrix_tmp ( left, n, mbasis, ndata, tdata, ydata, tval );

      if ( 0 < i & j == 0 )
        mark = '*';
      else
        mark = ' ';
      end

      fprintf ( 1, '%c  %12f  %12f\n', mark, tval, yval );

    end

  end
%
%  Second test
%
  beta1 = 1.0E+00;
  beta2 = 100.0E+00;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  BETA1 = %14f\n', beta1 );
  fprintf ( 1, '  BETA2 = %14f\n', beta2 );

  mbasis = basis_matrix_beta_uni ( beta1, beta2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '    TDATA, YDATA\n' );
  fprintf ( 1, '\n' );
  for i = 1 : ndata
    fprintf ( 1, '%12f  %12f\n', tdata(i), ydata(i) );
  end

  left = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, '    T, Spline(T)\n' );
  fprintf ( 1, '\n' );

  for i = 0 : ndata

    if ( i == 0 )
      tlo = tdata(1) - 0.5E+00 * ( tdata(2) - tdata(1) );
      thi = tdata(1);
    elseif ( i < ndata )
      tlo = tdata(i);
      thi = tdata(i+1);
    elseif ( ndata <= i )
      tlo = tdata(ndata);
      thi = tdata(ndata) + 0.5E+00 * ( tdata(ndata) - tdata(ndata-1) );
    end

    if ( i < ndata )
      jhi = nsample - 1;
    else
      jhi = nsample;
    end

    for j = 0 : jhi

      tval = ( ( nsample - j ) * tlo   ...
             + (           j ) * thi ) ...
             /   nsample;

      yval = basis_matrix_tmp ( left, n, mbasis, ndata, tdata, ydata, tval );

      if ( 0 < i & j == 0 )
        mark = '*';
      else
        mark = ' ';
      end

      fprintf ( 1, '%c  %12f  %12f\n', mark, tval, yval );

    end

  end
%
%  Third test
%
  beta1 = 100.0E+00;
  beta2 = 0.0E+00;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  BETA1 = %14f\n', beta1 );
  fprintf ( 1, '  BETA2 = %14f\n', beta2 );

  mbasis = basis_matrix_beta_uni ( beta1, beta2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '    TDATA, YDATA\n' );
  fprintf ( 1, '\n' );
  for i = 1 : ndata
    fprintf ( 1, '%12f  %12f\n', tdata(i), ydata(i) );
  end

  left = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, '    T, Spline(T)\n' );
  fprintf ( 1, '\n' );

  for i = 0 : ndata

    if ( i == 0 )
      tlo = tdata(1) - 0.5E+00 * ( tdata(2) - tdata(1) );
      thi = tdata(1);
    elseif ( i < ndata )
      tlo = tdata(i);
      thi = tdata(i+1);
    elseif ( ndata <= i )
      tlo = tdata(ndata);
      thi = tdata(ndata) + 0.5E+00 * ( tdata(ndata) - tdata(ndata-1) );
    end

    if ( i < ndata )
      jhi = nsample - 1;
    else
      jhi = nsample;
    end

    for j = 0 : jhi

      tval = ( ( nsample - j ) * tlo   ...
             + (           j ) * thi ) ...
             /   nsample;

      yval = basis_matrix_tmp ( left, n, mbasis, ndata, tdata, ydata, tval );

      if ( 0 < i & j == 0 )
        mark = '*';
      else
        mark = ' ';
      end

      fprintf ( 1, '%c  %12f  %12f\n', mark, tval, yval );

    end

  end

  return
end


text:
clear all; close all; clc

%% folder setup
outputDir = GetOutputDataDir;

%% init

hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);

%% run fish
range_fish = GetFishRange;%[1:3,5:18];%

M_stim = {'5','P','O','L','D','Y'};
M_stimname = {'phT-DF','phT','OMR','Loom','DF','Dot'};

M_pair_range = {[2,1],[2,1],[12,11],[14,15],[0,3],[13,0]}; % [3,2] for PT L vs R; [9,8] for OMR L vs R
%     leftstates_full = [2,12,14,22];
%     rightstates_full = [1,11,15,21];
%     leftstates_full = [2,12,14,22];
%     rightstates_full = [1,11,15,21];
%
% if fishset == 1,
%     States = [0,1,2,3];
%     singleNames = {'black','phototaxis R','phototaxis L','white'};%,...
% %         'L on','R on','L off','R off'};
% % elseif fishset == 2,
% %     States = [0,1,2,3,4,10,11,12];
% %     names = {'black','phototaxis left','phototaxis right','white','grey',...
% %         'OMR forward','OMR left','OMR right',...
% %         'left PT&OMR','right PT&OMR'};
% else%if fishset == 2,
%        States = [0,1,2,3,4,9,10,11,12,13,14,15,21,22,23];
%     singleNames = {'black','phototaxis R','phototaxis L','white','grey',...% 0-4
%         'OMR backward', 'OMR forward','OMR right','OMR left',... % 9-12
%         'Dot','looming L','looming R',... % 13-15
%         'red/blue R','red/blue L','red/red',...
%         };
% end
nStim = length(M_stim);
% for fictive channels (5 channels)
P_pval5= nan(18,5,nStim);
% for motor-seeds (L/R)
P_pval2= nan(18,2,nStim);

ClusterIDs = [1,1];%GetClusterIDs('all');
% M_stimname = {'OMR'};
for i_stim = 1:length(M_stim)
    [M_stimrange,stimrange] = GetStimRange(M_stim{i_stim});

    for i_fish = range_fish
        stimrange = M_stimrange{i_fish};

        tteststimrange = M_pair_range{i_stim};

        if ~isempty(stimrange)
            setappdata(hfig,'isMotorseed',0);
            LoadSingleFishDefault(i_fish,hfig,ClusterIDs,stimrange,0);
            stim = getappdata(hfig,'stim');
            behavior = getappdata(hfig,'behavior');

            %%
            samples = cell(1,length(tteststimrange));
            for i = 1:length(tteststimrange),
                IX = find(stim==tteststimrange(i));
                samples{i} = behavior(:,IX)';
            end

            [~, p] = ttest2(samples{1},samples{2});

            P_pval5(i_fish,:,i_stim) = p;

            %%
            setappdata(hfig,'isMotorseed',1);
            [~,~,behavior] = UpdateTimeIndex(hfig);

            samples = cell(1,length(tteststimrange));
            for i = 1:length(tteststimrange),
                IX = find(stim==tteststimrange(i));
                samples{i} = behavior(:,IX)';
            end

            [~, p] = ttest2(samples{1},samples{2});

            P_pval2(i_fish,:,i_stim) = p;

        end
    end
end

%%
figure('Position',[550,100,800,120]);
for i_stim = 2:6
    subplot(1,5,i_stim-1); hold on;
    y = P_pval2(:,:,i_stim);
    h = histogram(y(:),[0:0.05:1]);
    h.FaceColor = [0.5,0.5,0.5];
    set(gca,'xscale','log');
%     xlim([0,1]);
    ylim([0,30])
    plot([0.05,0.05],[0,30],'r--')
    title(M_stimname{i_stim})
    length(find(y(:)>1))
end

author: xiuyechen

repo: FishExplorer


text:
display('- - - - - - 0 - - - - - -')
display('Closing open figures and clearing variables...')
display('- - - - - - - - - - - - -')
clear all,close all,set(0,'DefaultFigureWindowStyle','docked')

path: github-repos/MATLAB/Sable-mcbench-benchmarks/mcbench-benchmarks-ba13b2f0296ef49491b95e3f984c7c41fccdb6d8/34985-monte-carlo-simulation-of-three-dimensional-grain-growth-code-version-no-1-basic/Monte Carlo Simulation Q-state Potts model 3D square-lattice - microstructure/InitialThings_3D_QPOTTS.m