"""
Simulation infrastructure with helpers and data structures for preprocessing of the given equations
and functions for postprocessing of simulation data.
"""
import warnings
from abc import ABCMeta, abstractmethod
from collections import OrderedDict
from copy import copy
from itertools import chain
import numpy as np
from scipy.integrate import ode
from scipy.interpolate import interp1d
from scipy.linalg import block_diag
from .core import (Domain, Parameters, Function, integrate_function,
calculate_scalar_product_matrix,
dot_product_l2, sanitize_input,
StackedBase, get_weight_transformation,
get_transformation_info,
EvalData, project_on_bases)
from .placeholder import (Scalars, TestFunction, Input, FieldVariable,
EquationTerm, get_common_target, get_common_form)
from .registry import get_base, register_base
__all__ = ["SimulationInput", "SimulationInputSum", "WeakFormulation",
"parse_weak_formulation",
"create_state_space", "StateSpace", "simulate_state_space",
"simulate_system", "simulate_systems",
"get_sim_result", "evaluate_approximation",
"parse_weak_formulations",
"get_sim_results", "set_dominant_labels", "CanonicalEquation",
"CanonicalForm"]
class EmptyInput(SimulationInput):
def __init__(self, dim):
SimulationInput.__init__(self)
self.dim = dim
def _calc_output(self, **kwargs):
return dict(output=np.zeros((len(np.atleast_1d(kwargs['time'])), self.dim)))
[docs]class StateSpace(object):
r"""
Wrapper class that represents the state space form of a dynamic system where
.. math::
\boldsymbol{\dot{x}}(t) &= \sum\limits_{k=0}^{L}\boldsymbol{A}_{k} \boldsymbol{x}^{p_k}(t)
+ \sum\limits_{j=0}^{V} \sum\limits_{k=0}^{L}\boldsymbol{B}_{j, k} \frac{\mathrm{d}^j u^{p_k}}{\mathrm{d}t^j}(t) \\
\boldsymbol{y}(t) &= \boldsymbol{C}\boldsymbol{x}(t) + \boldsymbol{D}u(t)
which has been approximated by projection on a base given by weight_label.
Args:
a_matrices (dict): State transition matrices :math:`\boldsymbol{A}_{p_k}`
for the corresponding powers of :math:`\boldsymbol{x}`
b_matrices (dict): Cascaded dictionary for the input matrices :math:`\boldsymbol{B}_{j, k}` in the sequence:
temporal derivative order, exponent .
input_handle: function handle, returning the system input :math:`u(t)`
c_matrix: :math:`\boldsymbol{C}`
d_matrix: :math:`\boldsymbol{D}`
"""
def __init__(self, a_matrices, b_matrices, base_lbl=None,
input_handle=None, c_matrix=None, d_matrix=None):
self.C = c_matrix
self.D = d_matrix
self.base_lbl = base_lbl
# mandatory
if isinstance(a_matrices, np.ndarray):
self.A = {1: a_matrices}
else:
self.A = a_matrices
if 0 not in self.A:
# this is the constant term (power 0) aka the f-vector
self.A[0] = np.zeros((self.A[1].shape[0],))
# optional
if isinstance(b_matrices, np.ndarray):
# fake import order and power for backward compatibility
self.B = {0: {1: b_matrices}}
else:
self.B = b_matrices
# TODO calculate available order
available_power = 1
if self.B is None:
self.B = {0: {available_power: np.zeros((self.A[available_power].shape[0], available_power))}}
if self.C is None:
self.C = np.zeros((available_power, self.A[available_power].shape[1]))
if self.D is None:
self.D = np.zeros((self.C.shape[0], np.atleast_2d(self.B[0][available_power]).T.shape[1]))
self.input = input_handle
if self.input is None:
self.input = EmptyInput(self.B[0][available_power].shape[1])
if not callable(self.input):
raise TypeError("input must be callable!")
# TODO export cython code?
[docs] def rhs(self, _t, _q):
r"""
Callback for the integration of the dynamic system, described by this object.
Args:
_t (float): timestamp
_q (array): weight vector
Returns:
(array): :math:`\boldsymbol{\dot{x}}(t)`
"""
q_t = self.A[0]
for p, a_mat in self.A.items():
[docs] q_t = q_t + a_mat @ np.power(_q, p)
# TODO make compliant with definition of temporal derived input
u = self.input(time=_t, weights=_q, weight_lbl=self.base_lbl)
for o, p_mats in self.B.items():
for p, b_mat in p_mats.items():
# q_t = q_t + (b_mat @ np.power(u, p)).flatten()
q_t = q_t + b_mat @ np.power(u, p)
return q_t
def simulate_system(weak_form, initial_states,
temporal_domain, spatial_domain,
derivative_orders=(0, 0), settings=None):
r"""
Convenience wrapper for :py:func:`.simulate_systems`.
Args:
weak_form (:py:class:`.WeakFormulation`): Weak formulation of the system
to simulate.
initial_states (numpy.ndarray): Array of core.Functions for
:math:`x(t=0, z), \dot{x}(t=0, z), \dotsc, x^{(n)}(t=0, z)`.
temporal_domain (:py:class:`.Domain`): Domain object holding information
for time evaluation.
spatial_domain (:py:class:`.Domain`): Domain object holding information
for spatial evaluation.
derivative_orders (tuple): tuples of derivative orders (time, spat) that
shall be evaluated additionally as values
settings: Integrator settings, see :py:func:`.simulate_state_space`.
"""
ics = sanitize_input(initial_states, Function)
initial_states = {weak_form.name: ics}
spatial_domains = {weak_form.name: spatial_domain}
derivative_orders = {weak_form.name: derivative_orders}
res = simulate_systems([weak_form], initial_states, temporal_domain, spatial_domains, derivative_orders, settings)
return res
[docs]def simulate_systems(weak_forms, initial_states, temporal_domain, spatial_domains, derivative_orders=None,
settings=None):
"""
Convenience wrapper that encapsulates the whole simulation process.
Args:
weak_forms ((list of) :py:class:`.WeakFormulation`): (list of) Weak
formulation(s) of the system(s) to simulate.
initial_states (dict, numpy.ndarray): Array of core.Functions for
:math:`x(t=0, z), \dot{x}(t=0, z), \dotsc, x^{(n)}(t=0, z)`.
temporal_domain (:py:class:`.Domain`): Domain object holding
information for time evaluation.
spatial_domains (dict): Dict with :py:class:`.Domain` objects holding
information for spatial evaluation.
derivative_orders (dict): Dict, containing tuples of derivative orders
(time, spat) that shall be evaluated additionally as values
settings: Integrator settings, see :py:func:`.simulate_state_space`.
Note:
The *name* attributes of the given weak forms must be unique!
Return:
list: List of :py:class:`.EvalData` objects, holding the results for the
FieldVariable and demanded derivatives.
"""
if derivative_orders is None:
derivative_orders = dict([(lbl, (0,0))for lbl in spatial_domains])
weak_forms = sanitize_input(weak_forms, WeakFormulation)
print("simulate systems: {}".format([f.name for f in weak_forms]))
print(">>> parse weak formulations")
canonical_equations = parse_weak_formulations(weak_forms)
print(">>> create state space system")
state_space_form = create_state_space(canonical_equations)
print(">>> derive initial conditions")
q0 = project_on_bases(initial_states, canonical_equations)
print(">>> perform time step integration")
sim_domain, q = simulate_state_space(state_space_form, q0, temporal_domain, settings=settings)
print(">>> perform postprocessing")
results = get_sim_results(sim_domain, spatial_domains, q, state_space_form, derivative_orders=derivative_orders)
print(">>> finished simulation")
return results
[docs]def get_sim_result(weight_lbl, q, temp_domain, spat_domain, temp_order, spat_order, name=""):
"""
Create handles and evaluate at given points.
Args:
weight_lbl (str): Label of Basis for reconstruction.
temp_order: Order or temporal derivatives to evaluate additionally.
spat_order: Order or spatial derivatives to evaluate additionally.
q: weights
spat_domain (:py:class:`.Domain`): Domain object providing values for
spatial evaluation.
temp_domain (:py:class:`.Domain`): Time steps on which rows of q are
given.
name (str): Name of the WeakForm, used to generate the data set.
"""
data = []
# temporal
ini_funcs = get_base(weight_lbl).fractions
for der_idx in range(temp_order + 1):
name = "{0}{1}".format(name, "_" + "".join(["d" for x in range(der_idx)] + ["t"]) if der_idx > 0 else "")
data.append(evaluate_approximation(weight_lbl, q[:, der_idx * ini_funcs.size:(der_idx + 1) * ini_funcs.size],
temp_domain, spat_domain, name=name))
# spatial (0th derivative is skipped since this is already handled above)
for der_idx in range(1, spat_order + 1):
name = "{0}{1}".format(name, "_" + "".join(["d" for x in range(der_idx)] + ["z"]) if der_idx > 0 else "")
data.append(
evaluate_approximation(weight_lbl, q[:, :ini_funcs.size], temp_domain, spat_domain, der_idx, name=name))
return data
[docs]def get_sim_results(temp_domain, spat_domains, weights, state_space, names=None,
derivative_orders=None):
"""
Convenience wrapper for :py:func:`.get_sim_result`.
Args:
temp_domain (:py:class:`.Domain`): Time domain
spat_domains (dict): Spatial domain from all subsystems which belongs to
*state_space* as values and name of the systems as keys.
weights (numpy.array): Weights gained through simulation. For example
with :py:func:`.simulate_state_space`.
state_space (:py:class:`.StateSpace`): Simulated state space instance.
names: List of names of the desired systems. If not given all available
subssystems will be processed.
derivative_orders (dict): Desired derivative orders.
Returns:
List of :py:class:`.EvalData` objects.
"""
ss_base = get_base(state_space.base_lbl)
if names is None:
# TODO: implement getter method in StackedBase or change function interface
if isinstance(ss_base, StackedBase):
labels = [lbl for lbl in ss_base._info.keys()]
names = [ss_base._info[lbl]["sys_name"] for lbl in labels]
else:
names = list(spat_domains)
labels = [state_space.base_lbl]
else:
if isinstance(ss_base, StackedBase):
labels = list()
for nm in names:
labels = [key for key, val in ss_base._info.items()
if val["sys_name"] is nm]
else:
labels = [state_space.base_lbl]
if derivative_orders is None:
derivative_orders = dict([(name, (0, 0)) for name in names])
results = []
for nm, lbl in zip(names, labels):
# if derivative_orders[n] is None derivatives of the
# corresponding variables are not provided
if derivative_orders[nm][0] is None:
derivative_orders[nm][0] = 0
if derivative_orders[nm][1] is None:
derivative_orders[nm][1] = 0
# acquire a transformation into the original weights
info = get_transformation_info(state_space.base_lbl,
lbl,
int(weights.shape[1] / ss_base.fractions.size) - 1,
derivative_orders[nm][0])
transformation = get_weight_transformation(info)
# project back
data = get_sim_result(info.dst_lbl,
np.apply_along_axis(transformation, 1, weights),
temp_domain,
spat_domains[nm],
info.dst_order,
derivative_orders[nm][1],
name=nm)
results += data
return results
[docs]class CanonicalEquation(object):
"""
Wrapper object, holding several entities of canonical forms for different
weight-sets that form an equation when summed up.
After instantiation, this object can be filled with information by passing
the corresponding coefficients to :py:meth:`.add_to`. When the parsing
process is completed and all coefficients have been collected, calling
:py:meth:`.finalize` is required to compute all necessary information for
further processing. When finalized, this object provides access to the
dominant form of this equation.
Args:
name (str): Unique identifier of this equation.
dominant_lbl (str): Label of the variable that dominates this equation.
"""
def __init__(self, name, dominant_lbl=None):
self.name = name
self.dominant_lbl = dominant_lbl
self.dynamic_forms = {}
self._static_form = CanonicalForm(self.name + "_static")
self._finalized = False
self._finalized_dynamic_forms = False
[docs] def add_to(self, weight_label, term, val, column=None):
"""
Add the provided *val* to the canonical form for *weight_label*,
see :py:meth:`.CanonicalForm.add_to` for further information.
Args:
weight_label (str): Basis to add onto.
term: Coefficient to add onto, see :py:func:`~CanonicalForm.add_to`.
val: Values to add.
column (int): passed to :py:func:`~CanonicalForm.add_to`.
"""
if self._finalized:
raise RuntimeError("Object has already been finalized, you are trying some nasty stuff there.")
if term["name"] in "fG":
# hold f and g vector separately
self._static_form.add_to(term, val, column)
return
if weight_label is None:
raise ValueError("weight_label can only be none if target is f or G.")
if weight_label not in list(self.dynamic_forms.keys()):
self.dynamic_forms[weight_label] = CanonicalForm("_".join([self.name + weight_label]))
self.dynamic_forms[weight_label].add_to(term, val)
[docs] def finalize(self):
"""
Finalize the Object.
After the complete formulation has been parsed and all terms have been
sorted into this Object via :py:meth:`.add_to` this function has to be
called to inform this object about it. Furthermore, the f and G parts of
the static_form will be copied to the dominant form for easier
state-space transformation.
Note:
This function must be called to use the :py:attr:`dominant_form`
attribute.
"""
if self.dominant_lbl is None:
raise ValueError("You have to set the dominant labels of the\n"
"canonical equation (weak form), for example\n"
"with pyinduct.simulation.set_dominant_labels().")
if not self._finalized_dynamic_forms:
self.finalize_dynamic_forms()
if self.dynamic_forms[self.dominant_lbl].singular:
raise ValueError("The form that has to be chosen is singular.")
# copy static terms to dominant form to transform them correctly
for letter in "fG":
if letter in self._static_form.matrices:
self.dynamic_forms[self.dominant_lbl].matrices.update({letter: self._static_form.matrices[letter]})
self._finalized = True
@property
def static_form(self):
"""
:py:class:`.WeakForm` that does not depend on any weights.
:return:
"""
return self._static_form
@property
def dominant_form(self):
"""
direct access to the dominant :py:class:`.CanonicalForm`.
Note:
:py:meth:`.finalize` must be called first.
Returns:
:py:class:`.CanonicalForm`: the dominant canonical form
"""
if self.dominant_lbl is None:
raise RuntimeError("Dominant label is not defined! Use for\n"
"expample pyinduct.simulation."
"set_dominant_label or set it manually.")
return self.dynamic_forms[self.dominant_lbl]
[docs] def get_static_terms(self):
"""
Return:
Terms that do not depend on a certain weight set.
"""
return self._static_form.get_terms()
[docs] def get_dynamic_terms(self):
"""
Return:
dict: Dictionary of terms for each weight set.
"""
return {label: val.get_terms() for label, val in self.dynamic_forms.items()}
@property
def input_function(self):
"""
The input handle for the equation.
"""
return self._static_form.input_function
@input_function.setter
def input_function(self, func):
self._static_form.input_function = func
[docs]def create_state_space(canonical_equations):
"""
Create a state-space system constituted by several
:py:class:`.CanonicalEquations` (created by
:py:func:`.parse_weak_formulation`)
Args:
canonical_equations: List of :py:class:`.CanonicalEquation`'s.
Raises:
ValueError: If compatibility criteria cannot be fulfilled
Return:
:py:class:`.StateSpace`: State-space representation of the approximated
system
"""
set_dominant_labels(canonical_equations)
if isinstance(canonical_equations, CanonicalEquation):
# backward compatibility
canonical_equations = [canonical_equations]
# check whether the formulations are compatible
for eq in canonical_equations:
for lbl, form in eq.dynamic_forms.items():
coupling_order = form.max_temp_order
# search corresponding dominant form in other equations
for _eq in canonical_equations:
# check uniqueness of name - dom_lbl mappings
if eq.name != _eq.name and eq.dominant_lbl == _eq.dominant_lbl:
raise ValueError("A dominant form has to be unique over all given Equations")
# identify coupling terms
if lbl == eq.dominant_lbl:
break
# identify corresponding dominant form
if _eq.dominant_lbl != lbl:
continue
dominant_order = _eq.dominant_form.max_temp_order
if dominant_order <= coupling_order:
# dominant order has to be at least one higher than
# the coupling order
raise ValueError("Formulations are not compatible")
# transform dominant forms into state-space representation
# and collect information
dominant_state_spaces = {}
state_space_props = Parameters(size=0,
parts=OrderedDict(),
powers=set(),
input_powers=set(),
dim_u=0,
input=None)
for eq in canonical_equations:
dom_lbl = eq.dominant_lbl
dom_form = eq.dominant_form
dom_ss = dom_form.convert_to_state_space()
dominant_state_spaces.update({dom_lbl: dom_ss})
# collect some information
state_space_props.parts[dom_lbl] = dict(start=copy(state_space_props.size),
orig_size=dom_form.dim_x,
size=dom_form.dim_xb,
order=dom_form.max_temp_order - 1,
sys_name=eq.name)
state_space_props.powers.update(dom_form.powers)
state_space_props.size += dom_form.dim_xb
state_space_props.dim_u = max(state_space_props.dim_u, dom_form.dim_u)
# update input handles
if state_space_props.input is None:
state_space_props.input = eq.input_function
else:
if eq.input_function is not None and state_space_props.input != eq.input_function:
raise ValueError("Only one input object allowed.")
# build new basis by concatenating the dominant bases of every equation
if len(canonical_equations) == 1:
new_name = next(iter(canonical_equations)).dominant_lbl
else:
members = state_space_props.parts.keys()
new_name = "_".join(members)
fracs = [frac for lbl in members for frac in get_base(lbl).fractions]
new_base = StackedBase(fracs, state_space_props.parts)
register_base(new_name, new_base)
# build new state transition matrices A_p_k for corresponding powers p_k of the state vector
a_matrices = {}
for p in state_space_props.powers:
a_mat = np.zeros((state_space_props.size, state_space_props.size))
for row_eq in canonical_equations:
row_dom_lbl = row_eq.dominant_lbl
row_dom_dim = state_space_props.parts[row_dom_lbl]["size"]
row_dom_trans_mat = row_eq.dominant_form.e_n_pb_inv
row_dom_sys_mat = dominant_state_spaces[row_dom_lbl].A.get(p, None)
row_idx = state_space_props.parts[row_dom_lbl]["start"]
for col_eq in canonical_equations:
col_dom_lbl = col_eq.dominant_lbl
# main diagonal
if col_eq.name == row_eq.name:
if row_dom_sys_mat is not None:
a_mat[row_idx:row_idx + row_dom_dim, row_idx:row_idx + row_dom_dim] = row_dom_sys_mat
continue
# coupling terms
if col_dom_lbl in row_eq.dynamic_forms:
for order, mats in row_eq.dynamic_forms[col_dom_lbl].matrices["E"].items():
orig_mat = mats.get(p, None)
if orig_mat is not None:
# transform matrix with row-transformation matrix and add to last "row"
# since it's not the dominant entry, revert sign change
def parse_weak_formulation(weak_form, finalize=False):
r"""
Parses a :py:class:`.WeakFormulation` that has been derived by projecting a
partial differential equation an a set of test-functions. Within this
process, the separating approximation
:math:`x^n(z, t) = \sum_{i=1}^n c_i^n(t) \varphi_i^n(z)` is plugged into
the equation and the separated spatial terms are evaluated, leading to a
ordinary equation system for the weights :math:`c_i^n(t)`.
Args:
weak_form: Weak formulation of the pde.
finalize (bool): Default: False. If you have already defined the
dominant labels of the weak formulations you can set this to True.
See :py:meth:`.CanonicalEquation.finalize`
Return:
:py:class:`.CanonicalEquation`: The spatially approximated equation in
a canonical form.
"""
if not isinstance(weak_form, WeakFormulation):
raise TypeError("Only able to parse WeakFormulation")
ce = CanonicalEquation(weak_form.name, weak_form.dominant_lbl)
# handle each term
for term in weak_form.terms:
# extract Placeholders
placeholders = dict(scalars=term.arg.get_arg_by_class(Scalars),
functions=term.arg.get_arg_by_class(TestFunction),
field_variables=term.arg.get_arg_by_class(FieldVariable),
inputs=term.arg.get_arg_by_class(Input))
# field variable terms: sort into E_np, E_n-1p, ..., E_0p
if placeholders["field_variables"]:
if len(placeholders["field_variables"]) != 1:
raise NotImplementedError
field_var = placeholders["field_variables"][0]
if not field_var.simulation_compliant:
raise ValueError("Shape- and test-function labels of FieldVariable must match for simulation purposes.")
temp_order = field_var.order[0]
exponent = field_var.data["exponent"]
term_info = dict(name="E", order=temp_order, exponent=exponent)
base = get_base(field_var.data["func_lbl"]).derive(field_var.order[1])
shape_funcs = base.raise_to(exponent)
if placeholders["inputs"]:
# essentially, this means that parts of the state-transition
# matrix will be time dependent
raise NotImplementedError
if placeholders["functions"]:
# is the integrand a product?
if len(placeholders["functions"]) != 1:
raise NotImplementedError
func = placeholders["functions"][0]
test_funcs = get_base(func.data["func_lbl"]).derive(func.order[1])
result = calculate_scalar_product_matrix(dot_product_l2,
test_funcs,
shape_funcs)
else:
# extract constant term and compute integral
# TODO this is a source of complex data, since integrate
# function will return complex dtype.
a = Scalars(np.atleast_2d(
[integrate_function(func, func.nonzero)[0]
for func in shape_funcs.fractions]))
if placeholders["scalars"]:
b = placeholders["scalars"][0]
result = _compute_product_of_scalars([a, b])
else:
result = a.data
ce.add_to(weight_label=field_var.data["weight_lbl"],
term=term_info,
val=result * term.scale)
continue
# TestFunctions or pre evaluated terms, those can end up in E, f or G
if placeholders["functions"]:
if not 1 <= len(placeholders["functions"]) <= 2:
raise NotImplementedError
func = placeholders["functions"][0]
test_funcs = get_base(func.data["func_lbl"]).derive(func.order[1]).fractions
if len(placeholders["functions"]) == 2:
# TODO this computation is nonsense. Result must be a vector containing int(tf1*tf2)
raise NotImplementedError
#
# func2 = placeholders["functions"][1]
# test_funcs2 = get_base(func2.data["func_lbl"], func2.order[2])
# result = calculate_scalar_product_matrix(dot_product_l2, test_funcs, test_funcs2)
# cf.add_to(("f", 0), result * term.scale)
# continue
if placeholders["scalars"]:
a = placeholders["scalars"][0]
b = Scalars(np.vstack([integrate_function(func, func.nonzero)[0] for func in test_funcs]))
result = _compute_product_of_scalars([a, b])
ce.add_to(weight_label=a.target_form,
term=get_common_target(placeholders["scalars"]),
val=result * term.scale)
continue
if placeholders["inputs"]:
if len(placeholders["inputs"]) != 1:
raise NotImplementedError
input_var = placeholders["inputs"][0]
input_func = input_var.data["input"]
input_index = input_var.data["index"]
input_exp = input_var.data["exponent"]
input_order = input_var.order[0]
term_info = dict(name="G", order=input_order, exponent=input_exp)
result = np.array([[integrate_function(func, func.nonzero)[0]] for func in test_funcs])
ce.add_to(weight_label=None, term=term_info, val=result * term.scale, column=input_index)
ce.input_function = input_func
continue
# pure scalar terms, sort into corresponding matrices
if placeholders["scalars"]:
result = _compute_product_of_scalars(placeholders["scalars"])
target = get_common_target(placeholders["scalars"])
target_form = get_common_form(placeholders)
if placeholders["inputs"]:
input_var = placeholders["inputs"][0]
input_func = input_var.data["input"]
input_index = input_var.data["index"]
input_exp = input_var.data["exponent"]
input_order = input_var.order[0]
term_info = dict(name="G", order=input_order, exponent=input_exp)
if input_order > 0:
# here we would need to provide derivative handles in the callable
raise NotImplementedError
if target["name"] == "E":
# this would mean that the input term should appear in a matrix like E1 or E2
# the result would be a time dependant sate transition matrix
raise NotImplementedError
ce.add_to(weight_label=None, term=term_info, val=result * term.scale, column=input_index)
ce.input_function = input_func
continue
ce.add_to(weight_label=target_form, term=target, val=result * term.scale)
continue
# inform object that the parsing process is complete
if finalize:
ce.finalize()
return ce
def _compute_product_of_scalars(scalars):
if len(scalars) > 2:
raise NotImplementedError
if len(scalars) == 1:
# simple scaling of all terms
res = scalars[0].data
elif scalars[0].data.shape == scalars[1].data.shape:
# element wise multiplication
res = np.prod(np.array([scalars[0].data, scalars[1].data]), axis=0)
elif scalars[0].data.shape == (1, 1) or scalars[1].data.shape == (1, 1):
# a lumped terms is present
res = scalars[0].data * scalars[1].data
else:
# dyadic product
try:
if scalars[0].data.shape[1] == 1:
[docs] res = scalars[0].data @ scalars[1].data
else:
res = scalars[1].data @ scalars[0].data
except ValueError as e:
raise ValueError("provided entries do not form a dyadic product")
return res
def simulate_state_space(state_space, initial_state, temp_domain, settings=None):
r"""
Wrapper to simulate a system given in state space form:
.. math:: \dot{q} = A_pq^p + A_{p-1}q^{p-1} + \dotsb + A_0q + Bu.
Args:
state_space (:py:class:`.StateSpace`): State space formulation of the
system.
initial_state: Initial state vector of the system.
temp_domain (:py:class:`.Domain`): Temporal domain object.
settings (dict): Parameters to pass to the :py:func:`set_integrator`
method of the :class:`scipy.ode` class, with the integrator name
included under the key :obj:`name`.
Return:
tuple: Time :py:class:`.Domain` object and weights matrix.
"""
if not isinstance(state_space, StateSpace):
raise TypeError
input_handle = state_space.input
if not isinstance(input_handle, SimulationInput):
raise TypeError("only simulation.SimulationInput supported.")
q = [initial_state]
t = [temp_domain[0]]
r = ode(state_space.rhs)
# TODO check for complex-valued matrices and use 'zvode'
if settings:
r.set_integrator(settings.pop("name"), **settings)
else:
# use some sane defaults
r.set_integrator(
"vode",
max_step=temp_domain.step,
method="adams",
nsteps=1e3
)
r.set_initial_value(q[0], t[0])
for t_step in temp_domain[1:]:
qn = r.integrate(t_step)
if not r.successful():
warnings.warn("*** Error: Simulation aborted at t={} ***".format(r.t))
break
t.append(r.t)
q.append(qn)
# create results
q = np.array(q)
return Domain(points=np.array(t), step=temp_domain.step), q
[docs]def evaluate_approximation(base_label, weights, temp_domain, spat_domain, spat_order=0, name=""):
"""
Evaluate an approximation given by weights and functions at the points given
in spatial and temporal steps.
Args:
weights: 2d np.ndarray where axis 1 is the weight index and axis 0 the
temporal index.
base_label (str): Functions to use for back-projection.
temp_domain (:py:class:`.Domain`): For steps to evaluate at.
spat_domain (:py:class:`.Domain`): For points to evaluate at (or in).
spat_order: Spatial derivative order to use.
name: Name to use.
Return:
:py:class:`.EvalData`
"""
funcs = get_base(base_label).derive(spat_order).fractions
if weights.shape[1] != funcs.shape[0]:
raise ValueError("weights (len={0}) have to fit provided functions (len={1})!".format(weights.shape[1],
funcs.size))
# evaluate shape functions at given points
shape_vals = np.array([func.evaluation_hint(spat_domain) for func in funcs])
def eval_spatially(weight_vector):
return np.real_if_close(np.dot(weight_vector, shape_vals), 1000)
data = np.apply_along_axis(eval_spatially, 1, weights)
return EvalData([temp_domain.points, spat_domain.points], data, name=name)
[docs]def set_dominant_labels(canonical_equations, finalize=True):
"""
Set the dominant label (*dominant_lbl*) member of all given canonical
equations and check if the problem formulation is valid (see background
section: http://pyinduct.readthedocs.io/en/latest/).
If the dominant label of one or more :py:class:`.CanonicalEquation`
is already defined, the function raise a UserWarning if the (pre)defined
dominant label(s) are not valid.
Args:
canonical_equations: List of :py:class:`.CanonicalEquation` instances.
finalize (bool): Finalize the equations? Default: True.
"""
if isinstance(canonical_equations, CanonicalEquation):
canonical_equations = [canonical_equations]
# collect all involved labels
labels = set(
chain(*[list(ce.dynamic_forms.keys()) for ce in canonical_equations]))
if len(labels) != len(canonical_equations):
raise ValueError("The N defined canonical equations (weak forms)\n"
"must hold exactly N different weight labels!\n"
"But your {} canonical equation(s) (weak form(s))\n"
"hold {} weight label(s)!"
"".format(len(canonical_equations),
len(labels)))
max_orders = dict()
for ce in canonical_equations:
ce.finalize_dynamic_forms()
for lbl in list(ce.dynamic_forms.keys()):
max_order = dict(
(("max_order", ce.dynamic_forms[lbl].max_temp_order),
("can_eqs", [ce])))
if lbl not in max_orders or \
max_orders[lbl]["max_order"] < max_order["max_order"]:
max_orders[lbl] = max_order
elif max_orders[lbl]["max_order"] == max_order["max_order"]:
max_orders[lbl]["can_eqs"].append(
max_order["can_eqs"][0])
non_valid1 = [(lbl, max_orders[lbl])
for lbl in labels if len(max_orders[lbl]["can_eqs"]) > 1]
if non_valid1:
raise ValueError("The highest time derivative from a certain weight\n"
"label may only occur in one canonical equation. But\n"
"each of the canonical equations {} holds the\n"
"weight label '{}' with order {} in time."
"".format(non_valid1[0][1]["can_eqs"][0].name,
non_valid1[0][0],
non_valid1[0][1]["max_order"]))
non_valid2 = [lbl for lbl in labels if max_orders[lbl]["max_order"] == 0]
if non_valid2:
raise ValueError("The defined problem leads to an differential\n"
"algebraic equation, since there is no time\n"
"derivative for the weights {}. Such problems are\n"
"not considered in pyinduct, yet."
"".format(non_valid2))
# set/check dominant labels
for lbl in labels:
pre_lbl = max_orders[lbl]["can_eqs"][0].dominant_lbl
max_orders[lbl]["can_eqs"][0].dominant_lbl = lbl
if pre_lbl is not None and pre_lbl != lbl:
warnings.warn("\n Predefined dominant label '{}' from\n"
"canonical equation / weak form '{}' not valid!\n"
"It will be overwritten with the label '{}'."
"".format(pre_lbl,
max_orders[lbl]["can_eqs"][0].name,
lbl),
UserWarning)
if finalize:
for ce in canonical_equations:
ce.finalize()