text:
function result=cosmo_dim_generalization_measure(ds,varargin)
% measure generalization across pairwise combinations over time (or any other dimension)
%
% result=cosmo_dim_generalization_measure(ds,varargin)
%
% Inputs:
%   ds                  dataset struct with d being a sample dimension, and
%                       with ds.sa.chunks==1 for samples to use for
%                       training and ds.sa.chunks==2 for those to use for
%                       testing. Other values for chunks are not allowed.
%   'measure',m         function handle to apply to combinations of samples
%                       in the input dataset ds, such as
%                       - @cosmo_correlation_measure
%                       - @cosmo_crossvalidation_measure
%                       - @cosmo_target_dsm_corr_measure
%   'dimension',d       dimension along which to generalize. Typically this
%                       will be 'time' for MEEG data
%   'radius',r          radius used for the d dimension. For example, when
%                       set to r=4 with d='time', then 4*2+1 time points
%                       are used to asses generalization, (except on the
%                       edges). Note that when using a radius>0, it is
%                       assumed that splits of the dataset by dimension d
%                       have corresponding elements in the same order
%                       (such as provided by cosmo_dim_transpose).
%   'nproc', np         Use np parallel threads. (Multiple threads may
%                       speed up computations). If parallel processing is
%                       not available, or if this option is not provided,
%                       then a single thread is used.
%   K,V                 any other key-value pairs necessary for the measure
%                       m, for example 'classifier' if
%                       m=@cosmo_crossvalidation_measure.
%
% Output:
%    result             dataset with ['train_' d] and ['test_' d] as sample
%                       dimensions, i.e. these are in ds.a.sdim.labels
%                       result.samples is Nx1, where N=K*J is the number of
%                       combinations of (1) the K points in ds with
%                       chunks==1 and different values in dimension d, and
%                       (2) the J points in ds with chunks==2 and different
%                       values in dimension d.
%
% Examples:
%     % Generalization over time
%     sz='big';
%     train_ds=cosmo_synthetic_dataset('type','timelock','size',sz,...
%                                              'nchunks',2,'seed',1);
%     test_ds=cosmo_synthetic_dataset('type','timelock','size',sz,...
%                                              'nchunks',3,'seed',2);
%     % set chunks
%     train_ds.sa.chunks(:)=1;
%     test_ds.sa.chunks(:)=2;
%     %
%     % construct the dataset
%     ds=cosmo_stack({train_ds, test_ds});
%     %
%     % make time a sample dimension
%     dim_label='time';
%     ds_time=cosmo_dim_transpose(ds,dim_label,1);
%     %
%     % set measure and its arguments
%     measure_args=struct();
%     %
%     % use correlation measure
%     measure_args.measure=@cosmo_correlation_measure;
%     % dimension of interest is 'time'
%     measure_args.dimension=dim_label;
%     %
%     % run time-by-time generalization analysis
%     dgm_ds=cosmo_dim_generalization_measure(ds_time,measure_args,...
%                                               'progress',false);
%     %
%     % the output has train_time and test_time as sample dimensions
%     cosmo_disp(dgm_ds.a)
%     %|| .sdim
%     %||   .labels
%     %||     { 'train_time'  'test_time' }
%     %||   .values
%     %||     { [  -0.2        [  -0.2
%     %||         -0.15          -0.15
%     %||          -0.1           -0.1
%     %||           :              :
%     %||             0              0
%     %||          0.05           0.05
%     %||           0.1 ]@7x1      0.1 ]@7x1 }
%
%
%     % Searchlight example
%     % (This example requires FieldTrip)
%     cosmo_skip_test_if_no_external('fieldtrip');
%     %
%     sz='big';
%     train_ds=cosmo_synthetic_dataset('type','timelock','size',sz,...
%                                              'nchunks',2,'seed',1);
%     test_ds=cosmo_synthetic_dataset('type','timelock','size',sz,...
%                                              'nchunks',3,'seed',2);
%     % set chunks
%     train_ds.sa.chunks(:)=1;
%     test_ds.sa.chunks(:)=2;
%     %
%     % construct the dataset
%     ds=cosmo_stack({train_ds, test_ds});
%     %
%     % make time a sample dimension
%     dim_label='time';
%     ds_time=cosmo_dim_transpose(ds,dim_label,1);
%     %
%     % set measure and its arguments
%     measure_args=struct();
%     %
%     % use correlation measure
%     measure_args.measure=@cosmo_correlation_measure;
%     % dimension of interest is 'time'
%     measure_args.dimension=dim_label;
%     %
%     % only to make this example run fast, most channels are eliminated
%     % (there is no other reason to do this step)
%     ds_time=cosmo_slice(ds_time,ds_time.fa.chan<=20,2);
%     ds_time=cosmo_dim_prune(ds_time);
%     %
%     % define neighborhood for channels
%     nbrhood=cosmo_meeg_chan_neighborhood(ds_time,...
%                                 'chantype','meg_combined_from_planar',...
%                                 'count',5,'label','dataset');
%     %
%     % run searchlight with generalization measure
%     measure=@cosmo_dim_generalization_measure;
%     dgm_sl_ds=cosmo_searchlight(ds_time,nbrhood,measure,measure_args,...
%                                                 'progress',false);
%     %
%     % the output has train_time and test_time as sample dimensions,
%     % and chan as feature dimension
%     cosmo_disp(dgm_sl_ds.a,'edgeitems',1)
%     %|| .fdim
%     %||   .labels
%     %||     { 'chan' }
%     %||   .values
%     %||     { { 'MEG0112+0113' ... 'MEG0712+0713'   }@1x7 }
%     %|| .meeg
%     %||   .samples_type
%     %||     'timelock'
%     %||   .samples_field
%     %||     'trial'
%     %||   .samples_label
%     %||     'rpt'
%     %|| .sdim
%     %||   .labels
%     %||     { 'train_time'  'test_time' }
%     %||   .values
%     %||     { [ -0.2        [ -0.2
%     %||           :             :
%     %||          0.1 ]@7x1     0.1 ]@7x1 }
%
%
% Notes:
%   - this function can be used together with searchlight
%   - to make a dimension d a sample dimension from a feature dimension
%     (usually necessary before running this function), or the other way
%     around (usually necessary after running this function), use
%     cosmo_dim_transpose.
%   - a 'partition' argument should not be provided, because this function
%     generates them itself. The partitions are generated so that there
%     is a single fold; samples with chunks==1 are always used for training
%     and those with chunks==2 are used for testing (e.g. when using
%     m=@cosmo_crossvalidation_measure). In the case of using
%     m=@cosmo_correlation_measure, this amounts to split-half
%     correlations.
%
% See also: cosmo_correlation_measure, cosmo_crossvalidation_measure
%           cosmo_target_dsm_corr_measure, cosmo_searchlight,
%           cosmo_dim_transpose
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    defaults=struct();
    defaults.radius=0;
    defaults.progress=1;
    defaults.check_partitions=false;
    defaults.nproc=1;
    opt=cosmo_structjoin(defaults,varargin);

    cosmo_check_dataset(ds);
    check_input(ds,opt);

    % get training and test set
    halves=split_dataset_in_train_and_test(ds);

    % split the data in two halves
    [train_values,train_splits]=split_half_by_dimension(halves{1},opt);
    [test_values,test_splits]=split_half_by_dimension(halves{2},opt);

    halves=[]; % let GC do its work

    % get number of processes available
    nproc_available=cosmo_parallel_get_nproc_available(opt);

    % Matlab needs newline character at progress message to show it in
    % parallel mode; Octave should not have newline character
    environment=cosmo_wtf('environment');
    progress_suffix=get_progress_suffix(environment);

    % split training data in multiple parts, so that each thread can do a
    % subset of all the work
    % set options for each worker process
    worker_opt_cell=cell(1,nproc_available);
    block_size=ceil(length(train_values)/nproc_available);
    first=1;
    for p=1:nproc_available
        last=min(first+block_size-1,length(train_values));
        block_idxs=first:last;

        worker_opt=struct();
        worker_opt.train_splits=train_splits(block_idxs);
        worker_opt.train_values=train_values(block_idxs);
        worker_opt.train_values_ori=train_values;
        worker_opt.train_values_idx=block_idxs;
        worker_opt.test_splits=test_splits;
        worker_opt.test_values=test_values;
        worker_opt.opt=opt;
        worker_opt.worker_id=p;
        worker_opt.nworkers=nproc_available;
        worker_opt.progress=opt.progress;
        worker_opt.progress_suffix=progress_suffix;
        worker_opt_cell{p}=worker_opt;
        first=last+1;
    end

    % Run process for each worker in parallel
    % Note that when using nproc=1, cosmo_parcellfun does actually not
    % use any parallellization; the result is a cell with a single element.
    result_map_cell=cosmo_parcellfun(opt.nproc,...
                                     @run_with_worker,...
                                    worker_opt_cell,...
                                    'UniformOutput',false);

    result=cosmo_stack(result_map_cell,1);
    cosmo_check_dataset(result);

function result=run_with_worker(worker_opt)
% run dimgen using the options in worker_opt

    train_splits = worker_opt.train_splits;
    train_values = worker_opt.train_values;
    train_values_ori = worker_opt.train_values_ori;
    train_values_idx = worker_opt.train_values_idx;
    test_splits = worker_opt.test_splits;
    test_values = worker_opt.test_values;
    opt = worker_opt.opt;
    worker_id=worker_opt.worker_id;
    nworkers=worker_opt.nworkers;
    progress=worker_opt.progress;
    progress_suffix=worker_opt.progress_suffix;

    % set partitions in case a crossvalidation or correlation measure is
    % used
    ntrain_elem=cellfun(@(x)size(x.samples,1),train_splits);
    ntest_elem=cellfun(@(x)size(x.samples,1),test_splits);
    opt.partitions=struct();
    opt.partitions.train_indices=cell(1);
    opt.partitions.test_indices=cell(1);

    % remove the dimension and measure arguments from the input
    dimension=opt.dimension;
    measure=opt.measure;

    opt=rmfield(opt,'dimension');
    opt=rmfield(opt,'measure');

    train_label=['train_' dimension];
    test_label=['test_' dimension];

    % see if progress has to be shown
    show_progress=~isempty(progress) && ...
                        progress && ...
                        worker_id==1;
    if show_progress
        prev_progress_msg='';
        clock_start=clock();
    end

    % allocate space for output
    ntrain=numel(train_values);
    ntest=numel(test_values);

    result_cell=cell(ntrain*ntest,1);

    % last non-empty row in result_cell
    pos=0;

    for k=1:ntrain
        % update partitions train set
        opt.partitions.train_indices{1}=1:ntrain_elem(k);
        for j=1:ntest
            % update partitions test set
            opt.partitions.test_indices{1}=ntrain_elem(k)+...
                                                (1:ntest_elem(j));
            % merge training and test dataset
            ds_merged=cosmo_stack({train_splits{k},test_splits{j}},...
                                                        1,1,false);

            opt.partitions=cosmo_balance_partitions(opt.partitions,...
                                            ds_merged,opt);

            % apply measure
            ds_result=measure(ds_merged,opt);

            % set dimension attributes
            nsamples=size(ds_result.samples,1);
            ds_result.sa.(train_label)=repmat(train_values_idx(k),...
                                            nsamples,1);
            ds_result.sa.(test_label)=repmat(j,nsamples,1);

            % store result
            pos=pos+1;
            result_cell{pos}=ds_result;
        end

        if show_progress
            if nworkers>1
                if k==ntrain
                    % other workers may be slower than first worker
                    msg=sprintf(['worker %d has completed; waiting for '...
                                    'other workers to finish...%s'],...
                                    worker_id, progress_suffix);
                else
                    % can only show progress from a single worker;
                    % therefore show progress of first worker
                    msg=sprintf('for worker %d / %d%s', worker_id, ...
                                    nworkers, progress_suffix);
                end
            else
                % no specific message
                msg='';
            end
            prev_progress_msg=cosmo_show_progress(clock_start, ...
                            k/ntrain, msg, prev_progress_msg);
        end

    end

    % merge results into a dataset
    result=cosmo_stack(result_cell,1,'drop_nonunique');
    if isfield(result,'sa') && isfield(result.sa,dimension)
        result.sa=rmfield(result.sa,dimension);
    end

    % set dimension attributes in the sample dimension
    result=add_sample_attr(result, {train_label;test_label},...
                                    {train_values_ori;test_values});

function check_input(ds,opt)
    % ensure input is kosher
    cosmo_isfield(opt,{'dimension','measure'},true);

    dimension=opt.dimension;
    dim_pos=cosmo_dim_find(ds,dimension,true);
    if dim_pos~=1
        error(['''%s'' must be a sample dimension (not a feature) '...
                    'dimension. To make ''%s'' a sample dimension in '...
                    'a dataset struct ds, use\n\n'...
                    '  cosmo_dim_transpose(ds,''%s'',1);'],...
                    dimension,dimension,dimension);
    end

    measure=opt.measure;
    if ~isa(measure,'function_handle')
        error('the ''measure'' argument must be a function handle');
    end

    if isfield(opt,'partitions')
        error(['the partitions argument is not allowed for this '...
                    'function, because it generates partitions itself.'...
                    'The dataset should have two chunks, with '...
                    'chunks set to 1 for the training set and '...
                    'set to 2 for the testing set']);
    end

function halves=split_dataset_in_train_and_test(ds)
    % return cell with {train_ds,test_ds}
    halves=cosmo_split(ds,'chunks',1);
    if numel(halves)~=2 || ...
            halves{1}.sa.chunks(1)~=1 || ...
            halves{2}.sa.chunks(1)~=2
        error(['chunks must be 1 (for the training set) or 2'...
                '(for the testing set)' ]);
    end

function [values,splits]=split_half_by_dimension(ds,opt)
    % split dataset by ds.a.(opt.dimension)
    dimension=opt.dimension;
    ds_pruned=cosmo_dim_prune(ds,'labels',{dimension},'dim',1);
    ds_tr=cosmo_dim_transpose(ds_pruned,dimension,2);

    nbrhood=cosmo_interval_neighborhood(ds_tr,dimension,opt);
    assert(isequal(nbrhood.a.fdim.labels,{dimension}));

    counts=cellfun(@numel,nbrhood.neighbors);

    keep_nbrs=find(counts==max(counts));

    % remove dimension information
    ds_tr=remove_fa_field(ds_tr,dimension);

    values=nbrhood.a.fdim.values{1}(keep_nbrs);
    sz=size(values);
    if sz(1)==1
        values=values';
    elseif sz(2)~=1
        error('dimension %s must be a row vector', dimension);
    end

    n=numel(keep_nbrs);
    splits=cell(n,1);

    for k=1:n
        idx=nbrhood.neighbors{keep_nbrs(k)};
        splits{k}=cosmo_slice(ds_tr,idx,2,false);
    end

function ds=add_sample_attr(ds, dim_labels, dim_values)
    if ~isfield(ds,'a') || ~isfield(ds.a,'sdim')
        ds.a.sdim=struct();
        ds.a.sdim.labels=cell(1,0);
        ds.a.sdim.values=cell(1,0);
    end

    ds.a.sdim.values=[ds.a.sdim.values(:); dim_values]';
    ds.a.sdim.labels=[ds.a.sdim.labels(:); dim_labels]';

function ds=remove_fa_field(ds,label)
    if isfield(ds.fa,label);
        ds.fa=rmfield(ds.fa,label);
    end

    [dim, index, attr_name, dim_name]=cosmo_dim_find(ds,...
                                                label,false);
    if ~isempty(dim)
        sfdim=ds.a.(dim_name);
        m=~cosmo_match(sfdim.labels,label);
        sfdim.values=sfdim.values(m);
        sfdim.labels=sfdim.labels(m);
        ds.a.(dim_name)=sfdim;
    end

function suffix=get_progress_suffix(environment)
    % Matlab needs newline character at progress message to show it in
    % parallel mode; Octave should not have newline character

    switch environment
        case 'matlab'
            suffix=sprintf('\n');
        case 'octave'
            suffix='';
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function y = gprnd(X, logtheta, covfunc)

if ~iscell(covfunc)
  covfunc = {covfunc};
end

K = regularize(feval(covfunc{:}, logtheta, X, X));% + 1e-5*diag(sparse(ones(cols(X),1)));

y = mymvnrnd(0, K);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function display_multifractal(S,t0)
	clrs = {'r','g','b','c','k','m'};

	S1 = [S{2}.signal{:}].';
	S2 = [S{3}.signal{:}].';

	subplot(121);
	plot(S1(:,t0));
	title('m=1');
	xlabel('j1');

	j1s = unique(S{3}.meta.scale(:,1)).';
	subplot(122);
	hold on
	for j1 = j1s
		plot(S2(S{3}.meta.scale(:,1)==j1,t0), ...
			clrs{mod(j1,length(clrs))+1});
	end
	hold off
	title('m=2');
	xlabel('j2-j1');
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function aedat = Reorientate(aedat, transpose, flipX, flipY)

%{
For polarity and frame data, if present, apply transpose, flipX and flipY
operations in that order according to the boolean inputs.
%}

dbstop if error

% ReorientatePolarity function depends on this step being carried out
% first
if transpose
    aedat.info.deviceAddressSpace = aedat.info.deviceAddressSpace([2 1]);
end

aedat = ReorientatePolarity(aedat, transpose, flipX, flipY);

aedat = ReorientateFrames(aedat, transpose, flipX, flipY);

author: panpanfei


text:
function x = mono_next_grlex ( m, x )

%*****************************************************************************80
%
%% MONO_NEXT_GRLEX: grlex next monomial.
%
%  Discussion:
%
%    Example:
%
%    M = 3
%
%    #  X(1)  X(2)  X(3)  Degree
%      +------------------------
%    1 |  0     0     0        0
%      |
%    2 |  0     0     1        1
%    3 |  0     1     0        1
%    4 |  1     0     0        1
%      |
%    5 |  0     0     2        2
%    6 |  0     1     1        2
%    7 |  0     2     0        2
%    8 |  1     0     1        2
%    9 |  1     1     0        2
%   10 |  2     0     0        2
%      |
%   11 |  0     0     3        3
%   12 |  0     1     2        3
%   13 |  0     2     1        3
%   14 |  0     3     0        3
%   15 |  1     0     2        3
%   16 |  1     1     1        3
%   17 |  1     2     0        3
%   18 |  2     0     1        3
%   19 |  2     1     0        3
%   20 |  3     0     0        3
%
%    Thanks to Stefan Klus for pointing out a discrepancy in a previous
%    version of this code, 05 February 2015.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 February 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer X(M), the current monomial.
%    The first item is X = [ 0, 0, ..., 0, 0 ].
%
%    Output, integer X(M), the next monomial.
%

%
%  Ensure that 1 <= M.
%
  if ( m < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MONO_NEXT_GRLEX - Fatal error!' );
    fprintf ( 1, '  M < 1\n' );
    error ( 'MONO_NEXT_GRLEX - Fatal error!' );
  end
%
%  Ensure that 0 <= XC(I).
%
  for i = 1 : m
    if ( x(i) < 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'MONO_NEXT_GRLEX - Fatal error!' );
      fprintf ( 1, '  X(I) < 0\n' );
      error ( 'MONO_NEXT_GRLEX - Fatal error!' );
    end
  end
%
%  Find I, the index of the rightmost nonzero entry of X.
%
  i = 0;
  for j = m : -1 : 1
    if ( 0 < x(j) )
      i = j;
      break
    end
  end
%
%  set T = X(I)
%  set X(I) to zero,
%  increase X(I-1) by 1,
%  increment X(M) by T-1.
%
  if ( i == 0 )
    x(m) = 1;
    return
  elseif ( i == 1 )
    t = x(1) + 1;
    im1 = m;
  elseif ( 1 < i )
    t = x(i);
    im1 = i - 1;
  end

  x(i) = 0;
  x(im1) = x(im1) + 1;
  x(m) = x(m) + t - 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function c4_normal_01_test ( )

%*****************************************************************************80
%
%% C4_NORMAL_01_TEST tests C4_NORMAL_01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 July 2006
%
%  Author:
%
%    John Burkardt
%
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C4_NORMAL_01_TEST\n' );
  fprintf ( 1, '  C4_NORMAL_01 computes pseudorandom complex values\n' );
  fprintf ( 1, '  normally distributed in the unit circle.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The initial seed is %d\n', seed );
  fprintf ( 1, '\n' );

  for i = 1 : 10
    [ r, seed ] = c4_normal_01 ( seed );
    fprintf ( 1, '  %6d  (%14f,  %14f)\n', i, real ( r ), imag ( r ) );
  end

  return
end


text:
function h = histp(x, xmin, xmax, nbins)
%HISTP	Histogram estimate of 1-dimensional probability distribution.
%
%	Description
%
%	HISTP(X, XMIN, XMAX, NBINS) takes a column vector X  of data values
%	and generates a normalized histogram plot of the  distribution. The
%	histogram has NBINS bins lying in the range XMIN to XMAX.
%
%	H = HISTP(...) returns a vector of patch handles.
%
%	See also
%	DEMGAUSS
%

%	Copyright (c) Ian T Nabney (1996-2001)

ndata = length(x);

bins = linspace(xmin, xmax, nbins);

binwidth = (xmax - xmin)/nbins;

num = hist(x, bins);

num = num/(ndata*binwidth);

h = bar(bins, num, 0.6);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function C = permutMat_yaw(tau_vec)
%permutMat_yaw takes tau_vec as an input and produces permutation matrix for
%yaw. Only initial conditions (yaw, yawdot, yaw2dot) and final (yawdot,
%yaw2dot) is being specified. Rest will be optimized with equality
%constraints.
r = 2; % up to 2nd derivative
n = 2*r+1; % 5th order
K = length(tau_vec);
C = eye((n+1)*K);
c = [];
idx = [];
for i=1:K-1
    if i == K-1
        ukn_row = (n+1)*i+1:(n+1)*i+4;  % Last term includes last angle
    else
        ukn_row = (n+1)*i+1:(n+1)*i+3;  % unknown rows: 7,8,9 , 13,14,15, ...
    end
    c = [c;C(ukn_row,:)];
    idx = [idx, ukn_row];
end
C(idx,:) = [];
C = [C;c];

author: yorgoon

repo: minimum-snap-geometric-control


text:
function [rad] = deg2rad(deg)
% Convert angle from degrees to radians.
rad = deg*pi/180;

author: Sable

repo: mcbench-benchmarks


text:
%% CONTRAST_BATCH_LOCAL uses the BATCH command to run the CONTRAST code locally.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 April 2010
%
%  Author:
%
%    John Burkardt
%
  clear

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CONTRAST_BATCH_LOCAL\n' );
  fprintf ( 1, '  Run CONTRAST_SCRIPT locally.\n' );
%
%  Read a black and white image and filter it locally.
%
  x = imread ( 'surfsup.tif' );
  xf = nlfilter ( x, [3,3], @contrast_enhance );
  xf = uint8 ( xf );
%
%  Set the number of workers.
%
  n = 4;
%
%  BATCH sends the script for execution.
%
  job = batch ( 'contrast_script', ...
    'Configuration', 'local', ...
    'CaptureDiary', true, ...
    'FileDependencies', { 'contrast_fun', 'contrast_enhance', 'surfsup.tif' }, ...
    'matlabpool', n );
%
%  WAIT pauses the MATLAB session til the job completes.
%
  wait ( job );
%
%  DIARY displays any messages printed during execution.
%
  diary ( job );
%
%  LOAD makes the script's workspace available.
%
%  y = output #1.
%
  load ( job );
%
%  DESTROY cleans up data about the job we no longer need.
%
  destroy ( job );
%
%  Display the original, single filtered, and SPMD filtered versions.
%
  figure ( 1 );

  subplot ( 1, 3, 1 );
  imshow ( x );
  title ( 'Original image' );

  subplot ( 1, 3, 2 );
  imshow( xf );
  title( 'Filtered on Client' );

  subplot ( 1, 3, 3 );
  imshow ( y );
  title ( 'Filtered on 4 SPMD Workers' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CONTRAST_BATCH_LOCAL\n' );
  fprintf ( 1, '  Normal end of execution.\n' );


text:
function rot = byEuler(varargin)
% define rotentations by Euler angles
%
% Syntax
%   rot = rotation.byEuler(phi1,Phi,phi2)          % Bunge convention
%   rot = rotation.byEuler(alpha,beta,gamma,'ZYZ') % Matthies convention
%   rot = rotation.byEuler(phi1,Phi,phi2,'Kocks')  % Kocks convention
%
% Input
%  phi1, Phi, phi2 - Euler angles in radiant
%
% Output
%  rot - @rotation
%
% Flags
%  Bunge, ZXZ -
%  ABG, Matthies, ZYZ -
%  Roe -
%  Kocks -
%  Canova -
%
% See also
% rotentation/rotentation rotentation/byMiller rotentation/byAxisAngle
% rotation/map

q = euler2quat(varargin{:});

rot = rotation(q);

author: mtex-toolbox

repo: mtex

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


text:
% STD_PLOTTF - plot ERSP/ITC images a component
%              or channel cluster in a STUDY. Also allows plotting scalp
%              maps.
% Usage:
%          >> std_plottf( times, freqs, data, 'key', 'val', ...)
% Inputs:
%  times - [vector] latencies in ms of the data points.
%  freqs - [vector] frequencies in Hz of the data points.
%  data  -  [cell array] mean data for each subject group and/or data
%           condition. For example, to plot mean ERPs from a STUDY
%           for epochs of 800 frames in two conditions from three groups
%           of 12 subjects:
%
%           >> data = { [800x12] [800x12] [800x12];... % 3 groups, cond 1
%                       [800x12] [800x12] [800x12] };  % 3 groups, cond 2
%           >> std_plottf(erp_ms,data);
%
%           By default, parametric statistics are computed across subjects
%           in the three groups. (group,condition) ERP averages are plotted.
%           See below and >> help statcond
%           for more information about the statistical computations.
%
% Optional display parameters:
%  'datatype'    - ['ersp'|'itc'] data type {default: 'ersp'}
%  'titles'      - [cell array of string] titles for each of the subplots.
%                  { default: none}
%
% Statistics options:
%  'groupstats'  - ['on'|'off'] Compute (or not) statistics across groups.
%                  {default: 'off'}
%  'condstats'   - ['on'|'off'] Compute (or not) statistics across groups.
%                  {default: 'off'}

%  'threshold'   - [NaN|real<<1] Significance threshold. NaN -> plot the
%                  p-values themselves on a different figure. When possible,
%                  significance regions are indicated below the data.
%                  {default: NaN}
%  'maskdata'    - ['on'|'off'] when threshold is non-NaN and not both
%                  condition and group statistics are computed, the user
%                  has the option to mask the data for significance.
%                  {default: 'off'}
%
% Other plotting options:
%  'plotmode'    - ['normal'|'condensed'] statistics plotting mode:
%                  'condensed' -> plot statistics under the curves
%                  (when possible); 'normal' -> plot them in separate
%                  axes {default: 'normal'}
%  'freqscale'   - ['log'|'linear'|'auto'] frequency plotting scale. This
%                  will only change the ordinate not interpolate the data.
%                  If you change this option blindly, your frequency scale
%                  might be inaccurate {default: 'auto'}
%
% ITC/ERSP image plotting options:
%  'tftopoopt'   - [cell array] TFTOPO plotting options (ERSP and ITC)
%  'caxis'       - [min max] color axis (ERSP, ITC, scalp maps)
%
% Scalp map plotting options:
%  'chanlocs'    - [struct] channel location structure
%
% Author: Arnaud Delorme, CERCO, CNRS, 2006-
%
% See also: POP_ERSPPARAMS, POP_ERPPARAMS, POP_SPECPARAMS, STATCOND

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

function [pgroup, pcond, pinter] = std_plottf(timevals, freqs, data, varargin)

pgroup = [];
pcond  = [];
pinter = [];
if nargin < 2
    help std_plottf;
    return;
end

opt = finputcheck( varargin, { 'titles'         'cell'   []              cellfun(@num2str, cell(20,20), 'uniformoutput', false);
                               'caxis'          'real'   []              [];
                               'ersplim'        'real'   []              []; % same as above
                               'itclim'         'real'   []              []; % same as above
                               'ylim'           'real'   []              [];
                               'tftopoopt'      'cell'   []              {};
                               'threshold'      'real'   []              NaN;
                               'unitx'          'string' []              'ms'; % just for titles
                               'unitcolor'      'string' {}              'dB';
                               'chanlocs'       'struct' []              struct('labels', {});
                               'freqscale'      'string' { 'log','linear','auto' }  'auto'; % note that paramsersp in std_erspplot contains the information as well
                               'effect'         'string' { 'main','marginal' }   'marginal';
                               'averagemode'    'string' { 'rms','ave' }   'rms';
                               'events'         'cell'   []              {};
                               'groupstats'     'cell'   []              {};
                               'condstats'      'cell'   []              {};
                               'interstats'     'cell'   []              {};
                               'maskdata'       'string' { 'on','off' }   'off';
                               'plottopo'       'string' { 'on','off' }   'off';
                               'datatype'       'string' { 'ersp','itc' 'erpim' }    'ersp';
                               'plotmode'       'string' { 'normal','condensed' }  'normal' }, 'std_plottf');
if ischar(opt), error(opt); end
if all(all(cellfun('size', data, 3)==1))               opt.singlesubject = 'on'; end

% remove empty entries
datapresent = ~cellfun(@isempty, data);
for c = size(data,1):-1:1, if sum(datapresent(c,:)) == 0, data(c,:) = []; opt.titles(c,:) = []; if ~isempty(opt.groupstats), opt.groupstats(c) = []; end; end; end
for g = size(data,2):-1:1, if sum(datapresent(:,g)) == 0, data(:,g) = []; opt.titles(:,g) = []; if ~isempty(opt.condstats ), opt.condstats( g) = []; end; end; end

if ~isempty(opt.groupstats) && ~isempty(opt.condstats) && strcmpi(opt.maskdata, 'on')
    disp('Cannot use ''maskdata'' option with both condition stat. and group stat. on');
    disp('Disabling statistics');
    opt.groupstats = {}; opt.condstats = {}; opt.maskdata = 'off';
end
if ~isempty(opt.ersplim), opt.caxis = opt.ersplim; end
if ~isempty(opt.itclim), opt.caxis = opt.itclim; end
onecol  = { 'b' 'b' 'b' 'b' 'b' 'b' 'b' 'b' 'b' 'b' };
manycol = { 'b' 'r' 'g' 'k' 'c' 'y' };

nc = size(data,1);
ng = size(data,2);
if nc >= ng, opt.transpose = 'on';
else         opt.transpose = 'off';
end

% test log frequencies
% --------------------
if length(freqs) > 2 && strcmpi(opt.freqscale, 'auto')
    midind  = floor(length(freqs)/2);
    if abs(freqs(midind)/freqs(end) - 1/2) < 0.1, opt.freqscale = 'linear';
    else                                          opt.freqscale = 'log';
    end
end

% condensed plot
% --------------
if strcmpi(opt.plotmode, 'condensed')
    meanplot = zeros(size(data{1},1), size(data{1},2));
    count = 0;
    for c = 1:nc
        for g = 1:ng
            if ~isempty(data{c,g})
                meanplot = meanplot + mean(data{c,g},3);
                count = count+1;
            end
        end
    end
    meanplot = meanplot/count;
    options = { 'chanlocs', opt.chanlocs, 'electrodes', 'off', 'cbar', 'on', ...
            'cmode', 'separate', opt.tftopoopt{:} };
    if strcmpi(opt.datatype, 'erpim'), options = { options{:} 'ylabel' 'Trials' }; end
    if strcmpi(opt.freqscale, 'log'), options = { options{:} 'logfreq', 'native' }; end
    tftopo( meanplot', timevals, freqs, 'title', opt.titles{1}, options{:});
    currentHangle = gca;
    if ~isempty( opt.caxis )
        caxis( currentHangle, opt.caxis )
    end
    cbar_standard(opt.datatype, ng, opt.unitcolor);
    axes(currentHangle);
    return;
end

% plotting parameters
% ------------------
if ng > 1 && ~isempty(opt.groupstats), addc = 1; else addc = 0; end
if nc > 1 && ~isempty(opt.condstats  ), addr = 1; else addr = 0; end

% compute significance mask
% --------------------------
pinterplot = {};
if strcmpi(opt.effect, 'marginal') || ng == 1 || nc == 1
    if ~isnan(opt.threshold) && ( ~isempty(opt.groupstats) || ~isempty(opt.condstats) )
        pcondplot  = opt.condstats;
        pgroupplot = opt.groupstats;
        maxplot = 1;
    else
        for ind = 1:length(opt.condstats),  pcondplot{ind}  = -log10(opt.condstats{ind}); end
        for ind = 1:length(opt.groupstats), pgroupplot{ind} = -log10(opt.groupstats{ind}); end
        maxplot = 3;
    end
elseif strcmpi(opt.effect, 'main') && ~isempty(opt.interstats)
    if ~isnan(opt.threshold) && ( ~isempty(opt.groupstats) || ~isempty(opt.condstats) )
        pcondplot  = { opt.interstats{2} };
        pgroupplot = { opt.interstats{1} };
        pinterplot = opt.interstats{3};
        maxplot = 1;
    else
        if ~isempty(opt.interstats{2}), pcondplot  = { -log10(opt.interstats{2}) }; end
        if ~isempty(opt.interstats{1}), pgroupplot = { -log10(opt.interstats{1}) }; end
        if ~isempty(opt.interstats{3}), pinterplot = -log10(opt.interstats{3}); end
        maxplot = 3;
    end
end

% -------------------------------
% masking for significance of not
% -------------------------------
statmask = 0;
if strcmpi(opt.maskdata, 'on') && ~isnan(opt.threshold) && ...
        (~isempty(opt.condstats) || ~isempty(opt.condstats))
    addc = 0; addr = 0; statmask = 1;
end

% -------------------------
% plot time/frequency image
% -------------------------
options = { 'chanlocs', opt.chanlocs, 'electrodes', 'off', 'cbar', 'off', ...
            'cmode', 'separate', opt.tftopoopt{:} };
if strcmpi(opt.freqscale, 'log'), options = { options{:} 'logfreq', 'native' }; end
if strcmpi(opt.datatype, 'erpim'), options = { options{:} 'ylabel' 'Trials' }; end

% adjust figure size
% ------------------
fig = figure('color', 'w');
pos = get(fig, 'position');
set(fig, 'position', [ pos(1)+15 pos(2)+15 pos(3)/2.5*(nc+addr), pos(4)/2*(ng+addc) ]);
pos = get(fig, 'position');
if strcmpi(opt.transpose, 'off'), set(gcf, 'position', [ pos(1) pos(2) pos(4) pos(3)]);
else                              set(gcf, 'position', pos);
end

% options
% -------
options = { 'limits' [NaN NaN NaN NaN opt.caxis] 'verbose' 'off' 'mode' opt.averagemode options{:} };

for c = 1:nc
    for g = 1:ng
        %hdl(c,g) = mysubplot(nc+addr, ng+addc, g + (c-1)*(ng+addc), opt.transpose);
        hdl(c,g) = mysubplot(nc+addr, ng+addc, c, g, opt.transpose);
        if ~isempty(data{c,g})
            if strcmpi(opt.plottopo, 'off')
                tmpplot = mean(data{c,g},3);
            else
                tmpplot = data{c,g};
                tmpplot = permute(tmpplot, [3 1 2 4]);
            end
            if ~isreal(tmpplot(1)), tmpplot = abs(tmpplot); end % comes second for processing single trials
            if statmask,
                if ~isempty(opt.condstats), tmpplot(find(pcondplot{g}(:) == 0)) = 0;
                else                        tmpplot(find(pgroupplot{c}(:) == 0)) = 0;
                end
            end
            if ~isempty(opt.events) && ~isempty(opt.events{c,g})
                 tmpevents = mean(opt.events{c,g},2);
            else tmpevents = [];
            end
            if strcmpi(opt.plottopo, 'on') && length(opt.chanlocs) > 1
                metaplottopo(tmpplot, 'chanlocs', opt.chanlocs, 'plotfunc', 'tftopo', 'squeeze', 'on', ...
                    'plotargs', { timevals, freqs, 'events', tmpevents, options{:} }, 'title', opt.titles{c,g});
            else
                tftopo( tmpplot, timevals, freqs, 'events', tmpevents, 'title', opt.titles{c,g}, options{:});
            end

            if c > 1
                ylabel('');
            end
        end

        % statistics across groups
        % -------------------------
        if strcmpi(opt.effect, 'marginal') || (strcmpi(opt.effect, 'main') && c == 1)
            if g == ng && ng > 1 && ~isempty(opt.groupstats) && ~isinf(pgroupplot{c}(1)) && ~statmask
                if strcmpi(opt.effect, 'main') && nc>1, centerc = nc/2-0.5; else centerc = 0; end
                hdl(c,g+1) = mysubplot(nc+addr, ng+addc, c+centerc, ng + 1, opt.transpose);
                pgroupplot{c}(pgroupplot{c}<0) = 0;
                tmpOptions = { 'limits' [nan nan nan nan -maxplot maxplot] options{3:end} };
                if strcmpi(opt.plottopo, 'on') && length(opt.chanlocs) > 1
                    metaplottopo(permute(pgroupplot{c}, [3 1 2]), 'chanlocs', opt.chanlocs, 'plotfunc', 'tftopo', 'squeeze', 'on', ...
                        'plotargs', { timevals, freqs, tmpOptions{:} }, 'title', opt.titles{c,g+1});
                else
                    tftopo( pgroupplot{c}, timevals, freqs, 'title', opt.titles{c,g+1}, tmpOptions{:});
                end
            end
        end
    end
end

for g = 1:ng
    % statistics across conditions
    % -----------------------------
    if strcmpi(opt.effect, 'marginal') || (strcmpi(opt.effect, 'main') && g == 1)
        if ~isempty(opt.condstats) && ~isinf(pcondplot{g}(1)) && ~statmask && nc > 1
            if strcmpi(opt.effect, 'main') && ng>1, centerg = ng/2-0.5; else centerg = 0; end
            hdl(nc+1,g) = mysubplot(nc+addr, ng+addc, nc+addr, g+centerg, opt.transpose);
            pcondplot{g}(pcondplot{g}<0) = 0;
            tmpOptions = { 'limits' [nan nan nan nan -maxplot maxplot] options{3:end} };
            if strcmpi(opt.plottopo, 'on') && length(opt.chanlocs) > 1
                metaplottopo(permute(pcondplot{g}, [3 1 2]), 'chanlocs', opt.chanlocs, 'plotfunc', 'tftopo', 'squeeze', 'on', ...
                    'plotargs', { timevals, freqs, tmpOptions{:} }, 'title', opt.titles{nc+1,g});
            else
                tftopo( pcondplot{g}, timevals, freqs, 'title', opt.titles{nc+1,g}, options{:});
            end
        end
    end
end

% statistics across group and conditions
% ---------------------------------------
if ~isempty(opt.groupstats) && ~isempty(opt.condstats) && ng > 1 && nc > 1 && ~isempty(pinterplot)
    hdl(nc+1,ng+1) = mysubplot(nc+addr, ng+addc, nc+addr, ng+1, opt.transpose);
    pinterplot(pinterplot<0) = 0;
    tftopo( pinterplot,  timevals, freqs, 'title', opt.titles{nc+1,ng+1}, options{:});
    caxis([-maxplot maxplot]);
    ylabel('');
end

% color bars
% ----------
axes(hdl(nc,ng));
cbar_standard(opt.datatype, ng, opt.unitcolor);
if isnan(opt.threshold) && (nc ~= size(hdl,1) || ng ~= size(hdl,2))
    ind = find(ishandle(hdl(end:-1:1)));
    axes(hdl(end-ind(1)+1));
    cbar_signif(ng, maxplot);
end

% mysubplot2 (allow to transpose if necessary)
% -------------------------------------------
function hdl = mysubplot(nr,nc,r,c,subplottype)

    cmargin = 0.2/nc;
    rmargin = 0.2/nr;
    if strcmpi(subplottype, 'transpose') || strcmpi(subplottype, 'on'),   hdl = subplot('position',[(r-1)/nr+rmargin (nc-c)/nc+cmargin 1/nr-2*rmargin 1/nc-2*cmargin]);
    elseif strcmpi(subplottype, 'normal') || strcmpi(subplottype, 'off'), hdl = subplot('position',[(c-1)/nc+cmargin (nr-r)/nr+rmargin 1/nc-2*cmargin 1/nr-2*rmargin]);
    elseif strcmpi(subplottype, 'noplot'), hdl = gca;
    else error('Unknown subplot type');
    end

% % mysubplot (allow to transpose if necessary)
% % -------------------------------------------
% function hdl = mysubplot(nr,nc,ind,transp);
%
%     r = ceil(ind/nc);
%     c = ind -(r-1)*nc;
%     if strcmpi(transp, 'on'), hdl = subplot(nc,nr,(c-1)*nr+r);
%     else                      hdl = subplot(nr,nc,(r-1)*nc+c);
%     end

% colorbar for ERSP and scalp plot
% --------------------------------
function cbar_standard(datatype, ng, unitcolor);
    pos = get(gca, 'position');
    tmpc = caxis;
    fact = fastif(ng == 1, 40, 20);
    tmp = axes('position', [ pos(1)+pos(3)+max(pos(3)/fact,0.006) pos(2) max(pos(3)/fact,0.01) pos(4) ]);
    set(gca, 'unit', 'normalized');
    if strcmpi(datatype, 'itc')
         cbar(tmp, 0, tmpc, 10); ylim([0.5 1]);
         title('ITC','fontsize',10,'fontweight','normal','interpreter','none');
    elseif strcmpi(datatype, 'erpim')
        cbar(tmp, 0, tmpc, 5);
    else
        cbar(tmp, 0, tmpc, 5);
        title(unitcolor);
    end

% colorbar for significance
% -------------------------
function cbar_signif(ng, maxplot);
    % Retrieving Defaults
    icadefs;

    pos = get(gca, 'position');
    tmpc = caxis;
    fact = fastif(ng == 1, 40, 20);
    tmp = axes('position', [ pos(1)+pos(3)+max(pos(3)/fact,0.006) pos(2) max(pos(3)/fact,0.01) pos(4) ]);
    map = colormap(DEFAULT_COLORMAP);
    n = size(map,1);
    cols = [ceil(n/2):n]';
    image([0 1],linspace(0,maxplot,length(cols)),[cols cols]);
    %cbar(tmp, 0, tmpc, 5);
    tick = linspace(0, maxplot, maxplot+1);
    set(gca, 'ytickmode', 'manual', 'YAxisLocation', 'right', 'xtick', [], ...
        'ytick', tick, 'yticklabel', round(10.^-tick*1000)/1000);
    xlabel('');
    colormap(DEFAULT_COLORMAP);

% rapid filtering for ERP
% -----------------------
function tmpdata2 = myfilt(tmpdata, lowpass, highpass, factor, filtertype)

    tmpdata2 = reshape(tmpdata, size(tmpdata,1), size(tmpdata,2)*size(tmpdata,3)*size(tmpdata,4));
    tmpdata2 = eegfiltfft(tmpdata2',lowpass, highpass, factor, filtertype)';
    tmpdata2 = reshape(tmpdata2, size(tmpdata,1), size(tmpdata,2), size(tmpdata,3), size(tmpdata,4));

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function f = p08_f ( m, n, x )

%*****************************************************************************80
%
%% P08_F evaluates the objective function for problem 08.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Marcin Molga, Czeslaw Smutnicki,
%    Test functions for optimization needs.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of arguments.
%
%    Input, real X(M,N), the arguments.
%
%    Output, real F(N), the function evaluated at the arguments.
%
  f = zeros ( n, 1 );

  y = r8vec_indicator ( m );
  y(1:m,1) = y(1:m,1) + 1.0;

  for j = 1 : n
    f(j) = sum ( abs ( x(1:m,j) ).^y(1:m) );
  end

  return
end


text:
function [ravenPath, prevDir] = findRAVENroot()
% findRAVENroot
%   Finds the root of the RAVEN directory, by searching for the path to
%   RAVEN2.png. Can also record the current directory, in case a function will
%   use the ravenPath to navigate to a precise folder, and it should return to
%   the previous directory afterwards. See e.g. optimizeProb calling glpk.

ST=dbstack('-completenames');
prevDir = pwd();
ravenPath = ST(strcmp({ST.name},'findRAVENroot')).file;
rootFound = 0;
while rootFound == 0
    isRoot = isfile(fullfile(ravenPath,'RAVEN2.png'));
    if isRoot
        rootFound = 1;
    else
        ravenPathOld = ravenPath;
        ravenPath = fileparts(ravenPath);
        if strcmp(ravenPathOld,ravenPath)
            error('Cannot find the RAVEN root directory. Make sure you have not removed the RAVEN2.png file from your RAVEN installation.')
        end
    end
end

author: SysBioChalmers

repo: RAVEN


text:
function [ft] = dm2ft(dm)
% Convert length from decimeters to feet.
% Chad A. Greene 2012
ft = dm*0.3280839895013;


text:
function [rc,fval,it] = ARfil(im,xm,rcinit,lagmax,p)

nobs = length(xm);

opties = optimset('Display','off','TolX',.001/sqrt(nobs),'TolFun',.0001);

[rc_tan,fval,exitflag,output]= fminunc('ARMA_MLfit',tan(.5*pi*rcinit),opties,im,xm,lagmax,p);
rc = 2/pi*atan(rc_tan);
it = output.iterations;


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ldbparam.m
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ldbparam.m
%
%% This function evaluate several room acoustic parameters, with the signal
%% processing done by Lundeby method.
%%
%% [s]=ldbparam(IR,fs,flag)
%%
%% The input is a room impulse response and its sampling frequency. The
%% flag variable indicates if the function should generate (1) or not (0)
%% the Schroeder decay plots.
%% The output it a text file with the value of several parameter for each
%% frequency band. If desired, returns a matrix with this values, where the
%% first line contains the frequency band central frequencies, instead of
%% text.

function [saida]=ldbparam(IR,fs,flag)

banda = filtros(IR,fs);
t = size(banda,2);
for n = 1:t
    s(1,n) = ceil(1000*2^(n-5));
    comeco = inicio(banda(:,n));
    fim = lundeby(banda(comeco:end,n),fs,flag);
    title(['Ponto de Cruzamento - Banda ',num2str(s(1,n))])
    if n == t-2
        title('Ponto de Cruzamento - Compensacao A ')
    elseif n == t-1
        title('Ponto de Cruzamento - Compensacao C ')
    elseif n == t
        title('Ponto de Cruzamento - Linear ')
    end

    aux = banda(comeco:fim,n).^2;
    [s(2,n),s(3,n),s(4,n),s(5,n),s(6,n)] = energeticos(aux,fs);
    [s(7,n),s(8,n),s(9,n),s(10,n)] = reverberacao(aux,fs,flag);
    title(['Curva de Decaimento - Banda ',num2str(s(1,n))])
    if n == t-2
        title('Curva de Decaimento - Compensacao A ')
    elseif n == t-1
        title('Curva de Decaimento - Compensacao C ')
    elseif n == t
        title('Curva de Decaimento - Linear ')
    end

end

if nargout == 1
    saida = s;
    saida(1,t-2) = (['A']);
    saida(1,t-1) = (['C']);
    saida(1,t) = (['L']);
else
    tabela(s,size(banda,2))
end


text:
function model = kpca_train(X,varargin)
% DESCRIPTION
% Kernel principal component analysis (KPCA)
%
%       model = kpca_train(X,varargin)
%
% INPUT
%   X            Training samples (N*d)
%                N: number of samples
%                d: number of features
%
% OUTPUT
%   model        KPCA model
%
%
% Created on 18th April 2019, by Kepeng Qiu.
%-------------------------------------------------------------%

% Default Parameters setting
options.sigma = 10;   % kernel width
options.dims  = 2;    % output dimension (dimensionality reduction)
                      % 2D data is easy to visualize
options.type  = 0;    % 0: dimensionality reduction or feature extraction
                      % 1: fault detection using KPCA
                      % 2: fault detection using dynamic KPCA(DKPCA)
options.beta  = 0.9;  % corresponding probabilities (for fault detection)
                      %
options.pcr   = 0.65; % principal contribution rate (for fault detection)
options.fd    = 0;    % 0: No fault diagnosis (fd)
                      % 1: fault diagnosis (fd)
options.lag   = 4;    % time lag (for DKPCA)
%

if rem(nargin-1,2)
    error('Parameters to kpca_train should be pairs')
end
numParameters = (nargin-1)/2;

for n =1:numParameters
    Parameters = varargin{(n-1)*2+1};
    value	= varargin{(n-1)*2+2};
    switch Parameters
            %
        case 'type'
            options.type = value;
            %
        case 'dims'
            options.dims = value;
            %
        case 'sigma'
            options.sigma = value;
            %
        case 'fd'
            options.fd = value;
            %
        case 'lag'
            options.lag = value;
            %
        case 'pcr'
            options.pcr = value;
            %
        case 'beta'
            options.beta = value;
    end
end

% number of training samples
L = size(X,1);

% DPCA
if options.type == 2
    %  Construct the augmented matrix
    X = constructAM(X,options.lag);
    L = size(X,1);
    model.lag = options.lag; % time lag
end

% Compute the kernel matrix
K = computeKM(X,X,options.sigma);

% Centralize the kernel matrix
unit = ones(L,L)/L;
K_c = K-unit*K-K*unit+unit*K*unit;

% Solve the eigenvalue problem
[V,D] = eigs(K_c/L);
% [V,D] = eig(K_c/L);
lambda = diag(D);

% Normalize the eigenvalue
V_s = V ./ sqrt(L*lambda)';
%  Lower version of MATLAB may report an error at line 70,
%  please replace it with the following code.
% ------------------------------
% for i = 1:size(lambda,1)
%     V_s(:,i) = V(:,i)/sqrt(L*lambda(i,1));
% end
% ------------------------------

% Compute the numbers of principal component
if options.type  == 1 || options.type  == 2  % fault detection
    dims = find(cumsum(lambda/sum(lambda)) >= ...
        options.pcr,1, 'first');
else % dimensionality reduction or feature extraction
    dims = options.dims;
end

% Extract the nonlinear component
mappedX  = K_c* V_s(:,1:dims) ;

% Store the results
model.mappedX =  mappedX ;
model.V_s = V_s;
model.lambda = lambda;
model.K_c = K_c;
model.L = L;
model.dims = dims;
model.X = X;
model.K = K;
model.unit = unit;
model.sigma = options.sigma;
model.diagnosis = options.fd;
model.type = options.type;

% Compute the threshold for fault detection
if options.type  == 1 || options.type  == 2
    model.beta = options.beta; % corresponding probabilities
    [SPE_limit,T2_limit] = comtupeLimit(model);
    model.SPE_limit = SPE_limit;
    model.T2_limit = T2_limit;
end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function ADimInfo = get_affine_matrix(this)
% Computes an affine transformation matrix based on the resolution and the
% first sampling point.
%
%   Y = MrDimInfo()
%   ADimInfo = Y.get_affine_matrix()
%
%   Note: The brackets after get_affine_matrix are necessary to prevent an
%   error in numArgumentsFromSubscript used for subsasign.
%
% This is a method of class MrDimInfo.
%
% DETAILS
%   - translation is the position of the first sampling point
%   - zoom is resolution or, if nan, samplingWidths
%
% IN
%
% OUT
%       Affine transformation matrix that describes the scaling
%       (resolution) and translation (given by the sampling points) that is
%       defined by the dimInfo.
%       Main use is within the geometry operations of MrImageGeometry.
%
% EXAMPLE
%
%   See also MrDimInfo MrImageGemeotry

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2018-10-19
% Copyright (C) 2018 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich

% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

% create template affine matrix (in case x, y, and z are not defined)
ADimInfo =  [1   0   0   0;
    0   1   0   0;
    0   0   1   0;
    0   0   0   1];

% check whether x, y, or z are specified
dimLabels = this.dimLabels;
hasX = ismember('x', dimLabels);
hasY = ismember('y', dimLabels);
hasZ = ismember('z', dimLabels);
if hasX, dimIndexX = this.get_dim_index('x'); end
if hasY, dimIndexY = this.get_dim_index('y'); end
if hasZ, dimIndexZ = this.get_dim_index('z'); end
% populate fields
% translation is the position of the first sampling point
% zoom is resolution or, if nan, samplingWidths
% x
if hasX
    ADimInfo(1,4) = this.samplingPoints{dimIndexX}(1);
    if ~isnan(this.resolutions(dimIndexX))
        ADimInfo(1,1) = this.resolutions(dimIndexX);
    elseif ~isnan(this.samplingWidths(dimIndexX))
        ADimInfo(1,1) = this.samplingWidths(dimIndexX);
    end
end

% y
if hasY
    ADimInfo(2,4) = this.samplingPoints{dimIndexY}(1);
    if ~isnan(this.resolutions(dimIndexY))
        ADimInfo(2,2) = this.resolutions(dimIndexY);
    elseif ~isnan(this.samplingWidths(dimIndexY))
        ADimInfo(2,2) = this.samplingWidths(dimIndexY);
    end
end

% z
if hasZ
    ADimInfo(3,4) = this.samplingPoints{dimIndexZ}(1);
    if ~isnan(this.resolutions(dimIndexZ))
        ADimInfo(3,3) = this.resolutions(dimIndexZ);
    elseif ~isnan(this.samplingWidths(dimIndexZ))
        ADimInfo(3,3) = this.samplingWidths(dimIndexZ);
    end
end
end

author: translationalneuromodeling

repo: tapas


text:
function [ path ] = full_path( db,enterie )
%FULL_PATH Summary of this function goes here
%   Detailed explanation goes here
   path=[db.path '/'];
   for ii=1:length(db.metaNames)
       path = [path enterie.(db.metaNames{ii}),'/'];
   end
   path =[path(1:end-1),'.',db.format];
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function node=orthdisk(c0,c1,r,ndiv)
%
% node=orthdisk(c0,c1,r,ndiv)
%
% Defining a 3D disk that is orthogonal to the vector c1-c0
%
% author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
%
% input:
%     c0: a 1x3 vector for the origin
%     c1: a 1x3 vector to define a direction vector c1-c0
%     r: the radius of the disk that is orthogonal to c1-c0, passing through c0
%     ndiv: division count to approximate a circle by a polygon, if ignored, ndiv=20
%
% output:
%     node: the 3D vertices of the disk
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

len=sqrt(sum((c0-c1).*(c0-c1)));
v0=c1-c0;

if(nargin<=2)
    r=1;
end
if(nargin<=3)
    ndiv=20;
end

dt=2*pi/ndiv;
theta=dt:dt:2*pi;
cx=r*cos(theta);
cy=r*sin(theta);
pp=[cx(:) cy(:) zeros(ndiv,1)];
node=rotatevec3d(pp,v0)+repmat(c0(:)',size(pp,1),1);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = gui_sim_eegdata(varargin)
%
% GUI_SIM_EEGDATA M-file for gui_sim_eegdata.fig
%      GUI_SIM_EEGDATA, by itself, creates a new GUI_SIM_EEGDATA or raises the existing
%      singleton*.
%
%      H = GUI_SIM_EEGDATA returns the handle to a new GUI_SIM_EEGDATA or the handle to
%      the existing singleton*.
%
%      GUI_SIM_EEGDATA('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in GUI_SIM_EEGDATA.M with the given input arguments.
%
%      GUI_SIM_EEGDATA('Property','Value',...) creates a new GUI_SIM_EEGDATA or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before gui_sim_eegdata_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to gui_sim_eegdata_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help gui_sim_eegdata

% Last Modified by GUIDE v2.5 13-Nov-2013 23:18:59

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @gui_sim_eegdata_OpeningFcn, ...
                   'gui_OutputFcn',  @gui_sim_eegdata_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before gui_sim_eegdata is made visible.
function gui_sim_eegdata_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to gui_sim_eegdata (see VARARGIN)

set(hObject,'name','Simulate EEG Data');

handles.output = hObject;

% set default termination behavior
handles.ExitButtonClicked = 'Cancel';

% set default EEGLAB background and text colors
%-----------------------------------------------
try, icadefs;
catch,
	GUIBACKCOLOR        =   [.8 .8 .8];
	GUIPOPBUTTONCOLOR   =   [.8 .8 .8];
	GUITEXTCOLOR        =   [0 0 0];
end;

allhandlers = hObject;

hh = findobj(allhandlers,'style', 'text');
%set(hh, 'BackgroundColor', get(hObject, 'color'), 'horizontalalignment', 'left');
set(hh, 'Backgroundcolor', GUIBACKCOLOR);
set(hh, 'foregroundcolor', GUITEXTCOLOR);
set(hObject, 'color',GUIBACKCOLOR );
% set(hh, 'horizontalalignment', g.horizontalalignment);

hh = findobj(allhandlers, 'style', 'edit');
set(hh, 'BackgroundColor', [1 1 1]); %, 'horizontalalignment', 'right');

hh =findobj(allhandlers, 'parent', hObject, 'style', 'pushbutton');
if ~strcmpi(computer, 'MAC') && ~strcmpi(computer, 'MACI') % this puts the wrong background on macs
    set(hh, 'backgroundcolor', GUIPOPBUTTONCOLOR);
    set(hh, 'foregroundcolor', GUITEXTCOLOR);
end;
hh =findobj(allhandlers, 'parent', hObject, 'style', 'popupmenu');
set(hh, 'backgroundcolor', GUIPOPBUTTONCOLOR);
set(hh, 'foregroundcolor', GUITEXTCOLOR);
hh =findobj(allhandlers, 'parent', hObject, 'style', 'checkbox');
set(hh, 'backgroundcolor', GUIBACKCOLOR);
set(hh, 'foregroundcolor', GUITEXTCOLOR);
hh =findobj(allhandlers, 'parent', hObject, 'style', 'listbox');
set(hh, 'backgroundcolor', GUIPOPBUTTONCOLOR);
set(hh, 'foregroundcolor', GUITEXTCOLOR);
hh =findobj(allhandlers, 'parent', hObject, 'style', 'radio');
set(hh, 'foregroundcolor', GUITEXTCOLOR);
set(hh, 'backgroundcolor', GUIPOPBUTTONCOLOR);
set(hObject, 'visible', 'on');

set(handles.pnlPropertyGrid,'backgroundcolor', GUIBACKCOLOR);
set(handles.pnlPropertyGrid,'foregroundcolor', GUITEXTCOLOR);
%-----------------------------------------------

drawnow

% render the PropertyGrid in the correct panel
handles.PropertyGridHandle = arg_guipanel( ...
                 handles.pnlPropertyGrid, ...
                'Function',@sim_eegdata, ...
                'params',{varargin{:}});

% Update handles structure
guidata(hObject, handles);

% Wait for user to click OK, Cancel or close figure
uiwait(handles.gui_sim_eegdata);

% --- Outputs from this function are returned to the command line.
function varargout = gui_sim_eegdata_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isempty(handles)
    % user closed the figure
    varargout = {[] hObject};
elseif strcmpi(handles.ExitButtonClicked,'OK')
    % user clicked OK
    % get PropertySpecification
    varargout = {handles.PropertyGridHandle handles.output};
else
    % user clicked cancel
    varargout = {[] handles.output};
end

try, close(hObject);
catch; end

function cmdCancel_Callback(hObject, eventdata, handles)

handles.ExitButtonClicked = 'Cancel';
guidata(hObject,handles);
uiresume(handles.gui_sim_eegdata)

function cmdOK_Callback(hObject, eventdata, handles)

handles.ExitButtonClicked ='OK';
guidata(hObject,handles);
uiresume(handles.gui_sim_eegdata);

function cmdHelp_Callback(hObject, eventdata, handles)
doc('sim_eegdata');

% --- Executes when gui_sim_eegdata is resized.
function gui_sim_eegdata_ResizeFcn(hObject, eventdata, handles)
% hObject    handle to gui_sim_eegdata (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

repo: REST


text:
function [Population,Range] = EnvironmentalSelection(Problem,Population,Range,N)
% The environmental selection of RSEA

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Non-dominated sorting
    [FrontNO,MaxFNO] = NDSort(Population.objs,N);
    Next             = find(FrontNO<=MaxFNO);

    %% Environmental selection
    if any(Range(1,:)==Range(2,:))
        Choose = LastSelection(Problem,Population(Next).objs,ismember(Next,find(FrontNO<MaxFNO)),N,ceil(sqrt(N)));
    else
        Choose = LastSelection(Problem,(Population(Next).objs-repmat(Range(1,:),length(Next),1))./repmat(Range(2,:)-Range(1,:),length(Next),1),ismember(Next,find(FrontNO<MaxFNO)),N,ceil(sqrt(N)));
    end
    Population = Population(Next(Choose));
	Range(1,:) = min([Range(1,:);Population.objs],[],1);
    Range(2,:) = max(Population(NDSort(Population.objs,1)==1).objs,[],1);
end

function Choose = LastSelection(Problem,PopObj,Choose,N,div)
% Select part of the solutions based on the radar grid

    %% Identify the extreme solutions
    [~,Extreme] = min(repmat(sqrt(sum(PopObj.^2,2)),1,size(PopObj,2)).*sqrt(1-(1-pdist2(PopObj,eye(size(PopObj,2)),'cosine')).^2),[],1); %Calculate the extreme points based on PBI
    Choose      = Choose | ismember(1:size(PopObj,1),Extreme);

    %% Calculate the convergence of each solution
	Con = sum(PopObj.^2,2).^0.5;
    Con = Con./max(Con);

    %% Calculate the radar grid of each solution
    [Site,RLoc] = RadarGrid(PopObj,div);
    RDis        = pdist2(RLoc,RLoc);
    RDis(logical(eye(length(RDis)))) = inf;
    CrowdG      = zeros(1,max(Site));
    temp        = tabulate(Site(Choose));
    CrowdG(temp(:,1)) = temp(:,2);

    %% Select N solutions one by one
    while sum(Choose) < N
        remainS  = find(~Choose);
        remainG  = unique(Site(remainS));
        bestG    = CrowdG(remainG) == min(CrowdG(remainG));
        current  = remainS(ismember(Site(remainS),remainG(bestG)));
        r        = 1-(Problem.FE/Problem.maxFE)^2;
        fitness  = Problem.M.*r.*Con(current) - min(RDis(current,Choose),[],2);
        [~,best] = min(fitness);
        Choose(current(best))       = true;
        CrowdG(Site(current(best))) = CrowdG(Site(current(best))) + 1;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function encseq=huffencode(huf,seq)
l=length(seq);
h=length(huf);
hufstr=struct2cell(huf);
hufstr=hufstr(1,1:h);
hufstr=cell2mat(hufstr);
encseq='';
 for i=1:l
     idx=find(seq(i)==hufstr);
     encseq=strcat(encseq,cell2mat(huf(idx).code));
 end
end


text:
function plot_his_map(groundTruthFile, expMapFile)
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     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/>.

    % load ground truth data
    [frameId, gt_x, gt_y, gt_z, gt_rx, gt_ry, gt_rz] = load_ground_truth_data(groundTruthFile);

    % load experience map data
    [expId, exp_x, exp_y, exp_z, exp_yaw, exp_pitch] = load_exp_map_data(expMapFile);

    hold on
    plot3(exp_x*1*(100), exp_y*100, exp_z*(-200), '.b');
    plot3((gt_x - gt_x(1))*100, (gt_y - gt_y(1))*100, (gt_z - gt_z(1))*(100), '.r');
    hold off;
    % grid on
    view(3)
    xl = xlabel('exp-x');
    yl = ylabel('exp-y');
    zlabel('exp-z');
    set(xl,'Rotation',15);
    set(yl,'Rotation',-30);
    title('3D Experience Map');
    %                     legend('Result','Truth' ,'1');
    % axis([-10 20 -10 20 -10 20]);
    %                     axis equal
    % axis on
    rotate3d on

end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
clear all
    % Generarea matricei A
A=[11,12,13,14;21,22,23,24;31,32,33,34;41,42,43,44]
    % Manipulari cu matricea A
B=A(2,2:4)
C=A(2,:)
D=A(1:2:3,:)
E(:,[3,5,7])=A(:,2:4)
F=A(:,[1 3 2 4])
    % Utilizarea filtrelor logice
Filt=zeros(size(A));
Filt(1,3)=1; Filt(2,4)=1; Filt(4,4)=1;
    % Transforarea matricei numerice in matrice logica
Filt=logical(Filt)
A_filtrat=A(Filt)'
    % Stergerea unei coloane din matricea A
A(:,2)=[ ]


text:
function [] = test_hexagram_shape

r1 = 2;
r2 = 1;

n = 6;
t = linspace(0, 2*pi, n+1);
x1 = r1 *[cos(t); sin(t) ];

t2 = t +pi /n;
x2 = r2 *[cos(t2); sin(t2) ];

m = 2 *n +1;
x(:, 1:2:m) = x1;
x(:, 2:2:m) = x2(:, 1:n);

ax = gca;
plotmd(ax, x)
axis(ax, 'equal')


text:
function train_timit_demo(context_win, hidden_units, num_layers, isdropout, isRNN, iscleanonly,...
    circular_step , isinputL1, MFCCorlogMelorSpectrum, framerate, pos_neg_r, outputnonlinear, opt, act, train_mode, const, const2, isGPU)
% Copyright (c) 2014-present University of Illinois at Urbana-Champaign
% All rights reserved.
%
% Developed by:     Po-Sen Huang, Paris Smaragdis
%                   Department of Electrical and Computer Engineering
%                   Department of Computer Science
%
%%
rand('state',0)
randn('state',0)
%%
%% setup paths for code. assumed this script runs in its own directory
baseDir= '../../';
codeDir = [baseDir,'codes', filesep];
minFuncDir = [baseDir, 'tools', filesep, 'minFunc_2012', filesep];

saveDir = [codeDir,filesep,'timit',...
           filesep,'discrim_joint_offset_all_results'];

addpath([baseDir, filesep,'tools', filesep,'labrosa']);
addpath([baseDir, filesep,'tools', filesep,'bss_eval']);

addpath(baseDir);
addpath(genpath(minFuncDir));
addpath([baseDir,filesep,'codes',filesep,'timit', filesep,'Data_with_dev']);
addpath(codeDir);
addpath([codeDir,'timit']);

CFGPath=[baseDir,'tools',filesep,'htk_features', filesep];
addpath(CFGPath);

%% setup network architecture
eI = [];
% 0- mfcc, 1- logmel, 2- spectrum
eI.MFCCorlogMelorSpectrum=MFCCorlogMelorSpectrum;
eI.CFGPath=CFGPath;
eI.seqLen = [1 10 25 50 100];
% eI.seqLen = [1];

eI.framerate=framerate;
if eI.framerate==64
    eI.winsize = 1024;    eI.nFFT = 1024;    eI.hop =eI.winsize/2;    eI.scf=1;%scf = 2/3;
else %32
    eI.winsize = 512;    eI.nFFT = 512;    eI.hop = eI.winsize/2;    eI.scf=1;%scf = 2/3;
end
winsize=eI.winsize; nFFT= eI.nFFT; hop= eI.hop; scf=eI.scf;
windows=sin(0:pi/winsize:pi-pi/winsize);

% single target or multiple targets
eI.cleanonly=iscleanonly;

% context window size of the input.
eI.num_contextwin = context_win;

% dimension of each input frame
if eI.MFCCorlogMelorSpectrum==0 %0 for mfcc 1 for logmel
    eI.featDim =39;
elseif eI.MFCCorlogMelorSpectrum==1 %0 for mfcc 1 for logmel
    eI.featDim =123;
else
    eI.featDim = (eI.nFFT/2+1);
end

eI.dropout = isdropout;

% weight tying in hidden layers
% if you want tied weights, must have odd number of *hidden* layers
eI.tieWeights = 0;

eI.const=const; eI.const2=const2;

hidden_units_set=[];

for il=1:num_layers
    hidden_units_set=[hidden_units_set, hidden_units];
end
% 2 hidden layers and output layer
if eI.cleanonly==1,
    eI.layerSizes = [hidden_units_set  nFFT/2+1 ];
else
    eI.layerSizes = [hidden_units_set  (nFFT/2+1)*2];
end
% highest hidden layer is temporal
eI.temporalLayer =isRNN;
% dim of network input at each timestep (final size after window & whiten)
eI.inputDim = eI.featDim * eI.num_contextwin;
% length of input sequence chunks.
% eI.seqLen = [1 10 25 50 100];
% activation function
switch act
    case 0
        eI.activationFn = 'logistic';
    case 1
        eI.activationFn = 'tanh';
    case 2
        eI.activationFn = 'RELU';
end

% temporal initialization type
eI.temporalInit = 'rand';
% weight norm penaly
eI.lambda = 0;
% file containing whitening matrices for outputs

eI.outputL1=0;

eI.inputL1=isinputL1;

eI.r=pos_neg_r;

eI.isdiscrim=2;

if opt==0,
    eI.opt='softlinear';
elseif opt==1,
    eI.opt='softabs';
elseif opt==2,
    eI.opt='softquad';
elseif opt==3,
    eI.opt='softabs_const';
elseif opt==4,
    eI.opt='softabs_kl_const';
end

eI.train_mode= train_mode;
eI.outputnonlinear=outputnonlinear;

%% setup weight caching
if isRNN,
    modelname=['model_RNN',num2str(isRNN)];
else
    modelname='model_DNN';
end
modelname=[modelname,'_win',num2str(context_win),'_h',num2str(hidden_units),'_l',num2str(num_layers)];
if iscleanonly, modelname=[modelname,'_cleanonly']; end
if isdropout,  modelname=[modelname,'_dropout'];    end
modelname=[modelname,['_r', num2str(eI.r)]];

modelname=[modelname,['_', num2str(eI.framerate),'ms']];
% modelname=[modelname,'_off',num2str(offset_step), '_snr', num2str(SNR_step)];
modelname=[modelname, '_', num2str(circular_step)];
eI.circular_step = circular_step;

modelname=[modelname,'_',eI.opt];
if outputnonlinear==0, modelname=[modelname,'_linearout']; end

modelname=[modelname, '_', eI.activationFn];

if eI.inputL1, modelname=[modelname, '_L',num2str(eI.inputL1)]; end

if eI.MFCCorlogMelorSpectrum==0
    modelname=[modelname,'_mfcc'];
elseif  eI.MFCCorlogMelorSpectrum==1
    modelname=[modelname,'_logmel'];
elseif  eI.MFCCorlogMelorSpectrum==2
    modelname=[modelname,'_spectrum'];
else
    modelname=[modelname,'_logpowspect'];
end

modelname= [modelname, '_trn', num2str(eI.train_mode)];

modelname=[modelname,'_c',num2str(const), '_c',num2str(const2)];

eI.modelname=modelname;
disp(modelname);

eI.saveDir = [saveDir, filesep, modelname, filesep];
if ~exist(eI.saveDir,'dir'), mkdir(eI.saveDir); end

%% initialize weights
[stack_i, W_t_i] = initialize_weights(eI);
[theta] = rnn_stack2params(stack_i, eI, W_t_i);

%% Directory of features
eI.featInBase =baseDir;

%% load data
eI.useCache = 0;

%% setup minFunc
options.Diagnostics = 'on';
options.Display = 'iter';
options.MaxIter = 800;
options.MaxFunEvals = 2500;
options.Corr = 50;
options.DerivativeCheck = 'off';
options.outputFcn = @save_callback_timit_general;

%% compute feature
SNRs=0;

[train1, fs, nbits]=wavread('female_train.wav');
[train2, fs, nbits]=wavread('male_train.wav');

maxLength=max([length(train1), length(train2)]);
train1(end+1:maxLength)=eps;
train2(end+1:maxLength)=eps;

train1=train1./sqrt(sum(train1.^2));
train2=train2./sqrt(sum(train2.^2));

eI.fs=fs;
%%
% chunk
[data_cell, targets_cell, mixture_spectrum]=formulate_data(train1, train2, eI, eI.train_mode); %0 -- chunk, 2--no chunk

  global SDR;
  SDR.deviter=0;   SDR.devmax=0;   SDR.testmax=0;

  eI.writewav=0;

if isGPU==1
  [theta,val]=minFunc(@drdae_discrim_joint_kl_obj_gpu, theta, options, eI, data_cell, targets_cell, mixture_spectrum, false, false);
else
  [theta,val]=minFunc(@drdae_discrim_joint_kl_obj, theta, options, eI, data_cell, targets_cell, mixture_spectrum, false, false);
end

  fprintf('%s\tdevmaxiter:\t%d\tdevSDR:\t%.3f\ttestSDR:\t%.3f\n',modelname, SDR.deviter, SDR.devmax, SDR.testmax);

return;

%% unit test

% context window size
context_win = 1;
% hidden units
hidden_units = 16;
num_layers = 1;
isdropout = 0;
% RNN temporal connection
isRNN = 2;
% One output source or two
iscleanonly = 0;
% Circular shift step
circular_step = 10000;
% normalize input as L1 norm = 1
isinputL1 = 0;
% 0: MFCC, 1: logmel, 2: spectra
MFCCorlogMelorSpectrum = 2;
% feature frame rate
framerate = 64;
% discriminative training gamma parameter
pos_neg_r = 0.05;
% Last layer - linear or nonlinear
outputnonlinear = 0;
% soft mask obj
softabs = 1;
% 0: logistic, 1: tanh, 2: RELU
act = 2;
% constant for avoiding numerical problems
const = 1e-10;
% constant for avoiding numerical problems
const2 = 0.001;
% 0: not using GPU, 1: using GPU
isGPU = 0;

train_mode = 0;
% 0:'softlinear',1:'softabs', 2:'softquad', 3:'softabs_const',
% 4:'softabs_kl_const'
opt = 1;

train_timit_demo(context_win, hidden_units, num_layers, isdropout, isRNN, iscleanonly,...
    circular_step , isinputL1, MFCCorlogMelorSpectrum, framerate, pos_neg_r, ...
    outputnonlinear, opt, act, train_mode, const, const2, isGPU)

author: posenhuang

repo: deeplearningsourceseparation


text:
% STD_READTOPO - returns the scalp map of a specified ICA component, assumed
%                  to have been saved in a Matlab file, [dataset_name].icatopo,
%                  in the same directory as the dataset file. If this file does
%                  not exist, use STD_TOPO to create it, else a pre-clustering
%                  function that calls it: POP_PRECLUST or EEG_PRECLUST.
% Usage:
%   >> [grid, y, x ] = std_readtopo(ALLEEG, setindx, component);
%   >> [grid, y, x ] = std_readtopo(ALLEEG, setindx, component, transform, mode);
%
% Inputs:
%   ALLEEG     - vector of EEG datasets (can also be one EEG set).
%                must contain the dataset of interest (see 'setindx' below).
%   setindx    - [integer] an index of an EEG dataset in the ALLEEG
%                structure, for which to get the component ERP.
%   component  - [integer] index of the component for which the scalp map
%                grid should be returned.
%   transform  - ['none'!'laplacian'|'gradient'] transform scalp map to
%                laplacian or gradient map. Default is 'none'.
%   mode       - ['2dmap'|'preclust'] return either a 2-D array for direct
%                plotting ('2dmap') or an array formatted for preclustering
%                with all the NaN values removed (ncomps x points). Default
%                is '2dmap' for 1 component and 'preclust' for several.
%
% Outputs:
%   grid      - square scalp-map color-value grid for the requested ICA component
%               in the specified dataset, an interpolated Cartesian grid as output
%               by TOPOPLOT.
%   y         - y-axis values for the interpolated grid
%   x         - x-axis values of the interpolated grid
%
%  See also  STD_TOPO, STD_PRECLUST
%
% Authors: Arnaud Delorme, Hilit Serby, SCCN, INC, UCSD, February, 2005

% Copyright (C) Hilit Serby, SCCN, INC, UCSD, October 11, 2004, hilit@sccn.ucsd.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.

function [X, yi, xi ] = std_readtopo(ALLEEG, abset, comps, option, mode)

X = [];
yi = [];
xi = [];
if nargin < 4
    option = 'none';
end
if nargin < 5
    mode = '2Dmap';
end
filename = correctfile(fullfile( ALLEEG(abset).filepath,[ ALLEEG(abset).filename(1:end-3) 'icatopo']),ALLEEG(abset).filepath);
tmpfile  = which(filename);
if ~isempty(tmpfile), filename = tmpfile; end

% 061411, 2:51pm
% Modified by Joaquin
% while getfield(dir(filename), 'bytes') < 1000
i = 1;
while getfield(dir(filename), 'bytes') < 5000
    topo = load( '-mat', filename);
    filename = correctfile(topo.file, ALLEEG(abset).filepath);
    tmpfile  = which(filename);
    if ~isempty(tmpfile), filename = tmpfile; end
    if(i>100)
        error('too many attempts to find valid icatopo');
    end
    i = i+1;
end

for k = 1:length(comps)

    if length(comps) < 3
        lastwarn('', '');
        try
            topo = load( '-mat', filename, ...
                         [ 'comp' int2str(comps(k)) '_grid'], ...
                         [ 'comp' int2str(comps(k)) '_x'], ...
                         [ 'comp' int2str(comps(k)) '_y'] );
        catch
            error( [ 'Cannot read file ''' filename '''' ]);
        end
        [~, warnId] = lastwarn();
        if ~isempty(warnId)
            error( 'Cannot find component %d in file %s', comps(k),  filename );
        end
    elseif k == 1
        try
            topo = load( '-mat', filename);
        catch
            error([ 'Missing scalp topography file - also necessary for ERP polarity' 10 'Try recomputing scalp topographies for components' ]);
        end
    end

    try
        tmp =  getfield(topo, [ 'comp' int2str(comps(k)) '_grid' ]);
    catch
        error([ 'Empty scalp topography file - also necessary for ERP polarity' 10 'Try recomputing scalp topographies for components' ]);
    end

    if strcmpi(option, 'gradient')
        [tmpx, tmpy]  = gradient(tmp); % Gradient
        tmp        = tmpx;
        tmp(:,:,2) = tmpy;
    elseif strcmpi(option, 'laplacian')
        tmp = del2(tmp); % Laplacian
    end

    if length(comps) > 1 || strcmpi(mode, 'preclust')
        tmp = tmp(find(~isnan(tmp))); % remove NaN for more than 1 component
    end
    if k == 1
        X = zeros([ length(comps) size(tmp) ]) ;
    end
    X(k,:,:,:) =  tmp;
    if k == 1
        yi   = getfield(topo, [ 'comp' int2str(comps(k)) '_y']);
        xi   = getfield(topo, [ 'comp' int2str(comps(k)) '_x']);
    end
end
X = squeeze(X);

return;

function filename = correctfile(filename, datasetpath)
    comp = computer;
    if filename(2) == ':' && ~strcmpi(comp(1:2), 'PC')
        filename = [filesep filename(4:end) ];
        filename(find(filename == '\')) = filesep;
    end

    if ~exist(filename)
        [tmpp tmpf ext] = fileparts(filename);
        if exist([tmpf ext])
            filename = [tmpf ext];
        else
            [tmpp2 tmpp1] = fileparts(tmpp);
            if exist(fullfile(tmpp1, [ tmpf ext ]))
                filename = fullfile(tmpp1, [ tmpf ext ]);
            else
                filename = fullfile(datasetpath, [ tmpf ext ]);
                if ~exist(filename)
                    error([ 'Cannot load file ''' [ tmpf ext ] '''' 10 'Go back and recompute the data file.' 10 'Note that plotting ICA component ERPs require' 10 'to precompute ICA topographies (see tutorial)']);
                end
            end
        end
    end;


text:
function tests = test_ft_preproc_slidingrange

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_preproc_slidingrange

if nargout
  % assume that this is called by RUNTESTS
  tests = functiontests(localfunctions);
else
  % assume that this is called from the command line
  fn = localfunctions;
  for i=1:numel(fn)
    feval(fn{i});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function testOptions(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nchan   = 8;
nsample = 1000;
dat     = randn(nchan, nsample) + 1;

result = [];
result{end+1} = ft_preproc_slidingrange(dat, 1,  false);
result{end+1} = ft_preproc_slidingrange(dat, 11, false);
result{end+1} = ft_preproc_slidingrange(dat, 31, false);
result{end+1} = ft_preproc_slidingrange(dat, 71, false);
result{end+1} = ft_preproc_slidingrange(dat, 11, true);
result{end+1} = ft_preproc_slidingrange(dat, 31, true);
result{end+1} = ft_preproc_slidingrange(dat, 71, true);

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequal(result{i}, result{j}), 'the results %d and %d should not be equal', i, j);
  end
end


text:
function r8vec_undex_test ( )

%*****************************************************************************80
%
%% R8VEC_UNDEX_TEST tests R8VEC_UNDEX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 October 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_UNDEX_TEST\n' );
  fprintf ( 1, '  R8VEC_UNDEX produces index vectors which create a sorted\n' );
  fprintf ( 1, '  list of the unique elements of an (unsorted) R8VEC, \n' );
  fprintf ( 1, '  and a map from the original vector to the (implicit) \n' );
  fprintf ( 1, '  vector of sorted unique elements.\n' );

  x_num = 9;
  x_val = [ 33.0, 55.0, 11.0, 11.0, 55.0, 33.0, 22.0, 22.0, 11.0 ];

  r8vec_print ( x_num, x_val, '  The vector X:' );

  tol = r8_epsilon ( );
  x_unique_num = r8vec_unique_count ( x_num, x_val, tol );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Tolerance for equality is %e\n', tol );
  fprintf ( 1, '  Number of unique entries in X is %d\n', x_unique_num );

  [ undx, xdnu ] = r8vec_undex ( x_num, x_val, x_unique_num, tol );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  UNDX can be used to list the unique elements of X\n' );
  fprintf ( 1, '  in sorted order.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I  UNDX   X(UNDX)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : x_unique_num
    fprintf ( 1, '  %4d  %4d  %8f\n', i, undx(i), x_val(undx(i)) );
  end

  xu_val(1:x_unique_num) = x_val(undx(1:x_unique_num));

  fprintf ( 1, '\n' );
  fprintf ( 1, '  UNDX can be used to created XU, a copy of X\n' );
  fprintf ( 1, '  containing only the unique elements, in sorted order.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I  UNDX     XU(I)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : x_unique_num
    fprintf ( 1, '  %4d  %4d  %8f\n', i, undx(i), xu_val(i) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  XDNU can be used to match each element of X with one of the\n' );
  fprintf ( 1, '  unique elements\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I  XDNU    X(I)       XU(XDNU(I))\n' );
  fprintf ( 1, '\n' );

  for i = 1 : x_num
    fprintf ( 1, '  %4d  %4d  %8f  %12f\n', i, xdnu(i), x_val(i), xu_val(xdnu(i)) );
  end

  return
end


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: various distances versus rotation angle
%
%   - data                 PETCT, Omega=(0,140)x(0,151), level=4:7, m=[128,128]
%   - viewer               viewImage2D
%   - interpolation        splineInter
%   - distance             {'SSD','NCC','MI','NGF'}
%   - transformation       rotation2D
% see also E7_extended
%==============================================================================

clear, close all, help(mfilename);

fprintf('%s\n','setup data, viewer, interpolator, transformation model, distance')
setup2DPETCTData; level = 6; omega = ML{level}.omega; m = ML{level}.m;

viewImage('reset',viewPara{:},'axis','off');
imgModel('reset','imgModel','splineInter','regularizer','moments','theta',1e0);
[T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega,'out',0);
center = (omega(2:2:end)-omega(1:2:end))'/2;
trafo('reset','trafo','rotation2D','c',center);
distance('reset','distance','SSD');
fprintf('%20s : %s\n','viewImage',viewImage);
fprintf('%20s : %s\n','imgModel',imgModel);
fprintf('%20s : %s\n','trafo',trafo);

xc  = getCellCenteredGrid(omega,m);
Rc = imgModel(R,omega,xc);

distances = {'SSD','NCC','MI','NGF'};

for k=1:length(distances),
  distance('reset','distance',distances{k});
  fprintf('%20s : %s\n','distance',distance);
  wc = linspace(-pi/2,pi/2,49);
  Dc = zeros(size(wc));

  % run the loop over all rotations
  for j = 1:length(wc),
    yc    = trafo(wc(j),xc);         % compute transformed grid
    Tc    = imgModel(T,omega,yc);    % compute transformed image
    Dc(j) = distance(Tc,Rc,omega,m); % compute distance

    % visualize
    if j == 1,
      th = [];
      FAIRfigure(k,'figname',mfilename); clf;
      subplot(1,3,1);  viewImage(Rc,omega,m);            th(1) = title('R');
      subplot(1,3,2);  vh = viewImage(Tc,omega,m);       th(2) = title('T(yc)');
      subplot(1,3,3);  ph = plot(wc(1),Dc(1),'r.','markersize',20);
      th(3) = title(sprintf('%s versus rotation',distance));
      axis([wc(1),wc(end),-inf,inf]); hold on; set(th,'fontsize',30);
      axis('auto y');
      FAIRpause;
    else
      set(vh,'cdata',reshape(Tc,m)')
      subplot(1,3,3); set(ph,'visible','off');
      plot(wc(1:j),Dc(1:j),'k-','linewidth',2);
      ph = plot(wc(j),Dc(j),'r.','markersize',20);
      drawnow,
      FAIRpause(1/100)
    end;
    fprintf('.'); if ~rem(j,50) || j == length(wc), fprintf('\n'); end;
  end;
  FAIRpause;
end;

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


text:
function regions = read_scg_regions(filePath, varargin)
opts.maxNumRegions = +inf ;
opts = vl_argparse(opts, varargin) ;

tmp = load(filePath) ;
regions.basis  = tmp.candidates_scg.superpixels ;
regions.labels = tmp.candidates_scg.labels(:)' ;
regions.scores = tmp.candidates_scg.scores(:)' ;

sareas = accumarray(regions.basis(:), 1) ;
for i=1:numel(regions.labels)
  regions.areas(i) = sum(sareas(regions.labels{i})) ;
end

% take only the top regions
[~,perm] = sort(regions.scores, 'descend') ;
perm = vl_colsubset(perm, opts.maxNumRegions, 'beginning') ;
regions.labels = regions.labels(perm) ;
regions.scores = regions.scores(perm) ;
regions.areas = regions.areas(perm) ;

author: mcimpoi

repo: deep-fbanks

save_path: github-repos/MATLAB/mcimpoi-deep-fbanks


text:
function [sts, vind] = checksubs_job(item, subs, dflag)

% function [sts vind] = checksubs_job(item, subs, dflag)
% Check whether a subscript reference is a valid reference in a job
% structure starting at item. subs(1) should have a subscript type of
% '.', and the subscript reference should be a tagname from item.val or
% item.values, depending on dflag.
%
% This function is identical for cfg_branch and cfg_(m)choice classes.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: checksubs_job.m 5678 2013-10-11 14:58:04Z volkmar $

rev = '$Rev: 5678 $'; %#ok

if strcmp(subs(1).type, '.')
    vind = find(strcmp(subs(1).subs, tagnames(item, dflag)),1);
    sts  = ~isempty(vind);
else
    sts = false;
    vind = [];
end

author: spm


text:
function [w]=ADMM_solve_w(params,use_sz,model_w,h_f)

    w = gpuArray(params.w_init*single(ones(use_sz)));
    q = w;
    m = w;

    mu    = 1;
    betha = 10;
    mumax = 10000;
    i = 1;
    params.admm_iterations=2;
    T = prod(use_sz);
    h=T*real(ifft2(h_f));
%     hw=bsxfun(@times,h,model_w);
    hw=h;
    Hh=sum(hw.^2,3);
%     Hh=sum(real(bsxfun(@times,h_f(:),conj(h_f(:)))));

    %   ADMM
    while (i <= params.admm_iterations)
        %   solve for w- please refer to the paper for more details
      %  w = (q-m)/(1+(params.admm_lambda1/mu)*Hh);
        w = bsxfun(@rdivide,(q-m),(1+(params.admm_lambda1/mu)*Hh));
        %   solve for q
        q=(params.admm_lambda2*model_w + mu*(w+m))/(params.admm_lambda2 + mu);

        %   update m
        m = m + (w - q);

        %   update mu- betha = 10.
        mu = min(betha * mu, mumax);
        i = i+1;

    end

end

author: Daikenan

repo: ASRCF

save_path: github-repos/MATLAB/Daikenan-ASRCF


text:
function out = forces_moments_wing(x, delta, wind, ~, p_drone, p_physics)

% FORCES_MOMENTS - Computes the forces and moments acting on the UAV.
%
% Inputs:
%   x     - state variables
%   delta - actuators values (angles,...)
%   wind  - wind parameters
%   P     - general parameters
% Outputs:
%   out
%       F     - forces
%       M     - moments
%       Va    - airspeed
%       alpha - angle of attack
%       beta  - sideslip angle
%       wind  - wind vector in the inertial frame
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Variables and Inputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% UAV velocity wrt inertial frame in body frame
v_xyz      = x(4:6);

% Euler angles
phi         = x(7);
theta       = x(8);
psi         = x(9);

% Rotation rates
p           = x(10);
q           = x(11);
r           = x(12);

% Actuators fixed-wing
de = delta(1); % delta elevator
da = delta(2); % delta aileron
dr = delta(3); % delta rudder
dt = delta(4); % delta thrust

Rbi = Rb2i(phi,theta,psi);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Wind computation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ws_ned  = wind(1:3);  % steady wind in NED frame
wg_xyz  = wind(4:6);  % gusts along body xyz axis

% Gust in NED frame = Rib * gust in Body frame
wg_ned = Rbi * wg_xyz;

% Wind in NED = steady in NED + gust in NED
w_ned = ws_ned + wg_ned;
wn = w_ned(1);
we = w_ned(2);
wd = w_ned(3);

%Steady wind in body (xyz) frame
Rib = Rbi';
ws_xyz = Rib * ws_ned;

% Wind in body frame
w_xyz = ws_xyz + wg_xyz;

% Compute air data wrt the inertial frame, in the body frame
va_xyz = v_xyz - w_xyz;
vax = va_xyz(1);
vay = va_xyz(2);
vaz = va_xyz(3);

% Airspeed norm (=V_inifity)
Va = norm(va_xyz);

if vax ~= 0
    alpha = atan2(vaz, vax);
else
    alpha = 0;
end

if Va ~=0
    beta = asin(vay / Va);
else
    beta = 0;
end

%   Relabel
ca = cos(alpha);
sa = sin(alpha);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Compute forces on fixed wing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   Aerodynamic coefficients in function of alpha
%       Lift and Drag Coefficients
CL = p_drone.CL0 + p_drone.CL_alpha * alpha;
CD = p_drone.CD0 + p_drone.CD_alpha * alpha;

%       Lift and Drag expressed in body frame
Cx     = - CD * ca       + CL * sa;
Cx_q   = - p_drone.CD_q * ca   + p_drone.CL_q  * sa;
Cx_de  = - p_drone.CD_de * ca  + p_drone.CL_de * sa;
Cz     = - CD * sa       - CL * ca;
Cz_q   = - p_drone.CD_q  * sa  - p_drone.CL_q  * ca;
Cz_de  = - p_drone.CD_de * sa  - p_drone.CL_de * ca;

% Weight
weight = p_drone.mass * p_physics.gravity * [-sin(theta);...
    cos(theta) * sin(phi);...
    cos(theta) * cos(phi)];

% Aerodynamic forces
p_dyn = 0.5 * p_physics.rho * Va^2;

if Va ~= 0
    kc = 0.5*p_drone.c/Va;
    kb = 0.5*p_drone.b/Va;
else
    kc = 0;
    kb = 0;
end

f_drag = p_dyn * p_drone.S_wing * ...
    [Cx + Cx_q*kc*q + Cx_de*de; ...
    p_drone.CY0 + p_drone.CY_beta*beta + p_drone.CY_p*kb*p + p_drone.CY_r*kb*r + p_drone.CY_da*da + p_drone.CY_dr*dr;...
    Cz + Cz_q*kc*q + Cz_de*de];

% Thrust force
f_thrust = 0.5*p_physics.rho*p_drone.S_prop*p_drone.C_prop*...
    [(p_drone.k_motor*dt)^2-Va^2; 0; 0];

% Total force = Weight + Aerodynamic forces + Thrust force
f_tot = weight + f_drag + f_thrust;

% Compute torques on drone
torque_aerodyn = p_dyn * p_drone.S_wing * ...
    [p_drone.b*(p_drone.Cl0 + p_drone.Cl_beta*beta + p_drone.Cl_p*kb*p + p_drone.Cl_r*kb*r + p_drone.Cl_da*da + p_drone.Cl_dr*dr);...
    p_drone.c*(p_drone.Cm0 + p_drone.Cm_alpha*alpha + p_drone.Cm_q*kc*q + p_drone.Cm_de*de);...
    p_drone.b*(p_drone.Cn0 + p_drone.Cn_beta*beta + p_drone.Cn_p*kb*p + p_drone.Cn_r*kb*r + p_drone.Cn_da*da + p_drone.Cn_dr*dr)];

% Thrust torque
%   Relabel thrust coefficients
torque_thrust = [-p_drone.k_TP*(p_drone.k_omega*dt)^2; 0; 0];

% Total torque
torque_tot = torque_aerodyn + torque_thrust;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Create output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

out = [f_tot; torque_tot; Va; alpha; beta; wn; we; wd];

end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab