text:
function [vert,conn,tria,tnum] = tridiv2(varargin)
%TRIDIV2 "quadtree" refinement for 2-simplex triangulations.
%   [VERT,EDGE,TRIA,TNUM] = TRIDIV2(VERT,EDGE,TRIA,TNUM) re-
%   turns a globally refined triangulation, in which all ed-
%   ges are bisected about their midpoints. Such refinement
%   splits each triangle into four new sub-triangles accord-
%   ing to a shape-preserving scheme.
%
%   VERT is a V-by-2 array of XY coordinates in the triangu-
%   lation, EDGE is an array of constrained edges, TRIA is a
%   T-by-3 array of triangles, and TNUM is a T-by-1 array of
%   part indices. Each row of TRIA and EDGE define an eleme-
%   nt. VERT(TRIA(II,1),:), VERT(TRIA(II,2),:) and VERT(TRIA
%   (II,3),:) are the coordinates of the II-TH triangle. The
%   edges in EDGE are defined in a similar manner. NUM is an
%   array of part indexing, such that TNUM(II) is the index
%   of the part in which the II-TH triangle resides.
%
%   [VERT,EDGE,TRIA,TNUM] = TRIDIV2(... ,TDIV) returns a se-
%   lectively refined mesh, where TDIV is a T-by-1 logical
%   array, with TDIV(KK) = TRUE if TRIA(KK,:) is to be refi-
%   ned. Such triangles are refined using the four-way split
%   described above. Additionally, a "halo" of adjacent tri-
%   angles are also refined to preseve compatibility of the
%   mesh. Such triangles are refined using a two-way bisect-
%   ion type scheme.

%   See also REFINE2, SMOOTH2

%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 29/01/2017

%---------------------------------------------- extract args
    vert = []; conn = []; tria = []; tnum = [] ;
    tdiv = [];

    if (nargin>=+1), vert = varargin{1}; end
    if (nargin>=+2), conn = varargin{2}; end
    if (nargin>=+3), tria = varargin{3}; end
    if (nargin>=+4), tnum = varargin{4}; end
    if (nargin>=+5), tdiv = varargin{5}; end

%---------------------------------------------- default arg.
    if (isempty(tnum))
        tnum = ones(size(tria,1),1) ;
    end
    if (isempty(tdiv))
        tdiv = true(size(tria,1),1) ;
    end

%---------------------------------------------- basic checks
    if ( ~isnumeric(vert) || ...
         ~isnumeric(conn) || ...
         ~isnumeric(tria) || ...
         ~isnumeric(tnum) || ...
         ~islogical(tdiv) )
        error('tridiv2:incorrectInputClass' , ...
            'Incorrect input class.') ;
    end

%---------------------------------------------- basic checks
    tnum = tnum(:) ; tdiv = tdiv(:) ;
    if (ndims(vert) ~= +2 || ...
        ndims(conn) ~= +2 || ...
        ndims(tria) ~= +2 || ...
        ndims(tnum) ~= +2 || ...
        ndims(tdiv) ~= +2 )
        error('tridiv2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end
    if (size(vert,2)~= +2 || ...
        size(conn,2)~= +2 || ...
        size(tria,2)~= +3 || ...
        size(tnum,2)~= +1 || ...
        size(tria,1)~= size(tnum,1) )
        error('tridiv2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

    nvrt = size(vert,1) ;

%---------------------------------------------- basic checks
    if (min(min(conn(:,1:2))) < +1 || ...
            max(max(conn(:,1:2))) > nvrt )
        error('tridiv2:invalidInputs', ...
            'Invalid EDGE input array.') ;
    end

    if (min(min(tria(:,1:3))) < +1 || ...
            max(max(tria(:,1:3))) > nvrt )
        error('tridiv2:invalidInputs', ...
            'Invalid TRIA input array.') ;
    end

%------------------------------ assemble extended adj. info.
   [edge,tria] = tricon2(tria,conn) ;

    ediv = false(size(edge,1),1) ;
    ediv(tria(tdiv,4:6))  = true ;

    snum = length(find(ediv));

    while (true)

    %-------------------------- tria's with >= 2 edge splits
        div3 = sum(double( ...
            ediv(tria(:,4:6))),2)>=2;

    %-------------------------- expand onto adj. edge splits
        ediv(tria(div3,4:6)) = true ;

        snew = length(find(ediv)) ;

        if (snew == snum), break; end

        snum = snew ;

    end

%------------------------------ tria's with == 1 edge splits
    div1 = sum( ...
    double( ediv(tria(:,4:6))),2)==1;

%------------------------------ indexing for mid-point vert.
    ivec = zeros(size(edge,1),1);
    ivec(ediv) = ...
        (1:snum)' + size(vert,1);

%------------------------------------------ update vert. set
    emid = vert(edge(ediv,1),:) ...
         + vert(edge(ediv,2),:) ;
    vert =[vert ; emid * 0.5] ;

%------------------------------------------ update edge. set
    [cvec,eloc] = ...
        setset2(conn,edge(ediv, 1:2)) ;

    epos = ivec(ediv);
    epos = epos(eloc(eloc>0)) ;

    conn = [conn(~cvec,1:2) ; ...
        conn(cvec,1) , epos ; ...
        conn(cvec,2) , epos ] ;

%------------------------------------ push 1-to-4 refinement
    tr41 = ones(length(find(div3)),3) ;
    tn41 = ones(length(find(div3)),1) ;
    tn41(:,1) = tnum(div3,1);
    tr41(:,1) = tria(div3,1);
    tr41(:,2) = ivec(tria(div3,4));
    tr41(:,3) = ivec(tria(div3,6));

    tr42 = ones(length(find(div3)),3) ;
    tn42 = ones(length(find(div3)),1) ;
    tn42(:,1) = tnum(div3,1);
    tr42(:,1) = ivec(tria(div3,4));
    tr42(:,2) = tria(div3,2);
    tr42(:,3) = ivec(tria(div3,5));

    tr43 = ones(length(find(div3)),3) ;
    tn43 = ones(length(find(div3)),1) ;
    tn43(:,1) = tnum(div3,1);
    tr43(:,1) = ivec(tria(div3,6));
    tr43(:,2) = ivec(tria(div3,5));
    tr43(:,3) = tria(div3,3);

    tr44 = ones(length(find(div3)),3) ;
    tn44 = ones(length(find(div3)),1) ;
    tn44(:,1) = tnum(div3,1);
    tr44(:,1) = ivec(tria(div3,6));
    tr44(:,2) = ivec(tria(div3,4));
    tr44(:,3) = ivec(tria(div3,5));

%----------------------- push 1-to-2 refinement about edge 1
    tvec = false(size(tria,1), 1) ;
    tvec(ediv(tria(:,4))&div1) = true ;

    tr21 = ones(length(find(tvec)),3) ;
    tn21 = ones(length(find(tvec)),1) ;
    tn21(:,1) = tnum(tvec,1);
    tr21(:,1) = ivec(tria(tvec,4));
    tr21(:,2) = tria(tvec,3);
    tr21(:,3) = tria(tvec,1);

    tr22 = ones(length(find(tvec)),3) ;
    tn22 = ones(length(find(tvec)),1) ;
    tn22(:,1) = tnum(tvec,1);
    tr22(:,1) = ivec(tria(tvec,4));
    tr22(:,2) = tria(tvec,2);
    tr22(:,3) = tria(tvec,3);

%----------------------- push 1-to-2 refinement about edge 2
    tvec = false(size(tria,1), 1) ;
    tvec(ediv(tria(:,5))&div1) = true ;

    tr23 = ones(length(find(tvec)),3) ;
    tn23 = ones(length(find(tvec)),1) ;
    tn23(:,1) = tnum(tvec,1);
    tr23(:,1) = ivec(tria(tvec,5));
    tr23(:,2) = tria(tvec,1);
    tr23(:,3) = tria(tvec,2);

    tr24 = ones(length(find(tvec)),3) ;
    tn24 = ones(length(find(tvec)),1) ;
    tn24(:,1) = tnum(tvec,1);
    tr24(:,1) = ivec(tria(tvec,5));
    tr24(:,2) = tria(tvec,3);
    tr24(:,3) = tria(tvec,1);

%----------------------- push 1-to-2 refinement about edge 3
    tvec = false(size(tria,1), 1) ;
    tvec(ediv(tria(:,6))&div1) = true ;

    tr25 = ones(length(find(tvec)),3) ;
    tn25 = ones(length(find(tvec)),1) ;
    tn25(:,1) = tnum(tvec,1);
    tr25(:,1) = ivec(tria(tvec,6));
    tr25(:,2) = tria(tvec,2);
    tr25(:,3) = tria(tvec,3);

    tr26 = ones(length(find(tvec)),3) ;
    tn26 = ones(length(find(tvec)),1) ;
    tn26(:,1) = tnum(tvec,1);
    tr26(:,1) = ivec(tria(tvec,6));
    tr26(:,2) = tria(tvec,1);
    tr26(:,3) = tria(tvec,2);

%------------------------------------------ update tria. set
    tria = [tria(~div1&~div3,1:3) ; ...
        tr41 ; tr42 ; tr43 ; tr44 ; ...
        tr21 ; tr22 ;
        tr23 ; tr24 ;
        tr25 ; tr26 ] ;

    tnum = [tnum(~div1&~div3,1:1) ; ...
        tn41 ; tn42 ; tn43 ; tn44 ; ...
        tn21 ; tn22 ;
        tn23 ; tn24 ;
        tn25 ; tn26 ] ;

end

author: CHLNDDEV


text:
function [xc,good,bad,type] = cornerfinder(xt,I,wintx,winty,wx2,wy2);

%[xc] = cornerfinder(xt,I);
%
%Finds the sub-pixel corners on the image I with initial guess xt
%xt and xc are 2xN matrices. The first component is the x coordinate
%(horizontal) and the second component is the y coordinate (vertical)
%
%Based on Harris corner finder method
%
%Finds corners to a precision below .1 pixel!
%Oct. 14th, 1997 - UPDATED to work with vertical and horizontal edges as well!!!
%Sept 1998 - UPDATED to handle diverged points: we keep the original points
%good is a binary vector indicating wether a feature point has been properly
%found.
%
%Add a zero zone of size wx2,wy2
%July 15th, 1999 - Bug on the mask building... fixed + change to Gaussian mask with higher
%resolution and larger number of iterations.

% California Institute of Technology
% (c) Jean-Yves Bouguet -- Oct. 14th, 1997

line_feat = 1; % set to 1 to allow for extraction of line features.

xt = xt';
xt = fliplr(xt);

if nargin < 4,
   winty = 5;
   if nargin < 3,
      wintx = 5;
   end;
end;

if nargin < 6,
   wx2 = -1;
   wy2 = -1;
end;

%mask = ones(2*wintx+1,2*winty+1);
mask = exp(-((-wintx:wintx)'/(wintx)).^2) * exp(-((-winty:winty)/(winty)).^2);

% another mask:
[X,Y] = meshgrid(-winty:winty,-wintx:wintx);
mask2 = X.^2 + Y.^2;
mask2(wintx+1,winty+1) = 1;
mask2 = 1./mask2;
%mask - mask2;

if (wx2>0) & (wy2>0),
   if ((wintx - wx2)>=2)&((winty - wy2)>=2),
      mask(wintx+1-wx2:wintx+1+wx2,winty+1-wy2:winty+1+wy2)= zeros(2*wx2+1,2*wy2+1);
   end;
end;

offx = [-wintx:wintx]'*ones(1,2*winty+1);
offy = ones(2*wintx+1,1)*[-winty:winty];

resolution = 0.005;

MaxIter = 10;

[nx,ny] = size(I);
N = size(xt,1);

xc = xt; % first guess... they don't move !!!

type = zeros(1,N);

for i=1:N,

   v_extra = resolution + 1; 		% just larger than resolution

   compt = 0; 				% no iteration yet

   while (norm(v_extra) > resolution) & (compt<MaxIter),

      cIx = xc(i,1); 			%
      cIy = xc(i,2); 			% Coords. of the point
      crIx = round(cIx); 		% on the initial image
      crIy = round(cIy); 		%
      itIx = cIx - crIx; 		% Coefficients
      itIy = cIy - crIy; 		% to compute
      if itIx > 0, 			% the sub pixel
	 vIx = [itIx 1-itIx 0]'; 	% accuracy.
      else
	 vIx = [0 1+itIx -itIx]';
      end;
      if itIy > 0,
	 vIy = [itIy 1-itIy 0];
      else
	 vIy = [0 1+itIy -itIy];
      end;

      % What if the sub image is not in?

      if (crIx-wintx-2 < 1), xmin=1; xmax = 2*wintx+5;
      elseif (crIx+wintx+2 > nx), xmax = nx; xmin = nx-2*wintx-4;
      else
	 xmin = crIx-wintx-2; xmax = crIx+wintx+2;
      end;

      if (crIy-winty-2 < 1), ymin=1; ymax = 2*winty+5;
      elseif (crIy+winty+2 > ny), ymax = ny; ymin = ny-2*winty-4;
      else
	 ymin = crIy-winty-2; ymax = crIy+winty+2;
      end;

      SI = I(xmin:xmax,ymin:ymax); % The necessary neighborhood
      SI = conv2(conv2(SI,vIx,'same'),vIy,'same');
      SI = SI(2:2*wintx+4,2:2*winty+4); % The subpixel interpolated neighborhood
      [gy,gx] = gradient(SI); 		% The gradient image
      gx = gx(2:2*wintx+2,2:2*winty+2); % extraction of the useful parts only
      gy = gy(2:2*wintx+2,2:2*winty+2); % of the gradients

      px = cIx + offx;
      py = cIy + offy;

      gxx = gx .* gx .* mask;
      gyy = gy .* gy .* mask;
      gxy = gx .* gy .* mask;

      bb = [sum(sum(gxx .* px + gxy .* py)); sum(sum(gxy .* px + gyy .* py))];

      a = sum(sum(gxx));
      b = sum(sum(gxy));
      c = sum(sum(gyy));

      dt = a*c - b^2;

      xc2 = [c*bb(1)-b*bb(2) a*bb(2)-b*bb(1)]/dt;

      %keyboard;

      if line_feat,

	 G = [a b;b c];
	 [U,S,V]  = svd(G);

	 %keyboard;

	 % If non-invertible, then project the point onto the edge orthogonal:

	 if (S(1,1)/S(2,2) > 50),
	    % projection operation:
	    xc2 = xc2 + sum((xc(i,:)-xc2).*(V(:,2)'))*V(:,2)';
	    type(i) = 1;
	 end;

      end;

      %keyboard;

%      G = [a b;b c];
%      [U,S,V]  = svd(G);

%      if S(1,1)/S(2,2) > 150,
%	 bb2 = U'*bb;
%	 xc2 = (V*[bb2(1)/S(1,1) ;0])';
%      else
%	 xc2 = [c*bb(1)-b*bb(2) a*bb(2)-b*bb(1)]/dt;
%      end;

      %if (abs(a)> 50*abs(c)),
%	 xc2 = [(c*bb(1)-b*bb(2))/dt xc(i,2)];
%      elseif (abs(c)> 50*abs(a))
%	 xc2 = [xc(i,1) (a*bb(2)-b*bb(1))/dt];
%      else
%	 xc2 = [c*bb(1)-b*bb(2) a*bb(2)-b*bb(1)]/dt;
%      end;

      %keyboard;

      v_extra = xc(i,:) - xc2;

      xc(i,:) = xc2;

%      keyboard;

      compt = compt + 1;

   end
end;

% check for points that diverge:

delta_x = xc(:,1) - xt(:,1);
delta_y = xc(:,2) - xt(:,2);

%keyboard;

bad = (abs(delta_x) > wintx) | (abs(delta_y) > winty);
good = ~bad;
in_bad = find(bad);

% For the diverged points, keep the original guesses:

xc(in_bad,:) = xt(in_bad,:);

xc = fliplr(xc);
xc = xc';

bad = bad';
good = good';

repo: ekfmonoslam


text:
function [logp, yhat, res] = tapas_softmax_mu3(r, infStates, ptrans)
% Calculates the log-probability of responses under the softmax model
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2017-2019 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Predictions or posteriors?
pop = 1; % Default: predictions
if r.c_obs.predorpost == 2
    pop = 3; % Alternative: posteriors
end

% Initialize returned log-probabilities, predictions,
% and residuals as NaNs so that NaN is returned for all
% irregualar trials
n = size(infStates,1);
logp = NaN(n,1);
yhat = NaN(n,1);
res  = NaN(n,1);

% Assumed structure of infStates:
% dim 1: time (ie, input sequence number)
% dim 2: HGF level
% dim 3: choice number
% dim 4: 1: muhat, 2: sahat, 3: mu, 4: sa

% Number of choices
nc = size(infStates,3);

% Belief trajectories at 1st level
states = squeeze(infStates(:,1,:,pop));

% Log-volatility trajectory
mu3 = squeeze(infStates(:,3,1,3));

% Responses
y = r.y(:,1);

% Weed irregular trials out from inferred states and responses
states(r.irr,:) = [];
mu3(r.irr) = [];
y(r.irr) = [];

% Inverse decision temperature
be = exp(-mu3);
be = repmat(be,1,nc);

% Partition functions
Z = sum(exp(be.*states),2);
Z = repmat(Z,1,nc);

% Softmax probabilities
prob = exp(be.*states)./Z;

% Extract probabilities of chosen options
probc = prob(sub2ind(size(prob), 1:length(y), y'));

% Calculate log-probabilities for non-irregular trials
reg = ~ismember(1:n,r.irr);
logp(reg) = log(probc);
yhat(reg) = probc;
res(reg) = -log(probc);

end

author: translationalneuromodeling

repo: tapas


text:
% Copyright (C) 2017-2018 Titus Cieslewski, RPG, University of Zurich,
%   Switzerland
%   You can contact the author at <titus at ifi dot uzh dot ch>
% Copyright (C) 2017-2018 Siddharth Choudhary, College of Computing,
%   Georgia Institute of Technology, Atlanta, GA, USA
% Copyright (C) 2017-2018 Davide Scaramuzza, RPG, University of Zurich,
%   Switzerland
%
% This file is part of dslam_open.
%
% dslam_open 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.
%
% dslam_open 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 dslam_open. If not, see <http://www.gnu.org/licenses/>.

function writeG2oConstraint(file_id, from_id, to_id, T_from_to, covariance)
% Adapted from Luca Carlone's
% https://bitbucket.org/lucacarlone/pgo3d-duality-opencode/src/ebb6e1b8cebaad7f2aaf581b1d0c0bad737faebb/lib/writeG2oDataset3D.m?at=master&fileviewer=file-view-default

dt = T_from_to(1:3, 4);
dx = dt(1); dy = dt(2); dz = dt(3);
dR = T_from_to(1:3, 1:3);
dq = rot2quat(dR);
assert(all(imag(dq)==0));
if norm(dq)>1e-3
    dq = dq/norm(dq);
else
    error('norm close to zero for unit quaternion (2)')
end
dqw = dq(1); dqx = dq(2); dqy = dq(3); dqz = dq(4);
I = covariance;

fprintf(file_id,'EDGE_SE3:QUAT %d %d   %f %f %f   %.7f %.7f %.7f %.7f   %f %f %f %f %f %f   %f %f %f %f %f   %f %f %f %f   %f %f %f   %f %f   %f\n', ...
    from_id, to_id, dx, dy, dz, dqx, dqy, dqz, dqw, ...
    I(1,1), I(1,2), I(1,3), I(1,4), I(1,5), I(1,6), ...
    I(2,2), I(2,3), I(2,4), I(2,5), I(2,6), ...
    I(3,3), I(3,4), I(3,5), I(3,6), ...
    I(4,4), I(4,5), I(4,6), ...
    I(5,5), I(5,6), ...
    I(6,6));

end

author: uzh-rpg


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function y = sprice_a(a, b, r, n, f, k, t,cp)
% sabr prices using the risk neutral density psabr
% and integrating this density with respect to the payoff

nl = length(k);
y = ones(1,nl);
if (cp == 1)
    % call
    for j = 1:nl
    F = @(x) (x-k(j)) .* psabr(a, b, r, n, f, x, t);
    y(j) = quad(F,k(j),1);
    end
else
    % put
    for j = 1:nl
    F = @(x) (k(j)-x) .* psabr(a, b, r, n, f, x, t);
    y(j) = quad(F,0.0001,k(j));
    end
end

end

author: Sable

repo: mcbench-benchmarks


text:
clear all

Nsamples = 41500;
EulerSaved = zeros(Nsamples, 3);

dt = 0.01;

for k=1:Nsamples
  [p q r] = GetGyro();
  [ax ay] = GetAccel();

  [phi_a theta_a] = EulerAccel(ax, ay);

  [phi theta psi] = EulerUKF([phi_a theta_a]', [p q r], dt);

  EulerSaved(k, :) = [ phi theta psi ];
end

PhiSaved   = EulerSaved(:, 1) * 180/pi;
ThetaSaved = EulerSaved(:, 2) * 180/pi;
PsiSaved   = EulerSaved(:, 3) * 180/pi;

t = 0:dt:Nsamples*dt-dt;

figure
plot(t, PhiSaved)

figure
plot(t, ThetaSaved)

figure
plot(t, PsiSaved)

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
% io_loadRFwaveform.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% [RF_struct]=io_loadRFwaveform(filename,type,off_res);
%
% DESCRIPTION:
% Initialize an RF pulse structure to contain an
% RF Pulse waveform as well as its accompanying header
% information.  This function finds the time-bandwidth
% product (tbw) and the time-w1 product (tw1) of the pulse,
% and stores this information in the header fields of the
% output RF structure.
%
% INPUTS:
% filename  = filename of RF pulse waveform text file.  Can be in Siemens
%             format (.pta), Varian/Agilent format (.RF), Bruker format
%             (.inv, .ref or .exc) or a plain text file (.txt, with two,
%             three or four columns (amplitude, phase, time-step
%             (optional),and gradient strength (optional).  Filename can
%             also be the name of a two, three or four column matlab vector
%             specifing the phase, amplitude, time vectors (optionally) and
%             gradient (optionally) of an RF pulse waveform.  Note about
%             Units: RF amplitude waveform [arbitrary units], B1 intensity
%             [kHz],  phase waveform [degrees], time step vector [arbitrary
%             units],  gradient waveform [G/cm].
% type      = Excitation ('exc'), Refocusing ('ref') or Inversion ('inv').
%             Alternatively, 'type' can specify the exact flip angle [in
%             degrees] of the pulse.  This can be useful if you want to use
%             a pulse that does not have an exact flip angle of 90 or 180
%             degrees (Thanks to Eric Plitman for help with this feature).
% off_res   = set to "1" if the centre of the RF pulse is not at 0 Hz.
%             Optional. Default=0.
%
% OUTPUTS:
% RF_struct = RF pulse waveform in FID-A rf pulse structure format.

function RF_struct=io_loadRFwaveform(filename,type,off_res);

if nargin<3
    off_res=0;
end

if isnumeric(type)
    %If 'type' was specified to be exactly 90 degrees, set the type to 'exc';
    if type==90;
        type='exc';
    end

    %If 'type' was specified to be exactly 180 degrees, set the type to 'inv'
    %(functionally there is no real difference between 'ref' and 'inv', so I
    %just chose 'Inv';
    if type==180;
        type='inv';
    end
end

%Now read in the waveform:
if isstr(filename)
    if exist(filename)
        if filename(end-3:end)=='.pta'
            disp('Siemens format .pta RF pulse file detected!! Loading waveform now.');
            rf=io_readpta(filename);
        elseif filename(end-2:end)=='.RF'
            disp('Varian/Agilent format .RF RF pulse file detected!! Loading waveform now.');
            rf=io_readRF(filename);
        elseif filename(end-3:end)=='.inv'
            disp('Bruker format .inv RF pulse file detected!! Loading waveform now.');
            rf=io_readRFBruk(filename);
        elseif filename(end-3:end)=='.rfc'
            disp('Bruker format .pta RF pulse file detected!! Loading waveform now.');
            rf=io_readRFBruk(filename);
        elseif filename(end-3:end)=='.exc'
            disp('Bruker format .exc RF pulse file detected!! Loading waveform now.');
            rf=io_readRFBruk(filename);
        elseif filename(end-3:end)=='.txt'
            disp('Basic .txt format RF pulse file detected!! Loading waveform now.');
            rf=io_readRFtxt(filename);
        else
            error('ERROR:  RF Pulse file not recognized.  Aborting!');
        end
    else
        error('ERROR:  File not found!  Aborting!');
    end
else
    if ismatrix(filename) && ndims(filename)==2 && ( 2 <= size(filename,2) <= 4 )
        disp('Input is an RF waveform already in the matlab workspace.  Loading waveform now.');
        if size(filename,2)==2
            rf=filename;
            rf(:,3)=ones(length(filename(:,1)),1);
        elseif size(filename,2)==3 || size(filename,2)==4
            rf=filename;
        end
    end
end

Tp=0.005;  %assume a 5 ms rf pulse;

%Frequency shifting the pulse to zero if off_res is 1. To be shifted back
%afterwards - PT,2021
if off_res
    %Initializations
    tmp_w1=0.03;
    tmp_min=1;

    %Looping through until a min of <-0.98 (inv/ref) or <0.2 (exc) is reached
    switch type
        case 'exc'
            min_floor = 0.2;
        case {'inv','ref'}
            min_floor = -0.98;
    end
    while tmp_min > min_floor
        [mv,sc]=bes(rf,Tp*1000,'f',tmp_w1,-5,5,1000); %running bes with bare settings
        [tmp_min,tmp_min_pos]=min(mv(3,:)); %finding the min and min position
        freq_shift=sc(tmp_min_pos)*1000;%the freq shift (in Hz) of the pulse @ 5ms
        tmp_w1=tmp_w1+0.02;
    end

    %frequency shifting the pulse to 0 Hz
    N=size(rf,1);
    dt=Tp/N;
    t=[0:dt:Tp-dt];
    phaseRamp=t*-freq_shift*360;
    rf(:,1)=rf(:,1)+phaseRamp';
end

%Find out if the pulse is phase modulated.  If it is not, then we can
%determine the time-w1 product of the pulse quite simply.  If it is phase
%modulated (adiabatic, etc) then the determination of the time-w1 product
%will need to me more interactive.
a=(round(rf(:,1))==180)|(round(rf(:,1))==0);

if sum(a)<length(rf(:,1))
    isPhsMod=true;
else
    isPhsMod=false;
end

%If there are any phase discontinuities in the phase function that are
%equal to a 360 degree jump we can remove these.  This will make it easier
%for rf_resample to do it's job later on:
jumps=diff(rf(:,1));
jumpsAbs=(abs(jumps)>355 & abs(jumps)<365);  %Assume jumps within this range are exactly = 360 degrees.
jumpIndex=find(jumpsAbs);
for n=1:length(jumpIndex)
    rf(jumpIndex(n)+1:end,1)=rf(jumpIndex(n)+1:end,1)-(360*(jumps(jumpIndex(n))/abs(jumps(jumpIndex(n)))));
end

%scale amplitude function so that maximum value is 1:
rf(:,2)=rf(:,2)./max(rf(:,2));

if ~isPhsMod
    %The pulse is not phase modulated, so we can calculate the w1max:
    %find the B1 max of the pulse in [kHz]:
    if isstr(type)
        if type=='exc'
            flipCyc=0.25; %90 degrees is 0.25 cycles;
        elseif type=='ref'
            flipCyc=0.5;  %180 degress is 0.5 cycles;
        elseif type=='inv'
            flipCyc=0.5;  %180 degrees is 0.5 cycles;
        end
    elseif isnumeric(type) %assume that a flip angle (in degrees) was given
        flipCyc=type/360;
    end

    intRF=sum(rf(:,2).*((-2*(rf(:,1)>179))+1))/length(rf(:,2));

    if intRF~=0
        w1max=flipCyc/(intRF*Tp); %w1max is in [Hz]
    else
        w1max=0;
    end
    tw1=Tp*w1max;
else
    %The pulse is phase modulated, so we will need to run some test to find
    %out the w1max;  To do this, we can plot Mz as a function of w1 and
    %find the value of w1 that results in the desired flip angle.
    [mv,sc]=bes(rf,Tp*1000,'b',0,0,5,40000);
    plot(sc,mv(3,:));
    xlabel('w1 (kHz)');
    ylabel('mz');
    w1max=input('Input desired w1max in kHz (for 5.00 ms pulse):  ');
    w1max=w1max*1000; %convert w1max to [Hz]
    tw1=Tp*w1max;
end

%now it's time to find out the time-bandwidth product:
%First make a high resolution plot the pulse profile over a wide bandwidth:
[mv,sc]=bes(rf,Tp*1000,'f',w1max/1000,-5,5,100000);
if isstr(type)
    if type=='exc'
        index=find(abs(mv(1,:)+1i*mv(2,:))>0.5);
        bw=sc(index(end))-sc(index(1))
        %plot(sc(index),mv(3,index),'.-',sc,mv(3,:));
    elseif type=='ref'
        index=find(mv(3,:)<0);
        bw=sc(index(end))-sc(index(1));
        %plot(sc(index),mv(3,index),'.-',sc,mv(3,:));
    elseif type=='inv'
        index=find(mv(3,:)<0);
        bw=sc(index(end))-sc(index(1));
        %plot(sc(index),mv(3,index),'.-',sc,mv(3,:));
    end
elseif isnumeric(type)
    mz=cos(type); %Find out the Mz value immediately following the pulse.
    thr=(1+mz)/2;  %Find out the Mz value mid-way between 1 and the mz (half-max):
    index=find(mv(3,:)<thr);  %Find the indices of the corresponding "full width"
    bw=sc(index(end))-sc(index(1));  %Now find the bandwidth at that point ("Full width at half max").
end

%Now make a very high resolution plot the pulse profile over a narrower bandwidth:
[mv,sc]=bes(rf,Tp*1000,'f',w1max/1000,-bw,bw,100000);
if isstr(type)
    if type=='exc'
        index=find(abs(mv(1,:)+1i*mv(2,:))>0.5);
        bw=sc(index(end))-sc(index(1))
        %plot(sc(index),mv(3,index),'.-',sc,mv(3,:));
    elseif type=='ref'
        index=find(mv(3,:)<0);
        bw=sc(index(end))-sc(index(1));
        %plot(sc(index),mv(3,index),'.-',sc,mv(3,:));
    elseif type=='inv'
        index=find(mv(3,:)<0);
        bw=sc(index(end))-sc(index(1));
        %plot(sc(index),mv(3,index),'.-',sc,mv(3,:));
    end
elseif isnumeric(type)
    mz=cos(type); %Find out the Mz value immediately following the pulse.
    thr=(1+mz)/2;  %Find out the Mz value mid-way between 1 and the mz (half-max):
    index=find(mv(3,:)<thr);  %Find the indices of the corresponding "full width"
    bw=sc(index(end))-sc(index(1));  %Now find the bandwidth at that point ("Full width at half max").
end

%finally, try to estimate the centre point of the RF pulse 'rfCentre'.
%This is a number between 0 and 1 that indicates where the peak of the RF
%amplitude waveform occurs.  If rfCentre = 0.5, then the pulse is symmetric.
%If rfCentre<0.5 then the peak occurs near the beginning of the pulse (i.e.
%a max phase pulse).  If rfCentre >0.5, then the peak occurs near the end
%of the pulse (i.e. min phase pulse).
maxIndex=find(rf(:,2)==max(rf(:,2)));
rfCentre=mean(maxIndex)/length(rf(:,2));

%Now store the output structure;
RF_struct.waveform=rf;
RF_struct.type=type;
RF_struct.f0=0;
if size(rf,2)>3 && any(rf(:,4))
    RF_struct.tbw='N/A - gradient modulated pulse';
    RF_struct.isGM=true;
    RF_struct.tthk=bw*Tp; %This is the time x sliceThickness product for
                         %gradient modulated pulses.  It is in units [cm.s]
else
    RF_struct.tbw=bw*Tp*1000;
    RF_struct.isGM=false;
    RF_struct.tthk='N/A - frequency selective pulse';
end
RF_struct.tw1=tw1;
RF_struct.rfCentre=rfCentre;

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
%% Steering Control Simple
% Steering Control of Autonomous Vehicles in Obstacle Avoidance Maneuvers.
%
% <<SteeringControlSimple.gif>>
%
% <html>
% <script src='https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'></script>
% </html>
%
%% Vehicle model
% *Nonlinear model*
%
% State vector
%
% \[ {\bf x} = \left[ \begin{array}{c} {\rm x}_1 \\ {\rm x}_2 \\ {\rm x}_3 \\ {\rm x}_4 \\ {\rm x}_5 \\ {\rm x}_6 \end{array} \right] = \left[ \begin{array}{c} x \\ y \\ \psi \\ v_{\rm T} \\ \alpha_{\rm T} \\ \dot{\psi} \end{array} \right] \]
%
% State equations
%
% \[ \dot{{\rm x}}_1 = {\rm x}_4 \cos \left( {\rm x}_3 + {\rm x}_5 \right) \]
%
% \[ \dot{{\rm x}}_2 = {\rm x}_4 \sin \left( {\rm x}_3 + {\rm x}_5 \right) \]
%
% \[ \dot{{\rm x}}_3 = {\rm x}_6 \]
%
% \[ \dot{{\rm x}}_4 = \frac{F_{y,{\rm F}} \sin \left( {\rm x}_5 - \delta \right) + F_{y,{\rm R}} \sin {\rm x}_5}{m_{T}} \]
%
% \[ \dot{{\rm x}}_5 = \frac{F_{y,{\rm F}} \cos \left( {\rm x}_5 - \delta \right) + F_{y,{\rm R}} \cos \alpha_{\rm T} - m_{T} {\rm x}_4 {\rm x}_6}{m_{T} {\rm x}_4} \]
%
% \[ \dot{{\rm x}}_6 = \frac{F_{y,{\rm F}} a \cos \delta - F_{y,{\rm R}} b}{I_{T}} \]
%
% Slip angles
%
% \[ \alpha_{\rm F} = \arctan \left( \frac{v_{\rm T} \sin \alpha_{\rm T} + a \dot{\psi}}{ v_{\rm T} \cos \alpha_{\rm T}} \right) - \delta \]
%
% \[ \alpha_{\rm R} = \arctan \left( \frac{v_{\rm T} \sin \alpha_{\rm T} - b \dot{\psi}}{ v_{\rm T} \cos \alpha_{\rm T}} \right) \]
%
%
% *Linear model*
%
% \[ \dot{x} = v_{\rm T} \]
%
% \[ \dot{y} = v_{{\rm T},0} \left( \psi + \alpha_{{\rm T}}\right) \]
%
% \[ \dot{\psi} = \dot{\psi} \]
%
% \[ \dot{v}_{\rm T} = 0 \]
%
% \[ \dot{\alpha}_{\rm T} = \frac{F_{y,{\rm F}} + F_{y,{\rm R}}}{m_{T} v_{{\rm T},0}} - \dot{\psi} \]
%
% \[ \ddot{\psi} = \frac{a F_{y,{\rm F}} -  b F_{y,{\rm R}}}{I_{T}} \]
%
% Neglecting equations of \(x\) and \(v_T\)
%
%
% \[ \left[ \begin{array}{c} \dot{y} \\ \dot{\psi} \\ \dot{\alpha}_T \\ \ddot{\psi} \end{array} \right] = \left[ \begin{array}{cccc} 0 & v_{T,0} & v_{T,0} & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & -\frac{K_F+K_R}{m_T v_{T,0}} & - \frac{m_T v_{T,0} + \frac{a K_F - b K_R}{v_{T,0}}}{m_T v_{T,0}} \\ 0 & 0 & - \frac{a K_F - b K_R}{I_T} & - \frac{a^2 K_F + b^2 K_R}{I_T v_{T,0}} \end{array} \right] \left[ \begin{array}{c} y \\ \psi \\ \alpha_T \\ \dot{\psi} \end{array} \right] + \left[ \begin{array}{c} 0 \\ 0 \\ \frac{K_F}{m_T v_{T,0}} \\ \frac{a K_F}{I_T}  \end{array} \right] \delta \]
%
% Slip angles
%
% \[ \alpha_{{\rm F},lin} = \alpha_{{\rm T}} + \frac{a}{v_{{\rm T},0}} \dot{\psi} - \delta \]
%
% \[ \alpha_{{\rm F},lin} = \alpha_{{\rm T}} - \frac{b}{v_{{\rm T},0}} \dot{\psi} \]
%
%% Tire model
%
% Typical characteristic curve and slip angle definition
%
% <<../../../docs/illustrations/CurvaCaracteristica.svg>>
%
% *Pacejka*
%
% \[ F_{y} = D \sin \left[ C \arctan{B \alpha - E( B \alpha -\arctan(B \alpha))} \right] \]
%
% *Linear*
%
% \[ F_ y = K \alpha \]
%

clear ; close all ; clc

slipAngle = (0:0.1:15)*pi/180;         % Slip angle [rad]

% Pacejka tire parameters
a0  = 1;
a1  = 0;
a2  = 800;
a3  = 10000;
a4  = 50;
a5  = 0;
a6  = 0;
a7  = -1;
a8  = 0;
a9  = 0;
a10 = 0;
a11 = 0;
a12 = 0;
a13 = 0;

TirePac = VehicleDynamicsLateral.TirePacejka();

Fz          = 4e+03;
camber      = 0;
TirePac.a0  = a0;
TirePac.a1  = a1;
TirePac.a2  = a2;
TirePac.a3  = a3;
TirePac.a4  = a4;
TirePac.a5  = a5;
TirePac.a6  = a6;
TirePac.a7  = a7;
TirePac.a8  = a8;
TirePac.a9  = a9;
TirePac.a10 = a10;
TirePac.a11 = a11;
TirePac.a12 = a12;
TirePac.a13 = a13;

muy0    = TirePac.a1 * Fz/1000 + TirePac.a2;
D       = muy0 * Fz/1000;
BCD     = TirePac.a3 * sin(2 * atan(Fz/1000/TirePac.a4))*(1-TirePac.a5 * abs(camber));

% Pneu linear equivalente

Ktire   = BCD * 180/pi;

TireLin     = VehicleDynamicsLateral.TireLinear();
TireLin.k   = Ktire;

% Lateral force
FyPac = TirePac.Characteristic(slipAngle, Fz, muy0/1000);
FyLin = TireLin.Characteristic(slipAngle);

% Graphics
g = VehicleDynamicsLateral.Graphics(TirePac);

%%
% Comparison of tire models

figure(1)
ax = gca;
set(ax, 'NextPlot', 'add', 'Box', 'on', 'XGrid', 'on', 'YGrid', 'on')
p = plot(slipAngle * 180/pi,-FyLin, 'Color', 'g', 'Marker', 's', 'MarkerFaceColor', 'g', 'MarkeredgeColor', 'k', 'MarkerSize', 7);
g.changeMarker(p, 10);
p = plot(slipAngle * 180/pi,-FyPac, 'Color', 'r', 'Marker', 'o', 'MarkerFaceColor', 'r', 'MarkeredgeColor', 'k', 'MarkerSize', 7);
g.changeMarker(p, 10);
xlabel('\(\alpha\) [grau]', 'Interpreter', 'Latex')
ylabel('\(F_y\) [N]', 'Interpreter', 'Latex')
l = legend('Linear', 'Pacejka');
set(l, 'Interpreter', 'Latex', 'Location', 'NorthWest')

%
%% Plant model
% Nonlinear vehicle + Pacejka tire

% Choosing vehicle
% System = VehicleDynamicsLateral.VehicleSimpleLinear();
VehiclePlant = VehicleDynamicsLateral.VehicleSimpleNonlinear();
% Defining vehicle parameters
VehiclePlant.mF0 = 700;
VehiclePlant.mR0 = 600;
VehiclePlant.IT = 10000;
VehiclePlant.lT = 3.5;
VehiclePlant.nF = 1;
VehiclePlant.nR = 1;
VehiclePlant.wT = 1.8;
VehiclePlant.muy = 1;
VehiclePlant.tire = TirePac;
VehiclePlant.deltaf = @ControlLaw;

disp(VehiclePlant)

% Choosing simulation
T       = 9;                     % Total simulation time             [s]
resol   = 500;                    % Resolution
TSPAN   = 0:T/resol:T;            % Time span                         [s]

simulator = VehicleDynamicsLateral.Simulator(VehiclePlant, TSPAN);
simulator.V0 = 16.7;

%% Controller design
%
% *Vehicle parameters*
mT  = 1300;
IT  = 10000;
a   = 1.6154;
b   = 1.8846;
vT0 = 16.7;
KF  = Ktire;
KR  = Ktire;

%%
% *Linear system*
%

A = [      0   vT0            vT0                         0                       ;...
           0    0              0                          1                       ;...
           0    0      -(KF+KR)/(mT*vT0)  -(mT*vT0+(a*KF-b*KR)/(vT0))/(mT*vT0)    ;...
           0    0      -(a*KF-b*KR)/IT    -(a^2*KF+b^2*KR)/(IT*vT0)               ];

B = [   0                  ;...
        0                  ;...
        KF/(mT*vT0)        ;...
        a*KF/IT            ];

C = [1 0 0 0];

%%
% A
disp(A)
%%
% B
disp(B)
%%
% C
disp(C)

%%
% *LQR design*
%

Q = [   0.3 0 0 0 ;...
        0   1 0 0 ;...
        0   0 1 0 ;...
        0   0 0 1 ];

%%
% Q

disp(Q)

R = 1;

%%
% R

disp(R)

Klqr = lqr(A,B,Q,R);

%%
% Klqr

disp(Klqr)

%%
% *Pole placement design*
%

polos = [-6 -6.3 -6.7 -7];

Kplace = place(A,B,polos);

%%
% Kplace

disp(Kplace)

%%
% Control law
%
% \[\delta = - {\bf K} {\bf z} + K_1 r \]
%
%% Double Lane Change Maneuver

% Simulation
simulator.Simulate();

g = VehicleDynamicsLateral.Graphics(simulator);
g.Frame();

    % Adding the double lane change track to the frame figure
    carWidth = 2;
    LaneOffset = 3.5;

    section1width = 1.1*carWidth + 0.25;
    section3width = 1.2*carWidth + 0.25;
    section5width = 1.3*carWidth + 0.25;

    section1Inf = -section1width/2;
    section1Sup = section1width/2;

    section3Inf = section1Inf+LaneOffset;
    section3Sup = section3Inf+section3width;
    section3Center = (section3Inf+section3Sup)/2;

    section5Inf = -section5width/2;
    section5Sup = section5width/2;

    % Section 1
    plot([0 15],[section1Inf section1Inf],'k')            % linha inferior
    plot([0 15],[section1Sup section1Sup],'k')            % linha superior
    plot([0 15],[0 0],'k--')                % linha central
    % Section 2
    plot([15 45],[0 section3Center],'k--')  % linha central
    % Section 3
    plot([45 70],[section3Inf section3Inf],'k')        % linha inferior
    plot([45 70],[section3Sup section3Sup],'k')        % linha superior
    plot([45 70],[section3Center section3Center],'k--')               % linha central
    % Section 4
    plot([70 95],[section3Center 0],'k--')
    % Section 5
    plot([95 130],[section5Inf section5Inf],'k')
    plot([95 130],[section5Sup section5Sup],'k')
    plot([95 130],[0 0],'k--')

g.Animation();
g.Animation();
% g.Animation('html/SteeringControlSimple');       % Uncomment to save animation gif

% Retrieving states
XT      = simulator.XT;
YT      = simulator.YT;
PSI     = simulator.PSI;
VEL     = simulator.VEL;
ALPHAT  = simulator.ALPHAT;
dPSI    = simulator.dPSI;

x = [YT PSI ALPHAT dPSI];

u = zeros(length(TSPAN),1);
output = zeros(length(TSPAN),1);

LateralDisp = 3.6;

for ii = 1:length(TSPAN)
    if XT(ii) <= 15
        r = 0;
    end
    if XT(ii) > 15 && XT(ii) <= 70
        r = LateralDisp;
    end
    if XT(ii) > 70
        r = 0;
    end

    u(ii) = - Kplace*x(ii,:)' + Kplace(1)*r;

    % Saturation
    if abs(u(ii)) < 42*pi/180
        output(ii) = u(ii);
    else
        output(ii) = sign(u(ii))*42*pi/180;
    end
end

% States
f = figure;
set(f,'PaperUnits','centimeters')
set(f,'PaperPosition',[0 0 8.9 5])
PaperPos = get(f,'PaperPosition');
set(f,'PaperSize',PaperPos(3:4))
hold on; box on; grid on
plot(TSPAN,YT,'r')
plot(TSPAN,PSI,'g')
plot(TSPAN,ALPHAT,'b')
plot(TSPAN,dPSI,'c')
xlabel('Time [s]')
ylabel('States')
l = legend('\(y\)','\(\psi\)','\(\alpha_T\)','\(\dot{\psi}\)');
set(l,'Interpreter','Latex','Location','NorthEast')

% Steering input
f = figure;
set(f,'PaperUnits','centimeters')
set(f,'PaperPosition',[0 0 8.9 3.5])
PaperPos = get(f,'PaperPosition');
set(f,'PaperSize',PaperPos(3:4))
hold on; box on; grid on
plot(TSPAN,output*180/pi,'k')
xlabel('Time [s]')
y = ylabel('\(\delta [deg]\)');
set(y,'Interpreter','Latex')

%% See Also
%
% <../../../index.html Home>
%

author: andresmendes

repo: Vehicle-Dynamics-Lateral


text:
function[b]=minmin(x)
%MINMIN  MINMIN(X)=MIN(X(ISFINITE(X)))
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2004--2012 J.M. Lilly --- type 'help jlab_license' for details
b=min(x(isfinite(x(:))));

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
classdef PoseGraph < handle
    %POSEGRAPH A class for doing pose graph optimization

    properties (SetAccess = private)
        node  % Pose nodes in graph
        edge  % Edge in graph
        H     % Information matrix
        b     % Information vector
    end  % properties set private

    properties (Dependent = true)
        n_node  % Number of nodes in graph
        n_edge  % Number of edges in graph
        pose    % Poses of all nodes
    end  % properties dependent

    methods

        function obj = PoseGraph()
            % Constructor of PoseGraph
            obj.node = PoseNode.empty;
            obj.edge = PoseEdge.empty;
        end

        function readGraph(obj, vfile, efile)
            % Reads graph from vertex and edge file

            % Try opening vertex file
            vfid = fopen(vfile);
            if (vfid < 0)
                fprintf('Fail to open %s\n.', vfile);
                return
            end
            vertices = fscanf(vfid, 'VERTEX2 %d %f %f %f\n', [4 Inf]);
            for i_node = 1:size(vertices,2)
                vi = vertices(:,i_node);
                id = vi(1) + 1;
                pose = vi(2:4);
                obj.node(i_node) = PoseNode(id, pose);
            end
            fclose(vfid);
            fprintf('Vertices loaded from: %s\n', vfile);
            % Try opening edge file
            efid = fopen(efile);
            if (efid < 0)
                fprintf('Fail to open %s\n.', vfile);
                return
            end
            edges = fscanf(efid,...
                'EDGE2 %d %d %f %f %f %f %f %f %f %f %f \n',[11 Inf]);
            for i_edge = 1:size(edges,2)
                ei = edges(:,i_edge);
                id_from   = ei(1) + 1;
                id_to     = ei(2) + 1;
                mean      = ei(3:5);
                infm      = zeros(3,3);
                infm(1,1) = ei(6,:);
                infm(2,1) = ei(7,:);
                infm(1,2) = ei(7,:);
                infm(2,2) = ei(8,:);
                infm(3,3) = ei(9,:);
                infm(1,3) = ei(10,:);
                infm(3,1) = ei(10,:);
                infm(3,2) = ei(11,:);
                infm(2,3) = ei(11,:);
                obj.edge(i_edge) = PoseEdge(id_from, id_to, mean, infm);
            end
            fclose(efid);
            fprintf('Edges loaded from: %s\n', vfile);
        end

        function plot(obj)
            % Plots pose graph
            obj.node.plot();
        end

        function optimize(obj, n_iter, vis)
            % Pose graph optimization
            if nargin < 3, vis = false; end
            if nargin < 2, n_iter = 1; end

            for i_iter = 1:n_iter
                fprintf('Pose Graph Optimization, Iteration %d.\n', i_iter);
                obj.iterate();
                fprintf('Iteration %d done.\n', i_iter);
                if vis
                    clf;
                    obj.plot();
                    title(sprintf('Iteration %d', i_iter));
                    drawnow;
                end
            end
        end

        function iterate(obj)
            % One iteration of pose graph optimization
            fprintf('Allocating Workspace.\n');
            % Create new H and b matrices each time
            obj.H = zeros(obj.n_node*3);   % 3n x 3n square matrix
            obj.b = zeros(obj.n_node*3,1); % 3n x 1  column vector

            fprintf('Linearizing.\n');
            obj.linearize();

            fprintf('Solving.\n');
            obj.solve();
        end

        function linearize(obj)
            % Linearize error functions and formulate a linear system
            for i_edge = 1:obj.n_edge
                ei = obj.edge(i_edge);
                % Get edge information
                i_node = ei.id_from;
                j_node = ei.id_to;
                T_z = v2t(ei.mean);
                omega = ei.infm;

                % Get node information
                v_i = obj.node(i_node).pose;
                v_j = obj.node(j_node).pose;
                i_ind = id2ind(i_node);
                j_ind = id2ind(j_node);

                T_i = v2t(v_i);
                T_j = v2t(v_j);
                R_i = T_i(1:2,1:2);
                R_z = T_z(1:2,1:2);

                si = sin(v_i(3));
                ci = cos(v_i(3));
                dR_i = [-si ci; -ci -si]';
                dt_ij = v_j(1:2) - v_i(1:2);

                % Caluclate jacobians
                A = [-R_z'*R_i' R_z'*dR_i'*dt_ij; 0 0 -1];
                B = [R_z'*R_i' [0;0]; 0 0 1];

                % Calculate error vector
                e = t2v(inv(T_z) * inv(T_i) * T_j);

                % Formulate blocks
                H_ii =  A' * omega * A;
                H_ij =  A' * omega * B;
                H_jj =  B' * omega * B;
                b_i  = -A' * omega * e;
                b_j  = -B' * omega * e;

                % Update H and b matrix
                obj.H(i_ind,i_ind) = obj.H(i_ind,i_ind) + H_ii;
                obj.H(i_ind,j_ind) = obj.H(i_ind,j_ind) + H_ij;
                obj.H(j_ind,i_ind) = obj.H(j_ind,i_ind) + H_ij';
                obj.H(j_ind,j_ind) = obj.H(j_ind,j_ind) + H_jj;
                obj.b(i_ind) = obj.b(i_ind) + b_i;
                obj.b(j_ind) = obj.b(j_ind) + b_j;
            end
        end

        function solve(obj)
            % Solves the linear system and update all pose node
            fprintf('Pose: %d, Edge: %d\n', obj.n_node, obj.n_edge);
            % The system (H b) is obtained only from relative constraints.
            % H is not full rank.
            % We solve this by anchoring the position of the 1st vertex
            % This can be expressed by adding teh equation
            % dx(1:3,1) = 0
            % which is equivalent to the following
            obj.H(1:3,1:3) = obj.H(1:3,1:3) + eye(3);
            H_sparse = sparse(obj.H);
            dx = H_sparse \ obj.b;
            dpose = reshape(dx, 3, obj.n_node);

            % Update node with solution
            for i_node = 1:obj.n_node
                obj.node(i_node).pose = obj.node(i_node).pose ...
                    + dpose(:,i_node);
            end
        end

        function n_node = get.n_node(obj)
            n_node = numel(obj.node);
        end

        function n_edge = get.n_edge(obj)
            n_edge = numel(obj.edge);
        end

        function pose = get.pose(obj)
            pose = [obj.node.pose];
        end

    end  % methods public

end  % classdef

function ind = id2ind(id)
%ID2IND converts id to indices in H and b
ind = (3*(id-1)+1):(3*id);
end

repo: graphslam


text:
%
%  Map a sequence of pen points into a Motor Program
%
%  Input
%   S: [ns x 1 nested cell array] of trajectories in motor space [n x 2]
%   I: [binary image]
%   fclassify: [function handle] fclassify([n x 2] trajectory, [scalar]
%                                   scale)
%               this function classifies trajectories are belong to a
%               specific library component
%
%  Output
%   M: MotorProgram object
%
function M = sequence_to_MP(S,I,fclassify,lib)

    assert(iscell(S));
    assert(iscell(S{1}));

    % find the new scale of primitives
    assert(exist('lib','var')>0);
    newscale = lib.newscale;
    ncpt = lib.ncpt;

%     % find the new scale of primitives
%     ps_clustering = defaultps_clustering;
%     newscale = ps_clustering.newscale_ss;
%     ncpt = ps_clustering.ncpt;
%     if exist('lib','var')
%        newscale = lib.newscale;
%        ncpt = lib.ncpt;
%     end

    % get the normalized trajectories, the stroke centers, scales, etc.
    verbose = false;
    [S_norm,~,S_scales] = normalize_dataset(S,newscale,verbose);

    % fit splines to the normalized trajectories, and assign them
    % indices in the library
    PM = defaultps;
    ps_bottomup = defaultps_bottomup;
    ns = length(S);
    S_splines = cell(ns,1);
    for sid=1:ns
        nsub = length(S{sid});
        S_splines{sid} = cell(nsub,1);
        for b=1:nsub
            S_splines{sid}{b} = fclassify(S_norm{sid}{b},S_scales{sid}{b});
        end
    end

    % make motor program structure
    M = MotorProgram(ns);
    M.I = I; % set the image
    M.parameters = PM;
    for sid=1:ns % for each stroke
       nsub = length(S{sid});
       M.S{sid}.ids = []; %zeros(nsub,1);
       M.S{sid}.invscales_token = zeros(nsub,1);
       M.S{sid}.shapes_token = zeros(ncpt,2,nsub);

       for b=1:nsub % for each sub-stroke
          M.S{sid}.ids(b,:) = S_splines{sid}{b}.indx;
          M.S{sid}.invscales_token(b) = 1 ./ S_scales{sid}{b}(1);
          M.S{sid}.shapes_token(:,:,b) = S_splines{sid}{b}.bspline;
       end
       M.S{sid}.invscales_type = M.S{sid}.invscales_token;

       M.S{sid}.pos_token = S{sid}{1}(1,:); % start position in real trajectory

    end

    % initialize the noise parameters
    M.epsilon = ps_bottomup.init_epsilon;
    M.blur_sigma = ps_bottomup.init_blur_sigma;

end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
function [route] = astar3(exbigraph,exbiloc,startnode,endnode,ts)
%% Written by Muhammet Balcilar, France
% all rights reverved

graph=exbigraph;
Loc=exbiloc;
n=size(graph,1);

curr=startnode;

Cost=[curr 0];
Hcost= sum(sum(abs(Loc(curr,:)-Loc(endnode,:))));

Route{1}=[curr];

Tabu=[curr];
cRoute=[curr];
cr=1;

i=1;
% if the endnode is found or all possibility is finished
while Hcost(cr)>0 %      sum(curr~=endnode)>0
    i=i+1;
    iter=i

    % find all possibility from current node
    %PosMove2=findpossmove(graph,curr,ts);
    PosMove=findpossmoveRec(graph,curr,ts);

    % find all possibility's cost from current to concerned node
    for j=1:size(PosMove,1)
        % manhattan distance for current and possible locs
        cost(j,1)=sum(sum( abs(Loc(curr,:)-Loc(PosMove(j,:),:))))+findindex(Cost,curr);
        % add time cost every step cost is just 1 move
        cost(j,1)=cost(j,1)+1;
        [val ind]=findindex(Cost,PosMove(j,:));
        if cost(j,1)<val
            Cost(ind,:)=[PosMove(j,:) cost(j,1)];
            Route{ind}=[cRoute; PosMove(j,:)];
            Hcost(ind,1)=sum(sum(abs(Loc(PosMove(j,:),:)-Loc(endnode,:))));
        end

    end

    % find minimum of normal cost plus heuristic cost of possible node
    tmp=Cost;
    tmp(:,end)=tmp(:,end)+Hcost;
    [mn cr]=finminexceptTabu(tmp,Tabu);

    if isinf(mn)
        break
    end
    % set minimum total cost's node as current node
    curr=Cost(cr,1:end-1);
    cRoute=Route{cr};
    % add current node to tabu list to prevent loops
    Tabu=[Tabu; curr];
end
% if current node is equal endnode means route is found
if Hcost(cr)==0
    route=cRoute;
else % i not means there is no way
    route=[];
end

author: balcilar

repo: Multi-Robot-Path-Planning-on-Graphs


text:
%kDeconvAxial 'Deconvolution using ML/EM algorithm '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros DeconvAxial.pane file
%
% Parameters:
% Integer: n 'Iterations', default: 60: 'Number of iterations until algorythm stops'
% Integer: rW 'Border (Pixels)', default: 4: 'a gaussian Window will be applied'
% Toggle: norichardson 'No Richardson', default: 0: 'do not use the richardson iteration, but the stadard ML-correction method instead'
% Toggle: noPSFnorm 'Do not normalize PSFs', default: 0: 'The integral intensity of the PSFs will then not be normalized'
% Double: mf 'Maximal Rel. Frequency X/Y', default: 0.8: 'Higher frequencies are suppressed in X/Y direction'
% Double: mzf 'Z', default: 0.4: 'Higher frequencies are suppressed in Z-direction'
% InputFile: i 'Measured Image', required: 'Measured confocal image'
% InputFile: p 'PSF', required: 'PSF needed for ML algorithm (same size as measured data)'
% InputFile: s 'Start File', optional: 'file to start iteration with'
% InputFile: Over 'Overrelaxation Table', optional: 'Table for Overrelaxation Factors'
% InputFile: patterns 'Pattern File', optional: 'File containing the different patterns (along elements)'
% Toggle: patsep 'Separable Patterns', default: 0: 'If selected, it is assumed, that the pattern can be described as a sum of separable components. In this case, the pattern-file has to consist of N planes (stacked along depth) for each element. The PSF elements are then assumed to contain the Z-modified PSFs.'
% Toggle: HolmesLiu 'Holmes & Liu', default: 0: 'If selected, the overrelaxation is choosen according to the Holmes&Liu paper. This specific implementation is however rather slow.'
% OutputFile: o 'Recon. Image', required: 'Resulting output data object'
% OutputFile: c 'Correction Image', required: 'Correction Image'
% OutputFile: f 'Fwd Projected File', optional: 'forward proj. guess'
% OutputFile: hf 'HF Content Output', optional: 'energy of highest 10% frequences'
% Toggle: background 'Estimate Background', default: 0: 'For patterned excitation data, a background (out-of-focus, not deconvolved) can be estimated.'
% OutputFile: bgFile 'Background', optional: 'Only for patterned excitation data, Estimate Background active'
%
% Example: [o, c, f, hf, bgFile] = kDeconvAxial({i, p, s, Over, patterns}, {'n',60;'rW',4;'norichardson',0;'noPSFnorm',0;'mf',0.8;'mzf',0.4;'i','';'p','';'s','';'Over','';'patterns','';'patsep',0;'HolmesLiu',0;'o','';'c','';'f','';'hf','';'background',0;'bgFile',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% DeconvAxial - Deconvolutes using many Psfs
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1996-2003, Rainer Heintzmann,  All rights reserved.
%

function varargout = kDeconvAxial(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kDeconvAxial(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'n', 60;'rW', 4;'norichardson', 0;'noPSFnorm', 0;'mf', 0.8;'mzf', 0.4;'i', '__input';'p', '__input';'s', '__input';'Over', '__input';'patterns', '__input';'patsep', 0;'HolmesLiu', 0;'o', '__output';'c', '__output';'f', '__output';'hf', '__output';'background', 0;'bgFile', '__output'};
maxval={1000,2,0,0,2,2,0,0,1,1,1,0,0,0,0,1,1,0,1};
minval={1,2,0,0,0,0,0,0,1,1,1,0,0,0,0,1,1,0,1};
istoggle=[0,1,1,1,1,1,0,0,1,1,1,1,1,0,0,1,1,1,1];
was_set=istoggle * 0;
paramtype={'Integer','Integer','Toggle','Toggle','Double','Double','InputFile','InputFile','InputFile','InputFile','InputFile','Toggle','Toggle','OutputFile','OutputFile','OutputFile','OutputFile','Toggle','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=2; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'axialrec"  -k'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [rfAmp,rfPhase,rfFreq,rfCoil,rfTime]=rfBIREF(p)
%create a BIREF adiabatic excitation rf pulse starting from tStart and ending at tEnd
%tStart rf start time
%tEnd rf end time
%dt rf sample time
%rfPhase rf phase
%rfFreq rf off-res freq

tStart=p.tStart;
tEnd=p.tEnd;
dt=p.dt;
MaxB1=p.MaxB1; % Maximum B1
MaxFreq=p.MaxFreq; % Frequency modulation amplitude
BIREFFlag=p.BIREFFlag;
rfCoil=p.CoilID;
Duplicates=max(1,p.Duplicates);
DupSpacing=max(0,p.DupSpacing);

tEnd=tEnd-tStart;
tStart=0; % time scale shift to 0

rfTime=linspace(tStart,tEnd,ceil((tEnd-tStart)/dt)+1);

switch BIREFFlag
    case 'BIREF-1'
        Zeta=pi/(tEnd-tStart);
        rfTime1=rfTime(rfTime<(tEnd-tStart)/2);
        rfAmp1=MaxB1*sin(Zeta*rfTime1); % rf frequency modulation
        rfTime2=rfTime(rfTime>=(tEnd-tStart)/2);
        rfAmp2=-MaxB1*sin(Zeta*rfTime2); % rf frequency modulation

        rfAmp=[rfAmp1 rfAmp2];
        rfFreq=MaxFreq*abs(cos(Zeta*rfTime));
        rfPhase=0*ones(size(rfTime)); % rf Phase
    case 'BIREF-2a'
        Zeta=pi/(tEnd-tStart);
        rfTime1=rfTime(rfTime<(tEnd-tStart)/2);
        rfFreq1=MaxFreq*sin(Zeta*rfTime1);
        rfTime2=rfTime(rfTime>=(tEnd-tStart)/2);
        rfFreq2=-MaxFreq*sin(Zeta*rfTime2);

        rfAmp=MaxB1*abs(cos(Zeta*rfTime));
        rfFreq=[rfFreq1 rfFreq2];
        rfPhase=0*ones(size(rfTime)); % rf Phase
    case 'BIREF-2b'
        Zeta=2*pi/(tEnd-tStart);
        rfTime1=rfTime(rfTime<(tEnd-tStart)/4);
        rfAmp1=MaxB1*abs(cos(Zeta*rfTime1));
        rfFreq1=MaxFreq*sin(Zeta*rfTime1);

        rfTime2=rfTime(rfTime>=(tEnd-tStart)/4 & rfTime<(tEnd-tStart)/2);
        rfAmp2=MaxB1*abs(cos(Zeta*rfTime2));
        rfFreq2=-MaxFreq*sin(Zeta*rfTime2);

        rfTime3=rfTime(rfTime>=(tEnd-tStart)/2 & rfTime<3*(tEnd-tStart)/4);
        rfAmp3=-MaxB1*cos(Zeta*rfTime3);
        rfFreq3=-MaxFreq*sin(Zeta*rfTime3);

        rfTime4=rfTime(rfTime>=3*(tEnd-tStart)/4);
        rfAmp4=-MaxB1*cos(Zeta*rfTime4);
        rfFreq4=MaxFreq*sin(Zeta*rfTime4);

        rfAmp=[rfAmp1 rfAmp2 rfAmp3 rfAmp4];
        rfFreq=[rfFreq1 rfFreq2 rfFreq3 rfFreq4];
        rfPhase=0*ones(size(rfTime)); % rf Phase
end

rfTime=rfTime+p.tStart; % time scale shift back
rfCoil=(rfCoil)*ones(size(rfTime));
rfAmp(1)=0;
rfAmp(end)=0;
rfFreq(1)=0;
rfFreq(end)=0;
rfPhase(1)=0;
rfPhase(end)=0;

% Create Duplicates
if Duplicates~=1 & DupSpacing ~=0
    rfAmp=repmat(rfAmp,[1 Duplicates]);
    rfFreq=repmat(rfFreq,[1 Duplicates]);
    rfPhase=repmat(rfPhase,[1 Duplicates]);
    rfCoil=repmat(rfCoil,[1 Duplicates]);
    TimeOffset = repmat(0:DupSpacing:(Duplicates-1)*DupSpacing,[length(rfTime) 1]);
    rfTime=repmat(rfTime,[1 Duplicates]) + (TimeOffset(:))';
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [ y, m, d, f, ierror ] = ymdf_check_common ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_CHECK_COMMON checks a Common YMDF date.
%
%  Discussion:
%
%    Certain simple errors in dates will be corrected, such as
%      "31 September 1996"
%    which will become
%      "1 October 1996".
%
%    The routine also knows that in the Common calendar, the dates
%    5 October 1582 through 14 October 1582 are illegal.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, integer M, integer D, real  F, the
%    YMDF date.
%
%    Output, integer Y, integer M, integer D, real F, the
%    YMDF date, which may be corrected if necessary and possible.
%
%    Output, integer IERROR, is 0 if the date is legal.
%
  ierror = 0;

  [ y, m, d, ierror ] = ymd_check_common ( y, m, d );

  if ( ierror ~= 0 )
    return
  end

  [ y, m, d, f ] = frac_borrow_common ( y, m, d, f );

  [ y, m, d, f ] = frac_carry_common ( y, m, d, f );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function dat = read_egis_data(filename, hdr, begtrial, endtrial, chanindx)

% READ_EGIS_DATA reads the data from an EGI EGIS format file
%
% Use as
%   dat = read_egis_data(filename, hdr, begtrial, endtrial, chanindx);
% where
%   filename       name of the input file
%   hdr            header structure, see FT_READ_HEADER
%   begtrial       first trial to read, mutually exclusive with begsample+endsample
%   endtrial       last trial to read,  mutually exclusive with begsample+endsample
%   chanindx       list with channel indices to read
%
% This function returns a 3-D matrix of size Nchans*Nsamples*Ntrials.
% Note that EGIS session files are defined as always being epoched.
% For session files the trials are organized with the members of each cell grouped
% together.  For average files the "trials" (subjects) are organized with the cells
% also grouped together (e.g., "cell1sub1, cell1sub2, ...).
%_______________________________________________________________________
%
%
% Modified from EGI's EGI Toolbox with permission 2007-06-28 Joseph Dien

% 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$

fh=fopen_or_error([filename],'r');
fclose(fh);

[fhdr,chdr,ename,cnames,fcom,ftext] = read_egis_header(filename);
fh=fopen_or_error([filename],'r');
fhdr(1)=fread(fh,1,'int32'); %BytOrd
[str,maxsize,cEndian]=computer;
if fhdr(1)==16909060
    if cEndian == 'B'
        endian = 'ieee-be';
    elseif cEndian == 'L'
        endian = 'ieee-le';
    end
elseif fhdr(1)==67305985
    if cEndian == 'B'
        endian = 'ieee-le';
    elseif cEndian == 'L'
        endian = 'ieee-be';
    end
else
    ft_error('This is not an EGIS average file.');
end

if fhdr(2) == -1
    fileType = 'ave';
elseif fhdr(2) == 3
    fileType = 'ses';
else
    ft_error('This is not an EGIS file.');
end

dat=zeros(hdr.nChans,hdr.nSamples,endtrial-begtrial+1);

%read to end of file header
status=fseek(fh,(130+(2*fhdr(18))+(4*fhdr(19))),'bof');
status=fseek(fh,2,'cof');
for loop=1:fhdr(18)
  temp=fread(fh,80,'uchar');
  theName=strtok(temp);
  theName=strtok(theName,char(0));
  cnames{loop}=deblank(char(theName))';
  status=fseek(fh,fhdr(24+(loop-1))-80,'cof');
end
fcom=fread(fh,fhdr(20),'uchar');
ftext=fread(fh,fhdr(21),'uchar');
fpad=fread(fh,fhdr(22),'uchar');
status=fseek(fh,-2,'cof');

%read to start of desired data
fseek(fh, ((begtrial-1)*hdr.nChans*hdr.nSamples*2), 'cof');

for segment=1:(endtrial-begtrial+1)
    dat(:,:,segment) = fread(fh, [hdr.nChans, hdr.nSamples],'int16',endian);
end
dat=dat(chanindx, :,:);

if fileType == 'ave'
    dat=dat/fhdr(12); %convert to microvolts
elseif fileType == 'ses'
    dat=dat/5; %convert to microvolts (EGIS sess files created by NetStation use a 5 bins per microvolt scaling)
end

fclose(fh);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function mpc = case15da
%CASE15DA  Power flow data for 15 bus distribution system from Das, et al
%   Please see CASEFORMAT for details on the case file format.
%
%   Data from ...
%       Das D, Kothari DP, Kalam A (1995) Simple and efficient method for load
%       flow solution of radial distribution networks. Int J Electr Power
%       Energy Syst 17:335-346. doi: 10.1016/0142-0615(95)00050-0
%       URL: https://doi.org/10.1016/0142-0615(95)00050-0

%% MATPOWER Case Format : Version 2
mpc.version = '2';

%%-----  Power Flow Data  -----%%
%% system MVA base
mpc.baseMVA = 1;

%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
mpc.bus = [ %% (Pd and Qd are specified in kW & kVAr here, converted to MW & MVAr below)
	1	3	0	0	0	0	1	1	0	11	1	1	1;
	2	1	44.1	44.991	0	0	1	1	0	11	1	1.1	0.9;
	3	1	70	71.4143	0	0	1	1	0	11	1	1.1	0.9;
	4	1	140	142.8286	0	0	1	1	0	11	1	1.1	0.9;
	5	1	44.1	44.991	0	0	1	1	0	11	1	1.1	0.9;
	6	1	140	142.8286	0	0	1	1	0	11	1	1.1	0.9;
	7	1	140	142.8286	0	0	1	1	0	11	1	1.1	0.9;
	8	1	70	71.4143	0	0	1	1	0	11	1	1.1	0.9;
	9	1	70	71.4143	0	0	1	1	0	11	1	1.1	0.9;
	10	1	44.1	44.991	0	0	1	1	0	11	1	1.1	0.9;
	11	1	140	142.8286	0	0	1	1	0	11	1	1.1	0.9;
	12	1	70	71.4143	0	0	1	1	0	11	1	1.1	0.9;
	13	1	44.1	44.991	0	0	1	1	0	11	1	1.1	0.9;
	14	1	70	71.4143	0	0	1	1	0	11	1	1.1	0.9;
	15	1	140	142.8286	0	0	1	1	0	11	1	1.1	0.9;
];

%% generator data
%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin	Pc1	Pc2	Qc1min	Qc1max	Qc2min	Qc2max	ramp_agc	ramp_10	ramp_30	ramp_q	apf
mpc.gen = [
	1	0	0	10	-10	1	100	1	10	0	0	0	0	0	0	0	0	0	0	0	0;
];

%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
mpc.branch = [  %% (r and x specified in ohms here, converted to p.u. below)
	1	2	1.35309	1.32349	0	0	0	0	0	0	1	-360	360;
	2	3	1.17024	1.14464	0	0	0	0	0	0	1	-360	360;
	3	4	0.84111	0.82271	0	0	0	0	0	0	1	-360	360;
	4	5	1.52348	1.0276	0	0	0	0	0	0	1	-360	360;
	2	9	2.01317	1.3579	0	0	0	0	0	0	1	-360	360;
	9	10	1.68671	1.1377	0	0	0	0	0	0	1	-360	360;
	2	6	2.55727	1.7249	0	0	0	0	0	0	1	-360	360;
	6	7	1.0882	0.734	0	0	0	0	0	0	1	-360	360;
	6	8	1.25143	0.8441	0	0	0	0	0	0	1	-360	360;
	3	11	1.79553	1.2111	0	0	0	0	0	0	1	-360	360;
	11	12	2.44845	1.6515	0	0	0	0	0	0	1	-360	360;
	12	13	2.01317	1.3579	0	0	0	0	0	0	1	-360	360;
	4	14	2.23081	1.5047	0	0	0	0	0	0	1	-360	360;
	4	15	1.19702	0.8074	0	0	0	0	0	0	1	-360	360;
];

%%-----  OPF Data  -----%%
%% generator cost data
%	1	startup	shutdown	n	x1	y1	...	xn	yn
%	2	startup	shutdown	n	c(n-1)	...	c0
mpc.gencost = [
	2	0	0	3	0	20	0;
];

%% convert branch impedances from Ohms to p.u.
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
Vbase = mpc.bus(1, BASE_KV) * 1e3;      %% in Volts
Sbase = mpc.baseMVA * 1e6;              %% in VA
mpc.branch(:, [BR_R BR_X]) = mpc.branch(:, [BR_R BR_X]) / (Vbase^2 / Sbase);

%% convert loads from kW to MW
mpc.bus(:, [PD, QD]) = mpc.bus(:, [PD, QD]) / 1e3;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


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

test_suite = functiontests(localfunctions);

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

p1 = [10 10];
p2 = [20 10];
p3 = [20 20];
p4 = [10 20];
square = [p1;p2;p3;p4];
testCase.assertEqual(polygonPoint(square, 0), p1, 'AbsTol', .01);
testCase.assertEqual(polygonPoint(square, 1), p2, 'AbsTol', .01);
testCase.assertEqual(polygonPoint(square, 2), p3, 'AbsTol', .01);
testCase.assertEqual(polygonPoint(square, 3), p4, 'AbsTol', .01);

testCase.assertEqual(polygonPoint(square, 4), p1);
testCase.assertEqual(polygonPoint(square, 3.5), centroid([p1;p4]), 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function aistats2012_slides_gp

randn('state', 1);

plot_gp(true, 0.2, 0.4, 100, false, ...
        '/home/jluttine/papers/aistats2012/slides/fig_short', ...
        8, 8);
plot_gp(true, 1.5, 0.4, 100, false, ...
        '/home/jluttine/papers/aistats2012/slides/fig_long', ...
        8, 8);
plot_gp(false, 1, 0.2, 6, 5, ...
        '/home/jluttine/papers/aistats2012/slides/fig_gp', ...
        6, 4);

function plot_gp(noisy, lengthscale, noise, N_all, samples, filename, width, height);
xmin = 0;
xmax = 10;

Nh = 1000;
xh = linspace(xmin, xmax, Nh);

x_all = linspace(xmin+0.5, xmax-0.5, N_all);
%x_all = unifrnd(xmin, xmax, 1, N_all);

% Generate data
covfunc_all = gp_cov_sum(gp_cov_se(sq_dist(x_all), ...
                                 'lengthscale', lengthscale), ...
                       gp_cov_scale(gp_cov_delta(N_all), ...
                                    'scale', noise));

K_y = covfunc_all([]);
L_y = chol(K_y + 1e-6*eye(N_all), 'lower');
y_all = L_y * randn(N_all,1);

for N = [N_all]
  x = x_all(:,1:N);
  y = y_all(1:N);

  covfunc_y = gp_cov_sum(gp_cov_se(sq_dist(x), ...
                                   'lengthscale', lengthscale), ...
                         gp_cov_scale(gp_cov_delta(N), ...
                                      'scale', noise));
  covfunc_yf = gp_cov_se(sq_dist(x,xh), ...
                         'lengthscale', lengthscale);
  if noisy
    covfunc_f = gp_cov_sum(gp_cov_se(zeros(Nh,1), ...
                                     'lengthscale', lengthscale), ...
                           gp_cov_scale(gp_cov_wrap(ones(Nh,1)), ...
                                        'scale', 0.3));
  else
    covfunc_f = gp_cov_se(zeros(Nh,1), ...
                          'lengthscale', lengthscale);
  end

  K_y = covfunc_y([]);
  K_yf = covfunc_yf([]);
  k_f = covfunc_f([]);

  [fh_mean, fh_var] = gp_predict(y, K_y, K_yf, k_f);
  figure
  errorplot(xh, fh_mean, sqrt(fh_var))
  hold on
  plot(x, y, 'r+')
  set(gca, 'xtick', [], 'ytick', []);
  set(gca, 'box', 'off');
  set_figure_size(width,height);
  print('-dpdf', [filename, '.pdf'])
  yl = get(gca, 'ylim');

  if samples
    covfunc_f = gp_cov_se(sq_dist(xh), ...
                          'lengthscale', lengthscale);
    K_f = covfunc_f([]);
    K_post = K_f - K_yf'*(K_y\K_yf) + 1e-6*eye(Nh);
    L_post = chol(K_post, 'lower');
    f = bsxfun(@plus, ...
               K_yf'*(K_y\y), ...
               L_post*randn(Nh,samples));
    figure
    plot(xh, f);
    hold on
    plot(x, y, 'k+');
    set(gca, 'xtick', [], 'ytick', []);
    set(gca, 'box', 'off');
    set(gca, 'ylim', yl);
    set_figure_size(width,height);
    print('-dpdf', [filename, '_samples.pdf'])
  end

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
% test101 Test querying strings from a table.

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 2006-04-13 Created (Johan Meijdam, VORtech)

function test101()
   OtStart();
   hquery = OtQueryNew();
   OtQuerySetSQL(hquery, ['SELECT name AS Hnames FROM "' testdir() 'name.DB" AS nametable WHERE name LIKE "H%"']);
   if OtQueryOpen(hquery) ~= 1
      error('OtQueryOpen failed')
   end
   disp(OtQueryFields(hquery));
   OtQueryFirst(hquery);
   while ~OtQueryEof(hquery)
      disp(OtQueryGetAll(hquery));
      OtQueryNext(hquery);
   end
   OtQueryClose(hquery);
   OtQueryFree(hquery);
   OtStop();
end


text:
function [y, pst] = tapas_ceode_gen_erp(P, M, U)
% [y, pst] = tapas_ceode_gen_erp(P, M, U)
%
% Compute predicted data based on a specified DCM, with variable integrator
% specification.
%
% Adapted from spm_gen_erp.m (original function lincense below).
%
%
% INPUT
%   P           struct          parameter structure
%   M           struct          model specification
%   U           struct          design specification
%
% OUTPUT
%   y           mat             Integrated activity in sensor space
%
% -------------------------------------------------------------------------
%
% Author: Dario Schöbi
% Created: 2020-08-10
% Copyright (C) 2020 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS ceode Toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------
%
% ORIGINAL FUNCTION LICENSE
%
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
%
% Karl Friston
% $Id: spm_gen_erp.m 5758 2013-11-20 21:04:01Z karl $
%--------------------------------------------------------------------------

% default inputs - one trial (no between-trial effects)
%--------------------------------------------------------------------------
if nargin < 3, U.X = sparse(1,0); end

% peristimulus time
%--------------------------------------------------------------------------
if nargout > 1
    pst = (1:M.ns)*U.dt - M.ons/1000;
end

% within-trial (exogenous) inputs
%==========================================================================
if ~isfield(U,'u')

    % peri-stimulus time inputs
    %----------------------------------------------------------------------
    U.u = feval(M.fu,(1:M.ns)*U.dt,P,M);

end

if isfield(M,'u')

    % remove M.u to preclude endogenous input
    %----------------------------------------------------------------------
    M = rmfield(M,'u');

end

% between-trial (experimental) inputs
%==========================================================================
if isfield(U,'X')
    X = U.X;
else
    X = sparse(1,0);
end

if ~size(X,1)
    X = sparse(1,0);
end

% cycle over trials
%==========================================================================
y      = cell(size(X,1),1);
for  c = 1:size(X,1)

    % condition-specific parameters
    %----------------------------------------------------------------------
    Q = spm_gen_Q(P, X(c, :));

    % solve for steady-state - for each condition
    %----------------------------------------------------------------------
    M.x  = spm_dcm_neural_x(Q,M);

    % integrate DCM - for this condition
    %----------------------------------------------------------------------
    y{c} = feval(M.int, Q, M, U);

end


text:
function [ xstar, seed ] = rk1_tv_step ( x, t, h, q, fv, gv, seed )

%*****************************************************************************80
%
%% RK1_TV_STEP takes one step of a stochastic Runge Kutta scheme.
%
%  Discussion:
%
%    The Runge-Kutta scheme is first-order, and suitable for time-varying
%    systems.
%
%    d/dx X(t,xsi) = F ( X(t,xsi), t ) + G ( X(t,xsi), t ) * w(t,xsi)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Jeremy Kasdin,
%    Runge-Kutta algorithm for the numerical integration of
%    stochastic differential equations,
%    Journal of Guidance, Control, and Dynamics,
%    Volume 18, Number 1, January-February 1995, pages 114-120.
%
%    Jeremy Kasdin,
%    Discrete Simulation of Colored Noise and Stochastic Processes
%    and 1/f^a Power Law Noise Generation,
%    Proceedings of the IEEE,
%    Volume 83, Number 5, 1995, pages 802-827.
%
%  Parameters:
%
%    Input, real X, the value at the current time.
%
%    Input, real T, the current time.
%
%    Input, real H, the time step.
%
%    Input, real Q, the spectral density of the input white noise.
%
%    Input, external real FV, the name of the deterministic
%    right hand side function.
%
%    Input, external real GV, the name of the stochastic
%    right hand side function.
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Output, real XSTAR, the value at time T+H.
%
  a21 = 1.0;

  q1 = 1.0;

  t1 = t;
  x1 = x;
  [ n1, seed ] = r8_normal_01 ( seed );
  w1 = n1 * sqrt ( q1 * q / h );
  k1 = h * fv ( t1, x1 ) + h * gv ( t1, x1 ) * w1;

  xstar = x1 + a21 * k1;

  return
end


text:
%DEMGMM3 Demonstrate density modelling with a Gaussian mixture model.
%
%	Description
%	 The problem consists of modelling data generated by a mixture of
%	three Gaussians in 2 dimensions with a mixture model using diagonal
%	covariance matrices.  The priors are 0.3, 0.5 and 0.2; the centres
%	are (2, 3.5), (0, 0) and (0,2); the covariances are all axis aligned
%	(0.16, 0.64), (0.25, 1) and the identity matrix. The first figure
%	contains a scatter plot of the data.
%
%	A Gaussian mixture model with three components is trained using EM.
%	The parameter vector is printed before training and after training.
%	The user should press any key to continue at these points.  The
%	parameter vector consists of priors (the column), and centres (given
%	as (x, y) pairs as the next two columns).  The diagonal entries of
%	the covariance matrices are printed separately.
%
%	The second figure is a 3 dimensional view of the density function,
%	while the third shows the axes of the 1-standard deviation circles
%	for the three components of the mixture model.
%
%	See also
%	GMM, GMMINIT, GMMEM, GMMPROB, GMMUNPAK
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Generate the data
ndata = 500;

% Fix the seeds for reproducible results
randn('state', 42);
rand('state', 42);
data = randn(ndata, 2);
prior = [0.3 0.5 0.2];
% Mixture model swaps clusters 1 and 3
datap = [0.2 0.5 0.3];
datac = [0 2; 0 0; 2 3.5];
datacov = [1 1;1 0.25; 0.4*0.4 0.8*0.8];
data1 = data(1:prior(1)*ndata,:);
data2 = data(prior(1)*ndata+1:(prior(2)+prior(1))*ndata, :);
data3 = data((prior(1)+prior(2))*ndata +1:ndata, :);

% First cluster has axis aligned variance and centre (2, 3.5)
data1(:, 1) = data1(:, 1)*0.4 + 2.0;
data1(:, 2) = data1(:, 2)*0.8 + 3.5;

% Second cluster has axis aligned variance and centre (0, 0)
data2(:,2) = data2(:, 2)*0.5;

% Third cluster is at (0,2) with identity matrix for covariance
data3 = data3 + repmat([0 2], prior(3)*ndata, 1);

% Put the dataset together again
data = [data1; data2; data3];

clc
disp('This demonstration illustrates the use of a Gaussian mixture model')
disp('with diagonal covariance matrices to approximate the unconditional')
disp('probability density of data in a two-dimensional space.')
disp('We begin by generating the data from a mixture of three Gaussians')
disp('with axis aligned covariance structure and plotting it.')
disp(' ')
disp('The first cluster has centre (0, 2).')
disp('The second cluster has centre (0, 0).')
disp('The third cluster has centre (2, 3.5).')
disp(' ')
disp('Press any key to continue')
pause

fh1 = figure;
plot(data(:, 1), data(:, 2), 'o')
set(gca, 'Box', 'on')

% Set up mixture model
ncentres = 3;
input_dim = 2;
mix = gmm(input_dim, ncentres, 'diag');

options = foptions;
options(14) = 5;	% Just use 5 iterations of k-means in initialisation
% Initialise the model parameters from the data
mix = gmminit(mix, data, options);

% Print out model
disp('The mixture model has three components and diagonal covariance')
disp('matrices.  The model parameters after initialisation using the')
disp('k-means algorithm are as follows')
disp('    Priors        Centres')
disp([mix.priors' mix.centres])
disp('Covariance diagonals are')
disp(mix.covars)
disp('Press any key to continue.')
pause

% Set up vector of options for EM trainer
options = zeros(1, 18);
options(1)  = 1;		% Prints out error values.
options(14) = 20;		% Number of iterations.

disp('We now train the model using the EM algorithm for 20 iterations.')
disp(' ')
disp('Press any key to continue.')
pause

[mix, options, errlog] = gmmem(mix, data, options);

% Print out model
disp(' ')
disp('The trained model has priors and centres:')
disp('    Priors        Centres')
disp([mix.priors' mix.centres])
disp('The data generator has priors and centres')
disp('    Priors        Centres')
disp([datap' datac])
disp('Model covariance diagonals are')
disp(mix.covars)
disp('Data generator covariance diagonals are')
disp(datacov)
disp('Note the close correspondence between these parameters and those')
disp('of the distribution used to generate the data.')
disp(' ')
disp('Press any key to continue.')
pause

clc
disp('We now plot the density given by the mixture model as a surface plot.')
disp(' ')
disp('Press any key to continue.')
pause

% Plot the result
x = -4.0:0.2:5.0;
y = -4.0:0.2:5.0;
[X, Y] = meshgrid(x,y);
X = X(:);
Y = Y(:);
grid = [X Y];
Z = gmmprob(mix, grid);
Z = reshape(Z, length(x), length(y));
c = mesh(x, y, Z);
hold on
title('Surface plot of probability density')
hold off
drawnow

clc
disp('The final plot shows the centres and widths, given by one standard')
disp('deviation, of the three components of the mixture model.  The axes')
disp('of the ellipses of constant density are shown.')
disp(' ')
disp('Press any key to continue.')
pause

% Try to calculate a sensible position for the second figure, below the first
fig1_pos = get(fh1, 'Position');
fig2_pos = fig1_pos;
fig2_pos(2) = fig2_pos(2) - fig1_pos(4);
fh2 = figure('Position', fig2_pos);

h = plot(data(:, 1), data(:, 2), 'bo');
hold on
axis('equal');
title('Plot of data and covariances')
for i = 1:ncentres
  v = [1 0];
  for j = 1:2
    start=mix.centres(i,:)-sqrt(mix.covars(i,:).*v);
    endpt=mix.centres(i,:)+sqrt(mix.covars(i,:).*v);
    linex = [start(1) endpt(1)];
    liney = [start(2) endpt(2)];
    line(linex, liney, 'Color', 'k', 'LineWidth', 3)
    v = [0 1];
  end
  % Plot ellipses of one standard deviation
  theta = 0:0.02:2*pi;
  x = sqrt(mix.covars(i,1))*cos(theta) + mix.centres(i,1);
  y = sqrt(mix.covars(i,2))*sin(theta) + mix.centres(i,2);
  plot(x, y, 'r-');
end
hold off

disp('Note how the data cluster positions and widths are captured by')
disp('the mixture model.')
disp(' ')
disp('Press any key to end.')
pause

close(fh1);
close(fh2);
clear all;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [ ptab, b, c, d, eps, ierror ] =  least_set ( ntab, xtab, ytab, ndeg )

%% LEAST_SET constructs the least squares polynomial approximation to data.
%
%  Discussion:
%
%    The least squares polynomial is not returned directly as a simple
%    polynomial.  Instead, it is represented in terms of a set of
%    orthogonal polynomials appopriate for the given data.  This makes
%    the computation more accurate, but means that the user can not
%    easily evaluate the computed polynomial.  Instead, the routine
%    LEAST_EVAL should be used to evaluate the least squares polynomial
%    at any point.  (However, the value of the least squares polynomial
%    at each of the data points is returned as part of this computation.)
%
%
%    A discrete unweighted inner product is used, so that
%
%      ( F(X), G(X) ) = sum ( 1 <= I <= NTAB ) F(XTAB(I)) * G(XTAB(I)).
%
%    The least squares polynomial is determined using a set of
%    orthogonal polynomials PHI.  These polynomials can be defined
%    recursively by:
%
%      PHI(0)(X) = 1
%      PHI(1)(X) = X - B(1)
%      PHI(I)(X) = ( X - B(I) ) * PHI(I-1)(X) - D(I) * PHI(I-2)(X)
%
%    The array B(1:NDEG) contains the values
%
%      B(I) = ( X*PHI(I-1), PHI(I-1) ) / ( PHI(I-1), PHI(I-1) )
%
%    The array D(2:NDEG) contains the values
%
%      D(I) = ( PHI(I-1), PHI(I-1) ) / ( PHI(I-2), PHI(I-2) )
%
%    Using this basis, the least squares polynomial can be represented as
%
%      P(X)(I) = sum ( 0 <= I <= NDEG ) C(I) * PHI(I)(X)
%
%    The array C(0:NDEG) contains the values
%
%      C(I) = ( YTAB(I), PHI(I) ) / ( PHI(I), PHI(I) )
%
%  Modified:
%
%    16 May 2004
%
%  Reference:
%
%    Gisela Engeln-Muellges and Frank Uhlig,
%    Numerical Algorithms with C, pages 191-193.
%    Springer, 1996.
%
%  Parameters:
%
%    Input, integer NTAB, the number of data points.
%
%    Input, real XTAB(NTAB), the X data.  The values in XTAB
%    should be distinct, and in increasing order.
%
%    Input, real YTAB(NTAB), the Y data values corresponding
%    to the X data in XTAB.
%
%    Input, integer NDEG, the degree of the polynomial which the
%    program is to use.  NDEG must be at least 1, and less than or
%    equal to NTAB-1.
%
%    Output, real PTAB(NTAB), the value of the least squares polynomial
%    at the points XTAB(1:NTAB).
%
%    Output, real B(1:NDEG), C(1:NDEG+1), D(1:NDEG-1), arrays containing
%    data about the polynomial.
%
%    Output, real EPS, the root-mean-square discrepancy of the
%    polynomial fit.
%
%    Output, integer IERROR, error flag.
%    zero, no error occurred;
%    nonzero, an error occurred, and the polynomial could not be computed.
%
  ierror = 0;
  C_OFFSET = 1;
  D_OFFSET = -1;
%
%  Check NDEG.
%
  if ( ndeg < 1 )
    ierror = 1;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LEAST_SET - Fatal error!\n' );
    fprintf ( 1, '  NDEG < 1.\n' );
    error ( 'LEAST_SET - Fatal error!' );
  end

  if ( ntab <= ndeg )
    ierror = 1;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LEAST_SET - Fatal error!\n' );
    fprintf ( 1, '  NTAB <= NDEG.\n' );
    error ( 'LEAST_SET - Fatal error!' );
  end
%
%  Check that the abscissas are strictly increasing.
%
  for ( i = 1 : ntab-1 )
    if ( xtab(i+1) <= xtab(i) )
      ierror = 1;
      fprintf ( 1, '\n' );
      fprintf ( 1, 'LEAST_SET - Fatal error!\n' );
      fprintf ( 1, '  XTAB must be strictly increasing, but\n' );
      fprintf ( 1, '  XTAB(%d) = %f\n', i,   xtab(i)   );
      fprintf ( 1, '  XTAB(%d) = %f\n', i+1, xtab(i+1) );
      error ( 'LEAST_SET - Fatal error!' );
    end
  end

  i0l1 = 0;
  i1l1 = ntab;
%
%  The polynomial is of degree at least 0.
%
  y_sum = sum ( ytab(1:ntab) );
  rn0 = ntab;
  c(0+C_OFFSET) = y_sum / ntab;

  ptab(1:ntab) = y_sum / ntab;

  if ( ndeg == 0 )
    eps = sum ( ( ptab(1:ntab) - ytab(1:ntab) ).^2 );
    eps = sqrt ( eps / ntab );
    b = [];
    d = [];
    return;
  end
%
%  The polynomial is of degree at least 1.
%
  b(1) = sum ( xtab(1:ntab) ) / ntab;

  s = 0.0E+00;
  sum2 = 0.0E+00;
  for ( i = 1 : ntab )
    ztab(i1l1+i) = xtab(i) - b(1);
    s = s + ztab(i1l1+i)^2;
    sum2 = sum2 + ztab(i1l1+i) * ( ytab(i) - ptab(i) );
  end

  rn1 = s;
  c(1+C_OFFSET) = sum2 / s;

  for ( i = 1 : ntab )
    ptab(i) = ptab(i) + c(1+C_OFFSET) * ztab(i1l1+i);
  end

  if ( ndeg == 1 )
    eps = sum ( ( ptab(1:ntab) - ytab(1:ntab) ).^2 );
    eps = sqrt ( eps / ntab );
    d = [];
    return;
  end

  ztable(1:ntab) = 1.0E+00;

  mdeg = 2;
  k = 2;

  while ( 1 )

    d(k+D_OFFSET) = rn1 / rn0;

    sum2 = 0.0E+00;
    for ( i = 1 : ntab )
      sum2 = sum2 + xtab(i) * ztab(i1l1+i) * ztab(i1l1+i);
    end

    b(k) = sum2 / rn1;

    s = 0.0E+00;
    sum2 = 0.0E+00;
    for ( i = 1 : ntab )
      ztab(i0l1+i) = ( xtab(i) - b(k) ) * ztab(i1l1+i) ...
        - d(k+D_OFFSET) * ztab(i0l1+i);
      s = s + ztab(i0l1+i) * ztab(i0l1+i);
      sum2 = sum2 + ztab(i0l1+i) * ( ytab(i) - ptab(i) );
    end

    rn0 = rn1;
    rn1 = s;

    c(k+C_OFFSET) = sum2 / rn1;

    it = i0l1;
    i0l1 = i1l1;
    i1l1 = it;

    for ( i = 1 : ntab )
      ptab(i) = ptab(i) + c(k+C_OFFSET) * ztab(i1l1+i);
    end

    if ( ndeg <= mdeg )
      break;
    end

    mdeg = mdeg + 1;
    k = k + 1;

  end
%
%  Compute the RMS error.
%
  eps = sum ( ( ptab(1:ntab) - ytab(1:ntab) ).^2 );
  eps = sqrt ( eps / ntab );

author: zouchuhang

repo: LayoutNet


text:
function [samples,hists] = SampleDGAnyMarginal(gammas,Lambda,supports,Nsamples)

% [samples,hists]=SampleDGAnyMarginal(gammas,Lambda,supports,Nsamples)
%   Generate samples for a Multivariate Discretized Gaussian with parameters
%   "gammas" and "Lambda" and "supports". The number of samples generated is "Nsamples"
%
%   input and output arguments are as described in "DGAnyMarginal"
%
% Usage:
%
% Code from the paper: 'Generating spike-trains with specified
% correlations', Macke et al., submitted to Neural Computation
%
% www.kyb.mpg.de/bethgegroup/code/efficientsampling

d=size(Lambda,1);

if isempty(supports)
    for k=1:d
        supports{k}=[0:numel(gammas{k})-1];
    end
end

cc=chol(Lambda);

B=randn(Nsamples,d)*cc;

for k=1:d
    [hists{k},dd]=histc(B(:,k),[-inf;gammas{k};inf]);
    hists{k}=hists{k}/Nsamples;
    samples(:,k)=supports{k}(dd);
    hists{k}=hists{k}(1:max(1,end-2));
end

repo: mcbench-benchmarks


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni

%Parseval's Theorem

%x(t)=exp(-t^2)

syms t w
x=exp(-t^2);
Et=int((abs(x))^2,t,-inf,inf) ;
eval(Et)
X=fourier(x,w);
Ew=(1/(2*pi))*int((abs(X))^2,w,-inf,inf);
eval(Ew)

%x(t)=exp(-t)u(t)
x=exp(-t) *heaviside(t);
Et=int((abs(x))^2,t,-inf,inf)
X=fourier(x,w);
Integ= int((abs(X))^2,w,-inf,inf);
Ew=(1/(2*pi))*Integ;
eval(Ew)


text:
function a = minij ( m, n )

%*****************************************************************************80
%
%% MINIJ returns the MINIJ matrix.
%
%  Formula:
%
%    A(I,J) = min ( I, J )
%
%  Example:
%
%    N = 5
%
%    1 1 1 1 1
%    1 2 2 2 2
%    1 2 3 3 3
%    1 2 3 4 4
%    1 2 3 4 5
%
%  Properties:
%
%    A is integral, therefore det ( A ) is integral, and
%    det ( A ) * inverse ( A ) is integral.
%
%    A is positive definite.
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    The inverse of A is tridiagonal.
%
%    The eigenvalues of A are
%
%      LAMBDA(I) = 0.5 / ( 1 - cos ( ( 2 * I - 1 ) * pi / ( 2 * N + 1 ) ) ),
%
%    For N = 12, the characteristic polynomial is
%      P(X) = X**12 - 78 X**11 + 1001 X**10 - 5005 X**9 + 12870 X**8
%        - 19448 X**7 + 18564 X**6 - 11628 X**5 + 4845 X**4 - 1330 X**3
%        + 231 X**2 - 23 X + 1.
%
%    (N+1)*ONES(N) - A also has a tridiagonal inverse.
%
%    Gregory and Karney consider the matrix defined by
%
%      B(I,J) = N + 1 - MAX(I,J)
%
%    which is equal to the MINIJ matrix, but with the rows and
%    columns reversed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Gregory, David Karney,
%    Example 3.12, Example 4.14,
%    A Collection of Matrices for Testing Computational Algorithms,
%    Wiley, 1969, page 41, page 74,
%    LC: QA263.G68.
%
%    Daniel Rutherford,
%    Some continuant determinants arising in physics and chemistry II,
%    Proceedings of the Royal Society Edinburgh,
%    Volume 63, A, 1952, pages 232-241.
%
%    John Todd,
%    Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
%    Academic Press, 1977, page 158.
%
%    Joan Westlake,
%    A Handbook of Numerical Matrix Inversion and Solution of
%    Linear Equations,
%    John Wiley, 1968,
%    ISBN13: 978-0471936756,
%    LC: QA263.W47.
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns
%    of the matrix.
%
%    Output, real A(M,N), the matrix.
%
  for i = 1 : m
    for j = 1 : n
      a(i,j) = min ( i, j );
    end
  end

  return
end


text:
%  Figure 10.11      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
%  fig10_11.m is a script to generate Fig. 10.11
%  the frequency response of the notch network
clf;
nnotch=[1/.81 0  1] ;
dnotch=[1/625 2/25  1];
% define frequency range
w=logspace(-1,1);
w(36)=1;
subplot(211)
w=logspace(-1,1);
w(24)=0.89;
% compute Bode
[magn phn]=bode(nnotch,dnotch,w);
% phn=php-360*ones(phn);
magn1=[magn, ones(size(magn))];
subplot(211); loglog(w,magn1); grid;
xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Fig. 10.11 Bode plot of a notch filter')
subplot(212); semilogx(w,phn); grid;
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');
%Bode grid
bodegrid


text:
% [R,Q] = vgg_rq(S)  Just like qr but the other way around.
%
% If [R,Q] = vgg_rq(X), then R is upper-triangular, Q is orthogonal, and X==R*Q.
% Moreover, if S is a real matrix, then det(Q)>0.

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

% By awf

function [U,Q] = rq(S)

S = S';
[Q,U] = qr(S(end:-1:1,end:-1:1));
Q = Q';
Q = Q(end:-1:1,end:-1:1);
U = U';
U = U(end:-1:1,end:-1:1);

if det(Q)<0
  U(:,1) = -U(:,1);
  Q(1,:) = -Q(1,:);
end

return

author: petercorke

repo: machinevision-toolbox-matlab


text:
function value = zsign1 ( z1, z2 )

%*****************************************************************************80
%
%% ZSIGN1 is a complex transfer-of-sign function.
%
%  Discussion:
%
%    The L1 norm is used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, complex Z1, Z2, the arguments.
%
%    Output, complex VALUE,  a complex value, with the magnitude of
%    Z1, and the argument of Z2.
%
  if ( zabs1 ( z2 ) == 0.0 )
    value = 0.0;
  else
    value = zabs1 ( z1 ) * ( z2 / zabs1 ( z2 ) );
  end

  return
end


text:
function [warnings] = addplotlystartup(startupPaths)
%[1]looks at statup.m files specified by the entries of startupPaths
%[2]appends the addplotly function to startup.m files (if not already present)
%[3]checks for other plotly addpath calls within any startup.m and outputs warning

%output warnings
warnings = cell(size(startupPaths));

for locs = 1:size(startupPaths,1);
    %addpath string for Plotly API
    addString = 'addpath(genpath(fullfile(matlabroot,''toolbox'',''plotly'')),''-end'');';
    %open current startup.m to read
    currentStartupID = fopen(startupPaths{locs},'r');
    if currentStartupID == -1, error('plotly:startupRead',...
            ['\n\nShoot! It looks like something went wrong reading the file: ' ...
            '\n' startupPaths{locs} '\n' ...
            '\nPlease contact your system admin or post a topic on https://community.plotly.com/c/api/matlab/ for more \ninformation. In the ' ...
            'mean time you can add the Plotly API \nto your search path manually whenever you need it! \n\n']);
    end
    %check for any instances of the addplotlyapi function
    startupScan = textscan(currentStartupID, '%s', 'delimiter', '\n', 'whitespace', '');
    startupLines = startupScan{1};
    Index = find(strcmp(startupLines,addString));
    otherPlotlyOccurrence = findstr(fileread(startupPaths{locs}),'plotly');
    %if addString is not in startup.m add it
    if(~any(Index));
        %reopen current startup.m with new permission
        currentStartupID = fopen(startupPaths{locs},'a+');
        if currentStartupID == -1, error('plotly:startupWrite',...
            ['\n\nShoot! It looks like something went wrong writing to the file: ' ...
            '\n\n' startupPaths{locs} '\n' ...
            '\nPlease contact your system admin or post a topic on https://community.plotly.com/c/api/matlab/ for more \ninformation. In the ' ...
            'mean time you can add the Plotly API \nto your search path manually whenever you need it! \n']);
        end
        fprintf(currentStartupID,['\n' addString]);
    end
    if(length(Index) ~= length(otherPlotlyOccurrence));
        warnings{locs} = ['\n[WARNING]: \n\nWe found an addpath specification for another version of Plotly at: ' ...
                                '\n\n' startupPaths{locs} '\n\nyou may be forcing MATLAB to look for an older version of Plotly!\n\n'];
    else
        warnings{locs} = '';
    end
    fclose(currentStartupID);
end
end

author: plotly


text:
%--- help for binary_operation ---
%
%  INTERNAL FUNCTION
%
%

author: jmaih