text:
function [ ] = montyhall(n,d,o)
%montyhall(n,d,o) Simulates n rounds of the Monty Hall problem and
%                 gives win/loss ratio
%
%                n   : Number of rounds to simulate (Any positive integer)
%                d   : Number of doors (Positive integer >= 3)
%                o=0 : Show final win/loss score and estimated probability
%                o=1 : Show win/loss tally chart and estimated probability
%                o=2 : Show results for each round
%
% All parameters are optional, by default they will take the values:
%  n=1, d=3, o=0.
%
% The Monty Hall problem is a probability puzzle based on the American
% television game show Let's Make a Deal. The name comes from the show's
% host, Monty Hall. The problem is also called the Monty Hall paradox,
% as it is a veridical paradox in that the result appears absurd but is
% demonstrably true.
% http://en.wikipedia.org/wiki/Monty_Hall_problem
%
%
%------------------------%

win=0;lose=0;tally=[0;0];Monty=0;
%-- Check for nonexistent variables and set to defaults
if exist('n')~=1, n=1; end
if exist('d')~=1 || d<3 || d~=round(d), d=3; end
if exist('o')~=1, o=0; end
%--

if o==2
    for i = 1:n
        %-- Picks a random door for the car to be behind
        Car = randi(d)-1;
        %-- Player picks a random door
        Choose = randi(d)-1;
        %-- Monty recursively picks a random door until he has chosen
        %   the door that neither has the car behind it or has been
        %   chosen by the player
        while (Monty==Car || Monty==Choose), Monty=randi(d)-1; end
        %-- If you originally choose the car, no matter what door Monty
        %   opens, you will lose by switching.
        %   If you originally choose a goat, Monty has to open a door to
        %   reveal the other goat, and switching will switch to the car
        %   meaning you win.
        %   Simply put, if you originally choose the car, switching wins
        %   and if you don't, switching loses.
        %   This is the basis of the score tally system
        if Choose==Car, lose=lose+1; else win=win+1; end
        %-- Display round results after game has been played.
        disp(['Round ',num2str(i)]);
        disp('Chosen door:');disp([zeros(1,Choose),1,zeros(1,(d-1)-Choose)]);
        if Car<Choose, disp('Monty opens:');disp([ones(1,Car),0,ones(1,Choose-Car-1),0,ones(1,(d-1)-Choose)]); end
        if Car>Choose, disp('Monty opens:');disp([ones(1,Choose),0,ones(1,Car-Choose-1),0,ones(1,(d-1)-Car)]); end
        if Car==Choose
            if Car<Monty, disp('Monty opens:');disp([ones(1,Car),0,ones(1,Monty-Car-1),0,ones(1,(d-1)-Monty)]); else disp('Monty opens:');disp([ones(1,Monty),0,ones(1,Car-Monty-1),0,ones(1,(d-1)-Car)]); end
        end
        disp('Car is behind:');disp([zeros(1,Car),1,zeros(1,(d-1)-Car)]);
        if Choose==Car, disp('Switching loses'); else disp('Switching wins'); end
        disp(' ');disp('   Win | Lose');disp([win,lose]);
        disp('---------------------');
    end
else
    for i = 1:n
        %-- When detail~=2, the next few lines simulate the game
        %   instead. This is a lot simpler as all unnecessary commands are
        %   removed. In the original block of code (where detail=2), the
        %   code can be summarised as;
        %   Pick 2 random numbers (Car and Choose), pick a random number
        %   not equal to those numbers (Monty) and finally check if the two
        %   original numbers are equal.
        %   The number that is chosen for Monty does not affect whether a
        %   win or a loss is tallied, so it is not needed for the simulation.
        Car = randi(d)-1;
        Choose = randi(d)-1;
        %-- Defining 'Monty' is not needed to decide a win or a loss.
        %   while (Monty==Car || Monty==Choose), Monty=randi(d)-1; end
        if Choose==Car, lose=lose+1; else win=win+1; end
        tally=[tally,[win;lose]];
    end
end
if o==1
    disp(['W: ',num2str(tally(1,1:n+1))]);
    disp(['L: ',num2str(tally(2,1:n+1))]);
else
    disp(['Won  ',num2str(win)]);
    disp(['Lost ',num2str(lose)]);
end
if n>=10, disp(['Estimated Chance Of Winning After Switch: ',num2str(win/n)]); end
end


text:
function [dpixc, dpixc_ind, blinkmat] = generatedataMultiPSF(size_vec, x_vec,y_vec,intensity_vec,psf_tmp,offset,Nt, probtrans)
% generate points

if size(intensity_vec,2)>size(intensity_vec,1); intensity_vec = intensity_vec'; end
if size(x_vec,2)>size(x_vec,1); x_vec = x_vec'; end
if size(y_vec,2)>size(y_vec,1); y_vec = y_vec'; end

nx = size_vec(2)-size_vec(1);
ny = size_vec(4)-size_vec(3);
N = length(x_vec); %number of points
center = [nx ny]/2;
center_round = round(center);

%
% intensity_mat = repmat(intensity_vec, 1,Nt);
% changemat = rand(N,Nt)<probtrans;
% statemat = mod(cumsum(changemat,2),2);
% initvec = rand(N,1)>0.5; %initial state of the blinkmat
% ivt = ~(initvec == statemat(:,1));
% statematinit = mod(statemat+repmat(ivt,1,Nt),2);
%
% % blinkmat_equal = rand(N, Nt);
% % blinkmat = blinkmat_equal .* intensity_mat; %different intensities...
% blinkmat = statematinit .* intensity_mat; %different intensities...
blinkmat = blinkmat_markov(N,Nt, intensity_vec, probtrans);
center_im = pixelize(center_round, 1, size_vec, nx, ny, [],0);
dpixc_ind = newimar(N);

shift_vec = [x_vec-center(1), y_vec-center(2)];
for ii=1:N
    sp = size(psf_tmp{ii});
    psf_tmp2 = double(psf_tmp{ii}/sum(psf_tmp{ii}(:))); %normalize to 1
    psf_tmp3 = padimage(psf_tmp2,fliplr([nx, ny])+2*sp);
    psf = psf_tmp3;
    psf_im = clip(dip_image(conv2(center_im,psf,'same')), 0, Inf);
    dpixc_ind{ii} = clip(shift(psf_im, shift_vec(ii,:)), 0, Inf);
end

dpixc_nonoise = array2im(dpixc_ind'*blinkmat + offset);
dpixc_dip = noise(dpixc_nonoise,'poisson');
dpixc = double(dpixc_dip);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function sitem = cfg2struct(item)

% function sitem = cfg2struct(item)
% Return a struct containing all fields of item plus a field type. This is
% the method suitable for cfg_choice and repeat classes. It descends down
% the values field to convert the cfg_items in this field into structs.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg2struct.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

% Get parent struct, re-classify as field 'type'
sitem = cfg2struct(item.cfg_item);
sitem.type = class(item);

% Need to cycle through added fields
fn = mysubs_fields;
for k = 1:numel(fn)
    sitem.(fn{k}) = item.(fn{k});
end;

% Treat values{:} fields
if numel(item.values) > 0
    for k = 1:numel(item.values)
        sitem.values{k} = cfg2struct(item.values{k});
    end;
end;

author: spm


text:
% FT_PREAMBLE_PROVENANCE is a helper script that records the time and memory at the
% start of the function. At the end of the function FT_POSTAMBLE_PROVENANCE will
% record and store the time and memory in the output configuration together with
% information about the environment, such as the name of the user and computer, the
% MATLAB and FieldTrip version, etc.
%
% FieldTrip also attempts to uniquely identify the input and the output data. The
% code that deals with tracking the input data structures is found in
% FT_PREAMBLE_LOADVAR. The code that deals with tracking the information about the
% output data structures is found in FT_POSTAMBLE_HISTORY.
%
% Use as
%   ft_preamble provenance
%   .... regular code goes here ...
%   ft_postamble provenance
%
% See also FT_PREAMBLE, FT_POSTAMBLE, FT_POSTAMBLE_PROVENANCE

% Copyright (C) 2011-2016, Robert Oostenveld, DCCN
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% Record the start time and memory. These are used by ft_postamble_callinfo, which
% stores them in the output cfg.callinfo.  In the mean time, they are stored in the
% function workspace, which is why they should have cryptical names to prevent any
% variable name clashes.

if (isfield(cfg, 'trackcallinfo') && ~istrue(cfg.trackcallinfo))
  % do not track the call information
  return
end

% add the user-specified cfg (before any defaults handling etc.) to the callinfo
% some fields are for internal use only and should not be stored
cfg.callinfo.usercfg = removefields(cfg, ignorefields('provenance'));

if isfield(cfg, 'trackdatainfo') && istrue(cfg.trackdatainfo)
  % compute the MD5 hash of each of the input arguments
  % temporarily remove the cfg field for getting the hash (creating a duplicate of the data, but still has the same mem ref, so no extra mem needed)
  if isequal(iW1aenge_preamble, {'varargin'})
    tmpargin = varargin;
  else
    isvar = cellfun(@(x) exist(x, 'var')==1, iW1aenge_preamble);
    tmpargin = cellfun(@eval, iW1aenge_preamble(isvar), 'UniformOutput', false);
    tmpargin( isvar) = tmpargin;
    tmpargin(~isvar) = {[]};
    clear isvar
  end
  cfg.callinfo.inputhash = cell(1,numel(tmpargin));
  for iargin = 1:numel(tmpargin)
    tmparg = tmpargin{iargin}; % can't get number of bytes with whos unless taken out of it's cell
    if isfield(tmparg,'cfg')
      tmparg = rmfield(tmparg, 'cfg');
    end
    % only calculate md5 when below 2^31 bytes (CalcMD5 can't handle larger input)
    bytenum = whos('tmparg');
    bytenum = bytenum.bytes;
    if bytenum<2^31
      try
        cfg.callinfo.inputhash{iargin} = ft_hash(tmparg);
      catch
        % the mxSerialize mex file is not available on all platforms, do not compute a hash
        % http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2452
      end
    else
      % the data is too large, do not compute a hash
    end
  end
  clear tmpargin iargin tmparg bytenum; % remove the extra references
end

stack = dbstack('-completenames');
% stack(1) is this script
% stack(2) is the calling ft_postamble function
% stack(3) is the main FieldTrip function that we are interested in
stack = stack(3);

% add information about the FieldTrip and MATLAB version used to the configuration
try
  cfg.callinfo.fieldtrip = ft_version();
catch
  cfg.callinfo.fieldtrip = 'unknown';
end
cfg.callinfo.matlab    = version();
cfg.callinfo.computer  = lower(computer); % for example maci64, glnx86, ...

% add information about the execution environment to the configuration
cfg.callinfo.hostname = gethostname();
cfg.callinfo.user     = getusername();
cfg.callinfo.pwd      = pwd;
cfg.callinfo.calltime = clock();

% add information about the function filename and revision to the configuration
cfg.version.name = stack.file;
clear stack

% the revision number is maintained by SVN in the ft_revision variable in the calling function
if ~exist('ft_revision', 'var')
  cfg.version.id   = 'unknown';
else
  cfg.version.id   = ft_revision;
end

ftohDiW7th_FuncTimer = tic();
ftohDiW7th_FuncMem   = memtic();


text:
%% Copyright (C) 2006 Sylvain Pelissier <sylvain.pelissier@gmail.com>
%% Copyright (C) 2015-2016 Colin B. Macdonald <cbm@m.fsf.org>
%%
%% This program is free software; you can redistribute it and/or modify it under
%% the terms of the GNU General Public License as published by the Free Software
%% Foundation; either version 3 of the License, or (at your option) any later
%% version.
%%
%% This program is distributed in the hope that it will be useful, but WITHOUT
%% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
%% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
%% details.
%%
%% You should have received a copy of the GNU General Public License along with
%% this program; if not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defun dirac (@var{x})
%% Compute the Dirac delta (generalized) function.
%%
%% The Dirac delta ``function'' is a generalized function (or distribution)
%% which is zero almost everywhere, except at the origin where it is
%% infinite.
%%
%% Examples:
%% @example
%% @group
%% dirac (0)
%%   @result{} Inf
%% dirac (1)
%%   @result{} 0
%% dirac ([-10 -1 0 1 inf])
%%   @result{} 0     0   Inf     0     0
%% @end group
%% @end example
%% @seealso{heaviside, @@sym/dirac}
%% @end defun

function y = dirac(x)
  if (nargin ~= 1)
    print_usage ();
  end

  if (~isreal (x))
    error ('dirac: X must not contain complex values');
  end

  y = zeros (size (x), class (x));
  y(x == 0) = Inf;

  y(isnan (x)) = NaN;
end

%!assert (isinf (dirac (0)))
%!assert (dirac (1) == 0)
%!assert (isnan (dirac (nan)))
%!assert (isequaln (dirac ([-1 1 0 eps inf -inf nan]), [0 0 inf 0 0 0 nan]))
%!error <complex values> dirac (1i)
%!assert (isa (dirac (single (0)), 'single'))

repo: octsympy


text:
function pdf = arcsin_pdf ( x, a )

%*****************************************************************************80
%
%% ARCSIN_PDF evaluates the Arcsin PDF.
%
%  Discussion:
%
%    The LOGISTIC EQUATION has the form:
%
%      X(N+1) = 4.0D+00 * LAMBDA * ( 1.0D+00 - X(N) ).
%
%    where 0 < LAMBDA <= 1.  This nonlinear difference equation maps
%    the unit interval into itself, and is a simple example of a system
%    exhibiting chaotic behavior.  Ulam and von Neumann studied the
%    logistic equation with LAMBDA = 1, and showed that iterates of the
%    function generated a sequence of pseudorandom numbers with
%    the Arcsin probability density function.
%
%    The derived sequence
%
%      Y(N) = ( 2 / PI ) * Arcsin ( SQRT ( X(N) ) )
%
%    is a pseudorandom sequence with the uniform probability density
%    function on [0,1].  For certain starting values, such as X(0) = 0, 0.75,
%    or 1.0D+00, the sequence degenerates into a constant sequence, and for
%    values very near these, the sequence takes a while before becoming
%    chaotic.
%
%    PDF(X) = 1 / ( PI * Sqrt ( A**2 - X**2 ) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Daniel Zwillinger and Stephen Kokoska,
%    CRC Standard Probability and Statistics Tables and Formulae,
%    Chapman and Hall/CRC, 2000, pages 114-115.
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%    -A < X < A.
%
%    Input, real A, the parameter of the CDF.
%    A must be positive.
%
%    Output, real PDF, the value of the PDF.
%
  if ( a <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'ARCSIN_PDF - Fatal error!\n' );
    fprintf ( 1, '  Parameter A must be positive.\n' );
    error ( 'ARCSIN_PDF - Fatal error!' );
  end

  if ( x <= -a | a <= x )
    pdf = 0.0;
  else
    pdf = 1.0 / ( pi * sqrt ( a * a - x * x ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function cnoise_test02 ( m, n, var, r, alpha )

%*****************************************************************************80
%
%% CNOISE_TEST02 tests F_ALPHA_TGAUSSIAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 April 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of realizations to create.
%
%    Input, integer N, the size of the noise vector.
%
%    Input, real VAR, the variance of the Gaussian distribution.
%
%    Input, real R, defines the truncated range of the distribution.
%
%    Input, real ALPHA, the value of ALPHA.
%
  scale = exp ( 0.5 * ( 1.0 - alpha ) * log ( n ) );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CNOISE_TEST02:\n' );
  fprintf ( 1, '  Test F_ALPHA_TGAUSSIAN.\n' );
  fprintf ( 1, '  Noise is generated from a truncated Gaussian distribution\n' );
  fprintf ( 1, '  with a specified variance and range.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of realizations M = %d\n', m );
  fprintf ( 1, '  Dimension of noise vector N = %d\n', n );
  fprintf ( 1, '  Noise exponent ALPHA = %f\n', alpha );
  fprintf ( 1, '  Variance of Gaussian distribution VAR = %f\n', var );
  fprintf ( 1, '  Truncated range = [ %f, %f ]\n', -r, r );
  fprintf ( 1, '  Scale factor SCALE = %f\n', scale );

  res = zeros ( m, 1 );

  tic
  for i = 1 : m
    x = f_alpha_tgaussian ( n, var, r, alpha );
    x = scale * x;
    res(i) = mean ( x );
  end
  result_time = toc;

  result_mean = mean ( res );
  result_std = std ( res );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Mean value of noise = %f\n', result_mean );
  fprintf ( 1, '  STD of noise        = %f\n', result_std );
  fprintf ( 1, '  Computation time    = %f seconds\n', result_time );

  return
end


text:
function area = polygonArea(poly, varargin)
%POLYGONAREA Compute the signed area of a polygon.
%
%   A = polygonArea(POINTS);
%   Compute area of a polygon defined by POINTS. POINTS is a N-by-2 array
%   of double containing coordinates of vertices.
%
%   Vertices of the polygon are supposed to be oriented Counter-Clockwise
%   (CCW). In this case, the signed area is positive.
%   If vertices are oriented Clockwise (CW), the signed area is negative.
%
%   If polygon is self-crossing, the result is undefined.
%
%   Examples
%     % compute area of a simple shape
%     poly = [10 10;30 10;30 20;10 20];
%     area = polygonArea(poly)
%     area =
%         200
%
%     % compute area of CW polygon
%     area2 = polygonArea(poly(end:-1:1, :))
%     area2 =
%         -200
%
%     % Computes area of a paper hen
%     x = [0 10 20  0 -10 -20 -10 -10  0];
%     y = [0  0 10 10  20  10  10  0 -10];
%     poly = [x' y'];
%     area = polygonArea(poly)
%     area =
%        400
%
%     % Area of unit square with 25% hole
%     pccw = [0 0; 1 0; 1 1; 0 1];
%     pcw = pccw([1 4 3 2], :) * .5 + .25;
%     polygonArea ([pccw; nan(1,2); pcw])
%     ans =
%        0.75
%
%   References
%   algo adapted from P. Bourke web page
%   http://paulbourke.net/geometry/polygonmesh/
%
%   See also
%   polygons2d, polygonCentroid, polygonSecondAreaMoments, triangleArea
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2004-05-05
% Copyright 2004-2022 INRA - TPV URPOI - BIA IMASTE

%% Process special cases

% in case of polygon sets, computes the sum of polygon areas
if iscell(poly)
    area = 0;
    for i = 1:length(poly)
        area = area + polygonArea(poly{i});
    end
    return;
end

% check there are enough points
if size(poly, 1) < 2
    area = 0;
    return;
end

% case of polygons with holes -> computes the sum of areas
if any(isnan(poly))
    area = sum(polygonArea(splitPolygons(poly)));
    return;
end

%% Process single polygons or single rings

% extract coordinates
if nargin == 1
    % polygon given as N-by-2 array
    px = poly(:, 1);
    py = poly(:, 2);

elseif nargin == 2
    % poylgon given as two N-by-1 arrays
    px = poly;
    py = varargin{1};
end

% indices of next vertices
N = length(px);
iNext = [2:N 1];

% compute area (vectorized version)
area = sum(px .* py(iNext) - px(iNext) .* py) / 2;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function test_bug1916

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_scalingfactor

if isempty(which('ft_scalingfactor'))
  % it is in the fieldtrip/private directory
  [ftver, ftpath] = ft_version;
  cd(fullfile(ftpath, 'private'));
end

old = 'mT/mm';
new = 'fT/km';

% the first time for any conversion takes 150 ms
stopwatch = tic;
factor = ft_scalingfactor(old, new);
toc(stopwatch)

n = 300;
old = repmat({old}, 1, n);
new = repmat({new}, 1, n);

% this used to work fine, because all inputs are the same
stopwatch = tic;
factor = cellfun(@ft_scalingfactor, old, new);
toc(stopwatch)

n = 100;
old = repmat({'T/cm' 'T/cm' 'T'}, 1, n);
new = repmat({'fT/m' 'fT/m' 'fT'}, 1, n);

% this used to be very slow, because of the A-A-B pattern
stopwatch = tic;
factor = cellfun(@ft_scalingfactor, old, new);
toc(stopwatch)

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Calculate the MFCC from filter bank
% Author: Xiao Xiong
% Created: 4 Feb 2005
% Last modified: 4 Feb 2005

function [feature] = fbank2mfcc(fbank,logE,DO_BLIND_EQUALIZATION);

global bias;
bias = zeros(12,1);

[N_vector, N_melBank]= size(fbank);

% for i=1:N_vector
%     % calculate the DCT of the vector. The mfcc sequences is
%     % c1,c2,c3...c12,c0
%     feature(i,1:13) = mydct(fbank(i,:), 13, N_melBank);
%     % do blind equalization to reduce convolutional distortion
%     if DO_BLIND_EQUALIZATION == 1
%         feature(i,1:12) = blind_equal(feature(i,1:12), logE(i));
%     end
% end

% Faster implementation
feature(:,1:13) = mydct(fbank, 13, N_melBank);

author: singaxiong

repo: SignalGraph


text:
% Computes the gradient using "finite differences" and gives
% us a numerical estimate of the gradient.
function numgrad = debug_numerical_gradient(gradient_step, theta)
    % The following code implements numerical gradient checking, and
    % returns the numerical gradient. It sets numgrad(i) to (a numerical
    % approximation of) the partial derivative of J with respect to the
    % i-th input argument, evaluated at theta. (i.e., numgrad(i) should
    % be the (approximately) the partial derivative of J with respect
    % to theta(i).)
    numgrad = zeros(size(theta));
    perturb = zeros(size(theta));

    e = 1e-4;
    for p = 1:numel(theta)
        % Set perturbation vector
        perturb(p) = e;
        [loss1 gradients1] = gradient_step(theta - perturb);
        [loss2 gradients2] = gradient_step(theta + perturb);

        % Compute numerical gradient.
        numgrad(p) = (loss2 - loss1) / (2 * e);
        perturb(p) = 0;
    end
end

author: trekhleb

repo: machine-learning-octave


text:
function [HDR] = sclose(HDR)
% SCLOSE closes the file with the handle HDR
% [HDR] = sclose(HDR)
%    HDR.FILE.status = -1 if file could not be closed.
%    HDR.FILE.status = 0 indicates the file has been closed.
%
% see also: SOPEN, SREAD, SSEEK, STELL, SCLOSE, SWRITE

% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 3
% of the License, or (at your option) any later version.

%	$Id$
%	(C) 1997-2005,2006,2008 by Alois Schloegl <a.schloegl@ieee.org>
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/

if (HDR.FILE.FID<0) || ~HDR.FILE.OPEN,
        HDR.FILE.status = -1;
        %fprintf(HDR.FILE.stderr,'Warning SCLOSE (%s): invalid handle\n',HDR.FileName);
end;

if HDR.FILE.OPEN >= 2,          % write-open of files
	% check file length - simple test for file integrity
        EndPos = ftell(HDR.FILE.FID);          % get file length
        % force file pointer to the end, otherwise Matlab 6.5 R13 on PCWIN
        status = fseek(HDR.FILE.FID,0,'eof');

        if strcmp(HDR.TYPE,'BKR');
                if HDR.NS<1,
                        fprintf(HDR.FILE.stderr,'Error SCLOSE BKR: number of channels (HDR.NS) must be larger than zero.\n');
                        return;
                end;
                if HDR.NRec<1,
                        fprintf(HDR.FILE.stderr,'Error SCLOSE BKR: number of blocks (HDR.NRec) must be larger than zero.\n');
                        return;
                end;
                % check file length and write Headerinfo.
                HDR.SPR = (EndPos-HDR.HeadLen)/(HDR.NRec*HDR.NS*2);
                if isnan(HDR.SPR), HDR.SPR=0; end;
                if HDR.FILE.OPEN==3;
                        if any(isnan([HDR.NRec,HDR.NS,HDR.SPR,HDR.DigMax,HDR.PhysMax,HDR.SampleRate])), 	% if any unknown, ...
                                fprintf(HDR.FILE.stderr,'Error SCLOSE BKR: some important header information is still undefined (i.e. NaN).\n');
                                fprintf(HDR.FILE.stderr,'\t HDR.NRec,HDR.NS,HDR.SPR,HDR.DigMax,HDR.PhysMax,HDR.SampleRate must be defined.\n');
                                fprintf(HDR.FILE.stderr,'\t Try again.\n');
                        end;

                        fclose(HDR.FILE.FID);
                        HDR.FILE.FID = fopen(HDR.FileName,'r+');

                        count=fwrite(HDR.FILE.FID,HDR.VERSION,'int16');	        % version number of header
                        count=fwrite(HDR.FILE.FID,HDR.NS,'int16');	        % number of channels
                        count=fwrite(HDR.FILE.FID,HDR.SampleRate,'int16');      % sampling rate
                        count=fwrite(HDR.FILE.FID,HDR.NRec,'int32');            % number of trials: 1 for untriggered data
                        count=fwrite(HDR.FILE.FID,HDR.SPR,'uint32');            % samples/trial/channel
                        count=fwrite(HDR.FILE.FID,HDR.PhysMax,'int16');		% Kalibrierspannung
                        count=fwrite(HDR.FILE.FID,HDR.DigMax, 'int16');		% Kalibrierwert
                        count=fwrite(HDR.FILE.FID,zeros(4,1),'uint8');
                        count=fwrite(HDR.FILE.FID,[HDR.Filter.LowPass,HDR.Filter.HighPass],'float');

                        fseek(HDR.FILE.FID,32,'bof');
                        HDR.FLAG.TRIGGERED   = HDR.NRec>1;	% Trigger Flag
                        count = fwrite(HDR.FILE.FID,HDR.FLAG.TRIGGERED,'int16');           % FLAG TRIGGERED
                end;

	elseif strcmp(HDR.TYPE,'EDF') || strcmp(HDR.TYPE,'BDF') || strcmp(HDR.TYPE,'GDF'),
         	tmp = floor((EndPos - HDR.HeadLen) / HDR.AS.bpb);  % calculate number of records
                if isnan(tmp)
                	tmp = 0;
                end;

                if (HDR.NRec~=tmp)
                       if ~any(HDR.FILE.PERMISSION=='z')
                               HDR.NRec=tmp;
                               fseek(HDR.FILE.FID,236,'bof');
                               if strcmp(HDR.TYPE,'GDF')
                                        c=fwrite(HDR.FILE.FID,[HDR.NRec,0],'int32');
                               else
                                        fprintf(HDR.FILE.FID,'%-8i',HDR.NRec);
                                end;
                        else    %% due to a limitation in zlib
                                fprintf(HDR.FILE.stderr,'ERROR SCLOSE: number-of-records-field (HDR.NRec) could not be updated in file %s.\n',HDR.FileName);
                        end;
                end;

                if strcmp(HDR.TYPE,'GDF') && isfield(HDR,'EVENT'),
                	if ~all([HDR.NS, HDR.NRec, HDR.AS.bpb]>0)
                		HDR.AS.EVENTTABLEPOS = HDR.HeadLen;
                	else
	                        HDR.AS.EVENTTABLEPOS = HDR.HeadLen+HDR.AS.bpb*HDR.NRec;
	                end;
                        if isfield(HDR.EVENT,'VAL')
                                % handling of sparse (non-equidistant) sampling values
                                if isempty(HDR.EVENT.DUR)
                                        HDR.EVENT.DUR = zeros(size(HDR.EVENT.TYP));
                                end;
                                if (HDR.VERSION<1.90)
                                        fprintf(2,'Warning SCLOSE: GDF v2.0 or higher required for storing sparse samples but version is only %4.2f.\n',HDR.VERSION);
                                        %ix0 = ~isnan(HDR.EVENT.VAL);
                                        %HDR.EVENT.POS = HDR.EVENT.POS(~ix0);
                                        %HDR.EVENT.TYP = HDR.EVENT.TYP(~ix0);
                                        %HDR.EVENT.CHN = HDR.EVENT.CHN(~ix0);
                                        %HDR.EVENT.DUR = HDR.EVENT.DUR(~ix0);
                                        %HDR.EVENT = rmfield(HDR.EVENT,'VAL');
                                end;
	                end;
                        if isfield(HDR.EVENT,'VAL')
                                ix  = find(~isnan(HDR.EVENT.VAL));
                                tmp = unique(HDR.EVENT.CHN(ix));

                                flag.invalid = 0;
                                if (~all(HDR.EVENT.CHN(ix) > 0))
                                        flag.invalid = 1;
                                        fprintf(2,'Warning SCLOSE: Sparse sampling values without valid channel. Samples are not stored.\n');
                                        tmp = tmp(tmp > 0);
                                end;
                                if any(HDR.GDFTYP(tmp) > 6)
                                        fprintf(2,'Warning SCLOSE: Sparse sampling values must be of type integer using no more than 32 bits. Data is converted to uint32.\n');
                                end;

                                if any(HDR.AS.SPR(tmp)) || flag.invalid,
                                        fprintf(2,'Warning SCLOSE: Sparse sampling value for non-sparse channels not allowed. The following channel(s) is(are) affected: ');
                                        fprintf(2,'%i', tmp(HDR.AS.SPR(tmp)>0) );
                                        fprintf(2,'.  Samples are not stored.\n')
                                        ix0 = ~isnan(HDR.EVENT.VAL) & ~(HDR.EVENT.CHN>0); % initialize index with logical(0)
                                        for k = 1:length(tmp),
                                                if (HDR.AS.SPR(tmp(k))>0),
                                                        ix0 = ix0 | (HDR.EVENT.CHN==tmp(k));
                                                end;
                                        end;
                                        HDR.EVENT.POS = HDR.EVENT.POS(~ix0);
                                        HDR.EVENT.TYP = HDR.EVENT.TYP(~ix0);
                                        HDR.EVENT.CHN = HDR.EVENT.CHN(~ix0);
                                        HDR.EVENT.DUR = HDR.EVENT.DUR(~ix0);
                                        HDR.EVENT.VAL = HDR.EVENT.VAL(~ix0);
                                        ix = ~isnan(HDR.EVENT.VAL);
                                end;

                                % prepare for storing
                                HDR.EVENT.TYP(ix) = hex2dec('7fff');
                                HDR.EVENT.DUR(ix) = HDR.EVENT.VAL(ix);
                        end;

                        len = [length(HDR.EVENT.POS),length(HDR.EVENT.TYP)];
                        EVENT.Version = 1;
                        if isfield(HDR.EVENT,'CHN') && isfield(HDR.EVENT,'DUR'),
                                if any(HDR.EVENT.CHN) || any(HDR.EVENT.DUR),
                                        EVENT.Version = 3;
					len = [len,length(HDR.EVENT.CHN),length(HDR.EVENT.DUR)];
                                end;
                        end;

                        if any(len~=len(1))
                                fprintf(HDR.FILE.stderr,'Error SCLOSE-GDF: cannot write Event table, file %s not closed.\n',HDR.FileName);
				return;
                        else
                                status  = fseek(HDR.FILE.FID,HDR.HeadLen+HDR.AS.bpb*HDR.NRec,'bof');
                                %status = fseek(HDR.FILE.FID,0,'eof');
                                if ftell(HDR.FILE.FID)~=HDR.AS.EVENTTABLEPOS,
                                        fprintf(HDR.FILE.stderr,'Warning SCLOSE-GDF: inconsistent GDF file\n');
                                end
				if ~isfield(HDR.EVENT,'SampleRate'),
					HDR.EVENT.SampleRate = HDR.SampleRate;
				end;

				% write eventtable info: Version, SampleRate & length
				if (HDR.VERSION<1.94)
					tmp = HDR.EVENT.SampleRate;
					tmp = [EVENT.Version,mod(tmp,256),floor(mod(tmp,65536)/256),floor(tmp/65536)];
	                                fwrite(HDR.FILE.FID,tmp,'uint8');
	                                fwrite(HDR.FILE.FID,length(HDR.EVENT.POS),'uint32');
				else
					tmp = length(HDR.EVENT.POS);
					tmp = [EVENT.Version,mod(tmp,256),floor(mod(tmp,65536)/256),floor(tmp/65536)];
	                                fwrite(HDR.FILE.FID,tmp,'uint8');
	                                fwrite(HDR.FILE.FID,HDR.EVENT.SampleRate,'float32');
				end;

				% write table
                                c1 = fwrite(HDR.FILE.FID,HDR.EVENT.POS,'uint32');
                                c2 = fwrite(HDR.FILE.FID,HDR.EVENT.TYP,'uint16');
                                c3 = length(HDR.EVENT.POS); c4 = c3;
                                if EVENT.Version==3;
                                        c3 = fwrite(HDR.FILE.FID,HDR.EVENT.CHN,'uint16');
                                        c4 = fwrite(HDR.FILE.FID,HDR.EVENT.DUR,'uint32');
                                end;
                                if any([c1,c2,c3,c4]~=length(HDR.EVENT.POS))
                                        fprintf(2,'\nError SCLOSE: writing of EVENTTABLE failed. File %s not closed.\n', HDR.FileName);
                                        return;
                                end
                        end;
                end;

        elseif strcmp(HDR.TYPE,'CFWB');
                tmp = (EndPos-HDR.HeadLen)/HDR.AS.bpb;
                if isnan(tmp), tmp=0; end;
                if (tmp~=HDR.NRec);
                        if ~any(HDR.FILE.PERMISSION=='z')
                                HDR.NRec = tmp;
                                fseek(HDR.FILE.FID,14*4,-1);
                                count = fwrite(HDR.FILE.FID,HDR.NRec,'int32');           % channels
                        else
                                fprintf(HDR.FILE.stderr,'ERROR SCLOSE (CFWB): number-of-samples-field (HDR.SPR) could not be updated in file %s.\n',HDR.FileName);
                        end;
                end;

        elseif strcmp(HDR.TYPE,'SND');
                tmp = (EndPos-HDR.HeadLen)/HDR.AS.bpb;
                if isnan(tmp), tmp=0; end;
                if (HDR.FILE.OPEN==3) && (tmp~=HDR.SPR);
                        if ~any(HDR.FILE.PERMISSION=='z')
                                HDR.SPR = tmp;
                                fseek(HDR.FILE.FID,8,'bof');
                                count = fwrite(HDR.FILE.FID,HDR.SPR*HDR.AS.bpb,'uint32');           % bytes
                                fseek(HDR.FILE.FID,20,-1);
                                count = fwrite(HDR.FILE.FID,HDR.NS,'uint32');           % channels
                        else
                                fprintf(HDR.FILE.stderr,'ERROR SCLOSE (SND): fields HDR.SPR and HDR.NS could not be updated in file %s.\n',HDR.FileName);
                        end;
		end;

        elseif strcmp(HDR.TYPE,'AIF');
                if HDR.FILE.OPEN==3;
                        fseek(HDR.FILE.FID,4,-1);
                        count = fwrite(HDR.FILE.FID,EndPos-8,'uint32');           % bytes
                        fseek(HDR.FILE.FID,HDR.WAV.posis(2),-1);
                        count = fwrite(HDR.FILE.FID,EndPos-4-HDR.WAV.posis(2),'uint32');           % channels
                end;

        elseif strcmp(HDR.TYPE,'FLT') ;
		% warning FLT/header might have changed in needs to be rewritten
        elseif strcmp(HDR.TYPE,'WAV') ;
                if HDR.FILE.OPEN==3;
                        fseek(HDR.FILE.FID,4,-1);
                        count = fwrite(HDR.FILE.FID,EndPos-16,'uint32');           % bytes
                        fseek(HDR.FILE.FID,HDR.WAV.posis(2),-1);
                        count = fwrite(HDR.FILE.FID,EndPos-4-HDR.WAV.posis(2),'uint32');           % channels
                end;
        end;
end;

if 0,
elseif strcmp(HDR.TYPE,'ZIP')
        [SUCCESS,MESSAGE,MESSAGEID] = rmdir(HDR.ZIP.TEMPDIR,'s');

elseif strcmp(HDR.TYPE,'FIF') && HDR.FILE.OPEN;
        global FLAG_NUMBER_OF_OPEN_FIF_FILES
        rawdata('close');
        HDR.FILE.OPEN = 0;
        HDR.FILE.status = 0;
        FLAG_NUMBER_OF_OPEN_FIF_FILES = FLAG_NUMBER_OF_OPEN_FIF_FILES-1;

elseif strmatch(HDR.TYPE,{'ADI','GTEC','LABVIEW','MAT','MAT4','MAT5','native','XML','XML-FDA','SierraECG'}),
        HDR.FILE.OPEN = 0;

elseif HDR.FILE.OPEN,
        HDR.FILE.OPEN = 0;
        HDR.FILE.status = fclose(HDR.FILE.FID);
end;


text:
function jed = transition_to_jed_jed ( )

%*****************************************************************************80
%
%% TRANSITION_TO_JED_JED returns the JED calendar transition as a JED.
%
%  Discussion:
%
%    In Scaliger's design of the JED, three cycles with different periods
%    began on JED = 0.  These three cycles coincide once more on the
%    transition day.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 December 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real JED, the Julian Ephemeris Date of the date.
%
  jed = 2913943.0;

  return
end


text:
function [] = getStateFiles(model_3d,statesDir,goodInds)

globals;
PASCAL_DIR = PASCAL_DIR;
if(nargin<3)
    goodInds = true(length(model_3d.c),1);
end
params = get_params();
viewClusters = cluster_views(model_3d,find(goodInds(:)),params.vis.numRotClusters,0);
%pBar =  TimedProgressBar( length(model_3d.c), 30, ...
%    'Generating state files ', ', completed ', 'Generated state files in ' );

parfor instance = 1:length(model_3d.c)

    if(~goodInds(instance))
        continue;
    end
    %% Read image
    if(strcmp(model_3d.voc_image_id{instance}(end-2:end),'jpg')~=0)
        fname = [PASCAL_DIR model_3d.voc_image_id{instance}];
    else
        fname = [PASCAL_DIR model_3d.voc_image_id{instance} '.jpg'];
    end
    im = imread(fname);

    %% Get Segmentation from pascal annotations
     %disp(['Instance ' num2str(instance)]);
     if(isfield(model_3d,'voc_rec_id'))
         voc_id = [model_3d.voc_image_id{instance} '_' num2str(model_3d.voc_rec_id(instance))];
     else
         voc_id = [model_3d.voc_image_id{instance} '_1'];
     end
    if(isfield(model_3d,'flip'))
        flip = model_3d.flip(instance);
    else
        flip = 0;
    end

    if(flip)
        im(:,:,1) = fliplr(im(:,:,1));im(:,:,2) = fliplr(im(:,:,2));im(:,:,3) = fliplr(im(:,:,3));
        voc_id = ['flip_' voc_id];
    end

    segs = model_3d.seg;
    if(iscell(segs))
        mask = segs{instance};
    else
        mask = roipoly(im,segs.poly_x{instance},segs.poly_y{instance});
    end

    [~,occlusionIDX] = bwdist(mask);
    %mask = roipoly(im,segs.poly_x{instance}{1},segs.poly_y{instance}{1});

    %% Prepare data for instance from model structure
    bbox = model_3d.bbox(instance,:); % Used this for bbox scaling!
    % Keep all the information in bbox while scaling
    bbox = [floor(bbox(1:2)) ceil(bbox(3:4))];
    sc = max(bbox(3),bbox(4))/model_3d.params.nrsfm.norm_dim;
    nImages = size(model_3d.points3,1)/3;
    tr = bbox(1,1:2)';
    vis_points = [model_3d.points3(instance,:);model_3d.points3(nImages+instance,:);...
        model_3d.points3(2*nImages+instance,:)];
    vis_points = sc*vis_points;
    vis_points(1:2,:) = vis_points(1:2,:)+repmat(tr,[1,size(model_3d.points3,2)]);
    %vis_points(3,:) = -vis_points(3,:);

    %% visualization

    %subplot(2,2,1), imagesc(im); rectangle('Position',bbox,'EdgeColor','r','LineWidth',1);axis equal;
    %hold on;plot(vis_points(1,:),vis_points(2,:),'r.','MarkerSize',20);
    %hold on; trimesh(f',vis_points(1,:),vis_points(2,:),'LineWidth',3);axis equal;
    %subplot(2,2,2), imagesc(dmap);colormap('jet');axis equal;
    %hold on; trimesh(f',vis_points(1,:),vis_points(2,:),'LineWidth',3);axis equal;
    %subplot(2,2,3), imagesc(mask);colormap('gray');axis equal;
    %subplot(2,2,4),trimesh(f,vis_points(1,:),vis_points(2,:));axis equal;
    %pause();close all;

    %% Setup state file output
    state=struct;
    %state.dmapCoarse = dmap;
    state.mask = mask;
    state.gtmask= mask;
    state.im=im;
    state.voc_id = voc_id;
    state.kps = vis_points;
    state.bbox = bbox;
    %state.cameraPos = cameraPos;
    %state.cameraF = cameraF;
    state.cameraRot = model_3d.rots{instance};
    state.cameraScale = model_3d.c(instance)*sc;
    state.def_weights = model_3d.def_weights(instance,:);
    state.translation = mean(state.kps(1:2,:),2);

    state.transform3d = giveTransformFunc3d(state.cameraRot,state.cameraScale,state.translation);
    state.invtransform3d = giveInvTransformFunc3d(state.cameraRot,state.cameraScale,state.translation);
    state.bdry = getBoundary(state.mask);
    state.occlusionIDX = occlusionIDX;
    state.viewClusters = viewClusters;
    state.globalID = instance;
    savefunc(fullfile(statesDir, [voc_id '.mat']), state);
    %pBar.progress();
end
%pBar.stop();
end

function savefunc(fname,state)
save(fname,'state');
end

function tfun3d = giveTransformFunc3d(R,c,tr)

function S3d = tf(S)
S3d = c*R*S';
tr3 = [tr;zeros(1,size(tr,2))];
S3d = bsxfun(@plus,S3d,tr3);
S3d = S3d';
end

tfun3d = @tf;
end

function tfun3d = giveInvTransformFunc3d(R,c,tr)

function S3d = invtf(S)
tr3 = [tr' zeros(size(tr,2),1)];
S3d = bsxfun(@minus,S,tr3);
S3d = (S3d*R)/c;
end

tfun3d = @invtf;
end

function bdry = getBoundary(mask)
f = [0 1 0;1 0 1; 0 1 0];
[I,J] = find(mask.*(conv2(double(mask),f,'same')<4));
bdry = [J,I];
end

repo: CategoryShapes


text:
function TessMat = in_tess_vtk(TessFile)
% IN_TESS_VTK: Read Visualization Toolkit .vtk mesh files
%
% USAGE:  TessMat = in_tess_vtk(TessFile);
%
% INPUT:
%     - TessFile : full path to a tesselation file
% OUTPUT:
%     - TessMat:  Brainstorm tesselation structure
% FORMAT:
%     ASCII file with the following structure:
%     - Header line 1: "# vtk DataFile Version ???"
%     - Header line 2: "surface file"
%     - Header line 3: "ASCII"
%     - Header line 4: "DATASET POLYDATA"
%     - Header line 5: "POINTS <nVertices>  float"
%     - Vertices     : nVertices lines ("x y z"), values in miliimeters
%     - Header line 6: "POLYGONS <nFaces> <nValues>"
%     - Faces        : nFaces lines ("nvert vertex1 vertex2 vertex3")
%
% SEE ALSO: in_tess

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2012

% Open tesselation file
fid = fopen(TessFile, 'r');
if fid < 0
    error(['Cannot open file ', TessFile])
end

% ===== READ HEADER =====
% Read 5 lines of header
h1 = fgetl(fid);
h2 = fgetl(fid);
h3 = fgetl(fid);
h4 = fgetl(fid);
h5 = fgetl(fid);
% Check file format
if isempty(strfind(h1, 'vtk'))
    error('Not a valid VTF ASCII mesh file.');
end
% Get number of vertices
nVertices = sscanf(h5, '%*s %d %*s', 1);

% ===== READ MESH =====
% Read vertices
Vertices = double(fscanf(fid, '%f', [3 nVertices])) / 1000;
% Go to next line
fgetl(fid);
% Read number of faces
h6 = fgetl(fid);
res = textscan(h6, '%s %d', 1);
if ~iscell(res{1}) || isempty(res{1}) || ~ischar(res{1}{1}) || ~isequal(res{1}{1}, 'POLYGONS')
    error(['Unsupported faces format. Only "POLYGONS" is supported.']);
end
nFaces = res{2};
% Read faces
Faces = double(fscanf(fid, '%f',[4 nFaces]) + 1);
Faces = Faces(2:4,:);
% Close file
fclose(fid);

%% ===== CONVERT IN BRAINSTORM FORMAT =====
TessMat.Vertices = Vertices';
TessMat.Faces    = Faces';

author: brainstorm-tools


text:
%PRTESTC Test routine for the PRTOOLS classifier
%
% This script tests a given, untrained classifier w, defined in the
% workspace, e.g. w = my_classifier. The goal is to find out whether
% w fulfills all the requirements of a PRTools classifier.
%

if exist('w') ~= 1 | ~ismapping(w) | ~isuntrained(w)
	error('No untrained classifier w found')
end

m = 50;
a = gendath([m,m]);
a = setprior(a,0);
[b,c] = gendat(a,0.5);
name = getname(w);
if (isempty(name))
	error('No name found for the untrained classifier.')
end

disp('')
disp(['     Testing Classifier ' name])
disp(['     ------------------ '])
v = b*w;
newfig(1,3);
scatterd(b);
plotc(v);
if (isempty(getname(v)))
	disp('No name found for the classifier.')
end
disp('Classification error for the Higleyman data: ')
disp(c*v*testc)

disp('Direct output: ')
disp(c(1:5,:)*v)

disp('Classifier output: ')
disp(c(1:5,:)*v*classc)

disp('Soft labels, classification error: ')
a = gendath([m,m],'soft');
a = setprior(a,0);
[b,c] = gendat(a,0.5);
disp(c*(b*w)*testc);

newfig(2,3);
learnsizes = [3,5,7,10,15];
e = cleval(w,b,learnsizes,2,c);
plote(e)
a = gendatm(repmat(20,1,8));
a = setprior(a,0);
[b,c] = gendat(a,0.5);
newfig(3,3)
scatterd(b)
u = b*w;
plotc(u)
disp('Classification error for a multi-class problem: ')
disp(c*u*testc)

if (isaffine(v))
	load nist16_38
	a = prdataset(a);
	v = a*w;
	newfig(6,3);
	show(v);
end

author: marianux

repo: ecg-kit

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


text:
%% createGoToWBC.m
%
%  Run this script once create the goToWholeBodyControllers.m script in your userpath.

clc
fprintf('\n## whole-body-controllers ##\n');
fprintf('\nCreating "goToWholeBodyControllers.m" in your userpath...\n\n');

% Path to the Matlab userpath
pathToUserpath        = userpath;
pathSeparatorLocation = strfind(pathToUserpath, pathsep);

if isempty(pathToUserpath)

    answer = input('Empty userpath. Do you want to reset the userpath? Y/N ','s');

    if strcmpi(answer,'Y')

       userpath('reset');
       disp('Resetting userpath..');
       pathToUserpath        = userpath;
       pathSeparatorLocation = strfind(pathToUserpath, pathsep);
    else
       error('Please set the userpath before running this script');
    end

elseif size(pathSeparatorLocation, 2) > 1

    answer = input('Multiple userpath. Do you want to reset the userpath? Y/N ','s');

    if strcmpi(answer,'Y')

       userpath('reset');
       disp('Resetting userpath..');
       pathToUserpath        = userpath;
       pathSeparatorLocation = strfind(pathToUserpath, pathsep);
    else
       error('Please set a single userpath before running this script');
    end
end

% check again the userpath
if isempty(pathToUserpath)

    error('userpath is still empty. Please set the userpath before running this script');

elseif size(pathSeparatorLocation, 2) > 1

    error('There are still multiple userpath. Please set a single userpath before running this script');
end

% save a script named "goToWholeBodyControllers" inside the pathdef folder,
% to facilitate the user to reach the WBC working folder

% get the WBC path
[~, currentFolderName, ~] = fileparts(pwd);
currentPath               = pwd;
wbcPath                   = currentPath(1:end-length(currentFolderName));

% note: exist returns 2 if a .m file with the specified name is found
if exist([pathToUserpath, filesep, 'goToWholeBodyControllers.m'],'file') == 2

    fprintf('\n')
    warning(['A goToWholeBodyControllers.m file exists already in your ',pathToUserpath, ' folder,' ...
             ' and therefore it won''t be created.'])
else
    % create the script and write inside the cd command
    fid = fopen([pathToUserpath, filesep, 'goToWholeBodyControllers.m'],'w');
    fprintf(fid, ['cd ', wbcPath, ';']);
    fclose(fid);

    fprintf('\n')
    fprintf(['A file called goToWholeBodyControllers.m has been created in your %s folder.\n', ...
             'This will help to quickly reach the WBC-project folder after ', ...
             'Matlab is launched.\n'], pathToUserpath);
end

author: robotology

repo: whole-body-controllers


text:
% Reflective multivariate slice sampling

function [func, fx] = mcmc_init_reflective(x_init, get_logpdf, get_dlogpdf, ...
                                           epsilon, L, varargin)

options = struct( ...
    'fx',    @func_x_default, ...
    'type',  'inside');
% $$$     'y',     @(x) x,  ...        % transform the variable
% $$$     'logdy', @(x) 0,  ...        % log-jacobi of the transformation
% $$$     'dy',    @(x) 1,  ...
% $$$     'ddy',   @(x) 0,  ...
% $$$     'f',     @(y, varargin) []); % function for making pre-evaluations for logpdf

% Parse arguments
[options, errmsg] = argparse( options, varargin{:} );
error(errmsg);

% Initialization
x_current = x_init;
[fx_current, dfx_current] = options.fx(x_current);
logpdf_current = get_logpdf(fx_current);

fx = fx_current;

x_proposal = x_current;

switch options.type

 case 'outside'
  func = @reflective_outside;

 case 'inside'
  func = @reflective_inside;

 case default
  error('Unknown type for slice sampling');

end

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

  function [x, fx] = reflective_inside(varargin)
  x_trajectory = x_current;
  fx_trajectory = fx_current;
  dfx_trajectory = dfx_current;
  logpdf_trajectory = logpdf_current;
  lp_trajectory = logpdf_current(varargin{:});

  % Momentum stuff
  step = exprnd(epsilon); % step length
  p = randn(size(x_current));

  %traj = [];

  % Simulate N steps
  for n=1:L

    if mod(n-1,ceil(L/10)) == 0
      % Update the slice
      lp_slice = lp_trajectory + log(rand());
    end

    % Move a step forward
    x_proposed = x_trajectory + step*p;
    % Check whether the new point is inside the slice
    [fx_proposed, dfx_proposed] = options.fx(x_proposed);
    logpdf_proposed = get_logpdf(fx_proposed);
    lp_proposed = logpdf_proposed(varargin{:});
    if lp_proposed <= lp_slice
      % Reflect based on the gradient at the inside point

      % Get the gradient
      dlogpdf = get_dlogpdf(dfx_trajectory);
      g = dlogpdf(varargin{:});

      % New momentum
      g2 = g'*g;
      if g2 > 0
        p = p - (2*(p'*g)/g2)*g;
      else
        % Random or opposite direction?
        p = -p;%orth(randn(numel(p))) * p;
      end

      % Check the reflected point is inside the slice
      x_reflect = x_trajectory + step*p;
      [fx_reflect, dfx_reflect] = options.fx(x_reflect);
      logpdf_reflect = get_logpdf(fx_reflect);
      lp_reflect = logpdf_reflect(varargin{:});
      if lp_reflect <= lp_slice
        disp('Reject trajectory')
        % TODO: Should I keep the state or ignore
        x = x_current;
        fx = fx_current;
        %x = nan;
        %fx = nan;
        return
      end

      %
      % TODO: HOW TO REQUEST THE GRADIENT?!?!?!!!!
      %

      % Verify opposite direction
      x_verify = x_trajectory - step*p;
      fx_verify = options.fx(x_verify);
      logpdf_verify = get_logpdf(fx_verify);
      lp_verify = logpdf_verify(varargin{:});
      if lp_verify > lp_slice
        disp('Reject trajectory')
        % TODO: Should I keep the state or ignore
        x = x_current;
        fx = fx_current;
        %x = nan;
        %fx = nan;
        return
      end

      % Accept the reflected step
      %disp('Reflect')
      x_trajectory = x_reflect;
      fx_trajectory = fx_reflect;
      dfx_trajectory = dfx_reflect;
      lp_trajectory = lp_reflect;
      logpdf_trajectory = logpdf_reflect;

    else

      % Accept the forward step
      x_trajectory = x_proposed;
      fx_trajectory = fx_proposed;
      dfx_trajectory = dfx_proposed;
      lp_trajectory = lp_proposed;
      logpdf_trajectory = logpdf_proposed;

    end
%    traj = [traj, x_trajectory];
  end

% $$$   figure
% $$$   plot(traj(1,:),traj(2,:), 'kx-');
% $$$   % Plot contours
% $$$   v = axis();
% $$$   vx = linspace(v(1), v(2), 20);
% $$$   vy = linspace(v(3), v(4), 10);
% $$$   [VY,VX] = meshgrid(vy,vx);
% $$$   V = [VX(:), VY(:)]';
% $$$   F = zeros(size(VX));
% $$$   for n=1:numel(F)
% $$$     logpdf = get_logpdf(V(:,n));
% $$$     F(n) = logpdf();
% $$$   end
% $$$   hold on
% $$$   contour(VX,VY,reshape(F,size(VX)));
% $$$   error('jou')

  disp('Accept trajectory')
  x_current = x_trajectory;
  fx_current = fx_trajectory;
  dfx_current = dfx_trajectory;
  logpdf_current = logpdf_trajectory;

  x = x_current;
  fx = fx_current;
  end

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

  function [y, fy] = reflective_outside(varargin)
  disp('Check gradient')
  mycheckgrad(@helper, x_current, 1e-6)
  error('Stopped because wanted to check gradients')
    function [f, df] = helper(x)
    [fx, dfx] = options.fx(x);
    logpdf = get_logpdf(fx);
    dlogpdf = get_dlogpdf(dfx);
    f = logpdf(varargin{:});
    df = dlogpdf(varargin{:});
    end
  end

end

function [fx, dfx] = func_x_default(x)
fx = x;
dfx = x;
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [deltaPoseNoisy] = poseSubNoisy(p_j , p_i, sigmaT, sigmaR, isUniform)
% pose_sub = (p_j , p_i) = inv(p_i) * p_j + noise
% returns the planar roto-translation that tranform pj in pi
% plus noise (sigmaT is the std of the noise on the cartesian components,
% sigmaR is the std of the noise added to the angle before wrapping)

if nargin < 5
   isUniform = 0; % Gaussian noise by default
end

if nargin == 2
  sigmaT = 0;
  sigmaR = 0;
end

p_j = p_j(:);
p_i = p_i(:);

rho_i = p_i(1:2);
rho_j = p_j(1:2);
th_i = p_i(3);
th_j = p_j(3);

R = rot2D(th_i);

delta_th = wrapToPi(th_j - th_i);

if(abs(delta_th) > pi)
    disp('Error in pose_sub')
    disp(delta_th)
end

deltaPose = [ R' * (rho_j   -  rho_i)
    delta_th  ];

if isUniform == 0 % Gaussian noise
    noise = [sigmaT*randn(2,1) ;  sigmaR*randn];
else % uniform noise in the interval [-sigma, + sigma]
    noise = [sigmaT*2*(rand(2,1)-0.5*ones(2,1)) ; sigmaR*2*(rand-0.5)];
end

deltaPoseNoisy = deltaPose + noise;
deltaPoseNoisy(3) = wrapToPi(deltaPoseNoisy(3)); % wrap angle

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function [energy3M,entropy3M,sumAvg3M,corr3M,invDiffMom3M,contrast3M,...
    clustShade3M,clustPromin3M,haralCorr3M] = textureByPatchAllOffsFromCtr(scanArray3M, nL, ...
    patchSizeV, flagv, hWait, minIntensity, maxIntensity)
% function [energy3M,entropy3M,sumAvg3M,corr3M,invDiffMom3M,contrast3M,...
%     clustShade3M,clustPromin3M,haralCorr3M] = textureByPatchAllOffsFromCtr(scanArray3M, nL, ...
%     patchSizeV, flagv, hWait, minIntensity, maxIntensity)
%
% This function calculates GLCM texture from patches w.r.t. their centers.
%
% INPUTS:
% scanArray3M - Scan matrix with NaN values for voxels that are not to be
% included in texture calculation.
% nL - Number of gray levels
% patchSizeV - 3-element vector for the patch radius. For example, [3 3 0]
% flagv - 9-element vector of booleans to toggle calculation of texture features.
% hWait - (optional) waitbar handle. Use NaN to turn off waitbar.
% minIntensity - (optional) minimum intensity for discretization.
% maxIntensity - (optional) maximum intensity for discretization.
%
% Example:
% scan3M(~mask3M) = NaN;
% flagV = zeros(1,9);
% flagV(2) = 1; % entropy only
% nL = 64;
% patchSizeV = [3 3 0];
% [~,entropy3M] = textureByPatchAllOffsFromCtr(scan3M, nL, ...
%     patchSizeV, offsetsM, flagV);
%
% APA, 11/07/2018

% Generate flags
if ~exist('flagv','var')
    flagv = ones(1,9);
elseif exist('flagv','var') && isempty(flagv)
    flagv = ones(1,9);
end

% Flag to draw waitbar
waitbarFlag = 0;
if exist('hWait','var') && ishandle(hWait)
    waitbarFlag = 1;
end

% Get indices of non-NaN voxels
calcIndM = ~isnan(scanArray3M);

slcWindow = 2 * patchSizeV(3) + 1;
rowWindow = 2 * patchSizeV(1) + 1;
colWindow = 2 * patchSizeV(2) + 1;

% Build distance matrices
numColsPad = floor(colWindow/2);
numRowsPad = floor(rowWindow/2);
numSlcsPad = floor(slcWindow/2);

% Get number of voxels per slice
[numRows, numCols, numSlices] = size(scanArray3M);
numVoxels = numRows*numCols;

% Quantize the image
if exist('minIntensity','var') && exist('maxIntensity','var')
    q = imquantize_cerr(scanArray3M,nL,minIntensity,maxIntensity);
else
    q = imquantize_cerr(scanArray3M,nL);
end

% Pad q, so that sliding window works also for the edge voxels
%scanArrayTmp3M = padarray(scanArray3M,[numRowsPad numColsPad
%numSlcsPad],NaN,'both'); % aa commented
if exist('padarray.m','file')
    q = padarray(q,[numRowsPad numColsPad numSlcsPad],NaN,'both');
else
    q = padarray_oct(q,[numRowsPad numColsPad numSlcsPad],NaN,'both');
end

nanFlag = 0; % the quantized image is always padded with NaN. Hence, always,
% nanFlag = 1.
if any(isnan(q(:)))
    nanFlag = 1;
    q(isnan(q)) = nL+1;
end
lq = nL + 1;

q = uint16(q); % q is the quantized image

% Create indices for 2D blocks
[m,n,~] = size(q);
m = uint32(m);
n = uint32(n);
colWindow = uint32(colWindow);
rowWindow = uint32(rowWindow);
slcWindow = uint32(slcWindow);

% Index calculation adapted from
% http://stackoverflow.com/questions/25449279/efficient-implementation-of-im2col-and-col2im

% Start indices for each block
start_ind = reshape(bsxfun(@plus,[1:m-rowWindow+1]',[0:n-colWindow]*m),[],1);

% Row indices
lin_row = permute(bsxfun(@plus,start_ind,[0:rowWindow-1])',[1 3 2]);

% Get linear indices based on row and col indices and get desired output
indM = reshape(bsxfun(@plus,lin_row,(0:colWindow-1)*m),rowWindow*colWindow,[]);

% Indices of last level to filter out
nanIndV = false([lq*lq,1]);
if nanFlag
    nanIndV([lq:lq:lq*lq-lq, lq*lq-lq:lq*lq]) = true;
end

% Build levels vector for mu, sig
levRowV = repmat(1:lq,[1 lq]);
levColV = repmat(1:lq,[lq 1]);
levColV = levColV(:)';

% Build list of indices for px and contrast calculation
numElems = nL*nL; % APA 7/13
indCtrstM = false(nL,numElems); % APA 7/13
indPxM = false(nL,numElems); % APA 7/13
indPxPlusYm = false(nL,numElems); % APA 7/13
for n=0:nL-1
    % indices for p(x-y), contrast
    indCtrstV = false(lq*lq,1);
    indCtrst1V = 1:lq-n;
    indCtrst2V = 1+n:lq;
    indCtrstTmpV = indCtrst1V + (indCtrst2V-1)*lq;
    indCtrstTmpV = [indCtrstTmpV indCtrst2V + (indCtrst1V-1)*lq];
    indCtrstV(indCtrstTmpV) = 1;
    indCtrstV(nanIndV) = [];
    % indCtrstC{n+1} = indCtrstV;
    indCtrstM(n+1,:) = indCtrstV;

    % indices for px
    indPxV = false(lq*lq,1);
    indPxV(lq*n+1:lq*(n+1)) = true;
    indPxV(nanIndV) = [];
    % indPxC{n+1} = indPxV;
    indPxM(n+1,:) = indPxV;

end
for n=1:2*nL
    % indices for p(x+y)
    indPxPlusYv = false(lq*lq,1);
    indPxPlusYv(levRowV + levColV == n) = 1;
    indPxPlusYv(nanIndV) = [];
    % indPxPlusYc{n} = indPxPlusYv;
    indPxPlusYm(n,:) = indPxPlusYv;
end

% Build linear indices column/row-wise for Symmetry
indRowV = zeros(1,lq*lq);
for i=1:lq
    indRowV((i-1)*lq+1:(i-1)*lq+lq) = i:lq:lq*lq;
end

% Filter out NaN levels
levRowV(nanIndV) = [];
levColV(nanIndV) = [];

% Initialize
energy3M = [];
entropy3M = [];
sumAvg3M = [];
corr3M = [];
invDiffMom3M = [];
contrast3M = [];
clustShade3M = [];
clustPromin3M = [];
haralCorr3M = [];
if flagv(1)
    energy3M = zeros([numRows, numCols, numSlices],'single');
end
if flagv(2)
    entropy3M = zeros([numRows, numCols, numSlices],'single');
end
if flagv(3)
    sumAvg3M = zeros([numRows, numCols, numSlices],'single');
end
if flagv(4)
    corr3M = zeros([numRows, numCols, numSlices],'single');
end
if flagv(5)
    invDiffMom3M = zeros([numRows, numCols, numSlices],'single');
end
if flagv(6)
    contrast3M = zeros([numRows, numCols, numSlices],'single');
end
if flagv(7)
    clustShade3M = zeros([numRows, numCols, numSlices],'single');
end
if flagv(8)
    clustPromin3M = zeros([numRows, numCols, numSlices],'single');
end
if flagv(9)
    haralCorr3M = zeros([numRows, numCols, numSlices],'single');
end

tic
% Iterate over slices. compute cooccurance for all patches per slice
for slcNum = 1:numSlices
    disp(['--- Texture Calculation for Slice # ', num2str(slcNum), ' ----'])
    if flagv(1), energyV = zeros(1,numVoxels,'single'); end
    if flagv(2), entropyV = zeros(1,numVoxels,'single'); end
    if flagv(3), sumAvgV = zeros(1,numVoxels,'single'); end
    if flagv(4), corrV = zeros(1,numVoxels,'single'); end
    if flagv(5), invDiffMomV = zeros(1,numVoxels,'single'); end
    if flagv(6), contrastV = zeros(1,numVoxels,'single'); end
    if flagv(7), clustShadeV = zeros(1,numVoxels,'single'); end
    if flagv(8), clustProminV = zeros(1,numVoxels,'single'); end
    if flagv(9), haralCorrV = zeros(1,numVoxels,'single'); end

    calcSlcIndV = calcIndM(:,:,slcNum);
    calcSlcIndV = calcSlcIndV(:);
    numCalcVoxs = sum(calcSlcIndV);
    %indSlcM = indM(:,calcSlcIndV);
    indSlcM = indM(:,calcSlcIndV);
    indCurrVox = ceil(size(indM,1)/2);
    indCurrVoxV = indSlcM(indCurrVox,:);
    %slcM = uint32(q(numRowsPad+(1:numRows),numColsPad+(1:numCols),slcNum));
    % List of Voxel numbers
    voxelNumsV = uint32(0:lq*lq:lq*lq*(numCalcVoxs-1));
    slcM = uint32(q(:,:,slcNum));
    voxelValsV = slcM(indCurrVoxV);
    slc2M = slcM(indSlcM);
    slc2M = bsxfun(@plus, slc2M, (voxelValsV-1)*lq + voxelNumsV);
    %slc2M = bsxfun(@plus,slc2M,voxelNumsV);
    cooccurPatchM = accumarray(slc2M(:) ...
        ,1, [lq*lq*numCalcVoxs,1],[],...
        [],true); % patch-wise cooccurance
    cooccurPatchM = reshape(cooccurPatchM,lq*lq,numCalcVoxs);
    cooccurPatchM = cooccurPatchM + cooccurPatchM(indRowV,:); % for symmetry
    cooccurPatchM(nanIndV,:) = [];
    cooccurPatchM = bsxfun(@rdivide,cooccurPatchM,sum(cooccurPatchM)+1e-5);
    cooccurPatchM = full(cooccurPatchM);

    % Angular Second Moment (Energy)
    if flagv(1)
        energyV(calcSlcIndV) = energyV(calcSlcIndV) + sum(cooccurPatchM.^2);
    end
    % Entropy
    if flagv(2)
        entropyV(calcSlcIndV) = entropyV(calcSlcIndV) - ...
            sum(cooccurPatchM.*log2(cooccurPatchM+1e-10));
    end

    %     % Contrast, inverse Difference Moment, sum avg
    px = indPxM * cooccurPatchM;
    if flagv(3) || flagv(5) || flagv(6)
        %pXminusY = indCtrstM * cooccurPatchM;
        %pXplusY = indPxPlusYm * cooccurPatchM;
        contrastV(calcSlcIndV) = contrastV(calcSlcIndV) + ...
            (0:nL-1).^2 * indCtrstM * cooccurPatchM;
        invDiffMomV(calcSlcIndV) = invDiffMomV(calcSlcIndV) + ...
            1./(1+(0:nL-1).^2) * indCtrstM * cooccurPatchM;
        sumAvgV(calcSlcIndV) = sumAvgV(calcSlcIndV) + ...
            (1:2*nL) * indPxPlusYm * cooccurPatchM;
    end

    % weighted pixel average (mu), weighted pixel variance (sig)
    mu = (1:nL) * px;
    sig = bsxfun(@minus,(1:nL)',mu);
    sig = sum(sig .*sig .* px, 1);

    if flagv(4) || flagv(7) || flagv(8)
        levIMinusMu = bsxfun(@minus,levRowV',mu);
        levJMinusMu = bsxfun(@minus,levColV',mu);
    end

    % Correlation
    if flagv(4)
        corrV(calcSlcIndV) = corrV(calcSlcIndV) + ...
            sum(levIMinusMu .* levJMinusMu  .* cooccurPatchM, 1) ...
            ./ (sig + 1e-10); % sig.^2 to match ITK results (ITK bug)
    end

    % Cluster Shade
    if flagv(7)
        clstrV = levIMinusMu + levJMinusMu;
        clustShadeV(calcSlcIndV) = clustShadeV(calcSlcIndV) + ...
            sum(clstrV.*clstrV.*clstrV .* cooccurPatchM, 1);
    end
    % Cluster Prominence
    if flagv(8)
        clstrV = levIMinusMu + levJMinusMu;
        clustProminV(calcSlcIndV) = clustProminV(calcSlcIndV) + ...
            sum(clstrV.*clstrV.*clstrV.*clstrV .* cooccurPatchM, 1);
    end

    % Haralick Correlation
    if flagv(9)
        % muX = mean(px,1);
        muX = 1/nL;
        % sigX = bsxfun(@minus,px,muX);
        sigX = px - muX;
        sigX = sum(sigX .*sigX, 1)/(nL);

        % % Knuth method for mean and standard deviation (like ITK)
        %              muX = px(1,:);
        %              muPrevX = muX;
        %              sigX = muX*0;
        %              for col = 2:size(px,1)
        %                  muX = muPrevX + (px(col,:) - muPrevX)/col;
        %                  sigX = sigX + (px(col,:)-muX).*(px(col,:)-muPrevX);
        %                  muPrevX = muX;
        %              end
        %              sigX = sigX/nL;

        haralCorrV(calcSlcIndV) = haralCorrV(calcSlcIndV) + ...
            (levRowV .* levColV * cooccurPatchM - ...
            muX .* muX) ./ (sigX + eps);   % (levRowV-1) .* (levColV-1) to match ITK? Bug?

    end

    if flagv(1)
        energy3M(:,:,slcNum) = reshape(energyV(:),[numRows, numCols]);
    end
    if flagv(2)
        entropy3M(:,:,slcNum) = reshape(entropyV(:),[numRows, numCols]);
    end
    if flagv(3)
        sumAvg3M(:,:,slcNum) = reshape(sumAvgV(:),[numRows, numCols]);
    end
    if flagv(4)
        corr3M(:,:,slcNum) = reshape(corrV(:),[numRows, numCols]);
    end
    if flagv(5)
        invDiffMom3M(:,:,slcNum) = reshape(invDiffMomV(:),[numRows, numCols]);
    end
    if flagv(6)
        contrast3M(:,:,slcNum) = reshape(contrastV(:),[numRows, numCols]);
    end
    if flagv(7)
        clustShade3M(:,:,slcNum) = reshape(clustShadeV(:),[numRows, numCols]);
    end
    if flagv(8)
        clustPromin3M(:,:,slcNum) = reshape(clustProminV(:),[numRows, numCols]);
    end
    if flagv(9)
        haralCorr3M(:,:,slcNum) = reshape(haralCorrV(:),[numRows, numCols]);
    end

    if waitbarFlag
        set(hWait, 'Vertices', [[0 0 slcNum/numSlices slcNum/numSlices]' [0 1 1 0]']);
        drawnow;
    end
end
toc

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
classdef IMMOEA_F8 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP for testing IM-MOEA

%------------------------------- Reference --------------------------------
% R. Cheng, Y. Jin, K. Narukawa, and B. Sendhoff, A multiobjective
% evolutionary algorithm using Gaussian process-based inverse modeling,
% IEEE Transactions on Evolutionary Computation, 2015, 19(6): 838-856.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 3;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            t = X(:,3:obj.D).^(1./(1+3*repmat(3:obj.D,size(X,1),1)/obj.D)) - repmat(X(:,1),1,obj.D-2);
            g = sum(t.^2,2);
            PopObj(:,1) = cos(pi/2*X(:,1)).*cos(pi/2*X(:,2)).*(1+g);
            PopObj(:,2) = cos(pi/2*X(:,1)).*sin(pi/2*X(:,2)).*(1+g);
            PopObj(:,3) = sin(pi/2*X(:,1)).*(1+g);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,3);
            R = R./repmat(sqrt(sum(R.^2,2)),1,3);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            a = linspace(0,pi/2,10)';
            R = {sin(a)*cos(a'),sin(a)*sin(a'),cos(a)*ones(size(a'))};
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [refl, lambdaNm] = hyperGetEnviSignature(filename)
% HYPERGETENVISIGNATURE Reads an ENVI hyperspectral reflectance signature
%   hyperGetEnviSignature reads the RIT Target Detection Blind Test
% signature files.
%
% Usage
%   [refl, lambdaNm] = hyperGetEnviSignature(filename)
%
% Input
%   filename - Filename of signature.
% Output
%   refl - Reflectance values [0, 1].
%   lambdaNm - corresponding wavelengths in nanometers

fid = fopen(filename);

for k=1:3
    dummy = fgetl(fid);
end

num = 1;
while 1
    tmp = fgetl(fid);
    if (tmp == -1), break, end;
    v = sscanf(tmp, '%f');
    refl(num) = v(2);
    lambdaNm(num) = v(1);
    num = num + 1;
end

refl = refl / 100;

author: davidkun

repo: HyperSpectralToolbox


text:
function [ Y ] = tapas_rdcm_reduce_zeros(X, Y)
% [ Y ] = tapas_rdcm_reduce_zeros(X, Y)
%
% If there are more zero-valued frequencies than informative ones,
% subsamples those frequencies to balance dataset
%
%   Input:
%   	X           - design matrix (predictors)
%       Y           - data
%
%   Output:
%       Y           - balanced data
%

% ----------------------------------------------------------------------
%
% Authors: Stefan Fraessle (stefanf@biomed.ee.ethz.ch), Ekaterina I. Lomakina
%
% Copyright (C) 2016-2022 Translational Neuromodeling Unit
%                         Institute for Biomedical Engineering
%                         University of Zurich & ETH Zurich
%
% This file is part of the TAPAS rDCM Toolbox, which is released under the
% terms of the GNU General Public License (GPL), version 3.0 or later. You
% can redistribute and/or modify the code under the terms of the GPL. For
% further see COPYING or <http://www.gnu.org/licenses/>.
%
% Please note that this toolbox is in an early stage of development. Changes
% are likely to occur in future releases.
%
% ----------------------------------------------------------------------

% get all indices
idx = 1:size(Y,1);

% data
data = sum(abs([Y X]),2);

% zero frequencies
idx_0 = idx(data == 0);

% number of zero frequencies
n0 = sum(data==0);

% number of non-zero and non NaN frequencies
n1 = sum(data>0);

% balance the data if there are too many zeros
if ( n0 > n1 )
    idx_del = [zeros(1,n1) ones(1,n0-n1)];
    idx_del = idx_del(randperm(n0))>0;
    Y(idx_0(idx_del),:) = NaN;
end

end

author: translationalneuromodeling

repo: tapas


text:
%% RISE Tutorial by Dr. Tao Zha
%% housekeeping
clear
close all
clc
%% Instructions
% - Please run this file block by block and make sure you read the
% comments in each block to understand what it does. If there is anything
% you do not understand, ask questions to the instructor or to your
% neighbor
% - to run a particular block, click on the block and then on your keyboard
% press CTRL+Enter
%% add the paths to RISE, the data and the models
setpaths=true;
if setpaths
    addpath Models % folder with the models
    addpath Data % folder containing the data
end
%% Bring in some data and transform them into RISE's time series format (ts)

tmp=load('Data/data_nk3eq_8501_1301');  %qdatae
dataList={
    'X','output gap  y_t (log GDP_t - log GDPPotential_t)'
    'PAI','PCE core inflation pi_t (log P_t - log P_{t-1})'
    'R','FFR R_t log(1+ffr/400)  (quarterly rate already)'
    };
mydata=struct();
startdate='1985Q1';
for id=1:size(dataList,1)
    % we just give the start date, RISE automatically understand that we
    % are dealing with quarterly data by the format startdate
    mydata.(dataList{id,1})=ts(startdate,... start date
        tmp.qdatae(:,id+1),... the data
        dataList{id,2});
end

%% plot your data, compute basic statistics and look at both carefully
varlist=fieldnames(mydata);
figure('name','US data')
nvars=numel(varlist);
for id=1:nvars
    subplot(nvars,1,id)
    dd=mydata.(varlist{id});
    plot(dd,'linewidth',2)
    title(mydata.(varlist{id}).varnames)
    fprintf('%s:: mean %0.3f  stdev %0.3f\n',mydata.(varlist{id}).varnames{1},mean(dd),std(dd));
end
[~,tmp]=sup_label(['US data ',mydata.(varlist{id}).start,':', mydata.(varlist{id}).finish],'t');
set(tmp,'fontsize',15)

%% Read the model(s)

model_names={'volatilityOnly','policyOnly','volatilityPolicySame',...
    'volatilityPolicyIndependent'};
nmodels=numel(model_names);
estim_models=cell(1,nmodels);

% rather than putting all the models in the same vector as we did earlier,
% we put them in a cell array. If we put them in the same vector and call
% the estimation function, RISE will think that we want to estimate a
% pareto-type of model. But this is not what we want to do and is probably
% beyond the scope of these lectures.

% we loop through the different models using the information in the labels
for imod=1:nmodels
    % replace "for" by "parfor" if you want to use parallel computation
    estim_models{imod}=rise(model_names{imod},... % name of the file to read
        'saveas',true,... % write the expanded model to disk with the default name
        'data',mydata... % we may assign the data now or later
        );
    % a model with multiple files inserted can be difficult to read. The
    % expanded model could be useful for understanding what RISE does and
    % for debugging purposes. The expanded model contains all the details
    % of the individual files (without the comments)
end

%% We estimate the models or filter them directly
close all,clc
% if we have the parallel computing toolbox, we can estimate all models in
% one go
filtration=cell{1,nmodels};
for imod=1:nmodels
    % replace "for" by "parfor" if you want to use parallel computation
    disp('*--------------------------------------------------------------*')
    disp(['*---------Estimation of ',model_names{imod},' model-----------*'])
    disp('*--------------------------------------------------------------*')
    [estim_models{imod},filtration{imod}]=estimate(estim_models{imod},'optimizer','fmincon');
end

%% plot the smoothed probabilities
% we plot the low response (coef_2) and the high volatility (vol_2) regimes
mystates={'coef_2','vol_2'};
mylabels={'low monetary policy response regime','High volatility regime'};
for imod=1:nmodels
    mytitle=['smoothed probabilities for ',model_names{imod},' model'];
    thisstates=mystates;
    thislabels=mylabels;
    discard=false(1,numel(thisstates));
    for ii=1:numel(thisstates)
        discard(ii)=~ismember(thisstates{ii},estim_models{imod}.markov_chains.state_names);
    end
    thisstates=thisstates(~discard);
    thislabels=thislabels(~discard);
    nstates=numel(thisstates);
    figure('name',mytitle)
    for istate=1:nstates
        subplot(nstates,1,istate)
        plot(filtration{imod}.smoothed_state_probabilities.(thisstates{istate}),...
            'linewidth',2)
        title([thislabels{istate},'(chain ',thisstates{istate}(1:end-2),' state ',thisstates{istate}(end),')'])
    end
    [junk,tmp]=sup_label(mytitle,'t');
    set(tmp,'fontsize',15)
    orient tall
end

author: jmaih


text:
function[dzdxout] = roiPooling_freeform_backward(dzdx, combineFgBox)
% [dzdxout] = roiPooling_freeform_backward(dzdx, combineFgBox)
%
% Freeform pooling backward pass.
%
% If specified, this function sums the gradients for the foreground
% mask and the entire box. Otherwise it just passed through the gradients.
%
% Copyright by Holger Caesar, 2015

if combineFgBox,
    % Sum the gradients from fg and entire box
    half = size(dzdx, 3) / 2;
    dzdxFg  = dzdx(:, :, 1:half, :);
    dzdxBox = dzdx(:, :, half+1:end, :);
    dzdxout = dzdxFg + dzdxBox;
else
    % Just pass through the gradients
    dzdxout = dzdx;
end;

author: nightrome

repo: matconvnet-calvin


text:
function circle_segment_test08 ( )

%*****************************************************************************80
%
%% CIRCLE_SEGMENT_TEST08 tests CIRCLE_SEGMENT_CONTAINS_POINT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 May 2013
%
%  Author:
%
%    John Burkardt
%
  n = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CIRCLE_SEGMENT_TEST08\n' );
  fprintf ( 1, '  CIRCLE_SEGMENT_CONTAINS_POINT reports whether\n' );
  fprintf ( 1, '  a circle segment contains a point.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Pick a circle segment at random.\n' );
  fprintf ( 1, '  Compute %d sample points in the surrounding box.\n', n );
  fprintf ( 1, '  Compare the area of the segment to the percentage of points\n' );
  fprintf ( 1, '  contained in the circle segment.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       N          Omega1          Omega2           Area         Estimate\n' );
  fprintf ( 1, '\n' );

  r = 1.0;
  c = [ 0.0; 0.0 ];

  for test = 1 : 5

    [ u1, seed ] = r8_uniform_01 ( seed );
    [ u2, seed ] = r8_uniform_01 ( seed );
    omega1 = 2.0 * pi * u1;
    omega2 = 2.0 * pi * u2;

    if ( omega2 < omega1 )
      omega2 = omega2 + 2.0 * pi;
    end

    [ xy, seed ] = r8mat_uniform_01 ( 2, n, seed );
    xy = 2.0 * xy - 1.0;

    inout = zeros ( n, 1 );
    for j = 1 : n
      inout(j) = circle_segment_contains_point ( r, c, omega1, omega2, xy(1:2,j) );
    end

    theta = circle_segment_angle_from_chord_angles ( omega1, omega2 );
    area = circle_segment_area_from_angle ( r, theta );
    area_est = 4.0 * sum ( inout(1:n) ) / n;

    fprintf ( 1, '  %6d  %14g  %14g  %14g  %14g\n', n, omega1, omega2, area, area_est );

  end

  return
end


text:
function [L, f, FAIL] = linop(N)
%LINOP   Convert a CHEBOP to a LINOP.
%   L = LINOP(N) converts a CHEBOP N to a linop L if N is a linear operator. If
%   N is not linear, an error message is returned.
%
%   [L, F] = LINOP(N) returns also the affine part F of the linear CHEBOP N such
%   that L*u + F(x) = N.op(x,u).
%
%   [L, F, FAIL] = LINOP(N) will prevent an error from being thrown if N is not
%   linear, but instead return FAIL = TRUE and L = []. If N is linear, FAIL =
%   FALSE.
%
% See also LINOP, CHEBOP/LINEARIZE, CHEBOP/ISLINEAR.

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

% Tell CHEBOP/LINEARIZE() to stop if it detects nonlinearity:
linCheck = 1;

% Linearize, thereby obtaining linearity information and a LINOP:
[L, f, isLinear] = linearize(N, N.init, [], linCheck);

% We need the entire operator (including BCs) to be linear:
FAIL = ~all(isLinear);

% Throw an error is the CHEBOP is nonlinear:
if ( FAIL && (nargout < 3) )
    error('CHEBFUN:CHEBOP:linop:nonlinear',...
        'This does not appear to be a linear operator.')
end

% All values of the LINOPCONSTRAINT stored in L will be of incorrect sign when
% returned from LINEARIZE(), if we want to use it for a LINOP backslash. This is
% because when problems are solved with LINOP backslash, the solution to the
% problem is the output itself, while in a Newton iteration, we have to add the
% output of the LINOP solution to the current guess. Thus, flip the signs of the
% values of L.constraint.
L.constraint = -L.constraint;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function Y_interp = interp_missing_data(Y)
% INTERP_MISSING_DATA - interpolate missing data using linear interpolation for each pixel
%
%  [Y_INTERP] = INTERP_MISSING_DATA(Y)
%
%  Given a matrix Y with possible NaN values, this function produces
%  a sparse matrix Y_INTERP that has the linearly interpolated values
%  of each pixel that exhibits a NaN value. The values are interpolated
%  over the last dimension of Y.
%
%  Example:
%     yy = [0 0 ; 0 0];
%     yy(:,:,2) = [0 NaN ; 0 2];
%     yy(:,:,3) = [0.1 0.2 ; 0 2.3];
%     Y_interp = interp_missing_data(yy);
%     % Y_interp is a sparse matrix :    (3,2)  0.1000
%
%  See also: INTERP1

sizY = size(Y);
dimY = length(sizY);
d = prod(sizY(1:dimY-1));
T = sizY(end);
mis_data = cell(d,1);

for i = 1:d
    [ii,jj,kk] = ind2sub(sizY(1:dimY-1),i);
    if dimY == 2
        ytemp = Y(i,:);
    elseif dimY == 3
        ytemp = squeeze(Y(ii,jj,:));
    elseif dimY == 4
        ytemp = squeeze(Y(ii,jj,kk,:));
    end
    f = isnan(ytemp(:));
    y_val = interp1(find(~f),ytemp(~f),find(f),'linear','extrap');
    mis_data{i} = [i*ones(length(y_val),1),find(f(:)),y_val(:)];
end

mis_data = cell2mat(mis_data);

Y_interp = sparse(mis_data(:,1),mis_data(:,2),mis_data(:,3),d,T);

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function offset = time2offset(time, fsample)

% TIME2OFFSET converts a time-axis of a trial into the offset in samples
% according to the definition from DEFINETRIAL
%
% Use as
%   [offset] = time2offset(time, fsample)
%
% The trialdefinition "trl" is an Nx3 matrix. The first column contains
% the sample-indices of the begin of the trial relative to the begin
% of the raw data , the second column contains the sample_indices of
% the end of the trials, and the third column contains the offset of
% the trigger with respect to the trial. An offset of 0 means that
% the first sample of the trial corresponds to the trigger. A positive
% offset indicates that the first sample is later than the triger, a
% negative offset indicates a trial beginning before the trigger.

% Copyright (C) 2005, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

offset = round(time(1)*fsample);


text:
function [coeff,latent] = pca_of_hog(sbin, n)

globals;
pascal_init;

try
  load([cachedir 'pca']);
catch
  ids = textread(sprintf(VOCopts.imgsetpath, 'trainval'), '%s');
  num = length(ids);
  if nargin > 1
    num = min(n, num);
  end
  X = zeros(31, 31);
  n = 0;
  for i = 1:num
    fprintf('pca: %d/%d\n', i, num);
    rec = PASreadrecord(sprintf(VOCopts.annopath, ids{i}));
    name = [VOCopts.datadir rec.imgname];
    im = color(imread(name));
    feat = features(resize(double(im), 0.25), sbin);
    % remove occlusion feature
    feat(:,:,32) = [];
    for x = 1:size(feat,2)
      for y = 1:size(feat,1);
        v = feat(y,x,:);
        X = X + v(:) * v(:)';
        n = n+1;
      end
    end
    feat = features(resize(double(im), 0.5), sbin);
    feat(:,:,32) = [];
    for x = 1:size(feat,2)
      for y = 1:size(feat,1);
        v = feat(y,x,:);
        X = X + v(:) * v(:)';
        n = n+1;
      end
    end
    feat = features(resize(double(im), 0.75), sbin);
    feat(:,:,32) = [];
    for x = 1:size(feat,2)
      for y = 1:size(feat,1);
        v = feat(y,x,:);
        X = X + v(:) * v(:)';
        n = n+1;
      end
    end
  end

  X = X/n;
  [coeff, latent] = pcacov(X);
  save([cachedir 'pca'], 'coeff', 'latent');
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function re =maxMed(img, sz)
img = double(img);
cx = (sz + 1)/2;
cy = (sz + 1)/2;
op1 = zeros(sz);
op2 = zeros(sz);
op1(:, cx) = 1;
op2(cy, :) = 1;
op3 = diag( ones(1, sz));
op4 = fliplr(op3);
re1 = ordfilt2(img, (sz +1)/2, op1, 'symmetric');
re2 = ordfilt2(img, (sz +1)/2, op2, 'symmetric');
re3 = ordfilt2(img, (sz +1)/2, op3, 'symmetric');
re4 = ordfilt2(img, (sz +1)/2, op4, 'symmetric');
back = max( cat(3, re1, re2, re3, re4), [], 3);
re = img - back;
end

author: daxjuanxiong

repo: infrared-small-target-detection


text:
function [ a_lu, pivot, info ] = r8gb_trf ( m, n, ml, mu, a )

%*****************************************************************************80
%
%% R8GB_TRF performs a LAPACK-style PLU factorization of a R8GB matrix.
%
%  Discussion:
%
%    The R8GB storage format is for an M by N banded matrix, with lower
%    bandwidth ML and upper bandwidth MU.  Storage includes room for ML
%    extra superdiagonals, which may be required to store nonzero entries
%    generated during Gaussian elimination.
%
%    The original M by N matrix is "collapsed" downward, so that diagonals
%    become rows of the storage array, while columns are preserved.  The
%    collapsed array is logically 2*ML+MU+1 by N.
%
%    This is a simplified, standalone version of the LAPACK
%    routine R8GBTRF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 March 2004
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Anderson, Bai, Bischof, Demmel, Dongarra, Du Croz, Greenbaum,
%    Hammarling, McKenney, Ostrouchov, Sorensen,
%    LAPACK User's Guide,
%    Second Edition,
%    SIAM, 1995.
%
%  Parameters:
%
%    Input, integer M, the number of rows of the matrix A.  0 <= M.
%
%    Input, integer N, the number of columns of the matrix A.  0 <= N.
%
%    Input, integer ML, the number of subdiagonals within the band of A.
%    0 <= ML.
%
%    Input, integer MU, the number of superdiagonals within the band of A.
%    0 <= MU.
%
%    Input, real A(2*ML+MU+1,N), the matrix A in band storage.
%
%    Output, real A_LU(2*ML+MU+1,N), information about the PLU factorization.
%
%    Output, integer PIVOT(min(M,N)), the pivot indices;
%    for 1 <= i <= min(M,N), row i of the matrix was interchanged with
%    row IPIV(i).
%
%    Output, integer INFO, error flag.
%    = 0: successful exit;
%    < 0: an input argument was illegal;
%    > 0: if INFO = +i, U(i,i) is exactly zero. The factorization
%         has been completed, but the factor U is exactly
%         singular, and division by zero will occur if it is used
%         to solve a system of equations.
%
  info = 0;
  a_lu(1:2*ml+mu+1,1:n) = a(1:2*ml+mu+1,1:n);
%
%  KV is the number of superdiagonals in the factor U, allowing for fill-in.
%
  kv = mu + ml;
%
%  Set fill-in elements in columns MU+2 to KV to zero.
%
  for j = mu + 2 : min ( kv, n )
    for i = kv - j + 2 : ml
      a_lu(i,j) = 0.0E+00;
    end
  end
%
%  JU is the index of the last column affected by the current stage
%  of the factorization.
%
  ju = 1;

  for j = 1 : min ( m, n )
%
%  Set the fill-in elements in column J+KV to zero.
%
    if ( j + kv <= n )
      a_lu(1:ml,j+kv) = 0.0;
    end
%
%  Find the pivot and test for singularity.
%  KM is the number of subdiagonal elements in the current column.
%
    km = min ( ml, m-j );

    piv = abs ( a_lu(kv+1,j) );
    jp = kv + 1;

    for i = kv + 2 : kv + km + 1
      if ( piv < abs ( a_lu(i,j) ) )
        piv = abs ( a_lu(i,j) );
        jp = i;
      end
    end

    jp = jp - kv;

    pivot(j) = jp + j - 1;

    if ( a_lu(kv+jp,j) ~= 0.0E+00 )

      ju = max ( ju, min ( j+mu+jp-1, n ) );
%
%  Apply interchange to columns J to JU.
%
      if ( jp ~= 1 )
        for i = 0 : ju - j
          t = a_lu(kv+jp-i,j+i);
          a_lu(kv+jp-i,j+i) = a_lu(kv+1-i,j+i);
          a_lu(kv+1-i,j+i) = t;
        end
      end
%
%  Compute the multipliers.
%
      if ( 0 < km )

        a_lu(kv+2:kv+km+1,j) = a_lu(kv+2:kv+km+1,j) / a_lu(kv+1,j);
%
%  Update the trailing submatrix within the band.
%
        if ( j < ju )

          for k = 1 : ju-j

            if ( a_lu(kv+1-k,j+k) ~= 0.0E+00 )

              for i = 1 : km
                a_lu(kv+i+1-k,j+k) = a_lu(kv+i+1-k,j+k) ...
                  - a_lu(kv+i+1,j) * a_lu(kv+1-k,j+k);
              end
            end
          end
        end
      end

    else
%
%  If pivot is zero, set INFO to the index of the pivot
%  unless a zero pivot has already been found.
%
      if ( info == 0 )
        info = j
      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function w = minus(w,q)
   %MINUS (-) Overloaded waveform subtraction    w - q
   %   w is an N-DIMENSIONAL waveform.
   %   if q is 1x1, it is subtracted from all data values of all waveforms
   %
   %   if q is numeric, and has the same size & shape as w, then each q is
   %     subtracted from the appropriate w  ie.  result(n) = w(n) + q(n)
   %
   %   if q is a vector, it should be of same length as waveform's data
   %     sample-by-sample subtracted from each waveform
   %
   %   if q is waveform, data from q are subtracted from w, assuming both
   %   waveforms have the same number of samples
   %
   %   If q is numeric, and w is of type waveform, then the operation
   %   q - w is invalid.  This is because of a possible ambiguous result type.
   %   Instead:
   %       for a numeric result, use the synatax:
   %              q  - double(w)
   %       for a waveform result, use the syntax:
   %              q + -w             or              -w + q
   %
   % See also: minus, wavform/plus, uminus, fix_data_length

   % AUTHOR: Celso Reyes, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   if ~isa(w,'waveform')
      %yipes.  we have a number minus a waveform
      errtext = [sprintf('subtracting a <waveform> from a <%s> leads to ambiguous answer type\n',class(w)), ...
         sprintf('  for a numerical result, try: <%s> - double(<waveform>)\n',class(w)), ...
         sprintf('  for a waveform result, try:  -<waveform> + <%s>\n',class(w))];
      error('Waveform:minus:invalidClass','Error attempting: <%s> - <%s>\n%s', class(w),class(q),errtext);
   end
   try
      w = w + (-q);
   catch caughtvalue
      %NOTE, it is important to keep this section in sync with waveform/plus
      switch caughtvalue.identifier
         case 'Waveform:plus:unknownOperation'
            wsize = num2str(size(w),'%dx'); wsize = wsize(1:end-1);
            qsize = num2str(size(q),'%dx'); qsize = qsize(1:end-1);
            error('Waveform:minus:unknownOperation',...
               ['unknown waveform subtraction operation:\n',...
               '< %s %s>  - <%s %s>'],wsize, class(w), qsize, class(q));
         case 'Waveform:plus:sizeMismatch'
            error('Waveform:minus:sizeMismatch',...
               ['error in dimensions: [NxM] - [MxN].\nOne possible'...
               ' fix would be to transpose ('') one of the terms '])
         case 'Waveform:plus:invalidClass'
            error('Waveform:minus:invalidClass',...
               'unknown subtraction operation: %s - %s', class(w), class(q));
         otherwise
            rethrow(caughtvalue);
      end
   end
   return;
end

author: geoscience-community-codes

repo: GISMO


text:
function [sicd_meta] = rgazcomp_file(input_filename_cphd, output_filename_sicd, varargin)
%RGAZCOMP_FILE Applies the most simple of IFPs (range-azimuth compression)
%on a file of motion-compensated phase history data in CPHD file format and
%writes to a SICD file with appropriate metadata
%
% RGAZCOMP_FILE(INPUT_FILENAME_CPHD, OUTPUT_FILENAME_SICD, 'PropertyName', PropertyValue, ...)
%
%       Property name     Description
%       resolution        Desired resolution 3dB IPR width
%                         ([range_resolution azimuth_resolution]) in
%                         meters.  The necessary pulses and samples for
%                         image formation at this resolution will be
%                         calculated (taken from middle of data.) If
%                         pulse_range and/or sample_range properties are
%                         defined, then this parameter is ignored.  Default
%                         is the full resolution that the collect supports.
%       pulse_range       Set of pulses to use for image formation. (For
%                         example, 1:1000.) Default is to calculate this
%                         from resolution property.
%       sample_range      Set of samples to use for image formation. (For
%                         example, 1:1000.) Default is to calculate this
%                         from resolution property.
%       channel           Channel to use from CPHD file.  Default is 1.
%       sample_rate       Image domain oversample ratios.  Default is 1.
%       max_block_size    When data can't be processed within memory, this
%                         is the size of the blocks (in bytes) to use
%                         process the data.
%       quiet             If false, this reports stats on collection and
%                         IFP parameters.  Default is true.
%
% Note for larger datasets, this routine executes "out of core" storing of
% intermediate results in temporary files, which should be automatically
% cleaned up upon completion of this function.  However, if function is
% stopped in the middle, either through an error or user intervention, some
% of the files may remain in MATLAB's temporary directory (tempdir).
%
% Authors: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Parse input parameters and read metadata
ph_reader = open_ph_reader(input_filename_cphd);
ifp_params = select_pulses_samples_cphd(ph_reader, varargin{:}); % Parse generic IFP parameters
p1 = inputParser;
p1.KeepUnmatched=true;
p1.addParamValue('max_block_size',[], @isscalar);
p1.addParamValue('sample_rate',1, @isscalar);
p1.FunctionName = mfilename;
p1.parse(varargin{:});
ifp_params.sample_rate = p1.Results.sample_rate;
max_block_size = p1.Results.max_block_size;
if isempty(max_block_size)
    if ispc
        [uv, sv] = memory;
        % max_block_size is just the size of the largest single array we
        % will handle at once.  We will need significantly more memory than
        % this.  A factor of 16 seems to be reasonable.  Its a nice round
        % power of two that seems to avoid swapping in limited testing.
        max_block_size = min(sv.PhysicalMemory.Available/16, ...
            uv.MaxPossibleArrayBytes);
    else % Can't gauge memory in UNIX, so just put arbitrary value
        max_block_size = 2^27; % About 100Meg
    end
end
cphd_meta = ph_reader.get_meta();
[ignore, nbdata] = ph_reader.read_cphd(ifp_params.pulse_range, [], 1);
if ~strcmpi(cphd_meta.CollectionID.RadarMode.ModeType,'SPOTLIGHT') || ...
        any(any(diff(nbdata.SRPPos)))  % Assure spotlight data
    error('RGAZCOMP_FILE:UNSUPPORTED_COLLECT_TYPE','Unsupported collection mode.  Only spotlight data is supported.');
end
if any(diff(nbdata.SC0)) || any(diff(nbdata.SCSS))
    error('RGAZCOMP_FILE:UNSUPPORTED_COLLECT_TYPE','Unsupported data type.  Must have constant Fx0/Fx_SS.');
end
if any(diff(diff(ifp_params.sample_range)))
    error('RGAZCOMP_FILE:UNSUPPORTED_COLLECT_TYPE','Selected samples must be regularly spaced.');
end
if ~strcmp(cphd_meta.Global.DomainType, 'FX')
    error('RGAZCOMP_FILE:UNSUPPORTED_FILE_TYPE','Unsupported CPHD type.  Currently only FX domain is supported.');
end

%% Compute some values we will need
% Number of pulses and samples that we will be processing in this function.
% Defined here purely for conciseness of notation.
num_pulses  = length(ifp_params.pulse_range);
num_samples = length(ifp_params.sample_range);

% Decide whether we can process in memory, or whether we need to use
% intermediate files.  We hope to avoid using intermediate files since file
% reads and writes take the great majority of the processing time here.  In
% one profiling test, file reads and writes took 98% of the time of this
% funcion.  For this reason, we have attempted to minimize the number of
% reads and writes, and we try to comment in the code sections that affect
% these file I/O delays.  Unfortunately virtual memory and swapping (which
% is non-optimal file I/O out of our control) is much worse than the file
% I/O explicitly handled in this function, so try to set max_block_size to
% something as large as possible without the operating system resorting to
% virtual memory.
data_element_size = 8; % In bytes.  Assumes single precision data
process_in_blocks = max_block_size<(num_pulses*num_samples*prod(ifp_params.sample_rate)*data_element_size);
if process_in_blocks
    warning('RGAZCOMP_FILE:FILE_PROCESSING','This data must be processed through intermediate files, rather than all in memory.  This could take a while...');
end

ifp_params.image_size_pixels = floor(ifp_params.sample_rate .* [num_samples num_pulses]);
sicd_meta = rgazcomp_sicd_meta(cphd_meta,nbdata,ifp_params);

%% Range FFT across samples
% The FFT will be done in separably in two steps, first in range and then
% in azimuth.

% Setup temporary file
% The range FFT stage includes an in-order read (CPHD is written in
% "pulse-major" order) and an in-order write (data interpolated in range).
% We have to "corner-turn" (or transpose) after the range FFT (as we do on
% all 4 steps), but it is typically faster to read from a file out-of-order
% than write out-of-order, so we will do the corner-turn in the read of the
% azimuth FFT, rather than the write of this range FFT stage.
if process_in_blocks
    filename_range_fft = tempname;
    tempsicdmeta.ImageData.NumRows   = num_pulses;
    tempsicdmeta.ImageData.NumCols   = ifp_params.image_size_pixels(1);
    writer_range_interp = SIOWriter(filename_range_fft, tempsicdmeta, 'include_sicd_metadata', false);
end

% Iterate through sets of pulses for range FFT
first_pulse_in_set=1; % Out of all pulses selected for processing
max_num_lines = floor(max_block_size/(num_samples*data_element_size));
wb_hand=waitbar(0,'Range FFT...');
while(first_pulse_in_set<num_pulses)
    waitbar((first_pulse_in_set-1)/num_pulses,wb_hand);
    pulse_indices = first_pulse_in_set:...
        min(num_pulses,first_pulse_in_set+max_num_lines-1);
    data_block = ph_reader.read_cphd(ifp_params.pulse_range(pulse_indices),...
        ifp_params.sample_range, ifp_params.channel);

    data_block = pfa_fft_zeropad_1d(data_block,ifp_params.sample_rate);

    if process_in_blocks
        writer_range_interp.write_chip(data_block,[1 first_pulse_in_set]); % Save range interpolated data in this block
    end
    first_pulse_in_set=pulse_indices(end)+1;
end
if process_in_blocks
    clear writer_range_interp % Close file
end
ph_reader.close();

%% Azimuth FFT across vectors
if process_in_blocks
    reader_range_fft = open_reader(filename_range_fft);
end
writer_final_sicd = SICDWriter(output_filename_sicd, sicd_meta); % Final image will be written here

% Iterate through blocks of data for azimuth FFT
first_line_in_block=1; % Out of all rows in the final image
max_num_lines = floor(max_block_size/(num_pulses*data_element_size));
waitbar(0,wb_hand,'Azimuth FFT...');
while(first_line_in_block<ifp_params.image_size_pixels(1))
    waitbar((first_line_in_block-1)/ifp_params.image_size_pixels(1),wb_hand);
    line_range = [first_line_in_block ...
        min(ifp_params.image_size_pixels(1),first_line_in_block+max_num_lines-1)];
    if process_in_blocks
        data_block = reader_range_fft.read_chip(line_range, [1 num_pulses]);
    end

    data_block = pfa_fft_zeropad_1d(data_block.',ifp_params.sample_rate);
    if isfield(sicd_meta,'SCPCOA') && isfield(sicd_meta.SCPCOA,'SideOfTrack') && ...
            isequal(sicd_meta.SCPCOA.SideOfTrack, 'R')
        data_block = flipud(data_block); % View from above
    end

    writer_final_sicd.write_chip(data_block,[1 first_line_in_block]);
    first_line_in_block=line_range(2)+1;
end
clear writer_final_sicd; % Close file
if process_in_blocks
    reader_range_fft.close(); % Close file
    delete(filename_range_fft); % Intermediate file no longer needed
end

close(wb_hand); % Finally, close waitbar

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint