text:
classdef (Abstract, CaseInsensitiveProperties = true) functionalPlot < ...
        iosr.statistics.statsPlot
%FUNCTIONALPLOT Abstract superclass for functional plots
%
%   As an abstract class, this class cannot be instantiated. It provides no
%   public methods. The class is a super class for:
%       - iosr.statistics.functionalBoxPlot
%       - iosr.statistics.functionalSpreadPlot
%
%   The class does, however, provide various properties to subclasses.
%
%   IOSR.STATISTICS.FUNCTIONALSPREADPLOT properties:
%       mainLineColor           - Specifies the color of the central line.
%                                 See 'addPrctilesLineColor' for details of
%                                 how to specify colors. The default is
%                                 'auto'.
%       mainLineStyle           - Specifies the style of the central line.
%                                 The property may be specified as a char
%                                 array (e.g. '-'), or as a length-P cell
%                                 vector of strings (if passing line styles
%                                 for each functional plot). The default is
%                                 '-'.
%       mainLineWidth           - Specifies the width of the central line.
%                                 The property may be specified as a
%                                 length-P numeric vector, or as a cell
%                                 vector (if passing line styles for each
%                                 functional plot). The default is 1.
%       outerLineColor          - Specifies the color of the outer lines.
%                                 See 'mainLineStyle' for details of how to
%                                 specify colors. The default is 'auto'.
%       outerLineStyle          - Specifies the style of the outer lines.
%                                 See 'mainLineStyle' for details of how to
%                                 specify styles. The default is '-'.
%       outerLineWidth          - Specifies the width of the outer line.
%                                 See 'mainLineWidth' for details of how to
%                                 specify widths. The default is 1.
%       outlierLineWidth        - Specifies the width of the outlier marker
%                                 line. See 'mainLineWidth' for details of
%                                 how to specify widths. The default is 1.
%       showOutliers            - Specifies whether to show outliers. The
%                                 property must be a logical value. The
%                                 default is false.
%       spreadAlpha             - Set the alpha of the shaded region(s).
%                                 The default is 0.5.
%       spreadBorderLineColor   - Set the color of the shaded regions'
%                                 / region's border. See
%                                 'addPrctilesLineColor' for details of how
%                                 to specify colors. The default is 'none'.
%       spreadBorderLineWidth   - Set the width of the shaded regions'
%                                 / region's border. The default is 1.
%       spreadColor             - Set the color of the shaded region(s).
%                                 See 'addPrctilesLineColor' for details of
%                                 how to specify colors. The default is
%                                 'auto'.
%       whiskers                - Specify whether the box is connected to
%                                 the outer line via a whisker. The line is
%                                 drawn with the same style as the outer
%                                 line. The default is false.
%
%   These properties can be referenced using dot notation - e.g. H.BOXCOLOR
%   where H is an instance of the FUNCTIONALSPREADPLOT object - or using
%   the SET and GET methods - e.g. GET(H,'BOXCOLOR'). Both methods are
%   case-insensitive.
%
%   Note that some handles will be empty unless the associated option is
%   turned on.
%
%   See also IOSR.STATISTICS.BOXPLOT, IOSR.STATISTICS.QUANTILE,
%   IOSR.STATISTICS.FUNCTIONALBOXPLOT,
%   IOSR.STATISTICS.FUNCTIONALSPREADPLOT, IOSR.STATISTICS.STATSPLOT.

    properties (AbortSet)
        mainLineColor = 'auto'          % The color of the central line.
        mainLineStyle = '-'             % The style of the central line.
        mainLineWidth = 2               % The width of the central line.
        outerLineColor = 'auto'         % The color of the outer lines.
        outerLineStyle = '-'            % The style of the outer lines.
        outerLineWidth = 0.5            % The width of the outer lines.
        outlierLineWidth = 1            % Width of the outlier marker edges.
        showOutliers = true             % Turn outliers on and off.
        spreadAlpha = 0.5               % The alpha of the shaded regions.
        spreadColor = 'auto'            % The color of the shaded regions.
        spreadBorderLineWidth = 1       % The width of the shaded region borders.
        spreadBorderLineColor = 'none'  % The color of the shaded region borders.
        whiskers = false                % Specify whether a line connects the box to the outer lines.
    end

    methods

        function set.showOutliers(obj, val)
            assert(numel(val)==1 && islogical(val), 'iosr:functionalPlot:invalidShowOutliers', ...
                '''SHOWOUTLIERS'' must be true or false.')
            obj.showOutliers = val;
            obj.draw();
        end

        % main lines

        function val = get.mainLineColor(obj)
            val = obj.parseColor(obj.mainLineColor);
        end

        function set.mainLineColor(obj,val)
            obj.mainLineColor = val;
            obj.draw();
        end

        function val = get.mainLineStyle(obj)
            val = obj.parseProps(obj.mainLineStyle, false);
        end

        function set.mainLineStyle(obj,val)
            obj.mainLineStyle = val;
            obj.draw();
        end

        function val = get.mainLineWidth(obj)
            val = obj.parseProps(obj.mainLineWidth, false);
        end

        function set.mainLineWidth(obj,val)
            obj.mainLineWidth = val;
            obj.draw();
        end

        % outer lines

        function val = get.outerLineColor(obj)
            val = obj.parseColor(obj.outerLineColor);
        end

        function set.outerLineColor(obj,val)
            obj.outerLineColor = val;
            obj.draw();
        end

        function val = get.outerLineStyle(obj)
            val = obj.parseProps(obj.outerLineStyle, false);
        end

        function set.outerLineStyle(obj,val)
            obj.outerLineStyle = val;
            obj.draw();
        end

        function val = get.outerLineWidth(obj)
            val = obj.parseProps(obj.outerLineWidth, false);
        end

        function set.outerLineWidth(obj,val)
            obj.outerLineWidth = val;
            obj.draw();
        end

        % outlier settings

        function val = get.outlierLineWidth(obj)
            val = obj.parseProps(obj.outlierLineWidth, false);
        end

        function set.outlierLineWidth(obj,val)
            obj.outlierLineWidth = val;
            obj.draw();
        end

        % set spread alpha

        function set.spreadAlpha(obj,val)
            assert(isnumeric(val) && isscalar(val), 'iosr:functionalPlot:invalidSpreadAlpha', ...
                '''SPREADALPHA'' must be a numeric scalar')
            obj.spreadAlpha = val;
            obj.draw();
        end

        % get/set spread line colors

        function val = get.spreadBorderLineColor(obj)
            val = obj.parseColor(obj.spreadBorderLineColor);
        end

        function set.spreadBorderLineColor(obj,val)
            obj.spreadBorderLineColor = val;
            obj.draw();
        end

        % set spread border line width

        function set.spreadBorderLineWidth(obj,val)
            assert(isnumeric(val) && isscalar(val), 'iosr:functionalPlot:invalidspreadBorderLineWidth', ...
                '''SPREADBORDERLINEWIDTH'' must be a numeric scalar')
            obj.spreadBorderLineWidth = val;
            obj.draw();
        end

        % set/get spread color

        function val = get.spreadColor(obj)
            val = obj.parseColor(obj.spreadColor);
        end

        function set.spreadColor(obj,val)
            obj.spreadColor = val;
            obj.draw();
        end

        % whisker

        function set.whiskers(obj,val)
            assert(islogical(val) && isscalar(val), 'iosr:functionalPlot:invalidWhiskers', ...
                '''whisker'' must be a scalar and logical');
            obj.whiskers = val;
            obj.draw();
        end

    end

    methods (Access = protected)

        % Draw the plot
        function draw(obj)

            if isfield(obj.handles,'axes')

                axes(obj.handles.axes);
                hold on;

                if length(obj.outDims) > 2
                    nlines = obj.outDims(3);
                else
                    nlines = 1;
                end
                for n = nlines:-1:1
                    % draw inner patch
                    obj.handles.spreads(n) = patch(...
                        [obj.x obj.x(end:-1:1)], ...
                        [obj.statistics.inner_u(:,:,n) obj.statistics.inner_l(:,end:-1:1,n)], ...
                        obj.spreadColor{n} , ...
                        'FaceAlpha', obj.spreadAlpha, ...
                        'EdgeColor', obj.spreadBorderLineColor{n}, ...
                        'LineWidth', obj.spreadBorderLineWidth ...
                    );
                end

                for n = nlines:-1:1

                    % draw outer limits
                    obj.handles.outer_u(n) = line(obj.x, obj.statistics.outer_u(:,:,n), ...
                        'linestyle', obj.outerLineStyle{n}, ...
                        'linewidth', obj.outerLineWidth{n}, ...
                        'color', obj.outerLineColor{n} ...
                    );
                    obj.handles.outer_l(n) = line(obj.x, obj.statistics.outer_l(:,:,n), ...
                        'linestyle', obj.outerLineStyle{n}, ...
                        'linewidth', obj.outerLineWidth{n}, ...
                        'color', obj.outerLineColor{n} ...
                    );

                    % draw centre line
                    obj.handles.main(n) = line(obj.x, obj.statistics.main(:,:,n), ...
                        'linestyle', obj.mainLineStyle{n}, ...
                        'linewidth', obj.mainLineWidth{n}, ...
                        'color', obj.mainLineColor{n} ...
                    );

                end

                % draw whiskers
                if obj.whiskers
                    x_range = max(obj.x) - min(obj.x);
                    whisker_xlim = [(min(obj.x) + (1/3)*x_range) (max(obj.x) - (1/3)*x_range)];
                    if nlines==1
                        whisker_x = mean(whisker_xlim);
                    else
                        whisker_x = linspace(whisker_xlim(1),whisker_xlim(2),nlines);
                    end
                    for n = nlines:-1:1
                        % interpolate to find y
                        [~,x_ix] = min(abs(obj.x-whisker_x(n)));
                        x_interp_ix = [max(1,x_ix-1) x_ix min(length(obj.x),x_ix+1)];
                        x_interp_ix = unique(x_interp_ix);
                        y_i_l = interp1(obj.x(x_interp_ix), obj.statistics.inner_l(:,x_interp_ix,n), whisker_x(n));
                        y_i_u = interp1(obj.x(x_interp_ix), obj.statistics.inner_u(:,x_interp_ix,n), whisker_x(n));
                        y_o_l = interp1(obj.x(x_interp_ix), obj.statistics.outer_l(:,x_interp_ix,n), whisker_x(n));
                        y_o_u = interp1(obj.x(x_interp_ix), obj.statistics.outer_u(:,x_interp_ix,n), whisker_x(n));

                        % plot
                        obj.handles.whiskers_l(n) = line([whisker_x(n) whisker_x(n)],[y_i_l y_o_l], ...
                            'linestyle', obj.outerLineStyle{n}, ...
                            'linewidth', obj.outerLineWidth{n}, ...
                            'color', obj.outerLineColor{n});
                        obj.handles.whiskers_u(n) = line([whisker_x(n) whisker_x(n)],[y_i_u y_o_u], ...
                            'linestyle', obj.outerLineStyle{n}, ...
                            'linewidth', obj.outerLineWidth{n}, ...
                            'color', obj.outerLineColor{n});
                    end
                else
                    try
                        delete(obj.handles.whiskers_l(n))
                        delete(obj.handles.whiskers_u(n))
                    catch
                    end
                end
            end

        end

    end

    methods (Access = protected)

        % Make properties in to a cell array
        function val = parseProps(obj, prop, addP)

            nLines = obj.outDims(3);
            if ~addP
                nAdd = 1;
            else
                nAdd = length(obj.addPrctiles);
            end
            val = cell(nLines, nAdd);
            if isnumeric(prop)
                cellprop = num2cell(prop);
            elseif ischar(prop)
                cellprop = cellstr(prop);
            else
                cellprop = prop;
            end
            try
            for n = 1:nLines
                if ~addP
                    val(n) = cellprop(mod(n-1, numel(cellprop)) + 1);
                else
                    for p = 1:nAdd
                        ixn = mod(n-1, size(cellprop, 1)) + 1;
                        ixp = mod(p-1, size(cellprop, 2)) + 1;
                        val(n,p) = cellprop(ixn, ixp);
                    end
                end
            end
            catch
                keyboard
            end

        end

        % Put colors in to a cell array
        function val = parseColor(obj, inColor)

            nLines = obj.outDims(3);
            val = cell(nLines, 1);
            if isnumeric(inColor)
                if size(inColor, 2) ~= 3
                    error('iosr:functionalPlot:colorInvalid','Color must be an N-by-3 array, where N is any positive integer.')
                end
                for n = 1:nLines
                    val(n,:) = {inColor(mod(n, size(inColor, 1)) + 1, :)};
                end
            elseif ischar(inColor)
                if strcmp(inColor, 'auto')
                    colors = lines(nLines);
                    for n = 1:nLines
                        val(n,:) = {colors(n, :)};
                    end
                else
                    colors = cellstr(inColor);
                    for n = 1:nLines
                        val(n,:) = colors(mod(n, size(inColor, 1)) + 1);
                    end
                end
            elseif isa(inColor,'function_handle')
                colors = inColor(nLines);
                for n = 1:nLines
                    val(n,:) = {colors(mod(n, numel(inColor)) + 1)};
                end
            elseif iscellstr(inColor)
                for n = 1:nLines
                    val(n,:) = inColor(mod(n, numel(inColor)) + 1);
                end
            else
                error('iosr:functionalPlot:invalidColorFormat','Invalid color format specified.')
            end

        end

    end

end

author: IoSR-Surrey

repo: MatlabToolbox


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
function kayttht
%KAYTTHT Unit test for function KAYTTH.

%	O. Lemoine - August 1996.

author: HeLiangHIT


text:
%Sensor Sensor superclass
%
% A superclass to represent robot navigation sensors.
%
% Methods::
%   plot        plot a line from robot to map feature
%   display     print the parameters in human readable form
%   char        convert to string
%
% Properties::
% robot   The Vehicle object on which the sensor is mounted
% map     The Map object representing the landmarks around the robot
%
% Reference::
%
%   Robotics, Vision & Control,
%   Peter Corke,
%   Springer 2011
%
% See also RangeBearing, EKF, Vehicle, Map.

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

classdef Sensor < handle
    % TODO, pose option, wrt vehicle

    properties
        robot
        map

        verbose

        ls
        animate     % animate sensor measurements
        interval    % measurement return subsample factor
        fail
        delay

    end

    methods

        function s = Sensor(robot, map, varargin)
        %Sensor.Sensor Sensor object constructor
        %
        % S = Sensor(VEHICLE, MAP, OPTIONS) is a sensor mounted on a vehicle described by the Vehicle class object
        % VEHICLE and observing landmarks in a map described by the Map class object MAP.
        %
        % Options::
        % 'animate'    animate the action of the laser scanner
        % 'ls',LS      laser scan lines drawn with style ls (default 'r-')
        % 'skip', I    return a valid reading on every I'th call
        % 'fail',T     sensor simulates failure between timesteps T=[TMIN,TMAX]
        %

            opt.skip = 1;
            opt.animate = false;
            opt.fail =  [];
            opt.ls = 'r-';
            opt.delay = 0.1;

            [opt,args] = tb_optparse(opt, varargin);

            s.interval = opt.skip;
            s.animate = opt.animate;

            s.robot = robot;
            s.map = map;
            s.verbose = false;
            s.fail = opt.fail;
            s.ls = opt.ls;
        end

        function plot(s, jf)
        %Sensor.plot Plot sensor reading
        %
        % S.plot(J) draws a line from the robot to the J'th map feature.
        %
        % Notes::
        % - The line is drawn using the linestyle given by the property ls
        % - There is a delay given by the property delay

            if isempty(s.ls)
                return;
            end

            h = findobj(gca, 'tag', 'sensor');
            if isempty(h)
                % no sensor line, create one
                h = plot(0, 0, s.ls, 'tag', 'sensor');
            end

            % there is a sensor line animate it

            xi = s.map.map(:,jf);
            set(h, 'XData', [s.robot.x(1), xi(1)], 'YData', [s.robot.x(2), xi(2)]);
            pause(s.delay);

            drawnow
        end

        function display(s)
            %Sensor.display Display status of sensor object
            %
            % S.display() displays the state of the sensor object in
            % human-readable form.
            %
            % Notes::
            % - This method is invoked implicitly at the command line when the result
            %   of an expression is a Sensor object and the command has no trailing
            %   semicolon.
            %
            % See also Sensor.char.
            loose = strcmp( get(0, 'FormatSpacing'), 'loose');
            if loose
                disp(' ');
            end
            disp([inputname(1), ' = '])
            disp( char(s) );
        end % display()

        function str = char(s)
            %Sensor.char Convert sensor parameters to a string
            %
            % s = S.char() is a string showing sensor parameters in
            % a compact human readable format.
            str = [class(s) ' sensor class:'];
            str = char(str, char(s.map));
        end

    end % method
end % classdef

author: Allopart

repo: rbpf-gmapping


text:
function [E,dE,f,g] = spm_DEM_eval(M,qu,qp)
% evaluates state equations and derivatives for DEM schemes
% FORMAT [E dE f g] = spm_DEM_eval(M,qu,qp)
%
% M  - model structure
% qu - conditional mode of states
%  qu.v{i} - casual states
%  qu.x(i) - hidden states
%  qu.y(i) - response
%  qu.u(i) - input
% qp - conditional density of parameters
%  qp.p{i} - parameter deviates for i-th level
%  qp.u(i) - basis set
%  qp.x(i) - expansion point ( = prior expectation)
%
% E  - generalised errors  (i.e.., y - g(x,v,P); x[1] - f(x,v,P))
%
% dE:
%  dE.du   - de[1:n]/du
%  dE.dy   - de[1:n]/dy[1:n]
%  dE.dc   - de[1:n]/dc[1:d]
%  dE.dp   - de[1:n]/dp
%  dE.dup  - d/dp[de[1:n]/du
%  dE.dpu  - d/du[de[1:n]/dp
%
% where u = x{1:d]; v[1:d]
%
% To accelerate computations one can specify the nature of the model using
% the field:
%
% M(1).E.linear = 0: full        - evaluates 1st and 2nd derivatives
% M(1).E.linear = 1: linear      - equations are linear in x and v
% M(1).E.linear = 2: bilinear    - equations are linear in x, v & x*v
% M(1).E.linear = 3: nonlinear   - equations are linear in x, v, x*v, & x*x
% M(1).E.linear = 4: full linear - evaluates 1st derivatives (for generalised
%                                  filtering, where parameters change)
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_DEM_eval.m 6270 2014-11-29 12:04:48Z karl $

% get dimensions
%==========================================================================
nl    = size(M,2);                       % number of levels
ne    = sum(spm_vec(M.l));               % number of e (errors)
nv    = sum(spm_vec(M.m));               % number of x (causal states)
nx    = sum(spm_vec(M.n));               % number of x (hidden states)
np    = sum(spm_vec(M.p));               % number of p (parameters)

% evaluate functions at each hierarchical level
%==========================================================================

% Get states {qu.v{1},qu.x{1}} in hierarchical form (v{i},x{i})
%--------------------------------------------------------------------------
v     = spm_unvec(qu.v{1},{M(1 + 1:end).v});
x     = spm_unvec(qu.x{1},{M(1:end - 1).x});
for i = 1:(nl - 1)
    p      = spm_unvec(spm_vec(M(i).pE) + qp.u{i}*qp.p{i},M(i).pE);
    f{i,1} = feval(M(i).f,x{i},v{i},p);
    g{i,1} = feval(M(i).g,x{i},v{i},p);
end

% Get Derivatives
%==========================================================================
persistent D
try
    method = M(1).E.linear;
catch
    method = 0;
end

switch method

    % get derivatives at each iteration of D-step - full evaluation
    %----------------------------------------------------------------------
    case{0}

        D     = spm_DEM_eval_diff(x,v,qp,M);

        % gradients w.r.t. states
        %------------------------------------------------------------------
        dedy  = D.dedy;
        dedc  = D.dedc;
        dfdy  = D.dfdy;
        dfdc  = D.dfdc;
        dgdx  = D.dgdx;
        dgdv  = D.dgdv;
        dfdv  = D.dfdv;
        dfdx  = D.dfdx;
        dgdxp = D.dgdxp;
        dfdxp = D.dfdxp;
        dgdvp = D.dgdvp;
        dfdvp = D.dfdvp;

        % gradients w.r.t. parameters
        %------------------------------------------------------------------
        dgdp  = D.dgdp;
        dfdp  = D.dfdp;

    % linear: assume equations are linear in x and v
    %----------------------------------------------------------------------
    case{1}

        % get derivatives and store expansion point (states)
        %------------------------------------------------------------------
        if isempty(D)

            D     = spm_DEM_eval_diff(x,v,qp,M);
            D.x   = x;
            D.v   = v;

            % gradients w.r.t. states
            %--------------------------------------------------------------
            dedy  = D.dedy;
            dedc  = D.dedc;
            dfdy  = D.dfdy;
            dfdc  = D.dfdc;
            dgdx  = D.dgdx;
            dgdv  = D.dgdv;
            dfdv  = D.dfdv;
            dfdx  = D.dfdx;

            % gradients w.r.t. parameters (state-dependent)
            %--------------------------------------------------------------
            dgdxp = D.dgdxp;
            dfdxp = D.dfdxp;
            dgdvp = D.dgdvp;
            dfdvp = D.dfdvp;

            % gradients w.r.t. parameters
            %--------------------------------------------------------------
            dgdp  = D.dgdp;
            dfdp  = D.dfdp;

        % linear expansion for derivatives w.r.t. parameters
        %------------------------------------------------------------------
        else

            % gradients w.r.t. states
            %--------------------------------------------------------------
            dedy  = D.dedy;
            dedc  = D.dedc;
            dfdy  = D.dfdy;
            dfdc  = D.dfdc;
            dgdx  = D.dgdx;
            dgdv  = D.dgdv;
            dfdv  = D.dfdv;
            dfdx  = D.dfdx;
            dgdxp = D.dgdxp;
            dfdxp = D.dfdxp;
            dgdvp = D.dgdvp;
            dfdvp = D.dfdvp;

            % gradients w.r.t. parameters
            %--------------------------------------------------------------
            dx    = spm_vec(qu.x{1}) - spm_vec(D.x);
            dv    = spm_vec(qu.v{1}) - spm_vec(D.v);
            dgdp  = D.dgdp;
            dfdp  = D.dfdp;
            for p = 1:np
                dgdp(:,p) = D.dgdp(:,p) + D.dgdxp{p}*dx + D.dgdvp{p}*dv;
                if nx
                    dfdp(:,p) = D.dfdp(:,p) + D.dfdxp{p}*dx + D.dfdvp{p}*dv;
                end
            end

        end

    % bilinear: assume equations are linear in x and v and x*v
    %----------------------------------------------------------------------
    case{2}

        % get derivatives and store expansion point (states)
        %------------------------------------------------------------------
        if isempty(D)

            % get high-order derivatives
            %--------------------------------------------------------------
            [Dv D] = spm_diff('spm_DEM_eval_diff',x,v,qp,M,2);

            for i = 1:nv, Dv{i} = spm_unvec(Dv{i},D); end
            D.x   = x;
            D.v   = v;
            D.Dv  = Dv;

            % gradients w.r.t. states
            %--------------------------------------------------------------
            dedy  = D.dedy;
            dedc  = D.dedc;
            dfdy  = D.dfdy;
            dfdc  = D.dfdc;
            dgdx  = D.dgdx;
            dgdv  = D.dgdv;
            dfdv  = D.dfdv;
            dfdx  = D.dfdx;
            dgdxp = D.dgdxp;
            dfdxp = D.dfdxp;
            dgdvp = D.dgdvp;
            dfdvp = D.dfdvp;

            % gradients w.r.t. parameters
            %--------------------------------------------------------------
            dgdp  = D.dgdp;
            dfdp  = D.dfdp;

        % linear expansion for derivatives w.r.t. parameters
        %------------------------------------------------------------------
        else

            % gradients w.r.t. causes and data
            %--------------------------------------------------------------
            dedy  = D.dedy;
            dedc  = D.dedc;
            dfdy  = D.dfdy;
            dfdc  = D.dfdc;

            % states (relative to expansion point)
            %--------------------------------------------------------------
            dv    = spm_vec(qu.v{1}) - spm_vec(D.v);

            % gradients w.r.t. states
            %--------------------------------------------------------------
            dgdx  = D.dgdx;
            dgdv  = D.dgdv;
            dfdx  = D.dfdx;
            dfdv  = D.dfdv;
            for i = 1:nv; dgdx = dgdx + D.Dv{i}.dgdx*dv(i); end
            for i = 1:nv; dgdv = dgdv + D.Dv{i}.dgdv*dv(i); end
            for i = 1:nv; dfdx = dfdx + D.Dv{i}.dfdx*dv(i); end
            for i = 1:nv; dfdv = dfdv + D.Dv{i}.dfdv*dv(i); end

            % second-order derivatives
            %--------------------------------------------------------------
            dgdxp = D.dgdxp;
            dgdvp = D.dgdvp;
            dfdxp = D.dfdxp;
            dfdvp = D.dfdvp;
            for p = 1:np
                for i = 1:nv; dgdxp{p} = dgdxp{p} + D.Dv{i}.dgdxp{p}*dv(i); end
                for i = 1:nv; dgdvp{p} = dgdvp{p} + D.Dv{i}.dgdvp{p}*dv(i); end
                for i = 1:nv; dfdxp{p} = dfdxp{p} + D.Dv{i}.dfdxp{p}*dv(i); end
                for i = 1:nv; dfdvp{p} = dfdvp{p} + D.Dv{i}.dfdvp{p}*dv(i); end
            end

            % gradients w.r.t. parameters
            %--------------------------------------------------------------
            dgdp  = D.dgdp;
            dfdp  = D.dfdp;
            for p = 1:np
                Dgdxp     = (D.dgdxp{p} + dgdxp{p})/2;
                Dgdvp     = (D.dgdvp{p} + dgdvp{p})/2;
                Dfdxp     = (D.dfdxp{p} + dfdxp{p})/2;
                Dfdvp     = (D.dfdvp{p} + dfdvp{p})/2;
                dgdp(:,p) = dgdp(:,p) + Dgdvp*dv;
                dfdp(:,p) = dfdp(:,p) + Dfdvp*dv;
            end

        end

    % nonlinear: assume equations are bilinear in x and v
    %----------------------------------------------------------------------
    case{3}

        % get derivatives and store expansion point (states)
        %------------------------------------------------------------------
        if isempty(D)

            % get high-order derivatives
            %--------------------------------------------------------------
            [Dx D] = spm_diff('spm_DEM_eval_diff',x,v,qp,M,1,'q');
            [Dv D] = spm_diff('spm_DEM_eval_diff',x,v,qp,M,2,'q');

            for i = 1:nx, Dx{i} = spm_unvec(Dx{i},D); end
            for i = 1:nv, Dv{i} = spm_unvec(Dv{i},D); end
            D.x   = x;
            D.v   = v;
            D.Dx  = Dx;
            D.Dv  = Dv;

            % gradients w.r.t. states
            %--------------------------------------------------------------
            dedy  = D.dedy;
            dedc  = D.dedc;
            dfdy  = D.dfdy;
            dfdc  = D.dfdc;
            dgdx  = D.dgdx;
            dgdv  = D.dgdv;
            dfdv  = D.dfdv;
            dfdx  = D.dfdx;
            dgdxp = D.dgdxp;
            dfdxp = D.dfdxp;
            dgdvp = D.dgdvp;
            dfdvp = D.dfdvp;

            % gradients w.r.t. parameters
            %--------------------------------------------------------------
            dgdp  = D.dgdp;
            dfdp  = D.dfdp;

        % linear expansion for derivatives w.r.t. parameters
        %------------------------------------------------------------------
        else

            % gradients w.r.t. causes and data
            %--------------------------------------------------------------
            dedy  = D.dedy;
            dedc  = D.dedc;
            dfdy  = D.dfdy;
            dfdc  = D.dfdc;

            % states (relative to expansion point)
            %--------------------------------------------------------------
            dx    = spm_vec(qu.x{1}) - spm_vec(D.x);
            dv    = spm_vec(qu.v{1}) - spm_vec(D.v);

            % gradients w.r.t. states
            %--------------------------------------------------------------
            dgdx  = D.dgdx;
            dgdv  = D.dgdv;
            dfdx  = D.dfdx;
            dfdv  = D.dfdv;
            for i = 1:nx; dgdx = dgdx + D.Dx{i}.dgdx*dx(i); end
            for i = 1:nv; dgdx = dgdx + D.Dv{i}.dgdx*dv(i); end
            for i = 1:nx; dgdv = dgdv + D.Dx{i}.dgdv*dx(i); end
            for i = 1:nv; dgdv = dgdv + D.Dv{i}.dgdv*dv(i); end
            for i = 1:nx; dfdx = dfdx + D.Dx{i}.dfdx*dx(i); end
            for i = 1:nv; dfdx = dfdx + D.Dv{i}.dfdx*dv(i); end
            for i = 1:nx; dfdv = dfdv + D.Dx{i}.dfdv*dx(i); end
            for i = 1:nv; dfdv = dfdv + D.Dv{i}.dfdv*dv(i); end

            % second-order derivatives
            %--------------------------------------------------------------
            dgdxp = D.dgdxp;
            dgdvp = D.dgdvp;
            dfdxp = D.dfdxp;
            dfdvp = D.dfdvp;
            for p = 1:np
                for i = 1:nx; dgdxp{p} = dgdxp{p} + D.Dx{i}.dgdxp{p}*dx(i); end
                for i = 1:nv; dgdxp{p} = dgdxp{p} + D.Dv{i}.dgdxp{p}*dv(i); end
                for i = 1:nx; dgdvp{p} = dgdvp{p} + D.Dx{i}.dgdvp{p}*dx(i); end
                for i = 1:nv; dgdvp{p} = dgdvp{p} + D.Dv{i}.dgdvp{p}*dv(i); end
                for i = 1:nx; dfdxp{p} = dfdxp{p} + D.Dx{i}.dfdxp{p}*dx(i); end
                for i = 1:nv; dfdxp{p} = dfdxp{p} + D.Dv{i}.dfdxp{p}*dv(i); end
                for i = 1:nx; dfdvp{p} = dfdvp{p} + D.Dx{i}.dfdvp{p}*dx(i); end
                for i = 1:nv; dfdvp{p} = dfdvp{p} + D.Dv{i}.dfdvp{p}*dv(i); end
            end

            % gradients w.r.t. parameters
            %--------------------------------------------------------------
            dgdp  = D.dgdp;
            dfdp  = D.dfdp;
            for p = 1:np
                Dgdxp     = (D.dgdxp{p} + dgdxp{p})/2;
                Dgdvp     = (D.dgdvp{p} + dgdvp{p})/2;
                Dfdxp     = (D.dfdxp{p} + dfdxp{p})/2;
                Dfdvp     = (D.dfdvp{p} + dfdvp{p})/2;
                dgdp(:,p) = dgdp(:,p) + Dgdxp*dx + Dgdvp*dv;
                dfdp(:,p) = dfdp(:,p) + Dfdxp*dx + Dfdvp*dv;
            end

        end

    % repeated evaluation of first order derivatives (for Laplace scheme)
    %----------------------------------------------------------------------
    case{4}

        % get derivatives and store expansion point (states)
        %------------------------------------------------------------------
        if isempty(D)

            D     = spm_DEM_eval_diff(x,v,qp,M);
            D.x   = x;
            D.v   = v;

            % gradients w.r.t. states
            %--------------------------------------------------------------
            dedy  = D.dedy;
            dedc  = D.dedc;
            dfdy  = D.dfdy;
            dfdc  = D.dfdc;
            dgdx  = D.dgdx;
            dgdv  = D.dgdv;
            dfdv  = D.dfdv;
            dfdx  = D.dfdx;

            % gradients w.r.t. parameters (state-dependent)
            %--------------------------------------------------------------
            dgdxp = D.dgdxp;
            dfdxp = D.dfdxp;
            dgdvp = D.dgdvp;
            dfdvp = D.dfdvp;

            % gradients w.r.t. parameters
            %--------------------------------------------------------------
            dgdp  = D.dgdp;
            dfdp  = D.dfdp;

        % re-evaluate first-order derivatives
        %------------------------------------------------------------------
        else

            % retain second-order gradients
            %--------------------------------------------------------------
            dgdxp = D.dgdxp;
            dfdxp = D.dfdxp;
            dgdvp = D.dgdvp;
            dfdvp = D.dfdvp;

            % re-evaluate first-order gradients
            %--------------------------------------------------------------
            D     = spm_DEM_eval_diff(x,v,qp,M,0);
            dedy  = D.dedy;
            dedc  = D.dedc;
            dfdy  = D.dfdy;
            dfdc  = D.dfdc;
            dgdx  = D.dgdx;
            dgdv  = D.dgdv;
            dfdv  = D.dfdv;
            dfdx  = D.dfdx;

            % replace second-order gradients
            %--------------------------------------------------------------
            D.dgdxp = dgdxp;
            D.dfdxp = dfdxp;
            D.dgdvp = dgdvp;
            D.dfdvp = dfdvp;

            % gradients w.r.t. parameters
            %--------------------------------------------------------------
            dx    = spm_vec(qu.x{1}) - spm_vec(x);
            dv    = spm_vec(qu.v{1}) - spm_vec(v);
            dgdp  = D.dgdp;
            dfdp  = D.dfdp;
            for p = 1:np
                dgdp(:,p) = D.dgdp(:,p) + D.dgdxp{p}*dx + D.dgdvp{p}*dv;
                if nx
                    dfdp(:,p) = D.dfdp(:,p) + D.dfdxp{p}*dx + D.dfdvp{p}*dv;
                end
            end

        end

    otherwise
        disp('Unknown method')

end

% order parameters (d = n = 1 for static models)
%--------------------------------------------------------------------------
d       = M(1).E.d + 1;                    % generalisation order of q(v)
n       = M(1).E.n + 1;                    % embedding order     (n >= d)

% Generalised prediction errors and derivatives
%==========================================================================
Ex      = cell(n,1);
Ev      = cell(n,1);
[Ex{:}] = deal(sparse(nx,1));
[Ev{:}] = deal(sparse(ne,1));

% prediction error (E) - causes
%--------------------------------------------------------------------------
for i = 1:n
    qu.y{i} = spm_vec(qu.y{i});
end
Ev{1} = [qu.y{1}; qu.v{1}] - [spm_vec(g); qu.u{1}];
for i = 2:n
    Ev{i} = dedy*qu.y{i} + dedc*qu.u{i} ...      % generalised response
          - dgdx*qu.x{i} - dgdv*qu.v{i};         % and prediction
end

% prediction error (E) - states
%--------------------------------------------------------------------------
try
    Ex{1} = qu.x{2} - spm_vec(f);
end
for i = 2:n - 1
    Ex{i} = qu.x{i + 1} ...                      % generalised motion
          - dfdx*qu.x{i} - dfdv*qu.v{i};         % and prediction
end

% error
%--------------------------------------------------------------------------
E     = spm_vec({Ev,Ex});

% Kronecker forms of derivatives for generalised motion
%==========================================================================
if nargout < 2, return, end

% dE.dp (parameters)
%--------------------------------------------------------------------------
dgdp  = {dgdp};
dfdp  = {dfdp};
for i = 2:n
    dgdp{i,1} = dgdp{1};
    dfdp{i,1} = dfdp{1};
    for p = 1:np
        dgdp{i,1}(:,p) = dgdxp{p}*qu.x{i} + dgdvp{p}*qu.v{i};
        dfdp{i,1}(:,p) = dfdxp{p}*qu.x{i} + dfdvp{p}*qu.v{i};
    end
end

% generalised temporal derivatives: dE.du (states)
%--------------------------------------------------------------------------
dedy  = kron(spm_speye(n,n),dedy);
dedc  = kron(spm_speye(n,d),dedc);
dfdy  = kron(spm_speye(n,n),dfdy);
dfdc  = kron(spm_speye(n,d),dfdc);
dgdx  = kron(spm_speye(n,n),dgdx);
dgdv  = kron(spm_speye(n,d),dgdv);
dfdv  = kron(spm_speye(n,d),dfdv);
dfdx  = kron(spm_speye(n,n),dfdx) - kron(spm_speye(n,n,1),speye(nx,nx));

% 1st error derivatives (states)
%--------------------------------------------------------------------------
dE.dy =  spm_cat({dedy; dfdy});
dE.dc =  spm_cat({dedc; dfdc});
dE.dp = -spm_cat({dgdp; dfdp});
dE.du = -spm_cat({dgdx, dgdv  ;
                  dfdx, dfdv});

% bilinear derivatives
%--------------------------------------------------------------------------
for i = 1:np
    dgdxp{i}  = kron(spm_speye(n,n),dgdxp{i});
    dfdxp{i}  = kron(spm_speye(n,n),dfdxp{i});
    dgdvp{i}  = kron(spm_speye(n,d),dgdvp{i});
    dfdvp{i}  = kron(spm_speye(n,d),dfdvp{i});
    dE.dup{i} = -spm_cat({dgdxp{i}, dgdvp{i};
                          dfdxp{i}, dfdvp{i}});
end
if np
    dE.dpu    =  spm_cell_swap(dE.dup);
else
    dE.dpu    =  {};
end

author: spm


text:
% Source Code Listings for the Calspan Version
% of the FlOl Engine Response Model
% Based on Report DNA-TR-93-124 by Adam J. Baran
% and Dr. Michael G. Dunn of Calspan Corporation
%
% AFRL/RQTI Sept 2012
%
% engine parameters:
% xn1 - fan speed (rpm)
% xn2 - core speed (rpm)
%    (fan speed correlation)
% xigv - fan inlet variable guide vane angle (deg.)
% psl4 - fan discharge static pressure (psig)
% pt25 - fan discharge total pressure (psig)
% xoil1 - oil pressure [xn1 correlation) (psig)
%    (core speed correlation)
% thrust - thrust (lbf)
% air - airflow (lbm/sec)
% vgv -compressor inlet variable guide vane angle (deg.)
% cdp - compressor discharge pressure Cpsig)
% pyro - pyrometer temperature (F)
% egt - exhaust gas temperature (F)
% fuel - fuel flow (Pph)
% xoil2c - oil pressure [xn2 correlation] (psig)
%
%       (subscripts)
% ()nom - nominal value
% ()del - Cin/de)crement over current interval
% ()sumc - sum of (in/de)crements over previous intervals
%        (prescripts)
% f() - function which generates the value of the variable
%
% Normal syntax to run is "Dunnv1('infile.dat','outfile.dat')" from matlab
% prompt. To build datafile on first run just run "Dunnv1"...
%% Main Routine
function Dunnv1(Infile,Outfile)
dbstop if error
% Check for whether a runtime datafile exists
if nargin == 0 ;            % no formatted input file exists
    defaultdata ;
    Infile = 'infile.dat';
    Outfile = 'outfile.dat';
elseif nargin == 1 ;        % formatted input file exists
    Outfile = 'outfile.dat';
end ;
STOR = [];
%% -------------------------------
% ROUTINE - READIN DATA FILE
%  -------------------------------
% readindata(Infile) ; Read in data from inputfile ; IN = fopen(Infile,'r') ;
%if (Infile == '');
%%%    Infile = 'infile.dat' ; %%% Remove this when converting to input
%end;
IN = fopen(Infile); %IN = fopen('infile.dat','w')
[pla] = textread(Infile, '%f',1,'headerlines',4) ;
[pamb] = textread(Infile, '%f',1,'headerlines',6) ;
[tamb] = textread(Infile, '%f',1,'headerlines',8) ;
[debug] = textread(Infile, '%u',1,'headerlines',10) ;
[nenc] = textread(Infile, '%u',1,'headerlines',12) ;
% Preallocate the arrays
% ib = zeros(nenc,1); blend = zeros(nenc,1); c = zeros(nenc,1); delt = zeros(nenc,1);
[ib blend c delt] =  textread(Infile,'%u %u %f %f','headerlines',14) ;
fclose('all') ;
%% Start writting out Header Block of Output File
%% -------------------------------
% ROUTINE - WRITEOUT RESULTS
%  -------------------------------
%if (Infile == '');
%%%    Outfile = 'outfile.dat' ; %%% Remove this when converting to input
%end;
OUT = fopen(Outfile,'w') ;
fprintf(OUT,'PEARL - F101 reponse model \n') ;
fprintf(OUT,'%s\n',datestr(now)) ;
fprintf(OUT,' \n') ; fprintf(OUT,' \n') ; fprintf(OUT,' \n') ;
fprintf(OUT,'inputs :  \n'); fprintf(OUT,' \n') ;
fprintf(OUT,'PLA (deg) from (16.5-80.0) =>  \n') ;
fprintf(OUT,'%-4.2e  \n',pla) ;
fprintf(OUT,'ambient pressure (psia) => ') ; fprintf(OUT,'%-4.2e  \n',pamb) ;
fprintf(OUT,'ambient temperature (F) => ') ; fprintf(OUT,'%-4.2e  \n',tamb) ;
fprintf(OUT,' \n') ;
fprintf(OUT,'#_of_encounters(-)(1-10)=> ') ; fprintf(OUT,'%-2i  \n',nenc) ;
fprintf(OUT,'enc#, blend (1=mpb,2=blend2), dust conc. (mg/m^3), duration (min) \n') ;
fprintf(OUT,' \n') ;
for i = 1:1:max(ib) ; % 'max(ib) is 'nenc'
    fprintf(OUT,'%-2i \t %-1i \t %-4.2e \t  %-4.2e \n', ib(i), blend(i),c(i),delt(i));
end ;
fprintf(OUT,' \n') ;
fprintf(OUT,'debug mode (0=off,1=on)  => ') ; fprintf(OUT,'%-2i  \n',debug) ;
fprintf(OUT,' \n') ; fprintf(OUT,' \n') ;
fprintf(OUT,'calculation begun   \n') ; fprintf(OUT,' \n'); fprintf(OUT,' \n');
fprintf(OUT,' i blend dust cnc. time int. AIR FLOW   ACC DUST    N2     CDP    FUEL FLOW   THRUST  \n');
fprintf(OUT,'         (mg/m^3)  (min)     lbm/sec)  (kg)        (rpm)  (psig) (lbm/hr)    (lbf)   \n');
fprintf(OUT,'+--+---+--------+--------+-------------+----------+-------+------+----------+-------+ \n');
%% START MAIN LOGIC HERE
% Initialize all Runtime Engine Parameters to zero
ienc = 0; none = 0; dust = 0.; icdpflag = 0; cdplim = 54.; accdust = 0. ;
ibug = 0; time = 0. ;
xn2nom = 0.; xn2del = 0.; xn2sum = 0.;            % Spool speeds
cdpnom = 0.; cdpdel = 0.; cdpsum = 0.;            % Compressor pressure
fuelnom = 0.; fueldel = 0.; fuelsum = 0.;         % Fuel mass flow
thrustnom = 0.; thrustdel = 0.; thrustsum = 0.;   % Engine Thrust
%%% F101 Specific DATA Here
const = 1.699e-06;                                % F101 engine correl. parameters for dust accum in F101
cdprate1 = 9.36e-03; cdprate2 = 1.59e-02;         % F101 engine correl. parameters for cdp w/ dust1 and 2
fuelrate1 = 1.22e-01; fuelrate2 = 3.32e-01;       % F101 engine correl. parameters for fuel mass flow w/ dust1 and 2
thrustrate1 = 1.41e-01; thrustrate2 = 3.33e-01;   % F101 engine correl. parameters for thrust w/ dust1 and 2
%%% F101 Specific DATA Here
% Calculate Ambient Density of Air (lbm/sec)
rhoamb = 144.*pamb/(53.3*(tamb+460.));            % Air equation of state
% Nominal clear air engine parameter values
xn2nom = FXN2NOM(pla) ;     % call function "figure nominal core speed from PLA"
xn2  = xn2nom ;             % set core speed from nominal core speed
air = FAIR(xn2) ;           % call function for "obtain air mass flow from core speed"
cdpnom = FCDPNOM(xn2) ;     % call function "get nominal compress discharge pressure from core speed"
cdp = cdpnom ;              % set comp. discharge pressure from nominal comp. discharge pressure
fuelnom = FFUELNOM(xn2) ;   % call function "get nominal fuel flow from core speed"
fuel = fuelnom ;            % set fuel flow from the nominal fuel flow
thrustnom = FTHRUSTNOM(xn2) ;
thrust = thrustnom ;
fprintf(OUT,'%2i\t%2i\t%4.1f\t%5.1f\t\t%4.1f\t\t%3.2e\t%6.1f\t%4.1f\t%6.1f\t%6.1f\n',0,0,0,0,air,accdust,xn2,cdp,fuel,thrust) ;
STORBASE = [0 0 0 0 0 air accdust xn2 cdp fuel thrust];
% LOOP over all dust clound encounters
for ienc = 1:1:nenc;
    %% Calculate Airflow, ingested dust over this interval and total accumulated dust
    % based on conditions at start of encounter
    air = FAIR(xn2);
    q = air/rhoamb ;
    dust = const*q*c(ienc)*delt(ienc) ; % Incremental Dust accumulation
    accdust = dust + accdust;           % Cumulative Dust accumulation
    %% Calculate Core Speed          %ibug %ib(ienc)%accdust %c(ienc)
    xn2rate = FXN2RATE(ibug,blend(ienc),accdust,c(ienc));
    xn2del = xn2rate*delt(ienc);
    xn2 = xn2nom - xn2del - xn2sum ;
    %% Calculate cdp
    cdpnom = FCDPNOM(xn2);
    if(blend(ienc) == 1);    % If dust blend #1 (mpb)
       cdprate = cdprate1 ;
    elseif(blend(ienc) == 2); % If dust blend #2 (blend-2)
       cdprate = cdprate2 ;
    end;
    cdpdel = cdprate*c(ienc)*delt(ienc);
    cdp = cdpnom + cdpdel + cdpsum ;
    %% Check for surge probability
    if((cdpdel+cdpsum >= cdplim)&&(icdpflag == 0));
       icdpflag = 1;
       cdpdelc = cdplim - cdpsum ;
       tau = cdpdelc/(cdprate*c(ienc));
       fprintf(OUT,'%s %3.1f  %s  %2i %s','*** engine surge at ',tau,'minutes into the',ienc,' -th encounter');
       fprintf(OUT,' \n') ;
    end
    %% Calculate Fuel Burn Rate
    fuelnom = FFUELNOM(xn2) ;
    if(blend(ienc) == 1);    % If dust blend #1 (mpb)
       fuelrate = fuelrate1 ;
    elseif(blend(ienc) == 2); % If dust blend #2 (blend-2)
       fuelrate = fuelrate2 ;
    end;
    fueldel = fuelrate*c(ienc)*delt(ienc) ;
    fuel = fuelnom + fueldel + fuelsum ;
    %% Calculate Thrust
    thrustnom = FTHRUSTNOM(xn2) ;
    if(blend(ienc) == 1);    % If dust blend #1 (mpb)
       thrustrate = thrustrate1 ;
    elseif(blend(ienc) == 2); % If dust blend #2 (blend-2)
       thrustrate = thrustrate2 ;
    end;
    thrustdel = thrustrate*c(ienc)*delt(ienc);
    thrust = thrustnom + thrustdel + thrustsum ;
    %% Updat the sum terms of each engine parameter
    xn2sum = xn2sum + xn2del ;
    cdpsum = cdpsum + cdpdel ;
    fuelsum = fuelsum + fueldel ;
    thrustsum = thrustsum + thrustdel ;
    %% Write and store out runtime results
    fprintf(OUT,'%2i\t%2i\t%4.1f\t%5.1f\t\t%4.1f\t\t%3.2e\t%6.1f\t%4.1f\t%6.1f\t%6.1f\n',ienc,blend(ienc),c(ienc),delt(ienc),air,accdust,xn2,cdp,fuel,thrust) ;
    time = time + delt(ienc);
    STOR(ienc,:) = [ienc time blend(ienc) c(ienc) delt(ienc) air accdust xn2 cdp fuel thrust];
end
STOR = [STORBASE ; STOR] ;
% Plot Out Trends
figure(1); scale = 0;
subplot(2,3,1); plot(STOR(:,2),STOR(:,4)); grid on ; axis([ 0 max(STOR(:,2)) 0 1.1.*max(STOR(:,4))]) ;
title('dust concentration at enc') ; xlabel('time (min)'); ylabel('concentration (mg/m^3)') ;
subplot(2,3,2); plot(STOR(:,2),STOR(:,6)); grid on ;
if (scale == 1); axis([ 0 max(STOR(:,2)) 0 1.1.*max(STOR(:,6))]) ; end;
title('air flow') ; xlabel('time (min)'); ylabel('air mass flow (lbm/sec)') ;
subplot(2,3,3); plot(STOR(:,2),STOR(:,8)); grid on ;
if (scale == 1); axis([ 0 max(STOR(:,2)) 0 1.1.*max(STOR(:,8))]) ; end;
title('fan speed') ; xlabel('time (min)'); ylabel('fan speed (rpm)') ;
subplot(2,3,4); plot(STOR(:,2),STOR(:,9)); grid on ;
if (scale == 1); axis([ 0 max(STOR(:,2)) 0 1.1.*max(STOR(:,9))]) ; end ;
title('compressor discharge pressure') ; xlabel('time (min)'); ylabel('pressure (psig)') ;
subplot(2,3,5); plot(STOR(:,2),STOR(:,10)); grid on ;
if (scale == 1); axis([ 0 max(STOR(:,2)) 0 1.1.*max(STOR(:,10))]) ; end;
title('fuel burn') ; xlabel('time (min)'); ylabel('fuel burn rate (pph)') ;
subplot(2,3,6); plot(STOR(:,2),STOR(:,11)); grid on ;
if (scale == 1); axis([ 0 max(STOR(:,2)) 0 1.1.*max(STOR(:,11))]) ; end;
title('thrust') ; xlabel('time (min)'); ylabel('thrust (lbf)') ;
%
fclose(OUT) ;
% END of MAIN Routine
end
%% Calculate Airflow as a function a fan of fan speed
% "x" is passed from "xn2" and "y" is passed from "air"
function [y] = FAIR(x);  %%% F101 Specific DATA Here
%%% F101 Specific DATA Here
a0= 440922.79744; a1= -174.0731208; a2= 0.027346652504;
a3= -2.1376216416e-06; a4= 8.3178976743e-11; a5= -1.2891058570e-15;
p = [440922.79744 -174.0731208 0.027346652504 -2.1376216416e-06 8.3178976743e-11];
%%% F101 Specific DATA Here
if(x<10348.);
    y=0.;
elseif(x>15000.);
    y=0.;
elseif((x>=10348.)&&(x<=15000.));
    y = a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4 + a5*x^5;
%    y = polyval(p,x)
end
end
%% Calculate the nominal fuel a fan of fan speed, x is passed from xn2
function [y] = FFUELNOM(x) ;%%% F101 Specific DATA Here
%%% F101 Specific DATA Here
a0 = -49816.415015; a1 = 15.16955213; a2 = -0.00090483247417;
a3 = -1.0788972042E-07; a4 = 1.3439057313e-11; a5 = -3.6898802700E-16;
%%% F101 Specific DATA Here
if(x<10348.);
    y = 0. ;
elseif(x>15000.);
    y = 0. ;
elseif((x>=10348.)&&(x<=15000.));
    y = a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4 + a5*x^5;
end
end
%% Function calculates the nominal thrust a fan of fan speed
function [y] = FTHRUSTNOM(x);%%% F101 Specific DATA Here
% 'x' is 'xn2' while 'y' is 'fthrustnom'
%%% F101 Specific DATA Here
a10=-90190833.084; a11=38206.13731; a12=-6.4537917919; a13=0.00054340227729;
a14=-2.2807645464e-08; a15=3.8182035691e-13;
a20=-26525834.184; a21=5596.7412494; a22=-0.24258263565; a23=-2.2544154617e-05;
a24=2.2276804666e-09; a25=-5.2033103514e-14;
%%% F101 Specific DATA Here
if(x<10543.) ;
    y = 0. ;
elseif(x>15000.) ;
    y = 0. ;
elseif((x>=10543.)&&(x<=13246.));
    y = a10 + a11*x + a12*x^2 + a13*x^3 + a14*x^4 + a15*x^5;
elseif((x>13426.)&&(x<=15000.));     % else
    y = a20 + a21*x + a22*x^2 + a23*x^3 + a24*x^4 + a25*x^5;
end
end
%%  Calculate Nominal value of Core Speed as function of power setting
function [y] = FXN2NOM(x);%%% F101 Specific DATA Here
% "x" is passed from the global variable "pla", and "y" is "xn2nom"
%%% F101 Specific DATA Here
a10=10353.360465; a11=0.18067978533; a20=10306.641571;
a21=220.99533596; a22=-21.326092808; a23=0.69563112177;
a24=-0.0087987179757; a25=3.8820512919e-05;
%%% F101 Specific DATA Here
if(x<16.5) ;
    y = 0. ;
elseif(x>80.);
    y = 0. ;
elseif((x>=16.5)&&(x<30.));
    y = a10 + a11*x ;
elseif((x>=30.)&&(x<=80.));
    y = a20 + a21*x + a22*x^2 + a23*x^3 + a24*x^4 + a25*x^5 ;
end
end
%% Calculate Compressor Discharge Pressure from core speed
function [y]= FCDPNOM(x);%%% F101 Specific DATA Here
% "x" is passed from "xn2" and "y" is passed at "cdpnom"
%%% F101 Specific DATA Here
a0 = 52767.656799 ; a1 = -20.435530217 ; a2 = 0.0031594211659 ;
a3 = -2.441489353e-07 ; a4 = 9.4426267256e-12 ; a5 = -1.4592076510e-16 ;
%%% F101 Specific DATA Here
if(x<10348.);
    y = 0. ;
elseif(x>15000.);
    y = 0. ;
elseif((x>=10348.)&&(x<=15000.));
    y = a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4 + a5*x^5 ;
end
end
%% Calculates core speed drop rate due to dust
function [y] = FXN2RATE(libug,lblend,lacc,lc) ;%%% F101 Specific DATA Here
%% where 'libug' is passed from 'ibug', is debug mode
% where 'li' is passed from 'i', is cloud encounter index counter
% where 'lacc' is passed from 'acc', is sum accumulation of dust
% where 'lc' is passed from 'c', concentration of the dust cloud
% where 'y' is passed back as 'fxn2rate', the core speed drop rate
%%% F101 Specific DATA Here
a110 = 0. ; a111 = 0.031091666884 ; a120 = 2.2386 ; a130 = 0.41616348802 ;
a131= 0.062687011283 ; a210= -0.023120229598 ; a211= 0.19569725331 ;
a212= -0.001954960762 ; a213= 6.9569031307e-06;
%%% F101 Specific DATA Here
if(lacc<=0.); y = 0.0; end;
if(lblend==1);                  % Model for the mpb blend
    if(lacc<72.);
        if((lc>=0.)&&(lc<72.));
            y = a110 + a111*lc ;
        elseif((lc>=72.)&&(lc<=293.));
            y = a120 ;
        elseif(lc>293);
            y = 1.0e10 ;
        end;
    elseif(lacc>=72.);
        if((lc>=0.)&&(lc<=282.));
            y = a130 + a131*lc ;
        elseif(lc>282.);
            y = 1.0e10;
        end;
    end;
elseif(lblend==2);              % Model for Blend #2
    if((lc>=0.)&&(lc<=295.));
        y = a210 + a211*lc + a212*lc^2 + a213*lc^3 ;
    elseif(lc>295.);
        y = 1.0e10;
    end;
end
end
%% DEFAULT DATA INPUT FILE
function defaultdata() ;
%% --------------------------------
%       DEFAULT DATA INPUT FILE
%  --------------------------------
% Write out a formatted 'infile.dat' for user editing if none exists in the
% directory where 'Dunnv1.m' is run.
% Defualt INPUT VALUES BLOCK
pla = 80.0 ;    % 65 percent throttle (from 16.5 to 80)
tamb = 59.0 ;   % Ambient temperature (in degF)
pamb =  14.7 ;  % Ambient pressure (in psia)
nenc = 10 ;     % No. of dust cloud encounters (No. of enc)
% Dust cloud encounter index numbers
ib = [0 1 2 3 4 5 6 7 8 9 10];
% CMAS blend (1 = mpb, 2 = blend2), as defined in DNA TR-92-121
blend = [0 2 2 2 2 2 2 2 2 2 2];
% dust conc. (mg/m^3)
c = [0. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100.];
% duration (min)
delt = [0.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0];
% debug mode {0=off, 1=on}
debug = 0;
% Write out these properties into default file in carefully formatted
% output
IN = fopen('infile.dat','w') ;
fprintf(IN,'PEARL - F101 reponse model \n') ;
fprintf(IN,' \n') ;
fprintf(IN,' \n') ;
fprintf(IN,'PLA (deg) from (16.5-80.0) =>  \n') ;
fprintf(IN,'%-4.2e  \n',pla) ;
fprintf(IN,'ambient pressure (psia) =>  \n') ;
fprintf(IN,'%-4.2e  \n',pamb) ;
fprintf(IN,'ambient temperature (F) =>  \n') ;
fprintf(IN,'%-4.2e  \n',tamb) ;
fprintf(IN,'debug mode (0=off,1=on)  =>  \n') ;
fprintf(IN,'%-2i  \n',debug) ;
fprintf(IN,'#_of_encounters(-)(1-10)=>  \n') ;
fprintf(IN,'%-2i  \n',nenc) ;
fprintf(IN,'enc#, blend (1=mpb,2=blend2), dust conc. (mg/m^3), duration (min) \n') ;
fprintf(IN,' \n') ;
for i = 1:1:nenc; fprintf(IN,'%-2i \t %-1i \t %-4.2e \t  %-4.2e \n', ib(i+1),blend(i+1),c(i+1),delt(i+1));end;
fclose(IN); % Close file
end

author: Sable

repo: mcbench-benchmarks


text:
function vF = sigma

vF = S2VectorFieldHandle(@(r) local(r));

  function v = local(r)
    [theta,rho] = polar(r); %#ok<POLAR>
    rot = rotation.byEuler(rho,theta,-rho,'ZYZ');
    v = rot .* xvector;
  end

end

author: mtex-toolbox

repo: mtex

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


text:
%calling function to search the optimal number of clusters
close all
clear all
path(path,'..\..\FUZZCLUST')

%the data
load motorcycle.txt
data.X=motorcycle(:,[1 2]);
[N,n]=size(data.X);
%data normalizaiton
data = clust_normalize(data,'range');

%parameters
ncmax=14; %maximal number of cluster
param.m=2;
param.e=1e-3;
%
ment=[];
figure(1)
for cln=2:ncmax
param.c=cln;
    param.ro = ones(1,param.c);
    result=GKclust(data,param);
    clf
    plot(data.X(:,1),data.X(:,2),'b.',result.cluster.v(:,1),result.cluster.v(:,2),'r*');
    hold on
    new.X=data.X;
    clusteval(new,result,param)
    %validation
    result=modvalidity(result,data,param);
    ment{cln}=result.validity;

end

PC=[];CE=[];SC=[];S=[];XB=[];DI=[];ADI=[];

    for i=2:ncmax
       PC=[PC ment{i}.PC];
       CE=[CE ment{i}.CE];
       SC=[SC ment{i}.SC];
       S=[S ment{i}.S];
       XB=[XB ment{i}.XB];
       DI=[DI ment{i}.DI];
       ADI=[ADI ment{i}.ADI];
   end
    figure(2)
    clf
    subplot(2,1,1), plot([2:ncmax],PC)
    title('Partition Coefficient (PC)')
    subplot(2,1,2), plot([2:ncmax],CE,'r')
    title('Classification Entropy (CE)')
    figure(3)
    subplot(3,1,1), plot([2:ncmax],SC,'g')
    title('Partition Index (SC)')
    subplot(3,1,2), plot([2:ncmax],S,'m')
    title('Separation Index (S)')
    subplot(3,1,3), plot([2:ncmax],XB)
    title('Xie and Beni Index (XB)')
    figure(4)
    subplot(2,1,1), plot([2:ncmax],DI)
    title('Dunn Index (DI)')
    subplot(2,1,2), plot([2:ncmax],ADI)
    title('Alternative Dunn Index (ADI)')


text:
function disp(X,name)
%DISP Command window display of a sparse tensor.
%
%   DISP(X) displays the tensor without printing its name.
%
%   DISP(X,NAME) displays the tensor with the given name.
%
%   See also SPTENSOR, SPTENSOR/DISPLAY.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

% Extract the number of nonzeros and number of dimensions
nz = nnz(X);

if ~exist('name','var')
    name = 'ans';
end

if (nz == 0)
    fprintf('%s is an all-zero sparse tensor of size %s\n',...
        name, tt_size2str(X.size));
    return;
else
    fprintf('%s is a sparse tensor of size %s with %d nonzeros\n',...
        name, tt_size2str(X.size), nz);
end

% Stop insane printouts
if (nz > 10000)
    r = input('Are you sure you want to print all nonzeros? (Y/N) ','s');
    if upper(r) ~= 'Y', return, end;
end

% preallocate
output = cell(nz,1);
%%
spc = floor(log10(max(double(X.subs),[],1)))+1;
if numel(spc) == 1
    fmt = ['\t(%' num2str(spc(1)) 'd)%s'];
else
    fmt = ['\t(%' num2str(spc(1)) 'd,'];
    for i = 2:numel(spc)-1
        fmt = [fmt '%' num2str(spc(i)) 'd,'];
    end
    fmt = [fmt '%' num2str(spc(end)) 'd)%s'];
end
%%
% Get values out so that they look nice
savefmt = get(0,'FormatSpacing');
format compact
S = evalc('disp(X.vals)');
set(0,'FormatSpacing',savefmt)
S = textscan(S,'%s','delimiter','\n','whitespace','');
S = S{1};
if ~isempty(strfind(S{1},'*'))
    fprintf('%s\n',S{1});
    S = S(2:end);
end
%%
for i = 1:nz
    output{i} = sprintf(fmt,X.subs(i,:),S{i});
end
fprintf('%s\n',output{:});

%     function y = fmt(s,v)
%         % nested function has access to n from disp function workspace
%         if n > 1
%             y = [sprintf('\t(') sprintf('%d,',s(1:n-1))...
%                 sprintf('%d) ',s(n)) sprintf('\t%g',v)];
%         else
%             y = sprintf('\t(%d) \t%f',s,v);
%         end
%     end

end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
%% Gaussian Mixture Model (GMM)
%
% Demonstrates EM clustering, and also compares againt K-means clustering.
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/python/gaussian_mix.py>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/python/kmeans.py>
%

function gaussian_mix_demo()
    %% Data
    % sample 2D points from a mixture distribution of K bivariate Gaussians
    K = 5;
    sz = 512;
    [pts, labels, mus, sigmas] = make_gaussian_mixture(K, sz);
    whos pts labels mus sigmas

    %%
    % draw points (color-coded) with the ground-truth mixtures
    clr = uint8(hsv(K) * 255);
    clr(:,4) = 0;
    img = zeros([sz sz 3], 'uint8');
    img0 = cv.circle(img, pts, 1, 'Colors',clr(labels,:), 'Thickness','Filled');
    for k=1:K
        img0 = draw_gaussian(img0, mus(k,:), sigmas{k}, [0 255 255]);
    end
    imshow(img0), title(sprintf('%d Samples, K=%d', numel(labels), K))

    %% EM
    % cluster using EM by fitting a Gaussian mixture model
    em = cv.EM();
    em.ClustersNumber = K;
    em.CovarianceMatrixType = 'Generic';
    tic
    [~,L] = em.trainEM(pts);
    toc

    %%
    % draw clustered points along with the estimated mixture model
    means = em.getMeans();
    covs = em.getCovs();
    img1 = cv.circle(img, pts, 1, 'Colors',clr(L+1,:), 'Thickness','Filled');
    for k=1:K
        img1 = draw_gaussian(img1, means(k,:), covs{k}, [0 255 255]);
    end
    figure, imshow(img1), title('Gaussian Mixtures')

    %% Kmeans
    % cluster using Kmeans
    crit = struct('type','Count+EPS', 'maxCount',30, 'epsilon',0.1);
    L = cv.kmeans(pts, K, 'Criteria',crit, 'Attempts',10);

    %%
    % draw points color-coded by assigned label
    img2 = cv.circle(img, pts, 1, 'Colors',clr(L+1,:), 'Thickness','Filled');
    figure, imshow(img2), title('Kmeans')
end

%% Helper functions

function [pts, labels, mus, sigmas] = make_gaussian_mixture(K, sz)
    %MAKE_GAUSSIAN_MIXTURE  Random points from Gaussian mixture distribution
    %
    %     [pts, labels, mus, sigmas] = make_gaussian_mixture(K, sz)
    %
    % ## Input
    % * __K__ number of components
    % * __sz__ image size, determines 2D points domain
    %
    % ## Output
    % * __pts__ matrix of 2D points
    % * __labels__ vector of corresponding component indices
    % * __mus__ K-by-2 mean of each component
    % * __sigmas__ cell array of length K, covariance of each component
    %
    % See also: gmdistribution.random
    %

    mus = zeros(K,2);
    sigmas = cell(K,1);
    pts = cell(K,1);
    labels = cell(K,1);
    for k=1:K
        mus(k,:) = (rand(1,2)*0.8 + 0.1) * sz;
        a = (rand(2,2) - 0.5) * sz * 0.1;
        sigmas{k} = (a.' * a) + sz*0.05*eye(2);
        n = 100 + randi([0 900]);
        pts{k} = my_mvnrnd(mus(k,:), sigmas{k}, n);
        labels{k} = k * ones(size(pts{k},1),1);
    end
    pts = single(cat(1, pts{:}));
    labels = int32(cat(1, labels{:}));
end

%%

function X = my_mvnrnd(mu, sigma, num)
    %MY_MVNRND  Random points from Gaussian distribution
    %
    %     X = my_mvnrnd(mu, sigma, num)
    %
    % ## Input
    % * __mu__ 1x2 mean vector
    % * __sigma__ 2x2 covariance matrix
    % * __num__ number of points to generate
    %
    % ## Output
    % * __X__ matrix of 2D points, num-by-2
    %
    % See also: mvnrnd, randn
    %

    if mexopencv.require('stats')
        X = mvnrnd(mu, sigma, num);
    else
        X = bsxfun(@plus, randn(num,numel(mu)) * cholcov(sigma), mu);
    end
end

%%

function img = draw_gaussian(img, mu, sigma, clr)
    %DRAW_GAUSSIAN  Draw a bivariate Gaussian
    %
    %     img = draw_gaussian(img, mu, sigma, clr)
    %
    % ## Input
    % * __img__ input image on which to draw
    % * __mu__ 1x2 mean vector
    % * __sigma__ 2x2 covariance matrix
    % * __clr__ color
    %
    % ## Output
    % * __img__ output image
    %
    % Represent a 2D Gaussian with an ellipse where:
    %
    % * its mean is the center of the ellipse
    % * its covariance matrix eigenvector corresponding to largest eigenvalue
    %   is the direction of the ellipse
    % * its (scaled) eigenvalues are the major/minor axis lengths of the
    %   ellipse
    %

    [w,u,~] = cv.SVD.Compute(sigma);
    ang = atan2(u(2,1), u(1,1)) * (180/pi);
    s = sqrt(w)*3;
    img = cv.ellipse(img, mu, s, 'Angle',ang, 'Color',clr, 'LineType','AA');
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function EmergencyTrajectory = calcEmergencyVelocityProfile_veldep(TargetTrajectory, ...
    drag_coefficient, roh_air, vehiclemass_kg, ...
    P_VDC_AccLimits_v_mps, P_VDC_AccLimits_ax_mps2, P_VDC_AccLimits_ay_mps2)

% Authors:       Alexander Wischnewski
%
% Description:
%   calculates a brake velocity profile to a given target trajectory to generate a reasonable
%   emergency backup

% Inputs:
%   TargetTrajectory    Trajectory to be checked (see tests for format)
%   drag_coefficient    Vehicle drag coefficient
%   roh_air             air density
%   vehiclemass_kg      vehicle mass
%   P_VDC_AccLimits_v_mps:  Acceleration limit map - velocity interpolation points
%   P_VDC_AccLimits_ax_mps2: Acceleration limit map - long. acc. limits
%   P_VDC_AccLimits_ay_mps2: Acceleration limit map - lat. acc. limits
%
% Outputs:
%   EmergencyTrajectory True if trajectory complies with all specs

% copy target trajectory
EmergencyTrajectory = TargetTrajectory;
% only if the trajectory is a driving trajectory
if(EmergencyTrajectory.v_mps(1) > 1)
    % calculate new velocity profile starting with the fifth element to ensure that the trajectory
    % does not begin prior to switching on it as this might cause serious issues
    for i = 6:1:50
        % distance to next point
        dS = double(EmergencyTrajectory.s_loc_m(i) - EmergencyTrajectory.s_loc_m(i-1));
        % update acceleration limits
        EmergencyTrajectory.ax_lim_mps2(i-1) = interp1(P_VDC_AccLimits_v_mps, P_VDC_AccLimits_ax_mps2, EmergencyTrajectory.v_mps(i-1));
        EmergencyTrajectory.ay_lim_mps2(i-1) = interp1(P_VDC_AccLimits_v_mps, P_VDC_AccLimits_ay_mps2, EmergencyTrajectory.v_mps(i-1));
        % calculate allowed longitudinal deceleration
        ax_mps2 = -(1 - abs(double(EmergencyTrajectory.kappa_radpm(i-1))*double(EmergencyTrajectory.v_mps(i-1)).^2/EmergencyTrajectory.ay_lim_mps2(i-1)))*EmergencyTrajectory.ax_lim_mps2(i-1);
        % add acceleration from drag
        ax_mps2 = ax_mps2 - 0.5*drag_coefficient*roh_air*EmergencyTrajectory.v_mps(i-1)^2/vehiclemass_kg;
        % calculate time needed to travel there using pq formula
        p = 2*EmergencyTrajectory.v_mps(i-1)/ax_mps2;
        q = -2*dS/ax_mps2;
        % check if the solution is still valid
        if((p/2)^2 > q && EmergencyTrajectory.v_mps(i-1) > 0.2)
            v_next = (-p/2 - sqrt((p/2)^2 - q))*ax_mps2 + EmergencyTrajectory.v_mps(i-1);
        else
            v_next = 0;
        end
        % set calculated speed and start next step with this speed
        EmergencyTrajectory.v_mps(i) = single(v_next);
        EmergencyTrajectory.ax_mps2(i-1) = ax_mps2;
    end
    % set final acceleration to value of before
    % use a lot less as this might lead to issues with the safety checks otherwise
    % speed should be zero here anyway
    EmergencyTrajectory.ax_mps2(50) = 0.5*EmergencyTrajectory.ax_mps2(49);
else
    % set speeds and acceleration to zero if the initial planning speed is too slow for driving
    EmergencyTrajectory.v_mps = zeros(50, 1);
    EmergencyTrajectory.ax_mps2 = zeros(50, 1);
end

author: TUMFTM


text:
function tDel = getTDel(planC)

 indexS = planC{end};
 tDelC = arrayfun(@(x) x.scanInfo(1).DICOMHeaders.TriggerTime,planC{indexS.scan},'un',0);
 nonZeroTDelV = [tDelC{:}]>0;
 n = find(nonZeroTDelV,1,'first');
 tDel = tDelC{n};

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [sets , scores] = mergeAll(I,labels,numlabels,iterationCount)

    graph = getLabelGraph(labels, numlabels+1);

    graphDistances= getGraphDistance(graph,numlabels);

    labelIndices = cell(1,numlabels);
    for i = 1:numlabels
        [row,col] = find(labels == i);
        labelIndices{1,i} = [row,col];
    end

    %GET GRADIENT OF IMAGE FOR COLOR-TEXTURE DISTANCE
    sigma = 0.5;

    Wx = floor((5/2)*sigma);
    if Wx < 1
      Wx = 1;
    end
    x = -Wx:Wx;

    % Evaluate 1D Gaussian filter (and its derivative).
    g = exp(-(x.^2)/(2*sigma^2));
    gp = -(x/sigma).*exp(-(x.^2)/(2*sigma^2));

    gradient = cell(2,3);

    gradient{1,1} = convolve2(convolve2(I(:,:,1),-gp,'same'),g','same');
    gradient{2,1} = convolve2(convolve2(I(:,:,1),g,'same'),-gp','same');

    gradient{1,2}= convolve2(convolve2(I(:,:,2),-gp,'same'),g','same');
    gradient{2,2} = convolve2(convolve2(I(:,:,2),g,'same'),-gp','same');

    gradient{1,3} = convolve2(convolve2(I(:,:,3),-gp,'same'),g','same');
    gradient{2,3} = convolve2(convolve2(I(:,:,3),g,'same'),-gp','same');

    irfx = gradient{1,1};
    irfy = gradient{2,1};

    igfx = gradient{1,2};
    igfy = gradient{2,2};

    ibfx = gradient{1,3};
    ibfy = gradient{2,3};

    %%%%%%%%

    orientations = cell(3,8);
    i = 1;
    for angle = 0:45:315
        orientations{1,i} = cos(angle*(pi/180))*irfx+sin(angle*(pi/180))*irfy;
        orientations{2,i} = cos(angle*(pi/180))*igfx+sin(angle*(pi/180))*igfy;
        orientations{3,i} = cos(angle*(pi/180))*ibfx+sin(angle*(pi/180))*ibfy;
        i = i + 1;
    end

    edges = -25:5:25;
    oHists = cell(numlabels,3,8);

    for l = 1:numlabels
        for o = 1:8
            for color = 1:3
                oHists{l,color, o} = histcounts(orientations{1,o}(labels == l),edges);
                oHists{l,color,o} = oHists{l,color, o} / sum(oHists{l,color, o});
            end
        end
    end

    ohists = cell(1,numlabels);

    for l = 1:numlabels
        ohists{1,numlabels} = zeros(24,10);
        for i = 1:8
            for j = 1:3
                for k = 1:10
                    ohists{1,l}((i-1) * 3 + j ,k) = oHists{l,j,i}(1,k);
                end
            end
        end
    end

    %GET COLOR HISTOGRAMS

    colorHists = cell(1,numlabels);

for l = 1:numlabels

    [row1,col1] = find(labels == l);
    edges = 0:1/20:1;

    rc1 = cat(2, row1,col1);

    ir1 = zeros(size(rc1,1) ,1);
    ig1 = zeros(size(rc1,1) ,1);
    ib1 = zeros(size(rc1,1), 1);

    for i = 1:size(rc1,1)
        ir1(i,1) = I(rc1(i,1),rc1(i,2) ,1);
        ig1(i,1) = I(rc1(i,1),rc1(i,2) ,2);
        ib1(i,1) = I(rc1(i,1),rc1(i,2) ,3);
    end

    r1 = histcounts(ir1,edges);
    r1= r1 / sum(r1);
    g1 = histcounts(ig1,edges);
    g1= g1 / sum(g1);
    b1 = histcounts(ib1,edges);
    b1= b1 / sum(b1);

    colorHists{1,l} = cat(1, r1 , g1 , b1);

    sets = cell(1,numlabels);

    for i = 1:numlabels
        sets{1,i} = i;
    end

end

 edgeImg = edge(rgb2gray(I),'Prewitt');
%[edgeImg, ~] = imgradient(rgb2gray(I),'prewitt');

scores = cell(1,2);
scoreCount = 1;

for i = 1:iterationCount
    [sets, lastMerged,~] = mergePixels(I,edgeImg, labels ,numlabels , graphDistances ,colorHists, ohists , sets , labelIndices);
    %score = scoreSet(edgeImg, labels, numlabels,lastMerged , labelIndices);
    score = getSophisticatedEdgeScore(edgeImg, labels, labelIndices, lastMerged);
    fprintf("Score: %f \n" , score);
    if(score > 0)
        scores{1,1}(1,scoreCount) = score;
        scores{1,2}{1,scoreCount} = lastMerged;
        scoreCount = scoreCount + 1;
    end
end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm

path: github-repos/MATLAB/zhangqianqianQQ-MachineVisionAlgorithm/MachineVisionAlgorithm-683338f6c3b1aab9fa2b80026915fe936aebf0ee/识别算法/Object Recognition based on super pixel/mergeAll.m


text:
%computes the novelty measure used by laroche
%> called by ::ComputeNoveltyFunction
%>
%> @param X: spectrogram (dimension FFTLength X Observations)
%> @param f_s: sample rate of audio data (unused)
%>
%> @retval d_lar novelty measure
% ======================================================================
function [d_lar] = NoveltyLaroche (X, f_s)

    % difference spectrum
    afDeltaX    = diff([sqrt(X(:,1)), sqrt(X)],1,2);

    % half-wave rectification
    afDeltaX(afDeltaX<0) = 0;

    % flux
    d_lar       = sum(afDeltaX)/size(X,1);
end

author: alexanderlerch

repo: ACA-Code


text:
function [ qy, qty, b, rsd, ab, info ] = dqrsl ( a, lda, n, k, qraux, y, job )

%*****************************************************************************80
%
%% DQRSL computes transformations, projections, and least squares solutions.
%
%  Discussion:
%
%    DQRSL requires the output of DQRDC.
%
%    For K <= min(N,P), let AK be the matrix
%
%      AK = ( A(JPVT(1)), A(JPVT(2)), ..., A(JPVT(K)) )
%
%    formed from columns JPVT(1), ..., JPVT(K) of the original
%    N by P matrix A that was input to DQRDC.  If no pivoting was
%    done, AK consists of the first K columns of A in their
%    original order.  DQRDC produces a factored orthogonal matrix Q
%    and an upper triangular matrix R such that
%
%      AK = Q * (R)
%               (0)
%
%    This information is contained in coded form in the arrays
%    A and QRAUX.
%
%    The parameters QY, QTY, B, RSD, and AB are not referenced
%    if their computation is not requested and in this case
%    can be replaced by dummy variables in the calling program.
%    To save storage, the user may in some cases use the same
%    array for different parameters in the calling sequence.  A
%    frequently occuring example is when one wishes to compute
%    any of B, RSD, or AB and does not need Y or QTY.  In this
%    case one may identify Y, QTY, and one of B, RSD, or AB, while
%    providing separate arrays for anything else that is to be
%    computed.
%
%    Thus the calling sequence
%
%      call dqrsl ( a, lda, n, k, qraux, y, dum, y, b, y, dum, 110, info )
%
%    will result in the computation of B and RSD, with RSD
%    overwriting Y.  More generally, each item in the following
%    list contains groups of permissible identifications for
%    a single calling sequence.
%
%      1. (Y,QTY,B) (RSD) (AB) (QY)
%
%      2. (Y,QTY,RSD) (B) (AB) (QY)
%
%      3. (Y,QTY,AB) (B) (RSD) (QY)
%
%      4. (Y,QY) (QTY,B) (RSD) (AB)
%
%      5. (Y,QY) (QTY,RSD) (B) (AB)
%
%      6. (Y,QY) (QTY,AB) (B) (RSD)
%
%    In any group the value returned in the array allocated to
%    the group corresponds to the last member of the group.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real A(LDA,P), contains the output of DQRDC.
%
%    Input, integer LDA, the leading dimension of the array A.
%
%    Input, integer N, the number of rows of the matrix AK.  It must
%    have the same value as N in DQRDC.
%
%    Input, integer K, the number of columns of the matrix AK.  K
%    must not be greater than min(N,P), where P is the same as in the
%    calling sequence to DQRDC.
%
%    Input, real QRAUX(P), the auxiliary output from DQRDC.
%
%    Input, real Y(N), a vector to be manipulated by DQRSL.
%
%    Input, integer JOB, specifies what is to be computed.  JOB has
%    the decimal expansion ABCDE, with the following meaning:
%      if A /= 0, compute QY.
%      if B /= 0, compute QTY.
%      if C /= 0, compute QTY and B.
%      if D /= 0, compute QTY and RSD.
%      if E /= 0, compute QTY and AB.
%    Note that a request to compute B, RSD, or AB automatically triggers
%    the computation of QTY, for which an array must be provided in the
%    calling sequence.
%
%    Output, real QY(N), contains Q * Y, if requested.
%
%    Output, real QTY(N), contains Q' * Y, if requested.
%
%    Output, real B(K), the solution of the least squares problem
%      minimize norm2 ( Y - AK * B),
%    if its computation has been requested.  Note that if pivoting was
%    requested in DQRDC, the J-th component of B will be associated with
%    column JPVT(J) of the original matrix A that was input into DQRDC.
%
%    Output, real RSD(N), the least squares residual Y - AK * B,
%    if its computation has been requested.  RSD is also the orthogonal
%    projection of Y onto the orthogonal complement of the column space
%    of AK.
%
%    Output, real AB(N), the least squares approximation Ak * B,
%    if its computation has been requested.  AB is also the orthogonal
%    projection of Y onto the column space of A.
%
%    Output, integer INFO, is zero unless the computation of B has
%    been requested and R is exactly singular.  In this case, INFO is the
%    index of the first zero diagonal element of R, and B is left unaltered.
%

%
%  Set info flag.
%
  info = 0;
%
%  Determine what is to be computed.
%
  cqy =  floor       ( job / 10000         ) ~= 0;
  cqty =         mod ( job,  10000 )         ~= 0;
  cb =   floor ( mod ( job,   1000 ) / 100 ) ~= 0;
  cr =   floor ( mod ( job,    100 ) /  10 ) ~= 0;
  cab =          mod ( job,     10 )       ~= 0;

  ju = min ( k, n - 1 );
%
%  Special action when N = 1.
%
  if ( ju == 0 )

    qy(1) = y(1);
    qty(1) = y(1);
    ab(1) = y(1);

    if ( a(1,1) == 0.0 )
      info = 1;
    else
      b(1) = y(1) / a(1,1);
    end

    rsd(1) = 0.0;

    return

  end
%
%  Set up to compute QY or QTY.
%
  qy(1:n) = y(1:n);
  qty(1:n) = y(1:n);
%
%  Compute QY.
%
  if ( cqy )

    for jj = 1 : ju

      j = ju - jj + 1;

      if ( qraux(j) ~= 0.0 )
        temp = a(j,j);
        a(j,j) = qraux(j);
        t = -ddot ( n-j+1, a(j:n,j), 1, qy(j:n), 1 ) / a(j,j);
        qy(j:n) = daxpy ( n-j+1, t, a(j:n,j)', 1, qy(j:n), 1 );
        a(j,j) = temp;
      end

    end

  end
%
%  Compute Q'*Y.
%
  if ( cqty )

    for j = 1 : ju
      if ( qraux(j) ~= 0.0 )
        temp = a(j,j);
        a(j,j) = qraux(j);
        t = -ddot ( n-j+1, a(j:n,j), 1, qty(j:n), 1 ) / a(j,j);
        qty(j:n) = daxpy ( n-j+1, t, a(j:n,j), 1, qty(j:n), 1 );
        a(j,j) = temp;
      end
    end

  end
%
%  Set up to compute B, RSD, or AB.
%
  b(1:k) = qty(1:k);
  ab(1:k) = qty(1:k);
  rsd(k+1:n) = qty(k+1:n);
  ab(k+1:n) = 0.0;
  rsd(1:k) = 0.0;
%
%  Compute B.
%
  if ( cb )

    for jj = 1 : k

      j = k - jj + 1;

      if ( a(j,j) == 0.0 )
        info = j;
        break
      end

      b(j) = b(j) / a(j,j);

      if ( j ~= 1 )
        t = -b(j);
        b(1:j-1) = daxpy ( j-1, t, a(1:j-1,j), 1, b(1:j-1), 1 );
      end

    end

  end
%
%  Compute RSD or AB as required.
%
  if ( cr | cab )

    for jj = 1 : ju

      j = ju - jj + 1;

      if ( qraux(j) ~= 0.0 )

        temp = a(j,j);
        a(j,j) = qraux(j);

        if ( cr )
          rsd(j:n) = -ddot ( n-j+1, a(j:n,j), 1, rsd(j:n), 1 ) / a(j,j);
          rsd(j:n) = daxpy ( n-j+1, t, a(j:n,j), 1, rsd(j:n), 1 );
        end

        if ( cab )
          t = -ddot ( n-j+1, a(j:n,j), 1, ab(j:n), 1 ) / a(j,j);
          ab(j:n) = daxpy ( n-j+1, t, a(j:n,j), 1, ab(j:n), 1 );
        end

        a(j,j) = temp;

      end
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m