text:
function resultGUI = matRad_calcDoseDirectMC(ct,stf,pln,cst,w,nHistories)
% matRad function to bypass dij calculation for MC dose calculation
% matRad dose calculation wrapper for MC dose calculation algorithms
% bypassing dij calculation for MC dose calculation algorithms.
%
% call
%   resultGUI = matRad_calcDoseDirecMC(ct,stf,pln,cst)
%   resultGUI = matRad_calcDoseDirecMC(ct,stf,pln,cst,w)
%   resultGUI = matRad_calcDoseDirectMC(ct,stf,pln,cst,nHistories)
%   resultGUI = matRad_calcDoseDirectMC(ct,stf,pln,cst,w,nHistories)
%
% input
%   ct:         ct cube
%   stf:        matRad steering information struct
%   pln:        matRad plan meta information struct
%   cst:        matRad cst struct
%   w:          (optional, if no weights available in stf): bixel weight
%               vector
%   nHistories: (optional) number of histories
%
% output
%   resultGUI:  matRad result struct
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2019 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

matRad_cfg =  MatRad_Config.instance();

calcDoseDirect = true;

if nargin < 6 || ~exist('nHistories')
  nHistories = matRad_cfg.propMC.direct_defaultHistories;
  matRad_cfg.dispInfo('Using default number of Histories: %d\n',nHistories);
end

% check if weight vector is available, either in function call or in stf - otherwise dose calculation not possible
if ~exist('w','var') && ~isfield([stf.ray],'weight')
     matRad_cfg.dispError('No weight vector available. Please provide w or add info to stf');
end

% copy bixel weight vector into stf struct
if exist('w','var')
    if sum([stf.totalNumOfBixels]) ~= numel(w)
        matRad_cfg.dispError('weighting does not match steering information');
    end
    counter = 0;
    for i = 1:size(stf,2)
        for j = 1:stf(i).numOfRays
            for k = 1:stf(i).numOfBixelsPerRay(j)
                counter = counter + 1;
                stf(i).ray(j).weight(k) = w(counter);
            end
        end
    end
else % weights need to be in stf!
    w = NaN*ones(sum([stf.totalNumOfBixels]),1);
    counter = 0;
    for i = 1:size(stf,2)
        for j = 1:stf(i).numOfRays
            for k = 1:stf(i).numOfBixelsPerRay(j)
                counter = counter + 1;
                w(counter) = stf(i).ray(j).weight(k);
            end
        end
    end
end

% dose calculation
if strcmp(pln.radiationMode,'protons')
  dij = matRad_calcParticleDoseMC(ct,stf,pln,cst,nHistories,calcDoseDirect);
else
    matRad_cfg.dispError('Forward MC only implemented for protons.');
end

% hack dij struct
dij.numOfBeams = 1;
dij.beamNum = 1;

% calculate cubes; use uniform weights here, weighting with actual fluence
% already performed in dij construction
resultGUI    = matRad_calcCubes(sum(w),dij);

% remember original fluence weights
resultGUI.w  = w;

repo: matRad


text:
function baseName=fileBaseName(fileName)

[~,baseName]=myfileparts(fileName);

end

author: kristinbranson

repo: JAABA


text:
function p = predictOneVsAll(all_theta, X)
%PREDICT Predict the label for a trained one-vs-all classifier. The labels
%are in the range 1..K, where K = size(all_theta, 1).
%  p = PREDICTONEVSALL(all_theta, X) will return a vector of predictions
%  for each example in the matrix X. Note that X contains the examples in
%  rows. all_theta is a matrix where the i-th row is a trained logistic
%  regression theta vector for the i-th class. You should set p to a vector
%  of values from 1..K (e.g., p = [1; 3; 1; 2] predicts classes 1, 3, 1, 2
%  for 4 examples)

m = size(X, 1);
num_labels = size(all_theta, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters (one-vs-all).
%               You should set p to a vector of predictions (from 1 to
%               num_labels).
%
% Hint: This code can be done all vectorized using the max function.
%       In particular, the max function can also return the index of the
%       max element, for more information see 'help max'. If your examples
%       are in rows, then, you can use max(A, [], 2) to obtain the max
%       for each row.
%

z = X * all_theta';

h = sigmoid(z);

[pval p] = max(h, [], 2);

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

end

author: fewtime

repo: ML

save_path: github-repos/MATLAB/fewtime-ML


text:
  function y = downsample3(x, m)
%|function y = downsample3(x, m)
%| downsample by averaging by integer factors
%| m can be a scalar (same factor for both dimensions)
%| or a 3-vector
%| in
%|	x	[nx ny nz]
%|	m	[1] or [1 3]
%| out
%|	y	[nx/m ny/m nz/m]
%| function modified to 3D space by: Taka Masuda

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

if ndims(x) == 2
	warn('2d case')
	if numel(m) == 1, m = [m m 1]; end
	if numel(m) ~= 3, fail 'm not length 3?', end
	y = downsample2(x, m(1:2));
return
end

if ndims(x) ~= 3
	fail('not 3d')
end

if numel(m) == 1
	m = m * ones(ndims(x),1);
end
if numel(m) ~= ndims(x), error 'bad m', end

% downsample along each dimension
y = downsample1(x, m(1));
y = downsample1(permute(y, [2 1 3]), m(2));
y = downsample1(permute(y, [3 2 1]), m(3)); % [3 1 2] order
y = permute(y, [2 3 1]);

% downsample1hide()
% 1d down sampling of each column
function y = downsample1hide(x, m)
'ok'
if m == 1, y = x; return; end
n1 = floor(size(x,1) / m);
n2 = size(x,2);
n3 = size(x,3);
y = zeros(n1,n2,n3);
for ii=0:m-1
	y = y + x(ii+[1:m:m*n1],:,:);
	ticker(mfilename, ii+1, m)
end
y = y / m;

% downsample3_test
function downsample3_test
x = [6 5 2];
x = reshape(2*[1:prod(x)], x);
down = 2;
y = downsample3(x, down);
if down == 1
	jf_equal(y, x)
else
	jf_equal(y, [39 63; 43 67; 47 71])
end

if has_aspire
	filex = [test_dir filesep 'testx.fld'];
	filey = [test_dir filesep 'testy.fld'];
	fld_write(filex, x)
	delete(filey)
	com = ['op sample3 mean ' filey ' ' filex ' %d %d %d'];
	com = sprintf(com, down, down, down);
	os_run(com)
	z = fld_read(filey);
	if ~isequal(y, z), error 'aspire/matlab mismatch', end
end

if 0 % big
	x = zeros(2.^[7 8 9]);
	cpu etic
	y = downsample3(x, 2);
	cpu etoc 'downsample3 time:'
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
load('history_734828928651');

n=1000;

r=hisrr(1,1:n,:);
hp=plot(r(1,:,1),r(1,:,2),'k.');
axis equal;
for lc=1:size(hisrr,2)/n
    r=hisrr(1,1+n*(lc-1):n*lc,:);
    %hp=plot(r(1,:,1),r(1,:,2),'k.');
    set(hp,'Xdata',r(1,:,1),'Ydata',r(1,:,2));
    drawnow;
end

author: Sable

repo: mcbench-benchmarks


text:
function L = kernPriorLogProb(kern)

% KERNPRIORLOGPROB Compute penalty terms associated with kernel priors.

% KERN

L = 0;
switch kern.type
 case {'cmpnd', 'multi', 'tensor'}
  for i = 1:length(kern.comp)
    L = L + kernPriorLogProb(kern.comp{i});
  end
 otherwise
  if isfield(kern, 'priors')
    fhandle = str2func([kern.type 'KernExtractParam']);
    params = fhandle(kern);
    if iscell(kern.priors),
      for i = 1:length(kern.priors)
	index = kern.priors{i}.index;
	L = L + priorLogProb(kern.priors{i}, params(index));
      end
    else
      for i = 1:length(kern.priors)
	index = kern.priors(i).index;
	L = L + priorLogProb(kern.priors(i), params(index));
      end
    end
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [spike] = ft_datatype_spike(spike, varargin)

% FT_DATATYPE_SPIKE describes the FieldTrip MATLAB structure for spike data
%
% Spike data is obtained using FT_READ_SPIKE to read files from a Plexon,
% Neuralynx or other animal electrophysiology data acquisition system. It
% is characterised as a sparse point-process, i.e. each neuronal firing is
% only represented as the time at which the firing happened. Optionally,
% the spike waveform can also be represented. Using this waveform, the
% neuronal firing events can be sorted into their single units.
%
% A required characteristic of the SPIKE structure is a cell-array with the
% label of the (single or multi) units.
%
%         label: {'unit1'  'unit2'  'unit3'}
%
% The fields of the SPIKE structure that contain the specific information
% per spike depends on the available information. A relevant distinction
% can be made between the representation of raw spikes that are not related
% to the temporal structure of the experimental design (i.e trials), and
% the data representation in which the spikes are related to the trial.
%
% For a continuous recording the SPIKE structure must contain a cell-array
% with the raw timestamps as recorded by the hardware system. As example,
% the original content of the .timestamp field can be
%
%         timestamp:  {[1x504 uint64]  [1x50 uint64]  [1x101 uint64]}
%
% An optional field that is typically obtained from the raw recording
% contains the waveforms for every unit and label as a cell-array. For
% example, the content of this field may be
%
%         waveform:   {[1x32x504 double] [1x32x50 double] [1x32x101 double]}
%
% If the data has been organised to reflect the temporal structure of the
% experiment (i.e. the trials), the SPIKE structure should contain a
% cell-array with the spike times relative to an experimental trigger. The
% FT_SPIKE_MAKETRIALS function can be used to reorganise the SPIKE
% structure such that the spike times are expressed relative to a trigger
% instead of relative to the acquisition devices internal timestamp clock.
% The time field then contains only those spikes that ocurred within one of
% the trials . The spike times are now expressed on seconds relative to the
% trigger.
%
%         time:       {[1x504 double] [1x50 double] [1x101 double]}
%
% In addition, for every spike we register in which trial the spike was
% recorded:
%
%         trial:      {[1x504 double] [1x50 double] [1x101 double]}
%
% To fully reconstruct the structure of the spike-train, it is required
% that the exact start- and end-point of the trial (in seconds) is
% represented. This is specified in a nTrials x 2 matrix.
%
%         trialtime:  [100x2 double]
%
% As an example, FT_SPIKE_MAKETRIALS could result in the following
% SPIKE structure that represents the spikes of three units that were
% observed in 100 trials:
%
%         label:           {'unit1'  'unit2'  'unit3'}
%         timestamp:       {[1x504 double] [1x50 double] [1x101 double]}
%         timestampdimord: '{chan}_spike'
%         time:            {[1x504 double] [1x50 double] [1x101 double]}
%         trial:           {[1x504 double] [1x50 double] [1x101 double]}
%         trialtime:       [100x2 double]
%         sampleinfo:      [100x2 double]
%         waveform:        {[1x32x504 double] [1x32x50 double] [1x32x101 double]}
%         waveformdimord: '{chan}_lead_time_spike'
%
% For analysing the relation between the spikes and the local field
% potential (e.g. phase-locking), the SPIKE structure can have additional
% fields such as fourierspctrm, lfplabel, freq and fourierspctrmdimord.
%
% For example, from the structure above we may obtain
%
%         label:               {'unit1'  'unit2'  'unit3'}
%         time:                {[1x504 double] [1x50 double] [1x101 double]}
%         trial:               {[1x504 double] [1x50 double] [1x101 double]}
%         trialtime:           [100x2 double]
%         timestamp:           {[1x504 double] [1x50 double] [1x101 double]}
%         timestampdimord:     '{chan}_spike'
%         waveform:            {[1x32x504 double] [1x32x50 double] [1x32x101 double]}
%         waveformdimord:      '{chan}_lead_time_spike'
%         fourierspctrm:       {504x2x20, 50x2x20, 101x2x20}
%         fourierspctrmdimord: '{chan}_spike_lfplabel_freq'
%         lfplabel:            {'lfpchan1', 'lfpchan2'}
%         freq:                [1x20 double]
%
% Required fields:
%   - label
%   - timestamp
%
% Optional fields:
%   - time, trial, trialtime
%   - timestampdimord
%   - unit, unitdimord
%   - waveform, waveformdimord
%   - fourierspctrm, fourierspctrmdimord, freq, lfplabel (these are extra outputs from FT_SPIKETRIGGEREDSPECTRUM and FT_SPIKE_TRIGGEREDSPECTRUM)
%   - hdr
%   - cfg
%
% Deprecated fields:
%   - origtime, origtrial
%
% Obsoleted fields:
%   - <unknown>
%
% Revision history:
%
% (2020/latest) Add an explicit xxxdimord for each of the known fields.
%
% (2012) Changed the dimensionality of the waveform to allow both
% stereotrode and tetrode data to be represented.
%
% (2011) Defined a consistent spike data representation that can
% also contain the Fourier spectrum and other fields. Use the xxxdimord
% to indicate the dimensions of the field.
%
% (2010) Introduced the time and the trialtime fields.
%
% (2007) Introduced the spike data representation.
%
% See also FT_DATATYPE, FT_DATATYPE_RAW, FT_DATATYPE_FREQ, FT_DATATYPE_TIMELOCK

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

% get the optional input arguments, which should be specified as key-value pairs
version = ft_getopt(varargin, 'version', 'latest');

if strcmp(version, 'latest')
  version = '2020';
end

if isempty(spike)
  return;
end

switch version
  case '2020'
    % ensure it is up to date with the previous version, which incoreporates changes from 2011 and 2012
    spike = ft_datatype_spike(spike, 'version', '2012');

    % remove the general dimord
    if isfield(spike, 'dimord')
      spike = rmfield(spike, 'dimord');
    end

    % add explicit dimords for the known fields
    if isfield(spike, 'timestamp') && ~isfield(spike, 'timestampdimord')
      spike.timestampdimord = '{chan}_spike';
    end
    if isfield(spike, 'unit') && ~isfield(spike, 'unitdimord')
      spike.unitdimord = '{chan}_spike';
    end
    if isfield(spike, 'waveform') && ~isfield(spike, 'waveformdimord')
      spike.waveformdimord = '{chan}_lead_time_spike';
    end
    if isfield(spike, 'fourierspctrm') && ~isfield(spike, 'fourierspctrmdimord')
      spike.fourierspctrmdimord = '{chan}_spike_lfplabel_freq';
    end

  case '2012'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if isfield(spike,'origtrial') && isfield(spike,'origtime')
      % this was the old spiketriggered spectrum output
      ft_warning('The spike datatype format you are using is depreciated. Converting to newer spike format');
      spike.trial = {spike.origtrial};
      spike       = rmfield(spike,'origtrial');
      spike.time  = {spike.origtime};
      spike       = rmfield(spike,'origtime');
      if ~isa(spike.fourierspctrm, 'cell')
        spike.fourierspctrm = {spike.fourierspctrm};
      end
      if ~isfield(spike, 'trialtime')
        % determine from the data itself
        ft_warning('Reconstructing the field trialtime from spike.origtime and spike.origtrial. This is not the original representation');
        tmax  = nanmax(spike.trial{1});
        tsmin = nanmin(spike.time{1});
        tsmax = nanmax(spike.time{1});
        spike.trialtime = [tsmin*ones(tmax,1) tsmax*ones(tmax,1)];
      end
      spike.lfplabel = spike.label; % in the old format, these were the lfp channels
      try
        spike.label    = spike.cfg.spikechannel;
      catch
        try
          spike.label = spike.spikechannel;
        catch
          spike.label = {'unit1'}; %default
        end
      end
      spike.dimord = '{chan}_spike_lfpchan_freq';
    end

    % fix the waveform dimensions
    if isfield(spike,'waveform')
      nUnits = length(spike.waveform);
      hasdat = false(1,nUnits);
      for iUnit = 1:nUnits
        hasdat(iUnit) = ~isempty(spike.waveform{iUnit});
      end

      if any(hasdat) %otherwise, ignore
        if ~isfield(spike, 'dimord')
          spike.dimord = '{chan}_lead_time_spike';
        end
        % fix the dimensions of the waveform dimord.
        for iUnit = 1:nUnits
          dim = size(spike.waveform{iUnit});
          if length(dim)==2 && ~isempty(spike.waveform{iUnit})
            nSpikes = length(spike.timestamp{iUnit}); % check what's the spike dimension from the timestamps
            spikedim = dim==nSpikes;
            if isempty(spikedim)
              ft_error('waveforms contains data but number of waveforms does not match number of spikes');
            end
            if spikedim==1
              spike.waveform{iUnit} = permute(spike.waveform{iUnit},[3 2 1]);
            else
              spike.waveform{iUnit} = permute(spike.waveform{iUnit},[3 1 2]);
            end

          elseif length(dim)==3 && ~isempty(spike.waveform{iUnit})
            nSpikes = length(spike.timestamp{iUnit}); % check what's the spike dimension from the timestamps
            spikedim = dim==nSpikes;
            % determine from the remaining dimensions which is the lead
            leaddim  = dim<6 & dim~=nSpikes;
            sampdim  = dim>=6 & dim~=nSpikes;
            if isempty(spikedim)
              ft_error('waveforms contains data but number of waveforms does not match number of spikes');
            end
            if sum(leaddim)~=1 || sum(sampdim)~=1, continue,end % in this case we do not know what to do
            if find(spikedim)~=3 && find(leaddim)~=1 && find(sampdim)~=2
              spike.waveform{iUnit} = permute(spike.waveform{iUnit}, [find(leaddim) find(sampdim) find(spikedim)]);
            end
          end
        end

      end

    end

    % ensure that we always deal with row vectors: for consistency of
    % representation
    if isfield(spike,'time')
      for iUnit = 1:length(spike.time)
        if size(spike.time{iUnit},2)==1
          spike.time{iUnit} = spike.time{iUnit}(:)';
        end
      end
    end

    if isfield(spike,'time')
      for iUnit = 1:length(spike.trial)
        if size(spike.trial{iUnit},2)==1
          spike.trial{iUnit} = spike.trial{iUnit}(:)';
        end
      end
    end

    if isfield(spike,'timestamp')
      for iUnit = 1:length(spike.timestamp)
        if size(spike.timestamp{iUnit},2)==1
          spike.timestamp{iUnit} = spike.timestamp{iUnit}(:)';
        end
      end
    end

    if isfield(spike,'unit')
      for iUnit = 1:length(spike.unit)
        if size(spike.unit{iUnit},2)==1
          spike.unit{iUnit} = spike.unit{iUnit}(:)';
        end
      end
    end

  otherwise
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ft_error('unsupported version "%s" for spike datatype', version);
end

spike = sortfields(spike);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s2 = sortfields(s1)
s2 = struct();
fn = sort(fieldnames(s1));
for i=1:numel(fn)
  s2.(fn{i}) = s1.(fn{i});
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%Starts from the current decoded state, takes input as minimum distance to
%reach that state and previous state And returns previous state and decoded
%information bit corresponding to that state.

function [prev_state,decoded_bit]=prev_stage(curr_state,distance_prev,metric)
    if(curr_state==1)
        if(distance_prev(1)+metric(1) <= distance_prev(3)+metric(5))
            prev_state=1;decoded_bit=0;
        else
            prev_state=3;decoded_bit=0;
        end
    end

    if(curr_state==2)
        if(distance_prev(1)+metric(2) <= distance_prev(3)+metric(6))
            prev_state=1;decoded_bit=1;
        else
            prev_state=3;decoded_bit=1;
        end
    end

    if(curr_state==3)
        if(distance_prev(2)+metric(3) <= distance_prev(4)+metric(7))
            prev_state=2;decoded_bit=0;
        else
            prev_state=4;decoded_bit=0;
        end
    end

    if(curr_state==4)
        if(distance_prev(2)+metric(4) <= distance_prev(4)+metric(8))
            prev_state=2;decoded_bit=1;
        else
            prev_state=4;decoded_bit=1;
        end
    end

end


text:
function out = matrix(L, varargin)
%   MATRIX(L, DIM) returns an DIM by (DIM+D) discretization of the linear
%   chebop L, where D is the differential order of L. If L is not linear
%   then an error is thrown.
%
%   If L has assigned boundary conditions then these are included in the
%   discretization, increasing the number of rows. If the domain of L has
%   breakpoints then DIM must be a vector specifying the discretization size
%   on each interval. In this case, appropriate continuity constraints are
%   included in the discretization, further increasing the number of rows
%   (by D times the number of breakpoints).
%
%   MATRIX(L, DIM, PREFS) allows additional preferences to be passed
%   via the CHEBOPPREF PREFS.
%
%   MATRIX(L, DIM, 'oldschool') forces the returned differentiation
%   matrices to be square rather than rectangular. See LINOP/MATRIX and
%   LINOP/FEVAL for further details.
%
%   MATRIX(L, DIM, 'oldschool', PREFS) allows, again, additional
%   preferences to be passed via the CHEBOPPREF PREFS.
%
%   For mathematical details of these matrices see Aurentz and Trefethen,
%   "Block operators and spectral discretizations," SIAM Review, 2017.
%   For the 'oldschool' discretizations, see Trefethen, Spectral Methods
%   in MATLAB, SIAM, 2000.
%
% See also FEVAL, LINOP/MATRIX, LINOP/FEVAL.

% Get the preferences if given:
isPrefGiven = 0;
for j = 1:nargin-1
    item = varargin{j};
    if ( isa(item,'cheboppref') )
        prefs = item;
        isPrefGiven = 1;
    end
end

% Otherwise, use default prefs
if ( ~isPrefGiven )
    prefs = cheboppref();
end

% Linearize:
[L_linop, ~, fail] = linop(L);
if ( fail )
    error('CHEBFUN:CHEBOP:matrix:nonlinear',...
        'Matrix expansion is only allowed for linear CHEBOP objects.')
end

% Determine the discretization:
prefs = determineDiscretization(L, L_linop, prefs);

% Derive continuity conditions in the case of breakpoints:
L_linop = deriveContinuity(L_linop);

% Call LINOP/FEVAL or LINOP/MATRIX:
if ( (numel(varargin) > 1) && strcmpi(varargin{2}, 'oldschool') )
    warnState = warning('off', 'CHEBFUN:LINOP:feval:deprecated');
    out = feval(L_linop, varargin{:});
    warning(warnState);
else
    % Add the preferences in vargarin to pass them to LINOP/MATRIX:
    if ( isPrefGiven )
        % If a CHEBOPPREF was passed to the method, it will have been at the
        % last position of varargin, indexed at nargin-1. Overwrite it with the
        % current PREFS.DISCRETIZATION, as the discretization might have changed
        % in the periodic case:
        varargin{nargin-1} = prefs;
    else
        % Otherwise, add the PREFS.DISCRETIZATION to VARARGIN, so that it can be
        % passed to the call to LINOP/MATRIX below.
        varargin{nargin} = prefs;
    end
    out = matrix(L_linop, varargin{:});
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


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

    enumeration
        Radial('Radial Vector');
        Velocity('Velocity Vector');
        AngularMomentum('Angular Momentum Vector');
        North('North Vector');
        East('East Vector');
        AscNode('Ascending Node Vector');
        Periapsis('Periapsis Vector');
    end

    properties
        name(1,:) char
    end

    methods
        function obj = VehicleStateVectorTypeEnum(name)
            obj.name = name;
        end
    end

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

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

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

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function drawSeg(h,x,c)

% DRAWSEG  Draw 2D or 3D segment with color.
%   DRAWSEG(H,S) Draws 2D or 3D segment S=[p1;p2] with handle H.
%
%   DRAWSEG(H,S,C) uses color C.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

if numel(x) == 4

    set(h,'xdata',x([1 3]),'ydata',x([2 4]),'vis','on')

elseif numel(x) == 6

    set(h,'xdata',x([1 4]),'ydata',x([2 5]),'zdata',x([3 6]),'vis','on')

else

    error('??? Size of vector ''x'' not correct.')

end

if nargin == 3
    set(h,'color',c)
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
% statcondfiledtrip()  - same as statcond except that it uses the fieldtrip
%                        statistical functions. This is useful to perform
%                        a wider variety of corrections for multiple
%                        comparisons for instance.
% Usage:
%          >> [stats, df, pvals, surrog] = statcond( data, 'key','val'... );
% Inputs:
%   data       = same as for statcond()
%
% Optional inputs:
%   'paired'   = ['on'|'off'] pair the data array {default: 'on' unless
%                the last dimension of data array is of different lengths}.
%   'mode'     = ['perm'|'param'] mode for computing the p-values:
%                 'param' = parametric testing (standard ANOVA or t-test);
%                 'perm' = non-parametric testing using surrogate data
%                  made by permuting the input data {default: 'perm'}
%   'naccu'    = this input is passed on as 'numrandomization' to Fieldtrip
%   'chanlocs' = EEGLAB channel location structure to perfom statistics
%                and cluster correction for multiple comparisons accross
%                channels (native fieldtrip entries may also be used).
%
% Fieldtrip options:
%   Any option to the freqanalysis, the statistics_montecarlo, the
%   statistics_analysis, statistics_stat, statistics_glm may be used
%   using 'key', val argument pairs.
%
% Outputs:
%   stats      = F- or T-value array of the same size as input data without
%                the last dimension. A T value is returned only when the data
%                includes exactly two conditions.
%   df         = degrees of freedom, a (2,1) vector, when F-values are returned
%   pvals      = array of p-values. Same size as input data without the last
%                data dimension. All returned p-values are two-tailed.
%   surrog     = surrogate data array (same size as input data with the last
%                dimension filled with a number ('naccu') of surrogate data sets.
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2005-
%         With thanks to Robert Oostenveld for fruitful discussions
%         and advice on this function.
%
% See also: freqanalysis(), statistics_montecarlol()

% Copyright (C) Arnaud Delorme
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [ ori_vals, df, pvals, surrogval ] = statcondfieldtrip( data, varargin );

    if nargin < 1
        help statcondfieldtrip;
        return;
    end;

    [g cfgparams] = finputcheck( varargin, { 'naccu'    'integer'   [1 Inf]             200;
                                             'mode'     'string'    { 'param' 'perm' 'bootstrap' }  'param';
                                             'chanlocs' 'struct'    { }   struct([]);
                                             'method'   'string'    { 'montecarlo' 'analytic' 'stat' 'glm' } 'analytic';
                                             'paired'   'string'    { 'on' 'off' }      'on' }, 'statcond', 'ignore');
    if isstr(g), error(g); end;

    if size(data,2) == 1, data  = transpose(data); end; % cell array transpose

    tmpsize   = size(data{1});

    % cfg configuration for Fieldtrip
    % -------------------------------
    cfg = struct(cfgparams{:});
    cfg.method      = g.method;
    if ~strcmpi(g.mode, 'param'), cfg.method = 'montecarlo'; end;
    if ~isempty(g.chanlocs)
        EEG = eeg_emptyset;
        EEG.chanlocs = g.chanlocs;
        EEG.nbchan   = length(g.chanlocs);
        EEG.data     = zeros(EEG.nbchan,100,1);
        EEG = eeg_checkset(EEG);
        tmpcfg = eeglab2fieldtrip(EEG, 'preprocessing', 'none');
        lay = prepare_layout(tmpcfg, tmpcfg);
        tmpcfg.layout        = lay;
        tmpcfg.neighbourdist = 30; %mean(abs([ EEG.chanlocs.X ]));
        cfg.neighbours    = neighbourselection(tmpcfg, []);
    else
        cfg.neighbours  = {};
    end;
    cfg.feedback    = 'no';
    cfg.ivar        = 1;
    cfg.alpha       = 0.05;
    cfg.numrandomization = g.naccu;

    % test if data can be paired
    % --------------------------
    if length(unique(cellfun('size', data, ndims(data{1}) ))) > 1
        g.paired = 'off';
    end;
    fprintf('%d x %d, ', size(data,1), size(data,2));
    if strcmpi(g.paired, 'on')
         fprintf('paired data, ');
    else fprintf('unpaired data, ');
    end;
    if size(data,1) == 1 & size(data,2) == 2
         fprintf('computing T values\n');
    else fprintf('computing F values\n');
    end;

    % output text
    % -----------
    if ~strcmpi(g.mode, 'param')
        fprintf('Accumulating (of %d):', g.naccu);
    end;

    if size(data,1) == 1, % only one row

        if size(data,2) == 2 & strcmpi(g.paired, 'on')

            % paired t-test (very fast)
            % -------------
            cfg.statistic   = 'depsamplesT';
            [newdata design1 design2 design3] = makefieldtripdata(data, 0, g.chanlocs);
            cfg.design      = [ design1; design3 ];
            cfg.uvar        = 2;

            cfg
            stat            = freqstatistics(cfg, newdata{:});
            ori_vals        = stat.stat;
            df              = [];
            pvals           = stat.prob;
            return;

        elseif size(data,2) == 2 & strcmpi(g.paired, 'off')

            % paired t-test (very fast)
            % -------------
            cfg.statistic   = 'indepsamplesT';
            [newdata design1] = makefieldtripdata(data, 0, g.chanlocs);
            cfg.design      = design1;
            stat            = freqstatistics(cfg, newdata{:});
            ori_vals        = stat.stat;
            df              = [];
            pvals           = stat.prob;
            return;

        elseif strcmpi(g.paired, 'on')

            % one-way ANOVA (paired) this is equivalent to unpaired t-test
            % -------------
            cfg.statistic   = 'depsamplesF';
            [newdata design1 design2 design3] = makefieldtripdata(data, 0, g.chanlocs);
            cfg.design      = [ design1; design3 ];
            cfg.uvar        = 2;
            stat            = freqstatistics(cfg, newdata{:});
            ori_vals        = stat.stat;
            df              = [];
            pvals           = stat.prob;
            return;

        else
            % one-way ANOVA (unpaired)
            % -------------
            cfg.statistic   = 'indepsamplesF';
            [newdata design1] = makefieldtripdata(data, 0, g.chanlocs);
            cfg.design      = [ design1 ];
            stat            = freqstatistics(cfg, newdata{:});
            ori_vals        = stat.stat;
            df              = [];
            pvals           = stat.prob;
            return;

        end;
    else
        if strcmpi(g.paired, 'on')

            % two-way ANOVA (paired)
            % -------------
            cfg.statistic   = 'anovan';
            [newdata design1 design2 design3] = makefieldtripdata(data, 0, g.chanlocs);
            cfg.design      = [ design1; design2; design3 ];
            cfg.effect      = 'X1*X2';
            cfg.ivar        = [1 2];
            cfg.uvar        = 3;
            stat            = freqstatistics(cfg, newdata{:});
            ori_vals        = stat.stat;
            df              = [];
            pvals           = stat.prob;
            return;

        else

            % two-way ANOVA (unpaired)
            % -------------
            cfg.statistic   = 'anovan';
            cfg.clustercritval = 4.5416; % 95 percentile of n =10000; a = { rand(n,10) rand(n,10); rand(n,10) rand(n,10) }; [F df p ] = statcondfieldtrip(a, 'paired', 'off');
            [newdata design1 design2] = makefieldtripdata(data, 0, g.chanlocs);
            if ~isempty(g.chanlocs)
                for index = 1:length(newdata)
                    newdata{index}.powspctrm = squeeze(newdata{index}.powspctrm);
                    newdata{index}.label     = { g.chanlocs.labels };
                    newdata{index}.freq      = 1;
                end;
            end;
            cfg
            newdata{1}
            cfg.design      = [ design1; design2 ];
            cfg.effect      = 'X1*X2';
            cfg.ivar        = [1 2];
            stat            = freqstatistics(cfg, newdata{:});
            ori_vals        = stat.stat;
            df              = [];
            pvals           = stat.prob;
            return;

        end;
    end;

function val = myndims(a)
    if ndims(a) > 2
        val = ndims(a);
    else
        if size(a,1) == 1,
            val = 2;
        elseif size(a,2) == 1,
            val = 1;
        else
            val = 2;
        end;
    end;

function [newdata, design1, design2, design3] = makefieldtripdata(data, chanflag, chanlocs);

    newdata = {};
    for i = 1:length(data(:))

        switch myndims(data{1})
          case 1,
            newdata{i}.powspctrm = data{i};

          case 2,
            if chanflag
                 newdata{i}.powspctrm = transpose(data{i});
            else newdata{i}.powspctrm = reshape(transpose(data{i}), size(data{i},2), 1, size(data{i},1));
            end;

          case 3,
            if chanflag
                 newdata{i}.powspctrm = permute(data{i}, [3 1 2]);
            else newdata{i}.powspctrm = permute(data{i}, [3 4 1 2]);
            end;

          case 4,
            if chanflag
                 newdata{i}.powspctrm = permute(data{i}, [4 1 2 3]);
            else newdata{i}.powspctrm = permute(data{i}, [3 5 1 2]);
            end;

        end;

        newdata{i}.dimord    = 'rpt_chan_freq_time';
        newdata{i}.label     = cell(1,size(newdata{i}.powspctrm,2));
        newdata{i}.label(:)  = { 'cz' };
        newdata{i}.freq      = [1:size(newdata{i}.powspctrm,3)];
        newdata{i}.time      = [1:size(newdata{i}.powspctrm,4)];

        % below in case channels are specified
        % not that statistics are done on time x frequencies or channels
        % so time x frequency x channels do not work yet here
        if ~isempty(chanlocs)
            newdata{i}.powspctrm = squeeze(newdata{i}.powspctrm);
            newdata{i}.label     = { chanlocs.labels };
            newdata{i}.freq      = 1;
            newdata{i}.time      = 1;
        end;

    end;

    design1 = [];
    design2 = [];
    design3 = [];
    for i = 1:size(data,2)
        for j = 1:size(data,1)
            ij = j+(i-1)*size(data,1);
            design1 = [ design1 ones(1, size(newdata{i}.powspctrm,1))*i ];
            design2 = [ design2 ones(1, size(newdata{i}.powspctrm,1))*j ];
            design3 = [ design3 [1:size(newdata{i}.powspctrm,1)] ];
        end;
    end;

author: PatternRecognition

repo: OpenBMI


text:
function p = part_table ( n )

%*****************************************************************************80
%
%% PART_TABLE tabulates the number of partitions of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the integer to be partitioned.
%    N must be positive.
%
%    Output, integer P(1:N+1), P(I+1) is the number of partitions of I.
%
  offset = 1;

  p = zeros ( n + 1, 1 );

  p(0+offset) = 1;
  p(1+offset) = 1;

  for i = 2 : n

    sign = 1;
    psum = 0;
    w = 1;
    j = 1;
    wprime = w + j;

    while ( w < n )

      if ( 0 <= i - w )
        if ( sign == 1 )
          psum = psum + p(i-w+offset);
        else
          psum = psum - p(i-w+offset);
        end
      end

      if ( wprime <= i )

        if ( sign == 1 )
          psum = psum + p(i-wprime+offset);
        else
          psum = psum - p(i-wprime+offset);
        end

      end

      w = w + 3 * j + 1;
      j = j + 1;
      wprime = w + j;
      sign = - sign;

    end

    p(i+offset) = psum;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Reference Frame Alignment
%
%%
% The most important difference between MTEX and many other EBSD software
% is that in MTEX the Euler angle reference is always the map reference
% frame. This mean the $x$ and $z$ axes of the map are exactly the rotation
% axes of the Euler angles.
%
% In case the map coordinates and the Euler angles in your data are with
% respect to different reference frames it is highly recommended to correct
% for this while importing the data into MTEX. This section explains in
% detail how to do this.
%
%% On Sreen Orientation of the EBSD Map
%
% Many people are concerned when the images produced by MTEX are not
% aligned exactly as they are in their commercial software. It is indeed
% very important to understand exactly the alignment of your data. However,
% the important point is not whether a map is upside down on your screen or
% not. The important point is how your map aligns with the specimen, as we
% want to use the map to describe properties of the specimen.
%
% There are basically two components in an EBSD data set that refer to the
% specimen reference frame: the spatial coordinates $x$, $y$ and the Euler
% angles $\phi_1$, $\Phi$, $\phi_2$. To explain the difference have a look
% at the EDAX export dialog
%
% <<edax_coordinate_systems.png>>
%
% Here we have the axes $x$ and $y$ which describe how the map coordinates
% needs to be interpreted and the axes $A_1$, $A_2$, $A_3$ which describe
% how the Euler angles, and in consequence, the pole figures needs to be
% interpreted. We see that in none of these settings the map reference
% system coincides with the Euler angle reference frame.
%
% This situation is not specific to EDAX but occurs as well with EBSD data
% from Oxford or Bruker, all of them using different reference system
% alignments. For that reason MTEX strongly recommends to transform the data
% such that both map coordinates and Euler angles refer to the same
% coordinate system.
%
% Doing this we have two choices:
%
% # transform everything to the reference system $x$, $y$ using the option
% |'convertEuler2SpatialReferenceFrame'|. This will keep the map
% coordinates while changing the Euler angles
% # transform everything to the reference system $A_1$, $A_2$, $A_3$ using
% the option |'convertSpatial2EulerReferenceFrame'|. This will keep the
% Euler angles while changing the map coordinates.
%
% In the case of EDAX data imported from an |*.ang| file we still need to
% specify the export option used by the EDAX software. This is done by the
% options |'setting 1'|, |'setting 2'|, |'setting 3'| or |'setting 4'|.
%
% Since setting 2 is default for most EDAX exports a typical command for
% importing data from an ang file would look like this

ebsd = EBSD.load([mtexEBSDPath filesep 'olivineopticalmap.ang'],...
  'convertEuler2SpatialReferenceFrame','setting 2')

plot(ebsd('olivine'),ebsd('olivine').orientations,'coordinates','on')

%%
% The plot does not yet fit the alignment of the map in the EDAX software
% as it plots the x-axis be default to east and the z-axis into the plane.
% This is only a plotting convention and can be set in MTEX by

plotx2east
plotzIntoPlane

plot(ebsd('olivine'),ebsd('olivine').orientations,'coordinates','on')

%%
% Other plotting conventions are |plotx2north|, |plotx2west|, |plotx2south|
% and |plotzOutOfPlane|. Note that these options only alter the orientation
% of the EBSD map and the pole figures on the screen but does not change
% any data.
%
%% Verify the reference system
% One way of verifying the reference systems is to visualize crystal shapes
% on top of the orientation map. To do this we proceed as follows

% reconstruct grains
grains = calcGrains(ebsd('indexed'));

% chose the correct crystal shape (cubic, hex are generic forms)
cS = crystalShape.olivine;

% select only large grains
largeGrains = grains(grains.grainSize>500)

% and plot the crystal shapes
hold on
plot(largeGrains,cS)
hold off

%%
% It may also be helpful to inspect pole figures

h = Miller({1,0,0},{0,1,0},{0,0,1},ebsd('O').CS);
plotPDF(ebsd('O').orientations,h,'contourf')

%%
% As pole figures display data relative to the specimen reference frame
% MTEX automatically aligns them on the screen exactly as the spatial map
% above, i.e., according to our last definition with x pointing towards
% east and y to the south.
%
%% Change the map reference system
% In order to manually change the map reference frame one may apply a
% rotation to the map coordinates only. E.g. to flip the map left to right
% while preserving the Euler angles one can do

rot = rotation.byAxisAngle(yvector,180*degree);
ebsd_rot = rotate(ebsd,rot,'keepEuler');

% reconstruct grains
grains = calcGrains(ebsd_rot('indexed'));

% select only large grains
largeGrains = grains(grains.grainSize>500);

plot(ebsd_rot('olivine'),ebsd_rot('olivine').orientations,'coordinates','on')

% and plot the crystal shapes
hold on
plot(largeGrains,cS)
hold off

%% Change the Euler angle reference system
% Analogously we may change the Euler angle reference frame while keeping
% the map coordinates

ebsd_rot = rotate(ebsd,rot,'keepXY');

% reconstruct grains
grains = calcGrains(ebsd_rot('indexed'));

% select only large grains
largeGrains = grains(grains.grainSize>500);

plot(ebsd_rot('olivine'),ebsd_rot('olivine').orientations,'coordinates','on')

% and plot the crystal shapes
hold on
plot(largeGrains,cS)
hold off

%% Changing both reference system simultaneously
%
% Sometimes it is necessary to relate the EBSD data to a different external
% reference frame, or to  change the external reference frame from one to
% the other, e.g. if one wants to concatenate several ebsd data sets where
% the mounting was not done in perfect coincidence. In these cases the data
% has to be rotated or shifted by the commands <EBSD.rotate.html rotate>
% and <EBSD.shift.html shift>. The following commands rotate both reference
% frames of the entire data set by 5 degree about the z-axis.

% define a rotation
rot = rotation.byAxisAngle(zvector,5*degree);

% rotate the EBSD data
ebsd_rot = rotate(ebsd,rot);

% reconstruct grains
grains = calcGrains(ebsd_rot('indexed'));

% select only large grains
largeGrains = grains(grains.grainSize>500);

plot(ebsd_rot('olivine'),ebsd_rot('olivine').orientations,'coordinates','on')

% and plot the crystal shapes
hold on
plot(largeGrains,cS)
hold off

author: mtex-toolbox

repo: mtex

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


text:
% std_editset() - modify a STUDY set structure.
%
% Usage:
%             >> [STUDY, ALLEEG] = std_editset(STUDY, ALLEEG, key1, val1, ...);
% Inputs:
%   STUDY      - EEGLAB STUDY set
%   ALLEEG     - vector of the EEG datasets included in the STUDY structure
%
% Optional inputs:
%   'commands' - {cell_array} change STUDY (see command description and
%                 example below.
%   'name'     - [string] specify a (mnemonic) name for the STUDY structure.
%                {default: ''}
%   'task'     - [string] attach a description of the experimental task(s)
%                performed by the STUDY subjects {default: ''}.
%   'filename' - [string] filename for the STUDY set.
%   'filepath' - [string] file path (directory/folder) in which the STUDY file
%                will be saved.
%   'addchannellabels' - ['on'|'off'] add channel labels ('1', '2', '3', ...)
%                to all datasets of a STUDY to ensure that all STUDY functions
%                will work {default: 'off'}
%   'notes'    - [string] notes about the experiment, the datasets, the STUDY,
%                or anything else to store with the STUDY itself {default: ''}.
%   'updatedat' - ['on'|'off'] update 'subject' 'session' 'condition' and/or
%                'group' fields of STUDY dataset(s).
%   'savedat'   - ['on'|'off'] re-save datasets
%   'inbrain'   - ['on'|'off'] select components for clustering from all STUDY
%                 datasets with equivalent dipoles located inside the brain volume.
%                 Dipoles are selected based on their residual variance and their
%                 location {default: 'off'}
%
% Each of the 'commands' (above) is a cell array composed of any of the following:
%   'index'     - [integer] modify dataset index.
%   'remove'    - [integer] remove dataset index.
%   'subject'   - [string] subject code.
%   'condition' - [string] dataset condition.
%   'session '  - [integer] dataset session number.
%   'group'     - [string] dataset group.
%   'load'      - [filename] load dataset from specified filename
%   'dipselect' - [float<1] select components for clustering from all STUDY
%                 datasets with dipole model residual var. below this value.
%   'inbrain'   - ['on'|'off'] same as above. This option may also be
%                 placed in the command list (preceeding the 'dipselect'
%                 option).
%
% Outputs:
%   STUDY      - a new STUDY set containing some or all of the datasets in ALLEEG,
%                plus additional information from the optional inputs above.
%   ALLEEG     - a vector of EEG datasets included in the STUDY structure
%
%  See also:  pop_createstudy(), std_loadalleeg(), pop_clust(), pop_preclust(),
%             eeg_preclust(), eeg_createdata()
%
% Authors: Arnaud Delorme, Hilit Serby, SCCN/INC/UCSD, October , 2004-

% Copyright (C) Arnaud Delorme & Scott Makeig, SCCN/INC/UCSD, October 11, 2004, smakeig@ucsd.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [STUDY, ALLEEG] = std_editset(STUDY, ALLEEG, varargin)

if (nargin < 3)
    help std_editset;
    return;
end;

% decode input parameters
% -----------------------
g = finputcheck(varargin, { 'updatedat' 'string'  { 'on' 'off' }  'off';
                            'name'      'string'  { }             '';
                            'task'      'string'  { }             '';
                            'notes'     'string'  { }             '';
                            'filename'  'string'  { }             '';
                            'filepath'  'string'  { }             '';
                            'resave'    'string'  { 'on' 'off' 'info' }  'off';
                            'savedat'   'string'  { 'on' 'off' }  'off';
                            'addchannellabels' 'string'  { 'on' 'off' }  'off';
                            'rmclust'   'string'  { 'on' 'off' }  'on';
                            'inbrain'   'string'  { 'on' 'off' }  'off';
                            'commands'  'cell'    {}              {} }, 'std_editset');
if isstr(g), error(g); end;

if ~isempty(g.name),  STUDY.name  = g.name; end
if ~isempty(g.task),  STUDY.task  = g.task; end
if ~isempty(g.notes), STUDY.notes = g.notes; end

% make one cell array with commands
% ---------------------------------
allcoms = {};
if ~isempty(g.commands)
    if iscell(g.commands{1})
        for k = 1:length(g.commands)
            % put index field first
            indindex   = strmatch('index', lower(g.commands{k}(1:2:end)));
            if ~isempty(indindex)
                 tmpcom = { 'index' g.commands{k}{2*(indindex-1)+1+1} g.commands{k}{:} };
            else tmpcom = g.commands{k};
            end;
            allcoms = { allcoms{:} tmpcom{:} };
        end;
    else
        allcoms = g.commands;
    end;
end;
g.commands = allcoms;

% add 'dipselect' command if 'inbrain' option is selected
% ---------------------------------
dipselectExists = false;
for k = 1:2:length(g.commands)
    if strcmp(g.commands{k},'dipselect')
       dipselectExists = true;
    end;
end;
if strcmp(g.inbrain,'on') && ~dipselectExists
    g.commands{length(g.commands)+1} = 'dipselect';
    g.commands{length(g.commands)+1} = 0.15;
end;

% copy values
% -----------
if ~isfield(STUDY, 'datasetinfo')
    for realindex = 1:length(ALLEEG)
        if ~isempty(ALLEEG(realindex).data)
            [tmppath tmpfile tmpext] = fileparts(  fullfile(ALLEEG(realindex).filepath, ALLEEG(realindex).filename) );
            STUDY.datasetinfo(realindex).filepath  = tmppath;
            STUDY.datasetinfo(realindex).filename  = [ tmpfile tmpext ];
            STUDY.datasetinfo(realindex).subject   = ALLEEG(realindex).subject;
            STUDY.datasetinfo(realindex).session   = ALLEEG(realindex).session;
            STUDY.datasetinfo(realindex).condition = ALLEEG(realindex).condition;
            STUDY.datasetinfo(realindex).group     = ALLEEG(realindex).group;
        end;
    end;
end;

% execute commands
% ----------------
currentind = 1;
rmlist = [];
for k = 1:2:length(g.commands)

    switch g.commands{k}
        case 'index'
            currentind = g.commands{k+1};
        case 'subject'
            STUDY.datasetinfo(currentind).subject = g.commands{k+1};
        case 'comps'
            STUDY.datasetinfo(currentind).comps = g.commands{k+1};
        case 'condition'
            STUDY.datasetinfo(currentind).condition = g.commands{k+1};
        case 'group'
            STUDY.datasetinfo(currentind).group   = g.commands{k+1};
        case 'session'
            STUDY.datasetinfo(currentind).session = g.commands{k+1};
        case 'session'
            STUDY.datasetinfo(currentind).session = g.commands{k+1};
        case 'remove'
            ALLEEG(end+2)            = ALLEEG(end);
            ALLEEG(g.commands{k+1})  = ALLEEG(end-1); % empty dataset
            ALLEEG(end-1:end)        = [];
            STUDY.datasetinfo(end+2)           = STUDY.datasetinfo(end);
            STUDY.datasetinfo(g.commands{k+1}) = STUDY.datasetinfo(end-1);
            STUDY.datasetinfo(end-1:end)       = [];
            if isfield(STUDY.datasetinfo, 'index')
                STUDY.datasetinfo = rmfield(STUDY.datasetinfo, 'index');
            end;
            STUDY.datasetinfo(1).index = [];
            STUDY.changrp = [];
        case 'return', return;
        case 'inbrain'
            g.inbrain = g.commands{k+1};
        case 'dipselect'
            STUDY = std_checkset(STUDY, ALLEEG);
            rv = g.commands{k+1};
            clusters = std_findsameica(ALLEEG);

            for cc = 1:length(clusters)

                idat = 0;
                for tmpi = 1:length(clusters{cc})
                    if isfield(ALLEEG(clusters{cc}(tmpi)).dipfit, 'model')
                        idat = clusters{cc}(tmpi);
                    end;
                end;

                indleft = [];
                if rv ~= 1
                    if idat ~= 0
                        if strcmp(g.inbrain,'on')
                            fprintf('Selecting dipoles with less than %%%2.1f residual variance and removing dipoles outside brain volume in dataset ''%s''\n', ...
                                100*rv, ALLEEG(idat).setname);
                            indleft = eeg_dipselect(ALLEEG(idat), rv*100,'inbrain');
                        else
                           fprintf('Selecting dipoles with less than %%%2.1f residual variance in dataset ''%s''\n', ...
                                100*rv, ALLEEG(idat).setname);
                            indleft = eeg_dipselect(ALLEEG(idat), rv*100,'rv');
                        end;
                    else
                        fprintf('No dipole information found in ''%s'' dataset, using all components\n', ALLEEG.setname)
                    end
                end;
                for tmpi = 1:length(clusters{cc})
                    STUDY.datasetinfo(clusters{cc}(tmpi)).comps = indleft;
                end;
            end;
            STUDY.cluster = [];
            STUDY = std_checkset(STUDY, ALLEEG); % recreate parent dataset

        case 'load'
            TMPEEG = std_loadalleeg( { g.commands{k+1} } );
            ALLEEG = eeg_store(ALLEEG, eeg_checkset(TMPEEG), currentind);
            ALLEEG(currentind).saved = 'yes';

            % update datasetinfo structure
            % ----------------------------
            [tmppath tmpfile tmpext] = fileparts( fullfile(ALLEEG(currentind).filepath, ...
                                                           ALLEEG(currentind).filename) );
            STUDY.datasetinfo(currentind).filepath  = tmppath;
            STUDY.datasetinfo(currentind).filename  = [ tmpfile tmpext ];
            STUDY.datasetinfo(currentind).subject   = ALLEEG(currentind).subject;
            STUDY.datasetinfo(currentind).session   = ALLEEG(currentind).session;
            STUDY.datasetinfo(currentind).condition = ALLEEG(currentind).condition;
            STUDY.datasetinfo(currentind).group     = ALLEEG(currentind).group;
            STUDY.datasetinfo(currentind).index     = currentind;
        otherwise, error(sprintf('Unknown command %s', g.commands{k}));
    end
end

% add channel labels automatically
% -------------------------------
if strcmpi(g.addchannellabels, 'on')
    disp('Generating channel labels for all datasets...');
    for currentind = 1:length(ALLEEG)
        for ind = 1:ALLEEG(currentind).nbchan
            ALLEEG(currentind).chanlocs(ind).labels = int2str(ind);
        end;
    end;
    ALLEEG(currentind).saved = 'no';
    g.savedat = 'on';
end;

% update ALLEEG structure?
% ------------------------
if strcmpi(g.updatedat, 'on')
    for currentind = 1:length(ALLEEG)
        if ~strcmpi(ALLEEG(currentind).subject,   STUDY.datasetinfo(currentind).subject)
            ALLEEG(currentind).subject          = STUDY.datasetinfo(currentind).subject;
            ALLEEG(currentind).saved            = 'no';
        end;
        if ~strcmpi(ALLEEG(currentind).condition, STUDY.datasetinfo(currentind).condition)
            ALLEEG(currentind).condition        = STUDY.datasetinfo(currentind).condition;
            ALLEEG(currentind).saved            = 'no';
        end;
        if ~isequal(ALLEEG(currentind).session, STUDY.datasetinfo(currentind).session)
            ALLEEG(currentind).session          = STUDY.datasetinfo(currentind).session;
            ALLEEG(currentind).saved            = 'no';
        end;
        if ~strcmpi(char(ALLEEG(currentind).group), char(STUDY.datasetinfo(currentind).group))
            ALLEEG(currentind).group            = STUDY.datasetinfo(currentind).group;
            ALLEEG(currentind).saved            = 'no';
        end;
    end;
end;

% remove empty datasets (cannot be done above because some empty datasets
% might not have been removed)
% ---------------------
[ ALLEEG STUDY.datasetinfo ] = removeempty(ALLEEG, STUDY.datasetinfo);

% save datasets if necessary
% --------------------------
if strcmpi(g.savedat, 'on')
    for index = 1:length(ALLEEG)
        if isempty(ALLEEG(index).filename)
            fprintf('Cannot resave ALLEEG(%d) because the dataset has no filename\n', index);
        else
            TMP = pop_saveset(ALLEEG(index), 'savemode', 'resave');
            ALLEEG = eeg_store(ALLEEG, TMP, index);
            ALLEEG(index).saved = 'yes';
        end;
    end;
end;

% remove cluster information if necessary
% ---------------------------------------
if strcmpi(g.rmclust, 'on')
    STUDY.cluster = [];
end;

% save study if necessary
% -----------------------
STUDY.changrp = [];
[STUDY ALLEEG] = std_checkset(STUDY, ALLEEG);
if ~isempty(g.filename),
    [STUDY.filepath STUDY.filename ext] = fileparts(fullfile( g.filepath, g.filename ));
    STUDY.filename = [ STUDY.filename ext ];
    g.resave = 'on';
end
if strcmpi(g.resave, 'on')
    STUDY = pop_savestudy(STUDY, ALLEEG, 'savemode', 'resave');
end;

% ---------------------
% remove empty elements
% ---------------------
function [ALLEEG, datasetinfo] = removeempty(ALLEEG, datasetinfo);

    rmindex = [];
    for index = 1:length(datasetinfo)
        if isempty(datasetinfo(index).subject) & isempty(ALLEEG(index).nbchan)
            rmindex = [ rmindex index ];
        end;
    end;
    datasetinfo(rmindex) = [];
    ALLEEG(rmindex)      = [];
    for index = 1:length(datasetinfo)
        datasetinfo(index).index = index;
    end;

    % remove empty ALLEEG structures
    % ------------------------------
    while length(ALLEEG) > length(datasetinfo)
       ALLEEG(end) = [];
    end;


text:
function zL = extrap(z, limtype);	% Slope-limited extrapolation (MUSCL)
					% Staggered scheme
% Function called: limiter

% Numbering for pressure nodes:
%   1 1 2 2                     J J
% |-o-|-o-|-o-|-o-|-o-|-o-|-o-|-o-|
% Numbering for momentum nodes:
% 1 1 2 2                        J-1  J
% |-o-|-o-|-o-|-o-|-o-|-o-|-o-|---o---|

J = length(z);
zL(1) = z(1);
for j = 2:J-1
  b     = (z(j+1) - z(j)); a = (z(j) - z(j-1));
  zL(j) = z(j) + 0.5*limiter(a, b, limtype)*(z(j+1)-z(j));
end
zL(J) = z(J);

repo: Aero-matlab


text:
function sys = IntDelayOp(m)
%
% INTDELAYOP: Integer Delay Operator. IntDelayOp delays the i-th channel of
% a digital LTI system by m(i) samples.
%
% Usage:    sys = IntDelayOp(m)
%
% INPUT
%       m is a vector of delays, m = (m0, m1, m2,..., mN)
%
% OUTPUT
%       sys: the delay operator in system format

den0 = [1 zeros(1,m(1)+1)];
num0 = [1 0];
sys = tf(num0, den0, -1);

for i = 2:length(m)

    deni = [1 zeros(1,m(i)+1)];
    numi = [1 0];
    sysi = tf(numi, deni, -1);

    sys = append(sys, sysi);
end


text:
function pass = test_battery( prefs )
% Check that chebop2 is working by using a battery of Laplace problems.
% Alex Townsend, March 2013.

if ( nargin < 1 )
    prefs = chebfunpref();
end
tol = 100*prefs.techPrefs.chebfuneps;

% Harmonic solution to the Laplace equation
N = chebop2(@(u) diff(u,2,1) + diff(u,2,2));
bdy = @(x,y) real(exp(x+1i*y));
N.lbc = @(y) bdy(-1,y); N.rbc = @(y) bdy(1,y);
N.dbc = @(x) bdy(x,-1); N.ubc = @(x) bdy(x,1);
u = N \ 0; exact = chebfun2(bdy);
pass(1) = ( norm( exact - u ) < tol );

% Harmonic solution to the Laplace equation
N = chebop2(@(u) diff(u,2,1) + diff(u,2,2));
bdy = @(x,y) real(exp(2*(x+1i*y)));
N.lbc = @(y) bdy(-1,y); N.rbc = @(y) bdy(1,y);
N.dbc = @(x) bdy(x,-1); N.ubc = @(x) bdy(x,1);
u = N \ 0; exact = chebfun2(bdy);
pass(2) = ( norm( exact - u ) < tol );

% Harmonic solution to the Laplace equation
N = chebop2(@(u) diff(u,2,1) + diff(u,2,2));
bdy = @(x,y) 10*real(exp(2*(x+1i*y)));
N.lbc = @(y) bdy(-1,y); N.rbc = @(y) bdy(1,y);
N.dbc = @(x) bdy(x,-1); N.ubc = @(x) bdy(x,1);
u = N \ 0; exact = chebfun2(bdy);
pass(3) = ( norm( exact - u ) < 10*tol );

% Harmonic solution to the Laplace equation
N = chebop2(@(u) diff(u,2,1) + diff(u,2,2));
bdy = @(x,y) real((x+1i*y).^2);
N.lbc = @(y) bdy(-1,y); N.rbc = @(y) bdy(1,y);
N.dbc = @(x) bdy(x,-1); N.ubc = @(x) bdy(x,1);
u = N \ 0; exact = chebfun2(bdy);
pass(4) = ( norm( exact - u ) < tol );

% Linearity check;
M = N;

N.lbc = @(x) (1+x).*(1-x);
N.rbc = 0; N.ubc = 0; N.dbc = 0;
u1 = N \ 0;

N = M; N.rbc = @(x) (1+x).*(1-x);
N.lbc = 0; N.ubc = 0; N.dbc = 0;
u2 = N \ 0;

N = M; N.ubc = @(x) (1+x).*(1-x);
N.rbc = 0; N.lbc = 0; N.dbc = 0;
u3 = N \ 0;

N = M; N.dbc = @(x) (1+x).*(1-x);
N.rbc = 0; N.ubc = 0; N.lbc = 0;
u4 = N \ 0;

N = M;
N.lbc = @(x) (1+x).*(1-x);
N.rbc = @(x) (1+x).*(1-x);
N.dbc = @(x) (1+x).*(1-x);
N.ubc = @(x) (1+x).*(1-x);
u = N \ 0;

[xx, yy] = chebfun2.chebpts2(100);
A = feval(u,xx,yy);
B = feval(u1,xx,yy) + feval(u2,xx,yy) + feval(u3,xx,yy) + feval(u4,xx,yy);
pass(5) = ( norm( A - B ) < 1e10*tol);

% Check we can use the notation lap(u) = div(grad(u))
N = chebop2(@(u) -divergence(gradient(u)) );
N.lbc = 0; N.rbc = 0; N.dbc = 0; N.ubc = 0;
u = N \ 1;
N = chebop2(@(u) -laplacian(u) );
N.lbc = 0; N.rbc = 0; N.dbc = 0; N.ubc = 0;
exact = N \ 1;
pass(6) = ( norm( u - exact ) < tol);

% check lap(u) exists:
N = chebop2(@(u) -lap(u) );
N.lbc = 0; N.rbc = 0; N.dbc = 0; N.ubc = 0;
exact = N \ 1;
pass(7) = ( norm( u - exact ) < tol);
end


text:
function Ch4Ex1_main()
%% SOS-based Policy Iteration for a car suspension systems
% The function is tested in MATLAB R2014b
% Copyright 2015 Yu Jiang
% Contact Yu Jiang (yu.jiang@nyu.edu)

% System requirements:
% - MATLAB (Manually Tested in MATLAB R2014b)
% - MATLAB Symbolic Toolbox
% - SDPT3-4.0
% - SISOTOOLS (free to download at http://www.cds.caltech.edu/sostools/)

% You can download tools.zip and run setuptools.m in the folder.

syms x1 x2 x3 x4 real
mb = 300;    % kg
mw = 60;     % kg
bs = 1000;   % N/m/s
ks = 16000 ; % N/m
kt = 190000; % N/m
kn = 0.1*ks;

% State matrices
A = [ 0 1 0 0;
    [-ks -bs ks bs]/mb ; ...
    0 0 0 1;
    [ks bs -ks-kt -bs]/mw];
B = [ 0 0; 0 10000/mb ; 0 0 ; [kt -10000]/mw];
B = B(:,2);

% LQR feedback gains for the linearized system
% Klqr = lqr(A,B,eye(4),1);
% f(x)
f = A*[x1;x2;x3;x4] + [0;-kn*(x1-x3)^3/mb;0;kn*(x1-x3)^3/mw];

% Polynomial Weighting functions
q0 = 100*x1^2+x2^2+x3^2+x4^2;
rx = 1;%+(x1^2+x2^2+x3^2+x4^2);
vars = [x1;x2;x3;x4];

% Initialize the SOSp
prog = sosprogram(vars);

% The Lyapunov function V(x)
[prog,V] = sospolyvar(prog,monomials([x1;x2;x3;x4],2:4),'wscoeff');

% Objective of the SOSp
myObj = int(int(int(int(V,-.5,.5),-10,10),-.5,.5),-10,10);

% Add Inequality constraint to assure V is positive definite
prog = sosineq(prog,V-0.0001*(x1^2+x2^2+x3^2+x4^2));

% Add Inequality constraint to assure stability and performance
expr = -[diff(V,x1) diff(V,x2) diff(V,x3) diff(V,x4)]*rx*f-rx*q0;
prog = sosineq(prog,expr);

% Solve the SOSp
prog = sossolve(prog);

% Obtain the Initial Lyapunov function
V0 = sosgetsol(prog,V);
V_old = V0;

% Initializing the old contol policy
u_prev = zeros(size(x1));

% Iteration
for i=1:10
    clear prog V
    %------------------------------ SOSp Start ----------------------------
    prog = sosprogram(vars);
    [prog,V] = sospolyvar(prog,monomials([x1;x2;x3;x4],2:4),'wscoeff');
    prog = sosineq(prog,V_old - V);
    prog = sosineq(prog, V);
    u = -1/2*B'*[diff(V_old,x1) diff(V_old,x2) diff(V_old,x3) diff(V_old,x4)].';
    qfcn =rx*q0 + u'*u;
    expr = -[diff(V,x1) diff(V,x2) diff(V,x3) diff(V,x4)]*(rx*f+B*u)-qfcn;
    prog = sosineq(prog,expr);
    prog = sossetobj(prog, myObj);
    prog = sossolve(prog);
    %------------------------------ SOSp End ------------------------------
    V_ = sosgetsol(prog,V);
    V_old = V_;
end

% Save the improved value function
Vnew = V_;

%% Post-processing results
x0 = [0,0,0,0];  % Iniial Condition
tIntv = [0 3];
[t1,y1] = ode23s(@(t,x) LocalSuspSys(t,x,u), [0 3], x0);
[t,y] = ode23s(@(t,x) LocalSuspSys(t,x,0), tIntv, x0);

%% Plot Results
figure(1)
subplot(411);
plot(t1,y1(:,1),t,y(:,1), 'r:','linewidth',2);
xlabel('time (sec)', 'FontSize',12)
ylabel('x_1','FontSize',12)
hl = legend('Improved performance', 'Uncontrolled performance');
set(hl, 'FontSize', 12');
subplot(412);
plot(t1,y1(:,2),t,y(:,2),'r:','linewidth',2);
xlabel('time (sec)', 'FontSize',12)
ylabel('x_2','FontSize',12)
hl = legend('Improved performance', 'Uncontrolled performance');
set(hl, 'FontSize', 12');
subplot(413);
plot(t1,y1(:,3),t,y(:,3), 'r:','linewidth',2);
xlabel('time (sec)', 'FontSize',12)
ylabel('x_3','FontSize',12)
hl = legend('Improved performance', 'Uncontrolled performance');
set(hl, 'FontSize', 12');
subplot(414);
plot(t1,y1(:,4),t,y(:,4),'r:','linewidth',2);
xlabel('time (sec)', 'FontSize',12)
ylabel('x_4','FontSize',12)
hl = legend('Improved performance', 'Uncontrolled performance');
set(hl, 'FontSize', 12');

figure(2)
xx1 = -.4:.04:.4;
xx2 = -5:0.5:5;
vn = zeros(length(xx1),length(xx2));
v1 = zeros(length(xx1),length(xx2));
un=vn;
ulqr = un;
kn=vn;
k1=v1;
x3=0;
x4=0;
for i=1:length(xx1)
    x1 = xx1(i);
    for j=1:length(xx2)
        x2 = xx2(j);
        vn(i,j)=eval(Vnew);
        v1(i,j)=eval(V0);
    end
end
surf(xx1,xx2,vn')
hold on
surf(xx1,xx2,v1')
hold off
xlabel('x_1', 'FontSize', 12)
ylabel('x_2', 'FontSize', 12)
view(gca,[-30.5 28]);
% Create textarrow
annotation(gcf,'textarrow',[0.210714285714286 0.174535137214669],...
    [0.895238095238095 0.631440045897884],'TextEdgeColor','none','FontSize',12,...
    'String',{'V_0(x1,x2,0,0)'});

% Create textarrow
annotation(gcf,'textarrow',[0.139285714285714 0.186735060271868],...
    [0.183333333333333 0.386454388984516],'TextEdgeColor','none','FontSize',12,...
    'String',{'V_{10}(x1,x2,0,0)'});
% export_fig Ex3_cost -pdf -transparent

end

%% LocalSuspSys
% Dynamics of the nonlinear suspension system
function dx = LocalSuspSys(t,x,u)
mb = 300;    % kg
mw = 60;     % kg
bs = 1000;   % N/m/s
ks = 16000 ; % N/m
kt = 190000; % N/m

[x1,x2,x3,x4] = deal(x(1),x(2),x(3),x(4));

% State matrices
A = [ 0 1 0 0;
    [-ks -bs ks bs]/mb ; ...
    0 0 0 1;
    [ks bs -ks-kt -bs]/mw];
B = [ 0; 10000/mb ; 0 ; -10000/mw];
B1 = [ 0; 0 ; 0 ; kt/mw];

if ~isdouble(u)
    u = eval(u);
end

if t <= 0.001
    r = 10;
else
    r = 0;
end

dx = A*x + B*u + B1*r;
end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function varargout = fevalm(varargin)
%FEVALM   Evaluate a DISKFUN in polar coordinates.
%
%   Z = FEVALM(F, THETA, R) returns a matrix of values Z of size
%   length(R)-by-length(THETA). (R,THETA) are polar coordinates for the
%   evaluation points on the disk.  They should be vectors of doubles.
%   Calling this function as above is equivalent to making a meshgrid of
%   the vectors THETA and R and then using FEVAL to evaluate at that grid.
%
% See also DISKFUN/FEVAL.

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

[varargout{1:nargout}] = fevalm@separableApprox(varargin{:});

end


text:
function [infixOut, notaVAR] = pref2inf(prefixIn)
%PREF2INF   Convert an expression on prefix form to infix form
%   [INFIXOUT, NOTAVAR] = PREF2INF(PREFIXIN) goes recursively through the
%   expression PREFIXIN, which is on prefix form. The output, INFIXOUT, is a
%   string, representing the expression on infix form. NOTAVAR is a string, that
%   corresponds to the second variable in the kernel for FRED and VOLT
%   operators.

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

prefCounter = 1;
NOTAVAROUT = [];
infixOut = getInfix();
notaVAR = NOTAVAROUT;
NOTAVAROUT = [];

function infixOut = getInfix()
%GETINFIX   Return the infix form.

next = char(prefixIn(prefCounter,2));

% Steep through the prefix expression, dealing with each different type of
% tokens separately.
if ( ~isempty(strmatch('OP', next)) )
    % Increase the counter.
    prefCounter = prefCounter + 1;

    % We are dealing with a binary operator, so obtain two infix expressions.
    exp1 = getInfix();
    exp2 = getInfix();

    % We now return different outputs depending on which operator we have.
    switch next
        case 'OP='
            infixOut = [exp1, '=', exp2];

        case 'OP+'
            if ( strcmp(exp1, '0') && strcmp(exp2, '0') )
                % Don't want to introduce unnecessary zeros.
                infixOut = '';
            elseif  ( strcmp(exp1, '0') )
                infixOut = exp2;
            elseif ( strcmp(exp2, '0') )
                infixOut = exp1;
            else
                infixOut = ['(', exp1, '+', exp2, ')'];
            end

        case 'OP-'
            if ( strcmp(exp1, '0') && strcmp(exp2, '0') )
                infixOut = '';
            elseif  ( strcmp(exp1, '0') )
                infixOut = ['-',exp2];
            elseif ( strcmp(exp2, '0') )
                infixOut = exp1;
            else
                infixOut = ['(', exp1, '-', exp2, ')'];
            end

        case 'OP*'
            % Check whether we have some special arguments, which we know we can
            % simplify.
            if ( strcmp(exp1, '1') && strcmp(exp2, '1') )
                infixOut = '1';
            elseif  ( strcmp(exp1, '1') )
                infixOut = exp2;
            elseif ( strcmp(exp2, '1') )
                infixOut = exp1;
            elseif  ( strcmp(exp1, '-1') )
                infixOut = ['-', exp2];
            elseif ( strcmp(exp2, '-1') )
                infixOut = ['-', exp1];
            elseif ( strcmp(exp1, '0') || strcmp(exp2, '0') )
                infixOut = '0';
            elseif ( strcmp(exp1, '-0') || strcmp(exp2, '-0') )
                infixOut = '0';
            else
                infixOut = ['(', exp1, '.*', exp2, ')'];
            end

        case 'OP/'
            if ( strcmp(exp2, '1') )
                % Division by 1 is trivial.
                infixOut = exp1;
            else
                infixOut = ['(', exp1, './', exp2, ')'];
            end

        case 'OP^'
            infixOut = [exp1, '.^(', exp2, ')'];

        case {'OP>', 'OP>=', 'OP<', 'OP<='}
            nextSym = next(3:end);
            infixOut = ['(', exp1, nextSym, exp2, ')'];
    end

elseif ( strcmp(next, 'FUNC1') )
    % A function which has one argument.
    nextFun = char(prefixIn(prefCounter, 1));
    prefCounter = prefCounter + 1;
    funcArg = getInfix();
    infixOut = [nextFun, '(', funcArg , ')'];

elseif ( strcmp(next, 'FUNC2') )
    % A function with two arguments.
    nextFun = char(prefixIn(prefCounter, 1));
    prefCounter = prefCounter + 1;
    % Obtain two infix expressions.
    funcArg1 = getInfix();
    funcArg2 = getInfix();

    % Some methods need a special treatment.
    if ( (strcmp(nextFun, 'diff') || strcmp(nextFun, 'cumsum')) && ...
            strcmp(funcArg2, '1') )
        infixOut = [nextFun, '(', funcArg1, ')'];

    elseif ( any(strcmp(nextFun, {'fred', 'volt'})) )
        % For fred() and volt(), we need to find what arguments denote the
        % kernel variable.
        [ignored1, yFredVar, ignored2, ignored3, xFredVar] = ...
            stringParser.lexer(funcArg1, 'bvp');
        anonStr = ['@(' xFredVar{1} ',' yFredVar{1} ')'];
        infixOut = [nextFun, '(', anonStr, funcArg1 , ',', funcArg2 ,  ')'];
        NOTAVAROUT = [NOTAVAROUT ;  yFredVar];

    else
        infixOut = [nextFun, '(', funcArg1 , ',', funcArg2 ,  ')'];
    end

elseif ( strcmp(next, 'FUNC3') )
    % Methods with three arguments.
    nextFun = char(prefixIn(prefCounter,1));
    prefCounter = prefCounter + 1;
    funcArg1 = getInfix();
    funcArg2 = getInfix();
    funcArg3 = getInfix();
    infixOut = [nextFun, '(', funcArg1 , ',', funcArg2, ',' , funcArg3 ')'];

elseif ( strcmp(next(1:end-1), 'DER') )
    % Differentiation via the ' symbol.
    prefCounter = prefCounter + 1;
    derivArg = getInfix();
    derivOrder = next(4:end);
    if ( ~strcmp(derivOrder, '1') )
        infixOut = ['diff(', derivArg, ',' derivOrder, ')'];
    else
        infixOut = ['diff(', derivArg, ')'];
    end

elseif ( ~isempty(strmatch('UN', next)) )
    % Unary operators
    nextUnary = char(prefixIn(prefCounter, 1));
    prefCounter = prefCounter + 1;
    unaryArg = getInfix();
    % Only care about unary -
    if ( strcmp(nextUnary, '-') )
        infixOut = [nextUnary, unaryArg];
    else
        infixOut = unaryArg;
    end

elseif ( strmatch(next, 'COMMA') )
    % The , symbol.
    prefCounter = prefCounter + 1;
    exp1 = getInfix();
    exp2 = getInfix();
    infixOut = ['(', exp1, ';', exp2, ')'];
else
    % We have arrived at a terminal token. Return it and increase the counter.
    infixOut = char(prefixIn(prefCounter,1));
    prefCounter = prefCounter + 1;
end

end

% Initialization for simplification.
newlen = length(infixOut);
len = inf;

% Do some simplifications on the returned string. In particular, we want to
% change -- to +, a -+ to -, a +- to -, and ++ to +.
while ( newlen ~= len )     % Loop until length doesn't change anymore.
    len = newlen;

    k = 1;
    while ( k < (numel(infixOut) - 1) )
        if ( strcmp(infixOut(k), '-') )
            if ( strcmp(infixOut(k+1), '-') )
                infixOut(k) = '+';
                infixOut(k+1) = [];
            elseif ( strcmp(infixOut(k+1), '+') )
                infixOut(k+1) = [];
            else
                k = k+1;
            end
        elseif ( strcmp(infixOut(k), '+') )
            if ( strcmp(infixOut(k+1), '-') )
                infixOut(k) = '-';
                infixOut(k+1) = [];
            elseif ( strcmp(infixOut(k+1), '+') )
                infixOut(k+1) = [];
            else
                k = k+1;
            end
        else
            k = k+1;
        end
    end

    newlen = length(infixOut);

end

% Clear global variables
prefixIn = [];
prefCounter = [];

end


text:
function [vx, vy] = vbeta2vxvy(v, beta)

vx = v*cos(beta);
vy = v*sin(beta);

author: TUMFTM


text:
%%*****************************************************************
%% schurmat_lblk: compute A*D*A'
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function [schur,UU,EE] = schurmat_lblk(blk,At,par,schur,UU,EE,p,dd);

   global idxdenAl

   iter = par.iter;
   n = sum(blk{p,2});

   if (iter==1)
      idxdenAl{p} = checkdense(At{p}');
   end
   ddsch = dd{p};
   if ~isempty(idxdenAl{p});
      idxden = idxdenAl{p};
      len = length(idxden);
      Ad = At{p}(idxden,:)' *spdiags(sqrt(ddsch(idxden)),0,len,len);
      UU = [UU, Ad];
      if isempty(EE)
         count = 0;
      else
         count = max(max(EE(:,1)),max(EE(:,2)));
      end
      tmp = count + [1:len]';
      EE = [EE; [tmp, tmp, -ones(len,1)] ];
      ddsch(idxden) = zeros(len,1);
   end
   schurtmp = At{p}' *spdiags(ddsch,0,n,n) *At{p};
   schur = schur + schurtmp;
%%*******************************************************************


text:
%
%
%this script produces the figures of the paper scattering_fractal_analysis
%september 2013
%Joan Bruna
%
%

close all;
%clear all;
%startup;

clear filt_opt;
clear scat_opt;
%set the maximum analysis scale
filt_opt.J = 14;
filt_opt.filter_format='fourier';
%size of realizations
N=2^(4+filt_opt.J);

%choose the wavelet
filt_opt.filter_type='selesnick_1d';
%filt_opt.filter_type='morlet_1d';
scat_opt.M=2;

scat_opt.path_margin= Inf;
[Wop, filters]=wavelet_factory_1d(N,filt_opt,scat_opt);

R=4;
options.Wop=Wop;
options.J=filt_opt.J;

poisson=@(alpha) (cumsum(double(rand(N,1) < N^(-alpha))));
dpoisson=@(alpha) ((double(rand(N,1) < N^(-alpha))));

[S{1},T{1},Tu{1},ex{1}] = scat_renorm_1d_scatnet(dpoisson, 0.6, options, R);

figure
plot_full_transfer(log2(T{1}));

oldopts.J=14;
oldopts.filters=selesnick_bis([N 1], oldopts);
oldopts.fullscatt=1;
oldopts.oversampling=oldopts.J;
[S{2},T{2},Tu{2},ex{2}] = scatt_renorm_estimation(dpoisson, 0.6, oldopts, R);

figure
plot_full_transfer(log2(T{2}));

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
%% Introduction to Genetic Algorithms
% Oren Rosen
% The MathWorks
% 11/7/2007
%
% This script is for a brief introduction to the MATLAB environment
% and to introduce the concept of Genetic Algorithms using Rastrigin's
% function.

%% Intro to MATLAB 1
x = 2;
y = 3;
z = ras([x,y])

%% Intro to MATLAB 2
x = -5:5

%% Define x and y vectors that sample the interval [-5,5]
x = -5:0.1:5;
y = -5:0.1:5;

%% Evaluate Rastrigin's function
% allcomb creates all combinations of elements in x,y
allcomb(x,y)

z = ras(allcomb(x,y));

%% Reshape resulting vector into a grid format
z = reshape(z,length(x),length(y));

%% Visualize the surface
surf(x,y,z);
contourf(x,y,z);


text:
function [x, i] = pcg_scs(A,b,x,M,rho_x,max_iters,tol,verbose)
% M is inverse preconditioner
r=b-(rho_x*x+A'*(A*x));
z = M.*r;
p = z;
ip = r'*z;
for i=1:max_iters
    Ap=rho_x*p + A'*(A*p);
    alpha= ip/(p'*Ap);
    x=x+alpha*p;
    r=r-alpha*Ap;
    if norm(r)<tol
        if verbose
            fprintf('CG took %i iterations to converge, resisdual %4f <= tolerance %4f\n',i,resid,tol)
        end
        return;
    end
    z = M.*r;
    ipold = ip;
    ip = z'*r;
    beta =  ip / ipold;
    p=z+beta*p;
end
if verbose
    fprintf('CG did not converge within %i iterations, resisdual %4f > tolerance %4f\n',max_iters,resid,tol)
end
end

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function Population = subNSGAII(Problem,Population,Operator,N)
% Sub-optimizer in LSMOF (NSGA-II)

%------------------------------- 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".
%--------------------------------------------------------------------------

% This function is written by Cheng He

    FrontNo  = NDSort(Population.objs,inf);
    CrowdDis = CrowdingDistance(Population.objs,FrontNo);
    if Operator == 1
        MatingPool = TournamentSelection(2,N,FrontNo,-CrowdDis);
        Offspring  = OperatorGA(Problem,Population(MatingPool));
    else
        MatingPool1 = TournamentSelection(2,N,FrontNo,-CrowdDis);
        MatingPool2 = TournamentSelection(2,N,FrontNo,-CrowdDis);
        MatingPool3 = TournamentSelection(2,N,FrontNo,-CrowdDis);
        Offspring   = OperatorDE(Problem,Population(MatingPool1),Population(MatingPool2),Population(MatingPool3));
    end
    Population = EnvironmentalSelection([Population,Offspring],N);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO