text:
%
% LibQPEP: A Library for Globally Optimal Solving Quadratic Pose Estimation Problems (QPEPs),
%          It also gives highly accurate uncertainty description of the solutions.
%
%
% Article:
%      Wu, J., Zheng, Y., Gao, Z., Jiang, Y., Hu, X., Zhu, Y., Jiao, J., Liu, M. (2020)
%           Quadratic Pose Estimation Problems: Globally Optimal Solutions,
%           Solvability/Observability Analysis and Uncertainty Description.
%           IEEE Transactions on Robotics.
%           https://doi.org/10.1109/TRO.2022.3155880
%
%
% Authors:      Jin Wu and Ming Liu
% Affiliation:  Hong Kong University of Science and Technology (HKUST)
% Emails:       jin_wu_uestc@hotmail.com; eelium@ust.hk
% Websites:     https://zarathustr.github.io
%               https://ram-lab.com
%
%
% syms_pnp.m: Generating symbolic expressions for solving PnP problem

clear all
close all
clc

addpath('func_files');
addpath('utils');

syms b1 b2 real
syms r1 r2 r3 real

image_pt = [b1, b2];
world_pt = [r1, r2, r3];

syms fx fy cx cy real
K = [fx, 0, 0;
      0, fy, 0;
     cx, cy, 1];

syms q0 q1 q2 q3 real
syms t1 t2 t3 lambda scale real
q = [q0; q1; q2; q3];
tt = [t1; t2; t3];
RR = q2R(q);
XX = [
    RR, tt;
    zeros(1, 3), 1];
J_pure = J_pnp_loss(image_pt, world_pt, K, RR, tt) * scale;
[coef_J_pure, mon_J_pure] = coeffs(J_pure, [q; tt]);
generateFuncFile(coef_J_pure, fullfile('func_files', 'coef_J_pure_pnp_func_new.m'), {image_pt, world_pt, [fx, fy, cx, cy], scale});
generateFuncFile(mon_J_pure, fullfile('func_files', 'mon_J_pure_pnp_func_new.m'), {q, tt});
J = J_pure - 0.5 * lambda * (q.' * q - 1);
x = [q; tt; lambda];
Jacob = jacobian(J, x);
Jacob = expand(Jacob.');

coef_len = length(coef_J_pure);
str = sprintf('coef_J = sym(''coef_J'', [1, %d]);', coef_len);
eval(str);
showSyms(symvar(coef_J));
J = coef_J * mon_J_pure.' - 0.5 * lambda * (q.' * q - 1);
x = [q; tt; lambda];
Jacob = jacobian(J, x);
Jacob = expand(Jacob.');

ts = solveSyms(Jacob(5 : 7), tt);
ss = Jacob(1 : 4);
ss = subs(ss, q0^3, (q0 * (1 - q1^2 - q2^2 - q3^2)));
ss = subs(ss, q0^2, (1 - q1^2 - q2^2 - q3^2));
Jacob(1 : 4) = ss;

ss = subs(ss, t1, ts.t1);
ss = subs(ss, t2, ts.t2);
ss = subs(ss, t3, ts.t3);

[coeft1, mont1] = coeffs(Jacob(5), tt);
g1_ = coeft1(1 : 3).';
g1 = vpa(zeros(3, 1));
for i = 1 : 3
    idx = sscanf(char(mont1(i)), 't%d');
    if(~isempty(idx))
        g1(idx) = g1_(i);
    end
end
[coeft2, mont2] = coeffs(Jacob(6), tt);
g2_ = coeft2(1 : 3).';
g2 = vpa(zeros(3, 1));
for i = 1 : 3
    idx = sscanf(char(mont2(i)), 't%d');
    if(~isempty(idx))
        g2(idx) = g2_(i);
    end
end
[coeft3, mont3] = coeffs(Jacob(7), tt);
g3_ = coeft3(1 : 3).';
g3 = vpa(zeros(3, 1));
for i = 1 : 3
    idx = sscanf(char(mont3(i)), 't%d');
    if(~isempty(idx))
        g3(idx) = g3_(i);
    end
end

[coeftq1, montq1] = coeffs(expand(Jacob(5) - g1.' * tt ), q);
generateFuncFile(coeftq1, fullfile('func_files', 'coeftq1_pnp_func_new.m'), {coef_J});
[coeftq2, montq2] = coeffs(expand(Jacob(6) - g2.' * tt ), q);
generateFuncFile(coeftq2, fullfile('func_files', 'coeftq2_pnp_func_new.m'), {coef_J});
[coeftq3, montq3] = coeffs(expand(Jacob(7) - g3.' * tt ), q);
generateFuncFile(coeftq3, fullfile('func_files', 'coeftq3_pnp_func_new.m'), {coef_J});
G = - [g1.'; g2.'; g3.'];
generateFuncFile(G, fullfile('func_files', 'G_pnp_func_new.m'), {coef_J});

pinvG = sym('pinvG', [3, 3]);
showSyms(symvar(pinvG));
coefs_tq = sym('coefs_tq', [3, 10]);
showSyms(symvar(coefs_tq));
ts = pinvG * coefs_tq * montq1.';
t1 = ts(1);
t2 = ts(2);
t3 = ts(3);
generateFuncFile(t1, fullfile('func_files', 't1_pnp_func_new.m'), {pinvG, coefs_tq, q});
generateFuncFile(t2, fullfile('func_files', 't2_pnp_func_new.m'), {pinvG, coefs_tq, q});
generateFuncFile(t3, fullfile('func_files', 't3_pnp_func_new.m'), {pinvG, coefs_tq, q});

[coef_Jacob1_qt, mon_Jacob1_qt] = coeffs(Jacob(1) + lambda * q0, [q; tt]);
generateFuncFile(coef_Jacob1_qt, fullfile('func_files', 'coef_Jacob1_qt_pnp_func_new.m'), {coef_J});
[coef_Jacob2_qt, mon_Jacob2_qt] = coeffs(Jacob(2) + lambda * q1, [q; tt]);
generateFuncFile(coef_Jacob2_qt, fullfile('func_files', 'coef_Jacob2_qt_pnp_func_new.m'), {coef_J});
[coef_Jacob3_qt, mon_Jacob3_qt] = coeffs(Jacob(3) + lambda * q2, [q; tt]);
generateFuncFile(coef_Jacob3_qt, fullfile('func_files', 'coef_Jacob3_qt_pnp_func_new.m'), {coef_J});
[coef_Jacob4_qt, mon_Jacob4_qt] = coeffs(Jacob(4) + lambda * q3, [q; tt]);
generateFuncFile(coef_Jacob4_qt, fullfile('func_files', 'coef_Jacob4_qt_pnp_func_new.m'), {coef_J});

coef_Jacob1_qt_syms = sym('coef_Jacob1_qt_syms', [1, 32]);
coef_Jacob2_qt_syms = sym('coef_Jacob2_qt_syms', [1, 32]);
coef_Jacob3_qt_syms = sym('coef_Jacob3_qt_syms', [1, 32]);
coef_Jacob4_qt_syms = sym('coef_Jacob4_qt_syms', [1, 32]);
showSyms(symvar(coef_Jacob1_qt_syms));
showSyms(symvar(coef_Jacob2_qt_syms));
showSyms(symvar(coef_Jacob3_qt_syms));
showSyms(symvar(coef_Jacob4_qt_syms));
coef_Jacob_qt_syms = [
    coef_Jacob1_qt_syms;
    coef_Jacob2_qt_syms;
    coef_Jacob3_qt_syms;
    coef_Jacob4_qt_syms;
    ];
Jacob_ = coef_Jacob_qt_syms * mon_Jacob1_qt.';

eqs = expand(eval([Jacob_; Jacob(8)]));

eqs_ = eqs;
eqs_(1 : 4) = eval(subs(eqs_(1 : 4), q0^2, (1 - q3 * q3 - q1 * q1 - q2 * q2)));
eqs_(1 : 4) = eval(subs(eqs_(1 : 4), q0^3, q0 * (1 - q3 * q3 - q1 * q1 - q2 * q2)));

f0 = eqs(1);
f1 = eqs(2);
f2 = eqs(3);
f3 = eqs(4);
[coef_f0_q, mon_f0_q] = coeffs(f0, q);
[coef_f1_q, mon_f1_q] = coeffs(f1, q);
[coef_f2_q, mon_f2_q] = coeffs(f2, q);
[coef_f3_q, mon_f3_q] = coeffs(f3, q);
coef_f0_q_sym = sym('coef_f0_q_sym', [1, 24]);
coef_f1_q_sym = sym('coef_f1_q_sym', [1, 24]);
coef_f2_q_sym = sym('coef_f2_q_sym', [1, 24]);
coef_f3_q_sym = sym('coef_f3_q_sym', [1, 24]);
showSyms(symvar(coef_f0_q_sym));
showSyms(symvar(coef_f1_q_sym));
showSyms(symvar(coef_f2_q_sym));
showSyms(symvar(coef_f3_q_sym));
f0 = coef_f0_q_sym * mon_f0_q.';
f1 = coef_f1_q_sym * mon_f1_q.';
f2 = coef_f2_q_sym * mon_f2_q.';
f3 = coef_f3_q_sym * mon_f3_q.';
coef_f_q_sym = [
    coef_f0_q;
    coef_f1_q;
    coef_f2_q;
    coef_f3_q;
    ];
generateFuncFile(coef_f_q_sym, fullfile('func_files', 'coef_f_q_sym_pnp_func_new.m'), {pinvG, coefs_tq, coef_Jacob_qt_syms});

eq = expand([
     f0 * q1 - f1 * q0;
     f0 * q2 - f2 * q0;
     f0 * q3 - f3 * q0;
     q0 * q0 + q1 * q1 + q2 * q2 + q3 * q3 - 1;
     ]);
generateFuncFile(eq, fullfile('func_files', 'eq_pnp_func_new.m'), {coef_f0_q_sym, coef_f1_q_sym, coef_f2_q_sym, coef_f3_q_sym, q});

eq_ = eq;
[coefeq1, moneq1] = coeffs(eq_(1), q);
[coefeq2, moneq2] = coeffs(eq_(2), q);
[coefeq3, moneq3] = coeffs(eq_(3), q);

v1 = moneq1.';
generateFuncFile(v1, fullfile('func_files', 'v1_func_pnp_new.m'), {q});
v2 = moneq2.';
generateFuncFile(v2, fullfile('func_files', 'v2_func_pnp_new.m'), {q});
v3 = moneq3.';
generateFuncFile(v3, fullfile('func_files', 'v3_func_pnp_new.m'), {q});

D = [
    coefeq1;
    coefeq2;
    coefeq3;
    ];
generateFuncFile(D, fullfile('func_files', 'D_func_pnp_new.m'), {coef_f0_q_sym, coef_f1_q_sym, coef_f2_q_sym, coef_f3_q_sym});

HH = jacobian(eq, q);
generateFuncFile(HH, fullfile('func_files', 'Jacob_pnp_func_new.m'), {coef_f0_q_sym, coef_f1_q_sym, coef_f2_q_sym, coef_f3_q_sym, q});
gradient = expand(HH.' * eq);

[coef1, mon1] = coeffs(eqs_(1), x);
[coef2, mon2] = coeffs(eqs_(2), x);
[coef3, mon3] = coeffs(eqs_(3), x);
[coef4, mon4] = coeffs(eqs_(4), x);

Q_sym = [
    coef1(7), coef1(14), coef1(18), coef1(20);
    coef2(7), coef2(14), coef2(18), coef2(20);
    coef3(7), coef3(14), coef3(18), coef3(20);
    coef4(7), coef4(14), coef4(18), coef4(20);
    ];
generateFuncFile(Q_sym, fullfile('func_files', 'Q_pnp_func_new.m'), {pinvG, coefs_tq, coef_Jacob_qt_syms});

res = expand(eqs_(1 : 4) - Q_sym * q);
H = jacobian(res, q);
h1 = H(:, 1);
h2 = H(:, 2);
h3 = H(:, 3);
h4 = H(:, 4);
P1 = jacobian(h1, q);
P2 = jacobian(h2, q);
P3 = jacobian(h3, q);
P4 = jacobian(h4, q);

for i = 1 : 4
    for j = 1 : 4
        str = sprintf('W%d%d = [jacobian(P%d(1, :), q%d).''; jacobian(P%d(2, :), q%d).''; jacobian(P%d(3, :), q%d).''; jacobian(P%d(4, :), q%d).''];', ...
                        i, j, i, j - 1, i, j - 1, i, j - 1, i, j - 1);
        eval(str);
    end
end

W1 = [W11, W12, W13, W14];
W2 = [W21, W22, W23, W24];
W3 = [W31, W32, W33, W34];
W4 = [W41, W42, W43, W44];
v = kron(q, q);
u = kron(v, q);
W = - [W1, W2, W3, W4] / 6;
generateFuncFile(W, fullfile('func_files', 'W_pnp_func_new.m'), {pinvG, coefs_tq, coef_Jacob_qt_syms});

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function value = mean(sF1)
% mean value of a spherical function
%

value = mean(sF1.values(:));

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [pTS,fmin]=grTravSale(C)
% Function [pTS,fmin]=grTravSale(C) solve the nonsymmetrical
% traveling salesman problem.
% Input parameter:
%   C(n,n) - matrix of distances between cities,
%     maybe, nonsymmetrical;
%     n - number of cities.
% Output parameters:
%   pTS(n) - the order of cities;
%   fmin - length of way.
% Uses the reduction to integer LP-problem:
% Look: Miller C.E., Tucker A. W., Zemlin R. A.
% Integer Programming Formulation of Traveling Salesman Problems.
% J.ACM, 1960, Vol.7, p. 326-329.
% Needed other products: MIQP.M.
% This software may be free downloaded from site:
% http://control.ee.ethz.ch/~hybrid/miqp/
% Author: Sergiy Iglin
% e-mail: siglin@yandex.ru
% personal page: http://iglin.exponenta.ru

% ============= Input data validation ==================
if nargin<1,
  error('There are no input data!')
end
if ~isnumeric(C),
  error('The array C must be numeric!')
end
if ~isreal(C),
  error('The array C must be real!')
end
s=size(C); % size of array C
if length(s)~=2,
  error('The array C must be 2D!')
end
if s(1)~=s(2),
  error('Matrix C must be square!')
end
if s(1)<3,
  error('Must be not less than 3 cities!')
end

% ============ Size of problem ====================
n=s(1); % number of vertexes
m=n*(n-1); % number of arrows

% ============ Parameters of integer LP problem ========
Aeq=[]; % for the matrix of the boundary equations
for k1=1:n,
  z1=zeros(n);
  z1(k1,:)=1;
  z2=[z1;eye(n)];
  Aeq=[Aeq z2([1:2*n-1],setdiff([1:n],k1))];
end
Aeq=[Aeq zeros(2*n-1,n-1)];
A=[]; % for the matrix of the boundary inequations
for k1=2:n,
  z1=[];
  for k2=1:n,
    z2=eye(n)*(n-1)*(k2==k1);
    z1=[z1 z2(setdiff([2:n],k1),setdiff([1:n],k2))];
  end
  z2=-eye(n);
  z2(:,k1)=z2(:,k1)+1;
  A=[A;[z1 z2(setdiff([2:n],k1),2:n)]];
end
beq=ones(2*n-1,1); % the right parts of the boundary equations
b=ones((n-1)*(n-2),1)*(n-2); % the right parts of the boundary inequations
C1=C'+diag(ones(1,n)*NaN);
C2=C1(:);
c=[C2(~isnan(C2));zeros(n-1,1)]; % the factors for objective function
vlb=[zeros(m,1);-inf*ones(n-1,1)]; % the lower bounds
vub=[ones(m,1);inf*ones(n-1,1)]; % the upper bounds
H=zeros(n^2-1); % Hessian

% ============= We solve the MIQP problem ==========
[xmin,fmin]=MIQP(H,c,A,b,Aeq,beq,[1:m],vlb,vub);

% ============= We return the results ==============
eik=round(xmin(1:m)); % the arrows of the way
e1=[zeros(1,n-1);reshape(eik,n,n-1)];
e2=[e1(:);0]; % we add zero to a diagonal
e3=(reshape(e2,n,n))'; % the matrix of the way
pTS=[1 find(e3(1,:))]; % we build the way
while pTS(end)>1, % we add the city to the way
  pTS=[pTS find(e3(pTS(end),:))];
end
return

author: HuangCongQing


text:
classdef QuadPlot < handle
    %QUADPLOT Visualization class for quad

    properties (SetAccess = public)
        k = 0;
        qn;             % quad number
        time = 0;       % time
        state;          % state
        des_state;      % desried state [x; y; z; xdot; ydot; zdot];
        rot;            % rotation matrix body to world

        color;          % color of quad
        wingspan;       % wingspan
        height;         % height of quad
        motor;          % motor position

        state_hist        % position history
        state_des_hist;   % desired position history
        time_hist;      % time history
        max_iter;       % max iteration
    end

    properties (SetAccess = private)
        h_3d
        h_m13;  % motor 1 and 3 handle
        h_m24;  % motor 2 and 4 handle
        h_qz;   % z axis of quad handle
        h_qn;   % quad number handle
        h_pos_hist;     % position history handle
        h_pos_des_hist; % desired position history handle
        text_dist;  % distance of quad number to quad
    end

    methods
        % Constructor
        function Q = QuadPlot(qn, state, wingspan, height, color, max_iter, h_3d)
            Q.qn = qn;
            Q.state = state;
            Q.wingspan = wingspan;
            Q.color = color;
            Q.height = height;
            Q.rot = QuatToRot(Q.state(7:10));
            Q.motor = quad_pos(Q.state(1:3), Q.rot, Q.wingspan, Q.height);
            Q.text_dist = Q.wingspan / 3;
            Q.des_state = Q.state(1:6);

            Q.max_iter = max_iter;
            Q.state_hist = zeros(6, max_iter);
            Q.state_des_hist = zeros(6, max_iter);
            Q.time_hist = zeros(1, max_iter);

            % Initialize plot handle
            if nargin < 7, h_3d = gca; end
            Q.h_3d = h_3d;
            hold(Q.h_3d, 'on')
            Q.h_pos_hist = plot3(Q.h_3d, Q.state(1), Q.state(2), Q.state(3), 'r.');
            Q.h_pos_des_hist = plot3(Q.h_3d, Q.des_state(1), Q.des_state(2), Q.des_state(3), 'b.');
            Q.h_m13 = plot3(Q.h_3d, ...
                Q.motor(1,[1 3]), ...
                Q.motor(2,[1 3]), ...
                Q.motor(3,[1 3]), ...
                '-ko', 'MarkerFaceColor', Q.color, 'MarkerSize', 5);
            Q.h_m24 = plot3(Q.h_3d, ...
                Q.motor(1,[2 4]), ...
                Q.motor(2,[2 4]), ...
                Q.motor(3,[2 4]), ...
                '-ko', 'MarkerFaceColor', Q.color, 'MarkerSize', 5);
            Q.h_qz = plot3(Q.h_3d, ...
                Q.motor(1,[5 6]), ...
                Q.motor(2,[5 6]), ...
                Q.motor(3,[5 6]), ...
                'Color', Q.color, 'LineWidth', 2);
            Q.h_qn = text(...
                Q.motor(1,5) + Q.text_dist, ...
                Q.motor(2,5) + Q.text_dist, ...
                Q.motor(3,5) + Q.text_dist, num2str(qn));
            hold(Q.h_3d, 'off')
        end

        % Update quad state
        function UpdateQuadState(Q, state, time)
            Q.state = state;
            Q.time = time;
            Q.rot = QuatToRot(state(7:10))'; % Q.rot needs to be body-to-world
        end

        % Update desired quad state
        function UpdateDesiredQuadState(Q, des_state)
            Q.des_state = des_state;
        end

        % Update quad history
        function UpdateQuadHist(Q)
            Q.k = Q.k + 1;
            Q.time_hist(Q.k) = Q.time;
            Q.state_hist(:,Q.k) = Q.state(1:6);
            Q.state_des_hist(:,Q.k) = Q.des_state(1:6);
        end

        % Update motor position
        function UpdateMotorPos(Q)
            Q.motor = quad_pos(Q.state(1:3), Q.rot, Q.wingspan, Q.height);
        end

        % Truncate history
        function TruncateHist(Q)
            Q.time_hist = Q.time_hist(1:Q.k);
            Q.state_hist = Q.state_hist(:, 1:Q.k);
            Q.state_des_hist = Q.state_des_hist(:, 1:Q.k);
        end

        % Update quad plot
        function UpdateQuadPlot(Q, state, des_state, time)
            Q.UpdateQuadState(state, time);
            Q.UpdateDesiredQuadState(des_state);
            Q.UpdateQuadHist();
            Q.UpdateMotorPos();
            set(Q.h_m13, ...
                'XData', Q.motor(1,[1 3]), ...
                'YData', Q.motor(2,[1 3]), ...
                'ZData', Q.motor(3,[1 3]));
            set(Q.h_m24, ...
                'XData', Q.motor(1,[2 4]), ...
                'YData', Q.motor(2,[2 4]), ...
                'ZData', Q.motor(3,[2 4]));
            set(Q.h_qz, ...
                'XData', Q.motor(1,[5 6]), ...
                'YData', Q.motor(2,[5 6]), ...
                'ZData', Q.motor(3,[5 6]))
            set(Q.h_qn, 'Position', ...
                [Q.motor(1,5) + Q.text_dist, ...
                Q.motor(2,5) + Q.text_dist, ...
                Q.motor(3,5) + Q.text_dist]);
            set(Q.h_pos_hist, ...
                'XData', Q.state_hist(1,1:Q.k), ...
                'YData', Q.state_hist(2,1:Q.k), ...
                'ZData', Q.state_hist(3,1:Q.k));
            set(Q.h_pos_des_hist, ...
                'XData', Q.state_des_hist(1,1:Q.k), ...
                'YData', Q.state_des_hist(2,1:Q.k), ...
                'ZData', Q.state_des_hist(3,1:Q.k));
            drawnow;
        end
    end

end

author: yrlu

repo: quadrotor

save_path: github-repos/MATLAB/yrlu-quadrotor


text:
function accuracy = evalAccuracyHungarian(label, labelTruth)
% Find segmentation accuracy.
% Input:
%   label: label given by algorithm(e.g. SSC).
%   labelTruth: ground truth label.

labelVal = unique(label);% e.g. [0 128 255];
labelTruthVal = unique(labelTruth);% e.g.[1 2 3 4]
n = length(labelVal);% e.g. 3
nTruth = length(labelTruthVal);% e.g. 4

if n > nTruth
    accuracy = evalAccuracyHungarian(labelTruth, label);
    return;
end

minBadLabelCount = +inf;
combMat = nchoosek(labelTruthVal, n);
costMat = zeros(n);
for ii = 1:size(combMat, 1)
    combVal = combMat(ii, :);% e.g. [1 3 4]

    for jj = 1:n
        for kk = 1:n
            costMat(jj, kk) = sum( labelTruth( label == labelVal(jj) ) ~= combVal(kk) );
        end
    end
    [~, badLabelCount] = munkres(costMat);
    if badLabelCount < minBadLabelCount
        minBadLabelCount = badLabelCount;
    end
end
accuracy = 1 - minBadLabelCount / length(label);

end

repo: BUPT-ICS-Courseware


text:
function [Ax,Ay,b,K] = convexhullConvex(varargin)
% Two upper bounds from tangents
% y > f(xL) + (x-xL)*df(xL)
% y > f(xU) + (x-xL)*df(xU)
% Upper bound from conneting extreme points
% y < f(xU)(x-xL)/(xU-xL) +  f(xL)(xU-x)/(xU-xL)
% can be wrtitten as
% Ax*x + Ay*y < b

if rem(nargin,3)
    error('The convex hull generator assumes n triplets (x,f,df).')
end
m = nargin/3;
x = [varargin{(1:m)}]';
f = [varargin{(1:m)+m}]';
df = [varargin{(1:m)+2*m}]';

if all(diff(x)>0)
    Ay = [-ones(m,1);1];
    b  = [-f + x.*df; -f(end)*x(1)/(x(end)-x(1)) +  f(1)*x(end)/(x(end)-x(1))];
    Ax  = [df;-f(end)/(x(end)-x(1)) + f(1)/(x(end)-x(1))];

     % Don't use ill-conditioned cuts
    if df(1)<-1000
        Ax(1)=[];
        Ay(1) = [];
        b(1) = [];
    end
    if df(end)>1000
        Ax(end)=[];
        Ay(end) = [];
        b(end) = [];
    end

else
    Ax = [];
    Ay = [];
    b = [];
end
K.f = 0;
K.l = length(b);


text:
function [ X ] = hyperFcls( M, U )
%HYPERFCLS Performs fully constrained least squares on pixels of M.
% hyperFcls performs fully constrained least squares of each pixel in M
% using the endmember signatures of U.  Fully constrained least squares
% is least squares with the abundance sum-to-one constraint (ASC) and the
% abundance nonnegative constraint (ANC).
%
% Usage
%   [ X ] = hyperFcls( M, U )
% Inputs
%   M - HSI data matrix (p x N)
%   U - Matrix of endmembers (p x q)
% Outputs
%   X - Abundance maps (q x N)
%
% References
%   "Fully Constrained Least-Squares Based Linear Unmixing." Daniel Heinz,
% Chein-I Chang, and Mark L.G. Althouse. IEEE. 1999.

if (ndims(U) ~= 2)
    error('M must be a p x q matrix.');
end

[p1, N] = size(M);
[p2, q] = size(U);
if (p1 ~= p2)
    error('M and U must have the same number of spectral bands.');
end

p = p1;
X = zeros(q, N);
Mbckp = U;
for n1 = 1:N
    count = q;
    done = 0;
    ref = 1:q;
    r = M(:, n1);
    U = Mbckp;
    while not(done)
        als_hat = inv(U.'*U)*U.'*r;
        s = inv(U.'*U)*ones(count, 1);

        % IEEE Magazine method (http://www.planetary.brown.edu/pdfs/3096.pdf)
        % Contains correction to sign.  Error in original paper.
        afcls_hat = als_hat - inv(U.'*U)*ones(count, 1)*inv(ones(1, count)*inv(U.'*U)*ones(count, 1))*(ones(1, count)*als_hat-1);

        % See if all components are positive.  If so, then stop.
        if (sum(afcls_hat>0) == count)
            alpha = zeros(q, 1);
            alpha(ref) = afcls_hat;
            break;
        end
        % Multiply negative elements by their counterpart in the s vector.
        % Find largest abs(a_ij, s_ij) and remove entry from alpha.
        idx = find(afcls_hat<0);
        afcls_hat(idx) = afcls_hat(idx) ./ s(idx);
        [val, maxIdx] = max(abs(afcls_hat(idx)));
        maxIdx = idx(maxIdx);
        alpha(maxIdx) = 0;
        keep = setdiff(1:size(U, 2), maxIdx);
        U = U(:, keep);
        count = count - 1;
        ref = ref(keep);
    end
    X(:, n1) = alpha;
end

return;

author: davidkun

repo: HyperSpectralToolbox


text:
function X=ctranspose(X)
%CTRANSPOSE (overloaded)

if isa(X,'blkvar')
    X = sdpvar(X);
end

n = X.dim(1);
m = X.dim(2);
ind = reshape(reshape(1:n*m,n,m)',n*m,1);
if isreal(X.basis)
    X.basis = X.basis(ind,:);
else
   X.basis = conj(X.basis(ind,:));
end
X.dim(1) = m;
X.dim(2) = n;
% Reset info about conic terms
X.conicinfo = [0 0];

% Flip noncommuting terms
% Get all the tables, and expand them so that they correspond to the same
% number of variables globally (nonCommutingTable is not up to date after a
% new commuting variables has been defined, to save flops)
nonCommutingTable         = yalmip('nonCommutingTable');
[monomtable,variabletype] = yalmip('monomtable');
if size(monomtable,1)>size(nonCommutingTable,1)
    nonCommutingTable((1+size(nonCommutingTable,1)):(size(monomtable,1)),1) = (1+size(nonCommutingTable,1)):(size(monomtable,1));
end
% Cast commutative variables as nc temporarily by adding them to the table
commuting = find(~any(nonCommutingTable,2));
nonCommutingTable(commuting,1) = commuting;

for i = 1:length(X.lmi_variables)
    if nonCommutingTable(X.lmi_variables(i),2)
        monoms = nonCommutingTable(X.lmi_variables(i),2:end);
        monoms(find(monoms)) = fliplr( monoms(find(monoms)));
        monoms = [nonCommutingTable(X.lmi_variables(i),1) monoms];
        old = findrows(nonCommutingTable,monoms);
        old = findrows(nonCommutingTable(:,2:end),monoms(2:end));
        if ~isempty(old)
            old = old(find((monoms(1) == nonCommutingTable(old,1)) | isnan(monoms(1)) & isnan(nonCommutingTable(old,1))));
        end
        if isempty(old)
            % Create a new monomial
            monomtable(end+1,end+1) = 0;
            variabletype(end+1) = variabletype(X.lmi_variables(i));
            nonCommutingTable = [nonCommutingTable;monoms];
            X.lmi_variables(i) = size(nonCommutingTable,1);
        else
            X.lmi_variables(i) = old;
        end
    end
end
% Fucked up order (lmi_variables should be sorted and unique)
if any(diff(X.lmi_variables)<0)
    [i,j]=sort(X.lmi_variables);
    X.basis = [X.basis(:,1) X.basis(:,j+1)];
    X.lmi_variables = X.lmi_variables(j);
end
[un_Z_vars2] = uniquestripped(X.lmi_variables);
if length(un_Z_vars2) < length(X.lmi_variables)
    [un_Z_vars,hh,jj] = unique(X.lmi_variables);
    if length(X.lmi_variables) ~=length(un_X_vars)
        X.basis = Z.basis*sparse([1 1+jj],[1 1+(1:length(jj))],ones(1,1+length(jj)))';
        X.lmi_variables = un_Z_vars;
    end
end
if size(monomtable,2) < size(monomtable,1)
    monomtable(size(monomtable,1),size(monomtable,1)) = 0;
end
yalmip('nonCommutingTable',nonCommutingTable);
yalmip('setmonomtable',monomtable,variabletype);


text:
function printLevelSetForMMG(fileName,ls)
fid = fopen(fileName,'w');

%HEADER
fprintf(fid,'%s\n\n','MeshVersionFormatted 2');
fprintf(fid,'%s\n\n',['Dimension ',num2str(dim)]);

% Write nodes
fprintf(fid,'SolAtVertices\n');
fprintf(fid,'[1 1]\n');
fprintf(fid,'\n');

N = size(ls,1);
fprintf(fid,'%d\n',N);
fprintf(fid,'%f %d\n',[ls(:,1),zeros(N,1)]');
fprintf(fid,'\n');

fprintf(fid,'End \n');
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


function refinement_model=train_refiner(imnames, region_meta_info, featdir, sptextdir, regspimgdir, sbddir, Wsz, categid)
feats=[];
labels=[];

regidsall=cell(numel(imnames),1);
movindall=cell(numel(imnames),1);
%collect training examples
for i=1:numel(imnames)
    %check if there is any groundtruth for this category
    if(~any(region_meta_info.gt{i}==categid))
        continue;
    end
    fprintf('Getting training regions %d\n', i);

    %find ground truths
    gtids=find(region_meta_info.gt{i}==categid);

    %find regions that overlap by more than 70\%
    [mov, movind]=max(region_meta_info.overlaps{i}(gtids,:),[],1);
    regids=find(mov>=0.7);
    movind=gtids(movind(regids));

    %index into features/sp
    nongtids=find(region_meta_info.gt{i}==0);
    regids=nongtids(regids);

    regidsall{i}=regids;
    movindall{i}=movind;
end
for i=1:numel(imnames)
    regids=regidsall{i};
    movind=movindall{i};
    if(isempty(regids))
        continue;
    end
    fprintf('Getting training features:%d\n', i);

    %load gt
    [cls, inst] = load_gt(sbddir, imnames{i});
    %load features
    d=load(fullfile(featdir, [imnames{i} '.mat']));
    f=d.feats(regids,:)';

    %load superpixels
    [sp, reg2sp]=read_sprep(fullfile(sptextdir, [imnames{i} '.txt']), fullfile(regspimgdir, [imnames{i} '.png']));

    %get boxes
    boxes=get_region_boxes(sp, reg2sp);
    spf=[];
    lab=[];
    %for each selected region
    for j=1:numel(regids)
        box=boxes(regids(j),:);
        box=expand_box(box, 16/(227-32));
        msk=reg2sp(:,regids(j));
        msk=msk(sp);
        clipped_msk=clip_img(msk, box);
        clipped_gt=clip_img(double(double(inst)==movind(j)),box);
        cl_msk_rsz=imresize(clipped_msk,Wsz);
        cl_gt_rsz=imresize(clipped_gt, Wsz);
        spf=[spf cl_msk_rsz(:)];
        lab=[lab cl_gt_rsz(:)>0.5];
    end

    f=[f; spf];
    if(isempty(feats))
        feats=zeros(size(f,1),10000);
        labels=zeros(10000,prod(Wsz));
        cnt=0;
    end
    feats(:,cnt+(1:size(f,2)))=f;
    labels(cnt+(1:size(f,2)),:)=double(lab');
    cnt=cnt+size(f,2);
end
feats=feats(:,1:cnt);
labels=labels(1:cnt,:);

%train svms
for i=1:prod(Wsz)
fprintf('Training coarse model for grid cell %d\n',i);
logisticmodel=liblinear_train(double(labels(:,i)), double(feats), '-s 0 -c 0.01 -B 1', 'col');
if(logisticmodel.Label(1)==0)
    W{i}=-logisticmodel.w;
else
    W{i}=logisticmodel.w;
end
end
W=cat(1, W{:});

%make all predictions
scr=bsxfun(@plus, W(:,1:end-1)*feats, W(:,end));
pred=1./(1+exp(-scr));

%train superpixel predictor
feats=zeros(2,1000000);
labels=zeros(1, 1000000);;
cnt=0;
cnt2=0;
for i=1:numel(imnames)
    regids=regidsall{i};
    movind=movindall{i};
    if(isempty(regids)) continue; end
    fprintf('Getting second stage feats:%d\n',i);
    %load superpixels
    [sp, reg2sp]=read_sprep(fullfile(sptextdir, [imnames{i} '.txt']), fullfile(regspimgdir, [imnames{i} '.png']));

    %get all coarse predictions
    predcurr=pred(:,cnt+(1:numel(regids)));
    cnt=cnt+numel(regids);

    boxes=get_region_boxes(sp, reg2sp);

    spareas=accumarray(sp(:),1);
    %apply it to each region
    for j=1:numel(regids)
        box=boxes(regids(j),:);
        box=expand_box(box, 16/(227-32));
        coarse_reg2sp=apply_mask(predcurr(:,j), box, sp,Wsz);
        orig_reg2sp=reg2sp(:,regids(j));

        %only choose large pixels and pixels for which prediction is large enough
        idx=find(spareas>100 & (coarse_reg2sp>0.2 | orig_reg2sp>0));

        %next get the true labels
        lab=reg2sp(:,movind(j));
        feats(1,cnt2+(1:numel(idx)))=coarse_reg2sp(idx);
        feats(2, cnt2+(1:numel(idx)))=orig_reg2sp(idx);
        labels(cnt2+(1:numel(idx)))=lab(idx);
        cnt2=cnt2+numel(idx);
    end
end
feats=feats(:,1:cnt2);
labels=labels(:,1:cnt2);

%learn second model
spmodel=liblinear_train(labels(:), feats, '-s 0 -c 1 -B 1', 'col');

refinement_model.W=W;
refinement_model.spmodel=spmodel;


text:
function order = circle_xy_size ( rule )

%*****************************************************************************80
%
%% CIRCLE_XY_SIZE sizes an XY quadrature rule inside the unit circle in 2D.
%
%  Integration region:
%
%    X*X + Y*Y <= 1.0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 April 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Frank Lether,
%    A Generalized Product Rule for the Circle,
%    SIAM Journal on Numerical Analysis,
%    Volume 8, Number 2, June 1971, pages 249-253.
%
%    Arthur Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971,
%    ISBN: 0130438936,
%    LC: QA311.S85.
%
%  Parameters:
%
%    Input, integer RULE, the rule desired.
%      1, 1 point 1-st degree;
%      2, 4 point 3-rd degree, Stroud S2:3-1;
%      3, 4 point 3-rd degree, Lether #1;
%      4, 4 point 3-rd degree, Stroud S2:3-2;
%      5, 5 point 3-rd degree;
%      6, 7 point 5-th degree;
%      7, 9 point 5-th degree;
%      8, 9 point 5-th degree, Lether #2;
%      9, 12 point 7-th degree;
%     10, 16 point 7-th degree, Lether #3;
%     11, 21 point 9-th degree, Stroud S2:9-3;
%     12, 25 point 9-th degree, Lether #4 (after correcting error);
%     13, 64 point 15-th degree Gauss product rule.
%
%    Output, integer ORDER, the order of the desired rule.
%
  if ( rule == 1 )

    order = 1;

  elseif ( rule == 2 )

    order = 4;

  elseif ( rule == 3 )

    order = 4;

  elseif ( rule == 4 )

    order = 4;

  elseif ( rule == 5 )

    order = 5;

  elseif ( rule == 6 )

    order = 7;

  elseif ( rule == 7 )

    order = 9;

  elseif ( rule == 8 )

    order = 9;

  elseif ( rule == 9 )

    order = 12;

  elseif ( rule == 10 )

    order = 16;

  elseif ( rule == 11 )

    order = 21;

  elseif ( rule == 12 )

    order = 25;

  elseif ( rule == 13 )

    order = 64;

  else

    order = -1;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CIRCLE_XY_SIZE - Fatal error!\n' );
    fprintf ( 1, '  There is no rule of index %d\n', rule );
    error ( 'CIRCLE_XY_SIZE - Fatal error!\n' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function s = spider_path(set)
% This function returns the path where spider is installed.
% If an arugment is provided than it adds all spider subdirectory paths
% to the Matlab search path.

subdirs=spider_subdirs;

% if nargin==0 %% can just return s and do not set paths
if nargin>0
    for i=1:length(subdirs)
        t = cellstr(subdirs{i});
        path(path,t{:});
    end
s=subdirs{1};

end

s=subdirs{1};
k = findstr(s,'spider');
s=s(1:k+5);

if nargin > 0
    path(path,s);
else
    subdirs=spider_subdirs;
% if nargin
    addpath(subdirs{:})
    addpath(s);

end

% show all default values for algorithms e.g C=Inf in svm, otherwise suppress
global display_tree_show_defaults;
display_tree_show_defaults=0;

% doesn't draw whole tree, only up to depth given
global display_tree_depth;
display_tree_depth=100;

% displays indexes e.g [1 3 4] if set to 1,
% displays last element of array e.g 4 of [1 3 4] if set to 2,
% or else displays a single integer index system for whole tree
%              (breadth-wise, leftmost is smallest value)
global display_tree_array_indexing;
display_tree_array_indexing=2;

% place index display at front or end of each line of output
global display_tree_index_at_front;
display_tree_index_at_front=1;

% draw brackets in display for object sets
global display_tree_show_brackets;
display_tree_show_brackets=1;

% number of characters of tab spacing for nested objects in tree
global display_tree_tab_step;
display_tree_tab_step=3;

% whether to allow recursive search for variables in statements
% such as:  b=algs({svm svm}); b.C=2
% (what you lose is it is much slower...)

global recursive_subsasgn_off;
recursive_subsasgn_off=0;
% else

if isunix
     s = [ s '/'];
 else
     s = [ s '\'];
 end

s(s=='\')='/';

%% Take into account weka.jar
wekajar=[matlabroot filesep 'java' filesep 'jar' filesep 'weka.jar'];
if( length(dir(wekajar))>0)
    javaaddpath(wekajar);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function table2 = i4mat_border_cut ( m, n, table )

%*****************************************************************************80
%
%% I4MAT_BORDER_CUT cuts the "border" of an I4MAT.
%
%  Discussion:
%
%    We suppose the input data gives values of a quantity on nodes
%    on a 2D grid, and we wish to create a new table corresponding only
%    to those nodes in the interior of the 2D grid.
%
%      0 0 0 0 0 0
%      0 * * * * 0    * * * *
%      0 * * * * 0 -> * * * *
%      0 * * * * 0    * * * *
%      0 0 0 0 0 0
%
%    The illustration suggests the situation in which a 5 by 6 array
%    is input, and a 3 by 4 array is to be output.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the table data.
%
%    Output, integer TABLE2(M-2,N-2), the new table data.
%
  if ( m <= 2 || n <= 2 )
    table2 = [];
    return
  end

  table2(1:m-2,1:n-2) = round ( table(2:m-1,2:n-1) );

  return
end


text:
function [Pt,a,b] = mci_lds_lat2par (P,M)
% Convert latent params to params
% FORMAT [Pt,a,b] = mci_lds_lat2par (P,M)
%
% P         Parameters (latent)
% M         model structure
%
% Pt        Parameters (transformed)
% a         diagonal values
% b         off-diagonal values
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: mci_lds_lat2par.m 6548 2015-09-11 12:39:47Z will $

% Diagonal entries
s=exp(P(1:M.d));
a=s*M.a_typical;

% Off-diagonal entries
b=P(M.d+1:end);

Pt=[a(:);b(:)];

author: spm


text:
function obj = updateSurfaceplot(obj, surfaceIndex)

%-AXIS INDEX-%
axIndex = obj.getAxisIndex(obj.State.Plot(surfaceIndex).AssociatedAxis);

%-CHECK FOR MULTIPLE AXES-%
[xsource, ysource] = findSourceAxis(obj,axIndex);

%-SURFACE DATA STRUCTURE- %
image_data = get(obj.State.Plot(surfaceIndex).Handle);

%-AXIS DATA-%
eval(['xaxis = obj.layout.xaxis' num2str(xsource) ';']);
eval(['yaxis = obj.layout.yaxis' num2str(ysource) ';']);

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

%-surface xaxis-%
obj.data{surfaceIndex}.xaxis = ['x' num2str(xsource)];

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

%-surface yaxis-%
obj.data{surfaceIndex}.yaxis = ['y' num2str(ysource)];

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

% check for 3D
if any(nonzeros(image_data.ZData))

    %-surface type-%
    obj.data{surfaceIndex}.type = 'surface';

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

    %-surface x-%
    obj.data{surfaceIndex}.x = image_data.XData;

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

    %-surface y-%
    obj.data{surfaceIndex}.y = image_data.YData;

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

    %-surface z-%
    obj.data{surfaceIndex}.z = image_data.ZData;

else

    %-surface type-%
    obj = updateImage(obj, surfaceIndex);

    %-surface x-%
    obj.data{surfaceIndex}.x = image_data.XData(1,:);

    %-surface y-%
    obj.data{surfaceIndex}.y = image_data.YData(:,1);
end

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

%-surface name-%
obj.data{surfaceIndex}.name = image_data.DisplayName;

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

%-surface showscale-%
obj.data{surfaceIndex}.showscale = false;

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

%-surface visible-%
obj.data{surfaceIndex}.visible = strcmp(image_data.Visible,'on');

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

leg = get(image_data.Annotation);
legInfo = get(leg.LegendInformation);

switch legInfo.IconDisplayStyle
    case 'on'
        showleg = true;
    case 'off'
        showleg = false;
end

obj.data{surfaceIndex}.showlegend = showleg;

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

end

author: goGPS-Project


text:
% [INPUT]
% f = A valid figure instance.
% t = A string representing the figure title.
%
% [OUTPUT]
% t = The handler of the figure title.

function t = figure_title(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
		ip.addRequired('f',@(x)validateattributes(x,{'matlab.ui.Figure'},{'scalar'}));
        ip.addRequired('s',@(x)validateattributes(x,{'char'},{'nonempty' 'size' [1 NaN]}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    f = ipr.f;
    s = ipr.s;

    nargoutchk(0,1);

    if (nargout == 0)
        figure_title_internal(f,s);
    else
        t = figure_title_internal(f,s);
    end

end

function t = figure_title_internal(f,s)

    f_font_size = get(f,'DefaultAxesFontSize') + 4;
    f_units = get(f,'Units');

    if (~strcmp(f_units,'pixels'))
        set(f,'Units','pixels');
        f_position = get(f,'Position');
        set(f,'Units',f_units);
    else
        f_position = get(f,'Position');
    end

    adjustment = ((f_font_size - 4) * 6.35) / f_position(4);
    title = NaN;
    y_max = 0;
    y_min = 1;

    handles = findobj(f,'Type','axes');
    handles_len = length(handles);
    handles_position = zeros(handles_len,4);

    for i = 1:handles_len
        h_current = handles(i);

        f_position = get(h_current,'Position');
        handles_position(i,:) = f_position;

        if (~strcmp(get(h_current,'Tag'),'FigureTitle'))
            f_y = f_position(2);
            f_height = f_position(4);

            if (f_y < y_min)
                y_min = f_y - (adjustment / 15);
            end

            if ((f_height + f_y) > y_max)
                y_max = f_height + f_y + (adjustment / 10);
            end
        else
            title = h_current;
        end
    end

    if (y_max > 0.92)
        scale = (0.92 - y_min) / (y_max - y_min);

        for i = 1:handles_len
            f_position = handles_position(i,:);
            f_position(2) = ((f_position(2) - y_min) * scale) + y_min;
            f_position(4) = (f_position(4) * scale) - ((1 - scale) * (adjustment / 15));

            set(handles(i),'Position',f_position);
        end
    end

    next_plot = get(f,'NextPlot');
    set(f,'NextPlot','add');

    if (ishghandle(title))
        delete(title);
    end

    axes('Position',[0 1 1 1],'Tag','FigureTitle','Visible','off');
    t_internal = text(0.5000,-0.0157,s,'FontSize',f_font_size,'HorizontalAlignment','center');

    set(f,'NextPlot',next_plot);

    axes(gca());

    if (nargout == 1)
        t = t_internal;
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
%SerialLink.itorque Inertia torque
%
% TAUI = R.itorque(Q, QDD) is the inertia force/torque vector (1xN) at the
% specified joint configuration Q (1xN) and acceleration QDD (1xN), and N
% is the number of robot joints. TAUI = INERTIA(Q)*QDD.
%
% If Q and QDD are matrices (KxN), each row is interpretted as a joint state
% vector, and the result is a matrix (KxN) where each row is the corresponding
% joint torques.
%
% Note::
% - If the robot model contains non-zero motor inertia then this will
%   included in the result.
%
% See also SerialLink.inertia, SerialLink.rne.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function it = itorque(robot, q, qdd)
	it = rne(robot, q, zeros(size(q)), qdd, [0;0;0]);

author: Allopart

repo: rbpf-gmapping


text:
function [kHz] = Hz2kHz(Hz)
% Convert frequency from hertz to kilohertz.
% Chad A. Greene 2012
kHz = Hz*1e-3;

author: Sable

repo: mcbench-benchmarks


text:
function k = plotDataTimeSeries(y, iS, iD, iDt, option)
if ~exist('option','var')
    option = '';
end
ml = y.GetMeasurementList('matrix');
k = find(ml(:,1)==iS & ml(:,2)==iD & ml(:,4)==iDt);
if ishandle(1)
    if ~strcmp(option,'keep') && ~strcmp(option,'keepopen') && ~strcmp(option,'open')
        close(1);
    end
end
figure
plot(y.dataTimeSeries(:,k))

author: BUNPC


text:
clear all, close all, clc
A = imread('../../CH01_SVD/DATA/dog.jpg');
B = rgb2gray(A);

%% Wavelet Compression
figure
[C,S] = wavedec2(B,4,'db1');

% Zero out all small coefficients and inverse transform
Csort = sort(abs(C(:)));
keepvec = [.1 .05 .01 .005 .002 .001];
for k=1:6
    figure
    keep = keepvec(k);
    % keep = 0.05;
    thresh = Csort(floor((1-keep)*length(Csort)));
    ind = abs(C)>thresh;
    Cfilt = C.*ind;

    % Plot Reconstruction
    Alow=uint8(waverec2(Cfilt,S,'db1'));
    imagesc(uint8(Alow))
    colormap gray

    colormap gray
    axis off
    axis tight
    set(gcf,'Position',[100 100 600 800])
    set(gcf,'PaperPositionMode','auto')
    %     print('-depsc2', '-loose', ['../figures/2DFFT_Compress',num2str(k)]);
        print('-dpng', '-loose', ['../figures/2DWavelet_Compress',num2str(k)]);

end

author: dynamicslab


text:
function geometry_test020 ( )

%*****************************************************************************80
%
%% TEST020 tests CUBE_SIZE_3D, CUBE_SHAPE_3D, SHAPE_PRINT_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST020\n' );
  fprintf ( 1, '  For the cube,\n' );
  fprintf ( 1, '  CUBE_SIZE_3D returns dimension information;\n' );
  fprintf ( 1, '  CUBE_SHAPE_3D returns face and order information.\n' );
  fprintf ( 1, '  SHAPE_PRINT_3D prints this information.\n' );

  [ point_num, edge_num, face_num, face_order_max ] = cube_size_3d ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of points =   %d\n', point_num );
  fprintf ( 1, '  Number of edges =    %d\n', edge_num );
  fprintf ( 1, '  Number of faces =    %d\n', face_num );
  fprintf ( 1, '  Maximum face order = %d\n', face_order_max );

  [ point_coord, face_order, face_point ] = cube_shape_3d ( ...
    point_num, face_num, face_order_max );

  shape_print_3d ( point_num, face_num, face_order_max, ...
    point_coord, face_order, face_point );

  return
end


text:
function legendre_polynomial_test04 ( )

%*****************************************************************************80
%
%% LEGENDRE_POLYNOMIAL_TEST04 tests P_QUADRATURE_RULE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 March 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LEGENDRE_POLYNOMIAL_TEST04:\n' );
  fprintf ( 1, '  P_QUADRATURE_RULE computes the quadrature rule\n' );
  fprintf ( 1, '  associated with P(n,x)\n' );

  n = 5;
  [ x, w ] = p_quadrature_rule ( n );

  r8vec2_print ( n, x, w,  '      X            W' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Use the quadrature rule to estimate:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Q = Integral ( -1 <= X < +1 ) X^E dx\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   E       Q_Estimate      Q_Exact\n' );
  fprintf ( 1, '\n' );

  for e = 0 : 2 * n - 1
    if ( e == 0 )
      f = ones ( n, 1 );
    else
      f = x .^ e;
    end
    q = w' * f;
    q_exact = p_integral ( e );
    fprintf ( 1, '  %2d  %14g  %14g\n', e, q, q_exact );
  end

  return
end


text:
addpath('./code/') ;
mydir = {'train_1','train_2','train_3','test_1','test_2','test_3','test_1_new','test_2_new','test_3_new'};
srateNew = 256;
windL = 30; % window length in s
windS = 30; % window shift in s

for i1 = 1: size(mydir,2)
    names = dir(['../data/' mydir{i1}]);
    for i2=3:size(names,1)
        fprintf('\n%d-%d',i1,i2-2);
        load(['../data/' mydir{i1} '/' names(i2).name]);
        dataStruct.data = dataStruct.data';

        % Highpass Filter coefficient
        Fchp = 0.5;
        alpha = (2*pi*Fchp)/double(dataStruct.iEEGsamplingRate);

        % 60Hz notch filter
        fs = double(dataStruct.iEEGsamplingRate);               %#sampling rate
        f0 = 60;                %#notch frequency
        fn = fs/2;              %#Nyquist frequency
        freqRatio = f0/fn;      %#ratio of notch freq. to Nyquist freq.
        notchWidth = 0.05;       %#width of the notch
        %Compute zeros
        myzeros = [exp( sqrt(-1)*pi*freqRatio ), exp( -sqrt(-1)*pi*freqRatio )];
        %Compute poles
        mypoles = (1-notchWidth) * myzeros;
        b = poly( myzeros ); %# Get moving average filter coefficients
        a = poly( mypoles ); %# Get autoregressive filter coefficients

        for ii=1:size(dataStruct.data,1)
            dataStruct.data(ii,:) = dataStruct.data(ii,:) - mean(dataStruct.data(ii,:)); % centering
            LPOutput = filter(b,a,dataStruct.data(ii,:));
            HPOutput = zeros(size(LPOutput));
            HPOutput(1) = LPOutput(1);
            for iIndex = 2:length(LPOutput)
                HPOutput(iIndex) = (1 - alpha)*HPOutput(iIndex-1) + LPOutput(iIndex) - LPOutput(iIndex-1);
            end
            dataStruct.data(ii,:) = HPOutput;
        end
        dataStruct.data = double(dataStruct.data(:,dataStruct.iEEGsamplingRate*4:end-dataStruct.iEEGsamplingRate*4)); % removing 4 first and last seconds due to filtering
        dataStruct.data = (resample(dataStruct.data',srateNew,double(dataStruct.iEEGsamplingRate)))';
        feat = NaN(size(dataStruct.data,1),fix(((length(dataStruct.data(1,:)))-windL*srateNew+windS*srateNew)/(windS*srateNew)),102);
        for ii=1:size(dataStruct.data,1)
            fprintf('.');
            dataCh = enframe(dataStruct.data(ii,:),windL*srateNew,windS*srateNew);
            for iii=1:size(dataCh,1)
                if range(dataCh(iii,:))<1000
                    feat(ii,iii,:) = SPC_extract_features(dataCh(iii,:)');
                end
            end
        end
        %seq = dataStruct.sequence;
        if exist(['./feat/' mydir{i1} '/']) ~= 7
            mkdir(['./feat/' mydir{i1} '/']);
        end
        save(['./feat/' mydir{i1} '/' names(i2).name(1:end-4) '_feat'],'feat');
    end
end

author: alexandrebarachant


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD 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.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

function [pX,pEx,bbP0] = tldGeneratePositiveData(tld,overlap,im0,p_par)

pX   = [];
pEx  = [];%zeros(prod(tld.patchsize),numWarps);

% Get closest bbox
[~,idxP] = max(overlap);
bbP0 =  tld.grid(1:4,idxP);

% Get overlapping bboxes
idxP = find(overlap > 0.6);
if length(idxP) > p_par.num_closest
    [~,sIdx] = sort(overlap(idxP),'descend');
    idxP = idxP(sIdx(1:p_par.num_closest));
end
bbP  = tld.grid(:,idxP);
if isempty(bbP), return; end

% Get hull
bbH  = bb_hull(bbP);
cols = bbH(1):bbH(3);
rows = bbH(2):bbH(4);

im1 = im0;
pEx = tldGetPattern(im1,bbP0,tld.model.patchsize);
if tld.model.fliplr
pEx = [pEx tldGetPattern(im1,bbP0,tld.model.patchsize,1)];
end
for i = 1:p_par.num_warps
    if i > 1
        randomize = rand; % Sets the internal randomizer to the same state
        %patch_input = img_patch(im0.input,bbH,randomize,p_par);
        patch_blur = img_patch(im0.blur,bbH,randomize,p_par);
        im1.blur(rows,cols) = patch_blur;
        %im1.input(rows,cols) = patch_input;
    end

    % Measures on blured image
    pX  = [pX fern(5,im1,idxP,0)];

    % Measures on input image
    %pEx(:,i) = tldGetPattern(im1,bbP0,tld.model.patchsize);
    %pEx = [pEx tldGetPattern(im1,tld.grid(1:4,idxP),tld.model.patchsize)];

end

repo: OpenTLD


text:
function mr = mrReadSSFP(pth, phaseCyclePath, nFrames, varargin);
%
%  mr = mrReadSSFP([non-phase-cycle directory], [phase-cycle directory], [nFrames], [options]);
%
% Read in a series of SSFP files, combining phase-cycle and non-phase cycle
% scans using the maximum intensity projection described in (Lee et al,
% Neuroimage, 2008).
%
% Note that this function reconstructs the P-files from K-space. To get the
% raw K-space data, use the function mrReadSSFPRaw.
%
%
% ras, 02/2009.
if notDefined('pth')
    pth = mrvSelectFile('r', '7', 'Select SSFP Raw (P*.7) file');
end

if notDefined('phaseCyclePath'), phaseCyclePath = '';   end

if notDefined('nFrames'),       nFrames = [];           end

%% params
combineMethod = 'sumofsquares';  % 'maxintensity' or 'sumofsquares'
kspaceFraction = 5/8;

% GUM: defaults are from SSFP pilot data
voxelSize = [1.5 1.5 2.5 2.16];

%% parse options
for ii = 1:2:length(varargin)
    eval( sprintf('%s = %s', varargin{ii}, num2str(varargin{ii+1})) );
end

%% get full path, and fileparts of the path:
pth = fullpath(pth);
[p f ext] = fileparts(pth);

%% initalize an empty mr struct:
mr = mrCreateEmpty;
mr.format = 'ssfp';
mr.name = [f ext];
mr.path = pth;

% set other fields
mr.voxelSize = voxelSize;
mr.dataUnits = 'Arbitrary';
mr.dimUnits = {'mm' 'mm' 'mm' 'sec'};

%% read the data
[R I] = ssfpReadData(pth, nFrames);

mr.data = I;
mr.dims = size(mr.data);
if length(mr.dims) < 4, mr.dims(4) = 1; end
mr.extent = mr.dims .* mr.voxelSize;
mr.dataRange = mrvMinmax(mr.data);

%% set other fields: coordinat spaces, etc.
mr.spaces = mrStandardSpaces(mr);

%% combine with phase-cycled data if needed
if ~isempty(phaseCyclePath)
    mr_cycle = mrReadSSFP(phaseCyclePath, '', nFrames);

    % note what data we're combining, and how:
    msg = sprintf('Combination of %s and %s using %s method.', mr.path, ...
                  mr_cycle.path, combineMethod);
    mr.comments = strvcat(mr.comments, msg);

    switch lower(combineMethod)
        case {'maxintensity' 'max' 'mip' 'maxintensityprojection'}
            mr = maxIntensityProjection(mr, mr_cycle);
        case {'sumofsquares' 'sos' 'sumsquared'}
            mr = sumOfSquares(mr, mr_cycle);
        otherwise
            error('Invalid combination method: %s.', combineMethod)
    end
    return
end

return
% /--------------------------------------------------------------------/ %

% /--------------------------------------------------------------------/ %
function mr = maxIntensityProjection(mr, mr2);
% Combine two mr data volumes using the maximum-intensity projection method
% described in Lee et al, NeuroImage, 2008). This puts the data in the
% first mr struct provided.
%
% ras 02/2009.

verbose = prefsVerboseCheck;
if verbose >= 1
    h = mrvWaitbar(0, 'Computing Max Intensity Projection');
end

% first, compute the mean images over time for each MR data set.
mu1 = nanmean(abs(mr.data), 4);
mu2 = nanmean(abs(mr2.data), 4);

% create a binary mask indicating which MR data set from which draw data
% for each voxel.
mask = zeros( size(mu1) );
for ii = 1:numel(mask)
    if mu1(ii) > mu2(ii)
        mask(ii) = 1;
    else
        mask(ii) = 2;
    end
end

% now, for those points for which the mask indicates to use mr2, copy over
% those data points into the mr1 data. (For the other points, where the
% mask indicates mr data 1, we already have the appropriate data loaded.
for t = 1:size(mr.data, 4)
    subvol = mr.data(:,:,:,t);
    subvol2 = mr2.data(:,:,:,t);

    subvol(mask==2) = subvol2(mask==2);

    mr.data(:,:,:,t) = subvol;

    if verbose >= 1, mrvWaitbar( t / size(mr.data, 4), h ); end
end

if verbose >= 1
    close(h);
end

return
% /--------------------------------------------------------------------/ %

% /--------------------------------------------------------------------/ %
function mr = sumOfSquares(mr, mr2);
% Combine two mr data volumes using the sum-of-squares method
% described in Lee et al, NeuroImage, 2008). This puts the data in the
% first mr struct provided.
%
% ras 02/2009.
mr.data = sqrt( (mr.data .^ 2) + (mr2.data .^ 2) );
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function bc = specific_bc(xbd,ybd)
%solution4_bc   Reference problem 1.4  boundary condition
%   bc = specific_bc(xbd,ybd);
%   input
%          xbd          x boundary coordinate vector
%          ybd          y boundary coordinate vector
%   IFISS function: DJS; 28 February 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
bc=(xbd.^2+ybd.^2).^(1/3) .*sin((pi+2*atan2(ybd,xbd))/3);
return


text: % labels01 = (labelsPM+1)/2; % maps -1->0, +1->1

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%   TTeMPS Toolbox.
%   Michael Steinlechner, 2013-2016
%   Questions and contact: michael.steinlechner@epfl.ch
%   BSD 2-clause license, see LICENSE.txt
function [X, residuum, cost, times] = alsLinsolve_fast( L, F, X, opts )

t_start = tic();
% set default opts
if ~exist( 'opts', 'var');       opts = struct();       end
if ~isfield( opts, 'nSweeps');   opts.nSweeps = 4;      end
if ~isfield( opts, 'solver');    opts.solver = 'pcg';   end

d = X.order;
n = X.size;

normF = norm(F);
g = apply(L, X) - F;
cost = cost_function_res( X, g );
residuum = norm( g ) / normF;
times = toc(t_start);

X = orthogonalize(X, 1);
for sweep = 1:opts.nSweeps
    % ====================================================================
    % LEFT-TO-RIGHT SWEEP
    % ====================================================================
    disp( ['STARTING SWEEP ', num2str(sweep), ' from left to right'] )
    disp( '===========================================================')
    for idx = 1:d-1
        disp( ['Current core: ', num2str(idx)] )

        Fi = contract( X, F, idx );
        sz = [X.rank(idx), X.size(idx), X.rank(idx+1)];

        if strcmpi( opts.solver, 'direct' )
            % if system very small
            Li = contract( L, X, idx );
            Ui = Li \ Fi(:);
            X.U{idx} = reshape( Ui, sz );

        elseif strcmpi( opts.solver, 'pcg' )

            [left, right] = Afun_prepare( L, X, idx );
            [B2, V, E] =  prepare_precond( L.L0, X, idx );

            Ui = pcg( @(y) Afun( L, y, idx, sz, left, right), ...
                     Fi(:), ...
                     1e-10, 1000, ...
                     @(y) apply_precond( B2, V, E, y, sz ), [],...
                     X.U{idx}(:) );

            X.U{idx} = reshape( Ui, sz );

        elseif strcmpi( opts.solver, 'diag' )
            X.U{idx} = solve_inner( L.L0, X, Fi, idx );

        else
            error( 'Unknown opts.solver type. Use either ''direct'', ''pcg'' (default) or ''diag''.' )
        end

        X = orth_at( X, idx, 'left', true );

        g = apply(L, X) - F;
        residuum = [residuum; norm( g ) / normF];
        cost = [cost; cost_function_res( X, g )];
        times = [times; toc(t_start)];
    end

    % ====================================================================
    % RIGHT-TO-LEFT
    % ====================================================================
    disp( 'Starting right-to-left half-sweep:')
    for idx = d:-1:2
        disp( ['Current core: ', num2str(idx)] )

        Fi = contract( X, F, idx );
        sz = [X.rank(idx), X.size(idx), X.rank(idx+1)];

        if strcmpi( opts.solver, 'direct' )
            % if system very small
            Li = contract( L, X, idx );
            Ui = Li \ Fi(:);
            X.U{idx} = reshape( Ui, sz );

        elseif strcmpi( opts.solver, 'pcg' )

            [left, right] = Afun_prepare( L, X, idx );
            [B2, V, E] =  prepare_precond( L.L0, X, idx );

            Ui = pcg( @(y) Afun( L, y, idx, sz, left, right), ...
                     Fi(:), ...
                     1e-10, 1000, ...
                     @(y) apply_precond( B2, V, E, y, sz ), [],...
                     X.U{idx}(:) );

            X.U{idx} = reshape( Ui, sz );

        elseif strcmpi( opts.solver, 'diag' )
            X.U{idx} = solve_inner( L.L0, X, Fi, idx );

        else
            error( 'Unknown opts.solver type. Use either ''direct'', ''pcg'' (default) or ''diag''.' )
        end

        X = orth_at( X, idx, 'right', true );

        g = apply(L, X) - F;
        residuum = [residuum; norm( g ) / normF];
        cost = [cost; cost_function_res( X, g )];
        times = [times; toc(t_start)];
    end

end

end

function res = cost_function( L, X, F )
res = 0.5*innerprod( X, apply(L, X) ) - innerprod( X, F );
end

function res = cost_function_res( X, res )
res = 0.5*innerprod( X, res );
end

function [left, right] = Afun_prepare( A, x, idx )
    y = A.apply(x);
    if idx == 1
        right = innerprod( x, y, 'RL', idx+1 );
        left = [];
    elseif idx == x.order
        left = innerprod( x, y, 'LR', idx-1 );
        right = [];
    else
        left = innerprod( x, y, 'LR', idx-1 );
        right = innerprod( x, y, 'RL', idx+1 );
    end
end

function res = Afun( A, U, idx, sz, left, right )

    V = reshape( U, sz );
    V = A.apply( V, idx );

    if idx == 1
        tmp = tensorprod_ttemps( V, right, 3 );
    elseif idx == A.order
        tmp = tensorprod_ttemps( V, left, 1 );
    else
        tmp = tensorprod_ttemps( V, right, 3);
        tmp = tensorprod_ttemps( tmp, left, 1);
    end

    res = tmp(:);
end

%function res = apply_local_precond( A, U, sz, expB)
%
%    p = size(U, 2);
%
%    x = reshape( U, [sz, p] );
%    res = zeros( [sz, p] );
%
%    for i = 1:size( expB, 1)
%        tmp = reshape( x, [sz(1), sz(2)*sz(3)*p] );
%        tmp = reshape( expB{1,i}*tmp, [sz(1), sz(2), sz(3), p] );
%
%        tmp = reshape( permute( tmp, [2 1 3 4] ), [sz(2), sz(1)*sz(3)*p] );
%        tmp = ipermute( reshape( expB{2,i}*tmp, [sz(2), sz(1), sz(3), p] ), [2 1 3 4] );
%
%        tmp = reshape( permute( tmp, [3 1 2 4] ), [sz(3), sz(1)*sz(2)*p] );
%        tmp = ipermute( reshape( expB{3,i}*tmp, [sz(3), sz(1), sz(2), p] ), [3 1 2 4] );
%
%        res = res + tmp;
%    end
%    res = reshape( res, [prod(sz), p] );
%
%end

function res = solve_inner( L0, X, Fi, idx )
    n = size(L0, 1);
    rl = X.rank(idx);
    rr = X.rank(idx+1);

    B1 = zeros( rl );
    % calculate B1 part:
    for i = 1:idx-1
        % apply L to the i'th core
        tmp = X;
        tmp.U{i} = tensorprod_ttemps( tmp.U{i}, L0, 2 );
        B1 = B1 + innerprod( X, tmp, 'LR', idx-1);
    end

    % calculate B2 part:
    B2 = L0;

    B3 = zeros( rr );
    % calculate B3 part:
    for i = idx+1:X.order
        tmp = X;
        tmp.U{i} = tensorprod_ttemps( tmp.U{i}, L0, 2 );
        B3 = B3 + innerprod( X, tmp, 'RL', idx+1);
    end

    [V,E] = eig( kron( eye(rr), B1 ) + kron( B3, eye(rl) ) );
    E = diag(E);

    rhs = matricize( Fi, 2 ) * V;
    Y = zeros(size(rhs));
    for i=1:length(E)
        Y(:,i) = (B2 + E(i)*speye(n)) \ rhs(:,i);
    end
    res = tensorize( Y*V', 2, [rl, n, rr] );
end

function [B2, V, E] = prepare_precond( L0, X, idx )
    n = size(L0, 1);
    rl = X.rank(idx);
    rr = X.rank(idx+1);

    B1 = zeros( rl );
    % calculate B1 part:
    for i = 1:idx-1
        % apply L to the i'th core
        tmp = X;
        tmp.U{i} = tensorprod_ttemps( tmp.U{i}, L0, 2 );
        B1 = B1 + innerprod( X, tmp, 'LR', idx-1);
    end

    % calculate B2 part:
    B2 = L0;

    B3 = zeros( rr );
    % calculate B3 part:
    for i = idx+1:X.order
        tmp = X;
        tmp.U{i} = tensorprod_ttemps( tmp.U{i}, L0, 2 );
        B3 = B3 + innerprod( X, tmp, 'RL', idx+1);
    end

    [V,E] = eig( kron( eye(rr), B1 ) + kron( B3, eye(rl) ) );
    E = diag(E);
end

function res = apply_precond( B2, V, E, rhs, sz )
    n = size(B2, 1);
    rhs = reshape( rhs, sz );
    rhs = matricize( rhs, 2 ) * V;
    Y = zeros(size(rhs));
    for i=1:length(E)
        Y(:,i) = (B2 + E(i)*speye(n)) \ rhs(:,i);
    end
    res = tensorize( Y*V', 2, sz );
    res = res(:);
end

author: NicolasBoumal

repo: manopt


text:
function [mean_err_rate, err_rates] = classif_mean_err_rate(labels,test_set,truth)

% Usage
%    [mean_err_rate,err_rates] =  classif_mean_err(labels,test_set,truth)
%
% Input
%    labels (int): The labels attributed to the testing instances.
%    test_set (int): The object indices of the testing instances.
%    truth: the actual labels of the testing instances
%
% Output
%    mean_err_rate (real): The mean error rate
%    recog_rate(real): array containing the individual error rates of
% each class.
%
% Description
%    This function computes the mean of the error rates of all the classes
% which are present in the dataset. We have mean_err_rate = 1 - mean_recog_rate
% where mean_recog_rate is the mean of all the recognition rate of all the
% classes which are present in the dataset and
% is a widely used measure of the performance of the classifier when the
% data set is highly unbalanced.
%
% See also
%    SVM_TEST, CLASSIF_ERR, CLASSIF_RECOG

if  isstruct(truth)
    src = truth;
    truth = [src.objects.class];
end

% Normally, the test_set contains samples of all classes
recog_rates = zeros(1,max(truth));
gdTruth = truth(:,test_set);

parfor k = 1:max(truth)

    mask = k == gdTruth;
    good_elts = find(labels==k & mask);

    mask1 = numel(find(mask==1));
    recog_rates(k)=numel(good_elts)/mask1;

end

err_rates = 1 - recog_rates;
mean_err_rate = 1 - mean(recog_rates);

end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
%% 1D Linear Advection Equation
%
% du/dt + a du/dx = 0
%
% Ref's:
%   Randall J. Leveque;
%   Finite Volume Methods for Hyperbolic Problems
%   Cambridge Press, 2002.
%
% Subroutine for solving using Upwind Method.
% by Manuel Diaz, manuel.ade'at'gmail.com
% Institute of Applied Mechanics, 2012.08.12

clear all; close all; clc;

%% Parameters
a = 0.5;
a_m = min(0,a);
a_p = max(0,a);
dx = 0.01;
cfl = 0.9;
dt = cfl*dx/abs(a);
t_end = 0.6;

%% Discretization of Domain
x = 1:dx:2;
t = 0:dt:t_end;

%% Initial Condition
n = length(x);
%u_0 = zeros(1,n);
u0(1:ceil(n/2)) = 1;
u0((ceil(n/2)+1):n) = 2;

%% Main Loop
% Initilize vector variables
u_next = zeros(size(u0)); % u^{n+1}
u_bar1 = zeros(size(u0)); % preallocate for MacCormack
u_bar2 = zeros(size(u0)); % preallocate for MacCormack

% Load Initial condition
u = u0;

% Main Loop
for k = t
    % plot every time step
    plot(x,u,'.')

    for j = 2:n-1
        % Single Sided Upwind
        %u_next(j) = u(j) - a*dt/dx*(u(j)-u(j-1)); % Upwind
        %u_next(j) = u(j) - a*dt/dx*(u(j+1)-u(j)); % Downwind

        % Double Sided Upwind
        u_next(j) = u(j) - a_p*dt/dx*(u(j)-u(j-1)) - a_m*dt/dx*(u(j+1)-u(j));

        % Comparing to Lax-Wendroff:
        %u_next(j) = u(j) - a*dt/(2*dx)*(u(j+1)-u(j-1)) ...
        %   +(dt^2/2)/(dx^2)*(a^2)*(u(j+1)-2*u(j)+u(j-1)); % LW

        % Comparing to MacCormack Method:
        % Predictor steps
        %u_bar1(j) = u(j) - a*dt/dx*(u(j+1)-u(j));
        %    u_bar1(1) = u_bar1(2); % Neumann BC in predictor step
        %u_bar2(j) = u_bar1(j) - a*dt/dx*(u_bar1(j)-u_bar1(j-1));
        % Corrector step
        %u_next(j) = 1/2*(u(j)+u_bar2(j));
    end
    % BC
    u_next(1) = u_next(2);      % Neumann BC
    u_next(n) = u_next(n-1);    % Neumann BC

    % UPDATE info
    u = u_next(1:n);

    % update figure
    pause(0.05); drawnow
end

%% Plot final Result
plot(x,u,'.')

repo: Aero-matlab


text:
function [fitness,erfc,sd_r,md] = compute_event_exceptionality(traces,N,robust_std)
%{
    Define a metric and order components according to the probabilty if some "exceptional events" (like a spike).
    Suvh probability is defined as the likeihood of observing the actual trace value over N samples given an estimated noise distribution.
    The function first estimates the noise distribution by considering the dispersion around the mode.
    This is done only using values lower than the mode. The estimation of the noise std is made robust by using the approximation std=iqr/1.349.
    Then, the probavility of having N consecutive eventsis estimated.
    This probability is used to order the components.

    Parameters:
    -----------

    traces: array
        Fluorescence traces

    Returns:
    --------

    fitness: array
        value estimate of the quality of components (the lesser the better)

    erfc: array
        probability at each time step of observing the N consequtive actual trace values given the distribution of noise

    noise_est: array
        the components ordered according to the fitness

%}

T=size(traces,2);

md = max(mode_robust(traces, 2),0);

if ~exist('N','var'); N=5; end
if ~exist('robust_std','var')
    robust_std = false;
end

ff1 = bsxfun(@minus,traces,md');
% only consider values under the mode to determine the noise standard deviation
ff1 = -ff1 .* (ff1 < 0);
if robust_std
    % compute 25 percentile
    ff1 = sort(ff1, 2);
    ff1(ff1 == 0) = nan;
    Ns = round(sum(ff1 > 0, 2) * .5);
    iqr_h = zeros(size(traces,1));
    idx = 1;
    for idx = 1:size(ff1,1)
        el = ff1(idx,:);
        iqr_h(idx) = ff1(idx, -Ns(idx));
    end
    % approximate standard deviation as iqr/1.349
    sd_r = 2 * iqr_h / 1.349;

else
    Ns = sum(ff1 > 0, 2);
    sd_r = sqrt(sum(ff1.^2, 2)./ Ns);
end
% compute z value
z = bsxfun(@times,bsxfun(@minus,traces,md'),1./(3 * sd_r));
% probability of observing values larger or equal to z given normal
% distribution with mean md and std sd_r
mu = 0;
sigma = 1;
pd = makedist('Normal',mu,sigma);
erf = 1 - cdf(pd,z);
% use logarithm so that multiplication becomes sum
erf = log(erf);
filt = ones(1,N);
% moving sum
erfc = conv2(1,filt,erf,'same');
erfc = erfc(:,1:T);

% select the maximum value of such probability for each trace
fitness = min(erfc, [], 2);

% ordered = np.argsort(fitness)

function outp = hsm_(data)
outp = [];
if numel(data) == 1
    outp = data(1);
elseif numel(data) == 2
    outp = mean(data(:));
elseif numel(data) == 3
    i1 = data(2) - data(1);
    i2 = data(3) - data(2);
    if i1 < i2
        outp = mean(data(1:2));
    elseif i2 > i1
        outp = mean(data(2:end));
    else
        outp = data(2);
    end
else

    % wMin = data[-1] - data[0]
    wMin = inf;
    N = idivide(int32(numel(data)),2,'floor') + mod(numel(data),2)-1;

    for i = 1:N
        w = data(i + N - 1) - data(i);
        if w < wMin
            wMin = w;
            j = i;
        end
    end
    outp = hsm_(data(j:j + N));
end

function dataMode = mode_robust(inputData, axis)
    %{
    Robust estimator of the mode of a data set using the half-sample mode.

    .. versionadded: 1.0.3
    """
    %}
    if ~isempty(axis)

        if axis == 2
            dataMode = arrayfun(@(n)  mode_robust(inputData(n,:),[]), 1:size(inputData,1));
        elseif axis == 1
            dataMode = arrayfun(@(n)  mode_robust(inputData(:,n),[]), 1:size(inputData,2));
        else
            error('axis can be 1 or two')
        end
    else
        % Create the function that we can use for the half-sample mode
        data = inputData(:);
        % The data need to be sorted for this to work
        data = sort(data);
        % Find the mode
        dataMode = hsm_(data);
    end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
state= acquire_bv(100, 'localhost');
while 1,
  [dmy,bn,mp,md,mt]= acquire_bv(state);
  for mm= 1:length(mt),
    %nowstr= datestr(now,31);
    %fprintf('marker %s at %s\n', md{mm}, nowstr);
    fprintf('marker %s at %d\n', md{mm}, bn+mp(mm));
  end
end
acquire_bv('close');

author: PatternRecognition

repo: OpenBMI


text:
function determ = pds_random_determinant ( n, key )

%*****************************************************************************80
%
%% PDS_RANDOM_DETERMINANT returns the determinant of the PDS_RANDOM matrix.
%
%  Discussion:
%
%    This routine will only work properly if the SAME value of SEED
%    is input that was input to PDS_RANDOM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, integer KEY, a positive value that selects the data.
%
%    Output, real DETERM, the determinant.
%
  seed = key;
  [ lambda, seed ] = r8vec_uniform_01 ( n, seed );

  determ = prod ( lambda(1:n) );

  return
end


text:
function [g1, g2] = lfmwhiteXwhiteKernGradient(lfmKern, whiteKern, t1, varargin)

% LFMWHITEXWHITEKERNGRADIENT Compute gradient between the LFM-WHITE and
% WHITE kernels.
% FORMAT
% DESC computes the gradient of an objective function with respect
% to cross kernel terms between LFM-WHITE and WHITE kernels for
% the multiple output kernel.
% ARG lfmKern : the kernel structure associated with the LFM-WHITE
% kernel.
% ARG whiteKern : the kernel structure associated with the WHITE
% kernel.
% ARG t1 : inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of objective function with respect to kernel
% parameters of LFM-WHITE kernel.
% RETURN g2 : gradient of objective function with respect to kernel
% parameters of WHITE kernel.
%
% FORMAT
% DESC computes the gradient of an objective function with respect
% to cross kernel terms between LFM-WHITE and WHITE kernels for
% the multiple output kernel.
% ARG lfmKern : the kernel structure associated with the LFM-WHITE
% kernel.
% ARG whiteKern : the kernel structure associated with the WHITE
% kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of objective function with respect to kernel
% parameters of LFM-WHITE kernel.
% RETURN g2 : gradient of objective function with respect to kernel
% parameters of WHITE kernel.
%
% SEEALSO : multiKernParamInit, multiKernCompute, lfmwhiteKernParamInit,
% whiteKernParamInit
%
% COPYRIGHT : David Luengo, 2009

% KERN

if nargin < 5
    t2 = t1;
else
    t2 = varargin{1};
end
covGrad = varargin{end};

if size(t1, 2) > 1 | size(t2, 2) > 1
  error('Input can only have one column');
end
if lfmKern.variance ~= whiteKern.variance
  error('Kernels cannot be cross combined if they have different variances.')
end

g1 = zeros(1,3);
g2 = 0; % The only parameter of the WHITE kernel (its variance) is already
        % accounted for in g1

mass = lfmKern.mass;
spring = lfmKern.spring;
damper = lfmKern.damper;
variance = lfmKern.variance;
sensitivity = lfmKern.sensitivity;
alpha = lfmKern.alpha;
omega = lfmKern.omega;

T1 = repmat(t1, 1, size(t2, 1));
T2 = repmat(t2.', size(t1, 1), 1);
deltaT = T1-T2;
indT = (T1 >= T2);
c = 1 / (mass * omega);
K0 = c * exp(-alpha*deltaT) .* indT;
K1 = K0 .* sin(omega*deltaT);
K2 = K0 .* cos(omega*deltaT);

gradMass = [1 0 0];
gradAlpha = [-damper/(2*mass^2) 1/(2*mass) 0];
c2 = sqrt(4*mass*spring-damper^2);
gradOmega = [(damper^2-2*mass*spring)/(2*c2*mass^2) -damper/(2*c2*mass) 1/c2];

% Gradient w.r.t. m_q
g1(1) = variance * sensitivity * sum(sum((gradOmega(1) * (T1-T2) .* K2 ...
    - gradAlpha(1) * (T1-T2) .* K1 ...
    - (gradMass(1)/mass + gradOmega(1)/omega) * K1) .* covGrad));

% Gradient w.r.t. D_q
g1(2) = variance * sensitivity * sum(sum((gradOmega(3) * (T1-T2) .* K2 ...
    - gradAlpha(3) * (T1-T2) .* K1 ...
    - (gradMass(3)/mass + gradOmega(3)/omega) * K1) .* covGrad));

% Gradient w.r.t. C_q
g1(3) = variance * sensitivity * sum(sum((gradOmega(2) * (T1-T2) .* K2 ...
    - gradAlpha(2) * (T1-T2) .* K1 ...
    - (gradMass(2)/mass + gradOmega(2)/omega) * K1) .* covGrad));

% Gradient w.r.t. sigma_r^2
g1(4) = sensitivity * sum(sum(K1 .* covGrad));

% Gradient w.r.t. S_{qr}
g1(5) = variance * sum(sum(K1 .* covGrad));

% Ensuring that g1 is real
g1 = real(g1);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [Zin,ip_it] = plain_foopsi(H,D,I_est,eps)

% solves argmin ||X-H||^2 subject to D*X>=0 with an interior point method
% using I_est as the initial value and eps as the initial barrier weight

ln = length(H);
step_back_frac = 0.5;
iter = 0;
if nargin == 2
    I_est = 1e-3*ones(ln,1);
    eps = 1;
end
Zin = I_est(:);

if nargin == 3
    eps = 1;
end
while eps>1e-5
    n = D*Zin;
    nnd = 10;
    E = norm(Zin-H)^2 - eps*sum(log(D*Zin));
    grad = 2*(Zin-H) - eps*D'*(n.^(-1));
    Hs = 2*speye(ln) + eps*D'*spdiags(n.^(-2),0,ln,ln)*D;
    while nnd/2>1
        iter = iter + 1;
        Z_dir = -Hs\grad;
        hit = -n./(D*Z_dir);
        if all(hit<0)
            s = 1;
        else
            s = min(1,.9*min(hit(hit>=0)));
        end
        E_new = E; s = s/step_back_frac;
        x_dg = grad'*Z_dir;
        while E_new > E + 0.25*s*x_dg
            s=s*step_back_frac;
            Z_new = Zin + s*Z_dir;
            n = D*Zin;
            E_new = norm(Z_new-H)^2 - eps*sum(log(D*Z_new));
        end
        %E = E_new;
        Zin = Zin + s*Z_dir;
        nnd = -x_dg;
        E = norm(Zin-H)^2 - eps*sum(log(D*Zin));
        n = D*Zin;
        grad = 2*(Zin-H) - eps*D'*(n.^(-1));
        Hs = 2*speye(ln) + eps*D'*spdiags(n.^(-2),0,ln,ln)*D;
        %disp(nnd)
    end
    eps = eps/10;
end
%fprintf('Interior point method converged after %i iterations \n',iter);
ip_it = iter;


text:
% [Sol] = homography2Motion(H)

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.

% Homography estimation between views of a plane
% decompose it into rotation and translation and plane normal
% p - retinal (calibrated) coordinates of the first view
% q - retinal (calibrated) coordinate in the second view
% Sol = [rotation, translation*1/d, plane normal] four possible solutions
% H - two view homography
%
% Algorith, 5.2 in Chapter 5, "An introduction to 3-D Vision"
% by Y. Ma, S. Soatto, J. Kosecka, S. Sastry (MASKS)
%
% Code distributed free for non-commercial use
% Copyright (c) MASKS, 2003
%
% Last modified 5/5/2003

function [Sol]  = homography2Motion(Hest);

 % decomposition of H into motion and structure in case of calibration case
 [u,s,v] = svd(Hest);
 H = Hest/s(2,2);
 [u,s,v] = svd(H'*H);

 if det(u) < 0 u = -u; end;

 s1 = s(1,1); s2 = s(2,2); s3 = s(3,3);
 v1 = u(:,1); v2 = u(:,2); v3 = u(:,3);
 u1 = (v1*sqrt(1-s3) + v3*sqrt(s1 -1))/sqrt(s1 - s3);
 u2 = (v1*sqrt(1-s3) - v3*sqrt(s1 -1))/sqrt(s1 - s3);

 U1 = [v2, u1, skew(v2)*u1];
 U2 = [v2, u2, skew(v2)*u2];
 W1 = [H*v2, H*u1, skew(H*v2)*H*u1];
 W2 = [H*v2, H*u2, skew(H*v2)*H*u2];

 N1 = skew(v2)*u1;
 N2 = skew(v2)*u2;

 R1 = W1*U1'
 R2 = W2*U2'

 Sol(:,:,1) = [W1*U1', (H - W1*U1')*N1, N1];
 Sol(:,:,2) = [W2*U2', (H - W2*U2')*N2, N2];
 Sol(:,:,3) = [W1*U1', -(H - W1*U1')*N1, -N1];
 Sol(:,:,4) = [W2*U2', -(H - W2*U2')*N2, -N2];

author: petercorke

repo: machinevision-toolbox-matlab


text:
function pack(descname)
%PACK Pack a descriptor to a zip file together with meta

% Copyright (C) 2016-2017 Karel Lenc
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

if nargin < 1, error('Descriptor not specified.'); end

descInfoFile = fullfile(hb_path, 'data', 'descriptors', descname, 'info.txt');
if ~exist(descInfoFile, 'file') ||...
    numel(utls.readfile(descInfoFile)) ~= 2
  fprintf('\nTo continue, we need to know few details about your submission:\n');
  submissisonName = input('Please enter submission name: ', 's');
  emailAddress = input('Please enter contact email: ', 's');
  ifd = fopen(descInfoFile, 'w');
  fprintf(ifd, '%s\n%s', submissisonName, emailAddress);
  fclose(ifd);
  fprintf('Submission info wrote to %s. Edit this file for changes.\n', ...
    descInfoFile);
end

zipFile = fullfile(hb_path(), [descname, '_descriptors.zip']);
fprintf('Packing the descriptors to %s.\n', zipFile);
zip(zipFile, fullfile(hb_path, 'data', 'descriptors', descname));

end

author: hpatches

repo: hpatches-benchmark


text:
report_this_filefun(mfilename('fullpath'));

str = [];
[newmatfile, newpath] = uiputfile([ hodo '*.dat'], 'Save As');

s = [a.Magnitude ];
le = a.Count;

fid = fopen([newpath newmatfile],'w') ;
fprintf(fid,'%5.0f %2.0f %2.0f\n',[floor(a(1,3)) a(1,4) a(1,5)]);
fprintf(fid,'%5.0f %2.0f %2.0f\n',[floor(a(le,3)) a(le,4) a(le,5)]);
fprintf(fid,'%4.1f \n',2.0);
fprintf(fid,'%4.1f \n',0.3);
fprintf(fid,'%4.1f\n',s');

fclose(fid);
clear s
cd kijko
!hn

return

%historic part

str = [];
[newmatfile, newpath] = uiputfile([ hodo '*.dat'], 'Save As');

s = [a.Date.Year a.Date.Month a.Date.Day a.Magnitude ];
le = a.Count;

fid = fopen([newpath newmatfile],'w') ;
fprintf(fid,'%5.0f %2.0f %2.0f\n',[floor(a(1,3))-1 a(1,4) a(1,5)]);
fprintf(fid,'%5.0f %2.0f %2.0f\n',[floor(a(le,3))+1 a(le,4) a(le,5)]);
fprintf(fid,'%5.2f \n',0.50);
fprintf(fid,'%5.0f %2.0f %2.0f %5.1f\n',s');

fclose(fid);

author: CelsoReyes


text:
% read_eeglabheader() - import EEGLAB dataset files
%
% Usage:
%   >> header = read_eeglabheader(filename);
%
% Inputs:
%   filename - [string] file name
%
% Outputs:
%   header   - FILEIO toolbox type structure
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, 2008-

%123456789012345678901234567890123456789012345678901234567890123456789012

% Copyright (C) 2008 Arnaud Delorme, SCCN, INC, UCSD, 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 header = read_eeglabheader(filename)

if nargin < 1
  help read_eeglabheader
  return
end

if ~isstruct(filename)
  load('-mat', filename, 'EEG');
else
  EEG = filename;
end

header.Fs          =  EEG.srate;
header.nChans      =  EEG.nbchan;
header.nSamples    =  EEG.pnts;
header.nSamplesPre = -EEG.xmin*EEG.srate;
header.nTrials     =  EEG.trials;
try
  header.label     = { EEG.chanlocs.labels }';
catch
  ft_warning('creating default channel names');
  for i=1:header.nChans
    header.label{i} = sprintf('chan%03d', i);
  end
end
ind = 1;
for i = 1:length( EEG.chanlocs )
  if isfield(EEG.chanlocs(i), 'X') && ~isempty(EEG.chanlocs(i).X)
    header.elec.label{ind, 1} = EEG.chanlocs(i).labels;
    % this channel has a position
    header.elec.elecpos(ind,1) = EEG.chanlocs(i).X;
    header.elec.elecpos(ind,2) = EEG.chanlocs(i).Y;
    header.elec.elecpos(ind,3) = EEG.chanlocs(i).Z;
    ind = ind+1;
  end
end

% remove data
% -----------
%if isfield(EEG, 'datfile')
%    if ~isempty(EEG.datfile)
%        EEG.data = EEG.datfile;
%    end;
%else
%    EEG.data = 'in set file';
%end;
EEG.icaact = [];

header.orig = EEG;

author: spm


text:
function [accel,C,S]=EGM08EarthAccel(rVec,accelOptions,M,TT1,TT2,effectsToInclude,EOP,C,S)
%%EGM08EARTHACCEL Get acceleration due to gravity from the Earth. This can
%            be a simple J2 model in a generic (not rigorously defined)
%            Earth-centered Earth-fixed (ECEF) or Earth-centered inertial
%            (ECI) coordinate system, or it can be more rigorously
%            defined in the International Terrestrial Reference System
%            (ITRS) or the Geocentric Celestial Reference System (GCRS),
%            accounting for effects such as polar motion, tides, and the
%            drift of the coefficients over time. The EGM2008 acceleration
%            model is used for parameters.
%
%INPUTS: rVec A 3XN (or 6XN) set of vectors of position (and velocity) in
%             meters in the same coordinate system as desired for the
%             output as specified by the parameter accelOptions. This is
%             the set of N positions where the acceleration due to the
%             Earth's gravity is desired. The velocity component is only
%             used if accelOptions specified an Earth-fixed coordinate
%             system. Specifically, if accelOptions is 0, 3, or 5. If a 3XN
%             vector is provided, but the velocity component is required,
%             the extra three elements are taken to be all zeros.
% accelOptions An optional parameter indicating the coordinate system and
%             options for the acceleration. Possible values are
%             0 (The default if omitted) The acceleration is found in a
%               generic ECEF coordinate system using either a Keplerian
%               model (if M=1) the J2 gravitational model (if M=2, the
%               default) with J2 given by the corresponding coefficient
%               from the zero-tide EGM2008 model without correcting for
%               polar motion, other tides, or the drift of the coefficients
%               over time. It is assumed that the rotation axis is aligned
%               with the z axis. Being an accelerating coordinate system,
%               Newtonian corrections for the Coriolis effect are used.
%               Values M>=3 are invalid. No Earth orientation parameters
%               are used. The time is also not used and neither is the
%               effectsToInclude nor C and S.
%             1 The same as 0 but in a generic Earth-centered inertial
%               coordinate system, meaning that the Coriolis effect is not
%               taken into account.
%             2 The acceleration is found using the first M elements of the
%               EGM2008 model in GCRS coordinates. The additional effects
%               accounted for (beyond just using a zero-tide model) are
%               determined by the effectsToInclude boolean vector. If C and
%               S are not provided, then the getEGMGravCoeffs function is
%               used. Earth orientation parameters are used, including the
%               length-of-day (LOD) parameter for the rotation rate of the
%               Earth.
%             3 The same as 2 but in ITRS coordinates, meaning that a
%               Newtonian Coriolis term is added to the coefficients,
%               because it is an accelerating coordinate system.
%           M The number of terms in the EGM2008 model to include. If
%             accelOptions=0 or accelOptions=1, then M chooses between
%             Keplerian, or J2 models. Otherwise, the number indicates
%             the highest order of the full EGM2008 model used. If omitted,
%             or an empty matrix is passed, a default value of 2 is used
%             for accelOptions=0 or 1 and a default of 3 is used for
%             accelOptions=2 or 3. If Inf or another number larger than the
%             highest coefficient order in the EGM2008 model is passed,
%             then the total number of coefficients in the EGM2008 model is
%             used.
%     TT1,TT2 Two parts of a Julian date given in terrestrial time (TT).
%             The units of the date are days. The full date is the sum of
%             both terms. The date is broken into two parts to provide more
%             bits of precision. It does not matter how the date is split.
%             The date is only used for accelOptions=2-3. Otherwise it can
%             be omitted or empty matrices can be passed.
% effectsToInclude A boolean vector indicating which effects are to be
%             included (if accelOptions>1). If an element of the vector is
%             1, that means that the effect will be taken into account. If
%             omitted, the default is [1;1;0;0;0]; The elements of the
%             vector are:
%             1) The drift of the low-order coefficients over time, as
%                given by the getGravCoeffOffset4Drift function.
%             2) The effects of polar motion on the coefficients, as given
%                by the getAdjustedGravCoeffs4PolarMotion function. They
%                are added after tidal effects (if any).
%             3) The effects of solid Earth tides as given by the
%                gravSolidTideOffset function.
%             4) The effects of pole tides as given by the
%                gravPoleTideOffset function.
%             5) The effects of ocean tides as given by the
%                gravOceanTideOffset function.
%         EOP A structure containing the Earth orientation parameters at
%             the given time. If omitted or an empty matrix is passed, the
%             values from the getEOP function are used. These are not used
%             if accelOptions=0 or 1. The elements of the structure are:
%             xpyp These are the polar motion coordinates in radians
%                  including the effects of tides and librations.
%             dXdY The celestial pole offsets with respect to the IAU
%                  2006/2000A precession/nutation model in radians.
%             deltaTTUT1 The difference between TT and UT1 in seconds
%              LOD The difference between the length of the day using TT
%                  and the length of the day in UT1. This is an
%                  instantaneous parameter (a derivative). The units are
%                  seconds.
%              C,S Optionally, the TIDE FREE EGM2008 coefficients with
%                  adjType=0 (the default) as obtained from the
%                  getEGMGravCoeffs function. Passing these saves a call to
%                  the getEGMGravCoeffs function and speeds up multiple
%                  calls to this function. These parameters are ignored if
%                  accelOptions=0 or 1.
%
%OUTPUTS: accel A 3XN matrix of the N accelerations due to gravity of the
%               Earth in meters per second squared in the specified
%               coordinate system (one for each input rVec).
%           C,S Values of the tide-free EGM2008 spherical harmonic
%               coefficients that can be passed to subsequent calls of this
%               function with the same options). Note that C and S are
%               modified (and then changed back), so care must be taken
%               with multithreading. Emptyt matrices are returned if
%               accelOptions=0 or accelOptions=1;
%
%For simple models, many parameters can be omitted.
%The J2 model in an ECEF system can be obtained using
%accel=EGM08EarthAccel(rVec)
%The function can be used for the Keplerian model in an ECEF system using
%accel=EGM08EarthAccel(rVec,0,1)
%or in an ECI system using
%accel=EGM08EarthAccel(rVec,1,1)
%The J2 model can be obtained in an ECI system using
%accel=EGM08EarthAccel(rVec,1,2)
%
%The J2 gravitational model is derived in Appendix E of [1].
%The Coriolis and centrifugal force models due to the use of non-inertial
%coordinate systems are derived from Appendix A of the same document.
%The various other effects are discussed in other sections.
%
%REFERENCES:
%[1] D. F. Crouse, "An Overview of Major Terrestrial, Celestial, and
%    Temporal Coordinate Systems for Target Tracking," Formal Report, Naval
%    Research Laboratory, no. NRL/FR/5344--16-10,279, 10 Aug. 2016, 173
%    pages.
%
%March 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%incorporating elements from a Coriolis correction by David Karnick.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(accelOptions))
    accelOptions=0;
end

if(nargin<3||isempty(M))
    if(accelOptions==0||accelOptions==1)
        M=2;
    else
        M=3;
    end
end
%If someone passes a number greater than the total number of coefficients,
%then just limit it to the total number in the EGM2008 model.
M=min(2190,M);

%Universal gravitational constant times the mass of the Earth
GM=Constants.EGM2008GM;

%Semi-major axis of the Earth
a=Constants.EGM2008SemiMajorAxis;

%The rotation rate of the Earth in radians per second.
omega=Constants.EGM2008EarthRotationRate;

%Check for the easy cases Keplerian or J2.
if(accelOptions==0||accelOptions==1)
    r=rVec(1:3,:);%Positions
    if(size(rVec,1)>3)
        v=rVec(4:6,:);%Velocities.
    else
        v=zeros(3,size(rVec,2));
    end

    rMag=sqrt(sum(r.*r,1));

    %Newton's law for a point or sphere.
    aNewton=-bsxfun(@times,(GM./rMag.^2),bsxfun(@rdivide,r,rMag));

    switch(M)
        case 1%Simple Keplerian dynamics.
            accel=aNewton;
        case 2%The J2 dynamic model
            %Putting this in explicitly saves a call to getEGMGravCoeffs
            %for just one parameter. This is the zero-tide value.
            C20Bar=-0.484169317366974e-03;
            C20=C20Bar*sqrt(5);
            J2=-C20;

            temp=5*r(3,:).^2./rMag.^2;
            %First-order oblateness correction.
            aOblateness=-bsxfun(@times,((3*GM*a.^2*J2./(2*rMag.^5))),[r(1,:).*(1-temp);
                                                                     r(2,:).*(1-temp);
                                                                     r(3,:).*(3-temp)]);

            accel=aNewton+aOblateness;
        otherwise
            error('Invalid value of M for the chosen accelOptions value')
    end

    %If Coriolis terms should be added due to the rotation of the
    %Earth.
    if(accelOptions==0)
        %The rotation vector for the Earth is
        %Omega=[0;0;1]*omega;
        %and the Coriolis and centrifugal forces to add to accel are
        %aCoriolis=-2*bsxfun(@cross,Omega,v);
        %aCentrifugal=-bsxfun(@cross,Omega,bsxfun(@cross,Omega,r));

        %However, it is faster to add the components directly:
        accel(1,:)=accel(1,:)+2*omega*v(2,:)+omega^2*r(1,:);
        accel(2,:)=accel(2,:)-2*omega*v(1,:)+omega^2*r(2,:);
    end
    C=[];
    S=[];
    return;
end

%Get the spherical harmonic coefficients up to the requested order if they
%are not provided.
if(nargin<8||isempty(C))
    isTideFree=true;
    modelType=0;
    [C,S]=getEGMGravCoeffs(M,isTideFree,modelType);
end

if(nargin<6||isempty(effectsToInclude))
    effectsToInclude=[1;1;0;0;0];
end

if(nargin<7||isempty(EOP))
    [JulUTC1,JulUTC2]=TT2UTC(TT1,TT2);
    [xpyp,dXdY,~,deltaTTUT1,LOD]=getEOP(JulUTC1,JulUTC2);
else
    xpyp=EOP.dxdy;
    dXdY=EOP.dXdY;
    deltaTTUT1=EOP.deltaTTUT1;
    LOD=EOP.LOD;
end

%The second one is added after all of the others.
numDelta=sum(effectsToInclude)-effectsToInclude(2);
deltaC=cell(numDelta,1);
deltaS=cell(numDelta,1);

curDelta=1;
if(effectsToInclude(1)~=false)
    %Get the offsets to C and S due to the drift of the coefficients.
    [deltaC{curDelta},deltaS{curDelta}]=getGravCoeffOffset4Drift(TT1,TT2,0);
    curDelta=curDelta+1;
end

if(effectsToInclude(3)~=false)
    %Compute the offsets due to solid Earth tides.

    [TDB1,TDB2]=TT2TDB(TT1,TT2);%Get approximate TDB.
    %Sun Position with respect to Earth.
    SunGCRSPosVel=readJPLEphem(TDB1,TDB2,11,3);
    rSunITRS=GCRS2ITRS(SunGCRSPosVel(1:3),TT1,TT2,deltaTTUT1,xpyp,dXdY);

    %Moon position with respect to Earth.
    MoonGCRSPosVel=readJPLEphem(TDB1,TDB2,10,3);
    rMoonITRS=GCRS2ITRS(MoonGCRSPosVel(1:3),TT1,TT2,deltaTTUT1,xpyp,dXdY);

    [deltaC{curDelta},deltaS{curDelta}]=gravSolidTideOffset(rMoonITRS,rSunITRS,TT1,TT2);
    curDelta=curDelta+1;
end

if(effectsToInclude(4)~=false)
    [deltaC{curDelta},deltaS{curDelta}]=gravPoleTideOffset(TT1,TT2,xpyp);
    curDelta=curDelta+1;
end

if(effectsToInclude(5)~=false)
    [deltaC{curDelta},deltaS{curDelta}]=gravOceanTideOffset(TT1,TT2);
    curDelta=curDelta+1;
end

%Find out the total number of coefficients in C and S that need to be
%changed --this is the maximum number of elements in the deltas.
numCoeffChanged=0;
for curDelta=1:numDelta
    curLength=length(deltaC{curDelta});
    if(curLength>numCoeffChanged)
        numCoeffChanged=curLength;
    end
end

if(length(C)<numCoeffChanged)
    numCoeffChanged=length(C);
end

%Now, save the elements in C and S, so that they can be restored for the
%return value when the function exits.
CSaved=C(1:numCoeffChanged);
SSaved=S(1:numCoeffChanged);

%Now, add in all of the effects, except the effects of polar motion on the
%coefficients.
for curDelta=1:numDelta
    numOffset=length(deltaC{curDelta});
    num2Change=min(numOffset,numCoeffChanged);

    C(1:num2Change)=C(1:num2Change)+deltaC{curDelta}(1:num2Change);
    S(1:num2Change)=S(1:num2Change)+deltaS{curDelta}(1:num2Change);
end

%Add in the effects of polar motion. This just changes the C21 and S21
%elements, which would otherwise be zero.
if(effectsToInclude(2)~=false)
    [~,~,C,S]=getAdjustedGravCoeffs4PolarMotion(C,S,TT1,TT2,true);
end

%Now, get the acceleration due to gravity in ITRS coordinates WITHOUT the
%Coriolis effect. This can be rotated to GCRS coordinates, if needed or
%have Coriolis terms added, if staying in this coordinate system. This
%requires getting the position in ITRS, spherical coordinates.

if(accelOptions~=3)%If r is not already in ITRS coordinates
   %Convert r (and v if present) into ITRS coordinates
   rVec=GCRS2ITRS(rVec,TT1,TT2,deltaTTUT1,xpyp,dXdY,LOD);
end

%Convert the position components into spherical coordinates for the
%spherHarmonicEval function.
r=rVec(1:3,:);%Positions
rSpher=Cart2Sphere(r);

[~,accel]=spherHarmonicEval(C,S,rSpher,a,GM);

%Now, if the acceleration is supposed to be in ITRS coordinates, then add
%the Coriolis terms. Otherwise, rotate it into GCRS coordinates.
if(accelOptions==3)
    %The output is in (approximate) ITRS coordinates, compute the Coriolis
    %terms.
    v=rVec(4:6,:);%The velocity in ITRS coordinates.

    %Get the axis of rotation in ITRS coordinates. This is the z-axis in
    %TIRS coordinates.
    omegaAxis=TIRS2ITRS([0;0;1],TT2,TT2,xpyp);

    %Modify the value of omega to deal with the LOD Earth orientation
    %parameter. This Equation is in Section IIIB of Crouse. The 86400 is
    %the number of seconds in a Julian day.
    omega=omega*(1-LOD/86400);

    %The rotation vector for the Earth
    Omega=omegaAxis*omega;

    aCoriolis=-2*bsxfun(@cross,Omega,v);
    aCentrifugal=-bsxfun(@cross,Omega,bsxfun(@cross,Omega,r));
    accel=accel+aCoriolis+aCentrifugal;
else
    %Rotate the output into GCRS coordinates. This is just a rotation (as
    %one would apply to position components). Thus, we do not want to pass
    %it as if it were a velocity input to this function.
    accel=ITRS2GCRS(accel,TT1,TT2,deltaTTUT1,xpyp,dXdY,LOD);
end

%Finally, restore the adjusted elements of C and S in case they are to be
%returned.
C(1:numCoeffChanged)=CSaved;
S(1:numCoeffChanged)=SSaved;
%Undo any changes from the getAdjustedGravCoeffs4PolarMotion function.
C(4)=0;
S(4)=0;

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
classdef(Abstract) AbstractCurveElement < matlab.mixin.SetGet & matlab.mixin.Copyable & matlab.mixin.Heterogeneous
    %AbstractCurveElement Summary of this class goes here
    %   Detailed explanation goes here

    properties
        indepVar(1,1) double %independent variable, such as indepVar
        depVar(1,1) double   %dependendent variable, such as throttleModifier
    end

    properties(Abstract, Constant)
        minIndepValue
        maxIndepValue

        minDepValue
        maxDepValue

        mandatoryIndepValues
    end

    methods(Abstract)
        indepVarName = getIndepVarName(obj)
        indepVarUnit = getIndepVarUnit(obj)
        depVarName = getDepVarName(obj)
        depVarUnit = getDepVarUnit(obj)
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot