Python scipy.optimize.fsolve() Examples

def _fitstart(self, data):
        g1 = _skew(data)
        g2 = _kurtosis(data)

        def func(x):
            a, b = x
            sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)
            ku = a**3 - a**2*(2*b-1) + b**2*(b+1) - 2*a*b*(b+2)
            ku /= a*b*(a+b+2)*(a+b+3)
            ku *= 6
            return [sk-g1, ku-g2]
        a, b = optimize.fsolve(func, (1.0, 1.0))
        return super(beta_gen, self)._fitstart(data, args=(a, b)) 

def _xinf_ND(xdot,x0,args=(),xddot=None,xtol=1.49012e-8):
    """Private function for wrapping the fsolving for x_infinity
    for a variable x in N dimensions"""
    try:
        result = fsolve(xdot,x0,args,fprime=xddot,xtol=xtol,full_output=1)
    except (ValueError, TypeError, OverflowError):
        xinf_val = NaN
    except:
        print "Error in fsolve:", sys.exc_info()[0], sys.exc_info()[1]
        xinf_val = NaN
    else:
        if result[2] in (1,2,3): #,4,5):
            # 4,5 means "not making good progress" (see fsolve docstring)
            xinf_val = float(result[0])
        else:
            xinf_val = NaN
    return xinf_val 

def calculate_psi_goal(psi_baseline, Au_baseline, Au_goal, deltaz,
                       c_baseline, c_goal):
    """Find the value for psi that has the same location w on the upper
    surface of the goal as psi_baseline on the upper surface of the
    baseline"""

    def integrand(psi_baseline, Au, deltaz, c):
        return c*np.sqrt(1 + dxi_u(psi_baseline, Au, deltaz/c)**2)

    def equation(psi_goal, L_baseline, Au_goal, deltaz, c):
        y, err = quad(integrand, 0, psi_goal, args=(Au_goal, deltaz, c))
        return y - L_baseline

    L_baseline, err = quad(integrand, 0, psi_baseline,
                           args=(Au_baseline, deltaz, c_baseline))
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        y = fsolve(equation, psi_baseline, args=(L_baseline, Au_goal, deltaz,
                                                 c_goal))
    return y[0] 

def computeShiftedGears(m, alpha, t1, t2, x1, x2):
    """Summary

    Args:
        m (float): common module of both gears [length]
        alpha (float): pressure-angle [rad]
        t1 (int): number of teeth of gear1
        t2 (int): number of teeth of gear2
        x1 (float): relative profile-shift of gear1
        x2 (float): relative profile-shift of gear2

    Returns:
        (float, float): distance between gears [length], pressure angle of the assembly [rad]
    """
    def inv(x): return np.tan(x) - x
    inv_alpha_w = inv(alpha) + 2 * np.tan(alpha) * (x1 + x2) / (t1 + t2)

    def root_inv(x): return inv(x) - inv_alpha_w
    alpha_w = opt.fsolve(root_inv, 0.)
    dist = m * (t1 + t2) / 2 * np.cos(alpha) / np.cos(alpha_w)
    return dist, alpha_w 

def _fitstart(self, data):
        g1 = _skew(data)
        g2 = _kurtosis(data)

        def func(x):
            a, b = x
            sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)
            ku = a**3 - a**2*(2*b-1) + b**2*(b+1) - 2*a*b*(b+2)
            ku /= a*b*(a+b+2)*(a+b+3)
            ku *= 6
            return [sk-g1, ku-g2]
        a, b = optimize.fsolve(func, (1.0, 1.0))
        return super(beta_gen, self)._fitstart(data, args=(a, b)) 

def Z_Hall_Yarborough(Tr, Pr):
    """Calculate gas compressibility factor using the correlation of Hall &
    Yarborough (1973)

    Parameters
    ------------
    Tr : float
        Reduced temperature [-]
    Pr : float
        Reduced pressure [-]

    Returns
    -------
    Z : float
        Gas compressibility factor [-]
    """
    # Check input in range of validity
    if Tr <= 1:
        warnings.warn("Using extrapolated values")

    X1 = -0.06125*Pr/Tr*exp(-1.2*(1-1/Tr)**2)
    X2 = 14.76/Tr - 9.76/Tr**2 + 4.58/Tr**3
    X3 = 90.7/Tr - 242.2/Tr**2 + 42.4/Tr**3
    X4 = 2.18+2.82/Tr

    def f(Y):
        return X1+(Y+Y**2+Y**3-Y**4)/(1-Y)**3-X2*Y**2+X3*Y**X4

    Yo = 0.0125*Pr/Tr*exp(-1.2*(1-1/Tr)**2)
    Y = fsolve(f, Yo, full_output=True)
    if Y[2] == 1 and abs(Y[1]["fvec"][0]) < 1e-5:
        Z = 0.06125*Pr/Tr/Y[0][0]*exp(-1.2*(1-1/Tr)**2)
    else:
        raise ValueError("Iteration not converge")

    return unidades.Dimensionless(Z) 

def test_solve_equilibrium_2():
    c = np.array([1.7e-03, 3.0e+06, 3.0e+06, 9.7e+07, 5.55e+09])
    stoich = (1, 1, 0, 0, -1)
    K = 55*1e-6

    def f(x):
        return prodpow(c+x*stoich, stoich) - K
    solution = solve_equilibrium(c, stoich, K)
    assert np.allclose(solution, c + stoich*fsolve(f, 0.1)) 

def test_solve_equilibrium_1():
    c = np.array((13., 11, 17))
    stoich = np.array((-2, 3, -4))
    K = 3.14

    def f(x):
        return (13 - 2*x)**-2 * (
            11 + 3*x)**3 * (17 - 4*x)**-4 - K
    assert np.allclose(solve_equilibrium(c, stoich, K),
                       c + stoich*fsolve(f, 3.48)) 

def _schedule(self, rho, v):
        if rho < 1e-100:
            epsilon = 0
        else:
            fun = lambda epsilon: pow(epsilon, 2) + v * pow(epsilon, 3 / 2) - pow(rho, 2)
            epsilon = optimize.fsolve(fun, rho / 2)

        return (epsilon) 

def _saturation(self, T):
        """Saturation calculation for two phase search"""
        rhoc = self._constants.get("rhoref", self.rhoc)
        Tc = self._constants.get("Tref", self.Tc)

        if T > Tc:
            T = Tc
        tau = Tc/T

        rhoLo = self._Liquid_Density(T)
        rhoGo = self._Vapor_Density(T)

        def f(parr):
            rhol, rhog = parr
            deltaL = rhol/rhoc
            deltaG = rhog/rhoc
            phirL = _phir(tau, deltaL, self._constants)
            phirG = _phir(tau, deltaG, self._constants)
            phirdL = _phird(tau, deltaL, self._constants)
            phirdG = _phird(tau, deltaG, self._constants)
            Jl = deltaL*(1+deltaL*phirdL)
            Jv = deltaG*(1+deltaG*phirdG)
            Kl = deltaL*phirdL+phirL+log(deltaL)
            Kv = deltaG*phirdG+phirG+log(deltaG)
            return Kv-Kl, Jv-Jl

        rhoL, rhoG = fsolve(f, [rhoLo, rhoGo])
        if rhoL == rhoG:
            Ps = self.Pc
        else:
            deltaL = rhoL/self.rhoc
            deltaG = rhoG/self.rhoc
            firL = _phir(tau, deltaL, self._constants)
            firG = _phir(tau, deltaG, self._constants)

            Ps = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(firL-firG+log(deltaL/deltaG))
        return rhoL, rhoG, Ps 

def _OsmoticPressure(T, P, S):
    """Procedure to calculate the osmotic pressure of seawater

    Parameters
    ----------
    T : float
        Tmperature, [K]
    P : float
        Pressure, [MPa]
    S : float
        Salinity, [kg/kg]

    Returns
    -------
    Posm : float
        Osmotic pressure, [MPa]

    References
    ----------
    IAPWS,  Advisory Note No. 5: Industrial Calculation of the Thermodynamic
    Properties of Seawater, http://www.iapws.org/relguide/Advise5.html, Eq 15
    """
    pw = _Region1(T, P)
    gw = pw["h"]-T*pw["s"]

    def f(Posm):
        pw2 = _Region1(T, P+Posm)
        gw2 = pw2["h"]-T*pw2["s"]
        ps = SeaWater._saline(T, P+Posm, S)
        return -ps["g"]+S*ps["gs"]-gw+gw2

    Posm = fsolve(f, 0)[0]
    return Posm 

def _Triple(S):
    """Procedure to calculate the triple point pressure and temperature for
    seawater

    Parameters
    ----------
    S : float
        Salinity, [kg/kg]

    Returns
    -------
    prop : dict
        Dictionary with the triple point properties:

            * Tt: Triple point temperature, [K]
            * Pt: Triple point pressure, [MPa]

    References
    ----------
    IAPWS,  Advisory Note No. 5: Industrial Calculation of the Thermodynamic
    Properties of Seawater, http://www.iapws.org/relguide/Advise5.html, Eq 7
    """
    def f(parr):
        T, P = parr
        pw = _Region1(T, P)
        gw = pw["h"]-T*pw["s"]

        pv = _Region2(T, P)
        gv = pv["h"]-T*pv["s"]

        gih = _Ice(T, P)["g"]
        ps = SeaWater._saline(T, P, S)

        return -ps["g"]+S*ps["gs"]-gw+gih, -ps["g"]+S*ps["gs"]-gw+gv

    Tt, Pt = fsolve(f, [273, 6e-4])

    prop = {}
    prop["Tt"] = Tt
    prop["Pt"] = Pt
    return prop 

def _frank(dependency, val):
    if(dependency=='kendall'):
        return fsolve(_frank_kendall_fopt, 1, args=(val))[0]
    elif(dependency=='spearman'):
        # TODO --  use function solvers in scipy to invert debye function for the closed form solution
        raise NotImplementedError('Spearmans Rho currently unsupported for Frank Copula family!')
    
    return r 

def equilibrium_pressure(minerals, stoichiometry, temperature, pressure_initial_guess=1.e5):
    """
    Given a list of minerals, their reaction stoichiometries
    and a temperature of interest, compute the
    equilibrium pressure of the reaction.

    Parameters
    ----------
    minerals : list of minerals
        List of minerals involved in the reaction.

    stoichiometry : list of floats
        Reaction stoichiometry for the minerals provided.
        Reactants and products should have the opposite signs [mol]

    temperature : float
        Temperature of interest [K]

    pressure_initial_guess : optional float
        Initial pressure guess [Pa]

    Returns
    -------
    pressure : float
        The equilibrium pressure of the reaction [Pa]
    """
    def eqm(P, T):
        gibbs = 0.
        for i, mineral in enumerate(minerals):
            mineral.set_state(P[0], T)
            gibbs = gibbs + mineral.gibbs * stoichiometry[i]
        return gibbs

    pressure = fsolve(eqm, [pressure_initial_guess], args=(temperature))[0]

    return pressure 

def nullcline_x(x): # nullcline_x

	def func(Ca, x):
		return x_inf(params['V_Ca']-Ca/params['K_c']/x)-x

	Ca = list(opt.fsolve(func, [0.17], args=(x[0],)))
	for i in xrange(1, x.size):
		Ca.append(opt.fsolve(func, [Ca[-1]], args=(x[i],))[0])

	return Ca 

def g_critical(I):

        k, E, m = params['k_0'], params['E_0'], params['m_0']

        def g_of_V(V):
                expfct = np.exp(-k*(V-V_0))
                return m*(1.-3.*V**2)-k*expfct/(1.+expfct)**2
        
        def to_minimize(V):
                return m*(V-V**3)+I+g_of_V(V)*(E-V)-nullcline_y(V, V_0, k)

        V_opt = opt.fsolve(to_minimize, x0=-1.)

        return g_of_V(V_opt) 

def Z_Elsharkawy(Tr, Pr):
    """Calculate gas compressibility factor using the Elsharkawy (2003)
    correlation

    Parameters
    ------------
    Tr : float
        Reduced temperature [-]
    Pr : float
        Reduced pressure [-]

    Returns
    -------
    Z : float
        Gas compressibility factor [-]
    """
    C1 = 0.3265 - 1.07/Tr - 0.5339/Tr**3 + 0.01569/Tr**4 - 0.05165/Tr**5
    C2 = 0.5475 - 0.7361/Tr + 0.1844/Tr**2
    C3 = 0.1056*(0.7361/Tr + 0.1844/Tr**2)

    # Eq 24
    def f(rho):
        C4 = 0.6134*(1+0.721*rho**2)*rho**2/Tr**3*3
        Z = 0.27*Pr/Tr/rho
        return 1 + C1*rho + C2*rho**2 - C3*rho**5 + C4*exp(-0.721*rho**2) - Z

    rho0 = 0.27*Pr/Tr
    rho = fsolve(f, rho0, full_output=True)
    if rho[2] == 1:
        Z = 0.27*Pr/Tr/rho[0]
    else:
        raise ValueError("Iteration not converge")

    return unidades.Dimensionless(Z) 

def _Tf(P, S):
    """Procedure to calculate the freezing temperature of seawater

    Parameters
    ----------
    P : float
        Pressure, [MPa]
    S : float
        Salinity, [kg/kg]

    Returns
    -------
    Tf : float
        Freezing temperature, [K]

    References
    ----------
    IAPWS,  Advisory Note No. 5: Industrial Calculation of the Thermodynamic
    Properties of Seawater, http://www.iapws.org/relguide/Advise5.html, Eq 12
    """
    def f(T):
        T = float(T)
        pw = _Region1(T, P)
        gw = pw["h"]-T*pw["s"]

        gih = _Ice(T, P)["g"]

        ps = SeaWater._saline(T, P, S)
        return -ps["g"]+S*ps["gs"]-gw+gih

    Tf = fsolve(f, 300)[0]
    return Tf 

def Z_Hall_Iglesias(Tr, Pr):
    """Calculate gas compressibility factor using the correlation of Hall-
    Iglesias (2007). This is a extension of Hall-Yarborough correlation to
    reduced values of temperature.

    Parameters
    ------------
    Tr : float
        Reduced temperature [-]
    Pr : float
        Reduced pressure [-]

    Returns
    -------
    Z : float
        Gas compressibility factor [-]
    """
    X1 = 14.54/Tr-8.23/Tr**2+3.39*Tr**3.5
    X2 = 90.7/Tr-242.2/Tr**2+42.4/Tr**3
    X3 = 1.18+2.82/Tr
    k1 = 1/(1.87+0.001*Tr**63)
    k2 = 171.8*Tr**13
    k3 = -3525*(1-exp(-219/Tr**63))

    def f(Z):
        Y = (0.06125*Pr*exp(-1.2*(1-1/Tr)**2))/Tr/Z
        return Z - (1+Y+Y**2-Y**3)/(1-Y)**3 - X1*Y + X2*Y**X3 + \
            k1*Y*exp(-k2*(Y-0.421)**2) + k3*Y**10*exp(-69279*(Y-0.374)**4)

    Z = fsolve(f, 0.5, full_output=True)
    if Z[2] != 1:
        raise ValueError("Iteration not converge")

    return unidades.Dimensionless(Z[0]) 

def main(self, solver_type='fsolve', iterations=10):

        # initial velocity profile
        self.boundary_conditions()

        for self.nx in range(1, self.gsimax):
            self.solution, self.solver_message = \
                self.solver(solver_type, iterations)
            self.shift_profiles() 

def _twb(self, P, phi_w, tdb):
        return None

    # def _V_Virial(self, T, P, phi_w):
#        """volumen por unidad de masa de aire seco"""
#        # FIXME: Don't work, for now use ideal gas equation
#        Bm, Cm=self.Virial(T, Xa)
#        vm=roots([1, -R_atml*T/self.P.atm, -R_atml*T*Bm/self.P.atm,-R_atml*T*Cm/self.P.atm])
#        if vm[0].imag==0.0:
#            v=vm[0].real
#        else:
#            v=vm[2].real
#        return unidades.SpecificVolume(v/self.aire.M/Xa)
#
#    def _h(self, Td, w):
#        """Enthalpy calculation procedure"""
#        cp_air = self.Air.Cp_Gas_DIPPR(Td)
#        cp_water = self.Water.Cp_Gas_DIPPR(Td)
#        h = (Td-273.15)*(cp_air+cp_water*w)
#        return unidades.Enthalpy(h, "kJkg")
#
#    def _HS(self, Td):
#        """Saturation humidity calculation procedure"""
#        pv = self._Ps(Td)
#        return self.Water.M*pv/(self.Air.M*(self.P-pv))
#
#    def _Ps(self, Td):
#        """Saturation pressure calculation procedure"""
#        if Td < 273.15:
#            return _Sublimation_Pressure(Td)
#        else:
#            return _PSat_T(Td)
#
#    def _Tw(self, Td, w):
#        Hv = self.agua.Hv_DIPPR(Td)
#        Cs = self.Calor_Especifico_Humedo(Td, w)
#
#        def f(Tw):
#            return self.Humedad_Absoluta(Tw)-w-Cs/Hv*(Td-Tw)
#        Tw = fsolve(f, Td)
#        return unidades.Temperature(Tw) 

def _tdp(self, P, phi_w, tdb):
        """Calculation of dew-point temperature"""
        if not phi_w:
            return None

        def f(t):
            f = self._f(t, P)
            pws = self._Pwsat(t)
            pw = phi_w*P
            return pw-f*pws

        tdp = fsolve(f, tdb)
        return tdp 

def _twb(tdb, W, P):
    """
    Calculate the wet-bulb temperature of a humid air as a perfet gas, as
    explain in [1]_, pag 1.9, inveted Eq 35-37

    Parameters
    ------------
    Tdb: float
        Dry bulb temperature, [K]
    W : float
       Humidity ratio, [kgw/kgda]
    P : float
        Pressure, [Pa]

    Returns
    -------
    Twb: float
        Wet bulb temperature, [K]
    """
    tdb_C = tdb - 273.15

    def f(twb):
        Pvs = _Psat(twb)
        Ws = 0.62198*Pvs/(P-Pvs)
        twb_C = twb - 273.15
        if tdb >= 0:
            w = ((2501.-2.326*twb_C)*Ws-1.006*(tdb_C-twb_C)) / \
                (2501.+1.86*tdb_C-4.186*twb_C)-W
        else:
            w = ((2830-0.24*twb_C)*Ws-1.006*(tdb-twb)) / \
                (2830+1.86*tdb_C-2.1*twb_C)
        return w

    twb = fsolve(f, tdb)[0]
    return twb


# Procedures only used in plot 

def calculo(self):
        entrada = self.kwargs["entrada"]
        self.Tout = unidades.Temperature(self.kwargs["Tout"])
        self.deltaP = unidades.DeltaP(self.kwargs["deltaP"])
        self.Hmax = unidades.Power(self.kwargs["Hmax"])
        if self.kwargs["eficiencia"]:
            eficiencia = self.kwargs["eficiencia"]
        else:
            eficiencia = 0.75

        if self.kwargs["poderCalorifico"]:
            poderCalorifico = self.kwargs["poderCalorifico"]
        else:
            poderCalorifico = 900

        salida = entrada.clone(T=self.Tout, P=entrada.P-self.deltaP)
        Ho = entrada.h
        H1 = salida.h
        Heat = unidades.Power(H1-Ho)

        if self.Hmax and Heat > self.Hmax:
            self.Heat = unidades.Power(self.Hmax)
            To = (entrada.T+self.Tout)/2
            T = fsolve(lambda T: entrada.clone(
                T=T, P=entrada.P-self.deltaP).h-Ho-self.Hmax, To)[0]
            self.salida = [entrada.clone(T=T, P=entrada.P-self.deltaP)]
        else:
            self.Heat = Heat
            self.salida = [salida]

        fuel = self.Heat.Btuh/poderCalorifico/eficiencia
        self.CombustibleRequerido = unidades.VolFlow(fuel, "ft3h")
        self.deltaT = unidades.DeltaT(self.salida[0].T-entrada.T)
        self.eficiencia = unidades.Dimensionless(eficiencia)
        self.poderCalorifico = unidades.Dimensionless(poderCalorifico)
        self.Tin = entrada.T
        self.Tout = self.salida[0].T 

def calculo(self):
        entrada = self.kwargs["entrada"]
        self.deltaP = unidades.DeltaP(self.kwargs["DeltaP"])
        self.HeatCalc = unidades.Power(self.kwargs["Heat"])
        if self.kwargs["Tout"]:
            Tout = unidades.Temperature(self.kwargs["Tout"])
        elif self.kwargs["DeltaT"]:
            Tout = unidades.Temperature(entrada.T+self.kwargs["DeltaT"])
        A = unidades.Area(self.kwargs["A"])
        U = unidades.HeatTransfCoef(self.kwargs["U"])
        Text = unidades.Temperature(self.kwargs["Text"])

        if self.modo == 1:
            self.salida = [entrada.clone(T=Tout, P=entrada.P-self.deltaP)]
            self.HeatCalc = unidades.Power(self.salida[0].h-entrada.h)
        else:
            if self.modo == 2:
                self.HeatCalc = unidades.Power(0)
            else:
                self.HeatCalc = unidades.Power(A*U*(Text-entrada.T))

            def f():
                output = entrada.clone(T=T, P=entrada.P-self.deltaP)
                return output.h-entrada.h-self.HeatCalc
            T = fsolve(f, entrada.T)[0]
            if T > max(Text, entrada.T) or T < min(Text, entrada.T):
                T = self.Text
            self.salida = [entrada.clone(T=T, P=entrada.P-self.deltaP)]

        self.Tin = entrada.T
        self.ToutCalc = self.salida[0].T
        self.deltaT = unidades.DeltaT(self.ToutCalc-entrada.T) 

def calculo(self):
        self.areaCalculada = Area(self.kwargs["area"])
        self.deltaP = Pressure(self.kwargs["deltaP"])

        if self.kwargs["potencialCarga"]:
            self.potencialCarga = PotencialElectric(
                self.kwargs["potencialCarga"])
        else:
            self.potencialCarga = PotencialElectric(24000)
        if self.kwargs["potencialDescarga"]:
            self.potencialDescarga = PotencialElectric(
                self.kwargs["potencialDescarga"])
        else:
            self.potencialDescarga = PotencialElectric(24000)
        if self.kwargs["epsilon"]:
            self.epsilon = Dimensionless(self.kwargs["epsilon"])
        else:
            self.epsilon = Dimensionless(4.)

        if self.kwargs["metodo"] == 1:
            def f(area):
                eta_i = self.calcularRendimientos_parciales(area)
                eta = self.calcularRendimiento(eta_i)
                return eta-self.kwargs["rendimientoAdmisible"]

            self.areaCalculada = Area(fsolve(f, 100)[0])

        eta_i = self.calcularRendimientos_parciales(self.areaCalculada)
        self.rendimiento_parcial = eta_i
        self.rendimiento = self.calcularRendimiento(eta_i)

        self.CalcularSalidas() 

def getTemperatureAtDensity(self, targetDensity, temperatureGuessInC):
        """Get the temperature at which the perturbed density occurs."""
        densFunc = (
            lambda temp: self.density(Tc=temp) - targetDensity
        )  # 0 at tempertature of targetDensity
        tAtTargetDensity = float(
            fsolve(densFunc, temperatureGuessInC)
        )  # is a numpy array if fsolve is called
        return tAtTargetDensity 

def _digammainv(y):
    # Inverse of the digamma function (real positive arguments only).
    # This function is used in the `fit` method of `gamma_gen`.
    # The function uses either optimize.fsolve or optimize.newton
    # to solve `sc.digamma(x) - y = 0`.  There is probably room for
    # improvement, but currently it works over a wide range of y:
    #    >>> y = 64*np.random.randn(1000000)
    #    >>> y.min(), y.max()
    #    (-311.43592651416662, 351.77388222276869)
    #    x = [_digammainv(t) for t in y]
    #    np.abs(sc.digamma(x) - y).max()
    #    1.1368683772161603e-13
    #
    _em = 0.5772156649015328606065120
    func = lambda x: sc.digamma(x) - y
    if y > -0.125:
        x0 = np.exp(y) + 0.5
        if y < 10:
            # Some experimentation shows that newton reliably converges
            # must faster than fsolve in this y range.  For larger y,
            # newton sometimes fails to converge.
            value = optimize.newton(func, x0, tol=1e-10)
            return value
    elif y > -3:
        x0 = np.exp(y/2.332) + 0.08661
    else:
        x0 = 1.0 / (-y - _em)

    value, info, ier, mesg = optimize.fsolve(func, x0, xtol=1e-11,
                                             full_output=True)
    if ier != 1:
        raise RuntimeError("_digammainv: fsolve failed, y = %r" % y)

    return value[0]


## Gamma (Use MATLAB and MATHEMATICA (b=theta=scale, a=alpha=shape) definition)

## gamma(a, loc, scale)  with a an integer is the Erlang distribution
## gamma(1, loc, scale)  is the Exponential distribution
## gamma(df/2, 0, 2) is the chi2 distribution with df degrees of freedom. 

def _fitstart(self, data):
        g1 = _skew(data)
        g2 = _kurtosis(data)

        def func(x):
            a, b = x
            sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)
            ku = a**3 - a**2*(2*b-1) + b**2*(b+1) - 2*a*b*(b+2)
            ku /= a*b*(a+b+2)*(a+b+3)
            ku *= 6
            return [sk-g1, ku-g2]
        a, b = optimize.fsolve(func, (1.0, 1.0))
        return super(beta_gen, self)._fitstart(data, args=(a, b)) 

def test_float32(self):
        func = lambda x: np.array([x[0] - 100, x[1] - 1000], dtype=np.float32)**2
        p = optimize.fsolve(func, np.array([1, 1], np.float32))
        assert_allclose(func(p), [0, 0], atol=1e-3)