text:
%ETS Elementary Transform Sequence class
%
% Manipulate a sequence (vector) of elementary transformations
%  - ETS.TX
%  - ETS.TY
%  - ETS.TZ
%  - ETS.RX
%  - ETS.RY
%  - ETS.RZ
%
% Methods::
%  ETS          Construct a sequence from string
%  isrot        True if rotational transform
%  istrans      True if translational transform
%  isjoint      Is ETS a function of qj
%  njoints      Maximum joint variable index
%  axis         Axis of translation or rotation
%  find         Find ETS that is a function of qj
%  subs         Substitute element of sequence
%-
%  eval         Evaluate ETS
%  jacobian     Compute Jacobian of ETS
%-
%  display      Display a sequence in human readable form
%  char         Convert sequence to a string
%
% Example::
%         ets = ETS('Rx(q1)Tx(a1)Ry(q2)Ty(a3)Rz(q3)Rx(pi/2)')
%         ets.eval([1 2 3]);
%
% Notes::
% - Still experimental
%
% See also trchain, trchain2.

% TODO:
%  - handle 2D case
%  - do DHFactor
%  - accept parameters from a passed struct rather as well as workspace

% 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

classdef ETS < handle

    properties
        type %  ETS.TX, ETS.TY, ETS.TZ, ETS.RX, ETS.RY, ETS.RZ

        val
        joint  % joint number, if a joint, else 0
        constant    % eg. 90, for angles
        sign
        symconstant  % eg. L1, for lengths

        % DH parameters
        theta
        D
        A
        alpha
        offset
        prismatic
    end

    properties (Constant)
        TX = 0;
        TY = 1;
        TZ = 2;
        RX = 3;
        RY = 4;
        RZ = 5;
        DH = 6;
        DHM = 7;

        names = {'Tx', 'Ty', 'Tz', 'Rx', 'Ry', 'Rz', 'DH', 'DHm' };
    end
    % Element.java
    %         // one of ETS.TX, ETS.TY ... ETS.RZ, DH_STANDARD/MODIFIED
    % 	int		type;
    %
    % 	// transform parameters, only one of these is set
    % 	String	var;        // eg. q1, for joint var ETS.TYpes
    % 	String	symconst;   // eg. L1, for lengths
    % 	int		constant;   // eg. 90, for angles
    %
    % 	// DH parameters, only set if ETS.TYpe is DH_STANDARD/MODIFIED
    % 	int 	theta,
    %             alpha;
    %     String  A,
    %             D;
    %     int     prismatic;
    %     int     offset;

    methods
        function ets = ETS(s, varargin)
            %ETS.ETS Construct elementary transform element or sequence
            %
            % e = ETS() is a new ETS object.
            %
            % e = ETS(t) is a clone of the ETS object t and all properties are copied.
            %
            % e = ETS(op, v) is a new ETS object of type op and value v.  OP can be any
            % of
            % 'Rx'       rotation about the x-axis
            % 'Ry'       rotation about the y-axis
            % 'Rz'       rotation about the z-axis
            % 'Tx'       translation along the x-axis
            % 'Ty'       translation along the y-axis
            % 'Tz'       translation along the z-axis
            % 'transl'   sequence of finite translations along the x-, y- and z-directions.
            % 'rpy'      sequence of finite rotations about the x-, y- and z-directions.
            %
            % e = ETS(str) is a sequence of ETS objects, each described by a
            % subexpression in the string STR.  Each subexpression comprises an
            % operation as per the table above followed by parentheses and a value.
            % For example:
            %
            %         ets = ETS('Rx(q1)Tx(a1)Ry(q2)Ty(a3)Rz(q3)Rx(pi/2)')

            if nargin == 0
                ets.joint = false;
                ets.constant = [];
                return;
            end

            if ischar(s) || isstring(s)
                s = string(s);
                if strfind(s, '(')
                    % string of tokens
                    tokens = regexp(s, '\s*(?<op>R.?|T.)\(\s*(?<arg>[^)]*)\s*\)\s*', 'names');

                    joint = 1;

                    ets = [];
                    for token = tokens

                        % get the argument for this transform element

                        % deal with case of symbolic arg or workspace expression
                        if token.arg.startsWith("q") || token.arg.startsWith("*")
                            x = ETS(token.op, "");
                            x.joint = joint;
                            joint = joint+1;
                        else
                            x = ETS(token.op, string(token.arg));
                            x.joint = 0;
                        end

                        ets = [ets x];  % append to the list
                    end
                else
                    % ETS('transl', [v])
                    % ETS('rpy', [v]);
                    switch s
                        case 'transl'
                            v = varargin{1};

                            ets = [];
                            if v(1) ~= 0
                                ets = [ets ETS('Tx', v(1))];
                            end
                            if v(2) ~= 0
                                ets = [ets ETS('Ty', v(2))];
                            end
                            if v(3) ~= 0
                                ets = [ets ETS('Tz', v(3))];
                            end
                        case 'rpy'
                            v = varargin{1};

                            ets = [];
                            if v(1) ~= 0
                                ets = [ets ETS('Rx', v(1))];
                            end
                            if v(2) ~= 0
                                ets = [ets ETS('Ry', v(2))];
                            end
                            if v(3) ~= 0
                                ets = [ets ETS('Rz', v(3))];
                            end

                        case 'Rx'
                            ets = ETS(ETS.RX, varargin{1});
                        case 'Ry'
                            ets = ETS(ETS.RY, varargin{1});
                        case 'Rz'
                            ets = ETS(ETS.RZ, varargin{1});
                        case 'Tx'
                            ets = ETS(ETS.TX, varargin{1});
                        case 'Ty'
                            ets = ETS(ETS.TY, varargin{1});
                        case 'Tz'
                            ets = ETS(ETS.TZ, varargin{1});

                        otherwise
                            error('RTB:trchain:badarg', 'unknown operator/option %s', s);
                    end
                end

            elseif isnumeric(s)
                ets = ETS();
                ets.type = s;
                ets.val = varargin{1};

            elseif isa(s, 'ETS')
                % clone the object
                ets.type = s.type;
                ets.joint = s.joint;
                ets.val = s.val;
                ets.symconstant = s.symconstant;

                if nargin > 1
                    ets.type = varargin{1};
                end

                if nargin > 2
                    if varargin{2} < 0
                        ets = ets.negate();
                    end
                end
            end
        end

        function s = op(ets)
            s = ets.names{ets.type+1};
        end

        function s = char(ets)
            s = '';
            for x=ets
                e =  x.op;

                if isjoint(x)
                    e = " " + e + "(*q" +  x.joint + x.val + ")";
                else
                    e = e + "(" + x.val + ")";
                end
                s = strcat(s, e);
            end
        end

        function display(ets)
            %ETS.display Display parameters
            %
            % ETS.display() displays the transform parameters in compact single line format.
            %
            % Notes::
            % - This method is invoked implicitly at the command line when the result
            %   of an expression is a Link object and the command has no trailing
            %   semicolon.
            %
            % See also Link.char, Link.dyn, SerialLink.showlink.
            loose = strcmp( get(0, 'FormatSpacing'), 'loose');
            if loose
                disp(' ');
            end
            disp([inputname(1), ' = '])
            disp( char(ets) );
        end % display()

%         function disp(ets)
%             disp( show(ets) );
%         end

        %-------------------------------------------------------------

        function T = eval(ets, q)

            T = eye(4,4);

            for x = ets
                % get the argument for this transform element
                if x.isjoint()
                    % from the passed in vector q

                    arg = q(x.joint);
                else            % or the workspace

                    try
                        arg = evalin('base', x.val);
                    catch
                        error('RTB:ETS:badarg', 'variable %s does not exist', x.val);
                    end
                end

                % now evaluate the element and update the transform chain
                switch x.type
                    case ETS.RX
                        T = T * trotx(arg, 'deg');
                    case ETS.RY
                        T = T * troty(arg, 'deg');
                    case ETS.RZ
                        T = T * trotz(arg, 'deg');
                    case ETS.TX
                        T = T * transl(arg, 0, 0);
                    case ETS.TY
                        T = T * transl(0, arg, 0);
                    case ETS.TZ
                        T = T * transl(0, 0, arg);
                    otherwise
                        error('RTB:ETS:badarg', 'unknown operator %s', x.op);
                end
            end

            T = trnorm(T);

            %             if isa(q, 'symfun')
            %                 T = formula(T);
            %             end
        end

        function J = jacobian(ets, q)
            n = ets.njoints();
            J = zeros(6, n);
            for j=1:n
                % find derivative with respect to q_j
                k = ets.find(j);
                deriv = ets(k).deriv(q(j));
                pre = eval( ets(1:k-1), q );
                post = eval( ets(k+1:end), q );
                Td = trnorm(pre)*deriv*trnorm(post);
                J(:,j) = [ Td(1:3,4); vex(t2r(Td)) ];
            end
        end

        function T = deriv(ets, q)
            switch ets.type
                case ETS.RX
                    T = r2t( skew([1 0 0])*rotx(q) );
                case ETS.RY
                    T = r2t( skew([0 1 0])*roty(q) );
                case ETS.RZ
                    T = r2t( skew([0 0 1])*rotz(q) );
                case ETS.TX
                    T = transl(1, 0, 0);
                case ETS.TY
                    T = transl(0, 1, 0);
                case ETS.TZ
                    T = transl(0, 0, 1);
            end
            T(4,4) = 0;
        end

%********************************************************************
%     D H F A C T O R
%********************************************************************

        function xo = subs(ets, k, x)
            fprintf('  Subs: %s := %s\n', char(ets(k)), char(x));
            xo = [ets(1:k-1) x ets(k+1:end)];
        end

        function k = find(ets, j)
            [~,k] = find([ets.joint] == j);
        end

        function n = njoints(ets)
            n = max([ets.joint]');
        end

        function b = isjoint(ets, j)
            b = [ets.joint] ~= 0;
        end

        function b = istrans(ets)
            b = ismember(ets.type, [ETS.TX ETS.TY ETS.TZ]);

        end

        function b = isrot(ets)
            b = ismember(ets.type, [ETS.RX ETS.RY ETS.RZ]);
        end

        function b = axis(ets)
            b = char( mod([ets.type],3) + 'x' );

        end

        function v = sametype(ets, next)
            v = ets.type == next.type;
        end

        function s = show(ets)
            for i=1:numel(ets)
                e = ets(i);
                s = sprintf('%s: type=%d, joint %d, val=%g. const=%s. sign=%d, sym=%s\n', ...
                    char(e), e.type, e.joint, e.val, e.constant, e.sign, e.symconstant);
                if nargout == 0
                    disp(s);
                end
            end
        end

                function negate(ets)
            ets.constant = -ets.constant;

            s = ets.symconst;

            % add initial sign char if none
            if s(1) ~= '+' && s(1) ~= '-'
                s = ['+' s];
            end

            % go through the string and flip all sign chars
            kp = strfind(s, '+');
            kn = strfind(s, '-');
            s(kp) = '-';
            s(kn) = '+';

            % if inital sign is + remove it
            if s(1) == '+';
                s = s(2:end);
            end

            ets.symconst = s;
        end

        % methods from Element.java
        %  *	public boolean istrans()
        %  *	public boolean isrot()
        %  *	public int axis()
        %  *	public boolean isjoint() {
        %  *	public boolean factorMatch(int dhWhich, int i, int verbose) {
        %  *	public void add(Element e) {
        %  *	public Element(int ETS.TYpe, int constant) {
        %  *	public Element(int ETS.TYpe)	// new of specified ETS.TYpe
        %  *
        %  *	public static String toString(Element [] e) {
        %  *	public String argString() {
        %  *	public String toString() {
        %  *
        %  * Constructors:
        %  *	Element(Element e) 	// clone of argument
        %  *	Element(Element e, int ETS.TYpe, int sign) // clone of argument with new ETS.TYpe
        %  *	Element(Element e, int ETS.TYpe) // clone of argument with new ETS.TYpe
        %  *	Element(String s)

        function v = plus(ets, e)
            if isstring(ets) && isa(e, 'ETS')
                % string concatenate
                v = ets + char(e);
            end
        end

        function symPlus(ets, e)
        end

        function new = merge(a, b)
            if isjoint(a)
                new = ETS(a);
                new.val = num2str(eval(new.val + "+" + b.val));
            elseif isjoint(b)
                new = ETS(b);
                new.val = num2str(eval(new.val + "+" + a.val));
            else
                % two constants
                new = ETS(a);
                new.val = num2str(eval(new.val + "+" + b.val));
                new.symconstant = new.symconstant + "+" +  b.symconstant;
                if eval(new.val) == 0
                    new = [];
                end
            end
            if isempty(new)
                fprintf('  Merging: %s %s -> (nil)\n', char(a), char(b));
            else
                fprintf('  Merging: %s %s -> %s\n', char(a), char(b), char(new));
            end
        end

        function b = swap(ets, next, dhWhich)

            b = false

            % don't swap if both are joint variables

            if ets.isjoint() && next.isjoint()
                return
            end

            switch (dhWhich)
                case 'standard'
                    % we want to sort terms into the order:	RZ	TX	TZ	RX

                    if  ((ets.type == ETS.TZ) && (next.type == ETS.TX)) || ...
                            ((ets.type == ETS.TX) && (next.type == ETS.RX) && next.isjoint()) || ...  % push constant translations through rotational joints of the same ETS.TYpe
                            ((ets.type == ETS.TY) && (next.type == ETS.RY)) && next.isjoint() || ...
                            ((ets.type == ETS.TZ) && (next.type == ETS.RZ)) && next.isjoint() || ...
                            (~ets.isjoint() && (this.type == ETS.RX) && (next.type == ETS.TX)) || ...
                            (~ets.isjoint() && (this.type == ETS.RY) && (next.type == ETS.TY)) ||....
                            (~ets.isjoint() && ~next.isjoint() && (this.type == ETS.TZ) && (next.type == ETS.RZ)) || ...
                            ((ets.type == ETS.TY) && (next.type == ETS.TZ)) || ...  % move ETS.TY terms to the right
                            ((ets.type == ETS.TY) && (next.type == ETS.TX))

                        fprintf(['Swap: ' char(ets) ' <-> ' char(next)] );
                        b = true;
                    end
                case 'modified'
                    if  ((ets.type == ETS.RX) && (next.type == ETS.TX)) || ...
                            ((ets.type == ETS.RY) && (next.type == ETS.TY)) || ...
                            ((ets.type == ETS.RZ) && (next.type == ETS.TZ)) || ...
                            ((ets.type == ETS.TZ) && (next.type == ETS.TX))

                        fprintf(['Swap: ' char(ets) ' <-> ' char(next)] );
                        b = true;
                    end
                otherwise
                    error('bad DH ETS.TYpe');
            end
            b = false;
        end

        %         	/**
        % 	 * Substitute this transform for a triple of transforms
        %      * that includes an ETS.RZ or ETS.TZ.
        %      *
        % 	 * @return	- null if no substituion required
        % 	 *			- array of Elements to substitute
        % 	 */
        function s = substituteToZ(ets, prev)

            switch ets.type

                case ETS.TX
                    s(1) = ETS(ETS.RY, 90);
                    s(2) = ETS(ets, ETS.TZ);
                    s(3) = ETS(ETS.RY, -90);

                case ETS.TY
                    % there are two options here
                    if prev.type == ETS.RZ
                        s(1) = ETS(ETS.RZ, 90);
                        s(2) = ETS(ets, ETS.TX);
                        s(3) = ETS(ETS.RZ, -90);
                    else
                        s(1) = ETS(ETS.RX, -90);
                        s(2) = ETS(ets, ETS.TZ);
                        s(3) = ETS(ETS.RX, 90);
                    end

                case ETS.RX
                    s(1) = ETS(ETS.RY, 90);
                    s(2) = ETS(ets, ETS.RZ);
                    s(3) = ETS(ETS.RY, -90);
                case ETS.RY
                    s(1) = ETS(ETS.RX, -90);
                    s(2) = ETS(ets, ETS.RZ);
                    s(3) = ETS(ETS.RX, 90);

                case ETS.RY
                    s(1) = ETS(ETS.RZ, 90);
                    s(2) = ETS(ets, ETS.RX);
                    s(3) = ETS(ETS.RZ, -90);
                otherwise
                    s = [];
            end
        end

        function s = substituteTY(ets, k)

            pref = -1;

            try
                if ets(k+1).axis == 'x'
                    pref = ETS.RX;
                end
            end

            try
                if ets(k+1).axis == 'z'
                    pref = ETS.RZ;
                end
            end

            assert(pref >= 0, 'substituteTY, can''t decide');

            if pref == ETS.RX
                s(1) = ETS(ETS.RX, -90);
                s(2) = ETS(ets(k), ETS.TZ);
                s(3) = ETS(ETS.RX, 90);
            else
                s(1) = ETS(ETS.RZ, 90);
                s(2) = ETS(ets(k), ETS.TX);
                s(3) = ETS(ETS.RZ, -90);
            end

        end

        function s = substituteRY(ets)

                s(1) = ETS(ETS.RZ, 90);
                s(2) = ETS(ets, ETS.RX);
                s(3) = ETS(ETS.RZ, -90);
        end

        % 	/**
        % 	 * Simple rewriting rule for adjacent transform pairs.  Attempt to
        % 	 * eliminate ETS.TY and ETS.RY.
        % 	 * @param	previous element in list
        % 	 * @return	- null if no substituion required
        % 	 *			- array of Elements to subsitute
        % 	 */
        function s = substituteY(this, prev, next)

            s = [];

            if (prev.isjoint() || e.isjoint())
                return
            end

            % note that if rotation is -90 we must make the displacement -ve */
            if ((prev.type == ETS.RX) && (this.type == ETS.TY))
                % RX.TY -> TZ.RX
                s(1) = ETS(this, ETS.TZ, prev.constant);
                s(2) = ETS(prev);
            elseif ((prev.type == ETS.RX) && (this.type == ETS.TZ))
                % RX.TZ -> TY.RX
                s(1) = ETS(this, ETS.TY, -prev.constant);
                s(2) = ETS(prev);
            elseif ((prev.type == ETS.RY) && (this.type == ETS.TX))
                % RY.TX-> TZ.RY
                s(1) = ETS(this, ETS.TZ, -prev.constant);
                s(2) = ETS(prev);
            elseif ((prev.type == ETS.RY) && (this.type == ETS.TZ))
                % RY.TZ-> TX.RY
                s(1) = ETS(this, ETS.TX, prev.constant);
                s(2) = ETS(prev);
            elseif ((prev.type == ETS.TY) && (this.type == ETS.RX))
                % TY.RX -> RX.TZ
                s(1) = ETS(this);
                s(2) = ETS(prev, ETS.TZ, -this.constant);
                %%return s;
                s = [];
            elseif ((prev.type == ETS.RY) && (this.type == ETS.RX))
                % RY(Q).RX -> RX.RZ(-Q)
                s(1) = ETS(this);
                s(2) = ETS(prev, ETS.RZ, -1);
            elseif ((prev.type == ETS.RX) && (this.type == ETS.RY))
                % RX.RY -> RZ.RX
                s(1) = ETS(this, ETS.RZ);
                s(2) = ETS(prev);
            elseif ((prev.type == ETS.RZ) && (this.type == ETS.RX))
                % RZ.RX -> RX.RY
                s(1) = ETS(this);
                s(2) = ETS(prev, ETS.RY);
                s = [];
            end
        end

        function out = mergeterms(ets)
            i = 1;
            while i < numel(ets)-1
                cur = ets(i);
                next = ets(i+1);
                if cur.sametype(next)
                    new = merge(cur, next);
                    if isempty(new)
                        % cancelation
                        ets(i:i+1) = [];
                    else
                        ets(i) = new;
                        ets(i+1) = [];
                    end
                else
                    i = i+ 1;
                end

            end
            out = ets;
        end

        function out = ordering(ets)
            for i=1:numel(ets)-1
                cur = ets(i);
                next = ets(i+1);

                % TZ RZ -> RZ TZ
                if cur.type == ETS.TZ && next.type == ETS.RZ && next.isjoint()
                    ets(i) = next;
                    ets(i+1) = cur;
                    fprintf('  Swap:  %s <--> %s\n', char(cur), char(next));
                end

                % RX TX -> TX RX
                if cur.type == ETS.RX && next.type == ETS.TX && ~next.isjoint()
                    ets(i) = next;
                    ets(i+1) = cur;
                    fprintf('  Swap:  %s <--> %s\n', char(cur), char(next));
                end
            end
            out = ets;
        end

        function out = simplify(ets)
            fprintf('--------- simplify\n');
            ets = mergeterms(ets);
            while 1
                len = numel(ets);

                ets = ordering(ets);
                ets = mergeterms(ets);

                if numel(ets) == len
                    break;
                end
            end

            out = ets;
        end

        function v = contains(ets, which)
            w = ETS.name2type(which);
            for e = ets
                v = e.type == w;
                if v
                    return
                end
            end
        end

        function dhfactor(ets)
            disp(char(ets))

            % find which joint to substitute
            fprintf('--------- align joints with z-axis\n');
            fixup = find( ets.isjoint & ([ets.type] ~= ETS.RZ) );

            for k=fliplr(fixup)   % work right to left so indices are not broken
                new = substituteToZ(ets(k));
                ets = subs(ets, k, new);
            end
            disp(char(ets))

            %
            ets = simplify(ets);
            disp(char(ets))

            % now eliminate all Ry
            if ets.contains('Ry')
                fprintf('--------- eliminate RY\n');

                fixup = find( [ets.type] == ETS.RY );

                for k=fliplr(fixup)   % work right to left so indices are not broken
                    new = substituteRY(ets(k));
                    ets = subs(ets, k, new);
                end

                disp(char(ets))

                ets = simplify(ets);
            end

            % now eliminate all Ty
            if ets.contains('Ty')
                fprintf('--------- eliminate TY\n');
                fixup = find( [ets.type] == ETS.TY );

                for k=fliplr(fixup)   % work right to left so indices are not broken
                    new = substituteTY(ets, k);
                    ets = subs(ets, k, new);
                end

                disp(char(ets))

                ets = simplify(ets);
            end

            ets = ordering(ets);

            disp(char(ets))

            factorize(ets)
        end

        function nfactors = factorize(ets, mdh)
            if nargin < 2
                mdh = false;
            end
            verbose = true;
            nfactors = 0;

            out = [];

            for i=1:length(ets)
                j = i;
                jvars = 0; match = 0;

                % scan next 4 terms
                if verbose
                    fprintf("start at " + ets(i));
                end

                match = 0; jvars = 0; j = i;
                for f=1:4
                    if j > length(ets)
                        break;
                    else
                        e = ets(j);
                    end

                    if e.factormatch(f, mdh)
                        j = j + 1;
                        match = match + 1;
                        if (e.isjoint())
                            jvars = jvars + 1;
                        end
                        if jvars > 1	% can only have 1 joint var per DH
                            break;
                        end
                    end
                end

                if match == 0 || jvars == 0
                    continue;		% no DH subexpression found, keep looking
                end

                if verbose
                    fprintf(" found subexpression " + match + " " + jvars + "\n");
                end

                if jvars == 0
                    continue;
                end

                first = i;
                last = j;
                if jvars > 1
                    last = last - 1;
                end

                [first last]

%                 dh = ETS();
%                 dh.type = ETS.DH;
%
%                 for j=first:last
%                     dh.add( this.get(i) );
%                     this.remove(i);
%                 end
%
%                 this.add(i, dh);
%                 nfactors = nfactors + 1;
%                 if verbose
%                     fprintf(" result: " + dh);
%                 end

            end
        end

        function match = factormatch(ets, i, mdh)

            if nargin < 2
                mdh = false;
            end

            if mdh
                dhFactors = [
                    ETS.RX 0
                    ETS.TX 0
                    ETS.RZ 1
                    ETS.RZ 1 ];
            else
                dhFactors = [
                    ETS.RZ 1
                    ETS.TX 0
                    ETS.TZ 1
                    ETS.RX 0 ];
            end

            match =	(ets.type == dhFactors(i,1)) && ...
            ~((dhFactors(i,2) == 0) && ets.isjoint());

            fprintf(" matching " + ets + " (i=" + i + ") " + " to " + ...
                op(dhFactors(i,1)) + "<" + dhFactors(i,2) + ">" + ...
                " -> " + match + "\n");
        end

        function new = 	add(e1, e2)
            assert(e1.type == DH || e1.type == DHM, 'wrong element type');

            fprintf("  adding: " + e1 + " += "  + e2);
            switch e2.type
                case RZ
                    if e2.isjoint()
                        e1.prismatic = 0;
                        e1.var = e1.var;
                        e1.offset = e1.constant;
                        e1.theta = 0.0;
                    else
                        e1.theta = e1.constant;
                    end
                case TX
                    e1.A = e2.symconst;
                case TZ
                    if e1.isjoint()
                        e1.prismatic = 1;
                        e1.var = e2.var;
                        e1.D = null;
                    else
                        e1.D = e2.symconst;
                    end
                case RX
                    e1.alpha = e2.constant;
                otherwise

                    error("cant factorize " + e);
            end
        end
    end

    methods (Static)
        function t = name2type(name)
            switch lower(name)
                case 'rx', t = ETS.RX;
                case 'ry', t = ETS.RY;
                case 'rz', t = ETS.RZ;
                case 'tx', t = ETS.TX;
                case 'ty', t = ETS.TY;
                case 'tz', t = ETS.TZ;
            end
        end
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function P = jacobiPD(n,a,b,x)
    P = 0;
    for  s = 0:n
        P = P + ...
            1/(factorial(s)*factorial(n+a-s)*factorial(b+s)*factorial(n-s)) * ...
            ((x-1)/2).^(n-s).*((x+1)/2).^s;
    end
    P = P*factorial(n+a) * factorial(n+b);
end

author: goGPS-Project


text:
function pageview(mode,cur_axs)
% PAGEVIEW View menu.
%          PAGEVIEW handles all the axis limit altering options under the
%          View pulldown menu.

% Author: Craig Borghesani
% Date: 8/10/94
% Revised: 10/18/94
% Copyright (c) 1999, Prentice-Hall

% obtain handle information
f = gcf;
ui_data = get(f,'userdata');
ui_han = ui_data{1};
last_lims = ui_han(36);

stat_bar = get(ui_han(43),'userdata');

if nargin == 1,
 cur_axs = get(ui_han(32),'userdata');
end

axs_data = get(cur_axs,'userdata');
axs_typ = axs_data(1);

if mode == 1, % Full

 if nargin == 1, set(stat_bar,'string','Setting axis limits to Full'); end

 cur_xlim = get(cur_axs,'xlim');
 cur_ylim = get(cur_axs,'ylim');
 set(last_lims,'userdata',[cur_axs,cur_xlim,cur_ylim]);

 lins = get(cur_axs,'children');
 cur_xlabel = get(cur_axs,'xlabel');
 cur_ylabel = get(cur_axs,'ylabel');
 cur_title = get(cur_axs,'title');
 lims = [1/eps,-1/eps,1/eps,-1/eps];
 for k = 1:length(lins),
  if strcmp(get(lins(k),'vis'),'on'),
   if strcmp(get(lins(k),'type'),'line'),
    xdata = get(lins(k),'xdata');  ydata = get(lins(k),'ydata');
   elseif ~any(lins(k) == [cur_xlabel,cur_ylabel,cur_title]),
    locdata = get(lins(k),'pos');
    xdata = locdata(1); ydata = locdata(2);
   else
    xdata = []; ydata = [];
   end
   if length(xdata) & (abs(xdata(1))<1000 & abs(ydata(1))<1000 & axs_typ==5) | axs_typ~=5,
    lims(1) = min([xdata(~isnan(xdata)),lims(1)]);
    lims(2) = max([xdata(~isnan(xdata)),lims(2)]);
    lims(3) = min([ydata(~isnan(ydata)),lims(3)]);
    lims(4) = max([ydata(~isnan(ydata)),lims(4)]);
   end
  end
 end

 if (lims(1) == 1/eps) | (lims(2) == -1/eps) | all(lims(1:2)==0),
  if any(axs_typ == [1,5]), lims(1) = -1; lims(2) = 1;
  elseif any(axs_typ==[2,3,6,7]), lims(1) = 0.01; lims(2) = 1000;
  elseif axs_typ == 4, lims(1) = -360; lims(2) = 0;
  elseif axs_typ == 8, lims(1) = 0; lims(2) = 10; end
 end

 if (lims(3) == 1/eps) | (lims(4) == -1/eps) | all(lims(3:4)==0),
  if any(axs_typ == [1,5]), lims(3) = -1; lims(4) = 1;
  elseif any(axs_typ == [2,4]), lims(3) = -20; lims(4) = 20;
  elseif axs_typ == 3, lims(3) = -360; lims(4) = 0;
  elseif axs_typ == 6, lims(3) = 0.1; lims(4) = 10;
  elseif axs_typ == 7, lims(3) = 0; lims(4) = 360;
  elseif axs_typ == 8, lims(3) = 0; lims(4) = 5; end
 end

 if any(axs_typ==[1,4,5]),
  if diff(lims(1:2))==0, lims(1:2) = [-abs(lims(1)),abs(lims(1))]; end
  if diff(lims(3:4))==0, lims(3:4) = [-abs(lims(3)),abs(lims(3))]; end
  lims(1:2) = lims(1:2) + diff(lims(1:2))*[-0.1,0.1];
  lims(3:4) = lims(3:4) + diff(lims(3:4))*[-0.1,0.1];
 elseif any(axs_typ == [2,3,7,8]),
  if diff(lims(3:4))==0, lims(3:4) = [lims(3)-10,lims(3)+10]; end
  lims(3:4) = lims(3:4) + diff(lims(3:4))*[-0.1,0.1];
 else
  loglims = log10(lims(3:4));
  loglims = loglims + diff(loglims)*[-0.1,0.1];
  lims(3:4) = 10 .^(loglims);
 end

 if axs_typ == 3,
  if lims(3) < -360, lims(3) = -360; end
  if lims(4) > 0, lims(4) = 0; end
 elseif axs_typ == 4,
  if lims(1) < -360, lims(1) = -360; end
  if lims(2) > 0, lims(2) = 0; end
 elseif axs_typ == 5,
  if lims(1) > 0, lims(1) = -1; end
  if lims(2) < 0, lims(2) = 1; end
 end

 if any(axs_typ == [1, 5]),
  set(cur_axs,{'xlim','ylim'},ctsqaxes(cur_axs,lims(1:2),lims(3:4)));
 else
  set(cur_axs,'xlim',lims(1:2),'ylim',lims(3:4));
 end

elseif any(mode == [2,3]), % In/Out

 if mode == 2,
  out = 0.15;
  if nargin == 1, set(stat_bar,'string','Zooming In'); end
 else
  out = -0.15;
  if nargin == 1, set(stat_bar,'string','Zooming Out'); end
 end

 cur_xlim = get(cur_axs,'xlim');
 cur_ylim = get(cur_axs,'ylim');
 set(last_lims,'userdata',[cur_axs,cur_xlim,cur_ylim]);

 if any(axs_typ == [1,4,5,8]),
  delx = diff(cur_xlim);
  new_xlim = cur_xlim+out*[delx,-delx];
  dely = diff(cur_ylim);
  new_ylim = cur_ylim+out*[dely,-dely];
 elseif any(axs_typ == [2,3,7]),
  loglims = log10(cur_xlim);
  delx = diff(loglims);
  loglims = loglims + out*[delx,-delx];
  new_xlim = 10 .^(loglims);
  dely = diff(cur_ylim);
  new_ylim = cur_ylim+out*[dely,-dely];
 else
  loglims = log10(cur_xlim);
  delx = diff(loglims);
  loglims = loglims + out*[delx,-delx];
  new_xlim = 10 .^(loglims);
  loglims = log10(cur_ylim);
  dely = diff(loglims);
  loglims = loglims + out*[dely,-dely];
  new_ylim = 10 .^(loglims);
 end

 if axs_typ == 3,
  if new_ylim(1) < -360, new_ylim(1) = -360; end
  if new_ylim(2) > 0, new_ylim(2) = 0; end
 elseif axs_typ == 4,
  if new_xlim(1) < -360, new_xlim(1) = -360; end
  if new_xlim(2) > 0, new_xlim(2) = 0; end
 end

 if any(axs_typ == [1, 5]),
  set(cur_axs,{'xlim','ylim'},ctsqaxes(cur_axs,new_xlim,new_ylim));
 else
  set(cur_axs,'xlim',new_xlim,'ylim',new_ylim);
 end

elseif mode == 4, % Zoom

 set(stat_bar,'string','Zoom Mode: press and drag mouse to define new limits');
 set(f,'pointer','crosshair');
 set(f,'windowbuttondownfcn','pageview(10)','windowbuttonupfcn','1;');

elseif mode == 10, % Zoom (RBBOX stage)

 cur_obj = get(f,'currentobject');
 par_cur = get(cur_obj,'parent');
 sel_axs = 0;
 if par_cur == f, % if the parent is the figure, we have an axis object
  sel_axs = cur_obj;
 elseif par_cur,  % if the parent is not the root, then it must be an axis
  sel_axs = par_cur;
 end

 if sel_axs,

  set(stat_bar,'string','Selecting region...');
  pt0 = get(sel_axs,'currentpoint');
  rbbox([get(f,'currentpoint'),0,0],get(f,'currentpoint'));
  drawnow;

  set(stat_bar,'string','Setting new axis limits');
  pt1 = get(sel_axs,'currentpoint');

  cur_xlim = get(sel_axs,'xlim');
  cur_ylim = get(sel_axs,'ylim');
  set(last_lims,'userdata',[sel_axs,cur_xlim,cur_ylim]);

  new_xlim = [min(pt0(1,1),pt1(1,1)), max(pt0(1,1),pt1(1,1))];
  new_ylim = [min(pt0(1,2),pt1(1,2)), max(pt0(1,2),pt1(1,2))];

  if axs_typ == 3,
   if new_ylim(1) < -360, new_ylim(1) = -360; end
   if new_ylim(2) > 0, new_ylim(2) = 0; end
  elseif axs_typ == 4,
   if new_xlim(1) < -360, new_xlim(1) = -360; end
   if new_xlim(2) > 0, new_xlim(2) = 0; end
  end

  set(f,'windowbuttondownfcn','pagemous','windowbuttonupfcn','',...
        'pointer','arrow');

  axs_data = get(sel_axs,'userdata');
  axs_typ = axs_data(1);
  if any(axs_typ == [1, 5]),
   set(sel_axs,{'xlim','ylim'},ctsqaxes(sel_axs,new_xlim,new_ylim));
  else
   set(sel_axs,'xlim',new_xlim,'ylim',new_ylim);
  end

 else

  set(stat_bar,'string','Need to be in an axis region to use Zoom');

 end

elseif mode == 5, % Last

 data = get(last_lims,'userdata');
 set(stat_bar,'string','Restoring Last axis limits');

 cur_xlim = get(cur_axs,'xlim');
 cur_ylim = get(cur_axs,'ylim');
 set(last_lims,'userdata',[cur_axs,cur_xlim,cur_ylim]);

 axs_data = get(data(1),'userdata');
 axs_typ = axs_data(1);
 if any(axs_typ == [1, 5]),
  set(data(1),{'xlim','ylim'},ctsqaxes(cur_axs,cur_xlim,cur_ylim));
 else
  set(data(1),'xlim',cur_xlim,'ylim',cur_ylim);
 end

elseif mode == 6, % grid on/off

 if strcmp(get(cur_axs,'xgrid'),'on'),
  set(stat_bar,'string','Turning grid Off');
  set(cur_axs,'xgrid','off','ygrid','off');
 else
  set(stat_bar,'string','Turning grid On');
  set(cur_axs,'xgrid','on','ygrid','on');
 end

end

author: Sable

repo: mcbench-benchmarks


text:
function zPoints=getEllipsePoints(z,A,gammaVal,numPoints,invertA)
%%GETELLIPSEPOINTS Generate numPoints around a 2D ellipse. These can be
%           used, for example, to plot the ellipse. A point zp on the
%           ellipse/ ellipsoid satisfies the equation
%           (zp-z)'*A*(zp-z)=gammaVal
%           or is invertA is true, the equation
%           (zp-z)'*inv(A)*(zp-z)=gammaVal
%           gammaVal determines what amount of probability is in the
%           region. If omitted, gammaVal=16.2235, which corresponds
%           approximately to the 99.97% confidence region, is used (The
%           value comes from ChiSquareD.invCDF(0.9997,2)).
%
%INPUTS: z A 2XN vector corresponding to the centers of the N ellipses for
%          which points should be obtained.
%        A A 2X2XN set of N positive definite matrices that specify the
%          size and shape of the ellipses, where a point zp is on the ith
%          ellipse if
%          (zp-z(:,i))'*A(:,:,i)*(zp-z(:,i))=gammaVal (if invertA is true,
%          then replace A(:,:,i) with inv(A(:,:,i))).
% gammaVal An optional parameter specifying the size of the ellipse/
%          ellipsoid. If omitted or an empty matrix is passed, then
%          gammaVal=16.2235 is used. This is approximately the value for a
%          99.97% confidence region if A are inverse covariance matrices of
%          a Gaussian distribution. gammaVal must be positive.
% numPoints An optional parameter specifying how many points should be
%          generated. The default if omitted or an empty matrix is passed
%          is 2000. If this number is not an even number, then the next
%          largest even number of points is generated.
%  invertA If this is true, then A is inverted before use. The default if
%          omitted or an empty matrix is passed is false.
%
%OUTPUTS: zPoints A 2XnumPointsXN set of the points for each of the
%                 ellipses.
%
%An eigendecomposition of a matrix A breaks it into parts V and D such that
%V*D*V'=A, where V is a rotation matrix and D is a diagonal matrix.
%Considering the problem here, this means that (zp-z)'*A*(zp-z)=
%(V'*zp-V'*z)'*D*(V'*zp-V'*z). If we substitute kp=V'*zp and k=V'*z, then
%the equation is (kp-k)'*D*(kp-k)=gammaVal, which is the equation for an
%ellipse where the axes are aligned with the coordinate axes. Thus, one can
%first find the points for an ellipse that is aligned with the coordinate
%axes (the kp points), and then rotate it back to the proper alignment (the
%zp points). However, it is simpler to first find an ellipse centered about
%zero. Thus, substitute l=kp-k and find the l points, then shift and rotate
%them to get the z-points.
%
%When considering an ellipse, suppose that l=[x;y] and D=diag(a,b). Then
%the ellipse equation is
%x^2*a+y^2*b=gammaVal
%Solving for y in terms of x, we have
%y=+/-sqrt((gammaVal-x^2*a)/b).
%The x values are limited so that the argument of the square root is
%positive. Thus the x values range from -sqrt(gammaVal/a) to
%+sqrt(gammaVal/a). Thus, a simple way to plot this centered, axis-aligned
%ellipse is to generate x-values in the valid range and then find the two
%sets of y values. This function does that and transforms the results to
%get the final result.
%
%EXAMPLE:
%Here, we generate points for two ellipses and then plot them.
% z=[[4;3],[3;4]];
% A=zeros(2,2,2);
% A(:,:,1)=[4,-2;
%           -2,5];
% A(:,:,2)=[5,-2;
%           -2,5];
% zPoints=getEllipsePoints(z,A,[],[],true);
% figure(1)
% clf
% hold on
% plot(zPoints(1,:,1),zPoints(2,:,1),'-r')
% plot(zPoints(1,:,2),zPoints(2,:,2),'-b')
%
%December 2020 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<5||isempty(invertA))
    invertA=false;
end

if(nargin<4||isempty(numPoints))
    numPoints=2000;
end

if(nargin<3||isempty(gammaVal))
    gammaVal=16.2235;
end

if(invertA)
    A=applyFunToEachMatrix(@inv,A);
end

N=size(z,2);

numPoints=ceil(numPoints/2);

zPoints=zeros(2,2*numPoints,N);
for curEllips=1:N
    %Perform an eigenvalue decomposition of A.
    [V,D]=eig(A(:,:,curEllips));
    %We will put the smallest eigenvalue at the start, if it is
    %not already there.
    [~,idx]=sort(diag(D),'ascend');
    V=V(idx,idx);
    D=D(idx,idx);
    a=D(1,1);
    b=D(2,2);

    xBound=sqrt(gammaVal/a);
    x=linspace(-xBound,xBound,numPoints);
    %The real command deals with possible finite precision issues.
    y=real(sqrt((gammaVal-x.^2*a)/b));
    %The centered, axis-aligned ellipse points.
    l=[x,fliplr(x);
       y,-fliplr(y)];
    zp=bsxfun(@plus,V*l,z(idx,curEllips));
    zPoints(:,:,curEllips)=zp(idx,:);%Undo the reordering.
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [L,perm,xsuper,split,tmpsiz] = symfctmex(adjncy, perm, cachsz) %#ok
% [L,perm,xsuper,split,tmpsiz] = symfctmex(X, perm, cachsz)
%
%   Computes sparse symbolic factor L, updated permutation PERM,
%   super-node partition XSUPER, and a splitting of supernodes
%   (SPLIT) to optimize use of the computer cache (assuming
%   CACHSZ*1024 byte available).   TMPSIZ is the amount of floating
%   point working storage that has to be allocated within blkfctmex.
%
%   Invokes ORNL block Cholesky library (Fortran).
%
% **********  INTERNAL FUNCTION OF CHOLTOOL  **********
%
% See also sedumi

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

%Indicate to the user Matlab cannot find the SeDuMi binaries
sedumi_binary_error();

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function test_failed = test_libltfat_dgtreal_long(varargin)
test_failed = 0;

fprintf(' ===============  %s ================ \n',upper(mfilename));

definput.flags.complexity={'double','single'};
[flags]=ltfatarghelper({},definput,varargin);
dataPtr = [flags.complexity, 'Ptr'];

[~,~,enuminfo]=libltfatprotofile;
phaseconv = enuminfo.ltfat_phaseconvention;

fftwflags = struct('FFTW_MEASURE',0,'FFTW_ESTIMATE',64,'FFTW_PATIENT',32,'FFTW_DESTROY_INPUT',1,...
    'FFTW_UNALIGNED',2,'FFTW_EXHAUSTIVE',8,'FFTW_PRESERVE_INPUT',16);

Larr  = [350 350   9   1];
aarr  = [ 10  10   9   1];
Marr  = [ 35  35   3   1];
Warr  = [  1   3   3   1];

for idx = 1:numel(Larr)
    L = Larr(idx);
    W = Warr(idx);
    a = aarr(idx);
    M = Marr(idx);
    M2 = floor(M/2) + 1;

    N = L/a;

    g = randn(L,1,flags.complexity);
    c = cast(randn(M2,N,W)+1i*randn(M2,N,W),flags.complexity);
    cout = complex2interleaved(c);
    f = randn(L,W,flags.complexity);

    fPtr = libpointer(dataPtr,f);
    gPtr = libpointer(dataPtr,g);
    coutPtr = libpointer(dataPtr,cout);

    truec = dgtreal(f,g,a,M);

    funname = makelibraryname('dgtreal_long',flags.complexity,0);
    status = calllib('libltfat',funname,fPtr,gPtr,L,W,a,M,phaseconv.LTFAT_FREQINV,coutPtr);

    res = norm(reshape(truec,M2,N*W) - interleaved2complex(coutPtr.Value),'fro');
    [test_failed,fail]=ltfatdiditfail(res+status,test_failed);
    fprintf(['DGTREAL FREQINV       L:%3i, W:%3i, a:%3i, M:%3i %s %s %s\n'],L,W,a,M,flags.complexity,ltfatstatusstring(status),fail);

    truec = dgtreal(f,g,a,M,'timeinv');
    status = calllib('libltfat',funname,fPtr,gPtr,L,W,a,M,phaseconv.LTFAT_TIMEINV,coutPtr);

    res = norm(reshape(truec,M2,N*W) - interleaved2complex(coutPtr.Value),'fro');
    [test_failed,fail]=ltfatdiditfail(res+status,test_failed);
    fprintf(['DGTREAL TIMEINV       L:%3i, W:%3i, a:%3i, M:%3i %s %s %s\n'],L,W,a,M,flags.complexity,ltfatstatusstring(status),fail);

    % With plan
    c = cast(randn(M2,N,W)+1i*randn(M2,N,W),flags.complexity);
    cout = complex2interleaved(c);
    coutPtr = libpointer(dataPtr,cout);

    plan = libpointer();
    funname = makelibraryname('dgtreal_long_init',flags.complexity,0);
    statusInit = calllib('libltfat',funname,gPtr,L,W,a,M,fPtr,coutPtr,phaseconv.LTFAT_FREQINV,fftwflags.FFTW_MEASURE,plan);

    funname = makelibraryname('dgtreal_long_execute',flags.complexity,0);
    statusExecute = calllib('libltfat',funname,plan);

    funname = makelibraryname('dgtreal_long_done',flags.complexity,0);
    statusDone = calllib('libltfat',funname,plan);

    truec = dgtreal(f,g,a,M);
    res = norm(reshape(truec,M2,N*W) - interleaved2complex(coutPtr.Value),'fro');
    [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
    fprintf(['DGTREAL FREQINV WP    L:%3i, W:%3i, a:%3i, M:%3i %s %s %s\n'],L,W,a,M,flags.complexity,ltfatstatusstring(status),fail);

    %%%%%%
    c = cast(randn(M2,N,W)+1i*randn(M2,N,W),flags.complexity);
    cout = complex2interleaved(c);
    coutPtr = libpointer(dataPtr,cout);

    plan = libpointer();
    funname = makelibraryname('dgtreal_long_init',flags.complexity,0);
    statusInit = calllib('libltfat',funname,gPtr,L,W,a,M,fPtr,coutPtr,phaseconv.LTFAT_TIMEINV,fftwflags.FFTW_MEASURE,plan);

    funname = makelibraryname('dgtreal_long_execute',flags.complexity,0);
    statusExecute = calllib('libltfat',funname,plan);

    funname = makelibraryname('dgtreal_long_done',flags.complexity,0);
    statusDone = calllib('libltfat',funname,plan);

    truec = dgtreal(f,g,a,M,'timeinv');
    res = norm(reshape(truec,M2,N*W) - interleaved2complex(coutPtr.Value),'fro');
    [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
    fprintf(['DGTREAL TIMEINV WP    L:%3i, W:%3i, a:%3i, M:%3i %s %s %s\n'],L,W,a,M,flags.complexity,ltfatstatusstring(status),fail);

    % With plan, new array
    c = cast(randn(M2,N,W)+1i*randn(M2,N,W),flags.complexity);
    cout = complex2interleaved(c);
    coutPtr = libpointer(dataPtr,cout);

    plan = libpointer();
    nullPtr = libpointer();
    funname = makelibraryname('dgtreal_long_init',flags.complexity,0);
    statusInit = calllib('libltfat',funname,gPtr,L,W,a,M,nullPtr,nullPtr,phaseconv.LTFAT_FREQINV,fftwflags.FFTW_ESTIMATE,plan);

    funname = makelibraryname('dgtreal_long_execute_newarray',flags.complexity,0);
    statusExecute = calllib('libltfat',funname,plan,fPtr,coutPtr);

    funname = makelibraryname('dgtreal_long_done',flags.complexity,0);
    statusDone = calllib('libltfat',funname,plan);

    truec = dgtreal(f,g,a,M);
    res = norm(reshape(truec,M2,N*W) - interleaved2complex(coutPtr.Value),'fro');
    [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
    fprintf(['DGTREAL FREQINV WP NA L:%3i, W:%3i, a:%3i, M:%3i %s %s %s\n'],L,W,a,M,flags.complexity,ltfatstatusstring(status),fail);

    %%%%%%
    c = cast(randn(M2,N,W)+1i*randn(M2,N,W),flags.complexity);
    cout = complex2interleaved(c);
    coutPtr = libpointer(dataPtr,cout);

    plan = libpointer();
    funname = makelibraryname('dgtreal_long_init',flags.complexity,0);
    statusInit = calllib('libltfat',funname,gPtr,L,W,a,M,nullPtr,nullPtr,phaseconv.LTFAT_TIMEINV,fftwflags.FFTW_ESTIMATE,plan);

    funname = makelibraryname('dgtreal_long_execute_newarray',flags.complexity,0);
    statusExecute = calllib('libltfat',funname,plan,fPtr,coutPtr);

    funname = makelibraryname('dgtreal_long_done',flags.complexity,0);
    statusDone = calllib('libltfat',funname,plan);

    truec = dgtreal(f,g,a,M,'timeinv');
    res = norm(reshape(truec,M2,N*W) - interleaved2complex(coutPtr.Value),'fro');
    [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
    fprintf(['DGTREAL TIMEINV WP NA L:%3i, W:%3i, a:%3i, M:%3i %s %s %s\n'],L,W,a,M,flags.complexity,ltfatstatusstring(status),fail);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [x, infos] = als_nmf(V, rank, in_options)
% Alternative least squares (ALS) for non-negative matrix factorization (NMF).
%
% The problem of interest is defined as
%
%       min || V - WH ||_F^2,
%       where
%       {V, W, H} > 0.
%
% Given a non-negative matrix V, factorized non-negative matrices {W, H} are calculated.
%
%
% Inputs:
%       V           : (m x n) non-negative matrix to factorize
%       rank        : rank
%       in_options
%           alg     : als: Alternative least squares (ALS)
%
%                   : hals: Hierarchical alternative least squares (Hierarchical ALS)
%                       Reference:
%                           Andrzej Cichocki and PHAN Anh-Huy,
%                           "Fast local algorithms for large scale nonnegative matrix and tensor factorizations,"
%                           IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences,
%                           vol. 92, no. 3, pp. 708-721, 2009.
%
%                   : acc_hals: Accelerated hierarchical alternative least squares (Accelerated HALS)
%                       Reference:
%                           N. Gillis and F. Glineur,
%                           "Accelerated Multiplicative Updates and hierarchical ALS Algorithms for Nonnegative
%                           Matrix Factorization,",
%                           Neural Computation 24 (4), pp. 1085-1105, 2012.
%                           See http://sites.google.com/site/nicolasgillis/.
%                           The corresponding code is originally created by the authors,
%                           Then, it is modifided by H.Kasai.
%
%
% Outputs:
%       x           : non-negative matrix solution, i.e., x.W: (m x rank), x.H: (rank x n)
%       infos       : log information
%           epoch   : iteration nuber
%           cost    : objective function value
%           optgap  : optimality gap
%           time    : elapsed time
%           grad_calc_count : number of sampled data elements (gradient calculations)
%
%
% This file is part of NMFLibrary
%
% Created by H.Kasai on Mar. 24, 2017
%
% Change log:
%
%       Oct. 27, 2017 (Hiroyuki Kasai): Fixed algorithm.
%
%       Apr. 22, 2019 (Hiroyuki Kasai): Fixed bugs.
%
%       May. 20, 2019 (Hiroyuki Kasai): Added initialization module.
%
%       Jun. 24, 2022 (Hiroyuki Kasai): Added momentum acceleration mode and mofified.
%
%       Jul. 12, 2022 (Hiroyuki Kasai): Modified code structures.
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.alg   = 'hals';
    local_options.sub_mode = 'std';
    local_options.alpha = 2;
    local_options.delta = 0.1;
    local_options.inner_max_epoch = 500;
    local_options.inner_max_epoch_parameter = 0.5;
    local_options.beta0 = 0.5;
    local_options.eta = 1.5;
    local_options.gammabeta = 1.01;
    local_options.gammabetabar = 1.005;
    local_options.momentum_h = 0;
    local_options.momentum_w = 0;
    local_options.scaling = true;
    local_options.warm_restart = false;

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % set paramters
    if ~strcmp(options.alg, 'als') && ~strcmp(options.alg, 'hals') ...
       && ~strcmp(options.alg, 'acc_hals')
        fprintf('Invalid algorithm: %s. Therfore, we use hals (i.e., Hierarchical ALS).\n', options.alg);
        options.alg = 'hals';
    end

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    W = init_factors.W;
    H = init_factors.H;

    % initialize
    method_name = sprintf('ALS (%s:%s)', options.alg, options.sub_mode);
    epoch = 0;
    grad_calc_count = 0;

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    % intialize for als
    if strcmp(options.alg, 'acc_hals')
        eit1 = cputime;
        VHt = V*H';
        HHt = H*H';

        scaling = sum(sum(VHt.*W))/sum(sum( HHt.*(W'*W) ));
        W = W * scaling;

        options_halsupdt = [];
    end

    if options.scaling
        [W, H] = normalize_WH(V, W, H, rank, 'type1');
    end

    [options, beta, betamax] = check_momemtum_setting(options);

    if options.warm_restart
        nV = norm(V, 'fro');
        rel_error = zeros(1, options.max_epoch);
        rel_error(1) = sqrt(nV^2 - 2*sum(sum(V * H' .* W)) + sum(sum( H * H' .* (W'*W)))) / nV;
    end
    W_prev = W;
    H_prev = H;

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(V, W, H, [], options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('%s: Epoch = 0000, cost = %.16e, optgap = %.4e\n', method_name, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        %% update H
        VtW = V'*W;
        WtW = W'*W;
        WtV = W' * V;

        if strcmp(options.alg, 'als')

            %H = (W*pinv(W'*W))' * V;
            H = WtW \ WtV;        % H = inv(W'*W) * W' * V;
            H = H .* (H>0);

        elseif strcmp(options.alg, 'hals')

            for k=1:rank
                tmp = (VtW(:,k)' - (WtW(:,k)' * H) + (WtW(k,k) * H(k,:))) / WtW(k,k);
                tmp(tmp<=eps) = eps;
                H(k,:) = tmp;
            end

        elseif strcmp(options.alg, 'acc_hals')

            eit1 = cputime;
            options_halsupdt.max_epoch = change_inner_max_epoch(V, W, options);
            H = HALSupdt(H, WtW, WtV, eit1, options.alpha, options.delta, options_halsupdt);

        end

        % perform momentum for H
        if strcmp(options.sub_mode, 'momentum')
            [H, H_tmp1, H_tmp2] = do_momentum_h(H, H_prev, beta, epoch, options);
        end

        %% update W
        VHt = V * H';
        HHt = H * H';

        if strcmp(options.alg, 'als')

            %W = ((inv(H*H')*H)*V')';
            W = VHt / HHt;        % W = V * H' * inv(H*H');
            W = (W>0) .* W;

            % normalize columns to unit
            W = W ./ (repmat(sum(W), m, 1)+eps);

        elseif strcmp(options.alg, 'hals')

            for k=1:rank
                tmp = (VHt(:,k) - (W * HHt(:,k)) + (W(:,k) * HHt(k,k))) / HHt(k,k);
                tmp(tmp<=eps) = eps;
                W(:,k) = tmp;
            end

        elseif strcmp(options.alg, 'acc_hals')

%            if epoch > 0 % Do not recompute A and B at first pass
                % Use actual computational time instead of estimates rhoU
                eit1 = cputime;
                eit1 = cputime-eit1;
%           end
            options_halsupdt.max_epoch = change_inner_max_epoch(V', H', options);
            W = HALSupdt(W', HHt',VHt', eit1, options.alpha, options.delta, options_halsupdt);
            W = W';

        end

        % perform momentum for W
        if strcmp(options.sub_mode, 'momentum')
            [W, H, W_tmp1] = do_momentum_w(W, W_prev, H, H_prev, H_tmp1, beta, epoch, options);
        end

        % perform warm_restart
        if options.warm_restart
            [W, H, W_prev, H_prev, rel_error, beta, betamax, options] = ...
                warm_restart(V, W, H, rank, W_prev, H_prev, W_tmp1, H_tmp1, H_tmp2, rel_error, beta, betamax, epoch, options);
        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % measure gradient calc count
        grad_calc_count = grad_calc_count + m*n;

        % update epoch
        epoch = epoch + 1;

        % store info
        infos = store_nmf_info(V, W, H, [], options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);

    end

    x.W = W;
    x.H = H;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function my_net_backward(net_info, work_info_batch, data_info, net_run_config)

    gpu_mode=net_run_config.use_gpu;

    if gpu_mode
      if ~net_info.ref.net_on_gpu
           my_move_net(net_info, 'gpu') ;
      end
      data_info.ref.output_info_layers{end}=...
          move_output_info_gpu(data_info.ref.output_info_layers{end});
    else
      assert(~net_info.ref.net_on_gpu);
      data_info.ref.output_info_layers{end}=...
          move_output_info_cpu(data_info.ref.output_info_layers{end});
    end

    do_backward(net_info, work_info_batch, data_info, net_run_config);

    if gpu_mode
        if ~net_info.ref.net_stay_on_gpu
            my_move_net(net_info, 'cpu') ;
        end
    end

end

function do_backward(net_info, work_info_batch, data_info, net_run_config)

one_optimizer_param=work_info_batch.ref.gen_optimizer_param_fn(work_info_batch, net_info);
net_info.ref.current_lr=one_optimizer_param.learning_rate;

bp_start_layer=net_info.ref.bp_start_layer;
layer_num = numel(net_info.ref.layers) ;
assert(bp_start_layer<=layer_num);

gpu_mode=net_run_config.use_gpu;

if gpu_mode && net_run_config.sync
    wait(gpuDevice) ;
end

if gpu_mode
    if ~net_info.ref.data_stay_on_gpu
        data_info.ref.output_info_layers{end}=move_output_info_gpu(...
            data_info.ref.output_info_layers{end});
    end
end

for layer_idx=layer_num:-1:bp_start_layer

    input_info=data_info.ref.output_info_layers{layer_idx};
    output_info=data_info.ref.output_info_layers{layer_idx+1};

    if gpu_mode && ~net_info.ref.data_stay_on_gpu
        input_info=move_output_info_gpu(input_info);
    end

    l = net_info.ref.layers{layer_idx} ;
    is_simple_layer= ~strcmp(l.type, 'my_custom');

    if is_simple_layer
        assert(~input_info.is_group_data);
        assert(~output_info.is_group_data);
        switch l.type
          case 'conv'

            one_dzdw=cell(2, 1);
            [input_info.dzdx, one_dzdw{1}, one_dzdw{2}] = ...
                vl_nnconv(input_info.x, l.filters, l.biases, ...
                          output_info.dzdx, ...
                          'pad', l.pad, 'stride', l.stride) ;
            input_info.dzdw=one_dzdw;

          case 'pool'
            input_info.dzdx = vl_nnpool(input_info.x, l.pool, output_info.dzdx, ...
              'pad', l.pad, 'stride', l.stride, 'method', l.method) ;
          case 'normalize'
            input_info.dzdx = vl_nnnormalize(input_info.x, l.param, output_info.dzdx) ;
          case 'softmax'
            input_info.dzdx = vl_nnsoftmax(input_info.x, output_info.dzdx) ;
          case 'relu'
            input_info.dzdx = vl_nnrelu(input_info.x, output_info.dzdx) ;
          case 'noffset'
            input_info.dzdx = vl_nnnoffset(input_info.x, l.param, output_info.dzdx) ;
          case 'dropout'
              input_info.dzdx = vl_nndropout(input_info.x, output_info.dzdx, 'mask', output_info.aux) ;
        end
    else
        if ~check_valid_net_output(output_info)
            break;
        end
        input_info = l.backward_fn(input_info, l, work_info_batch, output_info);
    end

    if gpu_mode && net_run_config.sync
      wait(gpuDevice) ;
    end

    input_info.bp_finished=true;
    data_info.ref.output_info_layers{layer_idx}=input_info;
    do_bp_update_one_layer(net_info, work_info_batch, input_info, layer_idx, one_optimizer_param);

    data_info.ref.output_info_layers{layer_idx+1}=[];

    if gpu_mode && net_run_config.sync
      wait(gpuDevice) ;
    end
end

end

function do_bp_update_one_layer(net_info, work_info_batch, input_info, layer_idx, one_optimizer_param)

    bp_start_layer=net_info.ref.bp_start_layer;
    if layer_idx<bp_start_layer
        return;
    end

      ly=net_info.ref.layers{layer_idx} ;

      if strcmp(ly.type, 'conv')

          lr=one_optimizer_param.learning_rate;
          momentum_param=one_optimizer_param.momentum;
          weightDecay_param=one_optimizer_param.weightDecay;

          ly.filtersMomentum = momentum_param * ly.filtersMomentum ...
              - weightDecay_param * ly.filtersWeightDecay ...
                  * lr * ly.filtersLearningRate * ly.filters ...
              - lr * ly.filtersLearningRate * input_info.dzdw{1} ;

          ly.biasesMomentum = momentum_param * ly.biasesMomentum ...
              - weightDecay_param * ly.biasesWeightDecay ...
                  * lr * ly.biasesLearningRate * ly.biases ...
              - lr * ly.biasesLearningRate * input_info.dzdw{2} ;

          ly.filters = ly.filters + ly.filtersMomentum ;
          ly.biases = ly.biases + ly.biasesMomentum ;
          net_info.ref.layers{layer_idx} = ly ;
      end

      if strcmp(ly.type, 'my_custom')
          layer_update_fn=ly.layer_update_fn;
          if ~isempty(layer_update_fn)

              update_info=[];
              update_info.input_info=input_info;
              update_info.work_info_batch=work_info_batch;

              ly=layer_update_fn('bp_update', ly, net_info, update_info);
              net_info.ref.layers{layer_idx} = ly ;
          end
      end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function kml_contour(lon, lat, z, varargin)
% KML_CONTOUR      overlay MATLAB contour lines onto Google Earth
%
% Syntax:
%     KML_CONTOUR(LON,LAT,Z) writes contour lines in the same format as
%     matlab's CONTOUR(LON,LAT,Z) or CONTOURC(LON,LAT,Z).
%     KML_CONTOUR(LON,LAT,Z,N) draws N contour lines, overriding the
%     automatic value
%     KML_CONTOUR(LON,LAT,Z,V) draws LENGTH(V) contour lines at the values
%     specified in the vector V
%     KML_CONTOUR(LON,LAT,Z,[v v]) computes a single contour at the level v
%
% Input:
%     LON: This can be either a matrix the same size as Z or a vector with
%     length the same as the number of columns in Z.
%     LAT: This can be either a matrix the same size as Z or a vector with
%     length the same as the number of rows in Z.
%     Z: Matrix of elevations
%
% Output:
%     This function creates a kml file called 'doc.kml' in the current
%     working directory
%

%
% Cameron Sparr - Nov. 10, 2011
% cameronsparr@gmail.com
%

    % STYLES:
    % Edit the width, color, and labelsize as you see fit.
    %       color: MATLAB color value (default is 'w')
    %       width: int or float (default is 1)
    %       labelsize: size of text contour labels (default is 0.9)
    color = ge_color('w');
    width = 1;
    labelsize = 0.9;

    % FEEL FREE TO PLAY AROUND WITH THE VALUES SPECIFIED BELOW FOR
    % 'labellimit', 'labelspace', and 'contourlimit'
    %
    % Limit to how many points on a contour line are required for the
    % function to place an altitude label:
    labellimit = round(sqrt(sqrt(numel(z))));
    % Spacing between altitude labels:
    labelspace = labellimit * 8;
    % Contour lines with length below the following limit will not be drawn
    contourlimit = round(labellimit / 3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % output kml file:
    kmlfile = fopen('doc.kml', 'w');
    if size(lon) == size(z)
        matlabel = true;
        [c, h] = contour(lon, lat, z, varargin{:});
        %close(gcf);
    else
        matlabel = false;
        c = contourc(lon, lat, z, varargin{:});
    end

    % specifies position in the 'c' matrix returned from MATLAB's CONTOUR
    % function.
    ind = 1;

    % begin writing kml file:
    kml_begin(kmlfile, labelsize, color, width);

    while ind < length(c)
        % altitude of current contour line
        zz = c(1,ind);
        % length of current contour line
        s = c(2,ind);

        % current contour line
        clon = c(1, ind+1 : ind+s);
        clat = c(2, ind+1 : ind+s);

        if ~matlabel
            % come up with my own contour labelling scheme...
            plon = c(1, ind+1);
            plat = c(2, ind+1);
            if s > labellimit
                place_label(kmlfile, plon, plat, zz);
            end

            for j = 1:numel(clon)
                if mod(j, labelspace) == 0
                    place_label(kmlfile, clon(j), clat(j), zz);
                end
            end
        end

        line_begin(kmlfile, zz);
        for i = 1:numel(clon)
            lon = num2str(clon(i), 8);
            lat = num2str(clat(i), 8);
            fprintf(kmlfile, [lon, ',', lat, ',',num2str(zz),' ']);
        end
        line_end(kmlfile);

        % move ind up to the next contour line.
        ind = ind + s + 1;
    end

    if matlabel
        disp('Using MATLAB native contour label positions');
        % use matlab positions for contour labels.
        lh = clabel(c,h);
        pos = get(lh, 'position');
        height = get(lh, 'UserData');
        close(gcf);
        for ii = 1:length(pos)
            mlon = pos{ii}(1);
            mlat = pos{ii}(2);
            zzz = height{ii}(1);
            place_label(kmlfile, mlon, mlat, zzz);
        end
    else
        disp('Using proprietary contour label positions');
    end

    kml_end(kmlfile);
end

function kml_begin(kmlfile, labelsize, color, width)
    fprintf(kmlfile, '<?xml version="1.0" encoding="UTF-8"?>\n');
    fprintf(kmlfile, '<kml xmlns="http://www.opengis.net/kml/2.2"');
    fprintf(kmlfile, ' xmlns:gx="http://www.google.com/kml/ext/2.2"');
    fprintf(kmlfile, ' xmlns:kml="http://www.opengis.net/kml/2.2"');
    fprintf(kmlfile, ' xmlns:atom="http://www.w3.org/2005/Atom">\n');
    fprintf(kmlfile, '<Document>\n');
    fprintf(kmlfile, '	<name>doc.kml</name>\n');

    fprintf(kmlfile, '	<Style id="sn_noicon">\n');
    fprintf(kmlfile, '      <IconStyle>\n');
    fprintf(kmlfile, '          <Icon>\n');
    fprintf(kmlfile, '          </Icon>\n');
    fprintf(kmlfile, '      </IconStyle>\n');
    fprintf(kmlfile, '      <LabelStyle>\n');
    fprintf(kmlfile, '          <scale>');
    fprintf(kmlfile, num2str(labelsize));
    fprintf(kmlfile, '</scale>\n');
    fprintf(kmlfile, '      </LabelStyle>\n');
    fprintf(kmlfile, '	</Style>\n');

    fprintf(kmlfile, '	<Style id="linestyle">\n');
    fprintf(kmlfile, '			<LineStyle>\n');
    fprintf(kmlfile, '				<color>#FF');
    fprintf(kmlfile, color);
    fprintf(kmlfile, '</color>\n');
    fprintf(kmlfile, '				<width>');
    fprintf(kmlfile, num2str(width));
    fprintf(kmlfile, '</width>\n');
    fprintf(kmlfile, '			</LineStyle>\n');
    fprintf(kmlfile, '	</Style>\n');
end

function kml_end(kmlfile)
    fprintf(kmlfile, '</Document>\n');
    fprintf(kmlfile, '</kml>\n');
    fclose(kmlfile);
end

function line_begin(kmlfile, zz)
    fprintf(kmlfile, '	<Placemark>\n');
    fprintf(kmlfile, ['		<name>', num2str(zz),'</name>\n']);
    fprintf(kmlfile, '		<styleUrl>#linestyle</styleUrl>\n');
    fprintf(kmlfile, '		<LineString>\n');
    fprintf(kmlfile, '			<tessellate>1</tessellate>\n');
    fprintf(kmlfile, '			<altitudeMode>clampToSeaFloor</altitudeMode>\n');
    fprintf(kmlfile, '			<gx:altitudeMode>clampToSeaFloor</gx:altitudeMode>\n');
    fprintf(kmlfile, '			<coordinates>\n');
    fprintf(kmlfile, '				');
end

function line_end(kmlfile)
    fprintf(kmlfile, '\n');
    fprintf(kmlfile, '			</coordinates>\n');
    fprintf(kmlfile, '		</LineString>\n');
    fprintf(kmlfile, '	</Placemark>\n');
end

function place_label(kmlfile, plon, plat, z)
    z = round(z);
    fprintf(kmlfile, '	<Placemark>\n');
    fprintf(kmlfile, ['		<name>', num2str(z),'</name>\n']);
    fprintf(kmlfile, '		<styleUrl>#sn_noicon</styleUrl>\n');
    fprintf(kmlfile, '		<Point>\n');
    fprintf(kmlfile, '			<altitudeMode>clampToSeaFloor</altitudeMode>\n');
    fprintf(kmlfile, '			<gx:altitudeMode>clampToSeaFloor</gx:altitudeMode>\n');
    fprintf(kmlfile, ['			<coordinates>',num2str(plon, 8),',',num2str(plat, 8),',0</coordinates>\n']);
    fprintf(kmlfile, '		</Point>\n');
    fprintf(kmlfile, '	</Placemark>\n');
end

function clrstr=ge_color(c,varargin)
%Jarrell Smith
%3/4/2008
    opacity=1;
    cspec=[0,0,0];

    nargchk(nargin,1,2);
    if nargin==2,
       mode='both';
       opacity=varargin{1};
       if length(opacity)>1 || ~isnumeric(opacity),
          error('Opacity must be numeric and length 1')
       elseif opacity>1 || opacity<0,
          error('Opacity must be between 0-1')
       end
    else
       mode='color';
    end
    if ischar(c), %process as color
       switch lower(c)
          case {'y','yellow'}
             cspec=[1,1,0];
          case {'m','magenta'}
             cspec=[1,0,1];
          case {'c','cyan'}
             cspec=[0,1,1];
          case {'r','red'}
             cspec=[1,0,0];
          case {'g','green'}
             cspec=[0,1,0];
          case {'b','blue'}
             cspec=[0,0,1];
          case {'w','white'}
             cspec=[1,1,1];
          case {'k','black'}
             cspec=[0,0,0];
          otherwise
             error('%s is an invalid Matlab ColorSpec.',c)
       end
    elseif isnumeric(c) && ndims(c)==2, %Determine if Color or Opacity
       if  all(size(c)==[1,1]), %Input is Opacity
          if c>1 || c<0
             error('Opacity must be scalar quantity between 0 to 1')
          end
          opacity=c;
          mode='opacity';
       %color
       elseif all(size(c)==[1,3]) %Input is Color
          if any(c<0|c>1)
             error('Numeric ColorSpec must be size [1,3] with values btw 0 to 1.')
          end
          cspec=c;
       else
          error('Incorrect size of first input argument.  Size must be [1,3] or [1,1].')
       end
    else
       error('Incorrect size of first input argument.  Size must be [1,3] or [1,1].')
    end
    opacity=round(opacity*255); %transparency (Matlab format->KML format)
    cspec=round(fliplr(cspec)*255); %color (Matlab format->KML format)
    switch mode
       case 'color'
          clrstr=sprintf('%s%s%s',dec2hex(cspec,2)');
       case 'opacity'
          clrstr=sprintf('%s',dec2hex(opacity,2));
       case 'both'
          clrstr=sprintf('%s%s%s%s',dec2hex(opacity,2),dec2hex(cspec,2)');
    end
end


text:
function r = logdet(M)
% Compute log(det(A)) without the usual numerical inaccuracies.

% Copyright (C) Christian Kothe, SCCN, 2011, christian@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify it under the terms of the GNU
% General Public License as published by the Free Software Foundation; either version 2 of the
% License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without
% even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License along with this program; if not,
% write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
% USA

r = 2 * sum(log(diag(chol(M))), 1);

repo: REST


text:
function showscene(cameras,voxels)
%SHOWSCENE: show a carve scene, including cameras, images and the model
%
%   SHOWSCENE(CAMERAS,VOXELS) shows the specified list of CAMERAS and the
%   current model as a surface fitted around VOXELS.
%
%   Example:
%   >> camera = loadcameradata(1);
%   >> camera.Silhouette = getsilhouette( camera.Image );
%   >> voxels = carve( makevoxels(50), camera );
%   >> showscene( camera, voxels );

%   Copyright 2005-2009 The MathWorks, Inc.
%   $Revision: 1.0 $    $Date: 2006/06/30 00:00:00 $

error( nargchk( 1, 2, nargin ) );

%% Plot each camera centre
set(gca,'DataAspectRatio',[1 1 1])
hold on

N = numel(cameras);
for ii=1:N
    spacecarving.showcamera(cameras(ii));
end
xlabel('X')
ylabel('Y')
zlabel('Z')

%% And show a surface around the object
if nargin>1 && ~isempty( voxels )
    spacecarving.showsurface( voxels );
end

%% Light the scene and adjust the view angle to make it a bit easier on
% the eye
view(3);
grid on;
light('Position',[0 0 1]);
axis tight


text:
function [x,tspan] = w20scalar(f,L,tspan,x0,Q,gaussian)
%% W20scalar - Weak order 2.0 Ito-Taylor
%
% Syntax:
%   [x,tspan] = w20scalar(f,L,tspan,x0,Q,gaussian)
%
% In:
%   f        - Drift function and derivatives, {f, df, ddf}
%   L        - Diffusion function and derivatives, {L, dL, ddL}
%   tspan    - Time steps to simulate, [t0,...,tend]
%   x0       - Initial condition
%   Q        - Spectral density (default: standard Brownian motion)
%   gaussian - Use Gaussian increments (default: false)
%
% Out:
%   x      - Solved values
%   tspan  - Time steps
%
% Description:
%   Integrates the system of stochatic differential equations
%     dx = f(x,t) dt + L(x,t) dbeta,  for x(0) = x0
%   over the time interval defined in tspan.
%
% Copyright:
%   2018 - Simo Särkkä and Arno Solin
%
% License:
%   This software is provided under the MIT License. See the accompanying
%   LICENSE file for details.

%%

  % Check if Q given
  if nargin<5 || isempty(Q), Q = eye(size(L(x0,tspan(1)),2)); end

  % Check the sort of increments to use
  if nargin<6 || isempty(gaussian), gaussian = false; end

  % Extract derivatives
  ddf = f{3};
  df = f{2};
  f = f{1};
  ddL = L{3};
  dL = L{2};
  L = L{1};

  % Cholesky factor of Q
  cQ = chol(Q,'lower');

  % Number of steps
  steps = numel(tspan);

  % Allocate space
  x = zeros(size(x0,1),steps);

  % Initial state
  x(:,1) = x0;

  % Pre-calculate random numbers
  RR = rand(size(Q,1),steps);
  R = (RR<1/6)*sqrt(3) - (RR>5/6)*sqrt(3);

  % Iterate
  for k=2:steps

    % Time discretization
    dt = tspan(k)-tspan(k-1);

    % Increment
    if gaussian
      db = sqrt(dt)*cQ*randn;
      dz = 1/2*db*dt;
    else
      db = sqrt(dt)*cQ*R(:,k);
      dz = 1/2*dt^(3/2)*R(:,k);
    end

    % Evaluate only once
    fx = f(x(:,k-1),[]);
    dfx = df(x(:,k-1),[]);
    ddfx = ddf(x(:,k-1),[]);
    Lx = L(x(:,k-1),[]);
    dLx = dL(x(:,k-1),[]);
    ddLx = ddL(x(:,k-1),[]);

    % Step
    x(:,k) = x(:,k-1) + ...
        fx*dt + ...
        Lx*db + ...
        1/2*Lx*dLx*(db^2 - dt) + ...
        Lx*df(x(:,k-1),[])*dz + ...
        1/2*(fx*dfx + 1/2*Lx^2*ddfx)*dt^2 + ...
        (fx*dLx - 1/2*Lx^2*ddLx)*(db*dt - dz);

  end

author: AaltoML

repo: SDE

save_path: github-repos/MATLAB/AaltoML-SDE


text:
% Calculate the log mel filter bank from waveform
% Author: Xiao Xiong
% Created: 4 Feb 2005
% Last modified: 4 Feb 2005

function [fbank] = wav2fbank(x, fs, frame_shift, nMel)
if nargin < 2
    fs = 8000;
end
if nargin < 3
    frame_shift = 0.01;
end
if nargin<4
    nMel = 23;
end

absX = wav2abs(x,fs,frame_shift);
fbank = log (abs2Mel(absX, fs, nMel));

author: singaxiong

repo: SignalGraph


text:
function regressedBox = BoxRegresssLog(box, regressionFactor)
% regressedBox = BoxRegresss(box, regressionFactor)
%
% Apply the regressionFactor to the box to get better bounding box:
% regressedBox(4) = boxMiddle + 1/2 boxWidth + regressionFactor(4) * 1/2 boxWidth
%                 = boxMiddle + (1+regressionFactor(4)) * 1/2 boxWidth
%
% box:              N x 4 vector with BB coordinates
% regressionFactor: N x 4 vector with regression factors
%
% regressedBox:     N x 4 vector with updated BB coordinates
%
% In log-space, like Girshick
%
% Jasper Uijlings - 2015

% regressedBox(N) = boxMiddle + (1+regressionFactor(N)) * 1/2 boxWidth
regressionFactor = 2 * (exp(regressionFactor / 6.537) - 0.5);

% Obtain middle
middleOdd = (box(:,1) + box(:,3)) / 2;
middleEven = (box(:,2) + box(:,4)) / 2;

% Apply factors with respect to middle.
% regress coord = middle       + distance to middle      .* F
regressedBox(:,4) = middleEven + (box(:,4) - middleEven) .* regressionFactor(:,4);
regressedBox(:,3) = middleOdd  + (box(:,3) - middleOdd)  .* regressionFactor(:,3);
regressedBox(:,2) = middleEven + (box(:,2) - middleEven) .* regressionFactor(:,2);
regressedBox(:,1) = middleOdd  + (box(:,1) - middleOdd)  .* regressionFactor(:,1);

author: nightrome

repo: matconvnet-calvin


text:
function hermite_polynomial_test11 ( p, e )

%*****************************************************************************80
%
%% HERMITE_POLYNOMIAL_TEST11 tests HEN_POWER_PRODUCT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer P, the maximum degree of the polynomial factors.
%
%    Input, integer E, the exponent of X in the integrand.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_POLYNOMIAL_TEST11\n' );
  fprintf ( 1, '  Compute a normalized probabilist''s Hermite polynomial power product table.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Tij = integral ( -oo < X < +oo ) X^E Hen(I,X) Hen(J,X) exp(-0.5*X*X) dx\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  where Hen(I,X) = normalized probabilist''s Hermite polynomial of degree I.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Maximum degree P = %d\n', p );
  fprintf ( 1, '  Exponent of X = %d\n', e );

  table = hen_power_product ( p, e );

  r8mat_print ( p + 1, p + 1, table, '  Power weighted table:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%   AUTHORSHIP
%   Primary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Math Lead & Secondary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%
function [jd,pane]=msgError(txt, pause, where, title)
if nargin<4
    title='Error...';
    if nargin<3
        where='center';
        if nargin<2
            pause=0;
        end
    end
end
[jd, pane]=msg(txt, pause, where, title, 'error.png');
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% Weighted random sampling.
%
% INPUTs: number of draws from a discrete distribution (n)
%         possible values to pick from, (P)
%         set of normalized weights/probabilities, (W)
% OUTPUTs: s - set of n numbers drawn from P
%              according to the weights in W
%
% GB: last updated, Oct 31 2012

function s = weightedRandomSample(n,P,W)

s = [];

if abs(sum(W)-1)>10**(-8);
    printf('The probabilities do not sum up to 1.\n');
    s = 'probabilities do not sum up to 1';
    return;
end

% divide the unit interval into |P| segments each with length W_i
unit = [0,W(1)];
for w=2:length(W)
  unit = [unit W(w)+unit(length(unit))];
end

% draw a random number in the unit interval uniformly - where does it fall?
while length(s)<n

  lb = max(find(unit<rand));    % rand is in [unit(lb), unit(lb+1)]
  s = [s P(lb)];                % pick P(lb)

end

author: aeolianine

repo: octave-networks-toolbox


text:
classdef MulConst < dagnn.ElementWise
% MulConst multiplies its input by its parameter.

    properties
    end

    methods
        function outputs = forward(obj, inputs, params)
            assert(numel(inputs) == 1, 'one input is needed');
            assert(numel(params) == 1, 'one param is needed');
            outputs{1} = mul_const(inputs{1}, params{1});
        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
            assert(numel(inputs) == 1, 'one input is needed');
            assert(numel(params) == 1, 'one param is needed');
            assert(numel(derOutputs) == 1, 'expect one gradient');
            [derInputs{1}, derParams{1}] = mul_const(...
                inputs{1}, params{1}, derOutputs{1});
        end

        function obj = MulConst(varargin)
            obj.load(varargin);
        end
    end
end

author: bertinetto

repo: cfnet

save_path: github-repos/MATLAB/bertinetto-cfnet


text:
%NAVIERSTOKES Navier-Stokes equations physics mode.

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function Cu = computeCu(nu)
    Cu_tmp1 = cell(nu,nu);
    Cu_tmp2 = cell(nu,nu);

    for i=1:nu
       for j=1:nu
          if i==j
              Cu_tmp1{i,j} = eye(2);
              Cu_tmp2{i,j} = -eye(2);
          else
              Cu_tmp1{i,j} = zeros(2,2);
              Cu_tmp2{i,j} = zeros(2,2);
          end
       end
    end

    Cu = [cell2mat(Cu_tmp1);cell2mat(Cu_tmp2)];
end

author: ccalas

repo: mpc

save_path: github-repos/MATLAB/ccalas-mpc


text:
function [vclassifier, hclassifier] = mcmcTrainSegmentClassifier2(features, labels, weights, maxdata, classparams)

if exist('classparams') && ~isempty(classparams)
    nnodes = classparams(1);
    ntrees = classparams(2);
    stopval = classparams(3);
else
    nnodes = 8;
    ntrees = 20;
    stopval = 0;
end

if ~exist('maxdata') || isempty(maxdata)
    maxdata = 25000;
end

if ~exist('weights') || isempty(weights)
    for f = 1:numel(labels)
        weights{f} = ones(size(labels{f})) / numel(labels{f});
    end
end

vlabels = labels;
hlabels = labels;
for f = 1:numel(labels)
    vlabels{f} = 1*(labels{f}==1) + 2*((labels{f}>1) & (labels{f}<7)) + ...
        3*(labels{f}==7);
    hlabels{f} = (labels{f}-1).*((labels{f}>1) & (labels{f}<7));
end

[vdata, vlab, vw] = formatData(features, vlabels, weights, maxdata);
[hdata, hlab, hw] = formatData(features, hlabels, weights, maxdata);

vnames = {'000', '090', 'sky'};
hnames = {'045', '090', '135', 'por', 'sol'};

%if size(vdata, 2) > 75
catids = []; % 75
%end
%num2str(find(std(vdata, 1)==0))
vclassifier =  train_boosted_dt_mc(vdata, catids, vnames(vlab)', ntrees, nnodes, stopval, vw, vnames);

hclassifier =  train_boosted_dt_mc(hdata, catids, hnames(hlab)', ntrees, nnodes, stopval, hw, hnames);

%% Reformat the input data to be used by classifier
function [data, lab, w] = formatData(features, labels, weights, maxdata)
% concatenate data

nimages = numel(features);

[tmp, nvars] = size(features{1});

% count segments
nseg = 0;
for f = 1:nimages
    nseg = nseg + sum(labels{f}>0);
end
%disp(num2str(nseg))

data = zeros(nseg, nvars);
lab = zeros(nseg, 1);
w = zeros(nseg, 1);

% concatenate data
vc = 0;
for f = 1:nimages
    ind = find(labels{f}>0);
    %disp(num2str([f size(labels{f}) size(features{f})]))
    data(vc+1:vc+numel(ind), :) = features{f}(ind, :);
    lab(vc+1:vc+numel(ind)) = labels{f}(ind);
    w(vc+1:vc+numel(ind)) = weights{f}(ind); % weight according to area in image
    vc = vc + numel(ind);
end
%disp([vc nseg])

if nseg > maxdata
    rind = randperm(nseg);
    rind = rind(1:maxdata);
    data = data(rind, :);
    lab = lab(rind);
    w = w(rind);
    nseg = maxdata;
end

w = w / sum(w);

author: zouchuhang

repo: LayoutNet


text:
function write_descriptors(path, descriptors)
% WRITE_DESCRIPTORS - Write the descriptors to a binary file.
%   path:
%       Path to descriptor file.
%   descriptors:
%       Floating point descriptors, where each row represents one
%       descriptor and each column represents one dimension of the
%       descriptor. For example, size(descriptors) = (1000, 128) contains
%       1000 descriptors of dimensionality 128. The data type of
%       the descriptors must be a single-precision real matrix.
%
% Copyright 2017: Johannes L. Schoenberger <jsch at inf.ethz.ch>

assert(isreal(descriptors) & isfloat(descriptors));

fid = fopen(path, 'w');
fwrite(fid, size(descriptors), 'int32');
fwrite(fid, descriptors', 'single');
fclose(fid);

end

author: ahojnnes

repo: local-feature-evaluation


text:
function pols = klegeypols ( x, y, n )

%*****************************************************************************80
%
%% KLEGEYPOLS evaluates scaled Legendre polynomials.
%
%  Discussion:
%
%    This routine evaluate a sequence of scaled Legendre polynomials
%    P_n(x/y) y^n, with the parameter y in [0,1].
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    13 June 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Input, real Y, the parameter.
%
%    Input, integer N, the highest degree to be evaluated.
%
%    Output, real POLS(N+1), the polynomial values.
%
  pols = zeros(n+1,1);

  pkp1 = 1.0;
  pols(1) = pkp1;
  if ( n == 0 )
    return
  end

  pk = pkp1;
  pkp1 = x;
  pols(2) = pkp1;
  if ( n == 1 )
    return
  end

  for k = 1 : n - 1
    pkm1 = pk;
    pk = pkp1;
    pkp1 = ( ( 2.0 * k + 1.0 ) * x * pk ...
      - k * pkm1 * y * y ) / ( k + 1.0 );
    pols(k+2) = pkp1;
  end

  return
end


text:
function pval = lme_mass_LR(statsfull,statsred,q)
% pval = lme_mass_LR(statsfull,statsred,q)
%
% Likelihood ratio test for the random effects. It can be used to test if a
% model with q+1 random effects is significantly better than a model with q
% random effects.
%
% Input
% statsfull: Structure array containing statistiscs for every voxel/vertex
% (see lme_FSfit for more details on these statistics) for the full
% model (the one with q+1 random effects).
% statsred: Structure array containing statistiscs for every voxel/vertex
% for the reduced model (the one with q random effects).
% q: Number of random effects in the reduced model.
%
% Output
% pval: P-value of the test at each voxel/vertex (based on a 50:50 mixture
% of chi-squared distributions with q and q+1 degrees of freedom).
%
% $Revision: 1.2 $  $Date: 2015/01/06 17:14:55 $
% Original Author: Jorge Luis Bernal Rusiel
% CVS Revision Info:
%    $Author: mreuter $
%    $Date: 2015/01/06 17:14:55 $
%    $Revision: 1.2 $
% References: Bernal-Rusiel J.L., Greve D.N., Reuter M., Fischl B., Sabuncu
% M.R., 2012. Statistical Analysis of Longitudinal Neuroimage Data with Linear
% Mixed Effects Models, NeuroImage, doi:10.1016/j.neuroimage.2012.10.065.
%

if nargin < 3
    error('Too few inputs');
end;

nv = length(statsfull);
pval = zeros(1,nv);
for i=1:nv
    lrstats = lme_LR(statsfull(i).lreml,statsred(i).lreml,q);
    pval(i) = lrstats.pval;
end;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [vertices, faces] = collapseEdgesWithManyFaces(vertices, faces, varargin)
%COLLAPSEEDGESWITHMANYFACES Remove mesh edges adjacent to more than two faces.
%
%   [V2, F2] = collapseEdgesWithManyFaces(V, F)
%   Count the number of faces adjacent to each edge, and collapse the edges
%   adjacent to more than two faces.
%
%
%   Example
%   collapseEdgesWithManyFaces
%
%   See also
%       trimMesh, isManifoldMesh

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2019-01-31, using Matlab 9.5.0.944444 (R2018b)
% Copyright 2019-2022 INRA - Cepia Software Platform

verbose = false;
while length(varargin) > 1 && ischar(varargin{1})
    name = varargin{1};
    if strcmpi(name, 'verbose')
        verbose = varargin{2};
    else
        error(['Unknown optional argument: ' name]);
    end
    varargin(1:2) = [];
end

while true
    % compute edge to vertex mapping
    edges = meshEdges(faces);

    % compute number of faces incident to each edge
    edgeFaces = trimeshEdgeFaces(faces);
    edgeFaceDegrees = sum(edgeFaces > 0, 2);

    inds = find(edgeFaceDegrees > 2);

    if isempty(inds)
        break;
    end

    edge = edges(inds(1), :);
    if verbose
        fprintf('remove edge with index %d: (%d, %d)\n', inds(1), edge);
    end
    [vertices, faces] = mergeMeshVertices(vertices, faces, edge);
end

% trim
[vertices, faces] = trimMesh(vertices, faces);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function map = preset_colormap(N,varargin)
% Wrapper for any COLORMAP function, to provide preset parameter values.
%
% (c) 2020-2022 Stephen Cobeldick
%
%%% Syntax:
% preset_colormap(@fun,p1,p2,...,pN) % store function and any parameters
% map = preset_colormap(N)           % generate colormap
%
%%% Examples %%%
%
% >> preset_colormap(@cubehelix,0.25,-0.67,1.5,1)
% >> colormap(preset_colormap)
% or
% >> preset_colormap(5)
% ans =
%          0         0         0
%     0.1055    0.2788    0.4895
%     0.1660    0.6705    0.4961
%     0.6463    0.8479    0.5076
%     1.0000    1.0000    1.0000
%
% >> preset_colormap(@brewermap,'PuOr')
% >> load topo
% >> load coast
% >> figure
% >> worldmap(topo, topolegend)
% >> contourfm(topo, topolegend);
% >> contourcmap('preset_colormap', 'Colorbar','on', 'Location','horizontal','TitleString','Contour Intervals in Meters');
% >> plotm(lat, long, 'k')
%
% See Also BREWERMAP CUBEHELIX COLORMAP CONTOURCMAP

persistent fnh arg
%
if nargin==0 % Default N same as MATLAB colormaps.
	N = cmDefaultN();
elseif nargin==1 && isnumeric(N)
	assert(isscalar(N) && isreal(N) && fix(N)==N,...
		'SC:preset_colormap:N:NotRealScalarNumeric',...
		'First input <N> must be a real scalar numeric of the colormap size.')
	N = double(N);
else % Store function handle and parameter values.
	assert(isa(N,'function_handle'),...
		'SC:preset_colormap:N:NotFunctionHandle',...
		'First input <N> must be a function handle to a colormap function.')
	fnh = N;
	arg = varargin;
	return
end
%
map = fnh(N,arg{:});
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%preset_colormap
function N = cmDefaultN()
% Get the colormap size from the current figure or default colormap.
try
	F = get(groot,'CurrentFigure');
catch %#ok<CTCH> pre HG2
	N = size(get(gcf,'colormap'),1);
	return
end
if isempty(F)
	N = size(get(groot,'DefaultFigureColormap'),1);
else
	N = size(F.Colormap,1);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%cmDefaultN

author: DrosteEffect

repo: BrewerMap


text:
% OtQueryFree Disconnect from an Omnitrans query.
%
% Usage:
%    OtQueryFree(hquery)
%
% Releases a query handle created with OtQueryNew().

% OtMatlab: Matlab interface to the Omnitrans transport planning software

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

% All rights reserved.

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

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

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

% 2006-02-28 Created (Johan Meijdam, VORtech)

function OtQueryFree(hquery)
   Ot('query.free', hquery);
end

author: Sable

repo: mcbench-benchmarks


text:
function overlayMaps(vw,map1,map2,anatFlag,map1Name,map2Name)
%
% overlayMaps(view,map1,map2,[anatFlag],[map1Name,map2Name]);
%
% Given two volumes containing values
% relevant for the current view (such
% as param maps, corAnal fields, etc),
% grab the relevant anatomical images from
% the view, take into account relevant
% clip values, zoom, and other factors, and
% overlay using the overlayVolumes GUI.
%
% map1Name and map2Name are optional strings
% specifying the names to use in the
% GUI for the two volumes.
%
%
% ras 04/02/2005
if ieNotDefined('map1Name')
    map1Name = inputname(2);
end

if ieNotDefined('map2Name')
    map2Name = inputname(3);
end

if ieNotDefined('anatFlag')
    anatFlag = 1;
end

if anatFlag==1
    % initialize anatomy
	viewType = viewGet(vw,'viewType');
	switch viewType
        case {'Inplane','Flat'},
            nSlices = numSlices(vw);
            for i = 1:nSlices
                bg(:,:,i) = recomputeAnatImage(vw,[],i);
            end

        case {'Volume','Gray'},
            bg = viewGet(vw,'anat');
            ori = getCurSliceOri(vw);
            if ori==1, % permute to axial view
                bg = permute(bg,[2 3 1]);
            elseif ori==2, % permute to coronal view
                bg = permute(bg,[1 3 2]);
            end
	end

    % resize overlay volumes to be same size as bg
    switch viewType
        case {'Inplane'},
            vs = viewGet(vw,'Size');
            for i = 1:vs(3)
                tmp1(:,:,i) = imresize(map1(:,:,i),[vs(1) vs(2)]);
                tmp2(:,:,i) = imresize(map2(:,:,i),[vs(1) vs(2)]);
            end
            map1 = tmp1; map2 = tmp2; clear tmp1 tmp2;

        case {'Volume','Gray','volume3View'},
            coords = canOri2CurOri(vw,vw.coords);
            ind = sub2ind(size(bg),coords(1,:),coords(2,:),coords(3,:));

            tmp1 = zeros(size(bg));
            tmp1(ind) = map1;
            map1 = tmp1;
            clear tmp1;

            tmp2 = zeros(size(bg));
            tmp2(ind) = map2;
            map2 = tmp2;
            clear tmp2;

            if isfield(vw.ui,'flipLR') & vw.ui.flipLR==1
                map1 = flipdim(map1,2);
                map2 = flipdim(map2,2);
                bg = flipdim(bg,2);
            end

        case {'Flat'},
            % threshold to see sulci/gyri clearly
            thresh = 0.5 .* max(bg(:));
            lightRng = [0.6 0.8] .* max(bg(:));
            darkRng = [0.2 0.4] .* max(bg(:));
            bg(bg < thresh) = normalize(bg(bg < thresh),darkRng(1),darkRng(2));
            bg(bg >= thresh) = normalize(bg(bg >= thresh),lightRng(1),lightRng(2));
    end

    % apply a zoom, if specified in ui
    ui = viewGet(vw,'ui');
    if isfield(ui,'zoom')
        switch viewType
            case {'Inplane','Flat'},
                zoom = round(ui.zoom);
                xrng = zoom(1):zoom(2);
                yrng = zoom(3):zoom(4);
                map1 = map1(yrng,xrng,:);
                map2 = map2(yrng,xrng,:);
                bg = bg(yrng,xrng,:);
            case {'Volume','Gray'},
                zoom = round(ui.zoom);
                ap = zoom(1,1):zoom(1,2);
                si = zoom(2,1):zoom(2,2);
                rl = zoom(3,1):zoom(3,2);
                map1 = map1(ap,si,rl);
                map2 = map2(ap,si,rl);
                bg = bg(ap,si,rl);
        end
    end

    bg = brighten(bg,0.5);

    overlayVolumes(map1,map2,bg,[],'map1Name',map1Name,'map2Name',map2Name);
else
    overlayVolumes(map1,map2,[],[],'map1Name',map1Name,'map2Name',map2Name);
end

return


text:
function trans = rotationOy(varargin)
%ROTATIONOY return 4x4 matrix of a rotation around y-axis
%
%   TRANS = rotationOy(THETA);
%   Returns the transform matrix corresponding to a rotation by the angle
%   THETA (in radians) around the Oy axis. A rotation by an angle of PI/2
%   would transform the vector [0 0 1] into the vector [1 0 0].
%
%   The returned matrix has the form:
%   [ cos(THETA)  0  sin(THETA)  0 ]
%   [    0        1       0      0 ]
%   [-sin(THETA)  0  cos(THETA)  0 ]
%   [    0        0       0      1 ]
%
%   TRANS = rotationOy(ORIGIN, THETA);
%   TRANS = rotationOy(X0, Y0, Z0, THETA);
%   Also specifies origin of rotation. The result is similar as performing
%   translation(-X0, -Y0, -Z0), rotation, and translation(X0, Y0, Z0).
%
%
%   See also:
%   transforms3d, transformPoint3d, rotationOx, rotationOz
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 18/02/2005.
%

%   HISTORY
%   2008/11/24 changed convention for angle
%   30/04/2009 deprecate: use createRotationOy instead

% deprecation warning
warning('geom3d:deprecated', ...
    [mfilename ' is deprecated, use ''createRotationOy'' instead']);

% call current implementation
trans = createRotationOy(varargin{:});

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips