text:
function [Q,R,S,U,P] = spm_MDP(MDP)
% solves the active inference problem for Markov decision processes
% FROMAT [Q,R,S,U,P] = spm_MDP(MDP)
%
% MDP.T           - process depth (the horizon)
% MDP.S(N,1)      - initial state
% MDP.B{M}(N,N)   - transition probabilities among hidden states (priors)
% MDP.C(N,1)      - terminal cost probabilities (prior N over hidden states)
% MDP.D(M,1)      - control probabilities (prior over M control states)
%
% optional:
%
% MDP.W           - log-precision of beliefs about transitions (default: 1)
% MDP.G{M}(N,N)   - transition probabilities used to generate outcomes
%                   (default: the prior transition probabilities)
% MDP.A(N,N)      - Likelihood of outcomes given hidden states
%                   (default: an identity mapping from states to outcomes)
% MDP.B{T,M}(N,N) - transition probabilities for each time point
% MDP.G{T,M}(N,N) - transition probabilities for each time point
%                   (default: MDP.B{T,M} = MDP.B{M})
%
% MDP.plot        -  swtich to suppress graphics
%
% produces:
%
% Q(N,K,T) - an array of conditional (posterior) expectations over N hidden
%            states and time 1,...,T at time 1,...,K
% R(M,K,T) - an array of conditional expectations over M control
%            states and time 1,...,T at time 1,...,K
% S(N,T)   - a sparse matrix of ones, encoding the state at time 1,...,T
% U(M,T)   - a sparse matrix of ones, encoding the action at time 1,...,T
% P(M,T)   - probabaility of emitting action 1,...,M at time 1,...,T
%
% This routine provides solutions of active inference (minimisation of
% variational free energy)using a generative model based upon a Markov
% decision process. This model and inference scheme is formulated
% in discrete space and time. This means that the generative model (and
% process) are  finite state  machines or hidden Markov models whose
% dynamics are given by transition probabilities among states. For
% simplicity, we assume an isomorphism between hidden states and outcomes,
% where the likelihood corresponds to a particular outcome conditioned upon
% hidden states. Similarly, for simplicity, this routine assumes that action
% and hidden controls are isomorphic. If the dynamics of transition
% probabilities of the true process are not provided, this routine will use
% the equivalent probabilities from the generative model.
%
% The transition probabilities are a cell array of probability transition
% matrices corresponding to each (discrete) the level of the control state.
%
% Mote that the conditional expectations are functions of time but also
% contain expectations about fictive states over time at each time point.
% To create time dependent transition probabilities, one can specify a
% function in place of the transition probabilities under different levels
% of control.
%
% partially observed Markov decision processes can be modelled by
% specifying a likelihood (as part of a generative model) and absorbing any
% probabilistic mapping between (isomorphic) hidden states and outcomes
% into the transition probabilities G.
%
% See also spm_MDP_game
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_MDP.m 6061 2014-06-21 09:02:42Z karl $

% set up and preliminaries
%==========================================================================

% plotting and precision defaults
%--------------------------------------------------------------------------
try PLOT   = MDP.plot;   catch, PLOT   = 1; end
try lambda = MDP.lambda; catch, lambda = 1; end
try W      = MDP.W;      catch, W      = 1; end

% get figure
%--------------------------------------------------------------------------
if PLOT, spm_figure('GetWin','MDP'); clf,   end

% generative model and initial states
%--------------------------------------------------------------------------
P0    = exp(-32);             % smallest probability
T     = MDP.T;                % process depth (the horizon)
Ns    = size(MDP.B{1},1);     % number of hidden states
Nb    = size(MDP.B,1);        % number of time-dependent probabilities
Nu    = size(MDP.B,2);        % number of hidden controls

% likelihood model (for a partially observed MDP implicit in G)
%--------------------------------------------------------------------------
try
    A = MDP.A + P0;
catch
    A = speye(Ns,Ns) + P0;
end
A     = A*diag(1./sum(A));

% transition probabilities (priors)
%--------------------------------------------------------------------------
for i = 1:T
    for j = 1:Nu
        if i == 1 || Nb == T
            B{i,j}   = MDP.B{i,j} + P0;
            H{i,j}   = B{i,j}';
            B{i,j}   = B{i,j}*diag(1./sum(B{i,j}));
            H{i,j}   = H{i,j}*diag(1./sum(H{i,j}));
            lnB{i,j} = log(B{i,j})*W;
            lnH{i,j} = log(H{i,j})*W;
        else
            B{i,j}   = B{1,j};
            lnB{i,j} = lnB{1,j};
            lnH{i,j} = lnH{1,j};
        end
    end
end

% terminal cost probabilities (priors)
%--------------------------------------------------------------------------
try
    C = spm_vec(MDP.C) + P0;
catch
    C = ones(Ns,1);
end

% control probabilities (priors)
%--------------------------------------------------------------------------
try
    D = spm_vec(MDP.D) + P0;
catch
    D = ones(Nu,1);
end

% state probabilities (priors)
%--------------------------------------------------------------------------
try
    E = spm_vec(MDP.E) + P0;
catch
    E = ones(Ns,1);
end

C     = C/sum(C);
D     = D/sum(D);
E     = E/sum(E);
lnA   = log(A);
lnD   = log(D)*W;
lnE   = log(E)*W;

% generative process (assume the true process is the same as the model)
%--------------------------------------------------------------------------
try
    G = MDP.G;
catch
    G = MDP.B;
end
Ng    = size(G,1);
for i = 1:T
    for j = 1:Nu
        if i == 1 || Ng == T
            G{i,j}   = G{i,j}*diag(1./sum(G{i,j}));
        else
            G{i,j}   = G{1,j};
        end
    end
end

% effector action (e) and sufficient statistics of proposal density
%--------------------------------------------------------------------------
a        = zeros(Ns,T);                          % state probability
b        = ones(Nu,T);                           % control probability
a(:,T)   = C;                                    % final state
b        = b*diag(1./sum(b));

% posterior expectations (states Q, control R)
%--------------------------------------------------------------------------
Q(:,1,:) = a;
R(:,1,:) = b;

%  initial state and posterior (states Q, control R) and action (E)
%--------------------------------------------------------------------------
s     = find(spm_vec(MDP.S));
S     = sparse(s,1,1,Ns,T);
U     = sparse(Nu,T);

% solve
%==========================================================================
for k = 1:(T - 1)

    % forward and backward passes at this time point
    %----------------------------------------------------------------------
    for i = 1:8
        for t = (T - 1):-1:k

            % get data likelihood if available at this time
            %--------------------------------------------------------------
            if t > k
                at = lnE;
            else
                at = lnA'*S(:,t) + lnE;
            end

            % and accumulate empirical priors
            %--------------------------------------------------------------
            for j = 1:Nu
                if t > 1
                    at  = at     + b(j,t - 1) *lnH{t,j}'*a(:,t - 1);
                end
                at      = at     + b(j,t    ) *lnB{t,j}'*a(:,t + 1);
                bt(j,1) = lnD(j) + a(:,t    )'*lnB{t,j}'*a(:,t + 1);
            end

            % update sufficient statistics of hidden states and control
            %--------------------------------------------------------------
            a(:,t) = spm_softmax(at);
            b(:,t) = spm_softmax(bt);

            % graphics to inspect update scheme
            %==============================================================
            if PLOT > 2

                % posterior beliefs about hidden states
                %----------------------------------------------------------
                spm_plot_states(a,b)

                % pause if requested
                %----------------------------------------------------------
                if PLOT > 3, pause, end

            end

        end

        % graphics to inspect update scheme
        %==================================================================
        if PLOT > 1

            % posterior beliefs about hidden states
            %--------------------------------------------------------------
            spm_plot_states(a,b)

        end
    end

    % sampling of next state (outcome)
    %======================================================================
    for i = 1:Nu
        F(i,1) = B{k,i}(:,s)'*lnA*a(:,k + 1);
    end

    % next action (the action and minimises expected free energy)
    %----------------------------------------------------------------------
    Pu       = spm_softmax(F,lambda);
    i        = find(rand < cumsum(Pu),1);

    % next state (assuming G mediates uncertainty modelled the likelihood)
    %----------------------------------------------------------------------
    Ps       = spm_softmax(G{k,i}(:,s),lambda);
    s        = find(rand < cumsum(Ps),1);

    % save action, state and posterior expectations (states Q, control R)
    %----------------------------------------------------------------------
    P(:,k)     = Pu;
    Q(:,k,:)   = a;
    R(:,k,:)   = b;
    S(s,k + 1) = 1;
    U(i,k)     = 1;

    % plot
    %======================================================================
    if PLOT > 0

        % posterior beliefs about hidden states
        %------------------------------------------------------------------
        spm_plot_states(a,b)

        % states sampled (outcome)
        %------------------------------------------------------------------
        subplot(2,2,3)
        if size(S,1) > 512
            spy(S,16)
        else
            imagesc(1 - S)
        end
        axis square
        title('Sampled state','FontSize',16)
        xlabel('time','FontSize',12)
        ylabel('action','FontSize',12)

        % action sampled (selected)
        %------------------------------------------------------------------
        subplot(2,2,4)
        if size(U,1) > 512
            spy(U,16)
        else
            imagesc(1 - U)
        end
        axis square
        title('Selected action','FontSize',16)
        xlabel('time','FontSize',12)
        ylabel('action','FontSize',12)
        drawnow

    end

end

function spm_plot_states(a,b)

% posterior beliefs about hidden states
%--------------------------------------------------------------------------
subplot(2,2,1)
if size(a,1) > 512
    spy(a > 1/32,16)
else
    imagesc(1 - a)
end
axis square
title('Expected hidden states','FontSize',16)
xlabel('time','FontSize',12)
ylabel('hidden state','FontSize',12)

% posterior beliefs about control states
%--------------------------------------------------------------------------
subplot(2,2,2)
if size(b,1) > 512
    spy(b > 1/8,16)
else
    imagesc(1 - b)
end
axis square
title('Expected control states','FontSize',16)
xlabel('time','FontSize',12)
ylabel('control state','FontSize',12)
drawnow

author: spm


text:
function s = maxitmess()
%MAXITMESS Create a standard error message when training reaches max. iterations.
%
%	Description
%	S = MAXITMESS returns a standard string that it used by training
%	algorithms when the maximum number of iterations (as specified in
%	OPTIONS(14) is reached.
%
%	See also
%	CONJGRAD, GLMTRAIN, GMMEM, GRADDESC, GTMEM, KMEANS, OLGD, QUASINEW, SCG
%

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

s = 'Maximum number of iterations has been exceeded';

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function pos = xyzhumanevaJoint2pos(joint)

% XYZHUMANEVAJOINT2POS
%
% COPYRIGHT : Carl Henrik Ek and Neil Lawrence, 2008

% MOCAP

pos = zeros(1,prod(size(joint)));

pos(1:3:end) = joint(:,1);
pos(2:3:end) = joint(:,2);
pos(3:3:end) = joint(:,3);

return

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef Blender < handle
    %BLENDER  Class for all image blenders
    %
    % See also: cv.Stitcher
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    methods
        function this = Blender(blenderType, varargin)
            %BLENDER  Constructor
            %
            %     obj = cv.Blender(blenderType)
            %     obj = cv.Blender(blenderType, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __blenderType__ image blender type. One of:
            %   * __NoBlender__ Simple blender which puts one image over
            %     another.
            %   * __FeatherBlender__ Simple blender which mixes images at its
            %     borders.
            %   * __MultiBandBlender__ Blender which uses multi-band blending
            %     algorithm (see [BA83]).
            %
            % ## Options
            % The following are options for the various algorithms:
            %
            % ### `FeatherBlender`
            % * __Sharpness__ default 0.02
            %
            % ### `MultiBandBlender`
            % * __TryGPU__ default false
            % * __NumBands__ default 5
            % * __WeightType__ One of:
            %   * __single__ (default)
            %   * __int16__
            %
            % ## References
            % [BA83]:
            % > Peter J Burt and Edward H Adelson.
            % > "A multiresolution spline with application to image mosaics".
            % > ACM Transactions on Graphics (TOG), 2(4):217-236, 1983.
            %
            % See also: cv.Blender.blend
            %
            this.id = Blender_(0, 'new', blenderType, varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.Blender
            %
            if isempty(this.id), return; end
            Blender_(this.id, 'delete');
        end

        function typename = typeid(this)
            %TYPEID  Name of the C++ type (RTTI)
            %
            %     typename = obj.typeid()
            %
            % ## Output
            % * __typename__ Name of C++ type
            %
            typename = Blender_(this.id, 'typeid');
        end
    end

    %% Blender
    methods
        function prepare(this, varargin)
            %PREPARE  Prepares the blender for blending
            %
            %     obj.prepare(corners, sizes)
            %     obj.prepare(dst_roi)
            %
            % ## Input
            % * __corners__ Source images top-left corners. Cell-array of 2D
            %   points `{[x,y],...}`.
            % * __sizes__ Source image sizes. Cell-array of sizes
            %   `{[w,h],...}`.
            % * **dst_roi** Destination ROI rectangle `[x,y,w,h]`, same as
            %   obtained by calling `cv.Blender.resultRoi(corners,sizes)`.
            %
            % See also: cv.Blender.feed
            %
            Blender_(this.id, 'prepare', varargin{:});
        end

        function feed(this, img, mask, tl)
            %FEED  Processes the image
            %
            %     obj.feed(img, mask, tl)
            %
            % ## Input
            % * __img__ Source image, 3-channel of type `int16`.
            % * __mask__ Source image mask, 1-channel of type `uint8`.
            % * __tl__ Source image top-left corners `[x,y]`.
            %
            % See also: cv.Blender.blend
            %
            Blender_(this.id, 'feed', img, mask, tl);
        end

        function [dst, dst_mask] = blend(this, varargin)
            %BLEND  Blends and returns the final pano
            %
            %     [dst, dst_mask] = obj.blend()
            %     [dst, dst_mask] = obj.blend('OptionName',optionValue, ...)
            %
            % ## Output
            % * __dst__ Final pano, of type `int16`.
            % * **dst_mask** Final pano mask, of type `uint8`.
            %
            % ## Options
            % * __Dst__
            % * __Mask__
            %
            % See also: cv.Blender.feed
            %
            [dst, dst_mask] = Blender_(this.id, 'blend', varargin{:});
        end
    end

    %% FeatherBlender
    methods
        function [weight_maps, dst_roi] = createWeightMaps(this, masks, corners)
            %CREATEWEIGHTMAPS  Creates weight maps for fixed set of source images by their masks and top-left corners
            %
            %     [weight_maps, dst_roi] = obj.createWeightMaps(masks, corners)
            %
            % ## Input
            % * __masks__ Cell-array of image masks.
            % * __corners__ Cell-array of top-left corners `{[x,y],...}`.
            %
            % ## Output
            % * **weight_maps** Cell-array of weight maps.
            % * **dst_roi** ROI rectangle `[x,y,w,h]`.
            %
            % Final image can be obtained by simple weighting of the source
            % images.
            %
            [weight_maps, dst_roi] = Blender_(this.id, 'createWeightMaps', masks, corners);
        end
    end

    %% Auxiliary methods
    methods (Static)
        function pyr = createLaplacePyr(img, num_levels, varargin)
            %CREATELAPLACEPYR  Create Laplace pyramid
            %
            %     pyr = cv.Blender.createLaplacePyr(img, num_levels)
            %     pyr = cv.Blender.createLaplacePyr(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __img__ Input image.
            % * **num_levels** Number of pyramid levels.
            %
            % ## Output
            % * __pyr__ Laplace pyramid, cell-array of matrices.
            %
            % ## Options
            % * __UseGPU__ (requires CUDA). default false
            %
            % See also: cv.Blender.restoreImageFromLaplacePyr, cv.buildPyramid
            %
            pyr = Blender_(0, 'createLaplacePyr', img, num_levels, varargin{:});
        end

        function img = restoreImageFromLaplacePyr(pyr, varargin)
            %RESTOREIMAGEFROMLAPLACEPYR  Restore source image from Laplace pyramid
            %
            %     img = cv.Blender.restoreImageFromLaplacePyr(pyr)
            %     img = cv.Blender.restoreImageFromLaplacePyr(pyr, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __pyr__ Laplace pyramid, cell-array of matrices.
            %
            % ## Output
            % * __img__ Output image.
            %
            % ## Options
            % * __UseGPU__ (requires CUDA). default false
            %
            % See also: cv.Blender.createLaplacePyr, cv.buildPyramid
            %
            img = Blender_(0, 'restoreImageFromLaplacePyr', pyr, varargin{:});
        end

        function [roi,success] = overlapRoi(tl1, tl2, sz1, sz2)
            %OVERLAPROI  Overlap ROI
            %
            %     [roi,success] = cv.Blender.overlapRoi(tl1, tl2, sz1, sz2)
            %
            % ## Input
            % * __tl1__ First top-left corner `[x1,y1]`.
            % * __tl2__ Second top-left corner `[x2,y2]`.
            % * __sz1__ First size `[w1,h1]`.
            % * __sz2__ Second size `[w2,h2]`.
            %
            % ## Output
            % * __roi__ Output rectangle `[x,y,w,h]`.
            % * __success__ Success flag.
            %
            % See also: cv.Blender.resultRoi
            %
            [roi,success] = Blender_(0, 'overlapRoi', tl1, tl2, sz1, sz2);
        end

        function roi = resultRoi(corners, sizes)
            %RESULTROI  Result ROI
            %
            %     roi = cv.Blender.resultRoi(corners, sizes)
            %
            % ## Input
            % * __corners__ Cell-array of top-left corners `{[x,y],...}`.
            % * __sizes__ Cell-array of corresponding sizes `{[w,h],...}`.
            %
            % ## Output
            % * __roi__ Output rectangle `[x,y,w,h]`.
            %
            % See also: cv.Blender.resultRoiIntersection
            %
            roi = Blender_(0, 'resultRoi', corners, sizes);
        end

        function roi = resultRoiIntersection(corners, sizes)
            %RESULTROIINTERSECTION  Result ROI intersection
            %
            %     roi = cv.Blender.resultRoiIntersection(corners, sizes)
            %
            % ## Input
            % * __corners__ Cell-array of top-left corners `{[x,y],...}`.
            % * __sizes__ Cell-array of corresponding sizes `{[w,h],...}`.
            %
            % ## Output
            % * __roi__ Output rectangle `[x,y,w,h]`.
            %
            % See also: cv.Blender.resultRoi
            %
            roi = Blender_(0, 'resultRoiIntersection', corners, sizes);
        end

        function tl = resultTl(corners)
            %RESULTTL  Result top-left corner
            %
            %     tl = cv.Blender.resultTl(corners)
            %
            % ## Input
            % * __corners__ Cell-array of top-left corners `{[x,y],...}`.
            %
            % ## Output
            % * __tl__ top-left corner `[x,y]`
            %
            % See also: cv.Blender.resultRoi
            %
            tl = Blender_(0, 'resultTl', corners);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [ comp ] = getHypInfo( hyps )
%GETHYPINFO Summary of this function goes here
%   Detailed explanation goes here
num_hyps = length(hyps);
comp = repmat(struct('objcenter',[],'objsize',[],'objtype',[]), num_hyps, 1 );
for i = 1:num_hyps
    if isempty(hyps{i})
        continue;
    end
    T = reshape(vertcat(hyps{i}.align)', 24, []);
    comp(i).objcenter = (T(1:3,:)+T(19:21,:))/2;
    comp(i).objcenter = comp(i).objcenter - ...
        repmat(comp(i).objcenter(:,1), 1, size(comp(i).objcenter, 2));
    comp(i).objsize = T(19:21,:)-T(1:3,:)+0.1;
    comp(i).objtype = horzcat(hyps{i}.objtype);
    comp(i).objaspectXY = comp(i).objsize(2,:)./comp(i).objsize(1,:);
    comp(i).objregionXY = comp(i).objsize(2,:).*comp(i).objsize(1,:);

    comp(i).objcenter = single(comp(i).objcenter);
    comp(i).objsize = single(comp(i).objsize);
    comp(i).objtype = single(comp(i).objtype);
    comp(i).objaspectXY = single(comp(i).objaspectXY);
    comp(i).objregionXY = single(comp(i).objregionXY);
end

end

author: zouchuhang

repo: LayoutNet


text:
function q = reshape(q,varargin)
% overloads reshape

q.a = reshape(q.a,varargin{:});
q.b = reshape(q.b,varargin{:});
q.c =  reshape(q.c,varargin{:});
q.d = reshape(q.d,varargin{:});

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
		% Determinarea domeniului de definitie al functiei
t=-10:10;
		% Apelarea functiei MATLAB def_funct
f=def_funct(t,0,3,5,5,1)

author: Sable

repo: mcbench-benchmarks


text:
function a = uplus(a)
%UPLUS        Gradient unary plus  + a
%

% written  10/16/98     S.M. Rump
% modified 12/18/02     S.M. Rump  performance improved
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

repo: OpenMAS


text:
function plotData(x, y)
%PLOTDATA Plots the data points x and y into a new figure
%   PLOTDATA(x,y) plots the data points and gives the figure axes labels of
%   population and profit.

% ====================== YOUR CODE HERE ======================
% Instructions: Plot the training data into a figure using the
%               "figure" and "plot" commands. Set the axes labels using
%               the "xlabel" and "ylabel" commands. Assume the
%               population and revenue data have been passed in
%               as the x and y arguments of this function.
%
% Hint: You can use the 'rx' option with plot to have the markers
%       appear as red crosses. Furthermore, you can make the
%       markers larger by using plot(..., 'rx', 'MarkerSize', 10);

figure; % open a new figure window

plot(x, y, 'rx', 'MarkerSize', 10);      % Plot the data
ylabel('Profit in $10,000s');            % Set the y-axis label
xlabel('Population of City in 10,000s'); % Set the x-axis label

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

end

repo: machine-learning-coursera


text:
function [R, scale]=ARFIT_arqr(v, p, mcor)
%ARQR	QR factorization for least squares estimation of AR model.
%
%  [R, SCALE]=ARQR(v,p,mcor) computes the QR factorization needed in
%  the least squares estimation of parameters of an AR(p) model. If
%  the input flag mcor equals one, a vector of intercept terms is
%  being fitted. If mcor equals zero, the process v is assumed to have
%  mean zero. The output argument R is the upper triangular matrix
%  appearing in the QR factorization of the AR model, and SCALE is a
%  vector of scaling factors used to regularize the QR factorization.
%
%  ARQR is called by ARFIT.
%
%  See also ARFIT.

%  Modified 29-Dec-99
%           24-Oct-10 Tim Mullen (added support for multiple realizatons)
%
%  Author: Tapio Schneider
%          tapio@gps.caltech.edu

  % n:   number of time steps (per realization)
  % m:   number of variables (dimension of state vectors)
  % ntr: number of realizations (trials)
  [n,m,ntr] = size(v);

  ne    = ntr*(n-p);            % number of block equations of size m
  np    = m*p+mcor;             % number of parameter vectors of size m

  % If the intercept vector w is to be fitted, least squares (LS)
  % estimation proceeds by solving the normal equations for the linear
  % regression model
  %
  %                  v(k,:)' = Aaug*u(k,:)' + noise(C)        (1)
  %
  % with Aaug=[w A] and `predictors'
  %
  %              u(k,:) = [1 v(k-1,:) ...  v(k-p,:)].         (2a)
  %
  % If the process mean is taken to be zero, the augmented coefficient
  % matrix is Aaug=A, and the regression model
  %
  %                u(k,:) = [v(k-1,:) ...  v(k-p,:)]          (2b)
  %
  % is fitted.
  % The number np is the dimension of the `predictors' u(k).
  %
  % If multiple realizations are given (ntr > 1), they are appended
  % as additional ntr-1 blocks of rows in the normal equations (1), and
  % the 'predictors' (2) correspondingly acquire additional row blocks.

  % Initialize the data matrix K (of which a QR factorization will be computed)
  K = zeros(ne,np+m);                 % initialize K
  if (mcor == 1)
    % first column of K consists of ones for estimation of intercept vector w
    K(:,1) = ones(ne,1);
  end

  % Assemble `predictors' u in K
  for itr=1:ntr
    for j=1:p
      K((n-p)*(itr-1) + 1 : (n-p)*itr, mcor+m*(j-1)+1 : mcor+m*j) = ...
          squeeze(v(p-j+1:n-j, :, itr));
    end
    % Add `observations' v (left hand side of regression model) to K
    K((n-p)*(itr-1) + 1 : (n-p)*itr, np+1 : np+m) = squeeze(v(p+1:n, :, itr));
  end

  % Compute regularized QR factorization of K: The regularization
  % parameter delta is chosen according to Higham's (1996) Theorem
  % 10.7 on the stability of a Cholesky factorization. Replace the
  % regularization parameter delta below by a parameter that depends
  % on the observational error if the observational error dominates
  % the rounding error (cf. Neumaier, A. and T. Schneider, 2001:
  % "Estimation of parameters and eigenmodes of multivariate
  % autoregressive models", ACM Trans. Math. Softw., 27, 27--57.).
  q     = np + m;             % number of columns of K
  delta = (q^2 + q + 1)*eps;  % Higham's choice for a Cholesky factorization
  scale = sqrt(delta)*sqrt(sum(K.^2));
  R     = triu(qr([K; diag(scale)]));

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% Compare the speeds of various inference engines on a coupled HMM

N = 2;
Q = 2;
rand('state', 0);
randn('state', 0);
discrete = 1;
if discrete
  Y = 2; % size of output alphabet
else
  Y = 1;
end
coupled = 1;
[bnet, onodes] = mk_chmm(N, Q, Y, discrete, coupled);
ss = N*2;

T = 3;

engine = {};
tic; engine{end+1} = jtree_dbn_inf_engine(bnet, 'observed', onodes);  toc
%tic; engine{end+1} = jtree_ndxSD_dbn_inf_engine(bnet, onodes);  toc
%tic; engine{end+1} = jtree_ndxB_dbn_inf_engine(bnet, onodes);  toc
engine{end+1} = hmm_inf_engine(bnet, onodes);
%engine{end+1} = dhmm_inf_engine(bnet, onodes);
tic; engine{end+1} = jtree_unrolled_dbn_inf_engine(bnet, T, onodes);  toc

%engine{end+1} = bk_inf_engine(bnet, 'ff', onodes);
%engine{end+1} = loopy_dbn_inf_engine(bnet, onodes);

exact = [1 2 3];

filter = 0;
single = 0;
maximize = 0;

[err, time, engine] = cmp_inference(bnet, onodes, engine, exact, T, filter, single, maximize);
%err = cmp_learning(bnet, onodes, engine, exact, T);


text:
function [h] = fdr(p, q)

% FDR false discovery rate
%
% Use as
%   h = fdr(p, q)
%
% This implements
%   Genovese CR, Lazar NA, Nichols T.
%   Thresholding of statistical maps in functional neuroimaging using the false discovery rate.
%   Neuroimage. 2002 Apr;15(4):870-8.
%
% There are two types of FDR correction (Benjamini-Hochberg & Benjamini-Yekutieli), of
% which the second is currently implemented.

% Copyright (C) 2005-2015, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% 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$

% convert the input into a row vector
dim = size(p);
p = reshape(p, 1, numel(p));

% sort the observed uncorrected probabilities
[ps, indx] = sort(p);

% count the number of voxels
V = length(p);

% compute the threshold probability for each voxel
pi = ((1:V)/V)  * q / c(V);

if any(ps<=pi)
  h = ps<=(max(ps(ps<=pi)));
else
  h = false(size(ps));
end

% undo the sorting
[dum, unsort] = sort(indx);
h = h(unsort);

% convert the output back into the original format
h = reshape(h, dim);

function s = c(V)
% See Genovese, Lazar and Holmes (2002) page 872, second column, first paragraph
if V<1000
  % compute it exactly
  s = sum(1./(1:V));
else
  % approximate it
  s = log(V) + 0.57721566490153286060651209008240243104215933593992359880576723488486772677766467093694706329174674951463144724980708248096050401448654283622417399764492353625350033374293733773767394279259525824709491600873520394816567;
end


text:
clear all;

fprintf('test Lasso weighted\n');
randn('seed',0);
% Data are generated
X=randn(64,10000);
X=X./repmat(sqrt(sum(X.^2)),[size(X,1) 1]);
D=randn(64,256);
D=D./repmat(sqrt(sum(D.^2)),[size(D,1) 1]);

% parameter of the optimization procedure are chosen
param.L=20; % not more than 20 non-zeros coefficients (default: min(size(D,1),size(D,2)))
param.lambda=0.15; % not more than 20 non-zeros coefficients
param.numThreads=8; % number of processors/cores to use; the default choice is -1
                    % and uses all the cores of the machine
param.mode=2;       % penalized formulation

W=rand(size(D,2),size(X,2));

tic
alpha=mexLassoWeighted(X,D,W,param);
t=toc;
toc

fprintf('%f signals processed per second\n',size(X,2)/t);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function parFormatDescription
% Description of .par file format for specifying event-related or block-
% design stimuli.
%
% The format is as follows: a tab-delimited ASCII text file with the
% following columns:
%
% (1) [\tab] (2) [\tab] (3) [\tab] (4) [\tab] (5) [\tab] (6) [\tab] (7)
%
% where: (1) is the onset time in seconds;
%        (2) is an integer specifying an event type, with 0 specifying
%            a 'baseline' or 'null' condition;
%        (3) is a text label specifying the name of a condition. (Only
%            the first occurrence of a condition is parsed to get the
%            name);
%        (4) specifies a color to associate w/ each condition, for use
%            in visualizing data that is sorted by condition. The color
%            is specified in matlab format, either a single letter
%            (e.g. 'w' for 'white', 'r' for 'red', 'm' for 'magenta')
%            or the format '[R G B]' where R, G, and B are the relative
%            intensities of the red, green and blue channels, each ranging
%            from 0-1 (see 'help plot');
%        (5) specifies a path to an image to associate w/ each condition,
%            for use in visualizing data that is sorted by condition;
%        (6) and (7) specify behavioral data, such as response times
%            or correct/incorrect responses.
%
%   All columns except for (1) and (2) are optional. If you don't have
%   information for a column, but would like to include it for a later
%   column, you can add the text '' or [] for that column (or add two
%   tabs).
%
%
%
% Importantly, .par files here have an implicit assumption that
% there's always _something_ going on during an experiment: events don't
% have a separately-specified duration, but are assumed to continue
% until the next event starts. Therefore, you may want to specify
% explicitly the end of a scan as a separate event (which can be
% assigned to a null condition or a negative number, if you like).

% if someone called this as a function, just give them the info:
p = fileparts(mfilename);
web(fullfile(p,'parFormatDescription.txt'));

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function filterS = runFilter(filterS, planDB)
%"runFilter"
%   Given a single filter struct and a planDB, return filterS updated with
%filterS.indicies containing indicies of plans that match the filter's
%parameters.  A filter struct filterS consists of the fields:
%        Filters   - the user set string "name" of the filter
%        fieldname - the field to be operated on
%        action    - the filtering action ('contains' 'exists' 'isempty')
%        regexp    - regular expression to apply to filter.
%        bool      - AND/OR
%        invert    - 0 or 1, 1 indicates NOT
%        indices   - indicies of plans that match the filter.
%
%   By JRA 12/15/03
%
%   filterS     :    filter struct
%   planDB      :    plan database
%
%   filterS     :    filter struct with indices updated
%
% Usage:
%   filterS = runFilter(filterS, planDB)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

fieldname = filterS.fieldname;
allFieldNames = {planDB.fieldIndex.fieldname};
isPlanC = [planDB.matFiles.isPlanC];
indices = repmat(logical(0), [1,length(planDB.matFiles)]);

%Get list of plans containing the requested field.
fIndex = find(strcmpi(allFieldNames, fieldname));
plansContaining = planDB.fieldIndex(fIndex).planIndices;

switch lower(filterS.action)
    case 'contains'
        for i = 1:length(plansContaining)
            data = getFieldContents(filterS.fieldname, planDB.matFiles(plansContaining(i)).extract);
            if matches(data, filterS.regexp)
                indices(plansContaining(i)) = 1;
            end
        end
        %     case 'does not contain'
        %         for i = 1:length(plansContaining)
        %             data = getFieldContents(filterS.fieldname, planDB.matFiles(plansContaining(i)).extract);
        %             if ~matches(data, filterS.regexp)
        %                indices(end+1)= plansContaining(i);
        %             end
        %         end
    case 'exists'
        indices(plansContaining) = logical(1);
        %     case 'does not exist'
        %         indices = isPlanC;
        %         indices(plansContaining) = 0;
        %         indices = find(indices);
    case 'is empty'
        for i = 1:length(plansContaining)
            data = getFieldContents(filterS.fieldname, planDB.matFiles(plansContaining(i)).extract);
            if isempty([data{:}])
                indices(plansContaining(i)) = 1;
            end
        end
        %     case 'is not empty'
        %         for i = 1:length(plansContaining)
        %             data = getFieldContents(filterS.fieldname, planDB.matFiles(plansContaining(i)).extract);
        %             if ~isempty([data{:}])
        %                indices(end+1)= plansContaining(i);
        %             end
        %         end
end

if filterS.invert
    indices = isPlanC & ~indices;
end
filterS.indices = find(indices);

end

function bool = matches(data, regexpS)
%Check of a piece of data matches a regexp. Passive, assumes false until
%proven otherwise.
bool = 0;
data = data(:);
for i=1:length(data)
    switch class(data{i})
        case {'cell', 'struct'}
            continue;
        case 'char'
            if regexpi(data{i}, regexpS)
                bool = 1;
            end
        otherwise
            if ~isnumeric(data{i})
                continue;
            else
                result = regexpi(num2str(data{i}), regexpS);
                for i=1:length(result)
                    if iscell(result)
                        if ~isempty(result{:})
                            bool = 1;
                        end
                    else
                        if ~isempty(result)
                            bool = 1;
                        end
                    end
                end
            end
    end
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [u,w,AE,AI,isExteriorTElem,isExteriorSElem] = interfacefittedPoisson(node,telem,selem,pde,interfaceEdge,bdEdge,option)
%% INTERFACEFITTEDPOISSON Poisson equation: P1 linear element.
%
%   u = INTERFACEPOISSON(node,telem,selem,pde,E) produces the linear finite element
%   approximation of the interface Poisson equation
%
%   input:
%       node, N*2 matrix, node(i,:) are the xy coordinates of i-th node;
%       telem, NTT*3 matrix, telem(j,:) are the three  global indices of the
%              vertices ofj-th triangle element;
%       selem,NTS*4 matrix, selem(j,:) are the four global indices of the
%             vertices of j-th quad element;
%
%

[N,Dim] = size(node);
Ndof = N;

%% Geometry structures of interface meshes

NTS = size(selem,1);
isExteriorSElem = false(NTS,1);
c = (node(selem(:,1),:) + node(selem(:,2),:) + node(selem(:,3),:)+node(selem(:,4),:))/4;
isExteriorSElem(pde.phi(c)>0) = true;

NTT = size(telem,1);
isExteriorTElem = false(NTT,1);
c = (node(telem(:,1),:) + node(telem(:,2),:) + node(telem(:,3),:))/3;
isExteriorTElem(pde.phi(c)>0) = true;

[sAE,sbE] = getstiffmatrixandrhsonquad(node,selem(isExteriorSElem,:));
[sAI,sbI] = getstiffmatrixandrhsonquad(node,selem(~isExteriorSElem,:));
[tAE,tbE] = getstiffmatrixandrhsontri(node,telem(isExteriorTElem,:));
[tAI,tbI] = getstiffmatrixandrhsontri(node,telem(~isExteriorTElem,:));

A = sAE + sAI + tAE + tAI;
b = sbE + sbI + tbE + tbI;

AI = sAI + tAI;
AE = sAE + tAE;

flux = getfluxconditiononinterface(node,interfaceEdge);
b = b - flux;

%% Extend w to the whole domain
isInterfaceNode = false(Ndof,1);
isInterfaceNode(interfaceEdge(:)) = true;
interfaceNode = find(isInterfaceNode);

isInNode = false(Ndof,1);
interiorSElem = selem(~isExteriorSElem,:);
interiorTElem = telem(~isExteriorTElem,:);
isInNode([interiorSElem(:);interiorTElem(:)]) = true;
inNode = find(isInNode & ~isInterfaceNode);

w = zeros(Ndof,1);
FI= zeros(Ndof,1);
w(interfaceNode) = pde.exactw(node(interfaceNode,:));
FI = FI - AI*w;
extensionoption.printlevel = 0;
w(inNode) = amg(AI(inNode, inNode),FI(inNode),extensionoption);

%% Dirichlet boundary condition
isBdNode = false(Ndof,1);
isBdNode(bdEdge(:)) = true;
bdNode = find(isBdNode);

u = zeros(Ndof,1);
u(bdNode) = pde.g_D(node(bdNode,:));
b = b - A*u+AI*w;
b(bdNode) = u(bdNode);

freeNode = find(~isBdNode);
bdidx = zeros(Ndof,1);
bdidx(bdNode) = 1;
Tbd = spdiags(bdidx,0,Ndof,Ndof);
T = spdiags(1-bdidx,0,Ndof,Ndof);
AD = T*A*T + Tbd;

%% Solve the system of linear equations
if isempty(freeNode), return; end
% Set up solver type
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if Ndof <= 1e3  % Direct solver for small size systems
        option.solver = 'direct';
    else            % MGCG  solver for large size systems
        option.solver = 'mg';
    end
end
solver = option.solver;
% solve
switch solver
    case 'direct'
        tic;
        u(freeNode) = AD(freeNode,freeNode)\b(freeNode);
        residual = norm(b - AD*u);
        info = struct('solverTime',toc,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'none'
        info = struct('solverTime',[],'itStep',0,'err',[],'flag',3,'stopErr',[]);
    case 'amg'
        option.solver = 'CG';
        [u(freeNode),info] = amg(AD(freeNode,freeNode),b(freeNode),option);
end

%%
eqn = struct('A',AD,'b',b,'freeNode',freeNode);

   %% Assemble stiffness matrix and rhs on quad mesh
    function [A,b] = getstiffmatrixandrhsonquad(node,elem)
        [NT,NV] = size(elem);
        % generate sparse pattern
        ii = zeros(10*NT,1); jj = zeros(10*NT,1);
        index = 0;
        for i = 1:4
            for j = i:4
                ii(index+1:index+NT) = double(elem(:,i));
                jj(index+1:index+NT) = double(elem(:,j));
                index = index + NT;
            end
        end
        % quadrature points
        if ~isfield(pde,'d'), pde.d = []; end
        if ~isfield(option,'dquadorder')
            option.dquadorder = 2;        % default order is exact for quadratic function
        end
        [pts, weight] = quadptsquad(option.dquadorder);
        nQuad = size(pts,1);
        % compute non-zeros
        sA = zeros(10*NT,nQuad);
        for p = 1:nQuad
            % Dphi at quadrature points
            [phi, Dphip, J] = quadbasis(node,elem,pts(p,:));
            index = 0;
            for i = 1:4
                for j = i:4
                    Aij = 0;
                    if isempty(pde.d) || isnumeric(pde.d)
                        Aij = Aij + weight(p)*dot(Dphip(:,:,i),Dphip(:,:,j),2);
                    else
                        pxy = zeros(NT, Dim);
                        for ip = 1:Dim
                            xi = node(:,ip);
                            pxy(:,ip) = xi(elem)*phi;
                        end
                        Aij = Aij + weight(p)*dot(Dphip(:,:,i),Dphip(:,:,j),2).*pde.d(pxy);
                    end
                    if ~isempty(pde.d) && isnumeric(pde.d) % d is piecewise constant
                        Aij = pde.d.*Aij;
                    end
                    Aij = Aij.*J;
                    sA(index+1:index+NT,p) = Aij;
                    index = index + NT;
                end
            end
        end
        sA = sum(sA,2);
        % assemble the matrix
        diagIdx = (ii == jj);   upperIdx = ~diagIdx;
        A = sparse(ii(diagIdx),jj(diagIdx),sA(diagIdx),Ndof,Ndof);
        AU = sparse(ii(upperIdx),jj(upperIdx),sA(upperIdx),Ndof,Ndof);
        A = A + AU + AU';
        clear Aij ii jj Dphip

        %% Assemble the right hand side
        b = zeros(Ndof,1);
        if ~isfield(option,'fquadorder')
            option.fquadorder = 3;   % default order
        end
        if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
            pde.f = [];
        end

        if ~isempty(pde.f)
            [pts, weight] = quadptsquad(option.fquadorder);
            nQuad = size(pts,1);
            bt = zeros(NT,NV);
            for p = 1:nQuad
                % quadrature points in the x-y coordinate
                [phi, ~, J] = quadbasis(node,elem, pts(p,:));
                pxy = zeros(NT, Dim);
                for i = 1:Dim
                    xi = node(:,i);
                    pxy(:,i) = xi(elem)*phi; % ? questionable
                end
                fp = pde.f(pxy);
                bt = bt + (weight(p)*fp.*J)*phi';
            end
            b = accumarray(elem(:),bt(:),[Ndof 1]);
        end

    end

    %% Assemble the stiffmatrix and right hand side on triangle mesh
    function [A,b] = getstiffmatrixandrhsontri(node,elem)
        NT = size(elem,1);
        % quadrature points
        center = (node(elem(:,1),:) + node(elem(:,2),:) + node(elem(:,3),:))/3;
       % Diffusion coefficient
        if isfield(pde,'d') && ~isempty(pde.d)
            if isnumeric(pde.d)
                K = pde.d;                   % d is an array
            else                            % d is a function
                K = pde.d(center);
            end
        else
            K = [];
        end
        [Dphi,area] = gradbasis(node,elem);
        A = sparse(Ndof,Ndof);
        for i = 1:3
            for j = i:3
                Aij = (Dphi(:,1,i).*Dphi(:,1,j)+Dphi(:,2,i).*Dphi(:,2,j)).*area;
                if ~isempty(K)
                    Aij = K.*Aij;
                end
                if (j==i)
                    A = A + sparse(elem(:,i),elem(:,j),Aij,Ndof,N);
                else
                    A = A + sparse([elem(:,i);elem(:,j)],[elem(:,j);elem(:,i)],[Aij; Aij],Ndof,Ndof);
                end
            end
        end

        b = zeros(Ndof,1);
        if ~isfield(option,'fquadorder')
            option.fquadorder = 3;   % default order
        end
        if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
            pde.f = [];
        end
        if ~isempty(pde.f)
            [lambda,weight] = quadpts(option.fquadorder);
            nQuad = size(lambda,1);
            ft = zeros(NT,3);
            for p = 1:nQuad
                % quadrature points in the x-y coordinate
                pxy = lambda(p,1)*node(elem(:,1),:) ...
                    + lambda(p,2)*node(elem(:,2),:) ...
                    + lambda(p,3)*node(elem(:,3),:);
                % function values at quadrature points
                fp = pde.f(pxy);
                % evaluate fp outside.
                for j = 1:3
                    ft(:,j) = ft(:,j) + lambda(p,j)*weight(p)*fp;
                end
            end
            ft = ft.*[area,area,area];
            b = accumarray(elem(:),ft(:),[Ndof 1]);
        end
    end

    %% Neumann boundary condition on interface edges
    function b = getfluxconditiononinterface(node,interfaceEdge)
        ve = node(interfaceEdge(:,1),:) - node(interfaceEdge(:,2),:);
        ve = [-ve(:,2), ve(:,1)];
        edgeLen = sqrt(sum(ve.^2,2));
        if ~isfield(option,'gNquadorder')
            option.gNquadorder = 5;   % default order
        end
        [lambda,weight] = quadpts1(option.gNquadorder);
        nQuad = length(weight);
        ge = zeros(size(interfaceEdge,1),2);
        for i = 1:nQuad
            pxy=lambda(i,1)*node(interfaceEdge(:,1),:)+lambda(i,2)*node(interfaceEdge(:,2),:);
            leftpt = pxy - ve/2;
            rightpt = pxy + ve/2;
            pxy = findintersectbisect(pde.phi,leftpt,rightpt);
            q = pde.exactq(pxy);
            ge(:,1)=ge(:,1)+weight(i)*lambda(i,1)*q;
            ge(:,2)=ge(:,2)+weight(i)*lambda(i,2)*q;
        end
        ge = ge.*[edgeLen,edgeLen];
        b = accumarray(interfaceEdge(:), ge(:),[Ndof,1]);
    end

end

repo: ifem


text:
classdef VideoMaker < handle

    properties (GetAccess = private, SetAccess = public)
        iterations
    end

    properties (Access = private)
        gidPath
        tclFileName
        filesFolder
        tclFileWriter
        tclTemplateName
        fieldName
        outputName
    end

    properties (Access = private)
        outPutNames
        fileName
        designVariableName
        tclTemplateNames
    end

    methods (Access = public)

        function makeVideo(obj)
            obj.makeDesignVariableVideo();
            obj.makeRegularizedDesignVariableVideo();
        end

    end

    methods (Access = public, Static)

        function obj = create(cParams)
            f = VideoMakerFactory();
            obj = f.create(cParams);
        end

    end

    methods (Access = protected)

        function init(obj,cParams)
            obj.fileName           = cParams.caseFileName;
            obj.designVariableName = cParams.designVarType;
            obj.tclTemplateNames   = cParams.tclTemplateNames;
            obj.outPutNames        = cParams.outPutNames;
            obj.createPaths();
            obj.createFolder();
        end

    end

    methods (Access = private)

        function createFolder(obj)
            if ~exist(obj.filesFolder,'dir')
                mkdir(obj.filesFolder);
            end
        end

        function createPaths(obj)
            %obj.gidPath = '/home/alex/GiDx64/14.1.9d';
            obj.gidPath = '/home/alex/GiDx64/gid-15.0.3';
            obj.filesFolder = fullfile(pwd,'Output',obj.fileName);
        end

        function makeDesignVariableVideo(obj)
            for ivideo = 1:numel(obj.tclTemplateNames)
                obj.tclTemplateName = obj.tclTemplateNames{ivideo};
                obj.fieldName = obj.designVariableName;
                obj.outputName = [obj.fileName,obj.outPutNames{ivideo}];
                obj.makeFieldVideo();
            end
        end

        function makeRegularizedDesignVariableVideo(obj)
            obj.fieldName = 'DensityGauss';
            obj.tclTemplateName = 'Make_Video_density';
            obj.outputName = [obj.fileName,obj.fieldName];
            obj.makeFieldVideo();
        end

        function makeFieldVideo(obj)
            obj.createTclFileName();
            obj.createTclFileWriter();
            obj.writeTclFile();
            obj.executeTclFile();
            obj.deleteTclFile();
        end

        function createTclFileName(obj)
            fName = 'tcl_gid.tcl';
            obj.tclFileName = fullfile(obj.filesFolder,fName);
        end

        function createTclFileWriter(obj)
            s.type            = obj.fieldName;
            s.tclFileName     = obj.tclFileName;
            s.filesFolder     = obj.filesFolder;
            s.outputName      = obj.outputName;
            s.iterations      = obj.iterations;
            s.fileName        = obj.fileName;
            s.tclTemplateName = obj.tclTemplateName;
            obj.tclFileWriter = TclFileWriter(s);
        end

        function writeTclFile(obj)
            obj.tclFileWriter.write();
        end

        function executeTclFile(obj)
            tFile = replace(obj.tclFileName,'\','\\');
            gFile = fullfile(obj.gidPath,'gid_offscreen');
            executingLine = ['"',gFile,'"',' -offscreen', ' -t ' ,'"source ',tFile,'"'];
            system(executingLine);
        end

        function deleteTclFile(obj)
            tFile = obj.tclFileName;
            if ispc
                system(['DEL ',tFile]);
            elseif isunix
                system(['rm ',tFile]);
            elseif ismac
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function err = labelic(subject,apslice)
% err = labelic(subject,apslice)
%
% Simple program that attempts to label the internal capsule based
% upon the automatic segmentation. It is very crude, to say the
% least, but it might serve as a good starting point.
%
% Uses SUBJECTS_DIR/subject/mri/aseg.mgz as input.
%
% Writes SUBJECTS_DIR/subject/mri/asegic.mgz as output using
% IC Labels are as found in FreeSurferColors.txt
%  155  Left-IntCapsule-Ant
%  156  Right-IntCapsule-Ant
%  157  Left-IntCapsule-Pos
%  158  Right-IntCapsule-Pos
%
% Anterior/Posterior is defined by apslice, ie, voxels anterior to
% coronal apslice are labeled anterior, etc. If unspecified, apslice
% defaults to 128.
%
% Also writes SUBJECTS_DIR/subject/mri/icmask.mgz as output. This
% is a binary mask of what the full IC label and may be helpful
% when editing the asegic.mgz.
%
%
% tkmedit subject orig.mgz -aux icmask.mgh -segmentation ./asegic.mgz
%
%

%
% labelic.m
%
% Original Author: Doug Greve
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:12 $
%    $Revision: 1.3 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

ndil = 6;

err = 1;
if(nargin < 1 | nargin > 2)
  fprintf('err = labelic(subject,apslice)\n');
  return;
end

if(~exist('apslice','var')) apslice = []; end
if(isempty(apslice)) apslice = 128; end

fprintf('Starting labelic\n');
fprintf('apslice = %d\n',apslice);

SUBJECTS_DIR = deblank(getenv('SUBJECTS_DIR'));

fspec = sprintf('%s/%s/mri/aseg.mgz',SUBJECTS_DIR,subject);
aseg = MRIread(fspec);
if(isempty(aseg)) return; end

% First, Create a binary mask of these structures
% Putamen
% Caudate
% Pallidum
% Lateral-Ventricle
% Thalamus-Proper
% VentralDC 28 60
m0 = (aseg.vol == 10 | aseg.vol == 11 | ...
      aseg.vol == 12 | aseg.vol == 13 | ...
      aseg.vol ==  4 | aseg.vol == 43 | ...
      aseg.vol == 49 | aseg.vol == 50 | ...
      aseg.vol == 51 | aseg.vol == 52 | ...
      aseg.vol == 28 | aseg.vol == 60);

% Dilate the mask. Needs to dilate enough to bridge the gap across
% the IC.
fprintf('Dilating %d \n',ndil);
tic;
md0 = fast_dilate(m0,ndil);
fprintf('  done %g\n',toc);

% Now erode by the same amount. Wont erode where mask bridged the
% IC gap.
fprintf('Eroding %d \n',ndil);
tic;
md = fast_dilate(md0,ndil,1);
fprintf('  done %g\n',toc);

wmleft  = (aseg.vol == 2);    % Left WM
icleft  = md & ~m0 & wmleft;  % Left IC
indicleft = find(icleft);
[r c s] = ind2sub(size(icleft),indicleft);
iant = find(s >  apslice);
ipos = find(s <= apslice);
indicleftant = indicleft(iant);
indicleftpos = indicleft(ipos);
aseg.vol(indicleftant) = 155;
aseg.vol(indicleftpos) = 157;

wmright = (aseg.vol == 41);    % Right WM
icright  = md & ~m0 & wmright; % Right IC
indicright = find(icright);
[r c s] = ind2sub(size(icright),indicright);
iant = find(s >  apslice);
ipos = find(s <= apslice);
indicrightant = indicright(iant);
indicrightpos = indicright(ipos);
aseg.vol(indicrightant) = 156;
aseg.vol(indicrightpos) = 158;

fspec = sprintf('%s/%s/mri/asegic.mgz',SUBJECTS_DIR,subject);
MRIwrite(aseg,fspec);

% Create the IC Mask
icmask = aseg;
icmask.vol = icleft | icright;
fspec = sprintf('%s/%s/mri/icmask.mgz',SUBJECTS_DIR,subject);
MRIwrite(icmask,fspec);

err = 0;

return;


text:
function lines = readLines( fileName )
%readLines reads the file line by line
%
% lines = readLines( fileName );
%
% Input:
%   fileName - string containing the full file name
%
% Output:
%   lines - cell array of lines contained in the file

fileID = fopen(fileName, 'r');
if fileID == -1
    error(['File ', fileName, ' can not be opened!']);
end
lines = textscan(fileID, '%s\n');
fclose(fileID);

lines = lines{1};

end

author: aosokin


text:
function x = r8sto_yw_sl ( n, b, x )

%*****************************************************************************80
%
%% R8STO_YW_SL solves the Yule-Walker equations for a R8STO matrix.
%
%  Discussion:
%
%    The R8STO storage format is used for a symmetric Toeplitz matrix.
%    It stores the N elements of the first row.
%
%    The matrix is also required to be positive definite.
%
%    This implementation of the algorithm assumes that the diagonal element
%    is 1.
%
%    The real symmetric Toeplitz matrix can be described by N numbers, which,
%    for convenience, we will label B(0:N-1).  We assume there is one more
%    number, B(N).  If we let A be the symmetric Toeplitz matrix whose first
%    row is B(0:N-1), then the Yule-Walker equations are:
%
%      A * X = -B(1:N)
%
%  Example:
%
%    To solve
%
%     1.0 0.5 0.2    x1   0.5
%     0.5 1.0 0.5 *  x2 = 0.2
%     0.2 0.5 1.0    x3   0.1
%
%    we input:
%
%      N = 3
%      B(1:3) = (/ 0.5, 0.2, 0.1 /)
%
%    with output:
%
%      X(1:3) = (/ -75, 12, -5 /) / 140
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gene Golub and Charles Van Loan,
%    Section 4.7.2, "Solving the Yule-Walker Equations",
%    Matrix Computations,
%    Third Edition,
%    Johns Hopkins, 1996.
%
%  Parameters:
%
%    Input, integer N, the order of the system.
%
%    Input, real B(N), defines the linear system.  The first entry of A is
%    a 1, followed by B(1) through B(N-1).  The right hand side of the
%    system is -B(1:N).
%
%    Output, real X(N), the solution of the linear system.
%
  x(1) = -b(1);
  beta = 1.0E+00;
  alpha = -b(1);

  for i = 1 : n-1
    beta = ( 1.0E+00 - alpha * alpha ) * beta;
    alpha = - ( b(i+1) + b(i:-1:1) * x(1:i)' ) / beta;
    x(1:i) = x(1:i) + alpha * x(i:-1:1);
    x(i+1) = alpha;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function  [C,phi,S12,confC,phierr,Cerr]=cohmathelper(J,err,Nsp)
% Helper function called by coherency matrix computations.
%
% Usage: [C,phi,S12,confC,phierr,Cerr]=cohmathelper(J,err,Nsp)
% Inputs:
% J   : Fourier transforms of data
% err : [0 p] or 0 for no errors; [1 p] for theoretical confidence level,
%       [2 p] for Jackknife (p - p value)
% Nsp : pass the number of spikes in each channel if finite size corrections are desired
%
% Outputs:
%
% C   : coherence
% phi : phase of coherency
% S12 : cross spectral matrix
% confC : confidence level for coherency - only for err(1)>=1
%       phierr - standard deviation for phi (note that the routine gives phierr as phierr(1,...)
%                and phierr(2,...) in order to incorporate Jackknife (eventually).
%                Currently phierr(1,...)=phierr(2,...). Note that phi + 2 phierr(1,...) and phi -2
%                phierr(2,...) will give 95% confidence bands for phi - only for err(1)>=1
% Cerr  : error bars for coherency (only for Jackknife estimates)-only for err(1)=2
%

errtype=err(1);
trialave=0;
[nf,K,Ch]=size(J);
clear K
confC=zeros(Ch,Ch);
C=zeros(nf,Ch,Ch);
S12=zeros(nf,Ch,Ch);
phi=zeros(nf,Ch,Ch);
phierr=zeros(2,nf,Ch,Ch);
if errtype==2; Cerr=zeros(2,nf,Ch,Ch);end;

for ch1=1:Ch;
     J1=squeeze(J(:,:,ch1));
     C(1:nf,ch1,ch1)=1;
     phi(1:nf,ch1,ch1)=0;
%      if errtype==2;
%           phierr(1:nf,ch1,ch1)=0;
%           Cerr(1:2,1:nf,ch1,ch1)=0;
%      elseif errtype==1
%            phierr(1:2,1:nf,ch1,ch1)=0;
%      end;
     s1=squeeze(mean(conj(J1).*J1,2));
     for ch2=1:ch1-1;
          J2=squeeze(J(:,:,ch2));
          s12=squeeze(mean(conj(J1).*J2,2));
          s2=squeeze(mean(conj(J2).*J2,2));
          C12=s12./sqrt(s1.*s2);
          C(:,ch1,ch2)=abs(C12);
          C(:,ch2,ch1)=C(:,ch1,ch2);
          phi(:,ch1,ch2)=angle(C12);
          phi(:,ch2,ch1)=phi(:,ch1,ch2);
          S12(:,ch1,ch2)=s12;
          S12(:,ch2,ch1)=S12(:,ch1,ch2);
          if errtype==2
             if nargin<3;
                 [conf,phie,Ce]=coherr(abs(C12),J1,J2,err,trialave);
             else
                 [conf,phie,Ce]=coherr(abs(C12),J1,J2,err,trialave,Nsp(ch1),Nsp(ch2));
             end
             confC(ch1,ch2)=conf;
             phierr(1,:,ch1,ch2)=phie;phierr(2,:,ch1,ch2)=phie;
             Cerr(1,:,ch1,ch2)=Ce(1,:);
             Cerr(2,:,ch1,ch2)=Ce(2,:);
             confC(ch2,ch1)=conf;
             phierr(1,:,ch2,ch1)=phie;phierr(2,:,ch2,ch1)=phie;
             Cerr(:,:,ch2,ch1)=Ce;
          elseif errtype==1
             if nargin<3;
                 [conf,phie]=coherr(abs(C12),J1,J2,err,trialave);
             else
                 [conf,phie]=coherr(abs(C12),J1,J2,err,trialave,Nsp(ch1),Nsp(ch2));
             end
             confC(ch1,ch2)=conf;
             phierr(1,:,ch1,ch2)=phie;phierr(2,:,ch1,ch2)=phie;
             confC(ch2,ch1)=conf;
             phierr(1,:,ch2,ch1)=phie;phierr(2,:,ch2,ch1)=phie;
          end;
     end;
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% STD_BUILDDESIGNMAT - Build the design matrix for an specific design
% specified in the structure ''design'' provided as input
%
% Usage:
%
% [tmpdmat,colLabels,catflag] = std_builddesignmat(design, trialinfo, 1)
%
% Inputs:
%  design      - Design structure as in the STUDY
%  trialinfo   - Structure of trial information. Each field should be a
%                cell array with one element for each trial
% expanding    - Expand the design matrix
% Optional inputs:
%
% Outputs:
%   tmpdmat    - Design matrix
%   colLabels  - Labels for each column of the design matrix
%   catflag    - Binary vector with dimension equal to the number of columns in the design matrix.
%                [0] mean a continuous regressor, [1] means a categotical one.
%
% See also: std_combtrialinfo , std_plodtmat
%
% Authors: Ramon Martinez-Cancino
%
% Copyright (C) 2015  Ramon Martinez-Cancino, UCSD, INC, SCCN
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [tmpdmat,colLabels,catflag] = std_builddesignmat(design, trialinfo, expanding)

if nargin < 3, expanding = 0; end
ntrials = length(trialinfo);
varindx = 1:length(design.variable);

% checking if 'group' var (temporal commit to detect group vars, right now only detecting variable 'group')
groupindx = find(strcmp({design.variable.label},'group'));
if ~isempty(groupindx)
    varindx(groupindx) = [];
end

tmpdmat   = NaN(ntrials,length(varindx));
catflag   = strcmp({design.variable.vartype}, 'categorical');
colLabels = {design.variable.label};

for i = 1 : length(varindx)
     % case for cont variables
    if ~catflag(i)
         [trialindsx, eventvals] = std_gettrialsind(trialinfo,design.variable(varindx(i)).label, '');
         if ~isempty(trialindsx)
             tmpdmat(trialindsx,i) = eventvals;
         end
    else % case for cat variables
        varvaluetmp = design.variable(varindx(i)).value;

        % Expanding cells cells
        c = 1; varlist = {}; varvallength = []; facval_indx = []; dmatval = 1; jcount = 1;
        for ivar = 1:length(varvaluetmp)

            % Length of varvaluetmp(i)
            if iscellstr(varvaluetmp(ivar))
                varvallength(ivar) = length(varvaluetmp(ivar)) ;
            elseif iscell(varvaluetmp(ivar))
                varvallength(ivar) = length(varvaluetmp{ivar}) ;
            else
                varvallength(ivar) = length(varvaluetmp{ivar}) ;
            end
            % Retrieving value and assigning index in design matrix
            if ~iscell(varvaluetmp{ivar})  && varvallength(ivar) == 1
                varlist{c} = varvaluetmp{ivar};
                varindxjoint{ivar} = c;
                facval_indx(c) = dmatval;
                c = c+1;
            else
                tmpcindx = [];
                for ival = 1: varvallength(ivar)
                    if iscellstr(varvaluetmp{ivar}(ival))
                        varlist{c} =  varvaluetmp{ivar}{ival};
                    else
                        varlist{c} =  varvaluetmp{ivar}(ival);
                    end
                    facval_indx(c) = dmatval;
                    tmpcindx       = cat(1,tmpcindx,c);
                    c = c+1;
                end
                varindxjoint{ivar} = tmpcindx;
            end
            dmatval = dmatval+1;
            %---
            for j = 1 :varvallength(ivar)
                %
                if iscellstr(varlist(varindxjoint{ivar}))
                    facval = varlist{varindxjoint{ivar}(j)};
                elseif iscell(varlist(varindxjoint{ivar}))
                    tmpval = varlist(varindxjoint{ivar});
                    if isnumeric(tmpval{j})
                        facval = tmpval{j};
                    else
                        facval = tmpval(j);
                    end
                elseif  isnumeric(varlist(varindxjoint{ivar}))
                    facval = varlist(varindxjoint{ivar}(j));
                elseif ischar(varlist(varindxjoint{ivar}))
                    %                     facval = varlist(varindxjoint{ivar}(j));
                end
                % Find indices of triasl for facval
                [trialindsx, eventvals] = std_gettrialsind(trialinfo,design.variable(varindx(i)).label, facval);

                % Populating the design matrix
                if ~isempty(trialindsx)
                    tmpdmat(trialindsx,i) = facval_indx(jcount);
                    jcount = jcount+1;
                end
            end
            %---
        end
    end
end
% -------------------------------------------------------------------------
% expand categ var
if expanding == 1
    % count number of columns
    nCols = 0;
    for iCol = 1:size(tmpdmat,2)
        if catflag(iCol)
             nCols = nCols+length(design.variable(varindx(iCol)).value);
        else
            nCols = nCols+1;
        end
    end

    tmpdmatExpanded = NaN(size(tmpdmat,1),nCols);
    countCol = 0;
    for iCol = 1:size(tmpdmat,2)
        if catflag(iCol)
            % get unique values for this given categ var
            uniqueVals = unique(tmpdmat(:,iCol));
            uniqueVals(isnan(uniqueVals)) = [];

            % scan unique values
            prevCol = countCol;
            for iUnique = 1:length(design.variable(varindx(iCol)).value)
                countCol    = countCol+1;
                try tmpval = uniqueVals(iUnique); catch, tmpval = NaN; end
                trialSelect = tmpdmat(:,iCol) == tmpval;
                otherCols = prevCol+setdiff(1:length(uniqueVals), iUnique);
                % putting 0 in other columns ensures we process the NaNs
                tmpdmatExpanded(trialSelect, countCol ) = 1;
                tmpdmatExpanded(trialSelect, otherCols) = 0;

                % get the label for that column
                if ischar(design.variable(varindx(iCol)).value{iUnique})
                    colLabels{countCol} = [design.variable(varindx(iCol)).label '-' design.variable(varindx(iCol)).value{iUnique}];
                elseif iscell(design.variable(varindx(iCol)).value{iUnique})
                    % Concat all vals
                    varnametmp = design.variable(varindx(iCol)).value{iUnique}{1};
                    for ivar = 2: length(design.variable(varindx(iCol)).value{iUnique})
                        varnametmp = [varnametmp '&' design.variable(varindx(iCol)).value{iUnique}{ivar}];
                    end
                    colLabels{countCol} = [design.variable(varindx(iCol)).label '-' varnametmp];
                elseif isnumeric(design.variable(varindx(iCol)).value{iUnique})
                    colLabels{countCol} = [design.variable(varindx(iCol)).label '-' int2str(design.variable(varindx(iCol)).value{iUnique})];
                end
            end
        else
            countCol = countCol+1;
            colLabels{countCol} = design.variable(varindx(iCol)).label;
            tmpdmatExpanded(:,countCol) = tmpdmat(:,iCol);
        end
    end

    tmpdmat = tmpdmatExpanded;
end

tmpdmat(:,end+1) = 1;
colLabels{numel(colLabels)+1}  = 'constant';

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
% The mapping from 61 phones set to 48 phonse set for TIMIT. The mapping
% table is as follows:
%
% cl vcl epi --> sil
% el --> l
% en --> n
% zh --> sh
% aa --> ao
% ix --> ih
% ax --> ah
%
% Author: Xiong Xiao, NTU, Singapore
% Date: 1 Jun 2016

function phone_seq_2 = TIMIT_map39to2(phone_seq_39)

found = strcmpi(phone_seq_39, 'sil');
phone_seq_2 = phone_seq_39;
idx = find(found==0);
for j=1:length(idx)
    phone_seq_2{idx(j)} = 'speech';
end
end

author: singaxiong

repo: SignalGraph


text:
function nav=decode_sp3b(nav,headinfo,fid)
global glc gls
peph=gls.peph;
NMAX=10000; nav.peph=repmat(gls.peph,NMAX,1);

if headinfo.type=='P'
    ns=headinfo.ns;
else
    ns=2*headinfo.ns;
end

%traversal by epoch
while ~feof(fid)

    line=fgetl(fid);
    if strcmp(line,'EOF'),continue;end

    if line(1)~='*'
        fprintf('sp3 invalid epoch\n');
        continue;
    end

    time=str2time(line(4:31));
    if strcmp(headinfo.tsys,'UTC'),time=utc2gpst(time);end
    peph.time=time;

    %initialize peph struct
    for i=1:glc.MAXSAT
        for j=1:4
            peph.pos(i,j)=0;
            peph.std(i,j)=0;
            peph.vel(i,j)=0;
            peph.vst(i,j)=0;
        end
        for j=1:3
            peph.cov(i,j)=0;
            peph.vco(i,j)=0;
        end
    end

    for i=1:ns
        line=fgetl(fid);
        if size(line,2)<4||(line(1)~='P'&&line(1)~='V'),continue;end

        if line(2)==' '
            sys=glc.SYS_GPS;
        else
            sys=code2sys(line(2));
        end
        prn=str2double(line(3:4));
        if sys==glc.SYS_QZS,prn=prn+192;end
        sat=satno(sys,prn);
        if ~sat,continue;end

        for j=1:4
            if j<=3
                m=2;k1=1000;base=headinfo.bfact(1);k2=10^-3;k3=0.1;k4=10^-7;
            else
                m=3;k1=10^-6;base=headinfo.bfact(2);k2=10^-12;k3=10^-10;k4=10^-16;
            end

            val=str2double(line(4+14*(j-1)+1:4+14*(j-1)+14));
            if size(line,2)>61
                std=str2double(line(61+3*(j-1)+1:61+3*(j-1)+m));
            else
                std=0;
            end

            if line(1)=='P'
                if val~=0&&abs(val-999999.999999)>=10^-6
                    peph.pos(sat,j)=val*k1; v=1;
                end
                if base>0&&std>0
                    peph.std(sat,j)=base^std*k2;
                end
            elseif v
                if val~=0&&abs(val-999999.999999)>=10^-6
                    peph.vel(sat,j)=val*k3;
                end
                if base>0&&std>0
                    peph.vst(sat,j)=base^std*k4;
                end
            end
        end
    end

    if v
        if nav.np>size(nav.peph,1)
            nav.peph(nav.np+1:nav.np+NMAX,1)=repmat(gls.peph,NMAX,1);
        end
        nav.peph(nav.np+1)=peph;
        nav.np=nav.np+1;
    end

end
fclose(fid);

if nav.np<size(nav.peph,1)
    nav.peph(nav.np+1:end,:)=[];
end

return

repo: GINav


text:
function exact = p51_exact ( )

%*****************************************************************************80
%
%% P51_EXACT returns the exact integral for problem 51.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  exact = - ( r8_ci ( 1.0 ) * sin ( 1.0 ) + ...
    ( 0.5 * pi  - r8_si ( 1.0 ) ) * cos ( 1.0 ) ) / pi;

  return
end


text:
function h = GB_plotboneVRML_New(bone_model,RT,faceColor,lightsOn,parentName,faceAlpha)

bone_conn1 = bone_model.conn;
bone_pts   = bone_model.pts;

if ~exist('lightsOn','var')
    lightsOn = 1;
end

if ~exist('faceColor','var')
    faceColor =  [0.9882 0.9490 0.9020];%[0.2 0.5 0.2];
end
if ~exist('RT','var')
    RT = eye(4,4);
end
if ~exist('faceAlpha','var')
    faceAlpha = 0.4;
end

new_pts = (RT*[bone_pts ones(length(bone_pts),1)]')';

bone_pts = new_pts;

% % % Check Connections
% % if min(min(bone_conn1(1,1:3))) == 1
% %     bone_conn1(:,1:3) = bone_conn1(:,1:3)-1;
% % end

if ~exist('parentName','var')
    h = trisurf(bone_conn1(:,1:3)+1, bone_pts(:,1),bone_pts(:,2),bone_pts(:,3),...
        'EdgeColor','none','LineStyle','none','FaceLighting','phong','FaceColor',faceColor,'FaceAlpha',faceAlpha); hold on
    if lightsOn ~= 0
        GB_giveMeSomeLights(lightsOn);
    end

else
    hold(parentName,'on');
    h = trisurf(bone_conn1(:,1:3)+1, bone_pts(:,1),bone_pts(:,2),bone_pts(:,3),...
        'Parent', parentName,'EdgeColor','none','LineStyle','none','FaceLighting','phong','FaceColor',faceColor,'FaceAlpha',faceAlpha);
    if lightsOn ~= 0
        GB_giveMeSomeLights(lightsOn,parentName);
    end

end

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function tw = trigauss ( n, alpha, beta )

%*****************************************************************************80
%
%% TRIGAUSS computes a trigonometric gaussian quadrature formula.
%
%  Discussion:
%
%    This function computes the N+1 angles and weights of a trigonometric
%    gaussian quadrature formula on [ALPHA,BETA], with
%    0 < BETA - ALPHA <= pi.
%
%    The formula integrates the canonical trigonometric basis with accuracy
%    from about 10^(-15) (for small omega) to about 10^(-13) (for omega --> pi)
%    up to N = 300.
%
%  Modified:
%
%    19 May 2013
%
%  Author:
%
%    Gaspare Da Fies, Alvise Sommariva, Marco Vianello
%
%  Parameters:
%
%    Input, integer N, the trigonometric degree of exactness.
%
%    Input, real ALPHA, BETA, the angular interval.
%    0 < BETA - ALPHA <= pi.
%
%    Output, real TW(N+1,2) array of angles and weights.
%
  np1 = n + 1;
%
%  Compute the half angle subtended by the circle segment.
%
  omega = ( beta - alpha ) / 2.0;
%
%  Compute the modified Chebyshev recursion coefficients.
%
  ab = r_subchebyshev ( np1, omega );
%
%  Compute the quadrature rule.
%
  xw_symm_eigw = SymmMw ( np1, ab );
  tw = trigauss_conversion ( xw_symm_eigw, omega );
%
%  Adjust angles from the reference system to the physical system.
%
  tw(:,1) = tw(:,1) + ( beta + alpha ) / 2.0;

  return
end
function ab = r_subchebyshev ( n, omega )

%*****************************************************************************80
%
%% R_SUBCHEBYSHEV computes the modified Chebyshev recursion coefficients.
%
%  Modified:
%
%    19 May 2013
%
%  Author:
%
%    Gerard Meurant, Alvise Sommariva
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real OMEGA, the arc angle.
%
%    Output, real AB(NN,2), the coefficients of the three term recursion.
%    The dimension NN is N + 1 if N is odd, or N if N is even.
%
  N = n;
  n = n - 1;

  if rem ( N, 2 ) == 1
    NN = N + 1;
    nn = n + 1;
  else
    NN = N;
    nn = n;
  end
%
%  Compute the moments.
%
  mom = fast_moments_computation ( omega, 2 * nn + 1 );
%
%  Recurrence coefficients of the monic Chebyshev polynomials.
%
  abm(:,1) = zeros ( 2 * nn + 1, 1 );
  abm(:,2) = 0.25 * ones ( 2 * nn + 1, 1 );
  abm(1,2) = pi;
  abm(2,2) = 0.5;
%
%  Recurrence coefficients for the monic orthogonal polynomials
%  with respect to the weight function
%    w(x) = 2 * sin ( omega / 2 ) / sqrt ( 1 - sin^2 ( omega / 2 ) * x^2 )
%  by the modified Chebyshev algorithm.
%
% ab = chebyshev ( NN + 1, mom, abm );
  ab = fast_chebyshev ( NN, mom, abm );

  return
end
function ab = fast_chebyshev ( N, mom, abm )

%*****************************************************************************80
%
%% SUBP_MOD_CHEBYSHEV carries out a modified Chebyshev algorithm.
%
%  Discussion:
%
%    This works only for the subperiodic weight function.
%
%    This is a simplified version of a routine by Walter Gautschi.
%
%  Modified:
%
%    19 May 2013
%
%  Parameters:
%
%    Input, integer N, ?
%
%    Input, real MOM(2*N), ?
%
%    Input, real ABM(?), ?
%
%    Output, real AB(N,2), ?
%
  ab = zeros(N,2);
  sig = zeros(N+1,2*N);

  ab(1,2) = mom(1);

  sig(1,1:2*N) = 0;
  sig(2,:) = mom(1:2*N);

  for n = 3:N+1
    for m = n-1:2*N-n+2
      sig(n,m) = sig(n-1,m+1) + abm(m,2) * sig(n-1,m-1) - ab(n-2,2) * sig(n-2,m);
    end

    ab(n-1,2) = sig(n,n-1) / sig(n-1,n-2);
  end

  return
end
function mom = fast_moments_computation ( omega, n )

%*****************************************************************************80
%
%% FAST_MOMENTS_COMPUTATION computes the moments.
%
%  Modified:
%
%    19 May 2013
%
%  Parameters:
%
%    Input, real OMEGA, the arc angle.
%
%    Input, integer N, the index of the highest moment to compute.
%
%    Output, real MOM(1,N+1), the 0-th through N-th moments.
%
  mom = zeros(1,n+1);
%
%  Set the first moment.
%
  mom(1) = 2.0 * omega;

  if ( 2 <= n )

    if ( omega <= 1/4*pi)
      l = 10;
    elseif ( omega <= 1/2*pi)
      l = 20;
    elseif ( omega <= 3/4*pi)
      l = 40;
    elseif ( omega == pi )
      l = 2*ceil(10*pi);
    else
      l = 2*ceil(10*pi/(pi-omega));
    end
%
%  Auxilliary vectors.
%
    temp=(2:2:n+2*l-2);
    temp2=temp.^2-1;
%
%  Diagonals.
%
    dl = 1/4 -1./(4*(temp-1));
    dc = 1/2 -1/sin(omega/2)^2 -1./(2*temp2);
    du = 1/4 +1./(4*(temp+1));

    d = 4*cos(omega/2)/sin(omega/2)./temp2';
    d(end) = d(end);
%
%  Solve the tridiagonal system.
%
    z = tridisolve ( dl(2:end), dc, du(1:end-1 ), d );
%
%  Set the odd moments.
%
    mom(3:2:n+1) = z(1:floor(n/2));

  end

  mom = mom';
%
%  Normalize.
%
  M = length ( mom );
  kk = 2.^(-((1:2:M)-2))';
  kk(1) = 1;
  v = ones(M,1);
  v(1:2:M) = kk;
  mom = v .* mom;

  return
end
function xw = SymmMw ( N, ab )

%*****************************************************************************80
%
%% SYMMMW computes a quadrature rule for a symmetric weight function.
%
%  Discussion:
%
%    This function uses the reduced matrix and eig and
%    computation of weights with the 3-term recurrence.
%
%  Modified:
%
%    19 May 2013
%
%  Author:
%
%    Gerard Meurant, Alvise Sommariva
%
%  Reference:
%
%    Gerard Meurant, Alvise Sommariva,
%    Fast variants of the Golub and Welsch algorithm for symmetric
%    weight functions,
%    Submitted, 2012.
%
%  Parameters:
%
%    Input, integer N, the cardinality of the rule.
%
%    Input, real AB(*,2), the 3-term recurrence for the orthogonal polynomials
%    same as in OPQ.  Note that AB(1,2) is the 0th moment.
%
%    Output, real XW(N,2), the nodes and weights of the quadrature rule.
%
  N0 = size ( ab, 1 );

  if ( N0 < N )
    error('SymmMw: input array ab is too short')
  end

  na = norm(ab(:,1));

  if na > 0
    error('SymmMw: the weight function must be symmetric')
  end
%
%  Computation of the reduced matrix in vectors (a,b)
%
  if mod(N,2) == 0
    even = 1;
    Nc = N / 2;
  else
    even = 0;
    Nc = fix(N / 2) +1;
  end

  absd = ab(:,2);
  absq = sqrt(absd);

  a = zeros(1,Nc);
  b = a;

  switch even
    case 1
        % N even
        a(1) = absd(2);
        b(1) = absq(2) * absq(3);

        k = (2:Nc-1);
        a(k) = absd(2*k-1) + absd(2*k);
        b(k) = absq(2*k) .* absq(2*k+1);
        a(Nc) = absd(N) + absd(N-1);
        start = 1;

       J = diag(a) + diag(b(1:Nc-1),1) + diag(b(1:Nc-1),-1);
       t = sort(eig(J));
       w = weights_3t(t',a,b);
%
% w are the squares of the first components
%
       w = w' / 2;
    case 0
        % N odd
        a(1) = absd(2);
        b(1) = absq(2) * absq(3);

        k = (2:Nc-1);
        a(k) = absd(2*k-1) + absd(2*k);
        b(k) = absq(2*k) .* absq(2*k+1);
        a(Nc) = absd(N);
        start = 2;
%
%  the first node must be zero
%
        J = diag(a) + diag(b(1:Nc-1),1) + diag(b(1:Nc-1),-1);
        t = sort(eig(J));
        t(1) = 0;
        w = weights_3t(t',a,b);
        w = [w(1); w(2:end)' / 2];
    otherwise
        error('this is not possible')
  end

  xwp = sqrt(t);

  xw(:,1) = [-xwp(end:-1:start,1); xwp];
  xw(:,2) = ab(1,2) * ([w(end:-1:start); w]);

  return
end
function tw = trigauss_conversion ( xw, omega )

%*****************************************************************************80
%
%% TRIGAUSS_CONVERSION ???
%
%  Modified:
%
%    19 May 2013
%
%  Author:
%
%    Gaspare Da Fies, Alvise Sommariva, Marco Vianello
%
%  Parameters:
%
%    Input, real XW(?,2), ?
%
%    Input, real OMEGA, the arc angle.
%
%    Output, real TW(?,2), the angles and weights for the trigonometic
%    quadrature rule.
%
  tw(:,1) = 2.0 * asin ( sin ( omega / 2.0 ) * xw(:,1) );
  tw(:,2) = xw(:,2);

  return
end
function w = weights_3t ( t, a, b )

%*****************************************************************************80
%
%% WEIGHTS_3T computes squares of the 1st components of eigenvectors.
%
%  Discussion:
%
%    The results are computed from the 3-term recurrence relation of
%    the orthogonal polynomials.
%
%  Modified:
%
%    19 May 2013
%
%  Author:
%
%    Gerard Meurant, Alvise Sommariva
%
%  Parameters:
%
%    Input, real T(N), the nodes.
%
%    Input, real A(N-1), B(N-1), the coefficients of the 3-term recurrence.
%
%    Output, real W(N), the squares of the first components of the eigenvectors.
%
  n = length ( t );

  P = zeros ( n, n );
  P(1,:) = ones ( 1, n );
  P(2,:) = ( t - a(1) ) / b(1);

  for k = 3 : n
    k1 = k - 1;
    k2 = k - 2;
    P(k,:) = ( ( t - a(k1) ) .* P(k1,:) - b(k2) * P(k2,:) ) / b(k1);
  end

  P2 = P .* P;

  w = 1.0 ./ sum ( P2 );

  return
end


text:
function [G, tt]=gsp_swiss_roll(N,rand_state,param)
%GSP_SWISS_ROLL Initialize a swiss roll graph
%   Usage:  G = gsp_swiss_roll(N,rand_state,param);
%
%   Input parameters:
%         N          : Number of vertices.
%         s          : sigma ( default: sqrt(2/N))
%         thresh     : threshold (default: 1e-6)
%         rand_state : rand seed (default: 0)
%   Output parameters:
%         G     : Graph structure.
%
%   'gsp_create_swiss_roll(N,s,thresh,rand_state)' initializes a graph
%   structure containing the swiss roll graph
%
%   Example:::
%
%          G = gsp_swiss_roll(200);
%          gsp_plot_graph(G);
%

% Author : David I Shuman, Nathanael Perraudin
% Test: test_graphs

if nargin<1
    N = 200;
end

if nargin<2
   rand_state = 0;
end

gsp_reset_seed(rand_state);

a = 1;   % swiss roll goes from a*pi to b*pi
b = 4;
y = rand(2,N);
% uniform distribution along the manifold (in data space)
tt = sqrt((b*b-a*a)*y(1,:)+a*a);
tt = pi*tt;
% now tt should go from a*pi to b*pi
height = y(2,:);
x = [tt.*cos(tt)/b^2; height; tt.*sin(tt)/b^2];

if nargin<3
    param = struct;
end

if ~isfield(param,'k'), param.k = 6; end
G = gsp_nn_graph(x',param);
G.map_coord = y';
G.type = 'Swiss Roll';

end

repo: gspbox


text:
classdef CF8 < PROBLEM
% <multi> <real> <large/none> <constrained>
% Constrained benchmark MOP

%------------------------------- Reference --------------------------------
% Q. Zhang, A. Zhou, S. Zhao, P. N. Suganthan, W. Liu, and S. Tiwari,
% Multiobjective optimization test instances for the CEC 2009 special
% session and competition, School of CS & EE, University of Essex, Working
% Report CES-487, 2009.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 3;
            if isempty(obj.D); obj.D = 10; end
            obj.lower    = [0,0,zeros(1,obj.D-2)-4];
            obj.upper    = [1,1,zeros(1,obj.D-2)+4];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values and constraint violations
        function Population = Evaluation(obj,varargin)
            X  = varargin{1};
            X  = max(min(X,repmat(obj.upper,size(X,1),1)),repmat(obj.lower,size(X,1),1));
            D  = size(X,2);
            J1 = 4 : 3 : D;
            J2 = 5 : 3 : D;
            J3 = 3 : 3 : D;
            Y  = X - 2*repmat(X(:,2),1,D).*sin(2*pi*repmat(X(:,1),1,D)+repmat(1:D,size(X,1),1)*pi/D);
            PopObj(:,1) = cos(0.5*X(:,1)*pi).*cos(0.5*X(:,2)*pi) + 2*mean(Y(:,J1).^2,2);
            PopObj(:,2) = cos(0.5*X(:,1)*pi).*sin(0.5*X(:,2)*pi) + 2*mean(Y(:,J2).^2,2);
            PopObj(:,3) = sin(0.5*X(:,1)*pi)                     + 2*mean(Y(:,J3).^2,2);
            PopCon      = 1 - (PopObj(:,1).^2+PopObj(:,2).^2)./(1-PopObj(:,3).^2) + 4*abs(sin(2*pi*((PopObj(:,1).^2-PopObj(:,2).^2)./(1-PopObj(:,3).^2)+1)));
            Population  = SOLUTION(X,PopObj,PopCon,varargin{2:end});
            obj.FE      = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            N      = ceil(N/5)*5;
            R      = zeros(N,3);
            R(:,3) = repmat(sin((0:1/(N/5-1):1).*pi/2)',5,1);
            for i = 0 : 4
                R(i*N/5+1:(i+1)*N/5,1) = sqrt(i/4*(1-R(i*N/5+1:(i+1)*N/5,3).^2));
            end
            R(:,2) = sqrt(max(1-R(:,1).^2-R(:,3).^2,0));
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [k, innerProd, arg, denom, numer, vec] = polyardKernCompute(kern, x, x2)

% POLYARDKERNCOMPUTE Compute the POLYARD kernel given the parameters and X.
% FORMAT
% DESC computes the kernel parameters for the automatic relevance determination polynomial
% kernel given inputs associated with rows and columns.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : the input matrix associated with the rows of the kernel.
% ARG x2 : the input matrix associated with the columns of the kernel.
% RETURN k : the kernel matrix computed at the given points.
%
% FORMAT
% DESC computes the kernel matrix for the automatic relevance determination polynomial
% kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : the kernel matrix computed at the given points.
%
% SEEALSO : polyardKernParamInit, kernCompute, kernCreate, polyardKernDiagCompute
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006

% KERN

scales = sparse(diag(sqrt(kern.inputScales)));
x = x*scales;

if nargin < 3
  innerProd = x*x';
  arg = innerProd*kern.weightVariance + kern.biasVariance;
  k = kern.variance*arg.^kern.degree;
else
  x2 = x2*scales;
  innerProd = x*x2';
  arg = innerProd*kern.weightVariance + kern.biasVariance;
  k = kern.variance*arg.^kern.degree;
end

author: SheffieldML

repo: GPmat


text:
function [ output_args ] = ESO_3rd( input_args )
%ESO_3RD Summary of this function goes here
%   Detailed explanation goes here

end

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
function K = covMaternard(d, hyp, x, z, i)

% Matern covariance function with nu = d/2 and with Automatic Relevance
% Determination (ARD) distance measure. For d=1 the function is also known as
% the exponential covariance function or the Ornstein-Uhlenbeck covariance
% in 1d. The covariance function is:
%
%   k(x^p,x^q) = sf^2 * f( sqrt(d)*r ) * exp(-sqrt(d)*r)
%
% with f(t)=1 for d=1, f(t)=1+t for d=3 and f(t)=1+t+t²/3 for d=5.
% Here r is the distance sqrt((x^p-x^q)'*inv(P)*(x^p-x^q)), where the P matrix
% is diagonal with ARD parameters ell_1^2,...,ell_D^2, where D is the dimension
% of the input space and sf2 is the signal variance. The hyperparameters are:
%
% hyp = [ log(ell_1)
%         log(ell_2)
%          ..
%         log(ell_D)
%         log(sf) ]
%
% Copyright (c) by Hannes Nickisch, 2013-10-13.
%
% See also COVFUNCTIONS.M.

if nargin<3, K = '(D+1)'; return; end              % report number of parameters
if nargin<4, z = []; end                                   % make sure, z exists
xeqz = isempty(z); dg = strcmp(z,'diag');                       % determine mode

[n,D] = size(x);
ell = exp(hyp(1:D));
sf2 = exp(2*hyp(D+1));
if all(d~=[1,3,5]), error('only 1, 3 and 5 allowed for d'), end         % degree

switch d
  case 1, f = @(t) 1;               df = @(t) 1./t;     % df(t) = (f(t)-f'(t))/t
  case 3, f = @(t) 1 + t;           df = @(t) 1;
  case 5, f = @(t) 1 + t.*(1+t/3);  df = @(t) (1+t)/3;
end
          m = @(t,f) f(t).*exp(-t); dm = @(t,f) df(t).*exp(-t);

% precompute distances
if dg                                                               % vector kxx
  K = zeros(size(x,1),1);
else
  if xeqz                                                 % symmetric matrix Kxx
    K = sq_dist(diag(sqrt(d)./ell)*x');
  else                                                   % cross covariances Kxz
    K = sq_dist(diag(sqrt(d)./ell)*x',diag(sqrt(d)./ell)*z');
  end
end

if nargin<5                                                        % covariances
  K = sf2*m(sqrt(K),f);
else                                                               % derivatives
  if i<=D                                               % length scale parameter
    if dg
      Ki = zeros(size(x,1),1);
    else
      if xeqz
        Ki = sq_dist(sqrt(d)/ell(i)*x(:,i)');
      else
        Ki = sq_dist(sqrt(d)/ell(i)*x(:,i)',sqrt(d)/ell(i)*z(:,i)');
      end
    end
    K = sf2*dm(sqrt(K),f).*Ki;
    K(Ki<1e-12) = 0;                                    % fix limit case for d=1
  elseif i==D+1                                            % magnitude parameter
    K = 2*sf2*m(sqrt(K),f);
  else
    error('Unknown hyperparameter')
  end
end


text:
%MEANC Mean combining classifier
%
%   W = MEANC(V)
%   W = V*MEANC
%
% INPUT
%   V    Set of classifiers (optional)
%
% OUTPUT
%   W    Mean combiner
%
% DESCRIPTION
% If V = [V1,V2,V3, ... ] is a set of classifiers trained on the same
% classes and W is the mean combiner: it selects the class with the mean of
% the outputs of the input classifiers. This might also be used as
% A*[V1,V2,V3]*MEANC in which A is a dataset to be classified.
%
% If it is desired to operate on posterior probabilities then the input
% classifiers should be extended like V = V*CLASSC;
%
% For affine mappings the coefficients may be averaged instead of the
% classifier results by using AVERAGEC.
%
% The base classifiers may be combined in a stacked way (operating in the
% same feature space by V = [V1,V2,V3, ... ] or in a parallel way
% (operating in different feature spaces) by V = [V1;V2;V3; ... ]
%
% EXAMPLES
% PREX_COMBINING
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, VOTEC, MAXC, MINC, MEDIANC, PRODC,
% AVERAGEC, STACKED, PARALLEL

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: meanc.m,v 1.3 2010/06/01 08:48:55 duin Exp $

function w = meanc(p1)

	type = 'mean';               % define the operation processed by FIXEDCC.
	name = 'Mean combiner';      % define the name of the combiner.

	% this is the general procedure for all possible
	% calls of fixed combiners handled by FIXEDCC
	if nargin == 0
		w = prmapping('fixedcc','combiner',{[],type,name});
	else
		w = fixedcc(p1,[],type,name);
	end

	if isa(w,'prmapping')
		w = setname(w,name);
	end

	return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function funhan = piecewise(varargin)
%PIECEWISE Piecewise-defined functions.
%
%  F = PIECEWISE(COND1,DEFN1,...CONDn,DEFNn,DEFAULT) returns a callable
%  function F that applies different definitions according to supplied
%  conditions. For a given X, F(X) will test COND1 and apply DEFN1 if true,
%  etc. If all conditions fail, then DEFAULT is applied. For any particular
%  input, the first condition to match is the only one applied.
%
%  Each condition should be either a:
%    * function handle evaluating to logical values, or
%    * vector [a b] representing membership in the half-open interval [a,b).
%  Each definition should be either a:
%    * function handle, or
%    * scalar value.
%  The DEFAULT definition is optional; if omitted, it will be set to NaN.
%
%  All function definitions can accept multiple input variables, but they
%  all must accept the same number as in the call to F. They also should
%  all be vectorized, returning results the same size and shape as their
%  inputs. Complex inputs will work if the definitions are set up
%  accordingly; in that case, intervals will be tested using only Re parts.
%
%  The special syntax F() displays all the conditions and definitions for F.
%
%  Examples:
%     f = piecewise([-1 1],@(x) 1+cos(pi*x),0);   % a "cosine bell"
%     ezplot(f,[-2 2])
%     g = piecewise(@(x) sin(x)<0.5,@sin,@(x) 1-sin(x));
%     ezplot(g,[-2*pi 2*pi])
%     h = piecewise(@(x,y) (x<0)|(y<0),@(x,y) sin(x-y));  % defined on L
%     ezsurf(h,[-1 1])
%     chi = piecewise(@(x,y,z) x.^2+y.^2+z.^2<1,1,0); % characteristic func
%     [ triplequad(chi,-1,1,-1,1,-1,1), 4/3*pi ]
%         ans =
%                4.1888    4.1888
%
%  See also FUNCTION_HANDLE.

% Copyright (c) 2007 by Toby Driscoll.
% Version 1, 6 Aug 2007.

% If an even number of inputs was given, no default value exists.
if rem(nargin,2)==0
  default = NaN;
else
  default = varargin{end};
end

% Number of condition/definition pairs.
numdefs = floor(nargin/2);
condn = varargin(1:2:2*numdefs);
defn = varargin(2:2:2*numdefs);

funhan = @piecewisefun;   % this is the return value

  % This is the defintion of the returned function.
  function f = piecewisefun(varargin)
  % Special syntax to display the function.
  if nargin==0
    fprintf('\nPiecewise function with conditions/definitions:\n')
    for k = 1:numdefs
      fprintf('  If %s : %s\n',funcornum2str(condn{k}),funcornum2str(defn{k}))
    end
    fprintf('  Otherwise : %s\n\n',funcornum2str(default))
    return
  end

  % Normal call
  f = zeros(size(varargin{1}));
  mask = false(size(f));
  for k = 1:numdefs
    % Determine the values satistfying condition k.
    if isa(condn{k},'function_handle')
      newpts = logical( condn{k}(varargin{:}) );
    else
      newpts = (varargin{1}>=condn{k}(1)) & (varargin{1}<condn{k}(2));
    end
    newpts = newpts & (~mask);
    % Evaluate definition k
    if isa(defn{k},'function_handle')
      % Tricky part: Extract subset of points from each variable.
      sref.type = '()';
      sref.subs = {newpts};
      extract = @(x) subsref(x,sref);
      x = cellfun(extract,varargin,'uniformoutput',false);
      f(newpts) = defn{k}(x{:});
    else
      f(newpts) = defn{k};  % scalar expansion
    end
    mask = mask | newpts;
  end

  % Default case
  if isa(default,'function_handle')
    sref.subs = {~mask};
    extract = @(x) subsref(x,sref);
    x = cellfun(extract,varargin,'uniformoutput',false);
    f(~mask) = default(x{:});
  else
    f(~mask) = default;  % scalar expansion
  end

  end

  % This subfunction converts a function or scalar to its native char form,
  % and a 2-vector into an interval notation.
  function s = funcornum2str(f)
    if isa(f,'function_handle')
      s = char(f);
    elseif numel(f)==1
      s = num2str(f);
    else
      s = [ 'in [' num2str(f(1)) ',' num2str(f(2))  ')' ];
    end
  end

end

repo: mcbench-benchmarks


text:
function [Ns ach_rcu ach_dt ach_gal conv normapx_val] = plot_all(delta, epsil);

if (nargin < 1) || isempty(delta)
	delta = 0.11;
end
if (nargin < 2) || isempty(epsil)
	epsil = 1e-3;
end

Ns = 50:50:2000;
Ns_cap = [Ns(1) Ns(end)];
cap = 1 + delta.*log2(delta) + (1-delta).*log2(1-delta);
Capr = cap + Ns_cap*0;

normapx_val = normapx(Ns, delta, epsil);
conv = converse(Ns, delta,epsil);
ach_gal = gallager_ach(Ns, delta, epsil);
ach_dt = dt_ach(Ns, delta, epsil);

% Compute RCU bound
ach_rcu = [];
for idx=1:length(Ns);
	n=Ns(idx);
	% provide good bracket for speeding up
	plow = max(ach_gal(idx),ach_dt(idx));
	pup = conv(idx);
	ach_rcu(idx) = rcu_ach(n, delta, epsil, plow, pup);
end;

%% normal approximation
figure;
plot(Ns_cap, Capr, 'r--', Ns, conv./Ns, 'r', Ns, normapx_val./Ns, 'k-', Ns, ach_rcu./Ns, 'b', Ns, ach_dt./Ns, 'b--', ...
	Ns, ach_gal./Ns, 'b-.');
xlabel('Blocklen, n'); ylabel('Rate, R'); ylim([0 1.05*cap]);
title(sprintf('Bounds for the BSC(%g), P_{e,max} = %g', delta, epsil));
legend('Capacity', 'Converse',  'Normal approximation', 'RCU achievability', 'DT achievability', 'Gallager achievability', 'Location', 'SouthEast');
grid on

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre