text:
function Population = subRVEA(Population,V,theta)
% The environmental selection of RVEA

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Cheng He

    PopObj = Population.objs;
    [N,M]  = size(PopObj);
    NV     = size(V,1);

    %% Translate the population
    PopObj = PopObj - repmat(min(PopObj,[],1),N,1);

    %% Calculate the degree of violation of each solution
    CV = sum(max(0,Population.cons),2);

    %% Calculate the smallest angle value between each vector and others
    cosine = 1 - pdist2(V,V,'cosine');
    cosine(logical(eye(length(cosine)))) = 0;
    gamma  = min(acos(cosine),[],2);

    %% Associate each solution to a reference vector
    Angle = acos(1-pdist2(PopObj,V,'cosine'));
    [~,associate] = min(Angle,[],2);

    %% Select one solution for each reference vector
    Next = zeros(1,NV);
    for i = unique(associate)'
        current1 = find(associate==i & CV==0);
        current2 = find(associate==i & CV~=0);
        if ~isempty(current1)
            % Calculate the APD value of each solution
            APD = (1+M*theta*Angle(current1,i)/gamma(i)).*sqrt(sum(PopObj(current1,:).^2,2));
            % Select the one with the minimum APD value
            [~,best] = min(APD);
            Next(i)  = current1(best);
        elseif ~isempty(current2)
            % Select the one with the minimum CV value
            [~,best] = min(CV(current2));
            Next(i)  = current2(best);
        end
    end
    % Population for next generation
    Population = Population(Next(Next~=0));
end


text:
function score = mutual_info_score(i,si,j,sj,data)
% G = mutual_info_score(i,si,j,sj,data)
% Only for tabular node which values are 1,2,...,size .
% si is size of node i, sj is size of node j.
% data(i,m) is the node i in the case m.
%
% Ref :
% C. Chow and C. Liu (1968). Approximating discrete probability distributions with dependence trees.
% IEEE Transactions on Information Theory, 14(3):462--467, May 1968.
%
% francois.olivier.c.h@gmail.com, philippe.leray@univ-nantes.fr, wangxiangyang@sjtu.edu.cn

[n N]=size(data);
Nj=hist(data(j,:),1:sj);
Ni=hist(data(i,:),1:si);
NiNj=Ni'*Nj;

for k=1:si
 ind=find(data(i,:)==k) ;
 Nij(k,:) = hist(data(j,ind),1:sj);
end

% sommons les valeurs non-infinies:
ind=find(NiNj~=0 & Nij~=0);
score=sum(sum(Nij(ind).*log(N*Nij(ind)./NiNj(ind))/N));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


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

    properties
        nonSeqEvts LaunchVehicleNonSeqEvents
        event LaunchVehicleEvent
    end

    methods
        function obj = NonSeqEvtsState(event, nonSeqEvts)
            obj.event = event;
            obj.nonSeqEvts = nonSeqEvts;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef crossvalidator
% CROSSVALIDATOR crossvalidation class.
%
%   DESCRIPTION
%   The crossvalidator class performs crossvalidation. It uses a
%   particular crossvalidation on a multivariate analysis specified in the
%   property mva.
%
%   One may use 'type' equal to:
%   - 'nfold' : n-fold cross-validation specified by 'folds' (default: 5)
%   - 'split' : data is split into train and test data; amount of training
%               data is given by 'proportion' (default: 0.75)
%   - 'loo'   : leave-one-out cross-validation; takes 1 trial per fold
%   - 'bloo'  : balanced loo cv; takes 1 trial from each class per fold.
%               Assumes that each class has the same nr of trials.
%
%   Alternatively one can specify the cell-arrays 'trainfolds' and
%   'testfolds' to contain for each fold the trial indices that need to be
%   used. If one of these is specified then the other one is automatically
%   filled using the complement of the trials. If these cell-arrays are
%   non-empty then the settings above are ignored.
%
%   In order to balance the occurrence of different classes one may set
%   'resample' equal to true (default: false). Resample will upsample less
%   occurring classes during training and downsample often occurring
%   classes during testing.
%
%   In case of memory problems, use 'compact' equal to true (default:
%   false). This will not store the trained multivariate methods but just
%   the classification outcomes and the parameters that are returned by the
%   model function of the multivariate methods.
%
%   EXAMPLE
%   X = rand(10,20); Y = [1 1 1 1 1 2 2 2 2 2]';
%   m = dml.crossvalidator('mva',dml.svm)
%   m = m.train(X,Y);
%   m.statistic('accuracy')
%
%   DEVELOPER
%   Marcel van Gerven (m.vangerven@donders.ru.nl)

    properties

      mva % multivariate analysis

      result % crossvalidation result

      design % the real outputs

      model % estimated model per fold

      type = 'nfold'; % 'type' : crossvalidation type 'nfold', 'split', 'loo'

      folds = 5; % 'folds' : number of folds when using nfold

      proportion = 0.75; % 'proportion' : proportion of used training data when using split

      resample = false; % 'resample' : whether or not to use resampling to balance classes

      % cell array indicating used trials per training fold
      % for multiple datasets this should be come a cell array of cell
      % arrays, e.g. trainfolds = { {f1 ... fn} {f1 ... fn} }
      trainfolds

      % cell array indicating used trials per test fold
      % for multiple datasets this should be come a cell array of cell
      % arrays, e.g. testfolds = { {f1 ... fn} {f1 ... fn} }
      testfolds

      stat = ''; % test statistic used to quantify performance of this cross-validator

      verbose = false; % generate output when true

      compact = false; % drop mva details when true

    end

    methods

      function obj = crossvalidator(varargin)

        % parse options
        for i=1:2:length(varargin)
          if ismember(varargin{i},fieldnames(obj))
            obj.(varargin{i}) = varargin{i+1};
          else
            error('unrecognized fieldname %s',varargin{i});
          end
        end

        if ~isa(obj.mva,'dml.analysis'), obj.mva = dml.analysis(obj.mva); end

      end

      function obj = train(obj,X,Y)

        % complete the folds when only train or test is specified
        if isempty(obj.trainfolds) && isempty(obj.testfolds)
          [obj.trainfolds,obj.testfolds] = obj.create_folds(Y);
        elseif isempty(obj.trainfolds)
          obj.trainfolds = obj.complement(Y,obj.testfolds);
        else
          obj.testfolds = obj.complement(Y,obj.trainfolds);
        end

        if iscell(Y)
          ndata = length(Y);
        else
          ndata = 1;
        end

        if ndata == 1
          nfolds = length(obj.trainfolds);
        else
          nfolds = length(obj.trainfolds{1});
        end

        obj.result = cell(nfolds,1);
        obj.design = cell(nfolds,1);
        obj.model = cell(nfolds,1);

        if ~iscell(obj.mva)
          obj.mva = repmat({obj.mva},[1 nfolds]);
        end

        for f=1:nfolds % iterate over folds

          if obj.verbose
            if ndata == 1
              fprintf('validating fold %d of %d for %d datasets\n',f,nfolds,ndata);
            else
              fprintf('validating fold %d of %d for %d datasets\n',f,nfolds,ndata);
            end
          end

          % construct X and Y for each fold
          if ndata == 1
            trainX = X(obj.trainfolds{f},:);
            testX = X(obj.testfolds{f},:);
            trainY = Y(obj.trainfolds{f},:);
            testY = Y(obj.testfolds{f},:);
          else
            trainX = cell([length(X) 1]);
            testX = cell([length(X) 1]);
            for c=1:length(X)
              trainX{c} = X{c}(obj.trainfolds{c}{f},:);
              testX{c} = X{c}(obj.testfolds{c}{f},:);
            end
            trainY = cell([length(Y) 1]);
            testY = cell([length(Y) 1]);
            for c=1:length(Y)
              trainY{c} = Y{c}(obj.trainfolds{c}{f},:);
              testY{c} = Y{c}(obj.testfolds{c}{f},:);
            end
          end

          tproc = obj.mva{f}.train(trainX,trainY);
          if ~isempty(testY)
            obj.result{f} = tproc.test(testX);
            obj.design{f} = testY;
          end
          obj.model{f} = tproc.model();

          if ~obj.compact, obj.mva{f} = tproc; end

          clear tproc;
          clear trainX;
          clear testX;
          clear trainY;
          clear testY;

        end

        % return unique model instead of cell array in case of one fold
        if length(obj.model)==1, obj.model = obj.model{1}; end

      end

      function [train,test] = create_folds(obj,Y)
      % get predefined sequence for train and test folds

      if strcmp(obj.type,'loo') || strcmp(obj.type,'bloo')

        if iscell(Y)
          obj.folds = min(cellfun(@(x)(size(x,1)),Y));
          if strcmp(obj.type,'bloo')
            nclasses = max(cellfun(@(x)(max(x)),Y));
            obj.folds = obj.folds/nclasses;
          end
        else
          obj.folds = size(Y,1);
          if strcmp(obj.type,'bloo')
            nclasses = max(Y);
            for k = 1:nclasses
              N(k,1) = sum(Y==k);
            end
            obj.folds = min(N);
            %obj.folds = obj.folds/nclasses;
          end
        end
        obj.type = 'nfold';

      end

      if strcmp(obj.type,'split'), obj.folds = 1; end

      if ~iscell(Y)

        test = create_test_folds(Y);

        train = obj.complement(Y,test);

        if obj.resample
          train = upsample(train,Y);
          test  = downsample(test,Y);
        end

      else

        train = cell(length(Y),1);
        test = cell(length(Y),1);

        for c=1:length(Y)

          test{c} = create_test_folds(Y{c});

          train{c} = obj.complement(Y{c},test{c});

          if obj.resample
            train{c} = upsample(train{c},Y{c});
            test{c}  = downsample(test{c},Y{c});
          end

        end

      end

        function y = create_test_folds(Y)

          % NOTE: Fixing of the random seed has been removed here in the fieldtrip
          % version since the randomseed is dealt with higher up in the code
          % hierarchy (i.e. specify cfg.randomseed in ft_timelockstatistics).

          y = cell(obj.folds,1);

          % determine indices of labeled (non-nan) datapoints
          % nan datapoints should end up in the training set
          labeled = find(any(~isnan(Y(1:size(Y,1),:)),2));
          if obj.verbose && numel(labeled)<size(Y,1)
            fprintf('adding unlabeled samples to training set\n');
          end

          % randomize labeled trials
          nsamples = size(labeled,1);
          idxs = labeled(randperm(nsamples));

          % make sure outcomes are evenly represented whenever possible
          [t,t,idx] = unique(Y,'rows');
          mx = max(idx);

          switch obj.type

            case 'split'

              if obj.verbose, fprintf('creating sample indices using %.0f%% of the data for testing\n',100*(1-obj.proportion)); end

              if mx == nsamples % unique labeled samples

                y{1} = idxs(1:round((1-obj.proportion)*nsamples));

              else

                % take labeled indices
                idx = idx(idxs);
                for j=1:mx
                  iidx = find(idx == j);
                  y{1} = [y{1}; idxs(iidx(1:floor((1-obj.proportion)*numel(iidx))))];
                end

                % add random instances to complete the number
                n = floor((1-obj.proportion) * numel(idxs)) - numel(y{1});
                if n>0
                  x = setdiff(idxs,y{1});
                  x = x(randperm(numel(x)));
                  y{1} = [y{1}; x(1:n)];
                end

              end

            case 'nfold'

              if obj.verbose, fprintf('creating sample indices using %d-fold cross-validation\n',obj.folds); end

              if mx == nsamples % unique samples

                for f=1:obj.folds
                  y{f} = idxs((floor((f-1)*(length(idxs)/obj.folds))+1):floor(f*(length(idxs)/obj.folds)));
                end

              else

                % take labeled indices
                idx = idx(idxs);

                f=1;
                for j=1:mx
                  iidx = find(idx == j);
                  for k=1:length(iidx)
                    y{f} = [y{f}; idxs(iidx(k))];
                    f = f+1; if f > obj.folds, f=1; end
                  end
                end

              end

            otherwise
              error('unrecognized option for type');
          end

        end

        function y = downsample(F,Y)
          % make sure all classes are represented equally in the
          % train samples by sampling with replacement and in the
          % test samples by using the minimum number of trials in a
          % condition. This ensures that training makes most use of
          % the data while testing has no inherent bias in the
          % data. We lose power by removing superfluous trials.

          y = cell(size(F));

          if obj.verbose
            fprintf('balancing test samples by removing superfluous trials\n');
          end

          for f=1:length(F)

            % get unique rows of Y
            [unq,t,idx] = unique(Y,'rows');

            idx = idx(F{f});

            minsmp = min(arrayfun(@(x)(sum(idx == x)),unq));

            if minsmp > 0

              tmp = [];
              for j=1:length(unq)

                iidx = find(idx == unq(j));

                if obj.verbose && length(iidx) - minsmp ~=0
                  fprintf('removing %d superfluous samples for label %d\n',length(iidx) - minsmp,unq(j));
                end

                tmp = [tmp; F{f}(iidx(1:minsmp))];

              end

              y{f} = tmp(randperm(length(tmp)));

            end

          end

        end

        function y = upsample(F,Y)
          % make sure all classes are represented equally in the
          % train samples by sampling with replacement and in the
          % test samples by using the minimum number of trials in a
          % condition. This ensures that training makes most use of
          % the data while testing has no inherent bias in the
          % data. We lose power by removing superfluous trials.

          y = cell(size(F));

          if obj.verbose
            fprintf('balancing training samples by sampling with replacement\n');
          end

          for f=1:length(F)

            % get unique rows of Y
            [unq,t,idx] = unique(Y,'rows');

            idx = idx(F{f});

            maxsmp = max(arrayfun(@(x)(sum(idx == x)),unq));

            tmp = F{f};
            for j=1:length(unq)

              iidx = find(idx == unq(j));

              if obj.verbose && maxsmp - length(iidx) ~= 0
                fprintf('drawing %d additional samples for label %d\n',maxsmp - length(iidx),unq(j));
              end

              % nameclash with fieldtrip_private!
              tmp = [tmp; randsample(F{f}(iidx),maxsmp - length(iidx),true)];

            end

            y{f} = tmp(randperm(length(tmp)));

          end

        end

      end

      function s = statistic(obj,stat)
        % report statistics on a trained crossvalidator object; also see
        % dml.statistic

        if nargin<2, stat = obj.stat; end

        if iscell(obj.design{1}) && iscell(obj.result{1})

          s = zeros(length(obj.design),1);
          for c=1:length(obj.design{1}) % iterate over datasets
            D = [];
            P = [];
            for f=1:length(obj.design) % iterate over folds
              D = cat(1,D,obj.design{f}{c});
              P = cat(1,P,obj.result{f}{c});
            end
            s(c) = dml.statistic(stat,D,P);
          end

        else

          D = cell2mat(obj.design);
          P = cell2mat(obj.result);
          s = dml.statistic(stat,D,P);

        end

      end

    end

    methods(Access=private,Static=true)

      function y = complement(Y,folds)
        % create complement of trials in each of the folds when only training
        % or test folds are specified

        if iscell(Y)
          ndata = length(Y);
        else
          ndata = 1;
        end

        y = cell(size(folds));
        for c=1:length(folds)
          if ndata == 1
            y{c} = setdiff(1:size(Y,1),folds{c})';
          else
            for f=1:length(folds{c})
              y{c}{f} = setdiff(1:size(Y{c},1),folds{c}{f})';
            end
          end
        end

      end

    end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef (Abstract) TrackInitiatorX < BaseX
% TrackInitiatorX Abstract class
%
% Summary of TrackInitiatorX:
% This is the base class for all TrackingX resamplers.
% Any custom defined Track Initiator should be derived from this TrackInitiatorX base class.
%
% TrackInitiatorX Properties:
%   None
%
% TrackInitiatorX Methods:
%   + TrackInitiatorX - Constructor method
%
% (+) denotes public properties/methods
%
% February 2018 Lyudmil Vladimirov, University of Liverpool.

    properties
    end

    methods
        function this = TrackInitiatorX(varargin)
        % TRACKINITIATORX Constructor method
        %
        % DESCRIPTION:
        % * TrackInitiatorX() returns a "TrackInitiatorX" object handle

            this@BaseX();
        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function params = setParams(params)

% HW params
params.Fs = 400;
% Number of channels
params.nChans = 16;

% Subjects to include
params.trainSubs = {'1', '2', '3'};
nSubs = numel(params.trainSubs);
params.testSubs = {'1', '2', '3'};

% OK Check params
% Marks files as OK in fileList if proportion of zeros is below this
% threshold
% Not used anymore
params.OKThresh = 0.5;

% Feature extraction params
params.HillsBands.Range = 1:47;

author: garethjns

repo: Kaggle-EEG


text:
function [Dict Coeff]=KSVD_Inpainting(Dict,blkMatrixIm,blkMask,sigma,rc_min,max_coeff,max_iter)
% Author: Chongyang
% Date: 2016/06/05
% KSVD for Inpainting
%
% sigma: residual error stopping criterion, normalized by signal norm  0.01
% rc_min: minimal residual correlation before stopping 0.01
% max_coeff: maximal number of non-zero coefficients for signal 10t

fprintf('KSVD begin .....\n')
tic
for KSVDiter=1:max_iter
    fprintf('iteration: %d/%d | time: %d \n',KSVDiter,max_iter,toc)
 	Coeff=OMP_Inpainting(Dict,blkMatrixIm.*blkMask,blkMask,sigma,rc_min,max_coeff); %256*146405
    for atom=1:1:size(Dict,2)
        Omega=find(abs(Coeff(atom,:))>0);
         % this atom will be useless
        if isempty(Omega)
            continue;
        end;
        CoeffM=Coeff;
        CoeffM(atom,:)=0;
        Err=blkMask(:,Omega).*(blkMatrixIm(:,Omega)-Dict*CoeffM(:,Omega));

        dd=diag(1./(blkMask(:,Omega)*(Coeff(atom,Omega).^2)'))*(Err*Coeff(atom,Omega)');
        Dict(:,atom)=dd/norm(dd);
        Coeff(atom,Omega)=(dd'*Err)./sum((blkMask(:,Omega).*(dd*ones(1,length(Omega)))).^2,1);

    end;
%     save Dict Dict;
%     save Coeff Coeff;
end;

author: chongyangtao

repo: Color-Image-Inpainting


text:
classdef TankMassConstraint < AbstractConstraint
    %TankMassConstraint Summary of this class goes here
    %   Detailed explanation goes here

    properties
        normFact = 1;
        tank LaunchVehicleTank
        event LaunchVehicleEvent
        eventNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        evalType(1,1) ConstraintEvalTypeEnum = ConstraintEvalTypeEnum.FixedBounds;
        stateCompType(1,1) ConstraintStateComparisonTypeEnum = ConstraintStateComparisonTypeEnum.Equals;
        stateCompEvent LaunchVehicleEvent
        stateCompNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;
    end

    methods
        function obj = TankMassConstraint(tank, event, lb, ub)
            obj.tank = tank;
            obj.event = event;
            obj.lb = lb;
            obj.ub = ub;

             obj.id = rand();
        end

        function [lb, ub] = getBounds(obj)
            lb = obj.lb;
            ub = obj.ub;
        end

        function [c, ceq, value, lwrBnd, uprBnd, type, eventNum, valueStateComp] = evalConstraint(obj, stateLog, celBodyData)
            type = obj.getConstraintType();

            switch obj.eventNode
                case ConstraintStateComparisonNodeEnum.FinalState
                    stateLogEntry = stateLog.getLastStateLogForEvent(obj.event);

                case ConstraintStateComparisonNodeEnum.InitialState
                    stateLogEntry = stateLog.getFirstStateLogForEvent(obj.event);

                otherwise
                    error('Unknown event node.');
            end

            value = lvd_TankMassTasks(stateLogEntry, 'tankMass', obj.tank);

            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                switch obj.stateCompNode
                    case ConstraintStateComparisonNodeEnum.FinalState
                        stateLogEntryStateComp = stateLog.getLastStateLogForEvent(obj.stateCompEvent);

                    case ConstraintStateComparisonNodeEnum.InitialState
                        stateLogEntryStateComp = stateLog.getFirstStateLogForEvent(obj.stateCompEvent);

                    otherwise
                        error('Unknown event node.');
                end

                valueStateComp = lvd_TankMassTasks(stateLogEntryStateComp, 'tankMass', obj.tank);
            else
                valueStateComp = NaN;
            end

            [c, ceq] = obj.computeCAndCeqValues(value, valueStateComp);

            lwrBnd = obj.lb;
            uprBnd = obj.ub;

            eventNum = obj.event.getEventNum();
        end

        function sF = getScaleFactor(obj)
            sF = obj.normFact;
        end

        function setScaleFactor(obj, sF)
            obj.normFact = sF;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = obj.tank == tank;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesEvent(obj, event)
            tf = obj.event == event;
            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                tf = tf || obj.stateCompEvent == event;
            end
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end

        function tf = usesExtremum(obj, extremum)
            tf = false;
        end

        function tf = usesGroundObj(obj, grdObj)
            tf = false;
        end

        function tf = canUseSparseOutput(obj)
            tf = true;
        end

        function event = getConstraintEvent(obj)
            event = obj.event;
        end

        function type = getConstraintType(obj)
            type = 'Tank Mass';
        end

        function [unit, lbLim, ubLim, usesLbUb, usesCelBody, usesRefSc] = getConstraintStaticDetails(obj)
            unit = 'mT';
            lbLim = 0;
            ubLim = Inf;
            usesLbUb = true;
            usesCelBody = false;
            usesRefSc = false;
        end

        function addConstraintTf = openEditConstraintUI(obj, lvdData)
            [~, tanks] = lvdData.launchVehicle.getTanksListBoxStr();
            if(numel(tanks) >= 1)
%                 addConstraintTf = lvd_EditTankConstraintGUI(obj, lvdData);

                output = AppDesignerGUIOutput({false});
                lvd_EditTankConstraintGUI_App(obj, lvdData, output);
                addConstraintTf = output.output{1};
            else
                errordlg('There are currently no tanks assigned to the launch vehicle in this scenario.  Add at least one tank first.');

                addConstraintTf = false;
            end
        end
    end

    methods(Static)
        function constraint = getDefaultConstraint(~, ~)
            constraint = TankMassConstraint(LaunchVehicleTank.empty(1,0), LaunchVehicleEvent.empty(1,0), 0, 0);
        end
    end
end


text:
function test_ft_mvaranalysis

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_mvaranalysis

fs = 500;
nchan = 32;
start_time = -1; %s
end_time = 2.5; %s
nsamples = (end_time - start_time) * fs + 1;

data = [];
data.time{1} = linspace(start_time, end_time, nsamples);
data.trial{1} = randn(nchan,nsamples);
data.label = cellstr(num2str((1:nchan).'));

cfg = [];
dataout = ft_mvaranalysis(cfg,data);


text:
function test_failed=test_filterbankscale()
test_failed = 0;

disp('----------FILTERBANKSCALE-----------');

 gr{1}=randn(20,1);
 gr{2}=randn(21,1);
 gr{3}=firfilter('hanning',19);
 gr{4}=firfilter('hanning',20);
 gr{5}=randn(4,1);
 gr{6}=firfilter('hanning',20,'causal');
 gr{7}=firfilter('hanning',20,'delay',13);
 gr{8}=firfilter('hanning',20,'delay',-13);
 gr{7}=firfilter('hanning',20,'delay',14);
 gr{8}=firfilter('hanning',20,'delay',-14);
 gr{9}=firfilter('hamming',19,.3);
 gr{10}=firfilter('hamming',19,.3,'real');
 gr{11}=blfilter('hanning',.19);
 gr{12}=blfilter('hanning',.2);
 gr{13}=blfilter('hanning',.132304);
 gr{14}=blfilter('hanning',.23,'delay',13);
 gr{15}=blfilter('hamming',.23,.3);
 gr{16}=blfilter('hamming',.23,.3,'real');
 % (almost) Allpass filter
 gr{17}=blfilter('hanning',2);
 gr{18}=blfilter('hanning',1);
 gr{19}=blfilter('hanning',2,1);
 gr{20}=blfilter('hanning',1,-1);
 gr{21}=blfilter('hamming',.1,-.3);
 gr{22}=blfilter('hanning',1.7);
 gr{23}=struct('H',randn(10,1),'L',100);
 gr{24}=crand(40,1);
 gr{25}=struct('h',crand(40,1));
 gr{26}=struct('h',crand(40,1),'realonly',1);

 Larr = [100, 211];
 norms = {'1','2','inf'};

 for L = Larr
     gr{23}=struct('H',randn(10,1),'L',L);
     for ii = 0:1
         if ii==0, freqstr = ''; freqflag = 'nofreq';
         else freqstr = 'FREQ'; freqflag = 'freq'; end

         for nId = 1:numel(norms)
             g = filterbankscale(gr,L,norms{nId},freqflag);
             if ii == 0
                 [~,n] = setnorm(ifft(filterbankfreqz(g,1,L)),norms{nId});
             else
                [~,n] = setnorm(filterbankfreqz(g,1,L),norms{nId});
             end

             res = sum(abs(n-1));
             [test_failed,fail]=ltfatdiditfail(res,test_failed);
             fprintf('NORM: %s %s L=%d %s\n',upper(norms{nId}),freqstr,L,fail);
         end
     end
 end

 scal = 0.1;
  for L = Larr
     gr{23}=struct('H',randn(10,1),'L',L);

         for nId = 1:numel(norms)
             grfreqz = filterbankfreqz(gr,1,L);
             g = filterbankscale(gr,scal);

             gfreqz = filterbankfreqz(g,1,L);

             res = sum(sum(scal*grfreqz-gfreqz));
             [test_failed,fail]=ltfatdiditfail(res,test_failed);
             fprintf('SCALE %f L=%d %s\n',scal,L,fail);

         end
 end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [hdl, inf] = dxAviOpen(fname);
	hdl = [];
	inf = [];
	if ~exist(fname,'file')
		fprintf(2,'dxAviOpen error, cannot find file %s\n', fname);
		return;
	end
	[hdl, t] = dxAviOpenMex(fname);
	if length(hdl) == 0;
		error(['Error opening video file ', fname]);
		return;
	end
	inf.Width = t(1);
	inf.Height = t(2);
	inf.NumFrames = t(3);
	inf.FrameRate = t(4);
	inf.TotalDuration = t(5);

author: Sable

repo: mcbench-benchmarks


text:
function test_bug2784

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_mvaranalysis

data = [];
data.trial{1} = randn(2,100);
data.trial{2} = randn(2,101);
data.time{1}  = (0:99)./100+0.001;
data.time{2}  = (0:100)./100;
data.label    = {'chan01';'chan02'};

try
  cfg = [];
  ft_mvaranalysis(cfg,data);
catch me
  if strcmp(me.message, 'time axes of all trials should be identical')
    fprintf('error caught by function\n');
  end
end

cfg.t_ftimwin = 0.5;
cfg.toi       = 0.25;
ft_mvaranalysis(cfg,data);

% load in the data provided by Tyler Grummett
filename = dccnpath('/home/common/matlab/fieldtrip/data/test/bug2784.mat');
load(filename);

cfg = [];
mdata = ft_mvaranalysis(cfg, temp);

% the issue Tyler had could be reproduced and was caused by a combination
% of things: the tfwin being 501 samples long, as well as the time axes of
% the trials being different. ft_mvaranalysis has been adjusted to more
% robustly deal with these cases.


text:
function BiObj = Estimation(PopObj,r)
% Estimate the proximity and crowding degree of each solution

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    N = size(PopObj,1);

    %% Proximity estimation
    fmax   = repmat(max(PopObj,[],1),N,1);
    fmin   = repmat(min(PopObj,[],1),N,1);
    PopObj = (PopObj-fmin)./(fmax-fmin);
    fpr    = sum(PopObj,2);

    %% Crowding degree estimation
    d     = pdist2(PopObj,PopObj);
    d(logical(eye(length(d)))) = inf;
    fprm  = repmat(fpr,1,N);
    case1 = d<r & fprm<=fprm';
    case2 = d<r & fprm>fprm';
    sh        = zeros(N);
    sh(case1) = (0.5*(1-d(case1)/r)).^2;
    sh(case2) = (1.5*(1-d(case2)/r)).^2;
    fcd   = sqrt(sum(sh,2));
    BiObj = [fpr,fcd];
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function f_x = ParFor2(in1)
%PARFOR2
%    F_X = PARFOR2(IN1)

%    This function was generated by the Symbolic Math Toolbox version 8.1.
%    16-Jul-2019 11:57:35

u = in1(:,1);
ux = in1(:,4);
uxx = in1(:,5);
uxxx = in1(:,6);
f_x = (u.*2.6068e4+ux.*7.4257e4+uxx.*2.6372e5+uxxx.*4.3251e5-u.*ux.*1.0e5+u.*uxx.*2.3836e4-u.*uxxx.*1.0264e6)./(u.*2.6078e4+2.6069e5);

author: dynamicslab

repo: SINDy-PI


text:
% @author: B. Schauerte
% @date:   2013
% @url:    http://cvhci.anthropomatik.kit.edu/~bschauer/

% create sub-folder
libs_folder = 'libs';
mkdir(libs_folder);

% install the spectral saliency toolbox
mkdir(libs_folder,'svst');
unzip('http://www.mathworks.com/matlabcentral/fileexchange/32455-spectral-phase-based-visual-saliency?download=true',fullfile(libs_folder,'svst',''));
cd(fullfile('libs','svst',''));
get_additional_files();
cd(fileparts(mfilename('fullpath')));

% install the progressbar package
mkdir(libs_folder,'prog');
unzip('http://www.mathworks.com/matlabcentral/fileexchange/6922-progressbar?download=true',fullfile(libs_folder,'prog',''));


text:
function [distance_precision, PASCAL_precision, average_center_location_error] = ...
    compute_performance_measures(positions, ground_truth, distance_precision_threshold, PASCAL_threshold)

% [distance_precision, PASCAL_precision, average_center_location_error] = ...
%    compute_performance_measures(positions, ground_truth, distance_precision_threshold, PASCAL_threshold)
%
% For the given tracker output positions and ground truth it computes the:
% * Distance Precision at the specified threshold (20 pixels as default if
% omitted)
% * PASCAL Precision at the specified threshold (0.5 as default if omitted)
% * Average Center Location error (CLE).
%
% The tracker positions and ground truth must be Nx4-matrices where N is
% the number of time steps in the tracking. Each row has to be on the form
% [c1, c2, s1, s2] where (c1, c2) is the center coordinate and s1 and s2
% are the size in the first and second dimension respectively (the order of
% x and y does not matter here).

if nargin < 3 || isempty(distance_precision_threshold)
    distance_precision_threshold = 20;
end
if nargin < 4 || isempty(PASCAL_threshold)
    PASCAL_threshold = 0.5;
end

if size(positions,1) ~= size(ground_truth,1),
    disp('Could not calculate precisions, because the number of ground')
    disp('truth frames does not match the number of tracked frames.')
    return
end

%calculate distances to ground truth over all frames
distances = sqrt((positions(:,1) - ground_truth(:,1)).^2 + ...
    (positions(:,2) - ground_truth(:,2)).^2);
distances(isnan(distances)) = [];

%calculate distance precision
distance_precision = nnz(distances < distance_precision_threshold) / numel(distances);

%calculate average center location error (CLE)
average_center_location_error = mean(distances);

%calculate the overlap in each dimension
overlap_height = min(positions(:,1) + positions(:,3)/2, ground_truth(:,1) + ground_truth(:,3)/2) ...
    - max(positions(:,1) - positions(:,3)/2, ground_truth(:,1) - ground_truth(:,3)/2);
overlap_width = min(positions(:,2) + positions(:,4)/2, ground_truth(:,2) + ground_truth(:,4)/2) ...
    - max(positions(:,2) - positions(:,4)/2, ground_truth(:,2) - ground_truth(:,4)/2);

% if no overlap, set to zero
overlap_height(overlap_height < 0) = 0;
overlap_width(overlap_width < 0) = 0;

% remove NaN values (should not exist any)
valid_ind = ~isnan(overlap_height) & ~isnan(overlap_width);

% calculate area
overlap_area = overlap_height(valid_ind) .* overlap_width(valid_ind);
tracked_area = positions(valid_ind,3) .* positions(valid_ind,4);
ground_truth_area = ground_truth(valid_ind,3) .* ground_truth(valid_ind,4);

% calculate PASCAL overlaps
overlaps = overlap_area ./ (tracked_area + ground_truth_area - overlap_area);
%-lx calculate the precision with different threshold(1-50)
pre_o=zeros(1,100);
pre_d=zeros(1,50);
for thre=1:100
pre_o(thre)=nnz(overlaps >=(thre/100))/numel(overlaps);
end

for thre=1:50
pre_d(thre)=nnz(distances < thre) / numel(distances);
end
% figure,plot([1:100]/100,pre_o);
PASCAL_precision=sum(pre_o)/100;
%save  'result.mat' pre_o pre_d;
%-lx 2015-7-1
% calculate PASCAL precision
%PASCAL_precision = nnz(overlaps >= PASCAL_threshold) / numel(overlaps);

end

author: XinLi-zn

repo: TADT

save_path: github-repos/MATLAB/XinLi-zn-TADT


text:
function CostLayer = F_mixture_mse(input_layers, CostLayer)
predicted{1} = input_layers{1}.a;
predicted{2} = input_layers{2}.a;
ref{1} = input_layers{3}.a;
ref{2} = input_layers{4}.a;

[D,T,N] = size(ref{1});

DEBUG = 0;

ref_idx = [];
for i=1:N

    diff1 = [predicted{1}(:,:,i)-ref{1}(:,:,i) predicted{2}(:,:,i)-ref{2}(:,:,i)];
    cost1 = 0.5/T * sum(sum( diff1 .* conj(diff1) ));   % support both real and complex numbers

    diff2 = [predicted{1}(:,:,i)-ref{2}(:,:,i) predicted{2}(:,:,i)-ref{1}(:,:,i)];
    cost2 = 0.5/T * sum(sum( diff2 .* conj(diff2) ));   % support both real and complex numbers

    if cost1>cost2
        ref_idx(:,i) = [2 1];   % record the pair information by remembering the reference index
        cost(i) = cost2;
    else
        ref_idx(:,i) = [1 2];
        cost(i) = cost1;
    end

    if DEBUG
        subplot(1,2,1); imagesc( [predicted{1}(1:257,:,i); ref{ref_idx(1,i)}(1:257,:,i)]); colorbar
        subplot(1,2,2); imagesc( [predicted{2}(1:257,:,i); ref{ref_idx(2,i)}(1:257,:,i)]); colorbar
%         subplot(2,2,1); imagesc(predicted{1}(1:257,:,i)); colorbar
%         subplot(2,2,2); imagesc(predicted{2}(1:257,:,i)); colorbar
%         subplot(2,2,3); imagesc(ref{ref_idx(1,i)}(1:257,:,i)); colorbar;
%         subplot(2,2,4); imagesc(ref{ref_idx(2,i)}(1:257,:,i)); colorbar
        pause
    end
end

CostLayer.a = mean(cost);
CostLayer.ref_idx = ref_idx;

end

author: singaxiong

repo: SignalGraph


text:
function a = r8vec_permute_cyclic ( n, k, a )

%*****************************************************************************80
%
%% R8VEC_PERMUTE_CYCLIC performs a cyclic permutation of an R8VEC.
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%    For 0 <= K < N, this function cyclically permutes the input vector
%    to have the form
%
%     ( A(K+1), A(K+2), ..., A(N), A(1), ..., A(K) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of objects.
%
%    Input, integer K, the increment used.
%
%    Input/output, real A(N), the array to be permuted.
%
  b = zeros ( size ( a ) );

  for i = 1 : n
    ipk = i4_wrap ( i + k, 1, n );
    b(i) = a(ipk);
  end

  a(1:n) = b(1:n);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function h = cbarDraw(cbar, parent);
%
% h = cbarDraw(cbar, [parent=gca]);
%
% Render a color bar, according to the settings specified in the cbar
% struct.
%
% These color bars differ from the built-in matlab color bar tools in a few
% basic respects:
%
%   * They draw colorbars based on any colormap the user inputs, rather
%   than the color map assigned to the current figure. This is useful for
%   illustrating color codes on true color images, or images with overlaid
%   responses (like fMRI activation maps)
%
%   * In addition to being able to render vertical and horizontal color bars,
%   these tools may also render 'color wheel' colorbars, showing e.g. polar
%   angle. Other possible renderings (like rings) may be added down the
%   line.
%
% cbar is a struct with the following fields:
%
%   cmap: color map (nColors x 3) for the color bar. (Columns ar [R G B],
%   from 0-255).
%
%   nColors: # of colors to use in the cmap.
%
%   clim: color limits (aka 'clip mode'), which determines primarily
%   the labeling of the color bar. Can be set to 'auto', in which case
%   the labeling will be from 1:nColors. Otherwise, will label according to
%   the clim values (format is [min max]).
%
%   colorWheel: use a color wheel instead of a bar (e.g., to show polar
%   angle for a polar angle map).
%
%   colorWheelStart: degrees clockwise from 12-o-clock which map to beginning of
%   color map.
%
%   colorWheelDirection: direction of the color wheel. Can be:
%   1 or 'clockwise' (equivalent); or, 2 or 'counterclockwise' (equiv.)
%
%   colorWheelExtent: degrees (1 - 360) subtended by the color map, for polar
%   angle maps.
%
% ras, 08/2006
if ~exist('cbar', 'var') | isempty(cbar), cbar = cbarDefault; end
if ~exist('parent', 'var') | isempty(parent), parent = gca; end

if cbar.colorWheel==1
    % draw a color wheel
    if isequal(lower(cbar.colorWheelDirection), 'clockwise')
        direction = 0;
    else
        direction = 1;
    end
    startAngle = cbar.colorWheelStart;
	cmap = cbar.cmap;
	p.doPlot = 0; p.trueColor = 1; p.background = get(gcf, 'Color');
	p.visualField = 'b'; % 'both' fields in cmapWedge
	h = image(cmapWedge(cmap, startAngle, direction, p), 'Parent', parent);
    axis(parent, 'image'); axis(parent, 'off');

%     AX = axis;
%     xText = AX(2) + .1*(AX(2)-AX(1));
%     yText = AX(3) + .5*(AX(4)-AX(3));
%     text(xText, yText, cbar.label, 'Parent', parent);

%     % correct funky aspect ratios (gum):
%     pos = get(parent, 'Position');
%     aspectRatio = pos(3) / pos(4);
%     if aspectRatio > 4 | aspectRatio < 1/4
%         set(parent, 'Position', [pos(1) 0 pos(3) 1]);
%     end

elseif isequal(cbar.direction, 'horiz');
    % draw a horizontal bar
    if isequal(cbar.clim, 'auto') | isempty(cbar.clim)
        clim = [1 cbar.nColors];
    else
        clim = cbar.clim;
    end
    img = ind2rgb(1:cbar.nColors, cbar.cmap);

	% show the image
	% (The Y values maintain an aspect ratio of 35:1)
	aspRatio = 1/15;
    h = image(clim, clim ./ aspRatio, img, 'Parent', parent);

% 	axis(parent, 'image');
% 	axis(parent, [clim clim .* aspRatio]);
    set(parent, 'YTick', []);
    cbar.label(cbar.label=='_') = '-';
    xlabel(parent, cbar.label, 'FontSize', cbar.fontsz, ...
		  'FontName', cbar.font, 'Interpreter', 'none');

else
    % draw a vertical bar
    if isequal(cbar.clim, 'auto')
        clim = [1 cbar.nColors];
    else
        clim = cbar.clim;
	end
    img = ind2rgb([cbar.nColors:-1:1]', cbar.cmap);

	% show the image
	% (The X values maintain an aspect ratio of 15:1)
	aspRatio = 1/15;
	h = image(clim ./ 20, fliplr(clim), img, 'Parent', parent);

% 	axis(parent, 'image');
	axis(parent, [clim .* aspRatio, clim]);
    set(parent, 'XTick', [], 'Ydir', 'normal');
	if checkfields(cbar, 'labelSide') & cbar.labelSide==1
		set(parent, 'YAxisLocation', 'right');
	end
    cbar.label(cbar.label=='_') = '-';
    ylabel(parent, cbar.label, 'FontSize', cbar.fontsz, ...
			'FontName', cbar.font, 'Interpreter', 'none');

end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [filename] = segFilename(present,uid,iid)
% function [filename] = segFileName(present,uid,iid)
%
% Return seg filename composed of components.
%
% INPUT
%	present		Presentation, one of {'gray','color'}.
%	uid		User ID.
%	iid		Image ID.
%
% OUTPUT
%	filename	Segmentation filename.
%
% David Martin <dmartin@eecs.berkeley.edu>
% January 2003

filename = fullfile(bsdsRoot,'human',present,sprintf('%d',uid),sprintf('%d.seg',iid));

author: Cloud-CV

repo: object-proposals


text:
function f = p09_f ( n, x )

%*****************************************************************************80
%
%% P09_F evaluates the objective function for problem 9.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the argument of the objective function.
%
%    Output, real F, the value of the objective function.
%
  ap = 0.00001;

  t1 = -1.0;
  t2 = 0.0;
  t3 = 0.0;
  d2 = 1.0;
  s2 = 0.0;

  for j = 1 : n
    t1 = t1 + ( n - j + 1 ) * x(j)^2;
    s1 = exp ( x(j) / 10.0 );
    if ( 1 < j )
      s3 = s1 + s2 - d2 * ( exp ( 0.1 ) + 1.0 );
      t2 = t2 + s3 * s3;
      t3 = t3 + ( s1 - 1.0 / exp ( 0.1 ) )^2;
    end
    s2 = s1;
    d2 = d2 * exp ( 0.1 );
  end

  f = ap * ( t2 + t3 ) + t1 * t1 + ( x(1) - 0.2 )^2;

  return
end


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcFreeROI(Temp,Event,handles)
global Figure_handles;
MU_main_handles=guidata(Figure_handles.MU_main);
handles = guidata(handles.MU_matrix_display);

MU_enable('off',[],handles);
ROI_h=imfreehand;
MU_enable('on',[],handles);

MU_main_handles.V.ROIs{end+1,1}='imfreehand';
MU_main_handles.V.ROIs{end,2}=ROI_h;
MU_main_handles.V.ROIs{end,3}=getPosition(ROI_h);
ROI_ind=length(MU_main_handles.V.ROIs(:,1));

p=round(getPosition(ROI_h));
if max(p(:,1))>handles.V.Column | max(p(:,2))>handles.V.Row | min(p(:,1))<1 | min(p(:,2))<1
    delete(ROI_h);
    errordlg('Out of range subscript.');
    return;
end
BW=createMask(ROI_h);
p=[min(p(:,1)) min(p(:,2)) max(p(:,1))-min(p(:,1)) max(p(:,2))-min(p(:,2))];
TTMatrix=handles.BMatrix(p(2):p(2)+p(4),p(1):p(1)+p(3));
TTMatrix=TTMatrix(double(BW(p(2):p(2)+p(4),p(1):p(1)+p(3)))~=0);
ROI_Stat_h=text(p(1)+p(3),p(2)+p(4),{[' ROI#: ' num2str(ROI_ind)]; ...
                                     [' mean: ' num2str(mean(double(TTMatrix(:))))]; ...
                                     [' sd:' num2str(std(double(TTMatrix(:))))]; ...
                                     [' sd(%):' num2str(abs(std(double(TTMatrix(:)))./mean(double(TTMatrix(:)))*100))]},...
                                     'FontSize',10,'Color','g');

handles.V.ROI=struct(...
                     'ROI_flag', 1,...
                     'ROI_mov',[],...  % ROI movement track
                     'ROI_Stat_h', ROI_Stat_h,... % ROI stats
                     'ROI_h', ROI_h ... % ROI handle
                     );
handles.ROIData=TTMatrix;
guidata(handles.MU_matrix_display, handles);
guidata(Figure_handles.MU_main,MU_main_handles);

addNewPositionCallback(ROI_h,@(p) MU_ROI_stat(p,ROI_h,ROI_ind,handles));
fcn=makeConstrainToRectFcn('imfreehand',[0.5 handles.V.Column+0.4],[0.5 handles.V.Row+0.4]);
setPositionConstraintFcn(ROI_h,fcn);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [R_ned_2_inert, ned_x, ned_y, ned_z] = computeNedFrame(ut, rVectECI, bodyInfo)
    %Source: https://en.wikipedia.org/wiki/North_east_down

    [phi, lambda, ~, ~, ~, ~, ~, ~, REci2Ecef] = getLatLongAltFromInertialVect(ut, rVectECI, bodyInfo);

    REcef2Ned = [-sin(phi)*cos(lambda), -sin(lambda), -cos(phi)*cos(lambda);
                 -sin(phi)*sin(lambda),  cos(lambda), -cos(phi)*sin(lambda);
                  cos(phi),              0,           -sin(phi)]';

    R_ned_2_inert = (REcef2Ned * REci2Ecef)';
    ned_x = R_ned_2_inert(:,1);
    ned_y = R_ned_2_inert(:,2);
    ned_z = R_ned_2_inert(:,3);
end


text:
function test_tbgppca_subset

data = tbload_temperature;
data = tbpreprocess(data);
data = tbsubset(data, 1:1:79, 1:10:89000);

data

% $$$ figure
% $$$ plot(data.observations');

inW = data.coordinates(:,1:2)';
inX = data.time';

[M,N] = size(data.observations);
D = 2;

% Covariance matrices
logthetaW = {[log(1);log(5e1)], [log(1);log(1e1)]}
logthetaX = {log(2000), log(100)};

% Covariance functions
gpcovW = {@gpcovScale, @(logtheta,x1,x2) gpcov(logtheta,x1,x2,@sqdistEarth)};
% $$$ pseudoW = cell(D,1);%zeros([dimW,Mp,D]);
% $$$ Mp = ceil(1*M);
% $$$ for d=1:D
% $$$   permM = randperm(M);
% $$$   pseudoW{d} = inW(:,permM(1:Mp));% + randn(2,Mp);
% $$$ end

gpcovX = @gpcov;

% Generate latent variables
for d=1:D
  covfunW{d} = gpcovW;
  covfunX{d} = gpcovX;
end

Y = data.observations;

% Learn GP PCA
maxiter = 3e1;
Q = vbgppcamv(Y,D,inW,inX,covfunW,logthetaW,covfunX,logthetaX, ...
              'maxiter',maxiter, 'pseudodensityx', 0.01, 'pseudodensityw', ...
              0.4, 'loglikelihood', true, 'updatehyper', 5, 'maxsearchx', ...
              3, 'maxsearchw', 50);

% $$$ % Learn other models
% $$$ if ppca
% $$$   Qppca = pca_full(Y,D,'maxiters',maxiter,'rotate2pca',true, ...
% $$$                    'algorithm','ppca');
% $$$   Yh_ppca = Qppca.A * Qppca.S + repmat(Qppca.Mu,1,N);
% $$$   testrmse_ppca = rmse(Ytest, Yh_ppca)
% $$$   noiselessrmse_ppca = rmse(Y, Yh_ppca)
% $$$ end
% $$$ if vbpca
% $$$   Qvbpca = vbpcamv(Ynm,M-1,'maxiters',maxiter);
% $$$   Yh_vbpca = Qvbpca.W * Qvbpca.X + repmat(Qvbpca.mu,1,N);
% $$$   testrmse_vbpca = rmse(Ytest, Yh_vbpca)
% $$$   noiselessrmse_vbpca = rmse(Y, Yh_vbpca)
% $$$ end

Yh_gppca = Q.W * Q.X;
rmse_gppca = rmse(Y, Yh_gppca)
%testrmse_gppca = rmse(Ytest, Yh_gppca)
%noiselessrmse_gppca = rmse(Y, Yh_gppca)

figure
subplot(5,1,1);
plot(Y');
title('Original data')
% $$$ subplot(5,1,2);
% $$$ plot(Ynm');
% $$$ title('Observed data')
subplot(5,1,3);
plot(Yh_gppca');
title('Reconstruction of GP VB PCA')
% $$$ if vbpca
% $$$   subplot(5,1,4);
% $$$   plot(Yh_vbpca');
% $$$   title('Reconstruction of VB PCA')
% $$$ end
% $$$ if ppca
% $$$   subplot(5,1,5);
% $$$   plot(Yh_ppca');
% $$$   title('Reconstruction of PPCA')
% $$$ end

vX = Q.varX;%diag(Q.CovX);
eX = 2*sqrt(vX);%sqrt( vX(reshape(1:(N*D), D, N)) );
vW = Q.varW;%diag(Q.CovW);
eW = 2*sqrt(vW);%sqrt( vW(reshape(1:(M*D), M, D)) );
tsgpplot(inX, Q.X', eX', 'pseudoinputs', {Q.pseudoX});

%tsgpplot(inW, Q.W, eW, 'pseudoinputs', {Q.pseudoW});
figure
for d=1:D
  subplot(D,1,d);
  mapproj('testbed');
  mapplot(inW(1,:),inW(2,:),'ro');
  hold on
  mapplot(Q.pseudoW{d}(1,:),Q.pseudoW{d}(2,:),'k+');
  mapcoast
end

figure
mapproj('testbed');
gpmapcolor(Q.pseudoW, Q.Wp, Q.CovWp, 22:0.1:27, 59:0.1:61, Q.logthetaW, covfunW);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
% This demonstrates the smoothing properties of the Kalman lagged filter.
% Hidden states follow a triangular wave, whose observation is perturbed
% with white noise. Critically, we render the inversion scheme blind during
% half a period of the oscillation. We then invert the model (under AR
% priors on hidden states), with and without large backward lag.

clear variables
close all

nt = 10;
eta = 0.*randn(1,3*nt);
x = 1:nt;
x = [x,fliplr(x),x] + eta;

e = randn(1,3*nt);
y = x+e;

theta = [];
phi = [];
alpha = 1/var(eta);
sigma = 1/var(e);

displaySimulations(y,x,eta,e);

f_fname = @f_AR;
g_fname = @g_Id;
dim.n_theta = size(theta,1);
dim.n_phi = size(phi,1);
dim.n = size(x,1);

options.isYout = zeros(1,3*nt);
options.isYout(nt:2*nt+1) = 1;
options.MaxIterInit = 0;
options.priors.a_alpha = 1;
options.priors.b_alpha = 1;

% VB-Kalman-filter
[p1,o1] = VBA_NLStateSpaceModel(y,[],f_fname,g_fname,dim,options);

% VB-Kalman-smoother (lag = size of blind window)
options.backwardLag = nt+1;
[p2,o2] = VBA_NLStateSpaceModel(y,[],f_fname,g_fname,dim,options);
set(gcf,'name',['Kalman lag = ',num2str(o2.options.backwardLag)])

VBA_ReDisplay(p1,o1,1);
set(gcf,'name',['Kalman lag = ',num2str(o1.options.backwardLag)])

author: MBB-team

repo: VBA-toolbox


text:
function [Dx, Vx, mean_suv, max_suv, min_suv, Slope]=analyze_ivh(structNum, scanSet, plot_flag)
% IVH analysis
global planC
indexS=planC{end};
optS.IVHBinWidth=0.05;
[scansV, volsV] = getIVH(structNum, scanSet, planC);
[scanBinsV, volsHistV] = doseHist(scansV, volsV, optS.IVHBinWidth);
cumVolsV = cumsum(volsHistV);
cumVols2V  = cumVolsV(end) - cumVolsV;  %cumVolsV is the cumulative volume lt that corresponding scan
% calculate stats
mean_suv= calc_meanDose(scanBinsV, volsHistV);
max_suv=  calc_maxDose(scanBinsV, volsHistV);
min_suv=  calc_minDose(scanBinsV, volsHistV);
D50 = calc_Dx(scanBinsV, volsHistV,50)

Slope=calc_Slope(scanBinsV, volsHistV, D50, 0);

range = max_suv-min_suv;
param = [10:10:90];
for i=1:length(param)
    per=param(i)*range/100+min_suv; % ith percentile
    Dx(i)= calc_Dx(scanBinsV, volsHistV,param(i));
    Vx(i)= calc_Vx(scanBinsV, volsHistV,per,1);
end

%including that scan bin.
% if plot_flag
% %     if ~isempty(planC{indexS.scan}(scanSet).scanOffset)
% %         h = plot([0, scanBinsV - planC{indexS.scan}(scanSet).scanOffset], [1, cumVols2V/cumVolsV(end)]);
% %     else
% %        h = plot([0, scanBinsV], [1, cumVols2V/cumVolsV(end)]);
% %     end
% end

return

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
 function scRay(point, LinCol)
%scRay : Draws a ray emanating from the origin towards a given point on smith chart
%
%  SYNOPSIS:
%     Draws a ray emanating from the origin towards a given point on smith chart.
%     One may need to do so when solving a transmission line problem using a smith
%     chart. This may be useful, because it also shows the phase angle and the
%     transmission line length on its outer extrimity which need not be read manualy
%     from the outer scales of the smith chart.
%
%
%  SYNTAX:
%     scRay(point, LinCol)
%
%  INPUT ARGUMENTS:
%     point  : coordinates of the given point [r x] normalized
%     LinCol : color of the ray
%
%  OUTPUT ARGUMENT:
%     none
%
%  EXAMPLE:
%         The Command sequence
%         scDraw;
%         scRay([2 3]);
%         will draw a blank smith chart and draw a line from smith chart
%         center [1 0] (corresponding to impedence Z_L) to the point [2 3]
%         in the smith chart corresponding to an impedance (2+j3)*Z_L. This
%         may be thought of locating an load impedance of (2+j3)*Z_L on the
%         smith chart
%
%
%     Mohammad Ashfaq - (31-05-2000)
%     Mohammad Ashfaq - (13-04-2006) Modified (example included)
%
 if nargin == 1
     LinCol='m';
 end
 r1 = point(1);
 x1 = point(2);
 [u1, v1] = scPOI(r1,x1);

 % MARK POINT
 plot(u1,v1,'r*')

 Theta    = atan2(v1,u1);

 plot([0 u1],[0 v1],LinCol);

 plot([1.35*cos(Theta) u1],[1.35*sin(Theta) v1],'k');

 h = text(1.4*cos(Theta),1.4*sin(Theta), '\theta=');
 set(h, 'Fontsize', 9, 'HorizontalAlignment','center','Color','r');
 h = text(1.4*cos(Theta)+0.08,1.4*sin(Theta), ['     ', num2str(Theta*180/pi,'%3.2f')]);
 set(h, 'Fontsize', 8, 'HorizontalAlignment','center','Color','k');

 h = text(1.4*cos(Theta),1.4*sin(Theta)-0.05, 'l / \lambda=');
 set(h, 'Fontsize', 9, 'HorizontalAlignment','center','Color','r');

 h = text(1.4*cos(Theta)+0.08,1.4*sin(Theta)-0.05, ['         ',num2str((1-Theta/pi)/4, '%0.3f')]);
 set(h, 'Fontsize', 8, 'HorizontalAlignment','center','Color','k');

if 0

 string1 = ['\theta=', num2str(Theta*180/pi,'%3.2f')];
 string2 = ['l/\lambda=', num2str((1-Theta/pi)/4, '%0.3f')];

 Thetad = Theta*180/pi;

 h1 = text(1.4*cos(Theta),1.4*sin(Theta)    , string1);
 h2 = text(1.4*cos(Theta),1.4*sin(Theta)-0.2, string2);
 set(h1, 'Fontsize', 8, 'HorizontalAlignment','right');
 set(h2, 'Fontsize', 8, 'HorizontalAlignment','right');

 if abs(Thetad)>90
    Thetad  = Thetad+180;
    set(h1,'rotation', Thetad, 'Fontsize', 8, 'HorizontalAlignment','right');
    set(h2,'rotation', Thetad, 'Fontsize', 8, 'HorizontalAlignment','right');
 else
    set(h1,'rotation', Thetad, 'Fontsize', 8);
    set(h2,'rotation', Thetad, 'Fontsize', 8);
 end
end


text:
%% t_mrdMBA
%
%  Show fibers using FP's mba code
%  Make sure AFQ is on your path
%  SPM needs to be there, too.
%  The repository can be downloaded from:
%     git clone
%
%%

remote.host = 'http://scarlet.stanford.edu/validation/MRI/VISTADATA';
remote.directory = fullfile('diffusion','sampleData','fibers');
rdata('cd',remote);
dirList = rdata('ls',remote,'pdb');                   % Match .pdb extension
%
oname = rdata('get',remote,dirList{2});
fname = [oname,'.pdb'];
copyfile(oname,fname);

%%  Load the fibers and display
fg = fgRead(fname);

% Render the Tract FA Profile for the left uncinate
% A small number of triangles (25 is the default).
nTriangles = 8;
[lgt , fSurf, fvc] = AFQ_RenderFibers(fg,'subdivs',nTriangles);

%%
[lgt , fSurf, fvc] = AFQ_RenderFibers(fg,...
    'subdivs',nTriangles,...
    'color',[.8 0 .8]);

%%
[lgt , fSurf, fvc] = AFQ_RenderFibers(fg,...
    'subdivs',nTriangles,...
    'color',[.6 .6 .8], ...
    'jittercolors',.4);
%%
[~, AFQdata] = AFQ_directories;
cortex = fullfile(AFQdata,'mesh','segmentation.nii.gz');
overlay = fullfile(AFQdata,'mesh','Left_Arcuate_Endpoints.nii.gz');
thresh = .01; % Threshold for the overlay image
crange = [.01 .8]; % Color range of the overlay image
% Render the cortical surface colored by the arcuate endpoint density
[p, msh, lightH] = AFQ_RenderCorticalSurface(cortex, 'overlay' , overlay, 'crange', crange, 'thresh', thresh)

% End


text:
function edgewt = SC_EdgeResponse(thickness)
% Ideal model: k =  EdgeCoef0, (0.4), EdgeRange = EdgeRangeMM, (80)
% i.e., y = k*x + 0.78
% Herein, we used an approximation method to model the edge response
% function
% Date: 2021-05-24

% Emperical Thickness Threshold for Edge Segementation
edgewt = double(imbinarize(thickness, 50));
tmpmask = edgewt;
h = fspecial('average', [25, 25]);
for ii = 1:5
    edgewt = imfilter(edgewt, h);
end
tmp = tmpmask.*edgewt;
tmp(tmp==0) = NaN;
tmp = rescale(tmp, 0.6, 1);
tmp(isnan(tmp)) = 0;
edgewt = tmp;

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function v = mono_value ( d, nx, f, x )

%*****************************************************************************80
%
%% MONO_VALUE evaluates a monomial.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 December 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer D, the spatial dimension.
%
%    Input, integer NX, the number of evaluation points.
%
%    Input, integer F(D), the exponents of the monomial.
%
%    Input, real X(D,NX), the coordinates of the evaluation points.
%
%    Output, real V(NX), the value of the monomial at X.
%
  v = ones ( nx, 1 );

  for i = 1 : d
    v(1:nx,1) = v(1:nx,1) .* x(i,1:nx)' .^ f(i);
  end

  return
end


text:
function [ rtn ] = TESTZERODAY( denum, testCase )
%TEST test the known zero-day values from JPL's ephemerides
  output = 1;
  rtn = 0;

  if nargin < 1
    denum = 405;
  end
  if isnumeric(denum)
    denum = sprintf('%d',denum);
  end
  if nargin < 2
    testCase = -1;
  end
  [ st, error ] = Ephem.readHeader( denum );
  if error == 0 && isfield(st,'JDEPOC') % contains zero data
    year = floor(st.JDEPOC/365.2425 - 4712.019); % possibly off by +1 year at New Years Day
    s = Ephem(denum);
    s.output = output;
    s.KM = 0;
    pos_err = 5.0e-8;
    vel_err = 7.0e-6;
    for i=1:2
      if i == 1
        fprintf(output,'Testing ascii ephemeris DE%s zero day\n',denum);
        fn = Ephem.asc_name(denum,year);
        [s,error] = s.openDEasc(denum, fn);
      elseif i == 2
        fprintf(output,'Testing binary ephemeris DE%s zero day\n',denum);
        fn = Ephem.eph_name(denum,year);
        [s,error] = s.openDEeph(denum, fn);
      end
      if error
        %if nargin > 1
        %  testCase.verifyEqual(error,0,'bad openDE');
        %end
        break;
      end
      target = Ephem.Mercury;
      [ target_pos, target_vel,~, error ] = s.state( st.JDEPOC, target );
      if ~error && isfield(st,'X1')
        check(st,'1',target_pos,target_vel,pos_err,vel_err,testCase);
      end
      target = Ephem.Venus;
      [ target_pos, target_vel,~, error ] = s.state( st.JDEPOC, target );
      if ~error && isfield(st,'X2')
        check(st,'2',target_pos,target_vel,pos_err,vel_err,testCase);
      end
      target = Ephem.Earth;
      [ target_pos, target_vel,~, error ] = s.state( st.JDEPOC, target );
      if ~error && isfield(st,'XB')
        check(st,'B',target_pos,target_vel,pos_err,vel_err,testCase);
      end
      target = Ephem.Mars;
      [ target_pos, target_vel,~, error ] = s.state( st.JDEPOC, target );
      if ~error && isfield(st,'X4')
        check(st,'4',target_pos,target_vel,pos_err,vel_err,testCase);
      end
      target = Ephem.Jupiter;
      [ target_pos, target_vel,~, error ] = s.state( st.JDEPOC, target );
      if ~error && isfield(st,'X5')
        check(st,'5',target_pos,target_vel,pos_err,vel_err,testCase);
      end
      target = Ephem.Saturn;
      [ target_pos, target_vel,~, error ] = s.state( st.JDEPOC, target );
      if ~error && isfield(st,'X6')
        check(st,'6',target_pos,target_vel,pos_err,vel_err,testCase);
      end
      target = Ephem.Uranus;
      [ target_pos, target_vel,~, error ] = s.state( st.JDEPOC, target );
      if ~error && isfield(st,'X7')
        check(st,'7',target_pos,target_vel,pos_err,vel_err,testCase);
      end
      target = Ephem.Neptune;
      [ target_pos, target_vel,~, error ] = s.state( st.JDEPOC, target );
      if ~error && isfield(st,'X8')
        check(st,'8',target_pos,target_vel,pos_err,vel_err,testCase);
      end
      target = Ephem.Pluto;
      [ target_pos, target_vel,~, error ] = s.state( st.JDEPOC, target );
      if ~error && isfield(st,'X9')
        check(st,'9',target_pos,target_vel,pos_err,vel_err,testCase);
      end
      target = Ephem.Moon;
      [ target_pos, target_vel,~, error ] = s.state( st.JDEPOC, target );
      if ~error && isfield(st,'XM')
        check(st,'M',target_pos,target_vel,pos_err,vel_err,testCase);
      end
      if isfield(st,'XS') && st.XS ~= 0 && st.YS ~= 0 && st.ZS ~= 0 % has Sun data
        target = Ephem.Sun;
        [ target_pos, target_vel,~, error ] = s.state( st.JDEPOC, target );
        check(st,'S',target_pos,target_vel,pos_err,vel_err,testCase);
      end
    end
  end
end
function rtn = check(st,n,pos,vel,pos_err,vel_err,testCase)
  rtn = true;
  output = 1;
  real_pos = [st.(sprintf('X%s',n)),st.(sprintf('Y%s',n)),st.(sprintf('Z%s',n))];
  real_vel = [st.(sprintf('XD%s',n)),st.(sprintf('YD%s',n)),st.(sprintf('ZD%s',n))];
  diff = (pos(1) - real_pos(1))/real_pos(1);
  if abs(diff) >= pos_err
    str = sprintf('X%c %g ~= %g (%g)\n',n,pos(1),real_pos(1),diff);
    if nargin > 6 && testCase ~= -1
      testCase.verifyEqual(pos(1),real_pos(1),'RelTol',pos_err,str);
    else
      fprintf(output,str);
    end
    rtn = false;
  end
  diff = (pos(2) - real_pos(2))/real_pos(2);
  if abs(diff) >= pos_err
    str = sprintf('Y%c %g ~= %g (%g)\n',n,pos(2),real_pos(2),diff);
    if nargin > 6 && testCase ~= -1
      testCase.verifyLessThan(abs(diff),pos_err,str);
    else
      fprintf(output,str);
    end
    rtn = false;
  end
  diff = (pos(3) - real_pos(3))/real_pos(3);
  if abs(diff) >= pos_err
    str = sprintf('Z%c %g ~= %g (%g)\n',n,pos(3),real_pos(3),diff);
    if nargin > 6 && testCase ~= -1
      testCase.verifyLessThan(abs(diff),pos_err,str);
    else
      fprintf(output,str);
    end
    rtn = false;
  end
  diff = (vel(1) - real_vel(1))/real_vel(1);
  if abs(diff) >= vel_err
    str = sprintf('XD%c %g ~= %g (%g)\n',n,vel(1),real_vel(1),diff);
    if nargin > 6 && testCase ~= -1
      testCase.verifyLessThan(abs(diff),vel_err,str);
    else
      fprintf(output,str);
    end
    rtn = false;
  end
  diff = (vel(2) - real_vel(2))/real_vel(2);
  if abs(diff) >= vel_err
    str = sprintf('YD%c %g ~= %g (%g)\n',n,vel(2),real_vel(2),diff);
    if nargin > 6 && testCase ~= -1
      testCase.verifyLessThan(abs(diff),vel_err,str);
    else
      fprintf(output,str);
    end
    rtn = false;
  end
  diff = (vel(3) - real_vel(3))/real_vel(3);
  if abs(diff) >= vel_err
    str = sprintf('ZD%c %g ~= %g (%g)\n',n,vel(3),real_vel(3),diff);
    if nargin > 6 && testCase ~= -1
      testCase.verifyLessThan(abs(diff),vel_err,str);
    else
      fprintf(output,str);
    end
    rtn = false;
  end
end

author: goGPS-Project