text:
function out = EN_DistributionEntropy(y,histOrKS,numBins,olremp)
% EN_DistributionEntropy    Distributional entropy.
%
% Estimates of entropy from the distribution of a data vector. The
% distribution is estimated either using a histogram with numBins bins, or as a
% kernel-smoothed distribution, using the ksdensity function from Matlab's
% Statistics Toolbox with width parameter, w (specified as the iunput numBins).
%
% An optional additional parameter can be used to remove a proportion of the
% most extreme positive and negative deviations from the mean as an initial
% pre-processing.
%
%---INPUTS:
%
% y, the input time series
%
% histOrKS: 'hist' for histogram, or 'ks' for ksdensity
%
% numBins: (*) (for 'hist'): an integer, uses a histogram with that many bins
%          (*) (for 'ks'): a positive real number, for the width parameter for
%                       ksdensity (can also be empty for default width
%                                       parameter, optimum for Gaussian)
%
% olremp [opt]: the proportion of outliers at both extremes to remove
%               (e.g., if olremp = 0.01; keeps only the middle 98% of data; 0
%               keeps all data. This parameter ought to be less than 0.5, which
%               keeps none of the data).
%               If olremp is specified, returns the difference in entropy from
%               removing the outliers.

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

doPlot = 0; % plot outputs to figure

% ------------------------------------------------------------------------------
%% Check inputs
% ------------------------------------------------------------------------------
if nargin < 2 || isempty(histOrKS)
    histOrKS = 'hist'; % use histogram by default
end
if nargin < 3 % (can be empty for default width for ksdensity)
    numBins = 10; % use 10 bins
end
if nargin < 4
    olremp = 0;
end

% ------------------------------------------------------------------------------
% (1) Remove outliers?
% ------------------------------------------------------------------------------
if olremp ~= 0
    yHat = y(y >= quantile(y,olremp) & y <= quantile(y,1-olremp));
    if isempty(yHat)
        % removed the entire time series?!
        % shouldn't be possible for good values of olremp with equality
        % in the above inequalities
        out = NaN; return
    else
        % Return the difference in entropy from removing outliers
        out = EN_DistributionEntropy(y,histOrKS,numBins) - ...
                EN_DistributionEntropy(yHat,histOrKS,numBins);
        return
    end
end

% ------------------------------------------------------------------------------
% (2) Form the histogram
% ------------------------------------------------------------------------------
switch histOrKS
case 'hist' % Use histogram to calculate pdf
    if isnumeric(numBins)
        [px,binEdges] = histcounts(y,numBins,'Normalization','probability');
    else
        [px,binEdges] = histcounts(y,'BinMethod',numBins,'Normalization','probability');
    end
    % Compute bin centers:
    xr = mean([binEdges(1:end-1); binEdges(2:end)]);
    % Compute bin widths:
    binWidths = diff(binEdges);

case 'ks' % Use ksdensity to calculate pdf
    if isempty(numBins)
        [px, xr] = ksdensity(y,'function','pdf'); % selects optimal width
    else
        [px, xr] = ksdensity(y,'width',numBins,'function','pdf'); % uses specified width
    end
    binWidths = ones(1,length(px))*(xr(2)-xr(1));

otherwise
    error('Unknown distribution method -- specify ''ks'' or ''hist''') % error; must specify 'ks' or 'hist'
end

if doPlot
    figure('color','w'); box('on');
    plot(xr,px,'k')
end

% ------------------------------------------------------------------------------
% (3) Compute the entropy sum and return it as output
% ------------------------------------------------------------------------------
% 0*log0 = 0:
out = -sum(px(px>0).*log(px(px>0)./binWidths(px>0)));

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function segment_length = p07_boundary_segment_length ( segment_index, h )

%*****************************************************************************80
%
%% P07_BOUNDARY_SEGMENT_LENGTH returns boundary segment lengths in problem 07.
%
%  Discussion:
%
%    No attempt has been made here to accurately compute a value of N
%    which would guarantee that the boundary would be divided into pieces
%    of length no more than H.  The curve is a little too complicated
%    to make this easy to do.
%
%    Moreover, the points that will be generated will only be equally
%    spaced in their X argument, not in their arc length.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer SEGMENT_INDEX, the index of one of the boundary segments.
%
%    Input, real H, the suggested spacing between points.
%
%    Output, integer SEGMENT_LENGTH, the number of points in the segment.
%
  if ( h <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P07_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );
    fprintf ( 1, '  Nonpositive H = %f\n', h );
    error ( 'P07_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );
  end

  if ( segment_index == 1 )

    n = round ( 10.0 * pi / h );
    n = max ( n, 13 );
    segment_length = n;

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P07_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );
    fprintf ( 1, '  Illegal SEGMENT_INDEX = %d\n', segment_index );
    error ( 'P07_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function qwgw_test07 ( )

%*****************************************************************************80
%
%% TEST07 tests QWGW for the Jacobi weight.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2014
%
%  Author:
%
%    John Burkardt
%

%
%  Set the quadrature interval and number of points.
%
  a = -1.0;
  b = +1.0;
  n = 5;
%
%  Set the weight function parameters.
%
  alpha = 0.25;
  beta = 0.75;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST07:\n' );
  fprintf ( 1, '  Compute points and weights for Gauss quadrature\n' );
  fprintf ( 1, '  with the Jacobi weight w(x) = (1-x^2)^alpha*(1+x)^beta\n' );
  fprintf ( 1, '  Order N = %d\n', n );
  fprintf ( 1, '  ALPHA = %g\n', alpha );
  fprintf ( 1, '  BETA =  %g\n', beta );
  fprintf ( 1, '  Interval = [%g,%g]\n', a, b );
%
%  Set the recursion coefficients.
%
  aj = zeros ( n, 1 );
  bj = zeros ( n, 1 );

  for j = 1 : n
    jr = j;
    aj(j) = ( beta - alpha ) * ( beta + alpha ) ...
      / ( alpha + beta + 2.0 * jr - 2.0 ) ...
      / ( alpha + beta + 2.0 * jr );
  end

  for j = 1 : n - 1
    jr = j;
    bj(j) = 4.0 * jr * ( alpha + jr ) * ( beta + jr ) ...
      * ( alpha + beta + jr ) ...
      / ( ( alpha + beta + 2.0 * jr )^2 - 1.0 ) ...
      / ( alpha + beta + 2.0 * jr )^2;
  end
  bj(n) = 0.0;

  bj(1:n) = sqrt ( bj(1:n) );

  mu0 = 2.0 ^ ( alpha + beta + 1.0 ) ...
    * gamma ( alpha + 1.0 ) * gamma ( beta + 1.0 ) ...
    / gamma ( alpha + beta + 2.0 );
%
%  Compute the points and weights.
%
  [ x, w ] = sgqf ( n, aj, bj, mu0 );

  r8vec_print ( n, x, '  Abscissas:' );
  r8vec_print ( n, w, '  Weights:' );

  return
end


text:
function PanoLabelIm(id, ANNO_ALL, config, options)
% load image
im_name = config.IMGLIST(id).name;
im = imread([config.folderName im_name '/' im_name '.jpg']);
[im_h, im_w, im_dim] = size(im);
im_o = im2double(im);
% alignment
if 1
imgSize = 320;
qError = 0.7;
global config;
config.NewProjFolder = './panoContext_code/';
config.lsdLocation = [config.NewProjFolder 'VpEstimation/lsd '];
config.bufferFolder = './Buffer_gt/';
if ~exist(config.bufferFolder, 'dir')
    mkdir(config.bufferFolder);
end
[ olines, vp, views, edges, panoEdge, score, angle] = panoEdgeDetection( im_o, imgSize, qError);
save(['./label_vp/pc/' im_name '.mat'], 'vp');
%keyboard
return
end
load(['./label_vp/pc/' im_name '.mat']);
[rotImg, R] = rotatePanorama(im_o, vp(3:-1:1,:));

if exist(['./label_cor/pc/' im_name '.mat'], 'file')
    gt = load(['./label_cor/pc/' im_name '.mat']);
    gt = gt.cor;
    keyboard
else

r_id = find([ANNO_ALL(id).ANNO3D.objects.type] == 3);
    pt = ANNO_ALL(id).ANNO3D.objects(r_id).points *ANNO_ALL(id).ANNO3D.R;
    pt = pt*R'; % algin pt also
    [ uv ] = xyz2uvN( pt ); % transform xyz to im coordinates
    [ uvcoord ] = uv2coords( uv, im_w, im_h );
%sort
[~,gt_id] = sort(uvcoord(:,1));
gt = uvcoord(gt_id,:);
for i = 1:2:size(gt,1)
    gt_t = gt(i:i+1,:);
    [~, gt_t_id] = sort(gt_t(:,2));
    gt(i:i+1,:) = gt_t(gt_t_id,:);
end
end
gt_all = [gt;
               gt(1,:);gt(3,:);gt(3,:);gt(5,:);
               gt(5,:);gt(7,:);gt(7,:);gt(1,:);
               gt(2,:);gt(4,:);gt(4,:);gt(6,:);
               gt(6,:);gt(8,:);gt(8,:);gt(2,:)];
    [ uv ] = coords2uv( gt_all, im_w, im_h );
    [ xyz ] = uv2xyzN( uv);
    [ lines ] = lineFromTwoPoint( xyz(1:2:end,:), xyz(2:2:end,:) );
    [ panoEdgeC2 ] = paintParameterLine2(lines, im_w, im_h, rotImg, 2 );

figure(1), axis off, hold off,
ha = tight_subplot(1, 1);
axes(ha(1)), imagesc(panoEdgeC2), axis image, hold on

if 1 % corner debug
disp('Click y if label this image');
[~,~,b] = ginput(1);
if b ~= 'y'
    return
end

% Allow user to input line segments; compute centers, directions, lengths
%disp('Draw Pano lines...')

% find concave lines
lines = [];
while 1
    disp(' ')
    disp('Draw concave lines or q to stop')
    [x1,y1,b] = ginput(1);
    if b=='q'
        break;
    end
    [x2,y2, ~] = ginput(1);
    % single line fitting
    vert = (x1+x2)/2;
    axes(ha(1)); h1 = plot([vert vert], [1 im_h], 'b');
    disp('Click y if accept this line')
    [~,~, b] = ginput(1);
    % save
    if b == 'y'
        lines(end+1,:) = [vert y1 1]; % 1 for concave signal
        lines(end+1,:) = [vert y2 1];
    else
        axes(ha(1)); delete(h1);
    end
end

% find convex lines
while 1
    disp(' ')
    disp('Draw convex lines or q to stop')
    [x1,y1,b] = ginput(1);
    if b=='q'
        break;
    end
    [x2,y2, ~] = ginput(1);
    % single line fitting
    vert = (x1+x2)/2;
    axes(ha(1)); h1 = plot([vert vert], [1 im_h], 'r');
    disp('Click y if accept this line')
    [~,~, b] = ginput(1);
    % save
    if b == 'y'
        lines(end+1,:) = [vert y1 2]; % 2 for convex signal
        lines(end+1,:) = [vert y2 2];
    else
        axes(ha(1)); delete(h1);
    end
end

% sort lines
[~,l_id] = sort(lines(:,1));
lines = lines(l_id,:);
for i = 2:2:size(lines,1)
    sub_l = lines(i-1:i,:);
    [~,l_id] = sort(sub_l,1);
    lines(i-1,:) = sub_l(l_id(1),:);
    lines(i,:) = sub_l(l_id(2),:);
end
end

% solve for pano
if size(lines,1) == 8
    [wall_d, x, f] = pano_line_solver(lines, im_w, options); % bos-shape case
else if size(lines,1) == 12
       [wall_d, x, f] = pano_line_solver_6(lines, im_w, options); % L-shape case
    else if size(lines,1) == 16
            [wall_d, x, f] = pano_line_solver_8(lines, im_w, options); % T-shape case
        end
    end
end

lines_ = lines;
% annotate box height
cor = [];
%c_h = 1.7;
flag = 0;
while 1
    disp('Draw corners or q to stop')
    [x1,y1,b] = ginput(1);
    axes(ha(1)); h1 = plot(x1, y1, 'xr');
    if b=='q'
        break;
    end
    [x2,y2] = ginput(1);
    axes(ha(1)); h3 = plot(x2, y2, 'xr');
    disp('Click y if accept this corner');
    [~,~,b] = ginput(1);
    if b == 'y'
        cor(end+1,:) = [x1 y1];
        cor(end+1,:) = [x2 y2];
        % find closet wall
        %keyboard
        vert = (x1+x2)/2;
        cor_d = abs(lines(:,1) - vert);
        [~, l_id] = min(cor_d);
        cor = [lines(l_id,1) cor(1,2); lines(l_id,1) cor(2,2)];
        % solve for all walls
        [~, cor_id] = sort(cor(:,2));
        cor = cor(cor_id,:);
        theta_y = pi*cor(2,2)/im_h-pi/2;
        d = abs(cot(theta_y));lines_ = lines;%c_h*abs(cot(theta_y));
        theta_y_ = pi/2 - pi*cor(1,2)/im_h;
        h = d*tan(theta_y_);%c_h + d*tan(theta_y_);
        d_m = d/wall_d((l_id+1)/2);
        line_ya = [];
        line_ya_ = [];
        for i = 1:2:size(lines,1)
            if i == l_id
                continue
            end
            line_d = d_m*wall_d((i+1)/2);
            line_theta_y = acot(line_d);%acot(line_d/c_h);
            line_y = (line_theta_y + pi/2)*im_h/pi;
            cor = [cor;lines(i,1) line_y];
            line_y_ = atan(h/line_d);%atan((h-c_h)/line_d);
            line_y_ = (pi/2-line_y_)*im_h/pi;
            cor = [cor;lines(i,1) line_y_ ];
            axes(ha(1)); h1 = plot([lines(i,1) lines(i,1)], [line_y line_y_], 'xr');
        end
        % sort cor
        [~, cor_id] = sort(cor(:,1));
        cor = cor(cor_id,:);
        for i = 1:2:size(cor,1)
                cor_sub = cor(i:i+1,:);
                [~, cor_id] = sort(cor_sub(:,2));
                cor(i,:) = cor_sub(cor_id(1),:);
                cor(i+1,:) = cor_sub(cor_id(2),:);
        end

        % visual full
        if size(cor,1) == 8
            cor_all = [cor;
                cor(1,:);cor(3,:);cor(3,:);cor(5,:);
                cor(5,:);cor(7,:);cor(7,:);cor(1,:);
                cor(2,:);cor(4,:);cor(4,:);cor(6,:);
                cor(6,:);cor(8,:);cor(8,:);cor(2,:)];
        else if size(cor,1) == 12
            cor_all = [cor;
                cor(1,:);cor(3,:);cor(3,:);cor(5,:);
                cor(5,:);cor(7,:);cor(7,:);cor(9,:);
                cor(9,:);cor(11,:);cor(11,:);cor(1,:);
                cor(2,:);cor(4,:);cor(4,:);cor(6,:);
                cor(6,:);cor(8,:);cor(8,:);cor(10,:);
                cor(10,:);cor(12,:);cor(12,:);cor(2,:)];
		else if size(cor,1) == 16
		    cor_all = [cor;
                	cor(1,:);cor(3,:);cor(3,:);cor(5,:);
                	cor(5,:);cor(7,:);cor(7,:);cor(9,:);
                	cor(9,:);cor(11,:);cor(11,:);cor(13,:);
                        cor(13,:);cor(15,:);cor(15,:);cor(1,:);
                	cor(2,:);cor(4,:);cor(4,:);cor(6,:);
                	cor(6,:);cor(8,:);cor(8,:);cor(10,:);
                	cor(10,:);cor(12,:);cor(12,:);cor(14,:);
			cor(14,:);cor(16,:);cor(16,:);cor(2,:)];
		end
            end
        end
        [ uv ] = coords2uv( cor_all, im_w, im_h );
        [ xyz ] = uv2xyzN( uv);
        [ lines ] = lineFromTwoPoint( xyz(1:2:end,:), xyz(2:2:end,:) );
        [ panoEdgeC ] = paintParameterLine2(lines, im_w, im_h, rotImg, 1 );
        axes(ha(1));imagesc(panoEdgeC); axis image; hold on;
        for i = 1:size(line,1)
            axes(ha(1)); plot([lines(i,1) lines(i,1)], [1 im_h], 'b');
            axes(ha(1)); plot([lines(i,1) lines(i,1)], [line_y line_y_], 'xr');
        end

        disp('Click y if accept this box');
        [~,~,b] = ginput(1);
        if b == 'y'
            flag = 1;
            break
        else
            axes(ha(1)); imagesc(rotImg), axis image, hold on
            for i = 1:size(lines,1)
                axes(ha(1)); plot([lines(i,1) lines(i,1)], [1 im_h], 'b');
            end
            cor = [];
        end
    else
        axes(ha(1)); delete(h1); delete(h3);
    end
    if flag
        break
    end
end

if 0
% sort cor
[~, cor_id] = sort(cor(:,1));
cor = cor(cor_id,:);
for i = 1:2:size(cor,1)
    cor_sub = cor(i:i+1,:);
    [~, cor_id] = sort(cor_sub(:,2));
    cor(i,:) = cor_sub(cor_id(1),:);
    cor(i+1,:) = cor_sub(cor_id(2),:);
end

% visual full
cor_all = [cor;
               cor(1,:);cor(3,:);cor(3,:);cor(5,:);
               cor(5,:);cor(7,:);cor(7,:);cor(1,:);
               cor(2,:);cor(4,:);cor(4,:);cor(6,:);
               cor(6,:);cor(8,:);cor(8,:);cor(2,:)];
[ uv ] = coords2uv( cor_all, im_w, im_h );
[ xyz ] = uv2xyzN( uv);
[ lines ] = lineFromTwoPoint( xyz(1:2:end,:), xyz(2:2:end,:) );
[ panoEdgeC ] = paintParameterLine2(lines, im_w, im_h, rotImg, 1 );
axes(ha(1));imagesc(panoEdgeC); axis image;
end

cor = [cor lines_(:,3)];

imwrite(panoEdgeC, ['./label_vis/pc/' im_name '.png'], 'png');
save(['./label_cor/pc/' im_name '.mat'], 'cor');

author: zouchuhang

repo: LayoutNet


text:
%CALIBRATIONMATRIXVALUES  Computes useful camera characteristics from the camera matrix
%
%     S = cv.calibrationMatrixValues(cameraMatrix, imageSize, apertureWidth, apertureHeight)
%
% ## Input
% * __cameraMatrix__ Input 3x3 camera matrix that can be estimated by
%   cv.calibrateCamera or cv.stereoCalibrate.
% * __imageSize__ Input image size `[w,h]` in pixels.
% * __apertureWidth__ Physical width in mm of the sensor.
% * __apertureHeight__ Physical height in mm of the sensor.
%
% ## Output
% * __S__ Struct with the following fields
%   * __fovx__ Output field of view in degrees along the horizontal sensor
%     axis.
%   * __fovy__ Output field of view in degrees along the vertical sensor axis.
%   * __focalLength__ Focal length of the lens in mm.
%   * __principalPoint__ Principal point `[cx,cy]` in mm.
%   * __aspectRatio__ Pixel aspect ratio `fy/fx`.
%
% The function computes various useful camera characteristics from the
% previously estimated camera matrix.
%
% ### Note
% Do keep in mind that the unity measure 'mm' stands for whatever unit of
% measure one chooses for the chessboard pitch (it can thus be any value).
%
% See also: cv.calibrateCamera, cv.stereoCalibrate
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
classdef TestRandMVNormal
    %TestRandMVNormal

    methods (Static)
        function test_1
            d = 3;
            mu = rand(1,d);
            sigma = eye(d);
            nsamples = 50;
            samples = cv.randMVNormal(mu, sigma, nsamples);
            validateattributes(samples, {'numeric'}, {'size',[nsamples d]});
            mean(samples);  % should be close to mu
            cov(samples);   % should be close to sigma
        end

        function test_error_argnum
            try
                cv.randMVNormal();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end


text:
function title = p23_title ( )

%*****************************************************************************80
%
%% P23_TITLE returns a title for problem 23.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'The De Jong Function F2';

  return
end


text:
% Sympoly demos

%% Various ways to create a sympoly

% A scalar (zero) sympoly
z = sympoly;

% Scalar sympolys 'x', 'y', 'u', 'v' created in the current workspace
sympoly x y u v

% A sympoly (identity matrix) array. The numeric element format is
% specified by the command window format style.
format short g
ayuh = sympoly(eye(3));

% Use deal to replicate a sympoly into several
[a,b] = deal(sympoly);

% Deal can also create a sympoly array
S(1:2) = deal(sympoly('x'));

% As can repmat
R = repmat(sympoly('x'),2,3);

whos

%% Arithmetic between sympolys, add, subtract, multiply, divide.

% add 1 to x
1 + x

%%

% double times a sympoly
2*y

%%

% subtraction, and a simple power
(x - y)^2

%%

% More complex expressions
(x - 2*y)^3/x + sqrt(y^3)

%% Synthetic division
[quotient,remainder] = syndivide(x^2+2*x-1,x+1)

%% Arrays of sympolys
[x , y ; 1 , x+y]

%%

% Arrays of sympolys
v = [1 x y x+y]

%% matrix multiplication
A = v*v'
B = v'*v
%% Selective extraction of terms
% The second term
terms(A,2)

%%
terms(A,x^2,'extract')

%%
% Delete a term
p = (1 + x^2 + x^7)^3
terms(p,x^2,'delete')

%% Selective deletion of terms
B = terms(A,x,'extract')

%%
% Operations on arrays
sympoly lambda
(rand(3) - lambda*eye(3))

%% Even eigenvalues, using det, then roots
roots(det(hilb(4) - lambda*eye(4)))

%% Sum on any dimension
sum(v,2)

%% And prod
prod(A(:))

%% Orthogonal polynomials from a variety of familes

% 3rd and 4th order Legendre polynomials
p3 = orthpoly(3,'legendre')
p4 = orthpoly(4,'legendre')

%%

% Orthogonal polynomials are orthogonal over the proper domain
defint(p3*p4,'x',[-1,1])

%%

% 2nd and 5th order Jacobi polynomials
p2 = orthpoly(2,'jacobi',2,3)
p5 = orthpoly(5,'jacobi',2,3)

%%

% Orthogonal polynomials are orthogonal over the proper domain.
% Numerical issures left this just eps shy from zero.
defint(p2*p5*(1-x)^2*(1+x)^5,'x',[-1,1])

%% Roots of the derivative of a sympoly
sort(roots(diff(orthpoly(6,'cheby2'))))

%% Error propagation through a sympoly

%  Given a unit Normal N(0,1) random variable, compute the
%  mean and variance of p(x) = 3*x + 2*x^2 - x^3

sympoly x
[polymean, polyvar] = polyerrorprop(3*x + 2*x^2 - x^3,'x',0,1)

%%

% Compute the mean and variance of x*y + 3*y^3, where x and y are
% respectively N(mux,sx^2), and N(muy,sy^2)

sympoly x y mux muy sx sy
[polymean,polyvar] = polyerrorprop(x*y+3*y^3,{'x' 'y'},[mux,muy],[sx,sy])

%% A simple construction for a Newton-Cotes integration rule
% Here I'll generate Simpson's 3/8 rule.
%
% <http://mathworld.wolfram.com/Simpsons38Rule.html |Simpson's 3/8 rule|>
%
M = vander(0:3);
sympoly x f0 f1 f2 f3
% an interpolating polynomial on this set of points
% { (0,f0), (1,f1), (2,f2), (3,f3) }
P_of_x = [x^3, x^2, x, 1]*pinv(M)*[f0;f1;f2;f3];

%% sympoly uses the command window format style to write out the coefficients
% Here, I'll force it to be in rational form
format rat

%%
% integrate the polynomial over its support
defint(P_of_x,'x',[0 3])

%%
% Or here, a 4 point open Newton-Cotes rule
M = vander(1:4);
sympoly x f1 f2 f3 f4
% an interpolating polynomial on this set of points
% { (1,f1), (2,f2), (3,f3) (4,f4) }
P_of_x = [x^3, x^2, x, 1]*pinv(M)*[f1;f2;f3;f4]

%%
% integrate the polynomial over the full domain of the rule
defint(P_of_x,'x',[0 5])

author: Sable

repo: mcbench-benchmarks


text:
function value = composite_abscissa ( order, i )

%*****************************************************************************80
%
%% COMPOSITE_ABSCISSA returns the I-th abscissa of a composite rule.
%
%  Discussion:
%
%    Our convention is that the abscissas are numbered from left to
%    right, that the interval is [-1,1], and the the abscissas are
%    evenly spaced, and that, except for the order 1 rule, the
%    endpoints are included.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ORDER, the order of the rule.
%
%    Input, integer I, the index of the desired abscissa.
%    1 <= I <= ORDER.
%
%    Output, real VALUE, the value of the I-th abscissa in the
%    rule of order ORDER.
%
  a = -1.0;
  b = +1.0;

  if ( order < 1 )
    value = - Inf;
  elseif ( i < 1 | order < i )
    value = - Inf;
  elseif ( order == 1 )
    value = 0.5 * ( a + b );
  else
    value = ( ( order - i     ) * a   ...
            + (         i - 1 ) * b ) ...
            / ( order     - 1 );
  end

  return
end


text:
fprintf(1,'Script that merges two "Cabib_Results.mat" data sets of the same camera into a single dataset\n')

dir;

cont = 1;
while cont
    data_set1 = input('Filename of the first dataset (with complete path if necessary): ','s');
    cont = ((exist(data_set1)~=2));
    if cont,
        fprintf(1,'File not found. Try again.\n');
    end;
end;
cont = 1;
while cont
    data_set2 = input('Filename of the second dataset (with complete path if necessary): ','s');
    cont = ((exist(data_set2)~=2));
    if cont,
        fprintf(1,'File not found. Try again.\n');
    end;
end;

load(data_set1); % part1\Calib_Results;

shift = n_ima;

load(data_set2); % part2\Calib_Results;

active_images2 = active_images;
n_ima2 = n_ima;

for kk = 1:n_ima

eval(['X_' num2str(kk+shift) ' = X_' num2str(kk) ';']);

eval(['dX_' num2str(kk+shift) ' = dX_' num2str(kk) ';']);
eval(['dY_' num2str(kk+shift) ' = dY_' num2str(kk) ';']);

eval(['wintx_' num2str(kk+shift) ' = wintx_' num2str(kk) ';']);
eval(['winty_' num2str(kk+shift) ' = winty_' num2str(kk) ';']);

eval(['x_' num2str(kk+shift) ' = x_' num2str(kk) ';']);
eval(['y_' num2str(kk+shift) ' = y_' num2str(kk) ';']);

eval(['n_sq_x_' num2str(kk+shift) ' = n_sq_x_' num2str(kk) ';']);
eval(['n_sq_y_' num2str(kk+shift) ' = n_sq_y_' num2str(kk) ';']);

eval(['omc_' num2str(kk+shift) ' = omc_' num2str(kk) ';']);
eval(['Tc_' num2str(kk+shift) ' = Tc_' num2str(kk) ';']);

end;

load(data_set1); % part1\Calib_Results;

n_ima = n_ima + n_ima2;
active_images = [active_images active_images2];

no_image = 1;

% Recompute the error (in the vector ex):
comp_error_calib;

fprintf('The two calibration datasets are now merged. You are now ready to run calibration. \n');

repo: ekfmonoslam


text:
function [decimaldates1,decimaldates2,stringdates1,stringdates2,stringdates3,Fstartlocation,Fendlocation]=gendates(names,lags,frequency,startdate,enddate,Fstartdate,Fenddate,Fcenddate,Fendsmpl,F,CF,favar)

% function [decimaldates1,decimaldates2,stringdates1,stringdates2,stringdates3,Fstartlocation,Fendlocation]=gendates(names,lags,frequency,startdate,enddate,Fstartdate,Fenddate,Fcenddate,Fendsmpl,F)
% generate cells of date strings and vector of dates converted into decimal numbers; used for plots
% inputs:  - vector 'decimaldates1': dates converted into decimal values, for the sample period
%          - vector 'decimaldates2': dates converted into decimal values, for the sample+forecasts period
%          - cell 'stringdates1': date strings for the sample period
%          - cell 'stringdates2': date strings for the sample+forecasts period
%          - cell 'stringdates3': date strings for the forecast evaluation period (i.e. period for which forecast is estimated and actual data exists)
%          - integer 'Fstartlocation': position of the forecast start date in stringdates2
%          - integer 'Fendlocation': position of the forecast end date in stringdates2
% outputs: - cell 'names': cell containing the excel spreadsheet labels (names and dates)
%          - integer 'lags': number of lags included in the model
%          - integer 'frequency': frequency of the data set
%          - string 'startdate': start date of the sample
%          - string 'enddate': end date of the sample
%          - string 'Fstartdate': start date of the forecasts
%          - string 'Fenddate': end date of the forecasts
%          - string 'Fcenddate': end date of the forecat evaluation (i.e. period for which forecast is estimated and actual data exists)
%          - integer 'Fendsmpl': 0-1 value to determine if forecasts must start after the final sample period
%          - integer 'F': 0-1 value to determine if forecasts must be estimated

% preliminary tasks
% define date strings
datestrings=names(2:end,1);

if favar.FAVAR==1 % in case we transform the data to first or second differences, we have a different startlocation
    if favar.transformation==1
        datestrings=names(1+favar.informationstartlocation:favar.informationendlocation_sub,:); %first row are labels
    end
end

% PHASE 1: CREATION OF DECIMAL DATES AND DATE STRINGS FOR THE ESTIMATION SAMPLE (DECIMALDATES1, STRINGDATES1)

% deal first with the data if it is yearly
if frequency==1
    % identify the number of years covered by the sample
    startyear=str2double(startdate(1,1:end-1));
    endyear=str2double(enddate(1,1:end-1));
    % create the decimal value vector
    decimaldates1=(startyear:endyear)';
    % convert into strings
    stringdates1=cellfun(@num2str,num2cell(decimaldates1),'UniformOutput',0);
    % add a 'y' character at the end
    for ii=1:size(stringdates1,1)
        stringdates1{ii,1}=[stringdates1{ii,1} 'y'];
    end
    % trim a number of initial periods equal to the number of lags, as these periods will be used to create initial conditions
    decimaldates1=decimaldates1(lags+1:end,:);
    stringdates1=stringdates1(lags+1:end,:);

    % deal now with the data if it is quarterly
elseif frequency==2
    % first identify the year and quarter of the initial date
    startyear=str2double(startdate(1,1:4));
    startquarter=str2double(startdate(1,6));
    % proceed similarly for the year and quarter of the final date
    endyear=str2double(enddate(1,1:4));
    endquarter=str2double(enddate(1,6));
    % initiate the decimal value vector and the string cell
    decimaldates1=[];
    stringdates1={};
    % create the decimal vector from start date up to the penultimate year
    year=startyear;
    quarter=startquarter;
    while year<=endyear-1
        while quarter<=4
            decimaldates1=[decimaldates1;year+(quarter-1)/4];
            temp=[num2str(year) 'q' num2str(quarter)];
            stringdates1{end+1,1}=temp;
            quarter=quarter+1;
        end
        quarter=1;
        year=year+1;
    end
    % complete with final year
    while quarter<=endquarter
        decimaldates1=[decimaldates1;year+(quarter-1)/4];
        temp=[num2str(year) 'q' num2str(quarter)];
        stringdates1{end+1,1}=temp;
        quarter=quarter+1;
    end
    % finally, trim a number of initial periods equal to the number of lags, as these periods will be used to create initial conditions
    decimaldates1=decimaldates1(lags+1:end,:);
    stringdates1=stringdates1(lags+1:end,:);

    % deal now with the data if it is monthly
elseif frequency==3
    % first identify the year and month of the initial date
    temp=startdate;
    temp(1,5)=' ';
    [startyear,startmonth]=strtok(temp);
    startyear=str2double(startyear);
    startmonth=str2double(startmonth);
    % proceed similarly for the year and month of the final date
    temp=enddate;
    temp(1,5)=' ';
    [endyear,endmonth]=strtok(temp);
    endyear=str2double(endyear);
    endmonth=str2double(endmonth);
    % initiate the decimal value vector and the string cell
    decimaldates1=[];
    stringdates1={};
    % create the decimal vector from start date up to the penultimate year
    year=startyear;
    month=startmonth;
    while year<=endyear-1
        while month<=12
            decimaldates1=[decimaldates1;year+(month-1)/12];
            temp=[num2str(year) 'm' num2str(month)];
            stringdates1{end+1,1}=temp;
            month=month+1;
        end
        month=1;
        year=year+1;
    end
    % complete with final year
    while month<=endmonth
        decimaldates1=[decimaldates1;year+(month-1)/12];
        temp=[num2str(year) 'm' num2str(month)];
        stringdates1{end+1,1}=temp;
        month=month+1;
    end
    % finally, trim a number of initial periods equal to the number of lags, as these periods will be used to create initial conditions
    decimaldates1=decimaldates1(lags+1:end,:);
    stringdates1=stringdates1(lags+1:end,:);

    % deal now with the data if it is weekly
elseif frequency==4
    % first identify the first and last years of the sample
    startyear=str2double(startdate(1,1:4));
    endyear=str2double(enddate(1,1:4));
    % calculate the total number of years over the period
    numyear=endyear-startyear+1;
    % create a vector of date strings for the sample period
    start=find(strcmp(datestrings,startdate));
    finish=find(strcmp(datestrings,enddate));
    stringdates1=datestrings(start:finish,1);
    % for each year in the sample, identify the earliest and latest week
    for ii=startyear:endyear
        % first identify which periods in the sample correspond to this year
        periods=strfind(stringdates1,num2str(ii));
        for jj=1:size(periods,1)
            if isempty(periods{jj,1})
                periods{jj,1}=0;
            end
        end
        periods=cell2mat(periods);
        % identify the position of the first week and the last week of this year
        first=min(find(periods==1));
        last=max(find(periods==1));
        % identify the week number to which these positions correspond
        temp=char(stringdates1(first,1));
        temp(1,5)=' ';
        [~,weeknum]=strtok(temp);
        week(ii-startyear+1,1)=str2double(weeknum);
        temp=char(stringdates1(last,1));
        temp(1,5)=' ';
        [~,weeknum]=strtok(temp);
        week(ii-startyear+1,2)=str2double(weeknum);
    end
    % now create the vector of decimal data
    decimaldates1=[];
    % compute until penultimate year
    for ii=1:numyear-1
        for jj=week(ii,1):week(ii,2)
            decimaldates1=[decimaldates1;(startyear+ii-1)+((jj-1)/week(ii,2))];
        end
    end
    % compute for last year (assuming a year of 52 weeks if the number of weeks is shorter)
    if week(numyear,2)<52
        total=52;
    else
        total=week(numyear,2);
    end
    for jj=week(numyear,1):week(numyear,2)
        decimaldates1=[decimaldates1;endyear+((jj-1)/total)];
    end
    % finally, trim a number of initial periods equal to the number of lags, as these periods will be used to create initial conditions
    decimaldates1=decimaldates1(lags+1:end,:);
    stringdates1=stringdates1(lags+1:end,:);

    % deal now with the data if it is daily
elseif frequency==5
    % first identify the first and last years of the sample
    startyear=str2double(startdate(1,1:4));
    endyear=str2double(enddate(1,1:4));
    % calculate the total number of years over the period
    numyear=endyear-startyear+1;
    % create a vector of date strings for the sample period
    start=find(strcmp(datestrings,startdate));
    finish=find(strcmp(datestrings,enddate));
    stringdates1=datestrings(start:finish,1);
    % for each year in the sample, identify the earliest and latest day
    for ii=startyear:endyear
        % first identify which periods in the sample correspond to this year
        periods=strfind(stringdates1,num2str(ii));
        for jj=1:size(periods,1)
            if isempty(periods{jj,1})
                periods{jj,1}=0;
            end
        end
        periods=cell2mat(periods);
        % identify the position of the first day and the last day of this year
        first=min(find(periods==1));
        last=max(find(periods==1));
        % indentify the day number to which these positions correspond
        temp=char(stringdates1(first,1));
        temp(1,5)=' ';
        [~,daynum]=strtok(temp);
        day(ii-startyear+1,1)=str2double(daynum);
        temp=char(stringdates1(last,1));
        temp(1,5)=' ';
        [~,daynum]=strtok(temp);
        day(ii-startyear+1,2)=str2double(daynum);
    end
    % now create the vector of decimal data
    decimaldates1=[];
    % compute until penultimate year
    for ii=1:numyear-1
        for jj=day(ii,1):day(ii,2)
            decimaldates1=[decimaldates1;(startyear+ii-1)+((jj-1)/day(ii,2))];
        end
    end
    % compute for last year (assuming a working year of 5 days a week, i.e. 261 opening days a year, if the total number of day is shorter)
    if day(numyear,2)<261
        total=261;
    else
        total=day(numyear,2);
    end
    for jj=day(numyear,1):day(numyear,2)
        decimaldates1=[decimaldates1;endyear+((jj-1)/total)];
    end
    % finally, trim a number of initial periods equal to the number of lags, as these periods will be used to create initial conditions
    decimaldates1=decimaldates1(lags+1:end,:);
    stringdates1=stringdates1(lags+1:end,:);

    % finally, if the data is undated
elseif frequency==6
    % identify the number of periods covered by the sample
    startperiod=str2double(startdate(1,1:end-1));
    endperiod=str2double(enddate(1,1:end-1));
    % create the decimal value vector
    decimaldates1=(startperiod:endperiod)';
    % convert into strings
    stringdates1=cellfun(@num2str,num2cell(decimaldates1),'UniformOutput',0);
    % add a 'u' character at the end
    for ii=1:size(stringdates1,1)
        stringdates1{ii,1}=[stringdates1{ii,1} 'u'];
    end
    % trim a number of initial periods equal to the number of lags, as these periods will be used to create initial conditions
    decimaldates1=decimaldates1(lags+1:end,:);
    stringdates1=stringdates1(lags+1:end,:);
end

% from now on, the code applies only if forecasts have been selected
% if not, simply return empty matrices

if F==0 && CF==0
    decimaldates2=[];
    stringdates2=[];
    Fstartlocation=[];
    Fendlocation=[];
    stringdates3=[];
else

    % PHASE 2: CREATION OF DECIMAL DATES FOR THE FORECAST PERIOD (DECIMALDATES2)

    % determine first whether the final forecast period is included or not in the sample
    % first deal with data if it is yearly
    if frequency==1
        % simply check wether the sample end year is anterior to the forecast end year
        if str2double(enddate(1,1:end-1))<str2double(Fenddate(1,1:end-1))
            included=0;
        else
            included=1;
        end
        % if data is quarterly
    elseif frequency==2
        % first identify the year and quarter of the final sample date
        smplendyear=str2double(enddate(1,1:4));
        smplendquarter=str2double(enddate(1,6));
        % similarly, identify the year and quarter of the final forecast date
        Fendyear=str2double(Fenddate(1,1:4));
        Fendquarter=str2double(Fenddate(1,6));
        if smplendyear<Fendyear
            included=0;
        elseif smplendyear==Fendyear && smplendquarter<Fendquarter
            included=0;
        else
            included=1;
        end
        % if data is monthly
    elseif frequency==3
        % first identify the year and month of the final sample date
        temp=enddate;
        temp(1,5)=' ';
        [smplendyear,smplendmonth]=strtok(temp);
        smplendyear=str2double(smplendyear);
        smplendmonth=str2double(smplendmonth);
        % identify the year and month of the final forecast date
        temp=Fenddate;
        temp(1,5)=' ';
        [Fendyear,Fendmonth]=strtok(temp);
        Fendyear=str2double(Fendyear);
        Fendmonth=str2double(Fendmonth);
        if smplendyear<Fendyear
            included=0;
        elseif smplendyear==Fendyear && smplendmonth<Fendmonth
            included=0;
        else
            included=1;
        end
        % if data is weekly
    elseif frequency==4
        % first identify the year and week of the final sample date
        temp=enddate;
        temp(1,5)=' ';
        [smplendyear,smplendweek]=strtok(temp);
        smplendyear=str2double(smplendyear);
        smplendweek=str2double(smplendweek);
        % identify the year and week of the final forecast date
        temp=Fenddate;
        temp(1,5)=' ';
        [Fendyear,Fendweek]=strtok(temp);
        Fendyear=str2double(Fendyear);
        Fendweek=str2double(Fendweek);
        if smplendyear<Fendyear
            included=0;
        elseif smplendyear==Fendyear && smplendweek<Fendweek
            included=0;
        else
            included=1;
        end
        % if data is daily
    elseif frequency==5
        % first identify the year and day of the final sample date
        temp=enddate;
        temp(1,5)=' ';
        [smplendyear,smplendday]=strtok(temp);
        smplendyear=str2double(smplendyear);
        smplendday=str2double(smplendday);
        % identify the year and day of the final forecast date
        temp=Fenddate;
        temp(1,5)=' ';
        [Fendyear,Fendday]=strtok(temp);
        Fendyear=str2double(Fendyear);
        Fendday=str2double(Fendday);
        if smplendyear<Fendyear
            included=0;
        elseif smplendyear==Fendyear && smplendday<Fendday
            included=0;
        else
            included=1;
        end
        % finally, if data is undated
    elseif frequency==6
        % simply check wether the sample end period is anterior to the forecast end period
        if str2double(enddate(1,1:end-1))<str2double(Fenddate(1,1:end-1))
            included=0;
        else
            included=1;
        end
    end

    % if the final forecast period is included in the sample, simply define decimaldates2 as decimaldates1
    if included==1
        decimaldates2=decimaldates1;
        stringdates2=stringdates1;

        % if the final forecast period lies beyond the end of the sample, define a new vector of decimal dates running from the sample start until the last forecast period
    elseif included==0
        % deal first with the data if it is yearly
        if frequency==1
            % simply complete decimaldates1 up to the final forecast period
            decimaldates2=[decimaldates1;(decimaldates1(end,1)+1:str2double(Fenddate(1,1:end-1)))'];
            stringdates2=cellfun(@num2str,num2cell(decimaldates2),'UniformOutput',0);
            % add a 'y' character at the end
            for ii=1:size(stringdates2,1)
                stringdates2{ii,1}=[stringdates2{ii,1} 'y'];
            end

            % if the data is quarterly
        elseif frequency==2
            % initiate the series
            decimaldates2=decimaldates1;
            stringdates2=stringdates1;
            % identify the year and quarter of the last sample period
            endyear=str2double(enddate(1,1:4));
            endquarter=str2double(enddate(1,6));
            % identify the year and quarter of the last forecast period
            Fendyear=str2double(Fenddate(1,1:4));
            Fendquarter=str2double(Fenddate(1,6));
            % advance sample end by one period
            year=endyear;
            quarter=endquarter;
            if quarter<4
                quarter=quarter+1;
            elseif quarter==4
                year=year+1;
                quarter=1;
            end
            while year<=Fendyear-1
                while quarter<=4
                    decimaldates2=[decimaldates2;year+(quarter-1)/4];
                    temp=[num2str(year) 'q' num2str(quarter)];
                    stringdates2{end+1,1}=temp;
                    quarter=quarter+1;
                end
                quarter=1;
                year=year+1;
            end
            % complete with final year
            while quarter<=Fendquarter
                decimaldates2=[decimaldates2;year+(quarter-1)/4];
                temp=[num2str(year) 'q' num2str(quarter)];
                stringdates2{end+1,1}=temp;
                quarter=quarter+1;
            end

            % if the data is monthly
        elseif frequency==3
            % initiate the series
            decimaldates2=decimaldates1;
            stringdates2=stringdates1;
            % identify the year and month of the last sample period
            temp=enddate;
            temp(1,5)=' ';
            [endyear,endmonth]=strtok(temp);
            endyear=str2double(endyear);
            endmonth=str2double(endmonth);
            % identify the year and month of the last forecast period
            temp=Fenddate;
            temp(1,5)=' ';
            [Fendyear,Fendmonth]=strtok(temp);
            Fendyear=str2double(Fendyear);
            Fendmonth=str2double(Fendmonth);
            % advance sample end by one period
            year=endyear;
            month=endmonth;
            if month<12
                month=month+1;
            elseif month==12
                year=year+1;
                month=1;
            end
            while year<=Fendyear-1
                while month<=12
                    decimaldates2=[decimaldates2;year+(month-1)/12];
                    temp=[num2str(year) 'm' num2str(month)];
                    stringdates2{end+1,1}=temp;
                    month=month+1;
                end
                month=1;
                year=year+1;
            end
            % complete with final year
            while month<=Fendmonth
                decimaldates2=[decimaldates2;year+(month-1)/12];
                temp=[num2str(year) 'm' num2str(month)];
                stringdates2{end+1,1}=temp;
                month=month+1;
            end

            % if the data is weekly
        elseif frequency==4
            % first identify the year and week of the final forecast date
            temp=Fenddate;
            temp(1,5)=' ';
            [Fendyear,Fendweek]=strtok(temp);
            Fendyear=str2double(Fendyear);
            Fendweek=str2double(Fendweek);
            % identify the year and week of the final data set date
            temp=datestrings{end,1};
            temp(1,5)=' ';
            [dataendyear,dataendweek]=strtok(temp);
            dataendyear=str2double(dataendyear);
            dataendweek=str2double(dataendweek);
            % now, there are two possibilities: either the forecast period is entirely included in the dataset, or it goes beyond the dataset
            % if entirely included in the dataset
            if Fendyear<dataendyear || (Fendyear==dataendyear && Fendweek<dataendweek)
                % then just copy the date strings from sample start to forecast end
                % find position for sample start
                start=find(cellfun(@isempty,strfind(datestrings,startdate))==0);
                % find position for forecast end
                finish=find(cellfun(@isempty,strfind(datestrings,Fenddate))==0);
                % create the strings
                stringdates2=datestrings(start:finish,1);
                % if the forecast periods goes beyond the data set
            else
                % find position for sample start
                start=find(cellfun(@isempty,strfind(datestrings,startdate))==0);
                % copy the dataset from sample start to its end
                stringdates2=datestrings(start:end,1);
                % advance the end of the dataset by one period
                if dataendweek<52
                    year=dataendyear;
                    week=dataendweek+1;
                else
                    year=dataendyear+1;
                    week=1;
                end
                % now complete stringdates2
                for ii=dataendyear:Fendyear-1
                    while week<=52
                        stringdates2{end+1,1}=[num2str(year) 'w' num2str(week)];
                        week=week+1;
                    end
                    week=1;
                    year=year+1;
                end
                % complete for the final year
                while week<=Fendweek
                    stringdates2{end+1,1}=[num2str(year) 'w' num2str(week)];
                    week=week+1;
                end
            end
            % finally, trim a number of initial conditions equal to the number of lags
            stringdates2=stringdates2(lags+1:end,1);
            % now that the strings are obtained, use them to generate the decimal dates
            % for each year in stringdates2, identify the number of weeks
            maxweek=zeros(size(stringdates2,1),1);
            for ii=str2double(startdate(1,1:4)):Fendyear-1
                periods=strfind(stringdates2,num2str(ii));
                for jj=1:size(periods,1)
                    if isempty(periods{jj,1})
                        periods{jj,1}=0;
                    end
                end
                periods=cell2mat(periods);
                location=max(find(periods==1));
                temp=char(stringdates2(location,1));
                temp(1,5)=' ';
                [~,weeknum]=strtok(temp);
                maxweek=maxweek+periods*str2double(weeknum);
            end
            % complete for the final year
            periods=strfind(stringdates2,num2str(Fendyear));
            for jj=1:size(periods,1)
                if isempty(periods{jj,1})
                    periods{jj,1}=0;
                end
            end
            periods=cell2mat(periods);
            maxweek=maxweek+periods*max(52,Fendweek);
            % finally, compute the decimal dates
            for ii=1:size(stringdates2,1)
                temp=stringdates2{ii,1};
                temp(1,5)=' ';
                [~,weeknum]=strtok(temp);
                decimaldates2(ii,1)=str2double(temp(1,1:4))+str2double(weeknum)/maxweek(ii,1);
            end

            % if the data is daily
        elseif frequency==5
            % first identify the year and day of the final forecast date
            temp=Fenddate;
            temp(1,5)=' ';
            [Fendyear,Fendday]=strtok(temp);
            Fendyear=str2double(Fendyear);
            Fendday=str2double(Fendday);
            % identify the year and day of the final data set date
            temp=datestrings{end,1};
            temp(1,5)=' ';
            [dataendyear,dataendday]=strtok(temp);
            dataendyear=str2double(dataendyear);
            dataendday=str2double(dataendday);
            % now, there are two possibilities: either the forecast period is entirely included in the dataset, or it goes beyond the dataset
            % if entirely included in the dataset
            if Fendyear<dataendyear || (Fendyear==dataendyear && Fendday<dataendday)
                % then just copy the date strings from sample start to forecast end
                % find position for sample start
                start=find(cellfun(@isempty,strfind(datestrings,startdate))==0);
                % find position for forecast end
                finish=find(cellfun(@isempty,strfind(datestrings,Fenddate))==0);
                % create the strings
                stringdates2=datestrings(start:finish,1);
                % if the forecast periods goes beyond the data set
            else
                % find position for sample start
                start=find(cellfun(@isempty,strfind(datestrings,startdate))==0);
                % copy the dataset from sample start to its end
                stringdates2=datestrings(start:end,1);
                % advance the day of the dataset by one period
                if dataendday<261
                    year=dataendyear;
                    day=dataendday+1;
                else
                    year=dataendyear+1;
                    day=1;
                end
                % now complete stringdates2
                for ii=dataendyear:Fendyear-1
                    while day<=261
                        stringdates2{end+1,1}=[num2str(year) 'd' num2str(day)];
                        day=day+1;
                    end
                    day=1;
                    year=year+1;
                end
                % complete for the final year
                while day<=Fendday
                    stringdates2{end+1,1}=[num2str(year) 'd' num2str(day)];
                    day=day+1;
                end
            end
            % finally, trim a number of initial conditions equal to the number of lags
            stringdates2=stringdates2(lags+1:end,1);
            % now that the strings are obtained, use them to generate the decimal dates
            % for each year in stringdates2, identify the number of day
            maxday=zeros(size(stringdates2,1),1);
            for ii=str2double(startdate(1,1:4)):Fendyear-1
                periods=strfind(stringdates2,num2str(ii));
                for jj=1:size(periods,1)
                    if isempty(periods{jj,1})
                        periods{jj,1}=0;
                    end
                end
                periods=cell2mat(periods);
                location=max(find(periods==1));
                temp=char(stringdates2(location,1));
                temp(1,5)=' ';
                [~,daynum]=strtok(temp);
                maxday=maxday+periods*str2double(daynum);
            end
            % complete for the final year
            periods=strfind(stringdates2,num2str(Fendyear));
            for jj=1:size(periods,1)
                if isempty(periods{jj,1})
                    periods{jj,1}=0;
                end
            end
            periods=cell2mat(periods);
            maxday=maxday+periods*max(261,Fendday);
            % finally, compute the decimal dates
            for ii=1:size(stringdates2,1)
                temp=stringdates2{ii,1};
                temp(1,5)=' ';
                [~,daynum]=strtok(temp);
                decimaldates2(ii,1)=str2double(temp(1,1:4))+str2double(daynum)/maxday(ii,1);
            end

            % finally, if the data is undated
        elseif frequency==6
            % simply complete decimaldates1 up to the final forecast period
            decimaldates2=[decimaldates1;(decimaldates1(end,1)+1:str2double(Fenddate(1,1:end-1)))'];
            stringdates2=cellfun(@num2str,num2cell(decimaldates2),'UniformOutput',0);
            % add a 'u' character at the end
            for ii=1:size(stringdates2,1)
                stringdates2{ii,1}=[stringdates2{ii,1} 'u'];
            end
        end
    end

    % PHASE 3: IDENTIFICATION OF THE POSITION OF THE INITIAL FORECAST PERIOD (IN TERMS OF THE VECTOR DECIMALDATES2)

    % identify the position of the start period for the forecasts
    % if the start period has been selected as the first period after the sample end, identifies it directly
    if Fendsmpl==1
        Fstartlocation=size(stringdates1,1)+1;
        % if the start period has not been selected as the first period after the sample end, it must be within the sample: look for it
    elseif Fendsmpl==0
        Fstartlocation=find(strcmp(stringdates1,Fstartdate));
    end

    % identify the position of the end period for the forecasts
    % if the end period is included in the sample, it can be identified directly
    if included==1
        Fendlocation=find(strcmp(stringdates1,Fenddate));
        % if the end period is beyond the sample, then its position is simply the last period of datestrings2
    elseif included==0
        Fendlocation=size(stringdates2,1);
    end

    % PHASE 3: CREATION OF DATE STRINGS FOR THE COMMON PERIOD

    % these date strings will be used only at one point: for the display of forecast evaluation
    % to identify the strings, simply copy from datestrings2 the dates over the common forecast period
    startfcperiod=find(strcmp(stringdates2,Fstartdate));
    endfcperiod=find(strcmp(stringdates2,Fcenddate));
    stringdates3=stringdates2(startfcperiod:endfcperiod,1);

end

author: european-central-bank

repo: BEAR-toolbox


text:
function [X,Y,Z] = cylinderWeights(r1,r2)
%CYLINDERWEIGHTS creates the innermost blue bars with heavy looking cylinders
%representing attatched weights.
t=0:0.01:1;
x=t;

for i=1:numel(t)
    if ((i>10 && i<20) || (i>80 && i<90))
        x(i)=r2;
    else
        x(i)=r1;
    end
end
[X,Y,Z]=cylinder(x,30);


text:
% Jiao Xianjun (putaoshu@msn.com; putaoshu@gmail.com)
% convert a sequence to a mat. each column shift previos to left by 1
% for example input sequence 1 2 3 4 5 6 7, and len_body=4
% output matrix:
% 1 2 3 4
% 2 3 4 5
% 3 4 5 6
% 4 5 6 7

% A script of project: https://github.com/JiaoXianjun/multi-rtl-sdr-calibration

function r = lin2col_shift_mat(s, len_body)

if length(s) < len_body
    disp('Length of input is too short!');
    r = -1;
    return;
end

len_tail = length(s) - len_body;

r = toeplitz(s, [s(1) zeros(1, len_tail)]);
r = r((len_tail+1):end, end:-1:1);

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function [P, beta] = x2p(X, u, tol)
%X2P Identifies appropriate sigma's to get kk NNs up to some tolerance
%
%   [P, beta] = x2p(xx, kk, tol)
%
% Identifies the required precision (= 1 / variance^2) to obtain a Gaussian
% kernel with a certain uncertainty for every datapoint. The desired
% uncertainty can be specified through the perplexity u (default = 15). The
% desired perplexity is obtained up to some tolerance that can be specified
% by tol (default = 1e-4).
% The function returns the final Gaussian kernel in P, as well as the
% employed precisions per instance in beta.
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    if ~exist('u', 'var') || isempty(u)
        u = 15;
    end
    if ~exist('tol', 'var') || isempty(tol)
        tol = 1e-4;
    end

    % Initialize some variables
    n = size(X, 1);                     % number of instances
    P = zeros(n, n);                    % empty probability matrix
    beta = ones(n, 1);                  % empty precision vector
    logU = log(u);                      % log of perplexity (= entropy)

    % Compute pairwise distances
%     disp('Computing pairwise distances...');
    D = squareform(pdist(X, 'euclidean') .^ 2);

    % Run over all datapoints
%     disp('Computing P-values...');
    for i=1:n

%         if ~rem(i, 500)
%             disp(['Computed P-values ' num2str(i) ' of ' num2str(n) ' datapoints...']);
%         end

        % Set minimum and maximum values for precision
        betamin = -Inf;
        betamax = Inf;

        % Compute the Gaussian kernel and entropy for the current precision
        Di = D(i, [1:i-1 i+1:end]);
        [H, thisP] = Hbeta(Di, beta(i));

        % Evaluate whether the perplexity is within tolerance
        Hdiff = H - logU;
        tries = 0;
        while abs(Hdiff) > tol && tries < 50

            % If not, increase or decrease precision
            if Hdiff > 0
                betamin = beta(i);
                if isinf(betamax)
                    beta(i) = beta(i) * 2;
                else
                    beta(i) = (beta(i) + betamax) / 2;
                end
            else
                betamax = beta(i);
                if isinf(betamin)
                    beta(i) = beta(i) / 2;
                else
                    beta(i) = (beta(i) + betamin) / 2;
                end
            end

            % Recompute the values
            [H, thisP] = Hbeta(Di, beta(i));
            Hdiff = H - logU;
            tries = tries + 1;
        end

        % Set the final row of P
        P(i, [1:i - 1, i + 1:end]) = thisP;
    end
%     disp(['Mean value of sigma: ' num2str(mean(sqrt(1 ./ beta)))]);
%     disp(['Minimum value of sigma: ' num2str(min(sqrt(1 ./ beta)))]);
%     disp(['Maximum value of sigma: ' num2str(max(sqrt(1 ./ beta)))]);
end

% Function that computes the Gaussian kernel values given a vector of
% squared Euclidean distances, and the precision of the Gaussian kernel.
% The function also computes the perplexity of the distribution.
function [H, P] = Hbeta(D, beta)
    P = exp(-D * beta);
    sumP = sum(P);
    H = log(sumP) + beta * sum(D .* P) / sumP;
    P = P / sumP;
end

repo: cluster


text:
function plist = partitions(total_sum,candidate_set,max_count,fixed_count)
% extracts the list of all partitions of a number as integer sums of a list of candidates
% usage: plist = partitions(total_sum,candidate_set)
% usage: plist = partitions(total_sum,candidate_set,max_count,fixed_count)
%
% PARTITIONS solves the money changing problem. E.g.,
% how can you make change for one dollar given coins
% of a given set of denominations. A good reference on
% the general problem is found here:
%
% http://en.wikipedia.org/wiki/Integer_partition
%
% PARTITIONS uses a recursive strategy to enumerate all
% possible partitions of the total_sum. This may be
% highly intensive for large sums or large sets of
% candidates.
%
% arguments: (input)
%  total_sum - scalar positive integer (to be partitioned)
%
%              BEWARE! a large total_sum can easily cause
%              stack problems. For example, the number of
%              partitions of 40 is 37338, a set that took 24
%              seconds to completely enumerate on my cpu.
%
%  candidate_set - (OPTIONAL) vector of (distinct) candidate
%              positive integers for the partitions.
%
%              Efficiency considerations force me to require
%              that the candidates be sorted in non-decreasing
%              order. An error is produced otherwise.
%
%              DEFAULT: candidate_set = 1:total_sum
%
%              BEWARE! large candidate sets can easily cause
%              stack problems
%
%  max_count - (OPTIONAL) the maximum quantity of any
%              candidate in the final sum.
%
%              max_count must be either a vector of the
%              same length as candidate_set, or a scalar
%              that applies to all elements in that set.
%
%              DEFAULT = floor(total_sum./candidate_set)
%
%  fixed_count - (OPTIONAL) Allows you to specify a fixed
%              number of terms in the partitioned sum.
%
%              fixed_count must be a positive integer if
%              supplied.
%
%              DEFAULT = []
%
% arguments: (output)
%  plist - array of partitions of total_sum. This is a list
%              of the quantity of each element such that
%              plist*candidate_set(:) yields total_sum
%
%
% example: Write 9 as an integer combination of the set [1 2 4 7]
%
%  partitions(9,[1 2 4 7])
%
% ans =
%    9     0     0     0
%    7     1     0     0
%    5     2     0     0
%    3     3     0     0
%    1     4     0     0
%    5     0     1     0
%    3     1     1     0
%    1     2     1     0
%    1     0     2     0
%    2     0     0     1
%    0     1     0     1
%
% Thus, we can write 9 = 9*1
% or 9 = 1*1 + 4*2
% or 9 = 1*2 + 1*7
% or any of 8 distinct other ways.
%
% There are 11 such ways to write 9 in terms of these
% candidates.
%
%
% example: Change a 1 dollar bill (100 cents) as an integer
%  combination of the set [1 5 10 25 50], using no more than
%  4 of any one coin denomination. Note that no pennies will
%  be allowed by the maximum constraint.
%
%  partitions(100,[1 5 10 25 50],4)
%
% ans =
%    0     4     3     2     0
%    0     2     4     2     0
%    0     3     1     3     0
%    0     1     2     3     0
%    0     0     0     4     0
%    0     4     3     0     1
%    0     2     4     0     1
%    0     3     1     1     1
%    0     1     2     1     1
%    0     0     0     2     1
%    0     0     0     0     2
%
% example: Write 13 as an integer combination of the set [2 4 6 8 10 12]
%  (Note that no such combination exists.)
%
%  partitions(13,[2 4 6 8 10 12])
%
% ans =
%   Empty matrix: 0-by-6
%
%
% example: find all possible ways to write 100 as a sum of EXACTLY 4
%  squares of the integers 1:9.
%
%  partitions(100,(1:9).^2,[],4)
% ans =
%   0    0    0    0    4    0    0    0    0
%   1    0    0    0    2    0    1    0    0
%   2    0    0    0    0    0    2    0    0
%   0    1    0    2    0    0    0    1    0
%   1    0    2    0    0    0    0    0    1
%
%
% Author: John D'Errico
% e-mail: woodchips@rochester.rr.com
% Release: 2
% Release date: 7/15/08

% default for candidate_set
if (nargin<2) || isempty(candidate_set)
  candidate_set = 1:total_sum;
end

% how many candidates are there
n = length(candidate_set);

% error checks
if any(candidate_set<=0)
  error('All members of candidate_set must be > 0')
end
% candidates must be sorted in increasng order
if any(diff(candidate_set)<0)
  error('Efficiency requires that candidate_set be sorted')
end

% check for a max_count. do we supply a default?
if (nargin<3) || isempty(max_count)
  % how high do we need look?
  max_count = floor(total_sum./candidate_set);
elseif length(max_count)==1
  % if a scalar was provided, then turn it into a vector
  max_count = repmat(max_count,1,n);
end

% check for a fixed_count
if (nargin<4) || isempty(fixed_count)
  fixed_count = [];
elseif (fixed_count<0) || (fixed_count~=round(fixed_count))
  error('fixed_count must be a positive integer if supplied')
end

% check for degenerate cases
if isempty(fixed_count)
  if total_sum == 0
    plist = zeros(1,n);
    return
  elseif (n == 0)
    plist = [];
    return
  elseif (n == 1)
    % only one element in the set. can we form
    % total_sum from it as an integer multiple?
    p = total_sum/candidate_set;
    if (p==fix(p)) && (p<=max_count)
      plist = p;
    else
      plist = [];
    end
    return
 end
else
  % there was a fixed_count supplied
  if (total_sum == 0) && (fixed_count == 0)
    plist = zeros(1,n);
    return
  elseif (n == 0) || (fixed_count <= 0)
    plist = [];
    return
  elseif (n==1)
    % there must be a non-zero fixed_count, since
    % we did not trip the last test. since there
    % is only one candidate in the set, will it work?
    if (fixed_count*candidate_set) == total_sum
      plist = fixed_count;
    else
      plist = [];
    end
    return
  end
end

% finally, we can do some work. start with the
% largest element and work backwards
m = max_count(end);
% do we need to back off on m?
c = candidate_set(end);
m = min([m,floor(total_sum/c),fixed_count]);

plist = zeros(0,n);
for i = 0:m
  temp = partitions(total_sum - i*c, ...
      candidate_set(1:(end-1)), ...
      max_count(1:(end-1)),fixed_count-i);
  plist = [plist;[temp,repmat(i,size(temp,1),1)]];  %#ok
end


text:
% Copyright (C) 2016, 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 speed_demo_with_movement_sinusoidal
close all

robot = load_robot('OMRON', 'ECOBRA600')

% example trajectories (10 seconds)
delta_t = 0.01; % s
t = 0:delta_t:20;
w1 = 0.5; %rad/s
w2 = 1.0; %rad/s
w3 = 1.5; %m/s

qt = [sin(w1*t); cos(w2*t);  0.1*(sin(w3*t)+0.3*cos(w3*t)); 0*t];

qds = [];
vs = [];

for i=1:length(t)-1
    % joint position
    q = qt(:, i);
    % joint speed
    qd = (qt(:, i+1) - qt(:, i))/delta_t;

    J = manipulator_jacobian(robot, q);
    % end effector's velocity
    v = J*qd;
    qds = [qds qd];
    vs = [vs v];
end
qds = [qds qd];
vs = [vs v];

figure, plot(t, qt)
title('Posición articular')
xlabel('tiempo (s)')
ylabel('q (rad)')
legend('q1 (rad)', 'q2 (rad)', 'q3 (m)')

figure, plot(t, qds)
title('Velocidad articular')
xlabel('tiempo (s)')
ylabel('qd (rad)')
legend('qd1 (rad/s)', 'qd2 (rad/s)', 'qd3 (m/s)')

figure, plot(t, vs)
title('Velocidad en el extremo')
xlabel('tiempo (s)')
ylabel('v (m/s)')
legend('vx m/s', 'vy m/s', 'vz m/s')

% animamos el resultado
%animate(robot, qt(:, 1:10:end))
% animamos el resultado con un vector de velocidad en el extremo del robot
for i=1:5:length(t)
    q = qt(:,i);
    v = vs(:,i);
    T = directkinematic(robot, q);

    %plot speed
    p0 = T(1:3,4);
    drawrobot3d(robot, q)
    draw_vector(v(1:3), p0, '       V', 3)
    pause(0.1)
end

repo: ARTE


text:
% Simple functino for visualizing the trajectory
% of a motor program M
function simpleViz(M)

    cell_traj = M.motor_warped;
    ns = numel(cell_traj);

    hold on
    for sid=1:ns
        stroke = cell_traj{sid};
        nsub = length(stroke);
        for b=1:nsub
            plot_traj(cell_traj{sid}{b});
        end
    end

    xlim([0 105]);
    ylim([-105 0]);
    box on
end

function plot_traj(T)
    plot(T(:,1),T(:,2),'k','LineWidth',2);
end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
function [ n_data, n, x, fx ] = p_polynomial_values ( n_data )

%*****************************************************************************80
%
%% P_POLYNOMIAL_VALUES returns values of the Legendre polynomials P(n,x).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, the order of the function.
%
%    Output, real X, the point where the function is evaluated.
%
%    Output, real FX, the value of the function.
%
  n_max = 22;

  fx_vec = [ ...
      0.1000000000000000E+01, ...
      0.2500000000000000E+00, ...
     -0.4062500000000000E+00, ...
     -0.3359375000000000E+00, ...
      0.1577148437500000E+00, ...
      0.3397216796875000E+00, ...
      0.2427673339843750E-01, ...
     -0.2799186706542969E+00, ...
     -0.1524540185928345E+00, ...
      0.1768244206905365E+00, ...
      0.2212002165615559E+00, ...
      0.0000000000000000E+00, ...
     -0.1475000000000000E+00, ...
     -0.2800000000000000E+00, ...
     -0.3825000000000000E+00, ...
     -0.4400000000000000E+00, ...
     -0.4375000000000000E+00, ...
     -0.3600000000000000E+00, ...
     -0.1925000000000000E+00, ...
      0.8000000000000000E-01, ...
      0.4725000000000000E+00, ...
      0.1000000000000000E+01 ];

  n_vec = [ ...
     0,  1,  2, ...
     3,  4,  5, ...
     6,  7,  8, ...
     9, 10,  3, ...
     3,  3,  3, ...
     3,  3,  3, ...
     3,  3,  3, ...
     3 ];

  x_vec = [ ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.00E+00, ...
     0.10E+00, ...
     0.20E+00, ...
     0.30E+00, ...
     0.40E+00, ...
     0.50E+00, ...
     0.60E+00, ...
     0.70E+00, ...
     0.80E+00, ...
     0.90E+00, ...
     1.00E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    x = 0.0;
    fx = 0.0;
  else
    n = n_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function result = blur(im,levels,filt)
% Blurs an image by blurring and subsampling repeatedly, followed by
% upsampling and blurring.
%
%      result=blur(im,[levels],[filt])
%
%  im     -  input image.
%  levels -  number of times to blur and subsample (default is 1).
%  filt   -  blurring 1d filter to be applied separably to the rows and
%            cols of im (default ='binom5').
%
% HISTORY:
%
% DJH '96
% update 12/97 to conform to Eero's updated pyrTools
%
% I don't think we should use this any more.  Why not just use Matlab's
% convolution or blurring or filtering? (BW)
% Because this is routine *much* faster than convolution. (RFD)
%
% 2007.05.03 RFD: fixed non-double class support.

% warning('This function calls very old mex files that might not work. Try imblur instead.');

if ~exist('levels','var'), levels=1; end
if ~exist('filt','var'),   filt = 'binom5'; end
if ischar(filt),           filt = namedFilter(filt); end

% ras 03/07: enforce double-precision data type
% (the compiled corrDn functions need double)
if ~isa(im, 'double'), type = class(im); im = double(im); end
tmp = blurDn(im,levels,filt);

% save upBlurDEBUG tmp levels filt
result = upBlur(tmp,levels,filt);

% Make sure its the same size as the input image
result = result((1:size(im,1)),(1:size(im,2)));

if exist('type', 'var')
	% we had a non-double matrix, convert back
	result = feval(type, result);
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft