text:
function [cleanedData, metToKeepSummary] = metsToKeep(dataOrg, metRep, calRelErr, param)
%
% removes the metabolomics measurements from the
% dataOrg table that have low measurement quality based on the relative SD
% of technical replicate sample (RSDqc provided in the metRep variable)
% and/or the average relative error (measured vs actual concentration) of
% the calibration line samples (RE provided in the calRelErr variable)
%
% USAGE:
%   [cleanedData, metToKeepSummary] = metsToKeep(dataOrg, metRep, calRelErr, param)
%
% INPUTS:
%  dataOrg:     A table with original metabolomics data in a long format
%               with information about measured samples, compounds,
%               retention times (RT), area(s), concentrations, etc.
%               the following column is required for the analysis:
%
%                     * .compound - compound (metabolite names) identical
%                                   to the names used in the metRep and
%                                   calRE variables
%  metRep:      A table with quality information based on the dataOrg
%               in a long format (output from mzQuality tool).
%               The following columns are required for the analysis:
%
%                     * .compound - compound (metabolite names) identical
%                                   to the names used in the dataOrg and
%                                   calRE variables
%                     * .RSDqc_* - one or more columns specifying the
%                                  relative standard deviation of the
%                                  repeated measurement of a (quality)
%                                  sample (per batch) per metabolite
%  calRelErr:   A table with quality information about the relative
%               error (RE) of the concentration estimation of the
%               calibration line samples based on the dataOrg
%               in a long format with information about measured samples,
%               compounds, retention times (RT), known concentrations,
%               estimated concentrations, and relative error etc.
%               The following columns are required for the analysis:
%
%                     * .compound - compound (metabolite names) identical
%                                   to the names used in the dataOrg and
%                                   calRE variables
%                     * .RE* - a columns specifying the relative error of
%                              the concentration estimation of the
%                              calibration line samples
%  param.tresholdRSDqc: the treshold value for the relative SD of the repeated
%                       measurement of sample (default = 25 (%); based on
%                       the based on the "Guidelines and considerations for
%                       the use of system suitability and quality control
%                       samples in mass spectrometry assays applied in
%                       untargeted clinical metabolomic studies")
%  param.tresholdCalRE: the treshold value for the relative error of the
%                       concentration estimation of the calibration line
%                       samples (default = 25 (%)
%
% OUTPUTS:
%  cleanedData:         table in the same format as dataOrg without the
%                       metabolites of low quality (above set tresholds)
%  metToKeepSummary:    table in the same format as metRep variable with an
%                       added columns:
%
%                     * .keep -     specifies whether a metabolite is to be
%                                   kept (1) or removed (0) from the further
%                                   analysis
%                     * .REcheck - specifies whether a metabolite passed
%                                   (1) or failed (0) the check based on
%                                   the relative error of the concentration
%                                   estimation of the calibration line sample
%                     * .RSDqcCheck - specifies whether a metabolite passed
%                                     (1) or failed (0) the check based on
%                                     the relative standard deviation of
%                                     the repeated measurement of a (QC) sample
%                     * .sumRE -    shows the average relative error of the
%                                   concentration estimation of the
%                                   calibration line sample per metabolite
%
% EXAMPLE:
%
% NOTE:
%
% Author(s): Agnieszka Wegrzyn (2021)

if exist('param', 'var')
    if isfield(param, "tresholdRSDqc")
        tresholdRSDqc = param.tresholdRSDqc;
    else
        tresholdRSDqc = 25;
    end
    if isfield(param, "tresholdCalRE")
        tresholdCalRE = param.tresholdCalRE;
    else
        tresholdCalRE = 25;
    end
else
    tresholdRSDqc = 25;
    tresholdCalRE = 25;
end

metToKeep = metRep;
calRE = calRelErr;
calRE.RE(isnan(calRE.RE)) = 0;
calRE.RE(isinf(calRE.RE)) = NaN;
metToKeep.REcheck = zeros(length(metToKeep.compound),1);
metToKeep.RSDqcCheck = zeros(length(metToKeep.compound),1);
for i=1:length(metToKeep.compound)
    if sum(ismember(calRE.compound,metToKeep.compound(i))) ~= 0
        avrCalRE = mean(calRE.RE(ismember(calRE.compound,metToKeep.compound(i))),'omitnan');
        sdCalRE = std(calRE.RE(ismember(calRE.compound,metToKeep.compound(i))),'omitnan');
        metToKeep.sumRE(i) = avrCalRE;
        %check if average + 1SD (to account for an outlier) of RE is below the set treshold
        if (avrCalRE + sdCalRE) > tresholdCalRE
            metToKeep.REcheck(i) = 0;
        else
            metToKeep.REcheck(i) = 1;
        end
    else
        metToKeep.sumRE(i) = NaN;
    end
    if all(table2array(metToKeep(i,contains(metToKeep.Properties.VariableNames, 'RSDqc_'))) <= tresholdRSDqc)
        metToKeep.RSDqcCheck(i) = 1;
    else
        metToKeep.RSDqcCheck(i) = 0;
    end
    if metToKeep.RSDqcCheck(i) == 1 && metToKeep.REcheck(i) == 1
        metToKeep.Var1(i) = 1;
    else
        metToKeep.Var1(i) = 0;
    end
end
metToKeepSummary = metToKeep;
metToKeepSummary.Properties.VariableNames(1) = "keep";
cleanedData = dataOrg;
cleanedData(contains(cleanedData.compound,metToKeepSummary.compound(metToKeepSummary.keep == 0)),:)= [];
cleanedData(:,1)= [];
disp(' ')
disp('--------------------------------------------------------------')
disp(['Number of metabolites in the dataset: ' num2str(numel(metToKeepSummary.compound))])
disp(['Number of metabolites left after quality check: ' num2str(sum(metToKeepSummary.keep))])
disp(' ')

author: opencobra

repo: cobratoolbox


text:
classdef spaces
% Utility functions to work with vector spaces

    methods(Static)

        function result  = insersection_space(A, B)
            % Let A and B represent vector spaces
            % The function returns basis for their intersection
            if size(A, 1) ~= size(B, 1)
                error('The ambient space must be same.');
            end

            % Orthogonalize the two bases
            A = orth(A);  % A is n x a
            B = orth(B);  % B is n x b
            % Count the rank of each of them
            % After orthogonalization number of columns may have reduced
            rank_a = size(A, 2);
            rank_b = size(B, 2);
            % Combine them
            C = [A B]; % C is n x (ra + rb)
            % Compute the null space of the whole thing
            D = null(C); % C * D = 0. D is (ra + rb) x  m
            % m is the nullity of C
            % Pick up the first rank_a rows from D
            D2 = D(1:rank_a, :); % ra * m
            % Multiply with A
            result = A * D2 ; % n * m.
            % We get the basis we were looking for
        end

        function result = orth_complement(A, B)
            % Find orthogonal complement of A in B.
            if size(A, 1) ~= size(B, 1)
                error('The ambient space must be same.');
            end
            % Orthogonalize A to make sure that it is full rank.
            A  = orth(A);
            rank_a = size(A,2);
            C = [A B];
            [Q, R] = qr(C, 0);
            % Leave the first rank_a columns from A
            result = Q(:, rank_a+ 1:end);
        end

        function result = principal_angles_orth_cos(A, B)
            % assumes that A and B are orthonormal bases
            % Compute the matrix of inner products of bases of A and B.
            M = A' * B;
            % Compute the SVD of M. The singular values are the cosine of principal angles
            result = svd(M);
            % make sure that result is bounded by 1.
            result = min(1, result);
        end

        function result = principal_angles_cos(A, B)
            % Finds the principal angles between the subspaces spanned by A and B
            % References
            % - http://in.mathworks.com/matlabcentral/newsreader/view_thread/284282
            if size(A, 1) ~= size(B, 1)
                error('The ambient space must be same.');
            end
            A = orth(A);
            B = orth(B);
            result = spx.la.spaces.principal_angles_orth_cos(A, B);
        end

        function result = principal_angles_radian(A, B)
            % Returns the principal angles in radians
            s = spx.la.spaces.principal_angles_cos(A, B);
            % Return the angles as cos inverse of singular values
            % ensure that cos-theta values are <= 1. Otherwise complex numbers may appear.
            result = acos(min(1, s));
        end

        function result = principal_angles_degree(A, B)
            % Returns the principal angles in degrees
            radians = spx.la.spaces.principal_angles_radian(A, B);
            % Return the angles as cos inverse of singular values
            result = rad2deg(radians);
        end

        function result = smallest_angle_cos(A, B)
            % Returns the smallest angle between two subspaces as cos(theta)
            s = spx.la.spaces.principal_angles_cos(A, B);
            result = s(1);
        end

        function result = smallest_angle_orth_cos(A, B)
            % Returns the smallest angle between two subspaces as cos(theta)
            % assumes A and B are orthonormal bases
            s = spx.la.spaces.principal_angles_orth_cos(A, B);
            result = s(1);
        end

        function result = smallest_angle_rad(A, B)
            % Returns the smallest angle between two subspaces in radians
            cos_theta = spx.la.spaces.smallest_angle_cos(A, B);
            result = acos(cos_theta);
        end

        function result = smallest_angle_deg(A, B)
            % Returns the smallest angle between two subspaces in degree
            theta = spx.la.spaces.smallest_angle_rad(A, B);
            result = rad2deg(theta);
        end

        function result = smallest_angles_cos(subspaces, d)
            % subspaces is either a cell array of bases or a concatenated matrix.
            % d is the dimension of each subspace [needed only if all bases are concatenated]
            if iscell(subspaces)
                bases = subspaces;
                % number of subspaces
                s = numel(subspaces);
                % the ambient dimension
                m = size(bases{1}, 1);
            else
                if nargin < 2
                    error('Dimension of each subspace must be specified.');
                end
                [m, n] = size(subspaces);
                if mod(n, d) ~= 0
                    error('n must be multiple of d');
                end
                % number of subspaces
                s = n /d;
                % create the cell array for bases
                bases = cell(s, 1);
                for i=0:s-1
                    %i-th subspace basis
                    si = subspaces(:, i*d + (1:d));
                    bases{i+1} = si;
                end
            end
            % Orthogonalize all subspaces
            for i=1:s
                bases{i} = orth(bases{i});
            end
            % The smallest angles result matrix
            result = eye(s);
            for i=1:s
                si = bases{i};
                for j=i+1:s
                    % j-th subspace
                    sj = bases{j};
                    result(i, j) = spx.la.spaces.smallest_angle_orth_cos(si, sj);
                    result(j, i) = result(i, j);
                end
            end
        end

        function result = smallest_angles_rad(subspaces, d)
            if nargin < 2
                % subspace dimensions are unspecified
                d = -1;
            end
            result = spx.la.spaces.smallest_angles_cos(subspaces, d);
            result = acos(result);
        end

        function result = smallest_angles_deg(subspaces, d)
            if nargin < 2
                % subspace dimensions are unspecified
                d = -1;
            end
            result = spx.la.spaces.smallest_angles_rad(subspaces, d);
            result = rad2deg(result);
        end

        function result = subspace_distance(A, B)
            % The distance between two subspaces based on Grassmannian
            % References
            % - http://math.stackexchange.com/questions/198111/distance-between-real-finite-dimensional-linear-subspaces
            if size(A, 1) ~= size(B, 1)
                error('The ambient space must be same.');
            end
            A = orth(A);
            B = orth(B);
            if size(A, 2) ~= size(B, 2)
                error('The two subspaces must be of same dimensions');
            end
            % Compute the projection matrices for the two spaces
            PA = A*A';
            PB = B*B';
            D = PA  - PB;
            % We return the operator norm of D as the distance between the two subspaces.
            result = norm(D);
        end

        function result = is_in_range_orth(v, U)
            % Returns if v is in the range of U where U is a unitary matrix
            nv = norm(v);
            if nv == 0
                % zero vector is always in the column space of U
                result = true;
                return;
            end
            % Compute the projection of v into the space spanned by U.
            pv = U (U' * v);
            % Compute the difference [projection to the orthogonal complement of U]
            d = v - pv;
            % Compute it's norm
            nd = norm(d);
            % Verify that the norm of difference is indeed very small
            result = nd <= 1e-6 * nv;
        end

        function result = is_in_range(v, A)
            % Returns if v is in the range of an arbitrary matrix A
            result = is_in_range_orth(v, orth(A));
        end

        function result = find_basis(A)
            % Returns a (not necessarily orthogonal) basis of A from columns of A
            [R, pivot_cols] = rref(A);
            [m, n] = size(A);
            i = 1;
            %pivot_cols = false(n, 1);
            %for j=1:n
            %    if R(i, j) == 1
            %        % Add the column containing the leading one
            %        pivot_cols(j) = true;
            %        % Move on to find the 1 in next row
            %        i = i + 1;
            %    end
            %end
            % Return a basis
            result = A (:, pivot_cols);
        end

        function [E, R] = elim(A)
            % References
            % - http://web.mit.edu/18.06/www/Course-Info/Mfiles/elim.m
            % Factorize: E R  = A
            % where E is a product of elementary matrices and
            % R is the row reduced echelon form
            [m, n] = size(A);
            I = eye(m);
            RE = rref([A I]);
            R = RE(:, 1:n);
            E = RE(:, (n+1):m+n);
        end

        function N = null_basis(A)
            % Returns a null basis for A
            % References
            % - http://web.mit.edu/18.06/www/Course-Info/Mfiles/nulbasis.m
            [m, n] = size(A);
            [R, pivot_cols] = rref(A, sqrt(eps));
            r = length(pivot_cols);
            % The columns of A which are not part of pivots
            freecol = 1:n;
            freecol(pivot_cols) = [];
            % Create space for storing the null basis
            N = zeros(n, n-r);
            N(freecol, : ) = eye(n-r);
            N(pivot_cols,  : ) = -R(1:r, freecol);
        end

        function [col_space, null_space, row_space, left_null_space] = four_bases(A)
            % Returns bases for the four spaces associated with the matrix A
            % These are not necessarily orthonormal bases
            [m, n] = size(A);
            [R, pivot_cols] = rref(A, sqrt(eps));
            rank_a = length(pivot_cols);
            % The first rank_a rows of the echelon matrix form the row space
            row_space = R(1:rank_a, :)';
            % The columns of A indexed by the pivot columns form the column space
            col_space = A(:, pivot_cols);

            % Computation of the null space
            % The columns of A which are not part of pivots
            freecol = 1:n;
            freecol(pivot_cols) = [];
            % Allocate memory for the null space basis
            null_space = zeros(n, n-r);
            null_space(freecol, : ) = eye(n-r);
            null_space(pivot_cols,  : ) = -R(1:r, freecol);

            left_null_space = E((r+1):m, :)';
        end

        function [col_space, null_space, row_space, left_null_space] = four_orth_bases(A)
            % Prepares the four bases using SVD
            [U, S, V] = svd(A);
            % Finding the rank of A
            s = diag(S);
            tol = max(size(A))*eps(max(s));
            r = sum(s > tol);
            col_space = U(:, 1:r);
            null_space = U(:, r+1:m);
            row_space = V(:, 1:r);
            left_null_space = V(:, r+1:n);
        end

        function Y = k_dim_to_n_dim(X, n, indices)
            % Maps the data in K dimensions to N-dimensions.
            [k, d] = size(X);
            if nargin < 2
                error('Target dimension must be specified');
            end
            if n < k
                error('n must be larger than k');
            end
            if nargin < 3
                indices = 1:k;
            end
            if ~isvector(indices)
                error('indices must be a vector.');
            end
            if numel(indices) > k
                error('Number of indices must be k');
            end
            if length(unique(indices))<length(indices)
                error('There must be exactly k unique entries in indices');
            end
            if max(indices) > n || min(indices) < 1
                error('Indices cannot point outside 1:n ');
            end
            Y = zeros(n, d);
            Y(indices, :) = X;
        end

        function [A, B, C] = three_spaces_at_angle(N, theta)
            if ~mod(N, 3) == 0
                error('N must be divisible by 3');
            end
            % First create two random orthonormal vectors
            X = orth(randn(N, 3));
            % Then tilt the second one w.r.t. first
            a1 = X(:, 1);
            a2 = X(:, 2);
            a3 = X(:, 3);
            p = cos(theta);
            % first vector for second space
            b1 = sqrt(1 - p^2) * a2 + p * a1;
            % first vector for third space
            c1_1 = p;
            c1_2 = p * (1 - p) / sqrt(1 - p^2);
            c1_3 = sqrt(1 - c1_1^2 - c1_2^2);
            c1 = c1_1 * a1 + c1_2 * a2 + c1_3 * a3;
            X = [a1 b1 c1];
            % Find the orthogonal complement of X
            [U S V] = svd(X);
            Y = U(:, 4:end);
            [~, n] = size(Y);
            % Distribute vectors from Y into A and B
            A = [a1 Y(:, 1:n/3)];
            B = [b1 Y(:, n/3 + (1:n/3))];
            C = [c1 Y(:, 2*n/3 + (1:n/3))];
        end

        function [A, B, C] = three_disjoint_spaces_at_angle(N, theta)
            X = eye(4);
            R1 = [cos(theta) -sin(theta); sin(theta) cos(theta)];
            a1 = [1 ; 0];
            a2 = R1*a1;
            a3 = R1*a2;
            theta2 = deg2rad(59);
            R2 = [cos(theta2) -sin(theta2); sin(theta2) cos(theta2)];
            b1 = [1 ; 0];
            b2 = R2*b1;
            b3 = R2*b2;
            z = zeros(2, 1);
            A = [a1 z; z b1];
            B = [a2 z; z b2];
            C = [a3 z; z b3];
            A =kron(A , eye(N / 2));
            B =kron(B , eye(N / 2));
            C =kron(C , eye(N / 2));
        end

        function describe_three_spaces(A, B, C)
            % Compute principal angles between the subspaces
            fprintf('Ranks: [A]: %d, [B]: %d, [A]: %d\n', ...
                rank(A), rank(B), rank(C));
            fprintf('Cols: [A]: %d, [B]: %d, [A]: %d\n', ...
                size(A, 2), size(B, 2), size(C, 2));
            fprintf('Ranks: [A B]: %d, [B C]: %d, [A C]: %d, \n', ...
                rank([A B]), rank([B C]), rank([A C]));
            D = [A B C];
            fprintf('Rank [A B C]: %d\n', rank(D));
            fprintf('Angle between A and B: %.4f deg\n', spx.la.spaces.smallest_angle_deg(A, B));
            fprintf('Angle between B and C: %.4f deg\n', spx.la.spaces.smallest_angle_deg(B, C));
            fprintf('Angle between A and C: %.4f deg\n', spx.la.spaces.smallest_angle_deg(A, C));
            fprintf('Column wise norms: \n');
            fprintf(' %.2f', spx.norm.norms_l2_cw(A));
            fprintf('\n');
            fprintf(' %.2f', spx.norm.norms_l2_cw(B));
            fprintf('\n');
            fprintf(' %.2f', spx.norm.norms_l2_cw(C));
            fprintf('\n');
        end

        function [A, B, C] = abc_spaces_junk_1(N, theta)
            if ~mod(N, 2) == 0
                error('N must be divisible by 2');
            end
            d = N / 2;
            % First create three random orthonormal vectors
            X = orth(randn(N, 2));
            % Then tilt the second one w.r.t. first
            a1 = X(:, 1);
            a2 = X(:, 2);
            p = cos(theta);
            % first vector for second space
            b1 = sqrt(1 - p^2) * a2 + p * a1;
            X = [a1 b1];
            % Find the orthogonal complement of X
            [U S V] = svd(X);
            Y = U(:, 3:end);
            [~, n] = size(Y);
            % Distribute vectors from Y into A and B
            A = [a1 Y(:, 1:n/2)];
            B = [b1 Y(:, n/2 + 1:end)];

            % choose a vector a3 which is orthogonal to A and b1
            X = [A B(:, 1:d-1)];
            [U S V] = svd(X);
            a3 = U(:, size(X, 2) + 1);
            % first vector for third space
            c1_1 = p;
            c1_2 = p * (1 - p) / sqrt(1 - p^2);
            c1_3 = sqrt(1 - c1_1^2 - c1_2^2);
            c1 = c1_1 * a1 + c1_2 * a2 + c1_3 * a3;

            X = [A c1];
            % Find the orthogonal complement of X
            [U S V] = svd(X);
            Y = U(:, size(X, 2)+1:end);
            C = [c1 Y];
        end

        function result = have_same_column_spans(A, B)
            % Checks if the column spans of two matrices are same.
            r1 = rank(A);
            r2 = rank(B);
            if r1 ~= r2
                result = false;
                return;
            end
            r3 = rank([A B]);
            if r3 ~= r1
                result = false;
                return;
            end
            result = true;
        end

    function result =  bases(X, counts)
        % bases for individual subspaces
        K = length(counts);
        % bases cell array
        result = cell(1, K);
        [start_indices, end_indices] = spx.cluster.start_end_indices(counts);
        for k=1:K
            ss = start_indices(k);
            ee = end_indices(k);
            XX = X(:, ss:ee);
            basis = orth(XX);
            result{k} = basis;
        end
    end

    end

end

author: indigits

repo: sparse-plex


text:
function update_temporal_parallel(obj, use_parallel, use_c_hat)
%% update the the temporal components for all neurons
% input:
%   use_parallel: boolean, do initialization in patch mode or not.
%       default(true); we recommend you to set it false only when you want to debug the code.
%   use_c_hat: use the previous estimation of C

%% Author: Pengcheng Zhou, Columbia University, 2017
%% email: zhoupc1988@gmail.com

%% process parameters

try
    % map data
    mat_data = obj.P.mat_data;

    % folders and files for saving the results
    log_file =  obj.P.log_file;
    flog = fopen(log_file, 'a');
    log_data = matfile(obj.P.log_data, 'Writable', true); %#ok<NASGU>

    % dimension of data
    dims = mat_data.dims;
    d1 = dims(1);
    d2 = dims(2);
    T = dims(3);
    obj.options.d1 = d1;
    obj.options.d2 = d2;

    % parameters for patching information
    patch_pos = mat_data.patch_pos;
    block_pos = mat_data.block_pos;

    % number of patches
    [nr_patch, nc_patch] = size(patch_pos);
catch
    error('No data file selected');
end
fprintf('\n-----------------UPDATE TEMPORAL---------------------------\n');

% frames to be loaded
frame_range = obj.frame_range;
T = diff(frame_range) + 1;

% threshold for detecting large residuals
thresh_outlier = obj.options.thresh_outlier;

% use parallel or not
if ~exist('use_parallel', 'var')||isempty(use_parallel)
    use_parallel = true; %don't save initialization procedure
end

% use previous estimation of C
if ~exist('use_c_hat', 'var') || isempty(use_c_hat)
    use_c_hat = true;
end
% options
options = obj.options;
bg_model = options.background_model;
maxIter = options.maxIter;

%% identify existing neurons within each patch
A = cell(nr_patch, nc_patch);
C = cell(nr_patch, nc_patch);
% if strcmpi(bg_model, 'ring')
    A_prev = A;
    C_prev = C;
% end
sn = cell(nr_patch, nc_patch);
ind_neurons = cell(nr_patch, nc_patch);

AA = cell(nr_patch, nc_patch);   % save the ai^T*ai for each neuron

for mpatch=1:(nr_patch*nc_patch)
    if strcmpi(bg_model, 'ring')
        tmp_block = block_pos{mpatch};
    else
        tmp_block = patch_pos{mpatch};
    end
    % find the neurons that are within the block
    mask = zeros(d1, d2);
    mask(tmp_block(1):tmp_block(2), tmp_block(3):tmp_block(4)) = 1;
    ind = (reshape(mask(:), 1, [])* obj.A>0);
    A{mpatch}= obj.A(logical(mask), ind);
    sn{mpatch} = obj.P.sn(logical(mask));
    C{mpatch} = obj.C(ind, :);
    ind_neurons{mpatch} = find(ind);    % indices of the neurons within each patch

    if strcmpi(bg_model, 'ring')
        ind = find(reshape(mask(:)==1, 1, [])* full(obj.A_prev)>0);
        A_prev{mpatch}= obj.A_prev((mask>0), ind);
        C_prev{mpatch} = obj.C_prev(ind, :);
    end
end
%% prepare for the variables for computing the background.
bg_model = obj.options.background_model;
bg_ssub = obj.options.bg_ssub;
W = obj.W;
b0 = obj.b0;
b = obj.b;
f = obj.f;

%% start updating temporal components
C_raw_new = C;
deconv_flag = obj.options.deconv_flag;
if options.deconv_flag
    deconv_options = obj.options.deconv_options;
else
    deconv_options = [];
end
if use_parallel
    parfor mpatch=1:(nr_patch*nc_patch)
        % no neurons within the patch
        [r, c] = ind2sub([nr_patch, nc_patch], mpatch);
        tmp_patch = patch_pos{mpatch};     %[r0, r1, c0, c1], patch location
        if strcmpi(bg_model, 'ring')
            tmp_block = block_pos{mpatch};
        else
            tmp_block = patch_pos{mpatch};
        end
        C_patch = C{mpatch};                % previous estimation of neural activity

        if isempty(C_patch)
            fprintf('Patch (%2d, %2d) is done. %2d X %2d patches in total. \n', r, c, nr_patch, nc_patch);
            continue;
        end
        A_patch = A{mpatch};

        % use ind_patch to indicate pixels within the patch
        ind_patch = false(diff(tmp_block(1:2))+1, diff(tmp_block(3:4))+1);
        ind_patch((tmp_patch(1):tmp_patch(2))-tmp_block(1)+1, (tmp_patch(3):tmp_patch(4))-tmp_block(3)+1) = true;

        % get data
        if strcmpi(bg_model, 'ring')
            % including areas outside of the patch for recorving background
            % in the ring model
            Ypatch = get_patch_data(mat_data, tmp_patch, frame_range, true);
        else
            Ypatch = get_patch_data(mat_data, tmp_patch, frame_range, false);
        end
        [nr_block, nc_block, ~] = size(Ypatch);

        % get background
        if strcmpi(bg_model, 'ring')
            A_patch_prev = A_prev{mpatch};
            C_patch_prev = C_prev{mpatch};
            W_ring = W{mpatch};
            b0_ring = b0{mpatch};
            Ypatch = reshape(Ypatch, [], T);
            tmp_Y = double(Ypatch)-A_patch_prev*C_patch_prev;
            if bg_ssub==1
                Ypatch = bsxfun(@minus, double(Ypatch(ind_patch,:))- W_ring*tmp_Y, b0_ring-W_ring*mean(tmp_Y, 2));
            else
                % get the dimension of the downsampled data
                [d1s, d2s] = size(imresize(zeros(nr_block, nc_block), 1/bg_ssub));
                % downsample data and reconstruct B^f
                temp = reshape(bsxfun(@minus, tmp_Y, mean(tmp_Y, 2)), nr_block, nc_block, []);
                temp = imresize(temp, 1./bg_ssub);
                Bf = reshape(W_ring*reshape(temp, [], T), d1s, d2s, T);
                Bf = imresize(Bf, [nr_block, nc_block]);
                Bf = reshape(Bf, [], T);

                Ypatch = bsxfun(@minus, double(Ypatch(ind_patch, :)) - Bf(ind_patch, :), b0_ring);
            end
        elseif strcmpi(bg_model, 'nmf')
            b_nmf = b{mpatch};
            f_nmf = f{mpatch};
            Ypatch = double(reshape(Ypatch, [], T))- b_nmf*f_nmf;
        else
            b_svd = b{mpatch};
            f_svd = f{mpatch};
            b0_svd = b0{mpatch};
            Ypatch = double(reshape(Ypatch, [], T)) - bsxfun(@plus, b_svd*f_svd, b0_svd);
        end

        % using HALS to update temporal components
        if ~use_c_hat
            [AA{mpatch}, C_raw_new{mpatch}] = fast_temporal(Ypatch, A_patch(ind_patch,:));
        else
            [~, C_raw_new{mpatch}] = HALS_temporal(Ypatch, A_patch(ind_patch,:), C_patch, maxIter, deconv_options);
            AA{mpatch}= sum(A_patch(ind_patch,:).^2, 1);
        end

        fprintf('Patch (%2d, %2d) is done. %2d X %2d patches in total. \n', r, c, nr_patch, nc_patch);
    end
else
    for mpatch=1:(nr_patch*nc_patch)
        % no neurons within the patch
        [r, c] = ind2sub([nr_patch, nc_patch], mpatch);
        tmp_patch = patch_pos{mpatch};     %[r0, r1, c0, c1], patch location
        if strcmpi(bg_model, 'ring')
            tmp_block = block_pos{mpatch};
        else
            tmp_block = patch_pos{mpatch};
        end

        C_patch = C{mpatch};                % previous estimation of neural activity

        if isempty(C_patch)
            fprintf('Patch (%2d, %2d) is done. %2d X %2d patches in total. \n', r, c, nr_patch, nc_patch);
            continue;
        end
        A_patch = A{mpatch};

        % use ind_patch to indicate pixels within the patch
        ind_patch = false(diff(tmp_block(1:2))+1, diff(tmp_block(3:4))+1);
        ind_patch((tmp_patch(1):tmp_patch(2))-tmp_block(1)+1, (tmp_patch(3):tmp_patch(4))-tmp_block(3)+1) = true;

        % get data
        if strcmpi(bg_model, 'ring')
            % including areas outside of the patch for recorving background
            % in the ring model
            Ypatch = get_patch_data(mat_data, tmp_patch, frame_range, true);
        else
            Ypatch = get_patch_data(mat_data, tmp_patch, frame_range, false);
        end
        [nr_block, nc_block, ~] = size(Ypatch);

        % get background
        if strcmpi(bg_model, 'ring')
            A_patch_prev = A_prev{mpatch};
            C_patch_prev = C_prev{mpatch};
            W_ring = W{mpatch};
            b0_ring = b0{mpatch};
            Ypatch = reshape(Ypatch, [], T);
            tmp_Y = double(Ypatch)-A_patch_prev*C_patch_prev;
            if bg_ssub==1
                Ypatch = bsxfun(@minus, double(Ypatch(ind_patch,:))- W_ring*tmp_Y, b0_ring-W_ring*mean(tmp_Y, 2));
            else
                % get the dimension of the downsampled data
                [d1s, d2s] = size(imresize(zeros(nr_block, nc_block), 1/bg_ssub));
                % downsample data and reconstruct B^f
                temp = reshape(bsxfun(@minus, tmp_Y, mean(tmp_Y, 2)), nr_block, nc_block, []);
                temp = imresize(temp, 1./bg_ssub);
                Bf = reshape(W_ring*reshape(temp, [], T), d1s, d2s, T);
                Bf = imresize(Bf, [nr_block, nc_block]);
                Bf = reshape(Bf, [], T);

                Ypatch = bsxfun(@minus, double(Ypatch(ind_patch, :)) - Bf(ind_patch, :), b0_ring);
            end
        elseif strcmpi(bg_model, 'nmf')
            b_nmf = b{mpatch};
            f_nmf = f{mpatch};
            Ypatch = double(reshape(Ypatch, [], T))- b_nmf*f_nmf;
        else
            b_svd = b{mpatch};
            f_svd = f{mpatch};
            b0_svd = b0{mpatch};
            Ypatch = double(reshape(Ypatch, [], T)) - bsxfun(@plus, b_svd*f_svd, b0_svd);
        end

        % using HALS to update temporal components
        if ~use_c_hat
            [AA{mpatch}, C_raw_new{mpatch}] = fast_temporal(Ypatch, A_patch(ind_patch,:));
        else
            [~, C_raw_new{mpatch}] = HALS_temporal(Ypatch, A_patch(ind_patch,:), C_patch, maxIter, deconv_options);
            AA{mpatch}= sum(A_patch(ind_patch,:).^2, 1);
        end

        fprintf('Patch (%2d, %2d) is done. %2d X %2d patches in total. \n', r, c, nr_patch, nc_patch);
    end
end
%% collect results
K = size(obj.C, 1);
C_new = zeros(K, T);
aa = zeros(K, 1);
fprintf('Collect results from all small patches...\n');
for mpatch=1:(nr_patch*nc_patch)
    C_raw_patch = C_raw_new{mpatch};
    ind_patch = ind_neurons{mpatch};
    aa_patch = AA{mpatch};
    for m=1:length(ind_patch)
        k = ind_patch(m);
        C_new(k, :) = C_new(k, :) + C_raw_patch(m,:) * aa_patch(m);
        aa(k) = aa(k) + aa_patch(m);
    end
end
aa(aa==0) = 1;
obj.C_raw = bsxfun(@times, C_new, 1./aa);
fprintf('Deconvolve and denoise all temporal traces again...\n');
if deconv_flag
    obj.C = obj.deconvTemporal();
else
    obj.C_raw = bsxfun(@minus, obj.C_raw, min(obj.C_raw,[],2));
    obj.C = obj.C_raw;
end
fprintf('Done!\n');

%% upadte b0
if strcmpi(bg_model, 'ring')
    fprintf('Update the constant baselines for all pixels..\n');
    obj.b0_new = cell2mat(obj.P.Ymean)-obj.reshape(obj.A*mean(obj.C,2), 2); -obj.reconstruct_b0();
    fprintf('Done!\n');
end

%% save the results to log
fprintf(flog, '[%s]\b', get_minute());
fprintf(flog, 'Finished updating temporal components.\n');
if obj.options.save_intermediate
    temporal.C_raw = obj.C_raw;
    temporal.ids = obj.ids;
    temporal.C = obj.C;
    temporal.S = obj.S;
    temporal.P.kernel_pars = obj.P.kernel_pars;
    temporal.b0_new = obj.b0_new;
    tmp_str = get_date();
    tmp_str=strrep(tmp_str, '-', '_');
    eval(sprintf('log_data.temporal_%s = temporal;', tmp_str));
    fprintf(flog, '\tThe results were saved as intermediate_results.temporal_%s\n\n', tmp_str);
end
fclose(flog);

function [aa, C_raw] = fast_temporal(Y, A)
%% estimate temporal components using the mean fluorescence
% input:
%   Y:  d*T, fluorescence data
%   A:  d*K, spatial components
% output:
%   aa: K*1, energy of each Ai
%   C_raw: K*T, the temporal components before thresholded or being
%   denoised.
% Author: Pengcheng Zhou, Columbia University, 2017
% zhoupc1988@gmail.com

%% options

%% initialization
% roughly initialize C
tmpA = bsxfun(@times, A, 1./max(A, [], 1));
ind_max = bsxfun(@ge, tmpA, 0.5); %max(tmpA, [], 2));
tmp_A = A.*double(ind_max);
aa = sum(tmp_A.^2, 1);
ind = (aa==0);
aa(ind) = inf;
C_raw = bsxfun(@times, tmp_A'*Y, 1./aa');
aa(ind) = 0;
% if any(ind)  % explain the residual using weak neurons
%     tmp_A = A(:, ind);
%     C_raw(ind, :) = (tmp_A'*tmp_A)\(tmp_A'*Y - tmp_A'*A*C_raw);
%     aa(ind) = sum(tmp_A.^2, 1);
% end

author: zhoupc


text:
function [grRules_new,genes,rxnGeneMat] = translateGrRules(grRules,targetFormat,origFormat,noMatch)
%translateGrRules  Translate grRules to other other gene ID types.
%
% translateGrRules converts grRules to a chosen target gene ID type(s).
% The original and target gene type(s) must be one of the following:
%
%       'ENSG'     Ensembl gene ID
%       'ENST'     Ensembl transcript ID
%       'ENSP'     Ensembl protein ID
%       'UniProt'  UniProt protein ID
%       'Symbol'   HUGO gene abbreviation
%       'Entrez'   Entrez (NCBI) gene ID
%
% The function will try to determine the original gene ID format, but this
% can be manually specified with the "origFormat" input. Furthermore, the
% user can supply a custom conversion key in the "targetFormat" input that
% uses gene IDs other than the above built-in ID types (see details below).
%
% NOTE: All trailing ".#" after Entrez or Ensembl ID (e.g., ENSG0000123.1)
%       will be REMOVED. The grRules may therefore change; e.g.:
%           '123.1 or 123.2 or 4365.1' --> '123 or 4365'
%
% NOTE: The function also "cleans" the grRules, meaning duplicate genes,
%       extra parentheses, etc. will be removed from each grRule. See
%       "cleanModelGeneRules.m" for more details.
%
%
% USAGE:
%
%   [grRules_new,genes,rxnGeneMat] = translateGrRules(grRules,targetFormat,origFormat,noMatch);
%
%
% INPUTS:
%
%   grRules     A cell vector of gene-reaction rules (e.g., model.grRules).
%
%   targetFormat    The desired output format of the grRules. This can be
%                   one or multiple gene ID types, and must be selected
%                   from the following: 'ENSG', 'ENSP', 'ENST', 'UniProt',
%                   'Symbol', 'Entrez'.
%                   DEFAULT: 'Symbol'.
%
%               *** NOTE: As a special case, instead of entering a target
%                   format or list of formats, this input can be replaced
%                   with a custom conversion key. This conversion key must
%                   be an Nx2 cell array, where the first column contains
%                   gene IDs corresponding to the current grRules, and the
%                   second column contains new gene IDs to which the rules
%                   will be translated.
%
%   origFormat  (Optional) A string specifying the original gene ID type
%               (e.g., 'UniProt'). If not provided, the function will try
%               to automatically determine the type of IDs present.
%
%   noMatch     What to do in case a gene in the original grRules does not
%               have a corresponding match in the new gene IDs to which the
%               grRules are to be translated. Choose from the following:
%
%               'delete'    (DEFAULT) Delete (exclude) the gene from the
%                           new rule.
%
%               'original'  keep the original gene ID. Note, this will
%                           result in mixed gene ID types in grRules and
%                           the corresponding list of genes, which may
%                           cause problems in other applications.
%
%
% OUTPUTS:
%
%   grRules_new  grRules with the gene IDs converted to the target format.
%                If targetFormat is more than one ID type, then grRules
%                will be returned as a structure, with a different field
%                for each gene ID type.
%
%   genes        A list of all genes found in the converted grRules.
%                If targetFormat is more than one ID type, then genes will
%                be returned as a structure, with a different field for
%                each gene ID type.
%
%   rxnGeneMat   A matrix associating reactions to genes, build from the
%                converted grRules and genes list. If targetFormat is more
%                than one ID type, then rxnGeneMat will be returned as a
%                structure, with a different field for each gene ID type.
%

% handle input arguments
if nargin < 4
    noMatch = 'delete';
elseif ~ismember(lower(noMatch),{'original','delete'})
    error('Valid inputs for noMatch are "original" or "delete".');
end

if nargin < 3
    gene_type_orig = {};
else
    gene_type_orig = origFormat;
end

custom_key = false;
if nargin < 2 || isempty(targetFormat)
    targetFormat = 'Symbol';
elseif ischar(targetFormat)
    % convert target format from string to cell
    targetFormat = {targetFormat};
elseif (size(targetFormat,2) == 2) && (size(targetFormat,1) > 10)  % check that it has some arbitrarly long length above 10
    % use custom conversion key
    custom_key = true;
    conv_key = targetFormat;
    targetFormat = {'NewFormat'};
    conv_key_head = {'OldFormat','NewFormat'};
    gene_type_orig = 'OldFormat';
else
    % ensure it is a column vector
    targetFormat = targetFormat(:);
end

% allow for non-exact matches to "Symbol" option
targetFormat(ismember(lower(targetFormat),{'name','names','symbol','symbols'})) = {'Symbol'};
if ismember(lower(gene_type_orig), {'name','names','symbol','symbols'})
    gene_type_orig = 'Symbol';
end

% get original list of genes from the grRules
genes_orig = getGenesFromGrRules(grRules);
if ismember('and',genes_orig) || ismember('or',genes_orig)
    error('Problem reading grRules. Verify that all "and" and "or" elements are lowercase and surrounded by spaces.');
end

% preprocess gene list and gene-reaction rules, if necessary
if any(ismember({'GPI','GAPDH'},genes_orig))  % need a more robust method to check if it's gene symbols
    % check if the list is gene symbols, if so, do not modify
    rules_orig = grRules;
    if ~custom_key && isempty(gene_type_orig)
        gene_type_orig = 'Symbol';
    end
else
    % remove ".#" from gene IDs (e.g. ENSG00000198888.2 -> ENSG00000198888, or 2597.1 -> 2597)
    genes_orig = regexprep(genes_orig,'\.\d+$','');
    rules_orig = regexprep(grRules,'\.\d+','');
end

% determine the original gene ID type, if not gene symbols
if isempty(gene_type_orig)
    if startsWith(genes_orig{1},{'ENSG','ENST','ENSP'})
        % ensembl ID type
        gene_type_orig = genes_orig{1}(1:4);
    elseif all(ismember(cellfun(@length,genes_orig),[6,10])) && all(cellfun(@(x) ~isempty(regexp(x,'^[A-Z]\d\w\w\w\d','once')),genes_orig))
        % Uniprot IDs are 6 or 10 characters, with format [A-Z][0-9][A-Z,0-9][A-Z,0-9][A-Z,0-9][0-9]...
        gene_type_orig = 'UniProt';
    elseif ~any(isnan(str2double(genes_orig)))
        % Entrez IDs are only digits
        gene_type_orig = 'Entrez';
    else
        error('Unable to recognize the type of gene IDs in the model.');
    end
end

% remove duplicate genes (if any)
genes_orig = unique(genes_orig,'stable');

if ~custom_key
    % import ID conversion key from the model "genes.tsv" annotation file

    % get the path
    [ST, I] = dbstack('-completenames');
    path = fileparts(ST(I).file);
    tmpfile = fullfile(path,'..','..','model','genes.tsv');

    % import as structure, convert to table, and extract header
    tmp = struct2table(importTsvFile(tmpfile));
    conv_key = table2array(tmp);
    conv_key_head = tmp.Properties.VariableNames;
    clear tmp

    % change header names to match contents of GENE_TYPES
    [hasMatch,ind] = ismember(conv_key_head,{'genes','geneENSTID',...
        'geneENSPID','geneUniProtID','geneSymbols','geneEntrezID'});
    type_abbrevs = {'ENSG';'ENST';'ENSP';'UniProt';'Symbol';'Entrez'};
    conv_key_head = type_abbrevs(ind(hasMatch));
end

% begin by "cleaning" the original grRules
rules_orig = cleanGrRules(rules_orig);

% convert rules to all other gene ID types
for i = 1:length(targetFormat)

    % extract portion of conversion key required; resulting variable will
    % be two columns of IDs, where column 1 is the original gene ID type,
    % and column 2 is the new gene ID type.
    [~,ind] = ismember([{gene_type_orig}, targetFormat(i)], conv_key_head);
    conv_key_sub = conv_key(:,ind);

    % remove rows with empty key entries or unneeded genes
    conv_key_sub(any(cellfun(@isempty,conv_key_sub),2),:) = [];
    conv_key_sub(~ismember(conv_key_sub(:,1),genes_orig),:) = [];

    % expand entries with multiple IDs (separated by semicolons) into
    % multiple rows with one ID per entry
    nested_key = cellfun(@(x) strip(strsplit(x, ';'))', conv_key_sub, 'UniformOutput', false);
    dup_col2 = arrayfun(@(i) repmat(nested_key(i,2), numel(nested_key{i,1}), 1), [1:size(nested_key,1)]', 'UniformOutput', false);
    semi_nested_key = [vertcat(nested_key{:,1}), vertcat(dup_col2{:})];
    dup_col1 = arrayfun(@(i) repmat(semi_nested_key(i,1), numel(semi_nested_key{i,2}), 1), [1:size(semi_nested_key,1)]', 'UniformOutput', false);
    conv_key_sub = [vertcat(dup_col1{:}), vertcat(semi_nested_key{:,2})];

    % only keep unique rows of conversion key
    [~,tmp] = ismember(conv_key_sub,unique(conv_key_sub));
    [~,ind] = unique(tmp,'rows');
    conv_key_sub = conv_key_sub(ind,:);

    if strcmpi(noMatch,'orig')
        % Identify genes that do not exist in the conversion key. For these
        % genes, add additional rows to the conversion key such that if an
        % original gene ID has no corresponding new gene ID, the original
        % gene ID will be used instead.
        unmatched = genes_orig(~ismember(genes_orig,conv_key_sub(:,1)));
        conv_key_sub = [conv_key_sub;repmat(unmatched,1,2)];
    end

    % define inline function to convert between gene IDs
    convertGene = @(g) strjoin(conv_key_sub(ismember(conv_key_sub(:,1),g),2), ' | ');

    % convert "and" to "&" and "or" to "|" (easier to work with)
    rules_new = rules_orig;
    rules_new = regexprep(rules_new, ' or ', ' | ');
    rules_new = regexprep(rules_new, ' and ', ' & ');

    % check if the target format is the original format
    if strcmp(targetFormat{i},gene_type_orig)
        % if current model rules are already in the target format, do not
        % convert gene IDs
    else
        % convert gene IDs
        % This next line identifies gene IDs as collections of characters that
        % are not spaces, parentheses, or the symbols & or |. It then replaces
        % those gene IDs with the new gene ID type, which calls the convertGene
        % inline function to retrieve.
        rules_new = regexprep(rules_new, '[^&|\(\) ]+', '(${convertGene($0)})');
    end

    % clean up rules (removes extra parentheses, repeated genes, etc.)
    rules_new = cleanGrRules(rules_new);

    % generate new rxnGeneMat and gene list based on converted grRules
    [genes.(targetFormat{i}),rxnGeneMat.(targetFormat{i})] = getGenesFromGrRules(rules_new);

    % restore "&" as "and" and "|" as "or"
    rules_new = regexprep(rules_new, ' \| ', ' or ');
    rules_new = regexprep(rules_new, ' & ', ' and ');

    % add new rules to model structure
    grRules_new.(targetFormat{i}) = rules_new;

end

% if only one target format is specified, structures are not necessary
if length(targetFormat) == 1
    grRules_new = grRules_new.(targetFormat{1});
    if ~isempty(genes)
        genes = genes.(targetFormat{1});
        rxnGeneMat = rxnGeneMat.(targetFormat{1});
    end
end

author: SysBioChalmers

repo: Human-GEM


text:
% APM Load Data File
%
% response = csv_load(server,app,filename)
%
% Function csv_load uploads the data file (csv) to the web-server
%   with the following arguments:
%
%   server = address of server
%      app = application name
% filename = data filename
%
% A response is returned indicating whether the file was
%   uploaded successfully to the server
%
function [response] = csv_load(server,app,filename)
   % load model
   fid=fopen(filename,'r');
   tline = [];
   while 1
      aline = fgets(fid);
      if ~ischar(aline), break, end
      % remove any double quote marks
      aline = [strrep(aline,'"',' ')];
      tline = [tline aline];
   end
   fclose(fid);

   % send to server once for every 2000 characters
   ts = size(tline,2);
   block = 2000;
   cycles = ceil(ts/block);
   for i = 1:cycles,
      if i<cycles,
         csv_block = ['csva ' tline((i-1)*block+1:i*block)];
      else
         csv_block = ['csv ' tline((i-1)*block+1:end)];
      end
      response = apm(server,app,csv_block);
   end
   response = 'Successfully loaded CSV file';

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


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)
    [theta] = trainLinearReg(X,y,lambda_vec(i));
    [error_train(i),~] = linearRegCostFunction(X,y,theta,0);
    [error_val(i),~] = linearRegCostFunction(Xval,yval,theta,0);
end;

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

end

repo: machine-learning-coursera


text:
function applyFcn= misc_getApplyFunc(model)
%applyFcn= misc_getApplyFunc(model)

if isstruct(model),
  func= misc_getFuncParam(model.classy);
else
  func= misc_getFuncParam(model);
end

trainFcnName= func2str(func);
if length(trainFcnName)<7 || ~strncmp('train_', trainFcnName, 6),
  error('names of classifier functions need to be prefixed by ''train_''.');
end

baseName= trainFcnName(7:end);
applyFcnName= ['apply_' baseName];
if exist(applyFcnName, 'file'),
  applyFcn= str2func(applyFcnName);
else
  applyFcn= @apply_separatingHyperplane;
end

author: bbci


text:
function [fs, bmg] = slfiltersize(fs0)
%SLFILTERSIZE Extracts information from filtersize
%
% $ Syntax $
%   - [fs, bmg] = slfiltersize(fs0)
%
% $ Arguments $
%   - fs0:      The input filter size
%   - fs:       The full filter size form
%   - bmg:      The boundary margins
%
% $ Description $
%   - [fs, bmg] = slfiltersize(fs0) restores the full form of the input
%     filtersize. In sltoolbox, filter size can be specified in either
%     of the following forms:
%     \*
%     \t    Table.  The forms of the filter size                    \\
%     \h     name      &           syntax                           \\
%            full      & [height, width, center_y, center_x]        \\
%            sizeonly  & [height, width]
%                        The center will be computed as:
%                        cy = floor((1 + h) / 2)
%                        cx = floor((1 + w) / 2)                    \\
%            lenonly   & [len]
%                        height = width = len
%     \*
%     bmg is the boundary margins in the form of
%     [top_margin, bottom_margin, left_margin, right_margin]
%
% $ History $
%   - Created by Dahua Lin, on Sep 1st, 2006
%

%% parse filter size

if ~isvector(fs0)
    error('sltoolbox:invalidarg', ...
        'fs0 should be a vector');
end

switch length(fs0)
    case 1
        h = fs0;
        w = fs0;
        cy = floor((1+h)/2);
        cx = floor((1+w)/2);
    case 2
        cencoords = floor((1 + fs0) / 2);
        h = fs0(1);
        w = fs0(2);
        cy = cencoords(1);
        cx = cencoords(2);
    case 4
        h = fs0(1);
        w = fs0(2);
        cy = fs0(3);
        cx = fs0(4);
    otherwise
        error('sltoolbox:sizmismatch', ...
            'The length of fs0 is illegal');
end

fs = [h, w, cy, cx];

%% compute boundary margins

if nargout >= 2
    tm = cy - 1;
    bm = h - cy;
    lm = cx - 1;
    rm = w - cx;
    bmg = [tm, bm, lm, rm];
end

author: lmthang


text:
%MPQ_POLY Polygon moments
%
% M = MPQ_POLY(V, P, Q) is the PQ'th moment of the polygon with vertices
% described by the columns of V.
%
% Notes::
% - The points must be sorted such that they follow the perimeter in
%   sequence (counter-clockwise).
% - If the points are clockwise the moments will all be negated, so centroids
%   will be still be correct.
% - If the first and last point in the list are the same, they are considered
%   to be a single vertex.
%
% See also MPQ, NPQ_POLY, UPQ_POLY, Polygon.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB 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.
%
% MVTB 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 MVTB.  If not, see <http://www.gnu.org/licenses/>.

function m = mpq(iv, p, q)
    if ~all(iv(:,1) == iv(:,end))
        %disp('closing the polygon')
        iv = [iv iv(:,1)];
    end
    [nr,n] = size(iv);
    if nr < 2,
        error('must be at least two rows of data')
    end
    x = iv(1,:);
    y = iv(2,:);

    m = 0.0;
    for l=1:n
        if l == 1
            dxl = x(l) - x(n);
            dyl = y(l) - y(n);
        else
            dxl = x(l) - x(l-1);
            dyl = y(l) - y(l-1);
        end
        Al = x(l)*dyl - y(l)*dxl;

        s = 0.0;
        for i=0:p
            for j=0:q
                s = s + (-1)^(i+j) * combin(p,i) * combin(q,j)/(i+j+1) * x(l)^(p-i)*y(l)^(q-j) * dxl^i * dyl^j;
            end
        end
        m = m + Al * s;
    end
    m = m / (p+q+2);

function c = combin(n, r)
%
% COMBIN(n,r)
%   compute number of combinations of size r from set n
%
    c = prod((n-r+1):n) / prod(1:r);

author: petercorke

repo: machinevision-toolbox-matlab


text:
function results = vl_test_ihashsum(varargin)
% VL_TEST_IHASHSUM
vl_test_init ;

function s = setup()
rand('state',0) ;
s.data = uint8(round(16*rand(2,100))) ;
sel = find(all(s.data==0)) ;
s.data(1,sel)=1 ;

function test_hash(s)
D = size(s.data,1) ;
K = 5 ;
h = zeros(1,K,'uint32') ;
id = zeros(D,K,'uint8');
next = zeros(1,K,'uint32') ;
[h,id,next] = vl_ihashsum(h,id,next,K,s.data) ;

sel = vl_ihashfind(id,next,K,s.data) ;
count = double(h(sel)) ;

[drop,i,j] = unique(s.data','rows') ;
for k=1:size(s.data,2)
  count_(k) = sum(j == j(k)) ;
end
vl_assert_equal(count,count_) ;

author: Cloud-CV

repo: object-proposals


text:
function fibonacci_spiral ( n )

%*****************************************************************************80
%
%% FIBONACCI_SPIRAL draws points on a Fibonacci spiral.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points to plot.
%    Default is 101.
%
  if ( nargin < 1 )
    n = 101;
  end
%
%  PHI is the golden ratio, the limit of the ratio of
%  successive Fibonacci numbers:
%
%    PHI = limit ( N->oo ) F(N+1)/F(N)
%
  phi = ( 1.0 + sqrt ( 5.0 ) ) / 2.0;
%
%  Allocate storage for the data.
%
  x = zeros ( n, 1 );
  y = zeros ( n, 1 );
%
%  Set the angle and radius of the first point.
%
  a = 0.0;
  r = 0.0;
%
%  Set the increments.
%
  da = 2.0 * pi * ( phi - 1.0 ) / phi;
  dr = 1.0;
%
%  Create a spiral in which the radius R and angle A both
%  increase by a constant increment,
%
  for i = 1 : n
    x(i) = r * cos ( a );
    y(i) = r * sin ( a );
    a = mod ( a + da, 2 * pi );
    r = r + dr;
  end
%
%  Display the data in a scatter plot.
%
  scatter ( x, y, 'b.' )
  axis equal
  title ( sprintf ( 'Fibonacci spiral, N = %d', n ) )

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function value = r8_cosd ( degrees )

%*****************************************************************************80
%
%% R8_COSD returns the cosine of an angle given in degrees.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real DEGREES, the angle in degrees.
%
%    Output, real VALUE, the cosine of the angle.
%
  radians = pi * ( degrees / 180.0 );

  value = cos ( radians );

  return
end


text:
function [distance_matrix, gram_matrix]= update_distance_matrix(distance_matrix, gram_matrix, gram_vector, new_sample_norm, id1, id2, w1, w2)
% Updates the distance matrix

% Normalise the weights so that they sum to one
alpha1 = w1/(w1+w2);
alpha2 = 1 - alpha1;

if id2 < 0
    norm_id1 = gram_matrix(id1, id1);

    % Update the gram matrix
    if alpha1 == 0
        % The new sample replaces an existing sample.
        gram_matrix(:,id1) = gram_vector;
        gram_matrix(id1,:) = gram_matrix(:,id1);
        gram_matrix(id1, id1) = new_sample_norm;
    elseif alpha2 == 0
        % The new sample is discared
    else
        % The new sample is merge with an existing sample
        gram_matrix(:,id1) = alpha1*gram_matrix(:,id1) + alpha2*gram_vector;
        gram_matrix(id1,:) = gram_matrix(:,id1);
        gram_matrix(id1, id1) = alpha1^2*norm_id1 + alpha2^2*new_sample_norm + 2*alpha1*alpha2*gram_vector(id1);
    end

    % Update distance matrix
    distance_matrix(:,id1) = max(gram_matrix(id1, id1) + diag(gram_matrix) - 2*gram_matrix(:,id1),0);
    distance_matrix(id1,:) = distance_matrix(:,id1) ;
    distance_matrix(id1,id1) = inf;
else
    if alpha1 == 0 || alpha2 == 0
        error('Error!');
    end

    % Two existing samples are merged and the new sample fills the empty
    % slot
    norm_id1 = gram_matrix(id1, id1);
    norm_id2 = gram_matrix(id2, id2);
    ip_id1_id2 = gram_matrix(id1,id2);

    % Handle the merge of existing samples
    gram_matrix(:,id1) = alpha1*gram_matrix(:,id1) + alpha2*gram_matrix(:,id2);
    gram_matrix(id1,:) = gram_matrix(:,id1);
    gram_matrix(id1, id1) = alpha1^2*norm_id1 + alpha2^2*norm_id2 + 2*alpha1*alpha2*ip_id1_id2;

    gram_vector(id1) = alpha1*gram_vector(id1) + alpha2*gram_vector(id2);

    % Handle the new sample
    gram_matrix(:,id2) = gram_vector;
    gram_matrix(id2,:) = gram_matrix(:,id2);
    gram_matrix(id2, id2) = new_sample_norm;

    % Update the distance matrix
    distance_matrix(:,id1) = max(gram_matrix(id1, id1) + diag(gram_matrix) - 2*gram_matrix(:,id1),0);
    distance_matrix(id1,:) = distance_matrix(:,id1) ;
    distance_matrix(id1,id1) = inf;

    distance_matrix(:,id2) = max(gram_matrix(id2, id2) + diag(gram_matrix) - 2*gram_matrix(:,id2),0);
    distance_matrix(id2,:) = distance_matrix(:,id2) ;
    distance_matrix(id2,id2) = inf;
end

author: martin-danelljan

repo: ECO


text:
function bench_ssc_mc_omp(md)

configs = cell(1, 4);
configs{1} = create_config([1 1 1 1 1 1 1 1 1 1], 4, '1-4');
configs{2} = create_config([2 1 1 1 1 1 1 1 1 1], 4, '2.1-4');
configs{3} = create_config([4 2 1 1 1 1 1 1 1 1], 4, '42.1-4');
configs{4} = create_config([2 2 2 2 2 2 2 2 2 2], 4, '2-4');
num_samples_per_digit_arr = [50 80 100 150 200 300 400];

nns = numel(num_samples_per_digit_arr);
nc = numel(configs);
for ns=1:nns
    num_samples_per_digit = num_samples_per_digit_arr(ns);
    for c=1:nc
        config = configs{c};
        solver_name = sprintf('ssc_mc_omp_%s', config.name);
        fprintf('Simulating: %s, for %d samples per digit\n', solver_name, num_samples_per_digit);
        result = check_digits(md, num_samples_per_digit, @ssc_mc_omp, solver_name, config);
        filepath = sprintf('bin/%s_digits=%d.mat', solver_name, num_samples_per_digit);
        save(filepath, 'result');
    end
end

end


text:
function [ALLEEG cfg] = pop_est_fitMVAR(ALLEEG,typeproc,varargin)
%
% Preprocess EEG dataset(s) for connectivity analysis. See [1] for
% mathematical details on preprocessing steps.
%
%
% Input:
%
%   ALLEEG:         Array of EEGLAB datasets to preprocess.
%   typeproc:       Reserved for future use. Use 0
%
% Optional:
%
%   <'Name',value> pairs as defined in pre_prepData()
%
% Output:
%
%   ALLEEG:         Prepocessed EEG structure(s)
%   cfg:            Argument specification structure.
%
%
% See Also: pre_prepData()
%
% References:
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual. Section 6.5.1
%   Available at: http://www.sccn.ucsd.edu/wiki/Sift
%
% Author: Tim Mullen 2009, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

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

if nargin<2
    typeproc = 0;
end

fcnName     = strrep(mfilename,'pop_','');
fcnHandle   = str2func(fcnName);

% check the dataset
res = hlp_checkeegset(ALLEEG,{'cat'});
if ~isempty(res)
    error(['SIFT:' fcnName],res{1});
end

if isfield(ALLEEG(1).CAT.configs,fcnName)
    % get default configuration (from prior use) and merge with varargin
    varargin = [hlp_struct2varargin(ALLEEG(1).CAT.configs.(fcnName)) varargin];
end

if strcmpi(typeproc,'nogui')
    % get the config from function
    cfg = arg_tovals(arg_report('rich',fcnHandle,[{'EEG',ALLEEG(1)},varargin]),false);
else
    % render the GUI
    [PGh figh] = feval(['gui_' fcnName],ALLEEG(1),varargin{:});

    if isempty(PGh)
        % user chose to cancel
        cfg = [];
        return;
    end

    % get the specification of the PropertyGrid
    ps = PGh.GetPropertySpecification;
    cfg = arg_tovals(ps,false);
end

drawnow;

if strcmpi(typeproc,'cfg_only')
    return;
end

% initialize progress bar
if cfg.verb==2 && length(ALLEEG)>1
    waitbarTitle = 'Fitting VAR Models';

    multiWaitbar(waitbarTitle,'Reset');
    multiWaitbar(waitbarTitle,'ResetCancel',true);
    multiWaitbar(waitbarTitle,...
                 'Color', [0.8 0.0 0.1],  ...
                 'CanCancel','on',        ...
                 'CancelFcn',@(a,b)disp('[Cancel requested. Please wait...]'));
end

% execute the low-level function
for cnd=1:length(ALLEEG)
    [ALLEEG(cnd).CAT.MODEL] = feval(fcnHandle,'EEG',ALLEEG(cnd),cfg);

    if ~isempty(cfg)
        % store the configuration structure
        ALLEEG(cnd).CAT.configs.(fcnName) = cfg;
    end

    if cfg.verb==2 && length(ALLEEG)>1
        % update waitbar
        drawnow;
        cancel = multiWaitbar(waitbarTitle,cnd/length(ALLEEG));
        if cancel && hlp_confirmWaitbarCancel(waitbarTitle)
            break;
        end
    end

end

% cleanup progress bar
if cfg.verb==2 && length(ALLEEG)>1
    multiWaitbar(waitbarTitle,'Close');
end

repo: REST


text:
function calpak_test375 ( )

%*****************************************************************************80
%
%% CALPAK_TEST375 tests MONTH_LENGTH_REPUBLICAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2013
%
%  Author:
%
%    John Burkardt
%
  n_test = 1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST375\n' );
  fprintf ( 1, '  For the Republican calendar:\n' );
  fprintf ( 1, '  MONTH_LENGTH_REPUBLICAN returns month lengths.\n' );

  y_test(1) = 4;

  for i_test = 1 : n_test

    y = y_test(i_test);
    sy = y_to_s_republican ( y );
    months = year_length_months_republican ( y );
    days = year_length_republican ( y );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %d\n', y );
    fprintf ( 1, '  %s\n', sy );
    fprintf ( 1, '  Year length in months = %d\n', months );
    fprintf ( 1, '  Year length in days = %d\n', days );
    fprintf ( 1, '\n' );

    for m = 1 : months
      month_name = month_to_month_name_republican ( m );
      fprintf ( 1, '  %10s  %2d\n', month_name, month_length_republican ( y, m ) );
    end

  end

  return
end


text:
function X=permute(X,p)
%PERMUTE (overloaded)

if length(X.dim) < length(p)
    X = ndsdpvar(X);
    X = permute(X,p);
    return
end
i = 1:prod(X.dim);
i = reshape(i,X.dim);
i = permute(i,p);
X.basis = X.basis(i,:);
X.dim = X.dim(p);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [po,freq] = spm_opm_psd(S)
% Compute PSD for OPM data(for checking noise floor)
% FORMAT D = spm_opm_psd(S)
%   S               - input structure
%  fields of S:
%   S.D             - SPM MEEG object                       - Default: no Default
%   S.triallength   - window size (ms)                      - Default: 1000
%   S.bc            - boolean to dc correct                 - Default: 0
%   S.channels      - channels to estimate PSD from         - Default: 'ALL'
%   S.plot          - boolean to plot or not                - Default: 0
%   S.units         - units of measurement                  - Default: 'fT'
%   S.constant      - constant line to draw as reference    - Default: 15
% Output:
%   psd             - power spectral density
%   f               - frequencies psd is sampled at
%__________________________________________________________________________
% Copyright (C) 2018 Wellcome Trust Centre for Neuroimaging

% Tim Tierney
% $Id: spm_opm_psd.m 7647 2019-07-25 13:58:53Z tim $

%-ArgCheck
%--------------------------------------------------------------------------

if ~isfield(S, 'units'),         S.units = 'fT'; end
if ~isfield(S, 'triallength'),   S.triallength = 1000; end
if ~isfield(S, 'constant'),      S.constant = 15; end
if ~isfield(S, 'bc'),            S.bc = 0; end
if ~isfield(S, 'channels'),      S.channels = 'ALL'; end
if ~isfield(S, 'plot'),          S.plot = 0; end
if ~isfield(S, 'D'),             error('D is required'); end

%-channel Selection
%--------------------------------------------------------------------------

    labs = S.channels;
    regex = {};

    for i = 1:length(labs)
        if isa(labs,'cell')
        regex{i} = ['regexp_(',labs{i},')'];
        else
            regex{i} = ['regexp_(',labs,')'];
        end
    end
    chans = [S.D.selectchannels(regex), indchantype(S.D,labs)];
%-epoch dataset
%--------------------------------------------------------------------------
if size(S.D,3)>1
    eD=S.D;
else
    nsamps = S.triallength/1000*S.D.fsample;
    beg = 1:nsamps:size(S.D,2);
    endsamp =  beg+(nsamps-1);
    inRange = ~(beg>size(S.D,2)|endsamp>size(S.D,2));
    eD = zeros(length(chans),nsamps,sum(inRange));
    for i =1:length(inRange)
        if(inRange(i))
            eD(:,:,i)=S.D(chans,beg(i):endsamp(i),:);
        end
    end

end
%- set window
%--------------------------------------------------------------------------
fs = S.D.fsample();
N =size(eD,2);
Nf= ceil((N+1)/2);
nepochs=size(eD,3);
pow = zeros(Nf,size(eD,1),nepochs);
wind  = window(@flattopwin,size(eD,2));
wind = repmat(wind,1,size(eD,1));

%- create PSD
%--------------------------------------------------------------------------

for j = 1:nepochs
    Btemp=eD(:,:,j)';
    Btemp = Btemp.*wind;
    mu=mean(Btemp);
    zf = bsxfun(@minus,Btemp,mu);
    if(S.bc)
        fzf = zf;
    else
        fzf=Btemp;
    end

    N= length(fzf);
    xdft = fft(fzf);
    xdft=xdft(1:floor(N/2+1),:);
    psdx = abs(xdft)./sqrt(N*fs);
    freq = 0:fs/size(fzf,1):fs/2;
    odd=mod(size(fzf,1),2)==1;
    if(odd)
        psdx(2:end) = sqrt(2)*psdx(2:end);
    else
        psdx(2:end-1) = sqrt(2)*psdx(2:end-1);
    end
    pow(:,:,j) =psdx;
end

%- plot
%--------------------------------------------------------------------------
po = median(pow(:,:,:),3);
if(S.plot)
    figure()
    semilogy(freq,po,'LineWidth',2);
    hold on
    xp2 =0:round(freq(end));
    yp2=ones(1,round(freq(end))+1)*S.constant;
    p2 =plot(xp2,yp2,'--k');
    p2.LineWidth=2;
    xlabel('Frequency (Hz)')
    labY = ['$$PSD (' S.units ' \sqrt[-1]{Hz}$$)'];
    ylabel(labY,'interpreter','latex')
    grid on
    ax = gca; % current axes
    ax.FontSize = 13;
    ax.TickLength = [0.02 0.02];
    fig= gcf;
    fig.Color=[1,1,1];

end

end

author: spm


text:
function [W] = amuse(X)
% BSS using eigenvalue value decomposition
% Program written by A. Cichocki and R. Szupiluk
%
% X [m x N] matrix of observed (measured) signals,
% W separating matrix,
% y estimated separated sources
% p time delay used in computation of covariance matrices
% optimal time-delay default p= 1
%
% First stage: Standard prewhitening

[m,N]=size(X);
if nargin==1,
 n=m; %
 end;

Rxx=(X*X')/N;

[Ux,Dx,Vx]=svd(Rxx);
 Dx=diag(Dx);
% n=xxx;
 if n<m, % under assumption of additive white noise and
        %when the number of sources are known or can a priori estimated
  Dx=Dx-real((mean(Dx(n+1:m))));
  Q= diag(real(sqrt(1./Dx(1:n))))*Ux(:,1:n)';
%
else    % under assumption of no additive noise and when the
        % number of sources is unknown
   n=max(find(Dx>1e-199)); %Detection the number of sources
   Q= diag(real(sqrt(1./Dx(1:n))))*Ux(:,1:n)';
end;
%
% else    %assumes no noise
%    Q=inv(sqrtm(Rxx));
% end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Second stage: Fast separation using sorting EVD
% notation the same as used in the Chapter 4
Xb=Q*X;
p=1;
% paramter p can take here value different than 1
% for example -1 or 2.
N=max(size(Xb));
Xb=Xb-kron(mean(Xb')',ones(1,N));

Rxbxbp=(Xb(:,1:N-1)*Xb(:,2:N)')/(N-1);
Rxbxbp= Rxbxbp+Rxbxbp';
[Vxb Dxb]=eig(Rxbxbp);
[D1 perm]=sort(diag(Dxb));
D1=flipud(D1);
Vxb=Vxb(:,flipud(perm));
W = Vxb'*Q;
%y = Vxb' * x1;

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function x = chi_cdf_inv ( cdf, a, b, c )

%*****************************************************************************80
%
%% CHI_CDF_INV inverts the Chi CDF.
%
%  Discussion:
%
%    A simple bisection method is used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%
%    Input, real A, B, C, the parameters of the PDF.
%    0 < B,
%    0 < C.
%
%    Output, real X, the corresponding argument of the CDF.
%
  it_max = 100;
  tol = 0.0001;

  if ( cdf <= 0.0 )
    x = a;
    return
  elseif ( 1.0 <= cdf )
    x = r8_huge ( );
    return
  end

  x1 = a;
  cdf1 = 0.0;

  x2 = a + 1.0;

  while ( 1 )

    cdf2 = chi_cdf ( x2, a, b, c );

    if ( cdf < cdf2 )
      break;
    end

    x2 = a + 2.0 * ( x2 - a );

  end
%
%  Now use bisection.
%
  it = 0;

  while ( 1 )

    it = it + 1;

    x3 = 0.5 * ( x1 + x2 );
    cdf3 = chi_cdf ( x3, a, b, c );

    if ( abs ( cdf3 - cdf ) < tol )
      x = x3;
      return
    end

    if ( it_max < it )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CHI_CDF_INV - Fatal error!\n' );
      fprintf ( 1, '  Iteration limit exceeded.\n' );
      error ( 'CHI_CDF_INV - Fatal error!' );
    end

    if ( ( cdf3 < cdf & cdf1 < cdf ) | ( cdf < cdf3 & cdf < cdf1 ) )
      x1 = x3;
      cdf1 = cdf3;
    else
      x2 = x3;
      cdf2 = cdf3;
    end

  end

  return
end


text:
function X = lleEmbed(Y, dims, neighbours)

% LLEEMBED Embed data set with LLE.

% MLTOOLS

% Wrapper for Sam Roweis' LLE code.

% Note LLE code uses the transpose of a design matrix.
if nargin < 3
  neighbours = 7;
end
if any(any(isnan(Y)))
  error('Cannot initialise gplvm using LLE when missing data is present.')
else
  X = lle(Y', neighbours, dims);
  X = X';
  % Rescale X so that variance is 1 and mean is zero.
  meanX = mean(X);
  X = X-ones(size(Y, 1), 1)*meanX;
  varX = var(X);
  X = X*diag(sqrt(1./varX));
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [medRV,minRV,medRVSS,minRVSS]=realized_min_med_variance(price,time,timeType,samplingType,samplingInterval,subsamples)
% Estimates truncated realized variance using the median and minimum method of Andersen, Dobrev and Schaumburg
%
% USAGE:
%   [MEDRV,MINRV] = realized_medrv(PRICE,TIME,TIMETYPE,SAMPLINGTYPE,SAMPLINGINTERVAL,SAMPLESPERBIN)
%   [MEDRV,MINRV,MEDRVSS,MINRVSS] = realized_medrv(PRICE,TIME,TIMETYPE,SAMPLINGTYPE,SAMPLINGINTERVAL,SAMPLESPERBIN,SUBSAMPLES)
%
% INPUTS:
%   PRICE            - m by 1 vector of high frequency prices
%   TIME             - m by 1 vector of times where TIME(i) corresponds to PRICE(i)
%   TIMETYPE         - String describing the way times are measured
%                       'wall'    24-hour clock of the form HHMMSS, e.g. 101543 or 153217
%                       'seconds' Time measured in seconds past midnight on the first day.
%                       'unit'  Unit normalized date format, e.g. .1, .234, .9
%                         Unit normalized times are more general than the other types and can be
%                         applied to data from more than one calendar day
%   SAMPLINGTYPE     - String describing the type of sampling to use when
%                        filtering PRICE
%                        'CalendarTime' - Sample in calendar time using observations separated by
%                          SAMPLINGINTERVAL seconds
%                        'CalendarUniform' - Sample in calendar time using SAMPLINGINTERVAL
%                          observations spread uniformly between TIME(1) and TIME(m)
%                        'BusinessTime' - Sample in business (tick) time using observation separated
%                          by SAMPLINGINTERVAL ticks
%                        'BusinessUniform' - Sample in business (tick) time using observations
%                          uniformly spaced in business time.
%                        'Fixed' - Sample at specific points in time. When using fixed,
%                          SAMPLINGINTERVAL must be a n by 1 vector of times with the same TIMETYPE
%                          as TIME (i.e. seconds if TIME is in seconds)
%   SAMPLINGINTERVAL  - Scalar integer or n by 1 vector whose meaning depends on SAMPLINGTYPE
%   SUBSAMPLES        - [OPTIONAL] Scalar integer indicating the number of subsample realized
%                         variance estimators to average with the original realized variance.
%                         Subsample realized variances are based on prices uniformly spaced between
%                         the times (Calendar sampling) or ticks (Business sampling).  SUBSAMPLES should
%                         be greater than 1, since 1 uses a single subsample and so is equivalent
%                         to the original RV will compute a subsample realized variance using the
%                         mid-point of the price sample points, 2 will use the original sample and
%                         the midpoint between the original sampling times, and so on.
%                         In general this number should be small so the subsample estimators will be "sparse".
%
% OUTPUTS:
%   MEDRV             - Median RV estimate
%   MINRV             - Minimum RV estimate
%   MEDRVSS           - Subsampled version of median RV estimate
%   MINRVSS           - Subsampled version of minimum RV estimate
%
% COMMENTS:
%
% EXAMPLE:
%  % Using all 1-minute returns between 9:30:00 and 16:00:00
%  fixedInterval = seconds2wall(wall2seconds(93000):60:wall2seconds(160000));
%  [MEDRV,MINRV] = realized_medrv(PRICE,TIME,'wall','Fixed',fixedInterval)
%
%  % Using all 5-minute returns between 9:30:00 and 16:00:00, but subsampling every minute
%  fixedInterval = seconds2wall(wall2seconds(93000):300:wall2seconds(160000));
%  [MEDRV,MINRV,MEDRVSS,MINRVSS] = realized_medrv(PRICE,TIME,'wall','Fixed',fixedInterval,5)
%
%  See also REALIZED_BIPOWER_VARIATION, REALIZED_KERNEL, REALIZED_QUANTILE_VARIANCE,
%  REALIZED_RANGE, REALIZED_PRICE_FILTER

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 10/24/2011

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<5 || nargin>6
    error('Five or Six inputs required.')
end
if size(price,2)>size(price,1)
    price=price';
end
if size(price,2)>1
    error('PRICE must be a m by 1 vector.')
end
if size(time,2)>size(time,1)
    time=time';
end
if any(diff(time)<0)
    error('TIME must be sorted and increasing')
end
if size(time,2)>1 || length(time)~=length(price)
    error('TIME must be a m by 1 vector.')
end
% Inserted to protect against inputing integer times
time = double(time);

timeType=lower(timeType);
if ~ismember(timeType,{'wall','seconds','unit'})
    error('TIMETYPE must be one of ''wall'', ''seconds'' or ''unit''.');
end
samplingType=lower(samplingType);
if ~ismember(samplingType,{'calendartime','calendaruniform','businesstime','businessuniform','fixed'})
    error('SAMPLINGTYPE must be one of ''CalendarTime'', ''CalendarUniform'', ''BusinessTime'', ''BusinessUniform'' or ''Fixed''.');
end

m=size(price,1);
t0=time(1);
tT=time(m);
if ismember(samplingType,{'calendartime','calendaruniform','businesstime','businessuniform'})
    % Must be a scalar integer if timeType is seconds or wall
    if ismember(timeType,{'wall','seconds'})
        if ~isscalar(samplingInterval) || floor(samplingInterval)~=samplingInterval || samplingInterval<1
            error('SAMPLINGINTERVAL must be a positive integer for the SAMPLINGTYPE selected when using ''wall'' or ''seconds'' as TIMETYPE.')
        end
    else
        if ~isscalar(samplingInterval) || samplingInterval<0
            error('SAMPLINGINTERVAL must be a positive value for the SAMPLINGTYPE selected when using ''unit'' as TIMETYPE.')
        end
    end
else
    if size(samplingInterval,2)>size(samplingInterval,1)
        samplingInterval=samplingInterval';
    end
    if ~(any(samplingInterval>=t0) && any(samplingInterval<=tT))
        error('At least one sampling interval must be between min(TIME) and max(TIME) when using ''Fixed'' as SAMPLINGTYPE.')
    end
    if any(diff(samplingInterval)<=0)
        error('When using ''Fixed'' as SAMPLINGTYPE the vector of sampling times in SAMPLINGINTERVAL must be sorted and strictly increasing.')
    end
end

if nargin<6 || isempty(subsamples)
    subsamples = 1;
end

if nargin==6 && ~isempty(subsamples)
    if ~isscalar(subsamples) || subsamples<0 || floor(subsamples)~=subsamples
        error('SUBSAMPLES must be a non-negative scalar.')
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

logPrice = log(price);
% Convert times to unit if not already unit
if ~strcmp(timeType,'unit')
    [time,~,~,samplingInterval]=realized_convert2unit(time,timeType,samplingType,samplingInterval);
end
% Since everything is converted to unit, set timeType to unit
timeType='unit';
% m is the length of price

% Filter prices and compute the RV
filteredLogPrice = realized_price_filter(logPrice,time,timeType,samplingType,samplingInterval);
% Construct the aggregated return series.  Remember will need to skip
returns = diff(filteredLogPrice);
m = length(returns);

returns2 = (sqrt(m)*returns).^2;
r2 = [returns2(1:m-2) returns2(2:m-1) returns2(3:m)];
medRV = mean(median(r2,2));
medRVScale = pi/(6-4*sqrt(3)+pi);
medRV = medRV*medRVScale ;

r2 = [returns2(1:m-1) returns2(2:m)];
minRV = mean(min(r2,[],2));
minRVscale = pi/(pi-2);
minRV = minRV*minRVscale;

subsampledLogPrices = realized_subsample(logPrice,time,timeType,samplingType,samplingInterval,subsamples);
medRVs = nan(m*subsamples,1);
minRVs = nan(m*subsamples,1);
medRVcount = 0;
minRVcount = 0;
for i=1:subsamples
    filteredLogPrice = subsampledLogPrices{i};
    returns = diff(filteredLogPrice);
    m = size(returns,1);
    returns2 = m*returns.^2;
    % MedRV
    r2 = [returns2(1:m-2) returns2(2:m-1) returns2(3:m)];
    n = size(r2,1);
    medRVs(medRVcount+(1:n)) = median(r2,2);
    medRVcount = medRVcount  + n;
    % MinRV
    r2 = [returns2(1:m-1) returns2(2:m)];
    n = size(r2,1);
    minRVs(minRVcount+(1:n)) = min(r2,[],2);
    minRVcount = minRVcount + n;
end
medRVSS = mean(medRVs(1:medRVcount));
medRVSS = medRVSS * medRVScale;
minRVSS = mean(minRVs(1:minRVcount));
minRVSS = minRVSS * minRVscale;

author: bashtage

repo: mfe-toolbox


text:
% Test file for SPINSPHERE:

function pass = test_spinsphere()

tol = 1e-2;

%% AC:

% Solve with DT and DT/2:
S = spinopsphere('AC'); S.tspan = S.tspan/10;
N = 128; dt = 1e-1;
u = spinsphere(S, N, dt, 'plot', 'off');
v = spinsphere(S, N, dt/2, 'plot', 'off');

% Compare:
dom = S.domain;
[xx, yy] = meshgrid(linspace(dom(1), dom(2), 50));
scale = max(max(abs(v(xx,yy))));
pass(1) = max(max(abs(u(xx,yy) - v(xx,yy))))/scale < tol;

%% GL:

% Solve with DT and DT/2:
S = spinopsphere('GL'); S.tspan = S.tspan/10;
S.init =  spherefun(@(x,y,z) cos(cosh(x.*z)-y));
N = 128; dt = 1e-1;
u = spinsphere(S, N, dt, 'plot', 'off');
v = spinsphere(S, N, dt/2, 'plot', 'off');

% Compare:
dom = S.domain;
[xx, yy] = meshgrid(linspace(dom(1), dom(2), 50));
scale = max(max(abs(v(xx,yy))));
pass(2) = max(max(abs(u(xx,yy) - v(xx,yy))))/scale < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function Problem = ITERmake_CON_MAT(Problem, Constraints, Objectives, Spec)

disp(' ')
nxny = input('Number of instances of swept variables [7]: ');
if isempty(nxny); nxny = 7; end
nxny = round(nxny);
if nxny < 3; return; end
disp(' ')

tic;
nPPreqs = Spec.nPPrequirements;
Constraints.current_active_cons = Constraints.active_cons(1:nPPreqs);
if length(Constraints.active_cons)>nPPreqs
    Constraints.current_active_nonlcons = Constraints.active_cons(nPPreqs+1:end);
    Constraints.activeNonLLabels = Constraints.conLabels(nPPreqs+1:end);
    Constraints.activeNonLLabels = Constraints.activeNonLLabels(Constraints.current_active_nonlcons);
else
    Constraints.current_active_nonlcons = [];
    Constraints.activeNonLLabels = [];
end

Constraints.activePPLabels = Constraints.conLabels(1:nPPreqs);
Constraints.activePPLabels = Constraints.activePPLabels(Constraints.current_active_cons);

Constraints.nactiveMar = nnz(Constraints.current_active_cons);
Constraints.nactivenonL = nnz(Constraints.current_active_nonlcons);

Constraints.colors = [0 0 1
    0 1 0
    1 0 0
    0 1 1
    1 0 1];
while length(Constraints.colors)<Constraints.nactiveMar+Constraints.nactivenonL
    Constraints.colors = [Constraints.colors;Constraints.colors];
end
Constraints.colors = Constraints.colors(1:Constraints.nactiveMar+Constraints.nactivenonL);

activeX = Problem.activeX;
x0 = Problem.x0(activeX); %the x0 of active values;
nvars = nnz(activeX);
lb = Problem.lb(activeX);
ub = Problem.ub(activeX);
XLabels=Problem.XLabels(activeX);

if nvars < 2;
    disp('Visualizer for 2 or more decision variables')
    beep; return
end

% first get values for baseline aircraft:
disp ('Evaluating baseline aircraft...')
[Constraints.C0,~,Constraints.MarStack0,...
    Constraints.NonLMarStack0] = EVALnonlcons(x0,Problem,Constraints,Spec);
[Objectives.ObjectiveValue0, Objectives.ObjVals0] = EVALobjective...
    (x0,Problem,Objectives,Spec,Constraints,'single');
nsubplots = nvars*(nvars-1)/2;
ticker = 0;
progressbar('Constraint matrix progress' , 'Subplot sweep')
Matrix(nvars,nvars).X = [];
for ii = 1:nvars;
    for jj = 1:nvars;
        if ii~=jj %make a plot from Matrix
            if jj > ii %make a carpet plot and store it
                action = ['Generating plot ' XLabels{ii} ' versus ' XLabels{jj} '...'];
                disp(action);
                [Matrix(ii,jj).X,Matrix(ii,jj).Y,Matrix(ii,jj).ObjStack,...
                    Matrix(ii,jj).MarStack,Matrix(ii,jj).feasible,...
                    Matrix(ii,jj).feasCompObj]=...
                    carpet_plot_vals(ii,jj,...
                    x0,lb,ub,Problem, Constraints, Objectives, Spec, nxny,...
                    ticker, nsubplots);
            ticker = ticker+1;
                progressbar(ticker/nsubplots);
            elseif jj < ii %just do a reflection
                Matrix(ii,jj).Y = Matrix(jj,ii).X;
                Matrix(ii,jj).X = Matrix(jj,ii).Y;
                Matrix(ii,jj).ObjStack = Matrix(jj,ii).ObjStack;
                Matrix(ii,jj).MarStack = Matrix(jj,ii).MarStack;
                Matrix(ii,jj).feasible = Matrix(jj,ii).feasible;
            end
        end
    end
end
progressbar(1);
save
disp('Workspace saved to matlab.mat')
Tstamp = num2str(fix(clock),'%02d');
fname = ['Autosaves\CarpetMatrixData' Tstamp '.mat'];
save(fname)
toc;
beep;pause(.5);beep;pause(.5);beep;pause(.5);
beep;pause(.1);beep;pause(.1);beep;
Haxes = CMATplot_the_carpet...
    (Matrix, Constraints,Objectives,x0,200,lb,ub,XLabels);
toc;
Problem = CMATpick_point(Problem,Constraints,Spec,Objectives,Haxes);
if get_yes_or_no('Save result to History.xls? [Y]/N: ',true)
    ITERsave_hist(Problem,Constraints,Objectives,Spec,'Visualizer');
end
end

function [XX,YY,ObjStack,MarsStack,feasible,feasCompObj]=...
    carpet_plot_vals(Xyind,Xxind,x0,lb,ub,...
    Problem, Constraints, Objectives, Spec, nxny,...
    ticker, nsubplots)
% nx = 9; ny = 8;
nx = nxny;
ny = nxny;

xx = linspace3(lb(Xxind),x0(Xxind),ub(Xxind),nx);
yy = linspace3(lb(Xyind),x0(Xyind),ub(Xyind),ny);

[YY,XX]=meshgrid(yy,xx);

ObjStack = NaN(nx,ny,1+Objectives.nObj);
MarsStack = NaN(nx,ny,Constraints.nactiveMar+Constraints.nactivenonL);

waitbar1tot = nx*ny-1;
waitbar1prog = 0;
progressbar([],0)
for ii = 1:length(xx)
    for jj = 1:length(yy)

        X = x0;
        X(Xxind) = xx(ii);
        X(Xyind) = yy(jj);
        if all(X==x0) %skip evaluating the seed point
            ObjStack(ii,jj,:) = ...
                [Objectives.ObjectiveValue0; Objectives.ObjVals0(:)];
            MarsStack(ii,jj,:) = -Constraints.C0(:);
        else
            C = EVALnonlcons(X,Problem,Constraints,Spec);
            MarsStack(ii,jj,:) = -C(:);
            [ObjectiveValue, ObjVals] = EVALobjective...
                (X,Problem,Objectives,Spec,Constraints,'single');
            ObjStack(ii,jj,:) = [ObjectiveValue; ObjVals(:)];
            waitbar1prog = waitbar1prog+1;
            frac2 = waitbar1prog/waitbar1tot;
            progressbar((ticker+frac2)/nsubplots,frac2)
        end
    end
end
feasible = all(MarsStack>=0,3);
feasCompObj = ObjStack(:,:,1);
feasCompObj(~feasible) = NaN;
end

author: Sable

repo: mcbench-benchmarks


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function result = xstar(ival, cp, a, b, iter, Grid_k, ...
    model, V, t, r, q, strike, varargin)

cfvals = exp(feval(@CF, model, pi * Grid_k * diag(1./ (b-a)), t,r,q,varargin{:}));

x = ival;
eps = 1e-6;

for n = 1:iter
    exp_t = exp( 1i * pi * Grid_k * diag((x - a) ./ (b - a)) );
    vec = real( cfvals .* exp_t ) .* V;
    vec(1,:) = 0.5*vec(1,:);

    g = (exp(-r * t) * sum(vec, 1)') ...
            - cp .* strike .* (exp(x) - 1);

    vec = imag(cfvals .* Grid_k .* exp_t) .* V;
    vec(1,:) = 0.5*vec(1,:);

    dg = - exp(-r * t) .* pi ./ (b - a) .* sum(vec, 1)' - cp .* strike .* exp(x);

    x = x - (g ./ dg);
    if abs(g) < eps
        break;
    end
end

result = x;

end


text:
function value = agud ( g )

%*****************************************************************************80
%
%% AGUD evaluates the inverse Gudermannian function.
%
%  Discussion:
%
%    The Gudermannian function relates the hyperbolic and trigonomentric
%    functions.  For any argument X, there is a corresponding value
%    G so that
%
%      SINH(X) = TAN(G).
%
%    This value G(X) is called the Gudermannian of X.  The inverse
%    Gudermannian function is given as input a value G and computes
%    the corresponding value X.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real G, the value of the Gudermannian.
%
%    Output, real VALUE, the argument of the Gudermannian.
%
  value = log ( tan ( 0.25 * pi + 0.5 * g ) );

  return
end


text:
function [q,qd,q2d] = fourier_series_traj(t,q0,A,B,w,N)
% -----------------------------------------------------------------------
% The function evaluates trancated Fourier series polynomial
% Inputs:
%   t - time instantes at which the polynomail should be evaluated
%   q0 - initial offset, should be vector (nx1)
%   A - sine coefficients
%   B - cosine coefficients
%   w - fundamental frequency
%   N - number of harmonics
% Outputs:
%   q - positions i.e. value of polynomial
%   qd - velocities i.e. derivative of polynomail
%   q2d - accelerations i.e. second derivative of polynomial
% -----------------------------------------------------------------------
q = q0;
qd = zeros(size(q0));
q2d = zeros(size(q0));
for k = 1:N
   q = q + A(:,k)/(w*k).*sin(w*k*t) - B(:,k)/(w*k).*cos(w*k*t);
   qd = qd + A(:,k).*cos(w*k*t) + B(:,k)*sin(w*k*t);
   q2d = q2d - A(:,k)*w*k.*sin(w*k*t) + B(:,k)*w*k*cos(w*k*t);
end

author: shamilmamedov


text:
function [Az El] = RaDec2AzEl(Ra,Dec,lat,lon,time)
% Programed by Darin C. Koblick 01/23/2010
%--------------------------------------------------------------------------
% External Function Call Sequence:
% [Az El] = RaDec2AzEl(0,0,0,-104,'1992/08/20 12:14:00')
%
% Worked Example: pg. 262 Vallado
%[Az El] = RaDec2AzEl(294.9891115,-20.8235624,39.007,-104.883,'1994/05/14 13:11:20.59856')
%[210.7514  23.9036] = RaDec2AzEl(294.9891115,-20.8235624,39.007,-104.883,'1994/05/14 13:11:20.59856')
%
% Worked Example: http://www.stargazing.net/kepler/altaz.html
% [Az El] = RaDec2AzEl(344.95,42.71667,52.5,-1.91667,'1997/03/14 19:00:00')
% [311.92258 22.40100] = RaDec2AzEl(344.95,42.71667,52.5,-1.91667,'1997/03/14 19:00:00')
%
% [Beta,el] = RaDec2AzEl(alpha_t,delta_t,phi,lamda,'yyyy/mm/dd hh:mm:ss')
%
% Function Description:
%--------------------------------------------------------------------------
% RaDec2AzEl will take the Right Ascension and Declination in the topocentric
% reference frame, site latitude and longitude as well as a time in GMT
% and output the Azimuth and Elevation in the local horizon
% reference frame.
%
% Inputs:                                                       Format:
%--------------------------------------------------------------------------
% Topocentric Right Ascension (Degrees)                         [N x 1]
% Topocentric Declination Angle (Degrees)                       [N x 1]
% Lat (Site Latitude in degrees -90:90 -> S(-) N(+))            [N x 1]
% Lon (Site Longitude in degrees -180:180 W(-) E(+))            [N x 1]
% UTC (Coordinated Universal Time YYYY/MM/DD hh:mm:ss)          [N x 1]
%
% Outputs:                                                      Format:
%--------------------------------------------------------------------------
% Local Azimuth Angle   (degrees)                               [N x 1]
% Local Elevation Angle (degrees)                               [N x 1]
%
%
% External Source References:
% Fundamentals of Astrodynamics and Applications
% D. Vallado, Second Edition
% Example 3-5. Finding Local Siderial Time (pg. 192)
% Algorithm 28: AzElToRaDec (pg. 259)
% -------------------------------------------------------------------------

%Example 3-5
[yyyy mm dd HH MM SS] = datevec(datenum(time,'yyyy/mm/dd HH:MM:SS'));
JD = juliandate(yyyy,mm,dd,HH,MM,SS);
T_UT1 = (JD-2451545)./36525;
ThetaGMST = 67310.54841 + (876600*3600 + 8640184.812866).*T_UT1 ...
+ .093104.*(T_UT1.^2) - (6.2*10^-6).*(T_UT1.^3);
ThetaGMST = mod((mod(ThetaGMST,86400*(ThetaGMST./abs(ThetaGMST)))/240),360);
ThetaLST = ThetaGMST + lon;

%Equation 4-11 (Define Siderial Time LHA)
LHA = mod(ThetaLST - Ra,360);

%Equation 4-12 (Elevation Deg)
El = asind(sind(lat).*sind(Dec)+cosd(lat).*cosd(Dec).*cosd(LHA));

%Equation 4-13 / 4-14 (Adaptation) (Azimuth Deg)
%Az = mod(atand(-(sind(LHA).*cosd(Dec)./(cosd(lat).*sind(Dec) - sind(lat).*cosd(Dec).*cosd(LHA)))),360);
Az = mod(atan2(-sind(LHA).*cosd(Dec)./cosd(El),...
    (sind(Dec)-sind(El).*sind(lat))./(cosd(El).*cosd(lat))).*(180/pi),360);

function jd = juliandate(year, month, day, hour, min, sec)
YearDur = 365.25;
for i = length(month):-1:1
    if (month(i)<=2)
        year(i)=year(i)-1;
        month(i)=month(i)+12;
    end
end
A = floor(YearDur*(year+4716));
B = floor(30.6001*(month+1));
C = 2;
D = floor(year/100);
E = floor(floor(year/100)*.25);
F = day-1524.5;
G = (hour+(min/60)+sec/3600)/24;
jd =A+B+C-D+E+F+G;


text:
function exactness_test08 ( )

%*****************************************************************************80
%
%% EXACTNESS_TEST08 tests Gauss-Chebyshev1 rules for the Chebyshev1 integral.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 May 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EXACTNESS_TEST08\n' );
  fprintf ( 1, '  Test Gauss-Chebyshev1 rules for the Chebyshev1 integral.\n' );
  fprintf ( 1, '  Density function rho(x) = 1/sqrt(1-x^2).\n' );
  fprintf ( 1, '  Region: -1 <= x <= +1.\n' );
  fprintf ( 1, '  Exactness: 2*N-1.\n' );

  for n = 1 : 5

    [ x, w ] = chebyshev1_set ( n );
    p_max = 2 * n;
    chebyshev1_exactness ( n, x, w, p_max );

  end

  return
end


text:
function test_failed=test_rangecompress
%TEST_RANGECOMPRESS Test range compression and expansion

test_failed=0;

disp(' ===============  TEST_RANGECOMPRESS ================');

x=tester_crand(5,11);

y=rangecompress(x,'mulaw');
x_r=rangeexpand(y,'mulaw');

res=norm(x-x_r,'fro');

[test_failed,fail]=ltfatdiditfail(res,test_failed);
fprintf(['RANGECOMPRESS MULAW %0.5g %s\n'],res,fail);

y=rangecompress(x,'alaw');
x_r=rangeexpand(y,'alaw');

res=norm(x-x_r,'fro');

[test_failed,fail]=ltfatdiditfail(res,test_failed);
fprintf(['RANGECOMPRESS  ALAW %0.5g %s\n'],res,fail);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [out1,out2] = spline_scale_function(varargin)
out2 = linspace(varargin{1:3});        % wavelet support.
out1 = spline_function(out2);