text:
%PGraph Graph class
%
%   g = PGraph()        create a 2D, planar embedded, directed graph
%   g = PGraph(n)       create an n-d, embedded, directed graph
%
% Provides support for graphs that:
%   - are directed
%   - are embedded in a coordinate system (2D or 3D)
%   - have multiple unconnected components
%   - have symmetric cost edges (A to B is same cost as B to A)
%   - have no loops (edges from A to A)
%
% Graph representation:
%   - vertices are represented by integer vertex ids (vid)
%   - edges are represented by integer edge ids (eid)
%   - each vertex can have arbitrary associated data
%   - each edge can have arbitrary associated data
%
% Methods::
%
% Constructing the graph::
%   g.add_node(coord)      add vertex
%   g.add_edge(v1, v2)     add edge fbetween vertices
%   g.setcost(e, c)        set cost for edge
%   g.setedata(e, u)       set user data for edge
%   g.setvdata(v, u)       set user data for vertex
%
% Modifying the graph::
%   g.clear()              remove all vertices and edges from the graph
%   g.delete_edge(e)       remove edge
%   g.delete_node(v)       remove vertex
%   g.setcoord(v)          set coordinate of vertex
%
% Information from graph::
%   g.about()                   summary information about node
%   g.component(v)              component id for vertex
%   g.componentnodes(c)         vertices in component
%   g.connectivity()            number of edges for all vertices
%   g.connectivity_in()         number of incoming edges for all vertices
%   g.connectivity_out()        number of outgoing edges for all vertices
%   g.coord(v)                  coordinate of vertex
%   g.cost(e)                   cost of edge
%   g.distance_metric(v1,v2)    distance between nodes
%   g.edata(e)                  get edge user data
%   g.edgedir(v1,v2)            direction of edge
%   g.edges(v)                  list of edges for vertex
%   g.edges_in(v)               list of edges into vertex
%   g.edges_out(v)              list of edges from vertex
%   g.lookup(name)              vertex from name
%   g.name(v)                   name of vertex
%   g.neighbours(v)             neighbours of vertex
%   g.neighbours_d(v)           neighbours of vertex and edge directions
%   g.neighbours_in(v)          neighbours with edges in
%   g.neighbours_out(v)         neighbours with edges out
%   g.samecomponent(v1,v2)      test if vertices in same component
%   g.vdata(v)                  vertex user data
%   g.vertices(e)               vertices for edge
%
% Display::
%
%   g.char()                   convert graph to string
%   g.display()                display summary of graph
%   g.highlight_node(v)        highlight vertex
%   g.highlight_edge(e)        highlight edge
%   g.highlight_component(c)   highlight all nodes in component
%   g.highlight_path(p)        highlight nodes and edge along path
%   g.pick(coord)              vertex closest to coord
%   g.plot()                   plot graph
%
%
% Matrix representations::
%   g.adjacency()          adjacency matrix
%   g.degree()             degree matrix
%   g.incidence()          incidence matrix
%   g.laplacian()          Laplacian  matrix
%
% Planning paths through the graph::
%   g.Astar(s, g)          shortest path from s to g
%   g.goal(v)              set goal vertex, and plan paths
%   g.path(v)              list of vertices from v to goal
%
% Graph and world points::
%   g.closest(coord)       vertex closest to coord
%   g.coord(v)             coordinate of vertex v
%   g.distance(v1, v2)     distance between v1 and v2
%   g.distances(coord)     return sorted distances from coord to all vertices
%
% Object properties (read only)::
%   g.n            number of vertices
%   g.ne           number of edges
%   g.nc           number of components
%
% Example::
%         g = PGraph();
%         g.add_node([1 2]');  % add node 1
%         g.add_node([3 4]');  % add node 1
%         g.add_node([1 3]');  % add node 1
%         g.add_edge(1, 2);    % add edge 1-2
%         g.add_edge(2, 3);    % add edge 2-3
%         g.add_edge(1, 3);    % add edge 1-3
%         g.plot()
%
% Notes::
% - Support for edge direction is quite simple.
% - The method distance_metric() could be redefined in a subclass.

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math

% Peter Corke 8/2009.

% TODO:
%       be able to delete nodes, must update connectivity

% update to use map.container class

classdef PGraph < matlab.mixin.Copyable

    properties (SetAccess=private, GetAccess=public)
        vertexlist      % vertex coordinates, columnwise, vertex number is the column number
        edgelist        % 2xNe matrix, each column is vertex index of edge start and end
        edgelen         % length (cost) of this edge

        labels          % label of each vertex (1xN)
        maxlabel        % set of all labels (1xNc)

        goaldist        % distance from goal, after planning

        vertexdata      % per vertex data, cell array
        edgedata        % per edge data, cell array
        ndims           % number of coordinate dimensions, height of vertices matrix
        verbose
        measure         % distance measure: 'Euclidean', 'SE2'
        dweight         % distance weighting for SE2 measure
        ncvalid
        names
    end

    properties (Dependent)
        n               % number of nodes/vertices
        ne              % number of edges
        nc              % number of components
    end

    methods

        function g = PGraph(ndims, varargin)
            %PGraph.PGraph Graph class constructor
            %
            % G=PGraph(D, OPTIONS) is a graph object embedded in D dimensions.
            %
            % Options::
            %  'distance',M   Use the distance metric M for path planning which is either
            %                 'Euclidean' (default) or 'SE2'.
            %  'verbose'      Specify verbose operation
            %
            % Notes::
            % - Number of dimensions is not limited to 2 or 3.
            % - The distance metric 'SE2' is the sum of the squares of the difference
            %   in position and angle modulo 2pi.
            % - To use a different distance metric create a subclass of PGraph and
            %   override the method distance_metric().

            if nargin < 1
                ndims = 2;  % planar by default
            elseif isa(ndims, 'PGraph')
                % do a deep copy
                g = ndims.copy();
                return
            end
            g.ndims = ndims;
            opt.distance = 'Euclidean';
            opt.dweight = 1;
            opt = tb_optparse(opt, varargin);

            g.clear();
            g.verbose = opt.verbose;
            g.measure = opt.distance;
            g.dweight = opt.dweight;
            g.vertexdata = {};
            g.edgedata = {};
            g.ncvalid = false;  % mark connectivity as suspect
            g.names = strings(0,0);
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% GRAPH MAINTENANCE
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        function v = add_node(g, coord, varargin)
            %PGraph.add_node Add a node
            %
            % V = G.add_node(X, OPTIONS) adds a node/vertex with coordinate X (Dx1) and
            % returns the integer node id V.
            %
            % Options:
            % 'name',N    Assign a string name N to this vertex
            % 'from',V    Create a directed edge from vertex V with cost equal to the distance between the vertices.
            % 'cost',C    If an edge is created use cost C
            %
            % Notes::
            % - Distance is computed according to the metric specified in the
            %   constructor.
            %
            % See also PGraph.add_edge, PGraph.data, PGraph.getdata.

            if length(coord) ~= g.ndims
                error('coordinate length different to graph coordinate dimensions');
            end

            opt.from = [];
            opt.name = [];
            opt.cost = NaN;

            opt = tb_optparse(opt, varargin);

            % append the coordinate as a column in the vertex matrix
            g.vertexlist = [g.vertexlist coord(:)];
            v = g.n;

            if g.verbose
                fprintf('add node (%d) = ', v);
                fprintf('%f ', coord);
                fprintf('\n');
            end

            % optionally add an edge
            if ~isempty(opt.from)
                if isnan(opt.cost)
                    opt.cost = g.distance(v, opt.from);
                end
                g.add_edge(opt.from, v, opt.cost);
            end

            if ~isempty(opt.name)
                g.names(v) = opt.name;
            end
            g.ncvalid = false;  % mark connectivity as suspect
        end

        function e = add_edge(g, v1, v2, d)
            %PGraph.add_edge Add an edge
            %
            % E = G.add_edge(V1, V2) adds a directed edge from vertex id V1 to vertex id V2, and
            % returns the edge id E.  The edge cost is the distance between the vertices.
            %
            % E = G.add_edge(V1, V2, C) as above but the edge cost is C.
            %
            % Notes::
            % - If V2 is a vector add edges from V1 to all elements of V2
            % - Distance is computed according to the metric specified in the
            %   constructor.
            %
            % See also PGraph.add_node, PGraph.edgedir.
            if g.verbose
                fprintf('add edge %d -> %d\n', v1, v2);
            end
            e = [];
            for vv=v2(:)'
                g.edgelist = [g.edgelist [v1; vv]];
                e = [e numcols(g.edgelist)];
                if (nargin < 4) || isempty(d)
                    d = g.distance(v1, vv);
                end
                g.edgelen = [g.edgelen d];
            end
            g.ncvalid = false;  % mark connectivity as suspect

        end

        function delete_node(g, vv)

            for v=sort(vv(:)', 2, 'descend')
                % remove all its edges
                el = g.edges(v);
                g.delete_edge(el);

                % remove the column from the vertex table
                g.vertexlist(:,v) = [];

                % now renumber all the edges that might have changed
                k = find(g.edgelist > v);
                g.edgelist(k) = g.edgelist(k) - 1;

                g.ncvalid = false;  % mark connectivity as suspect
            end
        end

        function delete_edge(g, e)
            g.edgelist(:,e) = [NaN; NaN];
            g.ncvalid = false;  % mark connectivity as suspect

        end

        function clear(g)
            %PGraph.clear Clear the graph
            %
            % G.clear() removes all vertices, edges and components.

            % set the orientation of the edge and vertex tables
            g.labels = zeros(1, 0);
            g.edgelist = zeros(2, 0);
            g.edgelen = zeros(1, 0);
            g.vertexlist = zeros(g.ndims, 0);
            g.ncvalid = false;  % mark connectivity as suspect
            g.maxlabel = 0;
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% GRAPH STRUCTURE
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % which edges contain v
        %  elist = g.edges(v)
        function e = edges(g, v)
            %PGraph.edges Find edges given vertex
            %
            % E = G.edges(V) is a vector containing the id of all edges connected to vertex id V.
            %
            % See also PGraph.edgedir.

            e = [find(g.edgelist(1,:) == v) find(g.edgelist(2,:) == v)];
        end

        function e = edges_in(g, v)
            %PGraph.edges Find edges given vertex
            %
            % E = G.edges(V) is a vector containing the id of all edges connected to vertex id V.
            %
            % See also PGraph.edgedir.
            e = find(g.edgelist(2,:) == v);
        end

        function e = edges_out(g, v)
            %PGraph.edges Find edges given vertex
            %
            % E = G.edges(V) is a vector containing the id of all edges connected to vertex id V.
            %
            % See also PGraph.edgedir.
            e = find(g.edgelist(1,:) == v);
        end

        function dir = edgedir(g, v1, v2)
            %PGraph.edgedir Find edge direction
            %
            % D = G.edgedir(V1, V2) is the direction of the edge from vertex id V1
            % to vertex id V2.
            %
            % If we add an edge from vertex 3 to vertex 4
            %         g.add_edge(3, 4)
            % then
            %         g.edgedir(3, 4)
            % is positive, and
            %         g.edgedir(4, 3)
            % is negative.
            %
            % See also PGraph.add_node, PGraph.add_edge.
            n = g.edges(v1);
            if any(ismember( g.edgelist(2, n), v2))
                dir = 1;
            elseif any(ismember( g.edgelist(1, n), v2))
                dir = -1;
            else
                dir = 0;
            end
        end

        function v = vertices(g, e)
            %PGraph.vertices Find vertices given edge
            %
            % V = G.vertices(E) return the id of the vertices that define edge E.
            v = g.edgelist(:,e);
        end

        function [n,c] = neighbours(g, v)
            %PGraph.neighbours Neighbours of a vertex
            %
            % N = G.neighbours(V) is a vector of ids for all vertices which are
            % directly connected neighbours of vertex V.
            %
            % [N,C] = G.neighbours(V) as above but also returns a vector C whose elements
            % are the edge costs of the paths corresponding to the vertex ids in N.
            e = g.edges(v);
            n = g.edgelist(:,e);
            n = n(:)';
            n(n==v) = [];   % remove references to self
            if nargout > 1
                c = g.cost(e);
            end
        end

        function [n,c] = neighbours_out(g, v)
            %PGraph.neighbours Outgoing neighbours of a vertex
            %
            % N = G.neighbours(V) is a vector of ids for all vertices which are
            % directly connected neighbours of vertex V.
            %
            % [N,C] = G.neighbours(V) as above but also returns a vector C whose elements
            % are the edge costs of the paths corresponding to the vertex ids in N.
            e = g.edges_out(v);
            n = g.edgelist(:,e);
            n = n(:)';
            n(n==v) = [];   % remove references to self
            if nargout > 1
                c = g.cost(e);
            end
        end

        function [n,c] = neighbours_in(g, v)
            %PGraph.neighbours Incoming neighbours of a vertex
            %
            % N = G.neighbours(V) is a vector of ids for all vertices which are
            % directly connected neighbours of vertex V.
            %
            % [N,C] = G.neighbours(V) as above but also returns a vector C whose elements
            % are the edge costs of the paths corresponding to the vertex ids in N.
            e = g.edges_in(v);
            n = g.edgelist(:,e);
            n = n(:)';
            n(n==v) = [];   % remove references to self
            if nargout > 1
                c = g.cost(e);
            end
        end

        function [n,c] = neighbours_d(g, v)
            %PGraph.neighbours_d Directed neighbours of a vertex
            %
            % N = G.neighbours_d(V) is a vector of ids for all vertices which are
            % directly connected neighbours of vertex V.  Elements are positive
            % if there is a link from V to the node (outgoing), and negative if the link
            % is from the node to V (incoming).
            %
            % [N,C] = G.neighbours_d(V) as above but also returns a vector C whose elements
            % are the edge costs of the paths corresponding to the vertex ids in N.
            e = g.edges(v);
            n = [-g.edgelist(1,e) g.edgelist(2,e)];
            n(abs(n)==v) = [];   % remove references to self
            if nargout > 1
                c = g.cost(e);
            end
        end

        function c = connectivity(g,nn)
            %PGraph.connectivity Node connectivity
            %
            % C = G.connectivity() is a vector (Nx1) with the number of edges per
            % vertex.
            %
            % The average vertex connectivity is
            %         mean(g.connectivity())
            %
            % and the minimum vertex connectivity is
            %         min(g.connectivity())

            if nargin == 1
                for k=1:g.n
                    c(k) = length(g.edges(k));
                end
            elseif nargin == 2
                c = [];
                for n=nn
                    c = [c length(g.edges(n))];
                end
            end
        end

        function c = connectivity_in(g,n )
            %PGraph.connectivity Graph connectivity
            %
            % C = G.connectivity() is a vector (Nx1) with the number of incoming edges per
            % vertex.
            %
            % The average vertex connectivity is
            %         mean(g.connectivity())
            %
            % and the minimum vertex connectivity is
            %         min(g.connectivity())

            if nargin == 1
                for k=1:g.n
                    c(k) = length(g.edges_in(k));
                end
            elseif nargin == 2
                c = length(g.edges_in(n));
            end
        end

        function c = connectivity_out(g,n )
            %PGraph.connectivity Graph connectivity
            %
            % C = G.connectivity() is a vector (Nx1) with the number of outgoing edges per
            % vertex.
            %
            % The average vertex connectivity is
            %         mean(g.connectivity())
            %
            % and the minimum vertex connectivity is
            %         min(g.connectivity())

            if nargin == 1
                for k=1:g.n
                    c(k) = length(g.edges_out(k));
                end
            elseif nargin == 2
                c = length(g.edges_out(n));
            end
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% NODE PROPERTIES
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        function p = coord(g, v)
            %PGraph.coord Coordinate of node
            %
            % X = G.coord(V) is the coordinate vector (Dx1) of vertex id V.

            if nargin < 2
                p = g.vertexlist;
            else
                p = g.vertexlist(:,v);
            end
        end

        function p = name(g, v)
            %PGraph.coord Name of node
            %
            % X = G.name(V) is the name (string) of vertex id V.

            if nargin < 2
                p = [g.names];
            else
                p = g.names(v);
            end
        end

        function p = lookup(g, name)
            p  = find( [g.names] == name );
        end

        function p = setcoord(g, v, c)
            %PGraph.coord Coordinate of node
            %
            % X = G.coord(V) is the coordinate vector (Dx1) of vertex id V.

            if nargin < 3
                if ~all(size(v) == size(g.vertexlist))
                    error('SMTB:PGraph:badarg', 'value must be size of vertex table');
                end
                g.vertexlist = v;
            else
                g.vertexlist(:,v) = c;
            end
                end

        function u = vdata(g, v)
            %PGraph.data Get user data for node
            %
            % U = G.data(V) gets the user data of vertex V which can be of any
            % type such as a number, struct, object or cell array.
            %
            % See also PGraph.setdata.
            u = g.vertexdata{v};
        end

        function u = setvdata(g, v, u)
            %PGraph.setdata Set user data for node
            %
            % G.setdata(V, U) sets the user data of vertex V to U which can be of any
            % type such as a number, struct, object or cell array.
            %
            % See also PGraph.data.

            g.vertexdata{v} = u;
        end

        function d = distance(g, v1, v2)
            %PGraph.distance Distance between vertices
            %
            % D = G.distance(V1, V2) is the geometric distance between
            % the vertices V1 and V2.
            %
            % See also PGraph.distances.

            d = g.distance_metric( g.vertexlist(:,v1), g.vertexlist(:,v2));

        end

        function [d,k] = distances(g, p)
            %PGraph.distances Distances from point to vertices
            %
            % D = G.distances(X) is a vector (1xN) of geometric distance from the point
            % X (Dx1) to every other vertex sorted into increasing order.
            %
            % [D,W] = G.distances(P) as above but also returns W (1xN) with the
            % corresponding vertex id.
            %
            % Notes::
            % - Distance is computed according to the metric specified in the
            %   constructor.
            %
            % See also PGraph.closest.

            d = g.distance_metric(p(:), g.vertexlist);
            [d,k] = sort(d, 'ascend');
        end

        function [c,dn] = closest(g, p, tol)
            %PGraph.closest Find closest vertex
            %
            % V = G.closest(X) is the vertex geometrically closest to coordinate X.
            %
            % [V,D] = G.closest(X) as above but also returns the distance D.
            %
            % See also PGraph.distances.
            d = g.distance_metric(p(:), g.vertexlist);
            [mn,c] = min(d);

            if nargin > 2 && mn > tol
                c = []; dn = [];
            end

            if nargout > 1
                dn = mn;
            end

        end

        function about(g, vv)
            if nargin < 2
                disp('pick a node using the mouse');
                vv = g.pick()
            end

            if ~g.ncvalid
                g.graphcolor();
            end

            for v=vv
                fprintf('Node %d #%d@ (', v, g.labels(v)); fprintf('%g ', g.coord(v)); fprintf(')\n');
                fprintf('  neighbours: ');
                fprintf('%d ', g.neighbours_in(v)); fprintf(' >-o-> ');
                fprintf('%d ', g.neighbours_out(v)); fprintf('\n');

                fprintf('  edges: ');
                fprintf('%d ', g.edges_in(v)); fprintf(' >-o-> ');
                fprintf('%d ', g.edges_out(v)); fprintf('\n');
            end
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% EDGE PROPERTIES
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        function d = cost(g, e)
            %PGraph.cost Cost of edge
            %
            % C = G.cost(E) is the cost of edge id E.
            d = g.edgelen(e);
        end

        function d = setcost(g, e, c)
            %PGraph.cost Set cost of edge
            %
            % G.setcost(E, C) set cost of edge id E to C.
            g.edgelen(e) = c;
        end

        function u = edata(g, e)
            %PGraph.data Get user data for node
            %
            % U = G.data(V) gets the user data of vertex V which can be of any
            % type such as a number, struct, object or cell array.
            %
            % See also PGraph.setdata.
            u = g.edgedata{e};
        end

        function u = setedata(g, e, u)
            %PGraph.setdata Set user data for node
            %
            % G.setdata(V, U) sets the user data of vertex V to U which can be of any
            % type such as a number, struct, object or cell array.
            %
            % See also PGraph.data.

            g.edgedata{e} = u;
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% GRAPH INFORMATION
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        function n = get.n(g)
            %Pgraph.n Number of vertices
            %
            % G.n is the number of vertices in the graph.
            %
            % See also PGraph.ne.
            n = numcols(g.vertexlist);
        end

        function ne = get.ne(g)
            %Pgraph.ne Number of edges
            %
            % G.ne is the number of edges in the graph.
            %
            % See also PGraph.n.
            ne = numcols(g.edgelist);
        end

        function ne = get.nc(g)
            %Pgraph.nc Number of components
            %
            % G.nc is the number of components in the graph.
            %
            % See also PGraph.component.

            if ~g.ncvalid
                g.graphcolor();
            end
            ne = g.maxlabel;
        end

        function display(g)
            %PGraph.display Display graph
            %
            % G.display() displays a compact human readable representation of the
            % state of the graph including the number of vertices, edges and components.
            %
            % See also PGraph.char.
            loose = strcmp( get(0, 'FormatSpacing'), 'loose');
            if loose
                disp(' ');
            end
            disp([inputname(1), ' = '])
            disp( char(g) );
        end % display()

        function s = char(g)
            %PGraph.char Convert graph to string
            %
            % S = G.char() is a compact human readable representation of the
            % state of the graph including the number of vertices, edges and components.

            s = '';
            s = char(s, sprintf('  %d dimensions', g.ndims));
            s = char(s, sprintf('  %d vertices', g.n));
            s = char(s, sprintf('  %d edges', numcols(g.edgelist)));
            s = char(s, sprintf('  %d components', g.nc));
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% GRAPH COMPONENTS
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        function graphcolor(g)
            % color the graph

            g.labels = repmat(NaN, 1, g.n);

            function colorComponent(g, v, l)
                g.labels(v) = l;
                for n = g.neighbours(v)
                    if isnan(g.labels(n))
                        colorComponent(g, n, l);
                    end
                end
            end

            for label = 1:g.n
                % find first vertex with no label
                v = find(isnan(g.labels));
                if isempty(v)
                    g.maxlabel = label-1;
                    break;
                end
                v = v(1);

                colorComponent(g, v, label);
            end
            g.ncvalid = true;

        end

        function c = component(g, v)
            %PGraph.component Graph component
            %
            % C = G.component(V) is the id of the graph component that contains vertex
            % V.
            c = g.labels(v);
        end

        function v = componentnodes(g, c)
            %PGraph.component Graph component
            %
            % C = G.component(V) are the ids of all vertices in the graph component V.
            v = find(g.labels == c);
        end

        function c = samecomponent(g, v1, v2)
            %PGraph.component Graph component
            %
            % C = G.component(V) is the id of the graph component that contains vertex
            % V.

            if ~g.ncvalid
                % lazy graph coloring
                g.graphcolor();
            end
            c = g.labels(v1) == g.labels(v2);
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% GRAPHICS
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        function plot(g, varargin)
            %PGraph.plot Plot the graph
            %
            % G.plot(OPT) plots the graph in the current figure.  Nodes
            % are shown as colored circles.
            %
            % Options::
            %  'labels'              Display vertex id (default false)
            %  'edges'               Display edges (default true)
            %  'edgelabels'          Display edge id (default false)
            %  'NodeSize',S          Size of vertex circle (default 8)
            %  'NodeFaceColor',C     Node circle color (default blue)
            %  'NodeEdgeColor',C     Node circle edge color (default blue)
            %  'NodeLabelSize',S     Node label text sizer (default 16)
            %  'NodeLabelColor',C    Node label text color (default blue)
            %  'EdgeColor',C         Edge color (default black)
            %  'EdgeLabelSize',S     Edge label text size (default black)
            %  'EdgeLabelColor',C    Edge label text color (default black)
            %  'componentcolor'      Node color is a function of graph component
            %  'only',N              Only show these nodes

            % show vertices
            holdon = ishold;
            hold on

            % parse options
            opt.componentcolor = false;
            opt.labels = false;
            opt.edges = true;
            opt.edgelabels = false;
            opt.NodeSize = 8;
            opt.NodeFaceColor = 'b';
            opt.NodeEdgeColor = 'b';
            opt.NodeLabelSize = 16;
            opt.NodeLabelColor = 'b';
            opt.EdgeColor = 'k';
            opt.EdgeLabelSize = 8;
            opt.EdgeLabelColor = 'k';
            opt.EdgeWidth = 0.5;
            opt.dims = g.ndims;
            opt.only = [1:g.n];

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

            % set default color if none specified
            if ~isempty(args)
                mcolor = args{1};
            else
                mcolor = 'b';
            end

            if opt.componentcolor

                colororder = get(gca, 'ColorOrder');

                % step through each component
                for c=1:g.nc
                    vertices = g.componentnodes(c);
                    if length(vertices) == 1
                        % singleton is grey
                        color = 0.8*[1 1 1];
                    else
                        % otherwise use next color from the axis color order
                        color = colororder(mod(c+1,numrows(colororder)-1)+1,:);
                    end

                                        % plot the edges
                    if opt.edges
                        coords = [];

                        for v = vertices
                            p0 = g.coord(v); p0 = p0(1:opt.dims);
                            for n = g.neighbours(v)
                                pn = g.coord(n); pn = pn(1:opt.dims);
                                coords = [coords; p0'; pn'; NaN*ones(1,g.ndims)];
                            end
                        end
                        if ~isempty(coords)
                            if opt.dims == 3
                                plot3(coords(1,:), coords(2,:), coords(3,:), 'Color', color, 'LineWidth', opt.EdgeWidth);
                            else
                                plot(coords(:,1), coords(:,2), 'Color', color, 'LineWidth', opt.EdgeWidth);
                            end
                        end
                    end

                    % plot the nodes
                    args = {'LineStyle', 'None', ...
                        'Marker', 'o', ...
                        'MarkerFaceColor', color, ...
                        'MarkerSize', opt.NodeSize, ...
                        'MarkerEdgeColor', color };

                    for v = vertices
                        if ~ismember(v, opt.only)
                            continue;
                        end
                        p = g.coord(v);
                        if opt.dims == 3
                            plot3(p(1), p(2), p(3), args{:});
                        else
                            plot(p(1), p(2), args{:});
                        end
                    end

                    if length(vertices) == 1
                        continue; % no edges to plot
                    end

                end
            else

                %% user selected colors

                % show edges
                if opt.edges
                for e=g.edgelist
                    v1 = g.vertexlist(:,e(1));
                    v2 = g.vertexlist(:,e(2));
                    if opt.dims == 3
                        plot3([v1(1) v2(1)], [v1(2) v2(2)], [v1(3) v2(3)], ...
                            'Color', opt.EdgeColor, 'LineWidth', opt.EdgeWidth);
                    else
                        plot([v1(1) v2(1)], [v1(2) v2(2)], ...
                            'Color', opt.EdgeColor, 'LineWidth', opt.EdgeWidth);
                    end
                end
                end

                % show the vertices as filled circles
                for i=1:g.n
                    % for each node
                        if ~ismember(i, opt.only)
                            continue;
                        end
                    args = {'LineStyle', 'None', ...
                        'Marker', 'o', ...
                        'MarkerFaceColor', opt.NodeFaceColor, ...
                        'MarkerSize', opt.NodeSize, ...
                        'MarkerEdgeColor', opt.NodeEdgeColor };
                    if opt.dims == 3
                        plot3(g.vertexlist(1,i), g.vertexlist(2,i), g.vertexlist(3,i), args{:});
                    else
                        plot(g.vertexlist(1,i), g.vertexlist(2,i), args{:});
                    end
                end
            end

            % show the edge labels
            if opt.edgelabels
                for i=1:g.ne
                    e = g.edgelist(:,i);
                    v1 = g.vertexlist(:,e(1));
                    v2 = g.vertexlist(:,e(2));

                    text('String', sprintf('  %g', g.cost(i)), ...
                        'Position', (v1 + v2)/2, ...
                        'HorizontalAlignment', 'left', ...
                        'VerticalAlignment', 'middle', ...
                        'FontUnits', 'pixels', ...
                        'FontSize', opt.EdgeLabelSize, ...
                        'Color', opt.EdgeLabelColor);
                end
            end
            % show the labels
            if opt.labels
                for i=1:g.n
                    text('String', sprintf('  %d', i), ...
                        'Position', g.vertexlist(:,i), ...
                        'HorizontalAlignment', 'left', ...
                        'VerticalAlignment', 'middle', ...
                        'FontUnits', 'pixels', ...
                        'FontSize', opt.NodeLabelSize, ...
                        'Color', opt.NodeLabelColor);
                end
            end
            if ~holdon
                hold off
            end
            grid on
        end

        function v = pick(g)
            %PGraph.pick Graphically select a vertex
            %
            % V = G.pick() is the id of the vertex closest to the point clicked
            % by the user on a plot of the graph.
            %
            % See also PGraph.plot.
            [x,y] = ginput(1);
            d = colnorm( bsxfun(@minus,[x; y], g.vertexlist(1:2,:)) );

            [~,v] = min(d);
        end

        function highlight_node(g, verts, varargin)
            %PGraph.highlight_node Highlight a node
            %
            % G.highlight_node(V, OPTIONS) highlights the vertex V with a yellow marker.
            % If V is a list of vertices then all are highlighted.
            %
            % Options::
            %  'NodeSize',S          Size of vertex circle (default 12)
            %  'NodeFaceColor',C     Node circle color (default yellow)
            %  'NodeEdgeColor',C     Node circle edge color (default blue)
            %
            % See also PGraph.highlight_edge, PGraph.highlight_path, PGraph.highlight_component.

            hold on

            % parse options
            opt.NodeSize = 12;
            opt.NodeFaceColor = 'y';
            opt.NodeEdgeColor = 'b';

            [opt,args] = tb_optparse(opt, varargin);
            markerprops = {'LineStyle', 'None', ...
                'Marker', 'o', ...
                'MarkerFaceColor', opt.NodeFaceColor, ...
                'MarkerSize', opt.NodeSize, ...
                'MarkerEdgeColor', opt.NodeEdgeColor };

            for v=verts
                if g.ndims == 3
                    plot3(g.vertexlist(1,v), g.vertexlist(2,v), g.vertexlist(3,v), ...
                        markerprops{:});
                else
                    plot(g.vertexlist(1,v), g.vertexlist(2,v), markerprops{:});
                end
            end
        end

        function highlight_component(g, c, varargin)
            %PGraph.highlight_component Highlight a graph component
            %
            % G.highlight_component(C, OPTIONS) highlights the vertices that belong to
            % graph component C.
            %
            % Options::
            %  'NodeSize',S          Size of vertex circle (default 12)
            %  'NodeFaceColor',C     Node circle color (default yellow)
            %  'NodeEdgeColor',C     Node circle edge color (default blue)
            %
            % See also PGraph.highlight_node, PGraph.highlight_edge, PGraph.highlight_component.
            nodes = find(g.labels == g.labelset(c));
            for v=nodes
                g.highlight_node(v, varargin{:});
            end
        end

        function highlight_edge(g, e, varargin)
            %PGraph.highlight_node Highlight a node
            %
            % G.highlight_edge(V1, V2) highlights the edge between vertices V1 and V2.
            %
            % G.highlight_edge(E) highlights the edge with id E.
            %
            % Options::
            % 'EdgeColor',C         Edge edge color (default black)
            % 'EdgeThickness',T     Edge thickness (default 1.5)
            %
            % See also PGraph.highlight_node, PGraph.highlight_path, PGraph.highlight_component.

            % parse options
            opt.EdgeColor = 'k';
            opt.EdgeThickness = 1.5;

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

            hold on
            if (length(args) > 0) && isnumeric(args{1})
                % highlight_edge(V1, V2)
                v1 = e;
                v2 = args{1};

                v1 = g.vertexlist(:,v1);
                v2 = g.vertexlist(:,v2);
            else
                % highlight_edge(E)
                e = g.edgelist(:,e);
                v1 = g.vertexlist(:,e(1));
                v2 = g.vertexlist(:,e(2));
            end

            % create the line properties for the edges
            lineprops = {
                'Color', opt.EdgeColor, ...
                'LineWidth', opt.EdgeThickness };

            if g.ndims == 3
                plot3([v1(1) v2(1)], [v1(2) v2(2)], [v1(3) v2(3)], lineprops{:});
            else
                plot([v1(1) v2(1)], [v1(2) v2(2)], lineprops{:});
            end
        end

        function highlight_path(g, path, varargin)
            %PGraph.highlight_path Highlight path
            %
            % G.highlight_path(P, OPTIONS) highlights the path defined by vector P
            % which is a list of vertex ids comprising the path.
            %
            % Options::
            %  'NodeSize',S          Size of vertex circle (default 12)
            %  'NodeFaceColor',C     Node circle color (default yellow)
            %  'NodeEdgeColor',C     Node circle edge color (default blue)
            %  'EdgeColor',C         Node circle edge color (default black)
            %  'EdgeThickness',T     Edge thickness (default 1.5)
            %
            % See also PGraph.highlight_node, PGraph.highlight_edge, PGraph.highlight_component.
            g.highlight_node(path, varargin{:});

            % highlight the edges
            for i=1:numel(path)-1
                v1 = path(i);
                v2 = path(i+1);
                g.highlight_edge(v1, v2, varargin{:});
            end
        end

        function dotfile(g, varargin)
            %Pgraph.dotfile Create a GraphViz dot file
            %
            % G.dotfile(filename, OPTIONS) creates the specified file which contains the
            % GraphViz code to represent the embedded graph.
            %
            % G.dotfile(OPTIONS) as above but outputs the code to the console.
            %
            % Options::
            %  'directed'    create a directed graph
            %
            % Notes::
            % - An undirected graph is default
            % - Use neato rather than dot to get the embedded layout

            opt.directed = false;

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

            if length(args) == 0
                fp = 1;
            else
                fp = fopen(args{1}, 'w');
            end

            if opt.directed
                fprintf(fp, 'digraph {\n');
            else
                fprintf(fp, 'graph {\n');
            end

            % add the nodes including name and position
            for i=1:g.n
                fprintf(fp, '"%s" [pos="%d,%d"]\n', g.names{i}, g.vertexlist(:,i));
            end

            % add the edges
            for i=1:g.ne
                edge = g.edgelist(:,i);
                if opt.directed
                    fprintf(fp, '"%s" -> "%s"\n', g.names{edge(1)}, g.names{edge(2)});
                else
                    fprintf(fp, '"%s" -- "%s"\n', g.names{edge(1)}, g.names{edge(2)});
                end
            end
            fprintf(fp, '}\n');
            if length(args) > 1
                fclose(fp);
            end
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% MATRIX REPRESENTATIONS
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        function L = laplacian(g)
            %Pgraph.laplacian Laplacian matrix of graph
            %
            % L = G.laplacian() is the Laplacian matrix (NxN) of the graph.
            %
            % Notes::
            % - L is always positive-semidefinite.
            % - L has at least one zero eigenvalue.
            % - The number of zero eigenvalues is the number of connected components
            %   in the graph.
            %
            % See also PGraph.adjacency, PGraph.incidence, PGraph.degree.

            L = g.degree() - (g.adjacency() > 0);
        end

        function D = degree(g)
            %Pgraph.degree Degree matrix of graph
            %
            % D = G.degree() is a diagonal matrix (NxN) where element D(i,i) is the number
            % of edges connected to vertex id i.
            %
            % See also PGraph.adjacency, PGraph.incidence, PGraph.laplacian.

            D = diag( g.connectivity() );
        end

        function A = adjacency(g)
            %Pgraph.adjacency Adjacency matrix of graph
            %
            % A = G.adjacency() is a matrix (NxN) where element A(i,j) is the cost
            % of moving from vertex i to vertex j.
            %
            % Notes::
            % - Matrix is symmetric.
            % - Eigenvalues of A are real and are known as the spectrum of the graph.
            % - The element A(I,J) can be considered the number of walks of one
            %   edge from vertex I to vertex J (either zero or one).  The element (I,J)
            %   of A^N are the number of walks of length N from vertex I to vertex J.
            %
            % See also PGraph.degree, PGraph.incidence, PGraph.laplacian.

            A = zeros(g.n, g.n);
            for i=1:g.n
                [n,c] = g.neighbours(i);
                for j=1:numel(n)
                    A(i,n(j)) = c(j);
                    A(n(j),i) = c(j);
                end
            end
        end

        function I = incidence(g)
            %Pgraph.degree Incidence matrix of graph
            %
            % IN = G.incidence() is a matrix (NxNE) where element IN(i,j) is
            % non-zero if vertex id i is connected to edge id j.
            %
            % See also PGraph.adjacency, PGraph.degree, PGraph.laplacian.
            I = zeros(g.n, g.ne);
            for i=1:g.n
                for n=g.edges(i)
                    I(i,n) = 1;
                end
            end
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% PATH PLANNING
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        function [path,cost] = Astar(g, vstart, vgoal, varargin)
            %PGraph.Astar path finding
            %
            % PATH = G.Astar(V1, V2) is the lowest cost path from vertex V1 to
            % vertex V2.  PATH is a list of vertices starting with V1 and ending
            % V2.
            %
            % [PATH,C] = G.Astar(V1, V2) as above but also returns the total cost
            % of traversing PATH.
            %
            % Notes::
            % - Uses the efficient A* search algorithm.
            % - The heuristic is the distance function selected in the constructor, it
            %   must be  admissible, meaning that it never overestimates the actual
            %   cost to get to the nearest goal node.
            %
            % References::
            % - Correction to "A Formal Basis for the Heuristic Determination of Minimum Cost Paths".
            %   Hart, P. E.; Nilsson, N. J.; Raphael, B.
            %   SIGART Newsletter 37: 28-29, 1972.
            %
            % See also PGraph.goal, PGraph.path.

            opt.directed = false;

            opt = tb_optparse(opt, varargin);

            % The set of vertices already evaluated.
            closedSet = [];

            % The set of tentative vertices to be evaluated, initially containing the start node
            openSet = [vstart];
            came_from(vstart) = 0;    % The map of navigated vertices.

            g_score(vstart) = 0;    % Cost from start along best known path.
            h_score(vstart) = g.distance(vstart, vgoal);
            % Estimated total cost from start to goal through y.
            f_score(vstart) = g_score(vstart) + h_score(vstart);

            while ~isempty(openSet)
                % current := the node in openSet having the lowest f_score[] value
                [~,k] = min(f_score(openSet));
                vcurrent = openSet(k);

                if vcurrent == vgoal
                    % we have arrived!
                    path = [];
                    p = vgoal;
                    while true
                        path = [p path];
                        p = came_from(p);
                        if p == 0
                            break;
                        end
                    end
                    if nargout > 1
                        cost = f_score(vgoal);
                    end
                    return
                end

                %remove current from openSet
                openSet = setdiff(openSet, vcurrent);
                %add current to closedSet
                closedSet = union(closedSet, vcurrent);

                if opt.directed
                    [neighbours,costs] = g.neighbours_out(vcurrent);

                else
                    [neighbours,costs] = g.neighbours(vcurrent);
                end

                for k=1:length(neighbours)

                    neighbour = neighbours(k);

                    if ismember(neighbour, closedSet)
                        continue;
                    end
                    tentative_g_score = g_score(vcurrent) + costs(k);

                    if ~ismember(neighbour, openSet)
                        %add neighbor to openSet
                        openSet = union(openSet, neighbour);
                        h_score(neighbour) = g.distance(neighbour, vgoal);
                        tentative_is_better = true;
                    elseif tentative_g_score < g_score(neighbour)
                        tentative_is_better = true;
                    else
                        tentative_is_better = false;
                    end
                    if tentative_is_better
                        % we found an edge that belongs to the path

                        % came_from is an array that records the path taken
                        %  - length g.n
                        %  -  came_from(A) = B means a path segment from A -> B
%                         if came_from(neighbour) > 0
%                             disp('problem')
%                         end
                        came_from(neighbour) = vcurrent;
                        g_score(neighbour) = tentative_g_score;
                        f_score(neighbour) = g_score(neighbour) + h_score(neighbour);
                    end
                end
            end
            path = [];
            cost = Inf;
        end

        function d = distance_metric(g, x1, x2)

            % distance between coordinates x1 and x2 using the relevant metric
            % x2 can be multiple points represented by multiple columns
            if isa(g.measure, 'function_handle')
                d = g.measure(x1(:), x2(:));
                if ~isscalar(d) && d <= 0
                    error('SMTB:PGraph:badresult', 'distance function must return a positive scalar');
                end
            else switch g.measure
                    case 'Euclidean'
                        d = colnorm( bsxfun(@minus, x1, x2) );

                    case 'SE2'
                        d = bsxfun(@minus, x1, x2) * g.dweight;
                        d(3,:) = angdiff(x1(3), x2(3,:));
                        d = colnorm( d );

                    case 'Lattice'
                        d = bsxfun(@minus, x1, x2) * g.dweight;
                        d(3,:) = angdiff(x1(3)*pi/2, x2(3,:)*pi/2);
                        d = colnorm( d );
                    otherwise
                        error('unknown distance measure', g.measure);
                end
            end
        end
    end %  methods

end % classdef

author: petercorke

repo: spatialmath-matlab


text:
classdef DropOut < dagnn.ElementWise
  properties
    rate = 0.5
    frozen = false
  end

  properties (Transient)
    mask
  end

  methods
    function outputs = forward(obj, inputs, params)
      if strcmp(obj.net.mode, 'test')
        outputs = inputs ;
        return ;
      end
      if obj.frozen & ~isempty(obj.mask)
        outputs{1} = vl_nndropout(inputs{1}, 'mask', obj.mask) ;
      else
        [outputs{1}, obj.mask] = vl_nndropout(inputs{1}, 'rate', obj.rate) ;
      end
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      if strcmp(obj.net.mode, 'test')
        derInputs = derOutputs ;
        derParams = {} ;
        return ;
      end
      derInputs{1} = vl_nndropout(inputs{1}, derOutputs{1}, 'mask', obj.mask) ;
      derParams = {} ;
    end

    % ---------------------------------------------------------------------
    function obj = DropOut(varargin)
      obj.load(varargin{:}) ;
    end

    function obj = reset(obj)
      reset@dagnn.ElementWise(obj) ;
      obj.mask = [] ;
      obj.frozen = false ;
    end
  end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function [varargout] = splitGradientAt(grad, timepoint, varargin)
%SplitGradient Splits a trapezoidal gradient into two extended trapezoids
%(currently shaped gradients) defined by the cut line.
%
%   [grads] = splitGradient(grad)
%   Returns the two gradient parts by cutting the original 'grad' at the
%   'timepoint' . For the input type 'trapezoid' the results are trtyurned
%   as extended trapezoids, for 'arb' as arbitrary gradient objects. The
%    delays in the individual gradient events are adapted such that
%   addGradients(...) produces an gradient equivalent to 'grad'.
%
%   See also  splitGradient makeExtendedTrapezoid makeTrapezoid
%             Sequence.addBlock  mr.opts
%
%   Maxim Zaitsev <maxim.zaitsev@uniklinik-freiburg.de>
%   Stefan Kroboth <stefan.kroboth@uniklinik-freiburg.de>

persistent parser

if isempty(parser)
    parser = inputParser;
    parser.FunctionName = 'splitGradientAt';
	parser.addRequired('grad', @isstruct);
    parser.addRequired('timepoint', @isnumeric);
    parser.addOptional('system', mr.opts(), @isstruct);
end
parse(parser, grad, timepoint, varargin{:});
opt = parser.Results;

gradRasterTime = opt.system.gradRasterTime;

% round the time point to the gradient raster;
timeindex = round(timepoint / gradRasterTime);
timepoint = timeindex * gradRasterTime;
timeindex = timeindex + 1; % convert to Matlab convention

ch = grad.channel;

if strcmp(grad.type, 'grad')
    % check if we have an arbitrary gradient or an exended trapezoid
    if abs(grad.tt(1)-0.5*gradRasterTime)<1e-10 && ...
       all(abs(grad.tt(2:end)-grad.tt(1:end-1)-gradRasterTime)<1e-10)
        % arbitrary gradient -- the most trivial conversion
        % if timepoint is out of range we have nothing to do
        if timeindex == 1 || timeindex >= length(grad.tt)
            varargout{1} = grad;
        else
            grad1=grad;
            grad2=grad;
            grad1.last=0.5*(grad.waveform(timeindex-1)+grad.waveform(timeindex)); % FIXME: retrive the double-sampling point (e.g. the corner of the trapezoid)
            grad2.first=grad1.last;
            grad2.delay=grad.delay + timepoint;
            grad1.tt=grad.tt(1:(timeindex-1));
            grad1.waveform=grad.waveform(1:(timeindex-1));
            grad2.tt=grad.tt(timeindex:end) - timepoint;
            grad2.waveform=grad.waveform(timeindex:end);

            if nargout==1
                varargout{1} = [grad1 grad2];
            else
                varargout{1} = grad1;
                varargout{2} = grad2;
            end
        end
        return; % early return
    else
        % we have an extended trapezoid -- excellent choice!
        times      = grad.tt;
        amplitudes = grad.waveform;
    end
elseif strcmp(grad.type, 'trap')
    grad.delay    = round(grad.delay   /gradRasterTime)*gradRasterTime; % MZ: was ceil
    grad.riseTime = round(grad.riseTime/gradRasterTime)*gradRasterTime; % MZ: was ceil
    grad.flatTime = round(grad.flatTime/gradRasterTime)*gradRasterTime; % MZ: was ceil
    grad.fallTime = round(grad.fallTime/gradRasterTime)*gradRasterTime; % MZ: was ceil

    % prepare the extended trapezoid structure
    if grad.flatTime == 0
        times      = [0 grad.riseTime  grad.riseTime+grad.fallTime];
        amplitudes = [0 grad.amplitude 0];
    else
        times      = [0 grad.riseTime  grad.riseTime+grad.flatTime grad.riseTime+grad.flatTime+grad.fallTime];
        amplitudes = [0 grad.amplitude grad.amplitude              0];
    end
else
    error('Splitting of unsupported event.');
end

% if the cutline is behind the gradient there is no second gradient to create
if timepoint >= grad.delay+times(end)
    error('trying to place the splitting time point after the end of the gradient');
end

% now we have everything in the extended trapezoid structure

% if the cutline goes through the delay we need special treatment
if timepoint < grad.delay
    times=[0 grad.delay+times];
    amplitudes = [0 amplitudes];
    grad.delay=0;
else
    timepoint = timepoint - grad.delay;
end

% sample at timepoint
amp_tp=interp1(times, amplitudes, timepoint, 'linear'); % MZ: interp1() is not OK here for the corner situation TODO: fixme! (e.g. by restoring the corners as done in waveforms_and_times())
% split the data
teps=1e-10; % we need this because of the rounding problems
times1 = [ times(times<timepoint-teps) timepoint ];
amplitudes1 = [ amplitudes(times<timepoint-teps) amp_tp ];
times2 = [ timepoint times(times>timepoint+teps) ] - timepoint;
amplitudes2 = [ amp_tp amplitudes(times>timepoint+teps) ];

% recreate gradients
grad1 = mr.makeExtendedTrapezoid(ch, opt.system, 'times', times1,...
                                  'amplitudes', amplitudes1, ...
                                  'skip_check', true);
grad1.delay = grad.delay;
grad2 = mr.makeExtendedTrapezoid(ch, opt.system, 'times', times2,...
                                  'amplitudes', amplitudes2, ...
                                  'skip_check', true);
grad2.delay = timepoint + grad.delay;

%grads = [grad1 grad2];
if nargout==1
    varargout{1} = [grad1 grad2];
else
    varargout{1} = grad1;
    varargout{2} = grad2;
end

end

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
clc;
clear;
close all;
addpath(genpath('image'));
addpath(genpath('whyte_code'));
addpath(genpath('cho_code'));
opts.prescale = 1; %%downsampling
opts.xk_iter = 5; %% the iterations
opts.gamma_correct = 1.0;
opts.k_thresh = 20;

% modify below line for the input dir for blurry image
inputdir = 'x2/';
% modify below line for the output dir for the kernels
outputdir = 'x2results/';

dirs = dir(inputdir);
L = length(dirs);
for l = 3:L
    picfiles = dir(strcat(inputdir, strcat(dirs(l).name, '/*.jpg')));
    LL = length(picfiles);
    for kk = 1:LL
      try
        opts.kernel_size = 25;  saturation = 0;
        lambda_dark = 4e-3; lambda_grad = 4e-3;opts.gamma_correct = 1.0;
        lambda_tv = 0.003; lambda_l0 = 5e-4; weight_ring = 1;
        img = imread([inputdir dirs(l).name '/' picfiles(kk).name]);
        y = img;
        isselect = 0; %false or true
        if isselect ==1
            figure, imshow(y);
            %tips = msgbox('Please choose the area for deblurring:');
            fprintf('Please choose the area for deblurring:\n');
            h = imrect;
            position = wait(h);
            close;
            B_patch = imcrop(y,position);
            y = (B_patch);
        else
            y = y;
        end
        if size(y,3)==3
            yg = im2double(rgb2gray(y));
        else
            yg = im2double(y);
        end
        [kernel, interim_latent] = blind_deconv(yg, lambda_dark, lambda_grad, opts);

        % for plotting
        k = kernel - min(kernel(:));
        k = k./max(k(:));

        % save kernel
      save(sprintf('%s%s%s%s%s',outputdir,dirs(l).name,'/',string(extractBetween(picfiles(kk).name,1,length(picfiles(kk).name)-4)),'_kernel.mat'), 'kernel');
        % uncomment the following line to generate png files for kernels
        % imwrite(k,sprintf('%s%s%s%s%s',outputdir,dirs(l).name,'/',string(extractBetween(picfiles(kk).name,1,length(picfiles(kk).name)-4)),'_kernel.png'));
      catch
        1
      end

    end
end

author: IVRL

repo: Kernel-Modeling-Super-Resolution


text:
function varargout = ft_plot_lay(varargin)

% This function is a backward compatibility wrapper for existing MATLAB scripts
% that call a function that is not part of the FieldTrip toolbox any more.
%
% Please update your code to make it future-proof.

oldname = mfilename;
newname = 'ft_plot_layout';

ft_warning('%s is only a backward compatibility wrapper, which will soon be removed. Please call %s instead.', upper(oldname), upper(newname));

[varargout{1:nargout}] = feval(newname, varargin{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function ct = matRad_calcHU(ct)
% matRad function to calculate Hounsfield units from a dicom ct
% that originally uses intensity values
%
% call
%   ct = matRad_calcHU(ct)
%
% input
%   ct: unprocessed dicom ct data which are stored as intensity values (IV)
%
%                      HU = IV * slope + intercept
%
% output
%   ct:                 ct struct with cube with HU
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1:ct.numOfCtScen
    ct.cubeHU{i} = double(ct.cubeIV{i}) * double(ct.dicomInfo.RescaleSlope) + double(ct.dicomInfo.RescaleIntercept);
end

ct = rmfield(ct,'cubeIV');

end

repo: matRad


text:
function pass = test_syntax( pref )
% Check the Chebfun2v constructor for different syntax.
% Alex Townsend, March 2013.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e5 * pref.cheb2Prefs.chebfun2eps;

for jj = 1 : 2

    f = @(x,y) jj*sin(x.*y);  % simple function.
    g = @(x,y ) jj*cos(x.*y);

    fd = diskfun(f);
    gd = diskfun(g);

    F1 = diskfunv( f, g);
    F2 = diskfunv( fd, gd);
    F3 =  [ fd; gd];

    pass(1, jj) = ( norm(F1 - F2) < tol );
    pass(2, jj) = ( norm(F2 - F3) < tol );

end
pass = pass(:)';

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [q num den]= get_q_RANSAC(N, N_I, k)

% [q num den] = get_q_RANSAC(N, N_I, k)
%
% DESC:
% calculates the probability q
%
% AUTHOR
% Marco Zuliani - marco.zuliani@gmail.com
%
% VERSION:
% 1.0
%
% INPUT:
% N                 = number of elements
% N_I               = number of inliers
% k                 = cardinality of the MSS
%
% OUTPUT:
% q                 = probability
% num, den          = q = num / den

if (k > N_I)
    error('k should be less or equal than N_I')
end;

if (N == N_I)
    q = 1;
    return;
end;

num = N_I:-1:N_I-k+1;
den = N:-1:N-k+1;
q = prod(num./den);

return

author: RANSAC

repo: RANSAC-Toolbox


text:
%JOY2TR Update transform from joystick
%
% T = JOY2TR(T, OPTIONS) updates the SE(3) homogeneous transform (4x4)
% according to spatial velocities sourced from a connected joystick device.
%
% Options::
% 'delay',D     Pause for D seconds after reading (default 0.1)
% 'scale',S     A 2-vector which scales joystick translational and
%               rotational to rates (default [0.5m/s, 0.25rad/s])
% 'world'       Joystick motion is in the world frame
% 'tool'        Joystick motion is in the tool frame (default)
% 'rotate',R    Index of the button used to enable rotation (default 7)
%
% Notes::
% - Joystick axes 0,1,3 map to X,Y,Z or R,P,Y motion.
% - A joystick button enables the mapping to translation OR rotation.
% - A 'delay' of zero means no pause
% - If 'delay' is non-zero 'scale' maps full scale to m/s or rad/s.
% - If 'delay' is zero 'scale' maps full scale to m/sample or rad/sample.
%
% See also joystick.

% Copyright (C) 1993-2017, 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

function Tout = joy2tr(T, varargin)

    % parse the options
    opt.delay = 0.1;
    opt.rotate = 7;
    opt.scale = [0.5 0.25];
    opt.frame = {'tool', 'world'};
    opt.min = 0.006;

    opt = tb_optparse(opt, varargin);

    % get the raw joystick output
    [j,b] = joystick();

    % update a 6-vector of translational and rotational velocities
    vel = zeros(1,6);
    if b(7) == 0
        % translation mode
        vel(1:2) = j(1:2);
        vel(3) = j(4);
    else
        % rotation mode
        vel(4:5) = j(1:2);
        vel(6) = j(4);
    end

    % values below threshold are set to zero
    vel(abs(vel) < opt.min) = 0;

    % apply scaling factors
    scale = opt.scale;

    if opt.delay > 0
        % normalize scaling factors by time
        pause(opt.delay);
        scale = scale * opt.delay;
    end

    vel(1:3) = vel(1:3) * scale(1);
    vel(4:6) = vel(4:6) * scale(2);

    % compute the incremental motion
    dT = transl(vel(1:3)) * rpy2tr(vel(4:6));

    % and apply it to the input transform
    switch opt.frame
        case 'tool'
            Tout = T * dT;
        case 'world'
            Tout = dT * T;
    end

    % normalize just to be safe
    Tout = trnorm(Tout);

repo: robotics-toolbox-matlab


text:
function [result]=sv_NodeCalcSchuster(params,mCatalog)
% function [result]=sv_NodeCalcSchuster(params,mCatalog)
% ----------------------------------------------------
% Function to calculate Schuster's test and fast Mc (not EMR method)
% estimates; calculates also uncertainties for Schuster's test
%
% Incoming variables:
% mCatalog     : current earthquake catalog
% params       : See sv_calcMc for parameters
%
% Outgoing variable:
% result.fMc_max    : magnitude of completeness using maximum cuvature approach
% result.fMc_90     : magnitude of completeness using GR-law fit at 90% goodness fit level
% result.fMc_95     : magnitude of completeness using GR-law fit at 95% goodness fit level
% result.fMc_com    : magnitude of completeness using GR-law fit at 90% goodness fit level,
%                     GR-law fit at 95% goodness fit level, and maximum cuvature approach
% result.fMcSch      : magnitude of completeness using Schusters method
% result.flogProbSch : log10(Probability) to obtain vector of length >= R from a random walkout
% result.v1Sigma     : 1Sigma (16 and 84 percentiles) for probability to obtain vector of length >= R
%                      from a random walkout from bootstrapping the catalog times
% result.flogv1Sigma : 1Sigma (16 and 84 percentiles) for log10(probability) to obtain vector of length >= R
%                      from a random walkout from bootstrapping the catalog times
% result.v1Srange    : 1Sigma range (84-16 percentile)
% result.flogv1Srange: 1Sigma log10 range (84-16 percentile)
%
% Author: J. Woessner
% woessner@seismo.ifg.ethz.ch
% last update: 17.04.03

% Init variable
result=[];

% Determine Mc by maximum curvature
nCalculateMC = 1;
result.fMc_max = calc_Mc(mCatalog, nCalculateMC);
if isempty(result.fMc_max)
    result.fMc_max = NaN;
end

% Determine Mc by goodness of fit 90%
nCalculateMC = 3;
result.fMc_90 = calc_Mc(mCatalog, nCalculateMC);
if isempty(result.fMc_90)
    result.fMc_90 = NaN;
end

% Determine Mc by goodness of fit 95%
nCalculateMC = 4;
result.fMc_95 = calc_Mc(mCatalog, nCalculateMC);
if isempty(result.fMc_95)
    result.fMc_95 = NaN;
end

% Determine Mc by best combination
nCalculateMC = 5;
result.fMc_com = calc_Mc(mCatalog, nCalculateMC);
if isempty(result.fMc_com)
    result.fMc_com = NaN;
end

% % Mc determination by modelling entire magnitude range using a NORMAL CDF
% [mResult result.fProbMcNorm result.fMcNorm fMu fSigma mDatPredBest vPredBest result.fBvalue] = calc_McCdfnormal(mCatalog, params.fBinning);
%
% if (isempty(result.fProbMcNorm) | isempty(result.fMcNorm))
%     result.fProbMcNorm = NaN;
%     result.fMcNorm = NaN;
% end
% if (isempty(result.fBvalue))
%     result.fBvalue = NaN;
% end

% Mc determination by Schuster
[mResult fMcSch,  fProbability] =  calc_SchusterMc(mCatalog,0.3);
[mWalkout fR95 fProbSch PHI,  R] = calc_Schusterwalk(mCatalog);
try
    result.fMcSch = fMcSch;
    result.flogProbSch = log10(fProbSch);
catch
    result.fMcSch = NaN;
    result.flogProbSch = NaN;
end

% Determine uncertainty of Schuster test
[vPerc v1Sigma,  fStdProb] = calc_BstSchuster(mCatalog,params.fBinning,params.fBstnum);
try
    result.v1Sigma = v1Sigma;
    result.flogv1Sigma = log10(v1Sigma);
    result.v1Srange = (v1Sigma(1,2) - v1Sigma(1,1));
    result.flogv1Srange = log10(result.v1Srange);
catch
    result.fv1Sigma= NaN;
    result.flogv1Sigma = [NaN NaN];
    result.v1Srange = NaN;
    result.flogv1Srange = NaN;
end

author: CelsoReyes


text:
function [ y, dzdw,dzdb,opts] = linear_layer( I,weight,bias,dzdy,opts)
%FAST_MLP_LAYER Summary of this function goes here
%   Detailed explanation goes here
%I: input_dim x batch_size

dzdw=[];
dzdb=[];

if exist('opts','var')&&isfield(opts,'use_nntoolbox')&&opts.use_nntoolbox==1
    [out,in]=size(weight);
    batch_size=size(I,2);
    %use nntoolbox
    if ~isfield(opts,'layer')||length(opts.layer)<opts.current_layer||~isfield(opts.layer{opts.current_layer},'fc_nntb')
        product_info=ver('nnet');
        opts.nnet_ver=str2double(product_info.Version);
        if opts.nnet_ver<11
            opts.layer{opts.current_layer}.fc_nntb=nnet.internal.cnn.layer.Convolution2D('conv2d_nntb', [1,1], in ,out, [1,1], [0,0]);
        else
            opts.layer{opts.current_layer}.fc_nntb=nnet.internal.cnn.layer.Convolution2D('conv2d_nntb', [1,1], in ,out, [1,1],'manual', [0,0,0,0]);
        end
        if opts.use_gpu
            opts.layer{opts.current_layer}.fc_nntb = setupForGPUPrediction(opts.layer{opts.current_layer}.fc_nntb);
        else
            opts.layer{opts.current_layer}.fc_nntb =setupForHostPrediction(opts.layer{opts.current_layer}.fc_nntb);
        end
    end
    fc_nntb=opts.layer{opts.current_layer}.fc_nntb;
    if ~isempty(bias)
        fc_nntb.Bias.Value=permute(bias(:),[3,4,1,2]);
    else
        fc_nntb.Bias.Value=zeros(1,1,out,1);
    end
    fc_nntb.Weights.Value=permute(weight,[4,3,2,1]);
    I=permute(I,[3,4,1,2]);

    if isempty(dzdy)
        y=fc_nntb.forward(I);
        y=permute(y,[3,4,1,2]);
    else

        dzdy=permute(dzdy,[3,4,1,2]);
        if opts.nnet_ver<11
            y = fc_nntb.backward(  I, [], dzdy, [] );
            gradients = fc_nntb.gradients(I, dzdy);
        else
            [y,gradients] = fc_nntb.backward( I, [], dzdy, [] );
        end

        dzdw=gradients{1}./batch_size;
        dzdb=gradients{2}./batch_size;
        y=permute(y,[3,4,1,2]);
        dzdw=permute(dzdw,[4,3,2,1]);
        dzdb=reshape(dzdb,size(bias));
    end

    return;
end

if isempty(dzdy)
    %forward mode

    y=weight*I;

    if ~isempty(bias)
        y=y+bias;
    end

else
    %backward mode

    y=weight'*dzdy;
    if ~isempty(bias)
        dzdb=mean(dzdy,2);%minibatch averaging
    end
    dzdy=permute(dzdy,[1,3,2]);
    I=permute(I,[3,1,2]);
    dzdw=dzdy.*I;
    dzdw=mean(dzdw,3);

end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
%% GMWB Script for GPUs
% This script is adapted from the code for the MathWorks webinar "Modeling
% Variable Annuities with MATLAB" found at:
%
% http://www.mathworks.com/matlabcentral/fileexchange/26960-modeling-variable-annuities-with-matlab
%
% It is designed to highlight some of the coding considerations that must
% be addressed if one wishes to move such an analysis from a CPU
% environment onto a GPU.
%
% Original: Yi Wang, MathWorks, 2010
% Adapted: Michael Weidman, Quantitative Support Services, Ltd., 2013

%% 1. Set Parameters and Load Data

annualFee               = 0.005;
riskFreeRate            = 0.05;
annualWithdrawalRate    = 0.07;
nYears                  = 20;

nTrials                 = 1e5;

% Investment Portfolio
Ticker = {'MMM', 'AA', 'AXP', 'T', 'BAC', 'BA', 'CAT', 'CVX', 'CSCO', 'KO'};

% Range for historical data
FromDate    = '01/01/2000';
ToDate      = '01/01/2010';
Period      = 'm'; % Monthly data

% Assume we're holding 10 shares of each stock in our portfolio
Holdings = 10*ones(length(Ticker), 1);

TimeSeries = getEquityData(Ticker, FromDate, ToDate, Period);

dispResults('header')

%% 2. Price GMWB I: Original Code

tic

assetPrice = simAssetPrice_orig(TimeSeries, Period, nYears, nTrials);

[AccountVal, PayoutVal, FeeVal] = calcValuePayoutAndFees_orig( ...
    Holdings, assetPrice, annualWithdrawalRate, annualFee);

toc1 = toc;

cost = pvvar(mean(PayoutVal, 2), riskFreeRate);
fee = pvvar(mean(FeeVal, 2), riskFreeRate);
probRuin = calcProbRuin(AccountVal);

dispResults('Original', cost, fee, probRuin, toc1)

%% 3. Price GMWB II: Vectorized Code, GPU

tic

assetPrice = simAssetPrice_orig(TimeSeries, Period, nYears, nTrials);

HoldingsG               = gpuArray(Holdings);
assetPriceG             = gpuArray(assetPrice);
annualWithdrawalRateG   = gpuArray(annualWithdrawalRate);
annualFeeG              = gpuArray(annualFee);

[AccountValG, PayoutValG, FeeValG] = calcValuePayoutAndFees_GPU( ...
    HoldingsG, assetPriceG, annualWithdrawalRateG, annualFeeG);

AccountVal  = gather(AccountValG);
PayoutVal   = gather(PayoutValG);
FeeVal      = gather(FeeValG);

toc1 = toc;

cost = pvvar(mean(PayoutVal, 2), riskFreeRate);
fee = pvvar(mean(FeeVal, 2), riskFreeRate);
probRuin = calcProbRuin(AccountVal);

dispResults('Vectorized, GPU', cost, fee, probRuin, toc1)

%% 4. Price GMWB III: Vectorized Code, CPU

tic

assetPrice = simAssetPrice_orig(TimeSeries, Period, nYears, nTrials);

[AccountVal, PayoutVal, FeeVal] = calcValuePayoutAndFees_VEC( ...
    Holdings, assetPrice, annualWithdrawalRate, annualFee);

toc1 = toc;

cost = pvvar(mean(PayoutVal, 2), riskFreeRate);
fee = pvvar(mean(FeeVal, 2), riskFreeRate);
probRuin = calcProbRuin(AccountVal);

dispResults('Vectorized, CPU', cost, fee, probRuin, toc1)

%% 5. Price GMWB IV: GPU for RNG, Vectorized CPU for rest

tic

assetPriceG = simAssetPrice_2GPU(TimeSeries, Period, nYears, nTrials);
assetPrice = gather(assetPriceG);

[AccountVal, PayoutVal, FeeVal] = calcValuePayoutAndFees_VEC( ...
    Holdings, assetPrice, annualWithdrawalRate, annualFee);

toc1 = toc;

cost = pvvar(mean(PayoutVal, 2), riskFreeRate);
fee = pvvar(mean(FeeVal, 2), riskFreeRate);
probRuin = calcProbRuin(AccountVal);

dispResults('Best: GPU & CPU', cost, fee, probRuin, toc1)

%% 6. Price GMWB V: CPU only with new RNG and Vectorized code
% This optional section confirms that the GPU + CPU method in the previous
% section is the fastest: Using the GPU to generate random numbers should
% be slightly faster that this (CPU-only) section.

tic

assetPrice = simAssetPrice_2(TimeSeries, Period, nYears, nTrials);

[AccountVal, PayoutVal, FeeVal] = calcValuePayoutAndFees_VEC( ...
    Holdings, assetPrice, annualWithdrawalRate, annualFee);

toc1 = toc;

cost = pvvar(mean(PayoutVal, 2), riskFreeRate);
fee = pvvar(mean(FeeVal, 2), riskFreeRate);
probRuin = calcProbRuin(AccountVal);

dispResults('2nd Best: CPU-only', cost, fee, probRuin, toc1)

author: Sable

repo: mcbench-benchmarks


text:
function results = eco_vot_deep(seq, res_path, bSaveImage, parameters)

params.use_gpu = false;
params.gpu_id = 1;

% Feature specific parameters
% hog_params.cell_size = 4;
% hog_params.compressed_dim = 10;

% grayscale_params.colorspace='gray';
% grayscale_params.cell_size = 1;

% cn_params.tablename = 'CNnorm';
% cn_params.useForGray = false;
% cn_params.cell_size = 4;
% cn_params.compressed_dim = 3;
%
% ic_params.tablename = 'intensityChannelNorm6';
% ic_params.useForColor = false;
% ic_params.cell_size = 4;
% ic_params.compressed_dim = 3;

% cnn_params.nn_name = 'imagenet-vgg-m-2048_new.mat'; % Name of the network
cnn_params.nn_name = 'UDT_Unsupervised.mat';

cnn_params.output_layer = [4];               % Which layers to use
cnn_params.downsample_factor = [4];           % How much to downsample each output layer
cnn_params.compressed_dim = [32];            % Compressed dimensionality of each output layer
cnn_params.input_size_mode = 'adaptive';        % How to choose the sample size
cnn_params.input_size_scale = 1;

% Which features to include
params.t_features = {
    struct('getFeature',@get_cnn_layers, 'fparams',cnn_params),...
    ...struct('getFeature',@get_colorspace, 'fparams',grayscale_params),...
    ...struct('getFeature',@get_fhog,'fparams',hog_params),...
    ...struct('getFeature',@get_table_feature, 'fparams',cn_params),...
    ...struct('getFeature',@get_table_feature, 'fparams',ic_params),...
};

% Global feature parameters1s
params.t_global.normalize_power = 2;    % Lp normalization with this p
params.t_global.normalize_size = true;  % Also normalize with respect to the spatial size of the feature
params.t_global.normalize_dim = true;   % Also normalize with respect to the dimensionality of the feature

% Image sample parameters
params.search_area_shape = 'square';    % The shape of the samples
params.search_area_scale = 4.0;         % The scaling of the target size to get the search area
params.min_image_sample_size = 200^2;   % Minimum area of image samples
params.max_image_sample_size = 250^2;   % Maximum area of image samples

% Detection parameters
params.refinement_iterations = 1;       % Number of iterations used to refine the resulting position in a frame
params.newton_iterations = 5;           % The number of Newton iterations used for optimizing the detection score
params.clamp_position = false;          % Clamp the target position to be inside the image

% Learning parameters
params.output_sigma_factor = 1/12;		% Label function sigma
params.learning_rate = 0.011;	 	    % Learning rate
params.nSamples = 50;                   % Maximum number of stored training samples
params.sample_replace_strategy = 'lowest_prior';    % Which sample to replace when the memory is full
params.lt_size = 0;                     % The size of the long-term memory (where all samples have equal weight)
params.train_gap = 5;                   % The number of intermediate frames with no training (0 corresponds to training every frame)
params.skip_after_frame = 10;           % After which frame number the sparse update scheme should start (1 is directly)
params.use_detection_sample = true;     % Use the sample that was extracted at the detection stage also for learning

% Factorized convolution parameters
params.use_projection_matrix = true;    % Use projection matrix, i.e. use the factorized convolution formulation
params.update_projection_matrix = true; % Whether the projection matrix should be optimized or not
params.proj_init_method = 'pca';        % Method for initializing the projection matrix
params.projection_reg = 2e-7;           % Regularization paremeter of the projection matrix

% Generative sample space model parameters
params.use_sample_merge = true;                 % Use the generative sample space model to merge samples
params.sample_merge_type = 'Merge';        % Strategy for updating the samples
params.distance_matrix_update_type = 'exact';  % Strategy for updating the distance matrix
params.neglect_higher_frequency = false;        % Neglect hiigher frequency components in the distance comparison for speed

% Conjugate Gradient parameters
params.CG_iter = 5;                     % The number of Conjugate Gradient iterations in each update after the first frame
params.init_CG_iter = 10*20;            % The total number of Conjugate Gradient iterations used in the first frame
params.init_GN_iter = 10;               % The number of Gauss-Newton iterations used in the first frame (only if the projection matrix is updated)
params.CG_use_FR = false;               % Use the Fletcher-Reeves (true) or Polak-Ribiere (false) formula in the Conjugate Gradient
params.CG_standard_alpha = true;        % Use the standard formula for computing the step length in Conjugate Gradient
params.CG_forgetting_rate = 60;	 	 	% Forgetting rate of the last conjugate direction
params.precond_data_param = 0.75;	 	% Weight of the data term in the preconditioner
params.precond_reg_param = 0.2;	 	    % Weight of the regularization term in the preconditioner
params.precond_proj_param = 40;         % Weight of the projection matrix part in the preconditioner

% Regularization window parameters
params.use_reg_window = true;           % Use spatial regularization or not
params.reg_window_min = 1e-4;			% The minimum value of the regularization window
params.reg_window_edge = 10e-3;         % The impact of the spatial regularization
params.reg_window_power = 2;            % The degree of the polynomial to use (e.g. 2 is a quadratic window)
params.reg_sparsity_threshold = 0.15;   % A relative threshold of which DFT coefficients that should be set to zero

% Interpolation parameters
params.interpolation_method = 'bicubic';    % The kind of interpolation kernel
params.interpolation_bicubic_a = -0.75;     % The parameter for the bicubic interpolation kernel
params.interpolation_centering = true;      % Center the kernel at the feature sample
params.interpolation_windowing = false;     % Do additional windowing on the Fourier coefficients of the kernel

% Scale parameters for the translation model
% Only used if: params.use_scale_filter = false

% params.number_of_scales = 3;            % Number of scales to run the detector
% params.scale_step = 1.02;               % The scale factor

% Scale filter parameters
% Only used if: params.use_scale_filter = true
params.use_scale_filter = true;        % Use the fDSST scale filter or not (for speed)
params.scale_sigma_factor = 1/16;       % Scale label function sigma
params.scale_learning_rate = 0.025;		% Scale filter learning rate
params.number_of_scales_filter = 17;    % Number of scales
params.number_of_interp_scales = 33;    % Number of interpolated scales
params.scale_model_factor = 1.0;        % Scaling of the scale model
params.scale_step_filter = 1.02;        % The scale factor for the scale filter
params.scale_model_max_area = 32*16;    % Maximume area for the scale sample patch
params.scale_feature = 'HOG4';          % Features for the scale filter (only HOG4 supported)
params.s_num_compressed_dim = 'MAX';    % Number of compressed feature dimensions in the scale filter
params.lambda = 1e-2;					% Scale filter regularization
params.do_poly_interp = true;           % Do 2nd order polynomial interpolation to obtain more accurate scale

% Other parameters
params.visualization = 0;               % Visualiza tracking and detection scores
params.debug = 0;                       % Do full debug visualization

% Initialize
params.seq = seq;

% Run tracker
results = tracker(params);

repo: UDT


text:
function SubOffspring = CreateOff(Problem,Population,SubPopulation,RMP)

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

    %% Parent selection
    Parent11 = [];
    Parent21 = [];
    Parent12 = [];
    Parent22 = [];
    for i = 1 : floor(length(Population)/2)
        P1  = Population(i);
        P2  = Population(i+floor(length(Population)/2));
        rmp = RMP(P1.dec(end),P2.dec(end));
        if (P1.dec(end) == P2.dec(end)) || (rand<rmp)
            Parent11 = [Parent11,P1];
            Parent21 = [Parent21,P2];
        else
            Parent12 = [Parent12,P1,P2];
            Parent22 = [Parent22,SubPopulation{P1.dec(end)}(randi(end)),SubPopulation{P2.dec(end)}(randi(end))];
        end
    end

    %% Offspring generation
    if ~isempty(Parent11)
        Parent1Dec     = Parent11.decs;
        Parent2Dec     = Parent21.decs;
        OffDec1        = OperatorGA(Problem,[Parent1Dec;Parent2Dec]);
        OffDec1(:,end) = [Parent1Dec(:,end);Parent2Dec(:,end)];
    else
        OffDec1 = [];
    end
    if ~isempty(Parent12)
        Parent1Dec     = Parent12.decs;
        Parent2Dec     = Parent22.decs;
        OffDec2        = OperatorGAhalf(Problem,[Parent1Dec;Parent2Dec]);
        OffDec2(:,end) = Parent1Dec(:,end);
    else
        OffDec2 = [];
    end
    SubOffspring = Divide(Problem.Evaluation([OffDec1;OffDec2]),length(Problem.SubD));
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%% ALL NONLINEAR VERSION
clc
clear
% Objective & Gradient
 fun = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
 grad = @(x)[2*x(1)-400*x(1)*(x(2)-x(1)^2)-2,200*x(2)-200*x(1)^2];

 % Constraints
A = [-1 1; 1 1];
rl = [-Inf;5];
ru = [-1;5];
lb = [0;0]; ub = [4;4];

% Nonlinear Constraint, Jacobian & Structure
 nlcon = @(x) A*x;
 nljac = @(x) sparse(A);
 jacstr = @() sparse(double(A~=0));

% Starting Guess
 x0 = [2;2];

% Build Function Structure
 funcs.objective = fun;
 funcs.gradient = grad;
 funcs.constraints = nlcon;
 funcs.jacobian = nljac;
 funcs.jacobianstructure = jacstr;

% Build Options Structure
 opts.lb = lb;
 opts.ub = ub;
 opts.cl = rl;
 opts.cu = ru;
 opts.ipopt.hessian_approximation = 'limited-memory';

% Call IPOPT
[x,output] = bonmin(x0,funcs,opts)

%% ALL LINEAR VERSION
clc
clear all
% Objective & Gradient
 fun = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
 grad = @(x)[2*x(1)-400*x(1)*(x(2)-x(1)^2)-2,200*x(2)-200*x(1)^2];

 % Constraints
A = [-1 1; 1 1];
rl = [-Inf;5];
ru = [-1;5];
lb = [0;0]; ub = [4;4];

% Starting Guess
 x0 = [2;2];

% Build Function Structure
 funcs.objective = fun;
 funcs.gradient = grad;

% Build Options Structure
 opts.lb = lb;
 opts.ub = ub;
 opts.A = sparse(A);
 opts.rl = rl;
 opts.ru = ru;
 opts.ipopt.hessian_approximation = 'limited-memory';
 opts.display = 2;

% Call IPOPT
[x,output] = bonmin(x0,funcs,opts)

%% NL INEQ, LIN EQ
clc
clear all
% Objective & Gradient
 fun = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
 grad = @(x)[2*x(1)-400*x(1)*(x(2)-x(1)^2)-2,200*x(2)-200*x(1)^2];

 % Constraints
Aeq = [1 1];
rl = 5;
ru = 5;
lb = [0;0]; ub = [4;4];

% Nonlinear Constraint, Jacobian & Structure
 A = [-1 1];
 nlcon = @(x) [A]*x;
 nljac = @(x) sparse([A]);
 jacstr = @() sparse(double([A]~=0));
 cl = -Inf;
 cu = -1;

% Starting Guess
 x0 = [2;2];

% Build Function Structure
 funcs.objective = fun;
 funcs.gradient = grad;
 funcs.constraints = nlcon;
 funcs.jacobian = nljac;
 funcs.jacobianstructure = jacstr;

% Build Options Structure
 opts.lb = lb;
 opts.ub = ub;
 opts.A = sparse(Aeq);
 opts.rl = rl;
 opts.ru = ru;
 opts.cl = cl;
 opts.cu = cu;
 opts.ipopt.hessian_approximation = 'limited-memory';
  opts.display = 2;

% Call IPOPT
% for i = 1:1000
    [x,output] = bonmin(x0,funcs,opts)
% end

%% NL EQ, LIN INEQ
clc
clear all
% Objective & Gradient
 fun = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
 grad = @(x)[2*x(1)-400*x(1)*(x(2)-x(1)^2)-2,200*x(2)-200*x(1)^2];

 % Constraints
A = [-1 1];
rl = -Inf;
ru = -1;
lb = [0;0]; ub = [4;4];

% Nonlinear Constraint, Jacobian & Structure
 Aeq = [1 1];
 nlcon = @(x) [Aeq]*x;
 nljac = @(x) sparse([Aeq]);
 jacstr = @() sparse(double([Aeq]~=0));
 cl = 5;
 cu = 5;

% Starting Guess
 x0 = [2;2];

% Build Function Structure
 funcs.objective = fun;
 funcs.gradient = grad;
 funcs.constraints = nlcon;
 funcs.jacobian = nljac;
 funcs.jacobianstructure = jacstr;

% Build Options Structure
 opts.lb = lb;
 opts.ub = ub;
 opts.A = sparse(A);
 opts.rl = rl;
 opts.ru = ru;
 opts.cl = cl;
 opts.cu = cu;
 opts.ipopt.hessian_approximation = 'limited-memory';
  opts.display = 1;

% Call IPOPT
% for i = 1:1000
    [x,output] = bonmin(x0,funcs,opts)
% end

%%
clc
%Objective
obj = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
% Constraints
A = [-1 1];
b = -1;
Aeq = [1 1];
beq = 5;
lb = [0;0]; ub = [4;4];
% Solve
x0 = [2;2];
opts = optiset('solver','bonmin','display','iter');
Opt = opti('obj',obj,'ndec',2,'bounds',lb,ub,'ineq',A,b,'eq',Aeq,beq,'ivars',2,'options',opts)
[x,fval,ef,info] = solve(Opt,x0)
%Plot
% plot(Opt,[],1)

%%
% Objective & Gradient
fun = @(x) -sum(x);
grad = @(x) [-1 -1 -1 -1];

% Linear constraints
n = 4;
A = randn(n,4);
rl = zeros(n,1);
ru = zeros(n,1);

% Nonlinear Constraint, Jacobian & Structure
nlcon = @(x) x(:).^4;
nljac = @(x) sparse(diag(4*x(:).^3));
jacstr = @() speye(4);
cl = [-inf;-inf;-inf;-inf];
cu = [1;1;1;1];

% Starting Guess
x0 = [0;0;0;0];
% Build Function Structure
funcs.objective = fun;
funcs.gradient = grad;
funcs.constraints = nlcon;
funcs.jacobian = nljac;
funcs.jacobianstructure = jacstr;

% Build Options Structure
lb = [-10;-10;-10;-10];
ub = [40;40;40;40];
opts.lb = lb;
opts.ub = ub;
opts.rl = rl;
opts.ru = ru;
opts.cl = cl;
opts.cu = cu;
opts.A = sparse(A);
opts.ipopt.hessian_approximation = 'limited-memory';
opts.display = 2;

% Call IPOPT
[x,output] = bonmin(x0,funcs,opts)

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function test_bug2069

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_appendfreq ft_freqstatistics ft_timelockstatistics

%% simulate some data

tl = cell(1,20);
tl_shuff = cell(1,20);
freq = cell(1,20);
freq_shuff = cell(1,20);

for k = 1:20
  % tl will hold timelocked structures with no channel order shuffling
  tl{k} = [];
  tl{k}.label = {'MLC11','MLC12','MLC13','MLC14','MLC15','MLC16','MLC17','MLC21','MLC22','MLC23','MLC24','MLC25','MLC31','MLC32','MLC41','MLC42','MLC51','MLC52','MLC53','MLC54','MLC55','MLC61','MLC62','MLC63','MLF11','MLF12','MLF13','MLF14','MLF21','MLF22','MLF23','MLF24','MLF25','MLF31','MLF32','MLF33','MLF34','MLF35','MLF41','MLF42','MLF43','MLF44','MLF45','MLF46','MLF51','MLF52','MLF53','MLF54','MLF55','MLF56','MLF61','MLF62','MLF63','MLF64','MLF65','MLF66','MLF67','MLO13','MLO14','MLO23','MLO24','MLO31','MLO32','MLO33','MLO34','MLO41','MLO42','MLO43','MLO44','MLO51','MLO52','MLO53','MLP11','MLP12','MLP21','MLP22','MLP23','MLP31','MLP32','MLP33','MLP34','MLP35','MLP41','MLP42','MLP43','MLP44','MLP45','MLP51','MLP52','MLP53','MLP54','MLP55','MLP56','MLP57','MLT11','MLT12','MLT13','MLT14','MLT15','MLT16','MLT21','MLT22','MLT23','MLT24','MLT25','MLT26','MLT27','MLT31','MLT32','MLT33','MLT34','MLT35','MLT36','MLT37','MLT41','MLT42','MLT43','MLT44','MLT45','MLT46','MLT47','MLT51','MLT52','MLT53','MLT54','MLT55','MLT56','MLT57','MRC11','MRC12','MRC13','MRC14','MRC15','MRC16','MRC17','MRC21','MRC22','MRC23','MRC24','MRC25','MRC31','MRC32','MRC41','MRC42','MRC51','MRC52','MRC53','MRC54','MRC55','MRC61','MRC62','MRC63','MRF11','MRF12','MRF13','MRF14','MRF21','MRF22','MRF23','MRF24','MRF25','MRF31','MRF32','MRF33','MRF34','MRF35','MRF41','MRF42','MRF43','MRF44','MRF45','MRF46','MRF51','MRF52','MRF53','MRF54','MRF55','MRF56','MRF61','MRF62','MRF63','MRF64','MRF65','MRF67','MRO11','MRO12','MRO13','MRO14','MRO21','MRO22','MRO23','MRO24','MRO31','MRO32','MRO33','MRO34','MRO41','MRO42','MRO43','MRO44','MRO51','MRO53','MRP11','MRP12','MRP21','MRP22','MRP23','MRP31','MRP32','MRP33','MRP34','MRP35','MRP41','MRP42','MRP43','MRP44','MRP45','MRP51','MRP52','MRP53','MRP54','MRP55','MRP56','MRP57','MRT11','MRT12','MRT13','MRT14','MRT15','MRT16','MRT21','MRT22','MRT23','MRT24','MRT25','MRT26','MRT27','MRT31','MRT32','MRT33','MRT34','MRT35','MRT36','MRT37','MRT41','MRT42','MRT43','MRT44','MRT45','MRT46','MRT47','MRT51','MRT52','MRT53','MRT54','MRT55','MRT56','MRT57','MZC01','MZC02','MZC03','MZC04','MZF01','MZF02','MZF03','MZO01','MZO03','MZP01','MLO11','MLO12','MLO21','MLO22','MZO02'};

  tl{k}.avg = randn(273,100);
  tl{k}.time = 0.01:0.01:1;
  tl{k}.dimord = 'chan_time';

  % shuffle the channel orders for tl_shuff
  perm = randperm(numel(tl{k}.label));
  tl_shuff{k} = tl{k};
  tl_shuff{k}.label = tl{k}.label(perm);
  tl_shuff{k}.avg = tl{k}.avg(perm,:);

  % take the same data as tl for a freq structure with singleton freq
  freq{k} = [];
  freq{k}.label = tl{k}.label;
  freq{k}.freq = 10;
  freq{k}.time = tl{k}.time;
  freq{k}.powspctrm = reshape(tl{k}.avg, [273 1 100]);
  freq{k}.dimord = 'chan_freq_time';

  % and also make a freq_shuff using the same shuffling as for tl_shuff
  freq_shuff{k} = freq{k};
  freq_shuff{k}.label = freq{k}.label(perm);
  freq_shuff{k}.powspctrm = freq{k}.powspctrm(perm,:,:);
end

% note: to visually verify that the shuffling was done correctly simply do
% cfg = [];
% cfg.layout = 'CTF275.lay';
% figure; ft_topoplotER(cfg, tl{1});
% figure; ft_topoplotER(cfg, tl_shuff{1});

%% stats cfg: common stuff for both methods
cfg = [];
cfg.method = 'montecarlo';
cfg.statistic = 'ft_statfun_depsamplesT';
cfg.correctm = 'cluster';
cfg.clusteralpha = 0.05;
cfg.clusterstatistic = 'maxsum';
cfg.minnbchan = 2;

c2 = [];
c2.method = 'template';
c2.template = 'CTF275_neighb.mat';
cfg.neighbours = ft_prepare_neighbours(c2);

cfg.tail = 0;
cfg.clustertail = 0;
cfg.alpha = 0.025;
cfg.numrandomization = 10;
cfg.clusterthreshold = 'nonparametric_common';

cfg.design = [
  ones(1,10) ones(1,10)*2;
  1:10 1:10
];
cfg.ivar  = 1;
cfg.uvar = 2;

%% do the stats

% timelockstatistics
stat_tl = ft_timelockstatistics(cfg, tl{:});
stat_tl_shuff = ft_timelockstatistics(cfg, tl_shuff{:});

% freqstatistics
cfg.avgoverfreq = 'yes';
cfg.frequency = [7 14];
stat_freq = ft_freqstatistics(cfg, freq{:});
stat_freq_shuff = ft_freqstatistics(cfg, freq_shuff{:});

%% diagnostics

% make channel order across structures consistent again
% note that we only fix the .stat field as that is the only one we test
% later on
[sel1,sel2] = match_str(stat_tl.label, stat_tl_shuff.label);
stat_tl_shuff.label = stat_tl_shuff.label(sel2);
stat_tl_shuff.stat = stat_tl_shuff.stat(sel2,:);
stat_tl_shuff.prob = stat_tl_shuff.prob(sel2,:);

[sel1,sel2] = match_str(stat_tl.label, stat_freq.label);
stat_freq.label = stat_freq.label(sel2);
stat_freq.stat = stat_freq.stat(sel2,:,:);
stat_freq.prob = stat_freq.prob(sel2,:,:);

[sel1,sel2] = match_str(stat_tl.label, stat_freq_shuff.label);
stat_freq_shuff.label = stat_freq_shuff.label(sel2);
stat_freq_shuff.stat = stat_freq_shuff.stat(sel2,:,:);
stat_freq_shuff.prob = stat_freq_shuff.prob(sel2,:,:);

fprintf('\ndiagnostic results:\n');
errflag = 0;

if all(stat_tl.stat(:) == stat_tl_shuff.stat(:))
  fprintf('for timelockstatistics, shuffling does not matter\n');
else
  fprintf('ERROR in timelockstatistics when channels are shuffled!\n');
  errflag = 1;
end

if all(stat_freq.stat(:) == stat_freq_shuff.stat(:))
  fprintf('for freqstatistics, shuffling does not matter\n');
else
  fprintf('ERROR in freqstatistics when channels are shuffled!\n');
  errflag = 1;
end

if all(stat_freq.stat(:) == stat_tl.stat(:))
  fprintf('with normal channel order, freqstats equal tlstats\n');
else
  fprintf('with normal channel order, freqstats DOES NOT equal tlstats\n');
  errflag = 1;
end

if all(stat_freq_shuff.stat(:) == stat_tl_shuff.stat(:))
  fprintf('with identically shuffled channel order, freqstats equal tlstats\n');
else
  fprintf('with identically shuffled channel order, freqstats DOES NOT equal tlstats\n');
  errflag = 1;
end

if errflag
  fprintf('\n');
  error('something bad happened, see above for details');
end