Fix merge

Author: Yaroslav-Muravev

epde/solver/base.py

@@ -0,0 +1,51 @@
+from abc import ABC, abstractmethod
+from typing import List, Union, Tuple, Dict, Any
+import numpy as np
+from epde.structure.main_structures import Equation, SoEq
+
+class BaseSolverAdapter(ABC):
+    """Базовый класс для всех адаптеров солверов."""
+
+    @abstractmethod
+    def solve(self, equation_or_system: Union[Equation, SoEq],
+              grids: List[np.ndarray],
+              data: Union[np.ndarray, List[np.ndarray]]) -> Tuple[List[np.ndarray], float]:
+        """
+        Решает уравнение/систему на заданных сетках.
+
+        Parameters
+        ----------
+        equation_or_system : Equation or SoEq
+            Одиночное уравнение или система.
+        grids : list of np.ndarray
+            Сетки координат (каждая размерность – отдельный массив).
+        data : np.ndarray or list of np.ndarray
+            Эталонные значения (для каждой переменной) для расчёта ошибки.
+
+        Returns
+        -------
+        solutions : list of np.ndarray
+            Список решений (по одному массиву на переменную).
+        loss : float
+            Числовая метрика ошибки (например, RMSE по всей сетке).
+        """
+        pass
+
+    @abstractmethod
+    def get_requirements(self) -> Dict[str, Any]:
+        """
+        Возвращает требования солвера к представлению уравнений.
+
+        Returns
+        -------
+        dict
+            Ключи:
+            - 'form': 'explicit_ode', 'pde_residual', 'system_of_odes'
+            - 'needs_initial_conditions': bool
+            - 'needs_boundary_conditions': bool
+            - 'supports_multisample': bool
+        """
+        pass
+
+    def supports_multisample(self) -> bool:
+        return self.get_requirements().get('supports_multisample', False)

epde/solver/classical_ode_adapter.py

@@ -0,0 +1,106 @@
+import numpy as np
+from scipy.integrate import solve_ivp
+from epde.solver.base import BaseSolverAdapter
+from epde.solver.ode_converter import ODEToFirstOrder
+from epde.structure.main_structures import Equation, SoEq
+import epde.globals as global_var
+import sympy as sp
+
+
+class ClassicalODEAdapter(BaseSolverAdapter):
+    def __init__(self, **config):
+        self.method = config.get('method', 'RK45')
+        self.rtol = config.get('rtol', 1e-6)
+        self.atol = config.get('atol', 1e-9)
+        self.rhs = config.get('rhs')
+        self.y0 = config.get('y0')
+        self._auto_rhs = self.rhs is None
+        self._converter = ODEToFirstOrder() if self._auto_rhs else None
+
+    def _get_order(self, equation):
+        if not self._auto_rhs:
+            return None
+        expr = self._converter._build_expression(equation, equation.main_var_to_explain)
+        u = sp.Function(equation.main_var_to_explain)(self._converter.t)
+        max_order = 0
+        for term in sp.preorder_traversal(expr):
+            if isinstance(term, sp.Derivative) and term.args[0] == u:
+                order = term.args[1][1]
+                if order > max_order:
+                    max_order = order
+        return max_order
+
+    def _get_initial_conditions(self, equation, data):
+        if self.y0 is not None:
+            return self.y0
+        order = self._get_order(equation)
+        if order == 1:
+            return [data[0]]
+        else:
+            raise ValueError(f"Для уравнения порядка {order} необходимо явно задать y0 в конфигурации.")
+
+    def solve(self, equation_or_system, grids, data):
+        print("[DEBUG] ClassicalODEAdapter.solve called")
+        print(f"[DEBUG] equation_or_system type: {type(equation_or_system)}")
+        print(f"[DEBUG] data type: {type(data)}")
+        if isinstance(data, (list, tuple)):
+            print(f"[DEBUG] data length: {len(data)}")
+            for i, d in enumerate(data):
+                print(f"[DEBUG] data[{i}].shape: {d.shape}")
+        else:
+            print(f"[DEBUG] data.shape: {data.shape}")
+
+        if len(grids) != 1:
+            raise ValueError("ClassicalODEAdapter работает только с 1D временной сеткой.")
+        mask = global_var.grid_cache.g_func_mask
+        t_full = grids[0].flatten()
+        t_masked = t_full[mask]
+        t_span = (t_masked.min(), t_masked.max())
+        t_eval = t_masked
+
+        # Определяем правую часть и начальные условия
+        if isinstance(equation_or_system, Equation):
+            if self.rhs is not None:
+                rhs = self.rhs
+                y0 = self.y0
+                if y0 is None:
+                    raise ValueError("Для явной rhs необходимо указать y0.")
+            else:
+                rhs = self._converter.equation_to_rhs(equation_or_system, equation_or_system.main_var_to_explain)
+                y0 = self._get_initial_conditions(equation_or_system, data)
+        elif isinstance(equation_or_system, SoEq):
+            # Для системы: rhs и y0 должны быть заданы явно
+            if self.rhs is None:
+                raise NotImplementedError(
+                    "Автоматическое преобразование систем ОДУ пока не поддерживается. Укажите rhs в конфигурации.")
+            rhs = self.rhs
+            y0 = self.y0
+            if y0 is None:
+                raise ValueError("Для системы необходимо указать y0 в конфигурации.")
+        else:
+            raise TypeError("Unsupported equation type")
+
+        sol = solve_ivp(rhs, t_span, y0, method=self.method,
+                        t_eval=t_eval, rtol=self.rtol, atol=self.atol)
+        if not sol.success:
+            raise RuntimeError(f"ODE solver failed: {sol.message}")
+
+        solutions = [sol.y[i] for i in range(sol.y.shape[0])]
+
+        print(f"[DEBUG] solutions length: {len(solutions)}")
+        for i, s in enumerate(solutions):
+            print(f"[DEBUG] solutions[{i}].shape: {s.shape}")
+        # data может быть списком массивов (для SoEq) или одним массивом (для Equation)
+        if isinstance(data, (list, tuple)):
+            loss = np.mean([np.sqrt(np.mean((solutions[i] - data[i]) ** 2)) for i in range(len(solutions))])
+        else:
+            loss = np.sqrt(np.mean((solutions[0] - data) ** 2))
+        return solutions, loss
+
+    def get_requirements(self):
+        return {
+            'form': 'explicit_ode',
+            'needs_initial_conditions': True,
+            'needs_boundary_conditions': False,
+            'supports_multisample': False
+        }

epde/solver/factory.py

@@ -0,0 +1,21 @@
+from epde.solver.base import BaseSolverAdapter
+from epde.solver.unified_deepxde_adapter import UnifiedDeepXDEAdapter
+from epde.solver.classical_ode_adapter import ClassicalODEAdapter
+
+class SolverFactory:
+    _registry = {
+        'deepxde': UnifiedDeepXDEAdapter,
+        'classical_ode': ClassicalODEAdapter,
+    }
+
+    @classmethod
+    def register(cls, name: str, adapter_class):
+        """Регистрирует новый тип солвера."""
+        cls._registry[name] = adapter_class
+
+    @classmethod
+    def create(cls, solver_type: str, **config) -> BaseSolverAdapter:
+        """Создаёт экземпляр адаптера по имени."""
+        if solver_type not in cls._registry:
+            raise ValueError(f"Неизвестный тип солвера: {solver_type}. Доступны: {list(cls._registry.keys())}")
+        return cls._registry[solver_type](**config)

epde/solver/ode_converter.py

@@ -0,0 +1,105 @@
+import sympy as sp
+import numpy as np
+from typing import Callable
+from epde.structure.main_structures import Equation
+
+class ODEToFirstOrder:
+    def __init__(self, t_symbol: str = 't'):
+        self.t = sp.Symbol(t_symbol, real=True)
+        self.debug = True  # можно выключить после отладки
+
+    def _build_expression(self, eq: Equation, var_name: str) -> sp.Expr:
+        u = sp.Function(var_name)(self.t)
+        expr = 0
+        if self.debug:
+            print(f"[DEBUG] Building expression for equation, target_idx={eq.target_idx}")
+        for term_idx, term in enumerate(eq.structure):
+            coeff = eq.weights_final[term_idx] if term_idx < len(eq.weights_final) else 1.0
+            term_expr = 1
+            if self.debug:
+                print(f"  Term {term_idx}: {term.name}, coeff={coeff}")
+            for factor in term.structure:
+                deriv_code = getattr(factor, "deriv_code", None)
+                is_deriv = factor.is_deriv and deriv_code is not None and len(deriv_code) > 0 and not all(v is None for v in deriv_code)
+                if self.debug:
+                    print(f"    Factor: is_deriv={factor.is_deriv}, variable={getattr(factor, 'variable', None)}, deriv_code={deriv_code}, params={factor.params}, is_deriv_flag={is_deriv}")
+                if is_deriv:
+                    # Настоящая производная (deriv_code содержит числа, не None)
+                    order = len(deriv_code)
+                    deriv = sp.Derivative(u, (self.t, order))
+                    power = factor.params[-1] if factor.params else 1.0
+                    term_expr *= deriv ** power
+                    if self.debug:
+                        print(f"      Derivative order {order}, power {power}, deriv={deriv}")
+                else:
+                    # Не производная – переменная или константа
+                    if hasattr(factor, 'variable') and factor.variable is not None:
+                        power = factor.params[-1] if factor.params else 1.0
+                        term_expr *= u ** power
+                        if self.debug:
+                            print(f"      Variable {factor.variable}, power {power}")
+                    else:
+                        if hasattr(factor, 'params') and len(factor.params) > 0:
+                            term_expr *= sp.Float(factor.params[-1])
+                            if self.debug:
+                                print(f"      Constant {factor.params[-1]}")
+                        else:
+                            if self.debug:
+                                print(f"      Unknown factor, ignored")
+            if term_idx == eq.target_idx:
+                expr -= coeff * term_expr
+                if self.debug:
+                    print(f"  Target term, adding -{coeff} * {term_expr}")
+            else:
+                expr += coeff * term_expr
+                if self.debug:
+                    print(f"  Adding +{coeff} * {term_expr}")
+        if self.debug:
+            print(f"  Total expression: {expr}")
+        return expr
+
+    def equation_to_rhs(self, eq: Equation, var_name: str) -> Callable:
+        expr = self._build_expression(eq, var_name)
+        u = sp.Function(var_name)(self.t)
+
+        derivs = [arg for arg in sp.preorder_traversal(expr) if isinstance(arg, sp.Derivative) and arg.args[0] == u]
+        if not derivs:
+            raise ValueError("Нет производных")
+        max_order = max(d.args[1][1] for d in derivs)
+        if self.debug:
+            print(f"[DEBUG] max_order = {max_order}")
+
+        u_deriv = sp.Derivative(u, (self.t, max_order))
+        sol = sp.solve(expr, u_deriv)
+        if self.debug:
+            print(f"[DEBUG] sol = {sol}")
+        if not sol:
+            raise ValueError(f"Не удалось выразить {u_deriv} из уравнения {expr}")
+        rhs_expr = sol[0]
+
+        y_sym = [sp.Symbol(f'{var_name}_{i}', real=True) for i in range(max_order)]
+        subs = {u: y_sym[0]}
+        for i in range(1, max_order):
+            subs[sp.Derivative(u, (self.t, i))] = y_sym[i]
+        rhs_expr = rhs_expr.subs(subs)
+        if self.debug:
+            print(f"[DEBUG] rhs_expr after substitution = {rhs_expr}")
+
+        rhs_expr = rhs_expr.replace(lambda x: isinstance(x, sp.Derivative), lambda x: 0)
+        rhs_expr = sp.simplify(rhs_expr)
+        if self.debug:
+            print(f"[DEBUG] rhs_expr after derivative replacement = {rhs_expr}")
+
+        derivatives = [y_sym[i+1] if i+1 < max_order else rhs_expr for i in range(max_order)]
+        if self.debug:
+            print(f"[DEBUG] derivatives = {derivatives}")
+
+        rhs_func_args = sp.lambdify([self.t] + y_sym, derivatives, modules='numpy')
+        def rhs_wrapper(t, y):
+            if self.debug:
+                print(f"[DEBUG] rhs_wrapper t={t}, y={y}")
+            res = rhs_func_args(t, *y)
+            if self.debug:
+                print(f"[DEBUG] rhs_wrapper res = {res}")
+            return np.array(res, dtype=float)
+        return rhs_wrapper

epde/solver/poka_tut_pobudem/__init__.py

@@ -0,0 +1,11 @@
+from epde.solver.base import BaseSolverAdapter
+from epde.solver.factory import SolverFactory
+from epde.solver.unified_deepxde_adapter import DeepXDEAdapter
+from epde.solver.classical_ode_adapter import ClassicalODEAdapter
+
+__all__ = [
+    'BaseSolverAdapter',
+    'SolverFactory',
+    'DeepXDEAdapter',
+    'ClassicalODEAdapter',
+]

epde/solver/poka_tut_pobudem/example.py

@@ -0,0 +1,95 @@
+import numpy as np
+import copy
+from scipy.integrate import solve_ivp
+from epde.interface.interface import EpdeSearch
+from epde.interface.equation_translator import translate_equation
+from epde.structure.main_structures import SoEq, Chromosome
+from epde.interface.token_family import TFPool
+from epde.operators.utils.default_parameter_loader import EvolutionaryParams
+from epde.solver.factory import SolverFactory
+
+def create_equation_from_str(eq_str, target_var, base_pool, all_vars):
+    families_copy = [copy.deepcopy(fam) for fam in base_pool.families]
+    for fam in families_copy:
+        if hasattr(fam, 'variable') and fam.variable is not None and fam.variable != target_var:
+            fam.status['demands_equation'] = False
+    temp_pool = TFPool(families_copy)
+    soeq = translate_equation(eq_str, temp_pool, all_vars=[target_var])
+    return soeq.vals[target_var]
+
+def lv_rhs(t, y):
+    u, v = y
+    alpha, beta, gamma, delta = 2/3, 4/3, 1.0, 1.0
+    du = alpha * u - beta * u * v
+    dv = delta * u * v - gamma * v
+    return [du, dv]
+
+t = np.linspace(0, 20, 200)
+sol_ref = solve_ivp(lv_rhs, (0,20), [1.0,1.0], t_eval=t, method='RK45', rtol=1e-6, atol=1e-9)
+exact_u, exact_v = sol_ref.y
+
+data_u = exact_u + 0.01 * np.random.normal(size=exact_u.shape)
+data_v = exact_v + 0.01 * np.random.normal(size=exact_v.shape)
+
+search = EpdeSearch(
+    use_solver=False,
+    multiobjective_mode=True,
+    coordinate_tensors=[t],
+    verbose_params={'show_iter_idx': False},
+    device='cpu'
+)
+search.set_preprocessor(default_preprocessor_type='FD', preprocessor_kwargs={})
+search.create_pool(data=[data_u, data_v], variable_names=['u', 'v'], max_deriv_order=1, additional_tokens=[])
+
+correct_eqs = [
+    '0.6666666666666666 * u{power: 1.0} + -1.3333333333333333 * u{power: 1.0} * v{power: 1.0} = du/dx0{power: 1.0}',
+    '1.0 * u{power: 1.0} * v{power: 1.0} + -1.0 * v{power: 1.0} = dv/dx0{power: 1.0}'
+]
+
+eq_u = create_equation_from_str(correct_eqs[0], 'u', search.pool, ['u', 'v'])
+eq_u.main_var_to_explain = 'u'
+eq_u.weights_internal = np.ones(len(eq_u.structure) - 1)
+eq_u.weights_internal_evald = True
+eq_u.weights_final_evald = True
+
+eq_v = create_equation_from_str(correct_eqs[1], 'v', search.pool, ['u', 'v'])
+eq_v.main_var_to_explain = 'v'
+eq_v.weights_internal = np.ones(len(eq_v.structure) - 1)
+eq_v.weights_internal_evald = True
+eq_v.weights_final_evald = True
+
+system = SoEq(search.pool, {})
+system.vals = Chromosome({'u': eq_u, 'v': eq_v}, {})
+system.moeadd_set = True
+
+def solve_with_solver(solver_type, solver_config, system, t, data):
+    adapter = SolverFactory.create(solver_type, **solver_config)
+    solutions, loss = adapter.solve(system, [t], data)
+    return solutions, loss
+
+
+print("=" * 50)
+print("Classical ODE solver (RK45)")
+print("=" * 50)
+
+solver_config_classical = {
+    "method": "RK45",
+    "rtol": 1e-6,
+    "atol": 1e-9,
+    "rhs": lv_rhs,
+    "y0": [1.0, 1.0]
+}
+solutions_cl, loss_cl = solve_with_solver("classical_ode", solver_config_classical, system, t, [data_u, data_v])
+print(f"Loss (RMSE): {loss_cl:.6f}")
+print(f"Max error u: {np.max(np.abs(solutions_cl[0] - exact_u)):.6f}")
+print(f"Max error v: {np.max(np.abs(solutions_cl[1] - exact_v)):.6f}")
+
+print("\n" + "=" * 50)
+print("DeepXDE (PINN) solver")
+print("=" * 50)
+
+solver_config_deepxde = EvolutionaryParams().get_default_params_for_operator('DeepXDEBasedFitness')
+solutions_dx, loss_dx = solve_with_solver("deepxde", solver_config_deepxde, system, t, [data_u, data_v])
+print(f"Loss (RMSE): {loss_dx:.6f}")
+print(f"Max error u: {np.max(np.abs(solutions_dx[0] - exact_u)):.6f}")
+print(f"Max error v: {np.max(np.abs(solutions_dx[1] - exact_v)):.6f}")