sim_lib_physics_study/
study.rs

1use sha2::{Digest, Sha256};
2use sim_kernel::{
3    CapabilityName, Consistency, ContentId, Cx, EvalFabric, EvalMode, EvalRequest, Expr, ShapeRef,
4    Symbol, Value,
5};
6use sim_lib_numbers_stats::{
7    CoverageEvidence, LatinHypercubePlan, SampleDesign, SobolPlan, SweepPlan,
8};
9pub use sim_lib_numbers_stats::{Scramble, UntestedRegion};
10use sim_lib_physics_influence::SelectionInput;
11use std::{
12    collections::{BTreeMap, BTreeSet},
13    sync::Arc,
14    time::Duration,
15};
16
17/// Stable identity used for reviewable plan components and execution records.
18pub type Identity = ContentId;
19
20/// Which disjoint statistical lane owns a point.
21#[derive(Clone, Copy, Debug, Eq, Hash, Ord, PartialEq, PartialOrd)]
22pub enum Partition {
23    Fit,
24    Selection,
25    Test,
26}
27
28/// Unit-cube sampling policy. Grid counts are per axis.
29#[derive(Clone, Debug, Eq, PartialEq)]
30pub enum SamplerPolicy {
31    Grid {
32        counts: Vec<usize>,
33    },
34    Latin {
35        points: usize,
36    },
37    Sobol {
38        points: usize,
39        skip: u64,
40        scramble: Scramble,
41    },
42}
43
44/// Semantic spacing applied after canonical unit-cube generation.
45#[derive(Clone, Copy, Debug, Eq, PartialEq)]
46pub enum Spacing {
47    Linear,
48    Logarithmic,
49}
50
51/// One exact injection, expressed in physical coordinates.
52#[derive(Clone)]
53pub struct BoundaryInjection {
54    pub label: String,
55    pub values: Vec<f64>,
56}
57
58/// A semantic quantity axis. `shape` is the canonical runtime Shape.
59#[derive(Clone)]
60pub struct ParameterAxis {
61    pub name: String,
62    pub quantity_shape: ShapeRef,
63    pub inclusive_bounds: (f64, f64),
64    pub spacing: Spacing,
65}
66
67/// Deterministic partition assignment. Every generated ordinal is assigned once.
68#[derive(Clone, Debug, Eq, PartialEq)]
69pub struct PartitionPlan {
70    pub fit: BTreeSet<usize>,
71    pub selection: BTreeSet<usize>,
72    pub test: BTreeSet<usize>,
73}
74
75/// Placement request remains data; EvalFabric resolves the provider late.
76#[derive(Clone, Debug, Eq, PartialEq)]
77pub struct PlacementRequest {
78    pub target: Symbol,
79    pub consistency: Consistency,
80    pub required_capabilities: Vec<CapabilityName>,
81}
82
83#[derive(Clone, Debug, Eq, PartialEq)]
84pub struct StudyLimits {
85    pub max_points: usize,
86    pub max_failures: usize,
87    pub deadline: Option<Duration>,
88}
89
90#[derive(Clone, Debug, Eq, PartialEq)]
91pub struct AuditPolicy {
92    pub energy_after_selection_only: bool,
93    pub record_residuals: bool,
94}
95
96#[derive(Clone, Debug, Eq, PartialEq)]
97pub struct InfluencePolicy {
98    pub require_clean_selection: bool,
99}
100
101/// All inputs to a study, immutable after validation.
102#[derive(Clone)]
103pub struct StudyPlan {
104    pub id: Identity,
105    pub model_id: Identity,
106    pub boundary_event_graph_id: Identity,
107    pub initial_state_id: Identity,
108    pub method_plan_ids: Vec<Identity>,
109    pub refinement_plan_ids: Vec<Identity>,
110    pub audit_policy: AuditPolicy,
111    pub influence_policy: InfluencePolicy,
112    pub outputs: Vec<String>,
113    pub limits: StudyLimits,
114    pub seed: u64,
115    pub placement: PlacementRequest,
116    pub axes: Vec<ParameterAxis>,
117    pub sampler: SamplerPolicy,
118    pub boundary_injections: Vec<BoundaryInjection>,
119    pub partitions: PartitionPlan,
120    pub untested_regions: Vec<UntestedRegion>,
121}
122
123#[derive(Clone, Debug, Eq, PartialEq)]
124pub enum PlanError {
125    Empty(&'static str),
126    InvalidAxis(String),
127    InvalidBoundary(String),
128    PartitionOverlap(usize),
129    MissingPartition(usize),
130    WorkLimit { requested: usize, limit: usize },
131    Sampler(String),
132}
133
134impl StudyPlan {
135    /// Validates the full immutable binding before sampling or execution.
136    pub fn validate(&self) -> Result<(), PlanError> {
137        if self.axes.is_empty() {
138            return Err(PlanError::Empty("axes"));
139        }
140        if self.outputs.is_empty() {
141            return Err(PlanError::Empty("outputs"));
142        }
143        if self.method_plan_ids.is_empty() {
144            return Err(PlanError::Empty("method plans"));
145        }
146        if self.refinement_plan_ids.is_empty() {
147            return Err(PlanError::Empty("refinement plans"));
148        }
149        if self.untested_regions.is_empty() {
150            return Err(PlanError::Empty("known untested regions"));
151        }
152        if !self.audit_policy.energy_after_selection_only {
153            return Err(PlanError::Empty("energy audit must be post-selection"));
154        }
155        for axis in &self.axes {
156            let (lo, hi) = axis.inclusive_bounds;
157            if axis.name.is_empty()
158                || !lo.is_finite()
159                || !hi.is_finite()
160                || lo > hi
161                || (axis.spacing == Spacing::Logarithmic && lo <= 0.0)
162            {
163                return Err(PlanError::InvalidAxis(axis.name.clone()));
164            }
165        }
166        for boundary in &self.boundary_injections {
167            if boundary.label.is_empty()
168                || boundary.values.len() != self.axes.len()
169                || boundary.values.iter().zip(&self.axes).any(|(v, a)| {
170                    !v.is_finite() || *v < a.inclusive_bounds.0 || *v > a.inclusive_bounds.1
171                })
172            {
173                return Err(PlanError::InvalidBoundary(boundary.label.clone()));
174            }
175        }
176        for index in self
177            .partitions
178            .fit
179            .iter()
180            .chain(&self.partitions.selection)
181            .chain(&self.partitions.test)
182        {
183            let lanes = usize::from(self.partitions.fit.contains(index))
184                + usize::from(self.partitions.selection.contains(index))
185                + usize::from(self.partitions.test.contains(index));
186            if lanes > 1 {
187                return Err(PlanError::PartitionOverlap(*index));
188            }
189        }
190        Ok(())
191    }
192
193    /// Generates points through the canonical stats designs and sweep evidence.
194    pub fn design(&self) -> Result<StudyDesign, PlanError> {
195        self.validate()?;
196        let dimensions = self.axes.len();
197        let base = match &self.sampler {
198            SamplerPolicy::Grid { counts } => {
199                grid_design(dimensions, counts, &self.untested_regions)?
200            }
201            SamplerPolicy::Latin { points } => LatinHypercubePlan {
202                dimensions,
203                points: *points,
204                seed: self.seed,
205                max_work: self.limits.max_points.saturating_mul(dimensions) as u64,
206                untested_regions: self.untested_regions.clone(),
207            }
208            .generate()
209            .map_err(|e| PlanError::Sampler(format!("{e:?}")))?,
210            SamplerPolicy::Sobol {
211                points,
212                skip,
213                scramble,
214            } => SobolPlan {
215                dimensions,
216                points: *points,
217                skip: *skip,
218                scramble: *scramble,
219                seed: self.seed,
220                max_work: self.limits.max_points.saturating_mul(dimensions) as u64,
221                untested_regions: self.untested_regions.clone(),
222            }
223            .generate()
224            .map_err(|e| PlanError::Sampler(format!("{e:?}")))?,
225        };
226        // Canonical stats owner computes exact duplicate evidence for all designs.
227        let mut sampled = SweepPlan {
228            inject_lower_boundary: false,
229            inject_upper_boundary: false,
230            untested_regions: vec![],
231        }
232        .apply(base);
233        let mut physical = sampled
234            .points
235            .iter()
236            .map(|p| scale_point(p, &self.axes))
237            .collect::<Vec<_>>();
238        let mut injected = Vec::new();
239        for boundary in &self.boundary_injections {
240            if let Some(existing) = physical
241                .iter()
242                .position(|p| exact_point_eq(p, &boundary.values))
243            {
244                injected.push((boundary.label.clone(), existing, true));
245            } else {
246                let index = physical.len();
247                physical.push(boundary.values.clone());
248                sampled
249                    .points
250                    .push(unscale_point(&boundary.values, &self.axes));
251                injected.push((boundary.label.clone(), index, false));
252            }
253        }
254        if physical.len() > self.limits.max_points {
255            return Err(PlanError::WorkLimit {
256                requested: physical.len(),
257                limit: self.limits.max_points,
258            });
259        }
260        sampled.coverage.boundary_injections = injected.iter().map(|(_, i, _)| *i).collect();
261        sampled.coverage.duplicates = exact_duplicates(&physical);
262        let mut points = Vec::with_capacity(physical.len());
263        for (ordinal, coordinates) in physical.into_iter().enumerate() {
264            let partition = partition_for(&self.partitions, ordinal)
265                .ok_or(PlanError::MissingPartition(ordinal))?;
266            points.push(StudyPoint {
267                id: point_id(&self.id, ordinal, &coordinates),
268                ordinal,
269                coordinates,
270                partition,
271            });
272        }
273        Ok(StudyDesign {
274            plan_id: self.id.clone(),
275            points,
276            coverage: sampled.coverage,
277            boundary_evidence: injected,
278        })
279    }
280}
281
282fn partition_for(plan: &PartitionPlan, i: usize) -> Option<Partition> {
283    if plan.fit.contains(&i) {
284        Some(Partition::Fit)
285    } else if plan.selection.contains(&i) {
286        Some(Partition::Selection)
287    } else if plan.test.contains(&i) {
288        Some(Partition::Test)
289    } else {
290        None
291    }
292}
293
294fn grid_design(
295    dimensions: usize,
296    counts: &[usize],
297    untested: &[UntestedRegion],
298) -> Result<SampleDesign, PlanError> {
299    if counts.len() != dimensions || counts.contains(&0) {
300        return Err(PlanError::Sampler(
301            "grid counts must match nonzero axes".into(),
302        ));
303    }
304    let count = counts
305        .iter()
306        .try_fold(1usize, |n, v| n.checked_mul(*v))
307        .ok_or(PlanError::Sampler("grid size overflow".into()))?;
308    let mut points = Vec::with_capacity(count);
309    for ordinal in 0..count {
310        let mut remainder = ordinal;
311        let mut point = Vec::with_capacity(dimensions);
312        for cells in counts {
313            let cell = remainder % *cells;
314            remainder /= *cells;
315            point.push(if *cells == 1 {
316                0.5
317            } else {
318                cell as f64 / (*cells - 1) as f64
319            });
320        }
321        points.push(point);
322    }
323    Ok(SampleDesign {
324        points,
325        coverage: CoverageEvidence {
326            sequence_identity: format!("cartesian-grid/inclusive-v1;counts={counts:?}"),
327            boundary_injections: vec![],
328            duplicates: vec![],
329            stratum_occupancy: counts.iter().map(|n| vec![1; *n]).collect(),
330            sampler: None,
331            untested_regions: untested.to_vec(),
332            work: count as u64,
333        },
334    })
335}
336
337fn scale_point(unit: &[f64], axes: &[ParameterAxis]) -> Vec<f64> {
338    unit.iter()
339        .zip(axes)
340        .map(|(u, a)| {
341            let (lo, hi) = a.inclusive_bounds;
342            match a.spacing {
343                Spacing::Linear => lo + u * (hi - lo),
344                Spacing::Logarithmic => (lo.ln() + u * (hi.ln() - lo.ln())).exp(),
345            }
346        })
347        .collect()
348}
349fn unscale_point(values: &[f64], axes: &[ParameterAxis]) -> Vec<f64> {
350    values
351        .iter()
352        .zip(axes)
353        .map(|(v, a)| {
354            let (lo, hi) = a.inclusive_bounds;
355            if lo == hi {
356                0.0
357            } else {
358                match a.spacing {
359                    Spacing::Linear => (v - lo) / (hi - lo),
360                    Spacing::Logarithmic => (v.ln() - lo.ln()) / (hi.ln() - lo.ln()),
361                }
362            }
363        })
364        .collect()
365}
366fn exact_point_eq(a: &[f64], b: &[f64]) -> bool {
367    a.len() == b.len() && a.iter().zip(b).all(|(x, y)| x.to_bits() == y.to_bits())
368}
369fn exact_duplicates(points: &[Vec<f64>]) -> Vec<(usize, usize)> {
370    let mut out = vec![];
371    for later in 0..points.len() {
372        if let Some(earlier) = (0..later).find(|i| exact_point_eq(&points[*i], &points[later])) {
373            out.push((later, earlier));
374        }
375    }
376    out
377}
378
379fn identity(parts: impl IntoIterator<Item = impl AsRef<[u8]>>) -> Identity {
380    let mut h = Sha256::new();
381    for part in parts {
382        let p = part.as_ref();
383        h.update((p.len() as u64).to_be_bytes());
384        h.update(p);
385    }
386    ContentId::from_bytes(Symbol::qualified("core", "sha256"), h.finalize().into())
387}
388fn point_id(plan: &Identity, ordinal: usize, values: &[f64]) -> Identity {
389    let mut parts = Vec::new();
390    parts.push(plan.bytes.to_vec());
391    parts.push((ordinal as u64).to_be_bytes().to_vec());
392    for value in values {
393        parts.push(value.to_bits().to_be_bytes().to_vec());
394    }
395    identity(parts)
396}
397
398/// Generated point set plus exact sampler and coverage evidence.
399pub struct StudyDesign {
400    pub plan_id: Identity,
401    pub points: Vec<StudyPoint>,
402    pub coverage: CoverageEvidence,
403    pub boundary_evidence: Vec<(String, usize, bool)>,
404}
405#[derive(Clone, Debug, PartialEq)]
406pub struct StudyPoint {
407    pub id: Identity,
408    pub ordinal: usize,
409    pub coordinates: Vec<f64>,
410    pub partition: Partition,
411}
412
413/// Caller-owned expression builder; study owns orchestration, not model encoding.
414pub trait PointExpression: Send + Sync {
415    fn expression(&self, plan: &StudyPlan, point: &StudyPoint) -> sim_kernel::Result<Expr>;
416}
417
418#[derive(Clone, Debug, Eq, PartialEq)]
419pub enum IncompleteOutcome {
420    Refused(String),
421    Interrupted(String),
422    Unknown(String),
423}
424
425/// One placement-transparent execution outcome.
426pub enum PointOutcome {
427    Complete(Value),
428    Incomplete(IncompleteOutcome),
429}
430
431pub struct PointResult {
432    pub id: Identity,
433    pub plan_id: Identity,
434    pub point_id: Identity,
435    pub reference_execution_id: Identity,
436    pub provider_id: Identity,
437    pub partition: Partition,
438    pub outcome: PointOutcome,
439}
440
441/// Executes points sequentially only to submit independent ordinary EvalFabric requests.
442/// Concurrency, transport, retry, and site choice remain exclusively fabric concerns.
443pub fn realize_sweep(
444    cx: &mut Cx,
445    fabric: &dyn EvalFabric,
446    plan: &StudyPlan,
447    design: &StudyDesign,
448    builder: &dyn PointExpression,
449    provider_id: Identity,
450) -> Vec<PointResult> {
451    let mut results = Vec::with_capacity(design.points.len());
452    for point in &design.points {
453        let execution_id = identity([
454            plan.id.bytes.as_slice(),
455            point.id.bytes.as_slice(),
456            provider_id.bytes.as_slice(),
457        ]);
458        let outcome = match builder.expression(plan, point) {
459            Err(error) => PointOutcome::Incomplete(IncompleteOutcome::Refused(error.to_string())),
460            Ok(expr) => match fabric.realize(
461                cx,
462                EvalRequest {
463                    expr,
464                    result_shape: None,
465                    required_capabilities: plan.placement.required_capabilities.clone(),
466                    deadline: plan.limits.deadline,
467                    consistency: plan.placement.consistency,
468                    mode: EvalMode::Eval,
469                    answer_limit: None,
470                    stream_buffer: None,
471                    stream: false,
472                    trace: true,
473                },
474            ) {
475                Ok(reply) => PointOutcome::Complete(reply.value),
476                Err(error) => PointOutcome::Incomplete(classify_failure(error.to_string())),
477            },
478        };
479        let result_id = identity([execution_id.bytes.as_slice(), outcome_tag(&outcome)]);
480        results.push(PointResult {
481            id: result_id,
482            plan_id: plan.id.clone(),
483            point_id: point.id.clone(),
484            reference_execution_id: execution_id,
485            provider_id: provider_id.clone(),
486            partition: point.partition,
487            outcome,
488        });
489    }
490    results
491}
492fn classify_failure(message: String) -> IncompleteOutcome {
493    let lower = message.to_ascii_lowercase();
494    if lower.contains("interrupt") || lower.contains("cancel") {
495        IncompleteOutcome::Interrupted(message)
496    } else if lower.contains("refus") || lower.contains("capab") || lower.contains("limit") {
497        IncompleteOutcome::Refused(message)
498    } else {
499        IncompleteOutcome::Unknown(message)
500    }
501}
502fn outcome_tag(outcome: &PointOutcome) -> &'static [u8] {
503    match outcome {
504        PointOutcome::Complete(_) => b"complete",
505        PointOutcome::Incomplete(IncompleteOutcome::Refused(_)) => b"refused",
506        PointOutcome::Incomplete(IncompleteOutcome::Interrupted(_)) => b"interrupted",
507        PointOutcome::Incomplete(IncompleteOutcome::Unknown(_)) => b"unknown",
508    }
509}
510
511/// Clean selection proof and candidates. There is deliberately no energy input.
512pub fn select<'a>(
513    proof: &SelectionInput,
514    candidates: impl IntoIterator<Item = &'a PointResult>,
515    score: impl Fn(&PointResult) -> Option<f64>,
516) -> Option<&'a PointResult> {
517    let _proof_identity = proof.proof_identity();
518    candidates
519        .into_iter()
520        .filter(|r| r.partition == Partition::Selection)
521        .filter_map(|r| score(r).filter(|v| v.is_finite()).map(|v| (r, v)))
522        .min_by(|a, b| a.1.total_cmp(&b.1))
523        .map(|(r, _)| r)
524}
525
526/// Post-selection observations. These can annotate but cannot choose a result.
527#[derive(Clone, Debug, Default, PartialEq)]
528pub struct Observations {
529    pub energy_store_change: Option<f64>,
530    pub unexplained_residual: Option<f64>,
531    pub sign_reversal: bool,
532    pub zero_crossing: bool,
533    pub event_transition: bool,
534    pub topology_transition: bool,
535    pub initial_state_sensitivity: Option<f64>,
536    pub boundary_value: bool,
537    pub fold: bool,
538    pub method_sensitivity: Option<f64>,
539    pub tolerance_sensitivity: Option<f64>,
540    pub precision_sensitivity: Option<f64>,
541}
542
543#[derive(Clone, Debug, Default)]
544pub struct ObservationInput {
545    pub energy_before: Option<f64>,
546    pub energy_after: Option<f64>,
547    pub residual: Option<f64>,
548    pub previous: Option<f64>,
549    pub current: Option<f64>,
550    pub event_changed: bool,
551    pub topology_changed: bool,
552    pub initial_delta: Option<f64>,
553    pub boundary_value: bool,
554    pub slope_before: Option<f64>,
555    pub slope_after: Option<f64>,
556    pub method_delta: Option<f64>,
557    pub tolerance_delta: Option<f64>,
558    pub precision_delta: Option<f64>,
559}
560pub fn detect(input: &ObservationInput) -> Observations {
561    let sign_reversal = matches!((input.previous,input.current),(Some(a),Some(b)) if a.signum()!=b.signum() && a!=0.0 && b!=0.0);
562    let zero_crossing = matches!((input.previous,input.current),(Some(a),Some(b)) if a==0.0 || b==0.0 || a.signum()!=b.signum());
563    let fold = matches!((input.slope_before,input.slope_after),(Some(a),Some(b)) if a.signum()!=b.signum());
564    Observations {
565        energy_store_change: input
566            .energy_before
567            .zip(input.energy_after)
568            .map(|(a, b)| b - a),
569        unexplained_residual: input.residual,
570        sign_reversal,
571        zero_crossing,
572        event_transition: input.event_changed,
573        topology_transition: input.topology_changed,
574        initial_state_sensitivity: input.initial_delta.map(f64::abs),
575        boundary_value: input.boundary_value,
576        fold,
577        method_sensitivity: input.method_delta.map(f64::abs),
578        tolerance_sensitivity: input.tolerance_delta.map(f64::abs),
579        precision_sensitivity: input.precision_delta.map(f64::abs),
580    }
581}
582
583/// Replay is exact only under the same reference execution identity.
584pub fn exact_replay(left: &PointResult, right: &PointResult) -> bool {
585    left.reference_execution_id == right.reference_execution_id && left.id == right.id
586}
587
588#[derive(Clone, Debug, PartialEq)]
589pub struct ProviderEnvelope {
590    pub provider_id: Identity,
591    pub absolute_tolerance: f64,
592    pub relative_tolerance: f64,
593}
594#[derive(Clone, Debug, PartialEq)]
595pub struct EnvelopeComparison {
596    pub equivalent: bool,
597    pub absolute_delta: f64,
598    pub allowed_delta: f64,
599}
600pub fn compare_provider_values(
601    left: f64,
602    right: f64,
603    left_provider: &ProviderEnvelope,
604    right_provider: &ProviderEnvelope,
605) -> EnvelopeComparison {
606    let absolute_delta = (right - left).abs();
607    let scale = left.abs().max(right.abs());
608    let allowed_delta = left_provider
609        .absolute_tolerance
610        .max(right_provider.absolute_tolerance)
611        + scale
612            * left_provider
613                .relative_tolerance
614                .max(right_provider.relative_tolerance);
615    EnvelopeComparison {
616        equivalent: absolute_delta <= allowed_delta,
617        absolute_delta,
618        allowed_delta,
619    }
620}
621
622/// Summary retains all partial outcomes and never manufactures retry safety.
623pub struct StudyResults {
624    pub plan_id: Identity,
625    pub results: Vec<PointResult>,
626    pub observations: BTreeMap<Identity, Observations>,
627    pub retry_safe: BTreeSet<Identity>,
628}
629impl StudyResults {
630    pub fn new(plan_id: Identity, results: Vec<PointResult>) -> Self {
631        Self {
632            plan_id,
633            results,
634            observations: BTreeMap::new(),
635            retry_safe: BTreeSet::new(),
636        }
637    }
638    pub fn annotate_selected(&mut self, result: &PointResult, observations: Observations) {
639        self.observations.insert(result.id.clone(), observations);
640    }
641}
642
643/// Runtime-facing loadable library identity and stable callable names.
644pub fn study_lib_symbol() -> Symbol {
645    Symbol::qualified("physics-study", "lib")
646}
647pub fn study_surface_symbols() -> [Symbol; 4] {
648    [
649        Symbol::qualified("physics-study", "design"),
650        Symbol::qualified("physics-study", "realize"),
651        Symbol::qualified("physics-study", "select"),
652        Symbol::qualified("physics-study", "compare-providers"),
653    ]
654}
655
656/// Shared fabric handle accepted by host adapters without exposing transport.
657pub type StudyFabric = Arc<dyn EvalFabric>;
sim_lib_physics_study/study.rs

sim_lib_physics_study/
study.rs
