This PR adds parameter optimization utilities to feos-core and exposes Regressor objects to Python.

Everything works, but there is still some WIP (what methods and properties to expose, naming of methods/modules).

Todos

• Currently, the parameter optimization has its own error. We should combine it with FeosError. Either as variant in FeosError or we bring the variants over to FeosError.
• There is probably some boilerplate that can be solved with macros for the Datasets.

Open Decisions

• The whole optimization is currently behind the rayon feature gate for convenience. Is there a case/need for optimization without rayon?
• Should we expose more features (Parameter covalence matrix, Fisher information matrix, leverage, optimality) or leave this to users (unscaled Jacobian and residuals can be accessed).
• I changed the ParametersAD trait quite a bit. Tried to resolve the "parameter order and indexing" synchronisation problem. We can discuss whether that change is appropriate and if we want a different interface.
• Also in the ParametersAD trait: should we add the option to blacklist parameters for AD? E.g. association sites (although there is currently an example that takes the derivative w.r.t those). I added the option - not sure we want to keep it. Currently there is the option but it not applied.
• Remove adaptive penalty for non-converged data points? Introduces a lot of noise and I am not sure it is any good.
• Cleanup: I added docstrings that should be moved to proper theory docs.

AI usage

• Lots of documentation and tests were generated by Claude. I changed a lot manually but we should discuss whether we are OK with LLM generated and manually checked/patched code and docstrings.
• Basically the whole Python bindings were written by Claude (single shot in about 40 seconds). I thoroughly checked the code and changed things I found weird but it's 95% Claude. I'd like to discuss whether that's something we are OK with and how/if we should go about this in future PRs.

Usage Python

import feos 
import pandas as pd

df_vp = pd.read_csv("vapor_pressure.csv")
df_rho = pd.read_csv("liquid_density.csv")

# Read data from np.arrays.
vp = feos.VaporPressureDataset(
    df_vp.temperature_k.values, 
    df_vp.vapor_pressure_pa.values
)
rho = feos.LiquidDensityDataset(
    df_rho.temperature_k.values, 
    df_rho.pressure_pa.values, 
    df_rho.liquid_density_molm3
)

# or from CSV:
vp = feos.VaporPressureDataset.from_csv("vapor_pressure.csv")

# Define regressor: model, data, initial parameters and what to fit
regressor = feos.PureRegressor(
    model=feos.EquationOfStateAD.PcSaftNonAssoc,
    datasets=[vp, rho],
    params={"m": 2.0, "sigma": 3.0, "epsilon_k": 300.0, "mu": 0.0},
    fit=["m", "sigma", "epsilon_k"],
)

# Run fit with custom configuration
result = regressor.fit(
    config=feos.FitConfig(
        patience=100,
        ftol=1e-12,
        xtol=1e-12,
        gtol=1e-12,
        stepbound=100.0,
        strategy=feos.NonConvergenceStrategy.ignore(), # set residual and grad to zero
    )
)

# Print results
print(f"  {result}")
print(f"  all parameters: {result.all_parameters()}")

# Evaluate all datasets and create data frames with original and predicted data.
ds = regressor.evaluate_datasets(result.all_parameters())

# ds is a list of dicts readily convertable to DataFrames
dfs = [pd.DataFrame(d) for d in ds]
for df in dfs:
    # generate plots or do something else

Output (250 vapor pressures and 250 liquid densities for hexane):

  FitResult(converged=true, n_eval=8, elapsed=2.8ms, aad=[vapor pressure: 6.65%, liquid density: 1.02%])
  all parameters: {'m': 2.9455258525936583, 'sigma': 3.8549088883947196, 'epsilon_k': 242.1332276164755, 'mu': 0.0} 

Usage Rust

use feos::core::parameter_optimization::*;
use feos::pcsaft::PcSaftPure;
use std::collections::HashMap;
use std::path::Path;

fn main() {
    let vp = VaporPressureDataset::from_csv(Path::new("vapor_pressure.csv"))
        .expect("failed to read vapor pressure CSV");
    let rho = LiquidDensityDataset::from_csv(Path::new("liquid_density.csv"))
        .expect("failed to read liquid density CSV");

    let datasets = vec![
        PureDataset::VaporPressure(vp),
        PureDataset::LiquidDensity(rho),
    ];

    let params = HashMap::from([
        ("m".to_string(), 2.0),
        ("sigma".to_string(), 3.0),
        ("epsilon_k".to_string(), 300.0),
        ("mu".to_string(), 0.0),
    ]);
    let fit = ["m", "sigma", "epsilon_k"];

    // Runtime static version
    let regressor = PureRegressor::<PcSaftPure<f64, 4>>::new(datasets, params, &fit)
        .expect("failed to build regressor");

    let config = FitConfig {
        patience: 100,
        ftol: 1e-12,
        xtol: 1e-12,
        gtol: 1e-12,
        stepbound: 100.0,
        strategy: NonConvergenceStrategy::Ignore,
        ..Default::default()
    };
    
    // this is the way LM does it - we could return a FitResult directly.
    let (regressor, report) = regressor.fit(config);
    let result = regressor.to_result(&report);

    println!("converged:          {}", result.converged);
    println!("reason:             {}", result.termination_reason);
    println!("n_evaluations:      {}", result.n_evaluations);
    println!("aad_per_dataset:    {:?}", result.aad_per_dataset);
    println!("all_parameters:     {:?}", result.all_parameters());

    // evaluate data at optimum...
    let ds_results = regressor.evaluate_datasets(&result.optimal_params);

generates this output:

converged:          true
reason:             Converged { ftol: true, xtol: false }
n_evaluations:      8
aad_per_dataset:    [Some(6.645721326146425), Some(1.0224205499118824)]
all_parameters:     {"m": 2.9455258525936583, "sigma": 3.8549088883947196, "epsilon_k": 242.1332276164755, "mu": 0.0}