Optimize FEBio simulations
==========================

Purpose
-------

The ``pyfebiopt.optimize`` module builds on XPLT reading. Optimization adapters
usually open the ``.xplt`` generated by a FEBio run, slice the variables needed
for one objective, and return arrays to the residual assembler.

Core classes
------------

``ExperimentSeries``
   Stores measured ``x`` and ``y`` arrays, plus optional weights. These are the
   reference data used to form residuals.

``Parameter`` and ``ParameterSpace``
   Define fitted variables, initial values, bounds, and active/inactive flags.
   ``ParameterSpace`` maps physical parameters :math:`\theta` to optimizer-space
   variables :math:`\phi` when reparameterization is enabled.

``FebTemplate``, ``ParameterBinding``, and ``EvaluationBinding``
   Connect Python parameters to XML nodes in a FEBio ``.feb`` file.
   ``ParameterBinding`` writes one parameter value directly. ``EvaluationBinding``
   computes an arbitrary value from the full parameter dictionary, useful for
   coupled values such as mirrored fiber directions.

``SimulationAdapter``
   Wraps a callable that reads one FEBio result file and returns simulated ``x``
   and ``y`` arrays. Adapters are where ``.xplt`` slicing and project-specific
   post-processing normally live.

``SimulationCase``
   Groups one template, one set of experiments, matching adapters, grid alignment
   policies, output folder names, and per-case FEBio runtime options.

``EngineOptions``
   Collects configuration for runners, storage, cleanup, Jacobians, logging,
   monitoring, and optimizer selection.

``Engine``
   Coordinates the full loop: render FEB files, launch FEBio, collect outputs,
   assemble residuals, log/report progress, and call the selected SciPy optimizer.

Assemble a script
-----------------

Minimal scripts follow the same shape even when the objective is complex:

.. code-block:: python

   import numpy as np

   from pyfebiopt.optimize import (
       CleanupOptions,
       Engine,
       EngineOptions,
       ExperimentSeries,
       FebTemplate,
       GridPolicyOptions,
       JacobianOptions,
       MonitorOptions,
       OptimizerOptions,
       ParameterBinding,
       ParameterSpace,
       RunnerOptions,
       SimulationAdapter,
       SimulationCase,
       StorageOptions,
   )
   from pyfebiopt.xplt import xplt

   def read_force_curve(xplt_path):
       xp = xplt(str(xplt_path))
       xp.readAllStates()
       force = xp.results["reaction forces"].nodeset("grip").comp("z")
       time = np.asarray(xp._time, dtype=float)
       return time, np.asarray(force, dtype=float)

   exp_data = np.loadtxt("data.txt")
   experiment = ExperimentSeries(x=exp_data[:, 0], y=exp_data[:, 1])

   parameters = ParameterSpace(xi=2.0)
   parameters.add_parameter(name="G", theta0=0.6, bounds=(0.0, 10.0))

   template = FebTemplate(
       "model_template.feb",
       bindings=[
           ParameterBinding(theta_name="G", xpath=".//Material/material[@id='1']/G"),
       ],
   )

   case = SimulationCase(
       template=template,
       subfolder="biaxial",
       experiments={"force": experiment},
       adapters={"force": SimulationAdapter(read_force_curve)},
       grids={"force": GridPolicyOptions(policy="exp_to_sim")},
       omp_threads=4,
   )

   options = EngineOptions(
       jacobian=JacobianOptions(enabled=True, perturbation=1e-3, parallel=True),
       cleanup=CleanupOptions(remove_previous=True, mode="all"),
       runner=RunnerOptions(jobs=2, command=("febio4", "-i")),
       storage=StorageOptions(mode="disk", root="runs"),
       monitor=MonitorOptions(enabled=True),
       optimizer=OptimizerOptions(
           name="least_squares",
           settings={"ftol": 1e-6, "xtol": 1e-6},
           reparametrize=True,
       ),
   )

   result = Engine(parameters, [case], options=options).run()
   print(result.theta)

Objective
---------

The default least-squares objective is assembled from all cases and experiments:

.. math::

   J(\theta) = \sum_{c \in \mathrm{cases}} \sum_{e \in \mathrm{experiments}}
               w_{c,e} \lVert f_{c,e}(\theta) - y_{c,e} \rVert_2^2

where :math:`f_{c,e}` is produced by the adapter, :math:`y_{c,e}` is the
experimental series, and :math:`w_{c,e}` is an optional weight.