use 64-bit safe types

This commit is contained in:
Zach Pincus
2011-11-07 15:51:39 -05:00
parent 4c88f94508
commit ea63e7d709
4 changed files with 141 additions and 123 deletions
+7 -7
View File
@@ -3,17 +3,17 @@ It contains the definitions of the mcp class, such that
other cython modules can "cimport mcp" and subclass it.
"""
cimport numpy as np
cimport heap
cimport numpy as np
# determine datatypes for mcp
ctypedef heap.BOOL_T BOOL_T
ctypedef unsigned char DIM_T
ctypedef np.float64_t FLOAT_T
ctypedef double FLOAT_C
cdef class MCP:
cdef heap.FastUpdateBinaryHeap costs_heap
cdef object costs_shape
cdef int dim
cdef DIM_T dim
cdef object flat_costs
cdef object flat_cumulative_costs
cdef object traceback_offsets
@@ -21,9 +21,9 @@ cdef class MCP:
cdef readonly object offsets
cdef object flat_offsets
cdef object offset_lengths
cdef int dirty
cdef int use_start_cost
cdef BOOL_T dirty
cdef BOOL_T use_start_cost
# if use_start_cost is true, the cost of the starting element is added to
# the cost of the path. Set to true by default in the base class...
cdef FLOAT_C _travel_cost(self, FLOAT_C old_cost, FLOAT_C new_cost, FLOAT_C offset_length)
cdef FLOAT_T _travel_cost(self, FLOAT_T old_cost, FLOAT_T new_cost, FLOAT_T offset_length)
+54 -40
View File
@@ -38,7 +38,19 @@ import numpy as np
cimport heap
import heap
FLOAT = np.float64
ctypedef np.int8_t OFFSET_T
OFFSET_D = np.int8
ctypedef np.int32_t FLAT_OFFSET_T
FLAT_OFFSET_D = np.int32
ctypedef np.int16_t OFFSETS_INDEX_T
OFFSETS_INDEX_D = np.int16
ctypedef np.int8_t EDGE_T
EDGE_D = np.int8
ctypedef np.intp_t INDEX_T
INDEX_D = np.intp
FLOAT_D = np.float64
def _get_edge_map(shape):
"""Return an array with edge points/lines/planes/hyperplanes marked.
@@ -53,7 +65,7 @@ def _get_edge_map(shape):
"""
d = len(shape)
edges = np.zeros(shape+(d,), order='F', dtype=np.int8)
edges = np.zeros(shape+(d,), order='F', dtype=EDGE_D)
for i in range(d):
slices = [slice(None)] * (d+1)
slices[d] = i
@@ -90,7 +102,7 @@ def _offset_edge_map(shape, offsets):
"""
d = len(shape)
edges = np.zeros(shape+(d,), order='F', dtype=np.int8)
edges = np.zeros(shape+(d,), order='F', dtype=EDGE_D)
offsets = np.asarray(offsets)
for i in range(d):
slices = [slice(None)] * (d+1)
@@ -244,16 +256,16 @@ cdef class MCP:
See class documentation.
"""
costs = np.asarray(costs)
if not np.can_cast(costs.dtype, FLOAT):
raise TypeError('cannot cast costs array to ' + str(FLOAT))
if not np.can_cast(costs.dtype, FLOAT_D):
raise TypeError('cannot cast costs array to ' + str(FLOAT_D))
# We use flat, fortran-style indexing here (could use C-style,
# but this is my code and I like fortran-style! Also, it's
# faster when working with image arrays, which are often
# already fortran-strided.)
self.flat_costs = costs.astype(FLOAT).flatten('F')
self.flat_costs = costs.astype(FLOAT_D).flatten('F')
size = self.flat_costs.shape[0]
self.flat_cumulative_costs = np.empty(size, dtype=FLOAT)
self.flat_cumulative_costs = np.empty(size, dtype=FLOAT_D)
self.flat_cumulative_costs.fill(np.inf)
self.dim = len(costs.shape)
self.costs_shape = costs.shape
@@ -263,7 +275,7 @@ cdef class MCP:
# This array stores, for each point, the index into the offset
# array (see below) that leads to that point from the
# predecessor point.
self.traceback_offsets = np.empty(size, dtype=np.int16)
self.traceback_offsets = np.empty(size, dtype=OFFSETS_INDEX_D)
self.traceback_offsets.fill(-1)
# The offsets are a list of relative offsets from a central
@@ -273,10 +285,10 @@ cdef class MCP:
# in the same way for flat indices to move to neighboring points.
if offsets is None:
offsets = make_offsets(self.dim, fully_connected)
self.offsets = np.array(offsets, dtype=np.int8)
self.offsets = np.array(offsets, dtype=OFFSET_D)
self.flat_offsets = np.array(
_ravel_index_fortran(self.offsets, self.costs_shape),
dtype=np.int32)
dtype=FLAT_OFFSET_D)
# Instead of unraveling each index during the pathfinding algorithm, we
# will use a pre-computed "edge map" that specifies for each dimension
@@ -292,7 +304,7 @@ cdef class MCP:
# The offset lengths are the distances traveled along each offset
self.offset_lengths = np.sqrt(
np.sum(self.offsets**2, axis=1)).astype(FLOAT)
np.sum(self.offsets**2, axis=1)).astype(FLOAT_D)
self.dirty = 0
self.use_start_cost = 1
@@ -306,8 +318,8 @@ cdef class MCP:
self.flat_cumulative_costs.fill(np.inf)
self.dirty = 0
cdef FLOAT_C _travel_cost(self, FLOAT_C old_cost,
FLOAT_C new_cost, FLOAT_C offset_length):
cdef FLOAT_T _travel_cost(self, FLOAT_T old_cost,
FLOAT_T new_cost, FLOAT_T offset_length):
return new_cost
@cython.boundscheck(False)
@@ -356,10 +368,10 @@ cdef class MCP:
"""
# basic variables to use for end-finding; also fix up the start and end
# lists
cdef int use_ends = 0
cdef int num_ends
cdef int all_ends = find_all_ends
cdef np.ndarray[np.uint32_t, ndim=1] flat_ends
cdef BOOL_T use_ends = 0
cdef INDEX_T num_ends
cdef BOOL_T all_ends = find_all_ends
cdef np.ndarray[INDEX_T, ndim=1] flat_ends
starts = _normalize_indices(starts, self.costs_shape)
if starts is None:
raise ValueError('start points must all be within the costs array')
@@ -371,7 +383,7 @@ cdef class MCP:
use_ends = 1
num_ends = len(ends)
flat_ends = np.array(_ravel_index_fortran(
ends, self.costs_shape), dtype=np.uint32)
ends, self.costs_shape), dtype=INDEX_D)
if self.dirty:
self._reset()
@@ -381,14 +393,14 @@ cdef class MCP:
cdef np.ndarray[FLOAT_T, ndim=1] flat_costs = self.flat_costs
cdef np.ndarray[FLOAT_T, ndim=1] flat_cumulative_costs = \
self.flat_cumulative_costs
cdef np.ndarray[np.int16_t, ndim=1] traceback_offsets = \
cdef np.ndarray[OFFSETS_INDEX_T, ndim=1] traceback_offsets = \
self.traceback_offsets
cdef np.ndarray[np.int8_t, ndim=2] flat_edge_map = self.flat_edge_map
cdef np.ndarray[np.int8_t, ndim=2] offsets = self.offsets
cdef np.ndarray[np.int32_t, ndim=1] flat_offsets = self.flat_offsets
cdef np.ndarray[EDGE_T, ndim=2] flat_edge_map = self.flat_edge_map
cdef np.ndarray[OFFSET_T, ndim=2] offsets = self.offsets
cdef np.ndarray[FLAT_OFFSET_T, ndim=1] flat_offsets = self.flat_offsets
cdef np.ndarray[FLOAT_T, ndim=1] offset_lengths = self.offset_lengths
cdef int dim = self.dim
cdef DIM_T dim = self.dim
cdef int num_offsets = len(flat_offsets)
# push each start point into the heap. Note that we use flat indexing!
@@ -398,14 +410,15 @@ cdef class MCP:
else:
costs_heap.push_fast(0, start)
cdef double cost, new_cost
cdef unsigned int index, new_index
cdef int is_at_edge, use_offset
cdef unsigned int d, i
cdef int offset, edge_val
cdef FLOAT_T cost, new_cost
cdef INDEX_T index, new_index
cdef BOOL_T is_at_edge, use_offset
cdef INDEX_T d, i
cdef OFFSET_T offset
cdef EDGE_T edge_val
cdef int num_ends_found = 0
cdef double inf = np.inf
cdef double travel_cost
cdef FLOAT_T inf = np.inf
cdef FLOAT_T travel_cost
while 1:
# Find the point with the minimum cost in the heap. Once
# popped, this point's minimum cost path has been found.
@@ -540,21 +553,22 @@ cdef class MCP:
'within the costs array')
traceback = [tuple(ends[0])]
cdef unsigned int flat_position =\
cdef INDEX_T flat_position =\
_ravel_index_fortran(ends, self.costs_shape)[0]
if self.flat_cumulative_costs[flat_position] == np.inf:
raise ValueError('no minimum-cost path was found '
'to the specified end point')
cdef np.ndarray[np.int32_t, ndim=1] position = \
np.array(ends[0], dtype=np.int32)
cdef np.ndarray[np.int16_t, ndim=1] traceback_offsets = \
cdef np.ndarray[INDEX_T, ndim=1] position = \
np.array(ends[0], dtype=INDEX_D)
cdef np.ndarray[OFFSETS_INDEX_T, ndim=1] traceback_offsets = \
self.traceback_offsets
cdef np.ndarray[np.int8_t, ndim=2] offsets = self.offsets
cdef np.ndarray[np.int32_t, ndim=1] flat_offsets = self.flat_offsets
cdef np.ndarray[OFFSET_T, ndim=2] offsets = self.offsets
cdef np.ndarray[FLAT_OFFSET_T, ndim=1] flat_offsets = self.flat_offsets
cdef unsigned int offset, d
cdef int dim = self.dim
cdef OFFSETS_INDEX_T offset
cdef DIM_T d
cdef DIM_T dim = self.dim
while 1:
offset = traceback_offsets[flat_position]
if offset == -1:
@@ -601,6 +615,6 @@ cdef class MCP_Geometric(MCP):
raise ValueError('all offset components must be 0, 1, or -1')
self.use_start_cost = 0
cdef FLOAT_C _travel_cost(self, FLOAT_C old_cost, FLOAT_C new_cost,
FLOAT_C offset_length):
cdef FLOAT_T _travel_cost(self, FLOAT_T old_cost, FLOAT_T new_cost,
FLOAT_T offset_length):
return offset_length * 0.5 * (old_cost + new_cost)
+18 -14
View File
@@ -7,29 +7,33 @@ value_of_fast()
# determine datatypes for heap
ctypedef double VALUE_T
ctypedef int REFERENCE_T
ctypedef Py_ssize_t REFERENCE_T
ctypedef REFERENCE_T INDEX_T
ctypedef unsigned char BOOL_T
ctypedef unsigned char LEVELS_T
cdef class BinaryHeap:
cdef readonly int count, levels, min_levels
cdef readonly INDEX_T count
cdef readonly LEVELS_T levels, min_levels
cdef VALUE_T *_values
cdef REFERENCE_T *_references
cdef int _popped_ref
cdef REFERENCE_T _popped_ref
cdef void _add_or_remove_level(self, int add_or_remove)
cdef void _add_or_remove_level(self, LEVELS_T add_or_remove)
cdef void _update(self)
cdef void _update_one(self, int i)
cdef void _remove(self, int i)
cdef void _update_one(self, INDEX_T i)
cdef void _remove(self, INDEX_T i)
cdef int push_fast(self, double value, int reference)
cdef INDEX_T push_fast(self, VALUE_T value, REFERENCE_T reference)
cdef VALUE_T pop_fast(self)
cdef class FastUpdateBinaryHeap(BinaryHeap):
cdef readonly int max_reference
cdef REFERENCE_T *_crossref
cdef int _invalid_ref
cdef int _pushed
cdef readonly REFERENCE_T max_reference
cdef INDEX_T *_crossref
cdef BOOL_T _invalid_ref
cdef BOOL_T _pushed
cdef VALUE_T value_of_fast(self, int reference)
cdef int push_if_lower_fast(self, double value, int reference)
cdef VALUE_T value_of_fast(self, REFERENCE_T reference)
cdef INDEX_T push_if_lower_fast(self, VALUE_T value,
REFERENCE_T reference)
+62 -62
View File
@@ -43,8 +43,7 @@ cdef extern from "pyport.h":
cdef VALUE_T inf = Py_HUGE_VAL
# this is handy
cdef inline int int_max(int a, int b): return a if a >= b else b
cdef inline int int_min(int a, int b): return a if a <= b else b
cdef inline INDEX_T index_min(INDEX_T a, INDEX_T b): return a if a <= b else b
cdef class BinaryHeap:
@@ -127,9 +126,9 @@ cdef class BinaryHeap:
#
# To calculate the capacity at a certain level:
# 2**l
def __cinit__(self, int initial_capacity=128, *args, **kws):
def __cinit__(self, INDEX_T initial_capacity=128, *args, **kws):
# calc levels from the default capacity
cdef int levels = 0
cdef LEVELS_T levels = 0
while 2**levels < initial_capacity:
levels += 1
# set levels
@@ -139,11 +138,11 @@ cdef class BinaryHeap:
self.count = 0
# allocate arrays
cdef int number = 2**self.levels
cdef INDEX_T number = 2**self.levels
self._values = <VALUE_T *>malloc( 2*number * sizeof(VALUE_T))
self._references = <REFERENCE_T *>malloc(number * sizeof(REFERENCE_T))
def __init__(self, int initial_capacity=128):
def __init__(self, INDEX_T initial_capacity=128):
"""__init__(initial_capacity=128)
Class constructor.
@@ -160,8 +159,8 @@ cdef class BinaryHeap:
Reset the heap to default, empty state.
"""
cdef int number = 2**self.levels
cdef int i
cdef INDEX_T number = 2**self.levels
cdef INDEX_T i
cdef VALUE_T *values = self._values
for i in range(number*2):
values[i] = inf
@@ -174,7 +173,6 @@ cdef class BinaryHeap:
free(self._references)
def __str__(self):
cdef int i0, i, n, level
s = ''
for level in range(1,self.levels+1):
i0 = 2**level-1 # LevelStart
@@ -187,15 +185,15 @@ cdef class BinaryHeap:
## C Maintanance methods
cdef void _add_or_remove_level(self, int add_or_remove):
cdef void _add_or_remove_level(self, LEVELS_T add_or_remove):
# init indexing ints
cdef int i, i1, i2, n
cdef INDEX_T i, i1, i2, n
# new amount of levels
cdef int new_levels = self.levels + add_or_remove
cdef LEVELS_T new_levels = self.levels + add_or_remove
# allocate new arrays
cdef int number = 2**new_levels
cdef INDEX_T number = 2**new_levels
cdef VALUE_T *values
cdef REFERENCE_T *references
values = <VALUE_T *>malloc(number*2 * sizeof(VALUE_T))
@@ -213,7 +211,7 @@ cdef class BinaryHeap:
if self.count:
i1 = 2**new_levels-1 # LevelStart
i2 = 2**self.levels-1 # LevelStart
n = int_min(2**new_levels, 2**self.levels)
n = index_min(2**new_levels, 2**self.levels)
for i in range(n):
values[i1+i] = old_values[i2+i]
for i in range(n):
@@ -238,7 +236,8 @@ cdef class BinaryHeap:
cdef VALUE_T *values = self._values
# Note that i represents an absolute index here
cdef int i0, i, ii, n, level
cdef INDEX_T i0, i, ii, n
cdef LEVELS_T level
# track tree
for level in range(self.levels,1,-1):
@@ -252,7 +251,7 @@ cdef class BinaryHeap:
values[ii] = values[i+1]
cdef void _update_one(self, int i):
cdef void _update_one(self, INDEX_T i):
"""Update the tree for one value."""
# shorter name for values
@@ -263,7 +262,8 @@ cdef class BinaryHeap:
i = i-1
# track tree
cdef int ii, level
cdef INDEX_T ii
cdef LEVELS_T level
for level in range(self.levels,1,-1):
ii = (i-1)//2 # CalcPrevAbs
@@ -279,18 +279,18 @@ cdef class BinaryHeap:
i = ii-1
cdef void _remove(self, int i1):
cdef void _remove(self, INDEX_T i1):
"""Remove a value from the heap. By index."""
cdef int levels = self.levels
cdef int count = self.count
cdef LEVELS_T levels = self.levels
cdef INDEX_T count = self.count
# get indices
cdef int i0 = (1 << levels) - 1 #2**self.levels - 1 # LevelStart
cdef int i2 = i0 + count - 1
cdef INDEX_T i0 = (1 << levels) - 1 #2**self.levels - 1 # LevelStart
cdef INDEX_T i2 = i0 + count - 1
# get relative indices
cdef int r1 = i1 - i0
cdef int r2 = count - 1
cdef INDEX_T r1 = i1 - i0
cdef INDEX_T r2 = count - 1
cdef VALUE_T *values = self._values
cdef REFERENCE_T *references = self._references
@@ -314,19 +314,19 @@ cdef class BinaryHeap:
## C Public methods
cdef int push_fast(self, double value, int reference):
cdef INDEX_T push_fast(self, VALUE_T value, REFERENCE_T reference):
"""The c-method for fast pushing.
Returns the index relative to the start of the last level in the heap."""
# We need to resize if currently it just fits.
cdef int levels = self.levels
cdef int count = self.count
cdef LEVELS_T levels = self.levels
cdef INDEX_T count = self.count
if count >= (1 << levels):#2**self.levels:
self._add_or_remove_level(+1)
levels += 1
# insert new value
cdef int i = ((1 << levels) - 1) + count # LevelStart + n
cdef INDEX_T i = ((1 << levels) - 1) + count # LevelStart + n
self._values[i] = value
self._references[count] = reference
@@ -347,9 +347,9 @@ cdef class BinaryHeap:
cdef VALUE_T *values = self._values
# init index. start at 1 because we start in level 1
cdef int level
cdef int i = 1
cdef int levels = self.levels
cdef LEVELS_T level
cdef INDEX_T i = 1
cdef LEVELS_T levels = self.levels
# search tree (using absolute indices)
for level in range(1, levels):
if values[i] <= values[i+1]:
@@ -364,8 +364,8 @@ cdef class BinaryHeap:
i = i+1
# get values
cdef int ir = i - ((1 << levels) - 1) #(2**self.levels-1) # LevelStart
cdef VALUE_T value = values[i]
cdef INDEX_T ir = i - ((1 << levels) - 1) #(2**self.levels-1) # LevelStart
cdef VALUE_T value = values[i]
self._popped_ref = self._references[ir]
# remove it
@@ -378,7 +378,7 @@ cdef class BinaryHeap:
## Python Public methods (that do not need to be VERY fast)
def push(self, double value, int reference=-1):
def push(self, VALUE_T value, REFERENCE_T reference=-1):
"""push(value, reference=-1)
Append a value to the heap, with optional reference.
@@ -416,10 +416,10 @@ cdef class BinaryHeap:
Get the values in the heap as a list.
"""
out = []
cdef int i, i0
cdef INDEX_T i, i0
i0 = 2**self.levels-1 # LevelStart
for i in range(self.count):
out.append( self._values[i0+i] )
out.append(self._values[i0+i])
return out
@@ -429,9 +429,9 @@ cdef class BinaryHeap:
Get the references in the heap as a list.
"""
out = []
cdef int i
cdef INDEX_T i
for i in range(self.count):
out.append( self._references[i] )
out.append(self._references[i])
return out
@@ -511,14 +511,14 @@ cdef class FastUpdateBinaryHeap(BinaryHeap):
is lower than the current value in the heap. This is again useful for
pathfinding algorithms.
"""
def __cinit__(self, int initial_capacity=128, max_reference=None):
def __cinit__(self, INDEX_T initial_capacity=128, max_reference=None):
if max_reference is None:
max_reference = initial_capacity - 1
self.max_reference = max_reference
self._crossref = <REFERENCE_T *>malloc((max_reference+1) * \
sizeof(REFERENCE_T))
self._crossref = <INDEX_T *>malloc((self.max_reference+1) * \
sizeof(INDEX_T))
def __init__(self, int initial_capacity=128, max_reference=None):
def __init__(self, INDEX_T initial_capacity=128, max_reference=None):
"""__init__(initial_capacity=128, max_reference=None)
Class constructor.
@@ -538,27 +538,27 @@ cdef class FastUpdateBinaryHeap(BinaryHeap):
"""
BinaryHeap.reset(self)
# set default values of crossrefs
cdef int i
cdef INDEX_T i
for i in range(self.max_reference+1):
self._crossref[i] = -1
cdef void _remove(self, int i1):
cdef void _remove(self, INDEX_T i1):
""" Remove a value from the heap. By index. """
cdef int levels = self.levels
cdef int count = self.count
cdef LEVELS_T levels = self.levels
cdef INDEX_T count = self.count
# get indices
cdef int i0 = (1 << levels) - 1 #2**self.levels - 1 # LevelStart
cdef int i2 = i0 + count - 1
cdef INDEX_T i0 = (1 << levels) - 1 #2**self.levels - 1 # LevelStart
cdef INDEX_T i2 = i0 + count - 1
# get relative indices
cdef int r1 = i1 - i0
cdef int r2 = count - 1
cdef INDEX_T r1 = i1 - i0
cdef INDEX_T r2 = count - 1
cdef VALUE_T *values = self._values
cdef REFERENCE_T *references = self._references
cdef REFERENCE_T *crossref = self._crossref
cdef INDEX_T *crossref = self._crossref
# update cross reference
crossref[references[r2]]=r1
@@ -574,14 +574,14 @@ cdef class FastUpdateBinaryHeap(BinaryHeap):
# update
self.count -= 1
count -= 1
if (levels>self.min_levels) & (count < (1 << (levels-2))):
if (levels > self.min_levels) & (count < (1 << (levels-2))):
self._add_or_remove_level(-1)
else:
self._update_one(i1)
self._update_one(i2)
cdef int push_fast(self, double value, int reference):
cdef INDEX_T push_fast(self, VALUE_T value, REFERENCE_T reference):
"""The c method for fast pushing.
If the reference is already present, will update its value, otherwise
@@ -593,11 +593,11 @@ cdef class FastUpdateBinaryHeap(BinaryHeap):
return -1
# init variable to store the index-in-the-heap
cdef int i
cdef INDEX_T i
# Reference is the index in the array where MCP is applied to.
# Find the index-in-the-heap using the crossref array.
cdef int ir = self._crossref[reference]
cdef INDEX_T ir = self._crossref[reference]
if ir != -1:
# update
@@ -611,7 +611,7 @@ cdef class FastUpdateBinaryHeap(BinaryHeap):
self._crossref[reference] = ir
return ir
cdef int push_if_lower_fast(self, double value, int reference):
cdef INDEX_T push_if_lower_fast(self, VALUE_T value, REFERENCE_T reference):
"""If the reference is already present, will update its value ONLY if
the new value is lower than the old one. If the reference is not
present, this append it. If a value was appended, self._pushed is
@@ -624,11 +624,11 @@ cdef class FastUpdateBinaryHeap(BinaryHeap):
return -1
# init variable to store the index-in-the-heap
cdef int i
cdef INDEX_T i
# Reference is the index in the array where MCP is applied to.
# Find the index-in-the-heap using the crossref array.
cdef int ir = self._crossref[reference]
cdef INDEX_T ir = self._crossref[reference]
cdef VALUE_T *values = self._values
self._pushed = 1
if ir != -1:
@@ -647,7 +647,7 @@ cdef class FastUpdateBinaryHeap(BinaryHeap):
return ir
cdef VALUE_T value_of_fast(self, int reference):
cdef VALUE_T value_of_fast(self, REFERENCE_T reference):
"""Return the value corresponding to the given reference. If inf
is returned, the reference may be invalid: check the _invaild_ref
field in this case."""
@@ -657,11 +657,11 @@ cdef class FastUpdateBinaryHeap(BinaryHeap):
return inf
# init variable to store the index-in-the-heap
cdef int i
cdef INDEX_T i
# Reference is the index in the array where MCP is applied to.
# Find the index-in-the-heap using the crossref array.
cdef int ir = self._crossref[reference]
cdef INDEX_T ir = self._crossref[reference]
self._invalid_ref = 0
if ir == -1:
self._invalid_ref = 1
@@ -750,7 +750,7 @@ cdef class FastUpdateBinaryHeap(BinaryHeap):
def cross_references(self):
"""Get the cross references in the heap as a list."""
out = []
cdef int i
cdef INDEX_T i
for i in range(self.max_reference+1):
out.append( self._crossref[i] )
return out