Skip to content
Permalink
master
Go to file
 
 
Cannot retrieve contributors at this time
4478 lines (3561 sloc) 136 KB
"""
Implementation of math operations on Array objects.
"""
import math
from collections import namedtuple
from enum import IntEnum
from functools import partial
import operator
import numpy as np
import llvmlite.llvmpy.core as lc
from numba import generated_jit
from numba.core import types, cgutils
from numba.core.extending import overload, overload_method, register_jitable
from numba.np.numpy_support import as_dtype, type_can_asarray
from numba.np.numpy_support import numpy_version
from numba.np.numpy_support import is_nonelike
from numba.core.imputils import (lower_builtin, impl_ret_borrowed,
impl_ret_new_ref, impl_ret_untracked)
from numba.core.typing import signature
from numba.np.arrayobj import make_array, load_item, store_item, _empty_nd_impl
from numba.np.linalg import ensure_blas
from numba.core.extending import intrinsic
from numba.core.errors import RequireLiteralValue, TypingError
def _check_blas():
# Checks if a BLAS is available so e.g. dot will work
try:
ensure_blas()
except ImportError:
return False
return True
_HAVE_BLAS = _check_blas()
@intrinsic
def _create_tuple_result_shape(tyctx, shape_list, shape_tuple):
"""
This routine converts shape list where the axis dimension has already
been popped to a tuple for indexing of the same size. The original shape
tuple is also required because it contains a length field at compile time
whereas the shape list does not.
"""
# The new tuple's size is one less than the original tuple since axis
# dimension removed.
nd = len(shape_tuple) - 1
# The return type of this intrinsic is an int tuple of length nd.
tupty = types.UniTuple(types.intp, nd)
# The function signature for this intrinsic.
function_sig = tupty(shape_list, shape_tuple)
def codegen(cgctx, builder, signature, args):
lltupty = cgctx.get_value_type(tupty)
# Create an empty int tuple.
tup = cgutils.get_null_value(lltupty)
# Get the shape list from the args and we don't need shape tuple.
[in_shape, _] = args
def array_indexer(a, i):
return a[i]
# loop to fill the tuple
for i in range(nd):
dataidx = cgctx.get_constant(types.intp, i)
# compile and call array_indexer
data = cgctx.compile_internal(builder, array_indexer,
types.intp(shape_list, types.intp),
[in_shape, dataidx])
tup = builder.insert_value(tup, data, i)
return tup
return function_sig, codegen
@intrinsic
def _gen_index_tuple(tyctx, shape_tuple, value, axis):
"""
Generates a tuple that can be used to index a specific slice from an
array for sum with axis. shape_tuple is the size of the dimensions of
the input array. 'value' is the value to put in the indexing tuple
in the axis dimension and 'axis' is that dimension. For this to work,
axis has to be a const.
"""
if not isinstance(axis, types.Literal):
raise RequireLiteralValue('axis argument must be a constant')
# Get the value of the axis constant.
axis_value = axis.literal_value
# The length of the indexing tuple to be output.
nd = len(shape_tuple)
# If the axis value is impossible for the given size array then
# just fake it like it was for axis 0. This will stop compile errors
# when it looks like it could be called from array_sum_axis but really
# can't because that routine checks the axis mismatch and raise an
# exception.
if axis_value >= nd:
axis_value = 0
# Calculate the type of the indexing tuple. All the non-axis
# dimensions have slice2 type and the axis dimension has int type.
before = axis_value
after = nd - before - 1
types_list = []
types_list += [types.slice2_type] * before
types_list += [types.intp]
types_list += [types.slice2_type] * after
# Creates the output type of the function.
tupty = types.Tuple(types_list)
# Defines the signature of the intrinsic.
function_sig = tupty(shape_tuple, value, axis)
def codegen(cgctx, builder, signature, args):
lltupty = cgctx.get_value_type(tupty)
# Create an empty indexing tuple.
tup = cgutils.get_null_value(lltupty)
# We only need value of the axis dimension here.
# The rest are constants defined above.
[_, value_arg, _] = args
def create_full_slice():
return slice(None, None)
# loop to fill the tuple with slice(None,None) before
# the axis dimension.
# compile and call create_full_slice
slice_data = cgctx.compile_internal(builder, create_full_slice,
types.slice2_type(),
[])
for i in range(0, axis_value):
tup = builder.insert_value(tup, slice_data, i)
# Add the axis dimension 'value'.
tup = builder.insert_value(tup, value_arg, axis_value)
# loop to fill the tuple with slice(None,None) after
# the axis dimension.
for i in range(axis_value + 1, nd):
tup = builder.insert_value(tup, slice_data, i)
return tup
return function_sig, codegen
#----------------------------------------------------------------------------
# Basic stats and aggregates
@lower_builtin(np.sum, types.Array)
@lower_builtin("array.sum", types.Array)
def array_sum(context, builder, sig, args):
zero = sig.return_type(0)
def array_sum_impl(arr):
c = zero
for v in np.nditer(arr):
c += v.item()
return c
res = context.compile_internal(builder, array_sum_impl, sig, args,
locals=dict(c=sig.return_type))
return impl_ret_borrowed(context, builder, sig.return_type, res)
@register_jitable
def _array_sum_axis_nop(arr, v):
return arr
def gen_sum_axis_impl(is_axis_const, const_axis_val, op, zero):
def inner(arr, axis):
"""
function that performs sums over one specific axis
The third parameter to gen_index_tuple that generates the indexing
tuples has to be a const so we can't just pass "axis" through since
that isn't const. We can check for specific values and have
different instances that do take consts. Supporting axis summation
only up to the fourth dimension for now.
typing/arraydecl.py:sum_expand defines the return type for sum with
axis. It is one dimension less than the input array.
"""
ndim = arr.ndim
if not is_axis_const:
# Catch where axis is negative or greater than 3.
if axis < 0 or axis > 3:
raise ValueError("Numba does not support sum with axis "
"parameter outside the range 0 to 3.")
# Catch the case where the user misspecifies the axis to be
# more than the number of the array's dimensions.
if axis >= ndim:
raise ValueError("axis is out of bounds for array")
# Convert the shape of the input array to a list.
ashape = list(arr.shape)
# Get the length of the axis dimension.
axis_len = ashape[axis]
# Remove the axis dimension from the list of dimensional lengths.
ashape.pop(axis)
# Convert this shape list back to a tuple using above intrinsic.
ashape_without_axis = _create_tuple_result_shape(ashape, arr.shape)
# Tuple needed here to create output array with correct size.
result = np.full(ashape_without_axis, zero, type(zero))
# Iterate through the axis dimension.
for axis_index in range(axis_len):
if is_axis_const:
# constant specialized version works for any valid axis value
index_tuple_generic = _gen_index_tuple(arr.shape, axis_index,
const_axis_val)
result += arr[index_tuple_generic]
else:
# Generate a tuple used to index the input array.
# The tuple is ":" in all dimensions except the axis
# dimension where it is "axis_index".
if axis == 0:
index_tuple1 = _gen_index_tuple(arr.shape, axis_index, 0)
result += arr[index_tuple1]
elif axis == 1:
index_tuple2 = _gen_index_tuple(arr.shape, axis_index, 1)
result += arr[index_tuple2]
elif axis == 2:
index_tuple3 = _gen_index_tuple(arr.shape, axis_index, 2)
result += arr[index_tuple3]
elif axis == 3:
index_tuple4 = _gen_index_tuple(arr.shape, axis_index, 3)
result += arr[index_tuple4]
return op(result, 0)
return inner
@lower_builtin(np.sum, types.Array, types.intp, types.DTypeSpec)
@lower_builtin(np.sum, types.Array, types.IntegerLiteral, types.DTypeSpec)
@lower_builtin("array.sum", types.Array, types.intp, types.DTypeSpec)
@lower_builtin("array.sum", types.Array, types.IntegerLiteral, types.DTypeSpec)
def array_sum_axis_dtype(context, builder, sig, args):
retty = sig.return_type
zero = getattr(retty, 'dtype', retty)(0)
# if the return is scalar in type then "take" the 0th element of the
# 0d array accumulator as the return value
if getattr(retty, 'ndim', None) is None:
op = np.take
else:
op = _array_sum_axis_nop
[ty_array, ty_axis, ty_dtype] = sig.args
is_axis_const = False
const_axis_val = 0
if isinstance(ty_axis, types.Literal):
# this special-cases for constant axis
const_axis_val = ty_axis.literal_value
# fix negative axis
if const_axis_val < 0:
const_axis_val = ty_array.ndim + const_axis_val
if const_axis_val < 0 or const_axis_val > ty_array.ndim:
raise ValueError("'axis' entry is out of bounds")
ty_axis = context.typing_context.resolve_value_type(const_axis_val)
axis_val = context.get_constant(ty_axis, const_axis_val)
# rewrite arguments
args = args[0], axis_val, args[2]
# rewrite sig
sig = sig.replace(args=[ty_array, ty_axis, ty_dtype])
is_axis_const = True
gen_impl = gen_sum_axis_impl(is_axis_const, const_axis_val, op, zero)
compiled = register_jitable(gen_impl)
def array_sum_impl_axis(arr, axis, dtype):
return compiled(arr, axis)
res = context.compile_internal(builder, array_sum_impl_axis, sig, args)
return impl_ret_new_ref(context, builder, sig.return_type, res)
@lower_builtin(np.sum, types.Array, types.DTypeSpec)
@lower_builtin("array.sum", types.Array, types.DTypeSpec)
def array_sum_dtype(context, builder, sig, args):
zero = sig.return_type(0)
def array_sum_impl(arr, dtype):
c = zero
for v in np.nditer(arr):
c += v.item()
return c
res = context.compile_internal(builder, array_sum_impl, sig, args,
locals=dict(c=sig.return_type))
return impl_ret_borrowed(context, builder, sig.return_type, res)
@lower_builtin(np.sum, types.Array, types.intp)
@lower_builtin(np.sum, types.Array, types.IntegerLiteral)
@lower_builtin("array.sum", types.Array, types.intp)
@lower_builtin("array.sum", types.Array, types.IntegerLiteral)
def array_sum_axis(context, builder, sig, args):
retty = sig.return_type
zero = getattr(retty, 'dtype', retty)(0)
# if the return is scalar in type then "take" the 0th element of the
# 0d array accumulator as the return value
if getattr(retty, 'ndim', None) is None:
op = np.take
else:
op = _array_sum_axis_nop
[ty_array, ty_axis] = sig.args
is_axis_const = False
const_axis_val = 0
if isinstance(ty_axis, types.Literal):
# this special-cases for constant axis
const_axis_val = ty_axis.literal_value
# fix negative axis
if const_axis_val < 0:
const_axis_val = ty_array.ndim + const_axis_val
if const_axis_val < 0 or const_axis_val > ty_array.ndim:
raise ValueError("'axis' entry is out of bounds")
ty_axis = context.typing_context.resolve_value_type(const_axis_val)
axis_val = context.get_constant(ty_axis, const_axis_val)
# rewrite arguments
args = args[0], axis_val
# rewrite sig
sig = sig.replace(args=[ty_array, ty_axis])
is_axis_const = True
gen_impl = gen_sum_axis_impl(is_axis_const, const_axis_val, op, zero)
compiled = register_jitable(gen_impl)
def array_sum_impl_axis(arr, axis):
return compiled(arr, axis)
res = context.compile_internal(builder, array_sum_impl_axis, sig, args)
return impl_ret_new_ref(context, builder, sig.return_type, res)
@lower_builtin(np.prod, types.Array)
@lower_builtin("array.prod", types.Array)
def array_prod(context, builder, sig, args):
def array_prod_impl(arr):
c = 1
for v in np.nditer(arr):
c *= v.item()
return c
res = context.compile_internal(builder, array_prod_impl, sig, args,
locals=dict(c=sig.return_type))
return impl_ret_borrowed(context, builder, sig.return_type, res)
@lower_builtin(np.cumsum, types.Array)
@lower_builtin("array.cumsum", types.Array)
def array_cumsum(context, builder, sig, args):
scalar_dtype = sig.return_type.dtype
dtype = as_dtype(scalar_dtype)
zero = scalar_dtype(0)
def array_cumsum_impl(arr):
out = np.empty(arr.size, dtype)
c = zero
for idx, v in enumerate(arr.flat):
c += v
out[idx] = c
return out
res = context.compile_internal(builder, array_cumsum_impl, sig, args,
locals=dict(c=scalar_dtype))
return impl_ret_new_ref(context, builder, sig.return_type, res)
@lower_builtin(np.cumprod, types.Array)
@lower_builtin("array.cumprod", types.Array)
def array_cumprod(context, builder, sig, args):
scalar_dtype = sig.return_type.dtype
dtype = as_dtype(scalar_dtype)
def array_cumprod_impl(arr):
out = np.empty(arr.size, dtype)
c = 1
for idx, v in enumerate(arr.flat):
c *= v
out[idx] = c
return out
res = context.compile_internal(builder, array_cumprod_impl, sig, args,
locals=dict(c=scalar_dtype))
return impl_ret_new_ref(context, builder, sig.return_type, res)
@lower_builtin(np.mean, types.Array)
@lower_builtin("array.mean", types.Array)
def array_mean(context, builder, sig, args):
zero = sig.return_type(0)
def array_mean_impl(arr):
# Can't use the naive `arr.sum() / arr.size`, as it would return
# a wrong result on integer sum overflow.
c = zero
for v in np.nditer(arr):
c += v.item()
return c / arr.size
res = context.compile_internal(builder, array_mean_impl, sig, args,
locals=dict(c=sig.return_type))
return impl_ret_untracked(context, builder, sig.return_type, res)
@lower_builtin(np.var, types.Array)
@lower_builtin("array.var", types.Array)
def array_var(context, builder, sig, args):
def array_var_impl(arr):
# Compute the mean
m = arr.mean()
# Compute the sum of square diffs
ssd = 0
for v in np.nditer(arr):
val = (v.item() - m)
ssd += np.real(val * np.conj(val))
return ssd / arr.size
res = context.compile_internal(builder, array_var_impl, sig, args)
return impl_ret_untracked(context, builder, sig.return_type, res)
@lower_builtin(np.std, types.Array)
@lower_builtin("array.std", types.Array)
def array_std(context, builder, sig, args):
def array_std_impl(arry):
return arry.var() ** 0.5
res = context.compile_internal(builder, array_std_impl, sig, args)
return impl_ret_untracked(context, builder, sig.return_type, res)
def zero_dim_msg(fn_name):
msg = ("zero-size array to reduction operation "
"{0} which has no identity".format(fn_name))
return msg
def _is_nat(x):
pass
@overload(_is_nat)
def ol_is_nat(x):
if numpy_version >= (1, 18):
return lambda x: np.isnat(x)
else:
nat = x('NaT')
return lambda x: x == nat
@lower_builtin(np.min, types.Array)
@lower_builtin("array.min", types.Array)
def array_min(context, builder, sig, args):
ty = sig.args[0].dtype
MSG = zero_dim_msg('minimum')
if isinstance(ty, (types.NPDatetime, types.NPTimedelta)):
# NP < 1.18: NaT is smaller than every other value, but it is
# ignored as far as min() is concerned.
# NP >= 1.18: NaT dominates like NaN
def array_min_impl(arry):
if arry.size == 0:
raise ValueError(MSG)
it = np.nditer(arry)
min_value = next(it).take(0)
if _is_nat(min_value):
return min_value
for view in it:
v = view.item()
if _is_nat(v):
if numpy_version >= (1, 18):
return v
else:
continue
if v < min_value:
min_value = v
return min_value
elif isinstance(ty, types.Complex):
def array_min_impl(arry):
if arry.size == 0:
raise ValueError(MSG)
it = np.nditer(arry)
min_value = next(it).take(0)
for view in it:
v = view.item()
if v.real < min_value.real:
min_value = v
elif v.real == min_value.real:
if v.imag < min_value.imag:
min_value = v
return min_value
elif isinstance(ty, types.Float):
def array_min_impl(arry):
if arry.size == 0:
raise ValueError(MSG)
it = np.nditer(arry)
min_value = next(it).take(0)
if np.isnan(min_value):
return min_value
for view in it:
v = view.item()
if np.isnan(v):
return v
if v < min_value:
min_value = v
return min_value
else:
def array_min_impl(arry):
if arry.size == 0:
raise ValueError(MSG)
it = np.nditer(arry)
min_value = next(it).take(0)
for view in it:
v = view.item()
if v < min_value:
min_value = v
return min_value
res = context.compile_internal(builder, array_min_impl, sig, args)
return impl_ret_borrowed(context, builder, sig.return_type, res)
@lower_builtin(np.max, types.Array)
@lower_builtin("array.max", types.Array)
def array_max(context, builder, sig, args):
ty = sig.args[0].dtype
MSG = zero_dim_msg('maximum')
if isinstance(ty, (types.NPDatetime, types.NPTimedelta)):
# NP < 1.18: NaT is smaller than every other value, but it is
# ignored as far as min() is concerned.
# NP >= 1.18: NaT dominates like NaN
def array_max_impl(arry):
if arry.size == 0:
raise ValueError(MSG)
it = np.nditer(arry)
max_value = next(it).take(0)
if _is_nat(max_value):
return max_value
for view in it:
v = view.item()
if _is_nat(v):
if numpy_version >= (1, 18):
return v
else:
continue
if v > max_value:
max_value = v
return max_value
elif isinstance(ty, types.Complex):
def array_max_impl(arry):
if arry.size == 0:
raise ValueError(MSG)
it = np.nditer(arry)
max_value = next(it).take(0)
for view in it:
v = view.item()
if v.real > max_value.real:
max_value = v
elif v.real == max_value.real:
if v.imag > max_value.imag:
max_value = v
return max_value
elif isinstance(ty, types.Float):
def array_max_impl(arry):
if arry.size == 0:
raise ValueError(MSG)
it = np.nditer(arry)
max_value = next(it).take(0)
if np.isnan(max_value):
return max_value
for view in it:
v = view.item()
if np.isnan(v):
return v
if v > max_value:
max_value = v
return max_value
else:
def array_max_impl(arry):
if arry.size == 0:
raise ValueError(MSG)
it = np.nditer(arry)
max_value = next(it).take(0)
for view in it:
v = view.item()
if v > max_value:
max_value = v
return max_value
res = context.compile_internal(builder, array_max_impl, sig, args)
return impl_ret_borrowed(context, builder, sig.return_type, res)
@lower_builtin(np.argmin, types.Array)
@lower_builtin("array.argmin", types.Array)
def array_argmin(context, builder, sig, args):
ty = sig.args[0].dtype
if (isinstance(ty, (types.NPDatetime, types.NPTimedelta))):
def array_argmin_impl(arry):
if arry.size == 0:
raise ValueError("attempt to get argmin of an empty sequence")
it = np.nditer(arry)
min_value = next(it).take(0)
min_idx = 0
if _is_nat(min_value):
return min_idx
idx = 1
for view in it:
v = view.item()
if _is_nat(v):
if numpy_version >= (1, 18):
return idx
else:
idx += 1
continue
if v < min_value:
min_value = v
min_idx = idx
idx += 1
return min_idx
elif isinstance(ty, types.Float):
def array_argmin_impl(arry):
if arry.size == 0:
raise ValueError("attempt to get argmin of an empty sequence")
for v in arry.flat:
min_value = v
min_idx = 0
break
if np.isnan(min_value):
return min_idx
idx = 0
for v in arry.flat:
if np.isnan(v):
return idx
if v < min_value:
min_value = v
min_idx = idx
idx += 1
return min_idx
else:
def array_argmin_impl(arry):
if arry.size == 0:
raise ValueError("attempt to get argmin of an empty sequence")
for v in arry.flat:
min_value = v
min_idx = 0
break
else:
raise RuntimeError('unreachable')
idx = 0
for v in arry.flat:
if v < min_value:
min_value = v
min_idx = idx
idx += 1
return min_idx
res = context.compile_internal(builder, array_argmin_impl, sig, args)
return impl_ret_untracked(context, builder, sig.return_type, res)
@lower_builtin(np.argmax, types.Array)
@lower_builtin("array.argmax", types.Array)
def array_argmax(context, builder, sig, args):
ty = sig.args[0].dtype
if (isinstance(ty, (types.NPDatetime, types.NPTimedelta))):
def array_argmax_impl(arry):
if arry.size == 0:
raise ValueError("attempt to get argmax of an empty sequence")
it = np.nditer(arry)
max_value = next(it).take(0)
max_idx = 0
if _is_nat(max_value):
return max_idx
idx = 1
for view in it:
v = view.item()
if _is_nat(v):
if numpy_version >= (1, 18):
return idx
else:
idx += 1
continue
if v > max_value:
max_value = v
max_idx = idx
idx += 1
return max_idx
elif isinstance(ty, types.Float):
def array_argmax_impl(arry):
if arry.size == 0:
raise ValueError("attempt to get argmax of an empty sequence")
for v in arry.flat:
max_value = v
max_idx = 0
break
if np.isnan(max_value):
return max_idx
idx = 0
for v in arry.flat:
if np.isnan(v):
return idx
if v > max_value:
max_value = v
max_idx = idx
idx += 1
return max_idx
else:
def array_argmax_impl(arry):
if arry.size == 0:
raise ValueError("attempt to get argmax of an empty sequence")
for v in arry.flat:
max_value = v
max_idx = 0
break
idx = 0
for v in arry.flat:
if v > max_value:
max_value = v
max_idx = idx
idx += 1
return max_idx
res = context.compile_internal(builder, array_argmax_impl, sig, args)
return impl_ret_untracked(context, builder, sig.return_type, res)
@overload(np.all)
@overload_method(types.Array, "all")
def np_all(a):
def flat_all(a):
for v in np.nditer(a):
if not v.item():
return False
return True
return flat_all
@overload(np.any)
@overload_method(types.Array, "any")
def np_any(a):
def flat_any(a):
for v in np.nditer(a):
if v.item():
return True
return False
return flat_any
def get_isnan(dtype):
"""
A generic isnan() function
"""
if isinstance(dtype, (types.Float, types.Complex)):
return np.isnan
else:
@register_jitable
def _trivial_isnan(x):
return False
return _trivial_isnan
@register_jitable
def less_than(a, b):
return a < b
@register_jitable
def greater_than(a, b):
return a > b
@register_jitable
def check_array(a):
if a.size == 0:
raise ValueError('zero-size array to reduction operation not possible')
def _check_is_integer(v, name):
if not isinstance(v, (int, types.Integer)):
raise TypingError('{} must be an integer'.format(name))
def nan_min_max_factory(comparison_op, is_complex_dtype):
if is_complex_dtype:
def impl(a):
arr = np.asarray(a)
check_array(arr)
it = np.nditer(arr)
return_val = next(it).take(0)
for view in it:
v = view.item()
if np.isnan(return_val.real) and not np.isnan(v.real):
return_val = v
else:
if comparison_op(v.real, return_val.real):
return_val = v
elif v.real == return_val.real:
if comparison_op(v.imag, return_val.imag):
return_val = v
return return_val
else:
def impl(a):
arr = np.asarray(a)
check_array(arr)
it = np.nditer(arr)
return_val = next(it).take(0)
for view in it:
v = view.item()
if not np.isnan(v):
if not comparison_op(return_val, v):
return_val = v
return return_val
return impl
real_nanmin = register_jitable(
nan_min_max_factory(less_than, is_complex_dtype=False)
)
real_nanmax = register_jitable(
nan_min_max_factory(greater_than, is_complex_dtype=False)
)
complex_nanmin = register_jitable(
nan_min_max_factory(less_than, is_complex_dtype=True)
)
complex_nanmax = register_jitable(
nan_min_max_factory(greater_than, is_complex_dtype=True)
)
@overload(np.nanmin)
def np_nanmin(a):
dt = determine_dtype(a)
if np.issubdtype(dt, np.complexfloating):
return complex_nanmin
else:
return real_nanmin
@overload(np.nanmax)
def np_nanmax(a):
dt = determine_dtype(a)
if np.issubdtype(dt, np.complexfloating):
return complex_nanmax
else:
return real_nanmax
@overload(np.nanmean)
def np_nanmean(a):
if not isinstance(a, types.Array):
return
isnan = get_isnan(a.dtype)
def nanmean_impl(a):
c = 0.0
count = 0
for view in np.nditer(a):
v = view.item()
if not isnan(v):
c += v.item()
count += 1
# np.divide() doesn't raise ZeroDivisionError
return np.divide(c, count)
return nanmean_impl
@overload(np.nanvar)
def np_nanvar(a):
if not isinstance(a, types.Array):
return
isnan = get_isnan(a.dtype)
def nanvar_impl(a):
# Compute the mean
m = np.nanmean(a)
# Compute the sum of square diffs
ssd = 0.0
count = 0
for view in np.nditer(a):
v = view.item()
if not isnan(v):
val = (v.item() - m)
ssd += np.real(val * np.conj(val))
count += 1
# np.divide() doesn't raise ZeroDivisionError
return np.divide(ssd, count)
return nanvar_impl
@overload(np.nanstd)
def np_nanstd(a):
if not isinstance(a, types.Array):
return
def nanstd_impl(a):
return np.nanvar(a) ** 0.5
return nanstd_impl
@overload(np.nansum)
def np_nansum(a):
if not isinstance(a, types.Array):
return
if isinstance(a.dtype, types.Integer):
retty = types.intp
else:
retty = a.dtype
zero = retty(0)
isnan = get_isnan(a.dtype)
def nansum_impl(a):
c = zero
for view in np.nditer(a):
v = view.item()
if not isnan(v):
c += v
return c
return nansum_impl
@overload(np.nanprod)
def np_nanprod(a):
if not isinstance(a, types.Array):
return
if isinstance(a.dtype, types.Integer):
retty = types.intp
else:
retty = a.dtype
one = retty(1)
isnan = get_isnan(a.dtype)
def nanprod_impl(a):
c = one
for view in np.nditer(a):
v = view.item()
if not isnan(v):
c *= v
return c
return nanprod_impl
@overload(np.nancumprod)
def np_nancumprod(a):
if not isinstance(a, types.Array):
return
if isinstance(a.dtype, (types.Boolean, types.Integer)):
# dtype cannot possibly contain NaN
return lambda a: np.cumprod(a)
else:
retty = a.dtype
is_nan = get_isnan(retty)
one = retty(1)
def nancumprod_impl(a):
out = np.empty(a.size, retty)
c = one
for idx, v in enumerate(a.flat):
if ~is_nan(v):
c *= v
out[idx] = c
return out
return nancumprod_impl
@overload(np.nancumsum)
def np_nancumsum(a):
if not isinstance(a, types.Array):
return
if isinstance(a.dtype, (types.Boolean, types.Integer)):
# dtype cannot possibly contain NaN
return lambda a: np.cumsum(a)
else:
retty = a.dtype
is_nan = get_isnan(retty)
zero = retty(0)
def nancumsum_impl(a):
out = np.empty(a.size, retty)
c = zero
for idx, v in enumerate(a.flat):
if ~is_nan(v):
c += v
out[idx] = c
return out
return nancumsum_impl
@register_jitable
def prepare_ptp_input(a):
arr = _asarray(a)
if len(arr) == 0:
raise ValueError('zero-size array reduction not possible')
else:
return arr
def _compute_current_val_impl_gen(op):
def _compute_current_val_impl(current_val, val):
if isinstance(current_val, types.Complex):
# The sort order for complex numbers is lexicographic. If both the
# real and imaginary parts are non-nan then the order is determined
# by the real parts except when they are equal, in which case the
# order is determined by the imaginary parts.
# https://github.com/numpy/numpy/blob/577a86e/numpy/core/fromnumeric.py#L874-L877 # noqa: E501
def impl(current_val, val):
if op(val.real, current_val.real):
return val
elif (val.real == current_val.real
and op(val.imag, current_val.imag)):
return val
return current_val
else:
def impl(current_val, val):
return val if op(val, current_val) else current_val
return impl
return _compute_current_val_impl
_compute_a_max = generated_jit(_compute_current_val_impl_gen(greater_than))
_compute_a_min = generated_jit(_compute_current_val_impl_gen(less_than))
@generated_jit
def _early_return(val):
UNUSED = 0
if isinstance(val, types.Complex):
def impl(val):
if np.isnan(val.real):
if np.isnan(val.imag):
return True, np.nan + np.nan * 1j
else:
return True, np.nan + 0j
else:
return False, UNUSED
elif isinstance(val, types.Float):
def impl(val):
if np.isnan(val):
return True, np.nan
else:
return False, UNUSED
else:
def impl(val):
return False, UNUSED
return impl
@overload_method(types.Array, 'ptp')
@overload(np.ptp)
def np_ptp(a):
if hasattr(a, 'dtype'):
if isinstance(a.dtype, types.Boolean):
raise TypingError("Boolean dtype is unsupported (as per NumPy)")
# Numpy raises a TypeError
def np_ptp_impl(a):
arr = prepare_ptp_input(a)
a_flat = arr.flat
a_min = a_flat[0]
a_max = a_flat[0]
for i in range(arr.size):
val = a_flat[i]
take_branch, retval = _early_return(val)
if take_branch:
return retval
a_max = _compute_a_max(a_max, val)
a_min = _compute_a_min(a_min, val)
return a_max - a_min
return np_ptp_impl
#----------------------------------------------------------------------------
# Median and partitioning
@register_jitable
def nan_aware_less_than(a, b):
if np.isnan(a):
return False
else:
if np.isnan(b):
return True
else:
return a < b
def _partition_factory(pivotimpl):
def _partition(A, low, high):
mid = (low + high) >> 1
# NOTE: the pattern of swaps below for the pivot choice and the
# partitioning gives good results (i.e. regular O(n log n))
# on sorted, reverse-sorted, and uniform arrays. Subtle changes
# risk breaking this property.
# Use median of three {low, middle, high} as the pivot
if pivotimpl(A[mid], A[low]):
A[low], A[mid] = A[mid], A[low]
if pivotimpl(A[high], A[mid]):
A[high], A[mid] = A[mid], A[high]
if pivotimpl(A[mid], A[low]):
A[low], A[mid] = A[mid], A[low]
pivot = A[mid]
A[high], A[mid] = A[mid], A[high]
i = low
j = high - 1
while True:
while i < high and pivotimpl(A[i], pivot):
i += 1
while j >= low and pivotimpl(pivot, A[j]):
j -= 1
if i >= j:
break
A[i], A[j] = A[j], A[i]
i += 1
j -= 1
# Put the pivot back in its final place (all items before `i`
# are smaller than the pivot, all items at/after `i` are larger)
A[i], A[high] = A[high], A[i]
return i
return _partition
_partition = register_jitable(_partition_factory(less_than))
_partition_w_nan = register_jitable(_partition_factory(nan_aware_less_than))
def _select_factory(partitionimpl):
def _select(arry, k, low, high):
"""
Select the k'th smallest element in array[low:high + 1].
"""
i = partitionimpl(arry, low, high)
while i != k:
if i < k:
low = i + 1
i = partitionimpl(arry, low, high)
else:
high = i - 1
i = partitionimpl(arry, low, high)
return arry[k]
return _select
_select = register_jitable(_select_factory(_partition))
_select_w_nan = register_jitable(_select_factory(_partition_w_nan))
@register_jitable
def _select_two(arry, k, low, high):
"""
Select the k'th and k+1'th smallest elements in array[low:high + 1].
This is significantly faster than doing two independent selections
for k and k+1.
"""
while True:
assert high > low # by construction
i = _partition(arry, low, high)
if i < k:
low = i + 1
elif i > k + 1:
high = i - 1
elif i == k:
_select(arry, k + 1, i + 1, high)
break
else: # i == k + 1
_select(arry, k, low, i - 1)
break
return arry[k], arry[k + 1]
@register_jitable
def _median_inner(temp_arry, n):
"""
The main logic of the median() call. *temp_arry* must be disposable,
as this function will mutate it.
"""
low = 0
high = n - 1
half = n >> 1
if n & 1 == 0:
a, b = _select_two(temp_arry, half - 1, low, high)
return (a + b) / 2
else:
return _select(temp_arry, half, low, high)
@overload(np.median)
def np_median(a):
if not isinstance(a, types.Array):
return
def median_impl(a):
# np.median() works on the flattened array, and we need a temporary
# workspace anyway
temp_arry = a.flatten()
n = temp_arry.shape[0]
return _median_inner(temp_arry, n)
return median_impl
@register_jitable
def _collect_percentiles_inner(a, q):
n = len(a)
if n == 1:
# single element array; output same for all percentiles
out = np.full(len(q), a[0], dtype=np.float64)
else:
out = np.empty(len(q), dtype=np.float64)
for i in range(len(q)):
percentile = q[i]
# bypass pivoting where requested percentile is 100
if percentile == 100:
val = np.max(a)
# heuristics to handle infinite values a la NumPy
if ~np.all(np.isfinite(a)):
if ~np.isfinite(val):
val = np.nan
# bypass pivoting where requested percentile is 0
elif percentile == 0:
val = np.min(a)
# convoluted heuristics to handle infinite values a la NumPy
if ~np.all(np.isfinite(a)):
num_pos_inf = np.sum(a == np.inf)
num_neg_inf = np.sum(a == -np.inf)
num_finite = n - (num_neg_inf + num_pos_inf)
if num_finite == 0:
val = np.nan
if num_pos_inf == 1 and n == 2:
val = np.nan
if num_neg_inf > 1:
val = np.nan
if num_finite == 1:
if num_pos_inf > 1:
if num_neg_inf != 1:
val = np.nan
else:
# linear interp between closest ranks
rank = 1 + (n - 1) * np.true_divide(percentile, 100.0)
f = math.floor(rank)
m = rank - f
lower, upper = _select_two(a, k=int(f - 1), low=0, high=(n - 1))
val = lower * (1 - m) + upper * m
out[i] = val
return out
@register_jitable
def _can_collect_percentiles(a, nan_mask, skip_nan):
if skip_nan:
a = a[~nan_mask]
if len(a) == 0:
return False # told to skip nan, but no elements remain
else:
if np.any(nan_mask):
return False # told *not* to skip nan, but nan encountered
if len(a) == 1: # single element array
val = a[0]
return np.isfinite(val) # can collect percentiles if element is finite
else:
return True
@register_jitable
def check_valid(q, q_upper_bound):
valid = True
# avoid expensive reductions where possible
if q.ndim == 1 and q.size < 10:
for i in range(q.size):
if q[i] < 0.0 or q[i] > q_upper_bound or np.isnan(q[i]):
valid = False
break
else:
if np.any(np.isnan(q)) or np.any(q < 0.0) or np.any(q > q_upper_bound):
valid = False
return valid
@register_jitable
def percentile_is_valid(q):
if not check_valid(q, q_upper_bound=100.0):
raise ValueError('Percentiles must be in the range [0, 100]')
@register_jitable
def quantile_is_valid(q):
if not check_valid(q, q_upper_bound=1.0):
raise ValueError('Quantiles must be in the range [0, 1]')
@register_jitable
def _collect_percentiles(a, q, check_q, factor, skip_nan):
q = np.asarray(q, dtype=np.float64).flatten()
check_q(q)
q = q * factor
temp_arry = np.asarray(a, dtype=np.float64).flatten()
nan_mask = np.isnan(temp_arry)
if _can_collect_percentiles(temp_arry, nan_mask, skip_nan):
temp_arry = temp_arry[~nan_mask]
out = _collect_percentiles_inner(temp_arry, q)
else:
out = np.full(len(q), np.nan)
return out
def _percentile_quantile_inner(a, q, skip_nan, factor, check_q):
"""
The underlying algorithm to find percentiles and quantiles
is the same, hence we converge onto the same code paths
in this inner function implementation
"""
dt = determine_dtype(a)
if np.issubdtype(dt, np.complexfloating):
raise TypingError('Not supported for complex dtype')
# this could be supported, but would require a
# lexicographic comparison
def np_percentile_q_scalar_impl(a, q):
return _collect_percentiles(a, q, check_q, factor, skip_nan)[0]
def np_percentile_impl(a, q):
return _collect_percentiles(a, q, check_q, factor, skip_nan)
if isinstance(q, (types.Number, types.Boolean)):
return np_percentile_q_scalar_impl
elif isinstance(q, types.Array) and q.ndim == 0:
return np_percentile_q_scalar_impl
else:
return np_percentile_impl
@overload(np.percentile)
def np_percentile(a, q):
return _percentile_quantile_inner(
a, q, skip_nan=False, factor=1.0, check_q=percentile_is_valid
)
@overload(np.nanpercentile)
def np_nanpercentile(a, q):
return _percentile_quantile_inner(
a, q, skip_nan=True, factor=1.0, check_q=percentile_is_valid
)
@overload(np.quantile)
def np_quantile(a, q):
return _percentile_quantile_inner(
a, q, skip_nan=False, factor=100.0, check_q=quantile_is_valid
)
@overload(np.nanquantile)
def np_nanquantile(a, q):
return _percentile_quantile_inner(
a, q, skip_nan=True, factor=100.0, check_q=quantile_is_valid
)
@overload(np.nanmedian)
def np_nanmedian(a):
if not isinstance(a, types.Array):
return
isnan = get_isnan(a.dtype)
def nanmedian_impl(a):
# Create a temporary workspace with only non-NaN values
temp_arry = np.empty(a.size, a.dtype)
n = 0
for view in np.nditer(a):
v = view.item()
if not isnan(v):
temp_arry[n] = v
n += 1
# all NaNs
if n == 0:
return np.nan
return _median_inner(temp_arry, n)
return nanmedian_impl
@register_jitable
def np_partition_impl_inner(a, kth_array):
# allocate and fill empty array rather than copy a and mutate in place
# as the latter approach fails to preserve strides
out = np.empty_like(a)
idx = np.ndindex(a.shape[:-1]) # Numpy default partition axis is -1
for s in idx:
arry = a[s].copy()
low = 0
high = len(arry) - 1
for kth in kth_array:
_select_w_nan(arry, kth, low, high)
low = kth # narrow span of subsequent partition
out[s] = arry
return out
@register_jitable
def valid_kths(a, kth):
"""
Returns a sorted, unique array of kth values which serve
as indexers for partitioning the input array, a.
If the absolute value of any of the provided values
is greater than a.shape[-1] an exception is raised since
we are partitioning along the last axis (per Numpy default
behaviour).
Values less than 0 are transformed to equivalent positive
index values.
"""
# cast boolean to int, where relevant
kth_array = _asarray(kth).astype(np.int64)
if kth_array.ndim != 1:
raise ValueError('kth must be scalar or 1-D')
# numpy raises ValueError: object too deep for desired array
if np.any(np.abs(kth_array) >= a.shape[-1]):
raise ValueError("kth out of bounds")
out = np.empty_like(kth_array)
for index, val in np.ndenumerate(kth_array):
if val < 0:
out[index] = val + a.shape[-1] # equivalent positive index
else:
out[index] = val
return np.unique(out)
@overload(np.partition)
def np_partition(a, kth):
if not isinstance(a, (types.Array, types.Sequence, types.Tuple)):
raise TypeError('The first argument must be an array-like')
if isinstance(a, types.Array) and a.ndim == 0:
raise TypeError('The first argument must be at least 1-D (found 0-D)')
kthdt = getattr(kth, 'dtype', kth)
if not isinstance(kthdt, (types.Boolean, types.Integer)):
# bool gets cast to int subsequently
raise TypeError('Partition index must be integer')
def np_partition_impl(a, kth):
a_tmp = _asarray(a)
if a_tmp.size == 0:
return a_tmp.copy()
else:
kth_array = valid_kths(a_tmp, kth)
return np_partition_impl_inner(a_tmp, kth_array)
return np_partition_impl
#----------------------------------------------------------------------------
# Building matrices
@register_jitable
def _tri_impl(N, M, k):
shape = max(0, N), max(0, M) # numpy floors each dimension at 0
out = np.empty(shape, dtype=np.float64) # numpy default dtype
for i in range(shape[0]):
m_max = min(max(0, i + k + 1), shape[1])
out[i, :m_max] = 1
out[i, m_max:] = 0
return out
@overload(np.tri)
def np_tri(N, M=None, k=0):
# we require k to be integer, unlike numpy
_check_is_integer(k, 'k')
def tri_impl(N, M=None, k=0):
if M is None:
M = N
return _tri_impl(N, M, k)
return tri_impl
@register_jitable
def _make_square(m):
"""
Takes a 1d array and tiles it to form a square matrix
- i.e. a facsimile of np.tile(m, (len(m), 1))
"""
assert m.ndim == 1
len_m = len(m)
out = np.empty((len_m, len_m), dtype=m.dtype)
for i in range(len_m):
out[i] = m
return out
@register_jitable
def np_tril_impl_2d(m, k=0):
mask = np.tri(m.shape[-2], M=m.shape[-1], k=k).astype(np.uint)
return np.where(mask, m, np.zeros_like(m, dtype=m.dtype))
@overload(np.tril)
def my_tril(m, k=0):
# we require k to be integer, unlike numpy
_check_is_integer(k, 'k')
def np_tril_impl_1d(m, k=0):
m_2d = _make_square(m)
return np_tril_impl_2d(m_2d, k)
def np_tril_impl_multi(m, k=0):
mask = np.tri(m.shape[-2], M=m.shape[-1], k=k).astype(np.uint)
idx = np.ndindex(m.shape[:-2])
z = np.empty_like(m)
zero_opt = np.zeros_like(mask, dtype=m.dtype)
for sel in idx:
z[sel] = np.where(mask, m[sel], zero_opt)
return z
if m.ndim == 1:
return np_tril_impl_1d
elif m.ndim == 2:
return np_tril_impl_2d
else:
return np_tril_impl_multi
@overload(np.tril_indices)
def np_tril_indices(n, k=0, m=None):
# we require integer arguments, unlike numpy
_check_is_integer(n, 'n')
_check_is_integer(k, 'k')
if not is_nonelike(m):
_check_is_integer(m, 'm')
def np_tril_indices_impl(n, k=0, m=None):
return np.nonzero(np.tri(n, m, k=k))
return np_tril_indices_impl
@overload(np.tril_indices_from)
def np_tril_indices_from(arr, k=0):
# we require k to be integer, unlike numpy
_check_is_integer(k, 'k')
if arr.ndim != 2:
raise TypingError("input array must be 2-d")
def np_tril_indices_from_impl(arr, k=0):
return np.tril_indices(arr.shape[0], k=k, m=arr.shape[1])
return np_tril_indices_from_impl
@register_jitable
def np_triu_impl_2d(m, k=0):
mask = np.tri(m.shape[-2], M=m.shape[-1], k=k - 1).astype(np.uint)
return np.where(mask, np.zeros_like(m, dtype=m.dtype), m)
@overload(np.triu)
def my_triu(m, k=0):
# we require k to be integer, unlike numpy
_check_is_integer(k, 'k')
def np_triu_impl_1d(m, k=0):
m_2d = _make_square(m)
return np_triu_impl_2d(m_2d, k)
def np_triu_impl_multi(m, k=0):
mask = np.tri(m.shape[-2], M=m.shape[-1], k=k - 1).astype(np.uint)
idx = np.ndindex(m.shape[:-2])
z = np.empty_like(m)
zero_opt = np.zeros_like(mask, dtype=m.dtype)
for sel in idx:
z[sel] = np.where(mask, zero_opt, m[sel])
return z
if m.ndim == 1:
return np_triu_impl_1d
elif m.ndim == 2:
return np_triu_impl_2d
else:
return np_triu_impl_multi
@overload(np.triu_indices)
def np_triu_indices(n, k=0, m=None):
# we require integer arguments, unlike numpy
_check_is_integer(n, 'n')
_check_is_integer(k, 'k')
if not is_nonelike(m):
_check_is_integer(m, 'm')
def np_triu_indices_impl(n, k=0, m=None):
return np.nonzero(1 - np.tri(n, m, k=k - 1))
return np_triu_indices_impl
@overload(np.triu_indices_from)
def np_triu_indices_from(arr, k=0):
# we require k to be integer, unlike numpy
_check_is_integer(k, 'k')
if arr.ndim != 2:
raise TypingError("input array must be 2-d")
def np_triu_indices_from_impl(arr, k=0):
return np.triu_indices(arr.shape[0], k=k, m=arr.shape[1])
return np_triu_indices_from_impl
def _prepare_array(arr):
pass
@overload(_prepare_array)
def _prepare_array_impl(arr):
if arr in (None, types.none):
return lambda arr: np.array(())
else:
return lambda arr: _asarray(arr).ravel()
def _dtype_of_compound(inobj):
obj = inobj
while True:
if isinstance(obj, (types.Number, types.Boolean)):
return as_dtype(obj)
l = getattr(obj, '__len__', None)
if l is not None and l() == 0: # empty tuple or similar
return np.float64
dt = getattr(obj, 'dtype', None)
if dt is None:
raise TypeError("type has no dtype attr")
if isinstance(obj, types.Sequence):
obj = obj.dtype
else:
return as_dtype(dt)
@overload(np.ediff1d)
def np_ediff1d(ary, to_end=None, to_begin=None):
if isinstance(ary, types.Array):
if isinstance(ary.dtype, types.Boolean):
raise TypeError("Boolean dtype is unsupported (as per NumPy)")
# Numpy tries to do this: return ary[1:] - ary[:-1] which
# results in a TypeError exception being raised
# since np 1.16 there are casting checks for to_end and to_begin to make
# sure they are compatible with the ary
if numpy_version >= (1, 16):
ary_dt = _dtype_of_compound(ary)
to_begin_dt = None
if not(is_nonelike(to_begin)):
to_begin_dt = _dtype_of_compound(to_begin)
to_end_dt = None
if not(is_nonelike(to_end)):
to_end_dt = _dtype_of_compound(to_end)
if to_begin_dt is not None and not np.can_cast(to_begin_dt, ary_dt):
msg = "dtype of to_begin must be compatible with input ary"
raise TypeError(msg)
if to_end_dt is not None and not np.can_cast(to_end_dt, ary_dt):
msg = "dtype of to_end must be compatible with input ary"
raise TypeError(msg)
def np_ediff1d_impl(ary, to_end=None, to_begin=None):
# transform each input into an equivalent 1d array
start = _prepare_array(to_begin)
mid = _prepare_array(ary)
end = _prepare_array(to_end)
out_dtype = mid.dtype
# output array dtype determined by ary dtype, per NumPy
# (for the most part); an exception to the rule is a zero length
# array-like, where NumPy falls back to np.float64; this behaviour
# is *not* replicated
if len(mid) > 0:
out = np.empty((len(start) + len(mid) + len(end) - 1),
dtype=out_dtype)
start_idx = len(start)
mid_idx = len(start) + len(mid) - 1
out[:start_idx] = start
out[start_idx:mid_idx] = np.diff(mid)
out[mid_idx:] = end
else:
out = np.empty((len(start) + len(end)), dtype=out_dtype)
start_idx = len(start)
out[:start_idx] = start
out[start_idx:] = end
return out
return np_ediff1d_impl
def _select_element(arr):
pass
@overload(_select_element)
def _select_element_impl(arr):
zerod = getattr(arr, 'ndim', None) == 0
if zerod:
def impl(arr):
x = np.array((1,), dtype=arr.dtype)
x[:] = arr
return x[0]
return impl
else:
def impl(arr):
return arr
return impl
def _get_d(dx, x):
pass
@overload(_get_d)
def get_d_impl(x, dx):
if is_nonelike(x):
def impl(x, dx):
return np.asarray(dx)
else:
def impl(x, dx):
return np.diff(np.asarray(x))
return impl
@overload(np.trapz)
def np_trapz(y, x=None, dx=1.0):
if isinstance(y, (types.Number, types.Boolean)):
raise TypingError('y cannot be a scalar')
elif isinstance(y, types.Array) and y.ndim == 0:
raise TypingError('y cannot be 0D')
# NumPy raises IndexError: list assignment index out of range
# inspired by:
# https://github.com/numpy/numpy/blob/7ee52003/numpy/lib/function_base.py#L4040-L4065 # noqa: E501
def impl(y, x=None, dx=1.0):
yarr = np.asarray(y)
d = _get_d(x, dx)
y_ave = (yarr[..., slice(1, None)] + yarr[..., slice(None, -1)]) / 2.0
ret = np.sum(d * y_ave, -1)
processed = _select_element(ret)
return processed
return impl
@register_jitable
def _np_vander(x, N, increasing, out):
"""
Generate an N-column Vandermonde matrix from a supplied 1-dimensional
array, x. Store results in an output matrix, out, which is assumed to
be of the required dtype.
Values are accumulated using np.multiply to match the floating point
precision behaviour of numpy.vander.
"""
m, n = out.shape
assert m == len(x)
assert n == N
if increasing:
for i in range(N):
if i == 0:
out[:, i] = 1
else:
out[:, i] = np.multiply(x, out[:, (i - 1)])
else:
for i in range(N - 1, -1, -1):
if i == N - 1:
out[:, i] = 1
else:
out[:, i] = np.multiply(x, out[:, (i + 1)])
@register_jitable
def _check_vander_params(x, N):
if x.ndim > 1:
raise ValueError('x must be a one-dimensional array or sequence.')
if N < 0:
raise ValueError('Negative dimensions are not allowed')
@overload(np.vander)
def np_vander(x, N=None, increasing=False):
if N not in (None, types.none):
if not isinstance(N, types.Integer):
raise TypingError('Second argument N must be None or an integer')
def np_vander_impl(x, N=None, increasing=False):
if N is None:
N = len(x)
_check_vander_params(x, N)
# allocate output matrix using dtype determined in closure
out = np.empty((len(x), int(N)), dtype=dtype)
_np_vander(x, N, increasing, out)
return out
def np_vander_seq_impl(x, N=None, increasing=False):
if N is None:
N = len(x)
x_arr = np.array(x)
_check_vander_params(x_arr, N)
# allocate output matrix using dtype inferred when x_arr was created
out = np.empty((len(x), int(N)), dtype=x_arr.dtype)