text:
classdef MrDimInfo < MrCopyData
    % Holds all dimensionality info (r.g. range/dimLabels/units) of multidimensional data
    %
    %   See also demo_dim_info MrImage.select

    % Author:   Saskia Bollmann & Lars Kasper
    % Created:  2016-01-23
    % Copyright (C) 2016 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/>.

    % the following properties can be fully derived from sampling points, ...
    % but are stored for convenience
    properties (Dependent)
        nDims;      % number of dimensions in dataset, default: 6

        nSamples;   % vector [1,nDims] of number of elements per dimension

        % vector [1,nDims] of resolutions for each dimension, ...
        % i.e. distance (in specified units) of adjacent elements,
        % NaN for non-equidistant spacing of elements
        % Example: TE = 2, 20, 35 ms => resolution = NaN
        resolutions;

        % vector [2,nDims] of [firstSample; lastSample] for each dimension ...
        ranges;

    end % properties (dependent)

    properties
        % cell(1,nDims) of string dimLabels for each dimension
        % default: {'x', 'y', 'z', 't', 'coil, 'echo'}
        dimLabels = {};

        % cell(1,nDims) of strings describing unit; '' for unit-free dims
        % default: {'mm', 'mm', 'mm', 's', '', 'ms'};
        units = {};

        % cell(1,nDims) of sampling position vectors for each dimension
        % Note: The entries for each dimension have to be numerical
        samplingPoints = {};

        % [1,nDims] vector of sampling width vectors for each dimension ...
        % Note: Typically, this will correspond to diff(samplingPoints);
        % However, if sampling does not cover the full interval between
        % consecutive points, it should be noted here
        % Example:  If z-coordinate has slice thickness and slice gap
        %           samplingWidths = sliceThickness
        %           resolutions = sliceThickness + sliceGap
        samplingWidths = [];

    end % properties

    methods

        function this = MrDimInfo(varargin)
            % Constructor of class, call via MrDimInfo('propertyName', propertyValue, ...
            % ...) syntax
            %               OR
            % MrDimInfo(fileName)
            %
            %               OR
            % MrDimInfo(fileArray) - additional dims are added via filename
            %
            %               OR
            % MrDimInfo(dimInfoStruct) - MrDimInfo object is costructed
            % from struct with same fields (especially for load and save)
            %
            % See also MrDimInfo.set_dims
            %
            % IN
            %   varargin    PropertyName/Value-pairs of MrDim-Info properties to be
            %               changed, e.g. resolutions, nSamples, units etc.
            %
            %   Properties:
            %
            %   'units'           cell(1,nDims) of strings describing unit; '' for unit-free dims
            %	'dimLabels'       cell(1,nDims) of string dimLabels for each changed dimension
            %
            %
            %   (1): 1st variant: explicit setting of sampling points for dimension(s)
            %
            %   'samplingPoints'  cell(1,nDims) of index vectors for each dimension
            %
            %
            %   (2)-(6): Other variants depend on setting some of the following parameters
            %
            %   'nSamples'              [1, nDims] number of samples per dimension
            %   'ranges'                [2, nDims] first and last sample per dimension
            %   'resolutions'           [1, nDims] resolution (in units), i.e. spacing
            %                           between sampling points per dimension
            %
            %   'arrayIndex'            index of samplingPoint to be set
            %   'samplingPoint'         value of sampling point at position arrayIndex
            %
            %   'firstSamplingPoint'    special case of samplingPoint, arrayIndex = 1 set
            %   'lastSamplingPoint'     special case of samplingPoint, arrayIndex = end set
            %
            %   Variants:
            %       (2) nSamples + ranges: sampling points computed automatically via
            %               samplingPoint(k) = ranges(1) + (ranges(2)-ranges(1))/nSamples*(k-1)
            %       (3) nSamples + resolutions + samplingPoint + arrayIndex:
            %               from samplingPoints(arrayIndex) = samplingPoints, others
            %               are constructed via
            %               [...    samplingPoint-resolution
            %                       samplingPoint
            %                       samplingPoint+resolution ...]
            %               until nSamples are created in total.
            %       (4) nSamples + resolutions + firstSamplingPoint:
            %               as (3), assuming arrayIndex = 1
            %       (5) nSamples + resolutions + lastSamplingPoint:
            %               as (3), assuming arrayIndex = end
            %       (6) nSamples Or resolution Or (samplingPoint+arrayIndex)
            %               missing input value taken from defaults
            %               nSamples        -> resolutions = 1, 1st sample 1
            %               resolutions     -> nSamples = 2, 1st sample 1
            %
            %           and always samplingWidths = resolution

            if nargin == 1 % single file, file array or struct is given
                if isstruct(varargin{1}) % is struct
                    this.update_properties_from(varargin{1});
                else % is file/file array
                    fileInput = varargin{1}; % extract input
                    isSingleFile = 0;
                    isFile = 1;
                    % check whether single filename or folder is given
                    if ischar(fileInput)
                        % determine whether file or folder
                        [~, ~, ext] = fileparts(fileInput);
                        if ~isempty(ext) % single file
                            fileName = fileInput;
                            isSingleFile = 1;
                        else
                            % or folder
                            fileInput = cellstr(spm_select('FPList',...
                                fileInput, '^*.*'));
                        end
                    elseif iscell(fileInput) && numel(fileInput) == 1
                        % fileArray with only one entry
                        fileName = fileInput{1}; % extract from cell array
                        isSingleFile = 1;
                    elseif isa(varargin{1}, 'MrAffineTransformation')
                        this.set_from_affine_geometry(varargin{1});
                        isFile = 0;
                    end

                    if isFile
                        if isSingleFile
                            % load from single file
                            this.load(fileName);
                        else
                            % fileArray is given
                            % load dimInfo from first file in file array
                            this.load(fileInput{1});
                            this.set_from_filenames(fileInput);
                        end
                    end
                end
            else % prop/value pairs given

                propertyNames = varargin(1:2:end);
                propertyValues = varargin(2:2:end);
                % Find nSamples property, and corresponding value to determine
                % dimension
                iArgNsamples = tapas_uniqc_find_string(propertyNames, 'nSamples');
                iArgSamplingPoints = tapas_uniqc_find_string(propertyNames, 'samplingPoints');
                iArgRanges = tapas_uniqc_find_string(propertyNames, 'ranges');

                hasNsamples = ~isempty(iArgNsamples);
                hasExplicitSamplingPoints = ~isempty(iArgSamplingPoints);
                hasRanges = ~isempty(iArgRanges);

                if hasExplicitSamplingPoints
                    nDims = numel(propertyValues{iArgSamplingPoints});
                elseif hasNsamples
                    % otherwise, allow empty constructor for copyobj-functionality
                    nDims = numel(propertyValues{iArgNsamples});
                elseif hasRanges
                    nDims = numel(propertyValues{iArgRanges})/2;
                else
                    % guessed number of update dimensions
                    % find shortest given input to dimInfo and determine
                    % dimensionality from that
                    nDims =[];
                    for p = 1:numel(propertyNames)
                        nDims(p) = numel(propertyValues{p});
                    end
                    nDims = min(nDims);

                    % populate nDims with non-samplingpoint parameters by adding dimensions
                    % , s.t. samplingPoints remains an empty cell for each
                    % dimension
                    % only samplingWidhts, units and dimLabels could be
                    % in varargin, so OK to add
                    this.add_dims(1:nDims, varargin{:});

                end

                % allows empty constructor for copyobj
                if ~isempty(nDims)
                    this.set_dims(1:nDims, varargin{:});
                end
            end

        end
        % NOTE: Most of the methods are saved in separate function.m-files in this folder;
        %       except: constructor, delete, set/get methods for properties.
        %

        function nDims = get.nDims(this)
            nDims = numel(this.nSamples);
        end

        function this = set.nDims(this, nDimsNew)
            % Sets number of dimensions of dimInfo, via add_dims or ...
            % remove_dims
            nDimsOld = this.nDims;

            if nDimsNew > nDimsOld
                this.add_dims((nDimsOld+1):nDimsNew);
            elseif nDimsNew < nDimsOld
                this.remove_dims((nDimsNew+1):nDimsOld);
            end
        end

        function this = set.nSamples(this, nSamplesNew)
            % Changes nSamples by keeping given resolution and adding samples ...
            % at end of samplingPoints-vectors
            % TODO: Can we do this via set_dims? Or is this a problem for
            % empty dimInfos which are just created?
            % No: more for when we don't have a defined resolution
            % (non-equidistant!)

            nSamplesOld = this.nSamples;
            if numel(nSamplesNew) ~= numel(nSamplesOld)
                error('tapas:uniqc:MrDimInfo:NDimsChanged', ...
                    'nDims cannot change via nSamples, use add_dims instead');
            end

            iChangedDims = find(nSamplesOld ~= nSamplesNew);

            for iDim = iChangedDims
                nOld = nSamplesOld(iDim);
                nNew = nSamplesNew(iDim);
                if nOld > nNew
                    this.samplingPoints{iDim}((nNew+1):end) = [];
                else
                    if nOld == 0 % no samples before, create all new
                        offset = 0;
                    else % start samples from last existing
                        offset = this.samplingPoints{iDim}(nOld);
                    end

                    % set default resolution to 1 for adding samples
                    res = this.resolutions(iDim);
                    if isnan(res)
                        res = 1;
                    end

                    this.samplingPoints{iDim}((nOld+1):nNew) = ...
                        offset + (1:(nNew-nOld))*res;
                end
            end

        end

        function nSamples = get.nSamples(this)
            if isempty(this.samplingPoints)
                nSamples = [];
            else
                nSamples = cell2mat(cellfun(@numel, this.samplingPoints, ...
                    'UniformOutput', false));
            end
        end

        function this = set.resolutions(this, resolutionsNew)
            % Changes resolutions by setting the new resolution via
            % changeDim, the first sample is preserved
            resolutionsOld = this.resolutions;

            if numel(resolutionsNew) ~= numel(resolutionsOld)
                error('tapas:uniqc:MrDimInfo:NDimsChanged', ...
                    'nDims cannot change via resolutions, use add_dims instead');
            end

            iChangedDims = find(~arrayfun(@isequaln, resolutionsOld, resolutionsNew));
            this.set_dims(iChangedDims, 'resolutions', resolutionsNew(iChangedDims));
        end

        function resolutions = get.resolutions(this)
            if isempty(this.samplingPoints)
                resolutions = [];
            else

                for iDim = 1:this.nDims
                    res = unique(diff(this.samplingPoints{iDim}));
                    switch numel(res)
                        case 0 % one element samplingPoints, take its value (?)
                            if ~isempty(this.samplingWidths) && ...
                                    numel(this.samplingWidths) >= iDim && ...
                                    ~isnan(this.samplingWidths(iDim))
                                resolutions(iDim) = this.samplingWidths(iDim);
                            else
                                resolutions(iDim) = NaN;
                            end
                            % resolutions(iDim) = this.samplingPoints{iDim};
                        case  1 % single resolution, return it
                            resolutions(iDim) = res;
                        otherwise % if no unique resolution,
                            % first check if within single floating precision,
                            % accept that as same!
                            if max(abs(diff(res))) < eps(single(1))
                                resolutions(iDim) = mean(res);
                            else
                                %  otherwise, really non-equidistant spacing, return NaN for this
                                %  dim
                                resolutions(iDim) = NaN;
                            end
                    end
                end
            end
        end

        function this = set.ranges(this, rangesNew)
            % Changes ranges by keeping given nSamples and adjusting ...
            % samplingPoints (i.e. spacing i.e. resolution)
            rangesOld = this.ranges;

            if numel(rangesNew) ~= numel(rangesOld)
                error('tapas:uniqc:MrDimInfo:NDimsChanged', ...
                    'nDims cannot change via ranges, use add_dims instead');
            end

            iChangedDims = union(find(rangesOld(1,:) ~= rangesNew(1,:)), ...
                find(rangesOld(2,:) ~= rangesNew(2,:)));

            this.set_dims(iChangedDims, 'ranges', rangesNew(:,iChangedDims), ...
                'nSamples', this.nSamples(iChangedDims));
        end

        function ranges = get.ranges(this)
            ranges = [first(this); last(this)];
        end

        function centerSamples = center(this, iDim)
            % return center sample for given dimension (i.e. index ceil(nSamples/2))
            if nargin < 2
                iDim = 1:this.nDims;
            end

            centerSamples = nan(1, numel(iDim));
            for d = iDim
                if ~isempty(this.samplingPoints{d})
                    centerSamples(find(d==iDim)) = this.samplingPoints{d}(...
                        ceil(this.nSamples(d)/2));
                end
            end

        end

        function firstSamples = first(this, iDim)
            if nargin < 2
                iDim = 1:this.nDims;
            end

            firstSamples = nan(1, numel(iDim));
            for d = iDim
                if ~isempty(this.samplingPoints{d})
                    firstSamples(find(d==iDim)) = this.samplingPoints{d}(1);
                end
            end

        end

        % return last sampling point for all or given dimensions
        function lastSamples = last(this, iDim)
            if nargin < 2
                iDim = 1:this.nDims;
            end

            lastSamples = nan(1,numel(iDim));
            for d = iDim
                if ~isempty(this.samplingPoints{d})
                    lastSamples(find(d==iDim)) = this.samplingPoints{d}(end);
                end
            end
        end

        function [iDim, isValidLabel] = get_dim_index(this, dimLabel, varargin)
            % return index of dimension(s) given by a dimLabel
            % IN
            %   dimLabel  dimension label string (or array of strings).
            %             or dimension number or cell of dim numbers (for
            %             compatibility)
            %   varargin
            %   'invert'    true or false (default)
            %               if true, all other indices not within dimLabel
            %               are returned
            %
            % OUT
            %   iDim            index of dimension with corresponding label
            %   isValidLabel    [nLabels,1] returns for each given label 1/0
            %                   i.e. whether it is indeed a label of dimInfo
            defaults.invert = false;
            args = tapas_uniqc_propval(varargin,defaults);
            if isnumeric(dimLabel) % (vector of) numbers
                iDim = dimLabel;
                % cell of numbers:
            elseif iscell(dimLabel) && isnumeric(dimLabel{1})
                iDim = cell2mat(dimLabel);
            else % string or cell of strings
                isExact = 1;
                iDim = tapas_uniqc_find_string(this.dimLabels, dimLabel, isExact);
                if iscell(iDim)
                    isValidLabel = ~cellfun(@isempty, iDim);
                    iDim = iDim(isValidLabel); % remove unfound dimensions;
                    iDim = cell2mat(iDim)';
                else
                    isValidLabel = ~isempty(iDim);
                end
            end
            if args.invert
                iDim = setdiff(1:this.nDims,iDim);
            end
        end

    end

end

author: translationalneuromodeling

repo: tapas


text:
function bestEpoch = findBestEpoch(expDir, varargin)
%FINDBESTEPOCH finds the best epoch of training
%   FINDBESTEPOCH(EXPDIR) evaluates the checkpoints
%   (the `net-epoch-%d.mat` files created during
%   training) in EXPDIR
%
%   FINDBESTEPOCH(..., 'option', value, ...) accepts the following
%   options:
%
%   `priorityMetric`:: 'classError'
%    Determines the highest priority metric by which to rank the
%    checkpoints.
%
%   `prune`:: false
%    Removes all saved checkpoints to save space except:
%
%       1. The checkpoint with the lowest validation error metric
%       2. The last checkpoint
%
% Copyright (C) 2017 Samuel Albanie
% Licensed under The MIT License [see LICENSE.md for details]

  opts.prune = false ;
  opts.priorityMetric = 'classError' ;
  opts = vl_argparse(opts, varargin) ;

  lastEpoch = findLastCheckpoint(expDir);
  if ~lastEpoch, return ; end % return if no checkpoints were found

  bestEpoch = findBestValCheckpoint(expDir, opts.priorityMetric);
  if opts.prune
    preciousEpochs = [bestEpoch lastEpoch];
    removeOtherCheckpoints(expDir, preciousEpochs);
    fprintf('----------------------- \n');
    fprintf('directory cleaned: %s\n', expDir);
    fprintf('----------------------- \n');
  end

% -------------------------------------------------------------------------
function removeOtherCheckpoints(expDir, preciousEpochs)
% -------------------------------------------------------------------------
  list = dir(fullfile(expDir, 'net-epoch-*.mat')) ;
  tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
  epochs = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
  targets = ~ismember(epochs, preciousEpochs);
  files = cellfun(@(x) fullfile(expDir, sprintf('net-epoch-%d.mat', x)), ...
          num2cell(epochs(targets)), 'UniformOutput', false);
  cellfun(@(x) delete(x), files)

% -------------------------------------------------------------------------
function bestEpoch = findBestValCheckpoint(expDir, priorityMetric)
% -------------------------------------------------------------------------
  lastEpoch = findLastCheckpoint(expDir) ;
  if strcmp(priorityMetric, 'last'), bestEpoch = lastEpoch ; return ; end
  % handle the different storage structures/error metrics
  path = fullfile(expDir, sprintf('net-epoch-%d.mat', lastEpoch)) ;
  try
    data = load(path) ;
  catch
		msg = 'checkopint at %s was malformed, trying agin in 10 secs....\n' ;
		warning(msg, path) ; pause(10) ; data = load(path) ;
  end
  if isfield(data, 'stats')
    valStats = data.stats.val;
  elseif isfield(data, 'info')
    valStats = data.info.val;
  elseif isfield(data, 'state')
    valStats = data.state.stats.val ;
  else
    error('storage structure not recognised');
  end

  ascending = {'mAP', 'accuracy'} ;
  descending = {'top1error', 'error', 'mbox_loss', 'class_loss'} ;

  % find best checkpoint according to the following priority
  metrics = [{priorityMetric} ascending descending] ;

  for i = 1:numel(metrics)
    if isfield(valStats, metrics{i})
      errorMetric = [valStats.(metrics{i})] ;
      selectedMetric = metrics{i} ;
      break ;
    end
  end

  assert(logical(exist('errorMetric', 'var')), 'error metrics not recognized') ;
  if ismember(selectedMetric, ascending)
    pick = @max ;
  else
    pick = @min ;
  end

  [~, bestEpoch] = pick(errorMetric);

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(expDir)
% -------------------------------------------------------------------------
  list = dir(fullfile(expDir, 'net-epoch-*.mat')) ;
  tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
  epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
  epoch = max([epoch 0]) ;

repo: MFT


text:
function [mm2] = cm22mm2(cm2)
% Convert area from square centimeters to square millimeters.
% Chad A. Greene 2012
mm2 = cm2*100;

author: Sable

repo: mcbench-benchmarks


text:
%% Copyright (C) 2015, 2016 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym ei (@var{x})
%% Symbolic exponential integral (Ei) function.
%%
%% Definition and example:
%% @example
%% @group
%% syms x
%% f = ei(x)
%%   @result{} f = (sym) Ei(x)
%% rewrite(f, 'Integral')         % doctest: +SKIP
%%   @result{} (sym)
%%       x
%%       ⌠
%%       ⎮   t
%%       ⎮  ℯ
%%       ⎮  ── dt
%%       ⎮  t
%%       ⌡
%%       -∞
%% @end group
%% @end example
%% (@strong{Note} rewriting as an integral is not yet supported.)
%%
%% Other examples:
%% @example
%% @group
%% diff(f)
%%   @result{} (sym)
%%        x
%%       ℯ
%%       ──
%%       x
%% @end group
%% @end example
%%
%% @seealso{@@sym/expint}
%% @end defmethod

function y = ei(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('Ei', x);
end

%!test
%! syms x
%! f = ei(sym(0));
%! assert (double(f) == -inf)

%!test
%! D = [1.895117816355937  4.954234356001890];
%! A = ei(sym([1 2]));
%! assert (all (abs(double(A) - D) < 1e-15))

repo: octsympy


text:
function model = mogOptimise(model, display, iters)

% MOGOPTIMISE Optimise an MOG model.
% FORMAT
% DESC optimises an mixtures of Gaussians model via the
% expectation maximisation algorithm.
% ARG model : the model to be optimised.
% RETURN model : the optimised model.
%
% SEEALSO : mmpcaCreate, modelOptimise
%
% COPYRIGHT : Neil D. Lawrence, 2006

% MLTOOLS

diffll = 1;
iter = 0;
ll = mogLowerBound(model);
while abs(diffll)>1e-6 & iter<iters
  iter = iter + 1;
  model =  mogEstep(model);
  if display > 1
    [ll, oldll] = boundCheck(model, ll, 'E-step');
  end
  model = mogUpdatePrior(model);
  if display > 1
    [ll, oldll] = boundCheck(model, ll, 'Prior Update');
  end
  model = mogUpdateMean(model);
  if display > 1
    [ll, oldll] = boundCheck(model, ll, 'Mean Update');
  end
  model = mogUpdateCovariance(model);
  if display > 1
    [ll, oldll] = boundCheck(model, ll, 'Covariance Update');
  end
  if display > 1
  else
    oldll = ll;
    ll = mogLowerBound(model);
  end
  diffll = ll -oldll;
  if display
    fprintf('Iteration %d log-likelihood: %2.6f\n', iter, ll)
  end
end

function [ll, oldll] = boundCheck(model, oldll, step)

% BOUNDCHECK Helper function for checking bound.

ll = mogLowerBound(model);
diffll = ll - oldll;
if ll -oldll < 0
  warning(['Log likelihood went down by ' num2str(diffll) 'in ' ...
           'step: ' step ' in mogOptimise'])
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function oldPath = ezPath(varargin)
%EZPATH add path more easier.
%USAGE
% ezPath folder1 folder2 ...

%SEE ALSO addpath rmpath savepath
%
%COPYRIGHT Zhenqiang Ying <yingzhenqiang-at-gmail-dot-com>

oldPath = path;

cfg.SavePath = true;
cfg.FullPath = true;

if isstruct(varargin{1})
	cfg = loaddefault(varargin{1}, cfg, true); % time-consuming
    varargin = varargin(2:end);
end

inputPath = strjoin(varargin, ';');
pathList = strsplit(inputPath, ';');

try
    caller = callerFile;
    defaultPath = caller.path;
catch
    defaultPath = cd;
end

if cfg.FullPath
    for pathStr = pathList
        [p, func] = strsplit(pathStr{1}, '>');p = p{1};
        if ~ismember(':', p) % is not fullpath
            p = [defaultPath filesep p];
        end
        if isempty(func) || ~exist(func{1}, 'file')
            if p(end) == '*'
                genpath(p(1:end-1))
                addpath(genpath(p(1:end-1)));
            else
                addpath(p);
            end
            %this.path{end+1} = p;
        end
    end
else
    addpath(inputPath);
end

if cfg.SavePath
	savepath
end

author: baidut

repo: BIMEF

save_path: github-repos/MATLAB/baidut-BIMEF


text:
function REP = REPSelection(PopObj,N,div)
% Select one of the particles in REP as the global best position for each
% particle

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

    NoP = size(PopObj,1);

    %% Calculate the grid location of each solution
    fmax = max(PopObj,[],1);
    fmin = min(PopObj,[],1);
    d    = (fmax-fmin)/div;
    fmin = repmat(fmin,NoP,1);
    d    = repmat(d,NoP,1);
    GLoc = floor((PopObj-fmin)./d);
    GLoc(GLoc>=div) = div - 1;
    GLoc(isnan(GLoc)) = 0;

    %% Detect the grid of each solution belongs to
    [~,~,Site] = unique(GLoc,'rows');

    %% Calculate the crowd degree of each grid
    CrowdG = hist(Site,1:max(Site));

    %% Roulette-wheel selection
    TheGrid = RouletteWheelSelection(N,CrowdG);
    REP     = zeros(1,N);
    for i = 1 : length(REP)
        InGrid = find(Site==TheGrid(i));
        Temp   = randi(length(InGrid));
        REP(i) = InGrid(Temp);
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% ctftopo
% EEGLAB script for plotting topo maps of CTF data
% Thomas Ferree @ UCSF
% Revised 2/17/2004

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

clear;

% toggle saving eeglab .set, 0 = no, 1 = yes
saveSet = 0;

% set machine: 1 = TiBoard, 2 = Glutamate/Acetylcholine/Brainmap
machine = 2;

% set input and output paths
if machine == 1
	pathdata = '/Users/tferree/Documents/Data/Attention/HighN/';
	pathersp = '/Users/tferree/Documents/Research/Attention/HighN/Numerical/';
	pathtemp = '/Users/tferree/Documents/Research/Attention/HighN/Temp/';
else
	addpath /home/tferree/EEGLAB42/;
	addpath /home/tferree/Attention/HighN/ctf2eeglab/;
	pathdata = '/home/cdale/HighN/';
	pathersp = '/home/tferree/Attention/HighN/Numerical/';
	pathtemp = '/home/tferree/Attention/HighN/Temp/';
end

% set subject and conditions

%subject = 'HM'
%date = '20040126';

%subject = 'GN'
%date = '20040109';

subject = 'CJ'
date = '20040116';

cue = 'CueL'

% load ERSP results
load([pathersp 'ERSP_' subject '_' cue '.mat']);

% call Darren's scripts
%ctffolder = [pathdata subject '_DNL_' date '_HiN_' cue 'InvNT' '-fEOG.ds'];
ctffolder = [pathdata subject '_epoch' filesep 'processed' filesep subject '_DNL_' date '_HiN_' cue 'InvNT' '-fEOG.ds'];

ctfchannels = 'meg';
ctftime = 'all';
ctftrials = [1];
ctf = ctf_read(ctffolder,ctfchannels,ctftime,ctftrials);
ctf2eeglab;

% make topo maps
freqs = allfreqs(1,:,1);
times = alltimes(1,:,1);
maxcolor = 5;

for plottime = 0:100:1800

    if plottime <= 9
        timestring = ['000' num2str(plottime)];
    elseif plottime >= 10 & plottime <= 99
        timestring = ['00' num2str(plottime)];
    elseif plottime >= 100 & plottime < 1000
        timestring = ['0' num2str(plottime)];
    else
        timestring = num2str(plottime);
    end

    timefreqlist = [plottime 5; plottime 7; plottime 9; plottime 11; plottime 13; plottime 15; plottime 17];
    h3 = figure;
    tftopo(allersp,times,freqs,'mode','ave','signifs', allerspboot,'limits',[min(times) max(times) 1 20 -maxcolor, maxcolor],'vert',[0 1000],'chanlocs',EEG.chanlocs,'timefreqs',timefreqlist,'electrodes','off','maplimits',[-maxcolor,maxcolor]);
    saveas(h3,[pathtemp 'Topo_' subject '_' cue '_' timestring 'msec.jpg'],'jpg');
    close(h3);

end

author: PatternRecognition

repo: OpenBMI


text:
function beta_nc_test01 ( )

%*****************************************************************************80
%
%% BETA_NC_TEST01 tests BETA_NONCENTRAL_CDF against tabulated values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2008
%
%  Author:
%
%    John Burkardt
%
  error_max = 1.0E-10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BETA_NC_TEST01:\n' );
  fprintf ( 1, '  Compare tabulated values of the noncentral\n' );
  fprintf ( 1, '  incomplete Beta Function against values\n' );
  fprintf ( 1, '  computed by BETA_NONCENTRAL_CDF.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1,'        A        B   LAMBDA      X      ' );
  fprintf ( 1,'  CDF                       CDF                     DIFF\n' );
  fprintf ( 1,'                                        ' );
  fprintf ( 1,' (tabulated)               (BETA_NC)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, a, b, lambda, x, fx ] = beta_noncentral_cdf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = beta_noncentral_cdf ( a, b, lambda, x, error_max );

    fprintf ( 1, '  %7.1f  %7.1f  %7.1f  %7.3f  %24.16e  %24.16e  %10.4e\n', ...
    a, b, lambda, x, fx, fx2, abs ( fx - fx2 ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% BOXPLOT Displays box plots of multiple data samples.
%    BOXPLOT(X) produces a box plot of the data in X.  If X is a matrix there
%    is one box per column, and if X is a vector there is just one box. On
%    each box, the central mark is the median, the edges of the box are the
%    25th and 75th percentiles, the whiskers extend to the most extreme
%    datapoints the algorithm considers to be not outliers, and the outliers
%    are plotted individually.
%
%    BOXPLOT(X,G) specifies one or more grouping variables G, producing a
%    separate box for each set of X values sharing the same G value or
%    values.  Grouping variables must have one row per element of X, or one
%    row per column of X. Specify a single grouping variable in G by using a
%    vector, a character array, a cell array of character vectors, a string
%    array, or a vector categorical array; specify multiple grouping
%    variables in G by using a cell array of these variable types, such as
%    {G1 G2 G3}, or by using a matrix.  If multiple grouping variables are
%    used, they must all be the same length.  Groups that contain a NaN or
%    an empty string ('') in a grouping variable are omitted, and are not
%    counted in the number of groups considered by other parameters.
%
%    By default, character and string grouping variables are sorted in the
%    order they initially appear in the data, categorical grouping variables
%    are sorted by the order of their levels, and numeric grouping variables
%    are sorted in numeric order.  To control the order of the groups,
%    you can either use categorical variables in G and specify the order of
%    their levels, or use the 'positions' argument.
%
%    BOXPLOT(AX, X, ...) produces a box plot in axes with handle AX.
%
%    BOXPLOT(..., 'PARAM1', val1, 'PARAM2', val2, ...) specifies optional
%    parameter name/value pairs.
%      'plotstyle'     'traditional' (default), or 'compact' to specify a
%                      box style designed for plots with many groups.  The
%                      plotstyle changes the defaults for some other
%                      parameters, as described below.
%
%      'boxstyle'      'outline' (default) to draw an unfilled box with
%                      dashed lines for whiskers, or 'filled' to draw a
%                      narrow filled box with solid lines for whiskers.
%      'colorgroup'    One or more grouping variables, of the same type as
%                      permitted for G, specifying that the box color should
%                      change when the specified variables change.  Default
%                      is [] for no box color change.
%      'colors'        Colors for boxes, specified as a single color (such
%                      as 'r' or [1 0 0]) or multiple colors (such as 'rgbm'
%                      or a three-column matrix of RGB values).  The sequence
%                      is replicated or truncated as required, so for example
%                      'rb' gives boxes that alternate in color.  Default
%                      when no 'colorgroup' is specified is to use the same
%                      color scheme for all boxes.  Default with
%                      'colorgroup' is a modified hsv colormap.
%      'datalim'       A two-element vector containing lower and upper limits,
%                      used by 'extrememode' to determine which points are
%                      extreme.  Default is [-Inf Inf].
%      'extrememode'   'clip' (default) to move data outside the 'datalim'
%                      limits to the limit, or 'compress' to distribute such
%                      points evenly in a region just outside the limit,
%                      retaining the relative order of the points.  A
%                      dotted line marks the limit if any points are outside
%                      it, and two gray lines mark the compression region if
%                      any points are compressed.  Values at +/-Inf can be
%                      clipped or compressed, but NaNs still do not appear
%                      on the plot.  Box notches are drawn to scale and may
%                      extend beyond the bounds if the median is inside the
%                      limit; they are not drawn if the median is outside
%                      the limits.
%      'factordirection' 'data' (default) to arrange the factors with the
%                      first value next to the origin, 'list' to arrange the
%                      factors left-to-right if on the x axis or top-to-
%                      bottom if on the y axis, or 'auto' to use 'data' for
%                      numeric grouping variables and 'list' for strings.
%      'fullfactors'   'off' (default) to have one group for each unique row
%                      of G, or 'on' to create a group for each possible
%                      combination of group variable values, including
%                      combinations that do not appear in the data.
%      'factorseparator' Specifies which factors should have their values
%                      separated by a grid line.  The value may be 'auto' or
%                      a vector of grouping variable numbers.  For example,
%                      [1 2] adds a separator line when the first or second
%                      grouping variable changes value.  'auto' is [] for
%                      one grouping variable and [1] for two or more
%                      grouping variables. Default is [].
%      'factorgap'     Specifies an extra gap to leave between boxes when
%                      the corresponding grouping factor changes value,
%                      expressed as a percentage of the width of the plot.
%                      For example, with [3 1], the gap is 3% of the width
%                      of the plot between groups with different values of
%                      the first grouping variable, and 1% between groups
%                      with the same value of the first grouping variable
%                      but different values for the second.  'auto'
%                      specifies that BOXPLOT should choose a gap
%                      automatically.  Default is [].
%      'grouporder'    Order of groups for plotting, specified as a cell
%                      array of strings.  With multiple grouping variables,
%                      separate values within each string with a comma.
%                      Using categorical arrays as grouping variables is an
%                      easier way to control the order of the boxes.
%      'jitter'        Maximum distance D to displace outliers along the
%                      factor axis by a uniform random amount, in order to
%                      make duplicate points visible.  D = 1 makes the
%                      jitter regions just touch between the closest
%                      adjacent groups.  The default is 0.
%      'labels'        Character array, cell array of strings, or numeric
%                      vector of box labels.  May have one label per group
%                      or per X value.  Multiple label variables may be
%                      specified via a numeric matrix or a cell array
%                      containing any of these types.
%      'labelorientation' 'horizontal' (default) for horizontal labels, or
%                      'inline' to draw the labels vertically when
%                      'orientation' has its default 'vertical' value.
%      'labelverbosity'  'all' (default) to display every label, 'minor' to
%                      display a label for a factor only when that factor
%                      has a different value from the previous group, or
%                      'majorminor' to display a label for a factor when
%                      that factor or any factor major to it has a
%                      different value from the previous group.
%      'medianstyle'   'line' (default) to draw a line for the median, or
%                      'target' to draw a black dot inside a white circle.
%      'notch'         'on' to draw comparison intervals using notches
%                      ('plotstyle' is 'traditional) or triangular markers
%                      ('plotstyle' is 'compact'), 'marker' to draw them
%                      using triangular markers, or 'off' (default) to omit
%                      them.  Two medians are significantly different at the
%                      5% level if their intervals do not overlap.  The
%                      interval endpoints are the extremes of the notches or
%                      the centers of the triangular markers.  When the
%                      sample size is small, notches may extend beyond the
%                      end of the box.
%      'orientation'   'vertical' (default) to plot X on the y axis, or
%                      'horizontal' to plot X on the x axis.
%      'outliersize'   Size of marker used for outliers, in points.
%                      Default is 6.
%      'positions'     Box positions specified as a numeric vector with one
%                      entry per group or X value (default 1:NGROUPS when the
%                      number of groups is NGROUPS).
%      'symbol'        Symbol and color to use for outliers, using the same
%                      values as the LineSpec parameter S in PLOT.  Default
%                      is 'r+'. If the symbol is omitted then the outliers
%                      are invisible; if the color is omitted then the
%                      outliers have the same color as their corresponding
%                      box.  Any line specification in S is ignored.
%      'whisker'       Maximum whisker length W.  Default is W=1.5.  Points
%                      are drawn as outliers if they are larger than
%                      Q3+W*(Q3-Q1) or smaller than Q1-W*(Q3-Q1), where Q1
%                      and Q3 are the 25th and 75th percentiles, respectively.
%                      The default value 1.5 corresponds to approximately +/-
%                      2.7 sigma and 99.3 coverage if the data are normally
%                      distributed.  The plotted whisker extends to the
%                      adjacent value, which is the most extreme data value
%                      that is not an outlier. Set 'whisker' to 0 to give no
%                      whiskers and to make every point outside of Q1 and Q3
%                      an outlier.
%      'widths'        A scalar or vector of box widths to use when the
%                      'boxstyle' is 'outline'.  The default is half of the
%                      minimum separation between boxes, which is .5 when
%                      the 'positions' argument takes its default value.
%                      The list of values is replicated or truncated as
%                      necessary.
%
%    When the 'plotstyle' parameter takes the value 'compact', then the
%    default values for other parameters are the following:
%        boxstyle - 'filled'            labelverbosity - 'majorminor'
%        factorgap - 'auto'             medianstyle - 'target'
%        factorseparator - 'auto'       outliersize - 4
%        jitter - 0.5                   symbol - 'o'
%        labelorientation - 'inline'
%
%    You can see the data values and group names by using the data cursor
%    tool, available from the figure window.  The data cursor shows the
%    original values of any points affected by the 'datalim' parameter.  You
%    can label the specific group to which an outlier belongs using the gname
%    function.
%
%    To modify the properties of box components, use findobj using tags to
%    find their handles as in one of the examples below.  The tag names
%    depend on the plotstyle and are:
%
%       all styles:  'Box', 'Outliers'
%       traditional: 'Median', 'Upper Whisker', 'Lower Whisker',
%                    'Upper Adjacent Value', 'Lower Adjacent Value',
%       compact:     'Whisker', 'MedianOuter', 'MedianInner'
%       when 'notch' is 'marker':
%                    'NotchLo', 'NotchHi'
%
%    Examples:
%       % Box plot of car gas mileage grouped by country
%       load carsmall
%       boxplot(MPG, Origin)
%       boxplot(MPG, Origin, 'sym','r*', 'colors',hsv(7))
%       boxplot(MPG, Origin, 'grouporder', ...
%                    {'France' 'Germany' 'Italy' 'Japan' 'Sweden' 'USA'})
%
%       % Plot by median gas mileage
%       [sortedMPG,sortedOrder] = sort(grpstats(MPG,Origin,@median));
%       pos(sortedOrder) = 1:6;
%       boxplot(MPG, Origin, 'positions', pos)
%
%       % Change some graphics properties
%       boxplot(chi2rnd(1,100,10)); % Generate box plot
%       h=findobj(gca,'tag','Outliers'); % Get handles for outlier lines.
%       set(h,'Marker','o'); % Change symbols for all the groups.
%       set(h(1),'MarkerEdgeColor','b'); % Change color for one group
%
%    See also ANOVA1, KRUSKALWALLIS, MULTCOMPARE.
%
%    Reference page in Doc Center
%       doc boxplot
%
%    Other functions named boxplot
%
%       ts/boxplot
%

author: jmaih


text:
function [M, RHS, Mx, My, Mz, RHSx, RHSy, RHSz] = convectionTvdTerm(u, phi, FL)
% This function uses the TVD scheme to discretize a
% convection term in the form $\grad (u \phi)$ where u is a face vactor
% It also returns the x, y, x parts of the matrix of coefficient.
%
% SYNOPSIS:
%   [M, RHS, Mx, My, Mz, RHSx, RHSy, RHSz] = convectionTvdTerm(u, phi, FL)
%
% PARAMETERS:
%   u  - velocity vector, FaceVariable
%   phi  - value of phi from the previous time step or iteration, CellVariable
%   FL  - Flux Limiter function
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

Mz=[];
d = u.domain.dimension;
switch d
    case 1
        [M, RHS] = convectionTvdTerm1D(u, phi, FL);
    case 1.5
        [M, RHS] = convectionTvdTermCylindrical1D(u, phi, FL);
    case 1.8
        [M, RHS] = convectionTvdTermSpherical1D(u, phi, FL);
    case 2
        [M, RHS, Mx, My, RHSx, RHSy] = convectionTvdTerm2D(u, phi, FL);
    case 2.5
        [M, RHS, Mx, My, RHSx, RHSy] = convectionTvdTermCylindrical2D(u, phi, FL);
    case 2.8
        [M, RHS, Mx, My, RHSx, RHSy] = ...
            convectionTvdTermRadial2D(u, phi, FL);
    case 3
        [M, RHS, Mx, My, Mz, RHSx, RHSy, RHSz] = convectionTvdTerm3D(u, phi, FL);
    case 3.2
        [M, RHS, Mx, My, Mz, RHSx, RHSy, RHSz] = ...
            convectionTvdTermCylindrical3D(u, phi, FL);
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function MS = createMeshCylindrical3D(varargin)
% MeshStructure = createMesh3D(Nr, Ntheta, Nz, Radius, theta, height)
% MeshStructure = createMesh3D(facelocationR, facelocationTheta, facelocationZ)
% creates a uniform 3D mesh:
% Nx is the number of cells in x (horizontal) direction
% Ny is the number of cells in y (vertical) direction
% Nz is the number of cells in z (perpendicular) direction
% Lx is the domain length in x direction
% Ly is the domain length in y direction
% Lz is the domain length in z direction
%
% SYNOPSIS:
%   MeshStructure = createMesh3D(Nx, Ny, Nz, Lx, Ly, Lz)
%
% PARAMETERS:
%   Nx: number of cells in the x direction
%   Lx: domain length in x direction
%   Ny: number of cells in the y direction
%   Ly: domain length in y direction
%   Nz: number of cells in the z direction
%   Lz: domain length in z direction
%
% RETURNS:
%   MeshStructure.
%                 dimensions=3 (3D problem)
%                 numbering: shows the indexes of cellsn from left to right
%                 and top to bottom and back to front
%                 cellsize: x, y, and z elements of the cell size =[Lx/Nx,
%                 Ly/Ny, Lz/Nz]
%                 cellcenters.x: location of each cell in the x direction
%                 cellcenters.y: location of each cell in the y direction
%                 cellcenters.z: location of each cell in the z direction
%                 facecenters.x: location of interface between cells in the
%                 x direction
%                 facecenters.y: location of interface between cells in the
%                 y direction
%                 facecenters.z: location of interface between cells in the
%                 z direction
%                 numberofcells: [Nx, Ny, Nz]
%
%
% EXAMPLE:
%   Nx = 2;
%   Lx = 1.0;
%   Ny = 3;
%   Ly = 2.0;
%   Nz = 4;
%   Lz = 3.0;
%   m = createMesh3D(Nx, Ny, Nz, Lx, Ly, Lz);
%   [X, Y, Z] = ndgrid(m.cellcenters.x, m.cellcenters.y, m.cellcenters.z);
%   [Xf, Yf, Zf] = ndgrid(m.facecenters.x, m.facecenters.y, m.facecenters.z);
%   plot3(X(:), Y(:), Z(:), 'or')
%   hold on;
%   plot3(Xf(:), Yf(:), Zf(:), '+b')
%   legend('cell centers', 'cell corners');
%
% SEE ALSO:
%     createMesh1D, createMesh2D, createMeshCylindrical1D, ...
%     createMeshCylindrical2D, createCellVariable, createFaceVariable

% Written by Ali A. Eftekhari
% See the license file

if nargin==6
  % uniform 1D mesh
  Nx=varargin{1};
  Ny=varargin{2};
  Nz=varargin{3};
  Width=varargin{4};
  Depth=varargin{6};
  Tetta=varargin{5};
  if Tetta>2*pi
      warning('Tetta is higher than 2*pi. It is scaled to 2*pi');
      Tetta = 2*pi;
  end
  MS=createMesh3D(Nx,Ny,Nz,Width, Tetta, Depth);
elseif nargin==3
  % nonuniform 1D mesh
  facelocationX=varargin{1};
  facelocationTheta=varargin{2};
  facelocationZ=varargin{3};
  if facelocationTheta(end)>2*pi
      facelocationTheta = facelocationTheta/facelocationTheta(end)*2.0*pi;
      warning('The domain size adjusted to match a maximum of 2*pi.')
  end
  MS=createMesh3D(facelocationX,facelocationTheta,facelocationZ);
end
MS.dimension=3.2;


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2014 National Renewable Energy Laboratory and National
% Technology & Engineering Solutions of Sandia, LLC (NTESS).
% Under the terms of Contract DE-NA0003525 with NTESS,
% the U.S. Government retains certain rights in this software.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

classdef waveClass<handle
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % The  ``waveClass`` creates a ``waves`` object saved to the MATLAB
    % workspace. The ``waveClass`` includes properties and methods used
    % to define WEC-Sim's wave input.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    properties (SetAccess = 'public', GetAccess = 'public')%input file
        bem             = struct(...                % (`structure`) Defines the BEM data implementation.
            'option',       'EqualEnergy',...       %
            'count',        [], ...                 %
            'frequency',    [], ...                 %
            'range',        [])                     % (`structure`) Defines the BEM data implementation. ``option`` (`string`) Method of frequency discretization for irregular waves, options include: ``'EqualEnergy'`` or ``'Traditional'``. Default = ``'EqualEnergy'``. ``count`` (`integer`) Number of interpolated wave frequencies, only used for ``irregular`` and ``spectrumImport``. Number of frequencies used varies depending on ``bem.option``, 1000 for ``'Traditional'``, and 500 for ``'EqualEnergy'`` and ``Imported``. Default = ``[]``. ``frequency`` (`vector`) Wave frequencies [rad/s] from BEM. Default = ``[]``. ``range`` (`2x1 vector`) Min and max wave frequency [rad/s], only used for ``irregular`` and ``spectrumImport``. If not specified, the BEM data frequency range is used. Default = ``[]``
        current         = struct(...                % (`structure`) Defines the current implementation.
            'option',       3,...                   %
            'depth',        0, ...                  %
            'direction',    0, ...                  %
            'speed',        0)                      % (`structure`) Defines the current implementation. ``option`` (`integer`) Define the sub-surface current model to be used in WEC-Sim, options include: ``0`` for depth-independent model, ``1`` for 1/7 power law variation with depth, ``2`` for linear variation with depth, or ``3`` for no current. Default = ``3``, ``depth`` (`float`) Current depth [m]. Define the depth over which the sub-surface current is modeled. Must be defined for options ``1`` and ``2``. The current is not calculated for any depths greater than the specified current depth. Default = ``0``, ``direction`` (`float`) Current direction [deg]. Surface current direction defined using WEC-Sim global coordinate system. Default = ``0``, ``speed``  (`float`) Current seed [m/s]. Surface current speed that is uniform along the water column. Default = ``0``
        direction       = 0;                        % (`float`) Incident wave direction(s) [deg]. Incident wave direction defined using WEC-Sim global coordinate system. Should be defined as a column vector for more than one wave direction. Default = ``0``
        elevationFile   = 'NOT DEFINED';            % (`string`) Data file that contains the times-series data file. Default = ``'NOT DEFINED'``
        gamma           = [];                       % (`float`) Defines gamma, only used for ``JS`` wave spectrum type. Default = ``[]``
        height          = 'NOT DEFINED';            % (`float`) Wave height [m]. Defined as wave height for ``regular``, or significant wave height for ``irregular``. Default =  ``'NOT DEFINED'``
        marker          = struct(...                % (`structure`) Defines the wave marker.
            'location',     [],...                  %
            'size',         10, ...                 %
            'style',        1)                      % (`structure`) Defines the wave marker. `loc` (`nx2 vector`) Marker [X,Y] locations [m]. Default = ``[]``. ``size`` (`float`) Marker size in Pixels. Default = ``10``. ``style`` Marker style, options include: ``1``: Sphere, ``2``: Cube, ``3``: Frame. Default = ``1``: Sphere
        period          = 'NOT DEFINED';            % (`float`) Wave period [s] . Defined as wave period for ``regular``, peak period for ``irregular``, or period of BEM data used for hydrodynamic coefficients for ``noWave``. Default = ``'NOT DEFINED'``
        phaseSeed       = 0;                        % (`integer`) Defines the random phase seed, only used for ``irregular`` and ``spectrumImport`` waves. Default = ``0``
        spectrumFile    = 'NOT DEFINED';            % (`string`) Data file that contains the spectrum data file.  Default = ``'NOT DEFINED'``
        spectrumType    = 'NOT DEFINED';            % (`string`) Specifies the wave spectrum type, options inlcude:``PM`` or ``JS``. Default = ``'NOT DEFINED'``
        viz             = struct(...                %
            'numPointsX',   50, ...                 %
            'numPointsY',   50 );                   % (`structure`) Defines visualization options, structure contains the fields ``numPointsX`` for the number of visualization points in x direction, and ``numPointsY`` for the number of visualization points in y direction.
        waterDepth      = [];                       % (`float`) Water depth [m]. Default to BEM water depth if not set.
        spread          = 1;                        % (`float`) Wave Spread probability associated with wave direction(s). Should be defined as a column vector for more than one wave direction. Default = ``1``
    end

    properties (SetAccess = 'private', GetAccess = 'public')%internal
        % The following properties are private, for internal use by WEC-Sim
        amplitude       = [];                       % Wave amplitude [m]. For regular waves or 2*(wave spectrum vector) for irregular waves
        deepWater       = [];                       % Deep water or not, depending on input from WAMIT, NEMOH and AQWA
        dOmega          = 0;                        % Frequency spacing [rad] for ``irregular`` waves.
        omega           = [];                       % Wave frequency (regular waves) or wave frequency vector (irregular waves), where omega = 2*pi/period [rad/s]
        phase           = 0;                        % (`float`) Wave phase [rad] . Only used for ``irregular`` waves.
        power           = [];                       % Wave Power Per Unit Wave Crest [W/m]
        spectrum        = [];                       % Wave Spectrum [m^2-s/rad] for ``Traditional``
        type            = 'NOT DEFINED';            % (`string`) Specifies the wave type, options include:``noWave``, ``noWaveCIC``, ``regular``, ``regularCIC``, ``irregular``, ``spectrumImport``, or ``elevationImport``. Default = ``'NOT DEFINED'``
        typeNum         = [];                       % Number to represent different type of waves
        waveAmpTime     = [];                       % Wave elevation time history [m] at the (0, 0, 0) origin
        waveAmpTimeViz  = [];                       % Wave elevation time history at marker locations specified by user [m]
        wavenumber      = [];                       % Wave Number, where wavenumber = 2*pi/wavelength [rad/m]
    end

    methods (Access = 'public')
        function obj = waveClass(type)
            % This method initializes the ``waveClass`` and creates a ``waves`` object.
            %
            % Parameters
            % ------------
            %     type : string
            %         String specifying the wave type, options include:
            %
            %           noWave
            %               No waves
            %
            %           noWaveCIC
            %               No waves but with the Convolution Integral Calculation to calculate radiation effects
            %
            %           regular
            %               Regular waves
            %
            %           regularCIC
            %               Regular waves with Convolution Integral Calculation to calculate radiation effects
            %
            %           irregular
            %           	Irregular Waves
            %
            %           spectrumImport
            %           	Irregular Waves with predefined phase
            %
            %           elevationImport
            %               Irregular Waves with predefined elevation
            %
            % Returns
            % ------------
            %     waves : obj
            %         waveClass object
            %

            obj.type = type;
            switch obj.type
                case {'noWave'}         % No Waves with Constant Hydrodynamic Coefficients
                    obj.typeNum = 0;
                case {'noWaveCIC'}      % No Waves w/Convolution Integral Calculation
                    obj.typeNum = 1;
                case {'regular'}        % Regular Waves with Constant Hydrodynamic Coefficients
                    obj.typeNum = 10;
                case {'regularCIC'}     % Regular Waves w/Convolution Integral Calculation
                    obj.typeNum = 11;
                case 'irregular'        % Irregular Waves with 'PM' or 'JS' wave spectrum
                    obj.typeNum = 20;
                case 'spectrumImport'   % Irregular waves with imported wave spectrum
                    obj.typeNum = 21;
                case 'elevationImport'        % Waves with imported wave elevation time-history
                    obj.typeNum = 30;
            end
        end

        function checkInputs(obj)
            % This method checks WEC-Sim user inputs and generates error messages if parameters are not properly defined.

            % check wave type
            types = {'noWave', 'noWaveCIC', 'regular', 'regularCIC', 'irregular', 'spectrumImport', 'elevationImport'};
            if sum(strcmp(types,obj.type)) ~= 1
                error(['Unexpected wave environment type setting, choose from: ' ...
                    '"noWave", "noWaveCIC", "regular", "regularCIC", "irregular", "spectrumImport", and "elevationImport".'])
            end
            % check 'waves.bem' fields
            if length(fieldnames(obj.bem)) ~=4
                error(['Unrecognized method, property, or field for class "waveClass", ' ...
                    '"waveClass.bem" structure must only include fields: "option", "count", "frequency", "range"']);
            end
            % check 'waves.current' fields
            if length(fieldnames(obj.current)) ~=4
                error(['Unrecognized method, property, or field for class "waveClass", ' ...
                    '"waveClass.current" structure must only include fields: "option", "depth", "direction", "speed"']);
            end
            % 'elevationFile' defined for 'elevationImport' case
            if strcmp(obj.type,'elevationImport')
                if strcmp(obj.elevationFile,'NOT DEFINED')
                    error('The "waves.elevationFile" must be defined when using the "elevationImport" wave type');
                end
                if ~strcmp(obj.period,'NOT DEFINED') || ~strcmp(obj.height,'NOT DEFINED')
                    warning('"waves.period" and "waves.height" are not used for "etaImport" wave types')
                end
            end
            % 'spectrumFile' defined for 'spectrumImport' case
            if strcmp(obj.type,'spectrumImport')
                if strcmp(obj.spectrumFile,'NOT DEFINED')
                    error('The "wave.spectrumFile" must be defined when using the "spectrumImport" wave type');
                end
                if ~strcmp(obj.period,'NOT DEFINED') || ~strcmp(obj.height,'NOT DEFINED')
                    warning('"waves.period" and "waves.height" are not used for "spectrumImport" wave types')
                end
            end

            % Check 'marker.location'
            if ~isempty(obj.marker.location) && ~ndims(obj.marker.location)==2
                error('The coordinates of the visualization markers should have an ordinate (y-coordinate) and an abscissa (x-coordinate)')
            end
            % Check wave spread
            if sum(obj.spread)~=1
                error('The wave spread should always sum to 1 to preserve spectrum/energy accuracy.')
            end

            % Check inputs based on type
            if strcmp(obj.type,'noWave') && strcmp(obj.period,'NOT DEFINED')
                error('"waves.period" must be defined for the hydrodynamic data period when using the "noWave" wave type');
            end
            if strcmp(obj.type,'noWaveCIC') && ~strcmp(obj.period,'NOT DEFINED')
                warning('"waves.period" is not used for the hydrodynamic data period when using the "noWaveCIC" wave type')
            end
            if strcmp(obj.type,'regular') || strcmp(obj.type,'regularCIC')
                if strcmp(obj.period,'NOT DEFINED') || strcmp(obj.height,'NOT DEFINED')
                    error('"waves.period" and "waves.height" need to be defined for "regular" and "regularCIC" wave types')
                end
                if ~strcmp(obj.spectrumType,'NOT DEFINED')
                    warning('"waves.spectrumType" is not used for the "regular" and "regularCIC" wave types')
                end
            end
            if strcmp(obj.type,'irregular')
                if strcmp(obj.spectrumType,'NOT DEFINED')
                    error('"waves.spectrumType" needs to be defined for "irregular" wave types')
                end
                if strcmp(obj.period,'NOT DEFINED') || strcmp(obj.height,'NOT DEFINED')
                    error('"waves.period" and "waves.height" need to be defined for "irregular" wave types')
                end
            end
        end

        function setup(obj,bemFreq,bemWaterDepth,rampTime,dt,maxIt,time,g,rho)
            % This method calculates WEC-Sim's wave properties
            % based on the specified wave type.

            if ~isempty(obj.marker.location)==1
                if ~width(obj.marker.location)==2
                    error('The coordinates of the visualization markers should have an ordinate (y-coordinate) and an abscissa (x-coordinate)')
                end
            end
            obj.bem.frequency = bemFreq;
            if isempty(obj.bem.range) && isempty(obj.bem.frequency)
                % No .h5 file and no frequency range defined --> error
                error('Must define frequency range in waves.bem.range when zero hydro bodies are used (no .h5 file).');
            elseif isempty(obj.bem.range) && ~isempty(obj.bem.frequency)
                % Use .h5 file range if not defined in input file
                obj.bem.range = obj.bem.frequency;
            elseif ~isempty(obj.bem.range) && ~isempty(obj.bem.frequency)
                % check that input file frequency range is not larger than
                % available BEM data
                if obj.bem.range(1) < min(obj.bem.frequency) || obj.bem.range(1) > max(obj.bem.frequency)
                    warning('Min frequency range outside BEM data, min frequency set to min BEM frequency')
                    obj.bem.range(1) = min(obj.bem.frequency);
                end
                if obj.bem.range(2) < min(obj.bem.frequency) || obj.bem.range(2) > max(obj.bem.frequency)
                    warning('Max frequency range outside BEM data, max frequency set to max BEM frequency')
                    obj.bem.range(2) = max(obj.bem.frequency);
                end
            end

            obj.setWaterDepth(bemWaterDepth);

            switch obj.type
                case {'noWave','noWaveCIC'}
                    if isempty(obj.omega) && strcmp(obj.period,'NOT DEFINED')
                        obj.omega = min(obj.bem.range);
                        obj.period = 2*pi/obj.omega;
                    elseif isempty(obj.omega)
                        obj.omega = 2*pi/obj.period;
                    else
                        obj.period = 2*pi/obj.omega;
                    end
                    obj.height = 0;
                    obj.amplitude = obj.height/2;
                    obj.wavenumber = calcWaveNumber(obj.omega,obj.waterDepth,g,obj.deepWater);
                    obj.waveElevNowave(time);
                case {'regular','regularCIC'}
                    if isempty(obj.omega) && strcmp(obj.period,'NOT DEFINED')
                        obj.omega = min(obj.bem.range);
                        obj.period = 2*pi/obj.omega;
                    elseif isempty(obj.omega)
                        obj.omega = 2*pi/obj.period;
                    else
                        obj.period = 2*pi/obj.omega;
                    end
                    obj.amplitude = obj.height/2;
                    obj.wavenumber = calcWaveNumber(obj.omega,obj.waterDepth,g,obj.deepWater);
                    obj.waveElevReg(rampTime, time);
                    obj.wavePowerReg(g,rho);
                case {'irregular','spectrumImport'}
                    if strcmp(obj.type,'spectrumImport')
                        obj.height = 0;
                        obj.period = 0;
                        obj.bem.option = 'Imported';
                        obj.spectrumType = 'spectrumImport';
                    end
                    minFrequency=min(obj.bem.range);
                    maxFrequency=max(obj.bem.range);
                    switch obj.bem.option
                        case {'Traditional'}
                            if isempty(obj.bem.count)
                                obj.bem.count = 1000;
                            end
                            obj.omega = (minFrequency:(maxFrequency-minFrequency)/(obj.bem.count-1):maxFrequency)';
                            obj.dOmega = ones(obj.bem.count,1).*(maxFrequency-minFrequency)./(obj.bem.count-1);
                        case {'EqualEnergy'}
                            bemCount_interp = 500000;
                            obj.omega = (minFrequency:(maxFrequency-minFrequency)/bemCount_interp:maxFrequency)';
                            obj.dOmega = mean(diff(obj.omega));
                            if isempty(obj.bem.count)
                                obj.bem.count = 500;
                            end
                        case {'Imported'}
                            data = importdata(obj.spectrumFile);
                            freqData = data(:,1);
                            freqLoc = freqData >= min(obj.bem.range)/2/pi & freqData <= max(obj.bem.range)/2/pi;
                            obj.omega    = freqData(freqLoc).*2.*pi;
                            obj.bem.count = length(obj.omega);
                            obj.dOmega(1,1)= obj.omega(2)-obj.omega(1);
                            obj.dOmega(2:obj.bem.count-1,1)=(obj.omega(3:end)-obj.omega(1:end-2))/2;
                            obj.dOmega(obj.bem.count,1)= obj.omega(end)-obj.omega(end-1);
                    end
                    obj.setWavePhase;
                    obj.irregWaveSpectrum(g,rho)
                    obj.wavenumber = calcWaveNumber(obj.omega,obj.waterDepth,g,obj.deepWater);
                    obj.waveElevIrreg(rampTime, time, obj.dOmega);
                case {'elevationImport'}    %  This does not account for wave direction
                    % Import 'elevationImport' time-series here and interpolate
                    data = importdata(obj.elevationFile) ;    % Import time-series
                    obj.waveElevUser(rampTime, maxIt, data, time);
            end
        end

        function listInfo(obj)
            % This method prints wave information to the MATLAB Command Window.

            fprintf('\nWave Environment: \n')
            switch obj.type
                case 'noWave'
                    fprintf('\tWave Type                            = No Wave (Constant Hydrodynamic Coefficients)\n')
                    fprintf('\tHydro Data Wave Period, period (sec)    	= %G\n',obj.period)
                case 'regular'
                    fprintf('\tWave Type                            = Regular Waves (Constant Hydrodynamic Coefficients)\n')
                    fprintf('\tWave Height, height (m)                   = %G\n',obj.height)
                    fprintf('\tWave Period, period (sec)                 = %G\n',obj.period)
                case 'noWaveCIC'
                    fprintf('\tWave Type                            = No Wave (Convolution Integral Calculation)\n')
                case 'regularCIC'
                    fprintf('\tWave Type                            = Regular Waves (Convolution Integral Calculation)\n')
                    fprintf('\tWave Height, height (m)                   = %G\n',obj.height)
                    fprintf('\tWave Period, period (sec)                 = %G\n',obj.period)
                case 'irregular'
                    if obj.phaseSeed == 0
                        fprintf('\tWave Type                            = Irregular Waves (Arbitrary Random Phase)\n')
                    else
                        fprintf('\tWave Type                            = Irregular Waves (Predefined Random Phase)\n')
                    end
                    obj.printWaveSpectrumType;
                    fprintf('\tSignificant Wave Height, Hs      (m) = %G\n',obj.height)
                    fprintf('\tPeak Wave Period, Tp           (sec) = %G\n',obj.period)
                case 'spectrumImport'
                    if size(importdata(obj.spectrumFile),2) == 3
                        fprintf('\tWave Type                            = Irregular waves with imported wave spectrum (Imported Phase)\n')
                    elseif obj.phaseSeed == 0
                        fprintf('\tWave Type                            = Irregular waves with imported wave spectrum (Random Phase)\n')
                    else
                        fprintf('\tWave Type                            = Irregular waves with imported wave spectrum (Seeded Phase)\n')
                    end
                    obj.printWaveSpectrumType;
                case 'elevationImport'
                    fprintf( '\tWave Type                           = Waves with imported wave elevation time-history\n')
                    fprintf(['\tWave Elevation Time-Series File    	= ' obj.elevationFile '  \n'])
            end
        end

        function calculateElevation(obj,rampTime,timeseries)
            % Calculates the wave elevation based on the wave type.
            % Used by postProcess.m
            %
            % Parameters
            % ------------
            %     waves : obj
            %         waveClass object
            %
            %     rampTime : float
            %         ramp time from the simulation class
            %
            %     timeseries : float array
            %         Array of all simulation output time steps
            %
            switch obj.type
                case {'noWave','noWaveCIC'}
                    obj.waveElevNowave(timeseries);

                case {'regular','regularCIC'}
                    obj.waveElevReg(rampTime, timeseries);

                case {'irregular','spectrumImport'}
                    obj.waveElevIrreg(rampTime, timeseries, obj.dOmega);

                case {'elevationImport'}
                    % Wave elevation is a necessary pre-processing step for
                    % the eta import case.
                    % Used by waveClass.setup()
            end
        end

        function Z = waveElevationGrid(obj, t, X, Y, it, g)
            % This method calculates wave elevation on a grid at a given
            % time, used by: :func:`write_paraview_wave`.
            %
            % Parameters
            % ------------
            %     waves: obj
            %         waveClass object
            %
            %     t : float
            %         the current time
            %
            %     X : matrix
            %         (m x n) matrix of X coordinates at which to calculate the wave elevation
            %
            %     Y : matrix
            %         (m x n) matrix of Y coordinates at which to calculate the wave elevation
            %
            %     it : time step iteration
            %
            %     g : gravitational acceleration constant from simulationClass
            %
            %
            % Returns
            % ---------
            %     Z : matrix
            %         (m x n) matrix of Z coordinates of the wave elevation
            %
            switch obj.type
                case {'noWave','noWaveCIC'}
                    Z = zeros (size (X));
                case {'regular', 'regularCIC'}
                    Xt = X*cos (obj.direction*pi/180) + Y * sin(obj.direction*pi/180);
                    Z = obj.amplitude * cos(-1 * obj.wavenumber * Xt  +  obj.omega * t);
                case {'irregular', 'spectrumImport'}
                    Z = zeros (size (X));
                    for idir=1:length(obj.direction)
                        Xt = X*cos(obj.direction(idir)*pi/180) + Y*sin(obj.direction(idir)*pi/180);
                        for iw = 1:length(obj.omega)
                            Z = Z + sqrt(obj.amplitude(iw)*obj.spread(idir).*obj.dOmega(iw)) * cos(-1*obj.wavenumber(iw)*Xt + obj.omega(iw)*t + obj.phase(iw,idir));
                        end
                    end
                case{'elevationImport'}
                    if it ==1
                        warning('Paraview wave surface discretization for qualitative purposes only.')
                    end
                    Z = zeros(size(X));
                    WaveEle = obj.waveAmpTime(:,2);
                    TimeWaveEle = obj.waveAmpTime(:,1);
                    L = length(TimeWaveEle);
                    Fs = 1/(TimeWaveEle(2)-TimeWaveEle(1));
                    YY = fft(WaveEle);
                    P2 = abs(YY/L);
                    P1 = P2(1:round(L/2)+1);
                    P1(2:end-1) = 2*P1(2:end-1);
                    ff = Fs*(0:round(L/2))/L;
                    PhasesEtaImp2 = angle(YY/L);
                    PhasesEtaImp1 = PhasesEtaImp2(1:round(L/2)+1);
                    kWaveEle = (2*pi.*ff).^2./g; % deep water wave approximation
                    Xt = X*cos(obj.direction*pi/180) + Y*sin(obj.direction*pi/180);
                    Zint=zeros(size(X));
                    for iw = 2:length(ff)
                        Zint = Zint + P1(iw) * cos(-1*kWaveEle(iw)*Xt + 2*pi*ff(iw).*t+ PhasesEtaImp1(iw));
                    end
                    Z = Zint;
            end
        end

        function plotElevation(obj,rampTime)
            % This method plots wave elevation time-history.
            %
            % Parameters
            % ------------
            %     rampTime : float, optional
            %         Specify wave ramp time to include in plot
            %
            % Returns
            % ------------
            %     figure : fig
            %         Plot of wave elevation versus time
            %

            arguments
                obj
                rampTime double {mustBeReal, mustBeNonNan, mustBeFinite} = 0
            end

            figure
            plot(obj.waveAmpTime(:,1),obj.waveAmpTime(:,2))
            title('Wave Surfave Elevation')
            if nargin==2
                hold on
                line([rampTime,rampTime],[1.5*min(obj.waveAmpTime(:,2)),1.5*max(obj.waveAmpTime(:,2))],'Color','k')
                title(['Wave Surface Elevation, Ramp Time ' num2str(rampTime) ' (s)'])
            end
            xlabel('Time (s)')
            ylabel('Elevation (m)')
        end

        function plotSpectrum(obj)
            % This method plots the wave spectrum.
            %
            % Returns
            % ------------
            %     figure : fig
            %         Plot of wave spectrum versus wave frequency
            %
            m0 = trapz(obj.omega,obj.spectrum);
            HsTest = 4*sqrt(m0);
            [~,I] = max(abs(obj.spectrum));
            wp = obj.omega(I);
            TpTest = 2*pi/wp;

            figure
            plot(obj.omega,obj.spectrum,'s-')
            hold on
            line([wp,wp],[0,max(obj.spectrum)],'Color','k')
            xlim([0 max(obj.omega)])
            title([obj.spectrumType, ' Spectrum, T_p= ' num2str(TpTest) ' [s], '  'H_m_0= ' num2str(HsTest), ' [m]'])
            if strcmp(obj.spectrumType,'JS') == 1
                title([obj.spectrumType, ' Spectrum, T_p= ' num2str(TpTest) ' [s], H_m_0= ' num2str(HsTest), ' [m], gamma = ' num2str(obj.gamma)])
            end
            xlabel('Frequency (rad/s)')
            ylabel('Spectrum (m^2-s/rad)');
        end

    end

    methods (Access = 'protected')
        function setWavePhase(obj)
            % Sets the irregular wave's random phase
            % used by: :meth:`waveClass.setup`.
            if obj.phaseSeed ~= 0
                rng(obj.phaseSeed);     % Phase seed = 1,2,3,...,etc
            else
                rng('shuffle');         % Phase seed shuffled
            end
            switch obj.bem.option
                case {'EqualEnergy','Traditional'}
                    obj.phase = 2*pi*rand(length(obj.direction),obj.bem.count);
                case {'Imported'}
                    data = importdata(obj.spectrumFile);
                    if size(data,2) == 3
                        freqData = data(:,1);
                        freqLoc = freqData>=min(obj.bem.range)/2/pi & freqData<=max(obj.bem.range)/2/pi;
                        phaseData = data(freqLoc,3);
                        obj.phase = phaseData';
                    else
                        obj.phase = 2*pi*rand(1,obj.bem.count);
                    end
            end
            obj.phase = obj.phase';
        end

        function setWaterDepth(obj,bemWaterDepth)
            % Set the water depth. If defined in input file, BEM depth is
            % not used. used by: :meth:`waveClass.setup`.
            if isempty(obj.waterDepth)
                if isempty(bemWaterDepth)
                    error('Must define water depth in waves.waterDepth when zero hydro bodies (no .h5 file).');
                elseif strcmp(bemWaterDepth,'infinite')
                    obj.deepWater = 1;
                    obj.waterDepth = 200;
                    fprintf('\tInfinite water depth specified in BEM and "waves.waterDepth" not specified in input file.\n')
                    fprintf('Set water depth to 200m for visualization.\n')
                else
                    obj.deepWater = 0;
                    obj.waterDepth = double(bemWaterDepth);
                end
            else
                if ~isempty(bemWaterDepth)
                    warning('Because water depth is specified in the wecSimInputFile, the water depth from the BEM data is ignored')
                end
            end
        end

        function waveElevNowave(obj,timeseries)
            % Set noWave elevation time-history
            % used by: :meth:`waveClass.setup`.
            % Used by postProcess for variable step solvers
            obj.waveAmpTime         = zeros(length(timeseries),2);
            obj.waveAmpTime(:,1)    = timeseries;

            if ~isempty(obj.marker.location)
                SZwaveAmpTimeViz = size(obj.marker.location);
                obj.waveAmpTimeViz = zeros(length(timeseries),SZwaveAmpTimeViz(1)+1);
                obj.waveAmpTimeViz(:,1) = timeseries;
            end
        end

        function waveElevReg(obj,rampTime,timeseries)
            % Calculate regular wave elevation time history
            % used by: :meth:`waveClass.setup`.
            % Used by postProcess for variable step solvers
            maxIt = length(timeseries);
            rampFunction = (1+cos(pi+pi*timeseries/rampTime))/2;
            rampFunction(timeseries>=rampTime) = 1;

            obj.waveAmpTime = zeros(maxIt,2);
            obj.waveAmpTime(:,1) = timeseries;
            obj.waveAmpTime(:,2) = rampFunction.*(obj.amplitude*cos(obj.omega*timeseries));

            % Wave Marker
            if ~isempty(obj.marker.location)==1
                if width(obj.marker.location)~=2
                    error('The coordinates of the visualization markers should have an ordinate (y-coordinate) and an abscissa (x-coordinate)')
                end
            end
            if ~isempty(obj.marker.location)
                SZwaveAmpTimeViz = size(obj.marker.location);
                obj.waveAmpTimeViz = zeros(maxIt,SZwaveAmpTimeViz(1)+1);
                for j = 1:SZwaveAmpTimeViz(1)
                    obj.waveAmpTimeViz(:,1) = timeseries;
                    obj.waveAmpTimeViz(:,j+1) = rampFunction.*obj.amplitude.*cos(obj.omega*timeseries - obj.wavenumber*(obj.marker.location(j,1).*cos(obj.direction*pi/180) + obj.marker.location(j,2).*sin(obj.direction*pi/180)));
                end
            end
        end

        function wavePowerReg(obj,g,rho)
            % Calculate wave power per unit wave crest for regular waves
            if obj.deepWater == 1
                % Deepwater Approximation
                obj.power = 1/(8*pi)*rho*g^(2)*(obj.amplitude).^(2).*obj.period;
            else
                % Full Wave Power Equation
                obj.power = rho*g*(obj.amplitude).^(2)/4*sqrt(g./obj.wavenumber.*tanh(obj.wavenumber.*obj.waterDepth))*(1+2*obj.wavenumber.*obj.waterDepth./sinh(2*obj.wavenumber.*obj.waterDepth));
            end
        end

        function irregWaveSpectrum(obj,g,rho)
            % Calculate wave spectrum vector (obj.amplitude)
            % Used by wavesIrreg (wavesIrreg used by: :meth:`waveClass.setup`.)
            frequency = obj.omega/(2*pi);
            Tp = obj.period;
            Hs = obj.height;
            switch obj.spectrumType
                case {'PM','JS'}
                    % Pierson-Moskowitz Spectrum from IEC TS 62600-2 ED2 Annex C.2 (2019)
                    bPM = (5/4)*(1/Tp)^(4);
                    aPM =  bPM*(Hs/2)^2;
                    fSpectrum  = (aPM*frequency.^(-5).*exp(-bPM*frequency.^(-4)));            % Wave Spectrum [m^2-s] for 'EqualEnergy'
                    if strcmp(obj.spectrumType,'JS')
                        % JONSWAP Spectrum from IEC TS 62600-2 ED2 Annex C.2 (2019)
                        fp = 1/Tp;
                        siga = 0.07;sigb = 0.09;                                    % cutoff frequencies for gamma function
                        [lind,~] = find(frequency<=fp);
                        [hind,~] = find(frequency>fp);
                        gammaAlpha = zeros(size(frequency));
                        if isempty(obj.gamma)
                            TpsqrtHs = Tp/sqrt(Hs);
                            if TpsqrtHs <= 3.6
                                obj.gamma = 5;
                            elseif TpsqrtHs > 5
                                obj.gamma = 1;
                            else
                                obj.gamma = exp(5.75 - 1.15*TpsqrtHs);
                            end
                        end
                        gammaAlpha(lind) = obj.gamma.^exp(-(frequency(lind)-fp).^2/(2*siga^2*fp^2));
                        gammaAlpha(hind) = obj.gamma.^exp(-(frequency(hind)-fp).^2/(2*sigb^2*fp^2));
                        C = 1 - 0.287*log(obj.gamma);
                        fSpectrum = C*fSpectrum.*gammaAlpha;                                % Wave Spectrum [m^2-s]
                    end
                    obj.spectrum = fSpectrum./(2*pi);                                       % Wave Spectrum [m^2-s/rad]
                case 'spectrumImport' % Imported Wave Spectrum
                    data = importdata(obj.spectrumFile);
                    freqData = data(:,1);
                    S_data = data(:,2);
                    freqLoc = freqData >= min(obj.bem.range)/2/pi & freqData <= max(obj.bem.range)/2/pi;
                    fSpectrum = S_data(freqLoc);                                    % Wave Spectrum [m^2-s] for 'EqualEnergy'
                    obj.spectrum = fSpectrum./(2*pi);                                       % Wave Spectrum [m^2-s/rad] for 'Traditional'
                    fprintf('\t"spectrumImport" uses the number of imported wave frequencies (not "Traditional" or "EqualEnergy")\n')
                case {'BS'}
                    error('Following IEC Standard, our Bretschneider Sprectrum (BS) option is exactly how the Pierson-Moskowitz (PM) Spectrum is defined. Please use PM instead');
            end
            % Power per Unit Wave Crest
            obj.wavenumber = calcWaveNumber(obj.omega,obj.waterDepth,g,obj.deepWater); %Calculate Wave Number for Larger Number of Frequencies Before Down Sampling in Equal Energy Method
            if obj.deepWater == 1
                % Deepwater Approximation
                obj.power = sum(1/2*rho*g^(2)*fSpectrum/(2*pi).*obj.dOmega./obj.omega);
            else
                % Full Wave Power Equation
                obj.power = sum((1/2)*rho*g.*fSpectrum/(2*pi).*obj.dOmega.*sqrt(g./obj.wavenumber.*tanh(obj.wavenumber.*obj.waterDepth)).*(1 + 2.*obj.wavenumber.*obj.waterDepth./sinh(2.*obj.wavenumber.*obj.waterDepth)));
            end
            %
            switch obj.bem.option
                case {'EqualEnergy'}
                    m0 = trapz(frequency,abs(fSpectrum));
                    numBins = obj.bem.count+1;
                    a_targ = m0/(numBins);
                    SfInt = cumtrapz(frequency,fSpectrum);
                    wn(1) = 1;
                    for kk = 1:numBins
                        jj = 1;
                        tmpa{kk}(1) = 0;
                        while tmpa{kk}(jj)-kk*a_targ < 0
                            tmpa{kk}(jj+1) = SfInt(wn(kk)+jj);
                            jj = jj+1;
                            if wn(kk)+jj >=length(fSpectrum)
                                break
                            end
                        end
                        [a_targ_real(kk),wna(kk)] = min(abs(tmpa{kk}-kk*a_targ));
                        wn(kk+1) = wna(kk)+wn(kk);
                        a_bins(kk) = trapz(frequency(wn(kk):wn(kk+1)),abs(fSpectrum(wn(kk):wn(kk+1))));
                    end
                    obj.omega = 2*pi*frequency(wn(2:end-1));
                    obj.dOmega = [obj.omega(1)-2*pi*frequency(wn(1)); diff(obj.omega)];
                    obj.spectrum = obj.spectrum(wn(2:end-1));                             % Wave Spectrum [m^2-s/rad]
            end
            obj.amplitude = 2 * obj.spectrum;                                              % Wave Amplitude [m]
        end

        function waveElevIrreg(obj,rampTime,timeseries,df)
            % Calculate irregular wave elevetaion time history
            % used by: :meth:`waveClass.setup`.
            % Used by postProcess for variable time step
            maxIt = length(timeseries);
            rampFunction = (1+cos(pi+pi*timeseries/rampTime))/2;
            rampFunction(timeseries>=rampTime) = 1;

            obj.waveAmpTime = zeros(maxIt,2);
            obj.waveAmpTime(:,1) = timeseries;

            %  Wave Markers
            if ~isempty(obj.marker.location)
                if width(obj.marker.location)~=2
                    error('The coordinates of the visualization markers should have an ordinate (y-coordinate) and an abscissa (x-coordinate)')
                end
            end
            if ~isempty(obj.marker.location)
               SZwaveAmpTimeViz = size(obj.marker.location);
               obj.waveAmpTimeViz = zeros(maxIt,SZwaveAmpTimeViz(1)+1);
               obj.waveAmpTimeViz(:,1) = timeseries;
            end

            % Calculate eta
            for i = 1:length(timeseries)
                tmp  = sqrt(obj.amplitude.*df*obj.spread);
                tmp1 = tmp.*real(exp(sqrt(-1).*(obj.omega.*timeseries(i) + obj.phase)));
                obj.waveAmpTime(i,2) = rampFunction(i)*sum(tmp1,'all');

                % Wave Markers
                if ~isempty(obj.marker.location)
                    for j = 1:SZwaveAmpTimeViz(1)
                        tmp14 = tmp.*real(exp(sqrt(-1).*(obj.omega.*timeseries(i) ...
                            - obj.wavenumber*(obj.marker.location(j,1).*cos(obj.direction*pi/180) ...
                            + obj.marker.location(j,2).*sin(obj.direction*pi/180)) + obj.phase)));
                        obj.waveAmpTimeViz(i,j+1) = rampFunction(i).*sum(tmp14,'all');
                    end
                end
            end
        end

        function waveElevUser(obj,rampTime,maxIt,data,time)
            % Calculate imported wave elevation time history
            % used by: :meth:`waveClass.setup`.
            rampFunction = (1+cos(pi+pi*time/rampTime))/2;
            rampFunction(time>=rampTime) = 1;

            obj.waveAmpTime = zeros(maxIt+1,2);
            data_t = data(:,1)';                    % Data Time [s]
            data_x = data(:,2)';                    % Wave Surface Elevation [m]
            obj.waveAmpTime(:,1) = time;
            obj.waveAmpTime(:,2) = rampFunction.*interp1(data_t,data_x,time);

            if any(~isempty(obj.marker.location))
                warning('Cannot use wave gauges or visualization markers with eta import. Gauges and markers removed.');
                obj.marker.location = [];
            end
        end

        function printWaveSpectrumType(obj)
            % Lists the wave spectrum type
            % Used by  :meth:`waveClass.listInfo`.
            if strcmp(obj.spectrumType,'JS')
                fprintf('\tSpectrum Type                        = JONSWAP \n')
            elseif strcmp(obj.spectrumType,'PM')
                fprintf('\tSpectrum Type                        = Pierson-Moskowitz  \n')
            elseif strcmp(obj.spectrumType,'spectrumImport')
                fprintf('\tSpectrum Type                        = Imported Spectrum \n')
            end
        end
    end
end

author: WEC-Sim

repo: WEC-Sim

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


text:
function cubePoissonP2
%% CUBEPOISSONP2 solves Poisson equation in a cube using quadratic element.
%
% Copyright (C) 2008 Long Chen. See COPYRIGHT.txt for details.

close all;
%% Parameters
maxIt = 3;
N = zeros(maxIt,1);
h = zeros(maxIt,1);

%% Generate an initial mesh
[node,elem] = cubemesh([-1,1,-1,1,-1,1],0.5);

%% Get the data of the pde
pde = sincosdata3;
% pde = polydata3;

%% Set up boundary condition
bdFlag = setboundary3(node,elem,'Dirichlet');
% bdFlag = setboundary3(node,elem,'Neumann');

%% Finite Element Method
errL2 = zeros(maxIt,1);
errH1 = zeros(maxIt,1);
erruIuh = zeros(maxIt,1);
for k=1:maxIt
    % refine grid
    [node,elem,bdFlag] = uniformrefine3(node,elem,bdFlag);
    % solve the equation
    [soln,eqn] = Poisson3P2(node,elem,bdFlag,pde);
    uh = soln.u;
    N(k) = length(uh);
    h(k) = 1./(size(node,1)^(1/3)-1);
    errH1(k) = getH1error3(node,elem,pde.Du,uh);
    errL2(k) = getL2error3(node,elem,pde.exactu,uh);
    uI = pde.exactu([node; (node(eqn.edge(:,1),:)+node(eqn.edge(:,2),:))/2]);
    erruIuh(k) = sqrt((uh-uI)'*eqn.Lap*(uh-uI));
    %error(k) = norm(u(1:N(k))-uI(1:N(k)));
end

%% Plot convergence rates
figure;
showrateh3(h,errH1,2,'-*', '|| Du - Du_h ||', ...
           h,errL2,2,'k-+', '|| u - u_h ||', ...
           h,erruIuh,2,'m-+','|| D u_I - D u_h ||');

fprintf('\n');
disp('Table: Error')
colname = {'#Dof','h','|| u-u_h ||','|| Du-Du_h ||','|| Du_I-Du_h ||'};
disptable(colname,N,[],h,'%0.3e',errL2,'%0.5e',errH1,'%0.5e',erruIuh,'%0.5e');

repo: ifem


text:
close all;
clear all;
clc;
rng('default');
% Create the directory for storing images
[status_code,message,message_id] = mkdir('bin');

% Signal space
N = 1000;
% Number of measurements
M = 200;
K = 50;

omp_success = 0;
gomp_success = 0;
for nt=1:100
    % Sensing matrix
    Phi = spx.dict.simple.gaussian_dict(M, N);
    % Construct the signal generator.
    gen  = spx.data.synthetic.SparseSignalGenerator(N, K);
    % Generate bi-uniform signals
    x = gen.gaussian;
    % Measurement vectors
    y = Phi.apply(x);

    % OMP solver instance
    solver = spx.pursuit.single.OrthogonalMatchingPursuit(Phi, K);
    % Solve the sparse recovery problem
    omp_result = solver.solve(y);
    % Solution vector
    z = omp_result.z;
    omp_stats = spx.commons.sparse.recovery_performance(Phi, K, y, x, z);
    %fprintf('OMP\n');
    %spx.commons.sparse.print_recovery_performance(omp_stats);

    % GOMP solver instance
    L = 2;
    options.ls_method = 'chol';
    z = spx.fast.gomp(double(Phi), y, K, L, 1e-6, options);
    gomp_stats = spx.commons.sparse.recovery_performance(Phi, K, y, x, z);
    %fprintf('GOMP\n');
    %spx.commons.sparse.print_recovery_performance(gomp_stats);

    fprintf('K=%d, Trial: %d, OMP: %s, GOMP: %s\n', ...
        K, nt, spx.io.true_false_short(omp_stats.success), ...
        spx.io.true_false_short(gomp_stats.success));
    omp_success = omp_success  + omp_stats.success;
    gomp_success = gomp_success  + gomp_stats.success;
end
fprintf('TOTAL: OMP : %d, GOMP: %d\n', omp_success, gomp_success);

author: indigits

repo: sparse-plex


text:
function varargout=fopenmkd(varargin)
%
% For calling details please see v_fopenmkd.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_fopenmkd.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: fopenmkd.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_fopenmkd(varargin{:});
else
    v_fopenmkd(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


%DAS_BF
clear;
close all;
load('Computed_RIRs.mat');

nChannels = size(m_pos,1);
dMic = m_pos(2, 2) - m_pos(1, 2);
c = 340;
delta = pdist2(s_pos(1,:), m_pos(2,:)) - pdist2(s_pos(1,:), m_pos(1,:));
DOA_est = acosd(delta/dMic);

%==================Generate array signals==================================%
speechfilename = {'6319-275224-0008.flac', '6319-275224-0011.flac'};
noisefilename = {'noise1.wav', 'noise2.wav'};

[source1, fs] = audioread(speechfilename{1});
[source2, ~] = audioread(noisefilename{1});
[noise, fs_n] = audioread(noisefilename{1});
noise = resample(noise,fs,fs_n);

n_f = fs * 10; %2 seoconds
source1 = source1(1:n_f);
source2 = source2(1:n_f) ./ 3;
noise = noise(1:n_f);

rir = RIR_sources(:,:,1);
speech1 = fftfilt(rir, source1).*30;

rir = RIR_sources(:,:,2);
speech2 = fftfilt(rir,source2).*30;

arraySignal = speech1 + repmat(noise, 1, 5);
%==========================================================================%

%=============================plot Beamformer==============================%
%h_mvdr = mvdr(d, Phi_u);
theta = linspace(0, 2*pi, 200); % incidence angle range
%plot
figure(1);
f = [10000];
hw = dsb(nChannels, dMic, DOA_est, f);        % delay-and-sum filter
p = plotBeamformer(nChannels, f, dMic, hw, theta);

% f = linspace(0, fs, 256);% freqeuncy vector
% hw = dsb(nChannels, dMic, DOA_est, f);        % delay-and-sum filter
% p = plotBeamformer(nChannels, f, dMic, hw, theta);
%==========================================================================%

%=============================apply Beamformer=============================%
x_das = applyBeamforming(arraySignal.', dMic, fs,  DOA_est, 'DSB', false);
audiowrite('output/x_dsb.wav', x_das, fs);

x_mvdr = applyBeamforming(arraySignal.', dMic, fs,  DOA_est, 'MVDR', false);
audiowrite('output/x_mvdr.wav', x_mvdr, fs);
%==========================================================================%

%====================Quality evaluation====================================%
%Compare
%[y,fs] = audioread('wav/speech1.wav');
%audiowrite('wav/speech1_16k.wav',y,16000);
%[y,fs] = audioread('output/x_dsb.wav');
%audiowrite('output/x_out_16k.wav',y,16000);
%%pesq = pesq('wav/speech1_16k.wav', 'wav/x_out_16k.wav')
%==========================================================================%