diff --git a/dev/vfm/call_vfm.py b/dev/vfm/call_vfm.py
index 8879176da..339a42849 100644
--- a/dev/vfm/call_vfm.py
+++ b/dev/vfm/call_vfm.py
@@ -13,7 +13,11 @@
     SpecimenGeometry,
 )
 from pyvale.vfm.hardening import LinearHardening
-from pyvale.vfm.identification import Identification, IdentificationPhase
+from pyvale.vfm.identification import run_identification
+from pyvale.vfm.identificationconfig import (
+    IdentificationConfig,
+    IdentificationPhase,
+)
 from pyvale.vfm.metricsbvf import SensitivityBasedVirtualFieldsMetric
 from pyvale.vfm.objectivefuncvector import VectorFirstResultPassthrough
 from pyvale.vfm.optimiserleastsquares import LeastSquares
@@ -21,7 +25,6 @@
     HomogeneousSpatialParameterisation,
 )
 from pyvale.vfm.spatialparamknown import KnownSpatialParameterisation
-from pyvale.vfm.vfm import run_identification
 
 inputs_path = Path(__file__).resolve().parent / "inputs"
 
@@ -100,7 +103,7 @@ def main():
         )
     ]
 
-    identification = Identification(
+    identification_config = IdentificationConfig(
         IsotropicVonMisesElastoplasticity(
             LinearHardening()
         ),
@@ -108,7 +111,7 @@ def main():
         phases
     )
 
-    vfm_result = run_identification(experiment_data, identification)
+    vfm_result = run_identification(experiment_data, identification_config)
     print(vfm_result)
 
 
diff --git a/src/pyvale/vfm/__init__.py b/src/pyvale/vfm/__init__.py
index 83356d6a5..771c2efaf 100644
--- a/src/pyvale/vfm/__init__.py
+++ b/src/pyvale/vfm/__init__.py
@@ -4,10 +4,10 @@
 # Copyright (C) 2025 The Computer Aided Validation Team
 #===============================================================================
 
-from .vfm import *
+from .identification import *
+from .identificationconfig import *
 
 from .experimentdata import *
-from .identification import *
 
 from .constlaw import *
 from .constlaws import *
diff --git a/src/pyvale/vfm/constlaw.py b/src/pyvale/vfm/constlaw.py
index 948d94f86..6483595a8 100644
--- a/src/pyvale/vfm/constlaw.py
+++ b/src/pyvale/vfm/constlaw.py
@@ -6,13 +6,51 @@
 
 
 class EIdentificationType(enum.Enum):
+    """
+    Identifies whether a constitutive law's parameters should be
+    identified with linear on nonlinear identification
+    """
+
     Linear = enum.auto()
+    """Identify parameters with linear identification"""
     Nonlinear = enum.auto()
+    """Identify parameters with nonlinear identification"""
 
 
 class IConstitutiveLaw(ABC):
+    """
+    Interface (abstract base class) for a constitutive law.
+
+    Provides the material model that relates strain to stress. Concrete
+    implementations define the specific constitutive equations and report
+    whether identification is linear or nonlinear
+    """
+
     @abstractmethod
     def get_identification_type(self) -> EIdentificationType:
+        """
+        Indicate whether this law is linear or nonlinear in its parameters.
+
+        Returns
+        -------
+        EIdentificationType
+            ``Linear`` or ``Nonlinear``
+        """
+        pass
+
+    @abstractmethod
+    def get_required_parameters(self) -> list[str]:
+        """
+        Return the list of required constitutive parameters for this law.
+
+        Concrete implementations should combine their own parameter names
+        with those from any nested hardening law
+
+        Returns
+        -------
+        list[str]
+            All parameter name strings this law requires
+        """
         pass
 
     @abstractmethod
@@ -21,4 +59,20 @@ def calculate_stress(
         strain: npt.NDArray[np.float64],
         constitutive_parameter_maps: dict[str, npt.NDArray[np.float64]],
     ) -> npt.NDArray[np.float64]:
+        """
+        Compute stress from the current strain and parameter maps.
+
+        Parameters
+        ----------
+        strain : npt.NDArray[np.float64]
+            Full-field strain history, shape ``(timesteps, components, y, x)``
+        constitutive_parameter_maps : dict[str, npt.NDArray[np.float64]]
+            Dictionary of current parameter values as 2D maps,
+            keyed by parameter name
+
+        Returns
+        -------
+        npt.NDArray[np.float64]
+            Stress field with the same shape as ``strain``
+        """
         pass
diff --git a/src/pyvale/vfm/constlaws.py b/src/pyvale/vfm/constlaws.py
index 79c005166..094cfd28e 100644
--- a/src/pyvale/vfm/constlaws.py
+++ b/src/pyvale/vfm/constlaws.py
@@ -38,6 +38,11 @@ def __init__(
     def get_identification_type(self) -> EIdentificationType:
         return EIdentificationType.Nonlinear
 
+    def get_required_parameters(self) -> list[str]:
+        params = [self.elastic_modulus_label, self.poissons_ratio_label]
+        params.extend(self.hardening_function.get_required_parameters())
+        return params
+
     def calculate_stress(
         self,
         strain: npt.NDArray[np.float64],
diff --git a/src/pyvale/vfm/constparam.py b/src/pyvale/vfm/constparam.py
index c5f29f437..8d5f4ddf5 100644
--- a/src/pyvale/vfm/constparam.py
+++ b/src/pyvale/vfm/constparam.py
@@ -6,19 +6,58 @@
 
 @dataclass(slots=True)
 class ConstitutiveParameter:
-    value: npt.NDArray[np.float64]
+    """
+    A single constitutive parameter with a spatially varying map of values and
+    associated value bounds.
+
+    The ``value`` may be supplied as a scalar (int or float) together with a
+    ``map_size`` to create a homogeneous 2D field, or as a full 2D array
+    directly
+    """
+
+    map: npt.NDArray[np.float64]
+    """Parameter map, shape ``(y, x)``"""
+
     lower_bound: float
+    """Lower bound for the parameter value"""
+
     upper_bound: float
+    """Upper bound for the parameter value"""
 
-    # TODO: enforce bounds on input value
-    # map_size must be in form (y, x)
     def __init__(
         self,
         value: int | float | npt.NDArray[np.float64],
         lower_bound: float,
         upper_bound: float,
-        map_size: npt.NDArray[np.float64] | None = None
+        map_size: npt.NDArray[np.float64] | None = None,
     ) -> None:
+        """
+        Parameters
+        ----------
+        value : int | float | npt.NDArray[np.float64]
+            Parameter value.  If scalar, ``map_size`` must also be provided
+            and the value is broadcast to create a homogeneous 2D field
+        lower_bound : float
+            Lower bound for the parameter
+        upper_bound : float
+            Upper bound for the parameter
+        map_size : npt.NDArray[np.float64] | None, optional
+            Shape ``(y, x)`` of the spatial parameterisation when ``value``
+            is a scalar.  Ignored when ``value`` is already an array
+
+        Raises
+        ------
+        ValueError
+            If ``lower_bound >= upper_bound``, or if any value in the
+            parameter map lies outside ``[lower_bound, upper_bound]``, or
+            if ``value`` is a scalar and ``map_size`` is ``None``
+        """
+        if lower_bound >= upper_bound:
+            raise ValueError(
+                f"lower_bound ({lower_bound}) must be less than "
+                f"upper_bound ({upper_bound})"
+            )
+
         if isinstance(value, (int, float)):
             if map_size is None:
                 raise ValueError(
@@ -26,11 +65,16 @@ def __init__(
                     "parameter value is int or float"
                 )
 
-            self.value = np.full((map_size[0], map_size[1]), value)
+            self.map = np.full((map_size[0], map_size[1]), value)
 
         else:
-            self.value = value
+            self.map = value
+
+        if np.any((self.map < lower_bound) | (self.map > upper_bound)):
+            raise ValueError(
+                f"parameter values must be within provided bounds "
+                f"[{lower_bound}, {upper_bound}]"
+            )
 
         self.lower_bound = lower_bound
         self.upper_bound = upper_bound
-
diff --git a/src/pyvale/vfm/dof.py b/src/pyvale/vfm/dof.py
index 1728bdd88..293fb808e 100644
--- a/src/pyvale/vfm/dof.py
+++ b/src/pyvale/vfm/dof.py
@@ -3,6 +3,14 @@
 
 @dataclass(slots=True)
 class DegreeOfFreedom:
+    """
+    A single scalar degree of freedom with bounds.
+
+    Used by the optimiser to explore the design space. Values are typically
+    normalised to ``[0, 1]`` during optimisation and denormalised to the
+    physical range ``[lower_bound, upper_bound]`` for evaluation
+    """
+
     value: float
     lower_bound: float
     upper_bound: float
diff --git a/src/pyvale/vfm/experimentdata.py b/src/pyvale/vfm/experimentdata.py
index 5d6f86b2f..5a146c15d 100644
--- a/src/pyvale/vfm/experimentdata.py
+++ b/src/pyvale/vfm/experimentdata.py
@@ -7,38 +7,95 @@
 
 @dataclass(slots=True)
 class SpecimenGeometry:
+    """
+    Physical geometry of the test specimen.
+
+    Stores the spatial coordinates, region-of-interest mask, thickness, and
+    per-point physical area of the DIC grid
+    """
+
     x: npt.NDArray[np.float64]
+    """
+    x-coordinates at each grid point, shape ``(y, x)`` (mm).
+
+    Always positive, increasing left to right (column index)
+    """
+
     y: npt.NDArray[np.float64]
-    # TODO: not sure exactly what this should be
-    region_of_interest: npt.NDArray[np.uint32]
+    """
+    y-coordinates at each grid point, shape ``(y, x)`` (mm).
+
+    Always positive, increasing top to bottom (row index)
+    """
+
+    region_of_interest: npt.NDArray[np.bool_]
+    """Boolean mask of valid analysis points, shape ``(y, x)``"""
+
     thickness: float
+    """Out-of-plane thickness of the specimen (mm)"""
+
     pixel_area: npt.NDArray[np.float64]
+    """Area per grid point, shape ``(y, x)`` (mm²)"""
 
 
 class EEdgeCondition(enum.Enum):
+    """Mechanical condition applied to an edge of the specimen"""
+
     Free = enum.auto()
+    """Unconstrained edge (stress-free)"""
+
     Fixed = enum.auto()
+    """Fully constrained edge (zero displacement)"""
+
     Traction = enum.auto()
+    """Edge with a known applied traction (force) applied"""
 
 
 @dataclass(slots=True)
 class Edge:
+    """Boundary condition for the two orthogonal directions on a single edge"""
+
     x: EEdgeCondition
+    """Condition in the global x-direction"""
+
     y: EEdgeCondition
+    """Condition in the global y-direction"""
 
 
 @dataclass(slots=True)
 class EdgeConditions:
+    """
+    Boundary conditions on the four edges of the specimen.
+
+    Edges are identified by the minimum/maximum coordinate value along each
+    axis
+    """
+
     min_x_edge: Edge
+    """Condition along the minimum x edge"""
+
     max_x_edge: Edge
+    """Condition along the maximum x edge"""
+
     min_y_edge: Edge
+    """Condition along the minimum y edge"""
+
     max_y_edge: Edge
+    """Condition along the maximum y edge"""
 
 
 @dataclass(slots=True)
 class BoundaryConditions:
+    """Combined kinematic and kinetic boundary conditions for the experiment"""
+
     edge_conditions: EdgeConditions
+    """Kinematic constraints on all four edges of the specimen"""
+
     force: npt.NDArray[np.float64]
+    """
+    Measured force history, shape ``(timesteps, 2)`` with columns
+    ``[Fx, Fy]`` (x-direction, y-direction)
+    """
 
 
 def _calculate_timestep_deltas(
@@ -54,11 +111,57 @@ def _calculate_timestep_deltas(
 
 @dataclass(slots=True)
 class ExperimentData:
+    """
+    Input data from a DIC experiment.
+
+    Stores the full-field strain history, specimen geometry, boundary
+    conditions, and temporal data needed to perform VFM identification.
+
+    Shape conventions
+    -----------------
+    strain                 (timesteps, components, y, x)
+    specimen_geometry:
+        x                  (y, x)
+        y                  (y, x)
+        pixel_area         (y, x)
+        region_of_interest (y, x)
+    boundary_conditions:
+        force              (timesteps, 2)  ``[Fx, Fy]``
+    timesteps              (timesteps,)
+
+    Coordinate system
+    -----------------
+    ``x`` increases left to right (column index)
+    ``y`` increases top to bottom (row index)
+    All coordinates are always positive, and start at 0.0
+
+    Notes
+    -----
+    ``delta_timesteps`` is computed automatically from ``timesteps`` on init
+    and is not user-supplied
+    """
+
     strain: npt.NDArray[np.float64]
+    """
+    Full-field strain history, shape ``(timesteps, components, y, x)``
+    where ``x`` increases left to right and ``y`` increases top to bottom.
+    Components are ordered as ``[xx, yy, xy]`` (normal x, normal y, shear xy)
+    """
+
     specimen_geometry: SpecimenGeometry
+    """Geometry of the specimen"""
+
     boundary_conditions: BoundaryConditions
+    """Kinematic and kinetic boundary conditions applied during the test"""
+
     timesteps: npt.NDArray[np.float64]
+    """Time value at each frame / load step, shape ``(timesteps,)``"""
+
     delta_timesteps: npt.NDArray[np.float64]
+    """
+    Time increment between consecutive frames (computed automatically),
+    shape ``(timesteps,)``
+    """
 
     def __init__(
         self,
diff --git a/src/pyvale/vfm/hardening.py b/src/pyvale/vfm/hardening.py
index 8a5d96322..addeb01f7 100644
--- a/src/pyvale/vfm/hardening.py
+++ b/src/pyvale/vfm/hardening.py
@@ -26,6 +26,19 @@ def hardening(
         """
         pass
 
+    @abstractmethod
+    def get_required_parameters(self) -> list[str]:
+        """
+        Return the list of constitutive parameters required for
+        this hardening law.
+
+        Returns
+        -------
+        list[str]
+            All parameter name strings this hardening law requires
+        """
+        pass
+
 
 @dataclass(slots=True)
 class LinearHardening(IHardeningFunction):
@@ -47,6 +60,9 @@ def __init__(
         else:
             self.hardening_modulus_label = "hardening_modulus"
 
+    def get_required_parameters(self) -> list[str]:
+        return [self.yield_strength_label, self.hardening_modulus_label]
+
     def hardening(
         self,
         constitutive_parameter_maps: dict[str, npt.NDArray[np.float64]],
@@ -94,6 +110,13 @@ def __init__(
         else:
             self.hardening_exponent_label = "hardening_exponent"
 
+    def get_required_parameters(self) -> list[str]:
+        return [
+            self.strength_coefficient_label,
+            self.strain_offset_label,
+            self.hardening_exponent_label,
+        ]
+
     def hardening(
         self,
         constitutive_parameter_maps: dict[str, npt.NDArray[np.float64]],
@@ -160,6 +183,14 @@ def __init__(
         else:
             self.rate_parameter_label = "rate_parameter"
 
+    def get_required_parameters(self) -> list[str]:
+        return [
+            self.yield_strength_label,
+            self.hardening_modulus_label,
+            self.saturation_stress_label,
+            self.rate_parameter_label,
+        ]
+
     def hardening(
         self,
         constitutive_parameter_maps: dict[str, npt.NDArray[np.float64]],
@@ -224,6 +255,13 @@ def __init__(
         else:
             self.hardening_exponent_label = "hardening_exponent"
 
+    def get_required_parameters(self) -> list[str]:
+        return [
+            self.yield_strength_label,
+            self.strength_coefficient_label,
+            self.hardening_exponent_label,
+        ]
+
     def hardening(
         self,
         constitutive_parameter_maps: dict[str, npt.NDArray[np.float64]],
diff --git a/src/pyvale/vfm/identification.py b/src/pyvale/vfm/identification.py
index 0377187b8..4305563f1 100644
--- a/src/pyvale/vfm/identification.py
+++ b/src/pyvale/vfm/identification.py
@@ -1,25 +1,48 @@
-from dataclasses import dataclass
+import numpy as np
 
-from pyvale.vfm.constlaw import IConstitutiveLaw
+from pyvale.vfm.constlaw import EIdentificationType
 from pyvale.vfm.constparam import ConstitutiveParameter
-from pyvale.vfm.metric import IMetric
-from pyvale.vfm.objectivefunc import IObjectiveFunction
-from pyvale.vfm.optimiser import IOptimiser
-from pyvale.vfm.spatialparam import ISpatialParameterisation
+from pyvale.vfm.experimentdata import ExperimentData
+from pyvale.vfm.validation import (
+    validate_experiment_data,
+    validate_identification_config,
+)
+from pyvale.vfm.identificationconfig import IdentificationConfig
 
 
-@dataclass(slots=True)
-class IdentificationPhase:
-    # Param name str must be same as in the input params
-    spatial_parameterisations: dict[str, ISpatialParameterisation]
-    metrics: list[IMetric]
-    objective_function: IObjectiveFunction
-    optimiser: IOptimiser
+def run_identification(
+    experiment_data: ExperimentData,
+    identification_config: IdentificationConfig
+) -> dict[str, ConstitutiveParameter]:
+    validate_experiment_data(experiment_data)
+    validate_identification_config(identification_config)
 
+    match identification_config.constitutive_law.get_identification_type():
+        # TODO: implement linear case
+        case EIdentificationType.Linear:
+            ...
+        case EIdentificationType.Nonlinear:
+            parameter_map_size = np.array(
+                experiment_data.specimen_geometry.x.shape,
+                dtype=np.uint32
+            )
 
-@dataclass(slots=True)
-class Identification():
-    constitutive_law: IConstitutiveLaw
-    parameters: dict[str, ConstitutiveParameter]
-    phases: list[IdentificationPhase]
+            for phase in identification_config.phases:
+                for param_name, sp in phase.spatial_parameterisations.items():
+                    sp.initialise_from_constitutive_parameter(identification_config.parameters[param_name])
 
+                optimised_spatial_parameterisations = phase.optimiser.optimise(
+                    identification_config.constitutive_law,
+                    parameter_map_size,
+                    phase.spatial_parameterisations,
+                    phase.metrics,
+                    phase.objective_function,
+                    experiment_data
+                )
+
+                for param_name, sp in optimised_spatial_parameterisations.items():
+                    identification_config.parameters[param_name].map = (
+                        sp.to_map(parameter_map_size)
+                    )
+
+    return identification_config.parameters
diff --git a/src/pyvale/vfm/identificationconfig.py b/src/pyvale/vfm/identificationconfig.py
new file mode 100644
index 000000000..2758069fa
--- /dev/null
+++ b/src/pyvale/vfm/identificationconfig.py
@@ -0,0 +1,53 @@
+from dataclasses import dataclass
+
+from pyvale.vfm.constlaw import IConstitutiveLaw
+from pyvale.vfm.constparam import ConstitutiveParameter
+from pyvale.vfm.metric import IMetric
+from pyvale.vfm.objectivefunc import IObjectiveFunction
+from pyvale.vfm.optimiser import IOptimiser
+from pyvale.vfm.spatialparam import ISpatialParameterisation
+
+
+@dataclass(slots=True)
+class IdentificationPhase:
+    """
+    A single identification phase.
+
+    Each phase defines its own spatial parameterisations, metrics, objective
+    function, and optimiser
+    """
+
+    spatial_parameterisations: dict[str, ISpatialParameterisation]
+    """
+    Mapping from constitutive parameter name to its spatial parameterisation
+    """
+
+    metrics: list[IMetric]
+    """Virtual-work metrics used to evaluate candidate stress fields"""
+
+    objective_function: IObjectiveFunction
+    """Scalar or vector objective that aggregates metric results"""
+
+    optimiser: IOptimiser
+    """Optimisation algorithm that drives the parameter search"""
+
+
+@dataclass(slots=True)
+class IdentificationConfig:
+    """
+    Complete configuration for a VFM identification run.
+
+    Combines the constitutive law to identify, initial guessed for its
+    parameters, and one or more identification phases that
+    are executed sequentially, where the output of one phase
+    becomes the initial guess for the next
+    """
+
+    constitutive_law: IConstitutiveLaw
+    """Constitutive model whose parameters are being identified"""
+
+    parameters: dict[str, ConstitutiveParameter]
+    """Initial guess, lower bound, and upper bound for each parameter"""
+
+    phases: list[IdentificationPhase]
+    """Identification phases that are executed in order"""
diff --git a/src/pyvale/vfm/metric.py b/src/pyvale/vfm/metric.py
index 4799934fa..8c08b152b 100644
--- a/src/pyvale/vfm/metric.py
+++ b/src/pyvale/vfm/metric.py
@@ -9,6 +9,15 @@
 
 
 class IMetric(ABC):
+    """
+    Interface (abstract base class) for a virtual-work metric.
+
+    A metric evaluates the discrepancy between a candidate stress field and
+    the measured boundary conditions using the virtual work principle.
+    Multiple metrics (e.g. for different virtual fields) may be
+    combined in an objective function
+    """
+
     @abstractmethod
     def evaluate(
         self,
@@ -16,6 +25,27 @@ def evaluate(
         constitutive_law: IConstitutiveLaw,
         parameter_map_size: npt.NDArray[np.uint32],
         spatial_parameterisations: dict[str, ISpatialParameterisation],
-        experiment_data: ExperimentData
+        experiment_data: ExperimentData,
     ) -> npt.NDArray[np.float64]:
+        """
+        Evaluate the metric for a given stress candidate.
+
+        Parameters
+        ----------
+        stress : npt.NDArray[np.float64]
+            Candidate stress field, shape ``(timesteps, components, y, x)``
+        constitutive_law : IConstitutiveLaw
+            Constitutive law used to produce the stress
+        parameter_map_size : npt.NDArray[np.uint32]
+            Spatial dimensions ``(y, x)`` of the parameter maps
+        spatial_parameterisations : dict[str, ISpatialParameterisation]
+            Current spatial parameterisations keyed by parameter name
+        experiment_data : ExperimentData
+            Measured DIC data
+
+        Returns
+        -------
+        npt.NDArray[np.float64]
+            Metric value(s) per timestep or per spatial point
+        """
         pass
diff --git a/src/pyvale/vfm/metricsbvf.py b/src/pyvale/vfm/metricsbvf.py
index 9c9e4fa73..9af35631e 100644
--- a/src/pyvale/vfm/metricsbvf.py
+++ b/src/pyvale/vfm/metricsbvf.py
@@ -45,7 +45,7 @@ def __init__(
         self,
         x: npt.NDArray[np.float64],
         y: npt.NDArray[np.float64],
-        region_of_interest: npt.NDArray[np.uint32],
+        region_of_interest: npt.NDArray[np.bool_],
         edge_conditions: EdgeConditions,
         mesh_size: npt.NDArray[np.uint32],
         # TODO: option to adjust fraction of largest timesteps used for calculating VF scaling factor
diff --git a/src/pyvale/vfm/normalisation.py b/src/pyvale/vfm/normalisation.py
index be2958d04..52c2c0e6e 100644
--- a/src/pyvale/vfm/normalisation.py
+++ b/src/pyvale/vfm/normalisation.py
@@ -3,20 +3,45 @@
 
 from pyvale.vfm.dof import DegreeOfFreedom
 
-#TODO implement normalisation for various scaling types (e.g. linear, log etc)
-#TODO depending on above, could vectorise functions 
 
 def normalise_degree_of_freedom(
-    degree_of_freedom: DegreeOfFreedom
+    degree_of_freedom: DegreeOfFreedom,
 ) -> float:
+    """
+    Normalise a single degree of freedom to ``[0, 1]``.
+
+    Parameters
+    ----------
+    degree_of_freedom : DegreeOfFreedom
+        DOF with ``value``, ``lower_bound``, and ``upper_bound``
+
+    Returns
+    -------
+    float
+        Normalised value in ``[0, 1]``
+    """
     return (
         (degree_of_freedom.value - degree_of_freedom.lower_bound)
         / (degree_of_freedom.upper_bound - degree_of_freedom.lower_bound)
     )
 
+
 def normalise_degrees_of_freedom(
-    degrees_of_freedom: list[DegreeOfFreedom]
+    degrees_of_freedom: list[DegreeOfFreedom],
 ) -> npt.NDArray[np.float64]:
+    """
+    Normalise a list of degrees of freedom to ``[0, 1]``.
+
+    Parameters
+    ----------
+    degrees_of_freedom : list[DegreeOfFreedom]
+        One or more DOFs to normalise
+
+    Returns
+    -------
+    npt.NDArray[np.float64]
+        1D array of normalised values
+    """
     normalised_dofs = []
 
     for dof in degrees_of_freedom:
@@ -24,16 +49,52 @@ def normalise_degrees_of_freedom(
 
     return np.array(normalised_dofs)
 
+
 def denormalise_degree_of_freedom(
     normalised_value: float,
     lower_bound: float,
-    upper_bound: float
+    upper_bound: float,
 ) -> float:
+    """
+    Reverse the normalisation from ``[0, 1]`` back to physical units.
+
+    Parameters
+    ----------
+    normalised_value : float
+        Value in ``[0, 1]``
+    lower_bound : float
+        Physical lower bound
+    upper_bound : float
+        Physical upper bound
+
+    Returns
+    -------
+    float
+        Denormalised value in ``[lower_bound, upper_bound]``.
+    """
     return ((upper_bound - lower_bound) * normalised_value) + lower_bound
 
+
 def denormalise_degrees_of_freedom(
     normalised_values: npt.NDArray[np.float64],
     lower_bounds: npt.NDArray[np.float64],
-    upper_bounds: npt.NDArray[np.float64]
+    upper_bounds: npt.NDArray[np.float64],
 ) -> npt.NDArray[np.float64]:
+    """
+    Reverse the normalisation for an array of values.
+
+    Parameters
+    ----------
+    normalised_values : npt.NDArray[np.float64]
+        Values in ``[0, 1]``
+    lower_bounds : npt.NDArray[np.float64]
+        Physical lower bounds per DOF
+    upper_bounds : npt.NDArray[np.float64]
+        Physical upper bounds per DOF
+
+    Returns
+    -------
+    npt.NDArray[np.float64]
+        1D array of denormalised values
+    """
     return ((upper_bounds - lower_bounds) * normalised_values) + lower_bounds
diff --git a/src/pyvale/vfm/objectivefunc.py b/src/pyvale/vfm/objectivefunc.py
index 745c44b87..e95342bc8 100644
--- a/src/pyvale/vfm/objectivefunc.py
+++ b/src/pyvale/vfm/objectivefunc.py
@@ -5,21 +5,61 @@
 
 
 class IScalarObjectiveFunction(ABC):
+    """
+    Interface (abstract base class) for a scalar objective function.
+
+    Aggregates a list of metric results into a single scalar value that the
+    optimiser minimises
+    """
+
     @abstractmethod
     def evaluate(
         self,
         metric_results: list[npt.NDArray[np.float64]],
     ) -> float:
+        """
+        Aggregate metric results into a scalar cost
+
+        Parameters
+        ----------
+        metric_results : list[npt.NDArray[np.float64]]
+            One array per metric, each with the metric's output
+
+        Returns
+        -------
+        float
+            Scalar objective value to minimise
+        """
         pass
 
 
 class IVectorObjectiveFunction(ABC):
+    """
+    Interface (abstract base class) for a vector objective function.
+
+    Aggregates metric results into a vector that the optimiser minimises
+    """
+
     @abstractmethod
     def evaluate(
         self,
         metric_results: list[npt.NDArray[np.float64]],
     ) -> npt.NDArray[np.float64]:
+        """
+        Aggregate metric results into a residual vector.
+
+        Parameters
+        ----------
+        metric_results : list[npt.NDArray[np.float64]]
+            One array per metric, each with the metric's output
+
+        Returns
+        -------
+        npt.NDArray[np.float64]
+            Residual vector for the optimiser
+        """
         pass
 
 
 IObjectiveFunction = IScalarObjectiveFunction | IVectorObjectiveFunction
+"""Union type for objective functions that produce a scalar or vector cost"""
diff --git a/src/pyvale/vfm/optimiser.py b/src/pyvale/vfm/optimiser.py
index b46b59772..641c46c9d 100644
--- a/src/pyvale/vfm/optimiser.py
+++ b/src/pyvale/vfm/optimiser.py
@@ -14,8 +14,26 @@
 
 
 class IOptimiser(ABC):
-    # Run a set of optimisation passes until a best guess is found
-    # TODO: figure out what to return for different optimisers
+    """
+    Interface (abstract base class) for a VFM optimisation algorithm.
+
+    Drives the parameter search by repeatedly evaluating candidate stress
+    fields, computing metric/objective values, and updating the spatial
+    parameterisations until convergence
+    """
+
+    @abstractmethod
+    def get_required_objective_function_type(self) -> type:
+        """
+        Return the required objective function type for this optimiser.
+
+        Returns
+        -------
+        type
+            ``IScalarObjectiveFunction`` or ``IVectorObjectiveFunction``
+        """
+        pass
+
     @abstractmethod
     def optimise(
         self,
@@ -26,12 +44,32 @@ def optimise(
         objective_function: IObjectiveFunction,
         experiment_data: ExperimentData,
     ) -> dict[str, ISpatialParameterisation]:
+        """
+        Run the optimisation loop for one identification phase.
+
+        Parameters
+        ----------
+        constitutive_law : IConstitutiveLaw
+            Constitutive model whose parameters are being identified
+        parameter_map_size : npt.NDArray[np.uint32]
+            Spatial dimensions ``(y, x)`` of the parameter maps
+        spatial_parameterisations : dict[str, ISpatialParameterisation]
+            Initial parameter distributions keyed by parameter name
+        metrics : list[IMetric]
+            Virtual-work metrics to evaluate candidate stress fields
+        objective_function : IObjectiveFunction
+            Aggregates metric results into the quantity to be minimised
+        experiment_data : ExperimentData
+            Measured DIC data
+
+        Returns
+        -------
+        dict[str, ISpatialParameterisation]
+            Optimised spatial parameterisations after convergence
+        """
         pass
 
 
-# TODO: in some optimisers will the input be a scalar not a vector?
-# TODO: should we constrain return type? or have 2 functions for scalar/vector?
-#   or have the scalar representation just be a 1 elem vector?
 def evaluate_candidate(
     degrees_of_freedom: npt.NDArray[np.float64],
     constitutive_law: IConstitutiveLaw,
@@ -41,9 +79,38 @@ def evaluate_candidate(
     objective_function: IObjectiveFunction,
     experiment_data: ExperimentData,
 ) -> float | npt.NDArray[np.float64]:
+    """
+    Evaluate one candidate point in the design space.
+
+    Unpacks the normalised degrees of freedom into spatial parameterisations,
+    computes the resulting stress via the constitutive law, evaluates all
+    metrics, and aggregates them with the objective function.
+
+    Parameters
+    ----------
+    degrees_of_freedom : npt.NDArray[np.float64]
+        Normalised spatial parameterisation variables
+    constitutive_law : IConstitutiveLaw
+        Constitutive model
+    parameter_map_size : npt.NDArray[np.uint32]
+        Spatial dimensions ``(y, x)`` of the parameter maps
+    spatial_parameterisations : dict[str, ISpatialParameterisation]
+        Reference spatial parameterisations (cloned internally)
+    metrics : list[IMetric]
+        Virtual-work metrics
+    objective_function : IObjectiveFunction
+        Scalar or vector objective
+    experiment_data : ExperimentData
+        Measured DIC data
+
+    Returns
+    -------
+    float | npt.NDArray[np.float64]
+        Scalar or vector objective value for the candidate
+    """
     updated_spatial_parameterisations = unpack_spatial_parameterisations(
         spatial_parameterisations,
-        degrees_of_freedom
+        degrees_of_freedom,
     )
 
     updated_constitutive_parameter_maps = {
@@ -52,7 +119,7 @@ def evaluate_candidate(
     }
 
     updated_stress = constitutive_law.calculate_stress(
-        experiment_data.strain, updated_constitutive_parameter_maps
+        experiment_data.strain, updated_constitutive_parameter_maps,
     )
 
     metric_results = []
@@ -63,9 +130,8 @@ def evaluate_candidate(
                 constitutive_law,
                 parameter_map_size,
                 updated_spatial_parameterisations,
-                experiment_data
-            )
+                experiment_data,
+            ),
         )
 
     return objective_function.evaluate(metric_results)
-
diff --git a/src/pyvale/vfm/optimiserleastsquares.py b/src/pyvale/vfm/optimiserleastsquares.py
index 9dc1ca318..ace724418 100644
--- a/src/pyvale/vfm/optimiserleastsquares.py
+++ b/src/pyvale/vfm/optimiserleastsquares.py
@@ -6,7 +6,10 @@
 from pyvale.vfm.experimentdata import ExperimentData
 from pyvale.vfm.metric import IMetric
 from pyvale.vfm.normalisation import normalise_degrees_of_freedom
-from pyvale.vfm.objectivefunc import IObjectiveFunction
+from pyvale.vfm.objectivefunc import (
+    IObjectiveFunction,
+    IVectorObjectiveFunction,
+)
 from pyvale.vfm.optimiser import (
     IOptimiser,
     evaluate_candidate,
@@ -25,6 +28,9 @@
 #   if we need these, should treat the below as a dataclass and
 #   take these options as inputs in construction
 class LeastSquares(IOptimiser):
+    def get_required_objective_function_type(self) -> type:
+        return IVectorObjectiveFunction
+
     def optimise(
         self,
         constitutive_law: IConstitutiveLaw,
diff --git a/src/pyvale/vfm/spatialparam.py b/src/pyvale/vfm/spatialparam.py
index 2e77c70d4..cc4a267fc 100644
--- a/src/pyvale/vfm/spatialparam.py
+++ b/src/pyvale/vfm/spatialparam.py
@@ -9,41 +9,94 @@
 from pyvale.vfm.dof import DegreeOfFreedom
 
 
-# In general, spatial parameterisations should start out empty,
-# then the update_from_constituitive_parameter fills
-# the dofs with values
 class ISpatialParameterisation(ABC):
+    """
+    Interface (abstract base class) for a spatial parameterisation.
+
+    Maps constitutive parameter values onto the 2D DIC grid. Concrete
+    implementations define how the parameter varies in space and how
+    those spatial degrees of freedom are collected, updated, and converted
+    back to maps
+    """
+
     @abstractmethod
     def get_num_degrees_of_freedom(self) -> int:
+        """
+        Number of adjustable degrees of freedom in this parameterisation.
+
+        Returns
+        -------
+        int
+            Count of degrees of freedom
+        """
         pass
 
     @abstractmethod
-    def update_from_constitutive_parameter(
+    def initialise_from_constitutive_parameter(
         self,
-        constitutive_parameter: ConstitutiveParameter
+        constitutive_parameter: ConstitutiveParameter,
     ) -> None:
+        """
+        Initialise this parameterisation from a constitutive parameter.
+
+        Populates the internal degrees of freedom using the initial value
+        and bounds stored in the parameter.
+
+        Parameters
+        ----------
+        constitutive_parameter : ConstitutiveParameter
+            Parameter with initial value, lower bound, and upper bound
+        """
         pass
 
     @abstractmethod
     def to_map(
         self,
-        size: npt.NDArray[np.uint32]
+        size: npt.NDArray[np.uint32],
     ) -> npt.NDArray[np.float64]:
+        """
+        Generate a 2D map from the current parameterisation.
+
+        Parameters
+        ----------
+        size : npt.NDArray[np.uint32]
+            Target spatial shape ``(y, x)``
+
+        Returns
+        -------
+        npt.NDArray[np.float64]
+            Parameter map with shape ``(y, x)``
+        """
         pass
 
     @abstractmethod
     def collect_degrees_of_freedom(
         self,
     ) -> list[DegreeOfFreedom]:
+        """
+        Export a copy of the degrees of freedom for the optimiser.
+
+        Returns
+        -------
+        list[DegreeOfFreedom]
+            Copies of each degree of freedom
+        """
         pass
 
-    # We assume the order of the list/array is
-    # the same as the order we provided when we collected the dofs
     @abstractmethod
     def update_from_degrees_of_freedom(
         self,
-        degrees_of_freedom: list[DegreeOfFreedom] | npt.NDArray[np.float64]
+        degrees_of_freedom: list[DegreeOfFreedom] | npt.NDArray[np.float64],
     ) -> None:
+        """
+        Update internal state from optimiser-supplied degrees of freedom.
+
+        Parameters
+        ----------
+        degrees_of_freedom : list[DegreeOfFreedom] | npt.NDArray[np.float64]
+            New values for each degree of freedom, in the same order as
+            returned by `collect_degrees_of_freedom()`.
+        """
         pass
 
 
@@ -51,6 +104,25 @@ def unpack_spatial_parameterisations(
     reference_spatial_parameterisations: dict[str, ISpatialParameterisation],
     normalised_degrees_of_freedom: npt.NDArray[np.float64],
 ) -> dict[str, ISpatialParameterisation]:
+    """
+    Create a copy of spatial parameterisations with updated degrees of freedom.
+
+    Denormalises the degrees of freedom vector, copies each reference
+    parameterisation, and applies the new degrees of freedom to the copy.
+
+    Parameters
+    ----------
+    reference_spatial_parameterisations : dict[str, ISpatialParameterisation]
+        Reference parameterisations keyed by parameter name
+    normalised_degrees_of_freedom : npt.NDArray[np.float64]
+        Normalised degrees of freedom
+
+    Returns
+    -------
+    dict[str, ISpatialParameterisation]
+        A copy of the reference spatial parameterisations with
+        updated degrees of freedom
+    """
     lower_bounds = []
     upper_bounds = []
 
@@ -62,7 +134,7 @@ def unpack_spatial_parameterisations(
     degrees_of_freedom = denormalise_degrees_of_freedom(
         normalised_degrees_of_freedom,
         np.array(lower_bounds),
-        np.array(upper_bounds)
+        np.array(upper_bounds),
     )
 
     unpacked_spatial_parameterisations = {}
diff --git a/src/pyvale/vfm/spatialparambasisfuncs.py b/src/pyvale/vfm/spatialparambasisfuncs.py
index b628e8fc5..78983e678 100644
--- a/src/pyvale/vfm/spatialparambasisfuncs.py
+++ b/src/pyvale/vfm/spatialparambasisfuncs.py
@@ -69,7 +69,6 @@ def collect_degrees_of_freedom(
 
         return dofs
 
-    # TODO: length check to match num dofs
     def update_from_degrees_of_freedom(
         self,
         degrees_of_freedom: list[DegreeOfFreedom] | npt.NDArray[np.float64]
@@ -159,7 +158,6 @@ def collect_degrees_of_freedom(
 
         return dofs
 
-    # TODO: length check to match num dofs
     def update_from_degrees_of_freedom(
         self,
         degrees_of_freedom: list[DegreeOfFreedom] | npt.NDArray[np.float64]
@@ -241,7 +239,7 @@ def get_num_degrees_of_freedom(self) -> int:
         return num_dofs
 
     # TODO: create our initial basis functions with fitting
-    def update_from_constitutive_parameter(
+    def initialise_from_constitutive_parameter(
         self,
         constitutive_parameter: ConstitutiveParameter
     ) -> None:
@@ -253,6 +251,12 @@ def to_map(
         self,
         size: npt.NDArray[np.uint32]
     ) -> npt.NDArray[np.float64]:
+        if not self.kernels:
+            raise RuntimeError(
+                "self.kernels is empty, initialise_from_constitutive_parameter"
+                "must be called before to_map"
+            )
+
         ...
         # map = np.zeros((size[0], size[1]))
 
@@ -307,6 +311,12 @@ def update_from_degrees_of_freedom(
         self,
         degrees_of_freedom: list[DegreeOfFreedom] | npt.NDArray[np.float64]
     ) -> None:
+        if len(degrees_of_freedom) != self.get_num_degrees_of_freedom():
+            raise ValueError(
+                f"expected {self.get_num_degrees_of_freedom()} degrees of "
+                f"freedom, got {len(degrees_of_freedom)}"
+            )
+
         index = 0
 
         if isinstance(self.floor, DegreeOfFreedom):
diff --git a/src/pyvale/vfm/spatialparamhomogeneous.py b/src/pyvale/vfm/spatialparamhomogeneous.py
index 160b6dbf9..879341a5d 100644
--- a/src/pyvale/vfm/spatialparamhomogeneous.py
+++ b/src/pyvale/vfm/spatialparamhomogeneous.py
@@ -9,10 +9,7 @@
 from pyvale.vfm.spatialparam import ISpatialParameterisation
 
 
-# TODO: How do I flag whether value should be a dof or fixed?
-#   Maybe a constructor with a is_value_fixed var which informs
-#   how we should construct
-@dataclass(slots=True, init=False)
+@dataclass(slots=True)
 class HomogeneousSpatialParameterisation(ISpatialParameterisation):
     value: float | DegreeOfFreedom | None = None
 
@@ -22,23 +19,29 @@ def get_num_degrees_of_freedom(self) -> int:
         else:
             return 0
 
-    def update_from_constitutive_parameter(
+    def initialise_from_constitutive_parameter(
         self,
         constitutive_parameter: ConstitutiveParameter
     ) -> None:
-        if isinstance(self.value, DegreeOfFreedom):
+        if self.value is None or isinstance(self.value, DegreeOfFreedom):
             self.value =  DegreeOfFreedom(
-                float(np.mean(constitutive_parameter.value)),
+                float(np.mean(constitutive_parameter.map)),
                 constitutive_parameter.lower_bound,
                 constitutive_parameter.upper_bound,
             )
         else:
-            self.value = float(np.mean(constitutive_parameter.value))
+            self.value = float(np.mean(constitutive_parameter.map))
 
     def to_map(
         self,
         size: npt.NDArray[np.uint32]
     ) -> npt.NDArray[np.float64]:
+        if self.value is None:
+            raise RuntimeError(
+                "self.value is None, initialise_from_constitutive_parameter"
+                "must be called before to_map"
+            )
+
         if isinstance(self.value, DegreeOfFreedom):
             value = self.value.value
         else:
@@ -54,11 +57,16 @@ def collect_degrees_of_freedom(
         else:
             return []
 
-    # TODO: length check to match num dofs
     def update_from_degrees_of_freedom(
         self,
         degrees_of_freedom: list[DegreeOfFreedom] | npt.NDArray[np.float64]
     ) -> None:
+        if len(degrees_of_freedom) != self.get_num_degrees_of_freedom():
+            raise ValueError(
+                f"expected {self.get_num_degrees_of_freedom()} degrees of "
+                f"freedom, got {len(degrees_of_freedom)}"
+            )
+
         if isinstance(degrees_of_freedom, list):
             if isinstance(self.value, DegreeOfFreedom):
                 self.value = degrees_of_freedom[0]
diff --git a/src/pyvale/vfm/spatialparamknown.py b/src/pyvale/vfm/spatialparamknown.py
index 0a8ee0b41..de428eb2f 100644
--- a/src/pyvale/vfm/spatialparamknown.py
+++ b/src/pyvale/vfm/spatialparamknown.py
@@ -15,17 +15,27 @@ class KnownSpatialParameterisation(ISpatialParameterisation):
     def get_num_degrees_of_freedom(self) -> int:
         return 0
 
-    def update_from_constitutive_parameter(
+    def initialise_from_constitutive_parameter(
         self,
         constitutive_parameter: ConstitutiveParameter
     ) -> None:
-        self.value = constitutive_parameter.value
+        self.value = constitutive_parameter.map
 
     def to_map(
         self,
         size: npt.NDArray[np.uint32]
     ) -> npt.NDArray[np.float64]:
-        # TODO: value error if param value not the right size
+        if self.value is None:
+            raise RuntimeError(
+                "self.value is None, initialise_from_constitutive_parameter"
+                "must be called before to_map"
+            )
+
+        if self.value.shape != (size[0], size[1]):
+            raise ValueError(
+                f"parameter map shape {self.value.shape} does not match "
+                f"expected size ({size[0]}, {size[1]})"
+            )
 
         return self.value
 
diff --git a/src/pyvale/vfm/validation.py b/src/pyvale/vfm/validation.py
new file mode 100644
index 000000000..690a533b3
--- /dev/null
+++ b/src/pyvale/vfm/validation.py
@@ -0,0 +1,212 @@
+import numpy as np
+
+from pyvale.vfm.experimentdata import ExperimentData
+from pyvale.vfm.identificationconfig import IdentificationConfig
+from pyvale.vfm.spatialparamknown import KnownSpatialParameterisation
+
+
+def validate_experiment_data(
+    experiment_data: ExperimentData
+) -> None:
+    geometry = experiment_data.specimen_geometry
+    boundary_conditions = experiment_data.boundary_conditions
+
+    # Shape and dtype checks
+    strain = experiment_data.strain
+    if strain.ndim != 4:
+        raise ValueError(
+            f"strain must be 4D (timesteps, components, y, x), "
+            f"got ndim={strain.ndim}"
+        )
+    if strain.shape[1] != 3:
+        raise ValueError(
+            f"strain must have exactly 3 components [xx, yy, xy], "
+            f"got {strain.shape[1]}"
+        )
+    if strain.dtype != np.float64:
+        raise ValueError(
+            f"strain must be float64, got {strain.dtype}"
+        )
+
+    timesteps = experiment_data.timesteps
+    if timesteps.ndim != 1:
+        raise ValueError(
+            f"timesteps must be 1D, got ndim={timesteps.ndim}"
+        )
+    if timesteps.dtype != np.float64:
+        raise ValueError(
+            f"timesteps must be float64, got {timesteps.dtype}"
+        )
+
+    force = boundary_conditions.force
+    if force.ndim != 2:
+        raise ValueError(
+            f"force must be 2D (timesteps, 2) with columns [Fx, Fy], "
+            f"got ndim={force.ndim}"
+        )
+    if force.dtype != np.float64:
+        raise ValueError(
+            f"force must be float64, got {force.dtype}"
+        )
+
+    for field_name, array in [
+        ("x", geometry.x),
+        ("y", geometry.y),
+        ("pixel_area", geometry.pixel_area),
+    ]:
+        if array.ndim != 2:
+            raise ValueError(
+                f"{field_name} must be 2D (y, x), got ndim={array.ndim}"
+            )
+        if array.dtype != np.float64:
+            raise ValueError(
+                f"{field_name} must be float64, got {array.dtype}"
+            )
+
+    region_of_interest = geometry.region_of_interest
+    if region_of_interest.ndim != 2:
+        raise ValueError(
+            f"region_of_interest must be 2D (y, x), "
+            f"got ndim={region_of_interest.ndim}"
+        )
+    if region_of_interest.dtype != np.bool_:
+        raise ValueError(
+            f"region_of_interest must be bool dtype, "
+            f"got {region_of_interest.dtype}"
+        )
+
+    # Cross-field dimension agreement
+    n_timesteps, _, n_y, n_x = strain.shape
+
+    if timesteps.shape[0] != n_timesteps:
+        raise ValueError(
+            f"timesteps length ({timesteps.shape[0]}) does not match "
+            f"strain timesteps ({n_timesteps})"
+        )
+    if force.shape[0] != n_timesteps:
+        raise ValueError(
+            f"force timesteps ({force.shape[0]}) does not match "
+            f"strain timesteps ({n_timesteps})"
+        )
+
+    if geometry.x.shape != (n_y, n_x):
+        raise ValueError(
+            f"x shape {geometry.x.shape} does not match "
+            f"strain spatial dims ({n_y}, {n_x})"
+        )
+    if geometry.y.shape != (n_y, n_x):
+        raise ValueError(
+            f"y shape {geometry.y.shape} does not match "
+            f"strain spatial dims ({n_y}, {n_x})"
+        )
+    if region_of_interest.shape != (n_y, n_x):
+        raise ValueError(
+            f"region_of_interest shape {region_of_interest.shape} "
+            f"does not match strain spatial dims ({n_y}, {n_x})"
+        )
+    if geometry.pixel_area.shape != (n_y, n_x):
+        raise ValueError(
+            f"pixel_area shape {geometry.pixel_area.shape} does not match "
+            f"strain spatial dims ({n_y}, {n_x})"
+        )
+
+    # Value constraints
+    if np.any(geometry.x < 0):
+        raise ValueError("x coordinates must be non-negative")
+    if np.any(geometry.y < 0):
+        raise ValueError("y coordinates must be non-negative")
+    if geometry.thickness <= 0:
+        raise ValueError(
+            f"thickness must be positive, got {geometry.thickness}"
+        )
+    if np.any(geometry.pixel_area <= 0):
+        raise ValueError("pixel_area must be positive everywhere")
+    if np.any(np.diff(timesteps) <= 0):
+        raise ValueError("timesteps must be strictly increasing")
+    if not np.all(np.isfinite(force)):
+        raise ValueError("force contains NaN or Inf values")
+    # NaN is only allowed where region_of_interest is False (outside the mask)
+    flat_roi = region_of_interest.ravel()
+    flat_strain = strain.reshape(strain.shape[0], strain.shape[1], -1)
+    if not np.all(np.isfinite(flat_strain[:, :, flat_roi])):
+        raise ValueError(
+            "strain contains NaN or Inf within the region of interest"
+        )
+
+
+def validate_identification_config(
+    config: IdentificationConfig
+) -> None:
+    # Structure checks
+    if not config.phases:
+        raise ValueError("identification must have at least one phase")
+
+    if not config.parameters:
+        raise ValueError("identification must have at least one parameter")
+
+    for i, phase in enumerate(config.phases):
+        if not phase.metrics:
+            raise ValueError(
+                f"phase {i} must have at least one metric"
+            )
+
+    # Constitutive-law parameter requirements
+    required = set(config.constitutive_law.get_required_parameters())
+    given = set(config.parameters.keys())
+
+    extra = given - required
+    missing = required - given
+
+    if extra:
+        raise ValueError(f"unexpected parameter(s): {extra}")
+    if missing:
+        raise ValueError(f"missing required parameter(s): {missing}")
+
+    # Cross-field consistency: parameter name agreement
+    param_names = set(config.parameters.keys())
+
+    for i, phase in enumerate(config.phases):
+        phase_param_names = set(phase.spatial_parameterisations.keys())
+
+        if phase_param_names != param_names:
+            missing = param_names - phase_param_names
+            extra = phase_param_names - param_names
+            parts = []
+            if missing:
+                parts.append(f"missing: {missing}")
+            if extra:
+                parts.append(f"unknown: {extra}")
+            raise ValueError(
+                f"phase {i} spatial parameterisations do not match "
+                f"config parameters; {'; '.join(parts)}"
+            )
+
+        required_type = phase.optimiser.get_required_objective_function_type()
+        if not isinstance(phase.objective_function, required_type):
+            raise ValueError(
+                f"phase {i}: optimiser requires {required_type.__name__}, "
+                f"got {type(phase.objective_function).__name__}"
+            )
+
+        if all(
+            isinstance(sp, KnownSpatialParameterisation)
+            for sp in phase.spatial_parameterisations.values()
+        ):
+            raise ValueError(
+                f"phase {i}: all spatial parameterisations are "
+                f"KnownSpatialParameterisation; at least one"
+                f"must be identifiable"
+            )
+
+    # Value constraints
+    for name, param in config.parameters.items():
+        if not np.isfinite(param.lower_bound):
+            raise ValueError(
+                f"parameter '{name}': lower_bound ({param.lower_bound}) "
+                f"must be finite"
+            )
+        if not np.isfinite(param.upper_bound):
+            raise ValueError(
+                f"parameter '{name}': upper_bound ({param.upper_bound}) "
+                f"must be finite"
+            )
diff --git a/src/pyvale/vfm/vfm.py b/src/pyvale/vfm/vfm.py
deleted file mode 100644
index 96869aec8..000000000
--- a/src/pyvale/vfm/vfm.py
+++ /dev/null
@@ -1,47 +0,0 @@
-import numpy as np
-
-from pyvale.vfm.constlaw import EIdentificationType
-from pyvale.vfm.constparam import ConstitutiveParameter
-from pyvale.vfm.experimentdata import ExperimentData
-from pyvale.vfm.identification import Identification
-
-
-# TODO: config validation
-#   - no forward referencing in phases list
-#   - individual weights cant be greater than 1.0 in total
-#   - sum of weights must be 1.0
-#   - optimiser is compatible with objective function
-# TODO: think about io, no pickling
-def run_identification(
-    experiment_data: ExperimentData,
-    identification: Identification
-) -> dict[str, ConstitutiveParameter]:
-    match identification.constitutive_law.get_identification_type():
-        # TODO: implement linear case
-        case EIdentificationType.Linear:
-            ...
-        case EIdentificationType.Nonlinear:
-            parameter_map_size = np.array(
-                experiment_data.specimen_geometry.x.shape,
-                dtype=np.uint32
-            )
-
-            for phase in identification.phases:
-                for param_name, sp in phase.spatial_parameterisations.items():
-                    sp.update_from_constitutive_parameter(identification.parameters[param_name])
-
-                optimised_spatial_parameterisations = phase.optimiser.optimise(
-                    identification.constitutive_law,
-                    parameter_map_size,
-                    phase.spatial_parameterisations,
-                    phase.metrics,
-                    phase.objective_function,
-                    experiment_data
-                )
-
-                for param_name, sp in optimised_spatial_parameterisations.items():
-                    identification.parameters[param_name].value = (
-                        sp.to_map(parameter_map_size)
-                    )
-
-    return identification.parameters
diff --git a/src/pyvale/vfm/vfmesh.py b/src/pyvale/vfm/vfmesh.py
index acb60ffb2..c2b422db1 100644
--- a/src/pyvale/vfm/vfmesh.py
+++ b/src/pyvale/vfm/vfmesh.py
@@ -1110,7 +1110,7 @@ def plot_virtual_fields_mesh(
 def generate_virtual_fields_mesh(
     x: npt.NDArray[np.float64],
     y: npt.NDArray[np.float64],
-    specimen_mask: npt.NDArray[np.uint32],
+    specimen_mask: npt.NDArray[np.bool_],
     edge_conditions: EdgeConditions,
     mesh_size: npt.NDArray[np.uint32],
     use_nlgeom: bool = False,
diff --git a/tests/vfm/test_vfm_end_to_end.py b/tests/vfm/test_vfm_end_to_end.py
index d722abd77..139875e81 100644
--- a/tests/vfm/test_vfm_end_to_end.py
+++ b/tests/vfm/test_vfm_end_to_end.py
@@ -18,7 +18,7 @@
     ExperimentData,
     SpecimenGeometry,
 )
-from pyvale.vfm.identification import Identification, IdentificationPhase
+from pyvale.vfm.identification import IdentificationConfig, IdentificationPhase
 from pyvale.vfm.metrics.virtual_fields.sensitivity_based_virtual_fields import (
     SensitivityBasedVirtualFieldsMetric,
 )
@@ -110,7 +110,7 @@ def test_end_to_end(experiment_data: ExperimentData) -> None:
         )
     ]
 
-    identification = Identification(
+    identification = IdentificationConfig(
         IsotropicVonMisesElastoplasticity(
             LinearHardening()
         ),