text:
function Population = LCSA_NSGAIIIEnvironmentalSelection(Population,N,Z,Zmin)
% -----------------------------------------------------------------------
%  Copyright (C) 2020 Heiner Zille
%
%  This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0
%  International License. (CC BY-NC-SA 4.0). To view a copy of this license,
%  visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or see the
%  pdf-file "License-CC-BY-NC-SA-4.0.pdf" that came with this code.
%
%  You are free to:
%  * Share ? copy and redistribute the material in any medium or format
%  * Adapt ? remix, transform, and build upon the material
%  Under the following terms:
%  * Attribution ? You must give appropriate credit, provide a link to the
%     license, and indicate if changes were made. You may do so in any reasonable
%     manner, but not in any way that suggests the licensor endorses you or your use.
%  * NonCommercial ? You may not use the material for commercial purposes.
%  * ShareAlike ? If you remix, transform, or build upon the material, you must
%    distribute your contributions under the same license as the original.
%  * No additional restrictions ? You may not apply legal terms or technological
%    measures that legally restrict others from doing anything the license permits.
%
%  Author of this Code:
%   Heiner Zille <heiner.zille@ovgu.de> or <heiner.zille@gmail.com>
%
%  This code is based on the following publications:
%
%  1) Heiner Zille
%     "Large-scale Multi-objective Optimisation: New Approaches and a Classification of the State-of-the-Art"
%     PhD Thesis, Otto von Guericke University Magdeburg, 2019
%     http://dx.doi.org/10.25673/32063
%
%  2) Heiner Zille and Sanaz Mostaghim
%     "Linear Search Mechanism for Multi- and Many-Objective Optimisation"
%     10th International Conference on Evolutionary Multi-Criterion Optimization (EMO 2019),
%        Lecture Notes in Computer Science, vol 11411.
%        Deb K. et al. (eds), Springer, Cham, East Lansing, Michigan, USA, March 2019
%     https://doi.org/10.1007/978-3-030-12598-1_32.
%
%  This file is intended to work with the PlatEMO framework version 2.5.
%  Date of publication of this code: 06.04.2020
%  Last Update of this code: 06.04.2020
%  A newer version of this algorithm may be available. Please contact the author
%  or see http://www.ci.ovgu.de/Research/Codes.html.
%
% The files may have been modified in Feb 2021 by the authors of the Platemo framework to work with the Platemo 3.0 release.
% -----------------------------------------------------------------------
% This file is derived from its original version containied in the PlatEMO
% framework.
% -----------------------------------------------------------------------

    if isempty(Zmin)
        Zmin = ones(1,size(Z,2));
    end

    %% Non-dominated sorting
    [FrontNo,MaxFNo] = NDSort(Population.objs,Population.cons,N);
    Next = FrontNo < MaxFNo;

    %% Select the solutions in the last front
    Last   = find(FrontNo==MaxFNo);
    Choose = LastSelection(Population(Next).objs,Population(Last).objs,N-sum(Next),Z,Zmin);
    Next(Last(Choose)) = true;
    % Population for next generation
    Population = Population(Next);
end

function Choose = LastSelection(PopObj1,PopObj2,K,Z,Zmin)
% Select part of the solutions in the last front

    PopObj = [PopObj1;PopObj2] - repmat(Zmin,size(PopObj1,1)+size(PopObj2,1),1);
    [N,M]  = size(PopObj);
    N1     = size(PopObj1,1);
    N2     = size(PopObj2,1);
    NZ     = size(Z,1);

    %% Normalization
    % Detect the extreme points
    Extreme = zeros(1,M);
    w       = zeros(M)+1e-6+eye(M);
    for i = 1 : M
        [~,Extreme(i)] = min(max(PopObj./repmat(w(i,:),N,1),[],2));
    end
    % Calculate the intercepts of the hyperplane constructed by the extreme
    % points and the axes
    Hyperplane = PopObj(Extreme,:)\ones(M,1);
    a = 1./Hyperplane;
    if any(isnan(a))
        a = max(PopObj,[],1)';
    end
    % Normalization
    PopObj = PopObj./repmat(a',N,1);

    %% Associate each solution with one reference point
    % Calculate the distance of each solution to each reference vector
    Cosine   = 1 - pdist2(PopObj,Z,'cosine');
    Distance = repmat(sqrt(sum(PopObj.^2,2)),1,NZ).*sqrt(1-Cosine.^2);
    % Associate each solution with its nearest reference point
    [d,pi] = min(Distance',[],1);

    %% Calculate the number of associated solutions except for the last front of each reference point
    rho = hist(pi(1:N1),1:NZ);

    %% Environmental selection
    Choose  = false(1,N2);
    Zchoose = true(1,NZ);
    % Select K solutions one by one
    while sum(Choose) < K
        % Select the least crowded reference point
        Temp = find(Zchoose);
        Jmin = find(rho(Temp)==min(rho(Temp)));
        j    = Temp(Jmin(randi(length(Jmin))));
        I    = find(Choose==0 & pi(N1+1:end)==j);
        % Then select one solution associated with this reference point
        if ~isempty(I)
            if rho(j) == 0
                [~,s] = min(d(N1+I));
            else
                s = randi(length(I));
            end
            Choose(I(s)) = true;
            rho(j) = rho(j) + 1;
        else
            Zchoose(j) = false;
        end
    end
end


text:
function b = dspsl ( ap, n, kpvt, b )

%*****************************************************************************80
%
%% DSPSL solves the real symmetric system factored by DSPFA.
%
%  Discussion:
%
%    To compute inverse(A) * C where C is a matrix with P columns:
%
%      call dspfa ( ap, n, kpvt, info )
%
%      if ( info /= 0 ) go to ...
%
%      do j = 1, p
%        call dspsl ( ap, n, kpvt, c(1,j) )
%      end do
%
%    A division by zero may occur if DSPCO has set RCOND == 0.0D+00
%    or DSPFA has set INFO /= 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real ( kind = 8 ) AP(N*(N+1)/2), the output from DSPFA.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer KPVT(N), the pivot vector from DSPFA.
%
%    Input, real B(N), the right hand side.
%
%    Output, real B(N), the solution.
%

%
%  Loop backward applying the transformations and D inverse to B.
%
  k = n;
  ik = floor ( ( n * ( n - 1 ) ) / 2 );

  while ( 0 < k )

    kk = ik + k;

    if ( 0 <= kpvt(k) )
%
%  1 x 1 pivot block.
%
      if ( k ~= 1 )

        kp = kpvt(k);
%
%  Interchange.
%
        if ( kp ~= k )
          temp = b(k);
          b(k) = b(kp);
          b(kp) = temp;
        end
%
%  Apply the transformation.
%
        b(1:k-1) = daxpy ( k-1, b(k), ap(ik+1:ik+k-1), 1, b(1:k-1), 1 );

      end
%
%  Apply D inverse.
%
      b(k) = b(k) / ap(kk);
      k = k - 1;
      ik = ik - k;

    else
%
%  2 x 2 pivot block.
%
      ikm1 = ik - ( k - 1 );

      if ( k ~= 2 )

        kp = abs ( kpvt(k) );
%
%  Interchange.
%
        if ( kp ~= k-1 )
          temp = b(k-1);
          b(k-1) = b(kp);
          b(kp) = temp;
        end
%
%  Apply the transformation.
%
        b(1:k-2) = daxpy ( k-2, b(k), ap(ik+1:ik+k-2), 1, b(1:k-2), 1 );
        b(1:k-2) = daxpy ( k-2, b(k-1), ap(ikm1+1:ikm1+k-2), 1, b(1:k-2), 1 );

      end
%
%  Apply D inverse.
%
      km1k = ik + k - 1;
      kk = ik + k;
      ak = ap(kk) / ap(km1k);
      km1km1 = ikm1 + k - 1;
      akm1 = ap(km1km1) / ap(km1k);
      bk = b(k) / ap(km1k);
      bkm1 = b(k-1) / ap(km1k);
      denom = ak * akm1 - 1.0;
      b(k) = ( akm1 * bk - bkm1 ) / denom;
      b(k-1) = ( ak * bkm1 - bk ) / denom;
      k = k - 2;
      ik = ik - ( k + 1 ) - k;

    end

  end
%
%  Loop forward applying the transformations.
%
  k = 1;
  ik = 0;

  while ( k <= n )

    if ( 0 <= kpvt(k) )
%
%  1 x 1 pivot block.
%
      if ( k ~= 1 )
%
%  Apply the transformation.
%
        b(k) = b(k) + ddot ( k-1, ap(ik+1:ik+k-1), 1, b(1:k-1), 1 );
        kp = kpvt(k);
%
%  Interchange.
%
        if ( kp ~= k )
          temp = b(k);
          b(k) = b(kp);
          b(kp) = temp;
        end

      end

      ik = ik + k;
      k = k + 1;

    else
%
%  2 x 2 pivot block.
%
      if ( k ~= 1 )
%
%  Apply the transformation.
%
        b(k) = b(k) + ddot ( k-1, ap(ik+1:ik+k-1), 1, b(1:k-1), 1 );
        ikp1 = ik + k;
        b(k+1) = b(k+1) + ddot ( k-1, ap(ikp1+1:ikp1+k-1), 1, b(1:k-1), 1 );
        kp = abs ( kpvt(k) );
%
%  Interchange.
%
        if ( kp ~= k )
          temp = b(k);
          b(k) = b(kp);
          b(kp) = temp;
        end

      end

      ik = ik + k + k + 1;
      k = k + 2;

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function turk = InitializeTurk(aws_access_key, aws_secret_key, sandbox)

if ~exist('sandbox','var')
    sandbox = true;
end

turk.aws_access_key = aws_access_key;
turk.aws_secret_key = aws_secret_key;

if sandbox
    turk.service_url='https://mechanicalturk.sandbox.amazonaws.com/?Service=AWSMechanicalTurkRequester';
else
    turk.service_url='https://mechanicalturk.amazonaws.com/?Service=AWSMechanicalTurkRequester';
end

% Operation List
disp('Operation List:');
disp('http://docs.amazonwebservices.com/AWSMechTurk/latest/AWSMturkAPI/ApiReference_OperationsArticle.html');

author: jianxiongxiao

repo: ProfXkit


text:
% Copyright (C) 2010 Joao Carreira
%
% This code is part of the extended implementation of the paper:
%
% J. Carreira, C. Sminchisescu, Constrained Parametric Min-Cuts for Automatic Object Segmentation, IEEE CVPR 2010
%
% Modified:
%   @authors:     Ahmad Humayun
%   @contact:     ahumayun@cc.gatech.edu
%   @affiliation: Georgia Institute of Technology
%   @date:        Fall 2013 - Summer 2014

%function pixel_ids = frame_pixel_ids(nrows, ncols, width, custom)
% creates a frame whose external border has nrows and ncols, and width = 'width'
function seeds = frame_pixel_ids(nrows, ncols, width, custom)
    assert(width>=1);

    boundary_pixels_horiz_top = false(nrows, ncols);
    boundary_pixels_horiz_top(1:width, 1:end) = true;
    boundary_pixels_horiz_bottom = false(nrows, ncols);
    boundary_pixels_horiz_bottom(end-width+1:end, 1:end) = true;
    boundary_pixels_vert_left = false(nrows, ncols);
    boundary_pixels_vert_left(1:end, 1:width) = true;
    boundary_pixels_vert_right = false(nrows, ncols);
    boundary_pixels_vert_right(1:end, end-width+1:end) = true;

    if(strcmp(custom, 'horiz'))
        seeds = boundary_pixels_horiz_bottom | boundary_pixels_horiz_top;
    elseif(strcmp(custom, 'down'))
        seeds = boundary_pixels_horiz_bottom;
    elseif(strcmp(custom, 'up'))
        seeds = boundary_pixels_horiz_top;
    elseif(strcmp(custom, 'vert'))
        seeds = boundary_pixels_vert_left | boundary_pixels_vert_right;
    elseif(strcmp(custom, 'all_but_down'))
        seeds = boundary_pixels_horiz_top | boundary_pixels_vert_left | ...
                boundary_pixels_vert_right;
    elseif(strcmp(custom, 'all')) % whole frame
        seeds = boundary_pixels_horiz_top | ...
                boundary_pixels_horiz_bottom | ...
                boundary_pixels_vert_left | boundary_pixels_vert_right;
    else
        error('no such option');
    end

    seeds = seeds(:);
end

author: Cloud-CV

repo: object-proposals


text:
figure;
% Load dataset
load('maneuvering_robot.mat');
truth = [TrueTrack.Trajectory.Vector];
obs = LinearGaussianX('NumMeasDims',2,'NumStateDims',4,'MeasurementErrVariance',100^2,'Mapping',[1 3]);
%measurement_model = RangeBearing2CartesianX('NumStateDims',4,'MeasurementErrVariance',[0.001,0.02],'Mapping',[1 3]);

%% Cartesian
% subplot(1,2,1);
hold on;
truth_100 = truth.*100;
%h1 = plot(ax(1),measurement(1,k),measurement(2,k),'k*','MarkerSize', 10);
measurement = obs.feval(truth_100, true);
meas = obs.finv([measurement(1,:);measurement(2,:)]);
h1 = plot(meas(obs.Mapping(1),:),meas(obs.Mapping(2),:),'r+','MarkerSize', 15);
h2 = plot(truth_100(1,1:end),truth_100(3,1:end),'k--','LineWidth',2);
h3 = plot(truth_100(1,1),truth_100(3,1),'ko','MarkerSize', 20, 'MarkerFaceColor','Green');
h4 = plot(truth_100(1,end),truth_100(3,end),'ko','MarkerSize', 20, 'MarkerFaceColor','Red');
h5 = plot(0,0,'rd','MarkerSize', 20, 'MarkerFaceColor','Blue');
legend([h1,h2,h3,h4,h5],'Measurements','Trajectory', 'Start', 'End','Radar','Location','southeast')
str = sprintf('Vessel trajectory');
title(str)
xlabel('X (m)')
ylabel('Y (m)')
axis([0,2500,0,1500])
box on

%% Polar
% subplot(1,2,2);
% truth2meas = obs.heval(truth_100(1:2,:));
% % [a,b]=obs.heval(truth_100(1:2,:));
% %h1 = polarplot(a(x_start),b(y_start),'ko','LineWidth',2,'MarkerSize',15,'MarkerFaceColor','w');hold on;
% polarplot(truth2meas(1,:),truth2meas(2,:),'k-','LineWidth',2,'MarkerSize',10,'MarkerFaceColor','w');hold on;
% thetalim([0 90])
% %h2 = polarplot(a(x_end),b(y_end),'k^','LineWidth',2,'MarkerSize',15,'MarkerFaceColor','w');hold on;
% h3 = polarplot(measurement(1,:),100*measurement(2,:),'+r','MarkerSize',15);
% pax = gca;
% pax.ThetaAxisUnits = 'radians';
% %rlabel('Range (m)');
% %thetalabel('Bearing (rad)');

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function colormap = BF_MakeBrightenedColorMap(rgb,numGrads)
% BF_MakeBrightenedColorMap  Given an rgb vector, progressively brightens into a colormap

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

if nargin < 2
    numGrads = 8;
end

colormap = zeros(numGrads,3);

brightenRange = linspace(0,1,numGrads);
for i = 1:numGrads
    colormap(i,:) = brighten(rgb,brightenRange(i));
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [MATreordered,MATindices,MATcost] = reorderMAT(MAT,H,cost)
%REORDERMAT         Reorder matrix for visualization
%
%   [MATreordered,MATindices,MATcost] = reorderMAT(MAT,H,cost);
%
%   This function reorders the connectivity matrix in order to place more
%   edges closer to the diagonal. This often helps in displaying community
%   structure, clusters, etc.
%
%   Inputs:     MAT,            connection matrix
%               H,              number of reordering attempts
%               cost,           'line' or 'circ', for shape of lattice
%                               (linear or ring lattice)
%
%               MATreordered    reordered connection matrix
%               MATindices      reordered indices
%               MATcost         cost of reordered matrix
%
%
%   Olaf Sporns, Indiana University

N = length(MAT);
diagMAT = diag(diag(MAT));
MAT = MAT-diagMAT;

% generate cost function
if strcmp(cost,'line')
    profil = fliplr(normpdf(1:N,0,N/2));
end;
if strcmp(cost,'circ')
    profil = fliplr(normpdf(1:N,N/2,N/4));
end;
COST = toeplitz(profil,profil);

% initialize lowCOST
lowMATcost = sum(sum(COST.*MAT));

% keep track of starting configuration
MATstart = MAT;
starta = 1:N;

% reorder
for h=1:H
    a = 1:N;
    % choose two positions at random and flip them
    r = randperm(N);
    a(r(1)) = r(2);
    a(r(2)) = r(1);
    MATcostnew = sum(sum(MAT(a,a).*COST));
    if (MATcostnew < lowMATcost)
        MAT = MAT(a,a);
        r2 = starta(r(2));
        r1 = starta(r(1));
        starta(r(1)) = r2;
        starta(r(2)) = r1;
        lowMATcost = MATcostnew;
    end;
end;	% h

MATreordered = MATstart(starta,starta) + diagMAT(starta,starta);
MATindices = starta;
MATcost = lowMATcost;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function varargout = gammainc(varargin)

if nargin ~= 2
    error('Not enough input arguments in gammainc.');
end

if isa(varargin{1},'sdpvar') && isa(varargin{2},'sdpvar')
    error('Only one argument in gammainc(X,A) can be an SDPVAR');
end

switch class(varargin{1})
    case 'double'

        if varargin{1} < 0
            error('X must be real and non-negative in gammainc(X,A)')
        end
        varargout{1} = InstantiateElementWise('gammainc_a',varargin{2:-1:1});

    case 'sdpvar'

        if varargin{2}<0
            error('A must be real and non-negative gammainc(X,A)');
        end
        varargout{1} = InstantiateElementWise('gammainc_x',varargin{:});

    otherwise
        error(['SDPVAR/' upper(mfilename) ' called with weird argument']);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
classdef LevelSetWithCircleInclusion < LevelSetSphereNdim

    methods (Access = public)

        function obj = LevelSetWithCircleInclusion(cParams)
            obj.fracRadius = cParams.fracRadius;
            obj.compute(cParams);
        end
    end

    methods (Access = protected)
        function computeLevelSetValue(obj)
            ls = 1 - obj.dist;
            obj.levelSet = ls;
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [DCM] = VBA_spm_sdcm_estimate2(P)
% Estimate parameters of a stochastic DCM for fMRI data (deconv HRF)
% FORMAT [DCM] = spm_sdcm_estimate2(DCM)
%
% DCM  - the DCM or its filename
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Jean Daunizeau

% load DCM structure
if ~nargin
    P = VBA_spm_select(1,'^DCM.*\.mat$','select DCM_???.mat');
    if isempty(P)
        return
    end
end
if isstruct(P)
    DCM = P;
    P   = ['DCM-' date];
else
    load(P)
end

disp(' ')
disp('-----------------------------')
disp('Inverting stochastic DCM...')

% Get data and remove confounds
y = DCM.Y.y';
X0 = DCM.Y.X0;
iX0 = pinv(X0'*X0)*X0';
nreg = size(y,1);
for i=1:nreg
    beta = iX0*y(i,:)';
    yc = X0*beta;
    y(i,:) = y(i,:) - yc';
end

% Unpack DCM specification
dt = DCM.U.dt;
uu = DCM.U.u';
[nu,nt] = size(uu);
[p,ny] = size(y);
A = DCM.a - eye(nreg);
for i=1:nu
    B{i} = DCM.b(:,:,i);
end
C = DCM.c;
if isfield(DCM,'d')
    for i=1:nreg
        try
            D{i} = DCM.d(:,:,i);
        catch
            D{i} = zeros(nreg,nreg);
        end
    end
else
    for i=1:nreg
        D{i} = zeros(nreg,nreg);
    end
end

% Prepare optional input for inversion routine
TR = DCM.Y.dt;
f_fname = @f_dcm4fmri;
g_fname = @g_Id;
microDT = 1e-1;             % micro-time resolution (in sec)
[inF] = prepare_dcm(A,B,C,D);
options.decim = max([1,floor(TR./microDT)]);
inF.deltat = TR./options.decim;
options.inF = inF;
options.microU = 1;
dim.n_theta = inF.indself;
dim.n_phi = 0;
dim.n = nreg;
dim.p = p;

% Build priors
alpha = 1e2/TR;
[priors] = getPriors(nreg,dim,alpha,ny);
options = getOptions(options,priors);

% Resample inputs on microtime integration grid
[u] = resampleU(uu,dt,nt,nu,ny,options);
u = [zeros(nu,1),u]; % add initial conditions

% Deconvolve fMRI data from HRF
[hrf] = VBA_spm_hrf(TR);
y0 = y;
y = [];
SNR = VBA_spm_getSNR(y0,uu,hrf);
for i=1:nreg
    try
        D = round(DCM.delays(i)./TR);
        y0i = [zeros(1,D),y0(i,1:end-D)];
    catch
        y0i = y0(i,:);
    end
    y(i,:) = VBA_spm_deconv(y0i,hrf,SNR);
end

% invert stochastic DCM using discrete-time VB routine:
[posterior,out] = VBA_NLStateSpaceModel(y,u,f_fname,g_fname,dim,options);

% extract relevant info and fills in DCM structure
[DCM] = exportDCMfromVBNLSS(posterior,out,DCM);

% change DCM output (account for HRF convolution)
for i=1:nreg
    y0i = conv(out.suffStat.gx(i,:),hrf);
    DCM.y(:,i) = y0i(1:end-length(hrf)+1)'; % predicted data
    DCM.R(:,i) = y0(i,:)' - DCM.y(:,i); % residuals of the model
end

disp('Inverting stochastic DCM... OK')
disp('-----------------------------')
disp(' ')

return

function [u] = resampleU(uu,dt,nt,nu,ny,options)
microDT = options.inF.deltat;
grid1 = 0:dt:dt*(nt-1);
grid2 = 0:microDT:microDT*options.decim*ny;
u = zeros(nu,ny);
for i=1:length(grid2)-1
    [tmp,ind1] = min(abs(grid2(i)-grid1));
    [tmp,ind2] = min(abs(grid2(i+1)-grid1));
    u(:,i) = mean(uu(:,ind1:ind2),2);
end
% [u,alpha] = spm_resample(full(uu),dt/microDT);

function [priors] = getPriors(nreg,dim,alpha,n_t)
priors.muX0 = 0.*ones(nreg,1);
priors.SigmaX0 = 0e-3*eye(nreg);
priors.muTheta = 0*ones(dim.n_theta,1);
priors.SigmaTheta = 1e-2*eye(dim.n_theta);
priors.muPhi = [];
priors.SigmaPhi = [];
% priors.SigmaTheta(options.inF.indself,options.inF.indself) = 0;

% NB on hyperpriors:
%   - fix state noise precision using high scale param
%   - use non-informative priors on the residual precision, with high
%   expectation.
% This is because of the first iteration of the hidden states posterior
% update, which has to deviate from the its prior predictive density (as
% derived from the deterministic inversion).
% The following iterations will then work with a realistic (expected)
% residual precision, and adapt.
SC = 1e0;
priors.a_alpha = SC*alpha;
priors.b_alpha = SC;
priors.a_sigma = 1e0;
priors.b_sigma = 1e-4;
for t = 1:n_t
    priors.iQx{t} = eye(dim.n);
    priors.iQy{t} = eye(dim.p);
end

function options = getOptions(options,priors)
options.priors = priors;
options.DisplayWin = 1;
options.GnFigs = 0;
options.gradF = 0;
options.updateHP = 1;
options.backwardLag = 2;
options.Laplace = 0;
% options.noSXi = 1;
% options.init0 = 0;
% options.embed = 0;

author: MBB-team

repo: VBA-toolbox


text:
function [q,r] = deconv(p1,p2)
%DECONV       Implements  p1 / p2  for univariate polynomials with remainder term
%
%   [q,r] = deconv(p1,p2)
%

% written  07/24/02     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if ~isa(p1,'polynom')
    p1 = polynom(p1,'');
  end

  if ~isa(p2,'polynom')
    p2 = polynom(p2,'');
  end

  if ( size(p1.e,2)>1 ) | ( size(p2.e,2)>1 )
    error('divison only for univariate polynomials')
  end

  % both p1 and p2 univariate polynomials
  if ~isequal(p1.v,p2.v) & ~isempty(p1.v) & ~isempty(p2.v)
    error('division only for univariate polynomials depending on the same variable')
  end

  isint = isa(p1.c,'intval') | isa(p2.c,'intval');

  if p2.e==0                           % p2 is constant
    q = p1;
    q.c = p1.c/p2.c;
    if nargout==2
      if isint
        r = polynom(intval(0),p1.v);
      else
        r = polynom(0,p1.v);
      end
    end
    if rndold
      setround(rndold)
    end
    return
  end

  if p1.e<p2.e                          % deg(p1) < deg(p2)
    q.e = 0;
    if isint
      q.c = intval(0);
    else
      q.c = 0;
    end
    q.v = p1.v;
    q = class(q,'polynom');
    if nargout==2
      r = p1;
    end
    if rndold
      setround(rndold)
    end
    return
  end

  q.e = max(p1.e-p2.e,0);
  if isint
    q.c = intval(zeros(1,q.e+1));
    p1.c = intval(p1.c);
    for i=1:(q.e+1)
      q.c(i) = p1.c(i)/p2.c(1);
      p1.c(i:(i+p2.e)) = p1.c(i:(i+p2.e)) - q.c(i)*p2.c;
    end
    if nargout==2
      np1 = length(p1.c);
      r.e = p2.e-1;
      r.c = p1.c((np1-p2.e+1):np1);
      r.v = p1.v;
      r = class(r,'polynom');
    end
  else
    if nargout==2
      [q.c,rc] = deconv(p1.c,p2.c);
      r.e = min(p2.e-1,length(rc)-1);
      r.c = rc(end-r.e:end);
      r.v = p1.v;
      r = class(r,'polynom');
    else
      q.c = deconv(p1.c,p2.c);
    end
  end
  q.v = p1.v;

  q = class(q,'polynom');

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function [empty,c] = emptyintersect(a,b)
%EMPTYINTERSECT    Compute intersection and check for empty components
%
%   [empty,c] = emptyintersect(a,b)
%
%Result c is that of intersect(a,b), and
%  empty(i) = 1     intersection of a(i) and b(i) is empty
%             0     intersection of a(i) and b(i) is not empty
%             NaN   at least one of a(i) and b(i) is NaN
%
%Input a and b must be both real or both complex
%

% written  04/22/09     S.M. Rump
% written  05/14/09     S.M. Rump  identical midpoints
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if prod(size(a))>1
    if prod(size(b))==1
      if issparse(a)
        b = b*spones(a);
      else
        b = b*ones(size(a));
      end
    end
  else
    if prod(size(b))>1
      if issparse(b)
        a = a*spones(b);
      else
        a = a*ones(size(b));
      end
    end
  end

  if ~isequal(size(a),size(b))
    error('intersect called with non-matching dimensions')
  end

  if ~isa(a,'intval')
    a = intval(a);
  end
  if ~isa(b,'intval')
    b = intval(b);
  end

  if a.complex & b.complex

    c.complex = 1;
    c.inf = [];
    c.sup = [];
    c.mid = b.mid;
    c.rad = b.rad;
    empty = zeros(size(b.mid));

    v = b.mid - intval(a.mid);         % connection of midpoints
    d2 = sqr(intval(real(v))) + sqr(intval(imag(v)));
    d = sqrt(d2);                      % inclusion of distance of midpoints
    arad = intval(a.rad);
    Delta = (arad+b.rad).*(arad-b.rad);
    wng = warning;
    warning off
    sumrad = arad+b.rad;
    index = ( d.inf > sumrad.sup );    % empty intersection
    if any(index(:))
      empty(index) = 1;
      c.mid(index) = NaN;
      c.rad(index) = NaN;
    end
    index0 = ( Delta.inf<=-d2 );
    if any(index0(:))                  % diameter of a in intersection
      c.mid(index0) = a.mid(index0);
      c.rad(index0) = a.rad(index0);
    end
    index1 = ( Delta.sup>=d2 );
    if any(index1(:))                  % diameter of b in intersection
      c.mid(index1) = b.mid(index1);
      c.rad(index1) = b.rad(index1);
    end
    index = ~( index | index0 | index1 );
    if any(index(:))
      x = ( 1 + (arad+b.rad).*(arad-b.rad)./d2 )/2;
      %VVVV  cmid = a.mid(index) + x(index).*v(index);
      %VVVV  c.mid(index) = mid(cmid);
      %VVVV  c.rad(index) = sup( rad(cmid) + sqrt(sqr(arad(index))-sqr(x(index).*d(index))) );
      s.type = '()'; s.subs = {index};
      cmid = a.mid(index) + subsref(x,s).*subsref(v,s);
      c.mid = subsasgn(c.mid,s,mid(cmid));
      c.rad = subsasgn(c.rad,s,sup( rad(cmid) + sqrt(sqr(subsref(arad,s))-sqr(subsref(x,s).*subsref(d,s))) ));
      %AAAA  Matlab bug fix
    end
    warning(wng);

    index = ( imag(c.rad)~=0 );
    if any(index(:))
      empty(index) = 1;
      c.mid(index) = NaN;
      c.rad(index) = NaN;
    end
    index = isnan(a.mid) | isnan(a.rad) | isnan(b.mid) | isnan(b.rad);
    if any(index(:))
      empty(index) = NaN;
      c.mid(index) = NaN;
      c.rad(index) = NaN;
    end

  elseif ~a.complex & ~b.complex    % two real intervals

    c.complex = 0;
    c.inf = max(a.inf,b.inf);
    c.sup = min(a.sup,b.sup);
    empty = zeros(size(a.inf));
    index = ( c.inf>c.sup );        % empty intersection
    if any(index(:))
      empty(index) = 1;
      c.inf(index) = NaN;
      c.sup(index) = NaN;
    end
    index = isnan(a.inf) | isnan(a.sup) | isnan(b.inf) | isnan(b.sup);
    if any(index(:))
      empty(index) = NaN;
      c.inf(index) = NaN;
      c.sup(index) = NaN;
    end
    c.mid = [];
    c.rad = [];

  else
    error('operands of intersect must be both real or both complex')
  end

  c = class(c,'intval');

  setround(rndold)


text:
%Script to ectract heterogenity metrics from MR images
%
%APA, 03/16/2010

global planC
indexS = planC{end};
structNum = 1;
disp('Heterogenity Metrics for MR')
[energy_40base,contrast_40base,Entropy_40base,Homogeneity_40base,standard_dev_40base,Ph_40base,slope_base] = getHaralicParams(structNum);
[jnk,jnk,raw] = xlsread('Cervix_MR_data.xls');
rowIndex = size(raw,1)+1;
dataM = [slope_base,energy_40base,contrast_40base,Entropy_40base,Homogeneity_40base,standard_dev_40base];
PtName = planC{indexS.scan}(1).scanInfo(1).patientName;
[SUCCESS,MESSAGE] = xlswrite('Cervix_MR_data.xls',{PtName},['A',num2str(rowIndex),':A',num2str(rowIndex)]);
if ~SUCCESS
    errordlg(MESSAGE.message,'Error Writing to Excel','modal')
    return
end
[SUCCESS,MESSAGE] = xlswrite('Cervix_MR_data.xls',{datestr(now)},['B',num2str(rowIndex),':B',num2str(rowIndex)]);
if ~SUCCESS
    errordlg(MESSAGE.message,'Error Writing to Excel','modal')
    return
end
acquisitionDate = datestr(datenum(num2str(planC{indexS.scan}(1).scanInfo(1).DICOMHeaders.AcquisitionDate),'yyyymmdd'));
[SUCCESS,MESSAGE] = xlswrite('Cervix_MR_data.xls',{acquisitionDate},['D',num2str(rowIndex),':D',num2str(rowIndex)]);
if ~SUCCESS
    errordlg(MESSAGE.message,'Error Writing to Excel','modal')
    return
end
[SUCCESS,MESSAGE] = xlswrite('Cervix_MR_data.xls',dataM,['E',num2str(rowIndex),':J',num2str(rowIndex)]);
if ~SUCCESS
    errordlg(MESSAGE.message,'Error Writing to Excel','modal')
    return
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% PURPOSE: subroutine for ploterpGUI.m
%          identifies minimum and maximum values of ERP amplitudes for Y scaling.
%
% FORMAT
%
%  [yylim, serror] = erpAutoYLim(ERP, binArray, chanArray, xxlim)
%
% INPUTS
%
% ERP         - ERPset
% binArray    - indices of bins from where to get the amplitude values
% chanArray   - indices of channels from where to get the amplitude values
% xxlim       - current scale for time ([min max] in ms)
%
%
% OUTPUT
%
% yylim       - range for Y scale
% serror      - error flag. 0 means no errors.
%
% *** This function is part of ERPLAB Toolbox ***
% Author: Javier Lopez-Calderon & Steven Luck
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2013

function [yylim, serror] = erpAutoYLim(ERP, binArray, chanArray, xxlim,blcorrdata)
%ams updated blcorrdata
serror = 0;
if nargin<1
        error('erpAutoYLim needs 1 input argument at least.')
end

datatype = checkdatatype(ERP);

if nargin<5
    blcorrdata = 'no';
end

if nargin<4
        if strcmpi(datatype, 'ERP')
                xxlim = [ERP.xmin ERP.xmax]*1000;
        else
                xxlim = [ERP.xmin ERP.xmax];
        end
end
if nargin<3
        chanArray = 1:ERP.nchan;
end
if nargin<2
        binArray = 1:ERP.nbin;
end
if isempty(binArray)
        binArray = 1:ERP.nbin;
end
if isempty(chanArray)
        chanArray = 1:ERP.nchan;
end
if isempty(xxlim)
    if strcmpi(datatype, 'ERP')
        if strcmp(datatype,'CSD')
            datatype = 'ERP';     % if CSD data, treat as ERP data here
        end
        xxlim = [ERP.xmin ERP.xmax]*1000;
    else
        xxlim = [ERP.xmin ERP.xmax];
    end
end
try
        nbin  = length(binArray);
        nchan = length(chanArray);
        fs    = ERP.srate;

        if strcmpi(datatype,'CSD')
            datatype = 'ERP';     % if CSD data, treat as ERP data here
        end

        if strcmpi(datatype, 'ERP')

            if xxlim(1)<round(ERP.xmin*1000)
                aux_xxlim(1) = round(ERP.xmin*1000);
            else
                aux_xxlim(1) = xxlim(1);
            end
            if xxlim(2)>round(ERP.xmax*1000)
                aux_xxlim(2) = round(ERP.xmax*1000);
            else
                aux_xxlim(2) = xxlim(2);
            end
        else  % fft
            if xxlim(1)<5
                aux_xxlim(1) = 5; % to avoid including the spectrum under 5 Hz in calculating Y auto (too big!)
            else
                aux_xxlim(1) = xxlim(1);
            end
            if xxlim(2)>round(fs/2);
                aux_xxlim(2) = round(fs/2);
            else
                aux_xxlim(2) = xxlim(2);
            end
        end

        %disp('start')
        %aux_xxlim
        %[ERP.xmin ERP.xmax]
        %fs
        [p1, p2, checkw] = window2sample(ERP, aux_xxlim(1:2) , fs, 'relaxed');
        %disp('end')

        if strcmpi(blcorrdata,'no')
            datresh = reshape(ERP.bindata(chanArray,p1:p2,binArray), 1, (p2-p1+1)*nbin*nchan);
            yymax   = max(datresh);
            yymin   = min(datresh);
        else
            datresh = reshape(blcorrdata(chanArray,p1:p2,binArray), 1, (p2-p1+1)*nbin*nchan);
            yymax   = max(datresh);
            yymin   = min(datresh);
        end

        if abs(yymax)<1 && abs(yymin)<1
            yylim(1:2) = [yymin*1.2 yymax*1.1]; % JLC. Mar 11, 2015
        else
            yylim(1:2) = round([yymin*1.2 yymax*1.1]); % JLC. Sept 26, 2012
        end

        %
        % in case of flatlined ERPs
        %
        if yylim(1)<1E-6 && yylim(2)<1E-6
                if strcmpi(datatype, 'ERP')
                        yylim(1:2) = [-1 1];
                        fprintf('WARNING: It seems like erpAutoYLim() found flatlined ERPs. So auto Y-limit was set to [-1 1].\n');
                else
                        yylim(1:2) = [0 1];
                        fprintf('WARNING: It seems like erpAutoYLim() found flatlined Spectrum. So auto Y-limit was set to [0 1].\n');
                end
        end
catch
        if strcmpi(datatype, 'ERP')
                yylim(1:2) = [-10 10];
        else
                yylim(1:2) = [0 1];
        end
        serror =1;
        fprintf('WARNING: ERPLAB could not find the auto Y limits for %s.\nPlease check your input parameters and waverforms.\n', ERP.erpname);
        fprintf('Default Y limit values were loaded.\n');
        return
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%	[velo2, tmax]= SINCHRONIZE(qini, qfinal, velocity) Finds a mean speed and the required
%   time to perform a movement between the joint coordinates qini and qfinal.
%   If the speed of each joint is different, the maximum time to perform the movement
%   by the slower joint is taken as a basis.
%
%   Inputs:
%		Qini: initial position in joint coordinates.
%		Qfinal: final position in joint coordinates.
%		Velocity: stores the maximum velocity of each joint.
%   Outputs:
%        velo2: new maximum speed for each joint.
%        tmax: time needed to perform the movement.
%
%	See also: MOVEJ, COMPUTE_JOINT_TRAJECTORY_INDEP
%
%   Author: Arturo Gil
%   Date:   29/03/2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function [actual_speed, maxtime]=synchronize(qini, qfinal, speed, accel)

tacel = speed./accel;

tcte = (abs(qfinal(:)-qini(:))-accel(:).*tacel.^2)./speed(:);

time_total = tcte + 2*tacel;

maxtime=max(time_total);

actual_speed = (qfinal(:)-qini(:)-accel(:).*tacel.^2)/maxtime(:);

repo: ARTE


text:
close all
zebra = imread('zebra.bmp');
figure; imagesc(zebra); title('original image');

[Zebra1] = Nonlinear_Diffusion(double(zebra),1,1e-2,1,20,0, 0);
title('total variation flow(p=1, tau=1), 20 time steps');
[Zebra2] = Nonlinear_Diffusion(double(zebra),1,1e-3,1.2,20, 0,0);
title('edge enhacing flow(p=1.2, tau=1), 20 time steps');
[Zebra3] = Nonlinear_Diffusion(double(zebra),3,1e-3,1.2,30, 0,0);
title('edge enhacing flow(p=1.2, tau=3), 30 time steps');
[Zebra4] = Nonlinear_Diffusion(double(zebra),3,1e-3,1.3,30, 0,0);
title('edge enhacing flow(p=1.3, tau=3), 30 time steps');

house = imread('house_noisy.bmp');
figure; imagesc(house); title('original image');

[House1] = Nonlinear_Diffusion(double(house),1,1e-3,1,50, 0,0);
title('total variation flow(p=1, tau=1) , 50 time steps');

[House2] = Nonlinear_Diffusion(double(house),5,1e-3,1,10, 0,0);
title('total variation flow(p=1, tau=5) , 10 time steps');

[House4] = Nonlinear_Diffusion(double(house),5,1e-3,1,20, 0,0.5);
title('total variation flow(p=1, tau=5), 20 time steps, regularization sigma=0.5');

[House5] = Nonlinear_Diffusion(double(house),5,1e-3,1,20, 0,1);
title('total variation flow(p=1, tau=5), 20 time steps, regularization sigma=1');

[House6] = Nonlinear_Diffusion(double(house),5,1e-3,1.2,40, 0,0);
title('edge ehnacing flow(p=1.2, tau=5), 40 time steps');

[House7] = Nonlinear_Diffusion(double(house),5,1e-3,1.4,60, 0,0.4);
title('edge ehnacing flow(p=1.4, tau=5), 60 time steps, regularization sigma=0.4');

figure;
subplot(2,2,1); imagesc(zebra); title('original image');
subplot(2,2,2); imagesc(Zebra4/255); title('filtered image');
subplot(2,2,3); imagesc(house); title('original image');
subplot(2,2,4); imagesc(House7/255); title('filtered image');

author: Sable

repo: mcbench-benchmarks


text:
%PRKMEANS PRTools k-means clustering
%
%   [LABELS,B] = PRKMEANS(A,K,MAXIT,INIT)
%
% INPUT
%  A       Matrix or dataset
%  K       Number of clusters to be found (optional; default: 2)
%  MAXIT   maximum number of iterations (optional; default: 50)
%  INIT    Labels for initialisation, or
%          'rand'     : take at random K objects as initial means, or
%          'kcentres' : use KCENTRES for initialisation (default)
%
% OUTPUT
%  LABELS  Cluster assignments, 1..K
%  B       Dataset with original data and labels LABELS:
%          B = PRDATASET(A,LABELS)
%
% DESCRIPTION
% K-means clustering of data vectors in A.
%
% SEE ALSO (PRTools Guide)
% DATASETS, HCLUST, KCENTRES, MODESEEK, EMCLUST, PRPROGRESS

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

function [assign,a] = prkmeans(a,varargin)

  [kmax,maxit,init] = setdefaults(varargin,2,50,'kcentres');
	n_ini = 100;				% Maximum size of subset to use for initialisation.

	% Create dataset with all equal labels and no priors.
	m = size(a,1);
	a = prdataset(a);
	islabtype(a,'crisp');
  a=set(a,'labels',ones(m,1),'lablist',[1:kmax]','prior',[]); % for speed

	n_ini = max(n_ini,kmax*5);  % initialisation needs sufficient samples
  text = sprintf('k-means clustering, %i iterations: ',maxit);
	prwaitbar(maxit,text);
	% Initialise by performing KCENTRES on...
	if (size(init,1) == 1) & strcmp(init,'kcentres') & (m > n_ini)
		%prwarning(2,'Initializing by performing KCENTRES on subset of %d samples.', n_ini);
		b = +gendat(a,n_ini); % ... a random subset of A.
		d = +distm(b);
		assign = kcentres(d,kmax,[]);
		bb = setprior(prdataset(b,assign),0);
		w = nmc(bb); % Initial partition W and assignments ASSIGN.
	elseif (size(init,1) == 1) & strcmp(init,'kcentres')
		%prwarning(2,'Initializing by performing KCENTRES on training set.');
		d = +distm(a);    % ... the entire set A.
		assign = kcentres(d,kmax,[]);
		aa = setprior(prdataset(a,assign),0);
		aa = setlablist(aa);
		w = nmc(aa); % mapping trained on the complete dataset
	elseif (size(init,1) == 1) & strcmp(init,'rand')
		%prwarning(2,'Initializing by randomly selected objects');
		R = randperm(m);
		w = nmc(prdataset(a(R(1:kmax),:),[1:kmax]')); % mapping trained on kmax random samples
	elseif (size(init,1) == m)
		assign = renumlab(init);
		kmax = max(assign);
		%prwarning(2,'Initializing by given labels, k = %i',kmax);
		w = nmc(prdataset(a,assign));
	else
		error('Wrong initialisation supplied')
	end
	assign = labeld(a*w);
	a = prdataset(a,assign);
	a = setprior(a,0);
	%tmp_assign = zeros(m,1); % Allocate temporary array.

	% Main loop, while assignments change
	it=1; % number of iterations
	ndif = 1;

	while (it<maxit) & (ndif > 0)
		prwaitbar(maxit,it,[text num2str(it)]);
		tmp_assign = assign;     % Remember previous assignments.
		a = setnlab(a,assign);
    a = setlablist(a);       % remove empty classes
    % disp([it classsizes(a)])
		w = a * nmc;   % Re-partition the space by assigning samples to nearest mean.
		[dummy,assign] = max(+a*w,[],2);  % Re-calculate assignments.
		it = it+1; % increase the iteration counter
		ndif = sum(tmp_assign ~= assign);
	end
	prwaitbar(0);

	if it>=maxit
		prwarning(1,['No convergence reached before the maximum number of %d iterations passed. ' ...
		'The last result was returned.'], maxit);
	end

return;

author: marianux

repo: ecg-kit

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


text:
function [x, y, h, hm, xm, ym] = generate_terrain(n, mesh_size, h0, r0, rr)

% [x, y, h, hm, xm, ym] = generate_terrain(n, mesh_size, h0, r0, rr)
%
% This function generates a series of points that approximate terrain
% according to a simple algorithm and very few parameters.
%
% Inputs:
% n         - Number of iterations of algorithm to perform, akin to the
%             order of magnitudes represented by the variation. Anything
%             beyond 7 will produce detail too fine to notice, while the
%             algorithm will require ~3 times longer for every iteration.
% mesh_size - Size of the output mesh (e.g., 512 for 512-by-512)
% h0        - Initial elevation
% r0        - Initial roughness (how much terrain can vary in a step)
% rr        - Roughness roughness (how much roughness can vary in a step)
%
% For an example, see terrain_generation_introduction.m or
% html/terrain_generation_introduction.html.
%
% Outputs:
% x, y, and h    - Vectors of points comprising terrain
% xm, ym, and hm - Meshes over landscape, useful for surf(...)
%
% Tucker McClure
% Copyright 2012, The MathWorks, Inc.

    % Set some defaults if the user didn't provide inputs.
    if nargin < 1, n         = 7;            end % Iterations
    if nargin < 2, mesh_size = 513;          end % Output mesh size
    if nargin < 3, h0        = 0.1 * rand(); end % Elevation
    if nargin < 4, r0        = 0.1 * rand(); end % Roughness
    if nargin < 5, rr        = 0.1 * rand(); end % Roughness roughness

    n0 = randi(5);     % Number of initial points
    m  = 3;            % How many points grow from each old point
    nf = n0 * (m+1)^n; % Total number of points

    % Create initial x, y, and height coordinates and roughness map.
    x = [randn(n0, 1);           zeros(nf-n0, 1)];
    y = [randn(n0, 1);           zeros(nf-n0, 1)];
    h = [r0 * randn(n0, 1) + h0; zeros(nf-n0, 1)];
    r = [rr * randn(n0, 1) + r0; zeros(nf-n0, 1)];

    % Create new points from old points n times.
    for k = 1:n

        % Calculate the new variance for the x, y random draws and for the
        % h, r random draws.
        dxy = 0.75^k;
        dh  = 0.5^k;

        % Number of new points to generate
        n_new = m * n0;

        % Parents for new points
        parents = reshape(repmat(1:n0, m, 1), [n_new, 1]);

        % Calculate indices for new and existing points.
        new = (n0+1):(n0+n_new);
        old = 1:n0;

        % Generate new x/y values.
        theta  = 2*pi * rand(n_new, 1);
        radius = dxy * (rand(n_new, 1) + 1);
        x(new) = x(parents) + radius .* cos(theta);
        y(new) = y(parents) + radius .* sin(theta);

        % Interpolate to find nominal new r and h values and add noise to
        % roughness and height maps.
        r(new) =   interpolate(x(old), y(old), r(old), x(new), y(new)) ...
                 + (dh * rr) .* randn(n_new, 1);
        h(new) =   interpolate(x(old), y(old), h(old), x(new), y(new)) ...
                 + (dh/dxy) * radius .* r(new) .* randn(n_new, 1);

        % Add up the new points.
        n0 = n_new + n0;

    end

    % Normalize the distribution of the points about the median.
    x = (x - median(x))/std(x);
    y = (y - median(y))/std(y);

    % If the user wants a mesh output, we can do that too. Create a mesh
    % over the significant part and interpolate over it.
    if nargout > 3 && ~isempty(mesh_size)
        [xm, ym] = meshgrid(linspace(-1, 1, mesh_size));
        hm = interpolate(x, y, h, xm, ym);
    end

end

% Perform our particular type of interpolation.
function vn = interpolate(x0, y0, v0, xn, yn)

    % Make an interpolator for height or roughness. We'll add safe far-away
    % points so we'll always be interpolating and not extrapolating.
    int = TriScatteredInterp([100*[-1 -1 1 1]'; x0], ...
                             [100*[-1 1 -1 1]'; y0], ...
                             [zeros(4, 1);      v0], 'linear'); %#ok<R%#ok<MSNU> EMFF1>

	% Perform the actual interpolation.
	vn = int(xn, yn);

end


text:
function varargout = Fading_BER_Curves(SNRs)
%ConvCode_BER_Curves Bit Error Rate plots for convolutional coded system.
%   h = ConvCode_BER_Curves plots an upper bound of the bit error rate
%   (BER) verus SNR per information bit (Eb/No) for a constraint length 7
%   (K=7), rate 1/2 code (R=1/2) over an AWGN channel with BPSK modulation.
%
%   The theoretical results are produced using the BERCODING function
%   from the Communications Toolbox, which uses expressions taken from:
%   [1] J. G. Proakis, Digital Communications, McGraw-Hill, 4th edition, 2001.

%   Written by Idin Motedayen-Aval
%   Applications Engineer
%   The MathWorks, Inc.
%   zq=[4 2 5 -15 -1 -3 24 -57 45 -12 19 -12 15 -8 3 -7 8 -69 53 12 -2];
%   char(filter(1,[1,-1],[105 zq])), clear zq

x = SNRs;    % Eb/No range

figure1 = figure;
k = 1;
for L = [1, 2, 3, 4, 6, 8];
    ber(k,:) = berfading(x,'psk',2,L);
    k = k+1;
end
ber(k,:) = berawgn(x,'psk',2,'nondiff'); % AWGN BER curve for reference

% Plot the results
line_h = semilogy(x,ber,':');
grid on
ylim([1e-006 1]);
xlim([min(SNRs) max(SNRs)]);
% legend show

% Create title
myT = sprintf('BER, Rayleigh Flat Fading Channel with diversity');
title(myT,'Interpreter','latex');
% Create xlabel
xlabel('SNR per transmitted bit, $^{E_b}/_{N_0}$ (dB)','Interpreter','latex');
% Create ylabel
ylabel('Bit Error Rate, BER');

% Create annotations

hold off

if nargout
    varargout{1} = figure1;
    if nargout > 1
        varargout{2} = line_h;
    end
end


text:
function[H]=loadH(pref)
% load the data from Excel
[H txt]=xlsread(pref.excelFile,'Long run prior');
end

author: european-central-bank

repo: BEAR-toolbox


text:
%%  Virtual Painter 1
%% discription
% Virtual Painter - this software allows you to draw a picture on the screen
% By following the red color
%
% Have the option to choose what color to paint by moving the red color to the desired color:
% Red, green, blue.
%
% To close the program, have received faint image on the screen
%
% Note: In order to achieve optimal results have run the program in the light of fluorescent lamp.
% % Operating Instructions
% To paint the color red object we vote to act as follows:
%% Play and exit:
%
% * To start the video: click on the Play button.
% * To exit from the function: Display to the camera very dark image ,
%    for example you can cover the camera with your hand,
%    Or turn off the light in the room.After doing so
%    The Camera and function should be closed automatically.
%% Notes
%  The identification of red objects in the picture depends on the lighting of the room.
%  Lighting Fluorescent light is the most appropriate lighting to this function to identify red objects.
%
%  first edited  in 8/06/2013 Saturday .
%  by Oren berkovich.

function virtualpaint()

%% creat figure and uicontrol
      figure('menubar','none')

    ply= uicontrol('style','togg','units','normalized',...
        'position',[0.1 0.1 0.1 0.1],'str','Play','callback',{@PlayvideoLive});

    uicontrol('style','text','str' ,'exit- cover the eye of the camera ',...
        'units','normalized','position',[0.2 0.1 0.5 0.05])

    %% main function
    function PlayvideoLive(~,~)

        set(ply,'enable','off');
        pause(0.5);

        try
    imaqreset
    vid = videoinput('winvideo', 1, 'YUY2_160x120');
    vid.ReturnedColorspace = 'rgb';
    set(vid,'framesperTrigger',10,'TriggerRepeat',Inf);

      color='y';
      x=[];
      y=[];
      Eraser='non';
      %
        tic

    start(vid);
        catch ex
            err=errordlg(ex.message);
            uiwait(err)
            close gcf
            return
        end

      while islogging(vid)

    tv=getdata(vid,1);

    r=tv(:,:,1,1);
    g=tv(:,:,2,1);
    b=tv(:,:,3,1);

    flushdata(vid);

    %% Find the differences of each object

rg=r-g;
gb=g-b;
rggb=rg-gb;

%% build the color filters

filter=rggb>26;
filter=uint8(filter);
F=filter.*rggb;

%% Remove objects less than the pixel size prescribing function

 bw=imfill(F,'holes');

 bw=bwareaopen(bw,50);
 %%

 %% Identify the position of the marked
 W=regionprops(bw,'Centroid');
 if isempty(W)==0
     %axis([0 160 0 120])

 x=horzcat(x,W(1).Centroid(1));
 y=horzcat(y,W(1).Centroid(2));

  elseif length(W)==2
     hold on
     line([W(1).Centroid(1) W(2).Centroid(1)],[W(1).Centroid(2) W(2).Centroid(2)],'color','r','linew',2)

 end
  %% erse mode
      if length(x)==100
          Eraser='era';
      end

      if eq(Eraser,'era')==1
           x=x(2:length(x));
           y=y(2:length(y));
           if  isempty(x)==1
               Eraser='non';
           end
      end
      %%
 hold on

 if isempty(x)==0
color=rgbcolor(x(length(x)),y(length(y)),color);
plot(x,y,'color',color,'linew',4)
 end

 %% Rebuild the RGB image, if the part has been lightened

    figure(gcf)
    hold off
    imshow(tv(:,:,:,1))

    state =mean2(r);
    if state<20&&toc>20;
       break
    end

      end
      % after you exit the loop
       stop(vid)
       imaqreset

       close gcf
       clc
       disp('good by')
    end

    end


text:
function [c ceq gradc gradceq] = hs71C(x)

c = -(prod(x)-25);
ceq = sum(x.^2)-40;

if(nargout > 2)
    gradc = -(prod(x)./x')';
    gradceq = 2*x;
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function spm_dem_cue_movie(DEM,q)
% creates a movie of cued pointing
% FORMAT spm_dem_cue_movie(DEM,q)
%
% DEM - DEM structure from reaching simulations
% q   - flag switching from true to perceived reaching
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dem_cue_movie.m 4230 2011-03-07 20:58:38Z karl $

% Dimensions
%--------------------------------------------------------------------------
N    = size(DEM.pU.x{1},2);
n    = size(DEM.pP.P{1},2);

% evaluate true location (targets)
%----------------------------------------------------------------------
for i = 1:N
    L(:,:,i) = DEM.pP.P{1};
end

if nargin > 1

    % evaluate perceived positions (motor plant)
    %----------------------------------------------------------------------
    x    = tan(DEM.qU.x{1}(1:2,:));                  % finger location
    c    = DEM.qU.x{1}(2 + (1:n),:);                 % target contrast

else

    % evaluate true positions (motor plant)
    %----------------------------------------------------------------------
    x    = tan(DEM.pU.x{1}(1:2,:));                  % finger location
    c    = DEM.pU.x{1}(2 + (1:n),:);                 % target contrast

end
c    = c - min(c(:)) + 1/32;
c    = c/max(c(:));
fin  = imread('pointfinger.jpg');

% movie
%--------------------------------------------------------------------------
s     = 2;
for i = 1:N

    cla
    axis image ij
    hold on

    % finger
    %----------------------------------------------------------------------
    imagesc(([-1 0] + .68)*s + x(1,i),([-1 0] + .96)*s + x(2,i),fin);
    hold on

    % trajectory
    %----------------------------------------------------------------------
    plot(x(1,1:i),x(2,1:i),'k:')

    % targets
    %----------------------------------------------------------------------
    for j = 1:n
        plot(L(1,j,i),L(2,j,i),'.','MarkerSize',64,'color',[c(j,i) (1 - c(j,i)) 0])
    end
    axis([-1 1 -1 1]*2)
    hold off
    drawnow

    % save
    %----------------------------------------------------------------------
    M(i) = getframe(gca);

end

% set ButtonDownFcn
%--------------------------------------------------------------------------
h = findobj(gca,'type','image');
set(h(1),'Userdata',{M,16})
set(h(1),'ButtonDownFcn','spm_DEM_ButtonDownFcn')

author: spm


text:
function fun = sym2fun(sobj,svar,var,skipSubs)
% SYM2FUN  Convert Symbolic Expression into Matlab Function Handle

if(nargin < 3 || isempty(var)), var = 'x'; end
if(nargin < 4 || isempty(skipSubs)), skipSubs = false; end

%Build cell array to store indexed vars
ivar = convIndex(svar);

%Subs out individual symbolic variables into our indexed list
if(~skipSubs)
    wstate = warning('off','symbolic:sym:sym:DeprecateExpressions');
    eq = subs(sobj,svar,ivar);
    warning(wstate);
else
    eq = sobj;
end

%Build a MATLAB function
str = char(eq);
%Remove 'matrix' if is found
ind = strfind(str,'matrix');
if(ind)
    str = str(ind+7:end-1);
end
%Ensure we have right size
[r,c] = size(eq);
if(c==1 && r > 1)
    str = regexprep(str,',',';'); %force column
elseif(c>1 && r>1)
    str = regexprep(str,'],','];'); %keep matrix
end
%Build function handle
fun = str2func(['@(' var ') ' str]);

%Convert index from m1 to m(1)
function [ivar,var] = convIndex(svar)
ivar = cell(size(svar));
%Find indexing number
for i = 1:length(svar)
    str = char(svar(i));
    for j = 1:length(str)
        if(~isletter(str(j)))
            break;
        end
    end
    if(i == 1) %assume first variable contains decision variable
        var = str(1:j-1);
    end
%   ivar{i} = sprintf('%s(%s)',str(1:j-1),str(j:end)); %this uses user index...
    ivar{i} = sprintf('%s(%d)',str(1:j-1),i); %this uses index based on array
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function c=comp_nonsepdgt_multi(f,g,a,M,lt)
%COMP_NONSEPDGT_MULTI  Compute Non-separable Discrete Gabor transform
%   Usage:  c=comp_nonsepdgt_multi(f,g,a,M,lt);
%
%   This is a computational subroutine, do not call it directly.

%   AUTHOR : Nicki Holighaus and Peter L. Søndergaard
%   TESTING: TEST_NONSEPDGT
%   REFERENCE: REF_NONSEPDGT

% Assert correct input.

L=size(f,1);
W=size(f,2);
N=L/a;

% ----- algorithm starts here, split into sub-lattices ---------------

c=zeros(M,N,W,assert_classname(f,g));

mwin=comp_nonsepwin2multi(g,a,M,lt,L);

% simple algorithm: split into sublattices

for ii=0:lt(2)-1
    c(:,ii+1:lt(2):end,:)=comp_dgt(f,mwin(:,ii+1),lt(2)*a,M,[0 1],0,0,0);
end;

% Phase factor correction
E = zeros(1,N,assert_classname(f,g));
for win=0:lt(2)-1
    for n=0:N/lt(2)-1
        E(win+n*lt(2)+1) = exp(-2*pi*i*a*n*rem(win*lt(1),lt(2))/M);
    end;
end;

c=bsxfun(@times,c,E);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function circle_integrals_test ( )

%*****************************************************************************80
%
%% CIRCLE_INTEGRALS_TEST tests the CIRCLE_INTEGRALS library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 January 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CIRCLE_INTEGRALS_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the CIRCLE_INTEGRALS library.\n' );

  circle_integrals_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CIRCLE_INTEGRALS_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function pass = test_toFileEIG(pref)
%TEST_TOFILEEIG  Test exporting all EIG demos to an m-file.
%
% This test only checks whether nothing breaks, it does not try to solve the
% problems.

% This test won't work if we can't write to the current directory.
[ignored, attr] = fileattrib(pwd);
if ( ~attr.UserWrite )
    warning('CHEBFUN:tests:chebgui:test_toFileEIG:perms', ...
        'Cannot write to chebfunroot/tests/chebgui/.  Bypassing test.');
    pass = true;
    return;
end

% Disable discretization warnings.
% TODO: Remove this when #1555 gets merged.
warnstate = warning('OFF', 'CHEBOPPREF:PARSEDISCRETIZATION');

% Find the folders which demos are stored in. The chebguiDemos folder lives in
% the trunk folder, find the path of the Chebfun trunk.
trunkPath = chebfunroot();

% Append directory information
bvppath = fullfile(trunkPath, 'chebguiDemos', 'eigdemos');

% Set up file to exporting to
tempPath = fullfile(trunkPath, 'tests', 'chebgui');
tempFileName = 'tempExporterTest.m';

% Obtain directory information:
D = dir(bvppath);

% Create a BVP exporter object:
exporter = chebguiExporter.constructor('eig');

% Wrap the test in a try-catch statement to ensure we don't files hanging around
% in the test folder:
try

    % Loop through the demos and export
    for demoCounter = 3:length(D)

        % The path to the .guifile
        demoPath = fullfile(bvppath, D(demoCounter).name);

        % Load demo
        cg = chebgui.demo2chebgui(demoPath);

        % Export the demo!
        toFile(exporter, cg, tempFileName, tempPath)
    end

    % Got here without crashing == success!
    pass = 1;

catch
    % Something went wrong == fail :(
    pass = 0;
end

% Delete the temporary file we wrote to:
delete(fullfile(tempPath, tempFileName))

% TODO: Remove this when #1555 gets merged.
warning(warnstate.state, 'CHEBOPPREF:PARSEDISCRETIZATION')
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function y = flipudf(x)
%FLIPUD Flip matrix in up/down direction.
%   FLIPUD(X) returns X with columns preserved and rows flipped
%   in the up/down direction.  For example,
%
%   X = 1 4      becomes  3 6
%       2 5               2 5
%       3 6               1 4
%
%   Class support for input X:
%      float: double, single
%
%   See also FLIPLR, ROT90, FLIPDIM.

%   Copyright 1984-2010 The MathWorks, Inc.
%   $Revision: 5.9.4.4 $  $Date: 2010/02/25 08:08:49 $

% if ~ismatrix(x)
%   error('MATLAB:flipud:SizeX', 'X must be a 2-D matrix.');
% end
y = x(end:-1:1,:);

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function [ val ] = check_dominance(vec1, vec2)
%   This procedure checks the dominance of two vectors,
%   where the vectors are in the form [obj,cv]
%   The non-domnination check will work like this --
%       1   if vec1 dominates vec2
%      -1   if vec2 dominates vec1
%       0   if both vec1 and vec2 are non-dominated
%
%   Please note: Vectorization on this function may not be too
%   efficient as it operates on only two vectors, and couple of
%   variables.

% split the obj and cv values
obj1 = vec1(1:end-1);
cv1 = vec1(end);
obj2 = vec2(1:end-1);
cv2 = vec2(end);

if(cv1 < 0 && cv2 < 0)
    if(cv1 > cv2)
        val = 1 ;
        return ;
    elseif(cv1 < cv2)
        val = -1 ;
        return ;
    else
        val = 0;
        return ;
    end
else
    if(cv1 < 0 && cv2 == 0)
        val = -1 ;
        return ;
    elseif(cv1 == 0 && cv2 < 0)
        val = 1  ;
        return ;
    else
        if(any(obj1 < obj2) && ~(any(obj1 > obj2)))
            val = 1 ;
        elseif(~(any(obj1 < obj2)) && any(obj1 > obj2))
            val = -1 ;
        else
            val = 0 ;
            return ;
        end
    end
end
% end of check_dominance
end

author: chudur-budur


text:
%inPoints = getPolygonGrid(xv,yv,ppa) returns points that are within a
%concave or convex polygon using the inpolygon function.

%xv and yv are columns representing the vertices of the polygon, as used in
%the Matlab function inpolygon

%ppa refers to the points per unit area you would like inside the polygon.
%Here unit area refers to a 1.0 X 1.0 square in the axes.

%Example:
% L = linspace(0,2.*pi,6); xv = cos(L)';yv = sin(L)'; %from the inpolygon documentation
% inPoints = getPolygonGrid(xv, yv, 10^5)
% plot(inPoints(:, 1),inPoints(:,2), '.k');

function [inPoints] = polygrid( xv, yv, ppa)

	N = sqrt(ppa);
%Find the bounding rectangle
	lower_x = min(xv);
	higher_x = max(xv);

	lower_y = min(yv);
	higher_y = max(yv);
%Create a grid of points within the bounding rectangle
	inc_x = 1/N;
	inc_y = 1/N;

	interval_x = lower_x:inc_x:higher_x;
	interval_y = lower_y:inc_y:higher_y;
	[bigGridX, bigGridY] = meshgrid(interval_x, interval_y);

%Filter grid to get only points in polygon
	in = inpolygon(bigGridX(:), bigGridY(:), xv, yv);
%Return the co-ordinates of the points that are in the polygon
	inPoints = [bigGridX(in), bigGridY(in)];

end


text:
% this function return the transformer's (input argument) loss
function FValue=SP_tlosscost1(InValue)

global KWTransformersTypes ShortCircuitLosses
Tag = find(KWTransformersTypes == InValue);
FValue = ShortCircuitLosses(1,Tag);


text:
function q = sphere_stereograph ( p )

%*****************************************************************************80
%
%% SPHERE_STEREOGRAPH computes the stereographic image of points on a sphere.
%
%  Discussion:
%
%    We start with a sphere of radius 1 and center (0,0,0).
%
%    The north pole N = (0,0,1) is the point of tangency to the sphere
%    of a plane, and the south pole S = (0,0,-1) is the focus for the
%    stereographic projection.
%
%    For any point P on the sphere, the stereographic projection Q of the
%    point is defined by drawing the line from S through P, and computing
%    Q as the intersection of this line with the plane.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 November 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    C F Marcus,
%    The stereographic projection in vector notation,
%    Mathematics Magazine,
%    Volume 39, Number 2, March 1966, pages 100-102.
%
%  Parameters:
%
%    Input, real P(3,N), a set of points on the unit sphere.
%
%    Output, real Q(3,N), the coordinates of the image points,
%
  n = size ( p, 2 );
  q = zeros ( 3, n );

  q(1,1:n) =  2.0 * p(1,1:n) ./ ( 1.0 + p(3,1:n) );
  q(2,1:n) =  2.0 * p(2,1:n) ./ ( 1.0 + p(3,1:n) );
  q(3,1:n) =  1.0;

  return
end


text:
function x = nco_abscissas_ab ( a, b, n )

%*****************************************************************************80
%
%% NCO_ABSCISSAS_AB computes the Newton Cotes Open abscissas for [A,B].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the endpoints of the interval.
%
%    Input, integer N, the order of the rule.
%
%    Output, real X(N), the abscissas.
%
  for i = 1 : n
    x(i) = ( ( n - i + 1 ) * a   ...
           + (     i     ) * b ) ...
           / ( n     + 1 );
  end

  return
end


text:
function [DCM] = spm_dcm_ssr_results(DCM,Action)
% Results for ERP Dynamic Causal Modeling (DCM)
% FORMAT spm_dcm_erp_results(DCM,'spectral data');
% FORMAT spm_dcm_erp_results(DCM,'Coupling (A)');
% FORMAT spm_dcm_erp_results(DCM,'Coupling (B)');
% FORMAT spm_dcm_erp_results(DCM,'Coupling (C)');
% FORMAT spm_dcm_erp_results(DCM,'trial-specific effects');
% FORMAT spm_dcm_erp_results(DCM,'Input');
% FORMAT spm_dcm_erp_results(DCM,'Cross-spectral density');
% FORMAT spm_dcm_erp_results(DCM,'Dipoles');
%
%___________________________________________________________________________
%
% DCM is a causal modelling procedure for dynamical systems in which
% causality is inherent in the differential equations that specify the model.
% The basic idea is to treat the system of interest, in this case the brain,
% as an input-state-output system.  By perturbing the system with known
% inputs, measured responses are used to estimate various parameters that
% govern the evolution of brain states.  Although there are no restrictions
% on the parameterisation of the model, a bilinear approximation affords a
% simple re-parameterisation in terms of effective connectivity.  This
% effective connectivity can be latent or intrinsic or, through bilinear
% terms, model input-dependent changes in effective connectivity.  Parameter
% estimation proceeds using fairly standard approaches to system
% identification that rest upon Bayesian inference.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_ssr_results.m 4096 2010-10-22 19:40:34Z karl $

% get figure handle
%--------------------------------------------------------------------------
Fgraph = spm_figure('GetWin','Graphics');
colormap(gray)
figure(Fgraph)
clf

% placespectral features in xY.y
%--------------------------------------------------------------------------
DCM.xY.y  = spm_cond_units(DCM.xY.csd,'csd');

% trial data
%--------------------------------------------------------------------------
xY  = DCM.xY;                   % data
nt  = length(xY.y);             % Nr trial types
nf  = size(xY.y{1},1);          % Nr frequency bins
nm  = size(xY.y{1},2);          % Nr spatial modes
Hz  = xY.Hz;                    % PST

% switch
%--------------------------------------------------------------------------
switch(lower(Action))

case{lower('spectral data')}

    % spm_dcm_ssr_results(DCM,'Data');
    %----------------------------------------------------------------------
    co = {'b', 'r', 'g', 'm', 'y', 'k', 'c'};
    Hz = xY.Hz;
    q  = max(spm_vec(xY.y));
    nm = min(nm,4);

    for k = 1:nt
        str{k} = sprintf('trial %i',k);
    end

    for i = 1:nm
        for j = i:nm

            % for each trial type
            %--------------------------------------------------------------
            subplot(nm,nm,(i - 1)*nm + j),cla
            for k = 1:nt
                plot(Hz,xY.y{k}(:,i,j),'color',co{k}), hold on
                set(gca,'YLim',[0 q])
            end
        end

        % spectral density
        %------------------------------------------------------------------
        subplot(2,2,3)
        for k = 1:nt
            plot(Hz,xY.y{k}(:,i,i),'color',co{i}), hold on
            set(gca,'YLim',[0 q])
        end
    end

    title('spectral density over modes')
    xlabel('Frequency (Hz)')
    ylabel('CSD')
    axis square
    return

end

% post inversion parameters
%--------------------------------------------------------------------------
nu  = length(DCM.B);          % Nr experimental inputs
ns  = size(DCM.A{1},2);       % Nr of sources

% switch
%--------------------------------------------------------------------------
switch(lower(Action))

case{lower('Coupling (A)')}

    % spm_dcm_ssr_results(DCM,'coupling (A)');
    %----------------------------------------------------------------------
    str = {'Forward','Backward','Lateral'};
    for  i = 1:3

        % images
        %------------------------------------------------------------------
        subplot(4,3,i)
        imagesc(exp(DCM.Ep.A{i}))
        title(str{i},'FontSize',10)
        set(gca,'YTick',[1:ns],'YTickLabel',DCM.Sname,'FontSize',8)
        set(gca,'XTick',[])
        xlabel('from','FontSize',8)
        ylabel('to','FontSize',8)
        axis square

        % table
        %------------------------------------------------------------------
        subplot(4,3,i + 3)
        text(0,1/2,num2str(full(exp(DCM.Ep.A{i})),' %.2f'),'FontSize',8)
        axis off,axis square

        % PPM
        %------------------------------------------------------------------
        subplot(4,3,i + 6)
        image(64*DCM.Pp.A{i})
        set(gca,'YTick',[1:ns],'YTickLabel',DCM.Sname,'FontSize',8)
        set(gca,'XTick',[])
        title('PPM')
        axis square

        % table
        %------------------------------------------------------------------
        subplot(4,3,i + 9)
        text(0,1/2,num2str(DCM.Pp.A{i},' %.2f'),'FontSize',8)
        axis off, axis square

    end

case{lower('Coupling (C)')}

    % spm_dcm_ssr_results(DCM,'coupling (C)');
    %----------------------------------------------------------------------

    % images
    %----------------------------------------------------------------------
    subplot(2,4,1)
    imagesc(exp(DCM.Ep.C))
    title('Factors','FontSize',10)
    set(gca,'XTick',[1:nu],'XTickLabel','Input','FontSize',8)
    set(gca,'YTick',[1:ns],'YTickLabel',DCM.Sname, 'FontSize',8)
    axis square

    % PPM
    %----------------------------------------------------------------------
    subplot(2,4,3)
    image(64*DCM.Pp.C)
    title('Factors','FontSize',10)
    set(gca,'XTick',[1:nu],'XTickLabel','Input','FontSize',8)
    set(gca,'YTick',[1:ns],'YTickLabel',DCM.Sname, 'FontSize',8)
    axis square
    title('PPM')

    % table
    %----------------------------------------------------------------------
    subplot(2,4,2)
    text(0,1/2,num2str(full(exp(DCM.Ep.C)),' %.2f'),'FontSize',8)
    axis off

    % table
    %----------------------------------------------------------------------
    subplot(2,4,4)
    text(0,1/2,num2str(DCM.Pp.C,' %.2f'),'FontSize',8)
    axis off

case{lower('Coupling (B)')}

    % spm_dcm_ssr_results(DCM,'coupling (B)');
    %----------------------------------------------------------------------
    for i = 1:nu

        % images
        %------------------------------------------------------------------
        subplot(4,nu,i)
        imagesc(exp(DCM.Ep.B{i}))
        title(DCM.xU.name{i},'FontSize',10)
        set(gca,'YTick',[1:ns],'YTickLabel',DCM.Sname,'FontSize',8)
        set(gca,'XTick',[])
        xlabel('from','FontSize',8)
        ylabel('to','FontSize',8)
        axis square

        % tables
        %------------------------------------------------------------------
        subplot(4,nu,i + nu)
        text(0,1/2,num2str(full(exp(DCM.Ep.B{i})),' %.2f'),'FontSize',8)
        axis off
        axis square

        % PPM
        %------------------------------------------------------------------
        subplot(4,nu,i + 2*nu)
        image(64*DCM.Pp.B{i})
        set(gca,'YTick',[1:ns],'YTickLabel',DCM.Sname,'FontSize',8)
        set(gca,'XTick',[])
        title('PPM')
        axis square

        % tables
        %------------------------------------------------------------------
        subplot(4,nu,i + 3*nu)
        text(0,1/2,num2str(DCM.Pp.B{i},' %.2f'),'FontSize',8)
        axis off
        axis square

    end

case{lower('trial-specific effects')}

    % spm_dcm_ssr_results(DCM,'trial-specific effects');
    %----------------------------------------------------------------------
    for i = 1:ns
        for j = 1:ns

            % ensure connection is enabled
            %--------------------------------------------------------------
            q     = 0;
            for k = 1:nu
                q = q | DCM.B{k}(i,j);
            end

            % plot trial-specific effects
            %--------------------------------------------------------------
            if q
                B     = zeros(nt,1);
                for k = 1:nu
                    B = B + DCM.xU.X(:,k)*DCM.Ep.B{k}(i,j);
                end

                subplot(ns,ns,(i - 1)*ns + j)
                bar(exp(B)*100,'c')
                title([DCM.Sname{j}, ' to ' DCM.Sname{i}],'FontSize',10)
                xlabel('trial',  'FontSize',8)
                ylabel('strength (%)','FontSize',8)
                set(gca,'XLim',[0 nt + 1])
                axis square

            end
        end
    end

case{lower('Input')}

    % spectrum of innovations or noise (Gu)
    %----------------------------------------------------------------------
    try
        Gu   = exp(DCM.Ep.a)*xY.Hz.^(-1)*2;    % spectral density of (AR) input
        Gu   = Gu + exp(DCM.Ep.b);             % spectral density of IID input
    catch
        Gu   = exp(DCM.Ep.a(1))*xY.Hz.^(-1);    % spectral density of (AR) input
        Gu   = Gu + exp(DCM.Ep.a(2));           % spectral density of IID input
    end

    % plot spectral density of innovations
    % ---------------------------------------------------------------------
    subplot(2,1,1)
    plot(xY.Hz,Gu)
    xlabel('frquency (Hz)')
    title('spectrum of innovations or noise')
    axis square, grid on

case{lower('Cross-spectral density')}

    % spm_dcm_ssr_results(DCM,'Cross-spectral density');
    %----------------------------------------------------------------------
    co = {'b', 'r', 'g', 'm', 'y', 'k', 'c'};
    Hz = xY.Hz;
    q  = max(spm_vec(DCM.Hc));
    nm = min(nm,4);

    tstr = {};
    mstr = {};
    for k = 1:nt
        tstr{end + 1} = sprintf('predicted: trial %i',k);
        tstr{end + 1} = sprintf('observed: trial %i',k);
    end
    for k = 1:nm
        mstr{end + 1} = sprintf('predicted: mode %i',k);
        mstr{end + 1} = sprintf('observed: mode %i',k);
    end

    for i = 1:nm
        for j = i:nm

            % for each trial type
            %--------------------------------------------------------------
            subplot(nm,nm,(i - 1)*nm + j),cla
            for k = 1:nt
                plot(Hz,DCM.Hc{k}(:,i,j),'color',co{k}), hold on
                plot(Hz,DCM.Hc{k}(:,i,j) + DCM.Rc{k}(:,i,j),':','color',co{k})
                set(gca,'YLim',[0 q])
            end
        end

        % legend
        %------------------------------------------------------------------
        if i == nm && j == nm
            legend(tstr)
        end

        % spectral density
        %------------------------------------------------------------------
        subplot(2,2,3)
        for k = 1:nt
            plot(Hz,DCM.Hc{k}(:,i,i),'color',co{i}), hold on
            plot(Hz,DCM.Hc{k}(:,i,i) + DCM.Rc{k}(:,i,i),':','color',co{i})
            set(gca,'YLim',[0 q])
        end
    end

    title({'Spectral density over modes';'(in channel-space)'},'FontSize',16)
    xlabel('Frequency (Hz)')
    ylabel('root CSD')
    axis square
    legend(mstr)

case{lower('Dipoles')}

    % return if LFP
    % ---------------------------------------------------------------------
    if strcmp(lower(DCM.xY.modality),'lfp')
        warndlg('There are no ECDs for these LFP data')
        return
    end

    % plot dipoles
    % ---------------------------------------------------------------------
    try
        P            = DCM.Ep;
        np           = size(P.L,2)/size(P.Lpos,2);
        sdip.n_seeds = 1;
        sdip.n_dip   = np*ns;
        sdip.Mtb     = 1;
        sdip.j{1}    = full(P.L);
        sdip.j{1}    = sdip.j{1}(:);
        sdip.loc{1}  = kron(ones(1,np),full(P.Lpos));
        spm_eeg_inv_ecd_DrawDip('Init', sdip)
    end

end


text:
function h = a2h(A);

% H = a2h(ABCD)
%
% ABCD to Hybrid-H transformation
% only for 2x2 matrices

d = A(1,1)*A(2,2) - A(1,2)*A(2,1);

while abs(A(2,2)) < 1e-8
    A(2,2) = A(2,2)*(1+rand*1e-8);
end;

h(1,1) = A(1,2)/A(2,2);
h(1,2) = d/A(2,2);
h(2,1) = -1/A(2,2);
h(2,2) = A(2,1)/A(2,2);


text:
%% getInnerVoxel
% Below is a demonstration of the features of the |getInnerVoxel| function

%%
clear; close all; clc;

%% Syntax
% |[indInternal]=getInnerVoxel(L,searchRadius,plotOn);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% reads multi-channel recording file to a matrix
% function .lfp] = function readmulti(fname,numchannel,chselect)
% last argument is optional (if omitted, it will read all the
% channels

function .lfp] = readmulti_ss(fname,numchannel,chselect)

if nargin == 2
  datafile = fopen(fname,'r');
 .lfp = fread(datafile,[numchannel,inf],'int16');
  fclose(datafile);
 .lfp =.lfp';
  return
end

if nargin == 3

  % the real buffer will be buffersize * numch * 2 bytes
  % (short = 2bytes)

  buffersize = 4096;

  % get file size, and calculate the number of samples per channel
  fileinfo = dir(fname);
  N_EL = ceil(fileinfo(1).bytes / 2 / numchannel);

  datafile = fopen(fname,'r');

  mmm = sprintf('%d elements',N_EL);
%  disp(mmm);

 .lfp=zeros(length(chselect),N_EL);
  N_EL=0;
  numelm=0;
  while ~feof(datafile),
    [data,count] = fread(datafile,[numchannel,buffersize],'int16');
    numelm = count/numchannel;
    if numelm>0 % Kenji modified 061009.Otherwise if numelm == 0 an error occur.
       .lfp(:,N_EL+1:N_EL+numelm) = data(chselect,:);
        N_EL = N_EL+numelm;
    end % Kenji modified 061009.Otherwise if numelm == 0 an error occur.

end
fclose(datafile);

end

eeg =.lfp';

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function lobatto_compute_test ( )

%*****************************************************************************80
%
%% LOBATTO_COMPUTE_TEST tests LOBATTO_COMPUTE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LOBATTO_COMPUTE_TEST\n' );
  fprintf ( 1, '  LOBATTO_COMPUTE computes a Lobatto rule;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '         I      X             W\n' );

  for n = 4 : 3 : 12

    [ x, w ] = lobatto_compute ( n );

    fprintf ( 1, '\n' );
    for i = 1 : n
      fprintf ( 1, '  %8d  %12f  %12f\n', i, x(i), w(i) );
    end

  end

  return
end


text:
function [imgOut, exposure_value] = SelectOverexposedTMO(img, percent)
%
%        [imgOut, exposure_value] = SelectOverexposedTMO(img, percent)
%
%
%        A simple TMO that generates an SDR image with a given percentage
%        of overexposed pixels.
%
%        Input:
%           -img: input HDR image
%           -percent: percentage of overexposed pixels.
%
%        Output:
%           -imgOut: a single exposure with the selected percentage of
%           overexposed pixels.
%
%
%
%     Copyright (C) 2022 Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

if ~exist('percent', 'var')
    percent = 0.1;
end

gamma_inv = 1.0 / 2.2;

    function err = residual(c)

        if c > 0.0
            tmp = ClampImg((img * c).^gamma_inv, 0.0, 1.0);
            L = lum(tmp);
            tmp_percent = length(find(L >= 0.95)) / numel(L);

            err = abs(tmp_percent - percent);
        else
            err = 1e9;
        end
    end

exposure_start = 1.0;
exposure_value = fminsearch(@residual, exposure_start, getOpts);

imgOut = ClampImg(img * exposure_value, 0.0, 1.0);

end

author: banterle


text:
function p = pyrGaussGen(img)
%
%
%        p = pyrGaussGen(img)
%
%
%        Input:
%           -img: an image
%
%        Output:
%           -p: a Gaussian pyramid of img
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

[r, c, ~] = size(img);
levels = floor(log2(min(r, c)));
list = [];
for i=1:levels
    %Detail layer
    ts   = struct('detail', img);
    list = [list, ts];

    %Next level
    img = pyrGaussGenAux(img);
end

%Base layer
p = struct('list', list, 'base', img);

end


text:
%{
 * Copyright (C) 2020-2030, The Regents of The University of Michigan.
 * All rights reserved.
 * This software was developed in the Biped Lab (https://www.biped.solutions/)
 * under the direction of Jessy Grizzle, grizzle@umich.edu. This software may
 * be available under alternative licensing terms; contact the address above.
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 * The views and conclusions contained in the software and documentation are those
 * of the authors and should not be interpreted as representing official policies,
 * either expressed or implied, of the Regents of The University of Michigan.
 *
 * AUTHOR: Bruce JK Huang (bjhuang[at]umich.edu)
 * WEBSITE: https://www.brucerobot.com/
%}

function point_cloud = loadPointCloud(path, name)
    if nargin == 1
        file = path;
    else
        file = path + name;
    end

    if ~isfile(file)
        disp('Invalid path when loading dataset')
        fprintf("path: %s\n", file);
        return
    else
        pc = load(file);
        point_cloud = pc.point_cloud; % [scan, point, [X, Y, Z, I, R]]
    end
end

author: UMich-BipedLab


text:
function variance = normal_truncated_a_variance ( mu, s, a )

%*****************************************************************************80
%
%% NORMAL_TRUNCATED_A_VARIANCE: variance of the lower truncated Normal PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real MU, S, the mean and standard deviation of the
%    parent Normal distribution.
%
%    Input, real A, the lower truncation limit.
%
%    Output, real VARIANCE, the variance of the PDF.
%
  alpha = ( a - mu ) / s;
% beta = Inf;

  alpha_pdf = normal_01_pdf ( alpha );

  alpha_cdf = normal_01_cdf ( alpha );

  variance = s * s * ( 1.0 ...
    + ( alpha * alpha_pdf ) / ( 1.0 - alpha_cdf ) ...
    - ( alpha_pdf / ( 1.0 - alpha_cdf ) ) ^ 2 );

  return
end


text:
% Copyright (C) 2017-2018 Titus Cieslewski, RPG, University of Zurich,
%   Switzerland
%   You can contact the author at <titus at ifi dot uzh dot ch>
% Copyright (C) 2017-2018 Siddharth Choudhary, College of Computing,
%   Georgia Institute of Technology, Atlanta, GA, USA
% Copyright (C) 2017-2018 Davide Scaramuzza, RPG, University of Zurich,
%   Switzerland
%
% This file is part of dslam_open.
%
% dslam_open is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% dslam_open 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 dslam_open. If not, see <http://www.gnu.org/licenses/>.

function quiver = posesToQuiver(T)
n = numel(T);
quiver = zeros(n, 6);
for i = 1:n
    quiver(i, :) = poseToQuiver(T{i});
end
end

author: uzh-rpg


text:
% SYNTAX:
% dc = hmrR_OD2Conc( dod, probe, ppf )
%
% UI NAME:
% OD_to_Conc
%
% DESCRIPTION:
% Convert OD to concentrations.
%
% INPUTS:
% dod: SNIRF.data container with the Change in OD tim course
% probe: SNIRF.probe container with the source/detector geometry
% ppf: Partial path length factors for each wavelength. This is a vector of
%      factors per wavelength.  Typical value is ~6 for each
%      wavelength if the absorption change is uniform over the volume of tissue measured.
%      To approximate the partial volume effect of a small localized absorption change
%      within an adult human head, this value could be as small as 0.1. Convention is
%      becoming to set ppf=1 and to not divide by the source-detector separation such that
%      the resultant "concentration" is in units of Molar mm (or Molar cm if those are the
%      spatial units). This is becoming wide spread in the literature but there is no
%      fixed citation. Use a value of 1 to choose this option.
%
% OUTPUTS:
% dc: SNIRF.data container with the concentration data
%
% USAGE OPTIONS:
% Delta_OD_to_Conc: dc = hmrR_OD2Conc( dod, probe, ppf )
%
% PARAMETERS:
% ppf: [1.0, 1.0]
%
% PREREQUISITES:
% Intensity_to_Delta_OD: dod = hmrR_Intensity2OD( intensity )
function dc = hmrR_OD2Conc( dod, probe, ppf )

dc = DataClass().empty();

for ii=1:length(dod)
    dc(ii) = DataClass();

    Lambda = probe.GetWls();
    SrcPos = probe.GetSrcPos();
    DetPos = probe.GetDetPos();
    nWav   = length(Lambda);
    ml     = dod(ii).GetMeasList();
    y      = dod(ii).GetDataTimeSeries();

    if length(ppf) < nWav
        warning('Length of ppf does not match the number of wavelengths. Falling back to ppf=1 for all wavelengths.');
        ppf = ones(1, nWav);
    elseif length(ppf) > nWav
        d = length(ppf)-nWav;
        ppf(end-d+1:end) = [];
    end

    if ~isempty(find(ppf==1))
        ppf = ones(size(ppf));
    end

    nTpts = size(y,1);

    e = GetExtinctions(Lambda);
    e = e(:,1:2) / 10; % convert from /cm to /mm
    einv = inv( e'*e )*e';

    lst = find( ml(:,4)==1 );
    y2 = zeros(nTpts, 3*length(lst));
    for idx=1:length(lst)
        k = 3*(idx-1)+1;
        idx1 = lst(idx);
        idx2 = find( ml(:,4)>1 & ml(:,1)==ml(idx1,1) & ml(:,2)==ml(idx1,2) );
        rho = norm(SrcPos(ml(idx1,1),:)-DetPos(ml(idx1,2),:));
        if ppf(1)~=1
            y2(:,k:k+1) = ( einv * (y(:,[idx1 idx2'])./(ones(nTpts,1)*rho*ppf))' )';
        else
            y2(:,k:k+1) = ( einv * (y(:,[idx1 idx2'])./(ones(nTpts,1)))' )';
        end
        y2(:,k+2) = y2(:,k) + y2(:,k+1);
        dc(ii).AddChannelHbO(ml(idx1,1), ml(idx1,2));
        dc(ii).AddChannelHbR(ml(idx1,1), ml(idx1,2));
        dc(ii).AddChannelHbT(ml(idx1,1), ml(idx1,2));
    end
    dc(ii).SetDataTimeSeries(y2);
    dc(ii).SetTime(dod(ii).GetTime());
end

author: BUNPC