text:
%   ****** See nozzle.m for instructions ******
%   solver.m

%   Solver for the macormack method
function [Q] = solver(Q)
global gamma R

[u,v,rho,p,e,T,ss,F,G] = flowvars(Q);

%   Take one MacCormack step
[Q] = mac(Q);

end

%%%%%%%%    flowvars    %%%%%%%%
function [u,v,rho,p,e,T,ss,F,G] = flowvars(Q)
global gamma R

%	Calculate the actual flow variables at each time step
rho = Q(:,:,1);
u = Q(:,:,2)./rho;
v = Q(:,:,3)./rho;
e = Q(:,:,4);
p = (gamma-1)*(e-(1/2)*rho.*(u.*u+v.*v));
T = p./(R*rho);
ss = sqrt(abs(gamma*R*T));

F(:,:,1) = rho.*u;
F(:,:,2) = rho.*u.*u + p;
F(:,:,3) = rho.*u.*v;
F(:,:,4) = (e+p).*u;

G(:,:,1) = rho.*v;
G(:,:,2) = rho.*u.*v;
G(:,:,3) = rho.*v.*v + p;
G(:,:,4) = (e+p).*v;

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%    mac    %%%%%%%%

function [Q] = mac(Q)
global gamma R

Q0 = Q;

%	Forward flux
[Qflux,dt] = flux_mc(Q,-1);
Qbar = Q - dt*Qflux;
Q = Qbar;

[Q] = boundary(Q);

%	Backward flux
[Qflux,dt] = flux_mc(Q,0);
Q = (1/2)*(Q0 + Qbar - dt*Qflux );

[Q] = boundary(Q);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%    flux_mc    %%%%%%%%

function [Qflux,dt] = flux_mc(Q,dd)
global x y Vol cfl

[u,v,rho,p,e,T,ss,F,G] = flowvars(Q);
nx = size(x,1);
ny = size(x,2);

a(1:nx+1,1:ny+1) = 0;
b(1:nx+1,1:ny+1) = 0;
c(1:nx+1,1:ny+1) = 0;
Qflux(1:nx+1,1:ny+1,1:4) = 0;

%	Get the fluxes
for i = 2: size(x,1)
    for j = 2: size(x,2)

	ii = i-1;
    jj = j-1;

    %	Right face
	sfpx = y(ii+1,jj+1)-y(ii+1,jj);
    sfpy = -( x(ii+1,jj+1)-x(ii+1,jj) );

	%	Left face
    sfmx = -( y(ii,jj+1) - y(ii,jj) );
    sfmy = ( x(ii,jj+1)-x(ii,jj) );

    %	Top face
    sgpx = -( y(ii+1,jj+1) - y(ii,jj+1) );
    sgpy = x(ii+1,jj+1) - x(ii,jj+1);

    %	Bottom face
    sgmx = ( y(ii+1,jj)-y(ii,jj) );
    sgmy = -( x(ii+1,jj) - x(ii,jj) );

	%	Get the flux
	Qflux(i,j,:) = ( F(i+1+dd,j,:)*sfpx + G(i+1+dd,j,:)*sfpy + ...
   		F(i+dd,j,:)*sfmx + G(i+dd,j,:)*sfmy + F(i,j+1+dd,:)*sgpx ...
        + G(i,j+1+dd,:)*sgpy + F(i,j+dd,:)*sgmx + G(i,j+dd,:)*sgmy );

	%	Normalize by Volume
    Qflux(i,j,:) = Qflux(i,j,:)./Vol(i,j);

	%	CFL terms
    a(i,j) = abs(u(i,j)*sfpx + v(i,j)*sfpy);
    b(i,j) = abs(u(i,j)*sgpx + v(i,j)*sgpy);
    c(i,j) = ss(i,j)*sqrt(abs( sfpx^2 + sfpy^2) ...
        + abs( sgpx^2 + sgpy^2) );

    end
end

dt = max(max((a+b+c)./Vol));
dt = cfl/dt;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%    boundary    %%%%%%%%

function [Q] = boundary(Q)
global x y gamma R P_amb T_c
[u,v,rho,p,e,T,ss,F,G] = flowvars(Q);

%   Problem boundary conditions here
nx = size(x,1);
ny = size(x,2);

%	Top Wall
p(:,ny+1) = p(:,ny);
v(:,ny+1) = 0;
u(:,ny+1) = 0;
rho(:,ny+1) = rho(:,ny);

%   Symmetry line
p(:,1) = p(:,2);
v(:,1) = -v(:,2);
u(:,1) = u(:,2);
rho(:,1) = rho(:,2);

%	Inflow-shouldn't change from initialization
u(1,:) = sqrt(gamma*R*T_c);
v(1,:) = v(2,:);

%	Out flow - set to upstream cells
u(nx+1,:) = u(nx,:);
v(nx+1,:) = v(nx,:);
p(nx+1,:) = p(nx,:);  %P_amb;
rho(nx+1,:) = rho(nx,:);

%	EOS
e = p/(gamma-1) + (1/2)*rho.*(u.*u + v.*v);

Q(:,:,1) = rho;
Q(:,:,2) = rho.*u;
Q(:,:,3) = rho.*v;
Q(:,:,4) = e;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: Sable


text:
function model = vargplvmAddDynamics(model, type, varargin)
% VARGPLVMADDDYNAMICS Add a dynamics structure to the model.
%
%	Description:
%
%	VARGPLVMADDDYNAMICS(MODEL, TYPE, ...) adds a dynamics structure to the
%	VARGPLVM.
%	 Arguments:
%	  MODEL - the model to add dynamics to.
%	  TYPE - the type of dynamics model to add in.
%	  ... - additional arguments to be passed on creation of the
%	   dynamics model.
%
% COPYRIGHT : Andreas C. Damianou, 2010-2011
% COPYRIGHT : Michalis K. Titsias, 2010-2011
% COPYRIGHT : Neil D. Lawrence, 2010-2011
%
%
%	See also
%	MODELCREATE
% VARGPLVM

type = [type 'Dynamics'];
model.dynamics = modelCreate(type, model.q, model.q, model.X, varargin{:}); % e.g. vargpTimeDynamicsCreate
params = vargplvmExtractParam(model);
model.dynamics.nParams = length(modelExtractParam(model.dynamics));
model.nParams = model.nParams + model.dynamics.nParams;
model = vargplvmExpandParam(model, params);

% This will be moved to "initDynamics"
%med = median(model.vardist.covars(:));
%if med < 0.3 || med > 0.7
%   warning('!!! Median value of variational covariances is %1.2f.\n', med);
%end

author: SheffieldML

repo: vargplvm


text:
function varargout = process_zscore_ab( varargin )
% PROCESS_ZSCORE: Compute Z-Score for a matrix A (normalization respect to a baseline).
%
% DESCRIPTION:  For each channel:
%     1) Compute mean m and variance v for baseline
%     2) For each time sample, subtract m and divide by v

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Z-score static (A=baseline) [DEPRECATED]';
    sProcess.FileTag     = 'zscore';
    sProcess.Category    = 'Filter2';
    sProcess.SubGroup    = 'Standardize';
    sProcess.Index       = 0;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/SourceEstimation#Z-score';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'results', 'timefreq', 'matrix'};
    sProcess.OutputTypes = {'data', 'results', 'timefreq', 'matrix'};
    sProcess.nInputs     = 2;
    sProcess.nMinFiles   = 1;
    % Default values for some options
    sProcess.isSourceAbsolute = 0;
    sProcess.processDim       = 1;    % Process channel by channel
    sProcess.isPaired         = 1;

    % Definition of the options
    sProcess.options.description.Comment = ['For each signal in input:
' ...
                                            '1) FilesA: Compute mean m and variance v for the baseline
' ...
                                            '2) FilesB: For each time sample, subtract m and divide by v
' ...
                                            'Z = (Data - m) / v

'];
    sProcess.options.description.Type    = 'label';
    % === Baseline time window
    sProcess.options.baseline.Comment = 'Baseline (Files A):';
    sProcess.options.baseline.Type    = 'baseline';
    sProcess.options.baseline.Value   = [];
    % === Sensor types
    sProcess.options.sensortypes.Comment = 'Sensor types or names (empty=all): ';
    sProcess.options.sensortypes.Type    = 'text';
    sProcess.options.sensortypes.Value   = 'MEG, EEG';
    sProcess.options.sensortypes.InputTypes = {'data'};
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    % Get frequency band
    if isfield(sProcess.options, 'baseline') && isfield(sProcess.options.baseline, 'Value') && iscell(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value{1})
        Time = sProcess.options.baseline.Value{1};
    else
        Time = [];
    end
    % Add frequency band
    if isempty(Time)
        Comment = 'Z-score normalization: [All file]';
    elseif any(abs(Time) > 2)
        Comment = sprintf('Z-score normalization (static): [%1.3fs,%1.3fs]', Time(1), Time(2));
    else
        Comment = sprintf('Z-score normalization (static): [%dms,%dms]', round(Time(1)*1000), round(Time(2)*1000));
    end
end

%% ===== RUN =====
function sInputB = Run(sProcess, sInputA, sInputB) %#ok<DEFNU>
    % Get options
    if isfield(sProcess.options, 'baseline') && isfield(sProcess.options.baseline, 'Value') && iscell(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value{1})
        BaselineBounds = sProcess.options.baseline.Value{1};
    else
        BaselineBounds = [];
    end
    % Get baseline indices
    if ~isempty(BaselineBounds)
        iBaseline = bst_closest(sProcess.options.baseline.Value{1}, sInputA.TimeVector);
        if (iBaseline(1) == iBaseline(2)) && any(iBaseline(1) == sInputA.TimeVector)
            error('Invalid baseline definition.');
        end
        iBaseline = iBaseline(1):iBaseline(2);
    % Get all file
    else
        iBaseline = 1:size(sInputA.A,2);
    end
    % Compute zscore
    sInputB.A = Compute(sInputA.A(:,iBaseline,:), sInputB.A);
    % Change DataType
    if ~strcmpi(sInputB.FileType, 'timefreq')
        sInputB.DataType = 'zscore';
    end
    % Default colormap
    if strcmpi(sInputB.FileType, 'results')
        sInputB.ColormapType = 'stat1';
        sInputB.Function = 'zscore';
    else
        sInputB.ColormapType = 'stat2';
    end
end

%% ===== COMPUTE =====
function B_data = Compute(A_baseline, B_data)
    disp('BST> process_zscore_ab.m is deprecated, use "Standardize > Baseline normalization" instead.');
    % Calculate mean and standard deviation
    [meanBaseline, stdBaseline] = process_zscore('ComputeStat', A_baseline);
    % Compute zscore
    B_data = bst_bsxfun(@minus, B_data, meanBaseline);
    B_data = bst_bsxfun(@rdivide, B_data, stdBaseline);
end

author: brainstorm-tools


text:
%
% Demo of how to load and use the dataset
%

% Parameters
nletter_per_alpha = 20; % number of letters of each alphabet to show

% Load the dataset
if ~exist('D','var')
    load('data_background','drawings','images','names','timing');
    D = Dataset(drawings,images,names,timing);
end

nrow = ceil(sqrt(nletter_per_alpha));
ncol = nrow;
nalpha = numel(D.images);

% for each alphabet
for a=1:nalpha
    msg = ['Alphabet ',D.names{a},' ',num2str(a),' of ',num2str(nalpha)];
    fprintf(1,[msg,'\n']);

    figure(1);
    clf
    for r=1:nletter_per_alpha % for each replication
       item = D.get('all',a,1,r); % get image 'r' of first character in alphabet 'a'

       % display the motor plotted on the image
       subplot(nrow,ncol,r);

       % switch between these two lines to show either image vs. motor data
       plot_motor_on_image(item.image,item.drawing,true,1);
%        plot_image_only(item.image);
    end

    % keyboard prompt
    str = input('  press <enter> to continue (''q'' to quit) : ','s');
    if str == 'q'
       break
    end
end

author: brendenlake

repo: omniglot


text:
%-----------------------------------------------------------------------
% Job configuration created by cfg_util (rev $Rev: 4899 $)
%-----------------------------------------------------------------------
matlabbatch{1}.menu_cfg.menu_entry.conf_files.type = 'cfg_files';
matlabbatch{1}.menu_cfg.menu_entry.conf_files.name = 'Directory';
matlabbatch{1}.menu_cfg.menu_entry.conf_files.tag = 'dir';
matlabbatch{1}.menu_cfg.menu_entry.conf_files.filter = 'dir';
matlabbatch{1}.menu_cfg.menu_entry.conf_files.ufilter = '.*';
matlabbatch{1}.menu_cfg.menu_entry.conf_files.dir = '';
matlabbatch{1}.menu_cfg.menu_entry.conf_files.num = [1 Inf];
matlabbatch{1}.menu_cfg.menu_entry.conf_files.check = [];
matlabbatch{1}.menu_cfg.menu_entry.conf_files.help = {'Directory to select files from.'};
matlabbatch{1}.menu_cfg.menu_entry.conf_files.def = [];
matlabbatch{2}.menu_cfg.menu_entry.conf_entry.type = 'cfg_entry';
matlabbatch{2}.menu_cfg.menu_entry.conf_entry.name = 'Filter';
matlabbatch{2}.menu_cfg.menu_entry.conf_entry.tag = 'filter';
matlabbatch{2}.menu_cfg.menu_entry.conf_entry.strtype = 's';
matlabbatch{2}.menu_cfg.menu_entry.conf_entry.extras = [];
matlabbatch{2}.menu_cfg.menu_entry.conf_entry.num = [1 Inf];
matlabbatch{2}.menu_cfg.menu_entry.conf_entry.check = [];
matlabbatch{2}.menu_cfg.menu_entry.conf_entry.help = {'A regular expression to filter files.'};
matlabbatch{2}.menu_cfg.menu_entry.conf_entry.def = [];
matlabbatch{3}.menu_cfg.menu_entry.conf_menu.type = 'cfg_menu';
matlabbatch{3}.menu_cfg.menu_entry.conf_menu.name = 'Descend into subdirectories';
matlabbatch{3}.menu_cfg.menu_entry.conf_menu.tag = 'rec';
matlabbatch{3}.menu_cfg.menu_entry.conf_menu.labels = {
                                                       'Yes'
                                                       'No'
                                                       }';
matlabbatch{3}.menu_cfg.menu_entry.conf_menu.values = {
                                                       'FPListRec'
                                                       'FPList'
                                                       }';
matlabbatch{3}.menu_cfg.menu_entry.conf_menu.check = [];
matlabbatch{3}.menu_cfg.menu_entry.conf_menu.help = {'Files can be selected from the specified directory only or from the specified directory and all its subdirectories.'};
matlabbatch{3}.menu_cfg.menu_entry.conf_menu.def = [];
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.type = 'cfg_exbranch';
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.name = 'File Selector (Batch Mode)';
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.tag = 'file_fplist';
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{1}(1) = cfg_dep;
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{1}(1).tname = 'Val Item';
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{1}(1).tgt_spec = {};
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{1}(1).sname = 'Directory (cfg_files)';
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{1}(1).src_exbranch = substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{1}(1).src_output = substruct('()',{1});
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{2}(1) = cfg_dep;
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{2}(1).tname = 'Val Item';
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{2}(1).tgt_spec = {};
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{2}(1).sname = 'Filter (cfg_entry)';
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{2}(1).src_exbranch = substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{2}(1).src_output = substruct('()',{1});
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{3}(1) = cfg_dep;
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{3}(1).tname = 'Val Item';
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{3}(1).tgt_spec = {};
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{3}(1).sname = 'Menu: Descend into subdirectories (cfg_menu)';
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{3}(1).src_exbranch = substruct('.','val', '{}',{3}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{3}(1).src_output = substruct('()',{1});
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.prog = @cfg_run_file_fplist;
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.vout = @cfg_vout_file_fplist;
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.check = [];
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.help = {'Select files from a directory using cfg_getfile(''FPList'',...).'};

author: spm


text:
function asa241_test01 ( )

%*****************************************************************************80
%
%% TEST01 tests R4_NORMAL_01_CDF_INVERSE, NORMAL_01_CDF_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 November 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01:\n' );
  fprintf ( 1, '  Let FX = NormalCDF ( X ).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  NORMAL_01_CDF_VALUES returns some values of ( X, FX ).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  R4_NORMAL_01_CDF_INVERSE takes FX and computes\n' );
  fprintf ( 1, '    an estimate X2, of the corresponding input argument,\n' );
  fprintf ( 1, '    accurate to about 7 decimal places.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X                       FX               X2\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = normal_01_cdf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    x2 = r4_normal_01_cdf_inverse ( fx );

    fprintf ( 1, '  %22.16e  %22.16e  %22.16e\n', x, fx, x2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function X = ttm(X,V,varargin)
%TTM Tensor times matrix for ktensor.
%
%   Y = TTM(X,A,N) computes the n-mode product of the ktensor X with a
%   matrix A; i.e., X x_N A.  The integer N specifies the dimension
%   (or mode) of X along which A should be multiplied.  If size(A) =
%   [J,I], then X must have size(X,N) = I.  The result will be a
%   ktensor of the same order and size as X except that size(Y,N) = J.
%
%   Y = TTM(X,{A,B,C,...}) computes the n-mode product of the ktensor
%   X with a sequence of matrices in the cell array.  The n-mode
%   products are computed sequentially along all dimensions (or modes)
%   of X. The cell array contains ndims(X) matrices.
%
%   Y = TTM(X,{A,B,C,...},DIMS) computes the sequence tensor-matrix
%   products along the dimensions specified by DIMS.
%
%   Y = TTM(...,'t') performs the same computations as above except
%   the matrices are transposed.
%
%   Examples
%   X = ktensor({rand(5,2),rand(3,2),rand(4,2),rand(2,2)});
%   A = rand(4,5); B = rand(4,3); C = rand(3,4); D = rand(3,2);
%   Y = ttm(X, A, 1)         %<-- computes X times A in mode-1
%   Y = ttm(X, {A,B,C,D}, 1) %<-- same as above
%   Y = ttm(X, A', 1, 't')   %<-- same as above
%   Y = ttm(X, {A,B,C,D}, [1 2 3 4]) %<-- 4-way multiply
%   Y = ttm(X, {D,C,B,A}, [4 3 2 1]) %<-- same as above
%   Y = ttm(X, {A,B,C,D})            %<-- same as above
%   Y = ttm(X, {A',B',C',D'}, 't')   %<-- same as above
%   Y = ttm(X, {C,D}, [3 4])     %<-- X times C in mode-3 & D in mode-4
%   Y = ttm(X, {A,B,C,D}, [3 4]) %<-- same as above
%   Y = ttm(X, {A,B,D}, [1 2 4])   %<-- 3-way multiply
%   Y = ttm(X, {A,B,C,D}, [1 2 4]) %<-- same as above
%   Y = ttm(X, {A,B,D}, -3)        %<-- same as above
%   Y = ttm(X, {A,B,C,D}, -3)      %<-- same as above
%
%   See also KTENSOR, KTENSOR/ARRANGE, TENSOR/TTM
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%%%%%%%%%%%%%%%%%%%%%%
%%% ERROR CHECKING %%%
%%%%%%%%%%%%%%%%%%%%%%

% Check the number of arguments
if (nargin < 2)
    error('TTM requires at least two arguments.');
end

% Check for transpose option
isTranspose = false;
if numel(varargin) > 0
  if isnumeric(varargin{1});
    dims = varargin{1};
  end
  isTranspose =  (ischar(varargin{end}) && (varargin{end} == 't'));
end

% Check for dims argument
if ~exist('dims','var')
    dims = [];
end

% Check that 2nd argument is cell array. If not, recall with V as a
% cell array with one element.
if ~iscell(V)
    X = ttm(X,{V},dims,varargin{end});
    return;
end

% Get sorted dims and index for multiplicands
[dims,vidx] = tt_dimscheck(dims,ndims(X),numel(V));

% Determine correct size index
if isTranspose
  j = 1;
else
  j = 2;
end

% Check that each multiplicand is the right size.
for i = 1:numel(dims)
    if (ndims(V) ~= 2) || (size(V{vidx(i)},j) ~= size(X,dims(i)))
disp(size(V{vidx(i)}))
disp(size(X))

        error('Multiplicand is wrong size');
    end
end

% Do the multiplications in the specified modes.
for i = 1:numel(dims)
  if isTranspose
    X.u{dims(i)} = V{vidx(i)}'* X.u{dims(i)};
  else
    X.u{dims(i)} = V{vidx(i)} * X.u{dims(i)};
  end
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [GAmp,GTime]=GyRadial(p)

global VCtl;
global VObj;
global VVar;

t1Start=p.t1Start;
t2Middle=p.t2Middle;
t3Start=p.t3Start;
tRamp=p.tRamp;
Gy1Sign=p.Gy1Sign;
Gy2Sign=p.Gy2Sign;
Gy3Sign=p.Gy3Sign;

% 2D radial encoding
FOV = VCtl.FOVFreq; % choose FOVFreq as real FOV
Res = VCtl.ResFreq; % choose ResFreq as real resolution

switch VCtl.R_AngPattern
    case 'Linear'
       eval(['R_AngRange=' VCtl.R_AngRange ';']);
       AngInc = R_AngRange / VCtl.R_SpokeNum;
    case 'Golden'
       AngInc = 111.246 * (pi/180); % Golden-angle sample
end

GyAmp=(1/FOV)/((VObj.Gyro/(2*pi))*(1/VCtl.BandWidth));
tHalf=1/(2*(VObj.Gyro/(2*pi))*GyAmp*(FOV/Res));
[GAmp1,GTime1]=StdTrap(t1Start-tRamp,          ...
                       t1Start+tHalf+tRamp,    ...
                       t1Start,                          ...
                       t1Start+tHalf,                    ...
                       GyAmp*sin(AngInc * (VVar.PhaseCount - 1))*Gy1Sign,2,2,2);
[GAmp2,GTime2]=StdTrap(t2Middle+VCtl.TEAnchorTime-tHalf-tRamp, ...
                       t2Middle+VCtl.TEAnchorTime+tHalf+tRamp, ...
                       t2Middle+VCtl.TEAnchorTime-tHalf,               ...
                       t2Middle+VCtl.TEAnchorTime+tHalf,               ...
                       GyAmp*sin(AngInc * (VVar.PhaseCount - 1))*Gy2Sign,2,2,2);
[GAmp3,GTime3]=StdTrap(t3Start-tRamp,            ...
                       t3Start+tHalf+tRamp,      ...
                       t3Start,                          ...
                       t3Start+tHalf,                    ...
                       GyAmp*sin(AngInc * (VVar.PhaseCount - 1))*Gy3Sign,2,2,2);

GAmp=GAmp2;
GTime=GTime2;

if Gy1Sign~=0
    GAmp=[GAmp1, GAmp];
    GTime=[GTime1, GTime];
end

if Gy3Sign~=0
    GAmp=[GAmp, GAmp3];
    GTime=[GTime, GTime3];
end

[GTime,m,n]=unique(GTime);
GAmp=GAmp(m);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcMeshc(Temp,Event,handles)
handles = guidata(handles.MU_matrix_display);

figure;
meshc(double(handles.BMatrix));
colormap(handles.V.Color_map);

end


text:
function digit = ch_to_digit ( c )

%*****************************************************************************80
%
%% CH_TO_DIGIT returns the integer value of a base 10 digit.
%
%  Example:
%
%     C   DIGIT
%    ---  -----
%    '0'    0
%    '1'    1
%    ...  ...
%    '9'    9
%    ' '    0
%    'X'   -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, the decimal digit, '0' through '9' or blank
%    are legal.
%
%    Output, integer DIGIT, the corresponding integer value.  If C was
%    'illegal', then DIGIT is -1.
%
  if ( '0' <= c & c <= '9' )

    digit = c - '0';

  elseif ( c == ' ' )

    digit = 0;

  else

    digit = -1;

  end

  return
end


text:
function tal = mrAnatGetTalairachDists
%
% tal = mrAnatGetTalairachDists
%
% Simple function to return the measured distances (in mm) from
% the original Talairach and Tournoux brain. The origin is assumed
% to be the ac. However, note that the ppc point is the distance
% from the pc- to get the most posterior distance from the ac, use
% ppc+acpc.
%
% HISTORY:
% 2005.04.29 RFD:wrote it.

tal.sac = 72;
tal.iac = -42;
tal.lac = -62;
tal.rac = 62;
tal.aac = 68;
tal.acpc = -24;
tal.ppc = -78;
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% Basis Pursuit with Douglas Rachford
% Test for DR algorithm for L1 minimization (BP).
% We do here a compressed sensing resolution
% (random matrix).

%%
% Add the toolbox.

addpath('../');
addpath('../toolbox/');

%%
% Dimensionality of the signal and number of measurements.

n = 200;
p = n/4;

%%
% Sensing matrix.

A = randn(p,n);

%%
% Measurements.

y = randn(p,1);

%%
% We aim at solving

%%
% |min_{A*x=y} norm(x,1)|

%%
% This can be rewriten as the minimization of |F(x)+G(x)|
% where |F=norm(x,1)| and |G=i_{A*x=y}| is the indicator function.

%%
% The proximity operator of the L1 norm is the soft thresholding.

ProxF = @(x,tau)perform_soft_thresholding(x, tau);

%%
% The proximity operator of the indicator of |A*x=y| is the orthogonal
% projection on A*x=y.

pA = A'*(A*A')^(-1);
ProxG = @(x,tau)x + pA*(y-A*x);

%%
% Create a function to record the values of F and the constraint at each iteration.

F = @(x)norm(x,1);
Constr = @(x)1/2*norm(y-A*x)^2;
options.report = @(x)struct('F', F(x), 'Constr', Constr(x));

%%
% Run the algorithm.

options.gamma = 5;
options.niter = 5000;
[x,R] = perform_dr(zeros(n,1), ProxF, ProxG, options);

%%
% Display the solution. At convergence, it should be of sparsity |p|.

clf;
plot(x);
axis tight;

%%
% Retrieve the F and constraint function values.

f = s2v(R,'F');
constr = s2v(R,'Constr');

%%
% Display.

clf;
subplot(2,1,1);
plot(f(2:end));
axis tight; title('Objective');
subplot(2,1,2);
plot(constr(2:end));
axis tight; title('Constraint');

author: gpeyre

repo: matlab-toolboxes


text:
function [update] = ...
    polynomial_expansion_registration3d(moving, fixed, varargin)
% POLYNOMIAL_EXPANSION_REGISTRATION3D Estimates a displacement field using polynomial expansion
%
% INPUT ARGUMENTS
% moving            - Moving image
% fixed             - Fixed image
%
% OPTIONAL INPUT ARGUMENTS
% 'signalModel'             - Local signal model to use when computing the
%                             polynomial expansion transformation
%                             'linear' (deafult), 'quadratic'
%
% 'transformationModel'     - Transformation model for estimating the
%                             displacement field
%                             translation, affine, non-rigid (default)
%
% 'multiModal'              - Set wheteher to perform multi-modal or
%                             uni-modal image registration
%                             false (default), true
%
% 'numberOfChannels'        - Number of channels to use in when computing
%                             the entropy (based on channel coding). This
%                             is only relevant if multiModal is set to
%                             true.
%                             Default value is 8
%
% OUTPUT ARGUMENTS
% update
%   displacementUpdate      - Estimated update field
%   certaintyUpdate         - Certainty related to the estimated update field
%   transformationMatrix    - Estimate transformation matrix (only if
%                             transformation model is set to translation or affine)

% Copyright (c) 2012 Daniel Forsberg
% danne.forsberg@outlook.com
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

%% Setup default parameters
% linear, quadratic
signalModel = 'linear';

% translation, affine, non-rigid
transformationModel = 'non-rigid';
multiModal = false;

% Only valid for multi-modal registration
numberOfChannels = 16;

% Only valid for non-rigid registration
sigma = 1.5;
alpha = 0.01;

% Overwrites default parameter
for k=1:2:length(varargin)
    eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

%% Perform polynomial expansion
if strcmp(signalModel,'quadratic')
    [A_moving, b_moving, c_moving] = make_Abc_fast(moving);
    [A_fixed, b_fixed, c_fixed] = make_Abc_fast(fixed);

    A = (A_moving + A_fixed)/2;
    delta_b = b_fixed - b_moving;
else
    [b_moving, c_moving] = make_bc_fast(moving);
    [b_fixed, c_fixed] = make_bc_fast(fixed);
    if multiModal
        b = b_fixed;

        [delta_c mask] = estimate_delta_c(c_fixed,c_moving,numberOfChannels);
        delta_c(mask ~= 1) = 0;
        mask = repmat(mask,[1 1 1 3]);
        b(mask ~= 1) = 0;
    else
        b = (b_moving + b_fixed)/2;

        delta_c = c_fixed - c_moving;
    end
end

if strcmp(signalModel,'quadratic')
    switch transformationModel
        case 'translation'
            [G, h] = build_G_h_quadratic3d(A, delta_b, transformationModel);
            d = G \ h;
            update.transformationMatrix = eye(4);
            update.transformationMatrix(1:3,4) = d;
        case {'rigid','affine'}
            [G, h] = build_G_h_quadratic3d(A, delta_b, transformationModel);
            p = G \ h;

            update.transformationMatrix = eye(4);
            update.transformationMatrix(1:3,1:4) = [1+p(1) p(2) p(3) p(10);...
                p(4) 1+p(5) p(6) p(11);...
                p(7) p(8) 1+p(9) p(12)];
        case 'non-rigid'
            A11 = A(:,:,:,5);
            A12 = 2*A(:,:,:,4);
            A13 = 2*A(:,:,:,8);
            A21 = 2*A(:,:,:,2);
            A22 = A(:,:,:,1);
            A23 = 2*A(:,:,:,7);
            A31 = 2*A(:,:,:,6);
            A32 = 2*A(:,:,:,3);
            A33 = A(:,:,:,9);
            B1 = delta_b(:,:,:,2);
            B2 = delta_b(:,:,:,1);
            B3 = delta_b(:,:,:,3);

            % Set the elements of the equation system
            G11 = A11.^2 + A21.^2 + A31.^2;
            G12 = A11.*A12 + A21.*A22 + A31.*A32;
            G13 = A11.*A13 + A21.*A23 + A31.*A33;
            G22 = A12.^2 + A22.^2 + A32.^2;
            G23 = A12.*A13 + A22.*A23 + A32.*A33;
            G33 = A13.^2 + A23.^2 + A33.^2;
            H1 = A11.*B1 + A21.*B2 + A31.*B3;
            H2 = A12.*B1 + A22.*B2 + A32.*B3;
            H3 = A13.*B1 + A23.*B2 + A33.*B3;

            % Smooth the elements of the equation system
            G11 = gauss_smoothing(G11, sigma);
            G12 = gauss_smoothing(G12, sigma);
            G13 = gauss_smoothing(G13, sigma);
            G22 = gauss_smoothing(G22, sigma);
            G23 = gauss_smoothing(G23, sigma);
            G33 = gauss_smoothing(G33, sigma);
            H1 = gauss_smoothing(H1, sigma);
            H2 = gauss_smoothing(H2, sigma);
            H3 = gauss_smoothing(H3, sigma);

            % Add alpha*identity matrix for stability
            scaleFactor = max([G11(:); G22(:); G33(:)]);
            G11 = G11 + alpha * scaleFactor;
            G22 = G22 + alpha * scaleFactor;
            G33 = G33 + alpha * scaleFactor;

            det = G33.*G12.^2 - 2*G12.*G13.*G23 + G22.*G13.^2 + G11.*G23.^2 - G11.*G22.*G33;

            % Estimate the displacement field
            update.displacement = cell(3,1);
            update.displacement{1} = removenan((H1.*(G23.^2 - G22.*G33) - ...
                H3.*(G12.*G23 - G13.*G22) - ...
                H2.*(G13.*G23 - G12.*G33))./(det + eps));
            update.displacement{2} = removenan((H2.*(G13.^2 - G11.*G33) - ...
                H3.*(G12.*G13 - G11.*G23) - ...
                H1.*(G13.*G23 - G12.*G33))./(det + eps));
            update.displacement{3} = removenan((H3.*(G12.^2 - G11.*G22) - ...
                H2.*(G12.*G13 - G11.*G23) - ...
                H1.*(G12.*G23 - G13.*G22))./(det + eps));

            % Estimate a certainty
            update.certainty = sqrt(b_fixed(:,:,:,2).^2 + b_fixed(:,:,:,1).^2 + b_fixed(:,:,:,3).^2);
    end
else
    switch transformationModel
        case 'translation'
            [G, h] = build_G_h_linear3d(b, delta_c, transformationModel);
            d = G \ h;
            update.transformationMatrix = eye(4);
            update.transformationMatrix(1:3,4) = d;
        case {'rigid','affine'}
            [G, h] = build_G_h_linear3d(b, delta_c, transformationModel);
            p = G \ h;

            update.transformationMatrix = eye(4);
            update.transformationMatrix(1:3,1:4) = [1+p(1) p(2) p(3) p(10);...
                p(4) 1+p(5) p(6) p(11);...
                p(7) p(8) 1+p(9) p(12)];
        case 'non-rigid'
            % Set the elements of the equation system
            % Set the elements of the equation system
            G11 = b(:,:,:,2) .* b(:,:,:,2);
            G12 = b(:,:,:,2) .* b(:,:,:,1);
            G13 = b(:,:,:,2) .* b(:,:,:,3);
            G22 = b(:,:,:,1) .* b(:,:,:,1);
            G23 = b(:,:,:,1) .* b(:,:,:,3);
            G33 = b(:,:,:,3) .* b(:,:,:,3);
            H1 = b(:,:,:,2) .* delta_c;
            H2 = b(:,:,:,1) .* delta_c;
            H3 = b(:,:,:,3) .* delta_c;

            % Smooth the elements of the equation system
            G11 = gauss_smoothing(G11, sigma);
            G12 = gauss_smoothing(G12, sigma);
            G13 = gauss_smoothing(G13, sigma);
            G22 = gauss_smoothing(G22, sigma);
            G23 = gauss_smoothing(G23, sigma);
            G33 = gauss_smoothing(G33, sigma);
            H1 = gauss_smoothing(H1, sigma);
            H2 = gauss_smoothing(H2, sigma);
            H3 = gauss_smoothing(H3, sigma);

            % Add alpha*identity matrix for stability
            scaleFactor = max([G11(:); G22(:); G33(:)]);
            G11 = G11 + alpha * scaleFactor;
            G22 = G22 + alpha * scaleFactor;
            G33 = G33 + alpha * scaleFactor;

            det = G33.*G12.^2 - 2*G12.*G13.*G23 + G22.*G13.^2 + G11.*G23.^2 - G11.*G22.*G33;

            % Estimate the displacement field
            update.displacement = cell(3,1);
            update.displacement{1} = removenan((H1.*(G23.^2 - G22.*G33) - ...
                H3.*(G12.*G23 - G13.*G22) - ...
                H2.*(G13.*G23 - G12.*G33))./(det + eps));
            update.displacement{2} = removenan((H2.*(G13.^2 - G11.*G33) - ...
                H3.*(G12.*G13 - G11.*G23) - ...
                H1.*(G13.*G23 - G12.*G33))./(det + eps));
            update.displacement{3} = removenan((H3.*(G12.^2 - G11.*G22) - ...
                H2.*(G12.*G13 - G11.*G23) - ...
                H1.*(G12.*G23 - G13.*G22))./(det + eps));

            % Estimate a certainty
            update.certainty = sqrt(b(:,:,:,2).^2 + b(:,:,:,1).^2 + b(:,:,:,3).^2);
    end
end

author: fordanic

repo: image-registration


text:
function logplot(x_orig, logscale)
% Normal plot with layered logscale plot
%   logplot(x_orig, logscale)
%     x_orig    Plotted data
%     logscale  Optional 2-element vector representing plot range
%               for cropping logarithmic plot area

% Copyright (C) 2004, 2005 DSS MATLAB package team (dss@cis.hut.fi).
% Distributed by Laboratory of Computer and Information Science,
% Helsinki University of Technology. http://www.cis.hut.fi/projects/dss/.
% $Id$

if size(x_orig,1)==1
  x=x_orig';
else
  x=x_orig;
end

if size(x,1)==0; return; end
i=[1:size(x,1)];

x(x==0)=NaN;
y=log(x);
y(isinf(y))=NaN;

if nargin<2
  offset = min(x)-min(y);
  scale = (max(x)-min(x))/(max(y)-min(y));
else
  offset = logscale(1)-min(y);
  scale = (logscale(2)-logscale(1))/(max(y)-min(y));
end
plot(i, x, '-', i, (y-min(y))*scale+min(x), ':');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function field_values = get_field(db, field_name)
    NUM_OF_ENTRIES = length(db);
    field_values = cell(NUM_OF_ENTRIES, 1);

    for k = 1:NUM_OF_ENTRIES
        metadata = getmeta(db, k);
        field_values{k} = metadata.(field_name);
    end;

    [field_values, ~] = unique(field_values);
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function dist = codedist_loss(C1,C2,Ylat,loss)
% compute the distance using a loss function metric between 'C1' and  'C2'
%
%    dist = codedist_loss(C1,C2,Ylatent)
%    dist = codedist_loss(C1,C2,Ylatent, loss_fct)
%
% 'C1' contains the result code,
% 'C2' contains the codebooks prototype. An infine small number 'eps'  represents the don't care
% 'loss_fct' is the loss function used in order to compute the loss
%            between 2 codewords. By default the sum of squares is used.
%            One can do 'winner-takes-all' decoding by using the
%            loss function 'max'.
%
% An example:
%  >> Ye = [1 1; 1 1; -1 -1; 1 -1];
%  >> codebook = [1 2 3];
%  >> old_codebook = [1 1 -1; 1 -1 -1];
%  >> code(Ye, codebook, [],old_codebook,'codedist_loss',{'Ylatent','mse'})
%
% To use this distance measure in LS-SVMlab, the following
% procedure is to be followed, assume input data 'X' and multiclass
% output 'Y'
%
% % encode for training
% >> model = initlssvm(X,Y,'classification',gam,sig2,'preprocess','RBF_kernel');
% >> model = changelssvm(model,'codetype','code_OneVsOne');
% >> model = trainlssvm(model);
%
% % decode for simulating
% >> [Yhamming, Ylatent] = simlssvm(model,Xt);
% >> model = changelssvm(model,'codedist_fct','codedist_loss');
% >> model = changelssvm(model,'codedist_args',{Ylatent,'sse'});
% >>   Yt  = simlssvm(model,Xt);
%
% see also:
%   bay_modoutClass, codedist_hamming, code_ECOC

% (c) SCD-KULeuven, rights & help @ http://www.esat.kuleuven.ac.be/sista/lssvmlab

if nargin<3,
  warning('The latent variables should be used, proceeding with the binary classifiers...');
  Ylat = C1;
end
eval('loss;','loss=''sse'';');

[nb,nbin] = size(Ylat);
[~,dim] = size(C2);
dist = zeros(nb,dim);

for d = 1:dim,
  for n= 1:nb,
    nondontcare = find(Ylat(n,:)~=eps & C2(:,d)'~=eps);
    dist(n,d) = feval(loss, Ylat(n,nondontcare),C2(nondontcare,d)');
  end
end
dist

function l = sse(X,Y)
l = sum(sum((X-Y).^2));

function l = winnertakesall(X,Y)
p = find(Y>0);
l = max(X(P));

author: peterhcharlton

repo: RRest


text:
function [C_smooth]=patchSmoothFaceMeasure(varargin)

% function [C_smooth]=patchSmoothFaceMeasure(F,V,C,smoothPar)

%% Parse input
switch nargin
    case 3
        F=varargin{1};
        V=varargin{2};
        C=varargin{3};
        smoothPar=[];
    case 4
        F=varargin{1};
        V=varargin{2};
        C=varargin{3};
        smoothPar=varargin{4};
end

smoothParDefault.lambda=0.5;
smoothParDefault.n=1;
smoothPar=structComplete(smoothPar,smoothParDefault,1);

%% Get connectivity array

[connectivityStruct]=patchConnectivity(F,V);
faceFaceConnectivity=connectivityStruct.face.face;

%%

nDims=size(C,2); %Number of dimensions
logicValid=faceFaceConnectivity>0;
C_smooth=C;
C_smooth_step=C;
for qIter=1:smoothPar.n
    %Loop for all dimensions
    for qDim=1:1:nDims
        Xp=NaN(size(C,1),size(faceFaceConnectivity,2));
        Xp(logicValid)=C_smooth(faceFaceConnectivity(logicValid),qDim);
        Xp=nanmean(Xp,2);
        C_smooth_step(:,qDim)=Xp;
    end
    C_smooth=((1-smoothPar.lambda).*C_smooth)+(smoothPar.lambda.*C_smooth_step);
end

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
%  INTERNAL FUNCTION
%
%  Note:
%       At present the restrictions are enforced only in the violating periods,
%       which implies that the contemporaneous shocks are altered. But one can
%       think of a scheme whereby only future shocks are changed, keeping the
%       original/actual shocks unchanged.
%
%

author: jmaih


text:
function [C, details] = spr_cvx(Y, options)
%%% validate arguments
if nargin < 2
    options = struct;
end
% verbosity
verbose = 0;
if isfield(options, 'verbose')
    verbose = options.verbose;
end
% noise factor
noise_factor = 0.01;
if isfield(options, 'noise_factor')
    noise_factor = options.noise_factor;
end
% affine
affine = false;
if isfield(options, 'affine')
    affine = options.affine;
end
% size of data-set
[M, S] = size(Y);
% iterate over signals
all_cols = 1:S;
% The representation matrix
C = zeros(S);
details = struct;
for s=all_cols
    if verbose > 0
        fprintf('.');
        if mod(s, 50) == 0
            fprintf('\n');
        end
    end
    % s-th data vector
    y = Y(:, s);
    % the signal dictionary to be used for reconstruction
    cols = all_cols ~= s; % S - 1 columns
    % dimensions are D x (S - 1)
    A = Y(:, cols);
    % threshold for the noise norm
    noise_norm_threshold = noise_factor * norm(y);
    if affine
        cvx_begin quiet;
            cvx_precision high;
            variable rep(S-1,1);
            minimize( norm(rep,1) );
            subject to
                norm(A * rep  - y) <= noise_norm_threshold;
                sum(rep) == 1;
        cvx_end;
    else
        cvx_begin quiet;
            cvx_precision high;
            variable rep(S-1,1);
            minimize( norm(rep,1) );
            subject to
                norm(A * rep  - y) <= noise_norm_threshold;
        cvx_end;
    end
    % Prepare the l1 solver
    % solver = spx.pursuit.single.BasisPursuit(A, x);
    % solver.Quiet = true;
    % lambda = 2.5;
    % Run the solver to obtain sparse representation
    % rep = solver.solve_l1_noise();
    % The obtained representation is in R^{S - 1} dimensions

    % Put back the representation
    C(cols, s) = rep;
end
if verbose > 0
    fprintf('\n');
end

end

author: indigits

repo: sparse-plex


text:
% POP_EDITEVENTVALS - Edit events contained in an EEG dataset structure.
%               If the dataset is the only input, a window pops up
%               allowing the user to insert the relevant parameter values.
%
% Usage: >> EEGOUT = pop_editeventvals( EEG, 'key1', value1, ...
%                                    'key2', value2, ... );
% Input:
%   EEG  - EEG dataset
%
% Optional inputs:
%   'sort'        - { field1 dir1 field2 dir2 } Sort events based on field1
%                   then on optional field2. Arg dir1 indicates the sort
%                   direction (0 = increasing, 1 = decreasing).
%   'changefield' - {num field value} Insert the given value into the specified
%                   field in event num. (Ex: {34 'latency' 320.4}). Note
%                   that for latency event latencies are in millisecond
%                   compared to the time locking event not in data samples.
%   'changeevent' - {num value1 value2 value3 ...} Change the values of
%                   all fields in event num.
%   'add','append','insert' - {num value1 value2 value3 ...} Insert event
%                   before or at event num, and assign value to structure
%                   fields. Note that the latency field must be in second
%                   and will be converted to data sample. Note also that
%                   the index of the event is often irrelevant, as events
%                   will be automatically resorted by latencies.
%   'delete'      - vector of indices of events to delete
%
% Outputs:
%   EEGOUT        - EEG dataset with the selected events only
%
% Ex:  EEG = pop_editeventvals(EEG,'changefield', { 1 'type' 'target'});
%        % set field type of event number 1 to 'target'
%
% Author: Arnaud Delorme & Hilit Serby, SCCN, UCSD, 15 March 2002
%
% See also: POP_SELECTEVENT, POP_IMPORTEVENT

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

% 03-16-02 text interface editing -sm & ad
% 03-18-02 automatic latency switching display (epoch/continuous) - ad & sm
% 03-18-02 debug soring order - ad
% 03-18-02 put latencies in ms - ad, lf & sm
% 03-29-02 debug latencies in ms - ad & sm
% 04-02-02 debugging test - ad & sm

function [EEG, com] = pop_editeventvals(EEG, varargin)

com ='';
if nargin < 1
   help pop_editeventvals;
   return;
end

if nargin >= 2 || ischar(EEG) % interpreting command from GUI or command line

    if ischar(EEG) % GUI
        gui    = 1;
        varargin = { EEG varargin{:} };

        % user data
        % ---------
        userdata  = get(gcf, 'userdata');
        EEG       = userdata{1};
        oldcom    = userdata{2};
        allfields = fieldnames(EEG.event);
        tmpind    = strmatch('urevent', allfields);
        allfields(tmpind) = [];

        % current event
        % -------------
        objevent  = findobj('parent', gcf, 'tag', 'numval');
        valnum    = str2double(get(objevent, 'string'));
        shift     = 0;

    else % command line

        gui = 0;
        if isempty(EEG.event)
            disp('Getevent: cannot deal with empty event structure');
            return;
        end

        allfields = fieldnames(EEG.event);
        tmpind = strmatch('urevent', allfields);
        allfields(tmpind) = [];

    end

    % retinterpret inputs (fix BUG 454)
    % --------------------------------
    if ~gui
        newvararg = {};
        for indfield = 1:2:length(varargin)

            com     = varargin{ indfield   };
            tmpargs = varargin{ indfield+1 };
            newvararg = { newvararg{:} com };

            if any(strcmpi({'add','insert','append'},com))
                evtind     = tmpargs{1};
                fields     = fieldnames(EEG.event);
                emptycells = cell(1,length(fields)-1);
                newvararg  = { newvararg{:}, { evtind emptycells{:} } };
                if strcmpi(com, 'append'), evtind = evtind+1; end

                for ind = 2:length( tmpargs )
                    if ~strcmpi(fields{ind-1}, 'urevent')
                        newvararg = { newvararg{:},'changefield',{ evtind fields{ind-1} tmpargs{ind} } };
                    end
                end
            else
                newvararg = { newvararg{:} tmpargs };
            end
        end
        varargin = newvararg;
    end

    % scan inputs
    % -----------
    for indfield = 1:2:length(varargin)
        if length(varargin) >= indfield+1
            tmparg = varargin{ indfield+1 };
        end

        switch lower(varargin{indfield})

     case 'goto' % ******************** GUI ONLY ***********************

      % shift time
      % ----------
      shift     = tmparg;
      valnum    = valnum + shift;
      if valnum < 1,                valnum = 1;                end
      if valnum > length(EEG.event), valnum = length(EEG.event); end
      set(objevent, 'string', num2str(valnum,5));

      % update fields
      % -------------
      for index = 1:length(allfields)

          enable = 'on';
          if isfield(EEG.event, 'type')
              if strcmpi(EEG.event(valnum).type, 'boundary'), enable = 'off'; end
          end
          if strcmp( allfields{index}, 'latency') && ~isempty(EEG.event(valnum).latency)
              if isfield(EEG.event, 'epoch')
                   value = eeg_point2lat( EEG.event(valnum).latency, EEG.event(valnum).epoch, ...
                                          EEG.srate, [EEG.xmin EEG.xmax]*1000, 1E-3);
              else value = (EEG.event(valnum).latency-1)/EEG.srate+EEG.xmin;
              end
          elseif strcmp( allfields{index}, 'duration') && ~isempty(EEG.event(valnum).duration)
              if isfield(EEG.event, 'epoch')
                   value = EEG.event(valnum).duration/EEG.srate*1000; % milliseconds
              else value = EEG.event(valnum).duration/EEG.srate;      % seconds
              end
          else
              value = getfield( EEG.event(valnum), allfields{index});
          end

          % update interface
          % ----------------
          tmpobj = findobj('parent', gcf, 'tag', allfields{index});
          if isnumeric(value)
              set(tmpobj, 'string', num2str(value,5), 'enable', enable);
          else % for strings and cells
              set(tmpobj, 'string', value, 'enable', enable);
          end

      end

      % update original
      % ---------------
      tmpobj = findobj('parent', gcf, 'tag', 'original');
      if isfield(EEG.event, 'urevent') && ~isempty(EEG.event(valnum).urevent) && EEG.event(valnum).urevent ~= valnum
           set(tmpobj, 'string', [ 'originally ' int2str(EEG.event(valnum).urevent)], ...
                       'horizontalalignment', 'center');
      else set(tmpobj, 'string', ' ');
      end
      return; % NO NEED TO SAVE ANYTHING

     case { 'append' 'insert' 'add' } % **********************************************************

      if gui
          shift     = tmparg; % shift is for adding before or after the event

          % add epoch number if data epoch
          % ------------------------------
          tmpcell    = cell(1,1+length(fieldnames(EEG.event)));
          tmpcell{1} = valnum;
          if EEG.trials > 1
              indepoch = strmatch('epoch', fieldnames(EEG.event), 'exact');
              if valnum > 1, tmpprevval = valnum-1;
              else           tmpprevval = valnum+1;
              end
              if tmpprevval <= length(EEG.event)
                  tmpcell{indepoch+1} = EEG.event(tmpprevval).epoch;
              end
          end

          % update commands
          % ---------------
          if shift
              oldcom     = { oldcom{:} 'append', tmpcell };
          else
              oldcom     = { oldcom{:} 'insert', tmpcell };
          end
      else
          if strcmpi(lower(varargin{indfield}), 'append') % not 'add' for backward compatibility
               shift = 1;
          else shift = 0;
          end
          valnum = tmparg{1};
      end

      % find ur index
      % -------------
      valnum    = valnum   + shift;
      if isfield(EEG.event, 'epoch'), curepoch = EEG.event(valnum).epoch; end
      if isfield(EEG, 'urevent') && isfield(EEG.event, 'urevent')
          % find non empty urvalnum
          urvalnum = [];
          count = 0;
          while isempty(urvalnum)
              tmpindex = mod(valnum+count-1, length(EEG.event)+1)+1;
              urvalnum = EEG.event(valnum+count).urevent;
              count = count+1;
          end
          if isfield(EEG.urevent, 'epoch'), urcurepoch = EEG.urevent(urvalnum).epoch; end
          urvalnum = urvalnum;
      end

      % update urevents
      % ---------------
      if isfield(EEG, 'urevent') && isfield(EEG.event, 'urevent')
          EEG.urevent(end+3)            = EEG.urevent(end);
          EEG.urevent(urvalnum+1:end-2) = EEG.urevent(urvalnum:end-3);
          EEG.urevent(urvalnum)         = EEG.urevent(end-1);
          EEG.urevent                   = EEG.urevent(1:end-2);
          if isfield(EEG.urevent, 'epoch'),  EEG.urevent(urvalnum).epoch = urcurepoch; end
      end

      % update events
      % -------------
      EEG.event(end+3)          = EEG.event(end);
      EEG.event(valnum+1:end-2) = EEG.event(valnum:end-3);
      EEG.event(valnum)         = EEG.event(end-1);
      EEG.event                 = EEG.event(1:end-2);
      if isfield(EEG.event, 'epoch'), EEG.event(valnum).epoch = curepoch; end;
      if isfield(EEG, 'urevent') && isfield(EEG.event, 'urevent')
          EEG.event(valnum).urevent = urvalnum;
          for index = valnum+1:length(EEG.event)
              EEG.event(index).urevent = EEG.event(index).urevent+1;
          end
      end

      % update type field
      % -----------------
      for tmpind = 1:length(allfields)
          EEG.event = checkconsistency(EEG.event, valnum, allfields{tmpind});
      end

     case 'delete' % **********************************************************

      if ~gui
          valnum = tmparg;
      end

      EEG.event(valnum) = [];

      if gui

          valnum           = min(valnum,length(EEG.event));
          set(objevent, 'string', num2str(valnum));

          % update commands
          % ---------------
          oldcom = { oldcom{:} 'delete', valnum };

      end

     case { 'assign' 'changefield' } % **********************************************************

      if gui % GUI case

          field    = tmparg;
          objfield = findobj('parent', gcf, 'tag', field);
          editval     = get(objfield, 'string');
          if ~isempty(editval) && ~isnan(str2double(editval)), editval = str2double(editval); end

          % update history
          % --------------
          oldcom = { oldcom{:},'changefield',{ valnum field editval }};

      else % command line case

          valnum  = tmparg{1};
          field   = tmparg{2};
          editval = tmparg{3};

      end

      % latency and duration case
      % -------------------------
      if strcmp( field, 'latency') && ~isempty(editval)
          if isfield(EEG.event, 'epoch')
               editval = eeg_lat2point( editval, EEG.event(valnum).epoch, ...
                                       EEG.srate, [EEG.xmin EEG.xmax]*1000, 1E-3);
          else editval = (editval- EEG.xmin)*EEG.srate+1;
          end
      end
      if strcmp( field, 'duration') && ~isempty(editval)
          if isfield(EEG.event, 'epoch')
               editval = editval/1000*EEG.srate; % milliseconds
          else editval = editval*EEG.srate;      % seconds
          end
      end

      % adapt to other formats
      % ----------------------
      EEG.event(valnum) = setfield(EEG.event(valnum), field, editval);
      EEG.event = checkconsistency(EEG.event, valnum, field);

      % update urevents
      % ---------------
      if isfield(EEG, 'urevent') && isfield(EEG.event, 'urevent')
          urvalnum  = EEG.event(valnum).urevent;

          % latency case
          % ------------
          if strcmp( field, 'latency') && ~isempty(editval)
              if isfield(EEG.urevent, 'epoch')
                  urepoch = EEG.urevent(urvalnum).epoch;

                  % find closest event latency
                  % --------------------------
                  if valnum<length(EEG.event)
                      if EEG.event(valnum+1).epoch == urepoch
                          urlatency = EEG.urevent(EEG.event(valnum+1).urevent).latency;
                          latency   = EEG.event(valnum+1).latency;
                      end
                  end
                  if valnum>1
                      if EEG.event(valnum-1).epoch == urepoch
                          urlatency = EEG.urevent(EEG.event(valnum-1).urevent).latency;
                          latency   = EEG.event(valnum-1).latency;
                      end
                  end

                  % update event
                  % ------------
                  if exist('urlatency') ~=1
                      disp('Urevent not updated: could not find other event in the epoch');
                  else
                      editval = urlatency - ( latency - editval ); % new latency value
                  end
              else
                  editval = eeg_urlatency(EEG.event, EEG.event(valnum).latency);
              end
          elseif strcmp( field, 'latency') % empty editval
              EEG.event(valnum).latency = NaN;
          end

          % duration case
          % ------------
          if strcmp( field, 'duration') && ~isempty(editval)
              if isfield(EEG.event, 'epoch')
                   editval = editval/1000*EEG.srate; % milliseconds -> point
              else editval = editval*EEG.srate;      % seconds -> point
              end
          end
          EEG.urevent = setfield(EEG.urevent, {urvalnum}, field, editval);
      end

     case 'sort' % **********************************************************

      if gui % retrieve data

          field1 = get(findobj('parent', gcf, 'tag', 'listbox1'), 'value');
          field2 = get(findobj('parent', gcf, 'tag', 'listbox2'), 'value');
          dir1   = get(findobj('parent', gcf, 'tag', 'order1'),   'value');
          dir2   = get(findobj('parent', gcf, 'tag', 'order2'),   'value');

          if field1 > 1, field1 = allfields{field1-1}; else return; end
          if field2 > 1, field1 = allfields{field2-1}; else field2 = []; end

          % update history
          % --------------
          oldcom = { oldcom{:},'sort',{ field1 dir1 field2 dir2 } };

      else % command line
          field1 = tmparg{1};
          if length(tmparg) < 2, dir1 = 0;
          else                   dir1 = tmparg{2};
          end
          if length(tmparg) < 3, field2 = [];
          else                   field2 = tmparg{3};
          end
          if length(tmparg) < 4, dir2 = 0;
          else                   dir2 = tmparg{4};
          end
      end

      % Toby edit 11/16/2005 This section is scrambling the eeg.event
      % fields. Requires further investigation.

      if ~isempty(field2)
          tmpevent = EEG.event;
          if ~ischar(EEG.event(1).(field2))
              tmparray = [ tmpevent.(field2) ];
          else tmparray = { tmpevent.(field2) };
          end
          % Commented out 11/18/2005, Toby
          % These lines were incorrectly sorting the event.latency field in
          % units of time (seconds) relevant to each event's relative
          % latency time as measured from the start of each epoch. It is
          % possible that there are occasions when it is desirable to do
          % so, but in the case of pop_mergset() it is not.

          %if strcmp(field2, 'latency') & EEG.trials > 1
          %    tmparray = eeg_point2lat(tmparray, {EEG.event.epoch}, EEG.srate, [EEG.xmin EEG.xmax], 1);
          %end
          [~, I] = mysort( tmparray );
          if dir2 == 1, I = I(end:-1:1); end
          events = EEG.event(I);
      else
          events = EEG.event;
      end
      tmpevent = EEG.event;
      if ~ischar(EEG.event(1).(field1))
           tmparray = [ tmpevent.(field1) ];
      else tmparray = { tmpevent.(field1) };
      end
      % Commented out 11/18/2005, Toby
      %if strcmp( field1, 'latency') & EEG.trials > 1
      %    tmparray = eeg_point2lat(tmparray, {events.epoch}, EEG.srate, [EEG.xmin EEG.xmax], 1);
      %end
      [~, I] = mysort( tmparray );
      if dir1 == 1, I = I(end:-1:1); end
      EEG.event = events(I);

      if gui
          % warn user
          % ---------
          warndlg2('Sorting done');
      else
          noeventcheck  = 1; % otherwise infinite recursion with eeg_checkset
      end

        end % end switch
    end % end loop

    % save userdata
    % -------------
    if gui
        userdata{1} = EEG;
        userdata{2} = oldcom;
        set(gcf, 'userdata', userdata);
        pop_editeventvals('goto', shift);
    else
        if ~exist('noeventcheck','var')
            EEG = eeg_checkset(EEG, 'eventconsistency');
            EEG = eeg_checkset(EEG, 'checkur');
        end
        if nargout > 1
            com = sprintf('EEG = pop_editeventvals(EEG, %s);', vararg2str(varargin));
        end
    end
    return;
end

% ----------------------
% graphic interface part
% ----------------------

if isempty(EEG.event)
    disp('Getevent: cannot deal with empty event structure');
    return;
end

allfields = fieldnames(EEG.event);
tmpind = strmatch('urevent', allfields);
allfields(tmpind) = [];

if nargin<2
    % add field values
    % ----------------
    geometry = { [2 0.5] };
    tmpstr = sprintf('Edit event field values (currently %d events)',length(EEG.event));
    uilist = { { 'Style', 'text', 'string', tmpstr, 'fontweight', 'bold'  } ...
               { 'Style', 'pushbutton', 'string', 'Delete event', 'callback', 'pop_editeventvals(''delete'');'  }};

    for index = 1:length(allfields)

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

        % input string
        % ------------
        if strcmp( allfields{index}, 'latency') || strcmp( allfields{index}, 'duration')
            if EEG.trials > 1
                 inputstr =  [ allfields{index} ' (ms)'];
            else inputstr =  [ allfields{index} ' (sec)'];
            end
		else inputstr =  allfields{index};
		end

		% callback for displaying help
		% ----------------------------
        if index <= length( EEG.eventdescription )
             tmptext = EEG.eventdescription{ index };
			 if ~isempty(tmptext)
				 if size(tmptext,1) > 15,    stringtext = [ tmptext(1,1:15) '...' ];
				 else                        stringtext = tmptext(1,:);
				 end
			 else stringtext = 'no-description'; tmptext = 'no-description';
			 end
        else stringtext = 'no-description'; tmptext = 'no-description';
        end
		cbbutton = ['questdlg2(' vararg2str(tmptext) ...
					',''Description of field ''''' allfields{index} ''''''', ''OK'', ''OK'');' ];

        % create control
        % --------------
        cbedit = [ 'pop_editeventvals(''assign'', ''' allfields{index} ''');' ];
		uilist   = { uilist{:}, { }, ...
					 { 'Style', 'pushbutton', 'string', inputstr, 'callback',cbbutton  }, ...
					 { 'Style', 'edit', 'tag', allfields{index}, 'string', '', 'callback', cbedit } ...
                     { } };
    end

    % add buttons
    % -----------
    geometry = { geometry{:} [1] [1.2 0.6 0.6 1 0.6 0.6 1.2] [1.2 0.6 0.6 1 0.6 0.6 1.2] [2 1 2] };

    tpappend = 'Append event after the current event';
    tpinsert = 'Insert event before the current event';
    tporigin = 'Original index of the event (in EEG.urevent table)';
    uilist   = { uilist{:}, ...
          { }, ...
          { },{ },{ }, {'Style', 'text', 'string', 'Event Num', 'fontweight', 'bold' }, { },{ },{ }, ...
          { 'Style', 'pushbutton', 'string', 'Insert event',  'callback', 'pop_editeventvals(''append'', 0);', 'tooltipstring', tpinsert }, ...
          { 'Style', 'pushbutton', 'string', '<<',            'callback', 'pop_editeventvals(''goto'', -10);' }, ...
          { 'Style', 'pushbutton', 'string', '<',             'callback', 'pop_editeventvals(''goto'', -1);' }, ...
          { 'Style', 'edit',       'string', '1',             'callback', 'pop_editeventvals(''goto'', 0);', 'tag', 'numval' }, ...
          { 'Style', 'pushbutton', 'string', '>',             'callback', 'pop_editeventvals(''goto'', 1);' }, ...
          { 'Style', 'pushbutton', 'string', '>>',            'callback', 'pop_editeventvals(''goto'', 10);' }, ...
          { 'Style', 'pushbutton', 'string', 'Append event',  'callback', 'pop_editeventvals(''append'', 1);', 'tooltipstring', tpappend }, ...
          { }, { 'Style', 'text',  'string', ' ', 'tag', 'original' 'horizontalalignment' 'center' 'tooltipstring' tporigin } { } };

    % add sorting options
    % -------------------
    listboxtext = 'No field selected';
    for index = 1:length(allfields)
         listboxtext = [ listboxtext '|' allfields{index} ];
    end

    userdata = { EEG {} };
    inputgui( geometry, uilist, 'pophelp(''pop_editeventvals'');', ...
                                  'Edit event values -- pop_editeventvals()', userdata, 'plot');
    pop_editeventvals('goto', 0);

    % wait for figure
    % ---------------
    fig = gcf;
    waitfor( findobj('parent', fig, 'tag', 'ok'), 'userdata');
    try, userdata = get(fig, 'userdata'); close(fig); % figure still exist ?
    catch, return; end

    % transfer events
    % ---------------
    if ~isempty(userdata{2})
        com = sprintf('EEG = pop_editeventvals(EEG,%s);', vararg2str(userdata{2}));
    end
    if isempty(findstr('''sort''', com))
        if ~isempty(userdata{2}) % some modification have been done
            EEG       = userdata{1};
            disp('Checking event consistency...');
            EEG = eeg_checkset(EEG, 'eventconsistency');
            EEG = eeg_checkset(EEG, 'checkur');
        end
    else
        com = '';
        disp('WARNING: all edits discarded because of event resorting. The EEGLAB event structure');
        disp('            must contain events sorted by latency (you may obtain an event structure');
        disp('            with resorted event by calling this function from the command line).');
    end
    return;

end
return;

% format the output field
% -----------------------
function strval = reformat( val, latencycondition, trialcondition, eventnum)
    if latencycondition
        if trialcondition
            strval = ['eeg_lat2point(' num2str(val) ', EEG.event(' int2str(eventnum) ').epoch, EEG.srate,[EEG.xmin EEG.xmax]*1000, 1E-3);' ];
        else
            strval = [ '(' num2str(val) '-EEG.xmin)*EEG.srate+1;' ];
        end
    else
        if ischar(val), strval = [ '''' val '''' ];
        else            strval = num2str(val);
        end
    end

% sort also empty values
% ----------------------
function [X, I] = mysort(tmparray)
    if iscell(tmparray)
        if all(cellfun('isreal', tmparray))
            tmpempty = cellfun('isempty', tmparray);
            tmparray(tmpempty) = { 0 };
            tmparray = [ tmparray{:} ];
        end
    end
    try
        [X, I] = sort(tmparray);
    catch
        disp('Sorting failed. Check that selected fields contain uniform value format.');
        X = tmparray;
        I = 1:length(X);
    end

% checkconsistency of new event
% -----------------------------
function eventtmp = checkconsistency(eventtmp, valnum, field)

    otherval = mod(valnum+1, length(eventtmp))+1;

    if ischar(getfield(eventtmp(valnum), field)) && ~ischar(getfield(eventtmp(otherval), field))
        eventtmp(valnum) = setfield(eventtmp(valnum), field, str2double(getfield(eventtmp(valnum), field)));
    end
    if ~ischar(getfield(eventtmp(valnum), field)) && ischar(getfield(eventtmp(otherval), field))
        eventtmp(valnum) = setfield(eventtmp(valnum), field, num2str(getfield(eventtmp(valnum), field)));
    end
    if strcmpi(field, 'latency') && isempty(getfield(eventtmp(valnum), field))
        eventtmp(valnum).latency = NaN;
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function write_serial_event(filename, event)

% WRITE_SERIAL_EVENT
%
% changed A.Hadjipapas 2010
%
% write to phyiscal serial port
% serial port on windows or linux platform

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

[port, opt] = filetype_check_uri(filename);
temp=instrfind; %% find isntruments including serial ports

%% check if a serial port with the specified name is already available
if isempty(temp)||~any(strcmpi(temp.Port,port))
  %% if not make one
  s=serial(port);
  %% if user has specified Baudrate, apply that, %%  FIX THIS: currently no
  %% automatic implementation for more options
  if ~isempty(opt)
    s.Baudrate =str2double(opt{2}); %% set Baudrate using optional parameters
  end
  fopen(s); %% open serial port
end

%% write the contents of the field event.value to the serial port as a string
if isfield(event,'value') && ~isempty(event.value)
  msg = event.value;
else
  msg = [];
end

%% write the message to the serial port if it is open
if strcmp(s.Status,'open')
  %% if numeric use fwrite otherwise fprintf
  if isnumeric(msg)
    fwrite(s,msg);
  end

  if ischar(msg)
    fprintf(s,'%s',msg);
  end
else
  error('could not write event to serial port');
end

%% this is the old code
%
%     % serial port on windows or linux platform
%     s = [];
%     [port, opt] = filetype_check_uri(filename);
%
%     % determine whether any serial port objects are already associated with the target serial port
%     temp = instrfind;
%     if isa(temp,'instrument')
%       % find all serial ports
%       i1 = strcmp('serial',{temp(:).Type});
%       if any(i1)
%         % find all serial ports whose name matches that of the specified port
%         i2 = strmatch(lower(['Serial-',port]),lower({temp(find(i1)==1).Name}));
%         % set s to the (first) matching port if present (and open if necessary)
%         if ~isempty(i2)
%           s = temp(i2(1));
%           if ~strcmp(s.Status,'open'), fopen(s); end;
%         end
%       end
%     end
%
%     % create, configure a serial port object if necessary and open the port
%     if ~isa(s,'serial')
%       s = serial(port);
%       if ~isempty(opt) && iscell(opt), s = set(s,opt); end;
%       fopen(s);
%     end
%
%     %     % convert the event structure into an appropriate message
%     %     if isfield(event,'type') && strcmp(event.type,'ctrlchar')
%     %       % use only a single control character
%     %       msg = char(event.value(1));
%     %     else
%     %       % convert the entire event structure into a message
%     %       msg = struct2msg(event);
%     %     end
%
%     % write the message to the serial port
%     if isa(s,'serial') && strcmp(s.Status,'open')
%       fwrite(s,char(event.value));
%     else
%       error('could not write event to serial port');
%     end

author: PatternRecognition

repo: OpenBMI


text:
clear all;
close all;
clc;
check_hopkins(@ssc_l1, 'ssc_l1');


text:
function [dpsi,deps]=nut_iau1980(t,f)

AS2R=pi/180/3600;
nut=[   0,   0,   0,   0,   1, -6798.4, -171996, -174.2, 92025,   8.9;
        0,   0,   2,  -2,   2,   182.6,  -13187,   -1.6,  5736,  -3.1;
        0,   0,   2,   0,   2,    13.7,   -2274,   -0.2,   977,  -0.5;
        0,   0,   0,   0,   2, -3399.2,    2062,    0.2,  -895,   0.5;
        0,  -1,   0,   0,   0,  -365.3,   -1426,    3.4,    54,  -0.1;
        1,   0,   0,   0,   0,    27.6,     712,    0.1,    -7,   0.0;
        0,   1,   2,  -2,   2,   121.7,    -517,    1.2,   224,  -0.6;
        0,   0,   2,   0,   1,    13.6,    -386,   -0.4,   200,   0.0;
        1,   0,   2,   0,   2,     9.1,    -301,    0.0,   129,  -0.1;
        0,  -1,   2,  -2,   2,   365.2,     217,   -0.5,   -95,   0.3;
        -1,   0,   0,   2,   0,    31.8,     158,    0.0,    -1,   0.0;
        0,   0,   2,  -2,   1,   177.8,     129,    0.1,   -70,   0.0;
        -1,   0,   2,   0,   2,    27.1,     123,    0.0,   -53,   0.0;
        1,   0,   0,   0,   1,    27.7,      63,    0.1,   -33,   0.0;
        0,   0,   0,   2,   0,    14.8,      63,    0.0,    -2,   0.0;
        -1,   0,   2,   2,   2,     9.6,     -59,    0.0,    26,   0.0;
        -1,   0,   0,   0,   1,   -27.4,     -58,   -0.1,    32,   0.0;
        1,   0,   2,   0,   1,     9.1,     -51,    0.0,    27,   0.0;
        -2,   0,   0,   2,   0,  -205.9,     -48,    0.0,     1,   0.0;
        -2,   0,   2,   0,   1,  1305.5,      46,    0.0,   -24,   0.0;
        0,   0,   2,   2,   2,     7.1,     -38,    0.0,    16,   0.0;
        2,   0,   2,   0,   2,     6.9,     -31,    0.0,    13,   0.0;
    2,   0,   0,   0,   0,    13.8,      29,    0.0,    -1,   0.0;
    1,   0,   2,  -2,   2,    23.9,      29,    0.0,   -12,   0.0;
    0,   0,   2,   0,   0,    13.6,      26,    0.0,    -1,   0.0;
    0,   0,   2,  -2,   0,   173.3,     -22,    0.0,     0,   0.0;
    -1,   0,   2,   0,   1,    27.0,      21,    0.0,   -10,   0.0;
    0,   2,   0,   0,   0,   182.6,      17,   -0.1,     0,   0.0;
    0,   2,   2,  -2,   2,    91.3,     -16,    0.1,     7,   0.0;
    -1,   0,   0,   2,   1,    32.0,      16,    0.0,    -8,   0.0;
    0,   1,   0,   0,   1,   386.0,     -15,    0.0,     9,   0.0;
    1,   0,   0,  -2,   1,   -31.7,     -13,    0.0,     7,   0.0;
    0,  -1,   0,   0,   1,  -346.6,     -12,    0.0,     6,   0.0;
    2,   0,  -2,   0,   0, -1095.2,      11,    0.0,     0,   0.0;
    -1,   0,   2,   2,   1,     9.5,     -10,    0.0,     5,   0.0;
    1,   0,   2,   2,   2,     5.6,      -8,    0.0,     3,   0.0;
    0,  -1,   2,   0,   2,    14.2,      -7,    0.0,     3,   0.0;
    0,   0,   2,   2,   1,     7.1,      -7,    0.0,     3,   0.0;
    1,   1,   0,  -2,   0,   -34.8,      -7,    0.0,     0,   0.0;
    0,   1,   2,   0,   2,    13.2,       7,    0.0,    -3,   0.0;
    -2,   0,   0,   2,   1,  -199.8,      -6,    0.0,     3,   0.0;
    0,   0,   0,   2,   1,    14.8,      -6,    0.0,     3,   0.0;
    2,   0,   2,  -2,   2,    12.8,       6,    0.0,    -3,   0.0;
    1,   0,   0,   2,   0,     9.6,       6,    0.0,     0,   0.0;
    1,   0,   2,  -2,   1,    23.9,       6,    0.0,    -3,   0.0;
    0,   0,   0,  -2,   1,   -14.7,      -5,    0.0,     3,   0.0;
    0,  -1,   2,  -2,   1,   346.6,      -5,    0.0,     3,   0.0;
    2,   0,   2,   0,   1,     6.9,      -5,    0.0,     3,   0.0;
    1,  -1,   0,   0,   0,    29.8,       5,    0.0,     0,   0.0;
    1,   0,   0,  -1,   0,   411.8,      -4,    0.0,     0,   0.0;
    0,   0,   0,   1,   0,    29.5,      -4,    0.0,     0,   0.0;
    0,   1,   0,  -2,   0,   -15.4,      -4,    0.0,     0,   0.0;
    1,   0,  -2,   0,   0,   -26.9,       4,    0.0,     0,   0.0;
    2,   0,   0,  -2,   1,   212.3,       4,    0.0,    -2,   0.0;
    0,   1,   2,  -2,   1,   119.6,       4,    0.0,    -2,   0.0;
    1,   1,   0,   0,   0,    25.6,      -3,    0.0,     0,   0.0;
    1,  -1,   0,  -1,   0, -3232.9,      -3,    0.0,     0,   0.0;
    -1,  -1,   2,   2,   2,     9.8,      -3,    0.0,     1,   0.0;
    0,  -1,   2,   2,   2,     7.2,      -3,    0.0,     1,   0.0;
    1,  -1,   2,   0,   2,     9.4,      -3,    0.0,     1,   0.0;
    3,   0,   2,   0,   2,     5.5,      -3,    0.0,     1,   0.0;
    -2,   0,   2,   0,   2,  1615.7,      -3,    0.0,     1,   0.0;
    1,   0,   2,   0,   0,     9.1,       3,    0.0,     0,   0.0;
    -1,   0,   2,   4,   2,     5.8,      -2,    0.0,     1,   0.0;
    1,   0,   0,   0,   2,    27.8,      -2,    0.0,     1,   0.0;
    -1,   0,   2,  -2,   1,   -32.6,      -2,    0.0,     1,   0.0;
    0,  -2,   2,  -2,   1,  6786.3,      -2,    0.0,     1,   0.0;
    -2,   0,   0,   0,   1,   -13.7,      -2,    0.0,     1,   0.0;
    2,   0,   0,   0,   1,    13.8,       2,    0.0,    -1,   0.0;
    3,   0,   0,   0,   0,     9.2,       2,    0.0,     0,   0.0;
    1,   1,   2,   0,   2,     8.9,       2,    0.0,    -1,   0.0;
    0,   0,   2,   1,   2,     9.3,       2,    0.0,    -1,   0.0;
    1,   0,   0,   2,   1,     9.6,      -1,    0.0,     0,   0.0;
    1,   0,   2,   2,   1,     5.6,      -1,    0.0,     1,   0.0;
    1,   1,   0,  -2,   1,   -34.7,      -1,    0.0,     0,   0.0;
    0,   1,   0,   2,   0,    14.2,      -1,    0.0,     0,   0.0;
    0,   1,   2,  -2,   0,   117.5,      -1,    0.0,     0,   0.0;
    0,   1,  -2,   2,   0,  -329.8,      -1,    0.0,     0,   0.0;
    1,   0,  -2,   2,   0,    23.8,      -1,    0.0,     0,   0.0;
    1,   0,  -2,  -2,   0,    -9.5,      -1,    0.0,     0,   0.0;
    1,   0,   2,  -2,   0,    32.8,      -1,    0.0,     0,   0.0;
    1,   0,   0,  -4,   0,   -10.1,      -1,    0.0,     0,   0.0;
    2,   0,   0,  -4,   0,   -15.9,      -1,    0.0,     0,   0.0;
    0,   0,   2,   4,   2,     4.8,      -1,    0.0,     0,   0.0;
    0,   0,   2,  -1,   2,    25.4,      -1,    0.0,     0,   0.0;
    -2,   0,   2,   4,   2,     7.3,      -1,    0.0,     1,   0.0;
    2,   0,   2,   2,   2,     4.7,      -1,    0.0,     0,   0.0;
    0,  -1,   2,   0,   1,    14.2,      -1,    0.0,     0,   0.0;
    0,   0,  -2,   0,   1,   -13.6,      -1,    0.0,     0,   0.0;
    0,   0,   4,  -2,   2,    12.7,       1,    0.0,     0,   0.0;
    0,   1,   0,   0,   2,   409.2,       1,    0.0,     0,   0.0;
    1,   1,   2,  -2,   2,    22.5,       1,    0.0,    -1,   0.0;
    3,   0,   2,  -2,   2,     8.7,       1,    0.0,     0,   0.0;
    -2,   0,   2,   2,   2,    14.6,       1,    0.0,    -1,   0.0;
    -1,   0,   0,   0,   2,   -27.3,       1,    0.0,    -1,   0.0;
    0,   0,  -2,   2,   1,  -169.0,       1,    0.0,     0,   0.0;
    0,   1,   2,   0,   1,    13.1,       1,    0.0,     0,   0.0;
    -1,   0,   4,   0,   2,     9.1,       1,    0.0,     0,   0.0;
    2,   1,   0,  -2,   0,   131.7,       1,    0.0,     0,   0.0;
    2,   0,   0,   2,   0,     7.1,       1,    0.0,     0,   0.0;
    2,   0,   2,  -2,   1,    12.8,       1,    0.0,    -1,   0.0;
    2,   0,  -2,   0,   1,  -943.2,       1,    0.0,     0,   0.0;
    1,  -1,   0,  -2,   0,   -29.3,       1,    0.0,     0,   0.0;
    -1,   0,   0,   1,   1,  -388.3,       1,    0.0,     0,   0.0;
    -1,  -1,   0,   2,   1,    35.0,       1,    0.0,     0,   0.0;
    0,   1,   0,   1,   0,    27.3,       1,    0.0,     0,   0.0];

dpsi=0;
deps=0;
for i=1:106
    ang=0;
    for j=1:5
        ang=ang+nut(i,j)*f(j);
    end
    dpsi=dpsi+(nut(i,7)+nut(i,8)*t)*sin(ang);
    deps=deps+(nut(i,9)+nut(i,10)*t)*cos(ang);
end
dpsi=dpsi*1e-4*AS2R;
deps=deps*1e-4*AS2R;

return

repo: GINav


text:
function lambda = lambda_measure ( dim_num, n, z )

%*****************************************************************************80
%
%% LAMBDA_MEASURE determines the pointset quality measure LAMBDA.
%
%  Discussion:
%
%    The LAMBDA measure of point distribution quality for a set Z of
%    N points in an DIM_NUM-dimensional region is defined as follows:
%
%    Let
%
%      GAMMA(I) = minimum ( 1 <= J <= N, I /= J ) distance ( Z(I), Z(J) )
%
%    and let
%
%      GAMMA_AVE = sum ( 1 <= I <= N ) GAMMA(I) / N
%
%    then
%
%      LAMBDA = sqrt ( sum ( 1 <= I <= N ) ( GAMMA(I) - GAMMA_AVE )**2 / N )
%        / GAMMA_AVE
%
%    An ideally regular mesh would have GAMMA(I) = GAMMA_AVE for all I,
%    so that LAMBDA would be 0.  Under this measure, the mesh with the
%    smaller value of LAMBDA is to be preferred.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2004
%
%  Author:
%
%    Max Gunzburger
%    John Burkardt
%
%  Reference:
%
%    Max Gunzburger and John Burkardt,
%    Uniformity Measures for Point Samples in Hypercubes.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real Z(DIM_NUM,N), the points.
%
%    Output, real LAMBDA, the LAMBDA quality measure.
%
%  Local parameters:
%
%    Local, real GAMMA_MAX, the maximum, over all points,
%    of the minimum distance to a distinct point.
%
%    Local, real GAMMA_MIN, the minimum, over all points,
%    of the minimum distance to a distinct point.
%

%
%  Take care of ridiculous cases.
%
  if ( n <= 1 )
    lambda = 0.0;
    return
  end
%
%  Compute the minimum spacing between distinct points of the set.
%
  gamma = pointset_spacing ( dim_num, n, z );
%
%  Average the minimum spacing.
%
  gamma_ave = sum ( gamma(1:n) ) / n;
%
%  Compute a weighted variance.
%
  if ( gamma_ave <= 0.0 )
    lambda = r8_huge ( );
  else
    lambda = sqrt ( sum ( ( gamma(1:n) - gamma_ave ).^2 ) / n ) / gamma_ave;
  end

  return
end


text:
classdef C10MOP8 < EvoXBenchProblem
% <multi> <real> <large/none> <expensive/none>

%------------------------------- Reference --------------------------------
% Z. Lu, R. Cheng, Y. Jin, K. C. Tan, and K. Deb, Neural architecture
% search as multiobjective optimization benchmarks: Problem formulation and
% performance assessment, IEEE Transactions on Evolutionary Computation,
% 2023.
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            config.name = 'darts';
            config.args.objs = 'err&params';
            config.args.normalized_objectives = false;
            obj.Setting@EvoXBenchProblem(config);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% MLTOOLSTOOLBOXES Load in the relevant toolboxes for the MLTOOLS.

importLatest('datasets');
importLatest('optimi');
importLatest('ndlutil');
importLatest('netlab');
importLatest('mvu');
importLatest('lmvu');
importLatest('lle');
importLatest('jdqr');
importLatest('isomap');
importLatest('sedumi');

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat