def compute_eccentric_anomaly(self, eccentricity, mean_anomaly):
        """compute eccentric anomaly, solve for Kepler Equation
            Parameter
            ----------
            eccentricity : array_like
                Eccentricity of binary system
            mean_anomaly : array_like
                Mean anomaly of the binary system
            Returns
            -------
            array_like
                The eccentric anomaly in radians, given a set of mean_anomalies
                in radians.
        """
        if hasattr(eccentricity,'unit'):
            e = np.longdouble(eccentricity).value
        else:
            e = eccentricity

        if any(e<0) or any(e>=1):
            raise ValueError('Eccentricity should be in the range of [0,1).')

        if hasattr(mean_anomaly,'unit'):
            ma = np.longdouble(mean_anomaly).value
        else:
            ma = mean_anomaly
        k = lambda E: E-e*np.sin(E)-ma   # Kepler Equation
        dk = lambda E: 1-e*np.cos(E)     # derivative Kepler Equation
        U = ma
        while(np.max(abs(k(U)))>5e-15):  # Newton-Raphson method
            U = U-k(U)/dk(U)
        return U*u.rad

    def __init__( self ):

        self.s = longdouble( 0 )
        self.m = longdouble( 0 )
        self.last_m = longdouble( 0 )
        self.n = ulonglong( 0 )
        self.is_started = False

def test_as_int():
    # Integer representation of number
    assert_equal(as_int(2.0), 2)
    assert_equal(as_int(-2.0), -2)
    assert_raises(FloatingError, as_int, 2.1)
    assert_raises(FloatingError, as_int, -2.1)
    assert_equal(as_int(2.1, False), 2)
    assert_equal(as_int(-2.1, False), -2)
    v = np.longdouble(2**64)
    assert_equal(as_int(v), 2**64)
    # Have all long doubles got 63+1 binary bits of precision?  Windows 32-bit
    # longdouble appears to have 52 bit precision, but we avoid that by checking
    # for known precisions that are less than that required
    try:
        nmant = type_info(np.longdouble)['nmant']
    except FloatingError:
        nmant = 63 # Unknown precision, let's hope it's at least 63
    v = np.longdouble(2) ** (nmant + 1) - 1
    assert_equal(as_int(v), 2**(nmant + 1) -1)
    # Check for predictable overflow
    nexp64 = floor_log2(type_info(np.float64)['max'])
    with np.errstate(over='ignore'):
        val = np.longdouble(2**nexp64) * 2  # outside float64 range
    assert_raises(OverflowError, as_int, val)
    assert_raises(OverflowError, as_int, -val)

    def compute_sentence_vector(self, sentence, model, w, vocab, zeros, word_list_to_exclude, settings):

        if len(sentence) == 1:
            tokens = sentence.split()  # tokenize by whitespace
        else:
            tokens = sentence  # sentence is already tokenized

        sentence_vector = np.longdouble(zeros)
        oov_count = 0

        for token in tokens:
            token = token.strip()

            try:
                if (settings['excludestopwords'] == '1') and token in word_list_to_exclude:
                    continue  # token somehow is not eligible to be used in our representation

                token_vector = self.get_word_vector(token, settings['method'], model, w, vocab)

            except Exception, e:
                # print(str(e))
                oov_count += 1
            else:
                token_vector = np.longdouble(np.asarray(token_vector))
                sentence_vector = sentence_vector + token_vector  # make vector summation for each token

 def __init__(self,):
     # Necessary parameters for all binary model
     self.binary_name = None
     self.param_default_value = {'PB':np.longdouble(10.0)*u.day,
                        'PBDOT':0.0*u.day/u.day,
                        'ECC': 0.9*u.Unit('') ,
                        'EDOT':0.0/u.second ,
                        'A1':10.0*ls,'A1DOT':0.0*ls/u.second,
                        'T0':np.longdouble(54000.0)*u.day,
                        'OM':10.0*u.deg,
                        'OMDOT':0.0*u.deg/u.year,
                        'XPBDOT':0.0*u.day/u.day,
                        'M2':0.0*u.M_sun,
                        'SINI':0*u.Unit(''),
                        'GAMMA':0*u.second, }
     # For Binary phase calculation
     self.param_default_value.update({'P0': 1.0*u.second,
                                      'P1': 0.0*u.second/u.second,
                                      'PEPOCH': np.longdouble(54000.0)*u.day
                                     })
     self.param_aliases = {'ECC':['E'],'EDOT':['ECCDOT'],
                           'A1DOT':['XDOT']}
     self.binary_params = list(self.param_default_value.keys())
     self.inter_vars = ['E','M','nu','ecc','omega','a1','TM2']
     self.binary_delay_funcs = []
     self.d_binarydelay_d_par_funcs = []

    def generate_affine_backtransformation(self):
        """ Generate synthetic examples and test them to determine transformation

        This is the key method!
        """
        if type(self.example) == FeatureVector:
            testsample = FeatureVector.replace_data(
                self.example, numpy.zeros(self.example.shape))
            self.offset = numpy.longdouble(self._execute(testsample))
            self.trafo = FeatureVector.replace_data(
                self.example, numpy.zeros(self.example.shape))
            for j in range(len(self.example.feature_names)):
                testsample = FeatureVector.replace_data(
                    self.example,
                    numpy.zeros(self.example.shape))
                testsample[0][j] = 1.0
                self.trafo[0][j] = \
                    numpy.longdouble(self._execute(testsample) - self.offset)
        elif type(self.example) == TimeSeries:
            testsample = TimeSeries.replace_data(
                self.example, numpy.zeros(self.example.shape))
            self.offset = numpy.longdouble(numpy.squeeze(
                self._execute(testsample)))
            self.trafo = TimeSeries.replace_data(
                self.example, numpy.zeros(self.example.shape))
            for i in range(self.example.shape[0]):
                for j in range(self.example.shape[1]):
                    testsample = TimeSeries.replace_data(
                        self.example, numpy.zeros_like(self.example))
                    testsample[i][j] = 1.0
                    self.trafo[i][j] = \
                        numpy.longdouble(numpy.squeeze(self._execute(testsample))
                                       - self.offset)

 def __init__(self):
     self.number = 10
     self.str = 'Test'
     self.list = [1,2,3]
     self.array = np.array([1,2,3])
     self.long_number = np.longdouble(1)
     self.long_array = np.longdouble([1,2,3])

def poincare_section_sets(datafile, colx, coly, colvx):
    datf = open(datafile, 'r')
    data = np.array([0, 0, 0], dtype=np.longdouble)
    data1 = np.array([0, 0, 0], dtype=np.longdouble)
    firstline = datf.readline()
    firstline = firstline.strip()
    columns = firstline.split()
    cols = [colx, coly, colvx]
    secpoint = [[], []]
    for i in range(3):
        data[i] = np.longdouble(columns[cols[i]])
    for line in datf:
        for i in range(3):
            data1[i] = data[i]
            line = line.strip()
            columns = line.split()
        for i in range(3):
            data[i] = np.longdouble(columns[cols[i]])
        if (data[1] > 0 and data1[1] < 0):
            secpoint[0].append((data1[0] + data[0]) / 2)
            secpoint[1].append((data1[2] + data[2]) / 2)
        elif (data[1] == 0 and data1[1] < 0):
            secpoint[0].append(data[0])
            secpoint[1].append(data[2])
    datf.close()
    return secpoint

    def d_phase_d_param(self, toas, delay, param):
        """ Return the derivative of phase with respect to the parameter.
        """
        # TODO need to do correct chain rule stuff wrt delay derivs, etc
        # Is it safe to assume that any param affecting delay only affects
        # phase indirectly (and vice-versa)??
        par = getattr(self, param)
        result = np.longdouble(np.zeros(len(toas))) * u.cycle/par.units
        param_phase_derivs = []
        phase_derivs = self.phase_deriv_funcs
        delay_derivs = self.delay_deriv_funcs
        if param in list(phase_derivs.keys()):
            for df in phase_derivs[param]:
                result += df(toas, param, delay).to(result.unit,
                            equivalencies=u.dimensionless_angles())
        else:
            # Apply chain rule for the parameters in the delay.
            # total_phase = Phase1(delay(param)) + Phase2(delay(param))
            # d_total_phase_d_param = d_Phase1/d_delay*d_delay/d_param +
            #                         d_Phase2/d_delay*d_delay/d_param
            #                       = (d_Phase1/d_delay + d_Phase2/d_delay) *
            #                         d_delay_d_param

            d_delay_d_p = self.d_delay_d_param(toas, param)
            dpdd_result = np.longdouble(np.zeros(len(toas))) * u.cycle/u.second
            for dpddf in self.d_phase_d_delay_funcs:
                dpdd_result += dpddf(toas, delay)
            result = dpdd_result * d_delay_d_p
        return result.to(result.unit, equivalencies=u.dimensionless_angles())

 def d_delay_d_DMs(self, toas, param_name, acc_delay=None): # NOTE we should have a better name for this.
     """Derivatives for constant DM
     """
     tbl = toas.table
     try:
         bfreq = self.barycentric_radio_freq(toas)
     except AttributeError:
         warn("Using topocentric frequency for dedispersion!")
         bfreq = tbl['freq']
     par = getattr(self, param_name)
     unit = par.units
     if param_name == 'DM':
         order = 0
     else:
         pn, idxf, idxv = split_prefixed_name(param_name)
         order = idxv
     dms = self.get_DM_terms()
     dm_terms = np.longdouble(np.zeros(len(dms)))
     dm_terms[order] = np.longdouble(1.0)
     if self.DMEPOCH.value is None:
         DMEPOCH = tbl['tdbld'][0]
     else:
         DMEPOCH = self.DMEPOCH.value
     dt = (tbl['tdbld'] - DMEPOCH) * u.day
     dt_value = (dt.to(u.yr)).value
     d_dm_d_dm_param = taylor_horner(dt_value, dm_terms)* (self.DM.units/par.units)
     return DMconst * d_dm_d_dm_param/ bfreq**2.0

def enn(x):
    enn = np.zeros_like(x, dtype=np.longdouble)
    enn[0] = np.longdouble(4.0 * w(0, d) * a0 ** 2)
    enn[1] = np.longdouble(4.0 * w(1, d) * a1 ** 2)
    for i in range(2, len(x)):
        enn[i] = np.longdouble(4.0 * w(1, d) * x[i] ** 2)
    return np.sum(enn, dtype=np.longdouble)

 def TestD_phase_d_toa(self):
     pint_d_phase_d_toa = self.modelB1855.d_phase_d_toa(self.toasB1855)
     mjd = np.array([np.longdouble(t.jd1 - ut.DJM0)+np.longdouble(t.jd2) for t in self.toasB1855.table['mjd']])
     tempo_d_phase_d_toa = self.plc.eval_spin_freq(mjd)
     diff = pint_d_phase_d_toa.value - tempo_d_phase_d_toa
     relative_diff = diff/tempo_d_phase_d_toa
     assert np.all(relative_diff < 1e-8), 'd_phae_d_toa test filed.'

 def test_ldouble_mapping(self):
     """ Test mapping for extended-precision """
     self.assertEqual(h5t.NATIVE_LDOUBLE.dtype, np.longdouble(1).dtype)
     if hasattr(np, 'float96'):
         self.assertEqual(h5t.py_create(np.dtype('float96')).dtype, np.longdouble(1).dtype)
     if hasattr(np, 'float128'):
         self.assertEqual(h5t.py_create(np.dtype('float128')).dtype, np.longdouble(1).dtype)

def ref_mjd(fits_file, hdu=1):
    """Read MJDREFF+ MJDREFI or, if failed, MJDREF, from the FITS header.

    Parameters
    ----------
    fits_file : str

    Returns
    -------
    mjdref : numpy.longdouble
        the reference MJD

    Other Parameters
    ----------------
    hdu : int
    """
    import collections

    if isinstance(fits_file, collections.Iterable) and\
            not is_string(fits_file):  # pragma: no cover
        fits_file = fits_file[0]
        logging.info("opening %s" % fits_file)

    try:
        ref_mjd_int = np.long(read_header_key(fits_file, 'MJDREFI'))
        ref_mjd_float = np.longdouble(read_header_key(fits_file, 'MJDREFF'))
        ref_mjd_val = ref_mjd_int + ref_mjd_float
    except:  # pragma: no cover
        ref_mjd_val = np.longdouble(read_header_key(fits_file, 'MJDREF'))
    return ref_mjd_val

	def __init__(self, delt_t, x=0.0, y=0.0, z=0.0,
				v_x=0.0, v_y=0.0, v_z=0.0, mass=0):
		"""Units:
		Position: km
		Volume: km/s
		Mass: kg
		Delta_t: days
		"""
		self.pos = np.array([np.longdouble(x),
							np.longdouble(y),
							np.longdouble(z)])

		self.vol = np.array([np.longdouble(v_x),
							np.longdouble(v_y),
							np.longdouble(v_z)])

		self.mass = mass

		self.del_t = np.longdouble(delt_t) * 86400

		# Initilize non input vars
		# Initilize history of positions
		self.hist = [[np.longdouble(x)],
					[np.longdouble(y)],
					[np.longdouble(z)]]
		# Initilize force array
		self.force = np.zeros(3, dtype=np.longdouble)
		# Initilize constants
		self.VELOCITY_FACTOR = self.del_t / self.mass / 1000
		self.FORCE_FACTOR = self.GRAVITATIONAL_CONSTANT * self.mass / 1000000

 def evalfreq(self,t):
     '''Return the freq at time t, computed with this polyco entry'''
     dt = (np.longdouble(t)-self.tmid.value)*np.longdouble(1440.0)
     s = np.longdouble(0.0)
     for i in range(1,self.ncoeff):
         s += np.longdouble(i) * self.coeffs[i] * dt**(i-1)
     freq = self.f0 + s/60.0
     return(freq)

 def __init__(self,toa,freq,error,site,flags=dict(),name="X"):
     e=toa.split(".")
     self.toa = (np.longdouble(e[0])+ np.longdouble("0."+e[1]))*u.day
     self.freq = freq.to(u.megahertz)
     self.site = site
     self.flags=flags
     self.error = error.to(u.microsecond)
     self.name=name

def gibb_roll(sides: int, biases: list):
    assert len(biases) == sides
    number = np.longdouble(random.uniform(0.0, np.sum(biases)))
    current = np.longdouble(0.0)
    for i, bias in enumerate(biases):
        current += bias
        if current >= number:
            return i

 def evalfreqderiv(self,t):
     """ Return the frequency derivative at time t."""
     dt = (np.longdouble(t)-self.tmid.value)*np.longdouble(1440.0)
     s = np.longdouble(0.0)
     for i in range(2,self.ncoeff):
         # Change to long double
         s += np.longdouble(i) * np.longdouble(i-1) * self.coeffs[i] * dt**(i-2)
     freqd = s/(60.0*60.0)
     return(freqd)

    def __init__(self, timfile):

        self.TOAs = []
        self.jump_parameters = []
        self.jump_labels = []
        self.TOA_uncertainties = []

        jump_count = 0
 
        for line in file(timfile):

            if ("MODE" in line):
                pass

            elif ("EFAC" in line):
                pass

            elif ("EQUAD" in line):
                pass

            elif ("FORMAT" in line):
                pass

            elif ("INCLUDE" in line):
                pass

            elif ("JUMP" in line):

                jump_count += 1

                if (jump_count % 2 != 0):
                    if (jump_count < 10):
                        self.jump_parameters.append('JUMP_000' + str(jump_count - (jump_count - 1) / 2))
                    if (jump_count >= 10 and jump_count < 100):
                        self.jump_parameters.append('JUMP_00' + str(jump_count - (jump_count - 1) / 2))

            elif ("INFO" in line):
                pass

            elif ("MODE" in line):
                pass

            else:
                
                elems = line.split()
                self.TOAs.append(np.longdouble(elems[0]))
                self.TOA_uncertainties.append(np.longdouble(elems[1]))
                
                if (jump_count % 2 == 0):
                    self.jump_labels.append('base')

                else:
                    self.jump_labels.append('JUMP_000' + str(jump_count - (jump_count - 1) / 2))

        self.TOAs = np.array(self.TOAs, dtype=np.longdouble)
        self.TOA_uncertainties = np.array(self.TOA_uncertainties, dtype=np.longdouble)