class Simulator(object):

    # List of all registered simulation modules. Note this are classes and we
    # need to instantiate them in __init__().
    _simulation_modules = []

    # Named tuple to be used in the recorder lambda function.
    _RecorderContext = namedtuple('RecorderContext', 'value time delta_time')

    @staticmethod
    @accepts((1, types.ClassType))
    def register_simulation_module(module_class):
        """
        register_simulation_module(module_class)

        Registers a simulation module class within the main simulator class.
        Note that you register the class, the simulator automatically
        instantiates an object of the class in its own constructor every time an
        instance of the simulator is created.

        :param module_class: The class to register as a simulation module.
        :type module_class: :class:`gridsim.core.AbstractSimulationModule`
        """
        Simulator._simulation_modules.append(module_class)

    @accepts((1, ExecutionManager))
    def __init__(self, execution_manager=DefaultExecutionManager()):
        """
        __init__(self)

        Gridsim main simulation class. Moves the simulation of all modules
        on in time and groups all simulation modules. The constructor creates
        automatically an instance of each registered simulation module. The fact
        to import a simulation module to the user's project registers the module
        automatically within the main simulation class. So only modules really
        used by the simulation are instantiated.
        """
        super(Simulator, self).__init__()

        self._execution_manager = execution_manager

        # Create an instance for each simulation module.
        self._modules = {}
        for module_class in Simulator._simulation_modules:
            module = module_class()
            module.simulator = self
            self._modules[module.attribute_name()] = module

        # Prepare an array for all recorders.
        self._time_recorders = []
        self._recorders = []
        self._recorderBindings = []

        # Initialize time to nothing.
        self.time = None

    @accepts((1, str))
    @returns(AbstractSimulationModule)
    def __getattr__(self, item):
        """
        __getattr__(self, item)

        If the main simulation object does not own an attribute with the given
        name <sim>.<name>, this method gets called. If there exists a simulation
        module with the given attribute name, we return the reference to that
        module. This enables the user to access very easily the different
        aspects of the simulation.

        :param item: The name of the attribute that does not exists per se which
            can be a potential simulation module.
        :type item: str
        
        :returns: :class:`AbstractSimulationModule` - The module corresponding to the
            item or None of no such module was found.
        """
        if item in self._modules.keys():
            return self._modules[item]
        else:
            raise AttributeError("The module "+str(item)+" is not loaded")

    @accepts(((1, 3, 6), (str, type(None))),((2, 4, 5), (int, type(None))),(7, (tuple, type(None))))
    @returns(list)
    def find(self, module=None, uid=None, friendly_name=None,
             element_class=None, instance_of=None,
             has_attribute=None, close_to=None):
        """
        find(self, module=None, uid=None, friendly_name=None, element_class=None, instance_of=None, has_attribute=None, close_to=None)

        Finds all :class:`.AbstractSimulationElement` derived
        objects matching the given criteria by searching on either the given
        Gridsim simulation module or by searching the whole simulation of the
        module was not specified. Note that the method returns always a list of
        element, even if only a single instance is found. All parameters are
        optional, if :func:`find()` will be called without any parameters,
        the list of all element in the actual simulation will be returned.

        :param module: The module to search for element.
        :type module: str
        
        :param uid: ID of the element. Note that these ID's are only unique for
            a given class inside a module, so you should never search just for
            an ID without specifying either the class or at least the module.
        :type uid: int
        
        :param friendly_name: The friendly name to search for.
        :type friendly_name: str
        
        :param element_class: The exact class the element has to be an instance 
            of. Superclasses are not considered.
        :type element_class: type
        
        :param instance_of: The object has to be an instance of the given class,
            whereas this can be the superclass of the object too.
        :type instance_of: type
        
        :param has_attribute: The object should have an attribute with the given
            name. This can be used in order to find all objects that have a 
            ``power`` attribute.
        :type has_attribute: str
        
        :param close_to: The object's position should be closer to the given one
            than the given radius. The parameters should be passed as tuple of 
            the position and the radius in meters [m]
        :type close_to: (Position, float)
        
        :return: List of :class:`.AbstractSimulationElement`
            matching the given criteria.

        *Example:*

        .. literalinclude:: ../../demo/find.py

        """
        elements = []

        # Get the complete list of element from either specific or all modules.
        if module is not None:
            if module in self._modules.keys():
                elements = self._modules[module].all_elements()
            else:
                return elements
        else:
            for module in self._modules.values():
                elements.extend(module.all_elements())

        # Only AbstractSimulationElements will be considered!
        elements = filter(
            lambda element: isinstance(element, AbstractSimulationElement),
            elements)

        # Apply the different filters.
        if uid is not None:
            elements = filter(lambda element: element.id == uid, elements)

        if friendly_name is not None:
            elements = filter(
                lambda element: element.friendly_name == friendly_name,
                elements)

        if element_class is not None:
            elements = filter(
                lambda element: element.__class__ == element_class, elements)

        if instance_of is not None:
            elements = filter(lambda element: isinstance(element, instance_of),
                              elements)

        if has_attribute is not None:
            elements = filter(lambda element: hasattr(element, has_attribute),
                              elements)

        if close_to is not None and len(close_to) is 2:
            position, radius = close_to
            if isinstance(position, Position) and isinstance(radius,
                                                             (float, int)):
                elements = filter(lambda element: element.position.distance_to(
                    position) <= radius, elements)

        return elements

    @units.wraps(None, (None, units.second))
    def reset(self, initial_time=0):
        """
        reset(self, initial_time=0*units.second)

        Resets (re-initializes) the simulation.

        :param initial_time: The initial time from which the simulation starts.
            Defaults to ``0``.
        :type initial_time: time, see :mod:`gridsim.unit`
        """
        self.time = initial_time

        for module in self._modules.values():
            module.reset()

        for recorder_binding in self._recorderBindings:
            recorder_binding.reset()

    @accepts((1, (int, float)))
    def _calculate(self, delta_time):
        for module in self._modules.values():
            module.calculate(self.time, delta_time)

    @accepts((1, (int, float)))
    def _update(self, delta_time):

        for module in self._modules.values():
            module.update(self.time, delta_time)

        for recorder in self._recorders:
            recorder.on_simulation_step(self.time)
        for recorder in self._time_recorders:
            recorder.on_simulation_step(self.time)

        for recorder_binding in self._recorderBindings:
            recorder_binding.update(self.time, delta_time)

    @accepts((1, (int, float)))
    def _step(self, delta_time):
        """
        _step(self, delta_time)

        Executes a single simulation step on all modules.

        :param delta_time: The delta time for the single step.
        :type delta_time: time, see :mod:`gridsim.unit`
        """
        delta_time_sec = units.value(delta_time, units.second)
        self._calculate(delta_time_sec)
        self.time += delta_time_sec
        self._update(delta_time_sec)

    def _end(self):

        """
        _end(self)

        Inform all module that the simulation ended

        """

        for module in self._modules.values():
            module.end(self.time)

    @units.wraps(None, (None, units.second))
    def step(self, delta_time):
        """
        step(self, delta_time)

        Executes a single simulation step on all modules.

        :param delta_time: The delta time for the single step.
        :type delta_time: time, see :mod:`gridsim.unit`
        """
        self._step(delta_time)

    @units.wraps(None, (None, units.second, units.second))
    def run(self, run_time, delta_time):
        """
        run(self, run_time, delta_time)

        Runs the simulation for a given time.

        The run_time parameter defines the duration the simulation has to run.
        If a simulation was already run before, the additional duration is
        run in addition.

        :param run_time: Total run time.
        :type run_time: time, see :mod:`gridsim.unit`
        
        :param delta_time: Time interval for the simulation.
        :type delta_time: time, see :mod:`gridsim.unit`
        """
        if self.time is None:
            self.reset()
        self._execution_manager.reset()

        self._execution_manager.preprocess()
        end_time = self.time + run_time
        self._update(delta_time)
        self._execution_manager.postprocess()

        while self.time < end_time:
            self._execution_manager.preprocess()
            self._step(delta_time)
            self._execution_manager.postprocess()

        if True: # from a cyberphysicalsystem need to end with a calculate 'Read'
            self._execution_manager.preprocess()
            self._calculate(delta_time)

        self._end()

    # Internal class. Holds a recorder and the binding of the recorder to an
    # attribute of an object.
    class _RecorderBinding(object):

        def __init__(self, recorder, subjects):
            super(Simulator._RecorderBinding, self).__init__()

            # Save the recorder instance and the subject's attribute.
            self._recorder = recorder
            self._subjects = subjects

        def reset(self):
            # Call the reset handler of the recorder.
            self._recorder.on_simulation_reset(
                [subject.friendly_name for subject in self._subjects])

        def update(self, time, delta_time):
            # This method is called by the simulation after each simulation step
            #   in order to update the recorder.
            for subject in self._subjects:

                value = getattr(subject, self._recorder.attribute_name)

                Simulator._RecorderContext(value, time, delta_time)

                self._recorder.on_observed_value(subject.friendly_name,
                                                 time, value)

    @accepts((1, Recorder), (2, (list, tuple, AbstractSimulationElement, types.NoneType)))
    @returns(Recorder)
    def record(self, recorder, subjects=None):
        """
        record(self, recorder, subjects, conversion=None)

        Adds a recorder to an attribute of an object. The recorder has to
        implement the :class:`Recorder` interface in order to
        receive the data of the given attribute after each simulation step.

        :param recorder: Reference to the recorder object. This object gets
            notified about the changes of the object's attribute at each
            simulation step.
        :type recorder: :class:`Recorder`
        
        :param subjects: The subjects of the recorder, in other words the 
            objects which is attributes has to be recorded. If `None`, the time
            is the only one value recorded, with :func:`on_simulation_step()`.
        :type subjects: list or tuple of :class:`.AbstractSimulationElement`

        :param conversion: Lambda function to convert the actual value taken
            from the attribute before recording. The lambda function gets a
            single parameter ``context`` which is a named tuple with the following
            element:
            
            **value**: The actual value just read from the simulation element's
            attribute.

            **time**: The actual simulation time.

            **delta_time**: The time step just simulated before updating the
            recorder. This can be handy in order to calculate for example the
            power of a delta-energy value.

        :type conversion: lambda context
        
        :returns: Reference to the recoder.
        """
        if subjects is None:
            self._time_recorders.append(recorder)
        else:
            if isinstance(subjects, AbstractSimulationElement):
                subjects = (subjects,)

            self._recorderBindings.append(self._RecorderBinding(recorder, subjects))

            if not recorder in self._recorders:
                self._recorders.append(recorder)
        return recorder class _Unit(object):
    def __init__(self):
    def Quantity(self):
    def Unit(self):
    def __call__(self, *args, **kwargs):
    def __getattr__(self, item):
    def to_si(self, measurement_or_unit):
    def value(self, measurement, unit=None):
    def unit(self, measurement):
    def dimension(self, measurement):
    def convert(self, m, u): class ElectricalPVBus(ElectricalBus):

    def __init__(self, friendly_name, position=Position()):
        """
        __init__(self, friendly_name, position=Position())

        Convenience class to define a PV Bus.

        .. seealso:: :class:`.ElectricalBus`

        :param friendly_name: Friendly name for the element.
            Should be unique within the simulation module.
        :type friendly_name: str

        :param position: Bus geographical position.
            Defaults to Position default value.
        :type position: :class:`.Position`
        """
        super(ElectricalPVBus, self)\
            .__init__(friendly_name, ElectricalBus.Type.PV_BUS, position) class ElectricalPQBus(ElectricalBus):

    def __init__(self, friendly_name, position=Position()):
        """
        __init__(self, friendly_name, position=Position())

        Convenience class to define a PQ Bus.

        .. seealso:: :class:`.ElectricalBus`

        :param friendly_name: Friendly name for the element.
            Should be unique within the simulation module.
        :type friendly_name: str

        :param position: Bus geographical position.
            Defaults to Position default value.
        :type position: :class:`.Position`
        """
        super(ElectricalPQBus, self)\
            .__init__(friendly_name, ElectricalBus.Type.PQ_BUS, position) class ElectricalTransmissionLine(AbstractElectricalTwoPort):

    @units.wraps(None, (None, None, units.metre, units.ohm, units.ohm, units.siemens))
    def __init__(self, friendly_name, length, X, R=0, B=0):
        """
        __init__(self, friendly_name, length, X, R=0, B=0)

        Class for representing a transmission line in an electrical network.
        Its parameters are linked to the
        transmission line PI model represented graphically below::

                                  R              jX
                            +-----------+   +-----------+
            o--->---o-------|           |---|           |-------o--->---o
                    |       +-----------+   +-----------+       |
                    |                                           |
                  -----  jB/2                                 -----  jB/2
                  -----                                       -----
                    |                                           |
                    |                                           |
                   ---                                         ---

        At initialization, in addition to the line ``friendly_name``,
        the line length, the line reactance (``X``), line resistance (``R``) and
        line charging (B) have to be given. ``R`` and ``B`` default to zero.

        :param friendly_name: Friendly name for the line. Should be unique
            within the simulation module, i.e. different for example from the
            friendly name of a bus.
        :type friendly_name: str

        :param length: Line length.
        :type length: length, see :mod:`gridsim.unit`

        :param X: Line reactance.
        :type X: ohm, see :mod:`gridsim.unit`

        :param R: Line resistance.
        :type R: ohm, see :mod:`gridsim.unit`

        :param B: Line charging.
        :type B: siemens, see :mod:`gridsim.unit`
        """

        # HACK: when object is constructed with *args or **kwargs
        if not isinstance(length, (int, float)):
            length = units.value(units.to_si(length))
        if not isinstance(X, (int, float)):
            X = units.value(units.to_si(X))
        if not isinstance(R, (int, float)):
            R = units.value(units.to_si(R))
        if not isinstance(B, (int, float)):
            B = units.value(units.to_si(B))

        # super constructor needs units as it is a "public" function
        super(ElectricalTransmissionLine, self).__init__(friendly_name, X*units.ohm, R*units.ohm)

        if length <= 0:
            raise RuntimeError('Length has to be a strictly positive number')
        if B < 0:
            raise RuntimeError('Line charging B can not be negative number')

        self.length = length
        """
        The transmission line length.
        """
        self.B = B
        """
        The transmission line charging.
        """ class ElectricalSlackBus(ElectricalBus):

    def __init__(self, friendly_name, position=Position()):
        """
        __init__(self, friendly_name, position=Position())

        Convenience class to define a Slack Bus.

        .. seealso:: :class:`.ElectricalBus`

        :param friendly_name: Friendly name for the element.
            Should be unique within the simulation module.
        :type friendly_name: str

        :param position: Bus geographical position.
            Defaults to Position default value.
        :type position: :class:`.Position`
        """
        super(ElectricalSlackBus, self)\
            .__init__(friendly_name, ElectricalBus.Type.SLACK_BUS, position) class ConstantElectricalCPSElement(AbstractElectricalCPSElement):

    @accepts((1, str))
    @units.wraps(None, (None, None, units.watt))
    def __init__(self, friendly_name, power):
        """
        __init__(self, friendly_name, power)

        This class provides the simplest Consuming-Producing-Storing element
        having a constant behavior.

        At initialization, the consumed or produced constant `power` has to
        be provided beside the element `friendly_name`. Power is positive, if
        consumed, negative, if produced. With the 'calculate' method the
        energy consumed or produced during the simulation step is calculated
        from the constant power value.

        :param friendly_name: Friendly name for the element. Should be unique
            within the simulation module.
        :type friendly_name: str

        :param power: The constant consumed (if positive) or produced
            (if negative) power.
        :type power: power, see :mod:`gridsim.unit`

        """
        # HACK: when object is constructed with *args or **kwargs
        if not isinstance(power, (int, float)):
            power = units.value(units.to_si(power))
        super(ConstantElectricalCPSElement, self).__init__(friendly_name)
        self.power = power

    @accepts(((1, 2), (int, float)))
    def calculate(self, time, delta_time):
        """
        calculate(self, time, delta_time)

        Calculates the element's the energy consumed or produced by the element
        during the simulation step.

        :param time: The actual time of the simulator in seconds.
        :type time: time, see :mod:`gridsim.unit`
        :param delta_time: The delta time for which the calculation has to be
            done in seconds.
        :type delta_time: time, , see :mod:`gridsim.unit`
        """
        self._internal_delta_energy = self.power * delta_time class ElectricalNetworkBranch(AbstractElectricalElement):

    @accepts((1, str),
             ((2, 3), ElectricalBus),
             (4, AbstractElectricalTwoPort))
    def __init__(self, friendly_name, from_bus, to_bus, two_port):
        """
        __init__(self, friendly_name, from_bus, to_bus, two_port)

        Class for a branch of an electrical network, i.e. connection between two
        buses (or nodes). It is oriented from one bus to the other. It is
        based on the general :class:`AbstractElectricalElement` class. At
        initialization, in addition to the 'friendly_name', the bus it is
        starting from, and the bus it is going to, have to be given, together
        with the two-port, e.g. transmission line, transformer,... it is made
        of.

        The chosen representation for branch electrical values is: active power
        ``Pi`` and reactive power ``Qi`` flowing into the branch at branch start,
        and active power Po and reactive power ``Qo`` flowing out of the branch at
        branch end. Their default values are None.

        :param friendly_name: Friendly name for the branch.
            Should be unique within the simulation module,
            i.e. different for example from the friendly name of a bus
        :type friendly_name: str

        :param from_bus: Electrical bus from which branch is starting.
        :type from_bus: :class:`.ElectricalBus`

        :param to_bus: Electrical bus to which branch is going.
        :type to_bus: :class:`.ElectricalBus`

        :param two_port: Electrical two-port on the branch,
            e.g. transmission line, transformer, ...
        :type two_port: :class:`.AbstractElectricalTwoPort`
        """
        if from_bus.id is None:
            raise RuntimeError('From_bus bus has not been added to simulator.')
        if to_bus.id is None:
            raise RuntimeError('To_bus bus has not been added to simulator.')

        super(ElectricalNetworkBranch, self).__init__(friendly_name)
        self._from_bus_id = from_bus.id
        self._to_bus_id = to_bus.id
        self._two_port = two_port
        self.Pij = None
        """
        Active power flowing into the branch from the from-bus terminal.
        """
        self.Qij = None
        """
        Reactive power flowing into the branch from the from-bus terminal.
        """
        self.Pji = None
        """
        Active power flowing into the branch from the to-bus terminal.
        """
        self.Qji = None
        """
        Reactive power flowing into the branch from the to-bus terminal.
        """

    @property
    def from_bus_id(self):
        """
        from_bus_id(self)

        Gets the id of the bus the branch is starting from.

        :returns: id of the bus the branch is starting from.
        :rtype: int
        """
        return self._from_bus_id

    @property
    def to_bus_id(self):
        """
        to_bus_id(self)

        Gets the id of the bus the branch is going to.

        :returns: id of the bus the branch is going to.
        :rtype: int
        """
        return self._to_bus_id

    def reset(self):
        """
        reset(self)

        Reset branch electrical values to their default : None
        """
        self.Pij = None
        self.Qij = None
        self.Pji = None
        self.Qji = None class PlotRecorder(Recorder, AttributesGetter):

    @accepts((1, str), ((2, 3), (units.Unit, type, type(None))))
    def __init__(self, attribute_name, x_unit=None, y_unit=None):
        """
        __init__(self, attribute_name, x_unit=None, y_unit=None)

        A :class:`PlotRecorder` can be used to record one or multiple signals
        and plot them at the end of the simulation either to a image file,
        add them to a PDF or show them on a window.

        :param attribute_name: The attribute name to observe. Every
            :class:`.AbstractSimulationElement` with a param of this name can
            be recorded
        :type attribute_name: str
        :param x_unit: The unit of the horizontal axis
        :type x_unit: int, float, type or unit see :mod:`gridsim.unit`
        :param y_unit: The unit of the vertical axis
        :type y_unit: int, float, type or unit see :mod:`gridsim.unit`

        If no information are give to the constructor of a :class:`.Recorder`,
        the data are not preprocessed.

        .. warning:: The preprocess can take a long time if a lot of data are
                     stored in a :class:`.PlotRecorder`.

        *Example:*

        .. literalinclude:: ../../demo/plotrecorder.py
            :linenos:

        * In this example, only new concept are detailed.
        * On line 46, we create a :class:`PlotRecorder` which recorder the
          ``temperature`` attribute in Kelvin.
        * On line 47, we use that recorder in order to record all element in
          the thermal module who have the attribute 'temperature' by using the
          powerful :func:`gridsim.simulation.Simulator.find()` method by asking
          to return a list of all element that have the attribute 'temperature'.
        * On line 51, we create a very similar recorder as before, but
          the attribute will be in celsius and the time in minutes.
        * On line 52, we recorder same objects as above,
        * On line 68, we create a first :class:`.FigureSaver` instance with the
          title 'Temperature (K)' which will use the data returned by the
          :class:`.Recorder` to :func:`save` the data in a figure.
        * On line 69, we also save the data of the other :class:`.Recorder`.
          The figures below show that the result is quite the same except that
          the y-coordinates have the right unit.
        * On line 72, we save the data using another saver: :class:`.CSVSaver`.

        The resulting 3 images of this simulation are show here:

        .. figure:: ../../demo/output/fig1.png
            :align: center

        .. figure:: ../../demo/output/fig2.png
            :align: center

        .. figure:: ../../demo/output/fig3.png
            :align: center

        """
        super(PlotRecorder, self).__init__(attribute_name, x_unit, y_unit)

        self._x = []
        self._y = {}

        # for optimisation when getting results
        self._res_x = []
        self._res_y = {}

    def on_simulation_reset(self, subjects):
        """
        on_simulation_reset(self, subjects)

        :class:`.Recorder` implementation. Prepares the data structures.

        .. seealso:: :func:`gridsim.simulation.Recorder.on_simulation_reset()`
                     for details.
        """
        for subject in subjects:
            self._y[subject] = []

    @accepts((1, (int, float)))
    def on_simulation_step(self, time):
        """
        on_simulation_step(self, time)

        :class:`.Recorder` implementation. Prepares the data structures.

        .. seealso:: :func:`gridsim.simulation.Recorder.on_simulation_step()`
                     for details.
        """
        self._x.append(time)

    @accepts((2, (int, float)), (3, (int, float, units.Quantity)))
    def on_observed_value(self, subject, time, value):
        """
        on_observed_value(self, subject, time, value)

        :class:`.Recorder` implementation. Saves the data in order to
        plot them afterwards.

        .. seealso:: :func:`gridsim.simulation.Recorder.on_observed_value()`
                     for details.
        """
        self._y[subject].append(value)

    def x_values(self):
        """
        x_values(self)

        Retrieves a list of time, on for each :func:`on_observed_value` and
        ordered following the call of :func:`gridsim.simulation.Simulator.step`.

        *Example*::

            sim = Simulator()

            kelvin = PlotRecorder('temperature')
            sim.record(kelvin, sim.thermal.find(has_attribute='temperature'))

            sim.step(17*units.second)
            sim.step(5*units.second)
            sim.step(42*units.second)

        In this example, this function should returns::

            [17, 5, 42]

        :return: time in second
        :rtype: list
        """
        if self._x_unit is None:
            self._res_x = self._x
        elif self._res_x is None or not len(self._res_x) == len(self._x):
            self._res_x = [units.value(units.convert(x, units.to_si(self._x_unit)), self._x_unit)
                           for x in self._x]
        return self._res_x

    @returns(str)
    def x_unit(self):
        """
        x_unit(self)

        Retrieves the unit of the time.

        :return: the time unit, see :mod:`gridsim.unit`
        :rtype: str
        """
        if type(self._x_unit) is units.Quantity:
            return units.unit(self._x_unit)
        else:
            return str(self._x_unit)

    def y_values(self):
        """
        y_values(self)

        Retrieves a map of the recorded data::

            {
                id1: [list of recorded values]
                id2: [list of recorded values]
            }

        with ``idX`` the ``subject`` given in parameter of :func:`on_observed_value`
        and :func:`on_simulation_reset`

        :return: a map associating id to recorded data
        :rtype: dict
        """
        if self._y_unit is None:
            self._res_y = self._y
        elif self._res_y is None or not len(self._res_y) == len(self._y):
                self._res_y = {}
                for key in self._y.keys():
                    if type(self._y_unit) is type:
                        self._res_y[key] = [self._y_unit(y) for y in self._y[key]]
                    else:
                        self._res_y[key] = [units.value(units.convert(y, units.to_si(self._y_unit)), self._y_unit)
                                            for y in self._y[key]]

        return self._res_y

    @returns(str)
    def y_unit(self):
        """
        y_unit(self)

        Retrieves the unit od the recorded data.

        :return: the type or the unit, see :mod:`gridsim.unit`
        :rtype: str
        """
        if type(self._y_unit) is units.Quantity:
            return units.unit(self._y_unit)
        else:
            return str(self._y_unit) class DirectLoadFlowCalculator(AbstractElectricalLoadFlowCalculator):

    @accepts(((1, 2), (int, float, type(None))))
    def __init__(self, s_base=None, v_base=None, is_PV=None, b=None, Yb=None):
        """
        This class implements the direct load flow method to solve the
        power-flow problem.

        .. seealso:: http://en.wikipedia.org/wiki/Power-flow_study#Power-flow_problem_formulation.

        The method is described for example in Hossein Seifi, Mohammad Sadegh
        Sepasian, Electric Power System Planning: Issues, Algorithms and
        Solutions, pp. 246-248

        At initialization the user has to give the reference power value
        `s_base` (all power values are then given relative to this reference
        value), the reference voltage value `v_base` (all voltage values are
        then given relative to this reference value), a boolean array `is_PV`
        specifying which one among the buses is a :class:`.ElectricalPVBus`
        (the bus with 1st position is slack, the others non
        :class:`.ElectricalPVBus` buses are :class:`.ElectricalPQBus` buses), an
        integer array `b` specifying for each branch the bus id it is starting
        from and the bus id it is going to, a complex array `Yb` specifying
        admittances of each network branch.
        """
        super(DirectLoadFlowCalculator, self).__init__()

        self._invBvq = None
        self._bA = None

        if s_base is not None and v_base is not None and is_PV is not None \
                and b is not None and Yb is not None:
            self.update(s_base, v_base, is_PV, b, Yb)

    @accepts(((1, 2), (int, float)))
    def update(self, s_base, v_base, is_PV, b, Yb):
        """
        update(self, s_base, v_base, is_PV, b, Yb)

        Updates values of the calculator.

        :param s_base: reference power value
        :type s_base: float
        :param v_base: reference voltage value
        :type v_base: float
        :param is_PV: N-long vector specifying which bus is of type
            :class:`.ElectricalPVBus`, where N is the number of buses including
            slack.
        :type is_PV: 1-dimensional numpy array of boolean
        :param b: Mx2 table containing for each branch the ids of start and end
            buses.
        :type b: 2-dimensional numpy array of int
        :param Yb: Mx4 table containing the admittances Yii, Yij, Yjj, and Yji
            of each branch.
        :type Yb: 2-dimensional numpy array of complex

        .. note:: variable names have been chosen to match the notation used in
                  Hossein Seifi, Mohammad Sadegh Sepasian, Electric Power System
                  Planning: Issues, Algorithms and Solutions, pp. 246-248
        """
        super(DirectLoadFlowCalculator, self).update(s_base, v_base,
                                                     is_PV, b, Yb)



        # compute matrix B from admittance matrix Y
        # off-diagonal element are equal to minus imaginary part of Y element
        B = -np.imag(self._Y)
        # diagonal element are equal to minus sum of off-diagonal element
        np.fill_diagonal(B, np.zeros(B.shape[0]))
        np.fill_diagonal(B, -B.sum(1))
        # compute inverse of B after removing first row and first column
        # (corresponding to slack bus)
        self._invBvq = np.linalg.inv(B[1:, 1:])

        # build bA matrix from branch susceptances as sparse matrix
        self._bA = lil_matrix((self._nBr, self._nBu))

        for i_branch in range(0, self._nBr):
            # this is minus the susceptance value
            mb = B[self._b[i_branch, 0], self._b[i_branch, 1]]
            self._bA[i_branch, self._b[i_branch, 0]] = -mb
            self._bA[i_branch, self._b[i_branch, 1]] = mb

        # change sparse matrix representation
        self._bA = self._bA.tocsr()

    @accepts((5, bool))
    def calculate(self, P, Q, V, Th, scaled):
        """
        calculate(self, P, Q, V, Th, scaled)

        Computes all bus electrical values determined by the network, i.e

        - takes as input active powers P of all non-slack buses.
        - outputs voltage angles of all buses.

        Reactive powers and voltage amplitudes are not used. However, to be
        compatible with the N-long vectors with the method of the base class,
        all N-long vectors P, Q, V, and Th have to be passed to this method.

        :param P: N-long vector of bus active powers, where N is the number of
            buses including slack.
        :type P: 1-dimensional numpy array of float
        :param Q: N-long vector of bus reactive powers, where N is the number of
            buses including slack.
        :type Q: 1-dimensional numpy array of float
        :param V: N-long vector of bus voltage amplitudes, where N is the number
            of buses including slack.
        :type V: 1-dimensional numpy array of float
        :param Th: N-long vector of bus voltage angles, where N is the number of
            buses including slack.
        :type Th: 1-dimensional numpy array of float
        :param scaled: specifies whether electrical input values are scaled or
            not
        :type scaled: boolean

        :return: modified [P, Q, V, Th]
        :rtype: a list of 4 element
        """
        # check input arguments and save them to internal
        # variables _P, _Q, and _V
        # variables _P, _Q, and _V
        self._read_calculate_args(P, Q, V, Th, scaled)
        # initialize voltage angles internal variable _Th to 0.0
        self._Th = np.zeros([self._nBu, 1])

        # compute slack active power (based on the assumption that there are no
        # branch losses)

        # update intern variable
        self._P[0] = -sum(self._P[1:])
        # update external variable
        if scaled:
            P[0] = self._P[0]
        else:
            P[0] = self.s_base * self._P[0]

        # TODO: decide if Q and V should be updated to be consistent
        # TODO (V should be set to ones)

        # vector of voltage angles
        # update intern variable
        self._Th = np.concatenate(([0.0], np.dot(self._invBvq, self._P[1:])))


        # return external variable
        self._calculate_done = True

        return [self._P, self._Q, self._V, self._Th]

    @accepts((1, bool))
    def get_branch_power_flows(self, scaled):
        """
        get_branch_power_flows(self, scaled)

        Computes all branch power flows. Cannot be called before method
        :func:`AbstractElectricalLoadFlowCalculator.calculate`.
        To be compatible with the generic method, returns a
        tuple made of 4 1-dimensional numpy arrays:

        - `Pij` containing branch input active powers,
        - `Qij` containing branch input reactive powers,
        - `Pji` containing branch output active powers,
        - `Qji` containing branch output reactive powers.

        Here, `Qij` and `Qji` are `None`. `Pji` is equal to `-Pij`.

        :param scaled: specifies whether output power flows have to be scaled
            or not
        :type scaled: boolean
        :returns: 4 M-long vector, where M is the number of branches

            :Pij: vector of active powers entering the branch from the from-bus
                  termination
            :Qij: None
            :Pji: vector of active powers entering the branch from the to-bus
                  termination
            :Qji: None

        :rtype: mixed tuple of None and one-dimensional numpy arrays of float
        """
        if self._P is None:
            # calculate has not been called
            raise RuntimeError('The calculate method has to be called first!')

        # branch active powers
        Pbr = self._bA.dot(self._Th)

        if not scaled:
            Pbr *= self.s_base

        return Pbr, None, -Pbr, None

    @accepts((1, bool))
    def get_branch_max_currents(self, scaled):
        """
        get_branch_max_currents(self, scaled)

        Returns active powers entering the branch from the from-bus termination.

        :param scaled: specifies whether output power flows have to be scaled
            or not
        :type scaled: boolean
        :returns: a M-long vector of active powers entering the branch from the
                  from-bus termination, where M is the number of branches

        .. warning:: This function is currently untested.
        """

        if self._P is None:
            # calculate has not been called
            raise RuntimeError('The calculate method has to be called first!')

        # it can be shown that in per unit the branch current magnitude is equal
        # to the branch active power (Vk~=1.0)

        # branch active powers
        [Pbr, _, _, _] = self.get_branch_power_flows(True)
        Ibr = Pbr  # FIXME Pij or Pji ???

        if not scaled:
            Ibr *= (self.s_base * self.v_sc) #

        return Ibr
import unittest
import numpy as np
from gridsim.simulation import Simulator
from gridsim.unit import units
from gridsim.electrical.network import ElectricalPVBus, ElectricalPQBus, \
    ElectricalTransmissionLine, ElectricalSlackBus
from gridsim.electrical.element import ConstantElectricalCPSElement
from gridsim.electrical.core import ElectricalNetworkBranch
from gridsim.recorder import PlotRecorder
from gridsim.electrical.loadflow import DirectLoadFlowCalculator
# This program checks whether the default DirectLoadFlowCalculator present in
# the gridsim.electrical module performs correctly. It computes the example
# given in Hossein Seifi, Mohammad Sadegh Sepasian, Electric Power System
# Planning: Issues, Algorithms and Solutions, pp. 247-248 and compare results
# with those of the reference.
#
# Example input:
#     _                                        _
#   /   \                                    /   \
#  |    |                                   |    |
#  \   /                                    \   /
#    +    -------------------------------     +
#    |    |                             |     |
#  --+----+--                         --+-----+--
#         |                             |
#         |                             |
#         |                             |
#         -------------     -------------
#                     |     |
#                   --+--+--+--
#                        |
#                        V
# Loads and generations
#-------------------------------------------------------------------
# Bus number    Bus type       PD (MW)      QD (MVAr)   PG (MW)
#-------------------------------------------------------------------
# 1             Slack           0           0           Unknown
# 2             PV              10          5           63
# 3             PQ              90          30          0
#-------------------------------------------------------------------
#
# Branches
#-------------------------------------------------------------------
# Line number   From bus        To bus      X (p.u.)    Rating (MVA)
#-------------------------------------------------------------------
# 1             1               2           0.0576      250
# 2             2               3           0.092       250
# 3             1               3           0.17        150
#-------------------------------------------------------------------






class TestEsimDLF(unittest.TestCase):

    def test_reference(self):

        # Initialize the Gridsim simulator and get a reference to
        # its electrical part
        sim = Simulator()
        esim = sim.electrical
        esim.load_flow_calculator = DirectLoadFlowCalculator()

        # network initialization
        #----------------------

        # add buses to simulator
        # slack bus has been automatically added
        bus1 = esim.bus('Slack Bus')
        bus2 = esim.add(ElectricalPVBus('Bus 2'))
        bus3 = esim.add(ElectricalPQBus('Bus 3'))

        # add branches to simulator
        # line length is arbitrarily set to 1.0
        bra1 = esim.connect('Branch 1-2', esim.bus('Slack Bus'),
                            esim.bus('Bus 2'),
                            ElectricalTransmissionLine('Line 1',
                            1.0*units.metre, 0.0576*units.ohm))
        # line length is arbitrarily set to 1.0
        bra2 = esim.connect('Branch 2-3', esim.bus('Bus 2'),
                            esim.bus('Bus 3'),
                            ElectricalTransmissionLine('Line 2',
                            1.0*units.metre, 0.092*units.ohm))
        # line length is arbitrarily set to 1.0
        bra3 = esim.connect('Branch 1-3', esim.bus('Slack Bus'),
                            esim.bus('Bus 3'),
                            ElectricalTransmissionLine('Line 3',
                            1.0*units.metre, 0.17*units.ohm))

        # input buses electrical values
        #------------------------------
        esim.attach('Bus 2', ConstantElectricalCPSElement('GD2', -.53*units.watt))

        esim.attach('Bus 3', ConstantElectricalCPSElement('GD3', .9*units.watt))