text:
classdef GrandlyUnifyView < handle
  properties
    model  % a ref for the model
    controller % a ref for the controller
      % this is necessary b/c we create buttons on-the-fly, and need
      % callbacks for them
    fig  % the figure handle
    projectFileLabelText
    projectFileText
    chooseProjectFileButton
    scoreFeatureLabelText
    scoreFeatureMatTexts=zeros(0,1)
    scoreFeatureJabTexts=zeros(0,1)
    chooseScoreFeatureJabButtons=zeros(0,1)
    classifierFileLabelText
    classifierFileText
    chooseClassifierFileButton
    gtExpDirsLabelText
    gtExpDirsListbox
    gtExpDirsAddButton
    gtExpDirsRemoveButton
    %quitButton
    convertButton
    spinDepth=0  % how many times spin() has been called, minus how many
                 % times unspin() has been called
  end
  methods
    function self=GrandlyUnifyView(model,controller)
      self.model=model;
      self.controller=controller;

      % Create the figure, make it invisible for now
      self.fig= ...
        figure('units','points', ...
               'resize','off', ...
               'numbertitle','off', ...
               'menubar','none', ...
               'name','Grandly Unify', ...
               'windowstyle','normal', ...
               'visible','off', ...
               'closeRequestFcn',@(g,e)(controller.quit()));

      % make the labels
      self.projectFileLabelText= ...
        uicontrol('parent',self.fig, ...
                  'units','points', ...
                  'style','text', ...
                  'backgroundcolor',get(self.fig,'color'), ...
                  'string','Project File ', ...
                  'HorizontalAlignment','right',...
                  'fontsize',12);

      self.classifierFileLabelText= ...
        uicontrol('parent',self.fig, ...
                  'units','points', ...
                  'style','text', ...
                  'backgroundcolor',get(self.fig,'color'), ...
                  'string','Classifier File ', ...
                  'HorizontalAlignment','right',...
                  'fontsize',12);

      self.gtExpDirsLabelText= ...
        uicontrol('parent',self.fig, ...
                  'units','points', ...
                  'style','text', ...
                  'backgroundcolor',get(self.fig,'color'), ...
                  'string',{'Ground-truth ', ...
                            'Experiment ', ...
                            'Directory Names '}, ...
                  'HorizontalAlignment','right',...
                  'fontsize',12);

      % Create the non-label controls
      self.projectFileText= ...
        uicontrol('parent',self.fig, ...
                  'style','text', ...
                  'backgroundcolor','w', ...
                  'HorizontalAlignment','left',...
                  'fontsize',12, ...
                  'units','points');
      self.classifierFileText= ...
        uicontrol('parent',self.fig, ...
                  'style','text', ...
                  'backgroundcolor','w', ...
                  'HorizontalAlignment','left',...
                  'fontsize',12, ...
                  'units','points');
      self.gtExpDirsListbox=...
        uicontrol('parent',self.fig, ...
                  'style','listbox', ...
                  'fontsize',12, ...
                  'units','points');
      self.chooseProjectFileButton= ...
        uicontrol('parent',self.fig, ...
                  'style','pushbutton', ...
                  'string','Choose', ...
                  'fontsize',12, ...
                  'units','points', ...
                  'callback',@(g,e)(controller.chooseProjectFileButtonPressed()));
      self.chooseClassifierFileButton= ...
        uicontrol('parent',self.fig, ...
                  'style','pushbutton', ...
                  'string','Choose', ...
                  'fontsize',12, ...
                  'units','points', ...
                  'callback',@(g,e)(controller.chooseClassifierFileButtonPressed()));
      self.gtExpDirsAddButton= ...
        uicontrol('parent',self.fig, ...
                  'style','pushbutton', ...
                  'string','Add', ...
                  'fontsize',12, ...
                  'units','points', ...
                  'callback',@(g,e)(controller.gtExpDirsAddButtonPressed()));
      self.gtExpDirsRemoveButton= ...
        uicontrol('parent',self.fig, ...
                  'style','pushbutton', ...
                  'string','Remove', ...
                  'fontsize',12, ...
                  'units','points', ...
                  'callback',@(g,e)(controller.gtExpDirsRemoveButtonPressed()));
      self.convertButton= ...
        uicontrol('parent',self.fig, ...
                  'style','pushbutton', ...
                  'string','Convert', ...
                  'fontsize',12, ...
                  'units','points', ...
                  'callback',@(g,e)(controller.convertButtonPressed()));

      % sync the view with the model
      self.update();

      % Get the screen dimensions
      savedUnits=get(0,'units');
      set(0,'units','points');
      screenPosition=get(0,'screensize');
      set(0,'units',savedUnits);

      % Figure the placement of the figure, want it centered relative to
      % primary screen
      pos=get(self.fig,'position');
      figureWidth=pos(3);
      figureHeight=pos(4);
      figureXOffset=screenPosition(1)+(screenPosition(3)-figureWidth )/2;
      figureYOffset=screenPosition(2)+(screenPosition(4)-figureHeight)/2;

      % Position the figure
      set(self.fig, ...
          'position',[figureXOffset figureYOffset figureWidth figureHeight]);

      % make it visible
      set(self.fig,'visible','on');
    end  % constructor method

    % ---------------------------------------------------------------------
    function spin(self)
      % Change the cursor to the watch cursor
      self.spinDepth=self.spinDepth+1;
      set(self.fig,'pointer','watch');
      drawnow('update');
    end

    % ---------------------------------------------------------------------
    function unspin(self)
      % Change the cursor back to normal, unless spin() has been called
      % multiple times without a corresponding unspin() calls.
      if self.spinDepth>0 ,
        self.spinDepth=self.spinDepth-1;
        drawnow('update');  % make sure the figure is as it should be before
                            % switching to the normal cursor
        set(self.fig,'pointer','arrow');
        drawnow('update');
      end
    end

    % ---------------------------------------------------------------------
    function update(self)
      % Update the view to match the model

      %
      % Do the layout
      %

      % Set figure layout parameters
      figureWidth=730+50;
      %figureHeightApprox=272;

      sideButtonWidth=80;  % choose buttons, add/remove buttons
      sideButtonHeight=26;
      sideButtonOffsetWidth=20;  % side buttons are this far to the right of
                                 % the center column
      smallSideButtonWidth=70;
      smallSideButtonHeight=24;
        % The 'Choose' buttons for the score feature jab file names are
        % small side buttons
      addRemoveInterButtonHeight=10;

      bottomButtonWidth=100;  % Convert and Quit are the bottom buttons
      bottomButtonHeight=30;
      %interBottomButtonWidth=20;

      centerColumnWidth=470+50;
      %fileNameTextHeight=20;
      gtExpDirsListboxHeight=100;

      labelOffsetWidth=20;  % Labels are this far to the left of the
                            % text/listbox that they label
      scoreFeatureJabIndentWidth=20;  % the names of the jab file names
                                      % are indented this much relative
                                      % to the score feature mat file names

      topPadHeight=22;
      underProjectPadHeight=10;
      interScoreFeaturePadHeight=5;
      belowScoreFeaturePadHeight=5;
      overClassifierPadHeight=10;
      classifierToGTPadHeight=21;
      bottomAreaHeight=70;
        % The "bottom area" is the rectangular area below the listbox.  It
        % is as wide as the figure, and it's bottom edge is aligned with
        % with bottom edge of the figure.  The row of bottom buttons is
        % centered within this area.

      % Want the top edge of the figure to stay fixed
      % The width of the figure is constant
      pos=get(self.fig,'position');
      figureTopYOffset=pos(2)+pos(4);
      figureXOffset=pos(1);

      % Get text extents
      extent=get(self.projectFileLabelText,'extent');
      extentSize=extent(3:4);
      projectFileLabelTextSize=extentSize;
      fileNameTextHeight=extentSize(2);  % height of a single line of text

      extent=get(self.classifierFileLabelText,'extent');
      extentSize=extent(3:4);
      %classifierFileLabelTextSize=extentSize+4;  % pad a little
      classifierFileLabelTextSize=extentSize;

      extent=get(self.gtExpDirsLabelText,'extent');
      extentSize=extent(3:4);
      %gtExpDirsLabelTextSize=extentSize+4;  % pad a little
      gtExpDirsLabelTextSize=extentSize;

      % Calculate overall figure height
      nScoreFeatures=length(self.model.scoreFeatureMatFileNames);
      scoreFeatureAreaHeight= ...
        2*nScoreFeatures*fileNameTextHeight + ...
        max(nScoreFeatures-1,0)*interScoreFeaturePadHeight + ...
        (nScoreFeatures>0)*belowScoreFeaturePadHeight;
      figureHeight=topPadHeight + ...
                   fileNameTextHeight + ...
                   underProjectPadHeight + ...
                   scoreFeatureAreaHeight + ...
                   overClassifierPadHeight + ...
                   fileNameTextHeight + ...
                   classifierToGTPadHeight + ...
                   gtExpDirsListboxHeight + ...
                   bottomAreaHeight;
      figureYOffset=figureTopYOffset-figureHeight;

      % Calculate positions
      projectCenterLineYOffset=figureHeight-topPadHeight-fileNameTextHeight/2;
      scoreFeatureAreaTopYOffset= ...
        projectCenterLineYOffset-fileNameTextHeight/2-underProjectPadHeight;
      classifierCenterLineYOffset= ...
        scoreFeatureAreaTopYOffset ...
          -scoreFeatureAreaHeight ...
          -overClassifierPadHeight ...
          -fileNameTextHeight/2;
      gtTopLineYOffset= ...
        classifierCenterLineYOffset ...
          -fileNameTextHeight/2 ...
          -classifierToGTPadHeight;

      centerColumnXOffset=(figureWidth-centerColumnWidth)/2;
      projectFileTextXOffset=centerColumnXOffset;
      classifierFileTextXOffset=centerColumnXOffset;
      gtExpDirsListboxXOffset=centerColumnXOffset;

      sideButtonXOffset=centerColumnXOffset+centerColumnWidth+sideButtonOffsetWidth;
      chooseProjectFileButtonXOffset=sideButtonXOffset;
      chooseClassifierFileButtonXOffset=sideButtonXOffset;
      gtExpDirsAddButtonXOffset=sideButtonXOffset;
      gtExpDirsRemoveButtonXOffset=sideButtonXOffset;

      convertButtonXOffset=(figureWidth-bottomButtonWidth)/2;
      projectFileLabelTextXOffset=centerColumnXOffset-labelOffsetWidth-projectFileLabelTextSize(1);
      classifierFileLabelTextXOffset=centerColumnXOffset-labelOffsetWidth-classifierFileLabelTextSize(1);
      gtExpDirsLabelTextXOffset=centerColumnXOffset-labelOffsetWidth-gtExpDirsLabelTextSize(1);

      projectFileLabelTextYOffset=projectCenterLineYOffset-projectFileLabelTextSize(2)/2;
      classifierFileLabelTextYOffset=classifierCenterLineYOffset-classifierFileLabelTextSize(2)/2;
      gtExpDirsLabelTextYOffset=gtTopLineYOffset-gtExpDirsLabelTextSize(2);

      projectFileTextYOffset=projectCenterLineYOffset-fileNameTextHeight/2;
      classifierFileTextYOffset=classifierCenterLineYOffset-fileNameTextHeight/2;
      gtExpDirsListboxYOffset=gtTopLineYOffset-gtExpDirsListboxHeight;

      chooseProjectFileButtonYOffset=projectCenterLineYOffset-sideButtonHeight/2;
      chooseClassifierFileButtonYOffset=classifierCenterLineYOffset-sideButtonHeight/2;
      gtExpDirsAddButtonYOffset=gtTopLineYOffset-sideButtonHeight;
      gtExpDirsRemoveButtonYOffset=gtTopLineYOffset-sideButtonHeight-addRemoveInterButtonHeight-sideButtonHeight;

      bottomButtonYOffset=(bottomAreaHeight-bottomButtonHeight)/2;
        % for both Quit and Convert buttons

      projectFileTextWidth=centerColumnWidth;
      classifierFileTextWidth=centerColumnWidth;
      gtExpDirsListboxWidth=centerColumnWidth;

      projectFileTextHeight=fileNameTextHeight;
      classifierFileTextHeight=fileNameTextHeight;

      % Position the figure
      set(self.fig, ...
          'position',[figureXOffset figureYOffset figureWidth figureHeight]);

      % Position the labels
      set(self.projectFileLabelText, ...
          'position',[projectFileLabelTextXOffset projectFileLabelTextYOffset projectFileLabelTextSize]);
      set(self.classifierFileLabelText, ...
          'position',[classifierFileLabelTextXOffset classifierFileLabelTextYOffset classifierFileLabelTextSize]);
      set(self.gtExpDirsLabelText, ...
          'position',[gtExpDirsLabelTextXOffset gtExpDirsLabelTextYOffset gtExpDirsLabelTextSize]);

      % Position everything else
      set(self.projectFileText, ...
          'position',[projectFileTextXOffset ...
                      projectFileTextYOffset ...
                      projectFileTextWidth ...
                      projectFileTextHeight]);
      set(self.classifierFileText , ...
          'position',[classifierFileTextXOffset ...
                      classifierFileTextYOffset ...
                      classifierFileTextWidth ...
                      classifierFileTextHeight]);
      set(self.gtExpDirsListbox, ...
          'position',[gtExpDirsListboxXOffset ...
                      gtExpDirsListboxYOffset ...
                      gtExpDirsListboxWidth ...
                      gtExpDirsListboxHeight]);
      set(self.chooseProjectFileButton, ...
          'position',[chooseProjectFileButtonXOffset ...
                      chooseProjectFileButtonYOffset ...
                      sideButtonWidth ...
                      sideButtonHeight]);
      set(self.chooseClassifierFileButton, ...
          'position',[chooseClassifierFileButtonXOffset ...
                      chooseClassifierFileButtonYOffset ...
                      sideButtonWidth ...
                      sideButtonHeight]);
      set(self.gtExpDirsAddButton, ...
          'position',[gtExpDirsAddButtonXOffset ...
                      gtExpDirsAddButtonYOffset ...
                      sideButtonWidth ...
                      sideButtonHeight]);
      set(self.gtExpDirsRemoveButton, ...
          'position',[gtExpDirsRemoveButtonXOffset ...
                      gtExpDirsRemoveButtonYOffset ...
                      sideButtonWidth ...
                      sideButtonHeight]);
      set(self.convertButton, ...
          'position',[convertButtonXOffset ...
                      bottomButtonYOffset ...
                      bottomButtonWidth ...
                      bottomButtonHeight]);

      %
      % Update the score feature controls
      %
      delete(self.scoreFeatureLabelText(ishandle(self.scoreFeatureLabelText)));
      delete(self.scoreFeatureMatTexts(ishandle(self.scoreFeatureMatTexts)));
      delete(self.scoreFeatureJabTexts(ishandle(self.scoreFeatureJabTexts)));
      delete(self.chooseScoreFeatureJabButtons(ishandle(self.chooseScoreFeatureJabButtons)));
      self.scoreFeatureLabelText=nan;
      self.scoreFeatureMatTexts=nan(nScoreFeatures,1);
      self.scoreFeatureJabTexts=nan(nScoreFeatures,1);
      self.chooseScoreFeatureJabButtons=nan(nScoreFeatures,1);
      if nScoreFeatures>0
        self.scoreFeatureLabelText= ...
          uicontrol('parent',self.fig, ...
                    'units','points', ...
                    'style','text', ...
                    'backgroundcolor',get(self.fig,'color'), ...
                    'string',{'Score Feature ', ...
                              'File Names '}, ...
                    'HorizontalAlignment','right',...
                    'fontsize',12, ...
                    'visible','off');
        extent=get(self.scoreFeatureLabelText,'extent');
        extentSize=extent(3:4);
        scoreFeatureLabelTextSize=extentSize;
        scoreFeatureLabelTextXOffset=centerColumnXOffset-labelOffsetWidth-scoreFeatureLabelTextSize(1);
        scoreFeatureLabelTextYOffset=scoreFeatureAreaTopYOffset-scoreFeatureLabelTextSize(2);
        set(self.scoreFeatureLabelText, ...
            'position',[scoreFeatureLabelTextXOffset scoreFeatureLabelTextYOffset scoreFeatureLabelTextSize], ...
            'visible','on');
      end
      for i=1:nScoreFeatures
        % Create the mat file name text object for this one
        xOffset=centerColumnXOffset;
        yOffset=scoreFeatureAreaTopYOffset-(i-1)*2*fileNameTextHeight-fileNameTextHeight;
        width=centerColumnWidth;
        height=fileNameTextHeight;
        self.scoreFeatureMatTexts(i)= ...
          uicontrol('parent',self.fig, ...
                    'units','points', ...
                    'style','text', ...
                    'backgroundcolor',get(self.fig,'color'), ...
                    'string',sprintf('%s:',self.model.scoreFeatureMatFileNames{i}), ...
                    'HorizontalAlignment','left',...
                    'fontsize',12, ...
                    'position',[xOffset yOffset width height]);
        % Create the jab file name text object for this one
        xOffset=centerColumnXOffset+scoreFeatureJabIndentWidth;
        yOffset=yOffset-fileNameTextHeight;
        width=centerColumnWidth-scoreFeatureJabIndentWidth;
        height=fileNameTextHeight;
        self.scoreFeatureJabTexts(i)= ...
          uicontrol('parent',self.fig, ...
                    'units','points', ...
                    'style','text', ...
                    'backgroundcolor',get(self.fig,'color'), ...
                    'HorizontalAlignment','right',...
                    'fontsize',12, ...
                    'position',[xOffset yOffset width height]);
        % Create the 'Choose' button for this one
        xOffset=centerColumnXOffset+centerColumnWidth+sideButtonOffsetWidth+(sideButtonWidth-smallSideButtonWidth)/2;
        yOffset=yOffset+fileNameTextHeight-sideButtonHeight/2+(sideButtonHeight-smallSideButtonHeight)/2;
        width=smallSideButtonWidth;
        height=smallSideButtonHeight;
        self.chooseScoreFeatureJabButtons(i)= ...
          uicontrol('parent',self.fig, ...
                    'style','pushbutton', ...
                    'string','Choose', ...
                    'fontsize',12, ...
                    'units','points', ...
                    'callback',@(g,e)(self.controller.chooseScoreFeatureJabFileButtonPressed(self.model.scoreFeatureMatFileNames{i})), ...
                    'position',[xOffset yOffset width height]);
      end

      %
      % Update the contents of UI elements
      %

      % Update the project file text object
      if self.model.isProjectFileSpecified() ,
        set(self.projectFileText, ...
            'String',self.model.projectFileName, ...
            'FontAngle','normal');
      else
        set(self.projectFileText, ...
            'String','Not specified', ...
            'FontAngle','italic');
      end

      % Update the score feature jab file text object(s)
      for i=1:nScoreFeatures
        if self.model.isScoreFeatureJabFileSpecified(i) ,
          set(self.scoreFeatureJabTexts(i), ...
              'String',[self.model.scoreFeatureJabFileNames{i} ' '], ...
              'FontAngle','normal');
        else
          set(self.scoreFeatureJabTexts(i), ...
              'String','Not specified ', ...
              'FontAngle','italic');
        end
      end

      % Update the classifier file text object
      if self.model.isClassifierFileSpecified() ,
        set(self.classifierFileText, ...
            'String',self.model.classifierFileName, ...
            'FontAngle','normal');
      else
        set(self.classifierFileText, ...
            'String','Not specified', ...
            'FontAngle','italic');
      end

      % Update the GT experiment listbox
      set(self.gtExpDirsListbox,'String',self.model.gtExpDirNames);
      set(self.gtExpDirsListbox,'Value',self.model.iCurrentGTExpDir);

      % Update enablement of Remove button
      set(self.gtExpDirsRemoveButton,'enable',offIff(isempty(self.model.gtExpDirNames)));

      % Update enablement of convert button
      set(self.convertButton,'enable',onIff(self.model.isProjectFileSpecified()));
    end  % method

    % ---------------------------------------------------------------------
    function close(self)
      if ishandle(self.fig)
        delete(self.fig);
        self.fig=[];
      end
    end  % method

  end
end

author: kristinbranson

repo: JAABA


text:
function r = invgamrand(a, b, varargin)
%INVGAMRAND Random matrices from inverse gamma distribution
%
%   R = INVGAMRAND(A,B) returns a matrix of random numbers chosen
%   from the inverse gamma distribution with parameters A and B.
%   The size of R is the common size of A and B if both are matrices.
%   If either parameter is a scalar, the size of R is the size of the other
%   parameter. Alternatively, R = GAMRND(A,B,M,N) returns an M by N matrix.
%
%   Note: Parameterization as in (Neal, 1996).
%      A is mean of the distribution
%      B is degrees of freedom
%
%	See also GAMRAND

% Copyright (c) 1999 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

error('No mex-file for this archtitecture. See Matlab help and convert.m in ./linuxCsource or ./winCsource for help.')

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function G=framematrix(F,L);
%FRAMEMATRIX  Frame synthesis operator matrix
%   Usage: G=framematrix(F,L);
%
%   `G=frsynmatrix(F,L)` returns the matrix representation *G* of the frame
%   synthesis operator for a frame *F* of length *L*. The frame object *F*
%   must have been created using |frame|.
%
%   The frame synthesis operator matrix contains all the frame atoms as
%   column vectors. It has dimensions $L \times Ncoef$, where $Ncoef$ is the
%   number of coefficients. The number of coefficients can be found as
%   `Ncoef=framered(F)*L`. This means that the frame matrix is usually
%   **very** large, and this routine should only be used for small values of
%   *L*.
%
%   The action of the frame analysis operator |frana| is equal to
%   multiplication with the Hermitean transpose of the frame
%   matrix. Consider the following simple example:::
%
%     L=200;
%     F=frame('dgt','gauss',10,20);
%     G=frsynmatrix(F,L);
%     testsig = randn(L,1);
%     res = frana(F,testsig)-G'*testsig;
%     norm(res)
%
%   See also: frame, frana, frsyn

warning(['LTFAT: FRAMEMATRIX has been deprecated and will be removed',...
         ' in the future releases, please use FRSYNMATRIX instead.']);

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

Lcheck=framelength(F,L);
if Lcheck~=L
    error('%s: Incompatible frame length.',upper(mfilename));
end;

if F.realinput

    %switch(F.type)
    %  case 'dgtreal'

    %  This code correctly reproduces the matrix represenation of the
    %  analysis operator, but not of the synthesis.
    %
    %    F2=frame('dgt',F.g,F.a,F.M);
    %    G2=frsynmatrix(F2,L);
    %    M2=floor(F.M/2)+1;
    %    N=L/F.a;
    %    G=zeros(L,M2*N);
    %    for n=0:N-1
    %        G(:,1+n*M2:(n+1)*M2)=G2(:,1+n*F.M:M2+n*F.M);
    %    end;

    %  otherwise
        error(['%s: The synthesis operator of real-valued-input frames does is ' ...
               'non-linear and does not have a matrix represenation.']);
        %end;
else

  % Generic code handles all frames where there are no extra coefficients
  % in the representation
  Ncoef = framered(F)*L;
  % sprintf for Octave compatibility
  assert(abs(Ncoef-round(Ncoef))<1e-3,sprintf('%s: There is a bug. Ncoef=%d should be an integer.',upper(mfilename),Ncoef));
  Ncoef=round(Ncoef);
  coef=eye(Ncoef);
  G = frsyn(F,coef);
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% By Roger Aarenstrup, roger.aarenstrup@mathworks.com
% 2006-08-18
%
% This is a state-space DC motor model of a
% Maxon RE25 10 Watt, precious metal brushes, 118743
%
% This model also have a weak connection to a load.
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%            The full dc motor model                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Vin is the input voltage to the motor
% i is the motor current
% th_m is the rotor angle, theta
% dth_m is the rotor angular velocity sometimes called omega
% th_l is the load angle
% dth_l is the load angular velocity

% Controller Sample Time
Ts = 100e-6;

% PARAMETERS DC MOTOR
Rm  = 2.06;             % Motor resistance (ohm)
Lm  = 0.000238;         % motor inductance (Henrys)
Kb = 1/((406*2*pi)/60); % Back EMF constant (Volt-sec/Rad)
Kt = 0.0235;            % Torque constand (Nm/A)
Jm = 1.07e-6;           % Rotor inertia (Kg m^2)
bm = 12e-7;             % MEchanical damping (linear model of
                        % friction: bm * dth)
% PARAMETERS LOAD
Jl = 10.07e-6;   % Load inertia (10 times the rotor)
bl = 12e-6;      % Load damping (friction)
Ks = 100;        % Spring constant for connection rotor/load
b = 0.0001;      % Spring damping for connection rotor/load

% SYSTEM MATRICES
%
% States:  [i dth_m th_m dth_l th_l]'
% Input:   Vin the motor voltage
% Outputs: same as states
%
Afull = [-Rm/Lm    -Kb/Lm       0       0       0;
         Kt/Jm   -(bm+b)/Jm   -Ks/Jm    b/Jm     Ks/Jm;
           0          1         0       0        0;
           0         b/Jl    Ks/Jl   -(b+bl)/Jl    -Ks/Jl;
           0          0         0       1        0];

Bfull = [1/Lm 0 0 0 0]';

Cfull = [0 1 0 0 0;
         0 0 1 0 0;
         0 0 0 1 0;
         0 0 0 0 1];

Dfull = [0 0 0 0]';

sys_full = ss(Afull, Bfull, Cfull, Dfull);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     The reduced dc motor model for current control      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% SYSTEM MATRICES
%
% States:  [dth_m th_m dth_l th_l]'
% Input:   I The current to the dc motor
% Outputs: same as states
%
Ared = [ -(bm+b)/Jm   -Ks/Jm    b/Jm     Ks/Jm;
              1         0       0        0;
             b/Jl    Ks/Jl   -(b+bl)/Jl    -Ks/Jl;
              0         0       1        0];

Bred = [Kt/Jm 0 0 0]';

Cred = eye(4);

Dred = [0 0 0 0]';

sys_red = ss(Ared, Bred, Cred, Dred);

% Discrete version of the model (as seen from the controller)
sys_red_d = c2d(sys_red, Ts, 'zoh');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        State-space controller design attempt 1          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The discrete poles of the closed loop system
dpole1 = exp(-2*pi*Ts*[20 22 24 26]);

% Calculate the control parameters
L1 = place(sys_red_d.a, sys_red_d.b, dpole1);

% Keep the static gain of the closed loop system to 1
F=feedback(sys_red_d, L1);   % The closed loop system, F.
Kstatic = freqresp(F, 0);    % Get the static gain of F.
Kstat = 1/Kstatic(4);        % It is the fourth output (load position)

% The above commands requires toolboxes, incase you don't have them
% you can simulate the system with a unit step is see the output static
% gain. Kstat will be the inverse of that.

% to verify the pole placement
%pzmap(F);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  State-space controller design attempt 2 - Integrator   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The discrete poles of the closed loop system
dpole2 = exp(-2*pi*Ts*[20 22 24 26 5]);

% We only consider one output here
Cred_one = [0 1 0 0];

% Calculate the control parameters
Aint = [sys_red_d.a [0 0 0 0]';
        -Cred_one 1];

Bint= [sys_red_d.b' 0]';

Cint = [Cred_one 0];

Dint = [0]';

sys_int_d2 = ss(Aint, Bint, Cint, Dint, Ts);

% Controllability VS Observability analysis
%ob = obsv(sys_int_d2.a, sys_int_d2.c);
%cr = ctrb(sys_int_d2.a, sys_int_d2.b);
%rank(ob)
%rank(cr)

[L2, a, m] = place(sys_int_d2.a, sys_int_d2.b, dpole2);

% Feed Forward gain
Kff = L2(5)/(dpole2(5) - 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  State-space controller design attempt 3 - Observer     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% These should be about twice as fast as the state feedback
% for the controller.
dpole3 = exp(-2*pi*Ts*[1400 1450 1500 1550]);
%dpole3 = [-0.3 -0.33 -0.35 -0.37];

% Recalculate the reduced model with one ouput
Cred2 = [0 1 0 0];

Dred2 = 0;

sys_red2 = ss(Ared, Bred, Cred2, Dred2);

% Discrete version of the model (as seen from the controller)
sys_red_d2 = c2d(sys_red2, Ts, 'zoh');

set(sys_red_d2, 'inputname', {'Um'}, ...
                'outputname', {'Enc_out'});

% Calculate the observer state feedback
K = place(sys_red_d2.a', sys_red_d2.c', dpole3);

% Put the observer together
Aobs = [sys_red_d2.a-K'*sys_red_d2.c];
Bobs = [sys_red_d2.b K'];
Cobs = eye(4);
Dobs = [0 0 0 0; 0 0 0 0]';

observer = ss(Aobs, Bobs, Cobs, Dobs, Ts, 'inputname', {'Ue', 'Enc_in'}, ...
              'outputname', {'dth_m', 'th_m', 'dth_l', 'th_l'});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  State-space controller design attempt 4 - The Servo Case  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Kinv = (Jl+Jm)/Kt;

dpole4 = exp(-2*pi*Ts*[200 204 220 224 300]);
[L3, a, m] = place(sys_int_d2.a, sys_int_d2.b, dpole4);

author: Sable

repo: mcbench-benchmarks


text:
function test_example_ica_eog

% MEM 8gb
% WALLTIME 00:30:00

%
%% Use independent component analysis (ICA) to remove EOG artifacts
%
%% # Description
%
% This script demonstrates how you can use ICA for cleaning the EOG artifacts from your MEG data. It consists of three steps:
%
% #  decomposition of the MEG data
% #  identifying the components that reflect eye artifacts
% #  removing those components and backprojecting the data
%
%% # Example dataset
%
% You can run the code below on your own data. Alternatively, try with the example MEG dataset _ArtifactMEG.ds_ (available from [ftp://ftp.fieldtriptoolbox.org/pub/fieldtrip/tutorial/ArtifactMEG.zip](ftp://ftp.fieldtriptoolbox.org/pub/fieldtrip/tutorial/ArtifactMEG.zip)). This dataset was acquired continuously with trials of 10 seconds. All figures in this example script are based on these data.
%
% To load this dataset into MATLAB and preprocess with FieldTrip, use:
%
% preprocessing of example dataset
cfg = [];
cfg.dataset            = dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/ArtifactMEG.ds');
cfg.trialdef.eventtype = 'trial';
cfg = ft_definetrial(cfg);

cfg.channel            = 'MEG';
cfg.continuous         = 'yes';
data = ft_preprocessing(cfg);

% downsample the data to speed up the next step
cfg = [];
cfg.resamplefs = 300;
cfg.detrend    = 'no';
data = ft_resampledata(cfg, data);

%% # ICA decomposition
%
% After reading in the preprocessed data into memory in FieldTrip format, you can continue with decomposing it in independent components.
%
% perform the independent component analysis (i.e., decompose the data)
cfg        = [];
cfg.method = 'runica'; % this is the default and uses the implementation from EEGLAB

comp = ft_componentanalysis(cfg, data);

% % Note that this is a time-consuming step. The output "comp" structure resembles the input raw data structure, i.e. it contains a time course for each component and each trial. Furthermore, it contains the spatial mixing matrix. In principle you can continue analyzing the data on the component level by doing
% %
%   cfg = [];
%   cfg = ...
%   freq = ft_freqanalysis(cfg, comp);
%
% % or
% %
%   cfg = [];
%   cfg = ...
%   timelock = ft_timelockanalysis(cfg, comp);

% but for this example we want to analyze the data eventually on the original channel level and only remove the components that represent the artifacts.
%
%% # Identify the artifacts
%
% plot the components for visual inspection
figure
cfg = [];
cfg.component = 1:20;       % specify the component(s) that should be plotted
cfg.layout    = 'CTF151.lay'; % specify the layout file that should be used for plotting
cfg.comment   = 'no';
ft_topoplotIC(cfg, comp)

% Make sure to plot and inspect all components. Write down the components that contain the eye artifacts. Very important is to know that on subsequent evaluations of the component decomposition result in components that **can have a different order**. That means that component numbers that you write down do not apply to another run of the ICA decomposition on the same data.
%
%
% The spatial topography of the components aids in interpreting whether a component represents activity from the cortex, or non-cortical physiological activity (muscle, eyes, heart) or even non-physiological activity (line noise and other environmental noise). If you are trained in this type of analysis, you can relatively easily spot the components that represent the eye movements: 9, 14 and 10.
%
% Besides the spatial topography you should inspect the time course of the components, which gives additional information on separating the cortical from the non-cortical contributions to the data.
%
% For further inspection of the time course of the components, use:
%
cfg = [];
cfg.layout = 'CTF151.lay'; % specify the layout file that should be used for plotting
cfg.viewmode = 'component';
ft_databrowser(cfg, comp)

% You can browse through the components and the trials. The EOG artifacts can be easily identified in the time course plots, see the figure below for an example.
%
%
%% # Remove the artifacts
%
% remove the bad components and backproject the data
cfg = [];
cfg.component = [9 10 14 24]; % to be removed component(s)
data = ft_rejectcomponent(cfg, comp, data);

% Compare the data before (red trace) and after (blue trace) the EOG removal - for example trial 4, channel MLF1
%


text:
function jed = epoch_to_jed_jelali ( )

%*****************************************************************************80
%
%% EPOCH_TO_JED_JELALI returns the epoch of the Jelali calendar as a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real JED, the Julian Ephemeris Date of the epoch.
%
  jed = 2114872.5;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [outexpdir] = ConvertMWTContour2Trx(experiment_name,rootoutputdir,varargin)

%% set parameters

% default parameters
rootdatadir = '/groups/branson/home/bransonk/behavioranalysis/code/Jdetect/data/larvae_mwt/rawdata';
contourext = 'outline';
trxfilestr = 'trx.mat';
perframedirstr = 'perframe';
MINDT = .001;
pxpermm = 11.111;
maxlength_mm = 5;
DEBUG = false;
makesoftlink = true;

[rootdatadir,contourext,trxfilestr,perframedirstr,DEBUG,maxlength_mm,pxpermm,makesoftlink] = ...
  myparse(varargin,...
  'rootdatadir',rootdatadir,...
  'contourext',contourext,...
  'trxfilestr',trxfilestr,...
  'perframedirstr',perframedirstr,...
  'debug',DEBUG,...
  'maxlength',maxlength_mm,...
  'pxpermm',pxpermm,...
  'makesoftlink',makesoftlink);

maxlength_px = ceil(maxlength_mm * pxpermm);

%% names of files

% find contour files for each trajectory
files = dir(fullfile(rootdatadir,[experiment_name,'.*.',contourext]));
contournames = {};
ids = [];
for i = 1:numel(files),
  match = regexp(files(i).name,['\.(\d+)\.',contourext,'$'],'tokens','once');
  if isempty(match),
    continue;
  end
  contournames{end+1} = fullfile(rootdatadir,files(i).name);
  ids(end+1) = str2double(match{1});
end

if isempty(contournames),
  error('Could not find any contours matching the pattern %s.*.%s',fullfile(rootdatadir,experiment_name),contourext);
end

ncontours = numel(contournames);

%% put all the contours into one trx file

% read everthing in from the contour files
alltimestamps = {};
trx = [];
arena = struct('x',nan,'y',nan,'r',nan);
for i = 1:ncontours,

  trk = struct('id',ids(i),'xcontour',{{}},'ycontour',{{}});
  timestamps = [];
  fid = fopen(contournames{i},'r');
  j = 1;
  while true,
    s = fgetl(fid);
    if ~ischar(s),
      break;
    end
    d = sscanf(s,'%f');
    if numel(d) < 3,
      warning('Contour line %s length < 3, skipping this line',s);
      continue;
    end
    if mod(numel(d),2) ~= 1,
      warning('Contour line %s length is not odd, skipping this line',s);
      continue;
    end
    d = d';
    timestamp = d(1);
    x = d(2:2:end-1);
    y = d(3:2:end);
    trk.xcontour{j} = x;
    trk.ycontour{j} = y;
    timestamps(j) = timestamp;
    j = j + 1;
  end
  fclose(fid);
  trk.arena = arena;
  alltimestamps{i} = timestamps;
  trx = structappend(trx,trk);
end

% make all the timestamps agree with each other
timestamps = union(alltimestamps{:});
dts = diff(timestamps);
% sanity check to look for rounding errors
if any(dts < MINDT),
  error('Rounding error check failed: there are timestamps < %f apart',MINDT);
end
t0s = cellfun(@(x) x(1), alltimestamps);
t1s = cellfun(@(x) x(end), alltimestamps);
fps = nanmean(dts);
for i = 1:ncontours,
  trx(i).firstframe = find(t0s(i) >= timestamps,1);
  trx(i).endframe = find(t1s(i) >= timestamps,1,'last');
  trx(i).off = 1 - trx(i).firstframe;
  % sanity check to make sure this matches the number of frames
  n = trx(i).endframe - trx(i).firstframe + 1;
  if numel(trx(i).xcontour) ~= n,
    error('Number of frames check failed: Number of frames of contour data does not match start and end timestamps');
  end
  trx(i).dt = dts(trx(i).firstframe:trx(i).endframe-1);
  trx(i).nframes = n;
  trx(i).fps = fps;
end

%% fit ellipses

[XGRID,YGRID] = meshgrid(0:maxlength_px-1,0:maxlength_px-1);

for i = 1:ncontours,

  trx(i).x = nan(1,trx(i).nframes);
  trx(i).y = nan(1,trx(i).nframes);
  trx(i).x_mm = nan(1,trx(i).nframes);
  trx(i).y_mm = nan(1,trx(i).nframes);
  trx(i).a = nan(1,trx(i).nframes);
  trx(i).b = nan(1,trx(i).nframes);
  trx(i).a_mm = nan(1,trx(i).nframes);
  trx(i).b_mm = nan(1,trx(i).nframes);
  trx(i).theta = nan(1,trx(i).nframes);
  trx(i).theta_mm = nan(1,trx(i).nframes);

  for t = 1:trx(i).nframes,

    % create an image
    xc = round(trx(i).xcontour{t}*pxpermm);
    yc = round(trx(i).ycontour{t}*pxpermm);
    minx = min(xc);
    maxx = max(xc);
    miny = min(yc);
    maxy = max(yc);
    nx = maxx-minx+1;
    ny = maxy-miny+1;
    bw = [inpolygon(XGRID(1:ny,1:nx),YGRID(1:ny,1:nx),[xc,xc(1)]-minx,[yc,yc(1)]-miny);false(maxlength_px-ny,nx)];

    % fit an ellipse
    x = XGRID(bw)+minx;
    y = YGRID(bw)+miny;
    mux = nanmean(x);
    muy = nanmean(y);
    S = cov([x,y],1);
    [a,b,theta] = cov2ell(S);
    % try and make theta be as continuous as possible
    if t > 1,
      dsame = abs(modrange(trx(i).theta(t-1)-theta,-pi,pi));
      dflip = abs(modrange(trx(i).theta(t-1)-theta+pi,-pi,pi));
      if dflip < dsame,
        theta = modrange(theta+pi,-pi,pi);
      end
    end

    % store quarter-major, quarter-minor axes
    a = a/2;
    b = b/2;
    mux_mm = mux / pxpermm;
    muy_mm = muy / pxpermm;
    a_mm = a / pxpermm;
    b_mm = b / pxpermm;
    trx(i).x(t) = mux_mm;
    trx(i).y(t) = muy_mm;
    trx(i).x_mm(t) = mux_mm;
    trx(i).y_mm(t) = muy_mm;
    trx(i).a(t) = a_mm;
    trx(i).b(t) = b_mm;
    trx(i).a_mm(t) = a_mm;
    trx(i).b_mm(t) = b_mm;
    trx(i).theta(t) = theta;
    trx(i).theta_mm(t) = theta;

    if DEBUG,
      hold off;
      plot(trx(i).xcontour{t},trx(i).ycontour{t},'k.-');
      hold on;
      drawellipse(trx(i).x(t),trx(i).y(t),trx(i).theta(t),trx(i).a(t)*2,trx(i).b(t)*2,'r');
      axis equal;
      title(num2str(t));
      drawnow;
    end

  end

end

%% plot the trajectories

if DEBUG,

nframes = numel(timestamps);
colors = jet(ncontours)*.7;
for t = 1:nframes,
  hold off;
  for i = 1:ncontours,
    if t > trx(i).endframe || t < trx(i).firstframe,
      continue;
    end
    j = t + trx(i).off;
    plot(trx(i).xcontour{j},trx(i).ycontour{j},'.-','color',colors(i,:));
    hold on;
    drawellipse(trx(i).x(j),trx(i).y(j),trx(i).theta(j),trx(i).a(j)*2,trx(i).b(j)*2,'color','r');
  end
  axis equal;
  drawnow;
end

end

%% create output directory

% make sure the root output directory exists
if ~exist(rootoutputdir,'dir'),
  mkdir(rootoutputdir);
end

% create the experiment directory
outexpdir = fullfile(rootoutputdir,experiment_name);
if ~exist(outexpdir,'dir'),
  mkdir(outexpdir);
end

%% create the trx file

outmatname = fullfile(outexpdir,trxfilestr);
save(outmatname,'trx');

%% copy over the contour files
for i = 1:ncontours,
  [~,basename] = myfileparts(contournames{i});
  outcontourname = fullfile(outexpdir,basename);
  if isunix && makesoftlink,
    cmd = sprintf('ln -s %s %s',contournames{i},outcontourname);
    unix(cmd);
  else
    [success,msg] = copyfile(contournames{i},outcontourname);
    if ~success,
      error('Error copying file %s to %s: %s',moviename,outmoviename,msg);
    end
  end
end

%% save the per-frame features
perframedir = fullfile(outexpdir,perframedirstr);
if ~exist(perframedir,'dir'),
  mkdir(perframedir);
end

% to do: save choreography's per-frame features


text:
function r = initialGuess(i_vars, stateGuess, algvarGuess)
ig = ocl.Variable.create(i_vars, 0);
ig.states.set(stateGuess);
ig.algvars.set(algvarGuess);
r = ig.value;

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
% Generate Figure 6 of the "Deep Scattering Spectrum" paper.

% Load the signal.
x = wavread('chord_signal.wav');

% Prepare the filters and scattering operators.
filt_opt.filter_type = {'gabor_1d','morlet_1d'};
filt_opt.B = [4 4];
filt_opt.Q = 4*filt_opt.B;
filt_opt.J = T_to_J([256 32768],filt_opt);

scat_opt.oversampling = 3;
scat_opt.M = 2;

Wop = wavelet_factory_1d(length(x), filt_opt, scat_opt);

% Compute scattering coefficients.
[S, U] = scat(x, Wop);

% Renormalize coefficients.
epsilon = 1e-3;
S = renorm_scat(S, epsilon);
S = renorm_1st(S, x, Wop{2}, epsilon);

S = resample_scat(S, 8, false);
U = resample_scat(U, 6, false);

% Compute the logarithm.
S = log_scat(S, 1e-3);
U = log_scat(U, 1e-3);

% Display scalogram U{2}, first-order coefficients S{2} and second-order
% coefficients S{3} for j1 = 58.
figure(6);
clf;
set(gcf,'Units','pixels','Position',[200 200 560 420]);
scattergram(U{2}, [], S{2}, [], S{3}, 58);
subplot(3,1,1);
set(gca,'XTick',[],'YTick',[]);
subplot(3,1,2);
set(gca,'XTick',[],'YTick',[]);
subplot(3,1,3);
set(gca,'XTick',[],'YTick',[]);
colormap(jet);

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [f,model] = penlab_callback_f(x,model)

global latest_x_f
global latest_df

x = x(:);
[f,latest_df] = fmincon_fun(x,model);
latest_x_f = x;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [stat, cfg] = statistics_wrapper(cfg, varargin)

% STATISTICS_WRAPPER performs the selection of the biological data for
% timelock, frequency or source data and sets up the design vector or
% matrix.
%
% The specific configuration options for selecting timelock, frequency
% of source data are described in FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS and
% FT_SOURCESTATISTICS, respectively.
%
% After selecting the data, control is passed on to a data-independent
% statistical subfunction that computes test statistics plus their associated
% significance probabilities and critical values under some null-hypothesis. The statistical
% subfunction that is called is FT_STATISTICS_xxx, where cfg.method='xxx'. At
% this moment, we have implemented two statistical subfunctions:
% FT_STATISTICS_ANALYTIC, which calculates analytic significance probabilities and critical
% values (exact or asymptotic), and FT_STATISTICS_MONTECARLO, which calculates
% Monte-Carlo approximations of the significance probabilities and critical values.
%
% The specific configuration options for the statistical test are
% described in FT_STATISTICS_xxx.

% This function depends on PREPARE_TIMEFREQ_DATA which has the following options:
% cfg.avgoverchan
% cfg.avgoverfreq
% cfg.avgovertime
% cfg.channel
% cfg.channelcmb
% cfg.datarepresentation (set in STATISTICS_WRAPPER cfg.datarepresentation = 'concatenated')
% cfg.frequency
% cfg.latency
% cfg.precision
% cfg.previous (set in STATISTICS_WRAPPER cfg.previous = [])
% cfg.version (id and name set in STATISTICS_WRAPPER)

% Copyright (C) 2005-2006, 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$

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed',     {'approach',   'method'});
cfg = ft_checkconfig(cfg, 'required',    {'method'});
cfg = ft_checkconfig(cfg, 'forbidden',   {'transform'});

% set the defaults
cfg.latency     = ft_getopt(cfg, 'latency', 'all');
cfg.frequency   = ft_getopt(cfg, 'frequency', 'all');
cfg.roi         = ft_getopt(cfg, 'roi',       []);
cfg.avgovertime = ft_getopt(cfg, 'avgovertime', 'no');
cfg.avgoverfreq = ft_getopt(cfg, 'avgoverfreq', 'no');
cfg.avgoverroi  = ft_getopt(cfg, 'avgoverroi',  'no');

% determine the type of the input and hence the output data
if ~exist('OCTAVE_VERSION', 'builtin')
  [s, i] = dbstack;
  if length(s)>1
    [caller_path, caller_name, caller_ext] = fileparts(s(2).name);
  else
    caller_path = '';
    caller_name = '';
    caller_ext  = '';
  end
  % evalin('caller', 'mfilename') does not work for Matlab 6.1 and 6.5
  istimelock = strcmp(caller_name,'ft_timelockstatistics');
  isfreq     = strcmp(caller_name,'ft_freqstatistics');
  issource   = strcmp(caller_name,'ft_sourcestatistics');
else
  % cannot determine the calling function in Octave, try looking at the
  % data instead
  istimelock  = isfield(varargin{1},'time') && ~isfield(varargin{1},'freq') && isfield(varargin{1},'avg');
  isfreq      = isfield(varargin{1},'time') && isfield(varargin{1},'freq');
  issource    = isfield(varargin{1},'pos') || isfield(varargin{1},'transform');
end

if (istimelock+isfreq+issource)~=1
  ft_error('Could not determine the type of the input data');
end

if istimelock || isfreq,
  % these defaults only apply to channel level data
  cfg.channel     = ft_getopt(cfg, 'channel',     'all');
  cfg.avgoverchan = ft_getopt(cfg, 'avgoverchan', 'no');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% collect the biological data (the dependent parameter)
%  and the experimental design (the independent parameter)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if issource

  % test that all source inputs have the same dimensions and are spatially aligned
  for i=2:length(varargin)
    if isfield(varargin{1}, 'dim') && (length(varargin{i}.dim)~=length(varargin{1}.dim) || ~all(varargin{i}.dim==varargin{1}.dim))
      ft_error('dimensions of the source reconstructions do not match, use FT_VOLUMENORMALISE first');
    end
    if isfield(varargin{1}, 'pos') && (length(varargin{i}.pos(:))~=length(varargin{1}.pos(:)) || ~all(varargin{i}.pos(:)==varargin{1}.pos(:)))
      ft_error('grid locations of the source reconstructions do not match, use FT_VOLUMENORMALISE first');
    end
    if isfield(varargin{1}, 'transform') && ~all(varargin{i}.transform(:)==varargin{1}.transform(:))
      ft_error('spatial coordinates of the source reconstructions do not match, use FT_VOLUMENORMALISE first');
    end
  end

  Nsource = length(varargin);
  Nvoxel  = length(varargin{1}.inside) + length(varargin{1}.outside);

  if ~isempty(cfg.roi)
    if ischar(cfg.roi)
      cfg.roi = {cfg.roi};
    end
    % the source representation should specify the position of each voxel in MNI coordinates
    x = varargin{1}.pos(:,1);  % this is from left (negative) to right (positive)
    % determine the mask to restrict the subsequent analysis
    % process each of the ROIs, and optionally also left and/or right separately
    roimask  = {};
    roilabel = {};
    for i=1:length(cfg.roi)
      if islogical(cfg.roi{i})
        tmp = cfg.roi{i};
      else
        tmpcfg.roi = cfg.roi{i};
        tmpcfg.inputcoord = cfg.inputcoord;
        tmpcfg.atlas = cfg.atlas;
        tmp = ft_volumelookup(tmpcfg, varargin{1});
      end
      if strcmp(cfg.avgoverroi, 'no') && ~isfield(cfg, 'hemisphere')
        % no reason to deal with separated left/right hemispheres
        cfg.hemisphere = 'combined';
      end

      if     strcmp(cfg.hemisphere, 'left')
        tmp(x>=0)    = 0;  % exclude the right hemisphere
        roimask{end+1}  = tmp;
        roilabel{end+1} = ['Left '  cfg.roi{i}];

      elseif strcmp(cfg.hemisphere, 'right')
        tmp(x<=0)    = 0;  % exclude the right hemisphere
        roimask{end+1}  = tmp;
        roilabel{end+1} = ['Right ' cfg.roi{i}];

      elseif strcmp(cfg.hemisphere, 'both')
        % deal separately with the voxels on the left and right side of the brain
        tmpL = tmp; tmpL(x>=0) = 0;  % exclude the right hemisphere
        tmpR = tmp; tmpR(x<=0) = 0;  % exclude the left hemisphere
        roimask{end+1}  = tmpL;
        roimask{end+1}  = tmpR;
        roilabel{end+1} = ['Left '  cfg.roi{i}];
        roilabel{end+1} = ['Right ' cfg.roi{i}];
        clear tmpL tmpR

      elseif strcmp(cfg.hemisphere, 'combined')
        % all voxels of the ROI can be combined
        roimask{end+1}  = tmp;
        if ischar(cfg.roi{i})
          roilabel{end+1} = cfg.roi{i};
        else
          roilabel{end+1} = ['ROI ' num2str(i)];
        end

      else
        ft_error('incorrect specification of cfg.hemisphere');
      end
      clear tmp
    end % for each roi

    % note that avgoverroi=yes is implemented differently at a later stage
    % avgoverroi=no is implemented using the inside/outside mask
    if strcmp(cfg.avgoverroi, 'no')
      for i=2:length(roimask)
        % combine them all in the first mask
        roimask{1} = roimask{1} | roimask{i};
      end
      roimask = roimask{1};  % only keep the combined mask
      % the source representation should have an inside and outside vector containing indices
      sel = find(~roimask);
      varargin{1}.inside  = setdiff(varargin{1}.inside, sel);
      varargin{1}.outside = union(varargin{1}.outside, sel);
      clear roimask roilabel
    end % if avgoverroi=no
  end

  % get the source parameter on which the statistic should be evaluated
  if strcmp(cfg.parameter, 'mom') && isfield(varargin{1}, 'avg') && isfield(varargin{1}.avg, 'csdlabel') && isfield(varargin{1}, 'cumtapcnt')
    [dat, cfg] = get_source_pcc_mom(cfg, varargin{:});
  elseif strcmp(cfg.parameter, 'mom') && isfield(varargin{1}, 'avg') && ~isfield(varargin{1}.avg, 'csdlabel')
    [dat, cfg] = get_source_lcmv_mom(cfg, varargin{:});
  elseif isfield(varargin{1}, 'trial')
    [dat, cfg] = get_source_trial(cfg, varargin{:});
  else
    [dat, cfg] = get_source_avg(cfg, varargin{:});
  end
  cfg.dimord = 'voxel';

  % note that avgoverroi=no is implemented differently at an earlier stage
  if strcmp(cfg.avgoverroi, 'yes')
    tmp = zeros(length(roimask), size(dat,2));
    for i=1:length(roimask)
      % the data only reflects those points that are inside the brain,
      % the atlas-based mask reflects points inside and outside the brain
      roi = roimask{i}(varargin{1}.inside);
      tmp(i,:) = mean(dat(roi,:), 1);
    end
    % replace the original data with the average over each ROI
    dat = tmp;
    clear tmp roi roimask
    % remember the ROIs
    cfg.dimord = 'roi';
  end

  % check whether the original input data contains a dim, which would allow
  % for 3D reshaping and clustering with bwlabeln
  if isfield(varargin{1}, 'transform') || ((isfield(varargin{1}, 'dim') && prod(varargin{1}.dim)==size(varargin{1}.pos,1)))
    cfg.connectivity = 'bwlabeln';
  end

elseif isfreq || istimelock
  % get the ERF/TFR data by means of PREPARE_TIMEFREQ_DATA
  cfg.datarepresentation = 'concatenated';
  [cfg, data] = prepare_timefreq_data(cfg, varargin{:});
  cfg = rmfield(cfg, 'datarepresentation');

  dim = size(data.biol);
  if length(dim)<3
    % seems to be singleton frequency and time dimension
    dim(3)=1;
    dim(4)=1;
  elseif length(dim)<4
    % seems to be singleton time dimension
    dim(4)=1;
  end
  cfg.dimord = 'chan_freq_time';

  % the dimension of the original data (excluding the replication dimension) has to be known for clustering
  cfg.dim = dim(2:end);
  % all dimensions have to be concatenated except the replication dimension and the data has to be transposed
  dat = transpose(reshape(data.biol, dim(1), prod(dim(2:end))));

  % remove to save some memory
  data.biol = [];

  % verify that neighbours are present in case needed (not needed when
  % averaging over channels)
  if ~(isfield(cfg, 'avgoverchan') && istrue(cfg.avgoverchan)) &&...
      ~isfield(cfg,'neighbours') && isfield(cfg, 'correctm') && strcmp(cfg.correctm, 'cluster')
    ft_error('You need to specify a neighbourstructure');
    %cfg.neighbours = ft_neighbourselection(cfg,varargin{1});
  end

end

% get the design from the information in cfg and data.
if ~isfield(cfg,'design')
  ft_warning('Please think about how you would create cfg.design.  Soon the call to prepare_design will be deprecated')
  cfg.design = data.design;
  [cfg] = prepare_design(cfg);
end

if size(cfg.design,2)~=size(dat,2)
  cfg.design = transpose(cfg.design);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute the statistic, using the data-independent statistical subfunction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% determine the function handle to the intermediate-level statistics function
if exist(['ft_statistics_' cfg.method], 'file')
  statmethod = str2func(['ft_statistics_' cfg.method]);
else
  ft_error('could not find the corresponding function for cfg.method="%s"\n', cfg.method);
end
fprintf('using "%s" for the statistical testing\n', func2str(statmethod));

% check that the design completely describes the data
if size(dat,2) ~= size(cfg.design,2)
  ft_error('the size of the design matrix (%d) does not match the number of observations in the data (%d)', size(cfg.design,2), size(dat,2));
end

% determine the number of output arguments
try
  % the nargout function in Matlab 6.5 and older does not work on function handles
  num = nargout(statmethod);
catch
  num = 1;
end

design = cfg.design;
cfg    = rmfield(cfg,'design'); % to not confuse lower level functions with both cfg.design and design input

% perform the statistical test
if strcmp(func2str(statmethod),'ft_statistics_montecarlo') % because ft_statistics_montecarlo (or to be precise, clusterstat) requires to know whether it is getting source data,
                                                        % the following (ugly) work around is necessary
  if num>1
    [stat, cfg] = statmethod(cfg, dat, design);
  else
    [stat] = statmethod(cfg, dat, design);
  end
else
  if num>1
    [stat, cfg] = statmethod(cfg, dat, design);
  else
    [stat] = statmethod(cfg, dat, design);
  end
end

if isstruct(stat)
  % the statistical output contains multiple elements, e.g. F-value, beta-weights and probability
  statfield = fieldnames(stat);
else
  % only the probability was returned as a single matrix, reformat into a structure
  dum = stat; stat = []; % this prevents a Matlab warning that appears from release 7.0.4 onwards
  stat.prob = dum;
  statfield = fieldnames(stat);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add descriptive information to the output and rehape into the input format
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if issource
   if ~isfield(varargin{1},'dim') % FIX ME; this is added temporarily (20110427) to cope with ft_sourceanalysis output not having a dim field since r3273
        varargin{1}.dim = [Nvoxel 1];
   end
  if isempty(cfg.roi) || strcmp(cfg.avgoverroi, 'no')
    % remember the definition of the volume, assume that they are identical for all input arguments
    try stat.dim       = varargin{1}.dim;        end
    try stat.xgrid     = varargin{1}.xgrid;      end
    try stat.ygrid     = varargin{1}.ygrid;      end
    try stat.zgrid     = varargin{1}.zgrid;      end
    try stat.inside    = varargin{1}.inside;     end
    try stat.outside   = varargin{1}.outside;    end
    try stat.pos       = varargin{1}.pos;        end
    try stat.transform = varargin{1}.transform;  end
    try stat.freq      = varargin{1}.freq;       end
    try stat.time      = varargin{1}.time;       end
  else
    stat.inside  = 1:length(roilabel);
    stat.outside = [];
    stat.label   = roilabel(:);
  end
  for i=1:length(statfield)
    tmp = getsubfield(stat, statfield{i});
    if isfield(varargin{1}, 'inside') && numel(tmp)==length(varargin{1}.inside)
      % the statistic was only computed on voxels that are inside the brain
      % sort the inside and outside voxels back into their original place
      if islogical(tmp)
        tmp(varargin{1}.inside)  = tmp;
        tmp(varargin{1}.outside) = false;
      else
        tmp(varargin{1}.inside)  = tmp;
        tmp(varargin{1}.outside) = nan;
      end
      extradim = 1;
    elseif isfield(varargin{1}, 'inside') && isfield(varargin{1}, 'freq') && isfield(varargin{1}, 'time') && numel(tmp)==length(varargin{1}.inside)*length(varargin{1}.freq)*length(varargin{1}.time)
      % the statistic is a higher dimensional matrix (here as a function of freq) computed only on the
      % inside voxels
      newtmp = zeros(length(varargin{1}.inside)+length(varargin{1}.outside), length(varargin{1}.freq), length(varargin{1}.time));
      if islogical(tmp)
        newtmp(varargin{1}.inside, :, :)  = reshape(tmp, length(varargin{1}.inside), length(varargin{1}.freq), []);
        newtmp(varargin{1}.outside, :, :) = false;
      else
        newtmp(varargin{1}.inside, :, :)  = reshape(tmp, length(varargin{1}.inside), length(varargin{1}.freq), []);
        newtmp(varargin{1}.outside, :, :) = nan;
      end
      tmp = newtmp; clear newtmp;
      extradim = length(varargin{1}.freq);
    elseif isfield(varargin{1}, 'inside') && isfield(varargin{1}, 'freq') && numel(tmp)==length(varargin{1}.inside)*length(varargin{1}.freq)
      % the statistic is a higher dimensional matrix (here as a function of freq) computed only on the
      % inside voxels
      newtmp = zeros(length(varargin{1}.inside)+length(varargin{1}.outside), length(varargin{1}.freq));
      if islogical(tmp)
        newtmp(varargin{1}.inside, :)  = reshape(tmp, length(varargin{1}.inside), []);
        newtmp(varargin{1}.outside, :) = false;
      else
        newtmp(varargin{1}.inside, :)  = reshape(tmp, length(varargin{1}.inside), []);
        newtmp(varargin{1}.outside, :) = nan;
      end
      tmp = newtmp; clear newtmp;
      extradim = length(varargin{1}.freq);
    elseif isfield(varargin{1}, 'inside') && isfield(varargin{1}, 'time') && numel(tmp)==length(varargin{1}.inside)*length(varargin{1}.time)
      % the statistic is a higher dimensional matrix (here as a function of time) computed only on the
      % inside voxels
      newtmp = zeros(length(varargin{1}.inside)+length(varargin{1}.outside), length(varargin{1}.time));
      if islogical(tmp)
        newtmp(varargin{1}.inside, :)  = reshape(tmp, length(varargin{1}.inside), []);
        newtmp(varargin{1}.outside, :) = false;
      else
        newtmp(varargin{1}.inside, :)  = reshape(tmp, length(varargin{1}.inside), []);
        newtmp(varargin{1}.outside, :) = nan;
      end
      tmp = newtmp; clear newtmp;
      extradim = length(varargin{1}.time);
    else
      extradim = 1;
    end
    if numel(tmp)==prod(varargin{1}.dim)
      % reshape the statistical volumes into the original format
      stat = setsubfield(stat, statfield{i}, reshape(tmp, varargin{1}.dim));
    else
      stat = setsubfield(stat, statfield{i}, tmp);
    end
  end
else
  haschan    = isfield(data, 'label');    % this one remains relevant, even after averaging over channels
  haschancmb = isfield(data, 'labelcmb'); % this one remains relevant, even after averaging over channels
  hasfreq = strcmp(cfg.avgoverfreq, 'no') && ~any(isnan(data.freq));
  hastime = strcmp(cfg.avgovertime, 'no') && ~any(isnan(data.time));
  stat.dimord = '';

  if haschan
    stat.dimord = [stat.dimord 'chan_'];
    stat.label  = data.label;
    chandim = dim(2);
  elseif haschancmb
    stat.dimord   = [stat.dimord 'chancmb_'];
    stat.labelcmb = data.labelcmb;
    chandim = dim(2);
  end

  if hasfreq
    stat.dimord = [stat.dimord 'freq_'];
    stat.freq   = data.freq;
    freqdim = dim(3);
  end

  if hastime
    stat.dimord = [stat.dimord 'time_'];
    stat.time   = data.time;
    timedim = dim(4);
  end

  if ~isempty(stat.dimord)
    % remove the last '_'
    stat.dimord = stat.dimord(1:(end-1));
  end

  for i=1:length(statfield)
    try
      % reshape the fields that have the same dimension as the input data
      if     strcmp(stat.dimord, 'chan')
        stat = setfield(stat, statfield{i}, reshape(getfield(stat, statfield{i}), [chandim 1]));
      elseif strcmp(stat.dimord, 'chan_time')
        stat = setfield(stat, statfield{i}, reshape(getfield(stat, statfield{i}), [chandim timedim]));
      elseif strcmp(stat.dimord, 'chan_freq')
        stat = setfield(stat, statfield{i}, reshape(getfield(stat, statfield{i}), [chandim freqdim]));
      elseif strcmp(stat.dimord, 'chan_freq_time')
        stat = setfield(stat, statfield{i}, reshape(getfield(stat, statfield{i}), [chandim freqdim timedim]));
      elseif strcmp(stat.dimord, 'chancmb_time')
        stat = setfield(stat, statfield{i}, reshape(getfield(stat, statfield{i}), [chandim timedim]));
      elseif strcmp(stat.dimord, 'chancmb_freq')
        stat = setfield(stat, statfield{i}, reshape(getfield(stat, statfield{i}), [chandim freqdim]));
      elseif strcmp(stat.dimord, 'chancmb_freq_time')
        stat = setfield(stat, statfield{i}, reshape(getfield(stat, statfield{i}), [chandim freqdim timedim]));
      elseif strcmp(stat.dimord, 'time')
        stat = setfield(stat, statfield{i}, reshape(getfield(stat, statfield{i}), [1 timedim]));
      elseif strcmp(stat.dimord, 'freq')
        stat = setfield(stat, statfield{i}, reshape(getfield(stat, statfield{i}), [1 freqdim]));
      elseif strcmp(stat.dimord, 'freq_time')
        stat = setfield(stat, statfield{i}, reshape(getfield(stat, statfield{i}), [freqdim timedim]));
      end
    end
  end
end

return % statistics_wrapper main()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION for extracting the data of interest
% data resembles PCC beamed source reconstruction, multiple trials are coded in mom
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dat, cfg] = get_source_pcc_mom(cfg, varargin)
Nsource = length(varargin);
Nvoxel  = length(varargin{1}.inside) + length(varargin{1}.outside);
Ninside = length(varargin{1}.inside);
dim     = varargin{1}.dim;
for i=1:Nsource
  ntrltap = sum(varargin{i}.cumtapcnt);
  dat{i}  = zeros(Ninside, ntrltap);
  for j=1:Ninside
    k = varargin{1}.inside(j);
    if ischar(varargin{i}.avg.csdlabel{1})
      % this represents the 'old-style' csdlabel, one list for all dipoles
      dipsel  = find(strcmp(varargin{i}.avg.csdlabel, 'scandip'));
    else
      % this represents the 'new-style' csdlabel, one list for each of the dipoles
      dipsel  = find(strcmp(varargin{i}.avg.csdlabel{k}, 'scandip'));
    end
    dat{i}(j,:) = reshape(varargin{i}.avg.mom{k}(dipsel,:), 1, ntrltap);
  end
end
% concatenate the data matrices of the individual input arguments
dat = cat(2, dat{:});
% remember the dimension of the source data
cfg.dim = dim;
% remember which voxels are inside the brain
cfg.inside = varargin{1}.inside;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION for extracting the data of interest
% data resemples LCMV beamed source reconstruction, mom contains timecourse
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dat, cfg] = get_source_lcmv_mom(cfg, varargin)
Nsource = length(varargin);
Nvoxel  = length(varargin{1}.inside) + length(varargin{1}.outside);
Ntime   = length(varargin{1}.avg.mom{varargin{1}.inside(1)});
Ninside = length(varargin{1}.inside);
dim     = [varargin{1}.dim Ntime];
dat     = zeros(Ninside*Ntime, Nsource);
for i=1:Nsource
  % collect the 4D data of this input argument
  tmp = nan(Ninside, Ntime);
  for j=1:Ninside
    k = varargin{1}.inside(j);
    tmp(j,:) = reshape(varargin{i}.avg.mom{k}, 1, dim(4));
  end
  dat(:,i) = tmp(:);
end
% remember the dimension of the source data
cfg.dim = dim;
% remember which voxels are inside the brain
cfg.inside = varargin{1}.inside;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION for extracting the data of interest
% data contains single-trial or single-subject source reconstructions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dat, cfg] = get_source_trial(cfg, varargin)
Nsource = length(varargin);
Nvoxel  = length(varargin{1}.inside) + length(varargin{1}.outside);

for i=1:Nsource
  Ntrial(i) = length(varargin{i}.trial);
end
k = 1;
for i=1:Nsource
  for j=1:Ntrial(i)
    tmp = getsubfield(varargin{i}.trial(j), cfg.parameter);
    if ~iscell(tmp),
      %dim = size(tmp);
      dim = [Nvoxel 1];
    else
      dim = [Nvoxel size(tmp{varargin{i}.inside(1)})];
    end
    if i==1 && j==1 && numel(tmp)~=Nvoxel,
      ft_warning('the input-data contains more entries than the number of voxels in the volume, the data will be concatenated');
      dat    = zeros(prod(dim), sum(Ntrial)); %FIXME this is old code should be removed
    elseif i==1 && j==1 && iscell(tmp),
      ft_warning('the input-data contains more entries than the number of voxels in the volume, the data will be concatenated');
      dat    = zeros(Nvoxel*numel(tmp{varargin{i}.inside(1)}), sum(Ntrial));
    elseif i==1 && j==1,
      dat = zeros(Nvoxel, sum(Ntrial));
    end
    if ~iscell(tmp),
      dat(:,k) = tmp(:);
    else
      Ninside   = length(varargin{i}.inside);
      %tmpvec    = (varargin{i}.inside-1)*prod(size(tmp{varargin{i}.inside(1)}));
      tmpvec    = varargin{i}.inside;
      insidevec = [];
      for m = 1:numel(tmp{varargin{i}.inside(1)})
        insidevec = [insidevec; tmpvec(:)+(m-1)*Nvoxel];
      end
      insidevec = insidevec(:)';
      tmpdat  = reshape(permute(cat(3,tmp{varargin{1}.inside}), [3 1 2]), [Ninside numel(tmp{varargin{1}.inside(1)})]);
      dat(insidevec, k) = tmpdat(:);
    end
    % add original dimensionality of the data to the configuration, is required for clustering
    %FIXME: this was obviously wrong, because often trial-data is one-dimensional, so no dimensional information is present
    %cfg.dim = dim(2:end);
    k = k+1;
  end
end
if isfield(varargin{1}, 'inside')
  fprintf('only selecting voxels inside the brain for statistics (%.1f%%)\n', 100*length(varargin{1}.inside)/prod(dim));
  for j=prod(dim(2:end)):-1:1
    dat((j-1).*dim(1) + varargin{1}.outside, :) = [];
  end
end
% remember the dimension of the source data
if isfield(varargin{1}, 'dim')
  cfg.dim = varargin{1}.dim;
else
  cfg.dim = size(varargin{1}.trial(1).(cfg.parameter));
  cfg.dim(1) = numel(varargin{1}.inside);
  cfg.insideorig = varargin{1}.inside;
  cfg.inside     = varargin{1}.inside;
end

%if ~isfield(cfg, 'dim')
  %warning('for clustering on trial-based data you explicitly have to specify cfg.dim');
%end
% remember which voxels are inside the brain
cfg.inside = varargin{1}.inside;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION for extracting the data of interest
% get the average source reconstructions, the repetitions are in multiple input arguments
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dat, cfg] = get_source_avg(cfg, varargin)
Nsource = length(varargin);
Nvoxel  = length(varargin{1}.inside) + length(varargin{1}.outside);
if isfield(varargin{1}, 'dim')
  dim    = varargin{1}.dim;
  inside = false(prod(dim),1);
else
  dim    = numel(varargin{1}.inside);
  inside = false(size(varargin{1}.pos,1),1);
end
inside(varargin{1}.inside) = true;

tmp     = getsubfield(varargin{1}, cfg.parameter);
if size(tmp,2)>1
  if isfield(varargin{1}, 'freq')
    Nvoxel = Nvoxel*numel(varargin{1}.freq);
    dim    = [dim numel(varargin{1}.freq)];
    inside = repmat(inside, [1 numel(varargin{1}.freq)]);
  end
  if isfield(varargin{1}, 'time')
    Nvoxel = Nvoxel*numel(varargin{1}.time);
    dim    = [dim numel(varargin{1}.time)];
    if isfield(varargin{1},'freq')
      inside = repmat(inside, [1 1 numel(varargin{1}.time)]);
    else
      inside = repmat(inside, [1 numel(varargin{1}.time)]);
    end
  end
end
dat = zeros(Nvoxel, Nsource);
for i=1:Nsource
  tmp = getsubfield(varargin{i}, cfg.parameter);
  dat(:,i) = tmp(:);
end
insideorig = find(inside(:,1));
inside = find(inside(:));
if isfield(varargin{1}, 'inside')
  if isfield(varargin{1}, 'pos')
    npos = size(varargin{1}.pos,1);
  elseif isfield(varargin{1}, 'dim')
    npos = prod(varargin{1}.dim);
  else
    npos = length(varargin{1}.inside); % uninformative, at least prevents crash
  end
  fprintf('only selecting voxels inside the brain for statistics (%.1f%%)\n', 100*length(varargin{1}.inside)/npos);
  dat = dat(inside,:);
end
% remember the dimension of the source data
cfg.dim = dim;
% remember which voxels are inside the brain
cfg.inside     = inside(:);
cfg.insideorig = insideorig;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION for creating a design matrix
% should be called in the code above, or in prepare_design
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [cfg] = get_source_design_pcc_mom(cfg, varargin)
% should be implemented

function [cfg] = get_source_design_lcmv_mom(cfg, varargin)
% should be implemented

function [cfg] = get_source_design_trial(cfg, varargin)
% should be implemented

function [cfg] = get_source_design_avg(cfg, varargin)
% should be implemented

author: spm


text:
function [Y,Error] = Sammon(X,d)
% Sammon mapping on X
% Usage: [Y,Error] = sammon(X,D)

% References:
% Sammon, John W. Jr., "A Nolinear Mapping for Data Structure Analysis",
% IEEE Transactions on Computers, vol. C-18, no. 5, pp 401-409, May 1969.

    %% Default setting
    maxHalves = 20;
    maxIter   = 500;
    theta     = 1e-9;

    %% Obtain distance matrix
    Dist = pdist2(X,X);

    %% Remaining inialisation
    N     = size(X,1);
    scale = 0.5/sum(Dist(:));
    Dist  = Dist + eye(N);
    Dinv  = 1./Dist;

    %% Initialized by PCA, map to d dimension
    [UU,DD] = svd(X);
    Y       = UU(:,1:d)*DD(1:d,1:d);
    %% Randomly initialize Y
    % Y      = randn(N,d);

    % Obtain error
    oneMatrix = ones(N,d);
    dist      = pdist2(Y,Y) + eye(N);
    dinv      = 1./dist;
    delta     = Dist - dist;
    Error     = sum(sum((delta.^2).*Dinv));

    %% Optimize process
    for i = 1 : maxIter
        % Compute gradient, Hessian and search direction
        delta    = dinv - Dinv;
        deltaOne = delta*oneMatrix;
        g        = delta*Y - Y.*deltaOne;
        dinv3    = dinv.^3;
        Y2       = Y.^2;
        H        = dinv3*Y2 - deltaOne - 2*Y.*(dinv3*Y) + Y2.*(dinv3*oneMatrix);
        s        = -g(:)./abs(H(:));
        yOld     = Y;

        % Use step-halving procedure to ensure progress id made
        for j = 1 : maxHalves
            Y(:)   = yOld(:) + s;
            dist   = pdist2(Y,Y) + eye(N);
            dinv   = 1./dist;
            delta  = Dist - dist;
            newError = sum(sum((delta.^2).*Dinv));
            if newError < Error
                break;
            else
                s = 0.5*s;
            end
        end

        % Evaluate termination criterion
        if abs((Error-newError)/Error) < theta
            break;
        end

        % Update error
        Error = newError;
    end
    Error = Error * scale;
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function RoPSs = RoPSFunc(mesh,neighborSize,binSize, rotaSize, LocalFrames)

%  Author: Yulan Guo {yulan.guo@nudt.edu.cn}
%  NUDT, China & CSSE, UWA, Australia
%
% This function takes a mesh and local reference frames (LRFs) for a set of
% keypoints as an input, and generate RoPS feature descriptors for the
% keypoints as an output. An illustration is shown in Fig4 in Ref.[1].
%
% Arguments : mesh - with vertices and faces
%                       neighborSize - the size of neighborhood to define a
%                                                  local surface for a selected keypoint
%                       binSize - number of bins for 2D plane partition
%                       rotaSize - number of rotations around each coordinate axis
%                       LocalFrames - local reference frames corresponding to all
%                                                keypoints of the input mesh

% Return :         RoPS - RoPS feature descriptors corresponding to all
%                                     keypoints
%
% Copyright : This code is written by Yulan Guo {yulan.guo@nudt.edu.cn}, NUDT.
%               The code may be used, modified and distributed for research purposes with
%              acknowledgement of the author and inclusion this copyright information.
%References:
%           [1] Yulan Guo, Ferdous Sohel, Mohammed Bennamoun, Min Lu, Jianwei Wan.
%           Rotational Projection Statistics for 3D Local Surface Description and Object Recognition.
%           Internation Journal of Computer Vision. 2013, 105 (1), 63-86
%           [2] Yulan Guo, Mohammed Bennamoun, Ferdous A Sohel, Min Lu, Jianwei Wan.
%           3D Object Recognition in Cluttered Scenes with Local Surface Features: A Survey.
%           IEEE Transactions on Pattern Analysis and Machine Intelligence,2014
%
% Disclaimer: This code is provided as is without any warranty.

%%%%%%%%%%%% parallel computing to accelerate.
p = gcp('nocreate'); % If no pool, do not create new one.
if isempty(p)   % if exists no parallel computing handle, create one.
    poolsize = parpool;
    % disp('Create a Parallel Computing Pool to accelarate process' );
else
    poolsize = p.NumWorkers;
    % disp(sprintf( 'Parpool has %d workers', poolsize) );
end

keypntIdx = mesh.keypntIdx;
% BucketSize = floor(length(mesh.vertices)/100);
% kdtreeVertices = KDTreeSearcher(mesh.vertices,'Distance','euclidean','BucketSize',BucketSize);
% BucketSize = floor(length(mesh.vertices)/100);
kdtreeVertices = KDTreeSearcher(mesh.vertices,'Distance','euclidean');

interval = pi/2/rotaSize;
parfor i = 1:length(keypntIdx)
    %obtain the neighboring point of a keypoint
    keypnt = mesh.vertices(keypntIdx(i),:);
    [neighborIdx2,neighborDis]= rangesearch(kdtreeVertices,keypnt,neighborSize);
    neighborIdx2 = cell2mat(neighborIdx2);
    neighborIdx = neighborIdx2(2:end);
    neighbNum = length(neighborIdx);
    if neighbNum<=1
        RoPSs{i,1} = ones(rotaSize*45,1)/sum(ones(rotaSize*45,1));
        continue;
    end

    %transfom the neighboring point to the local reference frame (LRF)
    rotation = LocalFrames{i};
    neighbor = [];
    for j=1:neighbNum
        neighbor(j,:) = (mesh.vertices(neighborIdx(j),:)-mesh.vertices(keypntIdx(i),:))*inv(rotation);
    end

    RoPS = [];
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %calculate the sub-feature of the keypoint along the Z axis
    for rotaIdx = 1:rotaSize
        rotaAngle = (rotaIdx-1)*interval + interval/2;
        R = [cos(rotaAngle) sin(rotaAngle) 0; -sin(rotaAngle) cos(rotaAngle) 0; 0 0 1]';
        rotaNeighbor = neighbor*R;
        %projection on the XY plane
        projNeighborXY = [rotaNeighbor(:,1),rotaNeighbor(:,2)];
        histTemp = subRoPSFunc(projNeighborXY,binSize);
        RoPS = [RoPS,histTemp];
        %projection on the XZ plane
        projNeighborXZ = [rotaNeighbor(:,1),rotaNeighbor(:,3)];
        histTemp = subRoPSFunc(projNeighborXZ,binSize);
        RoPS = [RoPS,histTemp];
        %projection on the YZ plane
        projNeighborYZ = [rotaNeighbor(:,2),rotaNeighbor(:,3)];
        histTemp = subRoPSFunc(projNeighborYZ,binSize);
        RoPS = [RoPS,histTemp];
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %calculate the sub-feature of the keypoint along the Y axis
    for rotaIdx = 1:rotaSize
        rotaAngle = (rotaIdx-1)*interval + interval/2;
        R = [cos(rotaAngle) 0 sin(rotaAngle); 0 1 0 ;-sin(rotaAngle) 0 cos(rotaAngle)]';
        rotaNeighbor = neighbor*R;
        %projection on the XY plane
        projNeighborXY = [rotaNeighbor(:,1),rotaNeighbor(:,2)];
        histTemp = subRoPSFunc(projNeighborXY,binSize);
        RoPS = [RoPS,histTemp];
        %projection on the XZ plane
        projNeighborXZ = [rotaNeighbor(:,1),rotaNeighbor(:,3)];
        histTemp = subRoPSFunc(projNeighborXZ,binSize);
        RoPS = [RoPS,histTemp];
        %projection on the YZ plane
        projNeighborYZ = [rotaNeighbor(:,2),rotaNeighbor(:,3)];
        histTemp = subRoPSFunc(projNeighborYZ,binSize);
        RoPS = [RoPS,histTemp];
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %calculate the sub-feature of the keypoint along the X axis
    for rotaIdx = 1:rotaSize
        rotaAngle = (rotaIdx-1)*interval + interval/2;
        R = [1,0,0; 0,cos(rotaAngle),sin(rotaAngle); 0, -sin(rotaAngle), cos(rotaAngle)]';
        rotaNeighbor = neighbor*R;
        %projection on the XY plane
        projNeighborXY = [rotaNeighbor(:,1),rotaNeighbor(:,2)];
        histTemp = subRoPSFunc(projNeighborXY,binSize);
        RoPS = [RoPS,histTemp];
        %projection on the XZ plane
        projNeighborXZ = [rotaNeighbor(:,1),rotaNeighbor(:,3)];
        histTemp = subRoPSFunc(projNeighborXZ,binSize);
        RoPS = [RoPS,histTemp];
        %projection on the YZ plane
        projNeighborYZ = [rotaNeighbor(:,2),rotaNeighbor(:,3)];
        histTemp = subRoPSFunc(projNeighborYZ,binSize);
        RoPS = [RoPS,histTemp];
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    RoPSs{i,1} = (RoPS)'/sum(RoPS);
end

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR


text:
function [ fx, fy ] = f06_f1 ( n, x, y )

%*****************************************************************************80
%
%% F06_F1 returns first derivatives of function 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N,1), Y(N,1), the evalution points.
%
%    Output, real FX(N,1), FY(N,1), the derivative values.
%
  fx = zeros ( n, 1 );
  fy = zeros ( n, 1 );

  t4(1:n,1) = 64.0 - 81.0 * ( ( x(1:n,1) - 0.5 ).^2 + ( y(1:n,1) - 0.5 ).^2 );

  i = ( 0.0 < t4(1:n,1) );

  fx(i,1) = - 9.0 * ( x(i,1) - 0.5 ) ./ sqrt ( t4(i,1) );
  fy(i,1) = - 9.0 * ( y(i,1) - 0.5 ) ./ sqrt ( t4(i,1) );

  return
end