# Source code for ahkab.trap

# -*- coding: iso-8859-1 -*-
# trap.py
# Trap DF (with a second order prediction)

# This file is part of the ahkab simulator.
#
# Ahkab is free software: you can redistribute it and/or modify
# the Free Software Foundation, version 2 of the License.
#
# Ahkab is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License v2
# along with ahkab.  If not, see <http://www.gnu.org/licenses/>.

"""
This file implements the Trapezoidal (TRAP) Differentiation Formula (DF) and a
second order prediction formula.

Module reference
----------------
"""

from __future__ import (unicode_literals, absolute_import,
division, print_function)

import numpy as np
from numpy.linalg import inv

from .py3compat import range_type

order = 2

[docs]def is_implicit():
"""Is this Differentiation Formula (DF) implicit?

**Returns:**

isit : boolean
In this case, that's True.
"""
return True

[docs]def get_required_values():
"""Get info regarding what values are needed by the DF

**Returns:**

tpl : tuple of tuples
A tuple of two tuples.

* The first tuple indicates what past values of the unknown are needed
for the DF.
* The second tuple indicates what past values of the unknown are needed
for the prediction method.

In particular, each of the sub-tuples is built this way:
::

(max_order_of_x, max_order_of_dx)

Where both the values are either int, or None.
If max_order_of_x is set to an integer value :math:k, the DF needs
all the :math:x_{n-i} values of x, for all :math:0 \\le i \\le k. In
the previous text, :math:x_{n-i} is the value the :math:x array
assumed :math:i steps before the one we are considering for the
derivative.

Similar considerations apply to max_order_of_dx, but regard rather
:math:dx_n/dt instead of :math:x_n.

If any of the values is set to None, it is to be assume that no
value is required.

The first array has to be used if no prediction is required, the second are the values needed for prediction.
"""
# we need x(n) and dx(n)/dt
return ((0, 0), (2, None))

[docs]def has_ff():
"""Has the method a Forward Formula for prediction?

**Returns:**

doesit : bool
In this particular case, this function always returns True.
"""
return True

[docs]def get_df(pv_array, suggested_step, predict=True):
"""Get the coefficients for DF and prediction formula

**Parameters:**

pv_array : sequence of sequences
Each element of pv_array must be of the form:

::

(time, x, derivate(x))

In particular, the :math:k element of pv_array contains the values
of:

* :math:t_{n-k} (the time),
* :math:x_{n-k},
* :math:dx_{n-k}/dt

evaluated :math:k time steps before the current one, labeled
:math:n+1.

How many samples are necessary is given by
:func:ahkab.trap.get_required_values.

Values that are not needed may be set to None, as they will be
disregarded.
suggested_step : float
The step that will be used for the current iteration, provided the error
will be deemed acceptable.
predict : bool, optional
Whether a prediction for :math:x_n is needed as well or not. Defaults
to True.

**Returns:**

ret : tuple
The return value has the form:

::

(C1, C0, x_lte_coeff, predict_x, predict_lte_coeff)

The derivative may be written as:

.. math::

d(x(n+1))/dt = C1 x(n+1) + C0

x_lte_coeff is the coefficient of the Local Truncation Error,
predict_x is the predicted value for :math:x and predict_lte_coeff
is the LTE coefficient for the prediction.

:raises ValueError: if the pv_array is malformed.
"""

# our method needs x(n) dx(n)/dt
if len(pv_array[0]) != 3:
raise ValueError('trap.get_df got a pv_array of wrong dimensions')
if pv_array[0][1] is None or pv_array[0][2] is None:
raise ValueError('trap.get_df got a pv_array that does not contain'
' required values')
if predict and (len(pv_array) < 2 or pv_array[0][1] is None or
pv_array[1][1] is None or pv_array[2][1] is None):
raise ValueError('trap.get_df got predict=True but pv_array that does'
' not contain the required values')

C1 = 2.0 / suggested_step
C0 = -1.0 * (2.0 / suggested_step * pv_array[0][1] + pv_array[0][2])
x_lte_coeff = 1.0 / 12 * suggested_step ** 3

if predict and len(pv_array) > 2 and pv_array[0][1] is not None and \
pv_array[1][1] is not None and pv_array[2][1] is not None:
A = np.zeros((len(pv_array[0]), len(pv_array[0])))
A[0, :] = 0
A[:, 0] = 1
for row in range_type(1, A.shape[0]):
for col in range_type(1, A.shape[0]):
A[row, col] = (pv_array[row][0] - pv_array[0][0]) ** col
Ainv = inv(A)

predict_x = np.zeros(pv_array[0][1].shape)
z = np.zeros((3, 1))
for var in range_type(pv_array[0][1].shape[0]):
for index in range(z.shape[0]):
z[index, 0] = pv_array[index][1][var, 0]
alpha = np.dot(Ainv, z)
predict_x[var, 0] = alpha[2, 0] * suggested_step**2 + \
alpha[1, 0] * suggested_step + alpha[0, 0]
predict_lte_coeff = -1.0 / 6.0 * suggested_step * \
(pv_array[0][0] + suggested_step - pv_array[1][0]) * \
(pv_array[0][0] + suggested_step - pv_array[2][0])
else:
predict_x, predict_lte_coeff = (None, None)
return C1, C0, x_lte_coeff, predict_x, predict_lte_coeff