text:
function varargout = hlp_nanocache(groupname,poolsize,func,varargin)
% Cache results of function calls.
% Results... = hlp_nanocache(Groupname,Poolsize,Function,Arguments...)
%
% The nanocache is useful when you want to cache the result of a computation for later reuse. This
% is done by comparing the given inputs with cached copies using isequal(), so you can only have a
% very small number of distinct inputs unless the inputs differ in their byte size (like data
% structures), in which case the nanocache is extremely efficient. You can also map them onto
% different groupnames (e.g., using a hash function), in which case the nanocache is also very
% efficient.
%
% In:
%   Groupname : A name that separates different groups of nanocache calls (e.g., name of the calling
%               function, or name of a key string in the Arguments). Can also be the same for every call.
%
%   Poolsize : Maximum number of entries to store per cache pool. Each group has one cache pool for
%              each Function and each length of Arguments (in bytes). The pool size limits how many
%              arguments of the same length the nanocache can hold for a given Function and Groupname.
%              Recommended number: 10.
%
%   Function : Function handle or name that should be called.
%
%   Arguments... : Arguments to pass to the function.
%
% Out:
%   Results... : Return values of the function for the given arguments, or cached copy.
%
% Limitations:
%   You cannot use the nanocache for functions whose output depends on anything other than the
%   Arguments (like global variables, files on disk, random numbers, and so on). However, you can make
%   each of these an extra argument to the function in order to use the nanocache.
%
%   If a function should give different results even through the arguments have the same value under
%   isequal, for instance when the type differs but the value is the same (like uint32(3) vs.
%   int32(3)), you need to make the type an extra argument to the function.
%
%   If the Arguments contain NaNs or anonymous functions the nanocache is not effective.
%
% Notes:
%   There is another function that supports more features to overcome most of the limitations of the
%   nanocache (but with extra overhead), called hlp_microcache.
%
%   To clear the nanocache, run the command "clear hlp_nanocache"
%
% Examples:
%   % if this line is executed for the first time, it is as slow as num2str(20.5)
%   str = hlp_nanocache('conversion',10,@num2str,20.5);
%
%   % if it is executed a second time, it is more than 60% faster than num2str(20.5)
%   str = hlp_nanocache('conversion',10,@num2str,20.5);
%
%                                   Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                   2013-10-16

% Copyright (C) Christian Kothe, SCCN, 2013, christian@sccn.ucsd.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

    persistent cache;

    % get the pool id based on the call signature (combination of nargout and input size)
    varinfo = whos('varargin');
    poolid = [char('b'+nargout) sprintf('%u',varinfo.bytes)];

    try
        % retrieve the pool of size-equivalent objects
        pool = cache.(groupname).(char(func)).(poolid);
        % search for the key in the pool
        for k=1:length(pool.keys)
            if isequal(varargin,pool.keys{k})
                varargout = pool.values{k};
                return;
            end
        end
    catch %#ok<CTCH>
        % create new pool
        pool = struct('keys',{{}},'values',{{}});
    end

    % did not find the entry: compute it
    [varargout{1:nargout}] = feval(func,varargin{:});

    % add to pool (max size per pool is limited)
    idx = 1+mod(length(pool.keys),poolsize);
    pool.keys{idx} = varargin;
    pool.values{idx} = varargout;
    cache.(groupname).(char(func)).(poolid) = pool;
end

repo: REST


text:
function [IM_result]=display_clust_image(IM,input_mask,fig)
% [IMseg]=display_clust_image(IM,input_mask,fig)
% display image segmentation results:
% input :
%          IM = input image
%          input_mask = segmentation mask
%          fig = matlab figure number for display
%
%  output:
%          IMseg = image with segments marked by different colors
%
%

if( nargin < 3 )
   fig = 0;
end

[height,width,colors] = size(IM);
segNum = max(max(input_mask));

if( size(IM,3) == 1 )
    IM = cat(3,IM,IM,IM);
end

map = [0 0 255; 150 0 255; 255 0 255; 255 0 0; 0 255 255; 0 255 87; 255 255 0; ...
        255 140 0; 128 0 0; 0 0 140; 255 196 125; 170 107 68]/255;
if( size(map,1)<segNum )
    map = hsv(segNum);
end
map_yiq = rgb2ntsc(map);

YIQ = rgb2ntsc(IM);
Y = YIQ(:,:,1);
I = YIQ(:,:,2);
Q = YIQ(:,:,3);
IM_result = ntsc2rgb(cat(3,Y,I,Q));

for i=1:segNum,
    mask = zeros(height,width);
    ind = find(input_mask==i);
    mask(ind) = 1;
    mask = medfilt2(mask);
    mask = medfilt2(mask);
    mask = medfilt2(mask);
    ind = find(mask);
    %%% mark segment boundary in white
    E = edge(mask);
    se = strel('disk',2);
    E = imdilate(E,se,'same');
    ind = find(E);
    %%% white
    Y(ind) = 1;
    I(ind) = 0;
    Q(ind) = 0;
    %%% green
    Y(ind) = 0.587;
    I(ind) = -0.2744;
    Q(ind) = -0.5299;
    %%% red
    Y(ind) = 0.2989;
    I(ind) = 0.5959;
    Q(ind) = 0.2115;
end
IM_result = ntsc2rgb(cat(3,Y,I,Q));
if( fig>0 )
    figure(fig);
    clf;
    imshow(IM_result);
end

author: jwyang


text:
function rename_brainvision_files(varargin)

    % RENAME_BRAINVISION_FILES renames a BrainVision EEG dataset, which consists of a vhdr header
    % file, vmrk marker file and a data file that usually has the extension dat, eeg or seg.
    %
    % Use as
    % rename_brainvision_files(oldname, newname, 'rmf', 'on')
    % where both the old and the new filename should be strings corresponding to the
    % header file, for instance: including the vhdr extension.
    % 'rmf' option indicates to remove old files and can be turned 'on' of 'off' (default)
    %
    % See also http://www.fieldtriptoolbox.org/ and https://sccn.ucsd.edu/wiki/EEGLAB for
    % open source software to process BrainVision EEG data.
    %
    % Robert Oostenveld https://gist.github.com/robertoostenveld/e31637a777c514bf1e86272e1092316e
    % Cyril Pernet - fixed few bugs here and there https://gist.github.com/CPernet/e037df46e064ca83a49fb4c595d4566a

    %% deal with inputs

    rmf = 'off'; % by default do not delete files
    if nargin >= 3
        if strcmpi(varargin{3}, 'rmf')
            if nargin == 3
                disp('no value associated to the remove file option, assuming off');
            else
                rmf = varargin{4};
            end
        else
            error(['unrecognized option argument ''' varargin{3} '''']);
        end
    end

    oldheaderfile = varargin{1};
    newheaderfile = varargin{2};
    clear varargin;

    % determine whether the file extensions should be in lower or upper case
    if ~isempty(regexp(newheaderfile, 'VHDR$', 'once'))
        switchcase = @upper;
    else
        switchcase = @lower;
    end

    % determine the filename without extension
    [~, f, ~] = fileparts(newheaderfile);

    %% do the renaming

    % deal with the header file
    assert(exist(oldheaderfile, 'file') ~= 0, 'the file %s does not exists', oldheaderfile);
    assert(exist(newheaderfile, 'file') == 0, 'the file %s already exists', newheaderfile);
    fid1 = fopen(oldheaderfile, 'r'); % read old
    fid2 = fopen(newheaderfile, 'w'); % write new

    while ~feof(fid1)
        line = fgetl(fid1);
        if ~isempty(regexp(line, '^MarkerFile', 'once'))
            [~, rem] = strtok(line, '=');
            oldmarkerfile = rem(2:end);
            [~, ~, x] = fileparts(oldmarkerfile);
            newmarkerfile = [f switchcase(x)];
            line = sprintf('MarkerFile=%s', newmarkerfile);
        elseif ~isempty(regexp(line, '^DataFile', 'once'))
            [~, rem] = strtok(line, '=');
            olddatafile = rem(2:end);
            [~, ~, x] = fileparts(olddatafile);
            newdatafile = [f switchcase(x)];
            line = sprintf('DataFile=%s', newdatafile);
        end
        fprintf(fid2, '%s\r\n', line);
    end
    fclose(fid1);
    fclose(fid2);

    % deal with the marker file
    if exist(oldmarkerfile, 'file') == 0
        [~, ~, ext] = fileparts(oldmarkerfile);
        [~, ff, ~] = fileparts(oldheaderfile); % re-reading as sometimes weird names comes up
        if exist([ff ext], 'file')
            oldmarkerfile = [ff ext];
        else
            error('the file %s does not exists', oldmarkerfile);
        end
    end
    assert(exist(newmarkerfile, 'file') == 0, 'the file %s already exists', newmarkerfile);
    fid1 = fopen(oldmarkerfile, 'r');
    fid2 = fopen(newmarkerfile, 'w');

    while ~feof(fid1)
        line = fgetl(fid1);
        if ~isempty(regexp(line, '^HeaderFile', 'once'))
            [~, rem] = strtok(line, '=');
            oldheaderfile = rem(2:end);
            [~, ~, x] = fileparts(oldheaderfile);
            newheaderfile = [f switchcase(x)];
            line = sprintf('HeaderFile=%s', newheaderfile);
        elseif ~isempty(regexp(line, '^DataFile', 'once'))
            [~, rem] = strtok(line, '=');
            olddatafile = rem(2:end);
            [~, ~, x] = fileparts(olddatafile);
            newdatafile = [f switchcase(x)];
            line = sprintf('DataFile=%s', newdatafile);
        end
        fprintf(fid2, '%s\r\n', line);
    end
    fclose(fid1);
    fclose(fid2);

    % deal with the data file
    if exist(olddatafile, 'file') == 0
        [~, ~, ext] = fileparts(olddatafile);
        [~, ff, ~] = fileparts(oldheaderfile); % re-reading as sometimes weird names comes up
        if exist([ff ext], 'file')
            olddatafile = [ff ext];
        else
            error('the file %s does not exists', oldmarkerfile);
        end
    end
    assert(exist(newdatafile, 'file') == 0, 'the file %s already exists', newdatafile);
    status = copyfile(olddatafile, newdatafile);
    if ~status
        error('failed to copy data from %s to %s', olddatafile, newdatafile);
    end

    %% delete old files *try* in case of user restriction
    if strcmpi(rmf, 'on')
        try
            delete(oldheaderfile);
        end
        try
            delete(oldmarkerfile);
        end
        try
            delete(olddatafile);
        end
    end

author: bids-standard

repo: bids-starter-kit


text:
function [sys,x0,str,ts] = nrotor_dynamics(t,x,u,flag, vehicle, x0, groundflag)
    % Flyer2dynamics lovingly coded by Paul Pounds, first coded 12/4/04
    % A simulation of idealised X-4 Flyer II flight dynamics.
    % version 2.0 2005 modified to be compatible with latest version of Matlab
    % version 3.0 2006 fixed rotation matrix problem
    % version 4.0 4/2/10, fixed rotor flapping rotation matrix bug, mirroring
    % version 5.0 8/8/11, simplified and restructured
   % version 6.0 25/10/13, fixed rotation matrix/inverse wronskian definitions, flapping cross-product bug

    % modified by PIC for N rotors
    global groundFlag

    warning off MATLAB:divideByZero

    % New in version 2:
    %   - Generalised rotor thrust model
    %   - Rotor flapping model
    %   - Frame aerodynamic drag model
    %   - Frame aerodynamic surfaces model
    %   - Internal motor model
    %   - Much coolage

    % Version 1.3
    %   - Rigid body dynamic model
    %   - Rotor gyroscopic model
    %   - External motor model

    %ARGUMENTS
    %   u       Reference inputs                1x4
    %   tele    Enable telemetry (1 or 0)       1x1
    %   crash   Enable crash detection (1 or 0) 1x1
    %   init    Initial conditions              1x12

    %INPUTS
    %   u      rotor speed                      1xN

    %CONTINUOUS STATES
    %   z      Position                         3x1   (x,y,z)
    %   v      Velocity                         3x1   (xd,yd,zd)
    %   n      Attitude                         3x1   (Y,P,R)
    %   o      Angular velocity                 3x1   (wx,wy,wz)
    %   w      Rotor angular velocity           4x1

    %CONTINUOUS STATE MATRIX MAPPING
    %   x = [z1 z2 z3 n1 n2 n3 z1 z2 z3 o1 o2 o3 w1 w2 w3 w4]

    %INITIAL CONDITIONS
    n0 = [0 0 0];               %   n0      Ang. position initial conditions    1x3
    v0 = [0 0 0];               %   v0      Velocity Initial conditions         1x3
    o0 = [0 0 0];               %   o0      Ang. velocity initial conditions    1x3
    init = [x0 n0 v0 o0];       % x0 is the passed initial position 1x3
    groundFlag = groundflag;;

    %CONTINUOUS STATE EQUATIONS
    %   z' = v
    %   v' = g*e3 - (1/m)*T*R*e3
    %   I*o' = -o X I*o + G + torq
    %   R = f(n)
    %   n' = inv(W)*o

    % basic sanity checks on number of rotors
    if ~isfield(vehicle, 'nrotors')
         error('RTB:nrotor_dynamics:bardarg', 'Number of rotors not specified in model structure')
    end
    if mod(vehicle.nrotors,2) == 1
        error('RTB:nrotor_dynamics:bardarg', 'Number of rotors must be even')
    end
    if vehicle.nrotors < 4
        error('RTB:nrotor_dynamics:bardarg', 'Number of rotors must be at least 4')
    end

    groundFlag = groundflag;

    % Dispatch the flag.
    switch flag
        case 0
            [sys,x0,str,ts]=mdlInitializeSizes(init, vehicle); % Initialization
        case 1
            sys = mdlDerivatives(t,x,u, vehicle); % Calculate derivatives
        case 3
            sys = mdlOutputs(t,x, vehicle); % Calculate outputs
        case { 2, 4, 9 } % Unused flags
            sys = [];
        otherwise
            error(['Unhandled flag = ',num2str(flag)]); % Error handling
    end
end % End of flyer2dynamics

%==============================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the
% S-function.
%==============================================================
%
function [sys,x0,str,ts] = mdlInitializeSizes(init, vehicle)
    %
    % Call simsizes for a sizes structure, fill it in and convert it
    % to a sizes array.
    %
    sizes = simsizes;
    sizes.NumContStates  = 12;
    sizes.NumDiscStates  = 0;
    sizes.NumOutputs     = 12;
    sizes.NumInputs      = vehicle.nrotors;
    sizes.DirFeedthrough = 0;
    sizes.NumSampleTimes = 1;
    sys = simsizes(sizes);
    %
    % Initialize the initial conditions.
    x0 = init;
    %
    % str is an empty matrix.
    str = [];
    %
    % Generic timesample
    ts = [0 0];

    if vehicle.verbose
        disp(sprintf('t\t\tz1\t\tz2\t\tz3\t\tn1\t\tn2\t\tn3\t\tv1\t\tv2\t\tv3\t\to1\t\to2\t\to3\t\tw1\t\tw2\t\tw3\t\tw4\t\tu1\t\tu2\t\tu3\t\tu4'))
    end
end % End of mdlInitializeSizes.

%==============================================================
% mdlDerivatives
% Calculate the state derivatives for the next timestep
%==============================================================
%
function sys = mdlDerivatives(t,x,u, quad)
    global a1s b1s

    for i = 1:quad.nrotors
        theta = (i-1)/quad.nrotors*2*pi;
        %   Di      Rotor hub displacements (1x3)
        % first rotor is on the x-axis, clockwise order looking down from above
        D(:,i) = [ quad.d*cos(theta); quad.d*sin(theta); quad.h];
    end

    %Body-fixed frame references
    e1 = [1;0;0];               %   ei      Body fixed frame references         3x1
    e2 = [0;1;0];
    e3 = [0;0;1];

    %EXTRACT ROTOR SPEED DEMANDS FROM U
    w = u(1:quad.nrotors);

    %EXTRACT STATES FROM X
    z = x(1:3);   % position in {W}
    n = x(4:6);   % RPY angles {W}
    v = x(7:9);   % velocity in {W}
    o = x(10:12); % angular velocity in {W}

    %PREPROCESS ROTATION AND WRONSKIAN MATRICIES
    phi = n(1);    % yaw
    the = n(2);    % pitch
    psi = n(3);    % roll

    % rotz(phi)*roty(the)*rotx(psi)
    R = [cos(the)*cos(phi) sin(psi)*sin(the)*cos(phi)-cos(psi)*sin(phi) cos(psi)*sin(the)*cos(phi)+sin(psi)*sin(phi);   %BBF > Inertial rotation matrix
         cos(the)*sin(phi) sin(psi)*sin(the)*sin(phi)+cos(psi)*cos(phi) cos(psi)*sin(the)*sin(phi)-sin(psi)*cos(phi);
         -sin(the)         sin(psi)*cos(the)                            cos(psi)*cos(the)];

    %Manual Construction
    %     Q3 = [cos(phi) -sin(phi) 0;sin(phi) cos(phi) 0;0 0 1];   % RZ %Rotation mappings
    %     Q2 = [cos(the) 0 sin(the);0 1 0;-sin(the) 0 cos(the)];   % RY
    %     Q1 = [1 0 0;0 cos(psi) -sin(psi);0 sin(psi) cos(psi)];   % RX
    %     R = Q3*Q2*Q1    %Rotation matrix
    %
    %    RZ * RY * RX
    iW = [0        sin(psi)          cos(psi);             %inverted Wronskian
          0        cos(psi)*cos(the) -sin(psi)*cos(the);
          cos(the) sin(psi)*sin(the) cos(psi)*sin(the)] / cos(the);

    %ROTOR MODEL

    for i=1:quad.nrotors %for each rotor
        %Relative motion

        Vr = cross(o,D(:,i)) + v;
        mu = sqrt(sum(Vr(1:2).^2)) / (abs(w(i))*quad.r);  %Magnitude of mu, planar components
        lc = Vr(3) / (abs(w(i))*quad.r);   %Non-dimensionalised normal inflow
        li = mu; %Non-dimensionalised induced velocity approximation
        alphas = atan2(lc,mu);
        j = atan2(Vr(2),Vr(1));  %Sideslip azimuth relative to e1 (zero over nose)
        J = [cos(j) -sin(j);
            sin(j) cos(j)];  %BBF > mu sideslip rotation matrix

        %Flapping
        beta = [((8/3*quad.theta0 + 2*quad.theta1)*mu - 2*(lc)*mu)/(1-mu^2/2); %Longitudinal flapping
            0;];%sign(w) * (4/3)*((Ct/sigma)*(2*mu*gamma/3/a)/(1+3*e/2/r) + li)/(1+mu^2/2)]; %Lattitudinal flapping (note sign)
        beta = J'*beta;  %Rotate the beta flapping angles to longitudinal and lateral coordinates.
        a1s(i) = beta(1) - 16/quad.gamma/abs(w(i)) * o(2);
        b1s(i) = beta(2) - 16/quad.gamma/abs(w(i)) * o(1);

        %Forces and torques
        T(:,i) = quad.Ct*quad.rho*quad.A*quad.r^2*w(i)^2 * [-cos(b1s(i))*sin(a1s(i)); sin(b1s(i));-cos(a1s(i))*cos(b1s(i))];   %Rotor thrust, linearised angle approximations
        Q(:,i) = -quad.Cq*quad.rho*quad.A*quad.r^3*w(i)*abs(w(i)) * e3;     %Rotor drag torque - note that this preserves w(i) direction sign
        tau(:,i) = cross(T(:,i),D(:,i));    %Torque due to rotor thrust
    end

    %RIGID BODY DYNAMIC MODEL
    dz = v;

    % vehicle can't fall below ground
    if groundFlag && (z(3) > 0)
        z(3) = 0;
        dz(3) = 0;
    end

    dn = iW*o;
    dv = quad.g*e3 + R*(1/quad.M)*sum(T,2);
    do = inv(quad.J)*(cross(-o,quad.J*o) + sum(tau,2) + sum(Q,2)); %row sum of torques
    sys = [dz;dn;dv;do];   %This is the state derivative vector
end % End of mdlDerivatives.

%==============================================================
% mdlOutputs
% Calculate the output vector for this timestep
%==============================================================
%
function sys = mdlOutputs(t,x, quad)
    %CRASH DETECTION
    if x(3)>0
        x(3) = 0;
        if x(6) > 0
            x(6) = 0;
        end
    end
    %if (x(3)>0)&(crash)
    %    error('CRASH!')
    %end

    %TELEMETRY
    if quad.verbose
        disp(sprintf('%0.3f\t',t,x))
    end

    % compute output vector as a function of state vector
    %   z      Position                         3x1   (x,y,z)
    %   v      Velocity                         3x1   (xd,yd,zd)
    %   n      Attitude                         3x1   (Y,P,R)
    %   o      Angular velocity                 3x1   (Yd,Pd,Rd)

    n = x(4:6);   % RPY angles
    phi = n(1);    % yaw
    the = n(2);    % pitch
    psi = n(3);    % roll

    % rotz(phi)*roty(the)*rotx(psi)
    R = [cos(the)*cos(phi) sin(psi)*sin(the)*cos(phi)-cos(psi)*sin(phi) cos(psi)*sin(the)*cos(phi)+sin(psi)*sin(phi);   %BBF > Inertial rotation matrix
         cos(the)*sin(phi) sin(psi)*sin(the)*sin(phi)+cos(psi)*cos(phi) cos(psi)*sin(the)*sin(phi)-sin(psi)*cos(phi);
         -sin(the)         sin(psi)*cos(the)                            cos(psi)*cos(the)];

    iW = [0        sin(psi)          cos(psi);             %inverted Wronskian
          0        cos(psi)*cos(the) -sin(psi)*cos(the);
          cos(the) sin(psi)*sin(the) cos(psi)*sin(the)] / cos(the);

    % return velocity in the body frame
    sys = [ x(1:6);
            inv(R)*x(7:9);   % translational velocity mapped to body frame
            iW*x(10:12)];    % RPY rates mapped to body frame
    %sys = [x(1:6); iW*x(7:9);  iW*x(10:12)];
    %sys = x;
end
% End of mdlOutputs.

author: petercorke

repo: robotics-toolbox-matlab


text:
function g = linearOutputGrad(model, X, dim)

% LINEAROUTPUTGRAD Evaluate derivatives of linear model outputs with respect to parameters.
%
% MODIFICATIONS : Carl Henrik Ek, 2009
%
% MLTOOLS

numData = size(X, 1);
if(nargin<=2)
  for i = 1:model.outputDim
    startZeros = zeros(numData, model.inputDim*(i - 1));
    finishZeros = zeros(numData, model.inputDim*(model.outputDim-i));
    startZeros2 = zeros(numData, (i - 1));
    finishZeros2 = zeros(numData, (model.outputDim-i));
    g(:, :, i) = [startZeros X finishZeros startZeros2 ones(numData, 1) finishZeros2];
  end
else
  g(:,:) = [X ones(numData,1)];
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function triangle_dunavant_rule_test04 ( )

%*****************************************************************************80
%
%% TRIANGLE_DUNAVANT_RULE_TEST04 tests DUNAVANT_RULE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 December 2006
%
%  Author:
%
%    John Burkardt
%
  rule_max = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGLE_DUNAVANT_RULE_TEST04\n' );
  fprintf ( 1, '  DUNAVANT_RULE returns the points and weights of\n' );
  fprintf ( 1, '  a Dunavant rule for the unit triangle.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This routine uses those rules to estimate the\n' );
  fprintf ( 1, '  integral of monomomials in the unit triangle.\n' );

  area = 0.5;

  for a = 0 : 10

    for b = 0 : 10 - a
%
%  Multiplying X^A * Y^B by COEF will give us an integrand
%  whose integral is exactly 1.  This makes the error calculations easy.
%
      coef = ( a + b + 2 ) * ( a + b + 1 );
      for i = 1 : b
        coef = coef * ( a + i ) / i;
      end

      fprintf ( 1, '\n' );
      fprintf ( 1, '  Integrate %f * X^%d * Y^%d\n', coef, a, b );
      fprintf ( 1, '\n' );
      fprintf ( 1, '      Rule       QUAD           ERROR\n' );
      fprintf ( 1, '\n' );

      for rule = 1 : rule_max

        order_num = dunavant_order_num ( rule );

        [ xy, w ] = dunavant_rule ( rule, order_num );

        quad = 0.0;

        for order = 1 : order_num

          x = xy(1,order);
          y = xy(2,order);

          if ( a == 0 & b == 0 )
            value = coef;
          elseif ( a == 0 & b ~= 0 )
            value = coef * y^b;
          elseif ( a ~= 0 & b == 0 )
            value = coef * x^a;
          elseif ( a ~= 0 & b ~= 0 )
            value = coef * x^a * y^b;
          end

          quad = quad + w(order) * value;

        end

        quad = area * quad;

        exact = 1.0;
        err = abs ( exact - quad );

        fprintf ( 1, '  %8d  %14f  %14f\n', rule, quad, err );

      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% computes registration offsets for data split into blocks
% loops over blocks and returns offsets dsall
function [dsall,ops1] = nonrigidOffsets(data, k, j, iplane0, ops, ops1)

nplanes = getOr(ops, {'nplanes'}, 1);
alignAcrossPlanes  = getOr(ops, {'alignAcrossPlanes'}, false);
planesToInterpolate = getOr(ops, {'planesToInterpolate'}, 1:nplanes);

nblocks = ops.numBlocks(1)*ops.numBlocks(2);
dsall = zeros(size(data,3), 2, nblocks);

for i = 1:numel(ops.planesToProcess)
    if alignAcrossPlanes && ismember(i, planesToInterpolate) % load frames of all planes
        indframes = 1:size(data,3);
    else
        ifr0 = iplane0(ops.planesToProcess(i));
        indframes = ifr0:ops.nplanes:size(data,3);
    end

    ds = zeros(numel(indframes), 2, nblocks,'double');
    Corr = zeros(numel(indframes), nblocks,'double');
    for ib = 1:nblocks
        % collect ds
        ops1{i}.mimg = ops1{i}.mimgB{ib};
	if ops.kriging
	  [ds(:,:,ib), Corr(:,ib)]  = ...
	      regoffKriging(data(ops1{i}.yBL{ib},ops1{i}.xBL{ib},indframes),ops1{i}, 0);
	else
	  [ds(:,:,ib), Corr(:,ib)]  = ...
	      regoffLinear(data(ops1{i}.yBL{ib},ops1{i}.xBL{ib},indframes),ops1{i},0);
	end
    end
    if j==1
        ds(1,:,:) = 0;
    end
    dsall(indframes,:,:)  = ds;
    ops1{i}.DS          = cat(1, ops1{i}.DS, ds);
    ops1{i}.CorrFrame   = cat(1, ops1{i}.CorrFrame, Corr);
    ops1{i}.Nframes(k)  = ops1{i}.Nframes(k) + length(indframes);
end

author: cortex-lab


text:
function jed = mjd_to_jed ( mjd )

%*****************************************************************************80
%
%% MJD_TO_JED converts a modified JED to a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real MJD, the modified Julian Ephemeris Date.
%
%    Output, real JED, the Julian Ephemeris Date.
%
  jed_epoch = epoch_to_jed_mjd ( );
  jed = mjd + jed_epoch;

  return
end


text:
function boundedIndices=boundsCheck3D(minCoords,maxCoords,inputCoords)
% function boundedIndices=boundsCheck(minCoords,maxCoords,inputCoords)
% PURPOSE:  returns the indices of those points in inputCoords
% That lie within the bounding box specified by the min and max coordinates
% Works with 3D coordinates.
% AUTHOR : Wade
% DATE : 020801
% Last: $Date: 2007/07/05 19:51:58 $
% Do some size checks
[nInputCoords,nInpDims]=size(inputCoords);
if (nInpDims~=3)
    error('In boundsCheck3D, inputCoords must be n*3');
end
if ((prod(size(minCoords))~=3) | (prod(size(maxCoords))~=3))
 error ('In boundsCheck3D, min and max coords must be 3*1 vectors');
end
% do bounds check
okPoints=(inputCoords(:,1)>minCoords(1)).*(inputCoords(:,2)>minCoords(2)).*(inputCoords(:,3)>minCoords(3));
okPoints=okPoints.*(inputCoords(:,1)<maxCoords(1)).*(inputCoords(:,2)<maxCoords(2)).*(inputCoords(:,3)<maxCoords(3));
boundedIndices=find(okPoints);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
clear all; close all; clc;

t=linspace(-pi,pi,10000);
delay=-pi/2;
delay2=0;
amp1=0.8;
amp2=1;
amp3=1;
slowness=150;

for ii=1:1:length(t)

    A=amp1*sin(amp1*2*pi*1*(ii/slowness));
    B=amp2*sin(amp2*2*pi*1*(ii/slowness)+delay);
    C=amp3*sin(amp3*2*pi*1*(ii/slowness)+delay2);
    sin_graph=amp1*sin(amp1*2*pi*1*(t+ii/slowness));
    cos_graph=amp2*sin(amp2*2*pi*1*(t+ii/slowness)+delay);
    third_graph=amp3*sin(amp3*2*pi*1*(t+ii/slowness)+delay2);

    subplot(411);
    plot(t,sin_graph);hold on;
    plot(0,A,'r*','Markersize',10);hold off;title('Signal A : X-Axis');
    xlim([-pi pi]);ylim([-1 1]);
    set(gcf,'Position',[250,100, 800, 800]);
    grid on;
    subplot(412);
    T=num2str(delay);
    Title_for_B=strcat('Signal B with Delay = ',T,' : Y-Axis');
    plot(t,cos_graph);hold on;title(Title_for_B);
    xlim([-pi pi]);ylim([-1 1]);
    plot(0,B,'r*','MarkerSize',10);hold off;
    grid on;
    subplot(413);
    T3=num2str(delay2);
    Title_for_C=strcat('Signal C with Delay = ',T3,' : Z-Axis');
    plot(t,third_graph);hold on;title(Title_for_C);
    xlim([-pi pi]);ylim([-1 1]);
    plot(0,C,'r*','MarkerSize',10);hold off;
    grid on;
    subplot(338);
%     T2=num2str(amp1/amp2);
%     Title_for_C=strcat('Ratio A/B = ',T2);
    plot3(sin_graph,cos_graph,third_graph);%title(Title_for_C);
    hold on;
    plot3(A,B,C,'r*','MarkerSize',10);hold off;%xlim([-1 1]); ylim([-1 1]);
    grid on;
%     drawnow;
    pause(0.00000001)

end

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function ispy(A,sizeMultiplier)
%ISPY Visualize sparsity pattern and the size of the elements of a matrix.
%
%   ISPY(S) plots the sparsity pattern of the matrix S and visualizes the
%   size and sign of its elements. Positive elements are plotted in
%   red and negative elements in blue. Zero elements are not plotted.
%   The size of an element is scaled logarithmically.
%
%   The imaginary part of any number is discarded.
%   NaNs are plotted as red stars on a yellow background.
%   Infs are plotted as blue stars on a green background.
%
%   If the matrix is very big, try zooming in to get a more accurate view.
%
%   If the visualisation should rely on coloring instead of sizes, or if the
%   matrix is really big, then you should consider using imagesc instead.
%
%   sizeMultiplier allows for manual tuning of the visualized sizes.
%
%   Examples:
%      ispy(gallery('wathen',15,15))
%      ispy(gallery('orthog',50,5))
%
%   See also SPY and IMAGESC.

%   -- v1.0 / 2013-08-20 / Samu Alanko, Patrick N. Raanes

% Input checks
if nargin<2
    sizeMultiplier = 1;
end
if nnz(A) == 0
    fprintf('Matrix is empty! Exiting. \n');
    return;
end

% Sizes
[M,N] = size(A);

% We neglect the complex parts of the elements
A = real(sparse(A));

% All nonzero elements
[Xnz,Ynz,Snz] = find(A);
% Finite elems
inds = isfinite(Snz);
X = Xnz(inds); Y = Ynz(inds); S = Snz(inds);
% Infinite elems
inds = isinf(Snz);
Xinf = Xnz(inds); Yinf = Ynz(inds);
% NaN elems
inds = isnan(Snz);
Xnan = Xnz(inds); Ynan = Ynz(inds);

% Get indices of positive and negative elements
posind = S>0;
negind = S<0;

% Set up figure
fh = gcf; % fig handle
ah = gca; % axes handle
xlim([0 N+1]); ylim([0 M+1]);

% Main routine
drawelems();

% Beautify
title(inputname(1));
box on;
set(ah, 'YDir', 'reverse')

% Draws the elements that are inside the current view (makes zoom fast)
function drawelems()

    % Get axes size and multiply the sizeMultiplier by the
    % ratio (axes area)/(standard axes area)
    [axW,axH] = getAxesSize(ah);
    sm = sizeMultiplier * axW*axH/1.2e5;

    % Elements included in current zoom
    xlims = get(ah,'YLim');
    ylims = get(ah,'XLim');
    width = abs(xlims(2)-xlims(1));

    % Select the elements that are inside the current view
    ninds = negind & X>=xlims(1) & X<=xlims(2) & Y>=ylims(1) & Y<=ylims(2);
    pinds = posind & X>=xlims(1) & X<=xlims(2) & Y>=ylims(1) & Y<=ylims(2);

    nn = sum(ninds); % num of neg elems
    np = sum(pinds); % num of pos elems

    Xp = [X(ninds); X(pinds)];
    Yp = [Y(ninds); Y(pinds)];
    Sp = [S(ninds); S(pinds)];

    % Markersizes - constants
    minmarkersize = 4;              % Minimum marker size
    maxmarkersize = 40*sm;          % Maximum marker size
    markersize    = 12*40*sm/width; % Average size

    % Markersize of elements
    meanval = mean(abs(Sp));
    Sp = markersize.*(log(abs(Sp)./meanval+1)/log(2));
    Sp = min(max(Sp,minmarkersize*ones(size(Sp))),maxmarkersize*ones(size(Sp)));
    Sp = markersize*Sp/mean(abs(Sp));

    % Plot
    cla; hold on;
    scatter(Yp(1:nn),    Xp(1:nn),    Sp(1:nn),    'b','filled'); % neg elems
    scatter(Yp(nn+1:end),Xp(nn+1:end),Sp(nn+1:end),'r','filled'); % pos elems

    % NaN and Inf
    if ~isempty(Xnan)
        scatter(Ynan,Xnan,1000*sm,'y.');
        scatter(Ynan,Xnan,80*sm,'r*');
    end
    if ~isempty(Xinf)
        scatter(Yinf,Xinf,1000*sm,'g.');
        scatter(Yinf,Xinf,80*sm,'b*');
    end

    hold off;

    % Attach callback functions. Also make them detachable if the figure
    % is reused. But if it's reused in resizing/zooming, then we must attach
    % the callbacks again. That's why the callbacks are inside drawelems().
    zh = zoom(fh);
    ph = pan(fh);
    set(zh,'ActionPostCallback',@zoomcallback);
    set(ph,'ActionPostCallback',@pancallback);
    set(fh,'ResizeFcn',@resizecallback);
    set(fh,'NextPlot','Replace');

end

% Callback functions
function zoomcallback(~,~)
    drawelems();
end

function pancallback(~,~)
    drawelems();
end

function resizecallback(~,~)
    drawelems();
end

end

% Get axes size in "points" units
function [w,h] = getAxesSize(ah)
oldUnits = get(ah,'Units');
set(ah, 'Units', 'points');
axPos = get(ah,'Position');
set(ah,'Units',oldUnits);
w = axPos(3);
h = axPos(4);
end

author: Sable

repo: mcbench-benchmarks


text:
%   CLASS QZSS_SS
% =========================================================================
%
% DESCRIPTION
%   container of QZSS Satellite System parameters
%
% REFERENCES
%   CRS parameters, according to each GNSS system CRS definition
%   (ICD document in brackets):
%
%   *_QZS --> GRS80    (IS-QZSS 1.8E)
%   Standard: http://qz-vision.jaxa.jp/USE/is-qzss/DOCS/IS-QZSS_18_E.pdf
%
%   Other useful links
%     - http://www.navipedia.net/index.php/QZSS_Signal_Plan
%     - Ellipsoid: http://www.unoosa.org/pdf/icg/2012/template/WGS_84.pdf
%       note that GM and OMEGAE_DOT are redefined in the standard IS-GPS200H (GPS is not using WGS_84 values)
%     - http://www.navipedia.net/index.php/Reference_Frames_in_GNSS
%     - http://gage6.upc.es/eknot/Professional_Training/PDF/Reference_Systems.pdf

%--------------------------------------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:        Andrea Gatti
%  Contributors:      Andrea Gatti, Giulio Tagliaferro ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

classdef QZSS_SS < Satellite_System
    properties (Constant, Access = 'public')
        SYS_EXT_NAME = 'QZSS';    % full name of the constellation
        SYS_NAME     = 'QZS';     % 3 characters name of the constellation, this "short name" is used as fields of the property list (struct) to identify a constellation
        SYS_C        = 'J';       % Satellite system (ss) character id

        % System frequencies as struct [MHz]
        F = struct('L1', 1575.420, ...
                   'L2', 1227.600, ...
                   'L5', 1176.450, ...
                   'LEX6', 1278.750)

        % Array of supported frequencies [MHz]
        F_VEC = struct2array(QZSS_SS.F) * 1e6;

        % Array of the corresponding wavelength - lambda => wavelengths
        L_VEC = 299792458 ./ QZSS_SS.F_VEC;

        N_SAT = 9;           % Maximum number of satellite in the constellation
        PRN = (1 : 9)';      % Satellites id numbers as defined in the constellation

        % CODE2DATA ftp://igs.org/pub/data/format/rinex303.pdf
        CODE_RIN3_ATTRIB  = {'ZXLSC F' 'XLS F', 'XIQ F', 'XLS F'}; % last letter of the observation code
        CODE_RIN3_DEFAULT_ATTRIB  = {'C' 'S' 'Q' 'S'}; % last letter of the observation code
        CODE_RIN3_2BAND  = '1256';                        % id for the freq as stored in F_VEC
        IONO_FREE_PREF  = ['12';'15';'16';'25';'26';'56']; % to be evaluated which combination is really better
    end

    properties (Constant, Access = 'private')
        % QZSS (GRS80) Ellipsoid semi-major axis [m]
        ELL_A = 6378137;
        % QZSS (GRS80) Ellipsoid flattening
        ELL_F = 1/298.257222101;
        % QZSS (GRS80) Ellipsoid Eccentricity^2
        ELL_E2 = (1 - (1 - QZSS_SS.ELL_F) ^ 2);
        % QZSS (GRS80) Ellipsoid Eccentricity
        ELL_E = sqrt(QZSS_SS.ELL_E2);
    end

    properties (Constant, Access = 'public')
        % Structure of orbital parameters (ellipsoid, GM, OMEGA_EARTH_DOT)
        ORBITAL_P = struct('GM', 3.986005e14, ...                  % Gravitational constant * (mass of Earth) [m^3/s^2]
                                    'OMEGAE_DOT', 7.2921151467e-5, ...      % Angular velocity of the Earth rotation [rad/s]
                                    'ELL',struct( ...                       % Ellipsoidal parameters QZSS (GRS80)
                                    'A', QZSS_SS.ELL_A, ...             % Ellipsoid semi-major axis [m]
                                    'F', QZSS_SS.ELL_F, ...             % Ellipsoid flattening
                                    'E', QZSS_SS.ELL_E, ...             % Eccentricity
                                    'E2', QZSS_SS.ELL_E2));             % Eccentricity^2
        ORBITAL_INC = 48;    % Orbital inclination
        ORBITAL_RADIUS  = 42164000 + 6378137; % Orbital radius
    end

    methods
        function this = QZSS_SS(offset)
            % Creator
            % SYNTAX: QZSS_SS(<offset>)
            if (nargin == 0)
                offset = 0;
            end
            this@Satellite_System(offset);
        end

        function copy = getCopy(this)
            % Get a copy of this
            copy = QZSS_SS(this.getOffset());
            copy.import(this);
        end
    end
end

author: goGPS-Project


text:
function u = minL2Potts2DADMM(f, gamma, varargin)
%minL2Potts2DADMM A ADMM splitting strategy for the two-dimensional vector valued L^2-Potts
%problem
%
%Description:
% Minimizes the (2D) Potts problem
%
%  \gamma \| u \|_0 + \| A u - f \|_p^p -> min
%
% using ADMM splitting an dynamic programming
%
%Reference:
% M. Storath, A. Weinmann
% Fast partitioning of vector-valued images",
% SIAM Journal on Imaging Sciences, 2014

% written by M. Storath
% $Date: 2014-06-30 11:26:34 +0200 (Mo, 30. Jun 2014) $	$Revision: 99 $

[m,n,~] = size(f);

% parse options
ip = inputParser;
addParamValue(ip, 'muInit', gamma*1e-2);
addParamValue(ip, 'muStep', 2);
addParamValue(ip, 'tol', 1e-10);
addParamValue(ip, 'isotropic', true);
addParamValue(ip, 'verbose', false);
addParamValue(ip, 'weights', ones(m,n));
addParamValue(ip, 'multiThreading', true);
addParamValue(ip, 'quantization', true);
addParamValue(ip, 'useADMM', true);
parse(ip, varargin{:});
par = ip.Results;

% check args
assert(par.muStep > 1, 'Variable muStep must be > 1.');
assert(all(par.weights(:) >= 0), 'Weights must be >= 0.');
assert(par.tol > 0, 'Stopping tolerance must be > 0.');
assert(par.muInit > 0, 'muInit must be > 0.');

% cast data to PLImage
plf = pottslab.PLImage(f);

% main program (calls Java routines)
if par.isotropic
    % near-isotropic discretization
    omega(1) = sqrt(2.0) - 1.0;
    omega(2) = 1.0 - sqrt(2.0)/2.0;
    % alternative neighborhood weights
    %omega(1) = (2 * sqrt(2.0) - 1.0)/3;
    %omega(2) = (2 - sqrt(2.0))/3;
    plu = pottslab.JavaTools.minL2PottsADMM8(plf, gamma, par.weights, par.muInit, par.muStep, par.tol, par.verbose, par.multiThreading, par.useADMM, omega);

else
    % anisotropic discretization
    plu = pottslab.JavaTools.minL2PottsADMM4(plf, gamma, par.weights, par.muInit, par.muStep, par.tol, par.verbose, par.multiThreading, par.useADMM);
end

% reshape the 1D array given by .toDouble()
u = reshape( plu.toDouble(), size(f) );

% to remove small remaining variations in result (algorithm works with floats)
if par.quantization
    u = round(u * 255)/255;
end

end

author: mstorath

repo: Pottslab

save_path: github-repos/MATLAB/mstorath-Pottslab


text:
function [varargout] = ft_notification(varargin)

% FT_NOTIFICATION works mostly like the WARNING and ERROR commands in MATLAB and
% is called by FT_ERROR, FT_WARNING, FT_NOTICE, FT_INFO and FT_DEBUG. Please note
% that you should not call this function directly.
%
% Some examples:
%  ft_info on
%  ft_info on msgId
%  ft_info off
%  ft_info off msgId
%  ft_info once
%  ft_info once msgId
%  ft_info on  backtrace
%  ft_info off backtrace
%  ft_info on  verbose
%  ft_info off verbose
%
%  ft_info query      % shows the status of all notifications
%  ft_info last       % shows the last notification
%  ft_info clear      % clears the status of all notifications
%  ft_info timeout 10 % sets the timeout (for 'once') to 10 seconds
%
% See also DEFAULTID, FT_ERROR, FT_WARNING, FT_NOTICE, FT_INFO, FT_DEBUG, ERROR, WARNING

% Copyright (C) 2012-2017, Robert Oostenveld, J?rn M. Horschig
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

global ft_default

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% determine who is calling

stack = dbstack('-completenames');
% st(1) is this function
% st(2) is the calling function
switch stack(2).name
  case 'ft_debug'
    level = 'debug';
  case 'ft_info'
    level = 'info';
  case 'ft_notice'
    level = 'notice';
  case 'ft_warning'
    level = 'warning';
  case 'ft_error'
    level = 'error';
  otherwise
    error('this function cannot be called from %s', stack(2).name);
end

% remove this function itself and the ft_xxx calling function
%stack = stack(3:end);
stack(1:2) = [];

% remove the non-FieldTrip functions from the path, these should not be part of the default message identifier
[v, p] = ft_version;
keep   = startsWith({stack.file}, p);
stack  = stack(keep);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% handle the defaults

if isfield(ft_default, 'notification') && isfield(ft_default.notification, level)
  % this would be error, warning, notice, info, debug
  s = ft_default.notification.(level);
else
  s = [];
end

% start with an empty state structure
if isempty(s)
  s = struct('identifier', {}, 'state', {}, 'timestamp', {});
end
ident = {s.identifier};

% set the default notification state
if ~any(strcmp(ident, 'all'))
  [s, ident] = setstate(s, 'all', 'on', ident);
end

% set the default backtrace state
defaultbacktrace = false;
if ~any(strcmp(ident, 'backtrace'))
  switch level
    case {'debug' 'info' 'notice'}
      [s, ident] = setstate(s, 'backtrace', 'off', ident);
    case 'warning'
      defaultbacktrace = true;
      t = warning('query', 'backtrace'); % get the default state
      [s, ident] = setstate(s, 'backtrace', t.state, ident);
    case 'error'
      [s, ident] = setstate(s, 'backtrace', 'on', ident);
  end % switch
end

% set the default verbose state
defaultverbose = false;
if ~any(strcmp(ident, 'verbose'))
  switch level
    case 'warning'
      defaultverbose = true;
      t = warning('query', 'verbose'); % get the default state
      [s, ident] = setstate(s, 'verbose', t.state, ident);
    otherwise
      [s, ident] = setstate(s, 'verbose', 'off', ident);
  end
end

% set the default timeout
if ~any(strcmp(ident, 'timeout'))
  [s, ident] = setstate(s, 'timeout', 60, ident);
end

% set the last notification to empty
if ~any(strcmp(ident, 'last'))
  state.message    = '';
  state.identifier = '';
  state.stack      = struct('file', {}, 'name', {}, 'line', {});
  [s, ident] = setstate(s, 'last', state, ident);
end

if strcmp(level, 'warning')
  ws = warning;
  % warnings should be on in general
  warning on
  % the backtrace is handled by this function
  warning off backtrace
  % the verbose message is handled by this function
  warning off verbose
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% set the state according to the input

if numel(varargin)==1 && (isstruct(varargin{1}) || isempty(varargin{1}))
  for i=1:numel(varargin{1})
    [s, ident] = setstate(s, varargin{1}(i).identifier, varargin{1}(i).state, ident);
  end
  ft_default.notification.(level) = s;
  return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% act according to the first input argument

if isempty(varargin)
  varargin{1} = 'query';
end

switch varargin{1}
  case 'on'
    if numel(varargin)>1
      msgId = varargin{2};
      % return the message state of this specific one
      varargout{1} = getreturnstate(s, msgId);
      % switch this specific item on
      s = setstate(s, msgId, 'on', ident);
      s = settimestamp(s, msgId, nan, ident);
      if strcmp(msgId, 'backtrace')
        defaultbacktrace = false;
      end
      if strcmp(msgId, 'verbose')
        defaultverbose = false;
      end
    else
      % return the message state of all
      varargout{1} = getreturnstate(s);
      % switch all on
      s = setstate(s, 'all', 'on', ident);
    end

  case 'off'
    if numel(varargin)>1
      msgId = varargin{2};
      % return the message state of this specific one
      varargout{1} = getreturnstate(s, msgId);
      % switch this specific item on
      s = setstate(s, msgId, 'off', ident);
      s = settimestamp(s, msgId, nan, ident);
      if strcmp(msgId, 'backtrace')
        defaultbacktrace = false;
      end
      if strcmp(msgId, 'verbose')
        defaultverbose = false;
      end
    else
      % return the message state of all
      varargout{1} = getreturnstate(s);
      % switch all off
      s = setstate(s, 'all', 'off', ident);
    end

  case 'once'
    if numel(varargin)>1
      msgId = varargin{2};
      % return the specific message state
      varargout{1} = getreturnstate(s, msgId);
      % switch a specific item to once
      s = setstate(s, msgId, 'once', ident);
    else
      % return the message state of all
      varargout{1} = getreturnstate(s);
      % switch all to once
      s = setstate(s, 'all', 'once', ident);
    end

  case 'timeout'
    % set the timeout, this is used for 'once'
    if ischar(varargin{2})
      s = setstate(s, 'timeout', str2double(varargin{2}), ident);
    else
      s = setstate(s, 'timeout', varargin{2}, ident);
    end

  case {'last' '-last'}
    % return the last notification
    varargout{1} = getstate(s, 'last', ident);

  case {'clear' '-clear'}
    % reset the notification system
    s = [];

  case {'query' '-query'}
    if numel(varargin)>1
      % select a specific item
      msgId = varargin{2};
      if ~any(strcmp(ident, msgId))
        error('Unknown setting or incorrect message identifier ''%s''.', msgId);
      end
      msgState = getstate(s, msgId, ident);
      if nargout
        varargout{1} = getreturnstate(s, msgId);
      elseif strcmp(msgId, 'verbose')
        if istrue(msgState)
          fprintf('%s output is verbose.\n', level);
        else
          fprintf('%s output is terse.\n', level);
        end
      elseif strcmp(msgId, 'backtrace')
        if istrue(msgState)
          fprintf('%s backtraces are enabled.\n', level);
        else
          fprintf('%s backtraces are disabled.\n', level);
        end
      else
        fprintf('The state of %s ''%s'' is ''%s''\n', level, msgId, msgState);
      end
    else
      % return all items
      r = getreturnstate(s);

      if nargout
        % return the state of all items
        varargout{1} = r;
      else
        % show the state of all items that are different from the default
        default = getstate(s, 'all', ident);
        fprintf('The default %s state is ''%s''.', level, default);
        r = r(~strcmp({r.state}, default));
        % don't show these
        r(strcmp('verbose',   {r.identifier})) = [];
        r(strcmp('backtrace', {r.identifier})) = [];
        if ~isempty(r)
          fprintf(' Items not set to the default are\n\n');
        end
        for i=1:numel(r)
          fprintf('  %4s  %s\n', r(i).state, r(i).identifier);
        end
        fprintf('\n');
      end
    end

  otherwise

    if nargout
      varargout{1} = getreturnstate(s);
    end

    % first input might be msgId
    if numel(varargin)>1 && ~isempty(regexp(varargin{1}, '^\w*:', 'once'))
      msgId = varargin{1};
      varargin = varargin(2:end); % shift them all by one
    else
      % use an automatically generated default identifier
      msgId = defaultId(stack);
    end

    % get the state for this notification, it will default to the 'all' state
    msgState = getstate(s, msgId, ident);

    % errors are always to be printed
    if strcmp(level, 'error')
      msgState = 'on';
    end

    if strcmp(msgState, 'once')
      timeout = getstate(s, 'timeout', ident);
      since   = elapsed(gettimestamp(s, msgId, ident));
      if (since>timeout)
        % the timeout has passed, update the timestamp and print the message
        s = settimestamp(s, msgId, tic, ident);
        msgState = 'on';
      else
        % the timeout has not yet passed, do not print the message
        msgState = 'off';
      end
    end

    if varargin{1}(end) =='\'
      % a lone trailing backslash, '\' , is not a valid control character
      varargin{1}(end+1) = '\';
    end

    % store the last notification
    state.message    = strtrim(sprintf(varargin{:})); % remove the trailing newline
    state.identifier = msgId;
    if ~isempty(stack)
      state.stack      = stack;
      s = setstate(s, 'last', state, ident);
    end

    if strcmp(msgState, 'on')

      if strcmp(level, 'error')
        % update the global variable, we won't return here after the error
        ft_default.notification.(level) = s;
        % the remainder is fully handled by the ERROR function
        if ~isempty(msgId)
          error(state);
        else
          error(rmfield(state, 'identifier'));
        end

      elseif strcmp(level, 'warning')
        if ~isempty(msgId)
          warning(msgId, varargin{:});
        else
          warning(varargin{:});
        end

      else
        % ensure there is a line end
        if isempty(regexp(varargin{1}, '\\n$', 'once')) % note the double \\
          varargin{1} = [varargin{1} '\n'];
        end
        fprintf(varargin{:});

      end % if level=error, warning or otherwise

      % decide whether the stack trace should be shown
      if istrue(getstate(s, 'backtrace', ident))
        for i=1:numel(stack)
          % show the deepest and lowest-level function first
          [p, f, x] = fileparts(stack(i).file);
          if isequal(f, stack(i).name)
            funname = stack(i).name;
          else
            funname = sprintf('%s>%s', f, stack(i).name);
          end
          filename = stack(i).file;
          if isempty(fileparts(filename))
            % it requires the full path
            filename = fullfile(pwd, filename);
          end
          if ft_platform_supports('html')
            fprintf(' In %s at line %d\n', filename, stack(i).line, funname, stack(i).line);
          else
            fprintf(' In ''%s'' at line %d\n', filename, stack(i).line);
          end
        end
        fprintf('\n');
      end

      % decide whether a verbose message should be shown
      if istrue(getstate(s, 'verbose', ident))
        if ~isempty(msgId)
          fprintf('Type "ft_%s off %s" to suppress this message.\n', level, msgId)
        else
          fprintf('Type "ft_%s off" to suppress this message.\n', level)
        end
      end

    end % if msgState is on

end % switch varargin{1}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% update the global variable
if defaultbacktrace && ~isempty(s)
  % do not keep the default, it will be determined again on the next call
  s = s(~strcmp({s.identifier}, 'backtrace'));
end
if defaultverbose && ~isempty(s)
  % do not keep the default, it will be determined again on the next call
  s = s(~strcmp({s.identifier}, 'verbose'));
end
ft_default.notification.(level) = s;

if strcmp(level, 'warning')
  % return to the original warning state
  warning(ws);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function r = getreturnstate(s, msgId)
if nargin<2
  r = s;
  ident = {r.identifier};
  % don't return these
  r(strcmp(ident, 'timeout')|strcmp(ident, 'last')) = [];
  % don't return the timestamps
  r = rmfield(r, 'timestamp');
else
  msgState = getstate(s, msgId, {s.identifier});
  r = struct('identifier', msgId, 'state', msgState);
end

function state = getstate(s, msgId, ident)
if nargin<3, ident = {s.identifier}; end
sel = find(strcmp(ident, msgId));
if numel(sel)==1
  state = s(sel).state;
  if isempty(state)
    state = getstate(s, 'all', ident);
  end
else
  state = getstate(s, 'all', ident);
end

function timestamp = gettimestamp(s, msgId, ident)
if nargin<3, ident = {s.identifier}; end
sel = find(strcmp(ident, msgId));
if numel(sel)==1
  timestamp = s(sel).timestamp;
else
  timestamp = nan;
end

function [s, ident] = setstate(s, msgId, state, ident)
if isempty(msgId), return; end % this happens from the command line
if nargin<4, ident = {s.identifier}; end
sel = find(strcmp(ident, msgId));
if numel(sel)==1
  s(sel).state = state;
else
  s(end+1).identifier = msgId;
  s(end  ).state      = state;
  s(end  ).timestamp  = nan;
  if nargout>1
    ident{end+1} = msgId;
  end
end

function s = settimestamp(s, msgId, timestamp, ident)
if isempty(msgId), return; end % this happens from the command line
if nargin<4, ident = {s.identifier}; end
sel = find(strcmp(ident, msgId));
if numel(sel)==1
  s(sel).timestamp = timestamp;
else
  s(end+1).identifier = msgId;
  s(end  ).state      = [];
  s(end  ).timestamp  = timestamp;
end

function t = elapsed(start)
if isempty(start) || isnan(start)
  t = inf;
else
  t = toc(start);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [rtnData, headData] = isfread (filename)
% isfread - will read the .ISF files produced by TEK TDS oscilloscopes
%
% useage:   [rtnData, headData] = isfread ('filename')
% where:    rtnData     - is a struct with the x and y data returned in
%                         rtnData.x and rtnData.y respectively
%           headData    - is a struct with all the header data returned in
%                         different fields
%           'filename'  - is a string with the name of the file to be
%                         extracted
%
% Example:  filename = 'TEK00000.ISF';
%           [data, header] = isfread (filename);
%           plot(data.x,data.y)
%
% The returned data is pre scaled using the information contained in the
% header of the ISF file.
%
% * Written 17/8/2004 by John Lipp - CCLRC Rutherford Appleton Laboratory *

fileID = fopen(filename,'r');

% read ASCII header
header_tmp = fread(fileID,267)';
header = char(header_tmp);

headData = parseHead(header);

% read binary data
inData = fread(fileID, headData.NR_PT, 'int16');

lowerXLimit = headData.XZERO;
upperXLimit = ((headData.NR_PT-1)* headData.XINCR  + headData.XZERO);
rtnData.x = [lowerXLimit:headData.XINCR:upperXLimit]';

rtnData.y = headData.YMULT*(inData-headData.YOFF);

% parseHead - used to parse and convert the header data into a more
% useful structure.
function headData = parseHead(header)

[headData.TYPE,rem] = strtok(header,':') ;

[headData.BYT_NR, rem] = getNextNum(rem);
[headData.BIT_NR, rem] = getNextNum(rem);
[headData.ENCDG,rem] = getNextStr(rem) ;
[headData.BN_FMT,rem] = getNextStr(rem) ;
[headData.BYT_OR,rem] = getNextStr(rem) ;
[headData.NR_PT, rem] = getNextNum(rem);
[headData.WFID,rem] = getNextStr(rem) ;
[headData.PT_FMT,rem] = getNextStr(rem) ;
[headData.XINCR, rem] = getNextNum(rem);
[headData.PT_OFF, rem] = getNextNum(rem);
[headData.XZERO, rem] = getNextNum(rem);
[headData.XUNIT,rem] = getNextStr(rem) ;
[headData.YMULT, rem] = getNextNum(rem);
[headData.YZERO, rem] = getNextNum(rem);
[headData.YOFF, rem] = getNextNum(rem);
[headData.YUNIT,rem] = getNextStr(rem) ;
[headData.CURVE,rem] = getNextStr(rem) ;

function [rtnNum, remStr] = getNextNum (string)
[junk,rem] = strtok(string,' ');
[tmp, remStr] = strtok(rem,';') ;
rtnNum = str2num(tmp);

function [rtnStr, remStr] = getNextStr(string)
[junk,rem] = strtok(string,' ');
rem = strtrim(rem);             % remove padding white space
[rtnStr,remStr] = strtok(rem,';') ;


text:
function isomapInfo = isomapDeconstruct(model)

% ISOMAPDECONSTRUCT break isomap in pieces for saving.
% FORMAT
% DESC takes an isomap model structure and breaks it into component
% parts for saving.
% ARG model : the model that needs to be saved.
% RETURN isomapInfo : a structure containing the other information
% from the isomap: what the sparse approximation is, what the inducing
% variables are.
%
% SEEALSO : isomapReconstruct
%
% COPYRIGHT : Neil D. Lawrence, 2009

% MLTOOLS

isomapInfo = model;
removeFields = {'Y'};

for i = 1:length(removeFields)
  if isfield(isomapInfo, removeFields{i})
    isomapInfo = rmfield(isomapInfo, removeFields{i});
  end
end


text:
%VideoCamera_fg Class to read from local video camera
%
% A concrete subclass of ImageSource that acquires images from a local
% camera using a simple open-source frame grabber interface.
%
% This class is not intended to be used directly, instead use the factory
% method VideoCamera.which will return an instance of this class if the interface
% is supported on your platform (Mac or Linux), for example
%
%         vid = VideoCamera.amera();
%
% Methods::
% grab    Aquire and return the next image
% size    Size of image
% close   Close the image source
% char    Convert the object parameters to human readable string
%
% See also ImageSource, AxisWebCamera, Movie.

% 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/>.

% mmread brings the whole movie into memory.  Not entirely sure what
% libavbin uses memory-wise, it takes a long time to "open" the file.

classdef VideoCamera_fg < ImageSource

    properties

        name;
    end

    properties(Constant=true)
        % set the operation codes for the camera driver
        CAMERA_OP_LIST = 1;
        CAMERA_OP_OPEN = 2;
        CAMERA_OP_CLOSE = 3;
        CAMERA_OP_START = 4;
        CAMERA_OP_STOP = 5;
        CAMERA_OP_ISRUNNING = 6;
        CAMERA_OP_GRAB = 7;
        CAMERA_OP_GET_PARAMS = 8;
        CAMERA_OP_SET_WIDTH_HEIGHT = 9;
        CAMERA_OP_GET_NAME = 10;
    end

    methods(Static)
        % class method to list the available video sources
        function list()
            framegrabber(VideoCamera_fg.CAMERA_OP_LIST);
        end
    end

    methods

        function m = VideoCamera_fg(varargin)
        %VideoCamera_fg.VideoCamera_fg Video camera constructor
        %
        % V = VideoCamera_fg.CAMERA, OPTIONS) is a VideoCamera_fg.object that acquires
        % images from the local video camera specified by the string CAMERA.
        %
        % If CAMERA is '?' a list of available cameras, and their
        % characteristics is displayed.
        %
        % Options::
        % 'uint8'          Return image with uint8 pixels (default)
        % 'float'          Return image with float pixels
        % 'double'         Return image with double precision pixels
        % 'grey'           Return greyscale image
        % 'gamma',G        Apply gamma correction with gamma=G
        % 'scale',S        Subsample the image by S in both directions.
        % 'resolution',S   Obtain an image of size S=[W H].
        % 'id',I           ID of camera
        %
        % Notes:
        % - The specified 'resolution' must match one that the camera is capable of,
        %   otherwise the result is not predictable.

            % invoke the superclass constructor and process common arguments
            m = m@ImageSource(varargin{:});

            % open the video source

            if length(varargin) == 0
                camera = -1;        % choose default
            else
                camera = varargin{1};
            end
            framegrabber(VideoCamera_fg.CAMERA_OP_OPEN, camera);

            for i=1:10
                k = framegrabber(VideoCamera_fg.CAMERA_OP_ISRUNNING);

                if k > 0
                    break;
                end
                pause(0.5);
            end

            % get the parameters of the video source
            sz = framegrabber(VideoCamera_fg.CAMERA_OP_GET_PARAMS);
            m.width = sz(1);
            m.height = sz(2);
            m.name = framegrabber(VideoCamera_fg.CAMERA_OP_GET_NAME);
        end

        function paramSet(v, a1)
        end

        % destructor
        function delete(m)
            framegrabber(VideoCamera_fg.CAMERA_OP_CLOSE);
        end

        function close(m)
        %VideoCamera_fg.close Close the image source
        %
        % V.close() closes the connection to the camera.
            framegrabber(VideoCamera_fg.CAMERA_OP_STOP);
            pause(2);
            framegrabber(VideoCamera_fg.CAMERA_OP_CLOSE);
        end

        function out = grab(m, opt)
        %VideoCamera_fg.grab Acquire image from the camera
        %
        % IM = V.grab() acquires an image from the camera.
        %
        % Notes::
        % - the function will block until the next frame is acquired.

            im = framegrabber(VideoCamera_fg.CAMERA_OP_GRAB);

            % apply options specified at construction time
            im = m.convert(im);
            m.width = numcols(im);
            m.height = numrows(im);

            if nargout > 0
                out = im;
            else
                idisp(im);
            end
        end

        function s = char(m)
        %VideoCamera_fg.char Convert to string
        %
        % V.char() is a string representing the state of the camera object in
        % human readable form.
            s = '';

            s = strvcat(s, sprintf('VideoCamera. %s %d x %d', m.name, m.width, m.height));

            % show constructor time options
            s2 = char@ImageSource(m);
            if ~isempty(s2)
                s = strvcat(s, strcat(' - ', s2));
            end
        end

    end
end

repo: machinevision-toolbox-matlab


text:
% ----- Global evaluation parameters -----
% Range of binning factors used for this evaluation.
binningFactors = [2, 3, 4];
% Range of compression levels (index 1 means uncoded; indices 2-5 mean
% coded with the corresponding QP setting using H.265/HEVC coder.
compressions   = [NaN, 10, 20, 30, 40];
% Range of number of frames for SR at the different binning factors.
numberOfFrames = [5, ...    % Number of frames for binning factor 2
                  11, ...   % Number of frames for binning factor 3
                  17];      % Number of frames for binning factor 4
% Index of the first reference frame for sliding window processing.
startReferenceFrame = 9;

% ----- Settings considered for this evaluation -----
% This is the index of the SR method that is evaluated.
sr_method           = 0:length(SRMethods);
% This is the index of the binning factor that is evaluated.
binning_val         = 1:3;
% This is the index of the compression setting.
compress_val        = 1;
% This is the index that indicates the number of input frames.
numberOfFrames_val  = 1;
% This is the index of the scene according to datasets.mat.
scenes_val          = 1:14;
% This is the index of the sliding window that is evaluated.
sliding_val         = 1;

author: thomas-koehler

repo: SupER


text:
function A = adjoint(X)
% ADJOINT Computes adjoint matrix
%
% A = ADJOINT(X)
%
% Brute-force implementation

[n,m] = size(X);
if n~=m
    error('Matrix must be square');
end

A = [];
if n == 1
    A = 1;
    return
end

% Ugly brute-force
for i = 1:n
    temp = [];
    noti = setdiff(1:n,i);
    for j = 1:n
        notj = setdiff(1:n,j);
        temp = [temp det(X(noti,notj))*((-1)^(i+j))];
    end
    A = [A;temp];
end
A = A';

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function inspect_bug2993

% WALLTIME 00:10:00
% MEM 1gb
% DEPENDENCY

cd(dccnpath('/home/common/matlab/fieldtrip/data/test/bug2993'));
elec_original = ft_read_sens('EEGMarkersLocalite.xml');

% elec_original = ft_determine_coordsys(elec_original);
elec_original.coordsys = 'lps';

% the electrodes should be rotated around the z-axis with 180 degrees to make them consistent
% with the MNI coordinate system (x=right, y=nose), which is expected by default in ft_prepare_layout

figure
ft_plot_sens(elec_original, 'label', 'on');
ft_plot_axes(elec_original);

%%
cfg = [];
cfg.elec = elec_original;
cfg.method = 'interactive';
elec_realigned1 = ft_electroderealign(cfg);

assert(~isequal(elec_realigned1.chanpos, elec_original.chanpos));

%%
cfg = [];
cfg.elec = elec_original;
cfg.method = 'interactive';
% this requires a headshape, but it can be empty
cfg.headshape = [];
cfg.headshape.pos  = zeros(0,3);
cfg.headshape.tri  = zeros(0,3);
cfg.headshape.unit = 'mm';
elec_realigned2 = ft_sensorrealign(cfg);

assert(~isequal(elec_realigned2.chanpos, elec_original.chanpos));

%%
cfg = [];
cfg.individual.elec = elec_original;
cfg.template.headshape = [];
cfg.template.headshape.pos  = randn(10,3);
cfg.template.headshape.tri  = zeros(0,3);
cfg.template.headshape.unit = 'mm';
cfg.method = 'interactive';
outcfg = ft_interactiverealign(cfg);

% this returns the transformation, not the updated electrode positions
assert(~isequal(outcfg.m, eye(4)));

%%

cfg = [];
cfg.elec = elec_original;
figure; ft_plot_layout(ft_prepare_layout(cfg));

cfg = [];
cfg.elec = elec_original;
cfg.rotate = 180;
figure; ft_plot_layout(ft_prepare_layout(cfg));

cfg = [];
cfg.elec = elec_realigned1;
figure; ft_plot_layout(ft_prepare_layout(cfg));


text:
function [Lmk,Obs] = reparametrizeLmk(Rob,Sen,Lmk,Obs,Opt)

%REPARAMETRIZELMK Reparametrize landmark.
%   [LMK,OBS] = REPARAMETRIZELMK(ROB,SEN,LMK,OBS,OPT) return the
%   reparametrized landmark and Obs structure.
%   The landmark is reparametrized only in the case of an inverse depth
%   point. A linearity test of cartesian point given inverse depth
%   point is performed.

%   Copyright 2009 David Marquez @ LAAS-CNRS.

global Map

switch Lmk.type

    case 'idpPnt'
        % we will convert from inverse-depth to euclidean.

        % Test for linearity:
        Ld = idpLinearTest(Rob,Sen,Lmk);

        if Ld < Opt.correct.linTestIdp

            % ranges
            ir = Lmk.state.r;  % idp
            er = ir(1:3);      % euclidean
            m  = Map.used;     % map

            % point coordinates
            idp     = Map.x(ir);   % idp
            [p,P_i] = idp2euc(idp);  % euclidean

            % map updates
            Map.x(er) = p;     % mean

            Map.P(er,m) = P_i * Map.P(ir,m); % co- and cross-variances
            Map.P(m,er) = Map.P(m,ir) * P_i';

            Map.used(Lmk.state.r(4:6)) = false; % used positions

            % Lmk and Obs updates
            Lmk.state.r = er;    % new range
            Lmk.type  = 'eucPnt'; % new type
            Obs.ltype = 'eucPnt'; % new type
        end

    case 'ahmPnt'
        % we will convert from anchored homogeneous to euclidean.

        % Test for linearity:
        Ld = ahmLinearTest(Rob,Sen,Lmk);

        if Ld < Opt.correct.linTestIdp

            % ranges
            ar  = Lmk.state.r;  % ahmPnt
            er  = ar(1:3);      % euclidean
            m   = Map.used;     % map

            % point coordinates
            ahm     = Map.x(ar);   % idp
            [p,P_a] = ahm2euc(ahm);  % euclidean

            % map updates
            Map.x(er)   = p;     % mean

            Map.P(er,m) = P_a * Map.P(ar,m); % co- and cross-variances
            Map.P(m,er) = Map.P(m,ar) * P_a';

            Map.used(Lmk.state.r(4:7)) = false; % used positions

            % Lmk and Obs updates
            Lmk.state.r = er;    % new range
            Lmk.type  = 'eucPnt'; % new type
            Obs.ltype = 'eucPnt'; % new type
        end

    case {'eucPnt'}
        % do nothing

    case {'hmgPnt','ahmPnt','plkLin','aplLin','idpLin','ahmLin','hmgLin'}
        % do nothing, by now <- probably add here something to do
        % Points should go to euclidean
        % Lines should go to some minimal representation (polar? 'plrLin')

    % case 'myLmk'
        % edit this 'myLmk' name to put your own landmark type
        % do something

    otherwise
        error('??? Unknown landmark type ''%s''.',Lmk.type)

end

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


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

% Vertices
vtx      = zeros(N,3);
vtx(:,1) = (-0.5*L:L/(N-1):0.5*L);

% Elements
elt = [(1:N-1) ; (2:N)]';

% Build mesh
mesh = msh(vtx,elt);
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [A, B] = mult_myForm(M, N, P, Q, C)
    % return myForm(M, N, C) * myForm(P, Q, C) = myForm(A, B, C);
    % -----------------------------------------------
    % Author: Tiep Vu, thv102@psu.edu, 4/8/2016
    %         (http://www.personal.psu.edu/thv102/)
    % -----------------------------------------------
    if nargin == 0
        d = 10;
        C = 10;
        M = rand(d, d);
        N = rand(d, d);
        P = rand(d, d);
        Q = rand(d, d);

    end
    %%
    A = M*P;
    B = M*Q + N*P + C*N*Q;
    % From here, we can see that if N == 0 => A = M*P; B = M*Q;
    % If Q == 0 => A = M*P; B = N*P.
    % toc;
    if nargin == 0
        M1 = myForm(A, B, C);
        M2 = myForm(M, N, C)*myForm(P, Q, C);
        norm(M1 - M2)
        pause
    end
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function [VV,FF,J] = extrude(V,F,varargin)
  % EXTRUDE Extrude a 2d mesh in the z direction by 1, connecting boundaries
  % apropriately
  %
  % [VV,FF] = extrude(V,F)
  % [VV,FF,J] = extrude(V,F,'ParmeterName',ParameterValue,...)
  %
  % Inputs:
  %   V  #V by 2 list of 2d vertex positions
  %   F  #F by 3 list of triangle indices into V
  %     or
  %   E  #E by 2 list of segment indices into V, CCW order
  %   Optional:
  %     'Cap' followed by {true} of false whether to put a *bottom cap*
  %     'Levels'  followed by number of "stacks" or "levels" {1}
  % Outputs:
  %   VV #VV by 3 list of 3d vertex positions
  %   FF  #FF by 3 list of triangle indices into VV
  %   J  #VV indices into V
  %

  cap = [];
  levels = 1;
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Cap','Levels'}, ...
    {'cap','levels'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end
  if isempty(cap)
    cap = size(F,2) > 2;
  end

  %assert(size(V,2) == 2, 'Vertices should be 2d');

  %% Copy top and bottom
  %  VV = [V 1+0*V(:,1);V 0*V(:,1)];
  switch size(F,2)
  case {3,4}
    % connect boundaries
    O = outline(F);
  case 2
    assert(~cap,'Cap should be false when input is curve');
    % Curve as input
    FF = [];
    O = F;
    % Use triangles
    F = zeros(0,3);
  end
  assert(size(unique(sort(O,2),'rows'),1) == size(O,1));

  % Rearrange vertices on bottom
  [VB,FB,OB,IM] = faces_first(V,F,O);
  n = size(V,1);
  JB = (1:n)';

  no = numel(unique(OB));

  VL = [ ...
      repmat(VB(1:max(OB(:)),:),[levels 1]) ...
        reshape(repmat(linspace(1-1/levels,0,levels),[no,1]),[no*levels 1])];
  JL = repmat((1:no)',levels,1);
  if cap
    VV = [VB ones(n,1);VL(1:no*(levels-1),:);VB zeros(n,1)];
    J = [JB;JL(1:no*(levels-1));JB];
  else
    VV = [VB ones(n,1);VL];
    J = [JB;JL];
  end
  FO = ones(2*size(OB,1)*levels,size(F,2));
  for level = 1:levels
    if level == 1
      off_t = 0;
    else
      off_t = n+(level-2)*no;
    end
    off_b = n+(level-1)*no;
    switch size(F,2)
    case 4
      FO((level-1)*size(OB,1)+(1:size(OB,1)),:) = [ ...
        off_t+OB(:,[2 1]) off_b+OB(:,1:2)];
    case 3
      FO((level-1)*2*size(OB,1)+(1:size(OB,1)*2),:) = [ ...
        off_b+OB(:,1) off_t+OB(:,[2 1]); ...
        off_t+OB(:,2) off_b+OB(:,1:2)];
    end
  end

  % remap so that faces on planes are consistent with input
  if cap
    FF = [FB;FO;n+(levels-1)*no+fliplr(FB)];
    RIM = 1:(n+(levels-1)*no+n);
    RIM(IM) = 1:n;
    RIM(n+IM) = n+(1:n);
  else
    FF = [FB;FO];
    RIM = 1:(n+levels*no);
    RIM(IM) = 1:n;
  end

  VV(RIM,:) = VV;
  J(RIM,:) = RIM(J);
  FF = RIM(FF);

end

author: alecjacobson

repo: gptoolbox


text:
function quad = fibonacci_lattice_q ( k, f )

%*****************************************************************************80
%
%% FIBONACCI_LATTICE_Q applies a Fibonacci lattice integration rule in 2D.
%
%  Discussion:
%
%    Because this is a standard lattice rule, it is really only suited
%    for functions which are periodic, of period 1, in both X and Y.
%
%    The related routines FIBONACCI_LATTICE_S and FIBONACCI_LATTICE_B
%    may be helpful in cases where the integrand does not satisfy this
%    requirement.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ian Sloan, Stephen Joe,
%    Lattice Methods for Multiple Integration,
%    Oxford, 1994,
%    ISBN: 0198534728,
%    LC: QA311.S56
%
%  Parameters:
%
%    Input, integer K, the index of the Fibonacci number to be used.
%    K must be at least 3.
%
%    Input, external real F, the name of the user-supplied routine
%    which evaluates the function, of the form:
%    function f ( dim_num, x )
%    integer dim_num
%    real f
%    real x(dim_num)
%    f = ...
%
%    Output, real QUAD, the estimated integral.
%
  dim_num = 2;

  quad = 0.0;

  m = fibonacci ( k );

  z(1) = 1;
  z(2) = fibonacci ( k - 1 );

  for j = 0 : m - 1
    x(1:dim_num) = mod ( j * z(1:dim_num) / m, 1.0 );
    quad = quad + f ( dim_num, x );
  end

  quad = quad / m;

  return
end


text:
% IALM + LMSVDS (Liu et al. 2012)
% process_video('RPCA', 'IALM_LMSVDS', 'dataset/demo.avi', 'output/demo_IALM_LMSVDS.avi');
[L,S] = inexact_alm_rpca_with_lmsvds(M);

author: andrewssobral

repo: lrslibrary


text:
classdef RWMOP1 < PROBLEM
% <multi> <real> <constrained>
% Pressure vessal problem

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

    methods
        %% Initialization
        function Setting(obj)
            obj.M        = 2;
            obj.D        = 4;
            obj.lower    = [0.51,0.51,10,10];
            obj.upper    = [99.49,99.49,200,200];
            obj.encoding = ones(1,obj.D);
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            x  = varargin{1};
            x1 = round(x(:,1));
            x2 = round(x(:,2));
            x3 = x(:,3);
            x4 = x(:,4);
            z1 = 0.0625*x1;
            z2 = 0.0625*x2;
            % Objective function
            f(:,1) = 1.7781.*z1.*x3.^2+0.6224.*z1.*x2.*x4+3.1661.*z1.^2.*x4+19.84.*z1.^2.*x3;
            f(:,2) = -pi.*x3.^2.*x4-(4/3).*pi.*x3.^3;
            % Constraints
            g(:,1) = 0.00954.*x3-z2;
            g(:,2) = 0.0193.*x3-z1;
            Population = SOLUTION(varargin{1},f,g,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
        %% Generate a point for hypervolume calculation
        function R = GetOptimum(obj,~)
            R = [3.5964885e+05  -7.3303829e+03];
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [OAoftrain,penalty] = svm_train_linear(traindata,trainlabel,rangeOfpenalty,incrOfpenalty)
% Using the training data and corresponding label to get the best penalty
% with respect to highest accuracy which is the result of 5-fold
% cross-validation(default).
% Using the grid search to look for the pair of parameters:(penalty)
% penalty lies in the range of rangeOfpenalty in the speed of incrOfpenalty

% svmtrain/svmpredict:
%       data:row represents # of samples, column represents # of features
%       label:column vector
% rangeOfpenalty:2-dims row vector
% incrOfpenalty:scaler on the basis of 2

i = 0;
row = size(rangeOfpenalty(1):incrOfpenalty:rangeOfpenalty(2),2);
OAoftrain = cell(1,row);
OAbest = 0;
penalty = 2^(rangeOfpenalty(1));
for iter_c = rangeOfpenalty(1):incrOfpenalty:rangeOfpenalty(2)
    i = i+1;
    cmd = ['-v 5 -t 0 -c ',sprintf('%.16f',2^(iter_c))];
    OAoftrain{i} = svmtrain(trainlabel,traindata,cmd);
    if OAoftrain{i} > OAbest
        penalty = 2^(iter_c);
        OAbest = OAoftrain{i};
    end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% SYNTAX:
% [data_dN, tInc, svs, nSV, tInc0] = hmrR_MotionCorrectPCArecurse(data_d, probe, mlActMan, mlActAuto, tIncMan, tMotion, tMask, STDEVthresh, AMPthresh, nSV, maxIter, turnon)
%
%
% UI NAME:
% Motion_Correct_PCA_Recurse
%
%
% DESCRIPTION:
% Identified motion artifacts in an input data matrix d. If any active
% data channel exhibits a signal change greater than STDEVthresh or
% AMPthresh, then a segment of data around that time point is marked as a
% motion artifact. Set maxIter=0 to skip this function.
%
%
% INPUTS:
% data_d: data matrix, timepoints x sd pairs
% sd: Source Detector Structure. The active data channels are indicated in
%     SD.MeasListAct.
% mlActMan: Cell array of vectors, one for each time base in data, specifying
%        active/inactive channels with 1 meaning active, 0 meaning inactive
% mlActAuto: Cell array of vectors, one for each time base in data, specifying
%        active/inactive channels with 1 meaning active, 0 meaning inactive
% tIncMan: Data that has been manually excluded. 0-excluded. 1-included.
%          Vector same length as d.
% tMotion: Check for signal change indicative of a motion artifact over
%     time range tMotion. Units of seconds. Typically tMotion=0.5.
% tMask: Mark data over +/- tMask seconds around the identified motion
%     artifact as a motion artifact. Units of seconds. Typicall tMask=1.
% STDEVthresh: If the signal d for any given active channel changes by more
%     that stdev_thresh * stdev(d) over the time interval tMotion, then
%     this time point is marked as a motion artifact. The standard deviation is
%     determined for each channel independently.
% AMPthresh: If the signal d for any given active channel changes by more
%     that AMPthresh over the time interval tMotion, then this time point
%     is marked as a motion artifact.
% nSV: This is the number of principal components to remove from the data.
%      If this number is less than 1, then the filter removes the first n
%      components of the data that removes a fraction of the variance
%      up to nSV. Yucel et al uses nSV=0.97.
% maxIter: Maximum number of iterations. Yucel et al uses maxIter=5;
%
%
% OUTPUTS:
% data_dN: This is the the motion corrected data.
% tInc: a vector of length time points with 1's indicating data included
%       and 0's indicate motion artifact AFTER correction of motion
%       artifacts
% svs: the singular values of the PCA for each iteration in each column
%      vector
% nSV: number of singular values removed from the data.
% tInc0: a vector of length time points with 1's indicating data included
%       and 0's indicate motion artifact BEFORE correction of motion
%       artifacts
%
% USAGE OPTIONS:
% Motion_Correct_PCA_Recurse:  [dod, tInc, svs, nSV, tInc0] = hmrR_MotionCorrectPCArecurse(dod, probe, mlActMan, mlActAuto, tIncMan, tMotion, tMask, STDEVthresh, AMPthresh, nSV, maxIter, turnon)
%
% PARAMETERS:
% tMotion: 0.5
% tMask: 1.0
% STDEVthresh: 20.0
% AMPthresh: 5.00
% nSV: 0.97
% maxIter: 5
% turnon: 1
%
% PREREQUISITES:
% Intensity_to_Delta_OD: dod = hmrR_Intensity2OD( intensity )
%

function [data_dN, tInc, svs, nSV, tInc0] = hmrR_MotionCorrectPCArecurse(data_d, probe, mlActMan,  mlActAuto, tIncMan, tMotion, tMask, STDEVthresh, AMPthresh, nSV, maxIter, turnon)

nBlks = length(data_d);
for iBlk=1:nBlks
    data_dN(iBlk) = DataClass(data_d(iBlk));
end

% Initialize output
tInc    = tIncMan;
svs     = cell(nBlks,1);
nSV     = repmat({nSV}, nBlks,1);
tInc0   = tIncMan;

% Check input args
if isempty(mlActMan)
    mlActMan = cell(nBlks,1);
end
if isempty(mlActAuto)
    mlActAuto = cell(nBlks,1);
end
if isempty(tIncMan)
    tIncMan = cell(nBlks,1);
end
if ~exist('turnon','var')
    turnon = 1;
end
if turnon==0
    return;
end

for iBlk=1:nBlks
    data_dN(iBlk) = DataClass(data_d(iBlk));

    tInc(iBlk) = hmrR_MotionArtifact(data_d(iBlk), probe, mlActMan(iBlk),  mlActAuto(iBlk), tIncMan(iBlk), tMotion, tMask, STDEVthresh, AMPthresh);

    tInc0{iBlk} = tInc{iBlk};

    ii=0;
    while length(find(tInc{iBlk}==0))>0 & ii<maxIter
        ii=ii+1;
        [data_dN(iBlk), svs_ii, nSV(iBlk)] = hmrR_MotionCorrectPCA(data_d(iBlk), mlActMan(iBlk), mlActAuto(iBlk), tIncMan(iBlk), tInc(iBlk), nSV(iBlk));
        tInc(iBlk) = hmrR_MotionArtifact(data_dN(iBlk), probe, mlActMan(iBlk), mlActAuto(iBlk), tIncMan(iBlk), tMotion, tMask, STDEVthresh, AMPthresh);
        data_d(iBlk).Copy(data_dN(iBlk));
        svs{iBlk}(:,ii) = svs_ii{1};
    end
end

author: BUNPC


text:
%% WGMGRATE Test convergence of multigrid methods for weak Galerkin method
%
% Reference
%
% An auxiliary space multigrid preconditioner for the weak Galerkin method.
% By Long Chen, Junping Wang, Yanqiu Wang, and Xiu Ye. Computers and
% Mathematics with Applications, 70(4):330?344 2015.
% doi:10.1016/j.camwa.2015.04.016

%% Options
clear variables;
close all;
option.maxIt = 4;
option.elemType = 'WG';
option.solver = 'mg';
option.smoothingstep = 2;
option.plotflag = 1;
option.rateflag = 0;
option.dispflag = 1;
colname = {'#Dof','Steps','Time'};

%% Example: circle mesh and Poisson equation
pde = sincosdata;
[node,elem] = circlemesh(0,0,1,0.25);
% showmesh(node,elem,'Facecolor','w');
bdFlag = setboundary(node,elem,'Dirichlet');
mesh = struct('node',node,'elem',elem,'bdFlag',bdFlag);

option.L0 = 1;
option.maxN = 3e5;
option.reducesystem = 0; % Solve the original system
[err,time,solver] = femPoisson(mesh,pde,option);
disptable(colname,solver.N,[],solver.itStep,[],solver.time,'%4.2g');

option.reducesystem = 1; % Solve the reduced system
[err,time,solver] = femPoisson(mesh,pde,option);
disptable(colname,solver.N,[],solver.itStep,[],solver.time,'%4.2g');

%% Example: square mesh and variable Poisson equation
[node,elem] = squaremesh([0,1,0,1],0.25);
% showmesh(node,elem,'Facecolor','w');
bdFlag = setboundary(node,elem,'Dirichlet');
mesh = struct('node',node,'elem',elem,'bdFlag',bdFlag);

option.L0 = 2;
option.dquadorder = 4;
pde = oscdiffdata;

option.reducesystem = 0; % Solve the original system
[err,time,solver] = femPoisson(mesh,pde,option);
disptable(colname,solver.N,[],solver.itStep,[],solver.time,'%4.2g');

option.reducesystem = 1; % Solve the reduced system
[err,time,solver] = femPoisson(mesh,pde,option);
disptable(colname,solver.N,[],solver.itStep,[],solver.time,'%4.2g');

%% Example: adaptive mesh and Poisson equation with less regularity
[node,elem] = squaremesh([-1,1,-1,1],1);
[node,elem] = delmesh(node,elem,'x>0 & y<0');
bdFlag = setboundary(node,elem,'Dirichlet');
mesh = struct('node',node,'elem',elem,'bdFlag',bdFlag);

pde = Lshapedata;
option.maxIt = 30;
option.maxN = 1e4;
option.printlevel = 0;
option.rateflag = 1;

option.reducesystem = 0; % Solve the original system
[err,time,solver] = afemPoisson(mesh,pde,option);
disptable(colname,solver.N,[],solver.itStep,[],solver.time,'%4.2g');

option.reducesystem = 1; % Solve the reduced system
[err,time,solver,eqn,node,elem] = afemPoisson(mesh,pde,option);
disptable(colname,solver.N,[],solver.itStep,[],solver.time,'%4.2g');
% figure; showmesh(node,elem,'Facecolor','w');

repo: ifem


text:
classdef SetSteeringModelAction < AbstractEventAction
    %SetSteeringModelAction Summary of this class goes here
    %   Detailed explanation goes here

    properties
        steeringModels SteeringModelsSet
    end

    %deprecated
    properties(Access=private)
        steeringModel(1,1) AbstractSteeringModel = RollPitchYawPolySteeringModel.getDefaultSteeringModel()
    end

    properties(Constant)
        emptyVarArr = AbstractOptimizationVariable.empty(0,1);
    end

    methods
        function obj = SetSteeringModelAction(steeringModel)
            obj.steeringModels = SteeringModelsSet();

            if(nargin > 0)
                obj.steeringModels.selectedModel = steeringModel;
            else
                obj.steeringModels.selectedModel = RollPitchYawPolySteeringModel.getDefaultSteeringModel();
            end

            obj.id = rand();
        end

%         function sm = get.steeringModel(obj)
%             sm = obj.steeringModel;
%             warning('Calling SetSteeringModelAction.steeringModel is deprecated.');
%         end
%
%         function set.steeringModel(obj,sm)
%             obj.steeringModel = sm;
%             warning('Calling SetThrottleModelAction.steeringModel is deprecated.');
%         end

        function newStateLogEntry = executeAction(obj, stateLogEntry)
            newStateLogEntry = stateLogEntry;
            newStateLogEntry.steeringModel = obj.steeringModels.selectedModel;
        end

        function initAction(obj, initialStateLogEntry)
            t0 = initialStateLogEntry.time;
            obj.steeringModels.selectedModel.setT0(t0);

            dcm = initialStateLogEntry.attitude.dcm;
            rVect = initialStateLogEntry.position;
            vVect = initialStateLogEntry.velocity;
            bodyInfo = initialStateLogEntry.centralBody;
            obj.steeringModels.selectedModel.setConstsFromDcmAndContinuitySettings(dcm, t0, rVect, vVect, bodyInfo);
        end

        function name = getName(obj)
            name = sprintf('Set Steering Model (%s)', obj.steeringModels.selectedModel.getTypeNameStr());
        end

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

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

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

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

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

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

        function tf = usesTankToTankConn(obj, tankToTank)
            tf = false;
        end

        function [tf, vars] = hasActiveOptimVar(obj)
            tf = false;
            vars = obj.emptyVarArr;

            optVar = obj.steeringModels.selectedModel.getExistingOptVar();
            if(not(isempty(optVar)))
                tf = any(optVar.getUseTfForVariable());
                vars(end+1) = optVar;
            end
        end
    end

    methods(Static)
        function addActionTf = openEditActionUI(action, lv)
            arguments
                action SetSteeringModelAction
                lv LaunchVehicle
            end

            addActionTf = action.steeringModels.openEditDialog(lv.lvdData, true);
        end

        function obj = loadobj(obj)
            arguments
                obj SetSteeringModelAction
            end

            if(isempty(obj.steeringModels))
                obj.steeringModels = SteeringModelsSet();
%                 obj.steeringModels.selectedModel = obj.steeringModel;

%             elseif(obj.steeringModels.selectedModel ~= obj.steeringModel)
%                 obj.steeringModels.selectedModel = obj.steeringModel;
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% This file is part of the following project:
% Oliver Parson, Siddhartha Ghosh, Mark Weal, Alex Rogers.
% Non-intrusive Load Monitoring using Prior Models of General Appliance Types.
% In: 26th AAAI Conference on Artificial Intelligence. Toronto, Canada. 2012.
% Code available for download: https://sites.google.com/site/oliparson/phd-work/research-files/aaai-2012-code.zip?attredirects=0
% Copyright: Oliver Parson et al., University of Southhampton, 2012.

% Modified by Romano Cicchetti, ETH Zurich, in the context of the NILM-Eval project

function [mpe] = my_viterbi_diff(bnet, evidence2, ignore_obs, lik_thres)

observed = ~isemptycell(evidence2);

T = length(evidence2);
mpe = evidence2(1:2,:);
d_states = bnet.node_sizes(bnet.dnodes_slice);
max_prob = zeros(d_states,T);
edges = zeros(d_states,T);

initial = struct(bnet.CPD{1});
emission2 = struct(bnet.CPD{2});
transition = struct(bnet.CPD{3});
emission = struct(bnet.CPD{4});

% prob of transitions
trans = transition.CPT;
diffMeans = permute(emission.mean, [2,3,1]);
diffCovs = sqrt(permute(emission.cov, [3,4,1,2]));
absoluteMeans = emission2.mean(:)';
absoluteCovs = emission2.cov(:)';

% max product forward pass
for t=1:T,
    % prob of states at t-1
    if t>1
        chain = max_prob(:,t-1);
    end

    % prob of diff emission
    emit = normpdf(evidence2{2,t}, diffMeans, diffCovs);

    % prob of absolute emission
    emit2 = normcdf(evidence2{3,t}, absoluteMeans, absoluteCovs);
    emit2 = emit2 / sum(emit2(:));

    % ignore observations of low probability
    if ignore_obs && sum(emit(:)) < lik_thres
        if t==1
            max_prob(:,t) = log(initial.CPT);
        else
            % product
            product = repmat(chain,[1,d_states]) + log(trans) + log(repmat(emit2,[d_states,1]));
            % max
            [max_prob(:,t), edges(:,t)] = max(product);
        end
    else
        if t==1
            max_prob(:,t) = log(initial.CPT);
        else
            %normalise emission probabilites
            emit = emit / sum(emit(:));
            % product
            product = repmat(chain,[1,d_states]) + log(trans) + log(emit) + log(repmat(emit2,[d_states,1]));
            % max
            [max_prob(:,t), edges(:,t)] = max(product);
        end
    end

    %normalise
    max_prob(:,t) = max_prob(:,t) - max(max_prob(:,t));

    a = max_prob(:,t);
    max_prob((isinf(a)),t) = min(a(~isinf(a))) - 100;

    if any(isinf(max_prob(:,t)))
        1;
    end

    % if observed then prob = 1
    if observed(1,t)
        max_prob(:,t) = 0;
        max_prob(evidence2{1,t},t) = 1;
    end
end

% backward pass
[~, mpe{1,T}] = max(max_prob(:,T));
for t=T:-1:2,
    mpe{1,t-1} = edges(mpe{1,t},t);
end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
 function s = sum(ob)
%function s = sum(ob)
%	"sum" method for Gblock object

s = sum(ob.base);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function nextpop = extractPop(opt, combinepop)
% Function: nextpop = extractPop(opt, combinepop)
% Description: Extract the best n individuals in 'combinepop'(population
%   size is 2n).
%
%         LSSSSWC, NWPU
%    Revision: 1.1  Data: 2011-07-12
%*************************************************************************

popsize = length(combinepop) / 2;
nextpop = combinepop(1:popsize);    %just for initializing

rankVector = vertcat(combinepop.rank);

n = 0;          % individuals number of next population
rank = 1;       % current rank number
idx = find(rankVector == rank);
numInd = length(idx);       % number of individuals in current front
while( n + numInd <= popsize )
    nextpop( n+1 : n+numInd ) = combinepop( idx );

    n = n + numInd;
    rank = rank + 1;

    idx = find(rankVector == rank);
    numInd = length(idx);
end

% If the number of individuals in the next front plus the number of individuals
% in the current front is greater than the population size, then select the
% best individuals by corwding distance(NSGA-II) or preference distance(R-NSGA-II).
if( n < popsize )
    if(~isempty(opt.refPoints))
        prefDistance = vertcat(combinepop(idx).prefDistance);
        prefDistance = [prefDistance, idx];
        prefDistance = sortrows( prefDistance, 1);
        idxSelect  = prefDistance( 1:popsize-n, 2);       % Select the individuals with smallest preference distance
        nextpop(n+1 : popsize) = combinepop(idxSelect);
    else
        distance = vertcat(combinepop(idx).distance);
        distance = [distance, idx];
        distance = flipud( sortrows( distance, 1) );      % Sort the individuals in descending order of crowding distance in the front.
        idxSelect  = distance( 1:popsize-n, 2);           % Select the (popsize-n) individuals with largest crowding distance.
        nextpop(n+1 : popsize) = combinepop(idxSelect);
    end
end

author: Eric-Bradford

repo: TS-EMO


text:
function  [fg] = nfgCreateConTrackPDB(phantomDir, fStepSizeMm, nMinNodes, nMaxNodes, bSanityCheck)
%Create ConTrack PDB
%
%   [fg] = nfgCreateSTTPDB(phantomDir, fStepSizeMm, nMinNodes, nMaxNodes,
%   bSanityCheck)
%
% AUTHORS:
% 2009.08.05 : AJS wrote it.
%
% NOTES:

if ieNotDefined('bSanityCheck'); bSanityCheck=0; end

% Directories
binDir = nfgGetName('binDir',phantomDir);
% Input Files
dtFile = nfgGetName('dtFile',phantomDir);
wmROIFile = nfgGetName('wmROIFile',phantomDir);
% Output Files
sttPDBFile = nfgGetName('sttPDBFile',phantomDir);
ctrparamsFile = nfgGetName('ctrparamsFile',phantomDir);
gold_sttPDBFile = nfgGetName('gold_sttPDBFile',phantomDir);
goldPDBFile = nfgGetName('goldPDBFile',phantomDir);

return;

author: vistalab

repo: vistasoft