text:
function [x,fval,gp] = gplite_fmin(gp,x0,maxflag)
%GPLITE_FMIN Find global minimum (or maximum) of GP.

if nargin < 2; x0 = 0; end
if nargin < 3 || isempty(maxflag); maxflag = 0; end

MaxBnd = 10;
hpd_frac = 0.5;
D = size(gp.X,2);
N0 = size(x0,1);
Nstarts = max(3,N0);

diam = max(gp.X) - min(gp.X);
LB = min(gp.X) - MaxBnd*diam;
UB = max(gp.X) + MaxBnd*diam;

% First, train GP
if ~isfield(gp,'post') || isempty(gp.post)
    % How many samples for the GP?
    if isfield(gp,'Ns') && ~isempty(gp.Ns); Ns_gp = gp.Ns; else; Ns_gp = 0; end
    options.Nopts = 1;  % Do only one optimization
    gp = gplite_train(...
        [],Ns_gp,gp.X,gp.y,gp.covfun,gp.meanfun,gp.noisefun,[],[],options);
end

% Start from the min (or max) of the training data
if maxflag
    [~,ord] = sort(gp.y,'descend');
else
    [~,ord] = sort(gp.y,'ascend');
end

% Take best for sure
X = gp.X(ord,:);
x0 = [x0; X(1,:)];
X(1,:) = [];

if Nstarts > N0+1
    Nx = size(X,1);
    N_hpd = ceil(Nx*hpd_frac);
    idx = randperm(N_hpd,min(Nstarts-N0,N_hpd));
    x0 = [x0; X(idx,:)];
end

N0 = size(x0,1);
x = zeros(N0,D);
f = zeros(N0,1);
opts = optimoptions('fmincon','GradObj','off','Display','off');
for i = 1:N0
    [x(i,:),f(i)] = fmincon(@(x) optfun(x,gp,maxflag),x0(i,:),[],[],[],[],LB,UB,[],opts);
end

[fval,idx] = min(f);
x = x(idx,:);

if maxflag; fval = -fval; end

end

function [f,df] = optfun(x,gp,maxflag)

if nargout > 1
    [f,df] = gplite_pred(gp,x);
else
    f = gplite_pred(gp,x);
end

if maxflag  % Want to find maximum, swap sign
    f = -f;
    if nargout > 1; df = -df; end
end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
classdef LogisticRegression < handle
    %LOGISTICREGRESSION  Logistic Regression classifier
    %
    % ## Logistic Regression
    %
    % ML implements logistic regression, which is a probabilistic
    % classification technique. Logistic Regression is a binary classification
    % algorithm which is closely related to Support Vector Machines (SVM).
    % Like SVM, Logistic Regression can be extended to work on multi-class
    % classification problems like digit recognition (i.e. recognizing digitis
    % like 0,1 2, 3,... from the given images). This version of Logistic
    % Regression supports both binary and multi-class classifications (for
    % multi-class it creates a multiple 2-class classifiers). In order to
    % train the logistic regression classifier, Batch Gradient Descent and
    % Mini-Batch Gradient Descent algorithms are used (see [BatchDesWiki]).
    % Logistic Regression is a discriminative classifier (see [LogRegTomMitch]
    % for more details). Logistic  Regression is implemented as a C++ class in
    % cv.LogisticRegression.
    %
    % In Logistic Regression, we try to optimize the training parameter
    % `theta` such that the hypothesis `0 <= h_theta(x) <= 1` is achieved. We
    % have `h_theta(x) = g(h_theta(x))` and `g(z)=1/(1+e^(-z))` as the
    % logistic or sigmoid function. The term "Logistic" in Logistic Regression
    % refers to this function. For given data of a binary classification
    % problem of classes 0 and 1, one can determine that the given data
    % instance belongs to class 1 if `h_theta(x) >= 0.5` or class 0 if
    % `h_theta(x) < 0.5`.
    %
    % In Logistic Regression, choosing the right parameters is of utmost
    % importance for reducing the training error and ensuring high training
    % accuracy:
    %
    % * The learning rate can be set with `LearningRate` property. It
    %   determines how fast we approach the solution. It is a positive real
    %   number.
    % * Optimization algorithms like Batch Gradient Descent and Mini-Batch
    %   Gradient Descent are supported in cv.LogisticRegression. It is
    %   important that we mention the number of iterations these optimization
    %   algorithms have to run. The number of iterations can be set with
    %   `Iterations` property. This parameter can be thought as number of
    %   steps taken and learning rate specifies if it is a long step or a
    %   short step. This and previous parameter define how fast we arrive at a
    %   possible solution.
    % * In order to compensate for overfitting regularization is performed,
    %   which can be enabled with `Regularization` property. One can specify
    %   what kind of regularization has to be performed by passing one of
    %   regularization kinds to this property.
    % * Logistic regression implementation provides a choice of 2 training
    %   methods with Batch Gradient Descent or the MiniBatch Gradient Descent.
    %   To specify this, set `TrainMethod` property as either `Batch` or
    %   `MiniBatch`. If training method is set to `MiniBatch`, the size of the
    %   mini batch has to be to a positive integer set with `MiniBatchSize`
    %   property.
    %
    % ## Example
    % A sample set of training parameters for the Logistic Regression
    % classifier can be initialized as follows:
    %
    %     lr = cv.LogisticRegression();
    %     lr.LearningRate = 0.05;
    %     lr.Iterations = 1000;
    %     lr.Regularization = 'L2';
    %     lr.TrainMethod = 'MiniBatch';
    %     lr.MiniBatchSize = 10;
    %
    % ## References
    % [LogRegWiki]:
    % > [Logistic regression](https://en.wikipedia.org/wiki/Logistic_regression)
    %
    % [BatchDesWiki]:
    % > [Gradient descent optimization](https://en.wikipedia.org/wiki/Gradient_descent_optimization)
    %
    % [LogRegTomMitch]:
    % > "Generative and Discriminative Classifiers: Naive Bayes and Logistic
    % > Regression" in Machine Learning, Tom Mitchell.
    % > [Machine Learning](http://www.cs.cmu.edu/~tom/NewChapters.html)
    %
    % [RenMalik2003]:
    % > "Learning a Classification Model for Segmentation". Proc. CVPR, Nice,
    % > France (2003).
    %
    % See also: cv.LogisticRegression.LogisticRegression, fitglm
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % The learning rate of the optimization algorithm.
        %
        % The higher the value, faster the rate and vice versa. If the value
        % is too high, the learning algorithm may overshoot the optimal
        % parameters and result in lower training accuracy. If the value is
        % too low, the learning algorithm converges towards the optimal
        % parameters very slowly. The value must a be a positive real number.
        % You can experiment with different values with small increments as in
        % 0.0001, 0.0003, 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, ... and select
        % the learning rate with less training error. Default 0.001
        LearningRate
        % Number of iterations.
        %
        % The number of iterations required for the learing algorithm
        % (Gradient Descent or Mini Batch Gradient Descent). It has to be a
        % positive integer. You can try different number of iterations like in
        % 100, 1000, 2000, 3000, 5000, 10000, .. so on. Default 1000
        Iterations
        % Kind of regularization to be applied.
        %
        % Default 'L2'. Possible values:
        %
        % * __Disable__ Regularization disabled.
        % * __L1__ L1 norm.
        % * __L2__ L2 norm.
        Regularization
        % Kind of training method used to train the classifier.
        %
        % Default 'Batch'. Possible values:
        %
        % * __Batch__ batch gradient descent.
        % * __MiniBatch__ Mini-Batch Gradient Descent. Set `MiniBatchSize` to
        %   a positive integer when using this method.
        TrainMethod
        % Number of training samples taken in each step of Mini-Batch Gradient
        % Descent.
        %
        % Will only be used if using `MiniBatch` training algorithm. It has to
        % take values less than the total number of training samples.
        % Default 1
        MiniBatchSize
        % Termination criteria of the training algorithm.
        %
        % A struct with the following fields is accepted:
        %
        % * __type__ one of 'Count', 'EPS', 'Count+EPS'. default 'Count+EPS'
        % * __maxCount__ maximum number of iterations. default `Iterations`
        % * __epsilon__ tolerance value. default `LearningRate`
        TermCriteria
    end

    %% Constructor/destructor
    methods
        function this = LogisticRegression(varargin)
            %LOGISTICREGRESSION  Creates/trains a logistic regression model
            %
            %     model = cv.LogisticRegression()
            %     model = cv.LogisticRegression(...)
            %
            % The first variant creates Logistic Regression model with default
            % parameters.
            %
            % The second variant accepts the same parameters as the train
            % method, in which case it forwards the call after construction.
            %
            % See also: cv.LogisticRegression, cv.LogisticRegression.train
            %
            this.id = LogisticRegression_(0, 'new');
            if nargin > 0
                this.train(varargin{:});
            end
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     model.delete()
            %
            % See also: cv.LogisticRegression
            %
            if isempty(this.id), return; end
            LogisticRegression_(this.id, 'delete');
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     model.clear()
            %
            % The method clear does the same job as the destructor: it
            % deallocates all the memory occupied by the class members. But
            % the object itself is not destructed and can be reused further.
            % This method is called from the destructor, from the `train` and
            % `load` methods, or even explicitly by the user.
            %
            % See also: cv.LogisticRegression.empty, cv.LogisticRegression.load
            %
            LogisticRegression_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Returns true if the algorithm is empty
            %
            %     b = model.empty()
            %
            % ## Output
            % * __b__ Returns true if the algorithm is empty (e.g. in the very
            %   beginning or after unsuccessful read).
            %
            % See also: cv.LogisticRegression.clear, cv.LogisticRegression.load
            %
            b = LogisticRegression_(this.id, 'empty');
        end

        function varargout = save(this, filename)
            %SAVE  Saves the algorithm parameters to a file or a string
            %
            %     model.save(filename)
            %     str = model.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to. In case of string
            %   output, only the filename extension is used to determine the
            %   output format (XML or YAML).
            %
            % ## Output
            % * __str__ optional output. If requested, the model is persisted
            %   to a string in memory instead of writing to disk.
            %
            % This method stores the complete model state to the specified
            % XML or YAML file (or to a string in memory, based on the number
            % of output arguments).
            %
            % See also: cv.LogisticRegression.load
            %
            [varargout{1:nargout}] = LogisticRegression_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     model.load(filename)
            %     model.load(str, 'FromString',true)
            %     model.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __filename__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model (switches
            %   between `Algorithm<T>::load()` and
            %   `Algorithm<T>::loadFromString()` C++ methods). default false
            %
            % This method loads the complete model state from the specified
            % XML or YAML file (either from disk or serialized string). The
            % previous model state is cleared.
            %
            % See also: cv.LogisticRegression.save
            %
            LogisticRegression_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = model.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.LogisticRegression.save, cv.LogisticRegression.load
            %
            name = LogisticRegression_(this.id, 'getDefaultName');
        end
    end

    %% StatModel
    methods
        function count = getVarCount(this)
            %GETVARCOUNT  Returns the number of variables in training samples
            %
            %     count = model.getVarCount()
            %
            % ## Output
            % * __count__ number of variables in training samples (plus one to
            %   account for the implicitly prepended bias/intercept term).
            %
            % See also: cv.LogisticRegression.train
            %
            count = LogisticRegression_(this.id, 'getVarCount');
        end

        function b = isTrained(this)
            %ISTRAINED  Returns true if the model is trained
            %
            %     b = model.isTrained()
            %
            % ## Output
            % * __b__ Returns true if the model is trained, false otherwise.
            %
            % See also: cv.LogisticRegression.empty, cv.LogisticRegression.train
            %
            b = LogisticRegression_(this.id, 'isTrained');
        end

        function b = isClassifier(this)
            %ISCLASSIFIER  Returns true if the model is a classifier
            %
            %     b = model.isClassifier()
            %
            % ## Output
            % * __b__ Returns true if the model is a classifier, false if the
            %   model is a regressor.
            %
            % Always true for logistic regression.
            %
            % See also: cv.LogisticRegression.isTrained
            %
            b = LogisticRegression_(this.id, 'isClassifier');
        end

        function status = train(this, samples, responses, varargin)
            %TRAIN  Trains the statistical model
            %
            %     status = model.train(samples, responses)
            %     status = model.train(csvFilename, [])
            %     [...] = model.train(..., 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __samples__ matrix of training samples. It should have
            %   `single` type. By default, each row represents a sample (see
            %   the `Layout` option).
            % * __responses__ matrix of associated responses. A vector of
            %   categorical labels, stored in an array of type `single`.
            % * __csvFilename__ The input CSV file name from which to load
            %   dataset. In this variant, you should set the second argument
            %   to an empty array.
            %
            % ## Output
            % * __status__ Success flag.
            %
            % ## Options
            % * __Data__ Training data options, specified as a cell array of
            %   key/value pairs of the form `{'key',val, ...}`. See below.
            % * __Flags__ The optional training flags, model-dependent.
            %   Not used. default 0
            %
            % ### Options for `Data` (first variant with samples and reponses)
            % * __Layout__ Sample types. Default 'Row'. One of:
            %   * __Row__ each training sample is a row of samples.
            %   * __Col__ each training sample occupies a column of samples.
            % * __VarIdx__ vector specifying which variables to use for
            %   training. It can be an integer vector (`int32`) containing
            %   0-based variable indices or logical vector (`uint8` or
            %   `logical`) containing a mask of active variables. Not set by
            %   default, which uses all variables in the input data.
            % * __SampleIdx__ vector specifying which samples to use for
            %   training. It can be an integer vector (`int32`) containing
            %   0-based sample indices or logical vector (`uint8` or
            %   `logical`) containing a mask of training samples of interest.
            %   Not set by default, which uses all samples in the input data.
            % * __SampleWeights__ optional floating-point vector with weights
            %   for each sample. Some samples may be more important than
            %   others for training. You may want to raise the weight of
            %   certain classes to find the right balance between hit-rate and
            %   false-alarm rate, and so on. Not set by default, which
            %   effectively assigns an equal weight of 1 for all samples.
            % * __VarType__ optional vector of type `uint8` and size
            %   `<num_of_vars_in_samples> + <num_of_vars_in_responses>`,
            %   containing types of each input and output variable. By default
            %   considers all variables as numerical (both input and output
            %   variables). In case there is only one output variable of
            %   integer type, it is considered categorical. You can also
            %   specify a cell-array of strings (or as one string of single
            %   characters, e.g 'NNNC'). Possible values:
            %   * __Numerical__, __N__ same as 'Ordered'
            %   * __Ordered__, __O__ ordered variables
            %   * __Categorical__, __C__ categorical variables
            % * __MissingMask__ Indicator mask for missing observation (not
            %   currently implemented). Not set by default
            % * __TrainTestSplitCount__ divides the dataset into train/test
            %   sets, by specifying number of samples to use for the test set.
            %   By default all samples are used for the training set.
            % * __TrainTestSplitRatio__ divides the dataset into train/test
            %   sets, by specifying ratio of samples to use for the test set.
            %   By default all samples are used for the training set.
            % * __TrainTestSplitShuffle__ when splitting dataset into
            %   train/test sets, specify whether to shuffle the samples.
            %   Otherwise samples are assigned sequentially (first train then
            %   test). default true
            %
            % ### Options for `Data` (second variant for loading CSV file)
            % * __HeaderLineCount__ The number of lines in the beginning to
            %   skip; besides the header, the function also skips empty lines
            %   and lines staring with '#'. default 1
            % * __ResponseStartIdx__ Index of the first output variable. If
            %   -1, the function considers the last variable as the response.
            %   If the dataset only contains input variables and no responses,
            %   use `ResponseStartIdx = -2` and `ResponseEndIdx = 0`, then the
            %   output variables vector will just contain zeros. default -1
            % * __ResponseEndIdx__ Index of the last output variable + 1. If
            %   -1, then there is single response variable at
            %   `ResponseStartIdx`. default -1
            % * __VarTypeSpec__ The optional text string that specifies the
            %   variables' types. It has the format
            %   `ord[n1-n2,n3,n4-n5,...]cat[n6,n7-n8,...]`. That is, variables
            %   from `n1` to `n2` (inclusive range), `n3`, `n4` to `n5` ...
            %   are considered ordered and `n6`, `n7` to `n8` ... are
            %   considered as categorical. The range
            %   `[n1..n2] + [n3] + [n4..n5] + ... + [n6] + [n7..n8]` should
            %   cover all the variables. If `VarTypeSpec` is not specified,
            %   then algorithm uses the following rules:
            %   * all input variables are considered ordered by default. If
            %     some column contains has non- numerical values, e.g.
            %     'apple', 'pear', 'apple', 'apple', 'mango', the
            %     corresponding variable is considered categorical.
            %   * if there are several output variables, they are all
            %     considered as ordered. Errors are reported when
            %     non-numerical values are used.
            %   * if there is a single output variable, then if its values are
            %     non-numerical or are all integers, then it's considered
            %     categorical. Otherwise, it's considered ordered.
            % * __Delimiter__ The character used to separate values in each
            %   line. default ','
            % * __Missing__ The character used to specify missing
            %   measurements. It should not be a digit. Although it's a
            %   non-numerical value, it surely does not affect the decision of
            %   of whether the variable ordered or categorical. default '?'
            % * __TrainTestSplitCount__ same as above.
            % * __TrainTestSplitRatio__ same as above.
            % * __TrainTestSplitShuffle__ same as above.
            %
            % The training algorithm uses one-vs-rest scheme for multiclass
            % case (i.e a binary problem is fit for each label).
            %
            % See also: cv.LogisticRegression.predict, cv.LogisticRegression.calcError
            %
            status = LogisticRegression_(this.id, 'train', samples, responses, varargin{:});
        end

        function [err,resp] = calcError(this, samples, responses, varargin)
            %CALCERROR  Computes error on the training or test dataset
            %
            %     err = model.calcError(samples, responses)
            %     err = model.calcError(csvFilename, [])
            %     [err,resp] = model.calcError(...)
            %     [...] = model.calcError(..., 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __samples__ See the train method.
            % * __responses__ See the train method.
            % * __csvFilename__ See the train method.
            %
            % ## Output
            % * __err__ computed error.
            % * __resp__ the optional output responses.
            %
            % ## Options
            % * __Data__ See the train method.
            % * __TestError__ if true, the error is computed over the test
            %   subset of the data, otherwise it's computed over the training
            %   subset of the data. Please note that if you loaded a
            %   completely different dataset to evaluate an already trained
            %   classifier, you will probably want not to set the test subset
            %   at all with `TrainTestSplitRatio` and specify
            %   `TestError=false`, so that the error is computed for the whole
            %   new set. Yes, this sounds a bit confusing. default false
            %
            % The method uses the predict method to compute the error. For
            % regression models the error is computed as RMS, for classifiers
            % as a percent of missclassified samples (0%-100%).
            %
            % See also: cv.LogisticRegression.train, cv.LogisticRegression.predict
            %
            [err,resp] = LogisticRegression_(this.id, 'calcError', samples, responses, varargin{:});
        end

        function [results,f] = predict(this, samples, varargin)
            %PREDICT  Predicts responses for input samples
            %
            %     [results,f] = model.predict(samples)
            %     [...] = model.predict(..., 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __samples__ The input data for the prediction algorithm. `MxN`
            %   matrix, where each row contains variables (features) of one
            %   object being classified. Should have `single` data type.
            %
            % ## Output
            % * __results__ Predicted labels as a column vector of `int32`
            %   type. In case `RawOutput` was set, it returns a `single`
            %   matrix of size `size(samples,1)-by-size(thetas,1)` which
            %   contains the raw output of the sigmoid function(s).
            % * __f__ The same as the response of the first sample, i.e
            %   `results(1)`.
            %
            % ## Options
            % * __Flags__ The optional predict flags, model-dependent. For
            %   convenience, you can set the individual flag options below,
            %   instead of directly setting bits here. default 0
            % * __RawOutput__ makes the method return the raw results (the
            %   value of the sigmoid function), not the class label.
            %   default false
            %
            % See also: cv.LogisticRegression.train, cv.LogisticRegression.calcError
            %
            [results,f] = LogisticRegression_(this.id, 'predict', samples, varargin{:});
        end
    end

    %% LogisticRegression
    methods
        function thetas = getLearntThetas(this)
            %GETLEARNTTHETAS  Returns the trained parameters
            %
            %     thetas = model.getLearntThetas()
            %
            % ## Output
            % * __thetas__ It returns learnt parameters of the Logistic
            %   Regression as a matrix of type `single` arranged across rows.
            %   For a two-class classifcation problem, it returns a single row
            %   matrix.
            %
            % `thetas` is a matrix of size `nclasses-by-model.getVarCount()`
            % if `nclasses>2`, otherwise `1-by-model.getVarCount()` if
            % `nclasses=2`.
            %
            % See also: cv.LogisticRegression.train
            %
            thetas = LogisticRegression_(this.id, 'get_learnt_thetas');
        end
    end

    %% Getters/Setters
    methods
        function value = get.Iterations(this)
            value = LogisticRegression_(this.id, 'get', 'Iterations');
        end
        function set.Iterations(this, value)
            LogisticRegression_(this.id, 'set', 'Iterations', value);
        end

        function value = get.LearningRate(this)
            value = LogisticRegression_(this.id, 'get', 'LearningRate');
        end
        function set.LearningRate(this, value)
            LogisticRegression_(this.id, 'set', 'LearningRate', value);
        end

        function value = get.MiniBatchSize(this)
            value = LogisticRegression_(this.id, 'get', 'MiniBatchSize');
        end
        function set.MiniBatchSize(this, value)
            LogisticRegression_(this.id, 'set', 'MiniBatchSize', value);
        end

        function value = get.Regularization(this)
            value = LogisticRegression_(this.id, 'get', 'Regularization');
        end
        function set.Regularization(this, value)
            LogisticRegression_(this.id, 'set', 'Regularization', value);
        end

        function value = get.TermCriteria(this)
            value = LogisticRegression_(this.id, 'get', 'TermCriteria');
        end
        function set.TermCriteria(this, value)
            LogisticRegression_(this.id, 'set', 'TermCriteria', value);
        end

        function value = get.TrainMethod(this)
            value = LogisticRegression_(this.id, 'get', 'TrainMethod');
        end
        function set.TrainMethod(this, value)
            LogisticRegression_(this.id, 'set', 'TrainMethod', value);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function slicefr(x,y,z,v,xs,ys,zs)
% This script is a routine built around Matlab's SLICE function.
% It does the same as SLICE, just, it plots a black 'frame' along the edges
% of each slice and a dotted line where two slices intersect.
% This works only of slices parallel to the co-ordinate axes.
%
% U. Theune, Geophysics, U of A
%

slice(x,y,z,v,xs,ys,zs)

hold on
axis tight
xl=get(gca,'xlim');
yl=get(gca,'ylim');
zl=get(gca,'zlim');

% Plot the frames
for i=1:length(xs)
    plot3([1 1 1 1 1]*xs(i),[yl(1) yl(1) yl(2) yl(2) yl(1)],[zl(1) zl(2) zl(2) zl(1) zl(1)],'k','linewidth',0.5)
end
for i=1:length(ys)
    plot3([xl(1) xl(1) xl(2) xl(2) xl(1)],[1 1 1 1 1]*ys(i),[zl(1) zl(2) zl(2) zl(1) zl(1)],'k','linewidth',0.5)
end
for i=1:length(zs)
    plot3([xl(1) xl(2) xl(2) xl(1) xl(1)],[yl(1) yl(1) yl(2) yl(2) yl(1)],[1 1 1 1 1]*zs(i),'k','linewidth',0.5)
end
% Plot the intersection lines
% - vertical
for i=1:length(xs)
    for j=1:length(ys)
        plot3([xs(i) xs(i)],[ys(j) ys(j)],[zl(1) zl(2)],':k','linewidth',0.25)
    end
end
% - and horizontal
for i=1:length(zs)
    for j=1:length(ys)
        plot3([xl(1) xl(2)],[ys(j) ys(j)],[zs(i) zs(i)],':k','linewidth',0.25)
    end
end
for i=1:length(zs)
    for j=1:length(xs)
        plot3([xs(j) xs(j)],[yl(1) yl(2)],[zs(i) zs(i)],':k','linewidth',0.25)
    end
end
hold off
axis normal

author: Sable

repo: mcbench-benchmarks


text:
function g = gauss(x, m, P)

% GAUSS Gaussian distribution value
%   GAUSS(X, M, P) gives the probability density at point X, following a
%   Gaussian or Normal distribution with mean M and covariances matrix P.
%
%   For X vector, the evaluation happens at a single point denoted by X,
%   and the result is a positive scalar in the range [0 , 1].
%
%   For X matrix, the evaluation happens at each row of X, and the result
%   is a row-vector of scalars, one for each column of X.

%   Copyright 2013 Joan Sola

d = size(x,1);
n = size(x,2);
MD2 = zeros(1,n);

if (size(m,1) == d) && (size(P,1) == d) && (size(P,2) == d)

    z = x - repmat(m, 1, n);
    for i = 1:n
        MD2(1,i) = z(:,i)' * P^-1 * z(:,i);
    end
    a = sqrt((2*pi)^d * det(P));
    g = exp(-0.5 * MD2) / a;

else
    error ('Input sizes don''t match')
end

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function varargout = abs(varargin)
%ABS Absolute value of a DISKFUN.
%   ABS(F) returns the absolute value of a DISKFUN. This function does not work
%   if the function passes through or becomes numerically close to zero.

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

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

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% Detect if the expected matte is a highly-transparent one
% This function implements the energy selection method described in
% Yagiz Aksoy, Tunc Ozan Aydin, Marc Pollefeys, "Designing Effective
% Inter-Pixel Information Flow for Natural Image Matting", CVPR, 2017.
% This is a very simple histogram-based classifier.

function ht = detectHighlyTransparent(image, trimap)

    image = reshape(im2double(image), [size(image, 1) * size(image, 2), 3]);
    trimap = im2double(trimap(:,:,1));

    fg = trimap > 0.8;
    bg = trimap < 0.2;
    unk = ~(fg | bg);
    fg = fg & imdilate(unk, ones(20));
    bg = bg & imdilate(unk, ones(20));

    fgi = image(fg, :);
    bgi = image(bg, :);
    uni = image(unk, :);

    fgh = [imhist(fgi(:, 1), 10); imhist(fgi(:, 3), 10); imhist(fgi(:, 3), 10);] / sum(fg(:));
    bgh = [imhist(bgi(:, 1), 10); imhist(bgi(:, 3), 10); imhist(bgi(:, 3), 10);] / sum(bg(:));
    unh = [imhist(uni(:, 1), 10); imhist(uni(:, 3), 10); imhist(uni(:, 3), 10);] / sum(unk(:));

    weights = ([fgh bgh]' * [fgh bgh]) \ ([fgh bgh]' * unh);
    recError = [fgh bgh] * weights - unh;
    recError = sqrt(sum(recError(:) .* recError(:))) / size(recError(:), 1);

    ht = recError > 0.0099;

end

author: yaksoy

repo: AffinityBasedMattingToolbox


text:
% © Rahul Kala, IIIT Allahabad, Creative Commons Attribution-ShareAlike 4.0 International License.
% The use of this code, its parts and all the materials in the text; creation of derivatives and their publication; and sharing the code publically is permitted without permission.

% Please cite the work in all materials as: R. Kala (2014) Code for Robot Path Planning using Probabilistic Roadmap (Octave), Indian Institute of Information Technology Allahabad, Available at: http://rkala.in/codes.html

function h=historic(a,b)
h = sqrt(sum((a-b).^2));

author: xuedidi


text:
function [ap, cmc] = compute_AP(good_image, junk_image, index)

cmc = zeros(length(index), 1);
ngood = length(good_image);

%remove junk_image
rows_junk = ismember(index,junk_image);
index(rows_junk) = [];

%find good_index
rows_good = find(ismember(index,good_image));
cmc(rows_good(1):end) = 1; %after first equals to one

ap = 0;
for i = 1:ngood
    d_recall = 1/ngood;
    precision = i/rows_good(i);
    if (rows_good(i)~=1)
        old_precision = (i-1)/(rows_good(i)-1);
    else
        old_precision = 1;
    end
    ap = ap + d_recall*(old_precision+precision)/2;
end

end

author: layumi

repo: Image-Text-Embedding


text:
function asa136_test ( )

%*****************************************************************************80
%
%% ASA136_TEST tests ASA136.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA136_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the ASA136 library.\n' );

  asa136_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA136_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [y, extra, invhess] = rbfevfwd(net, x, t, x_test, invhess)
%RBFEVFWD Forward propagation with evidence for RBF
%
%	Description
%	Y = RBFEVFWD(NET, X, T, X_TEST) takes a network data structure  NET
%	together with the input X and target T training data and input test
%	data X_TEST. It returns the normal forward propagation through the
%	network Y together with a matrix EXTRA which consists of error bars
%	(variance) for a regression problem or moderated outputs for a
%	classification problem.
%
%	The optional argument (and return value)  INVHESS is the inverse of
%	the network Hessian computed on the training data inputs and targets.
%	Passing it in avoids recomputing it, which can be a significant
%	saving for large training sets.
%
%	See also
%	FEVBAYES
%

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

y = rbffwd(net, x_test);
% RBF outputs must be linear, so just pass them twice (second copy is
% not used
if nargin == 4
  [extra, invhess] = fevbayes(net, y, y, x, t, x_test);
else
  [extra, invhess] = fevbayes(net, y, y, x, t, x_test, invhess);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
 function ob = Gtomo3(sys_type, mask, nx, ny, nz, varargin)
%function ob = Gtomo3(sys_type, mask, nx, ny, nz, options)
%|
%| Construct Gtomo3 object, for 3D tomographic image reconstruction,
%| especially for SPECT and perhaps PET.  (For CT use Gcone; see below.)
%| This object does 3d forward and backprojection using mex files.
%| See Gtomo3_test.m for example usage.
%|
%| Basically, you create an object calling:
%|		A = Gtomo3(sys_type, mask)
%| and then you can use it thereafter by typing commands like
%|		y = A * x;
%| which will auto-magically do the multiplication by calling a mex file.
%|
%| in
%|	sys_type	string	see 3D ASPIRE User's Guide
%|	mask	[nx ny nz] or string
%|		mask can be a filename, in which case the file is read.
%|		mask can be empty (for some system types, including 3s),
%|		in which case it is determined from backprojecting a "1" vector.
%|		Else, mask must be [nx ny nz] binary array of object support.
%|
%| options
%|	'nthread' 1, 2...	1, unless you have a multi-processor system
%|	'chat' 0|1		verbosity.  default: 1.
%|	'checkmask' 0|1		verify that the mask is correct.  default: 1.
%|	'view2d' 0|1		for 2z and 2dsc system types, use 2d views
%|				instead of sinograms, for use with OS.
%|				default: 0
%|	'permute213' 0|1	convert from 123 to 213 output order
%|				(e.g., to overcome vu vs st order of 3l@)
%|				default: 0
%|	'class' ''		default 'fatrix2'
%|	'scale' [1]		scale factor.  default '' (none)
%|
%| out
%|	ob [nd np]	np = sum(mask(:)), so it is already "masked"
%|
%| For more help, type 'f3d_mex' and see Gtomo3_test.m
%|
%| For X-ray CT, using Gcone() is recommended instead of Gtomo3().
%| Nevertheless, the following alternate usage is permitted for CT here:
%| function ob = Gtomo3(cg, ig, options)
%|
%| Copyright 01-04-22, Jeff Fessler, University of Michigan

if nargin == 1 && streq(sys_type, 'test'), Gtomo3_test, return, end
if nargin < 1, ir_usage, end

% handle new usage (for CT)
% (cg, ig, [options])
if isa(sys_type, 'strum')
	cg = sys_type;
	ig = mask;
	options = {};
	if isvar('nx')
		options = {nx};
		if isvar('ny')
			options = {nx, ny};
			if isvar('nz')
				options = {nx, ny, nz};
			end
		end
	end
	options = {options{:}, varargin{:}};

	sys_type = aspire_pair(cg, ig, 'system', '3l');
	ob = Gtomo3(sys_type, ig.mask, ig.nx, ig.ny, ig.nz, options{:});
return
end

arg.sys_type = sys_type;
arg.nx = nx;
arg.ny = ny;
arg.nz = nz;
arg.nxyz = nx*ny*nz;
arg.mask = logical(mask);

% defaults
arg.nthread = 1;
arg.chat = true;
arg.checkmask = false;
arg.scale = [];
arg.view2d = false;
arg.permute213 = false;
arg.class = 'fatrix2';
% arg.class = 'fatrix2';
arg.nthread_max = 8; % > 8 threads may crash matlab in Gtomo3_test for '3s'
		% due to f3d_mex bugs?  increase this at your own risk!

if ischar(arg.mask)
	arg.mask = fld_read(arg.mask); % read mask file
end

% options
arg = vararg_pair(arg, varargin);

if arg.nthread > arg.nthread_max
	warn('%d threads downgraded to %d, proceed with caution!', ...
		arg.nthread, arg.nthread_max)
	arg.nthread = arg.nthread_max;
end

if arg.view2d && ~streq(arg.sys_type, '2', 1)
	fail 'view2d option only relevant to 2dsc and 2z system types'
end
if arg.permute213 && ~streq(arg.sys_type, '3l', 2)
	warn 'permute213 designed only for "3l" system model!?'
end

% if no mask, then backproject to find it!
if isempty(arg.mask)
	arg.optr = f3d_mex('init', arg.sys_type, uint8(ones(nx,ny,nz)), ...
		int32(arg.nx), int32(arg.ny), int32(arg.nz), ...
		int32(arg.nthread), int32(arg.chat))
	warn('creating mask by summing will be a bit slow')
	t = double(f3d_mex('dims', arg.optr)); % data dimensions, matlab wants doubles
	arg.mask = f3d_mex('back', ones(t, 'single'), ...
			int32(0), int32(1), arg.optr, int32(0)) > 0;
	f3d_mex('free:1', arg.optr)
end

% initialize mex file
arg.optr = f3d_mex('init', arg.sys_type, uint8(arg.mask), ...
	int32(arg.nx), int32(arg.ny), int32(arg.nz), ...
	int32(arg.nthread), int32(arg.chat));
if arg.chat
%	printm('arg.optr = %lX', arg.optr)
end

arg.odim = Gtomo3_nn(arg);
arg.nd = prod(arg.odim);
arg.np = sum(arg.mask(:));
dim = [arg.nd arg.np]; % trick: make it masked by default!

% trick: handle 3l@ uv vs st mismatch (if requested) by permuting.
if arg.permute213
%	if ~isempty(after), fail('permute213 and scale not done'), end
%	after = @(x, is_transpose, istart, nblock) permute(x, [2 1 3]);
%	after = {'cascade_after', after};
	warn 'still testing permute213'
end

abs_arg = {'abs', @(ob) ob}; % aij nonnegative, except perhaps SPECT rotate
if streq(arg.sys_type, '3s@', 3)
	abs_arg = {}; % todo: check for spline_filter 3, else ok
end

switch arg.class
case 'fatrix2'
	forw = @(arg, x) Gtomo3_block_forw_array(arg, x, 1, 1);
	back = @(arg, y) Gtomo3_block_back_array(arg, y, 1, 1);
	forw_block = @(arg, x, iblock, nblock) ...
		Gtomo3_block_forw_array(arg, x, iblock, nblock);
	back_block = @(arg, y, iblock, nblock) ...
		Gtomo3_block_back_array(arg, y, iblock, nblock);

	ob = fatrix2('mask', arg.mask, 'arg', arg, abs_arg{:}, ...
		'does_many', 1, ...
		'forw', forw, 'back', back, ...
		'free', @Gtomo3_free, ...
		'forw_block', forw_block, 'back_block', back_block);
	if ~isempty(arg.scale)
		ob = arg.scale * ob;
		arg.scale = [];
	end

case 'Fatrix' % build Fatrix object
	if ~isempty(arg.scale)
		after = {'cascade_after', arg.scale};
	else
		after = {};
	end

	ob = Fatrix(dim, arg, 'caller', 'Gtomo3', abs_arg{:}, ...
		'forw', @Gtomo3_forw_Fatrix, ...
		'back', @Gtomo3_back_Fatrix, ...
		'free', @Gtomo3_free, ...
		'mtimes_block', @Gtomo3_mtimes_block_Fatrix, ...
		after{:});

otherwise
	fail('class %s', arg.class)
end

if arg.chat
	f3d_mex('show:1', arg.optr, int32(arg.chat)); % display header info
end

% check mask
if arg.checkmask
	printm 'Warn checking mask, which may be unnecessary'
	tmp = sum(ob).' > 0;
	tmp = embed(tmp, arg.mask);
	if ~isequal(tmp, arg.mask)
		if exist('imfill') == 2
			tmpfill = imfill(tmp);
			if ~isequal(tmpfill, arg.mask)
				if im
					im plc 1 2
					im(1, arg.mask, 'input mask')
					im(2, tmp, 'computed mask')
				end
				warn('mask-sum mismatch not corrected by fill')
			else
				warn('mask-sum mismatch; corrected by filling')
			end
			printm('Warn: this geometry may have poor sampling')
		else
			fail 'mask inconsistent with sum'
		end
	end
	printm('mask checking completed')
end

% Gtomo3_nn(): data dimensions
function nn = Gtomo3_nn(arg)
nn = double(f3d_mex('dims', arg.optr)); % silly matlab wants size to be double
if arg.view2d
	nn = nn([1 3 2]); % swap projection angle and slice index
end
if arg.permute213
	nn = nn([2 1 3]); % trick: swap vu to st for 3l@
end

% Gtomo3_forw_Fatrix(): y = G * x
% full projection
function y = Gtomo3_forw_Fatrix(arg, x)
y = Gtomo3_mtimes_block_Fatrix(arg, 0, x, 1, 1);

% Gtomo3_back_Fatrix(): x = G' * y
% full backprojection
function x = Gtomo3_back_Fatrix(arg, y)
x = Gtomo3_mtimes_block_Fatrix(arg, 1, y, 1, 1);

% Gtomo3_mtimes_block_Fatrix()
function y = Gtomo3_mtimes_block_Fatrix(arg, is_transpose, x, istart, nblock)
if is_transpose
	y = Gtomo3_block_back_Fatrix(arg, x, istart, nblock);
else
	y = Gtomo3_block_forw_Fatrix(arg, x, istart, nblock);
end

% Gtomo3_block_forw_Fatrix()
function y = Gtomo3_block_forw_Fatrix(arg, x, istart, nblock)

[x ei] = embed_in(x, arg.mask, arg.np);
y = Gtomo3_block_forw_array(arg, x, istart, nblock);
y = ei.shape(y);

% Gtomo3_block_back_Fatrix()
function x = Gtomo3_block_back_Fatrix(arg, y, istart, nblock)

na = arg.odim(end); % last index must be the one over which we subsetize
ia = istart:nblock:na;
[y eo] = embed_out(y, [arg.odim(1:end-1) length(ia)]); % [(M) *L]
x = Gtomo3_block_back_array(arg, y, istart, nblock);
x = eo.shape(x, arg.mask, arg.np);

% Gtomo3_free()
function Gtomo3_free(arg)
if arg.chat
	printm('freeing Gtomo3 object static memory')
	warn('further use of this Gtomo3 object is invalid')
end
f3d_mex('free:1', arg.optr);

%
% fatrix2 versions "_array"
%

% Gtomo3_block_forw_array()
% x: [nx ny nz *L]
% y: [odim(1:end-1) "na" *L] where "na" is # of angles in this block
function y = Gtomo3_block_forw_array(arg, x, istart, nblock)

for i4=1:size(x,4)
	y(:,:,:,i4) = f3d_mex('proj', single(x(:,:,:,i4)), ...
		int32(istart-1), int32(nblock), arg.optr, int32(arg.chat));
end
if arg.view2d % trick: for 2z and 2dsc
	y = permute(y, [1 3 2 4]);
end
if arg.permute213 % trick: for 3l@
	y = permute(y, [2 1 3 4]);
end
na = arg.odim(end); % last index must be the one over which we subsetize
ia = istart:nblock:na;
y = reshapee(y, prod(arg.odim(1:end-1)), na, []);
LL = size(y,3);
y = y(:,ia,:);
y = reshape(y, [arg.odim(1:end-1) length(ia) LL]); % [(M) *L]

% Gtomo3_block_back_array()
function x = Gtomo3_block_back_array(arg, y, istart, nblock)

na = arg.odim(end); % last index must be the one over which we subsetize
ia = istart:nblock:na;

n1 = prod(arg.odim(1:end-1));
y = reshapee(y, n1, length(ia), []);
yy = size(y); yy(2) = na; yy = zeros(yy);
yy(:,ia,:) = y;
dim = num2cell(arg.odim);
y = reshapee(yy, dim{:}, []);

if arg.permute213 % trick: for 3l@
	y = permute(y, [2 1 3 4]);
end
if arg.view2d % trick for 2d
	y = permute(y, [1 3 2 4]);
end
for i4 = 1:size(y,4)
	x(:,:,:,i4) = f3d_mex('back', single(y(:,:,:,i4)), ...
		int32(istart-1), int32(nblock), arg.optr, int32(arg.chat));
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [traj, infStates] = tapas_ehgf(r, p, varargin)
% Calculates the trajectories of the agent's representations under the HGF
%
% This function can be called in two ways:
%
% (1) tapas_ehgf(r, p)
%
%     where r is the structure generated by tapas_fitModel and p is the parameter vector in native space;
%
% (2) tapas_ehgf(r, ptrans, 'trans')
%
%     where r is the structure generated by tapas_fitModel, ptrans is the parameter vector in
%     transformed space, and 'trans' is a flag indicating this.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2020 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (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/>.

% Transform paramaters back to their native space if needed
if ~isempty(varargin) && strcmp(varargin{1},'trans');
    p = tapas_ehgf_transp(r, p);
end

% Number of levels
l = length(p)/5;

if l ~= floor(l)
    error('tapas:hgf:UndetNumLevels', 'Cannot determine number of levels');
end

% Unpack parameters
mu_0 = p(1:l);
sa_0 = p(l+1:2*l);
rho  = p(2*l+1:3*l);
ka   = p(3*l+1:4*l-1);
om   = p(4*l:5*l-2);
th   = exp(p(5*l-1));
al   = 1/p(5*l);

% Add dummy "zeroth" trial
u = [0; r.u(:,1)];

% Number of trials (including prior)
n = length(u);

% Assume that if u has more than one column, the last contains t
try
    if r.c_prc.irregular_intervals
        if size(u,2) > 1
            t = [0; r.u(:,end)];
        else
            error('tapas:hgf:InputSingleColumn', 'Input matrix must contain more than one column if irregular_intervals is set to true.');
        end
    else
        t = ones(n,1);
    end
catch
    if size(u,2) > 1
        t = [0; r.u(:,end)];
    else
        t = ones(n,1);
    end
end

% Initialize updated quantities

% Representations
mu = NaN(n,l);
pi = NaN(n,l);

% Other quantities
muhat = NaN(n,l);
pihat = NaN(n,l);
v     = NaN(n,l);
w     = NaN(n,l-1);
da    = NaN(n,l);
dau   = NaN(n,1);

% Representation priors
% Note: first entries of the other quantities remain
% NaN because they are undefined and are thrown away
% at the end; their presence simply leads to consistent
% trial indices.
mu(1,:) = mu_0;
pi(1,:) = 1./sa_0;

% Representation update loop
% Pass through trials
for k = 2:1:n
    if not(ismember(k-1, r.ign))

        %%%%%%%%%%%%%%%%%%%%%%
        % Effect of input u(k)
        %%%%%%%%%%%%%%%%%%%%%%

        % 1st level
        % ~~~~~~~~~
        % Prediction
        muhat(k,1) = mu(k-1,1) +t(k) *rho(1);

        % Precision of prediction
        pihat(k,1) = 1/(1/pi(k-1,1) +t(k) *exp(ka(1) *mu(k-1,2) +om(1)));

        % Input prediction error
        dau(k) = u(k) -muhat(k,1);

        % Updates
        pi(k,1) = pihat(k,1) +1/al;
        mu(k,1) = muhat(k,1) +1/pihat(k,1) *1/(1/pihat(k,1) +al) *dau(k);

        % Volatility prediction error
        da(k,1) = (1/pi(k,1) +(mu(k,1) -muhat(k,1))^2) *pihat(k,1) -1;

        if l > 2
            % Pass through higher levels
            % ~~~~~~~~~~~~~~~~~~~~~~~~~~
            for j = 2:l-1
                % Prediction
                muhat(k,j) = mu(k-1,j) +t(k) *rho(j);

                % Precision of prediction
                pihat(k,j) = 1/(1/pi(k-1,j) +t(k) *exp(ka(j) *mu(k-1,j+1) +om(j)));

                % Weighting factor
                v(k,j-1) = t(k) *exp(ka(j-1) *mu(k-1,j) +om(j-1));
                w(k,j-1) = v(k,j-1) *pihat(k,j-1);

                % Mean update
                mu(k,j) = muhat(k,j) +1/2 *1/pihat(k,j) *ka(j-1) *w(k,j-1) *da(k,j-1);

                % Ingredients of precision update which depend on the mean
                % update
                vv = t(k) *exp(ka(j-1) *mu(k,j) +om(j-1));
                pimhat = 1/(1/pi(k-1,j-1) +vv);
                ww = vv *pimhat;
                rr = (vv -1/pi(k-1,j-1)) *pimhat;
                dd = (1/pi(k,j-1) +(mu(k,j-1) -muhat(k,j-1))^2) *pimhat -1;

                % Precision update
                pi(k,j) = pihat(k,j) +max(0, 1/2 *ka(j-1)^2 *ww*(ww +rr*dd));

                % Volatility prediction error
                da(k,j) = (1/pi(k,j) +(mu(k,j) -muhat(k,j))^2) *pihat(k,j) -1;
            end
        end

        % Last level
        % ~~~~~~~~~~
        % Prediction
        muhat(k,l) = mu(k-1,l) +t(k) *rho(l);

        % Precision of prediction
        pihat(k,l) = 1/(1/pi(k-1,l) +t(k) *th);

        % Weighting factor
        v(k,l)   = t(k) *th;
        v(k,l-1) = t(k) *exp(ka(l-1) *mu(k-1,l) +om(l-1));
        w(k,l-1) = v(k,l-1) *pihat(k,l-1);

        % Mean update
        mu(k,l) = muhat(k,l) +1/2 *1/pihat(k,l) *ka(l-1) *w(k,l-1) *da(k,l-1);

        % Ingredients of the precision update which depend on the mean
        % update
        vv = t(k) *exp(ka(l-1) *mu(k,l) +om(l-1));
        pimhat = 1/(1/pi(k-1,l-1) +vv);
        ww = vv *pimhat;
        rr = (vv -1/pi(k-1,l-1)) *pimhat;
        dd = (1/pi(k,l-1) +(mu(k,l-1) -muhat(k,l-1))^2) *pimhat -1;

        pi(k,l) = pihat(k,l) +max(0, 1/2 *ka(l-1)^2 *ww*(ww +rr*dd));

        % Volatility prediction error
        da(k,l) = (1/pi(k,l) +(mu(k,l) -muhat(k,l))^2) *pihat(k,l) -1;
    else

        mu(k,:) = mu(k-1,:);
        pi(k,:) = pi(k-1,:);

        muhat(k,:) = muhat(k-1,:);
        pihat(k,:) = pihat(k-1,:);

        v(k,:)  = v(k-1,:);
        w(k,:)  = w(k-1,:);
        da(k,:) = da(k-1,:);

    end
end

% Remove representation priors
mu(1,:)  = [];
pi(1,:)  = [];

% Remove other dummy initial values
muhat(1,:) = [];
pihat(1,:) = [];
v(1,:)     = [];
w(1,:)     = [];
da(1,:)    = [];
dau(1)     = [];

% Create result data structure
traj = struct;

traj.mu     = mu;
traj.sa     = 1./pi;

traj.muhat  = muhat;
traj.sahat  = 1./pihat;

traj.v      = v;
traj.w      = w;
traj.da     = da;
traj.dau    = dau;

% Updates with respect to prediction
traj.ud = mu -muhat;

% Psi (precision weights on prediction errors)
psi        = NaN(n-1,l);
psi(:,1)   = 1./(al*pi(:,1));
psi(:,2:l) = pihat(:,1:l-1)./pi(:,2:l);
traj.psi   = psi;

% Epsilons (precision-weighted prediction errors)
epsi        = NaN(n-1,l);
epsi(:,1)   = psi(:,1) .*dau;
epsi(:,2:l) = psi(:,2:l) .*da(:,1:l-1);
traj.epsi   = epsi;

% Full learning rate (full weights on prediction errors)
wt        = NaN(n-1,l);
wt(:,1)   = psi(:,1);
wt(:,2:l) = 1/2 *(v(:,1:l-1) *diag(ka(1:l-1))) .*psi(:,2:l);
traj.wt   = wt;

% Create matrices for use by the observation model
infStates = NaN(n-1,l,4);
infStates(:,:,1) = traj.muhat;
infStates(:,:,2) = traj.sahat;
infStates(:,:,3) = traj.mu;
infStates(:,:,4) = traj.sa;

end

author: translationalneuromodeling

repo: tapas


text:
function module = noise_module_diagonal(varargin)

options = struct('init', struct(), ...
                 'prior', struct(), ...
                 'separation', 'rows', ...
                 'update_tau', 1);
[options, errmsg] = argparse(options, varargin{:});
error(errmsg);

init = struct('tau', 1);
[init, errmsg] = argparse(init, options.init);
error(errmsg);

prior = struct('a_tau', 1e-5, ...
               'b_tau', 1e-5);
[prior, errmsg] = argparse(prior, options.prior);
error(errmsg);

tau = [];
a_tau = [];
b_tau = [];
weights = [];

module.initialize = @initialize;
module.update = @update;
module.get_struct = @get_struct;

  function [Tau] = initialize(M,N)
  switch options.separation
   case 'rows'
    tau = nan(M,1);
    a_tau = nan(M,1);;
    b_tau = nan(M,1);
   case 'columns'
    tau = nan(1,N);
    a_tau = nan(1,N);
    b_tau = nan(1,N);
   otherwise
    error('Unknown separation');
  end

  tau(:) = init.tau(:);
  Tau = bsxfun(@plus, tau, zeros(M,N));
  end

  function S = get_struct()
  % TODO: not ready yet..
  S.module = 'noise_module_isotropic';
  S.tau = tau;
  S.a_tau = a_tau;
  S.b_tau = b_tau;
  S.prior = prior;
  S.init = init;
  S.options = options;
  end

  function [Tau, LogTau, KL_Tau] = update(iter, E2, Obs)

  % OBS is the number of observations for each element of E2.

  % <(Y-WX)^2>
  switch options.separation
   case 'rows'
    e2 = sum(E2,2);
    obs = sum(Obs,2);
   case 'columns'
    e2 = sum(E2,1);
    obs = sum(Obs,1);
   otherwise
    error('Unknown separation');
  end

  %
  % Update tau
  %

  if index_selected(iter, options.update_tau)

    a_tau(:) = prior.a_tau(:) + 0.5*obs(:);
    b_tau(:) = prior.b_tau(:) + 0.5*e2(:);

    tau = a_tau ./ b_tau;
    logtau = psi(a_tau) - log(b_tau);

    %
    % Compute Kullback-Leibler divergence KL(q(tau)||p(tau))
    %

    % <log p(Tau)>
    logpdf_p = sum(gamma_logpdf(tau(:), prior.a_tau(:), prior.b_tau(:), logtau(:)));

    % <log q(Tau)>
    logpdf_q = -sum(gamma_entropy(a_tau,b_tau));

    % KL(q||p)
    KL_Tau = logpdf_q - logpdf_p;

  else
    logtau = log(tau);
    KL_Tau = 0;
  end

  % Compute expectations of weighted Tau
  O = zeros(size(E2));
  Tau = bsxfun(@plus, tau, O);
  LogTau = bsxfun(@plus, logtau, O);

  end

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [vectorResult] = niftiCreateComplementString(vectorString)
%
%Take the complement of the vectorString input. So, if we receive AIR, we
%want to return ASR. That is, take the complement of each input with
%respect to the "correct" RPS format.

vectorResult = vectorString;

searchStringMap = containers.Map;

searchStringMap('L') = 'R';
searchStringMap('A') = 'P';
searchStringMap('I') = 'S';

for i = 1:numel(vectorResult)
    if searchStringMap.isKey(vectorResult(i))
        %We need to flip it!
        vectorResult(i) = searchStringMap(vectorResult(i));
    end %if
end %for

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [arrival,assoc,origin,W] = db_get_arrival_info(arid,dbname,varargin)
% [ARRIVAL, ASSOC, ORIGIN, W] = DB_GET_ARRIVAL_INFO(ARID_LIST, DATABASE_NAME) Returns pertinent
% fields from the arrival, assoc and origin tables based on input arrival id numbers (arid);
% ARRIVAL_LIST must be a vector of arrival id numbers (arid). The returned structures ARRIVAL contains
% fields pulled from each of the corresponding tables.
%
% Unlike db_get_origin_info, there is currently no method for entering arrivals by time.
% arid's must be used.

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

if isnumeric(arid)
	arrival.arid = reshape( arid , numel(arid), 1 );
else
    error('First input must be either arrival ids (arid)');
end

% Check for waveform
if numel(varargin) == 1
	W = varargin{1};
	if ~isa(W,'waveform')
		error('Third argument must be a WAVEFORM object');
	end
elseif numel(varargin) > 1
	error('Too many arguments');
end

% GET ORIGIN TIME
db = dbopen(dbname,'r');
db = dblookup_table(db,'arrival');
nlist = [];
for n = 1:length(arrival.arid);
	db1 = dbsubset(db,['arid==' num2str(arrival.arid(n))]);
	db2 = dblookup_table(db,'assoc');
	db1 = dbjoin(db1,db2);
	db2 = dblookup_table(db,'origin');
	db1 = dbjoin(db1,db2);
	[ origin.lat(n),origin.lon(n),origin.depth(n),origin.time_epoch(n),origin.ml(n),origin.orid(n)] = dbgetv(db1,'origin.lat','origin.lon','origin.depth','origin.time','origin.ml','origin.orid');
	[ arrival.sta{n},arrival.chan{n}] = dbgetv(db1,'arrival.sta','arrival.chan');
	[ arrival.time(n),arrival.iphase(n),arrival.deltim(n)] = dbgetv(db1,'arrival.time','arrival.iphase','arrival.deltim');
    [ assoc.delta(n),assoc.seaz(n),assoc.esaz(n),assoc.timeres(n) ] = dbgetv(db1,'assoc.delta','assoc.seaz','assoc.esaz','assoc.timeres');
end
dbclose(db);
origin.lat = reshape( origin.lat , numel(origin.lat), 1 );
origin.lon = reshape( origin.lon , numel(origin.lon), 1 );
origin.depth = reshape( origin.depth , numel(origin.depth), 1 );
origin.time_epoch = reshape( origin.time_epoch , numel(origin.time_epoch), 1 );
origin.time_matlab = datenum(strtime(origin.time_epoch));
origin.ml = reshape( origin.ml , numel(origin.ml), 1 );
origin.orid = reshape( origin.orid , numel(origin.orid), 1 );
%
arrival.time_epoch = reshape( arrival.time , numel(arrival.time), 1 );
arrival.time_matlab =  datenum(strtime(arrival.time_epoch));
arrival.iphase = reshape( arrival.iphase , numel(arrival.iphase), 1 );
arrival.deltim = reshape( arrival.deltim , numel(arrival.deltim), 1 );
arrival.sta = reshape( arrival.sta , numel(arrival.sta), 1 );
arrival.chan = reshape( arrival.chan , numel(arrival.chan), 1 );

%
assoc.delta = reshape( assoc.delta , numel(assoc.delta), 1 );
assoc.seaz = reshape( assoc.seaz , numel(assoc.seaz), 1 );
assoc.esaz = reshape( assoc.esaz , numel(assoc.esaz), 1 );
assoc.timeres = reshape( assoc.timeres , numel(assoc.timeres), 1 );
assoc.traveltime = arrival.time_matlab - origin.time_matlab;

if exist('W')
	for n = 1:numel(W)
		disp([ 'Retreiving origin, arrival and assoc info for orid: ' num2str(origin.orid(n)) ' ... '])
		W(n) = addfield(W(n),'ORIGIN_LAT',origin.lat(n));
		W(n) = addfield(W(n),'ORIGIN_LON',origin.lon(n));
		W(n) = addfield(W(n),'ORIGIN_DEPTH',origin.depth(n));
		W(n) = addfield(W(n),'ORIGIN_TIME_EPOCH',origin.time_epoch(n));
		W(n) = addfield(W(n),'ORIGIN_ML',origin.ml(n));
		W(n) = addfield(W(n),'ORIGIN_TIME_MATLAB',origin.time_matlab(n));
		W(n) = addfield(W(n),'ORIGIN_ORID',origin.orid(n));
		%
		W(n) = addfield(W(n),'ARRIVAL_TIME_EPOCH',arrival.time_epoch(n));
		W(n) = addfield(W(n),'ARRIVAL_TIME_MATLAB',arrival.time_matlab(n));
		W(n) = addfield(W(n),'ARRIVAL_IPHASE',arrival.iphase(n));
		W(n) = addfield(W(n),'ARRIVAL_DELTIME',arrival.deltime(n));
		%
		W(n) = addfield(W(n),'ASSOC_DELTA',assoc.delta(n));
		W(n) = addfield(W(n),'ASSOC_SEAZ',assoc.seaz(n));
		W(n) = addfield(W(n),'ASSOC_ESAZ',assoc.esaz(n));
		W(n) = addfield(W(n),'ASSOC_TIMERES',assoc.timeres(n));
	end
end

author: geoscience-community-codes

repo: GISMO


text:
function [g1, g2] = simwhiteXsimwhiteKernGradient(simKern1, simKern2, t1, varargin)

% SIMWHITEXSIMWHITEKERNGRADIENT Compute a cross gradient between two
% SIM-WHITE kernels.
% FORMAT
% DESC computes cross gradient of parameters of a cross kernel
% between two SIM-WHITE kernels for the multiple output kernel.
% ARG simKern1 : the kernel structure associated with the first SIM-WHITE
% kernel.
% ARG simKern2 : the kernel structure associated with the second SIM-WHITE
% kernel.
% ARG t1 : inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of the parameters of the first kernel, for
% ordering see simwhiteKernExtractParam.
% RETURN g2 : gradient of the parameters of the second kernel, for
% ordering see simwhiteKernExtractParam.
%
% FORMAT
% DESC computes cross kernel terms between two SIM-WHITE kernels for
% the multiple output kernel.
% ARG simKern1 : the kernel structure associated with the first SIM-WHITE
% kernel.
% ARG simKern2 : the kernel structure associated with the second SIM-WHITE
% kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of the parameters of the first kernel, for
% ordering see simwhiteKernExtractParam.
% RETURN g2 : gradient of the parameters of the second kernel, for
% ordering see simwhiteKernExtractParam.
%
% SEEALSO : multiKernParamInit, multiKernCompute, simwhiteKernParamInit,
% simwhiteKernExtractParam
%
% COPYRIGHT : David Luengo, 2009

% KERN

if nargin < 5
    t2 = t1;
else
    t2 = varargin{1};
end
covGrad = varargin{end};

if size(t1, 2) > 1 | size(t2, 2) > 1
  error('Input can only have one column');
end
if simKern1.variance ~= simKern2.variance
  error('Kernels cannot be cross combined if they have different variances.')
end

% Parameters of the kernels required in the computation
variance = simKern1.variance;
sensitivity1 = simKern1.sensitivity;
sensitivity2 = simKern2.sensitivity;
decay1 = simKern1.decay;
decay2 = simKern2.decay;

isStationary = (simKern1.isStationary == true) & (simKern2.isStationary == true);

% Auxiliary constants, vectors and matrices
g1 = zeros(1,3);
g2 = zeros(1,3);

T1 = repmat(t1, 1, size(t2, 1));
T2 = repmat(t2.', size(t1, 1), 1);
ind1 = (T1 < T2);
ind2 = ~ind1;
Dv = decay2 .* ind1 + decay1 .* ind2;
K1 = exp(-Dv.*abs(T1-T2));
if (isStationary == false)
    K2 = exp(-(decay1 * T1 + decay2 * T2));
end

if (isStationary == false)
    % Gradient w.r.t. the decays (D_p and D_q)
    g1(1) = (variance * sensitivity1 * sensitivity2 / (decay1 + decay2)) ...
        * sum(sum((-(K1 - K2)/(decay1 + decay2) ...
            - abs(T1 - T2) .* K1 .* ind2 + T1 .* K2) .* covGrad));
    g2(1) = (variance * sensitivity1 * sensitivity2 / (decay1 + decay2)) ...
        * sum(sum((-(K1 - K2)/(decay1 + decay2) ...
            - abs(T1 - T2) .* K1 .* ind1 + T2 .* K2) .* covGrad));
    % Gradient w.r.t. the variance of the driving process
    g1(2) = (sensitivity1 * sensitivity2 / (decay1 + decay2)) ...
        * sum(sum((K1-K2) .* covGrad));
    g2(2) = 0; % Otherwise it is counted twice
    % Gradient w.r.t. the sensitivities
    g1(3) = (variance * sensitivity2 / (decay1 + decay2)) ...
        * sum(sum((K1-K2) .* covGrad));
    g2(3) = (variance * sensitivity1 / (decay1 + decay2)) ...
        * sum(sum((K1-K2) .* covGrad));
else
    % Gradient w.r.t. the decays (D_p and D_q)
    g1(1) = - (variance * sensitivity1 * sensitivity2 / (decay1 + decay2)) ...
        * sum(sum((K1/(decay1 + decay2) + abs(T1 - T2) .* K1 .* ind2) .* covGrad));
    g2(1) = - (variance * sensitivity1 * sensitivity2 / (decay1 + decay2)) ...
        * sum(sum((K1/(decay1 + decay2) + abs(T1 - T2) .* K1 .* ind1) .* covGrad));
    % Gradient w.r.t. the variance of the driving process
    g1(2) = (sensitivity1 * sensitivity2 / (decay1 + decay2)) ...
        * sum(sum(K1 .* covGrad));
    g2(2) = 0; % Otherwise it is counted twice
    % Gradient w.r.t. the sensitivities
    g1(3) = (variance * sensitivity2 / (decay1 + decay2)) ...
        * sum(sum(K1 .* covGrad));
    g2(3) = (variance * sensitivity1 / (decay1 + decay2)) ...
        * sum(sum(K1 .* covGrad));
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ y, m ] = month_carry_eg_civil ( y, m )

%*****************************************************************************80
%
%% MONTH_CARRY_EG_CIVIL carries a year of months on the Egyptian Civil calendar.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, the year and month.
%
  months = year_length_months_eg_civil ( y );

  while ( 1 )

    if ( m <= months )
      break
    end

    m = m - months;
    y = y + 1;

  end

  return
end


text:
%  Y = typecast(X, type) converts a numeric value in X to the data type specified by type.
%*************************************************************************************
%
%  MATLAB (R) is a trademark of The Mathworks (R) Corporation
%
%  Function:    typecast
%  Filename:    typecast.m
%  Programmer:  James Tursa
%  Version:     3.00
%  Date:        March 17, 2011
%  Copyright:   (c) 2009, 2011 by James Tursa, All Rights Reserved
%
%  This code uses the BSD License:
%
%  Redistribution and use in source and binary forms, with or without
%  modification, are permitted provided that the following conditions are
%  met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
%  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%  POSSIBILITY OF SUCH DAMAGE.
%
% typecast is a mex function intended to mimic the MATLAB intrinsic typecast function
% for those users with older versions of MATLAB that do not have this intrinsic. Users
% of newer versions of MATLAB may be interested in this C-mex version to take advantage
% of the several extensions offered. This C-mex version of typecast differs from the
% intrinsic typecast in the following important aspects:
%
%                      Intrinsic typecast    C-mex typecast
%                      ------------------    --------------
% Type of copy:        Deep Data Copy        Shared Data Copy
% Allows complex:      No                    Yes
% Allows logical:      No                    Yes (cannot convert from complex)
% Allows char:         No                    Yes (cannot convert from complex)
% Allows non-vector:   No                    Yes
%
% Since this C-mex typecast produces a shared data copy of the original, it
% is more efficient than the MATLAB intrinsic typecast, which may be important
% if you are working with large variables. For non-vector inputs, the first
% non-singleton dimension must be compatible for the conversion.
%
% Building:
%
% typecast requires that a mex routine be built (one time only). This
% process is typically self-building the first time you call the function
% as long as you have the files typecast.m and typecast.c in the same
% directory somewhere on the MATLAB path. If you need to manually build
% the mex function, here are the commands:
%
% >> mex -setup
%   (then follow instructions to select a C or C++ compiler of your choice)
% >> mex typecast.c
%
% The usage is as follows (from The Mathworks website documentation):
%
% Syntax
%
% Y = typecast(X, type)
%
% Description
%
% Y = typecast(X, type) converts a value in X to the data type specified by type.
% Input X must be a full numeric or char or logical variable. The type input is a string
% set to one of the following: 'uint8', 'int8', 'uint16', 'int16', 'uint32', 'int32',
% 'uint64', 'int64', 'single', 'double', 'char' or 'logical'. typecast is different
% from the MATLAB cast function in that it does not alter the input data. typecast
% always returns the same number of bytes in the output Y as were in the input X.
% For example, casting the 16-bit integer 1000 to uint8 with typecast returns the full
% 16 bits in two 8-bit segments (3 and 232) thus keeping its original value
% (3*256 + 232 = 1000). The cast function, on the other hand, truncates the input value
% to 255.
%
% The output of typecast can be formatted differently depending on what system you use it on.
% Some computer systems store data starting with its most significant byte (an ordering
% called big-endian), while others start with the least significant byte (called little-endian).
%
% typecast issues an error if X contains fewer values than are needed to make an output value.
%
%*************************************************************************************

function varargout = typecast(varargin)
fname = 'typecast';
disp(' ');
disp(['Detected that the mex routine for ' fname ' is not yet built.']);
disp('Attempting to do so now ...');
disp(' ');
try
    mname = mfilename('fullpath');
catch
    mname = fname;
end
cname = [mname '.c'];
if( isempty(dir(cname)) )
    disp(['Cannot find the file ' fname '.c in the same directory as the']);
    disp(['file ' fname '.m. Please ensure that they are in the same']);
    disp('directory and try again. The following file was not found:');
    disp(' ');
    disp(cname);
    disp(' ');
    error(['Unable to compile ' fname '.c']);
else
    disp(['Found file ' fname '.c in ' cname]);
    disp(' ');
    disp('Now attempting to compile ...');
    disp('(If prompted, please press the Enter key and then select any C/C++');
    disp('compiler that is available, such as lcc.)');
    disp(' ');
    disp(['mex(''' cname ''',''-output'',''',mname,''')']);
    disp(' ');
    try
        mex(cname,'-output',mname);
        disp([ fname ' mex build completed ... you may now use ' fname '.']);
        disp(' ');
    catch
        disp(' ');
        error(['Unable to compile ' fname ' ... Contact author.']);
    end
    [varargout{1:nargout}] = typecast(varargin{:});
end
end


text:
function [MI, sequences, indexes, mini, sds, mn, sd] = motionCompDetectMotionMI(view,ROI,scans,plotMI,histogramClip)

%    [MI, sequences, indexes, mini, sds, mn, sd] = motionCompDetectMotionMI(view, [ROI], [scans], [plotMI], [histogramClip])
%
% gb 01/31/05
%
% Remove outliers within a scan and returns an array 'sequences' whose rows
% are the sequences between the outliers. The error computed is based on
% mutual information.
%
% Input arguments :
%       - view : current inplane
%       - ROI : Region of interest where the error has to be computed
%           Default = ones
%       - scan : scan to analyse.
%           Default = current scan
%       - plotMI : indicates if the user wants a plot of the Mutual Information.
%       Should be the value '1' if so.
%           Default = 1
%       - histogramClip : sets the clipping value of the histogram
%           Default = 0.1
%
% It returns :
%       - MI = Mutual Information between consecutive frames
%       - sequences = array representing the resulting sequences between the outliers
%       - indexes = indexes corresponding to the outliers
%       - mini = minimum of MI
%       - sds = array representing the standard deviations of each sequence
%       - mn = mean of MSE after having been clipped
%       - sd = standard deviation of MSE after having been clipped

global dataTYPES

% Initializes arguments and variables
if ieNotDefined('histogramClip')
    histogramClip = .1;
end

if ieNotDefined('scans')
    scans = selectScans(view);
end

if ieNotDefined('ROI')
    ROI = '';
end

curDataType = viewGet(view,'currentdatatype');
MI = [];
lastFrame = 0;
h = mrvWaitbar(0,'Computing Mutual Information...');

for scan = scans

    mrvWaitbar((scan - 1)/scans(end),h,['Computing Mutual Information for scan ' num2str(scan)]);

    clear tSeriesAllSlices;

	tSeriesAllSlices = motionCompLoadImages(view,scan);
	nFrames = size(tSeriesAllSlices,1);
    nSlices = size(tSeriesAllSlices,4);

    % Computes the Mutual Information when changing scans
    firstFrame = tSeriesAllSlices(1,:,:,:);
    if isequal(lastFrame,0)
        MI = [MI;0];
    else
        MI = [MI;motionCompMI(firstFrame,ROI,lastFrame)];
    end
    lastFrame = tSeriesAllSlices(nFrames,:,:,:);

    % Computes the Mutual Information
	MI = [MI;motionCompMI(tSeriesAllSlices,ROI)];

end
mrvWaitbar(1,h)

MI(1) = MI(2);

% Clips the histogram to calculate the mean and the standard deviation
% of MI without taking into account the outliers
histClip = mrAnatHistogramClip(MI,histogramClip,1,0);

% Computes the mean and the standard deviation of the clipped MI
mn = mean(histClip);
sd =  std(histClip);

% Finds the outliers. Are considered outliers all the values above the mean
% - two standard deviations
indexes = (find(MI < mn - 2*sd));
mini = min(MI);

% Finds the resulting sequences between the outliers.
sequences = horzcat([1;indexes + 2],[indexes - 2;nFrames - 1]);
sequences = sequences(find(sequences(:,2) - sequences(:,1) > 5),:);

% Computes all the standard deviations
sds = zeros(size(sequences,1),1);
for i = 1:size(sequences,1)
    sds(i) = std(MI(sequences(i,1):sequences(i,2)));
end

% Plots the curves
if ieNotDefined('plotMI')
    plotMI = 0;
end

if plotMI > 0
    figure
    title(['MI for ' view.sessionCode]);

    plot(MI)
    hold on
    plot(1:length(MI),mn*ones(1,length(MI)),'-.')
    plot([1 length(MI)],[mn - 2*sd, mn - 2*sd],'r')

    % Another measure can be used. Sometimes we can assume that the signal
    % is equally distributed around its mean. The explicit measure is :
    % mn - (max(MI) - mn)
    %
    % plot([1 nFrames],[2*mn - max(MI),2*mn - max(MI)],'black')

    for scanIndex = 1:(length(scans) - 1)
        plot([nFrames*scanIndex nFrames*scanIndex],[mini max(MI)],'-.','Color','m');
    end
end

close(h);


text:
%PLOT_POINT	Draw a feature point
%
% PLOT_POINT(P, OPTIONS) adds point markers to the current plot, where P (2xN)
% and each column is the point coordinate.
%
% Options::
%  'textcolor', colspec     Specify color of text
%  'textsize', size         Specify size of text
%  'bold'                   Text in bold font.
%  'printf', {fmt, data}    Label points according to printf format
%                           string and corresponding element of data
%  'sequence'               Label points sequentially
%
% Additional options are passed through to PLOT for creating the marker.
%
% Examples::
%   Simple point plot
%        P = rand(2,4);
%        plot_point(P);
%
%   Plot points with markers
%        plot_point(P, '*');
%
%   Plot points with square markers and labels 1 to 4
%        plot_point(P, 'sequence', 's');
%
%   Plot points with circles and annotations P1 to P4
%        data = [1 2 4 8];
%        plot_point(P, 'printf', {' P%d', data}, 'o');
%
% Notes::
% - The point(s) is added to the current plot.
% - 2D points only.
%
% See also PLOT, TEXT.

% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function plot_point(p, varargin)

    opt.textcolor = 'g';
    opt.textsize = 12;
    opt.printf = [];
    opt.sequence = false;
    opt.bold = false;

    [opt,arglist] = tb_optparse(opt, varargin);

    % default marker style
    if isempty(arglist)
        arglist = {'sb'};
    end

    % add stuff to pull .u and .v out of a vector of objects
    if ~isnumeric(p) && any(strcmp('u_', properties(p)))
        % p is an object with u_ and v_ properties
        p = [[p.u_]; [p.v_]];
    end

    holdon = ishold();
	hold on
	for i=1:numcols(p)
		plot(p(1,i), p(2,i), arglist{:});
        if opt.sequence
            show(p(:,i), '%d', i, opt);
        end

        if ~isempty(opt.printf)
            show(p(:,i), opt.printf{1}, opt.printf{2}(i), opt);
        end

	end
    if ~holdon
        hold off
    end
    figure(gcf)
end

function show(p, fmt, val, opt)
    if opt.bold
        fw = 'bold';
    else
        fw = 'normal';
    end
    text(p(1), p(2), sprintf([' ' fmt], val), ...
        'HorizontalAlignment', 'left', ...
        'VerticalAlignment', 'middle', ...
        'FontUnits', 'pixels', ...
        'FontSize', opt.textsize, ...
        'FontWeight', fw, ...
        'Color', opt.textcolor);
end

author: Allopart

repo: rbpf-gmapping


text:
function y = sdpfun(varargin)
%SDPFUN Obsolete - changed to blackbox

y = blackbox(varargin{:});

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% RMSEANALYSIS  RMSE analysis for slamtb.
%   RMSEANALYSIS is a script for evaluating average RMSE performance of
%   slamtb. It uses a slave version of SLAMTB, called SLAMTBSLAVE, and runs
%   it for N times with different kinds of landmarks.
%
%   Specify the types of landmarks in variable:  lmkTypes.
%   Specify the number of runs for each lmk in:  numRuns.
%   Specify the length of each run in         :  numFrames.
%   Specify the destination of log files in   :  logsDir.
%   Specify the rendering period in           :  rendPeriod.
%   Specify lmk parametrization with the flag :  reparametrize.
%   Specify N random seeds in                 :  randSeeds.
%
%   The result of this process is a set of log files. The contents of these
%   files can be plotted with NEESPLOTS.
%
%   See also SLAMTBSLAVE, NEESPLOTS, NEES, ROBOTNEES.

errorAnalysisFunction = @robotRmse;
logsDir    = '~/SLAM/logs/pose6d/RMSE/test/';
logFileHeader = 'frame, err: x y z r p y, std: x y z r p y';

% lmkTypes   = {'ahmPnt'};
% lmkTypes   = {'hmgPnt','idpPnt','ahmPnt'};
% lmkTypes   = {'hmgPnt'};
lmkTypes   = {'idpPnt'};
% lmkTypes   = {'idpPnt','ahmPnt'};
% lmkTypes   = {'plkLin','aplLin','hmgLin','ahmLin','idpLin'};
% lmkTypes   = {'aplLin','idpLin','ahmLin'};
% lmkTypes   = {'idpLin','ahmLin'};
% lmkTypes   = {'aplLin'};
% lmkTypes   = {'plkLin','aplLin'};

numRuns    = 25;
numFrames  = 400;

rendPeriod = 40;
reparametrize = false;
pixelNoise = 1;
noiseFactor = 1;

% randSeeds  = round(10000*rand(1,numRuns));

randSeeds = [8687 8440 3999 2599 8001 4314 9106 1818 2638 1455 1361 8693 5797 5499 1450 8530 6221 3510 5132 4018 7600 2399 1233 1839 2400];
% save [logsDir 'randSeeds.log'] randSeeds -ascii

for l = 1:numel(lmkTypes)
    lmkType = lmkTypes{l};
    for nRun = 1:numRuns
        disp(' ')
        disp('==============================')
        fprintf('Lmk type: %s -- Run #: %d\n',lmkType, nRun);
        disp('==============================')
        logFileName = [logsDir lmkType '-rmse-' num2str(nRun,'%02d') '.log'];
%         logFileName = [lmkType '-' num2str(nRun,'%02d') '.log']
        slamtbSlave;
    end
end

rmsePlots;

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE


text:
function [TorF, vstr, rdate] = have_feature_sdpt3()
%HAVE_FEATURE_SDPT3  Detect availability/version info for SDPT3
%
%   Feature detection function implementing 'sdpt3' tag for HAVE_FEATURE
%   to detect availability/version of SDPT3 SDP solver
%   (https://github.com/sqlp/sdpt3)
%
%   See also HAVE_FEATURE, SDPT3.

%   MP-Opt-Model
%   Copyright (c) 2004-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

TorF = exist('sdpt3','file') == 2;
vstr = '';
rdate = '';
if TorF && have_feature('evalc')
    str = evalc('help sdpt3');
    pat = 'version\s+([^\s]+).*Last Modified: ([^\n]+)\n';
    [s,e,tE,m,t] = regexp(str, pat);
    if ~isempty(t)
        vstr = t{1}{1};
        rdate = datestr(t{1}{2}, 'dd-mmm-yyyy');
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function nPhysioRegressors = tapas_physio_count_physio_regressors(physio)
% Returns number of physiological regressors created, given model
% specification;
% NOTE: only reproducible numbers (data-independent) are
% returned, i.e. session-specific movement spikes and %-variance explained
% PCA-components are not included
%
% IN
%   physio  physio-structure, See also tapas_physio_new
%
% OUT
% nPhysioRegressors     number of physiological regressors, e.g. motion,
%                       retroicor, noise_rois
%                       but: ignores
%
% EXAMPLE
%   tapas_physio_report_contrasts
%
%   See also

% Author: Lars Kasper
% Created: 2014-10-16
% Copyright (C) 2014 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public
% Licence (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/>.

model = physio.model;

nPhysioRegressors = 0;

if model.hrv.include
    nPhysioRegressors = nPhysioRegressors + numel(model.hrv.delays);
end

if model.rvt.include
    nPhysioRegressors = nPhysioRegressors + numel(model.rvt.delays);
end

if model.retroicor.include
    order = model.retroicor.order;
    nPhysioRegressors = nPhysioRegressors + ...
        2*order.c + ...
        2*order.r + ...
        4* order.cr;
end

if model.noise_rois.include
    % TODO: what if number of components implicit?...shall we save this?
    nPhysioRegressors = nPhysioRegressors + ...
        numel(model.noise_rois.roi_files)*ceil(model.noise_rois.n_components+1); % + 1 for mean
end

if model.movement.include
    % TODO: what about variable regressors, that should not be
    % concatenated, e.g. movement outlier censoring
    nPhysioRegressors = nPhysioRegressors + model.movement.order;
end

end


text:
%script to run the ring render gui:

clear
close all

ring_gui;

author: Sable


text:
%kkgeneric 'Generic Routine Interface (glue)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kgeneric.pane file
%
% Parameters:
% String: prefix 'Prefix', default: ' ': 'prefix string'
% InputFile: i 'Input Image ', optional: 'input image'
% OutputFile: o 'Output Image', optional: 'output image or file'
% String: iname1 'Input name #1 ', default: ' ': 'name of input argument #1'
% InputFile: i1 'Input Filename #1 ', optional: 'input filename #1'
% String: iname2 'Input name #2 ', default: ' ': 'name of input argument #2'
% InputFile: i2 'Input Filename #2 ', optional: 'input filename #2'
% String: iname3 'Input name #3 ', default: ' ': 'name of input argument #3'
% InputFile: i3 'Input Filename #3 ', optional: 'input filename #3'
% String: oname1 'Output name #1', default: ' ': 'name of output argument #1'
% OutputFile: o1 'Output Filename #1', optional: 'output filename #1'
% String: oname2 'Output name #2', default: ' ': 'name of output argument #2'
% OutputFile: o2 'Output Filename #2', optional: 'output filename #2'
% String: oname3 'Output name #3', default: ' ': 'name of output argument #3'
% OutputFile: o3 'Output Filename #3', optional: 'output filename #3'
% String: name1 'Arg name #1', default: ' ': 'Name of argument #1'
% String: val1 'Argument #1 value', default: ' ': 'argument #1 value'
% String: name2 'Arg name #2', default: ' ': 'Name of argument #2'
% String: val2 'Argument #2 value', default: ' ': 'argument #2 value'
% String: name3 'Arg name #3', default: ' ': 'Name of argument #3'
% String: val3 'Argument #3 value', default: ' ': 'argument #3 value'
% String: name4 'Arg name #4', default: ' ': 'Name of argument #4'
% String: val4 'Argument #4 value', default: ' ': 'argument #4 value'
% String: name5 'Arg name #5', default: ' ': 'Name of argument #5'
% String: val5 'Argument #5 value', default: ' ': 'argument #5 value'
% String: name6 'Arg name #6', default: ' ': 'Name of argument #6'
% String: val6 'Argument #6 value', default: ' ': 'argument #6 value'
% String: name7 'Arg name #7', default: ' ': 'Name of argument #7'
% String: val7 'Argument #7 value', default: ' ': 'argument #7 value'
% String: name8 'Arg name #8', default: ' ': 'Name of argument #8'
% String: val8 'Argument #8 value', default: ' ': 'argument #8 value'
% String: trailer 'Trailer', default: ' ': 'trailer string'
%
% Example: [o, o1, o2, o3] = kkgeneric({i, i1, i2, i3}, {'prefix',' ';'i','';'o','';'iname1',' ';'i1','';'iname2',' ';'i2','';'iname3',' ';'i3','';'oname1',' ';'o1','';'oname2',' ';'o2','';'oname3',' ';'o3','';'name1',' ';'val1',' ';'name2',' ';'val2',' ';'name3',' ';'val3',' ';'name4',' ';'val4',' ';'name5',' ';'val5',' ';'name6',' ';'val6',' ';'name7',' ';'val7',' ';'name8',' ';'val8',' ';'trailer',' '})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kgeneric - Generic Routine Interface (glue)
%
%  DESCRIPTION
% .I kgeneric
% is a programmable interface that allows a generalized Un*x command line
% to be executed by cantata.
%
% Entries on the pane are spit out to a string that is then system()'d
% off. The string is built up by appending the contents of the SELECTED
% entries together in scan-line order.
%
% Note: All entries are strings. Unfortunately, the forms don't
% seem to like multiple-worded strings, so you have to build up your
% command piece at a time. You CAN use a multiple-worded string in the
% prefix and trailer.
%
% Note: Anything entered into a FILENAME box will be treated as a file
% using the normal Khoros file access checking. This means the UI will
% try to open the filename for access. If it is a Un*x command or other
% object then things won't work.
%
% Note: There appears to be no way to protect special characters from
% the shell.
%
% Note: The oval button to the left of argument name selection allows
% you to put a minus sign ("-") immediately in front of the name itself
% in the command string that is executed. The parser won't let you
% put a minus sign in the string field itself - it will complain about
% unknown arguments. So if you want an argument name that looks like
%
% -blah
%
% then click the button to get the - to show and enter blah as the name.
% The actual command string executed will look like that above.
%
% Note: There is no job control from this routine. Hitting the
% kill button in cantata is not guaranteed to kill the command that was
% system()'d off. Additionally, status reporting is unsupported.
%
%
%
%  EXAMPLES
%
% To execute:
%
% ls -al > snarf
%
% using this routine, do:
%
% % kgeneric -prefix ls -iname1 al -iname2 > -o1 snarf -di1 1 -di2 0 -di3 0 -do1 0
% -do2 0 -do3 0 -dn1 0 -dn2 0 -dn3 0 -dn4 0 -dn5 0 -dn6 0 -dn7 0 -dn8 0
%
% and you should get the desired result. Note that this routine is not much good
% for use on the command line. It is for cantata.
%
%  "SEE ALSO"
% None. Yet.
%
%  RESTRICTIONS
% Kgeneric executes the command via a ksystem() call, and thus there is no
% reporting of errors from the command. If nothing happens when you run a command,
% try  executing the whole command from command line. If the command line fails,
% fix the problem before trying to use kgeneric.
%
% No data type or range checking is available when using this interface.
%
%  REFERENCES
% None.
%
%  COPYRIGHT
% Copyright (C) 1993 - 1996, Khoral Research, Inc.  All rights reserved.
%

function varargout = kkgeneric(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkgeneric(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'prefix', ' ';'i', '__input';'o', '__output';'iname1', ' ';'i1', '__input';'iname2', ' ';'i2', '__input';'iname3', ' ';'i3', '__input';'oname1', ' ';'o1', '__output';'oname2', ' ';'o2', '__output';'oname3', ' ';'o3', '__output';'name1', ' ';'val1', ' ';'name2', ' ';'val2', ' ';'name3', ' ';'val3', ' ';'name4', ' ';'val4', ' ';'name5', ' ';'val5', ' ';'name6', ' ';'val6', ' ';'name7', ' ';'val7', ' ';'name8', ' ';'val8', ' ';'trailer', ' '};
maxval={0,1,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
minval={0,1,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
istoggle=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1];
was_set=istoggle * 0;
paramtype={'String','InputFile','OutputFile','String','InputFile','String','InputFile','String','InputFile','String','OutputFile','String','OutputFile','String','OutputFile','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kgeneric"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [lambda_vec, error_train, error_val] = ...
    validationCurve(X, y, Xval, yval)
%VALIDATIONCURVE Generate the train and validation errors needed to
%plot a validation curve that we can use to select lambda
%   [lambda_vec, error_train, error_val] = ...
%       VALIDATIONCURVE(X, y, Xval, yval) returns the train
%       and validation errors (in error_train, error_val)
%       for different values of lambda. You are given the training set (X,
%       y) and validation set (Xval, yval).
%

% Selected values of lambda (you should not change this)
lambda_vec = [0 0.001 0.003 0.01 0.03 0.1 0.3 1 3 10]';

% You need to return these variables correctly.
error_train = zeros(length(lambda_vec), 1);
error_val = zeros(length(lambda_vec), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return training errors in
%               error_train and the validation errors in error_val. The
%               vector lambda_vec contains the different lambda parameters
%               to use for each calculation of the errors, i.e,
%               error_train(i), and error_val(i) should give
%               you the errors obtained after training with
%               lambda = lambda_vec(i)
%
% Note: You can loop over lambda_vec with the following:
%
%       for i = 1:length(lambda_vec)
%           lambda = lambda_vec(i);
%           % Compute train / val errors when training linear
%           % regression with regularization parameter lambda
%           % You should store the result in error_train(i)
%           % and error_val(i)
%           ....
%
%       end
%
%

for i = 1:length(lambda_vec)
    lambda = lambda_vec(i);
    theta = trainLinearReg(X, y, lambda);
    error_train(i) = linearRegCostFunction(X, y, theta, 0);
    error_val(i) = linearRegCostFunction(Xval, yval, theta, 0);
end

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

end

author: xjwhhh

repo: AndrewNgMachineLearning


text:
function inls = countCorrespondences(M, Btree, epsilon)

    [idx, dist] = knnsearch(Btree, M');
    in_idx = find(dist<=epsilon);

    inls = length(idx(in_idx));

end

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function pos2 = propmo (tjd1, pos1, vel1, tjd2)

% this function applies proper motion, including foreshortening
% effects, to a star's position.

% input

%  tjd1 = tdb julian date of first epoch

%  pos1 = position vector of star at first epoch

%  vel1 = velocity vector of star at first epoch

%  tjd2 = tdb julian date of second epoch

% output

%  pos2 = position vector of star at second epoch

% ported from NOVAS 3.0

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

for j = 1:1:3

    pos2(j) = pos1(j) + vel1(j) * (tjd2 - tjd1);

end

author: goGPS-Project