Merge branch 'DeepXDEBasedFitness' into main

Author: Yaroslav-Muravev

epde/operators/common/fitness.py

@@ -14,15 +14,10 @@
 from matplotlib import cm
 
 from epde.integrate import SolverAdapter
-# DeepXDEAdapter is imported lazily inside DeepXDEBasedFitness.apply() to
-# avoid triggering deepxde's import-time backend banner when no DeepXDE
-# solver is used (e.g. legacy L2/L2LR fitness paths).
+
 from epde.structure.main_structures import SoEq, Equation
 from epde.operators.utils.template import CompoundOperator
 import epde.globals as global_var
-from sklearn.linear_model import LinearRegression, Ridge
-from scipy.optimize import minimize
-from epde.supplementary import minmax_normalize
 from epde.supplementary import calculate_weights
 
 LOSS_NAN_VAL = 1e7
@@ -474,13 +469,18 @@ class DeepXDEBasedFitness(CompoundOperator):
 
     def __init__(self, param_keys: list):
         super().__init__(param_keys)
-        self.adapter = None
-
-    def set_adapter(self, config: dict = None, pretrained_net=None):
-        if self.adapter is None:
-            from epde.integrate.deepxde_integration import DeepXDEAdapter
-            cfg = self.params.get('deepxde_config', {}) if config is None else config
-            self.adapter = DeepXDEAdapter(pretrained_net=pretrained_net, **cfg)
+        self._solver = None
+
+    def _get_solver(self):
+        if self._solver is None:
+            from epde.solver.factory import SolverFactory
+            solver_type = self.params.get('solver_type', 'deepxde')
+            solver_config = self.params.get('solver_config', {})
+            # Обратная совместимость со старым параметром deepxde_config
+            if 'deepxde_config' in self.params and solver_type == 'deepxde':
+                solver_config = self.params['deepxde_config']
+            self._solver = SolverFactory.create(solver_type, **solver_config)
+        return self._solver
 
     def apply(self, objective, arguments: dict, force_out_of_place: bool = False):
         self_args, subop_args = self.parse_suboperator_args(arguments=arguments)
@@ -489,13 +489,7 @@ def apply(self, objective, arguments: dict, force_out_of_place: bool = False):
             self.suboperators['sparsity'].apply(objective, subop_args.get('sparsity', {}))
         self.suboperators['coeff_calc'].apply(objective, subop_args.get('coeff_calc', {}))
 
-        try:
-            pretrained_net = deepcopy(global_var.solution_guess_nn)
-        except:
-            pretrained_net = None
-        self.set_adapter(pretrained_net=pretrained_net)
-
-        keys, grids = global_var.grid_cache.get_all(mode='numpy')
+        keys, grids = global_var.grid_cache.get_all(mode='numpy', structural=True)
 
         if isinstance(objective, SoEq):
             data_list = []
@@ -507,15 +501,26 @@ def apply(self, objective, arguments: dict, force_out_of_place: bool = False):
             _, target, _ = objective.evaluate(normalize=False, return_val=False)
             data_list = [target.reshape(-1)]
 
+        solver = self._get_solver()
+
+        print(f"[DEBUG] objective type: {type(objective)}")
+        if isinstance(objective, SoEq):
+            print(f"[DEBUG] vars_to_describe: {objective.vars_to_describe}")
+            print(f"[DEBUG] data_list length: {len(data_list)}")
+            for i, d in enumerate(data_list):
+                print(f"[DEBUG] data_list[{i}].shape: {d.shape}")
+
         try:
-            solution_list, loss = self.adapter.solve(equation_or_system=objective,
-                                                     grids=grids,
-                                                     data=data_list)
+            solution_list, loss = solver.solve(equation_or_system=objective,
+                                               grids=grids,
+                                               data=data_list)
             if np.isnan(loss):
                 raise ValueError("NaN loss")
 
             if isinstance(objective, SoEq):
-                for idx, (var_name, eq) in enumerate({val: objective.vals[val] for val in objective.vars_to_describe}.items()):
+                # Перебираем уравнения в порядке vars_to_describe
+                for idx, var_name in enumerate(objective.vars_to_describe):
+                    eq = objective.vals[var_name]
                     err = self._compute_error(solution_list[idx], data_list[idx], eq)
                     if force_out_of_place:
                         pass
@@ -534,7 +539,9 @@ def apply(self, objective, arguments: dict, force_out_of_place: bool = False):
                     objective.fitness_calculated = True
                     self._compute_stability_for_equation(objective)
         except Exception as e:
-            print(f'[DeepXDEBasedFitness] DeepXDE solve failed: {e}')
+            print(f'[DeepXDEBasedFitness] Solver failed: {e}')
+            import traceback
+            traceback.print_exc()
             fitness_value = 1e7
             if force_out_of_place:
                 return fitness_value
@@ -566,7 +573,6 @@ def _compute_error(self, solution, data, eq):
         return err
 
     def _compute_stability_for_equation(self, eq: Equation):
-        # Повторно вычисляется evaluate
         _, target, features = eq.evaluate(normalize=False, return_val=False)
         data_shape = global_var.grid_cache.inner_shape
         self.get_g_fun_vals()

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