text:
%% Show results for system identification [training, validation]

clear ph
figure,box on, hold on,
ccolors = get(gca,'colororder');
plot([tA(1),tA(1)],[min(xB(:)) max([xA(:);xB(:)])],':','Color',[0.4,0.4,0.4],'LineWidth',1.5)
ylim([min(xB(:)) max([xA(:);xB(:)])])
plot(t,x(:,1),'Color',ccolors(1,:),'LineWidth',1);
plot(t,x(:,2),'Color',ccolors(2,:),'LineWidth',1);
ph(1) = plot([t;tA],[x(:,1);xA(:,1)],'Color',ccolors(1,:),'LineWidth',1);
ph(2) = plot([t;tA],[x(:,2);xA(:,2)],'Color',ccolors(2,:),'LineWidth',1);
ph(3) = plot([t;tA],[x(:,3);xA(:,3)],'Color',ccolors(3,:),'LineWidth',1);
ph(4) = plot(tB,xB(:,1),'-.','Color',ccolors(1,:)-[0 0.2 0.2],'LineWidth',2);
ph(5) = plot(tB,xB(:,2),'-.','Color',ccolors(2,:)-[0.1 0.2 0.09],'LineWidth',2);
ph(6) = plot(tB,xB(:,3),'-.','Color',ccolors(3,:)-[0.1 0.2 0.09],'LineWidth',2);
grid off
xlim([0 tB(end)])
xlabel('Time')
ylabel('xi')
t1 = text(tA(1)-5,40,'Training', 'FontSize',12);
t2 = text(1+tA(1),40,'Validation', 'FontSize',12);
lh = legend(ph([1,4]),'True',ModelName);
lh.Position = [lh.Position(1)-0.5,lh.Position(2)-0.15,lh.Position(3:4)];
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
print('-depsc2', '-loose','-cmyk', [figpath,'EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_Validation_OneFig.eps']);

delete(lh); delete(t1), delete(t2);
print('-depsc2', '-loose','-cmyk', [figpath,'EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_Validation_OneFig_noleg.eps']);

%% Actuation signal
clear ph
figure,box on, hold on,
ccolors = get(gca,'colororder');
plot([tA(1),tA(1)],[-15 260],':','Color',[0.4,0.4,0.4],'LineWidth',1.5)
plot(t,u,'-k','LineWidth',1);
plot(tv,uv,'-k','LineWidth',1);
grid off
ylim([min([u uv])+0.05*min([u uv]) max([u uv])+0.05*max([u uv])])
xlim([0 tB(end)])
xlabel('Time')
ylabel('Input')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
print('-depsc2', '-loose','-cmyk', [figpath,'EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_Actuation_OneFig.eps']);

author: eurika-kaiser

repo: SINDY-MPC


text:
% By Sherif A. Tawfik, Faculty of Engineering, Cairo University
function y=findzero(fun,x0,x1,varargin)
% fun is the function that we need to compute its root.
% x0 and x1 are two arguments to the function such that the root that we
% seek lie between them and the function has different sign at these two
% points
% varargin are the other arguments of the function

% the value of the function at the first point
f0=feval(fun,x0,varargin{:});

%the value of the function at the second point
f1=feval(fun,x1,varargin{:});

% check that the sign of the function at x0 and x1 is different. Otherwise
% report an error
if sign(f0)==sign(f1)
    error('the function at the two endpoints must be of opposite signs');
end

%find a closer point to the root using the method of chords. In the method
%of chords we simply connect the two points (x0, f(x0))  and (x1, f(x1))
%using a straight line and compute the intersection point with this line
%and the horizontal axis. This new point is closer to the desired root
x=x0 - f0 * ((x1-x0)/(f1-f0));

%evaluate the function at this new point
f=feval(fun,x,varargin{:});

% enter this root as long as the difference between the two points that
% sandwitch the desired root is larger than a certain threshold
while abs(f)>2^-52
    % we keep one of the two old points that has a different sign than the
    % new point and we overwrite the other old point with the new point
    if sign(f)==sign(f0)
        x0=x;
        f0=f;
    else
        x1=x;
        f1=f;
    end
    x=x0 - f0 * ((x1-x0)/(f1-f0));
    f=feval(fun,x,varargin{:});
end
% at the end of the loop we reach the root with the desired precision and
% it is given by x
y=x;

author: Sable

repo: mcbench-benchmarks


text:
function [uOutput] = import_harvard(nFunction, sFilename)
% [uOutput] = harvardcmt_imp(nFunction, sFilename);
% ----------------------------------------------------------
% Imports Harvard CMT data of a selfmade ASCII-file.
% Use 'getharvardcatalog' in the importfilters/harvard directory.
% See import_harvard_doc.html for more information.
%
% D. Schorlemmer; schorlemmer@sed.ethz.ch
%
% 17.05.2005

% Filter function switchyard
if nFunction == FilterOp.getDescription
  uOutput = 'Harvard CMT catalog';
elseif nFunction == FilterOp.getWebpage
  uOutput = 'import_harvard_doc.html';
elseif nFunction == FilterOp.importCatalog
  % Read formated data
  mData = textread(sFilename, '%s', 'delimiter', '\n', 'whitespace', '');
  % Create empty catalog
  uOutput = zeros(length(mData), 10);
  % Loop thru all lines of catalog and convert them
  for i = 1:length(mData)
    if rem(i,100) == 0 ; disp([ num2str(i) ' of ' num2str(length(mData)) ' events processed']); end
    try
      uOutput(i,1) = str2num(mData{i}(28:34));      % Longitude (PDE)
      uOutput(i,2) = str2num(mData{i}(21:25));      % Latitude  (PDE)
      uOutput(i,3) = str2num(mData{i}(1:2));        % Year
      if uOutput(i,3) < 76
        uOutput(i,3) = uOutput(i,3)+2000;
      else
        uOutput(i,3) = uOutput(i,3)+1900;
      end
      uOutput(i,4) = str2num(mData{i}(4:5));        % Month
      uOutput(i,5) = str2num(mData{i}(7:8));        % Day
      uOutput(i,6) = str2num(mData{i}(82:85));      % Magnitude Mw
      uOutput(i,7) = str2num(mData{i}(36:40));      % Depth (PDE)
      uOutput(i,8) = str2num(mData{i}(10:11));      % Hour
      uOutput(i,9) = str2num(mData{i}(13:14));      % Minute
      uOutput(i,10) = str2num(mData{i}(16:19));     % Second
      uOutput(i,11) = nan;                          % Reserved for cross-section values

      uOutput(i,12) = str2num(mData{i}(87:89));     % Strike (Plane 1)
      uOutput(i,13) = uOutput(i,12) + 90;           % Dip direction (Plane 1)
      uOutput(i,13) = mod(uOutput(i,13) + 360, 360);
      uOutput(i,14) = str2num(mData{i}(91:92));     % Dip (Plane 1)
      uOutput(i,15) = str2num(mData{i}(94:97));     % Rake (Plane 1)

      uOutput(i,16) = str2num(mData{i}(99:101));    % Strike (Plane 2)
      uOutput(i,17) = uOutput(i,16) + 90;           % Dip direction (Plane 2)
      uOutput(i,17) = mod(uOutput(i,17) + 360, 360);
      uOutput(i,18) = str2num(mData{i}(103:104));   % Dip (Plane 2)
      uOutput(i,19) = str2num(mData{i}(106:109));   % Rake (Plane 2)

      uOutput(i,20) = str2num(mData{i}(42:44));     % Magnitude mb
      uOutput(i,21) = str2num(mData{i}(46:48));     % Magnitude Ms
      uOutput(i,22) = str2num(mData{i}(50:55));     % Latitude (HAV)
      uOutput(i,23) = str2num(mData{i}(57:63));     % Longitude (HAV)
      uOutput(i,24) = str2num(mData{i}(65:69));     % Depth (HAV)
      uOutput(i,25) = str2num(mData{i}(72:75));     % Cen_time
      uOutput(i,26) = str2num(mData{i}(77:80));     % Half Duration
      % Create decimal year
      uOutput(i,3) = decyear([uOutput(i,3) uOutput(i,4) uOutput(i,5) uOutput(i,8) uOutput(i,9) uOutput(i,16)]);
    catch
      msg.dbfprintf('Import: Problem in line %d of %s. Line ignored.\n',i, sFilename);
    end
  end
end

author: CelsoReyes


text:
function [x,F,JAC,EXITFLAG] = NewtonRaphson(FUN,x,lambda, maxIter ,Display)
%NewtonRaphson solves a system of nonlinear equations via newton method
%
%   NewtonRaphson solves equations of the form:
%
%   F(X) = 0    where F and X may be scalars or vectors
%
%  NewtonRaphson implements the damped newton method with adaptive step
%  size. Theory and discussion can be found here:
%  http://forum.vis.ethz.ch/showthread.php?13434-Damped-Newton-method-Schrittweitensteuerung
%
%  The Jacobian is calculated numerically by a simple forward differential
%  method. Find more here:
%  http://en.wikipedia.org/wiki/Numerical_differentiation
%
%  The error estimation ist made by root-mean-square
%  http://en.wikipedia.org/wiki/Root_mean_square
%
%
%   x = NewtonRaphson(FUN,X0) starts at the initial guess X0 and tries to
%    solve the equations in FUN. FUN is a function handle and has to accept
%    input x and return a vector of equation values F evaluated at x.
%    Default values for solver and display setting.
%
%   x = NewtonRaphson(FUN,X0,lambda) starts at the initial guess X0 and tries to
%    solve the equations in FUN with user supplied initial relaxation factor. Might
%    be useful to increase solution speed.
%   default value: lambda = 0.1
%
%   x = NewtonRaphson(FUN,X0,lambda,maxIter) ...
%    User supplied maximum number of iterations
%    default value: maxIter = 100
%
%   x = NewtonRaphson(FUN,X0,lambda,maxIter,Display) ...
%    Display options:
%         'on' - full output during solution process
%         'off' - hide output
%    default value: maxIter = 'off'
%
%  OUTPUT Arguments:
%   x - Solution
%   F - Value at x
%   Jac - Jacobian at x
%   Exitflat: exit conditions of damped newton
%     1: all ok. Solution found
%     -1: no solution found
%     -2: FUN is not a function handle
%
%   To enter Demo Mode start without any arguments
%     NewtonRaphson();
%
%   Examples:
%
%  Find zero of function atan(x)
%
%   create function handle (here as an anonymous function)
%  F = @(x)atan(x);
%  x0 = 5; % start value /initial guess
%
%  x = NewtonRaphson(F,x0,1);
%
%  Find solution of following system of equations (same as in fsolve help)
%  F = @(x)[2*x(1) - x(2) - exp(-x(1));
%      -x(1) + 2*x(2) - exp(-x(2))];
%  x0 = [1;2];
%
%  x = NewtonRaphson(F,x0,0.1,100,'on');
%
%  Open Issues:
%   - At least some input error checks have to be done... I'm so lazy...
%   - correct some wording/grammar/spelling
%
% Author: Andi S. 2013

% Check Input arguments / find more at help nargin
if nargin < 5, Display = 'off'; end
if nargin < 4, maxIter = 100; end
if nargin < 3, lambda  = .1; end

if(nargin==0)
    disp('Demo-Mode on f(x) = atan(x)');
    FUN = @(x)atan(x);
    x = 5; % initial guess
    Display = 'on';
    % maximum number of iterations
    maxIter = 100;
    % Initial relaxation factor
    lambda = 1.;

end

% check for Input Arguments error
if(~isa(FUN,'function_handle'))
    disp('FUN must be a function handle.');
    x = 0; F = 0; JAC = 0; EXITFLAG = -2;
    return
end

% -------------------------------------------------------------------------
% Solver settings
% -------------------------------------------------------------------------
% Estimate problem size
N = length(x);

% Minimal relaxation factor
Lmin = 1e-10;

% Convergence Criteria
ConvCrit = 1e-6;

% Allocate Memory for Jacobian
JAC = zeros(N);

% -------------------------------------------------------------------------
% Solution
% -------------------------------------------------------------------------

% First Newton Step
% get first Values for x
F = FUN(x);

% numerical Jacobian
for i = 1:N
    xtemp = x;
    h = sqrt(eps)*x(i);
    xtemp(i) = xtemp(i)+h;
    JAC(i,:) = (FUN(xtemp)-F)./h;
end

% calculate step size
s = (JAC\-F);

% first newton step
x = x + lambda*s;
F = FUN(x);
err = sqrt(sum(F.^2)/N);

n = 0; % Iteration counter

if(strcmp(Display,'on'));
    fprintf('n: %i rlx fac: %e error: %e x: ',n,lambda,err);
    fprintf(repmat('%e ',1,N),x);
    fprintf('\n');
end

while( err > ConvCrit && n < maxIter)

    n = n+1;

    % calculate Jacobian
    for i = 1:N
        xtemp = x;
        h = sqrt(eps)*x(i);
        xtemp(i) = xtemp(i)+h;
        JAC(i,:) = (FUN(xtemp)-F)./h;
    end

    sTemp = (JAC\-F);

    % Find optimal step size

    if(~(max(abs(lambda * sTemp)) < max(abs(s)) ))

        % step size is not gettin' shorter. So we decrease relaxation factor
        while( ( max(abs(s)) < max(abs(lambda * sTemp))))

            % till step size decreases
            lambda = lambda / 2;
            if( lambda < Lmin )
                fprintf('Damping to strong. No convergence.\n');
                return;
            end
        end % step size is ok now
    end

    % next newton step
    s = sTemp;
    x = x + lambda*s;
    F = FUN(x);

    lambda = min( 1., lambda * 2);% increase relaxation factor

    err = sqrt(sum(F.^2)/N);

  if(strcmp(Display,'on'));
    fprintf('n: %i rlx fac: %e error: %e x: ',n,lambda,err);
    fprintf(repmat('%e ',1,N),x);
    fprintf('\n');
  end

end % End of while loop

% -------------------------------------------------------------------------
% Some Postprocessing
% -------------------------------------------------------------------------
if(n<maxIter)
    if(strcmp(Display,'on'));
        disp('Solver converged.');
    end
    EXITFLAG = 1;
else
    if(strcmp(Display,'on'));
        disp('No Solution found... Sorry.');
    end
    EXITFLAG = -1;
end

end


text:
function q = mat2quat(mat,varargin)
% converts direction cosine matrix to quaternion
%
% Syntax
%   q = mat2quat(mat)
%
% Input
%  mat - vector of matrixes
%
% Output
%  q - @quaternion
%
% See also
%
% quaternion_matrix Euler axis2quat hr2quat
%
% Description
% Wertz says to the algo similar to this with largest divisor
%       q4 = 1/2*sqrt((1+mat(1,1)+mat(2,2)+mat(3,3)));
% Eqn 12-14a - c
%  Quat(1) = (mat(2,3)-mat(3,2))/q4/4;
%  Quat(2) = (mat(3,1)-mat(1,3))/q4/4;
%  Quat(3) = (mat(1,2)-mat(2,1))/q4/4;
%  Quat(4) = q4

Quat=zeros(4,1);

%Compute absolute values of the four quaternions from diags of Eqn 12-13
%absQ=0.5*sqrt([1 -1 -1;-1 1 -1;-1 -1 1;1 1 1]*diag(mat)+1);

absQ(1,:) = 0.5 * sqrt(1+mat(1,1,:)+mat(2,2,:)+mat(3,3,:));
absQ(2,:) = 0.5 * sqrt(1-mat(1,1,:)-mat(2,2,:)+mat(3,3,:));
absQ(3,:) = 0.5 * sqrt(1+mat(1,1,:)-mat(2,2,:)-mat(3,3,:));
absQ(4,:) = 0.5 * sqrt(1-mat(1,1,:)+mat(2,2,:)-mat(3,3,:));

[~,ind]=max(absQ); % Select biggest for best accuracy

qind = ind == 1;
if any(qind)
  Quat(1,qind)=absQ(1,qind);
  Quat(2,qind)=squeeze((mat(2,3,qind)-mat(3,2,qind))).'.*0.25./absQ(1,qind);
  Quat(3,qind)=squeeze((mat(3,1,qind)-mat(1,3,qind))).'.*0.25./absQ(1,qind);
  Quat(4,qind)=squeeze((mat(1,2,qind)-mat(2,1,qind))).'.*0.25./absQ(1,qind);
end

qind = ind == 2;
if any(qind)
  Quat(1,qind)=squeeze(mat(1,2,qind)-mat(2,1,qind)).'.*0.25./absQ(2,qind);
  Quat(2,qind)=squeeze(mat(3,1,qind)+mat(1,3,qind)).'.*0.25./absQ(2,qind);
  Quat(3,qind)=squeeze(mat(3,2,qind)+mat(2,3,qind)).'.*0.25./absQ(2,qind);
  Quat(4,qind)=absQ(2,qind);
end

qind = ind == 3;
if any(qind)
  Quat(1,qind)=squeeze(mat(2,3,qind)-mat(3,2,qind)).'.*0.25./absQ(3,qind);
  Quat(2,qind)=absQ(3,qind);
  Quat(3,qind)=squeeze(mat(1,2,qind)+mat(2,1,qind)).'.*0.25./absQ(3,qind);
  Quat(4,qind)=squeeze(mat(3,1,qind)+mat(1,3,qind)).'.*0.25./absQ(3,qind);
end

qind = ind == 4;
if any(qind)
  Quat(1,qind)=squeeze(mat(3,1,qind)-mat(1,3,qind)).'.*0.25./absQ(4,qind);
  Quat(2,qind)=squeeze(mat(1,2,qind)+mat(2,1,qind)).'.*0.25./absQ(4,qind);
  Quat(3,qind)=absQ(4,qind);
  Quat(4,qind)=squeeze(mat(2,3,qind)+mat(3,2,qind)).'.*0.25./absQ(4,qind);
end

q = quaternion(real(Quat));
q = q./norm(q);
q = inv(q).';

author: mtex-toolbox

repo: mtex

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


text:
if isoctave
  status=2;
else
  status=0;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [ux, vx, wx] = spm_fx_cmc_tfm_gen(x,u,P,M,option)
% Generate pre synaptic signals for multimodal DCM for fMRI and M/EEG
% FORMAT [u,v,w]  = spm_fx_cmc_tfm_gen(x,u,P,M)
% FORMAT [u,v]    = spm_fx_cmc_tfm_gen(x,u,P,M)
% FORMAT [u]      = spm_fx_cmc_tfm_gen(x,u,P,M)
%
% Inputs:
% -------------------------------------------------------------------------
% x      - state vector
%   x(:,1) - voltage     (spiny stellate cells)
%   x(:,2) - conductance (spiny stellate cells)
%   x(:,3) - voltage     (superficial pyramidal cells)
%   x(:,4) - conductance (superficial pyramidal cells)
%   x(:,5) - current     (inhibitory interneurons)
%   x(:,6) - conductance (inhibitory interneurons)
%   x(:,7) - voltage     (deep pyramidal cells)
%   x(:,8) - conductance (deep pyramidal cells)
% P        - parameters of canonical micro circuits
% u        - exogenous input
% M        - neural-mass model structure
% option   - options array for calculation pre synaptic signals {1 x 4}:
%   option{1} - 'pre' (pre synaptic) or 'de' (decomposed into intrinsic
%               inhibitory, intrinsic excitatory and extrinsic excitatory)
%               NVC drive.
%   option{2} - 'd' (different) or 's' (same) parameters of neurovascular
%               scaling (this option is not used within this function).
%   option{3} - 'int' (only intrinsic neuronal signals are taken to account
%               for simulating presynaptic signals) or 'ext' (external
%               neuronal signals are additional included).
%   option{4} - EX, a 4x1 matrix  with either 0 or 1 elements (order as
%               follows: [x(:,1) x(:,3) x(:,5) x(:,7)]) to exclude or
%               include populations from calculation of pre synaptic signal
%
%  Examples of options               {'pre',   'd',    'int', EX},
%                                    {'pre',   's',    'int', EX},
%                                    {'pre',   'd',    'ext', EX},
%                                    {'pre',   's',    'ext', EX},
%                                    {'de',    's',    'int', EX},
%                                    {'de',    's',    'ext', EX},
%                                    {'de',    'd',    'int', EX},
%                                    {'de',    'd',    'ext', EX},
%
% Outputs:
% -------------------------------------------------------------------------
%  ux = spm_fx_cmc_tfm(x,u,P,M,option)
%  ux  -  simulated presynaptic signal (including or exclude distal regions)
%
% [ux,vx] = spm_fx_cmc_tfm_gen(x,u,P,M,option)
% ux  - intrinsic presynaptic input, (inhibitory)-without external input
% vx  - intrinsic presynaptic input (excitatory)-without external input
%
% [ux,vx,wx] = spm_fx_cmc_tfm_gen(x,u,P,M,,option)
% ux  - intrinsic presynaptic input, (inhibitory)
% vx  - intrinsic presynaptic input (excitatory)
% wx  - extrinsic presynaptic input
%--------------------------------------------------------------------------
% Prior fixed parameter scaling [Defaults]
%
% E  = (forward and backward) extrinsic rates
% G  = intrinsic rates
% D  = propagation delays (intrinsic, extrinsic)
% T  = synaptic time constants
% R  = slope of sigmoid activation function
%
%__________________________________________________________________________
% David O, Friston KJ (2003) A neural mass model for MEG/EEG: coupling and
% neuronal dynamics. NeuroImage 20: 1743-1755
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Amirhossein Jafarian, and Karl Friston
% $Id: spm_fx_cmc_tfm_gen.m 7735 2019-12-02 09:15:27Z peter $

persistent in1 in2 in3 ;

% get dimensions and configure state variables
%--------------------------------------------------------------------------
x  = spm_unvec(x,M.x);                % neuronal states
n  = size(x,1);                       % number of sources

% [default] fixed parameters
%--------------------------------------------------------------------------
E  = [2 1 1 1]*512;                   % extrinsic (forward and backward)
T  = [16 8 1 2]*16;                   % synaptic rate constants
R  = 1;                               % gain of activation function
B  = 0;                               % baseline firing

% Extrinsic connections
%--------------------------------------------------------------------------
% ss = spiny stellate
% sp = superficial pyramidal
% dp = deep pyramidal
% ii = inhibitory interneurons
%--------------------------------------------------------------------------
A{1} = exp(P.A{1})*E(1);              % forward  connections (sp -> ss)
A{2} = exp(P.A{2})*E(2);              % forward  connections (sp -> dp)
A{3} = exp(P.A{3})*E(3);              % backward connections (dp -> sp)
A{4} = exp(P.A{4})*E(4);              % backward connections (dp -> ii)

% input connections
%--------------------------------------------------------------------------
C    = exp(P.C);

% pre-synaptic inputs: s(V)
%--------------------------------------------------------------------------
ig   = 2*(1:4);                       % indices of conductance
iv   = ig - 1;                        % indices of voltage
V    = x(:,iv);                       % Voltage
F    = 1./(1 + exp(B - R*V));         % firing rate
S    = F - 1/(1 + exp(B));            % deviation

% input
%==========================================================================
if isfield(M,'u')

    % endogenous input
    %----------------------------------------------------------------------
    U   = u(:)*128;
    M.m = size(U,1);

else
    % exogenous input
    %----------------------------------------------------------------------
    U   = C*u(:)*2;
    M.m = size(C,2);

end
clear u

% time constants and intrinsic connections
%==========================================================================
T      = ones(n,1)*T;
i      = 1:size(P.T,2);
T(:,i) = T(:,i).*exp(P.T);

% extrinsic connections
%--------------------------------------------------------------------------
% forward  A{1}  sp -> ss
% forward  A{2}  sp -> dp
% backward A{3}  dp -> sp
% backward A{4}  dp -> ii
%--------------------------------------------------------------------------

% Neuronal states (deviations from baseline)
%--------------------------------------------------------------------------
%   x(:,1) - voltage     (spiny stellate cells)
%   x(:,2) - conductance (spiny stellate cells)but
%   x(:,3) - voltage     (superficial pyramidal cells)
%   x(:,4) - conductance (superficial pyramidal cells)
%   x(:,5) - current     (inhibitory interneurons)
%   x(:,6) - conductance (inhibitory interneurons)
%   x(:,7) - voltage     (deep pyramidal cells)
%   x(:,8) - conductance (deep pyramidal cells)
%--------------------------------------------------------------------------
%     ss sp ii dp   % intrinsic connections
%--------------------------------------------------------------------------
g  = [-8 -4  -4  0;  % ss
       4 -8  -2  0;  % sp
       4  2  -4 -2;  % ii
       0  1  -2 -4]; % dp

g  = g*256*exp(P.S);

% intrinsic connections to be optimised (only the first is modulated)
%--------------------------------------------------------------------------
G       = ones(n,1)*diag(g)';
i       = 1:size(P.G,2);
G(:,i)  = G(:,i).*exp(P.G);

% Modulatory effects of sp depolarisation on recurrent inhibition
%--------------------------------------------------------------------------
if isfield(P,'M')
    G(:,2) = G(:,2).*exp(-P.M*32*S(:,2));
end

% Motion of states: f(x)
%--------------------------------------------------------------------------
% Afferents
%==========================================================================
u      = zeros(n,4);                % intrinsic - inhibitory
v      = zeros(n,4);                % intrinsic - excitatory
w      = zeros(n,4);                % extrinsic - excitatory
pre_ext    = zeros(n,4);            % pre-synaptics with extrinsic input
pre_no_ext = zeros(n,4);            % pre-synaptics without extrinsic input
pincluded  = option{4};             % included populations (binary vector)

% Granular layer (excitatory interneurons): spiny stellate: Hidden causes
%--------------------------------------------------------------------------
u(:,1) = G(:,1).*S(:,1) + g(1,3)*S(:,3) + g(1,2)*S(:,2);
w(:,1) = A{1}*S(:,2) + U;
pre_no_ext(:,1) = u(:,1);
pre_ext(:,1)    = u(:,1) + w(:,1);

% Supra-granular layer (superficial pyramidal cells): Hidden causes - error
%--------------------------------------------------------------------------
u(:,2) = G(:,2).*S(:,2) + g(2,3)*S(:,3);
v(:,2) = g(2,1)*S(:,1);
w(:,2) = A{3}*S(:,4);
pre_no_ext(:,2) = u(:,2) + v(:,2);
pre_ext(:,2)    = u(:,2) + v(:,2) + w(:,2);

% Supra-granular layer (inhibitory interneurons): Hidden states - error
%--------------------------------------------------------------------------
u(:,3) = G(:,3).*S(:,3);
v(:,3) = g(3,1)*S(:,1) + g(3,4)*S(:,4) + g(3,2)*S(:,2);
w(:,3) = A{4}*S(:,4);
pre_no_ext(:,3) = u(:,3) + v(:,3);
pre_ext(:,3)    = u(:,3) + v(:,3) + w(:,3);

% Infra-granular layer (deep pyramidal cells): Hidden states
%--------------------------------------------------------------------------
u(:,4) = G(:,4).*S(:,4) + g(4,3)*S(:,3);
v(:,4) = g(4,2)*S(:,2);
w(:,4) = A{2}*S(:,2);
pre_no_ext(:,4) = u(:,4) + v(:,4);
pre_ext(:,4)    = u(:,4) + v(:,4) + w(:,4);

% Switch off unused connections
%--------------------------------------------------------------------------
u(:,~pincluded)=0;
v(:,~pincluded)=0;
w(:,~pincluded)=0;
pre_ext(:,~pincluded)=0;
pre_no_ext(:,~pincluded)=0;

% output
%--------------------------------------------------------------------------
t = 1 ;
if (strcmp(option{3}, 'ext')&& strcmp(option{1}, 'pre'))
    in1(:,:,t)  = pre_ext;
    ux(:,:,t)    = in1;
end
if(strcmp(option{3}, 'int') && strcmp(option{1}, 'pre'))
    in1(:,:,t)  = pre_no_ext;
    ux(:,:,t)    = in1;
end
if(strcmp(option{1}, 'de'))
    in1(:,:,t)  = u ;
    in2(:,:,t)  = v;
    in3(:,:,t)  = w;
    ux(:,:,t)  = in1;
    vx(:,:,t)  = in2;
    wx(:,:,t)  = in3;
end

author: spm


text:
function F = acsch(F, varargin)
%ACSCH   Inverse hyperbolic cosecant of a CHEBFUN.
%   ACSCH(F) computes the inverse hyperbolic cosecant of the CHEBFUN F.
%
%   ACSCH(F, PREF) does the same but uses the CHEBFUNPREF object PREF when
%   computing the composition.
%
% See also CSCH.

% Copyright 2017 by The University of Oxford and The Chebfun Developers. See
% http://www.chebfun.org/ for Chebfun information.

% Call the compose method:
F = compose(F, @acsch, varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ B ] = normalizerow( A ) %#codegen
  % NORMALIZEROW Normalize each row so that each row's l2 norm as a vector is 1
  %
  % [ B ] = normalizerow( A )
  %
  % Input:
  %  A  #A by D list of row vectors of dimension D
  % Output:
  %  B  #B by D list of normalized row vectors
  %
  % Copyright 2011, Alec Jacobson (jacobson@inf.ethz.ch), Daniele Panozzo
  %

  if issparse(A)
    % speed up (20x) for large sparse matrices
    B = bsxfun(@times,A,1./sqrt(sum(A.^2,2)));
  else
    % normrow will use robust norm
    B = A./normrow(A);
  end
end

author: alecjacobson

repo: gptoolbox


text:
% GETCHANLIST - Obtain indices of specified channel types.
%
% Usage:
%       >> chanlist = getchanlist(chanlocs, type)
%
% Inputs:
%     chanlocs - EEGLAB channel location structure
%     type     - [string] select channel of specified type
%                can enter a cell array to select several channel types
%
% Output:
%     chanlist - list of channel indices for the selected type(s)
%                sorted in increasing order
%
% Note: this function is not case sensitive
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, 2004
%
% See also: POP_CHANEDIT

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

function chanlist = getchanlist(chanlocs, type)

    if nargin < 1
        help getchanlist;
        return;
    end

    if nargin < 2 || ~isfield(chanlocs, 'type')
        chanlist = [1:length(chanlocs)];
        return;
    end

    % search channel types
    % --------------------
    if ischar(type), type = { type }; end
    type = lower(type);
    chanlist = [];
    alltypes = lower({ chanlocs.type });
    for index = 1:length(type)
        tmplist = strmatch(type{index}, alltypes, 'exact');
        if isempty(tmplist)
            fprintf('Warning: no channel of type ''%s'' found\n', type{index});
        end
        chanlist = [ chanlist tmplist' ];
    end
    chanlist = sort(chanlist);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [w, A, C, sbc, fpe, th]=ARFIT_arfit(v, pmin, pmax, selector, no_const)
%ARFIT	Stepwise least squares estimation of multivariate AR model.
%
%  [w,A,C,SBC,FPE,th]=ARFIT(v,pmin,pmax) produces estimates of the
%  parameters of an m-variate AR model of order p,
%
%      v(k,:)' = w' + A1*v(k-1,:)' +...+ Ap*v(k-p,:)' + noise(C),
%
%  where w is the (m x 1) intercept vector, A1, ..., Ap are (m x m)
%  coefficient matrices, and C is a (m x m) noise covariance
%  matrix. The estimated order p lies between pmin and pmax and is
%  chosen as the optimizer of Schwarz's Bayesian Criterion.
%
%  The input matrix v must contain the time series data, with
%  columns v(:,l) representing m variables l=1,...,m and rows
%  v(k,:) representing n observations at different (equally
%  spaced) times k=1,..,n. Optionally, v can have a third
%  dimension, in which case the matrices v(:,:, itr) represent
%  the realizations (e.g., measurement trials) itr=1,...,ntr of the
%  time series. ARFIT returns least squares estimates of the
%  intercept vector w, of the coefficient matrices A1,...,Ap (as
%  A=[A1 ... Ap]), and of the noise covariance matrix C.
%
%  As order selection criteria, ARFIT computes approximations to
%  Schwarz's Bayesian Criterion and to the logarithm of Akaike's Final
%  Prediction Error. The order selection criteria for models of order
%  pmin:pmax are returned as the vectors SBC and FPE.
%
%  The matrix th contains information needed for the computation of
%  confidence intervals. ARMODE and ARCONF require th as input
%  arguments.
%
%  If the optional argument SELECTOR is included in the function call,
%  as in ARFIT(v,pmin,pmax,SELECTOR), SELECTOR is used as the order
%  selection criterion in determining the optimum model order. The
%  three letter string SELECTOR must have one of the two values 'sbc'
%  or 'fpe'. (By default, Schwarz's criterion SBC is used.) If the
%  bounds pmin and pmax coincide, the order of the estimated model
%  is p=pmin=pmax.
%
%  If the function call contains the optional argument 'zero' as the
%  fourth or fifth argument, a model of the form
%
%         v(k,:)' = A1*v(k-1,:)' +...+ Ap*v(k-p,:)' + noise(C)
%
%  is fitted to the time series data. That is, the intercept vector w
%  is taken to be zero, which amounts to assuming that the AR(p)
%  process has zero mean.
%
%  Modified 14-Oct-00
%           24-Oct-10 Tim Mullen (added support for multiple realizatons)
%
%  Authors: Tapio Schneider
%           tapio@gps.caltech.edu
%
%           Arnold Neumaier
%           neum@cma.univie.ac.at

  % n:   number of time steps (per realization)
  % m:   number of variables (dimension of state vectors)
  % ntr: number of realizations (trials)
  [n,m,ntr]   = size(v);

  if (pmin ~= round(pmin) | pmax ~= round(pmax))
    error('Order must be integer.');
  end
  if (pmax < pmin)
    error('PMAX must be greater than or equal to PMIN.')
  end

  % set defaults and check for optional arguments
  if (nargin == 3)              % no optional arguments => set default values
    mcor       = 1;               % fit intercept vector
    selector   = 'sbc';	          % use SBC as order selection criterion
  elseif (nargin == 4)          % one optional argument
    if strcmp(selector, 'zero')
      mcor     = 0;               % no intercept vector to be fitted
      selector = 'sbc';	          % default order selection
    else
      mcor     = 1; 		  % fit intercept vector
    end
  elseif (nargin == 5)          % two optional arguments
    if strcmp(no_const, 'zero')
      mcor     = 0;               % no intercept vector to be fitted
    else
      error(['Bad argument. Usage: ', ...
	     '[w,A,C,SBC,FPE,th]=AR(v,pmin,pmax,SELECTOR,''zero'')'])
    end
  end

  ne  	= ntr*(n-pmax);         % number of block equations of size m
  npmax	= m*pmax+mcor;          % maximum number of parameter vectors of length m

  if (ne <= npmax)
    error('Time series (N = %u) too short.',size(v,1))
  end

  % compute QR factorization for model of order pmax
  [R, scale]   = ARFIT_arqr(v, pmax, mcor);

  % compute approximate order selection criteria for models
  % of order pmin:pmax
  [sbc, fpe]   = ARFIT_arord(R, m, mcor, ne, pmin, pmax);

  % get index iopt of order that minimizes the order selection
  % criterion specified by the variable selector
  [val, iopt]  = min(eval(selector));

  % select order of model
  popt         = pmin + iopt-1; % estimated optimum order
  np           = m*popt + mcor; % number of parameter vectors of length m

  % decompose R for the optimal model order popt according to
  %
  %   | R11  R12 |
  % R=|          |
  %   | 0    R22 |
  %
  R11   = R(1:np, 1:np);
  R12   = R(1:np, npmax+1:npmax+m);
  R22   = R(np+1:npmax+m, npmax+1:npmax+m);

  % get augmented parameter matrix Aaug=[w A] if mcor=1 and Aaug=A if mcor=0
  if (np > 0)
    if (mcor == 1)
      % improve condition of R11 by re-scaling first column
      con 	= max(scale(2:npmax+m)) / scale(1);
      R11(:,1)	= R11(:,1)*con;
    end;
    Aaug = (R11\R12)';

    %  return coefficient matrix A and intercept vector w separately
    if (mcor == 1)
      % intercept vector w is first column of Aaug, rest of Aaug is
      % coefficient matrix A
      w = Aaug(:,1)*con;        % undo condition-improving scaling
      A = Aaug(:,2:np);
    else
      % return an intercept vector of zeros
      w = zeros(m,1);
      A = Aaug;
    end
  else
    % no parameters have been estimated
    % => return only covariance matrix estimate and order selection
    % criteria for ``zeroth order model''
    w   = zeros(m,1);
    A   = [];
  end

  % return covariance matrix
  dof   = ne-np;                % number of block degrees of freedom
  C     = R22'*R22./dof;        % bias-corrected estimate of covariance matrix

  % for later computation of confidence intervals return in th:
  % (i)  the inverse of U=R11'*R11, which appears in the asymptotic
  %      covariance matrix of the least squares estimator
  % (ii) the number of degrees of freedom of the residual covariance matrix
  invR11 = inv(R11);
  if (mcor == 1)
    % undo condition improving scaling
    invR11(1, :) = invR11(1, :) * con;
  end
  Uinv   = invR11*invR11';
  th     = [dof zeros(1,size(Uinv,2)-1); Uinv];

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function RelPath = getRelativeFilePath(FullPath, RootPath, FlagRequireSubdir)
% getRelativeFilePath - Utility to return a relative file path within a root
%
% This function will find the relative file path, given a full absolute
% path and a root path
%
% Syntax:
%       RelPath = wt.utility.getRelativeFilePath(FullPath, RootPath)
%
% Inputs:
%       FullPath - the full absolute path to a file or folder
%       RootPath - the root folder to get a relative path for
%       FlagRequireSubdir - optional flag indicating whether FullPath must
%           be a subdirectory of RootPath [(true)|false]
%
% Outputs:
%       RelPath - the relative path
%
% Examples:
%
%     >> FullPath = 'C:\Program Files\MATLAB\R2016b\toolbox'
%     >> RootPath = 'C:\Program Files\MATLAB'
%     >> RelPath = wt.utility.getRelativeFilePath(FullPath, RootPath)
%
%     RelPath =
%          \R2016b\toolbox
%
% Notes:
%   If FullPath is not a subdirectory of RootPath and FlagRequireSubdir is
%   false, the path will contain parent directory separators "..\"
%

% Copyright 2016-2020 The MathWorks Inc.
% ---------------------------------------------------------------------

% Validate inputs
arguments
    FullPath (1,1) string
    RootPath (1,1) string
    FlagRequireSubdir (1,1) logical = true
end

% Is RootPath empty?
if ~strlength(FullPath)

    RelPath = RootPath;

elseif ~strlength(RootPath)

    RelPath = FullPath;

else

    % Remove trailing filesep
    if endsWith(RootPath,filesep)
        RootPath{1}(end) = '';
    end
    if endsWith(FullPath,filesep)
        FullPath{1}(end) = '';
    end

    % Split the paths apart
    RootParts = strsplit(RootPath,filesep);
    FullParts = strsplit(FullPath,filesep);

    % Find where the paths diverge
    idx = 1;
    SmallestPath = min(numel(RootParts), numel(FullParts));
    while idx<=SmallestPath && strcmpi(RootParts(idx),FullParts(idx))
        idx = idx+1;
    end

    % Is the specified path outside of the root directory?
    NumAbove = max(numel(RootParts) - idx + 1, 0);
    if FlagRequireSubdir && NumAbove>0
        error('The specified path:\n\t"%s"\nis not a subdirectory of the root path:\n\t"%s"',...
            char(FullPath), char(RootPath) );
    else
        % In case full path is above the RootPath, add ".." paths
        ParentPaths = string(repmat(['..' filesep],1,NumAbove));

        % Form the relative path
        RelPath = filesep + fullfile(ParentPaths, FullParts{idx:end});
    end

    % What if paths are still the same?
    if isempty(RelPath)
        RelPath = "." + filesep;
    end

end %if isempty(RootPath)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function pass = test_optimization(pref)
% Can we do global optimization over [0 1]^3?

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1000*pref.cheb3Prefs.chebfun3eps;

dom = [0 1 0 1 0 1];

Battery = {@(x,y,z) cos(pi*x.*y.*z),...
    @(x,y,z) cos(2*pi*x.*y.*z), ...
    @(x,y,z) cos(3*pi*x.*y.*z),...
    @(x,y,z) cos(4*pi*x.*y.*z),...
    @(x,y,z) cos(5*pi*x.*y.*z),...
    @(x,y,z) cos(6*pi*x.*y.*z),...
    @(x,y,z) cos(7*pi*x.*y.*z),...
    @(x,y,z) sin(pi*x.*y.*z),...
    @(x,y,z) cos(0*pi*(x-y-z).^2),...
    @(x,y,z) sin(x+y+z)
    };

maxima = [1
    1
    1
    1
    1
    1
    1
    1
    1
    1];

minima = [-1
    -1
    -1
    -1
    -1
    -1
    -1
     0
    +1
     0];

for j = 1:length(Battery)
    f = Battery{j};
    g = chebfun3(f, dom);
     [Y, ignored] = minandmax3(g);
    err(j) = norm(Y(1) - minima(j)) + norm(Y(2) - maxima(j));
end

%%
pass = err < tol;

end


text:
function [x0,LB,UB,PLB,PUB] = boundscheck_vbmc(x0,LB,UB,PLB,PUB,prnt)
%BOUNDSCHECK Initial check of bounds.

[N0,nvars] = size(x0);

% Expand scalar inputs to full vectors
if isscalar(LB); LB = LB*ones(1,nvars); end
if isscalar(UB); UB = UB*ones(1,nvars); end
if isscalar(PLB); PLB = PLB*ones(1,nvars); end
if isscalar(PUB); PUB = PUB*ones(1,nvars); end

if isempty(PLB) || isempty(PUB)
    if N0 > 1
        if prnt > 0
            fprintf('PLB and/or PUB not specified. Estimating plausible bounds from starting set X0...\n');
        end
        width = max(x0) - min(x0);
        if isempty(PLB)
            PLB = min(x0) - width/N0;
            PLB = max(PLB,LB);
        end
        if isempty(PUB)
            PUB = max(x0) + width/N0;
            PUB = min(PUB,UB);
        end

        idx = any(PLB == PUB);
        if any(idx)
            PLB(idx) = LB(idx);
            PUB(idx) = UB(idx);
            warning('vbmc:pbInitFailed', ...
                'Some plausible bounds could not be determined from starting set. Using hard upper/lower bounds for those instead.')
        end
    else
        warning('vbmc:pbUnspecified', ...
            'Plausible lower/upper bounds PLB and/or PUB not specified and X0 is not a valid starting set. Using hard upper/lower bounds instead.');
        if isempty(PLB); PLB = LB; end
        if isempty(PUB); PUB = UB; end
    end
end

% Test that all vectors have the same length
if any([numel(LB),numel(UB),numel(PLB),numel(PUB)] ~= nvars)
    error('All input vectors (X0, LB, UB, PLB, PUB), if specified, need to have the same size.');
end

% Test that all vectors are row vectors
if ~isvector(LB) || ~isvector(UB) || ~isvector(PLB) || ~isvector(PUB) ...
        || size(LB,1) ~= 1 || size(UB,1) ~= 1 || size(PLB,1) ~= 1 || size(PUB,1) ~= 1
    error('All input vectors LB, UB, PLB, PUB, if specified, should be row vectors.');
end

% Test that plausible bounds are finite
if ~all(isfinite(([PLB, PUB])))
    error('Plausible interval bounds PLB and PUB need to be finite.');
end

% Test that all vectors are real-valued
if ~isreal([x0(:)', LB, UB, PLB, PUB])
    error('All input vectors should be real-valued.');
end

% Fixed variables (all bounds equal) are not supported
fixidx = (LB == UB) & (UB == PLB) & (PLB == PUB);
if any(fixidx)
    error('vbmc:FixedVariables', ...
        'VBMC does not support fixed variables. Lower and upper bounds should be different.');
end

% Test that plausible bounds are different
if any(PLB == PUB)
    error('vbmc:MatchingPB', ...
        'For all variables, plausible lower and upper bounds need to be distinct.')
end

% Check that all X0 are inside the bounds
if any(any(bsxfun(@lt,x0,LB))) || any(any(bsxfun(@gt,x0,UB)))
    error('vbmc:InitialPointsNotInsideBounds', ...
        'The starting points X0 are not inside the provided hard bounds LB and UB.');
end

% Compute "effective" bounds (slightly inside provided hard bounds)
bounds_range = UB - LB;
bounds_range(isinf(bounds_range)) = 1e3;
scale_factor = 1e-3;
LB_eff = LB + scale_factor*bounds_range;
LB_eff(abs(LB) <= realmin) = scale_factor*bounds_range(abs(LB) <= realmin);
UB_eff = UB - scale_factor*bounds_range;
UB_eff(abs(UB) <= realmin) = -scale_factor*bounds_range(abs(UB) <= realmin);
LB_eff(isinf(LB)) = LB(isinf(LB));  % Infinities stay the same
UB_eff(isinf(UB)) = UB(isinf(UB));

if any(LB_eff >= UB_eff)
    error('vbmc:StrictBoundsTooClose', ...
        'Hard bounds LB and UB are numerically too close. Make them more separate.');
end

% Fix when provided X0 are almost on the bounds -- move them inside
if any(any(bsxfun(@le,x0,LB_eff))) || any(any(bsxfun(@ge,x0,UB_eff)))
    warning('vbmc:InitialPointsTooClosePB', ...
        'The starting points X0 are on or numerically too close to the hard bounds LB and UB. Moving the initial points more inside...');
    x0 = bsxfun(@max, bsxfun(@min,x0,UB_eff), LB_eff);
end

% Test order of bounds (permissive)
ordidx = LB <= PLB & PLB < PUB & PUB <= UB;
if any(~ordidx)
    error('vbmc:StrictBounds', ...
        'For each variable, hard and plausible bounds should respect the ordering LB < PLB < PUB < UB.');
end

% Test that plausible bounds are reasonably separated from hard bounds
ordidx = LB_eff < PLB & PUB < UB_eff;
if any(~ordidx)
    warning('vbmc:TooCloseBounds', ...
        'For each variable, hard and plausible bounds should not be too close. Moving plausible bounds.');
    PLB = max(PLB,LB_eff);
    PUB = min(PUB,UB_eff);
end

% Check that all X0 are inside the plausible bounds, move bounds otherwise
if any(any(bsxfun(@le,x0,PLB))) || any(any(bsxfun(@ge,x0,PUB)))
    warning('vbmc:InitialPointsOutsidePB', ...
        'The starting points X0 are not inside the provided plausible bounds PLB and PUB. Expanding the plausible bounds...');
    PLB = min(PLB,min(x0,[],1));
    PUB = max(PUB,max(x0,[],1));
end

% Test order of bounds
ordidx = LB < PLB & PLB < PUB & PUB < UB;
if any(~ordidx)
    error('vbmc:StrictBounds', ...
        'For each variable, hard and plausible bounds should respect the ordering LB < PLB < PUB < UB.');
end

% Test that variables are either bounded or unbounded (not half-bounded)
halfbnd = (isinf(LB) & isfinite(UB)) | (isfinite(LB) & isinf(UB));
if any(halfbnd)
    error('vbmc:HalfBounds', ...
        'Each variable needs to be unbounded or bounded. Variables bounded only below/above are not supported.');
end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function [amp,phi,yr]=HANTS(ni,nb,nf,y,ts,HiLo,low,high,fet,dod,delta)
% HANTS processing
%
% Wout Verhoef
% NLR, Remote Sensing Dept.
% June 1998
%
% Converted to MATLAB:
% Mohammad Abouali (2011)
%
% NOTE: This version is tested in MATLAB V2010b. In some older version you
% might get an error on line 117. Refer to the solution provided on that
% line.
%
% Modified:
%   Apply suppression of high amplitudes for near-singular case by
%	adding a number delta to the diagonal elements of matrix A,
%	except element (1,1), because the average should not be affected
%
%	Output of reconstructed time series in array yr June 2005
%
%   Change call and input arguments to accommodate a base period length (nb)
%   All frequencies from 1 (base period) until nf are included
%
% Parameters
%
% Inputs:
%   ni    = nr. of images (total number of actual samples of the time
%           series)
%   nb    = length of the base period, measured in virtual samples
%           (days, dekads, months, etc.)
%   nf    = number of frequencies to be considered above the zero frequency
%   y     = array of input sample values (e.g. NDVI values)
%   ts    = array of size ni of time sample indicators
%           (indicates virtual sample number relative to the base period);
%           numbers in array ts maybe greater than nb
%           If no aux file is used (no time samples), we assume ts(i)= i,
%           where i=1, ..., ni
%   HiLo  = 2-character string indicating rejection of high or low outliers
%   low   = valid range minimum
%   high  = valid range maximum (values outside the valid range are rejeced
%           right away)
%   fet   = fit error tolerance (points deviating more than fet from curve
%           fit are rejected)
%   dod   = degree of overdeterminedness (iteration stops if number of
%           points reaches the minimum required for curve fitting, plus
%           dod). This is a safety measure
%   delta = small positive number (e.g. 0.1) to suppress high amplitudes
%
% Outputs:
%
% amp   = returned array of amplitudes, first element is the average of
%         the curve
% phi   = returned array of phases, first element is zero
% yr	= array holding reconstructed time series

mat=zeros(min(2*nf+1,ni),ni,'single');
amp=zeros(nf+1,1,'single');
phi=zeros(nf+1,1,'single');
yr=zeros(ni,1,'double');

% if (Opt.FirstRun==true)
    sHiLo = 0;
    if (strcmp(HiLo,'Hi'))
        sHiLo =-1;
    end
    if (strcmp(HiLo,'Lo'))
        sHiLo = 1;
    end
    nr=min(2*nf+1,ni);
    noutmax=ni-nr-dod;
    dg=180.0/pi;
    mat(1,:)=1.0;

    ang=2.*pi*(0:nb-1)/nb;
    cs=cos(ang);
    sn=sin(ang);
%     Opt.FirstRun=false;
% end
i=1:nf;
for j=1:ni
    index=1+mod(i*(ts(j)-1),nb);
    mat(2*i  ,j)=cs(index);
    mat(2*i+1,j)=sn(index);
end

p=ones(ni,1);
p(or(y<low,y>high))=0;
nout=sum(p==0);

if (nout>noutmax)
%     disp('Not enough data points')
%     disp(['nout:' num2str(nout) ' ,noutmax:' num2str(noutmax)])
%     error('nout > noutmax')
    return
end

ready=false;
nloop=0;
nloopmax=ni;

while ((~ready)&&(nloop<nloopmax))
    nloop=nloop+1;
    za=mat*(p.*y);

    A=mat*diag(p)*mat';
    A=A+diag(ones(nr,1))*delta;
    A(1,1)=A(1,1)-delta;
    zr=A\za;

    yr=mat'*zr;
    diffVec=sHiLo*(yr-y);
    err=p.*diffVec;

	[~, rankVec]=sort(err,'ascend');
% The above line may not be recognized on some older MATLAB versions.
% Simply comment the above line and uncomment the line below.
%    [tmp, rankVec]=sort(err,'ascend');

    maxerr=diffVec(rankVec(ni));
    ready=(maxerr<=fet)||(nout==noutmax);

    if (~ready)
        i=ni;
        j=rankVec(i);
        while ( (p(j)*diffVec(j)>maxerr*0.5)&&(nout<noutmax) )
				p(j)=0;
				nout=nout+1;
				i=i-1;
				j=rank(i);
        end
    end
end

amp(1)=zr(1);
phi(1)=0.0;

zr(ni+1)=0.0;

i=2:2:nr;
ifr=(i+2)/2;
ra=zr(i);
rb=zr(i+1);
amp(ifr)=sqrt(ra.*ra+rb.*rb);
phase=atan2(rb,ra)*dg;
phase(phase<0)=phase(phase<0)+360;
phi(ifr)=phase;

end


text:
function value = r8_huge ( )

%*****************************************************************************80
%
%% R8_HUGE returns a "huge" real number.
%
%  Discussion:
%
%    MATLAB provides a built-in symbolic function realmax ( ) that returns
%    the largest real number.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VALUE, a huge number.
%
  value = 1.79769313486231571E+308;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [x,y,xy,mv,mp,bound] = mg_cavity_domain(nc,grid_type)
%mg_cavity_domain   square cavity Q2 grid generator for GMG
%   [x,y,xy,mv,mp,bound] = mg_cavity_domain(nc,grid_type)
%   input
%          nc           grid level indicator
%          grid_type    1 for uniform, 2 for stretched
%   output
%          grid-defining data
%
%   IFISS function: HCE; 24 January 2004.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage

%   modelled after IFISS function cavity_domain, adapted for GMG
n=2^nc; np=n/2; nq=n/4;
% y-direction
if grid_type==2
   hmax=nc/(2^(nc+1));
   x1=-1;x2=-2*hmax;x3=2*hmax;x4=1;nx1=2^(nc-1)-1;nx2=2;nx3=2^(nc-1)-1;
   y1=-1;y2=-2*hmax;y3=2*hmax;y4=1;ny1=2^(nc-1)-1;ny2=2;ny3=2^(nc-1)-1;
   y=subint(y1,y2,y3,y4,ny1,ny2,ny3);
   stretch=(y(3)-y(2))/(y(2)-y(1));
   x=y;
else
   yy=[1/np:1/np:1];
   ypos=[0,yy];
   yneg=-yy(length(yy):-1:1);
   y=[yneg,ypos]';
   x=y;
end
%
%% compute biquadratic element coordinates
nvtx=(n+1)*(n+1);
[X,Y]=meshgrid(x,y);
xx=reshape(X',nvtx,1);
yy=reshape(Y',nvtx,1);
xy=[xx(:),yy(:)];
%
kx = 1;
ky = 1;
mel=0;
for j=1:np
   for i=1:np
      mref=(n+1)*(ky-1)+kx;
      pref=(np+1)*(j-1)+i;
      mel=mel+1;
      nvv(1) = mref;
      nvv(2) = mref+2;
      nvv(3) = mref+2*n+4;
      nvv(4) = mref+2*n+2;
      nvv(5) = mref+1;
      nvv(6) = mref+n+3;
      nvv(7) = mref+2*n+3;
      nvv(8)=  mref+n+1;
      nvv(9)=  mref+n+2;
      npp(1) = pref;
      npp(2) = pref+1;
      npp(3) = pref+np+2;
      npp(4) = pref+np+1;
      mv(mel,1:9)=nvv(1:9);
      mp(mel,1:4)=npp(1:4);
      kx = kx + 2;
   end
   ky = ky + 2;
   kx = 1;
end

%% compute boundary vertices and edges
% four boundary edges
k1=find( xy(:,2)==-1  );
e1=[]; for k=1:mel, if any(mv(k,5)==k1), e1=[e1,k]; end, end
ef1=ones(size(e1));
%
k2=find( xy(:,1)==1 & xy(:,2)<1   & xy(:,2) >-1);
e2=[]; for k=1:mel, if any(mv(k,6)==k2), e2=[e2,k]; end, end
ef2=2*ones(size(e2));
%
k3=find( xy(:,2)==1  );
e3=[]; for k=1:mel, if any(mv(k,7)==k3), e3=[e3,k]; end, end
ef3=3*ones(size(e3));
%
k4=find( xy(:,1)==-1 & xy(:,2)<1   & xy(:,2) >-1);
e4=[]; for k=1:mel, if any(mv(k,8)==k4), e4=[e4,k]; end, end
ef4=4*ones(size(e4));
%
bound=sort([k1;k2;k3;k4]);
mbound=[e1',ef1';e2',ef2';e3',ef3';e4',ef4'];


text:
function [state, y] = ft_preproc_online_downsample_apply(state, x)

% FT_PREPROC_ONLINE_DOWNSAMPLE_APPLY passes a signal through the online downsampler
% and returns the downsampler state and the downsampled signal. The state keeps track
% of the number of samples to be skipped in the next call.
%
% Use as
%    [state, dat] = ft_preproc_online_downsample_apply(state, x)
% where
%   dat   = Nchan x Ntime
%   state = downsampler state, see FT_PREPROC_ONLINE_DOWNSAMPLE_INIT
%
% See also PREPROC

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

N = size(x,2);

% to get number K of sample we can write out, subtract the skipped samples,
% and then add maximum possible number of skip samples for next time (=state.factor-1)
K  = floor((N - state.numSkip + state.factor-1)/state.factor);

startIdx = 1+state.numSkip;
endIdx   = 1+state.numSkip + (K-1)*state.factor;

y = x(:,startIdx:state.factor:endIdx);
state.numSkip = state.factor-1-(N-endIdx); % for next time

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ x, status ] = p00_newton ( problem, option, nvar, x, par_index )

%*****************************************************************************80
%
%% P00_NEWTON applies Newton's method to an approximate root.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the problem index.
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Input, real X(NVAR), the starting point of Newton's method.
%
%    Input, integer PAR_INDEX, the index of the parameter to be held fixed.
%    This variable should be between 1 and NVAR.  However, the user can
%    set it to 0, indicating that the program should make an intelligent
%    choice for the index.
%
%    Output, real X(NVAR), an improved estimate of the root of F(X)=0.
%
%    Output, integer STATUS, the status of the iteration.
%    -3, the full number of steps was taken without convergence.
%        (however, the output X might be CLOSE to a good solution).
%    -2, the iteration seemed to be diverging, and was halted.
%    -1, the jacobian was singular, and the iteration was halted.
%     nonnegative, the convergence test was satisfied, and this is the
%        number of steps taken (possibly 0).
%
  COND_MAX = 1.0E+10;
  FX_ABS_TOL = 0.000001;
  IT_MAX = 20;
  VERBOSE = 0;

  x = x(:);

  if ( par_index < 1 | nvar < par_index )
    par_index = p00_par_index ( problem, option, nvar, x );
    if ( VERBOSE )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Iteration will hold index %d fixed.\n', par_index );
    end
  end
  par_value = x(par_index);

  if ( VERBOSE )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_NEWTON\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '      Step    ||F(X)||\n' );
    fprintf ( 1, '\n' );
  end

  for it = 0 : IT_MAX
%
%  Compute the function value.
%
    fx = p00_fun ( problem, option, nvar, x );
    fx(nvar) = x(par_index) - par_value;
%
%  Compute the norm of the function value.
%
    fx_max = max ( abs ( fx(1:nvar) ) );

    if ( VERBOSE )
      fprintf ( 1, '  %8d  %14e\n', it, fx_max );
    end

    if ( it == 0 )
      fx_max_init = fx_max;
    end
%
%  If the function norm is small enough, return.
%
    if ( abs ( fx_max ) < FX_ABS_TOL )
      status = it;
      break
    end
%
%  If the function norm seems to be exploding, halt.
%
    if ( 1000.0 * fx_max_init < abs ( fx_max ) )
      status = -2;
      break
    end

    if ( it == IT_MAX )
      status = -3;
      break
    end
%
%  Compute the jacobian.
%
    jac = p00_jac ( problem, option, nvar, x );

    jac(nvar,1:nvar) = 0.0;
    jac(nvar,par_index) = 1.0;
%
%  Solve the system JAC * DX = - FX
%
    if ( cond ( jac ) < COND_MAX )
      dx = - jac \ fx;
    else
      if ( VERBOSE )
        fprintf ( 1, '  (Using pseudoinverse.)\n' );
      end
      dx = - pinv ( jac ) * fx;
    end
%
%  Update X = X - DX.
%
    x(1:nvar) = x(1:nvar) + dx(1:nvar);

  end

  return
end


text:
function writePOLY_pyramid(filename,V,E,F,H)
  % WRITEPOLY_PYRAMID prints vertices, segments and facets to a .poly suitable
  % for use with PYRAMID
  %
  % writePOLY_pyramid(filename,V,E,F,H)
  %
  % Inputs:
  %   V  #V by dim=3 list of vertex positions
  %   E  #E by 2+boundary_markers list of segment indices, indexing V, and
  %     optional boundary markers
  %   F  #F struct containing polygon information arrays
  %     .facets  a #facets list of facets,  each facet is a polygon
  %       **NOTE: facets index E *not* V, contrary to typical (V,F) meshes and
  %       contrary to writePOLY_tetgen prototype
  %     .boundary_markers a #facets list of boundary_markers
  %     OR
  %   F  #F by constant-degree+boundary_markers  list of facets
  %
  % Example:
  %   % Mesh in (V,F)
  %   % Gather all edges
  %   E = [F(:,[2 3]);F(:,[3 1]);F(:,[1 2])];
  %   [~,IE,IuE] = unique(sort(E,2),'rows');
  %   % unique edges
  %   uE = E(IE,:);
  %   % reindex F into E
  %   FinE = reshape(IuE,size(F));
  %   Facets = [];
  %   Facets.facets = mat2cell(FinE,ones(size(FinE,1),1),[3]);
  %   Facets.boundary_marker = -ones(size(FinE,1),1);
  %   writePOLY_pyramid(path,V,uE,Facets,[]);
  %
  % See also: cdt, tetgen, writePOLY_triangle, writePOLY_tetgen
  %
  %

  % open file for writing
  poly_file_handle = fopen(filename,'w');

  dim = size(V,2);

  if dim ~= 3
    error('writePOLY_pyramid is for 3d meshes. Try writePOLY_triangle etc.');
  end

  % vertices section
  fprintf(poly_file_handle,'# Part 1 - node list\n');
  format = '%d %.17g %.17g %.17g\n';
  fprintf(poly_file_handle,'%d %d 0 0\n', size(V,1),size(V,2));
  if ~isempty(V)
    fprintf(poly_file_handle,format,[1:size(V,1);V']);
  end

  fprintf(poly_file_handle,'# Part 2 - segment list\n');
  fprintf(poly_file_handle,'%d %d\n',size(E,1),size(E,2)-2);
  format = ['%d %d %d' repmat(' %d',1,size(E,2)-2) '\n'];
  fprintf(poly_file_handle,format, [1:size(E,1);E']);

  fprintf(poly_file_handle,'# Part 2 - facet list\n');
  % for now, always include boundary markers
  % [num facets] [boundary markers]
  assert(isempty(F) || isempty(F.facets) || iscell(F.facets));
  if isempty(F)
    fprintf(poly_file_handle,'0\n');
  else
      fprintf(poly_file_handle,'%d %d\n',numel(F.facets),~isempty(F.boundary_marker));
      fs = cell2mat(cellfun(@size,F.facets,'UniformOutput',false));
      % Try to print all at once if facets are all the same size
      if ~isempty(fs) && all(fs(:,1) == 1) && all(fs(:,2) == fs(1,2))
          % build format
          fformat = [ ...
              ... % index and size
              '%d ' num2str(fs(1,2)) ...
              ... % facet indices into segments
              repmat(' %d',1,fs(1,2)) ...
              ... % boundary markers
              repmat(' %d',1,size(F.boundary_marker,2)) '\n'];
          % print all at once
          fprintf(poly_file_handle,fformat, ...
              [1:size(F.facets,1);[cell2mat(F.facets) F.boundary_marker]']);
      else
          % irregular face valences
          for f=1:numel(F.facets)
              % 1d list
              assert(size(F.facets{f},1)==1 || size(F.facets{f},2) == 1);
              fprintf('%d %d',f, numel(F.facets{f}));
              % print indices
              for p=1:numel(F.facets{f})
                  fprintf(poly_file_handle,' %d',F.facets{f}(p));
              end
              fprintf(poly_file_handle,'\n');
          end
      end
  end

  % [num holes]
  fprintf(poly_file_handle,'# Part 3 - hole list\n');
  fprintf(poly_file_handle,'%d\n',size(H,1));
  if ~isempty(H)
    assert(isempty(V) || size(H,2) == size(V,2));
    fprintf(poly_file_handle,'%0.17g %0.17g %0.17g\n',[1:size(H,1);H']);
  end

end


text:
classdef HaloCoordAdjustEnum < matlab.mixin.SetGet
    %HaloCoordAdjustEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        X('X_0', 1)
        Z('Z_0', 3);
    end

    properties
        name char = '';
        adjustInd(1,1) double = 0;
    end

    methods
        function obj = HaloCoordAdjustEnum(name, adjustInd)
            obj.name = name;
            obj.adjustInd = adjustInd;
        end
    end

    methods(Static)
        function [listBoxStr, enums] = getListBoxStr()
            m = enumeration('HaloCoordAdjustEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
            enums = m;
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('HaloCoordAdjustEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('HaloCoordAdjustEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end


text:
function [rVect, vVect] = getStateAtTime(bodyInfo, time, gmu)
%getStateAtTime Summary of this function goes here
%   Detailed explanation goes here
    if(isstruct(bodyInfo) || (isprop(bodyInfo,'propTypeEnum') && bodyInfo.propTypeEnum == BodyPropagationTypeEnum.TwoBody) || (numel(time) == 1 && time == bodyInfo.epoch))
        numTimes = length(time);

        oneArray = (zeros(1, numTimes)+1);

        sma = bodyInfo.sma * oneArray;
        ecc = bodyInfo.ecc * oneArray;
        inc = AngleZero2Pi(deg2rad(bodyInfo.inc)) * oneArray;
        raan = AngleZero2Pi(deg2rad(bodyInfo.raan)) * oneArray;
        argp = AngleZero2Pi(deg2rad(bodyInfo.arg)) * oneArray;
        M0 = deg2rad(bodyInfo.mean) * oneArray;

        n = computeMeanMotion(sma, gmu);
        deltaT = time - bodyInfo.epoch;
        M = (M0(:) + n(:).*deltaT(:))';
        tru = computeTrueAnomFromMean(M, ecc);

        if(length(tru) > 1)
            [rVect, vVect] = vect_getStatefromKepler(sma, ecc, inc, raan, argp, tru, gmu, true);
        else
            [rVect, vVect] = getStatefromKepler_Alg(sma, ecc, inc, raan, argp, tru, gmu);
%             [rVect, vVect] = vect_getStatefromKepler(sma, ecc, inc, raan, argp, tru, gmu, true);
        end

    elseif(bodyInfo.propTypeEnum == BodyPropagationTypeEnum.Numerical)
        [rVect, vVect] = bodyInfo.numIntStateCache.getCachedBodyStateAtTime(time);

    else
        error('Unknown celestial body prop sim type.');
    end
end


text:
function varargout = process_aec1( varargin )
% PROCESS_PLV1: Compute amplitude envelope correlation between one signal and all the others, in one file.

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Amplitude Envelope Correlation 1xN';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Connectivity';
    sProcess.Index       = 680;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Connectivity';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data',     'results',  'matrix'};
    sProcess.OutputTypes = {'timefreq', 'timefreq', 'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;

    % === CONNECT INPUT
    sProcess = process_corr1n('DefineConnectOptions', sProcess, 0);
    % === FREQ BANDS
    sProcess.options.freqbands.Comment = 'Frequency bands for the Hilbert transform:';
    sProcess.options.freqbands.Type    = 'groupbands';
    sProcess.options.freqbands.Value   = bst_get('DefaultFreqBands');
    % === Orthogonalize pairs of signals
    sProcess.options.isorth.Comment = 'Orthogonalize signal pairs before envelope computation';
    sProcess.options.isorth.Type    = 'checkbox';
    sProcess.options.isorth.Value   = 0;
    % === OUTPUT MODE
    sProcess.options.outputmode.Comment = {'Save individual results (one file per input file)', 'Concatenate input files before processing (one file)', 'Save average connectivity matrix (one file)'};
    sProcess.options.outputmode.Type    = 'radio';
    sProcess.options.outputmode.Value   = 1;
    sProcess.options.outputmode.Group   = 'output';
end

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

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputA) %#ok<DEFNU>
    % Input options
    OPTIONS = process_corr1n('GetConnectOptions', sProcess, sInputA);
    if isempty(OPTIONS)
        OutputFiles = {};
        return
    end

    OPTIONS.Method = 'aec';
    % Hilbert and frequency bands options
    OPTIONS.Freqs = sProcess.options.freqbands.Value;
    OPTIONS.isOrth = sProcess.options.isorth.Value;

    % Compute metric
    OutputFiles = bst_connectivity({sInputA.FileName}, {sInputA.FileName}, OPTIONS);
end

author: brainstorm-tools


function J=jacob_bvar(param)

% Filippo Ferroni, 6/1/2015
% Revised, 2/15/2017
% Revised, 3/21/2018

J=zeros(length(param),1);
for jj = 1 : length(param)
    J(jj)= bound0prime(param(jj));
end
J=diag(J);

%
% function y = bound01prime(x);
% y = exp(x)/(1+exp(x))^2;
%
function y = bound0prime(x);
y = exp(x);


text:
% pop_dipfit_nonlinear() - interactively do dipole fit of selected ICA components
%
% Usage:
%  >> EEGOUT = pop_dipfit_nonlinear( EEGIN )
%
% Inputs:
%   EEGIN       input dataset
%
% Outputs:
%   EEGOUT      output dataset
%
% Author: Robert Oostenveld, SMI/FCDC, Nijmegen 2003
%         Arnaud Delorme, SCCN, La Jolla 2003
%         Thanks to Nicolas Robitaille for his help on the CTF MEG
%         implementation

% SMI, University Aalborg, Denmark http://www.smi.auc.dk/
% FC Donders Centre, University Nijmegen, the Netherlands http://www.fcdonders.kun.nl/

% Copyright (C) 2003 Robert Oostenveld, SMI/FCDC roberto@smi.auc.dk
%
% 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [EEGOUT, com] = pop_dipfit_nonlinear( EEG, subfunction, parent, dipnum )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% the code for this interactive dialog has 4 major parts
% - draw the graphical user interface
% - synchronize the gui with the data
% - synchronize the data with the gui
% - execute the actual dipole analysis

% the subfunctions that perform handling of the gui are
% - dialog_selectcomponent
% - dialog_checkinput
% - dialog_setvalue
% - dialog_getvalue
% - dialog_plotmap
% - dialog_plotcomponent
% - dialog_flip
% the subfunctions that perform the fitting are
% - dipfit_position
% - dipfit_moment

if nargin<1
  help pop_dipfit_nonlinear;
  return
elseif nargin==1

  EEGOUT = EEG;
  com = '';

  if ~isfield(EEG, 'chanlocs')
    error('No electrodes present');
  end

  if ~isfield(EEG, 'icawinv')
    error('No ICA components to fit');
  end

  if ~isfield(EEG, 'dipfit')
    error('General dipolefit settings not specified');
  end

  if ~isfield(EEG.dipfit, 'vol') & ~isfield(EEG.dipfit, 'hdmfile')
    error('Dipolefit volume conductor model not specified');
  end

  % select all ICA components as 'fitable'
  select = 1:size(EEG.icawinv,2);
  if ~isfield(EEG.dipfit, 'current')
    % select the first component as the current component
    EEG.dipfit.current = 1;
  end

  % verify the presence of a dipole model
  if ~isfield(EEG.dipfit, 'model')
    % create empty dipole model for each component
    for i=select
      EEG.dipfit.model(i).posxyz = zeros(2,3);
      EEG.dipfit.model(i).momxyz = zeros(2,3);
      EEG.dipfit.model(i).rv     = 1;
      EEG.dipfit.model(i).select = [1];
    end
  end

  % verify the size of each dipole model
  for i=select
    if ~isfield(EEG.dipfit.model, 'posxyz') | length(EEG.dipfit.model) < i | isempty(EEG.dipfit.model(i).posxyz)
      % replace all empty dipole models with a two dipole model, of which one is active
      EEG.dipfit.model(i).select = [1];
      EEG.dipfit.model(i).rv = 1;
      EEG.dipfit.model(i).posxyz = zeros(2,3);
      EEG.dipfit.model(i).momxyz = zeros(2,3);
    elseif size(EEG.dipfit.model(i).posxyz,1)==1
      % replace all one dipole models with a two dipole model
      EEG.dipfit.model(i).select = [1];
      EEG.dipfit.model(i).posxyz = [EEG.dipfit.model(i).posxyz; [0 0 0]];
      EEG.dipfit.model(i).momxyz = [EEG.dipfit.model(i).momxyz; [0 0 0]];
    elseif size(EEG.dipfit.model(i).posxyz,1)>2
      % replace all more-than-two dipole models with a two dipole model
      warning('pruning dipole model to two dipoles');
      EEG.dipfit.model(i).select = [1];
      EEG.dipfit.model(i).posxyz = EEG.dipfit.model(i).posxyz(1:2,:);
      EEG.dipfit.model(i).momxyz = EEG.dipfit.model(i).momxyz(1:2,:);
    end
  end

  % default is not to use symmetry constraint
  constr = [];

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % construct the graphical user interface
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % define the callback functions for the interface elements
  cb_plotmap         = 'pop_dipfit_nonlinear(EEG, ''dialog_plotmap'', gcbf);';
  cb_selectcomponent = 'pop_dipfit_nonlinear(EEG, ''dialog_selectcomponent'', gcbf);';
  cb_checkinput      = 'pop_dipfit_nonlinear(EEG, ''dialog_checkinput'', gcbf);';
  cb_fitposition     = 'pop_dipfit_nonlinear(EEG, ''dialog_getvalue'', gcbf); pop_dipfit_nonlinear(EEG, ''dipfit_position'', gcbf); pop_dipfit_nonlinear(EEG, ''dialog_setvalue'', gcbf);';
  cb_fitmoment       = 'pop_dipfit_nonlinear(EEG, ''dialog_getvalue'', gcbf); pop_dipfit_nonlinear(EEG, ''dipfit_moment''  , gcbf); pop_dipfit_nonlinear(EEG, ''dialog_setvalue'', gcbf);';
  cb_close           = 'close(gcbf)';
  cb_help            = 'pophelp(''pop_dipfit_nonlinear'');';
  cb_ok              = 'uiresume(gcbf);';
  cb_plotdip         = 'pop_dipfit_nonlinear(EEG, ''dialog_plotcomponent'', gcbf);';
  cb_flip1           = 'pop_dipfit_nonlinear(EEG, ''dialog_flip'', gcbf, 1);';
  cb_flip2           = 'pop_dipfit_nonlinear(EEG, ''dialog_flip'', gcbf, 2);';
  cb_sym             = [ 'set(findobj(gcbf, ''tag'', ''dip2sel''), ''value'', 1);' cb_checkinput ];

  % vertical layout for each line
  geomvert  =   [1 1 1 1 1 1 1 1 1];

  % horizontal layout for each line
  geomhoriz = {
    [0.8 0.5 0.8 1 1]
    [1]
    [0.7  0.7     2 2 1]
    [0.7 0.5 0.2 2 2 1]
    [0.7 0.5 0.2 2 2 1]
    [1]
    [1 1 1]
    [1]
    [1 1 1]
  };

  % define each individual graphical user element
  elements  = { ...
      { 'style' 'text'       'string' 'Component to fit'    } ...
      { 'style' 'edit'       'string' 'dummy'    'tag' 'component' 'callback' cb_selectcomponent } ...
      { 'style' 'pushbutton' 'string' 'Plot map'                   'callback' cb_plotmap } ...
      { 'style' 'text'       'string' 'Residual variance = '                                             } ...
      { 'style' 'text'       'string' 'dummy'         'tag' 'relvar'                                  } ...
      { } ...
      { 'style' 'text'    'string' 'dipole'     } ...
      { 'style' 'text'    'string' 'fit'        } ...
      { 'style' 'text'    'string' 'position'   } ...
      { 'style' 'text'    'string' 'moment'     } ...
      { } ...
      ...
      { 'style' 'text'        'string' '#1' 'tag' 'dip1'                                 } ...
      { 'style' 'checkbox'    'string' ''   'tag' 'dip1sel'    'callback' cb_checkinput  } { } ...
      { 'style' 'edit'        'string' ''   'tag' 'dip1pos'    'callback' cb_checkinput  } ...
      { 'style' 'edit'        'string' ''   'tag' 'dip1mom'    'callback' cb_checkinput  } ...
      { 'style' 'pushbutton'  'string' 'Flip (in|out)'         'callback' cb_flip1       } ...
      ...
      { 'style' 'text'        'string' '#2' 'tag' 'dip2'                                 } ...
      { 'style' 'checkbox'    'string' ''   'tag' 'dip2sel'    'callback' cb_checkinput  } { } ...
      { 'style' 'edit'        'string' ''   'tag' 'dip2pos'    'callback' cb_checkinput  } ...
      { 'style' 'edit'        'string' ''   'tag' 'dip2mom'    'callback' cb_checkinput  } ...
      { 'style' 'pushbutton'  'string' 'Flip (in|out)'         'callback' cb_flip2       } ...
      ...
      { } { 'style' 'checkbox' 'string' 'Symmetry constrain for dipole #2' 'tag' 'dip2sym' 'callback' cb_sym  'value' 1 } ...
      { } { } { } ...
      { 'style' 'pushbutton'  'string' 'Fit dipole(s)'' position & moment' 'callback' cb_fitposition } ...
      { 'style' 'pushbutton'  'string' 'OR fit only dipole(s)'' moment'   'callback' cb_fitmoment   } ...
      { 'style' 'pushbutton'  'string' 'Plot dipole(s)'                  'callback' cb_plotdip     } ...
    };

  % add the cancel, help and ok buttons at the bottom

  geomvert  = [geomvert 1 1];

  geomhoriz = {geomhoriz{:} [1] [1 1 1]};

  elements  = { elements{:} ...
      { } ...
      { 'Style', 'pushbutton', 'string', 'Cancel', 'callback', cb_close } ...
      { 'Style', 'pushbutton', 'string', 'Help',   'callback', cb_help  } ...
      { 'Style', 'pushbutton', 'string', 'OK',     'callback', cb_ok    } ...
    };

  % activate the graphical interface
  supergui(0, geomhoriz, geomvert, elements{:});
  dlg = gcf;
  set(gcf, 'name', 'Manual dipole fit -- pop_dipfit_nonlinear()');
  set(gcf, 'userdata', EEG);
  pop_dipfit_nonlinear(EEG, 'dialog_setvalue', dlg);
  uiwait(dlg);
  if ishandle(dlg)
    pop_dipfit_nonlinear(EEG, 'dialog_getvalue', dlg);
    % FIXME, rv is undefined since the user may have changed dipole parameters
    % FIXME, see also dialog_getvalue subfucntion
    EEGOUT = get(dlg, 'userdata');
    close(dlg);
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % implement all subfunctions through a switch-yard
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif nargin>=3

  %disp(subfunction)
  EEG = get(parent, 'userdata');

  switch subfunction

    case 'dialog_selectcomponent'
      current = get(findobj(parent, 'tag', 'component'), 'string');
      current = str2num(current);
      current = current(1);
      current = min(current, size(EEG.icaweights,1));
      current = max(current, 1);
      set(findobj(parent, 'tag', 'component'), 'string', int2str(current));
      EEG.dipfit.current = current;
      % reassign the global EEG object back to the dialogs userdata
      set(parent, 'userdata', EEG);
      % redraw the dialog with the current model
      pop_dipfit_nonlinear(EEG, 'dialog_setvalue', parent);

    case 'dialog_plotmap'
      current = str2num(get(findobj(parent, 'tag', 'component'), 'string'));
      figure; pop_topoplot(EEG, 0, current, [ 'IC ' num2str(current) ], [1 1], 1);
      title([ 'IC ' int2str(current) ]);

    case 'dialog_plotcomponent'
      current = get(findobj(parent, 'tag', 'component'), 'string');
      EEG.dipfit.current = str2num(current);
      if ~isempty( EEG.dipfit.current )
        pop_dipplot(EEG, 'DIPFIT',  EEG.dipfit.current, 'normlen', 'on', 'mri', EEG.dipfit.mrifile);
      end;

    case 'dialog_checkinput'
      if get(findobj(parent, 'tag', 'dip1sel'), 'value') & ~get(findobj(parent, 'tag', 'dip1act'), 'value')
        set(findobj(parent, 'tag', 'dip1act'), 'value', 1);
      end
      if get(findobj(parent, 'tag', 'dip2sel'), 'value') & ~get(findobj(parent, 'tag', 'dip2act'), 'value')
        set(findobj(parent, 'tag', 'dip2act'), 'value', 1);
      end
      if ~all(size(str2num(get(findobj(parent, 'tag', 'dip1pos'), 'string')))==[1 3])
        set(findobj(parent, 'tag', 'dip1pos'), 'string', sprintf('%0.3f %0.3f %0.3f', EEG.dipfit.model(EEG.dipfit.current).posxyz(1,:)));
      else
        EEG.dipfit.model(EEG.dipfit.current).posxyz(1,:) = str2num(get(findobj(parent, 'tag', 'dip1pos'), 'string'));
      end
      if ~all(size(str2num(get(findobj(parent, 'tag', 'dip2pos'), 'string')))==[1 3])
        set(findobj(parent, 'tag', 'dip2pos'), 'string', sprintf('%0.3f %0.3f %0.3f', EEG.dipfit.model(EEG.dipfit.current).posxyz(2,:)));
      else
        EEG.dipfit.model(EEG.dipfit.current).posxyz(2,:) = str2num(get(findobj(parent, 'tag', 'dip2pos'), 'string'));
      end
      if ~all(size(str2num(get(findobj(parent, 'tag', 'dip1mom'), 'string')))==[1 3])
        set(findobj(parent, 'tag', 'dip1mom'), 'string', sprintf('%0.3f %0.3f %0.3f', EEG.dipfit.model(EEG.dipfit.current).momxyz(1,:)));
      else
        EEG.dipfit.model(EEG.dipfit.current).momxyz(1,:) = str2num(get(findobj(parent, 'tag', 'dip1mom'), 'string'));
      end
      if ~all(size(str2num(get(findobj(parent, 'tag', 'dip2mom'), 'string')))==[1 3])
        set(findobj(parent, 'tag', 'dip2mom'), 'string', sprintf('%0.3f %0.3f %0.3f', EEG.dipfit.model(EEG.dipfit.current).momxyz(2,:)));
      else
        EEG.dipfit.model(EEG.dipfit.current).momxyz(2,:) = str2num(get(findobj(parent, 'tag', 'dip2mom'), 'string'));
      end
      if get(findobj(parent, 'tag', 'dip2sel'), 'value') & get(findobj(parent, 'tag', 'dip2sym'), 'value') & ~get(findobj(parent, 'tag', 'dip1sel'), 'value')
        set(findobj(parent, 'tag', 'dip2sel'), 'value', 0);
      end
      set(parent, 'userdata', EEG);

    case 'dialog_setvalue'
      % synchronize the gui with the data
      set(findobj(parent, 'tag', 'component'), 'string', EEG.dipfit.current);
      set(findobj(parent, 'tag', 'relvar' ), 'string', sprintf('%0.2f%%', EEG.dipfit.model(EEG.dipfit.current).rv * 100));
      set(findobj(parent, 'tag', 'dip1sel'), 'value', ismember(1, EEG.dipfit.model(EEG.dipfit.current).select));
      set(findobj(parent, 'tag', 'dip2sel'), 'value', ismember(2, EEG.dipfit.model(EEG.dipfit.current).select));
      set(findobj(parent, 'tag', 'dip1pos'), 'string', sprintf('%0.3f %0.3f %0.3f', EEG.dipfit.model(EEG.dipfit.current).posxyz(1,:)));
      if strcmpi(EEG.dipfit.coordformat, 'CTF')
           set(findobj(parent, 'tag', 'dip1mom'), 'string', sprintf('%f %f %f', EEG.dipfit.model(EEG.dipfit.current).momxyz(1,:)));
      else set(findobj(parent, 'tag', 'dip1mom'), 'string', sprintf('%0.3f %0.3f %0.3f', EEG.dipfit.model(EEG.dipfit.current).momxyz(1,:)));
      end;
      Ndipoles = size(EEG.dipfit.model(EEG.dipfit.current).posxyz, 1);
      if Ndipoles>=2
          set(findobj(parent, 'tag', 'dip2pos'), 'string', sprintf('%0.3f %0.3f %0.3f', EEG.dipfit.model(EEG.dipfit.current).posxyz(2,:)));
          if strcmpi(EEG.dipfit.coordformat, 'CTF')
               set(findobj(parent, 'tag', 'dip2mom'), 'string', sprintf('%f %f %f', EEG.dipfit.model(EEG.dipfit.current).momxyz(2,:)));
          else set(findobj(parent, 'tag', 'dip2mom'), 'string', sprintf('%0.3f %0.3f %0.3f', EEG.dipfit.model(EEG.dipfit.current).momxyz(2,:)));
          end;
      end

    case 'dialog_getvalue'
      % synchronize the data with the gui
      if get(findobj(parent, 'tag', 'dip1sel'), 'value'); select = [1]; else select = []; end;
      if get(findobj(parent, 'tag', 'dip2sel'), 'value'); select = [select 2]; end;
      posxyz(1,:) = str2num(get(findobj(parent, 'tag', 'dip1pos'), 'string'));
      posxyz(2,:) = str2num(get(findobj(parent, 'tag', 'dip2pos'), 'string'));
      momxyz(1,:) = str2num(get(findobj(parent, 'tag', 'dip1mom'), 'string'));
      momxyz(2,:) = str2num(get(findobj(parent, 'tag', 'dip2mom'), 'string'));
      % assign the local values to the global EEG object
      EEG.dipfit.model(EEG.dipfit.current).posxyz = posxyz;
      EEG.dipfit.model(EEG.dipfit.current).momxyz = momxyz;
      EEG.dipfit.model(EEG.dipfit.current).select = select;
      % FIXME, rv is undefined after a manual change of parameters
      % FIXME, this should either be undated continuously or upon OK buttonpress
      % EEG.dipfit.model(EEG.dipfit.current).rv = nan;

      % reassign the global EEG object back to the dialogs userdata
      set(parent, 'userdata', EEG);

    case 'dialog_flip'
      % flip the orientation of the dipole
      current = EEG.dipfit.current;
      moment  = EEG.dipfit.model(current).momxyz;
      EEG.dipfit.model(current).momxyz(dipnum,:) = [ -moment(dipnum,1) -moment(dipnum,2) -moment(dipnum,3)];
      set(findobj(parent, 'tag', ['dip' int2str(dipnum) 'mom']), 'string', ...
        sprintf('%0.3f %0.3f %0.3f', EEG.dipfit.model(current).momxyz(dipnum,:)));
      set(parent, 'userdata', EEG);

    case {'dipfit_moment', 'dipfit_position'}
      % determine the selected dipoles and components
      current  = EEG.dipfit.current;
      select   = find([get(findobj(parent, 'tag', 'dip1sel'), 'value') get(findobj(parent, 'tag', 'dip2sel'), 'value')]);
      if isempty(select)
        warning('no dipoles selected for fitting');
        return
      end
      % remove the dipoles from the model that are not selected, but keep
      % the original dipole model (to keep the GUI consistent)
      model_before_fitting = EEG.dipfit.model(current);
      EEG.dipfit.model(current).posxyz = EEG.dipfit.model(current).posxyz(select,:);
      EEG.dipfit.model(current).momxyz = EEG.dipfit.model(current).momxyz(select,:);
      if strcmp(subfunction, 'dipfit_moment')
        % the default is 'yes' which should only be overruled for fitting dipole moment
        cfg.nonlinear  = 'no';
      end
      dipfitdefs;
      if get(findobj(parent, 'tag', 'dip2sym'), 'value') & get(findobj(parent, 'tag', 'dip2sel'), 'value')
          if strcmpi(EEG.dipfit.coordformat,'MNI')
              cfg.symmetry = 'x';
          else
              cfg.symmetry = 'y';
          end;
      else
          cfg.symmetry = [];
      end

      cfg.component  = current;
      % convert structure into list of input arguments
      arg = [fieldnames(cfg)' ; struct2cell(cfg)'];
      arg = arg(:)';

      % make a dialog to interrupt the fitting procedure
      fig = figure('visible', 'off');
      supergui( fig, {1 1}, [], ...
        {'style' 'text' 'string' 'Press button below to stop fitting' }, ...
        {'style' 'pushbutton' 'string' 'Interupt' 'callback' 'figure(gcbf); set(gcbf, ''tag'', ''stop'');' } );
      drawnow;
      % start the dipole fitting
      try
          warning backtrace off;
          EEG = dipfit_nonlinear(EEG, arg{:});
          warning backtrace on;
      catch,
          disp('Dipole localization failed');
      end;

      % should the following string be put into com? ->NOT SUPPORTED
      % --------------------------------------------------------
      com = sprintf('%s = dipfit_nonlinear(%s,%s)\n', inputname(1), inputname(1), vararg2str(arg));

      % this GUI always requires two sources in the dipole model
      % first put the original model back in and then replace the dipole parameters that have been fitted
      model_after_fitting = EEG.dipfit.model(current);
      newfields = fieldnames( EEG.dipfit.model );
      for index = 1:length(newfields)
      	  eval( ['EEG.dipfit.model(' int2str(current) ').' newfields{index} ' = model_after_fitting.' newfields{index} ';' ]);
      end;
      EEG.dipfit.model(current).posxyz(select,:) = model_after_fitting.posxyz;
      EEG.dipfit.model(current).momxyz(select,:) = model_after_fitting.momxyz;
      EEG.dipfit.model(current).rv               = model_after_fitting.rv;
      %EEG.dipfit.model(current).diffmap          = model_after_fitting.diffmap;

      % reassign the global EEG object back to the dialogs userdata
      set(parent, 'userdata', EEG);
      % close the interrupt dialog
      if ishandle(fig)
        close(fig);
      end

    otherwise
      error('unknown subfunction for pop_dipfit_nonlinear');
    end		% switch subfunction

  end		% if nargin

author: PatternRecognition

repo: OpenBMI