text:
function result = convolutionalPropagation(edges, edgeWeights, fixedVertices, ...
                    fixedDists,areaWeights,kernel,kernelTranspose,entropyLimit)

ne = size(edges,1);
nv = max(edges(:));
n = size(fixedDists,1);

% Initialize marginals
vv = ones(n,ne);
ww = ones(n,ne);
d = zeros(n,ne);

result = ones(n,nv);
result(:,fixedVertices) = fixedDists;

% For convenience, a binary array
isBoundary = zeros(nv,1);
isBoundary(fixedVertices) = 1;

% For each vertex compute his incoming and outgoing edges
incomingEdges = cell(nv,1);
outgoingEdges = cell(nv,1);

for i=1:ne
    e1 = edges(i,1);
    e2 = edges(i,2);
    outgoingEdges{e1} = [outgoingEdges{e1} i];
    incomingEdges{e2} = [incomingEdges{e2} i];
end

% No reason to project a distribution whose neighbors haven't been projected
queue = zeros(nv,1);
queue(fixedVertices) = 1;

projectEntropy = 0;

for j=1:1000
    toVisit = find(queue);

%     if mod(j,2)==0
%         toVisit = toVisit(end:-1:1); % reverse order each time for fun
%     end
%     toVisit = toVisit(randperm(length(toVisit))); % to avoid bias

    % Reset queue
    queue = zeros(nv,1);
    queue(fixedVertices) = 1;

    oldResult = result;

    for k=1:length(toVisit)
        v = toVisit(k);

        if isBoundary(v)
            p = result(:,v);

            for i=1:length(outgoingEdges{v}) %(v,w)
                e = outgoingEdges{v}(i);
                ww(:,e) = p ./ kernelTranspose(vv(:,e).*areaWeights);
            end

            for i=1:length(incomingEdges{v}) %(w,v)
                e = incomingEdges{v}(i);
                vv(:,e) = p ./ kernel(ww(:,e).*areaWeights);
            end
        else % not boundary
            omega = sum(edgeWeights(outgoingEdges{v}))+sum(edgeWeights(incomingEdges{v}));
            p = ones(n,1);

            for i=1:length(outgoingEdges{v}) %(v,w)
                e = outgoingEdges{v}(i);
                d(:,e) = ww(:,e).*kernelTranspose(vv(:,e).*areaWeights);
                d(:,e) = max(d(:,e),1e-10);
                p = p .* d(:,e).^(edgeWeights(e)/omega);
            end

            for i=1:length(incomingEdges{v}) %(w,v)
                e = incomingEdges{v}(i);
                d(:,e) = vv(:,e).*kernel(ww(:,e).*areaWeights);
                d(:,e) = max(d(:,e),0);
                p = p .* d(:,e).^(edgeWeights(e)/omega);
            end

            entropy = -sum(p.*log(p).*areaWeights);
            if nargin == 8 && entropy > entropyLimit && projectEntropy
                try % just ignore projection if it fails
                    fn = @(x) full(-sum(x*areaWeights.*((p.^x).*log(p))) - entropyLimit);
                    options = optimset('Display','none','tolfun',1e-6,'tolx',1e-6);
                    a = fzero(fn,[0 10],options);
                    p = p.^a;
                end
            end

            result(:,v) = p;

            for i=1:length(outgoingEdges{v}) %(v,w)
                e = outgoingEdges{v}(i);
                ww(:,e) = ww(:,e) .* p ./ d(:,e);
            end

            for i=1:length(incomingEdges{v}) %(w,v)
                e = incomingEdges{v}(i);
                vv(:,e) = vv(:,e) .* p ./ d(:,e);
            end
        end

        queue(edges(incomingEdges{v},:)) = 1;
        queue(edges(outgoingEdges{v},:)) = 1;
        queue(v) = 1;
    end

%     subplot(1,2,1);
%     imagesc(log(max(result,1e-20)));
%     colorbar;
%     title(sprintf('Iteration %d, log scale',j));
%     axis off;
%
%     subplot(1,2,2)
%     imagesc(result);
%     colorbar;
%     title(sprintf('Iteration %d, result',j));
%     axis off;
%
%     drawnow;

    change = sum(sum(bsxfun(@times,areaWeights,abs(result-oldResult))))/nv;
    change = full(change);

    fprintf('Iteration %d:  %g\n', j, change);

    if ~projectEntropy && change < 1e-3 && length(toVisit) == nv % was -3
        fprintf('Starting to project entropy.\n');
        projectEntropy = 1;
    elseif change < 1e-5 && j > 5 && projectEntropy
        return
    end
end

author: gpeyre


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Q = INVERSEKINEMATIC_IRB7600_500_230(robot, T)
%   Solves the inverse kinematic problem for the ABB IRB7600_500_2300 robot
%   where:
%   robot stores the robot parameters.
%   T is an homogeneous transform that specifies the position/orientation
%   of the end effector.
%
%   A call to Q=INVERSEKINEMATIC_IRB7600_500_2300 returns 8 possible solutions, thus,
%   Q is a 6x8 matrix where each column stores 6 feasible joint values.
%
%
%   Example code:
%
%   >>abb=load_robot('ABB', 'IRB7600_500_2300');
%   >>q = [0 0 0 0 0 0];
%   >>T = directkinematic(abb, q);
%   %Call the inversekinematic for this robot
%   >>qinv = inversekinematic(abb, T);
%
%   check that all of them are feasible solutions!
%   and every Ti equals T
%
%   for i=1:8,
%        Ti = directkinematic(abb, qinv(:,i))
%   end
%
%	See also DIRECTKINEMATIC.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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/>.
function [q] = inversekinematic_irb7600_500_230(robot, T)

%initialize q,
%eight possible solutions are generally feasible
q=zeros(6,8);

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%See geometry at the reference for this robot
L1=d(1);
L2=a(2);
L3=d(4);
A2=a(3);
L6=d(6);

A1 = a(1);

%T= [ nx ox ax Px;
%     ny oy ay Py;
%     nz oz az Pz];
Px=T(1,4);
Py=T(2,4);
Pz=T(3,4);

%Compute the position of the wrist, being W the Z component of the end effector's system
W = T(1:3,3);

% Pm: wrist position
Pm = [Px Py Pz]' - L6*W;

%first joint, two possible solutions admited:
% if q(1) is a solution, then q(1) + pi is also a solution
q1=atan2(Pm(2), Pm(1));

%solve for q2
q2_1=solve_for_theta2(robot, [q1 0 0 0 0 0 0], Pm);

q2_2=solve_for_theta2(robot, [q1+pi 0 0 0 0 0 0], Pm);

%solve for q3
q3_1=solve_for_theta3(robot, [q1 0 0 0 0 0 0], Pm);

q3_2=solve_for_theta3(robot, [q1+pi 0 0 0 0 0 0], Pm);

%Arrange solutions, there are 8 possible solutions so far.
% if q1 is a solution, q1* = q1 + pi is also a solution.
% For each (q1, q1*) there are two possible solutions
% for q2 and q3 (namely, elbow up and elbow up solutions)
% So far, we have 4 possible solutions. Howefer, for each triplet (theta1, theta2, theta3),
% there exist two more possible solutions for the last three joints, generally
% called wrist up and wrist down solutions. For this reason,
%the next matrix doubles each column. For each two columns, two different
%configurations for theta4, theta5 and theta6 will be computed. These
%configurations are generally referred as wrist up and wrist down solution
q = [q1         q1         q1        q1       q1+pi   q1+pi   q1+pi   q1+pi;
     q2_1(1)    q2_1(1)    q2_1(2)   q2_1(2)  q2_2(1) q2_2(1) q2_2(2) q2_2(2);
     q3_1(1)    q3_1(1)    q3_1(2)   q3_1(2)  q3_2(1) q3_2(1) q3_2(2) q3_2(2);
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0];

%leave only the real part of the solutions
q=real(q);

%Note that in this robot, the joint q3 has a non-simmetrical range. In this
%case, the joint ranges from 60 deg to -219 deg, thus, the typical normalizing
%step is avoided in this angle (the next line is commented). When solving
%for the orientation, the solutions are normalized to the [-pi, pi] range
%only for the theta4, theta5 and theta6 joints.

%normalize q to [-pi, pi]
q(1,:) = normalize(q(1,:));
q(2,:) = normalize(q(2,:));

% solve for the last three joints
% for any of the possible combinations (theta1, theta2, theta3)
for i=1:2:size(q,2),
    % use solve_spherical_wrist2 for the particular orientation
    % of the systems in this ABB robot
    % use either the geometric or algebraic method.
    % the function solve_spherical_wrist2 is used due to the relative
    % orientation of the last three DH reference systems.

    %use either one algebraic method or the geometric
    %qtemp = solve_spherical_wrist2(robot, q(:,i), T, 1, 'geometric'); %wrist up
    qtemp = solve_spherical_wrist2(robot, q(:,i), T, 1,'algebraic'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i)=qtemp;

    %qtemp = solve_spherical_wrist2(robot, q(:,i), T, -1, 'geometric'); %wrist down
    qtemp = solve_spherical_wrist2(robot, q(:,i), T, -1, 'algebraic'); %wrist down
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i+1)=qtemp;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for second joint theta2, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q2 = solve_for_theta2(robot, q, Pm)

%Evaluate the parameters
d = eval(robot.DH.d);
a = eval(robot.DH.a);

%See geometry
L2=a(2);
L3=d(4);
A2=a(3);

%compute L4
L4=sqrt(A2^2+L3^2);
%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

beta = atan2(-p1(2), p1(1));
gamma = (acos((L2^2+r^2-L4^2)/(2*r*L2)));

if ~isreal(gamma)
    disp('WARNING:inversekinematic_irb7600_500_230: the point is not reachable for this configuration, imaginary solutions');
%gamma = real(gamma);
end

%return two possible solutions
%elbow up and elbow down
%the order here is important and is coordinated with the function
%solve_for_theta3
q2(1) = pi/2 - beta - gamma; %elbow up
q2(2) = pi/2 - beta + gamma; %elbow down

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for third joint theta3, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q3 = solve_for_theta3(robot, q, Pm)

%Evaluate the parameters
d = eval(robot.DH.d);
a = eval(robot.DH.a);

%See geometry
L2=a(2);
L3=d(4);
A2=a(3);

%See geometry of the robot
L4=sqrt(A2^2+L3^2);
%the angle phi is fixed
phi= acos((A2^2-L3^2+L4^2)/(2*A2*L4));

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

eta = (acos((L2^2 + L4^2 - r^2)/(2*L2*L4)));

if ~isreal(eta)
   disp('WARNING:inversekinematic_irb7600_500_230: the point is not reachable for this configuration, imaginary solutions');
    %eta = real(eta);
end

%return two possible solutions
%elbow up and elbow down solutions
%the order here is important
q3(1) = pi - phi - eta;
q3(2) = pi - phi + eta;

repo: ARTE


text:
function varargout = drawArrow(varargin)
%DRAWARROW Draw an arrow on the current axis
%
%   drawArrow(x1, y1, x2, y2)
%   draws an arrow between the points (x1 y1) and (x2 y2).
%
%   drawArrow([x1 y1 x2 y2])
%   gives argument as a single array.
%
%   drawArrow(..., L, W)
%   specifies length and width of the arrow.
%
%   drawArrow(..., L, W, TYPE)
%   also specifies arrow type. TYPE can be one of the following :
%   0: draw only two strokes
%   1: fill a triangle
%   .5: draw a half arrow (try it to see ...)
%
%   Arguments can be single values or array of size N-by-1. In this case,
%   the function draws multiple arrows.
%
%   H = drawArrow(...)
%   return handle(s) to created arrow elements.
%   The handles are returned in a structure with the fields
%   'body', 'wing' and 'head' containing the handles to the different
%   parts of the arrow(s).
%
%   Example
%     t = linspace(0, 2*pi, 200);
%     figure; hold on;
%     plot(t, sin(t));
%     drawArrow([2 -1 pi 0], .1, .05, .5)
%
%   See also
%     drawEdge
%

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 11/11/2004 from drawEdge
%

%   HISTORY
%   2014-09-17 fix managment of handle values as suggested by Benoit Botton
%   2016-05-23 Improve codee and reduce calculations (by JuanPi Carbajal)

if isempty (varargin)
    error ('should specify at least one argument');
end

% parse arrow coordinates
var = varargin{1};
if size (var, 2) == 4
    x1 = var(:,1);
    y1 = var(:,2);
    x2 = var(:,3);
    y2 = var(:,4);
    varargin = varargin(2:end);

elseif length (varargin) > 3
    x1 = varargin{1};
    y1 = varargin{2};
    x2 = varargin{3};
    y2 = varargin{4};
    varargin = varargin(5:end);

else
    error ('MatGeom:drawArrow:invalidArgumentNumber', ...
        'wrong number of arguments, please read the doc');
end
N     = size (x1, 1);

% default values
l = 10  * ones (N, 1); % Body length
w = 5   * ones (N, 1); % Head width
r = 0.1 * ones (N, 1); % Head to body ratio
h = zeros (N, 1);      % Head type

if ~isempty (varargin)
    % Parse parameters
    k      = length (varargin);
    vartxt = 'lwrh';
    cmd    = ['%s = varargin{%d}; %s = %s(:);' ...
              'if length (%s) < N; %s = %s(1) * ones (N , 1); end'];
    for i = 1:k
        v = vartxt(i);
        eval (sprintf (cmd, v, i, v, v, v, v, v));
    end
end

hold on;
oldHold = ishold (gca);
if ~oldHold
    hold on;
end
axis equal;

% angle of the edge
theta = atan2 (y2-y1, x2-x1);

rl = r .* l;
rh = r .* h;
cT = cos (theta);
sT = sin (theta);
% point on the 'left'
xa1 = x2 - rl .* cT - w .* sT / 2;
ya1 = y2 - rl .* sT + w .* cT / 2;
% point on the 'right'
xa2 = x2 - rl .* cT + w .* sT / 2;
ya2 = y2 - rl .* sT - w .* cT / 2;
% point on the middle of the arrow
xa3 = x2 - rh .* cT;
ya3 = y2 - rh .* sT;

% draw main edge
tmp         = line ([x1.'; x2.'], [y1.'; y2.'], 'color', [0 0 1]);
handle.body = tmp;

% draw only 2 wings
ind = find (h == 0);
if ~isempty (ind)
    tmp              = line ([xa1(ind).'; x2(ind).'], [ya1(ind).'; y2(ind).'], ...
                             'color', [0 0 1]);
    handle.wing(:,1) = tmp;

    tmp              = line ([xa2(ind).'; x2(ind).'], [ya2(ind).'; y2(ind).'], ...
                             'color', [0 0 1]);
    handle.wing(:,2) = tmp;
end

% draw a full arrow
ind = find (h ~= 0);
if ~isempty (ind)
    tmp         = patch ([x2(ind) xa1(ind) xa3(ind) xa2(ind) x2(ind)].', ...
                         [y2(ind) ya1(ind) ya3(ind) ya2(ind) y2(ind)].', [0 0 1]);
    handle.head = tmp;
end

% format output arguments
if nargout > 0
    varargout{1} = handle;
end

if ~oldHold
    hold off;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [ nodeco, nodecotype, elnode, elnodetype, bdynde, bdyndetype ] = ...
  mesh2d_extract ( filename )

%*****************************************************************************80
%
%% MESH2D_EXTRACT extracts the information from the mesh file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILENAME, the name of the MAT file.
%
%    Output, real NODECO(NODE_NUM,2), the node coordinates.
%
%    Output, integer NODECOTYPE(NODE_NUM,1), the node type.
%
%    Output, integer ELNODE(ELEMENT_NUM,3), the node indices that form
%    each element.
%
%    Output, integer ELNODETYPE(ELEMENT_NUM,1), the element type.
%
%    Output, integer BDYNDE(BDY_NUM,2), pairs of nodes that form the boundary.
%
%    Output, integer BDYNDETYPE(BDY_NUM), the boundary condition type.
%
  input = fopen ( filename );

  if ( input < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MESH2D_EXTRACT - Error!\n' );
    fprintf ( 1, '  Could not open the file.\n' );
    error ( 'MESH2D_EXTRACT - Error!' );
  end

  line = fgets ( input );
  [ abc, count ] = sscanf ( line, '%d%d%d' );

  node_num = abc(1);
  element_num = abc(2);
  bdy_num = abc(3);

  nodeco = zeros ( node_num, 2 );
  nodecotype = zeros ( node_num, 1 );

  for node = 1 : node_num
    line = fgets ( input );
    [ abc, count ] = sscanf ( line, '%f%f%d' );
    nodeco(node,1) = abc(1);
    nodeco(node,2) = abc(2);
    nodecotype(node) = abc(3);
  end

  elnode = zeros ( element_num, 3 );
  elnodetype = zeros ( element_num,1 );

  for element = 1 : element_num
    line = fgets ( input );
    [ abcd, count ] = sscanf ( line, '%d%d%d%d' );
    elnode(element,1) = abcd(1);
    elnode(element,2) = abcd(2);
    elnode(element,3) = abcd(3);
    elnodetype(element) = abcd(4);
  end

  bdynde = zeros ( bdy_num, 2 );
  bdyndetype = zeros ( bdy_num, 1 );

  for bdy = 1 : bdynde
    line = fgets ( input );
    [ abc, count ] = sscanf ( line, '%d%d%d' );
    bdynde(node,1) = abc(1);
    bdynde(node,2) = abc(2);
    bdyndetype(node) = abc(3);
  end

  fclose ( input );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%Full Length simulation test
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 1/12/2011
classdef Simulation_FullLength_Test < TestCase
    methods
        function this = Simulation_FullLength_Test(name)
            this = this@TestCase(name);
        end
    end

    methods
        function testSimulation(~)
            runSimulation('outDir', 'output/runSimulationTests', 'logToDisk', true, 'logToDb', false);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [sY] = smooth2(Y,va)
% smooth 2D image with gaussian kernel
% function [sY] = smooth2(Y,va)
% IN:
%   - Y: the input image
%   - va: the smoothing value

if ~exist('va','var') || isempty(va)
    va = 0;
end

if va > 0
    nv = 4;
    w = zeros(2*nv+1);
    for i=-nv:nv
        for j=-nv:nv
            w(i+nv+1,j+nv+1) = exp(-0.5.*(i.^2+j.^2)/va);
        end
    end
    w = w./sum(w(:));
    sY = conv2(Y,w,'same');
else
    sY = Y;
end

author: MBB-team

repo: VBA-toolbox


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: computes Porous Slip Velocity based on Darcy's Law,
%
%           U_porous = -alpha <F_Lag,hat{n}> / | d vec{X}_Lag /ds |
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Por_Mat,nX,nY] = please_Compute_Porous_Slip_Velocity(ds,xLag,yLag,porous_info,F_Lag)

% xLag: vector of x-Pts associated w/ x-Lagrangian pts
% yLag: vector of y-Pts associated w/ y-Lagrangian pts
% porous_info:  col 1: lag-ids for porous media
%               col 2: x-Lag pts for lag-ids
%               col 3: y-Lag pts for lag-ids
%               col 4: porosity coefficient

% # of porous media pts.
Np = length( porous_info(:,1) );

% Initialize storage
%Por_X = zeros(length(xLag),1);
%Por_Y = Por_X;

% Compute Lagrangian Derivatives
[xL_s,yL_s] = give_Me_Lagrangian_Derivatives(ds,Np,porous_info);

% Compute Normal Vector (unit normals)
[nX,nY,sqrtNorm] = give_Me_Lagrangian_Normal_Vectors(xL_s,yL_s);

% Compute Porous Slip Velocity
Up_X = - ( porous_info(:,4) ) .* F_Lag( porous_info(:,1) ,1).*nX ./ sqrtNorm;
Up_Y = - ( porous_info(:,4) ) .* F_Lag( porous_info(:,1) ,2).*nY ./ sqrtNorm;
Por_Mat = [Up_X Up_Y];

% Store porous slip velocities in appropriate vector for adding to current Lagrangian Velocity Computation
%Por_X( porous_info(:,1) ) = Up_X;
%Por_Y( porous_info(:,1) ) = Up_Y;
%Por_Mat = [Por_X Por_Y];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: computes Lagrangian Derivatives
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xL_s,yL_s] = give_Me_Lagrangian_Derivatives(ds,Np,porous_info)

xL = porous_info(:,2);   % x-Values for Porous points
yL = porous_info(:,3);   % y-Values for Porous points
c = porous_info(:,5);    % choice of where it falls on stencil

xL_s = zeros(Np,1);
yL_s = zeros(Np,1);

for i=1:Np
   if c(i) == -2
       xL_s(i) = ( -25/12*xL(i) + 4*xL(i+1) - 3*xL(i+2) + 4/3*xL(i+3) - 1/4*xL(i+4) ) / ds;
       yL_s(i) = ( -25/12*yL(i) + 4*yL(i+1) - 3*yL(i+2) + 4/3*yL(i+3) - 1/4*yL(i+4) ) / ds;

   elseif c(i) == -1
       xL_s(i) = ( -0.25*xL(i-1) - 5/6*xL(i) + 1.5*xL(i+1) - 0.5*xL(i+2) + 1/12*xL(i+3) ) / ds;
       yL_s(i) = ( -0.25*yL(i-1) - 5/6*yL(i) + 1.5*yL(i+1) - 0.5*yL(i+2) + 1/12*yL(i+3) ) / ds;

   elseif c(i) == 0
       xL_s(i) = ( 1/12*xL(i-2) - 2/3*xL(i-1) + 2/3*xL(i+1) - 1/12*xL(i+2) ) / ds;
       yL_s(i) = ( 1/12*yL(i-2) - 2/3*yL(i-1) + 2/3*yL(i+1) - 1/12*yL(i+2) ) / ds;

   elseif c(i) == 1
       xL_s(i) = ( -1/12*xL(i-3) + 0.5*xL(i-2) - 1.5*xL(i-1) + 5/6*xL(i) + 0.25*xL(i+1) ) / ds;
       yL_s(i) = ( -1/12*yL(i-3) + 0.5*yL(i-2) - 1.5*yL(i-1) + 5/6*yL(i) + 0.25*yL(i+1) ) / ds;

   elseif c(i) == 2
       xL_s(i) = ( 0.25*xL(i-4) - 4/3*xL(i-3) + 3*xL(i-2) - 4*xL(i-1) + 25/12*xL(i) ) / ds;
       yL_s(i) = ( 0.25*yL(i-4) - 4/3*yL(i-3) + 3*yL(i-2) - 4*yL(i-1) + 25/12*yL(i) ) / ds;

   else
      fprintf('\n\n');
      error('Error: Bad file format inside your .porous file!\n');
   end
end

% IF CLOSED STRUCTURE!
% for i=1:Np
%     if i==1
%        xL_s(1) = ( xL(2) - xL(end) ) / (2*ds);
%        yL_s(1) = ( yL(2) - yL(end) ) / (2*ds);
%     elseif i<Np
%        xL_s(i) = ( xL(i+1) - xL(i-1) ) / (2*ds);
%        yL_s(i) = ( yL(i+1) - yL(i-1) ) / (2*ds);
%     else
%        xL_s(i) = ( xL(1) - xL(end-1) ) / (2*ds);
%        yL_s(i) = ( yL(1) - yL(end-1) ) / (2*ds);
%     end
% end

% Matrices for calculating derivatives
%AL2 = [1 1 1 1 1; 0 1 2 3 4; 0 1 4 9 16; 0 1 8 27 64; 0 1 16 3^4 4^4];
%AL1 = [1 1 1 1 1; -1 0 1 2 3; 1 0 1 4 9; -1 0 1 8 27; 1 0 1 16 81];
%A0 = [1 1 1 1 1; -2 -1 0 1 2; 4 1 0 1 4; -8 -1 0 1 8; 16 1 0 1 16];
%AR1 = [1 1 1 1 1; -3 -2 -1 0 1; 9 4 1 0 1; -27 -8 -1 0 1; 81 16 1 0 1];
%AR2 = [1 1 1 1 1; -4 -3 -2 -1 0; 16 9 4 1 0; -64 -27 -8 -1 0; 256 81 16 1 0];
%rhs = [0 1 0 0 0]';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% % FUNCTION: compute Porous Connections
% %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% function connects_Por = give_Me_Porous_Connections(ds,Np,porous_info)
%
% xL = porous_info(:,2);
% yL = porous_info(:,3);
%
% start = 1;
% last = Np;
% prev = last;
% for i=1:Np
%
%     x0 = xL(i); y0 = yL(i);        % central node value
%     xm1= xL(prev); ym1 = yL(prev); % 'previous' values
%     xp1 =xL(i+1); yp1=yL(i+1);     % 'next' values
%
%     distL = sqrt( (x0-xm1)^2 + (y0-ym1)^2 );
%     distR = sqrt( (x0-xp1)^2 + (y0-yp1)^2 );
%
%
% end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: computes Lagrangian UNIT Normal Vectors
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [nX,nY,sqrtN] = give_Me_Lagrangian_Normal_Vectors(xL_s,yL_s)

sqrtN = sqrt( (xL_s).^2 + (yL_s).^2 );

nX = ( yL_s ) ./ sqrtN;
nY = ( -xL_s) ./ sqrtN;

author: nickabattista


text:
% snirfDownsample()
%
% A utility that let's you pick fNIRS files in a file dialog box that
% are then decimated to the specified frequency using a low pass filter as
% implemented by the matlab 'downsample' command.
%
% Note that original copies of the fNIRS files are not preserved. It is
% advised that you create a back up of your fNIRS files in case you want to
% recover the higher sample rate date.
% Modified from resample code
% Fixed the problem with edges
% Fixed the problem with s vector - Meryem Oct 2016
% Add code to allow non-integer downsampling factor - Meryem Oct 2018
% Code refined for Homer3 and fNIRS files

% Homer3 Note
% Downsample tool does not work for a flat file structure

function snirfDownsample()

global maingui;
if isempty(maingui)
    [files, pathnm] = uigetfile( '*.snirf', 'Pick the .snirf file', 'multiselect','on');
else
    fpath = [maingui.dataTree.currElem.path maingui.dataTree.currElem.name];
    [files, pathnm] = uigetfile( '*.snirf', 'Pick files to downsample', fpath, 'multiselect','on');
end

if ~iscell(files)
    if files==0
        return;
    end
end

[~,name,~] = fileparts(files);

fsn = inputdlg( 'Decrease the sampling rate of a sequence by a factor of', 'Downsample SNIRF files', 1 );
if isempty(fsn)
    return
end
fsn = str2num(fsn{1});

wd = cd;
cd(pathnm)

if ~iscell(files)
    foo{1} = files;
    files = foo;
end

for iFile = 1:length(files)
    snirfData = SnirfClass(files{iFile});
%     varLst = whos();

%     fs = 1/(snirfData.data.time(2)-snirfData.data.time(1));

    if floor(fsn) == fsn % integer check

        snirfData.data.dataTimeSeries = downsample(snirfData.data.dataTimeSeries,fsn);
        if ~isempty(snirfData.aux)
            for iAux = 1:length(snirfData.aux)
                snirfData.aux(iAux).dataTimeSeries = downsample(snirfData.aux(iAux).dataTimeSeries,fsn);
                snirfData.aux(iAux).time = downsample(snirfData.aux(iAux).time,fsn);
            end
        end
        snirfData.data.time = downsample(snirfData.data.time,fsn);
%         s_sampled = zeros(size(snirfData.data.time,1),length(snirfData.aux));
%         for j=1:size(snirfData.stim,2)
%             lst_1 = snirfData.stim(j).data(:,3) == 1;
% %             lst_1 = round(lst_1/fsn);
%             s_sampled(lst_1,j) = 1;
%             lst_2 = snirfData.stim(j).data(:,3) == -1;
% %             lst_2 = round(lst_2/fsn);
%             s_sampled(lst_2,j) = -1;
%             snirfData.stim(j).data = s_sampled(:,j);
%         end
    else  % if downsample factor is not an integer (first upsample then downsample)
        t_new = linspace(1, size(snirfData.data.dataTimeSeries,1), 10*size(snirfData.data.dataTimeSeries,1));
        snirfData.data.dataTimeSeries = interp1(snirfData.data.dataTimeSeries, t_new);
        snirfData.data.dataTimeSeries = downsample(snirfData.data.dataTimeSeries,round(fsn*10));

        if ~isempty(snirfData.aux)
            for iAux = 1:length(snirfData.aux)
                snirfData.aux(iAux).dataTimeSeries = interp1(snirfData.aux(iAux).dataTimeSeries,t_new);
                snirfData.aux(iAux).dataTimeSeries = downsample(snirfData.aux(iAux).dataTimeSeries,round(fsn*10))';
                snirfData.aux(iAux).time = interp1(snirfData.aux(iAux).time,t_new);
                snirfData.aux(iAux).time = downsample(snirfData.aux(iAux).time,round(fsn*10))';
            end
        end

        snirfData.data.time = interp1(snirfData.data.time, t_new);
        snirfData.data.time = downsample(snirfData.data.time,round(fsn*10))';

%         s_sampled = zeros(size(snirfData.data.time,1),length(snirfData.aux));
%         for j=1:size(snirfData.stim,2)
%             lst_1 = snirfData.stim(j).data(:,3) == 1;
% %             lst_1 = round(lst_1/fsn);
%             s_sampled(lst_1,j) = 1;
%             lst_2 = snirfData.stim(j).data(:,3) == -1;
% %             lst_2 = round(lst_2/fsn);
%             s_sampled(lst_2,j) = -1;
%             snirfData.stim(j).data = s_sampled(:,j);
%         end

    end

    snirfName = sprintf([name '_downsample_' num2str(fsn) '.snirf']);
    snirfData.Save(snirfName);
    msgbox(['File created with name' snirfName]);
end

new_name = [files{1} '.orig'];
movefile (files{1}, new_name);

cd(wd);

author: BUNPC


text:
% ----------------------------------------------------------------------------
% function hfssCylinder(fid, Name, Axis, Center, Radius, Height, Units)
%
% Description :
% -------------
% Creates the VB script necessary to model a cylinder in HFSS.
%
% Parameters :
% ------------
% fid     - file identifier of the HFSS script file.
% Name    - name of the cylinder (in HFSS).
% Center  - center of the cylinder (specify as [x, y, z]). This is also the
%           starting point of the cylinder.
% Axis    - axis of the cylinder (specify as 'X', 'Y', or 'Z').
% Radius  - radius of the cylinder (scalar).
% Height  - height of the cylidner (from the point specified by Center).
% Units   - specify as 'in', 'mm', 'meter' or anything else defined in HFSS.
%
% Note :
% ------
%
% Example :
% ---------
% fid = fopen('myantenna.vbs', 'wt');
% ...
% hfssCylinder(fid, 'Cyl1', 'Z', [0, 0, 0], 0.1, 10, 'in');
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% This file is part of HFSS-MATLAB-API.
%
% HFSS-MATLAB-API 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.
%
% HFSS-MATLAB-API 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
% Foobar; if not, write to the Free Software Foundation, Inc., 59 Temple
% Place, Suite 330, Boston, MA  02111-1307  USA
%
% Copyright 2004, Vijay Ramasami (rvc@ku.edu)
% ----------------------------------------------------------------------------
function hfssCylinder(fid, Name, Axis, Center, Radius, Height, Units)

% Cylinder Parameters.
fprintf(fid, '\n');
fprintf(fid, 'oEditor.CreateCylinder _\n');
fprintf(fid, 'Array("NAME:CylinderParameters", _\n');
fprintf(fid, '"XCenter:=", "%f%s", _\n', Center(1), Units);
fprintf(fid, '"YCenter:=", "%f%s", _\n', Center(2), Units);
fprintf(fid, '"ZCenter:=", "%f%s", _\n', Center(3), Units);
fprintf(fid, '"Radius:=", "%f%s", _\n', Radius, Units);
fprintf(fid, '"Height:=", "%f%s", _\n', Height, Units);
fprintf(fid, '"WhichAxis:=", "%s"), _\n', upper(Axis));

% Cylinder Properties.
fprintf(fid, 'Array("NAME:Attributes", _\n');
fprintf(fid, '"Name:=", "%s", _\n', Name);
fprintf(fid, '"Flags:=", "", _\n');
fprintf(fid, '"Color:=", "(132 132 193)", _\n');
fprintf(fid, '"Transparency:=", 0, _\n');
fprintf(fid, '"PartCoordinateSystem:=", "Global", _\n');
fprintf(fid, '"MaterialName:=", "vacuum", _\n');
fprintf(fid, '"SolveInside:=", true)\n');
fprintf(fid, '\n');

author: yuip

repo: hfss-api

save_path: github-repos/MATLAB/yuip-hfss-api


text:
function [ ~ ] = intergrate(fname,fname_out,model,FBAsolution,listRxn)
%
% USAGE:
%
%    intergrate(fname, fname_out, model, FBAsolution, listRxn)
%
% INPUTS:
%    fname:               The name fo the file
%    fname_out:           The name of the output XML file name.
%    model:               A COBRA model structure
%    FBAsolution:         FBA results of a COBRA simulation by the
%                         `optimizeCbModel` function
%    listRxn:             A list of reaction IDs, after which the flux values are modified
%                         or added. P.s., the reaction IDs must be present in both the
%                         parsed CD model and the COBRA model structures. the reactions
%                         can be examined by `cmpR` function.
%
% OUTPUT:
%    fname_out:           A new XML file with the width of reaction links
%                         changed according to flux values calculated by COBRA
%
% EXAMPLE:
%
%    intergrate('anno_test.xml', 'anno_test_2_ecoli.xml', model_ecoli,solution, model_ecoli.rxns(:))
%
% .. Author: - Longfei Mao April, 2015

FBA_result=FBAsolution;
[parsed] = parseCD(fname);   % Parse a CellDesigner XML file into a CD Matlab structure
[ref,normalizedFlux] = addFlux(model,FBA_result,parsed,listRxn) % Overaly the flux distribution with the CD Matlab structure
[annotedText] = writeXML(fname,parsed,fname_out); % Write the CD Matlab strcuture to a new CD XML file.

author: opencobra

repo: cobratoolbox


text:
function y = cvx_isaffine( x, full ) %#ok
narginchk(1,2);
if nargin == 1,
    y = true;
else
    y = true( size( x ) );
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function rv = mptestver(varargin)
%MPTESTVER  Prints or returns MP-Test version info.
%   V = MPTESTVER returns the current MP-Test version numbers.
%   V = MPTESTVER('all') returns a struct with the fields Name, Version,
%   Release and Date (all char arrays). Calling MPTESTVER without assigning
%   the return value prints the version and release date of the current
%   installation of MP-Test.

%   MP-Test
%   Copyright (c) 2010-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Test.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mptest for more info.

v = struct( 'Name',     'MP-Test', ...
            'Version',  '8.0b1', ...
            'Release',  '', ...
            'Date',     '12-Dec-2022' );
if nargout > 0
    if nargin > 0
        rv = v;
    else
        rv = v.Version;
    end
else
    fprintf('%-22s Version %-10s %-11s\n', v.Name, v.Version, v.Date);
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function varargout = specgramscope(varargin)
% SPECGRAMSCOPE   Live updating spectrogram display
%
% STEP 1: Initialize the scope
% SPECGRAMSCOPE(FS,NFFT) initializes a spectrogram scope in the current axes.
% This spectrogram scope will compute and displays the NFFT-point FFT of a vector
% signal with sample rate FS Hz.  It will show the last 10 FFTs
%
% STEP 2: Update the scope
% SPECGRAMSCOPE(S) updates the spectrogram scope in the current axes with the
% FFT of vector S.  The scope should first be initialized as above with
% sample rate and FFT length.  If not, the sample rate will be 1 Hz and the FFT
% length will be the length of S.  Differences between the length of S and
% the specified FFT length are handled the same as MATLAB's built-in FFT
% function (i.e., zero-padding or truncation, as appropriate).
%
% SPECGRAMSCOPE(FS,NFFT,NTRACES) initializes a spectrogram scope in the
% current axes with NTRACES traces.  This is how many records to keep in
% time domain.  Default = 10
%
% SPECGRAMSCOPE(HAX, ...) defines the scope in specified axes HAX instead of GCA.  i.e.,
% SPECGRAMSCOPE(HAX,FS,NFFT) initializes axes HAX as a spectrogram scope, and
% SPECGRAMSCOPE(HAX,S) updates axes HAX with vector S.
%
% HAX = SPECGRAMSCOPE(...) returns a handle to the axes initialized by the
% spectrogram scope.  This is useful if you allow SPECGRAMSCOPE to create an
% axes for you, and want to be able to easily reference the axes for
% updates.  The surface created by SPECGRAMSCOPE all have the tag
% 'SpecgramScope'.  If you would like to manually modify the properties of
% these lines, their handles can be found by:
%
%        HAX = SPECGRAMSCOPE(...);
%        HSurf = findobj(HAX,'Tag','SpecgramScope');
%
% Example
%         %% Initialize data
%         Fs = 16384;
%         Nfft = 2048;
%         t = (0:1:Nfft-1)'/Fs;
%         fo = logspace(3.5,3.7);         % Range of fundamental frequencies
%         s1 = sin(2*pi*t*fo) + .1*rand(Nfft,length(fo));
%
%
%         %% Initialize scope
%         specgramscope(Fs,Nfft,30);
%         view([103 30])
%
%         %% Update scope
%         for ii = 1:length(fo)
%             specgramscope(s1(:,ii));
%             drawnow;pause(.01);
%         end;

%    Scott Hirsch 5-05.
%    shirsch@mathworks.com
%    Copyright 2004-2005 The MathWorks, Inc.

%% Parse input arguments
% Decision tree:
% + Initialize or update?
%   o If update -> OK
%   o If initialize -> Axes specified, or use GCA?

error(nargchk(1,4,nargin))
NTracesDefault = 10;

%% Initialize or update?
% If first or second input argument is not a scalar, it must be data - i.e. we are
% updating

if prod(size(varargin{1})) > 1 | prod(size(varargin{2})) > 1 % Update
    action = 'update';

    if nargin==1                % Use current axes
        hAxes = gca;
        data = varargin{1};
    else
        hAxes = varargin{1};    % Axes was specified
        data = varargin{2};
    end;

    % If the user has not initialized this scope, do it for them
    parms = getappdata(hAxes,'SpecgramScopeParameters');

    % Ensure that scope has been initialized
    if isempty(parms)
        % Use default values
        Fs = 1;
        data = rowmajor(data);
        Nfft = length(data);
        NTraces = NTracesDefault;       % Number of time histories
        feval(mfilename,hAxes,Fs,Nfft,NTraces);       % This recursive call will initialize the scope
        % Get the new parameter structure
        parms = getappdata(hAxes,'SpecgramScopeParameters');
    end;

else                                    % Initialize
    action = 'init';

    if ~isaxes(varargin{1})             % Easy mode, no handle passed in
        % Use current axes
        hAxes = gca;
        Fs = varargin{1};
        Nfft = varargin{2};
        if nargin==3
            NTraces = varargin{3};
        else
            NTraces = NTracesDefault;
        end;

    else                                % Expert mode, passed handle in
        hAxes = varargin{1};
        Fs = varargin{2};
        Nfft = varargin{3};
        if nargin==4
            NTraces = varargin{4};
        else
            NTraces = NTracesDefault;
        end;
    end;
end;

%% Dole out the work
%
switch action
    case 'init'     % Initialize

        % Build structure to internally pass information
        parms.Fs = Fs;                      % Sample Rate
        parms.NTraces = NTraces;            % Number of records in time
        parms.hAxes = hAxes;                % Handle to axes
        parms.Nfft = Nfft;                  % FFT Block size

        % Store parameter structure
        setappdata(hAxes,'SpecgramScopeParameters',parms);

        localInitScope(parms)               % Initialize scope

    case 'update'   % Update
        parms = getappdata(hAxes,'SpecgramScopeParameters');

        % Error checking
        % Ensure that scope has been initialized.  This shouldn't slip
        % through to here.
        if isempty(parms)
            error(['The spectrogram scope must first be initialized ' ...
                    'with the sample rate: specgramscope(hAxes,Fs)']);
        end;

        % Force data to be in columns.  Allow for multiple columns.  This will
        % error if data actually has more channels than samples.
        data = rowmajor(data);

        % Check that the number of columns corresponds to the number of lines
        nc = size(data,2);      % Number of columns
        if nc ~= 1
            error(['spectrogram scope requires a single column vector of data']);
        end;

        localUpdateScope(data,parms)            % Update the scope
end;

% Return appropriate output argument
if nargout
    varargout{1} = parms.hAxes;
end;

% ***********************************************************************
% Initialize the Scope
function localInitScope(parms)

% Set axes
f = (0:parms.Nfft/2-1)*parms.Fs/parms.Nfft;
f = f(:);

% Add surface
[X,Y] = meshgrid(-parms.NTraces+1:0,f);
parms.hSurf = surf(parms.hAxes,X,Y,NaN*ones(length(f),parms.NTraces), ...
    'Tag','SpecgramScope');

set(parms.hAxes, ...
    'XLim',[-parms.NTraces+1 0], ...
    'YLim',[0 f(end)]);%, ...
%     'YDir','reverse');  % Reverse frequency axes direction
shading(parms.hAxes,'interp');

setappdata(parms.hAxes,'SpecgramScopeParameters',parms);

%% Get handle to the figure
% Turn doublebuffer on to eliminate flickering
hFig = get(parms.hAxes,'Parent');

% In R14, it's possible that hFig would return a handle to a panel, not a
% figure
if ~strcmp(get(hFig,'Type'),'figure')
    hFig = get(hFig,'Parent');
end;

%%
% Label the plot.
% There's a bug in R13 when creating xlabel and ylabel with direct
% parenting - the alignment gets all messed up. Instead, make hAx current
% axes
ca = gca;
set(hFig,'CurrentAxes',parms.hAxes);
xlabel('History')
ylabel('Frequency (Hz)');
zlabel('Magnitude (dB)');
set(hFig,'CurrentAxes',ca);
view([103 30])

%%
% Turn doublebuffer on to eliminate flickering
set(hFig,'DoubleBuffer','on');

% ***********************************************************************
% Update the plot.
function localUpdateScope(data,parms)

[f,mag] = localfft(data,parms);

% Dynamically modify Magnitude axis as we go.  Expand, but don't shrink.
maxM=max(mag(:));
minM=min(mag(:));
yax2=get(parms.hAxes,'YLim');
if minM<yax2(1),
   yax2(1)=minM;
end
if maxM>yax2(2),
   yax2(2)=maxM;
end
set(parms.hAxes,'YLim',yax2)

% Update the plot
hSurf = parms.hSurf;
zd = get(hSurf,'ZData');
zd = [mag zd(:,1:end-1)];
set(hSurf,'ZData',zd,'CData',zd)

% set(parms.hLine, 'XData', f(:,1), 'YData', mag(:,1));
% set(parms.hLine, {'YData'}, Mag');

% Note: It looks like it's faster to update one line at a time
%  in a loop than to update with a cell array

% ***********************************************************************
% Calculate the fft of the data.
function [f, mag] = localfft(data,parms)

% Calculate the fft of the data.
xfft = 2/parms.Nfft*fft(data,parms.Nfft);

% Avoid taking the log of 0.
xfft(xfft == 0) = 1e-17;

% Compute magnitude, dB
mag = 20*log10(abs(xfft(1:parms.Nfft/2,:)));

f = (0:length(mag)-1)*parms.Fs/parms.Nfft;
f = f(:);

% ***********************************************************************
% Utility - isaxes
function truefalse = isaxes(h);
% ISAXES(H)  True if H is a handle to a valid axes

truefalse = 0;      % Start false
if ishandle(h)
    if strcmp('axes',get(h,'Type'))
        truefalse = 1;
    end;
end;

% ***********************************************************************
% Utility - rowmajor
function data = rowmajor(data);
% Force data to be row major. i.e. more rows than columns

[nr,nc] = size(data);
if nc>nr
    data = data';
    [nr,nc] = size(data);
end;

author: Sable

repo: mcbench-benchmarks


text:
function [parcel] = ft_sourceparcellate(cfg, source, parcellation)

% FT_SOURCEPARCELLATE combines the source-reconstruction parameters over the parcels, for
% example by averaging all the values in the anatomically or functionally labeled parcel.
%
% Use as
%    output = ft_sourceparcellate(cfg, source, parcellation)
% where the input source is a 2D surface-based or 3-D voxel-based source grid that was for
% example obtained from FT_SOURCEANALYSIS or FT_COMPUTE_LEADFIELD. The input parcellation is
% described in detail in FT_DATATYPE_PARCELLATION (2-D) or FT_DATATYPE_SEGMENTATION (3-D) and
% can be obtained from FT_READ_ATLAS or from a custom parcellation/segmentation for your
% individual subject. The output is a channel-based representation with the combined (e.g.
% averaged) representation of the source parameters per parcel.
%
% The configuration "cfg" is a structure that can contain the following fields
%   cfg.method       = string, method to combine the values, see below (default = 'mean')
%   cfg.parcellation = string, fieldname that contains the desired parcellation
%   cfg.parameter    = cell-array with strings, fields that should be parcellated (default = 'all')
%
% The values within a parcel or parcel-combination can be combined with different methods:
%   'mean'      compute the mean
%   'median'    compute the median (unsupported for fields that are represented in a cell-array)
%   'eig'       compute the largest eigenvector
%   'min'       take the minimal value
%   'max'       take the maximal value
%   'maxabs'    take the signed maxabs value
%   'std'       take the standard deviation
%
% See also FT_SOURCEANALYSIS, FT_DATATYPE_PARCELLATION, FT_DATATYPE_SEGMENTATION

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar source parcellation
ft_preamble provenance source parcellation

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% get the defaults
cfg.parcellation = ft_getopt(cfg, 'parcellation');
cfg.parameter    = ft_getopt(cfg, 'parameter', 'all');
cfg.method       = ft_getopt(cfg, 'method', 'mean'); % can be mean, min, max, svd
cfg.feedback     = ft_getopt(cfg, 'feedback', 'text');

% the data can be passed as input argument or can be read from disk
hasparcellation = exist('parcellation', 'var');

if ischar(cfg.parameter)
  cfg.parameter = {cfg.parameter};
end

if hasparcellation
  % the parcellation is specified as separate structure
else
  % the parcellation is represented in the source structure itself
  parcellation = source;
end

% keep the transformation matrix
if isfield(parcellation, 'transform')
  transform = parcellation.transform;
else
  transform = [];
end

% ensure it is a parcellation, not a segmentation
parcellation = ft_checkdata(parcellation, 'datatype', 'parcellation', 'parcellationstyle', 'indexed', 'hasunit', 'yes');

% keep the transformation matrix
if ~isempty(transform)
  parcellation.transform = transform;
end

% ensure it is a source, not a volume
source = ft_checkdata(source, 'datatype', 'source', 'insidestyle', 'logical', 'hasunit', 'yes');

% ensure that the source and the parcellation are anatomically consistent
if ~strcmp(source.unit, parcellation.unit)
  ft_error('the units of the source and parcellation structure are not consistent, please use FT_SOURCEINTERPOLATE');
end

% ensure that the source and the parcellation are anatomically consistent
tolerance = 0.1 * ft_scalingfactor('mm', source.unit);
if ~isalmostequal(source.pos, parcellation.pos, 'abstol', tolerance)
  ft_error('the positions of the source and parcellation structure are not consistent, please use FT_SOURCEINTERPOLATE');
end

if isempty(cfg.parcellation)
  % determine the first field that can be used for the parcellation
  fn = fieldnames(parcellation);
  for i=1:numel(fn)
    if isfield(parcellation, [fn{i} 'label'])
      ft_warning('using "%s" for the parcellation', fn{i});
      cfg.parcellation = fn{i};
      break
    end
  end
end

if isempty(cfg.parcellation)
  ft_error('you should specify the field containing the parcellation');
end

% determine the fields and corresponding dimords to work on
fn = fieldnames(source);
fn = setdiff(fn, {'pos', 'tri', 'inside', 'outside', 'time', 'freq', 'dim', 'transform', 'unit', 'coordsys', 'cfg', 'hdr'}); % remove fields that do not represent the data
fn = fn(cellfun(@isempty, regexp(fn, 'dimord'))); % remove dimord fields
fn = fn(cellfun(@isempty, regexp(fn, 'label'))); % remove label fields
dimord = cell(size(fn));
for i=1:numel(fn)
  dimord{i} = getdimord(source, fn{i});
end

if any(strcmp(cfg.parameter, 'all'))
  cfg.parameter = fn;
else
  [inside, i1, i2] = intersect(cfg.parameter, fn);
  [outside       ] = setdiff(cfg.parameter, fn);
  if ~isempty(outside)
    ft_warning('\nparameter "%s" cannot be parcellated', outside{:});
  end
  cfg.parameter = fn(i2);
  fn     = fn(i2);
  dimord = dimord(i2);
end

% although it is technically feasible, don't parcellate the parcellation itself
sel    = ~strcmp(cfg.parcellation, fn);
fn     = fn(sel);
dimord = dimord(sel);

if numel(fn)==0
  ft_error('there are no source parameters that can be parcellated');
end

% get the parcellation and the labels that go with it
tissue      = parcellation.(cfg.parcellation);
tissuelabel = parcellation.([cfg.parcellation 'label']);
ntissue     = length(tissuelabel);

if isfield(source, 'inside')
  % determine the conjunction of the parcellation and the inside source points
  n0 = numel(source.inside);
  n1 = sum(source.inside(:));
  n2 = sum(tissue(:)~=0);
  fprintf('there are in total %d positions, %d positions are inside the brain, %d positions have a label\n', n0, n1, n2);
  fprintf('%d of the positions inside the brain have a label\n',        sum(tissue(source.inside)~=0));
  fprintf('%d of the labeled positions are inside the brain\n',         sum(source.inside(tissue(:)~=0)));
  fprintf('%d of the positions inside the brain do not have a label\n', sum(tissue(source.inside)==0));
  % discard the positions outside the brain and the positions in the brain that do not have a label
  tissue(~source.inside) = 0;
end

% start preparing the output data structure
parcel       = keepfields(source, {'freq','time','cumtapcnt'});
parcel.label = tissuelabel;

for i=1:numel(fn)
  % parcellate each of the desired parameters
  dat = source.(fn{i});
  siz = getdimsiz(source, fn{i});
  siz(contains(tokenize(dimord{i},'_'),'pos')) = ntissue;

  if startsWith(dimord{i}, '{pos_pos}')
    fprintf('creating %d*%d parcel combinations for parameter %s by taking the %s\n', numel(tissuelabel), numel(tissuelabel), fn{i}, cfg.method);
    tmp = zeros(siz);
    ft_progress('init', cfg.feedback, 'computing parcellation');
    k = 0;
    K = numel(tissuelabel)^2;
    for j1=1:numel(tissuelabel)
      for j2=1:numel(tissuelabel)
        k = k + 1;
        ft_progress(k/K, 'computing parcellation for %s combined with %s', tissuelabel{j1}, tissuelabel{j2});
        switch cfg.method
          case 'mean'
            tmp(j1,j2,:,:) = cellmean(dat(tissue==j1,tissue==j2));
          case 'median'
            tmp(j1,j2,:,:) = cellmedian(dat(tissue==j1,tissue==j2));
          case 'min'
            tmp(j1,j2,:,:) = cellmin(dat(tissue==j1,tissue==j2));
          case 'max'
            tmp(j1,j2,:,:) = cellmax(dat(tissue==j1,tissue==j2));
          case 'eig'
            tmp(j1,j2,:,:) = celleig(dat(tissue==j1,tissue==j2));
          case 'std'
            tmp(j1,j2,:,:) = cellstd(dat(tissue==j1,tissue==j2));
          otherwise
            ft_error('method %s not implemented for %s', cfg.method, dimord{i});
        end % switch
      end % for j2
    end % for j1
    ft_progress('close');

  elseif startsWith(dimord{i}, '{pos}')
    fprintf('creating %d parcels for parameter %s by taking the %s\n', numel(tissuelabel), fn{i}, cfg.method);
    tmp = zeros(siz);
    ft_progress('init', cfg.feedback, 'computing parcellation');
    for j=1:numel(tissuelabel)
      ft_progress(j/numel(tissuelabel), 'computing parcellation for %s', tissuelabel{j});
      switch cfg.method
        case 'mean'
          tmp(j,:,:) = cellmean(dat(tissue==j));
        case 'median'
          tmp(j,:,:) = cellmedian(dat(tissue==j));
        case 'min'
          tmp(j,:,:) = cellmin(dat(tissue==j));
        case 'max'
          tmp(j,:,:) = cellmax(dat(tissue==j));
        case 'eig'
          tmp(j,:,:) = celleig(dat(tissue==j));
        case 'std'
          tmp(j,:,:) = cellstd(dat(tissue==j));
        otherwise
          ft_error('method %s not implemented for %s', cfg.method, dimord{i});
      end % switch
    end % for
    ft_progress('close');

  elseif startsWith(dimord{i}, 'pos_pos')
    fprintf('creating %d*%d parcel combinations for parameter %s by taking the %s\n', numel(tissuelabel), numel(tissuelabel), fn{i}, cfg.method);
    siz     = size(dat);
    siz(1)  = ntissue;
    siz(2)  = ntissue;
    tmp     = nan(siz);
    ft_progress('init', cfg.feedback, 'computing parcellation');
    k = 0;
    K = numel(tissuelabel)^2;
    for j1=1:numel(tissuelabel)
      for j2=1:numel(tissuelabel)
        k = k + 1;
        ft_progress(k/K, 'computing parcellation for %s combined with %s', tissuelabel{j1}, tissuelabel{j2});
        switch cfg.method
          case 'mean'
            tmp(j1,j2,:) = arraymean2(dat(tissue==j1,tissue==j2,:));
          case 'median'
            tmp(j1,j2,:) = arraymedian2(dat(tissue==j1,tissue==j2,:));
          case 'min'
            tmp(j1,j2,:) = arraymin2(dat(tissue==j1,tissue==j2,:));
          case 'max'
            tmp(j1,j2,:) = arraymax2(dat(tissue==j1,tissue==j2,:));
          case 'eig'
            tmp(j1,j2,:) = arrayeig2(dat(tissue==j1,tissue==j2,:));
          case 'maxabs'
            tmp(j1,j2,:) = arraymaxabs2(dat(tissue==j1,tissue==j2,:));
          case 'std'
            tmp(j1,j2,:) = arraystd2(dat(tissue==j1,tissue==j2,:));
          otherwise
            ft_error('method %s not implemented for %s', cfg.method, dimord{i});
        end % switch
      end % for j2
    end % for j1
    ft_progress('close');

  elseif startsWith(dimord{i}, 'pos')
    fprintf('creating %d parcels for %s by taking the %s\n', numel(tissuelabel), fn{i}, cfg.method);
    siz     = size(dat);
    siz(1)  = ntissue;
    tmp     = nan(siz);
    ft_progress('init', cfg.feedback, 'computing parcellation');
    for j=1:numel(tissuelabel)
      ft_progress(j/numel(tissuelabel), 'computing parcellation for %s', tissuelabel{j});
      switch cfg.method
        case 'mean'
          tmp(j,:) = arraymean1(dat(tissue==j,:));
        case 'mean_thresholded'
          cfg.mean = ft_getopt(cfg, 'mean', struct('threshold', []));
          if isempty(cfg.mean.threshold)
            ft_error('when cfg.method = ''mean_thresholded'', you should specify a cfg.mean.threshold');
          end
          if numel(cfg.mean.threshold)==size(dat,1)
            % assume one threshold per vertex
            threshold = cfg.mean.threshold(tissue==j,:);
          else
            threshold = cfg.mean.threshold;
          end
          tmp(j,:) = arraymean1(dat(tissue==j,:), threshold);
        case 'median'
          tmp(j,:) = arraymedian1(dat(tissue==j,:));
        case 'min'
          tmp(j,:) = arraymin1(dat(tissue==j,:));
        case 'max'
          tmp(j,:) = arraymax1(dat(tissue==j,:));
        case 'maxabs'
          tmp(j,:) = arraymaxabs1(dat(tissue==j,:));
        case 'eig'
          tmp(j,:) = arrayeig1(dat(tissue==j,:));
        case 'std'
          tmp(j,:)  = arraystd1(dat(tissue==j,:));
        otherwise
          ft_error('method %s not implemented for %s', cfg.method, dimord{i});
      end % switch
    end % for
    ft_progress('close');

  else
    ft_error('unsupported dimord %s', dimord{i})

  end % if pos, pos_pos, {pos}, etc.

  % update the dimord, use chan rather than pos
  % this makes it look just like timelock or freq data
  tok = tokenize(dimord{i}, '_');
  tok(strcmp(tok,  'pos' )) = {'chan'}; % replace pos by chan
  tok(strcmp(tok, '{pos}')) = {'chan'}; % replace pos by chan
  tok(strcmp(tok, '{pos'))  = {'chan'}; % replace pos by chan
  tok(strcmp(tok, 'pos}'))  = {'chan'}; % replace pos by chan

  % squeeze out any singleton oris
  siz  = [size(tmp) 1]; % add trailing singleton to be sure
  oris = contains(tok, 'ori') & siz(1:numel(tok))==1;
  siz(oris) = [];
  tmp = reshape(tmp, siz);
  tok(oris) = [];

  tmpdimord = sprintf('%s_', tok{:});
  tmpdimord = tmpdimord(1:end-1);         % exclude the last _

  % store the results in the output structure
  parcel.(fn{i})            = tmp;
  parcel.([fn{i} 'dimord']) = tmpdimord;

  % to avoid confusion
  clear dat tmp tmpdimord j j1 j2
end % for each of the fields that should be parcellated

% a brainordinate is a brain location that is specified by either a surface vertex (node) or a volume voxel
parcel.brainordinate = keepfields(parcellation, {'pos', 'tri', 'dim', 'transform'}); % keep the information about the geometry
fn = fieldnames(parcellation);
for i=1:numel(fn)
  if isfield(parcellation, [fn{i} 'label'])
    % keep each of the labeled fields from the parcellation
    parcel.brainordinate.( fn{i}         ) = parcellation.( fn{i}         );
    parcel.brainordinate.([fn{i} 'label']) = parcellation.([fn{i} 'label']);
  end
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   source parcellation
ft_postamble provenance parcel
ft_postamble history    parcel
ft_postamble savevar    parcel

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTIONS to compute something over the first dimension
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = arraymean1(x, threshold)

if nargin==1
  y = mean(x,1);
else
  if numel(threshold)==1
    % scalar comparison is possible
  elseif size(threshold,1) == size(x,1)
    % assume threshold to be column vector
    threshold = repmat(threshold, [1, size(x,2)]);
  end
  sel = sum(x>threshold,2);
  if ~isempty(sel)
    y = mean(x(sel>0,:),1);
  else
    y = nan+zeros(1,size(x,2));
  end
end

function y = arraymedian1(x)
if ~isempty(x)
  y = median(x,1);
else
  y = nan(1,size(x,2));
end

function y = arraymin1(x)
if ~isempty(x)
  y = min(x,[], 1);
else
  y = nan(1,size(x,2));
end

function y = arraymax1(x)
if ~isempty(x)
  y = max(x,[], 1);
else
  y = nan(1,size(x,2));
end

function y = arrayeig1(x)
if ~isempty(x)
  siz = size(x);
  x = reshape(x, siz(1), prod(siz(2:end)));
  [u, s, v] = svds(x, 1);         % x = u * s * v'
  y = s(1,1) * v(:,1);            % retain the largest eigenvector with appropriate scaling
  y = reshape(y, [siz(2:end) 1]); % size should have at least two elements
else
  siz = size(x);
  y   = nan([siz(2:end) 1]);
end

function y = arraymaxabs1(x)
if ~isempty(x)
  % take the value that is at max(abs(x))
  [dum,ix] = max(abs(x), [], 1);
  y        = x(ix);
else
  y = nan(1,size(x,2));
end

function y = arraystd1(x)
if ~isempty(x)
  y = std(x,0, 1);
else
  y = nan(1,size(x, 2));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTIONS to compute something over the first two dimensions
% all of these functions should be implemented the same
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = arraymean2(x)
siz = size(x);
x = reshape(x, [siz(1)*siz(2) siz(3:end) 1]); % simplify it into a single dimension
y = arraymean1(x);

function y = arraymedian2(x)
siz = size(x);
x = reshape(x, [siz(1)*siz(2) siz(3:end) 1]); % simplify it into a single dimension
y = arraymedian1(x);

function y = arraymin2(x)
siz = size(x);
x = reshape(x, [siz(1)*siz(2) siz(3:end) 1]); % simplify it into a single dimension
y = arraymin1(x);

function y = arraymax2(x)
siz = size(x);
x = reshape(x, [siz(1)*siz(2) siz(3:end) 1]); % simplify it into a single dimension
y = arraymax1(x);

function y = arrayeig2(x)
siz = size(x);
x = reshape(x, [siz(1)*siz(2) siz(3:end) 1]); % simplify it into a single dimension
y = arrayeig1(x);

function y = arraymaxabs2(x)
siz = size(x);
x = reshape(x, [siz(1)*siz(2) siz(3:end) 1]); % simplify it into a single dimension
y = arraymaxabs1(x);

function y = arraystd2(x)
siz = size(x);
x = reshape(x, [siz(1)*siz(2) siz(3:end) 1]); % simplify it into a single dimension
y = arraystd1(x);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTIONS for doing something over all elements of a cell-array
% add a singleton dimension, concatenate into an array, and do the computatioon
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = cellmean(x)
siz = size(x{1});
for i=1:numel(x)
  x{i} = reshape(x{i}, [1 siz]);
end
x = cat(1, x{:});
y = arraymean1(x);

function y = cellmedian(x)
siz = size(x{1});
for i=1:numel(x)
  x{i} = reshape(x{i}, [1 siz]);
end
x = cat(1, x{:});
y = arraymedian1(x);

function y = cellmin(x)
siz = size(x{1});
for i=1:numel(x)
  x{i} = reshape(x{i}, [1 siz]);
end
x = cat(1, x{:});
y = arraymin1(x);

function y = cellmax(x)
siz = size(x{1});
for i=1:numel(x)
  x{i} = reshape(x{i}, [1 siz]);
end
x = cat(1, x{:});
y = arraymax1(x);

function y = celleig(x)
siz = size(x{1});
for i=1:numel(x)
  x{i} = reshape(x{i}, [1 siz]);
end
x = cat(1, x{:});
y = arrayeig1(x);

function y = cellstd(x)
siz = size(x{1});
for i=1:numel(x)
  x{i} = reshape(x{i}, [1 siz]);
end
x = cat(1, x{:});
y = arraystd1(x);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef MimModelMap < CoreBaseClass

    properties (Access = private)
        Mim
        ModelCache
    end

    methods
        function obj = MimModelMap(mim)
            obj.Mim = mim;
            obj.clear();
        end

        function [value, hash] = getValue(obj, modelId)
            model = obj.ModelCache(modelId);
            [value, hash] = model.getOrRun();
        end

        function value = getDependentValue(obj, modelId, dependentModelId)
            % Gets a model value and adds the caller as a dependent

            model = obj.ModelCache(modelId);
            model.addDependentModel(obj.ModelCache(dependentModelId));
            value = model.getOrRun();
        end

        function setValue(obj, modelId, value)
            model = obj.ModelCache(modelId);
            model.setValue(value);
        end

        function addItem(obj, modelId, itemId)
            model = obj.ModelCache(modelId);
            model.addItem(modelId, itemId);
        end

        function removeItem(obj, modelId, itemId)
            model = obj.ModelCache(modelId);
            model.removeItem(itemId);
        end

        function autoUpdate(obj)
            % Triggers updating of invalid models that have AutoUpdate

            for model = obj.ModelCache.values()
                model{1}.update();
            end
        end

        function modelId = buildModelId(obj, modelClassName, parameters)
            % ToDo: Check parameters so as not to duplicate identical models

            modelConstructor = str2func(modelClassName);
            modelId = CoreSystemUtilities.GenerateUid();
            model = modelConstructor();
            model.initialise(modelId, parameters, MimModelCallback(modelId, obj, obj.Mim), true);
            obj.ModelCache(modelId) = model;
        end

        function clear(obj)
            obj.ModelCache = containers.Map();
            % Explicitly add in models with unique names
            model = MimSubjectList();
            model.initialise('MimSubjectList', {}, MimModelCallback('MimSubjectList', obj, obj.Mim), true);
            obj.ModelCache('MimSubjectList') = model;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
classdef StressAndMaxNormWithExponentPlotter < StressNormPlotter

    methods (Access = public)

        function obj = StressAndMaxNormWithExponentPlotter(d)
            obj.init(d);
        end

        function plot(obj)
            obj.fig = figure(1);
            obj.plotSmoothInfo();
            obj.plotNonSmoothInfo();
            obj.addLegend();
            obj.print();
        end

    end

    methods (Access = private)

        function plotSmoothInfo(obj)
            dB = obj.smoothDB;
            obj.plotStressNorm(dB);
            obj.plotStressMax(dB);
        end

        function plotNonSmoothInfo(obj)
            dB = obj.nonSmoothDB;
            obj.plotStressNorm(dB);
            obj.plotStressMax(dB);
        end

        function plotStressNorm(obj,dB)
            obj.x        = dB.pNorm;
            obj.y        = dB.invStressNorm;
            obj.lineInfo = [dB.color,'-+'];
            obj.appendPlot()
        end

        function plotStressMax(obj,dB)
            obj.x            = dB.pNorm;
            obj.y            = dB.invStressMax*ones(1,length(obj.x));
            obj.lineInfo     = [dB.color,'--'];
            obj.appendPlot()
        end

        function addLegend(obj)
            legStr = {'Smooth','NonSmooth','SmoothMaxNorm','NonSmoothMaxNorm'};
            lines = obj.plotLines;
            nlinesType = length(legStr);
            ilines = 1:obj.nExp:obj.nExp*nlinesType;
            legLines = [lines{ilines(1)} lines{ilines(2)} lines{ilines(3)} lines{ilines(4)}];
            legend(legLines,legStr);
        end

        function print(obj)
            sFig = obj.fig;
            ca = get(sFig,'CurrentAxes');
            xl = get(ca,'xlabel');
            set(xl,'string','Pnorm');
            figureName = 'StressNormWithPnorm';
            outPutFigName = fullfile(obj.outPutPath,figureName);
            pL   = obj.plotLines;
            fp   = plotPrinter(sFig,pL);
            fp.print(outPutFigName)
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%% script for making plots used to demonstrate accuracy of detectors

%% set up path

addpath ../perframe;
addpath ../perframe/compute_perframe_features;
addpath ../misc;
addpath ../filehandling;
addpath C:\Code\Ctrax\trunk\matlab\netlab;
outfigdir = 'C:/Code/Jdetect/figures/AccuracyOut';
if ~exist(outfigdir,'dir'),
  mkdir(outfigdir);
end

%% global parameters

% groundtruthing
ifpcolor = [0,0,.7];
fpcolor = [.3,.3,1];
ifncolor = [.7,0,0];
fncolor = [1,.3,.3];
ierrorcolor = [0,0,0];
errorcolor = [.3,.3,.3];
chancecolor = [.7,.7,.7];

TOUCH = 1;
CRABWALK = 2;
BACKUP = 3;
CHASE = 4;
WINGGROOMING = 5;
JUMP = 6;
STOP = 7;
RIGHTING = 8;
TAILPIVOTTURN = 9;
WALK = 10;

MOUSE_WALK = 11;
MOUSE_MALEFOLLOW = 12;

LARVA_HEADCAST = 13;
LARVA_CRAWL = 14;

%% TOUCH

behavior = 'Touch';
configfilename = '../perframe/params/JLabelParams_FlyBowl_Touch.xml';
roottraindir = 'C:/Data/JAABA/FlyBowl';
rootgroundtruthdir = 'C:/Data/JAABA/groundtruth_pBDPGAL4U_data';

train_experiment_names = {
  'GMR_71G01_AE_01_TrpA_Rig2Plate17BowlD_20110921T084818'
  'pBDPGAL4U_TrpA_Rig1Plate15BowlB_20110922T145928'
  'pBDPGAL4U_TrpA_Rig2Plate14BowlD_20110615T164545'
  'pBDPGAL4U_TrpA_Rig2Plate14BowlD_20110707T085804'};

groundtruth_experiment_names = {
  'pBDPGAL4U_TrpA_Rig1Plate10BowlB_20110609T091905'
  'pBDPGAL4U_TrpA_Rig1Plate15BowlD_20120210T152130'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlA_20120202T145723'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlC_20110811T133805'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlD_20111130T094443'};

train_expdirs = cellfun(@(x) fullfile(roottraindir,x),train_experiment_names,'UniformOutput',false);
groundtruth_expdirs = cellfun(@(x) fullfile(rootgroundtruthdir,x),groundtruth_experiment_names,'UniformOutput',false);

[hfigs,errordata,figfilenames,outmatfilename] = ...
  Compute_GroundTruthError_vs_TrainingSetSize(behavior,configfilename,train_expdirs,groundtruth_expdirs,...
  'hfigs',[3,4]+100*TOUCH,...
  'figpos',[95 550 1246 443],...
  'ifpcolor',ifpcolor,...
  'fpcolor',fpcolor,...
  'ifncolor',ifncolor,...
  'fncolor',fncolor,...
  'ierrorcolor',ierrorcolor,...
  'errorcolor',errorcolor,...
  'chancecolor',chancecolor,...
  'outfigdir',outfigdir);

%% CRABWALK

behavior = 'Crabwalk';
configfilename = '../perframe/params/JLabelParams_FlyBowl_Crabwalk.xml';

roottraindir = 'C:/Data/JAABA/FlyBowl';
rootgroundtruthdir = 'C:/Data/JAABA/groundtruth_pBDPGAL4U_data';

train_experiment_names = {
  'GMR_71G01_AE_01_TrpA_Rig2Plate17BowlD_20110921T084818'
  'pBDPGAL4U_TrpA_Rig1Plate15BowlB_20110922T145928'
  'pBDPGAL4U_TrpA_Rig2Plate14BowlD_20110615T164545'
  'pBDPGAL4U_TrpA_Rig2Plate14BowlD_20110707T085804'};

groundtruth_experiment_names = {
  'pBDPGAL4U_TrpA_Rig1Plate10BowlB_20110609T091905'
  'pBDPGAL4U_TrpA_Rig1Plate15BowlD_20120210T152130'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlA_20120202T145723'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlC_20110811T133805'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlD_20111130T094443'};

train_expdirs = cellfun(@(x) fullfile(roottraindir,x),train_experiment_names,'UniformOutput',false);
groundtruth_expdirs = cellfun(@(x) fullfile(rootgroundtruthdir,x),groundtruth_experiment_names,'UniformOutput',false);

try

[hfigs,errordata,figfilenames,outmatfilename] = ...
  Compute_GroundTruthError_vs_TrainingSetSize(behavior,configfilename,train_expdirs,groundtruth_expdirs,...
  'hfigs',[3,4]+100*CRABWALK,...
  'figpos',[95 550 1246 443],...
  'ifpcolor',ifpcolor,...
  'fpcolor',fpcolor,...
  'ifncolor',ifncolor,...
  'fncolor',fncolor,...
  'ierrorcolor',ierrorcolor,...
  'errorcolor',errorcolor,...
  'chancecolor',chancecolor,...
  'outfigdir',outfigdir);

catch ME
  warning('%s: %s',behavior,getReport(ME));
end

%% backup

behavior = 'Backup';
configfilename = '../perframe/params/JLabelParams_FlyBowl_Backup.xml';

roottraindir = 'C:/Data/JAABA/FlyBowl';
rootgroundtruthdir = 'C:/Data/JAABA/groundtruth_pBDPGAL4U_data';

train_experiment_names = {
  'GMR_71G01_AE_01_TrpA_Rig2Plate17BowlD_20110921T084818'
  'pBDPGAL4U_TrpA_Rig1Plate15BowlB_20110922T145928'
  'pBDPGAL4U_TrpA_Rig2Plate14BowlD_20110615T164545'
  'pBDPGAL4U_TrpA_Rig2Plate14BowlD_20110707T085804'};

groundtruth_experiment_names = {
  'pBDPGAL4U_TrpA_Rig1Plate10BowlB_20110609T091905'
  'pBDPGAL4U_TrpA_Rig1Plate15BowlD_20120210T152130'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlA_20120202T145723'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlC_20110811T133805'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlD_20111130T094443'};

train_expdirs = cellfun(@(x) fullfile(roottraindir,x),train_experiment_names,'UniformOutput',false);
groundtruth_expdirs = cellfun(@(x) fullfile(rootgroundtruthdir,x),groundtruth_experiment_names,'UniformOutput',false);

try
[hfigs,errordata,figfilenames,outmatfilename] = ...
  Compute_GroundTruthError_vs_TrainingSetSize(behavior,configfilename,train_expdirs,groundtruth_expdirs,...
  'hfigs',3:4+BACKUP*100,...
  'figpos',[95 550 1246 443],...
  'ifpcolor',ifpcolor,...
  'fpcolor',fpcolor,...
  'ifncolor',ifncolor,...
  'fncolor',fncolor,...
  'ierrorcolor',ierrorcolor,...
  'errorcolor',errorcolor,...
  'chancecolor',chancecolor,...
  'outfigdir',outfigdir);
catch ME
  warning('%s: %s',behavior,getReport(ME));
end

%% chase

behavior = 'Chase';
configfilename = '../perframe/params/JLabelParams_FlyBowl_Chase.xml';

roottraindir = 'C:/Data/JAABA/FlyBowl';
rootgroundtruthdir = 'C:/Data/JAABA/groundtruth_pBDPGAL4U_data';

train_experiment_names = {
  'GMR_14C07_AE_01_TrpA_Rig1Plate15BowlA_20120404T141155'
  'GMR_53B02_AE_01_TrpA_Rig1Plate15BowlC_20110930T140347'
  'GMR_71G01_AE_01_TrpA_Rig2Plate14BowlC_20110707T154934'
  'GMR_94B10_AE_01_TrpA_Rig1Plate15BowlC_20111007T155325'
  'pBDPGAL4U_TrpA_Rig1Plate15BowlB_20110922T145928'
  'pBDPGAL4U_TrpA_Rig2Plate14BowlD_20110615T164545'};

groundtruth_experiment_names = {
  'pBDPGAL4U_TrpA_Rig1Plate10BowlB_20110609T091905'
  'pBDPGAL4U_TrpA_Rig1Plate15BowlD_20120210T152130'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlA_20120202T145723'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlC_20110811T133805'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlD_20111130T094443'};

train_expdirs = cellfun(@(x) fullfile(roottraindir,x),train_experiment_names,'UniformOutput',false);
groundtruth_expdirs = cellfun(@(x) fullfile(rootgroundtruthdir,x),groundtruth_experiment_names,'UniformOutput',false);

try
  [hfigs,errordata,figfilenames,outmatfilename] = ...
  Compute_GroundTruthError_vs_TrainingSetSize(behavior,configfilename,train_expdirs,groundtruth_expdirs,...
  'hfigs',[3,4]+CHASE*100,...
  'figpos',[95 550 1246 443],...
  'ifpcolor',ifpcolor,...
  'fpcolor',fpcolor,...
  'ifncolor',ifncolor,...
  'fncolor',fncolor,...
  'ierrorcolor',ierrorcolor,...
  'errorcolor',errorcolor,...
  'chancecolor',chancecolor,...
  'outfigdir',outfigdir);
catch ME
  warning('%s: %s',behavior,getReport(ME));
end

%% wing grooming

behavior = 'WingGrooming';
configfilename = '../perframe/params/JLabelParams_FlyBowl_WingGrooming.xml';

roottraindir = 'C:/Data/JAABA/FlyBowl';
rootgroundtruthdir = 'C:/Data/JAABA/groundtruth_pBDPGAL4U_data';

train_experiment_names = {
  'GMR_71G01_AE_01_TrpA_Rig2Plate17BowlD_20110921T084818'
  'pBDPGAL4U_TrpA_Rig1Plate15BowlB_20110922T145928'
  'pBDPGAL4U_TrpA_Rig2Plate14BowlD_20110615T164545'
  'pBDPGAL4U_TrpA_Rig2Plate14BowlD_20110707T085804'};

groundtruth_experiment_names = {
  'pBDPGAL4U_TrpA_Rig1Plate10BowlB_20110609T091905'
  'pBDPGAL4U_TrpA_Rig1Plate15BowlD_20120210T152130'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlA_20120202T145723'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlC_20110811T133805'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlD_20111130T094443'};

train_expdirs = cellfun(@(x) fullfile(roottraindir,x),train_experiment_names,'UniformOutput',false);
groundtruth_expdirs = cellfun(@(x) fullfile(rootgroundtruthdir,x),groundtruth_experiment_names,'UniformOutput',false);

try
  [hfigs,errordata,figfilenames,outmatfilename] = ...
  Compute_GroundTruthError_vs_TrainingSetSize(behavior,configfilename,train_expdirs,groundtruth_expdirs,...
  'hfigs',[3,4]+WINGGROOMING*100,...
  'figpos',[95 550 1246 443],...
  'ifpcolor',ifpcolor,...
  'fpcolor',fpcolor,...
  'ifncolor',ifncolor,...
  'fncolor',fncolor,...
  'ierrorcolor',ierrorcolor,...
  'errorcolor',errorcolor,...
  'chancecolor',chancecolor,...
  'outfigdir',outfigdir);
catch ME
  warning('%s: %s',behavior,getReport(ME));
end

%% jump

behavior = 'Jump';
configfilename = '../perframe/params/JLabelParams_FlyBowl_Jump.xml';

roottraindir = 'C:/Data/JAABA/FlyBowl';
rootgroundtruthdir = 'C:/Data/JAABA/groundtruth_pBDPGAL4U_data';

train_experiment_names = {
'GMR_14C08_AE_01_TrpA_Rig1Plate15BowlB_20110914T113113'
'GMR_20B07_AE_01_TrpA_Rig2Plate17BowlC_20120216T143940'
'GMR_71G01_AE_01_TrpA_Rig2Plate14BowlC_20110707T154934'
'GMR_82D11_AE_01_TrpA_Rig1Plate15BowlB_20111028T093039'
'pBDPGAL4U_TrpA_Rig1Plate10BowlA_20110610T153218'
'pBDPGAL4U_TrpA_Rig2Plate17BowlA_20110929T143440'
};

groundtruth_experiment_names = {
  'pBDPGAL4U_TrpA_Rig1Plate10BowlB_20110609T091905'
  'pBDPGAL4U_TrpA_Rig1Plate15BowlD_20120210T152130'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlA_20120202T145723'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlC_20110811T133805'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlD_20111130T094443'};

train_expdirs = cellfun(@(x) fullfile(roottraindir,x),train_experiment_names,'UniformOutput',false);
groundtruth_expdirs = cellfun(@(x) fullfile(rootgroundtruthdir,x),groundtruth_experiment_names,'UniformOutput',false);

try
  [hfigs,errordata,figfilenames,outmatfilename] = ...
  Compute_GroundTruthError_vs_TrainingSetSize(behavior,configfilename,train_expdirs,groundtruth_expdirs,...
  'hfigs',[3,4]+JUMP*100,...
  'figpos',[95 550 1246 443],...
  'ifpcolor',ifpcolor,...
  'fpcolor',fpcolor,...
  'ifncolor',ifncolor,...
  'fncolor',fncolor,...
  'ierrorcolor',ierrorcolor,...
  'errorcolor',errorcolor,...
  'chancecolor',chancecolor,...
  'outfigdir',outfigdir);
catch ME
  warning('%s: %s',behavior,getReport(ME));
end

%% stop

behavior = 'Stop';
configfilename = '../perframe/params/JLabelParams_FlyBowl_Stop.xml';

roottraindir = 'C:/Data/JAABA/FlyBowl';
rootgroundtruthdir = 'C:/Data/JAABA/groundtruth_pBDPGAL4U_data';

% find train directories by loading in the classifier
train_experiment_names = {
  'pBDPGAL4U_TrpA_Rig2Plate17BowlB_20110805T101647'
};

% always doing groundtruthing on these
groundtruth_experiment_names = {
  'pBDPGAL4U_TrpA_Rig1Plate10BowlB_20110609T091905'
  'pBDPGAL4U_TrpA_Rig1Plate15BowlD_20120210T152130'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlA_20120202T145723'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlC_20110811T133805'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlD_20111130T094443'};

train_expdirs = cellfun(@(x) fullfile(roottraindir,x),train_experiment_names,'UniformOutput',false);
groundtruth_expdirs = cellfun(@(x) fullfile(rootgroundtruthdir,x),groundtruth_experiment_names,'UniformOutput',false);

try
  [hfigs,errordata,figfilenames,outmatfilename] = ...
  Compute_GroundTruthError_vs_TrainingSetSize(behavior,configfilename,train_expdirs,groundtruth_expdirs,...
  'hfigs',[3,4]+STOP*100,...
  'figpos',[95 550 1246 443],...
  'ifpcolor',ifpcolor,...
  'fpcolor',fpcolor,...
  'ifncolor',ifncolor,...
  'fncolor',fncolor,...
  'ierrorcolor',ierrorcolor,...
  'errorcolor',errorcolor,...
  'chancecolor',chancecolor,...
  'outfigdir',outfigdir);
catch ME
  warning('%s: %s',behavior,getReport(ME));
end

%% righting

behavior = 'Righting';
configfilename = '../perframe/params/JLabelParams_FlyBowl_Righting.xml';

roottraindir = 'C:/Data/JAABA/FlyBowl';
rootgroundtruthdir = 'C:/Data/JAABA/groundtruth_pBDPGAL4U_data';

% find train directories by loading in the classifier
train_experiment_names = {
  'pBDPGAL4U_TrpA_Rig1Plate10BowlA_20110610T153218'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlA_20110929T143440'
  'GMR_14C08_AE_01_TrpA_Rig1Plate15BowlB_20110914T113113'
  'GMR_71G01_AE_01_TrpA_Rig2Plate14BowlC_20110707T154934'
  'GMR_82D11_AE_01_TrpA_Rig1Plate15BowlB_20111028T093039'
};

% always doing groundtruthing on these
groundtruth_experiment_names = {
  'pBDPGAL4U_TrpA_Rig1Plate10BowlB_20110609T091905'
  'pBDPGAL4U_TrpA_Rig1Plate15BowlD_20120210T152130'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlA_20120202T145723'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlC_20110811T133805'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlD_20111130T094443'};

train_expdirs = cellfun(@(x) fullfile(roottraindir,x),train_experiment_names,'UniformOutput',false);
groundtruth_expdirs = cellfun(@(x) fullfile(rootgroundtruthdir,x),groundtruth_experiment_names,'UniformOutput',false);

try
  [hfigs,errordata,figfilenames,outmatfilename] = ...
  Compute_GroundTruthError_vs_TrainingSetSize(behavior,configfilename,train_expdirs,groundtruth_expdirs,...
  'hfigs',[3,4]+RIGHTING*100,...
  'figpos',[95 550 1246 443],...
  'ifpcolor',ifpcolor,...
  'fpcolor',fpcolor,...
  'ifncolor',ifncolor,...
  'fncolor',fncolor,...
  'ierrorcolor',ierrorcolor,...
  'errorcolor',errorcolor,...
  'chancecolor',chancecolor,...
  'outfigdir',outfigdir);
catch ME
  warning('%s: %s',behavior,getReport(ME));
end

%% tail pivot turn

behavior = 'TailPivotTurn';
configfilename = '../perframe/params/JLabelParams_FlyBowl_TailPivotTurn.xml';

roottraindir = 'C:/Data/JAABA/FlyBowl';
rootgroundtruthdir = 'C:/Data/JAABA/groundtruth_pBDPGAL4U_data';

% find train directories by loading in the classifier
train_experiment_names = {
  'GMR_82D11_AE_01_TrpA_Rig1Plate15BowlB_20111028T093039'
  'GMR_80A01_AE_01_TrpA_Rig2Plate14BowlD_20110408T140618'
  'pBDPGAL4U_TrpA_Rig1Plate15BowlB_20110922T145928'
  'pBDPGAL4U_TrpA_Rig2Plate14BowlD_20110615T164545'
  'FCF_attP2_1500062_None_Rig2Plate17BowlB_20120519T165213'
};

% always doing groundtruthing on these
groundtruth_experiment_names = {
  'pBDPGAL4U_TrpA_Rig1Plate10BowlB_20110609T091905'
  'pBDPGAL4U_TrpA_Rig1Plate15BowlD_20120210T152130'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlA_20120202T145723'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlC_20110811T133805'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlD_20111130T094443'};

train_expdirs = cellfun(@(x) fullfile(roottraindir,x),train_experiment_names,'UniformOutput',false);
groundtruth_expdirs = cellfun(@(x) fullfile(rootgroundtruthdir,x),groundtruth_experiment_names,'UniformOutput',false);

try
  [hfigs,errordata,figfilenames,outmatfilename] = ...
  Compute_GroundTruthError_vs_TrainingSetSize(behavior,configfilename,train_expdirs,groundtruth_expdirs,...
  'hfigs',[3,4]+TAILPIVOTTURN*100,...
  'figpos',[95 550 1246 443],...
  'ifpcolor',ifpcolor,...
  'fpcolor',fpcolor,...
  'ifncolor',ifncolor,...
  'fncolor',fncolor,...
  'ierrorcolor',ierrorcolor,...
  'errorcolor',errorcolor,...
  'chancecolor',chancecolor,...
  'outfigdir',outfigdir);
catch ME
  warning('%s: %s',behavior,getReport(ME));
end

%% walk

behavior = 'Walk';
configfilename = '../perframe/params/JLabelParams_FlyBowl_Walk.xml';

roottraindir = 'C:/Data/JAABA/FlyBowl';
rootgroundtruthdir = 'C:/Data/JAABA/groundtruth_pBDPGAL4U_data';

% find train directories by loading in the classifier
train_experiment_names = {
 'pBDPGAL4U_TrpA_Rig2Plate17BowlB_20110805T101647'
 'pBDPGAL4U_TrpA_Rig2Plate14BowlD_20110707T085804'
 'pBDPGAL4U_TrpA_Rig2Plate14BowlD_20110615T164545'
 'pBDPGAL4U_TrpA_Rig1Plate15BowlB_20110922T145928'
 'GMR_80A01_AE_01_TrpA_Rig2Plate14BowlD_20110408T140618'
 'GMR_82D11_AE_01_TrpA_Rig1Plate15BowlB_20111028T093039'
};

% always doing groundtruthing on these
groundtruth_experiment_names = {
  'pBDPGAL4U_TrpA_Rig1Plate10BowlB_20110609T091905'
  'pBDPGAL4U_TrpA_Rig1Plate15BowlD_20120210T152130'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlA_20120202T145723'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlC_20110811T133805'
  'pBDPGAL4U_TrpA_Rig2Plate17BowlD_20111130T094443'};

train_expdirs = cellfun(@(x) fullfile(roottraindir,x),train_experiment_names,'UniformOutput',false);
groundtruth_expdirs = cellfun(@(x) fullfile(rootgroundtruthdir,x),groundtruth_experiment_names,'UniformOutput',false);

try
  [hfigs,errordata,figfilenames,outmatfilename] = ...
  Compute_GroundTruthError_vs_TrainingSetSize(behavior,configfilename,train_expdirs,groundtruth_expdirs,...
  'hfigs',[3,4]+WALK*100,...
  'figpos',[95 550 1246 443],...
  'ifpcolor',ifpcolor,...
  'fpcolor',fpcolor,...
  'ifncolor',ifncolor,...
  'fncolor',fncolor,...
  'ierrorcolor',ierrorcolor,...
  'errorcolor',errorcolor,...
  'chancecolor',chancecolor,...
  'outfigdir',outfigdir);
catch ME
  warning('%s: %s',behavior,getReport(ME));
end

author: kristinbranson

repo: JAABA


text:
function [optarg] = ft_cfg2keyval(cfg)

% FT_CFG2KEYVAL converts between a structure and a cell-array with key-value
% pairs which can be used for optional input arguments.
%
% Use as
%   optarg = ft_cfg2keyval(cfg)
%
% See also FT_KEYVAL2CFG, FT_GETOPT

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

if ~isempty(cfg)
  optarg = [fieldnames(cfg) struct2cell(cfg)]';
  optarg = optarg(:)';
else
  optarg = {};
end


text:
function [varargout] = myProcessOptions(options,varargin)
% Similar to processOptions, but case insensitive and
%   using a struct instead of a variable length list

options = toUpper(options);

for i = 1:2:length(varargin)
    if isfield(options,upper(varargin{i}))
        v = getfield(options,upper(varargin{i}));
        if isempty(v)
            varargout{(i+1)/2}=varargin{i+1};
        else
            varargout{(i+1)/2}=v;
        end
    else
        varargout{(i+1)/2}=varargin{i+1};
    end
end

end

function [o] = toUpper(o)
if ~isempty(o)
    fn = fieldnames(o);
    for i = 1:length(fn)
        o = setfield(o,upper(fn{i}),getfield(o,fn{i}));
    end
end
end

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function [V,F,UV,C,N] = readOFF( filename )
  % READOFF reads an OFF file with vertex/face information
  %
  % [V,F,UV,C,N] = readOFF( filename )
  %
  % Input:
  %  filename  path to .off file
  % Outputs:
  %  V  #V by 3 list of vertices
  %  F  #F by 3 list of triangle indices
  %  UV  #V by 2 list of texture coordinates
  %  C  #V by 3 list of colors
  %  N  #V by 3 list of normals
  %
  % See also: load_mesh, readOBJfast, readOBJ

% (C) 2007 Denis Kovacs, NYU
%-------------------------------------------------------------------------

  V = [];
  F = [];
  UV = [];
  C = [];
  N = [];

  fp = fopen( filename, 'r' );
  OFFheader = upper(fscanf( fp, '%s\n', 1 ));
  if OFFheader(end-2:end) ~= 'OFF'
    warning('no OFF file!')
    fclose(fp);
    return;
  end
  OFFdim = 3;
  OFF_N = 0; OFF_C=0; OFF_ST=0;

  if find(OFFheader=='N') OFFdim = OFFdim+3; OFF_N=1; end
  if find(OFFheader=='C') OFFdim = OFFdim+3; OFF_C=1; end
  if find(OFFheader=='S') OFFdim = OFFdim+2; OFF_ST=1; end

  % eat any comments before
  line = eat_comments(fp,'#');

  d = sscanf( line, '%d', 3);
  nV = d(1); nF = d(2); nE = d(3);

  %disp(sprintf('  - Reading %d vertices', nV));

  switch OFFdim
      case  3; OFFV = textscan( fp, '%f %f %f', nV);
      case  5; OFFV = textscan( fp, '%f %f %f %f %f', nV);
      case  6; OFFV = textscan( fp, '%f %f %f %f %f %f', nV);
      case  7; OFFV = textscan( fp, '%f %f %f %f %f %f %f', nV);
      case  8; OFFV = textscan( fp, '%f %f %f %f %f %f %f %f', nV);
      case  9; OFFV = textscan( fp, '%f %f %f %f %f %f %f %f %f', nV);
      case 10; OFFV = textscan( fp, '%f %f %f %f %f %f %f %f %f %f', nV);
      case 11; OFFV = textscan( fp, '%f %f %f %f %f %f %f %f %f %f %f', nV);
      otherwise; fclose(fp); error('Unsupported number of vertex entries');
  end

  try
     OFFV = cell2mat(OFFV);
  end

  OFFdim = 1;
  V = OFFV(:,OFFdim:(OFFdim+2)); OFFdim = OFFdim + 3;
  if (OFF_N) N = OFFV(:,OFFdim:(OFFdim+2)); OFFdim = OFFdim + 3; end
  if (OFF_C) C = OFFV(:,OFFdim:(OFFdim+2)); OFFdim = OFFdim + 3; end
  if (OFF_ST) UV = OFFV(:,OFFdim:(OFFdim+1)); OFFdim = OFFdim + 2; end

  if (nF ~= 0)
    %disp(sprintf('  - Reading %d faces', nF));
    temp = textscan( fp, '%d %d %d %d %d %d %d %d %d %d %d', nF );
    sz = temp{1}(1);
    if all(sz == cell2mat(temp(1)))
      F = double (cell2mat( temp(2:sz+1 ))) +1;
    else
      warning('Trivially triangulating high degree facets');
      F = zeros(sum(temp{1}-2),3);
      fi = 1;
      for f = 1:size(temp{1},1)
        for j = 3:temp{1}(f);
          F(fi,:) = [temp{2}(f) temp{j}(f) temp{j+1}(f)]+1;
          fi = fi+1;
        end
      end
    end
  else
    F = [];
  end

  fclose(fp);

  %disp('  - done.');
end

author: alecjacobson

repo: gptoolbox


text:
function epoch = mydatestdi (epoch_std, year)
    [year_len, year_start] = mydatestd_aux(year);

    std_year_len = mydatestd_aux();

    epoch = epoch_std ./ std_year_len .* year_len + year_start;

    if any(epoch_std > 2.*std_year_len) ...
    || any(epoch_std <   -std_year_len)
        error('MATLAB:mydatestdi:outRange', ...
        'First argument is out of valid input range.');
    end

    year = repmat(year, length(epoch_std), 1);

    idx = (epoch_std > std_year_len);
    if any(idx)
        epoch(idx) = mydatestdi(epoch_std(idx)-std_year_len, year(idx)+1);
    end

    idx = (epoch_std < 0);
    if any(idx)
        epoch(idx) = mydatestdi(std_year_len-abs(epoch_std(idx)), year(idx)-1);
    end
end

%!test
%! % mydatestdi()
%! test('mydatestd_aux')


text:
function x = parametric(x)
%PARAMETRIC Defines a variable as parametric
%
%   F = PARAMETRIC(x) is used to describe the set of parametric variables
%   in a multi-parametric program, as an alternative to using the 4th input
%   in SOLVEMP
%
%   It can also be used to define the parametric variables in a SOS program
%
%   INPUT
%    x : SDPVAR object
%
%   OUTPUT
%    F : Constraint
%
%   EXAMPLE
%    F = parametric(x);               % Short notation
%
%   See also SOLVEMP, SOLVESOS, SDPVAR

x.typeflag = 13;
x = lmi(x);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [hFig, iDS, iFig] = view_channels(ChannelFile, Modality, isMarkers, isLabels, hFig, is3DElectrodes, selChannels)
% VIEW_CHANNELS: Display a channel file in all the associated 3DViz figure.
%
% USAGE: view_channels(ChannelFile, Modality, isMarker=1, isLabels=1, hFig=[], is3DElectrodes=0, selChannels=[])
%        view_channels(ChannelFile, Modality)
% OUTPUT:
%     - hFig : Matlab handle to the 3DViz figure that was created or updated
%     - iDS  : DataSet index in the GlobalData variable
%     - iFig : Indice of returned figure in the GlobalData(iDS).Figure array

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

global GlobalData;
% Parse inputs
if (nargin < 7) || isempty(selChannels)
    selChannels = [];
end
if (nargin < 6) || isempty(is3DElectrodes)
    is3DElectrodes = 0;
end
if (nargin < 5) || isempty(hFig)
    hFig = [];
end
if (nargin < 4) || isempty(isLabels) || isempty(isMarkers)
    isMarkers = 1;
    isLabels  = 1;
end
iDS  = [];
iFig = [];
% Get short path to this file
ChannelFile = file_short(ChannelFile);

% ===== GET/CREATE DATASET =====
if isempty(hFig)
    % Get Study that holds this ChannelFile
    sStudy = bst_get('ChannelFile', ChannelFile);
    % If this surface does not belong to any subject
    if isempty(sStudy)
        StudyFile = '';
        SubjectFile = '';
        %warning('This surface does not belong to any database subject.');
        % Check that the SurfaceFile really exist as an absolute file path
        if ~file_exist(ChannelFile)
            bst_error(['File not found : "', ChannelFile, '"'], 'Display surface');
            return
        end
    else
        StudyFile   = sStudy.FileName;
        SubjectFile = sStudy.BrainStormSubject;
        iDS = [];
        % Get GlobalData DataSet associated with channel file
        iDSChannel = bst_memory('GetDataSetChannel', ChannelFile);
        if ~isempty(iDSChannel) && ~isempty([GlobalData.DataSet(iDSChannel).Figure])
            iDS = iDSChannel;
        end
        iDSStudy   = [];
        iDSSubject = [];
        % ChannelFile not found in DataSets (or DataSet has no figure), try StudyFile
        if isempty(iDS)
            iDSStudy = bst_memory('GetDataSetStudy', StudyFile);
            if ~isempty(iDSStudy) && ~isempty([GlobalData.DataSet(iDSStudy).Figure])
                iDS = iDSStudy;
            end
        end
        % StudyFile not found in DataSets, try SubjectFile
        if isempty(iDS)
            iDSSubject = bst_memory('GetDataSetSubject', SubjectFile, 0);
            % Do not accept DataSet if an other ChannelFile is already attributed to the DataSet
            iDSSubject = iDSSubject(cellfun(@(c)isempty(c), {GlobalData.DataSet(iDSSubject).ChannelFile}));
            % Check if Datasets are acceptable
            if ~isempty(iDSSubject) && ~isempty([GlobalData.DataSet(iDSSubject).Figure])
                iDS = iDSSubject(1);
            end
        end
        % If no dataset with figure was selected, get the first dataset without figures
        if isempty(iDS)
            iDS = [iDSChannel iDSStudy iDSSubject];
            if ~isempty(iDS)
                iDS = iDS(1);
            end
        end
    end
    % If no existing dataset is found : create a new one
    if isempty(iDS)
       iDS = bst_memory('GetDataSetEmpty');
    end
    iDS = iDS(1);
    GlobalData.DataSet(iDS).SubjectFile = SubjectFile;
    GlobalData.DataSet(iDS).StudyFile   = StudyFile;

% ===== RE-USE EXISTING FIGURE =====
else
    % Get figure definition
    [hFig, iFig, iDS] = bst_figures('GetFigure', hFig);
end

% ===== LOAD CHANNEL FILE IF NOT DONE =====
% Also check if we are processing the same channel file, which can be
% different when doing real-time head position alignment (2 helmets from
% different datasets in the same figure).
if isempty(GlobalData.DataSet(iDS).ChannelFile) || ~strcmpi(GlobalData.DataSet(iDS).ChannelFile, ChannelFile) || ...
   (~GlobalData.DataSet(iDS).isChannelModified && isempty(GlobalData.DataSet(iDS).Channel))
    % Check if this channel file is already loaded and modified in another DataSet
    iDSother = setdiff(1:length(GlobalData.DataSet), iDS);
    if ~isempty(iDSother) && any([GlobalData.DataSet(iDSother).isChannelModified])
        % Ask user
        isSave = java_dialog('confirm', ...
            ['This channel file is being edited in another window.' 10 ...
             'Save the modifications so the new figure can show updated positions?'], 'Save modifications');
        % Force saving of the modifications
        if isSave
            bst_memory('SaveChannelFile', iDSother(1));
        end
    end
    % Load channel file
    ChannelMat = in_bst_channel(ChannelFile);
    % Copy information in memory
    GlobalData.DataSet(iDS).ChannelFile     = ChannelFile;
    GlobalData.DataSet(iDS).Channel         = ChannelMat.Channel;
    GlobalData.DataSet(iDS).MegRefCoef      = ChannelMat.MegRefCoef;
    GlobalData.DataSet(iDS).Projector       = ChannelMat.Projector;
    GlobalData.DataSet(iDS).Clusters        = ChannelMat.Clusters;
    GlobalData.DataSet(iDS).IntraElectrodes = ChannelMat.IntraElectrodes;
    GlobalData.DataSet(iDS).HeadPoints      = ChannelMat.HeadPoints;
end

% ===== CREATE 3DVIZ FIGURE =====
% Progress bar
isProgress = ~bst_progress('isVisible');
if isProgress
    bst_progress('start', 'View sensors', 'Loading data...');
end
if isempty(hFig)
    % Prepare FigureId structure
    FigureId = db_template('FigureId');
    FigureId.Type     = '3DViz';
    FigureId.SubType  = '';
    FigureId.Modality = Modality;
    % Create figure
    [hFig, iFig, isNewFig] = bst_figures('CreateFigure', iDS, FigureId);
    % If figure was not created: Display an error message and return
    if isempty(hFig)
        bst_error('Cannot create figure', '3D figure creation...', 0);
        return;
    end
    setappdata(hFig, 'StudyFile',    StudyFile);
    setappdata(hFig, 'SubjectFile',  SubjectFile);
else
    isNewFig = 0;
    % Get existing list of selected channels
    if ~isempty(GlobalData.DataSet(iDS).Figure(iFig).SelectedChannels)
        selChannels = GlobalData.DataSet(iDS).Figure(iFig).SelectedChannels;
    end
end
% Make sure that the Modality is saved
GlobalData.DataSet(iDS).Figure(iFig).Id.Modality = Modality;
% Set application data
setappdata(hFig, 'DataFile', '');
setappdata(hFig, 'isSensorsOnly', 1);

% ===== DISPLAY SENSORS =====
% Update figure selection
bst_figures('SetCurrentFigure', hFig, '3D');
% Get default selected channels for this modality
if isempty(selChannels)
    % Get selected channels
    selChannels = good_channel(GlobalData.DataSet(iDS).Channel, [], Modality);
    % Remove the channels with empty positions
    iNoLoc = find(cellfun(@isempty, {GlobalData.DataSet(iDS).Channel.Loc}));
    if ~isempty(iNoLoc)
        selChannels = setdiff(selChannels, iNoLoc);
    end
    if isempty(selChannels) && ~isempty(iNoLoc)
        error('None of the channels have positions defined.');
    end
end
% Set figure selected sensors
GlobalData.DataSet(iDS).Figure(iFig).SelectedChannels = selChannels;
% Display sensors
if is3DElectrodes
    figure_3d('PlotSensors3D', iDS, iFig);
    if isLabels
        figure_3d('ViewSensors', hFig, 1, 1, 0, Modality);
    end
    % Open iEEG tab
    if ismember(Modality, {'SEEG', 'ECOG', 'ECOG+SEEG'})
        gui_brainstorm('ShowToolTab', 'iEEG');
    end
elseif ismember(lower(Modality), {'ctf', 'vectorview306', '4d', 'kit', 'kriss', 'babymeg', 'ricoh'})
    figure_3d('PlotCoils', hFig, Modality, isMarkers);
elseif ismember(lower(Modality), {'nirs','nirs-brs'})
    figure_3d('PlotNirsCap', hFig, isMarkers);
else
    isMesh = ~isequal(Modality, 'SEEG');
    figure_3d('ViewSensors', hFig, isMarkers, isLabels, isMesh, Modality);
end
% Open Cluster tab
if ~isempty(GlobalData.DataSet(iDS).Clusters)
    gui_brainstorm('ShowToolTab', 'Cluster');
end
% Update lights
camlight(findobj(hFig, 'Tag', 'FrontLight'), 'headlight');
% Update figure name
bst_figures('UpdateFigureName', hFig);
% Camera basic orientation
if isNewFig
    figure_3d('SetStandardView', hFig, 'top');
end
% Show figure
set(hFig, 'Visible', 'on');
% Close progress bar
if isProgress
    bst_progress('stop');
end

author: brainstorm-tools