diff --git a/Cpufit/python/pycpufit/__init__.py b/Cpufit/python/pycpufit/__init__.py index eb008c8b..8b137891 100755 --- a/Cpufit/python/pycpufit/__init__.py +++ b/Cpufit/python/pycpufit/__init__.py @@ -1 +1 @@ -from cpufit import * + diff --git a/Gpufit/python/pygpufit/gpufit.py b/Gpufit/python/pygpufit/gpufit.py index ee018ac8..450a7b96 100644 --- a/Gpufit/python/pygpufit/gpufit.py +++ b/Gpufit/python/pygpufit/gpufit.py @@ -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) @@ -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 @@ -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 diff --git a/examples/python/gauss2d_cupy.py b/examples/python/gauss2d_cupy.py new file mode 100644 index 00000000..adad31b8 --- /dev/null +++ b/examples/python/gauss2d_cupy.py @@ -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])) diff --git a/examples/python/gauss2d_plot.py b/examples/python/gauss2d_plot.py index bcaa3ac8..c4201a3d 100644 --- a/examples/python/gauss2d_plot.py +++ b/examples/python/gauss2d_plot.py @@ -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)