text:
% MAIN  --  Quad-Rotor  --  Minimal-Acceleration trajectory
%
% Fin the minimal acceleration-squared trajectory to move the quad-rotor from an
% arbitrary state to the origin.
%
% NOTES:
%   X = [x;y;q] = [x pos, y pos, angle] = configuration
%  dX = [dx;dy;dq] = [x vel, y vel, angle rate] = rate
% ddX = [ddx;ddy;ddq] = acceleration
%
% PROBLEM:
%
% ddX = f(X,dX,u);     <-- dynamics
%
% cost = integral(  ddX^2  );     <-- cost function
%
% subject to:
%   X(0) = X0;
%   X(1) = XF;
%   dX(0) = dX0;
%   dX(1) = dXF;
%
% How to pose as a standard trajectory optimization problem?
%
% dX = V1;
% dV1 = f(X,V1,U1)
%
% V2 == V1;   % <-- Key line.
% dV2 = U2;
% cost = integral(  U2^2  );
%
% z = [X;V1;V2]
% u = [U1;U2]
%

clc; clear;

addpath ../../

% Dynamics paramters
p.g = 9.81; % (m/s^2) gravity
p.d = 0.3;  % (m) half-width of quad rotor
p.m = 0.2;  % (m) half-mass of the quad rotor

% Trajectory Parameters:
duration = 1;
uMax = 5*p.g*p.m;

% Initial State:
X0 = [1;0;0];   %  initial configuration
dX0 = zeros(3,1);  % initial rates
z0 = [X0; dX0; dX0];  % initial state

XF = [0;0;0];   % final configuration
dXF = zeros(3,1);  % final rates
zF = [XF; dXF; dXF];  % final state

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                     Set up function handles                             %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
w = 1./[1,1,1];  %weighting vector for path objective

problem.func.dynamics = @(t,z,u)( dynAcc(z,u,p) );
problem.func.pathObj = @(t,z,u)( pathObj(u,w) );  %accel-squared cost function
problem.func.pathCst = @(t,z,u)( pathCst(z) );

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                     Set up problem bounds                               %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

problem.bounds.initialTime.low = 0;
problem.bounds.initialTime.upp = 0;
problem.bounds.finalTime.low = duration;
problem.bounds.finalTime.upp = duration;

problem.bounds.initialState.low = z0;
problem.bounds.initialState.upp = z0;
problem.bounds.finalState.low = zF;
problem.bounds.finalState.upp = zF;

problem.bounds.control.low = [-uMax*[1;1];  -inf(3,1)];   %[torque, accel]
problem.bounds.control.upp = [uMax*[1;1];  inf(3,1)];

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                    Initial guess at trajectory                          %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

problem.guess.time = [0,duration];
problem.guess.state = [z0, zF];
problem.guess.control = [p.g*p.m*ones(2,2); zeros(3,2)];

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                         Solver options                                  %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

problem.options.nlpOpt = optimset(...
    'Display','iter',...
    'MaxFunEvals',1e5);

% problem.options.method = 'trapezoid';
% problem.options.trapezoid.nGrid = 40;
%
% problem.options.method = 'hermiteSimpson';
% problem.options.hermiteSimpson.nSegment = 10;

problem.options.method = 'chebyshev';
problem.options.chebyshev.nColPts = 15;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                            Solve!                                       %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

soln = optimTraj(problem);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                        Display Solution                                 %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

%%%% Unpack the simulation
t = linspace(soln(end).grid.time(1), soln(end).grid.time(end), 150);

z = soln(end).interp.state(t);
x = z(1,:);
y = z(2,:);
q = z(3,:);
dx = z(4,:);
dy = z(5,:);
dq = z(6,:);

X1 = z(1:3,:);
V1 = z(4:6,:);
V2 = z(7:9,:);

u = soln(end).interp.control(t);
u1 = u(1,:);
u2 = u(2,:);
A2 = u(3:5,:);

[dObj,uStar] = pathObj(u,w);

%%%% Plots:

%%%% Plots:

figure(2); clf;

subplot(2,2,1)
plot(t,X1);
legend('x','y','q')
title('configuration')

subplot(2,2,3)
plot(t,V2);
legend('x','y','q')
title('rates')

subplot(2,2,2)
plot(t,A2);
legend('x','y','q')
title('acceleration')

subplot(2,2,4); hold on;
plot(t,u1);  plot(t,u2);
title('actuators')
legend('u1','u2');

% Configuration trajectories
figure(1); clf;

subplot(2,2,1); hold on;
plot(t,x);
xlabel('t')
ylabel('x')
title('Minimum acceleration-squared trajectory')

subplot(2,2,2); hold on;
plot(t,y);
xlabel('t')
ylabel('y')

subplot(2,2,3); hold on;
plot(t,q);
xlabel('t')
ylabel('q')

subplot(2,2,4); hold on;
plot(t,u1);  plot(t,u2);
xlabel('t')
ylabel('u')
legend('u1','u2');

author: MatthewPeterKelly

repo: OptimTraj


text:
function [W, p, t] = build4(im, angles)
t_s = toc;
%Pad image
[im_pad,D] = impad(im,'diag',0);
sz = size(im);
m = sz(1); %rows
n = sz(2); %cols
n_proj = length(angles);
N = m*n; %number of unknown variables
M = D*n_proj; %number of equations
W_c = cell(n_proj,1);
%proj_info = zeros(M,2); %[theta No.]
p = zeros(M,1); %projection vector
ix = false(M,1);
rc = [m;n]/2+0.5;
%cnt = 0;
fprintf('Building Weight Matrix...\r')
% for kp = 1:n_proj
%     im_rot = imrotate(im_pad,-angles(kp),'bilinear','crop');
%     projmat(kp,:) = sum(im_rot,1);
%     t = -angles(kp)/180*pi;
%     R = [cos(t) -sin(t);sin(t) cos(t)];
%     %L = m*abs(sin(t))+n*abs(cos(t));
%     %offset = ceil((D-L)/2);
%     fprintf('\nAngle No.%d(%d Degree)\r',kp,angles(kp))
%     for kn = 1:N
%         [x,y] = ind2sub(sz,kn);
%         xy_rot = R*([x;y]-rc)+D/2+0.5;
%         idx = round(xy_rot(2));%#
%         %corresponding indice in W and p matrix
%         ixM = D*(kp-1)+idx;
%         W(ixM,kn) = 1;
%         %W(D*(kp-1)+idx,kn) = 1;
%         %W(cnt+idx-offset,kn) = 1;
%         %p(cnt+idx-offset) = projmat(kp,idx);%#
%         if ~ix(ixM)
%             p(ixM) = projmat(kp,idx);
%             ix(ixM) = true;
%         end
%         %p(D*(kp-1)+idx) = projmat(kp,idx);
%         %proj_info(cnt+idx-offset,:) = [kp,idx];
%     end
%     %cnt = cnt+(D-2*offset);
%     %fprintf('%d equations built.\r',D-2*offset)
% end
% %
[x,y] = ind2sub(sz,1:N);
%xy = [x;y]; %x,y coordinate
%x,y coordinate with respect to rotation center
xy_c = [x;y]-repmat(rc,1,N);
for kp = 1:n_proj
    try
        W_c{kp} = zeros(D,N); %weighting factor matrix
    catch expr
        fprintf(['Build W{%d}...\n' expr.message '\nGenerate a sparse.\r'],kp)
        W_c{kp} = sparse(D,N);
    end
    im_rot = imrotate(im_pad,-angles(kp),'bilinear','crop');
    pvec = sum(im_rot,1);
    t = -angles(kp)/180*pi;
    R = [cos(t) -sin(t);sin(t) cos(t)];
    fprintf('\nAngle No.%d(%d Degree)\r',kp,angles(kp))
    xy_rot = R*xy_c+D/2+0.5;
    idx = round(xy_rot(2,:));
    ixM = D*(kp-1)+idx; %corresponding indice in W and p matrix
    for kn = 1:N
        W_c{kp}(idx(kn),kn) = 1;
        p(ixM(kn)) = pvec(idx(kn));
        ix(ixM(kn)) = true;
    end
end
% %
%Convert cell to matrix
% try
%     W = zeros(M,N);
% catch expr
%     fprintf([expr.message '\nGenerate a sparse.\r'])
%     W = sparse(M,N);
% end
W = sparse(M,N);
for kp = 1:n_proj
    W((D*(kp-1)+(1:D)),:) = W_c{kp};
end

%Delete all-zero rows in W
fprintf('\nDelete all-zero rows in W...\r')
%ix = sum(W,2)==0;
W(~ix,:) = [];
p(~ix) = [];
n_eq = size(W,1);
fprintf('\nFinish building A. \n%d equations in total.\r',n_eq);
figure('name','Sparsity pattern of matrix W');
spy(W)
t_e = toc;
t = t_e-t_s;
fprintf('\nTotal: %f seconds\n\r',t);
end

author: Sable

repo: mcbench-benchmarks


text:
% Copyright 2017 Lime Microsystems Ltd.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%    http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%
% provide EVM for random sequences
% W-CDMA does not enforce equal amplitude signals, so we need AGC for each signal
% we assume we are dealing with BPSK or QAM, so rms amplitude should be unity
function evm=WCDMAULevm(x,qam)
  agc=1/sqrt(sum(x.*conj(x))/length(x));
  y=x*agc;
	evmref2=sum(abs(qam))/length(qam);
	[evm,loc]=min(abs(repmat(conj(qam'),1,length(y))-repmat(y,length(qam),1))); % decode QAM symbols
%	evm=evm./abs(qam(loc)); % some EVMs define relative to desired symbol
	evm/=evmref2; % some EVMs define relative to rms of QAM symbol table
end

author: myriadrf


text:
function [lin_restr,nonlin_restr,markov_chains]=create_restrictions_and_markov_chains0(markov_chains)
% create_restrictions_and_markov_chains0 -- creates restrictions for
% the constant-parameter SVAR model
%
% ::
%
%
%   [lin_restr,nonlin_restr,tpl]=create_restrictions_and_markov_chains0(tpl)
%
% Args:
%
%    - **markov_chains** [empty|struct]: structure of previously defined
%    markov chains
%
% Returns:
%    :
%
%    - **lin_restr** [cell]: cell array of restrictions (see below).
%
%    - **nonlin_restr** [cell]: cell array of inequality restrictions
%
%    - **markov_chains** [struct]: modified markov chains
%
% Note:
%
%    - The syntax to construct a restriction
%      --> ai(eqtn)
%      --> ai(eqtn,vbl)
%      --> ai(eqtn,vbl,chain_name,state)
%      --> a(eqtn)
%      --> a(eqtn,vbl)
%      --> a(eqtn,vbl,chain_name,state)
%      - **eqtn** [integer]: integer
%      - **vbl** [integer|char]: integer or variable name
%      - **i** [integer]: lag
%      - **chain_name** [char]: name of the markov chain
%      - **state** [integer]: state number
%
%    - The lag coefficients are labelled a0, a1, a2,...,ak, for a model with k
%    lags. Obviously, a0 denotes the contemporaneous coefficients.
%
%    - The constant terms labelled c_1_1, c_2_2,...,c_n_n, for a model with n
%    endogenous variables.
%
%    - The standard deviations labelled s_1_1, s_2_2,...,s_n_n, for a
%    model with n endogenous variables.
%
% Example:
%
%    See also:

if nargin==0||isempty(markov_chains)

    markov_chains=struct('name',{},...
    'states_expected_duration',{},...
    'controlled_parameters',{});

end

% syntax is coef(eqtn,vname,lag)
%-------------------------------
lin_restr={
    % first equation or "FFR" equation
    %----------------------------------
    'a1(1,pi)=0'
    'a2(1,pi)=0'
    'a1(1,ygap)=0'
    'a2(1,ygap)=0'
    'a2(1,FFR)=0'
    % second equation or "pi" equation
    %----------------------------------
    'a0(2,FFR)=0'
    'a1(2,FFR)=0'
    'a2(2,FFR)=0'
    'a1(2,ygap)=0'
    'a2(2,ygap)=0'
    % third equation or "ygap" equation
    %-----------------------------------
    'a1(3,FFR)=0'
    'a2(3,FFR)=0'
    'a1(3,pi)=0'
    'a2(3,pi)=0'
    'a0(3,pi)+a0(3,FFR)=0'
    };
nonlin_restr={
    'a0(3,FFR)>=0'
    'a1(1,FFR)>=0'
    'a1(1,FFR)<=1'
    };

end

author: jmaih


text:
%isiscalar - Test if parameter is a scalar (integer) satisfying an optional list of tests.
%
%  USAGE
%
%    test = isiscalar(x,test1,test2,...)
%
%    x              parameter to test
%    test1...       optional list of additional tests
%
%  EXAMPLES
%
%    % Test if x is a scalar (double)
%    isiscalar(x)
%
%    % Test if x is a strictly positive scalar (double)
%    isiscalar(x,'>0')
%
%    % Test if x is a scalar (double) included in [2,3]
%    isiscalar(x,'>=2','<=3')
%
%  NOTE
%
%    The tests ignore NaN, e.g. isiscalar(nan), isiscalar(nan,'>0') and isiscalar(nan,'<=0')
%    all return 1.
%
%  SEE ALSO
%
%    See also isdmatrix, isdvector, isdscalar, isimatrix, isivector, isstring,
%    islscalar, islvector, islmatrix.
%

% Copyright (C) 2010 by Michaël Zugaro
%
% 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.

function test = isiscalar(x,varargin)

% Check number of parameters
if nargin < 1,
  error('Incorrect number of parameters (type ''help isiscalar'' for details).');
end

% Test: double, scalar
test = isa(x,'double') & isscalar(x);
if ~test, return; end

% Test: integers?
test = test & round(x)==x;

% Optional tests
for i = 1:length(varargin),
	try
		if ~eval(['x' varargin{i} ';']), test = false; return; end
	catch err
		error(['Incorrect test ''' varargin{i} ''' (type ''help isiscalar'' for details).']);
	end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% main box sizes:
a=2;
b=1;

Nb=5; % number of boxes

% random boxes sizes:
% mab=mean([a b]);
% aa=0.05*mab+0.3*mab*rand(1,Nb);
% bb=0.05*mab+0.3*mab*rand(1,Nb);

mrg=0.95; % for margine
aa=[1    1     1     0.5  0.5];
bb=[1    0.25  0.25  0.5  0.5];
aa=aa*mrg;
bb=bb*mrg;

m2=min([aa bb]/2); % smallest half-size
AA=aa.*bb; % boxes areas

penalty=0.2*a*b;
nac=0.8; % negative area coefficient

N=5000; % population size
ng=20000; % number of generations
pmpe=0.1; % places exchange mutation probability
pmbj=0.02; % big gauss jump
pmsj=0.04; % small gauss jump
pmrr=0.1; % random rotation
pmvi=0.1; % random visible/invisible
pmne=0.2; % move to nearest adge

figure;
%ha1=axes;
ha1=subplot(2,1,1);
plot([0 a a 0 0], [0 0 b b 0],'b-');
xlim([-0.1*a 1.1*a]);
ylim([-0.1*b 1.1*b]);
set(ha1,'NextPlot','add');
ht=title(ha1,'start');
ha2=subplot(2,1,2);
drawnow;

set_cl; % set color table cl to plot boxes with different colors

% random initial population:
G=zeros(N,4*Nb);
Gch=zeros(N,4*Nb); % children
for Nc=1:N % for each individual
    G1=zeros(4,Nb); % one individual
    % G1(1,i)=1 if i-box is visible
    % G1(2,i)=1 if i-box is rotated at 90 degrees
    % G1(3,i) - x-coordinate of i-box center
    % G1(4,i) - y-coordinate of i-box center

    G1(1,:)=double(rand(1,Nb)<0.2);
    G1(2,:)=double(rand(1,Nb)<0.5);

    G1(3,:)=m2+(a-m2)*rand(1,Nb);
    G1(4,:)=m2+(b-m2)*rand(1,Nb);

    G(Nc,:)=(G1(:))'; % (G1(:))' converts matrix to row-vector
end

hi=imagesc(G,'parent',ha2);
drawnow;

Gpr1=zeros(4,Nb);
Gpr2=zeros(4,Nb); % two parents
Gch1=zeros(4,Nb);
Gch2=zeros(4,Nb); % two children
for ngc=1:ng % generations counting
    % find fitnesses:
    fitnesses=zeros(N,1);
    for Nc=1:N % for each individual
        G1(:)=(G(Nc,:))';
        vis=G1(1,:);
        ind=find(vis);
        L=length(ind);
        if L>0
            % only visible:
            rot=G1(2,ind);
            x=G1(3,ind);
            y=G1(4,ind);
            if L==1
                aaa=aa(ind);
                bbb=bb(ind);
                if rot
                    tmp=aaa;
                    aaa=bbb;
                    bbb=tmp;
                end
                A0=AA(ind); % box area
                x1=max([x-aaa/2  0]);
                y1=max([y-bbb/2  0]);
                x2=min([x+aaa/2  a]);
                y2=min([y+bbb/2  b]);
                % x1 - x2,  y1 - y2 is box (part of current box) that inside main box
                if (x1>=x2)||(y1>=y2)
                    A=0; % box that inside main box area
                else
                    A=(x2-x1)*(y2-y1); % box that inside main box area
                end
                %if A<A0 % if not fully inside main box
                if (aaa/2<=x)&&(x<=a-aaa/2)&&(bbb/2<=y)&&(y<=b-bbb/2) % if filly inside
                    fitness=A;
                else
                    fitness=A-nac*(A0-A)-penalty;
                end

            else
                fitness=0;
                ispen=false; % true if penality

                % check cross with main box:
                % add boxes arreas and strong subtract out areas:
                for n=1:L % for each box
                    ind1=ind(n);
                    aaa=aa(ind1);
                    bbb=bb(ind1);
                    if rot(n)
                        tmp=aaa;
                        aaa=bbb;
                        bbb=tmp;
                    end
                    A0=AA(ind1); % box area
                    x1=max([x(n)-aaa/2  0]);
                    y1=max([y(n)-bbb/2  0]);
                    x2=min([x(n)+aaa/2  a]);
                    y2=min([y(n)+bbb/2  b]);
                    % x1 - x2,  y1 - y2 is box (part of current box) that inside main box
                    if (x1>=x2)||(y1>=y2)
                        A=0; % box that inside main box area
                    else
                        A=(x2-x1)*(y2-y1); % box that inside main box area
                    end
                    %if A<A0 % if not fully inside main box
                        %fitness=fitness + A-nac*(A0-A);
                        %ispen=true; % penality
                    %else
                        %fitness=fitness + A;
                    %end

                    if (aaa/2<=x(n))&&(x(n)<=a-aaa/2)&&(bbb/2<=y(n))&&(y(n)<=b-bbb/2) % if filly inside
                        fitness=fitness + A;
                    else
                        fitness=fitness + A-nac*(A0-A);
                        ispen=true; % penality
                    end

                end

                % for each pair of boxes:
                for n1=1:L-1
                    ind1=ind(n1);
                    aaa1=aa(ind1);
                    bbb1=bb(ind1);
                    if rot(n1)
                        tmp=aaa1;
                        aaa1=bbb1;
                        bbb1=tmp;
                    end
                    A1=AA(ind1);
                    x1=x(n1);
                    y1=y(n1); % position of 1st box of pair
                    for n2=n1+1:L
                        ind2=ind(n2);
                        aaa2=aa(ind2);
                        bbb2=bb(ind2);
                        if rot(n2)
                            tmp=aaa2;
                            aaa2=bbb2;
                            bbb2=tmp;
                        end
                        A2=AA(ind2);
                        x2=x(n2);
                        y2=y(n2); % position of 2nd box of pair
                        dx=abs(x1-x2);
                        dy=abs(y1-y2); % distancies
                        a12=(aaa1/2+aaa2/2);
                        b12=(bbb1/2+bbb2/2);
                        if (dx<a12)&&(dy<b12) % if cross
                            ispen=true;
                            Ac=(a12-dx)*(b12-dy); % area of cross
                            fitness=fitness-Ac-Ac; % becuse area of n1 and n2 was added fully
                            fitness=fitness-2*nac*Ac;
                        end

                    end
                end

                if ispen
                    fitness=fitness-penalty;
                end

            end
        else
            fitness=0;
        end
        fitnesses(Nc)=fitness;
    end

    [fb bi]=max(fitnesses); % best

    % plot best:
    G1(:)=(G(bi,:))';
    Gb=G(bi,:); % best
    if mod(ngc,10)==0
        cla(ha1);
        Atmp=0;
        for Nbc=1:Nb
            vis1=G1(1,Nbc);
            if vis1
                rot1=G1(2,Nbc);
                aaa=aa(Nbc);
                bbb=bb(Nbc);
                if rot1
                    tmp=aaa;
                    aaa=bbb;
                    bbb=tmp;
                end
                x=G1(3,Nbc);
                y=G1(4,Nbc);
                plot([x-aaa/2  x+aaa/2  x+aaa/2  x-aaa/2  x-aaa/2],...
                     [y-bbb/2  y-bbb/2  y+bbb/2  y+bbb/2  y-bbb/2],...
                     '-','color',cl(Nbc,:),...
                     'parent',ha1);
                hold on;
                Atmp=Atmp+aaa*bbb;
            end
        end
        plot([0 a a 0 0], [0 0 b b 0],'b-','parent',ha1);
        xlim(ha1,[-0.1*a 1.1*a]);
        ylim(ha1,[-0.1*b 1.1*b]);

        set(hi,'Cdata',G);

        nvb=length(find(G1(1,:))); % number of visible boxes

        set(ht,'string',[' generation: ' num2str(ngc)  ', boxes: ' num2str(nvb) ', area: ' num2str(fb)]);

        drawnow;
    end

    % prepare for crossover, selection:
    fmn=min(fitnesses);
    fst=std(fitnesses);
    if fst<1e-7
        fst=1e-7;
    end
    fmn1=fmn-0.01*fst; % little low then minimum
    P=fitnesses-fmn1; % positive values
    p=P/sum(P); % probabilities
    ii=roulette_wheel_indexes(N,p);
    Gp=G(ii,:); % parents

    % crossover:
    for n=1:2:N
        pr1=Gp(n,:);
        pr2=Gp(n+1,:); % two parents
        % in matrix form:
        Gpr1(:)=pr1';
        Gpr2(:)=pr2';

        for Nbc=1:Nb

            % visibility:
            if rand<0.5
                Gch1(1,Nbc)=Gpr1(1,Nbc);
            else
                Gch1(1,Nbc)=Gpr2(1,Nbc);
            end
            if rand<0.5
                Gch2(1,Nbc)=Gpr1(1,Nbc);
            else
                Gch2(1,Nbc)=Gpr2(1,Nbc);
            end

            % rotation:
            if rand<0.5
                Gch1(2,Nbc)=Gpr1(2,Nbc);
            else
                Gch1(2,Nbc)=Gpr2(2,Nbc);
            end
            if rand<0.5
                Gch2(2,Nbc)=Gpr1(2,Nbc);
            else
                Gch2(2,Nbc)=Gpr2(2,Nbc);
            end

            % position:
            % child 1:
            %i3=ceil(3*rand);
            %i3=roulette_wheel_indexes(1,[0.2 0.4 0.4]);
            i3=1+ceil(2*rand);
            switch i3
                case 1 % get mean position
                    Gch1(3,Nbc)=(Gpr1(3,Nbc)+Gpr2(3,Nbc))/2;
                    Gch1(4,Nbc)=(Gpr1(4,Nbc)+Gpr2(4,Nbc))/2;
                case 2 %get position of parent 1
                    Gch1(3,Nbc)=Gpr1(3,Nbc);
                    Gch1(4,Nbc)=Gpr1(4,Nbc);
                case 3 %get position of parent 2
                    Gch1(3,Nbc)=Gpr2(3,Nbc);
                    Gch1(4,Nbc)=Gpr2(4,Nbc);
            end
            % child 2:
            %i3=ceil(3*rand);
            %i3=roulette_wheel_indexes(1,[0.2 0.4 0.4]);
            i3=1+ceil(2*rand);
            switch i3
                case 1 % get mean position
                    Gch2(3,Nbc)=(Gpr1(3,Nbc)+Gpr2(3,Nbc))/2;
                    Gch2(4,Nbc)=(Gpr1(4,Nbc)+Gpr2(4,Nbc))/2;
                case 2 %get position of parent 1
                    Gch2(3,Nbc)=Gpr1(3,Nbc);
                    Gch2(4,Nbc)=Gpr1(4,Nbc);
                case 3 %get position of parent 2
                    Gch2(3,Nbc)=Gpr2(3,Nbc);
                    Gch2(4,Nbc)=Gpr2(4,Nbc);
            end

        end
        ch1=(Gch1(:))';
        ch2=(Gch2(:))';
        Gch(n,:)=ch1;
        Gch(n+1,:)=ch2;

    end
    G=Gch; % now children

    % mutations:
    % places exchange
    for Nc=1:N % for each individual
        if rand<pmpe
            G1(:)=(G(Nc,:))';
            ir1=ceil(Nb*rand);
            ir2=ceil(Nb*rand);
            tmp1=G1(3:4,ir1);
            G1(3:4,ir1)=G1(3:4,ir2);
            G1(3:4,ir2)=tmp1;
            G(Nc,:)=(G1(:))';
        end
    end

    % big gauss jump:
    for Nc=1:N % for each individual
        if rand<pmbj
            G1(:)=(G(Nc,:))';
            ir=ceil(Nb*rand);
            G1(3:4,ir)=G1(3:4,ir)+[0.05*a*randn;
                                   0.05*b*randn];
            G(Nc,:)=(G1(:))';
        end
    end

    % small gauss jump:
    for Nc=1:N % for each individual
        if rand<pmsj
            G1(:)=(G(Nc,:))';
            ir=ceil(Nb*rand);
            G1(3:4,ir)=G1(3:4,ir)+[0.005*a*randn;
                                   0.005*b*randn];
            G(Nc,:)=(G1(:))';
        end
    end

    % random rotation:
    for Nc=1:N % for each individual
        if rand<pmrr
            G1(:)=(G(Nc,:))';
            ir=ceil(Nb*rand);
            G1(2,ir)=double(rand<0.5);
            G(Nc,:)=(G1(:))';
        end
    end

    % random visible/invisible:
    for Nc=1:N % for each individual
        if rand<pmvi
            G1(:)=(G(Nc,:))';
            ir=ceil(Nb*rand);
            G1(1,ir)=double(rand<0.5);
            G(Nc,:)=(G1(:))';
        end
    end

    % move to nearest edge:
    for Nc=1:N % for each individual
        if rand<pmne
            G1(:)=(G(Nc,:))';
            ir=ceil(Nb*rand); % random small box
            rv=find((G1(1,:))&((1:Nb)~=Nc)); % find rest visible
            if rand<0.5
                % to veritcile edge
                eax=[G1(3,rv)-aa(rv)/2  G1(3,rv)+aa(rv)/2  0  a]; % edge xs
                deax=[(G1(3,ir)-aa(ir)/2) - eax  (G1(3,ir)+aa(ir)/2) - eax]; % distancies
                [dmn indm]=min(abs(deax));
                G1(3,ir)=G1(3,ir)-deax(indm);
            else
                % to horizontal edge
                eay=[G1(4,rv)-bb(rv)/2  G1(4,rv)+bb(rv)/2  0  b]; % edge ys
                deay=[(G1(4,ir)-bb(ir)/2) - eay  (G1(4,ir)+bb(ir)/2) - eay]; % distancies
                [dmn indm]=min(abs(deay));
                G1(4,ir)=G1(4,ir)-deay(indm);
            end
        end
    end

    % ellitism:
    G(1,:)=Gb;

end


text:
function this = MrAffineTransformation_transformation(this)
% Unit test for MrAffineTransformation applying transformations
%
%   Y = MrUnitTest()
%   run(Y, 'MrAffineTransformation_transformation')
%
% This is a method of class MrUnitTest.
%
% IN
%
% OUT
%
% EXAMPLE
%   MrAffineTransformation_transformation
%
%   See also MrUnitTest

% Author:   Saskia Bollmann
% Created:  2018-01-16
% Copyright (C) 2018 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.
%

% construct MrAffineTransformation object from sampling points
affineTransformation = this.make_affineTransformation_reference(0);

% define expected solution
expSolution = affineTransformation.copyobj;

% create transformation matrix
transformationMatrix = [-0.2769 0.6948 0.4387 18.69; ...
    0.0462 -0.3171 0.3816 4.898; ...
    0.0971 0.9502 -0.7655 4.456; ...
    -0.8235 0.0344 0.7952 0.6463];
% apply transformation matrix
affineTransformation.apply_transformation(transformationMatrix);
% verify that affine transformation is applied
if affineTransformation.isequal(expSolution)
    this.assertFail('The transformation matrix has not been applied.');
end
% apply inverse transformation matrix
affineTransformation.apply_inverse_transformation(transformationMatrix);

% define actual solution
actSolution = affineTransformation;

% verify equality of expected and actual solution
% import matlab.unittests to apply tolerances for objects
import matlab.unittest.TestCase
import matlab.unittest.constraints.IsEqualTo
import matlab.unittest.constraints.AbsoluteTolerance
import matlab.unittest.constraints.PublicPropertyComparator

this.verifyThat(actSolution, IsEqualTo(expSolution,...
    'Within', AbsoluteTolerance(10e-7),...
    'Using', PublicPropertyComparator.supportingAllValues));

author: translationalneuromodeling

repo: tapas


text:
classdef InertialSideSlipAngleConstraint < AbstractConstraint
    %InertialSideSlipAngleConstraint Summary of this class goes here
    %   Detailed explanation goes here

    properties
        normFact = 1;
        event LaunchVehicleEvent
        eventNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        evalType(1,1) ConstraintEvalTypeEnum = ConstraintEvalTypeEnum.FixedBounds;
        stateCompType(1,1) ConstraintStateComparisonTypeEnum = ConstraintStateComparisonTypeEnum.Equals;
        stateCompEvent LaunchVehicleEvent
        stateCompNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;
    end

    methods
        function obj = InertialSideSlipAngleConstraint(event, lb, ub)
            obj.event = event;
            obj.lb = lb;
            obj.ub = ub;

             obj.id = rand();
        end

        function [lb, ub] = getBounds(obj)
            lb = obj.lb;
            ub = obj.ub;
        end

        function [c, ceq, value, lwrBnd, uprBnd, type, eventNum, valueStateComp] = evalConstraint(obj, stateLog, celBodyData)
            type = obj.getConstraintType();

            switch obj.eventNode
                case ConstraintStateComparisonNodeEnum.FinalState
                    stateLogEntry = stateLog.getLastStateLogForEvent(obj.event);

                case ConstraintStateComparisonNodeEnum.InitialState
                    stateLogEntry = stateLog.getFirstStateLogForEvent(obj.event);

                otherwise
                    error('Unknown event node.');
            end

            ut = stateLogEntry.time;
            rVect = stateLogEntry.position;
            vVect = stateLogEntry.velocity;
            bodyInfo = stateLogEntry.centralBody;

            [~,~,inertAngOfSideslip] = stateLogEntry.attitude.getInertialAeroAngles(ut, rVect, vVect, bodyInfo);
            value = rad2deg(inertAngOfSideslip);

            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                switch obj.stateCompNode
                    case ConstraintStateComparisonNodeEnum.FinalState
                        stateLogEntryStateComp = stateLog.getLastStateLogForEvent(obj.stateCompEvent).deepCopy();

                    case ConstraintStateComparisonNodeEnum.InitialState
                        stateLogEntryStateComp = stateLog.getFirstStateLogForEvent(obj.stateCompEvent).deepCopy();

                    otherwise
                        error('Unknown event node.');
                end

                cartElem = stateLogEntryStateComp.getCartesianElementSetRepresentation();
                cartElem = cartElem.convertToFrame(stateLogEntry.centralBody.getBodyCenteredInertialFrame());
                stateLogEntryStateComp.setCartesianElementSet(cartElem);

                ut = stateLogEntryStateComp.time;
                rVect = stateLogEntryStateComp.position;
                vVect = stateLogEntryStateComp.velocity;
                bodyInfo = stateLogEntryStateComp.centralBody;

                [~,~,inertAngOfSideslip] = stateLogEntryStateComp.attitude.getInertialAeroAngles(ut, rVect, vVect, bodyInfo);
                valueStateComp = rad2deg(inertAngOfSideslip);
            else
                valueStateComp = NaN;
            end

            [c, ceq] = obj.computeCAndCeqValues(value, valueStateComp);

            lwrBnd = obj.lb;
            uprBnd = obj.ub;

            eventNum = obj.event.getEventNum();
        end

        function sF = getScaleFactor(obj)
            sF = obj.normFact;
        end

        function setScaleFactor(obj, sF)
            obj.normFact = sF;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesEvent(obj, event)
            tf = obj.event == event;
            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                tf = tf || obj.stateCompEvent == event;
            end
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end

        function tf = usesExtremum(obj, extremum)
            tf = false;
        end

        function tf = canUseSparseOutput(obj)
            tf = true;
        end

        function event = getConstraintEvent(obj)
            event = obj.event;
        end

        function type = getConstraintType(obj)
            type = 'Inertial Side Slip Angle';
        end

%         function name = getName(obj)
%             name = sprintf('%s - Event %i', obj.getConstraintType(), obj.event.getEventNum());
%         end

        function [unit, lbLim, ubLim, usesLbUb, usesCelBody, usesRefSc] = getConstraintStaticDetails(obj)
            unit = 'deg';
            lbLim = -360;
            ubLim = 360;
            usesLbUb = true;
            usesCelBody = false;
            usesRefSc = false;
        end

        function addConstraintTf = openEditConstraintUI(obj, lvdData)
%             addConstraintTf = lvd_EditGenericMAConstraintGUI(obj, lvdData);

            output = AppDesignerGUIOutput({false});
            lvd_EditGenericMAConstraintGUI_App(obj, lvdData, output);
            addConstraintTf = output.output{1};
        end
    end

    methods(Static)
        function constraint = getDefaultConstraint(~, ~)
            constraint = InertialSideSlipAngleConstraint(LaunchVehicleEvent.empty(1,0),0,0);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function L = laplacian( F )
%LAPLACIAN Vector Laplacian of a DISKFUNV.
%   LAPLACIAN(F) returns a DISKFUNV representing the vector Laplacian of F.
%
% See also DISKFUN/LAPLACIAN

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Empty check:
if ( isempty( F ) )
    L = diskfunv();
    return
end

% The vector LAPLACIAN of a DISKFUNV is equal to F_xx + F_yy:
L = diff(F, 1, 2) + diff(F, 2, 2);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function a = i4_to_a ( i )

%*****************************************************************************80
%
%% I4_TO_A returns the I-th alphabetic character.
%
%  Example:
%
%    I   A
%
%   -8  ' '
%    0  ' '
%    1  'A'
%    2  'B'
%   ..
%   26  'Z'
%   27  'a'
%   52  'z'
%   53  ' '
%   99  ' '
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the index of the letter to be returned.
%    0 is a space;
%    1 through 26 requests 'A' through 'Z', (ASCII 65:90);
%    27 through 52 requests 'a' through 'z', (ASCII 97:122);
%
%    Output, character A, the requested alphabetic letter.
%
  if ( i <= 0 )
    a = ' ';
  elseif ( 1 <= i & i <= 26 )
    a = 'A' + i - 1;
  elseif ( 27 <= i & i <= 52 )
    a = 'a' + i - 27;
  else
    a = ' ';
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
a=VideoReader("train_vid.mp4");
for img=1:a.NumberOfFrames;
    filename=strcat('frame',num2str(img),'.jpg');
       b = read(a, img);
    imwrite(b,filename);
end

author: Aarchishya


text:
function jac = autoJac(fun,x0,varargin)
%AUTOJAC  Returns a Automatically Differentiated (AD) Jacobian
%
%   jac = autoJac(fun,x0) uses the user supplied Matlab function to
%   generate the Jacobian using automatic differentiation. The user
%   function must be a Matlab function and must not call external code or
%   class / toolbox functions. If your function breaks any of these
%   conditions consider using mklJac instead.
%
%   The underlying AD algorithm is adiff by William McIlhagga and its
%   documentation pdf is provided in the Differentiation folder. See the
%   BSD license below the code.

%   Copyright (C) 2011 Jonathan Currie (I2C2)

if(~isa(fun,'function_handle'))
    error('Fun should be a function handle!');
end

[~,jac] = adiffget(fun(adiff(x0),varargin{:}));

end

% Copyright (c) 2010, William McIlhagga
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
classdef spatiotemporal_datasource
   %UNTITLED2 Summary of this class goes here
   %   Detailed explanation goes here

   properties
   end

   methods
      function data = retrieve(obj, where_, from_ , until_)
         error('unimplemented base datasource class');
      end

   end

end

author: geoscience-community-codes

repo: GISMO


text:
%DRAWMATCHES  Draws the found matches of keypoints from two images
%
%     out = cv.drawMatches(im1, keypoints1, im2, keypoints2, matches1to2)
%     out = cv.drawMatches(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __im1__ First source image.
% * __keypoints1__ Keypoints from the first source image. A 1-by-N structure
%   array with the following fields:
%   * __pt__ coordinates of the keypoint `[x,y]`
%   * __size__ diameter of the meaningful keypoint neighborhood
%   * __angle__ computed orientation of the keypoint (-1 if not applicable).
%     Its possible values are in a range [0,360) degrees. It is measured
%     relative to image coordinate system (y-axis is directed downward), i.e
%     in clockwise.
%   * __response__ the response by which the most strong keypoints have been
%     selected. Can be used for further sorting or subsampling.
%   * __octave__ octave (pyramid layer) from which the keypoint has been
%     extracted.
%   * **class_id** object id that can be used to clustered keypoints by an
%     object they belong to.
% * __im2__ Second source image.
% * __keypoints2__ Keypoints from the second source image. Same format as
%   `keypoints1`.
% * __matches1to2__ Matches from the first image to the second one, which means
%   that `keypoints1(i)` has a corresponding point in `keypoints2(matches(i))`.
%   A 1-by-M structure array with the following fields:
%   * __queryIdx__ query descriptor index (zero-based index)
%   * __trainIdx__ train descriptor index (zero-based index)
%   * __imgIdx__ train image index (zero-based index)
%   * __distance__ distance between descriptors (scalar)
%
% ## Output
% * __out__ Output image. Its content depends on the option values defining
%   what is drawn in the output image. See possible options below. By default,
%   the two source images, matches, and single keypoints will be drawn. For
%   each keypoint, only the center point will be drawn (without a circle
%   around the keypoint with the keypoint size and orientation).
%
% ## Options
% * __MatchColor__ Color of matches (lines and connected keypoints). If all -1,
%   the color is generated randomly. default [-1,-1,-1,-1].
% * __SinglePointColor__ Color of single keypoints (circles), which means that
%   keypoints do not have the matches. If all -1, the color is generated
%   randomly. default [-1,-1,-1,-1].
% * __MatchesMask__ Mask determining which matches are drawn. If the mask is
%   empty, all matches are drawn. default empty.
% * __NotDrawSinglePoints__ Single keypoints will not be drawn. default false
% * __DrawRichKeypoints__ For each keypoint, the circle around keypoint with
%   keypoint size and orientation will be drawn. default false.
% * __OutImage__ If set, matches will be drawn on existing content of output
%   image, otherwise source image is used instead. Default not set
%
% This function draws matches of keypoints from two images in the output
% image. Match is a line connecting two keypoints (circles).
%
% In pseudo-code, the function works as follows:
%
%     for m=1:numel(matches1to2)
%         if ~MatchesMask(m), continue, end
%         kp1 = keypoints1(matches1to2(m).queryIdx + 1);
%         kp2 = keypoints2(matches1to2(m).trainIdx + 1);
%         draw_match(kp1, kp2);
%     end
%
% See also: cv.drawKeypoints, showMatchedFeatures
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [B,L] = ordered_outline(F)
  % ORDERED_OUTLINE Find outline (boundary) edges of mesh
  %
  % [B,L] = outline(F)
  %
  % Input:
  %   F  #F by 3 face list of indices
  % Outputs:
  %   B  #B by 1 list of outline edges
  %   L  #loops+1 by 1 list of boundary loop start indices into B, the last
  %     entries is (by tradition) always the numel of B + 1
  %
  % Example:
  %   V = [0 0; 1 0; 1 1 ; 0 1; 4 0; 4 4; 0 4];
  %   F = [1 2 3; 1 3 4; 5 6 7];
  %   [B,L] = ordered_outline(F);
  %   hold on
  %   for l = 1:(numel(L)-1)
  %     plot(V([B(L(l):(L(l+1)-1)) B(L(l))],1),V([B(L(l):(L(l+1)-1)) B(L(l))],2));
  %   end
  %   hold off
  %
  % See also: outline
  %

  if size(F,2) == 2
    O = F;
  else
    O = outline(F);
  end
  % determine uniqueness of indices
  [u,m,n] = unique(O(:),'rows');
  % determine counts for each unique edge
  counts = accumarray(n(:), 1);
  % all counts should be 2
  assert(all(counts == 2));

  % number of outline edges
  no = size(O,1);

  % list of boundary
  B = [];
  L = [];
  unseen = true(no,1);
  while any(unseen)
    L = [L numel(B)+1];
    % find first unseen index
    f = find(unseen,1);
    % append start vertex to boundary
    next = O(f,1);
    % keep track of start
    start = next;
    unseen(f) = false;
    B = [B next];
    next = O(f,2);
    % loop until back at start
    while next ~= start
      B = [B next];
      % try to find next as source
      f = find(unseen & O(:,1) == next);
      if isempty(f)
        % try to find next as dest
        f = find(unseen & O(:,2) == next);
        assert(~isempty(f));
        next = O(f,1);
      else
        next = O(f,2);
      end
      unseen(f) = false;
    end
  end
  L = [L numel(B)+1];

end

author: alecjacobson

repo: gptoolbox


text:
function [ a, sX ] = size_compat( sX, sY )
a = true;
switch class( sX ),
    case 'double',
        if isempty( sX ) || all( sX == 0 ),
            sX = sY;
        elseif isempty( sY ) || all( sY == 0 ),
        elseif ~isequal( sX, sY ),

            % Feb 29, 2012. Special case:
            %   One represents the size a x b x c, where c = 1
            %   The other is a x b (since Matlab often automatically squeezes
            %   3D arrays to 2D if the 3rd dimension is a singletone)
%             if (length(sX) >= 3 && length(sX) == length(sY)+1 && sX(end)==1) || ...
%                     (length(sY) >= 3 && length(sY) == length(sX)+1 && sY(end)==1)
                % do nothing
            % March 2012, a better fix (also due to Graham Coleman)
            if min( length(sX),length(sY) ) >= 2
                truncA  = cutTrailingOnes(sX);
                truncB  = cutTrailingOnes(sY);
                a = isequal( truncA, truncB );
            else
                a = false;
            end
        end
    case 'cell',
        if ~isa( sY, 'cell' ) || numel( sX ) ~= numel( sY ) || isa( sX{1}, 'function_handle' ) && ~isequal( sX, sY ),
            a = false;
        elseif isa( sX{1}, 'function_handle' ),
            a = isequal( sX, sY );
        else
            for k = 1 : numel( sX ),
                [ ta, sX{k} ] = size_compat( sX{k}, sY{k} );
                a = a && ta;
            end
        end
    otherwise,
        a = isequal( sX, sY );
end
if ~a,
    sX = [];
end

function y = cutTrailingOnes( x )
%cuts the final ones of a vector, and columnize

lastNotOne = find( x(:)~=1, 1, 'last' );

%do not cut before 2nd position
lastNotOne = max( 2, lastNotOne );

y = x(1:lastNotOne);

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function [er,output2,output3] = errorFunction( errFcn, varargin )
% err = errorFunction( errFcn, arguments )
%   just returns the evaluation of errFcn( arguments ).
%   This seems pointless, but the point is that this function
%   has memory of all the errors, so you look at this later.
%
% errHist = errorFunction()
%   will return the history of errors and also clear the memory
%
% ... = errorFunction( ticReference )
%   will allow the function to enable time logging
%   It will subtract off the time used for the error function
%
% [errHist, timeLog] = errorFunction( )
%   will output the time logging information.
%   The time log shows at what time, relative to ticReference,
%   the error was recorded
%
% The function uses persistent variables so don't forget to clear it
%   and set it with a new ticReference!
%
% Stephen Becker, March 6 2014 stephen.beckr@gmail.com
% March 11 2014, updated to subtract off the time involved
%   for computing the error. Also allows error function
%   to return a row vector instead of a scalar
% March 14 2014, if ticReference is really a cell array with
%   {ticReference,maxTime}, then this will throw an error
%   after reaching "maxTime" seconds

persistent erHist timeLog ticReference timeCorrection maxTime
if nargin <= 1
    if nargin == 1 && ~isempty(errFcn)
        if iscell( errFcn )
            ticReference = errFcn{1};
            maxTime      = errFcn{2};
        else
            maxTime     = Inf;
            ticReference = errFcn;
        end
    else
        ticReference = [];
    end
    er = erHist;
    if nargout >= 2, output2 = timeLog; end
    if nargout >= 3, output3 = timeCorrection; end
    erHist = [];
    timeLog= [];
    timeCorrection = [];
    return
end
if isempty( maxTime ), maxTime = Inf; end

tc1 = tic;
er = errFcn( varargin{:} );
erTime = toc(tc1);
erHist(end+1,:) = er;

% And mark at what time we took this error measurement
if ~isempty( ticReference )
    % And we also want to add "erTime" to ticReference
    % but this isn't straightforward since ticReference is some huge
    % uint64 and not in seconds. So instead we update a correction
    if isempty( timeCorrection), timeCorrection = 0; end
    timeCorrection = timeCorrection + erTime;
    timeLog(end+1) = toc( ticReference ) - timeCorrection;

    if timeLog(end) > maxTime
        error('errorFunction:timedOut','Reached maximum allowed time');
    end
end

author: andrewssobral

repo: lrslibrary


text:
function valid=is_valid_state(x,x_dot,theta,theta_dot)

one_degree=0.0174532;	% 2pi/360 */
six_degrees=0.1047192;
twelve_degrees=0.2094384;
fifty_degrees=0.87266;

valid=1;
if (x < -2.4 | x > 2.4  | theta < -twelve_degrees | theta > twelve_degrees)
    valid=-1;
end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function [I J] = itril(sz, k)
% function [I J] = itril(sz, k) % OR
% I = itril(sz, k)
%
% Return the subindices [I J] (or linear indices I if single output call)
% in the purpose of extracting an lower triangular part of the matrix of
% the size SZ. Input k is optional shifting. For k=0, extract from the main
% diagonal. For k>0 -> above the diagonal, k<0 -> below the diagonal
%
% This returnd same as [...] = find(tril(ones(sz),k))
% - Output is a column and sorted with respect to linear indice
% - No intermediate matrix is generated, that could be useful for large
%   size problem
% - Mathematically, A(itril(size(A)) is called (lower) "half-vectorization"
%   of A
%
% Example:
%
% A = [ 7     5     4
%       4     2     3
%       9     1     9
%       3     5     7 ]
%
% I = itril(size(A))  % gives [1 2 3 4 6 7 8 11 12]'
% A(I)                % gives [7 4 9 3 2 1 5  9  7]' OR A(tril(A)>0)
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% Date: 21/March/2009

if isscalar(sz)
    sz = [sz sz];
end
m=sz(1);
n=sz(2);

% Main diagonal by default
if nargin<2
    k=0;
end

nc = min(n,m+k); % number of columns of the triangular part
lo = max((1:nc).'-k,1); % lower row indice for each column
hi = m + zeros(nc,1); % upper row indice for each column

if isempty(lo)
    I = zeros(0,1);
    J = zeros(0,1);
else
    c=cumsum([0; hi-lo]+1); % cumsum of the length
    I = accumarray(c(1:end-1), (lo-[0; hi(1:end-1)]-1), ...
                   [c(end)-1 1]);
    I = cumsum(I+1); % row indice
    J = cumsum(accumarray(c,1));
    J = J(1:end-1); % column indice
end

if nargout<2
    % convert to linear indices
    I = sub2ind([m n], I, J);
end

end % itril


text:
function [m2] = km22m2(km2)
% Convert area from square kilometers to square meters.
% Chad A. Greene 2012
m2 = km2*1000000;


text:
function data_num = p06_data_num ( )

%*****************************************************************************80
%
%% P06_DATA_NUM returns the dimension of the data vector for problem 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer DATA_NUM, the dimension of the data vector.
%
  num_int = 5;

  data_num = 1 + num_int * ( num_int + 1 );

  return
end


text:
function centroids = kMeansInitCentroids(X, K)
%KMEANSINITCENTROIDS This function initializes K centroids that are to be
%used in K-Means on the dataset X
%   centroids = KMEANSINITCENTROIDS(X, K) returns K initial centroids to be
%   used with the K-Means on the dataset X
%

% You should return this values correctly
centroids = zeros(K, size(X, 2));

% ====================== YOUR CODE HERE ======================
% Instructions: You should set centroids to randomly chosen examples from
%               the dataset X
%

% Initialize the centroids to be random examples

% Randomly reorder the indices of examples
randidx = randperm(size(X, 1));
% Take the first K examples as centroids
centroids = X(randidx(1:K), :);

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

end

repo: machine-learning


text:
function OutMatrix = IncompleteBellPoly(Nin,Kin,DataList)
%INCOMPLETEBELLPOLY  An iterative method for computing the incomplete (also known as the second kind of) Bell polynomials.
%
%   Purpose:
%   Given a list of input values (x_{1},x_{2},...,x_{N}), the script returns a matrix of Bell polynomials B_{n,k} for
%   n=0,...,N and k=0,...,K.
%
%   Syntax:
%   OutMatrix = IncompleteBellPoly(Nin,Kin,DataList)
%   where
%   B_{n,k} = OutMatrix(n+1,k+1)
%   n=0,...,Nin
%   k=0,...,Kin (Kin<=Nin)
%   DataList = (x_{1},x_{2},...,x_{Nin})
%
%   Authors:
%   Patrick Kano and Moysey Brio
%   University of Arizona
%
%   Email:
%   brio@math.arizona.edu
%
%   Latest Modification Date:
%   April 3, 2007
%
%   Discussion:
%   Given Taylor expansion coefficients of a function g(t) {g_{0},g_{1},g_{2},...} with g_{0}=0,
%   B_{n,k}(g_{0},g_{1},...,g_{n-k+1}) is the nth Taylor coefficient of the kth derivative of g(t)/(k!) in terms of {g_{0},g_{1},g_{2},...}
%   \frac{1}{k!} g^{k}(t) = \sum_{n=0}^{\infty} B_{n,k} \frac{t^{n}}{n!}
%
%   The Bell polynomials can be computed efficiently by a recursion relation
%   B_{n,k} = \sum_{m=1}^{n-k+1} \binom{n-1}{m-1} g_{m} B_{n-m,k-1}
%   where
%   B_{0,0} = 1;
%   B_{n,0} = 0; for n=>1
%   B_{0,k} = 0; for k=>1
%
%   The coefficients can be stored in a lower triangular matrix.  The elements of the kth column are the Taylor coefficients
%   of the kth derivative of g(t)/k!.
%
%   In application, the polynomials arise in multiple contexts in combinatorics.  They also can be found in Riordan's
%   formulation of di Bruno's formula for computing an arbitrary order derivative of the composition of two functions
%   \frac{d^{n}}{dt^{n}} g(f(t)) = \sum_{k=0}^{n} g^{k}(f(t)) B_{n,k}(f^{1}(t),f^{2}(t),...,f^{n-k+1}(t))
%
%   Script Check:
%   If DataList=1 for all entries, B_{n,k} = S(n,k) = Stirling number of the second kind for (n,k)
%
%   Failure Return:
%   OutMatrix is undefined if the code fails.  An error statement is issued and the function exits.
%
%   References:
%   Ferrell S. Wheeler, Bell Polynomials, ACM SIGSAM Bulletin, vol. 21, issue 3, pp.44-53, 1987.
%
%   Warren P. Johnson, The curious history of Faa di Bruno's formula, The American mathematical monthly, vol. 109, no. 3, pp. 217-234, March 2002.
%
%   http://en.wikipedia.org/wiki/Bell_polynomials
%
%   Supported by: Grant # NSF ITR-0325097

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

%n - row index
%k - column index

%Check that sufficient input variables x_{1},x_{2},...,x_{N} have been given.
%Note: To compute B_{n,k} requires (x_{1},x_{2},...,x_{n-k+1}) input terms.
%However, we want B_{n=1...N,k=1...N), which includes B_{n=N,k=1}.
%To determine B_{N,1}, we thus require (x_{1},x_{2},...,x_{N}) input terms.

ListLength = length(DataList);

if(ListLength<Nin)
 ErrorStatement = sprintf('Insufficient number of user defined input values %d.  Computation requires %d.  Exiting...\n',ListLength,Nin);
 error(ErrorStatement);

elseif(Nin<0 || Kin<0)
 ErrorStatement = sprintf('N and K must be non-negative.  Exiting... \n');
 error(ErrorStatement);

elseif(Nin<Kin)
 ErrorStatement = sprintf('K must be less than or equal to N.  Exiting... \n');
 error(ErrorStatement);

elseif(Nin==0 && Kin==0)
 OutMatrix = [1];
 return;

elseif(Nin>0 && Kin==0)
 OutMatrix = zeros(Nin+1,1);
 OutMatrix(1,1) = 1;
 return;

else
 %Compute the recursion relationship
 %Note: To keep with Matlab's convention for starting indices at 1 we compute only the non-trivial elements of the matrix.

 %Initialize the Triangular Matrix
 Bm = zeros(Nin,Kin);

 %Starting values for the recurrence relations
 %Bm(0,0) = 1;
 %Bm(n,0) = 0; for n=>1
 %Bm(0,k) = 0; for k=>1
 %kidx=1
 Bm(1:Nin,1) = DataList(1:Nin);

 %Generate each required row of data
 for(nidx=2:Nin) %row, Note: Row 1 has only 1 non-zero element

  for(kidx=2:Kin) %column, Note: Column 1 is the given input array

   for(midx=1:nidx-kidx+1) %recursion
    Bm(nidx,kidx) = Bm(nidx,kidx) + nchoosek(nidx-1,midx-1)*DataList(midx)*Bm(nidx-midx,kidx-1);
   end

  end
 end

 %Append [1,0,...] to generate the complete output matrix
 OutMatrix = eye(Nin+1,Kin+1);
 OutMatrix(2:Nin+1,2:Kin+1) = Bm;
end

end %End of the Function


text:
% LOTKA-VOLTERRA system
% System identification: DelayDMDc

clear all, close all, clc
figpath = '../FIGURES/';
datapath = '../DATA/';
addpath('../utils');

ModelName = 'SparseModels';
%% Generate Data
InputSignalType = 'sphs';%prbs; chirp; noise; sine2; sphs; mixed
Ndelay = 1;%35;%35;
ONLY_TRAINING_LENGTH = 1;
getTrainingData

DataTrain.x = x;
DataTrain.t = t;
DataTrain.tspan = tspan;
DataTrain.u = u;
DataTrain.xmean = xmean;
xstd = std(DataTrain.x(:,1));

%% Parameters
close all
rng(0,'twister')

eta_vec = [0.01 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5];
N_ETA = length(eta_vec);

lambda_vec = [0,1e-5,1e-4,1e-3,1e-2,1e-1,1e-0,1e+1];
N_LAMBDA = length(lambda_vec);

model_collection = {'sDMDc','sEDMDc2', 'sEDMDc3','SINDYc'};
N_MODEL = length(model_collection);

N_REP = 30;

options_method.sparsify = 'ILSTH'; % 'LASSO'
options_method.usesine  = 0;

nVars = 2;
%% Initialization
TrainingError = inf*ones(N_ETA,N_LAMBDA,N_MODEL,N_REP);
ValidationError = inf*ones(N_ETA,N_LAMBDA,N_MODEL,N_REP);

PredictionResults_Training = zeros(1000,nVars,N_ETA,N_LAMBDA,N_MODEL,N_REP);
PredictionResults_Validation = zeros(1000,nVars,N_ETA,N_LAMBDA,N_MODEL,N_REP);

%% Model Identification

for iR = 1:N_REP
    counter = 0;

    tstart = tic;
    for iEta = 1:N_ETA

        % Add noise to data
        eps     = eta_vec(iEta)*xstd;
        x       = DataTrain.x + eps*randn(size(DataTrain.x));
        u       = DataTrain.u;
        M       = length(x);

        for iLambda = 1:N_LAMBDA
            lambda = lambda_vec(iLambda);

            %% Model Identification
            Hx  = getHankelMatrix_MV(x - repmat(DataTrain.xmean,[M 1]),1); % No time delay (=1)
            Hu  = getHankelMatrix_MV(u',1);

            % sparse DMDc
            options_method.order  = 1;
            yout = poolData(Hx',nVars,options_method.order,0);
            Y = yout(1:end-1,2:end)'; % remove constant
            Yp = yout(2:end,2:end)';
            U = Hu(1:end-1);
            sys_DMDc = SparseRegression(Y,Yp,U,dt, options_method, lambda);

            % sparse EDMDc
            options_method.order  = 2;
            yout = poolData(Hx',nVars,options_method.order,0);
            Y = yout(1:end-1,2:end)'; % remove constant
            Yp = yout(2:end,2:end)';
            U = Hu(1:end-1);
            sys_EDMDc2 = SparseRegression(Y,Yp,U,dt, options_method, lambda);

            % sparse EDMDc
            options_method.order  = 3;
            yout = poolData(Hx',nVars,options_method.order,0);
            Y = yout(1:end-1,2:end)'; % remove constant
            Yp = yout(2:end,2:end)';
            U = Hu(1:end-1);
            sys_EDMDc3 = SparseRegression(Y,Yp,U,dt, options_method, lambda);

            % SINDYc (not xmean-substracted)
            options_method.order  = 3;
            Xi = NonlinearSparseRegression(x,u',dt,options_method,lambda);

            %% Prediction on training data
            Nt = M-1;
            % sparse DMDc
            options_method.order  = 1;
            x0train = poolData(Hx(:,1)',nVars,options_method.order,0)';
            x0train = x0train(2:end);
            [xDMDc,~] = lsim(sys_DMDc,Hu,tspan(1:Nt),x0train);
            xDMDc = xDMDc(:,end-1:end);
            xDMDc = xDMDc + repmat(xmean,[Nt 1]);
            xPredTrain{1}.x = xDMDc;

            % sparse EDMDc
            options_method.order  = 2;
            x0train = poolData(Hx(:,1)',nVars,options_method.order,0)';
            x0train = x0train(2:end);
            [xEMD2,~] = lsim(sys_EDMDc2,Hu,tspan(1:Nt),x0train);
            xEMD2 = xEMD2(:,1:2);
            xEMD2 = xEMD2 + repmat(xmean,[Nt 1]);
            xPredTrain{2}.x = xEMD2;

            % sparse EDMDc
            options_method.order  = 3;
            x0train = poolData(Hx(:,1)',nVars,options_method.order,0)';
            x0train = x0train(2:end);
            [xEMD3,~] = lsim(sys_EDMDc3,Hu,tspan(1:Nt),x0train);
            xEMD3 = xEMD3(:,1:2);
            xEMD3 = xEMD3 + repmat(xmean,[Nt 1]);
            xPredTrain{3}.x = xEMD3;

            % SINDYc
            p.ahat = Xi(:,1:2);
            p.polyorder = options_method.order; p.usesine = options_method.usesine; p.dt = dt;
            [Ns,N] = size(Hx);
            xSINDYc = zeros(Ns,N); xSINDYc(:,1) = x(1,:)';
            for ct=1:N-1
                xSINDYc(:,ct+1) = rk4u(@sparseGalerkinControl_Discrete,xSINDYc(:,ct),DataTrain.u(ct),dt,1,[],p);
            end
            xSINDYc = xSINDYc';
            xPredTrain{4}.x = xSINDYc;

            % Store results
            for iM = 1:N_MODEL
                PredictionResults_Training(:,:,N_ETA,N_LAMBDA,iM,N_REP) = xPredTrain{iM}.x;
            end

            %             clear ph
            %             figure,box on,
            %             ccolors = get(gca,'colororder');
            %             ph(1) = plot(tspan(1:Nt),x(1:Nt,1),'-','Color','k','LineWidth',1); hold on
            %             plot(tspan(1:Nt),x(1:Nt,2),'-','Color','k','LineWidth',1);
            %             for iM = 1:N_MODEL
            %                 for iVar = 1:nVars
            %                     ph(iM+1) = plot(tspan(1:Nt),xPredTrain{iM}.x(:,iVar),'--','Color',ccolors(iM,:),'LineWidth',1); hold on
            %                 end
            %             end
            %             legend(ph,model_collection)

            for iM = 1:N_MODEL
                TrainingError(iEta,iLambda,iM,iR) = norm(x(1:Nt,:)-xPredTrain{iM}.x,'fro');
            end

            %disp(['Training Error: ', num2str(squeeze(TrainingError(iEta,iLambda,:,iR))')])
            %             disp(['                DMDc: ', num2str(DMDc_err)])
            %             disp(['          sparseDMDc: ', num2str(sparseDMDc_err)])
            %% Prediction on validation data
            Nt = M-1;
            % sparse DMDc
            options_method.order  = 1;
            x0train = poolData(Hx(:,1)',nVars,options_method.order,0)';
            x0train = x0train(2:end);
            [xDMDc,~] = lsim(sys_DMDc,Hu,tspan(1:Nt),x0train);
            xDMDc = xDMDc(:,end-1:end);
            xDMDc = xDMDc + repmat(xmean,[Nt 1]);
            xPredTrain{1}.x = xDMDc;

            % sparse EDMDc
            options_method.order  = 2;
            x0train = poolData(Hx(:,1)',nVars,options_method.order,0)';
            x0train = x0train(2:end);
            [xEMD2,~] = lsim(sys_EDMDc2,Hu,tspan(1:Nt),x0train);
            xEMD2 = xEMD2(:,1:2);
            xEMD2 = xEMD2 + repmat(xmean,[Nt 1]);
            xPredTrain{2}.x = xEMD2;

            % sparse EDMDc
            options_method.order  = 3;
            x0train = poolData(Hx(:,1)',nVars,options_method.order,0)';
            x0train = x0train(2:end);
            [xEMD3,~] = lsim(sys_EDMDc3,Hu,tspan(1:Nt),x0train);
            xEMD3 = xEMD3(:,1:2);
            xEMD3 = xEMD3 + repmat(xmean,[Nt 1]);
            xPredTrain{3}.x = xEMD3;

            % SINDYc
            p.ahat = Xi(:,1:2);
            p.polyorder = options_method.order; p.usesine = options_method.usesine; p.dt = dt;
            [Ns,N] = size(Hx);
            xSINDYc = zeros(Ns,N); xSINDYc(:,1) = x(1,:)';
            for ct=1:N-1
                xSINDYc(:,ct+1) = rk4u(@sparseGalerkinControl_Discrete,xSINDYc(:,ct),DataTrain.u(ct),dt,1,[],p);
            end
            xSINDYc = xSINDYc';
            xPredTrain{4}.x = xSINDYc;

            % Store results
            for iM = 1:N_MODEL
                PredictionResults_Validation(:,:,N_ETA,N_LAMBDA,iM,N_REP) = xPredTrain{iM}.x;
            end

            for iM = 1:N_MODEL
                ValidationError(iEta,iLambda,iM,iR) = norm(x(1:Nt,:)-xPredTrain{iM}.x,'fro');
            end

            counter = counter + 1;
        end
    end

    tstop = toc(tstart);
    disp(['STATUS (',num2str(iR),' of ', num2str(N_REP),'): ', num2str(100*counter/(N_ETA*N_LAMBDA)),'%', '[Time=',num2str(tstop),'s]'])
end

%% STATS
idxnan = isnan(TrainingError(:));
maxerr = max(TrainingError(:));
TrainingError(idxnan) = maxerr; % change later
ErrStats.mean = zeros(N_ETA,N_LAMBDA,N_MODEL);
ErrStats.median = zeros(N_ETA,N_LAMBDA,N_MODEL);
ErrStats.std = zeros(N_ETA,N_LAMBDA,N_MODEL);
ErrStats.min = zeros(N_ETA,N_LAMBDA,N_MODEL);
for iEta = 1:N_ETA
    for iLambda = 1:N_LAMBDA
        for iModel = 1:N_MODEL
            ErrStats.mean(iEta,iLambda,iModel) = mean(TrainingError(iEta,iLambda,iModel,:));
            ErrStats.median(iEta,iLambda,iModel) = median(TrainingError(iEta,iLambda,iModel,:));
            ErrStats.std(iEta,iLambda,iModel) = std(TrainingError(iEta,iLambda,iModel,:));
            ErrStats.min(iEta,iLambda,iModel) = min(TrainingError(iEta,iLambda,iModel,:));
        end
    end
end

%%
err_axis = [min(TrainingError(:)), 1e3];%max(TrainingError(:))
figure,
for iModel = 1:N_MODEL
    subplot(4,4,(iModel-1)*4+1)
    surf(lambda_vec,eta_vec,ErrStats.median(:,:,iModel)), shading interp
    xlabel('\lambda'), ylabel('\eta'), zlabel('median')
    %     zlim([0 1e6])
    caxis(err_axis);
    set(gca,'xscale','log')
    view(2)

    subplot(4,4,(iModel-1)*4+2)
    surf(lambda_vec,eta_vec,ErrStats.mean(:,:,iModel)) , shading interp
    xlabel('\lambda'), ylabel('\eta'), zlabel('mean')
    %     zlim([0 1e6])
    caxis(err_axis);
    set(gca,'xscale','log')
    view(2)

    subplot(4,4,(iModel-1)*4+3)
    surf(lambda_vec,eta_vec,ErrStats.std(:,:,iModel)) , shading interp
    xlabel('\lambda'), ylabel('\eta'), zlabel('std')
    %     zlim([0 1e6])
    caxis(err_axis);
    set(gca,'xscale','log')
    view(2)

    subplot(4,4,(iModel-1)*4+4)
    surf(lambda_vec,eta_vec,ErrStats.min(:,:,iModel)) , shading interp
    xlabel('\lambda'), ylabel('\eta'), zlabel('min')
    %     zlim([0 1e6])
    caxis(err_axis);
    set(gca,'xscale','log')
    view(2)
end
% for iEta = 1:N_ETA
%     for iLambda = 1:N_LAMBDA
%         plot3()
%     end
% end

%% TIme series stats
iEta = 1; iLambda = 1; iM = 4;
data = squeeze(PredictionResults_Training(:,:,iEta,iLambda,iM,N_REP));

xBmin = min(data(:,1:nVars,:),[],3);
xBmax = max(data(:,1:nVars,:),[],3);
clear ph
figure('visible','off'),box on, hold on,
ccolors = get(gca,'colororder');
plot([tB(1),tB(1)],[-25 65],':','Color',[0.4,0.4,0.4],'LineWidth',1.5)
plot([t(end),t(end)],[-25 65],':','Color',[0.4,0.4,0.4],'LineWidth',1.5)
ylim([-25 65])
t1 = text(5,55,'Training', 'FontSize',12);
t2 = text(5+tA(1),55,'Validation', 'FontSize',12);

X=[tB(2:end)',fliplr(tB(2:end)')];                %#create continuous x value array for plotting
Y=[xBmin(:,1)',flipud(xBmax(:,1))'];              %#create y values for out and then back
fillh1 = fill(X,Y,ccolors(1,:)-[0 0.2 0.2]);                  %#plot filled area
fillh1.EdgeColor = ccolors(1,:)-[0 0.2 0.2]; fillh1.FaceAlpha = 0.5;

X=[tB(2:end)',fliplr(tB(2:end)')];                %#create continuous x value array for plotting
Y=[xBmin(:,2)',flipud(xBmax(:,2))'];              %#create y values for out and then back
fillh = fill(X,Y,ccolors(2,:)-[0.1 0.2 0.09]);                  %#plot filled area
fillh.EdgeColor = ccolors(2,:)-[0.1 0.2 0.09]; fillh.FaceAlpha = 0.5;

X=[tB(2:end)',fliplr(tB(2:end)')];                %#create continuous x value array for plotting
Y=[xBmin(:,3)',flipud(xBmax(:,3))'];              %#create y values for out and then back
fillh = fill(X,Y,ccolors(3,:)-[0.1 0.2 0.09]);                  %#plot filled area
fillh.EdgeColor = ccolors(3,:)-[0.1 0.2 0.09]; fillh.FaceAlpha = 0.5;

%             plot(tB(2:end),xBmin(:,1),'-k','LineWidth',2)
%             plot(tB(2:end),xBmin(:,1),'--g','LineWidth',2)

if eps~=0
    ph(4) = plot(tspan,x(:,1),'-','Color',0.7*ones(1,3),'LineWidth',1); %ccolors(1,:)+[0.15 0.3 0.25]
    plot(tspan,x(:,2),'-','Color',0.7*ones(1,3),'LineWidth',1); %ccolors(2,:)+[0.15 0.3 0.25]
    plot(tspan,x(:,3),'-','Color',0.7*ones(1,3),'LineWidth',1);
    ph(1) = plot([DataTrain.t;tA],[DataTrain.x(:,1);xA(:,1)],'-','Color',ccolors(1,:),'LineWidth',1);
    ph(2) = plot([DataTrain.t;tA],[DataTrain.x(:,2);xA(:,2)],'-','Color',ccolors(2,:),'LineWidth',1);
    ph(3) = plot([DataTrain.t;tA],[DataTrain.x(:,3);xA(:,3)],'-','Color',ccolors(3,:),'LineWidth',1);
else
    ph(1) = plot([DataTrain.t;tA],[DataTrain.x(:,1);xA(:,1)],'-','Color',ccolors(1,:),'LineWidth',1);
    ph(2) = plot([DataTrain.t;tA],[DataTrain.x(:,2);xA(:,2)],'-','Color',ccolors(2,:),'LineWidth',1);
    ph(3) = plot([DataTrain.t;tA],[DataTrain.x(:,3);xA(:,3)],'-','Color',ccolors(3,:),'LineWidth',1);
    ph(4) = plot(t,x(:,1),'--','Color',[0 1 0],'LineWidth',1); % Training data
    plot(t,x(:,2),'--','Color',[0 1 0],'LineWidth',1);
    plot(t,x(:,3),'--','Color',[0 1 0],'LineWidth',1);
end

%             ph(4) = plot(tB,xB(:,1),'-.','Color',ccolors(1,:)-[0 0.2 0.2],'LineWidth',2);
%             ph(5) = plot(tB,xB(:,2),'-.','Color',ccolors(2,:)-[0.1 0.2 0.09],'LineWidth',2);
grid off
xlim([0 tv(end)])
xlabel('Time')
ylabel('xi')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto');

if eps~=0
    lh = legend([ph([1,4]),fillh1],'Truth','Training',ModelName,'Location','NorthWest');
    %             lh.Position = [lh.Position(1)+0.13,lh.Position(2)-0.2,lh.Position(3:4)];
    lh.Position = [lh.Position(1)+0.02,lh.Position(2)-0.06,lh.Position(3:4)];
else
    lh = legend(ph([1,4]),'Truth','Training',ModelName,'Location','NorthWest');
    %             lh.Position = [lh.Position(1)+0.13,lh.Position(2)-0.2,lh.Position(3:4)];
    lh.Position = [lh.Position(1)+0.02,lh.Position(2)-0.06,lh.Position(3:4)];
end

print('-depsc2', '-painters','-loose','-cmyk', [figpath,'EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_Validation_N',sprintf('%04g',Ntrain_vec(iN)),'_Eta',sprintf('%03g',100*eta_vec(iNoise)),'_STATS.eps']);

delete(lh), delete(t1), delete(t2)
%             print('-depsc2', '-painters','-loose','-cmyk', [figpath,'EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_Validation_N',sprintf('%04g',Ntrain_vec(iN)),'_Eta',sprintf('%03g',100*eta_vec(iNoise)),'_STATS_noleg.eps']);

return
%% Show validation
clear ph
figure,box on,
ccolors = get(gca,'colororder');
ph(1) = plot(tspan,x(:,1),'-','Color','k','LineWidth',1); hold on
ph(2) = plot(tspan,x(:,2),'-','Color','k','LineWidth',1);
ph(3) = plot(tspan(Ndelay:Nt+Ndelay-1),xDMDc(:,1),'--','Color','r','LineWidth',2);
ph(4) = plot(tspan(Ndelay:Nt+Ndelay-1),xDMDc(:,2),'--','Color','r','LineWidth',2);
ph(5) = plot(tspan(Ndelay:Nt+Ndelay-1),xDMDc_sparse(:,1),':','Color','b','LineWidth',2);
ph(6) = plot(tspan(Ndelay:Nt+Ndelay-1),xDMDc_sparse(:,2),':','Color','b','LineWidth',2);
xlim([0 100]), ylim([0,120]);
xlabel('Time')
ylabel('Population size')
% legend('Prey (True)','Predator (True)', 'Prey (DMDc)','Predator (DMDc)')
legend(ph([1,3,5]),'True','DMDc',ModelName)
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
print('-depsc2', '-loose', '-cmyk', [figpath,'EX_LOTKA_SI_',ModelName,'_',InputSignalType,'.eps']);

DMDc_err = norm(x(1:end-1,:)-xDMDc,'fro');
sparseDMDc_err = norm(x(1:end-1,:)-xDMDc_sparse,'fro');
disp(['Training Error: '])
disp(['                DMDc: ', num2str(DMDc_err)])
disp(['          sparseDMDc: ', num2str(sparseDMDc_err)])

%% Prediction
% Reference
tspanV   = [100:dt:200];
xA      = xv;
tA      = tv;

% Model DMDc
x0      = [x(end,1:2)];
Hunew   = [u(end),uv(1:end)];
[xB,tB] = lsim(sysmodel_DMDc,Hunew,tspanV,[x0-[xmean]]');
xB = xB + repmat(xmean,[length(tB) 1]);

% Model sparseDMDc / sparseEDMDc
if Nstates == 2
    [xB_sparse,tB] = lsim(sysmodel_sparseDMDc,Hunew,tspanV,[x0-[xmean]]');
    xB_sparse = xB_sparse(:,end-1:end);
    xB_sparse = xB_sparse + repmat(xmean,[length(tB) 1]);
elseif Nstates > 2
    x0dmdc = [x0(1)-xmean(1); x0(2)-xmean(2); (x0(1)-xmean(1)).^2; (x0(1)-xmean(1)).*(x0(2)-xmean(2)); (x0(2)-xmean(2)).^2];
    [xB_sparse,tB] = lsim(sysmodel_sparseDMDc,Hunew,tspanV,x0dmdc);
    xB_sparse = xB_sparse(:,1:2);
    xB_sparse = xB_sparse + repmat(xmean,[length(tB) 1]);
end

DMDc_err = norm(xA-xB(1:end-1,:),2);
sparseDMDc_err = norm(xA-xB_sparse(1:end-1,:),2);
disp(['Validation Error: '])
disp(['                DMDc: ', num2str(DMDc_err)])
disp(['          sparseDMDc: ', num2str(sparseDMDc_err)])
%% Show training and prediction
VIZ_SI_Validation

%% Save Data
Model.name = 'sparseDMDc';
Model.sys = sysmodel_DMDc;
Model.Ndelay = Ndelay;
Model.dt = dt;
save(fullfile(datapath,['EX_LOTKA_SI_',ModelName,'_',InputSignalType,'.mat']),'Model')

author: eurika-kaiser

repo: SINDY-MPC


text:
function [output,status] = urlread(urlChar,method,params)
%URLREAD Returns the contents of a URL as a string.
%   S = URLREAD('URL') reads the content at a URL into a string, S.  If the
%   server returns binary data, the string will contain garbage.
%
%   S = URLREAD('URL','method',PARAMS) passes information to the server as
%   part of the request.  The 'method' can be 'get', or 'post' and PARAMS is a
%   cell array of param/value pairs.
%
%   [S,STATUS] = URLREAD(...) catches any errors and returns 1 if the file
%   downloaded successfully and 0 otherwise.
%
%   Examples:
%   s = urlread('http://www.mathworks.com')
%   s = urlread('ftp://ftp.mathworks.com/README')
%   s = urlread(['file:///' fullfile(prefdir,'history.m')])
%
%   From behind a firewall, use the Preferences to set your proxy server.
%
%   See also URLWRITE.

%   Matthew J. Simoneau, 13-Nov-2001
%   Copyright 1984-2011 The MathWorks, Inc.
%   $Revision: 1.3.2.12 $ $Date: 2011/09/03 22:43:00 $

% This function requires Java.
if ~usejava('jvm')
   error(message('MATLAB:urlread:NoJvm'));
end

import com.mathworks.mlwidgets.io.InterruptibleStreamCopier;

% Be sure the proxy settings are set.
com.mathworks.mlwidgets.html.HTMLPrefs.setProxySettings

% Check number of inputs and outputs.
if ~ischar(urlChar)
    error('MATLAB:urlread:InvalidInput','The first input, the URL, must be a character array.');
end
if (nargin > 1) && ~strcmpi(method,'get') && ~strcmpi(method,'post')
    error('MATLAB:urlread:InvalidInput','Second argument must be either "get" or "post".');
end

% Do we want to throw errors or catch them?
if nargout == 2
    catchErrors = true;
else
    catchErrors = false;
end

% Set default outputs.
output = '';
status = 0;

% GET method.  Tack param/value to end of URL.
if (nargin > 1) && strcmpi(method,'get')
    if mod(length(params),2) == 1
        error('MATLAB:urlread:InvalidInput','Invalid parameter/value pair arguments.');
    end
    for i=1:2:length(params)
        if (i == 1), separator = '?'; else separator = '&'; end
        param = char(java.net.URLEncoder.encode(params{i}));
        value = char(java.net.URLEncoder.encode(params{i+1}));
        urlChar = [urlChar separator param '=' value];
    end
end

% Create a urlConnection.
[urlConnection,errorid,errormsg] = plugin_urlreadwrite(mfilename,urlChar);
if isempty(urlConnection)
    if catchErrors, return
    else error(errorid,errormsg);
    end
end
urlConnection.setReadTimeout(5000); % timeout in 5 seconds

% POST method.  Write param/values to server.
if (nargin > 1) && strcmpi(method,'post')
    try
        urlConnection.setDoOutput(true);
        urlConnection.setRequestProperty( ...
            'Content-Type','application/x-www-form-urlencoded');
        printStream = java.io.PrintStream(urlConnection.getOutputStream);
        for i=1:2:length(params)
            if (i > 1), printStream.print('&'); end
            param = char(java.net.URLEncoder.encode(params{i}));
            value = char(java.net.URLEncoder.encode(params{i+1}));
            printStream.print([param '=' value]);
        end
        printStream.close;
    catch
        if catchErrors, return
        else error('MATLAB:urlread:ConnectionFailed','Could not POST to URL.');
        end
    end
end

% Read the data from the connection.
try
    inputStream = urlConnection.getInputStream;
    byteArrayOutputStream = java.io.ByteArrayOutputStream;
    % This StreamCopier is unsupported and may change at any time.
    isc = InterruptibleStreamCopier.getInterruptibleStreamCopier;
    isc.copyStream(inputStream,byteArrayOutputStream);
    inputStream.close;
    byteArrayOutputStream.close;
    output = native2unicode(typecast(byteArrayOutputStream.toByteArray','uint8'),'UTF-8');
catch
    if catchErrors, return
    else error('MATLAB:urlread:ConnectionFailed','Error downloading URL. Your network connection may be down or your proxy settings improperly configured.');
    end
end

status = 1;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ value, ifault ] = betain ( x, p, q, beta )

%*****************************************************************************80
%
%% BETAIN computes the incomplete Beta function ratio.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 January 2008
%
%  Author:
%
%    Original FORTRAN77 version by KL Majumder, GP Bhattacharjee.
%    FORTRAN90 version by John Burkardt.
%
%  Reference:
%
%    KL Majumder, GP Bhattacharjee,
%    Algorithm AS 63:
%    The incomplete Beta Integral,
%    Applied Statistics,
%    Volume 22, Number 3, 1973, pages 409-411.
%
%  Parameters:
%
%    Input, real X, the argument, between 0 and 1.
%
%    Input, real P, Q, the parameters, which
%    must be positive.
%
%    Input, real BETA, the logarithm of the complete
%    beta function.
%
%    Output, real BETAIN, the value of the incomplete
%    Beta function ratio.
%
%    Output, integer IFAULT, error flag.
%    0, no error.
%    nonzero, an error occurred.
%
  acu = 0.1D-14;

  value = x;
  ifault = 0;
%
%  Check the input arguments.
%
  if ( p <= 0.0 | q <= 0.0 )
    ifault = 1;
    return
  end

  if ( x < 0.0 | 1.0 < x )
    ifault = 2;
    return
  end
%
%  Special cases.
%
  if ( x == 0.0 | x == 1.0 )
    return
  end
%
%  Change tail if necessary and determine S.
%
  psq = p + q;
  cx = 1.0 - x;

  if ( p < psq * x )
    xx = cx;
    cx = x;
    pp = q;
    qq = p;
    indx = 1;
  else
    xx = x;
    pp = p;
    qq = q;
    indx = 0;
  end

  term = 1.0;
  ai = 1.0;
  value = 1.0;
  ns = floor ( qq + cx * psq );
%
%  Use the Soper reduction formula.
%
  rx = xx / cx;
  temp = qq - ai;
  if ( ns == 0 )
    rx = xx;
  end

  while ( 1 )

    term = term * temp * rx / ( pp + ai );
    value = value + term;
    temp = abs ( term );

    if ( temp <= acu & temp <= acu * value )

      value = value * exp ( pp * log ( xx ) ...
      + ( qq - 1.0 ) * log ( cx ) - beta ) / pp;

      if ( indx )
        value = 1.0 - value;
      end

      break

    end

    ai = ai + 1.0;
    ns = ns - 1;

    if ( 0 <= ns )
      temp = qq - ai;
      if ( ns == 0 )
        rx = xx;
      end
    else
      temp = psq;
      psq = psq + 1.0;
    end

  end

  return
end


text:
function s1= struct_sloppycat(s1, s2, varargin)
% STRUCT_SLOPPYCAT - Copy fields
%
% Synopsis:
%   s1= file_loadMatlab(s1,s2,<OPT>)
%
%Arguments:
%   S1:     STRUCT - target struct in which fields from s2 are copied
%   S2:     STRUCT
%   OPT: PROPLIST - Structure or property/value list of optional properties:
%     'Dim' - INT (default 2): ?
%     'Keepfields' -
%     'Matchsize' -
%
% See also: struct_copyFields

props = {'Dim',         2           '!INT[1]';
         'Keepfields'   3           '!INT[1]|!CHAR(none first last all)';
         'Matchsize'    0           '!BOOL';
         };

if nargin==0
    s1=props;return
end

if isempty(s1),
  s1= s2;
  return;
elseif isempty(s2),
  return;
end

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);
misc_checkType(s1,'STRUCT');
misc_checkType(s2,'STRUCT');

flds1= fieldnames(s1);
flds2= fieldnames(s2);

switch(opt.Keepfields),
 case {0,'none'},
  flds= intersect(flds1, flds2,'legacy');
  s1= rmfield(s1, setdiff(flds1, flds,'legacy'));
  s2= rmfield(s2, setdiff(flds2, flds,'legacy'));
 case {1,'first'},
  flds= flds1;
  if opt.Matchsize,
    s2= struct_createfields(s2, setdiff(flds, flds2,'legacy'), 'matchsize',s1(1));
  else
    s2= struct_createfields(s2, setdiff(flds, flds2,'legacy'));
  end
 case {2,'last'},
  flds= flds2;
  if opt.Matchsize,
    s1= struct_createfields(s1, setdiff(flds, flds1,'legacy'), 'matchsize',s2(1));
  else
    s1= struct_createfields(s1, setdiff(flds, flds1,'legacy'));
  end
 case {3,'all'},
  flds= union(flds1, flds2,'legacy');
  if opt.Matchsize,
    s1= struct_createfields(s1, setdiff(flds, flds1,'legacy'), 'matchsize',s2(1));
    s2= struct_createfields(s2, setdiff(flds, flds2,'legacy'), 'matchsize',s1(1));
  else
    s1= struct_createfields(s1, setdiff(flds, flds1,'legacy'));
    s2= struct_createfields(s2, setdiff(flds, flds2,'legacy'));
  end
 otherwise
  error('unknown value for OPT.Keepfields');
end

s1= cat(opt.Dim, s1, s2);

author: bbci


text:
function Q = updateQTable(state, a, reward, Q, next_state, alpha, gamma)

	% Update values in the Q-table

    Q(state, a) = Q(state,a) + alpha*(reward + gamma* max(Q(next_state,:)) - Q(state,a));
end

author: kennydl

repo: Reinforcment-Learning-With-Q-Learning


text:
%% (Internal) Design a filter to downsample signals prior printing to a report
%
%   b = design_downsample_filter( downsample_factor )
%
% Arguments:
%
%      + downsample_factor: times to downsample signal
%
% Output:
%
%      + b: filter coefficients
%
% Example:
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 21/4/2015
% Copyright 2008-2015
%
function b = design_downsample_filter( downsample_factor )
% MATLAB Code
% Generated by MATLAB(R) 8.2 and the Signal Processing Toolbox 6.20.
% Generated on: 23-Jun-2014 10:45:52

% FIR least-squares Lowpass filter designed using the FIRLS function.

% All frequency values are normalized to 1.

N     = 60;   % Order
Fpass = 1/downsample_factor;  % Passband Frequency
Fstop = min(1, Fpass + 0.1);  % Stopband Frequency
Wpass = 1;    % Passband Weight
Wstop = 100;  % Stopband Weight

% Calculate the coefficients using the FIRLS function.
b  = firls(N, [0 Fpass Fstop 1], [1 1 0 0], [Wpass Wstop]);

% [EOF]

author: marianux

repo: ecg-kit

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


text:
function minc = spm_cfg_minc
% SPM Configuration file for MINC Import
%__________________________________________________________________________
% Copyright (C) 2005-2011 Wellcome Trust Centre for Neuroimaging

% $Id: spm_cfg_minc.m 6952 2016-11-25 16:03:13Z guillaume $

%--------------------------------------------------------------------------
% data MINC files
%--------------------------------------------------------------------------
data         = cfg_files;
data.tag     = 'data';
data.name    = 'MINC files';
data.help    = {'Select the MINC files to convert.'};
data.filter  = 'mnc';
data.ufilter = '.*';
data.num     = [1 Inf];

%--------------------------------------------------------------------------
% dtype Data Type
%--------------------------------------------------------------------------
dtype         = cfg_menu;
dtype.tag     = 'dtype';
dtype.name    = 'Data Type';
dtype.help    = {'Data-type of output images. Note that the number of bits used determines the accuracy, and the amount of disk space needed.'};
dtype.labels  = {
                'UINT8   - unsigned char'
                'INT16   - signed short'
                'INT32   - signed int'
                'FLOAT32 - single prec. float'
                'FLOAT64 - double prec. float'
}';
dtype.values = {spm_type('uint8') spm_type('int16') spm_type('int32') ...
                spm_type('float32') spm_type('float64')};
dtype.val    = {spm_type('int16')};

%--------------------------------------------------------------------------
% ext Output image format
%--------------------------------------------------------------------------
ext         = cfg_menu;
ext.tag     = 'ext';
ext.name    = 'Output image format';
ext.help    = {'Output files can be written as .img + .hdr, or the two can be combined into a .nii file.'};
ext.labels = {
              'Two file (img+hdr) NIfTI'
              'Single file (nii) NIfTI'
}';
ext.values = { 'img', 'nii' };
ext.def    = @(val)spm_get_defaults('images.format', val{:});

%--------------------------------------------------------------------------
% opts Options
%--------------------------------------------------------------------------
opts         = cfg_branch;
opts.tag     = 'opts';
opts.name    = 'Options';
opts.val     = {dtype ext};
opts.help    = {'Conversion options'};

%--------------------------------------------------------------------------
% minc MINC Import
%--------------------------------------------------------------------------
minc         = cfg_exbranch;
minc.tag     = 'minc';
minc.name    = 'MINC Import';
minc.val     = {data opts};
minc.help    = {
    'MINC Conversion.'
    'MINC is the image data format used for exchanging data within the ICBM community, and the format used by the MNI software tools. It is based on NetCDF. MINC is no longer supported for reading images into SPM, so MINC files need to be converted to NIFTI format in order to use them. See http://www.bic.mni.mcgill.ca/software/ for more information.'
    }';
minc.prog    = @spm_run_minc;
minc.vout    = @vout;

%==========================================================================
function out = spm_run_minc(job)
for i=1:numel(job.data)
    spm_mnc2nifti(job.data{i},job.opts);
end

out.files = cell(size(job.data));
for i=1:numel(job.data)
    out.files{i} = spm_file(job.data{i}, 'path',pwd, 'ext',job.opts.ext);
end

%==========================================================================
function dep = vout(job)
dep            = cfg_dep;
dep.sname      = 'Converted Images';
dep.src_output = substruct('.','files');
dep.tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});

author: spm


text:
function a = rss_mp(alg,hyper)

%=================================================================
% RSC_MP reduced set selection
%=================================================================
% a = rss_mp(alg,hyper)
% generates a rss object, using the matching pursuit selection method
%
% hyperparameters:
% child=svm         algorithm worked on
% max_loops=1e5     maximum number of basis functions
% tolerance=1e-5    tolerated loss in ||w-w*||^2
% backfit=1         backfit on every nth iteration
% backfit_at_start=100  always backfit for first e.g. 100 iterations
% dont_revisit=1    dont return to old basis function optimization, always get new one
% reoptimize_b=1    recalculate the threshold b0
% alpha_cutoff=0    throw away svs with abs(alpha)<n
% bal_w=0           treat multiple w's as equal by normalizing by length
% optimizer='iterative' iterative update of matrix inverse
%
% model:
% alpha         new alphas for rs-vectors
% Xsv           rs vectors
% b0            the threshold
%
% stats:
% w2=0          final value of ||w-w*||^2
% res=[]        results on a separate test set
% dtst=[]       separate test set
% test_on=0     iterations to test on
%
% methods:
% train         constructs a reduced set, returns trained rs-machine
% test          tests new rs-machine on supplied data
%
% example:
% d=gen(toy2d('2circles','l=100'));
% [r,a]=train(svm({kernel('rbf',1),'C=10000','alpha_cutoff=1e-2'}),d);
% [r,a2]=train(rss_mp(a,'tolerance=1e-2'),d);
% test(a2,d,loss)
%
%=================================================================
% author: goekhan bakir, jason weston
% reference: fast binary and multi-output rss, 2004
%=================================================================

  %hyperparams
  a.child=svm;
  a.max_loops=100000;% maximum number of basis functions
  a.tolerance=1e-5;
  a.backfit=1;      % backfit on every nth iteration
  a.backfit_at_start=100; % always backfit for first e.g. 100 iterations
  a.dont_revisit=1;  % dont return to old basis function optimization, always get new one
  a.reoptimize_b=1;
  a.alpha_cutoff=0;
  a.bal_w=0;         % treat multiple w's as equal by normalizing by length
  a.optimizer='iterative'; % iterative update of matrix inverse

  % model
  a.alpha=[];
  a.Xsv=[];
  a.b0=0;

  % training / testing statistics
  a.w2=0;  % final value of ||w*-w^2||^2
  a.res=[];% results on a separate test set
  a.dtst=[];% separate test set
  a.test_on=0; % iterations to test on

  if nargin==0
    a.child=svm;
  else
    a.child=alg; %% algorithm to use
  end

  p=algorithm('rss_mp');
  a= class(a,'rss_mp',p);
  a.algorithm.use_signed_output=0;

  if nargin==2
    eval_hyper;
  end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef FiveLinkAnimator < Animator.AbstractAnimator
    properties
        leg1Color = 'r';
        leg2Color = 'b';
        torsoColor = 'k';
        groundColor = 'g';
    end

    properties (Access = private)
        ground;

        pH_RKnee;
        pH_LKnee;
        pRKnee_RFoot;
        pLKnee_LFoot;
        pHT;

        text;

        starting_index;
        next_frame_time;

        H;

        q_all;
        t_all;
    end

    properties
        updateWorldPosition logical
    end

    methods
        function obj = FiveLinkAnimator(t, q, varargin)
            obj = obj@Animator.AbstractAnimator(); % Calling super constructor

            obj.q_all = q;
            obj.t_all = t;

            obj.startTime = t(1);
            obj.currentTime = obj.startTime;
            obj.endTime = t(end);
            obj.updateWorldPosition = false;

            if isempty(varargin)
                [terrain.Tx, terrain.Ty] = meshgrid(-10:1:10, -10:1:10);
                terrain.Tz = 0.*terrain.Tx;
            else
                terrain = varargin{1};
            end

            Rz = @(th) [cos(th), -sin(th), 0; sin(th), cos(th), 0; 0,0,1];
            Ry = @(th) [cos(th), 0, sin(th); 0, 1, 0; -sin(th), 0, cos(th)];
            Rx = @(th) [1,0,0; 0, cos(th), -sin(th); 0, sin(th), cos(th)];

            r = obj.q_all(1:3, end) - obj.q_all(1:3, 1);
            th = obj.q_all(4:6, end) - obj.q_all(4:6, 1);
            obj.H = [Rx(th(1))*Ry(th(2))*Rz(th(3)), r; 0,0,0,1];

            % Initialization
            q = obj.q_all(:,1);

            pH = [q(1);0;q(2)];
            pT = p_Torso(q);
            pRK = p_q2_right(q);
            pLK = p_q2_left(q);
            pR = p_RightToe(q);
            pL= p_LeftToe(q);

            % Define Terrain
            obj.ground = surf(terrain.Tx,terrain.Ty,terrain.Tz); hold on;

            % Define links
            obj.pH_RKnee = line([0,0],[pH(1) pRK(1)],[pH(3) pRK(3)]);
            obj.pH_LKnee = line([0,0],[pH(1) pLK(1)],[pH(3) pLK(3)]);
            obj.pRKnee_RFoot = line([0,0],[pRK(1) pR(1)],[pRK(3) pR(3)]);
            obj.pLKnee_LFoot = line([0,0],[pLK(1) pL(1)],[pLK(3) pL(3)]);
            obj.pHT = line([0,0],[pH(1) pT(1)],[pH(3) pT(3)]);

            set(obj.pH_RKnee,'LineWidth',3,'Color',obj.leg1Color);
            set(obj.pH_LKnee,'LineWidth',3,'Color',obj.leg2Color);
            set(obj.pRKnee_RFoot,'LineWidth',3,'Color',obj.leg1Color);
            set(obj.pLKnee_LFoot,'LineWidth',3,'Color',obj.leg2Color);
            set(obj.pHT,'LineWidth',3,'Color',obj.torsoColor);

        end

        function Draw(obj, t, x)
            q = x;

            pH = [q(1);0;q(2)];
            pT = p_Torso(q);
            pRK = p_q2_right(q);
            pLK = p_q2_left(q);
            pR = p_RightToe(q);
            pL= p_LeftToe(q);

            set(obj.pH_RKnee,'YData',[pH(1) pRK(1)],'ZData',[pH(3) pRK(3)], 'XData',[0 0]);
            set(obj.pH_LKnee,'YData',[pH(1) pLK(1)],'ZData',[pH(3) pLK(3)],'XData',[0 0]);
            set(obj.pRKnee_RFoot,'YData',[pRK(1) pR(1)],'ZData',[pRK(3) pR(3)],'XData',[0 0]);
            set(obj.pLKnee_LFoot,'YData',[pLK(1) pL(1)],'ZData',[pLK(3) pL(3)],'XData',[0 0]);
            set(obj.pHT,'YData',[pH(1) pT(1)],'ZData',[pH(3) pT(3)],'XData',[0 0]);

            delete(obj.text);
            obj.text = text(0,pH(1),pH(3)+1,['t = ',sprintf('%.2f',t)]);
            obj.text.FontSize = 14;
            obj.text.FontWeight = 'Bold';
            obj.text.Color = [0,0,0];
            %  set(obj.text);

            drawnow;
        end

        function x = GetData(obj, t)
            t_start = obj.t_all(1);
            t_end = obj.t_all(end);
            delta_t = t_end - t_start;

            val = 0;

            if t < t_start || t > t_end
                val = floor((t - t_start) / delta_t);
                t = t - val * delta_t;
            end

            if t < t_start
                t = t_start;
            elseif t > t_end
                t = t_end;
            end

            n = length(obj.t_all);
            x = obj.q_all(:, 1); % Default

            a = 1;
            b = n;

            while (a <= b)
                c = floor((a + b) / 2);

                if t < obj.t_all(c)
                    x = obj.q_all(:, c);
                    b = c - 1;
                elseif t > obj.t_all(c)
                    a = c + 1;
                else
                    x = obj.q_all(:, c);
                    break;
                end
            end

            delta_q = obj.q_all(1:6, end) - obj.q_all(1:6, 1);

            %T = obj.H(1:3,4);
            %roll = atan2(-obj.H(2,3),obj.H(3,3));
            %pitch = asin(obj.H(1,3));
            %yaw = atan2(-obj.H(1,2),obj.H(1,1));

            if obj.updateWorldPosition
                Rz = @(th) [cos(th), -sin(th), 0; sin(th), cos(th), 0; 0,0,1];
                Ry = @(th) [cos(th), 0, sin(th); 0, 1, 0; -sin(th), 0, cos(th)];
                Rx = @(th) [1,0,0; 0, cos(th), -sin(th); 0, sin(th), cos(th)];

                x_orig_init = obj.q_all(1:6, 1);
                x_current_init = obj.q_all(1:6, 1);
                if val > 0
                    for i = 1:val
                        x_end = obj.q_all(1:6, end);

                        r1 = x_current_init(1:3) - x_orig_init(1:3);
                        th1 = x_current_init(4:6) - x_orig_init(4:6);
                        H1 = [Rx(th1(1))*Ry(th1(2))*Rz(th1(3)), r1; 0,0,0,1];

                        r2 = x_end(1:3) - x_orig_init(1:3);
                        th2 = x_end(4:6) - x_orig_init(4:6);
                        H2 = [Rx(th2(1))*Ry(th2(2))*Rz(th2(3)), r2; 0,0,0,1];

                        H = H1*H2;
                        T = H(1:3,4);
                        roll = atan2(-H(2,3),H(3,3));
                        pitch = asin(H(1,3));
                        yaw = atan2(-H(1,2),H(1,1));

                        x_current_init = x_orig_init(1:6) + [T;roll;pitch;yaw];
                    end

                    x_current = x(1:6);

                    r1 = x_current_init(1:3) - x_orig_init(1:3);
                    th1 = x_current_init(4:6) - x_orig_init(4:6);
                    H1 = [Rx(th1(1))*Ry(th1(2))*Rz(th1(3)), r1; 0,0,0,1];

                    r2 = x_current(1:3) - x_orig_init(1:3);
                    th2 = x_current(4:6) - x_orig_init(4:6);
                    H2 = [Rx(th2(1))*Ry(th2(2))*Rz(th2(3)), r2; 0,0,0,1];

                    H = H1*H2;
                    T = H(1:3,4);
                    roll = atan2(-H(2,3),H(3,3));
                    pitch = asin(H(1,3));
                    yaw = atan2(-H(1,2),H(1,1));

                    x(1:6) = x_orig_init(1:6) + [T;roll;pitch;yaw];
                end
            end
        end

        function [center, radius, yaw] = GetCenter(obj, t, x)
            q = x;

            center = [0; q(1:2)];
            radius = 1.5;
            yaw = 0;
        end
    end
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev