text:
%% GENERIC object CLASS (objectDefinition.m) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This class is designed to define a generic object and import its
% variables into the simulation space for the purpose of multi-vehicle
% control simulation.

% Author: James A. Douthwaite 20/05/2016

classdef objectDefinition < handle
    %% OBJECT BASE CLASS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % This class contains the basic properties of a generic agent, neither
    % aerial or ground based.
    properties
        % SIMULATION IDENTIFIERS
        objectID = 0;
        name;
        radius = 1;                 % Characteristic radius
        % KINEMATIC PROPERTIES
        localState  = zeros(6,1);   % 6DOF state vector
        % GEOMETRIC PROPERTIES
        GEOMETRY = struct(...
            'vertices',[],...
            'faces',[],...
            'normals',[],...
            'centroid',zeros(3,1));
    end
    % Private (OpenMAS-facing) properties
    properties (Access = private)
        % GLOBAL PROPERTIES (Virtual (SIMULATION) data container)
        GLOBAL;
    end
    %% ///////////////////////// MAIN METHODS /////////////////////////////
    methods (Access = public)
        % Constructor
        function [this] = objectDefinition(varargin)
            % This function is designed to automatically apply a generic
            % naming convention to all objects generated prior to a
            % simulation.
            % INPUTS:
            % varargin - A cell array of property/value pairs
            % OUTPUTS:
            % this     - The generated object

            % Check all system paths are available
            this.GetPathDependencies();

            % Assign the objectID (if required)
            if this.objectID == 0                                           % Generate objectID if not allocated by a super class
                this.objectID = this.CreateobjectID();                       % Assign the number as an ID tag
                spawnNewobject = 1;
            else
                spawnNewobject = 0;
            end

            % Assign a name string
            if numel(this.name) < 1
                this.name = this.CreateName(this.objectID);                   % No input name string specified, use greek naming scheme
            end

            % Assign default
            this.GLOBAL   = this.CreateGLOBAL();    % The GLOBAL structure
            this.GEOMETRY = this.CreateGEOMETRY();  % Get the geometry of the object

            % //////////////// Check for user overrides ///////////////////
            [this] = this.ApplyUserOverrides(varargin); % Recursive overrides
            % /////////////////////////////////////////////////////////////

            % ////////// DISPLAY THE CURRENT object PROPERTIES ////////////
            if spawnNewobject
                fprintf('objectcount: %d\tname: %s\ttype: %s\n',this.objectID,this.name,class(this));
            end
        end
        % Setup (default 6-dof state)
        function [this] = setup(this,localXYZvelocity,localXYZrotations)   % [x y z phi theta psi]
            % This function is called in order to build the initial state
            % vector for the generic agent class 'objectDefinition'.
            % INPUTS:

            % ASSUMPTIONS:
            % - The designed state is described in the ENU axes.
            % - The state is of the form:
            % - [x y z phi theta psi]

            % Standard state vector is the 6DOF: x_k = [x y z phi theta psi]
            this.localState = [zeros(3,1);localXYZrotations];
        end
        % Main  (default 6-dof state)
        function [this] = main(this,ENV,varargin)
            % This is a generic process cycle of an object that accepts no
            % input commands/feedback and simply updates its states based
            % on its current attributes
            % INPUT:
            % varargin - Generic variable input container, containing the simulation timestep
            % OUTPUT:
            % this - The updated object

            % DETERMINE THE TIMESTEP
            if isstruct(ENV)
                dt = ENV.dt;
            else
                error('object TIME packet is invalid.');
            end

            % The state differential
            x_dot = [0;0;0;0;0;0];

            % MAINTAIN LOCAL INPUT CONDITIONS
            % Default behaviour of a system is to move with constant input
            % conditions.
            [dXdt]   = this.SingleIntegratorDynamics(this.localState,x_dot);
            newState = this.localState + dt*dXdt;

            % UPDATE THE CLASS GLOBAL PROPERTIES
            this = this.GlobalUpdate(dt,newState);
        end
    end

    %% /////////////////////// CHECK FUNCTIONS ////////////////////////////
    methods
        % Get the 'idle' status
        function [flag] = IsIdle(this)
            flag = this.GetGLOBAL('idleStatus');
        end
        % Get the dimensionality of the object
        function [flag] = Is3D(this)
            flag = this.GetGLOBAL('is3D');
        end
    end
    %% ////////////////////// GET/SET FUNCTIONS ///////////////////////////
    methods
        % Set the radius
        function set.radius(this,radius)
            this.radius = radius;               % Set the radius
            this.SetGLOBAL('radius',radius);    % Set the simulation radius
            % If agent has GEOMETRY
            if ~isempty(this.GEOMETRY.vertices)
                % Normalise the geometry
                this.GEOMETRY = OMAS_graphics.normalise(this.GEOMETRY);
                % Scale the geometry to the radius
                this.GEOMETRY = OMAS_graphics.scale(this.GEOMETRY,this.radius);
            end
        end
        % Get the GLOBAL structure
        function [param] = GetGLOBAL(this,label)
            % By default return the GLOBAL structure
            if nargin < 2
                param = this.GLOBAL;
                return
            end
            assert(ischar(label),"Expecting a string parameter label");
            % Return the associated parameter
            param = this.GLOBAL.(label);
        end
        % Set the virtual structure
        function [this]  = SetGLOBAL(this,label,value)
            % By default set the GLOBAL structure
            if nargin < 3
                % Input sanity check
                assert(isstruct(label),"The object's GLOBAL field must be a structure.");
                % Set the GLOBAL structure
                this.GLOBAL = label;
                return
            end
            assert(ischar(label),"Expecting a string parameter label");
            % Assign the virtual structure
            this.GLOBAL.(label) = value;
        end
    end
    %% ////////////////// BASIC PROPERTY ASSIGNMENTS //////////////////////
    methods
        % PARSE A GENERIC INPUT SET AGAINST A DEFAULT CONFIG STRUCTURE
        function [this]     = ApplyUserOverrides(this,pairArray)
            % This function is designed to parse a generic set of user
            % inputs and allow them to be compared to a default input
            % structure. This should be called in the class constructor

            % Input sanity check #1
            if nargin < 2 || numel(pairArray) == 0
                return
            end
            % Call the all parameter overrider
            [this] = GetParameterOverrides_recursive(this,pairArray);
        end
        % Get the object GEOMETRY structure
        function [GEOMETRY] = CreateGEOMETRY(this)
            % This function prepares object STL files for representation.
            % ASSUMPTIONS:
            % - The object stl file has the same name as the object being simulated.
            % - The STL is correctly rotated to match the axes of XYZ of the agents
            %   local frame.

            % This functions allow additional operations to the GEOMETRY
            % structure to happen other than sourcing the basic structure.

            % Call the generic object.stl geometry parser
            [GEOMETRY] = this.GetObjectGeometry(this);
        end
    end
    % Only the objectDefinition class has access
    methods (Static, Access = private)
        % Get the object GLOBAL structure
        function [GLOBAL]   = CreateGLOBAL()
            % This is function that assembles the template desciption of an
            % object

            % Define the GLOBAL structure
            GLOBAL = struct();
            GLOBAL.type = OMAS_objectType.misc;
            GLOBAL.hitBoxType = OMAS_hitBoxType.none;
            GLOBAL.radius = 0.5;                        % Diameter of 1m
            GLOBAL.colour = rand(1,3,'single');         % Virtual colour (for plotting)
            GLOBAL.symbol = 'square';                   % Representative symbol
            GLOBAL.position = [0;0;0];                  % Global Cartesian position
            GLOBAL.velocity = [0;0;0];                  % Global Cartesian velocity
            GLOBAL.quaternion = [1;0;0;0];              % Global quaternion pose
            GLOBAL.idleStatus = true;                   % object idle logical
            GLOBAL.is3D = true;                         % object operates in 3D logical
            GLOBAL.priorState = [];
        end
        % Get the object name
        function [namestr]  = CreateName(objectID)
            % No input name string specified, use greek naming scheme
            defaultID = { 'alpha', 'beta','gamma', 'delta','epsilon'  ...
                        ,  'zeta',  'eta','theta',  'iota','kappa'    ...
                        ,'lambda',   'mu',   'nu',    'xi','omicron'  ...
                        ,    'pi',  'rho','sigma',   'tau','upsilon'  ...
                        ,   'phi',  'chi',  'psi', 'omega'};
            % GENERATE DEFAULT NAME STRING
            defLength = length(defaultID);
            cycle = floor(double(objectID)/defLength);
            IDvalue = mod(objectID,defLength);
            if IDvalue == 0
                IDvalue = 24;
            end
            namestr = sprintf('%s%03d',defaultID{IDvalue},cycle+1);
        end
        % Get the object ID
        function [objectID] = CreateobjectID()
            % This function allocates a basic ID convention to the class
            % heirarchy to provide object automatic numeration.

            persistent objectCount;

            % DEFINE ID BASED ON EXISTING objectS
            if isempty(objectCount)
                % INITIAL object
                objectCount = 1;
            else
                % objectS EXIST
                objectCount = objectCount + 1;
            end
            % ALLOCATE object ID
            objectID = uint16(objectCount);
        end
    end

    %% STATE UPDATE METHODS
    methods (Static)
        % STATE UPDATE - TWO WHEELED DIFFERENTIAL ROBOT
        function [dXdt] = DifferencialDriveDynamics(X,d,a_left,a_right)
            % Becker, M. (2006). Obstacle avoidance procedure for mobile robots. ABCM Symposium Series
            % [x;y;theta;v;omega]

            % BUILD THE STATE DIFFERENTAL
            dXdt = zeros(5,1);
            dXdt(1) = X(4)*cos(X(3));             % x velocity
            dXdt(2) = X(4)*sin(X(3));             % y velocity
            dXdt(3) = X(5);                       % rotational rate
            dXdt(4) = (a_right + a_left)/2;       % Linear acceleration
            dXdt(5) = (a_right - a_left)/d;       % Angular acceleration
        end
        % STATE UPDATE - UNI-CYCLE ROBOT
        function [dXdt] = UnicycleDynamics(X,speed,headingRate)
            % This function provides the dynamics of a uni-cycle model
            % moving in on a 2D plane.

            % THE STATE UPDATE
            dXdt = zeros(3,1);
            dXdt(1) = speed*cos(X(3));
            dXdt(2) = speed*sin(X(3));
            dXdt(3) = headingRate;
        end
        % STATE UPDATE - DOUBLE INTEGRATOR
        function [dXdt] = DoubleIntegratorDynamics(X,U)
            % The relationship between the input and state vector is by
            % definition of a double integrator: dX(dX/dt)/dt = U

            % X(t) here is assumed to be of the form X(t) = [q;dq]

            qDim = length(X)/2;

            % dXdt is of the form dX(t) = [dq;ddq];
            dXdt = [X(qDim+1:2*qDim,1);     % Velocity states
                    U];                     % Acceleration states
        end
        % STATE UPDATE - SINGLE INTEGRATOR
        function [dXdt] = SingleIntegratorDynamics(X,U)
            % The relationship between the input and state vector is by
            % definition of a single integrator: dX/dt = U

            % Sanity check
            assert(length(X) == length(U),"Expecting the number of inputs to correspond to the state vector.");

            % X(t) here is assumed to be of the form X(t) = [q]

            % The change in state is literally the input
            dXdt = U;
        end
        % STATE UPDATE - SIMPLE EULER VELOCITIES
        function [dXdt] = SimpleDynamics(X,velocity_k_plus,omega_k_plus)
            % This function assumes that the velocity changes are
            % implemented this timestep directly, integration then occurs
            % including the updates from this timestep:
            % State vector is assumed:
            % 2D - [x y psi dx dy dpsi]'
            % 3D - [x y z phi theta psi dx dy dz dphi dtheta dpsi]'

            % CHECK EXPECTED INPUT DIMENSIONS
            Aflag = numel(velocity_k_plus) == 3 && numel(omega_k_plus) == 3; % IS 3D
            Bflag = numel(velocity_k_plus) == 2 && numel(omega_k_plus) == 1; % IS 2D
            assert((Aflag || Bflag) == 1,'Integrator inputs are of the wrong dimensions');

            % Either rotation occurs in one dimension (2D) or in three (3D)
            rotationDOF = numel(omega_k_plus); % Is heading a yaw (2D) or [phi,theta psi] (3D)
            if rotationDOF == 1
                linearIndices = 1:2;
                angularIndices = 3;
            else
                linearIndices = 1:3;
                angularIndices = 4:6;
            end

            % THE STATE DIFFERENTIAL
            dXdt = zeros(numel([linearIndices,angularIndices]),1);
            dXdt(linearIndices) = velocity_k_plus;
            dXdt(angularIndices) = omega_k_plus;
        end
    end
    %% GLOBAL REPRESENTATION UPDATE METHODS
    methods
        % INITIAL 3D STATE - [x y z phi theta psi dx dy dz dphi dtheta dpsi]
        function [this] = setup_3DVelocities(this,localXYZVelocity,localXYZrotations)
            % This function calculates the intial state for a generic
            % object.
            % The default state vector:
            % [x y z phi theta psi; dx dy dz dphi dtheta dpsi]

            % BUILD STATE VECTOR
            this.localState = zeros(12,1);
            this.localState(4:5) = localXYZrotations(1:2);
            %this.localState(4:6) = localXYZrotations;       % Represent yaw relative to initial position
            this.localState(7:9) = localXYZVelocity;
            % Define the
            this.SetGLOBAL('priorState',this.localState);
        end
        % GLOBAL UPDATE - STATE VECTOR DEFINED AS: [x_t;x_dot]'
        function [this] = GlobalUpdate_3DVelocities(this,dt,eulerState)
            % This function preforms the state update for a state vector
            % defined as [x y z phi theta psi dx dy dz dphi dtheta dpsi].

            positionIndices = 1:3;
            angleIndices = 4:6;

            % USE THE 'RATE' STATES DIRECTLY
            localLinearRates  = eulerState(positionIndices+6);
            localAngularRates = eulerState(angleIndices+6);
            % Get current properties
            p_k = this.GetGLOBAL('position');
            v_k = this.GetGLOBAL('velocity');
            q_k = this.GetGLOBAL('quaternion');

            % MAP THE LOCAL RATES TO GLOBAL RATES AND INTEGRATE QUATERNION
            % PREVIOUS QUATERNION % [ IF THE QUATERNION DEFINES GB ]

            % PREVIOUS ROTATION-MATRIX
            R_k_plus = quat2rotm(q_k');
            % THE GLOBAL AXIS RATES
            omega = R_k_plus'*localAngularRates;
            % UPDATE THE QUATERNION POSE
            q_k_plus = OMAS_geometry.integrateQuaternion(q_k,omega,dt);
            % REDEFINE THE ROTATION MATRIX
            R_k_plus = quat2rotm(q_k_plus');
            % MAP THE LOCAL VELOCITY TO THE GLOBAL AXES
            v_k_plus = R_k_plus'*localLinearRates;
            p_k_plus = p_k + dt*v_k;

            % ///////////////// REASSIGN K+1 PARAMETERS ///////////////////
            this = this.GlobalUpdate_direct(...
                p_k_plus,...    % Global position at k plius
                v_k_plus,...    % Global velocity at k plus
                q_k_plus);      % Quaternion at k plus
        end
        % GLOBAL UPDATE - EULER 6DOF(3DOF) [x y z phi theta psi],[x y psi]
        function [this] = GlobalUpdate(this,dt,eulerState)
            % This function computes the position, velocity and
            % quaternion in the global ENU frame from a 6DOF state vector
            % with euler rotations.

            % This function is intended for state vectors of the form:
            % - [x y z phi theta psi]
            % - [x y phi]

            % NOTATION INDICIES
            if this.Is3D
                positionIndices = 1:3;
                eulerIndices = 4:6;     % The rotations about X,Y,Z
            else
                positionIndices = 1:2;
                eulerIndices = 3;       % The rotations about Z
            end

            % Get current properties
            p_k = this.GetGLOBAL('position');
            v_k = this.GetGLOBAL('velocity');
            q_k = this.GetGLOBAL('quaternion');
            x_k = this.GetGLOBAL('priorState');

            % Populate the prior state if necessary
            if isempty(x_k)
                x_k = this.localState;
            end

            % Equivalent rates
            localLinearRates  = (eulerState(positionIndices) - x_k(positionIndices))/dt;
            localAngularRates = (eulerState(eulerIndices) - x_k(eulerIndices))/dt;

            % EULER ROTATIONS -> GLOBAL ROTATIONS
            if numel(localAngularRates) ~= 3
                localAngularRates = [0;0;1]*localAngularRates;
                localLinearRates  = [localLinearRates;0];
            end

            % MAP THE LOCAL RATES TO GLOBAL RATES AND INTEGRATE QUATERNION
            % PREVIOUS ROTATION-MATRIX
            R_k_plus = quat2rotm(q_k');
            % THE GLOBAL AXIS RATES
            omega = R_k_plus'*localAngularRates;
            % UPDATE THE QUATERNION POSE
            q_k_plus = OMAS_geometry.integrateQuaternion(q_k,omega,dt);
            % REDEFINE THE ROTATION MATRIX
            R_k_plus = quat2rotm(q_k_plus');
            % MAP THE LOCAL VELOCITY TO THE GLOBAL AXES
            velocity_k_plus = R_k_plus'*localLinearRates;
            position_k_plus = p_k + dt*v_k;
            % ///////////////// REASSIGN K+1 PARAMETERS ///////////////////
            this = this.GlobalUpdate_direct(...
                position_k_plus,... 	% Global position at k plius
                velocity_k_plus,...     % Global velocity at k plus
                q_k_plus);              % Quaternion at k plus
        end
        % GLOBAL UPDATE - EULER 6DOF(3DOF) (NO ROTATIONS)
        function [this] = GlobalUpdate_fixedFrame(this,dt,eulerState)
            % This function computes the position and velocity whilst
            % maintaining a fixed quaternion rotation (constant reference
            % orientation) in the global ENU frame from a 6DOF state vector
            % with euler rotations.

            % NOTATION INDICIES
            if this.Is3D
                positionIndices = 1:3;
            else
                positionIndices = 1:2;
            end

            % Get current properties
            position_k   = this.GetGLOBAL('position');
            velocity_k   = this.GetGLOBAL('velocity');
            quaternion_k = this.GetGLOBAL('quaternion');
            state_k      = this.GetGLOBAL('priorState');

            % Equivalent velocity
            velocity_k_plus = (eulerState(positionIndices) - state_k(positionIndices))/dt;

            % ROTATION RATES ABOUT THE GLOBAL AXES
            if numel(eulerState) ~= 6
                velocity_k_plus = [velocity_k_plus;0];
            end

            % NEW ROTATION MATRIX (G>B)
            R_k_plus = OMAS_geometry.quaternionToRotationMatrix(quaternion_k);
            % MAP THE LOCAL VELOCITY TO THE GLOBAL AXES
            velocity_k_plus = R_k_plus'*velocity_k_plus;
            position_k_plus = position_k + dt*velocity_k;

            % ///////////////// REASSIGN K+1 PARAMETERS ///////////////////
            [this] = this.GlobalUpdate_direct(...
                position_k_plus,...     % Global position at k plius
                velocity_k_plus,...     % Global velocity at k plus
                quaternion_k_plus,...   % Quaternion at k plus
                eulerState);            % The new state for reference
        end
        % GLOBAL UPDATE - DIRECTLY (global default)
        function [this] = GlobalUpdate_direct(this,p,v,q)
            % Under this notation, the state vector already contains the
            % global parameters of the object.
            % INPUTS:
            % globalPosition - 3D global cartesian position given in the ENU coordinate frame.
            % globalVelocity - 3D global cartesian velocity given in the ENU coordinate frame.
            % quaternion     - The new quaternion pose of body in the global frame.
            % R              - The rotation of the body
            % this.localState - The previous localState (independant of convention)
            % eulerState     - The new state as reported by the agent

            % Input sanity check
            assert(IsColumn(p,3),'Global position must be a 3D column vector [3x1].');
            assert(IsColumn(v,3),'Global velocity must be a 3D column vector [3x1].');
            assert(IsColumn(q,4),'Global pose must be a 4D quaternion vector [4x1].');
            assert(size(this.localState,2) == 1,'The length of the objects state update must match the its local state.');

            % ///////////////// REASSIGN K+1 PARAMETERS ///////////////////
            % Convert the quaternion to the equivalent rotation matrix
            R = OMAS_geometry.quaternionToRotationMatrix(q);
            % Assign the global parameters
            this.SetGLOBAL('position',p);                 	% Reassign the global position
            this.SetGLOBAL('velocity',v);                 	% Reassign the global velocity
            this.SetGLOBAL('quaternion',q);                	% Reassign the quaternion
            this.SetGLOBAL('R',R);                          % New rotation matrix
            this.SetGLOBAL('priorState',this.localState);  	% Record the previous state
        end
    end
    methods (Static)
        % Generic .STL to geometry parser
        function [GEOMETRY] = GetObjectGeometry(entity)
            % This function prepares object STL files for representation.
            % ASSUMPTIONS:
            % - The object stl file has the same name as the object being simulated.
            % - The STL is correctly rotated to match the axes of XYZ of the agents
            %   local frame.

            % ABSOLUTE PATH TO THE object STL LOCATION
            absPath = mfilename('fullpath');
            absPath = strrep(absPath,'objectDefinition','');               % Get the current m-files directory

            % GET THE AGENT SUPERCLASSES
            aliases = superclasses(entity);
            aliases = vertcat(class(entity),aliases);
            % MOVE THROUGH THE CLASS HEIRARCHY
            for aliasNo = 1:numel(aliases)
                % BUILD THE STL ASSOCIATION
                filename = strcat(absPath,char(aliases(aliasNo)),'.stl');  % Build the associated STL path
                % GET THE object PATCH FROM STL FILE (in the object directory)
                [GEOMETRY,getFlag] = OMAS_graphics.importStlFromFile(filename);
                if getFlag
                    break                                                  % Successful import
                end
            end

            % IF PATCH RETURNED, PROCESS FOR SIMULATION
            if isstruct(GEOMETRY)
                % Get the entities representative radius
                entityRadius = entity.GetGLOBAL('radius');
                % NORMALISE THE IMPORTED GEOMETRY
                GEOMETRY = OMAS_graphics.normalise(GEOMETRY);              % Normalise
                GEOMETRY = OMAS_graphics.scale(GEOMETRY,entityRadius);     % Scale
                GEOMETRY.normals = OMAS_graphics.normals(GEOMETRY);
                % ADD CENTROID
                GEOMETRY.centroid = zeros(1,3);                            % Assume vertices are relative to a centroid
            else
                % GEOMETRIC PROPERTIES
                GEOMETRY = struct('vertices',[],...
                                  'faces',[],...
                                  'normals',[],...
                                  'centroid',zeros(1,3));                  % If the object is something other than a point.
            end
            % Return patch, empty if not successful
        end
        % Dependancy check for all objects
        function GetPathDependencies()
            % This program preforms a check of the dependencies necessary
            % for the instantiation of object classes.

            % Get the path to the install directory
            repoString = mfilename('fullpath');
            repoString = repoString(1:(strfind(repoString,'objects')-1));
            % Add the path to the repo's environment
            addpath([repoString,'environment']);
            % Check all the OMAS_dependencies
            if ~OMAS_system.GetFileDependancies()
                error('Unable to get OpenMAS dependancies.');
            end
        end
    end
end


text:
function [avg_runtime, plug_events] = get_average_runtime_from_plug_data(appliance_name, training_days, setup)

    appliance_id = getApplianceID(appliance_name);
    global caching;
    if caching == 1
        if exist('cache_plugs_training.mat') == 2
            load('cache_plugs_training');
        else
            appliance_consumption_training = read_plug_data(setup.dataset, setup.household, appliance_id, training_days, setup.granularity);
            save('cache_plugs_training', 'appliance_consumption_training');
        end
    else
        appliance_consumption_training = read_plug_data(setup.dataset, setup.household, appliance_id, training_days, setup.granularity);
    end

    if strcmp(appliance_name, 'TV') == 1 ...
            || strcmp(appliance_name, 'Stereo') == 1 ...
            || strcmp(appliance_name, 'Fridge') == 1 ...
            || strcmp(appliance_name, 'Freezer') == 1

        on_events = [];
        off_events = [];
        cons = appliance_consumption_training;
        for i = 1:length(cons)-1
            % on events
            if cons(i) < 20 && cons(i+1) >= 20
                on_events = [on_events; i, i+1, 1, 1];
            elseif cons(i) >= 20 && cons(i+1) < 20
                if length(on_events) == 0
                    continue;
                end
                off_events = [off_events; i, i+1, -1, 1];
            end
        end

        if length(off_events) > length(on_events)
            off_events(end) = [];
        end

        avg_runtime = mean(off_events(:,1) - on_events(:,1));
        plug_events = {off_events, on_events};

    else
        error('Not implemented for appliance %s', appliance_name);
    end
end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function [D,C,K] = spm_dcm_KL(M)
% Computes the distance between two models based on prior responses
% FORMAT [D,C,K] = spm_dcm_KL(Mi,Mj)
%
% M{1:n}   - structure array of models
%
% D(n x n) - distance matrix (KL divergence)
% C{1:n}   - response covariances
% K{1:n}   - response means
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_KL.m 5219 2013-01-29 17:07:07Z spm $

% Volterra kernels
%==========================================================================

% time bins (if not specified)
%--------------------------------------------------------------------------
try
    dt = M{1}.dt;
    N  = M{1}.N;
catch

    % Bilinear representation
    %----------------------------------------------------------------------
    M0 = spm_bireduce(M{1},P{1});
    s  = real(eig(full(M0)));
    s  = max(s(find(s < 0)));
    N  = 32;
    dt = -4/(s*N);

end

% get covariances of prior responses
%==========================================================================
m     = length(M);
for i = 1:m

    % Get parameters (adding a little to prevent expansion around zero)
    %----------------------------------------------------------------------
    P      = M{i}.pE;
    pC     = M{i}.pC;
    P      = spm_vec(P) + sqrt(diag(pC))/8;
    P      = spm_unvec(P,M{i}.pE);

    % add a little to inputs (to prevent expansion around zero)
    %----------------------------------------------------------------------
    M{i}.u = ones(M{i}.m,1)/8;

    % get eigen-space of parameters for computational efficiency
    %----------------------------------------------------------------------
    V      = spm_svd(pC);
    pC     = V'*pC*V;

    % get derivative of kernels w.r.t. parameters
    %----------------------------------------------------------------------
    [dkdp,k] = spm_diff('spm_kernel',M{i},P,N,dt,2,{[],V});

    % prior mean and covariance of kernels
    %----------------------------------------------------------------------
    k     = spm_vec(k);
    dk    = sparse(length(k),length(dkdp));
    for j = 1:length(dkdp)
        dk(:,j) = spm_vec(dkdp{j});
    end
    K{i} = k;
    C{i} = dk*pC*dk';
    n(i) = length(pC);

end

% evaluate KL divergence
%--------------------------------------------------------------------------
for i = 1:m
    for j = 1:m
        Pi     = spm_pinv(C{i});
        k      = K{i} - K{j};
        d      = spm_logdet(C{i}) - spm_logdet(C{j}) + ...
                 trace(Pi*C{j}) + k'*Pi*k - n(i);
        D(i,j) = d/2;
    end
end

% ensure D is symmetric
%--------------------------------------------------------------------------
D   = (D + D')/2;

author: spm


text:
function varargout = lvd_EditCoordSysOriginRefFrameGUI(varargin)
% LVD_EDITCOORDSYSORIGINREFFRAMEGUI MATLAB code for lvd_EditCoordSysOriginRefFrameGUI.fig
%      LVD_EDITCOORDSYSORIGINREFFRAMEGUI, by itself, creates a new LVD_EDITCOORDSYSORIGINREFFRAMEGUI or raises the existing
%      singleton*.
%
%      H = LVD_EDITCOORDSYSORIGINREFFRAMEGUI returns the handle to a new LVD_EDITCOORDSYSORIGINREFFRAMEGUI or the handle to
%      the existing singleton*.
%
%      LVD_EDITCOORDSYSORIGINREFFRAMEGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in LVD_EDITCOORDSYSORIGINREFFRAMEGUI.M with the given input arguments.
%
%      LVD_EDITCOORDSYSORIGINREFFRAMEGUI('Property','Value',...) creates a new LVD_EDITCOORDSYSORIGINREFFRAMEGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before lvd_EditCoordSysOriginRefFrameGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to lvd_EditCoordSysOriginRefFrameGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help lvd_EditCoordSysOriginRefFrameGUI

% Last Modified by GUIDE v2.5 12-Jan-2021 18:10:17

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @lvd_EditCoordSysOriginRefFrameGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @lvd_EditCoordSysOriginRefFrameGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before lvd_EditCoordSysOriginRefFrameGUI is made visible.
function lvd_EditCoordSysOriginRefFrameGUI_OpeningFcn(hObject, eventdata, handles, varargin)
    % This function has no output args, see OutputFcn.
    % hObject    handle to figure
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    % varargin   command line arguments to lvd_EditCoordSysOriginRefFrameGUI (see VARARGIN)

    % Choose default command line output for lvd_EditCoordSysOriginRefFrameGUI
    handles.output = hObject;

    centerUIFigure(hObject);

    refFrame = varargin{1};
    setappdata(hObject, 'refFrame', refFrame);

    lvdData = varargin{2};
    setappdata(hObject, 'lvdData', lvdData);

	populateGUI(handles, lvdData, refFrame);

    % Update handles structure
    guidata(hObject, handles);

    % UIWAIT makes lvd_EditCoordSysOriginRefFrameGUI wait for user response (see UIRESUME)
    uiwait(handles.lvd_EditCoordSysOriginRefFrameGUI);

function populateGUI(handles, lvdData, refFrame)
    set(handles.refFrameNameText,'String',refFrame.getName());

    pointsListBoxStr = lvdData.geometry.points.getListboxStr();
    handles.frameOriginCombo.String = pointsListBoxStr;
    handles.frameOriginCombo.Value = lvdData.geometry.points.getIndsForPoints(refFrame.origin);

    coordSysListBoxStr = lvdData.geometry.coordSyses.getListboxStr();
    handles.coordSysCombo.String = coordSysListBoxStr;
    handles.coordSysCombo.Value = lvdData.geometry.coordSyses.getIndsForCoordSyses(refFrame.coordSys);

    handles.xColorCombo.String = ColorSpecEnum.getListboxStr();
    handles.xColorCombo.Value = ColorSpecEnum.getIndForName(refFrame.xAxisColor.name);

    handles.yColorCombo.String = ColorSpecEnum.getListboxStr();
    handles.yColorCombo.Value = ColorSpecEnum.getIndForName(refFrame.yAxisColor.name);

    handles.zColorCombo.String = ColorSpecEnum.getListboxStr();
    handles.zColorCombo.Value = ColorSpecEnum.getIndForName(refFrame.zAxisColor.name);

    handles.xLineStyleCombo.String = LineSpecEnum.getListboxStr();
    handles.xLineStyleCombo.Value = LineSpecEnum.getIndForName(refFrame.xAxisLineSpec.name);

    handles.yLineStyleCombo.String = LineSpecEnum.getListboxStr();
    handles.yLineStyleCombo.Value = LineSpecEnum.getIndForName(refFrame.yAxisLineSpec.name);

    handles.zLineStyleCombo.String = LineSpecEnum.getListboxStr();
    handles.zLineStyleCombo.Value = LineSpecEnum.getIndForName(refFrame.zAxisLineSpec.name);

    handles.xLineWidthText.String = fullAccNum2Str(refFrame.xAxisLineWidth);
    handles.yLineWidthText.String = fullAccNum2Str(refFrame.yAxisLineWidth);
    handles.zLineWidthText.String = fullAccNum2Str(refFrame.zAxisLineWidth);

    handles.scaleFactorText.String = fullAccNum2Str(refFrame.scaleFactor);

% --- Outputs from this function are returned to the command line.
function varargout = lvd_EditCoordSysOriginRefFrameGUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
    if(isempty(handles))
        varargout{1} = false;
    else
        lvdData = getappdata(hObject, 'lvdData');
        refFrame = getappdata(hObject, 'refFrame');

        refFrame.setName(handles.refFrameNameText.String);

        refFrame.origin = lvdData.geometry.points.getPointsForInds(handles.frameOriginCombo.Value);
        refFrame.coordSys = lvdData.geometry.coordSyses.getCoordSysesForInds(handles.coordSysCombo.Value);

        str = handles.xColorCombo.String{handles.xColorCombo.Value};
        refFrame.xAxisColor = ColorSpecEnum.getEnumForListboxStr(str);

        str = handles.xColorCombo.String{handles.yColorCombo.Value};
        refFrame.yAxisColor = ColorSpecEnum.getEnumForListboxStr(str);

        str = handles.xColorCombo.String{handles.zColorCombo.Value};
        refFrame.zAxisColor = ColorSpecEnum.getEnumForListboxStr(str);

        str = handles.xLineStyleCombo.String{handles.xLineStyleCombo.Value};
        refFrame.xAxisLineSpec = LineSpecEnum.getEnumForListboxStr(str);

        str = handles.xLineStyleCombo.String{handles.yLineStyleCombo.Value};
        refFrame.yAxisLineSpec = LineSpecEnum.getEnumForListboxStr(str);

        str = handles.xLineStyleCombo.String{handles.zLineStyleCombo.Value};
        refFrame.zAxisLineSpec = LineSpecEnum.getEnumForListboxStr(str);

        refFrame.xAxisLineWidth = str2double(handles.xLineWidthText.String);
        refFrame.yAxisLineWidth = str2double(handles.yLineWidthText.String);
        refFrame.zAxisLineWidth = str2double(handles.zLineWidthText.String);

        refFrame.scaleFactor = str2double(handles.scaleFactorText.String);

        varargout{1} = true;
        close(handles.lvd_EditCoordSysOriginRefFrameGUI);
    end

% --- Executes on button press in saveAndCloseButton.
function saveAndCloseButton_Callback(hObject, eventdata, handles)
% hObject    handle to saveAndCloseButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    errMsg = validateInputs(handles);

    if(isempty(errMsg))
        uiresume(handles.lvd_EditCoordSysOriginRefFrameGUI);
    else
        msgbox(errMsg,'Invalid Reference Frame Inputs','error');
    end

function errMsg = validateInputs(handles)
    errMsg = {};

    if(isempty(strtrim(handles.refFrameNameText.String)))
        errMsg{end+1} = 'Reference frame name must contain more than white space and must not be empty.';
    end

    val = str2double(get(handles.xLineWidthText,'String'));
    enteredStr = get(handles.xLineWidthText,'String');
    numberName = 'X Axis Line Width';
    lb = 0.01;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(val, numberName, lb, ub, isInt, errMsg, enteredStr);

    val = str2double(get(handles.yLineWidthText,'String'));
    enteredStr = get(handles.yLineWidthText,'String');
    numberName = 'Y Axis Line Width';
    lb = 0.01;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(val, numberName, lb, ub, isInt, errMsg, enteredStr);

    val = str2double(get(handles.zLineWidthText,'String'));
    enteredStr = get(handles.zLineWidthText,'String');
    numberName = 'Z Axis Line Width';
    lb = 0.01;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(val, numberName, lb, ub, isInt, errMsg, enteredStr);

    val = str2double(get(handles.scaleFactorText,'String'));
    enteredStr = get(handles.scaleFactorText,'String');
    numberName = 'Axes Scale Factor';
    lb = 0.001;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(val, numberName, lb, ub, isInt, errMsg, enteredStr);

% --- Executes on button press in cancelButton.
function cancelButton_Callback(hObject, eventdata, handles)
% hObject    handle to cancelButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    close(handles.lvd_EditCoordSysOriginRefFrameGUI);

function refFrameNameText_Callback(hObject, eventdata, handles)
% hObject    handle to refFrameNameText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of refFrameNameText as text
%        str2double(get(hObject,'String')) returns contents of refFrameNameText as a double

% --- Executes during object creation, after setting all properties.
function refFrameNameText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to refFrameNameText see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on key press with focus on lvd_EditCoordSysOriginRefFrameGUI or any of its controls.
function lvd_EditCoordSysOriginRefFrameGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_EditCoordSysOriginRefFrameGUI (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'enter'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'escape'
            close(handles.lvd_EditCoordSysOriginRefFrameGUI);
    end

% --- Executes on selection change in frameOriginCombo.
function frameOriginCombo_Callback(hObject, eventdata, handles)
% hObject    handle to frameOriginCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns frameOriginCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from frameOriginCombo

% --- Executes during object creation, after setting all properties.
function frameOriginCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to frameOriginCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in coordSysCombo.
function coordSysCombo_Callback(hObject, eventdata, handles)
% hObject    handle to coordSysCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns coordSysCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from coordSysCombo

% --- Executes during object creation, after setting all properties.
function coordSysCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to coordSysCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in xColorCombo.
function xColorCombo_Callback(hObject, eventdata, handles)
% hObject    handle to xColorCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns xColorCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from xColorCombo

% --- Executes during object creation, after setting all properties.
function xColorCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to xColorCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in xLineStyleCombo.
function xLineStyleCombo_Callback(hObject, eventdata, handles)
% hObject    handle to xLineStyleCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns xLineStyleCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from xLineStyleCombo

% --- Executes during object creation, after setting all properties.
function xLineStyleCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to xLineStyleCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function xLineWidthText_Callback(hObject, eventdata, handles)
% hObject    handle to xLineWidthText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of xLineWidthText as text
%        str2double(get(hObject,'String')) returns contents of xLineWidthText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

% --- Executes during object creation, after setting all properties.
function xLineWidthText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to xLineWidthText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in yColorCombo.
function yColorCombo_Callback(hObject, eventdata, handles)
% hObject    handle to yColorCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns yColorCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from yColorCombo

% --- Executes during object creation, after setting all properties.
function yColorCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to yColorCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in yLineStyleCombo.
function yLineStyleCombo_Callback(hObject, eventdata, handles)
% hObject    handle to yLineStyleCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns yLineStyleCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from yLineStyleCombo

% --- Executes during object creation, after setting all properties.
function yLineStyleCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to yLineStyleCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function yLineWidthText_Callback(hObject, eventdata, handles)
% hObject    handle to yLineWidthText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of yLineWidthText as text
%        str2double(get(hObject,'String')) returns contents of yLineWidthText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

% --- Executes during object creation, after setting all properties.
function yLineWidthText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to yLineWidthText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in zColorCombo.
function zColorCombo_Callback(hObject, eventdata, handles)
% hObject    handle to zColorCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns zColorCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from zColorCombo

% --- Executes during object creation, after setting all properties.
function zColorCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to zColorCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in zLineStyleCombo.
function zLineStyleCombo_Callback(hObject, eventdata, handles)
% hObject    handle to zLineStyleCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns zLineStyleCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from zLineStyleCombo

% --- Executes during object creation, after setting all properties.
function zLineStyleCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to zLineStyleCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function zLineWidthText_Callback(hObject, eventdata, handles)
% hObject    handle to zLineWidthText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of zLineWidthText as text
%        str2double(get(hObject,'String')) returns contents of zLineWidthText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

% --- Executes during object creation, after setting all properties.
function zLineWidthText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to zLineWidthText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function scaleFactorText_Callback(hObject, eventdata, handles)
% hObject    handle to scaleFactorText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of scaleFactorText as text
%        str2double(get(hObject,'String')) returns contents of scaleFactorText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

% --- Executes during object creation, after setting all properties.
function scaleFactorText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to scaleFactorText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [RefinedReactionsCarryingFlux, BlockedRefinedReactions] = testRefinedReactions(microbeID, BlockedRxns)
% Tests if reactions that were added through the refinement of gene
% annotations in the reconstructed organism can carry flux.
%
% INPUT
% microbeID         Microbe ID in data file with refined genome annotations
% BlockedRxns       Reactions in the COBRA model that cannot carry flux:
%                   output of identifyBlockedRxns function
%
% OUTPUT
% RefinedReactionsCarryingFlux     Reactions that were added through
%                                  refinement of the genome annotations
%                                  and can carry flux
%
% BlockedRefinedReactions          Reactions that were added through
%                                  refinement of the genome annotations
%                                  and cannot carry flux
%
%
% Almut Heinken, Nov 2019

BlockedRefinedReactions = {};
RefinedReactionsCarryingFlux = {};

if isfile('gapfilledGenomeAnnotation.txt')

    genomeAnnotation = readtable('gapfilledGenomeAnnotation.txt', 'ReadVariableNames', false, 'Delimiter', 'tab','TreatAsEmpty',['UND. -60001','UND. -2011','UND. -62011']);
    genomeAnnotation = table2cell(genomeAnnotation);

    % get the reactions that were added for the reconstruction
    findRxns=find(strcmp(microbeID,genomeAnnotation(:,1)));
    if ~isempty(findRxns)
        annRxns(:,1)=genomeAnnotation(findRxns(:,1),2);

        % find the overlap with blocked reactions in the model
        BlockedRefinedReactions = intersect(annRxns(:,1),BlockedRxns);
        RefinedReactionsCarryingFlux = setdiff(annRxns(:,1),BlockedRxns);
    end

end

end

author: opencobra

repo: cobratoolbox


text:
% Ot Matlab interface to OmniTRANS.
%
% Commands:
%    OtStart                     Initialize the Ot-interface.
%    OtStop                      Stop using the Ot-interface.
%
%    OtTableNew,OtQueryNew       Create a table/query handle.
%    OtTableFree,OtQueryFree     Release a table handle.
%
%    OtTableSetName              Set the path to a .db file.
%    OtQuerySetSql               Set an SQL statement for a query.
%
%    OtTableOpen,OtQueryOpen     Open an OmniTRANS table/SELECT-query.
%    OtTableClose,OtQueryClose   Close an OmniTRANS table/SELECT-query.
%    OtQueryExecute              Execute a non-SELECT SQL query.
%
%    OtTableFields,OtQueryFields             Get an array field names.
%    OtTableEmpty,OtQueryEmpty               Is the record set empty?
%    OtTableRecordCount,OtQueryRecordCount   Return the number of records.
%
%    OtTableFirst,OtQueryFirst   Move to the first record.
%    OtTableNext,OtQueryNext     Move to the next record.
%    OtTableEof,OtQueryEof       Check whether the last record was read.
%
%    OtTableLast,OtQueryLast     Move to the last record.
%    OtTablePrior,OtQueryPrior   Move to the previous record.
%    OtTableBof,OtQueryBof       Check whether the first record was read.
%
%    OtTableGet,OtQueryGet       Get contents of a field at the current record.
%    OtTableGetAll,OtQueryGetAll Get contents of all fields at the
%                                current record in a string array.
%
%    OtTableAppend               Put the table into append mode.
%    OtTableEdit                 Put the table into edit mode.
%    OtTableSet                  Modify the current record.
%    OtTableDelete               Delete the current record.
%    OtTablePost                 Post changes to the database table.
%
%    OtPmtu,OtPmturi                      Create a PMTU or PMTURI structure.
%    OtMatrixCubeNew,OtSkimCubeNew        Create a new matrix/skim cube handle.
%    OtMatrixCubeFree,OtSkimCubeFree      Release a matrix/skim cube handle.
%    OtMatrixCubeGet,OtSkimCubeGet        Obtain a matrix from a matrix/skim cube.
%    OtMatrixCubeSet,OtSkimCubeSet        Store a matrix in a matrix/skim cube.
%    OtMatrixCubeDelete,OtSkimCubeDelete  Delete a matrix in a matrix/skim cube.
%
% To be able to use this interface, ensure the OmniTRANS directory
% (with OtInterface.dll and all those BPL files) is on the system path.
%

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 2006-02-28 Created (Johan Meijdam, VORtech)
% 2010-01-22 Updated by Feike Brandt, OmniTRANS International B.V.

% Implementation: Ot.cpp

function [result] = Ot(command, arguments)
    error(['This is Ot.m. In case you read this error message,' ...
           'check whether there is an Ot.dll file in the same directory.']);

author: Sable

repo: mcbench-benchmarks


text:
function y = centroid(M, x)
% Attempts the computation of a centroid of a set of points on a manifold.
%
% function y = centroid(M, x)
%
% M is a structure representing a manifold.
% x is a cell of points on that manifold.

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

    % For now, just apply a few steps of gradient descent for Karcher means

    n = numel(x);

    problem.M = M;

    problem.cost = @cost;
    function val = cost(y)
        val = 0;
        for i = 1 : n
            val = val + M.dist(y, x{i})^2;
        end
        val = val/2;
    end

    problem.grad = @grad;
    function g = grad(y)
        g = M.zerovec(y);
        for i = 1 : n
            g = M.lincomb(y, 1, g, -1, M.log(y, x{i}));
        end
    end

    % This line can be uncommented to check that the gradient is indeed
    % correct. This should always be the case if the dist and the log
    % functions in the manifold are correct.
    % checkgradient(problem); pause;

    query = warning('query', 'manopt:getHessian:approx');
    warning('off', 'manopt:getHessian:approx');
    options.verbosity = 0;
    options.maxiter = 15;
    y = trustregions(problem, x{randi(n)}, options);
    warning(query.state, 'manopt:getHessian:approx');

end

author: NicolasBoumal

repo: manopt


text:
function varargout = help1(varargin)
% HELP1 M-file for help1.fig
%      HELP1, by itself, creates a new HELP1 or raises the existing
%      singleton*.
%
%      H = HELP1 returns the handle to a new HELP1 or the handle to
%      the existing singleton*.
%
%      HELP1('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in HELP1.M with the given input arguments.
%
%      HELP1('Property','Value',...) creates a new HELP1 or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before help1_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to help1_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help help1

% Last Modified by GUIDE v2.5 16-Jul-2010 19:35:59

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @help1_OpeningFcn, ...
                   'gui_OutputFcn',  @help1_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before help1 is made visible.
function help1_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to help1 (see VARARGIN)

% Choose default command line output for help1
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes help1 wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = help1_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes during object creation, after setting all properties.
function axes1_CreateFcn(hObject, eventdata, handles)
imshow(imread('telefono.jpg'))
% hObject    handle to axes1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate axes1


text:
function [varsigma] = ppcaPosteriorVar(model, X);

% PPCAPOSTERIORVAR Mean and variances of the posterior at points given by X.
% FORMAT
% DESC returns the posterior mean and variance for a given set of
% points.
% ARG model : the model for which the posterior will be computed.
% ARG x : the input positions for which the posterior will be
% computed.
% RETURN sigma : the variances of the posterior distributions.
%
% SEEALSO : ppcaCreate, ppcaPosteriorMeanVar
%
% COPYRIGHT : Neil D. Lawrence, 2008

% MLTOOLS

varsigma = repmat(1/model.beta, size(X, 1), 1);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ t_stop, y_stop ] = p38_stop ( neqn )

%*****************************************************************************80
%
%% P38_STOP returns the stopping point for problem p38.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_STOP, Y_STOP(NEQN), the final data.
%
  delta = p38_param ( 'GET', 'DELTA' );
  t_stop = 2.0 / delta;
  y_stop = zeros ( neqn, 1 );
  y_stop(1) = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%                           benchmarkSynthetic.m
%
% A general framework for benchmarking different phase retrieval algorithms
% using synthetic signals and either synthetic or real meaurement matrices.
%
% I/O
% Inputs
% xitem:  A string that describes the desired x-axis label of the plot that
%         that is produced by this benchmark.  Valid options are
%         ['m/n', 'snr','iterations', 'time'].  When using a 2D image with
%         Fourier measurements, then 'masks' should be used instead of
%         'm/n' to conntrol the number of Fourier masks.
% xvalues: A list of scalar values for which the performance of each
%          algorithm is measured.
% yitem:  A string to appear as the y-axis label.  Value should be
%         drawn from ['reconError', 'measurementError', 'correlation'];
% algorithms: a cell array of options structs,
%             where each struct is the same as the input parameter 'opts'
%             for solvePhaseRetrieval. See the example scripts for details.
% dataSet: The name of dataset used. Currently supported options are
%          ['1DGaussian', '2DImage', 'transmissionMatrix'].
%
% params: a struct of options containing the following fields:
%          verbose(boolean, default=false):
%                   If true, the result of each trial will be reported.
%          numTrials(integer, default=1):
%                   The number of trials each algorithm/dataset
%                   combination will run.
%          policy(string, default='median'):
%                   How to compute the final yvalue used for ploting from
%                   the values one gets by running numTrials trials. It
%                   currently supports
%                   ['median','mean','best','successRate'].
%          successConstant(real number,defualt=1e-5):
%                   If the yvalue of the current  trial is less than this,
%                   the trial will be counted as a success. This parameter
%                   will only be used when policy='successRate'.
%          maxTime(positive real number,default=120) :
%                   Max time allowed for a single algorithm.
%          recordSignals(boolean, default=false):
%                   Whether to record the recovered signal at each trial.
%
%
% Outputs
%          results : A 3D struct consisting of the errors(error
%                   metric is based on yitem chosen by the user) of all
%                   trials of algorithm/dataset combinations. Coordinates
%                   are (x-axis value, algorithm index, trial number).
% recoveredSignals: A 4D cell array consisting of the recovered signal at
%                   each trial for each algorithm. Coordinates are (x-axis
%                   value, algorithm index, current trial number, the index
%                   of the recovered signal).
%
% For more details, please look at the Phasepack User Guide.
%
%
% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

%% -----------------------------START----------------------------------

% Main benchmark interface
function [finalResults, results, recoveredSignals] = benchmarkSynthetic(xitem, xvalues, yitem, algorithms, dataSet, params)

% Check if the inputs are valid. Note: the fields of algorithms will be checked by
% solvePhaseRetrieval implicitly.
checkValidityOfInputs(xitem, xvalues, yitem, dataSet);

% If params is not provided, create it.
if ~exist('params', 'var')
    params = struct;
end

% Provide default value/validate the params provided by the user.
params = manageOptionsForBenchmark(dataSet, params);

% Check if the params are valid for the dataSet chosen by the user.
params = checkValidityOfParams(dataSet, params, xitem);

% Get the labels for x, y axis, numTrials, policy, successConstant,
% and recordSignals.
% For details of what these represent, see the header in this file or the User Guide.
numTrials = params.numTrials;
policy = params.policy;
successConstant = params.successConstant;
recordSignals = params.recordSignals;

% create struct to store results of each trial.
results = zeros(length(xvalues), length(algorithms), numTrials);
if recordSignals
    recoveredSignals = cell(length(xvalues), length(algorithms), numTrials);
end

fprintf('Benchmarking on dataset %s with %s as x axis and the %s\n %s of %d trials as y axis...\n\n',dataSet,xitem,policy,yitem,numTrials);

% Loop over the xvalues
for p=1:length(xvalues)
    fprintf('Running trails: %s=%g\n',xitem,xvalues(p))
    % Loop over the algorithms
    for k=1:length(algorithms)
        [opts,params] = setupTrialParameters(algorithms{k},xitem,xvalues(p),dataSet,params);
        fprintf('  %s:',algorithms{k}.algorithm);
        if numTrials==1
            fprintf('\n');
        end
        % Loop over the random trials
        for q=1:numTrials
            if numTrials>1 && params.verbose==0
                fprintf('*');
            end

            % Create a random test problem with the right dimensions and SNR
            [A, At, b0, xt, plotter, params] = createProblemData(dataSet, params);
            n = numel(xt);
            opts.xt = xt;

            % Call the solver specified in opts
            startTime = tic;
            [x, outs, opts] = solvePhaseRetrieval(A, At, b0, n, opts);
            elapsedTime = toc(startTime);

            % Update plot
            plotter(x);
            title(sprintf('%s (%s=%s)',opts.algorithm,xitem,num2str(xvalues(p))),'fontsize',16);
            drawnow();

            % Calculate the value the user is looking for.
            yvalue = evalY(yitem, x, xt, A, b0);
            % Store the value for this trial in a table.
            results(p, k, q) = yvalue;
            if recordSignals
                recoveredSignals{p,k,q} = x;
            end

            % Report results of a trial if verbose is true.
            if params.verbose
                reportResult(opts.initMethod, opts.algorithm, xitem, xvalues(p),...
                    yitem, yvalue, elapsedTime, q);
            end
        end
        fprintf('\n');
    end
end

% Get final results in order to plot a comparison graph for algorithms
% at different x values.
finalResults = getFinalResults(results, policy, successConstant, yitem);

% Plot the comparison of performance graph among different chosen algorithms and
% initial methods combinations.
plotComparison(xitem,yitem,xvalues,algorithms,dataSet,finalResults,policy);

end

%% Helper functions

% Set the parameters needed to solve a phase retrieval problem with the
% specified dataset, xtime, and yitem.
function [opts,params] = setupTrialParameters(opts,xitem,xval,dataSet,params)
switch lower(xitem)
    case 'iterations'
        opts.maxIters = xval;    % Update algorithm's max iterations.
        opts.tol = 1e-10;              % Set algorithm's tolerance to a very
        % small number in order to make it
        % run maxIters iterations.
        opts.maxTime = params.maxTime; % Set a time limit for each single trial.
        opts.isComplex = params.isComplex;
        opts.isNonNegativeOnly = params.isNonNegativeOnly;
    case 'm/n'
        if strcmp(lower(dataSet),'1dgaussian') || strcmp(lower(dataSet),'transmissionmatrix')
            params.m = round(params.n * xval); % Update m according m/n ratio.
        end
        opts.maxTime = params.maxTime; % Set a time limit for each single trial
    case 'snr'
        opts.maxTime = params.maxTime; % Set a time limit for each single trial.
        params.snr = xval;
    case 'time'
        opts.maxTime = xval;   % Set a time limit for each single trial.
        opts.tol = 1e-10;            % Set algorithm's tolerance to be very small
        opts.maxIters = 1e8;         % and algorithm's maxIters to be very large so it runs maxTime.
        opts.recordTimes = false;    % To save space, record nothing
        opts.recordResiduals = false;
    case 'masks'
        params.numMasks = xval;
        opts.maxTime = params.maxTime; % Set a time limit for each single trial
    case 'angle'
        opts.maxTime = params.maxTime; % Set a time limit for each single trial.
        opts.initMethod = 'angle';
        opts.initAngle = xval;
    otherwise
        % Raise an error if the given label for x axis is invalid.
        error(['invalid x label: ' xitem]);
end
end

% Run a specified algorithm on a specific dataset and get results
% Inputs:
% dataSet and params are as defined in benchmarkPhaseRetrieval
% opts: a struct consists of options specified by user for running the algorithm
%       on a specified dataset.
% Outputs:
% x:  n x 1 vector, estimation of xt given by the PR algorithm the subrountine invokes.
% xt: n x 1 vector, the real unknown signal.
% A:  m x n matrix/function handle(n x 1 -> m x 1).
% b0: m x 1 vector, the measurement.
% opts: a struct consists of options finally used for running the algorithm on a specified
%       dataset.
function [A, At, b0, xt, plotter, params] = createProblemData(dataSet, params)
    switch lower(dataSet)
        case '1dgaussian'
            [A, At, b0, xt, plotter] = experimentGaussian1D(params.n, params.m,...
                params.isComplex, params.isNonNegativeOnly);
        case '2dimage'
            [A, At, b0, xt, plotter] = experimentImage2D(params.numMasks, params.imagePath);
        case 'transmissionmatrix'
            [A, b0, xt, plotter] = experimentTransMatrixWithSynthSignal(params.n, params.m, params.A_cached);
            params.A_cached = A;
            At = [];
        otherwise
            error('unknown dataset: %s\n',dataSet);
    end

    % Add noise to achieve specified SNR
    if params.snr ~= inf
        noise = randn(params.m,1); % create noise
        noise = noise/norm(noise)*norm(b0)/params.snr;
        b0 = max(b0+noise,0);
    end

end

% Calculate how good the solution is. Use the metric specified by yitem.
function yvalue=evalY(yitem, x, xt, A, b0)
    %  Compute optimal rotation
    switch lower(yitem)
        case 'reconerror'
            % solve for least-squares solution:  alpha*x = xt
            alpha = (x(:)'*xt(:))/(x(:)'*x(:));
            x = alpha*x;
            yvalue = norm(xt(:)-x(:))/norm(xt(:));
        case 'measurementerror'
            % Transform A into function handle if A is a matrix
            if isnumeric(A)
                At = @(x) A'*x;
                A = @(x) A*x;
            end
            yvalue = norm(abs(A(x))-b0(:))/norm(b0(:));
        case 'correlation'
            yvalue = abs(x'*xt/norm(x)/norm(xt));
        otherwise
            error(['invalid y label: ' yitem]);
    end
end

% Get final results by averaging across all trials.
% The possible policies are the following:
% mean: take the mean of all trials.
% best: take the best of all trials. If yitem=='reconerror' or 'measurementerror',
%       min value will be taken; If yitem=='correlation', max value will be taken.
% median: take the median of all trials.
% successrate: A success rate will be calculated. If yitem=='reconerror' or
%              'measurementerror', it is the percentage of values that are
%              smaller than the successConstant. If yitem=='correlation', it is
%              the percentage of values that are larger than the successConstant.
% The input struct results has size length(xvalues) x length(algorithms) x numTrials
% The output struct finalResults has size length(xvalues) x length(algorithms)
function finalResults = getFinalResults(results, policy, successConstant, yitem)
    switch lower(policy)
        case 'mean'
            finalResults = mean(results,3);
        case 'best'
            switch lower(yitem)
                case {'reconerror', 'measurementerror'}
                    finalResults = min(results,[],3);
                case 'correlation'
                    finalResults = max(results,[],3);
                otherwise
                    error('invalid yitem: %s',yitem);
            end
        case 'median'
            finalResults = median(results,3);
        case 'successrate'
            switch lower(yitem)
                case {'reconerror', 'measurementerror'}
                    finalResults = mean(results<successConstant,3);
                case 'correlation'
                    finalResults = mean(results>successConstant,3);
                otherwise
                    error('invalid yitem: %s',yitem);
            end
        otherwise
            error('Invalid policy: %s', policy);
    end
end

% Plot a performance curve for each algorithm
function plotComparison(xitem,yitem,xvalues,algorithms,dataSet,finalResults,policy)
    algNames = {};
    % get the labels to appear in the legend
    for k=1:length(algorithms)
        if isfield(algorithms{k},'label') && numel(algorithms{k}.label)>0             % use user-specified label if available
            algNames{k} =  algorithms{k}.label;
        else
            algNames{k} = algorithms{k}.algorithm;  % otherwise use the algorithm name
        end
    end

    autoplot(xvalues,finalResults,algNames);

    title(char(strcat({xitem},{' VS '},{yitem},{' on '},{dataSet})));                         % title for plot
    xlabel(xitem);                                       % x-axis label
    ylabel(char(strcat({policy},{'('},{yitem},{')'})));  % y-axis label
    ax = gca;                                           % create an axes object
    ax.FontSize = 12;                                   % adjust font size on the axes
    legend('show', 'Location', 'northeastoutside');     % show the legend
end

% Report the result at the end of each trial if verbose is true.
function reportResult(initMethod, algorithm, xitem, xvalue, yitem, yvalue, time, currentTrialNum)
    if currentTrialNum==1
        fprintf('\n');
    end
    fprintf('Trial: %d, initial method: %s, algorithm: %s, %s: %d, %s: %f, time: %f\n',...
        currentTrialNum, initMethod, algorithm, xitem, xvalue, yitem, yvalue, time);
end

% Check if the inputs are valid.
function checkValidityOfInputs(xitem, xvalues, yitem, dataSet)
    assert(~isempty(xvalues), 'The list xvalues must not be empty.');

    % Check if yitem is a valid choice
    yitemList = {'reconerror', 'measurementerror', 'correlation'};
    checkIfInList('yitem', yitem, yitemList);

    % Check if dataSet chosen supports xitem
    xitem = lower(xitem);
    switch lower(dataSet)
        case '1dgaussian'
            supportedList = {'iterations','m/n','snr','time','angle'};
        case 'transmissionmatrix'
            supportedList = {'iterations','m/n','snr','time','angle'};
        case '2dimage'
            supportedList = {'masks','iterations','angle'};
        otherwise
            error('unknown dataset: %s\n',dataSet);
    end
    checkIfInList('xitem',xitem,supportedList,strcat('For ',lower(dataSet),', '));
end

% Check if the params are valid
function params = checkValidityOfParams(dataSet, params, xitem)
    switch lower(dataSet)
        case '1dgaussian'
            assert(isfield(params,'n'),'User must specify signal dimension in params.n');
            checkIfNumber('params.n',params.n);
            if ~strcmp(xitem,'m/n')
                assert(isfield(params,'m'), ['User must specify the number of measurements in params.m when using xlabel ',xitem]);
            end
        case '2dimage'
            assert(~strcmp(xitem,'m/n'),'xlabel m/n is not supported when using 2D images.  Instead use "masks"');
            if ~(xitem=='masks')
                assert(exist('params.numMasks'), ['User must specify the number of Fourier masks in params.numMasks when using xlabel ',xitem]);
            end
        case 'transmissionmatrix'
            params.A_cached = [];
            assert(isfield(params,'n'),'User must specify signal dimension in params.n.  Options are {256,1600,4096} for the tranmissionMatrix dataset.');
            assert(params.n==256 || params.n==1600 || params.n==4096,...
                ['Invalid choice (',num2str(params.n),') ofr params.n for the transmissionMatrix dataset. Valid choices are {256,1600,4096}.'])
        otherwise
            error('unknown dataset: %s\n',dataSet);
    end

    if ~isfield(params,'snr')
        params.snr = Inf;
    end

end

author: tomgoldstein

repo: phasepack-matlab


text:
filename='RVE_Square_Triangle';
ptype = 'MICRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'circleInclusion';
cost={'volume','perimeter'};
weights=[1 0.1];
constraint = {'enforceCh_CCstar_eq'};%inf & equality
optimizer = 'MMA'; incrementFactor = 1;
filterType = 'P1';

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

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

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

%Micro
epsilon_isotropy_initial=1e-1;
epsilon_isotropy_final = 1e-3;
selectiveC_Cstar = [1,1,1;
    1,1,1;
    1,1,1]; % 1 to select the component

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function IIDAnalysis(Dates,Data)
% this function performs simple invariance (i.i.d.) tests on a time series
% 1. it checks that the variables are identically distributed by looking at the
%    histogram of two subsamples
% 2. it checks that the variables are independent by looking at the 1-lag scatter plot
% under i.i.d. the location-dispersion ellipsoid should be a circle
% see "Risk and Asset Allocation"-Springer (2005), by A. Meucci

% plot time series

figure
h=plot(Dates,Data,'.');
xlim([Dates(1) Dates(end)])
datetick('x','mmmyy','keeplimits','keepticks');
grid on

% test "identically distributed hypothesis": split observations into two sub-samples and plot histogram
Sample_1=Data(1:round(length(Data)/2));
Sample_2=Data(round(length(Data)/2)+1:end);
num_bins_1=round(5*log(length(Sample_1)));
num_bins_2=round(5*log(length(Sample_2)));
X_lim=[min(Data)-.1*(max(Data)-min(Data)) max(Data)+.1*(max(Data)-min(Data))];
[n1,xout1]=hist(Sample_1,num_bins_1);
[n2,xout2]=hist(Sample_2,num_bins_2);

figure
subplot('Position',[0.03 .59 .44 .35])
h1=bar(xout1,n1,1);
set(h1,'FaceColor',[.7 .7 .7],'EdgeColor','k')
set(gca,'ytick',[],'xlim',X_lim,'ylim',[0 max(max(n1),max(n2))])
grid off

subplot('Position',[.53 .59 .44 .35])
h2=bar(xout2,n2,1);
set(h2,'FaceColor',[.7 .7 .7],'EdgeColor','k');
set(gca,'ytick',[],'xlim',X_lim,'ylim',[0 max(max(n1),max(n2))]);
grid off

% test "independently distributed hypothesis": scatter plot of observations at lagged times
subplot('Position',[.28 .01 .43 .43])
X=Data(1:end-1);
Y=Data(2:end);
h3=plot(X,Y,'.');
grid off
axis equal
set(gca,'xlim',X_lim,'ylim',X_lim);

m=mean([X Y])';
S=cov([X Y]);
TwoDimEllipsoid(m,S,2,0,0);


text:
%% Copyright 2012 The MathWorks, Inc.
% Create a detector object
faceDetector = vision.CascadeObjectDetector;

% Read input image
I = imread('visionteam.jpg');

% Detect faces
bbox = step(faceDetector, I);

% Create a shape inserter object to draw bounding boxes around detections
shapeInserter = vision.ShapeInserter('BorderColor','Custom','CustomBorderColor',[255 255 0]);

% Draw boxes around detected faces and display results
I_faces = step(shapeInserter, I, int32(bbox));
figure, imshow(I_faces), title('Detected faces');
%%
% Create a detector object
bodyDetector = vision.CascadeObjectDetector('UpperBody');
bodyDetector.MinSize = [60 60];
bodyDetector.ScaleFactor = 1.05;

% Read input image and detect upper body
bbox_body = step(bodyDetector, I);

% Draw bounding boxes
shapeInserter = vision.ShapeInserter('BorderColor','Custom','CustomBorderColor',[255 255 0]);
I_body = step(shapeInserter, I, int32(bbox_body));
figure, imshow(I_body);
%%
bbox_face = zeros(size(bbox_body));
for i=1:length(bbox_body)
    Icrop = imcrop(I,bbox_body(i,:));
    bbox = step(faceDetector,Icrop);
    bbox_face(i,:) = bbox + [bbox_body(i,1:2)-1 0 0];
end

I_faces2 = step(shapeInserter, I, int32(bbox_face));
figure, imshow(I_faces2);
%%
Icrop = imcrop(I,bbox_body(1,:));
figure;imshow(Icrop);
bbox = step(faceDetector,Icrop);
hold on;rectangle('Position',bbox,'EdgeColor','y');
%%
x = 5;
Irotate = imrotate(Icrop,x);
imshow(Irotate);
bbox = step(faceDetector, Irotate);
if bbox > 0
    hold on;rectangle('Position',bbox,'EdgeColor','y'); hold off;
end


text:
function sftpack_test03 ( )

%*****************************************************************************80
%
%% TEST03 tests R8VEC_SHT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2010
%
%  Author:
%
%    John Burkardt
%
  n = 17;
  alo = 0.0;
  ahi = 5.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST03\n' );
  fprintf ( 1, '  For real slow Hartley transforms,\n' );
  fprintf ( 1, '  R8VEC_SHT does a forward or backward transform.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of data items is N = %d\n', n );
%
%  Set the data values.
%
  seed = 123456789;

  [ c, seed ] = r8vec_uniform ( n, alo, ahi, seed );

  r8vec_print_part ( n, c, 10, '  The original data:' );
%
%  Compute the coefficients.
%
  d = r8vec_sht ( n, c );

  r8vec_print_part ( n, d, 10, '  The Hartley coefficients:' );
%
%  Now compute inverse transform of coefficients.  Should get back the
%  original data.

  e = r8vec_sht ( n, d );

  r8vec_print_part ( n, e, 10, '  The retrieved data:' );

  return
end


text:
function nc2010_vbpca_mnistexperiment(flatW, rotate, ncomps)

M = load('mnist_all');

n = 100;

Y = double(M.train5');
inds = randperm(cols(Y));
Y = Y(:,inds(1:n));
size(Y)

m = rows(Y);

startupdatehyper = 1
maxiters = 1e4;

randn('state', 1);
rand('state', 1);

% Divide into training (90%) and test sets (10%)
testset = (rand(size(Y)) < 0.1);
Ytrain = Y;
Ytest = Y;
Ytrain(testset) = nan;
Ytest(~testset) = nan;

%size(Ytrain)

if flatW
  error('Using flat W might be problematic because m>n ???');
  disp('Using non-hierarchical flat prior for W');
  stringW = 'flatW';
else
  disp('Using hierarchical prior for W');
  stringW = 'hierW';
end

filename = sprintf(['/home/jluttine/matlab/neurocomputing2010/' ...
                    'nc2010_vbpca_mnistexperiment_%s_ncomps=%d_rotate=%d'], ...
                   stringW, ncomps, rotate);

results = vbpcamv(Ytrain, ncomps, 'testset', Ytest, 'maxiters', 1e7, ...
                  'rotate', rotate, 'startupdatehyper', startupdatehyper, ...
                  'startrotate', 1, 'autosavetime', 1, ...
                  'autosavefile', filename, 'fixw', flatW);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [ node_num2, triangle_num2, edge_data ] = refine_size ( node_num1, ...
  triangle_num1, triangle_node1 )

%*****************************************************************************80
%
%% TRIANGULATION_ORDER3_REFINE_SIZE sizes a refined order 3 triangulation.
%
%  Discussion:
%
%    Given a triangle defined by nodes 1, 2, 3, we need to generate
%    nodes 12, 23, and 13, and create 4 new subtriangles, T1, T2, T3
%    and T4.
%
%    The task is more complicated by the fact that we are working with
%    a mesh of triangles, so that we want to create a node only once,
%    even though it may be shared by other triangles.
%
%          3
%         / \
%        /T3 \
%      13----23
%      / \T4 / \
%     /T1 \ /T2 \
%    1----12-----2
%
%    This routine simply determines the sizes of the resulting node
%    and triangle arrays.
%
%    The primary amount of work occurs in sorting a list of 3 * TRIANGLE_NUM
%    data items, one item for every edge of every triangle.  Each
%    data item records, for a given edge, the global indices
%    of the two endpoints, the local indices of the two endpoints,
%    and the index of the triangle.
%
%    Through careful sorting, it is possible to arrange this data in
%    a way that allows the proper generation of the interpolated nodes.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM1, the number of nodes in the original mesh.
%
%    Input, integer TRIANGLE_NUM1, the number of triangles in the
%    original mesh.
%
%    Input, integer TRIANGLE_NODE1(TRIANGLE_NUM1,3), the indices of the nodes
%    that form the triangles in the input mesh.
%
%    Output, integer NODE_NUM2, the number of nodes in the refined mesh.
%
%    Output, integer TRIANGLE_NUM2, the number of triangles in the
%    refined mesh.
%
%    Output, integer EDGE_DATA(3*TRIANGLE_NUM1,5), edge data that will
%    be needed by TRIANGULATION_ORDER3_REFINE_COMPUTE.
%
  edge_data = zeros(3*triangle_num1,5);
%
%  Step 1.
%  From the list of nodes for triangle T, of the form: (I,J,K)
%  construct the edge relations:
%
%    (I,J,1,2,T)
%    (I,K,1,3,T)
%    (J,K,2,3,T)
%
%  In order to make matching easier, we reorder each pair of nodes
%  into ascending order.
%
  row = 0;

  for triangle = 1 : triangle_num1

    i = triangle_node1(triangle,1);
    j = triangle_node1(triangle,2);
    k = triangle_node1(triangle,3);

    a = min ( i, j );
    b = max ( i, j );
    row = row + 1;
    edge_data(row,1:5) = [ a, b, 1, 2, triangle ];

    a = min ( j, k );
    b = max ( j, k );
    row = row + 1;
    edge_data(row,1:5) = [ a, b, 2, 3, triangle ];

    a = min ( k, i );
    b = max ( k, i );
    row = row + 1;
    edge_data(row,1:5) = [ a, b, 1, 3, triangle ];

  end
%
%  Step 2. Perform an ascending dictionary sort on the neighbor relations.
%  We only intend to sort on items in columns 1 and 2; the routine we call here
%  sorts on the full column but that won't hurt us.
%
%  What we need is to find all cases where triangles share an edge.
%  By sorting the columns of the EDGE_DATA array, we will put shared edges
%  next to each other.
%
  edge_data = sortrows ( edge_data );
% edge_data = i4col_sort_a ( 5, 3*triangle_num1, edge_data );
%
%  Step 3. All the triangles which share an edge show up as consecutive
%  columns with identical first two entries.  Figure out how many new
%  nodes there are, and allocate space for their coordinates.
%
  node_num2 = node_num1;

  n1_old = -1;
  n2_old = -1;

  for edge = 1 : 3 * triangle_num1
    n1 = edge_data(edge,1);
    n2 = edge_data(edge,2);
    if ( n1 ~= n1_old || n2 ~= n2_old )
      node_num2 = node_num2 + 1;
      n1_old = n1;
      n2_old = n2;
    end
  end

  triangle_num2 = 4 * triangle_num1;

  return
end


text:
%KMEANS PRTools k-means clustering, deprecated
%
%   [LABELS,B] = KMEANS(A,K,MAXIT,INIT)
%
% INPUT
%  A       Dataset
%  K       Number of clusters to be found (optional; default: 2)
%  MAXIT   maximum number of iterations (optional; default: 50)
%  INIT    Labels for initialisation, or
%          'rand'     : take at random K objects as initial means, or
%          'kcentres' : use KCENTRES for initialisation (default)
%
% OUTPUT
%  LABELS  Cluster assignments, 1..K
%  B       Dataset with original data and labels LABELS:
%          B = PRDATASET(A,LABELS)
%
% DESCRIPTION
% K-means clustering of data vectors in A. This routine calls PRKMEANS
%
% SEE ALSO (PRTools Guide)
% DATASETS, PRKMEANS, HCLUST, KCENTRES, MODESEEK, EMCLUST, PRPROGRESS

author: marianux

repo: ecg-kit

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


text:
function [spls, corresp, spls_adj, joints, root_id] = merge_nearby_joints(spls, corresp, spls_adj, joints, root_id)
while true
    for i=joints(:,1)'
        edges = find( spls_adj(i,:)==1 );
        edges = edges(edges~=i);
        tmp = ismember(edges, joints(:,1)');
        if sum(tmp)~=1, continue, end;

        j = edges(tmp);
        edges = find( spls_adj(j,:)==1 );
        edges = edges(edges~=j);
        tmp = ismember(edges, joints(:,1)');
        if sum(tmp)~=1, continue, end;

      % update the location
        spls( i,: ) = mean( spls( [i,j],: ) );
        spls( j,: ) = NaN;
        % update the correspondents
       corresp(corresp==j ) = i;

        % update the A matrix
        for k=1:size(spls_adj,1)
            if spls_adj( j,k ) == 1,
                spls_adj( i,k )=1;
                spls_adj( k,i)=1;
            end
        end
        % remove the row
        spls_adj( j,: ) = 0;
        spls_adj( :,j ) = 0;

        segments(j) = 0;
        j = find( joints(:,1)==j );
        if root_id == j
            root_id = find( joints(:,1)==i );
        end
        joints(j,:) = [];
%         break;
    end
    if i == joints(end,1), break, end;
end

%%
figure('Name','Merge nearby joints','NumberTitle','off');set(gcf,'color','white');view3d rot;
movegui('north');
plot_skeleton(spls, spls_adj);
axis off; axis equal; camorbit(0,0,'camera'); axis vis3d; view(-90,0);view3d rot;

author: taiya

repo: cloudcontr

save_path: github-repos/MATLAB/taiya-cloudcontr


text:
function varargout = drawDigraph(varargin)
%DRAWDIGRAPH Draw a directed graph, given as a set of vertices and edges.
%
%   drawDigraph(NODES1, NODES2, EDGES)
%   NODES1 are originating vertices
%   NODES2 are destination vertices
%   EDGES is an array, with first column containing index of origin vertex
%   (index in NODES1), and second column containing index of destination
%   vertex (index in NODES2).
%   Edges are drawn with arrows.
%
%   H = drawDigraph(...)
%   return handle to the set of edges.
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2004-08-17
% Copyright 2004-2022 INRA - TPV URPOI - BIA IMASTE

%% Initialisations

% check number of arguments
if nargin < 3
    help drawDigraph;
    return;
end

% initialisations
sn1 = 'bo';     % nodes are red circles
sn2 = 'ro';     % nodes are red circles

%% process input arguments

% First extract the graph structure
n1 = varargin{1};
n2 = varargin{2};
e  = varargin{3};
varargin = varargin(4:length(varargin));

% extract drawing style
if ~isempty(varargin)
    sn1 = varargin{1};
end

if length(varargin)>1
    sn2 = varargin{2};
end

%% main drawing processing

hold on;

nodes = [n1 ; n2];
e(:,2) = e(:,2)+length(n1);

% Draw a 2 dimensional directed graph

% Draw 2D Edges
%if ~strcmp(se, 'none') & size(e, 1)>0
%    he = plot([n1(e(:,1),1) n2(e(:,2),1)]', [n1(e(:,1),2) n2(e(:,2),2)]', se);
%end
he = drawDirectedEdges(nodes, e);

hn = [];
% Draw 2D nodes
if ~strcmp(sn1, 'none')
    hn = plot(n1(:,1), n1(:,2), sn1);
end

 % Draw 2D nodes
if ~strcmp(sn2, 'none')
    hn = plot(n2(:,1), n2(:,2), sn2);
end

%% format output arguments

if nargout==1
    varargout{1} = he;
end

if nargout==2
    varargout{1} = hn;
    varargout{2} = he;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom