text:
function MDP_DEM_Mixed_Models_Movement
% This demo illustrates a series of computational pathologies as elicited
% during a synthetic neurological examination. This focuses upon an
% examination of the biceps reflex, and a simple coordination task. Each of
% these are simulated for an arm that may move in three dimensions through
% internal and external rotation of the shoulder, flexion and extension of
% the shoulder, and flexion and extension of the elbow. These dynamics play
% out through an active Bayesian filtering scheme, where a series of
% attracting points draw the hand to different locations in 3D space. The
% selection of these attracting points involves a hierarhical Markov
% Decision Process, which identifies these sequences based upon the prior
% belief that (1) the sequence will minimise expected free energy and (2)
% the sequence is consistent with the trajectory predicted by the highest
% (contextual) level of the model.
%__________________________________________________________________________

rng default

% REFLEXES
%==========================================================================
% Demonstration 1 - Tendon reflexes (healthy, pyramidal lesion,
% cerebellar lesion)

% Specify generative model and process
%--------------------------------------------------------------------------

a = 2;                                      % length of upper arm
b = 1.5;                                    % length of forearm
c = [-1;-1;-1];                             % location of shoulder

DEM        = DEM_MDP_Movement_Disorders_GM; % Continuous state generative model
x = [0;-1;1;0;0;0];                         % Initial joint positions

DEM.G(1).g = @(x,v,a,P) Gg1Reflex(x,v,a,P); % Generative process allows for tendon tap
DEM.M(1).x = x;
DEM.G(1).x = x;
DEM.U = [MDP_DEM_MD_transform(x(1:3),a,b,c);0;0;0]*ones(1,64);
DEM.C = zeros(size(DEM.U));
DEM.C(3,10:12) = -0.25;                     % Transient stimulation of type II sensory afferents

% Solve model and plot healthy reflexes
%--------------------------------------------------------------------------
dem = spm_ADEM(DEM);

% Animate
spm_figure('GetWin','Figure 1'); clf
opengl('software')
for i = 1:size(dem.pU.x{1},2)
    MDP_DEM_MD_plot_arm(full(dem.pU.x{1}(1:3,i)),dem.M(1).pE.b1,dem.M(1).pE.b2,[-1;-1;-1]), hold on
    light
    view(10,0)
    zlim([-4 1])
    xlim([-1.2 2.1])
    drawnow
    hold off
    title('Normal')
end

clear dem

% Upper motor lesion
%--------------------------------------------------------------------------
DEM.M(1).V = exp(5);        % Overestimate precision

dem = spm_ADEM(DEM);

% Animate
spm_figure('GetWin','Figure 1'); clf
for i = 1:size(dem.pU.x{1},2)
    MDP_DEM_MD_plot_arm(full(dem.pU.x{1}(1:3,i)),dem.M(1).pE.b1,dem.M(1).pE.b2,[-1;-1;-1]), hold on
    light
    view(10,0)
    zlim([-4 1])
    xlim([-1.2 2.1])
    drawnow
    hold off
    title('Pyramidal')
end

clear dem

% Cerebellar
%--------------------------------------------------------------------------
DEM.M(1).V = exp(4);        % Restore sensory precision
DEM.M(1).E.s = 3/4;         % Overestimate smoothness

dem = spm_ADEM(DEM);

spm_figure('GetWin','Figure 1'); clf
for i = 1:size(dem.pU.x{1},2)
    MDP_DEM_MD_plot_arm(full(dem.pU.x{1}(1:3,i)),dem.M(1).pE.b1,dem.M(1).pE.b2,[-1;-1;-1]), hold on
    light
    view(10,0)
    zlim([-4 1])
    xlim([-1.2 2.1])
    drawnow
    hold off
    title('Cerebellar')
end

clear all

% COORDINATION
%==========================================================================
% Demonstration 2 - Mixed models (healthy)

% Specify discrete generatve model and mapping between continuous and
% discrete
%--------------------------------------------------------------------------

% Locations
%--------------------------------------------------------------------------

% 3D coordinates for three targets
L1  =   [   1;   -1; -3];
L2  =   [  -1;    1; -3];
L3  =   [   0;    0; -1];

% Intermediate (fictive) attracting points)
for i = 1:4
    L{i} = L1 + (i-1)*(L2 - L1)/4;
end
for i = 5:8
    L{i} = L2 + (i-5)*(L3 - L2)/4;
end
for i = 9:12
    L{i} = L3 + (i-9)*(L1 - L3)/4;
end

% MDP outcomes to DEM causes
%--------------------------------------------------------------------------
N  = 14;
nl = length(L);
nh = 3;

for i = 1:nh
    for j = 1:nl
        for k = 1:2
            c           = [L{j}; sparse(i,1,1,nh,1)];
            u           = [L{j}; sparse(i,1,1,nh,1)];
            demi.U{i,j,k} = u*ones(1,N);
            demi.C{i,j,k} = c*ones(1,N);
        end
    end
end

o     = [1 1 1];
O{1}  = spm_softmax(ones(nh,1));
O{2}  = spm_softmax(ones(nl,1));
O{3}  = spm_softmax(ones(2,1));

% generative model (Continuous)
%==========================================================================
DEM      = DEM_MDP_Movement_Disorders_GM;
DEM      = spm_MDP_DEM(DEM,demi,O,o);

% generative model (Discrete)
%==========================================================================
mdp      = MDP_DEM_Movement_Disorders_GM_1;
mdp.demi = demi;
mdp.DEM  = DEM;

mdp      = MDP_DEM_Movement_Disorders_GM_2(mdp);

% invert or solve
%--------------------------------------------------------------------------
rng default
MDP  = spm_MDP_check(mdp);
MDP  = spm_MDP_VB_X(MDP);

% illustrate belief updating - discrete
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 2'); clf
spm_MDP_VB_LFP(MDP,[],2);

% illustrate belief updating - continuous (last movement)
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 3'); clf;
spm_DEM_qU(MDP.mdp(end).dem(end).qU)

% illustrate movement
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 4'); clf

[Q, R, X, Y, Z] = MDP_DEM_Movement_Disorders_Animate(MDP,L);

spm_figure('GetWin','Figure 4'); clf
plot3(Q(1,:),Q(2,:),Q(3,:),'LineWidth',2,'Color','k'), hold on
for m = 1:length(MDP.mdp)
    for j = 1:length(MDP.mdp(m).dem)
        h = (m - 1)*length(MDP.mdp(m).dem) + j;
        H = length(MDP.mdp)*length(MDP.mdp(m).dem);
        surf(0.1*X + MDP.mdp(m).dem(j).qU.v{2}(1,i),0.1*Y + MDP.mdp(m).dem(j).qU.v{2}(2,i),0.1*Z + MDP.mdp(m).dem(j).qU.v{2}(3,2),'Facecolor','r','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'FaceAlpha',h/(1.5*H))
    end
end

light
axis equal
spm_figure('GetWin','Figure 5'); clf
subplot(5,1,1)
plot(0.05*(0:size(R,2)-1),R' - repmat(R(:,1)',size(R,2),1))
title('Normal')

% Demonstration 3 - Mixed models (Pyramidal versus extrapyramidal)
%--------------------------------------------------------------------------
% Pyramidal

mdp.MDP.DEM.M(1).V = exp(5); % Overestimate precision

% invert or solve
%--------------------------------------------------------------------------
rng default
MDP  = spm_MDP_check(mdp);
MDP  = spm_MDP_VB_X(MDP);

spm_figure('GetWin','Figure 6'); clf

[Q, R, X, Y, Z] = MDP_DEM_Movement_Disorders_Animate(MDP,L);

spm_figure('GetWin','Figure 6'); clf
plot3(Q(1,:),Q(2,:),Q(3,:),'LineWidth',2,'Color','k'), hold on
for m = 1:length(MDP.mdp)
    for j = 1:length(MDP.mdp(m).dem)
        h = (m - 1)*length(MDP.mdp(m).dem) + j;
        H = length(MDP.mdp)*length(MDP.mdp(m).dem);
        surf(0.1*X + MDP.mdp(m).dem(j).qU.v{2}(1,i),0.1*Y + MDP.mdp(m).dem(j).qU.v{2}(2,i),0.1*Z + MDP.mdp(m).dem(j).qU.v{2}(3,2),'Facecolor','r','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'FaceAlpha',h/(1.5*H))
    end
end
light
axis equal

spm_figure('GetWin','Figure 5');
subplot(5,1,2)
plot(0.05*(0:size(R,2)-1),R' - repmat(R(:,1)',size(R,2),1))
title('Pyramidal')

% Extrapyramidal
%--------------------------------------------------------------------------
mdp.MDP.DEM.M(1).V = exp(4);
mdp.MDP.beta = 50;

% invert or solve
%--------------------------------------------------------------------------
rng default
MDP  = spm_MDP_check(mdp);
MDP  = spm_MDP_VB_X(MDP);

% illustrate movement
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 8'); clf

[Q, R, X, Y, Z] = MDP_DEM_Movement_Disorders_Animate(MDP,L);

spm_figure('GetWin','Figure 8'); clf
plot3(Q(1,:),Q(2,:),Q(3,:),'LineWidth',2,'Color','k'), hold on
for m = 1:length(MDP.mdp)
    for j = 1:length(MDP.mdp(m).dem)
        h = (m - 1)*length(MDP.mdp(m).dem) + j;
        H = length(MDP.mdp)*length(MDP.mdp(m).dem);
        surf(0.1*X + MDP.mdp(m).dem(j).qU.v{2}(1,i),0.1*Y + MDP.mdp(m).dem(j).qU.v{2}(2,i),0.1*Z + MDP.mdp(m).dem(j).qU.v{2}(3,2),'Facecolor','r','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'FaceAlpha',h/(1.5*H))
    end
end
light
axis equal

spm_figure('GetWin','Figure 5');
subplot(5,1,3)
plot(0.05*(0:size(R,2)-1),R' - repmat(R(:,1)',size(R,2),1))
title('Extrapyramidal')

% Cerebellar
%--------------------------------------------------------------------------
mdp.MDP.beta = 16;
mdp.MDP.DEM.M(1).E.s = 3/4;

% invert or solve
%--------------------------------------------------------------------------
rng default
MDP  = spm_MDP_check(mdp);
MDP  = spm_MDP_VB_X(MDP);

% illustrate movement
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 9'); clf

[Q, R, X, Y, Z] = MDP_DEM_Movement_Disorders_Animate(MDP,L);

spm_figure('GetWin','Figure 9'); clf
plot3(Q(1,:),Q(2,:),Q(3,:),'LineWidth',2,'Color','k'), hold on
for m = 1:length(MDP.mdp)
    for j = 1:length(MDP.mdp(m).dem)
        h = (m - 1)*length(MDP.mdp(m).dem) + j;
        H = length(MDP.mdp)*length(MDP.mdp(m).dem);
        surf(0.1*X + MDP.mdp(m).dem(j).qU.v{2}(1,i),0.1*Y + MDP.mdp(m).dem(j).qU.v{2}(2,i),0.1*Z + MDP.mdp(m).dem(j).qU.v{2}(3,2),'Facecolor','r','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'FaceAlpha',h/(1.5*H))
    end
end
light
axis equal

spm_figure('GetWin','Figure 5');
subplot(5,1,4)
plot(0.05*(0:size(R,2)-1),R' - repmat(R(:,1)',size(R,2),1))
title('Cerebellar')

% Executive
%--------------------------------------------------------------------------
mdp.MDP.DEM.M(1).E.s = 1/2;
mdp.A{3} = ones(size(mdp.A{3})); % Disconnect higher level

% invert or solve
%--------------------------------------------------------------------------
rng default
MDP  = spm_MDP_check(mdp);
MDP  = spm_MDP_VB_X(MDP);

% illustrate movement
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 10'); clf

[Q, R, X, Y, Z] = MDP_DEM_Movement_Disorders_Animate(MDP,L);

spm_figure('GetWin','Figure 10'); clf
plot3(Q(1,:),Q(2,:),Q(3,:),'LineWidth',2,'Color','k'), hold on
for m = 1:length(MDP.mdp)
    for j = 1:length(MDP.mdp(m).dem)
        h = (m - 1)*length(MDP.mdp(m).dem) + j;
        H = length(MDP.mdp)*length(MDP.mdp(m).dem);
        surf(0.1*X + MDP.mdp(m).dem(j).qU.v{2}(1,i),0.1*Y + MDP.mdp(m).dem(j).qU.v{2}(2,i),0.1*Z + MDP.mdp(m).dem(j).qU.v{2}(3,2),'Facecolor','r','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'FaceAlpha',h/(1.5*H))
    end
end
light
axis equal

spm_figure('GetWin','Figure 5');
subplot(5,1,5)
plot(0.05*(0:size(R,2)-1),R' - repmat(R(:,1)',size(R,2),1))
title('Executive')

function MDP = MDP_DEM_Movement_Disorders_GM_1

% First level
%--------------------------------------------------------------------------

% Initial states
%--------------------------------------------------------------------------
D{1} = ones(3,1); % Location of visual cue
D{2} = zeros(12,1);  % Location of arm
D{2}(1) = 1;

% Likelihood
%--------------------------------------------------------------------------
for f1 = 1:numel(D{1})
    for f2 = 1:numel(D{2})
        A{1}(f1,f1,f2) = 1; % Vision (target)
        A{2}(f2,f1,f2) = 1; % Vision + Proprioception (arm)
        A{3}(2,f1,f2)  = (f1 == ((f2-1)/4)+1);  % At target
        A{3}(1,f1,f2)  = ~(f1 == ((f2-1)/4)+1); % Not at target
    end
end

% Transitions
%--------------------------------------------------------------------------
B{1} = eye(3);

for k = 1:length(D{2})
    for j = 1:length(D{2})
        if k == j
            B{2}(k,j,1) = 1; % Stay still
        elseif (k - j) == 1
            B{2}(k,j,2) = 1; % Move to next point
        elseif (j - k) == 1
            B{2}(k,j,3) = 1; % Move to previous point
        end
    end
end

B{2}(1,end,2) = 1;
B{2}(end,1,3) = 1;

% Preferences
%--------------------------------------------------------------------------
C{1} = zeros(3,1);
C{2} = zeros(12,1);
C{3} = [-3 3]';

% Policies
%--------------------------------------------------------------------------
V(:,:,1) = ones(8,3); % (time, policy, factor)
V(:,:,2) = [1 2 3;
    1 2 3;
    1 2 3;
    1 2 3;
    1 1 1;
    1 1 1;
    1 1 1;
    1 1 1];

E = [1 0.5 0.5]';

% First level MDP
%--------------------------------------------------------------------------
MDP.A = A;
MDP.B = B;
MDP.C = C;
MDP.D = D;
MDP.E = E;
MDP.V = V;
MDP.T = 9;
MDP.chi = -exp(64);
MDP.beta = 16;
% MDP.beta = 50;
MDP     = spm_MDP_check(MDP);

function MDP = MDP_DEM_Movement_Disorders_GM_2(mdp)

% First level
%--------------------------------------------------------------------------
MDP.MDP = mdp;

% Second level
%--------------------------------------------------------------------------

% Initial states
%--------------------------------------------------------------------------
D{1} = ones(3,1); % Target location
D{2} = zeros(9,1); % Initial locations x policies at lower level ({L1,P1},{L1,P2},{L1,P3},{L2,P1},...)
D{2}([1 2 3]) = 1;

% Likelihood
%--------------------------------------------------------------------------
for f1 = 1:length(D{1})
    for f2 = 1:length(D{2})
        A{1}(f1,f1,f2) = 1; % Target location 1st level

        if f2 < 4
            A{2}(11,f1,f2) = 0.03;  % Start location 1st level
            A{2}(12,f1,f2) = 0.12;
            A{2}(1,f1,f2)  = 0.7;
            A{2}(2,f1,f2)  = 0.12;
            A{2}(3,f1,f2)  = 0.03;

        elseif f2 <7
            A{2}(3,f1,f2) = 0.03;
            A{2}(4,f1,f2) = 0.12;
            A{2}(5,f1,f2)  = 0.7;
            A{2}(6,f1,f2)  = 0.12;
            A{2}(7,f1,f2)  = 0.03;
        else
            A{2}(7,f1,f2) = 0.03;
            A{2}(8,f1,f2) = 0.12;
            A{2}(9,f1,f2)  = 0.7;
            A{2}(10,f1,f2)  = 0.12;
            A{2}(11,f1,f2)  = 0.03;
        end

        A{3}(1,f1,f2)  = (rem(f2+2,3) == 0) + 0.5;  % Policy 1st level (slight bias towards 'stay still' policy)
        A{3}(2,f1,f2)  = (rem(f2+1,3) == 0);
        A{3}(3,f1,f2)  = (rem(f2,3) == 0);
    end
end

% Transitions
%--------------------------------------------------------------------------
B{1}(:,:,1) = ones(length(D{1}));
B{1}(:,:,2) = ones(length(D{1})); % Action that causes no change to prevent entering HMM mode
I = eye(3);
B{2} = [I, [I(3,:);I(1:2,:)],([I(3,:);I(1:2,:)])']/3;
B{2} = kron(B{2},ones(3,1));

% Preferences
%--------------------------------------------------------------------------
C{1} = zeros(3,1);
C{2} = zeros(12,1);
C{3} = zeros(3,1);

% Second level MDP
%--------------------------------------------------------------------------
MDP.A = A;
MDP.B = B;
MDP.C = C;
MDP.D = D;
MDP.T = 4;

MDP.link  = [1 0 0;
    0 1 0];
MDP.linkE = [0 0 1];

MDP     = spm_MDP_check(MDP);

function DEM = DEM_MDP_Movement_Disorders_GM
% hidden states
%--------------------------------------------------------------------------
x   = [0;-1;0.5;0;0;0]; % shoulder rotation, flexion, elbow flexion (position and velocity)

% and actions
%--------------------------------------------------------------------------
a   = [0;0;0];

% parameters
%--------------------------------------------------------------------------
P.kv = 2;              % damping
P.a = a;
P.b1 = 2;
P.b2 = 1.5;

M(1).pE = P;

% recognition model
%==========================================================================
M(1).E.s = 1/2;                               % smoothness: default 1/2 [3/4 for lesion]
M(1).E.n = 4;                                 % order of
M(1).E.d = 2;                                 % generalised motion

% level 1
%--------------------------------------------------------------------------
M(1).f  = @(x,v,P) Mf1(x,v,P);                % plant dynamics
M(1).g  = @(x,v,P) Mg1(x,v,P);                % prediction

M(1).x  = x;                                  % hidden states
M(1).V  = exp(4);                             % error precision (g)[4] [3/8 ?]
M(1).W  = exp(2);                             % error precision (f) [2]

% level 2:
%--------------------------------------------------------------------------
M(2).v  = [0;0;0;0;0;0];                      % priors (attracting point)
M(2).V  = exp(15);

% generative model
%==========================================================================

G(1).E.s = 1/2;                               % smoothness: default 1/2
G(1).E.n = 4;                                 % order of
G(1).E.d = 2;                                 % generalised motion

% first level
%--------------------------------------------------------------------------
G(1).f  = @(x,v,a,P) Gf1(x,v,a,P);
G(1).g  = @(x,v,a,P) Gg1(x,v,a,P);
G(1).x  = x;                                  % hidden states
G(1).V = exp(16);
% G(1).V  = exp([4*ones(1,6) 16*ones(1,6)]);    % error precision [16]
G(1).W  = exp(16);                            % error precision
G(1).pE = P;

% second level
%--------------------------------------------------------------------------
G(2).v  = [0;0;0;0;0;0];                      % exogenous forces
G(2).a  = spm_vec(a);                         % action forces
G(2).V  = exp(16);

DEM.G = G;
DEM.M = M;

% Equations for DEM model and process
%==========================================================================
function f = Mf1(x,v,P)
a = P.b1; % length of upper arm
b = P.b2; % length of forearm
c = [-1;-1;-1]; % location of shoulder

f = [x(4:6);MDP_DEM_MD_dthetadt(v(1:3),x,a,b,c) - P.kv*x(4:6)];

function g = Mg1(x,v,P)
a = P.b1; % length of upper arm
b = P.b2; % length of forearm
c = [-1;-1;-1]; % location of shoulder

g = [x;MDP_DEM_MD_transform(x(1:3),a,b,c);v(4:6)]; % proprioceptive, visual (hand), visual (target)

function f = Gf1(x,~,a,P)
f = [x(4:6);a  - P.kv*x(4:6)];

function g = Gg1(x,v,~,P)
a = P.b1; % length of upper arm
b = P.b2; % length of forearm
c = [-1;-1;-1]; % location of shoulder

g = [x;MDP_DEM_MD_transform(x(1:3),a,b,c);v(4:6)]; % proprioceptive, visual (hand), visual (target)

function g = Gg1Reflex(x,v,~,P)
a = P.b1; % length of upper arm
b = P.b2; % length of forearm
c = [-1;-1;-1]; % location of shoulder

g = [x + [zeros(3,1);v(1:3)];MDP_DEM_MD_transform(x(1:3),a,b,c);v(4:6)]; % proprioceptive, visual (hand), visual (target)

% Additional routines (plotting/coordinate transforms)
%==========================================================================
function y = MDP_DEM_MD_transform(x,a,b,c)
% converts angles to hand position in euclidean space
% x = angle of shoulder rotation, flexion, elbow flexion

y1 = c + [a*cos(x(2))*cos(x(1));
    a*cos(x(2))*sin(x(1));
    a*sin(x(2))];

y = y1 + [b*cos(x(3)+x(2))*cos(x(1));
    b*cos(x(3)+x(2))*sin(x(1));
    b*sin(x(3)+x(2))];

function MDP_DEM_MD_plot_arm(x,a,b,c)

y1 = c + [a*cos(x(2))*cos(x(1));
    a*cos(x(2))*sin(x(1));
    a*sin(x(2))];

y = y1 + [b*cos(x(3)+x(2))*cos(x(1));
    b*cos(x(3)+x(2))*sin(x(1));
    b*sin(x(3)+x(2))];

[S1, S2, S3] = sphere;

surf(S1*0.2+c(1),S2*0.2+c(2),S3*0.2+c(3),'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7), hold on
surf(S1*0.1+y1(1),S2*0.1+y1(2),S3*0.1+y1(3),'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7)
surf(S1*0.07+y(1),S2*0.07+y(2),S3*0.07+y(3),'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7)

plot3([c(1) y1(1)],[c(2) y1(2)],[c(3) y1(3)],'LineWidth',6,'Color',[0.1 0.1 0.8])
plot3([y1(1) y(1)],[y1(2) y(2)],[y1(3) y(3)],'LineWidth',5,'Color',[0.1 0.1 0.8])
plot3([y1(1) y(1)],[y1(2) y(2)],[y1(3)-0.1 y(3)-0.1],'LineWidth',4,'Color',[0.1 0.1 0.8])

MDP_DEM_Movements_Hand([pi,pi/2 + atan((y(3)-y1(3))/(sqrt((y(1) - y1(1))^2+(y(2) - y1(2))^2))),x(1)],0.25,y)
axis equal

function y = MDP_DEM_MD_dthetadt(v,x,a,b,c)

y1 = c + [a*cos(x(2))*cos(x(1));
    a*cos(x(2))*sin(x(1));
    a*sin(x(2))];

y2 = y1 + [b*cos(x(3)+x(2))*cos(x(1));
    b*cos(x(3)+x(2))*sin(x(1));
    b*sin(x(3)+x(2))];

y(1) = [v - y2]'*[-((a*cos(x(2)) + b*cos(x(3)+x(2))*sin(x(1))));
    ((a*cos(x(2)) + b*cos(x(3)+x(2))*cos(x(1))));
    0];

y(2) = [v - y2]'*[-((a*sin(x(2)) + b*sin(x(3) + x(2)))*cos(x(1)));
    -((a*sin(x(2)) + b*sin(x(3) + x(2))*sin(x(1))));
    (a*cos(x(2)) +(b*cos(x(2) + x(3))))];

y(3) = [v - y2]'*[-(a*sin(x(3) + x(2))*cos(x(1)));
    -(a*sin(x(3)+ x(2))*sin(x(1)));
    (b*cos(x(2) + x(3)))];

y = 0.05*y'; % 0.05

function MDP_DEM_Movements_Hand(a,s,l)
% angle    - a
% scale    - s
% location - l

% Carpals (as a single sphere)
%--------------------------------------------------------------------------
[X,Y,Z] = sphere;
Rx = makehgtform('xrotate',a(1));
Ry = makehgtform('yrotate',a(2));
Rz = makehgtform('zrotate',a(3));
T  = makehgtform('translate',l);
S  = makehgtform('scale',s);
q = hgtransform;

LW = 5;

% 2nd finger
%--------------------------------------------------------------------------

% MCP
%--------------------------------------------------------------------------
plot3([0 0],[0 0],[0 1.1],'b','LineWidth',LW,'Parent',q)
surf(0.1*X, 0.1*Y ,0.1*Z + 1.1,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)
% PIP
%--------------------------------------------------------------------------
plot3([0 0],[0 0.4],[1.1 1.4],'b','LineWidth',LW,'Parent',q)
surf(0.1*X, 0.1*Y + 0.4 ,0.1*Z + 1.4,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)
% MIP
%--------------------------------------------------------------------------
plot3([0 0],[0.4 0.7],[1.4 1.5],'b','LineWidth',LW,'Parent',q)
surf(0.1*X, 0.1*Y + 0.7 ,0.1*Z + 1.5,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)
% DIP
%--------------------------------------------------------------------------
plot3([0 0],[0.7 1],[1.5 1.6],'b','LineWidth',LW,'Parent',q)
surf(0.1*X, 0.1*Y + 1 ,0.1*Z + 1.6,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)

% 3rd finger
%--------------------------------------------------------------------------

% MCP
%--------------------------------------------------------------------------
plot3([-0.1 -0.2],[0 0],[0 1.05],'b','LineWidth',LW,'Parent',q)
surf(0.1*X - 0.2, 0.1*Y ,0.1*Z + 1.05,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)
% PIP
%--------------------------------------------------------------------------
plot3([-0.2 -0.2],[0 0.3],[1.05 1.4],'b','LineWidth',LW,'Parent',q)
surf(0.1*X - 0.2, 0.1*Y + 0.3 ,0.1*Z + 1.4,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)
% MIP
%--------------------------------------------------------------------------
plot3([-0.2 -0.2],[0.3 0.6],[1.4 1.5],'b','LineWidth',LW,'Parent',q)
surf(0.1*X - 0.2, 0.1*Y + 0.6 ,0.1*Z + 1.5,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)
% DIP
%--------------------------------------------------------------------------
plot3([-0.2 -0.2],[0.6 0.9],[1.5 1.6],'b','LineWidth',LW,'Parent',q)
surf(0.1*X - 0.2, 0.1*Y + 0.9 ,0.1*Z + 1.6,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)

% 4th finger
%--------------------------------------------------------------------------

% MCP
%--------------------------------------------------------------------------
plot3([-0.2 -0.4],[0 0],[0 1],'b','LineWidth',LW,'Parent',q)
surf(0.1*X - 0.4, 0.1*Y ,0.1*Z + 1,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)
% PIP
%--------------------------------------------------------------------------
plot3([-0.4 -0.4],[0 0.2],[1 1.4],'b','LineWidth',LW,'Parent',q)
surf(0.1*X - 0.4, 0.1*Y + 0.2 ,0.1*Z + 1.4,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)
% MIP
%--------------------------------------------------------------------------
plot3([-0.4 -0.4],[0.2 0.5],[1.4 1.5],'b','LineWidth',LW,'Parent',q)
surf(0.1*X - 0.4, 0.1*Y + 0.5 ,0.1*Z + 1.5,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)
% DIP
%--------------------------------------------------------------------------
plot3([-0.4 -0.4],[0.5 0.8],[1.5 1.6],'b','LineWidth',LW,'Parent',q)
surf(0.1*X - 0.4, 0.1*Y + 0.8 ,0.1*Z + 1.6,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)

% 1st finger
%--------------------------------------------------------------------------

% MCP
%--------------------------------------------------------------------------
plot3([0.1 0.2],[0 0],[0 1],'b','LineWidth',LW,'Parent',q)
surf(0.1*X + 0.2, 0.1*Y ,0.1*Z + 1,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)
% PIP
%--------------------------------------------------------------------------
plot3([0.2 0.2],[0 0.3],[1 1.4],'b','LineWidth',LW,'Parent',q)
surf(0.1*X + 0.2, 0.1*Y + 0.3 ,0.1*Z + 1.4,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)
% MIP
%--------------------------------------------------------------------------
plot3([0.2 0.2],[0.3 0.6],[1.4 1.5],'b','LineWidth',LW,'Parent',q)
surf(0.1*X + 0.2, 0.1*Y + 0.6 ,0.1*Z + 1.5,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)
% DIP
%--------------------------------------------------------------------------
plot3([0.2 0.2],[0.6 0.9],[1.5 1.6],'b','LineWidth',LW,'Parent',q)
surf(0.1*X + 0.2, 0.1*Y + 0.9 ,0.1*Z + 1.6,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)

% Thumb
%--------------------------------------------------------------------------

% MCP
%--------------------------------------------------------------------------
plot3([0.2 0.5],[0 0.05],[0 0.4],'b','LineWidth',LW,'Parent',q)
surf(0.1*X + 0.5, 0.1*Y +0.05,0.1*Z + 0.4,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)
% PIP
%--------------------------------------------------------------------------
plot3([0.5 0.5],[0.05 0.45],[0.4 0.6],'b','LineWidth',LW,'Parent',q)
surf(0.1*X + 0.5, 0.1*Y +0.45,0.1*Z + 0.6,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)
% DIP
%--------------------------------------------------------------------------
plot3([0.5 0.45],[0.45 0.6],[0.6 0.62],'b','LineWidth',LW,'Parent',q)
surf(0.1*X + 0.45, 0.1*Y +0.6,0.1*Z + 0.62,'Facecolor','b','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7,'Parent',q)

set(q,'Matrix',T*Rz*Rx*Ry*S);

function [Q, R, X, Y, Z] = MDP_DEM_Movement_Disorders_Animate(MDP,L)
opengl('software')
Q = [];
R = [];

for m = 1:length(MDP.mdp)
    for j = 1:length(MDP.mdp(m).dem)
        for i = 2:size(MDP.mdp(m).dem(j).pU.x{1},2)
            MDP_DEM_MD_plot_arm(full(MDP.mdp(m).dem(j).pU.x{1}(1:3,i)),MDP.mdp(m).dem(j).M(1).pE.b1,MDP.mdp(m).dem(j).M(1).pE.b2,[-1;-1;-1]), hold on
            [X,Y,Z] = sphere;
            surf(0.1*X + MDP.mdp(m).dem(j).qU.v{2}(1,i),0.1*Y + MDP.mdp(m).dem(j).qU.v{2}(2,i),0.1*Z + MDP.mdp(m).dem(j).qU.v{2}(3,i),'Facecolor','r','EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7)
            for k = 1:3
                surf(0.2*X + L{(k-1)*4+1}(1),0.2*Y + L{(k-1)*4+1}(2),0.2*Z + L{(k-1)*4+1}(3),'Facecolor',[1 1 1]-MDP.mdp(m).dem(j).Y(9+k,i),'EdgeColor','none','FaceLighting','gouraud','AmbientStrength',0.7)
            end
            axis equal
            light
            drawnow
            hold off

            Q(:,end+1) = MDP_DEM_MD_transform(full(MDP.mdp(m).dem(j).pU.x{1}(1:3,i)),MDP.mdp(m).dem(j).M(1).pE.b1,MDP.mdp(m).dem(j).M(1).pE.b2,[-1;-1;-1]);
        end
        R = [R full(MDP.mdp(m).dem(j).pU.x{1}(1:3,:))];
    end
end


text:
function [costs] = initCosts()
  costs.total = 0;
  costs.word = 0;
end

author: lmthang


text:
function ecc = getEccFromRpAndSma(Sma, Rp)
%getEccFromRpAndSma Summary of this function goes here
%   Detailed explanation goes here
    ecc = 1 - Rp/Sma;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ w, x ] = line_unit_o02 ( )

%*****************************************************************************80
%
%% LINE_UNIT_O02 returns a 2 point quadrature rule for the unit line.
%
%  Discussion:
%
%    The integration region is:
%
%    - 1.0 <= X <= 1.0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carlos Felippa,
%    A compendium of FEM integration formulas for symbolic work,
%    Engineering Computation,
%    Volume 21, Number 8, 2004, pages 867-890.
%
%  Parameters:
%
%    Output, real W(2), the weights.
%
%    Output, real X(2), the abscissas.
%
  w(1:2,1) = [ ...
    1.0000000000000000000, ...
    1.0000000000000000000 ];

  x(1:2,1) = [ ...
    -0.57735026918962576451, ...
     0.57735026918962576451 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% ASINH  Inverse hyperbolic sine.
%    ASINH(X) is the inverse hyperbolic sine of the elements of X.
%
%    See also SINH.
%
%    Reference page in Doc Center
%       doc asinh
%
%    Other functions named asinh
%
%       codistributed/asinh    gpuArray/asinh    sym/asinh    ts/asinh
%

author: jmaih


text:
function [cq mcq Xc]=getmeanquality(X,c,linkagetype)

    %point estimate cluster solution
    Y = pdist1(X);

    if exist('linkage.m', 'file') && exist('cluster.m', 'file')
        Z = linkage(Y,linkagetype);
        Xc = cluster(Z,c);
    else
        Z = linkage_t(Y,linkagetype);
        Xc = cluster_t(Z,c);
    end

    %convert X to X, binary input to doquality
    Xcx=makebinary(Xc)';

    %point estimate quality of each element;
    [cq center]=clustquality(Xcx,X);      %cq is cluster quality of each element

    %mean cq
    if all(cq == 1)
        mcq = 0;
    else
        mcq = mean( cq(cq~=1) );    %don't include single element clusters in your quality estimate
    end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [c,u,l]=bsds(bench,models,B,w,type,boot)
% Calculate Whites and Hansens p-vals for out-performance using unmodified data or studentized
% residuals,  the latter often providing better power, particularly when the losses functions are
% heteroskedastic
%
% USAGE:
%   [C] = bsds(BENCH,MODELS,B,W)
%   [C,U,L] = bsds(BENCH,MODELS,B,W,TYPE,BOOT)
%
% INPUTS:
%   BENCH  - Losses from the benchmark model
%   MODELS - Losses from each of the models used for comparison
%   B      - Number of Bootstrap replications
%   W      - Desired block length
%   TYPE   - String, either 'STANDARD' or 'STUDENTIZED'.  'STUDENTIZED' is the default, and
%              generally leads to better power.
%   BOOT   - [OPTIONAL] 'STATIONARY' or 'BLOCK'.  Stationary is used as the default.
%
% OUTPUTS:
%   C      - Consistent P-val(Hansen)
%   U      - Upper P-val(White) (Original RC P-vals)
%   L      - Lower P-val(Hansen)
%
% COMMENTS:
%   This version of the BSDS operates on quantities that should be 'bads', such as losses.  The null
%   hypothesis is that the average performance of  the benchmark is as small as the minimum average
%   performance across the models.  The alternative is that the minimum average loss across the
%   models is smaller than the the average performance of the benchmark.
%
%   If the quantities of interest are 'goods', such as returns, call bsds with
%   -1*BENCH and -1*MODELS
%
% EXAMPLES:
%   Standard Reality Check with 1000 bootstrap replications and a window size of 12
%       bench = randn(1000,1).^2;
%       models = randn(1000,100).^2;
%       [c,realityCheckPval] = bsds(bench, models, 1000, 12)
%   Standard Reality Check with 1000 bootstrap replications, a window size of 12 and a circular
%   block bootstrap
%       [c,realityCheckPval] = bsds(bench, models, 1000, 12, 'BLOCK')
%   Hansen's P-values
%       SPAPval = bsds(bench, models, 1000, 12)
%   Both Pvals on "goods"
%       bench = .01 + randn(1000,1);
%       models = randn(1000,100);
%       [SPAPval,realityCheckPval] = bsds(-bench, -models, 1000, 12)
%
% See also MCS

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 4/1/2007

%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<4 || nargin>6
    error('4 to inputs required')
end
if nargin == 4
    boot = 'STATIONARY';
    type = 'STUDENTIZED';
elseif nargin == 5
    boot = 'STATIONARY';
end
if isempty(type)
    type = 'STUDENTIZED';
end
if strcmpi(type,'STUDENTIZED')
    isStudentized = true;
else
    isStudentized = false;
end
% Get the length of the data
[tb,kb]=size(bench);
if kb>1
    error('BENCH must be a column vector')
end
if tb<2
    error('BENCH must have at least 2 observations.')
end
[t,k]=size(models);
if t~=tb
    error('BENCH and MODELS must have the same number of observations.')
end
if ~isscalar(B) || B<1 || floor(B)~=B
    error('B must be a positive scalar integer')
end
if ~isscalar(w) || w<1 || floor(w)~=w
    error('W must be a positive scalar integer')
end
boot = upper(boot);
if ~ismember(boot,{'STATIONARY','BLOCK'})
    error('BOOT must be either ''STATIONARY'' or ''BLOCK''.')
end
%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(boot,'BLOCK')
    bsdata=block_bootstrap((1:t)',B,w);
else
    bsdata=stationary_bootstrap((1:t)',B,w);
end

%OK now we have the bootstraps, what to do with them?
diffs=models-repmat(bench,1,k);

% First compute the boostarap sample averages, db*
% Second compute the variance estimate, omegak
% First the weghts
q=1/w;
i=1:t-1;
kappa=((t-i)./t).*(1-q).^i+i./t.*(1-q).^(t-i);
% Next compute the variances
vars=zeros(k,1)';
for i=1:k
    workdata = diffs(:,i)-mean(diffs(:,i));
    vars(i)= workdata'*workdata/t;
    for j=1:t-1
        vars(i) = vars(i) + 2*kappa(j)*(workdata(1:t-j)'*workdata(j+1:t))/t;
    end
end

% Aold is the original method to compute the truncation point
Aold=1/4*t^(0.25)*sqrt(vars/t);
mean(Aold);
% A new used the log(log(t)) rule
Anew = sqrt((vars/t)*2*log(log(t)));

% Only recenter if the average is reasonably small or the model is better
% (in which case mean(diffs) is negative).  If it is unreasonably large set
% the mean adjustment to 0
gc=mean(diffs).*(mean(diffs)<Anew);

% The lower assumes that every loss function that is worse than BM is
% unimportant for the asymptotic distribution, hence if its mean is
% less than 0, g=0.  This is different from the consistent where the
% threshold was it had to be greater than -A(i)
gl=min(0,mean(diffs));

%Then the upper, which assumes all models used are reasonably close to
%the benchmark that they coudl be better
gu=mean(diffs);

% Perf will hold the boostrapped statistics for B iterations
perfc=zeros(B,k);
perfl=zeros(B,k);
perfu=zeros(B,k);
if isStudentized
    stdDev = sqrt(vars);
else
    stdDev = ones(1,k);
end

for i=1:k
    workdata=diffs(:,i);
    % the i'th column of perf holds the B bootstrapped statistics
    mworkdata=mean(workdata(bsdata));
    perfc(:,i)=(mworkdata-gc(i))'/stdDev(i);
    perfl(:,i)=(mworkdata-gl(i))'/stdDev(i);
    perfu(:,i)=(mworkdata-gu(i))'/stdDev(i);
end
% Compute the test statistic
stat = min(mean(diffs)./stdDev);
% Compute the min in each row
perfc=min(perfc,[],2);
perfc=min(perfc,0);
perfl=min(perfl,[],2);
perfl=min(perfl,0);
perfu=min(perfu,[],2);
perfu=min(perfu,0);
% Count the number of time the min is below the statistic
c=mean(perfc<stat);
l=mean(perfl<stat);
u=mean(perfu<stat);

author: bashtage

repo: mfe-toolbox


text:
 function jf_slicer(data, varargin)
%function jf_slicer(data, [options])
%|
%| slice 3d data interactively (along 3rd dimension).
%| Uses scroll wheel to sweep through all slices.
%|
%| todo: hsv!
%|
%| in
%|	data	[nx ny nz]
%|
%| options
%|	clim	[1 2]		clim arg to im()
%|	iz	[1]		initial slice (default: nz/2+1)
%|
%| Jeff Fessler, University of Michigan

if ~nargin, ir_usage, end
if streq(data, 'test'), jf_slicer_test, return, end

arg.clim = [];
arg.iz = [];
arg = vararg_pair(arg, varargin);
if ~isreal(data)
	printm 'warning: taking abs of complex data'
	data = abs(data);
end

if isempty(arg.clim)
	arg.clim = double(minmax(data)');
end

nz = size(data, 3);

if isempty(arg.iz)
	iz = ceil((nz+1)/2);
else
	iz = arg.iz;
end

clamp = @(iz) max(min(iz, nz), 1);

stuff.data = data;
stuff.arg = arg;

%jf_slicer_call(iz, stuff)
%drawnow
jf_slicer_show

set(gcf, 'WindowScrollWheelFcn', @jf_slicer_scroll)
set(gcf, 'WindowKeyPressFcn', @jf_slicer_key)

%h = jf_add_slider('callback', @jf_slicer_call, 'data', {stuff}, ...
%	'min', 1, 'max', nz, 'sliderstep', [1 1]/(nz-1), 'value', iz);

function jf_slicer_scroll(src, event)
	iz = iz + event.VerticalScrollCount;
	iz = clamp(iz);
	jf_slicer_show
end % jf_slicer_scroll

function jf_slicer_key(src, event)
	key = get(gcf, 'CurrentCharacter');
	key = 0 + key;
	switch key
	case {28, 31, 45} % left arrow, down arrow, -
		iz = iz - 1;
	case {29, 30, 43} % right arrow, up arrow, +
		iz = iz + 1;
	otherwise
		if 49 <= key && key <= 57
			iz = key - 48; % 1-9
		end
	end
	iz = clamp(iz);
	pr iz
	jf_slicer_show
end % jf_slicer_key

function jf_slicer_show
	im(data(:,:,iz), arg.clim), cbar
	xlabelf('%d / %d', iz, nz)
end % jf_slicer_show

end % jf_slicer

function jf_slicer_call(iz, stuff)
persistent iz_save
if isempty(iz_save)
	iz_save = -1;
end
iz = round(iz);
if iz ~= iz_save
	iz_save = iz;
	arg = stuff.arg;
	im(stuff.data(:,:,iz), arg.clim), cbar
end

end % jf_slicer_call

function jf_slicer_test
data = reshape(1:7*8*9, 7, 8, 9);
im plc 1 2
im subplot 2
if im
	jf_slicer(data, 'clim', [0 500])
end
end % jf_slicer_test

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function basis_mn_t4_test ( )

%*****************************************************************************80
%
%% BASIS_MN_T4_TEST verifies BASIS_MN_T4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  node_num = 4;

  t = [ ...
    2.0, 0.0; ...
    4.0, 2.0; ...
    0.0, 4.0; ...
    2.0, 2.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BASIS_MN_T4_TEST:\n' );
  fprintf ( 1, '  Verify basis functions for element T4.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes = %d\n', node_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Physical Nodes:\n' );
  fprintf ( 1, '\n' );
  for j = 1 : node_num
    fprintf ( 1, '  %8d  %7f  %7f\n', j, t(1:2,j) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The basis function values at basis nodes\n' );
  fprintf ( 1, '  should form the identity matrix.\n' );
  fprintf ( 1, '\n' );

  [ phi, dphidx, dphidy ] = basis_mn_t4 ( t, node_num, t );

  for i = 1 : node_num
    for j = 1 : node_num
      fprintf ( 1, '  %7f', phi(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The X and Y derivatives should sum to 0.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      dPhidX sum    dPhidY sum\n' );
  fprintf ( 1, '\n' );

  for j = 1 : node_num

    sum_x = sum ( dphidx(1:node_num,j) );
    sum_y = sum ( dphidy(1:node_num,j) );
    fprintf ( 1, '  %14f  %14f\n', sum_x, sum_y );

  end

  return
end


text:
function combo_test37 ( )

%*****************************************************************************80
%
%% COMBO_TEST37 tests STIRLING_NUMBERS1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
  maxm = 6;
  maxn = 6;
  offset = 1;

  fprintf ( 1, ' \n' );
  fprintf ( 1, 'COMBO_TEST37\n' );
  fprintf ( 1, '  STIRLING_NUMBERS1 computes a table of Stirling\n' );
  fprintf ( 1, '  numbers of the first kind.\n' );

  s = stirling_numbers1 ( maxm, maxn );

  fprintf ( 1, ' \n' );
  fprintf ( 1, '    I S(I,0) S(I,1) S(I,2) S(I,3) S(I,4) S(I,5)\n' );
  fprintf ( 1, ' \n' );

  for i = 0 : maxm
    fprintf ( 1, '%5d', i );
    for j = 0 : maxn
      fprintf ( 1, '%5d', s(i+offset,j+offset) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function findedge(node,edge,range,varargin)
%% FINDEDGE highlights edges
%
%    FINDEDGE(node,edge,range) finds all elements whose indices are in the
%    range array by displaying these elements in yellow.
%
%    FINDEDGE(node,edge) finds all elements.
%
%    FINDEDGE(node,edge,'noindex') skip the display of indices.
%
%    FINDEDGE(node,edge,range,'param','value','param','value'...) allows
%    additional patch param/value pairs to highlight the elements.
%
% Example:
%     node = [1,0; 1,1; 0,1; -1,1; -1,0; -1,-1; 0,-1; 0,0];
%     elem = [1,2,8; 3,8,2; 8,3,5; 4,5,3; 7,8,6; 5,6,8];
%     subplot(1,2,1);
%     showmesh(node,elem);
%     T = auxstructure(elem);
%     findedge(node,T.edge,5,'index','color','r','MarkerSize',24);
%     subplot(1,2,2);
%     showmesh(node,elem);
%     findedge(node,T.edge);
%
%   See also findelem3, findnode3, findedge.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

% set up range
if (nargin<=2) || (isempty(range))
    range = (1:size(edge,1))';
end
if islogical(range)
    range = find(range);
end
if size(range,2)>size(range,1)
    range = range';
end
% draw edges in red
midEdge = (node(edge(range,1),:)+node(edge(range,2),:))/2;
hold on
if length(range) < size(edge,1)
    h = line([node(edge(range,1),1)'; node(edge(range,2),1)'],...
             [node(edge(range,1),2)'; node(edge(range,2),2)'],...
             'LineWidth',2,'Color','r');
 else
    h = plot(midEdge(:,1),midEdge(:,2),'r.','MarkerSize', 18);
end
% else
%     h = plot(midEdge(:,1),midEdge(:,2),'r.','MarkerSize', 18);
% end
if nargin > 3
    if strcmp(varargin{1},'noindex') || strcmp(varargin{1},'index') && length(varargin)>1
        startidx = 2;
    else
        startidx = 1;
    end
    set(h,varargin{startidx:end});
end
if (nargin <=3) || ~(strcmp(varargin{1},'noindex'))
    text(midEdge(:,1)+0.025,midEdge(:,2)-0.025,int2str(range), ...
         'FontSize',12,'FontWeight','bold','Color','k');
end

repo: Aero-matlab


text:
function include_params = include_codes(segm_params, filepath_params)
%INCLUDE_CODES is used for putting codes in path according to the parameter
% settings
%
% @authors:     Ahmad Humayun
% @contact:     ahumayun@cc.gatech.edu
% @affiliation: Georgia Institute of Technology
% @date:        Fall 2013 - Summer 2014

    extern_dir = filepath_params.extern_src_dir;
    code_root_dir = filepath_params.code_root_dir;
    extern_codes_dir = filepath_params.extern_codes_dir;
    utils_dir = filepath_params.utils_src_dir;

    include_params.vlfeat_dir = fullfile(extern_dir, 'toolboxes', ...
                                         'vlfeat');
    addpath(fullfile(include_params.vlfeat_dir, 'toolbox'));
%    vl_setup;

    include_params.pb_dir = fullfile(extern_dir, ...
                                     'segmentation', 'segbench');

    % includes from piotr's toolbox for sketchtokens
    piotr_toolbox_dir = fullfile(extern_dir, 'toolboxes', ...
                                 'piotr_toolbox');

    switch segm_params.boundaries_method
        case 'Gb'
            addpath(fullfile(include_params.pb_dir, 'lib', 'matlab'));
            include_params.gb_dir = fullfile(extern_dir, ...
                                             'segmentation', ...
                'boundaries--leordeanu_ECCV_2012_gb');
            addpath(include_params.gb_dir);
            addpath(fullfile(extern_codes_dir, 'extra_gb_code'));
        case 'GPb'
            include_params.globalpb_dir = fullfile(extern_dir, ...
                                                   'segmentation', ...
                'boundaries+segments--arbelaez_PAMI_2010_bsr');
            addpath(include_params.globalpb_dir);
            addpath(fullfile(include_params.globalpb_dir, 'lib'));
        case 'Pb'
            addpath(fullfile(include_params.pb_dir, 'lib', 'matlab'));
        case 'SketchTokens'
            include_params.sketchtokens_dir = fullfile(extern_dir, ...
                                                       'segmentation', ...
                'boundaries--lim_CVPR_2013_sketchtokens');
            addpath(include_params.sketchtokens_dir);

            addpath(fullfile(piotr_toolbox_dir, 'channels'));
        case 'StructEdges'
            include_params.structedges_dir = fullfile(extern_dir, ...
                                                      'segmentation', ...
                'boundaries--dollar_ICCV_2013_structedges');
            addpath(include_params.structedges_dir);

            addpath(fullfile(piotr_toolbox_dir, 'channels'));
        otherwise
    end

    % add parametric min-cut code to the path
    pmc_maxflow_method = segm_params.pmc_maxflow_method;
    if strcmp(pmc_maxflow_method, 'hochbaum')
        addpath(fullfile(extern_codes_dir, 'para_pseudoflow'));
    elseif strcmp(pmc_maxflow_method, 'kolmogorov')
        addpath(fullfile(extern_dir, ...
                'optimization/maxflow--boykov_PAMI_2004_maxflow/mex_maxflow'));
    elseif ~isempty(strfind(pmc_maxflow_method, 'nodynamic')) || ...
           ~isempty(strfind(pmc_maxflow_method, 'kohli')) || ...
           ~isempty(strfind(pmc_maxflow_method, 'multiseed'))
       addpath(fullfile(code_root_dir, 'boykov_maxflow'));
    else
    end

    addpath(filepath_params.fuxin_lib_dir);
    addpath(fullfile(filepath_params.fuxin_lib_dir, 'boosting'));

    % if using the seed sampler, include stuff from fuxin's library
%    if any(strcmpi(segm_params.graph_seed_gen_method, 'sp_seed_sampling'))
%        addpath(fullfile(filepath_params.fuxin_lib_dir, 'myqueue_1.1'));
%    end

   % only include vgg if planning to use FelzHutten superpixels
   if any(cellfun(@(x) any(strcmp(x, 'felzenszwalb_seeds_caller')), ...
          segm_params.graph_seed_params))
       addpath(fullfile(extern_dir, 'segmentation', 'imrender', 'vgg'));
   end

   addpath(fullfile(extern_dir, 'segmentation', ...
                    'stein_boundaryprocessing'));

   addpath(utils_dir);
   addpath(extern_dir);

%     % get export_fig in the path
%     addpath(fullfile(extern_dir, 'utils', 'export_fig'));
end

author: Cloud-CV

repo: object-proposals


text:
classdef RMMEDA_F8 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP for testing RM-MEDA

%------------------------------- Reference --------------------------------
% Q. Zhang, A. Zhou, and Y. Jin, RM-MEDA: A regularity model-based
% multiobjective estimation of distribution algorithm, IEEE Transactions on
% Evolutionary Computation, 2008, 12(1): 41-63.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 3;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            g = sum((X(:,3:end).^2-repmat(X(:,1),1,size(X,2)-2)).^2,2);
            PopObj(:,1) = cos(pi/2*X(:,1)).*cos(pi/2*X(:,2)).*(1+g);
            PopObj(:,2) = cos(pi/2*X(:,1)).*sin(pi/2*X(:,2)).*(1+g);
            PopObj(:,3) = sin(pi/2*X(:,1)).*(1+g);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,3);
            R = R./repmat(sqrt(sum(R.^2,2)),1,3);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            a = linspace(0,pi/2,10)';
            R = {sin(a)*cos(a'),sin(a)*sin(a'),cos(a)*ones(size(a'))};
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function aa = fitnesstra2(pop,para)
[pops,numvar]=size(pop);
for i=1:pops
    time=pop(i,5)+pop(i,6);
    kk=trajt(para,pop(i,:));
    tt=torque(kk);
    ft=ftorque(tt);
    poo=kk(1:2,:);
    fq=sum(sum(abs(diff(poo'))));
    pos=forkin(poo);xx=pos(1,:);yy=pos(2,:);
    x=diff(xx);
    y=diff(yy);
    dis=sqrt(x.^2+y.^2);
    fdis=sum(dis);
    a1=2;
    a2=1.8;
    a3=2;
    a4=1;
    f=a1*ft+a2*fq+a3*fdis+a4*time;
    aa(i)=1/f;
end

author: Sable

repo: mcbench-benchmarks


text:
clc; close all; clear all;
%%%%%%%%% this directory should be properly changed.
%%%%%% make sure that RANSAC and GP_toolbox had been added the path.
% addpath('GPSegmentation\Ransac\');
% addpath('iGPR/gpml-matlab-v4.0-2016-10-19');
% run('iGPR/gpml-matlab-v4.0-2016-10-19/setpath.m');
% load data_straightRoad.mat
load data_intersection.mat
cloud = pcdownsample(pointCloud(data'), 'gridAverage', 0.1);
data = cloud.Location';
Dist = sqrt(sum(data(1:2, :).^2));
data = data(:, Dist >= 3.0 & Dist <= 60.0);
params = [];
params.GapThr = [0.10 1.0];
params.RadArray = 0.0 : 2.0 : 80.0;
params.AngRes   = deg2rad(5.0);
params.GroundH  = -1.8;
params.IS_SHOW = 1;
tic
[GrdIdx, ObsIdx, UnkownIdx, GrdPts] = GPSegFun(data, params);
toc
%% Run DBSCAN Clustering Algorithm
obsData = data(:, ObsIdx);
X = obsData(1:2, :)';
epsilon = 1.0; % the unit is meter.
MinPts = 10;
tic
IDX = DBSCAN(X,epsilon,MinPts);
toc
%% Plot Results
figure;
hold on;
grid on;
axis equal;
xlabel('X/m');
ylabel('Y/m');
% PlotClusterinResult(X, IDX);
Colors = hsv(20);
nDim = 3;
Legends = {};
%%%%% IDX == 0 denotes noise.
K=max(IDX);
for i=1:K
    Tmp = obsData(:, IDX==i);
    if length(Tmp) <= 20 %%%% discarding small obstacle
        continue;
    end
    minRange = min(Tmp');
    maxRange = max(Tmp');
    xRange = [minRange(1) maxRange(1)];
    yRange = [minRange(2) maxRange(2)];
    zRange = [minRange(3) maxRange(3)];
    boxPts = [xRange([1 1 2 2 1]);
              yRange([1 2 2 1 1]);
              zRange([1 1 1 1 1])];
    id = mod(i, length(Colors));
    if id == 0
        id = length(Colors);
    end
    Color = Colors(id,:);
    plot3(boxPts(1, :), boxPts(2, :), boxPts(3, :), 'color', Color, 'linestyle', '-', ...
        'linewidth', 2, 'marker', 's');
    if ~isempty(Tmp)
        str = sprintf('Id=%02d, Len=%02d', i-1, length(Tmp));
        pt = Tmp(:, 1);
        text(pt(1), pt(2), pt(3), str, 'FontSize', 12, 'color', 'k' );
        pcshow(Tmp', Color);
    end
end
title(['DBSCAN Clustering (\epsilon = ' num2str(epsilon) ', MinPts = ' num2str(MinPts) ')']);

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR


text:
function params = crazyflie()
% crazyflie: physical parameters for the Crazyflie 2.0
%
% 2016 Bernd Pfrommer
%
% This function creates a struct with the basic parameters for the
% Crazyflie 2.0 quad rotor (without camera, but with about 5 vicon
% markers)
%
% Model assumptions based on physical measurements:
%
% motor + mount + vicon marker = mass point of 3g
% arm length of mass point: 0.046m from center
% battery pack + main board are combined into cuboid (mass 18g) of
% dimensions:
%
%   width  = 0.03m
%   depth  = 0.03m
%   height = 0.012m
%

m = 0.030;  % weight (in kg) with 5 vicon markers (each is about 0.25g)
g = 9.81;   % gravitational constant
I = [1.43e-5,   0,          0; % inertial tensor in m^2 kg
     0,         1.43e-5,    0;
     0,         0,          2.89e-5];
L = 0.046; % arm length in m

params.mass = m;
params.I    = I;
params.invI = inv(I);
params.grav = g;
params.arm_length = L;

params.maxangle = 40*pi/180; % you can specify the maximum commanded angle here
params.maxF     = 2.5*m*g;   % left these untouched from the nano plus
params.minF     = 0.05*m*g;  % left these untouched from the nano plus

% You can add any fields you want in params
% for example you can add your controller gains by
% params.k = 0, and they will be passed into controller.m

end

author: yrlu

repo: quadrotor

save_path: github-repos/MATLAB/yrlu-quadrotor


text:
% ====
% disp
% ====
%
% Description :
% -------------
%
% disp function for class inf_engine.
%
% Syntax :
% --------
%
% [] = disp( obj )
%
% Input(s) :
% ----------
%
% obj - class inf_engine
% An instance of the class inf_engine
%
% Output(s) :
% -----------
%
% Example(s) :
% ------------
%
% disp( obj );
%
% Reference(s) :
% --------------
%
% See also :
% ----------
%
% display
function [] = disp( obj )

disp( struct( obj ) );

% End of function

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function im1 = cropImage(im,center,scale)

w = 200*scale;
h = w;
x = center(1) - w/2;
y = center(2) - h/2;
bbox = [x,y,w,h];
padsize = round(bbox(3:4));
im = padarray(im,padsize,0);
im1 = imcrop(im,[bbox(1:2)+padsize,bbox(3:4)]);
im1 = imresize(im1,[200,200]);

author: geopavlakos


text:
function [id, gc_x, gc_y, gc_z] = load_gc_traj_data(gcTrajFile)
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

    %% Import data from text file.
    % Script for importing data from the following text file:

    %% Initialize variables.
    delimiter = ' ';
    startRow = 1;

    %% Format string for each line of text:
    %   column1: double (%d)
    %	column2: double (%f)
    %   column3: double (%f)
    %	column4: double (%f)
    %   column5: double (%f)
    %	column6: double (%f)
    %   column7: double (%f)
    % For more information, see the TEXTSCAN documentation.
    formatSpec = '%d%f%f%f%[^\n\r]';

    %% Open the text file.
    fileID = fopen(gcTrajFile,'r');

    %% Read columns of data according to format string.
    % This call is based on the structure of the file used to generate this
    % code. If an error occurs for a different file, try regenerating the code
    % from the Import Tool.
    dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'EmptyValue' ,NaN,'HeaderLines' ,startRow-1, 'ReturnOnError', false);

    %% Close the text file.
    fclose(fileID);

    %% Allocate imported array to column variable names
    id = dataArray{:, 1};
    gc_x = dataArray{:, 2};
    gc_y = dataArray{:, 3};
    gc_z = dataArray{:, 4};

    %% Clear temporary variables
    clearvars gcTrajFile delimiter startRow formatSpec fileID dataArray ans;
end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function [sFile, ChannelMat] = in_fopen_spm(DataFile)
% IN_FOPEN_SPM: Open a SPM .mat/.dat file.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2017

%% ===== READ HEADER =====
% Check if SPM is in the path
if ~exist('file_array', 'file')
    error('SPM must be in the Matlab path to use this feature.');
end
% Get the two input file names: .mat and .dat
[fPath, fBase, fExt] = bst_fileparts(DataFile);
MatFile = bst_fullfile(fPath, [fBase, '.mat']);
DatFile = bst_fullfile(fPath, [fBase, '.dat']);
% If one is missing: error
if ~file_exist(MatFile) || ~file_exist(DatFile)
    error('The two files .dat and .mat must be available in the same folder.');
end
% Read header
sMat = load(MatFile, 'D');
D = sMat.D;
nChannels = length(D.channels);

% Warning: Supporting only files with one epoch
if (length(D.trials) > 1)
    error(['Only continuous SPM files are currently supported. Files with multiple trials cannot be imported.' 10 ...
           'Please contact us through the Brainstorm user forum to request this feature.']);
end

%% ===== FILL STRUCTURE =====
% Initialize returned file structure
sFile = db_template('sfile');
% Add information read from header
sFile.byteorder    = 'l';
sFile.filename     = MatFile;
sFile.format       = 'SPM-DAT';
sFile.prop.sfreq   = double(D.Fsample);
sFile.prop.nAvg    = 1;
sFile.prop.times   = (round(D.timeOnset(1) .* sFile.prop.sfreq) + [0, (D.Nsamples - 1)]) ./ sFile.prop.sfreq;
sFile.channelflag  = ones(nChannels,1);
sFile.device       = 'SPM';
sFile.comment      = fBase;
if isa(D.data, 'file_array')
    sFile.header.file_array = D.data;
elseif isstruct(D.data) && isfield(D.data, 'y') && isa(D.data.y, 'file_array')
    sFile.header.file_array = D.data.y;
else
    error('Could not find the file_array object in the SPM structure.');
end
sFile.header.nChannels  = nChannels;
sFile.header.gain       = ones(nChannels,1);

%% ===== CHANNEL FILE =====
% Initialize structure
ChannelMat = db_template('ChannelMat');
ChannelMat.Comment = [sFile.device ' channels'];
ChannelMat.Channel = repmat(db_template('channeldesc'), [1, nChannels]);
% Loop on each channel
for i = 1:nChannels
    if (D.channels(i).bad)
        sFile.channelflag(i) = -1;
    end
    if iscell(D.channels(i).label) && ~isempty(D.channels(i).label)
        ChannelMat.Channel(i).Name = D.channels(i).label{1};
    elseif ischar(D.channels(i).label) && ~isempty(D.channels(i).label)
        ChannelMat.Channel(i).Name = D.channels(i).label;
    else
        disp(sprintf('BST> Warning: No information avaible for channel #%d.', i));
    end
    % Convert channel types
    switch upper(D.channels(i).type)
        case 'MEGPLANAR',   ChannelMat.Channel(i).Type = 'MEG GRAD';
        case 'MEGMAG',      ChannelMat.Channel(i).Type = 'MEG MAG';
        otherwise,          ChannelMat.Channel(i).Type = upper(D.channels(i).type);
    end
    % Channel gains
    if isfield(D.channels(i), 'units') && ~isempty(D.channels(i).units)
        switch (D.channels(i).units)
            case 'fT',        sFile.header.gain(i) = 1e-15;
            case 'fT/mm',     sFile.header.gain(i) = 1e-12;
            case 'mV',        sFile.header.gain(i) = 1e-3;
            case {'uV','?V'}, sFile.header.gain(i) = 1e-6;
            otherwise,        sFile.header.gain(i) = 1;
        end
    end
end
% Read detailed information from .meg and .eeg fields
ChannelMat = read_fieldtrip_chaninfo(ChannelMat, D.sensors);

% Convert head points
if isfield(D, 'fiducials') && isfield(D.fiducials, 'pnt') && isfield(D.fiducials, 'label')
    for i = 1:length(D.fiducials.label)
        ChannelMat.HeadPoints.Label = D.fiducials.label(:)';
        ChannelMat.HeadPoints.Type  = repmat({'EXTRA'}, size(ChannelMat.HeadPoints.Label));
        ChannelMat.HeadPoints.Loc   = bst_units_ui(D.fiducials.unit, D.fiducials.pnt');
    end
end
% Convert fiducials
if isfield(D, 'fiducials') && isfield(D.fiducials, 'fid') && isfield(D.fiducials.fid, 'label') && isfield(D.fiducials.fid, 'pnt')
    for i = 1:length(D.fiducials.fid.label)
        switch lower(D.fiducials.fid.label{i})
            case {'nas', 'nasion', 'nz', 'fidnas', 'fidnz'}  % NASION
                ChannelMat.SCS.NAS = bst_units_ui(D.fiducials.unit, D.fiducials.fid.pnt(i,:));
            case {'lpa', 'pal', 'og', 'left', 'fidt9', 'leftear'} % LEFT EAR
                ChannelMat.SCS.LPA = bst_units_ui(D.fiducials.unit, D.fiducials.fid.pnt(i,:));
            case {'rpa', 'par', 'od', 'right', 'fidt10', 'rightear'} % RIGHT EAR
                ChannelMat.SCS.RPA = bst_units_ui(D.fiducials.unit, D.fiducials.fid.pnt(i,:));
        end
    end
    % Force re-alignment on the new set of NAS/LPA/RPA
    if ~isempty(ChannelMat.SCS) && ~isempty(ChannelMat.SCS.NAS) && ~isempty(ChannelMat.SCS.LPA) && ~isempty(ChannelMat.SCS.RPA)
        ChannelMat = channel_detect_type(ChannelMat, 1, 0);
    end
end

%% ===== EVENTS =====
if isfield(D, 'trials') && isfield(D.trials, 'events') && isfield(D.trials.events, 'type')
    % Get all the event types
    evtList = {D.trials.events.type};
    % Events list
    [uniqueEvt, iUnique] = unique(evtList);
    uniqueEvt = evtList(sort(iUnique));
    % Initialize events list
    sFile.events = repmat(db_template('event'), 1, length(uniqueEvt));
    % Build events list
    for iEvt = 1:length(uniqueEvt)
        % Find all the occurrences of this event
        iOcc = find(strcmpi(uniqueEvt{iEvt}, evtList));
        % Concatenate all times
        t = [D.trials.events(iOcc).time];
        % If there is a duration: add it
        occDuration = [D.trials.events(iOcc).duration];
        if (length(occDuration) == length(t))
            t(2,:) = t(1,:) + occDuration;
        end
        % Set event
        sFile.events(iEvt).label   = strtrim(uniqueEvt{iEvt});
        sFile.events(iEvt).times   = t;
        sFile.events(iEvt).epochs  = 1 + 0*t(1,:);
        sFile.events(iEvt).select  = 1;
        sFile.events(iEvt).channels = [];
        sFile.events(iEvt).notes    = [];
    end
end

end

author: brainstorm-tools


text:
 function [xs, ni] = eml_sps(x, Gt, yi, ci, ri, niter, curv)
%function [xs, ni] = eml_sps(x, Gt, yi, ci, ri, niter, curv)
%	One iteration of the ML-SPS algorithm for emission Poisson problem
%	(separable paraboloidal surrogates)
%	model: Y_i ~ Poisson(c_i [G x]_i + r_i)		WITH r_i > 0 REQUIRED!
%	in:
%		Gt	transpose of system matrix
%		see em_fbp.m for model, G, yi, ci, ri
%	out:
%		x [np,niter]	updated image vectors each iteration
%
%	Copyright Mar 2000, Jeff Fessler, The University of Michigan

if nargin < 3, ir_usage, end

[nb, na] = size(yi);

if ~isvar('ci') || isempty(ci)
	ci = ones(size(yi));
end
if ~isvar('ri') || isempty(ri)
	ri = zeros(size(yi));
end
if ~isvar('niter') || isempty(niter)
	niter = 1;
end

eml_check(yi, ci, ri)

	gi = sum(Gt)';	% g_i = sum_j g_ij

xs = zeros(numel(x), niter);
xs(:,1) = x(:);
%
%	loop over iterations
%
for ii=2:niter
	li = reshape(Gt' * x(:), size(yi));	% l=G*x "line integrals"
	yb = ci .* li + ri;			% predicted measurement means

	%	curvatures
	ni = eml_curvature(yi, ci, ri, li, yb, curv);

	dothi = ci .* (yi ./ yb - 1);

	x = x + (Gt * dothi(:)) ./ (Gt * (gi .* ni(:)));
	x = max(x,0);

	xs(:,ii) = x;
end


text:
%  Multivariate Squarefree factorization that is capable of identifying
%  squarefree factors and multiplicities accurately even if the polynomial
%  is perturbed
%
%   Syntax:  >>[F,res,cond] = SquarefreeFactor(p,tol);
%
%   Input   p -- a multivariate polynomial (as a coeff. matrix) to be factored
%         tol -- (optional) backward error tolerance
%
%  Output   F -- cell of size kx2. For each j = 1, 2, ..., k, F{j,1}, F{j,2}
%                    are a squarefree factor and its multiplicity respectively
%         res -- residual
%        cond -- condition number of the factorization

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function newnode=sms(node,face,iter,alpha,method)
%
% newnode=sms(node,face,iter,useralpha,method)
%
% simplified version of surface mesh smoothing
%
% author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
% date: 2009/10/21
%
% input:
%    node:  node coordinates of a surface mesh
%    face:  face element list of the surface mesh
%    iter:  smoothing iteration number
%    alpha: scaler, smoothing parameter, v(k+1)=alpha*v(k)+(1-alpha)*mean(neighbors)
%    method: same as in smoothsurf, default is 'laplacianhc'
%
% output:
%    newnode: output, the smoothed node coordinates
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(nargin<5)
   method='laplacianhc';
end
if(nargin<4)
   if(nargin<3)
      iter=10;
   end
   alpha=0.5;
end

conn=meshconn(face,size(node,1));
newnode=smoothsurf(node(:,1:3),[],conn,iter,alpha,method,alpha);

repo: REST


text:
%% Display several images on the same figure
%  Changed: Dec 9th, 2011
%
function showimage(varargin)

    % Check parameters
    nb_args   = size(varargin,2);
    nb_images = nb_args;
    nb_cols   = 0;
    nb_rows   = 1;
    row       = 1;
    crange    = [0 1]; % default image intensities

    for i=1:nb_args
        if ischar(varargin{i})
            if isequal(varargin{i},'lim')
                lim       = varargin{i+1};
                nb_images = nb_images-2;
            elseif isequal(varargin{i},'nbcols')
                nb_cols   = varargin{i+1};
                nb_images = nb_images-2;
            elseif isequal(varargin{i},'nbrows')
                nb_rows   = varargin{i+1};
                nb_images = nb_images-2;
            elseif isequal(varargin{i},'row')
                row       = varargin{i+1};
                if row>nb_rows; nb_rows = row; end;
                nb_images = nb_images-2;
            elseif isequal(varargin{i},'caxis')
                crange    = varargin{i+1};
                nb_images = nb_images-2;
            else
                nb_images = nb_images-1;
            end
        end
    end

    if nb_cols==0; nb_cols = nb_images; end;

    % Display images
    iter_image = 1;
    for iter_arg=1:nb_args
        if ~ischar(varargin{iter_arg})
            I = varargin{iter_arg};

            % Use mid slice as image
            slice = ceil(size(I,3)/2);
            I = I(:,:,slice);

            subplot(nb_rows,nb_cols,(row-1)*nb_cols + iter_image);
            imagesc(I,crange);
            daspect([1 1 1]);
            if exist('lim'); axis([lim(3)-0.5 lim(4)+0.5 lim(1)-0.5 lim(2)+0.5]); end
            axis off;
            if iter_arg+1<=nb_args && ischar(varargin{iter_arg+1})
                title(varargin{iter_arg+1});
            end
            iter_image = iter_image+1;
        end
        if iter_image>nb_images
            break;
        end
    end

end


text:
function save_hdf5(key,outfilename)
% Saves the Key to an hdf5 file.
% Inputs:
%   key: An object of type Key.
%   outfilename: The name for the output hdf5 file.

assert(nargin==2)
assert(isa(key,'Key'))
assert(isa(outfilename,'char'))
assert(key.validate())

nummods = length(key.modelset);
numsegs = length(key.segset);

assert(size(key.tar,1)==nummods)
assert(size(key.tar,2)==numsegs)

assert(size(key.non,1)==nummods)
assert(size(key.non,2)==numsegs)

tar = int8(key.tar);
non = int8(key.non);

trialmask = tar - non;

hdf5write(outfilename,'/trial_mask',trialmask,'V71Dimensions',true);
hdf5write(outfilename,'/file_format','2.0','V71Dimensions',true,'WriteMode','append');
hdf5write(outfilename,'/ID/row_ids',key.modelset,'V71Dimensions',true,'WriteMode','append');
hdf5write(outfilename,'/ID/column_ids',key.segset,'V71Dimensions',true,'WriteMode','append');

author: nesl


text:
function [f,file] = stripfile(x,flag,sep)

% function [f,file] = stripfile(x,flag,sep)
%
% <x> is a string referring to a file (it is okay if the file
%   does not actually exist).  if <x> ends in /, we automatically
%   act as if that / does not exist.
% <flag> (optional) is whether to swap the output arguments.  default: 0.
% <sep> (optional) is the actual "/" character to use.  default is
%   the output of filesep.m.
%
% if <flag> is 0,
%   return <f> as the string but with the file name removed.
%   return <file> with the file name.
% if <flag> is 1, these arguments are swapped.
%
% history:
% 2014/07/14 - make more general by defaulting <sep> to filesep.m.
%
% example:
% isequal(stripfile('blah/temp.png',[],'/'),'blah/')
% isequal(stripfile('temp.png',[],'/'),'')
% isequal(stripfile('ok/blah/',1,'/'),'blah')

% input
if ~exist('flag','var') || isempty(flag)
  flag = 0;
end
if ~exist('sep','var') || isempty(sep)
  sep = filesep;
end

% find any /
locs = strfind(x,sep);

% if none, return ''
if isempty(locs)
  f = '';
  file = x;
% otherwise, return entire string up to the last /
else
  if locs(end)==length(x)  % ignore trailing /
    x = x(1:end-1);
    locs = locs(1:end-1);
  end
  if isempty(locs)
    f = '';
    file = x;
  else
    f = x(1:locs(end));
    file = x(locs(end)+1:end);
  end
end

% swap the output?
if flag
  [f,file] = swap(f,file);
end

author: cvnlab

repo: GLMsingle

save_path: github-repos/MATLAB/cvnlab-GLMsingle


text:
classdef RWMOP3 < PROBLEM
% <multi> <real> <constrained>
% Two bar truss design problem

%------------------------------- Reference --------------------------------
% A. Kumar, G. Wu, M. Ali, Q. Luo, R. Mallipeddi, P. Suganthan, and S. Das,
% A benchmark-suite of real-world constrained multi-objective optimization
% problems and some baseline results, Swarm and Evolutionary Computation,
% 2021, 67: 100961.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Initialization
        function Setting(obj)
            obj.M        = 2;
            obj.D        = 3;
            obj.lower    = [1e-5,1e-5,1];
            obj.upper    = [100, 100, 3];
            obj.encoding = ones(1,obj.D);
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            x  = varargin{1};
            x1 = x(:,1);
            x2 = x(:,2);
            x3 = x(:,3);
            % Objective function
            f(:,1) = x1.*(16+x3.^2).^(0.5)+x2.*(1+x3.^2).^(0.5);
            f(:,2) = (20.*(16+x3.^2).^(0.5))./(x3.*x1);
            % Constraints
            g(:,1) = f(:,1)-0.1;
            g(:,2) = f(:,2)-1e5;
            g(:,3) = (80.*(1+x3.^2).^(0.5))./(x3.*x2)-1e5;
            Population = SOLUTION(varargin{1},f,g,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
         %% Generate a point for hypervolume calculation
        function R = GetOptimum(obj,~)
            R = [1.0000000e-01   1.0000000e+05];
        end
    end
end


text:
function Results = MTrick(TrainX,TrainY,TestX,TestY,alpha,beta,numK,numCircle)

for id = 1:length(TrainY)
    if TrainY(id) == 2
        TrainY(id) = -1;
    end
end

for id = 1:length(TestY)
    if TestY(id) == 2
        TestY(id) = -1;
    end
end

G0 = [];
for i = 1:length(TrainY)
    if TrainY(i) == 1
        G0(i,1) = 1;
        G0(i,2) = 0;
    else
        G0(i,1) = 0;
        G0(i,2) = 1;
    end
end

TrainXY = scale_cols(TrainX,TrainY);
fprintf('......start to train logistic regression model1111.........\n');
w00 = zeros(size(TrainXY,1),1);
lambda = exp(linspace(-0.5,6,20));
wbest = [];
f1max = -inf;
for i = 1:length(lambda)
    w_0 = train_cg(TrainXY,w00,lambda(i));
    f1 = logProb(TrainXY,w_0);
    if f1 > f1max
        f1max = f1;
        wbest = w_0;
        se_lambda = lambda(i);
    end
end
% csvwrite(strcat('model_','test','.model'),wbest);
% wbest = load(strcat('model_','test','.model'));

ptemp = 1./(1 + exp(-wbest'*TrainX));
oriA = getResult(ptemp,TrainY);
fprintf('Test accuracy on source domain is :%g\n',oriA);
ptemp = 1./(1 + exp(-wbest'*TestX));
oriA = getResult(ptemp,TestY);
fprintf('Test accuracy on target domain is :%g\n',oriA);

fprintf('......start to learn PLSA model.........\n');
DataSetX = [TrainX TestX];
% set some variables
Learn.Verbosity = 1;
Learn.Max_Iterations = 50;
Learn.heldout = .1; % for tempered EM only, percentage of held out data
Learn.Min_Likelihood_Change = 1;
Learn.Folding_Iterations = 20; % for TEM only: number of fiolding
% in iterations
Learn.TEM = 0; %tempered or not tempered

[Pw_z,Pz_d,Pd,Li,perp,eta] = pLSA(DataSetX,[],numK,Learn);
pz = Pz_d*Pd';
pw = Pw_z*pz;
A = Pw_z;
for i = 1:size(Pw_z,1)
    A(i,:) = A(i,:).*pz';
end

for i = 1:size(Pw_z,2)
    for j = 1:length(A(:,i))
        if pw(j) > 0
            A(j,i) = A(j,i)./pw(j);
        else
            A(j,i) = 1/size(Pw_z,2);
        end
    end
end
pwz = A;
clear A;
% csvwrite(strcat('pwz_common.pwz'),pwz);
%
% pwz = load(strcat('pwz_common.pwz'));

Fs = pwz;
Ft = Fs;

Gs = G0;
Gt = [];
for i = 1:length(TestY)
    Gt(i,1) = ptemp(i);
    Gt(i,2) = 1 - ptemp(i);
end

Xs = TrainX;
Xt = TestX;
Xs = Xs/sum(sum(Xs));
Xt = Xt/sum(sum(Xt));

b = 1/(size(Gs,1));

S = ones(size(Fs,2),size(Gs,2));
for i = 1:size(S,1)
    S(i,:) = S(i,:)/sum(S(i,:));
end

fvalue = trace(Xs'*Xs-2*Xs'*Fs*S*Gs'+Gs*S'*Fs'*Fs*S*Gs')+alpha*b*trace(Gs*Gs'-2*Gs*G0'+G0*G0')+beta*trace(Xt'*Xt-2*Xt'*Ft*S*Gt'+Gt*S'*Ft'*Ft*S*Gt');
tempf = 0;
for circleID = 1:numCircle
    tempM = (Fs*S*Gs'*Gs*S');
    tempM1 = Xs*Gs*S';
    for i = 1:size(Fs,1)
        for j = 1:size(Fs,2)
            if tempM(i,j)~=0
                Fs(i,j) = Fs(i,j)*(tempM1(i,j)/tempM(i,j))^(0.5);
            else
                Fs(i,j) = 0;
            end
        end
    end
    for i = 1:size(Fs,1)
        if sum(Fs(i,:))~= 0
            Fs(i,:) = Fs(i,:)/sum(Fs(i,:));
        else
            for j = 1:size(Fs,2)
                Fs(i,j) = 1/(size(Fs,2));
            end
        end
    end
    tempM = (Gs*S'*Fs'*Fs*S+alpha*b*Gs);
    tempM1 = Xs'*Fs*S + alpha*b*G0;
    for i = 1:size(Gs,1)
        for j = 1:size(Gs,2)
            if tempM(i,j)~=0
                Gs(i,j) = Gs(i,j)*(tempM1(i,j)/tempM(i,j))^(0.5);
            else
                Gs(i,j) = 0;
            end
        end
    end
    for i = 1:size(Gs,1)
        if sum(Gs(i,:))~= 0
            Gs(i,:) = Gs(i,:)/sum(Gs(i,:));
        else
            for j = 1:size(Gs,2)
                Gs(i,j) = 1/(size(Gs,2));
            end
        end
    end

    tempM = (Ft*S*Gt'*Gt*S');
    tempM1 = Xt*Gt*S';
    for i = 1:size(Ft,1)
        for j = 1:size(Ft,2)
            if tempM(i,j)~=0
                Ft(i,j) = Ft(i,j)*(tempM1(i,j)/tempM(i,j))^(0.5);
            else
                Ft(i,j) =0;
            end
        end
    end
    for i = 1:size(Ft,1)
        if sum(Ft(i,:))~= 0
            Ft(i,:) = Ft(i,:)/sum(Ft(i,:));
        else
            for j = 1:size(Ft,2)
                Ft(i,j) = 1/(size(Ft,2));
            end
        end
    end

    tempM = (Gt*S'*Ft'*Ft*S);
    tempM1 = Xt'*Ft*S;
    for i = 1:size(Gt,1)
        for j = 1:size(Gt,2)
            if tempM(i,j)~=0
                Gt(i,j) = Gt(i,j)*(tempM1(i,j)/tempM(i,j))^(0.5);
            else
                Gt(i,j) = 0;
            end
        end
    end
    for i = 1:size(Gt,1)
        if sum(Gt(i,:))~= 0
            Gt(i,:) = Gt(i,:)/sum(Gt(i,:));
        else
            for j = 1:size(Gt,2)
                Gt(i,j) = 1/(size(Gt,2));
            end
        end
    end

    tempM = (Fs'*Fs*S*Gs'*Gs+beta*Ft'*Ft*S*Gt'*Gt);
    tempM1 = Fs'*Xs*Gs+beta*Ft'*Xt*Gt;
    for i = 1:size(S,1)
        for j = 1:size(S,2)
            if tempM(i,j)~=0
                S(i,j) = S(i,j)*(tempM1(i,j)/tempM(i,j))^(0.5);
            else
                S(i,j) = 0;
            end
        end
    end

    fvalue = trace(Xs'*Xs-2*Xs'*Fs*S*Gs'+Gs*S'*Fs'*Fs*S*Gs')+alpha*b*trace(Gs*Gs'-2*Gs*G0'+G0*G0')+beta*trace(Xt'*Xt-2*Xt'*Ft*S*Gt'+Gt*S'*Ft'*Ft*S*Gt');

    pp = [];
    for i = 1:length(TestY)
        if sum(Gt(i,:))~= 0
            pp(1,i) = Gt(i,1)/sum(Gt(i,:));
        else
            pp(1,i) = 0.5;
        end
    end
    Results(circleID) = getResult(pp,TestY);
    fprintf('the %g iteration is %g,the value of objective is %g\n',circleID,getResult(pp,TestY),fvalue);

    if circleID == 1
        tempf = fvalue;
    end
    if circleID > 1
        if abs(tempf - fvalue) < 10^(-11)
            break;
        end
        tempf = fvalue;
    end
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function p = towardsback(p0, vp, amount)
dir = vp{3} - p0;
dir = dir / norm(dir);
p = p0 + amount * dir;

author: zouchuhang

repo: LayoutNet


text:
% pop_export() - export EEG dataset
%
% Usage:
%   >> com = pop_export(EEG);   % a window pops up
%   >> com = pop_export(EEG, filename, 'key', 'val', ... );
%
% Inputs:
%   EEG            - eeglab dataset
%   filename       - file name
%
% Optional inputs:
%   'ica'          - ['on'|'off'] export ICA activities (or ERP). Default 'off'.
%   'time'         - ['on'|'off'] include time values. Default 'on'.
%   'timeunit'     - [float] time unit rel. to seconds. Default: 1E-3 = msec.
%   'elec'         - ['on'|'off'] include electrodes names or component numbers.
%                    Default 'on'.
%   'transpose'    - ['on'|'off'] 'off'-> electrode data = rows; 'on' -> electrode
%                    data = columns. Default 'off'.
%   'erp'          - ['on'|'off'] export ERP instead of raw data. Default 'off'.
%   'expr'         - [string] evaluate epxression on data. The expression must
%                    contain a variable 'x' representing the 2-D or 3-D data.
%   'precision'    - [float] number of significant digits in output. Default 4.
%
% Outputs:
%   com            - The expresion that execute this function. i.e. 'pop_export(MyEEG, 'ExpEEG.mat')'
%
% Note: tabulation are used as a delimiter.
%
% Author: Arnaud Delorme, CNL / Salk Institute, May 13, 2003

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

function com = pop_export(EEG, filename, varargin);

com = '';
if nargin < 1
    help pop_export;
    return;
end;

if nargin < 2
    commandload = [ '[filename, filepath] = uiputfile(''*'', ''Select a text file'');' ...
                    'if filename ~=0,' ...
                    '   set(findobj(''parent'', gcbf, ''tag'', ''tagedit''), ''string'', [ filepath filename ]);' ...
                    'end;' ...
                    'clear filename filepath tagtest;' ];

   uilist = { { 'style' 'text' 'string' 'Output file name' }, ...
              { 'style' 'edit' 'string' '' 'tag' 'tagedit' }, ...
              { 'style' 'pushbutton' 'string' 'Browse' 'callback' commandload },...
              { 'style' 'text' 'string' 'Export ICA activities instead of EEG data:' }, ...
              { 'style' 'checkbox' 'string' '' }, { },{ },  ...
              { 'style' 'text' 'string' 'Export ERP average instead of trials:' }, ...
              { 'style' 'checkbox' 'string' '' }, { }, { }, ...
              { 'style' 'text' 'string' 'Transpose matrix (elec -> rows):' }, ...
              { 'style' 'checkbox' 'string' '' }, { }, { }, ...
              { 'style' 'text' 'string' 'Export channel labels/component numbers:' }, ...
              { 'style' 'checkbox' 'string' '' 'value' 1 }, { }, { }, ...
              { 'style' 'text' 'string' 'Export time values:' }, ...
              { 'style' 'checkbox' 'string' '' 'value' 1 }, ...
              { 'style' 'text' 'string' 'Unit (re. sec)' }, { 'style' 'edit' 'string' '1E-3' }, ...
              { 'style' 'text' 'string' 'Number of significant digits to output:' }, ...
              { 'style' 'edit' 'string' '4' }, { }, { }, ...
              { 'style' 'text' 'string' 'Apply an expression to the output (see ''expr'' help):'} , ...
              { 'style' 'edit' 'string' '' } { } { } };
   bcheck = [1.7 0.25 0.7 0.6];
   uigeom = { [1 2 0.5] bcheck bcheck bcheck bcheck bcheck [3 0.7 0.8 0.6] [3 2 0.05 0.05] };
   result = inputgui( uigeom, uilist, 'pophelp(''pop_export'');', 'Export data - pop_export()' );
   if length( result ) == 0 return; end;

   % decode options
   % --------------
   if isempty(result{1}), error('File name required'); end;
   filename = result{1};
   options = {};
   if result{2},  options = { options{:} 'ica' 'on' }; end;
   if result{3},  options = { options{:} 'erp' 'on' }; end;
   if result{4},  options = { options{:} 'transpose' 'on' }; end;
   if ~result{5}, options = { options{:} 'elec' 'off' }; end;
   if ~result{6}, options = { options{:} 'time' 'off' }; end;
   if ~strcmpi(result{7}, '1E-3'), options = { options{:} 'timeunit' eval(result{7}) }; end;
   if ~strcmpi(result{8}, '4'),    options = { options{:} 'precision' eval(result{8}) }; end;
   if ~isempty(result{9}), options = { options{:} 'expr' result{9} }; end;
else
    options = varargin;
end;

% test inputs
% -----------
g = finputcheck(options, { ...
    'ica'       'string'    { 'on' 'off' }     'off';
    'time'      'string'    { 'on' 'off' }     'on';
    'timeunit'  'float'     [0 Inf]            1E-3;
    'elec'      'string'    { 'on' 'off' }     'on';
    'transpose' 'string'    { 'on' 'off' }     'off';
    'erp'       'string'    { 'on' 'off' }     'off';
    'precision' 'integer'   [0 Inf]            4;
    'expr'      'string'    []                 '' }, 'pop_export');
if isstr(g), error(g); end;

% select data
% ----------
if strcmpi(g.ica, 'on');
    eeglab_options; % changed from eeglaboptions 3/30/02 -sm
	if option_computeica
		x = EEG.icaact;
	else
        x = EEG.icaweights*EEG.icasphere*reshape(EEG.data(EEG.icachansind,:,:), length(EEG.icachansind), EEG.trials*EEG.pnts);
        x = reshape(x, size(x,1), EEG.pnts, EEG.trials);
    end;
else
    x = EEG.data;
end;

% select erp
% ----------
if strcmpi(g.erp, 'on');
    x = mean(x, 3);
else
    x = reshape(x, size(x,1), size(x,2)*size(x,3));
end;

% write data
% ----------
if ~isempty(g.expr)
    eval([ g.expr ';' ]);
end;

% add time axis
% -------------
if strcmpi(g.time, 'on');
    timeind = repmat( linspace(EEG.xmin, EEG.xmax, EEG.pnts)/g.timeunit, ...
                      [ 1 fastif(strcmpi(g.erp,'on'), 1, EEG.trials) ]);
    xx = zeros(size(x,1)+1, size(x,2));
    xx(1,:)     = timeind;
    xx(2:end,:) = x;
    x = xx; clear xx;
end

% transpose and write to disk
% ---------------------------
fid = fopen(filename, 'w');
if strcmpi(g.transpose, 'on');
    % writing electrodes
    % ------------------
    strprintf = '';
    for index = 1:size(x,1)
        if strcmpi(g.time, 'on'), tmpind = index-1;
        else                      tmpind = index;
        end;
        if strcmpi(g.elec, 'on')
            if tmpind > 0
                if ~isempty(EEG.chanlocs) & ~strcmpi(g.ica, 'on')
                    fprintf(fid, '%s\t', EEG.chanlocs(tmpind).labels);
                else fprintf(fid, '%d\t', tmpind);
                end;
            else
                fprintf(fid, ' \t');
            end;
        end;
        strprintf = [ strprintf '%.' num2str(g.precision) 'f\t' ];
    end;
    strprintf(end) = 'n';
    if strcmpi(g.elec, 'on'), fprintf(fid, '\n'); end;
    fprintf(fid, strprintf, x);
else
    % writing electrodes
    % ------------------
    for index = 1:size(x,1)
        if strcmpi(g.time, 'on'), tmpind = index-1;
        else                      tmpind = index;
        end;
        if strcmpi(g.elec, 'on')
            if tmpind > 0
                if ~isempty(EEG.chanlocs) & ~strcmpi(g.ica, 'on')
                    fprintf(fid,'%s\t', EEG.chanlocs(tmpind).labels);
                else fprintf(fid,'%d\t', tmpind);
                end;
            else
                fprintf(fid, ' \t');
            end;
        end;
        fprintf(fid,[ '%.' num2str(g.precision) 'f\t' ], x(index, :));
        fprintf(fid, '\n');
    end;
end;
fclose(fid);

com = sprintf('pop_export(%s,%s);', inputname(1), vararg2str({ filename, options{:} }));
return;

author: PatternRecognition

repo: OpenBMI


text:
function updateGeoScatter(obj,geoIndex)

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

    %-INTIALIZATIONS-%

    axIndex = obj.getAxisIndex(obj.State.Plot(geoIndex).AssociatedAxis);
    geoData = get(obj.State.Plot(geoIndex).Handle);
    axisData = geoData.Parent;
    [xSource, ~] = findSourceAxis(obj,axIndex);

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

    %-set trace-%
    if strcmpi(obj.PlotOptions.geoRenderType, 'geo')
        obj.data{geoIndex}.geo = sprintf('geo%d', xSource+1);
        obj.data{geoIndex}.type = 'scattergeo';

    elseif strcmpi(obj.PlotOptions.geoRenderType, 'mapbox')
        obj.data{geoIndex}.subplot = sprintf('mapbox%d', xSource+1);
        obj.data{geoIndex}.type = 'scattermapbox';
    end

    obj.data{geoIndex}.mode = 'markers+text';

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

    %-set trace data-%
    obj.data{geoIndex}.lat = geoData.LatitudeData;
    obj.data{geoIndex}.lon = geoData.LongitudeData;

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

    %-set trace marker-%
    marker = extractGeoMarker(geoData, axisData);

    if strcmpi(obj.PlotOptions.geoRenderType, 'mapbox')
        marker.allowoverlap = true;
        marker = rmfield(marker, 'symbol');

        if strcmp(marker.color, 'rgba(0,0,0,0)') && isfield(marker, 'line')
            marker.color = marker.line.color;
        end
    end

    obj.data{geoIndex}.marker = marker;

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

author: plotly


text:
function pco_camera_create_deffile()
 wfid = fopen('pco_camera_def.txt','wt');
 if(wfid == -1 )
  error('Unable to create file: pco_camera_def.txt');
 end

 get_defines('PCO_err.h',wfid);
 get_defines('sc2_defs.h',wfid);
 get_defines('sc2_ml_sdkstructures.h',wfid);
 get_defines('PCO_Recorder_Defines.h',wfid);

 fclose(wfid);
end

function get_defines(hfile,wfid)
n = 0;
fid = fopen(hfile);
if( fid == -1 )
    error(['Unable to open file: ' hfile]);
end
while( true )
    tline = fgetl(fid);
    if( tline == -1 )
     break;
    end
    if(numel(tline)<2)
%     disp(['blank line ',tline]);
     continue;
    end
    if((tline(1)=='/')&&(tline(2)=='/'))
%     disp(['found comment ',tline]);
     continue;
    end
    d = strfind(tline,'#define ');
    if( ~isempty(d) )
        str = strtrim(tline(d(1)+8:end));
%        disp(['string is ',str]);
        s = strfind(str,' ');
        if( ~isempty(s) )
            p = strfind(str(1:s(1)-1),'(');
            if( isempty(p) )
             q = strfind(str,'0x');
             r = strfind(str,'/');
             if( isempty(q) )
              str(s(1)) = '=';
              if( ~isempty(r) )
               str=str(1:r-1);
              end
              s=strfind(str,'.');
              if( ~isempty(s) )
               str=str(1:end-1);
              end
               fprintf(wfid,'%s\n',str);
%              disp(['string is ',str]);
%              eval(str);
               n = n + 1;
             else
              hstr=str(1:q-1);
              hstr=strtrim(hstr);
              hstr=strcat(hstr,'=');
              if( ~isempty(r) )
               numstr=str(q+2:r-1);
              else
               numstr=str(q+2:end);
              end
              numstr=strtrim(numstr);
              num=hex2dec(numstr);
              str=strcat(hstr,num2str(num));
              fprintf(wfid,'%s\n',str);

%              disp(['string is ',hstr,' numstr is ',numstr,' num is ',num2str(num)]);
              n = n + 1;
             end
            end
        end
    end
end
disp(['Number of variables defined in ',hfile,' = ' num2str(n)]);
fclose(fid);

end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function [seg, SEG_l] = ahmLinSegment(l,t)

% AHMLINSEGMENT  AHM line endpoints.
%   AHMLINSEGMENT(L,T) returns the 3D segment corresponding to the HMG line
%   L at abscissas T = [t1;t2].
%
%   [s,S_l] = AHMLINSEGMENT(...) returns the Jacobian wrt L.
%
%   See also AHMLINENDPOINTS, AHMLIN2SEG, AHMLIN2IDPPNTS.

%   Copyright 2009 Teresa Vidal.

% abscissas
t1 = t(1);
t2 = t(2);

% support 3d segment
[s, S_l] = ahmLin2seg(l);

% support points
p1 = s(1:3);
p2 = s(4:6);

P1_l = S_l(1:3,:);
P2_l = S_l(4:6,:);

% 3d endpoints
e1 = (1-t1)*p1 + t1*p2;
e2 = (1-t2)*p1 + t2*p2;

E1_p1 = 1-t1;
E1_p2 = t1;
E2_p1 = 1-t2;
E2_p2 = t2;

E1_l = E1_p1*P1_l + E1_p2*P2_l;
E2_l = E2_p1*P1_l + E2_p2*P2_l;

% 3d segment
seg   = [e1;e2];
SEG_l = [E1_l;E2_l];

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [] = twod_to_tecplot(x,e_conn,u,filename,title)
%------------------------------------------------------------------------------
%  twod_to_tecplot - function that accepts 2D finite element data and writes
%                    and input file to the graphics program Tecplot.
%
%  Usage:
%         [] = twod_to_tecplot(x,e_conn,u,filename,title)
%
%  Variables:     x
%                      spatial coordinates  (nodes x 2)
%                 e_conn
%                      element connectivity matrix  (3 nodes are used)
%                 u
%                      field variables
%                 filename
%                      name for the output tecplot file
%                 title
%                      title string required by tecplot
%------------------------------------------------------------------------------

[n_node,n_var] = size(u);
[n_elem,n_dof] = size(e_conn);

[fid] = fopen(filename,'w');

fprintf(fid,strcat(strcat('Title="',title),'"\n'));

temp_str = 'Variables="x","y"';
for i=1:n_var
  temp_str = strcat(strcat(temp_str,strcat(',"vec',strcat(num2str(i),'"'))));
end
temp_str = strcat(temp_str,'\n');
fprintf(fid,temp_str);

temp_str = 'Zone N = ';
temp_str = strcat( temp_str, num2str(n_node) );
temp_str = strcat( temp_str, ', E = ' );
temp_str = strcat( temp_str, num2str(n_elem) );
temp_str = strcat( temp_str, ', DATAPACKING = POINT, ZONETYPE = FETRIANGLE\n' );
fprintf(fid,temp_str);

print_str = ' %10.5g  %10.5g ';
for i=1:n_var
  print_str = strcat( print_str, ' %10.5g ' );
end
print_str = strcat( print_str, '\n' );

for i=1:n_node
  fprintf(fid,print_str,x(i,:),u(i,:));
end

print_str = ' %6d  %6d  %6d \n';
for i=1:n_elem
  fprintf(fid,print_str,e_conn(i,1:3));
end

fclose(fid);
% end function output_to_tecplot

author: johannesgerer

repo: jburkardt-m


text:
function output = callqpip(interfacedata)

% Retrieve needed data
options = interfacedata.options;
model = yalmip2quadprog(interfacedata);

if options.savedebug
    model.ops = options.qpip;
    save debugfile model
end

if options.showprogress;showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);end
solvertime = tic;
[x,flag,lm] = qpip(model.Q, model.c, model.A, model.b, model.Aeq, model.beq, model.lb, model.ub, options.verbose,options.qpip.mu,options.qpip.method);
solvertime = toc(solvertime);

% Internal format for duals
if ~isempty(lm)
    D_struc = [lm.equality;lm.inequality];
else
    D_struc = [];
end

% Check, currently not exhaustive...
switch flag
    case 0
        problem = 0;
    case 2
        problem = 1;
    otherwise
        problem = 1;
end

infostr = yalmiperror(problem,'QPIP');

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

% Save all data from the solver?
if options.savesolveroutput
    solveroutput.x = x;
    solveroutput.flag = flag;
    solveroutput.lm = lm;
else
    solveroutput = [];
end

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


text:
function [x, err] = rombergQuadrature(fun,tSpan,tol)
% [x, err] = rombergQuadrature(fun,tSpan,tol)
%
% Compute the integral(fun), over the domain tSpan, to an accuracy of tol,
% using Romberg quadrature. Fully vectorized.
%
% Good for high-accuracy quadrature over smooth vector functions.
%
% If necessary, this function will automatically sub-divide the interval to
% achieve the desired accuracy. This should only occur when fun is stiff or
% non-smooth.
%
% INPUTS:
%   fun = vector function to be integrated
%       dx = fun(t)
%           t = [1, nt] = time vector
%           dx = [nx, nt] = function value at each point in t
%   tSpan = [tLow, tUpp] = time span (domain) for integration
%   tol = [nx,1] = desired error tolerance along each dimension
%
% OUTPUT:
%   x = [nx,1] = integral along each dimension
%   err = [nx, 1] = error estimate along each dimension
%
% NOTES:
%   algorithm from:
%   http://www.math.usm.edu/lambers/mat460/fall09/lecture29.pdf
%

a = tSpan(1);
b = tSpan(2);
h = b-a;
f0 = fun(tSpan);
nx = size(f0,1);

if isscalar(tol)
    tol = tol*ones(nx,1);
end

nMin = 4;   % Complete at least this many iterations
nMax = 12;  % Maximum iteration count
nSubSegment = 4;  % Sub-divide segment if max iteration is reached

T = zeros(nMax,nMax,nx);  % (iteration, extrapolation, dimension)
Err = zeros(nMax,nx); %(extrapolation, dimension)

T(1,1,:) = 0.5*h*sum(f0,2);
h = h/2;

converged = false;
for j=2:nMax

    % Trapazoid method
    i = 1:2^(j-2);
    t = a+h*(2*i-1);
    f = fun(t);
    T(j,1,:) = 0.5*T(j-1,1,:) + reshape(h*sum(f,2),1,1,nx);

    % Richardson extrapolation
    for k=2:j
        T(j,k,:) = T(j,k-1,:) + ...
            ( 1/(4^(k-1)-1) )*( T(j,k-1,:) - T(j-1,k-1,:) );
    end
    h = h/2;

    % Check the error estimate for convergence
    Err(j,:) = T(j,k,:) - T(j-1,k-1,:);
    err = abs(Err(j,:))';
    if all( err < tol ) && j > nMin
        converged = true;
        break;
    end
end

x = reshape(T(j,j,:),nx,1);

if ~converged  %Then non-smooth problem. Recursively sub-divide interval.
    time = linspace(tSpan(1), tSpan(2), nSubSegment+1);
    x = zeros(nx,1);
    err = zeros(nx,1);
    for i=1:nSubSegment
        [xTmp, errTmp] = rombergQuadrature(fun,time([i, i+1]),tol);
        x = x + xTmp;
        err = max(err,errTmp);
    end
end

end

author: MatthewPeterKelly

repo: OptimTraj