text:
function sim = gaussianKernel(x1, x2, sigma)
%RBFKERNEL returns a radial basis function kernel between x1 and x2
%   sim = gaussianKernel(x1, x2) returns a gaussian kernel between x1 and x2
%   and returns the value in sim

% Ensure that x1 and x2 are column vectors
x1 = x1(:); x2 = x2(:);

% You need to return the following variables correctly.
sim = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the similarity between x1
%               and x2 computed using a Gaussian kernel with bandwidth
%               sigma
%
%

upper = sum((x1 .- x2) .^ 2);
sim = exp(-upper / (2 * sigma ^ 2));

% =============================================================

end

author: ecmadao

repo: Coding-Guide


text:
function [rtk,sat_,stat]=estpos(rtk,obs,nav,sv,opt)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%estimate reciever position and clock bias
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global  glc
NX=3+glc.NSYS; MAXITER=10; iter=1;
time=obs(1).time; xr0=[rtk.sol.pos';zeros(glc.NSYS,1)];

while iter<=MAXITER

    % compute residual,measurement model,weight matrix
    [v,H,P,vsat,azel,resp,nv,ns]=rescode(iter,obs,nav,sv,xr0,opt);
    sat_.vsat=vsat;
    sat_.azel=azel;
    sat_.resp=resp;

    if nv<NX, stat=0;return;end

    % least square algrithm
    [dx,Q]=least_square(v,H,P);

    xr0=xr0+dx;
    iter=iter+1;

    if dot(dx,dx)<1e-4
        rtk.sol.time=timeadd(obs(1).time,-xr0(4)/glc.CLIGHT);
        rtk.sol.ns  =ns;
        rtk.sol.ratio=0;
        rtk.sol.pos =xr0(1:3)';
        rtk.sol.vel =zeros(1,3);
        rtk.sol.posP(1)=Q(1,1);rtk.sol.posP(2)=Q(2,2);rtk.sol.posP(3)=Q(3,3);
        rtk.sol.posP(4)=Q(1,2);rtk.sol.posP(5)=Q(2,3);rtk.sol.posP(6)=Q(1,3);
        rtk.sol.dtr(1) =xr0(4)/glc.CLIGHT;
        rtk.sol.dtr(2) =xr0(5)/glc.CLIGHT;
        rtk.sol.dtr(3) =xr0(6)/glc.CLIGHT;
        rtk.sol.dtr(4) =xr0(7)/glc.CLIGHT;
        rtk.sol.dtr(5) =xr0(8)/glc.CLIGHT;

        % validate solution
        stat=valsol(time,v,P,vsat,azel,opt);
        if stat==1
            rtk.sol.stat=glc.SOLQ_SPP; return;
        end
        return;
    end

end

if iter>MAXITER
    stat=0;
    [week,sow]=time2gpst(time);
    fprintf('Warning:GPS week = %d sow = %.3f,SPP iteration divergent!\n',week,sow);
    return;
end

return

repo: GINav


text:
%SerialLink.ikine6s Analytical inverse kinematics
%
% Q = R.ikine6s(T) is the joint coordinates (1xN) corresponding to the robot
% end-effector pose T represented by an SE(3) homogenenous transform (4x4).  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)
%
% 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 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.

function theta = ikine6s(robot, T, 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

    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

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

    %% 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)
                warning('RTB:ikine6s:notreachable', 'point not reachable');
                theta = [NaN NaN NaN NaN NaN NaN];
                return
            end
            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);

        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
    % 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 = transl(0, 0, L(4).d);      % Tz(d4)
    Tt = transl(L(6).a, 0, L(6).d) * trotx(L(6).alpha);  % Tz(d6) Tx(a5) Rx(alpha6)

    R = inv(Td4) * inv(T13) * 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 i=1:robot.n   %#ok<*AGROW>
        theta(i) = theta(i) - L(i).offset;
    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).sigma] == 0) && ...
            (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).sigma] == 0);
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).sigma] == [0 0 1]);
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).sigma] == 0);
end

author: Allopart

repo: rbpf-gmapping


text:
function output = callmosek(model)

% Retrieve needed data
options = model.options;
F_struc = model.F_struc;
c       = model.c;
Q       = model.Q;
K       = model.K;
x0      = model.x0;
integer_variables = model.integer_variables;
binary_variables = model.binary_variables;
extended_variables = model.extended_variables;
ub      = model.ub;
lb      = model.lb;
mt      = model.monomtable;

% *********************************
% What type of variables do we have
% *********************************
model.linear_variables = find((sum(abs(mt),2)==1) & (any(mt==1,2)));
model.nonlinear_variables = setdiff((1:size(mt,1))',model.linear_variables);
model.sigmonial_variables = find(any(0>mt,2) | any(mt-fix(mt),2));

% Some meta-solver thinks we handle binaries
if ~isempty(model.binary_variables)
    integer_variables = union(model.integer_variables, model.binary_variables);
    if isempty(lb)
        lb = repmat(-inf,1,length(model.c));
    end
    lb(model.binary_variables) = max(lb(model.binary_variables),0);
    if isempty(ub)
        ub = repmat(inf,1,length(model.c));
    end
    ub(model.binary_variables) = min(ub(model.binary_variables),1);
    model.lb = lb;
    model.ub = ub;
    model.integer_variables = integer_variables;
end

% Some meta solvers might construct model with empty cones
if any(model.K.s) && any(model.K.s == 0)
    model.K.s(model.K.s==0)=[];
end
if any(model.K.q) && any(model.K.q == 0)
    model.K.q(model.K.q==0)=[];
end

if ~isempty(model.sigmonial_variables) | isequal(model.solver.version,'GEOMETRIC')
    [x,D_struc,problem,r,res,solvertime,prob] = call_mosek_geometric(model);
else
    [x,D_struc,problem,r,res,solvertime,prob] = call_mosek_lpqpsocpsdp(model);
end

% YALMIP has introduced internal variables for socp/exp cones etc
if length(x) > 0 && length(x) ~= length(model.c)
    x = x(1:length(model.c));
end

infostr = yalmiperror(problem,'MOSEK');

% Save all data sent to solver?
if options.savesolverinput
    solverinput.prob = prob;
    solverinput.param = options.mosek;
else
    solverinput = [];
end

% Save all data from the solver?
if options.savesolveroutput
    solveroutput.r = r;
    solveroutput.res = res;

else
    solveroutput = [];
end

% Standard interface
output = createOutputStructure(x,D_struc,[],problem,infostr,solverinput,solveroutput,solvertime);

function [x,D_struc,problem,r,res,solvertime,prob] = call_mosek_lpqpsocpsdp(model);

[model,output] = normalizeExponentialCone(model);
if output.problem
    problem = output.problem;
    x = [];
    D_struc = [];
    r = [];
    res = [];
    solvertime = 0;
    prob = [];
else
    if nnz(model.Q)==0 && isempty(model.integer_variables) && isempty(model.x0)
        % Standard cone problem which we can model by sending our standard dual
        % and then recover solution via Moseks dual
        [x,D_struc,problem,r,res,solvertime,prob] = call_mosek_dual(model);
    else
        % Integer conic program
        % Quadratic objective
        % Exponential cones
        [x,D_struc,problem,r,res,solvertime,prob] = call_mosek_primal(model);
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function bsv_test04 ( )

%*****************************************************************************80
%
%% BSV_TEST04 varies the location of the left boundary A in the Burgers equation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BSV_TEST04:\n' );
  fprintf ( 1, '  Solution of steady viscous Burgers equation.\n' );
  fprintf ( 1, '  Vary the left boundary location A around the value -1.\n' );

  a_test = [ -1.04, -1.02, -1.01, -1.005, -1.0, -0.995, -0.99, -0.98, -0.96 ];
  test_num = length ( a_test );
  b = +1.0;
  alpha = 1.0;
  beta = -1.0;
  nu = 0.1;
  n = 161;
  output = 0;

  u = zeros(n,test_num);

  fprintf ( 1, '\n' );

  for test = 1 : test_num

    a = a_test(test);

    fprintf ( 1, '  Using A = %g\n', a );

    u(:,test) = bsv ( a, b, alpha, beta, nu, n, output );

  end

  x = linspace ( a, b, n );

  figure ( 4 )
  clf
  hold on
  plot ( x, u, 'b-', 'LineWidth', 3 );
  plot ( [a,b], [0.0,0.0], 'r-', 'LineWidth', 2 )
  grid on
  title ( 'A = -1.04, -1.02, -1.01, -1.005, -1, -0.995, -0.99, -0.98, -0.96' )
  xlabel ( '<--- X --->' )
  ylabel ( '<---U(X) --->' )
  axis ( [ a, b, -1.5, alpha ] )
  hold off
  filename = 'bsv_test04.png';
  print ( '-dpng', filename )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saved plot to file "%s".\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function I = domIntegral4(data)
%+========================================================================+
%|                                                                        |
%|              OPENDOM - LIBRARY FOR NUMERICAL INTEGRATION               |
%|           openDom is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal & Francois Alouges (c) 2017-2018.          |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%|             francois.alouges@polytechnique.edu                         |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : domIntegral4.m                                |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Numerical integation with 4 input arguments   |
%|  `---'  |                                                              |
%+========================================================================+

%%% FINITE ELEMENT OPERATOR --> \int_{mesh(x)} psi(x)' f(x) psi(x) dx
if isa(data{1},'dom') && isa(data{2},'fem')
    % Domain with quadrature
    Xdom   = data{1};
    [X,Wx] = Xdom.qud;
    Nx     = size(X,1);
    Wx     = spdiags(Wx,0,Nx,Nx);

    % Integrated finite element matrix
    u  = data{2};
    Mu = u.uqm(Xdom);
    if iscell(Mu)
        Mu{1} = Mu{1}' * Wx;
        Mu{2} = Mu{2}' * Wx;
        Mu{3} = Mu{3}' * Wx;
    else
        Mu = Mu' * Wx;
    end

    % Function evaluation on quadrature
    F = data{3};
    if iscell(F)
        Fx{1} = spdiags(F{1}(X),0,Nx,Nx);
        Fx{2} = spdiags(F{2}(X),0,Nx,Nx);
        Fx{3} = spdiags(F{3}(X),0,Nx,Nx);
    else
        Fx = spdiags(F(X),0,Nx,Nx);
    end

    % Finite element matrix
    v  = data{4};
    Mv = v.uqm(Xdom);

    % Integration
    I = femMultiplyCell(Mu,Fx,Mv);

%%% BOUNDARY ELEMENT INTEGRATION --> \int_{mesh(y)} f(x,y) psi(y) dy
elseif isnumeric(data{1}) && isa(data{2},'dom')
    % Evaluation points
    X  = data{1};
    Nx = size(X,1);

    % Domain with quadrature
    Ydom   = data{2};
    [Y,Wy] = Ydom.qud;
    Ny     = size(Y,1);
    Wy     = spdiags(Wy,0,Ny,Ny);

    % Function evaluation on quadrature
    F = data{3};
    if iscell(F)
        Fxy = {zeros(Nx,Ny),zeros(Nx,Ny),zeros(Nx,Ny)};
        for i = 1:3
            for j = 1:Ny
                Fxy{i}(:,j) = F{i}(X,Y(j,:));
            end
        end
    else
        Fxy = zeros(Nx,Ny);
        for j = 1:Ny
            Fxy(:,j) = F(X,Y(j,:));
        end
    end

    % Integrated finite element matrix
    v  = data{4};
    Mv = v.uqm(Ydom);
    if iscell(Mv)
        Mv{1} = Wy * Mv{1};
        Mv{2} = Wy * Mv{2};
        Mv{3} = Wy * Mv{3};
    else
        Mv = Wy * Mv;
    end

    % Integration
    I = femMultiplyCell(Fxy,Mv);

%%% BOUNDARY ELEMENT INTEGRATION --> \int_{mesh(x)} psi(x)' f(x,y) dx
elseif isa(data{1},'dom') && isnumeric(data{2})
    % Domain with quadrature
    Xdom   = data{1};
    [X,Wx] = Xdom.qud;
    Nx     = size(X,1);
    Wx     = spdiags(Wx,0,Nx,Nx);

    % Evaluation points
    Y  = data{2};
    Ny = size(Y,1);

    % Integrated finite element matrix
    u  = data{3};
    Mu = u.uqm(Xdom);
    if iscell(Mu)
        Mu{1} = Mu{1}' * Wx;
        Mu{2} = Mu{2}' * Wx;
        Mu{3} = Mu{3}' * Wx;
    else
        Mu = Mu' * Wx;
    end

    % Function evaluation on quadrature
    F = data{4};
    if iscell(F)
        Fxy = {zeros(Nx,Ny),zeros(Nx,Ny),zeros(Nx,Ny)};
        for i = 1:3
            for j = 1:Ny
                Fxy{i}(:,j) = F{i}(X,Y(j,:));
            end
        end
    else
        Fxy = zeros(Nx,Ny);
        for j = 1:Ny
            Fxy(:,j) = F(X,Y(j,:));
        end
    end

    % Integration
    I = femMultiplyCell(Mu,Fxy);

else
    error('domIntegral4.m : unavailable case')
end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function test_suite = test_distancePoints
%TESTDISTANCEPOINTS  One-line description here, please.
%   output = test_distancePoints(input)
%
%   Example
%   testDistancePoints
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-04-22,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.
% Licensed under the terms of the LGPL, see the file "license.txt"

test_suite = functiontests(localfunctions);

function testSingleSingle(testCase) %#ok<*DEFNU>

pt1 = [10 10];
pt2 = [10 20];
pt3 = [20 20];

testCase.assertEqual(distancePoints(pt1, pt2), 10, 'AbsTol', .01);
testCase.assertEqual(distancePoints(pt2, pt3), 10, 'AbsTol', .01);
testCase.assertEqual(distancePoints(pt1, pt3), 10*sqrt(2), 'AbsTol', .01);

function testSingleSingleNorm1(testCase)
% test norm 1, equivalent to sum of absolute differences

pt1 = [10 10];
pt2 = [10 20];
pt3 = [20 20];

testCase.assertEqual(distancePoints(pt1, pt2, 1), 10, 'AbsTol', .01);
testCase.assertEqual(distancePoints(pt2, pt3, 1), 10, 'AbsTol', .01);
testCase.assertEqual(distancePoints(pt1, pt3, 1), 20, 'AbsTol', .01);

function testSingleSingleMaxNorm(testCase)

pt1 = [10 10];
pt2 = [10 20];
pt3 = [20 20];

testCase.assertEqual(distancePoints(pt1, pt2, inf), 10, 'AbsTol', .01);
testCase.assertEqual(distancePoints(pt2, pt3, inf), 10, 'AbsTol', .01);
testCase.assertEqual(distancePoints(pt1, pt3, inf), 10, 'AbsTol', .01);

function testSingleSingle3d(testCase)

pt1 = [10 10 10];
pt2 = [10 20 10];
pt3 = [20 20 10];
pt4 = [20 20 20];

testCase.assertEqual(distancePoints(pt1, pt2), 10, 'AbsTol', .01);
testCase.assertEqual(distancePoints(pt2, pt3), 10, 'AbsTol', .01);
testCase.assertEqual(distancePoints(pt1, pt3), 10*sqrt(2), 'AbsTol', .01);
testCase.assertEqual(distancePoints(pt1, pt4), 10*sqrt(3), 'AbsTol', .01);

function testSingleSingle3dNorm1(testCase)
% test norm 1, equivalent to sum of absolute differences

pt1 = [10 10 30];
pt2 = [10 20 30];
pt3 = [20 20 30];
pt4 = [20 20 40];

testCase.assertEqual(distancePoints(pt1, pt2, 1), 10, 'AbsTol', .01);
testCase.assertEqual(distancePoints(pt2, pt3, 1), 10, 'AbsTol', .01);
testCase.assertEqual(distancePoints(pt1, pt3, 1), 20, 'AbsTol', .01);
testCase.assertEqual(distancePoints(pt1, pt4, 1), 30, 'AbsTol', .01);

function testSingleSingleMaxNorm3d(testCase)

pt1 = [10 10 10];
pt2 = [10 20 10];
pt3 = [20 20 10];
pt4 = [20 20 20];

testCase.assertEqual(distancePoints(pt1, pt2, inf), 10, 'AbsTol', .01);
testCase.assertEqual(distancePoints(pt2, pt3, inf), 10, 'AbsTol', .01);
testCase.assertEqual(distancePoints(pt1, pt3, inf), 10, 'AbsTol', .01);
testCase.assertEqual(distancePoints(pt1, pt4, inf), 10, 'AbsTol', .01);

function testSingleArray(testCase)

pt1 = [10 10];
pt2 = [10 20];
pt3 = [20 20];

testCase.assertEqual(...
    distancePoints(pt1, [pt1; pt2; pt3]), ...
    [0 10 10*sqrt(2)], 'AbsTol', .01);

function testArrayArray(testCase)

pt1 = [10 10];
pt2 = [10 20];
pt3 = [20 20];
pt4 = [20 10];

array1 = [pt1;pt2;pt3];
array2 = [pt1;pt2;pt3;pt4];
res = [...
    0 10 10*sqrt(2) 10;...
    10 0 10 10*sqrt(2);...
    10*sqrt(2) 10 0 10];

testCase.assertEqual(distancePoints(array1, array2), res, 'AbsTol', .01);

function testArrayArrayDiag(testCase)

pt1 = [10 10];
pt2 = [10 20];
pt3 = [20 20];

array = [pt1;pt2;pt3];

testCase.assertEqual(...
    distancePoints(array, array, 'diag'), ...
    [0;0;0], 'AbsTol', .01);

function testArrayArray3dDiag(testCase)

pt1 = [10 10 30];
pt2 = [10 20 30];
pt3 = [10 20 40];

array1 = [pt1;pt2;pt3];
array2 = [pt2;pt3;pt1];

testCase.assertEqual(...
    distancePoints(array1, array2, 'diag'), ...
    [10;10;10*sqrt(2)], 'AbsTol', .01);

function testArrayArray3dNorm1Diag(testCase)

pt1 = [10 10 30];
pt2 = [10 20 30];
pt3 = [10 20 40];

array1 = [pt1;pt2;pt3];
array2 = [pt2;pt3;pt1];

testCase.assertEqual(...
    distancePoints(array1, array2, 1, 'diag'), ...
    [10;10;20], 'AbsTol', .01);

function testArrayArrayDiagMaxNorm(testCase)

pt1 = [10 10];
pt2 = [10 20];
pt3 = [20 20];

array1 = [pt1;pt2;pt3];
array2 = [pt2;pt3;pt1];

testCase.assertEqual(...
    distancePoints(array1, array2, inf, 'diag'), ...
    [10;10;10], 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function run_writelsl(varargin)
% Provide real-time BCI outputs to the lab streaming layer.
% run_writetcp(Model,SourceStream)
%
% This function runs in the background and processes data from some MATLAB stream (created with some
% other background input plugin, e.g., the BioSemi reader). The resulting estimates are offered to
% the lab streaming layer for use by other programs.
%
% In:
%   Model : A model data structure (as obtained from bci_train) based on which the predictor shall be
%           created; typically this is a model struct, but for convenience it can be a file name,
%           variable name in the base workspace, or a cell array of {file name, variable name} to
%           refer to a variable inside a .mat file. The model is not modified by this function.
%           (default: 'lastmodel')
%
%   SourceStreamNames : Optional names of stream data structures in the MATLAB base workspace to
%                       consider as possible data sources (previously created with onl_newstream);
%                       if a stream contains all channels that are needed by the predictor, or
%                       alternatively has the right number and type of channels it will be considered
%                       as a potential source stream unless ambiguous. (default: 'laststream')
%
%   LabStreamName : Name of the output stream in the lab streaming layer (default: 'BCI')
%
%   ChannelNames : labels for each output channel (if this is for classification and the output form
%                  is set to distribution, the channels correspond to the labels of the respective classes
%                  (default: {'class1','class2', ..., 'classN'} depending on the # of classes)
%
%   OutputForm : Output data form, see onl_predict (default: 'distribution')
%
%   UpdateFrequency : The rate at which new outputs will be computed. (default: 10)
%
%   PredictAt : Predict at markers. If nonempty, this is a cell array of online target markers relative
%               to which predictions shall be made. If empty, predictions are always made on the most
%               recently added sample. (default: {})
%
%   Verbose : Verbose output. If false, the console output of the online pipeline will be suppressed.
%             (default: false)
%
%   PredictorName : Name of the predictor to be created; a variable of this name will be created in
%                   the MATLAB base workspace to hold the predictor's state. If a variable with this
%                   name already exists it will be overridden. (default: 'lastpredictor')
%
% Examples:
%   % open a new BCILAB processing stream, using the previously learned predictive model 'mymodel',
%   % and reading from a previously opened input stream named 'mystream'. Name the stream 'BCI' and
%   % name the channels according to the classes whose probabilities they carry (here for some hypothetical
%   % error perception classifier)
%   run_writelsl('mymodel','mystream','BCI',{'error','no_error'})
%
%   % as before, but pass arguments by name
%   run_writelsl('Model','mymodel','SourceStream','mystream','ChannelNames',{'error','no_error'})
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2012-03-21

persistent lib;

% declare the name of this component (shown in the menu)
declare_properties('name','Lab streaming layer');

% define arguments
opts = arg_define(varargin, ...
    arg({'pred_model','Model'}, 'lastmodel', [], 'Predictive model. As obtained via bci_train or the Model Calibration dialog.','type','expression'), ...
    arg({'in_stream','SourceStreamNames','SourceStream'}, 'laststream',[],'Input Matlab stream name(s). Optional names of stream data structures in the MATLAB base workspace to consider as possible data sources (previously created with onl_newstream); if a stream contains all channels that are needed by the predictor, or alternatively has the right number and type of channels it will be considered as a potential source stream unless ambiguous.','typecheck',false,'shapecheck',false), ...
    arg({'out_stream','LabStreamName','Target'},'bci',[],'Name of the lab stream. This is the name under which the stream is provided to the lab streaming layer.'), ...
    arg({'channel_names','ChannelNames'},{'class1','class2'},[],'Output channel labels. These are the labels of the stream''s channels. In a typical classification setting each channel carries the probability for one of the possible classes.'), ...
    arg({'out_form','OutputForm','Form'},'expectation',{'expectation','distribution','mode','raw'},'Output form. Can be the expected value (posterior mean) of the target variable, or the distribution over possible target values (probabilities for each outcome, or parametric distribution), or the mode (most likely value) of the target variable.'), ...
    arg({'update_freq','UpdateFrequency'},10,[0 Inf],'Update frequency. The rate at which new outputs will be computed.'), ...
    arg({'predict_at','PredictAt'}, {},[],'Predict at markers. If nonempty, this is a cell array of online target markers relative to which predictions shall be made. If empty, predictions are always made on the most recently added sample.','type','expression'), ...
    arg({'verbose','Verbose'}, false,[],'Verbose output. If false, the console output of the online pipeline will be suppressed.'), ...
    arg({'source_id','SourceID'}, 'input_data',{'input_data','model'},'Use as source ID. This is the data that determines the source ID of the stream (if the stream is restarted, readers will continue reading from it if it has the same source ID). Can be input_data (use a hash of dataset ID + target markers used for training) or model (use all model parameters).'), ...
    arg({'pred_name','PredictorName'}, 'lastpredictor',[],'Name of new predictor. A variable of this name will be created in the MATLAB base workspace to hold the predictor''s state. If a variable with this name already exists it will be overridden.'));

% load the model
model = utl_loadmodel(opts.pred_model);

% check if channel labels make sense for the model
if strcmp(opts.out_form,'distribution')
    if isfield(model,'classes') && ~isempty(model.classes)
        if length(opts.channel_names) ~= length(model.classes)
            disp('The number of classes provided by the model does not match the number of provided channel names; falling back to default names.');
            opts.channel_names = cellfun(@(k)['class' num2str(k)],num2cell(1:length(model.classes),1),'UniformOutput',false);
        end
    end
else
    if isfield(model,'classes') && ~isempty(model.classes)
        if length(opts.channel_names) ~= 1
            disp('A classification model will produce just one channel if the output is not in distribution form, but a different number of channels was given. Falling back to the default channel label.');
            opts.channel_names = {'class'};
        end
    end
end

% try to calculate a UID for the stream
try
    if strcmp(opts.source_id,'model')
        uid = hlp_cryptohash({rmfield(model,'timestamp'),opts.predict_at,opts.in_stream,opts.out_stream});
    elseif strcmp(opts.source_id,'input_data')
        uid = hlp_cryptohash({model.source_data,opts.predict_at,opts.in_stream,opts.out_stream});
    else
        error('Unsupported SourceID option: %s',hlp_tostring(opts.source_id));
    end
catch e
    disp('Could not generate a unique ID for the predictive model; the BCI stream will not be recovered automatically after the provider system had a crash.');
    hlp_handleerror(e);
    uid = '';
end

% instantiate the library
disp('Opening the library...');
if isempty(lib)
    lib = lsl_loadlib(env_translatepath('dependencies:/liblsl-Matlab/bin')); end

% describe the stream
disp('Creating a new streaminfo...');
info = lsl_streaminfo(lib,opts.out_stream,'MentalState',length(opts.channel_names),opts.update_freq,'cf_float32',uid);
% ... including some meta-data
desc = info.desc();
channels = desc.append_child('channels');
for c=1:length(opts.channel_names)
    newchn = channels.append_child('channel');
    newchn.append_child_value('name',opts.channel_names{c});
    newchn.append_child_value('type',opts.out_form);
end

% create an outlet
outlet = lsl_outlet(info);

% start background writer job
onl_write_background( ...
    'ResultWriter',@(y)send_samples(outlet,y),...
    'MatlabStream',opts.in_stream, ...
    'Model',opts.pred_model, ...
    'OutputFormat',opts.out_form, ...
    'UpdateFrequency',opts.update_freq, ...
    'PredictorName',opts.pred_name, ...
    'PredictAt',opts.predict_at, ...
    'Verbose',opts.verbose, ...
    'StartDelay',0,...
    'EmptyResultValue',[]);

disp('Now writing...');

function send_samples(outlet,y)
if ~isempty(y)
    outlet.push_chunk(y'); end

repo: REST


text:
function engine = bk_inf_engine(bnet, varargin)
% BK_INF_ENGINE Boyen-Koller approximate inference algorithm for DBNs.
%
% In the BK algorithm, the belief state is represented as a product of marginals,
% even though the factors may not be independent.
%
% engine = bk_inf_engine(bnet, ...)
%
% The following optional arguments can be specified in the form of name/value pairs:
% [default value in brackets]
%
% clusters - if a cell array, clusters{i} specifies the terms in the i'th factor.
%          - 'exact' means create one cluster that contains all the nodes in a slice [exact]
%          - 'ff' means create one cluster per node (ff = fully factorised).
%
%
% For details, see
% - "Tractable Inference for Complex Stochastic Processes", X. Boyen and D. Koller, UAI 98.
% - "Approximate learning of dynamic models",  X. Boyen and D. Koller, NIPS 98.
% (The UAI98 paper discusses filtering and theory, and the NIPS98 paper discusses smoothing.)

ss = length(bnet.intra);
% set default params
clusters = 'exact';

if nargin >= 2
  args = varargin;
  nargs = length(args);
  for i=1:2:nargs
    switch args{i},
     case 'clusters',  clusters = args{i+1};
     otherwise, error(['unrecognized argument ' args{i}])
    end
  end
end

if strcmp(clusters, 'exact')
  %clusters = { compute_interface_nodes(bnet.intra, bnet.inter) };
  clusters = { 1:ss };
elseif strcmp(clusters, 'ff')
  clusters = num2cell(1:ss);
end

% We need to insert the prior on the clusters in slice 1,
% and extract the posterior on the clusters in slice 2.
C = length(clusters);
clusters2 = cell(1,2*C);
clusters2(1:C) = clusters;
for c=1:C
  clusters2{c+C} = clusters{c} + ss;
end

onodes = bnet.observed;
obs_nodes = [onodes(:) onodes(:)+ss];
engine.sub_engine = jtree_inf_engine(bnet, 'clusters', clusters2);

engine.clq_ass_to_cluster = zeros(C, 2);
for c=1:C
  engine.clq_ass_to_cluster(c,1) = clq_containing_nodes(engine.sub_engine, clusters{c});
  engine.clq_ass_to_cluster(c,2) = clq_containing_nodes(engine.sub_engine, clusters{c}+ss);
end
engine.clusters = clusters;

engine.clq_ass_to_node = zeros(ss, 2);
for i=1:ss
  engine.clq_ass_to_node(i, 1) = clq_containing_nodes(engine.sub_engine, i);
  engine.clq_ass_to_node(i, 2) = clq_containing_nodes(engine.sub_engine, i+ss);
end

% Also create an engine just for slice 1
bnet1 = mk_bnet(bnet.intra1, bnet.node_sizes_slice, 'discrete', myintersect(bnet.dnodes, 1:ss), ...
		'equiv_class', bnet.equiv_class(:,1), 'observed', onodes);
for i=1:max(bnet1.equiv_class)
  bnet1.CPD{i} = bnet.CPD{i};
end

engine.sub_engine1 = jtree_inf_engine(bnet1, 'clusters', clusters);

engine.clq_ass_to_cluster1 = zeros(1,C);
for c=1:C
  engine.clq_ass_to_cluster1(c) = clq_containing_nodes(engine.sub_engine1, clusters{c});
end

engine.clq_ass_to_node1 = zeros(1, ss);
for i=1:ss
  engine.clq_ass_to_node1(i) = clq_containing_nodes(engine.sub_engine1, i);
end

engine.clpot = []; % this is where we store the results between enter_evidence and marginal_nodes
engine.filter = [];
engine.maximize = [];
engine.T = [];

engine.bel = [];
engine.bel_clpot = [];
engine.slice1 = [];
%engine.pot_type = 'cg';
% hack for online inference so we can cope with hidden Gaussians and discrete
% it will not affect the pot type used in enter_evidence
engine.pot_type = determine_pot_type(bnet, onodes);

engine = class(engine, 'bk_inf_engine', inf_engine(bnet));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function f = make_gaussian_filter(sigma, cutoff)
%MAKE_GAUSSIAN_FILTER   Create Gaussian filter
%   MAKE_GAUSSSIAN_FILTER(SIGMA[, CUTOFF]) creates a Gaussian filter with
%   standard deviation SIGMA.  The filter mask is cut off at CUTOFF times
%   the standard deviation.  If CUTOFF is omitted, a default value of 3
%   is assumed.
%
%   See also FSPECIAL.
%
%   Author:  Stefan Roth, Department of Computer Science, TU Darmstadt
%   Contact: sroth@cs.tu-darmstadt.de
%   $Date$
%   $Revision$

% Copyright 2004-2007, Brown University, Providence, RI. USA
% Copyright 2007-2010 TU Darmstadt, Darmstadt, Germany.
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

  if (nargin < 2)
    cutoff = 3;
  end

  sz = ceil(cutoff * sigma);

  f = exp(-(-sz:sz).^2 / (2*sigma^2));
  f = f ./ sum(f);

author: kristinbranson

repo: JAABA


text:
function regresdiagplot(xdist,ydist,cutoffx,cutoffy,k,multi,attrib,labsd,labresd)

%REGRESDIAGPLOT makes a regression outlier map.
% The score distances (SD)/mahalanobis distances (MD)/robust distances (RD) and
% residual distances (ResD) calculated in a PCR or PLS analysis, or in a
% univariate or multivariate regression are plotted on the x and y-axis, respectively.
% The cutoff values marked by a red line indicate which observations are
% outlying with respect to the majority of the data. The terminology used is
% as indicated in the tabular:
%                         small SD/MD/RD  | large SD/MD/RD
%      ----------------------------------------------------------
%      small ResD/StdRes | regular point      | good leverage point
%      large ResD/StdRes | vertical outlier   | bad leverage point
%
% For more details see:
%    Hubert, M., Verboven, S. (2003),
%    "A robust PCR method for high-dimensional regressors,"
%    Journal of Chemometrics, 17, 438-452.
%
% Required input arguments:
%
%  xdist   : Score distances from a PCA analysis or
%            Mahalanobis distances from LS or MLR analysis or
%            Robust distances from LTS or MCDREG analysis,
%  ydist   : Residual distances of the PCR method or
%            Standardized residuals of one of the regression methods,
%  cutoffx : Cutoff value for the SD/MD/RD,
%  cutoffy : Cutoff value for the residual-distance,
%  k       : Number of  principal components used in the PCR/PLS method, or zero
%            if not available.
%  multi   : 0 for univariate analysis, 1 for multivariate analysis
%  attrib  : String identifying the method used =  'LS', 'MLR', 'LTS', 'MCDREG', 'RPCR',
%           'CPCR', 'RSIMPLS', 'CSIMPLS'
%
% Optional inputs:
%   labsd    : number of displayed points with largest distance on x-axes (default = 3)
%   labresd  : number of displayed points with largest distance on y-axes (default = 3)
%
%
% I/O: regresdiagplot(out.sd,out.rd,out.cutoff.sd,out.cutoff.rd,k,multi,attrib,labsd,labresd)
%
% Example: regresdiagplot(out.sd,out.rd,out.cutoff.sd,out.cutoff.rd,out.k,1,'RPCR',5,5)
%          regresdiagplot(out.md,out.stdres,out.cutoffmd,0,0,0,'LS',5,4)
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by Sabine Verboven
% Created 25 February 2002
% Revised : 12/02/2004
%

if nargin==8
    labresd=3;
end
if nargin==7
    labsd=3;
    labresd=3;
end
if nargin<7
    error('A required input variable is missing!')
end

if ischar(xdist) | ischar(ydist)
    mess=sprintf(['Warning: A singularity was detected in the analysis. ',...
            '\n No regression outlier map available.']);
    disp(mess)
    return
end
%all LTS-analysis in RPCR are intercept included!!!
if multi==1 %multivariate analysis
    x=xdist; %distances on x-axis: SD, MD or RD
    y=ydist; %residual distances
    quanty=cutoffy;
    quantx=cutoffx;
else  %univariate analysis
    quantx=cutoffx;
    quanty=sqrt(chi2inv(0.975,1));
    x=xdist;   %distances on x-axis: SD, MD, RD
    y=ydist;   %standardized residuals!!!
end

plot(x,y,'o');
set(gcf,'Name', 'Regression outlier map', 'NumberTitle', 'off');
hold on
xmin=0;%max(0,min(x));%-0.5;
xmax=max([quantx max(x)])+0.5;

if multi==1
    ymin=0;%max(0,min(y));%-1;
    if size(y,2)==1
        y=y';
    end
    ymax=max([y,quanty])+1;
    xlim([xmin,xmax]);
    ylim([ymin,ymax]);
    ylimi=get(gca,'Ylim');
    line([quantx,quantx],[ylimi(1),ylimi(2)],'Linestyle','-','Color','r');
    xlimi=get(gca,'Xlim');
    line([xlimi(1),xlimi(2)],[quanty,quanty],'Linestyle','-','Color','r');
    if labsd
        plotnumbers(x',y',labsd,labresd,2)%sort by  x and  y
    end
    ylabel('Residual distance')
else
    ymin=min([-4 min(y)])-0.5;
    ymax=max([4 max(y)])+0.5;
    xlim([xmin,xmax]);
    ylim([ymin,ymax]);
    ylimi=get(gca,'Ylim');
    ylimi=[min(ylimi(1),min([-4 min(y)]))-0.5,max(ylimi(2),max([4 max(y)]))+0.5];
    line([quantx,quantx],[ylimi(1),ylimi(2)],'Linestyle','-','Color','r');
    xlimi=get(gca,'Xlim');
    xlimi=[min(xlimi(1),min(x))-1,max(xlimi(2),max([quantx max(x)]))+1];
    line([xlimi(1),xlimi(2)],[quanty,quanty],'Linestyle','-','Color','r');
    line([xlimi(1),xlimi(2)],[-quanty,-quanty],'Linestyle','-','Color','r');
    if labsd
        plotnumbers(x',y',labsd,labresd,3)%sort by x and abs(y)
    end
    ylabel('Standardized residual')
end

box on

switch attrib
    case 'LS'
        xlabel('Mahalanobis distance');
        ylabel('Standardized LS residual');
    case 'MLR'
        xlabel('Mahalanobis distance');
    case 'LTS'
        xlabel('Robust distance computed by MCD');
        ylabel('Standardized LTS residual');
    case 'MCDREG'
        xlabel('Robust distance computed by MCD');
    case {'CPCR'}
        xlabel(['Score distance (',num2str(k),' LV)']);
    case {'RPCR'}
        xlabel(['Score distance (',num2str(k),' LV)']);
    case {'CSIMPLS'}
        xlabel(['Score distance (',num2str(k),' LV)']);
    case {'RSIMPLS'}
        xlabel(['Score distance (',num2str(k),' LV)']);

end
title(attrib)
hold off

author: marianux

repo: ecg-kit

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


text:
function title = p01_title ( title )

%*****************************************************************************80
%
%% P01_TITLE returns the title of problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'F(X) = SIN(X) - 0.5 * X';

  return
end


text:
function results = classify(data, net, labels, debug);
% CLASSIFY Classifies the given data using the given trained SFAM.
% RESULTS = CLASSIFY(DATA, NET, LABELS, DEBUG)
%	DATA is an M-by-D matrix where M is the number of samples and D is the size of the feature
%	space. NET is a previously trained SFAM network. LABELS is a M-vector containing the correct
%	labels for the data. If you don't have them, give it as an empty-vector [].
%	DEBUG is a scalar to control the verbosity of the program during training. If 0, nothing will
%	be printed, otherwise every DEBUG iterations an informatory line will be printed.
%
% Emre Akbas, May 2006
%

results = [];
hits=0;

tic;
for s=1:size(data,1)
    input = data(s,:);

    % Complement code input
    input = [input 1-input];

    % Compute the activation values for each prototype.
    activation = ones(1,length(net.weights));
    for i=1:length(net.weights)
	activation(i)  = sum(min(input,net.weights{i}))/...
		    (net.alpha + sum(net.weights{i}));
    end

    % Sort activation values
    [sortedActivations, sortedIndices] = sort(activation,'descend');

    % Iterate over the prototypes with decreasing activation-value
    results(s)=-1;
    for p=sortedIndices
	% Compute match of the current candidate prototype
	match = sum(min(input,net.weights{p}))/net.D;

	% Check resonance
	if match>=net.vigilance
	    results(s) = net.labels(p);

	    if ~isempty(labels)
		if labels(s)==results(s), hits = hits + 1; end;
	    end

	    break;
	end
    end

    if mod(s,debug)==0
	elapsed = toc;
	fprintf(1,'Tested %4dth sample. Hits so far: %3d which is %.3f%%.\tElapsed %.2f seconds.\n',s,hits,100*hits/s,elapsed);
	tic;
    end
end % samples loop

author: Sable

repo: mcbench-benchmarks


text:
function linpack_s_test11 ( )

%*****************************************************************************80
%
%% TEST11 tests SGEFA and SGESL.
%
%  Discussion:
%
%    In this example, we solve a relatively large linear system.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 100;
  lda = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST11\n' );
  fprintf ( 1, '  For a general matrix,\n' );
  fprintf ( 1, '  SGEFA computes the LU factors;\n' );
  fprintf ( 1, '  SGESL solves a factored linear system;\n' );
  fprintf ( 1, '  The matrix size is N = %d\n', n );
%
%  Assign values to the matrix A and the right hand side B.
%
%  The problem is just an enlarged version of the
%  problem for N = 5, which is:
%
%  Matrix A is ( n -1 -1 -1 -1)    Right hand side B is  (1)
%              (-1  n -1 -1 -1)                          (1)
%              (-1 -1  n -1 -1)                          (1)
%              (-1 -1 -1  n -1)                          (1)
%              (-1 -1 -1 -1  n)                          (1)
%
%  Solution is   (1)
%                (1)
%                (1)
%                (1)
%                (1)
%
  b(1:n) = 1.0;

  a(1:n,1:n) = -1.0;
  for i = 1 : n
    a(i,i) = n;
  end
%
%  Factor the matrix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Factor the matrix\n' );

  [ a, ipvt, info ] = sgefa ( a, lda, n );

  if ( info ~= 0 )
    fprintf ( 1, '  SGEFA returned an error flag INFO = %d\n', info );
    return
  end
%
%  Solve the linear system.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve the linear system.\n' );

  job = 0;
  b = sgesl ( a, lda, n, ipvt, b, job );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The first and last five entries of the solution:\n' );
  fprintf ( 1, '  (All of them should be 1.)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    if ( i <= 5 | n-5 < i )
      fprintf ( 1, '  %6d  %14f\n', i, b(i) );
    end
    if ( i == 5 )
      fprintf ( 1, '  ......  ..............\n' );
    end
  end

  return
end


text:
function [F_RF,F_BB]=OMP_Precoding()
global Ns Nt Nrf H Codebook_v
F_RF = [];
[~,~,V] = svd(H);
F_opt = V(:,1:Ns);
Fres = F_opt;
for i = 1:Nrf
    y = Codebook_v'*Fres;
    k = find(diag(y*y')==max(diag(y*y')));
    F_RF (:,i) = Codebook_v(:,k);
    F_BB = (F_RF' * F_RF)^(-1) * F_RF' *F_opt;
    Fres = (F_opt - F_RF * F_BB)/norm(F_opt - F_RF * F_BB,'fro');
end
F_BB = sqrt(Ns) * F_BB / norm(F_RF * F_BB,'fro');

author: Zzhaoxingyu

repo: hybrid-beamforming-for-three-scenes


text:
function [y]=attachToolToFlange(t_kuka,ts)
%% About
% Attach the tool specified by the transform (ts) into the flange

%% Areguments
% t_kuka: TCP/IP comunication object
% ts: transform vector [x,y,z,a,b,c] of TCP to the flange of the robot
% x,y,z are in mm
% a,b,c are XYZ fixed rotation angles in radians

% Copyright: Mohammad SAFEEA, 2018

    if max(max(size(ts==6)))
    else
        y=false;
        return;
    end
%%
    theCommand=['TFtrans']; % instruction part
    for i=1:6
        st=num2str(ts(i));
        theCommand=[theCommand,'_',st];
    end
    fprintf(t_kuka, theCommand);
    message=fgets(t_kuka);
    [ret]=checkAcknowledgment(message);
    y=ret;
    if(y==false)
        % turn off the server
        disp('Could not attach the tool to the Robot');
        net_turnOffServer(t_kuka);
    end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function s = getBehAssociation(Gamma,y,T,synch)
% Returns, for each time point, how well the state time courses predict
% the variable y across trials, which must have dimension (trials by 1)
%
% Author: Diego Vidaurre, OHBA, University of Oxford (2018)

if ~all(T(1)==T)
    error('Synchronisity can only be measured if trials have the same length');
end
N = length(T);
ttrial = size(Gamma,1)/N;
if nargin<4, synch = zeros(1,T(1)); end
K = size(Gamma,2);
Gamma = permute(reshape(Gamma,[ttrial N K]),[2 3 1]);
s = zeros(1,ttrial);
these = ~isnan(y);
y = y(these);
N = length(y);
y = zscore(y);
for t = 1:ttrial
   if synch(t)==1, continue; end
   x = [ones(N,1) Gamma(these,1:K-1,t)];
   beta = (x' * x + 1e-8 * eye(K) ) \ (x' *  y);
   s(t) = 1 - sum((y - x * beta).^2) / sum(y.^2);
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function gammaM = gammaScan3d(doseArray1, doseArray2, deltaXYZv, doseAgreement, distAgreement, maxDistance)
% function gammaM = gammaScan3d(doseArray1, doseArray2, deltaXYZv, doseAgreement, distAgreement, maxDistance)
%
% APA, 04/27/2012

deltaX = deltaXYZv(1);
deltaY = deltaXYZv(1);
deltaZ = deltaXYZv(1);
incrementRadius = min([deltaX deltaY deltaZ]);

% Initial gamma
gammaM = ((doseArray1-doseArray2).^2).^0.5/doseAgreement;
convergedM = false(size(gammaM));

% Calculate until 4 times the permissible distance to agreement.
if ~exist('maxDistance', 'var')
    maxDistance = distAgreement*1.5;
end
outerRadiusV = incrementRadius:incrementRadius:maxDistance;
numIters = length(outerRadiusV);

siz = size(gammaM);
minDoseDiffM = zeros(siz,'single');
maxDoseDiffM = minDoseDiffM;

% convergenceCountM = zeros(siz,'uint8');

% Update waitbar on gamma GUI
gammaGUIFig = findobj('tag','CERRgammaInputGUI');
if ~isempty(gammaGUIFig)
    ud = get(gammaGUIFig,'userdata');
    set(ud.wb.patch,'xData',[0 0 0 0])
end

for radNum = 1:numIters

    disp(['--- Gamma Calculation Iteraion ', num2str(radNum), ' ----'])

    % Create an ellipsoid ring neighborhood
    outerRadius = outerRadiusV(radNum);

    rowOutV = floor(outerRadius/deltaY);
    colOutV = floor(outerRadius/deltaX);
    slcOutV = floor(outerRadius/deltaZ);

    NHOOD_outer = createEllipsoidNHOOD(1:rowOutV,1:colOutV,1:slcOutV);

    % Compute Min and Max for the ellipsoid neighborhood
    [minLocalM, maxLocalM] = getMinMaxIM(doseArray2,NHOOD_outer);

    % Compute difference between dose1 and (min,max) for dose2
    minDoseDiffM(~convergedM) = doseArray1(~convergedM) - minLocalM(~convergedM) + doseAgreement; % D1-(min(D2)-5)
    maxDoseDiffM(~convergedM) = maxLocalM(~convergedM) + doseAgreement - doseArray1(~convergedM); % (max(D2)+5)-D1

    % If dose1 is contained within (min,max) for dose2, then it converged!
    newConvergedM = ~convergedM & minDoseDiffM >= 0 & maxDoseDiffM >= 0;

    % Compute gamma for voxels that converged
    %gammaM(newConvergedM) =  min(gammaM(newConvergedM), outerRadius/distAgreement);
    gammaM(newConvergedM) =  outerRadius;

    % Add newly converged voxels to the list
    convergedM = convergedM | newConvergedM;

    % Compute gamma for voxels that have not yet converged
    %gammaForNotConvergedM =  (min((minDoseDiffM(~convergedM)-doseAgreement).^2, (maxDoseDiffM(~convergedM)-doseAgreement).^2)./doseAgreement^2 + (outerRadius/distAgreement)^2).^0.5;
    %gammaM(~convergedM) = min(gammaM(~convergedM), gammaForNotConvergedM);

    % If gamma is less or equal to outerRadius/distAgreement, then the voxel conveged!
    %convergedM(~convergedM) = gammaM(~convergedM) <= outerRadius/distAgreement;

    if ~isempty(gammaGUIFig)
        set(ud.wb.patch,'xData',[0 0 radNum/numIters radNum/numIters])
        drawnow
    end

    if all(convergedM(:))
        disp('All converged!!!')
        if ~isempty(gammaGUIFig)
            set(ud.wb.patch,'xData',[0 0 1 1])
        end
        break
    end

end

gammaM(~convergedM) = NaN;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ z, stats ] = solve_cosamp( Phi, K, y, x)
%SOLVEBPIC Solves the CoSaMP problem
    [M, N] = size(Phi);
    epsilon = 1e-6;
    %We now setup our recovery program.
    solver = spx.pursuit.single.CoSaMP(Phi, K);
    result = solver.solve(y);
    % recovered representation. N length vector
    z = result.z;
    stats.recoveredRepresentationVector  = z;
    % recovery error vector. N length vector
    h = x - z;
    stats.recoveryErrorVector = h;
    % l_2 norm of reconstruction error
    stats.recoveryErrorL2Norm = norm(h);
    % The K non-zero coefficients in x (set of indices)
    stats.T0 = find (x ~= 0);
    % The portion of recovery error over T0 K length vector
    stats.recoveryErrorVectorT0 = h(stats.T0);
    % Positions of other places (set of indices)
    stats.T0C = setdiff(1:N , stats.T0);
    % Recovery error at T0C places [N - K] length vector
    hT0C = h(stats.T0C);
    stats.recoveryErrorVectorT0C = hT0C;
    % The K largest indices after T0 in recovery error (set of indices)
    stats.T1 = spx.commons.signals.largest_indices(hT0C, K);
    % The recovery error component over T1. [K] length vector.
    hT1 = h(stats.T1);
    stats.recoveryErrorVectorT1 = hT1;
    % Remaining indices [N - 2K] set of indices
    stats.TRest = setdiff(stats.T0C , stats.T1);
    % Recovery error over remaining indices [N - 2K] length vector
    hTRest = h(stats.TRest);
    stats.recoveryErrorVectorTRest = hTRest;
    % largest indices of the recovered vector
    stats.TT0 = spx.commons.signals.largest_indices(z, K);
    % Support Overlap
    stats.supportOverlap = intersect(stats.T0, stats.TT0);
    % Support recovery ratio
    stats.supportRecoveryRatio = numel(stats.supportOverlap) / K;
    % measurement error vector [M] length vector
    e = y - Phi * z;
    stats.measurementError = e;
    % Norm of measurement error.  This must be less than epsilon
    stats.measurementErrorL2Norm = norm(e);
    % Ratio between the norm of recovery error and measurement error
    stats.recoveryToMeasurementErrorNormRatio = ...
        stats.recoveryErrorL2Norm / stats.measurementErrorL2Norm;
    % whether we consider the process to be success or not.
    % We consider success only if the support has been recovered
    % completely.
    stats.success = numel(stats.supportOverlap) == K;
end

author: indigits

repo: sparse-plex


text:
function gui_check_range(edithandle,resetstr)
% check whether the the given edit field contains a range; warn and reset if not
if ~exist('resetstr','var')
    resetstr = ''; end

str = get(edithandle,'String');
if ~isempty(str) && numel(str2num(str)) ~= 2
    errordlg2('This field must be of the form [number, number]');
    set(edithandle,'String',resetstr);
end


text:
function pde = sincosdata3
%% SINCOSDATA3 trigonometric data for Poisson equation in 3-D
%
%     u = cos(pi*x)*cos(pi*y)*cos(pi*z);
%     f = - Delta u = 3*pi^2*cos(pi*x)*cos(pi*y)*cos(pi*z);
%     Du = (-pi*sin(pi*x)*cos(pi*y)*cos(pi*z),
%           -pi*cos(pi*x)*sin(pi*y)*cos(pi*z),
%           -pi*cos(pi*x)*cos(pi*y)*sin(pi*z));
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

pde = struct('f',@f,'exactu',@exactu,'g_D',@g_D,'Du',@Du,'g_N',@g_N, 'phi', @phi);

    % load data (right hand side function)
    function s = f(p)
    x = p(:,1); y = p(:,2); z = p(:,3);
    s = 3*pi^2*cos(pi*x).*cos(pi*y).*cos(pi*z);
    end
    % exact solution
    function s = exactu(p)
    x = p(:,1); y = p(:,2); z = p(:,3);
    s = cos(pi*x).*cos(pi*y).*cos(pi*z); % for neumann boundary condition, int_u =0
    end
    % Dirichlet boundary condition
    function s = g_D(p)
    s = exactu(p);
    end
    % Neumann boundary condtion
    function f = g_N(p)
    f = zeros(size(p,1),1);
    x = p(:,1); y = p(:,2); z = p(:,3);
    uprime = [-pi*sin(pi*x).*cos(pi*y).*cos(pi*z) ...
              -pi*cos(pi*x).*sin(pi*y).*cos(pi*z) ...
              -pi*cos(pi*x).*cos(pi*y).*sin(pi*z)];
    downbd = (abs(z)<eps); % n = (0,0,-1)
    f(downbd) = -uprime(downbd,3);
    upbd = (abs(z-1)<eps);% n = (0,0,1)
    f(upbd) = uprime(upbd,3);
    leftbd = (abs(x)<eps); % n = (-1,0,0)
    f(leftbd) = -uprime(leftbd,1);
    rightbd = (abs(x-1)<eps); % n = (1,0,0)
    f(rightbd) = uprime(rightbd,1);
    backbd = (abs(y-1)<eps); % n =(0,1,0)
    f(backbd) = uprime(backbd,2);
    frontbd = (abs(y)<eps); % n = (0,-1,0)
    f(frontbd) = -uprime(frontbd,2);
    end
    % Derivative of the exact solution
    function s = Du(p)
    x = p(:,1); y = p(:,2); z = p(:,3);
    s(:,1) = -pi*sin(pi*x).*cos(pi*y).*cos(pi*z);
    s(:,2) = -pi*cos(pi*x).*sin(pi*y).*cos(pi*z);
    s(:,3) = -pi*cos(pi*x).*cos(pi*y).*sin(pi*z);
    end
    function s = phi(p) % level set function
        x = p(:,1); y = p(:,2); z =p(:,3);
        s = x.^2 + y.^2 +z.^2 - 0.75^2;
    end
end

repo: ifem


function res=SODW(x,a,b,w);
% function res=SODW(x,a,b,w);
%
% Computes and sums over all outer products of the columns in x and
% weights them according to w.
%
% equivalent to:
%
% res=zeros(size(x,1));
% for i=1:n
%   res=res+w(i).*x(:,a(i))*x(:,b(i))';
% end;
%
%
% copyright 2005 by Kilian Q. Weinberger
% University of Pennsylvania
% kilianw@seas.upenn.edu
% ********************************************

if(min(w)<0) error('Weights must be non-negative in matlab version!\nPlease call "mex SODW.cpp" in the mexfunctions directory.\n');end;
[D,N]=size(x);
B=round(2500/D^2*1000000);
res=zeros(D^2,1);
sw=sqrt(w);
for i=1:B:length(a)
  BB=min(B,length(a)-i);
  Xa=mulh(x(:,a(i:i+BB)),sw(i:i+BB));
  Xb=mulh(x(:,b(i:i+BB)),sw(i:i+BB));
  XaXb=Xa*Xb';
  res=res+vec(Xa*Xa'+Xb*Xb'-XaXb-XaXb');

  if(i>1)   fprintf('.');end;
end;


text:
%kkexprn 'Output = Function(Input1, Input2)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kexprn.pane file
%
% Parameters:
% InputFile: i1 ' i1 ', required: 'input data object 1'
% InputFile: i2 ' i2 ', optional: 'input data object 2'
% InputFile: i3 ' i3 ', optional: 'input data object 3'
% InputFile: i4 ' i4 ', optional: 'input data object 4'
% InputFile: i5 ' i5 ', optional: 'input data object 5'
% InputFile: i6 ' i6 ', optional: 'input data object 6'
% String: expr ' ', default: ' ': 'arithmetic expression'
% OutputFile: o ' o ', required: 'resulting output data object'
%
% Example: o = kkexprn({i1, i2, i3, i4, i5, i6}, {'i1','';'i2','';'i3','';'i4','';'i5','';'i6','';'expr',' ';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kexprn - Output = Function(Input1, Input2)
%
%  DESCRIPTION
% .I kgexpr
% allows you to generate
% 1D, 2D, 3D, 4D, or 5D datasets of any data type.  For dataset generation,
% you must specify double precision "begin" and \fIend\fP values, as well
% as an integer "size" value, for each of width, height, depth, time, and
% elements, as relevant to the desired dimensionality of the dataset to be
% created.
% .PP
% The function that will define the dataset is specified as:
% f(w,h,d,t,e) = <expression>
% Where "w, h, d, t," and \fIe\fP are double precision values ranging from
% "begin" to \fIend\fP with the number of values to be generated determined by\fIsize\fP.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkexprn(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,..] = kkexprn(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'i2', '__input';'i3', '__input';'i4', '__input';'i5', '__input';'i6', '__input';'expr', ' ';'o', '__output'};
maxval={0,1,1,1,1,1,0,0};
minval={0,1,1,1,1,1,0,0};
istoggle=[0,1,1,1,1,1,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','InputFile','InputFile','InputFile','InputFile','String','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 'kexprn"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB