text:
% CONFIGSTATEMACHINE configures the state machine (type of demo, velocity
%                    of the demo, repeat movements, and so on).

%% --- Initialization ---

% If true, the robot CoM will follow a desired reference trajectory (COORDINATOR DEMO ONLY)
Config.LEFT_RIGHT_MOVEMENTS = false;

% If equal to one, the desired values of the center of mass are smoothed internally
Config.SMOOTH_COM_DES       = true;

% If equal to one, the desired values of the postural tasks are smoothed internally
Config.SMOOTH_JOINT_DES     = true;

% Joint torques saturation [Nm]
Sat.torque                  = 60;

% Joint torques rate of change saturation
Sat.uDotMax                 = 300;

% max unsigned difference between two consecutive (measured) joint positions,
% i.e. delta_qj = abs(qj(k) - qj(k-1))
Sat.maxJointsPositionDelta  = 15*pi/180; % [rad]

%% Regularization parameters
Reg.pinvDamp_baseVel        = 1e-7;
Reg.pinvDamp                = 1;
Reg.pinvTol                 = 1e-5;
Reg.KP_postural             = 0.1;
Reg.KD_postural             = 0;
Reg.HessianQP               = 1e-7;

%% State Machine configuration

% time between two yoga positions
StateMachine.joints_pauseBetweenYogaMoves = 3;

% contact forces threshold
StateMachine.wrench_thresholdContactOn    = 25;
StateMachine.wrench_thresholdContactOff   = 85;

% threshold on CoM and joints error
StateMachine.CoM_threshold                = 0.01;
StateMachine.joints_thresholdNotInContact = 5;
StateMachine.joints_thresholdInContact    = 50;

% initial state for state machine
StateMachine.initialState                 = 1;

% other configuration parameters for state machine
StateMachine.tBalancing                   = 1;
StateMachine.tBalancingBeforeYoga         = 1;
StateMachine.skipYoga                     = false;
StateMachine.demoOnlyBalancing            = false;
StateMachine.demoStartsOnRightSupport     = false; % If false, the Yoga demo is performed on the left foot first
StateMachine.yogaAlsoOnRightFoot          = false; % TO DO: yoga on both feet starting from right foot (not available for now)

%%%% List of possible "Yoga in loop" %%%%

% the robot will repeat the FULL DEMO (two feet balancing, yoga on left
% foot, back on two feet, yoga right foot, back on two feet). The demo is
% repeated until the user stops the Simulink model. This option is ignored
% if Sm.demoStartsOnRightSupport = true.
StateMachine.twoFeetYogaInLoop            = false;

% the robot will repeat the ONE FOOT yoga for the number of times the user
% specifies in the Sm.yogaCounter option. The robot WILL NOT go back to two
% feet balancing in between to consecutive yoga. WARNING: if the option
% Sm.yogaAlsoOnRightFoot is true, then the robot will repeat first the yoga
% on left foot for the number of times the user specifies in the Sm.yogaCounter,
% and then it will repeat the yoga on the right foot for the same number of times.
StateMachine.oneFootYogaInLoop            = false;
StateMachine.yogaCounter                  = 5;

%% Parameters for motors reflected inertia

% transmission ratio (1/N)
Config.Gamma                 = 0.01*eye(ROBOT_DOF);

% modify the value of the transmission ratio for the hip pitch.
% TODO: avoid to hard-code the joint numbering
Config.Gamma(end-5, end-5)   = 0.0067;
Config.Gamma(end-11,end-11)  = 0.0067;

% motors inertia (Kg*m^2)
legsMotors_I_m               = 0.0827*1e-4;
torsoPitchRollMotors_I_m     = 0.0827*1e-4;
torsoYawMotors_I_m           = 0.0585*1e-4;
armsMotors_I_m               = 0.0585*1e-4;

% add harmonic drives reflected inertia
if Config.INCLUDE_HARMONIC_DRIVE_INERTIA

    legsMotors_I_m           = legsMotors_I_m + 0.054*1e-4;
    torsoPitchRollMotors_I_m = torsoPitchRollMotors_I_m + 0.054*1e-4;
    torsoYawMotors_I_m       = torsoYawMotors_I_m + 0.021*1e-4;
    armsMotors_I_m           = armsMotors_I_m + 0.021*1e-4;
end

Config.I_m                   = diag([torsoPitchRollMotors_I_m*ones(2,1);
                                     torsoYawMotors_I_m;
                                     armsMotors_I_m*ones(8,1);
                                     legsMotors_I_m*ones(12,1)]);

% parameters for coupling matrices. Updated according to the wiki:
%
% http://wiki.icub.org/wiki/ICub_coupled_joints
%
% and corrected according to https://github.com/robotology/robots-configuration/issues/39
t            = 0.615;
r            = 0.022;
R            = 0.04;

% coupling matrices
T_LShoulder  = [-1  0  0;
                -1 -t  0;
                 0  t -t];

T_RShoulder  = [ 1  0  0;
                 1  t  0;
                 0 -t  t];

T_torso      = [ 0.5    -0.5     0;
                 0.5     0.5     0;
                 r/(2*R) r/(2*R) r/R];

if Config.INCLUDE_COUPLING

    Config.T = blkdiag(T_torso,T_LShoulder,1,T_RShoulder,1,eye(12));
else
    Config.T = eye(ROBOT_DOF);
end

% gain for feedforward term in joint torques calculation. Valid range: a
% value between 0 and 1
Config.K_ff  = 0;

% Config.USE_DES_JOINT_ACC_FOR_MOTORS_INERTIA if true, the desired joints
% accelerations are used for computing the feedforward term in joint
% torques calculations. Not effective if Config.K_ff = 0.
Config.USE_DES_JOINT_ACC_FOR_MOTORS_INERTIA = false;

%% Constraints for QP for balancing

% The friction cone is approximated by using linear interpolation of the circle.
% So, numberOfPoints defines the number of points used to interpolate the circle
% in each cicle's quadrant
numberOfPoints               = 4;
forceFrictionCoefficient     = 1/3;
torsionalFrictionCoefficient = 1/75;
fZmin                        = 10;

% physical size of the foot
feet_size                    = [-0.05  0.10;     % xMin, xMax
                                -0.025 0.025];   % yMin, yMax

% Compute contact constraints (friction cone, unilateral constraints)
[ConstraintsMatrix, bVectorConstraints] = wbc.computeRigidContactConstraints ...
    (forceFrictionCoefficient, numberOfPoints, torsionalFrictionCoefficient, feet_size, fZmin);

author: robotology

repo: whole-body-controllers


text:
function days = month_length_common ( y, m )

%*****************************************************************************80
%
%% MONTH_LENGTH_COMMON returns the number of days in a Common month.
%
%  Discussion:
%
%    The "common" calendar is meant to be the calendar which is Julian up to
%    day JED = 2299160, and Gregorian from day JED = 2299161 and after.
%
%    The routine knows that February has 28 days, except in leap years,
%    when it has 29.
%
%    In the Common calendar, October 1582 had only 21 days
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year in which the month occurred.
%
%    Input, integer M, the number of the month.
%
%    Output, integer DAYS, the number of days
%    in the month.
%
  mdays = [ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 ];
%
%  Check the input.
%
  [ y2, m2, ierror ] = ym_check_common ( y, m );

  if ( ierror ~= 0 )
    days = 0;
    return
  end
%
%  Take care of the special case.
%
  if ( y2 == 1582 )
    if ( m2 == 10 )
      days = 21;
      return
    end
  end
%
%  Get the number of days in the month.
%
  days = mdays(m2);
%
%  If necessary, add 1 day for February 29.
%
  if ( m2 == 2 && year_is_leap_common ( y2 ) )
    days = days + 1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function hermite_polynomial_test09 ( p, e )

%*****************************************************************************80
%
%% HERMITE_POLYNOMIAL_TEST09 tests HN_POWER_PRODUCT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer P, the maximum degree of the polynomial factors.
%
%    Input, integer E, the exponent of X in the integrand.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_POLYNOMIAL_TEST09\n' );
  fprintf ( 1, '  Compute a normalized physicist''s Hermite polynomial power product table.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Tij = integral ( -oo < X < +oo ) X^E Hn(I,X) Hn(J,X) exp(-X*X) dx\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  where Hn(I,X) = normalized physicist''s Hermite polynomial of degree I.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Maximum degree P = %d\n', p );
  fprintf ( 1, '  Exponent of X = %d\n', e );

  table = hn_power_product ( p, e );

  r8mat_print ( p + 1, p + 1, table, '  Power product table:' );

  return
end


text:
function centroid = polyhedron_centroid_3d ( coord, maxorder, face_num, node, ...
  node_num, order )

%*****************************************************************************80
%
%% POLYHEDRON_CENTROID_3D computes the centroid of a polyhedron in 3D.
%
%  Discussion:
%
%    The centroid can be computed as the volume-weighted average of
%    the centroids of the tetrahedra defined by choosing a point in
%    the interior of the polyhedron, and using as a base every triangle
%    created by triangulating the faces of the polyhedron.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real COORD(3,NODE_NUM), the vertices.
%    The vertices may be listed in any order.
%
%    Input, integer MAXORDER, the maximum number of vertices that make
%    up a face of the polyhedron.
%
%    Input, integer FACE_NUM, the number of faces of the polyhedron.
%
%    Input, integer NODE(FACE_NUM,MAXORDER).  Face I is defined by
%    the vertices NODE(I,1) through NODE(I,ORDER(I)).  These vertices
%    are listed in neighboring order.
%
%    Input, integer NODE_NUM, the number of points stored in COORD.
%
%    Input, integer ORDER(FACE_NUM), the number of vertices making up
%    each face.
%
%    Output, real CENTROID(3), the centroid of the polyhedron.
%
  dim_num = 3;
%
%  Compute a point in the interior.
%  We take the area-weighted centroid of each face.
%
  point(1:dim_num) = 0.0;
  area = 0.0;

  for face = 1 : face_num

    vert_num = order(face);

    v(1:dim_num,1:vert_num) = coord(1:dim_num,node(face,1:vert_num));

    [ polygon_area, normal ] = polygon_area_3d ( vert_num, v );

    polygon_centroid(1:dim_num) = polygon_centroid_3d ( vert_num, v );

    point(1:dim_num) = point(1:dim_num) + polygon_area * polygon_centroid(1:dim_num);

    area = area + polygon_area;

  end

  point(1:dim_num) = point(1:dim_num) / area;
%
%  Now triangulate each face.
%  For each triangle, consider the tetrahedron created by including POINT.
%
  centroid(1:dim_num) = 0.0;
  volume = 0.0;

  for face = 1 : face_num

    n3 = node(face,order(face));

    for vert = 1 : order(face) - 2

      n1 = node(face,vert);
      n2 = node(face,vert+1);

      tetra(1:dim_num,1:4) = [ ...
        coord(1:dim_num,n1)'; coord(1:dim_num,n2)'; coord(1:dim_num,n3)'; point(1:dim_num) ]';

      tetra_volume = tetrahedron_volume_3d ( tetra );

      tetra_centroid(1:dim_num) = tetrahedron_centroid_3d ( tetra );

      centroid(1:dim_num) = centroid(1:dim_num) + tetra_volume * tetra_centroid(1:dim_num);

      volume = volume + tetra_volume;

    end
  end

  centroid(1:dim_num) = centroid(1:dim_num) / volume;

  return
end


text:
% [PYR, INDICES, STEERMTX, HARMONICS] = buildSpyr(IM, HEIGHT, FILTFILE, EDGES)
%
% Construct a steerable pyramid on matrix IM.
%
% HEIGHT (optional) specifies the number of pyramid levels to build. Default
% is maxPyrHt(size(IM),size(FILT));
%
% FILTFILE (optional) should be a string referring to an m-file that
% returns the rfilters.  (examples: 'sp0Filters', 'sp1Filters',
% 'sp3Filters','sp5Filters'.  default = 'sp1Filters'). EDGES specifies
% edge-handling, and defaults to 'reflect1' (see corrDn).
%
% PYR is a vector containing the N pyramid subbands, ordered from fine
% to coarse.  INDICES is an Nx2 matrix containing the sizes of
% each subband.  This is compatible with the MatLab Wavelet toolbox.
% See the function STEER for a description of STEERMTX and HARMONICS.

% Eero Simoncelli, 6/96.

function [pyr,pind,steermtx,harmonics] = buildSpyr(im, ht, filtfile, edges)

%-----------------------------------------------------------------
%% DEFAULTS:

if (exist('filtfile') ~= 1)
  filtfile = 'sp1Filters';
end

if (exist('edges') ~= 1)
  edges= 'reflect1';
end

if (isstr(filtfile) & (exist(filtfile) == 2))
   [lo0filt,hi0filt,lofilt,bfilts,steermtx,harmonics] = eval(filtfile);
else
  fprintf(1,'\nUse buildSFpyr for pyramids with arbitrary numbers of orientation bands.\n');
  error('FILTFILE argument must be the name of an M-file containing SPYR filters.');
end

max_ht = maxPyrHt(size(im), size(lofilt,1));
if (exist('ht') ~= 1)
  ht = max_ht;
else
  if (ht > max_ht)
    error(sprintf('Cannot build pyramid higher than %d levels.',max_ht));
  end
end

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

hi0 = corrDn(im, hi0filt, edges);
lo0 = corrDn(im, lo0filt, edges);

[pyr,pind] = buildSpyrLevs(lo0, ht, lofilt, bfilts, edges);

pyr = [hi0(:) ; pyr];
pind = [size(hi0); pind];

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
cur_dir = pwd;
cd(fileparts(mfilename('fullpath')));

try
    fprintf('Downloading model_VGG16...\n');
    urlwrite('https://onedrive.live.com/download?resid=36FEC490FBC32F1A!114&authkey=!AE8uV9B07dREbhM&ithint=file%2czip', ...
        'model_VGG16.zip');

    fprintf('Unzipping...\n');
    unzip('model_VGG16.zip', '..');

    fprintf('Done.\n');
    delete('model_VGG16.zip');
catch
    fprintf('Error in downloading, please try links in README.md https://github.com/ShaoqingRen/faster_rcnn');
end

cd(cur_dir);

author: ShaoqingRen


text:
function [update_flag, latest_version, latest_release_date, local_version] = check_update(local_version)

    if nargin < 1
        local_version = nmflibrary_version();
    end

    update_flag = 0;

    site_url = 'https://raw.githubusercontent.com/hiroyuki-kasai/NMFLibrary/master/auxiliary/nmflibrary_version.m';
    filename_full_local = 'dl_nmflibrary_version.m';

    connection_available = check_network_connection();
    if connection_available
        websave(filename_full_local, site_url);

        [latest_version, latest_release_date] = dl_nmflibrary_version();

        if str2double(latest_version) > str2double(local_version)
            update_flag = 1;
            %fprintf('Latest: %s, Local: %s\n', latest_version, local_version);
        else
            %fprintf('This version is latest (%s)\n', local_version);
        end

        delete dl_nmflibrary_version.m
    else
        update_flag = -1;
        latest_version = [];
        latest_release_date = [];
    end

end

function available = check_network_connection()
    available = false;
    if ispc
        [~, b] = system('ping -n 1 www.google.com');
    elseif isunix
        [~, b] = system('ping -c 1 www.google.com');
    else
        error('how did you even get Matlab to work?')
    end

    n = strfind(b, 'cannot resolve');
    if n > 0
    else
        available = true;
    end
end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [Ain,Cin,bin,fin] = sparse_NMF_initialization(Y,K,options) %beta,eta,X0,err_thr,max_iter)

T = size(Y,ndims(Y));
if ~ismatrix(Y)
    Y = reshape(Y,numel(Y)/T,T);
end

defoptions = CNMFSetParms;
if nargin < 3 || isempty(options); options = defoptions; end
if ~isfield(options,'snmf_max_iter'); options.snmf_max_iter = defoptions.snmf_max_iter; end
    max_iter = options.snmf_max_iter;
if ~isfield(options,'err_thr'); options.err_thr = defoptions.err_thr; end
    err_thr = defoptions.err_thr;
if ~isfield(options,'eta'); options.eta = defoptions.eta; end
    eta = options.eta*max(Y(:))^2;
if ~isfield(options,'beta'); options.beta = defoptions.beta; end
    beta = options.beta;
if ~isfield(options,'nb'), options.nb = defoptions.nb; end
    nb = options.nb;

C = rand(K,T);

repeat = 1;
iter = 1;
obj_ = 1e-10;

v = ver;
flag_optim = any(strcmp('Optimization Toolbox', {v.Name})); % check if optimization toolbox is present
if flag_optim
    if verLessThan('optim','6.3')
        min_options = optimset('Algorithm','interior-point','GradObj','On','Display','Off');
    else
        min_options = optimoptions('fmincon','Algorithm','interior-point','GradObj','On','Display','Off');
    end
end

% remove median
medY = median(Y, 2);
Y = bsxfun(@minus, Y, medY);

while (iter <= max_iter) && repeat
    A = max((Y*C')*pinv(C*C'+beta*ones(K,K)),0);
    C = max((A'*A + eta*eye(K))\(A'*Y),0);

    ff = find(sum(C,2)==0);
    if ~isempty(ff)
        A(:,ff) = [];
        C(ff,:) = [];
        K = K - length(ff);
    end

    iter = iter + 1;
    if mod(iter,10) == 0;
        A = threshold_components(A,options);
        nC = sum(C.^2,2);
        AA = A'*A;
        mine = @(e) min_e(e,beta,eta,nC,A,AA);

        e = (eta*nC./(beta*A'*sum(A,2))).^(1/4);
        if flag_optim
            e = fmincon(mine,max(e,1e-4),[],[],[],[],1e-4*ones(K,1),[],[],min_options);
        end
        C = diag(e)\C;
        A = A*diag(e);
        fprintf('%i out of maximum %i iterations done \n',iter,max_iter);
    end
    obj = norm(Y - A*C,'fro')^2 + eta*norm(C,'fro')^2 + beta*norm(sum(A,2))^2;
    repeat = abs(obj - obj_) > err_thr*obj_;
    obj_ = obj;
end

fprintf('Algorithm converged after %i iterations. \n',iter-1);

[Ain,Cin] = order_components(A,C);
[bin,fin] = nnmf(max(Y - Ain*Cin + repmat(medY,1,T),0),nb);

    function [f,grad] = min_e(e,beta,eta,nC,A,AA)
        f = eta*norm(sqrt(nC)./e)^2 + beta*norm(A*e)^2;
        grad = -2*eta*nC./(e.^3) + 2*beta*AA*e;
    end

    function [A_or,C_or] = order_components(A,C)
        nA = sqrt(sum(A.^2));
        nr = length(nA);
        A = A/spdiags(nA(:),0,nr,nr);
        mA = max(A);
        C = spdiags(nA(:),0,nr,nr)*C;
        nC2 = sqrt(sum(C.^2,2));
        mC = max(C,[],2);
        [~,srt] = sort(mC.*mA'./nC2,'descend');
        A_or = A(:,srt);
        C_or = C(srt,:);
    end

end

author: zhoupc


text:
function [flat, allDist, sumDist, geodesicCoords, geodesicLayer] = measureCorticalDistance(flat, coords, verboseFlag, gray)
% USAGE: measureCorticalDistance(flat, [coords], [verboseFlag], [gray])
%
% AUTHOR:  Dougherty
% DATE:    2002.01.15
% PURPOSE:
%   Compute the shortest cortical manifold distance between
%   points. If coords is passed in, then the points are drawn
%   from there. Otherwise, the points are obtained via ginput.
%   (Note that coords must be a valid nx3 list of flat view coords,
%   such as would be obtained by 'getCurROIcoords(flatView)'.)
%
% HISTORY
%
% 7/16/02 djh, replaced mrSESSION.vAnatomyPath with global vANATOMYPATH
% 9/02/04 mms, inserted a warning not to use rotated flats.
% 9/22/04 rfd, now returns the geodesic path (as coords) and the gray layer
% of each point along the path.

global vANATOMYPATH;
mmPerPix = readVolAnatHeader(vANATOMYPATH);

% Check some flat and volume stuff here when you
% get a chance to fix this code up

% Get a gray structure because we need the gray nodes.
if notDefined('gray')
    gray = getSelectedGray;
    if isempty(gray)
        gray=initHiddenGray;
    end
end

if(~exist('coords','var') | isempty(coords))
    % Select flat figure and get a single point from the user
    figure(flat.ui.figNum)
    disp('Click left to add points, right to quit');
    button = 1;
    count = 0;
    w = 0.5;
    coords = [];
    z = viewGet(flat, 'Current Slice');
    while(button~=3)
        [x,y,button] = ginput(1);
        if(button==3)
            break;
        end
        count = count+1;
        coords = [coords,round([y;x;z])];
        h(count) = line([x-w,x-w,x+w,x+w,x-w],[y-w,y+w,y+w,y-w,y-w],'Color','w');
        if flat.rotateImageDegrees(coords(3,1))>0  %detects if the Flat your are currently selecting the points ist rotated
            errordlg('Flat must not be rotated','Warning!'); %Results would be wrong!
        end
    end
    % Delete the temporarily drawn squares
    for ii=1:length(h)
        delete(h(ii));
    end
    clear h;
end
if(~exist('verboseFlag','var') || isempty(verboseFlag))
    verboseFlag=1;
end

% the third coordinate is the 'slice', which, for flat views, means left or right hemisphere.
slice = coords(3,1);
if (slice==1)
    nodes = gray.allLeftNodes;
    edges = gray.allLeftEdges;
else
    nodes = gray.allRightNodes;
    edges = gray.allRightEdges;
end

% We loop for the number of line segements, which is the number of
% coords - 1.
geodesic = [];
for(ii=1:size(coords,2)-1)
    % get nearest flat coordinate (not all points on the flat correspond to flat coordinates)
    flatDistances = (flat.coords{slice}(1,:) - coords(1,ii)).^2 + ...
        (flat.coords{slice}(2,:) - coords(2,ii)).^2;
    % There is a one-to-many mapping of flatCoords to grayCoords, but we ignore that
    % here by using 'min', which will always reuturn one value, even if there are several
    % identical minima.
    % FIX THIS- we should always grab layer 1, or something more consistent
    % than relying on min's arbitrary sort.
    [val,startIndex] = min(flatDistances);
    % Do it again, for ii+1, to find the end point of this line segment.
    flatDistances = (flat.coords{slice}(1,:) - coords(1,ii+1)).^2 + ...
        (flat.coords{slice}(2,:) - coords(2,ii+1)).^2;
    [val,endIndex] = min(flatDistances);

    % Draw a line for each measured segment
    % we use the actual gray node coords rather than the ROI coords, so the use can see if
    % there is any non-trivial discrepancy.
    h(ii) = line([flat.coords{slice}(2,startIndex),flat.coords{slice}(2,endIndex)], ...
                 [flat.coords{slice}(1,startIndex),flat.coords{slice}(1,endIndex)], ...
                 'Color', 'r', 'LineWidth', 2);

    % Extract the gray node corresponding to the start
    startGrayNode = find(nodes(2,:) == flat.grayCoords{slice}(1,startIndex) & ...
                         nodes(1,:) == flat.grayCoords{slice}(2,startIndex) & ...
                         nodes(3,:) == flat.grayCoords{slice}(3,startIndex));
    endGrayNode = find(nodes(2,:) == flat.grayCoords{slice}(1,endIndex) & ...
                       nodes(1,:) == flat.grayCoords{slice}(2,endIndex) & ...
                       nodes(3,:) == flat.grayCoords{slice}(3,endIndex));

    % Catch errors. If we give mrManDist an empty startPoint array, it barfs.
    if(isempty(startGrayNode) | isempty(endGrayNode))
        myErrorDlg('No gray nodes were found for these coords!');
    end

    % Now, compute the manifold distance between these points.
    % mrManDist returns the distance to all other points from the given 'start' point.
    %allDist = mrManDist(nodes, edges, startGrayNode, mmPerPix, -1, 0);
    [allDist,nPts,lastPoint] = mrManDist(nodes, edges, startGrayNode, ...
                                      mmPerPix, -1,0);
    nextPoint = endGrayNode;
    while(nextPoint~=startGrayNode)
      geodesic(end+1) = lastPoint(nextPoint);
      nextPoint = lastPoint(nextPoint);
    end

    % We just want the distance from the start point to the end point, so we
    % pull that out by providing the index of the end point.
    dist(ii) = allDist(endGrayNode);
    if(verboseFlag>=0)
        disp(['Cortical distance of segment ',num2str(ii),': ',num2str(dist(ii)),' mm.']);
    end
end
sumDist=sum(dist);
if(verboseFlag>=0)
    disp(['Total cortical distance: ',num2str(sumDist),' mm.']);
end
if(verboseFlag>0)
    uiwait(msgbox(['Total cortical distance: ',num2str(sumDist),' mm.'], ...
        'Cortical Distance', 'modal'));
end

for ii=1:length(h)
    delete(h(ii))
end
geodesicCoords = nodes([2,1,3],geodesic);
geodesicLayer =  nodes(6,geodesic);
return


text:
function S = PMSAFinit(w0,mu,N,L,H,F,alpha,delta)

% PMSAFinit        Initialize Parameter Structure for the PMSAF Algorithm
%
%                  Psuedo-QMF CMFB is Used by default, H and F are used otherwise
% Arguments:
% w0               Coefficients of FIR filter at start
% mu               Step size
% N                Number of subbands
% L                length of analysis filter
% M                length of adaptive weight vector
% H                Analysis filter bank (optional), each column represents a filter
% F                Synthesis filter bank (optional), each column represents a filter
% alpha            Adjust scaling of tap weights
% delta            Small constant

if nargin > 4
    if (size(H,2)~=N)||(size(F,2)~=N)
        error('Columns of H (%d) or F (%d) not match with N = %d',size(H,2),size(F,2),N);
    end
else
% Defualt filter bank: Pseudo-QMF CMFB

    [hopt,passedge] = opt_filter(L-1,N); % Generate prototype lowpass filter
    [H,F] = make_bank(hopt,N);           % Generate filter banks using cosine modulation
    H = sqrt(N)*H';                      % Analysis section
    F = sqrt(N)*F';                      % Synthesis section
end

if nargin < 7 || isempty(alpha),         % Use default
    alpha = 0;                           % Try -0.5
end

if nargin < 8 || isempty(delta),
    delta = 1e-4;                        % Use default
end

% Assign structure fields

S.coeffs        = w0(:);                 % Convert to column vector of length M
S.step          = mu;                    % Step size
S.analysis      = H;                     % Analysis filter
S.synthesis     = F;                     % Synthesis filter
S.iter          = 0;                     % Iteration count
S.alpha         = alpha;                 % Tap weight scaling
S.delta         = ones(1,N)*delta;       % Small positive constant
S.AdaptStart    = L + length(w0);        % Running effect of analysis and adaptive filter,
                                         %   minimum L+M

author: CharlesThaCat

repo: acoustic-interference-cancellation

path: github-repos/MATLAB/CharlesThaCat-acoustic-interference-cancellation/acoustic-interference-cancellation-edb394499ea6f9c96445a3e9613bd64a854c289e/Subband processing/Book/After reading Chapter 5 -- Critically Sampled and Oversampled Subband Structures/PMSAFinit.m


text:
%gohome   positions the command prompt at top level directory
%   IFISS scriptfile: DJS; 24 May 2006.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage (see readme.m)
%
cd ( '/panfs/panasas1/users/jburkardt/public_html/m_src/toms866' );


text:
function test (varargin)
% TODO: in stack information, replace pointers to temporary files with
%       pointes to test section in the mfile being tested.
% TODO: add test case for get_stack_trace
% TODO: add test case for functions other than m-files that have
%       an accompanying homonymous m-file containing the test blocks.
% TODO: support subfunctions. (Nested functions will never be supported
%!      because they cannot stand alone.)
% TODO: is seems that a test block, at the end of a file, without a blank
%       line after it, is not supported.

    %%%%%%%%%%
    switch nargin
    case 0
        error('test:notEnoughInputs', ...
        'Not enough input arguments.');
    case 1
        mfunc = varargin{1};
        param = '';

        mfile = which(mfunc);
        if isempty(mfile)
            error('test:UndefinedFunction', ...
            'Function %s not found in path.', mfunc);
        end
        [mdir, mfunc, ext] = fileparts(mfile);
        if ~strcmp(ext, '.m')
            mfile = fullfile(mdir, [mfunc '.m']);
            if ~exist(mfile, 'file')
                error('test:FileNotFound',...
                'File not found: %s', mfile);
            end
        end

        %[tdir, tfunc] = fileparts(tempname);
        %tfile = fullfile(tdir, [tfunc '.m']);
        % under Vista, it's significantly faster to addpath
        % an empty directory rather than the standard temporary dir.
        tdir = tempname;
        tfunc = 'testit';
        tfile = fullfile(tdir, [tfunc '.m']);

        pdir = tdir;
        create_pdir = true;

        test_type = 'func';
    case 2
        mfunc = varargin{1};
        param = varargin{2};

        try
            eval(sprintf('%s;', param));
        catch
            param = sprintf('%s(1)', param);
        end
        mfunc2 = sprintf('%s(%s)', mfunc, param);
        mfile = which(mfunc2);
        if isempty(mfile)
            error('test:UndefinedFunction',...
            'Function %s not found in path.', mfunc2);
        end
        [mdir, mfunc, ext] = fileparts(mfile);
        if ~strcmp(ext, '.m')
            mfile = fullfile(mdir, [mfunc '.m']);
            if ~exist(mfile, 'file')
                error('test:FileNotFound',...
                'File not found: %s', mfile);
            end
        end

        %[tdir, tfunc] = fileparts(tempname);
        %tfile = fullfile(tdir, [tfunc '.m']);
        % under Vista, it's significantly faster to addpath
        % an empty directory rather than the standard temporary dir.
        tdir = tempname;
        tfunc = 'testit';
        tfile = fullfile(tdir, [tfunc '.m']);

        pdir = tdir;
        create_pdir = true;

        test_type = 'method';
    case 3
        mfunc = varargin{1};
        param = char(varargin{2});  % so that printf works when param = []
        mdir  = varargin{3};

        if ~isempty(param) && strcmp(mdir, 'private')
            mdir = [fileparts(which(param)) '\private'];
        end
        if ~(   strend([filesep 'private' filesep], mdir) ...
             || strend([filesep 'private']        , mdir) )
            error('test:DirEndInPrivate',...
            ['Third parameter %s should end in '...
            '"%cprivate".'], mdir, filesep);
        end

        mfile = fullfile(mdir, [mfunc '.m']);
        if ~exist(mfile, 'file')
            error('test:FileNotFound',...
            'File not found: %s', mfile);
        end

        tdir = dirup(mdir);
        create_pdir = false;
        while true
            [temp, tfunc] = fileparts(tempname); %#ok<ASGLU>
            tfile = fullfile(tdir, [tfunc '.m']);
            if ~exist(tfile, 'file'),  break;  end
        end

        if isempty(param)
            pdir = tdir;
            test_type = 'priv_func';
        else
            if ~strend([filesep '@' param], tdir)
                error('test:DirEndInClass',...
                ['Third parameter %s should end in '...
                '"%c@%s%cprivate".'], mdir, filesep, param, filesep);
            end
            pdir = dirup(tdir);
            test_type = 'priv_method';
        end
    otherwise
        error('test:tooManyInputs',...
        'Too many input parameters.');
    end
    %mdir, tdir, pdir  % DEBUG

    %%%%%%%%%%
    %mfile  % DEBUG
    %type(mfile)  % DEBUG
    [mfid, msg] = fopen (mfile, 'rt');
    if (mfid == -1)
        error('test:FileOpening', msg);
    end

    % extract blocks (test, error, shared) in m-file
    block_start = [];
    is_in_block = false;
    i = 1;  % # of %! lines
    j = 1;  % # of test or error blocks
    mline = fgets(mfid);
    while ~isequal(mline, -1)  % while not end of file
        if (   strncmp(mline, '%!test', 6) ...
            || strncmp(mline, '%!error', 7) ...
            || strncmp(mline, '%!shared', 8) )
            % start block
            is_in_block = true;
            block_start(j) = i;
            if strncmp(mline, '%!test', 6)
                block_type{j} = 'test';
            elseif strncmp(mline, '%!error', 7)
                block_type{j} = 'error';
            elseif strncmp(mline, '%!shared', 8)
                block_type{j} = 'shared';
            end
        elseif is_in_block && strncmp(mline, '%!', 2)
            % read block contents
            tline{i} = mline(4:end);
            i = i + 1;
        elseif is_in_block
            % end block
            block_end(j) = i-1;
            j = j + 1;
            is_in_block = false;
        end
        mline = fgets(mfid);
    end
    if is_in_block
        % end latest block
        block_end(j) = i-1;
    end
    fclose(mfid);

    if isempty(block_start) || all(strcmp(block_type, 'shared'))
        error('test:noBlock',...
        'Found no test or error section in file %s.', mfile);
    end
    myassert (length(block_start), length(block_end));
    myassert (length(block_start), length(block_type));

    %%%%%%%%%%

    %% create test file
    old_path = path;
    if create_pdir,  mkdir(pdir);  end
    addpath(pdir);
    %addpath(pdir, '-frozen');

    %% write and run each block
    for j=1:length(block_start)
        if strcmp(block_type{j}, 'shared'),  continue;  end

        [tfid, msg] = fopen(tfile, 'wt');
            if (tfid == -1),  error('test:FileOpening', msg);  end

        if strcmp(test_type, 'priv_method')
            fprintf(tfid, 'function %s (obj)\n', tfunc);
        else
            fprintf(tfid, 'function %s\n', tfunc);
        end

        found_mfunc_call = false;

        %% Write optional, previous, shared block.
        k = find (strcmp({block_type{1:j}}, 'shared'), 1, 'last');
        if ~isempty(k)
            myassert (strcmp(block_type{k}, 'shared'));
            fprintf(tfid, '%% shared\n');
            for i=block_start(k):block_end(k)
                fprintf(tfid, '%s', tline{i});
                found_mfunc_call = found_mfunc_call ...
                                   || is_function_called (tline{i}, mfunc);
            end
            fprintf(tfid, '\n');
        end

        %% Is it an error block?
        if strcmp(block_type{j}, 'error')
            fprintf(tfid, 'try\n');
            % In an error block, error is the expected behaviour.
            % We do we following:
            %     try,
            %         % INSERT USER CODE HERE
            %     catch,
            %         return;
            %     end
            %     error('test:error:fail', '');
            % If the user code generates an error, it is catched and
            % the test runs OK. If the user code does _not_ generate
            % an error(contrary to what was expected), then the test fails.
        end

        for i=block_start(j):block_end(j)
            fprintf(tfid, '%s', tline{i});
            found_mfunc_call = found_mfunc_call ...
                               || is_function_called (tline{i}, mfunc);
        end

        if (~found_mfunc_call)
            warning('test:noFuncCall', ...
                'Found no call to function %s itself in its %s block #%d.',...
                mfunc, block_type{j}, j);
        end

        if strcmp(block_type{j}, 'error')
            fprintf(tfid, 'catch\n');
            fprintf(tfid, '    return;\n');
            fprintf(tfid, 'end\n');
            fprintf(tfid, 'clear s;\n');
            fprintf(tfid, 's.identifier=''test:error:fail'';\n');
            fprintf(tfid, 's.message=''No error found in error section.'';\n');
            fprintf(tfid, 'error(s);\n');
        end

        fprintf(tfid, 'end\n');
        fclose(tfid);
        %type(tfile); % DEBUG

        % reload function definition:
        clear(tfile);
        if strcmp(test_type, 'priv_method')
            % we have added a new class methods; reload the class:
            clear(which(param));
        else
            clear(which(tfunc));
        end

        %% run test file
        try
            if strcmp(test_type, 'priv_method')
                eval(sprintf('%s(%s);', tfunc, param));
            else
                %which(tfunc),  type(which(tfunc))  % DEBUG
                feval(tfunc);
                %run([tdir filesep tfunc]);
            end
        catch %#ok<CTCH>
            s = lasterror; %#ok<LERR>
            idx = strcmp({s.stack.name}, mfunc);
            if strcmp(test_type, 'priv_func') && any(idx)
                s.stack(idx).file = mfile;  % it needs full path
            end
            type(which(tfunc));
            s.message = [s.message, get_stack_trace(s)];
            %if ~strcmp(s.identifier, 'MATLAB:UndefinedFunction')
                delete(tfile);
            %end
            path(old_path);
            if create_pdir,  rmdir(pdir);  end
            rethrow(s);
        end
    end

    %% remove test file
    if exist(tfile, 'file'),  delete(tfile);  end
    path(old_path);
    if create_pdir,  rmdir(pdir);  end
end

%!#test
%! % test()
%! myassert (false)

%!test
%! % test
%! myassert(true);
%! lasterr('', 'test:testing');

%!test
%! % test
%! s = lasterror;
%! myassert(s.identifier, 'test:testing');
%! lasterr('', '');

%!error
%! % test
%! myassert(false);

%!shared
%! var = pi;

%!test
%! % test
%! myassert (var == pi);
%! var = 2*pi;

%!test
%! % test
%! myassert (var == pi);

%!shared
%! var = 3*pi;

%!test
%! % test
%! myassert (var == 3*pi);

%!shared  % to clear previous %!shared

%!error
%! test;

%! % test
%! s = lasterror;
%! myassert(s.identifier, 'test:notEnoughInputs');

%!error
%! test(tempname);

%! % test
%! s = lasterror;
%! myassert(s.identifier, 'test:UndefinedFunction');

%!error
%! test(1, 2, 3, 4);

%! % test
%! s = lasterror;
%! myassert(s.identifier, 'test:tooManyInputs');

%!test
%! % It's not an error to have a single, empty test section.
%! func = tempname;
%! file = [func '.m'];
%!
%! fid = fopen(file, 'wt');
%! fprintf(fid, '%%!test\n');
%! fprintf(fid, '%%! %% %s\n', func);
%! fclose(fid);
%!
%! old_path = path;
%! addpath(fileparts(file));
%!
%! lasterr ('', '');
%! try,  test(func);  end
%! s = lasterror;
%! %type (file);  % DEBUG
%! delete (file);
%! path(old_path);
%! myassert (isempty(s.identifier));

%!test
%! % Issue an error if there is shared section but not test or error sections.
%! func = tempname;
%! file = [func '.m'];
%!
%! fid = fopen(file, 'wt');
%! fprintf(fid, '%%!shared\n');
%! fprintf(fid, '%%! \n');
%! fclose(fid);
%!
%! old_path = path;
%! addpath(fileparts(file));
%!
%! lasterr ('', '');
%! try,  test(func);  end
%! s = lasterror;
%! %type (file);  % DEBUG
%! delete (file);
%! path(old_path);
%! myassert (s.identifier, 'test:noBlock');

%!shared
%! % Here we generate the following directory-tree:
%! % dir/
%! %     @aclass/
%! %         aclass.m  % class constructor
%! %         afunc.m
%! % Then test() creates:
%! % dir/
%! %     atest.m
%!
%! class_code = {...
%!     'function answer = aclass (value)'
%!     '    if (nargin < 1),  value = [];  end'
%!     '    answer.value = value;'
%!     '    answer = class(answer, ''aclass'');'
%!     'end'
%! };
%! func_code = {...
%!     'function answer = afunc (obj)'
%!     '    answer = obj.value;'
%!     'end'
%!     ''
%!     '%!shared'
%!     '%! value = pi;'
%!     '%! obj = aclass(value);'
%!     ''
%!     '%!error'
%!     '%! % afunc'
%!     '%! myassert(false);'
%!     ''
%!     '%!test'
%!     '%! myassert(afunc(obj), value);'
%! };
%! %class_code, func_code  % DEBUG
%!
%! dir = tempname;
%! class_dir = fullfile(dir, '@aclass');
%! temp = mkdir(class_dir);
%!     myassert(temp ~= 0);
%!
%! class_file = fullfile(class_dir, 'aclass.m');
%! [class_fid, msg] = fopen (class_file, 'wt');
%!     if (class_fid == -1),  error(msg);  end
%! fprintf(class_fid, '%s\n', class_code{:});
%! temp = fclose(class_fid);
%!     myassert (temp ~= -1);
%!
%! func_dir = class_dir;
%! func_file = fullfile(func_dir, 'afunc.m');
%! func_fid = fopen (func_file, 'wt');
%!     myassert (func_fid ~= -1);
%! fprintf(func_fid, '%s\n', func_code{:});
%! temp = fclose(func_fid);
%!     myassert (temp ~= -1);
%!
%! %type(class_file);  type(func_file);  % DEBUG
%! %pause  % DEBUG
%!
%! old_path = path;
%! addpath (dir);  % NOT addpath (class_dir);

%!test
%! % Test a class method:
%! % test('afunc', 'aclass');
%!
%! try,
%!     test('afunc', 'aclass');
%! catch,
%!     s = lasterror;
%!     if ~strcmp(s.identifier, 'MATLAB:UndefinedFunction')
%!         path (old_path);
%!         [status, msg, msgid] = rmdir(dir, 's');
%!             if (status == 0),  error(msgid, msg);  end
%!     end
%!     rethrow (s);
%! end
%!
%! path (old_path);
%! [status, msg, msgid] = rmdir(dir, 's');
%!     if (status == 0),  error(msgid, msg);  end

%!test
%! % Test a class private method:
%! % test('afunc', 'aclass', 'adir');
%! %
%! % Here we generate the following directory-tree:
%! % dir/
%! %     @aclass/
%! %         aclass.m  % class constructor
%! %         private/
%! %             afunc.m
%! % Then test() creates:
%! % dir/
%! %     @aclass/
%! %         atest.m
%!
%! func_file_old = func_file;
%! func_dir = fullfile(class_dir, 'private');
%! temp = mkdir(func_dir);
%!     myassert(temp ~= 0);
%! func_file = fullfile(func_dir, 'afunc.m');
%! [status, msg, msgid] = movefile (func_file_old, func_file);
%!     if (status ~= 1),  error(msgid, msg);  end
%!
%! try,
%!     %keyboard  % DEBUG
%!     test('afunc', 'aclass', func_dir);
%! catch,
%!     s = lasterror;
%!     if ~strcmp(s.identifier, 'MATLAB:UndefinedFunction')
%!         path (old_path);
%!         [status, msg, msgid] = rmdir(dir, 's');
%!             if (status == 0),  error(msgid, msg);  end
%!     end
%!     rethrow (s);
%! end
%!
%! path (old_path);
%! [status, msg, msgid] = rmdir(dir, 's');
%!     if (status == 0),  error(msgid, msg);  end

%!test
%! % Test a class private method, letting test() guess adir:
%! % test('afunc', 'aclass', 'private');
%! %
%! % Here we generate the following directory-tree:
%! % dir/
%! %     @aclass/
%! %         aclass.m  % class constructor
%! %         private/
%! %             afunc.m
%! % Then test() creates:
%! % dir/
%! %     @aclass/
%! %         atest.m
%!
%! func_file_old = func_file;
%! func_dir = fullfile(class_dir, 'private');
%! temp = mkdir(func_dir);
%!     myassert(temp ~= 0);
%! func_file = fullfile(func_dir, 'afunc.m');
%! [status, msg, msgid] = movefile (func_file_old, func_file);
%!     if (status ~= 1),  error(msgid, msg);  end
%!
%! try,
%!     %keyboard  % DEBUG
%!     test('afunc', 'aclass', 'private');
%! catch,
%!     s = lasterror;
%!     if ~strcmp(s.identifier, 'MATLAB:UndefinedFunction')
%!         path (old_path);
%!         [status, msg, msgid] = rmdir(dir, 's');
%!             if (status == 0),  error(msgid, msg);  end
%!     end
%!     rethrow (s);
%! end
%!
%! path (old_path);
%! [status, msg, msgid] = rmdir(dir, 's');
%!     if (status == 0),  error(msgid, msg);  end

%!shared  % to clear previous %!shared

%!test
%! % Test a private function:
%! % test('afunc', [], 'adir');
%! %
%! % Here we generate the following directory-tree:
%! % dir/
%! %     private/
%! %         afunc.m
%! % Then test() creates:
%! % dir/
%! %     atest.m
%!
%! func_code = {...
%!     'function out = afunc (in)'
%!     '    if isequal(in, 0),  error(''afunc:any'', '''');  end'
%!     '    out = 1/in;'
%!     'end'
%!     ''
%!     '%!error'
%!     '%! afunc(0)';
%!     ''
%!     '%!test'
%!     '%! myassert(afunc(pi), 1/pi);'
%! };
%! %func_code  % DEBUG
%!
%! dir = tempname;
%! temp = mkdir(dir);
%!     myassert(temp ~= 0);
%! temp = mkdir(dir, 'private');
%!     myassert(temp ~= 0);
%!
%! func_dir = fullfile(dir, 'private');
%! func_file = fullfile(func_dir, 'afunc.m');
%! [func_fid, msg] = fopen (func_file, 'wt');
%!     if (func_fid == -1),  error(msg);  end
%! fprintf(func_fid, '%s\n', func_code{:});
%! temp = fclose(func_fid);
%!     myassert (temp ~= -1);
%! %type(func_file);  % DEBUG
%!
%! old_path = path;
%! addpath (dir);  % NOT addpath (func_dir);
%! test('afunc', [], func_dir);
%! path (old_path);
%!
%! [status, msg, msgid] = rmdir(dir, 's');
%!     if (status == 0),  error(msgid, msg);  end

%!#test
%! % Check if test files are being deleted,
%! % even in case of error.
%!
%! code = {...
%!     'function answer = func'
%!     '    persistent temp'
%!     '    if isempty(temp)'
%!     '        temp = which(get_caller_name);'
%!     '    end'
%!     '    answer = temp;'
%!     'end'
%! };
%! test1 = {...
%!     '%!test'
%!     '%! func;'
%!     '%! myassert(true);'
%! };
%! test2 = {...
%!     '%!test'
%!     '%! func;'
%!     '%! myassert(false);  % THROW ERROR'
%! };
%!
%! dir = tempname;
%! mkdir(dir);
%! file = fullfile(dir, 'func.m');
%!
%! old_path = path;
%! addpath(dir);
%!
%! for i=1:2
%!     switch i
%!     case 1
%!         testi = test1;
%!     case 2
%!         testi = test2;
%!     end
%!
%!     fid = fopen(file, 'wt');
%!     fprintf(fid, '%s\n', code{:});
%!     fprintf(fid, '%s\n', testi{:});
%!     fclose(fid);
%!
%!     clear('func');
%!     try,  test('func');  end
%!
%!     tfile = func;
%!     try,
%!         switch i
%!         case 1
%!             myassert (~exist(tfile, 'file'));
%!         case 2
%!             myassert (~exist(tfile, 'file'));
%!         end
%!     catch,
%!         s = lasterror;
%!         path(old_path);
%!         rmdir(dir, 's');
%!         rethrow (s);
%!     end
%! end
%! path(old_path);
%! rmdir(dir, 's');

%!test
%! % test ()
%! lastwarn('', '');  % clear last warning
%! warning('off', 'test:noFuncCall');

%!test
%! % generate warning

%!test
%! % test ()
%! [msg, msgid] = lastwarn;
%! myassert (msgid, 'test:noFuncCall');  % check warning
%! warning('on', 'test:noFuncCall');

%!test
%! % test ()
%! lastwarn('', '');  % clear last warning

%!shared
%! % test ()  % avoid test()'s warning

%!test
%! [msg, msgid] = lastwarn;
%! myassert (msgid, '');  % check warning is still clear

author: Sable

repo: mcbench-benchmarks


text:
function [out1, out2] = gpS00(X, P);
%%function [out1, out2] = gpS00(X, {input, target, test});

% gpS00: Gaussian process regression with "squared negative exponential"
% covariance function and independent Gaussian noise model. Two modes are
% possible: training and prediction: if no test data are given, the function
% returns minus the log likelihood and its partial derivatives with respect to
% the hyperparameters; this mode is used to fit the hyperparameters. If test
% data are given, then (marginal) Gaussian predictions are computed, whose mean
% and (noise free) variance are returned.
%
% usage: [fX dfX] = gpS00(X, input, target)
%    or: [mu S2]  = gpS00(X, input, target, test)
%
% where:
%
%   X      is a (column) vector (of size D+2) of hyperparameters
%   input  is a n by D matrix of training inputs
%   target is a (column) vector (of size n) of targets
%   test   is a nn by D matrix of test inputs
%   fX     is the returned value of minus log likelihood
%   dfX    is a (column) vector (of size D+2) of partial derivatives
%            of minus the log likelihood wrt each of the hyperparameters
%   mu     is a (column) vector (of size nn) of prediced means
%   S2     is a (column) vector (of size nn) of predicted variances
%
% where D is the dimension of the input. The form of the covariance function is
%
% C(x^p,x^q) = v^2 * exp[-(x^p - x^q)'*inv(P)*(x^p - x^q)/2]
%            + u^2 * delta_{p,q}
%
% where the first term is the squared negative exponential and the second term
% with the kronecker delta is the noise contribution. The P matrix is diagonal
% with "Automatic Relevance Determination" (ARD) or "input length scale"
% parameters w_1^2,...,w_D^2; The hyperparameter v is the "signal std dev" and
% u is the "noise std dev". All hyperparameters are collected in the vector X
% as follows:
%
% X = [ log(w_1)
%       log(w_2)
%        .
%       log(w_D)
%       log(v)
%       log(u) ]
%
% Note: the reason why the log of the parameters are used in X is that this
% often leads to a better conditioned (and unconstrained) optimization problem
% than using the raw hyperparameters themselves.
%
% This function can conveniently be used with the "minimize" function to train
% a Gaussian process:
%
% [X, fX, i] = minimize(X, 'gpS00', length, input, target)
%
% See also: minimize, hybrid
%
% (C) Copyright 1999 - 2003, Carl Edward Rasmussen (2003-08-08).

input=P{1};
target=P{2};
if length(P) == 3
	test=P{3};
	testmode=1;
else
	testmode=0;
end

[n, D] = size(input);         % number of examples and dimension of input space
input = input ./ repmat(exp(X(1:D))',n,1);

% first, we write out the covariance matrix Q

warning off;

Q = zeros(n,n);
for d = 1:D
  Q = Q + (repmat(input(:,d),1,n)-repmat(input(:,d)',n,1)).^2;
end
Q = exp(2*X(D+1))*exp(-0.5*Q);

if ~testmode   % if no test cases, we compute the negative log likelihood ...

  W = (Q+exp(2*X(D+2))*eye(n))\eye(n);               % W is inv (Q plus noise term)
  invQt = W*target;                               % don't compute determinant..
  [C,flag]=chol(Q+exp(2*X(D+2))*eye(n));
  logdetQ = 2*sum(log(diag(C)));        % ..directly
  out1 = 0.5*logdetQ + 0.5*target'*invQt + 0.5*n*log(2*pi);

  % ... and its partial derivatives

  out2 = zeros(D+2,1);                  % set the size of the derivative vector
  W = W-invQt*invQt';
  Q = W.*Q;
  for d = 1:D
    out2(d) = ...
            sum(sum(Q.*(repmat(input(:,d),1,n)-repmat(input(:,d)',n,1)).^2))/2;
  end
  out2(D+1) = sum(sum(Q));
  out2(D+2) = trace(W)*exp(2*X(D+2));

else                    % ... otherwise compute (marginal) test predictions ...

  [nn, D] = size(test);     % number of test cases and dimension of input space
  test = test ./ repmat(exp(X(1:D))',nn,1);

  a = zeros(n, nn);    % compute the covariance between training and test cases
  for d = 1:D
    a = a + (repmat(input(:,d),1,nn)-repmat(test(:,d)',n,1)).^2;
  end
  a = exp(2*X(D+1))*exp(-0.5*a);

  % ... write out the desired terms

  if nargout == 1
    out1 = a'*((Q+exp(2*X(D+2))*eye(n))\target);              % predicted means
  else
    invQ = inv(Q+exp(2*X(D+2))*eye(n));
    out1 = a'*(invQ*target);                                  % predicted means
    out2 = exp(2*X(D+1)) - sum(a.*(invQ*a),1)'; % predicted noise-free variance
  end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


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

%%% H-Matrix (recursion)
if (Mh.typ == 0)
    for i = 1:4
        Mh.chd{i} = hmxCopy(Mh.chd{i},tol);
        Mh.tol    = max(tol,Mh.tol);
    end
    Mh = hmxFusion(Mh);

%%% Compressed leaf
elseif (Mh.typ == 1)
    if (tol > Mh.tol)
        [A,B]  = hmxQRSVD(Mh.dat{1},Mh.dat{2},tol);
        Mh.dat = {A,B};
        Mh.tol = tol;
    end

%%% Full leaf
elseif (Mh.typ == 2)

%%% Unknown type
else
    error('hmxCopy.m : unavailable case')
end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function gmsh_io_test02 ( )

%*****************************************************************************80
%
%% GMSH_IO_TEST02 reads the example data from a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2014
%
%  Author:
%
%   John Burkardt
%
  gmsh_filename = 'example_2d.msh';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'GMSH_IO_TEST02:\n' );
  fprintf ( 1, '  Read data from a file.\n' );
%
%  Get the data size.
%
  [ node_num, m, element_num, element_order ] = ...
    gmsh_size_read ( gmsh_filename );
%
%  Print the sizes.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Node data read from file "%s"\n', gmsh_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes = %d\n', node_num );
  fprintf ( 1, '  Spatial dimension = %d\n', m );
  fprintf ( 1, '  Number of elements = %d\n', element_num );
  fprintf ( 1, '  Element order = %d\n', element_order );
%
%  Get the data.
%
  [ node_x, element_node ] = gmsh_data_read ( gmsh_filename, m, node_num, ...
    element_order, element_num );
%
%  Print some of the data.
%
  r8mat_transpose_print_some ( m, node_num, node_x, ...
    1, 1, m, 10, '  Coordinates for first 10 nodes:' );

  i4mat_transpose_print_some ( element_order, element_num, element_node, ...
    1, 1, element_order, 10, '  Connectivity for first 10 elements:' );

  return
end


text:
%YUV Class to read YUV4MPEG file
%
% A concrete subclass of ImageSource that returns images from a YUV4MPEG format
% uncompressed video file.
%
% Methods::
% grab    Aquire and return the next image
% size    Size of image
% close   Close the image source
% char    Convert the object parameters to human readable string
%
% Properties::
% curFrame        The index of the frame just read
%
% See also ImageSource, Video.
%
%
% SEE ALSO: Video

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

classdef YUV < ImageSource

    properties

        rate            % frame rate at which movie was capture

        skippedFrames

        curFrame
        skip

        fp
        hdr

    end

    methods

        function yuv = YUV(filename, varargin)
        %YUV.YUV YUV4MPEG sequence constructor
        %
        % Y = YUV(FILE, OPTIONS) is a YUV4MPEG object that returns frames
        % from the yuv4mpeg format file FILE.  This file contains uncompressed
        % color images in 4:2:0 format, with a full resolution luminance plane
        % followed by U and V planes at half resolution both directions.
        %
        % Options::
        % 'uint8'     Return image with uint8 pixels (default)
        % 'float'     Return image with float pixels
        % 'double'    Return image with double precision pixels
        % 'grey'      Return greyscale image
        % 'gamma',G   Apply gamma correction with gamma=G
        % 'scale',S   Subsample the image by S in both directions
        % 'skip',S    Read every S'th frame from the movie

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

            yuvm.curFrame = 1;
            yuv.skip = 1;

            yuv.fp = fopen(filename, 'r');

            hdr = fgets(yuv.fp);
            yuv.hdr = hdr;
            while length(hdr) > 1,
                [s, hdr] = strtok(hdr);
                switch s(1),
                case {'Y'}
                    if strcmp(s, 'YUV4MPEG2') == 0,
                        fclose(yuv.fp);
                        error('not a YUV4MPEG stream');
                    end
                case {'W'}
                    yuv.width = str2num(s(2:end));
                case {'H'}
                    yuv.height = str2num(s(2:end));
                otherwise
                    fprintf('found <%s>\n', s);
                end
            end
            yuv.curFrame = 0;
        end

        % destructor
        function delete(m)
            close(yuv.fp);
        end

        function close(m)
        %Movie.close Close the image source
        %
        % M.close() closes the connection to the movie.

            delete(m.movie);
        end

        function sz = size(yuv)
            sz = [yuv.width yuv.height];
        end

        function [o1,o2,o3] = grab(yuv, varargin)
        %Movie.grab Acquire next frame from movie
        %
        % IM = Y.grab(OPTIONS) is the next frame from the file.
        %
        % [Y,U,V] = Y.grab(OPTIONS) is the next frame from the file
        %
        % Options::
        % 'skip',S    Skip frames, and return current+S frame (default 1)
        % 'rgb'       Return as an RGB image, Y image is downsized by two (default).
        % 'rgb2'      Return as an RGB image, U and V images are upsized by two.
        % 'yuv'       Return Y, U and V images.
        %
        % Notes::
        % - If no output argument given the image is displayed using IDISP.
        % - For the 'yuv' option three output arguments must be given.

            opt.skip = yuv.skip;
            opt.mode = {'rgb', 'rgb2', 'yuv'};

            opt = tb_optparse(opt, varargin);

            % read next frame from the file
            [Y,U,V] = yuv.yuvread(opt.skip);
            if isempty(Y)
                out = [];
                return;
            end

            yuv.curFrame = yuv.curFrame + opt.skip;

            % apply options specified at construction time
            Y = yuv.convert(Y);
            U = yuv.convert(U);
            V = yuv.convert(V);

            % convert to RGB if required
            switch opt.mode
            case 'rgb'
                [R,G,B] = yuv2rgb(Y,U,V);
                o1 = cat(3, R, G, B);
            case 'rgb2'
                [R,G,B] = yuv2rgb2(Y,U,V);
                o1 = cat(3, R, G, B);
            case 'yuv'
                o1 = Y;
                o2 = U;
                o3 = V;
            end

        end

        function s = char(yuv)
        %Movie.char Convert to string
        %
        % M.char() is a string representing the state of the movie object in
        % human readable form.

            s = '';
            s = strvcat(s, sprintf('YUV file: %d x %d', yuv.width, yuv.height));
            s = strvcat(s, sprintf('cur frame %d (skip=%d)', yuv.curFrame, yuv.skip));
        end

        function paramSet(m, varargin)
            opt.skip = 1;

            disp(varargin)
            opt = tb_optparse(opt, varargin);
            m.skip = opt.skip;
        end

    end % methods

    methods(Access=protected)

        %YUVREAD	Read frame from a YUV4MPEG file
        %
        %	[y,u,v] = yuvread(yuv, skip)
        %	[y,u,v, h] = yuvread(yuv, skip)
        %
        %	Returns the Y, U and V components from the specified frame of
        %	YUV file.  Optionally returns the frame header h.

        function [Y,U,V, h] = yuvread(yuv, skip)

            while skip > 0
                % read and display the header
                hdr = fgets(yuv.fp);
                if isempty(hdr)
                    Y = [];
                    return;
                end
                fprintf('header: %s', hdr);

                % read the YUV data
                [Y,count] = fread(yuv.fp, yuv.width*yuv.height, 'uchar');
                if count ~= yuv.width*yuv.height
                    Y = [];
                    return;
                end
                [V,count] = fread(yuv.fp, yuv.width*yuv.height/4, 'uchar');
                if count ~= yuv.width*yuv.height/4
                    Y = [];
                    return;
                end
                [U,count] = fread(yuv.fp, yuv.width*yuv.height/4, 'uchar');
                if count ~= yuv.width*yuv.height/4
                    Y = [];
                    return;
                end

                skip = skip - 1;
            end

            Y = reshape(Y, yuv.width, yuv.height)';
            U = reshape(U, yuv.width/2, yuv.height/2)';
            V = reshape(V, yuv.width/2, yuv.height/2)';

            if nargout == 4
                h = hdr;
            end
        end

        %YUV2RGB	Convert YUV format to RGB
        %
        %	[r,g,b] = yuvread(y, u, v)
        %	rgb = yuvread(y, u, v)
        %
        %	Returns the equivalent RGB image from YUV components.  The Y image is
        %	halved in resolution.
        function [R,G,B] = yuv2rgb(y, u, v)

            % subsample, probably should smooth first...
            y = y(1:2:end, 1:2:end);

            Cr = u - 128;
            Cb = v - 128;

            % convert to RGB
            r = (y + 1.366*Cr - 0.002*Cb);
            g = (y - 0.700*Cr - 0.334*Cb);
            b = (y - 0.006*Cr + 1.732*Cb);

            % clip the values into range [0, 255]
            r = max(0, min(r, 255));
            g = max(0, min(g, 255));
            b = max(0, min(b, 255));

            if nargout == 1,
                R(:,:,1) = r;
                R(:,:,2) = g;
                R(:,:,3) = b;
            else
                R = r;
                G = g;
                B = b;
            end
        end

        %YUV2RGB	Convert YUV format to RGB
        %
        %	[r,g,b] = yuvread2(y, u, v)
        %	rgb = yuvread(y, u, v)
        %
        %	Returns the equivalent RGB image from YUV components.  The UV images are
        %	doubled in resolution so the resulting color image is original size.
        %
        function [R,G,B] = yuv2rgb2(y, u2, v2)

            % subsample, probably should smooth first...
            u = zeros(2*size(u2));
            v = zeros(2*size(v2));

            for i=1:4,
                u(1:2:end,1:2:end) = u2;
                v(1:2:end,1:2:end) = v2;

                u(2:2:end,1:2:end) = u2;
                v(2:2:end,1:2:end) = v2;

                u(1:2:end,2:2:end) = u2;
                v(1:2:end,2:2:end) = v2;

                u(2:2:end,2:2:end) = u2;
                v(2:2:end,2:2:end) = v2;
            end

            Cr = u - 128;
            Cb = v - 128;

            % convert to RGB
            r = (y + 1.366*Cr - 0.002*Cb);
            g = (y - 0.700*Cr - 0.334*Cb);
            b = (y - 0.006*Cr + 1.732*Cb);

            % clip the values into range [0, 255]
            r = max(0, min(r, 255));
            g = max(0, min(g, 255));
            b = max(0, min(b, 255));

            if nargout == 1,
                R(:,:,1) = r;
                R(:,:,2) = g;
                R(:,:,3) = b;
            else
                R = r;
                G = g;
                B = b;
            end
        end
    end % protected methods
end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function test_failed = test_libltfat_fftifftshift(varargin)
test_failed = 0;

fprintf(' ===============  %s ================ \n',upper(mfilename));

definput.flags.complexity={'double','single'};
[flags]=ltfatarghelper({},definput,varargin);
dataPtr = [flags.complexity, 'Ptr'];

Larr = [1,9,11,110,111];

for L = Larr

        z = cast((1:L)',flags.complexity);
        zi = complex2interleaved(fft(z));
        zout = randn(size(zi),flags.complexity);

        ziPtr = libpointer(dataPtr,zi);
        zoutPtr = libpointer(dataPtr,zout);

        trueres = ifftshift(z);

        funname = makelibraryname('fftifftshift',flags.complexity,1);
        status = calllib('libltfat',funname,ziPtr,L,zoutPtr);

        res = norm(trueres - real(ifft(interleaved2complex(zoutPtr.Value))));

        [test_failed,fail]=ltfatdiditfail(res+status,test_failed);
        fprintf(['FFTIFFTSHIFT OP L:%3i, %s %s %s\n'],L,flags.complexity,ltfatstatusstring(status),fail);

        status = calllib('libltfat', funname, ziPtr,L,ziPtr);

        res = norm(trueres - real(ifft(interleaved2complex(ziPtr.Value))));

        [test_failed,fail]=ltfatdiditfail(res+status,test_failed);
        fprintf(['FFTIFFTSHIFT IP L:%3i, %s %s %s\n'],L,flags.complexity,ltfatstatusstring(status),fail);

        z = cast((1:L)',flags.complexity);
        zi = complex2interleaved(fftreal(z));
        zout = randn(size(zi),flags.complexity);

        ziPtr = libpointer(dataPtr,zi);
        zoutPtr = libpointer(dataPtr,zout);

        trueres = ifftshift(z);

        funname = makelibraryname('fftrealifftshift',flags.complexity,1);
        status = calllib('libltfat',funname,ziPtr,L,zoutPtr);

        res = norm(trueres - ifftreal(interleaved2complex(zoutPtr.Value),L));

        [test_failed,fail]=ltfatdiditfail(res+status,test_failed);
        fprintf(['FFTREALIFFTSHIFT OP L:%3i, %s %s %s\n'],L,flags.complexity,ltfatstatusstring(status),fail);

        status = calllib('libltfat', funname, ziPtr,L,ziPtr);

        res = norm(trueres - ifftreal(interleaved2complex(ziPtr.Value),L));

        [test_failed,fail]=ltfatdiditfail(res+status,test_failed);
        fprintf(['FFTREALIFFTSHIFT IP L:%3i, %s %s %s\n'],L,flags.complexity,ltfatstatusstring(status),fail);

end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [sFile, ChannelMat, ImportOptions] = in_fopen_edf(DataFile, ImportOptions)
% IN_FOPEN_EDF: Open a BDF/EDF file (continuous recordings)
%
% USAGE:  [sFile, ChannelMat, ImportOptions] = in_fopen_edf(DataFile, ImportOptions)

% @=============================================================================
% 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, 2012-2021

% Parse inputs
if (nargin < 2) || isempty(ImportOptions)
    ImportOptions = db_template('ImportOptions');
end

%% ===== READ HEADER =====
% Open file
fid = fopen(DataFile, 'r', 'ieee-le');
if (fid == -1)
    error('Could not open file');
end
% Read all fields
hdr.version    = fread(fid, [1  8], 'uint8=>char');  % Version of this data format ('0       ' for EDF, [255 'BIOSEMI'] for BDF)
hdr.patient_id = fread(fid, [1 80], '*char');  % Local patient identification
hdr.rec_id     = fread(fid, [1 80], '*char');  % Local recording identification
hdr.startdate  = fread(fid, [1  8], '*char');  % Startdate of recording (dd.mm.yy)
hdr.starttime  = fread(fid, [1  8], '*char');  % Starttime of recording (hh.mm.ss)
hdr.hdrlen     = str2double(fread(fid, [1 8], '*char'));  % Number of bytes in header record
hdr.unknown1   = fread(fid, [1 44], '*char');             % Reserved ('24BIT' for BDF)
hdr.nrec       = str2double(fread(fid, [1 8], '*char'));  % Number of data records (-1 if unknown)
hdr.reclen     = str2double(fread(fid, [1 8], '*char'));  % Duration of a data record, in seconds
hdr.nsignal    = str2double(fread(fid, [1 4], '*char'));  % Number of signals in data record
% Check file integrity
if isnan(hdr.nsignal) || isempty(hdr.nsignal) || (hdr.nsignal ~= round(hdr.nsignal)) || (hdr.nsignal < 0)
    error('File header is corrupted.');
end
% Read values for each nsignal
for i = 1:hdr.nsignal
    hdr.signal(i).label = strtrim(fread(fid, [1 16], '*char'));
end
for i = 1:hdr.nsignal
    hdr.signal(i).type = strtrim(fread(fid, [1 80], '*char'));
end
for i = 1:hdr.nsignal
    hdr.signal(i).unit = strtrim(fread(fid, [1 8], '*char'));
end
for i = 1:hdr.nsignal
    hdr.signal(i).physical_min = str2double(fread(fid, [1 8], '*char'));
end
for i = 1:hdr.nsignal
    hdr.signal(i).physical_max = str2double(fread(fid, [1 8], '*char'));
end
for i = 1:hdr.nsignal
    hdr.signal(i).digital_min = str2double(fread(fid, [1 8], '*char'));
end
for i = 1:hdr.nsignal
    hdr.signal(i).digital_max = str2double(fread(fid, [1 8], '*char'));
end
for i = 1:hdr.nsignal
    hdr.signal(i).filters = strtrim(fread(fid, [1 80], '*char'));
end
for i = 1:hdr.nsignal
    hdr.signal(i).nsamples = str2num(fread(fid, [1 8], '*char'));
end
for i = 1:hdr.nsignal
    hdr.signal(i).unknown2 = fread(fid, [1 32], '*char');
end
% Unknown record size, determine correct nrec
if (hdr.nrec == -1)
    datapos = ftell(fid);
    fseek(fid, 0, 'eof');
    endpos = ftell(fid);
    hdr.nrec = floor((endpos - datapos) / (sum([hdr.signal.nsamples]) * 2));
end
% Close file
fclose(fid);

%% ===== RECONSTRUCT INFO =====
% Individual signal gain
for i = 1:hdr.nsignal
    % Interpet units
    switch (hdr.signal(i).unit)
        case 'mV',                        unit_gain = 1e3;
        case {'uV', char([166 204 86])},  unit_gain = 1e6;
        otherwise,                        unit_gain = 1;
    end
    % Check min/max values
    if isempty(hdr.signal(i).digital_min) || isnan(hdr.signal(i).digital_min)
        disp(['EDF> Warning: The digitial minimum is not set for channel "' hdr.signal(i).label '".']);
        hdr.signal(i).digital_min = -2^15;
    end
    if isempty(hdr.signal(i).digital_max) || isnan(hdr.signal(i).digital_max)
        disp(['EDF> Warning: The digitial maximum is not set for channel "' hdr.signal(i).label '".']);
        hdr.signal(i).digital_max = -2^15;
    end
    if isempty(hdr.signal(i).physical_min) || isnan(hdr.signal(i).physical_min)
        disp(['EDF> Warning: The physical minimum is not set for channel "' hdr.signal(i).label '".']);
        hdr.signal(i).physical_min = hdr.signal(i).digital_min;
    end
    if isempty(hdr.signal(i).physical_max) || isnan(hdr.signal(i).physical_max)
        disp(['EDF> Warning: The physical maximum is not set for channel "' hdr.signal(i).label '".']);
        hdr.signal(i).physical_max = hdr.signal(i).digital_max;
    end
    if (hdr.signal(i).physical_min >= hdr.signal(i).physical_max)
        disp(['EDF> Warning: Physical maximum larger than minimum for channel "' hdr.signal(i).label '".']);
        hdr.signal(i).physical_min = hdr.signal(i).digital_min;
        hdr.signal(i).physical_max = hdr.signal(i).digital_max;
    end
    % Calculate and save channel gain
    hdr.signal(i).gain   = unit_gain ./ (hdr.signal(i).physical_max - hdr.signal(i).physical_min) .* (hdr.signal(i).digital_max - hdr.signal(i).digital_min);
    hdr.signal(i).offset = hdr.signal(i).physical_min ./ unit_gain - hdr.signal(i).digital_min ./ hdr.signal(i).gain;
    % Error: The number of samples is not specified
    if isempty(hdr.signal(i).nsamples)
        % If it is not the first electrode: try to use the previous one
        if (i > 1)
            disp(['EDF> Warning: The number of samples is not specified for channel "' hdr.signal(i).label '".']);
            hdr.signal(i).nsamples = hdr.signal(i-1).nsamples;
        else
            error(['The number of samples is not specified for channel "' hdr.signal(i).label '".']);
        end
    end
    hdr.signal(i).sfreq = hdr.signal(i).nsamples ./ hdr.reclen;
end
% Find annotations channel
iAnnotChans = find(strcmpi({hdr.signal.label}, 'EDF Annotations'));  % Mutliple "EDF Annotation" channels allowed in EDF+
iStatusChan = find(strcmpi({hdr.signal.label}, 'Status'), 1);        % Only one "Status" channel allowed in BDF
iOtherChan = setdiff(1:hdr.nsignal, [iAnnotChans iStatusChan]);
% % Remove channels with lower sampling rates
% iIgnoreChan = find([hdr.signal(iOtherChan).sfreq] < max([hdr.signal(iOtherChan).sfreq]));    % Ignore all the channels with lower sampling rate
% if ~isempty(iIgnoreChan)
%     iOtherChan = setdiff(iOtherChan, iIgnoreChan);
% end
% Get all the other channels
if isempty(iOtherChan)
    error('This file does not contain any data channel.');
end
% Read events preferencially from the EDF Annotations track
if ~isempty(iAnnotChans)
    iEvtChans = iAnnotChans;
elseif ~isempty(iStatusChan)
    iEvtChans = iStatusChan;
else
    iEvtChans = [];
end
% % Detect channels with inconsistent sampling frenquency
% iErrChan = find([hdr.signal(iOtherChan).sfreq] ~= hdr.signal(iOtherChan(1)).sfreq);
% iErrChan = setdiff(iErrChan, iAnnotChans);
% if ~isempty(iErrChan)
%     error('Files with mixed sampling rates are not supported yet.');
% end
% Detect interrupted signals (time non-linear)
hdr.interrupted = ischar(hdr.unknown1) && (length(hdr.unknown1) >= 5) && isequal(hdr.unknown1(1:5), 'EDF+D');
if hdr.interrupted
    if ImportOptions.DisplayMessages
        [res, isCancel] = java_dialog('question', ...
            ['Interrupted EDF file ("EDF+D") detected. It is recommended to convert it' 10 ...
            'to a continuous ("EDF+C") file first. Do you want to continue reading this' 10 ...
            'file as continuous and attempt to fix the timing of event markers?' 10 ...
            'NOTE: This may not work as intended, use at your own risk!']);
        hdr.fixinterrupted = ~isCancel && strcmpi(res, 'yes');
    else
        hdr.fixinterrupted = 1;
    end
    if ~hdr.fixinterrupted
        warning(['Interrupted EDF file ("EDF+D"): requires conversion to "EDF+C". ' 10 ...
             'Brainstorm will read this file as a continuous file ("EDF+C"), the timing of the samples after the first discontinuity will be wrong.' 10 ...
             'This may not cause any major problem unless there are time markers in the file, they might be inaccurate in all the segments >= 2.']);
    end
else
    hdr.fixinterrupted = 0;
end

%% ===== CREATE BRAINSTORM SFILE STRUCTURE =====
% Initialize returned file structure
sFile = db_template('sfile');
% Add information read from header
sFile.byteorder  = 'l';
sFile.filename   = DataFile;
if (uint8(hdr.version(1)) == uint8(255))
    sFile.format = 'EEG-BDF';
    sFile.device = 'BDF';
else
    sFile.format = 'EEG-EDF';
    sFile.device = 'EDF';
end
sFile.header = hdr;
% Comment: short filename
[tmp__, sFile.comment, tmp__] = bst_fileparts(DataFile);
% No info on bad channels
sFile.channelflag = ones(hdr.nsignal,1);
% Acquisition date
sFile.acq_date = str_date(hdr.startdate);

%% ===== PROCESS CHANNEL NAMES/TYPES =====
% Try to split the channel names in "TYPE NAME"
SplitType = repmat({''}, 1, hdr.nsignal);
SplitName = repmat({''}, 1, hdr.nsignal);
for i = 1:hdr.nsignal
    if ~isempty(hdr.signal(i).label)
        % Removing trailing dots (eg. "Fc5." instead of "FC5", as in: https://www.physionet.org/pn4/eegmmidb/)
        if (hdr.signal(i).label(end) == '.') && (length(hdr.signal(i).label) > 1)
            hdr.signal(i).label(end) = [];
            if (hdr.signal(i).label(end) == '.') && (length(hdr.signal(i).label) > 1)
                hdr.signal(i).label(end) = [];
                if (hdr.signal(i).label(end) == '.') && (length(hdr.signal(i).label) > 1)
                    hdr.signal(i).label(end) = [];
                end
            end
        end
        % Remove extra spaces
        signalLabel = strrep(hdr.signal(i).label, ' - ', '-');
        % Find space chars (label format "Type Name")
        iSpace = find(signalLabel == ' ');
        % Only if there is one space only, and what is after the space is not only made of numbers
        if (length(iSpace) == 1) && (iSpace >= 3) && ~all(ismember(signalLabel(iSpace+1:end), '0123456789'))
            SplitName{i} = signalLabel(iSpace+1:end);
            SplitType{i} = signalLabel(1:iSpace-1);
        % Accept also 2 spaces
        elseif (length(iSpace) == 2) && (iSpace(1) >= 3) && ~all(ismember(signalLabel(iSpace(1)+1:end), ' 0123456789'))
            SplitName{i} = strrep(signalLabel(iSpace(1)+1:end), ' ', '_');
            SplitType{i} = signalLabel(1:iSpace(1)-1);
        end
    else
        hdr.signal(i).label = sprintf('E%d', i);
    end
end
% Remove the classification if it makes some names non unique
uniqueNames = unique(SplitName);
for i = 1:length(uniqueNames)
    if ~isempty(uniqueNames{i})
        iName = find(strcmpi(SplitName, uniqueNames{i}));
        if (length(iName) > 1)
            [SplitName{iName}] = deal('');
            [SplitType{iName}] = deal('');
        end
    end
end

%% ===== CREATE EMPTY CHANNEL FILE =====
ChannelMat = db_template('channelmat');
ChannelMat.Comment = [sFile.device ' channels'];
ChannelMat.Channel = repmat(db_template('channeldesc'), [1, hdr.nsignal]);
chRef = {};
% For each channel
for i = 1:hdr.nsignal
    % If is the annotation channel
    if ~isempty(iAnnotChans) && ismember(i, iAnnotChans)
        ChannelMat.Channel(i).Type = 'EDF';
        ChannelMat.Channel(i).Name = 'Annotations';
    elseif ~isempty(iStatusChan) && (i == iStatusChan)
        ChannelMat.Channel(i).Type = 'BDF';
        ChannelMat.Channel(i).Name = 'Status';
    % Regular channels
    else
        % If there is a pair name/type already detected
        if ~isempty(SplitName{i}) && ~isempty(SplitType{i})
            ChannelMat.Channel(i).Name = SplitName{i};
            ChannelMat.Channel(i).Type = SplitType{i};
        else
            % Channel name
            ChannelMat.Channel(i).Name = hdr.signal(i).label(hdr.signal(i).label ~= ' ');
            % Channel type
            if ~isempty(hdr.signal(i).type)
                if (length(hdr.signal(i).type) == 3)
                    ChannelMat.Channel(i).Type = hdr.signal(i).type(hdr.signal(i).type ~= ' ');
                elseif isequal(hdr.signal(i).type, 'Active Electrode') || isequal(hdr.signal(i).type, 'AgAgCl electrode')
                    ChannelMat.Channel(i).Type = 'EEG';
                else
                    ChannelMat.Channel(i).Type = 'Misc';
                end
            else
                ChannelMat.Channel(i).Type = 'EEG';
            end
        end
        % Extract reference name (at the end of the channel name, separated with a "-", eg. "-REF")
        iDash = find(ChannelMat.Channel(i).Name == '-');
        if ~isempty(iDash) && (iDash(end) < length(ChannelMat.Channel(i).Name))
            chRef{end+1} = ChannelMat.Channel(i).Name(iDash(end):end);
        end
    end
    ChannelMat.Channel(i).Loc     = [0; 0; 0];
    ChannelMat.Channel(i).Orient  = [];
    ChannelMat.Channel(i).Weight  = 1;
    % ChannelMat.Channel(i).Comment = hdr.signal(i).type;
end

% If the same reference is indicated for all the channels: remove it
if (length(chRef) >= 2)
    % Get the shortest reference tag
    lenRef = cellfun(@length, chRef);
    minLen = min(lenRef);
    % Check if all the ref names are equal (up to the max length - some might be cut because the channel name is too long)
    if all(cellfun(@(c)strcmpi(c(1:minLen), chRef{1}(1:minLen)), chRef))
        % Remove the reference tag from all the channel names
        for i = 1:length(ChannelMat.Channel)
            ChannelMat.Channel(i).Name = strrep(ChannelMat.Channel(i).Name, chRef{1}, '');
            ChannelMat.Channel(i).Name = strrep(ChannelMat.Channel(i).Name, chRef{1}(1:minLen), '');
        end
    end
end

% If there are only "Misc" and no "EEG" channels: rename to "EEG"
iMisc = find(strcmpi({ChannelMat.Channel.Type}, 'Misc'));
iEeg  = find(strcmpi({ChannelMat.Channel.Type}, 'EEG'));
if ~isempty(iMisc) && isempty(iEeg)
    [ChannelMat.Channel(iMisc).Type] = deal('EEG');
    iEeg = iMisc;
end

%% ===== DETECT MULTIPLE SAMPLING RATES =====
% Use the first "EEG" channel as the reference sampling rate (or the first channel if no "EEG" channels available)
if ~isempty(iEeg) && ismember(iEeg(1), iOtherChan)
    iChanFreqRef = iEeg(1);
else
    iChanFreqRef = iOtherChan(1);
end
% Mark as bad channels with sampling rates different from EEG
iChanWrongRate = find([sFile.header.signal.sfreq] ~= sFile.header.signal(iChanFreqRef).sfreq);
iChanWrongRate = intersect(iChanWrongRate, iOtherChan);
if ~isempty(iChanWrongRate)
    sFile.channelflag(iChanWrongRate) = -1;
    disp([sprintf('EDF> WARNING: Excluding channels with sampling rates other than %.3f Hz : ', hdr.signal(iChanFreqRef).sfreq), sprintf('%s ', ChannelMat.Channel(iChanWrongRate).Name)]);
end

% Consider that the sampling rate of the file is the sampling rate of the first signal
sFile.prop.sfreq = hdr.signal(iChanFreqRef).sfreq;
sFile.prop.times = [0, hdr.signal(iChanFreqRef).nsamples * hdr.nrec - 1] ./ sFile.prop.sfreq;
sFile.prop.nAvg  = 1;

%% ===== READ EDF ANNOTATION CHANNEL =====
if ~isempty(iEvtChans) % && ~isequal(ImportOptions.EventsMode, 'ignore')
    % Set reading options
    ImportOptions.ImportMode = 'Time';
    ImportOptions.UseSsp     = 0;
    ImportOptions.UseCtfComp = 0;
    % Read EDF annotations
    if strcmpi(sFile.format, 'EEG-EDF')
        evtList = {};
        % By default: the file starts at 0s (definitions useful in case the first records are missing)
        t0_file = 0;
        prev_rec = 0;
        prev_irec = 1;
        % In EDF+, the first annotation channel has epoch time stamps (EDF
        % calls epochs records).  So read all annotation channels per epoch.
        for irec = 1:hdr.nrec
            for ichan = 1:length(iEvtChans)
                bst_progress('text', sprintf('Reading annotations... [%d%%]', round((ichan + (irec-1)*length(iEvtChans))/length(iEvtChans)/hdr.nrec*100)));
                % Sample indices for the current epoch (=record)
                SampleBounds = [irec-1,irec] * sFile.header.signal(iEvtChans(ichan)).nsamples - [0,1];
                % Read record
                F = in_fread(sFile, ChannelMat, 1, SampleBounds, iEvtChans(ichan), ImportOptions);
                % Find all the TALs (Time-stamped Annotations Lists): separated with [20][0]
                % Possible configurations:
                %    Onset[21]Duration[20]Annot1[20]Annot2[20]..AnnotN[20][0]
                %    Onset[20]Annot1[20]Annot2[20]..AnnotN[20][0]
                %    Onset[20]Annot1[20][0]
                %    Onset[20][20][0]                       : First TAL of a record, with an empty annotation, indicating the offset from the beginning of the file
                %    Onset[20][20]Annot2[20]..AnnotN[20][0] : First TAL + extra annotations
                iSeparator = [-1, find((F(1:end-1) == 20) & (F(2:end) == 0))];
                % Loop on all the TALs
                for iAnnot = 1:length(iSeparator)-1
                    % Get annotation
                    strAnnot = char(F(iSeparator(iAnnot)+2:iSeparator(iAnnot+1)-1));
                    if isempty(strAnnot) || (strAnnot(1) == 20)
                        continue;
                    end
                    % Split in blocks with [20]
                    splitAnnot = str_split(strAnnot, 20, 0);
                    % The first TAL in a record should be indicating the timing of the first sample of the record (empty annotation)
                    if (iAnnot == 1) && ((length(splitAnnot) == 1) || isempty(splitAnnot{2})) && ~any(splitAnnot{1} == 21)
                        % Ignore if this is not the first channel
                        if (ichan > 1)
                            continue;
                        end
                        % Get record time stamp
                        t0_rec = str2double(splitAnnot{1});
                        if isempty(t0_rec) || isnan(t0_rec)
                            continue;
                        end
                        if (irec == 1)
                            t0_file = t0_rec;
                        % Find discontinuities equal or larger than 2 samples (1 is expected)
                        elseif abs(t0_rec - prev_rec - (irec - prev_irec) * hdr.reclen) >= (2 / sFile.prop.sfreq)
                            % Brainstorm fills partial/interrupted records with zeros
                            bstTime = prev_rec + hdr.reclen;
                            timeDiff = bstTime - t0_rec;
                            % If we want to fix timing, apply skip to initial timestamp
                            if hdr.fixinterrupted
                                t0_file = t0_file - timeDiff;
                            end
                            % Warn user of discontinuity
                            if timeDiff > 0
                                expectMsg = 'blank data';
                            else
                                expectMsg = 'skipped data';
                            end
                            startTime = min(t0_rec - t0_file - [0, timeDiff]); % before and after t0_file adjustment
                            endTime  = max(t0_rec - t0_file - [0, timeDiff]);
                            fprintf('WARNING: Found discontinuity between %.3fs and %.3fs (%d samples), expect %s in between.\n', startTime, endTime, round(abs(timeDiff).*sFile.prop.sfreq), expectMsg);
                            % Create event for users information
                            if timeDiff < 0
                                endTime = startTime; % no extent in this case, there is skipped time.
                            end
                            evtList(end+1,:) = {'EDF+D Discontinuity', [startTime; endTime]};
                        end
                        prev_rec = t0_rec;
                        prev_irec = irec;
                        % If there are extra annotations for this TAL: Create one event for each label in the TAL
                        for iLabel = 3:length(splitAnnot)
                            evtList(end+1,:) = {splitAnnot{iLabel}, (t0_rec-t0_file) + [0;0]};
                        end
                    % Regular TAL: indicating a marker
                    else
                        % Split time in onset/duration
                        splitTime = str_split(splitAnnot{1}, 21);
                        % Get time
                        t = str2double(splitTime{1});
                        if isempty(t) || isnan(t)
                            continue;
                        end
                        % Get duration
                        if (length(splitTime) > 1)
                            duration = str2double(splitTime{2});
                            % Exclude 1-sample long events
                            if isempty(duration) || isnan(duration) || (round(duration .* sFile.prop.sfreq) <= 1)
                                duration = 0;
                            end
                        else
                            duration = 0;
                        end
                        % Unnamed event
                        if (length(splitAnnot) == 1) || isempty(splitAnnot{2})
                            splitAnnot{2} = 'Unnamed';
                        end
                        % Create one event for each label in the TAL
                        for iLabel = 2:length(splitAnnot)
                            evtList(end+1,:) = {splitAnnot{iLabel}, (t-t0_file) + [0;duration]};
                        end
                    end
                end
            end
        end

        % If there are events: create a create an events structure
        if ~isempty(evtList)
            % Initialize events list
            sFile.events = repmat(db_template('event'), 0);
            % Events list
            [uniqueEvt, iUnique] = unique(evtList(:,1));
            uniqueEvt = evtList(sort(iUnique),1);
            % Build events list
            for iEvt = 1:length(uniqueEvt)
                % Find all the occurrences of this event
                iOcc = find(strcmpi(uniqueEvt{iEvt}, evtList(:,1)));
                % Concatenate all times
                t = [evtList{iOcc,2}];
                % If second row is equal to the first one (no extended events): delete it
                if all(t(1,:) == t(2,:))
                    t = t(1,:);
                end
                % Set event
                sFile.events(iEvt).label    = strtrim(uniqueEvt{iEvt});
                sFile.events(iEvt).times    = round(t .* sFile.prop.sfreq) ./ sFile.prop.sfreq;
                sFile.events(iEvt).epochs   = 1 + 0*t(1,:);
                sFile.events(iEvt).select   = 1;
                sFile.events(iEvt).channels = [];
                sFile.events(iEvt).notes    = [];
            end
        end

    % BDF Status line
    elseif strcmpi(sFile.format, 'EEG-BDF') && ~strcmpi(ImportOptions.EventsTrackMode, 'ignore')
        % Ask how to read the events
        [events, ImportOptions.EventsTrackMode] = process_evt_read('Compute', sFile, ChannelMat, ChannelMat.Channel(iEvtChans).Name, ImportOptions.EventsTrackMode);
        if isequal(events, -1)
            sFile = [];
            ChannelMat = [];
            return;
        end
        % Report the events in the file structure
        sFile.events = events;
        % Remove the 'Status: ' string in front of the events
        for i = 1:length(sFile.events)
            sFile.events(i).label = strrep(sFile.events(i).label, 'Status: ', '');
        end
        % Group events by time
        % sFile.events = process_evt_grouptime('Compute', sFile.events);
    end
end

author: brainstorm-tools


text:
if ~exist('namecore','var')
    namecore = 'D10ptBMrandPSFgauss_StretchCoef2_N10_offset100_nc';
end
chall=2*ones(length(iteration),length(ncvec),(ncvec(end)-1)^2);
for ii=1:length(ncvec)
    for jj=1:length(iteration)
%         r=load (['D10ptBMrandPSFgauss_StretchCoef2_N10_offset100_nc' num2str(ncvec(ii)) '/results_iter' num2str(iteration(jj)) '.mat']);
        r=load ([namecore num2str(ncvec(ii)) '/results_iter' num2str(iteration(jj)) '.mat']);
        if isfield(r.res,'a')
            r.res.h=r.res.a;
        end
        %     r=load (['results_updates_nmfclassic_nc' num2str(ncvec(ii)) '/results_iter' num2str(iteration) '.mat']);
        if isfield(r.peval,'data_dir')
            d = load(['~/' r.peval.data_path '/' r.peval.data_dir '/' r.peval.data_file]);
        elseif isfield (r.peval,'data_file')
            d=load(['~/' peval.data_path '/' peval.data_file]);
        else
            if exist('dpixc', 'var')
                d.dpixc = dpixc;
            else
                d = load ('dpixc');
            end
        end
        resid = (r.res.w*r.res.h - reshape(d.dpixc, r.peval.nx*r.peval.ny, r.peval.nt));
        % pixels correlations:
        cp = (corrcoef(resid'));
        cpmax(jj,ii) = max(cp(cp<1));
        cpmin(jj,ii) = min(cp(cp<1));
        cpmean(jj,ii) = mean(cp(cp<1));
        cpmeanpos(jj,ii) = mean(cp(and(cp>0,cp<1)));
        cpmeanneg(jj,ii) = mean(cp(cp<0));
        % intensity correlations:
        if size(r.res.h,1)>2
            ch = corr(r.res.h(1:end-1,:)');
            chall(jj,ii,1:numel(ch))=ch(:);
            chmax(jj,ii) = max(ch(ch<1));
            %     chmin(ii) = min(min(corr(r.res.h(1:end-1,:)')-eye(size(r.res.h,1)-1)));
            chmin(jj,ii) = min(ch(:));

            chmean(jj,ii)=mean(ch(ch<1));
            chmeanpos(jj,ii)=mean(ch(and(ch<1,ch>0)));
            chmeanneg(jj,ii)=mean(ch(ch<0));
        end
%         h=hinton(ch);
        ll(jj,ii)=loglikelihoodPoisson(reshape(d.dpixc,r.peval.nx*r.peval.ny,r.peval.nt),r.res.w*r.res.h);
        if isfield(r.res, 'lb')
            lb(jj,ii) = r.res.lb;
        end

    end
end

% figure
% plot(ncvec, chmax.*(abs(chmin)),'s-r')

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [Dlambda,volume,elemSign] = gradbasis3(node,elem)
%% GRADBASIS3 gradient of barycentric basis in 3-D.
%
% [Dlambda,volume,elemSign] = gradbasis3(node,elem) compute gradient of
% barycentric basis and volume of tetrahedron. The array volume is NT by 1
% and volume(t) is the volume of the t-th tetrahedron. Dlambda(1:NT, 1:3,
% 1:4) the first index is the label of tetrahedron, the second is the x-y-z
% coordinate, and the last one is the label of four indices of a tet. For
% example, Dlambda(t,:,1) is the gradient of lambda of the 1st index of the
% t-th tet. The elemSign array taking values 1 or -1 records the
% sign of the signed volume.
%
% See also gradbasis, gradu, gradu3
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

NT = size(elem,1);
face = uint32([elem(:,[2 4 3]);elem(:,[1 3 4]);elem(:,[1 4 2]);elem(:,[1 2 3])]);
v12 = node(face(:,2),:)-node(face(:,1),:);
v13 = node(face(:,3),:)-node(face(:,1),:);
normal = mycross(v12,v13,2);
v12 = v12(3*NT+1:4*NT,:);
v13 = v13(3*NT+1:4*NT,:);
v14 = node(elem(:,4),:)-node(elem(:,1),:);
volume = dot(mycross(v12,v13,2),v14,2)/6;
Dlambda = zeros(NT,3,4);
Dlambda(1:NT,:,1) = normal(1:NT,:)./[6*volume, 6*volume, 6*volume];
Dlambda(1:NT,:,2) = normal(NT+1:2*NT,:)./[6*volume, 6*volume, 6*volume];
Dlambda(1:NT,:,3) = normal(2*NT+1:3*NT,:)./[6*volume, 6*volume, 6*volume];
Dlambda(1:NT,:,4) = normal(3*NT+1:4*NT,:)./[6*volume, 6*volume, 6*volume];

% When the tetrahedron is not positive orientated, we reverse the sign of
% the volume. The sign of Dlambda is always right since signed volume is
% used in the computation.
idx = (volume<0);
volume(idx,:) = -volume(idx,:);
elemSign = ones(NT,1);
elemSign(idx) = -1;

repo: ifem


text:
function DEM = Readopentopo(varargin)

%READOPENTOPO read DEM using the opentopography.org API
%
% Syntax
%
%     DEM = readopentopo(pn,pv,...)
%
% Description
%
%     readopentopo reads DEMs from opentopography.org using the API
%     described on:
%     http://www.opentopography.org/developers
%     The DEM comes in geographic coordinates (WGS84) and should be
%     projected to a projected coordinate system (use reproject2utm) before
%     analysis in TopoToolbox.
%
% Input arguments
%
%     Parameter name values
%     'filename'       provide filename. By default, the function will save
%                      the DEM to a temporary file in the system's temporary
%                      folder.
%     'north'          northern boundary in geographic coordinates (WGS84)
%     'south'          southern boundary
%     'west'           western boundary
%     'east'           eastern boundary
%     'demtype'        The global raster dataset - SRTM GL3 (90m) is
%                      'SRTMGL3', SRTM GL1 (30m) is 'SRTMGL1', SRTM GL1
%                      (Ellipsoidal) is 'SRTMGL1_E', and ALOS World 3D 30m
%                      is 'AW3D30'
%     'deletefile'     'true' or false. True, if file should be deleted
%                      after it was downloaded and added to the workspace.
%
% Output arguments
%
%     DEM            Digital elevation model in geographic coordinates
%                    (GRIDobj)
%
% See also: GRIDobj, websave
%
% Reference: http://www.opentopography.org/developers
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 19. June, 2017

p = inputParser;
addParameter(p,'filename',[tempname '.tif']);
% addParameter(p,'interactive',false);
addParameter(p,'north',37.091337);
addParameter(p,'south',36.738884);
addParameter(p,'west',-120.168457);
addParameter(p,'east',-119.465576);
addParameter(p,'demtype','SRTMGL3');
addParameter(p,'deletefile',true);
parse(p,varargin{:});

demtype = validatestring(p.Results.demtype,{'SRTMGL3','SRTMGL1','SRTMGL1_E','AW3D30'},'readopentopo');

url = 'http://opentopo.sdsc.edu/otr/getdem';

% create output file
f = fullfile(p.Results.filename);

% save to drive
options = weboptions('Timeout',inf);

west = p.Results.west;
east = p.Results.east;
south = p.Results.south;
north = p.Results.north;

% if any([isempty(west) isempty(east) isempty(south) isempty(north)]) || p.Results.interactive;
%     wm = webmap;
%     % get dialog box
%     messagetext = ['Zoom and resize the webmap window to choose DEM extent. ' ...
%                          'Click the close button when you''re done.'];
%     d = waitdialog(messagetext);
%     uiwait(d);
%     [latlim,lonlim] = wmlimits(wm);
%     west = lonlim(1);
%     east = lonlim(2);
%     south = latlim(1);
%     north = latlim(2);
% end

websave(f,url,'west',west,...
              'east',east,...
              'north',north,...
              'south',south,...
              'outputFormat', 'GTiff', ...
              'demtype', demtype, ...
              options);

DEM      = GRIDobj(f);
DEM.name = demtype;

if p.Results.deletefile
    delete(f);
end
end

% function d = waitdialog(messagetext)
%     d = dialog('Position',[300 300 250 150],'Name','Choose rectangle region',...
%         'WindowStyle','normal');
%
%     txt = uicontrol('Parent',d,...
%                'Style','text',...
%                'Position',[20 80 210 40],...
%                'String',messagetext);
%
%     btn = uicontrol('Parent',d,...
%                'Position',[85 20 70 25],...
%                'String','Close',...
%                'Callback','delete(gcf)');
% end

author: GERSL

repo: Fmask

save_path: github-repos/MATLAB/GERSL-Fmask


text:
function [cache, place] = score_add_to_cache(cache,j,ps,score,scoring_fn)
% [cache place] = score_add_to_cache(cache,j,ps,score,scoring_fn)
%
% j is the son node,
% ps is the list of parents of j, for example [12 5 7],
% score is the score to add for this familly.
% scoring_fn is 'bic' or 'bayesian'.
%
% place = where the entry was add.
%
% example for 2 nodes with cache of size 5 :
%
% cache =
%   5   b        0      0      0  --> number of writing in cache +1 and b==1 iff the cache is full
%   0   0        1   -239.12   1  --> 1st familly in the cache (node 1 without parents) calculate with bic
%   0   0        2   -318.98   1
%   1   0        2   -189.23   2  --> 3rd familly in the cache (node 2 with 1 as parent) calculate with bayesian
%   0   1        1   -251.09   1
% .ps2bool.      j    score  1or2 --> new entry
%   |   |        |      |      |
%   |   |        |      |      |___> 1 for 'bic' or 2 for 'bayesian'
%   |   |        |      |__________> score of the familly
%   |   |        |_________________> son node of the familly
%   |   |__________________________> ==1 iff node 2 is parent of son node
%   |______________________________> ==1 iff node 1 is parent of son node
%
% If the cache is FULL then the new place is RanDoMly choose.
%
% francois.olivier.c.h@gmail.com

N=size(cache,2)-3;
place=0;

if ~isempty(find(ps==j))
  disp('This is a cyclic entry, nothing was done.');
elseif j>N | j<0
  disp('This entry is not valid, nothing was done.');
else

  switch scoring_fn
    case 'bic',
      fn=1;
    case 'bayesian',
      fn=2;
    otherwise,
      fn=3;
      %error(['unrecognized scoring fn ' scoring_fn]);
  end
  L=size(cache,1);

  if cache(1,2)==0
    place=cache(1,1);
  else
    place=ceil(rand(1)*(L-1))+1;
  end

  cache(place,:)=0
  cache(place,ps)=1;
  cache(place,N+1)=j;
  cache(place,N+2)=score;
  cache(place,N+3)=fn;
  cache(1,1)=place+1;
  if place==L | cache(1,2)~=0
    cache(1,2)=1
  end

end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function bool = is_array(x)
    % Check if x is regular dense array.
    bool = isa(x, 'double') || isa(x, 'single') || isa(x, 'gpuArray');
end

author: oseledets

repo: TT-Toolbox