Source code for simulai.utilities.opinf_deviation

# (C) Copyright IBM Corp. 2019, 2020, 2021, 2022.

#    Licensed under the Apache License, Version 2.0 (the "License");
#    you may not use this file except in compliance with the License.
#    You may obtain a copy of the License at

#           http://www.apache.org/licenses/LICENSE-2.0

#     Unless required by applicable law or agreed to in writing, software
#     distributed under the License is distributed on an "AS IS" BASIS,
#     WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#     See the License for the specific language governing permissions and
#     limitations under the License.

import os
import pickle
import numpy as np
import sympy as sp
from sympy import MatrixExpr

[docs]class OpInfDeviation: def __init__(self, A_hat: np.ndarray = None, H_hat: np.ndarray = None) -> None: """Evaluating the deviation evolution in an OpInf model :param A_hat: linear OpInf operator :type A_hat: np.ndarray :param H_hat: quadratic OpInf operator :type: np.ndarray :return: nothing """ self.A_hat = A_hat self.H_hat = H_hat self.n = A_hat.shape[0] epsilon = sp.MatrixSymbol('epsilon', self.n, 1) u = sp.MatrixSymbol('u', self.n, 1) # u*e.T ue = sp.MatMul(u, epsilon.T) # (u*e.T).T ue_T = ue.T # e*e.T ee = sp.MatMul(epsilon, epsilon.T) # U(u*e.T + (u*e.T).T) + e X e v_u = np.array(ue + ue_T + ee)[np.triu_indices(self.n)] v = sp.Matrix(v_u) H_symb = sp.Matrix(self.H_hat) A_symb = sp.Matrix(self.A_hat) v_jac = sp.Matrix.jacobian(H_symb @ v, epsilon) # A + H*(U(u*e.T + (u*e.T).T) + e X e) self.jac_expressions = A_symb + v_jac self.error_expressions = A_symb @ epsilon + H_symb @ v self.epsilon = epsilon self.u = u self.jac = None self.error = None self.compile()
[docs] def compile(self) -> None: self.jac = self.lambdify(expression=self.jac_expressions) self.error = self.lambdify(expression=self.error_expressions)
[docs] def lambdify(self, expression:MatrixExpr=None) -> callable: return sp.lambdify([self.epsilon, self.u], expression, 'numpy')
[docs] def eval_jacobian(self, u:np.ndarray=None, epsilon:np.ndarray=None) -> np.ndarray: """Evaluating error Jacobian :param u: reference solution :type u: np.ndarray :param epsilon: error associated to u :type epsilon: np.ndarray :return: error Jacobian :rtype: np.ndarray """ u = u.T epsilon = epsilon.T return self.jac(epsilon, u)
[docs] def eval_error(self, u:np.ndarray=None, epsilon:np.array=None) -> np.ndarray: """Evaluating error :param u: reference solution :type u: np.ndarray :param epsilon: error associated to u :type epsilon: np.ndarray :return: error :rtype: np.ndarray """ u = u.T epsilon = epsilon.T return self.error(epsilon, u)
[docs] def save(self, name:str=None, path:str=None) -> None: """Complete saving :param path: path to the saving directory :type path: str :param name: name for the model :type name: str :return: nothing """ blacklist = ['jac', 'error'] for item in blacklist: delattr(self, item) with open(os.path.join(path, name + '.pkl'), 'wb') as fp: pickle.dump(self, fp, protocol=4)