Coverage for chemreac/util/analysis : 100%

# -*- coding: utf-8 -*- 

""" 

chemreac.util.analysis 

---------------------- 

Functions to analyze numierc output from integration. 

""" 

from __future__ import (absolute_import, division, 

                        print_function, unicode_literals) 

import numpy as np 

def solver_linear_error(y, rtol, atol, logy=False, scale_err=1.0): 

    """ 

    Returns linear estimated error bounds from numerical integration 

    Parameters 

    ========== 

    y : array_like 

         Output from integration (before back-transformation is applied) 

    rtol : float 

         Relative tolerance 

    atol : float 

         Absolute tolerance 

    logy : bool 

         Is y from a run with logarithmic concentration? 

    scale_err : float 

         scale estimated error bounds (useful for testing) 

    Returns 

    ======= 

    Array of shape (2, len(y)) with rows corresponding to lower and 

    upper bounds around y. 

    .. note:: Assumes maximum mangitude of error be: \ 

    :math:`\\boldsymbol{e}_{max} = \|\\boldsymbol{y} \ 

    \\cdot \\mathrm{rtol}\| + \\mathrm{atol}` 

    """ 

    solver_err = scale_err*(np.abs(y*rtol) + atol) 

    if logy: 

        res = np.exp(y - solver_err), np.exp(y + solver_err) 

    else: 

        res = y - solver_err, y + solver_err 

    return np.array(res) 

def solver_linear_error_from_integration(integration, ti=slice(None), bi=0, 

                                         si=0, **kwargs): 

    """ 

    Convenience function wrapping :func:`solver_linear_error` 

    Parameters 

    ---------- 

    integration: Integration instance 

    """ 

    try: 

        atol_i = integration.info['atol'][si] 

    except TypeError: 

        atol_i = integration.info['atol'] 

    return integration.with_units( 

        solver_linear_error( 

            integration.yout[ti, bi, si], 

            integration.info['rtol'], 

            atol_i, 

            integration.rd.logy, 

            **kwargs), 

        'concentration') 

def suggest_t0(rd, y0, max_f=1.0): 

    """ 

    Suggests an appropriate initial time, 

    useful when logy==True and logt==True, 

    If suggested t0 > 1, 1 is returned. 

    Parameters 

    ========== 

    rd: ReactionDiffusion instance 

         System at hand 

    y0: sequence 

         initial values 

    max_f: float 

         upper bound of absolute value for largest element in for the 

         inital step. 

    """ 

    fout = rd.alloc_fout() 

    rd.f(0, y0, fout) 

    fout_maxabs = np.max(np.abs(fout)) 

    if fout_maxabs < max_f: 

        return 1.0 

    else: 

        return max_f/fout_maxabs