text:
%% Permutation test example
%
% A simple example of running a permutation test to determine the
% signifance of classification accuracies
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

%% Set the number of permutations
niter=1000;

%% Define dataset, classifier, partitioner
config=cosmo_config();
data_path=fullfile(config.tutorial_data_path,'ak6','s01');

data_fn=fullfile(data_path,'glm_T_stats_perrun.nii');
mask_fn=fullfile(data_path,'vt_mask.nii');
ds=cosmo_fmri_dataset(data_fn,'mask',mask_fn,...
                        'targets',repmat(1:6,1,10),...
                        'chunks',floor(((1:60)-1)/6)+1);

% remove constant features
ds=cosmo_remove_useless_data(ds);

% Only consider four classes (otherwise the classifier does extremily well)
ds=cosmo_slice(ds,ds.sa.targets<=4);

classifier=@cosmo_classify_nn;

% for more speed, just do odd-even partitioning
partitions=cosmo_oddeven_partitioner(ds);

%% compute classification accuracy of the original data
[pred, acc]=cosmo_crossvalidate(ds, classifier, partitions);

%% prepare for permutations
acc0=zeros(niter,1); % allocate space for permuted accuracies
ds0=ds; % make a copy of the dataset

%% for _niter_ iterations, reshuffle the labels and compute accuracy
% Use the helper function cosmo_randomize_targets
% >@@>
for k=1:niter
    ds0.sa.targets=cosmo_randomize_targets(ds);
    [foo, acc0(k)]=cosmo_crossvalidate(ds0, classifier, partitions);
end
% <@@<

p=sum(acc<acc0)/niter;
fprintf('%d permutations: accuracy=%.3f, p=%.4f\n', niter, acc, p);

bins=0:10/niter:1;
h=histc(acc0,bins);
bar(bins,h)
hold on
line([acc acc],[0,max(h)])
hold off
title(sprintf('acc=%.3f',acc))


text:
function y = logDirichlet(X, a)
% Compute log pdf of a Dirichlet distribution.
% Input:
%   X: d x n data matrix, each column sums to one (sum(X,1)==ones(1,n) && X>=0)
%   a: d x k parameter of Dirichlet
%   y: k x n probability density
% Output:
%   y: k x n probability density in logrithm scale y=log p(x)
% Written by Mo Chen (sth4nth@gmail.com).
X = bsxfun(@times,X,1./sum(X,1));
if size(a,1) == 1
    a = repmat(a,size(X,1),1);
end
c = gammaln(sum(a,1))-sum(gammaln(a),1);
g = (a-1)'*log(X);
y = bsxfun(@plus,g,c');

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function [ids, dis]= yael_nn(v, q, k, distype)
    assert( nargin<4 || distype==2 );
    if nargin<3, k= 1; end
    assert(k<=size(v,2));

    ids= zeros(k, size(q,2), 'int32');
    dis= zeros(k, size(q,2), 'single');

    for iVec= 1:size(q,2)
        ds= sum( bsxfun(@minus, v, q(:,iVec)).^2, 1 );
        [ds, inds]= sort(ds);
        dis(:,iVec)= ds(1:k);
        ids(:,iVec)= inds(1:k);
    end
end

author: Relja

repo: netvlad

save_path: github-repos/MATLAB/Relja-netvlad


text:
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%              COMMON ROBOT CONFIGURATION PARAMETERS                      %
%                                                                         %
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Init simulator core physics paramaters
physics_config.GRAVITY_ACC = [0,0,-9.81];
physics_config.TIME_STEP = Config.tStepSim;

% Robot configuration for WBToolbox
WBTConfigRobotSim = WBToolbox.Configuration;
WBTConfigRobotSim.RobotName = 'icubSim';
WBTConfigRobotSim.UrdfFile = 'model.urdf';
WBTConfigRobotSim.LocalName = 'WBT';
WBTConfigRobotSim.GravityVector = physics_config.GRAVITY_ACC;

% Controlboards and joints list. Each joint is associated to the corresponding controlboard
WBTConfigRobotSim.ControlBoardsNames = {'torso','left_arm','right_arm','left_leg','right_leg'}; %,'head'};
WBTConfigRobotSim.ControlledJoints = [];
numOfJointsForEachControlboard = [];

ControlBoards = struct();
ControlBoards.(WBTConfigRobotSim.ControlBoardsNames{1}) = {'torso_pitch','torso_roll','torso_yaw'};
ControlBoards.(WBTConfigRobotSim.ControlBoardsNames{2}) = {'l_shoulder_pitch','l_shoulder_roll','l_shoulder_yaw','l_elbow'};
ControlBoards.(WBTConfigRobotSim.ControlBoardsNames{3}) = {'r_shoulder_pitch','r_shoulder_roll','r_shoulder_yaw','r_elbow'};
ControlBoards.(WBTConfigRobotSim.ControlBoardsNames{4}) = {'l_hip_pitch','l_hip_roll','l_hip_yaw','l_knee','l_ankle_pitch','l_ankle_roll'};
ControlBoards.(WBTConfigRobotSim.ControlBoardsNames{5}) = {'r_hip_pitch','r_hip_roll','r_hip_yaw','r_knee','r_ankle_pitch','r_ankle_roll'};
% ControlBoards.(WBTConfigRobotSim.ControlBoardsNames{6}) = {'neck_pitch','neck_roll','neck_yaw'};

for n = 1:length(WBTConfigRobotSim.ControlBoardsNames)
    WBTConfigRobotSim.ControlledJoints = [WBTConfigRobotSim.ControlledJoints, ControlBoards.(WBTConfigRobotSim.ControlBoardsNames{n})];
    numOfJointsForEachControlboard = [numOfJointsForEachControlboard; length(ControlBoards.(WBTConfigRobotSim.ControlBoardsNames{n}))];
end

% structure used to configure the Robot class
%
robot_config.jointOrder = WBTConfigRobotSim.ControlledJoints;
robot_config.numOfJointsForEachControlboard = numOfJointsForEachControlboard;

% Note: Since iDynTree 3.0.0, if meshFilePrefix='', the standard iDynTree workflow of locating the
% mesh via the ExternalMesh.getFileLocationOnLocalFileSystem method is used. The iCub model meshes
% file tree is compatible with this workflow.
robot_config.meshFilePrefix = '';
robot_config.fileName = WBTConfigRobotSim.UrdfFile;
robot_config.N_DOF = numel(WBTConfigRobotSim.ControlledJoints);
robot_config.N_DOF_MATRIX = eye(robot_config.N_DOF);

% Initial condition of iCub and for the integrators.
initialConditions.base_position = [0; 0; 0.619];
initialConditions.orientation = diag([-1, -1, 1]);
initialConditions.w_H_b = mwbs.Utils.Rp2H(initialConditions.orientation, initialConditions.base_position);

% joint (inital) position
initialConditions.s = [
    0; 0; 0; ...
    -35.97; 29.97; 0.06; 50.00; ...
    -35.97; 29.97; 0.06; 50.00; ...
    10; 0; 0; -20; -10; 0; ...
    10; 0; 0; -20; -10; 0]*pi/180;

% velocty initial conditions
initialConditions.base_linear_velocity = [0; 0; 0];
initialConditions.base_angular_velocity = [0; 0; 0];
initialConditions.base_pose_dot = [initialConditions.base_linear_velocity; initialConditions.base_angular_velocity];
initialConditions.s_dot = zeros(robot_config.N_DOF, 1);

robot_config.initialConditions = initialConditions;

% Reflected inertia
robot_config.SIMULATE_MOTOR_REFLECTED_INERTIA = true;
INCLUDE_COUPLING = true;

% Robot frames list
FramesSim.BASE = 'root_link';
FramesSim.IMU = 'imu_frame';
FramesSim.LEFT_FOOT = 'l_sole';
FramesSim.RIGHT_FOOT = 'r_sole';
FramesSim.COM = 'com';

robot_config.robotFrames = FramesSim;

% structure used to configure the Contacts class
%

% foot print of the feet (iCub)
vertex = zeros(3, 4);
vertex(:, 1) = [-0.06; 0.04; 0];
vertex(:, 2) = [0.11; 0.04; 0];
vertex(:, 3) = [0.11; -0.035; 0];
vertex(:, 4) = [-0.06; -0.035; 0];

contact_config.foot_print = vertex;
contact_config.total_num_vertices = size(vertex,2)*2;

% friction coefficient for the feet
contact_config.friction_coefficient = 0.1;

%% Motors reflected inertia

% transmission ratio (1/N)
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
Gamma(end-5, end-5)   = 0.0067;
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
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;

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 INCLUDE_COUPLING
    T = blkdiag(T_torso,T_LShoulder,1,T_RShoulder,1,eye(12));
else
    T = eye(robot_config.N_DOF);
end

motorsReflectedInertia = wbc.computeMotorsReflectedInertia(Gamma,T,I_m);

% Joint friction

% === Mapping ===
jointDefaultOrder = {...
    'torso_pitch','torso_roll','torso_yaw', ...
    'l_shoulder_pitch','l_shoulder_roll','l_shoulder_yaw','l_elbow', ...
    'r_shoulder_pitch','r_shoulder_roll','r_shoulder_yaw','r_elbow', ...
    'l_hip_pitch','l_hip_roll','l_hip_yaw','l_knee','l_ankle_pitch','l_ankle_roll', ...
    'r_hip_pitch','r_hip_roll','r_hip_yaw','r_knee','r_ankle_pitch','r_ankle_roll'};

KvmechMappingTorso = 0.2*ones(3,1);
KvmechMappingLeftArm = 0.2*ones(4,1);
KvmechMappingRightArm = 0.2*ones(4,1);
KvmechMappingLeftLeg = [0.2 0.2 0.2 0.2 0.6 0.6]';
KvmechMappingRightLeg = [0.2 0.2 0.2 0.2 0.6 0.6]';

KvmechMapping = containers.Map(...
    jointDefaultOrder, ...
    [
    KvmechMappingTorso
    KvmechMappingLeftArm
    KvmechMappingRightArm
    KvmechMappingLeftLeg
    KvmechMappingRightLeg
    ]);

KvmechMat = diag(cell2mat(KvmechMapping.values(robot_config.jointOrder)));

jointFrictionMat = wbc.computeMotorsReflectedInertia(eye(robot_config.N_DOF),T,KvmechMat);

%% size of the square you see around the robot
visualizerAroundRobot = 1; % mt

clear ControlBoards numOfJointsForEachControlboard FramesSim initialConditions vertex

author: robotology

repo: whole-body-controllers


text:
function [headmodel, sens] = ft_prepare_vol_sens(headmodel, sens, varargin)

% FT_PREPARE_VOL_SENS does some bookkeeping to ensure that the volume
% conductor model and the sensor array are ready for subsequent forward
% leadfield computations. It takes care of some pre-computations that can
% be done efficiently prior to the leadfield calculations.
%
% Use as
%   [headmodel, sens] = ft_prepare_vol_sens(headmodel, sens, ...)
% with input arguments
%   headmodel = structure with volume conductor definition
%   sens      = structure with gradiometer or electrode definition
%
% The headmodel structure represents a volume conductor model of the head,
% its contents depend on the type of model. It is described in more detail
% in FT_DATATYPE_HEADMODEL. The sens structure represents a electrode or
% gradiometer array. It is described in more detail in FT_DATATYPE_SENS.
%
% Additional options should be specified in key-value pairs and can be
%   'channel'  = cell-array with strings (default = 'all')
%
% The detailed behavior of this function depends on whether the input
% consists of EEG or MEG and furthermoree depends on the type of volume
% conductor model:
% - in case of EEG single and concentric sphere models, the electrodes are
%   projected onto the skin surface.
% - in case of EEG boundary element models, the electrodes are projected on
%   the surface and a blilinear interpoaltion matrix from vertices to
%   electrodes is computed.
% - in case of MEG and a localspheres model, a local sphere is determined
%   for each coil in the gradiometer definition.
%  - in case of MEG with a singleshell Nolte model, the volume conduction
%    model is initialized
% In any case channel selection and reordering will be done. The channel
% order returned by this function corresponds to the order in the 'channel'
% option, or if not specified, to the order in the input sensor array.
%
% See also FT_COMPUTE_LEADFIELD, FT_READ_HEADMODEL, FT_READ_SENS

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

if iscell(headmodel) && iscell(sens)
  % this represents combined EEG, ECoG and/or MEG
  for i=1:numel(headmodel)
    [headmodel{i}, sens{i}] = ft_prepare_vol_sens(headmodel{i}, sens{i}, varargin{:});
  end
  return
end

% get the optional input arguments
% fileformat = ft_getopt(varargin, 'fileformat');
channel = ft_getopt(varargin, 'channel', sens.label);   % cell-array with channel labels, default is all

% ensure that the sensor description is up-to-date (Aug 2011)
sens = ft_datatype_sens(sens);

% this is to support volumes saved in mat-files, particularly interpolated
if ischar(headmodel)
  vpath     = fileparts(headmodel);   % remember the path to the file
  headmodel = ft_read_headmodel(headmodel); % replace the filename with the content of the file
end

% ensure that the volume conduction description is up-to-date (Jul 2012)
headmodel = ft_datatype_headmodel(headmodel);

% determine whether the input contains EEG or MEG sensors
iseeg = ft_senstype(sens, 'eeg');
ismeg = ft_senstype(sens, 'meg');

% determine the skin compartment
if ~isfield(headmodel, 'skin_surface')
  if isfield(headmodel, 'bnd')
    headmodel.skin_surface   = find_outermost_boundary(headmodel.bnd);
  elseif isfield(headmodel, 'r') && length(headmodel.r)<=4
    [dum, headmodel.skin_surface] = max(headmodel.r);
  end
end

% determine the inner_skull_surface compartment
if ~isfield(headmodel, 'inner_skull_surface')
  if isfield(headmodel, 'bnd')
    headmodel.inner_skull_surface  = find_innermost_boundary(headmodel.bnd);
  elseif isfield(headmodel, 'r') && length(headmodel.r)<=4
    [dum, headmodel.inner_skull_surface] = min(headmodel.r);
  end
end

% FT_HEADMODELTYPE to an empty struct won't work further down
if isempty(headmodel)
  headmodel = [];
end

% this makes them easier to recognise
sens.type       = ft_senstype(sens);
headmodel.type  = ft_headmodeltype(headmodel);

if isfield(headmodel, 'unit') && isfield(sens, 'unit') && ~strcmp(headmodel.unit, sens.unit)
  ft_error('inconsistency in the units of the volume conductor and the sensor array');
end

if ismeg && iseeg
  % this is something that could be implemented relatively easily
  ft_error('simultaneous EEG and MEG not yet supported');

elseif ~ismeg && ~iseeg
  ft_error('the input does not look like EEG, nor like MEG');

elseif ismeg

  % always ensure that there is a linear transfer matrix for combining the coils into gradiometers
  if ~isfield(sens, 'tra');
    Nchans = length(sens.label);
    Ncoils = size(sens.coilpos,1);
    if Nchans~=Ncoils
      ft_error('inconsistent number of channels and coils');
    end
    sens.tra = eye(Nchans, Ncoils);
  end

  if ~ft_headmodeltype(headmodel, 'localspheres')
    % select the desired channels from the gradiometer array
    [selchan, selsens] = match_str(channel, sens.label);
    % only keep the desired channels, order them according to the users specification
    try, sens.chantype = sens.chantype(selsens,:); end
    try, sens.chanunit = sens.chanunit(selsens,:); end
    try, sens.chanpos  = sens.chanpos (selsens,:); end
    try, sens.chanori  = sens.chanori (selsens,:); end
    sens.label    = sens.label(selsens);
    sens.tra      = sens.tra(selsens,:);
  else
    % for the localspheres model it is done further down
  end

  % remove the coils that do not contribute to any channel output
  selcoil      = any(sens.tra~=0,1);
  sens.coilpos = sens.coilpos(selcoil,:);
  sens.coilori = sens.coilori(selcoil,:);
  sens.tra     = sens.tra(:,selcoil);

  switch ft_headmodeltype(headmodel)
    case {'infinite' 'infinite_monopole' 'infinite_currentdipole' 'infinite_magneticdipole'}
      % nothing to do

    case 'singlesphere'
      % nothing to do

    case 'concentricspheres'
      % nothing to do

    case 'neuromag'
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % if the forward model is computed using the external Neuromag toolbox,
      % we have to add a selection of the channels so that the channels
      % in the forward model correspond with those in the data.
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      [selchan, selsens] = match_str(channel, sens.label);
      headmodel.chansel = selsens;

    case 'localspheres'
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % If the volume conduction model consists of multiple spheres then we
      % have to match the channels in the gradiometer array and the volume
      % conduction model.
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      % the initial localspheres volume conductor has a local sphere per
      % channel, whereas it should have a local sphere for each coil
      if size(headmodel.r,1)==size(sens.coilpos,1) && ~isfield(headmodel, 'label')
        % it appears that each coil already has a sphere, which suggests
        % that the volume conductor already has been prepared to match the
        % sensor array
        return
      elseif size(headmodel.r,1)==size(sens.coilpos,1) && isfield(headmodel, 'label')
        if ~isequal(headmodel.label(:), sens.label(:))
          % if only the order is different, it would be possible to reorder them
          ft_error('the coils in the volume conduction model do not correspond to the sensor array');
        else
          % the coil-specific spheres in the volume conductor should not have a label
          % because the label is already specified for the coils in the
          % sensor array
          headmodel = rmfield(headmodel, 'label');
        end
        return
      end

      % the CTF way of representing the headmodel is one-sphere-per-channel
      % whereas the FieldTrip way of doing the forward computation is one-sphere-per-coil
      Nchans   = size(sens.tra,1);
      Ncoils   = size(sens.tra,2);
      Nspheres = size(headmodel.label);

      if isfield(headmodel, 'orig')
        % these are present in a CTF *.hdm file
        singlesphere.o(1,1) = headmodel.orig.MEG_Sphere.ORIGIN_X;
        singlesphere.o(1,2) = headmodel.orig.MEG_Sphere.ORIGIN_Y;
        singlesphere.o(1,3) = headmodel.orig.MEG_Sphere.ORIGIN_Z;
        singlesphere.r      = headmodel.orig.MEG_Sphere.RADIUS;
        % ensure consistent units
        singlesphere = ft_convert_units(singlesphere, headmodel.unit);
        % determine the channels that do not have a corresponding sphere
        % and use the globally fitted single sphere for those
        missing = setdiff(sens.label, headmodel.label);
        if ~isempty(missing)
          ft_warning('using the global fitted single sphere for %d channels that do not have a local sphere', length(missing));
        end
        for i=1:length(missing)
          headmodel.label(end+1) = missing(i);
          headmodel.r(end+1,:)   = singlesphere.r;
          headmodel.o(end+1,:)   = singlesphere.o;
        end
      end

      % make a new structure that only holds the local spheres, one per coil
      localspheres = [];
      localspheres.type = headmodel.type;
      localspheres.unit = headmodel.unit;

      % for each coil in the MEG helmet, determine the corresponding channel and from that the corresponding local sphere
      for i=1:Ncoils
        coilindex = find(sens.tra(:,i)~=0); % to which channel does this coil belong
        if length(coilindex)>1
          % this indicates that there are multiple channels to which this coil contributes,
          % which happens if the sensor array represents a synthetic higher-order gradient.
          [dum, coilindex] = max(abs(sens.tra(:,i)));
        end

        coillabel = sens.label{coilindex};               % what is the label of this channel
        chanindex = find(strcmp(coillabel, headmodel.label));  % what is the index of this channel in the list of local spheres
        localspheres.r(i,:) = headmodel.r(chanindex);
        localspheres.o(i,:) = headmodel.o(chanindex,:);
      end
      headmodel = localspheres;

      % finally do the selection of channels and coils
      % order them according to the users specification
      [selchan, selsens] = match_str(channel, sens.label);

      % first only modify the linear combination of coils into channels
      try, sens.chantype = sens.chantype(selsens,:); end
      try, sens.chanunit = sens.chanunit(selsens,:); end
      try, sens.chanpos  = sens.chanpos (selsens,:); end
      try, sens.chanori  = sens.chanori (selsens,:); end
      sens.label   = sens.label(selsens);
      sens.tra     = sens.tra(selsens,:);
      % subsequently remove the coils that do not contribute to any sensor output
      selcoil      = find(sum(sens.tra,1)~=0);
      sens.coilpos = sens.coilpos(selcoil,:);
      sens.coilori = sens.coilori(selcoil,:);
      sens.tra     = sens.tra(:,selcoil);
      % make the same selection of coils in the localspheres model
      headmodel.r = headmodel.r(selcoil);
      headmodel.o = headmodel.o(selcoil,:);

    case 'singleshell'
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % if the forward model is computed using the code from Guido Nolte, we
      % have to initialize the volume model using the gradiometer coil
      % locations
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      % compute the surface normals for each vertex point
      if ~isfield(headmodel.bnd, 'nrm')
        fprintf('computing surface normals\n');
        headmodel.bnd.nrm = normals(headmodel.bnd.pos, headmodel.bnd.tri);
      end

      % estimate center and radius
      [center, radius] = fitsphere(headmodel.bnd.pos);

      % order of spherical spherical harmonics, for 'real' realistic volume conductors order=10 is o.k
      if isfield(headmodel, 'order')
        order = headmodel.order;
      else
        order = 10;
      end

      % initialize the forward calculation (only if  coils are available)
      if size(sens.coilpos,1)>0 && ~isfield(headmodel, 'forwpar')
        s = ft_scalingfactor(headmodel.unit, 'cm');
        headmodel.forwpar = meg_ini([s*headmodel.bnd.pos headmodel.bnd.nrm], s*center', order, [s*sens.coilpos sens.coilori]);
        headmodel.forwpar.scale = s;
      end

    case  'openmeeg'
        % don't do anything, h2em or h2mm generated later in ft_prepare_leadfield

    case 'simbio'
      ft_error('MEG not yet supported with simbio');

    otherwise
      ft_error('unsupported volume conductor model for MEG');
  end

elseif iseeg

  % the electrodes are used, the channel positions are not relevant any more
  % channel positinos need to be recomputed after projecting the electrodes on the skin
  if isfield(sens, 'chanpos'); sens = rmfield(sens, 'chanpos'); end

  % select the desired channels from the electrode array
  % order them according to the users specification
  [selchan, selsens] = match_str(channel, sens.label);
  Nchans = length(sens.label);

  sens.label     = sens.label(selsens);
  try, sens.chantype  = sens.chantype(selsens); end
  try, sens.chanunit  = sens.chanunit(selsens); end

  if isfield(sens, 'tra')
    % first only modify the linear combination of electrodes into channels
    sens.tra     = sens.tra(selsens,:);
    % subsequently remove the electrodes that do not contribute to any channel output
    selelec      = any(sens.tra~=0,1);
    sens.elecpos = sens.elecpos(selelec,:);
    sens.tra     = sens.tra(:,selelec);
  else
    % the electrodes and channels are identical
    sens.elecpos = sens.elecpos(selsens,:);
  end

  switch ft_headmodeltype(headmodel)
    case {'infinite' 'infinite_monopole' 'infinite_currentdipole'}
      % nothing to do

    case {'halfspace', 'halfspace_monopole'}
      % electrodes' all-to-all distances
      numelec = size(sens.elecpos,1);
      ref_el = sens.elecpos(1,:);
      md = dist( (sens.elecpos-repmat(ref_el,[numelec 1]))' );
      % take the min distance as reference
      md = min(md(1,2:end));
      pos = sens.elecpos;
      % scan the electrodes and reposition the ones which are in the
      % wrong halfspace (projected on the plane)... if not too far away!
      for i=1:size(pos,1)
        P = pos(i,:);
        is_in_empty = acos(dot(headmodel.ori,(P-headmodel.pos)./norm(P-headmodel.pos))) < pi/2;
        if is_in_empty
          dPplane = abs(dot(headmodel.ori, headmodel.pos-P, 2));
          if dPplane>md
            ft_error('Some electrodes are too distant from the plane: consider repositioning them')
          else
            % project point on plane
            Ppr = pointproj(P,[headmodel.pos headmodel.ori]);
            pos(i,:) = Ppr;
          end
        end
      end
      sens.elecpos = pos;

    case {'slab_monopole'}
      % electrodes' all-to-all distances
      numelc  = size(sens.elecpos,1);
      ref_elc = sens.elecpos(1,:);
      md  = dist( (sens.elecpos-repmat(ref_elc,[numelc 1]))' );
      % choose min distance between electrodes
      md  = min(md(1,2:end));
      pos = sens.elecpos;
      % looks for contacts outside the strip which are not too far away
      % and projects them on the nearest plane
      for i=1:size(pos,1)
        P = pos(i,:);
        instrip1 = acos(dot(headmodel.ori1,(P-headmodel.pos1)./norm(P-headmodel.pos1))) > pi/2;
        instrip2 = acos(dot(headmodel.ori2,(P-headmodel.pos2)./norm(P-headmodel.pos2))) > pi/2;
        is_in_empty = ~(instrip1&instrip2);
        if is_in_empty
          dPplane1 = abs(dot(headmodel.ori1, headmodel.pos1-P, 2));
          dPplane2 = abs(dot(headmodel.ori2, headmodel.pos2-P, 2));
          if dPplane1>md && dPplane2>md
            ft_error('Some electrodes are too distant from the planes: consider repositioning them')
          elseif dPplane2>dPplane1
            % project point on nearest plane
            Ppr = pointproj(P,[headmodel.pos1 headmodel.ori1]);
            pos(i,:) = Ppr;
          else
            % project point on nearest plane
            Ppr = pointproj(P,[headmodel.pos2 headmodel.ori2]);
            pos(i,:) = Ppr;
          end
        end
      end
      sens.elecpos = pos;

    case {'singlesphere', 'concentricspheres'}
      % ensure that the electrodes ly on the skin surface
      radius = max(headmodel.r);
      pos    = sens.elecpos;
      if isfield(headmodel, 'o')
        % shift the the centre of the sphere to the origin
        pos(:,1) = pos(:,1) - headmodel.o(1);
        pos(:,2) = pos(:,2) - headmodel.o(2);
        pos(:,3) = pos(:,3) - headmodel.o(3);
      end
      distance = sqrt(sum(pos.^2,2)); % to the center of the sphere
      if any((abs(distance-radius)/radius)>0.005)
        ft_warning('electrodes do not lie on skin surface -> using radial projection')
      end
      pos = pos * radius ./ [distance distance distance];
      if isfield(headmodel, 'o')
        % shift the center back to the original location
        pos(:,1) = pos(:,1) + headmodel.o(1);
        pos(:,2) = pos(:,2) + headmodel.o(2);
        pos(:,3) = pos(:,3) + headmodel.o(3);
      end
      sens.elecpos = pos;

    case {'bem', 'dipoli', 'asa', 'bemcp'}
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % do postprocessing of volume and electrodes in case of BEM model
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      % project the electrodes on the skin and determine the bilinear interpolation matrix
      if ~isfield(headmodel, 'tra') && (isfield(headmodel, 'mat') && ~isempty(headmodel.mat))
        % determine boundary corresponding with skin and inner_skull_surface
        if ~isfield(headmodel, 'skin_surface')
          headmodel.skin_surface = find_outermost_boundary(headmodel.bnd);
          fprintf('determining skin compartment (%d)\n', headmodel.skin_surface);
        end
        if ~isfield(headmodel, 'source')
          headmodel.source = find_innermost_boundary(headmodel.bnd);
          fprintf('determining source compartment (%d)\n', headmodel.source);
        end
        if size(headmodel.mat,1)~=size(headmodel.mat,2) && size(headmodel.mat,1)==length(sens.elecpos)
          fprintf('electrode transfer and system matrix were already combined\n');
        else
          fprintf('projecting electrodes on skin surface\n');
          % compute linear interpolation from triangle vertices towards electrodes
          [el, prj] = project_elec(sens.elecpos, headmodel.bnd(headmodel.skin_surface).pos, headmodel.bnd(headmodel.skin_surface).tri);
          tra       = transfer_elec(headmodel.bnd(headmodel.skin_surface).pos, headmodel.bnd(headmodel.skin_surface).tri, el);

          % replace the original electrode positions by the projected positions
          sens.elecpos = prj;

          if size(headmodel.mat,1)==size(headmodel.bnd(headmodel.skin_surface).pos,1)
            % construct the transfer from only the skin vertices towards electrodes
            interp = tra;
          else
            % construct the transfer from all vertices (also inner_skull_surface/outer_skull_surface) towards electrodes
            interp = [];
            for i=1:length(headmodel.bnd)
              if i==headmodel.skin_surface
                interp = [interp, tra];
              else
                interp = [interp, zeros(size(el,1), size(headmodel.bnd(i).pos,1))];
              end
            end
          end

          % incorporate the linear interpolation matrix and the system matrix into one matrix
          % this speeds up the subsequent repeated leadfield computations
          fprintf('combining electrode transfer and system matrix\n');

          % convert to sparse matrix to speed up the subsequent multiplication
          interp  = sparse(interp);
          headmodel.mat = interp * headmodel.mat;
          % ensure that the model potential will be average referenced
          avg = mean(headmodel.mat, 1);
          headmodel.mat = headmodel.mat - repmat(avg, size(headmodel.mat,1), 1);
        end
      end
    case  'openmeeg'
        % don't do anything, h2em or h2mm generated later in ft_prepare_leadfield

    case 'fns'
      if isfield(headmodel,'bnd')
        [el, prj] = project_elec(sens.elecpos, headmodel.bnd.pos, headmodel.bnd.tri);
        sens.tra = transfer_elec(headmodel.bnd.pos, headmodel.bnd.tri, el);
        % replace the original electrode positions by the projected positions
        sens.elecpos = prj;
      end

    case 'simbio'
      % check that the external toolbox is present
      ft_hastoolbox('simbio', 1);

      % extract the outer surface
      bnd = mesh2edge(headmodel);
      for j=1:length(sens.label)
        d = bsxfun(@minus, bnd.pos, sens.elecpos(j,:));
        [d, i] = min(sum(d.^2, 2));
        % replace the position of each electrode by the closest vertex
        sens.elecpos(j,:) = bnd.pos(i,:);
      end

      if (isfield(headmodel,'transfer') && isfield(headmodel,'elec'))
          if all(ismember(sens.label,headmodel.elec.label))
              [sensmember, senslocation] = ismember(sens.label,headmodel.elec.label);
              if (norm(sens.elecpos - headmodel.elec.elecpos(senslocation,:))<1e-8)
                  headmodel.transfer = headmodel.transfer(senslocation,:);
                  headmodel.elec = sens;
              else
                  ft_error('Electrode positions do not fit to the given transfer matrix!');
              end
          else
              ft_error('Transfer matrix does not fit the given set of electrodes!');
          end
      else
          headmodel.transfer = sb_transfer(headmodel,sens);
          headmodel.elec = sens;
      end

    case 'interpolate'
      % this is to allow moving leadfield files
      if ~exist(headmodel.filename{1}, 'file')
        for i = 1:length(headmodel.filename)
          [p, f, x] = fileparts(headmodel.filename{i});
          headmodel.filename{i} = fullfile(vpath, [f x]);
        end
      end

      matchlab = isequal(sens.label, headmodel.sens.label);
      matchpos = isequal(sens.elecpos, headmodel.sens.elecpos);
      matchtra = (~isfield(sens, 'tra') && ~isfield(headmodel.sens, 'tra')) || isequal(sens.tra, headmodel.sens.tra);

      if matchlab && matchpos && matchtra
        % the input sensor array matches precisely with the forward model
        % no further interpolation is needed
      else
        % interpolate the channels in the forward model to the desired channels
        filename = tempname;
        headmodel  = ft_headmodel_interpolate(filename, sens, headmodel);
        % update the sensor array with the one from the volume conductor
        sens = headmodel.sens;
      end % if recomputing interpolation

      % for the leadfield computations the @nifti object is used to map the image data into memory
      ft_hastoolbox('spm8up', 1);
      for i=1:length(headmodel.sens.label)
        % map each of the leadfield files into memory
        headmodel.chan{i} = nifti(headmodel.filename{i});
      end

    otherwise
      ft_error('unsupported volume conductor model for EEG');
  end

  % FIXME this needs careful thought to ensure that the average referencing which is now done here and there, and that the linear interpolation in case of BEM are all dealt with consistently
  % % always ensure that there is a linear transfer matrix for
  % % rereferencing the EEG potential
  % if ~isfield(sens, 'tra');
  %   sens.tra = eye(length(sens.label));
  % end

  % update the channel positions as the electrodes were projected to the skin surface
  [pos, ori, lab] = channelposition(sens);
  [selsens, selpos] = match_str(sens.label, lab);
  sens.chanpos = nan(length(sens.label),3);
  sens.chanpos(selsens,:) = pos(selpos,:);

end % if iseeg or ismeg

if isfield(sens, 'tra')
  if issparse(sens.tra) && size(sens.tra, 1)==1
    % this multiplication would result in a sparse leadfield, which is not what we want
    % the effect can be demonstrated as sparse(1)*rand(1,10), see also http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1169#c7
    sens.tra = full(sens.tra);
  elseif ~issparse(sens.tra) && size(sens.tra, 1)>1
    % the multiplication of the "sensor" leadfield (electrode or coil) with the tra matrix to get the "channel" leadfield
    % is faster for most cases if the pre-multiplying weighting matrix is made sparse
    sens.tra = sparse(sens.tra);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Ppr = pointproj(P,plane)
% projects a point on a plane
% plane(1:3) is a point on the plane
% plane(4:6) is the ori of the plane
Ppr  = [];
ori  = plane(4:6);
line = [P ori];
% get indices of line and plane which are parallel
par = abs(dot(plane(4:6), line(:,4:6), 2))<1e-14;
% difference between origins of plane and line
dp = plane(1:3) - line(:, 1:3);
% Divide only for non parallel vectors (DL)
t = dot(ori(~par,:), dp(~par,:), 2)./dot(ori(~par,:), line(~par,4:6), 2);
% compute coord of intersection point
Ppr(~par, :) = line(~par,1:3) + repmat(t,1,3).*line(~par,4:6);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function serves as a replacement for the dist function in the Neural
% Networks toolbox.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [d] = dist(x)
n = size(x,2);
d = zeros(n,n);
for i=1:n
  for j=(i+1):n
    d(i,j) = sqrt(sum((x(:,i)-x(:,j)).^2));
    d(j,i) = d(i,j);
  end
end

author: spm


text:
function prob_test1555 ( )

%*****************************************************************************80
%
%% TEST1555 tests VON_MISES_CDF, VON_MISES_CDF_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST1555:\n' );
  fprintf ( 1, '  VON_MISES_CDF evaluates the von Mises CDF.\n' );
  fprintf ( 1, '  VON_MISES_CDF_VALUES returns some exact values.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A is the dominant angle;\n' );
  fprintf ( 1, '  B is a measure of spread;\n' );
  fprintf ( 1, '  X is the angle;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      A     B         X   Exact F     Computed F\n' );
  fprintf ( 1, '\n' );
  n_data = 0;

  while ( 1 )

    [ n_data, a, b, x, fx ] = von_mises_cdf_values ( n_data );

    if ( n_data == 0 );
      break
    end

    fx2 = von_mises_cdf ( x, a, b );

    fprintf ( 1, '  %8f  %8f  %8f  %14f  %14f\n', a, b, x, fx, fx2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of the
%   IRB_4600.
%
%   Author: Arturo Gil. Universidad Miguel Hern�ndez de Elche.
%   email: arturo.gil@umh.es date:   09/01/2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

robot.name= 'abb_IRB_4600';

%Path where everything is stored for this robot
robot.path = 'robots/abb/IRB_4600';

robot.DH.theta= '[q(1) q(2)-pi/2 q(3) q(4) q(5) q(6)+pi]';
robot.DH.d='[0.495 0 0 0.960 0 0.135]';
robot.DH.a='[0.175 0.9 0.175 0 0 0]';
robot.DH.alpha= '[-pi/2 0 -pi/2 pi/2 -pi/2 0]';
robot.J=[];

robot.inversekinematic_fn = 'inversekinematic_IRB_4600(robot, T)';

%number of degrees of freedom
robot.DOF = 6;

%rotational: 0, translational: 1
robot.kind=['R' 'R' 'R' 'R' 'R' 'R'];

%minimum and maximum rotation angle in rad
robot.maxangle =[deg2rad(-180) deg2rad(180); %Axis 1, minimum, maximum
                deg2rad(-90) deg2rad(150); %Axis 2, minimum, maximum
                deg2rad(-180) deg2rad(75); %Axis 3
                deg2rad(-400) deg2rad(400); %Axis 4: Unlimited (400� default)
                deg2rad(-125) deg2rad(120); %Axis 5
                deg2rad(-400) deg2rad(400)]; %Axis 6: Unlimited (800� default)

%maximum absolute speed of each joint rad/s or m/s
robot.velmax = [deg2rad(200); %Axis 1, rad/s
                deg2rad(200); %Axis 2, rad/s
                deg2rad(260); %Axis 3, rad/s
                deg2rad(360); %Axis 4, rad/s
                deg2rad(360); %Axis 5, rad/s
                deg2rad(450)];%Axis 6, rad/s

robot.accelmax=robot.velmax/0.1; % 0.1 is here an acceleration time

            % end effectors maximum velocity
robot.linear_velmax = 2.5; %m/s

%base reference system
robot.T0 = eye(4);

%INITIALIZATION OF VARIABLES REQUIRED FOR THE SIMULATION
%position, velocity and acceleration
robot=init_sim_variables(robot);

% GRAPHICS
robot.graphical.has_graphics=1;
robot.graphical.color = [255 20 40]./255;
%for transparency
robot.graphical.draw_transparent=0;
%draw DH systems
robot.graphical.draw_axes=1;
%DH system length and Font size, standard is 1/10. Select 2/20, 3/30 for
%bigger robots
robot.graphical.axes_scale=1;
%adjust for a default view of the robot
robot.axis=[-0.75 0.75 -0.75 0.75 0 1.2];
%read graphics files
robot = read_graphics(robot);

%DYNAMICS
robot.has_dynamics=0;

repo: ARTE


text:
% COVAREP parallel feature formant extraction function
%
% Description
%  This function performs COVAREP_feature_formant_extraction
%  on all .wav files in the specified folder in parallel
%
% Input
%  in_dir [directory path] : Path to directory containing wav files to be analysed
%  options:     Struct with optional arguments (applied to all files in folder).
%
% Output
%  None
%
% 'options' can have the following fields:
%  feature_fs:  Sampling rate of the returned features (default: 0.01ms)
%  gender:      Used for vowelSpace (0: female (default); 1: male; 2: child)
%  features:    The list of features to be extracted (default: all).
%               Possible features are: 'f0', 'F', 'VUV', 'NAQ', 'QOQ',
%               'H1H2', 'PSP', 'MDQ', 'HRF', 'peakSlope', 'Rd', 'Rd_conf',
%               'MCEP', 'HMPDM', 'HMPDD', 'vowelSpace', 'VAD'
%  mcep_order:  Order of the MCEPs (default: 24)
%  hmpdm_order:	Order of the HMPDMs (default: 24)
%  hmpdd_order:	Order of the HMPDDs (default: 12)
%  channel:     If the audio file has multiple channels (default: 1)
%  save_mat:    Save the results as a MAT file (default: false)
%  save_csv:    Save the results as a CSV file (default: true)
%                   If 'save_mat' or 'save_csv' are strings, they are
%                   interpreted as paths to the file to be created.
% Author
%   Florian Helmhold <Florian@Helmhold.de> based on
%   COVAREP_feature_formant_extraction_perfile (Torsten Wörtwein <twoertwe@cs.cmu.edu>)
%   COVAREP_feature_extraction.m (John Kane <kanejo@.tcd.ie>)
%   COVAREP_formant_extraction.m (Stefan Scherer <scherer@ict.usc.edu>)

function COVAREP_feature_formant_extraction(in_dir, options)

%% Default settings
if nargin < 2
    options = struct; % Default empty options struct
end

%% Sanity check options
if isfield(options, 'save_mat') && ~islogical(options.save_mat)
    warning('Non boolean save_mat options will overwrite the output file.')
end

if isfield(options, 'save_csv') && ~islogical(options.save_csv)
    warning('Non boolean save_csv options will overwrite the output file.')
end

% Analysis settings
fileList=dir([in_dir filesep '*.wav']);
N=length(fileList);

if N==0
    disp('No wav files in inputted directory!!!')
end

if(hasParallelComputing)
    parfor n=1:N
        basename=regexp(fileList(n).name,'\.wav','split');
        basename=char(basename(1));
        disp(['Analysing file: ' basename])
        try
            COVAREP_feature_formant_extraction_perfile([in_dir filesep basename '.wav'], options)
            disp([basename ' successfully analysed'])
        catch err
            warning(['An error occurred while analysing ' basename ': ' getReport(err)])
        end
    end
else
    for n=1:N
        basename=regexp(fileList(n).name,'\.wav','split');
        basename=char(basename(1));
        disp(['Analysing file: ' basename])
        try
            COVAREP_feature_formant_extraction_perfile([in_dir filesep basename '.wav'], options)
            disp([basename ' successfully analysed'])
        catch err
            warning(['An error occurred while analysing ' basename ': ' getReport(err)])
        end
    end
end

%% Helper functions

function has_parallel = hasParallelComputing
% get all installed toolbox names
v = ver;
% collect the names in a cell array
[installedToolboxes{1:length(v)}] = deal(v.Name);

% check
has_parallel = all(ismember({'Parallel Computing Toolbox'},installedToolboxes));

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [ subG ] = gsp_subgraph( G,c )
%GSP_SUBGRAPH Create a subgraph from G
%   Usage: [ subG ] = gsp_subgraph( G,c )
%
%   Input parameters:
%       G   : Original graph
%       c   : Node to keep
%
%   Output parameters:
%       subG: Subgraph
%
%   This function create a subgraph from G taking only the node in c.

% Author: Nathanael Perraudin
% Date : Mai 2015

subG.W = G.W(c,c);
subG.N = length(c);
if isfield(G,'type')
    subG.type = ['sub-',G.type];
end

if isfield(G,'directed')
    subG.directed = G.directed;
end

if isfield(G,'lap_type')
    subG.lap_type = G.lap_type;
end

if isfield(G,'coords')
    subG.coords = G.coords(c,:);
end

if isfield(G,'limits')
    subG.limits = G.limits;
end

if isfield(G,'edge_width')
    subG.edge_width = G.edge_width;
end

if isfield(G,'edge_color')
    subG.edge_color = G.edge_color;
end

if isfield(G,'edge_style')
    subG.edge_style = G.edge_style;
end

if isfield(G,'vertex_size')
    subG.vertex_size = G.vertex_size;
end

if isfield(G,'vertex_color')
    subG.vertex_color = G.vertex_color;
end

if isfield(G,'vertex_edge_color')
    subG.vertex_edge_color = G.vertex_edge_color;
end

subG = gsp_graph_default_parameters(subG);

end

repo: gspbox


text:
function imdb = voc_get_database(vocDir, varargin)
opts.seed = 1 ;
opts = vl_argparse(opts, varargin) ;

assert(opts.seed == 1) ; % there is only one split

imdb.imageDir = fullfile(vocDir, 'JPEGImages') ;
imdb.maskDir = fullfile(vocDir, 'Masks'); % doesn't exist
imdb.classes.name={...
  'aeroplane'
  'bicycle'
  'bird'
  'boat'
  'bottle'
  'bus'
  'car'
  'cat'
  'chair'
  'cow'
  'diningtable'
  'dog'
  'horse'
  'motorbike'
  'person'
  'pottedplant'
  'sheep'
  'sofa'
  'train'
  'tvmonitor'};
numClass = length(imdb.classes.name);

% source images
imageFiles = dir(fullfile(imdb.imageDir, '*.jpg'));
imdb.images.name = {imageFiles.name};
numImages = length(imdb.images.name);

imdb.images.label = zeros(numClass, numImages);
imdb.images.label = zeros(1, numImages);
imdb.images.set = zeros(1, numImages);
imdb.images.id = 1:numImages;
imdb.images.vocid = cellfun(@(S) S(1:end-4), imdb.images.name, 'UniformOutput', false);
imageSets = {'train', 'val', 'test'};

% Loop over classes load labels
for c = 1:numClass,
  for s = 1:length(imageSets);
    imageSetPath = fullfile(vocDir, 'ImageSets', 'Main', sprintf('%s_%s.txt', imdb.classes.name{c}, imageSets{s}));
    [gtids,gt]=textread(imageSetPath,'%s %d');
    [membership, loc] = ismember(gtids, imdb.images.vocid);
    assert(all(membership));
    imdb.images.label(c, loc) = gt ;
  end
end

% Loop over images and record the imag sets
for s = 1:length(imageSets),
  imageSetPath = fullfile(vocDir, 'ImageSets', 'Main', sprintf('%s.txt',imageSets{s}));
  gtids = textread(imageSetPath,'%s');
  [membership, loc] = ismember(gtids, imdb.images.vocid);
  assert(all(membership));
  imdb.images.set(loc) = s;
end

% Remove images not part of train, val, test sets
valid = ismember(imdb.images.set, 1:length(imageSets));
imdb.images.name = imdb.images.name(imdb.images.id(valid));
imdb.images.id = 1:numel(imdb.images.name);
imdb.images.label = imdb.images.label(:, valid);
imdb.images.set = imdb.images.set(valid);
imdb.images.vocid = imdb.images.vocid(valid);

% Write out the segments
imdb.segments = imdb.images ;
imdb.segments.imageId = imdb.images.id ;
imdb.segments.difficult = false(1, numel(imdb.segments.id)) ;
% no masks

% make this compatible with the OS imdb
imdb.meta.classes = imdb.classes.name ;
imdb.meta.inUse = true(1,numel(imdb.meta.classes)) ;

author: mcimpoi

repo: deep-fbanks

save_path: github-repos/MATLAB/mcimpoi-deep-fbanks


text:
function g = gibbsKernGradient(kern, x, varargin)

% GIBBSKERNGRADIENT Gradient of GIBBS kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the
% Mark Gibbs's non-stationary
% kernel's parameters. As well as the kernel structure and the
% input positions, the user provides a matrix PARTIAL which gives
% the partial derivatives of the function with respect to the
% relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x : the input locations for which the gradients are being
% computed.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in X.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x1 : the input locations associated with the rows of the
% kernel matrix.
% ARG x2 : the input locations associated with the columns of the
% kernel matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as X1 and the same number of columns
% as X2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% SEEALSO gibbsKernParamInit, kernGradient, gibbsKernDiagGradient, kernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2006, 2009

% KERN

%  pretty(simple(diff((2*(l_i*l_j)/(l_i*l_i + l_j*l_j))^(d/2)*exp(-r*r/(l_i*l_i + l_j*l_j)),l_i)))

%                                         2
%       /    l_i l_j  \(1/2 d)           r               4        4        2  2
%   1/2 |2 -----------|        exp(- -----------) (-d l_i  + d l_j  + 4 l_i  r )
%       |     2      2|                 2      2
%       \  l_i  + l_j /              l_i  + l_j

%            /      2      2 2
%           /  ((l_i  + l_j )  l_i)
%          /
% >> pretty(simple(diff((2*(l_i*l_i)/(l_i*l_i + l_i*l_i))^(d/2)*exp(-r*r/(l_i*l_i + l_i*l_i)),l_i)))

%                                              2
%                                2            r
%                               r  exp(- 1/2 ----)
%                                               2
%                                            l_i
%                               ------------------
%                                         3
%                                      l_i
% >> pretty(simple(diff((2*(l_i*l_i)/(l_i*l_i + l_i*l_i))^(d/2)*exp(-r*r/(l_i*l_i + l_i*l_i)),l_i)))

fhandle = str2func([kern.lengthScaleTransform, 'Transform']);
g = zeros(1, kern.nParams);
% The last argument is covGrad
if length(varargin)<2
  [k, sk, n2, w2, l] = gibbsKernCompute(kern, x);
  gOut = modelOutputGrad(kern.lengthScaleFunc, x);
  gradFact = fhandle(l, 'gradfact');
  L1 = repmat(l, 1, size(l, 1));
  L2 = L1';
  covGrad = varargin{end};
  covGrad(1:size(covGrad, 1)+1:end) = 0;
  base = covGrad.*k;
  base2 = base.*(kern.inputDimension/2*(L2.^4 - L1.^4) + 2*L1.*L1.*n2)./(w2.*w2.*L1);
  for i = 1:size(g, 2)-1
    g(i) = g(i) + 2*sum(sum(base2.*repmat(gOut(:, i).*gradFact, 1, size(x, 1))));
  end
%/~
%  covGrad = varargin{end};
%  covGradDiag = diag(covGrad);
%  covGrad(1:size(covGrad, 1)+1:end) = 0;
%  base = covGradDiag.*diag(k)./(l.^3).*diag(n2);
%  for i = 1:size(g, 2) - 1
%    g(i) = g(i) + sum(base.*gOut(:, i));
%  end
%  base = covGrad.*k;
%  base = base.*(d/2*(L2.^4 - L1.^4) + 2*L1.*L1.*n2)./(w2.*w2.*L1);
%  for i = 1:size(g, 2)-1
%    g(i) = g(i) + sum(sum(base.*repmat(gOut(:, i), 1, size(x, 1))));
%  end
%~/
else
  [k, sk, n2, w2, l, l2] = gibbsKernCompute(kern, x, varargin{1});
  gOut = modelOutputGrad(kern.lengthScaleFunc, x);
  gradFact = fhandle(l, 'gradfact');
  gOut2 = modelOutputGrad(kern.lengthScaleFunc, varargin{1});
  gradFact2 = fhandle(l2, 'gradfact');
  L1 = repmat(l, 1, size(l2, 1));
  L2 = repmat(l2, 1, size(l, 1))';
  base = varargin{end}.*k;
  base2 = base.*(kern.inputDimension/2*(L2.^4 - L1.^4) + 2*L1.*L1.*n2)./(w2.*w2.*L1);
  for i = 1:size(g, 2)-1
    g(i) = g(i) + sum(sum(base2.*repmat(gOut(:, i).*gradFact, 1, size(varargin{1}, ...
                                                      1))));
  end
  base2 = base.*(kern.inputDimension/2*(L1.^4 - L2.^4) + 2*L2.*L2.* ...
                 n2)./(w2.*w2.*L2);
  for i = 1:size(g, 2)-1
    g(i) = g(i) + sum(sum(base2.*repmat(gOut2(:, i)'.*gradFact2', size(x, ...
                                                      1), 1)));
  end

end
g(end) =  sum(sum(varargin{end}.*sk));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function b = r83t_mv ( m, n, a, x )

%*****************************************************************************80
%
%% R83T_MV multiplies an R83T matrix times an R8VEC.
%
%  Discussion:
%
%    The R83T storage format is used for a tridiagonal matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,3), the matrix.
%
%    Input, real X(N), the vector to be multiplied by A.
%
%    Output, real B(M), the product A * x.
%
  b = zeros ( m, 1 );
  x = x(:);

  mn = min ( m, n );

  if ( n == 1 )
    b(1) = a(1,2) * x(1);
    if ( 1 < m )
      b(2) = a(2,1) * x(1);
    end
    return
  end

  b(1)      = a(1,2)       * x(1) ...
            + a(1,3)       * x(2);

  b(2:mn-1) = a(2:mn-1,1) .* x(1:mn-2) ...
            + a(2:mn-1,2) .* x(2:mn-1) ...
            + a(2:mn-1,3) .* x(3:mn);

  b(mn)     = a(mn,1)      * x(mn-1) ...
            + a(mn,2)      * x(mn);

  if ( n < m )
    b(mn+1) = b(mn+1) + a(mn+1,1) * x(mn);
  elseif ( m < n )
    b(mn) = b(mn) + a(mn,3) * x(mn+1);
  end

  return
end


text:
% Chebychev design of an FIR filter given a desired H(w)
% "Filter design" lecture notes (EE364) by S. Boyd
% (figures are generated)
%
% Designs an FIR filter given a desired frequency response H_des(w).
% The design is judged by the maximum absolute error (Chebychev norm).
% This is a convex problem (after sampling it can be formulated as an SOCP).
%
%   minimize   max |H(w) - H_des(w)|     for w in [0,pi]
%
% where H is the frequency response function and variable is h
% (the filter impulse response).
%
% Written for CVX by Almir Mutapcic 02/02/06

%********************************************************************
% problem specs
%********************************************************************
% number of FIR coefficients (including the zeroth one)
n = 20;

% rule-of-thumb frequency discretization (Cheney's Approx. Theory book)
m = 15*n;
w = linspace(0,pi,m)'; % omega

%********************************************************************
% construct the desired filter
%********************************************************************
% fractional delay
D = 8.25;            % delay value
Hdes = exp(-1j*D*w);  % desired frequency response

% Gaussian filter with linear phase (uncomment lines below for this design)
% var = 0.05;
% Hdes = 1/(sqrt(2*pi*var))*exp(-(w-pi/2).^2/(2*var));
% Hdes = Hdes.*exp(-j*n/2*w);

%*********************************************************************
% solve the minimax (Chebychev) design problem
%*********************************************************************
% A is the matrix used to compute the frequency response
% A(w,:) = [1 exp(-j*w) exp(-j*2*w) ... exp(-j*n*w)]
A = exp( -1j*kron(w,0:n-1) );

% optimal Chebyshev filter formulation
cvx_begin
  variable h(n,1)
  minimize( max( abs( A*h - Hdes ) ) )
cvx_end

% check if problem was successfully solved
disp(['Problem is ' cvx_status])
if ~strfind(cvx_status,'Solved')
  h = [];
end

%*********************************************************************
% plotting routines
%*********************************************************************
% plot the FIR impulse reponse
figure(1)
stem(0:n-1,h)
xlabel('n')
ylabel('h(n)')

% plot the frequency response
H = exp(-1j*kron(w,0:n-1))*h;
figure(2)
% magnitude
subplot(2,1,1);
plot(w,20*log10(abs(H)),w,20*log10(abs(Hdes)),'--')
xlabel('w')
ylabel('mag H in dB')
axis([0 pi -30 10])
legend('optimized','desired','Location','SouthEast')
% phase
subplot(2,1,2)
plot(w,angle(H))
axis([0,pi,-pi,pi])
xlabel('w'), ylabel('phase H(w)')

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [ a, b ] = p38_lim ( )

%*****************************************************************************80
%
%% P38_LIM returns the integration limits for problem 51.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = pi;

  return
end


text:
function node_order = tet_mesh_node_order ( tetra_order, tetra_num, ...
  tetra_node, node_num )

%*****************************************************************************80
%
%% TET_MESH_NODE_ORDER: determines the order of nodes.
%
%  Discussion:
%
%    The order of a node is the number of tetrahedrons that use that node
%    as a vertex.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer TETRA_ORDER, the order of the tetrahedrons.
%
%    Input, integer TETRA_NUM, the number of tetrahedrons.
%
%    Input, integer TETRA_NODE(TETRA_ORDER,TETRA_NUM), the nodes
%    that make up the tetrahedrons.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Output, integer NODE_ORDER(NODE_NUM), the order of each node.
%
  node_order(1:node_num) = 0;

  for tetra = 1 : tetra_num
    for i = 1 : tetra_order
      node = tetra_node(i,tetra);
      if ( node < 1 | node_num < node )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'TET_MESH_NODE_ORDER - Fatal error!\n' );
        fprintf ( 1, '  Illegal entry in TETRA_NODE.\n' );
        error ( 'TET_MESH_NODE_ORDER - Fatal error!' );
      else
        node_order(node) = node_order(node) + 1;
      end
    end
  end

  return
end


text:
function [varargout] = ft_qualitycheck(cfg)

% FT_QUALITYCHECK performs a quality inspection of a given MEG/EEG dataset,
% stores (.mat), and visualizes the result (.png and .pdf).
%
% This function segments the data into 10-second pieces and performs the
% following analyses:
%  1) reads the properties of the dataset
%  2) computes the headpositions and distance covered from recording onset (CTF only)
%  3) computes the mean, max, min, and range of the signal amplitude
%  4) detects trigger events
%  5) detects jump artifacts
%  6) computes the powerspectrum
%  7) estimates the low-frequency (<2 Hz) and line noise (~50 Hz)
%
% Use as
%   [info, timelock, freq, summary, headpos] = ft_qualitycheck(cfg)
% where info contains the dataset properties, timelock the timelocked data,
% freq the powerspectra, summary the mean descriptives, and headpos the
% headpositions throughout the recording
%
% The configuration should contain:
%   cfg.dataset = a string (e.g. 'dataset.ds')
%
% The following parameters can be used:
%   cfg.analyze   = string, 'yes' or 'no' to analyze the dataset (default = 'yes')
%   cfg.savemat   = string, 'yes' or 'no' to save the analysis (default = 'yes')
%   cfg.matfile   = string, filename (e.g. 'previousoutput.mat'), preferably in combination
%                    with analyze = 'no'
%   cfg.visualize = string, 'yes' or 'no' to visualize the analysis (default = 'yes')
%   cfg.saveplot  = string, 'yes' or 'no' to save the visualization (default = 'yes')
%   cfg.linefreq  = scalar, frequency of power line (default = 50)
%   cfg.plotunit  = scalar, the length of time to be plotted in one panel (default = 3600)
%
% See also FT_PREPROCESSING, FT_READ_HEADER, FT_READ_DATA, FT_READ_EVENT

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

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

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

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

% set the defaults
cfg.analyze   = ft_getopt(cfg, 'analyze',   'yes');
cfg.savemat   = ft_getopt(cfg, 'savemat',   'yes');
cfg.matfile   = ft_getopt(cfg, 'matfile',   []);
cfg.visualize = ft_getopt(cfg, 'visualize', 'yes');
cfg.saveplot  = ft_getopt(cfg, 'saveplot',  'yes');
cfg.linefreq  = ft_getopt(cfg, 'linefreq',  50);
cfg.plotunit  = ft_getopt(cfg, 'plotunit',  3600);
cfg.trl       = ft_getopt(cfg, 'trl', []);
%% ANALYSIS
if strcmp(cfg.analyze,'yes')
  tic

  % checks
  cfg   = ft_checkconfig(cfg, 'dataset2files', 'yes'); % translate into datafile+headerfile

  % these will be replaced by more appropriate values
  info.filename    = cfg.dataset;
  info.datasetname = 'unknown';
  info.starttime   = 'unknown';
  info.startdate   = 'unknown';
  info.stoptime    = 'unknown';
  info.stopdate    = 'unknown';

  % the exportname is also used in the cron job
  exportname = qualitycheck_exportname(cfg.dataset);
  [iseeg, ismeg, isctf, fltp] = filetyper(cfg.dataset);
  if isctf
    try
      % update the info fields
      info = read_ctf_hist(cfg.dataset);
    end
  end

  % add info
  info.event                  = ft_read_event(cfg.dataset);
  info.hdr                    = ft_read_header(cfg.dataset);
  info.filetype               = fltp;

  if isempty(cfg.trl)
    % trial definition
    cfgdef                      = [];
    cfgdef.dataset              = cfg.dataset;
    cfgdef.trialfun             = 'ft_trialfun_general';
    cfgdef.trialdef.triallength = 10;
    cfgdef.continuous           = 'yes';
    cfgdef                      = ft_definetrial(cfgdef);
    cfgdef.trl = cfgdef.trl(1:end-1, :); % remove last trial
  else
    % UNDOCUMENTED POSSIBILITY to add cfg.trl
    cfgdef.trl = cfg.trl;
  end
  ntrials                     = size(cfgdef.trl,1);
  timeunit                    = cfgdef.trialdef.triallength;

  % channelselection for jump detection (all) and for FFT (brain)
  if ismeg
    allchans                   = ft_channelselection({'MEG','MEGREF'}, info.hdr.label);
    chans                      = ft_channelselection('MEG', info.hdr.label); % brain
    allchanindx                = match_str(info.hdr.label, allchans);
    chanindx                   = match_str(chans, allchans);
    jumpthreshold              = 1e-10;
  elseif iseeg
    allchans                   = ft_channelselection('EEG', info.hdr.label);
    if isempty(allchans)
      % some EEG systems and data files use non-standard channel names that are not detected automatically
      ft_warning('no EEG channels detected, selecting all channels');
      allchans = info.hdr.label;
    end
    chans                      = allchans;  % brain
    allchanindx                = match_str(info.hdr.label, allchans);
    chanindx                   = match_str(chans, allchans);
    jumpthreshold              = 1e4;
  end

  % find headcoil channels
  if isctf % this fails for older CTF data sets
    Nx = strmatch('HLC0011', info.hdr.label); % x nasion coil
    Ny = strmatch('HLC0012', info.hdr.label); % y nasion
    Nz = strmatch('HLC0013', info.hdr.label); % z nasion
    Lx = strmatch('HLC0021', info.hdr.label); % x left coil
    Ly = strmatch('HLC0022', info.hdr.label); % y left
    Lz = strmatch('HLC0023', info.hdr.label); % z left
    Rx = strmatch('HLC0031', info.hdr.label); % x right coil
    Ry = strmatch('HLC0032', info.hdr.label); % y right
    Rz = strmatch('HLC0033', info.hdr.label); % z right
    headpos.dimord = 'chan_time';
    headpos.time   = (timeunit-timeunit/2:timeunit:timeunit*ntrials-timeunit/2);
    headpos.label  = {'Nx';'Ny';'Nz';'Lx';'Ly';'Lz';'Rx';'Ry';'Rz'};
    headpos.avg    = NaN(length(headpos.label), ntrials);
    headpos.grad   = info.hdr.grad;

    if numel(cat(1,Nx,Ny,Nz,Lx,Ly,Lz,Rx,Ry,Rz))==9
      hasheadpos = true;
    else
      hasheadpos = false;
    end

  end % if

  % analysis settings
  cfgredef             = [];
  cfgredef.length      = 1;
  cfgredef.overlap     = 0;

  cfgfreq              = [];
  cfgfreq.output       = 'pow';
  cfgfreq.channel      = allchans;
  cfgfreq.method       = 'mtmfft';
  cfgfreq.taper        = 'hanning';
  cfgfreq.keeptrials   = 'no';
  cfgfreq.foilim       = [0 min(info.hdr.Fs/2, 400)];

  % output variables
  timelock.dimord = 'chan_time';
  timelock.label  = allchans;
  timelock.time   = (timeunit-timeunit/2:timeunit:timeunit*ntrials-timeunit/2);
  timelock.avg    = NaN(length(allchans), ntrials); % updated in loop
  timelock.median = NaN(length(allchans), ntrials); % updated in loop
  timelock.jumps  = NaN(length(allchans), ntrials); % updated in loop
  timelock.range  = NaN(length(allchans), ntrials); % updated in loop
  timelock.min    = NaN(length(allchans), ntrials); % updated in loop
  timelock.max    = NaN(length(allchans), ntrials); % updated in loop

  freq.dimord     = 'chan_freq_time';
  freq.label      = allchans;
  freq.freq       = (cfgfreq.foilim(1):cfgfreq.foilim(2));
  freq.time       = (timeunit-timeunit/2:timeunit:timeunit*ntrials-timeunit/2);
  freq.powspctrm  = NaN(length(allchans), length(freq.freq), ntrials); % updated in loop

  summary.dimord  = 'chan_time';
  summary.time    = (timeunit-timeunit/2:timeunit:timeunit*ntrials-timeunit/2);
  summary.label   = {'Mean';'Median';'Min';'Max';'Range';'HmotionN';'HmotionL';'HmotionR';'LowFreqPower';'LineFreqPower';'Jumps'};
  summary.avg     = NaN(length(summary.label), ntrials); % updated in loop

  % try add gradiometer info
  if isfield(info.hdr, 'grad')
    timelock.grad = info.hdr.grad;
    freq.grad     = info.hdr.grad;
    summary.grad  = info.hdr.grad;
  end

  % process trial by trial
  for t = 1:ntrials
    fprintf('analyzing trial %s of %s \n', num2str(t), num2str(ntrials));

    % preprocess
    cfgpreproc     = cfgdef;
    cfgpreproc.trl = cfgdef.trl(t,:);
    data           = ft_preprocessing(cfgpreproc); clear cfgpreproc;

    % determine headposition
    if isctf && hasheadpos
      headpos.avg(1,t) = mean(data.trial{1,1}(Nx,:) * 100);  % meter to cm
      headpos.avg(2,t) = mean(data.trial{1,1}(Ny,:) * 100);
      headpos.avg(3,t) = mean(data.trial{1,1}(Nz,:) * 100);
      headpos.avg(4,t) = mean(data.trial{1,1}(Lx,:) * 100);
      headpos.avg(5,t) = mean(data.trial{1,1}(Ly,:) * 100);
      headpos.avg(6,t) = mean(data.trial{1,1}(Lz,:) * 100);
      headpos.avg(7,t) = mean(data.trial{1,1}(Rx,:) * 100);
      headpos.avg(8,t) = mean(data.trial{1,1}(Ry,:) * 100);
      headpos.avg(9,t) = mean(data.trial{1,1}(Rz,:) * 100);
    end

    % update values
    timelock.avg(:,t)    = mean(data.trial{1}(allchanindx,:),2);
    timelock.median(:,t) = median(data.trial{1}(allchanindx,:),2);
    timelock.range(:,t)  = max(data.trial{1}(allchanindx,:),[],2) - min(data.trial{1}(allchanindx,:),[],2);
    timelock.min(:,t)    = min(data.trial{1}(allchanindx,:),[],2);
    timelock.max(:,t)    = max(data.trial{1}(allchanindx,:),[],2);

    % detect jumps
    for c = 1:size(data.trial{1}(allchanindx,:),1)
      timelock.jumps(c,t) = length(find(diff(data.trial{1,1}(allchanindx(c),:)) > jumpthreshold));
    end

    % remove the grad and elec to make the subsequent code faster. See https://github.com/fieldtrip/fieldtrip/issues/1551
    data = removefields(data, {'grad', 'elec'});

    % FFT and noise estimation
    redef                 = ft_redefinetrial(cfgredef, data); clear data;
    FFT                   = ft_freqanalysis(cfgfreq, redef); clear redef;
    freq.powspctrm(:,:,t) = FFT.powspctrm;
    summary.avg(9,t)      = mean(mean(findpower(0,  2,  FFT, chanindx))); % Low Freq Power
    summary.avg(10,t)     = mean(mean(findpower(cfg.linefreq-1, cfg.linefreq+1, FFT, chanindx))); clear FFT; % Line Freq Power

    toc
  end % end of trial loop

  % determine headmotion: distance from initial trial (in cm)
  if isctf && hasheadpos
    summary.avg(6,:) = sqrt(sum((headpos.avg(1:3,:)-repmat(headpos.avg(1:3,1),1,size(headpos.avg,2))).^2,1)); % N
    summary.avg(7,:) = sqrt(sum((headpos.avg(4:6,:)-repmat(headpos.avg(4:6,1),1,size(headpos.avg,2))).^2,1)); % L
    summary.avg(8,:) = sqrt(sum((headpos.avg(7:9,:)-repmat(headpos.avg(7:9,1),1,size(headpos.avg,2))).^2,1)); % R
  end

  % summarize/mean and store variables of brain info only
  summary.avg(1,:)   = mean(timelock.avg(chanindx,:),1);
  summary.avg(2,:)   = mean(timelock.median(chanindx,:),1);
  summary.avg(3,:)   = mean(timelock.min(chanindx,:),1);
  summary.avg(4,:)   = mean(timelock.max(chanindx,:),1);
  summary.avg(5,:)   = mean(timelock.range(chanindx,:),1);
  summary.avg(11,:)  = mean(timelock.jumps(chanindx,:),1);

  % save to .mat
  if strcmp(cfg.savemat, 'yes')
    if isctf && hasheadpos
      headpos.cfg        = cfg;
      save(exportname, 'info','timelock','freq','summary','headpos');
    else
      save(exportname, 'info','timelock','freq','summary');
    end
  end

end % end of analysis

%% VISUALIZATION
if strcmp(cfg.visualize, 'yes')

  % load data
  if strcmp(cfg.analyze, 'no')
    if ~isempty(cfg.matfile)
      exportname = cfg.matfile;
    else
      exportname = qualitycheck_exportname(cfg.dataset);
    end
    fprintf('loading %s \n', exportname);
    load(exportname);
  end

  % determine number of 1-hour plots to be made
  nplots = ceil(length(freq.time)/(cfg.plotunit/10));

  % create GUI-like figure(s)
  for p = 1:nplots
    fprintf('visualizing %s of %s \n', num2str(p), num2str(nplots));
    toi = [p*cfg.plotunit-(cfg.plotunit-5) p*cfg.plotunit-5]; % select 1-hour chunks

    tmpcfg.latency = toi;
    temp_timelock  = ft_selectdata(tmpcfg, timelock);
    temp_freq      = ft_selectdata(tmpcfg, freq);
    temp_summary   = ft_selectdata(tmpcfg, summary);
    if exist('headpos','var')
      temp_headpos  = ft_selectdata(tmpcfg, headpos);
      draw_figure(info, temp_timelock, temp_freq, temp_summary, temp_headpos, toi);
      clear temp_timelock; clear temp_freq; clear temp_summary; clear temp_headpos; clear toi;
    else
      draw_figure(info, temp_timelock, temp_freq, temp_summary, toi);
      clear temp_timelock; clear temp_freq; clear temp_summary; clear toi;
    end

    % export to .PNG and .PDF
    if strcmp(cfg.saveplot, 'yes')
      [pathstr,name,extr] = fileparts(exportname);
      if p == 1
        exportfilename = name;
      else
        exportfilename = strcat(name,'_pt',num2str(p));
      end
      fprintf('exporting %s of %s \n', num2str(p), num2str(nplots));
      set(gcf, 'PaperType', 'a4');
      print(gcf, '-dpng', strcat(exportfilename,'.png'));
      orient landscape;
      print(gcf, '-dpdf', strcat(exportfilename,'.pdf'));
      close
    end
  end % end of nplots
end % end of visualization

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug

ft_postamble provenance
ft_postamble history timelock   % add the input cfg to multiple outputs
ft_postamble history freq       % add the input cfg to multiple outputs
ft_postamble history summary    % add the input cfg to multiple outputs

%% VARARGOUT
if nargout>0
  mOutputArgs{1} = info;
  mOutputArgs{2} = timelock;
  mOutputArgs{3} = freq;
  mOutputArgs{4} = summary;
  try
    mOutputArgs{5} = headpos;
  end
  [varargout{1:nargout}] = mOutputArgs{:};
  clearvars -except varargout
else
  clear
end

%%%%%%%%%%%%%%%%%%%%%%%% SUBFUNCTION %%%%%%%%%%%%%%%%%%%%%%%%%
function [x] = clipat(x, v, v2)
v = [v v2]; % clip between value v and v2
if length(v) == 1
  x(x>v) = v;
elseif length(v) == 2
  x(x<v(1)) = v(1);
  x(x>v(2)) = v(2);
end

%%%%%%%%%%%%%%%%%%%%%%%% SUBFUNCTION %%%%%%%%%%%%%%%%%%%%%%%%%
function [iseeg, ismeg, isctf, fltp] = filetyper(dataset)
fltp = ft_filetype(dataset);
iseeg = ft_filetype(dataset,'brainvision_eeg') | ...
  ft_filetype(dataset,'ns_eeg') | ...
  ft_filetype(dataset,'bci2000_dat') | ...
  ft_filetype(dataset,'neuroprax_eeg') | ...
  ft_filetype(dataset,'egi_sbin') | ...
  ft_filetype(dataset,'biosemi_bdf');
ismeg = ft_filetype(dataset,'ctf_ds') | ...
  ft_filetype(dataset,'4d') | ...
  ft_filetype(dataset,'neuromag_fif') | ...
  ft_filetype(dataset,'itab_raw');
isctf = ft_filetype(dataset, 'ctf_ds');
if ~ismeg && ~iseeg % if none found, try less strict checks
  [p, f, ext] = fileparts(dataset);
  if strcmp(ext, '.eeg')
    fltp = 'brainvision_eeg';
    iseeg = 1;
  elseif strcmp(ext, '.bdf')
    fltp = 'biosemi_bdf';
    iseeg = 1;
  elseif strcmp(ext, '.ds')
    fltp = 'ctf_ds';
    ismeg = 1;
  else % otherwise use eeg settings for stability reasons
    iseeg = 1;
  end
end

%%%%%%%%%%%%%%%%%%%%%%%% SUBFUNCTION %%%%%%%%%%%%%%%%%%%%%%%%%
function [power, freq] = findpower(low, high, freqinput, chans)
% replace value with the index of the nearest bin
xmin  = nearest(getsubfield(freqinput, 'freq'), low);
xmax  = nearest(getsubfield(freqinput, 'freq'), high);
% select the freq range
power = freqinput.powspctrm(chans, xmin:xmax);
freq  = freqinput.freq(:, xmin:xmax);

%%%%%%%%%%%%%%%%%%%%%%%% SUBFUNCTION %%%%%%%%%%%%%%%%%%%%%%%%%
function draw_figure(varargin)
% deal with input
if nargin == 6
  info     = varargin{1};
  timelock = varargin{2};
  freq     = varargin{3};
  summary  = varargin{4};
  headpos  = varargin{5};
  toi      = varargin{6};
elseif nargin == 5
  info     = varargin{1};
  timelock = varargin{2};
  freq     = varargin{3};
  summary  = varargin{4};
  toi      = varargin{5};
end

% determine whether it is EEG or MEG
if isfield(info, 'filename') % supported as of January 2018
  filename = info.filename;
elseif isfield(timelock.cfg, 'dataset')
  filename = timelock.cfg.dataset;
elseif isfield(info.hdr.orig, 'FileName')
  filename = info.hdr.orig.FileName;
elseif exist('headpos','var') && isfield(headpos.cfg.previous, 'dataset')
  filename = headpos.cfg.previous.dataset;
else
  error('could not determine the filename');
end
[iseeg, ismeg, isctf, fltp] = filetyper(filename);

if ismeg
  scaling = 1e15; % assuming data is in T and needs to become fT
  powscaling = scaling^2;
  ylab = 'fT';
elseif iseeg
  scaling = 1e0; % assuming data is in muV already
  powscaling = scaling^2;
  ylab = '\muV';
end

% PARENT FIGURE
h.MainFigure = figure(...
  'MenuBar','none',...
  'Name','ft_qualitycheck',...
  'Units','normalized',...
  'color','white',...
  'Position',[0.01 0.01 .99 .99]); % nearly fullscreen

if strcmp(info.startdate,'unknown')
  tmp = 'unknown';
else
  [d,w] = weekday(info.startdate);
  tmp = [w ' ' info.startdate];
end

h.MainText = uicontrol(...
  'Parent',h.MainFigure,...
  'Style','text',...
  'Units','normalized',...
  'FontSize',10,...
  'String',tmp,...
  'Backgroundcolor','white',...
  'Position',[.06 .96 .15 .02]);

h.MainText2 = uicontrol(...
  'Parent',h.MainFigure,...
  'Style','text',...
  'Units','normalized',...
  'FontSize',10,...
  'String','Jump artefacts',...
  'Backgroundcolor','white',...
  'Position',[.08 .46 .12 .02]);

h.MainText3 = uicontrol(...
  'Parent',h.MainFigure,...
  'Style','text',...
  'Units','normalized',...
  'FontSize',10,...
  'String','Mean powerspectrum',...
  'Backgroundcolor','white',...
  'Position',[.4 .3 .15 .02]);

h.MainText4 = uicontrol(...
  'Parent',h.MainFigure,...
  'Style','text',...
  'Units','normalized',...
  'FontSize',10,...
  'String','Timecourses',...
  'Backgroundcolor','white',...
  'Position',[.5 .96 .11 .02]);

h.MainText5 = uicontrol(...
  'Parent',h.MainFigure,...
  'Style','text',...
  'Units','normalized',...
  'FontSize',10,...
  'String','Events',...
  'Backgroundcolor','white',...
  'Position',[.81 .3 .06 .02]);

% HEADMOTION PANEL
h.HmotionPanel = uipanel(...
  'Parent',h.MainFigure,...
  'Units','normalized',...
  'Backgroundcolor','white',...
  'Position',[.01 .5 .25 .47]);

h.DataText = uicontrol(...
  'Parent',h.HmotionPanel,...
  'Style','text',...
  'Units','normalized',...
  'FontSize',10,...
  'String',info.datasetname,...
  'Backgroundcolor','white',...
  'Position',[.01 .85 .99 .1]);

h.TimeText = uicontrol(...
  'Parent',h.HmotionPanel,...
  'Style','text',...
  'Units','normalized',...
  'FontSize',10,...
  'String',[info.starttime ' - ' info.stoptime],...
  'Backgroundcolor','white',...
  'Position',[.01 .78 .99 .1]);

if ismeg
  allchans = ft_senslabel(ft_senstype(timelock));
  misschans = setdiff(ft_channelselection('MEG', info.hdr.label), allchans);
  nchans = numel(ft_channelselection('MEG', info.hdr.label));
else
  misschans = '';
  nchans = numel(ft_channelselection('EEG', info.hdr.label));
end

h.DataText2 = uicontrol(...
  'Parent',h.HmotionPanel,...
  'Style','text',...
  'Units','normalized',...
  'FontSize',10,...
  'String',[fltp ', fs: ' num2str(info.hdr.Fs) ', nchans: ' num2str(nchans)],...
  'Backgroundcolor','white',...
  'Position',[.01 .71 .99 .1]);

h.DataText3 = uicontrol(...
  'Parent',h.HmotionPanel,...
  'Style','text',...
  'Units','normalized',...
  'FontSize',10,...
  'String',['missing chans: ' misschans'],...
  'Backgroundcolor','white',...
  'Position',[.01 .64 .99 .1]);

% boxplot headmotion (*10; cm-> mm) per coil
if exist('headpos','var')
  h.HmotionAxes = axes(...
    'Parent',h.HmotionPanel,...
    'Units','normalized',...
    'color','white',...
    'Position',[.05 .08 .9 .52]);

  hmotions = ([summary.avg(8,:)' summary.avg(7,:)' summary.avg(6,:)'])*10;
  boxplot(h.HmotionAxes, hmotions, 'orientation', 'horizontal', 'notch', 'on');
  set(h.HmotionAxes,'YTick',1:3);
  set(h.HmotionAxes,'YTickLabel',{'R','L','N'});
  xlim(h.HmotionAxes, [0 10]);
  xlabel(h.HmotionAxes, 'Headmotion from start [mm]');
end

% TIMECOURSE PANEL
h.SignalPanel = uipanel(...
  'Parent',h.MainFigure,...
  'Units','normalized',...
  'Backgroundcolor','white',...
  'Position',[.28 .34 .71 .63]);

h.SignalAxes = axes(...
  'Parent',h.SignalPanel,...
  'Units','normalized',...
  'color','white',...
  'Position',[.08 .36 .89 .3]);

h.LinenoiseAxes = axes(...
  'Parent',h.SignalPanel,...
  'Units','normalized',...
  'color','white',...
  'Position',[.08 .23 .89 .1]);

h.LowfreqnoiseAxes = axes(...
  'Parent',h.SignalPanel,...
  'Units','normalized',...
  'color','white',...
  'Position',[.08 .1 .89 .1]);

% plot hmotion timecourses per coil (*10; cm-> mm)
if exist('headpos','var')
  h.HmotionTimecourseAxes = axes(...
    'Parent',h.SignalPanel,...
    'Units','normalized',...
    'color','white',...
    'Position',[.08 .73 .89 .22]);

  plot(h.HmotionTimecourseAxes, summary.time, clipat(summary.avg(6,:)*10, 0, 10), ...
    summary.time, clipat(summary.avg(7,:)*10, 0, 10), ...
    summary.time, clipat(summary.avg(8,:)*10, 0, 10), 'LineWidth',2);
  ylim(h.HmotionTimecourseAxes,[0 10]);
  ylabel(h.HmotionTimecourseAxes, 'Coil distance [mm]');
  xlim(h.HmotionTimecourseAxes,toi);
  grid(h.HmotionTimecourseAxes,'on');
  legend(h.HmotionTimecourseAxes, 'N','L','R');
end

% plot mean and range of the raw signal
plot(h.SignalAxes, summary.time, summary.avg(5,:)*scaling, summary.time, summary.avg(1,:)*scaling, 'LineWidth', 2);
set(h.SignalAxes,'Nextplot','add');
plot(h.SignalAxes, summary.time, summary.avg(3,:)*scaling, summary.time, summary.avg(4,:)*scaling, 'LineWidth', 1, 'Color', [255/255 127/255 39/255]);
grid(h.SignalAxes,'on');
ylabel(h.SignalAxes, ['Amplitude [' ylab ']']);
xlim(h.SignalAxes,toi);
legend(h.SignalAxes,'Range','Mean','Min','Max');
set(h.SignalAxes,'XTickLabel','');

% plot linenoise
semilogy(h.LinenoiseAxes, summary.time, clipat(summary.avg(10,:)*powscaling, 1e2, 1e4), 'LineWidth',2);
grid(h.LinenoiseAxes,'on');
legend(h.LinenoiseAxes, ['LineFreq [' ylab '^2/Hz]']);
set(h.LinenoiseAxes,'XTickLabel','');
xlim(h.LinenoiseAxes,toi);
ylim(h.LinenoiseAxes,[1e2 1e4]); % before april 28th this was 1e0 - 1e3

% plot lowfreqnoise
semilogy(h.LowfreqnoiseAxes, summary.time, clipat(summary.avg(9,:)*powscaling, 1e10, 1e12), 'LineWidth',2);
grid(h.LowfreqnoiseAxes,'on');
xlim(h.LowfreqnoiseAxes,toi);
ylim(h.LowfreqnoiseAxes,[1e10 1e12]);
legend(h.LowfreqnoiseAxes, ['LowFreq [' ylab '^2/Hz]']);
xlabel(h.LowfreqnoiseAxes, 'Time [seconds]'); % before april 28th this was 1e0 - 1e10

% EVENT PANEL
h.EventPanel = uipanel(...
  'Parent',h.MainFigure,...
  'Units','normalized',...
  'Backgroundcolor','white',...
  'Position',[.7 .01 .29 .3]);

% event details
eventtypes    = {};
eventtriggers = {};
eventvalues   = {};

if ~isempty(info.event)
  [a,b,c] = unique({info.event.type});
  for j=1:length(a)
    eventtypes{j,1}    = a{j};
    eventtriggers{j,1} = sum(c==j);
    eventvalues{j,1}   = length(unique([info.event(c==j).value]));
  end
end
if isempty(eventtypes)
  eventtypes{1,1} = 'no triggers found';
end

h.EventText = uicontrol(...
  'Parent',h.EventPanel,...
  'Style','text',...
  'Units','normalized',...
  'FontSize',10,...
  'String',['Types'; ' '; eventtypes],...
  'Backgroundcolor','white',...
  'Position',[.05 .05 .4 .85]);

h.EventText2 = uicontrol(...
  'Parent',h.EventPanel,...
  'Style','text',...
  'Units','normalized',...
  'FontSize',10,...
  'String',['Triggers'; ' '; eventtriggers],...
  'Backgroundcolor','white',...
  'Position',[.55 .05 .2 .85]);

h.EventText3 = uicontrol(...
  'Parent',h.EventPanel,...
  'Style','text',...
  'Units','normalized',...
  'FontSize',10,...
  'String',['Values'; ' '; eventvalues],...
  'Backgroundcolor','white',...
  'Position',[.8 .05 .15 .85]);

% POWERSPECTRUM PANEL
h.SpectrumPanel = uipanel(...
  'Parent',h.MainFigure,...
  'Units','normalized',...
  'Backgroundcolor','white',...
  'Position',[.28 .01 .4 .3]);

h.SpectrumAxes = axes(...
  'Parent',h.SpectrumPanel,...
  'Units','normalized',...
  'color','white',...
  'Position',[.15 .2 .8 .7]);

% plot powerspectrum
loglog(h.SpectrumAxes, freq.freq, mean(mean(freq.powspctrm,1),3)*powscaling,'r','LineWidth',2);
xlabel(h.SpectrumAxes, 'Frequency [Hz]');
ylabel(h.SpectrumAxes, ['Power [' ylab '^2/Hz]']);

% ARTEFACT PANEL
h.JumpPanel = uipanel(...
  'Parent',h.MainFigure,...
  'Units','normalized',...
  'Backgroundcolor','white',...
  'Position',[.01 .01 .25 .46]);

% jump details
jumpchans  = {};
jumpcounts = {};
[jumps,i] = find(timelock.jumps>0); % find all jumps
[a,b,c] = unique(jumps);
for j=1:length(a)
  jumpchans{j,1}  = timelock.label{a(j)};
  jumpcounts{j,1} = sum(c==j);
end
if isempty(jumpchans)
  jumpchans{1,1} = 'no jumps detected';
end

h.JumpText = uicontrol(...
  'Parent',h.JumpPanel,...
  'Style','text',...
  'Units','normalized',...
  'FontSize',10,...
  'String',jumpchans,...
  'Backgroundcolor','white',...
  'Position',[.15 .5 .25 .4]);

h.JumpText2 = uicontrol(...
  'Parent',h.JumpPanel,...
  'Style','text',...
  'Units','normalized',...
  'FontSize',10,...
  'String',jumpcounts,...
  'Backgroundcolor','white',...
  'Position',[.65 .5 .2 .4]);

% plot jumps on the dewar sensors
if ismeg
  h.TopoMEG = axes(...
    'Parent',h.JumpPanel,...
    'color','white',...
    'Units','normalized',...
    'Position',[0.4 0.05 0.55 0.4]);

  MEGchans                 = ft_channelselection('MEG', timelock.label);
  MEGchanindx              = match_str(timelock.label, MEGchans);
  cfgtopo                  = [];
  cfgtopo.marker           = 'off';
  cfgtopo.colorbar         = 'no';
  cfgtopo.comment          = 'no';
  cfgtopo.style            = 'blank';
  cfgtopo.layout           = ft_prepare_layout([], timelock);
  cfgtopo.highlight        = 'on';
  cfgtopo.highlightsymbol  = '.';
  cfgtopo.highlightsize    = 14;
  cfgtopo.highlightchannel = find(sum(timelock.jumps(MEGchanindx,:),2)>0);
  data.label               = MEGchans;
  data.powspctrm           = sum(timelock.jumps(MEGchanindx,:),2);
  data.dimord              = 'chan_freq';
  data.freq                = 1;
  axes(h.TopoMEG);
  ft_topoplotTFR(cfgtopo, data); % FIXME should this not be ft_topoplotER?
  clear data
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function PhiX = polynomial(X, degree, varargin)
    params = inputParser;
    params.addRequired('X');
    params.addRequired('degree');
    params.addParameter('bias', true);
    params.parse(X, degree, varargin{:});
    results = params.Results;
    bias = results.bias;
    [rows, dim] = size(X);
    if bias
        num_features = 1;
    else
        num_features = 0;
    end
    for deg=1:degree
        num_features = num_features + nchoosek(dim+deg-1, deg);
    end
    PhiX = zeros(rows, num_features);
    if bias
        PhiX(:, 1) = 1;
    end
    for r=1:rows
        x = X(r, :);
        if bias
            s = 2;
        else
            s = 1;
        end
        for deg=1:degree
            num_features = nchoosek(dim+deg-1, deg);
            e = s + num_features -1;
            % form all k degree combinations
            combs = spx.discrete.combinations_with_replacement(x, deg);
            features = prod(combs, 2)';
            PhiX(r, s:e) = features;
            s = s + num_features;
        end
    end
end

author: indigits

repo: sparse-plex


function [ price ] = Analytical_Price_Equity_Swap_Cap_ConstRate(r, M, dt, contract_param, sigma )
% Price the equity swap and cap/floor under constant interest rate assumptino

contract_type = contract_param.contract_type; % 1 = swap, 2 = cap/floor
K = contract_param.K;
c = contract_param.c;
T = dt*M;

gamma = (1 + r*dt + c);
mult = (1 - exp(-r*T))/(exp(r*dt) - 1);

if contract_type == 1  % Swap
    price = K*(exp(r*dt) - gamma)*mult;

elseif contract_type == 2  % cap/floor

    C = contract_param.C; % cap
    F = contract_param.F; % floor

    mu = (r - sigma*sigma/2)*dt;
    sigdt = sigma*sqrt(dt);
    a = (log(max(0, F + gamma)) - mu)/sigdt;
    b = (log(C + gamma) - mu)/sigdt;
    t1 = normcdf(b - sigdt) - normcdf(a - sigdt);
    t2 = normcdf(b) - normcdf(a);
    H = F*normcdf(a) + exp(r*dt)*t1 - gamma*t2 + C*(1-normcdf(b));
    price = K*H*mult;
else
    price = -123456789;
end

end


text:
function [bvecs,bvals] = dtiRawBuildBvecs(nVols, xform, gradsFile, bval, outBaseName, assetFlag)
%
% dtiRawBuildBvecs(nvols, xform, [gradsFile=uigetfile], [bval=0.8], ...
%                                 [outBaseName=[]], [assetFlag=false]);
%
% Creates "FSL-style" 'bvecs' and 'bvals' files given a Bammer/GE style
% 'grads' file of gradient directions. The bvecs are reoriented to image
% space given the xform (leave empty for no xform).
%
% xform should be the rotation component of the scanner-to-image xform.
% This transform should be in the NIFTI qto_ijk field, but sometimes it is
% set in sto_ijk. This transform maps scanner coordinates (x,y,z) to
% physical, image coordinates (i,j,k). In other words, it tells us how to
% rotate the scanner space to be aligned with the image space. The bvecs
% are usually specified in the scanner space coordinate frame, so we can
% use this transform to rotate them to the image frame. Note that some
% protocols might keep the directions in image (or 'logical') space, in
% which case the xform was applied at acquisition. For such protocols, you
% should not rotate the bvecs (xform=[] or xform=eye(3)).
%
% If xform is 4x4, it is assumed to be an affine xform and the rotation
% compoenent will be extracted.
%
% assetFlag is an idiosyncratic thing for the Hedehus/Bammer sequence used
% at Stanford. This sequence reorients the gradient directions specified in
% the dwepi.grads file to logical space rather than keeping the directions
% in scanner space. Thus, the bvecs do not need to be reoriented for
% oblique prescriptions as with some other DTI sequences. However, this
% sequence assumes that the 2nd column in dwepi.grads is the phase-encode
% dim. If your phase-enmcode is the usual '2', then this is fine. But, if
% you run ASSET and change the phase encode to L-R (dim 1), you need to
% swap the first and second columns of dwepi.grads. Also, there appears to
% be a flip in the phase-encode dim, so you also need to flip the sign on
% the phase-encode column.
%
% WEB RESOURCES:
%   mrvBrowseSVN('dtiRawBuildBvecs');
%   http://white.stanford.edu/newlm/index.php/DTI_Preprocessing
%
% EXAMPLE USAGE:
%   dwRaw = niftiRead('rawDti.nii.gz');
%   nVols = size(dwRaw.data,4);
%   xform = affineExtractRotation(dwRaw.qto_ijk);
%   [bvecs,bvals] = dtiRawBuildBvecs(nVols, xform, gradsFile, bval, 'rawDti')
%
%
% 2006.12.11 RFD: wrote it with BAW & AJS.
% 2007.01.29 RFD: reverted to using qto_ijk to get rotation matrix rather
% than computing it directly from the quaternion. I think it is safer to
% let the NIFTI code build the matrix for us. Also, we were occasionally
% getting imaginary values in our rotation matrix.
% 2007.03.21 RFD: changed default bval to 0.8 (consistent with our
% preference for msec/micrometers^2 units) and renamed the file to be more
% consistent with the other dti pre-processing functions.
% 2009.05.20 RFD: we now adjust bvalues to reflect gradient strength for
% bvecs with norm != 1.
%
% (C) Stanford University, VISTA
%

%% Check Inputs

if(~exist('outBaseName','var')), outBaseName = ''; end
if(nargout==0&&isempty(outBaseName))
    error('outBaseName must be specificed when nargout==0.');
end
if(~exist('gradsFile','var')||isempty(gradsFile))
  if(isunix)
    defaultDir = '/usr/local/dti/diffusion_grads/';
  else
    defaultDir = pwd;
  end
  [f,p] = uigetfile({'*.grads';'*.*'}, 'Select the GE grads file...',defaultDir);
  if(isnumeric(f)), error('User cancelled.'); end
  gradsFile = fullfile(p,f);
end

if(~exist('bval','var'))
  bval = 0.8;
  warning(sprintf('Using default b-value of %f msec/micrometers^2.',bval));
end

if(~exist('assetFlag','var')||isempty(assetFlag)), assetFlag = false; end

if(~exist('xform','var')); xform = []; end
if(isempty(xform))
  xform = eye(3);
elseif(size(xform,1)==4||size(xform,2)==4)
  xform = affineExtractRotation(xform);
end

%% Build the Bvecs/Bvals

% Read the grads file
grads = dlmread(gradsFile);

% We may need to transpose grads file
if size(grads,1)<size(grads,2)
    grads = grads';
end
if(size(grads,1)<nVols)
  % Assume that there are repeats and that the last repeat might not have
  % finished.
  grads = repmat(grads,ceil(nVols/size(grads,1)),1);
  grads = grads(1:nVols,:);
end

if(assetFlag)
    disp('Applying ASSET reorientation for Hedehus/Bammer sequence...');
    grads = grads*[0 -1 0; 1 0 0; 0 0 1];
end

% X-form the bvecs
bvecs       = xform*grads';
bvecNorm    = sqrt(sum(bvecs.^2));
nz          = bvecNorm~=0;
bvecs(:,nz) = bvecs(:,nz)./repmat(bvecNorm(nz),[3 1]);
bvals       = repmat(bval,1,size(bvecs,2));

% Zero-out non-dwi bvals
bvals(all(bvecs==0)) = 0;

% Scale the bvalues according to gradient magnitude. This assumes that the
% specified b-value is the max bvalue. Note that we only do this when all
% the bvec norms are 1 or less. If any of the bvec norms are >1, then we
% assume that the bvecs are just encoded sloppily and their norms do not
% really reflect the gradient amplitudes.
if(all(bvecNorm<=1))
    bvals = bvals.*bvecNorm.^2;
end

% Write out the bvals and bvecs
if(~isempty(outBaseName))
    dlmwrite([outBaseName '.bvecs'],bvecs,' ');
    dlmwrite([outBaseName '.bvals'],bvals,' ');
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [P loglikelihood] = LearnCPDsGivenGraph(dataset, G, labels)
%
% Inputs:
% dataset: N x 10 x 3, N poses represented by 10 parts in (y, x, alpha)
% G: graph parameterization as explained in PA description
% labels: N x 2 true class labels for the examples. labels(i,j)=1 if the
%         the ith example belongs to class j and 0 elsewhere
%
% Outputs:
% P: struct array parameters (explained in PA description)
% loglikelihood: log-likelihood of the data (scalar)
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

N = size(dataset, 1);
K = size(labels,2);

loglikelihood = 0;
P.c = zeros(1,K);

% estimate parameters
% fill in P.c, MLE for class probabilities
% fill in P.clg for each body part and each class
% choose the right parameterization based on G(i,1)
% compute the likelihood - you may want to use ComputeLogLikelihood.m
% you just implemented.
%%%%%%%%%%%%%%%%%%%%%%%%%
% YOUR CODE HERE
if size(size(G),2) == 2
	G1 = G; G2 = G;
else
	G1 = reshape(G(:,:,1),10,2);
	G2 = reshape(G(:,:,2),10,2);
end
human = dataset([labels(:,1)==1],:,:);
alien = dataset([labels(:,2)==1],:,:);
P.c = [size(human,1),size(alien,1)];
P.c /= N;
for i = 1:10
	if G1(i,1) == 0
		[P.clg(i).mu_y(1),P.clg(i).sigma_y(1)] = FitGaussianParameters(human(:,i,1));
		[P.clg(i).mu_x(1),P.clg(i).sigma_x(1)] = FitGaussianParameters(human(:,i,2));
		[P.clg(i).mu_angle(1),P.clg(i).sigma_angle(1)] = FitGaussianParameters(human(:,i,3));
	else
		% calc beta and sigma
		[Betay,P.clg(i).sigma_y(1)] = FitLinearGaussianParameters(human(:,i,1),reshape(human(:,G1(i,2),:),size(human,1),3));
		[Betax,P.clg(i).sigma_x(1)] = FitLinearGaussianParameters(human(:,i,2),reshape(human(:,G1(i,2),:),size(human,1),3));
		[Betaangle,P.clg(i).sigma_angle(1)] = FitLinearGaussianParameters(human(:,i,3),reshape(human(:,G1(i,2),:),size(human,1),3));
		P.clg(i).theta(1,:) = [Betay(4),Betay(1:3)',Betax(4),Betax(1:3)',Betaangle(4),Betaangle(1:3)'];

	end
	if G2(i,1) == 0
		[P.clg(i).mu_y(2),P.clg(i).sigma_y(2)] = FitGaussianParameters(alien(:,i,1));
		[P.clg(i).mu_x(2),P.clg(i).sigma_x(2)] = FitGaussianParameters(alien(:,i,2));
		[P.clg(i).mu_angle(2),P.clg(i).sigma_angle(2)] = FitGaussianParameters(alien(:,i,3));
	else
		%calc beta and sigma
		[Betay,P.clg(i).sigma_y(2)] = FitLinearGaussianParameters(alien(:,i,1),reshape(alien(:,G2(i,2),:),size(alien,1),3));
		[Betax,P.clg(i).sigma_x(2)] = FitLinearGaussianParameters(alien(:,i,2),reshape(alien(:,G2(i,2),:),size(alien,1),3));
		[Betaangle,P.clg(i).sigma_angle(2)] = FitLinearGaussianParameters(alien(:,i,3),reshape(alien(:,G2(i,2),:),size(alien,1),3));
		P.clg(i).theta(2,:) = [Betay(4),Betay(1:3)',Betax(4),Betax(1:3)',Betaangle(4),Betaangle(1:3)'];
	end
end
loglikelihood = ComputeLogLikelihood(P,G,dataset);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('log likelihood: %f\n', loglikelihood);

author: anhncs

repo: Probabilistic-Graphical-Models


text:
% RSADJUST - adjust l-values (Ramberg-Schmeiser distribution)
%                with respect to signal mean and variance
%
% Usage: p = rsadjust(l3, l4, m, var, skew)
%
% Input:
%   l3   - value lambda3 for Ramberg-Schmeiser distribution
%   l4   - value lambda4 for Ramberg-Schmeiser distribution
%   m    - mean of the signal distribution
%   var  - variance of the signal distribution
%   skew - skewness of the signal distribution (only the sign of
%          this parameter is used).
%
% Output:
%   l1  - value lambda3 for Ramberg-Schmeiser distribution
%   l2  - value lambda4 for Ramberg-Schmeiser distribution
%   l3  - value lambda3 for Ramberg-Schmeiser distribution (copy
%         from input)
%   l4  - value lambda4 for Ramberg-Schmeiser distribution (copy
%         from input)
%
% Author: Arnaud Delorme, SCCN, 2003
%
% See also: RSFIT, RSGET
%
% Reference: Ramberg, J.S., Tadikamalla, P.R., Dudewicz E.J., Mykkytka, E.F.
%            A probability distribution and its uses in fitting data.
%            Technimetrics, 1979, 21: 201-214.

% Copyright (C) 2003 Arnaud Delorme, SCCN, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [l1,l2,l3,l4] = rsadjust( l3, l4, mu, sigma2, m3);

    % swap l3 and l4 for negative skewness
    % ------------------------------------
    if m3 < 0
        ltmp = l4;
        l4   = l3;
        l3   = ltmp;
    end

    A = 1/(1 + l3) - 1/(1 + l4);
    B = 1/(1 + 2*l3) + 1/(1 + 2*l4) - 2*beta(1+l3, 1+l4);
    C = 1/(1 + 3*l3) - 1/(1 + 3*l4) ...
             - 3*beta(1+2*l3, 1+l4) + 3*beta(1+l3, 1+2*l4);

    % compute l2 (and its sign)
    % ------------------------
    l2 = sqrt( (B-A^2)/sigma2 );
    if m3 == 0, m3 = -0.000000000000001; end
    if (m3*(C - 2*A*B + 2*A^3)) < 0, l2 = -l2; end
    %l22 = ((C - 2*A*B + 2*A^3)/m3)^(1/3) % also equal to l2

    % compute l1
    % ----------
    l1 = mu - A/l2;

    return;

    % fitting table 1 of
    % ---------------
    [l1 l2]  = pdffitsolv2(-.0187,-.0388, 0, 1, 1)
    [l1 l2]  = pdffitsolv2(-.1359,-.1359, 0, 1, 0)

    % sign problem
    [l1 l2]  = pdffitsolv2(1.4501,1.4501, 0, 1)

    % numerical problem for l1
    [l1 l2]  = pdffitsolv2(-.00000407,-.001076, 0, 1, 2)
    [l1 l2]  = pdffitsolv2(0,-.000580, 0, 1, 2)

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab