########################################################################
##
## Copyright (C) 1995-2021 The Octave Project Developers
##
## See the file COPYRIGHT.md in the top-level directory of this
## distribution or .
##
## This file is part of Octave.
##
## Octave 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.
##
## Octave 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 Octave; see the file COPYING. If not, see
## .
##
########################################################################
## -*- texinfo -*-
## @deftypefn {} {@var{B} =} fillmissing (@var{A}, 'constant', @var{v})
## @deftypefnx {} {@var{B} =} fillmissing (@var{A}, @var{method})
## @deftypefnx {} {@var{B} =} fillmissing (@var{A}, @var{move_method}, @var{window_size})
## @deftypefnx {} {@var{B} =} fillmissing (@var{A}, @var{fill_function}, @var{window_size})
## @deftypefnx {} {@var{B} =} fillmissing (@dots{}, @var{dim})
## @deftypefnx {} {@var{B} =} fillmissing (@dots{}, @var{PropertyName}, @var{PropertyValue})
## @deftypefnx {} {[@var{B}, @var{idx}] =} fillmissing (@dots{})
##
## Replace missing entries of array @var{A} either with values in @var{v} or
## as determined by other specified methods. 'missing' values are determined
## by the data type of @var{A} as identified by the function @ref{ismissing},
## curently defined as :
##
## @itemize
## @item
## @code{NaN}: @code{single}, @code{double}
##
## @item
## @code{' '} (white space): @code{char}
##
## @item
## @code{@{''@}} (white space in cell): string cells.
## @end itemize
##
## @var{A} can be a numeric scalar or array, a character vector or array, or
## a cell array of character vectors (a.k.a. string cells).
##
## @var{v} can be a scalar or an array containing values for replacing the
## missing values in @var{A} with a compatible data type for isertion into
## @var{A}. The shape of @var{v} must be a scalar or an array with number
## of elements in @var{v} equal to the number of elements orthoganal to the
## operating dimension. E.g., if @code{size(@var{A})} = [3 5 4], operating
## along @code{dim} = 2 requires @var{v} to contain either 1 or 3x4=12
## elements.
##
## If requested, the optional output @var{idx} will contain a logical array
## the same shape as @var{A} indicating with 1's and 0's which locations in
## @var{A} were and were not filled.
##
## Alternate Input Arguments and Values:
## @itemize
## @item @var{method} - replace missing values with:
## @table @code
##
## @item next
## @itemx previous
## @itemx nearest
## next, previous, or nearest non-missing value (nearest defaults to next
## when equidistant as determined by 'SamplePoints'.)
##
## @item linear
## linear interpolation of neigboring, non-missing values
##
## @item spline
## piecewise cubic spline interpolation of neigboring, non-missing values
##
## @item pchip
## 'shape preserving' piecewise cubic spline interposaliton of neighboring,
## non-missing values
## @end table
##
## @item @var{move_method} - moving window calculated replacement values:
## @table @code
##
## @item movmean
## @itemx movmedian
## moving average or median using a window determined by @var{window_size}.
## @var{window_size} must be either a positive integer scalar value, or a
## two element positive integer vector. For scalar values, the window is
## centered on the missing element and contains @var{window_size}/2 elements
## on either side of the missing element. If @var{window_size} is even, it is
## centered on the current missing element and the previous element. If
## @var{window_size} is a vector, the first and second elements represent
## the number of elements backward and forward, respectively, from the current
## missing element to be included in the window.
## @end table
##
## @item @var{fill_function} - custom method specified as a function handle
## Supplied function must accept three inputs:
## @table @var
## @item xval -
## reference data values near the missing value
## @item xloc -
## location of the reference data
## @item missingloc -
## location of the missing data that needs to be filled
## @end table
##
## Supplied function must return a scalar or vector with the same number of
## elements as @var{missingloc}. The required @var{window_size} parameter
## must be either a positive integer scalar value, or a two element positive
## integer vector. @var{window_size} values are assumed to have the same
## scale as the sample points (and @var{xloc} and @var{missingloc}). For
## scalar values, the window is centered on the missing element and contains
## elements extending a distance of @var{window_size}/2 units to either side
## of the element gap. If @var{window_size} is a vector, the first and second
## elements represent the distance backward and forward, respectively, that
## the window extends from the current element gap to be included in the
## window.
##
## @item @var{dim} - specify a dimension for vector operations to operate
## along (default = first non-singeton dimension)
##
## @item @var{PropertyName}-@var{PropertyValue} pairs
## @table @code
## @item SamplePoints
## @var{PropertyValue} is a vector of sample point values representing the
## sorted and unique x-axis values of the data in @var{A}. If unspecified,
## the default is assumed to be the vector @var{[1 : size (A, dim)]}.
##
## @item EndValues
## Special handling method for missing values in first and last positions.
## @var{PropertyValue} can be a constant scalar value to use for filling, or
## one of the following options:
## @table @code
## @item extrap
## use the same process as @var{method}
## @item previous
## @itemx next
## @itemx nearest
## use the previous, next, or nearest non-missing value
## @item none
## do not fill the end values
## @end table
##
## @item MissingLocations
## @var{PropertyValue} must be a logical array the same size as @var{A}
## indicating locations of known missing data with a value of @code{true}.
## (cannot be combined with MaxGap)
##
## @item MaxGap
## @var{PropertyValue} is a numeric scalar indicating the maximum gap length
## to fill, and assumes the same distance scale as the sample points. Gap
## length is calculated by the difference in locations of the sample points
## on either side of the gap, and gaps larger than MaxGap are ignored by
## @var{fillmissing}. (cannot be combined with MissingLocations)
## @end table
## @end itemize
##
## Compatibility Notes:
## @itemize
## @item Numerical and logical inputs for @var{A} and @var{v} may be specified
## in any combination. The output will be the same class as @var{A}, with the
## @var{v} converted to that data type for filling. Only @code{single} and
## @code{double} have defined 'missing' values, so @var{A} of other data types
## will always output @var{B} = @var{A}.
##
## @item Method "makima" is not yet implemented in @code{interp1}, which is
## used by @code{fillmissing}.
## @end itemize
## @end deftypefn
##
## @seealso{ismissing, rmmissing, standardizeMissing}
function [A, idx_out] = fillmissing (A, varargin)
if (nargin < 2)|| (nargin > 12)
print_usage ();
endif
method = varargin{1};
if (ischar (method))
method = lower (method);
elseif (! is_function_handle (method))
error ("fillmissing: second input must be a string or function handle");
endif
input_class = class (A);
sz_A = size (A);
ndims_A = numel (sz_A);
samplepoints = [];
maxgap = [];
missing_locs = [];
endgap_method = [];
endgap_locs = [];
v = [];
endgap_val = [];
dim = [];
reshape_flag = false;
## process input arguments
if is_function_handle (method)
## error ("fillmissing: function handle methods not yet implemented");
## verify function handle and window
if ((nargin < 3) || ! isnumeric (varargin{2}) || ...
! any( numel (varargin{2})==[1 2]))
error (["fillmissing: fill function handle must be followed by ", ...
"a numeric scalar or two-element vector window size"]);
elseif (nargin (method) < 3)
error ("fillmissing: fill function must accept at least three inputs.");
endif
window_size = varargin{2};
next_varg = 3;
else
switch (method)
case {"previous", "next", "nearest"}
next_varg = 2;
case {"linear", "spline", "pchip", "makima"}
next_varg = 2;
##TODO - interpolant methods need to verify numeric only
case "constant"
if ((nargin < 3))
error (["fillmissing: 'constant' method must be followed by ", ...
"a numeric scalar or array"]);
endif
v = varargin{2}(:);
if isempty(v)
error ("fillmissing: the fill value cannot be an emtpy array");
endif
## type check v against A
if (iscellstr (A) && ischar (v) && ! iscellstr (v))
v = {v};
endif
if ((! isempty (v)) && ...
((isnumeric (A) && ! (isnumeric (v) || islogical (v))) || ...
(ischar (A) && ! ischar (v)) || ...
(iscellstr (A) && ! (iscellstr (v)))))
error ("fillmissing: fill value must have the same data type as 'A'");
endif
## v can't be size checked until after processing rest of inputs
next_varg = 3;
case {"movemean", "movemedian"}
if ((nargin < 3) || ! isnumeric (varargin{2}) || ...
! any( numel (varargin{2})==[1 2]))
error (["fillmissing: moving window method must be followed by ", ...
"a numeric scalar or two-element vector"]);
endif
window_size = varargin{2};
next_varg = 3;
otherwise
error ("fillmissing: unknown fill method '%s'", method);
endswitch
endif
## process any more parameters
if (next_varg < nargin)
#set dim. if specified, it is only numeric option that can appear next
if isnumeric (varargin{next_varg})
dim = varargin{next_varg};
if (!((isscalar (dim)) && (dim > 0)))
error ("fillmissing: DIM must be a positive scalar");
endif
next_varg++;
else
## default dim is first nonsingleton dimension of A
if isscalar (A)
dim = 1;
else
dim = find (sz_A > 1, 1, "first");
endif
endif
sz_A_dim = size (A, dim);
## process any remaining inputs, must be name-value pairs
while (next_varg < nargin)
propname = varargin{next_varg};
if (next_varg + 1 == nargin)
## must be at least one more input with 1st containing value
error ("fillmissing: properties must be given as name,value pairs");
else
propval = varargin{next_varg + 1};
next_varg = next_varg + 2;
if (! ischar (propname))
error ("invalid parameter name specified");
else
propname = lower (propname);
endif
## input validation for names and values
switch (propname)
case "samplepoints"
## val must be sorted, unique, numeric vector the same size
## as size(A,dim)
if (! (isnumeric (propval) && isvector (propval)
&& (numel (propval) == sz_A_dim) && issorted (propval)
&& (numel (propval) == numel (unique (propval)))))
error (["fillmissing: SamplePoints must be a sorted ", ...
"non-repeating, numeric vector with %d elements"], ...
sz_A_dim);
endif
samplepoints = propval(:);
case "endvalues"
## for numeric A, val must be numeric scalar, a numeric
## array with numel equal to the elements orthogonal to
## the dim or certain string methads. For non-numeric A,
## "constant" method is not valid.
if ischar (propval)
switch (lower (propval))
case {"extrap", "previous", "next", "nearest", "none"}
endgap_method = propval;
otherwise
error ("fillmissing: invalid EndValues method '%s'", propval);
endswitch
elseif (isnumeric (propval))
if (! (isnumeric (A) || islogical (A)))
error (["fillmissing: EndValues method 'constant' only ", ...
"valid for numeric arrays."]);
endif
endgap_method = 'constant';
endgap_val = propval;
else
error (["fillmissing: EndValues must be numeric or a ", ...
"valid method name"]);
endif
case "missinglocations"
## val must be logical array same size as A
if (! (islogical (propval) && isequal (sz_A, size (propval))))
error (["fillmissing: MissingLocations must be a logical ", ...
"array the same size as A"]);
endif
if (! isempty (maxgap))
error (["fillmissing: MissingLocations and MaxGap methods ", ...
"cannot be used simultaneously"]);
endif
missing_locs = propval;
case "maxgap"
## val must be positive numeric scalar
if (! (isnumeric (propval) && isscalar (propval) && (propval > 0)))
error ("fillmissing: MaxGap must be a positive numeric scalar");
endif
if (! isempty (missing_locs))
error (["fillmissing: MissingLocations and MaxGap methods ", ...
"cannot be used simultaneously"]);
endif
maxgap = propval;
otherwise
error ("invalid parameter name '%s'", propname);
endswitch
endif
endwhile
else
## no inputs after method
## set default dim
if isscalar (A)
dim = 1;
else
dim = find (sz_A > 1, 1, "first");
endif
sz_A_dim = size (A, dim);
endif
## reduce calls to size and avoid overruns checking sz_A for high dims
if dim > ndims_A
sz_A = [sz_A, ones(1, dim - numel(sz_A))];
ndims_A = numel (sz_A);
endif
## set defaults for any unspecified parameters
if (isempty (samplepoints))
samplepoints = [1 : sz_A_dim]';
endif
if isempty (missing_locs)
missing_locs = ismissing (A);
endif
## endvalues treated separately from interior missing_locs
if (isempty (endgap_method) || strcmp (endgap_method, "extrap"))
endgap_method = method;
if strcmp(endgap_method, "constant")
endgap_val = v;
endif
endif
## verify size of v and endgap_val for 'constant' methods, resize for A
orthogonal_size = [sz_A(1:dim-1),1,sz_A(dim+1:end)]; ## orthog. to dim size
numel_orthogonal = prod (orthogonal_size); ## numel perpen. to dim
if (strcmp (method, "constant") && (! isscalar (v)))
if (numel (v) != numel_orthogonal)
error (["fillmissing: fill value 'V' must be a scalar or a ", ...
"%d element array"], numel_orthogonal);
else
v = reshape (v, orthogonal_size);
endif
endif
if (strcmp (endgap_method, "constant") && (! isscalar (endgap_val)))
if (numel (endgap_val) != numel_orthogonal)
error (["fillmissing: EndValues must be a scalar or a ", ...
"%d element array"], numel_orthogonal);
else
endgap_val = reshape (endgap_val, orthogonal_size);
endif
endif
## simplify processing by temporarily permuting A so operation always on dim1
## revert permutation at the end
dim_idx_perm = [1 : ndims_A];
dim_idx_flip(1 : max(dim, ndims_A)) = {':'};
dim_idx_flip(1) = [sz_A_dim:-1:1];
if (dim != 1)
dim_idx_perm([1, dim]) = [dim, 1];
A = permute (A, dim_idx_perm);
sz_A([1, dim]) = sz_A([dim, 1]);
missing_locs = permute (missing_locs, dim_idx_perm);
reshape_flag = true;
orthogonal_size = [1, sz_A(2:end)];
if (! isempty (v))
v = permute (v, dim_idx_perm);
endif
if (! isempty (endgap_val))
endgap_val = permute (endgap_val, dim_idx_perm);
endif
endif
## precalculate fill data for several methods
zero_padding = zeros (orthogonal_size);
samplepoints_expand = samplepoints(:, ones (1, prod (sz_A(2:end))));
##find endgap locations
if (sz_A_dim < 3)
## all missing are endgaps
endgap_locs = missing_locs;
else
## use cumsums starting from first and last part in dim to find missing
## values in and adjacent to end locations.
endgap_locs = cumprod (missing_locs,1) | ...
(cumprod (missing_locs(dim_idx_flip{:}),1))(dim_idx_flip{:});
endif
## remove endgap_locs from missing_locs to avoid double processing
missing_locs(endgap_locs) = false;
## remove elements from missing and end location arrays if maxgap is specified
if (! isempty (maxgap))
## missing_locs: if samplepoints value diff on either side of missing
## elements is > maxgap, remove those values.
## for endgaps, use diff of inside and missing end samplepoint values
## and remove from endgaps
## First check gapsize of any interior missings in missing_locs
if (any (missing_locs(:)))
## locations in front of gaps
loc_first = [diff(missing_locs,1,1); zero_padding] == 1;
## locations in back of gaps
loc_last = diff ([zero_padding; missing_locs],1,1) == -1;
## value of samplepoints at front and back locs
samplepointvals_first = samplepoints_expand(loc_first);
samplepointvals_last = samplepoints_expand(loc_last);
## evaluate which gaps are too big to fill
gaps_to_remove = (samplepointvals_last - samplepointvals_first) > maxgap;
## convert those gaps into an array element list
idxs_to_remove = arrayfun ('colon', ...
((find (loc_first))(gaps_to_remove ) + 1), ...
((find (loc_last))(gaps_to_remove ) - 1), ...
"UniformOutput", false);
## remove those elements from missing_locs
missing_locs([idxs_to_remove{:}]) = false;
endif
##then do any endgaps
if (any (endgap_locs(:)))
## if any are all missing, remove for any value of maxgap
endgap_locs &= ! prod (endgap_locs, 1);
if ((sz_A_dim < 3) && (abs (samplepoints(2) - samplepoints(1)) > maxgap))
## shortcut - all missings are ends and exceed maxgap.
endgap_locs(:) = false;
else
## check gap size of front endgaps
##find loc element after gap
nextvals = sum (cumprod (endgap_locs,1)) + 1;
## compare diff between values at those points and at base with maxgap
ends_to_remove = abs (samplepoints(nextvals) - samplepoints(1)) ...
> maxgap;
## remove any with gap>maxgap
endgap_locs((cumprod (endgap_locs,1)) & ends_to_remove) = false;
## flip, repeat for back endgaps, then unflip and remove.
nextvals = sum (cumprod (endgap_locs(dim_idx_flip{:}),1)) + 1;
ends_to_remove = abs (samplepoints(end:-1:1)(nextvals) ...
- samplepoints(end)) > maxgap;
endgap_locs((cumprod (...
endgap_locs(dim_idx_flip{:}), 1)(dim_idx_flip{:})) & ...
ends_to_remove) = false;
endif
endif
endif
## fill the missing data
if is_function_handle (method)
%%dosomething
error ("fillmissing: function handle methods not yet implemented");
else
## process central missing values (gaps bound by two valid datapoints)
if (any (missing_locs(:)))
switch (method)
case "constant"
if (isscalar (v))
A(missing_locs) = v;
else
A(missing_locs) = (missing_locs .* v)(missing_locs);
endif
case {"previous", "next", "nearest", "linear"}
## interp1 should be able to provide all of these functions, but
## it cannot work with only 1 valid datapoint, but fillmissing needs
## to work for those. So must manually process these options.
## find element locations before and after each gap
loc_first = [diff(missing_locs, 1, 1); zero_padding] == 1;
loc_last = diff ([zero_padding; missing_locs], 1, 1) == -1;
gapsizes = find (loc_last) - find (loc_first) - 1;
## find prev, next, nearest fill values and replace in A.
switch (method)
case "previous"
fill_vals = repelems (A(loc_first), ...
[1:numel(gapsizes); gapsizes'])';
case "next"
fill_vals = repelems (A(loc_last), ...
[1:numel(gapsizes); gapsizes'])';
case {"nearest", "linear"}
## determine which missings go with values before or after
## gap based on samplevalue distance. (equal dist goes to after)
## find sample values before and after gaps
samplepointvals_first = samplepoints_expand(loc_first);
samplepointvals_last = samplepoints_expand(loc_last);
## build cell with linear indices of elements in each gap
gap_locations = arrayfun ('colon', ...
((find (loc_first)) + 1), ...
((find (loc_last)) - 1), ...
"UniformOutput", false);
## get sample values at those elements
[samplevals_at_gapvals, ~] = ind2sub (sz_A, [gap_locations{:}]);
samplevals_at_gapvals = samplepoints(samplevals_at_gapvals);
switch (method)
case "nearest"
gap_mids = (samplepointvals_last - ...
samplepointvals_first)/2 + samplepointvals_first;
## build comparison vector
sampleval_comp = repelems (gap_mids, ...
[1:numel(gapsizes);(gapsizes')])';
## determine which ones are replaced with values before or
## after to generate fill array
prev_fill = repelems (A(loc_first), ...
[1:numel(gapsizes); gapsizes'])' ...
.* (samplevals_at_gapvals < sampleval_comp);
next_fill = repelems (A(loc_last), ...
[1:numel(gapsizes); gapsizes'])' ...
.* (samplevals_at_gapvals >= sampleval_comp);
## build fill vector
fill_vals = next_fill + prev_fill;
case {"linear"}
## middle gaps already sorted to guarantee bounding
## point front and back
## expand first and last vectors for each gap point
Avalues_first = repelems (A(loc_first), ...
[1:numel(gapsizes); gapsizes'])';
Avalues_last = repelems (A(loc_last), ...
[1:numel(gapsizes); gapsizes'])';
## expand samplepoint values for interpolation x-values
samplepointvals_first = repelems (samplepointvals_first, ...
[1:numel(gapsizes);(gapsizes')])';
samplepointvals_last = repelems (samplepointvals_last, ...
[1:numel(gapsizes);(gapsizes')])';
## linearly interpolate
interp_slopes = (Avalues_last - Avalues_first) ...
./ (samplepointvals_last - samplepointvals_first);
fill_vals = interp_slopes .* (samplevals_at_gapvals ...
- samplepointvals_first) + Avalues_first;
endswitch
endswitch
## replace missings with appropriate fill values
A(missing_locs) = fill_vals;
case {"spline", "pchip", "makima"}
## use interp1 to perform more complex interpolations
## TODO: 'linear' is ~10-100x faster than this columnwise interp1
## approach. look to speed these up as well.
## previous processing guarantees missing_locs gaps bound by 2 points.
## determine columns needing operation to reduce empty operations
cols_to_use = any (missing_locs, 1);
## find logical data and empty location indices for each column to be
## used in columnwise call to interp1 (cast as columnwise cell arrays)
interp_locations = num2cell (...
(! (missing_locs | endgap_locs))(:, cols_to_use), 1);
interp_fill_locs = num2cell (...
(missing_locs & cols_to_use)(:, cols_to_use), 1);
## use locations to build cell arrays with sample and interpolant
## values to pass to interp1
interp_loc_ref = num2cell (struct ("type", "()", "subs", ...
num2cell (interp_locations)));
A_interpvalues = cellfun ('subsref', num2cell ( ...
A(:, cols_to_use), 1), interp_loc_ref, ...
"UniformOutput", false);
interp_samplevals = cellfun ('subsref', {samplepoints}, ...
interp_loc_ref, "UniformOutput", false);
interp_empty_samplevals = cellfun ('subsref', {samplepoints}, ...
num2cell (struct ("type", "()", "subs", ...
num2cell(interp_fill_locs))), "UniformOutput", false);
## generate fill_vals using interp1 for missing locations.
fill_vals = vertcat(cellfun ('interp1', interp_samplevals, ...
A_interpvalues, interp_empty_samplevals, {method}, ...
'UniformOutput', false){:});
## replace missings with appropriate fill values
A(missing_locs) = fill_vals;
endswitch
endif
## process endgaps
if (any (endgap_locs(:)))
switch (endgap_method)
case "none"
endgap_locs(:) = false;
case "constant"
if (isscalar (endgap_val))
A(endgap_locs) = endgap_val;
else
#reshape to size A trimmed, then do permute
A(endgap_locs) = (endgap_locs .* endgap_val)(endgap_locs);
endif
case {"previous", "next", "nearest", "linear", ...
"spline", "pchip", "makima"}
## all of these methods require sz_A_dim >= 2. shortcut path otherwise
if sz_A_dim < 2
endgap_locs(:) = false;
else
## index offset array for linear indexincg
offset = reshape ([0 : sz_A_dim : ...
sz_A_dim * (numel_orthogonal-1)], orthogonal_size);
switch (endgap_method)
case "previous"
## remove any points at front of array, includes all-NaN columns
endgap_locs (logical (cumprod (endgap_locs,1))) = false;
##find locations of the 'prev' value to use for filling
subsval_loc = [diff(endgap_locs, 1, 1); zero_padding] == 1;
##calculate the number of spots each 'prev' needs to fill
gapsizes = sz_A_dim - (find (subsval_loc) - ...
(subsval_loc .* offset)(subsval_loc));
## construct substitution value vector
fill_vals = repelems (A(subsval_loc), ...
[1:numel(gapsizes); gapsizes'])';
## replace missings with appropriate fill values
case "next"
##take any points at back of array and remove from endgap_locs
## includes any all-NaN columns
endgap_locs(logical (cumprod ( ...
endgap_locs(dim_idx_flip{:}), 1)(dim_idx_flip{:}))) ...
= false;
##find locations of the 'next' value to use for filling
subsval_loc = diff ([zero_padding; endgap_locs],1,1) == -1;
##calculate the number of spots each 'next' needs to fill
gapsizes = (find (subsval_loc) - ...
(subsval_loc .* offset)(subsval_loc)) - 1;
## construct substitution value vector
fill_vals = repelems (A(subsval_loc), ...
[1:numel(gapsizes); gapsizes'])';
case "nearest"
## remove any all-missing columns
endgap_locs &= (! prod (endgap_locs, 1));
## find front end info
front_gap_locs = logical (cumprod (endgap_locs, 1));
front_next_loc = diff ( ...
[zero_padding; front_gap_locs], 1, 1) == -1;
front_gapsizes = (find (front_next_loc) - ...
(front_next_loc .* offset)(front_next_loc)) - 1;
## find back end info
back_gap_locs = logical ( ...
cumprod (endgap_locs(dim_idx_flip{:}), 1)(dim_idx_flip{:}));
back_prev_loc = [diff(back_gap_locs, 1, 1); zero_padding] == 1;
back_gapsizes = sz_A_dim - (find (back_prev_loc) - ...
(back_prev_loc .* offset)(back_prev_loc));
## combine into fill variables
[~, fb_sort_idx] = sort ...
([find(front_gap_locs);find(back_gap_locs)]);
subsval_loc = [find(front_next_loc); find(back_prev_loc)];
gapsizes = [front_gapsizes; back_gapsizes];
## construct substitution value vector with sort order to mix
## fronts and backs in column order
fill_vals = (repelems (A(subsval_loc), ...
[1:numel(gapsizes); gapsizes'])')(fb_sort_idx);
case "linear"
## endgaps not guaranteed to have enough points to interpolate
cols_to_use = (sum (! (missing_locs | endgap_locs), 1) > 1) ...
& any (endgap_locs, 1);
interp_locs = ! (missing_locs | endgap_locs) & cols_to_use;
endgap_locs &= cols_to_use;
## process front endgaps
front_gap_locs = logical (cumprod (endgap_locs, 1));
fill_vals_front = [];
if (any (front_gap_locs(:)))
front_gapsizes = nonzeros (sum (front_gap_locs, 1));
## collect first data point after gap & expand to gapsize
front_interppoint_1 = repelems ( ...
find (diff ([zero_padding; front_gap_locs], 1, 1) == -1), ...
[1:numel(front_gapsizes); front_gapsizes'])';
## collect second data point after gap & expand to gapsize
front_interppoint_2 = repelems ( ...
find (diff ([zero_padding; ((cumsum (interp_locs, 1) .* ...
any (front_gap_locs, 1)) == 2)], 1, 1) == 1), ...
[1:numel(front_gapsizes); front_gapsizes'])';
front_interp_Avals = A([front_interppoint_1, ...
front_interppoint_2]);
front_interp_sampvals = samplepoints_expand( ...
[front_interppoint_1, front_interppoint_2]);
front_gap_loc_sampvals = samplepoints_expand(front_gap_locs);
## hack for vector automatic orientation forcing col vector
if (isvector (front_interp_Avals))
front_interp_Avals = (front_interp_Avals (:)).';
endif
if (isvector (front_interp_sampvals))
front_interp_sampvals = (front_interp_sampvals (:)).';
endif
## perform interpolation for every gap point
interp_slopes_front = diff (front_interp_Avals, 1, 2) ...
./ diff (front_interp_sampvals, 1, 2);
fill_vals_front = interp_slopes_front .* ...
(front_gap_loc_sampvals - ...
front_interp_sampvals(:,1)) + ...
front_interp_Avals(:,1);
endif
## process back endgaps
back_gap_locs = logical ( ...
cumprod (endgap_locs(dim_idx_flip{:}), 1)(dim_idx_flip{:}));
fill_vals_back = [];
if (any (back_gap_locs(:)))
back_gapsizes = ...
nonzeros (sum (back_gap_locs(dim_idx_flip{:}), 1));
## collect last data point before gap & expand to gapsize
back_interppoint_2 = repelems ( ...
find ([diff(back_gap_locs, 1, 1); zero_padding] == 1), ...
[1:numel(back_gapsizes); back_gapsizes'])';
## collect 2nd to last data point before gap & expand to gap
back_interppoint_1 = repelems ( ...
find ((diff ([zero_padding; ...
((cumsum (interp_locs(dim_idx_flip{:}), 1) .* ...
any (back_gap_locs, 1)) == 2)], ...
1, 1) == 1)(dim_idx_flip{:})), ...
[1:numel(back_gapsizes); back_gapsizes'])';
## build linear interpolant vectors
back_interp_Avals = A([back_interppoint_1, ...
back_interppoint_2]);
back_interp_sampvals = samplepoints_expand( ...
[back_interppoint_1, ...
back_interppoint_2]);
back_gap_loc_sampvals = samplepoints_expand(back_gap_locs);
## hack for vector automatic orientation forcing col vector
if (isvector (back_interp_Avals))
back_interp_Avals = (back_interp_Avals (:)).';
endif
if (isvector (back_interp_sampvals))
back_interp_sampvals = (back_interp_sampvals (:)).';
endif
## perform interpolation for every gap point
interp_slopes_back = diff (back_interp_Avals, 1, 2) ...
./ diff (back_interp_sampvals, 1, 2);
fill_vals_back = interp_slopes_back .* ...
(back_gap_loc_sampvals - ...
back_interp_sampvals(:,1)) + ...
back_interp_Avals(:,1);
endif
[~, fb_sort_idx] = sort ...
([find(front_gap_locs); find(back_gap_locs)]);
fill_vals = [fill_vals_front; fill_vals_back](fb_sort_idx);
case {"spline", "pchip", "makima"}
## use interp1 to perform more complex interpolations
## TODO: split out from columnwise interp1 approach to speed up
## endgap_locs not guaranteed to have 2 points.
## need to ignore columns with < 2 values, or with nothing to do
## find columns needing operation to reduce empty operations
cols_to_use = (sum (! (endgap_locs | missing_locs), 1) > 1) ...
& any (endgap_locs, 1);
## trim out unused cols from endgap_locs
endgap_locs &= cols_to_use;
## find logical data and empty location indices for each column
## to be used in columnwise call to interp1 (cast as columnwise
## cell arrays)
interp_locations = num2cell ...
((! (missing_locs | endgap_locs))(:, cols_to_use), 1);
interp_fill_locs = num2cell ...
(endgap_locs(:, cols_to_use), 1);
## use locations to build cell arrays with sample and
## interpolant values to pass to interp1
interp_loc_ref = num2cell (struct ("type", "()", "subs", ...
num2cell (interp_locations)));
A_interpvalues = cellfun ('subsref', num2cell (...
A(:, cols_to_use), 1), interp_loc_ref, ...
"UniformOutput", false);
interp_samplevals = cellfun ('subsref', {samplepoints}, ...
interp_loc_ref, "UniformOutput", false);
interp_empty_samplevals = cellfun ('subsref', ...
{samplepoints}, ...
num2cell (struct ("type", "()", "subs", ...
num2cell(interp_fill_locs))), ...
"UniformOutput", false);
## generate fill_vals using interp1 for missing locations.
fill_vals = vertcat(cellfun ('interp1', interp_samplevals, ...
A_interpvalues, interp_empty_samplevals, ...
{method}, {'extrap'}, 'UniformOutput', false){:});
otherwise
error ("EndValues method not yet implemented");
endswitch
## replace missings with appropriate fill values
A(endgap_locs) = fill_vals;
## some methods remove fill locations, only process if not empty
if (any (endgap_locs(:)))
## replace missings with appropriate fill values
A(endgap_locs) = fill_vals;
endif
endif
endswitch
endif
endif
if (nargout > 1)
idx_out = (missing_locs | endgap_locs);
endif
if (reshape_flag)
## revert permutation
A = permute (A, dim_idx_perm);
if (nargout > 1)
idx_out = permute (idx_out, dim_idx_perm);
endif
endif
endfunction
## numeric tests
%!assert (fillmissing ([1 2 3], "constant", 99), [1 2 3])
%!assert (fillmissing ([1 2 NaN], "constant", 99), [1 2 99])
%!assert (fillmissing ([NaN 2 NaN], "constant", 99), [99 2 99])
%!assert (fillmissing ([1 2 3]', "constant", 99), [1 2 3]')
%!assert (fillmissing ([1 2 NaN]', "constant", 99), [1 2 99]')
%!assert (fillmissing ([1 2 3; 4 5 6], "constant", 99), [1 2 3; 4 5 6])
%!assert (fillmissing ([1 2 NaN; 4 NaN 6], "constant", 99), [1 2 99; 4 99 6])
%!assert (fillmissing ([NaN 2 NaN; 4 NaN 6], "constant", [97, 98, 99]), [97 2 99; 4 98 6])
%!test
%! x = cat (3, [1, 2, NaN; 4, NaN, 6], [NaN, 2, 3; 4, 5, NaN]);
%! y = cat (3, [1, 2, 99; 4, 99, 6], [99, 2, 3; 4, 5, 99]);
%! assert (fillmissing (x, "constant", 99), y);
%! y = cat (3, [1, 2, 96; 4, 95, 6], [97, 2, 3; 4, 5, 99]);
%! assert (fillmissing (x, "constant", [94:99]), y);
%! assert (fillmissing (x, "constant", [94:99]'), y);
%! assert (fillmissing (x, "constant", permute ([94:99], [1 3 2])), y);
%! assert (fillmissing (x, "constant", [94, 96, 98; 95, 97, 99]), y);
%! assert (fillmissing (x, "constant", [94:99], 1), y);
%! y = cat (3, [1, 2, 96; 4, 97, 6], [98, 2, 3; 4, 5, 99]);
%! assert (fillmissing (x, "constant", [96:99], 2), y);
%! y = cat (3, [1, 2, 98; 4, 97, 6], [94, 2, 3; 4, 5, 99]);
%! assert (fillmissing (x, "constant", [94:99], 3), y);
%! y = cat (3, [1, 2, 92; 4, 91, 6], [94, 2, 3; 4, 5, 99]);
%! assert (fillmissing (x, "constant", [88:99], 99), y);
%!test
%! x = reshape([1:24],4,3,2);
%! x([1, 6, 7, 9, 12, 14, 16, 19, 22, 23]) = NaN;
%! y = x;
%! y([1,6,7,9 12, 14, 16, 19, 22, 23]) = [94 95 95 96 96 97 97 98 99 99];
%! assert (fillmissing (x, "constant", [94:99], 1), y);
%! y([1,6,7,9 12, 14, 16, 19, 22, 23]) = [92 93 94 92 95 97 99 98 97 98];
%! assert (fillmissing (x, "constant", [92:99], 2), y);
%! y([1,6,7,9 12, 14, 16, 19, 22, 23]) = [88 93 94 96 99 89 91 94 97 98];
%! assert (fillmissing (x, "constant", [88:99], 3), y);
%! y([1,6,7,9 12, 14, 16, 19, 22, 23]) = [76 81 82 84 87 89 91 94 97 98];
%! assert (fillmissing (x, "constant", [76:99], 99), y);
## tests with different endvalues behavior
%!assert (fillmissing ([1 2 3], "constant", 99, "endvalues", 88), [1 2 3])
%!assert (fillmissing ([1 NaN 3], "constant", 99, "endvalues", 88), [1 99 3])
%!assert (fillmissing ([1 2 NaN], "constant", 99, "endvalues", 88), [1 2 88])
%!assert (fillmissing ([NaN 2 3], "constant", 99, "endvalues", 88), [88 2 3])
%!assert (fillmissing ([NaN NaN 3], "constant", 99, "endvalues", 88), [88 88 3])
%!assert (fillmissing ([1 NaN NaN], "constant", 99, "endvalues", 88), [1 88 88])
%!assert (fillmissing ([NaN 2 NaN], "constant", 99, "endvalues", 88), [88 2 88])
%!assert (fillmissing ([NaN 2 NaN]', "constant", 99, "endvalues", 88), [88 2 88]')
%!assert (fillmissing ([1 NaN 3 NaN 5], "constant", 99, "endvalues", 88), [1 99 3 99 5])
%!assert (fillmissing ([1 NaN NaN NaN 5], "constant", 99, "endvalues", 88), [1 99 99 99 5])
%!assert (fillmissing ([NaN NaN NaN NaN 5], "constant", 99, "endvalues", 88), [88 88 88 88 5])
%!assert (fillmissing ([1 NaN 3 4 NaN], "constant", 99, "endvalues", 88), [1 99 3 4 88])
%!assert (fillmissing ([1 NaN 3 4 NaN], "constant", 99, 1, "endvalues", 88), [1 88 3 4 88])
%!assert (fillmissing ([1 NaN 3 4 NaN], "constant", 99, 1, "endvalues", "extrap"), [1 99 3 4 99])
%!test
%! x = reshape ([1:24], 3, 4, 2);
%! y = x;
%! x([1,2,5,6,8,10,13,16,18,19,20,21,22]) = NaN;
%! y([1,2,5,6,10,13,16,18,19,20,21,22])= 88; y([8])=99;
%! assert (fillmissing (x, "constant", 99, "endvalues", 88), y);
%! assert (fillmissing (x, "constant", 99, 1, "endvalues", 88), y);
%! y = x; y([1,2,5,8,10,13,16,19,22])= 88; y([6,18,20,21])=99;
%! assert (fillmissing (x, "constant", 99, 2, "endvalues", 88), y);
%! y(y==99) = 88;
%! assert (fillmissing (x, "constant", 99, 3, "endvalues", 88), y);
%! assert (fillmissing (x, "constant", 99, 4, "endvalues", 88), y);
%! assert (fillmissing (x, "constant", 99, 99, "endvalues", 88), y);
%! y([8]) = 94;
%! assert (fillmissing (x, "constant", [92:99], 1, "endvalues", 88), y);
%! y([6,8,18,20,21]) = [96,88,99,98,99];
%! assert (fillmissing (x, "constant", [94:99], 2, "endvalues", 88), y);
%! y = x; y(isnan(y)) = 88;
%! assert (fillmissing (x, "constant", [88:99], 3, "endvalues", 88), y);
%! y = x; y(isnan(y)) = [82,82,83,83,94,85,86,87,87,88,88,88,89];
%! assert (fillmissing (x, "constant", [92:99], 1, "endvalues", [82:89]), y);
%! y = x; y(isnan(y)) = [84,85,85,96,85,84,87,87,99,87,98,99,87];
%! assert (fillmissing (x, "constant", [94:99], 2, "endvalues", [84:89]), y);
%! y = x; y(isnan(y)) = [68,69,72,73,75,77,68,71,73,74,75,76,77];
%! assert (fillmissing (x, "constant", [88:99], 3, "endvalues", [68:79]), y);
%! assert (fillmissing (x, "constant", [88:93;94:99]', 3, "endvalues", [68:73;74:79]'), y)
%!test
%! x = reshape([1:24],4,3,2);
%! x([1, 6, 7, 9, 12, 14, 16, 19, 22, 23]) = NaN;
%! y = x;
%! y([1,6,7,9 12, 14, 16, 19, 22, 23]) = [94 95 95 96 96 97 97 98 99 99];
%! assert (fillmissing (x, "constant", [94:99], 1), y);
%! y([1,6,7,9 12, 14, 16, 19, 22, 23]) = [92 93 94 92 95 97 99 98 97 98];
%! assert (fillmissing (x, "constant", [92:99], 2), y);
%! y([1,6,7,9 12, 14, 16, 19, 22, 23]) = [88 93 94 96 99 89 91 94 97 98];
%! assert (fillmissing (x, "constant", [88:99], 3), y);
%! y([1,6,7,9 12, 14, 16, 19, 22, 23]) = [76 81 82 84 87 89 91 94 97 98];
%! assert (fillmissing (x, "constant", [76:99], 99), y);
## next/previous tests
%!assert (fillmissing ([1 2 3], "previous"), [1 2 3])
%!assert (fillmissing ([1 2 3], "next"), [1 2 3])
%!assert (fillmissing ([1 2 3]', "previous"), [1 2 3]')
%!assert (fillmissing ([1 2 3]', "next"), [1 2 3]')
%!assert (fillmissing ([1 2 NaN], "previous"), [1 2 2])
%!assert (fillmissing ([1 2 NaN], "next"), [1 2 NaN])
%!assert (fillmissing ([NaN 2 NaN], "previous"), [NaN 2 2])
%!assert (fillmissing ([NaN 2 NaN], "next"), [2 2 NaN])
%!assert (fillmissing ([1 NaN 3], "previous"), [1 1 3])
%!assert (fillmissing ([1 NaN 3], "next"), [1 3 3])
%!assert (fillmissing ([1 2 NaN; 4 NaN 6], "previous", 1), [1 2 NaN; 4 2 6])
%!assert (fillmissing ([1 2 NaN; 4 NaN 6], "previous", 2), [1 2 2; 4 4 6])
%!assert (fillmissing ([1 2 NaN; 4 NaN 6], "previous", 3), [1 2 NaN; 4 NaN 6])
%!assert (fillmissing ([1 2 NaN; 4 NaN 6], "next", 1), [1 2 6; 4 NaN 6])
%!assert (fillmissing ([1 2 NaN; 4 NaN 6], "next", 2), [1 2 NaN; 4 6 6])
%!assert (fillmissing ([1 2 NaN; 4 NaN 6], "next", 3), [1 2 NaN; 4 NaN 6])
%!test
%! x = reshape([1:24],4,3,2);
%! x([1, 6, 7, 9, 12, 14, 16, 19, 22, 23]) = NaN;
%! y = x;
%! y([1, 6, 7, 9, 14, 19, 22, 23]) = [2 8 8 10 15 20 24 24];
%! assert (fillmissing (x, "next", 1), y);
%! y = x;
%! y([1, 6, 7, 14, 16]) = [5, 10, 11, 18, 20];
%! assert (fillmissing (x, "next", 2), y);
%! y = x;
%! y([1, 6, 9, 12]) = [13 18 21 24];
%! assert (fillmissing (x, "next", 3), y);
%! assert (fillmissing (x, "next", 99), x);
%! y = x;
%! y([6, 7, 12, 14, 16, 19, 22, 23]) = [5 5 11 13 15 18 21 21];
%! assert (fillmissing (x, "previous", 1), y);
%! y = x;
%! y([6, 7, 9, 12, 19, 22, 23]) = [2 3 5 8 15 18 15];
%! assert (fillmissing (x, "previous", 2), y);
%! y = x;
%! y([14, 16, 22, 23]) = [2 4 10 11];
%! assert (fillmissing (x, "previous", 3), y);
%! assert (fillmissing (x, "previous", 99), x);
## next/previous tests with different endvalue behavior
%!assert (fillmissing ([1 2 3], "constant", 0, "endvalues", "previous"), [1 2 3])
%!assert (fillmissing ([1 2 3], "constant", 0, "endvalues", "next"), [1 2 3])
%!assert (fillmissing ([1 NaN 3], "constant", 0, "endvalues", "previous"), [1 0 3])
%!assert (fillmissing ([1 NaN 3], "constant", 0, "endvalues", "next"), [1 0 3])
%!assert (fillmissing ([1 2 NaN], "constant", 0, "endvalues", "previous"), [1 2 2])
%!assert (fillmissing ([1 2 NaN], "constant", 0, "endvalues", "next"), [1 2 NaN])
%!assert (fillmissing ([1 NaN NaN], "constant", 0, "endvalues", "previous"), [1 1 1])
%!assert (fillmissing ([1 NaN NaN], "constant", 0, "endvalues", "next"), [1 NaN NaN])
%!assert (fillmissing ([NaN 2 3], "constant", 0, "endvalues", "previous"), [NaN 2 3])
%!assert (fillmissing ([NaN 2 3], "constant", 0, "endvalues", "next"), [2 2 3])
%!assert (fillmissing ([NaN NaN 3], "constant", 0, "endvalues", "previous"), [NaN NaN 3])
%!assert (fillmissing ([NaN NaN 3], "constant", 0, "endvalues", "next"), [3 3 3])
%!assert (fillmissing ([NaN NaN NaN], "constant", 0, "endvalues", "previous"), [NaN NaN NaN])
%!assert (fillmissing ([NaN NaN NaN], "constant", 0, "endvalues", "next"), [NaN NaN NaN])
%!assert (fillmissing ([NaN 2 NaN 4 NaN], "constant", 0, "endvalues", "previous"), [NaN 2 0 4 4])
%!assert (fillmissing ([NaN 2 NaN 4 NaN], "constant", 0, "endvalues", "next"), [2 2 0 4 NaN])
%!assert (fillmissing ([NaN 2 NaN 4 NaN], "constant", 0, 1, "endvalues", "previous"), [NaN 2 NaN 4 NaN])
%!assert (fillmissing ([NaN 2 NaN 4 NaN], "constant", 0, 1, "endvalues", "next"), [NaN 2 NaN 4 NaN])
%!assert (fillmissing ([NaN 2 NaN 4 NaN], "constant", 0, 2, "endvalues", "previous"), [NaN 2 0 4 4])
%!assert (fillmissing ([NaN 2 NaN 4 NaN], "constant", 0, 2, "endvalues", "next"), [2 2 0 4 NaN])
%!assert (fillmissing ([NaN 2 NaN 4 NaN], "constant", 0, 3, "endvalues", "previous"), [NaN 2 NaN 4 NaN])
%!assert (fillmissing ([NaN 2 NaN 4 NaN], "constant", 0, 3, "endvalues", "next"), [NaN 2 NaN 4 NaN])
%!test
%! x = reshape ([1:24], 3, 4, 2);
%! x([1,2,5,6,8,10,13,16,18,19,20,21,22]) = NaN;
%! y = x;
%! y([5,6,8,18])=[4,4,0,17];
%! assert (fillmissing (x, "constant", 0, "endvalues", "previous"), y);
%! assert (fillmissing (x, "constant", 0, 1, "endvalues", "previous"), y);
%! y = x;
%! y([6,10,18,20,21])=[0,7,0,0,0];
%! assert (fillmissing (x, "constant", 0, 2, "endvalues", "previous"), y);
%! y = x;
%! y([16,19,21])=[4,7,9];
%! assert (fillmissing (x, "constant", 0, 3, "endvalues", "previous"), y);
%! assert (fillmissing (x, "constant", 0, 4, "endvalues", "previous"), x);
%! assert (fillmissing (x, "constant", 0, 99, "endvalues", "previous"), x);
%! y = x;
%! y([1,2,8,10,13,16,22])=[3,3,0,11,14,17,23];
%! assert (fillmissing (x, "constant", 0, "endvalues", "next"), y);
%! assert (fillmissing (x, "constant", 0, 1, "endvalues", "next"), y);
%! y = x;
%! y([1,2,5,6,8,18,20,21])=[4,11,11,0,11,0,0,0];
%! assert (fillmissing (x, "constant", 0, 2, "endvalues", "next"), y);
%! y = x;
%! y([2,5])=[14,17];
%! assert (fillmissing (x, "constant", 0, 3, "endvalues", "next"), y);
%! assert (fillmissing (x, "constant", 0, 4, "endvalues", "next"), x);
%! assert (fillmissing (x, "constant", 0, 99, "endvalues", "next"), x);
##tests for nearest
%!assert (fillmissing ([1 2 3], "nearest"), [1 2 3])
%!assert (fillmissing ([1 2 3]', "nearest"), [1 2 3]')
%!assert (fillmissing ([1 2 NaN], "nearest"), [1 2 2])
%!assert (fillmissing ([NaN 2 NaN], "nearest"), [2 2 2])
%!assert (fillmissing ([1 NaN 3], "nearest"), [1 3 3])
%!assert (fillmissing ([1 2 NaN; 4 NaN 6], "nearest", 1), [1 2 6; 4 2 6])
%!assert (fillmissing ([1 2 NaN; 4 NaN 6], "nearest", 2), [1 2 2; 4 6 6])
%!assert (fillmissing ([1 2 NaN; 4 NaN 6], "nearest", 3), [1 2 NaN; 4 NaN 6])
%!assert (fillmissing ([1 NaN 3 NaN 5], "nearest"), [1 3 3 5 5])
%!assert (fillmissing ([1 NaN 3 NaN 5], "nearest", "samplepoints", [0 1 2 3 4]), [1 3 3 5 5])
%!assert (fillmissing ([1 NaN 3 NaN 5], "nearest", "samplepoints", [0.5 1 2 3 5]), [1 1 3 3 5])
%!test
%! x = reshape([1:24],4,3,2);
%! x([1, 6, 7, 9, 12, 14, 16, 19, 22, 23]) = NaN;
%! y = x;
%! y([1, 6, 7, 9, 12, 14, 16, 19, 22, 23]) = [2 5 8 10 11 15 15 20 21 24];
%! assert (fillmissing (x, "nearest", 1), y);
%! y = x;
%! y([1, 6, 7, 9, 12, 14, 16, 19, 22, 23]) = [5 10 11 5 8 18 20 15 18 15];
%! assert (fillmissing (x, "nearest", 2), y);
%! y = x;
%! y([1, 6, 9, 12, 14, 16, 22, 23]) = [13 18 21 24 2 4 10 11];
%! assert (fillmissing (x, "nearest", 3), y);
%! assert (fillmissing (x, "nearest", 99), x);
##tests for nearest with diff endvalue behavior
%!assert (fillmissing ([1 2 3], "constant", 0, "endvalues", "nearest"), [1 2 3])
%!assert (fillmissing ([1 NaN 3], "constant", 0, "endvalues", "nearest"), [1 0 3])
%!assert (fillmissing ([1 2 NaN], "constant", 0, "endvalues", "nearest"), [1 2 2])
%!assert (fillmissing ([1 NaN NaN], "constant", 0, "endvalues", "nearest"), [1 1 1])
%!assert (fillmissing ([NaN 2 3], "constant", 0, "endvalues", "nearest"), [2 2 3])
%!assert (fillmissing ([NaN NaN 3], "constant", 0, "endvalues", "nearest"), [3 3 3])
%!assert (fillmissing ([NaN NaN NaN], "constant", 0, "endvalues", "nearest"), [NaN NaN NaN])
%!assert (fillmissing ([NaN 2 NaN 4 NaN], "constant", 0, "endvalues", "nearest"), [2 2 0 4 4])
%!assert (fillmissing ([NaN 2 NaN 4 NaN], "constant", 0, 1, "endvalues", "nearest"), [NaN 2 NaN 4 NaN])
%!assert (fillmissing ([NaN 2 NaN 4 NaN], "constant", 0, 2, "endvalues", "nearest"), [2 2 0 4 4])
%!assert (fillmissing ([NaN 2 NaN 4 NaN], "constant", 0, 3, "endvalues", "nearest"), [NaN 2 NaN 4 NaN])
%!test
%! x = reshape ([1:24], 3, 4, 2);
%! x([1,2,5,6,8,10,13,16,18,19,20,21,22]) = NaN;
%! y = x;
%! y([1,2,5,6,8,10,13,16,18,22])=[3 3 4 4 0 11 14 17 17 23];
%! assert (fillmissing (x, "constant", 0, "endvalues", "nearest"), y);
%! assert (fillmissing (x, "constant", 0, 1, "endvalues", "nearest"), y);
%! y = x;
%! y([1,2,5,6,8,10,18,20,21])=[4 11 11 0 11 7 0 0 0];
%! assert (fillmissing (x, "constant", 0, 2, "endvalues", "nearest"), y);
%! y = x;
%! y([2,5,16,19,21])=[14 17 4 7 9];
%! assert (fillmissing (x, "constant", 0, 3, "endvalues", "nearest"), y);
%! assert (fillmissing (x, "constant", 0, 99, "endvalues", "nearest"), x);
##tests for linear
%!assert (fillmissing ([1 2 3], "linear"), [1 2 3])
%!assert (fillmissing ([1 2 3]', "linear"), [1 2 3]')
%!assert (fillmissing ([1 2 NaN], "linear"), [1 2 3])
%!assert (fillmissing ([NaN 2 NaN], "linear"), [NaN 2 NaN])
%!assert (fillmissing ([1 NaN 3], "linear"), [1 2 3])
%!assert (fillmissing ([1 2 NaN; 4 NaN 6], "linear", 1), [1 2 NaN; 4 NaN 6])
%!assert (fillmissing ([1 2 NaN; 4 NaN 6], "linear", 2), [1 2 3; 4 5 6])
%!assert (fillmissing ([1 2 NaN; 4 NaN 6], "linear", 3), [1 2 NaN; 4 NaN 6])
%!assert (fillmissing ([1 NaN 3 NaN 5], "linear"), [1 2 3 4 5])
%!assert (fillmissing ([1 NaN 3 NaN 5], "linear", "samplepoints", [0 1 2 3 4]), [1 2 3 4 5])
%!assert (fillmissing ([1 NaN 3 NaN 5], "linear", "samplepoints", [0 1.5 2 5 14]), [1 2.5 3 3.5 5], eps)
%!test
%! x = reshape([1:24],4,3,2);
%! x([1, 6, 7, 9, 12, 14, 16, 19, 22, 23]) = NaN;
%! assert (fillmissing (x, "linear", 1), reshape([1:24],4,3,2));
%! y = reshape([1:24],4,3,2);
%! y([1 9 14 19 22 23]) = NaN;
%! assert (fillmissing (x, "linear", 2), y);
%! y = reshape([1:24],4,3,2);
%! y([1, 6, 7, 9, 12, 14, 16, 19, 22, 23]) = NaN;
%! assert (fillmissing (x, "linear", 3), y);
%! assert (fillmissing (x, "linear", 99), x);
##tests for linear with diff endvalue behavior
%!assert (fillmissing ([1 2 3], "linear", "endvalues", 0), [1 2 3])
%!assert (fillmissing ([1 NaN 3], "linear", "endvalues", 0), [1 2 3])
%!assert (fillmissing ([1 2 NaN], "linear", "endvalues", 0), [1 2 0])
%!assert (fillmissing ([1 NaN NaN], "linear", "endvalues", 0), [1 0 0])
%!assert (fillmissing ([NaN 2 3], "linear", "endvalues", 0), [0 2 3])
%!assert (fillmissing ([NaN NaN 3], "linear", "endvalues", 0), [0 0 3])
%!assert (fillmissing ([NaN NaN NaN], "linear", "endvalues", 0), [0 0 0])
%!assert (fillmissing ([NaN 2 NaN 4 NaN], "linear", "endvalues", 0), [0 2 3 4 0])
%!assert (fillmissing ([NaN 2 NaN 4 NaN], "linear", 1, "endvalues", 0), [0 2 0 4 0])
%!assert (fillmissing ([NaN 2 NaN 4 NaN], "linear", 2, "endvalues", 0), [0 2 3 4 0])
%!assert (fillmissing ([NaN 2 NaN 4 NaN], "linear", 3, "endvalues", 0), [0 2 0 4 0])
%!test
%! x = reshape ([1:24], 3, 4, 2);
%! x([1,2,5,6,8,10,13,16,18,19,20,21,22]) = NaN;
%! y = x;
%! y([1,2,5,6,10,13,16,18,19,20,21,22])=0; y(8)=8;
%! assert (fillmissing (x, "linear", "endvalues", 0), y);
%! assert (fillmissing (x, "linear", 1, "endvalues", 0), y);
%! y = x;
%! y([1,2,5,8,10,13,16,19,22])=0; y([6,18,20,21])=[6,18,20,21];
%! assert (fillmissing (x, "linear", 2, "endvalues", 0), y);
%! y = x;
%! y(isnan(y))=0;
%! assert (fillmissing (x, "linear", 3, "endvalues", 0), y);
%! assert (fillmissing (x, "linear", 99, "endvalues", 0), y);
##test other interpolants
%! x = reshape ([1:24], 3, 4, 2);
%! x([1,2,5,6,8,10,13,16,18,19,20,21,22]) = NaN;
%! y = x;
%! y([1,6,10,18,20,21]) = [2.5, 5, 8.5, 17.25, 21, 21.75];
%! assert (fillmissing (x, "linear", 2, "samplepoints", [2 4 8 10]), y, eps);
%! y([1,6,10,18,20,21]) = [2.5, 4.5, 8.5, 17.25, 21.5, 21.75];
%! assert (fillmissing (x, "spline", 2, "samplepoints", [2 4 8 10]), y, eps);
%! y([1,6,10,18,20,21]) = [2.5, 4.559386973180077, 8.5, 17.25, 21.440613026819925, 21.75];
%! assert (fillmissing (x, "pchip", 2, "samplepoints", [2 4 8 10]), y, 10*eps);
## known fail: makima method not yet implemented in interp1
%!test <60965>
%! x = reshape ([1:24], 3, 4, 2);
%! x([1,2,5,6,8,10,13,16,18,19,20,21,22]) = NaN;
%! y = x;
%! y([1,6,10,18,20,21]) = [2.5, 4.609523809523809, 8.5, 17.25, 21.390476190476186, 21.75];
%! assert (fillmissing (x, "makima", 2, "samplepoints", [2 4 8 10]), y, 10*eps);
##test maxgap for mid and end points
%!assert (fillmissing ([1 2 3], "constant", 0, "maxgap", 1), [1 2 3])
%!assert (fillmissing ([1 2 3], "constant", 0, "maxgap", 99), [1 2 3])
%!assert (fillmissing ([1 NaN 3], "constant", 0, "maxgap", 1), [1 NaN 3])
%!assert (fillmissing ([1 NaN 3], "constant", 0, "maxgap", 1.999), [1 NaN 3])
%!assert (fillmissing ([1 NaN 3], "constant", 0, "maxgap", 2), [1 0 3])
%!assert (fillmissing ([1 NaN NaN 4], "constant", 0, "maxgap", 2), [1 NaN NaN 4])
%!assert (fillmissing ([1 NaN NaN 4], "constant", 0, "maxgap", 3), [1 0 0 4])
%!assert (fillmissing ([1 NaN 3 NaN 5], "constant", 0, "maxgap", 2), [1 0 3 0 5])
%!assert (fillmissing ([NaN 2 NaN], "constant", 0, "maxgap", 0.999), [NaN 2 NaN])
%!assert (fillmissing ([NaN 2 NaN], "constant", 0, "maxgap", 1), [0 2 0])
%!assert (fillmissing ([NaN 2 NaN NaN], "constant", 0, "maxgap", 1), [0 2 NaN NaN])
%!assert (fillmissing ([NaN 2 NaN NaN], "constant", 0, "maxgap", 2), [0 2 0 0])
%!assert (fillmissing ([NaN NaN NaN], "constant", 0, "maxgap", 1), [NaN NaN NaN])
%!assert (fillmissing ([NaN NaN NaN], "constant", 0, "maxgap", 3), [NaN NaN NaN])
%!assert (fillmissing ([NaN NaN NaN], "constant", 0, "maxgap", 999), [NaN NaN NaN])
%!assert (fillmissing ([1 NaN 3 NaN 5], "constant", 0, "maxgap", 2, "samplepoints", [0 1 2 3 5]), [1 0 3 NaN 5])
%!assert (fillmissing ([1 NaN 3 NaN 5]', "constant", 0, "maxgap", 2, "samplepoints", [0 1 2 3 5]), [1 0 3 NaN 5]')
%!assert (fillmissing ([1 NaN 3 NaN 5], "constant", 0, "maxgap", 2, "samplepoints", [0 2 3 4 5]), [1 NaN 3 0 5])
%!assert (fillmissing ([1 NaN 3 NaN 5; 1 NaN 3 NaN 5], "constant", 0, 2, "maxgap", 2, "samplepoints", [0 2 3 4 5]), [1 NaN 3 0 5; 1 NaN 3 0 5])
%!test
%! x = cat (3, [1, 2, NaN; 4, NaN, NaN], [NaN, 2, 3; 4, 5, NaN]);
%! assert (fillmissing (x, "constant", 0, "maxgap", 0.1), x);
%! y = x;
%! y([4,7,12]) = 0;
%! assert (fillmissing (x, "constant", 0, "maxgap", 1), y);
%! assert (fillmissing (x, "constant", 0, 1, "maxgap", 1), y);
%! y = x;
%! y([5,7,12]) = 0;
%! assert (fillmissing (x, "constant", 0, 2, "maxgap", 1), y);
%! y = x;
%! y([4,5,7]) = 0;
%! assert (fillmissing (x, "constant", 0, 3, "maxgap", 1), y);
##2nd output, verify consistent with dim
%!test
%! x = cat (3, [1, 2, NaN; 4, NaN, NaN], [NaN, 2, 3; 4, 5, NaN]);
%! [~, idx] = fillmissing (x, "constant", 0, "maxgap", 1);
%! assert (idx, logical(cat(3,[0 0 0; 0 1 0],[1 0 0; 0 0 1])));
%! [~, idx] = fillmissing (x, "constant", 0, 1, "maxgap", 1);
%! assert (idx, logical(cat(3,[0 0 0; 0 1 0],[1 0 0; 0 0 1])));
%! [~, idx] = fillmissing (x, "constant", 0, 2, "maxgap", 1);
%! assert (idx, logical(cat(3,[0 0 1; 0 0 0],[1 0 0; 0 0 1])));
%! [~, idx] = fillmissing (x, "constant", 0, 3, "maxgap", 1);
%! assert (idx, logical(cat(3,[0 0 1; 0 1 0],[1 0 0; 0 0 0])));
##TEST Char and cellstr
%!assert (fillmissing (' foo bar ', "constant", 'X'), 'XfooXbarX')
%!assert (fillmissing ([' foo';'bar '], "constant", 'X'), ['Xfoo';'barX'])
%!assert (fillmissing ([' foo';'bar '], "next"), ['bfoo';'bar '])
%!assert (fillmissing ([' foo';'bar '], "next", 1), ['bfoo';'bar '])
%!assert (fillmissing ([' foo';'bar '], "previous"), [' foo';'baro'])
%!assert (fillmissing ([' foo';'bar '], "previous", 1), [' foo';'baro'])
%!assert (fillmissing ([' foo';'bar '], "next", 2), ['ffoo';'bar '])
%!assert (fillmissing ([' foo';'bar '], "previous", 2), [' foo';'barr'])
%!assert (fillmissing ([' foo';'bar '], "next", 3), [' foo';'bar '])
%!assert (fillmissing ([' foo';'bar '], "previous", 3), [' foo';'bar '])
%!assert (fillmissing ({'foo','bar'}, "constant", 'a'), {'foo','bar'})
%!assert (fillmissing ({'foo','bar'}, "constant", {'a'}), {'foo','bar'})
%!assert (fillmissing ({'foo', '', 'bar'}, "constant", 'a'), {'foo', 'a', 'bar'})
%!assert (fillmissing ({'foo', '', 'bar'}, "constant", {'a'}), {'foo', 'a', 'bar'})
%!assert (fillmissing ({'foo', '', 'bar'}, "previous"), {'foo', 'foo', 'bar'})
%!assert (fillmissing ({'foo', '', 'bar'}, "next"), {'foo', 'bar', 'bar'})
%!assert (fillmissing ({'foo', '', 'bar'}, "previous", 2), {'foo', 'foo', 'bar'})
%!assert (fillmissing ({'foo', '', 'bar'}, "next", 2), {'foo', 'bar', 'bar'})
%!assert (fillmissing ({'foo', '', 'bar'}, "previous", 1), {'foo', '', 'bar'})
%!assert (fillmissing ({'foo', '', 'bar'}, "next", 1), {'foo', '', 'bar'})
## Test input validation and error messages
%!error fillmissing ()
%!error fillmissing (1)
%!error fillmissing (1,2,3,4,5,6,7,8,9,10,11,12,13)
%!error fillmissing (1, 2)
%!error fillmissing (1, "foo")
%!error fillmissing (1, @(x) x, 1)
%!error fillmissing (1, @(x,y) x+y, 1)
%!error <'constant' method must be followed by> fillmissing (1, "constant")
%!error fillmissing (1, "constant", [])
%!error fillmissing (1, "constant", "a")
%!error fillmissing ("a", "constant", 1)
%!error fillmissing ("a", "constant", {"foo"})
%!error fillmissing ({"foo"}, "constant", 1)
%!error fillmissing (1, "movemean")
%!error fillmissing (1, "movemedian")
%!error fillmissing (1, "constant", 1, 0)
%!error fillmissing (1, "constant", 1, -1)
%!error fillmissing (1, "constant", 1, [1 2])
%!error fillmissing (1, "constant", 1, "samplepoints")
%!error fillmissing (1, "constant", 1, "foo")
%!error fillmissing (1, "constant", 1, 1, "foo")
%!error fillmissing (1, "constant", 1, 2, {1}, 4)
%!error fillmissing ([1 2 3], "constant", 1, 2, "samplepoints", [1 2])
%!error fillmissing ([1 2 3], "constant", 1, 2, "samplepoints", [3 1 2])
%!error fillmissing ([1 2 3], "constant", 1, 2, "samplepoints", [1 1 2])
%!error fillmissing ([1 2 3], "constant", 1, 2, "samplepoints", "abc")
%!error fillmissing ([1 2 3], "constant", 1, 2, "samplepoints", logical([1 1 1]))
%!error fillmissing ([1 2 3], "constant", 1, 1, "samplepoints", [1 2 3])
%!error fillmissing ('foo', "next", "endvalues", 1)
%!error fillmissing (1, "constant", 1, 1, "endvalues", "foo")
%!error fillmissing (1, "constant", 1, 1, "endvalues", "linear")
%!error fillmissing ([1 2 3], "constant", 1, 2, "endvalues", [1 2 3])
%!error fillmissing ([1 2 3], "constant", 1, 1, "endvalues", [1 2])
%!error fillmissing (randi(5,4,3,2), "constant", 1, 3, "endvalues", [1 2])
%!error fillmissing (1, "constant", 1, 1, "endvalues", {1})
%!error fillmissing (1, "constant", 1, 2, "foo", 4)
%!error fillmissing (1, "constant", 1, 2, "maxgap", 1, "missinglocations", false)
%!error fillmissing (1, "constant", 1, 2, "missinglocations", false, "maxgap", 1)
%!error fillmissing (1, "constant", 1, 2, "missinglocations", 1)
%!error fillmissing (1, "constant", 1, 2, "missinglocations", 'a')
%!error fillmissing (1, "constant", 1, 2, "missinglocations", [true false])
%!error fillmissing (1, "constant", 1, 2, "maxgap", true)
%!error fillmissing (1, "constant", 1, 2, "maxgap", 'a')
%!error fillmissing (1, "constant", 1, 2, "maxgap", [1 2])
%!error fillmissing (1, "constant", 1, 2, "maxgap", 0)
%!error fillmissing (1, "constant", 1, 2, "maxgap", -1)
%!error fillmissing ([1 2 3], "constant", [1 2 3])
%!error fillmissing ([1 2 3]', "constant", [1 2 3])
%!error fillmissing ([1 2 3]', "constant", [1 2 3], 1)
%!error fillmissing ([1 2 3], "constant", [1 2 3], 2)
%!error fillmissing (randi(5,4,3,2), "constant", [1 2], 1)
%!error fillmissing (randi(5,4,3,2), "constant", [1 2], 2)
%!error fillmissing (randi(5,4,3,2), "constant", [1 2], 3)
%!error fillmissing (1, @(x,y,z) x+y+z)