def problem_has_exact_solution(problem, checking_obj=None):
    """Convenience accessor for exact solution.

    Parameters
    ----------
    problem : :py:class:`.IProblem`
        The problem to check for an exact solution function.
    checking_obj : object
        *(optional)*
        The object calling this function for a meaningful error message.
        For debugging purposes only.

    Returns
    -------
    has_exact_solution : :py:class:`bool`
        :py:class:`True` if exact solution was given, :py:class:`False` otherwise

    Raises
    ------
    ValueError :
        If the given problem is not an instance of :py:class:`.IProblem`.
    """
    assert_is_instance(
        problem, IProblem, message="It needs to be a problem to have an exact solution.", checking_obj=checking_obj
    )
    return isinstance(problem, HasExactSolutionMixin)

 def time_interval(self, value):
     assert_is_instance(value, np.ndarray, descriptor="Time Interval", checking_obj=self)
     assert_condition(value.size == 2, ValueError,
                      message="Time Interval must have two values: NOT %d" % value.size)
     self.validate_time_interval(start=value[0], end=value[1])
     self.time_start = value[0]
     self.time_end = value[1]

    def __init__(self, restriction_stencil, decrease_in_points=None, center=None):
        assert_is_instance(restriction_stencil, np.ndarray, "Not an ndarray")
        self.rst_stencil = Stencil(restriction_stencil, center)


        if decrease_in_points is None:
            self.dip = np.asarray(self.rst_stencil.shape) - 1
        elif isinstance(decrease_in_points, np.ndarray) and \
                        decrease_in_points.ndim == restriction_stencil.ndim:
            self.dip = decrease_in_points
        elif isinstance(decrease_in_points, int):
            self.dip = np.asarray([decrease_in_points]*restriction_stencil.ndim)
        else:
            raise ValueError("Wrong decrease in points")
        self.dim = restriction_stencil.ndim

        if self.dim == 1:
            self.eval = self.evalA_1D
        elif self.dim == 2:
            self.eval = self.evalA_2D
        elif self.dim == 3:
            self.eval = self.evalA_3D
        else:
            raise NotImplementedError("More than 3 dimensions " +
                                      "are not implemented")

    def __init__(self, *args, **kwargs):
        """
        Parameters
        ----------
        exact_function : :py:class:`callable`
            *(optional)*
            If given initializes the problem with the exact solution function.
        """
        assert_is_instance(
            self,
            IProblem,
            descriptor="Problem needs to be a IProblem first: NOT %s" % class_name(self),
            checking_obj=self,
        )
        self._exact_function = None
        if "exact_function" in kwargs:
            self.exact_function = kwargs["exact_function"]

        self._strings["exact"] = None
        if "strings" in kwargs:
            if "exact" in kwargs["strings"]:
                assert_is_instance(
                    kwargs["strings"]["exact"],
                    str,
                    descriptor="String representation of Exact Function",
                    checking_obj=self,
                )
                self._strings["exact"] = kwargs["strings"]["exact"]

    def add_level_transition(self, transitioner, coarse_level, fine_level):
        """Adds specialized level transitioner for specified levels.

        Parameters
        ----------
        transitioner : :py:class:`.ILevelTransitionProvider`
            Special level transitioner for specified prolongation and restringation between given coarse and fine level.

        coarse_level : :py:class:`int`
            Coarse level of the transitioner.

        fine_level : :py:class:`int`
            Fine level of the transitioner.

        Raises
        ------
        ValueError
            if ``transitioner`` is not an :py:class:`.ILevelTransitionProvider`
        """
        assert_is_instance(transitioner, ILevelTransitionProvider, descriptor="Level Transitioner", checking_obj=self)

        # extend/initialize level_transition_provider map if necessary
        if coarse_level not in self._level_transitioners:
            self._level_transitioners[coarse_level] = {}

        self._level_transitioners[coarse_level][fine_level] = transitioner

    def __init__(self, *args, **kwargs):
        super(HeatEquation, self).__init__(*args, **kwargs)

        # HasExactSolutionMixin.__init__(self, *args, **kwargs)

        self._thermal_diffusivity = kwargs.get('thermal_diffusivity', 1.0)

        if self.time_start is None:
            self.time_start = 0.0
        if self.time_end is None:
            self.time_end = 1.0
        if self.initial_value is None:
            self.initial_value = np.zeros(self.dim_for_time_solver)

        if isinstance(self.thermal_diffusivity, complex):
            self.numeric_type = np.complex

        self._mg_stencil = kwargs.get('mg_stencil')
        self._mg_level = kwargs.get('mg_level')
        self._direct_solvers = {}

        if kwargs.get('delta_times_for_time_levels') is not None and self._mg_level is not None:
            assert_is_instance(kwargs['delta_times_for_time_levels'], (list, np.ndarray),
                               descriptor="Delta Times for Time Levels", checking_obj=self)
            for time_level in kwargs['delta_times_for_time_levels']:
                assert_is_instance(kwargs['delta_times_for_time_levels'][time_level], (list, np.ndarray),
                                   descriptor="Delta Times for Time Level %d" % time_level, checking_obj=self)
                for delta_time in kwargs['delta_times_for_time_levels'][time_level]:
                    self.initialize_direct_space_solver(time_level, delta_time, kwargs['mg_level'])

    def write_buffer(self, tag=None, **kwargs):
        """Writes data into this communicator's buffer

        Parameters
        ----------
        value :
            data values to be send to the next solver
        time_point : :py:class:`float`
            time point of the data values
        flag : :py:class:`.Message.SolverFlag`
            message flag

        Raises
        ------
        ValueError

            * if no arguments are given
            * if ``time_point`` is not a :py:class:`float`
        """
        assert_condition(len(kwargs) > 0, ValueError, "At least one argument must be given.", self)

        if tag not in self._buffer:
            self._buffer[tag] = Message()
            self.write_buffer(tag=tag, **kwargs)

        if "value" in kwargs:
            self._buffer[tag].value = deepcopy(kwargs["value"])

        if "time_point" in kwargs:
            assert_is_instance(kwargs["time_point"], float, descriptor="Time Point", checking_obj=self)
            self._buffer[tag].time_point = deepcopy(kwargs["time_point"])

        if "flag" in kwargs:
            self._buffer[tag].flag = deepcopy(kwargs["flag"])

    def __init__(self, stencil_list, center=None, *args, **kwargs):
        """init

        """
        # in the case of just one stencil
        if isinstance(stencil_list, np.ndarray) and stencil_list.ndim == 1:
            self.mode = "own"
            self.increase_of_points = 2
            self.stencil_list = []
            self.stencil_list.append(stencil_list)
            self.center = []
            if center is None:
                self.center.append(np.floor(np.asarray(self.stencil[0].shape)*0.5))
            elif isinstance(center, int) and center < self.stencil[0].size:
                self.center.append(center)
            else:
                raise ValueError("Wrong argument")
        # now the case of a stencil set
        elif isinstance(stencil_list, list):
            # check if each one is a nd.array
            for stencil in stencil_list:
                assert_is_instance(stencil, Stencil,
                                   "not real stencil", self)
            # each stencil
                self.mode = "list"
                self.stencil_list = stencil_list
                self.increase_of_points = len(stencil_list)
        if self.mode == "own":
            self.eval = self.eval_own
        elif self.mode == "list":
            self.eval = self.eval_list
        else:
            raise RuntimeError("Something went terribly wrong")
        super().__init__(*args, **kwargs)

def problem_has_direct_implicit(problem, checking_obj=None):
    """Convenience checker for existence of a direct implicit formulation of a problem.

    Parameters
    ----------
    problem : :py:class:`.IProblem`
        The problem to check for a direct implicit formulation.
    checking_obj : :py:class:`object`
        *(optional)*
        The object calling this function for a meaningful error message.
        For debugging purposes only.

    Returns
    -------
    has_direct_impl : :py:class:`bool`
        :py:class:`True` if exact solution was given, :py:class:`False` otherwise

    Raises
    ------
    ValueError :
        If the given problem is not an instance of :py:class:`.IProblem`.
    """
    assert_is_instance(problem, IProblem,
                       message="It needs to be a problem to have a direct implicit formula.", checking_obj=checking_obj)
    return isinstance(problem, HasDirectImplicitMixin)

    def evaluate(self, data, **kwargs):
        """Applies this integrator to given data in specified time interval.

        Parameters
        ----------
        data : :py:class:`numpy.ndarray`
            Data vector of the values at given time points.
            Its length must equal the number of integration nodes.
        time_start : :py:class:`float`
            *(optional)*
            Begining of the time interval to integrate over.
        time_end : :py:class:`float`
            *(optional)*
            End of the time interval to integrate over.

        Raises
        ------
        ValueError :

            * if ``data`` is not a :py:class:`numpy.ndarray`
            * if either ``time_start`` or ``time_end`` are not given
            * if ``time_start`` is larger or equals ``time_end``
        """
        assert_is_instance(data, np.ndarray, descriptor="Data to integrate", checking_obj=self)
        assert_condition("time_start" in kwargs or "time_end" in kwargs,
                         ValueError, message="Either start or end of time interval need to be given.",
                         checking_obj=self)
        assert_condition(kwargs["time_start"] < kwargs["time_end"],
                         ValueError,
                         message="Time interval need to be non-zero positive: [{:f}, {:f}]"
                                 .format(kwargs["time_start"], kwargs["time_end"]),
                         checking_obj=self)

    def add_coefficient(self, coefficient, power):
        """Adds or sets the coefficient :math:`c` of :math:`cx^p` for a specific :math:`p`.

        The polynomial gets automatically extended to hold the new coefficient in case it didn't included the specified
        power previously.
        Unset, but skipped powers have a coefficient of zero by default.

        Parameters
        ----------
        coefficient : :py:class:`float`
            Coefficient :math:`c` of :math:`cx^p`.
        power : :py:class:`int`
             Power :math:`p` of :math:`cx^p`.

        Examples
        --------
        >>> polyWeights = PolynomialWeightFunction()
        >>> # To set the coefficient of x^3 to 3.14 use:
        >>> polyWeights.add_coefficient(3.14, 3)
        >>> # Similar, to set the constant coefficient 42, e.i. 42*x^0, use:
        >>> polyWeights.add_coefficient(42, 0)
        """
        assert_is_instance(power, int, descriptor="Power", checking_obj=self)
        assert_condition(power >= 0, ValueError,
                         message="Power must be zero or positive: {:d}".format(power), checking_obj=self)

        if self._coefficients.size <= power + 1:
            self._coefficients = np.resize(self._coefficients, (power + 1))

        self._coefficients[power] = coefficient

 def value(self, value):
     assert_condition(not self.finalized, AttributeError,
                      message="Cannot change this solution data storage any more.", checking_obj=self)
     assert_is_instance(value, np.ndarray, descriptor="Values", checking_obj=self)
     self._dim = value.shape
     self._numeric_type = value.dtype
     self._data = value

    def construct_space_tensor(self, number_of_points_list, stencil=None):
        """Constructs the Spacetensor which is important for the evaluation in the case of Dirichlet boundary conditions

        Parameters
        ----------
        number_of_points_list : :py:class:`int` or :py:class:`numpy.ndarray`
            Number of points which will be distributed equiv-spaced on the grid
        """
        if isinstance(number_of_points_list, (int, float, complex)):
            assert_is_instance(stencil, Stencil, descriptor="Stencil", checking_obj=self)
            npoints = int(number_of_points_list)
            # LOG.debug("Your number %s was modified to %s" % (number_of_points_list, npoints))
            assert_condition(npoints > max(stencil.arr.shape), ValueError,
                             message="Not enough points for the stencil", checking_obj=self)
            npoints = np.asarray([npoints] * len(self.spacial_dim))
        elif isinstance(number_of_points_list, np.ndarray):
            assert_condition(len(number_of_points_list.shape) == 1 and number_of_points_list.size == len(self.spacial_dim),
                             ValueError, message="The number_of_points list is wrong", checking_obj=self)
            npoints = np.floor(number_of_points_list)
        else:
            raise ValueError("Wrong number of points list")

        # first we assign the memory using numpy
        # spt(npoints,dim)
        self._act_npoints = npoints
        lspc = []
        for i in range(len(self.spacial_dim)):
            lspc.append(np.linspace(self._geometry[i, 0], self._geometry[i, 1], npoints[i]))
        if len(self.spacial_dim) > 1:
            space_tensor = np.asarray(np.meshgrid(*lspc))
        else:
            space_tensor = np.linspace(self._geometry[0, 0], self._geometry[0, 1], npoints)

        return space_tensor

def find_root(fun, x0, method="hybr"):
    """Wrapper around SciPy's generic root finding algorithm to support complex numbers.

    SciPy's generic root finding algorithm (``scipy.optimize.root``) is not able to deal with functions returning
    and/or accepting arrays with complex numbers.

    This wrapped call will first convert all arrays of complex numbers into arrays of floats while splitting each
    complex number up into two floats.

    Parameters
    ----------
    fun : :py:class:`callable`
        Complex function to find the root of

    x0 : :py:class:`numpy.ndarray`
        Initial guess.

    method : :py:class:`str`
        Root finding method to be used.
        See ``scipy.optimize.root`` for details.

    Returns
    -------
    Same solution object as ``scipy.optimize.root`` but with ``x`` being converted back including complex numbers.

    Examples
    --------
        from pypint.plugins.implicit_solvers.find_root import find_root
        import numpy
        fun = lambda x: (-1.0 + 1.0j) * x
        sol = find_root(fun, numpy.array([0.0]))
    """
    assert_is_instance(x0, np.ndarray, descriptor="Initial Guess")
    assert_is_callable(fun, descriptor="Function to find root of")
    assert_is_instance(method, str, descriptor="Root finding method")

    _value_map = {}
    _transformed_size = 0
    _transform_necessary = False
    for i in range(0, x0.size):
        if isinstance(x0[i], complex):
            _value_map[i] = [_transformed_size, _transformed_size + 1]
            _transformed_size += 2
            _transform_necessary = True
        else:
            _value_map[i] = [_transformed_size]
            _transformed_size += 1

    if _transform_necessary:
        _wrapped_func = \
            lambda x_next: _transform_to_real(fun(_transform_to_complex(x_next, _value_map)),
                                              _value_map, _transformed_size)
        sol = root(fun=_wrapped_func, x0=_transform_to_real(x0, _value_map, _transformed_size), method=method)
    else:
        sol = root(fun=fun, x0=x0, method=method)

    if sol.success and _transform_necessary:
        sol.x = _transform_to_complex(sol.x, _value_map)
    return sol

 def implicit_solve(self, next_x, func, method="unused", **kwargs):
     """A solver for the implicit equations.
     """
     assert_is_instance(next_x, np.ndarray, descriptor="Initial Guess", checking_obj=self)
     assert_is_callable(func, descriptor="Function of RHS for Implicit Solver", checking_obj=self)
     sol = scop.newton_krylov(func, next_x.reshape(-1))
     assert_is_instance(sol, np.ndarray, descriptor="Solution", checking_obj=self)
     return sol.reshape(self.dim_for_time_solver)

    def _init_new_interval(self, start):
        """Initializes a new work interval

        Parameters
        ----------
        start : :py:class:`float`
            start point of new interval

        Returns
        -------
        has_work : :py:class:`bool`
            :py:class:`True` if new interval have been initialized;
            :py:class:`False` if no new interval have been initialized (i.e. new interval end would exceed end of time
            given by problem)
        """
        assert_is_instance(start, float, descriptor="Time Point", checking_obj=self)

        if start + self._dt > self.problem.time_end:
            return False

        if self.state and start == self.state.initial.time_point:
            return False

        self._init_new_state()

        # set width of current interval
        self.state.delta_interval = self._dt

        # compute time step and node distances
        self._deltas["t"] = self.state.delta_interval / self.num_time_steps  # width of a single time step (equidistant)

        # start time points of time steps
        self.__time_points["steps"] = np.linspace(start, start + self._dt, self.num_time_steps + 1)

        # initialize and transform integrator for time step width
        self._integrator.init(
            self.__nodes_type,
            self.__num_nodes,
            self.__weights_type,
            interval=np.array([self.__time_points["steps"][0], self.__time_points["steps"][1]], dtype=np.float),
        )

        self.__time_points["nodes"] = np.zeros((self.num_time_steps, self.num_nodes), dtype=np.float)
        _deltas_n = np.zeros(self.num_time_steps * (self.num_nodes - 1) + 1)

        # copy the node provider so we do not alter the integrator's one
        _nodes = deepcopy(self._integrator.nodes_type)
        for _t in range(0, self.num_time_steps):
            # transform Nodes (copy) onto new time step for retrieving actual integration nodes
            _nodes.interval = np.array([self.__time_points["steps"][_t], self.__time_points["steps"][_t + 1]])
            self.__time_points["nodes"][_t] = _nodes.nodes.copy()
            for _n in range(0, self.num_nodes - 1):
                _i = _t * (self.num_nodes - 1) + _n
                _deltas_n[_i + 1] = _nodes.nodes[_n + 1] - _nodes.nodes[_n]
        self._deltas["n"] = _deltas_n[1:].copy()

        return True

    def run(self, state, **kwargs):
        """Apply the solver core to the current state

        Parameters
        ----------
        state : :py:class:`.ISolverState`
            Current state of the solver.
        """
        assert_is_instance(state, ISolverState, descriptor="Solver's State", checking_obj=self)

    def compute_error(self, state, **kwargs):
        """Computes the error of the current state

        Parameters
        ----------
        state : :py:class:`.ISolverState`
            Current state of the solver.
        """
        assert_is_instance(state, ISolverState, descriptor="Solver's State", checking_obj=self)

 def __lt__(self, other):
     assert_is_instance(other, StepSolutionData,
                        message="Can not compare StepSolutionData with {}".format(class_name(other)),
                        checking_obj=self)
     return (
         self.dim == other.dim
         and self.numeric_type == other.numeric_type
         and self.time_point < other.time_point
     )

 def initial_value(self, initial_value):
     assert_is_instance(initial_value, np.ndarray, descriptor="Initial Value", checking_obj=self)
     assert_condition(
         initial_value.shape == self.dim_for_time_solver,
         ValueError,
         message="Initial Values shape must match problem DOFs: %s != %s"
         % (initial_value.shape, self.dim_for_time_solver),
         checking_obj=self,
     )
     self._initial_value = initial_value