text:
classdef InitialStateModel < matlab.mixin.SetGet
    %InitialStateModel Summary of this class goes here
    %   Detailed explanation goes here

    properties
        orbitModel(1,1) = GeographicElementSet.getDefaultElements();

        lvState LaunchVehicleState
        stageStates LaunchVehicleStageState

        aero(1,1) LaunchVehicleAeroState
        thirdBodyGravity(1,1) LaunchVehicle3BodyGravState

        steeringModels SteeringModelsSet
        throttleModels ThrottleModelsSet

        optVar InitialStateVariable
    end

    %deprecated
    properties(Access=private)
        steeringModel AbstractSteeringModel = RollPitchYawPolySteeringModel.getDefaultSteeringModel();
        throttleModel(1,1) AbstractThrottleModel = ThrottlePolyModel.getDefaultThrottleModel();
    end

    properties(Dependent)
        time(1,1) double
        centralBody(1,1) KSPTOT_BodyInfo
    end

    methods
        function obj = InitialStateModel()
            obj.steeringModels = SteeringModelsSet();
            obj.throttleModels = ThrottleModelsSet();
        end

        function time = get.time(obj)
            time = obj.orbitModel.time;
        end

        function set.time(obj, newTime)
            obj.orbitModel.time = newTime;
        end

        function time = get.centralBody(obj)
            time = obj.orbitModel.frame.getOriginBody();
        end

        function set.centralBody(obj, newCentralBody)
            if(isempty(obj.orbitModel.frame))
                obj.orbitModel.frame = newCentralBody.getBodyCenteredInertialFrame();
            end
            obj.orbitModel.frame.setOriginBody(newCentralBody);
        end

        function set.orbitModel(obj, newOrbitModel)
            if(isa(newOrbitModel,'AbstractOrbitStateModel'))
                if(isa(newOrbitModel,'BodyFixedOrbitStateModel'))
                    newFrame = obj.centralBody.getBodyFixedFrame(); %#ok<MCSUP>
                    elemSet = GeographicElementSet(obj.time, newOrbitModel.lat, newOrbitModel.long, newOrbitModel.alt, ...
                                                   newOrbitModel.vVectNEZ_az, newOrbitModel.vVectNEZ_el, newOrbitModel.vVectNEZ_mag, newFrame); %#ok<MCSUP>

                elseif(isa(newOrbitModel,'KeplerianOrbitStateModel'))
                    newFrame = obj.centralBody.getBodyCenteredInertialFrame(); %#ok<MCSUP>
                    elemSet = KeplerianElementSet(obj.time, newOrbitModel.sma, newOrbitModel.ecc, newOrbitModel.inc, ...
                                                  newOrbitModel.raan, newOrbitModel.arg, newOrbitModel.tru, newFrame); %#ok<MCSUP>

                elseif(isa(newOrbitModel,'CR3BPOrbitStateModel'))
                    error('No conversion available for CR3BPOrbitStateModel to new element set models.');
                end

                obj.orbitModel = elemSet;
            else
                obj.orbitModel = newOrbitModel;
            end
        end

        function cb = getCentralBodyForStateLog(obj)
            cb = obj.orbitModel.frame.getOriginBody();
        end

        function addStageState(obj, newStageState)
            obj.stageStates(end+1) = newStageState;
        end

        function removeStageStateForStage(obj, stage)
            stageStateInd = find([obj.stageStates.stage] == stage,1,'first');

            obj.stageStates(stageStateInd) = [];
        end

        function stateLogEntry = getInitialStateLogEntry(obj)
            stateLogEntry = LaunchVehicleStateLogEntry();

            ut = obj.time;
            stateLogEntry.time = ut;

            iFrame = obj.centralBody.getBodyCenteredInertialFrame();
            cartElemSet = obj.orbitModel.convertToFrame(iFrame).convertToCartesianElementSet();
            stateLogEntry.position = cartElemSet.rVect;
            stateLogEntry.velocity = cartElemSet.vVect;

            stateLogEntry.centralBody = obj.getCentralBodyForStateLog();
            stateLogEntry.lvState = obj.lvState.deepCopy();

            for(i=1:length(obj.stageStates))
                stateLogEntry.stageStates(i) = obj.stageStates(i).deepCopy(true, stateLogEntry.lvState);
            end

            tankStates = stateLogEntry.getAllTankStates();
            for(i=1:length(tankStates))
                tankState = tankStates(i);
                tankState.tankMass = tankState.tank.initialMass;
            end

            stateLogEntry.event = LaunchVehicleEvent.empty(0,1);
            stateLogEntry.aero = obj.aero.deepCopy();
            stateLogEntry.thirdBodyGravity = obj.thirdBodyGravity.copy();

            stopwatches = stateLogEntry.launchVehicle.stopwatches;
            for(i=1:length(stopwatches))
                stateLogEntry.stopwatchStates(end+1) = stopwatches(i).createInitialState();
            end

            extrema = stateLogEntry.launchVehicle.extrema;
            for(i=1:length(extrema))
                stateLogEntry.extremaStates(end+1) = extrema(i).createInitialState();
            end

            calcObjs = stateLogEntry.launchVehicle.calcObjs;
            for(i=1:length(calcObjs))
                stateLogEntry.calcObjStates(end+1) = calcObjs(i).createInitialState();
            end

            obj.steeringModels.selectedModel.setT0(obj.time);
            obj.throttleModels.selectedModel.setT0(obj.time);

            stateLogEntry.steeringModel = obj.steeringModels.selectedModel;
            stateLogEntry.throttleModel = obj.throttleModels.selectedModel;

            [~,sensors] = stateLogEntry.launchVehicle.lvdData.sensors.getListboxStr();
            for(i=1:length(sensors))
                stateLogEntry.sensorStates(end+1) = sensors(i).getInitialState();
            end
        end

        function optVar = getNewOptVar(obj)
            optVar = InitialStateVariable(obj);
        end

        function optVar = getExistingOptVar(obj)
            optVar = obj.optVar;
        end

        function vars = getAllOptVars(obj)
            vars = obj.getExistingOptVar();

            steeringVar = obj.steeringModels.selectedModel.getExistingOptVar();
            if(not(isempty(steeringVar)))
                vars(end+1) = obj.steeringModels.selectedModel.getExistingOptVar();
            end

            throttleVar = obj.throttleModels.selectedModel.getExistingOptVar();
            if(not(isempty(throttleVar)))
                vars(end+1) = throttleVar;
            end
        end

        function clearAllTankStatesAndRegenerate(obj)
            for(i=1:length(obj.stageStates))
                stgState = obj.stageStates(i);

                stgState.tankStates = LaunchVehicleTankState.empty(1,0);

                stage = stgState.stage;
                for(j=1:length(stage.tanks))
                    tank = stage.tanks(j);

                    newTankState = LaunchVehicleTankState(stgState);
                    newTankState.tank = tank;
                    newTankState.tankMass = tank.initialMass;

                    stgState.addTankState(newTankState);
                end
            end
        end

        function clearDuplicateEngineStates(obj)
            for(i=1:length(obj.stageStates))
                stgState = obj.stageStates(i);

                stage = stgState.stage;
                for(j=1:length(stage.engines))
                    engine = stage.engines(j);
                    engineStates = stgState.engineStates;

                    thisEngineStates = engineStates([engineStates.engine] == engine);
                    if(length(thisEngineStates) > 1)
                        thisEngineStateToSave = thisEngineStates(1);

                        engineStates(engineStates == engine) = LaunchVehicleEngineState.empty(1,0);
                        engineStates(end+1) = thisEngineStateToSave; %#ok<AGROW>
                    elseif(isempty(thisEngineStates))
                        newEngineState = LaunchVehicleEngineState(stgState);
                        newEngineState.engine = engine;
                        engineStates(end+1) = newEngineState; %#ok<AGROW>
                    end

%                     notThisEngineStates = engineStates([engineStates.engine] ~= engine);
%                     if(not(isempty(notThisEngineStates)))
%                         engineStates = setdiff(engineStates, notThisEngineStates);
%                     end

                    stgState.engineStates = engineStates;
                end
            end
        end

        function tf = isVarFromInitialState(obj, var)
            tf = obj.optVar == var;

            if(not(isempty(obj)) && not(isempty(obj.optVar)))
                tf = tf || obj.optVar.isVarContainedWithin(var);
            end

            if(not(isempty(obj.steeringModels.selectedModel.getExistingOptVar())))
                tf = tf || obj.steeringModels.selectedModel.getExistingOptVar() == var;
            end

            if(not(isempty(obj.throttleModels.selectedModel.getExistingOptVar())))
                tf = tf || obj.throttleModels.selectedModel.getExistingOptVar() == var;
            end
        end

        function setInitialStateFromStateLogEntry(obj, stateLogEntry)
            lvdData = stateLogEntry.lvdData;
            varSet = lvdData.optimizer.vars;

            %remove variables
            orbitVar = obj.optVar.orbitVar;
            varSet.removeVariable(orbitVar);

            initStateVar = obj.optVar;
            varSet.removeVariable(initStateVar);

            steerVar = obj.steeringModels.selectedModel.getExistingOptVar();
            varSet.removeVariable(steerVar);

            throttleVar = obj.throttleModels.selectedModel.getExistingOptVar();
            varSet.removeVariable(throttleVar);

            %set elements
            obj.orbitModel = stateLogEntry.getCartesianElementSetRepresentation();

            obj.lvState = stateLogEntry.lvState;
            obj.stageStates = stateLogEntry.stageStates;

            tankStates = stateLogEntry.getAllTankStates();
            for(i=1:length(tankStates))
                tankState = tankStates(i);
                tankState.tank.initialMass = tankState.tankMass;
            end

            obj.aero = stateLogEntry.aero;
            obj.thirdBodyGravity = stateLogEntry.thirdBodyGravity;

            oldSteerModelT0 = stateLogEntry.steeringModel.getT0();
            newSteerModelT0 = stateLogEntry.time;
            tOffsetDelta = newSteerModelT0 - oldSteerModelT0;
            obj.steeringModels.selectedModel = stateLogEntry.steeringModel;
            obj.steeringModels.selectedModel.setInitialAttitudeFromState(stateLogEntry, tOffsetDelta);

            oldThrottleModelT0 = stateLogEntry.throttleModel.getT0();
            newThrottleModelT0 = stateLogEntry.time;
            tOffsetDelta = newThrottleModelT0 - oldThrottleModelT0;
            obj.throttleModels.selectedModel = stateLogEntry.throttleModel;
            obj.throttleModels.selectedModel.setInitialThrottleFromState(stateLogEntry, tOffsetDelta);

            %clean up
            obj.clearDuplicateEngineStates();
        end
    end

    methods(Static)
        function stateLogModel = getDefaultInitialStateLogModelForLaunchVehicle(lv, bodyInfo)
            celBodyData = lv.lvdData.celBodyData;
            stateLogModel = InitialStateModel();

            ut = 0;
            bfFrame = bodyInfo.getBodyFixedFrame();
            geoElemSet = GeographicElementSet(ut, 0, 0, 0, 0, 0, 0, bfFrame);
            stateLogModel.orbitModel = geoElemSet;

            lvsState = LaunchVehicleState(lv);
            stateLogModel.lvState = lvsState;
            lvsState.holdDownEnabled = false;

            for(i=1:length(lv.engineTankConns)) %#ok<*NO4LP>
                e2TConnState = EngineToTankConnState(lv.engineTankConns(i));
                e2TConnState.active = true;
                lvsState.e2TConns(end+1) = e2TConnState;
            end

            stageStates = LaunchVehicleStageState.empty(1,0);
            for(i=1:length(lv.stages))
                stage = lv.stages(i);
                stgState = LaunchVehicleStageState(stage);
                stgState.active = true;

                engines = stage.engines;
                for(j=1:length(engines))
                    engine = engines(j);

                    engineState = LaunchVehicleEngineState(stgState);
                    engineState.engine = engine;
                    engineState.active = true;

                    stgState.engineStates(end+1) = engineState;
                end

                tanks = stage.tanks;
                for(j=1:length(tanks))
                    tank = tanks(j);

                    tankState = LaunchVehicleTankState(stgState);
                    tankState.tank = tank;
                    tankState.tankMass = tank.initialMass;

                    stgState.tankStates(end+1) = tankState;
                end

                stageStates(end+1) = stgState; %#ok<AGROW>
            end
            stateLogModel.stageStates = stageStates;

            aeroState = LaunchVehicleAeroState();
            stateLogModel.aero = aeroState;

            grav3Body = LaunchVehicle3BodyGravState();
            grav3Body.celBodyData = celBodyData;
            stateLogModel.thirdBodyGravity = grav3Body;

            rpyModel = RollPitchYawPolySteeringModel.getDefaultSteeringModel();
            stateLogModel.steeringModel = rpyModel;

            throtModel = ThrottlePolyModel.getDefaultThrottleModel();
            stateLogModel.throttleModel = throtModel;
        end

        function obj = loadobj(obj)
            arguments
                obj InitialStateModel
            end

            if(isempty(obj.steeringModels))
                obj.steeringModels = SteeringModelsSet();
                obj.steeringModels.selectedModel = obj.steeringModel;

            elseif(obj.steeringModels.selectedModel ~= obj.steeringModel)
                obj.steeringModels.selectedModel = obj.steeringModel;
            end

            if(isempty(obj.throttleModels))
                obj.throttleModels = ThrottleModelsSet();
                obj.throttleModels.selectedModel = obj.throttleModel;

            elseif(obj.throttleModels.selectedModel ~= obj.throttleModel)
                obj.throttleModels.selectedModel = obj.throttleModel;
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [dynpcm2] = mbar2dynpcm2(mbar)
% Convert pressure from millibars to dyne per centimeter squared.
% Chad Greene 2012
dynpcm2 = mbar*1000.00;

author: Sable

repo: mcbench-benchmarks


text:
function [ iq_signal ] = read_file_iq( filename )
%read_file_tdoa reads a file with IQ data (e.g. from rtl_sdr)

    disp('read_file_iq');

    % daten vom 1. RX
    disp(['IQ read from data file = ' filename]);
    fileID = fopen(filename);
    a = fread(fileID);
    fclose(fileID);

    inphase1 = a(1:2:end) -128;
    quadrature1 = a(2:2:end) -128;
    disp(['successfully read ' int2str(length(inphase1)) ' samples']);

    % complex representation
    iq_signal = inphase1 + 1i.*quadrature1;

end

repo: tdoa-evaluation-rtlsdr


text:
function r = mescatter(f,E)
  N = size(f,1);
  A = 2*(N+E);

  r = zeros(A,size(f,2));

  ridx = -N-E:-N-1;  fidx = N-E:N-1;
  r(mod(ridx,A)+1,:) = -f(mod(fidx,N)+1,:);
  ridx = -N+1:-1;  fidx = N-1:-1:1;
  r(mod(ridx,A)+1,:) = f(mod(fidx,N)+1,:);
  ridx = 0;  fidx = 0;
  r(mod(ridx,A)+1,:) = f(mod(fidx,N)+1,:) * sqrt(2);
  ridx = 1:N-1;  fidx = 1:N-1;
  r(mod(ridx,A)+1,:) = f(mod(fidx,N)+1,:);
  ridx = N+1:N+E-1;  fidx = N-1:-1:N-E+1;
  r(mod(ridx,A)+1,:) = -f(mod(fidx,N)+1,:);
  ridx = 0;

author: thomaskuestner


text:
function [tf, tf_lambda, tf_U] = isequal(A,B)
%ISEQUAL True if each component of two ktensor's is numerically equal.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

tf = false;
tf_lambda = false;
tf_U = false;

if ~isa(B,'ktensor')
    return;
end

tf_lambda = isequal(A.lambda, B.lambda);
if ncomponents(A) == ncomponents(B)
    tf_U = cellfun(@isequal, A.u, B.u);
end
tf = tf_lambda & all(tf_U);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
% --------------------------------------------------------
% ObjectNet3D
% Copyright (c) 2016 CVGL Stanford University
% Licensed under The MIT License [see LICENSE for details]
% Written by Yu Xiang
% --------------------------------------------------------
% display CAD model and anchor points
% example: draw_cad('car', 1)
function draw_cad(cls, index)

opt = globals();

% load mat file
filename = fullfile(opt.root, sprintf('CAD/mat/%s.mat', cls));
object = load(filename);
cads = object.(cls);
cad = cads(index);

% display mesh
trimesh(cad.faces, cad.vertices(:,1), cad.vertices(:,2), cad.vertices(:,3), 'EdgeColor', 'b');
axis equal;
hold on;

% display anchor points
pnames = cad.pnames;
for i = 1:numel(pnames)
    X = cad.(pnames{i});
    if isempty(X) == 0
        plot3(X(1), X(2), X(3), 'ro', 'LineWidth', 5);
    end
end

author: yuxng


text:
function [result, t] = linfactor (arg1, arg2)
%LINFACTOR factorize a matrix, or use the factors to solve Ax=b.
% Uses LU or CHOL to factorize A, or uses a previously computed factorization to
% solve a linear system.  This function automatically selects an LU or Cholesky
% factorization, depending on the matrix.  A better method would be for you to
% select it yourself.  Note that mldivide uses a faster method for detecting
% whether or not A is a candidate for sparse Cholesky factorization (see spsym
% in the CHOLMOD package, for example).
%
% Example:
%   F = linfactor (A) ;     % factorizes A into the object F
%   x = linfactor (F,b) ;   % uses F to solve Ax=b
%   norm (A*x-b)
%
% A second output is the time taken by the method, ignoring the overhead of
% determining which method to use.  This makes for a fairer comparison between
% methods, since normally the user will know if the matrix is supposed to be
% symmetric positive definite or not, and whether or not the matrix is sparse.
% Also, the overhead here is much higher than mldivide or spsym.
%
% This function has its limitations:
%
% (1) determining whether or not the matrix is symmetric via nnz(A-A') is slow.
%     mldivide (and spsym in CHOLMOD) do it much faster.
%
% (2) MATLAB really needs a sparse linsolve.  See cs_lsolve, cs_ltsolve, and
%     cs_usolve in CSparse, for example.
%
% (3) this function really needs to be written as a mexFunction.
%
% (4) the full power of mldivide is not brought to bear.  For example, UMFPACK
%     is not very fast for sparse tridiagonal matrices.  It's about a factor of
%     four slower than a specialized tridiagonal solver as used in mldivide.
%
% (5) permuting a sparse vector or matrix is slower in MATLAB than it should be;
%     a built-in linfactor would reduce this overhead.
%
% (6) mldivide when using UMFPACK uses relaxed partial pivoting and then
%     iterative refinement.  This leads to sparser LU factors, and typically
%     accurate results.  linfactor uses sparse LU without iterative refinement.
%
% The primary purpose of this function is to answer The Perennially Asked
% Question (or The PAQ for short (*)):  "Why not use x=inv(A)*b to solve Ax=b?
% How do I use LU or CHOL to solve Ax=b?"  The full answer is below.  The short
% answer to The PAQ (*) is "PAQ=LU ... ;-) ... never EVER use inv(A) to solve
% Ax=b."
%
% The secondary purpose of this function is to provide a prototype for some of
% the functionality of a true MATLAB built-in linfactor function.
%
% Finally, the third purpose of this function is that you might find it actually
% useful for production use, since its syntax is simpler than factorizing the
% matrix yourself and then using the factors to solve the system.
%
% See also lu, chol, mldivide, linsolve, umfpack, cholmod.
%
% Oh, did I tell you never to use inv(A) to solve Ax=b?
%
% Requires MATLAB 7.3 (R2006b) or later.

% Copyright 2007, Timothy A. Davis, University of Florida
% VERSION 1.1.0, Nov 1, 2007

if (nargin < 1 | nargin > 2 | nargout > 2)          %#ok
    error ('Usage: F=linfactor(A) or x=linfactor(F,b)') ;
end

if (nargin == 1)

    %---------------------------------------------------------------------------
    % F = linfactor (A) ;
    %---------------------------------------------------------------------------

    A = arg1 ;
    [m n] = size (A) ;
    if (m ~= n)
        error ('linfactor: A must be square') ;
    end

    if (issparse (A))

        % try sparse Cholesky (CHOLMOD): L*L' = P*A*P'
        if (nnz (A-A') == 0 & all (diag (A) > 0))   %#ok
            try
                tic ;
                [L, g, PT] = chol (A, 'lower') ;
                t = toc ;
                if (g == 0)
                    result.L = L ;
                    result.LT = L' ;    % takes more memory, but solve is faster
                    result.P = PT' ;    % ditto.  Need a sparse linsolve here...
                    result.PT = PT ;
                    result.kind = 'sparse Cholesky: L*L'' = P*A*P''' ;
                    result.code = 0 ;
                    return
                end
            catch
		% matrix is symmetric, but not positive definite
		% (or we ran out of memory)
            end
        end

        % try sparse LU (UMFPACK, with row scaling): L*U = P*(R\A)*Q
        tic ;
        [L, U, P, Q, R] = lu (A) ;
        t = toc ;
        result.L = L ;
        result.U = U ;
        result.P = P ;
        result.Q = Q ;
        result.R = R ;
        result.kind = 'sparse LU: L*U = P*(R\A)*Q where R is diagonal' ;
        result.code = 1 ;

    else

        % try dense Cholesky (LAPACK): L*L' = A
        if (nnz (A-A') == 0 & all (diag (A) > 0))                           %#ok
            try
                tic ;
                L = chol (A, 'lower') ;
                t = toc ;
                result.L = L ;
                result.kind = 'dense Cholesky: L*L'' = A' ;
                result.code = 2 ;
                return
            catch
		% matrix is symmetric, but not positive definite
		% (or we ran out of memory)
            end
        end

        % try dense LU (LAPACK): L*U = A(p,:)
        tic ;
        [L, U, p] = lu (A, 'vector') ;
        t = toc ;
        result.L = L ;
        result.U = U ;
        result.p = p ;
        result.kind = 'dense LU: L*U = A(p,:)' ;
        result.code = 3 ;

    end

else

    %---------------------------------------------------------------------------
    % x = linfactor (F,b)
    %---------------------------------------------------------------------------

    F = arg1 ;
    b = arg2 ;

    if (F.code == 0)

        % sparse Cholesky: MATLAB could use a sparse linsolve here ...
        tic ;
        result = F.PT * (F.LT \ (F.L \ (F.P * b))) ;
        t = toc ;

    elseif (F.code == 1)

        % sparse LU: MATLAB could use a sparse linsolve here too ...
        tic ;
        result = F.Q * (F.U \ (F.L \ (F.P * (F.R \ b)))) ;
        t = toc ;

    elseif (F.code == 2)

        % dense Cholesky: result = F.L' \ (F.L \ b) ;
        lower.LT = true ;
        upper.LT = true ;
        upper.TRANSA = true ;
        tic ;
        result = linsolve (F.L, linsolve (F.L, b, lower), upper) ;
        t = toc ;

    elseif (F.code == 3)

        % dense LU: result = F.U \ (F.L \ b (F.p,:)) ;
        lower.LT = true ;
        upper.UT = true ;
        tic ;
        result = linsolve (F.U, linsolve (F.L, b (F.p,:), lower), upper) ;
        t = toc ;
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function product_rule ( list_filename, quad_filename )

%*****************************************************************************80
%
%% MAIN is the main program for PRODUCT_RULE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 May 2007
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n');
  fprintf ( 1, 'PRODUCT_RULE\n');
  fprintf ( 1, '  MATLAB version\n');
  fprintf ( 1, '\n');
  fprintf ( 1, '  Create a multidimensional product rule\n');
  fprintf ( 1, '  as a product of distinct 1D integration rules.\n');
%
%  Get the list filename.
%
  if ( nargin < 1 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'PRODUCT_RULE:\n' );
    list_filename = input ( '  Enter the name of the file listing the factors.' );

  end
%
%  Get the product file prefix.
%
  if ( nargin < 2 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'PRODUCT_RULE:\n' );
    quad_filename = input ( '  Enter the product file prefix to use.' );

  end
%
%  Count the items in the list file.
%
  list_num = file_row_count ( list_filename );
%
%  Determine the spatial dimension and number of points in the product.
%
  dim_num = list_num;
  point_num = product_rule_size ( list_filename, list_num );
%
%  Allocate the product items.
%
  x(1:dim_num,1:point_num) = 0.0;
  w(1:point_num) = 1.0;
  r(1:dim_num,1:2) = 0.0;

  list_unit = fopen ( list_filename, 'rt' );

  if ( list_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PRODUCT_RULE - Fatal error!\n' );
    fprintf ( 1, '  Nonzero value of IOS while opening list file.\n' );
    error ( 'PRODUCT_RULE - Fatal error!' );
  end
%
%  Read the factor information and apply it.
%
  for dim = 1 : dim_num

    quad_1d_filename = fgetl ( list_unit );

    quad_x_1d_filename = strcat ( quad_1d_filename, '_x.txt' );
    quad_w_1d_filename = strcat ( quad_1d_filename, '_w.txt' );
    quad_r_1d_filename = strcat ( quad_1d_filename, '_r.txt' );
%
%  Read the X file.
%
    [ dim_num_1d, point_num_1d ] = r8mat_header_read ( quad_x_1d_filename );

    if ( dim_num_1d ~= 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PRODUCT_RULE - Fatal error!\n' );
      fprintf ( 1, '  The 1D quadrature abscissa file should have exactly\n' );
      fprintf ( 1, '  one value on each line.\n' );
      error ( 'PRODUCT_RULE - Fatal error!' );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Number of points in 1D rule = %d\n', point_num_1d );

    x_1d = r8mat_data_read ( quad_x_1d_filename, dim_num_1d, point_num_1d );
%
%  Read the W file.
%
    [ dim_num_1d, point_num_1d2 ] = r8mat_header_read ( quad_w_1d_filename );

    if ( dim_num_1d ~= 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PRODUCT_RULE - Fatal error!\n' );
      fprintf ( 1, '  The 1D quadrature weight file should have exactly\n' );
      fprintf ( 1, '  one value on each line.\n' );
      error ( 'PRODUCT_RULE - Fatal error!' );
    end

    if ( point_num_1d2 ~= point_num_1d )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PRODUCT_RULE - Fatal error!\n' );
      fprintf ( 1, '  The 1D quadrature weight file should have exactly\n' );
      fprintf ( 1, '  the same number of lines as the abscissa file.\n' );
      error ( 'PRODUCT_RULE - Fatal error!' );
    end

    w_1d = r8mat_data_read ( quad_w_1d_filename, dim_num_1d, point_num_1d );
%
%  Read the R file.
%
    [ dim_num_1d, point_num_1d2 ] = r8mat_header_read ( quad_r_1d_filename );

    if ( dim_num_1d ~= 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PRODUCT_RULE - Fatal error!\n' );
      fprintf ( 1, '  The 1D quadrature region file should have exactly\n' );
      fprintf ( 1, '  one value on each line.\n' );
      error ( 'PRODUCT_RULE - Fatal error!' );
    end

    if ( point_num_1d2 ~= 2 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PRODUCT_RULE - Fatal error!\n' );
      fprintf ( 1, '  The 1D quadrature region file should have exactly\n' );
      fprintf ( 1, '  two lines.\n' );
      error ( 'PRODUCT_RULE - Fatal error!' );
    end

    r_1d = r8mat_data_read ( quad_r_1d_filename, 1, 2 );
%
%  Update the X, W, and R of the product rule.
%
    x = r8vec_direct_product ( dim, point_num_1d, x_1d, dim_num, point_num, x );

    w = r8vec_direct_product2 ( dim, point_num_1d, w_1d, dim_num, point_num, w );

    r(dim,1) = r_1d(1);
    r(dim,2) = r_1d(2);

  end

  fclose ( list_unit );
%
%  Write the product rule.
%
  quad_x_filename = strcat ( quad_filename, '_x.txt' );
  quad_w_filename = strcat ( quad_filename, '_w.txt' );
  quad_r_filename = strcat ( quad_filename, '_r.txt' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Creating product quadrature rule X file = "%s".\n', ...
    quad_x_filename );

  r8mat_write ( quad_x_filename, dim_num, point_num, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Creating product quadrature rule W file = "%s".\n', ...
    quad_w_filename );

  r8mat_write ( quad_w_filename, 1, point_num, w );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Creating product quadrature rule R file = "%s".\n', ...
    quad_r_filename );

  r8mat_write ( quad_r_filename, dim_num, 2, r );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PRODUCT_RULE:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function column_num = file_column_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_COLUMN_COUNT counts the columns in the first line of a file.
%
%  Discussion:
%
%    The file is assumed to be a simple text file.
%
%    Most lines of the file are presumed to consist of COLUMN_NUM words,
%    separated by spaces.  There may also be some blank lines, and some
%    comment lines, which have a "#" in column 1.
%
%    The routine tries to find the first non-comment non-blank line and
%    counts the number of words in that line.
%
%    If all lines are blanks or comments, it goes back and tries to analyze
%    a comment line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the file.
%
%    Output, integer COLUMN_NUM, the number of columns in the file.
%
  FALSE = 0;
  TRUE = 1;
%
%  Open the file.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_COLUMN_COUNT - Error!' );
  end
%
%  Read one line, but skip blank lines and comment lines.
%  Use FGETL so we drop the newline character!
%
  got_one = FALSE;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

    elseif ( line(1) == '#' )

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

  if ( got_one == FALSE )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Warning!\n' );
    fprintf ( 1, '  The file does not seem to contain any data.\n' );
    column_num = -1;
    return;
  end

  column_num = s_word_count ( line );

  return
end
function row_num = file_row_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_ROW_COUNT counts the number of row records in a file.
%
%  Discussion:
%
%    Each input line is a "RECORD".
%
%    The records are divided into three groups:
%
%    * BLANK LINES (nothing but blanks)
%    * COMMENT LINES (begin with a '#')
%    * DATA RECORDS (anything else)
%
%    The value returned by the function is the number of data records.
%
%    By the way, if the MATLAB routine FGETS is used, instead of
%    FGETL, then the variable LINE will include line termination
%    characters, which means that a blank line would not actually
%    have zero characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Output, integer ROW_NUM, the number of rows found.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_ROW_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_ROW_COUNT - Error!' );
  end

  blank_num = 0;
  comment_num = 0;
  row_num = 0;

  record_num = 0;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    record_num = record_num + 1;
    record_length = s_len_trim ( line );

    if ( record_length <= 0 )
      blank_num = blank_num + 1;
    elseif ( line(1) == '#' )
      comment_num = comment_num + 1;
    else
      row_num = row_num + 1;
    end

  end

  fclose ( input_unit );

  return
end
function point_num = product_rule_size ( list_filename, list_num )

%*****************************************************************************80
%
%% PRODUCT_RULE_SIZE returns the size of a product rule of distinct factors.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 May 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string LIST_FILENAME, a file containing a list
%    of prefixes defining quadrature rules.
%
%    Input, integer LIST_NUM, the number of prefixes in the file.
%
%    Output, integer POINT_NUM, the number of points in the product rule.
%
  point_num = 1;

  list_unit = fopen ( list_filename, 'rt' );

  if ( list_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PRODUCT_RULE_SIZE - Fatal error!\n' );
    fprintf ( 1, '  Nonzero value of IOS while opening list file.\n' );
    error ( 'PRODUCT_RULE_SIZE - Fatal error!' );
  end

  for list = 1 : list_num

    quad_1d_filename = fgetl ( list_unit );

    quad_x_1d_filename = strcat ( quad_1d_filename, '_x.txt' );

    [ dim_num_1d, point_num_1d ] = r8mat_header_read ( quad_x_1d_filename );

    point_num = point_num * point_num_1d;

  end

  fclose ( list_unit );

  return
end
function table = r8mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% R8MAT_DATA_READ reads data from an R8MAT file.
%
%  Discussion:
%
%    An R8MAT is an array of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns of data.
%
%    Output, real TABLE(M,N), the point coordinates.
%
  table = zeros ( m, n );
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %f' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the file.\n' );
    error ( 'R8MAT_DATA_READ - Error!' );
  end

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = r8mat_header_read ( input_filename )

%*****************************************************************************80
%
%% R8MAT_HEADER_READ reads the header from an R8MAT file.
%
%  Discussion:
%
%    An R8MAT is an array of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
%  For smaller data files, and less precision, try:
%
%     fprintf ( output_unit, '  %14.6f', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %24.16f', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function x = r8vec_direct_product ( factor_index, factor_order, ...
  factor_value, factor_num, point_num, x )

%*****************************************************************************80
%
%% R8VEC_DIRECT_PRODUCT creates a direct product of R8VEC's.
%
%  Discussion:
%
%    To explain what is going on here, suppose we had to construct
%    a multidimensional quadrature rule as the product of K rules
%    for 1D quadrature.
%
%    The product rule will be represented as a list of points and weights.
%
%    The J-th item in the product rule will be associated with
%      item J1 of 1D rule 1,
%      item J2 of 1D rule 2,
%      ...,
%      item JK of 1D rule K.
%
%    In particular,
%      X(J) = ( X(1,J1), X(2,J2), ..., X(K,JK))
%    and
%      W(J) = W(1,J1) * W(2,J2) * ... * W(K,JK)
%
%    So we can construct the quadrature rule if we can properly
%    distribute the information in the 1D quadrature rules.
%
%    This routine carries out that task.
%
%    Another way to do this would be to compute, one by one, the
%    set of all possible indices (J1,J2,...,JK), and then index
%    the appropriate information.  An advantage of the method shown
%    here is that you can process the K-th set of information and
%    then discard it.
%
%  Example:
%
%    Rule 1:
%      Order = 4
%      X(1:4) = ( 1, 2, 3, 4 )
%
%    Rule 2:
%      Order = 3
%      X(1:3) = ( 10, 20, 30 )
%
%    Rule 3:
%      Order = 2
%      X(1:2) = ( 100, 200 )
%
%    Product Rule:
%      Order = 24
%      X(1:24) =
%        ( 1, 10, 100 )
%        ( 2, 10, 100 )
%        ( 3, 10, 100 )
%        ( 4, 10, 100 )
%        ( 1, 20, 100 )
%        ( 2, 20, 100 )
%        ( 3, 20, 100 )
%        ( 4, 20, 100 )
%        ( 1, 30, 100 )
%        ( 2, 30, 100 )
%        ( 3, 30, 100 )
%        ( 4, 30, 100 )
%        ( 1, 10, 200 )
%        ( 2, 10, 200 )
%        ( 3, 10, 200 )
%        ( 4, 10, 200 )
%        ( 1, 20, 200 )
%        ( 2, 20, 200 )
%        ( 3, 20, 200 )
%        ( 4, 20, 200 )
%        ( 1, 30, 200 )
%        ( 2, 30, 200 )
%        ( 3, 30, 200 )
%        ( 4, 30, 200 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer FACTOR_INDEX, the index of the factor being processed.
%    The first factor processed must be factor 1%
%
%    Input, integer FACTOR_ORDER, the order of the factor.
%
%    Input, real FACTOR_VALUE(FACTOR_ORDER), the factor values
%    for factor FACTOR_INDEX.
%
%    Input, integer FACTOR_NUM, the number of factors.
%
%    Input, integer POINT_NUM, the number of elements in the direct product.
%
%    Input, real X(FACTOR_NUM,POINT_NUM), the elements of the
%    direct product, which are built up gradually.
%
%    Output, real X(FACTOR_NUM,POINT_NUM), the elements of the
%    direct product, updated by the latest factor.
%
%  Local Parameters:
%
%    Local, integer START, the first location of a block of values to set.
%
%    Local, integer CONTIG, the number of consecutive values to set.
%
%    Local, integer SKIP, the distance from the current value of START
%    to the next location of a block of values to set.
%
%    Local, integer REP, the number of blocks of values to set.
%
  persistent contig;
  persistent rep;
  persistent skip;

  if ( factor_index == 1 )
    contig = 1;
    skip = 1;
    rep = point_num;
    x(1:factor_num,1:point_num) = 0.0;
  end

  rep = rep / factor_order;
  skip = skip * factor_order;

  for j = 1 : factor_order

    start = 1 + ( j - 1 ) * contig;

    for k = 1 : rep
      x(factor_index,start:start+contig-1) = factor_value(j);
      start = start + skip;
    end

  end

  contig = contig * factor_order;

  return
end
function w = r8vec_direct_product2 ( factor_index, factor_order, ...
  factor_value, factor_num, point_num, w )

%*****************************************************************************80
%
%% R8VEC_DIRECT_PRODUCT2 creates a direct product of R8VEC's.
%
%  Discussion:
%
%    To explain what is going on here, suppose we had to construct
%    a multidimensional quadrature rule as the product of K rules
%    for 1D quadrature.
%
%    The product rule will be represented as a list of points and weights.
%
%    The J-th item in the product rule will be associated with
%      item J1 of 1D rule 1,
%      item J2 of 1D rule 2,
%      ...,
%      item JK of 1D rule K.
%
%    In particular,
%      X(J) = ( X(1,J1), X(2,J2), ..., X(K,JK))
%    and
%      W(J) = W(1,J1) * W(2,J2) * ... * W(K,JK)
%
%    So we can construct the quadrature rule if we can properly
%    distribute the information in the 1D quadrature rules.
%
%    This routine carries out that task for the weights W.
%
%    Another way to do this would be to compute, one by one, the
%    set of all possible indices (J1,J2,...,JK), and then index
%    the appropriate information.  An advantage of the method shown
%    here is that you can process the K-th set of information and
%    then discard it.
%
%  Example:
%
%    Rule 1:
%      Order = 4
%      W(1:4) = ( 2, 3, 5, 7 )
%
%    Rule 2:
%      Order = 3
%      W(1:3) = ( 11, 13, 17 )
%
%    Rule 3:
%      Order = 2
%      W(1:2) = ( 19, 23 )
%
%    Product Rule:
%      Order = 24
%      W(1:24) =
%        ( 2 * 11 * 19 )
%        ( 3 * 11 * 19 )
%        ( 4 * 11 * 19 )
%        ( 7 * 11 * 19 )
%        ( 2 * 13 * 19 )
%        ( 3 * 13 * 19 )
%        ( 5 * 13 * 19 )
%        ( 7 * 13 * 19 )
%        ( 2 * 17 * 19 )
%        ( 3 * 17 * 19 )
%        ( 5 * 17 * 19 )
%        ( 7 * 17 * 19 )
%        ( 2 * 11 * 23 )
%        ( 3 * 11 * 23 )
%        ( 5 * 11 * 23 )
%        ( 7 * 11 * 23 )
%        ( 2 * 13 * 23 )
%        ( 3 * 13 * 23 )
%        ( 5 * 13 * 23 )
%        ( 7 * 13 * 23 )
%        ( 2 * 17 * 23 )
%        ( 3 * 17 * 23 )
%        ( 5 * 17 * 23 )
%        ( 7 * 17 * 23 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer FACTOR_INDEX, the index of the factor being processed.
%    The first factor processed must be factor 1.
%
%    Input, integer FACTOR_ORDER, the order of the factor.
%
%    Input, real FACTOR_VALUE(FACTOR_ORDER), the factor values for
%    factor FACTOR_INDEX.
%
%    Input, integer FACTOR_NUM, the number of factors.
%
%    Input, integer POINT_NUM, the number of elements in the direct product.
%
%    Output, real W(POINT_NUM), the elements of the
%    direct product, updated by the latest factor.
%
%  Local Parameters:
%
%    Local, integer START, the first location of a block of values to set.
%
%    Local, integer CONTIG, the number of consecutive values to set.
%
%    Local, integer SKIP, the distance from the current value of START
%    to the next location of a block of values to set.
%
%    Local, integer REP, the number of blocks of values to set.
%
  persistent contig;
  persistent rep;
  persistent skip;

  if ( factor_index == 1 )
    contig = 1;
    skip = 1;
    rep = point_num;
    w(1:point_num) = 1.0;
  end

  rep = rep / factor_order;
  skip = skip * factor_order;

  for j = 1 : factor_order

    start = 1 + ( j - 1 ) * contig;

    for k = 1 : rep
      w(start:start+contig-1) = w(start:start+contig-1) * factor_value(j);
      start = start + skip;
    end

  end

  contig = contig * factor_order;

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
%% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LEN, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
function word_num = s_word_count ( s )

%*****************************************************************************80
%
%% S_WORD_COUNT counts the number of "words" in a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be examined.
%
%    Output, integer WORD_NUM, the number of "words" in the string.
%    Words are presumed to be separated by one or more blanks.
%
  FALSE = 0;
  TRUE = 1;

  word_num = 0;
  s_length = length ( s );

  if ( s_length <= 0 )
    return;
  end

  blank = TRUE;

  for i = 1 : s_length

    if ( s(i) == ' ' )
      blank = TRUE;
    elseif ( blank == TRUE )
      word_num = word_num + 1;
      blank = FALSE;
    end

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function hFit = dfaddparamfit(hFit, fitname, distname, dsname, fitframe, exclname, useestimated, fixedvals)
%DFADDPARAMFIT Add parametric fit in dfittool

%   $Revision: 1.1.6.10 $  $Date: 2004/01/24 09:35:10 $
%   Copyright 2003-2004 The MathWorks, Inc.

badfit = false;    % badfit=true means fit failed or not attempted
covok = true;      % covariance calculation can be done
if isempty(hFit)
    newfit = true;
    hFit = stats.dffit(fitname, fitframe);
else
    newfit = false;
end
listeners = hFit.listeners;
set(listeners, 'Enabled', 'off');

frameContents = fitframe.getContentPane;
componentsVector = get(frameContents,'Components');
fittingPanel = componentsVector(1);

% Get data set to fit
ds=find(getdsdb,'name',dsname);
hFit.distname = distname;
hFit.dataset  = dsname;
hFit.fittype = 'param';
hFit.dshandle = ds;

% Store some GUI values in fit
hFit.pfixedtext = fixedvals;
hFit.pestimated = useestimated;

% Extract data from this data set
alpha = 0.05;
hExcl = dfgetexclusionrule(exclname);
[x, cens, freq] = getincludeddata(ds,hExcl);

% Get information about the requested distribution
dist = dfgetdistributions(distname);
if length(dist)~=1 || isempty(x)
   if length(dist)~=1
      emsg = 'Bad distribution name.';
   else
      emsg = 'No data remaining after exclusion rule applied.';
   end
   wmsg = '';
   badfit = true;
end
if length(dist)==1
   hFit.enablebounds = dist.hasconfbounds;
end

% Perform the fit
lasterr('');
lastwarn('');
ws = warning('off');
if badfit
   p = [];
else
   try
      nparams = length(dist.pnames);
      if dist.censoring
         censargs = {cens freq};
      else
         if ~isempty(cens) && any(cens)
            error('stats:dfaddparamfit:NoCensoring',...
                  'Censoring not allowed with the %s distribution', distname);
         elseif ~isempty(freq) && any(freq~=1)
            x = expandInput(x,freq);
            freq = [];
         end
         censargs = {};
      end

      % How many output variables will this return?
      if dist.paramvec
         nparamvars = 1;
      else
         nparamvars = nparams;
      end

      fixedparams = cell(0,1);
      if any(~useestimated)
         for j=1:nparams
            if ~useestimated(j)
               txt = deblank(fixedvals{j});
               if isempty(txt)
                  error('stats:dfaddparamfit:BadParam',...
                        'Invalid value for parameter %s', dist.pnames{j});
               end
               num = str2double(txt);
               if ~isfinite(num)
                  error('stats:dfaddparamfit:BadParam',...
                        'Invalid value for parameter %s', dist.pnames{j});
               end
               fixedparams{length(fixedparams)+1} = num;
            end
         end
      end

      % Set up a cell array to receive outputs, then do the fit
      pcell = cell(nparamvars,1);
      [pcell{:}] = feval(dist.fitfunc, x, fixedparams{:}, alpha, censargs{:});

      % Extract results into a single vector
      if dist.paramvec
         p = pcell{1};
      else
         p = [pcell{:}];
      end
   catch
      p = [];
   end
end
warning(ws);

if ~badfit
   if ~isempty(lastwarn)
      wmsg = sprintf('Warning:  %s',lastwarn);
   else
      wmsg = '';
   end
   emsg = lasterr;
   newmsg = '';
   if any(~isfinite(p))
      newmsg = 'Fit produced infinite parameter estimates.';
   elseif numel(p)~=numel(dist.pnames) || ~isnumeric(p)
      newmsg = 'Fit function returned bad parameter values';
   end
   if ~isempty(newmsg)
      badfit = true;
      emsg = combinemsg(emsg,newmsg);
   end

   % Any type of failure so far makes the covariance calculation questionable
   if ~isempty(wmsg) || ~isempty(emsg)
      covok = false;
   end
end

% Try to get a likelihood value
if isempty(p)
   pcov = [];
   nlogl = NaN;
else
   try
      if ~isempty(dist.likefunc)
         if covok
            [nlogl,pcov] = feval(dist.likefunc, p, x, censargs{:});
         else
            nlogl = feval(dist.likefunc, p, x, censargs{:});
            pcov = [];
         end
      else
         pcov = [];
         nlogl = localnlogl(num2cell(p),dist.pdffunc,dist.cdffunc,x,cens,freq);
      end
      newmsg = '';
   catch
      newmsg = lasterr;
   end
   if isempty(newmsg) && (~isnumeric(nlogl) || ~isscalar(nlogl))
      newmsg = 'Result must be a numeric scalar';
   end
   if isnan(nlogl)
      nlogl = NaN;      % explicitly set to real nan to remove imaginary part
   end
   if ~isempty(newmsg);
      pcov = [];
      nlogl = NaN;
      wmsg = combinemsg(wmsg,...
                        sprintf('Error while evaluating likelihood:\n%s',...
                                newmsg));
   end

end

% Get the range over which to show the fit
dffig = dfgetset('dffig');
ax = findall(dffig,'Type','axes','Tag','main');
xlim = get(ax,'XLim');

% Create a fit object using the information we calculated
if badfit
   resultsText = emsg;
else
   try
      hFit = storefitresults(hFit, dist, p, pcov, nlogl, xlim, hExcl, exclname);
      resultsText = getresults(hFit);
   catch
      resultsText = lasterr;
      badfit = true;
   end
end

resultsText = combinemsg(wmsg,resultsText);

% Show results
hFit.resultstext = resultsText;
fittingPanel.setResults(resultsText)

if ~isempty(hFit)
   if ~newfit && ~(hFit.isgood == ~badfit)
   		com.mathworks.toolbox.stats.FitsManager.getFitsManager.fitIsGoodChanged(java(hFit), ~badfit);
   end
   hFit.isgood = ~badfit;
   if newfit
	  hFit.plot = 1;
      % Add to fit array
      connect(hFit,getfitdb,'up');
   end
end

if hFit.plot
   % Determine if bounds can be shown
   if ~dist.hasconfbounds
      hFit.showbounds = false;
   end

   % Update plotted curve
   updateplot(hFit);

   % Update plot limits
   dfswitchyard('dfupdatexlim');
   dfswitchyard('dfupdateylim');
end

set(listeners, 'Enabled', 'on');

if ~newfit
   com.mathworks.toolbox.stats.FitsManager.getFitsManager.fitChanged(...
       java(hFit),fitname,fitname);
end

% Display a more prominent warning outside the results text
if ~badfit && ~isempty(wmsg)
   warndlg(wmsg,'Distribution Fitting Warning','modal');
end

% ----------------------------------------------
function hFit = storefitresults(hFit, dist, p, pcov, nlogl, xlim, hExcl, exclname)
% Update its properties
hFit.distspec = dist;
hFit.params   = p;
hFit.pcov     = pcov;
hFit.pfixed   = false(size(p));
hFit.loglik = -nlogl;
hFit.support = dist.support;
hFit.exclusionrule = hExcl;
hFit.exclusionrulename = exclname;

hFit.xlim = xlim;
setftype(hFit,dfgetset('ftype'));

% ---------------------------------------------
function nlogl = localnlogl(p, pdf, cdf, x, cens, freq)
% Calculate negative log likelihood

% Handle defaults for option inputs
if isempty(cens)
   cens = false(size(x));
else
   cens = (cens == 1);
end
if isempty(freq)
   freq = ones(size(x));
end

% Compute for uncensored observations
nlogl = - sum(freq(~cens) .* log(feval(pdf, x(~cens), p{:})));

% Add component for censored observations
if any(cens)
   nlogl = nlogl - sum(freq(cens) .* log(1-feval(cdf, x(cens), p{:})));
end

% -----------------------------------
function msg = combinemsg(msg,newmsg)
%COMBINEMSG Combine multiple messages into a single message
if isempty(msg)
   msg = newmsg;
elseif ~isempty(newmsg)
   msg = sprintf('%s\n\n%s',msg,newmsg);
end

% -----------------------------------------
function expanded = expandInput(input,freq)
%EXPANDDATA Expand out an input vector using element frequencies.
if ~isequal(size(input),size(freq))
    error('stats:dfaddparamfit:InputSizeMismatch',...
          'Input argument sizes must match.');
end

% Remove points that have zero frequency
t = (freq == 0);
if any(t)
   input(t) = [];
   freq(t) = [];
end

% Expand the remainder
i = cumsum(freq);
j = zeros(1, i(end));
j(i(1:end-1)+1) = 1;
j(1) = 1;
expanded = input(cumsum(j));

author: zouchuhang

repo: LayoutNet


text:
function pass = test_battery(pref)
% A large battery of functions.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e8 * pref.cheb3Prefs.chebfun3eps;

Battery = {@(x,y,z) cos(pi*x.*y.*z),...                                        % 1
    @(x,y,z) cos(2*pi*x.*y.*z),...                                             % 2
    @(x,y,z) cos(3*pi*x.*y.*z),...                                             % 3
    @(x,y,z) cos(4*pi*x.*y.*z),...                                             % 4
    @(x,y,z) cos(5*pi*x.*y.*z),...                                             % 5
    @(x,y,z) cos(6*pi*x.*y.*z),...                                             % 6
    @(x,y,z) cos(7*pi*x.*y.*z),...                                             % 7
    @(x,y,z) sin(pi*x.*y.*z),...                                               % 8
    @(x,y,z) sin(8*pi*x.*(1-x).*y.*(1-y).*z.*(1-z)),...                        % 9
    @(x,y,z) sin(8*pi*x.*(1-x).*y.*(1-y).*z.*(1-z).*(x-y-z).^2),...            % 10
    @(x,y,z) cos(0*pi*(x-y-z).^2),...                                          % 11
    @(x,y,z) cos(pi*(x-y-z).^2),...                                            % 12
    @(x,y,z) cos(2*pi*(x-y-z).^2),...                                          % 13
    @(x,y,z) exp(sin(4*pi./(1+x)).*sin(4*pi./(1+y)).*sin(4*pi./(1+z))),...     % 14
    @(x,y,z) log(1+x.*y.*z),...                                                % 15
    @(x,y,z) cos(pi*x.*sin(pi*y)) + cos(pi*y.*z.*sin(pi*x.*z)),...             % 16
    @(x,y,z) cos(2*pi*x.*sin(pi*y).*cos(pi*z)) + cos(2*pi*y.*sin(pi*x)),...    % 17
    @(x,y,z) (1-x.*y.*z)./(1+x.^2+y.^2),...                                    % 18
    @(x,y,z) cos(pi*x.*y.^2.*z.^3).*cos(pi*z.*y.^2.*x.^3),...                  % 19
    @(x,y,z) cos(2*pi*x.*y.^2.*z.^3).*cos(2*pi*y.*x.^2.*z.^3),...              % 20
    @(x,y,z) cos(3*pi*x.*y.^2.*(1+z)).*cos(3*pi*y.*x.^2.*(1+z)),...            % 21
    @(x,y,z) (x-y)./(3-x.^2+z.^2),...                                          % 22
    @(x,y,z) exp(-y.*x.^2.*z) + exp(-x.*y.^2.*z),...                           % 23
    @(x,y,z) exp((1-x.^2.*z)./(1+y.^2.*z)) + exp((1-y.^2.*z)./(1+x.^2.*z)),... % 24
    @(x,y,z) 10.^(-x.*y.*z),...                                                % 25
    @(x,y,z) 10.^(-10*x.*y.*z),...                                             % 26
    @(x,y,z) sin(x+y+z),...                                                    % 27
    @(x,y,z) exp(x+y+z).*cos(x+y+z),...                                       % 28
    @(x,y,z) cos(pi/4+x+y+z),...                                               % 29
    @(x,y,z) cos(pi/3+5*x+5*y+5*z),...                                        % 30
    @(x,y,z) (x+y+z).^12,...                                                   % 31
    @(x,y,z) exp(x.^2.*y.^2.*z.^2).*cos(x+y+z),...                             % 32
    @(x,y,z) sin(x.*y.*z)+x+y+z,...                                            % 33
    @(x,y,z) cos(100*x)-cos(100*y)-cos(100*z),...                              % 34
    @(x,y,z) x-y+z,...                                                         % 35
    @(x,y,z) -sin(x).*(sin(x.^2/pi).^2 - sin(y).*sin(y.^2/pi).^2) - sin(z).*sin(z.^2/pi).^2,...                  % 36
    @(x,y,z) -x.*sin(sqrt(abs(0.2+x))) - y.*sin(sqrt(abs(0.3+y))) - z.*sin(sqrt(abs(0.1+z))),...                % 37
    @(x,y,z) (x.^2+y.^2+z.^2)/4000 - cos(x).*cos(y/sqrt(2)) .* cos(z/sqrt(3))+1,...                             % 38
    @(x,y,z) exp(1./(1+(((x+1)/2).*((y+1)/2).*(z+1)/2).^2)),...                % 39
    @(x,y,z) 2*exp(exp((x+1)/2)),...                                           % 40
    };

% Exact values for definite integrals computed with MATLAB's symbolic
% toolbox.
exact = [
    0.8461106227350253, ...
    0.6052157118483288, ...
    0.4697661749231274, ...
    0.3886397051915788, ...
    0.3330920107414529, ...
    0.2928314960173565, ...
    0.2619772790103158, ...
    0.3226674912855009, ...
    0.1153949594968427, ...
    0.0463743346003271, ...
    1.0, ...
    0.3854484437166375, ...
    0.2595975105701650, ...
    1.0700418045865364, ...
    0.1103040719136995, ...
    1.1941278093664705, ...
    0.4916086077355209, ...
    0.5741043338997425, ...
    0.9381453489074456, ...
    0.8830048694732021, ...
    0.3411136775418946, ...
    0.0101828389788635, ...
    1.8529116285910792, ...
    4.2715143503929390, ...
    0.7859343211742496, ...
    0.3352208605577725, ...
    0.8793549306454008, ...
    -0.801590008944990, ...
    -0.577703137905825, ...
    -0.008766419956637, ...
    5220.0021978021978, ...
    0.0366931467095498, ...
    1.6224340287967378, ...
    0.0050636564110975, ...
    0.5, ...
    -0.022434416365347, ...
    -1.181086511303566, ...
    0.2694080277004217, ...
    2.3340368233541618, ...
    17.816513659679896, ...
    ];

dom = [0 1 0 1 0 1];
for jj=1:1:length(Battery)
    % Test construction from the function handle:
    g = chebfun3(Battery{jj}, dom);
    pass(5*jj-4) = abs(sum3(g) - exact(jj)) < tol;

    % Test construction from the Chebfun3 object:
    % 1: Without specifying the variable to be separated in Phase 1:
    h = chebfun3(@(x,y,z) g(x,y,z), dom);
    pass(5*jj-3) = abs(sum3(h) - exact(jj)) < tol;

    % 2: Different variables specified to be separated in Phase 1. The
    % goal is to check that constructor has no issues with any choice of
    % variable chosen for the first separation.
    h1 = chebfun3(@(x,y,z) g(x,y,z), dom, 'fiberDim', 1);
    pass(5*jj-2) = abs(sum3(h1) - exact(jj)) < tol;

    h2 = chebfun3(@(x,y,z) g(x,y,z), dom, 'fiberDim', 2);
    pass(5*jj-1) = abs(sum3(h2) - exact(jj)) < tol;

    h3 = chebfun3(@(x,y,z) g(x,y,z), dom, 'fiberDim', 3);
    pass(5*jj) = abs(sum3(h3) - exact(jj)) < tol;

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef prtPreProcLogDiscPostPlsda < prtPreProcClass
    % prtPreProcLogDisc   Histogram equalization processing
    %
    %   LOGDISC = prtPreProcLogDisc creates a logistic discriminant pre
    %   processing object. A prtPreProcLogDisc object processes the input data
    %   so that each feature dimension is scaled between 0 and 1 to best
    %   match the data set class labels.
    %
    %   prtPreProcLogDisc has no user settable properties.
    %
    %   A prtPreProcLogDisc object also inherits all properties and
    %   functions from the prtAction class
    %
    %   Example:
    %
    %   dataSet = prtDataGenUnimodal;     % Load a data set
    %   logDisc = prtPreProcLogDisc;      % Create a pre processing object
    %
    %   logDisc = logDisc.train(dataSet);  % Train
    %   dataSetNew = logDisc.run(dataSet); % Run
    %
    %   % Plot
    %   subplot(2,1,1); plot(dataSet);
    %   title('Original Data');
    %   subplot(2,1,2); plot(dataSetNew);
    %   title('LogDisc Data');
    %
    %   See Also: prtPreProc, prtPreProcPca, prtPreProcPls,
    %   prtPreProcHistEq, prtPreProcZeroMeanColumns, prtPreProcLda,
    %   prtPreProcZeroMeanRows, prtPreProcLogDisc, prtPreProcZmuv,
    %   prtPreProcMinMaxRows

    properties (SetAccess=private)
        name = 'Logistic Discriminant' % 'Logistic Discriminant'
        nameAbbreviation = 'LogDisc' % LogDisc
    end

    properties (SetAccess=private, Hidden = true)
        % General Classifier Properties
        logDiscWeights = [];
        logDiscMeans = [];
    end

    methods
        function Obj = prtPreProcLogDiscPostPlsda(varargin)
            % Allow for string, value pairs
            % There are no user settable options though.
            Obj = prtUtilAssignStringValuePairs(Obj,varargin{:});
        end
    end

    methods (Access = protected, Hidden = true)

        function Obj = trainAction(Obj,DataSet)

			y = DataSet.getTargetsAsBinaryMatrix;
			LogDisc = prtClassLogisticDiscriminant;
			for iFeature = 1:DataSet.nFeatures
				cLogDisc = LogDisc.train(setTargets(DataSet.retainFeatures(iFeature),y(:,iFeature)));
				Obj.logDiscMeans(iFeature) = cLogDisc.w(1);
				Obj.logDiscWeights(iFeature) = cLogDisc.w(2);
			end

        end

        function DataSet = runAction(Obj,DataSet)
			sigmaFn = @(x) 1./(1 + exp(-x));
			for iFeature = 1:length(Obj.logDiscWeights)
                DataSet.X(:,iFeature) = sigmaFn(DataSet.X(:,iFeature)*Obj.logDiscWeights(iFeature) + Obj.logDiscMeans(iFeature));
			end

			if length(Obj.logDiscWeights) > 1
				DataSet.X = bsxfun(@rdivide, DataSet.X, sum(DataSet.X,2));
			else
				DataSet.X = 1-DataSet.X;
			end

        end

		function X = runActionFast(Obj,X)
			for iFeature = 1:length(Obj.logDiscWeights)
				X(:,iFeature) = 1./(1 + exp(- (X(:,iFeature)*Obj.logDiscWeights(iFeature) + Obj.logDiscMeans(iFeature))));
			end
			if size(X,2) > 1
				X = bsxfun(@rdivide, X, sum(X,2));
			else
				X = 1-X;
			end

        end
    end

end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function v = calcSymAxis(S2F,varargin)

v = equispacedS2Grid('resolution',5*degree,'antipodal');

err = nan(size(v));

for i = 1:length(v)
  err(i) = norm(S2F - S2F.symmetrise(v(i)));
end

[~,id] = min(err);

% refine
v2 = equispacedS2Grid('resolution',0.5*degree,'maxTheta',5*degree,'center',v(id));

err = nan(size(v2));
for i = 1:length(v2)
  err(i) = norm(S2F - S2F.symmetrise(v2(i)));
end

[~,id] = min(err);
v = v2(id);

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function lambda = fibonacci2_eigenvalues ( n )

%*****************************************************************************80
%
%% FIBONACCI2_EIGENVALUES returns the eigenvalues of the FIBONACCI2 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 May 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );

  if ( n == 1 )

    lambda(1,1) = 0.0;

  else

    phi = 0.5 * ( 1.0 + sqrt ( 5.0 ) );

    phi(1)          = phi;
    lambda(2:n-1,1) = 1.0;
    lambda(n,1)     = 1.0 - phi;

  end

  return
end


text:
function ye = fL(t,D,nalpha,nbeta,u,h)
alpha = nalpha*pi;
beta = nbeta*pi;
ye = h*cos(beta*u*t) + exp(-D*alpha*alpha*t)*cos(alpha*u*t);

repo: Aero-matlab


text:
function thickness = SC_ThicknessEstimator(blk, page, sccalib, step_du, step_dv)
%% Estimate Water-Equivalent Thickness (2D)
% Reference: Improved scatter correction using adaptive scatter kernel superposition
% Input:
%               Blk(u,v): Total intensity, i.e., I_0
%               BlkAirNorm:
%               Prm(u,v): Primary intensity, i.e., I_p
%               AirNorm: Primary intensity airnorm chamber value
%               sccalib: Scatter Calibration Structure
% Output:
%               thickness(u,v): Estimated object thickness, i.e., tau(x,y) in Reference
% Date: 2021-05-05
% Author: Yi Du (yi.du@hotmail.com)

% mu H2O = 0.02 /mm
muH2O = str2double(sccalib.CalibrationResults.Globals.muH2O.Text);

% unit mm
tmp = blk./page;
tmp(tmp<0)=0.0001;
thickness = log(tmp) /muH2O;

% fill holes by interpolation
thickness(thickness<0) = NaN;

thickness = single(inpaint_nans(double(thickness), 2));

%% Smooth the estimated thickness
% thickness(vv, uu, ntheta)
thickness = SC_SmoothThickness(thickness, sccalib, step_du, step_dv);

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function Gc=leadlagc(G,Wc,Gam_c,Kv,key)
G=tf(G); [Gai,Pha]=bode(G,Wc);
Phi_c=sin((Gam_c-Pha-180)*pi/180);
den=G.den{1}; a=den(length(den):-1:1);
ii=find(abs(a)<=0); num=G.num{1};
G_n=num(length(num));
if length(ii)>0
   if ii(1)>1, a=a(ii(1)+1);
   else, a=a(ii(1)+1); end
else, a=a(1); end;
alpha=sqrt((1-Phi_c)/(1+Phi_c));
Zc=alpha*Wc; Pc=Wc/alpha;
Kc=sqrt((Wc*Wc+Pc*Pc)/(Wc*Wc+Zc*Zc))/Gai;
K1=G_n*Kc*alpha/a;
if nargin==4, key=1;
   if Phi_c<0, key=2;
   else, if K1<Kv, key=3; end, end
end
switch key
   case 1, Gc=tf([1 Zc]*Kc,[1 Pc]);
   case 2
      Kc=1/Gai; K1=G_n*Kc/a;
      Gc=tf([1 0.1*Wc],[1 K1*Gcn(2)/Kv]);
   case 3
      Zc2=Wc*0.1; Pc2=K1*Zc2/Kv;
      Gcn=Kc*conv([1 Zc],[1,Zc2]);
      Gcd=conv([1 Pc],[1,Pc2]);
      Gc=tf(Gcn,Gcd);
end


text:
function [ out ] = MSE( pic1,pic2 )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
e=0;
[m n]=size(pic1);

for i=1:m
    for j=1:n
        e=e+double((pic1(i,j)-pic2(i,j))^2);
    end
end
out=e/(m*n);
end


text:
% minimum distance from point on ellipse of fly to nose of any other fly
function [data,units] = compute_dell2nose(trx,n)

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

for i1 = 1:nflies,
  fly1 = flies(i1);
  % access closestfly to ensure that dell2nose is computed
  trx(fly1).closestfly_ell2nose;
  data{i1} = trx(fly1).dell2nose;
end
units = parseunits('mm');

author: kristinbranson

repo: JAABA


text:
function fiff_write_double_complex_matrix(fid,kind,mat)
%
% fiff_write_double_complex_matrix(fid,kind,mat)
%
% Writes a double-precision complex matrix tag
%
%     fid           An open fif file descriptor
%     kind          The tag kind
%     mat           The data matrix
%

%
%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.1  2006/09/23 14:43:56  msh
%   Added routines for writing complex and double complex matrices.
%   Added routine for writing double-precision real matrix.
%
%

me='MNE:fiff_write_double_complex_matrix';

if nargin ~= 3
        error(me,'Incorrect number of arguments');
end

FIFFT_COMPLEX_DOUBLE = 21;
FIFFT_MATRIX  = bitshift(1,30);
FIFFT_MATRIX_COMPLEX_DOUBLE = bitor(FIFFT_COMPLEX_DOUBLE,FIFFT_MATRIX);
FIFFV_NEXT_SEQ=0;

datasize = 2*8*numel(mat) + 4*3;

count = fwrite(fid,int32(kind),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(FIFFT_MATRIX_COMPLEX_DOUBLE),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(datasize),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(FIFFV_NEXT_SEQ),'int32');
if count ~= 1
   error(me,'write failed');
end
nrow = size(mat,1);
ncol = size(mat,2);
for j = 1:nrow
   for k = 1:ncol
      count = fwrite(fid,real(mat(j,k)),'double');
      if count ~= 1
          error(me,'write failed');
      end
      count = fwrite(fid,imag(mat(j,k)),'double');
      if count ~= 1
          error(me,'write failed');
      end
   end
end
dims(1) = size(mat,2);
dims(2) = size(mat,1);
dims(3) = 2;
count = fwrite(fid,int32(dims),'int32');
if count ~= 3
    error(me,'write failed');
end

return;


text:
% TEST_openLoopDynamics.m
%
% This script performs some basic checks on the equations of motion.
%
% For example, the total energy should be constant to the tolerance of the
% integrator if the applied torque (u) is zero.
%
% If m1 >> m2, then q should behave like a simple pendulum
%

clc; clear;

%%%% Set up the simulation
z0 = [
    0.0;   %horizontal position
    (pi/180)*80;  %pendulum angle (wrt gravity)
    0.3;   %horizontal velocity
    0.5];  %pendulum angular rate

tSpan = [0,1.5];

p.m1 = 1.0;  % (kg) Cart mass
p.m2 = 0.3;  % (kg) pole mass
p.g = 9.81;  % (m/s^2) gravity
p.l = 0.5;   % (m) pendulum (pole) length

%%%% Function Handles
ctrlFun = @(z)( zeros(size(z(1,:))) );  %Passive controller for now
dynFun = @(t,z)( cartPoleDynamics(z, ctrlFun(z), p) );

%%%% Simulate the system!
options = odeset(...
    'RelTol',1e-8, ...
    'AbsTol',1e-8);
sol = ode45(dynFun, tSpan, z0, options);

%%%% Unpack the simulation
t = linspace(tSpan(1), tSpan(2), 200);
z = deval(sol,t);
u = ctrlFun(z);

%%%% Plots:
figure(1); clf;
plotCartPole(t,z,u,p);  %Plots state, control, and energy vs time.

%%%% Draw Trajectory:
[p1,p2] = cartPoleKinematics(z,p);

figure(2); clf;
nFrame = 5;  %Number of frames to draw
drawCartPoleTraj(t,p1,p2,nFrame);

author: MatthewPeterKelly

repo: dscTutorials


text:
function H=calcSpherHessian(x,systemType,useHalfRange,lTx,lRx,M)
%%CALCSPHERHESSIAN Calculate the Hessian matrix (a matrix of second partial
%          derivatives) for a monostatic or bistatic spherical measurement,
%          ignoring atmospheric effects, with respect to 3D Cartesian
%          position.
%
%INPUTS: x The 3X1 position of the target in Cartesian coordinates in the
%          order [x;y;z].
% systemType An optional parameter specifying the axis from which the
%          angles are measured in radians. Possible values are
%          0 (The default if omitted) Azimuth is measured
%            counterclockwise from the x-axis in the x-y plane. Elevation
%            is measured up from the x-y plane (towards the z-axis). This
%            is consistent with common spherical coordinate systems for
%            specifying longitude (azimuth) and geocentric latitude
%            (elevation).
%          1 Azimuth is measured counterclockwise from the z-axis in the
%            z-x plane. Elevation is measured up from the z-x plane
%            (towards the y-axis). This is consistent with some spherical
%            coordinate systems that use the z axis as the boresight
%            direction of the radar.
%          2 This is the same as 0 except instead of being given
%            elevation, one desires the angle away from the z-axis, which
%            is (pi/2-elevation).
%          3 This is the same as 0 except azimuth is measured clockwise
%            from the y-axis in the x-y plane instead of counterclockwise
%            from the x-axis. This coordinate system often arises when
%            given "bearings" in a local East-North-Up coordinate system,
%            where the bearing directions are measured East of North.
% useHalfRange An optional boolean value specifying whether the bistatic
%          (round-trip) range value has been divided by two. This normally
%          comes up when operating in monostatic mode (the most common
%          type of spherical coordinate system), so that the range
%          reported is a one-way range (or just half a bistatic range).
%          The default if this parameter is not provided is false if lTx
%          is provided and true if it is omitted (monostatic).
%      lTx The 3X1 transmitter position in the global coordinate system
%          with [x;y;z] components. If omitted or an empty matrix is
%          passed, then a vector of zeros is used.
%      lRx The 3X1 receiver position in the global coordinate system
%          with [x;y;z] components. If omitted or an empty matrix is
%          passed, then a vector of zeros is used.
%        M A rotation matrix from the global Coordinate system to the
%          orientation of the coordinate system at the receiver. This is
%          only necessary if UV direction components are desired. If
%          omitted, it is assumed to be the identity matrix.
%
%OUTPUTS: H The 3X3X3 Hessian matrix, where H(:,:,1) is the Hessian for the
%           range component, H(:,:,2) is the Hessian for the azimuth
%           component, and H(:,:,3) is the Hessian for the elevation
%           component. The ordering of the derivatives in each matrix is:
%                  [d^2/(dxdx), d^2/(dxdy), d^2/(dxdz);
%                   d^2/(dydx), d^2/(dydy), d^2/(dydz);
%                   d^2/(dzdx), d^2/(dzdy), d^2/(dzdz)];
%          note that each matrix is symmetric (i.e.
%                   d^2/(dydx)=d^2/(dxdy) ).
%
%This function just calls rangeHessian and spherAngHessian.
%
%June 2017 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<6||isempty(M))
   M=eye(3,3);
end

if((nargin<4||isempty(lTx))&&(nargin<3||isempty(useHalfRange)))
    useHalfRange=true;
elseif(nargin<3||isempty(useHalfRange))
    useHalfRange=false;
end

if(nargin<5||isempty(lRx))
    lRx=zeros(3,1);
end

if(nargin<4||isempty(lTx))
    lTx=zeros(3,1);
end

if(nargin<2||isempty(systemType))
    systemType=0;
end

H=zeros(3,3,3);
H(:,:,1)=rangeHessian(x,useHalfRange,lTx,lRx);
H(:,:,2:3)=spherAngHessian(x,systemType,lRx,M);
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function p = approxeq_pot(A, B, tol)

if nargin < 3, tol = 1e-3; end

p = approxeq(A.p, B.p, tol) & approxeq(A.u, B.u, tol);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% DEMSWISSROLLVARGPLVM1 Run variational GPLVM on swiss roll data.

% VARGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'swissRoll';
experimentNo = 1;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = vargplvmOptions('dtcvar');
options.kern = {'rbfard2', 'bias', 'white'};
options.numActive = 100;
%options.initX = 'ppca';
%options.scale2var1 = 1; % scale data to have variance 1
options.tieParam = 'tied';

options.optimiser = 'scg';
latentDim = 3;
d = size(Y, 2);

% demo using the variational inference method for the gplvm model
model = vargplvmCreate(latentDim, d, Y, options);

model = vargplvmParamInit(model, model.m, model.X);
%model.vardist.covars = 10*ones(model.N,model.q) + 0.001*randn(model.N,model.q);

% Optimise the model.
iters = 15; % Default: 1500
display = 1;

model = vargplvmOptimise(model, display, iters);

% Save the results.
modelWriteResult(model, dataSetName, experimentNo);

if exist('printDiagram') & printDiagram
   % order wrt to the inputScales
  mm = vargplvmReduceModel(model,2);
  lvmPrintPlot(mm, lbls, dataSetName, experimentNo);
end

% Load the results and display dynamically.
lvmResultsDynamic(model.type, dataSetName, experimentNo, 'plot3', model.y)

lvmScatterPlotColor(mm, model.y(:, 2)); %%% ?

% See also: plotting with scatter3

author: SheffieldML

repo: vargplvm


text:
%% This file will generate the guess of the iSINDYc left hand side.
% Last Updated: 2019/05/08
% Coded By: K
%
% Note that this file is modified. When the iter = 6, we will add dx6*x6^4
% to the left hand side guess.
function [P_Data,P_sym]=GuessLib(X,dX,iter,u,Highest_Poly_Order,Highest_Trig_Order,Highest_U_Order)
%% First get the size of the X matrix, determin the data length and the number of variables we have.
[Data_Length,Variable_Number]=size(X);
[~,Variable_Number_dX]=size(dX);
[~,Variable_Number_u]=size(u);
% Also create the symbolic variable
Symbol=sym('z',[Variable_Number,1]);
Symbol_dX=sym('dz',[Variable_Number,1]);
Symbol_u=sym('u',[Variable_Number_u,1]);

%% Now according the Highest Polynomial Order entered, we will calculate the data matrix.
Data=[];
Index=1;

%% First calculate the polynomial term

%Order Zero:
Data(:,Index)=ones(Data_Length,1);
Sym_Struct{1,Index}=sym(1);

%Order One:
if Highest_Poly_Order>=1
    for i=1:Variable_Number
        Index=Index+1;
        Data(:,Index)=X(:,i);
        Sym_Struct{1,Index}=Symbol(i,1);
    end
end

%Order Two:
if Highest_Poly_Order>=2
    for i=1:Variable_Number
        for j=i:Variable_Number
            Index=Index+1;
            Data(:,Index)=X(:,i).*X(:,j);
            Sym_Struct{1,Index}=Symbol(i,1)*Symbol(j,1);
        end
    end
end

%Order Three:
if Highest_Poly_Order>=3
    for i=1:Variable_Number
        for j=1:Variable_Number
            for k=1:Variable_Number
                Index=Index+1;
                Data(:,Index)=X(:,i).*X(:,j).*X(:,k);
                Sym_Struct{1,Index}=Symbol(i,1)*Symbol(j,1)*Symbol(k,1);
            end
        end
    end
end

%Order Four:
if Highest_Poly_Order>=4
    for i=1:Variable_Number
        for j=1:Variable_Number
            for k=1:Variable_Number
                for p=1:Variable_Number
                    Index=Index+1;
                    Data(:,Index)=X(:,i).*X(:,j).*X(:,k).*X(:,p);
                    Sym_Struct{1,Index}=Symbol(i,1)*Symbol(j,1)*Symbol(k,1)*Symbol(p,1);
                end
            end
        end
    end
end

%% Then add the Trigonometric Function in the output data:

%Order One:
if Highest_Trig_Order>=1
    for i=1:Variable_Number
        Index=Index+1;
        Data(:,Index)=sin(X(:,i));
        Sym_Struct{1,Index}=sin(Symbol(i,1));
    end
    for i=1:Variable_Number
        Index=Index+1;
        Data(:,Index)=cos(X(:,i));
        Sym_Struct{1,Index}=cos(Symbol(i,1));
    end
end

%Order Two:
if Highest_Trig_Order>=2
    %
    for i=1:Variable_Number
        for j=i+1:Variable_Number
            Index=Index+1;
            Data(:,Index)=sin(X(:,i)-X(:,j)).^2;
            Sym_Struct{1,Index}=sin(Symbol(i,1)-Symbol(j,1))^2;
        end
    end
    for i=1:Variable_Number
        for j=i+1:Variable_Number
            Index=Index+1;
            Data(:,Index)=cos(X(:,i)-X(:,j)).^2;
            Sym_Struct{1,Index}=cos(Symbol(i,1)-Symbol(j,1))^2;
        end
    end
end

%Order Three:
if Highest_Trig_Order>=3
    %
    for i=1:Variable_Number
        for j=1:Variable_Number
            Index=Index+1;
            Data(:,Index)=sin(X(:,i)+X(:,j));
            Sym_Struct{1,Index}=sin(Symbol(i,1)+Symbol(j,1));
        end
    end
    for i=1:Variable_Number
        for j=1:Variable_Number
            Index=Index+1;
            Data(:,Index)=cos(X(:,i)+X(:,j));
            Sym_Struct{1,Index}=cos(Symbol(i,1)+Symbol(j,1));
        end
    end
    %
    for i=1:Variable_Number
        for j=i+1:Variable_Number
            Index=Index+1;
            Data(:,Index)=sin(X(:,i)-X(:,j));
            Sym_Struct{1,Index}=sin(Symbol(i,1)-Symbol(j,1));
        end
    end
    for i=1:Variable_Number
        for j=i+1:Variable_Number
            Index=Index+1;
            Data(:,Index)=cos(X(:,i)-X(:,j));
            Sym_Struct{1,Index}=cos(Symbol(i,1)-Symbol(j,1));
        end
    end
end

%Order Four:
if Highest_Trig_Order>=4
    for i=1:Variable_Number
        for j=i+1:Variable_Number
            Index=Index+1;
            Data(:,Index)=sin(X(:,i)-X(:,j)).*X(:,i).^2;
            Sym_Struct{1,Index}=sin(Symbol(i,1)-Symbol(j,1))*Symbol(i,1)^2;
        end
    end
    %
    for i=1:Variable_Number
        for j=i+1:Variable_Number
            Index=Index+1;
            Data(:,Index)=cos(X(:,i)-X(:,j)).*X(:,i).^2;
            Sym_Struct{1,Index}=cos(Symbol(i,1)-Symbol(j,1))*Symbol(i,1)^2;
        end
    end
    %
    for i=1:Variable_Number
        for j=1:Variable_Number
            Index=Index+1;
            Data(:,Index)=sin(X(:,i)-2*X(:,j));
            Sym_Struct{1,Index}=sin(Symbol(i,1)-2*Symbol(j,1));
        end
    end
    %
    for i=1:Variable_Number
        for j=1:Variable_Number
            Index=Index+1;
            Data(:,Index)=cos(X(:,i)-2*X(:,j));
            Sym_Struct{1,Index}=cos(Symbol(i,1)-2*Symbol(j,1));
        end
    end
    %
    for i=1:Variable_Number
        for j=i:Variable_Number
            Index=Index+1;
            Data(:,Index)=sin(2*X(:,i)-2*X(:,j)).*X(:,i).^2;
            Sym_Struct{1,Index}=sin(2*Symbol(i,1)-2*Symbol(j,1))*Symbol(i,1)^2;
        end
    end
    %
    for i=1:Variable_Number
        for j=i:Variable_Number
            Index=Index+1;
            Data(:,Index)=cos(2*X(:,i)-2*X(:,j)).*X(:,i).^2;
            Sym_Struct{1,Index}=cos(2*Symbol(i,1)-2*Symbol(j,1))*Symbol(i,1)^2;
        end
    end
end

% Add the guess x6^4 manually for the selected states.
if iter==1 || iter==2 || iter==6
    Index=Index+1;
    Data(:,Index)=X(:,6).^4;
    Sym_Struct{1,Index}=Symbol(6,1)^4;
end

pin=Index;
j=0;
%% Frome here, we add the right hand side to our data

%If the Highest_U_Order is zero, we assume that our system does not have
%control input, vice versa.
if Highest_U_Order==0
    for k=1:pin
        j=j+1;
        P_Data(:,j)=dX(:,1).*Data(:,k);
        P_sym{1,j}=Symbol_dX(iter,1)*(Sym_Struct{1,k});
    end
else
    for i=1:Variable_Number_u
        for k=1:pin
            j=j+1;
            P_Data(:,j)=u(:,1).*Data(:,k);
            P_sym{1,j}=Symbol_u(1,1)*(Sym_Struct{1,k});
        end
    end
end

author: dynamicslab

repo: SINDy-PI


text:
% cmp = MS_complexitybs(x,n);
%
% calculate the Lempel-Ziv complexity of the n-symbol stream x.
%
% x\in{0,1,2,...,(n-1)}
%
% cmp is the normalised complexity, that is the number of distinct
% symbol sequences in x, divided by the expected number of distinct
% symbols for a noise sequence.
%
% Algorithm is implemented in MS_complexitybs.c
%
% Michael Small
% michael.small@uwa.edu.au, http://school.maths.uwa.edu.au/~small/
% 24/9/03
% For further details, please see M. Small. Applied Nonlinear Time Series
% Analysis: Applications in Physics, Physiology and Finance. Nonlinear Science
% Series A, vol. 52. World Scientific, 2005. (ISBN 981-256-117-X) and the
% references therein.

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function dist_signed = line_exp_point_dist_signed_2d ( p1, p2, p )

%*****************************************************************************80
%
%% LINE_EXP_POINT_DIST_SIGNED_2D: signed distance ( explicit line, point ) in 2D.
%
%  Discussion:
%
%    The explicit form of a line in 2D is:
%
%      ( P1, P2 ) = ( (X1,Y1), (X2,Y2) ).
%
%    The signed distance has two interesting properties:
%
%    *  The absolute value of the signed distance is the
%        usual (Euclidean) distance.
%
%    *  Points with signed distance 0 lie on the line,
%       points with a negative signed distance lie on one side
%         of the line,
%       points with a positive signed distance lie on the
%         other side of the line.
%
%    Assuming that C is nonnegative, then if a point is a positive
%    distance away from the line, it is on the same side of the
%    line as the point (0,0), and if it is a negative distance
%    from the line, it is on the opposite side from (0,0).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(2,1), P2(2,1), two points on the line.
%
%    Input, real P(2,1), the point whose signed distance is desired.
%
%    Output, real DIST_SIGNED, the signed distance from the
%    point to the line.
%

%
%  If the explicit line degenerates to a point, the computation is easy.
%
  if ( p1(1:2,1) == p2(1:2,1) )

    dist_signed = sqrt ( sum ( ( p1(1:2,1) - p(1:2,1) ).^2 ) );
%
%  Convert the explicit line to the implicit form A * X + B * Y + C = 0.
%  This makes the computation of the signed distance  to (X,Y) easy.
%
  else

    a = p2(2,1) - p1(2,1);
    b = p1(1,1) - p2(1,1);
    c = p2(1,1) * p1(2,1) - p1(1,1) * p2(2,1);

    dist_signed = ( a * p(1,1) + b * p(2,1) + c ) / sqrt ( a * a + b * b );

  end

  return
end


text:
function [xi,w]=thirdOrderSpherShellCubPoints(numDim,rMin)
%%THIRDORDERSPHERSHELLCUBPOINTS Obtain third order cubature points and
%               weight for integrating spherical region surrounding the
%               origin from radius rMin to radius 1, where rMin>0.
%
%INPUTS: numDim An integer specifying the dimensionality of the points to
%               be generated.
%          rMin The minimum radius of the spherical shell. 0<rMin<1.
%
%OUTPUTS: xi A numDim X numCubaturePoints matrix containing the cubature
%            points. (Each "point" is a vector)
%          w A numCubaturePoints X 1 vector of the weights associated with
%            the cubature points.
%
%This implements formula S_n^{shell} 3-1 in [1], pg. 293, 2*numDim points.
%
%REFERENCES:
%[1] A.H. Stroud, Approximate Calculation of Multiple Integrals. Cliffs,
%    NJ: Prentice-Hall, Inc., 1971.
%
%November 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

rho=rMin;
n=numDim;
V=(2*pi^(n/2)/(n*gamma(n/2)))*(1-rho^n);

r=sqrt(n/(n+2)*((1-rho^(n+2))/(1-rho^n)));
B=V/(2*n);

xi=fullSymPerms([r;zeros(n-1,1)]);
w=B*ones(2*n,1);

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
% Test file for conformal.m.

function pass = test_conformal(pref)

if ( nargin < 1 )
    pref = chebfunpref();
end

% Save the current warning state and then clear it
[lastmsg, lastid] = lastwarn();
lastwarn('');

circle = chebfun('exp(pi*1i*x)','trig');
C = real(circle) + .8i*imag(circle);

% compare the two algorithms on a simple problem
ctr = .1+.1i;
[f,finv] = conformal(C,ctr);
[f2,f2inv] = conformal(C,ctr,'poly');
z = .2-.1i;
err = abs(f(z)-f2(z));
errinv = abs(finv(z)-f2inv(z));
pass(1) = (err < 1e-4) && (errinv < 1e-4);

% snowflake
C = chebfun('exp(pi*1i*t)*(1+.2*cos(6*pi*t))','trig');
[f,finv] = conformal(C);
Z = .5*exp(1i*pi*(1:100)'/100);
err = norm(Z- finv(f(Z)));
pass(2) = (err < 1e-4);

% tanh function
ff = @(z) atanh(2*z)/1.2;
ffinv = @(w) tanh(1.2*w)/2;
C = ffinv(circle);
[f,finv] = conformal(C);
z = -.1i;
pass(3) = abs(f(z)-ff(z)) < 1e-4;
w = .2;
pass(4) = abs(finv(w)-ffinv(w)) < 1e-4;

end


text:
%% TMP

addpath('./libqp/matlab/');

%%
% Example: Training two-class SVM classifier with L2-regularized bias term.

% load training and testing data
load('riply_dataset','trn_X','trn_y','tst_X','tst_y');

% ensure that the labels are +1/-1
trn_y(find(trn_y~=1)) = -1;
tst_y(find(tst_y~=1)) = -1;

% input arguments
lambda_bmrm = 1e-2;        % regularization parameter
options.verb = 1;     % will display progress info

nExamples = length(trn_y);
NUM_THREADS = 2;
% TMP
try
   matlabpool('close');
catch ME2
   idSegLast = regexp(ME2.identifier, '(?<=:)\w+$', 'match');
   if (strcmp(idSegLast, 'OpenConnection'))
       fprintf('Cannot close matlabpool. Skipping...\n');
   end;
end;
fprintf('Starting MATLAB Pool...\n');
try
   matlabpool('open', 'local', NUM_THREADS);
catch ME1
   idSegLast = regexp(ME1.identifier, '(?<=:)\w+$', 'match');
   if (strcmp(idSegLast, 'OpenConnection'))
       fprintf('Already connected to matlab pool. Skipping...\n');
   end;
end
%
idx = randperm(nExamples);
data = cell(NUM_THREADS,1);
to = 0;
for i=1:NUM_THREADS
    from=to+1;
    to = round(i*nExamples/NUM_THREADS);
    data{i}.X0 = 1;
    data{i}.X = trn_X(:,idx(from:to));
    data{i}.y = trn_y(idx(from:to));
end

% lambda = lambda_bmrm*NUM_THREADS;
lambda = lambda_bmrm;

% call BMRM solver
%[W,stat]= parbmrm(data,@risk_svm,lambda,options);
[W,stat]= parbmrm(data,@risk_svm_normalized,lambda,options);

if data{1}.X0 == 0
    W0 = 0;
else
    W0 = W(end);
    W = W(1:end-1);
end

% compute training and testing error
ypred = sign(W'*trn_X + W0);
trn_error = sum(ypred(:)~=trn_y(:))/length(trn_y)

ypred = sign(W'*tst_X + W0);
tst_error = sum(ypred(:)~=tst_y(:))/length(tst_y)

% visualize bounary of the trained liner classification rule
figure;
ppatterns(trn_X,trn_y);
pline(W,W0);

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
% PRTBRVDISCRETEHIERARCHY - PRT BRV Discrete hierarchical model structure
%   Has parameters that specify a dirichlet density

classdef prtBrvDiscreteHierarchy

    properties
        lambda
    end
    methods
        function self = prtBrvDiscreteHierarchy(varargin)
            if nargin < 1
                return
            end

            self = defaultParameters(self,varargin{1});
        end

        function self = defaultParameters(self, nDimensions)
            self.lambda = ones(1,nDimensions)/nDimensions;
            %self.lambda = ones(1,nDimensions)*nDimensions;
        end

        function tf = isValid(self)
            tf = ~isempty(self.lambda);
        end
    end
end


text:
function display(SO3F,varargin)
% standard output

if check_option(varargin,'skipHeader')
  disp('  SO3-function component');
else
  displayClass(SO3F,inputname(1),[],'moreInfo',symChar(SO3F),varargin{:});
end

if length(SO3F) > 1, disp(['  size: ' size2str(SO3F)]); end

if SO3F.antipodal, disp('  antipodal: true'); end

disp(['  weight: ' xnum2str(mean(SO3F,'resolution',5*degree))])

disp(' ')

end


text:
function b = r8bto_vxm ( m, l, a, x )

%*****************************************************************************80
%
%% R8BTO_VXM multiplies a vector by a R8BTO matrix.
%
%  Discussion:
%
%    The R8BTO storage format is for a block Toeplitz matrix. The matrix
%    can be regarded as an L by L array of blocks, each of size M by M.
%    The full matrix has order N = M * L.  The L by L matrix is Toeplitz,
%    that is, along its diagonal, the blocks repeat.
%
%    Storage for the matrix consists of the L blocks of the first row,
%    followed by the L-1 blocks of the first column (skipping the first row).
%    These items are stored in the natural way in an (M,M,2*L-1) array.
%
%  Example:
%
%    M = 2, L = 3
%
%    1 2 | 3 4 | 5 6
%    5 5 | 6 6 | 7 7
%    ----+-----+-----
%    7 8 | 1 2 | 3 4
%    8 8 | 5 5 | 6 6
%    ----+-----+-----
%    9 0 | 7 8 | 1 2
%    9 9 | 8 8 | 5 5
%
%    X = (/ 1, 2, 3, 4, 5, 6 /)
%
%    B = (/ 163, 122, 121, 130, 87, 96 /)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the order of the blocks of the matrix A.
%
%    Input, integer L, the number of blocks in a row or column of A.
%
%    Input, real A(M,M,2*L-1), the R8BTO matrix.
%
%    Input, real X(M*L), the vector to be multiplied.
%
%    Output, real B(M*L), the product X * A.
%

%
%  Construct the right hand side by blocks.
%
  for i = 1 : l

    b(1:m,i) = 0.0;

    for j = 1 : i
      b(1:m,i) = b(1:m,i) + a(1:m,1:m,i+1-j)' * x(1:m,j);
    end

    for j = i+1 : l
      b(1:m,i) = b(1:m,i) + a(1:m,1:m,l+j-i)' * x(1:m,j);
    end

  end

  return
end


text:
function a = uplus(a)
%UPLUS        Implements  +a  for intervals
%

% written  10/16/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged, improved performance
%

repo: OpenMAS


text:
function Population = EnvironmentalSelection(Population,R,N)
% The environmental selection of RPEA

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

    %% Calculate the Tchebychev distances
    Distance = TchebychevDistance(Population.objs,R);

    %% Environmental selection
    RemainP = 1 : length(Population);
    RemainR = 1 : size(R,1);
    while length(RemainP) > length(Population)-N
        if isempty(RemainR)
            RemainR = 1 : size(R,1);
        end
        [temp,imin] = min(Distance(RemainP,RemainR),[],1);
        [~,jmin]    = min(temp);
        imin        = imin(jmin);
        RemainP(imin) = [];
        RemainR(jmin) = [];
    end
    Population = Population(setdiff(1:length(Population),RemainP));
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
classdef (Sealed) MimGuiSingleton < handle
    % MimGuiSingleton. The singleton used by all instances of MimGui
    %
    %     You should not use this class within your own code. It is intended to
    %     be used internally within the TD MIM Toolkit.
    %
    %     MimGuiSingleton ensures there is only a single instance for writing
    %     settings to disk, allowing multiple GUI windows to exist simultaneously.
    %
    %     MimGuiSingleton is a singleton. It cannot be created using the
    %     constructor; instead call MimGuiSingleton.GetGuiSingleton;
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Access = private)
        Settings
    end

    methods (Static)
        function gui_singleton = GetGuiSingleton(app_def, reporting)
            persistent GuiSingleton
            if isempty(GuiSingleton) || ~isvalid(GuiSingleton)
                GuiSingleton = MimGuiSingleton(app_def, reporting);
            end
            gui_singleton = GuiSingleton;
        end
    end

    methods
        function settings = GetSettings(obj)
            settings = obj.Settings;
        end
    end

    methods (Access = private)
        function obj = MimGuiSingleton(app_def, reporting)
            obj.Settings = PTKSettings.LoadSettings(app_def, reporting);
        end
    end

end

author: tomdoel

repo: pulmonarytoolkit


text:
function h = max2(f, g, dims)
%MAX2   Maximum value of a CHEBFUN3 in two directions.
%   MAX2(F) returns a 1D CHEBFUN representing the maximum of the CHEBFUN3
%   object F along the y and z directions, i.e,
%                                          MAX2(F) = @(z) max(F( :, :, z)).
%
%   MAX2(F, [], dims) returns a CHEBFUN representing the maximum of F along
%   the dimensions DIMS, where DIMS = [1, 2] means along the x and y
%   directions, etc.
%
%   WARNING: This function is not always accurate to the expected precision.
%
%   For the global maximum use MAX3.
%
% See also CHEBFUN3/MAX and CHEBFUN3/MAX3.

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

% Empty check:
if ( isempty(f) )
    error('CHEBFUN:CHEBFUN3:max2:input', 'CHEBFUN3 is empty');
end

% Default to max2 of one chebfun3:
if ( nargin < 2 )
    g = [];
end

% Default to maximum along the x and y directions:
if ( nargin < 3 )
    dims = [1, 2];
end

% Do not allow max(F, G):
if ( nargin > 1 && ~isempty(g) )
    error('CHEBFUN:CHEBFUN3:max2:twoCHEBFUN3Inputs', ...
        'Unable to maximize two CHEBFUN3 objects.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We have no idea how to achieve this in an efficient way. This
% is an attempt to return a result, but typically it won't be accurate to
% more than 4-5 digits.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dom = f.domain;
n = 129;
if ( all(dims == [1, 2]) || all(dims == [2, 1]) )
    vals = sample(f, n, n, n);
    temp = chebfun3.unfold(vals, [3]);
    h = chebfun(max(temp, [], 2), dom(5:6), 'splitting', 'on');
    h = simplify(h);
elseif ( all(dims == [1, 3]) || all(dims == [3, 1]) )
    vals = sample(f, n, n, n);
    temp = chebfun3.unfold(vals, [2]);
    h = chebfun(max(temp, [], 2), dom(3:4), 'splitting', 'on');
    h = simplify(h);
elseif ( all(dims == [2, 3]) || all(dims == [3, 2]) )
    vals = sample(f, n, n, n);
    temp = chebfun3.unfold(vals, [1]);
    h = chebfun(max(temp, [], 2), dom(1:2), 'splitting', 'on');
    h = simplify(h);
elseif ( dims == 0 )
    error('CHEBFUN:CHEBFUN3:max2:dims', ...
        'Dimension arguments must be two positive integer scalars within indexing range.')
else
   % return the CHEBFUN3. This is analogous to that MAX() command in
   % MATLAB.
   h = f;
end

end