Skip to content

Add support for cupy. #132

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
syamajala wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Add support for cupy. #132

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Cpufit/python/pycpufit/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1 +1 @@
from cpufit import *

211 changes: 210 additions & 1 deletion Gpufit/python/pygpufit/gpufit.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@

import os
import time
from ctypes import cdll, POINTER, byref, c_int, c_float, c_char, c_char_p, c_size_t
from ctypes import cdll, POINTER, byref, c_int, c_float, c_char, c_char_p, c_size_t, cast
import numpy as np

# define library loader (actual loading is lazy)
Expand All @@ -30,6 +30,12 @@
POINTER(c_float), POINTER(c_int), c_float, c_int, POINTER(c_int), c_int, c_size_t,
POINTER(c_char), POINTER(c_float), POINTER(c_int), POINTER(c_float), POINTER(c_int)]

gpufit_func_cuda = lib.gpufit_constrained_cuda_interface
gpufit_func_cuda.restype = c_int
gpufit_func_cuda.argtypes = [c_size_t, c_size_t, POINTER(c_float), POINTER(c_float), c_int, c_float, c_int,
POINTER(c_int), POINTER(c_float), POINTER(c_int), c_int, c_size_t, POINTER(c_char),
POINTER(c_float), POINTER(c_int), POINTER(c_float), POINTER(c_int)]

# gpufit_get_last_error function in the dll
error_func = lib.gpufit_get_last_error
error_func.restype = c_char_p
Expand Down Expand Up @@ -320,3 +326,206 @@ def get_cuda_version():
driver_version = (driver_version // 1000, driver_version % 1000 // 10)

return runtime_version, driver_version

try:
import cupy as cp

def fit_gpu(data, weights, model_id, initial_parameters, tolerance=None, max_number_iterations=None, \
parameters_to_fit=None, estimator_id=None, user_info=None):
"""
Calls the C interface gpufit_cuda_interface function in the library.
(see also http://gpufit.readthedocs.io/en/latest/bindings.html#python)

All 2D NumPy arrays must be in row-major order (standard in NumPy), i.e. array.flags.C_CONTIGUOUS must be True
(see also https://docs.scipy.org/doc/numpy/reference/arrays.ndarray.html#internal-memory-layout-of-an-ndarray)

:param data: The data - 2D NumPy array of dimension [number_fits, number_points] and data type np.float32
:param weights: The weights - 2D NumPy array of the same dimension and data type as parameter data or None (no weights available)
:param model_id: The model ID
:param initial_parameters: Initial values for parameters - NumPy array of dimension [number_fits, number_parameters] and data type np.float32
:param tolerance: The fit tolerance or None (will use default value)
:param max_number_iterations: The maximal number of iterations or None (will use default value)
:param parameters_to_fit: Which parameters to fit - NumPy array of length number_parameters and type np.int32 or None (will fit all parameters)
:param estimator_id: The Estimator ID or None (will use default values)
:param user_info: User info - NumPy array of type np.char or None (no user info available)
:return: parameters, states, chi_squares, number_iterations, execution_time
"""

# call fit_constrained without any constraints
return fit_constrained_gpu(data, weights, model_id, initial_parameters, tolerance=tolerance,
max_number_iterations=max_number_iterations,
parameters_to_fit=parameters_to_fit,
estimator_id=estimator_id, user_info=user_info)


def fit_constrained_gpu(data, weights, model_id, initial_parameters, constraints=None, constraint_types=None,
tolerance=None, max_number_iterations=None, \
parameters_to_fit=None, estimator_id=None, user_info=None):
"""
Calls the C interface gpufit_cuda_interface function in the library.
(see also http://gpufit.readthedocs.io/en/latest/bindings.html#python)

All 2D NumPy arrays must be in row-major order (standard in NumPy), i.e. array.flags.C_CONTIGUOUS must be True
(see also https://docs.scipy.org/doc/numpy/reference/arrays.ndarray.html#internal-memory-layout-of-an-ndarray)

:param data: The data - 2D NumPy array of dimension [number_fits, number_points] and data type np.float32
:param weights: The weights - 2D NumPy array of the same dimension and data type as parameter data or None (no weights available)
:param model_id: The model ID
:param initial_parameters: Initial values for parameters - NumPy array of dimension [number_fits, number_parameters] and data type np.float32
:param constraints: Constraint bounds intervals - NumPy array of dimension [number_fits, 2*number_parameters] and data type np.float32
:param constraint_types: Types of constraints for all parameters (including fixed parameters) - NumPy array of length number_parameters and type np.int32 or None (means no constraints) with values from class ConstraintType
:param tolerance: The fit tolerance or None (will use default value)
:param max_number_iterations: The maximal number of iterations or None (will use default value)
:param parameters_to_fit: Which parameters to fit - NumPy array of length number_parameters and type np.int32 or None (will fit all parameters)
:param estimator_id: The Estimator ID or None (will use default values)
:param user_info: User info - NumPy array of type np.char or None (no user info available)
:return: parameters, states, chi_squares, number_iterations, execution_time
"""
# check all 2D NumPy arrays for row-major memory layout (otherwise interpretation of order of dimensions fails)
if not data.flags.c_contiguous:
raise RuntimeError('Memory layout of data array mismatch.')

if weights is not None and not weights.flags.c_contiguous:
raise RuntimeError('Memory layout of weights array mismatch.')

if not initial_parameters.flags.c_contiguous:
raise RuntimeError('Memory layout of initial_parameters array mismatch.')

# size check: data is 2D and read number of points and fits
if data.ndim != 2:
raise RuntimeError('data is not two-dimensional')
number_points = data.shape[1]
number_fits = data.shape[0]

# size check: consistency with weights (if given)
if weights is not None and data.shape != weights.shape:
raise RuntimeError('dimension mismatch between data and weights')
# the unequal operator checks, type, length and content (https://docs.python.org/3.7/reference/expressions.html#value-comparisons)

# size check: initial parameters is 2D and read number of parameters
if initial_parameters.ndim != 2:
raise RuntimeError('initial_parameters is not two-dimensional')
number_parameters = initial_parameters.shape[1]
if initial_parameters.shape[0] != number_fits:
raise RuntimeError('dimension mismatch in number of fits between data and initial_parameters')

# size check: constraints is 2D and number of fits, 2x number of parameters if given
if constraints is not None:
if constraints.ndim != 2:
raise RuntimeError('constraints not two-dimensional')
if constraints.shape != (number_fits, 2*number_parameters):
raise RuntimeError('constraints array has invalid shape')

# size check: constraint_types has certain length (if given)
if constraint_types is not None and constraint_types.shape[0] != number_parameters:
raise RuntimeError('constraint_types should have length of number of parameters')

# size check: consistency with parameters_to_fit (if given)
if parameters_to_fit is not None and parameters_to_fit.shape[0] != number_parameters:
raise RuntimeError(
'dimension mismatch in number of parameters between initial_parameters and parameters_to_fit')

# default value constraint types
if constraint_types is None:
constraint_types = cp.full(number_parameters, ConstraintType.FREE, dtype=cp.int32)

# default value: tolerance
if tolerance is None:
tolerance = 1e-4

# default value: max_number_iterations
if max_number_iterations is None:
max_number_iterations = 25

# default value: estimator ID
if estimator_id is None:
estimator_id = EstimatorID.LSE

# default value: parameters_to_fit
if parameters_to_fit is None:
parameters_to_fit = np.ones(number_parameters, dtype=cp.int32)
# now only weights and user_info could be not given

# type check: data, weights (if given), initial_parameters, constraints (if given) are all cp.float32
if data.dtype != cp.float32:
raise RuntimeError('type of data is not cp.float32')
if weights is not None and weights.dtype != cp.float32:
raise RuntimeError('type of weights is not cp.float32')
if initial_parameters.dtype != cp.float32:
raise RuntimeError('type of initial_parameters is not cp.float32')
if constraints is not None and constraints.dtype != cp.float32:
raise RuntimeError('type of constraints is not cp.float32')

# type check: parameters_to_fit, constraint_types is cp.int32
if parameters_to_fit.dtype != cp.int32:
raise RuntimeError('type of parameters_to_fit is not cp.int32')
if constraint_types.dtype != cp.int32:
raise RuntimeError('type of constraint_types is not cp.int32')

# type check: valid model, estimator id, constraint_types
if not _valid_id(ModelID, model_id):
raise RuntimeError('Invalid model ID, use an attribute of ModelID')
if not _valid_id(EstimatorID, estimator_id):
raise RuntimeError('Invalid estimator ID, use an attribute of EstimatorID')
if not all(_valid_id(ConstraintType, constraint_type) for constraint_type in constraint_types):
raise RuntimeError('Invalid constraint type, use an attribute of ConstraintType')

# we don't check type of user_info, but we extract the size in bytes of it
if user_info is not None:
user_info_size = user_info.nbytes
else:
user_info_size = 0

# pre-allocate output variables
# parameters = cp.zeros((number_fits, number_parameters), dtype=cp.float32)
states = cp.zeros(number_fits, dtype=cp.int32)
chi_squares = cp.zeros(number_fits, dtype=cp.float32)
number_iterations = cp.zeros(number_fits, dtype=cp.int32)

# conversion to ctypes types for optional C interface parameters using NULL pointer (None) as default argument
if weights is not None:
weights_p = cast(weights.data.ptr, gpufit_func_cuda.argtypes[3])
else:
weights_p = None
if constraints is not None:
constraints_p = cast(constraints.data.ptr, gpufit_func_cuda.argtypes[8])
else:
constraints_p = None
if user_info is not None:
user_info_p = cast(user_info.data.ptr, gpufit_func_cuda.argtypes[12])
else:
user_info_p = None

# call into the library (measure time)
t0 = time.perf_counter()
status = gpufit_func_cuda(
gpufit_func_cuda.argtypes[0](number_fits),
gpufit_func_cuda.argtypes[1](number_points),
cast(data.data.ptr, gpufit_func_cuda.argtypes[2]),
weights_p,
gpufit_func_cuda.argtypes[4](model_id),
gpufit_func_cuda.argtypes[5](tolerance),
gpufit_func_cuda.argtypes[6](max_number_iterations),
parameters_to_fit.ctypes.data_as(gpufit_func_cuda.argtypes[7]),
constraints_p,
cast(constraint_types.data.ptr, gpufit_func_cuda.argtypes[9]),
gpufit_func_cuda.argtypes[10](estimator_id),
gpufit_func_cuda.argtypes[11](user_info_size),
user_info_p,
cast(initial_parameters.data.ptr, gpufit_func_cuda.argtypes[13]),
cast(states.data.ptr, gpufit_func_cuda.argtypes[14]),
cast(chi_squares.data.ptr, gpufit_func_cuda.argtypes[15]),
cast(number_iterations.data.ptr, gpufit_func_cuda.argtypes[16]))
t1 = time.perf_counter()

# check status
if status != Status.Ok:
# get error from last error and raise runtime error
error_message = error_func()
raise RuntimeError('status = {}, message = {}'.format(status, error_message))

# return output values
return initial_parameters, states, chi_squares, number_iterations, t1 - t0

except ImportError as e:
raise e
119 changes: 119 additions & 0 deletions examples/python/gauss2d_cupy.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,119 @@
"""
Example of the Python binding of the Gpufit library which implements
Levenberg Marquardt curve fitting in CUDA
https://github.com/gpufit/Gpufit
http://gpufit.readthedocs.io/en/latest/bindings.html#python

Multiple fits of a 2D Gaussian peak function with Poisson distributed noise
This example additionally requires numpy.
"""

import cupy as np
import pygpufit.gpufit as gf


def generate_gauss_2d(p, xi, yi):
"""
Generates a 2D Gaussian peak.
http://gpufit.readthedocs.io/en/latest/api.html#gauss-2d

:param p: Parameters (amplitude, x,y center position, width, offset)
:param xi: x positions
:param yi: y positions
:return: The Gaussian 2D peak.
"""

arg = -(np.square(xi - p[1]) + np.square(yi - p[2])) / (2 * p[3] * p[3])
y = p[0] * np.exp(arg) + p[4]

return y


if __name__ == '__main__':

# cuda available checks
print('CUDA available: {}'.format(gf.cuda_available()))
if not gf.cuda_available():
raise RuntimeError(gf.get_last_error())
print('CUDA versions runtime: {}, driver: {}'.format(*gf.get_cuda_version()))

# number of fits and fit points
number_fits = 10000
size_x = 12
number_points = size_x * size_x
number_parameters = 5

# set input arguments

# true parameters
true_parameters = np.array((10, 5.5, 5.5, 3, 10), dtype=np.float32)

# initialize random number generator
np.random.seed(0)

# initial parameters (relative randomized, positions relative to width)
initial_parameters = np.tile(true_parameters, (number_fits, 1))
initial_parameters[:, (1, 2)] += true_parameters[3] * (-0.2 + 0.4 * np.random.rand(number_fits, 2))
initial_parameters[:, (0, 3, 4)] *= 0.8 + 0.4 * np.random.rand(number_fits, 3)

# generate x and y values
g = np.arange(size_x)
yi, xi = np.meshgrid(g, g, indexing='ij')
xi = xi.astype(np.float32)
yi = yi.astype(np.float32)

# generate data
data = generate_gauss_2d(true_parameters, xi, yi)
data = np.reshape(data, (1, number_points))
data = np.tile(data, (number_fits, 1))

# add Poisson noise
data = np.random.poisson(data)
data = data.astype(np.float32)

# tolerance
tolerance = 0.0001

# maximum number of iterations
max_number_iterations = 20

# estimator ID
estimator_id = gf.EstimatorID.MLE

# model ID
model_id = gf.ModelID.GAUSS_2D

# run Gpufit
parameters, states, chi_squares, number_iterations, execution_time = gf.fit_gpu(data, None, model_id,
initial_parameters,
tolerance,
max_number_iterations, None,
estimator_id, None)

# print fit results
converged = states == 0
print('*Gpufit*')

# print summary
print('\nmodel ID: {}'.format(model_id))
print('number of fits: {}'.format(number_fits))
print('fit size: {} x {}'.format(size_x, size_x))
print('mean chi_square: {:.2f}'.format(np.mean(chi_squares[converged])))
print('iterations: {:.2f}'.format(np.mean(number_iterations[converged])))
print('time: {:.2f} s'.format(execution_time))

# get fit states
number_converged = np.sum(converged)
print('\nratio converged {:6.2f} %'.format(number_converged / number_fits * 100))
print('ratio max it. exceeded {:6.2f} %'.format(np.sum(states == 1) / number_fits * 100))
print('ratio singular hessian {:6.2f} %'.format(np.sum(states == 2) / number_fits * 100))
print('ratio neg curvature MLE {:6.2f} %'.format(np.sum(states == 3) / number_fits * 100))

# mean, std of fitted parameters
converged_parameters = parameters[converged, :]
converged_parameters_mean = np.mean(converged_parameters, axis=0)
converged_parameters_std = np.std(converged_parameters, axis=0)
print('\nparameters of 2D Gaussian peak')
for i in range(number_parameters):
print('p{} true {:6.2f} mean {:6.2f} std {:6.2f}'.format(i, true_parameters[i], converged_parameters_mean[i],
converged_parameters_std[i]))
2 changes: 1 addition & 1 deletion examples/python/gauss2d_plot.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ def gaussians_2d(x, y, p):
model_id = gf.ModelID.GAUSS_2D

# loop over different number of fits
n_fits_all = np.around(np.logspace(2, 6, 20)).astype(np.int)
n_fits_all = np.around(np.logspace(2, 6, 20)).astype(int)

# generate x and y values
g = np.arange(size_x)
Expand Down