Python scipy.optimize.brenth() Examples

def flash_TVF_zs(self, T, VF, zs):
        assert 0 <= VF <= 1
        Psats = self._Psats(T)
        # handle one component
        if self.N == 1:
            return 'l/g', [1.0], [1.0], VF, Psats[0]
        
        if VF == 0:
            P = bubble_at_T(zs, Psats)
        elif VF == 1:
            P = dew_at_T(zs, Psats)
        else:
            P = brenth(self._T_VF_err, min(Psats)*(1+1E-7), max(Psats)*(1-1E-7), args=(VF, zs, Psats))
        Ks = [K_value(P=P, Psat=Psat) for Psat in Psats]
        V_over_F, xs, ys = flash_inner_loop(zs=zs, Ks=Ks)
        return 'l/g', xs, ys, V_over_F, P 

def P_dew_at_T(self, T, zs, Psats=None):
        Psats = self._Psats(Psats, T)
        Pmax = self.P_bubble_at_T(T, zs, Psats)
        diff = 1E-7
        # EOSs do not solve at very low pressure
        if self.use_phis:
            Pmin = max(Pmax*diff, 1)
        else:
            Pmin = Pmax*diff
        P_dew = brenth(self._T_VF_err, Pmin, Pmax, args=(T, zs, Psats, Pmax, 1))
        self.__TVF_solve_cache = None
        return P_dew
#        try:
#            return brent(self._dew_P_UNIFAC_err, args=(T, zs, Psats, Pmax), brack=(Pmax*diff, Pmax*(1-diff), Pmax))
#        except:
#        return golden(self._dew_P_UNIFAC_err, args=(T, zs, Psats, Pmax), brack=(Pmax, Pmax*(1-diff)))
# 

def flash_TVF_zs(self, T, VF, zs):
        assert 0 <= VF <= 1
        Psats = self._Psats(T=T)
        Pbubble = self.P_bubble_at_T(T=T, zs=zs, Psats=Psats)
        if VF == 0:
            P = Pbubble
        else:
            diff = 1E-7
            Pmax = Pbubble
            # EOSs do not solve at very low pressure
            if self.use_phis:
                Pmin = max(Pmax*diff, 1)
            else:
                Pmin = Pmax*diff

            P = brenth(self._T_VF_err, Pmin, Pmax, args=(T, zs, Psats, Pmax, VF))
            self.__TVF_solve_cache = None
#            P = brenth(self._T_VF_err, Pdew, Pbubble, args=(T, VF, zs, Psats))
        V_over_F, xs, ys = self._flash_sequential_substitution_TP(T=T, P=P, zs=zs, Psats=Psats)
        return 'l/g', xs, ys, V_over_F, P 

def ci_corr(self, sig=.05, upper_bound=None, lower_bound=None):
        """
        Returns the confidence intervals for the correlation coefficient

        Parameters
        ----------
        sig : float
            The significance level.  Default is .05

        upper_bound : float
            Maximum value the upper confidence limit can be.
            Default is  99% confidence limit assuming normality.

        lower_bound : float
            Minimum value the lower condidence limit can be.
            Default is 99% confidence limit assuming normality.

        Returns
        -------
        interval : tuple
            Confidence interval for the correlation

        """
        endog = self.endog
        nobs = self.nobs
        self.r0 = chi2.ppf(1 - sig, 1)
        point_est = np.corrcoef(endog[:, 0], endog[:, 1])[0, 1]
        if upper_bound is None:
            upper_bound = min(.999, point_est + \
                          2.5 * ((1. - point_est ** 2.) / \
                          (nobs - 2.)) ** .5)

        if lower_bound is None:
            lower_bound = max(- .999, point_est - \
                          2.5 * (np.sqrt((1. - point_est ** 2.) / \
                          (nobs - 2.))))

        llim = optimize.brenth(self._ci_limits_corr, lower_bound, point_est)
        ulim = optimize.brenth(self._ci_limits_corr, point_est, upper_bound)
        return llim, ulim 

def ci_corr(self, sig=.05, upper_bound=None, lower_bound=None):
        """
        Returns the confidence intervals for the correlation coefficient

        Parameters
        ----------
        sig : float
            The significance level.  Default is .05

        upper_bound : float
            Maximum value the upper confidence limit can be.
            Default is  99% confidence limit assuming normality.

        lower_bound : float
            Minimum value the lower condidence limit can be.
            Default is 99% confidence limit assuming normality.

        Returns
        -------
        interval : tuple
            Confidence interval for the correlation

        """
        endog = self.endog
        nobs = self.nobs
        self.r0 = chi2.ppf(1 - sig, 1)
        point_est = np.corrcoef(endog[:, 0], endog[:, 1])[0, 1]
        if upper_bound is None:
            upper_bound = min(.999, point_est + \
                          2.5 * ((1. - point_est ** 2.) / \
                          (nobs - 2.)) ** .5)

        if lower_bound is None:
            lower_bound = max(- .999, point_est - \
                          2.5 * (np.sqrt((1. - point_est ** 2.) / \
                          (nobs - 2.))))

        llim = optimize.brenth(self._ci_limits_corr, lower_bound, point_est)
        ulim = optimize.brenth(self._ci_limits_corr, point_est, upper_bound)
        return llim, ulim 

def _Tsats(self, P):
        Tsats = []
        for i in self.VaporPressures:
            try: 
                Tsats.append(i.solve_prop(P))
            except:
                error = lambda T: i.extrapolate_tabular(T) - P
                Tsats.append(brenth(error, i.Tmax, i.Tmax*5))
        return Tsats 

def flash_PVF_zs(self, P, VF, zs):
        assert 0 <= VF <= 1
        Tsats = self._Tsats(P)
        # handle one component
        if self.N == 1:
            return 'l/g', [1.0], [1.0], VF, Tsats[0]
        
        T = brenth(self._P_VF_err, min(Tsats)*(1+1E-7), max(Tsats)*(1-1E-7), args=(P, VF, zs))
        Psats = self._Psats(T)
        Ks = [K_value(P=P, Psat=Psat) for Psat in Psats]
        V_over_F, xs, ys = flash_inner_loop(zs=zs, Ks=Ks)
        return 'l/g', xs, ys, V_over_F, T 

def flash_PVF_zs(self, P, VF, zs):
        try:
            # In some caases, will find a false root - resort to iterations which
            # are always between Pdew and Pbubble if this happens
            T = brenth(self._P_VF_err_2, min(self.Tms), min(self.Tcs), args=(P, VF, zs), maxiter=500)
            V_over_F, xs, ys = self._flash_sequential_substitution_TP(T=T, P=P, zs=zs)
            assert abs(V_over_F-VF) < 1E-6
        except:
            T = ridder(self._P_VF_err, min(self.Tms), min(self.Tcs), args=(P, VF, zs))
            V_over_F, xs, ys = self._flash_sequential_substitution_TP(T=T, P=P, zs=zs)
        return 'l/g', xs, ys, V_over_F, T 

def solve_prop(self, goal, reset_method=True):
        r'''Method to solve for the temperature at which a property is at a
        specified value. `T_dependent_property` is used to calculate the value
        of the property as a function of temperature; if `reset_method` is True,
        the best method is used at each temperature as the solver seeks a
        solution. This slows the solution moderately.

        Checks the given property value with `test_property_validity` first
        and raises an exception if it is not valid. Requires that Tmin and
        Tmax have been set to know what range to search within.

        Search is performed with the brenth solver from SciPy.

        Parameters
        ----------
        goal : float
            Propoerty value desired, [`units`]
        reset_method : bool
            Whether or not to reset the method as the solver searches

        Returns
        -------
        T : float
            Temperature at which the property is the specified value [K]
        '''
        if self.Tmin is None or self.Tmax is None:
            raise Exception('Both a minimum and a maximum value are not present indicating there is not enough data for temperature dependency.')
        if not self.test_property_validity(goal):
            raise Exception('Input property is not considered plausible; no method would calculate it.')

        def error(T):
            if reset_method:
                self.method = None
            return self.T_dependent_property(T) - goal
        try:
            return brenth(error, self.Tmin, self.Tmax)
        except ValueError:
            raise Exception('To within the implemented temperature range, it is not possible to calculate the desired value.') 

def find_steady_state(self, a, b, method='brentq', **kwargs):
        """
        Compute the equilibrium value of capital stock (per unit effective
        labor).

        Parameters
        ----------
        a : float
            One end of the bracketing interval [a,b].
        b : float
            The other end of the bracketing interval [a,b]
        method : str (default=`brentq`)
            Method to use when computing the steady state. Supported methods
            are `bisect`, `brenth`, `brentq`, `ridder`. See `scipy.optimize`
            for more details (including references).
        kwargs : optional
            Additional keyword arguments. Keyword arguments are method specific
            see `scipy.optimize` for details.

        Returns
        -------
        x0 : float
            Zero of `f` between `a` and `b`.
        r : RootResults (present if ``full_output = True``)
            Object containing information about the convergence. In particular,
            ``r.converged`` is True if the routine converged.

        """
        if method == 'bisect':
            result = optimize.bisect(self.evaluate_k_dot, a, b, **kwargs)
        elif method == 'brenth':
            result = optimize.brenth(self.evaluate_k_dot, a, b, **kwargs)
        elif method == 'brentq':
            result = optimize.brentq(self.evaluate_k_dot, a, b, **kwargs)
        elif method == 'ridder':
            result = optimize.ridder(self.evaluate_k_dot, a, b, **kwargs)
        else:
            mesg = ("Method must be one of : 'bisect', 'brenth', 'brentq', " +
                    "or 'ridder'.")
            raise ValueError(mesg)

        return result 

def find_steady_state(self, a, b, method='brentq', **kwargs):
        """
        Compute the equilibrium value of capital stock (per unit
        effective labor).

        Parameters
        ----------
        a : float
            One end of the bracketing interval [a,b].
        b : float
            The other end of the bracketing interval [a,b]
        method : str (default=`brentq`)
            Method to use when computing the steady state. Supported
            methods are `bisect`, `brenth`, `brentq`, `ridder`. See
            `scipy.optimize` for more details (including references).
        kwargs : optional
            Additional keyword arguments. Keyword arguments are method
            specific see `scipy.optimize` for details.

        Returns
        -------
        x0 : float
            Zero of `f` between `a` and `b`.
        r : RootResults (present if ``full_output = True``)
            Object containing information about the convergence. In
            particular, ``r.converged`` is True if the routine
            converged.

        """
        if method == 'bisect':
            result = optimize.bisect(self.evaluate_k_dot, a, b, **kwargs)
        elif method == 'brenth':
            result = optimize.brenth(self.evaluate_k_dot, a, b, **kwargs)
        elif method == 'brentq':
            result = optimize.brentq(self.evaluate_k_dot, a, b, **kwargs)
        elif method == 'ridder':
            result = optimize.ridder(self.evaluate_k_dot, a, b, **kwargs)
        else:
            mesg = ("Method must be one of : 'bisect', 'brenth', 'brentq', " +
                    "or 'ridder'.")
            raise ValueError(mesg)

        return result 

def flash_PH_zs_bounded(self, P, Hm, zs, T_low=None, T_high=None, 
                            Hm_low=None, Hm_high=None):
        # Begin the search at half the lowest chemical's melting point
        if T_low is None:
            T_low = min(self.Tms)/2 
                
        # Cap the T high search at 8x the highest critical point
        # (will not work well for helium, etc.)
        if T_high is None:
            max_Tc = max(self.Tcs)
            if max_Tc < 100:
                T_high = 4000
            else:
                T_high = max_Tc*8
    
        temp_pkg_cache = []
        def PH_error(T, P, zs, H_goal):
            if not temp_pkg_cache:
                temp_pkg = self.to(T=T, P=P, zs=zs)
                temp_pkg_cache.append(temp_pkg)
            else:
                temp_pkg = temp_pkg_cache[0]
                temp_pkg.flash(T=T, P=P, zs=zs)
            temp_pkg._post_flash()
            return temp_pkg.Hm - H_goal
    
        try:
            T_goal = brenth(PH_error, T_low, T_high, args=(P, zs, Hm))
            return T_goal

        except ValueError:
            if Hm_low is None:
                pkg_low = self.to(T=T_low, P=P, zs=zs)
                pkg_low._post_flash()
                Hm_low = pkg_low.Hm
            if Hm < Hm_low:
                raise ValueError('The requested molar enthalpy cannot be found'
                                 ' with this bounded solver because the lower '
                                 'temperature bound %g K has an enthalpy (%g '
                                 'J/mol) higher than that requested (%g J/mol)' %(
                                                             T_low, Hm_low, Hm))
            if Hm_high is None:
                pkg_high = self.to(T=T_high, P=P, zs=zs)
                pkg_high._post_flash()
                Hm_high = pkg_high.Hm
            if Hm > Hm_high:
                raise ValueError('The requested molar enthalpy cannot be found'
                                 ' with this bounded solver because the upper '
                                 'temperature bound %g K has an enthalpy (%g '
                                 'J/mol) lower than that requested (%g J/mol)' %(
                                                             T_high, Hm_high, Hm)) 

def flash_PS_zs_bounded(self, P, Sm, zs, T_low=None, T_high=None, 
                            Sm_low=None, Sm_high=None):
        # Begin the search at half the lowest chemical's melting point
        if T_low is None:
            T_low = min(self.Tms)/2 
                
        # Cap the T high search at 8x the highest critical point
        # (will not work well for helium, etc.)
        if T_high is None:
            max_Tc = max(self.Tcs)
            if max_Tc < 100:
                T_high = 4000
            else:
                T_high = max_Tc*8
    
        temp_pkg_cache = []
        def PS_error(T, P, zs, S_goal):
            if not temp_pkg_cache:
                temp_pkg = self.to(T=T, P=P, zs=zs)
                temp_pkg_cache.append(temp_pkg)
            else:
                temp_pkg = temp_pkg_cache[0]
                temp_pkg.flash(T=T, P=P, zs=zs)
            temp_pkg._post_flash()
            return temp_pkg.Sm - S_goal
    
        try:
            T_goal = brenth(PS_error, T_low, T_high, args=(P, zs, Sm))
            return T_goal

        except ValueError:
            if Sm_low is None:
                pkg_low = self.to(T=T_low, P=P, zs=zs)
                pkg_low._post_flash()
                Sm_low = pkg_low.Sm
            if Sm < Sm_low:
                raise ValueError('The requested molar entropy cannot be found'
                                 ' with this bounded solver because the lower '
                                 'temperature bound %g K has an entropy (%g '
                                 'J/mol/K) higher than that requested (%g J/mol/K)' %(
                                                             T_low, Sm_low, Sm))
            if Sm_high is None:
                pkg_high = self.to(T=T_high, P=P, zs=zs)
                pkg_high._post_flash()
                Sm_high = pkg_high.Sm
            if Sm > Sm_high:
                raise ValueError('The requested molar entropy cannot be found'
                                 ' with this bounded solver because the upper '
                                 'temperature bound %g K has an entropy (%g '
                                 'J/mol/K) lower than that requested (%g J/mol/K)' %(
                                                             T_high, Sm_high, Sm))