Repository: flimlib/flimlib (language: python)

Issue:
Python bindings; CI and versioning
@facetorched's [work](https://github.com/facetorched/flimlib/tree/python) on Python bindings is going well, and we should be ready to create a PR very soon. So I wanted to update you (esp. @paulbarber and @ctrueden) on what we're cooking up and also start a discussion on some logistical issues (mostly CI).

The expectation is that one will be able to do something like:

```bash
pip install flimlib  # Doesn't work yet
```

```python
import flimlib
result = flimlib.GCI_marquardt_fitting_engine(dt, photoncounts, [0.0, A_guess, tau_guess],
        fitfunc=flimlib.GCI_multiexp_tau)
print(result.chisq)
```

The bindings (hand-written using ctypes) are thin (but well-type-checked) wrappers around the C functions, with some Python niceness, such as many parameters being optional with reasonable defaults, and not needing to explicitly specify array lengths. Also, all array inputs and outputs can be NumPy arrays.

Only a subset of the RLD, LMA, and phasor functions (namely, the `*_fitting_engine` functions and their associated utility functions) are wrapped so far; more can be added later.

The build uses Python's standard `setuptools`, and bypasses both Maven and CMake. While this technically duplicates the C build, it ends up being much simpler for many detailed reasons. Adding a new `.c` file to flimlib does not require updating the Python build files, as `src/main/c/**/*.c` are automatically built. There are also unit(-ish) tests intended to test the Python bindings.

So, here are two things I wanted to ask about, planning slightly ahead:

- Seeing that build requirements are quite different for Python vs Java (e.g., we need a binary for each version of Python 3.6-3.9, and the Linux builds need to happen within a specific Docker container), I intend to help @facetorched set up CI in a way that is largely orthogonal to the existing CI for Java and C. Would you be okay if we did this using GitHub Actions? (Mostly because I have much more [experience](https://github.com/micro-manager/pymmcore) with it compared to Travis and AppVeyor.) Of course I am willing to maintain any such build for the foreseeable future.

- Versioning. It is simplest (and easy) for the Python build to read the version in `pom.xml`; we would convert `2.1.2-SNAPSHOT` to `2.1.2.dev0` to follow Python convensions. The question is, are we okay with a single version number for C/Java/Python packages, or do we want some decoupling? Assuming this is intended to be a SemVer, would you be okay if we occasionally incremented the minor version just because the Python bindings, say, gained more wrappers? I like the simplicity of just one version number for all, but just checking. In theory, we could also add a 4th number that indicates the Python bindings version.
python interface?
Any plans to add a python interface to call flimlib functions from python? Or any suggestions on how to do this? Or is there an equivalent package to flimlib in python that you are aware of?


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 36681da..6de9838 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,3 +8,11 @@ src/main/c/doxygen.exe
 *.log
 /bin/
 /Debug/
+dist/
+venv/
+__pycache__/
+*.egg-info
+_flimlib.*.pyd
+_flimlib.*.so
+.tox/
+.vscode/
\ No newline at end of file
diff --git a/MANIFEST.in b/MANIFEST.in
new file mode 100644
index 0000000..e7815dc
--- /dev/null
+++ b/MANIFEST.in
@@ -0,0 +1,1 @@
+recursive-include src/main/c *.h *.def
\ No newline at end of file
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..b5a3c46
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,6 @@
+[build-system]
+requires = [
+    "setuptools>=42",
+    "wheel"
+]
+build-backend = "setuptools.build_meta"
\ No newline at end of file
diff --git a/setup.cfg b/setup.cfg
new file mode 100644
index 0000000..9c2a052
--- /dev/null
+++ b/setup.cfg
@@ -0,0 +1,23 @@
+[metadata]
+name = flimlib
+version = 0.0.1.dev0
+author = Paul Barber, Curtis Rueden, Kevin Tan
+description = A python wrapper for flimlib's C functions
+long_description = file: README.md
+long_description_content_type = text/markdown
+url = https://github.com/flimlib/flimlib
+project_urls =
+    Bug Tracker = https://github.com/flimlib/flimlib/issues
+classifiers =
+    Programming Language :: Python :: 3
+    License :: OSI Approved :: GNU General Public License v3 (GPLv3)
+    Operating System :: OS Independent
+
+[options]
+package_dir =
+    = src/main/python
+packages = find:
+python_requires = >=3.6
+
+[options.packages.find]
+where = src/main/python
\ No newline at end of file
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000..c7105a3
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,65 @@
+import glob
+import os.path
+import setuptools
+import sys
+
+from distutils.command.build_ext import build_ext
+from distutils.unixccompiler import UnixCCompiler
+
+# Workaround for setuptools/distutils limitations described in:
+# https://github.com/pypa/setuptools/issues/1192 and
+# https://github.com/pypa/setuptools/issues/1732
+# In our case, we need -std=c++14 (or c++11) for macOS clang++, but setuptools
+# uses 'clang $CFLAGS' for both C and C++ sources. Passing a C++ std to the C
+# compiler is an error, so we need a mechanism to pass a flag only when
+# compiling C++ sources. (Credit to https://stackoverflow.com/a/65865696 for
+# this idea.)
+
+unix_cxx_flags = '-std=c++14' # We could also read os.environ.get('CXXFLAGS')
+
+class UnixCCxxCompiler(UnixCCompiler):
+    def _compile(self, obj, src, ext, cc_args, extra_postargs, pp_opts):
+        is_cxx_source = os.path.splitext(src)[-1] in ('.cpp', '.cxx', '.cc')
+        my_cc_args = ([unix_cxx_flags] + cc_args if is_cxx_source else cc_args)
+        super()._compile(obj, src, ext, my_cc_args, extra_postargs, pp_opts)
+
+class build_ext_c_cxx(build_ext):
+    def build_extension(self, ext):
+        if self.compiler.compiler_type == 'unix':
+            compiler = UnixCCxxCompiler()
+            for attr, value in self.compiler.__dict__.items():
+                setattr(compiler, attr, value)
+            self.compiler = compiler
+        super().build_extension(ext)
+
+
+# normpath cleans up slashes on Windows
+flimlib_c_sources = [os.path.normpath(p) for p in
+                     glob.glob('src/main/c/**/*.c', recursive=True)]
+flimlib_cpp_sources = [os.path.normpath(p) for p in
+                       glob.glob('src/main/c/**/*.cpp', recursive=True)]
+
+flimlib_msvc_def = 'src/main/c/flimlib.def'
+
+module_source = 'src/main/python/flimlib/flimlib_dummy.c'
+
+c_sources = flimlib_c_sources + flimlib_cpp_sources + [module_source]
+
+# Technically we should make conditional on compiler, not OS, but we only
+# support MSVC for now.
+link_args = ['/DEF:' + flimlib_msvc_def] if sys.platform == 'win32' else []
+
+
+flimlib_ext = setuptools.Extension(
+    'flimlib._flimlib',
+    sources=c_sources,
+    extra_link_args=link_args,
+)
+
+setuptools.setup(
+    install_requires=["numpy>=1.12.0"],
+    ext_modules=[flimlib_ext],
+    cmdclass={
+        'build_ext': build_ext_c_cxx,
+    },
+)
diff --git a/src/main/python/flimlib/__init__.py b/src/main/python/flimlib/__init__.py
new file mode 100644
index 0000000..3f9ecbc
--- /dev/null
+++ b/src/main/python/flimlib/__init__.py
@@ -0,0 +1,3 @@
+__all__ = ['cfunc', '_flimlib_dummy']
+
+from .cfunc import *
\ No newline at end of file
diff --git a/src/main/python/flimlib/cfunc.py b/src/main/python/flimlib/cfunc.py
new file mode 100644
index 0000000..b67c4f6
--- /dev/null
+++ b/src/main/python/flimlib/cfunc.py
@@ -0,0 +1,492 @@
+import ctypes
+import glob
+import os
+import warnings
+from collections import namedtuple
+
+import numpy as np
+
+folder = os.path.dirname(os.path.abspath(__file__))
+dll_candidates = (glob.glob(os.path.join(folder, '_flimlib.*.pyd')) +
+                  glob.glob(os.path.join(folder, '_flimlib.*.so')))
+if len(dll_candidates) > 1:
+    raise RuntimeError("More than one flimlib extension found???")
+if not dll_candidates:
+    raise RuntimeError("flimlib extension missing")
+dll_path = dll_candidates[0]
+
+# 0x8 = LOAD_WITH_ALTERED_SEARCH_PATH, allowing absolute path loading
+_flimlib = ctypes.CDLL(dll_path, winmode=0x8)
+
+def _prep_sig(noise_type, sig, len):
+    """
+    Helper function to make sure noise_type is valid and sig is the correct type 
+    for the noise_type given
+    """
+    if not noise_type in _noise_types.keys():
+        raise ValueError(
+            "invalid noise type. The valid types are: ", _noise_types.keys)
+    elif noise_type == 'NOISE_GIVEN':
+        sig = np.asarray(sig, dtype=np.float32)
+        if sig.shape != (len,):
+            raise ValueError("incorrect shape of sig")
+        sig = np.ctypeslib.as_ctypes(sig)
+    elif noise_type == 'NOISE_CONST':
+        sig = float(np.asarray(sig))  # convert to float
+        sig = ctypes.c_float(sig)
+    elif sig is not None:
+        message = "Expected sig=None for noise type " + str(noise_type) + ". The given value of sig will be ignored"
+        warnings.warn(message)
+        return None
+    return sig
+
+def _prep_instr_ninstr(instr):
+    """
+    Helper function to handle instr input and generate a value for ninstr
+    """
+    ninstr = 0
+    if instr is not None:
+        instr = np.asarray(instr, dtype=np.float32)
+        if instr.ndim != 1:
+            raise ValueError("instr must be 1 dimensional")
+        ninstr = instr.shape[0] # presumably shorter than photon_count
+        instr = np.ctypeslib.as_ctypes(instr)
+    return instr, ninstr
+
+def _prep_common_fit_params(photon_count):
+    """
+    Helper function to generate common flimlib inputs given the shape of photon_count
+    """
+    photon_count = np.asarray(photon_count, dtype=np.float32)
+    if photon_count.ndim != 1:
+        raise ValueError("photon_count must be a 1 dimensional")
+    
+    fit_start = 0
+    fit_end = photon_count.shape[0] # should it be len or len-1?
+
+    photon_count = np.ctypeslib.as_ctypes(photon_count)
+    fitted = np.ctypeslib.as_ctypes(np.empty(fit_end, dtype=np.float32))
+    residuals = np.ctypeslib.as_ctypes(np.empty(fit_end, dtype=np.float32))
+
+    return photon_count, fit_start, fit_end, fitted, residuals
+
+_GCI_ecf_matrix = _flimlib.GCI_ecf_matrix
+_GCI_ecf_matrix.argtypes = [ctypes.c_int, ctypes.c_int]
+_GCI_ecf_matrix.restype = ctypes.POINTER(ctypes.POINTER(ctypes.c_float))
+
+_GCI_ecf_free_matrix = _flimlib.GCI_ecf_free_matrix
+_GCI_ecf_free_matrix.argtypes = [ctypes.POINTER(ctypes.POINTER(ctypes.c_float))]
+
+
+class _EcfMatrix:
+    """
+    A class used to represent an allocated matrix used by flimlib
+
+    ...
+
+    Attributes
+    ----------
+    matrix : LP_LP_c_float
+        a ctypes object representing a float** matrix
+    nrows : int
+        the number of rows
+    ncols : int
+        the number of columns
+
+    Methods
+    -------
+    asarray()
+        Returns a numpy.ndarray copy of the matrix
+    """
+    def __init__(self, nrows, ncols):
+        self.matrix = _GCI_ecf_matrix(nrows, ncols)
+        self.nrows = nrows
+        self.ncols = ncols
+
+    def __del__(self):
+        _GCI_ecf_free_matrix(self.matrix)
+
+    def asarray(self):
+        """Returns a numpy.ndarray copy of the matrix"""
+        # find the address of the start of the pointers at the beginning of the matrix
+        # the first pointer can be found by indexing
+        ptr_addr = ctypes.addressof(self.matrix[0])
+        # perform pointer arithmatic to skip the pointers at the beginning
+        # this gives us the address of the beginning of the data we are interested in
+        data_addr = ptr_addr + self.nrows*ctypes.sizeof(ctypes.c_void_p)
+        # apparently we must first cast to c_void_p to properly access the address that we calculated
+        data_void_ptr = ctypes.c_void_p(data_addr)
+        # cast the void pointer into the desired type (a normal 2d float array)
+        data_ptr = ctypes.cast(data_void_ptr, ctypes.POINTER(ctypes.c_float))
+
+        return np.ctypeslib.as_array(data_ptr, shape=(self.nrows, self.ncols)).copy()
+    
+    def fill(self):
+        """for debug purposes only"""
+        for r in range(self.nrows):
+            for c in range(self.ncols):
+                self.matrix[r][c] = r*self.ncols+c # enumerate
+
+# Noise types used by flimlib
+# 4 and 5 should raise an error for GCI_triple_integral_fitting_engine
+_noise_types = {'NOISE_CONST': 0, 'NOISE_GIVEN': 1, 'NOISE_POISSON_DATA': 2,
+                'NOISE_POISSON_FIT': 3, 'NOISE_GAUSSIAN_FIT': 4, 'NOISE_MLE': 5}
+
+# restrain types used by flimlib
+_restrain_types = {'ECF_RESTRAIN_DEFAULT': 0, 'ECF_RESTRAIN_USER': 1}
+
+_GCI_triple_integral_fitting_engine = _flimlib.GCI_triple_integral_fitting_engine
+TripleIntegralResult = namedtuple('TripleIntegralResult', 'tries Z A tau fitted residuals chisq')
+_GCI_triple_integral_fitting_engine.argtypes = [
+    ctypes.c_float,                 # float xincr
+    ctypes.POINTER(ctypes.c_float), # float y[]
+    ctypes.c_int,                   # int fit_start
+    ctypes.c_int,                   # int fit_end
+    ctypes.POINTER(ctypes.c_float), # float instr[]
+    ctypes.c_int,                   # int ninstr
+    ctypes.c_int,                   # noise_type noise
+    ctypes.POINTER(ctypes.c_float), # float sig[]
+    ctypes.POINTER(ctypes.c_float), # float *Z
+    ctypes.POINTER(ctypes.c_float), # float *A
+    ctypes.POINTER(ctypes.c_float), # float *tau
+    ctypes.POINTER(ctypes.c_float), # float *fitted
+    ctypes.POINTER(ctypes.c_float), # float *residuals
+    ctypes.POINTER(ctypes.c_float), # float *chisq
+    ctypes.c_float                  # float chisq_target
+    ]
+
+def GCI_triple_integral_fitting_engine(period, photon_count,
+                                       instr=None, noise_type='NOISE_POISSON_FIT', sig=None,
+                                       chisq_target=1.1):
+    """
+    Performs an exponential fit on the data using Rapid Lifetime Determination
+
+    Parameters
+    ----------
+    period : float
+        The time between samples in `photon_count`
+    photon_count : array_like
+        The data to be fit. the length of this array determines the length of the fit
+    instr : {None, array_like}, optional
+        The instrument response. If `instr` is None, no instrument response will be used
+        (default is None)
+    noise_type : str, optional
+        The noise type to use. Valid values are: 'NOISE_CONST', 'NOISE_GIVEN', 
+        'NOISE_POISSON_DATA', 'NOISE_POISSON_FIT' (default is 'NOISE_POISSON_FIT')
+    sig : {None, float, array_like}, optional
+        The noise. A float array the same length as `photon_count` if `noise_type` 
+        is 'NOISE_GIVEN', a float if `noise_type` is 'NOISE_CONST' and None otherwise 
+        (default is None)
+    chisq_target : float, optional
+        The target chi squared value for the fit (default is 1.1)
+
+    Returns
+    -------
+    TripleIntegralResult
+        A namedtuple containing values in order: tries, Z, A, tau, fitted, residuals, chisq
+
+    Raises
+    ------
+    ValueError
+        If noise type passed is not implemented yet
+    """
+    period = ctypes.c_float(period)
+
+    photon_count, fit_start, fit_end, fitted, residuals = _prep_common_fit_params(photon_count)
+
+    instr, ninstr = _prep_instr_ninstr(instr)
+
+    # this lack of implementation is unique to triple integral
+    if noise_type == 'NOISE_GAUSSIAN_FIT' or noise_type == 'NOISE_MLE':
+        raise ValueError(
+            "Noise types 'NOISE_GAUSSIAN_FIT' and 'NOISE_MLE' are currently unimplemented")
+
+    sig = _prep_sig(noise_type, sig, np.asarray(photon_count).shape[0])
+
+    chisq_target = ctypes.c_float(chisq_target)
+
+    Z, A, tau, chisq = ctypes.c_float(), ctypes.c_float(
+    ), ctypes.c_float(), ctypes.c_float()  # output values
+
+    tries = _GCI_triple_integral_fitting_engine(
+        period, photon_count, fit_start, fit_end, instr, ninstr, 
+        _noise_types[noise_type], sig, Z, A, tau, fitted, 
+        residuals, chisq, chisq_target)
+    
+    return TripleIntegralResult(tries, Z.value, A.value, tau.value, 
+        np.asarray(fitted), np.asarray(residuals), chisq.value)
+
+
+_GCI_Phasor = _flimlib.GCI_Phasor # C function
+PhasorResult = namedtuple('PhasorResult', 'error_code Z u v taup taum tau fitted residuals chisq')
+_GCI_Phasor.argtypes = [
+    ctypes.c_float,                 # float xincr
+    ctypes.POINTER(ctypes.c_float), # float y[]
+    ctypes.c_int,                   # int fit_start
+    ctypes.c_int,                   # int fit_end
+    ctypes.POINTER(ctypes.c_float), # float *Z
+    ctypes.POINTER(ctypes.c_float), # float *u
+    ctypes.POINTER(ctypes.c_float), # float *v
+    ctypes.POINTER(ctypes.c_float), # float *taup
+    ctypes.POINTER(ctypes.c_float), # float *taum
+    ctypes.POINTER(ctypes.c_float), # float *tau
+    ctypes.POINTER(ctypes.c_float), # float *fitted
+    ctypes.POINTER(ctypes.c_float), # float *residuals
+    ctypes.POINTER(ctypes.c_float)  # float *chisq
+    ]
+
+def GCI_Phasor(period, photon_count):
+    """
+    Performs an exponential fit on the data using Phasors
+
+    Parameters
+    ----------
+    period : float
+        The time between samples in `photon_count`
+    photon_count : array_like
+        The data to be fit. the length of this array determines the length of the fit
+
+    Returns
+    -------
+    TripleIntegralResult
+        A namedtuple containing values in order: error_code, Z, u, v, taup, taum, tau, fitted, residuals, chisq
+    """
+    period = ctypes.c_float(period)
+
+    photon_count, fit_start, fit_end, fitted, residuals = _prep_common_fit_params(photon_count)
+    
+    Z, u, v, taup, taum, tau, chisq = (ctypes.c_float(), 
+    ctypes.c_float(), ctypes.c_float(), ctypes.c_float(), ctypes.c_float(), 
+    ctypes.c_float(), ctypes.c_float())  # output values
+
+    error_code = _GCI_Phasor(period, photon_count, fit_start, fit_end, Z, u, v, taup, taum, tau, fitted, residuals, chisq)
+    return PhasorResult(error_code, Z.value, u.value, v.value, taup.value, taum.value, 
+        tau.value, np.asarray(fitted), np.asarray(residuals), chisq.value)
+
+
+class FitFunc:
+    """
+    A class used to represent a fitfunc used by flimlib
+
+    ...
+
+    Attributes
+    ----------
+    python_func : function
+        A python function passed by the user to wrap as a ctypes function pointer
+    c_func : ctypes._CFuncPtr
+        A ctypes function pointer for the fitfunc for flimlib to use
+    nparam_predicate : function
+        A python function that returns True if the correct number of parameters
+        has been passed for the given fitfunc
+
+    Methods
+    -------
+    get_c_func(nparam_in)
+        Returns ctypes function pointer if `nparam_in` is valid
+    """
+    def __init__(self, python_func, nparam_predicate=None):
+        """
+        Parameters
+        ----------
+        python_func : {function, ctypes._CFuncPtr)
+            A python function passed by the user to wrap as a ctypes function pointer. 
+            If a ctypes function pointer is passed, it will be used directly instead
+        nparam_predicate : {None, function}, optional
+            A python function that returns True if the correct number of parameters
+            has been passed for the given fitfunc (default is None)
+        """
+        if isinstance(python_func, ctypes._CFuncPtr):
+            # a ctypes c function was passed instead of a python function. We just use it directly
+            self.c_func = python_func
+        else:
+            self.python_func = python_func
+            CFUNC = ctypes.CFUNCTYPE(None, ctypes.c_float, ctypes.POINTER(ctypes.c_float), 
+                ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_float), ctypes.c_int)
+            self.c_func = CFUNC(self._wrapped_python_func)
+        self.nparam_predicate=nparam_predicate
+
+    def _wrapped_python_func(self, x, param, y, dy_dparam, nparam):
+        param_in = np.ctypeslib.as_array(param, shape=(nparam,))
+        y_out, dy_dparam_out = self.python_func(x, param_in)
+        y.contents.value = ctypes.c_float(y_out).value
+        dy_dparam = np.ctypeslib.as_array(dy_dparam,shape=(nparam,))
+        dy_dparam[:] = np.asarray(dy_dparam_out, dtype=np.float32)
+
+    def get_c_func(self, nparam_in):
+        if self.nparam_predicate is None or self.nparam_predicate(nparam_in):
+            return self.c_func
+        else:
+            raise TypeError("Incorrect size of param")
+
+
+def _multiexp_predicate(n_param):
+    # 3 or greater and odd
+    return n_param >= 3 and n_param % 2 == 1
+
+def _stretchedexp_predicate(n_param):
+    return n_param == 4
+
+GCI_multiexp_tau = FitFunc(_flimlib.GCI_multiexp_tau, nparam_predicate=_multiexp_predicate)
+GCI_multiexp_lambda = FitFunc(_flimlib.GCI_multiexp_lambda, nparam_predicate=_multiexp_predicate)
+GCI_stretchedexp = FitFunc(_flimlib.GCI_stretchedexp, nparam_predicate=_stretchedexp_predicate)
+
+_GCI_set_restrain_limits = _flimlib.GCI_set_restrain_limits
+_GCI_set_restrain_limits.argtypes= [
+    ctypes.POINTER(ctypes.c_int),   # int restrain[]
+    ctypes.c_int,                   # int nparam
+    ctypes.POINTER(ctypes.c_float), # float minval[]
+    ctypes.POINTER(ctypes.c_float)  # float maxval[]
+]
+
+def GCI_set_restrain_limits(restrain, minval, maxval):
+    """
+    Sets global constraints on fit parameters
+
+    Parameters
+    ----------
+    restrain : array_like
+        Array of 1s and 0s. 1 indicates that the parameter will be constrained, 
+        0 indicates it will not be constrained
+    minval : array_like
+        The minimum acceptable value of each parameter
+    maxval : array_like
+        The maximum acceptable value of each parameter
+
+    Returns
+    -------
+    None
+    """
+    minval = np.asarray(minval,dtype=np.float32)
+    maxval = np.asarray(maxval,dtype=np.float32)
+    restrain = np.asarray(restrain,dtype=np.intc)
+
+    if not(restrain.ndim == minval.ndim == maxval.ndim == 1):
+        raise ValueError("restrain, minval and maxval must be 1 dimentional!")
+    if not(restrain.shape == minval.shape == maxval.shape):
+        raise ValueError("restrain, minval and maxval must have the same shape!")
+
+    nparam = restrain.shape[0]
+    minval = np.ctypeslib.as_ctypes(minval)
+    maxval = np.ctypeslib.as_ctypes(maxval)
+    restrain = np.ctypeslib.as_ctypes(restrain)
+
+    error_code = _GCI_set_restrain_limits(restrain, nparam, minval, maxval)
+
+    if(error_code == -1):
+        raise ValueError("invalid size of restrain")
+    if(error_code == -2):
+        raise ValueError('maxval must be element-wise greater than minval')
+
+_GCI_marquardt_fitting_engine = _flimlib.GCI_marquardt_fitting_engine # C function
+MarquardtResult = namedtuple('MarquardtResult', 'error_code param fitted residuals chisq covar alpha erraxes')
+_GCI_marquardt_fitting_engine.argtypes= [
+    ctypes.c_float,                                 # float xincr
+    ctypes.POINTER(ctypes.c_float),                 # float *trans
+    ctypes.c_int,                                   # int ndata
+    ctypes.c_int,                                   # int fit_start
+    ctypes.c_int,                                   # int fit_end
+    ctypes.POINTER(ctypes.c_float),                 # float instr[]
+    ctypes.c_int,                                   # int ninstr
+    ctypes.c_int,                                   # noise_type noise
+    ctypes.POINTER(ctypes.c_float),                 # float sig[]
+    ctypes.POINTER(ctypes.c_float),                 # float param[]
+    ctypes.POINTER(ctypes.c_int),                   # int paramfree[]
+    ctypes.c_int,                                   # int nparam 
+    ctypes.c_int,                                   # restrain_type restrain
+    ctypes.c_void_p,                                # void (*fitfunc)(float, float [], float *, float [], int)
+    ctypes.POINTER(ctypes.c_float),                 # float *fitted
+    ctypes.POINTER(ctypes.c_float),                 # float *residuals
+    ctypes.POINTER(ctypes.c_float),                 # float *chisq
+    ctypes.POINTER(ctypes.POINTER(ctypes.c_float)), # float **covar
+    ctypes.POINTER(ctypes.POINTER(ctypes.c_float)), # float **alpha
+    ctypes.POINTER(ctypes.POINTER(ctypes.c_float)), # float **erraxes
+    ctypes.c_float,                                 # float chisq_target
+    ctypes.c_float,                                 # float chisq_delta
+    ctypes.c_int,                                   # int chisq_percent
+    ]
+
+def GCI_marquardt_fitting_engine(period, photon_count, param, paramfree=None, restrain_type='ECF_RESTRAIN_DEFAULT',
+                                       fitfunc=GCI_multiexp_tau, instr=None, noise_type='NOISE_POISSON_FIT', sig=None,
+                                       chisq_target=1.1, chisq_delta=1E-5, chisq_percent=95):
+    """
+    Performs an exponential fit on the data using the Levenberg–Marquardt Algorithm
+
+    Parameters
+    ----------
+    period : float
+        The time between samples in `photon_count`
+    photon_count : array_like
+        The data to be fit. the length of this array determines the length of the fit
+    param : array_like
+        The initial parameters for the fit. The "guess"
+    paramfree : array_like, optional
+        Array of 1s and 0s. 1 indicates that the parameter is free, 0 indicates it is fixed.
+        If None is passed, all parameters will be free. (default is None)
+    restrain_type : str, optional
+        restrain type to use. Valid values are 'ECF_RESTRAIN_DEFAULT' and 'ECF_RESTRAIN_USER' (default is 'ECF_RESTRAIN_DEFAULT')
+    fitfunc : FitFunc, optional
+        A FitFunc object that contains the fit function to be used in the fit. (default is GCI_multiexp_tau)
+    instr : {None, array_like}, optional
+        The instrument response. If `instr` is None, no instrument response will be used
+        (default is None)
+    noise_type : str, optional
+        The noise type to use. Valid values are: 'NOISE_CONST', 'NOISE_GIVEN', 
+        'NOISE_POISSON_DATA', 'NOISE_POISSON_FIT' (default is 'NOISE_POISSON_FIT')
+    sig : {None, float, array_like}, optional
+        The noise. A float array the same length as `photon_count` if `noise_type` 
+        is 'NOISE_GIVEN', a float if `noise_type` is 'NOISE_CONST' and None otherwise 
+        (default is None)
+    chisq_target : float, optional
+        The target chi squared value for the fit (default is 1.1)
+    chisq_delta : float, optional
+        An individual fit will continue if the chisq value changes by more then this amount (default is 1E-5)
+    chisq_percent : int, optional
+        Defines the confidence interval when calculating the error axes (default is 95)
+
+    Returns
+    -------
+    MarquardtResult
+        A namedtuple containing values in order: error_code, param, fitted, residuals, chisq, covar, alpha, erraxes
+    """
+    period = ctypes.c_float(period)
+
+    photon_count, fit_start, fit_end, fitted, residuals = _prep_common_fit_params(photon_count)
+
+    ndata = np.asarray(photon_count).shape[0]
+
+    instr, ninstr = _prep_instr_ninstr(instr)
+
+    sig = _prep_sig(noise_type, sig, np.asarray(photon_count).shape[0])
+
+    param = np.asarray(param,dtype=np.float32)
+    nparam = param.shape[0]
+    param = np.ctypeslib.as_ctypes(param)
+
+    if paramfree is None:
+        paramfree = np.ones(nparam,dtype=np.intc) # default all parameters are free
+    else:
+        paramfree = np.asarray(paramfree,dtype=np.intc) # TODO should I check to make sure it's ones and zeros?
+    paramfree = np.ctypeslib.as_ctypes(paramfree)
+
+    chisq = ctypes.c_float()
+
+    chisq_target = ctypes.c_float(chisq_target)
+    chisq_delta = ctypes.c_float(chisq_delta)
+    chisq_percent = ctypes.c_int(chisq_percent)
+
+    # allocate the specific style of 2d array used by flimlib using EcfMatrix class
+    covar = _EcfMatrix(nparam, nparam)
+    alpha = _EcfMatrix(nparam, nparam)
+    erraxes = _EcfMatrix(nparam, nparam)
+
+    # cast ctypes function pointer to c_void_p
+    fitfunc = ctypes.cast(fitfunc.get_c_func(nparam), ctypes.c_void_p)
+
+    error_code = _GCI_marquardt_fitting_engine(period, photon_count, ndata, fit_start, fit_end, 
+        instr, ninstr, _noise_types[noise_type], sig, param, paramfree, nparam, _restrain_types[restrain_type],
+        fitfunc, fitted, residuals, chisq,
+        covar.matrix, alpha.matrix, erraxes.matrix, chisq_target, chisq_delta, chisq_percent)
+    return MarquardtResult(error_code, np.asarray(param), np.asarray(fitted), 
+        np.asarray(residuals), chisq.value, covar.asarray(), alpha.asarray(), erraxes.asarray())
+
diff --git a/src/main/python/flimlib/flimlib_dummy.c b/src/main/python/flimlib/flimlib_dummy.c
new file mode 100644
index 0000000..bcee81e
--- /dev/null
+++ b/src/main/python/flimlib/flimlib_dummy.c
@@ -0,0 +1,28 @@
+#define PY_SSIZE_T_CLEAN
+#include <Python.h>
+
+/*
+ * This is an empty Python module that allows us to use setuptools' build
+ * system to build flimlib as a Python extension, thereby automatically taking
+ * care of the platform and ABI tags in the DLL and wheel filenames.
+ * Actual access to fliblib only happens through ctypes.
+ */
+
+static PyModuleDef module_def = {
+    PyModuleDef_HEAD_INIT,
+    .m_name = "_flimlib",
+    .m_doc = "Dummy module for flimlib C library",
+    .m_size = -1,
+};
+
+PyMODINIT_FUNC
+PyInit__flimlib(void)
+{
+    PyObject *m;
+
+    m = PyModule_Create(&module_def);
+    if (m == NULL)
+        return NULL;
+
+    return m;
+}
\ No newline at end of file
diff --git a/tox.ini b/tox.ini
new file mode 100644
index 0000000..6a802f2
--- /dev/null
+++ b/tox.ini
@@ -0,0 +1,7 @@
+[tox]
+envlist = py3
+
+[testenv]
+deps = pytest
+changedir = src/test/python
+commands = pytest --basetemp="{envtmpdir}" {posargs}


Test patch:
diff --git a/src/test/python/test_flimlib.py b/src/test/python/test_flimlib.py
new file mode 100644
index 0000000..5626f13
--- /dev/null
+++ b/src/test/python/test_flimlib.py
@@ -0,0 +1,196 @@
+#unit tests
+import unittest
+import flimlib
+import numpy as np
+import math
+
+samples = 256
+a_in = 10.0
+tau_in = 2.0
+period = 0.04
+time = np.linspace(0,(samples-1)*period,samples,dtype=np.float32)
+photon_count32 = a_in * np.exp(-time/tau_in)
+photon_count64 = np.asarray(photon_count32,dtype=np.float64)
+photon_count2d = np.zeros([256,256])
+
+linear_const = 1
+photon_count_linear = linear_const * time
+
+class Test3Integral(unittest.TestCase):
+    def test_output_margin(self):
+        result = flimlib.GCI_triple_integral_fitting_engine(
+            period, photon_count32)
+        
+        # test to an integer precision since the triple integral fit isn't as rigorous
+        self.assertAlmostEqual(result.A, a_in, 0)
+        self.assertAlmostEqual(result.tau, tau_in, 0)
+        self.assertAlmostEqual(result.Z, 0.0, 0)
+        self.assertTrue(result.tries>0)
+        self.assertAlmostEqual(result.fitted[0], a_in, 0)
+        self.assertAlmostEqual(result.residuals[0], 0.0, 0)
+        self.assertAlmostEqual(result.residuals[-1], 0.0, 0)
+
+
+    def test_photon_count(self):
+        flimlib.GCI_triple_integral_fitting_engine(period, photon_count64) #passing float 64 should be fine
+        with self.assertRaises(ValueError):
+            flimlib.GCI_triple_integral_fitting_engine(period, photon_count2d) #passing 2d array will fail
+        with self.assertRaises(ValueError):
+            flimlib.GCI_triple_integral_fitting_engine(period, 5.7) #must be 1d array
+        flimlib.GCI_triple_integral_fitting_engine(period, [1,2])
+
+    def test_instr(self):
+        flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, instr=[1,2,3,4,5])
+        flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, instr=["42"])
+        flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, instr=list(range(300))) # no errors??? ask mark/jenu. ok we should actually just check that it's the right length
+        with self.assertRaises(ValueError):
+            flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, instr=5)
+        with self.assertRaises(ValueError):
+            flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, instr="foo")
+        with self.assertRaises(TypeError):
+            flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, instr={"foo": "bar"})
+
+    def test_noise_type_input(self):
+        with self.assertRaises(ValueError):
+            flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='foo')
+        with self.assertRaises(TypeError):
+            flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type={"foo": "bar"})
+
+    def test_unused_noise_types(self):
+        with self.assertRaises(ValueError):
+            flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='NOISE_GAUSSIAN_FIT')
+        with self.assertRaises(ValueError):
+            flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='NOISE_MLE')
+
+    def test_noise_const(self):
+        flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='NOISE_CONST',sig=2)
+        flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='NOISE_CONST',sig=[5])
+        with self.assertRaises(TypeError):
+            flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='NOISE_CONST',sig=[1,2,3,4,5])
+        with self.assertRaises(ValueError):
+            flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='NOISE_CONST',sig='foo')
+    
+    def test_noise_given(self):
+        flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='NOISE_GIVEN',sig=list(range(samples)))
+        with self.assertRaises(ValueError):
+            flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='NOISE_GIVEN',sig=2)
+        with self.assertRaises(TypeError):
+            flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='NOISE_GIVEN',sig={"foo": "bar"})
+        with self.assertRaises(ValueError):
+            flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='NOISE_GIVEN',sig='foo')
+        with self.assertRaises(ValueError):
+            flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='NOISE_GIVEN',sig=[1,2,3,4,5])
+        
+    def test_noise_const_and_given(self):
+        #the result should be the same if noise given is constant
+        result_const = flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='NOISE_CONST',sig=1.0)
+        result_noise_given = flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='NOISE_GIVEN',sig=np.ones(samples))
+        self.assertAlmostEqual(result_const.chisq,result_noise_given.chisq)
+        self.assertAlmostEqual(result_const.A,result_noise_given.A)
+        self.assertAlmostEqual(result_const.Z,result_noise_given.Z)
+        self.assertAlmostEqual(result_const.tau,result_noise_given.tau)
+
+    def test_noise_poisson_data(self):
+        flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='NOISE_POISSON_DATA')
+        # this should print a warning
+        flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='NOISE_POISSON_DATA',sig=2)
+
+    def test_noise_poisson_fit(self):
+        flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='NOISE_POISSON_FIT')
+        # this should print a warning
+        flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, noise_type='NOISE_POISSON_FIT',sig=2)
+    
+    def test_chisq_target(self):
+        result = flimlib.GCI_triple_integral_fitting_engine(period, photon_count32, chisq_target=0.0) #impossible target causes it to run more than once and give up
+        self.assertTrue(result.tries>1)
+
+class TestGCIPhasor(unittest.TestCase):
+    def test_output_margin(self):
+        result = flimlib.GCI_Phasor(period, photon_count32)
+        self.assertEqual(result.error_code, 0)
+        self.assertAlmostEqual(result.fitted[0], a_in, 1)
+        self.assertAlmostEqual(result.residuals[0], 0.0, 1)
+
+
+def dummy_exp_tau(x, param):
+    y = param[1]*np.exp(-x / param[2])
+    #y = param[0] + param[1]*np.exp(-x / param[2])
+    dy_dparam = [1,0,0]
+    dy_dparam[1] = np.exp(-x/param[2])
+    dy_dparam[2] = param[1]*x*np.exp(-x/param[2])/(param[2]**2)
+    return y, dy_dparam
+
+def dummy_exp_tau_predicate(n_param):
+    return n_param == 3
+
+# linear fit!
+def dummy_linear(x, param):
+    y = param[1]*x
+    #y = param[0] + param[1]*x
+    dy_dparam = [1,0]
+    dy_dparam[1] = x
+    return y, dy_dparam
+
+def dummy_linear_predicate(n_param):
+    return n_param == 2
+
+class TestGCIMarquardt(unittest.TestCase):
+    def test_output_margin(self):
+        param_in = [0,a_in+1,tau_in+1] # slight offset to detect if the fitting works!
+        result = flimlib.GCI_marquardt_fitting_engine(period, photon_count32, param_in)
+        self.assertAlmostEqual(0.0,result.param[0],1)
+        self.assertAlmostEqual(a_in,result.param[1],1)
+        self.assertAlmostEqual(tau_in,result.param[2],1)
+
+    def test_paramfree(self):
+        param_in = [0,a_in+1,tau_in+1] # slight offset to detect if the fitting works!
+        result = flimlib.GCI_marquardt_fitting_engine(period, photon_count32, 
+            param_in, paramfree=[1,0,1])
+        # second parameter should have been held fixed!
+        self.assertAlmostEqual(a_in+1, result.param[1],1)
+
+    def test_restraintype(self):
+        with self.assertRaises(ValueError):
+            flimlib.GCI_set_restrain_limits([0,0],[0],[0])
+        with self.assertRaises(ValueError):
+            flimlib.GCI_set_restrain_limits([[0,0],[0,0]],[0,0],[0,0])
+        param_in = [0,a_in+1,tau_in+1]
+        flimlib.GCI_set_restrain_limits([0,1,0],[0,0,0],[0,a_in-1,0])
+        result = flimlib.GCI_marquardt_fitting_engine(period, photon_count32, param_in, restrain_type='ECF_RESTRAIN_USER')
+        self.assertTrue(result.param[1] <= a_in-1)
+    
+    def test_multiexp_lambda(self):
+        lambda_in = 1/tau_in # lambda is the decay rate!
+        param_in = [0,a_in+1,lambda_in+1] # slight offset to detect if the fitting works!
+        result = flimlib.GCI_marquardt_fitting_engine(period, photon_count32, param_in, fitfunc=flimlib.GCI_multiexp_lambda)
+
+        self.assertAlmostEqual(a_in, result.param[1], 1)
+        self.assertAlmostEqual(lambda_in, result.param[2],1)
+    
+    def test_user_defined_fitfunc(self):
+        param_in = [0.0,a_in+1,tau_in+1] # slight offset to detect if the fitting works!
+        fitfunc_in = flimlib.FitFunc(dummy_exp_tau, nparam_predicate=dummy_exp_tau_predicate)
+        result = flimlib.GCI_marquardt_fitting_engine(period, photon_count32, param_in, fitfunc=fitfunc_in)
+        self.assertAlmostEqual(0.0,result.param[0],1)
+        self.assertAlmostEqual(a_in,result.param[1],1)
+        self.assertAlmostEqual(tau_in,result.param[2],1)
+
+        # linear fit! (did not work with the first param being nonzero)
+        param_in = [0, linear_const+1]
+        fitfunc_in = flimlib.FitFunc(dummy_linear, nparam_predicate=dummy_linear_predicate)
+        result = flimlib.GCI_marquardt_fitting_engine(period, photon_count_linear, param_in, fitfunc=fitfunc_in, noise_type='NOISE_CONST',sig=1.0)
+        self.assertAlmostEqual(0.0,result.param[0],1)
+        self.assertAlmostEqual(linear_const,result.param[1],1)
+
+    def test_nparam(self):
+        with self.assertRaises(TypeError):
+            param_in = [a_in+1,tau_in+1] # forgot the first parameter!
+            flimlib.GCI_marquardt_fitting_engine(period, photon_count_linear, param_in)
+
+        # any odd number of parameters greater or equal to 3 should work!
+        param_in = [0, a_in+1,tau_in+1, 1, 1] # pass 5 parameters
+        flimlib.GCI_marquardt_fitting_engine(period, photon_count_linear, param_in)
+        
+
+if __name__ == '__main__':
+    unittest.main()
\ No newline at end of file