text:
%SerialLink.ikine6s Analytical inverse kinematics
%
% Q = R.ikine(T) are the joint coordinates (1xN) corresponding to the robot
% end-effector pose T which is an SE3 object or homogenenous transform
% matrix (4x4), and N is the number of robot joints.  This is a analytic
% solution for a 6-axis robot with a spherical wrist (the most common form
% for industrial robot arms).
%
% If T represents a trajectory (4x4xM) then the inverse kinematics is
% computed for all M poses resulting in Q (MxN) with each row representing
% the joint angles at the corresponding pose.
%
% Q = R.IKINE6S(T, CONFIG) as above but specifies the configuration of the arm in
% the form of a string containing one or more of the configuration codes:
%
% 'l'   arm to the left (default)
% 'r'   arm to the right
% 'u'   elbow up (default)
% 'd'   elbow down
% 'n'   wrist not flipped (default)
% 'f'   wrist flipped (rotated by 180 deg)
%
% Trajectory operation::
%
% In all cases if T is a vector of SE3 objects (1xM) or a homogeneous
% transform sequence (4x4xM) then R.ikcon() returns the joint coordinates
% corresponding to each of the transforms in the sequence.
%
% Notes::
% - Treats a number of specific cases:
%   - Robot with no shoulder offset
%   - Robot with a shoulder offset (has lefty/righty configuration)
%   - Robot with a shoulder offset and a prismatic third joint (like Stanford arm)
%   - The Puma 560 arms with shoulder and elbow offsets (4 lengths parameters)
%   - The Kuka KR5 with many offsets (7 length parameters)
% - The inverse kinematics for the various cases determined using ikine_sym.
% - The inverse kinematic solution is generally not unique, and
%   depends on the configuration string.
% - Joint offsets, if defined, are added to the inverse kinematics to
%   generate Q.
% - Only applicable for standard Denavit-Hartenberg parameters
%
% Reference::
% - Inverse kinematics for a PUMA 560,
%   Paul and Zhang,
%   The International Journal of Robotics Research,
%   Vol. 5, No. 2, Summer 1986, p. 32-44
%
% Author::
% - The Puma560 case: Robert Biro with Gary Von McMurray,
%   GTRI/ATRP/IIMB, Georgia Institute of Technology, 2/13/95
% - Kuka KR5 case: Gautam Sinha,
%   Autobirdz Systems Pvt. Ltd.,  SIDBI Office,
%   Indian Institute of Technology Kanpur, Kanpur, Uttar Pradesh.
%
% See also SerialLink.fkine, SerialLink.ikine, SerialLink.ikine_sym.

function thetavec = ikine6s(robot, TT, varargin)

    if robot.mdh ~= 0
        error('RTB:ikine:notsupported','Solution only applicable for standard DH conventions');
    end

    if robot.n ~= 6
        error('RTB:ikine:notsupported','Solution only applicable for 6-axis robot');
    end

    %
    %     % recurse over all poses in a trajectory
    %     if ndims(T) == 3
    %         theta = zeros(size(T,3),robot.n);
    %         for k=1:size(T,3)
    %             theta(k,:) = ikine6s(robot, T(:,:,k), varargin{:});
    %         end
    %         return;
    %     end
    %
    %
    %     if ~ishomog(T)
    %         error('RTB:ikine:badarg', 'T is not a homog xform');
    %     end

    TT = SE3(TT);

    L = robot.links;

    if ~robot.isspherical()
        error('RTB:ikine:notsupported', 'wrist is not spherical');
    end

    % The configuration parameter determines what n1,n2,n4 values are used
    % and how many solutions are determined which have values of -1 or +1.

    if nargin < 3
        configuration = '';
    else
        configuration = lower(varargin{1});
    end

    % default configuration

    sol = [1 1 1];  % left, up, noflip

    for c=configuration
        switch c
            case 'l'
                sol(1) = 1;
            case 'r'
                sol(1) = 2;
            case 'u'
                sol(2) = 1;
            case 'd'
                sol(2) = 2;
            case 'n'
                sol(3) = 1;
            case 'f'
                sol(3) = 2;
        end
    end

    % determine the arm structure and the relevant solution to use
    if isempty(robot.ikineType)
        if is_simple(L)
            robot.ikineType = 'nooffset';
        elseif is_puma(L)
            robot.ikineType = 'puma';
        elseif is_offset(L)
            robot.ikineType = 'offset';
        elseif is_rrp(L)
            robot.ikineType = 'rrp';
        else
            error('RTB:ikine6s:badarg', 'This kinematic structure not supported');
        end
    end

    for k=1:length(TT)

        % undo base and tool transformations
        T = inv(robot.base) * TT(k) * inv(robot.tool);

        % drop back to matrix form
        T = T.T;

        %% now solve for the first 3 joints, based on position of the spherical wrist centre

        switch robot.ikineType
            case 'puma'
                % Puma model with shoulder and elbow offsets
                %
                % - Inverse kinematics for a PUMA 560,
                %   Paul and Zhang,
                %   The International Journal of Robotics Research,
                %   Vol. 5, No. 2, Summer 1986, p. 32-44
                %
                % Author::
                % Robert Biro with Gary Von McMurray,
                % GTRI/ATRP/IIMB,
                % Georgia Institute of Technology
                % 2/13/95

                a2 = L(2).a;
                a3 = L(3).a;
                d1 = L(1).d;
                d3 = L(3).d;
                d4 = L(4).d;

                % The following parameters are extracted from the Homogeneous
                % Transformation as defined in equation 1, p. 34

                Ox = T(1,2);
                Oy = T(2,2);
                Oz = T(3,2);

                Ax = T(1,3);
                Ay = T(2,3);
                Az = T(3,3);

                Px = T(1,4);
                Py = T(2,4);
                Pz = T(3,4) - d1;

                %
                % Solve for theta(1)
                %
                % r is defined in equation 38, p. 39.
                % theta(1) uses equations 40 and 41, p.39,
                % based on the configuration parameter n1
                %

                r = sqrt(Px^2 + Py^2);
                if sol(1) == 1
                    theta(1) = atan2(Py,Px) + pi - asin(d3/r);
                else
                    theta(1) = atan2(Py,Px) + asin(d3/r);
                end

                %
                % Solve for theta(2)
                %
                % V114 is defined in equation 43, p.39.
                % r is defined in equation 47, p.39.
                % Psi is defined in equation 49, p.40.
                % theta(2) uses equations 50 and 51, p.40, based on the configuration
                % parameter n2
                %
                if sol(2) == 1
                    n2 = -1;
                else
                    n2 = 1;
                end
                if sol(1) == 2
                    n2 = -n2;
                end

                V114 = Px*cos(theta(1)) + Py*sin(theta(1));
                r = sqrt(V114^2 + Pz^2);
                Psi = acos((a2^2-d4^2-a3^2+V114^2+Pz^2)/(2.0*a2*r));
                if ~isreal(Psi)
                    theta = [];
                else

                    theta(2) = atan2(Pz,V114) + n2*Psi;

                    %
                    % Solve for theta(3)
                    %
                    % theta(3) uses equation 57, p. 40.
                    %
                    num = cos(theta(2))*V114+sin(theta(2))*Pz-a2;
                    den = cos(theta(2))*Pz - sin(theta(2))*V114;
                    theta(3) = atan2(a3,d4) - atan2(num, den);
                end
            case 'nooffset'
                a2 = L(2).a;
                a3 = L(3).a;
                d1 = L(1).d;

                px = T(1,4); py = T(2,4); pz = T(3,4);

                %%% autogenerated code
                if L(1).alpha < 0
                    if sol(1) == 1
                        q(1) = angle(-px-py*1i);
                    else
                        q(1) = angle(px+py*1i);
                    end
                    S1 = sin(q(1));
                    C1 = cos(q(1));

                    if sol(2) == 1
                        q(2) = -angle(a2*d1*-2.0+a2*pz*2.0-C1*a2*px*2.0i-S1*a2*py*2.0i)+angle(d1*pz*2.0i-a2^2*1i+a3^2*1i-d1^2*1i-pz^2*1i-C1^2*px^2*1i-sqrt((a2*d1*2.0-a2*pz*2.0)^2+(C1*a2*px*2.0+S1*a2*py*2.0)^2-(d1*pz*-2.0+a2^2-a3^2+d1^2+pz^2+C1^2*px^2+S1^2*py^2+C1*S1*px*py*2.0)^2)-S1^2*py^2*1i-C1*S1*px*py*2.0i);
                    else
                        q(2) = -angle(a2*d1*-2.0+a2*pz*2.0-C1*a2*px*2.0i-S1*a2*py*2.0i)+angle(d1*pz*2.0i-a2^2*1i+a3^2*1i-d1^2*1i-pz^2*1i-C1^2*px^2*1i+sqrt((a2*d1*2.0-a2*pz*2.0)^2+(C1*a2*px*2.0+S1*a2*py*2.0)^2-(d1*pz*-2.0+a2^2-a3^2+d1^2+pz^2+C1^2*px^2+S1^2*py^2+C1*S1*px*py*2.0)^2)-S1^2*py^2*1i-C1*S1*px*py*2.0i);
                    end
                    S2 = sin(q(2));
                    C2 = cos(q(2));

                    if sol(3) == 1
                        q(3) = -angle(a2*-1i+C2*d1-C2*pz+S2*d1*1i-S2*pz*1i+C1*C2*px*1i-C1*S2*px+C2*S1*py*1i-S1*S2*py)+angle(a3*1i-sqrt((-a2+S2*d1-S2*pz+C1*C2*px+C2*S1*py)^2+(-C2*d1+C2*pz+C1*S2*px+S1*S2*py)^2-a3^2));
                    else
                        q(3) = -angle(a2*-1i+C2*d1-C2*pz+S2*d1*1i-S2*pz*1i+C1*C2*px*1i-C1*S2*px+C2*S1*py*1i-S1*S2*py)+angle(a3*1i+sqrt((-a2+S2*d1-S2*pz+C1*C2*px+C2*S1*py)^2+(-C2*d1+C2*pz+C1*S2*px+S1*S2*py)^2-a3^2));
                    end
                else
                    if sol(1) == 1
                        q(1) = angle(px+py*1i);
                    else
                        q(1) = angle(-px-py*1i);
                    end
                    S1 = sin(q(1));
                    C1 = cos(q(1));

                    if sol(2) == 1
                        q(2) = -angle(a2*d1*2.0-a2*pz*2.0-C1*a2*px*2.0i-S1*a2*py*2.0i)+angle(d1*pz*2.0i-a2^2*1i+a3^2*1i-d1^2*1i-pz^2*1i-C1^2*px^2*1i-sqrt((a2*d1*2.0-a2*pz*2.0)^2+(C1*a2*px*2.0+S1*a2*py*2.0)^2-(d1*pz*-2.0+a2^2-a3^2+d1^2+pz^2+C1^2*px^2+S1^2*py^2+C1*S1*px*py*2.0)^2)-S1^2*py^2*1i-C1*S1*px*py*2.0i);
                    else
                        q(2) = -angle(a2*d1*2.0-a2*pz*2.0-C1*a2*px*2.0i-S1*a2*py*2.0i)+angle(d1*pz*2.0i-a2^2*1i+a3^2*1i-d1^2*1i-pz^2*1i-C1^2*px^2*1i+sqrt((a2*d1*2.0-a2*pz*2.0)^2+(C1*a2*px*2.0+S1*a2*py*2.0)^2-(d1*pz*-2.0+a2^2-a3^2+d1^2+pz^2+C1^2*px^2+S1^2*py^2+C1*S1*px*py*2.0)^2)-S1^2*py^2*1i-C1*S1*px*py*2.0i);
                    end
                    S2 = sin(q(2));
                    C2 = cos(q(2));

                    if sol(3) == 1
                        q(3) = -angle(a2*-1i-C2*d1+C2*pz-S2*d1*1i+S2*pz*1i+C1*C2*px*1i-C1*S2*px+C2*S1*py*1i-S1*S2*py)+angle(a3*1i-sqrt((-a2-S2*d1+S2*pz+C1*C2*px+C2*S1*py)^2+(C2*d1-C2*pz+C1*S2*px+S1*S2*py)^2-a3^2));
                    else
                        q(3) = -angle(a2*-1i-C2*d1+C2*pz-S2*d1*1i+S2*pz*1i+C1*C2*px*1i-C1*S2*px+C2*S1*py*1i-S1*S2*py)+angle(a3*1i+sqrt((-a2-S2*d1+S2*pz+C1*C2*px+C2*S1*py)^2+(C2*d1-C2*pz+C1*S2*px+S1*S2*py)^2-a3^2));
                    end
                end

                theta(1:3) = q;

            case'offset'
                % general case with 6 length parameters
                a1 = L(1).a;
                a2 = L(2).a;
                a3 = L(3).a;
                d1 = L(1).d;
                d2 = L(2).d;
                d3 = L(3).d;

                px = T(1,4); py = T(2,4); pz = T(3,4);

                %%% autogenerated code
                if L(1).alpha < 0

                    if sol(1) == 1
                        q(1) = -angle(-px+py*1i)+angle(d2*1i+d3*1i-sqrt(d2*d3*-2.0-d2^2-d3^2+px^2+py^2));
                    else
                        q(1) = angle(d2*1i+d3*1i+sqrt(d2*d3*-2.0-d2^2-d3^2+px^2+py^2))-angle(-px+py*1i);
                    end
                    S1 = sin(q(1));
                    C1 = cos(q(1));

                    if sol(2) == 1
                        q(2) = angle(d1*pz*2.0i-sqrt(d1*pz^3*4.0+d1^3*pz*4.0-a1^4-a2^4-a3^4-d1^4-py^4-pz^4-C1^4*px^4+C1^2*py^4*2.0-C1^4*py^4+a1^2*a2^2*2.0+a1^2*a3^2*2.0+a2^2*a3^2*2.0-a1^2*d1^2*2.0+a2^2*d1^2*2.0+a3^2*d1^2*2.0-a1^2*py^2*6.0-a2^2*py^2*2.0+a3^2*py^2*2.0-a1^2*pz^2*2.0+a2^2*pz^2*2.0+a3^2*pz^2*2.0-d1^2*py^2*2.0-d1^2*pz^2*6.0-py^2*pz^2*2.0+d1*py^2*pz*4.0+C1^3*a1*px^3*4.0+S1^3*a1*py^3*4.0-C1^2*a1^2*px^2*6.0+C1^2*a2^2*px^2*2.0+C1^2*a3^2*px^2*2.0+C1^2*a1^2*py^2*6.0+C1^2*a2^2*py^2*1.0e1-C1^2*a3^2*py^2*2.0-C1^4*a2^2*py^2*1.2e1+C1^6*a2^2*py^2*4.0-C1^2*d1^2*px^2*2.0+C1^2*d1^2*py^2*2.0-C1^2*px^2*py^2*6.0+C1^4*px^2*py^2*6.0-C1^2*px^2*pz^2*2.0+C1^2*py^2*pz^2*2.0+S1^6*a2^2*py^2*4.0+C1*a1^3*px*4.0+S1*a1^3*py*4.0+a1^2*d1*pz*4.0-a2^2*d1*pz*4.0-a3^2*d1*pz*4.0-C1*a1*a2^2*px*4.0-C1*a1*a3^2*px*4.0-C1*S1*px^3*py*4.0+C1*a1*d1^2*px*4.0+C1*a1*px*py^2*1.2e1+C1*a1*px*pz^2*4.0+S1*a1*a2^2*py*4.0-S1*a1*a3^2*py*4.0+S1*a1*d1^2*py*4.0+S1*a1*px^2*py*4.0+S1*a1*py*pz^2*4.0-C1*S1^3*px*py^3*4.0+C1*S1^3*px^3*py*4.0-C1^3*a1*px*py^2*1.2e1-S1^3*a1*a2^2*py*8.0+C1^2*d1*px^2*pz*4.0-C1^2*d1*py^2*pz*4.0-S1^3*a1*px^2*py*4.0-C1*a1*d1*px*pz*8.0-S1*a1*d1*py*pz*8.0-C1*S1*a1^2*px*py*1.2e1-C1*S1*a2^2*px*py*4.0+C1*S1*a3^2*px*py*4.0-C1*S1*d1^2*px*py*4.0-C1*S1*px*py*pz^2*4.0-C1^2*S1*a1*a2^2*py*8.0+C1^2*S1*a1*px^2*py*8.0+C1*S1^3*a2^2*px*py*8.0+C1^3*S1*a2^2*px*py*8.0+C1*S1*d1*px*py*pz*8.0)-a1^2*1i-a2^2*1i+a3^2*1i-d1^2*1i-py^2*1i-pz^2*1i-C1^2*px^2*1i+C1^2*py^2*1i+C1*a1*px*2.0i+S1*a1*py*2.0i-C1*S1*px*py*2.0i)-angle(-a2*(a1*-1i+d1-pz+C1*px*1i+S1^3*py*1i+C1^2*S1*py*1i));
                    else
                        q(2) = angle(d1*pz*2.0i+sqrt(d1*pz^3*4.0+d1^3*pz*4.0-a1^4-a2^4-a3^4-d1^4-py^4-pz^4-C1^4*px^4+C1^2*py^4*2.0-C1^4*py^4+a1^2*a2^2*2.0+a1^2*a3^2*2.0+a2^2*a3^2*2.0-a1^2*d1^2*2.0+a2^2*d1^2*2.0+a3^2*d1^2*2.0-a1^2*py^2*6.0-a2^2*py^2*2.0+a3^2*py^2*2.0-a1^2*pz^2*2.0+a2^2*pz^2*2.0+a3^2*pz^2*2.0-d1^2*py^2*2.0-d1^2*pz^2*6.0-py^2*pz^2*2.0+d1*py^2*pz*4.0+C1^3*a1*px^3*4.0+S1^3*a1*py^3*4.0-C1^2*a1^2*px^2*6.0+C1^2*a2^2*px^2*2.0+C1^2*a3^2*px^2*2.0+C1^2*a1^2*py^2*6.0+C1^2*a2^2*py^2*1.0e1-C1^2*a3^2*py^2*2.0-C1^4*a2^2*py^2*1.2e1+C1^6*a2^2*py^2*4.0-C1^2*d1^2*px^2*2.0+C1^2*d1^2*py^2*2.0-C1^2*px^2*py^2*6.0+C1^4*px^2*py^2*6.0-C1^2*px^2*pz^2*2.0+C1^2*py^2*pz^2*2.0+S1^6*a2^2*py^2*4.0+C1*a1^3*px*4.0+S1*a1^3*py*4.0+a1^2*d1*pz*4.0-a2^2*d1*pz*4.0-a3^2*d1*pz*4.0-C1*a1*a2^2*px*4.0-C1*a1*a3^2*px*4.0-C1*S1*px^3*py*4.0+C1*a1*d1^2*px*4.0+C1*a1*px*py^2*1.2e1+C1*a1*px*pz^2*4.0+S1*a1*a2^2*py*4.0-S1*a1*a3^2*py*4.0+S1*a1*d1^2*py*4.0+S1*a1*px^2*py*4.0+S1*a1*py*pz^2*4.0-C1*S1^3*px*py^3*4.0+C1*S1^3*px^3*py*4.0-C1^3*a1*px*py^2*1.2e1-S1^3*a1*a2^2*py*8.0+C1^2*d1*px^2*pz*4.0-C1^2*d1*py^2*pz*4.0-S1^3*a1*px^2*py*4.0-C1*a1*d1*px*pz*8.0-S1*a1*d1*py*pz*8.0-C1*S1*a1^2*px*py*1.2e1-C1*S1*a2^2*px*py*4.0+C1*S1*a3^2*px*py*4.0-C1*S1*d1^2*px*py*4.0-C1*S1*px*py*pz^2*4.0-C1^2*S1*a1*a2^2*py*8.0+C1^2*S1*a1*px^2*py*8.0+C1*S1^3*a2^2*px*py*8.0+C1^3*S1*a2^2*px*py*8.0+C1*S1*d1*px*py*pz*8.0)-a1^2*1i-a2^2*1i+a3^2*1i-d1^2*1i-py^2*1i-pz^2*1i-C1^2*px^2*1i+C1^2*py^2*1i+C1*a1*px*2.0i+S1*a1*py*2.0i-C1*S1*px*py*2.0i)-angle(-a2*(a1*-1i+d1-pz+C1*px*1i+S1^3*py*1i+C1^2*S1*py*1i));
                    end
                    S2 = sin(q(2));
                    C2 = cos(q(2));

                    if sol(3) == 1
                        q(3) = angle(a3*1i-sqrt(d1*pz*-2.0+a1^2+a2^2-a3^2+d1^2+py^2+pz^2+C1^2*px^2-C1^2*py^2+C2*a1*a2*2.0-C1*a1*px*2.0-S2*a2*d1*2.0-S1*a1*py*2.0+S2*a2*pz*2.0-C1*C2*a2*px*2.0-C2*S1*a2*py*2.0+C1*S1*px*py*2.0))-angle(a2*-1i-C2*a1*1i+C2*d1-C2*pz+S2*a1+S2*d1*1i-S2*pz*1i+C1*C2*px*1i-C1*S2*px+C2*S1*py*1i-S1*S2*py);
                    else
                        q(3) = -angle(a2*-1i-C2*a1*1i+C2*d1-C2*pz+S2*a1+S2*d1*1i-S2*pz*1i+C1*C2*px*1i-C1*S2*px+C2*S1*py*1i-S1*S2*py)+angle(a3*1i+sqrt(d1*pz*-2.0+a1^2+a2^2-a3^2+d1^2+py^2+pz^2+C1^2*px^2-C1^2*py^2+C2*a1*a2*2.0-C1*a1*px*2.0-S2*a2*d1*2.0-S1*a1*py*2.0+S2*a2*pz*2.0-C1*C2*a2*px*2.0-C2*S1*a2*py*2.0+C1*S1*px*py*2.0));
                    end
                else
                    if sol(1) == 1
                        q(1) = -angle(px-py*1i)+angle(d2*1i+d3*1i-sqrt(d2*d3*-2.0-d2^2-d3^2+px^2+py^2));
                    else
                        q(1) = -angle(px-py*1i)+angle(d2*1i+d3*1i+sqrt(d2*d3*-2.0-d2^2-d3^2+px^2+py^2));
                    end
                    S1 = sin(q(1));
                    C1 = cos(q(1));

                    if sol(2) == 1
                        q(2) = angle(d1*pz*2.0i-sqrt(d1*pz^3*4.0+d1^3*pz*4.0-a1^4-a2^4-a3^4-d1^4-py^4-pz^4-C1^4*px^4+C1^2*py^4*2.0-C1^4*py^4+a1^2*a2^2*2.0+a1^2*a3^2*2.0+a2^2*a3^2*2.0-a1^2*d1^2*2.0+a2^2*d1^2*2.0+a3^2*d1^2*2.0-a1^2*py^2*6.0-a2^2*py^2*2.0+a3^2*py^2*2.0-a1^2*pz^2*2.0+a2^2*pz^2*2.0+a3^2*pz^2*2.0-d1^2*py^2*2.0-d1^2*pz^2*6.0-py^2*pz^2*2.0+d1*py^2*pz*4.0+C1^3*a1*px^3*4.0+S1^3*a1*py^3*4.0-C1^2*a1^2*px^2*6.0+C1^2*a2^2*px^2*2.0+C1^2*a3^2*px^2*2.0+C1^2*a1^2*py^2*6.0+C1^2*a2^2*py^2*1.0e1-C1^2*a3^2*py^2*2.0-C1^4*a2^2*py^2*1.2e1+C1^6*a2^2*py^2*4.0-C1^2*d1^2*px^2*2.0+C1^2*d1^2*py^2*2.0-C1^2*px^2*py^2*6.0+C1^4*px^2*py^2*6.0-C1^2*px^2*pz^2*2.0+C1^2*py^2*pz^2*2.0+S1^6*a2^2*py^2*4.0+C1*a1^3*px*4.0+S1*a1^3*py*4.0+a1^2*d1*pz*4.0-a2^2*d1*pz*4.0-a3^2*d1*pz*4.0-C1*a1*a2^2*px*4.0-C1*a1*a3^2*px*4.0-C1*S1*px^3*py*4.0+C1*a1*d1^2*px*4.0+C1*a1*px*py^2*1.2e1+C1*a1*px*pz^2*4.0+S1*a1*a2^2*py*4.0-S1*a1*a3^2*py*4.0+S1*a1*d1^2*py*4.0+S1*a1*px^2*py*4.0+S1*a1*py*pz^2*4.0-C1*S1^3*px*py^3*4.0+C1*S1^3*px^3*py*4.0-C1^3*a1*px*py^2*1.2e1-S1^3*a1*a2^2*py*8.0+C1^2*d1*px^2*pz*4.0-C1^2*d1*py^2*pz*4.0-S1^3*a1*px^2*py*4.0-C1*a1*d1*px*pz*8.0-S1*a1*d1*py*pz*8.0-C1*S1*a1^2*px*py*1.2e1-C1*S1*a2^2*px*py*4.0+C1*S1*a3^2*px*py*4.0-C1*S1*d1^2*px*py*4.0-C1*S1*px*py*pz^2*4.0-C1^2*S1*a1*a2^2*py*8.0+C1^2*S1*a1*px^2*py*8.0+C1*S1^3*a2^2*px*py*8.0+C1^3*S1*a2^2*px*py*8.0+C1*S1*d1*px*py*pz*8.0)-a1^2*1i-a2^2*1i+a3^2*1i-d1^2*1i-py^2*1i-pz^2*1i-C1^2*px^2*1i+C1^2*py^2*1i+C1*a1*px*2.0i+S1*a1*py*2.0i-C1*S1*px*py*2.0i)-angle(-a2*(a1*-1i-d1+pz+C1*px*1i+S1^3*py*1i+C1^2*S1*py*1i));
                    else
                        q(2) = angle(d1*pz*2.0i+sqrt(d1*pz^3*4.0+d1^3*pz*4.0-a1^4-a2^4-a3^4-d1^4-py^4-pz^4-C1^4*px^4+C1^2*py^4*2.0-C1^4*py^4+a1^2*a2^2*2.0+a1^2*a3^2*2.0+a2^2*a3^2*2.0-a1^2*d1^2*2.0+a2^2*d1^2*2.0+a3^2*d1^2*2.0-a1^2*py^2*6.0-a2^2*py^2*2.0+a3^2*py^2*2.0-a1^2*pz^2*2.0+a2^2*pz^2*2.0+a3^2*pz^2*2.0-d1^2*py^2*2.0-d1^2*pz^2*6.0-py^2*pz^2*2.0+d1*py^2*pz*4.0+C1^3*a1*px^3*4.0+S1^3*a1*py^3*4.0-C1^2*a1^2*px^2*6.0+C1^2*a2^2*px^2*2.0+C1^2*a3^2*px^2*2.0+C1^2*a1^2*py^2*6.0+C1^2*a2^2*py^2*1.0e1-C1^2*a3^2*py^2*2.0-C1^4*a2^2*py^2*1.2e1+C1^6*a2^2*py^2*4.0-C1^2*d1^2*px^2*2.0+C1^2*d1^2*py^2*2.0-C1^2*px^2*py^2*6.0+C1^4*px^2*py^2*6.0-C1^2*px^2*pz^2*2.0+C1^2*py^2*pz^2*2.0+S1^6*a2^2*py^2*4.0+C1*a1^3*px*4.0+S1*a1^3*py*4.0+a1^2*d1*pz*4.0-a2^2*d1*pz*4.0-a3^2*d1*pz*4.0-C1*a1*a2^2*px*4.0-C1*a1*a3^2*px*4.0-C1*S1*px^3*py*4.0+C1*a1*d1^2*px*4.0+C1*a1*px*py^2*1.2e1+C1*a1*px*pz^2*4.0+S1*a1*a2^2*py*4.0-S1*a1*a3^2*py*4.0+S1*a1*d1^2*py*4.0+S1*a1*px^2*py*4.0+S1*a1*py*pz^2*4.0-C1*S1^3*px*py^3*4.0+C1*S1^3*px^3*py*4.0-C1^3*a1*px*py^2*1.2e1-S1^3*a1*a2^2*py*8.0+C1^2*d1*px^2*pz*4.0-C1^2*d1*py^2*pz*4.0-S1^3*a1*px^2*py*4.0-C1*a1*d1*px*pz*8.0-S1*a1*d1*py*pz*8.0-C1*S1*a1^2*px*py*1.2e1-C1*S1*a2^2*px*py*4.0+C1*S1*a3^2*px*py*4.0-C1*S1*d1^2*px*py*4.0-C1*S1*px*py*pz^2*4.0-C1^2*S1*a1*a2^2*py*8.0+C1^2*S1*a1*px^2*py*8.0+C1*S1^3*a2^2*px*py*8.0+C1^3*S1*a2^2*px*py*8.0+C1*S1*d1*px*py*pz*8.0)-a1^2*1i-a2^2*1i+a3^2*1i-d1^2*1i-py^2*1i-pz^2*1i-C1^2*px^2*1i+C1^2*py^2*1i+C1*a1*px*2.0i+S1*a1*py*2.0i-C1*S1*px*py*2.0i)-angle(-a2*(a1*-1i-d1+pz+C1*px*1i+S1^3*py*1i+C1^2*S1*py*1i));
                    end
                    S2 = sin(q(2));
                    C2 = cos(q(2));

                    if sol(3) == 1
                        q(3) = angle(a3*1i-sqrt(d1*pz*-2.0+a1^2+a2^2-a3^2+d1^2+py^2+pz^2+C1^2*px^2-C1^2*py^2+C2*a1*a2*2.0-C1*a1*px*2.0+S2*a2*d1*2.0-S1*a1*py*2.0-S2*a2*pz*2.0-C1*C2*a2*px*2.0-C2*S1*a2*py*2.0+C1*S1*px*py*2.0))-angle(a2*-1i-C2*a1*1i-C2*d1+C2*pz+S2*a1-S2*d1*1i+S2*pz*1i+C1*C2*px*1i-C1*S2*px+C2*S1*py*1i-S1*S2*py);
                    else
                        q(3) = -angle(a2*-1i-C2*a1*1i-C2*d1+C2*pz+S2*a1-S2*d1*1i+S2*pz*1i+C1*C2*px*1i-C1*S2*px+C2*S1*py*1i-S1*S2*py)+angle(a3*1i+sqrt(d1*pz*-2.0+a1^2+a2^2-a3^2+d1^2+py^2+pz^2+C1^2*px^2-C1^2*py^2+C2*a1*a2*2.0-C1*a1*px*2.0+S2*a2*d1*2.0-S1*a1*py*2.0-S2*a2*pz*2.0-C1*C2*a2*px*2.0-C2*S1*a2*py*2.0+C1*S1*px*py*2.0));
                    end
                end

                theta(1:3) = q;

            case 'rrp'
                % RRP (Stanford arm like)

                px = T(1,4); py = T(2,4); pz = T(3,4);
                d1 = L(1).d;
                d2 = L(2).d;

                %%% autogenerated code
                if L(1).alpha < 0
                    if sol(1) == 1
                        q(1) = -angle(-px+py*1i)+angle(d2*1i-sqrt(-d2^2+px^2+py^2));
                    else
                        q(1) = angle(d2*1i+sqrt(-d2^2+px^2+py^2))-angle(-px+py*1i);
                    end
                    S1 = sin(q(1));
                    C1 = cos(q(1));

                    if sol(2) == 1
                        q(2) = angle(d1-pz-C1*px*1i-S1*py*1i);
                    else
                        q(2) = angle(-d1+pz+C1*px*1i+S1*py*1i);
                    end
                    S2 = sin(q(2));
                    C2 = cos(q(2));

                    q(3) = -C2*d1+C2*pz+C1*S2*px+S1*S2*py;

                else
                    if sol(1) == 1
                        q(1) = -angle(px-py*1i)+angle(d2*1i-sqrt(-d2^2+px^2+py^2));
                    else
                        q(1) = -angle(px-py*1i)+angle(d2*1i+sqrt(-d2^2+px^2+py^2));
                    end
                    S1 = sin(q(1));
                    C1 = cos(q(1));

                    if sol(2) == 1
                        q(2) = angle(-d1+pz-C1*px*1i-S1*py*1i);
                    else
                        q(2) = angle(d1-pz+C1*px*1i+S1*py*1i);
                    end
                    S2 = sin(q(2));
                    C2 = cos(q(2));

                    q(3) = -C2*d1+C2*pz-C1*S2*px-S1*S2*py;
                end
                theta(1:3) = q;

            case 'kr5'
                %Given function will calculate inverse kinematics for KUKA KR5 robot

                % Equations are calculated and implemented by
                % Gautam Sinha
                % Autobirdz Systems Pvt. Ltd.
                % SIDBI Office,
                % Indian Institute of Technology Kanpur, Kanpur, Uttar Pradesh
                % 208016
                % India
                %email- gautam.sinha705@gmail.com

                % get the a1, a2 and a3-- link lenghts for link no 1,2,3
                L = robot.links;
                a1 = L(1).a;
                a2 = L(2).a;
                a3 = L(3).a;

                % Check wether wrist is spherical or not
                if ~robot.isspherical()
                    error('wrist is not spherical')
                end

                % get d1,d2,d3,d4---- Link offsets for link no 1,2,3,4
                d1 = L(1).d;
                d2 = L(2).d;
                d3 = L(3).d;
                d4 = L(4).d;

                % Get the parameters from transformation matrix
                Ox = T(1,2);
                Oy = T(2,2);
                Oz = T(3,2);

                Ax = T(1,3);
                Ay = T(2,3);
                Az = T(3,3);

                Px = T(1,4);
                Py = T(2,4);
                Pz = T(3,4);

                % Set the parameters n1, n2 and n3 to get required configuration from
                % solution
                n1 = -1;   % 'l'
                n2 = -1;   % 'u'
                n4 = -1;   % 'n'
                if ~isempty(strfind(configuration, 'l'))
                    n1 = -1;
                end
                if ~isempty(strfind(configuration, 'r'))
                    n1 = 1;
                end
                if ~isempty(strfind(configuration, 'u'))
                    if n1 == 1
                        n2 = 1;
                    else
                        n2 = -1;
                    end
                end
                if ~isempty(strfind(configuration, 'd'))
                    if n1 == 1
                        n2 = -1;
                    else
                        n2 = 1;
                    end
                end
                if ~isempty(strfind(configuration, 'n'))
                    n4 = 1;
                end
                if ~isempty(strfind(configuration, 'f'))
                    n4 = -1;
                end

                % Calculation for theta(1)
                r=sqrt(Px^2+Py^2);

                if (n1 == 1)
                    theta(1)= atan2(Py,Px) + asin((d2-d3)/r);
                else
                    theta(1)= atan2(Py,Px)+ pi - asin((d2-d3)/r);
                end

                % Calculation for theta(2)
                X= Px*cos(theta(1)) + Py*sin(theta(1)) - a1;
                r=sqrt(X^2 + (Pz-d1)^2);
                Psi = acos((a2^2-d4^2-a3^2+X^2+(Pz-d1)^2)/(2.0*a2*r));

                if ~isreal(Psi)
                    warning('RTB:ikine6s:notreachable', 'point not reachable');
                    theta = [NaN NaN NaN NaN NaN NaN];
                    return
                end

                theta(2) = atan2((Pz-d1),X) + n2*Psi;

                % Calculation for theta(3)
                Nu = cos(theta(2))*X + sin(theta(2))*(Pz-d1) - a2;
                Du = sin(theta(2))*X - cos(theta(2))*(Pz-d1);
                theta(3) = atan2(a3,d4) - atan2(Nu, Du);

                % Calculation for theta(4)
                Y = cos(theta(1))*Ax + sin(theta(1))*Ay;
                M2 = sin(theta(1))*Ax - cos(theta(1))*Ay ;
                M1 =  ( cos(theta(2)-theta(3)) )*Y + ( sin(theta(2)-theta(3)) )*Az;
                theta(4) = atan2(n4*M2,n4*M1);

                % Calculation for theta(5)
                Nu =  -cos(theta(4))*M1 - M2*sin(theta(4));
                M3 =  -Az*( cos(theta(2)-theta(3)) ) + Y*( sin(theta(2)-theta(3)) );
                theta(5) = atan2(Nu,M3);

                % Calculation for theta(6)
                Z = cos(theta(1))*Ox + sin(theta(1))*Oy;
                L2 = sin(theta(1))*Ox - cos(theta(1))*Oy;
                L1 = Z*( cos(theta(2)-theta(3) )) + Oz*( sin(theta(2)-theta(3)));
                L3 = Z*( sin(theta(2)-theta(3) )) - Oz*( cos(theta(2)-theta(3)));
                A1 = L1*cos(theta(4)) + L2*sin(theta(4));
                A3 = L1*sin(theta(4)) - L2*cos(theta(4));
                Nu =  -A1*cos(theta(5)) - L3*sin(theta(5));
                Du =  -A3;
                theta(6) = atan2(Nu,Du);

            otherwise
                error('RTB:ikine6s:badarg', 'Unknown solution type [%s]', robot.ikineType);
        end

        if ~isempty(theta)
            % Solve for the wrist rotation

            % we need to account for some random translations between the first and last 3
            % joints (d4) and also d6,a6,alpha6 in the final frame.

            T13 = robot.A(1:3, theta(1:3));  % transform of first 3 joints

            % T = T13 * Tz(d4) * R * Tz(d6) Tx(a5)
            Td4 = SE3(0, 0, L(4).d);      % Tz(d4)
            Tt = SE3(L(6).a, 0, L(6).d) * SE3.Rx(L(6).alpha);  % Tz(d6) Tx(a5) Rx(alpha6)

            R = inv(Td4) * inv(T13) * SE3(T) * inv(Tt);

            % the spherical wrist implements Euler angles

            if sol(3) == 1
                theta(4:6) = tr2eul(R, 'flip');
            else
                theta(4:6) = tr2eul(R);

            end
            if L(4).alpha > 0
                theta(5) = -theta(5);
            end

            % remove the link offset angles
            for j=1:robot.n   %#ok<*AGROW>
                theta(j) = theta(j) - L(j).offset;
            end

            % stack the rows
            thetavec(k,:) = theta;
        else
            warning('RTB:ikine6s:notreachable', 'point not reachable');
            thetavec(k,:) = [NaN NaN NaN NaN NaN NaN];
        end
    end
end

% predicates to determine which kinematic solution to use
function s = is_simple(L)
    alpha = [-pi/2 0 pi/2];
    s =     all([L(2:3).d] == 0) && ...
        (all([L(1:3).alpha] == alpha) || all([L(1:3).alpha] == -alpha)) && ...
        all([L(1:3).isrevolute] == 1) && ...
        (L(1).a == 0);
end

function s = is_offset(L)
    alpha = [-pi/2 0 pi/2];
    s =     (all([L(1:3).alpha] == alpha) || all([L(1:3).alpha] == -alpha)) && ...
        all([L(1:3).isrevolute] == 1);
end

function s = is_rrp(L)
    alpha = [-pi/2 pi/2 0];
    s =     all([L(2:3).a] == 0) && ...
        (all([L(1:3).alpha] == alpha) || all([L(1:3).alpha] == -alpha)) && ...
        all([L(1:3).isrevolute] == [1 1 0]);
end

function s = is_puma(L)
    alpha = [pi/2 0 -pi/2];
    s =     (L(2).d == 0) && (L(1).a == 0) && ...
        (L(3).d ~= 0) && (L(3).a ~= 0) && ...
        all([L(1:3).alpha] == alpha) && ...
        all([L(1:3).isrevolute] == 1);
end


text:
function siz = size(a,idx)
%SIZE Return size of sptenmat.
%
%   D = SIZE(T) returns the size of the tensor.
%
%   I = size(T,DIM) returns the sizes of the dimensions specified by DIM.
%
%   See also SPTENMAT.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

if isempty(a.tsize)
    siz = [];
    return;
end

m = prod(a.tsize(a.rdims));
n = prod(a.tsize(a.cdims));
siz = [m n];

if  exist('idx','var')
    siz = siz(idx);
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
%SHOWFIGS Show all figures on the screen
%
%  SHOWFIGS(K)
%
% Use K figures on a row

function showfigs(k)

h = sort(get(0,'children'));  % handles for all figures
n = length(h);                % number of figure
if nargin == 0
	k = ceil(sqrt(n));          % figures to be shown
end
s = 0.95/k;   % screen stitch
r = 0.93;     % image size reduction
t = 0.055;    % top gap
b = 0.005;    % border gap
fig = 0;
set(0,'units','pixels');
monpos = get(0,'monitorposition');
monpos = monpos(1,:);
for i=1:k
	for j=1:k
		fig = fig+1;
		if fig > n, break; end
		set(h(fig),'units','pixels','position',[(j-1)*s+b,(1-t)-i*s,s*r,s*r]*monpos(4));
		figure(h(fig));
	end
end
for j=n:-1:1, figure(h(j)); end

author: marianux

repo: ecg-kit

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


text:
% Usage: compute_rmse( Yrec, dataset, remove_val )
%        compute_rmse( name_of_file_with_reconstructions )
function [ out, lat, lon, time ] = metoffice_compute_rmse( Yrec, dataset, remove_val )

if isstr(Yrec)
    setup = get_setup( Yrec )

    % Load reconstructions
    datapath = metoffice_get_path();
    load( [ datapath, Yrec ], 'Yrec' );
    %load( [ datapath '/RecTest/RESULTS/' Yrec ], 'Yrec' );
else
    setup.dataset = dataset;
    setup.remove_val = remove_val;
    setup
end

% Remove values below -1.8
Yrec( Yrec < -1.8 ) = -1.8;

[ true_sst, Itest, Itrain, time, train_sst ] = metoffice_get_testdata( setup.dataset );
%[ true_sst, Itest, Itrain, time, train_sst ] = load_testdata( setup.dataset );

if setup.remove_val
    Ival = get_valset( setup.remove_val, time, Itrain );
    Itrain = Itrain & ~Ival;
end

% Load lat and lon to compute the area weights
datapath = metoffice_get_path();
load( [ datapath, dataset ], 'lat', 'lon' )
%load( [ datapath '/RecTest/DATASETS/' dataset ], 'lat', 'lon' )
lts = repmat( lat, 1, length(lon) );
weights = cosd( lts(:) );
clear lts lat lon

err = true_sst - Yrec;
out.train = get_rmse( err, Itrain, weights );
out.test = get_rmse( err, Itest, weights );
if setup.remove_val
    out.val = get_rmse( err, Ival, weights );
end

%out.Yrec = Yrec;

out.trdata = get_rmse( train_sst-Yrec, Itrain, weights );

%out.test.err = err;
%out.test.err(~Itest) = NaN;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = get_rmse( err, I, weights )

err2 = err.^2;
err2 = bsxfun( @times, weights, err2 );
Iw = bsxfun( @times, weights, I );

err2(~I) = 0;
out.rmse_fld = sqrt( sum(err2,2)./sum(I,2) );

out.rmse_ts = sqrt( sum(err2,1)./sum(I,1) );

out.rmse = sqrt( sum( err2(:) )/sum(I(:)) );

if 1
    out.rmse2 = sum(err2,1)./sum(Iw,1);
    out.rmse2( isnan(out.rmse2) ) = [];
    out.rmse2 = sqrt( mean( out.rmse2 ) );
end

out.err_fld = sqrt(sum(err2,2));
out.err_ts = sqrt(sum(err2,1));

% Smoothed time series: take 5 year moving average
err2_ts = sum(err2,1);
nobs_ts = sum(I,1);

nfilt = 5*12;
T = size(err2,2);
err2_sm = conv( ones(1,nfilt), err2_ts );
err2_sm = err2_sm( nfilt/2+(1:T) );
nobs_sm = conv( ones(1,nfilt), nobs_ts );
nobs_sm = nobs_sm( nfilt/2+(1:T) );

out.rmse_smts = sqrt( err2_sm./nobs_sm );

out.err_smts = sqrt(err2_sm);

% Median
out.medwe = median( err(I) );
out.medwe_fld = NaN*zeros(size(err,1),1);
for i = 1:size(err,1)
    out.medwe_fld(i) = median( err(i, I(i,:) ) );
end

out.medwe_ts = NaN*zeros(1,size(err,2));
for i = 1:size(err,2)
    out.medwe_ts(i) = median( err( I(:,i), i ) );
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
%kvspatial 'Compute Spatial Features Using NxM Window (K1)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vspatial.pane file
%
% Parameters:
% InputFile: i 'Input Image ', required: 'input image filename'
% Integer: w 'Window width ', default: 3: 'window width argument'
% Integer: h 'Window height', default: 3: 'window height argument'
% OutputFile: o 'Output Image', required: 'output image filename'
%
% Example: o = kvspatial(i, {'i','';'w',3;'h',3;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vspatial - Compute Spatial Features Using NxM Window  (K1)
%
%  DESCRIPTION
% .I vspatial
% performs a spatial feature extraction on an image using a selected
% statistical operation.  The statistical operation is computed, using
% overlapping windows, over the entire image.  The center pixel is updated
% in each window based on the calculated statistics of that window.
%
% The input image, which must be of data storage type BYTE, corresponds
% to the -i argument.  The output image, which corresponds to the -o
% argument, is stored as data type FLOAT.  The output image may be a
% multiband image depending on how many spatial operators were selected.
% The resulting image will have a border of pixels of value 0, since the
% window operation does not calculate values for border pixels.
%
% The width and height of the window is specified by the -w and -h window
% width and height arguments.  The width and height arguments cannot be an
% even number.  The window dimensions must be odd (ie. 3x3, 5x5, 7x7, etc.),
% since the center pixel is replaced with the computed value of the window
% data.  The default window size is 3x3.
%
% Each pixel in the image is updated with a computed value except the border
% pixels.  A 3x3 window will result in a border of one pixel, a 5x5 window
% will result in a border of two pixels, etc.
%
% There are six statistical operations that may be selected.  Any combination
% of the statistical operations may be selected by toggling the desired
% argument on (ie. argument followed by a 1).  The default for each argument
% is unselected (ie. argument followed by a 0), except -m (mean), which is
% selected.  The six possible statistical operations and their arguments are:
% .DS
%   -m  --  calculates the mean of the image.
%   -v  --  calculates the variance of the image.
%   -c  --  calculates the contrast of the image.
%   -s  --  calculates the angular second moment of the image.
%   -e  --  calculates the entropy of the image.
%   -d  --  calculates the dispersion of the image.
% .DE
% The resulting feature vector may contain from one to six bands of data,
% depending upon the number of arguments selected.
%
% All input images must be of data storage type BYTE.
%
% All output images are of data storage type FLOAT.
%
%
%
%  EXAMPLES
% vspatial -i input_image.xv -w 5 -h 5 -m 0 -v 1 -e 1 -o output_image
%
% the input image must of type BYTE; the selected width and height of the
% window is 5 x 5, which will result in a border width of two pixels in the
% output image.  The mean spatial operator is unselected (default is selected),
% and the variance and entropy operations are selected, which will result in
% an output image (feature vector) of 2 bands, based on the operations selected.
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
% All input images must be of data storage type BYTE.
%
% All output images are of data storage type FLOAT.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kvspatial(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kvspatial(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'w', 3;'h', 3;'o', '__output'};
maxval={0,64,64,0};
minval={0,1,1,0};
istoggle=[0,1,1,0];
was_set=istoggle * 0;
paramtype={'InputFile','Integer','Integer','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'vspatial"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function prediction = BasicCifarNet(varargin)
%BASICCIFARNET Returns a simple network for CIFAR10
%   M = models.BasicCifarNet() returns a model inspired by:
%
%     LeCun et al., "Gradient-Based Learning Applied to Document
%     Recognition", Proceedings of the IEEE, 1998.
%
%   models.BasicCifarNet(..., 'option', value, ...) accepts the following
%   options:
%
%   `input`:: default input
%     Specifies an input (images) layer for the network. If unspecified, a
%     new one is created.
%
%   `numClasses`:: 10
%     Number of output classes.
%
%   `batchNorm`:: true
%     Whether to use batch normalization.
%
%   Any other options will be passed to models.ConvBlock(), and can be used
%   to change the activation function, weight initialization, etc.
%
%   Suggested SGD training options are also returned in the struct M.meta.

  % parse options. unknown arguments will be passed to ConvBlock (e.g.
  % activation).
  opts.input = Input('name', 'images', 'gpu', true) ;  % default input layer
  opts.numClasses = 10 ;
  opts.batchNorm = true ;  % whether to use batch normalization
  [opts, convBlockArgs] = vl_argparse(opts, varargin, 'nonrecursive') ;

  % get conv block generator with the given options. default activation is
  % ReLU, with pre-activation batch normalization (can be overriden).
  conv = models.ConvBlock('batchNorm', opts.batchNorm, ...
    'preActivationBatchNorm', true, convBlockArgs{:}) ;

  % build network
  images = opts.input ;

  x = conv(images, 'size', [5, 5, 3, 32], 'pad', 2, 'weightScale', 0.01) ;
  x = vl_nnpool(x, 3, 'stride', 2, 'method', 'max', 'pad', 1) ;

  x = conv(x, 'size', [5, 5, 32, 32], 'pad', 2, 'weightScale', 0.05) ;
  x = vl_nnpool(x, 3, 'stride', 2, 'method', 'avg', 'pad', 1) ;

  x = conv(x, 'size', [5, 5, 32, 64], 'pad', 2, 'weightScale', 0.05) ;
  x = vl_nnpool(x, 3, 'stride', 2, 'method', 'avg', 'pad', 1) ;

  x = conv(x, 'size', [4, 4, 64, 64], 'weightScale', 0.05) ;

  prediction = conv(x, 'size', [1, 1, 64, opts.numClasses], 'weightScale', 0.05, ...
    'batchNorm', false, 'activation', 'none') ;

  % default training options for this network
  defaults.batchSize = 128 ;
  defaults.weightDecay = 0.0005 ;
  if ~opts.batchNorm
    defaults.learningRate = 0.01 ;
    defaults.numEpochs = 100 ;
  else
    defaults.learningRate = 0.1 ;
    defaults.numEpochs = 40 ;
  end
  defaults.imageSize = [32, 32, 3] ;
  prediction.meta = defaults ;

end

repo: MFT


text:
function f = addtheme(f, theme)

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

    %-validate theme name-%

    themePath = 'plotly/themes';
    S = dir(themePath);

    if isempty(S)
        paths = split(path, ':');

        for p = 1:length(paths)
            if ~isempty(strfind(paths{p}, themePath))
                themePath = paths{p};
                break;
            end
        end

        S = dir(themePath);
    end

    N = {S.name};

    if ~any(strcmp(N,strcat(theme, '.json'))) == 1
        ME = MException('MyComponent:noSuchVariable',...
            [strcat('\n<strong>', theme,...
            '</strong> is not a supported themes.'),...
            ' Please choose one of these theme names:\n\n',...
            strrep(strrep([S.name], '...', ''), '.json', ' | ')]);
        throw(ME)
    end

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

    %-add theme to figure-%

    fname = sprintf('%s/%s.json', themePath, theme);
    % fname = strcat('plotly/themes/', theme, '.json');
    fid = fopen(fname);
    raw = fread(fid,inf);
    str = char(raw');
    fclose(fid);
    theme_template = jsondecode(str);

    f.layout.template = theme_template;

    if isfield(f.layout.template.layout, 'paper_bgcolor')
        disp(strcat('layout.bg_paper:::',...
            f.layout.template.layout.paper_bgcolor))
    end

    if isfield(f.layout.template.layout, 'plot_bgcolor')
        disp(strcat('layout.plot_bgcolor:::',...
            f.layout.template.layout.plot_bgcolor))
    end

    %-------------------------------------------------------------------------%
end

author: plotly


text:
clc;

S = [22  30];
bb = [1 1 15 15]';

grid = bb_scan(bb,S)
%  grid = grid(:,1:1000);
%  o = bb_overlap(grid,grid,2);

repo: OpenTLD


text:
%% mergeImageData
% Below is a demonstration of the features of the |mergeImageData| function

%%
clear; close all; clc;

%% Syntax
% |[N12,G12]=mergeImageData(N1,G1,N2,G2,nFac,plotOn);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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 3 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, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [z,pv1,pv2,pv3]=compare_bino_prob(X,Y)
% [z,pv1,pv2,pv3]=compare_bino_prob(X,Y)
% Compares probabilities of two binomial samples X, Y
% Let X and Y be sets of 0 or 1, results from a Bernoully experiments
% with probability p1 for set X and p2 for set Y.
% Test the null hypothesis
%       H_0: p1 = p2
%    vs alternatives
%       H_1: p1 > p2
%       H_2: p1 < p2
%       H_3: p1 != p2
%
% Input:
%       X - row of 0 or 1
%       Y - row of 0 or 1
%
% Output:
%       z   - normal statistics N(0,1) for H_0
%       pv1 - p-value for H_1
%       pv2 - p-value for H_2
%       pv3 - p-value for H_3

% Dimiar Atanasov (2008)
% datanasov@nbu.bg

if size(X,1) > 1 || size(Y,1) > 1
    error('Sets should be rows');
end;

n_x1 = size( find(X == 1), 1 );
n_x0 = size( find(X == 0), 1 );

n_y1 = size( find(Y == 1), 2);
n_y0 = size( find(Y == 0), 2);

n_x = n_x1 + n_x0;
n_y = n_y1 + n_y0;

n_1 = n_x1 + n_y1;
n_0 = n_x0 + n_y0;

n = n_0 + n_1;

h_x = n_x1 / n_x;
h_y = n_y1 / n_y;

h = (n_x1 + n_y1)/(n_x + n_y);

if (n_1^2  / n < 5) || (n_0^2 / n < 5) || ( n_1*n_0 / n < 5)
    disp('Missing asymptotic behaviour!!!');
end;

s = h*(1-h)*(1/n_x + 1/n_y);

z = (h_x - h_y) / sqrt(s);

pv1 = 1 - normcdf(z);
pv2 = normcdf(z);
pv3 = (1 - normcdf( abs(z) ))/2;

author: Sable

repo: mcbench-benchmarks


text:
function t_opf_dc_cplex(quiet)
%T_OPF_DC_CPLEX  Tests for DC optimal power flow using CPLEX solver.

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

if nargin < 1
    quiet = 0;
end

algs = [1; 2; 3; 4; 5; 6];
alg_names = {
    'primal simplex',
    'dual simplex',
    'network simplex',
    'barrier',
    'sifting',
    'concurrent'
};
num_tests = 43 * length(algs);

t_begin(num_tests, quiet);

[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;

casefile = 't_case9_opf';
if quiet
    verbose = 0;
else
    verbose = 0;
end

mpopt = mpoption('out.all', 0, 'verbose', verbose);
mpopt = mpoption(mpopt, 'opf.dc.solver', 'CPLEX');

%% run DC OPF
if have_feature('cplex')
  for k = 1:length(algs)
    mpopt = mpoption(mpopt, 'cplex.lpmethod', algs(k), 'cplex.qpmethod', algs(k));
    t0 = sprintf('DC OPF (CPLEX %s): ', alg_names{k});

    %% set up indices
    ib_data     = [1:BUS_AREA BASE_KV:VMIN];
    ib_voltage  = [VM VA];
    ib_lam      = [LAM_P LAM_Q];
    ib_mu       = [MU_VMAX MU_VMIN];
    ig_data     = [GEN_BUS QMAX QMIN MBASE:APF];
    ig_disp     = [PG QG VG];
    ig_mu       = (MU_PMAX:MU_QMIN);
    ibr_data    = (1:ANGMAX);
    ibr_flow    = (PF:QT);
    ibr_mu      = [MU_SF MU_ST];
    ibr_angmu   = [MU_ANGMIN MU_ANGMAX];

    %% get solved DC power flow case from MAT-file
    load soln9_dcopf;       %% defines bus_soln, gen_soln, branch_soln, f_soln

    %% run OPF
    t = t0;
    [baseMVA, bus, gen, gencost, branch, f, success, et] = rundcopf(casefile, mpopt);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %%-----  test OPF with angle difference limits  -----
    t = [t0 'w/angle diff lims : '];
    mpc = loadcase(casefile);
    mpc.branch(4, ANGMAX) = 3;
    mpc.branch(7, ANGMIN) = -4.5;
    r = rundcopf(mpc, mpopt);
    [bus, gen, branch, f, success] = deal(r.bus, r.gen, r.branch, r.f, r.success);
    t_ok(success, [t 'success']);
    t_is(   f, 6456.7213, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,PG        ),    [99.98497;89.35133;125.66371], 4, [t 'gen dispatch']);
    t_is(branch(:,ibr_data  ), mpc.branch(:,ibr_data   ), 10, [t 'branch data']);
    e = zeros(size(branch, 1), 1);
    e(4) = 297.83776;
    e(7) = -26.94788;
    t_is(branch(:,MU_ANGMAX )-branch(:,MU_ANGMIN ), e, 4, [t 'branch ang diff mu']);

    t = [t0 'w/ignored angle diff lims : '];
    mpopt1 = mpoption(mpopt, 'opf.ignore_angle_lim', 1);
    r = rundcopf(mpc, mpopt1);
    [bus, gen, branch, f, success] = deal(r.bus, r.gen, r.branch, r.f, r.success);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), mpc.branch(:,ibr_data   ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %%-----  run OPF with extra linear user constraints & costs  -----
    %% two new z variables
    %%      0 <= z1, P3 - P1 <= z1
    %%      0 <= z2, P3 - P2 <= z2
    %% with A and N sized for DC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1;1;2;2;2],[10;12;13;12;11;14],[-1;1;-1;1;-1;-1],2,14);
    mpc.u = [0; 0];
    mpc.l = [-Inf; -Inf];
    mpc.zl = [0; 0];

    mpc.N = sparse([1;2], [13;14], [1;1], 2, 14);   %% new z variables only
    mpc.fparm = ones(2,1) * [1 0 0 1];              %% w = r = z
    mpc.H = sparse(2,2);                            %% no quadratic term
    mpc.Cw = [1000;1];

    t = [t0 'w/extra constraints & costs 1 : '];
    [r, success] = rundcopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    t_is(r.gen(1, PG), 116.15974, 5, [t 'Pg1 = 116.15974']);
    t_is(r.gen(3, PG), 116.15974, 5, [t 'Pg3 = 116.15974']);
    t_is(r.var.val.z, [0; 0.3348], 4, [t 'user vars']);
    t_is(r.cost.usr, 0.3348, 4, [t 'user costs']);

    %% with A and N sized for AC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1;1;2;2;2],[19;21;25;21;20;26],[-1;1;-1;1;-1;-1],2,26);
    mpc.u = [0; 0];
    mpc.l = [-Inf; -Inf];
    mpc.zl = [0; 0];

    mpc.N = sparse([1;2], [25;26], [1;1], 2, 26);   %% new z variables only
    mpc.fparm = ones(2,1) * [1 0 0 1];              %% w = r = z
    mpc.H = sparse(2,2);                            %% no quadratic term
    mpc.Cw = [1000;1];

    t = [t0 'w/extra constraints & costs 2 : '];
    [r, success] = rundcopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    t_is(r.gen(1, PG), 116.15974, 5, [t 'Pg1 = 116.15974']);
    t_is(r.gen(3, PG), 116.15974, 5, [t 'Pg3 = 116.15974']);
    t_is(r.var.val.z, [0; 0.3348], 4, [t 'user vars']);
    t_is(r.cost.usr, 0.3348, 4, [t 'user costs']);

    t = [t0 'infeasible : '];
    %% with A and N sized for DC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1], [10;11], [1;1], 1, 14);   %% Pg1 + Pg2
    mpc.u = Inf;
    mpc.l = 600;
    [r, success] = rundcopf(mpc, mpopt);
    t_ok(~success, [t 'no success']);

    %% OPF with all buses isolated
    t = [t0 'all buses isolated : '];
    mpc = loadcase(casefile);
    mpc.bus(:, BUS_TYPE) = NONE;
    try
        r = rundcopf(mpc, mpopt);
        t_is(r.success, 0, 12, [t 'success = 0']);
    catch
        t_ok(0, [t 'unexpected fatal error']);
    end
  end
else
    t_skip(num_tests, 'CPLEX not available');
end

t_end;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function rawdata = ReadRawDataFile(datafile)

metadatafns = {'mouse','exp','id','upper'};
metadataisnum = [true,false,false,false];

datafnsspecial = {'auto_Grab_success','auto_Grab_successtype'};

fid = fopen(datafile,'r');

% read headers
l = fgetl(fid);
headers = regexp(l,',','split');
idxdata = find(~cellfun(@isempty,regexp(headers,'^auto','once')) & ...
  ~ismember(headers,datafnsspecial));% & ...
  %cellfun(@isempty,regexp(headers,'successfultrial')));
metadatai = nan(1,numel(metadatafns));
for i = 1:numel(metadatafns),
  metadatai(i) = find(strcmpi(headers,metadatafns{i}));
end

successi = find(strcmpi(headers,'auto_Grab_success'),1);
successtypei = find(strcmpi(headers,'auto_Grab_successtype'),1);

datei = find(strcmpi(headers,'date'),1);
parentdiri = find(strcmpi(headers,'upper'),1);
expi = find(strcmpi(headers,'exp'),1);

datafns = headers(idxdata);

rawdata = [];

while true,
  l = fgetl(fid);
  if ~ischar(l),
    break;
  end
  if isempty(l),
    continue;
  end
  s = regexp(l,',','split');

  rawdatacurr = struct;

  m = regexp(s{datei},'^(\d+)([^\d]*)$','once','tokens');
  rawdatacurr.date = m{1};
  rawdatacurr.iscno = strcmpi(m{2},'CNO');
  m = regexp(s{parentdiri},'[\\/]([^\\/]+)$','tokens','once');
  rawdatacurr.session = m{1};
  m = regexp(s{expi},'_[vV](\d+)$','tokens','once');
  rawdatacurr.trial = str2double(m{1});

  for i = 1:numel(metadatafns),
    fn = metadatafns{i};
    rawdatacurr.(fn) = s{metadatai(i)};
    if metadataisnum(i),
      rawdatacurr.(fn) = str2double(rawdatacurr.(fn));
    end
  end

  rawdatacurr.success = find(strcmpi(s{successi},{'false','true'}))-1;
  rawdatacurr.successtype = s{successtypei};

  for i = 1:numel(datafns),
    fn = datafns{i};
    rawdatacurr.(fn) = str2double(s{idxdata(i)});
  end

  if isfield(rawdatacurr,'exp') && ~isfield(rawdatacurr,'expfull'),
    rawdatacurr.expfull = rawdatacurr.exp;
  end

  rawdata = structappend(rawdata,rawdatacurr);

end

fclose(fid);

author: kristinbranson

repo: JAABA


text:
%% Copyright (C) 2014-2017, 2019 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defop  Method   @@sym {vertcat} {(@var{x}, @var{y}, @dots{})}
%% @defopx Operator @@sym {[@var{x}; @var{y}; @dots{}]} {}
%% Vertically concatentate symbolic arrays.
%%
%% Example:
%% @example
%% @group
%% A = sym([1 2])
%%   @result{} A = (sym) [1  2]  (1×2 matrix)
%% [A; A; 2*A]
%%   @result{} (sym 3×2 matrix)
%%       ⎡1  2⎤
%%       ⎢    ⎥
%%       ⎢1  2⎥
%%       ⎢    ⎥
%%       ⎣2  4⎦
%% @end group
%% @end example
%% @seealso{@@sym/horzcat, @@sym/cat}
%% @end defop

function h = vertcat(varargin)

  % special case for 0x0 but other empties should be checked for
  % compatibilty
  cmd = {
          '_proc = []'
          'for i in _ins:'
          '    if i is None or not i.is_Matrix:'
          '        _proc.append(sp.Matrix([[i]]))'
          '    else:'
          '        if i.shape == (0, 0):'
          '            pass'
          '        else:'
          '            _proc.append(i)'
          'return sp.MatrixBase.vstack(*_proc),'
          };

  for i = 1:nargin
    varargin{i} = sym(varargin{i});
  end
  h = pycall_sympy__ (cmd, varargin{:});

end

%!test
%! % basic
%! syms x
%! A = [x; x];
%! B = vertcat(x, x);
%! C = vertcat(x, x, x);
%! assert (isa (A, 'sym'))
%! assert (isa (B, 'sym'))
%! assert (isa (C, 'sym'))
%! assert (isequal (size(A), [2 1]))
%! assert (isequal (size(B), [2 1]))
%! assert (isequal (size(C), [3 1]))

%!test
%! % basic, part 2
%! syms x
%! A = [x; 1];
%! B = [1; x];
%! C = [1; 2; x];
%! assert (isa (A, 'sym'))
%! assert (isa (B, 'sym'))
%! assert (isa (C, 'sym'))
%! assert (isequal (size(A), [2 1]))
%! assert (isequal (size(B), [2 1]))
%! assert (isequal (size(C), [3 1]))

%!test
%! % column vectors
%! a = [sym(1); 2];
%! b = [sym(3); 4];
%! assert (isequal ( [a;b] , [1; 2; 3; 4]  ))
%! assert (isequal ( [a;b;a] , [1; 2; 3; 4; 1; 2]  ))

%!test
%! % row vectors
%! a = [sym(1) 2];
%! b = [sym(3) 4];
%! assert (isequal ( [a;b] , [1 2; 3 4]  ))
%! assert (isequal ( [a;b;a] , [1 2; 3 4; 1 2]  ))

%!test
%! % row vector, other row
%! a = [sym(1) 2];
%! assert (isequal ( [a; [sym(3) 4]] , [1 2; 3 4]  ))

%!test
%! % empty vectors
%! v = [sym(1) sym(2)];
%! a = [v; []];
%! assert (isequal (a, v))
%! a = [[]; v; []];
%! assert (isequal (a, v))
%! a = [v; []; []];
%! assert (isequal (a, v))

%!xtest
%! % FIXME: is this Octave bug? worth worrying about
%! syms x
%! a = [x; [] []];
%! assert (isequal (a, x))

%!test
%! % more empty vectors
%! v = [sym(1) sym(2)];
%! q = sym(ones(0, 2));
%! assert (isequal ([v; q], v))

%!error <ShapeError>
%! v = [sym(1) sym(2)];
%! q = sym(ones(0, 3));
%! w = vertcat(v, q);

%!test
%! % Octave 3.6 bug: should pass on 3.8.1 and matlab
%! a = [sym(1) 2];
%! assert (isequal ( [a; [3 4]] , [1 2; 3 4]  ))
%! assert (isequal ( [a; sym(3) 4] , [1 2; 3 4]  ))
%! % more examples
%! syms x
%! [x [x x]; x x x];
%! [[x x] x; x x x];
%! [[x x] x; [x x] x];
%! [x x x; [x x] x];

%!test
%! % issue #700
%! A = sym ([1 2]);
%! B = simplify (A);
%! assert (isequal ([B; A], [A; B]))

repo: octsympy


text:
% Fast community finding algorithm by M. Newman
% Source: "Fast algorithm for detecting community
%                   structure in networks", Mark Newman
%
% INPUTs: adjacency matrix, nxn
% OUTPUTs: sequential group (cluster) formation,
%          modularity metric for each cluster breakdown
%
% Other functions used: numEdges.m
% Note: To save the plot generated in this routine:
%       uncomment "print newmanCommFast_example.pdf"
%
% GB: last updated, Oct 12 2012

function [groups_hist,Q]=newmanCommFast(adj)

groups={};
groups_hist={};  % a list of groups at each step
Q=[];
n=size(adj,1);   % number of vertices
groups_hist{1}={};

for i=1:n
  groups{i}=i;   % each node starts out as one community
  groups_hist{1}{i}=i;
end

Q(1) = Qfn(groups,adj);

for s=1:n-1 % all possible iterations

  dQ=zeros(length(groups));

  for i=1:length(groups)
    for j=i+1:length(groups)
      group_i=groups{i};
      group_j=groups{j};

      dQ(i,j)=0; % default

      % check if connected
      connected=0;
      if sum(sum(adj([group_i group_j],[group_i group_j])))>sum(sum(adj(group_i,group_i)))+sum(sum(adj(group_j,group_j)))
	connected=1;
      end

      if connected && not(isempty(group_i)) && not(isempty(group_j))
	% join groups i and j?
	dQ(i,j)=deltaQ(groups,adj,i,j);
      end

    end
  end

  % pick max nonzero dQ
  dQ(find(dQ==0))=-Inf;
  [minv,minij]=max(dQ(:));

  [mini,minj]=ind2sub([size(dQ,1),size(dQ,1)],minij);   % i and j corresponding to min dQ
  groups{mini}=sort([groups{mini} groups{minj}]);
  groups{minj}=groups{length(groups)}; % swap minj with last group
  groups=groups(1:length(groups)-1);

  groups_hist{length(groups_hist)+1}=groups; % save current snapshot
  Q(length(Q)+1) = Qfn(groups,adj);

end % end of all iterations

num_modules=[];
for g=1:length(groups_hist)
  num_modules=[num_modules length(groups_hist{g})];
end

set(gcf,'Color',[1 1 1],'Colormap',jet);
plot(num_modules,Q,'ko-')
xlabel('number of modules / communities')
ylabel('modularity metric, Q')

% print newmanCommFast_example.pdf

function dQ=deltaQ(groups,adj,i,j)

% computing the delta Q between two community configurations
% dQ = 2(e_ij - ai*aj) or (Q_groups_(i,j)_merged - Q_original)
% inputs: current community assignments: groups, adjacency matrix, i,j to be joined
% outputs: delta Q value (real number)

% $$$ Q1=Qfn(groups,adj);
% $$$ groups{i}=[groups{i} groups{j}];
% $$$ groups{j}=groups{length(groups)};
% $$$ groups=groups(1:length(groups)-1);
% $$$ Q2=Qfn(groups,adj);
% $$$ dQ = Q2-Q1;

% alternative dQ = 2(e_ij - ai*aj) from paper;
nedges=numEdges(adj); % compute the total number of edges
e_ii = numEdges( adj(groups{i},groups{i}) ) / nedges;
e_jj = numEdges( adj(groups{j},groups{j}) ) / nedges;
e_ij = numEdges( adj([groups{i} groups{j}],[groups{i} groups{j}]) ) / nedges - e_ii - e_jj;

a_i=sum(sum(adj(groups{i},:)))/nedges-e_ii;
a_j=sum(sum(adj(groups{j},:)))/nedges-e_jj;

dQ = 2*(e_ij-a_i*a_j);

function Q=Qfn(modules,adj) % ....... same as modularityMetric.m ........

% Computing the modularity for the final module break-down
% Defined as: Q=sum_over_modules_i (eii-ai^2) (eq 5) in Newman and Girvan.
% eij = fraction of edges that connect community i to community j
% ai=sum_j (eij)

nedges=numEdges(adj); % compute the total number of edges

Q = 0;
for m=1:length(modules)

  e_mm=numEdges(adj(modules{m},modules{m}))/nedges;
  a_m=sum(sum(adj(modules{m},:)))/nedges - e_mm;  % counting e_mm only once

  Q = Q + (e_mm - a_m^2);
end

author: aeolianine

repo: octave-networks-toolbox


text:
function h=Iabais(x)
L=length(x);
k=1;
i=1;
while i<=L/2
   h(i)=x(k)*256+x(k+1);
     k=k+2;
   i=i+1;
end;


text:
function clean = escapechars(dirty)
    % escape \ and % chars for sprintf
    clean = strrep(dirty,'\', '\\');
    clean = strrep(clean,'%', '%%');
end


text:
function idxCst = addCst(obj, cst2add)
% ADDCST Add constant(s) to problem

noCst2add = numel(cst2add); % no. of constants to add

% Add constants with 'direct indexing' (faster than other methods)
% Indices
if isempty(obj.cst)
    idxCst2add = 1:noCst2add;
else
    noCst = numel(obj.cst); % actual no. of constants
    idxCst2add = noCst+1:noCst+noCst2add;
end

% Add
obj.cst(idxCst2add,1) = cst2add;

% Indices of constants
idxCst = uint32(idxCst2add)';

end

author: pglira


text:
function [re_image] = imresize_gpu(image,re_param)
%This file implements the imresize function running on GPU based on
%vl_nnbilinearsampler.
% Input:
% image         -      image to be processed, gpuArray
% re_param      -   if its size is one, it is the ratio that to be resized
%               -   if it contains two numbers, it is the resized size.
% Output:
% re_image  -      the resized image.
% By Xin Li, 2018-6-5, UCM

if numel(re_param)==1
    sz = size(image);
    output_sz = round(sz(1:2)*re_param);
else
    output_sz = re_param;
end

yi = linspace(-1, 1, output_sz(1));
xi = linspace(-1, 1, output_sz(2));
[xx,yy] = meshgrid(xi,yi);
yyxx = single([yy(:), xx(:)]') ;
g = reshape(yyxx, 2, output_sz(1), output_sz(2), []);
re_image = vl_nnbilinearsampler(image, gpuArray(g));

end

author: XinLi-zn

repo: TADT

save_path: github-repos/MATLAB/XinLi-zn-TADT


text:
function b_n = beamWeightsDifferential2Spherical(a_n)
%BEAMWEIGHTSDIFFERENTIAL2SPHERICAL Transform differential to spherical coefficients.
%
%   For a set of a_n coefficients describing a pattern of a differential
%   array of the form
%   D(theta)=Sum_{n=1}^{N} a_n * cos^n(theta),
%   generate the beamweights for the same pattern in the SHD. Since the
%   pattern is axisymmetric only the N+1 coefficients of m=0 are returned.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% BEAMWEIGHTSDIFFERENTIAL2SPHERICAL.M - 10/4/2013
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% order N (or polynomial degree)
N = length(a_n)-1;

% build Legendre polynomial matrix
P_N = zeros(N+1);
for n = 0:N
    P_N(1:n+1,n+1) = returnLegePolyCoeffs(n);
end
% build the normalisation matrix
w = sqrt((2*(0:N)+1)/(4*pi));
W = diag(w);
W_inv = diag(1./w);
% convert coefficients
b_n = W_inv * inv(P_N) * a_n;

author: polarch

repo: Spherical-Array-Processing


text:
filename='Bridge_triangle_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'SLERP'; incrementFactor = 1;
filterType = 'P1';

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

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

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [D,S] = dictlearn(varargin)

% References:
%     "Online Learning for Matrix Factorization and Sparse Coding"
%     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro
%     arXiv:0908.0050
%
%     "Online Dictionary Learning for Sparse Coding"
%     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro
%     ICML 2009.

args = arg_define(varargin, ...
    arg_nogui({'X','Data'},[],[]), ...
    arg({'whiten','Whiten'}, true, [], 'Whiten the data.'), ...
    arg({'rescale','Rescale'}, true, [], 'Rescale the data.'), ...
    arg({'shuffle','Shuffle'}, true, [], 'Shuffle the data.'), ...
    arg({'haveMemory','HaveMemory'}, true, [], 'Have lots of memory.'), ...
    arg({'D','InitialGuess'},[],[],'Initial guess.'), ...
    arg({'K','NumComponents'},256,[], 'Number of components.'), ...
    arg({'iter','MaxIterations'},500,[], 'Maximum number of iterations.'), ...
    arg({'mode','ObjectiveFunction'},'(1/2)||x_i-Dalpha_i||_2^2 + lambda||alpha_i||_1 + lambda_2||alpha_i||_2^2',{'(1/2)||x_i-Dalpha_i||_2^2  s.t. ||alpha_i||_1 <= lambda','||alpha_i||_1  s.t. ||x_i-Dalpha_i||_2^2 <= lambda','(1/2)||x_i-Dalpha_i||_2^2 + lambda||alpha_i||_1 + lambda_2||alpha_i||_2^2','(1/2)||x_i-Dalpha_i||_2^2  s.t. ||alpha_i||_0 <= lambda','||alpha_i||_0  s.t.  ||x_i-Dalpha_i||_2^2 <= lambda','0.5||x_i-Dalpha_i||_2^2 +lambda||alpha_i||_0'},'Objective function to impose.'), ...
    arg({'modeD','DictionaryConstraints'},'||d_j||_2^2 <= 1',{'||d_j||_2^2 <= 1','||d_j||_2^2 + gamma1||d_j||_1 <= 1','||d_j||_2^2 + gamma1||d_j||_1 + gamma2 FL(d_j) <= 1','(1-gamma1)||d_j||_2^2 + gamma1||d_j||_1 <= 1'},'Constraints to impose on dictionary.'), ...
    arg({'modeParam','OptimizationMode'},'parameter-free (ICML)',{'parameter-free (ICML)','using rho (arXiv)','exponential decay weights (rho as memory factor, arXiv)'},'Optimization mode to use.'), ...
    arg({'lambda','Lambda'},1,[], 'Main penalty parameter.'), ...
    arg({'lambda2','Lambda2'},1,[], 'Second penalty parameter.'), ...
    arg({'posAlpha','PositiveActivations'},false,[], 'Enforce positive activations. Not compatible with learning modes 3 and 4.'), ...
    arg({'posD','PositiveDictionary'},false,[], 'Enforce positive dictionary. Not compatible with constraint set 2.'), ...
    arg({'gamma1','GammaParam1'},[],[], 'Gamma parameter 1 for dictionary constraint set.'), ...
    arg({'gamma2','GammaParam2'},[],[], 'Gamma parameter 2 for dictionary constraint set.'), ...
    arg({'batchsize','MinibatchSize'},512,[], 'Size of the mini-batches to use.'), ...
    arg({'iter_updateD','NumBCDIterations'},1,[], 'Number of block-coordinate descent iterations. For the dictionary update step.'), ...
    arg({'rho','TuningParameterRho'},[],[], 'Rho tuning parameter. For online opt mode.'), ...
    arg({'t0','TuningParameterT0'},[],[], 'T0 tuning parameter. For online opt mode.'), ...
    arg({'clean','PruneUnused'},true,[], 'Prune unused elements.'), ...
    arg({'verbose','Verbose'},true,[], 'Verbose output.'), ...
    arg({'numThreads','NumThreads'},16,[], 'Number of threads.'), ...
    arg_nogui({'chanlocs','Chanlocs'},16,[], 'Channel locations.'));

% fix up arguments
if isempty(args.D)
    args = rmfield(args,'D'); end
if isempty(args.gamma1)
    args = rmfield(args,'gamma1'); end
if isempty(args.gamma2)
    args = rmfield(args,'gamma2'); end
if isempty(args.rho)
    args = rmfield(args,'rho'); end
args.modeParam = hlp_rewrite(args.modeParam,'parameter-free (ICML)',0, 'using rho (arXiv)',1, ...
    'exponential decay weights (rho as memory factor, arXiv)',3);
args.mode = hlp_rewrite(args.mode, ...
    '(1/2)||x_i-Dalpha_i||_2^2  s.t. ||alpha_i||_1 <= lambda',0, ...
    '||alpha_i||_1  s.t. ||x_i-Dalpha_i||_2^2 <= lambda',1, ...
    '(1/2)||x_i-Dalpha_i||_2^2 + lambda||alpha_i||_1 + lambda_2||alpha_i||_2^2',2, ...
    '(1/2)||x_i-Dalpha_i||_2^2  s.t. ||alpha_i||_0 <= lambda',3,...
    '||alpha_i||_0  s.t.  ||x_i-Dalpha_i||_2^2 <= lambda',4,...
    '0.5||x_i-Dalpha_i||_2^2 +lambda||alpha_i||_0',5);
args.modeD = hlp_rewrite(args.modeD, ...
    '||d_j||_2^2 <= 1',0, ...
    '||d_j||_2^2 + gamma1||d_j||_1 <= 1',1,...
    '||d_j||_2^2 + gamma1||d_j||_1 + gamma2 FL(d_j) <= 1',2, ...
    '(1-gamma1)||d_j||_2^2 + gamma1||d_j||_1 <= 1',3);

X = args.X;

%args.whiten = false;

% Whiten and rescale the data
if args.whiten
    sphere = inv(real(sqrtm(cov_robust(X'))));
else
    sphere = eye(size(X,1));
end
X = sphere*X;

if args.rescale
    m = sqrt(sum(X.^2) + (1e-8));
    X = bsxfunwrap(@rdivide,X,m);
end

if args.shuffle
    idx = shuffle(1:size(X,2));
    X = X(:,idx);
end

%% Run the dictionary learning
if args.haveMemory
    D = mexTrainDL_Memory(X,args);
else
    D = mexTrainDL(X,args);
end
S = sphere;
%Df = inv(S)*D;

repo: REST


text:
function err = save_nifti(hdr,niftifile)
% err = save_nifti(hdr,niftifile)
%
% Pixel data should be in hdr.vol
%
% Handles data structures with more than 32k cols by setting
% hdr.dim(2) = -1 and hdr.glmin = ncols. This is FreeSurfer specific,
% for handling surfaces. The exception to this is when the total
% number of spatial voxels equals 163842, then the volume is
% reshaped to 27307x1x6xnframes. This is for handling the 7th
% order icosahedron used by FS group analysis.
%

%
% save_nifti.m
%
% Original Author: Doug Greve
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

err = 1;
if(nargin ~= 2)
  fprintf('err = save_nifti(hdr,niftifile)\n');
  return;
end

ext = niftifile((strlen(niftifile)-2):strlen(niftifile));
if(strcmpi(ext,'.gz'))
  gzip_needed = 1;
  niftifile = niftifile(1:strlen(niftifile)-3);
  %fprintf('First, saving to %s before compressing\n',niftifile);
else
  gzip_needed = 0;
end

% Check for ico7
sz = size(hdr.vol);
if(sz(1) == 163842)
  %fprintf('save_nifti: ico7 reshaping\n');
  dim = [27307 1 6 size(hdr.vol,4)];
  hdr.vol = reshape(hdr.vol, dim);
end

fp = fopen(niftifile,'w');
if(fp == -1)
  fprintf('ERROR: could not open %s\n',niftifile);
  return;
end

hdr.data_type = [hdr.data_type(:)' repmat(' ',[1 10])];
hdr.data_type = hdr.data_type(1:10);

hdr.db_name = [hdr.db_name(:)' repmat(' ',[1 18])];
hdr.db_name = hdr.db_name(1:18);

hdr.dim = ones(1,8);
if size(hdr.vol,4)>1
  hdr.dim(1) = 4;
else
  hdr.dim(1) = 3;
end
hdr.dim(2) = size(hdr.vol,1);
hdr.dim(3) = size(hdr.vol,2);
hdr.dim(4) = size(hdr.vol,3);
hdr.dim(5) = size(hdr.vol,4);

% This is to accomodate structures with more than 32k cols
% FreeSurfer specific. See also mriio.c.
if(hdr.dim(2) > 2^15)
  hdr.glmin = hdr.dim(2);
  hdr.dim(2) = -1;
end

hdr.pixdim = [hdr.pixdim(:)' repmat(0,[1 8])];
hdr.pixdim = hdr.pixdim(1:8);

hdr.descrip = [hdr.descrip(:)' repmat(' ',[1 80])];
hdr.descrip = hdr.descrip(1:80);

hdr.aux_file = [hdr.aux_file(:)' repmat(' ',[1 24])];
hdr.aux_file = hdr.aux_file(1:24);

hdr.intent_name = [hdr.intent_name(:)' repmat(' ',[1 16])];
hdr.intent_name = hdr.intent_name(1:16);

hdr.magic = [hdr.magic(:)' repmat(0,[1 4])];
hdr.magic = hdr.magic(1:4);

if(isempty(hdr.regular))  hdr.regular  = ' '; end
if(isempty(hdr.dim_info)) hdr.dim_info = ' '; end
if(isempty(hdr.slice_code)) hdr.slice_code = ' '; end
if(isempty(hdr.xyzt_units)) hdr.xyzt_units = ' '; end % should be err

hdr.vox_offset = 352; % not 348
fwrite(fp,348,'int');
fwrite(fp,hdr.data_type,    'char');
fwrite(fp,hdr.db_name,      'char');
fwrite(fp,hdr.extents,      'int');
fwrite(fp,hdr.session_error,'short');
fwrite(fp,hdr.regular,      'char');
fwrite(fp,hdr.dim_info,     'char');
fwrite(fp,hdr.dim,          'short');
fwrite(fp,hdr.intent_p1,    'float');
fwrite(fp,hdr.intent_p2,    'float');
fwrite(fp,hdr.intent_p3,    'float');
fwrite(fp,hdr.intent_code,  'short');
fwrite(fp,hdr.datatype,     'short');
fwrite(fp,hdr.bitpix,       'short');
fwrite(fp,hdr.slice_start,  'short');
fwrite(fp,hdr.pixdim,       'float');
fwrite(fp,hdr.vox_offset,   'float');
fwrite(fp,hdr.scl_slope,    'float');
fwrite(fp,hdr.scl_inter,    'float');
fwrite(fp,hdr.slice_end,    'short');
fwrite(fp,hdr.slice_code,   'char');
fwrite(fp,hdr.xyzt_units,   'char');
fwrite(fp,hdr.cal_max,      'float');
fwrite(fp,hdr.cal_min,      'float');
fwrite(fp,hdr.slice_duration,'float');
fwrite(fp,hdr.toffset,       'float');
fwrite(fp,hdr.glmax,         'int');
fwrite(fp,hdr.glmin,         'int');
fwrite(fp,hdr.descrip,       'char');
fwrite(fp,hdr.aux_file,      'char');
fwrite(fp,hdr.qform_code,    'short');
fwrite(fp,hdr.sform_code,    'short');
fwrite(fp,hdr.quatern_b,     'float');
fwrite(fp,hdr.quatern_c,     'float');
fwrite(fp,hdr.quatern_d,     'float');
fwrite(fp,hdr.quatern_x,     'float');
fwrite(fp,hdr.quatern_y,     'float');
fwrite(fp,hdr.quatern_z,     'float');
fwrite(fp,hdr.srow_x,        'float');
fwrite(fp,hdr.srow_y,        'float');
fwrite(fp,hdr.srow_z,        'float');
fwrite(fp,hdr.intent_name,   'char');
fwrite(fp,hdr.magic,         'char');

% Pad to get to 352 bytes (header size is 348)
fwrite(fp,0,'char');
fwrite(fp,0,'char');
fwrite(fp,0,'char');
fwrite(fp,0,'char');

npix = prod(size(hdr.vol));
switch(hdr.datatype)
 case   2, nitemswritten = fwrite(fp,hdr.vol,'uchar'); % dont use char
 case   4, nitemswritten = fwrite(fp,hdr.vol,'short');
 case   8, nitemswritten = fwrite(fp,hdr.vol,'int');
 case  16, nitemswritten = fwrite(fp,hdr.vol,'float');
 case  64, nitemswritten = fwrite(fp,hdr.vol,'double');
 case 256, nitemswritten = fwrite(fp,hdr.vol,'int8');
 case 512, nitemswritten = fwrite(fp,hdr.vol,'ushort');
 case 768, nitemswritten = fwrite(fp,hdr.vol,'uint');
 otherwise,
  fprintf('WARNING: data type %d not supported, but writing as float',...
	  hdr.datatype);
  nitemswritten = fwrite(fp,hdr.vol,'float');
  return;
end
fclose(fp);

if(npix ~= nitemswritten)
  fprintf('ERROR: tried to write %d, but only wrote %d',npix,nitemswritten);
  return;
end

if(gzip_needed)
  cmd = sprintf('gzip -f %s', niftifile);
  %fprintf('Compressing with\n');
  %fprintf('   %s\n',cmd);
  unix(cmd);
end

err = 0;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ xyz, line_pointer, line_data ] = xyzl_example ( point_num, ...
  line_num, line_data_num )

%*****************************************************************************80
%
%% XYZL_EXAMPLE sets data suitable for a pair of XYZ and XYZL files.
%
%  Discussion:
%
%    There are 8 points.
%    There are 6 lines.
%    There are 18 line items.
%
%       8------7
%      /|     /|
%     / |    / |
%    5------6  |
%    |  4---|--3
%    | /    | /
%    |/     |/
%    1------2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, integer LINE_NUM, the number of lines.
%
%    Input, integer  LINE_DATA_NUM, the number of line items.
%
%    Output, real XY(2,POINT_NUM), the point coordinates.
%
%    Output, integer LINE_POINTER(LINE_NUM+1), pointers to the
%    first line item for each line.
%
%    Output, integer LINE_DATA(LINE_DATA_NUM), indices
%    of points that form lines.
%
  xyz(1:3,1:point_num) = [ ...
     0.0,  0.0,  0.0; ...
     1.0,  0.0,  0.0; ...
     1.0,  1.0,  0.0; ...
     0.0,  1.0,  0.0; ...
     0.0,  0.0,  1.0; ...
     1.0,  0.0,  1.0; ...
     1.0,  1.0,  1.0; ...
     0.0,  1.0,  1.0 ]';

  line_pointer(1:line_num+1) = [ 1, 6, 11, 13, 15, 17, 19 ];

  line_data(1:line_data_num) = [ ...
     1,  2,  3,  4,  1, ...
     5,  6,  7,  8,  5, ...
     1,  5, ...
     2,  6, ...
     3,  7, ...
     4,  8 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Run setupThermoModel for the E. coli metabolic reconstruction iAF1260

% Requires an installation of:
% OpenBabel

%% Add required fields and directories to path
initVonBertalanffy

%todo, clean up this so it is not needed
global CBTDIR
pth=which('initCobraToolbox.m');
CBTDIR = pth(1:end-(length('initCobraToolbox.m')+1));

cd([CBTDIR filesep 'test' filesep 'testVonBertalanffy'])

%% Configure inputs
load iAF1260
if model.S(952,350)==0
    model.S(952,350)=1; % One reaction needing mass balancing in iAF1260
end
model.metCharges(strcmp('asntrna[c]',model.mets))=0; % One reaction needing charge balancing

molfileDir = 'iAF1260Molfiles'; % Directory containing molfiles

cid = []; % KEGG Compound identifiers. Not required since molfile directory is specified.

T = 310.15; % Temperature in Kelvin
cellCompartments = ['c'; 'e'; 'p']; % Cell compartment identifiers
ph = [7.7; 7.7; 7.7]; % Compartment specific pH
is = [0.25; 0.25; 0.25]; % Compartment specific ionic strength in mol/L
chi = [0; 90; 90]; % Compartment specific electrical potential relative to cytosol in mV

xmin = 1e-5*ones(size(model.mets)); % Lower bounds on metabolite concentrations in mol/L
xmax = 0.02*ones(size(model.mets)); % Upper bounds on metabolite concentrations in mol/L

confidenceLevel = 0.95; % Confidence level for estimated standard transformed reaction Gibbs energies. Used to quantitatively assign reaction directionality.

%% Call setupThermoModel
modelT = setupThermoModel(model,molfileDir,cid,T,cellCompartments,ph,is,chi,xmin,xmax,confidenceLevel);

save('iAF1260Thermo_test.mat', 'modelT', '-v7');

%% Compare test results to expected results
clear all;

old = load('iAF1260Thermo.mat');
new = load('iAF1260Thermo_test.mat');

% Check for differences in estimated standard transformed Gibbs energies of formation
fig = figure(1);
subplot(1,3,1);
rmse1 = sqrt(mean( (new.modelT.DfGt0 - old.modelT.DfGt0).^2 ));
fprintf('RMSE difference between the old and new DfGt0: %g\n', rmse1);
cdfplot(abs((new.modelT.DfGt0 - old.modelT.DfGt0)));
xlabel('|D_f G^{\prime\circ}(new) - D_f G^{\prime\circ}(old)|');
title(['\Delta_f G^{\prime\circ} RMSE = ' sprintf('%g', rmse1)]);

% Check for differences in estimated standard transformed reaction Gibbs energies
subplot(1,3,2);
rmse2 = sqrt(mean( (new.modelT.DrGt0 - old.modelT.DrGt0).^2 ));
fprintf('RMSE difference between the old and new DrGt0: %g\n', rmse2);
cdfplot(abs((new.modelT.DrGt0 - old.modelT.DrGt0)));
xlabel('|D_r G^{\prime\circ}(new) - D_r G^{\prime\circ}(old)|');
title(['\Delta_r G^{\prime\circ} RMSE = ' sprintf('%g', rmse2)]);

% Check for differences in uncertainty levels - indicative of differences in coverage
subplot(1,3,3);
rmse3 = sqrt(mean( (new.modelT.uf - old.modelT.uf).^2 ));
fprintf('RMSE difference between the old and new uf: %g\n', rmse3);
cdfplot(abs((new.modelT.uf - old.modelT.uf)));
xlabel('|U_f (new) - U_f (old)|');
title(['U_f RMSE = ' sprintf('%g', rmse3)]);

print(fig, 'iAF1260_compare.eps', '-deps');

author: opencobra

repo: cobratoolbox


text:
%function f = s6990029.m
%This function was completed on Tuesday, 16 November 2002
%as a part of a course, in order to form a good-interfaced program
%for producing random variables from the negative binomial distribution.
%
%Made by P.B.(6990029)
function s6990029

fprintf('In a few seconds a pop-up screen will prompt you to choose from one\nof 4 methods, in order to generate random variables from the Negative Binomial Distribution.');
fprintf('\n\nLOADING('),
for j=1:30
fprintf('---')
pause(0.06990029)
end
fprintf(')')
fprintf('\n')

K=menu('Choose a method for generating n random variables from the Neg.Binomial :','1) Inverse transform Method','2) Convolution Method (using Bernoulli D.)','3) Convolution Method (using Geometrical D.)','4) Acceptance Rejection Method','**************************Which is faster=better***********************','Exit');

%Method 1

if (K==1)

disp('You have selected the Inverse Transform method. Please wait...');
pause(2);
clc;clear;
N = input('Please give me the size of the sample : ');
R = input('Type the value of the parameter of successes R : ');
P = input('Now give the probability of success P : ');
disp('Please standby...');

%set up the storage space and the clock
tic
X = zeros(1,N);
x =1:N;

for x =1:N
pr(x)  = NBINCDF(x,R,P);
end

pr(N)=1;    %epeidh en vgainei panta 1 to teleutaio stoixeio tou pinaka ths a8roistikhs ...kammia fora... dhmiourgei provlhmata!!!

for i=1:N
 s=1;
u= rand;
while u > pr(s)
    s=s+1;
end
X(i) = s;
end

%stop the clock
t=toc;
fp=fopen('c:\data1.txt','w');
fprintf(fp,' %3.0f \n',X);
fclose(fp);
disp('Your sample has been saved in data1.txt');
fprintf('The time needed, depending on your machine'' speed \n and the complexity of the algorithm, was: ');disp(t)
pause(2);

load handel;
sound(y,Fs);

%the empirical histogram

[n,h]= hist(X,10);
n=n/(h(2)-h(1))/N;
bar(h,n,1,'w')
hold on

%theoritical curve

for x =1:N
pr(x)  = NBINPDF(x,R,P);
end
x =1:N;
plot(x,pr,'k');
hold off
title('Empirical & Theoritical pdf of the Neg. Binomial Distribution');
xlabel('X');
ylabel('P(X)');

%Method 2

elseif (K==2)

fprintf('You have selected the Convolution method \n(using as a prime several Bernoulli variables). \n\nPlease wait...');
pause(4);
clc;clear;
N = input('Please give me the size of the sample : ');
R = input('Type the value of the parameter of successes R : ');
P = input('Now give the probability of the bernoulli trials P : ');
disp('Please standby...');

%set up the storage space and the clock
tic

X = zeros(1,N);
for j=1:N
h=0;
i=0;
while h < R
u=rand;
i=i+1;
if  u <= P
h = h+1;
end
end
X(j)=i-R;               % duskolo shmeio...orismos:...o arithmos twn dokimwn mexri R epituxies!
end

t=toc;
fp=fopen('c:\data2.txt','w');
fprintf(fp,' %3.0f \n',X);
fclose(fp);
disp('Your sample has been saved in data2.txt');
fprintf('The time needed, depending on your machine'' speed \n and the complexity of the algorithm, was: ');disp(t)
pause(2);

load handel;
sound(y,Fs);

%the empirical histogram

[n,h]= hist(X,10);
n=n/(h(2)-h(1))/N;
bar(h,n,1,'w')
hold on

%theoritical curve

for x =1:N
pr(x)  = NBINPDF(x,R,P);
end
x =1:N;
plot(x,pr,'k');
hold off
title('Empirical & Theoritical pdf of the Neg. Binomial Distribution');
xlabel('X');
ylabel('P(X)');

%METHOD 3

elseif (K==3)

fprintf('You have selected the Convolution method \n(using as a prime several Geometrical variables). \n\nPlease wait...');
pause(4);
clc;clear;
N = input('Please give me the size of the sample : ');
R = input('Type the value of the parameter of successes R : ');
P = input('Now give the probability of the bernoulli trials P : ');
disp('Please standby...');

%set up the storage space and the clock
tic

X = zeros(1,N);

G=geornd(P,N,R);

X=G*ones(R,1);

%stop the clock
t=toc;
fp=fopen('c:\data3.txt','w');
fprintf(fp,' %3.0f \n',X);
fclose(fp);
disp('Your sample has been saved in data3.txt');
fprintf('The time needed, depending on your machine'' speed \n and the complexity of the algorithm, was: ');disp(t)
pause(2);

load handel;
sound(y,Fs);

%the empirical histogram

[n,h]= hist(X,10);
n=n/(h(2)-h(1))/N;
bar(h,n,1,'w')
hold on

%theoritical curve

for x =1:N
pr(x)  = NBINPDF(x,R,P);
end
x =1:N;
plot(x,pr,'k');
hold off
title('Empirical & Theoritical pdf of the Neg. Binomial Distribution');
xlabel('X');
ylabel('P(X)');

%METHOD 4

elseif (K==4)

fprintf('You have selected the Acceptance-Rejection method). \n\nPlease wait...');
pause(4);
clc;clear;
N = input('Please give me the size of the sample : ');
R = input('Type the value of the parameter of successes R : ');
P = input('Now give the probability of the bernoulli trials P : ');
disp('Please standby...');

%set up the storage space and the clock
tic

for x =1:N
pr(x)  = NBINPDF(x,R,P);
end

c=max(pr)*N;

irv=1;
while irv<=N
    y = unidrnd(N);
    u = rand(1);
    if u <= (pr(y))/c ;
        X(irv) =  y;
        irv = irv +1;
    end
end

%stop the clock
t=toc;
%save to disk
fp=fopen('c:\data4.txt','w');
fprintf(fp,' %3.0f \n',X);
fclose(fp);
disp('Your sample has been saved in data4.txt');
fprintf('The time needed, depending on your machine'' speed \n and the complexity of the algorithm, was: ');disp(t)
pause(2);

load handel;
sound(y,Fs);

%the empirical histogram

[n,h]= hist(X,10);
n=n/(h(2)-h(1))/N;
bar(h,n,1,'w')
hold on

%theoritical curve

for x =1:N
pr(x)  = NBINPDF(x,R,P);
end
x =1:N;
plot(x,pr,'k');
hold off
title('Empirical & Theoritical pdf of the Neg. Binomial Distribution');
xlabel('X');
ylabel('P(X)');

%Test for overall performance within the range of the two parameters (R,P)

%#########################################################################################

elseif (K==5)

    clc;
    fprintf('You have selected to test between the three first  methods.\nFor the 4th one the envelope function selected was rather unsuitable since it was\nthe discrete uniform.\nWith an other function as an envelope it would be very hard to test\nthrough a wide range of (R,P)\nsince the Neg. Binomial''s shape would change.\n\nPlease wait. It will take some time ...\n');
pause(5);

clear;
N = input('Please give me the size of the sample. It would be most preferable to choose sample bigger than 10000! : ');

fprintf('Please go drink a coffee...\n')

    q=1;

    for R=10:N
    for P=0.3:0.2:0.7
        %1st
tic
X = zeros(1,N);
x =1:N;

for x =1:N
pr(x)  = NBINCDF(x,R,P);
end

pr(N)=1;    %epeidh en vgainei panta 1 to teleutaio stoixeio tou pinaka ths a8roistikhs ...kammia fora... dhmiourgei provlhmata!!!

for i=1:N
 s=1;
u= rand;
while u> pr(s)
    s=s+1;
end
    X(i) = s;
end
%stop the clock
t1=toc;

%2nd

tic

X = zeros(1,N);
for j=1:N
h=0;
i=0;
while h < R
u=rand;
i=i+1;
if  u <= P
h = h+1;
end
end
X(j)=i-R;               % duskolo shmeio...orismos:...o arithmos twn dokimwn mexri R epituxies!
end

t2=toc;

%3rd

tic

X = zeros(1,N);

G=geornd(P,N,R);

X=G*ones(R,1);

%stop the clock
t3=toc;

q=q+1;

T(q,:)=[t1 t2 t3];

end
end

boxplot(T);
anova1(T);

fprintf('\n\n\n\n\nOne would choose the method which uses the sum of the geometrical variables\nsince the Anova table also recommends this method!\nHave a nice day.')

else
 disp('Have a nice Day!!!');

%break

end


text:
classdef SPX_ShortTermFourierTransform < handle

properties(SetAccess=private)
    % Length of block
    BlockLength
    % Number of samples in overlaps
    OverlappedSamplesCount
    % FFT length
    FFTLength
    % Window type
    WindowType
    % Window
    Window
end

methods
    function self  = SPX_ShortTermFourierTransform(...
        blockLength, overlappedSamples, fftLength, windowType)
        self.BlockLength  = blockLength;
        self.OverlappedSamplesCount = overlappedSamples;
        self.FFTLength = fftLength;
        self.WindowType = windowType;
        if strcmp(windowType, 'hamming')
            self.Window = hamming(blockLength)';
        elseif strcmp(windowType,'blackman')
            self.Window = blackman(blockLength)';
        else
            error('Unsupported window type');
        end
    end

    function [ spec, timePoints]  = run(self, signal)
        blockLength = self.BlockLength;
        overlappedSamples = self.OverlappedSamplesCount;
        N = length(signal);
        NBlocks = spx.dsp.overlapped_blocks_count(N, blockLength...
            , overlappedSamples);
        w = self.Window;
        shift = blockLength - overlappedSamples;
        fftLength = self.FFTLength;
        start = 0;
        for ii=1:NBlocks
            % Let's identify the start and end of the block
            range = start + (1:blockLength);
            % Let's read off the block and apply windowing
            x = signal(range) .* w;
            % Identify normalized start time for this block
            timePoints(ii) = start +1;
            % Compute fft of this
            Xf  = fft(x, fftLength);
            spec(:, ii) = Xf;
            start = start + shift;
        end
    end
end

end

author: indigits

repo: sparse-plex