text:
classdef matRad_MinMaxDose < DoseConstraints.matRad_DoseConstraint
    % matRad_MinMaxDose Implements a MinMaxDose constraint
    %   See matRad_DoseConstraint for interface description
    %
    % use log sum exp approximation, see appendix A in
    % http://scitation.aip.org/content/aapm/journal/medphys/41/8/10.1118/1.4883837
    %
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    % Copyright 2020 the matRad development team.
    %
    % This file is part of the matRad project. It is subject to the license
    % terms in the LICENSE file found in the top-level directory of this
    % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
    % of the matRad project, including this file, may be copied, modified,
    % propagated, or distributed except according to the terms contained in the
    % LICENSE file.
    %
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    properties (Constant)
        name = 'Min/Max dose constraint';
        parameterNames = {'d^{min}', 'd^{max}','method'};
        parameterTypes = {'dose','dose',{'approx','voxelwise'}};
    end

    properties
        parameters = {0,30,1};
        epsilon = 1e-3; %slack parameter for the logistic approximation
    end

    methods
        function obj = matRad_MinMaxDose(minDose,maxDose,method)

            %If we have a struct in first argument
            if nargin == 1 && isstruct(minDose)
                inputStruct = minDose;
                initFromStruct = true;
            else
                initFromStruct = false;
                inputStruct = [];
            end

            %Call Superclass Constructor (for struct initialization)
            obj@DoseConstraints.matRad_DoseConstraint(inputStruct);

            %now handle initialization from other parameters
            if ~initFromStruct

                if nargin < 3 || ~ischar(method)
                    method = 'approx';
                end

                methodIx = find(strcmp(method,obj.parameterTypes{3}));

                if isempty(methodIx) || numel(methodIx) > 1
                    methodIx = 1;
                    msg = ['Dose Constraint method can only be ', strjoin(obj.parameterTypes{3},' or '), '! Using method ''', obj.parameterTypes{3}{methodIx}, '''.'];
                    warning(msg);
                end

                obj.parameters{3} = methodIx;

                if nargin >= 2 && isscalar(maxDose)
                    obj.parameters{2} = maxDose;
                end

                if nargin >= 1 && isscalar(minDose)
                    obj.parameters{1} = minDose;
                end
            end
        end

        %Overloads the struct function to add constraint specific
        %parameters
        function s = struct(obj)
            s = struct@DoseConstraints.matRad_DoseConstraint(obj);
            s.epsilon = obj.epsilon;
        end

        function cu = upperBounds(obj,n)
            switch obj.parameters{3}
                case 1 %logsumexp approx
                    %Validate parameters
                    if obj.parameters{1} <= 0 && isinf(obj.parameters{2}) %Constraint doesn't make sense (min = 0 & max = Inf)
                        cu = [];
                    elseif obj.parameters{2} == Inf %Only min dose
                        cu = Inf;
                    elseif obj.parameters{1} <= 0 %Only max dose
                        cu = obj.parameters{2};
                    else %both are set sensible
                        cu = [Inf; obj.parameters{2}];
                    end
                case 2 %voxelwise
                    cu = obj.parameters{2}*ones(n,1);
                otherwise
                    error(['Min/max dose constraint evaluation method not known!']);
            end
            %cu = [Inf; obj.parameters{2}];
        end
        function cl = lowerBounds(obj,n)
            switch obj.parameters{3}
                case 1 %logsumexp approx
                    if obj.parameters{1} <= 0 && isinf(obj.parameters{2})
                        cl = [];
                    elseif obj.parameters{2} == Inf
                        cl = obj.parameters{1};
                    elseif obj.parameters{1} <= 0
                        cl = 0;
                    else
                        cl = [obj.parameters{1}; 0];
                    end
                case 2
                    cl = obj.parameters{1}*ones(n,1);
                otherwise
                    matRad_cfg = MatRad_Config.instance();
                    matRad_cfg.dispError('Min/max dose constraint evaluation method not known!');
            end
        end

        function jStruct = getDoseConstraintJacobianStructure(obj,n)
            switch obj.parameters{3}
                case 1 %logsumexp approx
                    %Validate parameters
                    if obj.parameters{1} <= 0 && isinf(obj.parameters{2}) %Constraint doesn't make sense (min = 0 & max = Inf)
                        jStruct = ones(n,0);
                    elseif obj.parameters{1} > 0 && isfinite(obj.parameters{2}) %both are set sensible
                        jStruct = ones(n,2);
                    else %Only min or max dose
                        jStruct = ones(n,1);
                    end
                    %jStruct = ones(n,2);
                case 2
                    jStruct = speye(n);
                otherwise
                    matRad_cfg = MatRad_Config.instance();
                    matRad_cfg.dispError('Min/max dose constraint evaluation method not known!');
            end

        end

        %% Calculates the Constraint Function value
        function cDose = computeDoseConstraintFunction(obj,dose)
            %cDose(2) = dose_max + obj.epsilon * log( sum(exp((dose - dose_max)/obj.epsilon)));
            %cDose(1) = dose_min - obj.epsilon * log( sum(exp((dose_min - dose)/obj.epsilon)));
            switch obj.parameters{3}
                case 1 %logsumexp approx
                    cDose = obj.computeDoseConstraintFunctionLogSumExp(dose);
                case 2
                    cDose = obj.computeDoseConstraintFunctionVoxelwise(dose);
                otherwise
                    matRad_cfg = MatRad_Config.instance();
                    matRad_cfg.dispError('Min/max dose constraint evaluation method not known!');
            end
        end

        %% Calculates the Constraint jacobian
        function cDoseJacob  = computeDoseConstraintJacobian(obj,dose)
            switch obj.parameters{3}
                case 1 %logsumexp approx
                    cDoseJacob = obj.computeDoseConstraintJacobianLogSumExp(dose);
                case 2
                    cDoseJacob = obj.computeDoseConstraintJacobianVoxelwise(dose);
                otherwise
                    matRad_cfg = MatRad_Config.instance();
                    matRad_cfg.dispError('Min/max dose constraint evaluation method not known!');
            end
        end
    end

    methods (Access = private)
        % LogSumExp Approximation
        function cDose = computeDoseConstraintFunctionLogSumExp(obj,dose)
            dose_min = min(dose);
            dose_max = max(dose);

            %Validate parameters
            if obj.parameters{1} <= 0 && isinf(obj.parameters{2}) %Constraint doesn't make sense (min = 0 & max = Inf)
                cDose = [];
            elseif obj.parameters{2} == Inf %Only min dose
                cDose = dose_min - obj.epsilon * log( sum(exp((dose_min - dose)/obj.epsilon)));
            elseif obj.parameters{1} <= 0 %Only max dose
                cDose = dose_max + obj.epsilon * log( sum(exp((dose - dose_max)/obj.epsilon)));
            else %both are set sensible
                cDose(2,1) = dose_max + obj.epsilon * log( sum(exp((dose - dose_max)/obj.epsilon)));
                cDose(1,1) = dose_min - obj.epsilon * log( sum(exp((dose_min - dose)/obj.epsilon)));
            end

        end
        function cDoseJacob  = computeDoseConstraintJacobianLogSumExp(obj,dose)
            %Validate parameters
            if obj.parameters{1} <= 0 && isinf(obj.parameters{2}) %Constraint doesn't make sense (min = 0 & max = Inf)
                cDoseJacob = [];
            elseif obj.parameters{2} == Inf %Only min dose
                cDoseJacob(:,1) = exp( (min(dose)-dose)/obj.epsilon );
                cDoseJacob(:,1) = cDoseJacob(:,1)/sum(cDoseJacob(:,1));
            elseif obj.parameters{1} <= 0 %Only max dose
                cDoseJacob(:,1) = exp( (dose-max(dose))/obj.epsilon );
                cDoseJacob(:,1) = cDoseJacob(:,1)/sum(cDoseJacob(:,1));
            else %both are set sensible
                cDoseJacob(:,1) = exp( (min(dose)-dose)/obj.epsilon );
                cDoseJacob(:,1) = cDoseJacob(:,1)/sum(cDoseJacob(:,1));

                cDoseJacob(:,2) = exp( (dose-max(dose))/obj.epsilon );
                cDoseJacob(:,2) = cDoseJacob(:,2)/sum(cDoseJacob(:,2));
            end

        end

        %Exact voxel-wise
        function cDose = computeDoseConstraintFunctionVoxelwise(obj,dose)
            cDose = dose;
        end
        function cDoseJacob  = computeDoseConstraintJacobianVoxelwise(obj,dose)
            cDoseJacob = speye(numel(dose),numel(dose));
        end
    end

end

repo: matRad


text:
function [ vX, mX ] = SolveLsL1ComplexAdmm( mA, vB, paramLambda, numIterations )
% ----------------------------------------------------------------------------------------------- %
%[ vX, mX ] = SolveLsL1ComplexSubGrad( mA, vB, lambdaFctr, numIterations )
% Solves the 0.5 * || A x - b ||_2 + \lambda || x ||_1 problem using ADMM
% Method. The model allows A, b and x to be Complex.
% Input:
%   - mA                -   Model Matrix.
%                           The model matrix.
%                           Structure: Matrix (m X n).
%                           Type: 'Single' / 'Double' (Complex).
%                           Range: (-inf, inf).
%   - vB                -   Input Vector.
%                           The model known data.
%                           Structure: Vector (m X 1).
%                           Type: 'Single' / 'Double' (Complex).
%                           Range: (-inf, inf).
%   - paramLambda       -   Parameter Lambda.
%                           The L1 Regularization parameter.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: (0, inf).
%   - numIterations     -   Number of Iterations.
%                           Number of iterations of the algorithm.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range {1, 2, ...}.
% Output:
%   - vX                -   Output Vector.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References
%   1.  Wikipedia ADMM Method - https://en.wikipedia.org/wiki/Augmented_Lagrangian_method#Alternating_direction_method_of_multipliers.
% Remarks:
%   1.  This implementation assumes the Rho Factor is 1.
% Known Issues:
%   1.  A
% TODO:
%   1.  B
% Release Notes:
%   -   1.0.000     07/11/2016
%       *   First realease version.
% ----------------------------------------------------------------------------------------------- %

numRows = size(mA, 1);
numCols = size(mA, 2);

vX  = pinv(mA) * vB; %<! Dealing with "Fat Matrix"
vU = zeros([numCols, 1]);
vV = zeros([numCols, 1]);

mAA = mA' * mA;
vAb = mA' * vB;
mI = eye(numCols);

mAAI = mAA + mI;

mL = chol(mAAI, 'lower');
mU = mL';

sL.LT = true();
sU.UT = true();

mX = zeros([numCols, numIterations]);
mX(:, 1) = vX;

for ii = 2:numIterations
    % vX = mAAI \ (vAb + vV - vU);
    vX = linsolve(mU, linsolve(mL, vAb + vV - vU, sL), sU); %<! https://www.mathworks.com/help/matlab/ref/linsolve.html
    vV = ProxL1(vX + vU, paramLambda);
    vU = vU + vX - vV;

    mX(:, ii) = vX;
end

end

function [ vX ] = ProxL1( vX, lambdaFactor )

% Soft Thresholding - Complex Domain -> Keep Phase, Soft Threshold the
% Modulus

% vXAbs   = abs(vX);
% vXPhase = angle(vX);
%
% vX = max(vXAbs - lambdaFactor, 0) .* exp(1i * vXPhase);

vXAbs = abs(vX);

vX              = (vX ./ vXAbs) .* max((vXAbs - lambdaFactor), 0);
vX(vXAbs == 0)  = 0;

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function model=addRxnsGenesMets(model,sourceModel,rxns,addGene,rxnNote,confidence)
% addRxnsGenesMets
%   Copies reactions from a source model to a new model, including
%   (new) metabolites and genes
%
%   model           draft model where reactions should be copied to
%   sourceModel     model where reactions and metabolites are sourced from
%   rxns            cell array with reaction IDs (from source model). Can also
%                   be string if only one reaction is added
%   addGene         three options:
%                   false   no genes are annotated to the new reactions
%                   true    grRules ared copied from the sourceModel and
%                           new genes are added when required
%                   string or cell array
%                           new grRules are specified as string or cell
%                           array, and any new genes are added when
%                           required
%                   (opt, default false)
%   rxnNote         cell array with strings explaining why reactions were copied
%                   to the model, to be included as newModel.rxnNotes. Can also
%                   be string if same rxnNotes should be added for each new
%                   reaction, or only one reaction is to be added (opt, default
%                   'Added via addRxnsAndMets()')
%   confidence      integer specifying confidence score for all reactions.
%                   4:  biochemical data: direct evidence from enzymes
%                       assays
%                   3:  genetic data: knockout/-in or overexpression
%                       analysis
%                   2:  physiological data: indirect evidence, e.g.
%                       secretion products or defined medium requirement
%                       sequence data: genome annotation
%                   1:  modeling data: required for functional model,
%                       hypothetical reaction
%                   0:  no evidence
%                   following doi:10.1038/nprot.2009.203 (opt, default 0)
%
%   newModel        an updated model structure
%
% 	This function only works if the draft model and source model follow
%	the same metabolite and compartment naming convention. Metabolites are
%	only matched by metaboliteName[compartment]. Useful if one wants to copy
%	additional reactions from source to draft after getModelFromHomology was
%	used involving the same models.
%
%   Usage: newModel=addRxnsGenesMets(model,sourceModel,rxns,addGene,rxnNote,confidence)

if nargin<6
    confidence=0;
end
rxns=convertCharArray(rxns);
if nargin<5
    rxnNote={'Added via addRxnsGenesMets()'};
else
    rxnNote=convertCharArray(rxnNote);
end
if numel(rxnNote)==1 && numel(rxns)>1
    rxnNoteArray=cell(1,numel(rxns));
    rxnNoteArray(:)=rxnNote;
    rxnNote=rxnNoteArray;
end
if nargin<4
    addGene=false;
end

% Obtain indexes of reactions in source model
[notNewRxn,oldRxn]=ismember(rxns,model.rxns);
rxns=rxns(~notNewRxn);
if isempty(rxns)
    error('All reactions are already in the model.');
elseif ~isempty(notNewRxn)
    fprintf('The following reactions were already present in the model and will not be added:\n')
    fprintf([strjoin(model.rxns(oldRxn(find(oldRxn))),'\n') '\n'])
end

[~, rxnIdx]=ismember(rxns,sourceModel.rxns); % Get rxnIDs
if any(rxnIdx==0)
    error('Not all reaction IDs could be found in the source model')
end

% Add new metabolites
metIdx=find(any(sourceModel.S(:,rxnIdx),2)); % Get metabolite IDs
% Many of the metabolites in are already in the draft model, so only add
% the new metabolites

% Match by metNames[metComps]. First make these structures for each model.
metCompsN =cellstr(num2str(model.metComps));
map=containers.Map(cellstr(num2str(transpose([1:length(model.comps)]))),model.comps);
metCompsN = map.values(metCompsN);
metCompsN = strcat(model.metNames,'[',metCompsN,']');

sourcemetCompsN=cellstr(num2str(sourceModel.metComps));
map=containers.Map(cellstr(num2str(transpose([1:length(sourceModel.comps)]))),sourceModel.comps);
sourcemetCompsN = map.values(sourcemetCompsN);
sourcemetCompsN=strcat(sourceModel.metNames,'[',sourcemetCompsN,']');

newMetCompsN=sourcemetCompsN(metIdx);
notNewMet=newMetCompsN(ismember(newMetCompsN,metCompsN));

if ~isempty(notNewMet)
    fprintf('\nThe following metabolites were already present in the model and will not be added:\n')
    fprintf([strjoin(transpose(notNewMet),'\n') '\n'])
end

metIdx=metIdx(~ismember(sourcemetCompsN(metIdx),metCompsN));

if ~isempty(metIdx)
    fprintf('\nThe following metabolites will be added to the model:\n')
    fprintf([strjoin(transpose(sourcemetCompsN(metIdx)),'\n') '\n'])

    if isfield(sourceModel,'mets')
        metsToAdd.mets=sourceModel.mets(metIdx);
    end
    if isfield(sourceModel,'metNames')
        metsToAdd.metNames=sourceModel.metNames(metIdx);
    end
    if isfield(sourceModel,'metFormulas')
        metsToAdd.metFormulas=sourceModel.metFormulas(metIdx);
    end
    if isfield(sourceModel,'metCharge')
        metsToAdd.metCharge=sourceModel.metCharge(metIdx);
    end
    if isfield(sourceModel,'metMiriams')
        metsToAdd.metMiriams=sourceModel.metMiriams(metIdx);
    end
    if isfield(sourceModel,'metFormulas')
        metsToAdd.metFormulas=sourceModel.metFormulas(metIdx);
    end
    if isfield(sourceModel,'inchis')
        metsToAdd.inchis=sourceModel.inchis(metIdx);
    end
    if isfield(sourceModel,'metSmiles')
        metsToAdd.metSmiles=sourceModel.metSmiles(metIdx);
    end

    metsToAdd.compartments=strtrim(cellstr(num2str(sourceModel.metComps(metIdx)))); % Convert from compartment string to compartment number
    [~,idx]=ismember(metsToAdd.compartments,strsplit(num2str(1:length(sourceModel.comps)))); % Match compartment number to compartment abbreviation
    metsToAdd.compartments=sourceModel.comps(idx); % Fill in compartment abbreviations

    model=addMets(model,metsToAdd);
end
fprintf('\nNumber of metabolites added to the model:\n')
fprintf([num2str(numel(metIdx)),'\n'])

% Add new genes
if ~islogical(addGene) | addGene ~= false
    if ischar(addGene)
        rxnToAdd.grRules={addGene};
    elseif iscell(addGene)
        rxnToAdd.grRules=addGene(~notNewRxn);
    else
        rxnToAdd.grRules=sourceModel.grRules(rxnIdx); % Get the relevant grRules
    end
end
% Add new reactions
rxnToAdd.equations=constructEquations(sourceModel,rxnIdx);
rxnToAdd.rxnNames=sourceModel.rxnNames(rxnIdx);
rxnToAdd.rxns=sourceModel.rxns(rxnIdx);
rxnToAdd.lb=sourceModel.lb(rxnIdx);
rxnToAdd.ub=sourceModel.ub(rxnIdx);
rxnToAdd.rxnNotes(:)=rxnNote(~notNewRxn);
rxnToAdd.rxnConfidenceScores=NaN(1,numel(rxnToAdd.rxns));
if ~isnumeric(confidence)
    EM='confidence score must be numeric';
    dispEM(EM, true);
end
rxnToAdd.rxnConfidenceScores(:)=confidence;
if isfield(sourceModel,'subSystems')
    rxnToAdd.subSystems=sourceModel.subSystems(rxnIdx);
end
if isfield(sourceModel,'eccodes')
    rxnToAdd.eccodes=sourceModel.eccodes(rxnIdx);
end
model=addRxns(model,rxnToAdd,3,'',false,true);

fprintf('\nNumber of reactions added to the model:\n')
fprintf([num2str(numel(rxnIdx)),'\n'])
end

author: SysBioChalmers

repo: RAVEN


text:
function GE_dumpHeader(IFileName, OutFile)
%
% GE_dumpHeader(IFileName,OutFile)
% Dumps much of the header from a GE lx2 or 5.X file to
% OutFile.  If OutFile is not specified, the output goes
% to stdout.
%
% Souheil J. Inati, PhD
% Dartmouth College
% May 2000
% souheil.inati@dartmouth.edu
%

%%%% Call GE_readHeader %%%%
[su_hdr,ex_hdr,se_hdr,im_hdr,pix_hdr] = GE_readHeader(IFileName);

% Open the output file
if nargin == 2
  outid = fopen(OutFile,'w');
else
  outid = 1;  % stdout
end

%%%% Exam Header %%%%
fprintf(outid,'Exam Header Section:\n');
fprintf(outid,GE_dumpExamHeader(ex_hdr));
fprintf(outid,'\n\n');

%%%% Series Header %%%%
fprintf(outid,'Series Header Section:\n');
fprintf(outid,GE_dumpSeriesHeader(se_hdr));
fprintf(outid,'\n\n');

%%%% Image Header %%%%
fprintf(outid,'Image Header Section:\n');
fprintf(outid,GE_dumpImageHeader(im_hdr));

% Close the File if not stdout
if outid ~= 1
  fclose(outid);
end

return

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Fig. 5.24  Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
% script for Figure 5.24, Lead design.
np=1;
dp=[1 1 0];
nc=[1 2];
dc=[1 10];
nol=conv(np,nc);
dol=conv(dp,dc);
rlocus(nol,dol)
axis([-19 1 -7.5 7.5])
title('Figure 5.24  Root locus for lead design')
hold on
r=roots([1 11 80 140]);
plot(r,'*')
z=0:.1:.9;
 wn=2:2:19;
 sgrid(z, wn)
 hold off

author: Sable

repo: mcbench-benchmarks


text:
function [K, Kc] = sdlfmXsdrbfKernCompute(sdlfmKern, sdrbfKern, t1, t2, type)

% SDLFMXSDRBFKERNCOMPUTE Cross kernel between a SDLFM and a SDRBF kernels.
% FORMAT
% DESC computes cross kernel terms between a switching dynamical
% LFM kernels and a switching dynamical RBF kernel for the multiple output
% kernel.
% ARG sdlfmKern : the kernel structure associated with the SDLFM
% kernel.
% ARG sdrbfKern : the kernel structure associated with the SDRBF
% kernel.
% ARG t : inputs for which kernel is to be computed.
% RETURN K : block of values from kernel matrix. K is a cell where the
% number of cells is related to the number of latent forces per interval
% RETURN Kc : the kernel matrices grouped in cells. Suitable form for the
% sparse approximations
%
% FORMAT
% DESC computes cross kernel terms between a SDLFM and a SDRBF kernels for
% the multiple output kernel.
% ARG sdlfmKern : the kernel structure associated with the SDLFM
% kernel.
% ARG sdrbfKern : the kernel structure associated with the SDRBF
% kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% RETURN K : block of values from kernel matrix. K is a cell where the
% number of cells is related to the number of latent forces per interval
% RETURN Kc : the kernel matrices grouped in cells. Suitable form for the
% sparse approximations
%
% FORMAT
% DESC computes cross kernel terms a SDLFM and a SDRBF kernels for
% the multiple output kernel. The SDLFM kernel can correspond to Position
% (default), Velocity or Acceleration. The type of kernel to be computed
% is specified in 'type'.
% ARG sdlfmKern : the kernel structure associated with the SDLFM
% kernel.
% ARG sdrbfKern : the kernel structure associated with the SDRBF
% kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG type : specifies the type of kerne to be computed
% RETURN K : block of values from kernel matrix. K is a cell where the
% number of cells is related to the number of latent forces per interval
% RETURN Kc : the kernel matrices grouped in cells. Suitable form for the
% sparse approximations
%
% SEEALSO : sdlfmKernParamInit, sdlfmKernCompute, sdlfmKernParamInit
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if nargin < 5
    type = 'Pos';
    if nargin < 4
        t2 = t1;
    end
end

if size(t1, 2) > 1 || size(t2, 2) > 1
    error('Input can only have one column');
end

% compareInverseWidth = sdlfmKern.inverseWidth == sdrbfKern.inverseWidth;
% if sum(sum(compareInverseWidth))~=(sdlfmKern.nIntervals*sdlfmKern.nlfPerInt)
%     error('Kernels cannot be cross combined if they have different inverse widths.')
% end
% compareSwitchingTimes = sdlfmKern.switchingTimes == sdrbfKern.switchingTimes;
% if sum(sum(compareSwitchingTimes))~=sdlfmKern.nIntervals
%     error('Kernels cannot be cross combined if they have different switching points.')
% end

switch type
    case 'Pos'
        fhandle = 'sdlfmXsdrbfKernComputeBlock';
    case 'Vel'
        fhandle = 'sdlfmvXsdrbfKernComputeBlock';
    case 'Accel'
        fhandle = 'sdlfmaXsdrbfKernComputeBlock';
end

fhandle = str2func(fhandle);

%Form the basic kernels
lfmKern = struct();
rbfKern = struct();

% Create structures that will make easy the computation of the kernels
%spVector = [sdlfmKern.switchingTimes t1(end)+0.1]; % For the last interval include the last point

spVector = [cumsum(sdlfmKern.switchingTimes) t1(end)+50];

dim1 = zeros(1, sdlfmKern.nIntervals); dim2 = zeros(1, sdrbfKern.nIntervals);

for i =1:sdlfmKern.nIntervals
    for j=1:sdlfmKern.nlfPerInt
        % Create the appropriate set of kernel structures
        lfmKern(i,j).mass = sdlfmKern.mass;
        lfmKern(i,j).spring = sdlfmKern.spring;
        lfmKern(i,j).damper = sdlfmKern.damper;
        lfmKern(i,j).inverseWidth = sdlfmKern.inverseWidth(j,i);
        lfmKern(i,j).sensitivity = sdlfmKern.sensitivity(j,i);
        rbfKern(i,j).variance = sdrbfKern.variance;
        rbfKern(i,j).inverseWidth = sdrbfKern.inverseWidth(j,i);
        lfmKern(i,j).limit = spVector(i+1) - spVector(i);
        rbfKern(i,j).limit = spVector(i+1) - spVector(i);
        lfmKern(i,j).isNormalised = sdlfmKern.isNormalised;
        rbfKern(i,j).isNormalised = sdrbfKern.isNormalised;
    end
    newt1 = t1(t1> spVector(i) & t1<spVector(i+1));
    newt2 = t2(t2> spVector(i) & t2<spVector(i+1));
    dim1(i) = length(newt1);
    dim2(i) = length(newt2);
end

if sum(dim1)~=length(t1) || sum(dim2)~=length(t2)
    error('A problem with the dimensions of the switching intervals occured')
end

% Compute some necessary constants

generalConst = sdlfmKernComputeConstant(sdlfmKern.nIntervals, ...
    lfmKern(1,1), lfmKern(1,1), spVector);

K = cell(sdrbfKern.nlfPerInt,1);
Kc = cell(sdrbfKern.nlfPerInt,1);

for q=1:sdrbfKern.nlfPerInt
    startValOne = 1;
    endValOne   = 0;
    tK = zeros(sum(dim1), sum(dim2));
    for i=1:sdlfmKern.nIntervals
        endValOne = endValOne + dim1(i);
        startValThree = 1;
        endValThree = 0;
        for j=1:i
            if i>j
                lfmKernLocal = lfmKern(j,q);
                rbfKernLocal = rbfKern(j,q);
            else
                lfmKernLocal = lfmKern(i,q);
                rbfKernLocal = rbfKern(i,q);
            end
            endValThree = endValThree + dim2(j);
            tK(startValOne:endValOne, startValThree:endValThree) = fhandle(lfmKernLocal, ...
                rbfKernLocal, t1(startValOne:endValOne) - spVector(i), ...
                t2(startValThree:endValThree) - spVector(j), i, j, generalConst);
            startValThree = endValThree + 1;
        end
        startValOne = endValOne + 1;
    end
    tK = real(tK);
    K{q} = tK;
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function divdif_test08 ( )

%*****************************************************************************80
%
%% DIVDIF_TEST17 tests NCO_RULE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 June 2011
%
%  Author:
%
%    John Burkardt
%
  norder = 8;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'DIVDIF_TEST17\n' );
  fprintf ( 1, '  NCO_RULE computes open Newton Cotes formulas.\n' );
  fprintf ( 1, '\n' );

  [ xtab, weight ] = nco_rule ( norder );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Newton-Cotes Open Quadrature Rule:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      Abscissa       Weight\n' );
  fprintf ( 1, '\n' );

  for i = 1 : norder
    fprintf ( 1, '  %3d  %14f  %14f\n', i, xtab(i), weight(i) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [dat] = read_edf(filename, hdr, begsample, endsample, chanindx)

% READ_EDF reads specified samples from an EDF continous datafile
% It neglects all trial boundaries as if the data was acquired in
% non-continous mode.
%
% Use as
%   [hdr] = read_edf(filename);
% where
%    filename        name of the datafile, including the .bdf extension
% This returns a header structure with the following elements
%   hdr.Fs           sampling frequency
%   hdr.nChans       number of channels
%   hdr.nSamples     number of samples per trial
%   hdr.nSamplesPre  number of pre-trigger samples in each trial
%   hdr.nTrials      number of trials
%   hdr.label        cell-array with labels of each channel
%   hdr.orig         detailled EDF header information
%
% Or use as
%   [dat] = read_edf(filename, hdr, begsample, endsample, chanindx);
% where
%    filename        name of the datafile, including the .bdf extension
%    hdr             header structure, see above
%    begsample       index of the first sample to read
%    endsample       index of the last sample to read
%    chanindx        index of channels to read (optional, default is all)
% This returns a Nchans X Nsamples data matrix

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

needhdr = (nargin==1);
needevt = (nargin==2);
needdat = (nargin>3);

if needhdr
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % read the header, this code is from EEGLAB's openbdf
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  FILENAME = filename;

  % defines Seperator for Subdirectories
  SLASH='/';
  BSLASH=char(92);

  cname=computer;
  if cname(1:2)=='PC' SLASH=BSLASH; end;

  fid=fopen(FILENAME,'r','ieee-le');
  if fid<0
    fprintf(2,['Error LOADEDF: File ' FILENAME ' not found\n']);
    return;
  end;

  EDF.FILE.FID=fid;
  EDF.FILE.OPEN = 1;
  EDF.FileName = FILENAME;

  PPos=min([max(find(FILENAME=='.')) length(FILENAME)+1]);
  SPos=max([0 find((FILENAME=='/') | (FILENAME==BSLASH))]);
  EDF.FILE.Ext = FILENAME(PPos+1:length(FILENAME));
  EDF.FILE.Name = FILENAME(SPos+1:PPos-1);
  if SPos==0
    EDF.FILE.Path = pwd;
  else
    EDF.FILE.Path = FILENAME(1:SPos-1);
  end;
  EDF.FileName = [EDF.FILE.Path SLASH EDF.FILE.Name '.' EDF.FILE.Ext];

  H1=char(fread(EDF.FILE.FID,256,'char')');     %
  EDF.VERSION=H1(1:8);                          % 8 Byte  Versionsnummer
  %if 0 fprintf(2,'LOADEDF: WARNING  Version EDF Format %i',ver); end;
  EDF.PID = deblank(H1(9:88));                  % 80 Byte local patient identification
  EDF.RID = deblank(H1(89:168));                % 80 Byte local recording identification
  %EDF.H.StartDate = H1(169:176);               % 8 Byte
  %EDF.H.StartTime = H1(177:184);               % 8 Byte
  EDF.T0=[str2num(H1(168+[7 8])) str2num(H1(168+[4 5])) str2num(H1(168+[1 2])) str2num(H1(168+[9 10])) str2num(H1(168+[12 13])) str2num(H1(168+[15 16])) ];

  % Y2K compatibility until year 2090
  if EDF.VERSION(1)=='0'
    if EDF.T0(1) < 91
      EDF.T0(1)=2000+EDF.T0(1);
    else
      EDF.T0(1)=1900+EDF.T0(1);
    end;
  else ;
    % in a future version, this is hopefully not needed
  end;

  EDF.HeadLen = str2num(H1(185:192));  % 8 Byte  Length of Header
  % reserved = H1(193:236);            % 44 Byte
  EDF.NRec = str2num(H1(237:244));     % 8 Byte  # of data records
  EDF.Dur = str2num(H1(245:252));      % 8 Byte  # duration of data record in sec
  EDF.NS = str2num(H1(253:256));       % 8 Byte  # of signals

  EDF.Label = char(fread(EDF.FILE.FID,[16,EDF.NS],'char')');
  EDF.Transducer = char(fread(EDF.FILE.FID,[80,EDF.NS],'char')');
  EDF.PhysDim = char(fread(EDF.FILE.FID,[8,EDF.NS],'char')');

  EDF.PhysMin= str2num(char(fread(EDF.FILE.FID,[8,EDF.NS],'char')'));
  EDF.PhysMax= str2num(char(fread(EDF.FILE.FID,[8,EDF.NS],'char')'));
  EDF.DigMin = str2num(char(fread(EDF.FILE.FID,[8,EDF.NS],'char')'));
  EDF.DigMax = str2num(char(fread(EDF.FILE.FID,[8,EDF.NS],'char')'));

  % check validity of DigMin and DigMax
  if (length(EDF.DigMin) ~= EDF.NS)
    fprintf(2,'Warning OPENEDF: Failing Digital Minimum\n');
    EDF.DigMin = -(2^15)*ones(EDF.NS,1);
  end
  if (length(EDF.DigMax) ~= EDF.NS)
    fprintf(2,'Warning OPENEDF: Failing Digital Maximum\n');
    EDF.DigMax = (2^15-1)*ones(EDF.NS,1);
  end
  if (any(EDF.DigMin >= EDF.DigMax))
    fprintf(2,'Warning OPENEDF: Digital Minimum larger than Maximum\n');
  end
  % check validity of PhysMin and PhysMax
  if (length(EDF.PhysMin) ~= EDF.NS)
    fprintf(2,'Warning OPENEDF: Failing Physical Minimum\n');
    EDF.PhysMin = EDF.DigMin;
  end
  if (length(EDF.PhysMax) ~= EDF.NS)
    fprintf(2,'Warning OPENEDF: Failing Physical Maximum\n');
    EDF.PhysMax = EDF.DigMax;
  end
  if (any(EDF.PhysMin >= EDF.PhysMax))
    fprintf(2,'Warning OPENEDF: Physical Minimum larger than Maximum\n');
    EDF.PhysMin = EDF.DigMin;
    EDF.PhysMax = EDF.DigMax;
  end
  EDF.PreFilt= char(fread(EDF.FILE.FID,[80,EDF.NS],'char')');   %
  tmp = fread(EDF.FILE.FID,[8,EDF.NS],'char')'; %   samples per data record
  EDF.SPR = str2num(char(tmp));               % samples per data record

  fseek(EDF.FILE.FID,32*EDF.NS,0);

  EDF.Cal = (EDF.PhysMax-EDF.PhysMin)./(EDF.DigMax-EDF.DigMin);
  EDF.Off = EDF.PhysMin - EDF.Cal .* EDF.DigMin;
  tmp = find(EDF.Cal < 0);
  EDF.Cal(tmp) = ones(size(tmp));
  EDF.Off(tmp) = zeros(size(tmp));

  EDF.Calib=[EDF.Off';(diag(EDF.Cal))];
  %EDF.Calib=sparse(diag([1; EDF.Cal]));
  %EDF.Calib(1,2:EDF.NS+1)=EDF.Off';

  EDF.SampleRate = EDF.SPR / EDF.Dur;

  EDF.FILE.POS = ftell(EDF.FILE.FID);
  if EDF.NRec == -1                            % unknown record size, determine correct NRec
    fseek(EDF.FILE.FID, 0, 'eof');
    endpos = ftell(EDF.FILE.FID);
    EDF.NRec = floor((endpos - EDF.FILE.POS) / (sum(EDF.SPR) * 2));
    fseek(EDF.FILE.FID, EDF.FILE.POS, 'bof');
    H1(237:244)=sprintf('%-8i',EDF.NRec);      % write number of records
  end;

  EDF.Chan_Select=(EDF.SPR==max(EDF.SPR));
  for k=1:EDF.NS
    if EDF.Chan_Select(k)
      EDF.ChanTyp(k)='N';
    else
      EDF.ChanTyp(k)=' ';
    end;
    if findstr(upper(EDF.Label(k,:)),'ECG')
      EDF.ChanTyp(k)='C';
    elseif findstr(upper(EDF.Label(k,:)),'EKG')
      EDF.ChanTyp(k)='C';
    elseif findstr(upper(EDF.Label(k,:)),'EEG')
      EDF.ChanTyp(k)='E';
    elseif findstr(upper(EDF.Label(k,:)),'EOG')
      EDF.ChanTyp(k)='O';
    elseif findstr(upper(EDF.Label(k,:)),'EMG')
      EDF.ChanTyp(k)='M';
    end;
  end;

  EDF.AS.spb = sum(EDF.SPR);    % Samples per Block
  bi=[0;cumsum(EDF.SPR)];

  idx=[];idx2=[];
  for k=1:EDF.NS,
    idx2=[idx2, (k-1)*max(EDF.SPR)+(1:EDF.SPR(k))];
  end;
  maxspr=max(EDF.SPR);
  idx3=zeros(EDF.NS*maxspr,1);
  for k=1:EDF.NS, idx3(maxspr*(k-1)+(1:maxspr))=bi(k)+ceil((1:maxspr)'/maxspr*EDF.SPR(k));end;

  %EDF.AS.bi=bi;
  EDF.AS.IDX2=idx2;
  %EDF.AS.IDX3=idx3;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % convert the header to Fieldtrip-style
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if all(EDF.SampleRate==EDF.SampleRate(1))
    hdr.Fs           = EDF.SampleRate(1);
    hdr.nChans       = EDF.NS;
    hdr.label        = cellstr(EDF.Label);
    % it is continuous data, therefore append all records in one trial
    hdr.nSamples     = EDF.Dur * EDF.SampleRate(1);
    hdr.nSamplesPre  = 0;
    hdr.nTrials      = EDF.NRec;
    hdr.orig         = EDF;

  elseif all(EDF.SampleRate(1:end-1)==EDF.SampleRate(1))
    % only the last channel has a deviant sampling frequency
    % this is the case for EGI recorded datasets that have been converted
    % to EDF+, in which case the annotation channel is the last
    chansel = find(EDF.SampleRate==EDF.SampleRate(1));
    % continue with the subset of channels that has a consistent sampling frequency
    hdr.Fs           = EDF.SampleRate(chansel(1));
    hdr.nChans       = length(chansel);
    warning('Skipping "%s" as continuous data channel because of inconsistent sampling frequency (%g Hz)', deblank(EDF.Label(end,:)), EDF.SampleRate(end));
    hdr.label        = cellstr(EDF.Label);
    hdr.label        = hdr.label(chansel);
    % it is continuous data, therefore append all records in one trial
    hdr.nSamples     = EDF.Dur * EDF.SampleRate(chansel(1));
    hdr.nSamplesPre  = 0;
    hdr.nTrials      = EDF.NRec;
    hdr.orig         = EDF;
    % this will be used on subsequent reading of data
    hdr.orig.chansel    = chansel;
    hdr.orig.annotation = find(strcmp(cellstr(hdr.orig.Label), 'EDF Annotations'));

  else
    % select the sampling rate that results in the most channels
    [a, b, c] = unique(EDF.SampleRate);
    for i=1:length(a)
      chancount(i) = sum(c==i);
    end
    [dum, indx] = max(chancount);
    chansel = find(EDF.SampleRate == a(indx));

    % continue with the subset of channels that has a consistent sampling frequency
    hdr.Fs           = EDF.SampleRate(chansel(1));
    hdr.nChans       = length(chansel);
    hdr.label        = cellstr(EDF.Label);
    hdr.label        = hdr.label(chansel);
    % it is continuous data, therefore append all records in one trial
    hdr.nSamples     = EDF.Dur * EDF.SampleRate(chansel(1));
    hdr.nSamplesPre  = 0;
    hdr.nTrials      = EDF.NRec;
    hdr.orig         = EDF;
    % this will be used on subsequent reading of data
    hdr.orig.chansel    = chansel;
    hdr.orig.annotation = find(strcmp(cellstr(hdr.orig.Label), 'EDF Annotations'));

    warning('channels with different sampling rate not supported, using a subselection of %d channels at %f Hz', length(hdr.label), hdr.Fs);
  end

  % return the header
  dat = hdr;

elseif needdat || needevt
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % read the data
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % retrieve the original header
  EDF = hdr.orig;

  % determine whether a subset of channels should be used
  % which is the case if channels have a variable sampling frequency
  variableFs = isfield(EDF, 'chansel');

  if variableFs && needevt
    % read the annotation channel, not the data channels
    EDF.chansel = EDF.annotation;
    begsample = 1;
    endsample = EDF.SampleRate(end)*EDF.NRec*EDF.Dur;
  end

  if variableFs
    epochlength = EDF.Dur * EDF.SampleRate(EDF.chansel(1));   % in samples for the selected channel
    blocksize   = sum(EDF.Dur * EDF.SampleRate);              % in samples for all channels
    chanoffset  = EDF.Dur * EDF.SampleRate;
    chanoffset  = cumsum([0; chanoffset(1:end-1)]);
    % use a subset of channels
    nchans = length(EDF.chansel);
  else
    epochlength = EDF.Dur * EDF.SampleRate(1);                % in samples for a single channel
    blocksize   = sum(EDF.Dur * EDF.SampleRate);              % in samples for all channels
    % use all channels
    nchans = EDF.NS;
  end

  % determine the trial containing the begin and end sample
  begepoch    = floor((begsample-1)/epochlength) + 1;
  endepoch    = floor((endsample-1)/epochlength) + 1;
  nepochs     = endepoch - begepoch + 1;

  if nargin<5
    chanindx = 1:nchans;
  end

  % allocate memory to hold the data
  dat = zeros(length(chanindx),nepochs*epochlength);

  % read and concatenate all required data epochs
  for i=begepoch:endepoch
    if variableFs
      % only a subset of channels with consistent sampling frequency is read
      offset = EDF.HeadLen + (i-1)*blocksize*2; % in bytes
      % read the complete data block
      buf = readLowLevel(filename, offset, blocksize); % see below in subfunction
      for j=1:length(chanindx)
        % cut out the part that corresponds with a single channel
        dat(j,((i-begepoch)*epochlength+1):((i-begepoch+1)*epochlength)) = buf((1:epochlength) + chanoffset(EDF.chansel(chanindx(j))));
      end

    elseif length(chanindx)==1
      % this is more efficient if only one channel has to be read, e.g. the status channel
      offset = EDF.HeadLen + (i-1)*blocksize*2; % in bytes
      offset = offset + (chanindx-1)*epochlength*2;
      % read the data for a single channel
      buf = readLowLevel(filename, offset, epochlength); % see below in subfunction
      dat(:,((i-begepoch)*epochlength+1):((i-begepoch+1)*epochlength)) = buf;

    else
      % read the data from all channels, subsequently select the desired channels
      offset = EDF.HeadLen + (i-1)*blocksize*2; % in bytes
      % read the complete data block
      buf = readLowLevel(filename, offset, blocksize); % see below in subfunction
      buf = reshape(buf, epochlength, nchans);
      dat(:,((i-begepoch)*epochlength+1):((i-begepoch+1)*epochlength)) = buf(:,chanindx)';
    end
  end

  % select the desired samples
  begsample = begsample - (begepoch-1)*epochlength;  % correct for the number of bytes that were skipped
  endsample = endsample - (begepoch-1)*epochlength;  % correct for the number of bytes that were skipped
  dat = dat(:, begsample:endsample);

  % Calibrate the data
  if variableFs
    % using a sparse matrix speeds up the multiplication
    calib = sparse(diag(EDF.Cal(EDF.chansel(chanindx))));
    dat   = full(calib * dat);
  elseif length(chanindx)==1
    % in case of one channel the calibration would result in a sparse array
    calib = EDF.Cal(chanindx);
    dat   = calib * dat;
  else
    % using a sparse matrix speeds up the multiplication
    calib = sparse(diag(EDF.Cal(chanindx)));
    dat   = calib * dat;
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION for reading the 16 bit values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function buf = readLowLevel(filename, offset, numwords)
if offset < 2*1024^2
  % use the external mex file, only works for <2GB
  buf = read_16bit(filename, offset, numwords);
else
  % use plain matlab, thanks to Philip van der Broek
  fp = fopen(filename,'r','ieee-le');
  status = fseek(fp, offset, 'bof');
  if status
    error(['failed seeking ' filename]);
  end
  [buf,num] = fread(fp,numwords,'bit16=>double');
  fclose(fp);
  if (num<numwords)
    error(['failed reading ' filename]);
    return
  end
end

author: PatternRecognition

repo: OpenBMI


text:
% Wrapper around Matlab's histc function with some modifications
%
% This function assigns each input value a bin number according to specified
% limits of values and width of each bin.
% Vector of bin edges is created out of specified limits and bin width. Value
% x(i) is assigned to bin k if edges(k) <= x(i) < edges(k+1). So far this is
% exactly what Matlab's function histc does. The difference comes in the last bin.
% Condition of assignment for the last bin is edges(end-1) <= x(i) <= edges(end).
%
%  USAGE
%   [bins, nBins, edges] = helpers.bin(x, limits, binWidth)
%   x           1D vector of values that should be binned.
%   limits      [min_x max_x] limit values that define bins. These need not be equal to
%               min(x) and max(x). Values of x that equal to min_x or max_x are all
%               included in the final histogram.
%   binWidth    width of each individual bin. Used to create edges vector.
%   bins        Vector of assigned bin numbers to x. length(bins) == length(x).
%               Values outside of limits are assigned to NaN.
%   nBins       Number of bins.
%   edges       Edge values of bins.
%
function [bins, nBins, edges] = bin(x, limits, binWidth)
    nBins = ceil((limits(2) - limits(1)) / binWidth);

    edges = limits(1):binWidth:limits(2);
    if length(edges) == nBins
        edges(end+1) = limits(2);
    end

    ind = (x == edges(end)); % find indices of values that according to mathematics should be
                                 % assigned to bin nBins+1
    [~, bins] = histc(x, edges);
    bins(ind) = nBins;
    bins(bins > nBins) = nan;
    bins(bins == 0) = nan;
end

author: kavli-ntnu


text:
function test_bug2570

% WALLTIME 00:10:00
% MEM 1gb
% DEPENDENCY ft_apply_montage ft_scalingfactor

montage = [];
montage.tra = 1e6;
montage.labelold = {'Cz'};
montage.labelnew = {'Cz'};
montage.chanunitold = {'V'};
montage.chanunitnew = {'uV'};

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

output = ft_apply_montage(montage, montage);
assert(isequal(output, montage)); % it should not have changed

output = ft_apply_montage(montage, montage, 'inverse', 1);
assert(isequal(output.tra, 1)); % this should be an identity transform
assert(isequal(output.chanunitold, output.chanunitnew));

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

data_V = [];
data_V.label = {'Cz'};
data_V.time = {1:10};
data_V.trial = {ones(1,10)};
data_V.chanunit = {'V'};

data_mV = [];
data_mV.label = {'Cz'};
data_mV.time = {1:10};
data_mV.trial = {ones(1,10)*1000};
data_mV.chanunit = {'mV'};

data_uV = [];
data_uV.label = {'Cz'};
data_uV.time = {1:10};
data_uV.trial = {ones(1,10)*1000000};
data_uV.chanunit = {'uV'};

output1 = ft_apply_montage(data_V,  montage);
output2 = ft_apply_montage(data_mV, montage);
output3 = ft_apply_montage(data_uV, montage);

assert(isequal(output1, output2));
assert(isequal(output1, output3));
assert(isequal(output1, data_uV)); % they should all be in uV

output1 = ft_apply_montage(data_V,  montage, 'inverse', 'yes');
output2 = ft_apply_montage(data_mV, montage, 'inverse', 'yes');
output3 = ft_apply_montage(data_uV, montage, 'inverse', 'yes');

assert(isequal(output1, output2));
assert(isequal(output1, output3));
assert(isequal(output1, data_V)); % they should all be in V

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

data_T = [];
data_T.label = {'Cz'};
data_T.time = {1:10};
data_T.trial = {ones(1,10)};
data_T.chanunit = {'T'};

try
  ft_apply_montage(data_T,  montage);
  ok = false;
catch
  ok = true;
end
if ~ok
  error('this should have given an error, T cannot be converted to V');
end


text:
function [eigvector, eigvalue, bSuccess] = OLPP(W, options, data)
% OLPP: Orthogonal Locality Preserving Projections
%
%       [eigvector, eigvalue, bSuccess] = OLPP(W, options, data)
%
%             Input:
%               data       - Data matrix. Each row vector of fea is a data point.
%               W       - Affinity matrix. You can either call "constructW"
%                         to construct the W, or construct it by yourself.
%               options - Struct value in Matlab. The fields in options
%                         that can be set:
%
%                         Please see OLGE.m for other options.
%
%             Output:
%               eigvector - Each column is an embedding function, for a new
%                           data point (row vector) x,  y = x*eigvector
%                           will be the embedding result of x.
%               eigvalue  - The sorted eigvalue of OLPP eigen-problem.
%
%               bSuccess  - 0 or 1. Indicates whether the OLPP calcuation
%                           is successful. (OLPP needs matrix inverse,
%                           which will lead to eigen-decompose a
%                           non-symmetrical matrix. The caculation precsion
%                           of malab sometimes will cause imaginary numbers
%                           in eigenvectors. It seems that the caculation
%                           precsion of matlab is a little bit random, you
%                           can try again if not successful. More robust
%                           and efficient algorithms are welcome!)
%
%                         Please see OLGE.m for other options.
%
%
%
%
%    Examples:
%
%       fea = rand(50,70);
%       options = [];
%       options.Metric = 'Euclidean';
%       options.NeighborMode = 'KNN';
%       options.k = 5;
%       options.WeightMode = 'HeatKernel';
%       options.t = 1;
%       W = constructW(fea,options);
%       options.PCARatio = 0.99
%       options.ReducedDim = 5;
%       bSuccess = 0
%       while ~bSuccess
%           [eigvector, eigvalue, bSuccess] = OLPP(W, options, fea);
%       end
%       Y = fea*eigvector;
%
%       fea = rand(50,70);
%       gnd = [ones(10,1);ones(15,1)*2;ones(10,1)*3;ones(15,1)*4];
%       options = [];
%       options.Metric = 'Euclidean';
%       options.NeighborMode = 'Supervised';
%       options.gnd = gnd;
%       options.bLDA = 1;
%       W = constructW(fea,options);
%       options.PCARatio = 1;
%       options.ReducedDim = 5;
%       while ~bSuccess
%           [eigvector, eigvalue, bSuccess] = OLPP(W, options, fea);
%       end
%       Y = fea*eigvector;
%
%
%
% See also constructW, LPP, LGE, OLGE.
%
%Reference:
%
%   Deng Cai and Xiaofei He, "Orthogonal Locality Preserving Indexing"
%   The 28th Annual International ACM SIGIR Conference (SIGIR'2005),
%   Salvador, Brazil, Aug. 2005.
%
%   Deng Cai, Xiaofei He, Jiawei Han and Hong-Jiang Zhang, "Orthogonal
%   Laplacianfaces for Face Recognition". IEEE Transactions on Image
%   Processing, vol. 15, no. 11, pp. 3608-3614, November, 2006.
%
%    Written by Deng Cai (dengcai2 AT cs.uiuc.edu), August/2004, Feb/2006,
%                                                   Mar/2007, May/2007

if (~exist('options','var'))
   options = [];
end

[nSmp] = size(data,1);
if size(W,1) ~= nSmp
    error('W and data mismatch!');
end

D = full(sum(W,2));

if isfield(options,'Regu') && options.Regu
    options.ReguAlpha = options.ReguAlpha*sum(D)/length(D);
end

D = sparse(1:nSmp,1:nSmp,D,nSmp,nSmp);

%==========================
% If data is too large, the following centering codes can be commented
%==========================
if isfield(options,'keepMean') && options.keepMean
else
    if issparse(data)
        data = full(data);
    end
    sampleMean = mean(data);
    data = (data - repmat(sampleMean,nSmp,1));
end
%==========================

[eigvector, eigvalue, bSuccess] = OLGE(W, D, options, data);

eigIdx = find(eigvalue < 1e-3);
eigvalue (eigIdx) = [];
eigvector(:,eigIdx) = [];

author: ZJULearning

repo: MatlabFunc


text:
function write(obj,filename,create_signature)
%WRITE Write sequence to file.
%   WRITE(seqObj, filename) Write the sequence data to the given
%   filename using the open file format for MR sequences.
%
%   Examples:
%   Write the sequence file to the my_sequences directory
%
%       write(seqObj,'my_sequences/gre.seq')
%
% See also  read

if (nargin<3)
    create_signature=true;
end

fid=fopen(filename, 'w');
assert(fid ~= -1, 'Cannot open file: %s', filename);
fprintf(fid, '# Pulseq sequence file\n');
fprintf(fid, '# Created by MATLAB mr toolbox\n\n');

fprintf(fid, '[VERSION]\n');
fprintf(fid, 'major %s\n', num2str(obj.version_major));
fprintf(fid, 'minor %s\n', num2str(obj.version_minor));
fprintf(fid, 'revision %s\n', num2str(obj.version_revision));
fprintf(fid, '\n');

if ~isempty(obj.definitions)
    fprintf(fid, '[DEFINITIONS]\n');
    keys = obj.definitions.keys;
    values = obj.definitions.values;
    for i=1:length(keys)
        fprintf(fid, '%s ', keys{i});
        if (ischar(values{i}))
            fprintf(fid, '%s ', values{i});
        else
            fprintf(fid, '%.9g ', values{i});
        end
        fprintf(fid, '\n');
    end
	fprintf(fid, '\n');
end

fprintf(fid, '# Format of blocks:\n');
fprintf(fid, '# NUM DUR RF  GX  GY  GZ  ADC  EXT\n');
fprintf(fid, '[BLOCKS]\n');
idFormatWidth = length(num2str(length(obj.blockEvents)));
idFormatStr = ['%' num2str(idFormatWidth) 'd'];
for i = 1:length(obj.blockEvents)
    %fprintf(fid,[idFormatStr ' %2d %2d %3d %3d %3d %2d 0\n'],[i obj.blockEvents(i,:)]);
    %fprintf(fid,[idFormatStr ' %2d %2d %3d %3d %3d %2d 0\n'],[i obj.blockEvents{i}]);
    bd=obj.blockDurations(i)/obj.blockDurationRaster;
    bdr=round(bd);
    assert(abs(bdr-bd)<1e-6); % this may still trigger false alarms for very long delays due to the limited accuracy of the double
    fprintf(fid,[idFormatStr ' %3d %3d %3d %3d %3d %2d %2d\n'], ...
            [i bdr obj.blockEvents{i}(2:end)]);
end
fprintf(fid, '\n');

if ~isempty(obj.rfLibrary.keys)
    fprintf(fid, '# Format of RF events:\n');
    fprintf(fid, '# id amplitude mag_id phase_id time_shape_id delay freq phase\n');
    fprintf(fid, '# ..        Hz   ....     ....          ....    us   Hz   rad\n');
    fprintf(fid, '[RF]\n');
    keys = obj.rfLibrary.keys;
    for k = keys
        libData1 = obj.rfLibrary.data(k).array(1:4);
        libData2 = obj.rfLibrary.data(k).array(6:7);
        delay = round(obj.rfLibrary.data(k).array(5)/obj.rfRasterTime)*obj.rfRasterTime*1e6; % a bit of a hack: round the delay
        fprintf(fid, '%d %12g %d %d %d %g %g %g\n', [k libData1 delay libData2]);
    end
    fprintf(fid, '\n');
end

arbGradMask = obj.gradLibrary.type == 'g';
trapGradMask = obj.gradLibrary.type == 't';

if any(arbGradMask)
    fprintf(fid, '# Format of arbitrary gradients:\n');
    fprintf(fid, '#   time_shape_id of 0 means default timing (stepping with grad_raster starting at 1/2 of grad_raster)\n');
    fprintf(fid, '# id amplitude amp_shape_id time_shape_id delay\n'); % do we need delay ???
    fprintf(fid, '# ..      Hz/m       ..         ..          us\n');
    fprintf(fid, '[GRADIENTS]\n');
    keys = obj.gradLibrary.keys;
    for k = keys(arbGradMask)
        fprintf(fid, '%d %12g %d %d %d\n', ...
                [k obj.gradLibrary.data(k).array(1:3) ...
                 round(obj.gradLibrary.data(k).array(4)*1e6)]);
    end
    fprintf(fid, '\n');
end

if any(trapGradMask)
    fprintf(fid, '# Format of trapezoid gradients:\n');
    fprintf(fid, '# id amplitude rise flat fall delay\n');
    fprintf(fid, '# ..      Hz/m   us   us   us    us\n');
    fprintf(fid, '[TRAP]\n');
    keys = obj.gradLibrary.keys;
    for k = keys(trapGradMask)
        data = obj.gradLibrary.data(k).array;
        data(2:end) = round(1e6*data(2:end));
        fprintf(fid, '%2d %12g %3d %4d %3d %3d\n', [k data]);
    end
    fprintf(fid, '\n');
end

if ~isempty(obj.adcLibrary.keys)
    fprintf(fid, '# Format of ADC events:\n');
    fprintf(fid, '# id num dwell delay freq phase\n');
    fprintf(fid, '# ..  ..    ns    us   Hz   rad\n');
    fprintf(fid, '[ADC]\n');
    keys = obj.adcLibrary.keys;
    for k = keys
        data = obj.adcLibrary.data(k).array(1:5).*[1 1e9 1e6 1 1];
        fprintf(fid, '%d %d %.0f %.0f %g %g\n', [k data]);
    end
    fprintf(fid, '\n');
end

%if ~isempty(obj.delayLibrary.keys)
%    fprintf(fid, '# Format of delays:\n');
%    fprintf(fid, '# id delay (us)\n');
%    fprintf(fid, '[DELAYS]\n');
%    keys = obj.delayLibrary.keys;
%    for k = keys
%        fprintf(fid, '%d %d\n', ...
%                [k round(1e6*obj.delayLibrary.data(k).array)]);
%    end
%    fprintf(fid, '\n');
%end

if ~isempty(obj.extensionLibrary.keys)
    fprintf(fid, '# Format of extension lists:\n');
    fprintf(fid, '# id type ref next_id\n');
    fprintf(fid, '# next_id of 0 terminates the list\n');
    fprintf(fid, '# Extension list is followed by extension specifications\n');
    fprintf(fid, '[EXTENSIONS]\n');
    keys = obj.extensionLibrary.keys;
    for k = keys
        fprintf(fid, '%d %d %d %d\n', ...
                [k round(obj.extensionLibrary.data(k).array)]);
    end
    fprintf(fid, '\n');
end

if ~isempty(obj.trigLibrary.keys)
    fprintf(fid, '# Extension specification for digital output and input triggers:\n');
    fprintf(fid, '# id type channel delay (us) duration (us)\n');
%     fprintf(fid, 'extension TRIGGERS 1\n'); % fixme: extension ID 1 is hardcoded here for triggers
    fprintf(fid, ['extension TRIGGERS ',num2str(obj.getExtensionTypeID('TRIGGERS')),'\n']);

    keys = obj.trigLibrary.keys;
    for k = keys
        fprintf(fid, '%d %d %d %d %d\n', ...
                [k round(obj.trigLibrary.data(k).array.*[1 1 1e6 1e6])]);
    end
    fprintf(fid, '\n');
end

if ~isempty(obj.labelsetLibrary.keys)
    lbls=mr.getSupportedLabels();

    fprintf(fid, '# Extension specification for setting labels:\n');
    fprintf(fid, '# id set labelstring\n');
    tid=obj.getExtensionTypeID('LABELSET');
    fprintf(fid, ['extension LABELSET ',num2str(tid),'\n']);
    keys = obj.labelsetLibrary.keys;
    for k = keys
            fprintf(fid, '%d %d %s\n', ...
                k, obj.labelsetLibrary.data(k).array(1),lbls{obj.labelsetLibrary.data(k).array(2)});
    end
    fprintf(fid, '\n');

    fprintf(fid, '# Extension specification for setting labels:\n');
    fprintf(fid, '# id set labelstring\n');
    tid=obj.getExtensionTypeID('LABELINC');
    fprintf(fid, ['extension LABELINC ',num2str(tid),'\n']);
    lbls=mr.getSupportedLabels();
    keys = obj.labelincLibrary.keys;
    for k = keys
            fprintf(fid, '%d %d %s\n', ...
                k, obj.labelincLibrary.data(k).array(1),lbls{obj.labelincLibrary.data(k).array(2)});
    end
    fprintf(fid, '\n');
end

if ~isempty(obj.shapeLibrary.keys)
    fprintf(fid, '# Sequence Shapes\n');
    fprintf(fid, '[SHAPES]\n\n');
    keys = obj.shapeLibrary.keys;
    for k = keys
        shape_dat = obj.shapeLibrary.data(k).array;
        fprintf(fid, 'shape_id %d\n', k);
        fprintf(fid, 'num_samples %d\n', shape_dat(1));
        fprintf(fid, '%.9g\n', shape_dat(2:end));
        fprintf(fid, '\n');
    end
end

fclose(fid);

if create_signature
    % sign the file (this version with re/loading the file is a factor 2 faster than the sprintf() based one that kept a memory-copy of the data written)

    % re-open and read in the file
    fid=fopen(filename, 'r');
    buf=fread(fid);
    fclose(fid);

    % calculate the digest
    md5hash=md5_java(buf);
    %fprintf('%s\n',md5hash);

    % store the signature in the object
    obj.signatureType='md5';
    obj.signatureFile='text';
    obj.signatureValue=md5hash;

    % re-open the file for appending
    fid=fopen(filename, 'a');
    fprintf(fid, '\n[SIGNATURE]\n'); % the preceding new line BELONGS to the signature (and needs to be sripped away to recalculate the signature)
    fprintf(fid, '# This is the hash of the Pulseq file, calculated right before the [SIGNATURE] section was added\n');
    fprintf(fid, '# It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE]\n');
    fprintf(fid, '# The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be sripped away for recalculating/verification)\n');
    fprintf(fid, 'Type md5\n');
    fprintf(fid, 'Hash %s\n', md5hash);
    fclose(fid);
end

end

function out=md5_java(buf)
    import java.security.*;
    import java.math.*;
    import java.lang.String;

    md = MessageDigest.getInstance('MD5');
    hash = md.digest(double(buf));
    bi = BigInteger(1, hash);

    out=char(String.format('%032x', bi));
end

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function f = mckinnon ( x )

%*****************************************************************************80
%
%% MCKINNON computes the McKinnon function.
%
%  Discussion:
%
%    This function has a global minimizer:
%
%      X* = ( 0.0, -0.5 ), F(X*) = -0.25.
%
%    There are three parameters, TAU, THETA and PHI.
%
%    1 < TAU, then F is strictly convex.
%             and F has continuous first derivatives.
%    2 < TAU, then F has continuous second derivatives.
%    3 < TAU, then F has continuous third derivatives.
%
%    However, this function can cause the Nelder-Mead optimization
%    algorithm to "converge" to a point which is not the minimizer
%    of the function F.
%
%    Sample parameter values which cause problems for Nelder-Mead
%    include:
%
%      TAU = 1, THETA = 15, PHI =  10;
%      TAU = 2, THETA =  6, PHI =  60;
%      TAU = 3, THETA =  6, PHI = 400;
%
%    To get the bad behavior, we also assume the initial simplex has the form
%
%      X1 = (0,0),
%      X2 = (1,1),
%      X3 = (A,B),
%
%    where
%
%      A = (1+sqrt(33))/8 =  0.84307...
%      B = (1-sqrt(33))/8 = -0.59307...
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ken McKinnon,
%    Convergence of the Nelder-Mead simplex method to a nonstationary point,
%    SIAM Journal on Optimization,
%    Volume 9, Number 1, 1998, pages 148-158.
%
%  Parameters:
%
%    Input, real X(2), the argument of the function.
%
%    Output, real F, the value of the function at X.
%
  if ( length ( x ) ~= 2 )
    error ( 'Error: function expects a two dimensional input\n' );
  end

  tau = 2.0;
  theta = 6.0;
  phi = 60.0;

  if ( x(1) <= 0.0 )
    f = theta * phi * abs ( x(1) ).^tau + x(2) * ( 1.0 + x(2) );
  else
    f = theta       *       x(1).^tau   + x(2) * ( 1.0 + x(2) );
  end

  return
end


text:
% GPCG algorithm adapted for NMF
function X = gpcg_nmf(Y,A,X,MaxIter,CostFun,Alpha)
% Our implementation of GPCG;

[R,T] = size(X); [M,T] = size(Y);
beta = 0.5; mu = 0.01; MaxSearchIter = 10;
Diff_DF_max = 0; gamma_GP = 0.01;

if CostFun == 1
   H = spalloc(R*T,R*T,T*R^2);
   Hx = A'*A;
   H = kron(speye(T),Hx); % Hessian
end

for k = 1:MaxIter

% Step 1
    Xp = X;
    if CostFun ~= 1
       H = CostFunHess(Y,A,X,Alpha); % Hessian
    end
    P = -CostFunGrad(Y,A,X,CostFun,Alpha); % Gradient Matrix
    p = P(:); % Vectorization

% Step 2
     eta = (p'*p)/(p'*(H*p));
     for i = 0:MaxSearchIter % Armijo rule
         eta = eta*beta.^i;
         X = max(X + eta*P,0);
         if (CostFunEval(Y,A,X,CostFun,Alpha) - CostFunEval(Y,A,Xp,CostFun,Alpha)) <=  (norm(X - Xp,'fro')*mu/eta)
             break;
         end
     end

% Step 3
  Z = zeros(size(X));
  Z(X > eps) = 1;
  z = vec(Z);

% Step 4
  Pc = CostFunGrad(Y,A,X,CostFun,Alpha); % Gradient Matrix
  pc = Pc(:); % Vectorization

% Step 5
  pR = z.*pc; % Reduced gradient
  if CostFun ~= 1
     H = CostFunHess(Y,A,X,Alpha); % Hessian
  end

HR = repmat(z,1,R*T).*H.*repmat(z',R*T,1) + speye(R*T) - spdiags(z,0,R*T,R*T); % reduced Hessian

% Step 6
 [pc,flag,relres,iter,resvec] = pcg(HR,-pR);
 P = reshape(pc,R,T); % Matricization

 % Step 7
     DF = CostFunEval(Y,A,X,CostFun,Alpha);
     eta = 1;
     for j = 0:MaxSearchIter % Armijo rule
         eta = eta*beta.^j;
         X = max(X + eta*P,0);
         if (CostFunEval(Y,A,X,CostFun,Alpha) < DF)
             break;
         end
     end

   Fn(k) = CostFunEval(Y,A,X,CostFun,Alpha);

 % Stopping criterion
   DF_old = norm(Y - A*Xp,'fro');  DF = norm(Y - A*X,'fro');
   Diff_DF = DF_old - DF;  Diff_DF_max = max(Diff_DF,Diff_DF_max);
   if (Diff_DF <= gamma_GP*Diff_DF_max) & (k > 1)
       break;
   end

end

% Cost Function
function F = CostFunEval(Y,A,X,CostFun,Alpha);

 switch CostFun

     case 1 % Eucliden distance

         F = norm(Y - A*X,'fro');

     case 2 % Alpha divergence

         Z = A*X + 1E2*eps;
         Y = Y + 1E2*eps;
         if Alpha == 1 % KL divergence
            F = sum(sum(Y.*log(Y./Z) + Z - Y));
         elseif Alpha == 0 % Dual KL divergence
            F = sum(sum(Z.*log(Z./(Y+eps)) + Y - Z));
         else % Alpha-divergence
            F = (1/Alpha)*sum(sum(Y.*( (Y./Z).^(Alpha - 1) - 1)/(Alpha - 1) + Z - Y));
         end
 end

% Gradient of Cost Function
function G = CostFunGrad(Y,A,X,CostFun,Alpha);

 switch CostFun

     case 1 % Eucliden distance

         G = A'*(A*X - Y);

     case 2 % Alpha divergence

         Z = A*X+1E2*eps;
         if ~Alpha % Dual KL divergence
            G = A'*log(Z./(Y + 1E2*eps));
         else
            G = (1/Alpha)*A'*(1 - ((Y+1E2*eps)./Z).^Alpha);
         end
  end

% Hessian of Cost Function
function H = CostFunHess(Y,A,X,Alpha);

[R,T] = size(X);
M = size(Y,1);
H = spalloc(R*T,R*T,T*R^2);
Z = A*X+1E2*eps;
if ~Alpha % Dual KL divergence
    Zx = 1./Z + 1E2*eps;
else
    Zx = ((Y+1E2*eps).^Alpha)./(Z.^(Alpha + 1));
end
 for t = 1:T
     H(((t-1)*R+1):t*R,((t-1)*R+1):t*R) = A'*repmat(Zx(:,t),1,M)*A;
 end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function sim = gaussianKernel(x1, x2, sigma)
%RBFKERNEL returns a radial basis function kernel between x1 and x2
%   sim = gaussianKernel(x1, x2) returns a gaussian kernel between x1 and x2
%   and returns the value in sim

% Ensure that x1 and x2 are column vectors
x1 = x1(:); x2 = x2(:);

% You need to return the following variables correctly.
sim = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the similarity between x1
%               and x2 computed using a Gaussian kernel with bandwidth
%               sigma
%
%

sim = exp((- sum((x1 - x2) .^ 2)) / (2 * (sigma ^ 2)));

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

end

author: UtkarshPathrabe


text:
function [lat,lon,h]=xyz2ell3(X,Y,Z,a,b,e2)
% XYZ2ELL3  Converts cartesian coordinates to ellipsoidal.
%   Uses direct algorithm in B.R. Bowring, "The accuracy of
%   geodetic latitude and height equations", Survey
%   Review, v28 #218, October 1985, pp.202-206.  Vectorized.
%   See also XYZ2ELL, XYZ2ELL2.
% Version: 2011-02-19
% Useage:  [lat,lon,h]=xyz2ell3(X,Y,Z,a,b,e2)
%          [lat,lon,h]=xyz2ell3(X,Y,Z)
% Input:   X \
%          Y  > vectors of cartesian coordinates in CT system (m)
%          Z /
%          a   - ref. ellipsoid major semi-axis (m); default GRS80
%          b   - ref. ellipsoid minor semi-axis (m); default GRS80
%          e2  - ref. ellipsoid eccentricity squared; default GRS80
% Output:  lat - vector of ellipsoidal latitudes (radians)
%          lon - vector of ellipsoidal longitudes (radians)
%          h   - vector of ellipsoidal heights (m)

% Copyright (c) 2011, Michael R. Craymer
% All rights reserved.
% Email: mike@craymer.com

if nargin ~= 3 & nargin ~= 6
  warning('Incorrect number of input arguments');
  return
end
if nargin == 3
  [a,b,e2]=refell('grs80');
end

lon=atan2(Y,X);
e=e2*(a/b)^2;
p=sqrt(X.*X+Y.*Y);
r=sqrt(p.*p+Z.*Z);
u=atan(b.*Z.*(1+e.*b./r)./(a.*p));
lat=atan((Z+e.*b.*sin(u).^3)./(p-e2.*a.*cos(u).^3));
v=a./sqrt(1-e2.*sin(lat).^2);
h=p.*cos(lat)+Z.*sin(lat)-a*a./v;


text:
function stackorder_menu(parent)
    % stackorder_menu context menu for sort order

    %FIXME doesn't seem to attach to all items
    c=uicontextmenu('Tag','StackOrderContext');
    set(parent,'UIContextMenu',c);
    uimenu(c,'Label','top'    , 'MenuSelectedFcn',@(s,~)setstack(s,'top'));
    uimenu(c,'Label','up'     , 'MenuSelectedFcn',@(s,~)setstack(s,'up'));
    uimenu(c,'Label','down'   , 'MenuSelectedFcn',@(s,~)setstack(s,'down'));
    uimenu(c,'Label','bottom' , 'MenuSelectedFcn',@(s,~)setstack(s,'bottom'));

    function setstack(source,val)
        uistack(source.Parent,val);
    end
end

author: CelsoReyes


text:
function new_name = fixIndex(name,destination_folderpath)
    index = 1;
    list = updateList(destination_folderpath);

    try_name = assembleName(name,index);

    for i = 1:length(list)
        if strcmpi(try_name,list(i).name)
            index = index + 1;
            try_name = assembleName(name,index);
        end
    end
    new_name = try_name;
end

function new_name = assembleName(name,index)
    us = strfind(name,'_');
    dots = strfind(name,'.');
    name(us(end)+1:dots(end)-1)=[];
    new_name = [name(1:us(end)), num2str(index), name(us(end)+1:end)];
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function neighbours = fixneighbours(neighbours)
% This function converts the old format of the neighbourstructure into the
% new format - although it just works as a wrapper
%
% See also FT_NEIGHBOURSELECTION

neighbours = cellStruct2StructCell(neighbours);

end


text:
function triangle_wandzura_rule_test04 ( )

%*****************************************************************************80
%
%% TRIANGLE_WANDZURA_RULE_TEST04 tests WANDZURA_RULE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 December 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGLE_WANDZURA_RULE_TEST04\n' );
  fprintf ( 1, '  WANDZURA_RULE returns the points and weights of\n' );
  fprintf ( 1, '  a Wandzura rule for the unit triangle.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This routine uses those rules to estimate the\n' );
  fprintf ( 1, '  integral of monomomials in the unit triangle.\n' );

  rule_num = wandzura_rule_num ( );

  area = 0.5;

  for a = 0 : 10

    for b = 0 : 10 - a
%
%  Multiplying X^A * Y^B by COEF will give us an integrand
%  whose integral is exactly 1.  This makes the error calculations easy.
%
      coef = ( a + b + 2 ) * ( a + b + 1 );
      for i = 1 : b
        coef = coef * ( a + i ) / i;
      end

      fprintf ( 1, '\n' );
      fprintf ( 1, '  Integrate %f * X^%d * Y^%d\n', coef, a, b );
      fprintf ( 1, '\n' );
      fprintf ( 1, '      Rule       QUAD           ERROR\n' );
      fprintf ( 1, '\n' );

      for rule = 1 : rule_num

        order_num = wandzura_order_num ( rule );

        [ xy, w ] = wandzura_rule ( rule, order_num );

        quad = 0.0;

        for order = 1 : order_num

          x = xy(1,order);
          y = xy(2,order);

          if ( a == 0 & b == 0 )
            value = coef;
          elseif ( a == 0 & b ~= 0 )
            value = coef * y^b;
          elseif ( a ~= 0 & b == 0 )
            value = coef * x^a;
          elseif ( a ~= 0 & b ~= 0 )
            value = coef * x^a * y^b;
          end

          quad = quad + w(order) * value;

        end

        quad = area * quad;

        exact = 1.0;
        err = abs ( exact - quad );

        fprintf ( 1, '  %8d  %14f  %14f\n', rule, quad, err );

      end

    end

  end

  return
end


text:
%% ---------------Truncated Amplitude Spectral Initializer-----------------

% Intializer proposed for Truncated Amplitude Flow as given in Algorithm 1
% of the Truncated Amplitude Flow (TAF) paper. For certain definitions and
% descriptions, user may need to refer to equations (14, 17) and Algorithm
% box 1 for details.
%
% The authors of the paper propose two different initialization schemes
% (this function implements the second approach, described below).
% In the first approach, one throws out the measurements of large
% magnitude, and only keeps the remaining measurement vectors that produce
% small measurements.  These vectors are nearly orthogonal to the signal,
% and so we approximate the signal by finding a vector that is
% un-correlated with these measurement vectors. This requires finding the
% smallest eigenvalue of a large matrix.
%    In the second method, the authors propose to throw out the smallest
% measurements, and find the signal most correlated with the remaining
% measurement vectors (which produce large measurement) by finding the
% largest eigenvalue of a matrix.
%    For certain isometric measurement matrices, these appraoches are both
% equivalent.  However, for non-isometric matrices they differ.  The
% authors of the paper claim that it is difficult to find a small
% eigenvector, and so they adopt the second method that requires a large
% eigenvector.  This second method, which requires the leading (largest)
% eigenvalue/vector, is implemented here.
%   Note:  in practice, small eigenvalues can be computed efficiently using
% an Arnoldi method, and this often leads to better initializers.  For an
% implementation using the smallest eigenvalue/vector, see initNull.m.
%
%  See the script 'testInitOrthogonal.m' for an example of proper usage of
%  this function.

% PAPER TITLE:
%              Solving Systems of Random Quadratic Equations via Truncated
%              Amplitude Flow.

% ARXIV LINK:
%              https://arxiv.org/pdf/1605.08285.pdf

% INPUTS:
%         A:   Function handle/numerical matrix for data matrix A. The rows
%              of this matrix are the measurement vectors that produce
%              amplitude measurements '\psi'.
%         At:  Function handle/numerical matrix for A transpose.
%         b0:  Observed data vector consisting of amplitude measurements
%              generated from b0 = |A*x|. We assign it to 'psi' to be
%              consistent with the notation in the paper.
%         n:   Length of unknown signal to be recovered by phase retrieval.

% OUPTUT :
%         x0:  The initial vector to be used by any solver.

% Note:        When a function handle is used, the value of 'n' (the length
%              of the unknown signal) and 'At' (a function handle for the
%              adjoint of 'A') must be supplied. When 'A' is numeric, the
%              values of 'At' and 'n' are ignored and inferred from the
%              arguments

% DESCRIPTION:
%              Random vectors in high dimensions are almost always
%              orthogonal to each other. Using this very intuitive result,
%              this method computes an initial guess that is orthogonal to
%              as many measurement vectors as possible. The orthogonal
%              promoting initializer solves the limitations of the spectral
%              methods, namely heavy tailed distributions due to 4th moment
%              generating functions by using a truncation step.
%              Specifically, it removes measurement vectors that are not
%              orthogonal to the initial random guess. We form the matrix Y
%              = (1/card_I) * S_bar^T* S_bar using the truncated vectors
%              and compute the leading eigenvector. S_bar is the complement
%              of the set S. Refer to the paper for a description of S.

% METHOD:
%         1.) Find the set 'I_Bar'. I is set of indices of |I| smallest
%             measurement vectors as defined in equation(14) in the paper,
%             where |I| is the cardinality of I.I_bar is the complement of
%             I. This is done because by using I, we are required to find
%             the smallest eignevalue of a matrix formed from the data
%             vectors. However, finding the smallest eigenvalue is often
%             intractible in practice. So we find the largest eigenvalue
%             instead.

%         2.) Using a mask R, form the matrix Y = (1/card_I) * S_bar^T *
%             S_bar where S_bar is the complement of the S which is defined
%             in equation (17) in the Truncated Amplitude Paper.

%         3.) Compute the leading eigenvector of Y (computed in previous
%             step) and scale it according to the norm of x as described in
%             Step 3, Algorithm 1 of the paper.
%
% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

%% -----------------------------START----------------------------------

function [x0] = initAmplitude(A,At,b0,n,verbose)

psi = b0; % To be consistent with the notation used in the paper

% If A is a matrix, infer n and At from A
if isnumeric(A)
    n = size(A, 2);
    % Transform matrix into function form
    At = @(x) A' * x;
    A = @(x) A * x;
end

m = length(psi);    % Number of measurements.

if ~exist('verbose','var') || verbose
fprintf(['Estimating signal of length %d using an orthogonal ',...
        'initializer with %d measurements...\n'],n,m);
end

% Cardinality of I. I is the set that contains the indices of the
% truncated vectors. Namely, it removes measurement vectors that are
% not orthogonal to the initial random guess
card_I = ceil(m/6);

% STEP 1: Construct the set I of indices. We approximate by assuming that
% the norm of each row of A is same.
[~, index_array] = sort(psi, 'descend');  % Sort the data vector
ind = index_array(1: card_I); % Truncate the indexes.

% STEP 2: Form Y
R = zeros(m, 1);
% Defining the mask for truncation
R(ind) = 1;
% Forming the truncated matrix Y according to equation (17) in referenced paper.
Y = @(x) 1/card_I * At(R.*A(x));

% STEP 3: Use eigs to compute leading eigenvector of Y (Y is computed in
% previous step)
opts = struct;
opts.isreal = false; % Create opts struct for eigs
[V, ~] = eigs(Y, n, 1, 'lr', opts);
% Scale the norm to match that of x
AV = abs(A(V));
alpha = (AV'*psi) / ( AV'*AV );
x0 = V*alpha;

if ~exist('verbose','var') || verbose
    fprintf('Initialization finished.\n');
end

end

author: tomgoldstein

repo: phasepack-matlab


text:
function [cMap]=turbo(varargin)

% function [cMap]=turbo(n)
% ------------------------------------------------------------------------
% Creates the colormap data for the turbo colormap.
%
% This colormap is a MATLAB implementation of the matplotlib turbo colormap
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2020/08/25
%------------------------------------------------------------------------

%%
switch nargin
    case 0
        n=250;
    case 1
        n=varargin{1};
end

%%

[cMap]=matplotlibColormap('turbo',n);

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [ y2, m2, d2, f2 ] = ymdf_next_republican ( y1, m1, d1, f1 )

%*****************************************************************************80
%
%% YMDF_NEXT_REPUBLICAN returns the Republican YMD date of the next day.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, M1, D1, real F1,
%    the YMDF date.
%
%    Output, integer Y2, M2, D2, real F2,
%    tomorrow's YMDF date.
%
  y2 = y1;
  m2 = m1;
  d2 = d1 + 1;
  f2 = f1;

  [ y2, m2, d2 ] = day_carry_republican ( y2, m2, d2 );

  return
end


text:
function E = read_buffer_offline_events(eventfile, hdr)
% function E = read_buffer_offline_events(eventfile, header)
%
% This function reads FCDC buffer-type events from a binary file.

% (C) 2010 S. Klanke

type = {
  'char'
  'uint8'
  'uint16'
  'uint32'
  'uint64'
  'int8'
  'int16'
  'int32'
  'int64'
  'single'
  'double'
};

wordsize = {
  1 % 'char'
  1 % 'uint8'
  2 % 'uint16'
  4 % 'uint32'
  8 % 'uint64'
  1 % 'int8'
  2 % 'int16'
  4 % 'int32'
  8 % 'int64'
  4 % 'single'
  8 % 'double'
};

F = fopen_or_error(eventfile,'rb',hdr.orig.endianness);
numEvt = 0;

fseek(F, 0, 'eof');
fileSize = ftell(F);
fseek(F, 0, 'bof');

maxEvt = floor(fileSize / 32);
if maxEvt == 0
   E = [];
   return;
end

E = repmat(struct('type', 0, 'sample', 0, 'duration', 0, 'offset', 0, 'value', 0), [maxEvt,1]);

num=0;

while ~feof(F)
  try
    type_type   = fread(F, 1, 'uint32');
    type_numel  = fread(F, 1, 'uint32');
    value_type  = fread(F, 1, 'uint32');
    value_numel = fread(F, 1, 'uint32');
    sample      = fread(F, 1, 'int32');
    offset      = fread(F, 1, 'int32');
    duration    = fread(F, 1, 'int32');
    bufsize     = fread(F, 1, 'uint32');
    buf         = fread(F, bufsize, 'uint8=>uint8');
  catch
    break;
  end

  type_size  = wordsize{type_type+1}*type_numel;
  value_size = wordsize{value_type+1}*value_numel;

  if type_size + value_size > bufsize
     ft_warning('Invalid event definition - skipping remaining contents');
     break;
  end

  num=num+1;
  E(num).sample = sample + 1; % offset are 1-based in Matlab
  E(num).offset = offset;
  E(num).duration = duration;
  if type_type==0
    E(num).type = char(buf(1:type_size))';
  else
    E(num).type = typecast(buf(1:type_size), type{type_type+1});
  end
  if value_type==0
    E(num).value = char(buf((type_size+1):(type_size+value_size)))';
  else
    E(num).value = typecast(buf((type_size+1):(type_size+value_size)), type{value_type+1});
  end
end

E = E(1:num);


text:
function [str, tag, cind, ccnt] = gencode_item(item, tag, tagctx, stoptag, tropts)

% function [str, tag, cind, ccnt] = gencode_item(item, tag, tagctx, stoptag, tropts)
% Generate code to recreate a cfg_(m)choice item. This code does not deal with
% arrays of cfg_items, such a configuration should not exist with the
% current definition of a configuration tree.
%
% Traversal options
% struct with fields
% stopspec - match spec to stop forced setting of eflag
% dflag    - if set to true, don't create code for .val children (code
%            for .val field is created)
% clvl     - current level in tree
% mlvl     - maximum level to force settings - range 1 (top level only) to
%            Inf (all levels)
% cnt      - item count - used for unique tags
% mcnt     - (not evaluated here)
%
% This function is identical for cfg_choice and cfg_mchoice classes.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

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

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

%% Parent object
% Generate generic object
itropts = tropts;
istoptag = stoptag;
if tropts.dflag
    % don't descend into .val tree of cfg_item
    itropts.clvl = 1;
    itropts.mlvl = 1;
    istoptag     = '';
end;
[str, tag, cind, ccnt] = gencode_item(item.cfg_item, tag, tagctx, istoptag, itropts);
tagctx = [tagctx {tag}];
% Check whether to generate code - ccnt == 0 means that generic object did
% not return code
if (tropts.clvl > tropts.mlvl || (~isempty(tropts.stopspec) && match(item, tropts.stopspec))) || ccnt == 0
    str = {};
    cind = [];
    ccnt = 0;
    return;
end;
% Reclassify generic object
str{cind} = sprintf('%s         = %s;', tag, class(item));
%% Values
% Generate values field
if numel(item.values) > 0
    % Traverse values{:} tree
    cstr = {};
    % Update clvl
    ctropts = tropts;
    ctropts.clvl = ctropts.clvl + 1;
    ctropts.cnt  = ctropts.cnt + ccnt;
    ctag = cell(size(item.values));
    for k = 1:numel(item.values)
        % tags are used as variable names and need to be unique in the
        % context of this .values tag. This includes the item's tag itself
        % and the tags of its immediate children.
        ctag{k} = genvarname(subsref(item.values{k}, substruct('.','tag')), ...
                             tagctx);
        [ccstr, ctag{k}, ccind, cccnt] = gencode_item(item.values{k}, ctag{k}, ...
                                              tagctx, stoptag, ctropts);
        tagctx = [tagctx ctag(k)];
        if ~isempty(ccstr)
            % Child has returned code
            cstr = [cstr(:)' ccstr(:)'];
            ctropts.cnt = ctropts.cnt + cccnt;
            ccnt = ccnt + cccnt;
        end;
    end;
    % Update position of class definition
    cind = cind+numel(cstr);
    % Prepend code of children
    str = [cstr(:)' str(:)'];
    str{end+1} = sprintf('%s.values  = {%s};', tag, sprintf('%s ', ctag{:}));
end;

author: spm


text:
function x = true(x)
% TRUE Constrains a variable to be positive
%
% TRUE(x) returns the constraint (x>=0.5). It is assumed that x is
% binary. The reason (x>=0.5) is used instead of x==1 is that some big-M
% modelling turns out to be less sensitive to numerical issues in this
% form. Once the model enters the solver, it will be trivially presolved
% anyway.
%
% For safety, it is advised to use the TRUE operator when working with
% logic constraints, instead of relying on the automatic constraints used
% by YALMIP (expression generated using AND and OR are automatically
% assumed to be constrained to be true.
%
% (a|b) is automatically changed to (TRUE(a|b)). To constrain a to be true,
% the user has to explicitely use (TRUE(a)).
%
%   See also SDPVAR/FALSE, SDPVAR/AND, SDPVAR/OR, SDPVAR/NOT, BINVAR, BINARY

x = (x>=.5);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [ n_data, x, fx ] = psi_values ( n_data )

%*****************************************************************************80
%
%% PSI_VALUES returns some values of the Psi or Digamma function.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      PolyGamma[x]
%
%    or
%
%      PolyGamma[0,x]
%
%    PSI(X) = d ln ( Gamma ( X ) ) / d X = Gamma'(X) / Gamma(X)
%
%    PSI(1) = -Euler's constant.
%
%    PSI(X+1) = PSI(X) + 1 / X.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 June 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 20;

  fx_vec = [ ...
    -10.42375494041108E+00, ...
     -5.289039896592188E+00, ...
     -3.502524222200133E+00, ...
     -2.561384544585116E+00, ...
     -1.963510026021423E+00, ...
     -1.540619213893190E+00, ...
     -1.220023553697935E+00, ...
     -0.9650085667061385E+00, ...
     -0.7549269499470514E+00, ...
     -0.5772156649015329E+00, ...
     -0.4237549404110768E+00, ...
     -0.2890398965921883E+00, ...
     -0.1691908888667997E+00, ...
     -0.6138454458511615E-01, ...
      0.3648997397857652E-01, ...
      0.1260474527734763E+00, ...
      0.2085478748734940E+00, ...
      0.2849914332938615E+00, ...
      0.3561841611640597E+00, ...
      0.4227843350984671E+00 ];

  x_vec = [ ...
     0.1E+00, ...
     0.2E+00, ...
     0.3E+00, ...
     0.4E+00, ...
     0.5E+00, ...
     0.6E+00, ...
     0.7E+00, ...
     0.8E+00, ...
     0.9E+00, ...
     1.0E+00, ...
     1.1E+00, ...
     1.2E+00, ...
     1.3E+00, ...
     1.4E+00, ...
     1.5E+00, ...
     1.6E+00, ...
     1.7E+00, ...
     1.8E+00, ...
     1.9E+00, ...
     2.0E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function sites = dbget_closest_sites(lon, lat, distkm, sitesdb, maxsta, snum, enum, chanmatch)
%DBGET_CLOSEST_SITES Get closest sites from an Antelop database to a given [lon, lat]
% sites = DBGET_CLOSEST_SITES(lon, lat, distkm, sitesdb, maxsta, snum, enum, chanmatch)
% load sites within distkm km of a given (LON, LAT)
% from site database sitesdb
%
% The sites structure returned has the followed fields:
%	channeltag - net.sta.loc.chan a ChannelTag object
%	longitude	- sites longitude
%	latitude	- sites latitude
%	elev	- sites elevation
%	distance	- distance (in km) from input (LON, LAT)
%
% Example: Get sites within 100km of (-150.0, 60.5) that were
% operational in 1990:
% sites = DBGET_CLOSEST_SITES(-150.0, 60.5, 100.0, '/avort/oprun/dbmaster/master_stations', 10, datenum(1990,1,1), datenum(1991,1,1));

% AUTHOR: Glenn Thompson, UAF-GI
% $Date: $
% $Revision: -1 $
debug.printfunctionstack('>')
warning(sprintf('%s assumes every station has a calibration - those without are missed',mfilename))
if ~exist('maxsta', 'var')
	maxsta = 999;
end
% if ~exist(sprintf('%s.site',sitesdb))
%     sitesdb = input('Please enter sites database path', 's')
% end
debug.print_debug(1,sprintf('sites db is %s',sitesdb))
dbptr = dbopen(sitesdb, 'r');

% Filter the site table
dbptr_site = dblookup_table(dbptr, 'site');
nrecs = dbquery(dbptr_site, 'dbRECORD_COUNT');
debug.print_debug(2,sprintf('Site table has %d records', nrecs));
dbptr_site = dbsubset(dbptr_site, sprintf('distance(lon, lat, %.4f, %.4f)<%.4f',lon,lat,km2deg(distkm)));
nrecs = dbquery(dbptr_site, 'dbRECORD_COUNT');
debug.print_debug(2,sprintf('After distance subset to lat %f and lon %f: %d records', lat, lon, nrecs));

if ~exist('snum', 'var')
    % No start time given, so assume we just want sites that exist today.
    % Remove any sites that have been decommissioned
    dbptr_site = dbsubset(dbptr_site, sprintf('offdate == NULL'));
else
    % Remove any sites that were decommissioned before the start time
    debug.print_debug(2,sprintf('offdate == NULL || offdate > %s',datenum2julday(snum)));
    dbptr_site = dbsubset(dbptr_site, sprintf('offdate == NULL || offdate > %s',datenum2julday(snum)));
end
% Remove any sites that were installed after the end time (this may remove
% some sites that exist today)
if exist('enum', 'var')
    debug.print_debug(2, sprintf('ondate  < %s',datenum2julday(enum)));
    dbptr_site = dbsubset(dbptr_site, sprintf('ondate  < %s',datenum2julday(enum)));
end
nrecs = dbquery(dbptr_site, 'dbRECORD_COUNT');
debug.print_debug(2,sprintf('After time subset: %d records', nrecs));
if nrecs==0
sites=struct();
return
end
dbgetv(dbptr_site, 'sta')

% Filter the sitechan table
dbptr_sitechan = dblookup_table(dbptr, 'sitechan');
nrecs = dbquery(dbptr_sitechan, 'dbRECORD_COUNT');
debug.print_debug(2,sprintf('sitechan has %d records', nrecs));

if exist('chanmatch','var')
    dbptr_sitechan = dbsubset(dbptr_sitechan, chanmatch);
    nrecs = dbquery(dbptr_sitechan, 'dbRECORD_COUNT');
    debug.print_debug(2,sprintf('After chan subset: %d records', nrecs));
end

if ~exist('snum', 'var')
    % No start time given, so assume we just want sites that exist today.
    % Remove any sites that have been decommissioned
    dbptr_sitechan = dbsubset(dbptr_sitechan, sprintf('offdate == NULL'));
else
    % Remove any sites that were decommissioned before the start time
    debug.print_debug(2,sprintf('offdate == NULL || offdate > %s',datenum2julday(snum)));
    dbptr_sitechan = dbsubset(dbptr_sitechan, sprintf('offdate == NULL || offdate > %s',datenum2julday(snum)));
end
% Remove any sites that were installed after the end time (this may remove
% some sites that exist today)
if exist('enum', 'var')
    debug.print_debug(2,sprintf('ondate  < %s',datenum2julday(enum)));
    dbptr_sitechan = dbsubset(dbptr_sitechan, sprintf('ondate  < %s',datenum2julday(enum)));
end
nrecs = dbquery(dbptr_sitechan, 'dbRECORD_COUNT');
debug.print_debug(2,sprintf('After time subset: %d records', nrecs));
dbgetv(dbptr_sitechan, 'sta')

% Join site and sitechan
dbptr_sitechan = dbjoin(dbptr_site, dbptr_sitechan);
nrecs = dbquery(dbptr_sitechan, 'dbRECORD_COUNT');
debug.print_debug(2,sprintf('After join site-sitechan %d records', nrecs));
dbgetv(dbptr_sitechan, 'sta')

% Read snetsta to get network code
dbptr_snetsta = dblookup_table(dbptr, 'snetsta');
nrecs = dbquery(dbptr_snetsta, 'dbRECORD_COUNT');
debug.print_debug(2,sprintf('snetsta has %d records', nrecs));
dbptr_snetsta = dbjoin(dbptr_sitechan, dbptr_snetsta);
nrecs = dbquery(dbptr_snetsta, 'dbRECORD_COUNT');
debug.print_debug(2,sprintf('After join site-sitechan-snetsta: %d records', nrecs));
snet = dbgetv(dbptr_snetsta, 'snet');
sta = dbgetv(dbptr_snetsta, 'sta');
stanetmap = containers.Map();
for c=1:numel(sta)
    stanetmap(sta{c}) = snet{c};
end

% If there is a calibration table, get the calib value and units
dbptr_calibration = dblookup_table(dbptr, 'calibration');
nrecs = dbquery(dbptr_calibration, 'dbRECORD_COUNT');
debug.print_debug(2,sprintf('calibration has %d records', nrecs));
%dbptr_calibration = dbjoin(dbptr_sitechan, dbptr_calibration,
%{'sta';'chan'}, {'sta';'chan'}); % danger with this is no subset on ondate/offdate
dbptr_calibration = dbjoin(dbptr_sitechan, dbptr_calibration);
dbptr_calibration = dbsubset(dbptr_calibration, 'sitechan.sta == calibration.sta && sitechan.chan == calibration.chan');
nrecs = dbquery(dbptr_calibration, 'dbRECORD_COUNT');
debug.print_debug(2,sprintf('After join site-sitechan-calibration: %d records', nrecs));

% Read data from final view
if nrecs == 0
    sites = [];
    return
end
% net = dbgetv(db4, 'snet');
% if ~iscell(net)
%     net = {net};
% end

latitude = dbgetv(dbptr_calibration, 'lat');
longitude = dbgetv(dbptr_calibration, 'lon');
elev = dbgetv(dbptr_calibration, 'elev');
staname = dbgetv(dbptr_calibration, 'sta');
if ~iscell(staname)
	staname = {staname};
end
channame = dbgetv(dbptr_calibration, 'chan');
if ~iscell(channame)
	channame = {channame};
end
ondate = dbgetv(dbptr_calibration, 'sitechan.ondate');
offdate = dbgetv(dbptr_calibration, 'sitechan.offdate');
calib = dbgetv(dbptr_calibration, 'calibration.calib');
units = dbgetv(dbptr_calibration, 'calibration.units');
dbclose(dbptr);

% Reformat data into sites return structure
numsites = numel(latitude);
for c=1:numsites
    yyyy = floor(ondate(c)/1000);
    if yyyy>1900
        jjj = ondate(c)-yyyy*1000;
        ondnum(c) = datenum(yyyy,1,jjj);
    else
        ondnum(c)=-Inf;
    end
end

for c=1:numsites
    yyyy = floor(offdate(c)/1000);
    if yyyy>1900
        jjj = offdate(c)-yyyy*1000;
        offdnum(c) = datenum(yyyy,1,jjj);
    else
       offdnum(c) = Inf;
    end
end

for c=1:numsites
    stadist(c) = deg2km(distance(lat, lon, latitude(c), longitude(c)));
end

% order the sites by distance
[y,i]=sort(stadist);
c=1;
while ((c<=numsites) && (stadist(i(c)) < distkm))
	%sites(c).name = staname{i(c)};
	%sites(c).channel = channame{i(c)};
	%sites(c).scnl = scnlobject(sites(c).name, sites(c).channel, net{i(c)});
    sta = staname{i(c)};
    chan = channame{i(c)};
    net = stanetmap(sta);
    sites(c).channeltag = ChannelTag(net, sta, '', chan);
    sites(c).string = string(sites(c).channeltag);
	sites(c).longitude = longitude(i(c));
	sites(c).latitude = latitude(i(c));
	sites(c).elev = elev(i(c));
	sites(c).distance = stadist(i(c));
    sites(c).ondnum = ondnum(i(c));
    sites(c).offdnum = offdnum(i(c));
    sites(c).calib = calib(i(c));
    sites(c).units = units(i(c));
	c = c + 1;
end

% remove any duplicate sites
[~,j]=unique({sites.string});
sites = sites(j);

% limit the number of sites
numsites = min([maxsta numel(sites)]);
sites = sites(1:numsites);

debug.printfunctionstack('<')

author: geoscience-community-codes

repo: GISMO


text:
% laplac2d() - generate a 2 dimensional gaussian matrice
%
% Usage :
%    >> [ gaussmatrix ] = laplac2d( rows, columns, sigma, ...
%                                       meanR, meanC, cut)
%
% Example :
%   >> laplac2d( 5, 5)
%
% Inputs:
%   rows    - number of rows
%   columns - number of columns
%   sigma   - standart deviation (default: rows/5)
%   meanR   - mean for rows (default: center of the row)
%   meanC   - mean for columns (default: center of the column)
%   cut	    - percentage (0->1) of the maximum value for removing
%             values in the matrix (default: 0)
%
% Note: this function implements a simple laplacian for exploratory
%       research. For a more rigorous validated approach use the freely
%       available Current Source Density Matlab toolbox.
%
% See also: eeg_laplac()
%
% Author: Arnaud Delorme, CNL, Salk Institute, 2001

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

function mat = laplac2d( sizeX, sizeY, sigma, meanX, meanY, cut);

if nargin < 2
	help laplac2d
	return;
end;
if nargin < 3
	sigma = sizeX/5;
end;
if nargin < 4
	meanX = (sizeX+1)/2;
end;
if nargin < 5
	meanY = (sizeY+1)/2;
end;
if nargin < 6
	cut = 0;
end;

X = linspace(1, sizeX, sizeX)'* ones(1,sizeY);
Y = ones(1,sizeX)'   		  * linspace(1, sizeY, sizeY);
%[-sizeX/2:sizeX/2]'*ones(1,sizeX+1);
%Y = ones(1,sizeY+1)'   *[-sizeY/2:sizeY/2];

r2 = (X-meanX).*(X-meanX) + (Y-meanY).*(Y-meanY);
sigma2 = sigma*sigma;

mat = - exp(-0.5*r2/sigma2) .* ((r2 - sigma2)/(sigma2*sigma2));
% zeros crossing at r = -/+ sigma;
% mat = r2;
if cut > 0
	maximun = max(max(mat))*cut;
	I = find(mat < maximun);
	mat(I) = 0;
end;

return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function c = correlation_cubic ( n, rho, rho0 )

%*****************************************************************************80
%
%% CORRELATION_CUBIC evaluates the cubic correlation function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Petter Abrahamsen,
%    A Review of Gaussian Random Fields and Correlation Functions,
%    Norwegian Computing Center, 1997.
%
%  Parameters:
%
%    Input, integer N, the number of arguments.
%
%    Input, real RHO(N,1), the arguments.
%
%    Input, real RHO0, the correlation length.
%
%    Output, real C(N,1), the correlations.
%
  rho = rho ( : );

  rhohat = min ( abs ( rho ) / rho0, 1.0 );

  c = 1.0 - 7.0 * rhohat.^2 + 8.75 * rhohat.^3 - 3.5 * rhohat.^5 + 0.75 * rhohat.^7;

  return
end


text:
function v = indexBy(idx, arr)
% VAL = INDEXBY(INDEX, ARRAY)
%   Get the INDEX'th of ARRAY regardless of if it is a cell or not.
%

  if iscell(arr)
    v = arr{idx};
  else
    v = arr(idx);
  end
end

author: Sable


text:
function [wstilda wstilda2,wstilda3]=fbssmoother(a,b,xf,wftilda)
[n,T]=size(xf);
%Initialize weights
wstilda(:,T)=wftilda(:,T);
for i=1:n
ws3(i,T)=wftilda(i,T)*weight(xf(i,T),xf(:,T),wstilda(i,T),wftilda(:,T),i,a,b);
end
wstilda3(:,T)=ws3(:,T)/sum(ws3(:,T));

for t=T-1:-1:1
%Compute smoother weights p(x(t))|Y)
for i=1:n
ws(i,t)=wftilda(i,t)*weight(xf(:,t+1),xf(:,t),wstilda(:,t+1),wftilda(:,t),i,a,b);
end
wstilda(:,t)=ws(:,t)/sum(ws(:,t));
%Compute joint smoother weight p(x(t),x(t+1)|Y) and p(x(t),x(t)|Y)
for i=1:n
ws2(i,t)=wftilda(i,t)*weight(xf(i,t+1),xf(:,t),wstilda(i,t+1),wftilda(:,t),i,a,b);
ws3(i,t)=wftilda(i,t)*weight(xf(i,t),xf(:,t),wstilda(i,t),wftilda(:,t),i,a,b);
end
wstilda2(:,t)=ws2(:,t)/sum(ws2(:,t));
wstilda3(:,t)=ws3(:,t)/sum(ws3(:,t));
end