text:
function [seg,glo]=snrseg(s,r,fs,m,tf)
%SNRSEG Measure segmental and global SNR [SEG,GLO]=(S,R,FS,M,TF)
%
%Usage: (1) seg=snrseg(s,r,fs);                  % s & r are noisy and clean signal
%       (2) seg=snrseg(s,r,fs,'wz');             % no VAD or inerpolation used ['Vq' is default]
%       (3) [seg,glo]=snrseg(s,r,fs,'Vq',0.03);  % 30 ms frames
%
% Inputs:    s  test signal
%            r  reference signal
%           fs  sample frequency (Hz)
%            m  mode [default = 'V']
%                 w = No VAD - use whole file
%                 v = use sohn VAD to discard silent portions
%                 V = use P.56-based VAD to discard silent portions [default]
%                 a = A-weight the signals
%                 b = weight signals by BS-468
%                 q = use quadratic interpolation to remove delays +- 1 sample
%                 z = do not do any alignment
%                 p = plot results
%           tf  frame increment [0.01]
%
% Outputs: seg = Segmental SNR in dB
%          glo = Global SNR in dB (typically 7 dB greater than SNR-seg)
%
% This function compares a noisy signal, S, with a clean reference, R, and
% computes the segemntal signal-to-noise ratio (SNR) in dB. The signals,
% which must be of the same length, are split into non-overlapping frames
% of length TF (default 10 ms) and the SNR of each frame in dB is calculated.
% The segmental SNR is the average of these values, i.e.
%         SEG = mean(10*log10(sum(Ri^2)/sum((Si-Ri)^2))
% where the mean is over frames and the sum runs over one particular frame.
% Two optional modifications can be made to this basic formula:
%
%    (a) Frames are excluded if there is no significant energy in the R
%        signal. The idea is to limit the calculation to frames in which
%        speech is active. By default, the voicebox function "activlev" is
%        used to detect the inactive frames (the 'V' mode option).
%
%    (b) In each frame independently, the reference signal is shifted by up
%        to +- 1 sample to find the alignment than minimizes the noise
%        component (S-R)^2. This shifting accounts for small misalignments
%        and/or sample frequency differences between the two signals. For
%        larger shifts, you can use the voicebox function "sigalign".
%        Accurate alignemnt is especially important at high SNR values.
%
% If no M argument is specified, both these modifications will be applied;
% this is equivalent to specifying M='Vq'.

% Bugs/suggestions
% (1) Optionally restrict the bandwidth to the smaller of the two
%     bandwidths either with an extra parameter or automatically determined

%      Copyright (C) Mike Brookes 2011
%      Version: $Id: snrseg.m 2953 2013-05-02 12:51:26Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<4 || ~ischar(m)
    m='Vq';
end
if nargin<5 || ~numel(tf)
    tf=0.01; % default frame length is 10 ms
end
snmax=100;  % clipping limit for SNR

% filter the input signals if required

if any(m=='a')  % A-weighting
    [b,a]=stdspectrum(2,'z',fs);
    s=filter(b,a,s);
    r=filter(b,a,r);
elseif any(m=='b') %  BS-468 weighting
    [b,a]=stdspectrum(8,'z',fs);
    s=filter(b,a,s);
    r=filter(b,a,r);
end

mq=~any(m=='z');
nr=min(length(r), length(s));
kf=round(tf*fs); % length of frame in samples
ifr=kf+mq:kf:nr-mq; % ending sample of each frame
ifl=ifr(end);
nf=numel(ifr);
rf=sum(reshape(r(mq+1:ifl).^2,kf,nf),1);
ef=sum(reshape((s(mq+1:ifl)-r(mq+1:ifl)).^2,kf,nf),1);
if mq
    efm=sum(reshape((s(3:ifl+1)-r(2:ifl)).^2,kf,nf),1);
    efp=sum(reshape((s(1:ifl-1)-r(2:ifl)).^2,kf,nf),1);
    efa=0.5*(efp+efm)-ef;
    efb=0.5*(efp-efm);
    efmk=(abs(efb)<2*efa) & (efa>0); % mask for frames with a valid minimum
    if any(efmk)
        ef(efmk)=ef(efmk)-0.25*efb(efmk).^2./efa(efmk);
    end
    ef=min(min(ef,efm),efp);
end

em=ef==0; % mask for zero noise frames
rm=rf==0; % mask for zero reference frames
snf=10*log10((rf+rm)./(ef+em));
snf(rm)=-snmax;
snf(em)=snmax;

% select the frames to include

if any(m=='w')
    vf=true(1,nf); % include all frames
elseif any(m=='v');
    vs=vadsohn(r,fs,'na');
    nvs=length(vs);
    [vss,vix]=sort([ifr'; vs(:,2)]);
    vjx=zeros(nvs+nf,5);
    vjx(vix,1)=(1:nvs+nf)'; % sorted position
    vjx(1:nf,2)=vjx(1:nf,1)-(1:nf)'; % prev VAD frame end (or 0 or nvs+1 if none)
    vjx(nf+1:end,2)=vjx(nf+1:end,1)-(1:nvs)'; % prev snr frame end (or 0 or nvs+1 if none)
    dvs=[vss(1)-mq; vss(2:end)-vss(1:end-1)];  % number of samples from previous frame boundary
    vjx(:,3)=dvs(vjx(:,1)); % number of samples from previous frame boundary
    vjx(1:nf,4)=vs(min(1+vjx(1:nf,2),nvs),3); % VAD result for samples between prev frame boundary and this one
    vjx(nf+1:end,4)=vs(:,3); % VAD result for samples between prev frame boundary and this one
    vjx(1:nf,5)=1:nf; % SNR frame to accumulte into
    vjx(vjx(nf+1:end,2)>=nf,3)=0;  % zap any VAD frame beyond the last snr fram
    vjx(nf+1:end,5)=min(vjx(nf+1:end,2)+1,nf); % SNR frame to accumulate into
    vf=full(sparse(1,vjx(:,5),vjx(:,3).*vjx(:,4),1,nf))>kf/2; % accumulate into SNR frames and compare with threshold
else  % default is 'V'
    [lev,af,fso,vad]=activlev(r,fs);    % do VAD on reference signal
    vf=sum(reshape(vad(mq+1:ifl),kf,nf),1)>kf/2; % find frames that are mostly active
end
seg=mean(snf(vf));
glo=10*log10(sum(rf(vf))/sum(ef(vf)));

if ~nargout || any (m=='p')
    subplot(311);
    plot((1:length(s))/fs,s);
    ylabel('Signal');
    title(sprintf('SNR = %.1f dB, SNR_{seg} = %.1f dB',glo,seg));
    axh(1)=gca;
    subplot(312);
    plot((1:length(r))/fs,r);
    ylabel('Reference');
    axh(2)=gca;
    subplot(313);
    snv=snf;
    snv(~vf)=NaN;
    snu=snf;
    snu(vf>0)=NaN;
    plot([1 nr]/fs,[glo seg; glo seg],':k',((1:nf)*kf+(1-kf)/2)/fs,snv,'-b',((1:nf)*kf+(1-kf)/2)/fs,snu,'-r');
    ylabel('Frame SNR');
    xlabel('Time (s)');
    axh(3)=gca;
    linkaxes(axh,'x');
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function cnndataout = cnntrain(cnndatain)
%CNNTRAIN Training of convolutional neural network.
%   CNNDATAOUT = CNNTRAIN(CNNDATAIN) Training of convolutional neural
%   network. Type
%
%     >> help cnninfo
%
%   at the prompt for a detailed description of all parameters.
%
%   NOTE: The custoim folder FCNNtoolbox must be in the MATLAB path.
%
%   TRAINING STARTING "FRESH." The inputs to cnntrain are the result of
%   running function cnninit: [cnn,fcnn] = cnninit(cnnparam), whose
%   outputs are the initialized cnn and fcnn. Additional inputs are as
%   listed below.
%
%   cnndatain.cnn = cnn
%   cnndatain.fcnn = fcnnn
%   cnndatain.Images
%   cnndatain.R
%   cnndatain.Epochs
%   cnndatain.MiniBatchSize
%   cnndatain.MSE (set to [] initially)
%   cnndatain.SmoothMSE (set to [] initially)
%
%   TRAINING IN MULTISTAGES.
%   When training takes a long time, it is helpful to train for a
%   manageable number of epochs and then repeat the process, but starting
%   with the last results of training. For this, we make the following
%   modifications in cnndatain before starting the next training session
%   if we want to do another training session without making any changes
%   to any of the cnndatain fields.
%
%   cnndatain = cnndataout;
%
%   If we want to make changes to any of the fields of structure
%   cnndatain, we make following the preceding line of code. For
%   example, to change the number of epochs and continue training, we
%   write:
%
%   cnndatain = cnndataout;
%   cnndatain.Epochs = some new value;
%
%   The other fields of cnndatain are not changed.
%
%   TRAINING STARTING WITH PRETRAINED CNN/FCNN If, instead of starting
%   with a freshly initialized CNN and FCNN, we want to start training
%   with previously-trained CNN and FCNN, we do it as follows. Suppose
%   that the saved the MATLAB workspace after some training session as
%   SavedCNNFCNN_200Epochs.mat. To continue training where the previous
%   training left off, we write:
%
%   load SavedCNNFCNN_200Epochs.mat
%   cnndatain = cnndataout;
%
%   The file SavedCNNFCNN_200Epochs.mat carries with it all the
%   parameters of the last application of cnninit and cnntrain from
%   which they were saved. Any changes in cnndatain, like the number of
%   Epochs, are implemented as explained before.
%
%   The other fields of cnndatain are not changed.
%
%
%   OUTPUTS
%   The outputs of this function in structure cnndataout, with fields:
%
%   cnndataout.cnn
%    The complete, trained cnn.
%
%   cnndataout.fcnn
%    The complete, trained fully-connected net
%
%   cnndataout.MSE
%    A vector containing the training MSE error values as a function of
%    epoch.
%
%   cnndataout.SmoothMSE
%    A smoothed version of cnndataout.MSE
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% SETUP

% Check input image for size M-by-M-by-depth-by-NI, where M is the
% (square) image size, depth is the image depth (e.g., 1 for grayscale
% images and 3 for RGB images), and NI is the number of images.
if (ndims(cnndatain.Images) ~= 4)
   error('Images have to be of sizeM-by-M-by-depth-by-NI')
end

% cnn and fcnn.
cnn = cnndatain.cnn;
fcnn = cnndatain.fcnn;

% Number of cnn layers.
Lc = numel([cnn.NumLayerKernels]);
% Number of fcnn layers.
L = numel([fcnn.NumNodes]);

% Pixel values must be in the range [0,1].
if (max(cnndatain.Images(:)) > 1) || (min(cnndatain.Images(:)) < 0)
   error('Image pixels must be in the range [0,1]')
end

% The input image array of size M-by-M-by-depth-by-NI, where M is the
% (square) image size, depth is the image depth (e.g., 1 for
% grayscale images and 3 for RGB images), and NI is the number of
% images.
NI = size(cnndatain.Images,4);

% Minibatch size. Defaults to the number of patterns.
if ~isfield(cnndatain,'MiniBatchSize')
   cnndatain.MiniBatchSize = size(cnndatain.R,2);
end
MBS = cnndatain.MiniBatchSize;

% The num of images divided by the minibatch size must be an integer.
if rem(NI,MBS) ~= 0
   error('(num images)/(minibatch size) must be an integer')
end
batchReps = NI/MBS;

% Number of training epochs.
NE = cnndatain.Epochs;

% Training errors.
MSE = cnndatain.MSE;
SmoothMSE = cnndatain.SmoothMSE;

% TRAIN THE NETWORK.

h = waitbar(0,'Training the CNN . . . .');
for t = 1:NE
   % Randomize the order of the input images.
   r = randperm(NI);
   for j = 1:batchReps
      % Extract a random imagebatch of MBS images from the input images.
      minibatch = cnndatain.Images(:,:,:,r((j - 1)*MBS + 1 : j*MBS));

      % Extract the corresponding columns from input.R.
      batchR = cnndatain.R(:,r((j - 1)*MBS + 1 : j*MBS));

      % Feedforward through the CNN.
      cnn = cnnff(cnn,minibatch);

      % Feedforward through the FCNN.
      X = cnn(Lc).OutputVectors;
      fcnn = fcnnff(fcnn,X);

      % Backprop through fully-connected net.
      R = batchR;
      fcnn = fcnnbp(fcnn,R);

      % Backprop through cnn.
      delback = fcnn(1).D; % del vector fedback from the fcnn.
      cnn = cnnbp(cnn,delback);

      % Compute average gradients for the minibatch.
      cnn = cnngradients(cnn);

      % Update the cnn weights.
      cnn = cnnupdateweights(cnn);

      % Update the fcnn weights.
      fcnn = fcnnupdateweights(fcnn);

      % Output error.
      fcnnerror = 0.5*(fcnn(L).A - R);

      % Mean-squared error (MSE) for the current minibatch.
      currentMSE = sum(sum((fcnnerror).^2))/size(fcnnerror,2);

      % Update MSE and SmoothMSE.
      if isempty(MSE) % Initialize at beginning.
         MSE(1) = currentMSE;
         SmoothMSE(1) = currentMSE;
      end
      MSE(end + 1) = currentMSE;
      % For smoothing, you could also use MATLAB function smoothdata,
      % but it tends to be much more aggresive that the following
      % "lowpass" filter which is more local and results in a smoothed
      % function whose shape better preserves the nuances of the raw
      % data. If you prefer to use funcdtion smoothdata you can apply ot
      % externally to function cnndataout.MSE
      SmoothMSE(end + 1) = 0.9*SmoothMSE(end) + 0.1*currentMSE;

      waitbar(t/NE)

   end

end
close(h)

% OUTPUT THE TRAINED CNN, FCNN, AND THE MSE ERROR VECTORS.
% Save original inputs.
cnndataout = cnndatain;
% Update the parameters that changed in this function.
cnndataout.cnn = cnn;
cnndataout.fcnn = fcnn;
cnndataout.MSE = MSE;
cnndataout.SmoothMSE = SmoothMSE;

author: dipum

repo: dipum-toolbox

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


text:
function Z = spinv(A, B)
% SPINV    Evaluate the sparse inverse matrix
%
% z = sinv(A)  returns the elements of inv(A)_ij, for which A_ij
%      is different from zero.
%
% z = sinv(LD, 1)  returns the elements of inv(A)_ij, for which A_ij
%     is different from zero, and where LD is the LDL cholesky
%     decomposition of A. LD has to be in the form returned by ldlchol
%     in SuiteSparse by Tim Davis.
%
%   Note! If z = sinv(LD, 1) is used LD must not be modified in Matlab
%   after ldlchol. Matlab destroys the symbolic sparsity structure in
%   the Cholesky decomposition, which is needed in the spinv
%   algorithm. If LD is modified the worst scenario is memory corruption.
%
%
%     For details, see:
%     Jarno Vanhatalo and Aki Vehtari (2008). Modelling local and
%     global phenomena with sparse Gaussian processes. Proceedings of
%     the 24th Conference on Uncertainty in Artificial Intelligence
%

% Copyright (c) 2008-2010      Jarno Vanhatalo

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

    n = size(A,1);

    if nargin == 1
        [LD, p, q] = ldlchol(A);
    else
        LD = A;
    end

    [I,J,ld] = find(LD);
    temp = [I(:) J(:) ; J(:) I(:)];
    temp = sortrows(unique(temp,'rows'),2);
    Iz = temp(:,1); Jz = temp(:,2);

    % Find the column starting points
    a1=zeros(n,1);
    a2 = cumsum(histc(J,1:n));
    a1(1) = 1; a1(2:end) = a2(1:end-1) + 1;
    az1=zeros(n,1);
    az2 = cumsum(histc(Jz,1:n));
    az1(1) = 1; az1(2:end) = az2(1:end-1) + 1;

    for j=1:n
        indaz{j} = az1(j):az2(j);
        indIz{j} = Iz(indaz{j})';
    end

    % Evaluate the sparse inverse
    z = zeros(size(Iz));
    z(end) = 1./ld(end);
    % Allocate memory
    cindit=zeros(n,1);
    for jj = n-1:-1:1
        fil = ld(a1(jj)+1:a1(jj+1)-1);
        fi = I(a1(jj)+1:a1(jj+1)-1);
        lfi = length(fi);
        Zt = zeros(lfi,lfi);
        indz = cumsum(histc(indIz{jj},[0 ; fi]));
        indz = az1(jj) + indz(1:end-1);

        i4=0;
        for i1 = 1:lfi
            cind1=indaz{fi(i1)};
            Icind1=indIz{fi(i1)};
            indfi = lfi;
            i2=length(Icind1);
            go = true;
            while go
                if Icind1(i2)==jj  % Find the indeces for the jj'th rows in fi columns
                    i4=i4+1;
                    cindit(i4)=cind1(i2);
                    go = false;
                end
                if indfi >= 1 && fi(indfi) == Icind1(i2) % Find the indeces for the fi'th rows in i2'nd columns
                    Zt(indfi,i1) = z(cind1(i2));
                    indfi = indfi-1;
                end
                i2 = i2-1;
            end
        end
        % remove extras
        cindi=cindit(1:i4);

        zij = -fil'*Zt;
        z(cindi) = zij;
        z(indz) = zij;
        zij = 1./ld(a1(jj)) - fil'*z(indz);
        z(az1(jj)-1+find(indIz{jj}==jj,1)) = zij;
    end

    Z = sparse(Iz,Jz,z);

    if nargin == 1
        r(q) = 1:n;
        Z = Z(r,r);
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function fileStr = file2string(fname)

fileStr = '';
fid = fopen(fname,'r');
tline = fgetl(fid);
while ischar(tline)
    fileStr = [fileStr ' ' tline];
    tline = fgetl(fid);
end
fclose(fid);

author: jianxiongxiao

repo: ProfXkit


text:
function gT = simwhiteKernGradX(kern, t1, t2)

% SIMWHITEKERNGRADX Gradient of SIM-WHITE kernel with respect to a point t.
% FORMAT
% DESC computes the gradient of the SIM-White (Single Input Motif - White)
% kernel with respect to the input positions.
% ARG kern : kernel structure for which gradients are being computed.
% ARG t1 : locations against which gradients are being computed.
% RETURN gT : the returned gradients. The gradients are returned in
% a matrix which is numData x numInputs x numData. Where numData is
% the number of data points and numInputs is the number of input
% dimensions in t1 (currently always one).
%
% FORMAT
% DESC computes the gradident of the SIM-White (Single Input Motif - White)
% kernel with respect to the input positions where both the row
% positions and column positions are provided separately.
% ARG kern : kernel structure for which gradients are being
% computed.
% ARG t1 : row locations against which gradients are being computed.
% ARG t2 : column locations against which gradients are being computed.
% RETURN gT : the returned gradients. The gradients are returned in
% a matrix which is numData2 x numInputs x numData1. Where numData1 is
% the number of data points in t1, numData2 is the number of data
% points in t2 and numInputs is the number of input dimensions in t1
% (currently always one).
%
% SEEALSO simwhiteKernParamInit, kernGradX, simwhiteKernDiagGradX
%
% COPYRIGHT : David Luengo, 2009

% KERN

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

% Initialisation of the gradient matrix
gT = zeros(size(t1, 1), 1, size(t2, 1));

% Parameters of the kernel required in the computation
variance = kern.variance;
decay = kern.decay;
sensitivity = kern.sensitivity;
isStationary = kern.isStationary;

c = 0.5 * variance * (sensitivity^2);
if (isStationary == false)
    for i = size(t1, 1)
        gT(i, 1, :) =  - sign(t1(i)-t2) .* exp(-decay*abs(t1(i)-t2)) ...
                + exp(-decay*(t1(i)+t2));
    end
else
    for i = size(t1, 1)
        gT(i, 1, :) =  - sign(t1(i)-t2) .* exp(-decay*abs(t1(i)-t2));
    end
end
gT = c * gT;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function quadrule_test07 ( )

%*****************************************************************************80
%
%% TEST07 tests FEJER1_RULE_COMPUTE and FEJER1_RULE_SET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 March 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST07\n' );
  fprintf ( 1, '  FEJER1_RULE_COMPUTE computes\n' );
  fprintf ( 1, '  a Fejer type 1 quadrature rule.\n' );
  fprintf ( 1, '  FEJER1_RULE_SET sets\n' );
  fprintf ( 1, '  a Fejer type 1 quadrature rule.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compare:\n' );
  fprintf ( 1, '    (X1,W1) from FEJER1_RULE_SET\n' );
  fprintf ( 1, '    (X2,W2) from FEJER1_RULE_COMPUTE\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '     Order        W1            W2            X1            X2\n' );
  fprintf ( 1, '\n' );

  for order = 1 : 3 : 10

    [ x1, w1 ] = fejer1_rule_set ( order );
    [ x2, w2 ] = fejer1_rule_compute ( order );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %8d\n', order );

    for i = 1 : order
      fprintf ( 1, '            %12f  %12f  %12f  %12f\n', ...
        w1(i), w2(i), x1(i), x2(i) );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function example4 ( )

%*****************************************************************************80
%
%% EXAMPLE4 uses BVP4C to solve the EXAMPLE4 problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2013
%
%  Author:
%
%    Original MATLAB version by Shampine, Kierzenka, Reichelt.
%    This version by John Burkardt.
%
%  Reference:
%
%    Lawrence Shampine, Jacek Kierzenka, Mark Reichelt,
%    Solving boundary value problems for ordinary differential equations
%    in MATLAB with bvp4c.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EXAMPLE4:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Use BVP4C to solve a boundary value problem with\n' );
  fprintf ( 1, '  a periodic solution of unknown period P.\n' );
  fprintf ( 1, '  We solve this on [0,1], with P an additional unknown.\n' );
  fprintf ( 1, '  y'' = 3 ( y + z - y^3/3 - 1.3\n' );
  fprintf ( 1, '  z'' = - ( y - 0.7 + 0.8 * z ) / 3\n' );
  fprintf ( 1, '  y(0) = y(P), z(0) = z(P)\n' );
  fprintf ( 1, '  Use initial guess y = sin ( 2 pi x ), z = cos ( 2 pi x ).\n' );
%
%  Set SOLINIT, the structure defining the initial guess.
%
  x_init = linspace ( 0.0, 1.0, 5 );
  p_init = 2.0 * pi;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial estimate for P is %g\n', p_init );
  solinit = bvpinit ( x_init, @example4_init, p_init );
%
%  Have BVP4C solve the problem.
%
  sol = bvp4c ( @example4_ode, @example4_bc, solinit );
%
%  Retrieve the updated value of the period.
%
  p = sol.parameters;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Period P = %g\n', p );
%
%  Use DEVAL to evaluate the solution.
%
  x = linspace ( 0.0, 1.0, 101 );
  px = p * x;
  sol = deval ( sol, x );
%
%  Evaluate the initial guess on the initial grid [0,1].
%
  y_init = example4_init ( x );
%
%  Display a plot of Y versus initial guess.
%
  plot ( px, sol(1,:), 'r-', ...
         px, y_init(1,:), 'g-', 'Linewidth', 2 );
  xlabel ( '<--- X --->', 'Fontsize', 16 );
  ylabel ( '<--- Y --->', 'Fontsize', 16 );
  title ( 'EXAMPLE4: Y(red) versus initial guess (green)', 'Fontsize', 16 )
  grid on
  filename = 'example4.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saving plot file as "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EXAMPLE4:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function y = example4_init ( x )

%*****************************************************************************80
%
%% EXAMPLE4_INIT evaluates the initial guess for the solution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the point at which the ODE is to be evaluated.
%
%    Output, real Y(M,1), the value of the initial guess.
%
  y(1,:) = sin ( 2.0 * pi * x );
  y(2,:) = cos ( 2.0 * pi * x );

  return
end
function dydx = example4_ode ( x, y, p )

%*****************************************************************************80
%
%% EXAMPLE4_ODE evaluates the right hand side of the ODE.
%
%  Discussion:
%
%    We assume that the differential equation has been rewritten as a
%    system of first order equations of the form
%
%      dydx = f(x,y)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the point at which the ODE is to be evaluated.
%
%    Input, real Y(M), the value of the solution at X.
%
%    Input, real P, the estimate for the period.
%
%    Output, real DYDX(M), the value of the right hand side given X and Y.
%
  dydx(1,1) = p * 3.0 * ( y(1) + y(2) - ( y(1) )^3 / 3.0 - 1.3 );
  dydx(2,1) =  - p * ( y(1) - 0.7 + 0.8 * y(2) ) / 3.0;

  return
end
function bc = example4_bc ( ya, yb, p )

%*****************************************************************************80
%
%% EXAMPLE4_BC evaluates the boundary conditions.
%
%  Discussion:
%
%    The third boundary condition applies the phase condition, and eliminates
%    the spurious solution y'(x) = 0, and solutions with p = 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real YA(M), YB(M), the solution value at the left and right endpoints.
%
%    Input, real P, the estimate for the period.
%
%    Output, real BC(M), the value of the boundary conditions.
%
  bc(1,1) = ya(1) - yb(1);
  bc(2,1) = ya(2) - yb(2);
  bc(3,1) = p * ( ya(1) - 0.7 + 0.8 * ya(2) ) / 3.0 - 1.0;

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function nc=normcor(a,b)
m=size(a);
n=size(b);
a=im2double(a);
b=im2double(b);
a1=mean2(a);
b1=mean2(b);
c1=0;c2=0;
for i=1:m
    for j=1:n
      c1=c1+(a(i,j)-a1);
      c2=c2+(b(i,j)-b1);
      num=c1*c2;
      c3=(c1^2)*(c2^2);
      dem=sqrt(c3);
    end
end
nc=num/dem;

author: Sable

repo: mcbench-benchmarks


text:
function [J,lambdaOpt,T] = ridgeSVD(Y,Ut, s2,V,nlambda,plotGCV,verb)
%[J,lambdaOpt,T] = ridgeSVD(Y,Ut, s2,V,nlambda,plotGCV)
%
% Estimates a ridge regression model, also know as Tikhonov regularization,
% or minimum norm with L2 prior (or Loreta in the EEG inverse solution literature).
% For an implementation of sLORETA model see the function inverseSolutionLoreta.
%
% Y: measurements (Nsensors X 1)
% Ut, s2,V are defined as the SVD decomposition of the standardized lead field matrix
% nlambda: maximum size of the grid for the hyperparameter lambda, default: 100
% plotGCV: plot the GCV curve (true/false), default: false
% Jest: estimated parapeters
% T: estimated inverse operatormaximum size of the grid for the hyperparameter lambda, default: 100
%
% Jest = argmin(J) ||Y-K*J||^2 + lambda*||L*J||^2 == argmin(J) ||Y-K/L*Jst||^2 + lambda*||I||^2, s.t. J = L/Jst
% and lambda > 0
%
% This code is based on a previous implementation used in Valdes-Hernandez
% et al. (2009), written by Alejandro Ojeda and Pedro Valdez-Hernandez at
% the Cuban Neuroscience Center in 2009.
%
% Author: Alejandro Ojeda, SCCN/INC/UCSD, Jul-2012
%
% References:
%   Pedro A. Valdés-Hernández, Alejandro Ojeda, Eduardo Martínez-Montes, Agustín
%       Lage-Castellanos, Trinidad Virués-Alba, Lourdes Valdés-Urrutia, Pedro A.
%       Valdes-Sosa, 2009. White matter White matter architecture rather than
%       cortical surface area correlates with the EEG alpha rhythm. NeuroImage 49
%       (2010) 2328–2339

if nargin < 4, error('Not enough input arguments.'); end
if nargin < 5 || isempty(nlambda) nlambda = 100; end
if nargin < 6 || isempty(plotGCV) plotGCV = false; end
if nargin < 7 || isempty(verb) verb = false; end

n = size(Ut,1);
p = size(V,1);
s = sqrt(s2);
UtY = Ut*Y;

tol = max([n p])*eps(max(s));
lambda = logspace(log10(tol),log10(max(s)),nlambda);
gcv = zeros(nlambda,1);

beta2 = mean(norms(Y).^2 - norms(UtY).^2);
[n,m] = size(Ut);
delta0 = 0;
  if (m > n && beta2 > 0)
      if verb
        fprintf('m>n criterion met\n');
        fprintf('m=%d, n=%d\n',m,n);
      end
      delta0 = beta2;
  end
for it=1:nlambda
    gcv(it) = gcvfun2(lambda(it),s2,UtY,delta0,m-n);
end

loc = getMinima(gcv);
if verb
    fprintf('GCV min found at loc=%d | lambda=%0.5g\n',loc,lambda(loc)); end
if isempty(loc),
    fprintf('GCV did not find a minimum.\n');
    loc = length(lambda);
end
loc = loc(end);
lambdaOpt = lambda(loc);

T = V*diag(s./(s2+lambdaOpt.^2))*Ut;
J = T*Y;                            % J = (K'*K+lambda*L'*L)\K'*Y

% J = bsxfun(@minus,J,median(J));
% J = bsxfun(@rdivide,J,(std(J)+eps));

if plotGCV
    figure;
    semilogx(lambda,gcv)
    xlabel('log-lambda');
    ylabel('GCV');
    hold on;
    plot(lambdaOpt,gcv(loc),'rx','linewidth',2)
    grid on;
end

%---
function indmin = getMinima(x)
fminor = diff(x)>=0;
fminor = ~fminor(1:end-1) & fminor(2:end);
fminor = [0; fminor; 0];
indmin = find(fminor);

function G = gcvfun2(lambda,s2,beta,delta0,mn,dsvd)

% Auxiliary routine for gcv.  PCH, IMM, Feb. 24, 2008.

% Note: f = 1 - filter-factors.
if (nargin==5)
   f = (lambda^2)./(s2 + lambda^2);
else
   f = lambda./(s2 + lambda);
end
G = mean((norms(bsxfun(@times,f,beta)).^2 + delta0)/(mn + sum(f))^2);

repo: REST


text:
% Template Matlab script to create a BIDS compatible:
%
%   sub-01_ses-01_recording-ManualFullExample_blood.tsv
%   sub-01_ses-01_recording-ManualFullExample_blood.json
%
% This example lists all REQUIRED, RECOMMENDED and OPTIONAL fields.

% Writing json files relies on the JSONio library
% https://github.com/bids-standard/bids-matlab
%
% Make sure it is in the matab/octave path
try
    bids.bids_matlab_version;
catch
    warning('%s\n%s\n%s\n%s', ...
            'Writing the JSON file seems to have failed.', ...
            'Make sure that the following library is in the matlab/octave path:', ...
            'https://github.com/bids-standard/bids-matlab');
end

clear;

this_dir = fileparts(mfilename('fullpath'));
root_dir = fullfile(this_dir, '..', filesep, '..');

project = 'templates';

sub_label = '01';
ses_label = '01';
recording_label = 'ManualFullExample';

name_spec.modality = 'pet';
name_spec.suffix = 'blood';
name_spec.ext = '.tsv';
name_spec.entities = struct('sub', sub_label, ...
                            'ses', ses_label, ...
                            'recording', recording_label);

% using the 'use_schema', true
% ensures that the entities will be in the correct order
bids_file = bids.File(name_spec, 'use_schema', true);

% Contrust the fullpath version of the filename
tsv_file_name = fullfile(root_dir, project, bids_file.bids_path, bids_file.filename);
json_file_name = fullfile(root_dir, project, bids_file.bids_path, bids_file.json_filename);

%% Write TSV
tsv.time = 0;
tsv.plasma_radioactivity = 0;
tsv.metabolite_parent_fraction = 0;
tsv.metabolite_polar_fraction = 0;
tsv.hplc_recovery_fractions = 0;
tsv.whole_blood_radioactivity = 0;

bids.util.tsvwrite(tsv_file_name, tsv);

%% Write JSON
json.PlasmaAvail = '';
json.MetaboliteAvail = '';
json.WholeBloodAvail = '';
json.DispersionCorrected = '';
json.WithdrawalRate = '';
json.TubingType = '';
json.TubingLength = '';
json.DispersionConstant = '';
json.Haematocrit = '';
json.BloodDensity = '';
json.PlasmaFreeFraction = '';
json.PlasmaFreeFractionMethod = '';
json.MetaboliteMethod = '';
json.MetaboliteRecoveryCorrectionApplied = '';

% get the definition of several columns,
% use the bids.Schema class from bids matlab
schema = bids.Schema;

columns = {'time'
           'plasma_radioactivity'
           'metabolite_parent_fraction'
           'metabolite_polar_fraction'
           'hplc_recovery_fractions'
           'whole_blood_radioactivity'};

for i = 1:numel(columns)
    def = schema.get_definition(columns{i});
    json.(columns{i}).Description = def.description;
    if isfield(def, 'unit')
        json.(columns{i}).Unit = def.unit;
    end
end

json.plasma_radioactivity.unit = 'kBq/ml';
json.whole_blood_radioactivity.unit = 'kBq/ml';
json.metabolite_parent_fraction.unit = 'unitless';
json.metabolite_polar_fraction.unit = 'unitless';

%% Write JSON
% Make sure the directory exists
bids.util.mkdir(fileparts(json_file_name));
bids.util.jsonencode(json_file_name, json);

author: bids-standard

repo: bids-starter-kit


text:
% composit.m        September 2, 2013

% sun-synchronous, repeating ground track,
% frozen orbit design - Kozai J2 perturbations

% Orbital Mechanics with MATLAB

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

clear all;

global j2 j3 mu req omega argper fi xldot xqp

% astrodynamic and utility constants

j2 = 0.00108263;
j3 = -0.00000254;
mu = 398600.4415;
req = 6378.14;
omega = 0.000072921151467;

pi2 = 2.0 * pi;
dtr = pi / 180.0;
rtd = 180.0 / pi;

x = zeros(3, 1);

clc; home;

fprintf('\n                    program composit');

fprintf('\n\n < sun-synchronous, repeating ground track, frozen orbits >\n\n');

% request initial guesses

while(1)

   fprintf('\nplease input an initial guess for the semimajor axis (kilometers)\n');

   x(1) = input('? ');

   if (x(1) > 0.0)
      break;
   end

end

while(1)

   fprintf('\nplease input an initial guess for the eccentricity (non-dimensional)');
   fprintf('\n(0 <= eccentricity < 1)\n');

   x(2) = input('? ');

   if (x(2) >= 0.0 && x(2) < 1.0)
      break;
   end

end

while(1)

    fprintf('\nplease input an initial guess for the inclination (degrees)');
    fprintf('\n(90 < inclination <= 180)\n');

    x(3) = input('? ');

    if (x(3) > 90.0 && x(3) <= 180.0)
       break;
    end

end

x(3) = x(3) * dtr;

while(1)

   fprintf('\nplease input the number of integer orbits in the repeat cycle\n');

   norbits = input('? ');

   if (norbits > 1)
      break;
   end

end

while(1)

   fprintf('\nplease input the number of integer days in the repeat cycle\n');

   ndays = input('? ');

   if (ndays > 0)
      break;
   end

end

% compute repetition factor

xqp = norbits / ndays;

% fundamental interval

fi = pi2 / xqp;

% required nodal regression rate

xldot = (360.0 / 365.25) * dtr / 86400.0;

argper = 90.0 * dtr;

% solve system of nonlinear equations

n = 3;

maxiter = 200;

[xf, niter, icheck] = snle ('compfunc', x, n, maxiter);

sma = xf(1);

ecc = xf(2);

inc = rtd * xf(3);

% keplerian period

tkepler = pi2 * sqrt(sma ^ 3 / mu);

tnode = pi2 * (1.0 / xqp) / (omega - xldot);

% print results

clc; home;

fprintf('\n                    program composit');

fprintf('\n\n < sun-synchronous, repeating ground track, frozen orbits >\n\n');

fprintf('mean semimajor axis                %12.4f  kilometers \n\n', sma);

fprintf('mean eccentricity                  %12.10f \n\n', ecc);

fprintf('mean orbital inclination           %12.4f  degrees \n\n', inc);

fprintf('mean argument of perigee           %12.4f  degrees \n\n', argper*rtd);

fprintf('keplerian period                   %12.4f  minutes \n\n', tkepler/60);

fprintf('nodal period                       %12.4f  minutes \n\n', tnode/60);

fprintf('number of orbits in repeat cycle   %12.4f  \n\n', norbits);

fprintf('number of days in repeat cycle     %12.4f  \n\n', ndays);

fprintf('ground trace repetition factor     %12.4f  \n\n', xqp);


text:
% SYNTAX:
% [s,tRange] = hmrR_StimRejection_Nirs(t,s,tIncAuto,tIncMan,tRange)
%
% UI NAME:
% Stim_Exclude
%
% DESCRIPTION:
% Excludes stims that fall within the time points identified as
% motion artifacts from HRF calculation.
%
%
% INPUT:
% t:        the time vector (#time points x 1)
% s:        s matrix (#time points x #conditions) containing 1 for
%           each time point and condition that has a stimulus and
%           zeros otherwise.
% tIncAuto: time points (#time points x 1) identified as motion
%           artifacts by processing stream.
% tIncMan:  time points (#time points x 1) identified as motion
%           artifacts by user.
% tRange:   an array of 2 numbers [t1 t2] specifying how many
%           seconds surrounding motion artifacts, tIncMan and tIncAuto,
%           to consider as excluded data and therefore exclude any stimuli
%           which fall within those buffers.
%           Typically values are t1=-2 and t2 equal to the stimulus
%           duration.
%
% OUTPUT:
% s:        s matrix (#time points x #conditions) containing 1 for
%           each time point and condition that has a stimulus that is
%           included in the HRF calculation, -1 for a stimulus that is
%           excluded automatically in the processing stream, -2
%           for each stimulus excluded by a manually set patch and
%           zeros otherwise.
% tRange:   same tRange array as in the input
%
% USAGE OPTIONS:
% Stim_Exclude: [s,tRange] = hmrR_StimRejection_Nirs(t,s,tIncAuto,tIncMan,tRange)
%
% PARAMETERS:
% tRange: [-5.0, 10.0]
%

function [s,tRange] = hmrR_StimRejection_Nirs(t,s,tIncAuto,tIncMan,tRange)

if isempty(tIncAuto)
    tIncAuto = cell(1,1);
end
if isempty(tIncMan)
    tIncMan = cell(1,1);
end
tIncMan = tIncMan{1};

dt = (t(end)-t(1))/length(t);
tRangeIdx = [floor(tRange(1)/dt):ceil(tRange(2)/dt)];

smax = max(s,[],2);
lstS = find(smax==1);
for iS = 1:length(lstS)
    lst = round(min(max(lstS(iS) + tRangeIdx,1),length(t)));
    if ~isempty(tIncAuto) && min(tIncAuto(lst))==0
        s(lstS(iS),:) = -2*abs(s(lstS(iS),:));
    end
    if ~isempty(tIncMan) && min(tIncMan(lst))==0
        s(lstS(iS),:) = -1*abs(s(lstS(iS),:));
    end
end

author: BUNPC


text:
function a = hilbert_inverse ( n )

%*****************************************************************************80
%
%% HILBERT_INVERSE returns the inverse of the Hilbert matrix.
%
%  Formula:
%
%    A(I,J) =  (-1)**(I+J) * (N+I-1)! * (N+J-1)! /
%           [ (I+J-1) * ((I-1)!*(J-1)!)**2 * (N-I)! * (N-J)! ]
%
%  Example:
%
%    N = 5
%
%       25    -300     1050    -1400     630
%     -300    4800   -18900    26880  -12600
%     1050  -18900    79380  -117600   56700
%    -1400   26880  -117600   179200  -88200
%      630  -12600    56700   -88200   44100
%
%  Properties:
%
%    A is symmetric.
%
%    Because A is symmetric, it is normal, so diagonalizable.
%
%    A is almost impossible to compute accurately by general routines
%    that compute the inverse.
%
%    A is integral.
%
%    The sum of the entries of A is N**2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the inverse Hilbert matrix.
%

%
%  Set the (1,1) entry.
%
  a(1,1) = n * n;
%
%  Define Row 1, Column J by recursion on Row 1 Column J-1
%
  i = 1;
  for j = 2 : n
    a(i,j) = -a(i,j-1) * ( ( n + j - 1 ) * ( i + j - 2 ) * ...
      ( n + 1 - j ) ) /( ( i + j - 1 ) * ( j - 1 ) * ( j - 1 ) );
  end
%
%  Define Row I by recursion on row I-1
%
  for i = 2 : n
    for j = 1 : n

      a(i,j) = -a(i-1,j) * ( ( n + i - 1 ) * ( i + j - 2 ) * ( n + 1 - i ) ) ...
        / ( (i+j-1) * ( i - 1 ) * ( i - 1 ) );

    end
  end

  return
end


text:
function [ a, b ] = p25_lim ( )

%*****************************************************************************80
%
%% P25_LIM returns the integration limits for problem 25.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 1.0;

  return
end


text:
% Mathematics Q2791227
% https://math.stackexchange.com/questions/2301266
% Proximal Operator of Huber Loss Function (For L1 Regularized Huber Loss)
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     20/03/2020
%   *   First release.

%% General Parameters

subStreamNumberDefault = 0; %<! Set to 0 for Random

run('InitScript.m');

figureIdx           = 0; %<! Continue from Question 1
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Parameters

numElements = 25;

outOfSetThr     = 1e-6;
outOfSetCost    = 1e9;

paramLambda = rand(1);
paramDelta = 3;

%% Load / Generate Data

vY      = 2 * randn(numElements, 1);
hObjFun = @(vX) 0.5 * sum((vX - vY) .^ 2) + (paramLambda * HuberLoss(vX, paramDelta));
% From https://math.stackexchange.com/a/1650535
% This is the Proximal Operator for the case paramDelta = 1.
hProxHuberLoss1 = @(vX, paramLambda) vX - ((paramLambda * vX) ./ (max(abs(vX), paramLambda + 1)));
% Supporting any Huber Loss Function by Scaling of the Prox
% HuberLoss(vY, paramDelta) = paramDelta * paramDelta * HuberLoss(vY / paramDelta, 1)
hProxHuberLoss = @(vX, paramDelta, paramLambda) paramDelta * hProxHuberLoss1(vX / paramDelta, paramLambda);

%% Solution by CVX

solverString = 'CVX';

tic();

cvx_begin('quiet')
    % cvx_precision('best');
    variable vX(numElements, 1);
    minimize( 0.5 * sum_square(vX - vY) + (0.5 * paramLambda * sum(huber(vX, paramDelta))) );
cvx_end

toc();

disp([' ']);
disp([solverString, ' Solution Summary']);
disp(['The ', solverString, ' Solver Status - ', cvx_status]);
% disp(['The Optimal Value Is Given By - ', num2str(cvx_optval)]);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(vX))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);

%% Solution by Analytic

solverString = 'Analytic';

tic();

% From https://math.stackexchange.com/a/1650535/33
vX = hProxHuberLoss(vY, paramDelta, paramLambda);

toc();

disp([' ']);
disp([solverString, ' Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(vX))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);

%% Solution by Analytic

solverString = 'Analytic';

tic();

% By Boyd's Book
vX = ProxHuberLossBoyd(vY, paramDelta, paramLambda);

toc();

disp([' ']);
disp([solverString, ' Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(vX))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function [cmps2] = Gal2cmps2(Gal)
% Convert acceleration from galileos to centimeters per second-squared
% Chad A. Greene 2012
cmps2 = Gal;
end


text:
classdef prtRegressRvm < prtRegress
    % prtRegressRvm  Relevance vector machine regression object
    %
    %   REGRESS = prtRegressRvm returns a prtRegressRvm object
    %
    %    REGRESS = prtRegressRVM(PROPERTY1, VALUE1, ...) constructs a
    %    prtRegressRvm object REGRESS with properties as specified by
    %    PROPERTY/VALUE pairs.
    %
    %    A prtRegressRvm object inherits all properties from the prtRegress
    %    class. In addition, it has the following properties:
    %
    %   SetAccess = public:
    %    kernels            - A cell array of prtKernel objects specifying
    %                         the kernels to use
    %    verbosePlot        - Flag indicating whether or not to plot during
    %                         training
    %    verboseText        - Flag indicating whether or not to display
    %                         a message during training
    %
    %   SetAccess = private/protected:
    %    learningConverged  - Flag indicating if the training converged
    %    learningResults    - Struct with information about the convergence
    %    beta               - The weights on each of the kernel elements;
    %                         learned during training
    %    Sigma              - The learned covariance
    %    sparseBeta         - The weights on the retained kernel elements;
    %                         learned durning training
    %    sparseKernels      - The retained kernels
    %
    %   This code is based on:
    %       Michael E Tipping, Sparse bayesian learning and the relevance
    %   vector machine, The Journal of Machine Learning Research, Vol 1.
    %
    %   Also see http://en.wikipedia.org/wiki/Relevance_vector_machine
    %
    %   A prtRegressionRvm object inherits the PLOT method from the
    %   prtRegress object, and the TRAIN, RUN, CROSSVALIDATE and KFOLDS
    %   methods from the prtAction object.
    %
    %   Example:
    %
    %   dataSet = prtDataGenNoisySinc;           % Load a prtDataRegress
    %   dataSet.plot;                    % Display data
    %   reg = prtRegressRvm;             % Create a prtRegressRvm object
    %   reg = reg.train(dataSet);        % Train the prtRegressRvm object
    %   reg.plot();                      % Plot the resulting curve
    %   dataSetOut = reg.run(dataSet);   % Run the regressor on the data
    %   hold on;
    %   plot(dataSet.getX,dataSetOut.getX,'c.') % Plot, overlaying the
    %                                           % fitted points with the
    %                                           % curve and original data
    %   legend('Regression curve','Original Points','Kernel Locations Used',0)
    %
    %
    %   See also prtRegress, prtRegressGP, prtRegressLslr

    properties (SetAccess=private)
        name = 'Relevance Vector Machine'  % Relevance Vector Machine
        nameAbbreviation = 'RVM'           % RVM
    end

    properties
        kernels = prtKernelDc & prtKernelRbfNdimensionScale;

        verbosePlot = false;   % Whether or not to plot during training
        verboseText = false;   % Whether or not to plot during training
    end

    properties (Hidden = true)
        learningMaxIterations = 1000;  % Maximum number of iteratoins
        learningConvergedTolerance = 1e-6;
        learningRelevantTolerance = 1e-3;
    end
    properties (SetAccess = 'protected',GetAccess = 'public')
        learningConverged = [];% Whether or not the training converged

        beta = [];      % Estimated in training
        Sigma = [];     % Estimated in training
        sigma2 = [];    % Estimated in training

        sparseBeta = [];% Estimated in training
        sparseKernels = {};% Estimated in training
    end
    methods
         % Allow for string, value pairs
        function Obj = prtRegressRvm(varargin)
            Obj = prtUtilAssignStringValuePairs(Obj,varargin{:});
        end

        function Obj = set.kernels(Obj,val)
            if ~isa(val,'prtKernel')
                error('prt:prtRegressRvm:kernels','kernels must be a prtKernel');
            end

            Obj.kernels = val;
        end

        function Obj = set.verbosePlot(Obj,val)
            if ~prtUtilIsLogicalScalar(val)
                error('prt:prtRegressRvm:verbosePlot','verbosePlot must be a logical value or a positive integer');
            end
            Obj.verbosePlot = val;
        end

        function Obj = set.verboseText(Obj,val)
            if ~prtUtilIsLogicalScalar(val)
                error('prt:prtRegressRvm:verboseText','verboseText must be a logical value or a positive integer');
            end
            Obj.verboseText = val;
        end
    end

    methods (Access = protected, Hidden = true)

        function Obj = trainAction(Obj,DataSet)
            %Rvm = trainAction(Rvm,DataSet) (Private; see prtClass\train)
            %   Implements Jefferey's prior based training of a relevance
            %   vector machine.  The Rvm output from this function contains
            %   fields "sparseBeta" and "sparseKernels"

            warningState = warning;

            if DataSet.nTargetDimensions ~= 1
                error('prt:prtRegressRvm:tooManyTargets','prtRegressRvm can only operate on single target data.');
            end

            y = DataSet.getTargets(1:DataSet.nObservations,1);

            localKernels = Obj.kernels.train(DataSet);
            gram = localKernels.run_OutputDoubleArray(DataSet);
            nBasis = size(gram,2);

            if Obj.verboseText
                fprintf('RVM training with %d possible vectors.\n', nBasis);
            end

            Obj.beta = zeros(nBasis,1);

            relevantIndices = true(nBasis,1); % Everybody!

            alpha = ones(nBasis,1); % Initialize

            Obj.sigma2 = var(y); % A descent guess to start with

            for iteration = 1:Obj.learningMaxIterations
                % Given currenet relevant stuff find the weight mean and
                % covariance
                if ~any(relevantIndices)
                    break;
                end
                cPhi = gram(:,relevantIndices);
                A = diag(alpha(relevantIndices));

                sigma2Inv = (Obj.sigma2^-1);

                SigmaInv = A + sigma2Inv*(cPhi'*cPhi);
                Obj.Sigma = inv(SigmaInv);
                mu = sigma2Inv*(SigmaInv\(cPhi'*y));

                % Find the current prediction
                yHat = cPhi*mu;

                % Update A
                logAlphaOld = log(alpha(relevantIndices));

                cG = 1 - alpha(relevantIndices).*diag(Obj.Sigma);
                alpha(relevantIndices) = cG./(mu.^2);

                % Update sigma2
                Obj.sigma2 = norm(y-yHat)./(length(yHat) - sum(cG));

                Obj.beta = zeros(nBasis,1);
                Obj.beta(relevantIndices) = mu;

                %check tolerance for basis removal
                TOL = abs(log(alpha(relevantIndices))-logAlphaOld);
                TOL(isnan(TOL)) = 0; % inf-inf = nan

                if Obj.verboseText
                    fprintf('\t Iteration %d: %d RV''s, Convergence tolerance: %g \n',iteration, sum(relevantIndices), max(TOL));
                end

                if all(TOL < Obj.learningConvergedTolerance) && iteration > 1
                    Obj.learningConverged = true;

                    if Obj.verboseText
                        fprintf('Convergence reached. Exiting...\n\n');
                    end

                    break;
                end
                % We didn't break so we can contiue on
                % Select relevant stuff
                newRelevantIndices = alpha < 1./Obj.learningRelevantTolerance;

                if ~mod(iteration,Obj.verbosePlot)
                    if DataSet.nFeatures == 1
                        Obj.verboseIterationPlot(DataSet,relevantIndices);
                    elseif iteration == 1
                        warning('prt:prtRegressRvm','Learning iteration plot can only be produced for training Datasets with 1 feature.');
                    end
                end

                relevantIndices = newRelevantIndices;
            end

            if Obj.verboseText && iteration == Obj.learningMaxIterations
                fprintf('Exiting...Convergence not reached before the maximum allowed iterations was reached.\n\n');
            end

            % Make sparse represenation
            Obj.sparseBeta = Obj.beta(relevantIndices,1);
            Obj.sparseKernels = localKernels.retainKernelDimensions(relevantIndices);

            % Very bad training
            if isempty(Obj.sparseBeta)
                warning('prt:prtClassRvm:NoRelevantFeatures','No relevant features were found during training.');
            end

            % Reset warning
            warning(warningState);

        end

        function DataSetOut = runAction(Obj,DataSet)

            if isempty(Obj.sparseBeta)
                DataSetOut = DataSet.setObservations(zeros(DataSet.nObservations,Obj.dataSetSummary.nTargetDimensions));
                return
            end

            memChunkSize = 1000; % Should this be moved somewhere?
            n = DataSet.nObservations;

            DataSetOut = prtDataSetRegress(zeros(n,1));
            for i = 1:memChunkSize:n;
                cI = i:min(i+memChunkSize,n);
                cDataSet = prtDataSetRegress(DataSet.X(cI,:));
                gram = Obj.sparseKernels.run_OutputDoubleArray(cDataSet);

                DataSetOut = DataSetOut.setObservations(gram*Obj.sparseBeta, cI);
            end
        end
    end

    methods
        function varargout = plot(Obj)

            HandleStructure = plot@prtRegress(Obj);

            holdState = get(gca,'nextPlot');

            % Plot the kernels
            hold on
            %This only plots the x-coordinate... which is weird, but it's
            %all we can do right now b/c kernels don't know about
            %targets...
            Obj.sparseKernels.plot();
            set(gca, 'nextPlot', holdState);

            varargout = {};
            if nargout > 0
                varargout = {HandleStructure};
            end
        end
    end

    methods (Access=protected,Hidden=true)
        function Obj = verboseIterationPlot(Obj,DataSet,relevantIndices)
            DsSummary = DataSet.summarize;

            [linGrid, gridSize] = prtPlotUtilGenerateGrid(DsSummary.lowerBounds, DsSummary.upperBounds, Obj.plotOptions);

            trainedKernel = train(Obj.kernels, DataSet);
            trainedKernel = trainedKernel.retainKernelDimensions(relevantIndices);
            cPhi = trainedKernel.run_OutputDoubleArray(prtDataSetClass(linGrid));

            yHat = reshape(cPhi*Obj.beta(relevantIndices),gridSize);

            colors = Obj.plotOptions.colorsFunction(Obj.dataSetSummary.nTargetDimensions);
            lineWidth = Obj.plotOptions.lineWidth;
            plot(linGrid,yHat,'color',colors(1,:),'lineWidth',lineWidth);
            hold on
            plot(DataSet);
            plot(trainedKernel);
            hold off
            drawnow;
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


%% BVAR tutorial: Bayesian Dynamic Factor model
% Author:   Filippo Ferroni and  Fabio Canova
% Date:     09/14/2021

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Estimation of a  static and  dynamic  factor  model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; close all;  clc;

addpath ../../cmintools/
addpath ../../bvartools/

rng('default');
rng(999);

% generate factors/observables

% Setting the parameters
% AR1 factors
Phi    = [0.75 0.0;0.1 0.8];
% Cov factors innovation
Sigma  = [1 -0.2; -0.2 2];
Q      =  chol(Sigma,'lower');
% factor loadings
Lambda = [1 0;
          0.3 1;
          0.5 -1.5;
          -1.1 -0.5;
          1.5 1.5;
      5*randn(45,size(Phi,1))];
% persistence of idisyncratic errors
rho    = 0;
% standard deviation of idisyncratic errors
sig    = 1;

% preallocating memory
% sample length
T      = 301;
f = zeros(T,size(Phi,1));
y = zeros(T,size(Lambda,1));
e = zeros(T,size(Lambda,1));
i = zeros(T,size(Phi,1));

for t = 1 : T-1
    i(t+1,:) = randn(size(f,2),1)';
    f(t+1,:) = (Phi*f(t,:)' + Q*i(t+1,:)')';
    e(t+1,:) = (rho*e(t,:)' + sig*randn(size(Lambda,1),1))';
    y(t+1,:) = (Lambda*f(t+1,:)' + e(t+1,:)')';
end
f(1,:) = [];
y(1,:) = [];

subplot(2,1,1)
plot(y)
subplot(2,1,2)
plot(f)

% true IRFs
hor = 24;
ir_true=zeros(size(y,2), hor,size(f,2));
true    = iresponse(Phi,eye(size(f,2)),hor,Q);
Lam_    = repmat(Lambda,1,1,size(f,2));
if exist('pagemtimes','builtin') == 5
    ir_true = pagemtimes(Lam_,true);
else
    for ff = 1 : size(f,2)
        ir_true(:,:,ff) = Lambda * true(:,:,ff);
    end
end
pause;

%% principal components

nfac   = 2;
transf = 0;

[~,pc,~,egv,~] = pc_T(y,nfac,transf);

% scree plot
figure,
plot(egv(1:10),'b','Linewidth',2);
grid on
title('Eigenvalues')
% STR_RECAP = [ dirname '/screeplot'];
% saveas(gcf,STR_RECAP,'fig');
% savefigure_pdf([STR_RECAP '.pdf']);

%% static factor model

% 2 static factors
nfac = size(f,2);
lags = 0;
% priors
options.priors.F.Lambda.mean = 0;
options.priors.F.Lambda.cov  = 6;
% estimation command
[BSFM] = bdfm_(y,lags,nfac,options);
% consider a subset of draws
%index = 5000:20:size(BSFM.Phi_draws,3);
index = 1:1:size(BSFM.Phi_draws,3);
% plot estimated and true factors
figure,
for gg=1:nfac
    subplot(nfac,1,gg)
    % draws of factors
    plot(squeeze(BSFM.f_draws(:,gg,index)),'Color',[0.7 0.7 0.7])
    hold on
    % true factor
    f_mean = mean(BSFM.f_draws(:,gg,index),3);
    sign_ = sign(corr(f_mean,f(:,gg)));
    plot(sign_ * f(:,gg),'b','Linewidth',2)
    hold on
    plot(f_mean,'k','Linewidth',2)
end
% save plot
dirname  = 'factor_plt';
mkdir(dirname)
STR_RECAP = [ dirname '/sfm'];
saveas(gcf,STR_RECAP,'fig');
if strcmp(version('-release'),'2022b') == 0
	savefigure_pdf([STR_RECAP '.pdf']);
end
pause;

%% dynamic factor model
clear options
nfac   = size(f,2);
lags   = round(size(Phi,2)/nfac);

% Priors options
% factors priors
options.priors.F.Phi.mean    = Phi;
options.priors.F.Phi.cov     = 10 * eye(size(Phi,1));
options.priors.F.Sigma.scale = Sigma;
options.priors.F.Sigma.df    = 4;
options.priors.F.Lambda.mean = 0;
options.priors.F.Lambda.cov  = 6;
% else Jeffrey prior on PhiSigma: options.priors.name = 'Jeffrey';
% idyosincratic error priors
options.priors.G.Sigma.scale = sig;
options.priors.G.Sigma.df    = 4;
% IRF options
% IV identification
options.proxy    = i(lags+2:end,1);
% Sign identification
options.signs{1} = 'y(4,1:3,1)<0'; %
options.signs{2} = 'y(5,1:3,1)>0'; %
% estimation command
[BDFM] = bdfm_(y,lags,nfac,options);
% consider a subset of draws
% index = 5000:20:size(BDFM.Phi_draws,3);
index = 1:1:size(BDFM.Phi_draws,3);
% plot estimated and true factors
figure,
for gg=1:nfac
    subplot(nfac,1,gg)
    plot(squeeze(BDFM.f_draws(:,gg,index)),'Color',[0.7 0.7 0.7])
    hold on
    f_mean = mean(BDFM.f_draws(:,gg,index),3);
    sign_ = sign(corr(f_mean,f(:,gg)));
    % true factor
    plot(sign_ * f(:,gg),'b','Linewidth',2)
    hold on
    plot(f_mean,'k','Linewidth',2)
end
STR_RECAP = [ dirname '/dfm'];
saveas(gcf,STR_RECAP,'fig');
if strcmp(version('-release'),'2022b') == 0
	savefigure_pdf([STR_RECAP '.pdf']);
end

%%

Phi_m = mean(BDFM.Phi_draws(:,:,index),3);
Sig_m = mean(BDFM.Sigma_draws(:,:,index),3);
Ptt  = lyapunov_symm(Phi_m, Sig_m);
f_m  = mean(BDFM.f_draws(:,:,index),3);
figure,
for gg=1:nfac
    subplot(nfac,1,gg)
    plot([f_m(:,gg)-2*Ptt(gg,gg) f_m(:,gg)+2*Ptt(gg,gg)],'Color',[0.7 0.7 0.7],'Linewidth',2)
    hold on
    plot(f(:,gg),'b','Linewidth',2)
    hold on
    plot(f_m(:,gg),'k','Linewidth',2)
end

Phi0 = Phi';
jj = 0;
figure('Name','Phi')
for gg=1: nfac
    for hh =1 : nfac
        jj = jj+1;
        subplot(nfac,nfac,jj)
        histogram(squeeze(BDFM.Phi_draws(gg,hh,index)),40)
        hold on
        plot([Phi0(gg,hh) Phi0(gg,hh)],[0 50],'r','Linewidth',2)
        xlim([-1 1])
    end
end
%

figure('Name','Sigma')
Sigma0 = Q*Q';
jj = 0;
for gg=1:2
    for hh =1 :2
        jj = jj+1;
        subplot(nfac,nfac,jj)
        histogram(squeeze(BDFM.Sigma_draws(gg,hh,index)),40)
        hold on
        plot([Sigma0(gg,hh) Sigma0(gg,hh)],[0 50],'r','Linewidth',2)
    end
end
%
% figure('Name','Lambda')
% Lambda0 = Lambda;
% jj = 0;
% for gg=1:10
%     for hh =1 :2
%         jj = jj+1;
%         subplot(5,4,jj)
%         histogram(squeeze(BDFM.lambda_draws(gg,hh,index)),40)
%         hold on
%         plot([Lambda0(gg,hh) Lambda0(gg,hh)],[0 30],'r','Linewidth',2)
%     end
% end
%
% figure('Name','sigma_g')
% sig0 = sig*ones(size(Lambda,1),1);
% for jj=1:10
%     subplot(3,4,jj)
%     histogram(squeeze(BDFM.sigma_draws(jj,index)),40)
%     hold on
%     plot([sig0(jj) sig0(jj)],[0 30],'r','Linewidth',2)
% end

%% Plot IRFs

index_var          = [1 2 4 5 15 20];
index_sho          = 1;

% some options:
% add the true IRF
options.add_irfs   = ir_true(index_var,:,index_sho);
% additional 90% HPD set
options.conf_sig_2 = 0.9;
% additional 90% HPD set
options.nplots = [2 3];
% variables names for the plots
options.varnames      = {'Var1','Var2','Var4','Var5','Var15','Var20'};
% name of the directory where the figure is saved
options.saveas_dir    = './factor_plt';
% name of the figure to save
options.saveas_strng  = 'sign';
% sign restricted IRF
irfs_to_plot       = BDFM.irsign_draws(index_var,:,index_sho,:);
plot_irfs_(irfs_to_plot,options)

% IV IRF
irfs_to_plot_iv     = BDFM.irproxy_draws(index_var,:,index_sho,:);
% name of the figure to save
options.saveas_strng  = 'iv';
plot_irfs_(irfs_to_plot_iv,options)


text:
function vel = veldim(vel,from,to)
%VELDIM  Convert velocity units
%
%   velOut = VELDIM(velIn, FROM, TO) converts velIn from the units
%   specified by the string FROM to the units specified by the string
%   TO.  FROM and TO are case-insensitive, and may equal any of the
%   following:
%
%        'meters per second', 'mps' or 'm/s'
%        'feet per second', 'ftps' or 'ft/s'   <== U.S. survey feet
%        'kilometers per hour', 'kmph', 'km/hr', or 'km/h'
%        'knots', 'kts', 'nm/hr' or 'nm/h'
%        'miles per hour', 'mph', 'mi/hr', 'mi/h', sm/hr' or 'sm/h' <== statute miles
%
%   Exercise caution with 'feet' and 'miles'
%   ----------------------------------------
%   VELDIM interprets 'feet per second', 'ftps', and 'ft/s' as U.S. survey
%   feet per second, and does not support international feet at all.
%
%   By definition, one international foot is exactly 0.3048 meters and
%   one U.S. survey foot is exactly 1200/3937 meters.  For many
%   applications, the difference is significant.
%
%   Likewise, VELDIM interprets 'miles per hour', 'mph', 'mi/h', and 'sm/h'
%   as statute miles per hour (also known as U.S. survey miles), and does
%   not support international miles per hour at all.  By definition, one
%   international mile is 5280 international feet and one statute mile is
%   5280 survey feet.
%
%   See also  MPS2FTPS, MPS2KMPH,  MPS2KTS,  MPS2MPH,
%             FTPS2MPS, FTPS2KMPH, FTPS2KTS, FTPS2MPH,
%             KMPH2MPS, KMPH2FTPS, KMPH2KTS, KMPH2MPH,
%             KTS2MPS,  KTS2FTPS,  KTS2KMPH, KTS2MPH,
%             MPH2MPS,  MPH2FTPS,  MPH2KMPH, MPH2KTS

error(nargchk(3, 4, nargin, 'struct'))

% Warn and convert to real if VEL is complex.
vel = ignoreComplex(vel, mfilename, 'VEL');

% Convert units only if there's something to change.
if ~strcmp(from, to)
    vel = applyconversion(vel, from, to);
end

function vel = applyconversion(vel, from, to)

from = lower(from);
to = lower(to);

toIsSupported   = true;
fromIsSupported = true;

switch from
    case {'meters per second','mps','m/s'}
        switch to
            case {'feet per second','ftps','ft/s'}
                vel = mps2ftps(vel);
            case {'kilometers per hour','kmph','km/h','km/hr'}
                vel = mps2kmph(vel);
            case {'knots','kts','nm/h','nm/hr'}
                vel = mps2kts(vel);
            case {'miles per hour','mph','mi/hr','mi/h','sm/hr','sm/h'}
                vel = mps2mph(vel);
            otherwise
                toIsSupported = false;
        end
    case {'feet per second','ftps','ft/s'}
        switch to
            case {'meters per second','mps','m/s'}
                vel = ftps2mps(vel);
            case {'kilometers per hour','kmph','km/h','km/hr'}
                vel = ftps2kmph(vel);
            case {'knots','kts','nm/h','nm/hr'}
                vel = ftps2kts(vel);
            case {'miles per hour','mph','mi/hr','mi/h','sm/hr','sm/h'}
                vel = ftps2mph(vel);
            otherwise
                toIsSupported = false;
        end
    case {'kilometers per hour','kmph','km/h','km/hr'}
        switch to
            case {'feet per second','ftps','ft/s'}
                vel = kmph2ftps(vel);
            case {'meters per second','mps','m/s'}
                vel = kmph2mps(vel);
            case {'knots','kts','nm/h','nm/hr'}
                vel = kmph2kts(vel);
            case {'miles per hour','mph','mi/hr','mi/h','sm/hr','sm/h'}
                vel = kmph2mph(vel);
            otherwise
                toIsSupported = false;
        end
    case {'knots','kts','nm/h','nm/hr'}
        switch to
            case {'feet per second','ftps','ft/s'}
                vel = kts2ftps(vel);
            case {'kilometers per hour','kmph','km/h','km/hr'}
                vel = kts2kmph(vel);
            case {'meters per second','mps','m/s'}
                vel = kts2mps(vel);
            case {'miles per hour','mph','mi/hr','mi/h','sm/hr','sm/h'}
                vel = kts2mph(vel);
            otherwise
                toIsSupported = false;
        end
    case {'miles per hour','mph','mi/hr','mi/h','sm/hr','sm/h'}
        switch to
            case {'feet per second','ftps','ft/s'}
                vel = mph2ftps(vel);
            case {'kilometers per hour','kmph','km/h','km/hr'}
                vel = mph2kmph(vel);
            case {'knots','kts','nm/h','nm/hr'}
                vel = mph2kts(vel);
            case {'meters per second','mps','m/s'}
                vel = mph2mps(vel);
            otherwise
                toIsSupported = false;
        end

    otherwise
        fromIsSupported = false;
end

assert(toIsSupported, 'map:distdim:UnsupportedToUnits', ...
    'Unsupported ''TO'' units: %s.', to)

assert(fromIsSupported, 'map:distdim:UnsupportedFromUnits', ...
    'Unsupported ''FROM'' units: %s.', from)


text:
seed = 1;
rand('state', seed);
randn('state', seed);

N = 4;
dag = zeros(N,N);
C = 1; S = 2; R = 3; W = 4;
dag(C,[R S]) = 1;
dag(R,W) = 1;
dag(S,W)=1;

false = 1; true = 2;
ns = 2*ones(1,N); % binary nodes

bnet = mk_bnet(dag, ns);
if 0
  bnet.CPD{C} = tabular_CPD(bnet, C, [0.5 0.5]);
  bnet.CPD{R} = tabular_CPD(bnet, R, [0.8 0.2 0.2 0.8]);
  bnet.CPD{S} = tabular_CPD(bnet, S, [0.5 0.9 0.5 0.1]);
  bnet.CPD{W} = tabular_CPD(bnet, W, [1 0.1 0.1 0.01 0 0.9 0.9 0.99]);
else
  for i=1:N, bnet.CPD{i} = tabular_CPD(bnet, i); end
end

evidence = cell(1,N);
onodes = [1 3];
data = sample_bnet(bnet);
evidence(onodes) = data(onodes);

clear engine;
engine{1} = belprop_inf_engine(bnet);
engine{2} = jtree_inf_engine(bnet);
engine{3} = global_joint_inf_engine(bnet);
engine{4} = var_elim_inf_engine(bnet);
E = length(engine);

clear mpe;
for e=1:E
  mpe{e} = find_mpe(engine{e}, evidence);
end
for e=2:E
  assert(isequal(mpe{1}, mpe{e}))
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [volOut,axisLabels]=changeAnatomyAxes(inVolume,axisTransform)
% function [volOut,axisLabels]=changeAnatomyAxes(inVolume,axisLabels)
% Rearrange the axes of a medical imaging brain volume.
% INPUT:
% inVolume:          3D input volume.
% axisTransform: This defines the way the axes will be re-arranged.
%              : It is a 3*1 array which must contain a permutation of the digits [1,2,3]
% Basically, we call 'permute' with this array (after doing some flips).

% Example:
% We have a volume y*x*z in size (so size(inVolume)=[y,x,z]
% The axes are y: axial, x: coronal, z: sagittal.
% In other words, if we take an image slice a=inVolume(:,:,fix(z/2));
% then 'a' is a sagittal view of the midline with the nose facing left, and dorsal pointing up.
% Now we want to rearrange the volume so that the same slice (:,:,fix(z/2)) will give us an axial view with the
% nose pointing up and the right ear on the right: y coronal, x sagittal, z axial
% So axisTransform is
% [2 3 1]
% In other words, the first axis in the original image becomes the 3rd axis in the target
% The second axis becomes the first and the third becomes the second.
% To flip the image so that left and right ears were reversed, we would make the second entry -ve: [2 -3 1]

% ARW 03.12.02

% Check number of input argumnents and size of inVolume

if (nargin~=2)
    error('This routine needs 2 input arguments');
end

volSize=size(inVolume);

axisTransform=axisTransform(:);

if (length(axisTransform)~=length(volSize))
    error('axisTransform matrix must be a vector with the same number of entries as there are dimensions in inVolume');
end

% Do flipdim operations
[flippedDims]=find(axisTransform<0);

for t=flippedDims
    flipdim(inVolume,abs(axisTransform(t)));
end

% Do the permutation
volOut=permute(inVolume,abs(axisTransform));

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function realimag(c)
%REALIMAG     Display real and imaginary part of interval Taylor separately
%
%   realimag(c)
%

% written  05/21/09     S.M. Rump
%

  loose = strcmp(get(0,'FormatSpacing'),'loose');

  name = inputname(1);
  if isempty(name)                    % happens for display(taylorinit(random))
    name = 'ans';
  end

  if isreal(c.t)
    display(c,name)
  else
    if loose, disp(' '); end
    display(real(c),['real(' name ')'])
    if loose, disp(' '); end

    if loose, disp(' '); end
    display(imag(c),['imag(' name ')'])
    if loose, disp(' '); end
  end

repo: OpenMAS


text:
function i4_seed_advance_test ( )

%*****************************************************************************80
%
%% I4_SEED_ADVANCE_TEST tests I4_SEED_ADVANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 May 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_SEED_ADVANCE_TEST\n' );
  fprintf ( 1, '  I4_SEED_ADVANCE advances the seed.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Step        SEED input       SEED output\n' );
  fprintf ( 1, '\n' );

  seed_new = 12345;

  for step = 1 : 10

    seed = seed_new;
    seed_new = i4_seed_advance ( seed );

    fprintf ( 1, '  %4d  %16d  %16d\n', step, seed, seed_new );

  end

  return
end


text:
function res = bf_sources_grid(BF, S)
% Generate beamforming grid
% Copyright (C) 2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: bf_sources_grid.m 7703 2019-11-22 12:06:29Z guillaume $

%--------------------------------------------------------------------------
if nargin == 0
    resolution = cfg_entry;
    resolution.tag = 'resolution';
    resolution.name = 'Grid resolution';
    resolution.strtype = 'n';
    resolution.num = [1 1];
    resolution.val = {5};
    resolution.help = {'Select the resolution of the grid (in mm)'};

    space = cfg_menu;
    space.tag = 'space';
    space.name = 'Coordinate system';
    space.help = {'Select the coordinate system in which the grid should be generated'};
    space.labels = {'MNI template', 'MNI-aligned', 'Head', 'Native'};
    space.values = {'MNI template', 'MNI-aligned', 'Head', 'Native'};
    space.val = {'MNI template'};

    constrain = cfg_menu;
    constrain.tag = 'constrain';
    constrain.name = 'Coordinate sources to';
    constrain.help = {'The boundary to which the grid is confined'};
    constrain.labels = {'Inner skull', 'Scalp'};
    constrain.values = {'iskull', 'scalp'};
    constrain.val = {'iskull'};

    grid = cfg_branch;
    grid.tag = 'grid';
    grid.name = 'Grid';
    grid.val = {resolution, space, constrain};

    res = grid;

    return
elseif nargin < 2
    error('Two input arguments are required');
end

if ~isfield(S, 'constrain')
    S.constrain = 'iskull';
end

switch S.constrain
    case 'iskull'
        constraint = export(gifti(BF.data.mesh.tess_iskull), 'ft');
    case 'scalp'
        constraint = export(gifti(BF.data.mesh.tess_scalp), 'ft');
end

M1 = BF.data.transforms.toNative;

switch S.space
    case 'MNI template'
        M1 = BF.data.transforms.toMNI/M1;
        M2 = inv(BF.data.transforms.toMNI);
    case 'MNI-aligned'
        M1 = BF.data.transforms.toMNI_aligned/M1;
        M2 = inv(BF.data.transforms.toMNI_aligned);
    case 'Head'
        M1 = BF.data.transforms.toHead/M1;
        M2 = inv(BF.data.transforms.toHead);
    case 'Native'
        M2 = inv(M1);
        M1 = eye(4);
end

constraint = ft_convert_units(ft_transform_geometry(M1, constraint));

mn = min(constraint.pnt);
mx = max(constraint.pnt);

resolution = S.resolution;

if isequal(constraint.unit, 'm')
    resolution = 1e-3*resolution;
end

% If zero is inside the brain, make sure grid points fall on multiples of
% resolution to ease simulating data from points on the grid
if mn(1)<0 && mx(1)>0
    grid.xgrid = [fliplr(0:-resolution:mn(1)) resolution:resolution:mx(1)];
else
    grid.xgrid = mn(1):resolution:mx(1);
end

if mn(2)<0 && mx(2)>0
    grid.ygrid = [fliplr(0:-resolution:mn(2)) resolution:resolution:mx(2)];
else
    grid.ygrid = mn(2):resolution:mx(2);
end

if mn(3)<0 && mx(3)>0
    grid.zgrid = [fliplr(0:-resolution:mn(3)) resolution:resolution:mx(3)];
else
    grid.zgrid = mn(3):resolution:mx(3);
end

grid.dim   = [length(grid.xgrid) length(grid.ygrid) length(grid.zgrid)];
[X, Y, Z]  = ndgrid(grid.xgrid, grid.ygrid, grid.zgrid);

pos   = [X(:) Y(:) Z(:)];

inside = ft_inside_headmodel(pos, struct('bnd', constraint));

pos    = spm_eeg_inv_transform_points(M2, pos);

grid.allpos  = pos;
grid.inside  = find(inside);
grid.outside = find(~inside);

grid.pos     = pos(inside, :);

res = grid;

author: spm


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Code to run part-based RCNNs for fine-grained detection %
%%%% Usage: put pretrained deep model paths into cnn_models %%
%%%% Also define model definition file %%%%%%%%%%%%%%%%%%%%%%%
%%%% Change all the paths to your path %%%%%%%%%%%%%%%%%%%%%%%
%%%% Written by Ning Zhang %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

init;

results = run_classification(cnn_models, model_def);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [networks, rxnNumGenes] = GPR2models(metabolic_model, selected_rxns, separate_transcript, numWorkers, printLevel)
% Each GPR rule is converted into a network where the reaction and genes
% involved are interconnected.
%
% USAGE:
%
%    [networks, rxnNumGenes] = GPR2models(metabolic_model, selected_rxns, separate_transcript, numWorkers, printLevel)
%
% INPUTS:
%    metabolic_model:        Metabolic model structure (COBRA Toolbox format)
%    selected_rxns:          Index array which indicates selected reactions to
%                            calculate the network.
%    separate_transcript:    Character used to discriminate
%                            different isoforms of a gene. Default ''.
%
% OPTIONAL INPUTS:
%    numWorkers:        Maximum number of workers
%                       * 0 - maximum provided by the system (automatic)
%                       (default). If parallel pool is active, numWorkers
%                       is defined by the system, otherwise is 1.
%                       * 1 - sequential
%                       * 2+ - parallel
%    printLevel:        show the reactions created in models.
%                       * 0 - shows nothing
%                       * 1 - shows progress by reactions (default)
%                       * 2+ - shows everything (reaction and network generation)
%
% OUTPUT:
%    networks:          Cell array which contains all the networks.
%    rxnNumGenes:       Array with the number of genes present in each GPR
%                       rule for each reaction.
%
% .. Authors:
%       - Inigo Apaolaza, Aug 2017, University of Navarra, TECNUN School of Engineering.
%       - Luis V. Valcarcel, Aug 2017, University of Navarra, TECNUN School of Engineering.
%       - Francisco J. Planes, Aug 2017, University of Navarra, TECNUN School of Engineering.
%       - Inigo Apaolaza, April 2018, University of Navarra, TECNUN School of Engineering.

p = inputParser;
% check required arguments
addRequired(p, 'metabolic_model');
addOptional(p, 'selected_rxns', [], @isnumeric);
addOptional(p, 'separate_transcript', '', @(x)ischar(x)); % default is empty
addOptional(p, 'numWorkers', 0, @(x)isnumeric(x)&&isscalar(x)); % Default is gpc('nocreate')
addOptional(p, 'printLevel', 1, @(x)isnumeric(x)&&isscalar(x));
% extract variables from parser
parse(p, metabolic_model, selected_rxns, separate_transcript, numWorkers, printLevel);
metabolic_model = p.Results.metabolic_model;
selected_rxns = p.Results.selected_rxns;
separate_transcript = p.Results.separate_transcript;
numWorkers = p.Results.numWorkers;
printLevel = p.Results.printLevel;

% fill with all reactions
if (isempty(selected_rxns))
    selected_rxns = 1:length(metabolic_model.rxns);
end

% Step 1: check the fields of the input metabolic model
if ~isfield(metabolic_model,'genes')
    error('The genes of the model are not defined');
end
if ~isfield(metabolic_model,'rxns')
    error('The reactions of the model are not defined');
end
if ~isfield(metabolic_model,'grRules') || ~isfield(metabolic_model,'rules') ||~isfield(metabolic_model,'rxnGeneMat')
    error('The GPR rules of the model are not defined');
end

%rxnGeneMat is a required field for this function, so if it does not exist,
%build it.
if ~isfield(metabolic_model,'rxnGeneMat')
    metabolic_model = buildRxnGeneMat(metabolic_model);
end

% Generate a cell array to store the resulting models
networks = cell(size(selected_rxns));

% Examine if there is an existing pool
p = gcp('nocreate'); % If no pool, do not create new one.
if isempty(p)
    poolsize = 0;               % Single core
elseif numWorkers == 0
    poolsize = p.NumWorkers;    % Multi core, all cores in the PC
elseif numWorkers == 1
    poolsize = 0;               % Single core
else
    p = parpool(numWorkers);
    poolsize = p.NumWorkers;    % Multi core, limited by user
end

% Step 2: Create models
if printLevel > 0 && length(selected_rxns)>1
    disp('Calculating Networks for GPR rules...');
end
parfor (i=1:length(selected_rxns),poolsize)
    if printLevel > 0
%         clc
        disp([num2str(i),' of ', num2str(length(selected_rxns)) ,' rxns']);
    end
    RXN = metabolic_model.rxns{selected_rxns(i)};
    % Create empty model
    model = createModel();
    % include metabolic reaction as objective metabolite
    model.mets{1} = RXN;
    % Add objective function
    model.rxns{1} = RXN;
    model.S(1,1) = -1;

    % Check if there is any gene related to this reaction and add them
    genes = metabolic_model.genes(metabolic_model.rxnGeneMat(selected_rxns(i),:)>0);
    if length(genes)>1
        model.mets = vertcat(RXN, genes); % Include genes in the model
%         model = addMetabolite(model,genes);
        model.genes = genes;
        fp = FormulaParser();
        head = fp.parseFormula(metabolic_model.rules{selected_rxns(i)});
        model = modelParser(model,metabolic_model.genes,head,[RXN,'_N_1'],RXN, printLevel);
        model = reduceModel(model);
    elseif length(genes)==1
        model.mets = vertcat(RXN, genes); % Include genes in the model
        model.genes = genes;
        model.rxns(2) = genes;
        model.S(1:2,2) = [+1 -1];
    else
        model.rxns{2} = 'NO_GENES';
        model.S(1,1:2) = [-1 +1];
    end

    % Include Exchange reactions
    [~, gene_idx] = ismember(genes, model.mets);
    demand_idx = (1:length(genes)) + length(model.rxns);
    % compulsitory fields
    model.rxns(demand_idx) = strcat('DM_',genes);
    model.rxns = model.rxns';
    model.S(gene_idx,demand_idx) = eye(length(genes));

    % converge isoforms (if neccesary)
    if ~isempty(separate_transcript) && length(genes)>1
        model = convergeTranscripts2Gene(model, separate_transcript);
    end

    % Generate rest of fields
    idx_rxns = (1:length(model.rxns))';
    idx_mets = (1:length(model.mets))';
    model.b(idx_mets) = 0;
    model.csense(idx_mets) = 'E';
    model.c(idx_rxns) = 0;
    model.c(1) = 1;
    model.lb(idx_rxns) = 0;
    model.ub(idx_rxns) = inf;
    model.ub(1) = 1000;
    model.rules(demand_idx) = {''};
    model.b = model.b';
    model.c = model.c';
    model.lb = model.lb';
    model.ub = model.ub';
    model.rules = model.rules';

    % Store the model
    networks{i} = model;
end

% Calculate the number of genes per GPR rule
rxnNumGenes = zeros(size(selected_rxns));
for i = 1:length(selected_rxns)
    rxn = selected_rxns(i);
    idx_genes = find(metabolic_model.rxnGeneMat(rxn,:));
    rxn_isoforms = metabolic_model.genes(idx_genes);
    rxn_genes = unique(strtok(rxn_isoforms,separate_transcript));
    rxnNumGenes(i) = length(rxn_genes);
end

end

function model = modelParser(model,genes,head,node_parent,parent,printLevel)
% This function look into the GPRparser and builds auxiliary nodes for each
% layer of the GPR rules.
% For example:
%     g1 & g2 & g3 is one layer in the model,   RXN = g1 + g2 + g3
%     (g1 | g2) & g3 is a two layer model:      RXN = (g1 | g2) + g3
%                                               (g1 | g2) = g1
%                                               (g1 | g2) = g2

%  Name nodes according to position and layer inside the GPR rule
node = strcat({[node_parent,'_']},cellfun(@num2str,num2cell(1:length(head.children)),'UniformOutput', false));

% model = addMetabolite(model, node);
model.mets = vertcat(model.mets,node');
[~, idx_new_mets] = ismember(node, model.mets);

idx_parent_mets = find(strcmp(parent, model.mets));

%  Expand model with every node, using OR/AND rules
if isa(head,'OrNode')
    idx_new_rxn = length(model.rxns) + (1:length(head.children));
    model.rxns(idx_new_rxn) = node;
    model.S(idx_new_mets,idx_new_rxn) = -eye(length(head.children));
    model.S(idx_parent_mets,idx_new_rxn) = +1;
else
    idx_new_rxn = length(model.rxns)+1;
    model.rxns{idx_new_rxn} = ['AND_' parent];
    model.S(idx_new_mets,idx_new_rxn) = -1;
    model.S(idx_parent_mets,idx_new_rxn) = +1;
end

for i=1:length(head.children)
    child = head.children(i);
    % check if the child is final layer
    if isa(child,'LiteralNode')
        gen = genes(str2double(child.id));
        idx_new_rxn = length(model.rxns)+1;
        idx_met_gen = strcmp(gen,model.mets);
        idx_met_node = strcmp(node(i),model.mets);
        model.rxns{idx_new_rxn} = [gen{1},'_',node{i}];
        model.S(idx_met_node,idx_new_rxn) = +1;
        model.S(idx_met_gen,idx_new_rxn) = -1;
    else
        model = modelParser(model,genes,child,node{i},node{i},printLevel);
    end
end

end

function modelOut = reduceModel(model)
% As the algorithm generates redundant nodes, this function is used to
% remove these nodes:
% Ej:       gen -> node -> rule  == gen -> rule

% This is done for every gene of the model
genes = model.genes;
for gen=1:length(genes)
    gen_idx = strcmp(model.mets,genes{gen}); % find idx of gene
    rxn_gen_out = find(model.S(gen_idx,:)<0);   % find reactions that consume gene
    % for each reaction, find consumption node
    [node_idx,~] = find(model.S(:,rxn_gen_out)>0);
    for i = 1:length(node_idx)
        % for each reaction, find consumption node
        node_name = model.mets(node_idx(i));
        % find reactions out of this node
        [~, rxn_idx] = find((model.S(node_idx(i),:)<0));
        rxn_name = model.rxns(rxn_idx);
        % in this reaction, change (node -> rule) for (gene -> rule)
        model = changeRxnMets(model, node_name, genes(gen), rxn_name);
    end
    % remove stoichimetri of intermediate nodes
    model.S(node_idx,:)=0;
    % remove output of genes to intermediate nodes
    model.S(:,rxn_gen_out)=0;
end

% remove unused reactions and unused metabolites
idx_rxns = sum(model.S~=0,1)>0;
idx_mets = sum(model.S~=0,2)>0;
modelOut.S = model.S(idx_mets,idx_rxns);
modelOut.rxns = model.rxns(idx_rxns);
modelOut.mets = model.mets(idx_mets);

end

function model = convergeTranscripts2Gene(model, separate_transcripts)
% Some COBRA models have transcripts instead of genes. Most of the
% transcripts have the same functionallity, so we can converge them as the
% same gene.
% Ej:       gene 10005.1
%           gene 10005.2        ==>     gene 10005
%           gene 10005.3
%
% This is done for every gene of the model at exchange reaction level

% seach demand reactions
idx_DM_rxns = startsWith(model.rxns,'DM_');
DM_rxns = model.rxns(idx_DM_rxns);
% obtain isoforms and genes
transcripts = regexprep(DM_rxns,'DM_','');
genes = strtok(transcripts,separate_transcripts);

% Remove all the previously existing exchange reactions
model.S = model.S(:,~idx_DM_rxns);
model.rxns = model.rxns(~idx_DM_rxns);

% Divide into genes with one transcript and two or more transcripts
x = tabulate(genes);
genes_unique = x(:,1);
genes_one_transcrits = x(cell2mat(x(:,2))==1,1);
genes_more_transcripts = x(cell2mat(x(:,2))>1,1);

% For genes with one transcript, just change the name
[~,idx_genes_one_transcript] = ismember(genes_one_transcrits, strtok( model.mets,separate_transcripts));
model.mets(idx_genes_one_transcript) = genes_one_transcrits;

% For genes with more than one transcript, converge them into one
% metabolite
% Generate new metabolites
model = addMetabolite(model, genes_more_transcripts);
% Converge genes
for i = 1:length(genes_more_transcripts)
    gen_idx = find(strcmp(genes_more_transcripts{i}, model.mets));
    transcripts_idx = find(strcmp(genes_more_transcripts{i}, strtok( model.mets,separate_transcripts)));
    transcripts_idx = setdiff(transcripts_idx,gen_idx);
    % new reactions to converge genes
    idx_new_rxns = length(model.rxns)+1 : length(model.rxns) + length(transcripts_idx);
    model.S(transcripts_idx,idx_new_rxns) = +eye(length(transcripts_idx)); % produce transcript
    model.S(gen_idx,idx_new_rxns) = -1; % consume gene
    model.rxns(idx_new_rxns) = strcat('Con_',genes_more_transcripts{i},'_',model.mets(transcripts_idx));
end

% Generate new demand reactions for all genes
[~, gene_idx] = ismember(genes_unique, model.mets);
demand_idx = (1:length(genes_unique)) + length(model.rxns);
model.rxns(demand_idx) = strcat('DM_',genes_unique);
model.S(gene_idx,demand_idx) = eye(length(genes_unique));

end

author: opencobra

repo: cobratoolbox


text:
%% Phase Calculation
% Generates the phase signal for Kalman Filtering
% Inputs
% peaks             QRS peak locations
% length            Length of data for phase generation
%
%
% Fetal Extraction Toolbox, version 1.0, February 2014
% Released under the GNU General Public License
%
% Copyright (C) 2014  Fernando Andreotti
% Dresden University of Technology, Institute of Biomedical Engineering
% fernando.andreotti@mailbox.tu-dresden.de
%
% Last updated : 24-07-2014
%
% Based on: Synthetic ECG model error
% Open Source ECG Toolbox, version 2.0, March 2008
% Released under the GNU General Public License
% Copyright (C) 2008  Reza Sameni
% Sharif University of Technology, Tehran, Iran -- LIS-INPG, Grenoble, France
% reza.sameni@gmail.com

% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version.
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details.

% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version.
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details.
%
function phase = FECGx_kf_PhaseCalc(peaks,lengthx)
% Based on Sameni's method for phase calculation
phase = zeros(1,lengthx);
m = diff(peaks);            % gets distance between peaks
% first interval
% dealing with borders (first and last peaks may not be full waves)
% uses second interval as reference
L = peaks(1);   %length of first interval
if isempty(m) % only ONE peak was detected
    phase(1:lengthx) = linspace(-2*pi,2*pi,lengthx);
else
    phase(1:L) = linspace(2*pi-L*2*pi/m(1),2*pi,L);
    % beats in the middle
    for i = 1:length(peaks)-1;  % generate phases between 0 and 2pi for almos all peaks
        phase(peaks(i):peaks(i+1)) = linspace(0,2*pi,m(i)+1);
    end                                         % 2pi is overlapped by 0 on every loop
    % last interval
    % uses second last interval as reference
    L = length(phase)-peaks(end);   %length of last interval
    phase(peaks(end):end) = linspace(0,L*2*pi/m(end),L+1);
end
phase = mod(phase,2*pi);
phase(find(phase>pi)) = phase(find(phase>pi))- 2*pi;

author: fernandoandreotti


text:
% Function to calculate Threshold for BayesShrink

function threshold=bayes(X,sigmahat)

len=length(X);
sigmay2=sum(X.^2)/len;
sigmax=sqrt(max(sigmay2-sigmahat^2,0));
if sigmax==0 threshold=max(abs(X));
else threshold=sigmahat^2/sigmax;
end


text:
function fibonacci_spiral_connected ( n )

%*****************************************************************************80
%
%% FIBONACCI_SPIRAL_CONNECTED draws connected points on a Fibonacci spiral.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points to plot.
%    Default is 101.
%
  if ( nargin < 1 )
    n = 101;
  end
%
%  PHI is the golden ratio, the limit of the ratio of
%  successive Fibonacci numbers:
%
%    PHI = limit ( N->oo ) F(N+1)/F(N)
%
  phi = ( 1.0 + sqrt ( 5.0 ) ) / 2.0;
%
%  Allocate storage for the point data.
%
  x = zeros ( n, 1 );
  y = zeros ( n, 1 );
%
%  Set the angle and radius of the first point.
%
  a = 0.0;
  r = 0.0;
%
%  Set the increments.
%
  da = 2.0 * pi * ( phi - 1.0 ) / phi;
  dr = 1.0;
%
%  Create a spiral in which the radius R and angle A both
%  increase by a constant increment,
%
  for i = 1 : n
    x(i) = r * cos ( a );
    y(i) = r * sin ( a );
    a = mod ( a + da, 2 * pi );
    r = r + dr;
  end
%
%  SCALE controls how many steps we take between the actual points.
%  A value of 5 is enough to see the basic spiral that connects the points.
%  A vale of 10 would make a smoother spiral.
%
  scale = 5.0;
%
%  Allocate storage for the intermediate data.
%
  n2 = scale * ( n - 1 ) + 1;
  x2 = zeros ( n2, 1 );
  y2 = zeros ( n2, 1 );
%
%  Set the angle and radius of the first point.
%
  a = 0.0;
  r = 0.0;
%
%  Set the increments.
%
  da = 2.0 * pi * ( phi - 1.0 ) / phi;
  dr = 1.0;

  da = da / scale;
  dr = dr / scale;
%
%  Create a spiral in which the radius R and angle A both
%  increase by a constant increment,
%
  for i = 1 : n2
    x2(i) = r * cos ( a );
    y2(i) = r * sin ( a );
    a = mod ( a + da, 2 * pi );
    r = r + dr;
  end
%
%  Display the points,
%  and use the intermediate points to draw lines that display the spiral.
%
  figure
  scatter ( x, y, 'b.' );
  hold on
  plot ( x2, y2, 'r-' )
  axis equal
  title ( sprintf ( 'Connected Fibonacci spiral, N = %d', n ) )
  hold off

  return
end


text:
function C = mtimes(A, B)
%*   Multiplication of VALSDISCRETIZATION discretizations.

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

C = A*B;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function value = r4_li ( x )

%*****************************************************************************80
%
%% R4_LI evaluates the logarithmic integral for an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the logarithmic integral evaluated at X.
%
  persistent sqeps

  if ( isempty ( sqeps ) )
    sqeps = sqrt ( r4_mach ( 3 ) );
  end

  if ( x < 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_LI - Fatal error!\n' );
    fprintf ( 1, '  Function undefined for X <= 0.\n' );
    error ( 'R4_LI - Fatal error!' )
  end

  if ( x == 0.0 )
    value = 0.0;
    return
  end

  if ( x == 1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_LI - Fatal error!\n' );
    fprintf ( 1, '  Function undefined for X = 1.\n' );
    error ( 'R4_LI - Fatal error!' )
  end

  if ( abs ( 1.0 - x ) < sqeps )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_LI - Warning!\n' );
    fprintf ( 1, '  Answer less than half precision.\n' );
    fprintf ( 1, '  X is too close to 1.\n' );
  end

  value = r4_ei ( log ( x ) );

  return
end


text:
%this shows an example of using the "get axes" information from graph_picker

im=imread('example_graph.jpg');%load image
%display('identify reference points in graph picker and ''get axes''');
%[tmp]=graph_picker(im);

figure
image(im);%this is the same image that was used for graph_picker
axis off
h_imax=axes;%get handle to current axes
X=[2.722513e-01 8.272251e-01 1.204188e+00 1.612565e+00 2.146597e+00 2.774869e+00 3.591623e+00 4.219895e+00 ];%sample data
Y=[4.821429e+00 8.928571e+00 1.696429e+01 2.767857e+01 3.142857e+01 2.946429e+01 3.232143e+01 3.732143e+01 ];
plot(X,Y,'ob');
axis([-2.129e+00 6.055e+00 -1.191e+01 5.744e+01]);%this comes from graph_picker
set(h_imax,'color','none');%turn off background color
grid on


text:
% concatenate cell array of 3D tensors
% assume that the cells are MxNxP_i tensors where M and N are shared by all
% cells, while P_i can be different.
% Generate a 3D tensor of MxNxsum(P_i).
%
function output = cell2mat_tensor3D(data, padnumber)
if nargin<2
    padnumber = 0;
end

precision = 'single';
if ~isempty(data) && ~isempty(data{1})
    precision = class(data{1}(1));
end

nCell = length(data);
for i=1:nCell
    [M(i),N(i),P(i)] = size(data{i});
    if M(i)*N(i)==0
        P(i)=0;
    end
end

if sum(abs(N-N(1)))==0  % if all cells have the same number of columns
    N = N(1);
else        % if some cells have less columns than others, zero pad them
    Nmax = max(N);
    for j=1:nCell
        data{j}(:,N(j)+1:Nmax,:) = padnumber;
    end
    N = max(N);
end

if sum(abs(P-1))==0     % all are matrices, not tensors
    output = cell2mat(data);
    output = reshape(output, max(M),N,nCell);
else
    output = zeros(max(M),N, sum(P), precision);
    for j=1:nCell
        idx1 = sum(P(1:j-1))+1;
        idx2 = sum(P(1:j));
        if idx1>idx2; continue; end
        [d1,d2,d3] = size(data{j});
        [n1,n2,n3] = size(output(:,:, idx1:idx2 ));

        if d1~=n1 || d2~=n2 || d3~=n3 %M(j)~=max(M) || idx2-idx1+1 ~= P(j)
            msg = sprintf('Warning: data{j} has a dimension of %d - %d - %d, while output(:,:,idx1:idx2) has a dimension of %d - %d - %d, %s\n', d1,d2,d3,n1,n2,n3, datestr(now));
            my_append('cell2mat_tensor3D_log.txt', msg);
            continue;
        end
        output(:,:, idx1:idx2 ) = data{j};
    end
end

end

author: singaxiong

repo: SignalGraph


text:
function [ctf] = ctf_read_badchannels(folder,ctf);

% ctf_read_badchannels - read BadChannels file from a CTF .ds folder
%
% [ctf] = ctf_read_badchannels(folder,ctf)
%
% This function reads the BadChannels text file in a .ds folder.
%
% The BadChannels file is a text file that contains the channel names of
% any bad channels.  The channel names are assumed to be listed in a single
% column, one channel per line.
%
% INPUTS
%
% If you do not wish to specify an input option, use [], but keep the order
% of the input options as above.  Only specify as many input options as
% required.  With no input options, the function will prompt for a folder
% and read the BadChannels file (if it exists).
%
% folder - the directory of the .ds data set to read.  By
% default, a gui prompts for the folder.
%
% ctf - a struct with folder, setup, sensor and data fields.
%
% OUTPUTS
%
% ctf.sensor.bad - a string cell array of channel labels
%
%      <>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> %
%      <                                                       > %
%      <                      DISCLAIMER:                      > %
%      <                                                       > %
%      < THIS PROGRAM IS INTENDED FOR RESEARCH PURPOSES ONLY.  > %
%      < THIS PROGRAM IS IN NO WAY INTENDED FOR CLINICAL OR    > %
%      <                     OFFICIAL USE.                     > %
%      <                                                       > %
%      <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<> %
%

% $Revision: 1.1 $ $Date: 2009-01-30 03:49:27 $

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

% Modified: 01/2005, Darren.Weber_at_radiology.ucsf.edu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%-----------------------------------------------
% ensure we have the folder

if exist('folder','var'),
  if exist('ctf','var'),
    ctf = ctf_folder(folder,ctf);
  else
    ctf = ctf_folder(folder);
  end
else
  if exist('ctf','var'),
    ctf = ctf_folder([],ctf);
  else
    ctf = ctf_folder;
  end
end

%--------------------------------------------------------------
ver = '$Revision: 1.1 $ $Date: 2009-01-30 03:49:27 $';
fprintf('\nCTF_READ_BADCHANNELS [v %s]\n',ver(11:15)); tic;

%----------------------------------------------------------------
% read the file

[folderPath,folderName,folderExt] = fileparts(ctf.folder);
BadChannelsFile = findBadChannelsFile( ctf.folder );

ctf.sensor.bad = textread(BadChannelsFile,'%s','delimiter','\n');

t = toc; fprintf('...done (%6.2f sec)\n\n',t);

return

% -------------------------------------------------------
function BadChannelsFile = findBadChannelsFile( folder )

filename = dir([ folder filesep 'BadChannels' ]);
if isempty(filename)
    filename = dir([ folder filesep 'badchannels' ]);
end;

if isempty(filename)
    error('No BadChannels file in selected folder');
else
    BadChannelsFile = [ folder filesep filename.name ];
end;

return

author: PatternRecognition

repo: OpenBMI