<commit_before># -*- coding: UTF-8 -*-
"""
:mod:`System` -- Systems
============================

Examples of ODE systems.

Collection of subclasses of the class :class:`odelab.system.System`.
The interface needed depends on the solver. For example, exponential solvers will need a `linear` method,
dae solvers will need a `constraint`, or `multi-dynamics` method. See the documentation of the :doc:`scheme`
section for more information.

.. module :: system
	:synopsis: Examples of ODE systems.
.. moduleauthor :: <NAME> <<EMAIL>>

"""
from __future__ import division

import numpy as np
from numpy import array

import scipy.sparse as sparse

import odelab.scheme.rungekutta as rk

class System(object):
	"""
	General System class to define a simple dynamics given by a right-hand side.
	"""
	def __init__(self, f=None):
		if f is not None:
			self.f = f

	def label(self, component):
		return '%d' % component

	def preprocess(self, u0):
		return u0

	def postprocess(self, u1):
		return u1

	def velocity(self,u):
		"""select the velocity component"""

	def position(self, u):
		"""select the position component"""

	def lag(self, u):
		"""select the lagrangian component"""

	def vel_lag_split(self, vl):
		"""splits between velocity and lagrangian"""

	def vel_lag_stack(self,v,l):
		"""stacks together the velocity and lagrangian parts"""

	def force(self, u):
		"""compute the force as a function of `u`"""

	def codistribution(self, u):
		"""compute the codistribution matrix at `u`"""


class JayExample(System):
	r"""
	 The example in [jay06]_ §5. This is a test to check implementation of the
	 SRK-DAE2 methods given in [jay06]_. We want to compare our results to
	 [jay06]_ Fig. 1.

	 The exact solution to this problem is known as is

.. math::

	y1(t) = \ee^t\\
	y2(t) = \ee^{-2t}\\
	z1(t) = \ee^{2t}

We will compute the global error at :math:`t=1` at plot this relative to the
stepsize :math:`h`. This is what is done in [jay06]_ Fig.1.


:References:

.. [jay06] Jay - *Specialized Runge-Kutta methods for index $ 2$ differential-algebraic equations.* Math. Comput. 75, No. 254, 641-654 (2006). :doi:`10.1090/S0025-5718-05-01809-0`
	"""


	def multi_dynamics(self, tu):
		y1,y2,z,t = tu
		return {
			rk.LobattoIIIA: array([y2 - 2*y1**2*y2, -y1**2]),
			rk.LobattoIIIB: array([y1*y2**2*z**2, np.exp(-t)*z - y1]),
			rk.LobattoIIIC: array([-y2**2*z, -3*y2**2*z]),
			rk.LobattoIIICs: array([2*y1*y2**2 - 2*np.exp(-2*t)*y1*y2, z]),
			rk.LobattoIIID: array([2*y2**2*z**2, y1**2*y2**2])
			}

	def constraint(self, tu):
		return array([tu[0]**2*tu[1] -  1])

	def state(self, u):
		return u[0:2]

	def lag(self, u):
		return u[2:3]

	@classmethod
	def exact(self, t, u0):
		if np.allclose(u0[:2], array([1.,1.])):
			return array([np.exp(t), np.exp(-2*t), np.exp(2*t)])
		raise ValueError("Exact solution not defined")

	def label(self, component):
		return ['y1','y2','z'][component]

class GraphSystem(System):
	r"""
	Trivial semi-explicit index 2 DAE of the form:

.. math::

		x' = 1\\
		y' = λ\\
		y  = f(x)
	"""

	def state(self, u):
		return u[:2]

	def lag(self, u):
		return u[-1:]

	def multi_dynamics(self, t, u):
		x,y = self.state(u)
		return {
			rk.LobattoIIIA: np.zeros_like(self.state(u)),
			rk.LobattoIIIB: array([np.ones_like(x), self.lag(u)[0]]),
			}

	def constraint(self, t, u):
		x,y = self.state(u)
		return array([y - self.f(x)])

	def hidden_error(self, t, u):
		return self.lag(u)[0] - self.f.der(t,self.state(u)[0])

class ODESystem(System):
	"""
	Simple wrapper to transform an ODE into a semi-explicit DAE.
	"""
	def __init__(self, f, RK_class=rk.LobattoIIIA):
		self.f = f
		self.RK_class = RK_class

	def state(self, u):
		return u[:1]

	def lag(self,u):
		return u[1:] # should be empty
<commit_msg>Fix the GraphSystem DAE system<commit_after># -*- coding: UTF-8 -*-
"""
:mod:`System` -- Systems
============================

Examples of ODE systems.

Collection of subclasses of the class :class:`odelab.system.System`.
The interface needed depends on the solver. For example, exponential solvers will need a `linear` method,
dae solvers will need a `constraint`, or `multi-dynamics` method. See the documentation of the :doc:`scheme`
section for more information.

.. module :: system
	:synopsis: Examples of ODE systems.
.. moduleauthor :: <NAME> <<EMAIL>>

"""
from __future__ import division

import numpy as np
from numpy import array

import scipy.sparse as sparse

import odelab.scheme.rungekutta as rk

class System(object):
	"""
	General System class to define a simple dynamics given by a right-hand side.
	"""
	def __init__(self, f=None):
		if f is not None:
			self.f = f

	def label(self, component):
		return '%d' % component

	def preprocess(self, u0):
		return u0

	def postprocess(self, u1):
		return u1

	def velocity(self,u):
		"""select the velocity component"""

	def position(self, u):
		"""select the position component"""

	def lag(self, u):
		"""select the lagrangian component"""

	def vel_lag_split(self, vl):
		"""splits between velocity and lagrangian"""

	def vel_lag_stack(self,v,l):
		"""stacks together the velocity and lagrangian parts"""

	def force(self, u):
		"""compute the force as a function of `u`"""

	def codistribution(self, u):
		"""compute the codistribution matrix at `u`"""


class JayExample(System):
	r"""
	 The example in [jay06]_ §5. This is a test to check implementation of the
	 SRK-DAE2 methods given in [jay06]_. We want to compare our results to
	 [jay06]_ Fig. 1.

	 The exact solution to this problem is known as is

.. math::

	y1(t) = \ee^t\\
	y2(t) = \ee^{-2t}\\
	z1(t) = \ee^{2t}

We will compute the global error at :math:`t=1` at plot this relative to the
stepsize :math:`h`. This is what is done in [jay06]_ Fig.1.


:References:

.. [jay06] Jay - *Specialized Runge-Kutta methods for index $ 2$ differential-algebraic equations.* Math. Comput. 75, No. 254, 641-654 (2006). :doi:`10.1090/S0025-5718-05-01809-0`
	"""


	def multi_dynamics(self, tu):
		y1,y2,z,t = tu
		return {
			rk.LobattoIIIA: array([y2 - 2*y1**2*y2, -y1**2]),
			rk.LobattoIIIB: array([y1*y2**2*z**2, np.exp(-t)*z - y1]),
			rk.LobattoIIIC: array([-y2**2*z, -3*y2**2*z]),
			rk.LobattoIIICs: array([2*y1*y2**2 - 2*np.exp(-2*t)*y1*y2, z]),
			rk.LobattoIIID: array([2*y2**2*z**2, y1**2*y2**2])
			}

	def constraint(self, tu):
		return array([tu[0]**2*tu[1] -  1])

	def state(self, u):
		return u[0:2]

	def lag(self, u):
		return u[2:3]

	@classmethod
	def exact(self, t, u0):
		if np.allclose(u0[:2], array([1.,1.])):
			return array([np.exp(t), np.exp(-2*t), np.exp(2*t)])
		raise ValueError("Exact solution not defined")

	def label(self, component):
		return ['y1','y2','z'][component]

class GraphSystem(System):
	r"""
	Trivial semi-explicit index 2 DAE of the form:

.. math::

		x' = 1\\
		y' = λ\\
		y  = f(x)
	"""

	def state(self, u):
		return u[:2]

	def lag(self, u):
		return u[2:3]

	def multi_dynamics(self, ut):
		x,y = self.state(ut)
		return {
			rk.RadauIIA: np.zeros_like(self.state(ut)),
			rk.RadauIIA: array([np.ones_like(x), self.lag(ut)[0]]),
			}

	def constraint(self, ut):
		x,y = self.state(ut)
		return array([y - self.f(x)])

	def hidden_error(self, t, u):
		return self.lag(u)[0] - self.f.der(t,self.state(u)[0])

	def exact(self, t, u0):
		x = u0[0] + t
		return array([x, self.f(x), self.f.der(x)])

class ODESystem(System):
	"""
	Simple wrapper to transform an ODE into a semi-explicit DAE.
	"""
	def __init__(self, f, RK_class=rk.LobattoIIIA):
		self.f = f
		self.RK_class = RK_class

	def state(self, u):
		return u[:1]

	def lag(self,u):
		return u[1:] # should be empty

{'fraction_non_alphanumeric': 0.17304908540597447, 'fraction_numerical': 0.03316225819186735, 'mean_word_length': 3.840764331210191, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7272620', 'n_tokens_mistral': 3214, 'n_tokens_neox': 2891, 'n_words': 855}