diff --git a/CMakeLists.txt b/CMakeLists.txt index e6278cd..1543461 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -182,6 +182,13 @@ if( PYTHONINTERP_FOUND ) add_python_launcher( RUN_PYTHON "${CMAKE_SOURCE_DIR}/Gpufit/python" "${Python_WORKING_DIRECTORY}" + ) + add_custom_target( RUN_PYTHON_CPUFIT ) + set_property( TARGET RUN_PYTHON_CPUFIT PROPERTY FOLDER CMakePredefinedTargets ) + add_dependencies( RUN_PYTHON_CPUFIT Cpufit ) + add_python_launcher( RUN_PYTHON_CPUFIT + "${CMAKE_SOURCE_DIR}/Cpufit/python" + "${Python_WORKING_DIRECTORY}" ) endif() endif() diff --git a/Cpufit/CMakeLists.txt b/Cpufit/CMakeLists.txt index cd4a7d6..af2038d 100644 --- a/Cpufit/CMakeLists.txt +++ b/Cpufit/CMakeLists.txt @@ -31,3 +31,4 @@ set_target_properties( Cpufit #install( TARGETS Cpufit RUNTIME DESTINATION bin ) add_subdirectory( matlab ) +add_subdirectory( python ) diff --git a/Cpufit/cpufit.h b/Cpufit/cpufit.h index db92aa7..e0156c2 100644 --- a/Cpufit/cpufit.h +++ b/Cpufit/cpufit.h @@ -9,6 +9,10 @@ #define VISIBLE #endif +#ifdef __APPLE__ +#define VISIBLE __attribute__((visibility("default"))) +#endif + #include #include "../Gpufit/constants.h" #include "../Gpufit/definitions.h" diff --git a/Cpufit/python/CMakeLists.txt b/Cpufit/python/CMakeLists.txt new file mode 100755 index 0000000..b65b27a --- /dev/null +++ b/Cpufit/python/CMakeLists.txt @@ -0,0 +1,58 @@ + +# Python + +# Python package + +set( build_directory "${CMAKE_BINARY_DIR}/${CMAKE_CFG_INTDIR}/pyCpufit" ) +set( setup_files + "${CMAKE_CURRENT_SOURCE_DIR}/README.txt" + "${CMAKE_CURRENT_SOURCE_DIR}/setup.py" + "${CMAKE_CURRENT_SOURCE_DIR}/setup.cfg" +) +set( module_directory "${build_directory}/pycpufit" ) +set( module_files + "${CMAKE_CURRENT_SOURCE_DIR}/pycpufit/__init__.py" + "${CMAKE_CURRENT_SOURCE_DIR}/pycpufit/cpufit.py" + "${CMAKE_CURRENT_SOURCE_DIR}/pycpufit/version.py" +) + + +set( binary $ ) + + +add_custom_target( PYTHON_PACKAGE_CPUFIT + COMMAND ${CMAKE_COMMAND} -E + remove_directory ${build_directory} + COMMAND ${CMAKE_COMMAND} -E + make_directory ${build_directory} + COMMAND ${CMAKE_COMMAND} -E + copy_if_different ${setup_files} ${build_directory} + COMMAND ${CMAKE_COMMAND} -E + make_directory ${module_directory} + COMMAND ${CMAKE_COMMAND} -E + copy_if_different ${module_files} ${module_directory} + COMMAND ${CMAKE_COMMAND} -E + copy_if_different ${binary} ${module_directory} +) + +set_property( TARGET PYTHON_PACKAGE_CPUFIT PROPERTY FOLDER CMakePredefinedTargets ) +add_dependencies( PYTHON_PACKAGE_CPUFIT Cpufit ) + +if( NOT PYTHONINTERP_FOUND ) + message( STATUS "Python NOT found - skipping creation of Python wheel!" ) + return() +endif() + +# Python wheel (output name is incorrect, requires plattform tag, see packaging) + +add_custom_target( PYTHON_WHEEL_CPUFIT ALL + COMMAND ${CMAKE_COMMAND} -E + chdir ${build_directory} "${PYTHON_EXECUTABLE}" setup.py clean --all + COMMAND ${CMAKE_COMMAND} -E + chdir ${build_directory} "${PYTHON_EXECUTABLE}" setup.py bdist_wheel + COMMENT "Preparing Python Wheel" +) +set_property( TARGET PYTHON_WHEEL_CPUFIT PROPERTY FOLDER CMakePredefinedTargets ) +add_dependencies( PYTHON_WHEEL_CPUFIT PYTHON_PACKAGE_CPUFIT ) + +# add launcher to Python package diff --git a/Cpufit/python/README.txt b/Cpufit/python/README.txt new file mode 100755 index 0000000..de3db9a --- /dev/null +++ b/Cpufit/python/README.txt @@ -0,0 +1,22 @@ +Python binding for the [Cpufit library](https://github.com/gpufit/Gpufit) which implements Levenberg Marquardt curve fitting in CUDA + +Requirements + +- Windows +- Python 2 or 3 with NumPy + +Installation + +Currently the wheel file has to be installed locally. + +If NumPy is not yet installed, install it using pip from the command line + +pip install numpy + +Then install pyCpufit from the local folder via: + +pip install --no-index --find-links=LocalPathToWheelFile pyCpufit + +Examples + +See project-root/examples/python folder. **Cpufit** examples begin with `cpufit_`. diff --git a/Cpufit/python/pycpufit/__init__.py b/Cpufit/python/pycpufit/__init__.py new file mode 100755 index 0000000..eb008c8 --- /dev/null +++ b/Cpufit/python/pycpufit/__init__.py @@ -0,0 +1 @@ +from cpufit import * diff --git a/Cpufit/python/pycpufit/cpufit.py b/Cpufit/python/pycpufit/cpufit.py new file mode 100755 index 0000000..8e79bc0 --- /dev/null +++ b/Cpufit/python/pycpufit/cpufit.py @@ -0,0 +1,297 @@ +""" +Program written from scratch following example of pyGPUfit by Jakub Dokulil (jakub.dokulil@imp.ac.at). + +For GPUfit see https://github.com/gpufit/Gpufit, http://gpufit.readthedocs.io/en/latest/bindings.html#python + +""" +# %% + +import os +import time +from ctypes import cdll, POINTER, byref, c_int, c_float, c_char, c_char_p, c_size_t +import numpy as np + +# %% +package_dir = os.path.dirname(os.path.realpath(__file__)) + +if os.name == 'nt': + lib_path = os.path.join(package_dir, 'cpufit.dll') # library name on Windows +elif os.name == 'posix': + lib_path = os.path.join(package_dir, 'libCpufit.so') # library name on Unix +else: + raise RuntimeError('OS {} not supported by pyGpufit.'.format(os.name)) + +lib = cdll.LoadLibrary(lib_path) + +# %% +cpufit_func = lib.cpufit + +# %% + +cpufit_func.restype = c_int +cpufit_func.argtypes = [c_size_t, c_size_t, POINTER(c_float), POINTER(c_float), c_int, POINTER(c_float), + 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)] + +# int cpufit +# ( +# std::size_t n_fits, - 1 +# std::size_t n_points, - number of points to be fited +# REAL * data, - flattened image with data +# REAL * weights, - wheights (see alesandros example) +# int model_id, - we want 2 GAUSS_2D_ELLIPTIC +# REAL * initial_parameters, - [p0, p1, ..., p5] +# REAL tolerance, - lets set 1e-3 +# int max_n_iterations, 1200 +# int * parameters_to_fit, - [1, 1, 1, 1, 1, 1] +# int estimator_id, - 0 +# std::size_t user_info_size, - 0 +# char * user_info, - NULL +# REAL * output_parameters, +# int * output_states, +# REAL * output_chi_squares, +# int * output_n_iterations +# ) + +error_func = lib.cpufit_get_last_error +error_func.restype = c_char_p +error_func.argtypes = None + +class ModelID: + GAUSS_1D = 0 + GAUSS_2D = 1 + GAUSS_2D_ELLIPTIC = 2 + GAUSS_2D_ROTATED = 3 + CAUCHY_2D_ELLIPTIC = 4 + LINEAR_1D = 5 + FLETCHER_POWELL = 6 + BROWN_DENNIS = 7 + SPLINE_1D = 8 + SPLINE_2D = 9 + SPLINE_3D = 10 + SPLINE_3D_MULTICHANNEL = 11 + SPLINE_3D_PHASE_MULTICHANNEL = 12 + +class EstimatorID: + LSE = 0 + MLE = 1 + +class ConstraintType: + FREE = 0 + LOWER = 1 + UPPER = 2 + LOWER_UPPER = 3 + +class Status: + Ok = 0 + Error = 1 + +def _valid_id(cls, id): + properties = [key for key in cls.__dict__.keys() if not key.startswith('__')] + values = [cls.__dict__[key] for key in properties] + return id in values + +def fit(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 fit function in the library. + (see also ) + + 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(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(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 fit 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 = np.full(number_parameters, ConstraintType.FREE, dtype=np.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=np.int32) + + # now only weights and user_info could be not given + + # type check: data, weights (if given), initial_parameters, constraints (if given) are all np.float32 + if data.dtype != np.float32: + raise RuntimeError('type of data is not np.float32') + if weights is not None and weights.dtype != np.float32: + raise RuntimeError('type of weights is not np.float32') + if initial_parameters.dtype != np.float32: + raise RuntimeError('type of initial_parameters is not np.float32') + if constraints is not None and constraints.dtype != np.float32: + raise RuntimeError('type of constraints is not np.float32') + + # type check: parameters_to_fit, constraint_types is np.int32 + if parameters_to_fit.dtype != np.int32: + raise RuntimeError('type of parameters_to_fit is not np.int32') + if constraint_types.dtype != np.int32: + raise RuntimeError('type of constraint_types is not np.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 = np.zeros((number_fits, number_parameters), dtype=np.float32) + states = np.zeros(number_fits, dtype=np.int32) + chi_squares = np.zeros(number_fits, dtype=np.float32) + number_iterations = np.zeros(number_fits, dtype=np.int32) + + # conversion to ctypes types for optional C interface parameters using NULL pointer (None) as default argument + if weights is not None: + weights_p = weights.ctypes.data_as(cpufit_func.argtypes[3]) + else: + weights_p = None + if constraints is not None: + constraints_p = constraints.ctypes.data_as(cpufit_func.argtypes[6]) + else: + constraints_p = None + if user_info is not None: + user_info_p = user_info.ctypes.data_as(cpufit_func.argtypes[13]) + else: + user_info_p = None + + # call into the library (measure time) + t0 = time.perf_counter() + status = cpufit_func( + cpufit_func.argtypes[0](number_fits), \ + cpufit_func.argtypes[1](number_points), \ + data.ctypes.data_as(cpufit_func.argtypes[2]), \ + weights_p, \ + cpufit_func.argtypes[4](model_id), \ + initial_parameters.ctypes.data_as(cpufit_func.argtypes[5]), \ + # constraints_p, \ + # constraint_types.ctypes.data_as(cpufit_func.argtypes[7]), \ + cpufit_func.argtypes[6](tolerance), \ + cpufit_func.argtypes[7](max_number_iterations), \ + parameters_to_fit.ctypes.data_as(cpufit_func.argtypes[8]), \ + cpufit_func.argtypes[9](estimator_id), \ + cpufit_func.argtypes[10](user_info_size), \ + user_info_p, \ + parameters.ctypes.data_as(cpufit_func.argtypes[12]), \ + states.ctypes.data_as(cpufit_func.argtypes[13]), \ + chi_squares.ctypes.data_as(cpufit_func.argtypes[14]), \ + number_iterations.ctypes.data_as(cpufit_func.argtypes[15])) + 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 parameters, states, chi_squares, number_iterations, t1 - t0 + +def get_last_error(): + """ + :return: Error message of last error. + """ + return error_func() diff --git a/Cpufit/python/pycpufit/version.py b/Cpufit/python/pycpufit/version.py new file mode 100644 index 0000000..2032e84 --- /dev/null +++ b/Cpufit/python/pycpufit/version.py @@ -0,0 +1,8 @@ +""" +Current version (short and full names). +""" + +__version_short__ = '1.2' +__version_full__ = '1.2.0' + +__version__ = __version_full__ \ No newline at end of file diff --git a/Cpufit/python/requirements.txt b/Cpufit/python/requirements.txt new file mode 100755 index 0000000..b316bf2 --- /dev/null +++ b/Cpufit/python/requirements.txt @@ -0,0 +1 @@ +NumPy>=1.8 \ No newline at end of file diff --git a/Cpufit/python/setup.cfg b/Cpufit/python/setup.cfg new file mode 100755 index 0000000..5eab412 --- /dev/null +++ b/Cpufit/python/setup.cfg @@ -0,0 +1,2 @@ +[bdist_wheel] +universal=1 diff --git a/Cpufit/python/setup.py b/Cpufit/python/setup.py new file mode 100755 index 0000000..6896398 --- /dev/null +++ b/Cpufit/python/setup.py @@ -0,0 +1,48 @@ +""" + setup script for pyCpufit + + TODO get version, get meaningful email +""" + +from setuptools import setup, find_packages +import os +from io import open # to have encoding as parameter of open on Python >=2.6 +import pycpufit.version as vs + +if os.name == 'nt': + lib_ext = '.dll' # library name extension on Windows +elif os.name == 'posix': + lib_ext = '.so' # library name extensions on Unix +else: + raise RuntimeError('OS {} not supported'.format(os.name)) + +HERE = os.path.abspath(os.path.dirname(__file__)) + +CLASSIFIERS = ['Development Status :: 5 - Production/Stable', + 'Intended Audience :: End Users/Desktop', + 'Operating System :: Microsoft :: Windows', + 'Topic :: Scientific/Engineering', + 'Topic :: Software Development :: Libraries'] + +def get_long_description(): + """ + Get the long description from the README file. + """ + with open(os.path.join(HERE, 'README.txt'), encoding='utf-8') as f: + return f.read() + +if __name__ == "__main__": + setup(name='pyCpufit', + version=vs.__version__, + description='Levenberg Marquardt curve fitting on the CPU', + long_description=get_long_description(), + url='https://github.com/gpufit/Gpufit', + author='M. Bates, A. Przybylski, B. Thiel, J. Dokulil, and J. Keller-Findeisen', + author_email='a@b.c', + license='MIT license', + classifiers=[], + keywords='Levenberg Marquardt, curve fitting', + packages=find_packages(where=HERE), + package_data={'pycpufit': ['*{}'.format(lib_ext)]}, + install_requires=['NumPy>=1.8'], + zip_safe=False) diff --git a/Cpufit/python/tests/run_tests.py b/Cpufit/python/tests/run_tests.py new file mode 100755 index 0000000..a2ea8dd --- /dev/null +++ b/Cpufit/python/tests/run_tests.py @@ -0,0 +1,19 @@ +""" +Discovers all tests and runs them. Assumes that initially the working directory is test. +""" + +import sys +import unittest + +if __name__ == '__main__': + + loader = unittest.defaultTestLoader + + tests = loader.discover('.') + + runner = unittest.TextTestRunner() + + results = runner.run(tests) + + # return number of failures + sys.exit(len(results.failures)) \ No newline at end of file diff --git a/Cpufit/python/tests/test_gaussian_fit_1d.py b/Cpufit/python/tests/test_gaussian_fit_1d.py new file mode 100755 index 0000000..e92b9de --- /dev/null +++ b/Cpufit/python/tests/test_gaussian_fit_1d.py @@ -0,0 +1,77 @@ +""" + Equivalent to https://github.com/gpufit/Gpufit/blob/master/Gpufit/tests/Gauss_Fit_1D.cpp +""" + +import unittest +import numpy as np +import pycpufit.cpufit as cf + +def generate_gauss_1d(parameters, x): + """ + Generates a 1D Gaussian curve. + + :param parameters: The parameters (a, x0, s, b) + :param x: The x values + :return: A 1D Gaussian curve. + """ + + a = parameters[0] + x0 = parameters[1] + s = parameters[2] + b = parameters[3] + + y = a * np.exp(-np.square(x - x0) / (2 * s**2)) + b + + return y + +class Test(unittest.TestCase): + + def test_gaussian_fit_1d(self): + # constants + n_fits = 1 + n_points = 5 + n_parameter = 4 # model will be GAUSS_1D + + # true parameters + true_parameters = np.array((4, 2, 0.5, 1), dtype=np.float32) + + # generate data + data = np.empty((n_fits, n_points), dtype=np.float32) + x = np.arange(n_points, dtype=np.float32) + data[0, :] = generate_gauss_1d(true_parameters, x) + + # tolerance + tolerance = 0.001 + + # max_n_iterations + max_n_iterations = 10 + + # model id + model_id = cf.ModelID.GAUSS_1D + + # initial parameters + initial_parameters = np.empty((n_fits, n_parameter), dtype=np.float32) + initial_parameters[0, :] = (2, 1.5, 0.3, 0) + + # call to cpufit + parameters, states, chi_squares, number_iterations, execution_time = cf.fit(data, None, model_id, + initial_parameters, tolerance, \ + max_n_iterations, None, None, None) + + # print results + for i in range(n_parameter): + print(' p{} true {} fit {}'.format(i, true_parameters[i], parameters[0, i])) + print('fit state : {}'.format(states)) + print('chi square: {}'.format(chi_squares)) + print('iterations: {}'.format(number_iterations)) + print('time: {} s'.format(execution_time)) + + assert (chi_squares < 1e-6) + assert (states == 0) + assert (number_iterations <= max_n_iterations) + for i in range(n_parameter): + assert (abs(true_parameters[i] - parameters[0, i]) < 1e-6) + +if __name__ == '__main__': + + unittest.main() diff --git a/Cpufit/python/tests/test_linear_regression.py b/Cpufit/python/tests/test_linear_regression.py new file mode 100755 index 0000000..60b3527 --- /dev/null +++ b/Cpufit/python/tests/test_linear_regression.py @@ -0,0 +1,61 @@ +""" + Equivalent to https://github.com/gpufit/Gpufit/blob/master/Gpufit/tests/Linear_Fit_1D.cpp +""" + +import unittest +import numpy as np +import pycpufit.cpufit as cf + +class Test(unittest.TestCase): + + def test_gaussian_fit_1d(self): + # constants + n_fits = 1 + n_points = 2 + n_parameter = 2 + + # true parameters + true_parameters = np.array((0, 1), dtype=np.float32) + + # data values + data = np.empty((n_fits, n_points), dtype=np.float32) + data[0, :] = (0, 1) + + # max number iterations + max_number_iterations = 10 + + # initial parameters + initial_parameters = np.empty((n_fits, n_parameter), dtype=np.float32) + initial_parameters[0, :] = (0, 0) + + # model id + model_id = cf.ModelID.LINEAR_1D + + # tolerance + tolerance = 0.001 + + # user info + user_info = np.array((0, 1), dtype=np.float32) + + # call to cpufit + parameters, states, chi_squares, number_iterations, execution_time = cf.fit(data, None, model_id, + initial_parameters, tolerance, \ + None, None, None, user_info) + + # print results + for i in range(n_parameter): + print(' p{} true {} fit {}'.format(i, true_parameters[i], parameters[0, i])) + print('fit state : {}'.format(states)) + print('chi square: {}'.format(chi_squares)) + print('iterations: {}'.format(number_iterations)) + print('time: {} s'.format(execution_time)) + + assert (chi_squares < 1e-6) + assert (states == 0) + assert (number_iterations <= max_number_iterations) + for i in range(n_parameter): + assert (abs(true_parameters[i] - parameters[0, i]) < 1e-6) + +if __name__ == '__main__': + + unittest.main() diff --git a/docs/installation.rst b/docs/installation.rst index edbc47d..506a3a8 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -233,6 +233,8 @@ The unit tests can be executed by building the target "RUN_TESTS" or by starting the created executables in the output directory from the command line. +To install the python gpufit package go to "Release" folder and here you find Gpufit resp Cpufit python packages. + Compiling Gpufit on Linux ------------------------- diff --git a/examples/python/cpufit_gauss2d.py b/examples/python/cpufit_gauss2d.py new file mode 100644 index 0000000..ab0be48 --- /dev/null +++ b/examples/python/cpufit_gauss2d.py @@ -0,0 +1,113 @@ +""" +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 numpy as np +import cpufit as cf + +# %% +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__': + + # 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 = cf.EstimatorID.MLE + + # model ID + model_id = cf.ModelID.GAUSS_2D + + # run Gpufit + parameters, states, chi_squares, number_iterations, execution_time = cf.fit(data, None, model_id, + initial_parameters, + tolerance, max_number_iterations, None, + estimator_id, None) + + # print fit results + converged = states == 0 + print('*Cpufit*') + + # 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/cpufit_gauss2d_plot.py b/examples/python/cpufit_gauss2d_plot.py new file mode 100644 index 0000000..50bb40c --- /dev/null +++ b/examples/python/cpufit_gauss2d_plot.py @@ -0,0 +1,122 @@ +""" +Example of the Python binding of the Gpufit library which implements +Levenberg Marquardt curve fitting in CUDA +https://github.com/gpufit/Gpufit + +Multiple fits of a 2D Gaussian peak function with Poisson distributed noise +repeated for a different total number of fits each time and plotting the results +http://gpufit.readthedocs.io/en/latest/bindings.html#python + +This example additionally requires numpy (http://www.numpy.org/) and matplotlib (http://matplotlib.org/). +""" + +import numpy as np +import matplotlib.pyplot as plt +import pycpufit as cf + + +def gaussians_2d(x, y, p): + """ + Generates many 2D Gaussians peaks for a given set of parameters + """ + + n_fits = p.shape[0] + + y = np.zeros((n_fits, x.shape[0], x.shape[1]), dtype=np.float32) + + # loop over each fit + for i in range(n_fits): + pi = p[i, :] + arg = -(np.square(xi - pi[1]) + np.square(yi - pi[2])) / (2 * pi[3] * pi[3]) + y[i, :, :] = pi[0] * np.exp(arg) + pi[4] + + return y + + +if __name__ == '__main__': + + # # cuda available checks + # print('CUDA available: {}'.format(cf.cuda_available())) + # if not cf.cuda_available(): + # raise RuntimeError(cf.get_last_error()) + # print('CUDA versions runtime: {}, driver: {}'.format(*cf.get_cuda_version())) + + # number of fit points + size_x = 5 + number_points = size_x * size_x + + # set input arguments + + # true parameters + mean_true_parameters = np.array((100, 2, 2, 1, 10), dtype=np.float32) + + # average noise level + average_noise_level = 10 + + # initialize random number generator + np.random.seed(0) + + # tolerance + tolerance = 0.0001 + + # maximum number of iterations + max_number_iterations = 10 + + # model ID + model_id = cf.ModelID.GAUSS_2D + + # loop over different number of fits + n_fits_all = np.around(np.logspace(2, 6, 20)).astype('int') + + # 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) + + # loop + speed = np.zeros(n_fits_all.size) + for i in range(n_fits_all.size): + n_fits = n_fits_all[i] + + # vary positions of 2D Gaussian peaks slightly + test_parameters = np.tile(mean_true_parameters, (n_fits, 1)) + test_parameters[:, (1, 2)] += mean_true_parameters[3] * (-0.2 + 0.4 * np.random.rand(n_fits, 2)) + + # generate data + data = gaussians_2d(xi, yi, test_parameters) + data = np.reshape(data, (n_fits, number_points)) + + # add noise + data += np.random.normal(scale=average_noise_level, size=data.shape) + + # initial parameters (randomized relative (to width for position)) + initial_parameters = np.tile(mean_true_parameters, (n_fits, 1)) + initial_parameters[:, (1, 2)] += mean_true_parameters[3] * (-0.2 + 0.4 * np.random.rand(n_fits, 2)) + initial_parameters[:, (0, 3, 4)] *= 0.8 + 0.4 * np.random.rand(n_fits, 3) + + # run Gpufit + parameters, states, chi_squares, number_iterations, execution_time = cf.fit(data, None, model_id, + initial_parameters, tolerance, + max_number_iterations) + + # analyze result + converged = states == 0 + speed[i] = n_fits / execution_time + precision_x0 = np.std(parameters[converged, 1] - test_parameters[converged, 1], axis=0, dtype=np.float64) + + # display result + '{} fits '.format(n_fits) + print('{:7} fits iterations: {:6.2f} | time: {:6.3f} s | speed: {:8.0f} fits/s' \ + .format(n_fits, np.mean(number_iterations[converged]), execution_time, speed[i])) + +# plot +plt.semilogx(n_fits_all, speed, 'bo-') +plt.grid(True) +plt.xlabel('number of fits per function call') +plt.ylabel('fits per second') +plt.legend(['Gpufit'], loc='upper left') +ax = plt.gca() +ax.set_xlim(n_fits_all[0], n_fits_all[-1]) + +plt.show()