text:
function S = globMatrix3DStiff(fun,mesh,fem1,fem2)

%% USAGE: generate stiffness global matrix on a 3D mesh
%
% INPUTS:
% fun --- coefficient function
% mesh --- a struct data contains very rich mesh information.
% fem1 --- global DoF for test function space
% fem2 --- global DoF for trial function space
%
% OUTPUTS:
% [IN JN XN] --- triplets of the sparse matrix from regular elements.
% [II JI XI] --- triplets of the sparse matrix from interface elements.

% Last Modified: 08/07/2020 by Xu Zhang

%% 0. Initializaiton
if strcmp(fem1.type,'P1')||strcmp(fem1.type,'DGP1')||strcmp(fem1.type,'CR')
    feEvalBas1 = @evalP1Bas3D;
elseif strcmp(fem1.type,'P2')||strcmp(fem1.type,'DGP2')
    feEvalBas1 = @evalP2Bas3D;
end

if strcmp(fem2.type,'P1')||strcmp(fem2.type,'DGP1')||strcmp(fem2.type,'CR')
    feEvalBas2 = @evalP1Bas3D;
elseif strcmp(fem2.type,'P2')||strcmp(fem2.type,'DGP2')
    feEvalBas2 = @evalP2Bas3D;
end

dof1 = fem1.ldof; dof2 = fem2.ldof; nloc = dof1*dof2;
nt = length(mesh.t);
A = fem1.area; gx = fem1.gx; gy = fem1.gy; gz = fem1.gz; gw = fem1.gw;
X = zeros(nloc*nt, 1);

coef = feval(fun,gx,gy,gz);
Ibasx = cell(dof1,1); Ibasy = cell(dof1,1); Ibasz = cell(dof1,1);
Jbasx = cell(dof2,1); Jbasy = cell(dof2,1); Jbasz = cell(dof2,1);
for i = 1:dof1
    Ibasx{i} = feEvalBas1(fem1.bas(:,:,i), gx, gy, gz, [1,0,0]);
    Ibasy{i} = feEvalBas1(fem1.bas(:,:,i), gx, gy, gz, [0,1,0]);
    Ibasz{i} = feEvalBas1(fem1.bas(:,:,i), gx, gy, gz, [0,0,1]);
end
for j = 1:dof2
    Jbasx{j} = feEvalBas2(fem2.bas(:,:,j), gx, gy, gz, [1,0,0]);
    Jbasy{j} = feEvalBas2(fem2.bas(:,:,j), gx, gy, gz, [0,1,0]);
    Jbasz{j} = feEvalBas2(fem2.bas(:,:,j), gx, gy, gz, [0,0,1]);
end

I = reshape(repmat(fem1.t,4,1),nloc*nt,1);
J = repmat(reshape(fem2.t,dof2*nt,1),4,1);
ind = 0;
for i = 1:dof1
    for j = 1:dof2
        X(ind+1:ind+nt) = A.*(sum(((Ibasx{i}.*(coef.*Jbasx{j})).*gw'),2) + ...
            sum(((Ibasy{i}.*(coef.*Jbasy{j})).*gw'),2) + ...
            sum(((Ibasz{i}.*(coef.*Jbasz{j})).*gw'),2));
        ind = ind + nt;
    end
end
ID = find(X~=0);
S = sparse(I(ID),J(ID),X(ID),size(fem1.p,1),size(fem2.p,1));

repo: ifem


text:
function [CnoCom, Com] = PurgeComments (C, ch)
%
%   [CnoCom, Com] = PurgeComments (C, [ch])
%
%Purpose:
%   purges all matlab coments from character array
%
%
%Input Parameters:
%   C : character array
%   ch : is the character indicating a comment line, default is '%'
%
%Output Parameters:
%   err : 0 No Problem
%       : 1 Mucho Problems
%   CnoCom : comment purged caharacter array
%   Com : Comment array
%
%More Info :
%
%  SkipMatlabHelp
%  FindChar
%	NextString
%
%   C = sprintf (...
%  '%%Hello \n%%Couci Couci \nAbc = 00\n%%Wald\njon\n%%ert\nMlk\n%%brt%%crt\nflp\n%%kr\n\nko=6\n%%bizz\n%%Bozz\n\nfly')
%   [CnoCom] = PurgeComments (C)
%
%
%     Author : Ziad Saad
%     Date : Sat Mar 27 15:42:11 CST 1999

%Define the function name for easy referencing
FuncName = 'PurgeComments';

if (nargin == 1), ch = '%'; end

%Debug Flag
DBG = 1;

%initailize return variables
err = 1;
CnoCom = [];
Com = '';

%skip the help lines (top comments) if there are any
if (C(1) == ch),
	[C] = SkipMatlabHelp (C);
end

[err, Loc_1] = FindChar (C, ch);

if (~isempty(Loc_1)),
	cnt = 1;
	CnoCom = C;
	nC = length(C);
	%fprintf ('C purged of heading is :\n');
	%fprintf ('%c', C);
	FirstChunk = C(1:Loc_1(1)-1);
	CnoCom(1:length(FirstChunk)) = FirstChunk;
	nCnoCom = length(FirstChunk);
	%fprintf ('CnoCom at first chunk : \n');
	%fprintf ('%c', CnoCom(1:nCnoCom));

	while (cnt < length(Loc_1)),
		%skip that line
		[err,Cnext, cend] = NextString (C, 'NewLine', Loc_1(cnt));
	   Com = [Com Cnext];
		%make sure next line is not a comment
		if (cend+1 == Loc_1(cnt+1)), lop = 1; else lop = 0; end
		while (lop), %looks like it is a comment, cycle through block
			% fprintf (1,'Looping, cnt = %g/%g\n', cnt, length(Loc_1));
			cnt = cnt+1;
			[err,Cnext, cend] = NextString (C, 'NewLine', Loc_1(cnt));
			Com = [Com Cnext];
         lop = 0;
			if (cnt < length(Loc_1))
				if (cend+1 == Loc_1(cnt+1)),
					lop = 1;
				end
			end
		end

		%copy the next good piece,
		%if you did not run out of cnt in the looping condition above
		if (cnt < length(Loc_1)),
			GoodChunk = C(cend+1:Loc_1(cnt+1)-1);
			CnoCom(1+nCnoCom:nCnoCom+length(GoodChunk)) = GoodChunk;
			nCnoCom = nCnoCom+length(GoodChunk);
			%fprintf ('CnoCom at middle chunk (cnt = %g) : \n', cnt);
			%fprintf ('%c', CnoCom(1:nCnoCom));

			%increment and go on
			cnt = cnt + 1;
		end
	end

		%append last piece if any
		if (Loc_1(cnt) < nC),
			[err,Cnext, cend] = NextString (C, 'NewLine', Loc_1(cnt));
			Com = [Com Cnext];
         LastChunk = C(cend+1:nC);
			CnoCom(1+nCnoCom:nCnoCom+length(LastChunk)) = LastChunk;
			nCnoCom = nCnoCom+length(LastChunk);
			%fprintf ('Last Chunk : \n');
			%fprintf ('%c',LastChunk);
		end

		CnoCom = CnoCom(1:nCnoCom);

		%fprintf ('\n***\nFinal CnoCom\n***\n');
		%fprintf ('%c', CnoCom(1:nCnoCom));
else %no comments
	CnoCom = C;
end
%Com
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function PoiseuilleFlow( varargin )
%POISOILLEFLOW Summary of this function goes here
%   Detailed explanation goes here
close all; % closes all figures
addpath('../Functions');

h=4; % height of the channel
dPdx=-0.01; % pressure gradient in x direction
mu=0.001; % dynamic viscosity

% function for the analytical solution
u_analFunc=@(y) -0.5*h^2/mu*dPdx*(1-(y/h).^2);

% create mesh
m=createMesh2D(20,20,2*h,2*h);

%% Boundary Conditions
U_BC = createBC(m); % all Neumann boundary condition structure

% Assign boundary conditions
U_BC.left.a(:) = 1; U_BC.left.b(:)=0; U_BC.left.c(:)=0; % Dirichlet for the left boundary
U_BC.right.a(:) = 1; U_BC.right.b(:)=0; U_BC.right.c(:)=0; % Dirichlet for the left boundary
U_BC.top.a(:) = 0; U_BC.top.b(:)=1; U_BC.top.c(:)=0; % Dirichlet for the left boundary
U_BC.bottom.a(:) = 0; U_BC.bottom.b(:)=1; U_BC.bottom.c(:)=0; % Dirichlet for the left boundary

V_BC = createBC(m); % all Neumann boundary condition structure

% Assign boundary conditions
V_BC.top.a(:) = 0; V_BC.top.b(:)=1; V_BC.top.c(:)=0; % Dirichlet for the left boundary
V_BC.bottom.a(:) = 0; V_BC.bottom.b(:)=1; V_BC.bottom.c(:)=0; % Dirichlet for the left boundary
V_BC.left.a(:) = 0; V_BC.left.b(:)=1; V_BC.left.c(:)=0; % Dirichlet for the left boundary
V_BC.right.a(:) = 0; V_BC.right.b(:)=1; V_BC.right.c(:)=0; % Dirichlet for the left boundary

BC_pressureCorrection = createBC(m);
BC_pressureCorrection.left.a(:) = 1; BC_pressureCorrection.left.b(:)=0; BC_pressureCorrection.left.c(:)=0;
BC_pressureCorrection.right.a(:) = 1; BC_pressureCorrection.right.b(:)=0; BC_pressureCorrection.right.c(:)=0;
BC_pressureCorrection.bottom.a(:) = 1; BC_pressureCorrection.bottom.b(:)=0; BC_pressureCorrection.bottom.c(:)=0;
BC_pressureCorrection.bottom.a(end/2) = 0; BC_pressureCorrection.bottom.b(end/2)=1; BC_pressureCorrection.bottom.c(end/2)=1;
BC_pressureCorrection.top.a(:) = 1; BC_pressureCorrection.top.b(:)=0; BC_pressureCorrection.top.c(:)=0;

% constant pressure - so there is no gradient
p_cellVar=createCellVariable(m,1,BC_pressureCorrection);

mu_cellVar=createCellVariable(m,0.001);
mu_faceVar=createFaceVariable(m,0.001);

rho=1000;
rho_cellVar=createCellVariable(m,rho);
rho_faceVar=createFaceVariable(m,rho);

U_cellVar=createCellVariable(m,0,U_BC);
V_cellVar=createCellVariable(m,0,V_BC);
faceVelocity=createFaceVariable(m,[0 0]);

pGradMat = gradientCellTerm(p_cellVar);
helpVarx = createCellVariable(m,dPdx);
helpVary = createCellVariable(m,0);
pGradMat.xvalue=helpVarx.value(2:end-1,2:end-1);
pGradMat.yvalue=helpVary.value(2:end-1,2:end-1);

% solve momentum equation Momentum equation
[U_cellVar,~,~,~]=MomentumEq(1,m,pGradMat,mu_faceVar,faceVelocity,U_cellVar,V_cellVar,U_BC,V_BC,rho,'SIMPLE');

y_anal=linspace(0,h,100);
u_anal=u_analFunc(y_anal);

figure;
plot(u_anal,y_anal); % plot analytical solution
hold on;
plot(U_cellVar.value(3,2:end-1),U_cellVar.domain.cellcenters.y-h,'o');
hold off;
figure;
visualizeCells(U_cellVar);

end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
% IMAGE_COMPLEX Return a color image that for a complex 2d signal.
%
% Usage
%	x_rgb = IMAGE_COMPLEX(x)
%
% Input
%    x (complex 2d matrix): the signals to be displayed
%
% Output
%    x_rgb (numeric): a color image of the signal.
%
% Description
%    Compute a color image corresponding to a complex 2d signal.
%    The saturation corresponds to the complex modulus while the hue
%    corresponds to the complex argument.
%
% See also
%    DISPLAY_FILTER_BANK_2D_COMPLEX, IMAGE_COMPLEX_CELL
function x_rgb = image_complex(x)
	x = x / max(abs(x(:)) + 1e-10);
	x_hsv(:,:,1) = (pi + 1e-5 + angle(x)) / (2*pi);
	x_hsv(:,:,2) = abs(x);
	x_hsv(:,:,3) = 1;
	x_rgb = hsv2rgb(x_hsv);
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
% SCRIPT TEST FOR THE KINEMATIC PROBLEM FOR SERIAL ROBOTS

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE 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.
%
% ARTE 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 ARTE.  If not, see <http://www.gnu.org/licenses/>.

close all

fprintf('\nTHE DEMO PRESENTS THE DIRECT AND INVERSE KINEMATIC PROBLEM')

%there are eight possible solutions for the inverse kinematic problem for most of these robots
n_solutions = 8;

%Try different configurations beware that, depending on the robot's topology
%not all the eight possible solutions will be feasible for an antropomorphic 6R robot.
%q=[0.5 -0.4 -0.2 0.1 0.1 0.1]

q=[0.2 -0.2 0.3 0.1 0 0.1]

%load robot parameters. You can try different robots%
robot=load_robot('ABB', 'IRB140'); n_solutions = 8;
%robot=load_robot('KUKA', 'KR60_3'); n_solutions = 8;
%robot=load_robot('ABB', 'IRB120'); n_solutions = 8;
%robot=load_robot('MOTOMAN', 'MH12'); n_solutions = 8;

%adjust 3D view as desired
adjust_view(robot)

%there are just 2 solutions for these robots and 4 DOF
%q = [pi/2 0.2 0.8 pi/4]
%robot=load_robot('kuka', 'KR5_scara_R350_Z200'); n_solutions = 2;
%robot=load_robot('example', 'scara'); n_solutions = 2;
%robot=load_robot('example', '2dofplanar'); n_solutions = 2;

%draw the robot
drawrobot3d(robot, q)

%Now compute direct kinematics for this position q
T = directkinematic(robot, q)

%Set to zero if you want to see the robot transparent
robot.graphical.draw_transparent=0;

%Set to one if you want to see the DH axes
%robot.graphical.draw_axes=1;

%Call the inversekinematic for this robot. All the possible solutions are
%stored at qinv. At least, one of the possible solutions should match q
qinv = inversekinematic(robot, T)

fprintf('\nNOW WE CAN REPRESENT THE DIFFERENT SOLUTIONS TO ACHIEVE THE SAME POSITION AND ORIENTATION\n')
fprintf('\nNote that some solutions may not be feasible since some joints may be out of range.\n')
correct=zeros(1,n_solutions);
%check that all of them are possible solutions!
for i=1:size(qinv,2),

    Ti = directkinematic(robot, qinv(:,i)) %Ti is constant for the different solutions

    % Note that all the solutions may not be feasible. Some of the joints may
    % be out of range. You can test this situation with test_joints
    test_joints(robot, qinv(:,i));

    %now draw the robot to see the solution
    drawrobot3d(robot, qinv(:,i))

    pause(1);

    k=sum(sum((T-Ti).^2));
    if k < 0.01 % a simple threshold to find differences in the solution
        correct(1,i)= 1;
    else
        correct(1,i)= 0; %uncorrect solution
        fprintf('\nERROR: One of the solutions seems to be uncorrect. Sum of errors: %f', i, k);
    end
end

fprintf('\n************** RESULTS **************')

%Display a message if any of the solutions is not correct
if sum(correct)==n_solutions
    fprintf('\nTEST 1--> OK: Every solution in qinv yields the same position/orientation T');
else
    fprintf('\nTEST 1--> ERROR: One or more of the solutions seem to be uncorrect.');
end

%Now, test if any of the solutions in qinv matches q
%find the solution that matches the initial q
%delta is just a squared sum of errors at each of the columns of the matrix
%which store the different solutions of qinv
delta=(repmat(q',[1 n_solutions])-qinv).^2;
i=find(sum(delta,1) < 0.01);
if ~isempty(i)
    fprintf('\nTEST 2--> OK!: Found a matching solution for the initial q.\n');
    solution=qinv(:,i)
else
    error_test2=1
    fprintf('\nTEST 2--> ERROR: Did not find a matching solution for the initial q.');
end

fprintf('\n************** ****** **************\n')

repo: ARTE


text:
function [X, fX, i] = fmincg(f, X, options, P1, P2, P3, P4, P5)
% Minimize a continuous differentialble multivariate function. Starting point
% is given by "X" (D by 1), and the function named in the string "f", must
% return a function value and a vector of partial derivatives. The Polack-
% Ribiere flavour of conjugate gradients is used to compute search directions,
% and a line search using quadratic and cubic polynomial approximations and the
% Wolfe-Powell stopping criteria is used together with the slope ratio method
% for guessing initial step sizes. Additionally a bunch of checks are made to
% make sure that exploration is taking place and that extrapolation will not
% be unboundedly large. The "length" gives the length of the run: if it is
% positive, it gives the maximum number of line searches, if negative its
% absolute gives the maximum allowed number of function evaluations. You can
% (optionally) give "length" a second component, which will indicate the
% reduction in function value to be expected in the first line-search (defaults
% to 1.0). The function returns when either its length is up, or if no further
% progress can be made (ie, we are at a minimum, or so close that due to
% numerical problems, we cannot get any closer). If the function terminates
% within a few iterations, it could be an indication that the function value
% and derivatives are not consistent (ie, there may be a bug in the
% implementation of your "f" function). The function returns the found
% solution "X", a vector of function values "fX" indicating the progress made
% and "i" the number of iterations (line searches or function evaluations,
% depending on the sign of "length") used.
%
% Usage: [X, fX, i] = fmincg(f, X, options, P1, P2, P3, P4, P5)
%
% See also: checkgrad
%
% Copyright (C) 2001 and 2002 by Carl Edward Rasmussen. Date 2002-02-13
%
%
% (C) Copyright 1999, 2000 & 2001, Carl Edward Rasmussen
%
% Permission is granted for anyone to copy, use, or modify these
% programs and accompanying documents for purposes of research or
% education, provided this copyright notice is retained, and note is
% made of any changes that have been made.
%
% These programs and documents are distributed without any warranty,
% express or implied.  As the programs were written for research
% purposes only, they have not been tested to the degree that would be
% advisable in any important application.  All use of these programs is
% entirely at the user's own risk.
%
% [ml-class] Changes Made:
% 1) Function name and argument specifications
% 2) Output display
%

% Read options
if exist('options', 'var') && ~isempty(options) && isfield(options, 'MaxIter')
    length = options.MaxIter;
else
    length = 100;
end

RHO = 0.01;                            % a bunch of constants for line searches
SIG = 0.5;       % RHO and SIG are the constants in the Wolfe-Powell conditions
INT = 0.1;    % don't reevaluate within 0.1 of the limit of the current bracket
EXT = 3.0;                    % extrapolate maximum 3 times the current bracket
MAX = 20;                         % max 20 function evaluations per line search
RATIO = 100;                                      % maximum allowed slope ratio

argstr = ['feval(f, X'];                      % compose string used to call function
for i = 1:(nargin - 3)
  argstr = [argstr, ',P', int2str(i)];
end
argstr = [argstr, ')'];

if max(size(length)) == 2, red=length(2); length=length(1); else red=1; end
S=['Iteration '];

i = 0;                                            % zero the run length counter
ls_failed = 0;                             % no previous line search has failed
fX = [];
[f1 df1] = eval(argstr);                      % get function value and gradient
i = i + (length<0);                                            % count epochs?!
s = -df1;                                        % search direction is steepest
d1 = -s'*s;                                                 % this is the slope
z1 = red/(1-d1);                                  % initial step is red/(|s|+1)

while i < abs(length)                                      % while not finished
  i = i + (length>0);                                      % count iterations?!

  X0 = X; f0 = f1; df0 = df1;                   % make a copy of current values
  X = X + z1*s;                                             % begin line search
  [f2 df2] = eval(argstr);
  i = i + (length<0);                                          % count epochs?!
  d2 = df2'*s;
  f3 = f1; d3 = d1; z3 = -z1;             % initialize point 3 equal to point 1
  if length>0, M = MAX; else M = min(MAX, -length-i); end
  success = 0; limit = -1;                     % initialize quanteties
  while 1
    while ((f2 > f1+z1*RHO*d1) || (d2 > -SIG*d1)) && (M > 0)
      limit = z1;                                         % tighten the bracket
      if f2 > f1
        z2 = z3 - (0.5*d3*z3*z3)/(d3*z3+f2-f3);                 % quadratic fit
      else
        A = 6*(f2-f3)/z3+3*(d2+d3);                                 % cubic fit
        B = 3*(f3-f2)-z3*(d3+2*d2);
        z2 = (sqrt(B*B-A*d2*z3*z3)-B)/A;       % numerical error possible - ok!
      end
      if isnan(z2) || isinf(z2)
        z2 = z3/2;                  % if we had a numerical problem then bisect
      end
      z2 = max(min(z2, INT*z3),(1-INT)*z3);  % don't accept too close to limits
      z1 = z1 + z2;                                           % update the step
      X = X + z2*s;
      [f2 df2] = eval(argstr);
      M = M - 1; i = i + (length<0);                           % count epochs?!
      d2 = df2'*s;
      z3 = z3-z2;                    % z3 is now relative to the location of z2
    end
    if f2 > f1+z1*RHO*d1 || d2 > -SIG*d1
      break;                                                % this is a failure
    elseif d2 > SIG*d1
      success = 1; break;                                             % success
    elseif M == 0
      break;                                                          % failure
    end
    A = 6*(f2-f3)/z3+3*(d2+d3);                      % make cubic extrapolation
    B = 3*(f3-f2)-z3*(d3+2*d2);
    z2 = -d2*z3*z3/(B+sqrt(B*B-A*d2*z3*z3));        % num. error possible - ok!
    if ~isreal(z2) || isnan(z2) || isinf(z2) || z2 < 0 % num prob or wrong sign?
      if limit < -0.5                               % if we have no upper limit
        z2 = z1 * (EXT-1);                 % the extrapolate the maximum amount
      else
        z2 = (limit-z1)/2;                                   % otherwise bisect
      end
    elseif (limit > -0.5) && (z2+z1 > limit)         % extraplation beyond max?
      z2 = (limit-z1)/2;                                               % bisect
    elseif (limit < -0.5) && (z2+z1 > z1*EXT)       % extrapolation beyond limit
      z2 = z1*(EXT-1.0);                           % set to extrapolation limit
    elseif z2 < -z3*INT
      z2 = -z3*INT;
    elseif (limit > -0.5) && (z2 < (limit-z1)*(1.0-INT))  % too close to limit?
      z2 = (limit-z1)*(1.0-INT);
    end
    f3 = f2; d3 = d2; z3 = -z2;                  % set point 3 equal to point 2
    z1 = z1 + z2; X = X + z2*s;                      % update current estimates
    [f2 df2] = eval(argstr);
    M = M - 1; i = i + (length<0);                             % count epochs?!
    d2 = df2'*s;
  end                                                      % end of line search

  if success                                         % if line search succeeded
    f1 = f2; fX = [fX' f1]';
    fprintf('%s %4i | Cost: %4.6e\r', S, i, f1);
    s = (df2'*df2-df1'*df2)/(df1'*df1)*s - df2;      % Polack-Ribiere direction
    tmp = df1; df1 = df2; df2 = tmp;                         % swap derivatives
    d2 = df1'*s;
    if d2 > 0                                      % new slope must be negative
      s = -df1;                              % otherwise use steepest direction
      d2 = -s'*s;
    end
    z1 = z1 * min(RATIO, d1/(d2-realmin));          % slope ratio but max RATIO
    d1 = d2;
    ls_failed = 0;                              % this line search did not fail
  else
    X = X0; f1 = f0; df1 = df0;  % restore point from before failed line search
    if ls_failed || i > abs(length)          % line search failed twice in a row
      break;                             % or we ran out of time, so we give up
    end
    tmp = df1; df1 = df2; df2 = tmp;                         % swap derivatives
    s = -df1;                                                    % try steepest
    d1 = -s'*s;
    z1 = 1/(1-d1);
    ls_failed = 1;                                    % this line search failed
  end
  if exist('OCTAVE_VERSION')
    fflush(stdout);
  end
end
fprintf('\n');

author: Ayatans

repo: Machine-Learning-homework


text:
function af = naca4gen(iaf)
%
% "naca4gen" Generates the NACA 4 digit airfoil coordinates with desired no.
% of panels (line elements) on it.
%      Author : Divahar Jayaraman (j.divahar@yahoo.com)
%
% INPUTS-------------------------------------------------------------------
%       iaf.designation = NACA 4 digit iaf.designation (eg. '2412') - STRING !
%                 iaf.n = no of panels (line elements) PER SIDE (upper/lower)
% iaf.HalfCosineSpacing = 1 for "half cosine x-spacing"
%                       = 0 to give "uniform x-spacing"
%          iaf.wantFile = 1 for creating airfoil data file (eg. 'naca2412.dat')
%                       = 0 to suppress writing into a file
%       iaf.datFilePath = Path where the data  file has to be created
%                         (eg. 'af_data_folder/naca4digitAF/')
%                         use only forward slash '/' (Just for OS portability)
%
% OUTPUTS------------------------------------------------------------------
% Data:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
%       af.x = x cordinate (nx1 array)
%       af.z = z cordinate (nx1 array)
%      af.xU = x cordinate of upper surface (nx1 array)
%      af.zU = z cordinate of upper surface (nx1 array)
%      af.xL = x cordinate of lower surface (nx1 array)
%      af.zL = z cordinate of lower surface (nx1 array)
%      af.xC = x cordinate of camber line (nx1 array)
%      af.zC = z cordinate of camber line (nx1 array)
%    af.name = Name of the airfoil
%  af.header = Airfoil name ; No of panels ; Type of spacing
%              (eg. 'NACA4412 : [50 panels,Uniform x-spacing]')
%
%
% File:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
% First line : Header eg. 'NACA4412 : [50 panels,Half cosine x-spacing]'
% Subsequent lines : (2*iaf.n+1) rows of x and z values
%
% Typical Inputs:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
% iaf.designation='2312';
% iaf.n=56;
% iaf.HalfCosineSpacing=1;
% iaf.wantFile=1;
% iaf.datFilePath='./'; % Current folder
% iaf.is_finiteTE=0;

% % [[Calculating key parameters-----------------------------------------]]
t=str2num(iaf.designation(3:4))/100;
m=str2num(iaf.designation(1))/100;
p=str2num(iaf.designation(2))/10;

a0= 0.2969;
a1=-0.1260;
a2=-0.3516;
a3= 0.2843;

if iaf.is_finiteTE ==1
    a4=-0.1015; % For finite thick TE
else
    a4=-0.1036;  % For zero thick TE
end

% % [[Giving x-spacing---------------------------------------------------]]
if iaf.HalfCosineSpacing==1
    beta=linspace(0,pi,iaf.n+1)';
    x=(0.5*(1-cos(beta))); % Half cosine based spacing
    iaf.header=['NACA' iaf.designation ' : [' num2str(2*iaf.n) 'panels,Half cosine x-spacing]'];
else
    x=linspace(0,1,iaf.n+1)';
    iaf.header=['NACA' iaf.designation ' : [' num2str(2*iaf.n) 'panels,Uniform x-spacing]'];
end

yt=(t/0.2)*(a0*sqrt(x)+a1*x+a2*x.^2+a3*x.^3+a4*x.^4);

xc1=x(find(x<=p));
xc2=x(find(x>p));
xc=[xc1 ; xc2];

if p==0
    xu=x;
    yu=yt;

    xl=x;
    yl=-yt;

    zc=zeros(size(xc));
else
    yc1=(m/p^2)*(2*p*xc1-xc1.^2);
    yc2=(m/(1-p)^2)*((1-2*p)+2*p*xc2-xc2.^2);
    zc=[yc1 ; yc2];

    dyc1_dx=(m/p^2)*(2*p-2*xc1);
    dyc2_dx=(m/(1-p)^2)*(2*p-2*xc2);
    dyc_dx=[dyc1_dx ; dyc2_dx];
    theta=atan(dyc_dx);

    xu=x-yt.*sin(theta);
    yu=zc+yt.*cos(theta);

    xl=x+yt.*sin(theta);
    yl=zc-yt.*cos(theta);
end
af.name=['NACA ' iaf.designation];

af.x=[flipud(xu) ; xl(2:end)];
af.z=[flipud(yu) ; yl(2:end)];

indx1=1:min( find(af.x==min(af.x)) );  % Upper surface indices
indx2=min( find(af.x==min(af.x)) ):length(af.x); % Lower surface indices
af.xU=af.x(indx1); % Upper Surface x
af.zU=af.z(indx1); % Upper Surface z
af.xL=af.x(indx2); % Lower Surface x
af.zL=af.z(indx2); % Lower Surface z

af.xC=xc;
af.zC=zc;

lecirFactor=0.8;
af.rLE=0.5*(a0*t/0.2)^2;

le_offs=0.5/100;
dyc_dx_le=(m/p^2)*( 2*p-2*le_offs );
theta_le=atan(dyc_dx_le);
af.xLEcenter=af.rLE*cos(theta_le);
af.yLEcenter=af.rLE*sin(theta_le);

% % [[Writing iaf data into file------------------------------------------]]
if iaf.wantFile==1
    F1=iaf.header;
    F2=num2str([af.x af.z]);
    F=strvcat(F1,F2);
    fileName=[iaf.datFilePath 'naca' iaf.designation '.dat'];
    dlmwrite(fileName,F,'delimiter','')
end

author: Sable

repo: mcbench-benchmarks


text:
function [bbci, data]= bbci_calibrate_ERPSpeller(bbci, data)
%BBCI_CALIBRATE_ERP_SPELLER - Calibrate online system for ERP Spellers
%
%This function is called by bbci_calibrate
%(if BBCI.calibate.fcn is set to @bbci_calibrate_ERPSpeller).
%Via BBCI.calibrate.settings, the details can be specified, se below.
%
%Synopsis:
% [BBCI, DATA]= bbci_calibrate_ERPSpeller(BBCI, DATA)
%
%Arguments:
%  BBCI -  the field 'calibrate' holds parameters specific to ERP Spellers.
%  DATA -  holds the calibration data
%
%Output:
%  BBCI - Updated BBCI structure in which all necessary fields for
%     online operation are set, see bbci_apply_structures.
%  DATA - As input but added some information of the analysis that might
%     be reused in a second run
%
%BBCI.calibrate.settings may include the following parameters:
%  ...to.be.filled.in...
%  'cue_markers': Markers that trigger the classification of an ERP,
%          default: [11:16,21:26,31:36,41:46]

% 11-2011 Benjamin Blankertz

default_grd= sprintf(['scale,FC3,FC1,FCz,FC2,FC4,legend\n' ...
                     'C5,C3,C1,Cz,C2,C4,C6\n' ...
                     'CP5,CP3,CP1,CPz,CP2,CP4,CP6\n' ...
                     'P7,P5,P3,Pz,P4,P6,P8\n' ...
                     'PO7,PO3,O1,Oz,O2,PO4,PO8']);

props= {'ref_ival'      [-200 0]                       '!DOUBLE[1 2]'
        'disp_ival'     [-200 800]                     '!DOUBLE[1 2]'
        'clab'          {'not','E*'}                   'CHAR|CELL{CHAR}'
				'clab_rereference'            ''               'CHAR|CELL{CHAR}'
        'cfy_clab'      {'not','E*','Fp*','A*'}        'CHAR|CELL{CHAR}'
        'cfy_ival'                    'auto'           'CHAR(auto)|DOUBLE[- 2]'
        'cfy_ival_pick_peak'          [100 700]        'DOUBLE[1 2]'
        'band'                        []               'DOUBLE[1 2]|DOUBLE[1]'
        'control_per_stimulus'        0                'BOOL'
        'model'             @train_RLDAshrink          'FUNC|CELL'
				'whitening'                   0                'BOOL'
        'nSequences'                  5                'INT'
        'nClasses'                    0                'INT'
        'cue_markers'       [11:16,21:26,31:36,41:46]  '!DOUBLE[1 -]'
        'cfy_maxival'                 []               'DOUBLE[1 2]'
        'reject_artifacts'            1                'BOOL'
        'reject_artifacts_opts'       {}               'PROPLIST'
        'reject_channels'             1                'BOOL'
        'reject_eyemovements'         0                'BOOL'
        'reject_eyemovements_crit'    []               'STRUCT'
        'grd'                         default_grd      'CHAR'
        'clab_erp'                    {'Cz','PO7'}     'CHAR|CELL{CHAR}'
        'clab_rsq'                    {}               'CHAR|CELL{CHAR}'
        'target_dist'                 0                'INT'
        'mrk2feedback_fcn'   @(x)(1+mod(x-1,10))       'FUNC'
        'create_figs'                 1                'BOOL'
       };
[opt, isdefault]= opt_setDefaults('bbci.calibrate.settings', props);

nClassesGuess= length(unique(opt.mrk2feedback_fcn(opt.cue_markers)));
[opt, isdefault]= ...
  opt_overrideIfDefault(opt, isdefault, 'nClasses', nClassesGuess);
[opt, isdefault]= ...
  opt_overrideIfDefault(opt, isdefault, 'clab_rsq', opt.clab_erp);

default_crit_ival= [100 800];
default_crit_ival(1)= max(opt.disp_ival(1), default_crit_ival(1));
default_crit_ival(2)= min(opt.disp_ival(2), default_crit_ival(2));
default_crit= {'maxmin'  100                        'DOUBLE[1]'
               'clab'    {{'EOG*','F9,10','Fp*'}}   'CHAR|CELL{CHAR}'
               'ival'    default_crit_ival          'DOUBLE[1 2]'
              };
[opt, isdefault]= ...
  opt_overrideIfDefault(opt, isdefault, 'reject_eyemovements_crit', ...
                        opt_propspecToStruct(default_crit));

% add possibly missing fields to opt.reject_eyemovements_crit
opt.reject_eyemovements_crit= ...
  opt_setDefaults(opt.reject_eyemovements_crit, default_crit);

[opt, isdefault]= ...
  opt_overrideIfDefault(opt, isdefault, ...
                        'clab_rsq', opt.clab_erp);

% store chosen default settings back in bbci variable
bbci.calibrate.settings= opt;

opt_scalp_erp= defopt_scalp_erp('colorOrder', [0.9 0 0.9; 0.4 0.4 0.4]);
opt_scalp_r= defopt_scalp_r;

mnt= mnt_setGrid(data.mnt, opt.grd);

if ~data.isnew && isfield(data, 'result'),
  previous= data.result;
else
  previous= struct;
end

if ~isempty(opt.band),
  [filt_b,filt_a]= cheby2(5, 20, opt.band/data.cnt.fs*2);
  data.cnt = proc_filt(data.cnt, filt_b, filt_a);
end

BC_result= [];
BC_result.mrk= data.mrk;

BC_result.clab= data.cnt.clab(util_chanind(data.cnt, opt.clab));
BC_result.cfy_clab= data.cnt.clab(util_chanind(data.cnt, opt.cfy_clab));
% The following figures are always generated
if opt.create_figs,
    data.figure_handles= [1 2 4 6];
end

if opt.target_dist,
  [BC_result.mrk, BC_result.selected_trials]= ...
      mrk_selectTargetDist(BC_result.mrk, opt.target_dist);
else
  BC_result.selected_trials= NaN;
end

if ~isempty(opt.clab_rereference),
	[data.cnt, dmy]= proc_commonAverageReference(data.cnt, opt.clab_rereference, '*');
	[dmy, A_reref]= proc_commonAverageReference(proc_selectChannels(data.cnt,BC_result.cfy_clab), opt.clab_rereference, '*');
	A_reref(:,	sum(A_reref==0,1) == size(A_reref,1)) = []; %removes the zero-channels
end

%% --- Artifact rejection (trials and/or channels) based on variance criterion
%
flds= {'reject_artifacts', 'reject_channels', ...
       'reject_artifacts_opts', 'target_dist', 'clab'};
if data.isnew || ~isfield(data, 'previous_settings') || ...
      ~struct_areFieldsEqual(opt, data.previous_settings, flds),
  BC_result.rejected_trials= NaN;
  BC_result.rejected_clab= NaN;
  if opt.reject_artifacts | opt.reject_channels,
    if opt.create_figs,
      fig_state= fig_set(5, 'Hide',1, 'Name','Artifact rejection');
      data.figure_handles(end+1)= 5;
    end
    [mk_clean , rClab, rTrials]= ...
        reject_varEventsAndChannels(data.cnt, BC_result.mrk, ...
                                    opt.disp_ival, ...
                                    'CLab',BC_result.clab, ...
                                    'DoMultipass', 1, ...
                                    'Verbose', 1, ...
                                    'Visualize', opt.create_figs, ...
                                    opt.reject_artifacts_opts{:});
    if opt.create_figs, fig_publish(fig_state); end
    if opt.reject_artifacts,
      bbci_log_write(data, 'Rejected: %d trial(s).', length(rTrials));
      BC_result.mrk= mrk_selectEvents(BC_result.mrk, 'not',rTrials);
      BC_result.rejected_trials= rTrials;
    end
    if opt.reject_channels,
      bbci_log_write(data, 'Rejected channels: <%s>', str_vec2str(rClab));
      BC_result.rejected_clab= rClab;
      if iscell(BC_result.rejected_clab),   %% that means rejected_clab is not NaN
          cidx= find(ismember(BC_result.clab, BC_result.rejected_clab));
          BC_result.clab(cidx)= [];
      end
    end
  else
    % Avoid confusion with old figure from previous run
    fig_closeIfExists(5);
  end
else
  result_flds= {'rejected_trials', 'rejected_clab', 'clab'};
  BC_result= struct_copyFields(BC_result, previous, result_flds);
end

%% --- Segmentation and baselining ---
%
epo= proc_segmentation(data.cnt, BC_result.mrk, opt.disp_ival, ...
                       'clab', BC_result.clab);

%% Rejection of eyemovements based on max-min criterium
%
if opt.reject_eyemovements && opt.reject_eyemovements_crit.maxmin>0,
  [epo, iArte]= ...
      proc_rejectArtifactsMaxMin(epo, ...
                                 opt.reject_eyemovements_crit.maxmin, ...
                                 'clab', opt.reject_eyemovements_crit.clab, ...
                                 'ival', opt.reject_eyemovements_crit.ival, ...
                                 'verbose',1);
  BC_result.mrk= mrk_selectEvents(BC_result.mrk, 'not',iArte);
  BC_result.eyemovement_trials= iArte;
else
  BC_result.eyemovement_trials= NaN;
end

epo= proc_baseline(epo, opt.ref_ival, 'channelwise', 1);

%% --- Plot r^2 matrix and select intervals if requested ---
%
epo_r= proc_rSquareSigned(proc_selectChannels(epo,BC_result.cfy_clab));
%epo_r= proc_rSquareSigned(epo);
epo_r.className= {'sgn r^2 ( T , NT )'};  %% just make it shorter
if opt.create_figs,
  fig_state= fig_set(6, 'Hide',1, 'Name','r^2 Matrix');
end
if isempty(opt.cfy_ival) || isequal(opt.cfy_ival, 'auto'),
  [BC_result.cfy_ival, nfo]= ...
      procutil_selectTimeIntervals(epo_r, 'Visualize', opt.create_figs, ...
                            'VisuScalps', opt.create_figs, ...
                            'Sort', 1, ...
                            'IntersampleTiming', 1, ...
                            'IvalPickPeak', opt.cfy_ival_pick_peak, ...
                            'IvalMax', opt.cfy_maxival);
  bbci_log_write(data.log.fid, 'Selected time intervals:');
  bbci_log_write(data.log.fid, '  [%g %g] ms', BC_result.cfy_ival');
else
  BC_result.cfy_ival= opt.cfy_ival;
  if opt.create_figs,
    plot_scoreMatrix(epo_r, BC_result.cfy_ival);
  end
end
if opt.create_figs,
  fig_publish(fig_state);
end

%% --- Visualize ERPs ---
%
if opt.create_figs,
  fig_state= fig_set(1, 'Hide',1, 'Name','ERP - grid plot');
  H= grid_plot(epo, mnt, defopt_erps, 'ColorOrder',opt_scalp_erp.ColorOrder);
  %grid_markInterval(BC_result.cfy_ival);
  if isfield(H, 'scale'),
    grid_addBars(epo_r, 'HScale',H.scale);
  else
    grid_addBars(epo_r);
  end
  fig_publish(fig_state);

  fig_state= fig_set(2, 'Hide',1, 'Name','ERP - scalp maps');
  H= plot_scalpEvolutionPlusChannel(epo, mnt, opt.clab_erp, ...
                                    BC_result.cfy_ival, ...
                                    opt_scalp_erp);
  grid_addBars(epo_r);
  fig_publish(fig_state);

  if isempty(opt.clab_rsq) || isequal(opt.clab_rsq,'auto'),
    opt.clab_rsq= unique_unsort({nfo.peak_clab}, 4);
  end
  fig_state= fig_set(4, 'Hide',1, 'Name','ERP - r^2 scalp maps');
  plot_scalpEvolutionPlusChannel(epo_r, mnt, opt.clab_rsq, ...
                                 BC_result.cfy_ival, ...
                                 opt_scalp_r);
  fig_publish(fig_state);
end
clear epo*

%% --- Feature extraction ---
%
BC_result.ref_ival= opt.ref_ival;
BC_result.nSequences= opt.nSequences; %% Future extension: select by simulation

bbci.signal.clab= BC_result.cfy_clab;

cnt_processed = proc_selectChannels(data.cnt,bbci.signal.clab);

if ~isempty(opt.clab_rereference),
	cnt_processed = proc_linearDerivation(cnt_processed, A_reref);
	bbci.signal.proc=  {{@online_linearDerivation, A_reref}};
else
	bbci.signal.proc= {};
end

if opt.whitening
	C= cov(cnt_processed.x);
	[V,D]= eig(C);
	A_whitening= V*diag(1./sqrt(diag(D)))*V';
    cnt_processed= proc_linearDerivation(cnt_processed, A_whitening);
	bbci.signal.proc= cat(2, bbci.signal.proc, {{@online_linearDerivation, A_whitening}});
end

if ~isempty(opt.band)  %set filter if required!
  online_fs= data.cnt.fs;   %extract sampling rate from bbci-strukt
  [filt_b,filt_a]= cheby2(5, 20, opt.band/online_fs*2);
  bbci.signal.proc= cat(2, bbci.signal.proc, {{@online_filt, filt_b, filt_a}});
end

if isempty(opt.ref_ival),     % no baselining
  bbci.feature.proc= {{@proc_jumpingMeans, BC_result.cfy_ival}};
  min_t = min(BC_result.cfy_ival(:));
else
  bbci.feature.proc= {{@proc_baseline, BC_result.ref_ival, 'beginning_exact'}, ...
                      {@proc_jumpingMeans, BC_result.cfy_ival}};
  min_t =  min(BC_result.ref_ival(1), min(BC_result.cfy_ival(:)));
end
bbci.feature.ival= [min_t max(BC_result.cfy_ival(:))];

if opt.control_per_stimulus,
  bbci.control.fcn= @bbci_control_ERPSpellerBinary;
else
  bbci.control.fcn= @bbci_control_ERPSpeller;
end
bbci.control.param= {struct('nClasses',         opt.nClasses, ...
                            'nSequences',       BC_result.nSequences, ...
                            'mrk2feedback_fcn', opt.mrk2feedback_fcn)};
bbci.control.condition.marker= opt.cue_markers;
if ~isfield(bbci, 'quit_condition') || ~isfield(bbci.quit_condition, 'marker'),
  bbci.quit_condition.marker = 255;
end

fv= proc_segmentation(cnt_processed, BC_result.mrk, bbci.feature.ival);
fv= bbci_calibrate_evalFeature(fv, bbci.feature);

bbci.classifier.C= trainClassifier(fv, opt.model);

opt_xv= struct('SampleFcn',    {{@sample_chronKFold, 10}}, ...
               'LossFcn',      @loss_rocArea);

loss= crossvalidation(fv, opt.model, opt_xv);
bbci_log_write(data.log.fid, 'ROC-loss: %.1f%%', 100*loss);
%[loss, dum, outTe]= crossvalidation(fv, opt.model, opt_xv);
%me= val_confusionMatrix(fv, outTe, 'Mode','normalized');
%bbci_log_write(data.log.fid, ['ROC-loss: %.1f%%  |  ' ...
%                    'Correct Hits: %.1f%%,  Correct Miss: %.1f%%'], ...
%               100*loss, 100*me(1,1), 100*me(2,2));

data.feature= fv;
data.result= BC_result;
data.result.rocloss = 100*loss;
%data.result.me = me;

author: bbci


text:
function handle = progressbar( handle,increment,string,titlestr )
%
% progressbar - shows a progress bar dialog based on the function "waitbar"
%
% Format: handle = progressbar( handle,increment [,string,titlestr] )
%
% Input:    handle      - handle to current progress bar, [] for a new one
%           increment   - a fraction (0..1) to increment by.
%                         (-1) signals the function to remove the handle from
%                         the persistant list and close the progressbar
%           string      - a string to be replaced in the progress bar (optional)
%           titlestr    - a title for the dialog box (optional)
%
% Output:   handle      - a graphic handle of the dialog box
%
%
% NOTE:     this function uses a persistant list of handles and thier values.
%           therefore,  to delete a progressbar, please call the function with:
%               progressbar( handle,-1 );
%
%           an "abort" button is formed inside the progressbar, if the calling process
%           uses the persistent function "gui_active". when the "abort" button is pressed,
%           the Callback function "gui_active" changes it's value to zero, which enables
%           the application to interactively stop execution
%
% Example:  gui_active(1);      % will add an abort button
%           h           = progressbar( [],0,'my test' );
%           max_count   = 1e+3;
%           for idx = 1:max_count
%               fprintf( '%d\n',idx )';
%               h = progressbar( h,1/max_count );
%               if ~gui_active
%                   break;
%               end
%           end
%           progressbar( h,-1 );
%

persistent handle_list counter_list last_handle last_idx;

% initialize
% =============
call_flag = min( nargin,4 );

% analyze input and decide what to do
% ====================================
if isempty( handle )            % create a new dialog
    counter_list(end+1) = 0;
    last_idx            = length( counter_list );
    switch call_flag
    case 2, last_handle = waitbar( increment,'Please Wait...' );
    case 3, last_handle = waitbar( increment,string );
    case 4, last_handle = waitbar( increment,string,'Name',titlestr );
    end
    handle_list(end+1)  = last_handle;
    handle              = last_handle;
    check_position( handle_list );      % so that the figures don't hide each other
    if (gui_active)
        add_button( last_handle );      % add the abort button if the state of the gui_active is set
    end

elseif ( increment == -1 )      % delete correct handle from the list
    last_handle             = handle;
    last_idx                = find( handle_list == handle );
    handle_list( last_idx ) = [];
    counter_list( last_idx )= [];
    if ishandle( last_handle )      % close the figure, if it's open
        close( last_handle );       % since user can close it by him self
    end
    last_handle             = [];

elseif (handle == last_handle)  % update last dialog
    counter_list(last_idx)  = counter_list(last_idx) + increment;
    if ishandle( handle )       % nobody killed my figure
        switch call_flag
        case 2, waitbar( counter_list(last_idx),handle );
        case 3, waitbar( counter_list(last_idx),handle,string );
        case 4, waitbar( counter_list(last_idx),handle,string,'Name',titlestr );
        end
    else                        % somebody killed my figure -> so I create it again
        switch call_flag
        case 2, handle = waitbar( counter_list(last_idx),'Please Wait...' );
        case 3, handle = waitbar( counter_list(last_idx),string );
        case 4, handle = waitbar( counter_list(last_idx),string,'Name',titlestr );
        end
        handle_list(last_idx)   = handle;
        last_handle             = handle;
        check_position( handle_list );      % so that the figures don't hide each other
        if (gui_active)
            add_button( last_handle );      % add the abort button if the state of the gui_active is set
        end
    end
else                            % find the handle inside the list
    last_handle = handle;
    last_idx    = find( handle_list == handle );
    if ~isempty( last_idx )
        counter_list(last_idx)  = counter_list(last_idx) + increment;
        switch call_flag
        case 2, waitbar( counter_list(last_idx),last_handle );
        case 3, waitbar( counter_list(last_idx),last_handle,string );
        case 4, waitbar( counter_list(last_idx),last_handle,string,'Name',titlestr );
        end
    end
end

% update display after all
% ==========================
drawnow;

% =======================================================================================
%                               Inner Function Implementation
% =======================================================================================

function add_button( fig_handle )
%
% adds the abort button to the waitbar
%

% collect handles and set control units to pixels
axes_handle     = get( fig_handle,'currentaxes' );
last_fig_units  = get( fig_handle,'units' );
last_axes_units = get( axes_handle,'units' );
set( fig_handle,'units','pixels' );
set( axes_handle,'units','pixels' );

% read controls position
fig_position    = get( fig_handle,'position' );
axes_position   = get( axes_handle,'position' );
fig_width       = fig_position(3);
fig_height      = fig_position(4);
axes_xcoord     = axes_position(1);
axes_ycoord     = axes_position(2);
axes_width      = axes_position(3);
axes_height     = axes_position(4);

% load the button icon and create the button
load( 'gauge_abort_icon' );
button_width    = ButtonSize16x16(1)+2;
button_height   = ButtonSize16x16(2)+2;
button_margin   = 10;
button_xcoord   = (fig_width + axes_width + axes_xcoord - button_width)/2 - button_margin;
button_ycoord   = (axes_height - button_height)/2 + axes_ycoord;
button_handle   = uicontrol( 'Parent',fig_handle,'units','pixels',...
    'Position',[ button_xcoord,button_ycoord,button_width,button_height ],...
    'Callback','gui_active(0);progressbar(get(gcbo,''parent''),-1);close(get(gcbo,''parent''));',...
    'CData',Icon16x16 );

% resize axis to accommodate the button, and restore axes and figure units back to previous
axes_position(3) = axes_width - button_width - button_margin;
set( axes_handle,'position',axes_position );
set( fig_handle,'units',last_fig_units );
set( axes_handle,'units',last_axes_units );

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

function check_position( handle_list )
%
% makes sure that the progressbar does not hide it's nested progressbars
%

% process only if there is more than one progress bar on the screen
if (length(handle_list)>1)
    y_increment = 70;                   % pixels
    x_increment = 30;                   % pixels

    % change units to pixels
    screen_units    = get( 0,'units' );
    last_fig_units  = get( handle_list(end-1),'units' );
    cur_fig_units   = get( handle_list(end),'units' );
    set( 0,'units','pixels' );
    set( handle_list(end-1),'units','pixels' );
    set( handle_list(end),'units','pixels' );

    % get positions, and calc new position for progress bar
    screen_size     = get( 0,'screensize' );
    last_position   = get( handle_list(end-1),'position' );
    cur_position    = get( handle_list(end),'position' );
    cur_position(1) = last_position(1) + x_increment;
    cur_position(2) = last_position(2) - y_increment;

    % check that we don't go outside the screen
    if (cur_position(1)+cur_position(3)>screen_size(3))
        cur_position(1) = x_increment;
    end
    if (cur_position(2)<screen_size(1))
        cur_position(2) = screen_size(4) - y_increment - cur_position(4);
    end

    % store new position and restore units
    set( handle_list(end),'position',cur_position );
    set( 0,'units',screen_units );
    set( handle_list(end-1),'units',last_fig_units );
    set( handle_list(end),'units',cur_fig_units );
end


text:
% OZMOVIE  This program creates a movie of the Newton-Krylov solution of
% the Ornstein-Zernike equations.
%
function ozmovie
global L U rho
n=257;
epsilon=.1; sigma=2.0; rho=.2; beta=10; L=9;
dx=L/(n-1); r=0:dx:L; r=r';
%
% U is a global variable used in the function evaluation. See the text.
%
U=elj(r,sigma,epsilon,beta);
%
tol=[1.d-8,1.d-8];
x=zeros(2*n,1);
t1=cputime;
parms=[40,80,-.1];
[sol, it_hist, ierr,oz_hist] = nsoli(x,'oz',tol,parms);
%
% Build the movie of the iterations in h. Put the
% iteration counter on the screen.
%
[hr,hc]=size(oz_hist); n=hr/2;
for im=1:hc
    axis([0 9 -4 4]);
    h=oz_hist(1:n,im);
    L=9; dx=L/(n-1); r=0:dx:L; r=r';
    plot(r,h,'-'); axis([0 9 -2 3]); text(3,-1,num2str(im));
    xlabel('r'); ylabel('h');
    mov(im)=getframe;
end
%
% Showtime!
%
movie(mov,1,1);
%
%
%
function u=elj(r,sigma,epsilon,beta)
n2=length(r);
ra=r(2:n2);
r12=(sigma./ra).^12; r6=(sigma./ra).^6;
ua=exp(-4*beta*epsilon*(r12-r6));
u=[0,ua']';

repo: Aero-matlab


text:
function DAt = getDAtm(A,Ablkjc,dense,DAtdenq,d,K)
% DAt = getDAtm(A,Ablkjc,dense,DAtdenq,d,K)
%
% GETDATM  Computes d[k]'*Aj[k] for each lorentz block k and constraint j.
%
% ******************** INTERNAL FUNCTION OF SEDUMI ********************
%
% See also sedumi, getada2.

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%
nq = length(K.q);
DAt.q = extractA(A,Ablkjc,1,2,K.mainblks(1),K.mainblks(2));
if nq > 0
    DAt.q = spdiags(d.q1,0,nq,nq) * DAt.q;
    DAt.q = DAt.q + ddot(d.q2, A, K.qblkstart, Ablkjc);
end
DAt.denq = adendotd(dense,d,DAt.q(dense.q,:)',DAtdenq,K.qblkstart);
if ~isempty(dense.q)
    DAt.q(dense.q,:) = 0.0;
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
%% Test of different Fourier boundaries detection strategies

clear all

%% User setup
% Choose the signal you want to analyze
% (sig1,sig2,sig3,sig4=ECG,sig5=seismic,lena,textures)
signal = 'sig4';

% Choose the wanted preprocessing (none,plaw,poly,morpho)
params.preproc = 'none';
params.degree=5; % degree for the polynomial interpolation

% Choose the wanted detection method (locmax,locmaxmin,ftc)
params.method = 'locmaxmin';
params.N = 6; % maximum number of band for the locmaxmin method
params.completion = 0;

% Perform the detection on the log spectrum instead the spectrum
params.log=0;

%% Load signals
switch lower(signal)
    case 'sig1'
        load('sig1.mat');
        t=0:1/length(f):1-1/length(f);
    case 'sig2'
        load('sig2.mat');
        t=0:1/length(f):1-1/length(f);
    case 'sig3'
        load('sig3.mat');
        t=0:1/length(f):1-1/length(f);
    case 'sig4'
        load('sig4.mat');
        t=0:length(f)-1;
    case 'sig5'
        load('seismic.mat');
        f=f(10000:20000); %sub portion of the signal used in the paper
        t=0:length(f)-1;
    case 'lena'
        load lena
        l=round(size(f,2)/2);
        imR=[f(:,(l-1:-1:1)) f f(:,(end:-1:end-l+1))];
        fftim=fft(imR');
        ff=abs(sum(abs(fftim),2)/size(fftim,2));
    case 'textures'
        load('texture.mat');
        l=round(size(f,2)/2);
        imR=[f(:,(l-1:-1:1)) f f(:,(end:-1:end-l+1))];
        fftim=fft(imR');
        ff=abs(sum(abs(fftim),2)/size(fftim,2));
end

if (~strcmp(signal,'lena')) && (~strcmp(signal,'textures'))
    % We extend the signal by miroring to deal with the boundaries
    l=round(length(f)/2);
    f=[f(l-1:-1:1);f;f(end:-1:end-l+1)];

    % We compute the Fourier transform of f
    ff=abs(fft(f));
end

%% Perform the detection and plot the detected boundaries
boundaries = EWT_Boundaries_Detect(ff,params);
boundaries = boundaries*pi/round(length(ff)/2);
Show_EWT_Boundaries(abs(fft(f)),boundaries,10);


text:
% function dodSpline = hmrR_MotionCorrectSpline_Nirs(dod, t, SD, tInc, p)
%
% UI NAME:
% Spline_Motion_Correction
%
% dodSpline = hmrR_MotionCorrectSpline_Nirs(dod, t, SD, tInc, p)
% Perform a cubic spline correction of the motion artifacts identified in
% tInc. The algorithm follows the procedure describe by
% Scholkmann et al., Physiol. Meas. 31, 649-662 (2010). Set p = -1 to skip
% this function.
%
% INPUTS:
% dod:   delta_OD
% t:     time vector
% SD:    SD structure
% tInc:  Matrix of included time points (1=included; 0=not included (movement)
%        The matrix is #time points x #channels and usually comes from
%        hmrMotionArtifactByChannel()
% p:     Parameter p used in the spline interpolation. The value
%        recommended in the literature is 0.99. Use -1 if you want to skip this
%        motion correction.
%
% OUTPUTS:
% dodSpline:  dod after spline interpolation correction, same size as dod
%            (Channels that are not in the active ml remain unchanged)
%
% LOG:
% created 01-26-2012, J. Selb
%
% TO DO:
%

function dodSpline = hmrR_MotionCorrectSpline_Nirs(dod, t, SD, tInc, p, turnon)

if exist('turnon')
   if turnon==0
       dodSpline = dod;
   return;
   end
end

% if p outside its authorized range, set to 0.99
if p>1 || p<0
    display('Parameter has to be between 0 and 1. Returning with no correction');
    dodSpline = dod;
    return;
end

fs = 1/mean(t(2:end)-t(1:end-1));

% window widths limits for computing the mean in the segment shifts
dtShort = 0.3;  % seconds
dtLong  = 3;    % seconds

ml = SD.MeasList;
mlAct = SD.MeasListAct; % prune bad channels

lstAct = find(mlAct==1);
dodSpline = dod;
t = t(:);  % needs to be a column vector

for ii = 1:length(lstAct)

    idx_ch = lstAct(ii);

    lstMA = find(tInc(:,idx_ch)==0);   % sublist of motion artifact segments

    if ~isempty(lstMA)

        % Find indexes of starts and ends of MA segments
        lstMs = find(diff(tInc(:,idx_ch))==-1);   % starting indexes of mvt segments
        lstMf = find(diff(tInc(:,idx_ch))==1);    % ending indexes of mvt segments

        % Case where there's a single MA segment, that either starts at the
        % beginning or ends at the end of the total time duration
        if isempty(lstMf)
            lstMf = size(tInc,1);
        end
        if isempty(lstMs)
            lstMs = 1;
        end
        % If any MA segment either starts at the beginning or
        % ends at the end of the total time duration
        if lstMs(1)>lstMf(1)
            lstMs = [1;lstMs];
        end
        if lstMs(end)>lstMf(end)
            lstMf(end+1,1) = size(tInc,1);
        end

        lstMl = lstMf-lstMs;    % lengths of MA segments
        nbMA = length(lstMl);   % number of MA segments

        % Do the spline interpolation on each MA segment
        % only include channels in the active meas list

        for jj = 1:nbMA
            lst = lstMs(jj):(lstMf(jj)-1);
            % spline interp
            SplInterp = csaps(t(lst)', dod(lst,idx_ch)', p, t(lst)')';
            % corrected signal = original signal - spline interpolation
            dodSpline(lst,idx_ch) = dod(lst,idx_ch) - SplInterp;
        end

        % Reconstruction of the whole time series (shift each segment)

        %% First MA segment: shift to the previous noMA segment if it exists,
        % to the next noMA segment otherwise
        lst = (lstMs(1)):(lstMf(1)-1);
        SegCurrLength = lstMl(1);
        if SegCurrLength < dtShort*fs
            windCurr = SegCurrLength;
        elseif SegCurrLength < dtLong*fs
            windCurr = floor(dtShort*fs);
        else
            windCurr = floor(SegCurrLength/10);
        end

        if lstMs(1)>1
            SegPrevLength = length(1:(lstMs(1)-1));
            if SegPrevLength < dtShort*fs
                windPrev = SegPrevLength;
            elseif SegPrevLength < dtLong*fs
                windPrev = floor(dtShort*fs);
            else
                windPrev = floor(SegPrevLength/10);
            end
            meanPrev = mean(dodSpline(lst(1)-windPrev:(lst(1)-1), idx_ch));
            meanCurr = mean(dodSpline(lst(1):(lst(1)+windCurr-1), idx_ch));
            dodSpline(lst,idx_ch) = dodSpline(lst,idx_ch) - meanCurr + meanPrev;

        else
            if length(lstMs)>1
                SegNextLength = length(lstMf(1):(lstMs(2)));
            else
                SegNextLength = length(lstMf(1):size(tInc,1));
            end
            if SegNextLength < dtShort*fs
                windNext = SegNextLength;
            elseif SegNextLength < dtLong*fs
                windNext = floor(dtShort*fs);
            else
                windNext = floor(SegNextLength/10);
            end
            meanCurr = mean(dodSpline((lst(end)-windCurr):(lst(end)-1),  idx_ch));
            meanNext = mean(dodSpline((lst(end)+1):(lst(end)+windNext), idx_ch));
            dodSpline(lst,idx_ch) = dodSpline(lst,idx_ch) - meanCurr + meanNext;
        end

        %% Intermediate segments
        for kk=1:(nbMA-1)
            % no motion
            lst = lstMf(kk):(lstMs(kk+1)-1);
            SegPrevLength = lstMl(kk);
            SegCurrLength = length(lst);
            if SegPrevLength < dtShort*fs
                windPrev = SegPrevLength;
            elseif SegPrevLength < dtLong*fs
                windPrev = floor(dtShort*fs);
            else
                windPrev = floor(SegPrevLength/10);
            end
            if SegCurrLength < dtShort*fs
                windCurr = SegCurrLength;
            elseif SegCurrLength < dtLong*fs
                windCurr = floor(dtShort*fs);
            else
                windCurr = floor(SegCurrLength/10);
            end
            meanPrev = mean(dodSpline((lst(1)-windPrev):(lst(1)-1), idx_ch));
            meanCurr = mean(dod(lst(1):(lst(1)+windCurr-1), idx_ch));

            dodSpline(lst,idx_ch) = dod(lst,idx_ch) - meanCurr + meanPrev;

            % motion
            lst = (lstMs(kk+1)):(lstMf(kk+1)-1);
            SegPrevLength = SegCurrLength;
            SegCurrLength = lstMl(kk+1);
            if SegPrevLength < dtShort*fs
                windPrev = SegPrevLength;
            elseif SegPrevLength < dtLong*fs
                windPrev = floor(dtShort*fs);
            else
                windPrev = floor(SegPrevLength/10);
            end
            if SegCurrLength < dtShort*fs
                windCurr = SegCurrLength;
            elseif SegCurrLength < dtLong*fs
                windCurr = floor(dtShort*fs);
            else
                windCurr = floor(SegCurrLength/10);
            end
            meanPrev = mean(dodSpline((lst(1)-windPrev):(lst(1)-1), idx_ch));
            meanCurr = mean(dodSpline(lst(1):(lst(1)+windCurr-1), idx_ch));

            dodSpline(lst,idx_ch) = dodSpline(lst,idx_ch) - meanCurr + meanPrev;
        end

        %% Last not MA segment
        if lstMf(end)<length(t)
            lst = (lstMf(end)-1):length(t);
            SegPrevLength = lstMl(end);
            SegCurrLength = length(lst);
            if SegPrevLength < dtShort*fs
                windPrev = SegPrevLength;
            elseif SegPrevLength < dtLong*fs
                windPrev = floor(dtShort*fs);
            else
                windPrev = floor(SegPrevLength/10);
            end
            if SegCurrLength < dtShort*fs
                windCurr = SegCurrLength;
            elseif SegCurrLength < dtLong*fs
                windCurr = floor(dtShort*fs);
            else
                windCurr = floor(SegCurrLength/10);
            end
            meanPrev = mean(dodSpline((lst(1)-windPrev):(lst(1)-1), idx_ch));
            meanCurr = mean(dod(lst(1):(lst(1)+windCurr-1), idx_ch));

            dodSpline(lst,idx_ch) = dod(lst,idx_ch) - meanCurr + meanPrev;
        end

    %else
     %   dodSpline(:,i_ch) = dod;
    end
end

author: BUNPC


text:
%  internal function

%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)

function [hdr, filetype, fileprefix, machine] = load_nii_hdr(fileprefix)

   if ~exist('fileprefix','var'),
      error('Usage: [hdr, filetype, fileprefix, machine] = load_nii_hdr(filename)');
   end

   machine = 'ieee-le';
   new_ext = 0;

   if findstr('.nii',fileprefix) & strcmp(fileprefix(end-3:end), '.nii')
      new_ext = 1;
      fileprefix(end-3:end)='';
   end

   if findstr('.hdr',fileprefix) & strcmp(fileprefix(end-3:end), '.hdr')
      fileprefix(end-3:end)='';
   end

   if findstr('.img',fileprefix) & strcmp(fileprefix(end-3:end), '.img')
      fileprefix(end-3:end)='';
   end

   if new_ext
      fn = sprintf('%s.nii',fileprefix);

      if ~exist(fn)
         msg = sprintf('Cannot find file "%s.nii".', fileprefix);
         error(msg);
      end
   else
      fn = sprintf('%s.hdr',fileprefix);

      if ~exist(fn)
         msg = sprintf('Cannot find file "%s.hdr".', fileprefix);
         error(msg);
      end
   end

   fid = fopen(fn,'r',machine);

   if fid < 0,
      msg = sprintf('Cannot open file %s.',fn);
      error(msg);
   else
      fseek(fid,0,'bof');

      if fread(fid,1,'int32') == 348
         hdr = read_header(fid);
         fclose(fid);
      else
         fclose(fid);

         %  first try reading the opposite endian to 'machine'
         %
         switch machine,
         case 'ieee-le', machine = 'ieee-be';
         case 'ieee-be', machine = 'ieee-le';
         end

         fid = fopen(fn,'r',machine);

         if fid < 0,
            msg = sprintf('Cannot open file %s.',fn);
            error(msg);
         else
            fseek(fid,0,'bof');

            if fread(fid,1,'int32') ~= 348

               %  Now throw an error
               %
               msg = sprintf('File "%s" is corrupted.',fn);
               error(msg);
            end

            hdr = read_header(fid);
            fclose(fid);
         end
      end
   end

   if strcmp(hdr.hist.magic, 'n+1')
      filetype = 2;
   elseif strcmp(hdr.hist.magic, 'ni1')
      filetype = 1;
   else
      filetype = 0;
   end

   return					% load_nii_hdr

%---------------------------------------------------------------------
function [ dsr ] = read_header(fid)

        %  Original header structures
	%  struct dsr
	%       {
	%       struct header_key hk;            /*   0 +  40       */
	%       struct image_dimension dime;     /*  40 + 108       */
	%       struct data_history hist;        /* 148 + 200       */
	%       };                               /* total= 348 bytes*/

    dsr.hk   = header_key(fid);
    dsr.dime = image_dimension(fid);
    dsr.hist = data_history(fid);

    %  For Analyze data format
    %
    if ~strcmp(dsr.hist.magic, 'n+1') & ~strcmp(dsr.hist.magic, 'ni1')
        dsr.hist.qform_code = 0;
        dsr.hist.sform_code = 0;
    end

    return					% read_header

%---------------------------------------------------------------------
function [ hk ] = header_key(fid)

    fseek(fid,0,'bof');

	%  Original header structures
	%  struct header_key                     /* header key      */
	%       {                                /* off + size      */
	%       int sizeof_hdr                   /*  0 +  4         */
	%       char data_type[10];              /*  4 + 10         */
	%       char db_name[18];                /* 14 + 18         */
	%       int extents;                     /* 32 +  4         */
	%       short int session_error;         /* 36 +  2         */
	%       char regular;                    /* 38 +  1         */
	%       char dim_info;   % char hkey_un0;        /* 39 +  1 */
	%       };                               /* total=40 bytes  */
	%
	% int sizeof_header   Should be 348.
	% char regular        Must be 'r' to indicate that all images and
	%                     volumes are the same size.

    v6 = version;
    if str2num(v6(1))<6
       directchar = '*char';
    else
       directchar = 'uchar=>char';
    end

    hk.sizeof_hdr    = fread(fid, 1,'int32')';	% should be 348!
    hk.data_type     = deblank(fread(fid,10,directchar)');
    hk.db_name       = deblank(fread(fid,18,directchar)');
    hk.extents       = fread(fid, 1,'int32')';
    hk.session_error = fread(fid, 1,'int16')';
    hk.regular       = fread(fid, 1,directchar)';
    hk.dim_info      = fread(fid, 1,'uchar')';

    return					% header_key

%---------------------------------------------------------------------
function [ dime ] = image_dimension(fid)

	%  Original header structures
	%  struct image_dimension
	%       {                                /* off + size      */
	%       short int dim[8];                /* 0 + 16          */
        %       /*
        %           dim[0]      Number of dimensions in database; usually 4.
        %           dim[1]      Image X dimension;  number of *pixels* in an image row.
        %           dim[2]      Image Y dimension;  number of *pixel rows* in slice.
        %           dim[3]      Volume Z dimension; number of *slices* in a volume.
        %           dim[4]      Time points; number of volumes in database
        %       */
	%       float intent_p1;   % char vox_units[4];   /* 16 + 4       */
	%       float intent_p2;   % char cal_units[8];   /* 20 + 4       */
	%       float intent_p3;   % char cal_units[8];   /* 24 + 4       */
	%       short int intent_code;   % short int unused1;   /* 28 + 2 */
	%       short int datatype;              /* 30 + 2          */
	%       short int bitpix;                /* 32 + 2          */
	%       short int slice_start;   % short int dim_un0;   /* 34 + 2 */
	%       float pixdim[8];                 /* 36 + 32         */
	%	/*
	%		pixdim[] specifies the voxel dimensions:
	%		pixdim[1] - voxel width, mm
	%		pixdim[2] - voxel height, mm
	%		pixdim[3] - slice thickness, mm
	%		pixdim[4] - volume timing, in msec
	%					..etc
	%	*/
	%       float vox_offset;                /* 68 + 4          */
	%       float scl_slope;   % float roi_scale;     /* 72 + 4 */
	%       float scl_inter;   % float funused1;      /* 76 + 4 */
	%       short slice_end;   % float funused2;      /* 80 + 2 */
	%       char slice_code;   % float funused2;      /* 82 + 1 */
	%       char xyzt_units;   % float funused2;      /* 83 + 1 */
	%       float cal_max;                   /* 84 + 4          */
	%       float cal_min;                   /* 88 + 4          */
	%       float slice_duration;   % int compressed; /* 92 + 4 */
	%       float toffset;   % int verified;          /* 96 + 4 */
	%       int glmax;                       /* 100 + 4         */
	%       int glmin;                       /* 104 + 4         */
	%       };                               /* total=108 bytes */

    dime.dim        = fread(fid,8,'int16')';
    dime.intent_p1  = fread(fid,1,'float32')';
    dime.intent_p2  = fread(fid,1,'float32')';
    dime.intent_p3  = fread(fid,1,'float32')';
    dime.intent_code = fread(fid,1,'int16')';
    dime.datatype   = fread(fid,1,'int16')';
    dime.bitpix     = fread(fid,1,'int16')';
    dime.slice_start = fread(fid,1,'int16')';
    dime.pixdim     = fread(fid,8,'float32')';
    dime.vox_offset = fread(fid,1,'float32')';
    dime.scl_slope  = fread(fid,1,'float32')';
    dime.scl_inter  = fread(fid,1,'float32')';
    dime.slice_end  = fread(fid,1,'int16')';
    dime.slice_code = fread(fid,1,'uchar')';
    dime.xyzt_units = fread(fid,1,'uchar')';
    dime.cal_max    = fread(fid,1,'float32')';
    dime.cal_min    = fread(fid,1,'float32')';
    dime.slice_duration = fread(fid,1,'float32')';
    dime.toffset    = fread(fid,1,'float32')';
    dime.glmax      = fread(fid,1,'int32')';
    dime.glmin      = fread(fid,1,'int32')';

    return					% image_dimension

%---------------------------------------------------------------------
function [ hist ] = data_history(fid)

	%  Original header structures
	%  struct data_history
	%       {                                /* off + size      */
	%       char descrip[80];                /* 0 + 80          */
	%       char aux_file[24];               /* 80 + 24         */
	%       short int qform_code;            /* 104 + 2         */
	%       short int sform_code;            /* 106 + 2         */
	%       float quatern_b;                 /* 108 + 4         */
	%       float quatern_c;                 /* 112 + 4         */
	%       float quatern_d;                 /* 116 + 4         */
	%       float qoffset_x;                 /* 120 + 4         */
	%       float qoffset_y;                 /* 124 + 4         */
	%       float qoffset_z;                 /* 128 + 4         */
	%       float srow_x[4];                 /* 132 + 16        */
	%       float srow_y[4];                 /* 148 + 16        */
	%       float srow_z[4];                 /* 164 + 16        */
	%       char intent_name[16];            /* 180 + 16        */
	%       char magic[4];   % int smin;     /* 196 + 4         */
	%       };                               /* total=200 bytes */

    v6 = version;
    if str2num(v6(1))<6
       directchar = '*char';
    else
       directchar = 'uchar=>char';
    end

    hist.descrip     = deblank(fread(fid,80,directchar)');
    hist.aux_file    = deblank(fread(fid,24,directchar)');
    hist.qform_code  = fread(fid,1,'int16')';
    hist.sform_code  = fread(fid,1,'int16')';
    hist.quatern_b   = fread(fid,1,'float32')';
    hist.quatern_c   = fread(fid,1,'float32')';
    hist.quatern_d   = fread(fid,1,'float32')';
    hist.qoffset_x   = fread(fid,1,'float32')';
    hist.qoffset_y   = fread(fid,1,'float32')';
    hist.qoffset_z   = fread(fid,1,'float32')';
    hist.srow_x      = fread(fid,4,'float32')';
    hist.srow_y      = fread(fid,4,'float32')';
    hist.srow_z      = fread(fid,4,'float32')';
    hist.intent_name = deblank(fread(fid,16,directchar)');
    hist.magic       = deblank(fread(fid,4,directchar)');

    fseek(fid,253,'bof');
    hist.originator  = fread(fid, 5,'int16')';

    return					% data_history

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [net, info] = wsddn_train(varargin)
% @author: Hakan Bilen
% wsddn_train: training script for WSDDN

opts.dataDir = fullfile(vl_rootnn, 'data') ;
opts.expDir = fullfile(vl_rootnn, 'exp') ;
opts.imdbPath = fullfile(vl_rootnn, 'data', 'imdbs', 'imdb-eb.mat');
opts.modelPath = fullfile(vl_rootnn, 'models', 'imagenet-vgg-f.mat') ;
opts.proposalType = 'eb' ;
opts.proposalDir = fullfile(vl_rootnn, 'data', 'EdgeBoxes') ;

opts.addBiasSamples = 1; % add Box Scores
opts.addLossSmooth  = 1; % add Spatial Regulariser
opts.softmaxTempCls = 1; % softmax temp for cls
opts.softmaxTempDet = 2; % softmax temp for det
opts.maxScale = 2000 ;

% if you have limited gpu memory (<6gb), you can change the next 2 params
opts.maxNumProposals = inf; % limit number (eg 1500)
opts.imageScales = [480,576,688,864,1200]; % scales
opts.minBoxSize = 20; % minimum bounding box size
opts.train.gpus = [] ;
opts.train.continue = true ;
opts.train.prefetch = true ;
opts.train.learningRate = 1e-5 * [ones(1,10) 0.1*ones(1,10)] ;
opts.train.weightDecay = 0.0005;
opts.train.numEpochs = 20;
opts.train.derOutputs = {'objective', 1} ;

opts.numFetchThreads = 1 ;
opts = vl_argparse(opts, varargin) ;

display(opts);

opts.train.batchSize = 1 ;
opts.train.expDir = opts.expDir ;
opts.train.numEpochs = numel(opts.train.learningRate) ;
%% -------------------------------------------------------------------------
%                                                   Database initialization
% -------------------------------------------------------------------------
fprintf('loading imdb...');
if exist(opts.imdbPath,'file')==2
  imdb = load(opts.imdbPath) ;
else
  if strcmp(opts.proposalType,'ssw')
    imdb = setup_voc07_ssw('dataDir',opts.dataDir, ...
      'proposalDir',opts.proposalDir,'loadTest',1);
  elseif strcmp(opts.proposalType,'eb')
    imdb = setup_voc07_eb('dataDir',opts.dataDir, ...
      'proposalDir',opts.proposalDir,'loadTest',1);
  else
    error('undefined proposal type %s\n',opts.proposalType)
  end

  imdbFolder = fileparts(opts.imdbPath);

  if ~exist(imdbFolder,'dir')
    mkdir(imdbFolder);
  end
  save(opts.imdbPath,'-struct', 'imdb', '-v7.3');
end

fprintf('done\n');

imdb = fixBBoxes(imdb, opts.minBoxSize, opts.maxNumProposals);

% use train + val for training
imdb.images.set(imdb.images.set == 2) = 1;
trainIdx = find(imdb.images.set == 1);

%% Compute image statistics (mean, RGB covariances, etc.)
imageStatsPath = fullfile(opts.dataDir, 'imageStats.mat') ;
if exist(imageStatsPath,'file')
  load(imageStatsPath, 'averageImage', 'rgbMean', 'rgbCovariance') ;
else

  images = imdb.images.name(imdb.images.set == 1) ;
  images = strcat([imdb.imageDir filesep],images) ;

  [averageImage, rgbMean, rgbCovariance] = getImageStats(images, ...
    'imageSize', [256 256], ...
    'numThreads', opts.numFetchThreads, ...
    'gpus', opts.train.gpus) ;
  save(imageStatsPath, 'averageImage', 'rgbMean', 'rgbCovariance') ;
end
[v,d] = eig(rgbCovariance) ;
rgbDeviation = v*sqrt(d) ;
clear v d ;

%% ------------------------------------------------------------------------
%                                                    Network initialization
% -------------------------------------------------------------------------
nopts.addBiasSamples = opts.addBiasSamples; % add Box Scores (only with Edge Boxes)
nopts.addLossSmooth  = opts.addLossSmooth; % add Spatial Regulariser
nopts.softmaxTempCls = opts.softmaxTempCls; % softmax temp for cls
nopts.softmaxTempDet = opts.softmaxTempDet; % softmax temp for det

nopts.averageImage = reshape(rgbMean,[1 1 3]) ;
% nopts.rgbVariance = 0.1 * rgbDeviation ;
nopts.rgbVariance = [] ;
nopts.numClasses = numel(imdb.classes.name) ;
nopts.classNames = imdb.classes.name ;

if ~exist(opts.modelPath,'file')
  [pname,fname,ext]  = fileparts(opts.modelPath) ;
  if ~exist(pname,'dir')
    mkdir(pname) ;
  end
  fprintf('Downloading %s to %s\n', [fname ext], pname) ;
  urlwrite(sprintf('http://www.vlfeat.org/matconvnet/models/%s',[fname ext]),...
    opts.modelPath) ;
end

net = load(opts.modelPath);
net = wsddn_init(net,nopts);

if nopts.addLossSmooth
  opts.train.derOutputs = {'objective', 1, 'lossTopB', 1e-4} ;
end

if ~exist(opts.expDir,'dir')
  mkdir(opts.expDir) ;
end

%% -------------------------------------------------------------------------
%                                                   Database stats
% -------------------------------------------------------------------------
bopts = net.meta.normalization;
net.meta.augmentation.jitterBrightness = 0 ;
% bopts.interpolation = 'bilinear';
bopts.jitterBrightness = net.meta.augmentation.jitterBrightness ;
bopts.imageScales = opts.imageScales;
bopts.numThreads = opts.numFetchThreads;
bopts.addLossSmooth = opts.addLossSmooth;
bopts.addBiasSamples = opts.addBiasSamples;
bopts.maxScale = opts.maxScale ;
bopts.vgg16 = any(arrayfun(@(a) strcmp(a.name, 'relu5_1'), net.layers)==1) ;
%% -------------------------------------------------------------------
%                                                                Train
% --------------------------------------------------------------------
% avoid test data
valIdx = find(imdb.images.set == 3);
valIdx = valIdx(1:5:end) ;
% valIdx = [];

%%
bopts.useGpu = numel(opts.train.gpus) >  0 ;
bopts.prefetch = opts.train.prefetch;

info = cnn_train_dag(net, imdb, @(i,b) ...
  getBatch(bopts,i,b), ...
  opts.train, 'train', trainIdx, ...
  'val', valIdx) ;

%% -------------------------------------------------------------------
%                                                       Deploy network
% --------------------------------------------------------------------
if ~exist(fullfile(opts.expDir,'net.mat'),'file')
  removeLoss = {'dagnn.Loss','dagnn.DropOut'};
  for i=1:numel(removeLoss)
    dagRemoveLayersOfType(net,removeLoss{i}) ;
  end

  net.mode = 'test' ;
  net_ = net ;
  net = net_.saveobj() ;
  save(fullfile(opts.expDir,'net.mat'), '-struct','net');
end
% --------------------------------------------------------------------
function inputs = getBatch(opts, imdb, batch)
% --------------------------------------------------------------------
if isempty(batch)
  inputs = {'input', [], 'label', [], 'rois', [], 'ids', []};
  return;
end

opts.scale = opts.imageScales(randi(numel(opts.imageScales)));
opts.flip = randi(2,numel(batch),1)-1; % random flip
is_vgg16 = opts.vgg16 ;
opts = rmfield(opts,'vgg16') ;

images = strcat([imdb.imageDir filesep], imdb.images.name(batch)) ;
opts.prefetch = (nargout == 0);

[im,rois] = wsddn_get_batch(images, imdb, batch, opts);

if nargout>0
  rois = single(rois') ;
  labels = imdb.images.label(:,batch) ;
  labels = reshape(labels,[1 1 size(labels,1) numel(batch)]);

  if opts.useGpu > 0
    im = gpuArray(im) ;
    rois = gpuArray(rois) ;
  end

  if ~isempty(rois)
   rois = rois([1 3 2 5 4],:) ;
  end

  ss = [16 16] ;

  if is_vgg16
    o0 = 8.5 ;
    o1 = 9.5 ;
  else
    o0 = 18 ;
    o1 = 9.5 ;
  end

  rois = [ rois(1,:); ...
    floor((rois(2,:) - o0 + o1) / ss(1) + 0.5) + 1;
    floor((rois(3,:) - o0 + o1) / ss(2) + 0.5) + 1;
    ceil((rois(4,:) - o0 - o1) / ss(1) - 0.5) + 1;
    ceil((rois(5,:) - o0 - o1) / ss(2) - 0.5) + 1];

  inputs = {'input', im, 'label', labels, 'rois', rois, 'ids', batch} ;

  if opts.addLossSmooth
    inputs{end+1} = 'boxes' ;
    inputs{end+1} = imdb.images.boxes{batch} ;
  end

  if opts.addBiasSamples==1
    boxScore = reshape(imdb.images.boxScores{batch},[1 1 1 numel(imdb.images.boxScores{batch})]);
    inputs{end+1} = 'boxScore';
    inputs{end+1} = boxScore ;
  end
end

% -------------------------------------------------------------------------
function imdb = fixBBoxes(imdb, minSize, maxNum)
% -------------------------------------------------------------------------
for i=1:numel(imdb.images.name)
  bbox = imdb.images.boxes{i};
  % remove small bbox
  isGood = (bbox(:,3)>=bbox(:,1)+minSize) & (bbox(:,4)>=bbox(:,2)+minSize);
  bbox = bbox(isGood,:);
  % remove duplicate ones
  [dummy, uniqueIdx] = unique(bbox, 'rows', 'first');
  uniqueIdx = sort(uniqueIdx);
  bbox = bbox(uniqueIdx,:);
  % limit number for training
  if imdb.images.set(i)~=3
    nB = min(size(bbox,1),maxNum);
  else
    nB = size(bbox,1);
  end

  if isfield(imdb.images,'boxScores')
    imdb.images.boxScores{i} = imdb.images.boxScores{i}(isGood);
    imdb.images.boxScores{i} = imdb.images.boxScores{i}(uniqueIdx);
    imdb.images.boxScores{i} = imdb.images.boxScores{i}(1:nB);
  end
  imdb.images.boxes{i} = bbox(1:nB,:);
  %   [h,w,~] = size(imdb.images.data{i});
  %   imdb.images.boxes{i} = [1 1 h w];

end

% -------------------------------------------------------------------------
function layers = dagFindLayersOfType(net, type)
% -------------------------------------------------------------------------
layers = [] ;
for l = 1:numel(net.layers)
  if isa(net.layers(l).block, type)
    layers{1,end+1} = net.layers(l).name ;
  end
end
% -------------------------------------------------------------------------
function dagRemoveLayersOfType(net, type)
% -------------------------------------------------------------------------
names = dagFindLayersOfType(net, type) ;
for i = 1:numel(names)
  layer = net.layers(net.getLayerIndex(names{i})) ;
  net.removeLayer(names{i}) ;
  net.renameVar(layer.outputs{1}, layer.inputs{1}, 'quiet', true) ;
end

author: hbilen

repo: WSDDN

save_path: github-repos/MATLAB/hbilen-WSDDN


text:
function removeCDReactions(fileName,rxnRemoveList,printLevel)
%removes a list of reactions from a cell designer map, also removes
%correspinding species and species aliases if necessary
%
% INPUT
%   xmlStruct:      Structure obtained from the "xml2struct" function.
%                   To be kept for the conversion back to an XML file
%                   of the structure.
%
%   map:            Matlab structure of the map containing all the
%                   relevant fields usable for checking and correction.
%
%   rxnRemoveList:  Cell array of reaction abbreviation to be removed
%
%   printLevel:     {0,(1)}
%
% OUTPUT
%   xmlStructOut:   Structure for the conversion back to an XML file
%                   of the structure.
%
%   mapOut:         Matlab structure of the smaller map containing all the
%                   relevant fields usable for checking and correction.
%
% rxnNotInMap: boolean vector the length of rxnRemoveList
%                        indicating reactions that could not be found in the map
%
% Ronan Fleming 2020

if ~exist('printLevel','var')
    printLevel = 1;
end

if 0
    xmlStruct = xml2struct(fileName);

    nMapReactions = length(xmlStruct.sbml.model.listOfReactions.reaction);

    mapRxns = cell(nMapReactions,1);
    for i=1:nMapReactions
        mapRxns{i}= xmlStruct.sbml.model.listOfReactions.reaction{i}.annotation.celldesigner_colon_extension.celldesigner_colon_name.Text;
    end

else
    %use the map structure to be able to delete isolated species and
    %species aliases
    [xmlStruct, map] = transformXML2Map(fileName);
    mapRxns = map.rxnName;
end

[mapRxnsToRemove,LOCB]=ismember(mapRxns,rxnRemoveList);

rxnNotInMap=~ismember((1:length(rxnRemoveList))',LOCB);
if any(rxnNotInMap)
    if printLevel>0
        ind = find(rxnNotInMap);
        for i=1:length(ind)
            fprintf('%s%s\n',rxnRemoveList{ind(i)}, ' not present in the map')
        end
    end
end

[s,~]= size(map.sID);
[m,~]=size(map.sAlias);

%use the map mapping matrices to determine which species and species
%aliases to remove
% map.sID	s x r	logical	Logical matrix with rows = speciesID and columns = reactionsID
mapSpeciesToRemove = getCorrespondingRows(map.sID, true(s,1), mapRxnsToRemove, 'exclusive');

% map.sAlias	m x r	logical	Logical matrix with rows = speciesAlias and columns = reactionsID
mapMolsToRemove = getCorrespondingRows(map.sAlias, true(m,1), mapRxnsToRemove, 'exclusive');

%remove reactions
xmlStruct.sbml.model.listOfReactions.reaction = xmlStruct.sbml.model.listOfReactions.reaction(~mapRxnsToRemove);

%remove species
xmlStruct.sbml.model.listOfSpecies.species = xmlStruct.sbml.model.listOfSpecies.species(~mapSpeciesToRemove);

%remove species alias
xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias =...
    xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias(~mapMolsToRemove);

%write out the reduced xml file
struct2xml(xmlStruct, [fileName(1:end-4) '_subset.xml']);

author: opencobra

repo: cobratoolbox


text:
function unimplemented_error(varargin)
    % provide error dialog (if interactive) and error message for unimplemented features
    %
    % see also under_construction

    beep;
    if ZmapDialog.Data.Interactive
        errordlg('Under Construction','Not reimplemented yet');
    end
    error('ZMAP:UNIMPLEMENTED','Under Construction : Not yet implemented / reimplemented');
end

author: CelsoReyes


text:
function Kinematics = kinematics(States)
% function Kinematics = kinematics(States)
%
% Computer Generated File -- DO NOT EDIT
%
% This function was created by the function Write_Power()
% 13-Dec-2013 11:41:00
%
%  ARGUMENTS:
%   States = [Nstate x Ntime] matrix of states
%
%  RETURNS:
%   Kinematics = a struct with fields:
%      th1 = angle of leg 1
%      th2 = angle of leg 2
%      L1 = length of leg 1
%      L2 = length of leg 2
%      dth1 = (d/dt) angle of leg 1
%      dth2 = (d/dt) angle of leg 2
%      dL1 = (d/dt) length of leg 1
%      dL2 = (d/dt) length of leg 2
%
%
% Matthew Kelly
% Cornell University
%

% See also DERIVE_EOM
%
x0 = States(:,1); % (m) Hip horizontal position
y0 = States(:,2); % (m) Hip vertical position
x1 = States(:,3); % (m) Foot One horizontal position
y1 = States(:,4); % (m) Foot One vertical position
x2 = States(:,5); % (m) Foot Two horizontal position
y2 = States(:,6); % (m) Foot Two vertical position
dx0 = States(:,7); % (m/s) Hip horizontal velocity
dy0 = States(:,8); % (m/s) Hip vertical velocity
dx1 = States(:,9); % (m/s) Foot One horizontal velocity
dy1 = States(:,10); % (m/s) Foot One vertical velocity
dx2 = States(:,11); % (m/s) Foot Two horizontal velocity
dy2 = States(:,12); % (m/s) Foot Two vertical velocity

% Commonly used expressions
L1 = ((x0 - x1).^2 + (y0 - y1).^2).^(1./2);
L2 = ((x0 - x2).^2 + (y0 - y2).^2).^(1./2);
% dL1 = ((dx0 - dx1)*(x0 - x1) + (dy0 - dy1)*(y0 - y1))/L1;
% dL2 = ((dx0 - dx2)*(x0 - x2) + (dy0 - dy2)*(y0 - y2))/L2;
% th1 = atan2(x1 - x0, y1 - y0);
% th2 = atan2(x2 - x0, y2 - y0);
% dth1 = -((dy0 - dy1)*(x0 - x1) - (dx0 - dx1)*(y0 - y1))/L1^2;
% dth2 = -((dy0 - dy2)*(x0 - x2) - (dx0 - dx2)*(y0 - y2))/L2^2;

Kinematics.L1 = L1;
Kinematics.L2 = L2;
Kinematics.dL1 = ((dx0 - dx1).*(x0 - x1) + (dy0 - dy1).*(y0 - y1))./L1;
Kinematics.dL2 = ((dx0 - dx2).*(x0 - x2) + (dy0 - dy2).*(y0 - y2))./L2;
Kinematics.th1 = atan2(x1 - x0, y1 - y0);
Kinematics.th2 = atan2(x2 - x0, y2 - y0);
Kinematics.dth1 = -((dy0 - dy1).*(x0 - x1) - (dx0 - dx1).*(y0 - y1))./L1.^2;
Kinematics.dth2 = -((dy0 - dy2).*(x0 - x2) - (dx0 - dx2).*(y0 - y2))./L2.^2;

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [GAmp,GTime]=GxFSE(p)

global VCtl;
global VObj;

t1Start=p.t1Start;
t2Middle=p.t2Middle;
Gx1Sign=p.Gx1Sign;
Gx2Sign=p.Gx2Sign;

GAmp=[];
GTime=[];
GxAmp=(1/VCtl.FOVFreq)/((VObj.Gyro/(2*pi))*(1/VCtl.BandWidth));
tHalf=1/(2*(VObj.Gyro/(2*pi))*GxAmp*VCtl.RFreq);
tRamp=GxAmp/VCtl.MaxSlewRate;

% prephasing
[GAmp1,GTime1]=StdTrap(t1Start-tRamp,            ...
                       t1Start+tHalf+tRamp,      ...
                       t1Start,                          ...
                       t1Start+tHalf,                    ...
                       GxAmp*Gx1Sign,2,2,2);

GAmp=[GAmp GAmp1];
GTime=[GTime GTime1];

% frequency encoding
TimeOffset = (t2Middle+VCtl.TEAnchorTime)-floor(VCtl.FSE_ETL/2)*VCtl.FSE_ESP;
for i = 1: VCtl.FSE_ETL
    [GAmpt,GTimet]=StdTrap(TimeOffset + (i-1)*VCtl.FSE_ESP -tHalf -tRamp, ...
                           TimeOffset + (i-1)*VCtl.FSE_ESP +tHalf +tRamp, ...
                           TimeOffset + (i-1)*VCtl.FSE_ESP -tHalf, ...
                           TimeOffset + (i-1)*VCtl.FSE_ESP +tHalf, ...
                           GxAmp*Gx2Sign,2,2,2);
     GAmp=[GAmp GAmpt];
     GTime=[GTime GTimet];

end

[GTime,m,n]=unique(GTime);
GAmp=GAmp(m);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
% This is the demo script for ASRCF

clear;
clc;
close all;
setup_path();

% get seq information
base_path  = './seq';
video_path = [base_path '/' 'FaceOcc1'];
[seq, ground_truth] = load_video_info(video_path,'FaceOcc1');
seq.startFrame = 1;
seq.endFrame = seq.len;
seq.ground_truth=ground_truth;

% Run ASRCF- main function
results = run_ASRCF(seq, video_path);

author: Daikenan

repo: ASRCF

save_path: github-repos/MATLAB/Daikenan-ASRCF


text:
function p = new_parpool(number)

    if ~exist('number', 'var')
        number = cpu_cores();
    end

    if ~isempty(gcp('nocreate'))
        delete(gcp);
    end
    p = parpool(number);
end

repo: R-FCN


text:
function my_call(value)
fprintf('%d \n', value); % might check here whether this is a label at all
                         % since we get here by clicking on ANY text in figure

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [Center,R] = adaptiveDivision(PopObj,K)
% Reference point adaption

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    [N,M] = size(PopObj);

    if K == 1
        Center = zeros(1,M);
        R      = inf;
        subNum = [1,N];
    else
       %% Detect the number of subregion
        % Calculate the distance between each solution
        fmin     = min(PopObj,[],1);
        fmax     = max(PopObj,[],1);
        PopObj   = (PopObj-repmat(fmin,N,1))./repmat(fmax-fmin,N,1);
        Distance = pdist2(PopObj,PopObj);
        Distance(logical(eye(N))) = inf;
        radius   = max(min(Distance));
        % Detect subregion(s)
        Transformation = zeros(N,1);
        Remain         = find(Transformation==0);
        RegionID       = 1;
        while ~isempty(Remain)
            seeds      = find(~Transformation,1);
            Transformation(seeds) = RegionID;
            Remain     = find(Transformation==0);
            while true
                neighbors = sum(Distance(seeds,Remain)<=radius,1);
                seeds     = Remain(neighbors>=1);
                Transformation(seeds) = RegionID;
                Remain    = find(Transformation==0);
                if sum(neighbors)==0
                    break;
                end
            end
            RegionID = RegionID + 1;
        end
        %% Region division
        % Count the number of subregions of the true PF
        TrueNum = length(unique(Transformation));

        % Calculate the center point of each subregion
        Center = zeros(TrueNum,M);
        R      = ones(TrueNum,1);
        for i = 1 : TrueNum
            current     = Transformation==i;
            Center(i,:) = mean(PopObj(current,:));
            R(i)        = max(pdist2(PopObj(current,:),Center(i,:)));
        end

        % Select K points
        subNum  = tabulate(Transformation);
        subNum  = subNum(:,1:end-1);
        if TrueNum > K
            % Merging small subregions
            while sum(subNum(:,2)~=inf) > K
                [~,I] = min(subNum(:,2));
                Center(I,:) = inf(1,M);
                subNum(I,2) = inf;
                R(I)        = -inf;
                current     = find(Transformation == I);
                [~,T]       = min(pdist2(PopObj(current,:),Center),[],2);
                Transformation(current) = T;

                % Update reference point
                Idx = find(subNum(:,2)~=inf);
                for k =  1 : length(Idx)
                    Center(Idx(k),:) = mean(PopObj(Transformation == Idx(k),:));
                    R(Idx(k))        = max(pdist2(PopObj(Transformation == Idx(k),:),Center(Idx(k),:)))/sqrt(M-1);
                end
            end
        elseif TrueNum < K
            % Splite large subregions
            while sum(subNum(:,2)~=-inf) < K
                [~,I] = max(subNum(:,2));
                Center(I,:) = -inf(1,M);
                subNum(I,2) = -inf;
                R(I)        = -inf;
                current     = find(Transformation == I);
                [~,T1]      = max(pdist2(PopObj(current,:),PopObj(current(randi(length(current))),:)),[],1);
                [~,T2]      = max(pdist2(PopObj(current,:),PopObj(current(T1),:)),[],1);
                [~,T]       = min(pdist2(PopObj(current,:),PopObj(current([T1,T2],:),:)),[],2);
                ExistNum    = length(subNum(:,1));
                Transformation(current) = T + ExistNum;

                % Update reference point
                Center(ExistNum+1,:) = mean(PopObj(Transformation==ExistNum+1,:));
                Center(ExistNum+2,:) = mean(PopObj(Transformation==ExistNum+2,:));
                [R(ExistNum+1),R(ExistNum+2)] = deal(0.5*pdist2(Center(ExistNum+1,:),Center(ExistNum+2,:)));
                subNum(end+1,:) = [ExistNum+1,sum(T==1)];
                subNum(end+1,:) = [ExistNum+2,sum(T==2)];
            end
        end
    end

    % Select reference point
    select = abs(subNum(:,2)) ~= inf;
    Center = Center(select,:);
    R      = R(select);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function allMets = evaluateTracking(allSeq, resDir, dataDir)
%% evaluate CLEAR MOT and other metrics
% concatenate ALL sequences and evaluate as one!
%
% SETUP:
%
% define directories for tracking results...
% resDir = fullfile('res','data',filesep);
% ... and the actual sequences
% dataDir = fullfile('..','data','2DMOT2015','train',filesep);
%
%

fprintf('Sequences: \n');
disp(allSeq')

% concat gtInfo
gtInfo=[];
gtInfo.Xi=[];
allFgt=zeros(1,length(allSeq));

% Find out the length of each sequence
% and concatenate ground truth
gtInfoSingle=[];
seqCnt=0;
for s=allSeq
    seqCnt=seqCnt+1;
    seqName = char(s);
    seqFolder= [dataDir,seqName,filesep];

    assert(isdir(seqFolder),'Sequence folder %s missing',seqFolder);

    gtFile = fullfile(dataDir,seqName,'gt','gt.txt');
    gtI = convertTXTToStruct(gtFile,seqFolder);

    [Fgt,Ngt] = size(gtInfo.Xi);
    [FgtI,NgtI] = size(gtI.Xi);
    newFgt = Fgt+1:Fgt+FgtI;
    newNgt = Ngt+1:Ngt+NgtI;

    gtInfo.Xi(newFgt,newNgt) = gtI.Xi;
    gtInfo.Yi(newFgt,newNgt) = gtI.Yi;
    gtInfo.W(newFgt,newNgt) = gtI.W;
    gtInfo.H(newFgt,newNgt) = gtI.H;

    gtInfoSingle(seqCnt).wc=0;

    % fill in world coordinates if they exist
    if isfield(gtI,'Xgp') && isfield(gtI,'Ygp')
        gtInfo.Xgp(newFgt,newNgt) = gtI.Xgp;
        gtInfo.Ygp(newFgt,newNgt) = gtI.Ygp;
        gtInfoSingle(seqCnt).wc=1;
    end

    allFgt(seqCnt) = FgtI;

    gtInfoSingle(seqCnt).gtInfo=gtI;

end
gtInfo.frameNums=1:size(gtInfo.Xi,1);

allMets=[];

mcnt=1;

fprintf('Evaluating ... \n');

clear stInfo
stInfo.Xi=[];

evalMethod=1;

% flags for entire benchmark
% if one seq missing, evaluation impossible
eval2D=1;
eval3D=1;

seqCnt=0;

% iterate over each sequence
for s=allSeq

    seqCnt=seqCnt+1;
    seqName = char(s);

    fprintf('\t... %s\n',seqName);

    % if a result is missing, we cannot evaluate this tracker
    resFile = fullfile(resDir,[seqName '.txt']);
    if ~exist(resFile,'file')
        fprintf('WARNING: result for %s not available!\n',seqName);
        eval2D=0;
        eval3D=0;
        continue;
    end

    stI = convertTXTToStruct(resFile);
%     stI.Xi(find(stI.Xi(:)))=-1;
    % check if bounding boxes available in solution
    imCoord=1;
    if all(stI.Xi(find(stI.Xi(:)))==-1)
        imCoord=0;
    end

    worldCoordST=0; % state
    if isfield(stI,'Xgp') && isfield(stI,'Ygp')
        worldCoordST=1;
    end

    [FI,NI] = size(stI.Xi);

    % if stateInfo shorter, pad with zeros
    % GT and result must be equal length
    if FI<allFgt(seqCnt)
        missingFrames = FI+1:allFgt(seqCnt);
        stI.Xi(missingFrames,:)=0;
        stI.Yi(missingFrames,:)=0;
        stI.W(missingFrames,:)=0;
        stI.H(missingFrames,:)=0;
        if worldCoordST
            stI.Xgp(missingFrames,:)=0; stI.Ygp(missingFrames,:)=0;
        end
        [FI,NI] = size(stI.Xi);

    end

    % get result for one sequence only

    [mets, mInf]=CLEAR_MOT_HUN(gtInfoSingle(seqCnt).gtInfo,stI);

    allMets(mcnt).mets2d(seqCnt).name=seqName;
    allMets(mcnt).mets2d(seqCnt).m=mets;

    allMets(mcnt).mets3d(seqCnt).name=seqName;
    allMets(mcnt).mets3d(seqCnt).m=zeros(1,length(mets));

    if imCoord
        fprintf('*** 2D (Bounding Box overlap) ***\n'); printMetrics(mets); fprintf('\n');
    else
        fprintf('*** Bounding boxes not available ***\n\n');
        eval2D=0;
    end

    % if world coordinates available, evaluate in 3D
    if  gtInfoSingle(seqCnt).wc &&  worldCoordST
        evopt.eval3d=1;evopt.td=1;
        [mets, mInf]=CLEAR_MOT_HUN(gtInfoSingle(seqCnt).gtInfo,stI,evopt);
            allMets(mcnt).mets3d(seqCnt).m=mets;
        fprintf('*** 3D (in world coordinates) ***\n'); printMetrics(mets); fprintf('\n');
    else
        eval3D=0;
    end

    [F,N] = size(stInfo.Xi);
    newF = F+1:F+FI;
    newN = N+1:N+NI;

    % concat result
    stInfo.Xi(newF,newN) = stI.Xi;
    stInfo.Yi(newF,newN) = stI.Yi;
    stInfo.W(newF,newN) = stI.W;
    stInfo.H(newF,newN) = stI.H;
    if isfield(stI,'Xgp') && isfield(stI,'Ygp')
        stInfo.Xgp(newF,newN) = stI.Xgp;stInfo.Ygp(newF,newN) = stI.Ygp;
    end

end
stInfo.frameNums=1:size(stInfo.Xi,1);

if eval2D
    fprintf('\n');
    fprintf(' ********************* Your Benchmark Results (2D) ***********************\n');

    [m2d, mInf]=CLEAR_MOT_HUN(gtInfo,stInfo);
    allMets.bmark2d=m2d;

    evalFile = fullfile(resDir, 'eval2D.txt');

    printMetrics(m2d);
    dlmwrite(evalFile,m2d);
end

if eval3D
    fprintf('\n');
    fprintf(' ********************* Your Benchmark Results (3D) ***********************\n');

    evopt.eval3d=1;evopt.td=1;

    [m3d, mInf]=CLEAR_MOT_HUN(gtInfo,stInfo,evopt);
    allMets.bmark3d=m3d;

    evalFile = fullfile(resDir, 'eval3D.txt');

    printMetrics(m3d);
    dlmwrite(evalFile,m3d);
end
if ~eval2D && ~eval3D
    fprintf('ERROR: results cannot be evaluated\n');
end

author: yuxng


text:
function [phi_g, ip, x_g,  w_g] = boundary(i_bnd, x_local, n_gauss)

% Compute values of 2D quadratic bases at 1D Gauss points along a boundary
% Also return a list of (3) indices for the non-zero basis functions,
% the co-ords of the Gauss points and the 1D integration weights
%
% i_bnd specifies the triangle boundary of interest
%       1 - along s   = 0 boundary   - points [1 4 2]
%       2 - along r+s = 1 boundary   - points [2 5 3]
%       3 - along r   = 0 boundary   - points [3 6 1]
% x_local 6 x 2 array of nodal coordinate points (x, y)
% n_gauss  the number of Gauss points for the 1D integration

% phi_g is a n_gauss array x 3 - values of the 3 non-zero quadratic basis
%                                functions at the n_gauss 1D integration points
% ip  array of indices for non-zero quadratic functions
% x_g   is an n_gauss x 2 array of coordinates for the 1D Gauss points
% w_g   is an n_gauss array of integration weights

%% extract the Gauss weights and points on [- 1, 1];
     [r,  w_g] = oned_gauss(n_gauss); % 1 \le n_gauss \le 6

     r1 = (1+r)/2; % map to [0, 1]
     r1 = r1(:);   % make sure it's a column

%% which boundary is requested
     switch i_bnd
         case {1}
             ip = [ 1 4 2];
         case {2}
             ip = [ 2 5 3];
         case {3}
             ip = [ 3 6 1];
         otherwise
             disp(['Illegal i_bnd = ', int2str(i_bnd)])
             return
     end

%% Coordinates of the Gauss points
     x_g = (1-r1)*x_local(ip(1),:) ...
          +   r1 *x_local(ip(3),:);

%% Scale the weights by the (unsigned) length of the interval
     djac  =  norm(x_local(ip(3), :) - x_local(ip(1), :))/2;
     w_g   =  w_g*djac;

%% Evaluate the test functions

      p1 = 1 -3*r1 + 2*r1.^2; % p1(r = 0)  = 1
      p2 =    4*r1 - 4*r1.^2; % p2(r = 1/2)= 1
      p3 =   -  r1 + 2*r1.^2; % p3(r = 1)  = 1

      phi_g = [p1(:)  p2(:) p3(:) ];

  return
end

author: johannesgerer

repo: jburkardt-m