text:
function [ALLEEG cfg handles] = pop_vis_TimeFreqGrid(ALLEEG,typeproc,varargin)
%
% Create a Time-Frequency Grid from a connectivity matrix. For details on
% the Interactive Time-Frequency Grid see [1].
% This function generates an interactive GUI. Default GUI options can be
% pre-specified as <name, value> pairs as listed in vis_TimeFreqGrid().
% Alternately, an arg configuration structure can be provided in lieu of
% name, value pairs (see Example below).
%
% Inputs:
%
%       ALLEEG:     Array of 1 or 2 EEGLAB datasets. If 2 sets, this function
%                   plots difference (and Conn must also be 1x2 array of
%                   structs)
%       Conn:       SIFT Connectivity Structure
%
% Optional:
%
%       <Name, Value> pairs as defined in vis_TimeFreqGrid()
%
% Outputs:
%
%       cfg:                          Argument specification struct. Can be
%                                     supplied in lieu of <name, value>
%                                     argument pairs to exactly reconstruct
%                                     TF-grid.
%       figureHandles:                Handles to figures.
%
% Example:  >> cfg = pop_vis_TimeFreqGrid(ALLEEG);    % now select some options from GUI and click "OK"
%           >> vis_TimeFreqGrid(ALLEEG,cfg);          % TF-grid will now be recreated from command-line
%                                                     % exactly as specified in GUI
%
% See Also: vis_TimeFreqGrid()
%
% References:
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual. Chapter 6.
%   Available at: http://www.sccn.ucsd.edu/wiki/Sift
%
% Author: Tim Mullen, 2010, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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 3 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.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

handles = [];

if length(varargin) == 1 && isempty(varargin{1})
    varargin = {};
end

if nargin<2
    typeproc = 0;
end

fcnName     = strrep(mfilename,'pop_','');
fcnHandle   = str2func(fcnName);

% check the dataset
res = hlp_checkeegset(ALLEEG,{'conn'});
if ~isempty(res)
    error(['SIFT:' fcnName],res{1});
end

% check if we have statistics
hasStats = isempty(hlp_checkeegset(ALLEEG,{'stats'}));

% pass in stats as separate option
if hasStats, varargin = [varargin 'Stats',ALLEEG(1).CAT.Stats]; end

if isfield(ALLEEG(1).CAT.configs,fcnName)
    % get default configuration (from prior use) and merge with varargin
    varargin = [hlp_struct2varargin(ALLEEG(1).CAT.configs.(fcnName)) varargin];
end

for cnd=1:length(ALLEEG);
    % make node labels a row array
    ALLEEG(cnd).CAT.curComponentNames = ALLEEG(cnd).CAT.curComponentNames(:)';
end

% extract the Connectivity structures
for cnd=1:length(ALLEEG);
    Conn(cnd) = ALLEEG(cnd).CAT.Conn;
end

if strcmpi(typeproc,'nogui')
    % get the config from function
    cfg = arg_tovals(arg_report('rich',fcnHandle,[{'EEG',ALLEEG,'Conn',Conn},varargin]),false);
else
    % render the GUI
    [PGh figh] = feval(['gui_' fcnName],ALLEEG,Conn,varargin{:});

    if isempty(PGh)
        % user chose to cancel
        cfg = [];
        return;
    end

    % get the specification of the PropertyGrid
    ps = PGh.GetPropertySpecification;
    cfg = arg_tovals(ps,false);
end

drawnow;

if strcmpi(typeproc,'cfg_only')
    return;
end

% execute the low-level function
if length(ALLEEG)==2 && ~cfg.plotCondDiff.arg_selection
    % there is more than one condition and user wants to
    % plot each condition sequentially
    for cnd=1:length(ALLEEG)
        if hasStats
            handles{cnd} = vis_TimeFreqGrid('ALLEEG',ALLEEG(cnd),'Conn',Conn(cnd),'Stats',ALLEEG(cnd).CAT.Stats,cfg);
        else
            handles{cnd} = vis_TimeFreqGrid('ALLEEG',ALLEEG(cnd),'Conn',Conn(cnd),cfg);
        end
        if ~isempty(cfg)
            % store the configuration structure
            ALLEEG(cnd).CAT.configs.(fcnName) = cfg;
        end
    end
else
    % either user wants to plot condition difference
    % or there is only one condition
    if hasStats
            handles = vis_TimeFreqGrid('ALLEEG',ALLEEG,'Conn',Conn,'Stats',ALLEEG(1).CAT.Stats,cfg);
        else
            handles = vis_TimeFreqGrid('ALLEEG',ALLEEG,'Conn',Conn,cfg);
    end
    if ~isempty(cfg) && length(ALLEEG)==1
        for k=1:length(ALLEEG)
            % store the configuration structure
            ALLEEG(k).CAT.configs.(fcnName) = cfg;
        end
    end
end

function x = hasStatsField(EEG)
    x = isfield(EEG.CAT,'Stats');

repo: REST


text:
% std_findsameica() - find groups of datasets with identical ICA decomposiotions
%
% Usage:
%        >> clusters = std_findsameica(ALLEEG);
%        >> clusters = std_findsameica(ALLEEG,icathreshold);
% Inputs:
%   ALLEEG           - a vector of loaded EEG dataset structures of all sets in the STUDY set.
%   icathreshold     - Threshold to compare icaweights
%
% Outputs:
%   cluster - cell array of groups of datasets
%
% Authors:  Arnaud Delorme, SCCN, INC, UCSD, July 2009-

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, arno@sccn.ucsd.edu
%
% 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.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% Coding notes: Useful information on functions and global variables used.

function cluster = std_findsameica(ALLEEG, varargin);

% 6/2/2014 Ramon : Allow ica threshold as input.
if nargin == 1
    icathreshold = 2e-5;
elseif nargin == 2;
    icathreshold = varargin{1};
end

cluster = { [1] };
for index = 2:length(ALLEEG)

    found = 0;
    for c = 1:length(cluster)
        if all(size(ALLEEG(cluster{c}(1)).icaweights) == size(ALLEEG(index).icaweights))
            %if isequal(ALLEEG(cluster{c}(1)).icaweights, ALLEEG(index).icaweights)
            if sum(sum(abs(ALLEEG(cluster{c}(1)).icaweights-ALLEEG(index).icaweights))) < icathreshold
                cluster{c}(end+1) = index;
                found = 1;
                break;
            end;
        end;
    end;
    if ~found
        cluster{end+1} = index;
    end;
end;


text:
function [ lat, long ] = xy2latlong( x, y, ref_lat, ref_long )
% converts cartesian coordinates x,y (in km) to lat/long
% ref_lat and ref_long is the geodetic reference point for plane approximation of the earth surface

earth_circumf = 40074;

lat = (y * 360 / earth_circumf) + ref_lat;
long = ( (x*360) / (earth_circumf * cos(ref_lat*pi/180)) ) + ref_long;

end

repo: tdoa-evaluation-rtlsdr


text:
function X_cache=fmridesign(frametimes,slicetimes, ...
                            events,S,exclude,hrf_parameters)

%FMRIDESIGN
%
% Produces a set of design matrices, one for each slice, for fmristat.
% With just the frametimes, it gives the hemodynamic response function.
%
% X_CACHE = FMRIDESIGN( FRAME_TIMES [, SLICE_TIMES [, EVENTS , [S ,
%                       [, EXCLUDE [, HRF_PARAMETERS ]]]]] )
%
% FRAME_TIMES is a row vector of frame acquisition times in seconds.
%
% SLICE_TIMES is a row vector of relative slice acquisition times,
% i.e. absolute acquisition time of a slice is FRAME_TIMES + SLICE_TIMES.
% Default is 0.
%
% EVENTS is a matrix whose rows are events and whose columns are:
% 1. id - an integer from 1:(number of events) to identify event type;
% 2. times - start of event, synchronised with frame and slice times;
% 3. durations (optional - default is 0) - duration of event;
% 4. heights (optional - default is 1) - height of response for event.
% For each event type, the response is a box function starting at the event
% times, with the specified durations and heights, convolved with the
% hemodynamic response function (see below). If the duration is zero, the
% response is the hemodynamic response function whose integral is
% the specified height - useful for `instantaneous' stimuli such as visual
% stimuli. The response is then subsampled at the appropriate frame and slice
% times to create a design matrix for each slice, whose columns correspond
% to the event id number. EVENT_TIMES=[] will ignore event times and just
% use the stimulus design matrix S (see next). Default is [1 0].
%
% S: Events can also be supplied by a stimulus design matrix,
% whose rows are the frames, and column are the event types. Events
% are created for each column, beginning at the frame time for each row
% of S, with a duration equal to the time to the next frame, and a height
% equal to the value of S for that row and column. Note that a
% constant term is not usually required, since it is removed by the
% polynomial trend terms provided N_POLY>=0. Note that all values for
% all frames must be supplied, because smoothing and lagging by the
% hemodynamic resonse is done BEFORE excluding time points by EXCLUDE.
% Default is [].
%
% EXCLUDE is a list of frames that should be excluded from the
% analysis. This must be used with Siemens EPI scans to remove the
% first few frames, which do not represent steady-state images.
% Default is [].
%
% HRF_PARAMETERS is a matrix whose rows are 6 parameters for the
% hemodynamic response function, one row for each event type and column
% of S (if there is just one row, this is repeated as necessary).
% The hrf is modeled as the difference of two
% gamma density functions (Glover, NeuroImage, 9:416-429).
% The components of HRF_PARAMETERS are:
% 1. PEAK1: time to the peak of the first gamma density;
% 2. FWHM1: approximate FWHM of the first gamma density;
% 3. PEAK2: time to the peak of the second gamma density;
% 4. FWHM2: approximate FWHM of the second gamma density;
% 5. DIP: coefficient of the second gamma density;
%    Final hrf is:   gamma1/max(gamma1)-DIP*gamma2/max(gamma2)
%    scaled so that its total integral is 1.
% 6. FIT_SCALE: 1 - fit the time scale of the hrf by convolving its
%    derivative with the specified column of the design matrix, to create an
%    additional column for the design matrix. Dividing the effect
%    of this column by the effect of the hrf itself estimates the
%    scale shift. 0 ignores this option.
% If PEAK1=0 then there is no smoothing of that event type with the hrf.
% Default is: [5.4 5.2 10.8 7.35 0.35 0] chosen by Glover (1999) for
% an auditory stimulus.
%
% X_CACHE: A cache of the design matrices; rows are the non-excluded frames,
% columns are all the regressor variables, with slices running slowest.

%############################################################################
% COPYRIGHT:   Copyright 2000 K.J. Worsley and C. Liao,
%              Department of Mathematics and Statistics,
%              McConnell Brain Imaging Center,
%              Montreal Neurological Institute,
%              McGill University, Montreal, Quebec, Canada.
%              worsley@math.mcgill.ca, liao@math.mcgill.ca
%
%              Permission to use, copy, modify, and distribute this
%              software and its documentation for any purpose and without
%              fee is hereby granted, provided that the above copyright
%              notice appear in all copies.  The author and McGill University
%              make no representations about the suitability of this
%              software for any purpose.  It is provided "as is" without
%              express or implied warranty.
%############################################################################

% Defaults:

if nargin < 2
   slicetimes=0
end
if nargin < 3
   events=[1 0]
end
if nargin < 4
   S=[]
end
if nargin < 5
   exclude=[]
end
if nargin < 6
   hrf_parameters=[5.4 5.2 10.8 7.35 0.35 0]
end

numframes=length(frametimes);
numslices=length(slicetimes);

% Keep time points that are not excluded:

allpts = 1:numframes;
allpts(exclude) = zeros(1,length(exclude));
keep = allpts( find( allpts ) );
n=length(keep);
scantimes=frametimes(keep);

if ~isempty(events)
   numevents=size(events,1);
   eventid=events(:,1);
   numeventypes=max(eventid);
   eventime=events(:,2);
   if size(events,2)>=3
      duration=events(:,3);
   else
      duration=zeros(numevents,1);
   end
   if size(events,2)>=4
      height=events(:,4);
   else
      height=ones(numevents,1);
   end
   mineventime=min(eventime);
   maxeventime=max(eventime+duration);
else
   numeventypes=0;
   mineventime=Inf;
   maxeventime=-Inf;
end

if ~isempty(S)
   numcolS=size(S,2);
else
   numcolS=0;
end

% Set up response matrix:

dt=0.02;
  startime=min(mineventime,min(frametimes)+min([slicetimes 0]));
finishtime=max(maxeventime,max(frametimes)+max([slicetimes 0]));
numtimes=ceil((finishtime-startime)/dt)+1;
numresponses=numeventypes+numcolS;
response=zeros(numtimes,numresponses);

if ~isempty(events)
   height=height./(1+(duration==0)*(dt-1));
   for k=1:numevents
      type=eventid(k);
      n1=ceil((eventime(k)-startime)/dt)+1;
      n2=ceil((eventime(k)+duration(k)-startime)/dt)+(duration(k)==0);
      if n2>=n1
         response(n1:n2,type)=response(n1:n2,type)+height(k)*ones(n2-n1+1,1);
      end
   end
end

if ~isempty(S)
   for j=1:numcolS
      for i=find(S(:,j)')
         n1=ceil((frametimes(i)-startime)/dt)+1;
         if i<numframes
            n2=ceil((frametimes(i+1)-startime)/dt);
         else
            n2=numtimes;
         end
         if n2>=n1
            response(n1:n2,numeventypes+j)= ...
               response(n1:n2,numeventypes+j)+S(i,j)*ones(n2-n1+1,1);
         end
      end
   end
end

% Set hrf parameters:

numscale=0;
for k=1:numresponses
   if k<=size(hrf_parameters,1)
      if hrf_parameters(k,1)>0
         peak1=hrf_parameters(k,1);
         fwhm1=hrf_parameters(k,2);
         peak2=hrf_parameters(k,3);
         fwhm2=hrf_parameters(k,4);
         dip=hrf_parameters(k,5);
         alpha1=peak1^2/fwhm1^2*8*log(2);
         alpha2=peak2^2/fwhm2^2*8*log(2);
         beta1=fwhm1^2/peak1/8/log(2);
         beta2=fwhm2^2/peak2/8/log(2);

         numlags=ceil(max(peak1+2*fwhm1,peak2+2*fwhm2)/dt);
         time=(0:(numlags-1))'*dt;
         gamma1=(time/peak1).^alpha1.*exp(-(time-peak1)./beta1);
         gamma2=(time/peak2).^alpha2.*exp(-(time-peak2)./beta2);
         hrf=gamma1-dip*gamma2;
         sumhrf=sum(hrf);
         hrf=hrf/sumhrf;
         if hrf_parameters(k,6)==1
            fit_scale=1;
            d_hrf=((time/beta1-alpha1-1).*gamma1- ...
               dip*(time/beta2-alpha2-1).*gamma2);
            d_hrf=d_hrf/sumhrf;
         else
            fit_scale=0;
         end
      else
         fit_scale=0;
         hrf=1;
      end
   end
   eventmatrix(:,k)=conv2(response(:,k),hrf);
   if fit_scale==1
      numscale=numscale+1;
      eventmatrix(:,numresponses+numscale)=conv2(response(:,k),d_hrf);
   end
end

% Make all the design matrices for each slice:

numcolX=numresponses+numscale;
X_cache=zeros(n,numcolX*numslices);
for slice = 1:numslices
   subtime=floor((scantimes+slicetimes(slice)-startime)/dt)+1;
   X_cache(:,(1:numcolX)+(slice-1)*numcolX)=eventmatrix(subtime,:);
end

% End.

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function tP = subsasgn(tP,s,b)
% overloads subsasgn

if ~isa(tP,'triplePointList'), tP = b([]); end

switch s(1).type

  case '()'

    if numel(s)>1, b =  subsasgn(subsref(tP,s(1)),s(2:end),b); end

    s(1).subs = {subsind(tP,s(1).subs)};

    if isempty(b)

      s(1).subs = {~s(1).subs{1}};
      tP = subsref(tP,s(1));

    else

      tP = subsasgn@dynProp(tP,s(1),b);
      tP.misrotation = subsasgn(tP.misrotation,s(1),b.misrotation);

      id1 = subsasgn(tP.grainId(:,1),s(1),b.grainId(:,1));
      id2 = subsasgn(tP.grainId(:,2),s(1),b.grainId(:,2));
      tP.grainId = [id1(:),id2(:)];

      id1 = subsasgn(tP.phaseId(:,1),s(1),b.phaseId(:,1));
      id2 = subsasgn(tP.phaseId(:,2),s(1),b.phaseId(:,2));
      tP.phaseId = [id1(:),id2(:)];

      id1 = subsasgn(tP.F(:,1),s(1),b.F(:,1));
      id2 = subsasgn(tP.F(:,2),s(1),b.F(:,2));
      tP.F = [id1(:),id2(:)];

      tP.V = b.V;
      tP.CSList = b.CSList;
      tP.phaseMap = b.phaseMap;

    end

  otherwise

    if tP.isProperty(s(1).subs) % otherwise a property

      tP = subsasgn@dynProp(tP,s,b);

    else

      tP = builtin('subsasgn',tP,s,b);

    end

end
end

author: mtex-toolbox

repo: mtex

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


text:
function u = fem2d_bvp_quadratic ( nx, ny, a, c, f, x, y )

%*****************************************************************************80
%
%% FEM2D_BVP_QUADRATIC solves boundary value problem on a rectangle.
%
%  Discussion:
%
%    The program uses the finite element method, with piecewise quadratic basis
%    functions to solve a 2D boundary value problem over a rectangle.
%
%    The following differential equation is imposed inside the region:
%
%      - d/dx a(x,y) du/dx - d/dy a(x,y) du/dy + c(x,y) * u(x,y) = f(x,y)
%
%    where a(x,y), c(x,y), and f(x,y) are given functions.
%
%    On the boundary, the solution is constrained to have the value 0.
%
%    The finite element method will use a regular grid of NX nodes in X, and
%    NY nodes in Y.  Both NX and NY must be odd.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NX, NY, the number of X and Y grid values.
%    NX and NY must be odd and at least 3.
%
%    Input, function A(X,Y), evaluates a(x,y);
%
%    Input, function C(X,Y), evaluates c(x,y);
%
%    Input, function F(X,Y), evaluates f(x,y);
%
%    Input, real X(NX), Y(NY), the mesh points.
%
%    Output, real U(NX,NY), the finite element coefficients, which are also
%    the value of the computed solution at the mesh points.
%

%
%  Quadrature definitions.
%
  quad_num = 3;
  abscissa(1) = -0.774596669241483377035853079956;
  abscissa(2) = 0.000000000000000000000000000000;
  abscissa(3) = 0.774596669241483377035853079956;
  weight(1) = 0.555555555555555555555555555556;
  weight(2) = 0.888888888888888888888888888889;
  weight(3) = 0.555555555555555555555555555556;
%
%  Make room for the matrix A and right hand side b.
%
  mn = nx * ny;
  A = zeros ( mn, mn );
  b = zeros ( mn, 1 );
%
%  Compute the matrix entries by integrating over each element.
%
  ex_num = ( nx - 1 ) / 2;
  ey_num = ( ny - 1 ) / 2;

  for ex = 1 : ex_num

    w = 2 * ex - 1;
    cc = 2 * ex;
    e = 2 * ex + 1;

    xx(1) = x(w);
    xx(2) = x(cc);
    xx(3) = x(e);

    for ey = 1 : ey_num

      s = 2 * ey - 1;
      mm = 2 * ey;
      n = 2 * ey + 1;

      yy(1) = y(s);
      yy(2) = y(mm);
      yy(3) = y(n);
%
%  Node indices
%
%  7  8  9   wn cn en
%  4  5  6   wm cm em
%  1  2  3   ws cs es
%
      node(1) = ( 2 * ey - 2 ) * nx + ( ex - 1 ) * 2 + 1;
      node(2) = ( 2 * ey - 2 ) * nx + ( ex - 1 ) * 2 + 2;
      node(3) = ( 2 * ey - 2 ) * nx + ( ex - 1 ) * 2 + 3;
      node(4) = ( 2 * ey - 1 ) * nx + ( ex - 1 ) * 2 + 1;
      node(5) = ( 2 * ey - 1 ) * nx + ( ex - 1 ) * 2 + 2;
      node(6) = ( 2 * ey - 1 ) * nx + ( ex - 1 ) * 2 + 3;
      node(7) = ( 2 * ey     ) * nx + ( ex - 1 ) * 2 + 1;
      node(8) = ( 2 * ey     ) * nx + ( ex - 1 ) * 2 + 2;
      node(9) = ( 2 * ey     ) * nx + ( ex - 1 ) * 2 + 3;

      for qx = 1 : quad_num

        xq = ( ( 1.0 - abscissa(qx) ) * xx(1)   ...
             + ( 1.0 + abscissa(qx) ) * xx(3) ) ...
               / 2.0;

        for qy = 1 : quad_num

          yq = ( ( 1.0 - abscissa(qy) ) * yy(1)   ...
               + ( 1.0 + abscissa(qy) ) * yy(3) ) ...
                 / 2.0;

          wq = weight(qx) * ( xx(3) - xx(1) ) / 2.0 ...
             * weight(qy) * ( yy(3) - yy(1) ) / 2.0;
%
%  Need to AUTOMATE THIS PROCEDURE.
%
%         v(1) =  ( xq - x(2) ) / ( x(1) - x(2) ) ...
%               * ( xq - x(3) ) / ( x(1) - x(3) ) ...
%               * ( yq - y(2) ) / ( y(1) - y(2) ) ...
%               * ( yq - y(3) ) / ( y(1) - y(3) );

%         vx(1) =           1.0 / ( x(1) - x(2) )  ...
%               * ( xq - x(3) ) / ( x(1) - x(3) ) ...
%               * ( yq - y(2) ) / ( y(1) - y(2) ) ...
%               * ( yq - y(3) ) / ( y(1) - y(3) ) ...
%               + ( xq - x(2) ) / ( x(1) - x(2) ) ...
%               *           1.0 / ( x(1) - x(3) ) ...
%               * ( yq - y(2) ) / ( y(1) - y(2) ) ...
%               * ( yq - y(3) ) / ( y(1) - y(3) );

%         vy(1) = ( xq - x(2) ) / ( x(1) - x(2) ) ...
%               * ( xq - x(3) ) / ( x(1) - x(3) ) ...
%               *          1.0  / ( y(1) - y(2) ) ...
%               * ( yq - y(3) ) / ( y(1) - y(3) ) ...
%               + ( xq - x(2) ) / ( x(1) - x(2) ) ...
%               * ( xq - x(3) ) / ( x(1) - x(3) ) ...
%               * ( yq - y(2) ) / ( y(1) - y(2) ) ...
%               *           1.0 / ( y(1) - y(3) );

          v = ones ( 9, 1 );
          vx = zeros ( 9, 1 );
          vy = zeros ( 9, 1 );

          k = 0;

          for jl = 1 : 3
            for il = 1 : 3

              k = k + 1;

              for il2 = 1 : 3
                if ( il2 ~= il )
                  v(k) = v(k) * ( xq - xx(il2) ) / ( xx(il) - xx(il2) );
                  t = 1.0 / ( xx(il) - xx(il2 ) );
                  for il3 = 1 : 3
                    if ( il3 ~= il & il3 ~= il2 )
                      t = t * ( xq - xx(il3) ) / ( xx(il) - xx(il3) );
                    end
                  end
                  for jl2 = 1 : 3
                    if ( jl2 ~= jl )
                      t = t * ( yq - yy(jl2) ) / ( yy(jl) - yy(jl2) );
                    end
                  end
                  vx(k) = vx(k) + t;
                end
              end

              for jl2 = 1 : 3
                if ( jl2 ~= jl )
                  v(k) = v(k) * ( yq - yy(jl2) ) / ( yy(jl) - yy(jl2) );
                  t = 1.0 / ( yy(jl) - yy(jl2 ) );
                  for il2 = 1 : 3
                    if ( il2 ~= il )
                      t = t * ( xq - xx(il2) ) / ( xx(il) - xx(il2) );
                    end
                  end
                  for jl3 = 1 : 3
                    if ( jl3 ~= jl & jl3 ~= jl2 )
                      t = t * ( yq - yy(jl3) ) / ( yy(jl) - yy(jl3) );
                    end
                  end
                  vy(k) = vy(k) + t;
                end
              end

            end
          end

          aq = a ( xq, yq );
          cq = c ( xq, yq );
          fq = f ( xq, yq );

          for i = 1 : 9
            ii = node(i);
            for j = 1 : 9
              jj = node(j);
              A(ii,jj) = A(ii,jj) + wq * ( vx(i) * aq * vx(j) ...
                                         + vy(i) * aq * vy(j) ...
                                         + v(i)  * cq * v(j) );
            end
            b(ii) = b(ii) + wq * ( v(i) * fq );
          end

        end
      end
    end
  end
%
%  Where a node is on the boundary,
%  replace the finite element equation by a boundary condition.
%
  k = 0;
  for y = 1 : ny
    for x = 1 : nx
      k = k + 1;
      if ( x == 1 | x == nx | y == 1 | y == ny )
        A(k,1:mn) = 0.0;
        A(1:mn,k) = 0.0;
        A(k,k) = 1.0;
        b(k) = 0.0;
      end
    end
  end

  if ( 0 )
    spy ( A );
    pause
  end
%
%  Solve the linear system.
%
  u = A \ b;
%
%  Make the vector U into a matrix.
%  Hope that of the various possible reorderings, you and MATLAB agree.
%
%    U(matrix)   U(vector)
%   ----------  ----------
%    U11 U12  =>  U11
%    U21 U22      U21
%                 U12
%                 U22
%
  u = reshape ( u, nx, ny );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function p = presolve_upforce(p)

p.upForce = {};
if p.K.l == 0
    return
end
% Extract linear inequalities
ff = p.F_struc(startofLPCone(p.K):startofLPCone(p.K)+p.K.l-1,:);
% Search for y <= x1 + x2 +...
candidates = find(ff(:,1)==0);
for i = 1:length(candidates)
    row = p.F_struc(candidates(i),2:end);
    xx = find(row==1);
    yy = find(row==-1);
    if length(yy)==1 && length(xx)+1==nnz(row)
        % Yes, this matches y>=x1+x2+x3+...
        if all(ismember(find(row),p.binary_variables))
            p.upForce{end+1}.forcing = yy;
            p.upForce{end}.forced = xx;
        end
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [ err ] = baselines( Xtrain, Ytrain, Xtest, Ytest, lambda, opts)
% Inputs
% Xtrain: input training data
% Ytrain: output training data
% Xtest: input test data
% Ytest: output test data
% lambda: regularization parameter (for global and local models)
% opts:
%   opts.avg: whether to compute avg or total error
%   opts.obj: 'R' for regression, 'C' for classification
%   opts.type:
%       consant - for each task, the model is the average (R) or mode (C) of the labels
%       global - concatenate data from all tasks, fit one global model
%       local - fit m separate models

% Output
% Average or total RMSE (R) or classification error (C) across tasks

%% set variables
m = length(Xtrain); % # of tasks
d = size(Xtrain{1}, 2); % # of features

%% compute baselines
switch opts.type
    case 'constant'
        % predict the mean output (R) or mode class (C) from training
        if(opts.obj == 'R') % regression
            if(opts.avg)
                errs = zeros(m, 1);
                for t=1:m
                    pred = mean(Ytrain{t});
                    errs(t) = sqrt(mean((Ytest{t} - pred).^2));
                end
                err = mean(errs);
            else
                Y_hat = cell(m,1);
                for t=1:m
                    Y_hat{t} = repmat(mean(Ytrain{t}), length(Ytest{t}), 1);
                end
                Y = cell2mat(Ytest(:));
                Y_hat = cell2mat(Y_hat);
                err = sqrt(mean((Y - Y_hat).^2));
            end
        else % classification
            if(opts.avg)
                errs = zeros(m, 1);
                for t=1:m
                    pred = mode(Ytrain{t});
                    errs(t) = mean(pred ~= Ytest{t});
                end
                err = mean(errs);
            else
                pred = mode(opts.Ytrainactual);
                err = sum(pred ~= opts.Yactual) / length(opts.Yactual);
            end
        end

    case 'global'
        % compute one global model for the data
        allX = cat(1, Xtrain{:});
        allY = cat(1, Ytrain{:});
        allXtest = cat(1, Xtest{:});
        allYtest = cat(1, Ytest{:});
        if(opts.obj == 'R') % regression
            w = inv(allX' * allX + lambda * eye(d)) * allX' * allY;
            err = sqrt(mean((allYtest - allXtest * w).^2));
        else % classification
            w = simple_svm(allX, allY, lambda, opts);
            if(opts.avg)
                errs = zeros(m, 1);
                for t=1:m
                    predvals = sign(Xtest{t} * w);
                    errs(t) = mean(predvals ~= Ytest{t});
                end
                err = mean(errs);
            else
                predvals = sign(allXtest * w);
                err = mean(allYtest ~= predvals);
            end
        end

    case 'local'
        % compute m completely separate local models
        if(opts.obj == 'R') % regression
            if(opts.avg)
                errs = zeros(m,1);
                for t=1:m
                    wt = inv(Xtrain{t}' * Xtrain{t} + lambda * eye(d)) * Xtrain{t}' * Ytrain{t};
                    errs(t) = sqrt(mean((Ytest{t} - Xtest{t} * wt).^2));
                end
                % compute average rmse
                err = mean(errs);
            else
                Y_hat = cell(m,1);
                for t=1:m
                    wt = inv(Xtrain{t}' * Xtrain{t} + lambda * eye(d)) * Xtrain{t}' * Ytrain{t};
                    Y_hat{t} = Xtest{t} * wt;
                end
                % compute total rmse
                Y = cell2mat(Ytest(:));
                Y_hat = cell2mat(Y_hat);
                size(Y_hat)
                size(Y)
                err = sqrt(mean((Y - Y_hat).^2));
            end
        else % classification
            Y_hat = cell(m,1);
            for t=1:m
                wt = simple_svm(Xtrain{t}, Ytrain{t}, lambda, opts);
                Y_hat{t} = sign(Xtest{t} * wt);
            end
            if(opts.avg)
                errs = zeros(m,1);
                for t=1:m
                    errs(t) = mean(Ytest{t} ~= Y_hat{t});
                end
                err = mean(errs);
            else
                Y = cell2mat(Ytest);
                Y_hat = cell2mat(Y_hat);
                err = mean(Y ~= Y_hat);
            end
        end
end

end

author: gingsmith

repo: fmtl

save_path: github-repos/MATLAB/gingsmith-fmtl


text:
function [P, J] = regionGrowing(cIM, initPos, thresVal, maxDist, tfMean, tfFillHoles, tfSimplify)
% REGIONGROWING Region growing algorithm for 2D/3D grayscale images
%
% Syntax:
%   P = regionGrowing();
%   P = regionGrowing(cIM);
%   P = regionGrowing(cIM, initPos)
%   P = regionGrowing(..., thresVal, maxDist, tfMean, tfFillHoles, tfSimpl)
%   [P, J] = regionGrowing(...);
%
% Inputs:
%          cIM: 2D/3D grayscale matrix                      {current image}
%      initPos: Coordinates for initial seed position     {ginput position}
%     thresVal: Absolute threshold level to be included     {5% of max-min}
%      maxDist: Maximum distance to the initial position in [px]      {Inf}
%       tfMean: Updates the initial value to the region mean (slow) {false}
%  tfFillHoles: Fills enclosed holes in the binary mask              {true}
%   tfSimplify: Reduces the number of vertices {true, if dpsimplify exists}
%
% Outputs:
%   P: VxN array (with V number of vertices, N number of dimensions)
%      P is the enclosing polygon for all associated pixel/voxel
%   J: Binary mask (with the same size as the input image) indicating
%      1 (true) for associated pixel/voxel and 0 (false) for outside
%
% Examples:
%   % 2D Example
%   load example
%   figure, imshow(cIM, [0 1500]), hold all
%   poly = regionGrowing(cIM, [], 300); % click somewhere inside the lungs
%   plot(poly(:,1), poly(:,2), 'LineWidth', 2)
%
%   % 3D Example
%   load mri
%   poly = regionGrowing(squeeze(D), [66,55,13], 60, Inf, [], true, false);
%   plot3(poly(:,1), poly(:,2), poly(:,3), 'x', 'LineWidth', 2)
%
% Requirements:
%   TheMathWorks Image Processing Toolbox for bwboundaries() and axes2pix()
%   Optional: Line Simplification by Wolfgang Schwanghart to reduce the
%             number of polygon vertices (see the MATLAB FileExchange)
%
% Remarks:
%   The queue is not preallocated and the region mean computation is slow.
%   I haven't implemented a preallocation nor a queue counter yet for the
%   sake of clarity, however this would be of course more efficient.
%
% Author:
%   Daniel Kellner, 2011, braggpeaks{}googlemail.com
%   History: v1.00: 2011/08/14

% error checking on input arguments
if nargin > 7
    error('Wrong number of input arguments!')
end

if ~exist('cIM', 'var')
    himage = findobj('Type', 'image');
    if isempty(himage) || length(himage) > 1
        error('Please define one of the current images!')
    end

    cIM = get(himage, 'CData');
end

if ~exist('initPos', 'var') || isempty(initPos)
    himage = findobj('Type', 'image');
    if isempty(himage)
        himage = imshow(cIM, []);
    end

    % graphical user input for the initial position
    p = ginput(1);

    % get the pixel position concerning to the current axes coordinates
    initPos(1) = round(axes2pix(size(cIM, 2), get(himage, 'XData'), p(2)));
    initPos(2) = round(axes2pix(size(cIM, 1), get(himage, 'YData'), p(1)));
end

if ~exist('thresVal', 'var') || isempty(thresVal)
    thresVal = double((max(cIM(:)) - min(cIM(:)))) * 0.05;
end

if ~exist('maxDist', 'var') || isempty(maxDist)
    maxDist = Inf;
end

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

if ~exist('tfFillHoles', 'var')
    tfFillHoles = true;
end

if isequal(ndims(cIM), 2)
    initPos(3) = 1;
elseif isequal(ndims(cIM),1) || ndims(cIM) > 3
    error('There are only 2D images and 3D image sets allowed!')
end

[nRow, nCol, nSli] = size(cIM);

if initPos(1) < 1 || initPos(2) < 1 ||...
   initPos(1) > nRow || initPos(2) > nCol
    error('Initial position out of bounds, please try again!')
end

if thresVal < 0 || maxDist < 0
    error('Threshold and maximum distance values must be positive!')
end

if ~isempty(which('dpsimplify.m'))
    if ~exist('tfSimplify', 'var')
        tfSimplify = true;
    end
    simplifyTolerance = 1;
else
    tfSimplify = false;
end

% initial pixel value
regVal = double(cIM(initPos(1), initPos(2), initPos(3)));

% text output with initial parameters
disp(['RegionGrowing Opening: Initial position (' num2str(initPos(1))...
      '|' num2str(initPos(2)) '|' num2str(initPos(3)) ') with '...
      num2str(regVal) ' as initial pixel value!'])

% preallocate array
J = false(nRow, nCol, nSli);

% add the initial pixel to the queue
queue = [initPos(1), initPos(2), initPos(3)];

%%% START OF REGION GROWING ALGORITHM
while size(queue, 1)
  % the first queue position determines the new values
  xv = queue(1,1);
  yv = queue(1,2);
  zv = queue(1,3);

  % .. and delete the first queue position
  queue(1,:) = [];

  % check the neighbors for the current position
  for i = -1:1
    for j = -1:1
      for k = -1:1

        if xv+i > 0  &&  xv+i <= nRow &&...          % within the x-bounds?
           yv+j > 0  &&  yv+j <= nCol &&...          % within the y-bounds?
           zv+k > 0  &&  zv+k <= nSli &&...          % within the z-bounds?
           any([i, j, k])       &&...      % i/j/k of (0/0/0) is redundant!
           ~J(xv+i, yv+j, zv+k) &&...          % pixelposition already set?
           sqrt( (xv+i-initPos(1))^2 +...
                 (yv+j-initPos(2))^2 +...
                 (zv+k-initPos(3))^2 ) < maxDist &&...   % within distance?
           cIM(xv+i, yv+j, zv+k) <= (regVal + thresVal) &&...% within range
           cIM(xv+i, yv+j, zv+k) >= (regVal - thresVal) % of the threshold?

           % current pixel is true, if all properties are fullfilled
           J(xv+i, yv+j, zv+k) = true;

           % add the current pixel to the computation queue (recursive)
           queue(end+1,:) = [xv+i, yv+j, zv+k];

           if tfMean
               regVal = mean(mean(cIM(J > 0))); % --> slow!
           end

        end
      end
    end
  end
end
%%% END OF REGION GROWING ALGORITHM

% loop through each slice, fill holes and extract the polygon vertices
P = [];
for cSli = 1:nSli
    if ~any(J(:,:,cSli))
        continue
    end

	% use bwboundaries() to extract the enclosing polygon
    if tfFillHoles
        % fill the holes inside the mask
        J(:,:,cSli) = imfill(J(:,:,cSli), 'holes');
        B = bwboundaries(J(:,:,cSli), 8, 'noholes');
    else
        B = bwboundaries(J(:,:,cSli));
    end

	newVertices = [B{1}(:,2), B{1}(:,1)];

    % simplify the polygon via Line Simplification
    if tfSimplify
        newVertices = dpsimplify(newVertices, simplifyTolerance);
    end

    % number of new vertices to be added
    nNew = size(newVertices, 1);

    % append the new vertices to the existing polygon matrix
    if isequal(nSli, 1) % 2D
        P(end+1:end+nNew, :) = newVertices;
    else                % 3D
        P(end+1:end+nNew, :) = [newVertices, repmat(cSli, nNew, 1)];
    end
end

% text output with final number of vertices
disp(['RegionGrowing Ending: Found ' num2str(length(find(J)))...
      ' pixels within the threshold range (' num2str(size(P, 1))...
      ' polygon vertices)!'])

author: Sable

repo: mcbench-benchmarks


text:
function A = warmUpExercise()
%WARMUPEXERCISE Example function in octave
%   A = WARMUPEXERCISE() is an example function that returns the 5x5 identity matrix

A = [];
% ============= YOUR CODE HERE ==============
% Instructions: Return the 5x5 identity matrix
%               In octave, we return values by defining which variables
%               represent the return values (at the top of the file)
%               and then set them accordingly.

A = eye(5)    % 5x5 identity matrix

% ===========================================

end

repo: machine-learning-coursera


text:
function trans = scaling3d(varargin)
%SCALING3D return 4x4 matrix of a 3D scaling
%
%   TRANS = scaling3d(S);
%   returns the scaling transform corresponding to a scaling factor S in
%   each direction. S can be a scalar, or a 1x3 vector containing the
%   scaling factor in each direction.
%
%   TRANS = scaling3d(SX, SY, SZ);
%   returns the scaling transform corresponding to a different scaling
%   factor in each direction.
%
%   The returned matrix has the form :
%   [SX  0  0  0]
%   [ 0 SY  0  0]
%   [ 0  0 SZ  0]
%   [ 0  0  0  0]
%
%
%   See also:
%   transforms3d, transformPoint3d
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 20/04/2006.
%

%   HISTORY
%   25/11/2008 rename from scale3d to scaling3d
%   HISTORY
%   30/04/2009 deprecate: use createScaling3d instead

% deprecation warning
warning('geom3d:deprecated', ...
    [mfilename ' is deprecated, use ''createScaling3d'' instead']);

% call current implementation
trans = createScaling3d(varargin{:});

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [ DKI, K_ADC, K_AKC, Kv1, Kv2, Kv3, I_DKI ] = fit_DKI_model( im, B_mat, thresh_val, method)
% FIT_DT    Fits the diffusion kurtosis model voxelwise to an image
%           S = S0 exp(-(bD + 1/6 b^2 * K))
%
%
% Inputs:
%
%   IM is the input image, of any dimensionality, with the diffusion
%   scans in the last dimension
%
%   B is the b-matrix, of size [3 3 N]. If the measured b matrix isn't
%   available, it can be approximated by:
%   b(:,:,n)=b_values(n)*unit_vectors(n,:)'*unit_vectors(n,:);
%
%   THRESH_VAL is a threshold for skipping the FA, ADC and vector field
%   computation in voxels x where im(x) < thresh_val. THRESH_VAL can
%   alternatively be a mask of voxels to process.
%
%   METHOD is a string, and can take the following values:
%       'linear' fits the tensor in the log domain (default)
%       'nonlinear' fits the tensor in the signal domain (recommended for
%       high b-values or when SNR is low)
%
% Outputs:
%
%   DKI is the array of coefficients

% Author: Darryl McClymont <darryl.mcclymont@gmail.com>
% Copyright © 2014-2015 University of Oxford
% Version: 0.1.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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 3 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
%narginchk(2, 5);
%nargoutchk(0, 5);

sz = size(im);

if (nargin < 3) || isempty(thresh_val)
    thresh_val = -inf;
end

if (nargin < 4) || isempty(method)% default method is linear least squares
    method = 'linear';
end

% check that a valid method has been given
if ~max([strcmp(method, 'linear'), strcmp(method, 'nonlinear')])
    disp('Unrecognised method, performing fast linear fitting')
end

% handle a vector with the wrong orientation
if (length(sz) == 2) && (sz(2) == 1)
    im = im';
    sz = size(im);
end

% reshape image
I_vector = reshape(im, [prod(sz(1:end-1)), sz(end)]);

% handle mask or threshold
if isscalar(thresh_val)
    thresh_val = I_vector(:,1) > thresh_val;
else
    thresh_val = logical(thresh_val);
end

% keep only the masked voxels
I_vector = double(I_vector(thresh_val(:), :));

[~, b0] = max(mean(I_vector, 1));

% normalise and take the log of the image to linearise the equation
imlog = log(abs(bsxfun(@rdivide, I_vector, I_vector(:,b0))));

% get the bvalue
bvalue = squeeze(B_mat(1,1,:) + B_mat(2,2,:) + B_mat(3,3,:))';

% normalise the B matrix
B_mat_norm = B_mat ./ repmat(permute(bvalue + eps, [3,1,2]), [3 3 1]);

% define A matrix

A_D = squeeze([B_mat_norm(1,1,:), B_mat_norm(2,2,:), B_mat_norm(3,3,:), 2*B_mat_norm(1,2,:), 2*B_mat_norm(1,3,:), 2*B_mat_norm(2,3,:)])';
% A_D in the paper is not scaled by the b value

A_K = squeeze([B_mat_norm(1,1,:).^2, B_mat_norm(2,2,:).^2, B_mat_norm(3,3,:).^2, ...
    4*B_mat_norm(1,1,:).*B_mat_norm(1,2,:), ...
    4*B_mat_norm(1,1,:).*B_mat_norm(1,3,:), ...
    4*B_mat_norm(1,2,:).*B_mat_norm(2,2,:), ...
    4*B_mat_norm(1,3,:).*B_mat_norm(3,3,:), ...
    4*B_mat_norm(2,2,:).*B_mat_norm(2,3,:), ...
    4*B_mat_norm(2,3,:).*B_mat_norm(3,3,:), ...
    2*B_mat_norm(1,1,:).*B_mat_norm(2,2,:) + 4*B_mat_norm(1,2,:).*B_mat_norm(1,2,:), ...
    2*B_mat_norm(1,1,:).*B_mat_norm(3,3,:) + 4*B_mat_norm(1,3,:).*B_mat_norm(1,3,:), ...
    2*B_mat_norm(2,2,:).*B_mat_norm(3,3,:) + 4*B_mat_norm(2,3,:).*B_mat_norm(2,3,:), ...
    8*B_mat_norm(1,2,:).*B_mat_norm(1,3,:) + 4*B_mat_norm(1,1,:).*B_mat_norm(2,3,:), ...
    8*B_mat_norm(1,2,:).*B_mat_norm(2,3,:) + 4*B_mat_norm(2,2,:).*B_mat_norm(1,3,:), ...
    8*B_mat_norm(1,3,:).*B_mat_norm(2,3,:) + 4*B_mat_norm(3,3,:).*B_mat_norm(1,2,:)])';
% this takes care of cross terms by saying that (for example) 6 x^2 y^2
% = 2 (x^2) (y^2) + 4 (xy) (xy)
% so it uses these elements of the B matrix.

A = [-bsxfun(@times, bvalue', A_D), 1/6 * bsxfun(@times, bvalue'.^2, A_K)];

% only keep b values less than 5000 in the linear fit
lin_bvals = bvalue < 5000;

% weighted fit
param.rician = 0;
param.unique_weights = 0;
param.verbose = 1;
M = weighted_linear_fit(imlog(:,lin_bvals), A(lin_bvals,:)', @(z)exp(z), param);
%M = weighted_linear_fit(imlog, A', @(z)exp(z), param);

%M = (pinv(A) * imlog')';
%M = (pinv(A' * (W.^2) * A) * A' * (W.^2) * (imlog'))';
% weighted version of M = (pinv(A) * imlog')';

% add the S0 column
M = [I_vector(:,1), M];

M_lin = M;

% perform nonlinear fitting (if required)
if strcmp(method, 'nonlinear')

    % get error
    I_recon = bsxfun(@times, M(:,1) , exp(M(:,2:end) * A'));
    rmse_lin = sqrt(mean((I_vector(:) - I_recon(:)).^2));
    %disp(['RMSE from quick fit = ' num2str(rmse_lin)])

    options = optimoptions('lsqcurvefit','Jacobian','on', 'DerivativeCheck', 'off', ...
        'display', 'off', 'MaxIter', 40);

    lb = [0 0 0 0 -1 -1 -1 ... % [S0 DTI
            0 0 0 -1 -1 -1 -1 -1 -1 0 0 0 -1 -1 -1] * 3E-3; % DKI ]

    ub = [inf, ones(1,21) * 3E-3]; % [S0 DTI DKI]

    % constrain initial guess between lower and upper bounds
    M = bsxfun(@max, M, lb);
    M = bsxfun(@min, M, ub);
    M(isnan(M(:))) = 0;
    M(isinf(M(:))) = 0;

    if (size(I_vector,1) > 1000) && (matlabpool('size') == 0)
        disp('Open matlabpool for parallel processing')
    end

    A_hat = -A';

    parfor i = 1:size(I_vector,1)

        M(i,:) = lsqcurvefit(@DKI_model, M(i,:), A_hat, I_vector(i,:), lb, ub, options);

    end

    I_recon = bsxfun(@times, M(:,1) , exp(M(:,2:end) * A'));
    rmse = sqrt(mean((I_vector(:) - I_recon(:)).^2));

    bad_fits = rmse > rmse_lin;
    %disp(['Number of bad fits: ' num2str(sum(bad_fits(:)))])
    M(bad_fits,:) = M_lin(bad_fits,:);

    %disp(['RMSE from slow fit = ' num2str(rmse)])

end

if nargout == 1
    % quick hack because reshape sz needs to be at least length 2
    if length(sz) == 2
        sz = [sz(1), 1, sz(2)];
    end

    % fill regions outside mask with zeros, and reshape
    DKI = zeros(prod(sz(1:end-1)), 22);
    DKI(thresh_val(:), :) = M;
    DKI = reshape(DKI, [sz(1:end-1), 22]);

    return
end

% compute the returned signal
I_DKI = bsxfun(@times, M(:,1), exp(M(:, 2:end) * A'));

K_v1 = zeros(1, size(M,1));
K_v2 = zeros(1, size(M,1));
K_v3 = zeros(1, size(M,1));
ADC = zeros(1, size(M,1));
AKC = zeros(1, size(M,1));

[x, y, z] = psphere(200);
x = x(1:end/2); y = y(1:end/2); z = z(1:end/2);

% Compute the kurtosis (mean and directional)
for v = 1:size(M,1)

    % reshape the coefficients

    % S0 = M(i,1);

    D_hat_i = M(v,2:7);

    % A_D = xx, yy, zz, xy, xz, yz

    % rearrange for convenience
    D_hat = D_hat_i([1 4 5 2 6 3]);

    D33 = [D_hat(1), D_hat(2), D_hat(3);
           D_hat(2), D_hat(4), D_hat(5);
           D_hat(3), D_hat(5), D_hat(6)];

    [vec, val] = eig(D33);
    % sort
    [~, ix] = sort(sum(val,1), 'descend');
    vec = vec(:, ix);
    val = sum(val(:, ix));

    % eigenvalues should be positive
    % taking the magnitude is a very bad way of handling them
    val = abs(val);

    ADC(v) = sum(val) / 3;

    D2i = M(v, 8:end);

    % A_K = xxxx, yyyy, zzzz, 1 2 3
    %       xxxy, xxxz, xyyy, 4 5 6
    %       xzzz, yyyz, yzzz, 7 8 9
    %       xxyy, xxzz, yyzz, 10 11 12
    %       xxyz, xyyz, xyzz ,13 14 15

    % re-arrange for convenience
    D2 = D2i([1, 4, 5, 10, 13, 11, 6, 14, 15, 7, 2, 8, 12, 9, 3]);

    % kurtosis tensor (includes ADC^2 internally)
    K = [ D2(01), D2(02), D2(03), D2(04), D2(05), D2(06); % xx
          D2(02), D2(04), D2(05), D2(07), D2(08), D2(09); % xy
          D2(03), D2(05), D2(06), D2(08), D2(09), D2(10); % xz
          D2(04), D2(07), D2(08), D2(11), D2(12), D2(13); % yy
          D2(05), D2(08), D2(09), D2(12), D2(13), D2(14); % yz
          D2(06), D2(09), D2(10), D2(13), D2(14), D2(15)];% zz
        %   xx      xy      xz      yy      yz      zz

    % Parallel kurtosis (kurtosis in the direction of highest diffusivity)
    v1 = vec(:, 1);

    % Diffusion in this direction
    D_v1 = v1' * D33 * v1; % = eigenvalue 1

    % kurtosis in this direction
    v1_sq = v1 * v1';
    v1_sq_vect = [v1_sq(1), 2*v1_sq(2), 2*v1_sq(3), v1_sq(5), 2*v1_sq(6), v1_sq(9)];
    K_v1(v) = (v1_sq_vect * K * v1_sq_vect') / (D_v1^2);

    % Parallel kurtosis (v2)
    v2 = vec(:, 2);

    % Diffusion in this direction
    D_v2 = v2' * D33 * v2; % = eigenvalue 2

    % kurtosis in this direction
    v2_sq = v2 * v2';
    v2_sq_vect = [v2_sq(1), 2*v2_sq(2), 2*v2_sq(3), v2_sq(5), 2*v2_sq(6), v2_sq(9)];
    K_v2(v) = (v2_sq_vect * K * v2_sq_vect') / (D_v2^2);

    % Parallel kurtosis (v3)
    v3 = vec(:, 3);

    % Diffusion in this direction
    D_v3 = v3' * D33 * v3; % = eigenvalue 3

    % kurtosis in this direction
    v3_sq = v3 * v3';
    v3_sq_vect = [v3_sq(1), 2*v3_sq(2), 2*v3_sq(3), v3_sq(5), 2*v3_sq(6), v3_sq(9)];
    K_v3(v) = (v3_sq_vect * K * v3_sq_vect') / (D_v3^2);

    % mean kurtosis
    % there are clever ways of getting this, but the easiest is to just
    % average the kurtosis over a bunch of uniformly spaced directions
    K_vi = zeros(1,length(x));
    for i = 1:length(x)
        vi = [x(i), y(i), z(i)]';
        D_vi = vi' * D33 * vi;
        vi_sq = vi * vi';
        vi_sq_vect = [vi_sq(1), 2*vi_sq(2), 2*vi_sq(3), vi_sq(5), 2*vi_sq(6), vi_sq(9)];
        K_vi(i) = (vi_sq_vect * K * vi_sq_vect') / (D_vi^2);
    end

    AKC(v) = mean(K_vi);

end

% quick hack because reshape sz needs to be at least length 2
if length(sz) == 2
    sz = [sz(1), 1, sz(2)];
end

% fill regions outside mask with zeros, and reshape
DKI = zeros(prod(sz(1:end-1)), 22);
DKI(thresh_val(:), :) = M;
DKI = reshape(DKI, [sz(1:end-1), 22]);

IDKI = zeros(prod(sz(1:end-1)), sz(end));
IDKI(thresh_val(:), :) = I_DKI;
I_DKI = reshape(IDKI, sz);

K_ADC = zeros(size(thresh_val));
K_ADC(thresh_val) = ADC;
K_ADC = reshape(K_ADC, sz(1:end-1));

K_AKC = zeros(size(thresh_val));
K_AKC(thresh_val) = AKC;
K_AKC = reshape(K_AKC, sz(1:end-1));

Kv1 = zeros(size(thresh_val));
Kv1(thresh_val) = K_v1;
Kv1 = reshape(Kv1, sz(1:end-1));

Kv2 = zeros(size(thresh_val));
Kv2(thresh_val) = K_v2;
Kv2 = reshape(Kv2, sz(1:end-1));

Kv3 = zeros(size(thresh_val));
Kv3(thresh_val) = K_v3;
Kv3 = reshape(Kv3, sz(1:end-1));

end

function [ F, J ] = DKI_model( x, Bv )

    F = x(1) * exp(-(x(2:end) * Bv));

    if nargout > 1 % Jacobian
        J = zeros(length(F), length(x));

        J(:,1) = exp(-(x(2:end) * Bv));

        J(:,2:end) = -bsxfun(@times, Bv, F)';

    end

end

%% This converts the 3x3 kurtosis tensor into a 3x3x3x3 kurtosis tensor
% It is not needed right now, but good to have for the future.

%     K_index =[01, 02, 03, 04, 05, 06; % xx
%               02, 04, 05, 07, 08, 09; % xy
%               03, 05, 06, 08, 09, 10; % xz
%               04, 07, 08, 11, 12, 13; % yy
%               05, 08, 09, 12, 13, 14; % yz
%               06, 09, 10, 13, 14, 15];% zz
%             % xx  xy  xz  yy  yz  zz
%
%     W = zeros([3 3 3 3]);
%
%     for i = 1:3
%         for j = 1:3
%              for k = 1:3
%                  for l = 1:3
%                     switch i
%                         case 1
%                             r1 = [1 2 3];
%                         case 2
%                             r1 = [2 4 5];
%                         case 3
%                             r1 = [3 5 6];
%                     end
%                     switch j
%                         case 1
%                             r2 = [1 2 3];
%                         case 2
%                             r2 = [2 4 5];
%                         case 3
%                             r2 = [3 5 6];
%                     end
%
%                     if i == j
%                         switch i
%                             case 1
%                                 r =1;
%                             case 2
%                                 r = 4;
%                             case 3
%                                 r = 6;
%                         end
%                     else
%                         r = intersect(r1, r2);
%                     end
%
%
%                     switch k
%                         case 1
%                             c1 = [1 2 3];
%                         case 2
%                             c1 = [2 4 5];
%                         case 3
%                             c1 = [3 5 6];
%                     end
%                     switch l
%                         case 1
%                             c2 = [1 2 3];
%                         case 2
%                             c2 = [2 4 5];
%                         case 3
%                             c2 = [3 5 6];
%                     end
%
%                     if k == l
%                         switch k
%                             case 1
%                                 c =1;
%                             case 2
%                                 c = 4;
%                             case 3
%                                 c = 6;
%                         end
%                     else
%                         c = intersect(c1, c2);
%                     end
%
%                     D2_index = K_index(r,c);
%
%                     W(i,j,k,l) = D2(D2_index);
%                  end
%              end
%         end
%     end
%     v1 = vec(:, 1);
%     % Diffusion in this direction
%     D_v1 = v1' * D33 * v1;
%     K_v1_vox = 0;
%     for i = 1:3
%         for j = 1:3
%             for k = 1:3
%                 for l = 1:3
%                     K_v1_vox = K_v1_vox + v1(i) * v1(j) * v1(k) * v1(l) * W(i,j,k,l);
%                 end
%             end
%         end
%     end
%     K_v1(v) = K_v1_vox / D_v1^2;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% Change to your downloaded location
clear
addpath('C:\liblinear\matlab')
addpath('../training_code/');
addpath('../utilities/');
addpath('../../data extraction/');

%% load shared definitions and AU data
all_dataset_aus = [1, 2, 4, 5, 6, 7, 9, 10, 12, 14, 15, 17, 20, 23, 25, 26, 28, 45];

% Set up the hyperparameters to be validated
hyperparams.c = 10.^(-7:0.5:1);
hyperparams.e = 10.^(-3);

hyperparams.validate_params = {'c', 'e'};

% Set the training function
svm_train = @svm_train_linear;

% Set the test function (the first output will be used for validation)
svm_test = @svm_test_linear;

BP4D_aus = [1, 2, 4, 6, 7, 10, 12, 14, 15, 17, 23];
SEMAINE_aus = [2, 12, 17, 25, 28, 45];
FERA2011_aus = [1, 2, 4, 6, 7, 10, 12, 15, 17, 25, 26];
UNBC_aus = [6, 7, 9, 12, 25, 26];
DISFA_aus = [1, 2, 4, 5, 6, 9, 12, 15, 17, 20, 25, 26];
Bosphorus_aus = [1, 2, 4, 5, 6, 7, 9, 10, 12, 14, 15, 17, 20, 23, 25, 26, 45];

for a=1:numel(all_dataset_aus)

    au = all_dataset_aus(a);

    train_samples = [];
    valid_samples = [];

    train_labels = [];
    valid_labels = [];

    % Keeping track which validation sample came from which dataset, as we
    % will be validating based on which hyperparam leads to best
    % performance on all datasets (mean F1)
    eval_ids = [];
    dataset_ids = {};

    if(~isempty(find(BP4D_aus == au, 1)))
        op = cd('../BP4D/');
        rest_aus = setdiff(BP4D_aus, au);
        shared_defs;
        % load the training and testing data for the current fold
        [train_samples_bp4d, train_labels_bp4d, valid_samples_bp4d, valid_labels_bp4d, ~, PC, means, scaling] = Prepare_HOG_AU_data_generic_dynamic(train_recs, devel_recs, au, BP4D_dir, hog_data_dir);

        train_samples = cat(1, train_samples, train_samples_bp4d);
        valid_samples = cat(1, valid_samples, valid_samples_bp4d);

        train_labels = cat(1, train_labels, train_labels_bp4d);
        valid_labels = cat(1, valid_labels, valid_labels_bp4d);

        if(isempty(eval_ids))
           eval_ids = ones(size(valid_labels_bp4d,1), 1);
        end

        clear 'train_samples_bp4d' 'train_labels_bp4d' 'valid_samples_bp4d' 'valid_labels_bp4d'

        dataset_ids = cat(1, dataset_ids, {'BP4D'});

        cd(op);

    end

    if(~isempty(find(FERA2011_aus == au, 1)))
        op = cd('../FERA2011/');
        rest_aus = setdiff(FERA2011_aus, au);
        shared_defs;
        all_recs = cat(2, train_recs, devel_recs);
        [users_train, users_valid_fera] = get_balanced_fold(FERA2011_dir, all_recs, au, 1/3, 1);

        % load the training and testing data for the current fold
        [train_samples_f2011, train_labels_f2011, valid_samples_f2011, valid_labels_f2011, ~, PC, means, scaling] = Prepare_HOG_AU_data_generic_dynamic(users_train, users_valid_fera, au, rest_aus, FERA2011_dir, features_dir);

        train_samples = cat(1, train_samples, train_samples_f2011);
        valid_samples = cat(1, valid_samples, valid_samples_f2011);

        train_labels = cat(1, train_labels, train_labels_f2011);
        valid_labels = cat(1, valid_labels, valid_labels_f2011);

        if(isempty(eval_ids))
           eval_ids = ones(size(valid_labels_f2011,1), 1);
        else
           eval_ids = cat(1, eval_ids, (eval_ids(end)+1)*ones(size(valid_labels_f2011,1), 1));
        end

        clear 'train_samples_f2011' 'train_labels_f2011' 'valid_samples_f2011' 'valid_labels_bp4d'

        dataset_ids = cat(1, dataset_ids, {'FERA2011'});

        cd(op);

    end

    if(~isempty(find(SEMAINE_aus == au, 1)))
        op = cd('../SEMAINE/');
        rest_aus = setdiff(SEMAINE_aus, au);
        shared_defs;
        % load the training and testing data for the current fold
        [train_samples_sem, train_labels_sem, valid_samples_sem, valid_labels_sem, ~, PC, means, scaling] = Prepare_HOG_AU_data_generic_dynamic(train_recs, devel_recs, au, rest_aus, SEMAINE_dir, hog_data_dir);

        train_samples = cat(1, train_samples, train_samples_sem);
        valid_samples = cat(1, valid_samples, valid_samples_sem);

        train_labels = cat(1, train_labels, train_labels_sem);
        valid_labels = cat(1, valid_labels, valid_labels_sem);

        if(isempty(eval_ids))
           eval_ids = ones(size(valid_labels_sem,1), 1);
        else
           eval_ids = cat(1, eval_ids, (eval_ids(end)+1)*ones(size(valid_labels_sem,1), 1));
        end

        clear 'train_samples_sem' 'train_labels_sem' 'valid_samples_sem' 'valid_labels_sem'
        dataset_ids = cat(1, dataset_ids, {'SEMAINE'});
        cd(op);

    end

    if(~isempty(find(Bosphorus_aus == au, 1)))
        op = cd('../Bosphorus/');
        rest_aus = setdiff(DISFA_aus, au);
        shared_defs;

        % make sure validation data's labels are balanced
        [users_train, users_valid_bosph] = get_balanced_fold(Bosphorus_dir, all_recs, au, 1/3, 1);

        % need to split the rest
        [train_samples_bosph, train_labels_bosph, valid_samples_bosph, valid_labels_bosph, ~, PC, means, scaling, valid_ids, valid_success] = Prepare_HOG_AU_data_dynamic(users_train, users_valid_bosph, au, rest_aus, Bosphorus_dir, hog_data_dir);

        train_labels_bosph(train_labels_bosph > 1) = 1;
        valid_labels_bosph(valid_labels_bosph > 1) = 1;

        % As there are only a few Bosphorus images (it being an image dataset)
        % we oversample it quite a bit
        train_samples_bosph = repmat(train_samples_bosph, 10, 1);
        train_labels_bosph = repmat(train_labels_bosph, 10, 1);

        train_samples = cat(1, train_samples, train_samples_bosph);
        valid_samples = cat(1, valid_samples, valid_samples_bosph);

        train_labels = cat(1, train_labels, train_labels_bosph);
        valid_labels = cat(1, valid_labels, valid_labels_bosph);

        if(isempty(eval_ids))
           eval_ids = ones(size(valid_labels_bosph, 1), 1);
        else
           eval_ids = cat(1, eval_ids, (eval_ids(end)+1)*ones(size(valid_labels_bosph, 1), 1));
        end

        clear 'train_samples_bosph' 'train_labels_bosph' 'valid_samples_bosph' 'valid_labels_bosph'
        dataset_ids = cat(1, dataset_ids, {'Bosphorus'});
        cd(op);

    end

    if(~isempty(find(DISFA_aus == au, 1)))
        op = cd('../DISFA/');
        rest_aus = setdiff(DISFA_aus, au);
        shared_defs;

        % make sure validation data's labels are balanced
        [users_train, users_valid_disfa] = get_balanced_fold(DISFA_dir, users, au, 1/3, 1);

        % need to split the rest
        [train_samples_disf, train_labels_disf, valid_samples_disf, valid_labels_disf, ~, PC, means, scaling, valid_ids, valid_success] = Prepare_HOG_AU_data_generic_dynamic(users_train, users_valid_disfa, au, rest_aus, DISFA_dir, hog_data_dir);

        train_labels_disf(train_labels_disf > 1) = 1;
        valid_labels_disf(valid_labels_disf > 1) = 1;

        train_samples = cat(1, train_samples, train_samples_disf);
        valid_samples = cat(1, valid_samples, valid_samples_disf);

        train_labels = cat(1, train_labels, train_labels_disf);
        valid_labels = cat(1, valid_labels, valid_labels_disf);

        if(isempty(eval_ids))
           eval_ids = ones(size(valid_labels_disf,1), 1);
        else
           eval_ids = cat(1, eval_ids, (eval_ids(end)+1)*ones(size(valid_labels_disf,1), 1));
        end

        clear 'train_samples_disf' 'train_labels_disf' 'valid_samples_disf' 'valid_labels_disf'
        dataset_ids = cat(1, dataset_ids, {'DISFA'});
        cd(op);

    end

    if(~isempty(find(UNBC_aus == au, 1)))
        op = cd('../UNBC/');
        rest_aus = setdiff(UNBC_aus, au);
        shared_defs;

        [users_train, users_valid_unbc] = get_balanced_fold(UNBC_dir, all_recs, au, 1/3, 1);

        % load the training and testing data for the current fold
        [train_samples_unbc, train_labels_unbc, valid_samples_unbc, valid_labels_unbc, ~, PC, means, scaling] = Prepare_HOG_AU_data_dynamic(users_train, users_valid_unbc, au, rest_aus, UNBC_dir, features_dir);

        % Binarizing the data
        train_labels_unbc(train_labels_unbc > 1) = 1;
        valid_labels_unbc(valid_labels_unbc > 1) = 1;

        train_samples = cat(1, train_samples, train_samples_unbc);
        valid_samples = cat(1, valid_samples, valid_samples_unbc);

        train_labels = cat(1, train_labels, train_labels_unbc);
        valid_labels = cat(1, valid_labels, valid_labels_unbc);

        if(isempty(eval_ids))
           eval_ids = ones(size(valid_labels_unbc,1), 1);
        else
           eval_ids = cat(1, eval_ids, (eval_ids(end)+1)*ones(size(valid_labels_unbc,1), 1));
        end

        clear 'train_samples_unbc' 'train_labels_unbc' 'valid_samples_unbc' 'valid_labels_unbc'
        dataset_ids = cat(1, dataset_ids, {'UNBC'});
        cd(op);

    end

    train_samples = sparse(train_samples);
    dataset_ids
    %% Cross-validate here
    hyperparams.eval_ids = eval_ids;
    [ best_params, ~ ] = validate_grid_search_no_par(svm_train, svm_test, false, train_samples, train_labels, valid_samples, valid_labels, hyperparams);
    model = svm_train(train_labels, train_samples, best_params);

    model.eval_ids = eval_ids;

    [~, predictions_all] = svm_test(valid_labels, valid_samples, model);

    name = sprintf('mat_models/AU_%d_dynamic.mat', au);

    eval_ids_uq = unique(eval_ids)';
    F1s = [];
    for i=eval_ids_uq
        F1s = cat(2, F1s, compute_F1(valid_labels(eval_ids==i), predictions_all(eval_ids==i)));
    end

    meanF1 = mean(F1s);

    save(name, 'model', 'F1s', 'meanF1', 'predictions_all', 'valid_labels', 'eval_ids', 'dataset_ids');

    % Write out the model
    name = sprintf('classifiers/AU_%d_dynamic.dat', au);

    pos_lbl = model.Label(1);
    neg_lbl = model.Label(2);

    w = model.w(1:end-1)';
    b = model.w(end);

    svs = bsxfun(@times, PC, 1./scaling') * w;

    write_lin_dyn_svm(name, means, svs, b, pos_lbl, neg_lbl);

    clear 'train_samples' 'train_labels' 'valid_samples' 'valid_labels'
end

author: TadasBaltrusaitis

repo: OpenFace


text:
%%
% This routine etsimates the b-value of a curve automatically
%  The b-valkue curve is differenciated and the point
%  of maximum curvature marked. The b-value will be calculated
%  using this point and the point half way toward the high
%  magnitude end of the b-value curve.
%
%  THIS IS RUN AUTOMATICALLY FROM MIN_OMORI1
%
%%

global cluscat mess bfig backcat xt3 bvalsum3  bval aw bw t1 t2 t3 t4;
global  ttcat n teb t0b cua b1 n1 b2 n2  ew si  S mrt bvalsumhold b;
global selt magco bvml avml bvls avls bv;
global hndl2 inpr1;
think

%report_this_filefun(mfilename('fullpath'));

maxmag = ceil(10*max(newt2.Magnitude))/10;
mima = min(newt2.Magnitude);
if mima > 0 ; mima = 0 ; end

% number of mag units
nmagu = (maxmag*10)+1;

bval = zeros(1,nmagu);
bvalsum = zeros(1,nmagu);
bvalsum3 = zeros(1,nmagu);

%%
%
% bval contains the number of events in each bin
% bvalsum is the cum. sum in each bin
% bval2 is number events in each bin, in reverse order
% bvalsum3 is reverse order cum. sum.
% xt3 is the step in magnitude for the bins == .1
%
%%

[bval,xt2] = hist(newt2.Magnitude,(mima:0.1:maxmag));
bvalsum = cumsum(bval); % N for M <=
bval2 = bval(length(bval):-1:1);
bvalsum3 = cumsum(bval(length(bval):-1:1));    % N for M >= (counted backwards)
xt3 = (maxmag:-0.1:mima);

backg_ab = log10(bvalsum3);

%%
% Estimate the b value -- based on one of 5 methods
%
% calculates max likelihood b value(bvml) && WLS(bvls)
%
%%

Nmin = 10;
bvs=newt2;
b=newt2;

% set DEFAULT to run using best of option
inpr1 = 5;

%% enough events??
if length(bvs) >= Nmin

    %%
    % calculation based on best combination of 90% and 95% probability -- default
    %%

    %            if inpr1 == 5
    mcperc_ca3;
    if isnan(Mc95) == 0
        magco = Mc95;
    elseif isnan(Mc90) == 0
        magco = Mc90;
    else
        [bv magco stan av me mer me2,  pr] =  bvalca3(bvs,1,1);
    end
    l = bvs(:,6) >= magco-0.05;
    if length(bvs(l)) >= Nmin
        %              [bvls magco0 stanls avls me mer me2,  pr] =  bvalca3(bvs(l,:),2,2);
        [mea bvml stanml avml ] =  bmemag(b(l,:));
    else
        bv = nan; bv2 = nan; magco = nan; av = nan; av2 = nan;
    end

end

author: CelsoReyes


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Eugenio Alcala Baselga
% Date: 02/06/2018
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

classdef automatic_dynamic_control
    %...
    properties (Constant)
        V_vec       = [2 18];
        Steer_vec   = [deg2rad(-25) deg2rad(25)];
        Alpha_vec   = [-0.1 0.1];

        sigma_vec   = [deg2rad(-25)+deg2rad(30) deg2rad(25)+deg2rad(30)];

        steer_min   = -deg2rad(25);
        steer_max   = deg2rad(25);
        sigma_min   = -deg2rad(25) + deg2rad(30);
        sigma_max   = deg2rad(25) + deg2rad(30);

%         Force_max   = 3500;
%         Force_min   = 0;

        Ts          = 0.001;
        gamma       = 0.001;

        Q       = [ 0.05 0 0 0 0 0 0;         %v
                    0 0.01 0 0 0 0 0;         %alpha
                    0 0 0.01 0 0 0 0;         %w
                    0 0 0 0.01 0 0 0;         %Fxr
                    0 0 0 0 10000 0 0;       %delta
                    0 0 0 0 0 90000 0;        %w integral
                    0 0 0 0 0    0 0];       %v integral
        R       = [0.01  0;
                   0    10];

    end

    methods
    end

end


text:
% Demo for DCM for fMRI (distributed responses)
% This demo inverts the DCM for fMRI model, which contains the ballon model
% as a generalized observation function (not affected by stochastic
% innovations).

close all
clear variables

%-----------------------------------------------------------
%-------------- DCM model specification --------------------

%--- Basic settings
TR = 1e0;                      % sampling period (in sec)
n_t = round(1e2/TR);            % number of time samples
dtU = round(2/TR)+1;           % input-on time interval
t0U = round(10/TR)+1;
microDT = 1e-1;               % micro-time resolution (in sec)
homogeneous = 1;              % params of g(x) homogeneous accross regions
reduced_f = 1;
alpha   = 2e2/TR;              % state noise precision
sigma   = 1e0;              % measurement noise precision

%--- Input
u       = zeros(2,n_t);
u(1,t0U:t0U+dtU) = 1;
u(1,4*t0U:4*t0U+dtU) = 1;
u(2,4*t0U:4*t0U+5*dtU) = 1;

%--- DCM structure
% invariant effective connectivity
A = [0 1 1
     1 0 1
     0 1 0];
nreg = size(A,1);
% modulatory effects
B{1} = zeros(nreg,nreg);
B{2} = [0 0 0
        1 0 0
        0 0 0];
% input-state coupling
C = [1 0
     0 0
     0 0];
% gating (nonlinear) effects
D{1} = [0 0 0
        0 0 0
        0 1 0];
D{2} = zeros(nreg,nreg);
D{3} = zeros(nreg,nreg);

%--- Build options/dim structures for model inversion
disp('Extracting HRF parameters...')
[thetaHRF,phiHRF]   = get_HRFparams(TR,microDT,homogeneous);
disp('Done.')
thetaHRF = 0.*thetaHRF;
phiHRF = 0.*phiHRF;
f_fname = @f_DCMwHRF;
g_fname = @g_HRF_distributed;
[options] = prepare_fullDCM(A,B,C,D,TR,microDT);
options.inG.n_phi = 4; % number of spatial modes
options.inG.B = VBA_sigmoid(randn(8,options.inG.n_phi)); % spatial modes
options.inG.ind_hrf = 1:2*nreg;
options.inG.n_reg = nreg;
for i=1:nreg
    options.inG.ind3{i} = options.inG.ind_hrf(end)+1+(i-1)*4:...
        options.inG.ind_hrf(end)+i*4;
end
options.MaxIter = 8; % to quicken inversion
dim.n_theta         = options.inF.ind5(end);
dim.n_phi           = (2+options.inG.n_phi).*nreg;
dim.n               = 5*nreg;
dim.p               = 8;
options.dim=dim;

%--- Build priors for model inversion
indHemo = options.inF.indself+1:dim.n_theta;
priors.muX0 = kron(ones(nreg,1),[0;0;0;0;0]);
priors.SigmaX0 = 0e-1*speye(5*nreg);
priors.muTheta = 1e-1*ones(dim.n_theta,1);
priors.muTheta(options.inF.indself) = -0;
priors.SigmaTheta = 1e-0*eye(dim.n_theta);
priors.SigmaTheta(options.inF.indself,options.inF.indself) = 0;
if reduced_f
    % fix some HRF params to their default values
    priors.muTheta(options.inF.ind1) = thetaHRF(1);
    priors.muTheta(options.inF.ind3) = thetaHRF(3);
    priors.muTheta(options.inF.ind5) = thetaHRF(6);
    priors.SigmaTheta(options.inF.ind1,options.inF.ind1) = 0;
    priors.SigmaTheta(options.inF.ind3,options.inF.ind3) = 0;
    priors.SigmaTheta(options.inF.ind5,options.inF.ind5) = 0;
end
priors.muPhi = 1e-1*ones(dim.n_phi,1);
priors.SigmaPhi = 1e-0*eye(dim.n_phi);
priors.SigmaPhi(options.inG.ind_hrf,options.inG.ind_hrf) = 0;

% NB on hyperpriors:
%   - fix state noise precision using high scale param
%   - use non-informative priors on the residual precision, with high
%   expectation.
% This is because of the first iteration of the hidden states posterior
% update, which has to deviate from the its prior predictive density (as
% derived from the deterministic inversion).
% The following iterations will then work with a realistic (expected)
% residual precision, and adapt.

SC = 1e2;
priors.a_alpha = Inf;%SC*alpha;
priors.b_alpha = 0;%SC;
priors.a_sigma = 1e0;
priors.b_sigma = 1e-4;

% State noise time-dependent covariance structure.
% NB: ratio of neural versus hemodynamic precision
for t = 1:n_t
    dq = 1e4*ones(dim.n,1);
    dq(options.inF.n5) = 1;
    priors.iQx{t} = diag(dq);
end
options.priors = priors;
options.backwardLag = 10;
options.gradF = 0;
options.updateHP = 1;
options.init0 = 0;

%-----------------------------------------------------------
%----------- simulated times series specification ----------

%--- simulated evolution parameters: neuronal level
% A matrix
t_A = exp([ -0.5
            -0.5
            -0.5
            -2.5
            -1.5 ]);
% self-inhibition gain
t_Aself = -0;
% B matrices
t_B{1} = [];
t_B{2} = exp([ -0.5 ]);
% C matrix
t_C = exp([ +0.1 ]);
% D matrices
t_D{1} = -exp([ -1.5 ]);
t_D{2} = [];
t_D{3} = [];

%--- simulated evolution parameters: hemodynamic level
t_E0 = thetaHRF(1)*ones(nreg,1);       % HbO2 extraction fraction gain
t_tau0 = thetaHRF(2)*ones(nreg,1);     % mean blood transit time gain
t_kaf = thetaHRF(3)*ones(nreg,1);      % vasodilatory signal feedback regulation
t_kas = thetaHRF(4)*ones(nreg,1);      % vasodilatory signal decay gain
t_alpha = thetaHRF(6)*ones(nreg,1);    % vessel stifness gain

%--- simulated observation parameters
if ~homogeneous
    p_E0 = phiHRF(1)*ones(nreg,1);       % HbO2 extraction fraction gain
    p_epsilon = phiHRF(2)*ones(nreg,1);  % ratio of intra- and extravascular signal
else
    p_E0 = phiHRF(1);
    p_epsilon = phiHRF(2);
end
p_w = repmat([1;2;-1;-2],nreg,1);   % spatial pattern parameters

%--- Recollect paramters for simulated data
nu = size(u,1);
theta = zeros(dim.n_theta,1);
theta(options.inF.indA) = t_A;
for i=1:nu
    theta(options.inF.indB{i}) = t_B{i};
end
theta(options.inF.indC) = t_C;
for i=1:nreg
    theta(options.inF.indD{i}) = t_D{i};
end
theta(options.inF.indself) = t_Aself;
theta(options.inF.ind1) = t_E0;
theta(options.inF.ind2) = t_tau0;
theta(options.inF.ind3) = t_kaf;
theta(options.inF.ind4) = t_kas;
theta(options.inF.ind5) = t_alpha;
phi(options.inG.ind1) = p_E0;
phi(options.inG.ind2) = p_epsilon;
phi = [ phi(:) ; p_w ];

% Build time series of hidden states and observations
[y,x,x0,eta,e] = VBA_simulate (n_t,f_fname,g_fname,theta,phi,u,alpha,sigma,options);

% display time series of hidden states and observations
displaySimulations(y,x,eta,e);
% disp('--paused--')
% pause

%-----------------------------------------------------------
%------------------- model inversion -----------------------

% Call inversion routine
[posterior,out] = VBA_NLStateSpaceModel(y,u,f_fname,g_fname,dim,options);

% Display results
displayResults(posterior,out,y,x,x0,theta,phi,alpha,sigma);

% Make predictions
try
    options = out.options;
    [xs,ys,xhat,vx,yhat,vy] = VBA_comparePredictions(...
        n_t,theta,phi,u,alpha,sigma,options,posterior,dim);
catch
    disp('------!!Unable to form predictions!!------')
end

author: MBB-team

repo: VBA-toolbox


text:
function [bbs, labels] = acfDetectBadacost( I, detector, fileName)
% Run aggregate channel features object detector on given image(s).
%
% The input 'I' can either be a single image (or filename) or a cell array
% of images (or filenames). In the first case, the return is a set of bbs
% where each row has the format [x y w h score] and score is the confidence
% of detection. If the input is a cell array, the output is a cell array
% where each element is a set of bbs in the form above (in this case a
% parfor loop is used to speed execution). If 'fileName' is specified, the
% bbs are saved to a comma separated text file and the output is set to
% bbs=1. If saving detections for multiple images the output is stored in
% the format [imgId x y w h score] and imgId is a one-indexed image id.
%
% A cell of detectors trained with the same channels can be specified,
% detected bbs from each detector are concatenated. If using multiple
% detectors and opts.pNms.separate=1 then each bb has a sixth element
% bbType=j, where j is the j-th detector, see bbNms.m for details.
%
% USAGE
%  bbs = acfDetectBadacostTrees( I, detector, [fileName] )
%
% INPUTS
%  I          - input image(s) of filename(s) of input image(s)
%  detector   - detector(s) trained via acfTrainBadacostTrees
%  fileName   - [] target filename (if specified return is 1)
%
% OUTPUTS
%  bbs        - [nx5] array of bounding boxes or cell array of bbsc
%  labels     - [nx1] array of subclass labels (for example orientation
%                     class in car detection
%
% EXAMPLE
%
% See also acfTrain, acfModify, bbGt>loadAll, bbNms
%
% Piotr's Computer Vision Matlab Toolbox      Version 3.40
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

% run detector on every image
if(nargin<3), fileName=''; end; multiple=iscell(I);
if(nargin<4), aRatio=[]; end;
if(~isempty(fileName) && exist(fileName,'file')), bbs=1; return; end
if(~multiple), [bbs, labels]=acfDetectImgBadacost(I,detector); else
  n=length(I); bbs=cell(n,1); labels=cell(n,1);
  parfor i=1:n, [bbs{i}, labels{i}]=acfDetectImgBadacost(I{i},detector); disp(i); end
end

% write results to disk if fileName specified
if(isempty(fileName)), return; end
d=fileparts(fileName); if(~isempty(d)&&~exist(d,'dir')), mkdir(d); end
if( multiple ) % add image index to each bb and flatten result
  for i=1:n, bbs{i}=[ones(size(bbs{i},1),1)*i bbs{i}]; end
  are_empty = cellfun(@(x) ((size(x,1) == 0) || size(x,2)==0),bbs);
  bbs(are_empty) = [];
  are_empty = cellfun(@(x) ((size(x,1) == 0) || size(x,2)==0),labels);
  labels(are_empty) = [];
  bbs=cell2mat(bbs);
  labels=cell2mat(labels);
end
dlmwrite(fileName,bbs); bbs=1;

end

function [bbs, labels] = acfDetectImgBadacost( I, detector ) %, aRatio )
% Run trained sliding-window object detector on given image.
Ds=detector; if(~iscell(Ds)), Ds={Ds}; end; nDs=length(Ds);
opts=Ds{1}.opts; pPyramid=opts.pPyramid; pNms=opts.pNms;
imreadf=opts.imreadf; imreadp=opts.imreadp;
shrink=pPyramid.pChns.shrink; pad=pPyramid.pad;
separate=nDs>1 && isfield(pNms,'separate') && pNms.separate;
% read image and compute features (including optionally applying filters)
if(all(ischar(I))), I=feval(imreadf,I,imreadp{:}); end
P=chnsPyramid(I,pPyramid); bbs=cell(P.nScales,nDs);
if(isfield(opts,'filters') && ~isempty(opts.filters)), shrink=shrink*2;
  for i=1:P.nScales, fs=opts.filters; C=repmat(P.data{i},[1 1 size(fs,4)]);
    for j=1:size(C,3), C(:,:,j)=conv2(C(:,:,j),fs(:,:,j),'same'); end
    P.data{i}=imResample(C,.5);
  end
end
% apply sliding window classifiers
for i=1:P.nScales
  for j=1:nDs, opts=Ds{j}.opts;
    modelDsPad=opts.modelDsPad; modelDs=opts.modelDs;

    % -- Use (with changed sign) min positive class cost minus negative class cost as score
    % WARNING: Change to accomodate with theory. We transpose Cprime because acfDetectBadacostTrees1
    %          internally transposes the cost matrix that is passed to it and we want that the
    %          actual matrix used is Cprime, transpose(transpose(Cprime)) => Cprime (2016/11)
    clf_aux = Ds{j}.clf;
    clf_aux.Cprime = clf_aux.Cprime';
    [bb, labels] = acfDetectBadacostTrees1(P.data{i},clf_aux,shrink,...
        modelDsPad(1),modelDsPad(2),opts.stride,opts.cascThr);
    % End WARNING.

    shift=(modelDsPad-modelDs)/2-pad;
    bb(:,1)=(bb(:,1)+shift(2))/P.scaleshw(i,2);
    bb(:,2)=(bb(:,2)+shift(1))/P.scaleshw(i,1);
    bb(:,3)=modelDs(2)/P.scales(i);
    bb(:,4)=modelDs(1)/P.scales(i);
    bb(:,6)=labels(:);  % We make the label of the subclass the 6th argument.
    if(separate), bb(:,7)=j; end; bbs{i,j}=bb;
  end;
end; bbs=cat(1,bbs{:});

if(isfield(Ds{1}.clf, 'aRatio'))
  % JMBUENA: Change bounding boxes to class specific bounding box. For
  % example in the KITTI benchmark we have cars. At each orientation the
  % bounding box of the car has a specific aspect ratio (a size view of the
  % car is rectangular and a frontal car is squared).
  use_fixed_width = 0;
  if isfield(Ds{1}.clf, 'aRatioFixedWidth')
    use_fixed_width = Ds{1}.clf.aRatioFixedWidth;
  end
  bbs = correctToClassSpecificBbs(bbs, Ds{1}.clf.aRatio, use_fixed_width, 1);
end

if(~isempty(pNms)), bbs=bbNms(bbs,pNms); end
end

author: MengyangPu

repo: EDTER

save_path: github-repos/MATLAB/MengyangPu-EDTER


text:
function [edgeS, maskDown3D] = getSurface(structNumV, marginV, xyDownsampleIndex, planC)
%"getSurface"
%   Get all the surface points of a a composite structure defined by the
%   structures given in structNumV, with margins epanded in 3D according to
%   marginV.  The iV, jV, kV index vectors give all the surface points with
%   respect to the uniformized CT scan.
%
%JOD, 14 Nov 03.
%JRA, 27 Mar 05.
%
%Usage:
%   [iEdgeV,jEdgeV,kEdgeV] = getSurface(structNumV,marginV, xyDownsampleIndex)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

if ~exist('planC')
    global planC
end
indexS = planC{end};

global stateS

if any(marginV ~= marginV(1))
  error('Currently only supports the same PB margin around each target.')
end
margin = marginV(1);

scanNum = getStructureAssociatedScan(structNumV(1));
CTUniformInfoS = planC{indexS.scan}(scanNum).uniformScanInfo;

xOffset = CTUniformInfoS.xOffset;
yOffset = CTUniformInfoS.yOffset;

sliceThickness = CTUniformInfoS.sliceThickness;

delta_xy = CTUniformInfoS.grid1Units;

%-----------build composite target volume---------------------%

mask3D = getUniformStr(structNumV(1));

SZ=size(mask3D);

maskDown3D = logical(getDownsample3(mask3D, xyDownsampleIndex, 1));

clear mask3D

S = size(maskDown3D);

len = length(structNumV);

for i = 2 : len

  maskSingle = getUniformStr(structNumV(i));

  mask3D = maskSingle; clear maskSingle;

  tmp3D = logical(getDownsample3(mask3D, xyDownsampleIndex, 1));

  clear mask3D

  maskDown3D = tmp3D | maskDown3D;

end

% II.  Get surface points from target volume

edge = [0 0 0; 0 1 0; 0 0 0];
edge(:,:,2) = [0 1 0; 1 1 1; 0 1 0];
edge(:,:,3) = [0 0 0; 0 1 0; 0 0 0];
edge = edge/7;

surfPoints = getSurfacePoints(maskDown3D);
edge3D = repmat(logical(0), size(maskDown3D));
for i=1:size(surfPoints,1)
   edge3D(surfPoints(i,1),surfPoints(i,2), surfPoints(i,3)) = 1;
end

%Expand margin using convolution
%Create margin ball:

c1 = ceil(margin/delta_xy);
c2 = ceil(margin/delta_xy);
c3 = ceil(margin/sliceThickness);

[uM,vM,wM] = meshgrid(- c1 : c1, -c2 : c2, - c3 : c3);

xM = uM * delta_xy;
yM = vM * delta_xy;
zM = wM * sliceThickness;

rM = (xM.^2 + yM.^2 + zM.^2).^0.5;

ball = [rM <= margin];

[iBallV,jBallV,kBallV] = find3d(ball);

sR = size(rM);

deltaV = (sR - 1)/2 +1;

onesV = repmat(logical(1), [1,length(iBallV)]);

[iV,jV,kV] = find3d(edge3D);

sV = size(maskDown3D);

ind_surfV = sub2ind(sV,iV,jV,kV);

ball_offsetV = (iBallV - deltaV(1)) + sV(1) * (jBallV - deltaV(2)) + sV(1) * sV(2) * (kBallV - deltaV(3));

for i = 1 : length(ind_surfV) %put ones in

  total_indV = ind_surfV(i) + ball_offsetV;

  total_indV = clip(total_indV,1,prod(sV),'limits');

  maskDown3D(total_indV) = onesV;

end

%Find new edge
surfPoints = getSurfacePoints(maskDown3D);

edgeS.rows   = surfPoints(:,1);%iV;
edgeS.cols   = surfPoints(:,2);%jV;
edgeS.slices = surfPoints(:,3);%kV;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function varargout = process_ft_mtmconvol( varargin )
% PROCESS_FT_MTMCONVOL: Call FieldTrip function ft_mtmconvol.
%
% REFERENCES:
%     - http://www.fieldtriptoolbox.org/tutorial/timefrequencyanalysis

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'FieldTrip: ft_mtmconvol (Multitaper)';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Frequency';
    sProcess.Index       = 506;
    sProcess.Description = 'http://www.fieldtriptoolbox.org/tutorial/timefrequencyanalysis';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'import', 'data'};
    sProcess.OutputTypes = {'import', 'data'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'results', 'matrix'};
    sProcess.OutputTypes = {'timefreq', 'timefreq', 'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Options: Time window
    sProcess.options.timewindow.Comment = 'Time window:';
    sProcess.options.timewindow.Type    = 'timewindow';
    sProcess.options.timewindow.Value   = [];
    sProcess.options.timewindow.Group   = 'input';
    % Options: Sensor types
    sProcess.options.sensortypes.Comment = 'Sensor types or names (empty=all): ';
    sProcess.options.sensortypes.Type    = 'text';
    sProcess.options.sensortypes.Value   = 'MEG, EEG';
    sProcess.options.sensortypes.InputTypes = {'data'};
    sProcess.options.sensortypes.Group   = 'input';
    % Options: Scouts
    sProcess.options.clusters.Comment = '';
    sProcess.options.clusters.Type    = 'scout_confirm';
    sProcess.options.clusters.Value   = {};
    sProcess.options.clusters.InputTypes = {'results'};
    sProcess.options.clusters.Group   = 'input';
    % Options: Scout function
    sProcess.options.scoutfunc.Comment    = {'Mean', 'Max', 'PCA', 'Std', 'All', 'Scout function:'};
    sProcess.options.scoutfunc.Type       = 'radio_line';
    sProcess.options.scoutfunc.Value      = 1;
    sProcess.options.scoutfunc.InputTypes = {'results'};
    sProcess.options.scoutfunc.Group   = 'input';

    % Options: Taper
    sProcess.options.mt_taper.Comment = 'Taper: ';
    sProcess.options.mt_taper.Type    = 'combobox_label';
    sProcess.options.mt_taper.Value   = {'dpss', {'dpss', 'hanning', 'rectwin', 'sine'; ...
                                                  'dpss', 'hanning', 'rectwin', 'sine'}};
    % Options: Frequencies
    sProcess.options.mt_frequencies.Comment = 'Frequencies (start:step:stop): ';
    sProcess.options.mt_frequencies.Type    = 'text';
    sProcess.options.mt_frequencies.Value   = '1:2:120';
    % Options: Frequency resolution
    sProcess.options.mt_freqmod.Comment = 'Modulation factor: ';
    sProcess.options.mt_freqmod.Type    = 'value';
    sProcess.options.mt_freqmod.Value   = {10, ' ', 0};
    % Options: Time resolution
    sProcess.options.mt_timeres.Comment = 'Time resolution: ';
    sProcess.options.mt_timeres.Type    = 'value';
    sProcess.options.mt_timeres.Value   = {1, 'ms', []};
    % Options: Time step
    sProcess.options.mt_timestep.Comment = 'Time step: ';
    sProcess.options.mt_timestep.Type    = 'value';
    sProcess.options.mt_timestep.Value   = {0.1, 'ms', []};
    % === CORRESPONDANCE WITH FIELDTRIP
    sProcess.options.mt_label.Comment = ['
' ...
        'Correspondance with FieldTrip inputs:
' ...
        ' - timeoi = tw(1)+tr/2 : time_step : tw(2)-tr/2-1/sfreq
' ...
        '        tr=time_resolution, tw=time_window
' ...
        ' - tapsmofrq = frequencies / modulation_factor
' ...
        ' - timwin = repmat(time_resolution, 1, length(frequencies))'];
    sProcess.options.mt_label.Type    = 'label';

    % === MEASURE
    sProcess.options.measure.Comment = {'Power', 'Magnitude', 'Measure: '; ...
                                        'power', 'magnitude', ''};
    sProcess.options.measure.Type    = 'radio_linelabel';
    sProcess.options.measure.Value   = 'power';
    sProcess.options.measure.Group   = 'output';
    % === AVERAGE OUTPUT FILES
    sProcess.options.avgoutput.Comment = 'Save average power across files
(do not save one new file per input file)';
    sProcess.options.avgoutput.Type    = 'checkbox';
    sProcess.options.avgoutput.Value   = 1;
    sProcess.options.avgoutput.Group   = 'output';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    OutputFiles = {};
    % Initialize FieldTrip or SPM: ft_specest_mtmconvol is in both
    [isInstalled, errMsg] = bst_plugin('InstallMultipleChoice', {'fieldtrip', 'spm12'});
    if ~isInstalled
        bst_report('Error', sProcess, [], errMsg);
        return;
    end
    bst_plugin('SetProgressLogo', 'fieldtrip');
    % Check for Signal Processing Toolbox when using DPSS
    if iscell(sProcess.options.mt_taper.Value)
        mt_taper = sProcess.options.mt_taper.Value{1};
    else
        mt_taper = sProcess.options.mt_taper.Value;
    end
    if strcmpi(mt_taper, 'dpss') && ~exist('dpss', 'file')
        bst_report('Error', sProcess, [], 'The option "dpss" requires the Signal Processing Toolbox.');
        return;
    end
    % Call TIME-FREQ process
    OutputFiles = process_timefreq('Run', sProcess, sInputs);
    % Remove logo
    bst_plugin('SetProgressLogo', []);
end

author: brainstorm-tools


text:
function [h, ax]= visutil_addTitle(tit, VERTICAL, shift)
%[h, ax]= addTitle(tit, <vert=0, shift=1>);

if ~exist('shift', 'var'), shift=1; end
if ~exist('VERTICAL', 'var'),
  VERTICAL= 0;
elseif ischar(VERTICAL),
  VERTICAL= strmatch(VERTICAL, {'horizontal','vertical'})-1;
end

if VERTICAL,
  borderSize= 0.06;
  posShrink= [1-borderSize 1 1-borderSize 1];
  posShift= [1 0 0 0];
else
  borderSize= 0.1;
  posShrink= [1 1-borderSize 1 1-borderSize];
  posShift= [0 0 0 0];
end

if shift,
  hc= get(gcf, 'children');
  for hi= 1:length(hc),
    if strcmp(get(hc(hi), 'type'), 'axes'),
      pos= get(hc(hi), 'position');
      set(hc(hi), 'position', posShift - (posShift-pos) .* posShrink);
    end
  end
end

ax= visutil_getBackgroundAxis;
prop= {'visible', 'on', ...
       'horizontalAlignment', 'center', ...
       'verticalAlignment', 'middle'};
if VERTICAL,
  prop= {prop{:}, 'rotation', 90};
  h= visutil_textFit(borderSize/2, 0.5, tit, [borderSize 0.9], prop{:});
else
  h= visutil_textFit(0.5, 1-borderSize/2, tit, [0.8 borderSize], prop{:});
end

author: bbci


text:
function quadrule_test39 ( )

%*****************************************************************************80
%
%% QUADRULE_TEST39 tests NCO_SET and SUM_SUB.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2009
%
%  Author:
%
%    John Burkardt
%
  order_max = 9;

  nfunc = func_set ( 'COUNT', 'DUMMY' );

  a = 0.0;
  b = 1.0;

  nsub = 1;
  xlo = -1.0;
  xhi = +1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUADRULE_TEST39\n' );
  fprintf ( 1, '  NCO_SET sets up an open Newton-Cotes rule;\n' );
  fprintf ( 1, '  SUM_SUB carries it out.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Integration interval is [%f, %f]\n', a, b );
  fprintf ( 1, '  Number of subintervals is %d\n', nsub );
  fprintf ( 1, '  Quadrature order will vary.\n' );
  fprintf ( 1, '  Integrand will vary.\n' );
  fprintf ( 1, '\n' );

  for ilo = 1 : 5 : nfunc

    ihi = min ( ilo + 4, nfunc );

    fprintf ( 1, '\n' );
    fprintf ( 1, '    ' );
    for i = ilo : ihi
      fprintf ( '%14s', fname(i) );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '\n' );

    for norder = 1 : order_max

      if ( norder == 8 )
        continue
      end

      fprintf ( 1, '  %2d', norder );

      for i = ilo : ihi

        func_set ( 'SET', i );

        [ xtab, weight ] = nco_set ( norder );

        result(i) = sum_sub ( @func, a, b, nsub, norder, xlo, xhi, xtab, ...
          weight );

        fprintf ( 1, '  %12f', result(i) );

      end

      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
function scores = ras(pop)
% shortcut to rastriginsfcn()

scores = 10.0 * size(pop,2) + sum(pop .^2 - 10.0 * cos(2 * pi .* pop),2);


text:
function [centers,mincenter,mindist,q2,quality] = kmeansFast(data,initcenters,method)
% output: final centers
% input: data points and initial centers
% if initcenters is a number k, create k centers and start with these
% otherwise, use centers given as input
% method = 0: unoptimized, using n by k matrix of distances O(nk) space
%          1: vectorized, using only O(n+k) space
%          2: like 1, in addition using distance inequalities (default)
% URL: http://cseweb.ucsd.edu/~elkan/fastkmeans.html

tic
if nargin < 3 method = 2; end
[n,dim] = size(data);

if max(size(initcenters)) == 1
    k = initcenters;
    [centers, mincenter, mindist, lower, computed] = anchors(mean(data),k,data);
    total = computed;
    skipestep = 1;
else
    centers = initcenters;
    mincenter = zeros(n,1);
    total = 0;
    skipestep = 0;
    [k,dim2] = size(centers);
    if dim ~= dim2 error('dim(data) ~= dim(centers)'); end;
end

nchanged = n;
iteration = 0;
oldmincenter = zeros(n,1);

while nchanged > 0
    % do one E step, then one M step
    computed = 0;

    if method == 0 && ~skipestep
        for i = 1:n
            for j = 1:k
                distmat(i,j) = calcdist(data(i,:),centers(j,:));
            end
        end
        [mindist,mincenter] = min(distmat,[],2);
        computed = k*n;

    elseif (method == 1 || (method == 2 && iteration == 0)) && ~skipestep
        mindist = Inf*ones(n,1);
        lower = zeros(n,k);
        for j = 1:k
           jdist = calcdist(data,centers(j,:));
           lower(:,j) = jdist;
           track = find(jdist < mindist);
           mindist(track) = jdist(track);
           mincenter(track) = j;
        end
        computed = k*n;

    elseif method == 2 && ~skipestep
        computed = 0;

% for each center, nndist is half the distance to the nearest center
% if d(x,center) < nndist then x cannot belong to any other center
% mindist is an upper bound on the distance of each point to its nearest center

        nndist = min(centdist,[],2);
% the following usually is not faster
%        ldist = min(lower,[],2);
%        mobile = find(mindist > max(nndist(mincenter),ldist));
        mobile = find(mindist > nndist(mincenter));

% recompute distances for point i and center j
%       only if j can possibly be the new nearest center
% for speed, the first check has been optimized by modifying centdist
% swapping the order of the checks is slower for data with natural clusters

        mdm = mindist(mobile);
        mcm = mincenter(mobile);

        for j = 1:k
% the following is incorrect: for j = unique(mcm)'
            track = find(mdm > centdist(mcm,j));
            if isempty(track) continue; end
            alt = find(mdm(track) > lower(mobile(track),j));
            if isempty(alt) continue; end
            track1 = mobile(track(alt));

% calculate exact distances to the mincenter
% recalculate separately for each jj to avoid copying too much of data
% redo may be empty, but we don't need to check this
            redo = find(~recalculated(track1));
            redo = track1(redo);
            c = mincenter(redo);
            computed = computed + size(redo,1);
            for jj = unique(c)'
                rp = redo(find(c == jj));
                udist = calcdist(data(rp,:),centers(jj,:));
                lower(rp,jj) = udist;
                mindist(rp) = udist;
            end
            recalculated(redo) = 1;

            track2 = find(mindist(track1) > centdist(mincenter(track1),j));
            track1 = track1(track2);
            if isempty(track1) continue; end

            % calculate exact distances to center j
            track4 = find(lower(track1,j) < mindist(track1));
            if isempty(track4) continue; end
            track5 = track1(track4);
            jdist = calcdist(data(track5,:),centers(j,:));
            computed = computed + size(track5,1);
            lower(track5,j) = jdist;

            % find which points really are assigned to center j
            track2 = find(jdist < mindist(track5));
            track3 = track5(track2);
            mindist(track3) = jdist(track2);
            mincenter(track3) = j;
        end % for j=1:k
    end % if method

    oldcenters = centers;

% M step: recalculate the means for each cluster
% if a cluster is empty, its mean is left unchanged
% we minimize computations for clusters with little changed membership

    diff = find(mincenter ~= oldmincenter);
    diffj = unique([mincenter(diff);oldmincenter(diff)])';
    diffj = diffj(diffj > 0);

    if size(diff,1) < n/3 && iteration > 0
         for j = diffj
            plus = find(mincenter(diff) == j);
            minus = find(oldmincenter(diff) == j);
            oldpop = pop(j);
            pop(j) = pop(j) + size(plus,1) - size(minus,1);
            if pop(j) == 0 continue; end
            centers(j,:) = (centers(j,:)*oldpop + sum(data(diff(plus),:),1) - sum(data(diff(minus),:),1))/pop(j);
        end
    else
        for j = diffj
            track = find(mincenter == j);
            pop(j) = size(track,1);
            if pop(j) == 0 continue; end
% it's correct to have mean(data(track,:),1) but this can make answer worse!
            centers(j,:) = mean(data(track,:),1);
        end
    end

    if method == 2
        for j = diffj
            offset = calcdist(centers(j,:),oldcenters(j,:));
            computed = computed + 1;
            if offset == 0 continue; end
            track = find(mincenter == j);
            mindist(track) = mindist(track) + offset;
            lower(:,j) = max(lower(:,j) - offset,0);
        end

% compute distance between each pair of centers
% modify centdist to make "find" using it faster
        recalculated = zeros(n,1);
        realdist = alldist(centers);
        centdist = 0.5*realdist + diag(Inf*ones(k,1));
        computed = computed + k + k*(k-1)/2;
    end

    nchanged = size(diff,1) + skipestep;
    iteration = iteration+1;
    skipestep = 0;
    oldmincenter = mincenter;

%   difference = max(max(abs(oldcenters - centers)));
%   [iteration toc nchanged computed size(diffj,2)]
    %[iteration toc nchanged computed]
    total = total + computed;
end % while nchanged > 0

udist = calcdist(data,centers(mincenter,:));
quality = mean(udist);
q2 = mean(udist.^2);
%[iteration toc quality q2 total]

author: adikhosla

repo: feature-extraction


text:
function audio = mergeAudio(data, window)

pad = [];
if window(1) <= 1
    pad=zeros(abs(window(1)), 1);
    window(1) = 1;
end
if isstring(data) |ischar(data)
    audio = audioread(data, window);
else
   audio = data;
end
audio = [mean(audio - mean(audio,1) ,2)]; % Take the mean of the audio channels
audio = int16(audio * 32767); % Convert to int16

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak