text:
function dot = permutation_puzzle ( n )

%*****************************************************************************80
%
%% PERMUTATION_PUZZLE exhibits a peculiar property of permutations.
%
%  Discussion:
%
%    Let P1 and P2 be arbitrary permutations of the integers from 1 to 10.
%
%    Pair corresponding entries of the permutations to create the set of 10
%    points (P1(1),P2(1)) through (P1(10),P2(10)).  Now subtract 5, that is,
%    N/2, from the X and Y coordinates of each point so that their average value
%    is zero.  Then apply a rotation of 45 degrees to each vector.
%
%    Now instead of looking at 10 vectors (X(I),Y(I)), each of length 2, consider
%    the data as two vectors X and Y, each of length 10.  If the instructions have
%    been followed correctly, then X are Y are perpendicular to each other in
%    10-dimensional space.
%
%    In fct, this seems to be true for all the permutations you try.
%    You can even set P1 and P2 equal.
%    You can try other values of N.
%
%    It does not work for most other rotation angles.  For example, a 30 degree
%    rotation will not work.
%
%    Once you've picked P1 and P2, you can take an arbitrary vector Q,
%    set Q1 = Q(P1) and Q2 = Q(P2) and carry out the transformation on Q1 and
%    Q2 instead of on P1 and P2, and it will still work.
%
%    The question is, can you show that all these statements are true, and
%    can you explain why?
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the vector length.
%
%    Output, real DOT, the dot product of the two vectors.
%
  p1 = randperm ( n );
  p2 = randperm ( n );

  xy = [ p1; p2 ];

  xy = xy - ( n / 2 );

  angle = pi / 4;

  A = [ cos ( angle ), - sin ( angle ); ...
        sin ( angle ),   cos ( angle ) ];

  xy = A * xy;

  dot = xy(1,1:n) * xy(2,1:n)';

  return
end


text:
% Copyright 2019 Jonas Koenemann, Moritz Diehl, University of Freiburg
% Redistribution is permitted under the 3-Clause BSD License terms. Please
% ensure the above copyright notice is visible in any derived work.
%

function [nFails] = run()
  nFails = ocl.tests.run.core();

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function x = subsasgn( x, S, y )
error( 'Subscripted assignment not allowed for these objects.' );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function out = SY_LocalDistributions(y,numSegs,eachOrPar,numPoints)
% SY_LocalDistributions  Compares the distribution in consecutive time-series segments
%
% Returns the sum of differences between each kernel-smoothed distributions
%
%---INPUTS:
%
% y, the input time series
%
% numSegs, the number of segments to break the time series into
%
% eachOrPar, (i) 'par': compares each local distribution to the parent (full time
%                       series) distribution
%            (ii) 'each': compare each local distribution to all other local
%                         distributions
%
% numPoints, number of points to compute the distribution across (in each local
%          segments) [200 by default]
%
% The operation behaves in one of two modes: each compares the distribution in
% each segment to that in every other segment, and par compares each
% distribution to the so-called 'parent' distribution, that of the full signal.
%
%---OUTPUTS: measures of the sum of absolute deviations between distributions
% across the different pairwise comparisons.

% ------------------------------------------------------------------------------
% 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/>.
% ------------------------------------------------------------------------------

% Plot outputs?
doPlot = false;

% ------------------------------------------------------------------------------
% Check inputs:
% ------------------------------------------------------------------------------
if nargin < 2 || isempty(numSegs) % number of segments
    numSegs = 5;
end
if nargin < 3 || isempty(eachOrPar)
    eachOrPar = 'par'; % compare each subsection to full (parent) distribution
end
if nargin < 4 || isempty(numPoints)
    % number of points to compute the distribution across
    numPoints = 200; % 200 by default
end

% ------------------------------------------------------------------------------
% Preliminaries
% ------------------------------------------------------------------------------
N = length(y); % Length of the time series (number of samples)
lseg = floor(N/numSegs);
dns = zeros(numPoints,numSegs);
r = linspace(min(y),max(y),numPoints); % Make range of ksdensity uniform across all subsegments

% ------------------------------------------------------------------------------
% Compute the kernel-smoothed distribution in all numSegs segments of the time series
% ------------------------------------------------------------------------------
for i = 1:numSegs
    dns(:,i) = ksdensity(y((i-1)*lseg+1:i*lseg),r,'function','pdf');
end

if doPlot
    figure('color','w')
    plot(dns,'k')
end

% ------------------------------------------------------------------------------
% Compare the local distributions
% ------------------------------------------------------------------------------
switch eachOrPar
    case {'par','parent'}
        % Compares each subdistribtuion to the parent (full signal) distribution
        pardn = ksdensity(y,r,'function','pdf');
        divs = zeros(numSegs,1);
        for i = 1:numSegs
            divs(i) = sum(abs(dns(:,i)-pardn')); % each is just divergence to parent
        end
        if doPlot
            hold on; plot(pardn,'r','LineWidth',2); hold off
        end
    case 'each'
        % Compares each subdistribtuion to the parent (full signal) distribution
        if numSegs == 2 % output is just an integer: only two distributions to compare
            out = sum(abs(dns(:,1)-dns(:,2)));
            return
        end

        % numSegs > 2: need to compare a number of different distributions against each other
        diffmat = NaN * ones(numSegs); % store pairwise differences
                                    % start as NaN to easily get upper triangle later
        for i = 1:numSegs
            for j = 1:numSegs
                if j > i
                    diffmat(i,j) = sum(abs(dns(:,i)-dns(:,j))); % store sum of absolute differences
                end
            end
        end

        divs = diffmat(~isnan(diffmat)); % (the upper triangle of diffmat)
                                         % set of divergences in all pairs of segments of the time series
        % divs = diffmat(diffmat > 0); % a set of non-zero divergences in all pairs of segments of the time series
        % if isempty(divs);
        %     fprintf(1,'That''s strange -- no changes in distribution??! This must be a really strange time series.\n');
        %     out = NaN; return
        % end
    otherwise
        error('Unknown method ''%s'', should be ''each'' or ''par''',eachOrPar);
end

%-------------------------------------------------------------------------------
% Return basic statistics on differences in distributions in different
% segments of the time series
out.meandiv = mean(divs);
out.mediandiv = median(divs);
out.mindiv = min(divs);
out.maxdiv = max(divs);
out.stddiv = std(divs);

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [rq,rqSS,diagnostics] = realized_quantile_variance(price,time,timeType,samplingType,samplingInterval,quantiles,blockSize,symmetric,overlap,subsamples)
% Computes the Quantile Realized Variance of Christensen, Oomen and Podolskij (2008), the MinRV and
% MedRV estimator of Andersen, Dobrev and Shaumberg and the general symmetrized version suggested by
% Sheppard in a discussion of COP
%
% USAGE:
%   [RQ] = realized_quantile_variance(PRICE,TIME,TIMETYPE,SAMPLINGTYPE,SAMPLINGINTERVAL,QUANTILES,BLOCKSIZE)
%   [RQ,RQSS,DIAGNOSTICS] = realized_quantile_variance(PRICE,TIME,TIMETYPE,SAMPLINGTYPE,SAMPLINGINTERVAL,QUANTILES,BLOCKSIZE,SYMMETRIC,OVERLAP,SUBSAMPLES)
%
% INPUTS:
%   PRICE            - m by 1 vector of high frequency prices
%   TIME             - m by 1 vector of times where TIME(i) corresponds to PRICE(i)
%   TIMETYPE         - String describing the way times are measured
%                       'wall'    24-hour clock of the form HHMMSS, e.g. 101543 or 153217
%                       'seconds' Time measured in seconds past midnight on the first day.
%                       'unit'  Unit normalized date format, e.g. .1, .234, .9
%                         Unit normalized times are more general than the other types and can be
%                         applied to data from more than one calendar day
%   SAMPLINGTYPE     - String describing the type of sampling to use when
%                        filtering PRICE
%                        'CalendarTime' - Sample in calendar time using observations separated by
%                          SAMPLINGINTERVAL seconds
%                        'CalendarUniform' - Sample in calendar time using SAMPLINGINTERVAL
%                          observations spread uniformly between TIME(1) and TIME(m)
%                        'BusinessTime' - Sample in business (tick) time using observation separated
%                          by SAMPLINGINTERVALticks
%                        'BusinessUniform' - Sample in business (tick) time using observations
%                          uniformly spaced in business time.
%                        'Fixed' - Sample at specific points in time. When using fixed,
%                          SAMPLINGINTERVAL must be a n by 1 vector of times with the same TIMETYPE
%                          as TIME (i.e. seconds if TIME is in seconds)
%   SAMPLINGINTERVAL - Scalar integer or n by 1 vector whose meaning depends on the selected SAMPLINGTYPE
%   QUANTILES        - k by 1 vector of quantile values to use when computing RQ.
%                        When SYMMETRIC = false, 0.5<QUANTILES<=1.
%                        When SYMMETRIC = true, 0<QUANTILES<=1
%                        In either case, QUANTILES * BLOCKSIZE must be an integer.
%                        The simplest method to set the quantiles is to determine which of the
%                        ordered returns should be used in the QRV, and then to set
%                        QUANTILES = ORDER/BLOCKSIZE.  For example, when SYMETRIC = FALSE, and
%                        BLOCKSIZE = 20, QUANTILES = [13 15 18]/20 will use the returns in position
%                        13, 15 and 18 (as well as 8 6 and 3).  When SYMMETRIC = TRUE, QUANTILES =
%                        [5 8 10 13 15 18]/20 will use the absolute value of returns in positions 5,
%                        8, 10, 13, 15 and 18.  For a fixed block size, you should generally more
%                        quantiles when SYMMETRIC = TRUE.
%   BLOCKSIZE        - Number of returns to use in each block when computing the quantiles. NOTE: If
%                        OVERLAP = FALSE, then the number of returns produced by filtering according
%                        to SAMPLINGTYPE and SAMPLINGINTERVAL must be an integer multiple of
%                        BLOCKSIZE. Otherwise there is no constraint.
%   SYMMETRIC        - [OPTIONAL] Boolean indicating whether the base the estimator on the absolute value
%                        of returns (symmetric) of on the non-adjusted returns. Default value is
%                        FALSE (as used by COP).
%   OVERLAP          - [OPTIONAL] Boolean indicating whether to use all overlapping blocks (TRUE),
%                        or only non-overlapping blocks. Default value is TRUE.
%   SUBSAMPLES       - [OPTIONAL] Scalar integer indicating the number of subsample realized
%                        variance estimators to average with the original realized variance.
%                        Subsample realized variances are based on prices uniformly spaced between
%                        the times (Calendar sampling) or ticks (Business sampling).  SUBSAMPLES=1
%                        will compute a subsample realized variance using the mid-point of the
%                        price sample points, 2 will use 1/3 and 2/3, and so on. In general this
%                        number should be small so the subsample estimators will be "sparse".
%
% OUTPUTS:
%   RQ               - Quantile realized variance vector (k by 1) corresponding to the values in QUANTILES
%   RQSS             - Subsample based version of RQ.  RQSS = RQ if SUBSAMPLES = 0 or omitted.
%   DIAGNOSTICS      - Structure with fields:
%                        RQINDIV   - k by 1 vector of RQ estimates
%                        RQSSINDIV - k by 1 vector of subsample RQ estimates
%                        RQWEIGHTS - Weights used in combining the individual RQ estimates.  RQ
%                                      weights are computed from RQCOV as the minimum variance
%                                      combination
%                        RQCOV     - k by k non-scaled (integrated quarticity) covariance matrix.
%
% COMMENTS:
%   The cases where SAMPLESPERBIN is in {4 5 6 10 13 15 18 20 25 26 30 36 39 50 60 65 72 75 100 144}
%   (non-symmetric) or {2 3 4 5 6 8 9 10 12 13 15 18 20 24 25 26 30 36 39 40 50 60 65 72 75 100 144}
%   (symmetric) the scales and non-scaled covariances have been pre-computed using Monte Carlo
%   integration using 100,000,000 simulations. Choosing another number for SAMPLESPERBIN requires a
%   run of realized_quantile_variance_scale which can be slow if SAMPLESPERBIN is large. If using
%   another value, 1,000,000 simulations will be used in computing the weights. If repeatedly using
%   another value, pre-computing the appropriate scales and non-scaled covariance using
%   realzed_quantile_weight_simulation is recommended.
%
% EXAMPLES:
%   % Estimate the RQ from prices available from 9:30 to 16:00 using 1 minute returns, 30 samples per
%   % bin and quantiles [18 22 27]/30 ([0.6 0.733 0.9])
%   sampleTimes = seconds2wall(wall2seconds(93000):60:wall2seconds(160000))
%   RQ=realized_quantile_variance(price,time,'wall','Fixed',sampleTimes,[18 22 27]/30,30);
%
%   % The same but using 4 subsamples
%   RQ=realized_quantile_variance(price,time,'wall','Fixed',sampleTimes,[18 22 27]/30,30,[],[],4);
%
%   % Estimate the RQ from prices available from 9:30 to 16:00 using 1 minute returns, 30 samples per
%   % bin, the symmetric verion and quantiles [5 10 15 18 22 27]/30 ([0.15 0.33 0.5 0.6 0.733 0.9])
%   RQ=realized_quantile_variance(price,time,'wall','Fixed',sampleTimes,[5 10 15 18 22 27]/30,30,true);
%
%   % Estimate the MedRV of ABS
%   RQ=realized_quantile_variance(price,time,'wall','Fixed',sampleTimes,2/3,3,true);
%
%   % Estimate the MinRV of ABS
%   RQ=realized_quantile_variance(price,time,'wall','Fixed',sampleTimes,1/2,2,true);
%
%   % Estimate a MedRV-like RQV only using blocks of size 5
%   RQ=realized_quantile_variance(price,time,'wall','Fixed',sampleTimes,3/5,5,true);
%
%  See also REALIZED_MIN_MED_VARIANCE, REALIZED_KERNEL, REALIZED_VARIANCE,
%  REALIZED_BIPOWER_VARIATION, REALIZED_THRESHOLD_VARIANCE

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 5/1/2008

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% InputChecking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<7 || nargin>10
    error('Seven to ten inputs required.')
end
switch nargin
    case 7
        symmetric = false;
        overlap = true;
        subsamples = 1;
    case 8
        overlap = true;
        subsamples = 1;
    case 9
        subsamples = 1;
end

if size(price,2)>size(price,1)
    price=price';
end
if size(price,2)>1
    error('PRICE must be a m by 1 vector.')
end
if size(time,2)>size(time,1)
    time=time';
end
if any(diff(time)<0)
    error('TIME must be sorted and increasing')
end
if size(time,2)>1 || length(time)~=length(price)
    error('TIME must be a m by 1 vector.')
end
% Inserted to protect against inputing integer times
time = double(time);

timeType=lower(timeType);
if ~ismember(timeType,{'wall','seconds','unit'})
    error('TIMETYPE must be one of ''wall'', ''seconds'' or ''unit''.');
end
samplingType=lower(samplingType);
if ~ismember(samplingType,{'calendartime','calendaruniform','businesstime','businessuniform','fixed'})
    error('SAMPLINGTYPE must be one of ''CalendarTime'', ''CalendarUniform'', ''BusinessTime'', ''BusinessUniform'' or ''Fixed''.');
end

m=size(price,1);
t0=time(1);
tT=time(m);
if ismember(samplingType,{'calendartime','calendaruniform','businesstime','businessuniform'})
    % Must be a scalar integer
    if ~isscalar(samplingInterval) || floor(samplingInterval)~=samplingInterval || samplingInterval<1
        error('SAMPLINGINTERVAL must be a positive integer for the SAMPLINGTYPE selected.')
    end
else
    if size(samplingInterval,2)>size(samplingInterval,1)
        samplingInterval=samplingInterval';
    end
    if ~(any(samplingInterval>=t0) && any(samplingInterval<=tT))
        error('At least one sampling interval must be between min(TIME) and max(TIME) when using ''Fixed'' as SAMPLINGTYPE.')
    end
    if any(diff(samplingInterval)<=0)
        error('When using ''Fixed'' as SAMPLINGTYPE the vector of sampling times in SAMPLINGINTERVAL must be sorted and strictly increasing.')
    end
end

% SAMPLERPERBIN, positive integer >=2
if blockSize<1 || floor(blockSize)~=blockSize
    error('SAMPLESPERBIN must be a positive integer greater than or equal to 2 (and usually at least 20)')
end

% SYMMETRIC
if isempty(symmetric)
    symmetric = false;
end
if ~isscalar(symmetric)
    error('SYMMETRIC must be a scalar logical value.')
end
symmetric = logical(symmetric);

% QUANTILES
if symmetric
    if any(quantiles<=0) || any(quantiles>1) || length(quantiles)>blockSize || length(unique(quantiles*blockSize)) ~= length(quantiles)
        error('QUANTILES must be unique, greater than 0 and less than or equal to 1, and the number of quantiles must be smaller than SAMPLESPERBIN when SYMMETRIC = true.')
    end
else
    if any(quantiles<=0.5) || any(quantiles>1) || length(quantiles)>(blockSize/2) || length(unique(quantiles*blockSize)) ~= length(quantiles)
        error('QUANTILES must be unique, greater than 0.5 and less than or equal to 1, and the number of quantiles must be smaller than 0.5*SAMPLESPERBIN when SYMMETRIC = false.')
    end
end

% OVERLAP
if isempty(overlap)
    overlap = true;
end
if ~isscalar(overlap)
    error('OVERLAP must be a scalar logical value.')
end
overlap = logical(overlap);

% SUBSAMPLES
if isempty(subsamples)
    subsamples = 1;
end
if ~isscalar(subsamples) || subsamples<0 || floor(subsamples)~=subsamples
    error('SUBSAMPLES must be a non-negative integer.')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 1. Compute returns
logPrice = log(price);
filteredLogPrice = realized_price_filter(logPrice,time,timeType,samplingType,samplingInterval);
returns = diff(filteredLogPrice);
% Check that the number of returns is compatible with SAMPLESPERBIN
n=length(returns);
if ~overlap
    if floor(n/blockSize)~=(n/blockSize)
        error(['The number of returns computed from the prices returned from realized_price_filter(PRICE,TIME,TIMETYPE,SAMPLINGTYPE,SAMPLINGINTERVAL) must be an integer multiple of SAMPLESPERBIN.  The number of return produced is ' num2str(n)]);
    end
end

if overlap
    binStart = 1:n-blockSize+1;
else
    binStart = 1:blockSize:n;
end
returns = returns * sqrt(n);
if symmetric
    returns = abs(returns);
end
q = length(quantiles);
rqindiv = zeros(length(binStart),q);

% Compute the indices to use
if symmetric
    indices = round(quantiles*blockSize);
else
    upperIndices = round(quantiles*blockSize);
    lowerIndices = round((1-quantiles)*blockSize+1);
end
% Loop over the bins

count = 1;
for j=binStart
    binreturns = returns(j:j+blockSize-1);
    binreturns = sort(binreturns)';
    if symmetric
        rqindiv(count,:)=binreturns(indices).^2;
    else
        rqindiv(count,:)=binreturns(lowerIndices).^2+binreturns(upperIndices).^2;
    end
    count = count + 1;
end
rqindiv = mean(rqindiv);

[rq,rqindiv,rqcov,rqweights] = realized_quantile_variance_core(rqindiv,blockSize,quantiles,symmetric);

diagnostics.rqindiv = rqindiv;
diagnostics.rqcov = rqcov;
diagnostics.rqweights = rqweights;

subsampledLogPrice = realized_subsample(logPrice,time,timeType,samplingType,samplingInterval,subsamples);
rqindiv = nan(subsamples * length(binStart),q);
count = 1;
for i=1:subsamples
    filteredLogPrice = subsampledLogPrice{i};
    returns = diff(filteredLogPrice);
    n = length(returns);
    returns = returns * sqrt(n);
    if symmetric
        returns = abs(returns);
    end
    binStart = 1:n-blockSize+1;
    for j=binStart
        binreturns = returns(j:j+blockSize-1);
        binreturns = sort(binreturns)';
         if symmetric
            rqindiv(count,:)=binreturns(indices).^2;
         else
            rqindiv(count,:)=binreturns(lowerIndices).^2+binreturns(upperIndices).^2;
        end
        count = count + 1;
    end
end
rqindiv = rqindiv(1:count-1,:);
rqindiv = mean(rqindiv);

[rqSS,rqindivSS] = realized_quantile_variance_core(rqindiv,blockSize,quantiles,symmetric);
diagnostics.rqindivSS = rqindivSS;

end

function [rq,rqindiv,rqcov,rqweights] = realized_quantile_variance_core(rqindiv,blockSize,quantiles,symmetric)

scales = load('realized_quantile_scales');
if symmetric
    simulationBlockSize = scales.symmetricSimulationSamplerperbin;
    simulationCovariance = scales.symmetricExpectedCovariance;
    simulationExpectedQuantiles = scales.symmetricExpectedQuantiles;
    simulationQuantile = scales.symmetricSimulationQuantile;
else
    simulationBlockSize = scales.asymmetricSimulationSamplerperbin;
    simulationCovariance = scales.asymmetricExpectedCovariance;
    simulationExpectedQuantiles = scales.asymmetricExpectedQuantiles;
    simulationQuantile = scales.asymmetricSimulationQuantile;

end
simulationBlockSize = cell2mat(simulationBlockSize);
if ismember(blockSize,simulationBlockSize )
    % If available load and use the pre-computed value

    [~,pl]=ismember(blockSize,simulationBlockSize);
    thisScale = simulationExpectedQuantiles{pl};
    thisQuantile = simulationQuantile{pl};
    thisCovariance = simulationCovariance{pl};
    q = length(quantiles);
    indicator = zeros(1,q);
    for i = 1:q
        [~,indicator(i)] = min(abs(thisQuantile-quantiles(i)));
    end
    rqExpectedSquaredQuantileValue = thisScale(indicator);
    % Find index
    rqcov     = thisCovariance(indicator,indicator);
    rqCovInv = rqcov \ eye(q);
    rqweights = rqCovInv*ones(q,1)/(ones(1,q)*rqCovInv*ones(q,1));
else
    % Need to simulate using 1,000,000 simulations
    if blockSize>100
        warning('oxfordRealized:realizedQuantileVariance','Computing the scales needed.  Since SAMPLESPERBIN is very large this may take a long time.  \nConsider pre-computing this value, especially if using this value of SAMPLESPERBIN many times.')
    else
        warning('oxfordRealized:realizedQuantileVariance','Computing the scales needed.  \nConsider pre-computing this value, especially if using this value of SAMPLESPERBIN many times.')
    end
    [rqExpectedSquaredQuantileValue,rqweights,rqcov]=realized_quantile_variance_scale(adjSamplesPerBin,quantiles,1000000);
end

rqindiv  = rqindiv ./ rqExpectedSquaredQuantileValue;
rq = rqweights'*rqindiv';
end

author: bashtage

repo: mfe-toolbox


text:
function [segments, nseg] = detectVoiced(x,fs,t)

%
% function [segments, fs] = detectVoiced(wavFileName)
%
% Theodoros Giannakopoulos
% http://www.di.uoa.gr/~tyiannak
%
% (c) 2010
%
% This function implements a simple voice detector. The algorithm is
% described in more detail, in the readme.pdf file
%
% ARGUMENTS:
%  - wavFileName: the path of the wav file to be analyzed
%  - t: if provided, the detected voiced segments are played and some
%  intermediate results are also ploted
%
% RETURNS:
%  - segments: a cell array of M elements. M is the total number of
%  detected segments. Each element of the cell array is a vector of audio
%  samples of the respective segment.
%  - fs: the sampling frequency of the audio signal
%  - nseg: number of voiced segments
%
% EXECUTION EXAMPLE:
%
% [segments, fs, nseg] = detectVoiced(audioread('example.wav'),8000,1);
%

% Convert mono to stereo
if (size(x, 2)==2)
	x = mean(x')';
end

% Window length and step (in seconds):
win = 0.050; % originally 0.050
step = 0.020; % originally 0.050

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  THRESHOLD ESTIMATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%

Weight = 5; % used in the threshold estimation method, originally 5

% Compute short-time energy and spectral centroid of the signal:
Eor = silence_removal.ShortTimeEnergy(x, win*fs, step*fs);
Cor = silence_removal.SpectralCentroid(x, win*fs, step*fs, fs);

% Apply median filtering in the feature sequences (twice), using 5 windows:
% (i.e., 250 mseconds)
NUM_WINDOWS_STE = 5;
NUM_WINDOWS_SC = 5;
E = medfilt1(Eor, NUM_WINDOWS_STE); E = medfilt1(E, NUM_WINDOWS_STE);
C = medfilt1(Cor, NUM_WINDOWS_SC); C = medfilt1(C, NUM_WINDOWS_SC);

% Get the average values of the smoothed feature sequences:
E_mean = mean(E);
Z_mean = mean(C);

% Find energy threshold:
[HistE, X_E] = hist(E, round(length(E) / 10));  % histogram computation
[MaximaE, countMaximaE] = silence_removal.findMaxima(HistE, 3); % find the local maxima of the histogram
if (size(MaximaE,2)>=2) % if at least two local maxima have been found in the histogram:
    T_E = (Weight*X_E(MaximaE(1,1))+X_E(MaximaE(1,2))) / (Weight+1); % ... then compute the threshold as the weighted average between the two first histogram's local maxima.
else
    T_E = E_mean / 1.5;
end

% Find spectral centroid threshold:
[HistC, X_C] = hist(C, round(length(C) / 10));
[MaximaC, countMaximaC] = silence_removal.findMaxima(HistC, 3);
if (size(MaximaC,2)>=2)
    T_C = (Weight*X_C(MaximaC(1,1))+X_C(MaximaC(1,2))) / (Weight+1);
else
    T_C = Z_mean / 2;
end

% Thresholding:
Flags1 = (E>=T_E);
Flags2 = (C>=T_C);
% Flags3 = (E == max(E));
flags = Flags1 & Flags2; % & Flags3;

VERBOSE = nargin==3 && logical(t) == true;

if VERBOSE % plot results:
	clf;
	subplot(3,1,1); plot(Eor, 'g'); hold on; plot(E, 'c'); legend({'Short time energy (original)', 'Short time energy (filtered)'});
    L = line([0 length(E)],[T_E T_E]); set(L,'Color',[0 0 0]); set(L, 'LineWidth', 2);
    axis([0 length(Eor) min(Eor) max(Eor)]);

    subplot(3,1,2); plot(Cor, 'g'); hold on; plot(C, 'c'); legend({'Spectral Centroid (original)', 'Spectral Centroid (filtered)'});
	L = line([0 length(C)],[T_C T_C]); set(L,'Color',[0 0 0]); set(L, 'LineWidth', 2);
    axis([0 length(Cor) min(Cor) max(Cor)]);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  SPEECH SEGMENTS DETECTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
count = 1;
WIN = 5;
Limits = [];
while (count < length(flags)) % while there are windows to be processed:
	% initilize:
	countTemp = 1;
	% while flags=1:
    while ((flags(count)==1) && (count < length(flags)))
		if (countTemp==1) % if this is the first of the current speech segment:
            Limit1 = round((count-WIN)*step*fs)+1; % set start limit:
            if Limit1 < 1
                Limit1 = 1;
            end
		end
		count = count + 1; 		% increase overall counter
		countTemp = countTemp + 1;	% increase counter of the CURRENT speech segment
    end

    if countTemp > 1 % if at least one segment has been found in the current loop:
        Limit2 = round((count+WIN)*step*fs);			% set end counter
        if (Limit2>length(x))
            Limit2 = length(x);
        end
        Limits(end+1, 1) = Limit1;
        Limits(end,   2) = Limit2;
    end
    count = count + 1; % increase overall counter
end

%%%%%%%%%%%%%%%%%%%%%%%
% POST - PROCESS      %
%%%%%%%%%%%%%%%%%%%%%%%

% A. MERGE OVERLAPPING SEGMENTS:
RUN = 1;
while (RUN==1)
    RUN = 0;
    for i=1:size(Limits,1)-1 % for each segment
        if (Limits(i,2)>=Limits(i+1,1))
            RUN = 1;
            Limits(i,2) = Limits(i+1,2);
            Limits(i+1,:) = [];
            break;
        end
    end
end

% B. Get final segments:
segments = {};
for i=1:size(Limits,1)
    segments{end+1} = x(Limits(i,1):Limits(i,2));
end

nseg = length(segments);

if VERBOSE
    subplot(3,1,3);
    % Plot results and play segments:
    time = 0:1/fs:(length(x)-1) / fs;
    for i=1:length(segments)
        hold off;
        P1 = plot(time, x); set(P1, 'Color', [0.7 0.7 0.7]);
        hold on;
        for j=1:nseg
            if (i~=j)
                timeTemp = Limits(j,1)/fs:1/fs:Limits(j,2)/fs;
                P = plot(timeTemp, segments{j});
                set(P, 'Color', [0.4 0.1 0.1]);
            end
        end
        timeTemp = Limits(i,1)/fs:1/fs:Limits(i,2)/fs;
        P = plot(timeTemp, segments{i});
        set(P, 'Color', [0.9 0.0 0.0]);
        axis([0 time(end) min(x) max(x)]);
        sound(segments{i}, fs);
%         clc;
        fprintf('Playing segment %d of %d. Press any key to continue...', i, length(segments));
        pause
    end
    clc
    hold off;
    P1 = plot(time, x); set(P1, 'Color', [0.7 0.7 0.7]);
    hold on;
    for i=1:nseg
        for j=1:nseg
            if (i~=j)
                timeTemp = Limits(j,1)/fs:1/fs:Limits(j,2)/fs;
                P = plot(timeTemp, segments{j});
                set(P, 'Color', [0.4 0.1 0.1]);
            end
        end
        axis([0 time(end) min(x) max(x)]);
    end
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function cS = plus(cS,v)

% crystal shape should be first argument
if ~isa(cS,'crystalShape'), [cS,v] = deal(v,cS); end

% shift should be vector3d
if isa(v,'double')
  if size(v,2)==2, v = [v,ones(size(v,1),1)]; end
  v = vector3d(v.');
else
  v = v(:).';
end

% extent shift by number of vertices of the crystal
v = repmat(v,size(cS.V,1),1);

% extend the crystals by the number of shifts
if size(cS.V,2) == 1
  cS = repmat(cS,size(v,2));
end

% shift vertices
cS.V = cS.V + v;

author: mtex-toolbox

repo: mtex

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


text:
function spm_opt_bfun

d       = round(get(gca,'CurrentPoint'));
i       = d(1);
pE      = get(gcf,'UserData');
P       = spm_fieldindices(pE,i);
j       = spm_fieldindices(pE,P);
j       = find(j == i);
[i j k] = ind2sub(size(spm_cat(getfield(pE,P))),j);

fprintf('\n%s(%i,%i,%i)\n',P,i,j,k);

return

author: spm


text:
function CPDpot = convert_dbn_CPDs_to_tables1(bnet, evidence)
% CONVERT_DBN_CPDS_TO_TABLES Convert CPDs of (possibly instantiated) DBN nodes to tables
% CPDpot = convert_dbn_CPDs_to_tables(bnet, evidence)
%
% CPDpot{n,t} is a table containing P(n,t|pa(n,t), ev)
% All hidden nodes are assumed to be discrete
% We assume the observed nodes are the same in every slice
%
% Evaluating the conditional likelihood of the evidence can be very slow,
% so we take pains to vectorize where possible, i.e., we try to avoid
% calling convert_to_table

[ss T] = size(evidence);
%obs_bitv = ~isemptycell(evidence(:));
obs_bitv = zeros(1, 2*ss);
obs_bitv(bnet.observed) = 1;
obs_bitv(bnet.observed+ss) = 1;

ns = bnet.node_sizes(:);
CPDpot = cell(ss,T);

for n=1:ss
  % slice 1
  t = 1;
  ps = parents(bnet.dag, n);
  e = bnet.equiv_class(n, 1);
  if ~any(obs_bitv(ps))
    CPDpot{n,t} = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, evidence{n,t});
  else
    CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps n], evidence(:,1));
  end

  % slices 2..T
  debug = 1;
  if ~obs_bitv(n)
    CPDpot = helper_hidden_child(bnet, evidence, n, CPDpot, obs_bitv, debug);
  else
    CPDpot = helper_obs_child(bnet, evidence, n, CPDpot, obs_bitv, debug);
  end
end

if 0
CPDpot2 = convert_dbn_CPDs_to_tables_slow(bnet, evidence);
for t=1:T
  for n=1:ss
    if ~approxeq(CPDpot{n,t}, CPDpot2{n,t})
      fprintf('CPDpot n=%d, t=%d\n',n,t);
      keyboard
    end
  end
end
end

% special cases: c=child, p=parents, d=discrete, h=hidden, 1=1slice
% if c=h=1 then c=d=1, since hidden nodes must be discrete
% c=h c=d p=h p=d p=1 method
% ---------------------------
% 1   1   1   1   -   replicate CPT
% 0   1   1   1   1   dhmm
% 0   0   1   1   1   ghmm
% -   1   -   1   -   evaluate CPT on evidence
% other               loop

%%%%%%%
function CPDpot = helper_hidden_child(bnet, evidence, n, CPDpot, obs_bitv, debug)

[ss T] = size(evidence);
self = n+ss;
ps = parents(bnet.dag, self);
e = bnet.equiv_class(n, 2);
ns = bnet.node_sizes(:);
if ~any(obs_bitv(ps)) % all parents are hidden (hence discrete)
  if debug, fprintf('node %d is hidden, all ps are hidden\n', n); end
  if myismember(n, bnet.dnodes)
    %CPT = CPD_to_CPT(bnet.CPD{e});
    %CPT = reshape(CPT, [prod(ns(ps)) ns(self)]);
    CPT = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, []);
    CPDpot(n,2:T) = num2cell(repmat(CPT, [1 1 T-1]), [1 2]);
  else
    error(['hidden cts node disallowed'])
  end
else % some parents are observed - slow
  if mysubset(ps, bnet.dnodes) % all parents are discrete
    % given CPT(p1, p2, p3, p4, c), where p1,p3 are observed
    % we create CPT([p2 p4 c], [p1 p3]).
    % We then convert all observed p1,p3 into indices ndx
    % and return CPT(:, ndx)
    CPT = CPD_to_CPT(bnet.CPD{e});
    domain = [ps self];
    % if dom is [3 7 8] and 3,8 are observed, odom_rel = [1 3], hdom_rel = 2,
    % odom = [3 8], hdom = 7
    odom_rel = find(obs_bitv(domain));
    hdom_rel = find(~obs_bitv(domain));
    odom = domain(odom_rel);
    hdom = domain(hdom_rel);
    CPT = permute(CPT, [hdom_rel odom_rel]);
    CPT = reshape(CPT, prod(ns(hdom)), prod(ns(odom)));
    parents_in_same_slice = all(ps > ss);
    if parents_in_same_slice
      if debug, fprintf('node %d is hidden, some ps are obs, all ps discrete, 1 slice\n', n); end
      data = cell2num(evidence(odom-ss,2:T)); %data(i,t) = val of i'th obs parent at t+1
    else
      if debug, fprintf('node %d is hidden, some ps are obs, all ps discrete, 2 slice\n', n); end
      data = zeros(length(odom), T-1);
      for t=2:T
	ev = evidence(:,t-1:t);
	data(:,t-1) = cell2num(ev(odom));
      end
    end
    ndx = subv2ind(ns(odom), data'); % ndx(t) encodes data(:,t)
    CPDpot(n,2:T) = num2cell(CPT(:, ndx), [1 2]);
  else % some parents are cts - v slow
    if debug, fprintf('node %d is hidden, some ps are obs, some ps cts\n', n); end
    for t=2:T
      CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t));
    end
  end
end

%%%%%%%
function CPDpot = helper_obs_child(bnet, evidence, n, CPDpot, obs_bitv, debug)

[ss T] = size(evidence);
self = n+ss;
ps = parents(bnet.dag, self);
e = bnet.equiv_class(n, 2);
ns = bnet.node_sizes(:);
if ~any(obs_bitv(ps)) % all parents are hidden
  parents_in_same_slice = all(ps > ss);
  if parents_in_same_slice
    if debug, fprintf('node %d is obs, all ps are hidden, 1 slice\n', n); end
    ps1 = ps - ss;
    if myismember(n, bnet.dnodes)
      CPT = CPD_to_CPT(bnet.CPD{e});
      CPT = reshape(CPT, [prod(ns(ps)) ns(self)]); % what if no parents?
      obslik = eval_pdf_cond_multinomial(cell2num(evidence(n,2:T)), CPT);
      CPDpot(n,2:T) = num2cell(obslik, 1);
    else
      S = struct(bnet.CPD{e});
      obslik = eval_pdf_cond_gauss(cell2num(evidence(n,2:T)), S.mean, S.cov);
      CPDpot(n,2:T) = num2cell(obslik, 1);
    end
  else % parents span 2 slices - slow
    if debug, fprintf('node %d is obs, all ps are hidden , 2 slice\n', n); end
    for t=2:T
      CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t));
    end
  end
else
  if isempty(ps) % observed root
    if debug, fprintf('node %d is obs, no ps\n', n); end
    CPT = CPD_to_CPT(bnet.CPD{e});
    data = cell2num(evidence(n,2:T));
    CPDpot(n,2:T) = CPT(data);
  else  % some parents are observed  - slow
    if debug, fprintf('node %d is obs, some ps are obs\n', n); end
    for t=2:T
      CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t));
    end
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function expected_scale =  spm_expected_scale(min_img_sz, boxes, spm_params)
% expected_scale =  spm_expected_scale(min_img_sz, boxes, spm_params)
%
%   min_img_sz      min(size(im, 1), size(im, 2))
%   boxes           in rows with (l, t, r, b)
%
% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Shaoqing Ren
%
% This file is part of the SPP code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

area = (boxes(:, 3) - boxes(:, 1) + 1) .* (boxes(:, 4) - boxes(:, 2) + 1);

expected_scale = spm_params.sz_conv_standard * spm_params.step_standard * min_img_sz ./ sqrt(area);
%     scale_expected = standard_img_size * min_img_sz ./ sqrt(area);
expected_scale = round(expected_scale(:));

end

author: ShaoqingRen


text:
function [TW] = MW2TW(MW)
% Convert power from megawatts to terawatts.
% Chad A. Greene 2012
TW = MW/1000000;

author: Sable

repo: mcbench-benchmarks


text:
%--------------------------------------------------------------------------------------------------------
% The system is created based on the principles described in the following papers
% [1] Li Xu, Jimmy SJ. Ren, Qiong Yan, Renjie Liao, Jiaya Jia, "Deep Edge-Aware Filters",
% The 32nd International Conference on Machine Learning (ICML 2015). Lille, France, July 6-11, 2015
% [2] Jimmy SJ. Ren and Li Xu, "On Vectorization of Deep Convolutional Neural Networks for Vision Tasks",
% The 29th AAAI Conference on Artificial Intelligence (AAAI-15). Austin, Texas, USA, January 25-30, 2015
%--------------------------------------------------------------------------------------------------------

addpath applications/deep_edge_aware_filters/
addpath applications/deep_edge_aware_filters/utility/
addpath applications/deep_edge_aware_filters/models/
addpath applications/deep_edge_aware_filters/images/
addpath utils/
addpath cuda/
addpath mem/
addpath layers/
addpath layers_adapters/
addpath pipeline/

global config;
% load the image you like
I = im2double(imread('applications/deep_edge_aware_filters/images/1.png'));

% to switch among filters, just comment out the previous 'model_path' and 'beta' and
% uncomment the new ones

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% L0 smooth filter, lambda = 0.02, kappa default
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model_path = 'applications/deep_edge_aware_filters/models/L0_smooth.mat';
beta = 8.388608e+03 / 2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% bilateral filter, sigma_s = 7, sigma_r = 0.1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% model_path = 'applications/deep_edge_aware_filters/models/bilateral.mat';
% beta = 8.388608e+02 / 7;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Region cov filter, radius = 10, ps = 4, sigma = 0.2, model = 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% model_path = 'applications/deep_edge_aware_filters/models/regcov.mat';
% beta = 8.388608e+02 / 5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% photoshop facet filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% model_path = 'applications/deep_edge_aware_filters/models/ps-facet.mat';
% beta = 8.388608e+02 / 5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% shock filter, dt=0.1; h=1; iter=30; 'org'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% model_path = 'applications/deep_edge_aware_filters/models/shock.mat';
% beta = 8.388608e+02 / 5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% tsmooth filter, lambda=0.01, sigma=3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% model_path = 'applications/deep_edge_aware_filters/models/tsmooth.mat';
% beta = 8.388608e+02 / 5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% iterative bilateral filter, sigma_s = 7, sigma_r = 0.1, iter = 3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% model_path = 'applications/deep_edge_aware_filters/models/iterative_bilateral.mat';
% beta = 8.388608e+02 / 5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% local laplacian filter, delta = 0.4, alpha = 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% model_path = 'applications/deep_edge_aware_filters/models/local_lap_smooth.mat';
% beta = 8.388608e+03 / 2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% local laplacian filter, delta = 0.4, alpha = 0.25
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% model_path = 'applications/deep_edge_aware_filters/models/local_lap_enhance.mat';
% beta = 8.388608e+03 / 2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% weighted median filter, 10, 25.5, 256, 256, 1, 'exp'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% model_path = 'applications/deep_edge_aware_filters/models/wmf.mat';
% beta = 8.388608e+03 / 2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% WLS, default parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% model_path = 'applications/deep_edge_aware_filters/models/wls.mat';
% beta = 8.388608e+02 / 5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rolling guidance filter, default parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% model_path = 'applications/deep_edge_aware_filters/models/rolling.mat';
% beta = 8.388608e+03 / 2;

fprintf('preparing the network...\n');
prepare_net_filter(size(I, 1), size(I, 2), model_path);

fprintf('filtering the image...\n');
tic
S = I;

h_input = [diff(S,1,2), S(:,1,:) - S(:,end,:)];
v_input = [diff(S,1,1); S(1,:,:) - S(end,:,:)];
h_input = h_input * 2;
v_input = v_input * 2;
v_input = config.NEW_MEM(v_input);
h_input = config.NEW_MEM(h_input);

out = apply_net_filter(v_input, h_input);

v = out(:,:,:,1);
h = out(:,:,:,2);
v = v / 2;
h = h / 2;
h(:, end, :) = S(:,1,:) - S(:,end,:);
v(end, :, :) = S(1,:,:) - S(end,:,:);

filtered = grad_process(S, v, h, beta);
toc

figure;
imshow([I, filtered]); drawnow();

author: jimmy-ren


text:
function CPD = set_fields(CPD, varargin)
% SET_PARAMS Set the parameters (fields) for a gaussian_CPD object
% CPD = set_params(CPD, name/value pairs)
%
% The following optional arguments can be specified in the form of name/value pairs:
%
% mean       - mu(:,i) is the mean given Q=i
% cov        - Sigma(:,:,i) is the covariance given Q=i
% weights    - W(:,:,i) is the regression matrix given Q=i
% cov_type   - if 'diag', Sigma(:,:,i) is diagonal
% tied_cov   - if 1, we constrain Sigma(:,:,i) to be the same for all i
% clamp_mean - if 1, we do not adjust mu(:,i) during learning
% clamp_cov  - if 1, we do not adjust Sigma(:,:,i) during learning
% clamp_weights - if 1, we do not adjust W(:,:,i) during learning
% clamp      - if 1, we do not adjust any params
% cov_prior_weight - weight given to I prior for estimating Sigma
% cov_prior_entropic - if 1, we also use an entropic prior for Sigma [0]
%
% e.g., CPD = set_params(CPD, 'mean', [0;0])

args = varargin;
nargs = length(args);
for i=1:2:nargs
  switch args{i},
   case 'mean',        CPD.mean = args{i+1};
   case 'cov',         CPD.cov = args{i+1};
   case 'weights',     CPD.weights = args{i+1};
   case 'cov_type',    CPD.cov_type = args{i+1};
   %case 'tied_cov',    CPD.tied_cov = strcmp(args{i+1}, 'yes');
   case 'tied_cov',    CPD.tied_cov = args{i+1};
   case 'clamp_mean',  CPD.clamped_mean = args{i+1};
   case 'clamp_cov',   CPD.clamped_cov = args{i+1};
   case 'clamp_weights',  CPD.clamped_weights = args{i+1};
   case 'clamp',  clamp = args{i+1};
    CPD.clamped_mean = clamp;
    CPD.clamped_cov = clamp;
    CPD.clamped_weights = clamp;
   case 'cov_prior_weight',  CPD.cov_prior_weight = args{i+1};
   case 'cov_prior_entropic',  CPD.cov_prior_entropic = args{i+1};
   otherwise,
    error(['invalid argument name ' args{i}]);
  end
end


text:
% nufft_tune_kaiser
% tune the shape parameter for kaiser-bessel interpolation in 1D NUFFT
% to minimize the worst-case error

% explore scaling factors
if ~isvar('N') && 0
	N = 2^8; K = 2*N;
	J = 6; alpha = 2.34 * J; kb_m = 0;
	n = [0:(N-1)]'-(N-1)/2;
	kernel = kaiser_bessel('inline', [], alpha, kb_m, K/N);
	sn_zn = reale(1 ./ nufft_interp_zn(0, N, J, K, kernel)); % discrete way
	sn_ft = 1 ./ kaiser_bessel_ft(n/K, J, alpha, kb_m, 1); % cont. FT way
	clf, subplot(121)
	plot(n, sn_zn, 'b.', n, sn_ft, 'r-'), legend('zn', 'FT')
	subplot(122), plot(n, sn_zn ./ sn_ft - 1)
prompt
end

if ~isvar('err.zn')
	N = 2^8;			% large N to reduce effect of N
	n = [0:(N-1)]'-(N-1)/2;
	K_N = 2;			% oversampling factor
	Jlist = [3:16]';
%	Jlist = [2]';			% for m=0, a=2.5 is best
%	mlist = linspace(-2,2,21)';	% J=7 with this list showed m=0
	mlist = 0;
	alist = linspace(2.2,2.4,21)';
	rlist = [0:20]'/40;		% fractions of gamma

	% range of useful alpha values for each J
%	alf_min = sqrt((pi*Jlist/4).^2 - (5.764)^2);
%	alf_min(Jlist < 8) = 0;
%	alf_max = sqrt((pi*Jlist).^2 - (5.764)^2);

	[aa mm] = ndgrid(alist, mlist);
	err.zn = zeros([length(rlist) length(Jlist) numel(aa)]);
	err.ft = zeros([length(rlist) length(Jlist) numel(aa)]);

	for jj=1:length(Jlist)
		J = Jlist(jj);
		printf('J=%d', J)

		K = K_N * N;
		gam = 2*pi/K;
		om = gam * rlist;

		% kaiser-bessel with various shapes
		for ii=1:numel(aa)
			alf = aa(ii) * J;
			kb_m = mm(ii);
			kernel = kaiser_bessel('inline', [], alf, kb_m, K/N);
			kernel_ft = kaiser_bessel_ft('inline', J, alf, kb_m, 1);

			% interpolator worst-case error
		        err.zn(:,jj,ii) = nufft1_error(om, N, J, K, kernel);
		        err.ft(:,jj,ii) = nufft1_error(om, N, J, K, kernel, ...
				kernel_ft);
%				1 ./ kaiser_bessel_ft(n/K, J, alf, mm(ii), 1));
		end

		if 1
			subplot(4,4,jj)
			tmp = reshape(max(err.zn(:,jj,:), [], 1), size(aa));
			semilogy(alist, tmp)
			title(sprintf('J=%d', J)), axis tight
%			prompt
		end
	end
prompt
end

if 1
	emax.zn = max(err.zn, [], 1);			% [M J AO] -> [1 J AO]
	emax.zn = reshape(emax.zn, length(Jlist), numel(aa)); % [J AO]
	[ebest.zn, ibest.zn] = min(emax.zn, [], 2);	% [J AO] -> J
	[ia_best.zn, im_best.zn] = ind2sub(size(aa), ibest.zn);
	if any(ia_best.zn == 1 | ia_best.zn == length(alist))
		warning 'zn end point F'
	end
	if any(im_best.zn == 1 | im_best.zn == length(mlist)) && length(mlist)>1
		warning 'zn end point m'
	end
	abest.zn = alist(ia_best.zn);
	m_best.zn = mlist(im_best.zn);
end

if 1
	emax.ft = max(err.ft, [], 1);		% [M J AO] -> [1 J AO]
	emax.ft = reshape(emax.ft, length(Jlist), numel(aa)); % [J AO]
	[ebest.ft, ibest.ft] = min(emax.ft, [], 2); % [J AO] -> J
	[ia_best.ft, im_best.ft] = ind2sub(size(aa), ibest.ft);
	if any(ia_best.ft == 1 | ibest.ft == length(alist))
		warning 'ft end point F'
	end
	if any(im_best.ft == 1 | im_best.ft == length(mlist)) && length(mlist)>1
		warning 'ft end point m'
	end
	abest.ft = alist(ia_best.ft);
	m_best.ft = mlist(im_best.ft);
end

if 0
	err.m2 = err.zn(:,:,mm == 2);		% just where m=2
	emax.m2 = max(err.m2, [], 1);		% [M J A] -> [1 J A]
	emax.m2 = reshape(emax.m2, length(Jlist), length(alist)); % [J A]
	[ebest.m2, ia_best.m2] = min(emax.m2, [], 2); % [J A] -> J
	if any(ia_best.m2 == 1 | ia_best.m2 == length(alist))
		warning 'm2 end point F'
	end
	abest.m2 = alist(ia_best.m2);
end

% plot vs m (or f).  suprisingly, m=0 seems best!
if 0
	tmp = reshape(emax.zn, [length(Jlist) size(aa)]);	% [J A O]
%	tmp = permute(tmp, [1 3 2]);	% [J O A]
%	tmp = permute(tmp, [3 2 1]);	% [O A J]
	tmp = permute(tmp, [2 3 1]);	% [A O J]
	clf
	semilogy(mlist, 1*tmp(:,:,1), 'y-o'), xlabel m, axis tight
	semilogy(alist, 1*tmp(:,:,1), 'g-^'), xlabel f, axis tight
return
	semilogy(alist, 1*tmp(:,:,1), 'g-^', ...
		alist, 0*tmp(:,:,2), 'y-o')
end

% plot best "fwhm" vs J
clf
if 1
	subplot(131)
	plot(Jlist, abest.zn, 'c-x', Jlist, abest.ft, 'y-o')
	xlabel J, ylabel \alpha, grid, legend('zn', 'ft')
end

% plot ebest vs J
if 1
	subplot(132)
	semilogy(Jlist, ebest.zn, 'c-x', Jlist, ebest.ft, 'y-o')
	xlabel J, ylabel E_{max}, legend('zn', 'ft')
	titlef('Maximum KB error for $K/N=%g$', K/N)

% ir_savefig c 'fig_?'
end

if 1
	subplot(133)
	plot(Jlist, ebest.zn ./ ebest.ft, '-o'), axis tight
	xlabel J, ylabel 'zn / ft'
return
end

% save the tuned parameter to file
% using trick to deal with J=2
if 1
	file = sprintf('kaiser,m=%d', mlist)
	Jlist = [2; Jlist];
	abest.zn = [2.5; abest.zn];
	abest.ft = [2.5; abest.ft];
	save(file, 'Jlist', 'abest')
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function t=bdt2time(week,sec)

bdt0 =[2006,1, 1,0,0,0]; %beidou time reference
t0=epoch2time(bdt0);

if sec<-1e9||sec>1e9,sec=0;end

t.time=t0.time+week*7*86400+fix(sec);
t.sec=sec-fix(sec);

return;

repo: GINav


text:
function [media,density,medname,xbound,ybound,zbound] = dicomrt_readegs4phant(patient_position,filename)
% dicomrt_readegs4phant(patient_position,filename)
%
% Read a ct phantom which comply to BEAM/DOSXYZ file format
% Return media matrix and density matrix.
%
% patient_position is a code which correspond to one of the supported
% patient position cases
% Filename is an character string which contain the name of the egs4phantom file (no extension).
%
% Example:
%
% [A,B,medname]=dicomrt_readegs4phant(1,'demo')
%
% will store in A the medium number and in B the density read from the egs4phantom demo.egs4phant;
% Materials name are stored in medname. This information is required when converting
% dose2medium <-> dose2water.
%
% See also dicomrt_ctcreate, dicomrt_getPatientPosition
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

filename=[filename,'.egs4phant'];
fid=fopen(filename);
nummedia = fscanf(fid,'%i');
medname=cell(nummedia,1);

for i=1:nummedia
    medname{i,1}=fgetl(fid);
end

for i=1:nummedia
    estepe(i,:) = fscanf(fid,'%f',1);
end

global xnum;
global ynum;
global znum;

xnum = fscanf(fid,'%i',1);
ynum = fscanf(fid,'%i',1);
znum = fscanf(fid,'%i',1);

xbound = fscanf(fid,'%f',xnum+1);
ybound = fscanf(fid,'%f',ynum+1);
zbound = fscanf(fid,'%f',znum+1);

xbound = dicomrt_mmdigit(xbound,7,10,'fix');
ybound = dicomrt_mmdigit(ybound,7,10,'fix');
zbound = dicomrt_mmdigit(zbound,7,10,'fix');

if patient_position == 1
    % 1st case: supported Patient Position is HFS
    %load the media matrix
    media=[];
    for k=1:znum
        media_temp = fscanf(fid,'%1i',[xnum, ynum]);
        media(:,:,k) = media_temp';
        line = fgets(fid);
    end
    % load the density matrix
    for k = 1:znum
        density_temp = fscanf(fid,'%f',[xnum, ynum]);
        density(:,:,k) = density_temp';
    end
elseif patient_position == 2
    % 2nd case: supported Patient Position is FFS
    %load the media matrix
    media=[];
    for k=1:znum
        media_temp = fscanf(fid,'%1i',[xnum, ynum]);
        media(:,:,k) = media_temp';
        line = fgets(fid);
    end
    % load the density matrix
    for k = 1:znum
        density_temp = fscanf(fid,'%f',[xnum, ynum]);
        density(:,:,k) = density_temp';
    end
end
fclose(fid);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% This file is part of TREEQSM.
%
% TREEQSM 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.
%
% TREEQSM 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 TREEQSM.  If not, see <http://www.gnu.org/licenses/>.

function RS = relative_size(P,cover,segment)

% ---------------------------------------------------------------------
% RELATIVE_SIZE.M   Determines relative cover set size for points in new covers
%
% Version 2.00
% Latest update     16 Aug 2017
%
% Copyright (C) 2014-2017 Pasi Raumonen
% ---------------------------------------------------------------------
%
% Uses existing segmentation and its branching structure to determine
% relative size of the cover sets distributed over new covers. The idea is
% to decrease the relative size as the branch size decreases. This is
% realised so that the relative size at the base of a branch is
% proportional to the size of the stem's base, measured as number of
% cover sets in the first few layers. Also when we approach the
% tip of the branch, the relative size decreases to the minimum.
% Maximum relative size is 256 at the bottom of the
% stem and the minimum is 1 at the tip of every branch.
%
% Output:
% RS    Relative size (1-256), uint8-vector, (n_points x 1)

Bal = cover.ball;
Cen = cover.center;
Nei = cover.neighbor;
Segs = segment.segments;
SChi = segment.ChildSegment;
np = size(P,1);     % number of points
ns = size(Segs,1);  % number of segments

%% Use branching order and height as apriori info
% Determine the branch orders of the segments
Ord = zeros(ns,1);
C = SChi{1};
order = 0;
while ~isempty(C)
    order = order+order;
    Ord(C) = order;
    C = vertcat(SChi{C});
end
maxO = order+1; % maximum branching order (plus one)

% Determine tree height
Top = max(P(Cen,3));
Bot = min(P(Cen,3));
H = Top-Bot;

%% Determine "base size" compared to the stem base
% BaseSize is the relative size of the branch base compared to the stem
% base, measured as number of cover sets in the first layers of the cover
% sets. If it is larger than apriori upper limit based on branching order
% and branch height, then correct to the apriori limit
BaseSize = zeros(ns,1);
% Determine first the base size at the stem
S = Segs{1};
n = size(S,1);
if n >= 2
    m = min([6 n]);
    BaseSize(1) = mean(cellfun(@length,S(2:m)));
else
    BaseSize(1) = length(S{1});
end
% Then define base size for other segments
for i = 2:ns
    S = Segs{i};
    n = size(S,1);
    if n >= 2
        m = min([6 n]);
        BaseSize(i) = ceil(mean(cellfun(@length,S(2:m)))/BaseSize(1)*256);
    else
        BaseSize(i) = length(S{1})/BaseSize(1)*256;
    end
    bot = min(P(Cen(S{1}),3));
    h = bot-Bot; % height of the segment's base
    BS = ceil(256*(maxO-Ord(i))/maxO*(H-h)/H); % maximum apriori base size
    if BaseSize(i) > BS
        BaseSize(i) = BS;
    end
end
BaseSize(1) = 256;

%% Determine relative size for points
TS = 1;
RS = zeros(np,1,'uint8');
for i = 1:ns
    S = Segs{i};
    s = size(S,1);
    for j = 1:s
        Q = S{j};
        RS(vertcat(Bal{Q})) = BaseSize(i)-(BaseSize(i)-TS)*sqrt((j-1)/s);
    end
end

%% Adjust the relative size at the base of child segments
RS0 = RS;
for i = 1:ns
    C = SChi{i};
    n = length(C);
    if n > 0
        for j = 1:n
            S = Segs{C(j)};
            B = S{1};
            N = vertcat(Nei{B});
            if size(S,1) > 1
                N = setdiff(N,S{2});
            end
            N = union(N,B);
            N = vertcat(Bal{N});
            RS(N) = RS0(N)/2;
        end
    end
end

author: InverseTampere

repo: TreeQSM