def halton_seq(self,output_array, start_dim, seed):
        array_1d_double = npct.ndpointer(dtype=numpy.double, ndim=1, flags='CONTIGUOUS')
        array_2d_double = npct.ndpointer(dtype=numpy.double, ndim=2, flags='CONTIGUOUS')

        self.so.halton_seq.argtypes = [c_int, c_int, c_int, c_int, array_2d_double]
        self.so.halton_seq.restype = c_int
        return self.so.halton_seq(len(output_array), len(output_array[0]), start_dim, seed, output_array)

    def solve(self, b, verbose = None):
        b_shp = b.shape
        if b.ndim == 2 and b.shape[1] == 1:
            b = b.ravel()
            nrhs = 1
        
        elif b.ndim == 2 and b.shape[1] != 1:
            nrhs = b.shape[1]
            b = b.ravel(order='F')
        
        else:
            b = b.ravel()
            nrhs = 1

        if self._is_complex:
            data_type = np.complex128
            if b.dtype != np.complex128:
                b = b.astype(np.complex128, copy=False)
        else:
            data_type = np.float64
            if b.dtype != np.float64:
                b = b.astype(np.float64, copy=False)

        # Create solution array (x) and pointers to x and b
        if self._is_complex:
            x = np.zeros(b.shape, dtype=np.complex128, order='C')
        else:
            x = np.zeros(b.shape, dtype=np.float64, order='C')
        
        np_x = x.ctypes.data_as(ctypeslib.ndpointer(data_type, ndim=1, flags='C')) 
        np_b = b.ctypes.data_as(ctypeslib.ndpointer(data_type, ndim=1, flags='C'))
        
        error = np.zeros(1,dtype=np.int32)
        np_error = error.ctypes.data_as(ctypeslib.ndpointer(np.int32, ndim=1, flags='C')) 
        #Call solver
        _solve_start = time.time()
        pardiso(self._np_pt, byref(c_int(1)), byref(c_int(1)), byref(c_int(self._mtype)),
            byref(c_int(33)), byref(c_int(self._dim)), self._data, self._indptr, self._indices, 
            self._np_perm, byref(c_int(nrhs)), self._np_iparm, byref(c_int(0)), np_b,
            np_x, np_error)
        self._solve_time = time.time() -_solve_start
        if error[0] != 0:
            raise Exception(pardiso_error_msgs[str(error[0])])
        
        if verbose:
            print('Solution Stage')
            print('--------------')
            print('Solution time:                  ',round(self._solve_time,4))
            print('Solution memory (Mb):           ',round(self._iparm[16]/1024.,4))
            print('Number of iterative refinements:',self._iparm[6])
            print('Total memory (Mb):              ',round(sum(self._iparm[15:17])/1024.,4))
            print()
        
        # Return solution vector x
        if nrhs==1:
            if x.shape != b_shp:
                x = np.reshape(x, b_shp)
            return x
        else:
            return np.reshape(x, b_shp, order='F')

def loadPupilFitC():
    pupil_fit = loadclib.loadCLibrary("pupil_fit")

    # From sa_library/multi_fit.c
    pupil_fit.mFitGetFitImage.argtypes = [ctypes.c_void_p,
                                          ndpointer(dtype=numpy.float64)]

    pupil_fit.mFitGetNError.argtypes = [ctypes.c_void_p]
    pupil_fit.mFitGetNError.restype = ctypes.c_int
    
    pupil_fit.mFitGetPeakPropertyDouble.argtypes = [ctypes.c_void_p,
                                                    ndpointer(dtype=numpy.float64),
                                                    ctypes.c_char_p]
    
    pupil_fit.mFitGetPeakPropertyInt.argtypes = [ctypes.c_void_p,
                                                 ndpointer(dtype=numpy.int32),
                                                 ctypes.c_char_p]

    pupil_fit.mFitGetResidual.argtypes = [ctypes.c_void_p,
                                          ndpointer(dtype=numpy.float64)]
    
    pupil_fit.mFitGetUnconverged.argtypes = [ctypes.c_void_p]
    pupil_fit.mFitGetUnconverged.restype = ctypes.c_int

    pupil_fit.mFitIterateLM.argtypes = [ctypes.c_void_p]

    pupil_fit.mFitNewBackground.argtypes = [ctypes.c_void_p,
                                            ndpointer(dtype=numpy.float64)]    
    
    pupil_fit.mFitNewImage.argtypes = [ctypes.c_void_p,
                                       ndpointer(dtype=numpy.float64)]

    pupil_fit.mFitRemoveErrorPeaks.argtypes = [ctypes.c_void_p]

    pupil_fit.mFitRemoveRunningPeaks.argtypes = [ctypes.c_void_p]

    pupil_fit.mFitSetPeakStatus.argtypes = [ctypes.c_void_p,
                                            ndpointer(dtype=numpy.int32)]       
    
    # From pupilfn/pupil_fit.c
    pupil_fit.pfitCleanup.argtypes = [ctypes.c_void_p]

    pupil_fit.pfitInitialize.argtypes = [ctypes.c_void_p,
                                         ndpointer(dtype=numpy.float64),
                                         ndpointer(dtype=numpy.float64),
                                         ctypes.c_double,
                                         ctypes.c_int,
                                         ctypes.c_int]
    pupil_fit.pfitInitialize.restype = ctypes.POINTER(daoFitC.fitData)

    pupil_fit.pfitNewPeaks.argtypes = [ctypes.c_void_p,
                                       ndpointer(dtype=numpy.float64),
                                       ctypes.c_char_p,
                                       ctypes.c_int]

    pupil_fit.pfitSetZRange.argtypes = [ctypes.c_void_p,
                                        ctypes.c_double,
                                        ctypes.c_double]

    return pupil_fit

    def __init__(self, scheme="midpoint", solver="hll"):
        from ctypes import CDLL, POINTER, CFUNCTYPE, Structure, c_double, c_int
        from numpy import float64, int32
        from numpy.ctypeslib import ndpointer
        from os.path import abspath, dirname
        from os import popen

        dbl_arr = ndpointer(dtype=float64, flags=("C_CONTIGUOUS", "WRITEABLE"))
        dbl_vec = ndpointer(dtype=float64, flags=("C_CONTIGUOUS"))
        int_arr = ndpointer(dtype=int32, flags=("C_CONTIGUOUS", "WRITEABLE"))
        int_vec = ndpointer(dtype=int32, flags=("C_CONTIGUOUS"))

        lib_home = dirname(abspath(popen("find . -name *.so").readline()))
        clib = CDLL(lib_home + "/" + self.libname + ".so")

        self.schemes = ["fwd_euler", "midpoint", "RK3", "ctu_hancock"]
        self.ghost_cells = dict(zip(self.schemes, (2, 4, 6, 4)))
        self.advance = {}

        assert scheme in self.schemes
        assert solver in self.solvers

        for sname in self.schemes:
            self.advance[sname] = getattr(clib, "advance_state_" + sname)
            self.advance[sname].argtypes = [dbl_arr, c_double]

        clib.integrate_init.argtypes = [int_vec, dbl_vec, c_int]
        clib.integrate_free.argtypes = []
        clib.get_failure_mask.argtypes = [int_arr]
        clib.set_riemann_solver.argtypes = [c_int]

        self.clib = clib
        self.scheme = scheme
        self.NumGhostCells = self.ghost_cells[self.scheme]
        self.solver = solver

def setCInterface(homotopy_lib):
    global directory
    global homotopy

    if sys.platform == "win32":
        homotopy = cdll.LoadLibrary(directory + homotopy_lib + ".dll")
    else:
        homotopy = cdll.LoadLibrary(directory + homotopy_lib + ".so")

    l1flt_size = homotopy.getL1FLTSize()
    if l1flt_size == 4:
        float_type = c_float
        float_array_type = numpy.float32
    elif l1flt_size == 8:
        float_type = c_double
        float_array_type = numpy.float64

    homotopy.getXVector.argtypes = [ndpointer(dtype=float_array_type)]
    homotopy.initialize.argtypes = [ndpointer(dtype=float_array_type), c_int, c_int, c_int, c_int, c_int]
    homotopy.l2Error.argtypes = [ndpointer(dtype=float_array_type)]
    homotopy.l2Error.restype = float_type
    homotopy.newYVector.argtypes = [ndpointer(dtype=float_array_type)]
    homotopy.solve.argtypes = [float_type, c_int]
    homotopy.solve.restype = float_type

    if hasattr(homotopy, "getVisited"):
        homotopy.getVisited.argtypes = [ndpointer(dtype=numpy.int32)]

    if hasattr(homotopy, "initializeGPU"):
        homotopy.initializeGPU.argtypes = [c_char_p, c_int, c_int, c_int]

def resize(img,scale):
    """
        downsample img to scale 
    """
    sdims=img.shape
    datatype=c_double
    if img.dtype!=datatype:
       print "Error the image must be of doubles!"
       raise RuntimeError

    if scale>1.0:
       print "Invalid scaling factor!"
       raise RuntimeError  
    
    img = asfortranarray(img,c_double) # make array continguous
    
    try:
        mresize = ctypeslib.load_library("libresize.so",".") 
    except:
        print "Unable to load resize library"
        raise RuntimeError
        
    #use two times the 1d resize to get a 2d resize
    fresize = mresize.resize1dtran
    fresize.restype = None
    fresize.argtypes = [ ctypeslib.ndpointer(dtype=datatype, ndim=3), c_int,ctypeslib.ndpointer(dtype=datatype, ndim=3), c_int, c_int , c_int ]
    ddims = [int(round(sdims[0]*scale)),int(round(sdims[1]*scale)),sdims[2]];
    mxdst = zeros((ddims), dtype=datatype)
    tmp = zeros((ddims[0],sdims[1],sdims[2]), dtype=datatype)
    img1=img
    t1=time()
    fresize(img1, sdims[0], tmp, ddims[0], sdims[1], sdims[2]);
    fresize(tmp, sdims[1], mxdst, ddims[1], ddims[0], sdims[2]);
    t2=time()
    return mxdst.reshape(ddims[2],ddims[1],ddims[0]).T

 def test_dtype(self):
     dt = np.intc
     p = ndpointer(dtype=dt)
     self.assertTrue(p.from_param(np.array([1], dt)))
     dt = '<i4'
     p = ndpointer(dtype=dt)
     self.assertTrue(p.from_param(np.array([1], dt)))
     dt = np.dtype('>i4')
     p = ndpointer(dtype=dt)
     p.from_param(np.array([1], dt))
     self.assertRaises(TypeError, p.from_param,
                       np.array([1], dt.newbyteorder('swap')))
     dtnames = ['x', 'y']
     dtformats = [np.intc, np.float64]
     dtdescr = {'names' : dtnames, 'formats' : dtformats}
     dt = np.dtype(dtdescr)
     p = ndpointer(dtype=dt)
     self.assertTrue(p.from_param(np.zeros((10,), dt)))
     samedt = np.dtype(dtdescr)
     p = ndpointer(dtype=samedt)
     self.assertTrue(p.from_param(np.zeros((10,), dt)))
     dt2 = np.dtype(dtdescr, align=True)
     if dt.itemsize != dt2.itemsize:
         self.assertRaises(TypeError, p.from_param, np.zeros((10,), dt2))
     else:
         self.assertTrue(p.from_param(np.zeros((10,), dt2)))

def kmeans(X, nclst, maxiter=0, numruns=1):
    """Wrapper for Peter Gehlers accelerated MPI-Kmeans routine."""

    mpikmeanslib = N.ctypeslib.load_library("libmpikmeans.so", ".")
    mpikmeanslib.kmeans.restype = c_double
    mpikmeanslib.kmeans.argtypes = [
        ndpointer(dtype=c_double, ndim=1, flags="C_CONTIGUOUS"),
        ndpointer(dtype=c_double, ndim=1, flags="C_CONTIGUOUS"),
        ndpointer(dtype=c_uint, ndim=1, flags="C_CONTIGUOUS"),
        c_uint,
        c_uint,
        c_uint,
        c_uint,
        c_uint,
    ]

    npts, dim = X.shape
    assignments = empty((npts), c_uint)

    #    bestSSE=N.Inf
    #    bestassignments=empty( (npts), c_uint)
    Xvec = array(reshape(X, (-1,)), c_double, order="C")
    permutation = N.random.permutation(range(npts))  # randomize order of points
    CX = array(X[permutation[:nclst], :], c_double, order="C").flatten()
    print CX
    SSE = mpikmeanslib.kmeans(CX, Xvec, assignments, dim, npts, min(nclst, npts), maxiter, numruns)
    return reshape(CX, (nclst, dim)), SSE, (assignments + 1)

def pmc(ei,ej,nnodes,nnedges): #ei, ej is edge list whose index starts from 0
    degrees = np.zeros(nnodes,dtype = np.int32)
    new_ei = []
    new_ej = []
    for i in range(nnedges):
        degrees[ei[i]] += 1
        if ej[i] <= ei[i] + 1:
            new_ei.append(ei[i])
            new_ej.append(ej[i])
    maxd = max(degrees)
    offset = 0
    new_ei = np.array(new_ei,dtype = np.int32)
    new_ej = np.array(new_ej,dtype = np.int32)
    outsize = maxd
    output = np.zeros(maxd,dtype = np.int32)
    lib = ctypes.cdll.LoadLibrary("libpmc.dylib")
    fun = lib.max_clique
    #call C function
    fun.restype = np.int32
    fun.argtypes = [ctypes.c_int32,ndpointer(ctypes.c_int32, flags="C_CONTIGUOUS"),
                  ndpointer(ctypes.c_int32, flags="C_CONTIGUOUS"),ctypes.c_int32,
                  ctypes.c_int32,ndpointer(ctypes.c_int32, flags="C_CONTIGUOUS")]
    clique_size = fun(len(new_ei),new_ei,new_ej,offset,outsize,output)
    max_clique = np.empty(clique_size,dtype = np.int32)
    max_clique[:]=[output[i] for i in range(clique_size)]

    return max_clique

def loadXDRLibrary(LoadDirectFromMSMBuilder=True):
    global _xdrlib

    xdr_library_path = find_library("xdrfile")
    if xdr_library_path and not LoadDirectFromMSMBuilder:
        _xdrlib = CDLL(xdr_library_path)
    elif LoadDirectFromMSMBuilder==True:
        MSMBuilderPath=imp.find_module("msmbuilder")[1]
        _xdrlib=CDLL(MSMBuilderPath+"/libxdrfile.so")
    else:
        # Lutz: find_library does not look in LD_LIBRARY_PATH on linux machines (see http://bugs.python.org/issue2936), even though cdll.LoadLibrary does
        # as a workaround, we try to load the shared object file manually
        if sys.platform.startswith("linux"):
            try:
                _xdrlib = CDLL("libxdrfile.so")
            except:
                raise RuntimeError("Unable to find xdrfile library. Make sure that it is compiled as a shared library, and that its location is known to the dynamic linker (e.g., set LD_LIBRARY_PATH on Linux).")
        else:
            raise RuntimeError("Unable to find xdrfile library. Make sure that it is compiled as a shared library, and that its location is known to the dynamic linker (e.g., set DYLD_LIBRARY_PATH on Mac OS X, LD_LIBRARY_PATH on Linux).")

    # declare interface for functions in the xdr library to insure some amount of type safety
    # for use of numpy arrays in ctypes, see http://thread.gmane.org/gmane.comp.python.numeric.general/7418/focus=7418
    _xdrlib.read_xtc_natoms.argtypes = [c_char_p, POINTER(c_int)]
    _xdrlib.read_trr_natoms.argtypes = [c_char_p, POINTER(c_int)]
    _xdrlib.xdrfile_open.argtypes = [c_char_p, c_char_p]
    _xdrlib.xdrfile_open.restype = c_void_p
    _xdrlib.xdrfile_close.argtypes = [c_void_p]
    _xdrlib.read_xtc.argtypes = [c_void_p, c_int, POINTER(c_int), POINTER(c_float), ndpointer(dtype="single",shape=(3,3),flags="C_CONTIGUOUS"), ndpointer(dtype="single",ndim=2,flags="C_CONTIGUOUS"), POINTER(c_float)]
    _xdrlib.read_trr.argtypes = [c_void_p, c_int, POINTER(c_int), POINTER(c_float),POINTER(c_float), ndpointer(dtype="single",shape=(3,3),flags="C_CONTIGUOUS"), ndpointer(dtype="single",ndim=2,flags="C_CONTIGUOUS"),ndpointer(dtype="single",ndim=2,flags="C_CONTIGUOUS"),ndpointer(dtype="single",ndim=2,flags="C_CONTIGUOUS")]
    _xdrlib.write_xtc.argtypes = [c_void_p, c_int, c_int, c_float, ndpointer(dtype="single",shape=(3,3),flags="C_CONTIGUOUS"), ndpointer(dtype="single",ndim=2,flags="C_CONTIGUOUS"), c_float]

    def __init__(self, rho, sigma, scorebook):
        """
        :param rho: gap opening penalty
        :param sigma: gap extension penalty
        :param scorebook: see `from_submatr` method
        """
        self.rho = rho
        self.sigma = sigma
        self.scorebook = scorebook

        arrflages = str('C')

        self.lib = ctypes.cdll.LoadLibrary(locate_lib('align.so'))
        self.build = self.lib.build
        self.build.argtypes = [
            ndpointer(self.score_type, flags=arrflages),
            self.index_type,   # isize
            self.index_type,   # jsize
            self.score_type,   # rho
            self.score_type,   # sigma
            self.short_type,  # local_
        ]

        self.backtrack = self.lib.backtrack
        self.backtrack.rettype = self.index_type
        self.backtrack.argtypes = [
            ndpointer(self.score_type, flags=arrflages),
            self.index_type,   # isize
            self.index_type,   # jsize
            self.index_type,   # istart
            self.index_type,   # jstart
            ndpointer(self.index_type, flags=arrflages),    # iarr index
            ndpointer(self.index_type, flags=arrflages),    # jarr index
            self.short_type,  # global_
        ]

def setCInterface(homotopy_lib):
    global homotopy

    homotopy = loadclib.loadCLibrary(homotopy_lib)

    l1flt_size = homotopy.getL1FLTSize()
    if(l1flt_size == 4):
        float_type = c_float
        float_array_type = numpy.float32
    elif(l1flt_size == 8):
        float_type = c_double
        float_array_type = numpy.float64

    homotopy.getXVector.argtypes = [ndpointer(dtype=float_array_type)]
    homotopy.initialize.argtypes = [ndpointer(dtype=float_array_type),
                                    c_int,
                                    c_int,
                                    c_int,
                                    c_int,
                                    c_int]
    homotopy.l2Error.argtypes = [ndpointer(dtype=float_array_type)]
    homotopy.l2Error.restype = float_type
    homotopy.newYVector.argtypes = [ndpointer(dtype=float_array_type)]
    homotopy.solve.argtypes = [float_type, c_int]
    homotopy.solve.restype = float_type

    if(hasattr(homotopy, "getVisited")):
        homotopy.getVisited.argtypes = [ndpointer(dtype=numpy.int32)]

    if(hasattr(homotopy, "initializeGPU")):
        homotopy.initializeGPU.argtypes = [c_char_p,
                                           c_int,
                                           c_int,
                                           c_int]

    def __init__(self, nx, ny):
        self.nx = nx
        self.ny = ny

        self.nx_p = byref(c_int(nx))
        self.ny_p = byref(c_int(ny))

        # Load the library using numpy
        libm = npct.load_library('my_math', './')
        modname = 'my_math'
        my_cos = getattr(libm, '__{}_MOD_{}'.format(modname, 'my_cos'))
        my_cos_2d = getattr(libm, '__{}_MOD_{}'.format(modname, 'my_cos_2d'))

        # Set the argument and return type
        c_int_p = POINTER(c_int)

        arr_1d_f8 = npct.ndpointer(ndim=1, dtype='f8')
        my_cos.argtypes = (c_int_p, arr_1d_f8, arr_1d_f8)
        my_cos.restype = None

        arr_2d_f8 = npct.ndpointer(ndim=2, dtype='f8', flags='F')
        my_cos_2d.argtypes = (c_int_p, c_int_p, arr_2d_f8, arr_2d_f8)
        my_cos_2d.restype = None

        # Set to public
        self.libm = libm
        self.my_cos = my_cos
        self.my_cos_2d = my_cos_2d

def does_mapping_exist(v, this_type, atom_pos, atomtype, eps):
    """No XML documentation summary.
    """
    libpath = path.join(path.dirname(__file__), "ftypes.celib.so")
    method = static_symbol("symmetry_module_c", "does_mapping_exist_c", libpath, True)
    method.argtypes = [ndpointer(dtype=float, ndim=1, shape=(3,), flags="F"), c_int_p, 
                       ndpointer(dtype=float, ndim=2, flags="F"), c_int_p, c_int_p, 
                       ndpointer(dtype=int, ndim=1, flags="F"), c_int_p, c_bool_p, 
                       c_double_p]

    v_a = require(v, float, "F")
    this_type_c = c_int(this_type)
    atom_pos_a = require(atom_pos, float, "F")
    atom_pos_0 = c_int(atom_pos_a.shape[0])
    atom_pos_1 = c_int(atom_pos_a.shape[1])
    atomtype_a = require(atomtype, int, "F")
    atomtype_0 = c_int(atomtype_a.shape[0])
    mapped_c = c_bool(False)
    eps_c = c_double(eps)
    method(v_a, byref(this_type_c), atom_pos_a, byref(atom_pos_0), byref(atom_pos_1), 
           atomtype_a, byref(atomtype_0), byref(mapped_c), byref(eps_c))

    result = FtypesResult("symmetry_module", "does_mapping_exist", "Subroutine")
    result.add("mapped", mapped_c.value)

    return result

 def test_dtype(self):
     dt = np.intc
     p = ndpointer(dtype=dt)
     self.assertTrue(p.from_param(np.array([1], dt)))
     dt = "<i4"
     p = ndpointer(dtype=dt)
     self.assertTrue(p.from_param(np.array([1], dt)))
     dt = np.dtype(">i4")
     p = ndpointer(dtype=dt)
     p.from_param(np.array([1], dt))
     self.assertRaises(TypeError, p.from_param, np.array([1], dt.newbyteorder("swap")))
     dtnames = ["x", "y"]
     dtformats = [np.intc, np.float64]
     dtdescr = {"names": dtnames, "formats": dtformats}
     dt = np.dtype(dtdescr)
     p = ndpointer(dtype=dt)
     self.assertTrue(p.from_param(np.zeros((10,), dt)))
     samedt = np.dtype(dtdescr)
     p = ndpointer(dtype=samedt)
     self.assertTrue(p.from_param(np.zeros((10,), dt)))
     dt2 = np.dtype(dtdescr, align=True)
     if dt.itemsize != dt2.itemsize:
         self.assertRaises(TypeError, p.from_param, np.zeros((10,), dt2))
     else:
         self.assertTrue(p.from_param(np.zeros((10,), dt2)))

def actionAngleTorus_Freqs_c(pot,jr,jphi,jz,
                             tol=0.003):
    """
    NAME:
       actionAngleTorus_Freqs_c
    PURPOSE:
       compute frequencies on a single torus
    INPUT:
       pot - Potential object or list thereof
       jr - radial action (scalar)
       jphi - azimuthal action (scalar)
       jz - vertical action (scalar)
       tol= (0.003) goal for |dJ|/|J| along the torus
    OUTPUT:
       (Omegar,Omegaphi,Omegaz,flag)
    HISTORY:
       2015-08-05/07 - Written - Bovy (UofT)
    """
    #Parse the potential
    npot, pot_type, pot_args= _parse_pot(pot,potfortorus=True)

    #Set up result
    Omegar= numpy.empty(1)
    Omegaphi= numpy.empty(1)
    Omegaz= numpy.empty(1)
    flag= ctypes.c_int(0)

    #Set up the C code
    ndarrayFlags= ('C_CONTIGUOUS','WRITEABLE')
    actionAngleTorus_FreqsFunc= _lib.actionAngleTorus_Freqs
    actionAngleTorus_FreqsFunc.argtypes=\
        [ctypes.c_double,
         ctypes.c_double,
         ctypes.c_double,
         ctypes.c_int,
         ndpointer(dtype=numpy.int32,flags=ndarrayFlags),
         ndpointer(dtype=numpy.float64,flags=ndarrayFlags),
         ctypes.c_double,
         ndpointer(dtype=numpy.float64,flags=ndarrayFlags),
         ndpointer(dtype=numpy.float64,flags=ndarrayFlags),
         ndpointer(dtype=numpy.float64,flags=ndarrayFlags),
         ctypes.POINTER(ctypes.c_int)]

    #Array requirements
    Omegar= numpy.require(Omegar,dtype=numpy.float64,requirements=['C','W'])
    Omegaphi= numpy.require(Omegaphi,dtype=numpy.float64,requirements=['C','W'])
    Omegaz= numpy.require(Omegaz,dtype=numpy.float64,requirements=['C','W'])
    
    #Run the C code
    actionAngleTorus_FreqsFunc(ctypes.c_double(jr),
                               ctypes.c_double(jphi),
                               ctypes.c_double(jz),
                               ctypes.c_int(npot),
                               pot_type,
                               pot_args,
                               ctypes.c_double(tol),
                               Omegar,Omegaphi,Omegaz,
                               ctypes.byref(flag))

    return (Omegar[0],Omegaphi[0],Omegaz[0],flag.value)

def actionAngleStaeckel_calcu0(E,Lz,pot,delta):
    """
    NAME:
       actionAngleStaeckel_calcu0
    PURPOSE:
       Use C to calculate u0 in the Staeckel approximation
    INPUT:
       E, Lz - energy and angular momentum
       pot - Potential or list of such instances
       delta - focal length of prolate spheroidal coordinates
    OUTPUT:
       (u0,err)
       u0 : array, shape (len(E))
       err - non-zero if error occured
    HISTORY:
       2012-12-03 - Written - Bovy (IAS)
    """
    #Parse the potential
    npot, pot_type, pot_args= _parse_pot(pot,potforactions=True)

    #Set up result arrays
    u0= numpy.empty(len(E))
    err= ctypes.c_int(0)

    #Set up the C code
    ndarrayFlags= ('C_CONTIGUOUS','WRITEABLE')
    actionAngleStaeckel_actionsFunc= _lib.calcu0
    actionAngleStaeckel_actionsFunc.argtypes= [ctypes.c_int,
                               ndpointer(dtype=numpy.float64,flags=ndarrayFlags),
                               ndpointer(dtype=numpy.float64,flags=ndarrayFlags),
                               ctypes.c_int,
                               ndpointer(dtype=numpy.int32,flags=ndarrayFlags),
                               ndpointer(dtype=numpy.float64,flags=ndarrayFlags),
                               ctypes.c_double,
                               ndpointer(dtype=numpy.float64,flags=ndarrayFlags),
                               ctypes.POINTER(ctypes.c_int)]

    #Array requirements, first store old order
    f_cont= [E.flags['F_CONTIGUOUS'],
             Lz.flags['F_CONTIGUOUS']]
    E= numpy.require(E,dtype=numpy.float64,requirements=['C','W'])
    Lz= numpy.require(Lz,dtype=numpy.float64,requirements=['C','W'])
    u0= numpy.require(u0,dtype=numpy.float64,requirements=['C','W'])

    #Run the C code
    actionAngleStaeckel_actionsFunc(len(E),
                                    E,
                                    Lz,
                                    ctypes.c_int(npot),
                                    pot_type,
                                    pot_args,
                                    ctypes.c_double(delta),
                                    u0,
                                    ctypes.byref(err))

    #Reset input arrays
    if f_cont[0]: E= numpy.asfortranarray(E)
    if f_cont[1]: Lz= numpy.asfortranarray(Lz)

    return (u0,err.value)

    def __init__(self, method='full', sampling=(-1,-1), rank=5, reg=0.01, tol=1E-6, iters=500, verbose=True):

        modes = {'full':0 , 'subsample': 1}
        if method not in modes:
            raise ValueError("'method' must be one of" + modes.keys())
        self.method = method
        self._mode = modes[method]
        if method == 'subsample' and -1 in sampling:
            raise ValueError("'method' is set to 'subsample' but 'sampling' is not set.")
        self.sampling = sampling
        self.rank = rank
        self.reg = reg
        self.tol = tol
        self.iters = iters
        self.verbose = verbose
        libpath = os.path.dirname(os.path.abspath(__file__)) + '/librosl.so.0.2'
        self._pyrosl = ctypes.cdll.LoadLibrary(libpath).pyROSL
        self._pyrosl.restype = ctypes.c_int
        self._pyrosl.argtypes = [
                           ndpointer(ctypes.c_double, flags="F_CONTIGUOUS"),
                           ndpointer(ctypes.c_double, flags="F_CONTIGUOUS"),
                           ndpointer(ctypes.c_double, flags="F_CONTIGUOUS"),
                           ndpointer(ctypes.c_double, flags="F_CONTIGUOUS"),
                           ctypes.c_int, ctypes.c_int,
                           ctypes.c_int, ctypes.c_double,
                           ctypes.c_double, ctypes.c_int,
                           ctypes.c_int, ctypes.c_int,
                           ctypes.c_int, ctypes.c_bool]
        self.components_ = None

def test_populate_z_shells(clib, formal_integral_model, p):
    '''
    Test the case where p > r[0]
    '''
    func = clib.populate_z
    func.restype = ctypes.c_int64
    func.argtypes = [
            ctypes.POINTER(StorageModel),   # storage
            c_double,                       # p
            ndpointer(dtype=np.float64),    # oz
            ndpointer(dtype=np.int64)       # oshell_id
            ]

    size = formal_integral_model.no_of_shells_i
    r_inner = as_array(formal_integral_model.r_inner_i, (size,))
    r_outer = as_array(formal_integral_model.r_outer_i, (size,))

    p = r_inner[0] + (r_outer[-1] - r_inner[0]) * p
    idx = np.searchsorted(r_outer, p, side='right')

    oz = np.zeros(size * 2)
    oshell_id = np.zeros_like(oz, dtype=np.int64)

    offset = size - idx

    expected_N = (offset) * 2
    expected_oz = np.zeros_like(oz)
    expected_oshell_id = np.zeros_like(oshell_id)

    # Calculated way to determine which shells get hit
    expected_oshell_id[:expected_N] = np.abs(
            np.arange(0.5, expected_N, 1) - offset) - 0.5 + idx

    expected_oz[0:offset] = 1 + calculate_z(
            r_outer[np.arange(size, idx, -1) - 1],
            p)
    expected_oz[offset:expected_N] = 1 - calculate_z(
            r_outer[np.arange(idx, size, 1)],
            p)

    N = func(
            formal_integral_model,
            p,
            oz,
            oshell_id
            )

    assert N == expected_N

    ntest.assert_allclose(
            oshell_id,
            expected_oshell_id
            )

    ntest.assert_allclose(
            oz,
            expected_oz,
            atol=1e-5
            )

 def __init__(self, dbname, host, user, password):
     self.connect = "dbname=%s host=%s user=%s password=%s"%(dbname, host, user, password) #TODO Harcode
     
     self.array_1d_double = npct.ndpointer(dtype=np.double, ndim=1, flags='CONTIGUOUS')
     self.array_1d_int = npct.ndpointer(dtype=np.int64, ndim=1, flags='CONTIGUOUS')
     self.array_2d_double = npct.ndpointer(dtype=np.double, ndim=2, flags='CONTIGUOUS')
     #TODO path harcode
     self.libcd = npct.load_library("cfsfdp", "/home/sergio/iibm/workspace/NeuroDB/NeuroDB/cfunctions/cfsfdp")