dfsc API

dfsc is the only public library name for the current differentiable fractional scientific-computing artifact.

The first mature component is a differentiable Mittag-Leffler spectral layer, exposed together with workflow helpers and component discovery.

Component Registry

import dfsc

print(dfsc.component_summary())
print(dfsc.implemented_components())
print(dfsc.algorithm_registry())

The registry separates implemented components from planned extension points, so future variable-order or distributed-order components are visible without being claimed as current capabilities.

Applicability Contract

import dfsc

report = dfsc.mlsl_applicability_report(layer.eigenvalues, layer.eigenvectors)
print(report.to_dict())

The report checks the current MLSL assumptions: retained diagonal spectral representation, finite real spectra, matching eigenvector dimensions, and the real non-positive Mittag-Leffler argument regime generated by non-negative times and non-negative Laplacian eigenvalues.

Core Layer

from dfsc import MittagLefflerSpectralLayer

layer = MittagLefflerSpectralLayer(eigenvalues, eigenvectors, terms=100)
u = layer(u0, times, alpha, beta=beta)

The layer evaluates

u(t) = Phi diag(E_alpha(-c^2 lambda_n^(beta/2) t^alpha)) Phi.T u0

where alpha and optional beta may require gradients.

Diagnostic Evaluation

evaluation = dfsc.evaluate_mittag_leffler(alpha, z, method="auto")
values = evaluation.values
print(evaluation.diagnostics())
print(evaluation.reliability.to_dict())

The diagnostic result reports the selected series/hybrid path, branch counts, finiteness, and disagreement with a richer embedded truncation. The disagreement is deliberately not presented as a rigorous error bound.

One restricted exception is available. For a real non-positive argument, 0 < alpha <= 1, and an alternating series whose omitted terms have entered their decreasing regime, alternating_series_remainder_bound returns the first omitted term as a certified series-truncation bound. Returning None means that the certificate is unavailable.

bound = dfsc.alternating_series_remainder_bound(alpha, z, terms=32)
budget = dfsc.ErrorBudget(rtol=1e-8, atol=1e-10)
report = dfsc.compose_error_budget_report(
    budget,
    reference_norm=float(values.norm()),
    evaluator_estimate=None if bound is None else float(bound.max()),
    evaluator_rigorous=bound is not None,
)

Unassessed Krylov or projection terms remain None; they are never silently counted as zero and therefore prevent a global-certification claim. Supplying a projection estimate does not make it rigorous unless the caller also sets projection_rigorous=True from an independently justified bound.

Set strict=True when an optimization or production workflow should reject inputs outside the documented reliability domain instead of accepting a merely finite tensor. Solutions returned by dfsc.solve expose the same conservative contract through solution.reliability and solution.summary().

Constructors

from dfsc import MLSLConfig, build_dirichlet_mlsl_1d, build_dirichlet_mlsl_2d

x, layer_1d = build_dirichlet_mlsl_1d(
    num_points=128,
    num_modes=24,
    config=MLSLConfig(terms=100),
)

coords, layer_2d = build_dirichlet_mlsl_2d(
    num_points_1d=32,
    num_modes_1d=8,
    config=MLSLConfig(terms=100),
)

Additional 1D boundary constructors are available:

from dfsc import (
    build_mixed_mlsl_1d,
    build_neumann_mlsl_1d,
    build_periodic_mlsl_1d,
)

2D tensor-product constructors are available for Dirichlet, Neumann, periodic, and mixed settings:

from dfsc import (
    build_dirichlet_mlsl_2d,
    build_mixed_mlsl_2d,
    build_neumann_mlsl_2d,
    build_periodic_mlsl_2d,
)

For broader negative-real spectral regimes:

x, layer = build_dirichlet_mlsl_1d(
    num_points=128,
    num_modes=24,
    config=MLSLConfig.stable(terms=120),
)

Operator and graph adapters are available when the user already has a finite dimensional self-adjoint positive semidefinite operator:

from dfsc import MLSLConfig, build_graph_mlsl, build_operator_mlsl

layer_from_matrix = build_operator_mlsl(operator_matrix, num_modes=64, config=MLSLConfig.stable())
layer_from_graph = build_graph_mlsl(adjacency_matrix, num_modes=64, config=MLSLConfig.stable())

These adapters extend dfsc to graph and unstructured-discretization workflows after discretization, but they do not solve mesh generation, finite-element assembly, or non-self-adjoint operator analysis.

Assembled finite-element-style stiffness and mass matrices are accepted through the generalized self-adjoint problem:

layer = dfsc.build_generalized_operator_mlsl(K, M, num_modes=64)
problem = dfsc.GeneralizedOperatorSpectralProblem(
    stiffness=K, mass=M, u0=u0, times=times, alpha=alpha
)
solution = dfsc.solve(problem)

The implementation solves K phi = lambda M phi and uses M phi for modal projection. Sparse assembly and mesh generation remain external responsibilities.

Problem--Algorithm--Solve Interface

problem = dfsc.FractionalSpectralProblem(
    layer=layer,
    u0=u0,
    times=times,
    alpha=alpha,
    beta=beta,
)
solution = dfsc.solve(problem)

The solve interface provides a compact library-level contract: problem objects hold data and parameters, algorithm objects declare the chosen dfsc method, and Solution stores values, times, metadata, and diagnostics. This follows the problem/algorithm separation common in mature scientific-computing ecosystems while remaining scoped to Mittag-Leffler spectral dynamics.

AutoDFSC estimates the largest real spectral argument using detached problem metadata and selects the direct or stable evaluator. Explicit MLSLDirect and MLSLStable objects now reconfigure the prebuilt layer rather than only naming the returned solution.

For OperatorSpectralProblem, AutoDFSC selects MLSLKrylov when the matrix dimension exceeds dense_eigh_limit and no explicit retained mode count was requested:

policy = dfsc.AutoDFSC(dense_eigh_limit=256, krylov_dimension=48)
solution = dfsc.solve(operator_problem, policy)

MLSLKrylov uses fully reorthogonalized Lanczos projection and diagonalizes only the reduced tridiagonal matrix. It supports leading batches of initial states and preserves gradients with respect to the initial state and fractional orders along the realized Lanczos path. Its embedded m versus m-8 disagreement is a convergence diagnostic, not a certified error bound. The the tensor path accepts dense or PyTorch sparse symmetric positive-semidefinite matrices.

Matrix-free operators use an explicit mathematical contract:

operator = dfsc.SelfAdjointLinearOperator(
    size=n,
    matvec=matvec,
    dtype=torch.float64,
    device="cuda",
    symmetric=True,
    positive_semidefinite=True,
)
problem = dfsc.LinearOperatorSpectralProblem(operator, u0, times, alpha)
solution = dfsc.solve(problem)

The callable must preserve vector shape, dtype, and device. Its self-adjoint and positive-semidefinite flags are user assertions: a finite Lanczos run cannot prove these global properties. dfsc checks sparse tensor symmetry and rejects significantly negative reduced Ritz values, but those checks are not a formal PSD certificate. Gradients can pass through differentiable tensors captured by the matvec callable.

For repeated time/order queries with a fixed operator and initial-state batch, prepare the Lanczos spaces once:

prepared = dfsc.prepare_lanczos_basis(A, u0_batch, krylov_dimension=32)
prediction = dfsc.apply_prepared_lanczos_basis(prepared, times, alpha, beta=beta)

The query retains gradients with respect to alpha and beta. Rebuild the basis when the operator or initial state changes.

Local inverse-problem diagnostics are available separately:

report = dfsc.local_identifiability(loss_fn, fitted_parameters, noise_variance=sigma2)

This Hessian-based report is local and model-conditional. It does not prove global identifiability.

History-Aware Fallback

problem = dfsc.CaputoL1Problem(
    operator=A,
    u0=u0,
    alpha=alpha,
    final_time=1.0,
    num_steps=200,
    forcing=forcing,  # optional (..., num_steps + 1, state_size)
)
solution = dfsc.solve(problem)

This first fallback solves constant-order linear Caputo systems with 0 < alpha < 1 using a differentiable implicit L1 scheme on a uniform grid. It stores the full history and is not a fast-memory or adaptive solver.

FFT History Operator

For a trajectory already sampled on a uniform grid, dfsc can evaluate all Caputo-L1 derivatives by direct or FFT convolution:

derivatives, diagnostics = dfsc.caputo_l1_history(
    trajectory,
    alpha=alpha,
    final_time=1.0,
    method="auto",
)

The default time axis is 0 for a scalar trajectory (T+1,) and -2 for batched state trajectories (..., T+1, state_size). method="auto" uses the direct reference implementation for short histories and FFT convolution for long histories. Both paths preserve gradients with respect to the trajectory and scalar alpha for 0 < alpha < 1.

The same capability is available through the solve interface:

problem = dfsc.CaputoHistoryProblem(trajectory, alpha, final_time=1.0)
solution = dfsc.solve(problem)

This is an offline, full-trajectory residual operator with quasi-linear FFT work. It still stores the supplied trajectory and does not replace the sequential implicit CaputoL1Problem solver. It is not an adaptive CQ method, an SOE compressor, or an online fast-memory time stepper.

Complex and General Operators

The controlled complex scalar evaluator accepts moderate complex arguments:

evaluation = dfsc.evaluate_complex_mittag_leffler(alpha, complex_z)
print(evaluation.diagnostics())

Its default validated region is |z| <= 4; embedded truncation disagreement is reported but is not a rigorous error bound. Inputs outside this radius are rejected unless the caller explicitly opts into an unvalidated raw series.

General real, complex, non-self-adjoint, or non-normal operators use Arnoldi:

problem = dfsc.GeneralOperatorProblem(A, u0, times, alpha)
solution = dfsc.solve(problem, dfsc.MLSLArnoldi())

The reduced Hessenberg matrix function is evaluated directly by a matrix power series, so this path does not require diagonalizability. Diagnostics report the effective Arnoldi dimension, breakdowns, reduced radius, an empirical non-normality measure, and disagreement with a smaller subspace. The default path rejects reduced radii above 4; large complex sectors, highly non-normal transient regimes, and pseudospectral error certification remain unsupported.

Current Numerical Scope

  • method="series" is the default and supports explicit custom backward.
  • MLSLConfig.stable() uses method="hybrid" for non-positive real arguments in longer-time or stiffer regimes. The backward pass is differentiable through PyTorch autograd over the active series/asymptotic branches.
  • The built-in bases cover 1D/2D Dirichlet, Neumann, periodic, and mixed cases, with dense operator, graph, and generalized stiffness/mass adapters.
  • The primitive is history-free at evaluation time: it maps u0, query times, and fractional parameters directly to states without storing all past steps.

Expected Use

Use MLSL as a physics-structured component inside inverse problems, operator learning pipelines, hybrid neural models, and differentiable simulators where fractional dynamics are present or hypothesized.

Workflow Helpers

Bounded trainable fractional orders:

from dfsc import make_trainable_orders

orders = make_trainable_orders(alpha_init=0.9, beta_init=1.4)
loss = model_loss(orders.alpha, orders.beta)
loss.backward()

Hybrid residual composition:

from dfsc import HybridResidualModel

model = HybridResidualModel(backbone=mlsl_layer, residual_head=neural_head)
u_pred = model(u0, times, alpha, beta=beta)

Forced Dynamics

The first forced extension is available as:

from dfsc import ForcedMittagLefflerSpectralLayer

forced_layer = ForcedMittagLefflerSpectralLayer(base_layer, forcing_terms=100)
u = forced_layer(u0, times, alpha, forcing_values, forcing_times, beta=beta)

It uses a midpoint quadrature approximation of the Duhamel-type convolution kernel involving E_{alpha,alpha}.

The current implementation uses a hybrid two-parameter Mittag-Leffler evaluator for non-positive real spectral arguments and has manufactured-solution validation in a multi-mode spectral space.

For float32 or lower precision, the hybrid evaluator automatically uses a more conservative series/asymptotic switching threshold to avoid high-order series overflow on GPU.

Semilinear Mild Solutions

problem = dfsc.SemilinearSpectralProblem(
    layer=layer,
    u0=u0,
    times=times,
    alpha=alpha,
    beta=beta,
    nonlinearity=lambda state: -gamma * state.pow(3),
)
solution = dfsc.solve(problem, dfsc.MLSLPicard())

The solver evaluates the Duhamel term by midpoint quadrature and iterates the mild equation to a relative fixed-point tolerance. Solution.retcode is "maxiters" when the tolerance is not reached. This implementation covers one unbatched state on increasing query times beginning at zero; it does not claim global convergence for arbitrary nonlinearities.

Experimental Order Wrappers

from dfsc import DistributedOrderMLSL, VariableOrderMLSL

vo_layer = VariableOrderMLSL(base_layer, order_fn)
do_layer = DistributedOrderMLSL(base_layer, alpha_nodes, trainable_weights=True)

VariableOrderMLSL evaluates the backbone independently with one alpha per query time. DistributedOrderMLSL forms a differentiable weighted mixture over retained alpha nodes. Both wrappers are intended for early ablations and interface validation; they are not advertised as general variable-order or distributed-order fractional solvers.

Public Benchmark Data Contract

from dfsc import benchmark_targets, validate_dataset_manifest

print(benchmark_targets())
ok, missing = validate_dataset_manifest(manifest)

dfsc does not bundle or invent real physical datasets. The benchmark contract requires dataset name, domain, task, source, license, citation, splits, and tensor-file mapping before a dataset can be treated as a reproducible benchmark.

Domain Application Templates

The domain layer packages an existing dfsc problem with a recommended algorithm, assumptions, limitations, and differentiable-parameter list:

case = dfsc.advection_diffusion_case(
    initial=lambda x: torch.sin(2 * torch.pi * x),
    times=torch.tensor([0.0, 0.01]),
    alpha=torch.tensor(0.9, requires_grad=True),
    diffusivity=torch.tensor(0.02, requires_grad=True),
    velocity=torch.tensor(0.15, requires_grad=True),
    num_points=16,
    arnoldi_dimension=16,
)
solution = case.solve()

Available constructors are anomalous_diffusion_case, assembled_relaxation_case, network_diffusion_case, and advection_diffusion_case. These are tested translations into existing dfsc algorithms, not independent domain solvers. application_catalog() reports the validated assumptions and remaining limitations.