def test_evolve_constant_voltage(self):
     ptree = PropertyTree()
     ptree.put_string('mode', 'constant_voltage')
     ptree.put_double('voltage', 2.1)
     evolve_one_time_step = TimeEvolution.factory(ptree)
     evolve_one_time_step(device, 0.1)
     self.assertEqual(device.get_voltage(), 2.1)

    def test_builders(self):
        for AbstractClass in [Observer, Observable]:
            # AbstractClass takes a PropertyTree as argument.
            self.assertRaises(TypeError, AbstractClass)

            # The PropertyTree must specify what concrete class derived from
            # AbstractClass to instantiate.
            ptree = PropertyTree()
            self.assertRaises(RuntimeError, AbstractClass, ptree)

            # The derived concrete class must be registerd in the dictionary
            # that holds the builders.
            ptree.put_string('type', 'Invalid')
            self.assertRaises(KeyError, AbstractClass, ptree)

            # Now declare a concrete class.
            class ConcreteClass(AbstractClass):
                def __new__(cls, *args, **kwargs):
                    return object.__new__(ConcreteClass)
                def __init__(*args, **kwargs):
                    pass
            # Here is how to register a derived concrete class to the base abstract class.
            AbstractClass._builders['ConcreteClass'] = ConcreteClass

            # Now instantiation works.
            ptree.put_string('type', 'ConcreteClass')
            o = AbstractClass(ptree)

            # Also can build directly from derived class.
            o = ConcreteClass()

            # Remove from the dictionary.
            del AbstractClass._builders['ConcreteClass']
            self.assertRaises(KeyError, AbstractClass, ptree)

 def test_evolve_constant_load(self):
     ptree = PropertyTree()
     ptree.put_string('mode', 'constant_load')
     ptree.put_double('load', 120)
     evolve_one_time_step = TimeEvolution.factory(ptree)
     evolve_one_time_step(device, 0.1)
     self.assertAlmostEqual(device.get_voltage()/device.get_current(), -120)

    def test_abstract_class(self):
        # Declare a concrete Experiment
        class DummyExperiment(Experiment):

            def __new__(cls, *args, **kwargs):
                return object.__new__(DummyExperiment)

            def __init__(self, ptree):
                Experiment.__init__(self)
        # Do not forget to register it to the builders dictionary.
        Observable._builders['Dummy'] = DummyExperiment

        # Construct directly via DummyExperiment with a PropertyTree as a
        # positional arguemnt
        ptree = PropertyTree()
        dummy = DummyExperiment(ptree)

        # ... or directly via Experiment by specifying the ``type`` of
        # Experiment.
        ptree.put_string('type', 'Dummy')
        dummy = Experiment(ptree)

        # The method run() must be overloaded.
        self.assertRaises(RuntimeError, dummy.run, None)

        # Override the method run().
        def run(self, device):
            pass
        DummyExperiment.run = run

        # Now calling it without raising an error.
        dummy.run(None)

 def test_hold(self):
     ptree = PropertyTree()
     ptree.put_string('mode', 'hold')
     evolve_one_time_step = TimeEvolution.factory(ptree)
     device.evolve_one_time_step_constant_voltage(0.1, 1.4)
     evolve_one_time_step(device, 0.1)
     self.assertEqual(device.get_voltage(), 1.4)

 def test_evolve_constant_current(self):
     ptree = PropertyTree()
     ptree.put_string('mode', 'constant_current')
     ptree.put_double('current', 100e-3)
     evolve_one_time_step = TimeEvolution.factory(ptree)
     evolve_one_time_step(device, 0.1)
     self.assertEqual(device.get_current(), 100e-3)

 def test_time_limit(self):
     ptree = PropertyTree()
     ptree.put_string('end_criterion', 'time')
     ptree.put_double('duration', 15)
     time_limit = EndCriterion.factory(ptree)
     time_limit.reset(0.0, device)
     self.assertFalse(time_limit.check(2.0, device))
     self.assertTrue(time_limit.check(15.0, device))
     self.assertTrue(time_limit.check(60.0, device))

 def testParallelRC(self):
     # make parallel RC equivalent circuit
     device_database = PropertyTree()
     device_database.put_string('type', 'ParallelRC')
     device_database.put_double('series_resistance', R)
     device_database.put_double('parallel_resistance', R_L)
     device_database.put_double('capacitance', C)
     device = EnergyStorageDevice(device_database, MPI.COMM_WORLD)
     # setup experiment and measure
     eis_database = setup_expertiment()
     spectrum_data = measure_impedance_spectrum(device, eis_database)
     # extract data
     f = spectrum_data['frequency']
     Z_computed = spectrum_data['impedance']
     M_computed = 20*log10(absolute(Z_computed))
     P_computed = angle(Z_computed)*180/pi
     # compute the exact solution
     Z_exact = R+R_L/(1+1j*R_L*C*2*pi*f)
     M_exact = 20*log10(absolute(Z_exact))
     P_exact = angle(Z_exact)*180/pi
     # ensure the error is small
     max_phase_error_in_degree = linalg.norm(P_computed-P_exact, inf)
     max_magniture_error_in_decibel = linalg.norm(M_computed-M_exact, inf)
     print('max_phase_error_in_degree = {0}'.format(max_phase_error_in_degree))
     print('max_magniture_error_in_decibel = {0}'.format(max_magniture_error_in_decibel))
     self.assertLessEqual(max_phase_error_in_degree, 1)
     self.assertLessEqual(max_magniture_error_in_decibel, 0.2)

    def test_export_eclab_ascii_format(self):
        # define dummy experiment
        # it is quicker than building an actual EIS experiment
        class DummyExperiment(Experiment):

            def __new__(cls, *args, **kwargs):
                return object.__new__(DummyExperiment)

            def __init__(self, ptree):
                Experiment.__init__(self)
        dummy = DummyExperiment(PropertyTree())
        # produce dummy data for the experiment
        # here just a circle on the complex plane
        n = 10
        f = ones(n, dtype=float)
        Z = ones(n, dtype=complex)
        for i in range(n):
            f[i] = 10**(i / (n - 1))
            Z[i] = cos(2 * pi * i / (n - 1)) + 1j * sin(2 * pi * i / (n - 1))
        dummy._data['frequency'] = f
        dummy._data['impedance'] = Z
        # need a supercapacitor here to make sure method inspect() is kept in
        # sync with the EC-Lab headers
        ptree = PropertyTree()
        ptree.parse_info('super_capacitor.info')
        super_capacitor = EnergyStorageDevice(ptree)
        dummy._extra_data = super_capacitor.inspect()

        # export the data to ECLab format
        eclab = ECLabAsciiFile('untitled.mpt')
        eclab.update(dummy)

        # check that all lines end up with Windows-style line break '/r/n'
        # file need to be open in byte mode or the line ending will be
        # converted to '\n'...
        # also check that the number of lines in the headers has been computed
        # correctly and that the last one contains the column headers
        with open('untitled.mpt', mode='rb') as fin:
            lines = fin.readlines()
            for line in lines:
                self.assertNotEqual(line.find(b'\r\n'), -1)
                self.assertNotEqual(line.find(b'\r\n'), len(line) - 4)
            header_lines = int(lines[1].split(
                b':')[1].lstrip(b'').rstrip(b'\r\n'))
            self.assertEqual(
                header_lines,
                len(eclab._unformated_headers)
            )
            self.assertEqual(lines[header_lines - 1].find(b'freq/Hz'), 0)

        # check Nyquist plot does not throw
        nyquist = NyquistPlot('nyquist.png')
        nyquist.update(dummy)

 def test_constant_current_charge_for_given_time(self):
     ptree = PropertyTree()
     ptree.put_string('mode', 'constant_current')
     ptree.put_double('current', 5e-3)
     ptree.put_string('end_criterion', 'time')
     ptree.put_double('duration', 15.0)
     ptree.put_double('time_step', 0.1)
     stage = Stage(ptree)
     data = initialize_data()
     steps = stage.run(device, data)
     self.assertEqual(steps, 150)
     self.assertEqual(steps, len(data['time']))
     self.assertAlmostEqual(data['time'][-1], 15.0)
     self.assertAlmostEqual(data['current'][-1], 5e-3)

 def testRetrieveData(self):
     try:
         from h5py import File
     except ImportError:
         print('module h5py not found')
         return
     device_database = PropertyTree()
     device_database.put_string('type', 'SeriesRC')
     device_database.put_double('series_resistance', R)
     device_database.put_double('capacitance', C)
     device = EnergyStorageDevice(device_database, MPI.COMM_WORLD)
     eis_database  = setup_expertiment()
     eis_database.put_int('steps_per_decade', 1)
     eis_database.put_int('steps_per_cycle', 64)
     eis_database.put_int('cycles', 2)
     eis_database.put_int('ignore_cycles', 1)
     fout = File('trash.hdf5', 'w')
     spectrum_data = measure_impedance_spectrum(device, eis_database, fout)
     fout.close()
     fin = File('trash.hdf5', 'r')
     retrieved_data = retrieve_impedance_spectrum(fin)
     fin.close()
     print(spectrum_data['impedance']-retrieved_data['impedance'])
     print(retrieved_data)
     self.assertEqual(linalg.norm(spectrum_data['frequency'] -
                                  retrieved_data['frequency'], inf), 0.0)
     # not sure why we don't get equality for the impedance
     self.assertLess(linalg.norm(spectrum_data['impedance'] -
                                 retrieved_data['impedance'], inf), 1e-10)

 def test_setup_frequency_range(self):
     ptree = PropertyTree()
     ptree.put_string('type', 'ElectrochemicalImpedanceSpectroscopy')
     # specify the upper and lower bounds of the range
     # the number of points per decades controls the spacing on the log
     # scale
     ptree.put_double('frequency_upper_limit', 1e+2)
     ptree.put_double('frequency_lower_limit', 1e-1)
     ptree.put_int('steps_per_decade', 3)
     eis = Experiment(ptree)
     print(eis._frequencies)
     f = eis._frequencies
     self.assertEqual(len(f), 10)
     self.assertAlmostEqual(f[0], 1e+2)
     self.assertAlmostEqual(f[3], 1e+1)
     self.assertAlmostEqual(f[9], 1e-1)
     # or directly specify the frequencies
     frequencies = [3, 2e3, 0.1]
     eis = Experiment(ptree, frequencies)
     self.assertTrue(all(equal(frequencies, eis._frequencies)))

 def test_voltage_limit(self):
     ptree = PropertyTree()
     ptree.put_double('voltage_limit', 1.7)
     # upper limit
     ptree.put_string('end_criterion', 'voltage_greater_than')
     voltage_limit = EndCriterion.factory(ptree)
     voltage_limit.reset(5.0, device)
     device.evolve_one_time_step_constant_voltage(0.2, 1.3)
     self.assertFalse(voltage_limit.check(0.0, device))
     self.assertFalse(voltage_limit.check(60.0, device))
     device.evolve_one_time_step_constant_voltage(0.2, 1.7)
     self.assertTrue(voltage_limit.check(45.0, device))
     device.evolve_one_time_step_constant_voltage(0.2, 2.1)
     self.assertTrue(voltage_limit.check(45.0, device))
     # lower limit
     ptree.put_string('end_criterion', 'voltage_less_than')
     voltage_limit = EndCriterion.factory(ptree)
     voltage_limit.reset(0.0, device)
     device.evolve_one_time_step_constant_voltage(0.2, 1.3)
     self.assertTrue(voltage_limit.check(0.0, device))
     device.evolve_one_time_step_constant_voltage(0.2, 1.7)
     self.assertTrue(voltage_limit.check(45.0, device))
     device.evolve_one_time_step_constant_voltage(0.2, 2.1)
     self.assertFalse(voltage_limit.check(45.0, device))

 def test_fourier_analysis(self):
     ptree = PropertyTree()
     ptree.put_int('steps_per_cycle', 3)
     ptree.put_int('cycles', 1)
     ptree.put_int('ignore_cycles', 0)
     ptree.put_string('harmonics', '1')
     # uninitialized data
     data = {}
     self.assertRaises(KeyError, fourier_analysis, data, ptree)
     # empty data
     data = initialize_data()
     self.assertRaises(IndexError, fourier_analysis, data, ptree)
     # bad data
     data['time'] = array([1, 2, 3], dtype=float)
     data['current'] = array([4, 5, 6], dtype=float)
     data['voltage'] = array([7, 8], dtype=float)
     self.assertRaises(AssertionError, fourier_analysis, data, ptree)
     # poor data (size not a power of 2)
     data['voltage'] = array([7, 8, 9], dtype=float)
     with catch_warnings():
         simplefilter("error")
         self.assertRaises(RuntimeWarning, fourier_analysis, data, ptree)
     # data unchanged after analyze
     dummy = array([1, 2, 3, 4, 5, 6, 7, 8], dtype=float)
     data['time'] = dummy
     data['current'] = dummy
     data['voltage'] = dummy
     # ptree needs to be updated
     self.assertRaises(AssertionError, fourier_analysis, data, ptree)
     ptree.put_int('steps_per_cycle', 4)
     ptree.put_int('cycles', 2)
     ptree.put_int('ignore_cycles', 0)
     fourier_analysis(data, ptree)
     self.assertTrue(all(equal(data['time'], dummy)))
     self.assertTrue(all(equal(data['current'], dummy)))
     self.assertTrue(all(equal(data['voltage'], dummy)))

    def test_retrieve_data(self):
        ptree = PropertyTree()
        ptree.put_string('type', 'SeriesRC')
        ptree.put_double('series_resistance', 50e-3)
        ptree.put_double('capacitance', 3)
        device = EnergyStorageDevice(ptree)

        ptree = PropertyTree()
        ptree.put_string('type', 'RagoneAnalysis')
        ptree.put_double('discharge_power_lower_limit', 1e-1)
        ptree.put_double('discharge_power_upper_limit', 1e+1)
        ptree.put_int('steps_per_decade', 1)
        ptree.put_double('initial_voltage', 2.1)
        ptree.put_double('final_voltage', 0.7)
        ptree.put_double('time_step', 1.5)
        ptree.put_int('min_steps_per_discharge', 20)
        ptree.put_int('max_steps_per_discharge', 30)
        ragone = Experiment(ptree)

        with File('trash.hdf5', 'w') as fout:
            ragone.run(device, fout)
        performance_data = ragone._data

        fin = File('trash.hdf5', 'r')
        retrieved_data = retrieve_performance_data(fin)
        fin.close()
        # a few digits are lost when power is converted to string
        self.assertLess(linalg.norm(performance_data['power'] -
                                    retrieved_data['power'], inf), 1e-12)
        self.assertEqual(linalg.norm(performance_data['energy'] -
                                     retrieved_data['energy'], inf), 0.0)

        # TODO: probably want to move this into its own test
        ragoneplot = RagonePlot("ragone.png")
        ragoneplot.update(ragone)

        # check reset reinitialize the time step and empty the data
        ragone.reset()
        self.assertEqual(ragone._ptree.get_double('time_step'), 1.5)
        self.assertFalse(ragone._data['power'])
        self.assertFalse(ragone._data['energy'])

 def test_verification_with_equivalent_circuit(self):
     R = 50e-3   # ohm
     R_L = 500   # ohm
     C = 3       # farad
     # setup EIS experiment
     ptree = PropertyTree()
     ptree.put_string('type', 'ElectrochemicalImpedanceSpectroscopy')
     ptree.put_double('frequency_upper_limit', 1e+4)
     ptree.put_double('frequency_lower_limit', 1e-6)
     ptree.put_int('steps_per_decade', 3)
     ptree.put_int('steps_per_cycle', 1024)
     ptree.put_int('cycles', 2)
     ptree.put_int('ignore_cycles', 1)
     ptree.put_double('dc_voltage', 0)
     ptree.put_string('harmonics', '3')
     ptree.put_string('amplitudes', '5e-3')
     ptree.put_string('phases', '0')
     eis = Experiment(ptree)
     # setup equivalent circuit database
     device_database = PropertyTree()
     device_database.put_double('series_resistance', R)
     device_database.put_double('parallel_resistance', R_L)
     device_database.put_double('capacitance', C)
     # analytical solutions
     Z = {}
     Z['SeriesRC'] = lambda f: R + 1 / (1j * C * 2 * pi * f)
     Z['ParallelRC'] = lambda f: R + R_L / (1 + 1j * R_L * C * 2 * pi * f)
     for device_type in ['SeriesRC', 'ParallelRC']:
         # create a device
         device_database.put_string('type', device_type)
         device = EnergyStorageDevice(device_database)
         # setup experiment and measure
         eis.reset()
         eis.run(device)
         f = eis._data['frequency']
         Z_computed = eis._data['impedance']
         # compute the exact solution
         Z_exact = Z[device_type](f)
         # ensure the error is small
         max_phase_error_in_degree = linalg.norm(
             angle(Z_computed) * 180 / pi - angle(Z_exact) * 180 / pi,
             inf)
         max_magniture_error_in_decibel = linalg.norm(
             20 * log10(absolute(Z_exact)) - 20 *
             log10(absolute(Z_computed)),
             inf)
         print(device_type)
         print(
             '-- max_phase_error_in_degree = {0}'.format(max_phase_error_in_degree))
         print(
             '-- max_magniture_error_in_decibel = {0}'.format(max_magniture_error_in_decibel))
         self.assertLessEqual(max_phase_error_in_degree, 1)
         self.assertLessEqual(max_magniture_error_in_decibel, 0.2)

def run():
    # parse uq database
    input_database=PropertyTree()
    input_database.parse_xml('uq.xml')
    uq_database=input_database.get_child('uq')

    # declare parameters
    params=uq_database.get_int('params')
    parameter_list=[]
    for p in range(params):
        parameter_database=uq_database.get_child('param_'+str(p))
        distribution_type=parameter_database.get_string('distribution_type')
        parameter_name=parameter_database.get_string('name')
        if distribution_type=='uniform':
            parameter_range=parameter_database.get_array_double('range')
            parameter_list.append(UniformParameter('param_'+str(p),
                                                   parameter_name,
                                                   min=parameter_range[0],
                                                   max=parameter_range[1]))
        elif distribution_type=='normal':
            parameter_mean=parameter_database.get_double('mean')
            parameter_standard_deviation=parameter_database.get_double('standard_deviation')
            parameter_list.append(NormalParameter('param_'+str(p),
                                                  parameter_name,
                                                  mean=parameter_mean,
                                                  dev=parameter_standard_deviation))
        else:
            raise RuntimeError('invalid distribution type '+distribution_type+' for param_'+str(p))

    # create a host
    host_database=uq_database.get_child('host')
    host_type=host_database.get_string('type')
    if host_type=="Interactive":
        host=InteractiveHost(cpus=host_database.get_int_with_default_value('cpus',1),
                             cpus_per_node=host_database.get_int_with_default_value('cpus_per_node',0))
    elif host_type=="PBS":
        host=PBSHost(host_database.get_string('env'),
                     cpus=host_database.get_int_with_default_value('cpus',0),
                     cpus_per_node=host_database.get_int_with_default_value('cpus_per_node',0),
                     qname=host_database.get_string_with_default_value('qname','standby'),
                     walltime=host_database.get_string_with_default_value('walltime','1:00:00'),
                     modules=host_database.get_string_with_default_value('modules',''),
                     pack=host_database.get_int_with_default_value('pack',1),
                     qlimit=host_database.get_int_with_default_value('qlimit',200))
    else:
        raise RuntimeError('invalid host type '+host_type)

    # pick UQ method
    method=uq_database.get_string('method')
    if method=='SmolyakSparseGrid':
        level=uq_database.get_int('level')
        uq=Smolyak(parameter_list,level=level)
    elif method=='MonteCarlo':
        samples=uq_database.get_int('samples')
        uq=MonteCarlo(parameter_list,num=samples)
    elif method=='LatinHypercubeSampling':
        samples=uq_database.get_int('samples')
        uq=LHS(parameter_list,num=samples)
    else:
        raise RuntimeError('invalid UQ method '+method)

    # make a test program
    test_program_database=uq_database.get_child('test_program')
    description=test_program_database.get_string('description')
    executable_name=test_program_database.get_string('executable')
    for p in range(params):
        executable_name+=' --param_'+str(p)+' $param_'+str(p)
    prog=TestProgram(exe=executable_name,
                     desc=description)

    # run
    return Sweep(uq,host,prog)

 def test_charge_constant_voltage(self):
     ptree = PropertyTree()
     ptree.put_string('charge_mode', 'constant_voltage')
     ptree.put_double('charge_voltage', 1.4)
     ptree.put_string('charge_stop_at_1', 'current_less_than')
     ptree.put_double('charge_current_limit', 1e-6)
     ptree.put_string('charge_stop_at_2', 'time')
     ptree.put_double('charge_max_duration', 60)
     ptree.put_double('time_step', 0.2)
     charge = Charge(ptree)
     data = initialize_data()
     charge.run(device, data)
     self.assertTrue(data['time'][-1] >= 60 or
                     abs(data['current'][-1]) <= 1e-6)
     self.assertAlmostEqual(data['voltage'][-1], 1.4)

 def test_charge_constant_current(self):
     ptree = PropertyTree()
     ptree.put_string('charge_mode', 'constant_current')
     ptree.put_double('charge_current', 10e-3)
     ptree.put_string('charge_stop_at_1', 'voltage_greater_than')
     ptree.put_double('charge_voltage_limit', 1.4)
     ptree.put_double('time_step', 0.2)
     charge = Charge(ptree)
     data = initialize_data()
     charge.run(device, data)
     self.assertAlmostEqual(data['current'][0], 10e-3)
     self.assertAlmostEqual(data['current'][-1], 10e-3)
     self.assertGreaterEqual(data['voltage'][-1], 1.4)
     self.assertAlmostEqual(data['time'][1] - data['time'][0], 0.2)

 def test_time_steps(self):
     ptree = PropertyTree()
     ptree.put_int('stages', 2)
     ptree.put_int('cycles', 1)
     ptree.put_double('time_step', 1.0)
     ptree.put_string('stage_0.mode', 'hold')
     ptree.put_string('stage_0.end_criterion', 'time')
     ptree.put_double('stage_0.duration', 2.0)
     ptree.put_string('stage_1.mode', 'rest')
     ptree.put_string('stage_1.end_criterion', 'time')
     ptree.put_double('stage_1.duration', 1.0)
     ptree.put_double('stage_1.time_step', 0.1)
     multi = MultiStage(ptree)
     data = initialize_data()
     steps = multi.run(device, data)
     self.assertEqual(steps, 12)
     self.assertEqual(steps, len(data['time']))
     self.assertAlmostEqual(data['time'][5], 2.4)
     self.assertAlmostEqual(data['voltage'][0], data['voltage'][1])
     self.assertAlmostEqual(data['current'][3], 0.0)