text:
% JFE 2003 Universal option valuation using quadrature methods
% Example: an European call option
%
% To use the program to files are needed : ECall_QUAD.m
%
clear all;
% option's contract parameters
 T = 0.5;                     % Maturity
 S0 = 100;                    % Initial stock price
 E = 105;                      % Strike price
 r = 0.06;                    % risk free interest rate
 sig = 0.2;                   % volatility
 D = 0;                       % continuous dividend yield
%%%%%%%%%%%%%%%%%
[OptionValue] = ECall_QUAD(S0,T,E,r,sig,D)
%Black_Scholes_call = blsprice(S0, E, r, T, sig)


text:
function gamefaqs(varargin)

	query = reshape(transpose([strvcat(varargin{:}),[repmat('+',size(strvcat(varargin{:}),1)-1,1);' ']]),1,[]);

	url = ['''http://www.gamefaqs.com/search/index.html?game=',query,''''];

	url(url==' ') = [];

	eval(['web ', url ,' -browser'])

end

author: Sable

repo: mcbench-benchmarks


text:
function value = r4_cosh ( x )

%*****************************************************************************80
%
%% R4_COSH evaluates the hyperbolic cosine of an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the hyperbolic cosine of X.
%
  persistent ymax

  if ( isempty ( ymax ) )
    ymax = 1.0 / sqrt ( r4_mach ( 3 ) );
  end

  y = exp ( abs ( x ) );

  value = 0.5 * y;

  if ( y < ymax )
    value = 0.5 * ( y + 1.0 / y );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% %Example of user input MATLAB file for post processing
% close all
%
% %Plot waves
% waves.plotElevation(simu.rampTime);
% try
%     waves.plotSpectrum();
% catch
% end
%
% % Plot RY forces for body 1
% plotForces(output,1,5)
%
% %Plot RY response for body 1
% output.plotResponse(1,5);
%
% % Plot x forces for body 2
% plotForces(output,2,1)

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function quad = lattice ( dim_num, m, z, f )

%*****************************************************************************80
%
%% LATTICE applies a lattice integration rule.
%
%  Discussion:
%
%    Because this is a standard lattice rule, it is really only suited
%    for functions which are periodic, of period 1, in both X and Y.
%
%    For a suitable F, and a given value of M (the number of lattice points),
%    the performance of the routine is affected by the choice of the
%    generator vector Z.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ian Sloan, Stephen Joe,
%    Lattice Methods for Multiple Integration,
%    Oxford, 1994,
%    ISBN: 0198534728,
%    LC: QA311.S56
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer M, the order (number of points) of the rule.
%
%    Input, integer Z(DIM_NUM), the generator vector.  Typically, the elements
%    of Z satisfy 1 <= Z(1:DIM_NUM) < M, and are relatively prime to M.
%    This is easy to guarantee if M is itself a prime number.
%
%    Input, external real F, the name of the user-supplied routine
%    which evaluates the function, of the form:
%    function f ( dim_num, x )
%    integer dim_num
%    real f
%    real x(dim_num)
%    f = ...
%
%    Output, real QUAD, the estimated integral.
%
  quad = 0.0;

  for j = 0 : m - 1
    x(1:dim_num) = mod ( j * z(1:dim_num) / m, 1.0 );
    quad = quad + f ( dim_num, x );
  end

  quad = quad / m;

  return
end


text:
function calpak_test0125 ( )

%*****************************************************************************80
%
%% CALPAK_TEST0125 tests JED_TO_YMDF_ALEXANDRIAN and YMDF_TO_JED_ALEXANDRIAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST0125\n' );
  fprintf ( 1, '  For the Alexandrian calendar:\n' );
  fprintf ( 1, '  JED_TO_YMDF_ALEXANDRIAN: JED -> YMDF.\n' );
  fprintf ( 1, '  YMDF_TO_JED_ALEXANDRIAN: YMDF -> JED.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  JED (in)    YMDF               JED (out)\n' );
  fprintf ( 1, '\n' );

  jed_epoch = epoch_to_jed_alexandrian ( );

  i = 0;

  while ( 1 )

    i = i + 1;
    jed1 = jed_test ( i );

    if ( jed1 < 0.0 )
      break
    end

    if ( jed_epoch <= jed1 )

      [ y2, m2, d2, f2 ] = jed_to_ymdf_alexandrian ( jed1 );

      s2 = ymd_to_s_alexandrian ( y2, m2, d2 );

      jed3 = ymdf_to_jed_alexandrian ( y2, m2, d2, f2 );

      fprintf ( 1, '  %11.2f     %25s  %11.2f\n', jed1, s2, jed3 );

    end

  end

  return
end


text:
function scale = p13_scale ( neqn )

%*****************************************************************************80
%
%% P13_SCALE returns scale factors for problem p13.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real SCALE(NEQN), the scaling factors.
%
  scale(1:neqn,1) = 1.0;

  return
end


text:
function a = subsasgn(a,s,value)

Q = 'a'; Q2='a';

%% ---- preprocess adding children calls each time there isn't one!!! --
%%  to minimic f.child{1}.child{1}  by just doing f{1}{1} !!

t=s; j=1; isfile=isobject(a); lasttext=0;
for i=1:length(s)

 if ischar(s(i).subs)  %% check if last thing accessed is an object
    isfile = logical(myexist(s(i).subs,'file')==2);
 end

  if equal(s(i).type,'{}') | equal(s(i).type,'()')  %% time for child!
    before=i-1; if i==1 before=1; end; %% look at s(before)
    if ~equal(s(before).subs,'child') & isfile
      % if it exists as file is hopefully an object, isnt child add one!
      t(j)=s(1); t(j).subs='child'; t(j).type='.';
      j=j+1;
    end
  end

  if ischar(s(i).subs)
    lasttext=j;
  end

  t(j)=s(i); j=j+1;
end
s=t;
%%-------------------------------------------------------------

for i=1:length(s)
  if equal(s(i).type, '.')
    Q=[Q, '.', s(i).subs];
  elseif equal(s(i).type, '()')
    tmp = length(s(i).subs);
    Q=[Q, '('];
    for j=1:(tmp-1)
      Q=[Q, 's(', num2str(i), ').subs{', num2str(j), '}, '];
    end;
    Q=[Q, 's(', num2str(i), ').subs{', num2str(tmp), '})'];

  elseif equal(s(i).type, '{}')
    tmp = length(s(i).subs);
    Q=[Q, '{'];
    for j=1:(tmp-1)
      Q=[Q, 's(', num2str(i), ').subs{', num2str(j), '}, '];
    end;
    Q=[Q, 's(', num2str(i), ').subs{', num2str(tmp), '}}'];
  end ;

  if i<lasttext Q2=Q; end;
end;

%% ---- if member doesn't exist search in the children,algorithm.. --------
%%            --- assign value to all possible hits -------

%if ~isfile b=eval(Q2);if isfield(struct(b),name) isfile=1; end; end;
global allow_recursive_subsasgn;

if ~isfile & allow_recursive_subsasgn
  name=s(lasttext).subs;  % this is the member we are searching for
  Qend=Q(length(Q2)+length(name)+2:length(Q));

  stack=[]; stack{1}=Q2;
  while length(stack)>0

      cur=stack{1}; stop=0; b=eval(cur);   % eval top line
      if isfield(struct(b),name)
	Q=[stack{1} '.' name Qend];
  	eval([Q ' = value;']) ;  %% evaluate final answer
	if strcmp(name,'child') stop=1; end; %% shouldnt
                                             %recursively add
                                             %children - it loops
      end; % check for  member

      if ~isfield(struct(b),'alias') & ~stop % as long as not an @algorithm
	ali=b.algorithm.alias;
	for i=1:2:length(ali)
	  if strcmp(name,ali{i})
	    if isfield(struct(b),ali{i+1}) %% found it!
	      newname=ali{i+1};
	      Q=[stack{1} '.' newname Qend];
	      eval([Q ' = value;']) ;  %% evaluate final answer
	    end
	  end
	end
      end

      if isfield(struct(b),'algorithm') & ~stop % check algorithm path
	stack{length(stack)+1}=[ cur '.algorithm'];
      end
      if isfield(struct(b),'child') & ~stop % check child path
	if ~isempty(b.child)
	  if ~iscell(b.child)
	    stack{length(stack)+1}=[ cur '.child'];
	  else
	    for i=1:length(b.child)
	      stack{length(stack)+1}=[ cur '.child{' num2str(i) '}'];
	    end
	  end
	end
      end
      if length(stack)>0
	 stack=stack(2:length(stack));
      end
  end
else
  %% simple evaluation.. don't search for multiple versions..
  Q = [Q, ' = value;'];
  eval(Q) ;  %% evaluate final answer
end

%%-------------------------------------------------------------

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


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

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

  return
end


text:
classdef (Abstract) TagGeneratorX < BaseX
% TagGeneratorX abstract class
%
% Summary of TagGeneratorX:
%   This is the base class for all Tag Generators

    methods (Abstract)
        generate(this); % Generate new tag(s)
%         forget(this);   % Remove tag(s) from list of allocated tags
    end

end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
% Test script for solving the 2D advection
Globals2D;

% Set polynomial order to use
N = 7;

% Read and initiate circular mesh
filename = 'circA01.neu';
[Nv, VX, VY, K, EToV, BCType] = MeshReaderGambitBC2D(filename);

StartUp2D;
BuildBCMaps2D;

% Push all boundary faces to unit cylinder
[k,f] = find(BCType);
curved = sort(unique(k));
MakeCylinder2D([k,f], 1, 0, 0);

% Set initial conditions
% First 6 modes of eigenmodes with 6 azimuthal periods
alpha = [9.936109524217684,13.589290170541217,17.003819667816014,...
         20.320789213566506,23.586084435581391,26.820151983411403];

% choose radial mode
alpha0 = alpha(2);
theta = atan2(y,x);
rad   = sqrt(x.^2+y.^2);

Ez = besselj(6, alpha0*rad).*cos(6*theta);
Hx = zeros(Np, K); Hy = zeros(Np, K);

% Solve Problem for exactly one period
FinalTime = .5;
[Hx,Hy,Ez,time] = MaxwellCurved2D(Hx,Hy,Ez,FinalTime);

exactEz = besselj(6, alpha0*rad).*cos(6*theta)*cos(alpha0*time(end));
maxabserror = max(max(abs(Ez-exactEz)))

repo: Aero-matlab


text:
classdef testHomogenizationLaminateForAnisotropicTensors < ...
         testShowingError

     properties (Access = protected)
         testName = 'HomogenizationLaminateForAnisotropicTensors';
         tol = 1e-6;
     end

     properties (Access = private)
         stiffTensor
         weakTensor
         ChForIso
         ChForAni
         laminateDirection
         theta
     end

     methods (Access = public)

         function obj = testHomogenizationLaminateForAnisotropicTensors()
             obj.init()
             obj.computeHomogenizerForIsotropicMaterials()
             obj.computeHomogenizerForAnisotropicMaterials()
         end

     end

     methods (Access = protected)

         function init(obj)
             obj.createTensors()
             obj.createLaminateDirection()
             obj.theta = 0.8;
         end

        function createLaminateDirection(obj)
            d = [0 1 0];
            dir = Vector3D;
            dir.setValue(d);
            dir.normalize()
            obj.laminateDirection = dir;
        end

         function createTensors(obj)
            E1 = 1;
            E0 = 1e-3;
            nu1 = 1/3;
            nu0 = 1/3;
            obj.stiffTensor = IsotropicConstitutiveTensor(E1,nu1);
            obj.weakTensor  = IsotropicConstitutiveTensor(E0,nu0);
         end

         function computeHomogenizerForIsotropicMaterials(obj)
            C0 = obj.weakTensor;
            C1 = obj.stiffTensor;
            dir{1} = obj.laminateDirection;
            m1 = 1;
            SeqHomog = VoigtPlaneStressHomogHomogenizer(C0,C1,dir,m1,obj.theta);
            obj.ChForIso  = SeqHomog.getPlaneStressHomogenizedTensor();
         end

         function computeHomogenizerForAnisotropicMaterials(obj)
            C0 = obj.weakTensor;
            C1 = obj.stiffTensor;
            dir = obj.laminateDirection;
            Lam = AnisotropicLaminateHomogenizer(C0,C1,dir,obj.theta);
            obj.ChForAni = Lam.getHomogenizedTensor();
         end

         function computeError(obj)
             Ch4Iso = obj.ChForIso.getValue();
             Ch4Ani = obj.ChForAni.getValue();
             obj.error = norm(Ch4Iso - Ch4Ani);
         end

     end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
classdef la

methods(Static)

function no = nonorthogonality(U)
    % Computes the deviation from orthogonality of a basis
    [m, k] = size(U);
    gap = U' * U - eye(k);
    no = norm(gap);
end

end % methods

end % classdef

author: indigits

repo: sparse-plex


text:
  function [grad, out, com] = ir_hct_rgrad(x, varargin)
%|function [grad, out, com] = ir_hct_rgrad(x, varargin)
%|
%| compute 3D regularizer gradient using hct binary
%| for UM testing only
%|
%| in
%|	x	in xyz order
%| option
%|	(many - see below)
%|
%| Notes: because hct2 uses 'zxy' order, it is essential to use Reg1
%| with zxy order also, to ensure gradient matches over entire object.
%| If Reg1 uses the usual xyz order, then slices [1 end] do not match.
%|
%| Jeff Fessler, 2012-06-15

if nargin == 1 && streq(x, 'test'), ir_hct_rgrad_test, return, end
if nargin < 2, ir_usage(), end

arg.R = [];
arg.bij_center = '';
arg.bij_above = '';
arg.log2reg_x = 0; % 0 by default because beta built into R.beta
arg.log2reg_z = []; % ""
arg.kappa_min = 0;
arg.clean = true; % remove files after?
arg.dir = ''; % work directory
arg.file_x = 'ir-hct-rgrad-x.fld'; % work file for x
arg.file_grad = 'ir-hct-rgrad-grad.fld'; % work file for grad
arg.file_mask = 'ir-hct-rgrad-mask.fld'; % work file for mask
arg.file_kappa = 'ir-hct-rgrad-kappa.fld'; % work file for kappa
arg.chat = 0;

arg = vararg_pair(arg, varargin);
if isempty(arg.dir)
	arg.dir = test_dir;
end

R = arg.R;

if isempty(arg.log2reg_z)
	arg.log2reg_z = arg.log2reg_x;
end

arg.file_x = [arg.dir arg.file_x];
arg.file_grad = [arg.dir arg.file_grad];
arg.file_mask = [arg.dir arg.file_mask];
arg.file_kappa = [arg.dir arg.file_kappa];

pn = jf_protected_names;
if ~pn.has_hct2
	fail 'need hct2'
end

if ~streq(R.type_penal, 'zxy') || ~R.offsets_is_zxy
	fail 'zxy required to ensure match'
end

% write files
tox = @(z) permute(z, [2 3 1]);
fld_write(arg.file_x, x)
fld_write(arg.file_mask, tox(R.mask), 'type', 'byte')
fld_write(arg.file_kappa, tox(R.wt.kappa))

if exist(arg.file_grad, 'file')
	eval(['!/bin/rm ' arg.file_grad])
end

if ~isempty(R.pot_params)
	delta = R.pot_params(1);
	qgg_q = R.pot_params(2);
else
	delta = 0;
	qgg_q = 0;
end

if isempty(arg.bij_center)
	if isfield(R, 'beta') && numel(R.beta) > 1
		beta = R.beta;
		displace_zxy = penalty_displace(R.offsets, R.dim);
		oz = displace_zxy(:,1);
		ox = displace_zxy(:,2);
		oy = displace_zxy(:,3);
		bij = zeros(3,3,3);
		ii = sub2ind([3 3 3], 2+ox, 2+oy, 2+oz);
		bij(ii) = beta;
		ii = sub2ind([3 3 3], 2-ox, 2-oy, 2-oz);
		bij(ii) = beta;
		bij = bij / 2^arg.log2reg_x; % scale to match ge style in hct2
		if arg.log2reg_x ~= arg.log2reg_z, fail 'not done', end
		jf_equal(bij(:,:,1), bij(:,:,3)) % verify symmetry
		b = @(i,j) bij(1+i, 1+j, 1); % above
		arg.bij_above = sprintf('%g,%g,%g,%g,%g,%g,%g,%g,%g', ...
                        b(0,0), b(0,1), b(0,2), ...
                        b(1,0), b(1,1), b(1,2), ...
                        b(2,0), b(2,1), b(2,2));
		b = @(i,j) bij(1+i, 1+j, 2); % center
		arg.bij_center = sprintf('%g,%g,%g,%g,%g,%g,%g,%g,%g', ...
                        b(0,0), b(0,1), b(0,2), ...
                        b(1,0), b(1,1), b(1,2), ...
                        b(2,0), b(2,1), b(2,2));
	else
		arg.bij_center = '-';
		arg.bij_above = '-';
	end
end

% trick due to zxy order:
nx = R.dim(2);
ny = R.dim(3);
nz = R.dim(1);

com = sprintf(['hct2 what rgrad chat %d ' ...
	'ns 8 nt 16 na 12 ', ... % fake
	'nx %d ny %d nz %d ' ...
	'bij_center %s bij_above %s ' ...
	'distance_power %g ', ...
	'pot %s delta %g qgg_q %g use_reg_zxy 1 ' ...
	'log2reg_x %.6f log2reg_z %.6f kappa_min %g ' ...
	'file_kappa %s file_mask %s file_init %s file_xh %s'], ...
	...
	arg.chat, nx, ny, nz, ...
	arg.bij_center, arg.bij_above, ...
	R.distance_power, ...
	R.pot_type, delta, qgg_q, ...
	arg.log2reg_x, arg.log2reg_z, arg.kappa_min, ...
	arg.file_kappa, arg.file_mask, arg.file_x, arg.file_grad);

out = os_run(com); % run hct2

grad = fld_read(arg.file_grad);

if arg.clean
	clean = ['!/bin/rm ' arg.file_x ' ' arg.file_grad ' ' arg.file_mask ' ' arg.file_kappa];
	eval(clean)
end

% ir_hct_rgrad_test()
function ir_hct_rgrad_test

rng(0)
kappa = 2 + rand(18,20,10);
%kappa = ones(18,20,10);
kappa([1 end],:,:) = 0; % zero border required for zxy
kappa(:,[1 end],:) = 0;

toz = @(x) permute(x, [3 1 2]);
tox = @(z) permute(z, [2 3 1]);
R = Reg1(toz(kappa), 'beta', 4, 'nthread', jf('ncore'), ...
	'type_penal', 'zxy', ...
	'offsets', '3d:26', ...
	'type_wt', 'fly', ...
	'pot_arg', {'qgg2', 10, 1.2}, ...
	'distance_power', 0);
Rx = Reg1(kappa, 'beta', 4, 'nthread', jf('ncore'), ...
	'offsets', '3d:26', ...
	'type_wt', 'fly', ...
	'pot_arg', {'qgg2', 10, 1.2}, ...
	'distance_power', 0);
%x = zeros(size(kappa));
%x(end/2,end/2,end/2) = 0;
%x(end/2,end/4,1) = 1;
rng(0)
x = randn(size(kappa)) .* (kappa ~= 0);
if 0 % kludge to make penalty value match
	x(:,:,1) = 0;
	x(:,:,end) = 0;
end
mat = R.cgrad(R, toz(x));
mat = tox(mat);

xlim = [-1 1] * 3;
clim = [-1 1] * 1000;
im plc 2 2

im(1, x, xlim)
im(mat, clim)

printm 'calling hct'
%tic
[hct, out] = ir_hct_rgrad(x, 'R', R, 'log2reg_x', 4, 'chat', 99);
%toc
printm(['hct2 output: \n ' out])
pr Rx.penal(Rx, x);

im(hct, clim)
im(hct - mat)
equivs(hct, mat, 'thresh', 2e-6)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [ c1, c2, c3, r1, r2, r3 ] = year_to_scaliger_common ( y )

%*****************************************************************************80
%
%% YEAR_TO_SCALIGER_COMMON converts a Common year to its Scaliger indices.
%
%  Discussion:
%
%    The year 4713 BCE was chosen by Joseph Scaliger for the start of
%    his Julian Ephemeris Date system, because three cycles coincided
%    in that year, the 28 year Julian calendar cycle, the 19 year Metonic
%    cycle, and the 15 year Roman Indiction cycle.  Thus, the year
%    4713 BCE has Scaliger index (1,1,1).  Each subsequent year has a distinct
%    set of Scaliger indices until 7980 years later, when the year
%    3266 CE will again have the Scaliger index (1,1,1).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the Common year.
%
%    Output, integer C1, C2, C3, the number of completed
%    Julian, Metonic and Indiction cycles.
%
%    Output, integer R1, R2, R3, the Julian, Metonic and
%    Indiction cycle numbers that make up the Scaliger index.
%
  if ( y == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'YEAR_TO_SCALIGER_COMMON - Fatal error!\n' );
    fprintf ( 1, '  Illegal input Y = 0.\n' );
    error ( 'YEAR_TO_SCALIGER_COMMON - Fatal error!' );
  end
%
%  Adjust for missing year 0.
%
  if ( y < 0 )
    y2 = y + 1;
  else
    y2 = y;
  end
%
%  Now shift so 4713 BC becomes the year 1.
%
  y2 = y2 + 4713;

  c1 = floor ( ( y2 - 1 ) / 28 );
  c2 = floor ( ( y2 - 1 ) / 19 );
  c3 = floor ( ( y2 - 1 ) / 15 );

  r1 = i4_wrap ( y2, 1, 28 );
  r2 = i4_wrap ( y2, 1, 19 );
  r3 = i4_wrap ( y2, 1, 15 );

  return
end


text:
function bst_headtracking(varargin)
% BST_HEADTRACKING: Displays a subject's head position in realtime; used
% for quality control before recording MEG.
%
% USAGE:    bst_headtracking()              Defaults: isRealtimeAlign=0, hostIP='localhost' and TCPIP=1972
%           bst_headtracking(isRealtimeAlign)
%           bst_headtracking(isRealtimeAlign, hostIP, TCPIP)
%
% Inputs:   isRealtimeAlign = [0,1], 1 turns on realtime alignment with saved headposition
%           hostIP  = IP address of host computer
%           TCPIP   = TCP/IP port of host computer
%
% @=============================================================================
% 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: Elizabeth Bock & Francois Tadel, 2012-2013

global isSaveAlignChannelFile
%% ===== DEFAULT INPUTS ====

% defaults
isRealtimeAlign = 0;
hostIP  = 'localhost';    % IP address of host computer
TCPIP   = 1972;             % TPC/IP port of host computer
PosFile = [];
if nargin == 1
    if strcmp(varargin{1},'RealtimeAlign')
        isRealtimeAlign = 1;
    end
elseif nargin == 2
    isRealtimeAlign = varargin{1};
    hostIP  = varargin{2};
elseif nargin == 3
    isRealtimeAlign = varargin{1};
    hostIP  = varargin{2};
    TCPIP   = varargin{3};
elseif nargin == 4
    isRealtimeAlign = varargin{1};
    hostIP  = varargin{2};
    TCPIP   = varargin{3};
    PosFile = varargin{4};
end

%% ===== CONFIGURATION ====
% % User Directory
% user_dir = bst_get('UserDir');

% Database
ProtocolName  =             'HeadTracking';
SubjectName   =             'HeadMaster';
ConditionHeadPoints =       'HeadPoints'; % Used to warp the head surface
ConditionChan =             'SensorPositions'; % Used to update sensor locations in real-time
ConditionRealtimeAlign =    'RealtimeAlign'; % Used for realtime alignment of previous head position
% Brainstorm
bst_dir =       bst_get('BrainstormHomeDir');
% bst_db_dir =    bst_get('BrainstormDbDir');
tmp_dir =       bst_get('BrainstormTmpDir');

% Initialize FieldTrip
[isInstalled, errMsg, PlugFt] = bst_plugin('Install', 'fieldtrip');
if ~isInstalled
    return;
end
ft_dir = bst_fileparts(which(PlugFt.TestFile));
% Add fieldtrip buffer to path
ft_rtbuffer = bst_fullfile(ft_dir, 'realtime', 'src', 'buffer', 'matlab');
ft_io = bst_fullfile(ft_dir, 'fileio');
if exist(ft_rtbuffer, 'dir') && exist(ft_io, 'dir')
    addpath(ft_rtbuffer);
    addpath(ft_io);
else
    error('Cannot find the FieldTrip buffer and/or io directories');
end
% Newer Fieldtrip version requires multithreading dll.
ft_pthreadlib = bst_fullfile(ft_dir, 'realtime', 'src', 'external', 'pthreads-win64', 'lib');
if ispc && ~exist(bst_fullfile(ft_rtbuffer, 'pthreadGC2-w64.dll'), 'file') && ...
        exist(bst_fullfile(ft_pthreadlib, 'pthreadGC2-w64.dll'), 'file')
    file_copy(bst_fullfile(ft_pthreadlib, 'pthreadGC2-w64.dll'), ...
        bst_fullfile(ft_rtbuffer, 'pthreadGC2-w64.dll'));
end

%% ===== PREPARE DATABASE =====
% Get protocol
iProtocol = bst_get('Protocol', ProtocolName);
% If the protocol doesn't exist yet
if isempty(iProtocol)
    % Create new protocol
    iProtocol = gui_brainstorm('CreateProtocol', ProtocolName, 1, 0); %#ok<NASGU>
else
    % Set as current protocol
    gui_brainstorm('SetCurrentProtocol', iProtocol);
end

%% ===== PREPARE SUBJECT =====
% Get subject
[sSubject, iSubject] = bst_get('Subject', SubjectName);
% If warping: we need to recreate the subject
isWarp = 0;
% Ask subject if the anatomy will be warped
res = java_dialog('question', ...
    'Do you want to warp the anatomy to subject headpoints?', ...
    'Add Headpoints', [], {'Yes', 'No', 'Cancel'});
% User cancelled operation
if isempty(res) || strcmpi(res, 'Cancel')
    return
end

% Warp
if strcmpi(res, 'Yes')
    isWarp = 1;
end

% Create if subject doesnt exist
if isempty(iSubject)
    [sSubject, iSubject] = db_add_subject(SubjectName, [], 1, 0);
end
% Update subject structure
sSubject = bst_get('Subject', iSubject);

if isWarp
    % Delete all the tess and anatomy files in the subject that are not
    % default
    for ii = 1:length(sSubject.Anatomy)
        if isempty(strfind(sSubject.Anatomy(ii).FileName, bst_get('DirDefaultSubject')))
            file_delete(file_fullpath(sSubject.Anatomy(ii).FileName), 1);
        end
    end
    for ii = 1:length(sSubject.Surface)
        if isempty(strfind(sSubject.Surface(ii).FileName, bst_get('DirDefaultSubject')))
            file_delete(file_fullpath(sSubject.Surface(ii).FileName), 1);
        end
    end
end
% Reload subject
db_reload_subjects(iSubject);
% Update subject structure
sSubject = bst_get('Subject', iSubject);

%% ===== PREPARE CONDITION: HEADPOINTS =====
% Get condition
[sStudy, iStudy] = bst_get('StudyWithCondition', [SubjectName '/' ConditionHeadPoints]);
% If warping: we need a channel file in HeadPoints, this will be populated
% with the head points measured and saved in a .pos file
if isWarp
    % Create if condition doesnt exist
    if isempty(sStudy)
        iStudy = db_add_condition(SubjectName, ConditionHeadPoints);
        sStudy = bst_get('Study', iStudy);
    end

    % Copy default channel file to this condition
    DefChannelFile = bst_fullfile(bst_dir, 'defaults', 'meg', 'channel_ctf_default.mat');
    file_copy(DefChannelFile, bst_fileparts(file_fullpath(sStudy.FileName)));
    % Reload condition
    db_reload_studies(iStudy);

    % Get updated study definition
    sStudy = bst_get('Study', iStudy);

    % Update the channel file with the head points
    if isempty(PosFile)
        LastUsedDirs = bst_get('LastUsedDirs');
        PosFile = java_getfile( 'open', 'Select POS file...', LastUsedDirs.ImportChannel, 'single', 'files', ...
            {{'*.pos'}, 'POS files', 'POLHEMUS'}, 0);
        if isempty(PosFile)
            return;
        end
    end
    % Read POS file and channel file
    HeadMat = in_channel_pos(PosFile);
    HPChannelFile = file_fullpath(sStudy.Channel.FileName);
    ChannelMat = in_bst_channel(HPChannelFile);

    % Copy head points
    ChannelMat.HeadPoints = HeadMat.HeadPoints;
    % Force re-alignment on the new set of NAS/LPA/RPA (switch from CTF coil-based to SCS anatomical-based coordinate system)
    ChannelMat = channel_detect_type(ChannelMat, 1, 0);
    save(HPChannelFile, '-struct', 'ChannelMat');

else
    HPChannelFile = file_fullpath(sStudy.Channel.FileName);
    ChannelMat = in_bst_channel(HPChannelFile);
end

% Get the transformation for HPI head coordinates (POS file) to Brainstorm
iTrans = find(~cellfun(@isempty,strfind(ChannelMat.TransfMegLabels, 'Native=>Brainstorm/CTF')));
if isempty(iTrans)
    bst_error('No SCS transformation in the channel file')
    return;
end
%ChannelMat.TransfMeg{end+1} = [ChannelMat.SCS.R, ChannelMat.SCS.T; 0 0 0 1];
trans = ChannelMat.TransfMeg{iTrans};
R = trans(1:3, 1:3);
T = trans(1:3, 4); %in meters

%% ===== WARP =====
if isWarp
    % Warp surface for new head points
    % bst_warp_prepare(ChannelFile, Options)
    %    Options     : Structure of the options
    %         |- tolerance    : Percentage of outliers head points, ignored in the calulation of the deformation.
    %         |                 Set to more than 0 when you know your head points have some outliers.
    %         |                 If not specified: asked to the user (default
    %         |- isInterp     : If 0, do not do a full interpolation (default: 1)
    %         |- isScaleOnly  : If 1, do not perform the full deformation but only a linear scaling in the three directions (default: 0)
    %         |- isSurfaceOnly: If 1, do not warp/scale the MRI (default: 0)
    %         |- isInteractive: If 0, do not ask anything to the user (default: 1)
    Options.tolerance = 0.02;
    Options.isInterp = [];
    Options.isScaleOnly = 0;
    Options.isSurfaceOnly = 1;
    Options.isInteractive = 1;
    hFig = bst_warp_prepare(HPChannelFile, Options);

    % Close figure
    close(hFig);
    % Update subject structure
    sSubject = bst_get('Subject', iSubject);
end

%% ===== PREPARE CONDITION: CHANNEL POSITIONS =====
% Get condition
[sStudyChan, iStudyChan] = bst_get('StudyWithCondition', [SubjectName '/' ConditionChan]);
% Create if condition doesnt exist, this will be populated with sensor
% positions measured from real-time res4 file
if isempty(sStudyChan)
    iStudyChan = db_add_condition(SubjectName, ConditionChan);
    sStudyChan = bst_get('Study', iStudyChan);
end

%% ===== INITIALIZE FIELDTRIP BUFFER =====
% find and kill any old matlab processes
pid = feature('getpid');
[tmp,tasks] = system('tasklist');
pat = '\s+';
str=tasks;
s = regexp(str, pat, 'split');
iMatlab = find(~cellfun(@isempty, strfind(s, 'MATLAB.exe')));

if numel(iMatlab) > 1
    % kill the extra matlab(s)
    disp('An old MATLAB process is still running, stopping now...');
    temp = str2double(s(iMatlab+1));
    iKill = find(temp ~= pid);
    for ii=1:length(iKill)
        [status,result] = system(['Taskkill /PID ' num2str(temp(iKill(ii))) ' /F']);
        if status
            disp(result);
            disp('The process could not be stopped.  Please stop it manually');
        else
            disp(result);
        end
    end
end

% Initialize the buffer
try
    disp('BST> Initializing FieldTrip buffer...');
    buffer('tcpserver', 'init', hostIP, TCPIP);
catch
    disp('BST> Warning: FieldTrip buffer is already initialized.');
    buffer('flush_hdr', [], hostIP, TCPIP);
end
% disp('BST> Remember to type ''clear buffer'' once done with head tracking.'); % Doesn't work: freezes.
% Waiting for acquisition to start
disp('BST> Waiting for acquisition to start...');
while (1) % To fix: This loop can make Matlab crash if wait is too long.
    try
        % [nsamples, nevents, timeout]
        numbers = buffer('wait_dat', [1 -1 1000], hostIP, TCPIP);
        disp('BST> Acquisition started.');
        break;
    catch
        pause(1);
    end
end

%% ===== READING RES4 INFO =====
% Reading header
hdr = buffer('get_hdr', [], hostIP, TCPIP);

% Write .res4 file
res4_file = bst_fullfile(tmp_dir, 'temp.res4');
fid = fopen(res4_file, 'w', 'l');
fwrite(fid, hdr.ctf_res4, 'uint8');
fclose(fid);
% Write empty .meg4
meg4_file = bst_fullfile(tmp_dir, 'temp.meg4');
fid = fopen(meg4_file, 'w', 'l');
fclose(fid);
% Reading structured res4
SensorPositionMat = in_channel_ctf(res4_file);
% Add channel file to the SensorPositions study
SensorPositionFile = db_set_channel(iStudyChan, SensorPositionMat, 2, 2);

%% ===== HEAD TRACKING =====
% Display subject's head
hFig = view_surface(sSubject.Surface(sSubject.iScalp).FileName);
% Set view from the left
figure_3d('SetStandardView', hFig, 'front');
% Check for RealtimeAlign
if isRealtimeAlign
    % Create save button
    btn = uicontrol('Style', 'pushbutton', 'String', 'Save',...
        'Position', [20 20 50 20],...
        'Callback', @SaveAlignCallback);
    [sStudyAlign, iStudyAlign] = bst_get('StudyWithCondition', [SubjectName '/' ConditionRealtimeAlign]);
else
    sStudyAlign = [];
end
ColorTable = [1,0,0; 0,1,0; 0,0,1];
colorInd = 2;
if ~isempty(sStudyAlign)
    % Display CTF helmet
    view_helmet_local(sStudyAlign.Channel.FileName, hFig); % _local
    % Get the helmet patch
    hHelmetPatch = findobj(hFig, 'Tag', 'HelmetPatch');
    color = ColorTable(mod(colorInd-1,size(ColorTable,1))+1,:);
    set(hHelmetPatch, 'FaceColor', color, 'FaceAlpha', .3, 'SpecularStrength', 0, ...
                            'EdgeColor', color, 'EdgeAlpha', .2, 'LineWidth', 1, ...
                            'Marker', 'none', 'Tag', 'MultipleSensorsPatches');
    AlignVertices = get(hHelmetPatch, 'Vertices');
end

% Display current position helmet
colorInd = 1;
% Display CTF helmet
view_helmet_local(SensorPositionFile, hFig);
% Get the helmet patch
hHelmetPatch = findobj(hFig, 'Tag', 'HelmetPatch');
color = ColorTable(mod(colorInd-1,size(ColorTable,1))+1,:);
set(hHelmetPatch, 'FaceColor', color, 'FaceAlpha', .3, 'SpecularStrength', 0, ...
                        'EdgeColor', color, 'EdgeAlpha', .2, 'LineWidth', 1, ...
                        'Marker', 'none', 'Tag', 'MultipleSensorsPatches', ...
                        'Visible', 'on');
% Get XYZ coordinates of the helmet patch object
InitVertices = get(hHelmetPatch, 'Vertices');

% Loop to update positions
while (1)
    % Number of samples to read
    nSamples = 300;
    % Read the last nSamples fiducial positions
    hdr = buffer('get_hdr', [], hostIP, TCPIP);
    currentSample = hdr.nsamples;
    if currentSample < nSamples-1
        % Wait
        pause(.2);
        continue;
    end
    dat = buffer('get_dat', [currentSample-nSamples-1,currentSample-1], hostIP, TCPIP);
    % Average in time
    buf = mean(double(dat.buf), 2);
    Fid = [buf(1:3), buf(4:6), buf(7:9)] ./ 1e3; % convert to mm
    % Get fiducial positions
    sMri.SCS.NAS = Fid(:,1)';
    sMri.SCS.LPA = Fid(:,2)';
    sMri.SCS.RPA = Fid(:,3)';
    % Compute transformation Dewar=>Native (CTF head coils)
    transfSCS = cs_compute(sMri, 'scs'); % NAS, LPA and RPA in mm
    sMri.SCS.R = transfSCS.R;
    sMri.SCS.T = transfSCS.T;
    sMri.SCS.Origin = transfSCS.Origin;
    % Apply transformation to helmet vertices.  sMri is in mm and vertices
    % are in m (as it should be according to cs_convert code).
    Vertices = cs_convert(sMri, 'mri', 'scs', InitVertices);
    % Convert HPI coordinates to Brainstorm coordinates (based on cardinal points)
    Vertices = bst_bsxfun(@plus, R * Vertices', T)';
    % Stop if the window was closed
    if ~ishandle(hHelmetPatch)
        break;
    end
    % Update helmet patch
    set(hHelmetPatch, 'Vertices', Vertices);

    % check difference between current and align position
    if ~isempty(sStudyAlign)
        maxDiff = (max(max(100*abs(AlignVertices - Vertices)/AlignVertices)));

        if maxDiff < 10
            colorInd = 2;
            color = ColorTable(mod(colorInd-1,size(ColorTable,1))+1,:);
            set(hHelmetPatch, 'FaceColor', color, 'FaceAlpha', .3, 'SpecularStrength', 0, ...
                                'EdgeColor', color, 'EdgeAlpha', .2, 'LineWidth', 1, ...
                                'Marker', 'none', 'Tag', 'MultipleSensorsPatches');
        else
            colorInd = 1;
            color = ColorTable(mod(colorInd-1,size(ColorTable,1))+1,:);
            set(hHelmetPatch, 'FaceColor', color, 'FaceAlpha', .3, 'SpecularStrength', 0, ...
                                'EdgeColor', color, 'EdgeAlpha', .2, 'LineWidth', 1, ...
                                'Marker', 'none', 'Tag', 'MultipleSensorsPatches');
        end
    end

    if isSaveAlignChannelFile
        % Check for RealtimeAlign condition
        [tmp, iStudyAlign] = bst_get('StudyWithCondition', [SubjectName '/' ConditionRealtimeAlign]);
        if isempty(iStudyAlign)
            % create the condition and channel file with the current sensor positions for subsequent alignment
            iStudyAlign = db_add_condition(SubjectName, ConditionRealtimeAlign);
        end
        SensorPositionMat.SCS = sMri.SCS; % in mm
        db_set_channel(iStudyAlign, SensorPositionMat, 2, 0);
        SaveAlignChannelFile(HPChannelFile, iStudyAlign, SensorPositionMat);
    end
    % Wait
    pause(.2);
end
end
%% SaveAlignCallback
function SaveAlignCallback(source,callbackdata)
    global isSaveAlignChannelFile
    isSaveAlignChannelFile = 1;
end
%% Head Localization
function SaveAlignChannelFile(HPChannelFile, iStudyAlign, ChannelMat)
    global isSaveAlignChannelFile
    % ChannelMat.SCS already contains the transformation Dewar=>Native (CTF
    % head coils), in mm

    % Get transformation Native=>Brainstorm/CTF (anatomical NAS/LPA/RPA)
    % from head points file.
    HPChannelMat = in_bst_channel(HPChannelFile);
    iTrans = find(~cellfun(@isempty,strfind(HPChannelMat.TransfMegLabels, 'Native=>Brainstorm/CTF')));
    if isempty(iTrans)
        bst_error('No SCS transformation in the head points channel file')
        return;
    end
    % Get the translation and rotation from the HeadPoints tranformation
    trans = HPChannelMat.TransfMeg{iTrans};
    anatR = trans(1:3, 1:3);
    anatT = trans(1:3, 4) * 1000; % convert from m to mm

    % Combine the tranformations
    transfAnat = [anatR, anatT; 0 0 0 1]*[ChannelMat.SCS.R, ChannelMat.SCS.T; 0 0 0 1]; % in mm

    % Update the ChannelMat structure
    ChannelMat.SCS.R = transfAnat(1:3, 1:3);
    ChannelMat.SCS.T = transfAnat(1:3, 4); % in mm

    % Process each sensor
    for i = 1:length(ChannelMat.Channel)
        if ~isempty(ChannelMat.Channel(i).Loc)
            % Converts the locations. ChannelMat.SCS is in mm and locations
            % are in m (as it should be according to cs_convert code).
            ChannelMat.Channel(i).Loc = cs_convert(ChannelMat, 'mri', 'scs', ChannelMat.Channel(i).Loc')';
        end
    end

    % Also transform the fiducial positions, temporarily converted to m.
    ChannelMat.SCS.NAS = cs_convert(ChannelMat, 'mri', 'scs', ChannelMat.SCS.NAS ./ 1000) .* 1000; % still in mm
    ChannelMat.SCS.LPA = cs_convert(ChannelMat, 'mri', 'scs', ChannelMat.SCS.LPA ./ 1000) .* 1000;
    ChannelMat.SCS.RPA = cs_convert(ChannelMat, 'mri', 'scs', ChannelMat.SCS.RPA ./ 1000) .* 1000;

    % Update the list of transformations.  Here the translation must be in
    % m.
    if isempty(ChannelMat.TransfMegLabels)
        iTrans = [];
    else
        iTrans = find(~cellfun(@isempty,strfind(ChannelMat.TransfMegLabels, 'Native=>Brainstorm/CTF')));
    end

    if isempty(iTrans)
        ChannelMat.TransfMeg{end+1} = [ChannelMat.SCS.R, ChannelMat.SCS.T ./ 1000; 0 0 0 1];
        ChannelMat.TransfMegLabels{end+1} = 'Native=>Brainstorm/CTF';
        ChannelMat.TransfEeg{end+1} = [ChannelMat.SCS.R, ChannelMat.SCS.T ./ 1000; 0 0 0 1];
        ChannelMat.TransfEegLabels{end+1} = 'Native=>Brainstorm/CTF';
    else
        ChannelMat.TransfMeg{iTrans} = [ChannelMat.SCS.R, ChannelMat.SCS.T ./ 1000; 0 0 0 1];
        ChannelMat.TransfMegLabels{iTrans} = 'Native=>Brainstorm/CTF';
        ChannelMat.TransfEeg{iTrans} = [ChannelMat.SCS.R, ChannelMat.SCS.T ./ 1000; 0 0 0 1];
        ChannelMat.TransfEegLabels{iTrans} = 'Native=>Brainstorm/CTF';
    end

    % Save new channel file to the target studies.
    db_set_channel(iStudyAlign, ChannelMat, 2, 0);
    isSaveAlignChannelFile = 0;

end

function view_helmet_local(ChannelFile, hFig)
% Local trimmed down version of view_helmet, such that it also uses
% channel_tesselate_local which was modified to work in dewar coordinates.
% The only change, at the end of this function, is that the sensor patch is
% deleted instead of making it invisible.

Modality = 'MEG';

% hFig = view_channels(ChannelFile, Modality, 1, 0, hFig);
  [hFig] = bst_figures('GetFigure', hFig);
  %   iDS = bst_memory('GetDataSetChannel', ChannelFile);
  %   if isempty(hFig)
  %       % Prepare FigureId structure
  %       FigureId = db_template('FigureId');
  %       FigureId.Type     = '3DViz';
  %       FigureId.SubType  = '';
  %       FigureId.Modality = Modality;
  %       % Create figure
  %       [hFig, iFig, isNewFig] = bst_figures('CreateFigure', iDS, FigureId);
  %       % If figure was not created: Display an error message and return
  %       if isempty(hFig)
  %           bst_error('Cannot create figure', '3D figure creation...', 0);
  %           return;
  %       end
  %   end

%   figure_3d('ViewSensors', hFig, 1, 0, 1, Modality);
    % Load channel file
    ChannelMat = in_bst_channel(ChannelFile);
    Channel = ChannelMat.Channel;
    selChan = good_channel(Channel, [], Modality);
    % Get sensors positions

%     [tmp, markersLocs] = GetChannelPositions(iDS, selChan);
    % Initialize returned variables
    vertices    = zeros(3,0);
    % Get device type
    Device = bst_get('ChannelDevice', ChannelFile);
    % Get selected channels
    Channel = Channel(selChan);
    % Find magnetometers
    if strcmpi(Device, 'Vectorview306')
        iMag = good_channel(Channel, [], 'MEG MAG');
    end
    % Loop on all the sensors
    for i = 1:length(Channel)
        % If position is not defined
        if isempty(Channel(i).Loc)
            Channel(i).Loc = [0;0;0];
        end
        % Get number of integration points or coils
        nIntegPoints = size(Channel(i).Loc, 2);
        % Switch depending on the device
        switch (Device)
            case {'CTF', '4D', 'KIT', 'RICOH'}
                if (nIntegPoints >= 4)
                    vertices    = [vertices,    Channel(i).Loc(:,1:4)]; %#ok<*AGROW>
                else
                    vertices    = [vertices,    Channel(i).Loc];
                end
            case 'KRISS'
                if (nIntegPoints >= 4)
                    vertices    = [vertices,    Channel(i).Loc(:,1:4)];
                else
                    vertices    = [vertices,    Channel(i).Loc(:,1)];
                end
            case 'Vectorview306'
                if isempty(iMag) || ismember(i, iMag)
                    vertices = [vertices, Channel(i).Loc];
                end
            case 'BabySQUID'
                vertices    = [vertices,    Channel(i).Loc(:,1)];
            case 'BabyMEG'
                vertices    = [vertices,    Channel(i).Loc];
            case {'NIRS-BRS', 'NIRS'}
                Factor = 1;
                % Position of the channel: mid-way between source and detector, organized in layers by wavelength
                vertices    = [vertices,    mean(Channel(i).Loc,2) .* Factor];
            otherwise
                vertices    = [vertices,    Channel(i).Loc];
        end
    end
    vertices    = double(vertices');
% End of content from GetChannelPositions subfunction.

    % Put focus on target figure
    hAxes = findobj(hFig, '-depth', 1, 'Tag', 'Axes3D');

%     PlotSensorsNet(hAxes, markersLocs, 1, isMesh, markersOrient);
      faces = channel_tesselate_local( vertices );
      % === DISPLAY PATCH ===
      % Display faces / edges / vertices
      FaceColor = [.7 .7 .5];
      EdgeColor = [.4 .4 .3];
      LineWidth = 1;
      % Create sensors patch
      hNet = patch(...
        'Vertices',         vertices, ...
        'Faces',            faces, ...
        'FaceVertexCData',  repmat([1 1 1], [length(vertices), 1]), ...
        'Parent',           hAxes, ...
        'Marker',           'o', ...
        'LineWidth',        LineWidth, ...
        'FaceColor',        FaceColor, ...
        'FaceAlpha',        1, ...
        'EdgeColor',        EdgeColor, ...
        'EdgeAlpha',        1, ...
        'MarkerEdgeColor',  [.4 .4 .3], ...
        'MarkerFaceColor',  'flat', ...
        'MarkerSize',       6, ...
        'BackfaceLighting', 'lit', ...
        'AmbientStrength',  0.5, ...
        'DiffuseStrength',  0.5, ...
        'SpecularStrength', 0.2, ...
        'SpecularExponent', 1, ...
        'SpecularColorReflectance', 0.5, ...
        'FaceLighting',     'gouraud', ...
        'EdgeLighting',     'gouraud', ...
        'Tag',              'SensorsPatch');
% End of content from figure_3d.

  % Update lights
  camlight(findobj(hFig, 'Tag', 'FrontLight'), 'headlight');
  % Update figure name
  bst_figures('UpdateFigureName', hFig);
  % Camera basic orientation
  %   if isNewFig
  %       figure_3d('SetStandardView', hFig, 'top');
  %   end
  % Show figure
  set(hFig, 'Visible', 'on');
% End of content from view_channels.

% Get sensors patch
hSensorsPatch = findobj(hFig, 'Tag', 'SensorsPatch');
if isempty(hSensorsPatch)
    return
elseif numel(hSensorsPatch) > 1
    hSensorsPatch = hSensorsPatch(1);
end
% Get sensors positions
vert = get(hSensorsPatch, 'Vertices');

% ===== CREATE HELMET SURFACE =====
% Get the acquisition device
Device = bst_get('ChannelDevice', ChannelFile);
% Distance sensors/helmet
switch (Device)
    case 'Vectorview306'
        dist = .019;
    case 'CTF'
        dist = .015;
    case '4D'
        dist = .015;
    case 'KIT'
        dist = .020;
    case 'KRISS'
        dist = .025;
    case 'BabyMEG'
        dist = .008;
    case 'RICOH'
        dist = .020;
    otherwise
        dist = 0;
end
% Shrink sensor patch to create the inner helmet surface
if (dist > 0)
    center = mean(vert);
    vert = bst_bsxfun(@minus, vert, center);
    [th,phi,r] = cart2sph(vert(:,1),vert(:,2),vert(:,3));
    [vert(:,1),vert(:,2),vert(:,3)] = sph2cart(th, phi, r - dist);
    vert = bst_bsxfun(@plus, vert, center);
end

% ===== DISPLAY HELMET SURFACE =====
% Copy sensor patch object
hHelmetPatch = copyobj(hSensorsPatch, get(hSensorsPatch,'Parent'));
% Make the sensor patch invisible
% set(hSensorsPatch, 'Visible', 'off');
delete(hSensorsPatch);
% Set patch properties
set(hHelmetPatch, 'Vertices',   vert, ...
                   'LineWidth',  1, ...
                   'EdgeColor',  [.5 .5 .5], ...
                   'EdgeAlpha',  1, ...
                   'FaceColor',  'y', ...
                   'FaceAlpha',  .3, ...
                   'Marker',     'none', ...
                   'Tag',        'HelmetPatch');

end

function Faces = channel_tesselate_local( Vertices, isPerimThresh )
% Modified copy of channel_tesselate to work in dewar coordinates:
% 2d projection changed to "centered - max(z)" coordinates and
% thresholdPerim increased to 8 std instead of 6 to avoid helmet holes.

% CHANNEL_TESSELATE: Tesselate a set of EEG or MEG sensors, for display purpose only.
%
% USAGE:  Faces = channel_tesselate( Vertices, isPerimThresh=1 )
%
% INPUT:
%    - Vertices      : [Nx3], set of 3D points (MEG or EEG sensors)
%    - isPerimThresh : If 1, remove the Faces that are too big
% OUTPUT:
%    - Faces    : [Mx3], result of the tesselation

% Parse inputs
if (nargin < 2) || isempty(isPerimThresh)
    isPerimThresh = 1;
end

% === TESSELATE ===
% Compute best fitting sphere
bfs_center = bst_bfs(Vertices)';
% Center Vertices on BFS center
coordC = bst_bsxfun(@minus, Vertices, bfs_center);
    % 2D Projection
    % Use centered points in case points are not in SCS coordinates, which
    % is the case e.g. when doing real-time head position display (dewar
    % coordinates).
    [X,Y] = bst_project_2d(coordC(:,1), coordC(:,2), coordC(:,3)-max(coordC(:,3)), '2dcap');
% Normalize coordinates
coordC = bst_bsxfun(@rdivide, coordC, sqrt(sum(coordC.^2,2)));
% Tesselation of the sensor array
Faces = convhulln(coordC);

% === REMOVE UNNECESSARY TRIANGLES ===
% For instance: the holes for the ears on high-density EEG caps
if isPerimThresh
    % Get border of the representation
    border = convhull(X,Y);
    % Keep Faces inside the border
    iInside = find(~(ismember(Faces(:,1),border) & ismember(Faces(:,2),border)& ismember(Faces(:,3),border)));
    %Faces   = Faces(iInside, :);

    % Compute perimeter
    triPerimeter = tess_perimeter(Vertices, Faces);
    % Threshold values
    thresholdPerim = mean(triPerimeter(iInside)) + 8 * std(triPerimeter(iInside));
    % Apply threshold
    iFacesOk = intersect(find(triPerimeter <= thresholdPerim), iInside);
    % Find Vertices that are not in the Faces matrix
    iVertNotInFaces = setdiff(1:length(Vertices), unique(Faces(:)));
    if ~isempty(iVertNotInFaces)
        disp(['CHANNEL_TESSELATE> WARNING: Some sensors are not in the Faces list: ' sprintf('%d ', iVertNotInFaces)]);
    end
    % Loop until all the Vertices are visible
    isMissing = 1;
    while isMissing
        % List all the Vertices ignored by the reduced mesh
        iVertOk = unique(reshape(Faces(iFacesOk,:),[],1));
        iVertMissing = setdiff(1:length(Vertices), iVertOk);
        iVertMissing = setdiff(iVertMissing, iVertNotInFaces);
        % If all the Vertices are included, next step
        if isempty(iVertMissing)
            isMissing = 0;
        else
            % Find Faces connected to the first missing vertex
            iFacesAdd = find(any(Faces == iVertMissing(1), 2));
            % From the potential candidate Faces, keep the one that has the smaller perimeter
            [minP, iMinP] = min(triPerimeter(iFacesAdd));
            % Add the smallest face to the list of Faces we keep
            iFacesOk(end+1) = iFacesAdd(iMinP);
        end
    end
    % Remove the Faces
    Faces = Faces(sort(iFacesOk),:);
end

end

author: brainstorm-tools


text:
function bivand2_test ( )

%*****************************************************************************80
%
%% BIVAND2_TEST tests BIVAND2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 May 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BIVAND2_TEST:\n' );
  fprintf ( 1, '  Compute a bidimensional Vandermonde matrix\n' );
  fprintf ( 1, '  associated with the product polynomials of\n' );
  fprintf ( 1, '  maximum degree less than N.\n' );

  n = 3;
  nn = n^2;

  alpha = [ 1.0; 2.0; 3.0 ];
  beta = [ 10.0; 20.0; 30.0 ];

  r8vec_print ( n, alpha, '  Vandermonde vector ALPHA:' );
  r8vec_print ( n, beta, '  Vandermonde vector BETA:' );

  a = bivand2 ( n, alpha, beta );

  r8mat_print ( nn, nn, a, '  Bidimensional Vandermonde matrix:' );

  return
end


text:
function [ptrace, ptrace_n,...
    idx_sort1, idx_sort2, idx_match1, idx_match2] = set_matching_figure(cgk1, cgk2, flag_figure)

% assume cgk1,2 already taken intersection
if ~exist('flag_figure', 'var')
    flag_figure = 1;
end

% set matching/covering index
M = set_match_count(cgk1, cgk2);
[n1, n2] = size(M);
din = sum(M,2);
din(din==0) = 1;
%         dout = sum(M,1);
%         dout(dout==0) = 1;
M1 = M./ repmat(din, [1, n2]);
% M1b = M./ repmat(dout, [n1, 1]);

% iterative sorting...
Msort = M;
idx1 = 1:n1;
idx2 = 1:n2;
flag_converge = 0;
for k = 1: 1000
    [~, idsort2] = sortMn_column(Msort);
    Msort = Msort(:, idsort2);
    idx2 = idx2(idsort2);
    [~, idsort1] = sortMn_column(Msort');
    Msort = Msort(idsort1, :);
    idx1 = idx1(idsort1);

    if sum(abs(idsort2-(1:n2)))<=2 && sum(abs(idsort1-(1:n1))) <= 2
        flag_converge = 1;
        break;
    end
end
if flag_converge == 0
    display(' sorting does not converge');
end

M2 = M1(idx1, idx2);
M2n = M(idx1, idx2);

if flag_figure == 1
    fig1 = figure;
    subplot(121);
    imagesc(M2);
    colormap(bluewhitered);
end

% matching
[assignment, ~] = munkres(-M2n); % based on count
idx3 = 1:n1;
M3 = M2(idx3(assignment>0), assignment(assignment>0));
M3n = M2n(idx3(assignment>0), assignment(assignment>0));

ptrace = trace(M3)/sum(sum(M3));
ptrace_n = trace(M3n)/sum(sum(M3n));

if flag_figure == 1
    figure(fig1);
    subplot(122);
    imagesc(M3);
    colormap(bluewhitered);
    axis equal
    axis tight
    % title([num2str(ptrace), ', ',...
    %     num2str(ptrace_n), ', ', num2str(size(M3,1))])
end

idx_sort1 = idx1;
idx_sort2 = idx2;
idx_match1 = idx1(idx3(assignment>0));
idx_match2 = idx2(assignment(assignment>0));
end

author: xiuyechen

repo: FishExplorer


text:
%FEAT2OBJ Transform feature images to object images in dataset
%
%   B = FEAT2OBJ(A)
%
% INPUT
%   A     Dataset with object images, possible with multiple bands
%
% OUTPUT
%   B     Dataset with features images
%
% SEE ALSO (PRTools Guide)
% DATASETS, IM2OBJ, IM2FEAT, DATA2IM, OBJ2FEAT

function b = feat2obj(a)

  isdataset(a)
  isfeatim(a);
  im = data2im(a);
  b = im2feat(im);

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function minLoc = overlaysurface(fitnessfcn,options)
% Overlays surface and marks known minimum for demo cases of pso. Called by
% PSOPLOTSWARMSURF.
%
% See also: PSODEMO, PSOPLOTSWARMSURF

xBound = options.PopInitRange(:,1) ;
yBound = options.PopInitRange(:,2) ;
[XX,YY] = meshgrid(xBound(1):(xBound(2)-xBound(1))/50:xBound(2),...
    yBound(1):(yBound(2)-yBound(1))/50:yBound(2)) ;
ZZ = zeros(size(XX)) ;
clear xBound yBound

for i = 1:size(XX,1)
    for j = 1:size(XX,2)
        ZZ(i,j) = fitnessfcn([XX(i,j) YY(i,j)]) ;
    end
end

if strcmpi(options.DemoMode,'pretty') || strcmpi(options.DemoMode,'on')
    surface(XX,YY,ZZ,'LineStyle','none',...
        'FaceAlpha',0.4,...
        'FaceLighting','gouraud',...
        'FaceColor','interp')
    flagfill = [1 0 0] ;
    flaglines = [1 0 0] ;
    set(gcf,'Colormap',1-colormap('cool'))
else
    surface(XX,YY,ZZ,'LineStyle','-',...
        'FaceColor','none',...
        'EdgeColor',[1 0.5 0.5])
    flagfill = 'none' ;
    flaglines = [0 0.5 0] ;
%     set(gca,'Color','none')
%     set(gcf,'Colormap',colormap('summer'))
end
axis tight

% Mark known minimum location, if it exists
if ~isempty(options.KnownMin)
    xmin = options.KnownMin(1) ;
    ymin = options.KnownMin(2) ;
    zminmax = get(gca,'ZLim') ;
    minLoc = line([xmin xmin],[ymin ymin],zminmax,...
        'Color',flaglines,...
        'LineWidth',1) ;
    aspr = get(gca,'DataAspectRatio') ;
    flagtip = 0.125*diff(zminmax)*aspr(1)/aspr(3) ;
    patch([xmin xmin xmin + flagtip],...
        [ymin ymin ymin],...
        [zminmax(2),...
         zminmax(2) + 0.1*diff(zminmax),...
         zminmax(2) + 0.05*diff(zminmax)],[1 0 0],...
         'FaceColor',flagfill,...
         'LineWidth',1,...
         'EdgeColor',flaglines)
     clear xmin ymin
end

axis tight
camlight
set(gca,...
    'XLimMode','manual',...
    'YLimMode','manual',...
    'ZLimMode','manual',...
    'PlotBoxAspectRatio',[1 1 1])

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function objNew = map2ecefTrafo(obj, mstruct, varargin)
% MAP2ECEFTRAFO Transform exterior orientation from ecef to map coordinates.

% Input parsing ----------------------------------------------------------------

p = inputParser;
p.addRequired('mstruct');
p.parse(mstruct);
p = p.Results;
% Clear required inputs to avoid confusion
clear mstruct

% Start ------------------------------------------------------------------------

% procHierarchy = {'IMG' 'MAP2ECEFTRAFO'};
% msg('S', procHierarchy);
% msg('I', procHierarchy, sprintf('Image label = ''%s''', obj.label));

% Conversion of position to ecef coordinates -----------------------------------

X0Map = obj.X0;
Y0Map = obj.Y0;
Z0Map = obj.Z0;

[lat, lon, hEll] = minvtran(p.mstruct, X0Map, Y0Map, Z0Map); % lat, lon in degrees!
[X0Ecef, Y0Ecef, Z0Ecef] = geodetic2ecef(lat*pi/180, lon*pi/180, hEll, p.mstruct.geoid); % lat, lon in radian!

obj.X0 = X0Ecef;
obj.Y0 = Y0Ecef;
obj.Z0 = Z0Ecef;

% Conversion of angles to ecef coordinates -------------------------------------

RMap = obj.R;

xUnitVectorMap = RMap(:,1);
yUnitVectorMap = RMap(:,2);
zUnitVectorMap = RMap(:,3);

coordSysAxesMap = [X0Map                   Y0Map                   Z0Map                     % origin
                   X0Map+xUnitVectorMap(1) Y0Map+xUnitVectorMap(2) Z0Map+xUnitVectorMap(3)   % endpoint x unit vector
                   X0Map+yUnitVectorMap(1) Y0Map+yUnitVectorMap(2) Z0Map+yUnitVectorMap(3)   % endpoint y unit vector
                   X0Map+zUnitVectorMap(1) Y0Map+zUnitVectorMap(2) Z0Map+zUnitVectorMap(3)]; % endpoint z unit vector

[lat, lon, hEll] = minvtran(p.mstruct, coordSysAxesMap(:,1), coordSysAxesMap(:,2), coordSysAxesMap(:,3)); % lat, lon in degrees!
[coordSysAxesEcef(:,1), coordSysAxesEcef(:,2), coordSysAxesEcef(:,3)] = geodetic2ecef(lat*pi/180, lon*pi/180, hEll, p.mstruct.geoid); % lat, lon in radian!

xUnitVectorECEF = [coordSysAxesEcef(2,:) - coordSysAxesEcef(1,:)]';
yUnitVectorECEF = [coordSysAxesEcef(3,:) - coordSysAxesEcef(1,:)]';
zUnitVectorECEF = [coordSysAxesEcef(4,:) - coordSysAxesEcef(1,:)]';

RECEF = [xUnitVectorECEF yUnitVectorECEF zUnitVectorECEF];

[ome, phi, kap] = R2opk(RECEF);

obj.ome = ome;
obj.phi = phi;
obj.kap = kap;

objNew = obj;

% End --------------------------------------------------------------------------

% msg('E', procHierarchy);

end

author: pglira


text:
% Test file for chebtech/flipud.m

function pass = test_flipud(pref)

if ( nargin < 1 )
    pref = chebtech.techPref();
end

for n = 1:2
    if ( n == 1 )
        testclass = chebtech1();
    else
        testclass = chebtech2();
    end

    % Try some standard calls to flipud
    f = testclass.make(@(x) sin(x+.5), [], pref);
    g = testclass.make(@(x) sin(-x+.5), [], pref);
    h = flipud(f);
    pass(n, 1) = norm(g.coeffs - h.coeffs, inf) < 10*vscale(h)*eps;

    f = testclass.make(@(x) [sin(x+.5), exp(x)], [], pref);
    g = testclass.make(@(x) [sin(-x+.5), exp(-x)], [], pref);
    h = flipud(f);
    pass(n, 2) = norm(g.coeffs - h.coeffs, inf) < 10*max(vscale(h)*eps);

    f = testclass.make(@(x) sin(1i*x+.5), [], pref);
    g = testclass.make(@(x) sin(-1i*x+.5), [], pref);
    h = flipud(f);
    pass(n, 3) = norm(g.coeffs - h.coeffs, inf) < 10*vscale(h)*eps;

    f = testclass.make(@(x) [sin(x+.5), exp(1i*x)], [], pref);
    g = testclass.make(@(x) [sin(-x+.5), exp(-1i*x)], [], pref);
    h = flipud(f);
    pass(n, 4) = norm(g.coeffs - h.coeffs, inf) < 10*max(vscale(h)*eps);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function lambda = tri_upper_eigenvalues ( n )

%*****************************************************************************80
%
%% TRI_UPPER_EIGENVALUES returns the eigenvalues of the TRI_UPPER matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );

  lambda(1:n,1) = 1.0;

  return
end


text:
function [center, U, obj_fcn] = my_fcm(data, cluster_n, options)
%FCM Data set clustering using fuzzy c-means clustering.
%
%   [CENTER, U, OBJ_FCN] = FCM(DATA, N_CLUSTER) finds N_CLUSTER number of
%   clusters in the data set DATA. DATA is size M-by-N, where M is the number of
%   data points and N is the number of coordinates for each data point. The
%   coordinates for each cluster center are returned in the rows of the matrix
%   CENTER. The membership function matrix U contains the grade of membership of
%   each DATA point in each cluster. The values 0 and 1 indicate no membership
%   and full membership respectively. Grades between 0 and 1 indicate that the
%   data point has partial membership in a cluster. At each iteration, an
%   objective function is minimized to find the best location for the clusters
%   and its values are returned in OBJ_FCN.
%
%   [CENTER, ...] = FCM(DATA,N_CLUSTER,OPTIONS) specifies a vector of options
%   for the clustering process:
%       OPTIONS(1): exponent for the matrix U             (default: 2.0)
%       OPTIONS(2): maximum number of iterations          (default: 100)
%       OPTIONS(3): minimum amount of improvement         (default: 1e-5)
%       OPTIONS(4): info display during iteration         (default: 1)
%   The clustering process stops when the maximum number of iterations
%   is reached, or when the objective function improvement between two
%   consecutive iterations is less than the minimum amount of improvement
%   specified. Use NaN to select the default value.
%
%   Example
%       data = rand(100,2);
%       [center,U,obj_fcn] = fcm(data,2);
%       plot(data(:,1), data(:,2),'o');
%       hold on;
%       maxU = max(U);
%       % Find the data points with highest grade of membership in cluster 1
%       index1 = find(U(1,:) == maxU);
%       % Find the data points with highest grade of membership in cluster 2
%       index2 = find(U(2,:) == maxU);
%       line(data(index1,1),data(index1,2),'marker','*','color','g');
%       line(data(index2,1),data(index2,2),'marker','*','color','r');
%       % Plot the cluster centers
%       plot([center([1 2],1)],[center([1 2],2)],'*','color','k')
%       hold off;
%
%   See also FCMDEMO, INITFCM, IRISFCM, DISTFCM, STEPFCM.

%   Roger Jang, 12-13-94, N. Hickey 04-16-01
%   Copyright 1994-2002 The MathWorks, Inc.
%   $Revision: 1.13 $  $Date: 2002/04/14 22:20:38 $

if nargin ~= 2 & nargin ~= 3,
	error('Too many or too few input arguments!');
end

data_n = size(data, 1);
in_n = size(data, 2);

% Change the following to set default options
default_options = [2;	% exponent for the partition matrix U
		100;	% max. number of iteration
		1e-5;	% min. amount of improvement
		1];	% info display during iteration

if nargin == 2,
	options = default_options;
else
	% If "options" is not fully specified, pad it with default values.
	if length(options) < 4,
		tmp = default_options;
		tmp(1:length(options)) = options;
		options = tmp;
	end
	% If some entries of "options" are nan's, replace them with defaults.
	nan_index = find(isnan(options)==1);
	options(nan_index) = default_options(nan_index);
	if options(1) <= 1,
		error('The exponent should be greater than 1!');
	end
end

expo = options(1);		% Exponent for U
max_iter = options(2);		% Max. iteration
min_impro = options(3);		% Min. improvement
display = options(4);		% Display info or not

obj_fcn = zeros(max_iter, 1);	% Array for objective function

U = my_initfcm(cluster_n, data_n);	% changed by Romano Cicchetti in order to obtain reproducable results
% Main loop
for i = 1:max_iter,
	[U, center, obj_fcn(i)] = stepfcm(data, U, cluster_n, expo);
	if display,
		fprintf('Iteration count = %d, obj. fcn = %f\n', i, obj_fcn(i));
	end
	% check termination condition
	if i > 1,
		if abs(obj_fcn(i) - obj_fcn(i-1)) < min_impro, break; end,
	end
end

iter_n = i;	% Actual number of iterations
obj_fcn(iter_n+1:max_iter) = [];

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function [ handle ] = ml_plot_cv_grid_states(stats,options)
%ML_PLOT_CV_GRID_STATES Plots the results of grid search K-fold Cross
% Validation
%
%   input -----------------------------------------------------------------
%
%       o stats     : struct,   multi-layer structure, see output of
%                               ml_get_cv_grid_states.m
%
%       o options   : struct,   plot options.
%
%               options.title = 'title_figure';
%               options.param_names = ['C', 'sigma'];
%
%   output ----------------------------------------------------------------
%
%       o handle : figure handle
%

if ~isfield(options,'title'),options.title = 'add a title in options.title'; end
if ~isfield(options,'log_grid'),options.log_grid = 0'; end
if ~isfield(options,'svm_metrics'),options.svm_metrics = 0'; end

title_name = options.title;

[P, N] = size(stats.train.acc.mean');

if P > 1 && N > 1
        handle = figure('Color', [1 1 1],  'Position', [0, 1000, 1295, 455]);
        hold on;
        colormap hot;
        x = options.param_ranges(2,:);
        y = options.param_ranges(1,:);

        if (options.svm_metrics == 1)
            subplot(2, 3, 1)
        else
            subplot(1, 2, 1)
        end
        z = stats.train.acc.mean;
        contourf(x,y,z)
        if (options.log_grid ==1)
            set(gca,'xscale','log')
            set(gca,'yscale','log')
            x_range = options.param_ranges(2,:);
            y_range = options.param_ranges(1,:);
            set(gca, 'XTick', options.param_ranges(2,:))
            set(gca,'XTickLabel', cellstr(num2str(x_range(:), '%4.2f')))
            set(gca, 'YTick', options.param_ranges(1,:))
            set(gca,'YTickLabel', cellstr(num2str(y_range(:), '%4.2f')))
        end
        title('Train Accuracy (Mean)','FontSize',14, 'FontWeight','Normal')
        xlabel(options.param_names(2), 'FontSize',14, 'FontWeight','Normal')
        ylabel(options.param_names(1), 'FontSize',14, 'FontWeight','Normal')

        colorbar
        grid off
        axis square

        if (options.svm_metrics == 1)
            subplot(2, 3, 2)
        else
            subplot(1, 2, 2)
        end
        z = stats.train.fmeasure.mean;
        contourf(x,y,z)

        if (options.log_grid ==1)
            set(gca,'xscale','log')
            set(gca,'yscale','log')
            x_range = options.param_ranges(2,:);
            y_range = options.param_ranges(1,:);
            set(gca, 'XTick', options.param_ranges(2,:))
            set(gca,'XTickLabel', cellstr(num2str(x_range(:), '%4.2f')))
            set(gca, 'YTick', options.param_ranges(1,:))
            set(gca,'YTickLabel', cellstr(num2str(y_range(:), '%4.2f')))
        end

        title('Train F-measure (Mean)','FontSize',14, 'FontWeight','Normal')
        xlabel(options.param_names(2), 'FontSize',14, 'FontWeight','Normal')
        ylabel(options.param_names(1), 'FontSize',14, 'FontWeight','Normal')

        colorbar
        grid off
        axis square

        if (options.svm_metrics == 1)
            subplot(2, 3, 3)
            z = stats.train.fpr.mean;
            contourf(x,y,z)
            if (options.log_grid ==1)
                set(gca,'xscale','log')
                set(gca,'yscale','log')
                x_range = options.param_ranges(2,:);
                y_range = options.param_ranges(1,:);
                set(gca, 'XTick', options.param_ranges(2,:))
                set(gca,'XTickLabel', cellstr(num2str(x_range(:), '%4.2f')))
                set(gca, 'YTick', options.param_ranges(1,:))
                set(gca,'YTickLabel', cellstr(num2str(y_range(:), '%4.2f')))
            end
            title('Mean Train FPR (Fall-out)', 'FontSize',14, 'FontWeight','Normal')
            xlabel(options.param_names(2), 'FontSize',14, 'FontWeight','Normal')
            ylabel(options.param_names(1), 'FontSize',14, 'FontWeight','Normal')

            colorbar
            grid off
            axis square

            subplot(2, 3, 4)
            z = stats.train.tnr.mean;
            contourf(x,y,z)
            if (options.log_grid ==1)
                set(gca,'xscale','log')
                set(gca,'yscale','log')
                x_range = options.param_ranges(2,:);
                y_range = options.param_ranges(1,:);
                set(gca, 'XTick', options.param_ranges(2,:))
                set(gca,'XTickLabel', cellstr(num2str(x_range(:), '%4.2f')))
                set(gca, 'YTick', options.param_ranges(1,:))
                set(gca,'YTickLabel', cellstr(num2str(y_range(:), '%4.2f')))
            end
            title('Mean Train TNR (Specificity)', 'FontSize',14, 'FontWeight','Normal')
            xlabel(options.param_names(2), 'FontSize',14, 'FontWeight','Normal')
            ylabel(options.param_names(1), 'FontSize',14, 'FontWeight','Normal')

            colorbar
            grid off
            axis square

            subplot(2, 3, 5)
            z = stats.model.ratioSV.mean;
            contourf(x,y,z)
            if (options.log_grid ==1)
                set(gca,'xscale','log')
                set(gca,'yscale','log')
                x_range = options.param_ranges(2,:);
                y_range = options.param_ranges(1,:);
                set(gca, 'XTick', options.param_ranges(2,:))
                set(gca,'XTickLabel', cellstr(num2str(x_range(:), '%4.2f')))
                set(gca, 'YTick', options.param_ranges(1,:))
                set(gca,'YTickLabel', cellstr(num2str(y_range(:), '%4.2f')))
            end
            title('% of SV/M Datapoints', 'FontSize',14, 'FontWeight','Normal')
            xlabel(options.param_names(2), 'FontSize',14, 'FontWeight','Normal')
            ylabel(options.param_names(1), 'FontSize',14, 'FontWeight','Normal')

            colorbar
            grid off
            axis square

            subplot(2, 3, 6)
            z = stats.model.boundSV.mean;
            contourf(x,y,z)
            if (options.log_grid ==1)
                set(gca,'xscale','log')
                set(gca,'yscale','log')
            end
            title('% of Bounded SVs / Total SVs')
            xlabel(options.param_names(2), 'FontSize',14, 'FontWeight','Normal')
            ylabel(options.param_names(1), 'FontSize',14, 'FontWeight','Normal')

            colorbar
            grid off
            axis square
        end

        if isfield(stats,'test')

            handle2 =  figure('Color', [1 1 1],  'Position', [0, 125,  1295, 455]);
            hold on;
            colormap hot;
            x = options.param_ranges(2,:);
            y = options.param_ranges(1,:);

            subplot(2, 2, 1)
            z = stats.test.acc.mean;
            contourf(x,y,z)
            if (options.log_grid ==1)
                set(gca,'xscale','log')
                set(gca,'yscale','log')
                x_range = options.param_ranges(2,:);
                y_range = options.param_ranges(1,:);
                set(gca, 'XTick', options.param_ranges(2,:))
                set(gca,'XTickLabel', cellstr(num2str(x_range(:), '%4.2f')))
                set(gca, 'YTick', options.param_ranges(1,:))
                set(gca,'YTickLabel', cellstr(num2str(y_range(:), '%4.2f')))
            end
            title('Test Accuracy (Mean)' , 'FontSize',14, 'FontWeight','Normal')
            xlabel(options.param_names(2), 'FontSize',14, 'FontWeight','Normal')
            ylabel(options.param_names(1), 'FontSize',14, 'FontWeight','Normal')

            colorbar
            grid off
            axis square

            subplot(2, 2, 2)
            z = stats.test.fmeasure.mean;
            contourf(x,y,z)
            if (options.log_grid ==1)
                set(gca,'xscale','log')
                set(gca,'yscale','log')
                x_range = options.param_ranges(2,:);
                y_range = options.param_ranges(1,:);
                set(gca, 'XTick', options.param_ranges(2,:))
                set(gca,'XTickLabel', cellstr(num2str(x_range(:), '%4.2f')))
                set(gca, 'YTick', options.param_ranges(1,:))
                set(gca,'YTickLabel', cellstr(num2str(y_range(:), '%4.2f')))
            end

            title('Test F-measure (Mean)', 'FontSize',14, 'FontWeight','Normal')
            xlabel(options.param_names(2), 'FontSize',14, 'FontWeight','Normal')
            ylabel(options.param_names(1), 'FontSize',14, 'FontWeight','Normal')

            colorbar
            grid off
            axis square

                subplot(2, 2, 3)
                z = stats.test.fpr.mean;
                contourf(x,y,z)
                if (options.log_grid ==1)
                    set(gca,'xscale','log')
                    set(gca,'yscale','log')
                    x_range = options.param_ranges(2,:);
                    y_range = options.param_ranges(1,:);
                    set(gca, 'XTick', options.param_ranges(2,:))
                    set(gca,'XTickLabel', cellstr(num2str(x_range(:), '%4.2f')))
                    set(gca, 'YTick', options.param_ranges(1,:))
                    set(gca,'YTickLabel', cellstr(num2str(y_range(:), '%4.2f')))
                end
                title('Test FPR (Fall-out)', 'FontSize',14, 'FontWeight','Normal')
                xlabel(options.param_names(2), 'FontSize',14, 'FontWeight','Normal')
                ylabel(options.param_names(1), 'FontSize',14, 'FontWeight','Normal')

                colorbar
                grid off
               axis square

                subplot(2, 2, 4)
                z = stats.test.tnr.mean;
                contourf(x,y,z)
                if (options.log_grid ==1)
                    set(gca,'xscale','log')
                    set(gca,'yscale','log')
                    x_range = options.param_ranges(2,:);
                    y_range = options.param_ranges(1,:);
                    set(gca, 'XTick', options.param_ranges(2,:))
                    set(gca,'XTickLabel', cellstr(num2str(x_range(:), '%4.2f')))
                    set(gca, 'YTick', options.param_ranges(1,:))
                    set(gca,'YTickLabel', cellstr(num2str(y_range(:), '%4.2f')))
                end
                title('Test TNR (Specificity)', 'FontSize',14, 'FontWeight','Normal')
                xlabel(options.param_names(2), 'FontSize',14, 'FontWeight','Normal')
                ylabel(options.param_names(1), 'FontSize',14, 'FontWeight','Normal')

                colorbar
                grid off
                axis square
        end

else
        handle = figure('Color', [1 1 1]);

        x_index = options.param_ranges;
        subplot(1,2,1)
        errorbar(x_index,stats.test.acc.mean,stats.test.acc.std,'-rs'); hold on;
        errorbar(x_index,stats.train.acc.mean,stats.train.acc.std,'-gs');

        xlabel(options.param_names, 'FontSize',14, 'FontWeight','Normal');
        ylabel('Accuracy','FontSize',14);
        xlim([x_index(1) x_index(end)])

        if (options.log_grid ==1)
                set(gca,'xscale','log')
                x_range = options.param_ranges;
                set(gca, 'XTick', options.param_ranges)
                set(gca,'XTickLabel', cellstr(num2str(x_range(:), '%4.2f')))
        end

        ylim([min(stats.test.acc.mean)*3/4 max(stats.train.acc.mean)*1.1])
        legend({'Test','Train'},'FontSize',14);
        title(title_name,'FontSize',14, 'FontWeight','Normal');
        axis square
        box on;
        grid on;

        subplot(1,2,2)

        errorbar(x_index,stats.test.fmeasure.mean,stats.test.fmeasure.std,'-rs'); hold on;
        errorbar(x_index,stats.train.fmeasure.mean,stats.train.fmeasure.std,'-gs');

        xlabel(options.param_names, 'FontSize',14, 'FontWeight','Normal');
        ylabel('F-Measure','FontSize',14);
        xlim([x_index(1) x_index(end)])
        if (options.log_grid ==1)
            set(gca,'xscale','log')
            x_range = options.param_ranges;
            set(gca, 'XTick', options.param_ranges)
            set(gca,'XTickLabel', cellstr(num2str(x_range(:), '%4.2f')))
        end

        ylim([min(stats.test.fmeasure.mean)*3/4 max(stats.train.fmeasure.mean)*1.1])
        legend({'Test','Train'},'FontSize',14);
        title(title_name,'FontSize',14, 'FontWeight','Normal');
        axis square
        box on;
        grid on;

end

end

author: epfl-lasa


text:
function G = divgrad(F)
%DIVGRAD   Laplacian of a CHEBFUN2V.
%   F = DIVGRAD(F) returns the Laplacian of a CHEBFUN2V i.e.,
%       divgrad(F) = F(1)_xx + F(2)_yy
%
% This command is not defined for a chebfun2v with 3 components.
%
% Also see CHEBFUN2V/LAP.

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

nComponents = F.nComponents;
if ( nComponents > 2 )
    error('CHEBFUN:CHEBFUN2V:divgrad:components',...
        'Command is not defined for CHEBFUN2V objects with >2 components.')
end

Fc = F.components;
G = diff(Fc{1}, [2,0]) + diff(Fc{2}, [0,2]);

end


text:
% Small technical example.
% Shows how to call for the multiresolution decomposition of an image
% and shows the detailed bookkeeping.
%
disp('Small technical example.');
disp('Shows how to call for the multiresolution decomposition of an image');
disp('and shows the detailed bookkeeping.');
disp('FOR MORE INFORMATION:  help QLiftDec2');
disp(' ');
disp('See also the report http://repository.cwi.nl:8888/cwi_repository/docs/IV/04/04178D.pdf');
disp('Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>');
disp(' (C) 1998-2006 Stichting CWI, Amsterdam, The Netherlands');
disp(' ');
%---PARAMETERS-----------------------------------------------------------------
% How to execute, set parameters
N = 6;                   %  maximum level (even number) in lifting scheme
filtername = 'Neville4';
%
%---INSERT YOUR IMAGE HERE-----------------------------------------------------
if exist('imread','file') == 2
  Orig = double(imread('zenithgray.TIF','tiff'));
else
  load zenithgray; Orig = zenithgray; clear zenithgray;
end
%
disp([' Dimensions of original      ' int2str( size(Orig) )]);
%---DECOMPOSITION--------------------------------------------------------------
disp([' Filter type is ' filtername]);
[C,S] = QLiftDec2(Orig,N,filtername);
disp([' Dimensions of decomposition ' int2str( size(C) )]);
disp(' Contents of bookkeeping ');
disp(int2str(S));
%


text:
classdef net_model_acp < mp.net_model_ac & mp.form_acp

%   MATPOWER
%   Copyright (c) 2019-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

    properties
        va = [];
        vm = [];
    end

    methods
        function obj = net_model_acp()
            obj@mp.net_model_ac();
            obj.element_classes = ...
                { @mp.nme_bus_acp, @mp.nme_gen_acp, @mp.nme_load_acp, ...
                    @mp.nme_branch_acp, @mp.nme_shunt_acp };

            %% Due to a bug related to inheritance in constructors in
            %% Octave 5.2 and earlier (https://savannah.gnu.org/bugs/?52614),
            %% INIT_SET_TYPES() cannot be called directly in the
            %% MP_IDX_MANAGER constructor, as desired.
            %%
            %% WORKAROUND:  INIT_SET_TYPES() is called explicitly as needed
            %%              (if obj.node is empty) in BUILD() and DISPLAY(),
            %%              after object construction, but before object use.
        end

        function obj = def_set_types(obj)
            def_set_types@mp.net_model_ac(obj);     %% call parent first
            obj.set_types.va = 'VOLTAGE ANG VARS (va)';
            obj.set_types.vm = 'VOLTAGE MAG VARS (vm)';
        end

        function va = initial_voltage_angle(obj, idx)
            va = obj.params_var('va');
            if nargin > 1 && ~isempty(idx)
                va = va(idx);
            end
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function offset = find_offset_multi_source(samples1, samples2)
    [c, lags] = xcorr(samples1, samples2);
    [~, ind] = max(c);
    offset = lags(ind) + 1; % add 1 for Matlab indexing
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function stroud_test19 ( )

%*****************************************************************************80
%
%% TEST19 tests CUBE_UNIT_3D, QMULT_3D, RECTANGLE_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  global FUNC_3D_INDEX;

  a1 = -1.0;
  b1 = +1.0;

  a(1) = -1.0;
  a(2) = -1.0;
  a(3) = -1.0;
  b(1) = 1.0;
  b(2) = 1.0;
  b(3) = 1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST19\n' );
  fprintf ( 1, '  CUBE_UNIT_3D approximates integrals\n' );
  fprintf ( 1, '    in the unit cube in 3D.\n' );
  fprintf ( 1, '  QMULT_3D approximates triple integrals.\n' );
  fprintf ( 1, '  RECTANGLE_3D approximates integrals\n' );
  fprintf ( 1, '    in a rectangular block.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '    F(X)    CUBE_UNIT_3D  QMULT_3D      RECTANGLE_3D\n' );
  fprintf ( 1, '\n' );

  num = function_3d_num ( );

  for i = 1 : num

    FUNC_3D_INDEX = i;

    result1 = cube_unit_3d ( 'function_3d' );
    result2 = qmult_3d ( 'function_3d', a1, b1, 'fu18', 'fl18', 'fu28', 'fl28' );
    result3 = rectangle_3d ( 'function_3d', a, b );

    fname = function_3d_name ( i );

    fprintf ( 1, '  %s  %12f  %12f  %12f\n', ...
      fname, result1, result2, result3 );

  end

  return
end

author: johannesgerer


text:
clear all; close all;
I=imread('circbw.tif');
se=strel('disk', 3);
J=imdilate(I, se);
a1=bwarea(I)
a2=bwarea(J)
(a2-a1)/a1
figure;
subplot(121);  imshow(I);
subplot(122);  imshow(J);

author: HuangCongQing


text:
function f=midi2frq(n,s)
%MIDI2FRQ	Convert musical note numbers to frequencies F=(N,S)
%		s is:	'e' equal tempered (default)
%			'p' pythagorean scale
%			'j' just intonation
%
% notes are numbered in semitones with middle C being 60
% On the equal tempered scale, note 69 (the A above middle C)
% has a frequency of 440 Hz.
%
% see FRQ2NOTE for the inverse transform

% Pythagorean
%     sharps 1 2187/2048 9/8 19683/16384 81/64 4/3 729/512  3/2 6561/4096 27/16 59049/32768 243/128 2
%     flats  1 256/243   9/8 32/27       81/64 4/3 1024/729 3/2 128/81    27/16 16/9        243/128 2
%
% Just Intonation
%     sharps 1 25/24 9/8 75/64 5/4 4/3 45/32  3/2 25/16 5/3 225/128 15/8 2
%     flats  1 16/15 9/8 6/5   5/4 4/3 108/75 3/2 8/5   5/3 18/10   15/8 2

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: midi2frq.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin > 1
  if s(1)=='p'
    r=[256/243 9/8 32/27 81/64 4/3 729/512 3/2 128/81 27/16 16/9 243/128];
  elseif s(1)=='j'
    r=[16/15 9/8 6/5 5/4 4/3 36/25 3/2 8/5 5/3 9/5 15/8];
  else
    r=0;
  end
  if r(1)
    c=[0 0 12*log(r)/log(2)-(1:11) 0];
    nm=mod(n,12);
    na=floor(nm);
    nb=nm-na;
    f=440*exp((n+c(na+2).*(1-nb)+c(na+3).*nb-69)*log(2)/12);
  else
    f=440*exp((n-69)*log(2)/12);
  end
else
  f=440*exp((n-69)*log(2)/12);
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function [patchHandles textHandles] = hlp_vrect(x, varargin)
% Tim Mullen, 2011, SCCN/INC/UCSD

g = finputcheck(varargin, ...
    {'yscale'           'real'      []           1; ...         % vertical scaling (relative to ylimits)
    'dock'              'string'    {'bottom','top','none'} 'top'; ...    % whether to 'dock' the patch at the top or bottom of the axis. If 'none' then scale is location on yaxis of [bottom top] of rect
    'label'             ''          []              {}; ...
    'axesHandle'        'real'      []              gca; ...
    'textPosition'      'real'      [0 1]           [0.5 0.02]; ...
    'textProperties'    'cell'      []              {}; ...
    'patchProperties'   'cell'      []              {'FaceColor',[0.7 0.7 1],'FaceAlpha',0.5,'EdgeColor',[0.2 0.2 0.2],'EdgeAlpha',0.5} ...
    },'hlp_vrect','error','quiet');

if ischar(g)
    help hlp_vrect;
    error(g);
end

pp = hlp_varargin2struct(g.patchProperties);
if ~isfield(pp,'FaceColor');
    FaceColor = [0.7 0.7 1];
else
    FaceColor = pp.FaceColor;
end
clear pp;

xsize = size(x);
if( xsize(2)~=2 )
    error('x must have two columns');
end

if ~iscell(g.label)
    g.label = {g.label};
end

if length(g.label)~=xsize(1)
    g.label = repmat(g.label,xsize(1),1);
end

if size(g.textPosition,2) ~= 2
    error('textPosition must have two columns [ (x,y) position ]');
end

if size(g.textPosition,1)~=xsize(1)
    g.textPosition=repmat(g.textPosition,xsize(1),1);
end

textPositionX = g.textPosition(:,1);
textPositionY = g.textPosition(:,2);

holdState = ishold(g.axesHandle);
hold(g.axesHandle, 'on');

if strcmpi(g.dock,'none')
    yLimits = g.yscale;
else
    yl = get(g.axesHandle,'ylim');             % Row vector
    yLimits = yl.*g.yscale;
end

if any(g.yscale<1)
    switch g.dock
        case 'bottom'
            yLimits = yLimits-(yLimits(1)-yl(1))+0.001*diff(yLimits);
        case 'top'
            yLimits = yLimits+(yl(2)-yLimits(2))-0.001*diff(yLimits);
    end
end

% bring desired axes into focus
set(gcf,'CurrentAxes',g.axesHandle);

patchHandles = zeros(1,xsize(1));
textHandles = [];

for k=1:xsize(1)
    % for each patch
    patchHandles(k) = patch([x(k,1) x(k,1) x(k,2) x(k,2)], [yLimits yLimits(end:-1:1)],FaceColor,g.patchProperties{:});

    if( ~isempty(g.label) )
%         xLimits = get(g.axesHandle,'xlim');
%         xLowerLimit = xLimits(1);
%         xUpperLimit = xLimits(2);
%         xRange      = xUpperLimit - xLowerLimit;
        xPosition   = x(k,1) + textPositionX(k,:)*(x(k,2)-x(k,1));

        yUpperLimit = yLimits(2);
        yLowerLimit = yLimits(1);
        yRange      = yUpperLimit - yLowerLimit;
        Yposition   = yLowerLimit + textPositionY(k,:)*yRange;
%         Yposition   = repmat(yPosition, length(x), 1);

        % textHandle is a vector correspond to the line
        textHandles(k) = text(xPosition, Yposition, g.label{k}, 'Parent', g.axesHandle,g.textProperties{:});
        % center text
        extent = get(textHandles(k),'Extent');
        set(textHandles(k),'Position',[xPosition-extent(3)/2 Yposition]);

    end
end

if( holdState==false )
    hold(g.axesHandle, 'off');
end

repo: REST


text:
function varargout = process_canoltymap2( varargin )
% This function generates Canolty like maps (Science 2006, figure 1) for the input signal.

% @=============================================================================
% 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: Esther Florin, Sylvain Baillet, 2011-2013
%          Francois Tadel, 2013-2014

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Canolty maps (FileB=MaxPAC)';
    sProcess.Category    = 'File2';
    sProcess.SubGroup    = 'Frequency';
    sProcess.Index       = 661;
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'raw',      'data',     'results',  'matrix'};
    sProcess.OutputTypes = {'timefreq', 'timefreq', 'timefreq', 'timefreq'};
    sProcess.nInputs     = 2;
    sProcess.nMinFiles   = 1;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Resting#Canolty_maps';
    % ==== INPUT ====
    sProcess.options.label_in.Comment = 'Input options:';
    sProcess.options.label_in.Type    = 'label';
    % === TIME WINDOW
    sProcess.options.timewindow.Comment = 'Time window: ';
    sProcess.options.timewindow.Type    = 'timewindow';
    sProcess.options.timewindow.Value   = [];
    % === SENSOR SELECTION
    sProcess.options.target_data.Comment    = 'Sensor types or names (empty=all): ';
    sProcess.options.target_data.Type       = 'text';
    sProcess.options.target_data.Value      = 'MEG, EEG';
    sProcess.options.target_data.InputTypes = {'data', 'raw'};
    % === SCOUTS SELECTION
    sProcess.options.scouts.Comment    = 'Use scouts';
    sProcess.options.scouts.Type       = 'scout_confirm';
    sProcess.options.scouts.Value      = {};
    sProcess.options.scouts.InputTypes = {'results'};
    % === 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'};
    % === SCOUT TIME ===
    sProcess.options.scouttime.Comment    = {'Before', 'After', 'When to apply the scout function:'};
    sProcess.options.scouttime.Type       = 'radio_line';
    sProcess.options.scouttime.Value      = 1;
    sProcess.options.scouttime.InputTypes = {'results'};
    % === ROW NAMES
    sProcess.options.target_tf.Comment    = 'Row names or indices (empty=all): ';
    sProcess.options.target_tf.Type       = 'text';
    sProcess.options.target_tf.Value      = '';
    sProcess.options.target_tf.InputTypes = {'timefreq', 'matrix'};

    % ==== ESTIMATOR ====
    sProcess.options.label_method.Comment = '
Estimator options:';
    sProcess.options.label_method.Type    = 'label';
    % === EPOCH TIME
    sProcess.options.epochtime.Comment = 'Epoch time: ';
    sProcess.options.epochtime.Type    = 'range';
    sProcess.options.epochtime.Value   = {[-0.5, 0.5], 'ms', []};
    % === MAX_BLOCK_SIZE
    sProcess.options.max_block_size.Comment = 'Number of signals to process at once: ';
    sProcess.options.max_block_size.Type    = 'value';
    sProcess.options.max_block_size.Value   = {100, ' ', 0};

    % ==== OUTPUT ====
    sProcess.options.label_out.Comment = '
Output configuration:';
    sProcess.options.label_out.Type    = 'label';
    % === SAVE AVERAGED LOW-FREQ SIGNALS
    sProcess.options.save_erp.Comment = 'Save averaged low frequency signals';
    sProcess.options.save_erp.Type    = 'checkbox';
    sProcess.options.save_erp.Value   = 1;
end

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

%% ===== RUN =====
function OutputFile = Run(sProcess, sInputA, sInputB) %#ok<DEFNU>
    % Load the optimal low-frequency values from the MaxPAC file in InputB
    DataMat = load(file_fullpath(sInputB.FileName), 'sPAC');
    if ~isfield(DataMat, 'sPAC') || isempty(DataMat.sPAC) || ~isfield(DataMat.sPAC, 'NestingFreq') || isempty(DataMat.sPAC.NestingFreq)
        bst_report('Error', sProcess, sInputB, 'Invalid MaxPAC file in FilesB.');
        OutputFile = {};
        return;
    end
    % Set option from the file
    sProcess.options.lowfreq.Value{1} = DataMat.sPAC.NestingFreq;
    % Run the megPAC process
    OutputFile = process_canoltymap('Run', sProcess, sInputA);
end


text:
function [gew,pve,H] = spm_csd2gew(csd,Hz,u)
% Convert cross sspectral density to Geweke Granger causality
% FORMAT [gew,pve,H] = spm_csd2gew(csd,Hz)
%
% ccf  (N,m,m)   - cross covariance functions
% Hz   (n x 1)   - vector of frequencies (Hz)
% u    (1)       - regularizer (default: 1);
%
% gwe  (N,m,m)   - Geweke's frequency domain Granger causality
% pve  (N,m,m)   - proportion of variance explained
% H    (N,m,m)   - transfer function matrix
%
% This routine uses the Wilson-Burg algorithm to perform spectral matrix
% factorisation. The minimum phase factor is then used to form the noise
% covariance (covariance of the innovations) and implicitly derive the
% transfer functions (and spectral Granger causality).
%
% See also:
%  spm_ccf2csd.m, spm_ccf2mar, spm_csd2ccf.m, spm_csd2mar.m, spm_mar2csd.m,
%  spm_csd2coh.m, spm_dcm_mtf.m, spm_Q.m, spm_mar.m and spm_mar_spectral.m
%
%__________________________________________________________________________
% Copyright (C) 2014 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_csd2gew.m 5908 2014-03-05 20:31:57Z karl $

% preliminaries
%--------------------------------------------------------------------------
if nargin < 3
    try
        [gew,pve,H] = spm_csd2gew(csd,Hz,2);
        return
    catch
        [gew,pve,H] = spm_csd2gew(csd,Hz,16);
        return
    end
end

% pad spectrum if necessary
%--------------------------------------------------------------------------
iw    = 1 + round(Hz/(Hz(2) - Hz(1)));
n     = size(csd,2);
nw    = 257;
is    = ceil(nw/2):nw;

% Wilson-Burg algorithm
%==========================================================================

% initialise transfer function
%--------------------------------------------------------------------------
H     = zeros(nw,n,n);
P     = zeros(nw,n,n);
e     = norm(squeeze(max(csd,[],1)))/128;
E     = eye(n,n)*e;

P(iw,:,:) = csd;
for i = 1:n
    P(:,i,i) = P(:,i,i) + e;
    H(:,i,i) = sqrt(P(:,i,i));
end

% iterate until convergence: solve for H*H' = P
%--------------------------------------------------------------------------
for t = 1:128

    % compute left-hand side (deconvolution)
    %----------------------------------------------------------------------
    for w = 1:nw
        S        = squeeze(H(w,:,:)) + E;
        A(w,:,:) = (eye(n,n) + S*S'\squeeze(P(w,:,:)));
    end

    % retain causal signal and half zero lag
    %----------------------------------------------------------------------
    S           = ifft(A);
    S(is,:,:)   = 0;
    S(1,:,:)    = S(1,:,:)/2;
    A           = fft(S);

    % recover next update (convolution)
    %----------------------------------------------------------------------
    nrm   = zeros(nw,1);
    for w = 1:nw
        U        = squeeze(A(w,:,:));
        H(w,:,:) = squeeze(H(w,:,:))*U^(1/u);
        nrm(w)   = nrm(w) + norm(eye(n,n) - U,'inf');
    end

    % break if convergence
    %----------------------------------------------------------------------
    nrm = mean(nrm);
    if nrm < 1e-6, break, end
    if nrm > 8,   return, end

end

% transfer function and noise covariance
%==========================================================================

% get noise covariance
%--------------------------------------------------------------------------
S     = ifft(H);
R     = squeeze(S(1,:,:));
C     = real(R*R');
c     = sqrtm(C);

% recover transfer function
%--------------------------------------------------------------------------
for w = 1:nw
    H(w,:,:) = squeeze(H(w,:,:))/c;
    P(w,:,:) = squeeze(H(w,:,:))*C*squeeze(H(w,:,:))';
end

% Geweke Granger Causality in the Frequency domain
%--------------------------------------------------------------------------
pve   = zeros(nw,n,n);
gew   = zeros(nw,n,n);
for j = 1:n
    for k = 1:n
        rkj        = C(j,j) - (C(j,k)^2)/C(k,k);
        sk         = abs(P(:,k,k));
        hkj        = abs(H(:,k,j)).^2;
        pve(:,k,j) = rkj*hkj./sk;
        gew(:,k,j) = -log(1 - pve(:,k,j));
    end
end

% return  specified frequencies
%--------------------------------------------------------------------------
gew = gew(iw,:,:);
pve = pve(iw,:,:);
H   = H(iw,:,:);

author: spm


text:
function g = ngaussNoiseGradientParam(noise, mu, varsigma, y)

% NGAUSSNOISEGRADIENTPARAM Gradient of NGAUSS noise's parameters.
% FORMAT
% DESC computes the gradient of the log Z of the noiseless Gaussian noise model with respect to the of functions with respect to the
% noiseless Gaussian
% noise's parameters.
% ARG noise : the noise structure for which the gradients are being
% computed.
% ARG mu : the input means for which the gradients are being computed.
% ARG varSigma : the input variances for which the gradients are being computed.
% ARG y : the target values for the noise model.
% RETURN g : gradients of the log Z with respect to
% the noise parameters. The ordering of the vector should match
% that provided by the function noiseExtractParam.
%
%
% SEEALSO ngaussNoiseParamInit, ngaussnoiseGradVals, noiseGradientParam
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% NOISE

D = size(y, 2);
u = zeros(size(y));

for i = 1:D
  mu(:, i) = mu(:, i) + noise.bias(i);
end

u = y - mu;
nu = 1./(varsigma+noise.sigma2);
u = u.*nu;
gbias = sum(u, 1);
g = [gbias];

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%Extremely similar to the false position method. The main difference is
%that the secant method does not actually have a defined interval where
%the root lies on. It converges faster than the false position method,
%but it is not always guaranteed to converge.

%INPUTS:
%Function handle f
%x1 = a
%x2 = b
%maximum tolerated error

%OUTPUTS:
%An approximated value for the root of f.

%Written by MatteoRaso

function y = secant(f, a, b, error)
  x = [a, b];
  n = 2;
  while abs(f(x(n))) > error
    x(n + 1) = -f(x(n)) * (x(n) - x(n - 1)) / (f(x(n)) - f(x(n - 1))) + x(n);
    n = n + 1;
    disp(f(x(n)))
  endwhile
  A = ["The root is approximately ", num2str(x(n))];
  disp(A)
  y = x(n);
endfunction

author: TheAlgorithms

repo: MATLAB-Octave


text:
% Test file for ADCHEBFUN Bessel functions

function pass = test_bessel

% List of airy functions to test, with different number of inputs
besselj1 = @(f) besselj(1, f);
besselj25 = @(f) besselj(2.5, f);

funcList = {besselj1, besselj25};

% Call the ADCHEBFUN testUnary() method to do the tests.
pass = adchebfun.testUnary(funcList);

end


text:
% Make an HMM with autoregressive Gaussian observations (switching AR model)
%   X1 -> X2
%   |     |
%   v     v
%   Y1 -> Y2

seed = 0;
rand('state', seed);
randn('state', seed);

intra = zeros(2);
intra(1,2) = 1;
inter = zeros(2);
inter(1,1) = 1;
inter(2,2) = 1;
n = 2;

Q = 2; % num hidden states
O = 2; % size of observed vector

ns = [Q O];
dnodes = 1;
onodes = [2];
bnet = mk_dbn(intra, inter, ns, 'discrete', dnodes, 'observed', onodes);

bnet.CPD{1} = tabular_CPD(bnet, 1);
bnet.CPD{2} = gaussian_CPD(bnet, 2);
bnet.CPD{3} = tabular_CPD(bnet, 3);
bnet.CPD{4} = gaussian_CPD(bnet, 4);

T = 10; % fixed length sequences

engine = {};
%engine{end+1} = hmm_inf_engine(bnet);
engine{end+1} = jtree_unrolled_dbn_inf_engine(bnet, T);
%engine{end+1} = smoother_engine(hmm_2TBN_inf_engine(bnet));
%engine{end+1} = smoother_engine(jtree_2TBN_inf_engine(bnet));

inf_time = cmp_inference_dbn(bnet, engine, T, 'check_ll',1);
learning_time = cmp_learning_dbn(bnet, engine, T, 'check_ll', 1);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function juldaystr = datenum2julday(dnum)
[yyyy, mm, dd, hh, mi, ss] = datevec(dnum);
dnum_jan1 = datenum(yyyy, 1, 1);
dnum_diff = dnum - dnum_jan1;
julday = ceil(dnum_diff + eps); % eps turns datenum(1997,1,1) into day 1 rather than day 0
juldaystr = sprintf('%4d%03d', yyyy, julday);

author: geoscience-community-codes

repo: GISMO