text:
%% DEBUG SIMULATION SETUP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; close all;

% ADD THE PROGRAM PATHS
addpath('environment');
addpath('objects');
addpath('toolboxes');
addpath('scenarios');

% CODE DEBUG COMMANDS %
% profile on
% profile viewer

fprintf('[SETUP]\tInitialising example script.\n');

%% INITIALISE ANY TOOLBOXES
% IntLab();   % Load Intlab
% OMAS_objectDiagnostics()

%% SIMULATION PARAMETERS
[~, userdir]   = system('echo %USERPROFILE%');  % Get desktop path
sim_outputPath = strcat(userdir,'\desktop\openmas-data');
sim_vebosity   = 1;
sim_warningDistance = 2;
sim_maxDuration = 15;
sim_timeStep    = 0.1;                        % Nominal (0.25s)
sim_idleTimeOut = 5*sim_timeStep;

sim_publishFigures = false;
% sim_publishFigures = true;
sim_figureSet = {'all'};
% sim_figureSet = {'events','plan','inputs','isometric','gif'};
% sim_figureSet = {'plan','inputs','isometric','gif'};
% sim_figureSet = {'isometric','gif'};

%% SCENARIO PARAMETERS
sim_agentNumber     = 5;
sim_agentRadius     = 0.5;
sim_agentOrbit      = 5;
sim_agentVelocity   = 2;
sim_adjacencyMatrix = double(~eye(sim_agentNumber));
sim_waypointOrbit   = 10;
sim_waypointRadius  = 0.1;
sim_offsetAngle     = pi/4;
sim_obstacleNumber  = 3;
sim_noiseSigma      = 0.2;
sim_plotScenario    = true;

%% INITIALISE AGENTS
fprintf('[SETUP]\tAssigning agent definitions:\n');
for index = 1:sim_agentNumber
% BASIC CLASSES
%     agentIndex{index} = objectDefinition('radius',sim_agentRadius);
%     agentIndex{index} = agent('radius',sim_agentRadius);
%     agentIndex{index} = agent_test('radius',sim_agentRadius);
%     agentIndex{index} = agent_2D();
%     agentIndex{index} = agent_2D_test('radius',sim_agentRadius);

    agentIndex{index} = agent_example('radius',sim_agentRadius);

% QUADCOPTER DYNAMICS
%     agentIndex{index} = quadcopter_legacy();
%     agentIndex{index} = quadcopter();
%     agentIndex{index} = quadcopter_formation('adjacencyMatrix',sim_adjacencyMatrix);

% ARdrone DYNAMICS
%     agentIndex{index} = ARdrone_prev();
%     agentIndex{index} = ARdrone('radius',sim_agentRadius);
%     agentIndex{index} = ARdrone_LQR();
%     agentIndex{index} = ARdrone_MPC();

%     agentIndex{index} = fixedWing();
%     agentIndex{index} = globalHawk();
%     agentIndex{index} = boeing737();
%     agentIndex{index} = A10();

%     agentIndex{index} = planetoid();
%     agentIndex{index} = earth();
%     agentIndex{index} = moon();
%     agentIndex{index} = ISS();

% FORMATION CONTROL
%     agentIndex{index} = agent_formation('adjacencyMatrix',sim_adjacencyMatrix);
%     agentIndex{index} = agent_formation_boids();
%     agentIndex{index} = agent_formation_VO();
%     agentIndex{index} = agent_formation_RVO();
%     agentIndex{index} = agent_formation_HRVO();
%     agentIndex{index} = agent_2D_formation_VO('adjacencyMatrix',sim_adjacencyMatrix);
%     agentIndex{index} = agent_2D_formation_RVO();
%     agentIndex{index} = agent_2D_formation_HRVO();
%     agentIndex{index} = agent_2D_formation_ORCA();

% VECTOR SHARING;
%     agentIndex{index} = agent_vectorSharing('radius',sim_agentRadius,'detectionRadius',25);
%     agentIndex{index} = agent_2D_vectorSharing('radius',sim_agentRadius);

% INTERVAL AVOIDANCE
%     agentIndex{index} = agent_interval();
%     agentIndex{index} = agent_IA('radius',sim_agentRadius);
%     agentIndex{index} = agent_2D_IA('radius',sim_agentRadius);

% VELOCITY OBSTACLE METHODS
%     agentIndex{index} = agent_VO('radius',sim_agentRadius);
%     agentIndex{index} = agent_RVO('radius',sim_agentRadius);
%     agentIndex{index} = agent_HRVO('radius',sim_agentRadius);
%     agentIndex{index} = agent_2D_VO('radius',sim_agentRadius);
%     agentIndex{index} = agent_2D_RVO('radius',sim_agentRadius);
%     agentIndex{index} = agent_2D_HRVO('radius',sim_agentRadius);
%     agentIndex{index} = agent_2D_ORCA('radius',sim_agentRadius);
%     agentIndex{index} = agent_2D_VO_withComplex('radius',sim_agentRadius);

% OBSTACLES
%     agentIndex{index} = obstacle();
%     agentIndex{index} = obstacle_cuboid();
%     agentIndex{index} = obstacle_spheroid();
end

for index = 1:sim_obstacleNumber
	% OBSTACLES
%     obstacleIndex{index} = obstacle();
%     obstacleIndex{index} = obstacle_cuboid();
%     obstacleIndex{index} = obstacle_spheroid();
end

%% PLACE AGENT OBJECTS IN PRE-DEFINED SCENARIO
% FORMATION CONTROL TESTS
% [ objectIndex ] = GetScenario_corridor('agents',agentIndex,'adjacencyMatrix',sim_adjacencyMatrix,'plot',sim_plotScenario);
% [ objectIndex ] = GetScenario_formation_split('agents',agentIndex,'agentSpacing',4,'adjacencyMatrix',sim_adjacencyMatrix,'plot',sim_plotScenario,'noiseFactor',sim_noiseSigma);
% [ objectIndex ] = GetScenario_formation_fourObstacles('agents',agentIndex,'obstacleRadius',2,'plot',sim_plotScenario);

% % OBSTACLE TESTS
% [ objectIndex ] = GetScenario_fourCuboidObstacles('agents',agentIndex,'plot',sim_plotScenario);
% [ objectIndex ] = GetScenario_obstacleTrack('agents',agentIndex,'obstacles',obstacleIndex,'plot',sim_plotScenario);

% % AGENT TESTS
% [ objectIndex ] = GetScenario_twoLines('agents',agentIndex,'agentVelocity',sim_agentVelocity,'padding',5,'agentSeparation',sim_agentOrbit,'waypointSeparation',sim_waypointOrbit,'agentVelocity',sim_agentVelocity,'plot',sim_plotScenario);
[ objectIndex ] = GetScenario_concentricRing('agents',agentIndex,'agentOrbit',sim_agentOrbit,'agentVelocity',sim_agentVelocity,'waypointOrbit',sim_waypointOrbit,'waypointRadius',sim_waypointRadius,'plot',sim_plotScenario,'noiseFactor',sim_noiseSigma);
% [ objectIndex ] = GetScenario_concentricRing('agents',agentIndex,'agentOrbit',sim_agentOrbit,'agentVelocity',sim_agentVelocity,'waypointOrbit',sim_waypointOrbit,'offsetAngle',pi/2,'plot',sim_plotScenario,'noiseFactor',sim_noiseSigma);
% [ objectIndex ] = GetScenario_concentricSphere('agents',agentIndex,'agentOrbit',sim_agentOrbit,'agentVelocity',sim_agentVelocity,'waypointOrbit',sim_waypointOrbit,'plot',sim_plotScenario);
% [ objectIndex ] = GetScenario_concentricAngle('agents',agentIndex,'agentOrbit',sim_agentOrbit,'agentVelocity',sim_agentVelocity,'waypointOrbit',sim_waypointOrbit,'angle',sim_offsetAngle,'plot',sim_plotScenario);

% RANDOM TESTS
% [ objectIndex ] = GetScenario_random('objects',agentIndex,'plot',sim_plotScenario);
% [ objectIndex ] = GetScenario_randomNormal('objects',agentIndex,'plot',sim_plotScenario);
% [ objectIndex ] = GetScenario_randomUniform('objects',agentIndex,'is3D',false,'plot',sim_plotScenario);

% WAYPOINT TESTS
% [ objectIndex ] = GetScenario_waypointCurve('agents',agentIndex,'agentVelocity',sim_agentVelocity,'plot',sim_plotScenario);
% [ objectIndex ] = GetScenario_waypoint_90Degrees('agents',agentIndex,'agentVelocity',sim_agentVelocity,'plot',sim_plotScenario);

% % RL TESTS
% [ objectIndex ] = GetScenario_earthOrbit('plot',sim_plotScenario);

%% %%%%%% INITIALISE THE SIMULATION WITH THE OBJECT INDEX %%%%%%%%%%%%%%%%%
[DATA,META] = OMAS_initialise('objects',objectIndex,...
                             'duration',sim_maxDuration,...
                                   'dt',sim_timeStep,...
                          'idleTimeOut',sim_idleTimeOut,...
                              'figures',sim_figureSet,...
                      'warningDistance',sim_warningDistance,...
                            'verbosity',sim_vebosity,...
                           'outputPath',sim_outputPath,...
                          'publishMode',sim_publishFigures);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clearvars -except DATA META
load(strcat(META.outputPath,'META.mat'));
load(strcat(META.outputPath,'EVENTS.mat'));
load(strcat(META.outputPath,'OBJECTS.mat'));

%% EXTRA GRAPHS (DEBUG)
% objectDATA = DATA.objectIndex{1,1}.DATA;
% figure()
% grid on;
% plot(META.TIME.timeVector,objectDATA.inputs)
% xlabel('time (s)');
% legend('Location','southeast');

repo: OpenMAS


text:
function [accuracy,predictlabel,elapse] = CSRKDApredict(fea, gnd, model)
% SRKDApredict: Spectral Regression Kernel Discriminant Analysis Prediction
%               SRKDApredict use SRKDA as a classifier. It used the nearest
%               center rule in the SRKDA subspace for classification.
%
%       [predictlabel,accuracy,elapse] = SRKDApredict(fea, gnd, model);
%
%             Input:
%
%               fea     - data matrix. Each row is a data point.
%               gnd     - Label vector of fea.
%             model     - model trained by SRKDAtrain.m
%
%             Output:
%
%            accuracy   - classification accuracy
%         predictlabel  - predict label for fea
%             elapse    - running time.
%
%    Examples:
%
%
% See also SRKDAtrain, KSR, KSR_caller
%
%Reference:
%
%   [1] Deng Cai, Xiaofei He, and Jiawei Han. "Speed Up Kernel Discriminant
%   Analysis", The VLDB Journal, vol. 20, no. 1, pp. 21-33, January, 2011.
%
%   [2] Deng Cai, Xiaofei He and Jiawei Han, "SRDA: An Efficient Algorithm for
%   Large Scale Discriminant Analysis" IEEE Transactions on Knowledge and
%   Data Engineering, vol. 20, no. 1, pp. 1-12, January, 2008.
%
%   [3] V. Sindhwani, P. Niyogi, M. Belkin, "Beyond  the  Point  Cloud:  from
%   Transductive  to  Semi-supervised  Learning", ICML 2005.
%
%   version 2.0 --December/2011
%   version 1.0 --May/2006
%
%   Written by Deng Cai (dengcai AT gmail.com)
%

MAX_MATRIX_SIZE = 8000; % You can change this number based on your memory.

nTrain = size(model.Landmark,1);
nTest = size(fea,1);
nBlock = ceil(MAX_MATRIX_SIZE*MAX_MATRIX_SIZE/nTrain);
Embed_Test = zeros(nTest,size(model.projection,2));
for i = 1:ceil(nTest/nBlock)
    if i == ceil(nTest/nBlock)
        smpIdx = (i-1)*nBlock+1:nTest;
    else
        smpIdx = (i-1)*nBlock+1:i*nBlock;
    end
    KTest= constructKernel(fea(smpIdx,:),model.Landmark,model.options);
    Embed_Test(smpIdx,:) = KTest*model.projection;
    clear KTest;
end

D = EuDist2(Embed_Test,model.ClassCenter,0);
[dump, idx] = min(D,[],2);
predictlabel = model.ClassLabel(idx);

accuracy = 1 - length(find(predictlabel-gnd))/nTest;

author: ZJULearning

repo: MatlabFunc


text:
classdef ZmapFigureFunction < ZmapFunction
    % ZMAPFIGUREFUNCTION is an interactive figure that is reused
    %    - can be refreshed
    %    - has menu items
    properties(Constant,Abstract)
        FigTag % tag for this figure
    end

    properties
        FigureDetails; % cell containing figure parameters (other than tag)
    end

    properties(Dependent)
        hFig % handle to this figure
    end

    methods(Abstract)
        CreateMenu(obj)
        ClearFigure(obj)
    end

    methods
        function fig = get.hFig(obj)
            fig=findobj('Tag',obj.FigTag,'-and','Type','Figure');
        end

        function doIt(obj)
           CreateFigure(obj);
            % do it is redefined here
        end

        function CreateFigure(obj)
            fig = obj.hFig;
            if isempty(obj.hFig)
                fig=figure('Tag',obj.FigTag,obj.figureDetails{:});
                obj.CreateMenu();
                assert(~isempty(fig) && isvalid(fig),'problem creating figure');

                obj.Calculate();
                obj.plot();
            else
                assert(~isempty(fig) && isvalid(fig),'problem creating figure');
                % duplicated in case a refresh/clearplot is necessary first
                obj.Calculate();
                obj.plot();

            end
        end
    end
end

author: CelsoReyes


text:
function sliceThicknessV = deduceSliceWidths(planC)
%function sliceThicknessV = deduceSliceWidths(planC)
%
%This function deduces slice thicknesses when the treatment planning system
%provides only z Values.
%
%
%Created:  30 Apr 03, JOD.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};

zValuesV = [planC{indexS.scan}.scanInfo(:).zValue];

sliceThicknessV = ones(size(zValuesV)) * nan;

for i = 2 : length(zValuesV) - 1

  nextDelta = abs(zValuesV(i+1) - zValuesV(i));

  sliceThicknessV(i) = nextDelta;

end

sliceThicknessV(1) = 2 * (abs(zValuesV(2) - zValuesV(1)) - 0.5 * sliceThicknessV(2));

sliceThicknessV(end) = 2 * (abs(zValuesV(end) - zValuesV(end - 1)) - 0.5 * sliceThicknessV(end - 1));

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ i, j ] = i4_to_triangle ( k )

%*****************************************************************************80
%
%% I4_TO_TRIANGLE converts an integer to triangular coordinates.
%
%  Discussion:
%
%    Triangular coordinates are handy when storing a naturally triangular
%    array (such as the lower half of a matrix) in a linear array.
%
%    Thus, for example, we might consider storing
%
%    (1,1)
%    (2,1) (2,2)
%    (3,1) (3,2) (3,3)
%    (4,1) (4,2) (4,3) (4,4)
%
%    as the linear array
%
%    (1,1) (2,1) (2,2) (3,1) (3,2) (3,3) (4,1) (4,2) (4,3) (4,4)
%
%    Here, the quantities in parenthesis represent the natural row and
%    column indices of a single number when stored in a rectangular array.
%
%    In this routine, we are given the location K of an item in the
%    linear array, and wish to determine the row I and column J
%    of the item when stored in the triangular array.
%
%  First Values:
%
%     K  I  J
%
%     0  0  0
%     1  1  1
%     2  2  1
%     3  2  2
%     4  3  1
%     5  3  2
%     6  3  3
%     7  4  1
%     8  4  2
%     9  4  3
%    10  4  4
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer K, the linear index of the (I,J) element, which
%    must be nonnegative.
%
%    Output, integer I, J, the row and column indices.
%
  if ( k < 0 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4_TO_TRIANGLE - Fatal error!\n' );
    fprintf ( 1, '  K < 0.\n' );
    fprintf ( 1, '  K = %d\n', k );
    error ( 'I4_TO_TRIANGLE - Fatal error!' );

  elseif ( k == 0 )

    i = 0;
    j = 0;
    return

  end

  i = floor ( sqrt ( 2 * k ) );

  if ( i * i + i < 2 * k )
    i = i + 1;
  end

  j = k - ( i * ( i - 1 ) ) / 2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function fx = p11_fun ( n, x )

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

  return
end


text:
function value = index3_row ( i_min, i, i_max, j_min, j, j_max, ...
  k_min, k, k_max, index_min )

%*****************************************************************************80
%
%% INDEX3_ROW indexes a 3D array by rows.
%
%  Discussion:
%
%    When we say "by rows", we really just mean that entries of the array are
%    indexed starting at entry (I_MIN,J_MIN,K_MIN), and the increasing the LAST
%    index first, then the next-to-the-last, and so on.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 April 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I_MIN, I, I_MAX, for row indices,
%    the minimum, the index, and the maximum.
%
%    Input, integer J_MIN, J, J_MAX, for column indices,
%    the minimum, the index, and the maximum.
%
%    Input, integer K_MIN, K, K_MAX, for plane indices,
%    the minimum, the index, and the maximum.
%
%    Input, integer INDEX_MIN, the index of (I_MIN,J_MIN,K_MIN).
%    Typically, this is 0 or 1.
%
%    Output, integer VALUE, the index of element (I,J,K).
%
  value = index_min ...
             + ( k - k_min ) ...
             + ( j - j_min ) * ( k_max + 1 - k_min ) ...
             + ( i - i_min ) * ( j_max + 1 - j_min ) * ( k_max + 1 - k_min );

  return
end


text:
% KERNEL = binomialFilter(size)
%
% Returns a vector of binomial coefficients of order (size-1) .

% Eero Simoncelli, 2/97.

function [kernel] = binomialFilter(sz)

if (sz < 2)
  error('size argument must be larger than 1');
end

kernel = [0.5 0.5]';

for n=1:sz-2
  kernel = conv([0.5 0.5]', kernel);
end

repo: SelfExSR


text:
classdef C45Classifier < Algorithm

    properties (Access = public)
        shouldTrain = false;
        inc_node;
        region;
    end

    properties (Access = private)
        classifier;
    end

    methods (Access = public)

        function obj = C45Classifier()
            obj.name = 'Ensemble';
            obj.inputPort = DataType.kTable;
            obj.outputPort = DataType.kTable;
        end

        function dataOut = compute(obj,data)
            if obj.shouldTrain
                obj.train(data);
                dataOut = [];
            else
                dataOut = obj.test(data);
            end
        end

        function D = train(obj,table)
            train_features = table.features;
            train_targets = table.labels;

            [Ni, M]		= size(train_features);
            obj.inc_node    = obj.inc_node*M/100;
            Nu          = 10;

            %For the decision obj.region
            N           = obj.region(5);
            mx          = ones(N,1) * linspace (obj.region(1),obj.region(2),N);
            my          = linspace (obj.region(3),obj.region(4),N)' * ones(1,N);
            flatxy      = [mx(:), my(:)]';

            %Find which of the input features are discrete, and discretisize the corresponding
            %dimension on the decision obj.region
            discrete_dim = zeros(1,Ni);
            for i = 1:Ni
                Nb = length(unique(train_features(i,:)));
                if (Nb <= Nu)
                    %This is a discrete feature
                    discrete_dim(i)	= Nb;
                    [~, flatxy(i,:)]	= high_histogram(flatxy(i,:), Nb);
                end
            end

            %Build the tree recursively
            disp('Building tree')
            obj.classifier        = make_tree(train_features, train_targets, obj.inc_node, discrete_dim, max(discrete_dim), 0);
            %{
            %Make the decision obj.region according to the tree
            disp('Building decision surface using the tree')
            targets		= use_tree(flatxy, 1:N^2, obj.classifier, discrete_dim, unique(train_targets));

            D   		= reshape(targets,N,N);
            %}
        end

        function labels = test(obj,table)
            labels = use_tree(table.features,1:N^2, obj.classifier, discrete_dim, unique(train_targets));
        end

        function metrics = computeMetrics(obj,table)
            flops = timeit(@()obj.test(table)) / Constants.kReferenceComputingTime;
            memory = Helper.ComputeObjectSize(obj.classifier);
            outputSize = table.height * Constants.kClassificationResultBytes;
            metrics = Metric(flops,memory,outputSize);
        end

        function str = toString(obj)
            str = sprintf('%s_%d',obj.name,obj.ensembleMethod);
        end

    end

    methods (Access = private)
        function targets = use_tree(features, indices, tree, discrete_dim, Uc)

            targets = zeros(1, size(features,2));

            if (tree.dim == 0)
                %Reached the end of the tree
                targets(indices) = tree.child;
                return;
            end

            %This is not the last level of the tree, so:
            %First, find the dimension we are to work on
            dim = tree.dim;
            dims= 1:size(features,1);

            %And classify according to it
            if (discrete_dim(dim) == 0)
                %Continuous feature
                in				= indices(find(features(dim, indices) <= tree.split_loc));
                targets		= targets + use_tree(features(dims, :), in, tree.child(1), discrete_dim(dims), Uc);
                in				= indices(find(features(dim, indices) >  tree.split_loc));
                targets		= targets + use_tree(features(dims, :), in, tree.child(2), discrete_dim(dims), Uc);
            else
                %Discrete feature
                Uf				= unique(features(dim,:));
                for i = 1:length(Uf)
                    in   	   = indices(find(features(dim, indices) == Uf(i)));
                    targets	= targets + use_tree(features(dims, :), in, tree.child(i), discrete_dim(dims), Uc);
                end
            end
        end

        function tree = make_tree(features, targets, discrete_dim, maxNbin, base)

            %Build a tree recursively

            [Ni, L]    					= size(features);
            Uc         					= unique(targets);
            tree.dim						= 0;
            %tree.child(1:maxNbin)	= zeros(1,maxNbin);
            tree.split_loc				= inf;

            if isempty(features)
                return;
            end

            %When to stop: If the dimension is one or the number of examples is small
            if ((obj.inc_node > L) | (L == 1) | (length(Uc) == 1))
                H					= hist(targets, length(Uc));
                [~, largest] 	= max(H);
                tree.child	 	= Uc(largest);
                return;
            end

            %Compute the node's I
            for i = 1:length(Uc)
                Pnode(i) = length(find(targets == Uc(i))) / L;
            end
            Inode = -sum(Pnode.*log(Pnode)/log(2));

            %For each dimension, compute the gain ratio impurity
            %This is done separately for discrete and continuous features
            delta_Ib    = zeros(1, Ni);
            split_loc	= ones(1, Ni)*inf;

            for i = 1:Ni
                data	= features(i,:);
                Nbins	= length(unique(data));
                if (discrete_dim(i))
                    %This is a discrete feature
                    P	= zeros(length(Uc), Nbins);
                    for j = 1:length(Uc)
                        for k = 1:Nbins
                            indices 	= find((targets == Uc(j)) & (features(i,:) == k));
                            P(j,k) 	= length(indices);
                        end
                    end
                    Pk          = sum(P);
                    P           = P/L;
                    Pk          = Pk/sum(Pk);
                    info        = sum(-P.*log(eps+P)/log(2));
                    delta_Ib(i) = (Inode-sum(Pk.*info))/-sum(Pk.*log(eps+Pk)/log(2));
                else
                    %This is a continuous feature
                    P	= zeros(length(Uc), 2);

                    %Sort the features
                    [sorted_data, indices] = sort(data);
                    sorted_targets = targets(indices);

                    %Calculate the information for each possible split
                    I	= zeros(1, L-1);
                    for j = 1:L-1
                        for k =1:length(Uc)
                            P(k,1) = length(find(sorted_targets(1:j) 		== Uc(k)));
                            P(k,2) = length(find(sorted_targets(j+1:end) == Uc(k)));
                        end
                        Ps		= sum(P)/L;
                        P		= P/L;
                        info	= sum(-P.*log(eps+P)/log(2));
                        I(j)	= Inode - sum(info.*Ps);
                    end
                    [delta_Ib(i), s] = max(I);
                    split_loc(i) = sorted_data(s);
                end
            end

            %Find the dimension minimizing delta_Ib
            [~, dim] = max(delta_Ib);
            dims		= 1:Ni;
            tree.dim = dim;

            %Split along the 'dim' dimension
            Nf		= unique(features(dim,:));
            Nbins	= length(Nf);
            if (discrete_dim(dim))
                %Discrete feature
                for i = 1:Nbins
                    indices    		= find(features(dim, :) == Nf(i));
                    tree.child(i)	= make_tree(features(dims, indices), targets(indices), discrete_dim(dims), maxNbin, base);
                end
            else
                %Continuous feature
                tree.split_loc		= split_loc(dim);
                indices1		   	= find(features(dim,:) <= split_loc(dim));
                indices2	   		= find(features(dim,:) > split_loc(dim));
                tree.child(1)		= make_tree(features(dims, indices1), targets(indices1), obj.inc_node, discrete_dim(dims), maxNbin);
                tree.child(2)		= make_tree(features(dims, indices2), targets(indices2), obj.inc_node, discrete_dim(dims), maxNbin);
            end

        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function A = sledges2adjmat(n, nt, edges, varargin)
%SLEDGES2ADJMAT Creates an adjacency matrix from edge set
%
% $ Syntax $
%   - A = sledges2adjmat(n, nt, edges, ...)
%
% $ Arguments $
%   - n:            The number of (source) nodes
%   - nt:           The number of (target) nodes
%   - edges:        The matrix of edge set
%
% $ Description $
%   - A = sledges2adjmat(n, nt, edges) creates an adjacency matrix
%     from the edge set. You can specify the following properties:
%       - 'valtype':        the value type of the target matrix
%                           - 'auto': if has value, make numeric matrix
%                                     if no value, make logical matrix
%                           - 'logical': make logical matrix always
%                           - 'numeric': make numeric matrix always
%                           (default = 'auto')
%       - 'sparse':         whether to create a sparse matrix
%                           (default = true)
%       - 'preprune':       whether to prune the edges first
%                           (default = false)
%       - 'prunemethod':    the method used to prune the edge set
%                           (default = [], means using default method)
%                           refer to slpruneedgeset for the specification
%                           of the prune methods.
%       - 'sym':            whether to create symmetric graph
%       - 'symmethod':      the method to symmetrize the graph
%                           (default = [], means using default method)
%                           refer to slsymedgeset for the specification.
%
% $ Remarks $
%   - The property sym can only be true when n == nt.
%
%   - It is an integrated wrapper for slmakeadjmat, slsymgraph and
%     slpruneedgeset.
%
% $ History $
%   - Created by Dahua Lin, on Sep 9, 2006
%

%% parse and verify input arguments

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

if ~isempty(edges)
    ncols = size(edges, 2);
    if ndims(edges) ~= 2 || (ncols ~= 2 && ncols ~= 3)
        error('sltoolbox:invalidarg', ...
            'The edges should be a 2D matrix with two or three columns');
    end
end

opts.valtype = 'auto';
opts.sparse = true;
opts.preprune = false;
opts.prunemethod = [];
opts.sym = false;
opts.symmethod = [];
opts = slparseprops(opts, varargin{:});

if opts.sym
    if n ~= nt
        error('sltoolbox:rterror', ...
            'The sym can only be true when n == nt');
    end
end

switch opts.valtype
    case 'auto'
        islogic = (ncols == 2);
    case 'numeric'
        islogic = false;
    case 'logical'
        islogic = true;
    otherwise
        error('sltoolbox:invalidarg', ...
        'Invalid value type of adjacency matrix: %s', opts.valtype);
end

%% main skeleton

% prune
if opts.preprune
    edges = slpruneedgeset(n, nt, edges, opts.prunemethod);
end

% make adjmat

A = slmakeadjmat(n, nt, edges, [], islogic, opts.sparse);

% symmetrize
if opts.sym
    A = slsymgraph(A, opts.symmethod);
end

author: lmthang


text:
function World = doMapping(World, sen, r_new, scan_data, handles)

    %%%%%%%%%%%%%%%%%%%%%%%%% MAPPING SETTINGS %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Amount by which to increment/reduce cell value when obstacle is detected:
    map_augment = 20;
    % Amount by which to reduce cell value when no obstacle is detected:
    map_reduce = 20;
    map_max = 255;          % Maximum value that a grid cell can take
    map_min = 0;            % Minimum value that a grid cell can take
    map_isOccupied = 150;   % Value at which cell is in occupied state
    map_isEmpty = 50;       % Value at which cell is in empty state
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    scan_polar = getMeasurement(scan_data);
    % Determine the region in the scan that was detected as being emtpy
    empty_region_pol = ones(2, size(scan_polar, 2) * ...
        (round((sen.range - sen.range_min)/World.map_res) + 1)) * -1;
    position = 1;
    for i = 1:size(scan_polar, 2)
        empty_region_r = sen.range_min : World.map_res : scan_polar(1, i) ...
            - World.map_res;
        empty_region_a = repmat(scan_polar(2, i), 1, length(empty_region_r));
        empty_region_pol(:, position : position + length(empty_region_r) - 1) ...
            = [empty_region_r; empty_region_a];
        position = position + length(empty_region_r);
    end
    empty_region_pol = empty_region_pol(:, empty_region_pol(1, :) ~= -1);
    % Map the empty region of the scan onto the global frame
    empty_region = invScanPoint(r_new, empty_region_pol);

    scan_global = transToGlobal(r_new, scan_data);
    gridmap_occupied = zeros(size(World.gridmap));
    gridmap_empty = zeros(size(World.gridmap));
    scan_global = interp1(World.map_vals, World.map_vals, scan_global, 'nearest');
    empty_region = interp1(World.map_vals, World.map_vals, empty_region, 'nearest');

    for i = 1:length(World.map_vals)
        for j = 1:length(World.map_vals)
            if ~isempty(scan_global)
                gridmap_occupied(i, j) = sum(scan_global(1,:) == World.map_vals(j) & ...
                    scan_global(2,:) == World.map_vals(i));
            end
            gridmap_empty(i, j) = sum(empty_region(1,:) == World.map_vals(j) & ...
                empty_region(2,:) == World.map_vals(i));
        end
    end

    World.gridmap_counter = World.gridmap_counter ...
        + map_augment * gridmap_occupied ...
        - map_reduce * gridmap_empty;

    % Set the gridmap counter such that the max value is 255 and min value
    % is 0:

    World.gridmap_counter(World.gridmap_counter > map_max) = map_max;
    World.gridmap_counter(World.gridmap_counter < map_min) = map_min;

    gridmap_occupied = World.gridmap_counter > map_isOccupied;
    gridmap_empty = World.gridmap_counter < map_isEmpty;
    World.gridmap = ones(size(World.gridmap_counter)) * 0.5;
    World.gridmap(gridmap_occupied) = 1;
    World.gridmap(gridmap_empty) = 0;

    World.gridmap_greyscale = mat2gray(World.gridmap_counter, [map_min map_max]);

    imagesc(flipud(World.gridmap_greyscale), 'parent', handles.map1, [0 1]);
    set(handles.map1, 'XTick', [], 'YTick', []);
    imagesc(flipud(World.gridmap), 'parent', handles.map2, [0 1]);
    set(handles.map2, 'XTick', [], 'YTick', []);

end

author: jaijuneja

repo: ekf-slam-matlab


text:
% Tests for chebfun plotting functions.
function pass = test_plot_xylim(pref)

tol = 1e-4;
% Create a figure, and make it invisible. Need to do this a number of time
% throughout the test.
hfig = figure;
set(hfig,'visible','off')
%% Finite functions on unbounded domains

dom1 = [0 pi];
pass1 = [];
x = chebfun(@(x) x, dom1);
plot(x)
pass1(length(pass1) + 1) = ( norm(dom1 - get(gca,'xlim')) < tol);

plot(0.62*sin(x))
yl = get(gca, 'ylim');
pass1(length(pass1) + 1) = ( norm(0.62 - yl(2)) > 0.05);

plot(0.62*[sin(x) 0*x -sin(x)])
yl = get(gca, 'ylim');
pass1(length(pass1) + 1) = ( norm(0.62*[-1 1] - yl) > 0.05);
pass1(length(pass1) + 1) = strcmp(get(gca,'ylimmode'), 'auto');

plot(sin(x))
hold on
plot(-sin(x))
% yl = get(gca,'ylim');
% pass1(length(pass1) + 1) = ( norm(yl - [-1 1]) < tol );
pass1(length(pass1) + 1) = strcmp(get(gca,'ylimmode'), 'auto');
hold off

plot(.62*sin(x))
hold on
plot(-.62*sin(x))
yl = get(gca,'ylim');
pass1(length(pass1) + 1) = (norm(0.62*[-1 1] - yl) > 0.05 );
pass1(length(pass1) + 1) = strcmp(get(gca,'ylimmode'), 'auto');
hold off

%% Unbounded functions

dom2 = [-inf 0];
pass2 = [];
f = chebfun(@(x) exp(x), dom2);
plot(f)
pass2(length(pass2) + 1) =  strcmp(get(gca,'xlimmode'), 'manual');
pass2(length(pass2) + 1) =  strcmp(get(gca,'ylimmode'), 'auto');
xl = get(gca,'xlim');
pass2(length(pass2) + 1)  = ( norm(xl - [-10 0]) < tol );

g = chebfun(@(x) -0.62*exp(x), dom2);
plot(g)
pass2(length(pass2) + 1) =  strcmp(get(gca,'ylimmode'), 'auto');
yl = get(gca,'ylim');
pass2(length(pass2) + 1)  = ( norm(yl(1) - 0.62) > 0.05 );

dom3 = [-20 20];
h = chebfun(@(x) cos(x), dom3);
hold on
plot(h,'r')
xl = get(gca,'xlim');
pass2(length(pass2) + 1) = ( norm(xl - dom3) < tol );
hold off

plot(h,'g')
hold on
plot(g)
xl = get(gca,'xlim');
pass2(length(pass2) + 1)  = ( norm(xl - dom3) < tol );
hold off

dom4 = [-2 2];
h = chebfun(@(x) cos(x), dom4);
plot(g)
hold on
plot(h,'r')
xl = get(gca,'xlim');
pass2(length(pass2) + 1) = ( norm(xl - [-10 2]) < tol );
hold off

plot(h,'g')
hold on
plot(g)
xl = get(gca,'xlim');
pass2(length(pass2) + 1)  = ( norm(xl - [-10 2]) < tol );
hold off
%% Functions that blow-up
pass3 = [];
f1 = 11.3*sin(x);
f2 = 1./x;

plot(f2)
pass3(length(pass3) + 1) =  strcmp(get(gca,'xlimmode'), 'manual');
pass3(length(pass3) + 1) =  strcmp(get(gca,'ylimmode'), 'manual');
xl = get(gca, 'xlim');
yl = get(gca, 'ylim');
% Check that we obtain reasonable ylimits
pass3(length(pass3) + 1)  = ( norm(yl(1) - 1/xl(2)) < tol && ( yl(2) < 10 ) );
hold on
plot(f1,'r')
% Check that we obtain reasonable ylimits
yl = get(gca, 'ylim');
pass3(length(pass3) + 1)  = ( norm(yl(1)) < tol && ( yl(2) > 10 ) );
hold off

% Do the plotting in reverse order
plot(f1,'g')
ylOld = get(gca, 'ylim');
hold on
plot(f2)
% Check that we obtain reasonable ylimits
ylNew = get(gca, 'ylim');
pass3(length(pass3) + 1) = ( (norm(ylOld(1) - ylNew(1)) < tol) && ...
    (ylNew(2) > 10) );
hold off

% Check x-lim and y-lim symmetry of x:
f = chebfun(@(x) x, [-inf inf]);
plot(f)
xl = get(gca, 'xlim');
yl = get(gca, 'ylim');
pass3(length(pass3) + 1) = (abs(sum(xl)) < 1e-10) && (abs(sum(yl)) < 1e-10);

%%
pass = [pass1, pass2, pass3];

close(hfig);

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function tf = areTheseToolboxesInstalled(requiredToolboxes)
%ARETHESETOOLBOXESINSTALLED takes a cell array of toolbox names and checks whether they are currently installed
% SYNOPSIS tf = areTheseToolboxesInstalled(requiredToolboxes)
%
% INPUT requiredToolboxes: cell array with toolbox names to test for. Eg.
%        {'MATLAB','Image Processing Toolbox'}
%
% OUTPUT tf: true or false if the required toolboxes are installed or not
%--------------------------------------------------------------------------
% http://stackoverflow.com/a/2061075/1485872
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri
%--------------------------------------------------------------------------

% get all installed toolbox names
v = ver;
% collect the names in a cell array
[installedToolboxes{1:length(v)}] = deal(v.Name);

% check
tf = all(ismember(requiredToolboxes,installedToolboxes));

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function VRes = spm_write_residuals(SPM,Ic)
% Write residual images
% FORMAT Vres = spm_write_residuals(SPM,Ic)
% SPM    - structure containing generic analysis details
% Ic     - contrast index used to adjust data (0:   no adjustment)
%                                             (NaN: adjust for everything)
%
% VRes   - struct array of residual image handles
%__________________________________________________________________________
% Copyright (C) 2012-2013 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_write_residuals.m 6656 2015-12-24 16:49:52Z guillaume $

%-Get SPM.mat
%--------------------------------------------------------------------------
 if ~nargin || isempty(SPM)
     [SPM,sts] = spm_select(1,'^SPM\.mat$','Select SPM.mat');
     if ~sts, VRes = ''; return; end
 end

if ~isstruct(SPM)
    swd = spm_file(SPM,'fpath');
    try
        load(fullfile(swd,'SPM.mat'));
        SPM.swd = swd;
    catch
        error(['Cannot read ' fullfile(swd,'SPM.mat')]);
    end
end

try, SPM.swd; catch, SPM.swd = pwd; end
cwd = pwd; cd(SPM.swd);

%-Get contrast used to adjust data
%--------------------------------------------------------------------------
if nargin<2 || isempty(Ic)
    q(1)   = 0;
    Con    = {'<don''t adjust>'};
    q(2)   = NaN;
    Con{2} = '<adjust for everything>';
    for i = 1:length(SPM.xCon)
        if strcmp(SPM.xCon(i).STAT,'F')
            q(end + 1) = i;
            Con{end + 1} = SPM.xCon(i).name;
        end
    end
    i  = spm_input('adjust data for (select contrast)','!+1','m',Con);
    Ic = q(i);
end

%-Compute and write residuals
%==========================================================================
spm('Pointer','Watch')
M   = SPM.xY.VY(1).mat;
DIM = SPM.xY.VY(1).dim(1:min(numel(SPM.xY.VY(1).dim),3));
[nScan, nBeta] = size(SPM.xX.X);

if spm_mesh_detect(SPM.xY.VY)
    file_ext = '.gii';
else
    file_ext = spm_file_ext;
end

%-Initialise residual images
%--------------------------------------------------------------------------
VRes(1:nScan) = deal(struct(...
    'fname',    [],...
    'dim',      DIM,...
    'dt',       [spm_type('float64') spm_platform('bigend')],...
    'mat',      M,...
    'pinfo',    [1 0 0]',...
    'descrip',  'Residuals'));
for i = 1:nScan
    VRes(i).fname   = [sprintf('Res_%04d', i) file_ext];
    VRes(i).descrip = sprintf('Residuals (%04d)', i);
end
VRes = spm_data_hdr_write(VRes);

%-Loop over chunks
%--------------------------------------------------------------------------
chunksize = floor(spm_get_defaults('stats.maxmem') / 8 / nScan);
nbchunks  = ceil(prod(DIM) / chunksize);
chunks    = min(cumsum([1 repmat(chunksize,1,nbchunks)]),prod(DIM)+1);

spm_progress_bar('Init',nbchunks,'Writing residuals','Chunks');

for i=1:nbchunks
    chunk = chunks(i):chunks(i+1)-1;

    %-Get mask
    %----------------------------------------------------------------------
    m = spm_data_read(SPM.VM,chunk) > 0;
    m = m(:)';

    %-Get raw data, whiten and filter
    %----------------------------------------------------------------------
    y = zeros(nScan,numel(chunk));
    for j=1:nScan
        y(j,:) = spm_data_read(SPM.xY.VY(j),chunk);
    end
    y(:,~m) = [];

    y = spm_filter(SPM.xX.K,SPM.xX.W*y);

    if Ic ~= 0

        %-Parameter estimates: beta = xX.pKX*xX.K*y
        %------------------------------------------------------------------
        beta = zeros(nBeta,numel(chunk));
        for j=1:nBeta
            beta(j,:) = spm_data_read(SPM.Vbeta(j),chunk);
        end
        beta(:,~m) = [];

        %-Subtract Y0 = XO*beta,  Y = Yc + Y0 + e
        %------------------------------------------------------------------
        if ~isnan(Ic)
            y = y - spm_FcUtil('Y0',SPM.xCon(Ic),SPM.xX.xKXs,beta);
        else
            y = y - SPM.xX.xKXs.X * beta;
        end

    end

    %-Write residuals
    %----------------------------------------------------------------------
    yy = NaN(numel(chunk),1);
    for j=1:nScan
        yy(m)   = y(j,:);
        VRes(j) = spm_data_write(VRes(j), yy, chunk);
    end

    spm_progress_bar('Set',i)
end

cd(cwd);
spm_progress_bar('Clear')
spm('Pointer','Arrow')

author: spm


text:
function [Mean_accuracy, Variance,p1,m1,Means,Vars,Maxsf,Minsf]=simo2(data,v,c, measure, p, m, N,rn,pl)
%Inputs:
% data: put your data file in matrix form, rows indicates samples and
% columns features of the data.
% v:   how many columns of features you are using e.g. [1:4] means you will
% use first four columns as features of your data.
% c: column where you class labels are. They should be in numerical form.
% measure: Which quantifier in OWA operators you want to use.
% 1=Basic RIM quantifier
% 2=Quadratic quantifier
% 3=Exponential quantifier
% 4=Trigonometric quantifier
% 5=O'Hagans method.
% p: parameter in generalized Lukasiewics similarity, can be studied as a range of
% parameter values and then given as a vector i.e. p=[0.1:0.1:5].
% alpha: Alpha value in owa operators
% Can be given as vector i.e. m=[0.1:0.1:5].
% N: how many times data is divided randomly into testing set and learning
% set.
% rn: portion of data which is used in learning set. default rn=1/2 (data
% is divided in half; half is used in testing set and half in learning set.
%pl: do you want to plot how the parameter changes in p and m changes the
% mean classification accuracies and variances.
%
%OUTPUTS:
% Mean_accuracy: Mean classification accuracy with best parameter values in
% p and m:
% p1 and m1: best parameter values w.r.t. mean classification accuracy.
% Variance: Variances with best parameter values

sv_name = 'results2'; % mat file where results are stored ('' = "no storing")

[data, lc, cs] = init_data(data,v,c); % data initialization

w_opt = 0; % Do we use weight optimization. Not implemented in this version.

%Initializations
fitness = zeros(1,N);
fitness_id = zeros(1,N);
fitness_dif = zeros(1,N);
Means = zeros(length(p),length(m),length(rn));
Vars = zeros(length(p),length(m),length(rn));
Maxsf = zeros(length(p),length(m),length(rn));
Minsf = zeros(length(p),length(m),length(rn));

Means_fit_dif = zeros(length(p),length(m),length(rn));
Vars_fit_dir = zeros(length(p),length(m),length(rn));
Ideal_var = zeros(length(p),length(m), length(lc),length(rn));

for n = 1 : length(rn)
    rn_ideal = ones(1,length(lc))*rn(n);
    for j = 1:length(m)
        for i = 1:length(p)
            y = [p(i), m(j),measure]; % p and m values and similarity measure
            for k = 1 : N
                ideal_ind = [];
                for l = 1 :length(lc)
                    temp = randperm(lc(l))-1;
                    ideal_ind = [ideal_ind, cs(l)+temp(1:floor(lc(l)*rn_ideal(l)))]; % learning set indexes
                    data_ind = setxor([1:size(data,1)], ideal_ind); % testing set indexes
                end
                ideals(:,:,k) = idealvectors(data(ideal_ind,:), y); % idealvectors
                if w_opt == 0
                    [fitness(k), class, Simil] = calcfit(data(data_ind,:), ideals(:,:,k), y);
                end
            end
            Means(i,j,n) = mean(fitness);
            Vars(i,j,n) = var(fitness);
            Maxsf(i,j,n) = max(fitness);
            Minsf(i,j,n) = min(fitness);
            fitness=[];
        end
    end
    tmp=max(max(Means));
    [p1,m1]=find(tmp==Means);
    Mean_accuracy=Means(p1,m1);
    Variance=Vars(p1,m1);
    if pl==1
    [X,Y] = meshgrid(m,p);
    figure
    subplot(2,2,2);
    surfc(X,Y,Vars(:,:,n))
    title('Variances','FontSize',15)
    xlabel('\alpha-values')
    ylabel('p-values')
    zlabel('Variance')
    subplot(2,2,1)
    surfc(X,Y,Means(:,:,n))
    title('Mean classification accuracies','FontSize',15)
    xlabel('\alpha-values')
    ylabel('p-values')
    zlabel('Classification accuracy')
     subplot(2,2,3)
     surfc(X,Y,Maxsf(:,:,n))
     title('Maximum accuracies','FontSize',15)
     xlabel('\alpha-values')
     ylabel('p-values')
     zlabel('Classification accuracy')
     subplot(2,2,4)
     surfc(X,Y,Minsf(:,:,n))
     title('Minimum accuracies','FontSize',15)
     xlabel('\alpha-values')
     ylabel('p-values')
     zlabel('Classification accuracy')
    end
    clear Y
end
if length(sv_name) ~= 0
    save(sv_name)
end

author: Sable

repo: mcbench-benchmarks


text:
function [sourcemodel] = lbex(cfg, sourcemodel)

% This function will add the field "subspace" to the sourcemodel definition.
%
% The subspace projection is based on the LBEX (local basis expansion)
% method.

% set the defaults
cfg.lbex       = ft_getopt(cfg, 'lbex',       3); % this is a distance, in units of sourcemodel.pos
cfg.lbexeigtol = ft_getopt(cfg, 'lbexeigtol', 1000*eps);
cfg.feedback   = ft_getopt(cfg, 'feedback',   'text');
cfg.keep       = ft_getopt(cfg, 'keep',       'all');

if isequal(cfg.keep, 'all')
  cfg.keep = sourcemodel.inside;
else
  if ~islogical(cfg.keep)
    keep = false(size(sourcemodel.pos,1),1);
    keep(cfg.keep) = true;
    cfg.keep = keep;
  end
end
assert(isequal(numel(cfg.keep), numel(sourcemodel.inside)));

cfg.keep = cfg.keep(:) & sourcemodel.inside(:);

Ndipoles = size(sourcemodel.pos,1);
Ninside  = sum(cfg.keep);
inside   = find(cfg.keep);

% concatenate the leadfield of all dipoles that are inside the brain into one large matrix
lfa = cat(2, sourcemodel.leadfield{:});

% do the computations on the svd basis to avoid numerical issues
[U,S,V] = svd(lfa, 'econ');
diagS   = diag(S);
Tol     = 1e-12;
sel     = find(diagS>Tol.*diagS(1));
P       = diag(1./sqrt(diag(S(sel,sel))))*U(:,sel)'; % prewhitening matrix
lfa     = P*lfa;

% covariance of all leadfields
Ca = lfa * lfa';

sourcemodel.subspace = cell(1,size(sourcemodel.pos,1));
ft_progress('init', cfg.feedback, 'computing lbex');
for dipindx=1:Ninside
  % renumber the loop-index variable to make it easier to print the progress bar
  i = inside(dipindx);

  % compute the distance from this dipole to each other dipole
  dist = sqrt(sum((sourcemodel.pos-repmat(sourcemodel.pos(i,:), [Ndipoles 1])).^2, 2));

  % define the region of interest around this dipole
  sel  = dist<=cfg.lbex & sourcemodel.inside;
  Nsel = sum(sel);

  % concatenate the leadfield of all dipoles that are inside the ROI into one matrix
  lfr = P*cat(2,sourcemodel.leadfield{sel(:)'});

  % covariance of leadfields of dipoles inside the ROI
  Cr = lfr * lfr';

  % The eigenvalue problem is to determine the nontrivial solutions of the equation
  %   A*x = l*x
  % The generalized eigenvalue decomposition solves
  %   A*x = l*B*x
  % If B is non-singular, the problem could be solved by reducing it into a standard eigenvalue problem
  %   inv(B)*A*x = l*x
  % See http://www.mathworks.com/access/helpdesk/help/techdoc/ref/eig.html

  % compute eigenspace decomposition, THIS IS NUMERICALLY UNSTABLE
  [v, d] = eig(Cr, Ca);
  % compute eigenspace decomposition, THIS ALSO DOES NOT SOLVE THE PROBLEM
  % [v, d] = eig(pinv(Ca, cfg.lbexeigtol)*Cr);

  % select the eigenvectors with a non-zero eigenvalue
  dd  = diag(d);
  dd  = dd./max(dd);
  sel = dd>cfg.lbexeigtol;
  Nsel2 = sum(sel);

  ft_progress(dipindx/Ninside, 'computing lbex %d/%d, number of dipoles in ROI=%d, subspace dimension=%d\n', dipindx, Ninside, Nsel, Nsel2);

  % remember the subspace projection matrix
  sourcemodel.subspace{inside(dipindx)} = flip(v(:, sel)',1)*P;
end
ft_progress('close');

if ~isequal(cfg.keep(:),sourcemodel.inside(:))
  sourcemodel.inside = cfg.keep;
  sourcemodel.leadfield(~cfg.keep) = {[]};
  sourcemodel.subspace(~cfg.keep)  = {[]};
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function r = overlap_ratio(rect1, rect2)
% OVERLAP_RATIO
% Compute the overlap ratio between two rectangles
%
% Hyeonseob Nam, 2015
%

inter_area = rectint(rect1,rect2);
union_area = rect1(:,3).*rect1(:,4) + rect2(:,3).*rect2(:,4) - inter_area;

r = inter_area./union_area;
end

repo: MFT


text:
function pde = TorusTime(mum,mup,sigm,sigp,epsm,epsp,omega,x0,y0,z0,r1,r2)
%% USAGE: polynomial solution for Poisson equation
%  Last Modified: 02/21/2020 by Xu Zhang

%% PDE Structure
pde = struct('intf',@intf,...
    'f1',@f1,'f2',@f2,'f3',@f3,...
    'fm1',@fm1,'fm2',@fm2,'fm3',@fm3,...
    'fp1',@fp1,'fp2',@fp2,'fp3',@fp3,...
    'Mu',@Mu,'Mum',@Mum,'Mup',@Mup,'one',@one,...
    'Epslon',@Epslon,'Epslonm',@Epslonm,'Epslonp',@Epslonp,...
    'Sig',@Sig,'Sigm',@Sigm,'Sigp',@Sigp);

pde.mum = mum;
pde.mup = mup;
pde.sigm = sigm;
pde.sigp = sigp;
pde.epsm = epsm;
pde.epsp = epsp;
%% interface function
    function u = intf(x,y,z)
        u = (sqrt((x-x0).^2+(y-y0).^2)-r2).^2 + (z-z0).^2 - r1^2;
    end

%% exact solution

%% right hand side function
    function u = f1(x,y,z,t)
        u = fm1(x,y,z,t);
        id = intf(x,y,z) > 0;
        u(id) = fp1(x(id),y(id),z(id),t);
    end
    function u = f2(x,y,z,t)
        u = fm2(x,y,z,t);
        id = intf(x,y,z) > 0;
        u(id) = fp2(x(id),y(id),z(id),t);
    end
    function u = f3(x,y,z,t)
        u = fm3(x,y,z,t);
        id = intf(x,y,z) > 0;
        u(id) = fp3(x(id),y(id),z(id),t);
    end

    function u = fm1(x,y,z,t)
        u = -2*pi*omega*sin(myacos(x,y))*cos(2*pi*omega*t);
    end
    function u = fm2(x,y,z,t)
        u = 2*pi*omega*cos(myacos(x,y))*cos(2*pi*omega*t);
    end
    function u = fm3(x,y,z,t)
        u = zeros(size(x));
    end
    function u = fp1(x,y,z,t)
        u = zeros(size(x));
    end
    function u = fp2(x,y,z,t)
        u = zeros(size(x));
    end
    function u = fp3(x,y,z,t)
        u = zeros(size(x));
    end
    function theta = myacos(x,y)
        theta = (2*pi*(1-sign(y))/2 + sign(y).*acos(x./sqrt(x.^2+y.^2))).*abs(sign(y))+...
            (1-abs(sign(y))).*((1-x)/2*pi);
    end

%% Diffusion coefficient function
    function u = Mu(x,y,z)
        u = Mum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Mup(x(id),y(id),z(id));
    end
    function u = Mum(x,y,z)
        u = mum^(-1)*ones(size(x));
    end
    function u = Mup(x,y,z)
        u = mup^(-1)*ones(size(x));
    end

%% Mass coefficient function
    function u = Epslon(x,y,z)
        u = Epslonm(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Epslonp(x(id),y(id),z(id));
    end
    function u = Epslonm(x,y,z)
        u = epsm*ones(size(x));
    end
    function u = Epslonp(x,y,z)
        u = epsp*ones(size(x));
    end

    function u = Sig(x,y,z)
        u = Sigm(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Sigp(x(id),y(id),z(id));
    end
    function u = Sigm(x,y,z)
        u = sigm*ones(size(x));
    end
    function u = Sigp(x,y,z)
        u = sigp*ones(size(x));
    end

%% Other function
    function u = one(x,y,z)
        u = ones(size(x));
    end
end

repo: ifem


text:
function [Q2,Pk2]=EkfFilter(Q1,ImuData,t,Vm,Pk1)
%EKF filter to Gyro atitude with Accelerate & Magnetic
% derivation <A Double-Stage Kalman Filter for Orientation Tracking
%            with an Integrated Processor in 9-D IMU>
% author  Zhang Xin

if isempty(Pk1)
   Pk1=[ 0.1 , 0.01 , 0.01 , 0.01;...
           0.01 , 0.1 , 0.01 , 0.01;...
           0.01 , 0.01 , 0.1 , 0.01;...
           0.01 , 0.01 , 0.01 , 0.1  ]*0.001;

end
para1=0.0001;%0.0001;
para2=0.1;
para3=0.5;
Qk=para1*eye(4);
Rk1=para2*eye(3);
Rk2=para3*eye(3);

wx=ImuData(1,5)*t;
wy=ImuData(1,6)*t;
wz=ImuData(1,7)*t;

norm_a=norm(ImuData(1,2:4));
norm_g=norm(ImuData(1,5:7));

Ak=[ 1    , -wx/2 , -wy/2 , -wz/2 ;...
     wx/2 ,   1   ,  wz/2 , -wy/2 ;...
     wy/2 , -wz/2 ,   1   ,  wx/2 ;...
     wz/2 ,  wy/2 , -wx/2 ,   1   ];

Qp=Ak*Q1';     %Q_predict

Qp=Qp/norm(Qp);

P_k=Ak*Pk1*Ak'+Qk;

if abs(norm_a-9.8)<2 && norm_g< 2
    %R = quatern2rotMat(Qp);
    R2=[2*Qp(2)*Qp(3)+2*Qp(1)*Qp(4);...
        2*Qp(1)^2+2*Qp(3)^2-1;...
        2*Qp(3)*Qp(4)-2*Qp(1)*Qp(2)];

    R3=[2*Qp(2)*Qp(4)-2*Qp(1)*Qp(3);...
        2*Qp(3)*Qp(4)+2*Qp(1)*Qp(2);...
        2*Qp(1)^2+2*Qp(4)^2-1 ];

    J2=2*[ Qp(4), Qp(3), Qp(2) , Qp(1);...
           Qp(1),-Qp(2), Qp(3) ,-Qp(4);...
          -Qp(2),-Qp(1), Qp(4) , Qp(3)] ;

    J3=2*[-Qp(3), Qp(4),-Qp(1) , Qp(2);...
           Qp(2), Qp(1), Qp(4) , Qp(3);...
           Qp(1),-Qp(2),-Qp(3) , Qp(4)] ;

    h1=R3;                               %acc in sensor fixed frame
    h2=Vm(2)*R2+Vm(3)*R3;          %mag in sensor fixed frame

    Hk1=J3;
    Hk2=Vm(2)*J2+Vm(3)*J3;

    Kk1=P_k*Hk1'*inv(Hk1*P_k*Hk1'+Rk1);

    qe1=Kk1*(ImuData(1,2:4)'/norm_a-h1);         % Acc difference value between real & prediction from attitude

    Pk_1=(eye(4)-Kk1*Hk1)*P_k;

    Kk2=P_k*Hk2'*inv(Hk2*P_k*Hk2'+Rk2);

    mag=ImuData(1,8:10);
    mag=mag/norm(mag);

    qe2=Kk2*(mag'-h2);       % mag difference value bwtween real & prediction from attitude
    Pk2=(eye(4)-Kk2*Hk2)*Pk_1;

else
    qe1=[0;0;0;0];
    Pk_1=P_k;
    qe2=[0;0;0;0];
    Pk2=Pk_1;

end

% Kk2=P_k*Hk2'*inv(Hk2*P_k*Hk2'+Rk2);
%
% mag=ImuData(1,8:10);
% mag=mag/norm(mag);
%
% qe2=Kk2*(mag'-h2);       % mag difference value bwtween real & prediction from attitude

qt=Qp+qe1+qe2;

Q2=qt'/norm(qt);

if Q2(1)<0
    Q2=-Q2;
end

%Pk2=(eye(4)-Kk2*Hk2)*Pk_1;

end

author: shenshikexmu

repo: IMUCalibration-Gesture


text:
function y = colsum(x)

y = sum(x,2);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [inputMatrix, coords] = cropMatrix(inputMatrix,varargin)
	% Crops a matrix either by removing rows or adding NaNs to where data was previously.
	% Biafra Ahanonu
	% 2014.01.23 [16:06:01]
	% inputs
		% inputMatrix - a [m n p] matrix of any class type
	% outputs
		% inputMatrix - cropped or NaN'd matrix, same name to reduce memory usage

	[inputMatrix, coords] = ciapkg.image.cropMatrix(inputMatrix,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
classdef ParamMPCGCBF < handle
    properties
        horizon
        gamma
        P % weight for terminal cost (CLF)
        xWeight
        uWeight
    end
    methods
        function self = ParamMPCGCBF(horizon, gamma, P, x_weight, u_weight)
            self.horizon = horizon;
            self.gamma = gamma;
            self.P = P;
            self.xWeight = x_weight;
            self.uWeight = u_weight;
        end
    end
end

author: HybridRobotics

repo: NMPC-DCLF-DCBF


text:
%Generate 2D Gabor Patch
%
% GABOR
%   Generates a spatially oriented sinusoidal grating multiplied by a gaussian
%   window. Gabor functions have been used to model simple cell receptive fields
%   and are commonly used filters for edge detection.
%
% USAGE
%   [x y F] = gabor(varargin)
%
%   Returns:
%       x: grid of x-values
%       y: grid of y-values
%       F: gabor function evaluated at (x,y) points
%
%   Parameters (Default):
%       theta (2*pi*rand): orientation of the gabor patch
%       lambda (20): spatial wavelength
%       Sigma (10): standard deviation of gaussian window
%       width (256): width of generated image
%       height (256): height of generated image
%       px (rand): horizontal center of gabor patch, relative to the image width (must be between 0 and 1)
%       py (rand): vertical center of gabor patch, relative to the image height (must be between 0 and 1)
%
% EXAMPLE
%   [x y F] = gabor;
%   pcolor(x,y,F); axis image;
%   shading('interp'); colormap copper;
%
% VERSION 1.0, Thu Jul 12 09:47:52 2012
%
% AUTHOR: Niru Maheswaranathan
%         nirum@stanford.edu

function [x y F] = gabor(varargin)

    % Parse Input
    p = inputParser;
    addParamValue(p,'theta',2*pi*rand,@isnumeric);
    addParamValue(p,'lambda',20,@isnumeric);
    addParamValue(p,'Sigma',10,@isnumeric);
    addParamValue(p,'width',256,@isnumeric);
    addParamValue(p,'height',256,@isnumeric);
    addParamValue(p,'px',rand*0.8 + 0.1,@isnumeric);
    addParamValue(p,'py',rand*0.8 + 0.1,@isnumeric);
    p.KeepUnmatched = true;
    parse(p,varargin{:});

    % Generate mesh
    [x y] = meshgrid(1:p.Results.width, 1:p.Results.height);

    % Center of gaussian window
    cx = p.Results.px*p.Results.width;
    cy = p.Results.py*p.Results.height;

    % Orientation
    x_theta=(x-cx)*cos(p.Results.theta)+(y-cy)*sin(p.Results.theta);
    y_theta=-(x-cx)*sin(p.Results.theta)+(y-cy)*cos(p.Results.theta);

    % Generate gabor
    F = exp(-.5*(x_theta.^2/p.Results.Sigma^2+y_theta.^2/p.Results.Sigma^2)).*cos(2*pi/p.Results.lambda*x_theta);

end


text:
function lines2d(varargin)
%LINES2D  Description of functions operating on planar lines.
%
%   The term 'line' refers to a planar straight line, which is an unbounded
%   curve. Line segments defined between 2 points, which are bounded, are
%   called 'edge', and are presented in file 'edges2d'.
%
%   A straight line is defined by a point (its origin), and a vector (its
%   direction). The parameters are bundled into a 1-by-4 row vector:
%   LINE = [x0 y0 dx dy];
%
%   A line contains all points (x,y) such that:
%       x = x0 + t*dx
%       y = y0 + t*dy;
%   for all t between -infinity and +infinity.
%
%   See also
%   points2d, vectors2d, edges2d, rays2d
%   createLine, cartesianLine, medianLine, edgeToLine, lineToEdge
%   orthogonalLine, parallelLine, bisector, radicalAxis
%   lineAngle, linePosition, projPointOnLine
%   isPointOnLine, distancePointLine, isLeftOriented
%   intersectLines, intersectLineEdge, clipLine
%   reverseLine, transformLine, drawLine
%   lineFit

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2008-10-13, using Matlab 7.4.0.287 (R2007a)
% Copyright 2008-2022 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas

help('lines2d');

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% caffe.reset_all();
if ~exist('imageList','var')
    BLUFR_config = 'D:\face project\BLUFR\config\lfw\blufr_lfw_config.mat';
    load(BLUFR_config);
end;
lfw_root = 'E:\datasets\lfw-aligned-wuxiang';
caffe.set_mode_gpu();
gpu_id = 0;  % we will use the first gpu in this demo
caffe.set_device(gpu_id);

% ROIx = 19:82;
% ROIy = 19:82;
ROIx = 9:136;
ROIy = 9:136;
feature_dim = 256;
mean_value = 0;
scale = 1 / 255;
batch_size = 50;
total_length = length(imageList);
total_iter = ceil(total_length / batch_size);

% ROIx = 1:64;
% ROIy = 1:64;
height = length(ROIx);
width = length(ROIy);

% net = caffe.Net('D:\face project\face_verification_experiment\proto\LightenedCNN_B_deploy.prototxt','D:\face project\experiment\fine-tune\wuxiang\face_train_test_iter_14000.caffemodel', 'test');
net = caffe.Net('D:\face project\face_verification_experiment\proto\LightenedCNN_B_deploy.prototxt','D:\face project\experiment\fine-tune\wuxiang\LightenedCNN_B.caffemodel', 'test');
Descriptors = zeros(total_iter * batch_size, feature_dim);
for i = 1 : total_iter
    disp([i total_iter]);
    J = zeros(height,width,1,batch_size,'single');
    for j = 1 : batch_size
        if (i-1)*batch_size+j > total_length
            break;
        end;
        file_name = imageList{(i-1)*batch_size+j};
        find_zero = strfind(file_name,'0');
        folder = file_name(1:find_zero(1)-2);

        I = imread(fullfile(lfw_root, folder, file_name));
        I = rgb2gray(I)';
        I = I(ROIx,ROIy,:);
        I = single(I) - mean_value;
        J(:,:,:,j) = I*scale;
    end;
    f1 = net.forward({J});
    f1 = f1{1};
    Descriptors((i-1)*batch_size+1:i*batch_size,:) = reshape(f1,[feature_dim,batch_size])';
end;
Descriptors = Descriptors(1:total_length,:);
save('lfw_feature.mat','Descriptors');

author: happynear

repo: FaceVerification


text:
function r8poly_print ( m, a, title )

%*****************************************************************************80
%
%% R8POLY_PRINT prints out a polynomial.
%
%  Discussion:
%
%    The power sum form is:
%
%      p(x) = a(1) + a(2) * x + ... + a(m-1) * x^(m-1) + a(m) * x^(m-1)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the nominal degree of the polynomial.
%
%    Input, real A(1:M+1), the polynomial coefficients.
%    A(1) is the constant term and
%    A(M+1) is the coefficient of X^(M).
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );

  if ( a(m+1) < 0.0 )
    plus_minus = '-';
  else
    plus_minus = ' ';
  end

  mag = abs ( a(m+1) );

  if ( 3 <= m + 1 )
    fprintf ( 1, '  p(x) = %c%14g * x^%d\n', plus_minus, mag, m );
  elseif ( m + 1 == 2 )
    fprintf ( 1, '  p(x) = %c%14g * x\n', plus_minus, mag );
  elseif ( m + 1 == 1 )
    fprintf ( 1, '  p(x) = %c%14g\n', plus_minus, mag );
  end

  for i = m : -1 : 1

    if ( a(i) < 0.0 )
      plus_minus = '-';
    else
      plus_minus = '+';
    end

    mag = abs ( a(i) );

    if ( mag ~= 0.0 )

      if ( 3 <= i )
        fprintf ( 1, '         %c%14g * x^%d\n', plus_minus, mag, i-1 );
      elseif ( i == 2 )
        fprintf ( 1, '         %c%14g * x\n', plus_minus, mag );
      elseif ( i == 1 )
        fprintf ( 1, '         %c%14g\n', plus_minus, mag );
      end
    end

  end

  return
end


text:
% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin
% if you have a gaming joystick connected to your computer you use the
% sticks to control the position and orientation of a coordinate frame.

T = eye(4,4);
h = trplot(T, 'axis', [-5 5 -5 5 -5 5]);

while true
    T = joy2tr(T, 'tool');
    trprint(T, 'fmt', '%.1f')
    trplot(h, T);
end

author: Allopart

repo: rbpf-gmapping


text:
function test_bug2834

% MEM 2gb
% WALLTIME 00:20:00
% DEPENDENCY ft_sourceanalysis

% get some data
filename = dccnpath('/home/common/matlab/fieldtrip/data/test/latest/freq/meg/freq_mtmconvol_powandcsd_ctf151.mat');
load(filename);

filename = dccnpath('/home/common/matlab/fieldtrip/data/test/latest/vol/Subject01vol_singleshell.mat');
load(filename);

% the issue seems to be due to the fact that ft_selectdata with cfg.channel
% does not automatically subselect the corresponding channelcmbs

cfg = [];
cfg.channel = {'MEG' '-MLO21' '-MLF21'};
freq2 = ft_selectdata(cfg, freq);

% this works from the point of view of ft_selectdata
tmp = repmat(ft_channelselection(cfg.channel, freq.label), [1 149]);
tmp2 = tmp';

cfg = [];
cfg.channel = {'MEG' '-MLO21' '-MLF21'};
cfg.channelcmb = [tmp(tril(ones(149),-1)>0) tmp2(tril(ones(149),-1)>0)];
freq3 = ft_selectdata(cfg, freq);

cfg = [];
cfg.headmodel = vol;
cfg.sourcemodel.resolution = 1;
cfg.method = 'dics';
cfg.frequency = 14;
cfg.latency   = 0.5;
source = ft_sourceanalysis(cfg, freq3);


text:
%% NLS1 Original
clc
%Function
fun = @(x,xdata) x(1)*exp(x(2)*xdata);
grad = @(x,xdata) [ exp(x(2).*xdata), x(1).*xdata.*exp(x(2).*xdata)];
%Fitting Data
xdata = [0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3];
ydata = [455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5];
%Setup Options
opts = optiset('solver','matlab','display','iter');
%Build & Solve
Opt = opti('fun',fun,'grad',grad,'data',xdata,ydata,'ndec',2,'options',opts)
x0 = [100; -1]; % Starting guess
[x,fval,exitflag,info] = solve(Opt,x0)
%Plot
confidence(Opt)
plot(Opt)

%% NLS1 Weighted [no grad]
clc
%Function
fun = @(x,xdata) x(1)*exp(x(2)*xdata);
grad = @(x,xdata) [ exp(x(2).*xdata), x(1).*xdata.*exp(x(2).*xdata)];
%Fitting Data
xdata = [0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3];
ydata = [455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5];
%Fitting Weights
weights = ones(size(ydata)); weights(end) = 1e3;
%Setup Options
opts = optiset('solver','lmder','display','iter');
%Build & Solve
Opt = opti('fun',fun,'grad',[],'data',xdata,ydata,'weights',weights,'ndec',2,'options',opts)
x0 = [100; -1]; % Starting guess
[x,fval,exitflag,info] = solve(Opt,x0)
%Plot
[cInt,stats] = confidence(Opt)
plot(Opt)

%% NLS1 Weighted no xdata [no grad]
clc
%Fitting Data
xdata = [0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3]';
ydata = [455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5];
%Function
fun = @(x) x(1)*exp(x(2)*xdata);
grad = @(x) [ exp(x(2).*xdata), x(1).*xdata.*exp(x(2).*xdata)];
%Fitting Weights
weights = ones(size(ydata)); weights(end) = 1e3;
%Setup Options
opts = optiset('solver','lmder','display','iter');
%Build & Solve
Opt = opti('fun',fun,'grad',[],'ydata',ydata,'weights',weights,'ndec',2,'options',opts)
x0 = [100; -1]; % Starting guess
[x,fval,exitflag,info] = solve(Opt,x0)
%Plot
plot(Opt)

%% NLS1 Weighted [with grad]
clc
%Function
fun = @(x,xdata) x(1)*exp(x(2)*xdata);
grad = @(x,xdata) [ exp(x(2).*xdata), x(1).*xdata.*exp(x(2).*xdata)];
%Fitting Data
xdata = [0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3];
ydata = [455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5];
%Fitting Weights
weights = ones(size(ydata)); weights(end) = 1e3;
%Setup Options
opts = optiset('solver','lmder','display','iter');
%Build & Solve
Opt = opti('fun',fun,'grad',grad,'data',xdata,ydata,'weights',weights,'ndec',2,'options',opts)
x0 = [100; -1]; % Starting guess
[x,fval,exitflag,info] = solve(Opt,x0)
%Plot
plot(Opt)

%% NLS1 Weighted no xdata [with grad]
clc
%Fitting Data
xdata = [0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3]';
ydata = [455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5];
%Function
fun = @(x) x(1)*exp(x(2)*xdata);
grad = @(x) [ exp(x(2).*xdata), x(1).*xdata.*exp(x(2).*xdata)];
%Fitting Weights
weights = ones(size(ydata)); weights(end) = 1e3;
%Setup Options
opts = optiset('solver','lmder','display','iter');
%Build & Solve
Opt = opti('fun',fun,'grad',grad,'ydata',ydata,'weights',weights,'ndec',2,'options',opts)
x0 = [100; -1]; % Starting guess
[x,fval,exitflag,info] = solve(Opt,x0)
%Plot
plot(Opt)

%% NLS1 Weighted as NLP [no grad], MANUAL WEIGHTS
clc
%Fitting Data
xdata = [0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3]';
%Fitting Weights
weights = ones(size(xdata)); weights(end) = 1e3;
ydata = [455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5]'.*weights;
%Function
fun = @(x,xdata) x(1)*exp(x(2)*xdata).*weights;
grad = @(x,xdata) [ exp(x(2).*xdata), x(1).*xdata.*exp(x(2).*xdata)];
%Setup Options
opts = optiset('solver','ipopt','display','iter');
%Build & Solve
Opt = opti('fun',fun,'grad',[],'data',xdata,ydata,'weights',[],'ndec',2,'options',opts)
x0 = [100; -1]; % Starting guess
[x,fval,exitflag,info] = solve(Opt,x0)
%Plot
plot(Opt)

%% NLS1 Weighted as NLP [no grad], AUTO WEIGHTS
clc
%Fitting Data
xdata = [0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3]';
ydata = [455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5]';
%Fitting Weights
weights = ones(size(xdata)); weights(end) = 1e3;
%Function
fun = @(x,xdata) x(1)*exp(x(2)*xdata);
grad = @(x,xdata) [ exp(x(2).*xdata), x(1).*xdata.*exp(x(2).*xdata)];
%Setup Options
opts = optiset('solver','ipopt','display','iter');
%Build & Solve
Opt = opti('fun',fun,'grad',[],'data',xdata,ydata,'weights',weights,'ndec',2,'options',opts)
x0 = [100; -1]; % Starting guess
[x,fval,exitflag,info] = solve(Opt,x0)
%Plot
plot(Opt)

%% NLS1 Weighted as NLP no xdata [no grad]
clc
%Fitting Data
xdata = [0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3]';
ydata = [455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5];
%Function
fun = @(x) x(1)*exp(x(2)*xdata);
grad = @(x) [ exp(x(2).*xdata), x(1).*xdata.*exp(x(2).*xdata)];
%Fitting Weights
weights = ones(size(ydata)); weights(end) = 1e3;
%Setup Options
opts = optiset('solver','ipopt','display','iter');
%Build & Solve
Opt = opti('fun',fun,'grad',[],'ydata',ydata,'weights',weights,'ndec',2,'options',opts)
x0 = [100; -1]; % Starting guess
[x,fval,exitflag,info] = solve(Opt,x0)
%Plot
plot(Opt)

%% NLS1 Weighted as NLP [WITH grad], MANUAL WEIGHTS
clc
%Fitting Data
xdata = [0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3]';
%Fitting Weights
weights = ones(size(xdata)); weights(end) = 1e3;
ydata = [455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5]'.*weights;
%Function
fun = @(x,xdata) x(1)*exp(x(2)*xdata).*weights;
grad = @(x,xdata) [ exp(x(2).*xdata).*weights, x(1).*xdata.*exp(x(2).*xdata).*weights];
%Setup Options
opts = optiset('solver','ipopt','display','iter');
%Build & Solve
Opt = opti('fun',fun,'grad',grad,'data',xdata,ydata,'weights',[],'ndec',2,'options',opts)
x0 = [100; -1]; % Starting guess
[x,fval,exitflag,info] = solve(Opt,x0)
%Plot
plot(Opt)

%% NLS1 Weighted as NLP [WITH grad], AUTO WEIGHTS
clc
%Fitting Data
xdata = [0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3]';
ydata = [455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5]';
%Fitting Weights
weights = ones(size(xdata)); weights(end) = 1e3;
%Function
fun = @(x,xdata) x(1)*exp(x(2)*xdata);
grad = @(x,xdata) [ exp(x(2).*xdata), x(1).*xdata.*exp(x(2).*xdata)];
%Setup Options
opts = optiset('solver','ipopt','display','iter');
%Build & Solve
Opt = opti('fun',fun,'grad',grad,'data',xdata,ydata,'weights',weights,'ndec',2,'options',opts)
x0 = [100; -1]; % Starting guess
[x,fval,exitflag,info] = solve(Opt,x0)
%Plot
plot(Opt)

%% NLS1 Weighted as NLP no xdata [WITH grad]
clc
%Fitting Data
xdata = [0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3]';
ydata = [455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5];
%Function
fun = @(x) x(1)*exp(x(2)*xdata);
grad = @(x) [ exp(x(2).*xdata), x(1).*xdata.*exp(x(2).*xdata)];
%Fitting Weights
weights = ones(size(ydata)); weights(end) = 1e3;
%Setup Options
opts = optiset('solver','ipopt','display','iter');
%Build & Solve
Opt = opti('fun',fun,'grad',grad,'ydata',ydata,'weights',weights,'ndec',2,'options',opts)
x0 = [100; -1]; % Starting guess
[x,fval,exitflag,info] = solve(Opt,x0)
%Plot
plot(Opt)

%% NLS2
clc
%Function
i = (1:40)';
fun = @(x) x(1)*exp(-x(2)*i) + x(3);
%Fitting Data
ydata=[5.8728, 5.4948, 5.0081, 4.5929, 4.3574, 4.1198, 3.6843, 3.3642, 2.9742, 3.0237, 2.7002, 2.8781,...
       2.5144, 2.4432, 2.2894, 2.0938, 1.9265, 2.1271, 1.8387, 1.7791, 1.6686, 1.6232, 1.571, 1.6057,...
       1.3825, 1.5087, 1.3624, 1.4206, 1.2097, 1.3129, 1.131, 1.306, 1.2008, 1.3469, 1.1837, 1.2102,...
       0.96518, 1.2129, 1.2003, 1.0743];
%Setup Options
opts = optiset('solver','mkltrnls','display','iter');
%Build & Solve
x0=[1.0; 0.0; 0.0];
Opt = opti('fun',fun,'ydata',ydata,'ndec',3,'options',opts)
[x,fval,exitflag,info] = solve(Opt,x0)

plot(Opt)

%% NLS2 w Weights
clc
%Function
i = (1:40)';
fun = @(x) x(1)*exp(-x(2)*i) + x(3);
%Fitting Data
ydata=[5.8728, 5.4948, 5.0081, 4.5929, 4.3574, 4.1198, 3.6843, 3.3642, 2.9742, 3.0237, 2.7002, 2.8781,...
       2.5144, 2.4432, 2.2894, 2.0938, 1.9265, 2.1271, 1.8387, 1.7791, 1.6686, 1.6232, 1.571, 1.6057,...
       1.3825, 1.5087, 1.3624, 1.4206, 1.2097, 1.3129, 1.131, 1.306, 1.2008, 1.3469, 1.1837, 1.2102,...
       0.96518, 1.2129, 1.2003, 1.0743];
%Weighting
wts = ones(size(ydata)); wts([1 end]) = 1e3;
%Build & Solve
x0=[1.0; 0.0; 0.0];
Opt = opti('fun',fun,'ydata',ydata,'x0',x0);
OptW = opti('fun',fun,'ydata',ydata,'weights',wts,'x0',x0);
x = solve(Opt)
xw = solve(OptW)

plot(i,ydata,'ko',i(end),ydata(end),'ksq',i,fun(xw),'r*-',i,fun(x),'m*-')
xlim([length(i)*0.45 length(i)*1.01]); xlabel('i'); ylabel('y');
legend('Original Data','Target Point','Weighted Fit','Standard Fit');
title('NLS Curve Fit - Comparison of Weighted and Un-weighted');

%% NLS3 (Modified NLP - HS76)
clc
%Function
fun = @(x) [x(1);
            sqrt(0.5)*x(2);
            x(3);
            sqrt(0.5)*x(4);];
%Fitting Data
ydata = [0.0, 0.0, 0.0, 0.0];
%Constraints
lb = [0.0, 0.0, 0.0, 0.0];
ub = [inf, inf, inf, inf];
A = -[-1.0, -2.0, -1.0, -1.0;
     -3.0, -1.0, -2.0, 1.0];
b = -[-5.0, -0.4]';
Aeq = [0.0, 1.0, 4.0, 0.0];
beq = 1.5;
%Setup Options
opts = optiset('solver','levmar','display','iter');
%Build & Solve
x0 = [0.5, 0.5, 0.5, 0.5];
Opt = opti('fun',fun,'ydata',ydata,'ineq',A,b,'eq',Aeq,beq,'bounds',lb,ub,'options',opts)
[x,fval,exitflag,info] = solve(Opt,x0)

%% NLS Example Prob
clc
%Get Problem
prob = nls_prob(19);
opts = optiset('display','iter');
%Build OPTI Object
Opt = opti(prob,opts);
%Solve
[x,fval,exitflag,info] = solve(Opt)
%Plot
plot(Opt,[],1)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef PTKEditedResultPatch < PTKPatch
% PTKPatch. Class for sharing image edits with other PTK clients
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %
    %
    %     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
        Schema = 1
        PatchType = 'EditedResult'
        SeriesUid
        PluginName
        EditedResult
    end

end

author: tomdoel

repo: pulmonarytoolkit


text:
function result = combine_structs(args,structmask)
% Turn a list of name-value pairs/structs into a single struct.
% Struct = combine_structs(Arguments,StructMask)
%
% In:
%   Arguments : cell array of function arguments in form of name-value pairs, possibly interleaved
%               with structs {'name',value,STRUCT,'name',value,'name',value, ...}
%
%   StructMask : a bitmask that encodes at what positions Arguments contains structs; optional.
%
% Out:
%   Struct : A struct with values assigned to fields
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2013-11-03

if nargin < 2
    structmask = cellfun('isclass',args,'struct'); end

if isscalar(structmask)
    % inputs has already a single struct
    result = args{1};
else
    % obtain flat NVP list
    nvps = flatten_structs(args,structmask);
    % split names/values and find indices of the last assignment for each name
    nvps = reshape(nvps,2,[]);
    [s,indices] = sort(nvps(1,:));
    indices(strcmp(s(1:end-1),s(2:end))) = [];
    % build the struct
    result = cell2struct(nvps(2,indices),nvps(1,indices),2);
end

repo: REST


text:
%  Figure 6.68      Feedback Control of Dynamic Systems, 6e
%                   Franklin, Powell, Emami
%

clear all
%close all;
clf

k=10;
num=k;
den=[1 1 0];
num=conv(num,[10 1]);
den=conv(den,[100 1]);
numcl=[1 0.1];
dencl=[1 1.01 1.01 0.1];
t=0:.2:50;
y=step(numcl,dencl,t);
%subplot(2,1,1)
plot(t,y);
xlabel('Time (sec)');
ylabel('y');
title('Fig. 6.68 Step response of lag-compensation design.');
nicegrid;


text:
function [returnVal, txtReturnVal] = lvmThreeDPlot(X, lbl, symbol)

% LVMTHREEDPLOT Helper function for plotting the labels in 3-D.
% FORMAT
% DESC helper function for plotting an embedding in 3-D with symbols.
% ARG X : the data to plot.
% ARG lbl : the labels of the data point.
% ARG symbol : the symbols to use for the different labels.
%
% SEEALSO : lvmScatterPlot, lvmVisualise
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006, 2008

% MLTOOLS
if nargin < 2
  lbl = [];
end
if(strcmp(lbl, 'connect'))
  connect = true;
  lbl = [];
else
  connect = false;
end

if iscell(lbl)
  lblStr = lbl{2};
  lbl = lbl{1};
  labelsString = true;
else
  labelsString = false;
end

if nargin < 3
  if isempty(lbl)
    symbol = getSymbols(1);
  else
    symbol = getSymbols(size(lbl,2));
  end
end
axisHand = gca;
returnVal = [];
textReturnVal = [];
nextPlot = get(axisHand, 'nextplot');
for i = 1:size(X, 1)
  if i == 2
    set(axisHand, 'nextplot', 'add');
  end
  if ~isempty(lbl)
    labelNo = find(lbl(i, :));
  else
    labelNo = 1;
  end
  try
    returnVal = [returnVal; plot3(X(i, 1), X(i, 2), X(i, 3), symbol{labelNo})];
    if labelsString
      textReturnVal = [textReturnVal; text(X(i, 1), X(i, 2), X(i, 3), lblStr{i})];
    end
    if connect
      if i>1
        line([X(i-1, 1) X(i, 1)], [X(i-1, 2), X(i, 2)], [X(i-1, 3), X(i, 3)]);
      end
    end
  catch
    if strcmp(lasterr, 'Index exceeds matrix dimensions.')
      error(['Only ' num2str(length(symbol)) ' labels supported (it''s easy to add more!)'])
    end
  end
end
set(axisHand, 'nextplot', nextPlot);
set(returnVal, 'markersize', 10);
set(returnVal, 'linewidth', 2);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
clear all, close all, clc

s = tf('s');

G = (s+1)/(s-2);
Gtrue = (s+.9)/(s-1.9);

K = 1/G;

L = K*Gtrue;

margin(L)

CL = feedback(L,1);

author: dynamicslab


text:
function [minFlux, maxFlux, minFD, maxFD, GRvector, result, LP] = SteadyComFVA(modelCom, options, varargin)
% Flux variability analysis for community model at community steady-state for a range of growth rates.
% The function is capable of saving intermediate results and continuing from previous results
% if the file path is given in `options.saveFVA`. It also allows switch from single thread to parallel
% computation from intermediate results (but not the other way round).
%
% USAGE:
%    [minFlux, maxFlux, minFD, maxFD, GRvector, result, LP] = SteadyComFVA(modelCom, options, parameter, 'param1', value1, 'param2', value2, ...)
%
% INPUT:
%    modelCom:     A community COBRA model structure with the following fields (created using `createMultipleSpeciesModel`)
%                  (the first 5 fields are required, at least one of the last two is needed. Can be obtained using `getMultiSpecisModelId`):
%
%                    * S - Stoichiometric matrix
%                    * b - Right hand side
%                    * c - Objective coefficients
%                    * lb - Lower bounds
%                    * ub - Upper bounds
%                    * infoCom - structure containing community reaction info
%                    * indCom - the index structure corresponding to `infoCom`
%
% OPTIONAL INPUTS:
%    options:    struct with the following possible fields:
%
%                  * optGRpercent - A vector of percentages. Perform FVA at these percents of max. growth rate (Default = [99.99])
%                  * optBMpercent - Only consider solutions that yield at least a certain percentage of the optimal biomass (Default = 99.99)
%                  * rxnNameList - List of reactions (index row vector or subset of `*.rxns`) for which FVA is performed.
%                    (Default = biomass reaction of each species)
%                    Or a :math:`(N_{rxns} + N_{organism}) x K` matrix for FVA of `K` linear combinations of fluxes and/or abundances
%                    e.g., `[1; -2; 0]` for finding the max/min of :math:`1 v_1 - 2 v_2 + 0 v_3`
%                  * rxnFluxList - List of reactions (index vector or subset of `*.rxns`) whose fluxes are
%                    also returned along with the FVA result of each entry in `rxnNameList`
%                    (Default = biomass reaction of each species)
%                  * GRmax - maximum growth rate of the model (default to be found `SteadyCom.m`)
%                    (the two parameters below are usually determined by solving the problem during the program.
%                    Provide them only if you want to constrain the total biomass to a particular value)
%                  * BMmaxLB - lower bound for the total biomass (default 1)
%                  * BMmaxUB - upper bound for the total biomass (other parameters below)
%                  * saveFVA - If non-empty, become the filename to save the FVA results
%                    (default empty, not saving)
%                  * saveFre - save frequency. Save every `(#rxns for FVA) * saveFre` (default 0.1)
%                  * threads - for parallelization: > 1 for explicitly stating the no. of threads used,
%                    0 or -1 for using all available threads. Default 1.
%                    (Requires Matlab parallel toolbox)
%                  * verbFlag - Verbose output. 1 to have waitbar, >1 to have stepwise output (default 3)
%                  * loadModel - (`ibm_cplex` only) String of filename to be loaded. If non-empty, load the cplex
%                    model ('loadModel.mps'), basis ('loadModel.bas') and parameters ('loadModel.prm').
%                    (May add also other parameters in `SteadyCom` for calculating the maximum growth rate.)
%
%    parameter:  structure for solver-specific parameters.
%                'param1', value1, ...   name-value pairs for `solveCobraLP` parameters. See `solveCobraLP` for details
%
% OUTPUTS:
%    minFlux:    Minimum flux for each reaction
%    maxFlux:    Maximum flux for each reaction
%
% OPTIONAL OUTPUTS:
%    minFD:      :math:`rxnFluxList * rxnNameList` matrix containing the fluxes in `options.rxnFluxList`
%                corresponding to minimizing each reaction in `options.rxnNameList`
%    maxFD:      :math:`rxnFluxList * rxnNameList` matrix containing the fluxes in `options.rxnFluxList`
%                corresponding to maximizing each reaction in `options.rxnNameList`
%    GRvector:   a vector of growth rates at which FVA has been performed
%    result:     result structure from `SteadyCom`
%    LP:         `LP` problem structure (`Cplex LP` object for `ibm_cplex`)

[modelCom, ibm_cplex, feasTol, solverParams, parameters] = SteadyComSubroutines('initialize', modelCom, varargin{:});
% Initialization above
if nargin < 2 || isempty(options)
    options = struct();
end
% handle solveCobraLP name-value arguments that are specially treated in SteadyCom functions
[options, varargin] = SteadyComSubroutines('solveCobraLP_arg', options, parameters, varargin);

% get SteadyCom paramters. If a required parameter is in options, get its value, else equal to the
% default value in SteadyComSubroutines('getParams') if there is. Otherwise an empty matrix.
[GRmax, optGRpercent, rxnNameList, rxnFluxList, ...
    GRfx, BMmaxLB, BMmaxUB, ...
    verbFlag, loadModel, saveFVA, threads] = SteadyComSubroutines('getParams',  ...
    {'GRmax', 'optGRpercent', 'rxnNameList', 'rxnFluxList',...
    'GRfx','BMmaxLB','BMmaxUB', ...
    'verbFlag', 'loadModel','saveFVA','threads'}, ...
    options, modelCom);

[m, n] = size(modelCom.S);  % model size
nSp = numel(modelCom.indCom.spBm);  % number of organisms
nRxnSp = sum(modelCom.indCom.rxnSps > 0);  % number of organism-specific rxns

if ischar(rxnNameList)
    rxnNameList = {rxnNameList};
end
if iscell(rxnNameList) || (min(size(rxnNameList)) == 1 && size(rxnNameList, 1) < n)
    nRxnFVA = numel(rxnNameList);
else
    nRxnFVA = size(rxnNameList, 2);
end
if ischar(rxnFluxList)
    rxnFluxList = {rxnFluxList};
end

addRow = false;
GRgiven = false;
if isempty(GRmax)
    % get maximum growth rate
    [sol, result, LP] = SteadyCom(modelCom, options, varargin{:});
    if strcmp(result.stat,'infeasible')
        % infeasible model
        warning('The model is infeasible.');
        [minFlux, maxFlux] = deal(NaN(nRxnFVA, 1));
        [minFD, maxFD] = deal(NaN(numel(rxnFluxList), nRxnFVA));
        GRvector = NaN(numel(optGRpercent), 1);
        return
    end
    GRmax = result.GRmax;
    if ibm_cplex
        idRow = size(LP.Model.A, 1);  % row that constrains the total biomass
    else
        idRow = size(LP.A, 1);  % row that constrains the total biomass
    end
else
    % if GRmax is given, BMmaxLB and BMmaxUB should be included in options in this case to ensure feasibility
    if ibm_cplex && ~isempty(loadModel)
        % load Cplex model if loadModel is given
        LP = Cplex('SteadyComFVA');
        LP.readModel([loadModel '.mps']);
        LP.readBasis([loadModel '.bas']);
        LP.readParam([loadModel '.prm']);
        LP.DisplayFunc = [];
        fprintf('Load model ''%s'' successfully.\n', loadModel);
        addRow = true;
        if size(LP.Model.A, 1) > m + 2 * nRxnSp + nSp
            %try to find the row that constrains total biomass
            [ynRow, idRow] = ismember(sparse(ones(nSp, 1), (n + 1) : (n + nSp), ones(nSp, 1), 1, n + nSp),...
                LP.Model.A((m + 2 * nRxnSp + nSp + 1):end, 1:(n + nSp)),'rows');
            if ynRow
                idRow = m + 2 * nRxnSp + nSp + idRow;
            end
            addRow = ~ynRow;
        end
    else
        % get the LP structure using SteadyCom
        options2 = options;
        options2.LPonly = true;
        [~, ~, LP] = SteadyCom(modelCom, options2, varargin{:});
        addRow = true;  % no constraint on total biomass using LPonly option
    end
    result = struct('GRmax',GRmax,'vBM',[],'BM',[],'Ut',[],'Ex',[],'flux',[],'iter0',[],'iter',[],'stat','optimal');
    GRgiven = true;
end
if addRow
    % add a row for constraining the sum of biomass if not exist
    if ibm_cplex
        LP.addRows(BMmaxLB, sparse(ones(1, nSp), (n + 1) : (n + nSp), ones(1, nSp), 1, size(LP.Model.A, 2)), BMmaxUB, 'UnityBiomass');
        idRow = size(LP.Model.A, 1);
    else
        LP.A = [LP.A; sparse([ones(nSp, 1); 2 * ones(nSp, 1)], repmat((n + 1):(n + nSp), 1, 2),...
            ones(nSp * 2, 1), 2, size(LP.A, 2))];
        LP.b = [LP.b; BMmaxUB; BMmaxLB];
        LP.csense = [LP.csense, 'LG'];
        idRow = size(LP.A, 1);
    end
else
    % using BMmaxLB and BMmaxUB stored in the LP if not given in options
    if ~isfield(options,'BMmaxLB')  % take from LP if not supplied
        if ibm_cplex
            BMmaxLB = LP.Model.lhs(idRow);
        else
            BMmaxLB = LP.b(idRow);
        end
    end
    if ~isfield(options,'BMmaxUB')  % take from LP if not supplied
        if ibm_cplex
            BMmaxUB = LP.Model.rhs(idRow);
        else
            BMmaxUB = LP.b(idRow - 1);
        end
    end
    % not allow the max. biomass to exceed the one at max growth rate,
    % can happen if optBMpercent < 100. May dismiss this constraint or
    % manually supply BMmaxUB in the options if sum of biomass should be variable
    if ibm_cplex
        LP.Model.lhs(idRow) = BMmaxLB;
        LP.Model.rhs(idRow) = BMmaxUB;
    else
        LP.b(idRow) = BMmaxLB;
        LP.b(idRow - 1) = BMmaxUB;
    end
end
if ibm_cplex
    LP = setCplexParam(LP, solverParams);  % set Cplex parameters
    % update the LP to ensure the current growth rate is constrained
    LP.Model.A = SteadyComSubroutines('updateLPcom', modelCom, GRmax, GRfx, [], LP.Model.A, []);
    LP.Model.sense = 'minimize';
    LP.Model.obj(:) = 0;
    LP.solve();
    dev = checkSolFeas(LP);
else
    % update the LP to ensure the current growth rate is constrained
    LP.A = SteadyComSubroutines('updateLPcom', modelCom, GRmax, GRfx, [], LP.A, []);
    LP.c(:) = 0;
    LP.osense = 1;
    sol = solveCobraLP(LP, varargin{:});
    dev = checkSolFeas(LP, sol);
    if isfield(sol, 'basis')
        LP.basis = sol.basis;  % reuse basis
    end
end
% check and adjust for feasibility
% (LP from SteadyCom should pass this automatically as the row has been added in SteadyComCplex)
kBMadjust = 0;
while ~(dev <= feasTol) && kBMadjust < 10
    kBMadjust = kBMadjust + 1;
    % relax the required sum of biomass
    if ibm_cplex
        LP.Model.lhs(idRow) = BMmaxLB * (1 - feasTol/(11 - kBMadjust));
        LP.solve();
        dev = checkSolFeas(LP);
    else
        LP.b(idRow) = BMmaxLB * (1 - feasTol/(11 - kBMadjust));
        sol = solveCobraLP(LP, varargin{:});
        dev = checkSolFeas(LP, sol);
        if isfield(sol, 'basis')
            LP.basis = sol.basis;  % reuse basis
        end
    end
    if verbFlag
        fprintf('BMmax adjustment: %d\n', kBMadjust);
    end
end
if ~(dev <= feasTol)  % dev can be NaN, which still means infeasibility. So use ~(dev <= feasTol) instead of dev > feasTol
    warning('Model not feasible.')
    [minFlux, maxFlux] = deal(NaN(nRxnFVA, numel(optGRpercent)));
    [minFD, maxFD] = deal(NaN(numel(rxnFluxList), nRxnFVA, numel(optGRpercent)));
    GRvector = NaN(numel(optGRpercent), 1);
    result.stat = 'infeasible';
    return
end
if GRgiven
    % assign result structure if in the rare case of given growth rate
    if ibm_cplex
        flux = LP.Solution.x;
    else
        flux = sol.full;
    end
    result.vBM = flux(modelCom.indCom.spBm);
    result.BM = flux(n+1:n+nSp);
    % two different types of indexing
    if size(modelCom.indCom.EXcom, 2) == 2
        % uptake and excretion reactions separated
        result.Ut = flux(modelCom.indCom.EXcom(:,1));
        result.Ex = flux(modelCom.indCom.EXcom(:,2));
    else
        % uptake and excretion in one exchange reaction
        [result.Ut, result.Ex] = deal(flux(modelCom.indCom.EXcom(:,1)));
        result.Ut(result.Ut > 0) = 0;
        result.Ut = -result.Ut;
        result.Ex(result.Ex < 0) = 0;
    end
    result.flux = flux(1:n);
end
% update BMmaxLB and BMmaxUB for FVA at each given growth rate
if ~isfield(options, 'BMmaxLB')
    if ibm_cplex
        options.BMmaxLB = LP.Model.lhs(idRow);
    else
        options.BMmaxLB = LP.b(idRow);
    end
end
if ~isfield(options, 'BMmaxUB')
    if ibm_cplex
        options.BMmaxUB = LP.Model.rhs(idRow);
    else
        options.BMmaxUB = LP.b(idRow - 1);
    end
end

GRvector = GRmax * optGRpercent / 100;
if ~isempty(saveFVA)
    % decide number of digits in the save name
    if numel(optGRpercent) == 1
        kDisp = 2;
    else
        d = min(GRvector(2:end) - GRvector(1:end-1));
        if d < 1
            kDisp = abs(floor(log10(abs(d))));
        else
            kDisp = 0;
        end
    end
    directory = strsplit(saveFVA,filesep);
    if numel(directory) > 1
        % not saving in the current directory. Create the directory.
        directory = strjoin([{pwd}, directory(1:end-1)],filesep);
        if ~exist(directory, 'dir')
            mkdir(directory);
        end
    end
    if ibm_cplex
        LPmodel = LP.Model;  % the Cplex dynamic object is not good for saving
        LPstart = LP.Start;
        save([saveFVA '_model.mat'], 'LPmodel', 'LPstart', 'options', 'solverParams', 'parameters');
        clear LPmodel LPstart
    else
        save([saveFVA '_model.mat'], 'LP', 'options', 'solverParams', 'parameters');
    end
end

[minFlux, maxFlux] = deal(zeros(nRxnFVA, numel(GRvector)));
[minFD, maxFD] = deal(zeros(numel(rxnFluxList), nRxnFVA, numel(GRvector)));

%parallel computation
try
    p = gcp('nocreate');
    if isempty(p)
        if threads > 1
            %given explicit no. of threads
            parpool(ceil(threads));
        elseif threads ~= 1
            %default max no. of threads (input 0 or -1 etc)
            parpool;
        end
    end
catch
    %No parallel pool existent
end
%perform FVA at each growth rate
for j = 1:numel(GRvector)
    optionsJ = options;
    optionsJ.GR = GRvector(j);
    if ~isempty(saveFVA)
        optionsJ.saveFVA = sprintf(['%s_GR%.' num2str(kDisp) 'f'], saveFVA, GRvector(j));
    end
    [minFluxJ,maxFluxJ,minFDj,maxFDj,LP] = SteadyComFVAgr(modelCom,optionsJ, LP, varargin{:});
    minFlux(:, j) = minFluxJ;
    maxFlux(:, j) = maxFluxJ;
    minFD(:,:, j) = minFDj;
    maxFD(:,:, j) = maxFDj;
end

end

author: opencobra

repo: cobratoolbox


text:
classdef xt_reserves < mp.extension

%   MATPOWER
%   Copyright (c) 2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end     %% properties

    methods
        function dmc_elements = dmc_element_classes(obj, dmc_class, fmt, mpopt)
            switch fmt
                case 'mpc2'
                    dmc_elements = { @mp.dmce_reserve_gen_mpc2, ...
                                     @mp.dmce_reserve_zone_mpc2 };
                otherwise
                    dmc_elements = {};      %% no modifications
            end
        end

        function dm_elements = dm_element_classes(obj, dm_class, task_tag, mpopt)
            switch task_tag
                case 'OPF'
                    dm_elements = { @mp.dme_reserve_gen, ...
                                    @mp.dme_reserve_zone };
                otherwise
                    dm_elements = {};       %% no modifications
            end
        end

        function mm_elements = mm_element_classes(obj, mm_class, task_tag, mpopt)
            switch task_tag
                case {'OPF'}
                    mm_elements = { @mp.mme_reserve_gen, ...
                                    @mp.mme_reserve_zone };
                otherwise
                    mm_elements = {};       %% no modifications
            end
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [estBias, estVar] = testSpSegmentationEstimator(imsegs, spdata, maps)

% UNFINISHED
[data, lab] = formatData(spdata, {imsegs(:).labels});

vconf = test_boosted_dt_mc(vclassifier, data);
vconf = 1 ./ (1+exp(-vconf));
vconf = vconf ./ repmat(sum(vconf, 2), 1, size(vconf, 2));

hconf = test_boosted_dt_mc(hclassifier, data);
hconf = 1 ./ (1+exp(-hconf));
hconf = hconf ./ repmat(sum(hconf, 2), 1, size(hconf, 2));

spEst = [vconf(:, 1) repmat(vconf(:, 2), 1, 5).*hconf2 vconf(:, 3)];
spTrue = lab;

nmaps = 0;
for f = 1:numel(maps)
    nmaps = nmaps + numel(maps{f});
end

pEst = zeros(nmaps, 1);
pTrue = zeros(nmaps, 1);
%
% True = zeros(size(valEst));
% for i = 1:numel(lab)
%     for k =
%
% spvlabels = {imsegs(:).vert_labels};
% sphlabels = {imsegs(:).horz_labels};
%
% % test vertical
% [vdata, vlab] = formatData(spfeatures, spvlabels);
% vconf = test_boosted_dt_mc(vclassifier, vdata);
% vconf = 1 ./ (1+exp(-vconf));
% vconf = vconf ./ repmat(sum(vconf, 2), 1, size(vconf, 2));
% [tmp, vmax] = max(vconf, [], 2);
% verror = mean(vmax~=vlab)
%
% vconferr = 0;
% for k = 1:numel(vlab)
%     vconferr = vconferr + 1 - vconf(k, vlab(k));
% end
% vconferr = vconferr / numel(vlab)
%
% % test subclass
% [hdata, hlab] = formatData(spfeatures, sphlabels);
% hconf = test_boosted_dt_mc(hclassifier, hdata);
% hconf = 1 ./ (1+exp(-hconf));
% hconf = hconf ./ repmat(sum(hconf, 2), 1, size(hconf, 2));
% [tmp, hmax] = max(hconf, [], 2);
% herror = mean(hmax~=hlab)
%
% hconferr = 0;
% for k = 1:numel(hlab)
%     hconferr = hconferr + 1 - hconf(k, hlab(k));
% end
% hconferr = hconferr / numel(hlab)
%
% % get same_label error
% %[serror, sconferr] = sameLabelError(spfeatures, spvlabels, sphlabels, ...
% %    {imsegs(:).adjmat}, vclassifier, hclassifier)
%
% errors.verr = verror;
% errors.vcerr = vconferr;
% errors.herr = herror;
% errors.hcerr = hconferr;
% %errors.serr = serror;
% %errors.scerr = sconferr;
%
% [hdata2, hlab2] = formatData(spfeatures, spvlabels);
% hconf2 = test_boosted_dt_mc(hclassifier, hdata2);
% hconf2 = 1 ./ (1+exp(-hconf2));
% hconf2 = hconf2 ./ repmat(sum(hconf2, 2), 1, size(hconf2, 2));
%
% lab = (vlab==1)*1 + (vlab==2)*(hlab2+1) + (vlab==3)*7;
% pEst = [vconf(:, 1) repmat(vconf(:, 2), 1, 5).*hconf2 vconf(:, 3)];
% pTrue = repmat((1:7), numel(lab), 1).*repmat(lab, 1, 7);
%
% size(pEst)
% size(pTrue)
%
% [pMean, pVal, estBias] = evaluateProbabilityEstimate(pEst(:), pTrue(:), 0.05)
% pErrors.pMean = pMean;
% pErrors.pVal = pVal;
% pErrors.estBias = estBias;
%
%
%

%
%
%
% pVal = [step/2:step:(1-step/2)];
%
% nsteps = numel(pVal);
%
% estBias = 0;
%
% for v = 1:numel(pVal)
%     v1 = pVal(v)-step/2;
%     v2 = min(v1+step, 1);
%
%     ind = find((pEst >= v1) & (pEst <= v2));
%
%     pMean(v) = mean(pTrue(ind));
%     estBias = estBias + numel(ind)/numel(pEst)*(pMean(v)-pVal(v));
% end

author: zouchuhang

repo: LayoutNet