text:
function msg = mpc_init(s,c)
    apm(s,c,'clear all');

    % load model and data
    apm_load(s,c,'control.apm');
    csv_load(s,c,'control.csv');

    % configure MV / CV
    apm_info(s,c,'MV','Q1');
    apm_info(s,c,'MV','Q2');
    apm_info(s,c,'CV','TC1');
    apm_info(s,c,'CV','TC2');

    % dynamic control
    apm_option(s,c,'apm.imode',6);
    apm_option(s,c,'apm.solver',3);
    apm_option(s,c,'apm.hist_hor',600);

    % tune MV
    % delta MV movement penalty
    apm_option(s,c,'Q1.dcost',1.0e-4);
    apm_option(s,c,'Q2.dcost',1.0e-4);
    % penalize energy use
    apm_option(s,c,'Q1.cost',0.0);
    apm_option(s,c,'Q2.cost',0.0);
    % limit MV movement each cycle
    apm_option(s,c,'Q1.dmax',50);
    apm_option(s,c,'Q2.dmax',50);
    % MV limits
    apm_option(s,c,'Q1.upper',100);
    apm_option(s,c,'Q1.lower',0);
    apm_option(s,c,'Q2.upper',100);
    apm_option(s,c,'Q2.lower',0);

    % tune CV
    % how fast to reach setpoint
    apm_option(s,c,'TC1.tau',10);
    apm_option(s,c,'TC2.tau',10);
    % trajectory type
    apm_option(s,c,'TC1.tr_init',2);
    apm_option(s,c,'TC2.tr_init',2);

    % let optimizer use MV
    apm_option(s,c,'Q1.status',1);
    apm_option(s,c,'Q2.status',1);
    % include CV in objective function
    apm_option(s,c,'TC1.status',1);
    apm_option(s,c,'TC2.status',1);

    % feedback status (whether we have measurements)
    apm_option(s,c,'Q1.fstatus',0);
    apm_option(s,c,'Q2.fstatus',0);
    apm_option(s,c,'TC1.fstatus',1);
    apm_option(s,c,'TC2.fstatus',1);

    % web-viewer option, update every second
    apm_option(s,c,'apm.web_plot_freq',3);

    msg = 'initialization complete';
end

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function engine = hmm_2TBN_inf_engine(bnet, varargin)
% HMM_2TBN_INF_ENGINE Inference engine for DBNs which uses the forwards-backwards algorithm.
% engine = hmm_2TBN_inf_engine(bnet, ...)
%
% The DBN is converted to an HMM with a single meganode, but the observed nodes remain factored.
% This can be faster than jtree if the num. hidden nodes is low, because of lower constant factors.
%
% All hidden nodes must be discrete.
% All observed nodes are assumed to be leaves.
% The parents of each observed leaf are assumed to be a subset of the hidden nodes within the same slice.
% The only exception is if bnet is an AR-HMM, where the parents are assumed to be self in the
% previous slice (continuous), plus all the discrete nodes in the current slice.

%% Optional arguments
%% ndx_type - 'B', 'D', or 'SD', used in marginal_family [ 'SD' ]

ndx_type = 'SD';
ss = bnet.nnodes_per_slice;

% parse optional params
args = varargin;
nargs = length(args);
if nargs > 0
  for i=1:2:nargs
    switch args{i},
     %case 'ndx_type', ndx_type = args{i+1};
     otherwise,
      error(['invalid argument name ' args{i}]);
    end
  end
end

% Stuff to do with speeding up marginal_family
%engine.ndx_type = ndx_type;

[int, engine.persist, engine.transient] = compute_interface_nodes(bnet.intra, bnet.inter);
engine.persist_bitv = zeros(1, ss);
engine.persist_bitv(engine.persist) = 1;

ns = bnet.node_sizes(:);
ns(bnet.observed) = 1;
ns(bnet.observed+ss) = 1;
engine.eff_node_sizes = ns;

% for n=1:ss
%   dom = 1:(2*ss); % domain of xi(:,:,1)
%   fam = family(bnet.dag, n+ss);
%   engine.marg_fam2_ndx_id(n) = add_ndx(dom, fam, ns, ndx_type);

%   dom = 1:ss; % domain of gamma(:,:,1)
%   fam = family(bnet.dag, n);
%   engine.marg_fam1_ndx_id(n) = add_ndx(dom, fam, ns, ndx_type);

%   engine.marg_singleton_ndx_id(n) = add_ndx(dom, n, ns, ndx_type);
% end

for o=bnet.observed(:)'
  %if bnet.equiv_class(o,1) ~= bnet.equiv_class(o,2)
  %  error(['observed node ' num2str(o) ' is not tied'])
  %end
  cs = children(bnet.dag, o);
  if ~isempty(cs)
    error(['observed node ' num2str(o) ' is not allowed children'])
  end
end

[engine.startprob, engine.transprob, engine.obsprob] = dbn_to_hmm(bnet);

% This is where we will store the results between enter_evidence and marginal_nodes
engine.one_slice_marginal = [];
engine.two_slice_marginal = [];

ss = length(bnet.intra);
engine.maximize = [];
engine.evidence = [];
engine.node_sizes = [];

% avoid the need to do bnet_from_engine, which is slow
engine.slice_size = ss;
engine.parents = bnet.parents;

engine.bel = [];
engine = class(engine, 'hmm_2TBN_inf_engine', inf_engine(bnet));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function slices_fig_h = cluster_orthviews_showcenters(cl, myview, overlay, xhairs, sliceorder, bgcolor, varargin)
% :Usage:
% ::
%
%    slices_fig_h = cluster_orthviews_showcenters(cl, myview, [overlay], [xhairs], [order slices flag], [background color], [var args])
%
%    cluster_orthviews_showcenters(cl, 'coronal');
%    cluster_orthviews_showcenters(cl, 'sagittal');
%    cluster_orthviews_showcenters(cl, 'axial');
%
% used in cluster_orthviews_classes
%
% ..
%    tor wager, aug 2006
% ..

    if nargin < 3 || isempty(overlay), overlay = which('scalped_single_subj_T1.img'); end
    if nargin < 4, xhairs = 1; end
    if nargin < 5, sliceorder = 0; end

    whichorth = [];
    for i = 1:length(varargin)
       if ischar(varargin{i})
            switch varargin{i}

                case 'whichorth', whichorth = varargin{i+1};

                otherwise, warning('scn_tools:badInput', ['Unknown input string option:' varargin{i}]);
            end
        end
    end

    % set up orthviews to be OK
    scn_export_spm_window('setup', overlay);
%     set(get(findobj('Tag','Graphics'), 'Children'), 'Color', 'white');
    set(findobj('Type', 'axes', 'Parent', findobj('Tag','Graphics')), 'Color', 'white');
    if xhairs, spm_orthviews('Xhairs', 'on'); end

    % Set window to be equal/constant
    [volInf_tmp, dat] = iimg_read_img(overlay, 2);
    dat = dat(volInf_tmp.wh_inmask);
    spm_orthviews('Window', whichorth, [prctile(dat, 2) prctile(dat, 98)]);

    myviews = {'axial' 'coronal' 'sagittal'};   % for selecting SPM window
    whview = find(strcmp(myviews, myview));
    if isempty(whview), error('myview must be axial, coronal, or sagittal.'); end

    if sliceorder
        % order slices from negative to positive along this dimension
        myviews2 = {'sagittal' 'coronal' 'axial' };  % for selecting coord
        whcoord = strmatch(myview, myviews2) ;
        cen = cat(1, cl.mm_center);
        [cen, wh] = unique(cen(:, whcoord));
        % remove duplicates and sort
        cl = cl(wh);

        % get text string base
        mystr = {'x = ' 'y = ' 'z = '};
        textbase = mystr{whcoord};
    end

    % get optimal number of axes
    num_clusters = length(cl);
    rc = ceil(sqrt(num_clusters));

    if isempty(whichorth)
        axh = get_orth_axishandles;
    else
        axh = get_orth_axishandles_whichorth(whichorth);
    end

    axh = axh(whview);

    slices_fig_h = figure; %create_figure(myview);
    set(slices_fig_h, 'Color', 'k');

    for i = 1:num_clusters
        newax(i) = subplot(rc, rc, i);
        axis off;
     end

    for i = 1:num_clusters
        %cluster_orthviews(cl{i}, colors2(classes(i)), 'overlay', overlay);
        spm_orthviews('Reposition', cl(i).mm_center);

        copyobj(get(axh, 'Children'), newax(i));
        axes(newax(i));
        axis image

        h = findobj(newax(i), 'Type', 'text');
        h= h(1); % kludgy fix for multiple matches ***

        if sliceorder
            % try to set a reasonable font size
            pos = get(newax(i), 'Position');
            height = pos(3);
            fs = round(height * 70);
            set(h, 'FontSize', fs)
            % set position of text (move down and right)
            pos = get(h, 'Position');
            pos(2) = 0; %pos(2) - fs./2;
            pos(1) = 0;
            set(h, 'Position', pos)
            % set text string based on cen
            set(h, 'String', [textbase num2str(cen(i))]);
        elseif ishandle(h)
            delete(h);
        end
    end

    % try to set a reasonable enlargement factor
    n = 1 + .15 * log(num_clusters ./ 2);
    n = max(n, 1); n = min(n, 3);
    enlarge_axes(gcf, n)

    % set background color to print as whatever it is on screen
    set(gcf,'InvertHardcopy', 'off');
end

function axish = get_orth_axishandles

    % Get figure and axis handles
    fh = findobj('Tag', 'Graphics');
    ch = get(fh, 'Children');
    for i= 1:length(ch)
        mytype = get(ch(i), 'Type');
        wh(i) = strcmp(mytype, 'axes');
    end
    axish = ch(find(wh));

    if isempty(axish), error('SPM figure orthviews do not exist'); end

    % get which axis is which
    for i = 1:length(axish)
        poss(i, :) = get(axish(i), 'Position');
    end

    % get rid of extra axes we may have created in the 4th quadrant
    other_axes = find(any(poss(:, 1) > .45 & poss(:, 2) < .2, 2));
    axish(other_axes) = [];
    poss(other_axes, :) = [];

    % sort into order:  axial, coronal, saggital
    ssum = sum(poss(:, 1:2), 2);
    [ssum, ind] = sort(ssum);
    axish = axish(ind);

end

function axish = get_orth_axishandles_whichorth(whichorth)
% Get the axis handles for the current orthviews
global st
for i = 1:length(st.vols), wh(i) = ~isempty(st.vols{i}); end
wh = find(wh); wh = wh(whichorth);
axish = cat(1, st.vols{wh}.ax{:});
axish = sort(cat(1, axish(:).ax));

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function determ = frank_determinant ( n )

%*****************************************************************************80
%
%% FRANK_DETERMINANT returns the determinant of the FRANK matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real DETERM, the determinant.
%
  determ = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


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

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

close all

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

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

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

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

%adjust 3D view as desired
adjust_view(robot)

%there are just 2 solutions for these robots and 4 DOF
%q = [pi/2 0.2 0.8 pi/4]
%q = [-pi/4 pi/2 0.5 pi]
%robot=load_robot('kuka', 'KR5_scara_R350_Z200'); n_solutions = 2;
%robot=load_robot('example', 'scara'); n_solutions = 2;
%robot=load_robot('example', '2dofplanar'); n_solutions = 2;
%robot=load_robot('example', '3dofplanar'); n_solutions = 2;
%robot=load_robot('example', 'prismatic');n_solutions = 1; %just one possible solutions for this case

%draw the robot
drawrobot3d(robot, q)

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

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

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

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

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

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

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

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

    pause(1);

    k=sum(sum((T-Ti).^2));
    if k < 0.01 % a simple threshold to find differences in the solution
        correct(1,i)= 1;
    else
        correct(1,i)= 0; %uncorrect solution
        fprintf('\nERROR: One of the solutions seems to be uncorrect. Sum of errors: %f', i, k);
    end
end
%Display a message if any of the solutions is not correct
if sum(correct)==n_solutions
    fprintf('\nOK: Every solution in qinv yields the same position/orientation T');
else
    fprintf('\nERROR: One or more of the solutions seems to be uncorrect.');
end

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

repo: ARTE


text:
function value = r8mat_min ( m, n, a )

%*****************************************************************************80
%
%% R8MAT_MIN returns the minimum entry of an M by N R8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in A.
%
%    Input, integer N, the number of columns in A.
%
%    Input, real A(M,N), the matrix.
%
%    Output, real VALUE, the minimum entry of A.
%
  value = min ( min ( a(1:m,1:n) ) );

  return
end


text:
function variance = chi_square_variance ( a )

%*****************************************************************************80
%
%% CHI_SQUARE_VARIANCE returns the variance of the central Chi squared PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, the parameter of the distribution.
%    1 <= A.
%
%    Output, real VARIANCE, the variance of the PDF.
%
  variance = 2.0 * a;

  return
end


text:
function hsvargplvmClassVisualise(call)

% LVMCLASSVISUALISE Callback function for visualising data.
% FORMAT
% DESC contains the callback functions for visualizing points from the
% latent space in the higher dimension space.
% ARG call : either 'click', 'move', 'toggleDynamics',
% 'dynamicsSliderChange'
%
% COPYRIGHT : Neil D. Lawrence, 2003, 2009
%
% SEEALSO : lvmResultsDynamic

% MLTOOLS

global visualiseInfo

% In which output modality do we predict? If there is only one modality,
% this should be set to either "-1" or 1:model.layer{1}.M
visIndex = visualiseInfo.model.vis.index;
% From which input layer do we predict all the way until the leaves?
visLayer = visualiseInfo.model.vis.layer;

switch call
 case 'click'
  [x, y]  = localCheckPointPosition(visualiseInfo);
  if ~isempty(x)
    visualiseInfo.latentPos(visualiseInfo.dim1) = x;
    visualiseInfo.latentPos(visualiseInfo.dim2) = y;
  end
  visualiseInfo.clicked = ~visualiseInfo.clicked;
  if isfield(visualiseInfo.model, 'dynamics') & ~isempty(visualiseInfo.model.dynamics)
    if visualiseInfo.runDynamics
      visualiseInfo.dynamicsRunning = 1;
      fhandle = str2func([visualiseInfo.model.type 'DynamicsRun']);
      feval(fhandle);
      visualiseInfo.dynamicsRunning = 0;
    end
  else
    visualiseInfo.dynamicsRunning = 0;
  end
 case 'move'
  if visualiseInfo.clicked & ~visualiseInfo.runDynamics
    [x, y]  = localCheckPointPosition(visualiseInfo);
    if ~isempty(x)
      % This should be changed to a model specific visualisation.
      visualiseInfo.latentPos(visualiseInfo.dim1) = x;
      visualiseInfo.latentPos(visualiseInfo.dim2) = y;
      set(visualiseInfo.latentHandle, 'xdata', x, 'ydata', y);
      fhandle = str2func([visualiseInfo.model.type 'PosteriorMeanVar']);
      % In contrast to the other posteriorMeanVar functions, this one
      % has outputs that act as inputs in the next layer. In the current implementation, only
      % one output can be visualised. This has to be instructed in a
      % special field "vis" of the complete model.
      visIndex = visualiseInfo.model.vis.index;
      visLayer = visualiseInfo.model.vis.layer;
      if visualiseInfo.showVariance
        [mu, varsigma] = fhandle(visualiseInfo.model, visualiseInfo.latentPos, [], visLayer, 1, visIndex);
      else
          mu = fhandle(visualiseInfo.model, visualiseInfo.latentPos, [], visLayer, 1, visIndex);
          varsigma = zeros(size(mu));
      end
      if isfield(visualiseInfo.model, 'noise')
        Y = noiseOut(visualiseInfo.model.noise, mu, varsigma);
      else
        Y = mu;
      end
       % Some skeletons need to have their data remapped
      if isfield(visualiseInfo.model, 'Ytochannels') && visualiseInfo.model.Ytochannels
        Y = Ytochannels(Y);
      end
      visualiseInfo.visualiseModify(visualiseInfo.visHandle, ...
                                    Y, visualiseInfo.varargin{:});
    end
  end
 case 'toggleDynamics'
  visualiseInfo.runDynamics = ~visualiseInfo.runDynamics;
  set(visualiseInfo.dynamicsRadio, 'value', visualiseInfo.runDynamics);

 case 'dynamicsSliderChange'
  X = modelOut(visualiseInfo.model.dynamics, get(visualiseInfo.dynamicsSlider, 'value'));
  x = X(1);
  y = X(2);
  visualiseInfo.latentPos(visualiseInfo.dim1) = x;
  visualiseInfo.latentPos(visualiseInfo.dim2) = y;
  set(visualiseInfo.latentHandle, 'xdata', x, 'ydata', y);
  fhandle = str2func([visualiseInfo.model.type 'PosteriorMeanVar']);
  [mu, varsigma] = fhandle(visualiseInfo.model, visualiseInfo.latentPos, [], visLayer, 1, visIndex);
  if isfield(visualiseInfo.model, 'noise')
    Y = noiseOut(visualiseInfo.model.noise, mu, varsigma);
  else
    Y = mu;
  end
  % Some skeletons need to have their data remapped
  if isfield(visualiseInfo.model, 'Ytochannels') && visualiseInfo.model.Ytochannels
      Y = Ytochannels(Y);
  end
  visualiseInfo.visualiseModify(visualiseInfo.visHandle, ...
                                Y, visualiseInfo.varargin{:});

 case 'latentSliderChange'
  counter = 0;
  for i = size(visualiseInfo.latentPos, 2):-1:1
    % Go through setting latent positions to slider positions.
    if i ~= visualiseInfo.dim1 && i ~= visualiseInfo.dim2
      counter = counter + 1;
      visualiseInfo.latentPos(i) = get(visualiseInfo.latentSlider(counter), 'value');
      set(visualiseInfo.sliderTextVal(counter), 'string', num2str(visualiseInfo.latentPos(i)));
    end
  end
  if visualiseInfo.showVariance
    lvmSetPlot;
  else
      lvmSetPlotNoVar('hsvargplvmClassVisualise');
  end
  fhandle = str2func([visualiseInfo.model.type 'PosteriorMeanVar']);
  if visualiseInfo.showVariance
    [mu, varsigma] = fhandle(visualiseInfo.model, visualiseInfo.latentPos, [], visLayer, 1, visIndex);
  else
      mu = fhandle(visualiseInfo.model, visualiseInfo.latentPos, [], visLayer, 1, visIndex);
      varsigma = zeros(size(mu));
  end
  if isfield(visualiseInfo.model, 'noise')
    Y = noiseOut(visualiseInfo.model.noise, mu, varsigma);
  else
    Y = mu;
  end
  % Some skeletons need to have their data remapped
  if isfield(visualiseInfo.model, 'Ytochannels') && visualiseInfo.model.Ytochannels
      Y = Ytochannels(Y);
  end
  visualiseInfo.visualiseModify(visualiseInfo.visHandle, ...
                                    Y, visualiseInfo.varargin{:});

  visualiseInfo.latentHandle = line(visualiseInfo.latentPos(visualiseInfo.dim1), ...
                                    visualiseInfo.latentPos(visualiseInfo.dim2), ...
                                    'markersize', 20, 'color', [0.5 0.5 0.5], ...
                                    'marker', '.', 'visible', 'on', ...
                                    'erasemode', 'xor');

 case 'updateLatentRepresentation'
  visualiseInfo.dim1 = get(visualiseInfo.xDimension, 'value');
  visualiseInfo.dim2 = get(visualiseInfo.yDimension, 'value');
  if visualiseInfo.showVariance
    lvmSetPlot;
  else
     lvmSetPlotNoVar('hsvargplvmClassVisualise');
  end

  visualiseInfo.latentHandle = line(visualiseInfo.latentPos(visualiseInfo.dim1), ...
                                    visualiseInfo.latentPos(visualiseInfo.dim2), ...
                                    'markersize', 20, 'color', [0.5 0.5 0.5], ...
                                    'marker', '.', 'visible', 'on', ...
                                    'erasemode', 'xor');

end

function point = localGetNormCursorPoint(figHandle)

point = get(figHandle, 'currentPoint');
figPos = get(figHandle, 'Position');
% Normalise the point of the curstor
point(1) = point(1)/figPos(3);
point(2) = point(2)/figPos(4);

function [x, y] = localGetNormAxesPoint(point, axesHandle)

position = get(axesHandle, 'Position');
x = (point(1) - position(1))/position(3);
y = (point(2) - position(2))/position(4);
lim = get(axesHandle, 'XLim');
x = x*(lim(2) - lim(1));
x = x + lim(1);
lim = get(axesHandle, 'YLim');
y = y*(lim(2) - lim(1));
y = y + lim(1);

function [x, y] = localCheckPointPosition(visualiseInfo)

% Get the point of the cursor
point = localGetNormCursorPoint(gcf);

% get the position of the axes
position = get(visualiseInfo.plotAxes, 'Position');

% Check if the pointer is in the axes
if point(1) > position(1) ...
      & point(1) < position(1) + position(3) ...
      & point(2) > position(2) ...
      & point(2) < position(2) + position(4);

  % Rescale the point according to the axes
  [x y] = localGetNormAxesPoint(point, visualiseInfo.plotAxes);

  % Find the nearest point
else
  % Return nothing
  x = [];
  y = [];
end

author: SheffieldML

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


text:
% deprecated function to install plugin - use plugin_menu instead

% Copyright (C) 2012- Arnaud Delorme
%
% 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 restartEeglabFlag = plugin_extract(type, pluginlist, page)

disp('The function plugin_extract has been deprecated and calls the plugin_menu function instead');
plugin_menu(pluginlist);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function value = dzasum ( n, x, incx )

%*****************************************************************************80
%
%% DZASUM takes the sum of the absolute values of a complex vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for FORTRAN usage,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, pages 308-323, 1979.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, complex X(*), the vector.
%
%    Input, integer INCX, the increment between successive entries of X.
%
%    Output, real VALUE, the sum of the absolute values.
%
  value = sum ( abs ( real ( x(1:incx:1+(n-1)*incx) ) ) ...
              + abs ( imag ( x(1:incx:1+(n-1)*incx) ) ) );

  return
end


text:
% A test script on image sequences, making sure grayscale and 16 bit
% sequences work

clear

% The location executable will depend on the OS
if(isunix)
    executable = '"../../build/bin/FeatureExtraction"';
else
    executable = '"../../x64/Release/FeatureExtraction.exe"';
end

% Input file
in_dirs = {'../../samples/image_sequence',...
    '../../samples/image_sequence_gray', ...
    '../../samples/image_sequence_16bit'};

% Where to store the output
output_dir = './processed_features/';

for i=1:numel(in_dirs)
    command = sprintf('%s -fdir "%s" -out_dir "%s" -verbose', executable, in_dirs{i}, output_dir);

    if(isunix)
        unix(command);
    else
        dos(command);
    end
end

author: TadasBaltrusaitis

repo: OpenFace


text:
% [succeeded,outmatname,trx] = cadabra2ctrax(featname,roiname,moviename,[outmatname],[doflipup],[dofliplr],[rot])
% inputs cadabra mat file output names, loads in the trajectories, and
% converts to ctrax trx format. these trx will be saved to the file
% outmatname (user will be prompted for this name if none given).
% inputs:
% featname: [optional] "feat" output by cadabra. if not input, user will be
% prompted for file.
% roiname: [optional] "roi" output by cadabra. if not input, user will be
% prompted for file.
% moviename: name of movie tracked
% outmatname: [optional] mat file to save trx trajectories to. user will be
% prompted for file name if none given
% doflipud: [optional] whether to flip the trajectories vertically across
% the middle of the video. Default: false.
% dofliplr: [optional] whether to flip the trajectories horizontally across
% the center of the video. Default: false.
% rot: [optional] degrees to rotate the trajectories by around the center
% of the video. Default: 0.

function [succeeded,outmatname,trx] = cadabra2ctrax(featname,roiname,moviename,outmatname,doflipud,dofliplr,rot,movieinfo)

% initialize outputs
succeeded = false;
if ~exist('outmatname','var'),
  outmatname = '';
end
trx = [];

if ~exist('featname','var') || isempty(featname),
  helpmsg = 'Choose CADABRA feat mat file';
  [featname,featpath] = uigetfilehelp('*.mat','Choose feat file','','helpmsg',helpmsg);
  if ~ischar(featname),
    return;
  end
  featname = fullfile(featpath,featname);
end

if ~exist('roiname','var') || isempty(roiname),
  helpmsg = sprintf('Choose CADABRA roi mat file corresponding to feat file %s',featname);
  roiname = strrep(featname,'feat','roi');
  [roiname,roipath] = uigetfilehelp('*.mat','Choose feat file',roiname,'helpmsg',helpmsg);
  if ~ischar(roiname),
    return;
  end
  roiname = fullfile(roipath,roiname);
end

if ~exist('moviename','var') || isempty(moviename),
  helpmsg = sprintf('Choose movie file corresponding to feat file %s',featname);
  moviename = strrep(featname,'_feat.mat','.avi');
  [moviename,moviepath] = uigetfilehelp('*.avi','Choose movie file',moviename,'helpmsg',helpmsg);
  if ~ischar(moviename),
    return;
  end
  moviename = fullfile(moviepath,moviename);
end

if ~exist('doflipud','var'),
  doflipud = false;
end

if ~exist('dofliplr','var'),
  dofliplr = false;
end

if ~exist('rot','var'),
  rot = 0;
end

if rot ~= 0,
  phi = rot*pi/180;
  R = [cos(rot),sin(rot);-sin(rot),cos(rot)];
end

% load in data
feat = load(featname);
roi = load(roiname);

if ~exist('movieinfo','var') || isempty(movieinfo),

  % read fps, frame size, moviename is either the name of the movie or an mmreader object
  if ~ischar(moviename),
    readerobj = moviename;
    moviename = fullfile(get(readerobj,'Path'),get(readerobj,'Name'));
  else
    readerobj = mmreader(moviename);
  end
  fps = get(readerobj,'FrameRate')+0;
  movieheight = get(readerobj,'Height')+0;
  moviewidth = get(readerobj,'Width')+0;
else
  fps = movieinfo.fps;
  movieheight = movieinfo.height;
  moviewidth = movieinfo.width;
end
movieheight_mm = movieheight * roi.scale.y;
moviewidth_mm = moviewidth * roi.scale.x;

% scale, pxpermm
pxpermm = 1/mean([roi.scale.x,roi.scale.y]);

% allocate
c = cell(1,2);
trx = struct('x',c,'y',c,'theta',c,'a',c,'b',c,...
  'id',c,'moviename',c,'firstframe',c,'arena',c,...
  'nframes',c,'endframe',c,'pxpermm',c,'fps',c,'x_mm',c,...
  'y_mm',c,'a_mm',c,'b_mm',c,'dt',c);
obj = [feat.fly_feat.obj1,feat.fly_feat.obj2];
firstframeoff = feat.fly_feat.frame(1) - 1;
arena = struct('x',nan,'y',nan,'r',nan);
for fly = 1:2,

  % frames for which fly is tracked

  % (x,y) = (0,0) for untracked frames
  badframes = obj(fly).pos_x == 0 & obj(fly).pos_y == 0;
  lastframe = find(~badframes,1,'last');
  firstframe = find(~badframes,1);
  if isempty(lastframe),
    firstframe = 1;
    lastframe = 0;
  end
  nframes = lastframe - firstframe + 1;

  % allocate
  trx(fly).x = nan(1,nframes);
  trx(fly).y = nan(1,nframes);
  trx(fly).theta = nan(1,nframes);
  trx(fly).a = nan(1,nframes);
  trx(fly).b = nan(1,nframes);
  trx(fly).xwingl = nan(1,nframes);
  trx(fly).ywingl = nan(1,nframes);
  trx(fly).xwingr = nan(1,nframes);
  trx(fly).ywingr = nan(1,nframes);
  trx(fly).x_mm = nan(1,nframes);
  trx(fly).y_mm = nan(1,nframes);
  trx(fly).a_mm = nan(1,nframes);
  trx(fly).b_mm = nan(1,nframes);
  trx(fly).dt = nan(1,nframes-1);

  % store parameters
  trx(fly).moviename = moviename;
  trx(fly).firstframe = firstframe + firstframeoff;
  trx(fly).endframe = lastframe + firstframeoff;
  trx(fly).nframes = nframes;
  trx(fly).pxpermm = pxpermm;
  trx(fly).fps = fps;
  trx(fly).id = fly;
  trx(fly).arena = arena;
  %trx(fly).f2i = @(f) f - trx(fly).firstframe + 1;

  % store data
  idx = feat.fly_feat.frame(firstframe:lastframe)-trx(fly).firstframe + 1;
  % we use the quarter major, minor axis length
  trx(fly).a_mm(idx) = obj(fly).FLength/4;
  % maybe area is major/2 * minor/2 * pi, store quarter minor axis length
  trx(fly).b_mm(idx) = (obj(fly).FArea./(obj(fly).FLength/2)/pi)/2;
  % convert degrees to radians
  trx(fly).theta(idx) = obj(fly).headdir*pi/180;

  % convert mm to px, incorporate offset
  trx(fly).x(idx) = obj(fly).pos_x*pxpermm + roi.ROI.cols(1) - 1;
  trx(fly).y(idx) = obj(fly).pos_y*pxpermm + roi.ROI.rows(1) - 1 - 1;
  trx(fly).x_mm = trx(fly).x / pxpermm;
  trx(fly).y_mm = trx(fly).y / pxpermm;
  trx(fly).a = trx(fly).a_mm*pxpermm;
  trx(fly).b = trx(fly).b_mm*pxpermm;

  phil = obj(fly).phil*pi/180;
  wingll = obj(fly).wingll*pxpermm;
  xwl = trx(fly).x(idx) + wingll.*cos(-phil-trx(fly).theta(idx)-pi);
  ywl = trx(fly).y(idx) + wingll.*sin(-phil-trx(fly).theta(idx));
  phir = obj(fly).phir*pi/180;
  winglr = obj(fly).winglr*pxpermm;
  xwr = trx(fly).x(idx) + winglr.*cos(phir-trx(fly).theta(idx)-pi);
  ywr = trx(fly).y(idx) + winglr.*sin(phir-trx(fly).theta(idx));
  trx(fly).xwingl(idx) = xwl;
  trx(fly).ywingl(idx) = ywl;
  trx(fly).xwingr(idx) = xwr;
  trx(fly).ywingr(idx) = ywr;
  trx(fly).xwingl_mm = trx(fly).xwingl / pxpermm;
  trx(fly).ywingl_mm = trx(fly).ywingl / pxpermm;
  trx(fly).xwingr_mm = trx(fly).xwingr / pxpermm;
  trx(fly).ywingr_mm = trx(fly).ywingr / pxpermm;

  trx(fly).dt(idx(1:end-1)) = diff(feat.fly_feat.time(firstframe:lastframe));

  % flipud if necessary
  if doflipud,
    trx(fly).y_mm = movieheight_mm - trx(fly).y_mm;
    trx(fly).y = movieheight - trx(fly).y;
    trx(fly).theta = -trx(fly).theta;
  end

  % fliplr if necessary
  if dofliplr,
    trx(fly).x_mm = moviewidth_mm - trx(fly).x_mm;
    trx(fly).x = moviewidth - trx(fly).x;
    trx(fly).theta = modrange(pi-trx(fly).theta,-pi,pi);
  end

  % rotate if necessary
  if rot ~= 0,
    trx(fly).x_mm = (trx(fly).x_mm - moviewidth_mm/2)*R + moviewidth_mm/2;
    trx(fly).y_mm = (trx(fly).y_mm - movieheight_mm/2)*R + movieheight_mm/2;
    trx(fly).x = (trx(fly).x - moviewidth/2)*R + moviewidth/2;
    trx(fly).y = (trx(fly).y - movieheight/2)*R + movieheight/2;
    trx(fly).theta = modrange(trx(fly).theta + phi,-pi,pi);
  end

end

if ~exist('outmatname','var') || isempty(outmatname),
  [pathstr,name] = fileparts(featname);
  outname = strrep(name,'_feat','');
  outmatname = fullfile(pathstr,[outname,'_trx.mat']);
  helpmsg = {};
  helpmsg{1} = 'Choose mat file to save Ctrax version of trajectories loaded from:';
  helpmsg{2} = sprintf('CADABRA feat file: %s',featname);
  helpmsg{3} = sprintf('CADABRA roi file: %s',roiname);
  [outmatname,outmatpath] = uiputfilehelp('*.mat',sprintf('Save ctrax version of %s',name),outmatname,'helpmsg',helpmsg);
  if ~ischar(outmatname),
    outmatname = '';
    return;
  end
  outmatname = fullfile(outmatpath,outmatname);
  fprintf('Saving to %s...\n',outmatname);
end

save(outmatname,'trx');
succeeded = true;

author: kristinbranson

repo: JAABA


text:
%% Configuration
%
%%
% MTEX is higly customizable. On a local level this can be done by passing
% <GeneralConceptsOptions options> to MTEX functions. On a global level
% this is achieved by editing the file <matlab:edit('mtex_settings.m')
% mtex_settings.m> via the command
%
%   edit mtex_settings
%
% There the following behaviour can be customized
%
% * the orientation of the x and y spatial coordinate axes
% * the alignment of the a and b crystallographic unit cell axes
% * default font size
% * default figure size
% * show or hide coordinates / micronbar on EBSD maps
% * default pole figure annotations
% * default Euler angle convention
% * default color map
% * file extensions associated with EBSD and pole figure files
% * path to CIF (Crystallographic Information Framework) and other files
% * the default maximum iteration depth of the function
% <PoleFigure.calcODF.html calcODF>
% * the amount of available memory
%
% and many more.
%
% You can use
%
%   getMTEXpref
%
% to display all mtex_settings.

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function test_suite=test_config
% tests for cosmo_config
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function s=randstr()
    s=char(ceil(rand(1,20)*20+64));

function test_config_unable_to_find_file()
    fn=tempname();
    assertExceptionThrown(@()cosmo_config(fn),'');

function test_config_unable_to_open_file()
    dirname=fileparts(mfilename('fullpath'));
    assertExceptionThrown(@()cosmo_config(dirname),'');

function test_config_reading_empty()
    helper_write_read_config(false);

function test_config_reading_with_paths()
    helper_write_read_config(true);

function s=helper_generate_config(n_keys,include_keys)
    s=struct();
    for k=1:n_keys
        if k<=numel(include_keys)
            key=include_keys{k};
            value=pwd();
        else
            key=randstr();
            value=randstr();
        end
        s.(key)=value;
    end

function tmp_fn=helper_write_config(c, suffix)
    tmp_fn=tempname();
    fid=fopen(tmp_fn,'w');
    file_closer=onCleanup(@()fclose(fid));

    key_value_pairs_cell=cellfun(@(key)sprintf('%s=%s\n',key,c.(key)),...
                                fieldnames(c),...
                                'UniformOutput',false);
    fprintf(fid,'%s',key_value_pairs_cell{:});
    fprintf(fid,'%s',suffix);

function test_error_when_quote_in_paths()
    orig_warning_state=cosmo_warning();
    warning_state_resetter=onCleanup(@()cosmo_warning(orig_warning_state));
    cosmo_warning('off');

    for with_path_suffix=[true false]
        key=randstr();
        if with_path_suffix
            key=[key '_path'];
        end

        around_chars={'','''','"'};
        for k=1:numel(around_chars)
            c=around_chars{k};

            value=randstr();
            config=struct();
            config.(key)=[c value c];

            fn=helper_write_config(config,'');
            cleaner=onCleanup(@()delete(fn));

            f_handle=@()cosmo_config(fn);
            if numel(c)>0
                assertExceptionThrown(f_handle,'');
            else
                f_handle();
            end

            clear cleaner
        end
    end

function helper_write_read_config(include_path_settings, config)
    orig_warning_state=cosmo_warning();
    warning_state_resetter=onCleanup(@()cosmo_warning(orig_warning_state));

    cosmo_warning('reset');
    cosmo_warning('off');

    n_keys=ceil(rand()*5+5);

    if include_path_settings
        include_keys={'tutorial_data_path','output_data_path'};
    else
        include_keys={};
    end

    % generate and write config
    if nargin<=2
        config=helper_generate_config(n_keys,include_keys);
    end

    tmp_fn=helper_write_config(config, '# this is a comment');
    file_deleter=onCleanup(@()delete(tmp_fn));

    % read configuration
    [c,read_tmp_fn]=cosmo_config(tmp_fn);
    assertEqual(c,config);
    assertEqual(read_tmp_fn,tmp_fn);

    w=cosmo_warning();
    if include_path_settings
        % no warnign shown
        assert(isempty(w.shown_warnings));
    else
        % warning must have been shown if not include path settings
        assert(numel(w.shown_warnings)>0);
        assert(iscellstr(w.shown_warnings))
    end

    % add one key-value pair
    c2=c;
    key=randstr();
    c2.(key)=randstr();

    % write and read new configuration
    cosmo_config(tmp_fn,c2);
    c3=cosmo_config(tmp_fn);

    % should be different with missing key-value pair
    assert(~isequal(c,c3));

    % should be the same
    assertEqual(c2,c3);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function view = LoadLaminarIndices(view)

% view = LoadLaminarIndices(view);
%
% Load the laminar indices, if they are available as a file in the current
% view directory. If not, offer to calculate them. For backward style
% conformance, this takes an optional volume view input, and returns the
% volume with laminar coordinate indices field "laminarIndices" loaded.
% However, if the view is not supplied, the currently selected VOLUME is
% operated upon and returned. Returns an empty variable if something goes
% wrong, e.g., the user % decides not to calculate the indices, which takes
% some time.
%
% Ress, 10/05

mrGlobals

if ieNotDefined('view')
  selectedVOLUME = viewSelected('Volume');
  view = VOLUME{selectedVOLUME};
end
if ~strcmp(view.viewType, 'Gray'), view = switch2Gray(view); end
fName = fullfile(viewDir(view), 'laminarIndices.mat');
if exist(fName, 'file')
  data = load(fName);
  VOLUME{selectedVOLUME}.laminarIndices = data.laminarIndices;
else
  yn = questdlg('No laminar coordinates file. Calculate coordinates?');
  if strcmp(yn, 'Yes'), MapLaminarIndices, end
end

if isfield(VOLUME{selectedVOLUME}, 'laminarIndices')
  view = VOLUME{selectedVOLUME};
else
  view = [];
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [IF_phase_unwrap, GDD_phase, MGD_phase]=phase_feature_extraction(sp_complex, sp_delay)

% This function is using for different kinds of phase feature extraction,
% include instantaneous frequency deviation (IF), group delay deviation
% (GDD), modified group delay (MGD), phase distortion (PD), phase
% distortion standard deviation (PDD) and relative phase shift (RPS).
%
% input:
%     sp_complex:   [n*d] n is frame number (time), d is the feature dimension (frequency)
%                   STFT of target waveform x(n) which contain both magnitude and phase information
%     sp_delay:     [n*d] n is frame number, d is the feature dimension
%                   STFT of modified target waveform n*x(n)
%
% output:
%                   [n*d] n is frame number (time), d is the feature dimension (frequency)
%     IF_phase:     instantaneous frequency deviation
%     GDD_phase:    group delay deviation
%     MGD_phase:    modified group delay
%     PD_phase:     phase distortion
%     PDD_phase:    phase distortion standard deviation
%     RPS_phase:    relative phase shift
%

is_plot = 1;

phase = angle(sp_complex);
phase_unwrap = unwrap(phase);

% [~, IF_phase] = gradient(phase); % deviation along time aixs
[~, IF_phase_unwrap] = gradient(phase_unwrap); % deviation along time aixs

GDD_phase = group_delay_feature(sp_complex); % deviation along frequency aixs

MGD_phase = modified_group_delay_raw(sp_complex, sp_delay);

if is_plot
    figure
    subplot(5,1,1)
    imagesc(phase');colormap jet
    title('original phase')

    subplot(5,1,2)
    imagesc(phase_unwrap');colormap jet
    title('original unwraped phase')

    subplot(5,1,3)
    imagesc(IF_phase_unwrap');colormap jet
    title('instantaneous frequency deviation')

    subplot(5,1,4)
    gdd = GDD_phase.^ 0.4;
    imagesc(real(gdd'));colormap jet
    title('group delay deviation')

    subplot(5,1,5)
    imagesc(MGD_phase');colormap jet
    title('modified group delay')
end

% RPS_phase = relative_phase_shift(phase);
%
% PD_phase = phase_distortion(phase);

% PDD_phase = phase_distortion_std_deviation(phase)

author: singaxiong

repo: SignalGraph


text:
function sqi = flatsqi(signal)
%flatSQI Flatline SQI
%
%
% Input:
%   signal:         single channel (F)ECG [1xN double]
%  MIN_AMP:         if the median of the filtered ECG is inferior to MINAMP
%                   then it is likely to be a flatline note the importance of
%                   the units here for the ECG (mV)
%
% Output:
%   sqi:            resulting xSQI for segment
%
% Fetal Extraction Toolbox, version 1.0, February 2014
% Released under the GNU General Public License
%
% Copyright (C) 2014 Fernando Andreotti
% Dresden University of Technology, Institute of Biomedical Engineering
% fernando.andreotti@mailbox.tu-dresden.de
%
% Last updated : 30-06-2016
%
%
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version.
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details.

MIN_AMP = 0.1;
% == Flatline detection
if (sum(abs(signal-median(signal))>MIN_AMP)/length(signal))<0.05
    % this is a flat line
    sqi = 0;
else
    sqi = 1;
end
% if 20% of the samples (or more) have an absolute amplitude which is higher
% than MIN_AMP then we are good to go.

author: fernandoandreotti


text:
function h=m_ruler(posx,posy,varargin)
% M_RULER Draws a distance scalebar for a map
%   M_RULER([X1 X2],Y1) draws a horizontal scale bar between the
%   normalized coordinates (X1,Y1) and (X2,Y1) where both X/Y are
%   in the range 0 to 1.
%   M_RULER(X1,[Y1 Y2]) draws a vertical scalebar
%   M_RULER(...,NINTS) draws the scalebar with NINTS intervals
%   if NINTS is a scalar (default 4).  Distances of each interval are
%   chosen to be 'nice'. If NINTS is a vector it is understood to be
%   the distances to be used (in meters)
%
%   M_RULER(....,'parameter','value',...) lets you specify
%   extra parameter/value pairs in the usual handle-graphics way.
%   'color' and 'fontsize' are probably the most useful, 'tickdir'
%   'in' and 'out' chooses between different styles.
%
%   Probably BEST to call this AFTER M_GRID otherwise placement might
%   seem a bit odd.
%
% WARNING - the scalebar is probably not useful for any global
%           (i.e. whole-world) or even a significant-part-of-the-globe
%           map, but I won't stop you using it. Caveat user!

% R. Pawlowicz rich@eos.ubc.ca  8/Nov/2006
% 7/Dec/11 - Octave 3.2.3 compatibility
%   Nov/17 - some changes in look for new matlab graphics
%   Jan/18 - patch colours must be double not logical! (no longer the same)

global MAP_PROJECTION

if isempty(MAP_PROJECTION)
   disp('No Map Projection initialized - call M_PROJ first!');
   return;
end

nints=4;
fixticks=0;

if length(varargin)>0
 if isnumeric(varargin{1})
   nints=varargin{1};
   varargin=varargin(2:end);
   fixticks=1;
 end
end

gcolor='k';
glinestyle='-';
glinewidth=3;
gfontsize=get(gca,'fontsize');
gfontname=get(gca,'fontname');
gticklen=get(gca,'ticklength'); gticklen=gticklen(1);
gtickdir=get(gca,'tickdir');

if MAP_PROJECTION.newgraphics
 gticklen=gticklen*3;
end

% Parse parameter list for options. I really should do some
% error checking here, but...

k=1;
while k<=length(varargin)
  switch lower(varargin{k}(1:3))
    case 'col'
      gcolor=varargin{k+1};
    case 'lin'
      switch lower(varargin{k}(1:5))
         case 'linew'
           glinewidth=varargin{k+1};
         case 'lines'
           glinestyle=varargin{k+1};
      end
    case 'fon'
       switch lower(varargin{k}(1:5))
         case 'fonts'
           gfontsize=varargin{k+1};
         case 'fontn'
           gfontname=varargin{k+1};
       end
    case 'tic'
       switch lower(varargin{k}(1:5))
         case 'tickl'
           gticklen=varargin{k+1}(1);
         case 'tickd'
           gtickdir=varargin{k+1};
       end
    case {'get','usa'}
      disp('      ''ticklength'',value');
      disp('      ''tickdir'',( ''in'' | ''out'' )');
      disp('      ''color'',colorspec');
      disp('      ''linewidth'', value');
      disp('      ''linestyle'', ( linespec | ''none'' )');
      disp('      ''fontsize'',value');
      disp('      ''fontname'',name');
       return;
    case 'set'
      disp(['      ticklength = ' num2str(gticklen)]);
      disp(['      tickdir = ' gtickdir]);
      disp(['      color = ' gcolor]);
      disp(['      linewidth = ' num2str(glinewidth)]);
      disp(['      linestyle = ' glinestyle]);
      disp(['      fontsize = ' num2str(gfontsize)]);
      disp(['      fontname = ' gfontname]);
       return;
      otherwise
       disp([' Unknown option: ' varargin{k}]);
  end
  k=k+2;
end

% Need earth radius, in m.
erad=6378137; %m (from WGS-84)

if  ( length(posx)==2 && length(posy)==1)
  posy=[posy posy];
  horiz=1;
elseif  ( length(posx)==1 && length(posy)==2)
  posx=[posx posx];
  horiz=0;
end

xlm=get(gca,'xlim');
ylm=get(gca,'ylim');

% Get into screen coords
posx=xlm(1) + posx*diff(xlm);
posy=ylm(1) + posy*diff(ylm);

if diff(xlm)>10 % we are probably already in meters, i.e. UTM
   scfac=1;
else
  scfac=erad;
end

distance=(diff(posx)+diff(posy))*scfac;

niceints=[ 1       2               5         	     ...
           10      20      25      50       75       ...
	       100     200     250     500      750      ...
	       1000    2000    4000    5000              ...
	       10000   20000   25000   50000    75000    ...
	       100000  200000  250000  500000   750000   ...
	       1000000 2000000 2500000 50000000 7500000  ];

if length(nints)==1

  exactint=distance/(nints);
  [dun,I]=min(abs(niceints-exactint));

  if ~fixticks
    nints=fix(distance/niceints(I));
  end
  dist=[0:nints]*niceints(I);

else
  dist=(nints-nints(1));
  nints=length(dist)-1;
end

if max(log10(dist(2:end)))>=3
   numfac=1000;
   units=' km';
else
   numfac=1;
   units= ' m';
end

if horiz

  if strcmp(gtickdir,'in')

    line(posx(1)+[0 dist(end)/scfac],posy(1)+[0 0],'color',gcolor,'linewidth',glinewidth,'linestyle',glinestyle,...
	 'clipping','off','tag','m_ruler_x');

    XX=posx(1)+[dist;dist]/scfac;
    YY=posy(1)+diff(ylm)*gticklen*[-1;1]*ones(1,nints+1);
    if MAP_PROJECTION.IsOctave
      for k=1:size(XX,2)
        line(XX(:,k),YY(:,k),...
	 'color',gcolor,'linewidth',glinewidth/3,'linestyle',glinestyle,...
	 'clipping','off','tag','m_ruler_y');
      end
    else
      line(XX,YY,...
	 'color',gcolor,'linewidth',glinewidth/3,'linestyle',glinestyle,...
	 'clipping','off','tag','m_ruler_y');
    end
  else

    patch(posx(1)+[dist(1:end-1);dist(1:end-1);dist(2:end);dist(2:end)]/scfac,...
           posy(1)+diff(ylm)*gticklen*[-1;1;1;-1]*ones(1,nints),...
	   repmat(MAP_PROJECTION.LARGVAL  ,4,nints),...
           reshape(double(rem(rem(0:nints*3-1,nints),2)==0),1,nints,3),...
	 'linewidth',glinewidth/3,'linestyle',glinestyle,...
         'clipping','off','tag','m_ruler');

  end

  if nints>1
    for k=1:nints+1
      text(posx(1)+dist(k)/scfac,posy(1)-diff(ylm)*gticklen,num2str(dist(k)/numfac), ...
          'fontsize',gfontsize,'fontname',gfontname,...
         'verticalalignment','top','horizontalalignment','center','tag','m_ruler_label');
    end
  %  text(posx(1)+dist(nints+1)/scfac,posy(1)-diff(ylm)*gticklen*2,sprintf('%d km',dist(end)/numfac),...
  %        'fontsize',gfontsize,'fontname',gfontname,...
  %        'verticalalignment','top','horizontalalignment','center','tag','m_ruler_label');
    text(posx(1)+diff(posx)/2,posy(1)+diff(ylm)*gticklen,units,...
          'fontsize',gfontsize,'fontname',gfontname,...
          'verticalalignment','bottom','horizontalalignment','center','tag','m_ruler_label');
  else
    text(posx(1)+mean(dist)/scfac,posy(1)-diff(ylm)*gticklen*2,[num2str(dist(end)/numfac) units],...
          'fontsize',gfontsize,'fontname',gfontname,...
          'verticalalignment','top','horizontalalignment','center','tag','m_ruler_label');
  end

else

  if strcmp(gtickdir,'in')

    line(posx(1)+[0 0],posy(1)+[0 dist(end)/scfac],'color',gcolor,'linewidth',glinewidth,'linestyle',glinestyle,...
	 'clipping','off','tag','m_ruler_x');

    XX=posx(1)+diff(xlm)*gticklen*[-1;1]*ones(1,nints+1);
    YY=posy(1)+[dist;dist]/scfac;

    if MAP_PROJECTION.IsOctave
      for k=1:size(XX,2)
       line(XX(:,k),YY(:,k),...
	 'color',gcolor,'linewidth',glinewidth/3,'linestyle',glinestyle,...
	 'clipping','off','tag','m_ruler_y');
      end
    else
       line(XX,YY,...
	 'color',gcolor,'linewidth',glinewidth/3,'linestyle',glinestyle,...
	 'clipping','off','tag','m_ruler_y');
    end
  else

    patch(posx(1)+diff(xlm)*gticklen*[-1;1;1;-1]*ones(1,nints),...
          posy(1)+[dist(1:end-1);dist(1:end-1);dist(2:end);dist(2:end)]/scfac,...
	   repmat(MAP_PROJECTION.LARGVAL  ,4,nints),...
           reshape(double(rem(rem(0:nints*3-1,nints),2)==0),1,nints,3),...
	 'linewidth',glinewidth/3,'linestyle',glinestyle,...
         'clipping','off','tag','m_ruler');

  end

  if nints>1
    for k=1:nints
      text(posx(1)+diff(xlm)*gticklen*2,posy(1)+dist(k)/scfac,num2str(dist(k)/numfac), ...
          'fontsize',gfontsize,'fontname',gfontname,...
         'verticalalignment','middle','horizontalalignment','left','tag','m_ruler_label');
    end
    text(posx(1)+diff(xlm)*gticklen*2,posy(1)+dist(nints+1)/scfac,[num2str(dist(end)/numfac) units],...
          'fontsize',gfontsize,'fontname',gfontname,...
          'verticalalignment','middle','horizontalalignment','left','tag','m_ruler_label');
  else
    text(posx(1)+diff(xlm)*gticklen*2,posy(1)+mean(dist)/scfac,[num2str(dist(end)/numfac) units],...
          'fontsize',gfontsize,'fontname',gfontname,...
          'verticalalignment','middle','horizontalalignment','left','tag','m_ruler_label');
  end
end

if nargout==0
    clear h
end

author: goGPS-Project


text:
function [EEG, globalTrendOut] = removeGlobalTrend(EEG, globalTrendIn)
% Perform detrending or high pass filtering to remove low frequency
%
% Usage:
%   EEG = removeTrend(EEG)
%   [EEG, detrendOut] = removeTrend(EEG, detrendIn)
%
% Input:
%   EEG               Structure that requires .data and .srate fields
%   detrendIn         Input structure with fields described below
%
% Structure parameters (detrendIn):
%   detrendChannels   Vector of channels to detrend or filter
%                     (default is all channels)
%   detrendType       One of the strings 'high pass', 'linear', or 'none'
%                     indicating type of detrending (default is'linear')
%   detrendCutoff     Detrend or high pass cutoff (default is 1 Hz)
%   detrendStepSize   Seconds for detrend window slide (default is 0.02 s)
%
% Output:
%   EEG               Revised EEG structure channels detrended or filtered
%   detrendOut        Structure with the following items described below:
%
% Structure parameters (detrendOut):
%   detrendChannels   Vector of detrended or filtered channels
%   detrendType       Type of detrending or filtering
%   detrendCutoff     High pass cutoff or detrend window (default is 1 Hz)
%   detrendStepSize   Seconds for detrend window slide (default is 0.02 s)
%   detrendCommand    String version of detrending command
%
% Implementation notes:
%   1) High pass filtering is done with EEGLAB pop_eegfiltnew FIR filter
%   2) Detrending is done with the chronux_2 runline command
%   3) The EEG.data array will be converted to double regardless of type
%
%% Check the parameters
if nargin < 1 || ~isstruct(EEG)
    error('removeGlobalTrend:NotEnoughArguments', 'first argument must be a structure');
elseif nargin < 2 || ~exist('globalTrendIn', 'var') || isempty(globalTrendIn)
    globalTrendIn = struct();
end
if ~isstruct(globalTrendIn)
    error('removeGlobalTrend:NoData', 'second argument must be a structure')
end

defaults = getPrepDefaults(EEG, 'globaltrend');
globalTrendOut = struct('globalTrendChannels', [], 'doLocal', [], ...
    'localCutoff', [], 'localStepSize', [], 'linearFit', [], ...
     'channelCorrelations', []);
[globalTrendOut, errors] = checkPrepDefaults(globalTrendIn, globalTrendOut, defaults);
if ~isempty(errors)
    error('removeGlobalTrend:BadParameters', ['|' sprintf('%s|', errors{:})]);
end
%% Detrend the data either using high pass or linear detrending
if globalTrendOut.doLocal
    params = struct('detrendChannels', globalTrendOut.globalTrendChannels, ...
                    'detrendType', 'linear', ...
                    'detrendCutoff', globalTrendOut.localCutoff, ...
                    'detrendStepSize', globalTrendOut.localStepSize);

    [EEG, trend] = removeTrend(EEG, params);
    globalTrendOut.localTrend = trend;
end

channels = 1:size(EEG.data, 1);
globalTrendOut.globalTrendChannels = channels;
t = 0:(size(EEG.data, 2) - 1);
t = t/EEG.srate;
data = EEG.data(channels, :);
myPval = zeros(size(data, 1), 2);
myCorr = zeros(size(data, 1), 1);
parfor k = 1:size(data, 1)
  myPval(k, :) = polyfit(t, data(k, :), 1);
  myCorr(k) = corr(t', (data(k, :))');
  data(k, :) = data(k, :) - polyval(myPval(k, :), t);
end
globalTrendOut.linearFit = myPval;
globalTrendOut.channelCorrelations = myCorr;
EEG.data(channels, :) = data;

author: VisLab

repo: EEG-Clean-Tools


text:
% STD_SPEC - Returns the data or ICA component spectra for a dataset. Updates the EEG structure
%              in the Matlab environment and in the .set file as well. Saves the spectra
%              in a file.
% Usage:
%           >> [spec freqs] = std_spec(EEG, 'key', 'val', ...);
%
%              Computes the mean spectra of the data channels or activities of specified
%              components of the supplied dataset. The spectra are saved in a Matlab file.
%              If such a file already exists, loads the spectral information from this file.
%              Options (below) specify which components to use, and the desired frequency
%              range. There is also an option to specify other SPECTOPO input variables
%              (see >> help spectopo for details).
%
%              Returns the removed mean spectra of the selected ICA components in the
%              requested frequency range. If the spectra were computed previously but a
%              different frequency range is selected, there is an overwrite option.
%              so. The function will load previously computed log spectra, if any, and
%              will remove the mean from the requested frequency range. The frequencies
%              vector is also returned.
% Inputs:
%   EEG - a loaded epoched EEG dataset structure.
%
% Optional inputs:
%   'components' - [numeric vector] components of the EEG structure for which
%                  activation spectrum will be computed. Note that because
%                  computation of ERP is so fast, all components spectrum are
%                  computed and saved. Only selected component
%                  are returned by the function to Matlab
%                  {default|[] -> all}
%   'channels'   - [cell array] channels of the EEG structure for which
%                  activation spectrum will be computed. Note that because
%                  computation of ERP is so fast, all channels spectrum are
%                  computed and saved. Only selected channels
%                  are returned by the function to Matlab
%                  {default|[] -> none}
%   'recompute'  - ['on'|'off'] force recomputing ERP file even if it is
%                  already on disk.
%   'trialindices' - [cell array] indices of trials for each dataset.
%                  Default is all trials.
%   'recompute'  - ['on'|'off'] force recomputing data file even if it is
%                  already on disk.
%   'rmcomps'    - [integer array] remove artifactual components (this entry
%                  is ignored when plotting components). This entry contains
%                  the indices of the components to be removed. Default is none.
%   'interp'     - [struct] channel location structure containing electrode
%                  to interpolate ((this entry is ignored when plotting
%                  components). Default is no interpolation.
%   'output'     - ['power'|'fft'] compute power of keep single complex
%                  'fft' estimate. Default is 'power'.
%   'fileout'    - [string] name of the file to save on disk. The default
%                  is the same name (with a different extension) as the
%                  dataset given as input.
%  'savetrials'  - ['on'|'off'] save single-trials ERSP. Requires a lot of disk
%                  space (dataset space on disk times 10) but allow for refined
%                  single-trial statistics.
%
% spectrum specific optional inputs:
%   'specmode'   - ['psd'|'fft'|'pburg'|'pmtm'] method to compute spectral
%                  decomposition. 'psd' uses the spectopo function (optional
%                  parameters to this function may be given as input). 'fft'
%                  uses a simple fft on each trial. For continuous data
%                  data trials are extracted automatically (see 'epochlim'
%                  and 'epochrecur' below). Two experimental modes are
%                  'pmtm' and 'pbug' which use multitaper and the Burg
%                  method to compute spectrum respectively. NOTE THAT SOME
%                  OF THESE OPTIONS REQUIRE THE SIGNAL PROCESSING TOOLBOX.
%   'epochlim'   - [min max] for FFT on continuous data, extract data
%                  epochs with specific epoch limits in seconds (see also
%                  'epochrecur' below). Default is [0 1].
%   'epochrecur' - [float] for FFT on continuous data, set the automatic
%                  epoch extraction recurence interval (default is 0.5 second).
%   'timerange'  - [min max] use data within a specific time range before
%                  computing the data spectrum. For instance, for evoked
%                  data trials, it is recommended to use the baseline time
%                  period.
%   'logtrials'  - ['on'|'off'] compute single-trial log transform before
%                  averaging them. Default is 'off' for 'psd' specmode and
%                  'on' for 'fft' specmode. Ignored when output is set to
%                  'fft'.
%   'continuous' - ['on'|'off'] force epoch data to be treated as
%                  continuous so small data epochs can be extracted for the
%                  'fft' specmode option. Default is 'off'.
%   'freqrange'  - [minhz maxhz] frequency range (in Hz) within which to
%                  return the spectrum {default|[]: [0 sample rate/2]}.
%                  Note that this does not affect the spectrum computed on
%                  disk, only the data returned by this function as output.
%   'nw'         - [integer] number of tapers for the 'pmtm' spectral
%                  method. Default is 4.
%   'burgorder'  - [integet] order for the Burg spectral method.
%
% Changes between EEGLAB 13 and later EEGLAB versions:
% For the 'specmode' option 'fft', EEGLAB 14 and later version detrend the
% data and apply hamming taper to it. EEGLAB 13 and earlier remove the
% baseline from the data and apply a hamming taper but only for continuous data.
%
% Other optional spectral parameters:
%   All optional parameters to the spectopo function may be provided to this
%   function as well (requires the 'specmode' option above to be set to
%   'psd').
%
% Outputs:
%   spec      - the mean spectra (in dB) of the requested ICA components in the selected
%               frequency range (with the mean of each spectrum removed).
%   freqs     - a vector of frequencies at which the spectra have been computed.
%
% Files output or overwritten for ICA:
%               [dataset_filename].icaspec,   % raw spectrum of ICA components
% Files output or overwritten for data:
%               [dataset_filename].datspec,
%
% See also  SPECTOPO, STD_ERP, STD_ERSP, STD_MAP, STD_PRECLUST
%
% Authors: Arnaud Delorme, SCCN, INC, UCSD, January, 2005

% Defunct:      0 -> if frequency range is different from saved spectra, ask via a
%                    pop-up window whether to keep existing spectra or to overwrite them.

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, October 11, 2004, arno@sccn.ucsd.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 [X, f, overwrt] = std_spec(EEG, varargin)

overwrt = 1; % deprecated
if nargin < 1
    help std_spec;
    return;
end

% decode inputs
% -------------
if ~isempty(varargin)
    if ~ischar(varargin{1})
        varargin = { varargin{:} [] [] };
        if all(varargin{1} > 0)
            options = { 'components' varargin{1} 'freqrange' varargin{2} };
        else
            options = { 'channels' -varargin{1} 'freqrange' varargin{2} };
        end
    else
        options = varargin;
    end
else
    options = varargin;
end

[g spec_opt] = finputcheck(options, { 'components' 'integer' []         [];
                                      'channels'   'cell'    {}         {};
                                      'timerange'  'float'   []         [];
                                      'specmode'   'string'  {'fft','psd','pmtm','pburg'} 'psd';
                                      'recompute'  'string'  { 'on','off' } 'off';
                                      'savetrials' 'string'  { 'on','off' } 'off';
                                      'continuous' 'string'  { 'on','off' } 'off';
                                      'logtrials'  'string'  { 'on','off' 'notset' } 'notset';
                                      'output'     'string'  { 'power','fft' } 'power';
                                      'savefile'   'string'  { 'on','off' } 'on';
                                      'epochlim'   'real'    []         [0 1];
                                      'trialindices' { 'integer','cell' } []         [];
                                      'epochrecur' 'real'    []         0.5;
                                      'rmcomps'    'cell'    []         cell(1,length(EEG));
                                      'nw'         'float'   []         4;
                                      'fileout'    'string'  []         '';
                                      'trialinfo'  'struct'  []         struct([]);
                                      'burgorder'  'integer' []         20;
                                      'interp'     'struct'  { }        struct([]);
                                      'nfft'       'integer' []         [];
                                      'freqrange'  'real'    []         [] }, 'std_spec', 'ignore');
if ischar(g), error(g); end
if isempty(g.trialindices), g.trialindices = cell(1, length(EEG)); end
if ~iscell(g.trialindices), g.trialindices = { g.trialindices }; end
if ~strcmpi(g.specmode, 'fft') && strcmpi(g.output, 'ftt'), error('FFT option only valid when computing FFT'); end
if isfield(EEG,'icaweights')
   numc = size(EEG(1).icaweights,1);
else
   error('EEG.icaweights not found');
end
if isempty(g.components)
    g.components = 1:numc;
end

EEG_etc = [];

% filename
% --------
if isempty(g.fileout), g.fileout = fullfile(EEG(1).filepath, EEG(1).filename(1:end-4)); end
if ~isempty(g.channels)
    filename = [ g.fileout '.datspec'];
    prefix = 'chan';
else
    filename = [ g.fileout '.icaspec'];
    prefix = 'comp';
end

% SPEC information found in datasets
% ---------------------------------
if exist(filename) && strcmpi(g.recompute, 'off')

    fprintf('File "%s" found on disk, no need to recompute\n', filename);
    if strcmpi(prefix, 'comp')
        [X,tmp,f] = std_readfile(filename, 'components', g.components, 'freqlimits', g.freqrange, 'measure', 'spec');
    else
        [X,tmp,f] = std_readfile(filename, 'channels', g.channels, 'freqlimits', g.freqrange, 'measure', 'spec');
    end
    if ~isempty(X), return; end
end

% No SPEC information found
% -------------------------
if isempty(g.channels)
     [X,boundaries]  = eeg_getdatact(EEG, 'component', [1:size(EEG(1).icaweights,1)], 'trialindices', g.trialindices );
else [X,boundaries]  = eeg_getdatact(EEG, 'trialindices', g.trialindices, 'rmcomps', g.rmcomps, 'interp', g.interp);
end
if ~isempty(boundaries) && boundaries(end) ~= size(X,2), boundaries = [boundaries size(X,2)]; end

% get specific time range for epoched and continuous data
% -------------------------------------------------------
oritrials = EEG.trials;
if ~isempty(g.timerange)
    if oritrials > 1
        timebef  = find(EEG(1).times >= g.timerange(1) & EEG(1).times < g.timerange(2) );
        X        = X(:,timebef,:);
        EEG(1).pnts = length(timebef);
    else
        disp('warning: ''timerange'' option cannot be used with continuous data');
    end
end

% extract epochs if necessary
% ---------------------------
if all([ EEG.trials] == 1) || strcmpi(g.continuous, 'on')
    epochCount  = 1;
    sampleCount = 1;
    for iEEG = 1:length(EEG)
        TMP = EEG(1);
        TMP.data = X;
        TMP.icaweights = [];
        TMP.icasphere  = [];
        TMP.icawinv    = [];
        TMP.icaact     = [];
        TMP.icachansind = [];
        TMP.trials = size(TMP.data,3);
        TMP.pnts   = size(TMP.data,2);
        TMP.event  = [];
        TMP.urevent  = [];
        TMP.epoch  = [];
        TMP.chanlocs = [];
        for index = 1:length(boundaries)
            TMP.event(index).type = 'boundary';
            TMP.event(index).latency = boundaries(index);
        end
        TMP = eeg_checkset(TMP);
        if TMP.trials > 1
            % epoch data - need to re-extract data
            TMP = pop_select(TMP, 'trial', [epochCount:(epochCount+EEG(iEEG).trials-1)]);
            epochCount = epochCount+EEG(iEEG).trials;
            TMP = eeg_epoch2continuous(TMP);
        else
            % continuous data - need to re-extract data
            TMP = pop_select(TMP, 'point', [sampleCount:(sampleCount+EEG(iEEG).pnts-1)]);
            sampleCount = sampleCount+EEG(iEEG).pnts;
        end
        TMP = eeg_regepochs(TMP, g.epochrecur, g.epochlim);
        disp('Warning: continuous data, extracting 1-second epochs');
        if iEEG == 1,
            XX = TMP.data;
            newTrialInfo = g.trialinfo(iEEG);
            newTrialInfo(1:size(TMP.data,3)) = g.trialinfo(iEEG);
        else
            XX(:,:,end+1:end+size(TMP.data,3)) = TMP.data;
            newTrialInfo(end+1:end+size(TMP.data,3)) = g.trialinfo(iEEG);
        end
    end
    g.trialinfo = newTrialInfo;
    X = XX;
end

% compute spectral decomposition
% ------------------------------
if strcmpi(g.logtrials, 'notset'), if strcmpi(g.specmode, 'fft') g.logtrials = 'on'; else g.logtrials = 'off'; end; end
if strcmpi(g.logtrials, 'on'), datatype = 'SPECTRUMLOG'; else datatype = 'SPECTRUMABS'; end
if strcmpi(g.specmode, 'psd')
    if strcmpi(g.savetrials, 'on') || strcmpi(g.logtrials, 'on')
        if all([ EEG.trials] == 1) || strcmpi(g.continuous, 'on')
            if isequal(g.epochlim, [0 1])
                fprintf('Spectopo(psd): randomly extracted epochs are only 1 seconds. PSD is better suited for longer epochs.\n');
            end
        end
        fprintf('Computing spectopo (psd) across trials: ');
        for iTrial = 1:size(X,3)
            [tmp, f] = spectopo(X(:,:,iTrial), size(X,2), EEG(1).srate, 'plot', 'off', 'boundaries', boundaries, 'nfft', g.nfft, 'verbose', 'off', spec_opt{:});
            if iTrial == 1
                XX = zeros(size(tmp,1), size(tmp,2), size(X,3));
            end
            XX(:,:,iTrial) = tmp;
            %if iTrial == 1 && size(X,3) > 1, XX(:,:,size(X,3)) = 0; end
            if mod(iTrial,10) == 0, fprintf('%d ', iTrial); end
        end
        fprintf('\n');
        if strcmpi(g.logtrials, 'off')
             X = 10.^(XX/10);
        else X = XX;
        end
        if strcmpi(g.savetrials, 'off')
            X = mean(X,3);
        end
    else
        [X, f] = spectopo(X, size(X,2), EEG(1).srate, 'plot', 'off', 'boundaries', boundaries, 'nfft', g.nfft, 'verbose', 'off', spec_opt{:});
        X = 10.^(X/10);
    end
elseif strcmpi(g.specmode, 'pmtm')
    if strcmpi(g.logtrials, 'on')
        error('Log trials option cannot be used in conjunction with the PMTM option');
    end
    if all([ EEG.trials ] == 1) && ~isempty(boundaries), disp('Warning: multitaper does not take into account boundaries in continuous data (use ''psd'' method instead)'); end
    fprintf('Computing spectrum using multitaper method:');
    for cind = 1:size(X,1)
        fprintf('.');
        for tind = 1:size(X,3)
            [tmpdat f] = pmtm(X(cind,:,tind), g.nw, g.nfft, EEG.srate);
            if cind == 1 && tind == 1
                X2 = zeros(size(X,1), length(tmpdat), size(X,3));
            end
            X2(cind,:,tind) = tmpdat;
        end
    end
    fprintf('\n');
    X = X2;
    if strcmpi(g.savetrials, 'off'), X = mean(X,3); end
elseif strcmpi(g.specmode, 'pburg')
    if strcmpi(g.logtrials, 'on')
        error('Log trials option cannot be used in conjunction with the PBURB option');
    end
    fprintf('Computing spectrum using Burg method:');
    if all([ EEG.trials ] == 1) && ~isempty(boundaries), disp('Warning: pburg does not take into account boundaries in continuous data (use ''psd'' method instead)'); end
    for cind = 1:size(X,1)
        fprintf('.');
        for tind = 1:size(X,3)
            [tmpdat f] = pburg(X(cind,:,tind), g.burgorder, g.nfft, EEG.srate);
            if cind == 1 && tind == 1
                X2 = zeros(size(X,1), length(tmpdat), size(X,3));
            end
            X2(cind,:,tind) = tmpdat;
        end
    end
    fprintf('\n');
    X = X2;
    if strcmpi(g.savetrials, 'off'), X = mean(X,3); end
else % fft mode
    %
    if size(X,3) > 1
        for iTrial = 1:size(X,3)
            X(:,:,iTrial) = detrend(X(:,:,iTrial)')';
        end
    else
        X = detrend(X')';
    end
    try
        X = bsxfun(@times, X, hamming(size(X,2))'); % apply hamming window even for data trials (not the case in EEGLAB 13)
    catch
        X = bsxfun(@times, X, hamming2(size(X,2))');
    end
    disp('Warning: std_spec function computation has changed since version 13 (see help message)');
    %end
    % if all([ EEG.trials ] == 1) && ~isempty(boundaries), disp('Warning: fft does not take into account boundaries in continuous data (use ''psd'' method instead)'); end
    tmp   = fft(X, g.nfft, 2);
    f     = linspace(0, EEG(1).srate/2, floor(size(tmp,2)/2));
    f     = f(2:end); % remove DC (match the output of PSD)
    tmp   = tmp(:,2:floor(size(tmp,2)/2),:);

    % To compute spectral density (but still need FFT correction
    %     dens  = f(3)-f(2)
    %     tmp   = tmp(:,2:floor(size(tmp,2)/2),:)/dens;

    if strcmpi(g.output, 'power')
        X     = tmp.*conj(tmp);
        if strcmpi(g.logtrials, 'on'),  X = 10*log10(X); end
    else
        X = tmp;
        datatype = 'SPECTRUMFFT';
    end
end
eeglab_options;
if option_single
    X = single(X);
end

% Save SPECs in file (all components or channels)
% -----------------------------------------------
fileNames = computeFullFileName( { EEG.filepath }, { EEG.filename });
if strcmpi(g.savefile, 'on')
    options = { options{:} spec_opt{:} 'timerange' g.timerange 'nfft' g.nfft 'specmode' g.specmode };
    if strcmpi(prefix, 'comp')
        savetofile( filename, f, X, 'comp', 1:size(X,1), options, {}, fileNames, g.trialindices, datatype , g.trialinfo);
    else
        if ~isempty(g.interp)
            savetofile( filename, f, X, 'chan', 1:size(X,1), options, { g.interp.labels }, fileNames, g.trialindices, datatype , g.trialinfo);
        else
            tmpchanlocs = EEG(1).chanlocs;
            savetofile( filename, f, X, 'chan', 1:size(X,1), options, { tmpchanlocs.labels }, fileNames, g.trialindices, datatype , g.trialinfo);
        end
    end
end
return;

% compute full file names
% -----------------------
function res = computeFullFileName(filePaths, fileNames);
for index = 1:length(fileNames)
    res{index} = fullfile(filePaths{index}, fileNames{index});
end

% -------------------------------------
% saving SPEC information to Matlab file
% -------------------------------------
function savetofile(filename, f, X, prefix, comps, params, labels, dataFiles, dataTrials, datatype , trialInfo);

    disp([ 'Saving SPECTRAL file ''' filename '''' ]);
    allspec = [];
    for k = 1:length(comps)
        allspec = setfield( allspec, [ prefix int2str(comps(k)) ], squeeze(X(k,:,:)));
    end
    if nargin > 6 && ~isempty(labels)
        allspec.labels = labels;
    end
    allspec.freqs      = f;
    allspec.parameters = params;
    allspec.datatype   = datatype;
    allspec.datafiles   = dataFiles;
    allspec.datatrials  = dataTrials;
    allspec.trialinfo   = trialInfo;
    allspec.average_spec = mean(X,1);
    std_savedat(filename, allspec);

% -------------------------------------
% Adapted from Octave version
% -------------------------------------
function c = hamming2(m, opt)

if (nargin < 1 || nargin > 2)
    help hamming;
    return;
end

if (~(isscalar (m) && (m == fix (m)) && (m > 0)))
    error ('hamming: M must be a positive integer');
end

N = m - 1;
if (nargin == 2)
    switch (opt)
        case 'periodic'
            N = m;
        case 'symmetric'
            % Default option, same as no option specified.
        otherwise
            error ('hamming: window type must be either "periodic" or "symmetric"');
    end
end

if (m == 1)
    c = 1;
else
    m = m - 1;
    c = 0.54 - 0.46 * cos (2 * pi * (0 : m)' / N);
end


text:
function [out] = grayero(im,se)

if (~isa(se,'strel'))
    se=strel(se);
end
seq=getsequence(se);
lunghezza=size(seq,1);

out=gerod(im,getnhood(seq(1)),2);
for ii=2:lunghezza
    out=gerod(out,getnhood(seq(ii)),2);
end

author: Sable

repo: mcbench-benchmarks


text:
classdef matRad_SquaredUnderdosing < DoseObjectives.matRad_DoseObjective
% matRad_SquaredUnderdosing Implements a penalized squared underdosing objective
%   See matRad_DoseObjective for interface description
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2020 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    properties (Constant)
        name = 'Squared Underdosing';
        parameterNames = {'d^{min}'};
        parameterTypes = {'dose'};
    end

    properties
        parameters = {60};
        penalty = 1;
    end

    methods
        function obj = matRad_SquaredUnderdosing(penalty,dMin)
            %If we have a struct in first argument
            if nargin == 1 && isstruct(penalty)
                inputStruct = penalty;
                initFromStruct = true;
            else
                initFromStruct = false;
                inputStruct = [];
            end

            %Call Superclass Constructor (for struct initialization)
            obj@DoseObjectives.matRad_DoseObjective(inputStruct);

            %now handle initialization from other parameters
            if ~initFromStruct
                if nargin == 2 && isscalar(dMin)
                    obj.parameters{1} = dMin;
                end

                if nargin >= 1 && isscalar(penalty)
                    obj.penalty = penalty;
                end
            end
        end

        %% Calculates the Objective Function value
        function fDose = computeDoseObjectiveFunction(obj,dose)
            % overdose : dose minus prefered dose
            underdose = dose - obj.parameters{1};

            % apply positive operator
            underdose(underdose>0) = 0;

            % claculate objective function
            fDose = obj.penalty/numel(dose) * (underdose'*underdose);
        end

        %% Calculates the Objective Function gradient
        function fDoseGrad   = computeDoseObjectiveGradient(obj,dose)
            % overdose : dose minus prefered dose
            underdose = dose - obj.parameters{1};

            % apply positive operator
            underdose(underdose>0) = 0;

            % calculate delta
            fDoseGrad = 2 * obj.penalty/numel(dose) * underdose;
        end
    end

end

repo: matRad


text:
function test_ft_plot_dipole

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_plot_dipole

figure
ft_plot_dipole([0 0 0], [1 2 3])

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function h = hs71H(x,lambda)

h = [ 2*x(4)             x(4)   x(4)   2*x(1)+x(2)+x(3);
      x(4)               0      0   x(1);
      x(4)               0      0   x(1);
      2*x(1)+x(2)+x(3)  x(1)  x(1)  0 ];

h = h + lambda.ineqnonlin*-[    0      x(3)*x(4) x(2)*x(4) x(2)*x(3);
                            x(3)*x(4)     0     x(1)*x(4) x(1)*x(3);
                            x(2)*x(4) x(1)*x(4)     0     x(1)*x(2);
                            x(2)*x(3) x(1)*x(3) x(1)*x(2)     0  ];

h = h + lambda.eqnonlin * diag([2 2 2 2]);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function J=vl_imwhiten(I,alpha,cutoff)
% VL_IMWHITEN  Whiten image
%   J = VL_IMWHITEN(I,ALPHA) approximatively whitens the power spectrum
%   of the natural image I. The algorithm assumes that the modulus of
%   the spectrum decays as 1/f^ALPHA (f is the frequency).
%
%   VL_IMWHITEN(I) uses ALPHA=1 (a typical value for natural images).
%
%   VL_IMWHITEN(I,ALPHA,CUTOFF) also applies a low-pass filter with
%   cutoff frequency equal to CUTOFF x FN, where FN is the Nyquist
%   frequency (half of the sampling frequency).
%
%   See also: VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

if ~exist('alpha'),  alpha = 1 ; end
if ~exist('cutoff'), cutoff = [] ; end

[M,N]=size(I) ;

% Frequency domain
fn = 0.5 ; % Nyquist freq (=1/2T, T=1)
fx_range=linspace(-fn, fn, N) ;
fy_range=linspace(-fn, fn, M) ;
[fx fy]=meshgrid(fx_range, fy_range) ;

% Whitening filter
rho=sqrt(fx.*fx+fy.*fy);
filt=rho.^alpha ;

% Low-pass filter
if ~isempty(cutoff)
  fcut = cutoff * fn ;
  filt = filt .* exp(-(rho/fcut).^4);
  %filt = filt .* exp( - 0.5 * (rho / fcut) .^ 2);
end

% Apply filter
J = real(ifft2(fft2(I).*fftshift(filt))) ;

author: Cloud-CV

repo: object-proposals


text:
function value = mono_422_3d ( n, x )

%*****************************************************************************80
%
%% MONO_422_3D evaluates X**4 Y**2 Z**2.
%
%  Modified:
%
%    25 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Input, real X(N), the point at which the monomial is to be evaluated.
%
%    Output, real VALUE, the value of the monomial.
%
  value = x(1)^4 * x(2)^2 * x(3)^2;

  return
end


text:
function gd=gabdual(g,a,M,varargin)
%GABDUAL  Canonical dual window of Gabor frame
%   Usage:  gd=gabdual(g,a,M);
%           gd=gabdual(g,a,M,L);
%           gd=gabdual(g,a,M,'lt',lt);
%
%   Input parameters:
%         g     : Gabor window.
%         a     : Length of time shift.
%         M     : Number of channels.
%         L     : Length of window. (optional)
%         lt    : Lattice type (for non-separable lattices).
%   Output parameters:
%         gd : Canonical dual window.
%
%   `gabdual(g,a,M)` computes the canonical dual window of the discrete Gabor
%   frame with window *g* and parameters *a*, *M*.
%
%   The window *g* may be a vector of numerical values, a text string or a
%   cell array. See the help of |gabwin| for more details.
%
%   If the length of *g* is equal to *M*, then the input window is assumed
%   to be an FIR window. In this case, the canonical dual window also has
%   length of *M*. Otherwise the smallest possible transform length is chosen
%   as the window length.
%
%   `gabdual(g,a,M,L)` returns a window that is the dual window for a system
%   of length *L*. Unless the dual window is a FIR window, the dual window
%   will have length *L*.
%
%   `gabdual(g,a,M,'lt',lt)` does the same for a non-separable lattice
%   specified by *lt*. Please see the help of |matrix2latticetype| for a
%   precise description of the parameter *lt*.
%
%   If $a>M$ then the dual window of the Gabor Riesz sequence with window
%   *g* and parameters *a* and *M* will be calculated.
%
%   Examples:
%   ---------
%
%   The following example shows the canonical dual window of the Gaussian
%   window:::
%
%     a=20;
%     M=30;
%     L=300;
%     g=pgauss(L,a*M/L);
%     gd=gabdual(g,a,M);
%
%     % Simple plot in the time-domain
%     figure(1);
%     plot(gd);
%
%     % Frequency domain
%     figure(2);
%     magresp(gd,'dynrange',100);
%
%   See also:  gabtight, gabwin, fir2long, dgt

%   AUTHOR : Peter L. Søndergaard.
%   TESTING: TEST_DGT
%   REFERENCE: REF_GABDUAL.

%% ---------- Assert correct input.

if nargin<3
  error('%s: Too few input parameters.',upper(mfilename));
end;

definput.keyvals.L=[];
definput.keyvals.lt=[0 1];
definput.keyvals.nsalg=0;
[flags,kv,L]=ltfatarghelper({'L'},definput,varargin);

%% ------ step 2: Verify a, M and L
if isempty(L)
    if isnumeric(g)
        % Use the window length
        Ls=length(g);
    else
        % Use the smallest possible length
        Ls=1;
    end;

    % ----- step 2b : Verify a, M and get L from the window length ----------
    L=dgtlength(Ls,a,M,kv.lt);

else

    % ----- step 2a : Verify a, M and get L

    Luser=dgtlength(L,a,M,kv.lt);
    if Luser~=L
        error(['%s: Incorrect transform length L=%i specified. Next valid length ' ...
               'is L=%i. See the help of DGTLENGTH for the requirements.'],...
              upper(mfilename),L,Luser);
    end;

end;

%% ----- step 3 : Determine the window

[g,info]=gabwin(g,a,M,L,kv.lt,'callfun',upper(mfilename));

if L<info.gl
  error('%s: Window is too long.',upper(mfilename));
end;

R=size(g,2);
% -------- Are we in the Riesz sequence of in the frame case

scale=1;
if a>M*R
  % Handle the Riesz basis (dual lattice) case.
  % Swap a and M, and scale differently.
  scale=a/M;
  tmp=a;
  a=M;
  M=tmp;
end;

% -------- Compute -------------

if kv.lt(2)==1
    % Rectangular case
    if (info.gl<=M) && (R==1)

        % Diagonal of the frame operator
        d = gabframediag(g,a,M,L);
        gd=g./long2fir(d,info.gl);

    else

        % Long window case

        % Just in case, otherwise the call is harmless.
        g=fir2long(g,L);

        gd=comp_gabdual_long(g,a,M)*scale;

    end;

else
    % Non-separable case
    g=fir2long(g,L);

    if (kv.nsalg==1) || (kv.nsalg==0 && kv.lt(2)<=2)

        mwin=comp_nonsepwin2multi(g,a,M,kv.lt,L);

        gdfull=comp_gabdual_long(mwin,a*kv.lt(2),M)*scale;

        % We need just the first vector
        gd=gdfull(:,1);

    else

        [s0,s1,br] = shearfind(L,a,M,kv.lt);

        if s1 ~= 0
            p1 = comp_pchirp(L,s1);
            g = p1.*g;
        end

        b=L/M;
        Mr = L/br;
        ar = a*b/br;

        if s0 == 0
            gd=comp_gabdual_long(g,ar,Mr);
        else
            p0=comp_pchirp(L,-s0);
            g = p0.*fft(g);
            gd=comp_gabdual_long(g,L/Mr,L/ar)*L;
            gd = ifft(conj(p0).*gd);
        end

        if s1 ~= 0
            gd = conj(p1).*gd;
        end

    end;

    if (info.gl<=M) && (R==1)
        gd=long2fir(gd,M);
    end;

end;

% --------- post process result -------

if isreal(g) && (kv.lt(2)==1 || kv.lt(2)==2)
  % If g is real and the lattice is either rectangular or quinqux, then
  % the output is known to be real.
  gd=real(gd);
end;

if info.wasrow
  gd=gd.';
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat