from __future__ import division
import unittest
import pytest
import numpy as np
from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials
from openmdao.api import IndepVarComp, Group, Problem
from openconcept.utilities.math.integrals import Integrator

class IntegratorTestGroup(Group):
    """An OpenMDAO group to test the every-node integrator component

    Options
    -------
    num_nodes : int
        Number of analysis points per segment
    quantity_units : str
        The units of the integral quantity q (NOT the rate)
    diff_units : str
        The units of the integrand (none by default)
    integrator : str
        Which integration scheme to use (default 'simpson')
    """

    def initialize(self):
        self.options.declare('quantity_units',default=None, desc="Units of the quantity being differentiated")
        self.options.declare('diff_units',default=None, desc="Units of the differential")
        self.options.declare('rate_units',default=None, desc="Units of the rate")
        self.options.declare('num_nodes',default=11, desc="Number of nodes per segment")
        self.options.declare('integrator',default='simpson', desc="Which simpson integrator to use")
        self.options.declare('time_setup',default='dt')
        self.options.declare('second_integrand',default=False)
        self.options.declare('zero_start', default=False)
        self.options.declare('final_only', default=False)
        self.options.declare('test_auto_names', default=False)
        self.options.declare('val', default=0.0)

    def setup(self):
        quantity_units = self.options['quantity_units']
        diff_units = self.options['diff_units']
        rate_units = self.options['rate_units']
        num_nodes = self.options['num_nodes']
        integrator_option = self.options['integrator']
        time_setup = self.options['time_setup']
        second_integrand = self.options['second_integrand']
        zero_start = self.options['zero_start']
        final_only = self.options['final_only']
        test_auto_names = self.options['test_auto_names']
        val = self.options['val']
        num_nodes = num_nodes

        iv = self.add_subsystem('iv', IndepVarComp())
        integ = Integrator(diff_units=diff_units, num_nodes=num_nodes, method=integrator_option, 
                              time_setup=time_setup)
        if not test_auto_names:
            integ.add_integrand('q', rate_name='dqdt', start_name='q_initial', end_name='q_final',
                                units=quantity_units, rate_units=rate_units, zero_start=zero_start,
                                final_only=final_only, val=val)
        else:
            integ.add_integrand('q', units=quantity_units, rate_units=rate_units, zero_start=zero_start,
                                final_only=final_only)
        if second_integrand:
            integ.add_integrand('q2', rate_name='dq2dt', start_name='q2_initial', end_name='q2_final', units='kJ')
            iv.add_output('rate_to_integrate_2', val=np.ones((num_nodes,)), units='kW')
            iv.add_output('initial_value_2', val=0., units='kJ')
            self.connect('iv.rate_to_integrate_2', 'integral.dq2dt')
            self.connect('iv.initial_value_2', 'integral.q2_initial')
        self.add_subsystem('integral', integ)

        if rate_units and quantity_units:
            # overdetermined and possibly inconsistent
            pass
        elif not rate_units and not quantity_units:
            if diff_units:
                rate_units =  '(' + diff_units +')** -1'
        elif not rate_units:
            # solve for rate_units in terms of quantity_units
            if not diff_units:
                rate_units = quantity_units
            else:
                rate_units = '('+quantity_units+') / (' + diff_units +')'
        elif not quantity_units:
            # solve for quantity units in terms of rate units
            if not diff_units:
                quantity_units = rate_units
            else:
                quantity_units = '('+rate_units+')*('+diff_units+')'

        iv.add_output('rate_to_integrate', val=np.ones((num_nodes,)), units=rate_units)
        iv.add_output('initial_value', val=0, units=quantity_units)

        if not test_auto_names:
            self.connect('iv.rate_to_integrate','integral.dqdt')
        else:
            self.connect('iv.rate_to_integrate','integral.q_rate')
        if not zero_start:
            self.connect('iv.initial_value', 'integral.q_initial')

        if time_setup == 'dt':
            iv.add_output('dt', val=1, units=diff_units)
            self.connect('iv.dt', 'integral.dt')
        elif time_setup == 'duration':
            iv.add_output('duration', val=1*(num_nodes-1), units=diff_units)
            self.connect('iv.duration', 'integral.duration')
        elif time_setup == 'bounds':
            iv.add_output('t_initial', val=2, units=diff_units)
            iv.add_output('t_final', val=2 + 1*(num_nodes-1), units=diff_units)
            self.connect('iv.t_initial','integral.t_initial')
            self.connect('iv.t_final','integral.t_final')

class IntegratorCommonTestCases(object):
    """
    A common set of test cases for the integrator component
    """

    def test_uniform_no_units(self):
        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator))
        prob.setup(check=True, force_alloc_complex=True)
        prob.run_model()
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        assert_near_equal(prob['integral.q'], np.linspace(0, nn_tot-1, nn_tot), tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units=None), nn_tot-1, tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_linear_no_units(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = x
        f = x ** 2 / 2

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob['integral.q'], f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_quadratic_no_units(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob['integral.q'], f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_machine_zero_rate(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 0.0 * x
        f = 0.0 * x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob['integral.q'], f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_auto_names(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(test_auto_names=True, num_nodes=self.num_nodes, integrator=self.integrator))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob['integral.q'], f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_quadratic_qty_units(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     quantity_units='kg', diff_units='s'))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_quadratic_qty_units_nonzero_start(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x + 25.2

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     quantity_units='kg', diff_units='s'))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob['iv.initial_value'] = 25.2
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_quadratic_qty_units_zero_start(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     quantity_units='kg', diff_units='s',
                                                     zero_start=True))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)
        with self.assertRaises(KeyError):
            # ensure the input hasn't been created
            prob['integral.q_initial'] = -1.0

    def test_quadratic_qty_units_final_only(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     quantity_units='kg', diff_units='s',
                                                     final_only=True))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        with self.assertRaises(KeyError):
            # this output shouldn't exist
            assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_quadratic_rate_units(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     rate_units='kg/s', diff_units='s'))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_quadratic_multiple_integrands(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime1 = 4 * x **2 - 8*x + 5
        f1 = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x
        fprime2 = -2 * x **2 + 10.5*x - 2
        f2 = -2 * x ** 3 / 3 + 10.5*x**2 /2 - 2*x
        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     rate_units='kg/s', diff_units='s', second_integrand=True))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime1
        prob['iv.rate_to_integrate_2'] = fprime2
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units='kg'), f1, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f1[-1], tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q2', units='kJ'), f2, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q2_final', units='kJ'), f2[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)


    def test_quadratic_both_units_correct(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x
        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                                rate_units='kg/s', quantity_units='kg', diff_units='s'))
        with self.assertRaises(ValueError) as cm:
            prob.setup(check=True)
            
        msg = ('Specify either quantity units or rate units, but not both')
        self.assertEqual(str(cm.exception), msg)

    def test_quadratic_duration(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     quantity_units='kg', diff_units='s',time_setup='duration'))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_quadratic_bounds(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     quantity_units='kg', diff_units='s',time_setup='bounds'))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_machine_zero_bounds(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 0.0 * x
        f = 0.0*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     quantity_units='kg', diff_units='s',time_setup='bounds'))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_quadratic_no_rate_units(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     diff_units='s'))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units=None), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

class SimpsonIntegrator5PtTestCases(unittest.TestCase, IntegratorCommonTestCases):
    """
    Only run the common test cases using second order accuracy because
    it cannot differentiate the quartic accurately
    """
    def __init__(self, *args, **kwargs):
        self.num_nodes = 11
        self.integrator = 'simpson'
        super(SimpsonIntegrator5PtTestCases, self).__init__(*args, **kwargs)

class SimpsonIntegrator3PtTestCases(unittest.TestCase, IntegratorCommonTestCases):
    """
    Only run the common test cases using second order accuracy because
    it cannot differentiate the quartic accurately
    """
    def __init__(self, *args, **kwargs):
        self.num_nodes = 7
        self.integrator = 'simpson'
        super(SimpsonIntegrator3PtTestCases, self).__init__(*args, **kwargs)

@pytest.mark.filterwarnings("ignore:divide by zero")
@pytest.mark.filterwarnings("ignore:invalid value")
class SimpsonIntegrator1PtTestCases(unittest.TestCase, IntegratorCommonTestCases):
    """
    Only run the common test cases using second order accuracy because
    it cannot differentiate the quartic accurately
    """
    def __init__(self, *args, **kwargs):
        self.num_nodes = 1
        self.integrator = 'simpson'
        super(SimpsonIntegrator1PtTestCases, self).__init__(*args, **kwargs)

class BDFIntegrator5PtTestCases(unittest.TestCase, IntegratorCommonTestCases):
    """
    Only run the common test cases using second order accuracy because
    it cannot differentiate the quartic accurately
    """
    def __init__(self, *args, **kwargs):
        self.num_nodes = 11
        self.integrator = 'bdf3'
        super(BDFIntegrator5PtTestCases, self).__init__(*args, **kwargs)

class BDFIntegrator3PtTestCases(unittest.TestCase, IntegratorCommonTestCases):
    """
    Only run the common test cases using second order accuracy because
    it cannot differentiate the quartic accurately
    """
    def __init__(self, *args, **kwargs):
        self.num_nodes = 7
        self.integrator = 'bdf3'
        super(BDFIntegrator3PtTestCases, self).__init__(*args, **kwargs)

class EdgeCaseTestCases(unittest.TestCase):
    def test_quadratic_even_num_nodes(self):
        num_nodes = 10
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x
        with self.assertRaises(ValueError) as cm:
            prob = Problem(IntegratorTestGroup(num_nodes=num_nodes, integrator='simpson',
                                                                    rate_units='kg/s', diff_units='s'))
            prob.setup()
            
        msg = ('num_nodes is ' +str(num_nodes) + ' and must be odd')
        self.assertEqual(str(cm.exception), msg)

    def test_default_value_scalar(self):
        num_nodes = self.num_nodes = 11
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator='simpson',
                                                     quantity_units='kg', diff_units='s',time_setup='duration',
                                                     val=5.0))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        # DO NOT RUN THE MODEL
        assert_near_equal(prob.get_val('integral.q', units='kg'), 5.0*np.ones((num_nodes,)), tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), 5.0, tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_default_value_vector(self):
        num_nodes = self.num_nodes = 11
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator='simpson',
                                                     quantity_units='kg', diff_units='s',time_setup='duration',
                                                     val=5.0*np.linspace(0.0, 1.0, num_nodes)))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        # DO NOT RUN THE MODEL
        assert_near_equal(prob.get_val('integral.q', units='kg'), 5.0*np.linspace(0.0, 1.0, num_nodes), tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), 5.0, tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

from __future__ import division
import unittest
import pytest
import numpy as np
from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials
from openmdao.api import IndepVarComp, Group, Problem
from openconcept.utilities.math.integrals import Integrator

class IntegratorTestGroup(Group):
    """An OpenMDAO group to test the every-node integrator component

    Options
    -------
    num_nodes : int
        Number of analysis points per segment
    quantity_units : str
        The units of the integral quantity q (NOT the rate)
    diff_units : str
        The units of the integrand (none by default)
    integrator : str
        Which integration scheme to use (default 'simpson')
    """

    def initialize(self):
        self.options.declare('quantity_units',default=None, desc="Units of the quantity being differentiated")
        self.options.declare('diff_units',default=None, desc="Units of the differential")
        self.options.declare('rate_units',default=None, desc="Units of the rate")
        self.options.declare('num_nodes',default=11, desc="Number of nodes per segment")
        self.options.declare('integrator',default='simpson', desc="Which simpson integrator to use")
        self.options.declare('time_setup',default='dt')
        self.options.declare('second_integrand',default=False)
        self.options.declare('zero_start', default=False)
        self.options.declare('final_only', default=False)
        self.options.declare('test_auto_names', default=False)
        self.options.declare('val', default=0.0)

    def setup(self):
        quantity_units = self.options['quantity_units']
        diff_units = self.options['diff_units']
        rate_units = self.options['rate_units']
        num_nodes = self.options['num_nodes']
        integrator_option = self.options['integrator']
        time_setup = self.options['time_setup']
        second_integrand = self.options['second_integrand']
        zero_start = self.options['zero_start']
        final_only = self.options['final_only']
        test_auto_names = self.options['test_auto_names']
        val = self.options['val']
        num_nodes = num_nodes

        iv = self.add_subsystem('iv', IndepVarComp())
        integ = Integrator(diff_units=diff_units, num_nodes=num_nodes, method=integrator_option, 
                              time_setup=time_setup)
        if not test_auto_names:
            integ.add_integrand('q', rate_name='dqdt', start_name='q_initial', end_name='q_final',
                                units=quantity_units, rate_units=rate_units, zero_start=zero_start,
                                final_only=final_only, val=val)
        else:
            integ.add_integrand('q', units=quantity_units, rate_units=rate_units, zero_start=zero_start,
                                final_only=final_only)
        if second_integrand:
            integ.add_integrand('q2', rate_name='dq2dt', start_name='q2_initial', end_name='q2_final', units='kJ')
            iv.add_output('rate_to_integrate_2', val=np.ones((num_nodes,)), units='kW')
            iv.add_output('initial_value_2', val=0., units='kJ')
            self.connect('iv.rate_to_integrate_2', 'integral.dq2dt')
            self.connect('iv.initial_value_2', 'integral.q2_initial')
        self.add_subsystem('integral', integ)

        if rate_units and quantity_units:
            # overdetermined and possibly inconsistent
            pass
        elif not rate_units and not quantity_units:
            if diff_units:
                rate_units =  '(' + diff_units +')** -1'
        elif not rate_units:
            # solve for rate_units in terms of quantity_units
            if not diff_units:
                rate_units = quantity_units
            else:
                rate_units = '('+quantity_units+') / (' + diff_units +')'
        elif not quantity_units:
            # solve for quantity units in terms of rate units
            if not diff_units:
                quantity_units = rate_units
            else:
                quantity_units = '('+rate_units+')*('+diff_units+')'

        iv.add_output('rate_to_integrate', val=np.ones((num_nodes,)), units=rate_units)
        iv.add_output('initial_value', val=0, units=quantity_units)

        if not test_auto_names:
            self.connect('iv.rate_to_integrate','integral.dqdt')
        else:
            self.connect('iv.rate_to_integrate','integral.q_rate')
        if not zero_start:
            self.connect('iv.initial_value', 'integral.q_initial')

        if time_setup == 'dt':
            iv.add_output('dt', val=1, units=diff_units)
            self.connect('iv.dt', 'integral.dt')
        elif time_setup == 'duration':
            iv.add_output('duration', val=1*(num_nodes-1), units=diff_units)
            self.connect('iv.duration', 'integral.duration')
        elif time_setup == 'bounds':
            iv.add_output('t_initial', val=2, units=diff_units)
            iv.add_output('t_final', val=2 + 1*(num_nodes-1), units=diff_units)
            self.connect('iv.t_initial','integral.t_initial')
            self.connect('iv.t_final','integral.t_final')

class IntegratorCommonTestCases(object):
    """
    A common set of test cases for the integrator component
    """

    def test_uniform_no_units(self):
        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator))
        prob.setup(check=True, force_alloc_complex=True)
        prob.run_model()
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        assert_near_equal(prob['integral.q'], np.linspace(0, nn_tot-1, nn_tot), tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units=None), nn_tot-1, tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_linear_no_units(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = x
        f = x ** 2 / 2

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob['integral.q'], f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_quadratic_no_units(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob['integral.q'], f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_machine_zero_rate(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 0.0 * x
        f = 0.0 * x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob['integral.q'], f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_auto_names(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(test_auto_names=True, num_nodes=self.num_nodes, integrator=self.integrator))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob['integral.q'], f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_quadratic_qty_units(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     quantity_units='kg', diff_units='s'))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_quadratic_qty_units_nonzero_start(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x + 25.2

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     quantity_units='kg', diff_units='s'))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob['iv.initial_value'] = 25.2
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_quadratic_qty_units_zero_start(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     quantity_units='kg', diff_units='s',
                                                     zero_start=True))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)
        with self.assertRaises(KeyError):
            # ensure the input hasn't been created
            prob['integral.q_initial'] = -1.0

    def test_quadratic_qty_units_final_only(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     quantity_units='kg', diff_units='s',
                                                     final_only=True))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        with self.assertRaises(KeyError):
            # this output shouldn't exist
            assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_quadratic_rate_units(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     rate_units='kg/s', diff_units='s'))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_quadratic_multiple_integrands(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime1 = 4 * x **2 - 8*x + 5
        f1 = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x
        fprime2 = -2 * x **2 + 10.5*x - 2
        f2 = -2 * x ** 3 / 3 + 10.5*x**2 /2 - 2*x
        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     rate_units='kg/s', diff_units='s', second_integrand=True))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime1
        prob['iv.rate_to_integrate_2'] = fprime2
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units='kg'), f1, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f1[-1], tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q2', units='kJ'), f2, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q2_final', units='kJ'), f2[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)


    def test_quadratic_both_units_correct(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x
        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                                rate_units='kg/s', quantity_units='kg', diff_units='s'))
        with self.assertRaises(ValueError) as cm:
            prob.setup(check=True)
            
        msg = ('Specify either quantity units or rate units, but not both')
        self.assertEqual(str(cm.exception), msg)

    def test_quadratic_duration(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     quantity_units='kg', diff_units='s',time_setup='duration'))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_quadratic_bounds(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     quantity_units='kg', diff_units='s',time_setup='bounds'))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_machine_zero_bounds(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 0.0 * x
        f = 0.0*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     quantity_units='kg', diff_units='s',time_setup='bounds'))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_quadratic_no_rate_units(self):
        num_nodes = self.num_nodes
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator,
                                                     diff_units='s'))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        prob.run_model()
        assert_near_equal(prob.get_val('integral.q', units=None), f, tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

class SimpsonIntegrator5PtTestCases(unittest.TestCase, IntegratorCommonTestCases):
    """
    Only run the common test cases using second order accuracy because
    it cannot differentiate the quartic accurately
    """
    def __init__(self, *args, **kwargs):
        self.num_nodes = 11
        self.integrator = 'simpson'
        super(SimpsonIntegrator5PtTestCases, self).__init__(*args, **kwargs)

class SimpsonIntegrator3PtTestCases(unittest.TestCase, IntegratorCommonTestCases):
    """
    Only run the common test cases using second order accuracy because
    it cannot differentiate the quartic accurately
    """
    def __init__(self, *args, **kwargs):
        self.num_nodes = 7
        self.integrator = 'simpson'
        super(SimpsonIntegrator3PtTestCases, self).__init__(*args, **kwargs)

@pytest.mark.filterwarnings("ignore:divide by zero")
@pytest.mark.filterwarnings("ignore:invalid value")
class SimpsonIntegrator1PtTestCases(unittest.TestCase, IntegratorCommonTestCases):
    """
    Only run the common test cases using second order accuracy because
    it cannot differentiate the quartic accurately
    """
    def __init__(self, *args, **kwargs):
        self.num_nodes = 1
        self.integrator = 'simpson'
        super(SimpsonIntegrator1PtTestCases, self).__init__(*args, **kwargs)

class BDFIntegrator5PtTestCases(unittest.TestCase, IntegratorCommonTestCases):
    """
    Only run the common test cases using second order accuracy because
    it cannot differentiate the quartic accurately
    """
    def __init__(self, *args, **kwargs):
        self.num_nodes = 11
        self.integrator = 'bdf3'
        super(BDFIntegrator5PtTestCases, self).__init__(*args, **kwargs)

class BDFIntegrator3PtTestCases(unittest.TestCase, IntegratorCommonTestCases):
    """
    Only run the common test cases using second order accuracy because
    it cannot differentiate the quartic accurately
    """
    def __init__(self, *args, **kwargs):
        self.num_nodes = 7
        self.integrator = 'bdf3'
        super(BDFIntegrator3PtTestCases, self).__init__(*args, **kwargs)

class EdgeCaseTestCases(unittest.TestCase):
    def test_quadratic_even_num_nodes(self):
        num_nodes = 10
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x
        with self.assertRaises(ValueError) as cm:
            prob = Problem(IntegratorTestGroup(num_nodes=num_nodes, integrator='simpson',
                                                                    rate_units='kg/s', diff_units='s'))
            prob.setup()
            
        msg = ('num_nodes is ' +str(num_nodes) + ' and must be odd')
        self.assertEqual(str(cm.exception), msg)

    def test_default_value_scalar(self):
        num_nodes = self.num_nodes = 11
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator='simpson',
                                                     quantity_units='kg', diff_units='s',time_setup='duration',
                                                     val=5.0))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        # DO NOT RUN THE MODEL
        assert_near_equal(prob.get_val('integral.q', units='kg'), 5.0*np.ones((num_nodes,)), tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), 5.0, tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

    def test_default_value_vector(self):
        num_nodes = self.num_nodes = 11
        nn_tot = num_nodes
        x = np.linspace(0, nn_tot-1, nn_tot)
        fprime = 4 * x **2 - 8*x + 5
        f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x

        prob = Problem(IntegratorTestGroup(num_nodes=self.num_nodes, integrator='simpson',
                                                     quantity_units='kg', diff_units='s',time_setup='duration',
                                                     val=5.0*np.linspace(0.0, 1.0, num_nodes)))
        prob.setup(check=True, force_alloc_complex=True)
        prob['iv.rate_to_integrate'] = fprime
        # DO NOT RUN THE MODEL
        assert_near_equal(prob.get_val('integral.q', units='kg'), 5.0*np.linspace(0.0, 1.0, num_nodes), tolerance=1e-14)
        assert_near_equal(prob.get_val('integral.q_final', units='kg'), 5.0, tolerance=1e-14)
        partials = prob.check_partials(method='cs',compact_print=True)
        assert_check_partials(partials, atol=1e-8, rtol=1e0)

An OpenMDAO group to test the every-node integrator component      Options     -------     num_nodes : int         Number of analysis points per segment     quantity_units : str         The units of the integral quantity q (NOT the rate)     diff_units : str         The units of the integrand (none by default)     integrator : str         Which integration scheme to use (default 'simpson') # overdetermined and possibly inconsistent # solve for rate_units in terms of quantity_units # solve for quantity units in terms of rate units A common set of test cases for the integrator component # ensure the input hasn't been created # this output shouldn't exist Only run the common test cases using second order accuracy because     it cannot differentiate the quartic accurately Only run the common test cases using second order accuracy because     it cannot differentiate the quartic accurately Only run the common test cases using second order accuracy because     it cannot differentiate the quartic accurately Only run the common test cases using second order accuracy because     it cannot differentiate the quartic accurately Only run the common test cases using second order accuracy because     it cannot differentiate the quartic accurately # DO NOT RUN THE MODEL # DO NOT RUN THE MODEL