🎼 Fugue User Guide

A production-ready, monadic probabilistic programming library for Rust

Write elegant probabilistic programs by composing Model values in direct style; execute them with pluggable interpreters and state-of-the-art inference algorithms.

Supported Rust: 1.70+ • Platforms: Linux / macOS / Windows • Crate: fugue-ppl on crates.io

About Fugue

• 🧩 Monadic PPL: Compose probabilistic programs using pure functional abstractions
• 🔒 Type-Safe Distributions: 10+ built-in probability distributions with natural return types
• 📊 Multiple Inference Methods: MCMC, SMC, Variational Inference, ABC
• 🔍 Comprehensive Diagnostics: R-hat convergence, effective sample size, validation
• 🚀 Production Ready: Numerically stable algorithms with memory optimization
• ✨ Ergonomic Macros: Do-notation (prob!), vectorization (plate!), addressing (addr!)

Installation

[dependencies]
fugue-ppl = "0.1.0"

Getting Started with Fugue

Welcome to Fugue, a type-safe probabilistic programming library for Rust! This guide will get you building Bayesian models in just 15-20 minutes.

Learning Path

We recommend following this path for the best learning experience:

flowchart LR
    A[Installation<br/>2 min] --> B[Your First Model<br/>5 min]
    B --> C[Understanding Models<br/>8 min]
    C --> D[Running Inference<br/>5 min]
    D --> E[Complete Tutorials<br/>45-60 min each]

Quick Start

If you're impatient and want to see Fugue in action immediately:

# Create a new project
cargo new my_bayesian_project
cd my_bayesian_project

# Add Fugue
cargo add fugue-ppl rand

# Copy our "Hello, Probabilistic World!" example into src/main.rs
# (See Installation section)

# Run it!
cargo run

What Makes Fugue Different?

🔒 Type Safety First

Unlike other probabilistic programming libraries, Fugue preserves natural types:

// In Fugue ✅
let coin: bool = sample(addr!("coin"), Bernoulli::new(0.5).unwrap());
let count: u64 = sample(addr!("events"), Poisson::new(3.0).unwrap());
let category: usize = sample(addr!("choice"), Categorical::uniform(5).unwrap());

// Other PPLs ❌
let coin: f64 = sample("coin", Bernoulli(0.5));  // Returns 0.0 or 1.0
let count: f64 = sample("events", Poisson(3.0)); // Need to cast to int
// let category: f64 = sample("choice", Categorical([...])); // Risky indexing

🚀 Zero-Cost Abstractions

Models compile to efficient code with no runtime overhead.

🧰 Composable Architecture

Separate model specification from execution strategy through handlers.

📊 Production Ready

Built-in diagnostics, memory optimization, and error handling.

Architecture Overview

Fugue's modular design separates concerns cleanly:

graph TB
    subgraph "Your Code"
        M[Model Definition]
        D[Data & Observations]
    end

    subgraph "Core System"
        C[Distributions & Types]
        H[Handlers & Interpreters]
        T[Traces & Memory]
    end

    subgraph "Inference Engines"
        MCMC[MCMC Sampling]
        SMC[Particle Filtering]
        VI[Variational Inference]
        ABC[ABC Methods]
    end

    M --> C
    D --> C
    C --> H
    H --> T
    T --> MCMC
    T --> SMC
    T --> VI
    T --> ABC

The Big Picture

Probabilistic Programming lets you:

1. Model uncertainty and relationships in data
2. Condition on observations to learn parameters
3. Infer posterior distributions and make predictions
4. Quantify uncertainty in your conclusions

Fugue makes this safe, fast, and composable in Rust.

Next Steps

Ready to dive in?

After completing Getting Started, explore:

• Complete Tutorials - End-to-end projects with real applications
• How-To Guides - Specific techniques and best practices
• API Documentation - Comprehensive technical reference

Prerequisites: Basic Rust knowledge (variables, functions, cargo commands)

Installation

Getting Fugue set up in your Rust project takes just 2 minutes. Let's get you running!

# Install Rust via rustup
curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh

# Update to latest stable
rustup update stable

Adding Fugue to Your Project

New Project

cargo new my_probabilistic_project
cd my_probabilistic_project

Add Fugue to your Cargo.toml:

[dependencies]
fugue-ppl = "0.1.0"
rand = "0.8"  # For random number generation

Existing Project

Add Fugue to your existing Cargo.toml:

[dependencies]
fugue-ppl = "0.1.0"
rand = "0.8"

Or use cargo add:

cargo add fugue-ppl rand

Verification: "Hello, Probabilistic World!"

Let's verify your installation with a simple example that showcases Fugue's type safety.

Create or replace src/main.rs:

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn main() {
    println!("🎲 Hello, Probabilistic World!");

    // Create a simple model: flip a biased coin
    let coin_model = sample(addr!("coin"), Bernoulli::new(0.7).unwrap());

    // Run the model with a seeded RNG for reproducible results
    let mut rng = StdRng::seed_from_u64(42);
    let (is_heads, trace) = runtime::handler::run(
        runtime::interpreters::PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        coin_model,
    );

    // Print the result - notice it's a bool, not a float!
    let result = if is_heads { "Heads" } else { "Tails" };
    println!("🪙 Coin flip result: {}", result);
    println!("📊 Log probability: {:.4}", trace.total_log_weight());

    // Demonstrate type safety with different distributions
    let mut rng = StdRng::seed_from_u64(123);

    // Count events - returns u64 directly
    let (event_count, _) = runtime::handler::run(
        runtime::interpreters::PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        sample(addr!("events"), Poisson::new(3.5).unwrap()),
    );
    println!("🎯 Event count: {} (type: u64)", event_count);

    // Choose category - returns usize for safe indexing
    let options = vec!["red", "green", "blue"];
    let (category_idx, _) = runtime::handler::run(
        runtime::interpreters::PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        sample(addr!("color"), Categorical::uniform(3).unwrap()),
    );
    println!("🎨 Chosen color: {} (safe indexing!)", options[category_idx]);

    println!("✅ Fugue is working correctly!");
}

Run it to verify everything works:

cargo run

You should see output like:

🎲 Hello, Probabilistic World!
🪙 Coin flip result: Heads
📊 Log probability: -0.3567
🎯 Event count: 4 (type: u64)
🎨 Chosen color: blue (safe indexing!)
✅ Fugue is working correctly!

IDE Setup

VS Code

Install the rust-analyzer extension for the best development experience:

1. Open VS Code
2. Go to Extensions (Ctrl+Shift+X)
3. Search for "rust-analyzer"
4. Install the official rust-analyzer extension

Other IDEs

• IntelliJ/CLion: Install the Rust plugin
• Vim/Neovim: Use coc-rust-analyzer or native LSP
• Emacs: Use lsp-mode with rust-analyzer

Optional: Running Examples

Fugue comes with comprehensive examples to explore:

# Clone the repository to access examples
git clone https://github.com/your-org/fugue-ppl
cd fugue-ppl

# List available examples
ls examples/

# Run a simple example
cargo run --example gaussian_mean -- --obs 2.5 --seed 42

# Try a more complex one
cargo run --example working_with_distributions

Troubleshooting

Common Issues

Build fails with dependency errors:

# Make sure you're using Rust 1.70+
rustc --version

# Update your dependencies
cargo update

Examples don't run:

# Make sure you're in the project root directory
pwd

# Check example names
ls examples/

IDE doesn't provide completions:

• Make sure rust-analyzer is installed and running
• Try restarting your IDE after installing dependencies
• Check that your Cargo.toml has the correct dependencies

Getting Help

If you encounter issues:

1. Check the GitHub Issues
2. Review the examples for working code
3. Read the API documentation

Next Steps

Installation complete! 🎉

Ready to build your first probabilistic model? → Your First Model

Want to explore examples first? → Complete Tutorials

Time: ~2 minutes • Next: Your First Model

Your First Model

Now that Fugue is installed, let's build your first probabilistic model step by step. We'll start simple and gradually introduce the key concepts.

Step 1: The Simplest Model

Let's start with the simplest possible model - one that always returns the same value:

use fugue::*;

fn constant_model() -> Model<f64> {
    pure(42.0)
}

This model always returns 42.0. The pure function creates a deterministic Model<f64>.

Key insight: Models are descriptions of computations, not the computations themselves.

Step 2: Adding Randomness

Now let's add some randomness by sampling from a probability distribution:

use fugue::*;

fn random_model() -> Model<f64> {
    sample(addr!("x"), Normal::new(0.0, 1.0).unwrap())
}

Step 3: Running Your Model

To actually get values from your model, you need to "run" it with a handler:

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn main() {
    let model = sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
    
    // Create a seeded random number generator
    let mut rng = StdRng::seed_from_u64(42);
    
    // Run the model with PriorHandler (forward sampling)
    let (value, trace) = runtime::handler::run(
        runtime::interpreters::PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        model,
    );
    
    println!("Sampled value: {:.4}", value);
    println!("Log probability: {:.4}", trace.total_log_weight());
}

Running this outputs:

Sampled value: 1.0175
Log probability: -0.9189

Step 4: Type Safety in Action

Fugue's type safety really shines with discrete distributions:

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn type_safe_examples() {
    let mut rng = StdRng::seed_from_u64(42);
    
    // Flip a coin - returns bool directly!
    let (is_heads, _) = runtime::handler::run(
        runtime::interpreters::PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        sample(addr!("coin"), Bernoulli::new(0.6).unwrap()),
    );
    
    // Natural boolean usage - no comparisons needed!
    let outcome = if is_heads { "Heads" } else { "Tails" };
    println!("Coin flip: {}", outcome);
    
    // Count events - returns u64 directly!
    let (count, _) = runtime::handler::run(
        runtime::interpreters::PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        sample(addr!("events"), Poisson::new(3.0).unwrap()),
    );
    
    println!("Event count: {} (no casting needed!)", count);
    
    // Choose category - returns usize for safe indexing!
    let colors = vec!["red", "green", "blue", "yellow"];
    let (idx, _) = runtime::handler::run(
        runtime::interpreters::PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        sample(addr!("color"), Categorical::uniform(4).unwrap()),
    );
    
    println!("Chosen color: {}", colors[idx]); // Safe indexing!
}

# Other PPLs - everything returns float
coin = sample("coin", Bernoulli(0.6))  # Returns 0.0 or 1.0 
if coin == 1.0:  # Need comparison ❌
    ...

count = sample("events", Poisson(3.0))  # Returns float
count_int = int(count)  # Need casting ❌

idx = sample("color", Categorical([...]))  # Returns float  
colors[int(idx)]  # Risky casting and indexing ❌

Step 5: Your First Bayesian Model

Now let's create a simple Bayesian model that learns from data:

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn estimate_mean(observation: f64) -> Model<f64> {
    sample(addr!("mu"), Normal::new(0.0, 2.0).unwrap())  // Prior belief
        .bind(move |mu| {
            observe(addr!("y"), Normal::new(mu, 1.0).unwrap(), observation)  // Likelihood
                .map(move |_| mu)  // Return the parameter
        })
}

fn main() {
    let observation = 3.0;  // We observed a value of 3.0
    let model = estimate_mean(observation);
    
    let mut rng = StdRng::seed_from_u64(42);
    let (estimated_mu, trace) = runtime::handler::run(
        runtime::interpreters::PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        model,
    );
    
    println!("Observation: {}", observation);
    println!("Estimated mean: {:.4}", estimated_mu);
    println!("Log probability: {:.4}", trace.total_log_weight());
}

Understanding Model Composition

Fugue models compose using two key operations:

map - Transform Values

let doubled = sample(addr!("x"), Normal::new(0.0, 1.0).unwrap())
    .map(|x| x * 2.0);  // Apply function to the result

bind - Dependent Computations

let dependent = sample(addr!("x"), Normal::new(0.0, 1.0).unwrap())
    .bind(|x| sample(addr!("y"), Normal::new(x, 0.5).unwrap()));  // y depends on x

Key Takeaways

After working through these examples, you should understand:

✅ Models are values: Model<T> represents a probabilistic computation
✅ Safe constructors: Distributions use .new().unwrap() for parameter validation
✅ Type safety: Distributions return natural types (bool, u64, f64)
✅ Addressing: addr!("name") gives names to random variables
✅ Execution: Models need handlers to run and produce values
✅ Composition: Use map and bind to build complex models from simple parts

What's Next?

You can now build and run basic probabilistic models! 🎉

Continue your journey:

Time: ~5 minutes • Next: Understanding Models

Understanding Models

Now that you can build basic models, let's understand the key concepts that make Fugue powerful. This will give you the mental framework to build sophisticated probabilistic programs.

The Big Picture: Models vs Execution

One of Fugue's key insights is separating model specification from execution:

graph LR
    subgraph "Your Code"
        M[Model Definition<br/>What to compute]
    end
    
    subgraph "Runtime System"  
        H[Handler<br/>How to execute]
        T[Trace<br/>What happened]
    end
    
    M --> H
    H --> T

Why this matters: The same model can be executed in different ways:

• Forward sampling (generate data from priors)
• Conditioning (inference given observations)
• Replay (MCMC proposals)
• Scoring (compute probabilities)

Addresses: The Key to Advanced Inference

Every sample and observe site needs a unique address:

use fugue::*;

// Good addressing ✅
let model = sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap())
    .bind(|mu| sample(addr!("sigma"), LogNormal::new(0.0, 0.5).unwrap())
        .bind(move |sigma| {
            observe(addr!("y1"), Normal::new(mu, sigma).unwrap(), 2.1);
            observe(addr!("y2"), Normal::new(mu, sigma).unwrap(), 1.9);
            pure((mu, sigma))
        }));

Addressing Patterns

Simple names for scalar parameters:

let mu = sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap());

Indexed addresses for collections:

for i in 0..10 {
    let x_i = sample(addr!("x", i), Normal::new(mu, 1.0).unwrap());
}

Scoped addresses for hierarchical models:

let encoder_z = sample(scoped_addr!("encoder", "z"), dist);
let decoder_z = sample(scoped_addr!("decoder", "z"), dist);

// ❌ DON'T: Random or non-deterministic addresses  
let addr = format!("param_{}", rng.gen::<u32>());  // NEVER!

// ❌ DON'T: Reuse addresses for different purposes
sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
sample(addr!("x"), Bernoulli::new(0.5).unwrap());  // Collision!

// ❌ DON'T: Missing addresses in loops
for i in 0..10 {
    sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());  // All same address!
}

Model Composition: The Monadic Structure

Fugue models follow a monadic pattern that makes complex models composable:

Three Fundamental Operations

graph TB
    subgraph "Model Building Blocks"
        P[pure: A → Model⟨A⟩<br/>Lift values into models]
        M[map: Model⟨A⟩ → ⟨A → B⟩ → Model⟨B⟩<br/>Transform outputs]  
        B[bind: Model⟨A⟩ → ⟨A → Model⟨B⟩⟩ → Model⟨B⟩<br/>Chain dependent computations]
    end

pure - Lift Values

Use pure to inject deterministic values into the probabilistic context:

use fugue::*;

// Lift a constant
let constant = pure(42.0);

// Lift computed values
let computed = pure(data.iter().sum::<f64>() / data.len() as f64);

map - Transform Outputs

Use map when you want to transform the result without adding randomness:

use fugue::*;

// Transform a single sample
let squared = sample(addr!("x"), Normal::new(0.0, 1.0).unwrap())
    .map(|x| x * x);

// Combine multiple values  
let sum = zip(
    sample(addr!("a"), Normal::new(0.0, 1.0).unwrap()),
    sample(addr!("b"), Normal::new(0.0, 1.0).unwrap())
).map(|(a, b)| a + b);

bind - Chain Dependent Computations

Use bind when the next random choice depends on a previous one:

use fugue::*;

// Dependent sampling: variance depends on mean
let hierarchical = sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap())
    .bind(|mu| sample(addr!("sigma"), LogNormal::new(mu.abs().ln(), 0.1).unwrap())
        .bind(move |sigma| sample(addr!("y"), Normal::new(mu, sigma).unwrap())));

// Conditional branching: choice affects distribution  
let mixture = sample(addr!("component"), Bernoulli::new(0.5).unwrap())
    .bind(|component| {
        if component {
            sample(addr!("value"), Normal::new(-2.0, 1.0).unwrap())
        } else {
            sample(addr!("value"), Normal::new(2.0, 1.0).unwrap())  
        }
    });

Advanced Composition Patterns

Building Collections

use fugue::*;

// Fixed-size collection
let samples = traverse_vec((0..10).collect(), |i| {
    sample(addr!("x", i), Normal::new(0.0, 1.0).unwrap())
});

// Data-driven collection
let observations = traverse_vec(data, |datum| {
    observe(addr!("y", datum.id), Normal::new(datum.mu, 1.0).unwrap(), datum.value)
});

Conditional Models

use fugue::*;

fn model_selection(data: &[f64]) -> Model<String> {
    sample(addr!("use_robust"), Bernoulli::new(0.2).unwrap())
        .bind(|use_robust| {
            if use_robust {
                // Robust model with t-distribution
                sample(addr!("df"), LogNormal::new(1.0, 0.5).unwrap())
                    .bind(|df| {
                        // Hypothetical t-distribution sampling
                        pure("robust".to_string())
                    })
            } else {
                // Standard normal model
                sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap())
                    .map(|_| "normal".to_string())
            }
        })
}

Hierarchical Structure

use fugue::*;

fn hierarchical_model(groups: Vec<Vec<f64>>) -> Model<Vec<f64>> {
    // Global hyperparameters
    sample(addr!("global_mu"), Normal::new(0.0, 2.0).unwrap())
        .bind(|global_mu| {
            sample(addr!("global_sigma"), LogNormal::new(0.0, 0.5).unwrap())
                .bind(move |global_sigma| {
                    // Group-level parameters  
                    traverse_vec(groups.into_iter().enumerate().collect(), move |(g, data)| {
                        sample(addr!("group_mu", g), Normal::new(global_mu, global_sigma).unwrap())
                            .bind(move |group_mu| {
                                // Individual observations
                                traverse_vec(data.into_iter().enumerate().collect(), move |(i, y)| {
                                    observe(addr!("y", g, i), Normal::new(group_mu, 1.0).unwrap(), y)
                                }).map(move |_| group_mu)
                            })
                    })
                })
        })
}

Mental Models for Success

Think in Terms of Generative Stories

Ask yourself: "How could this data have been generated?"

// Story: "Each person has a skill level, and their performance 
//         on each task reflects that skill plus task-specific noise"

let model = sample(addr!("skill"), Normal::new(100.0, 15.0).unwrap())  // Person's skill
    .bind(|skill| {
        traverse_vec(tasks, move |task| {
            let difficulty = task.difficulty;
            let expected_score = skill - difficulty;
            observe(addr!("score", task.id), Normal::new(expected_score, 5.0).unwrap(), task.actual_score)
        }).map(move |_| skill)
    });

Separate What from How

• What: Model describes relationships and distributions
• How: Handler determines execution strategy (sampling, inference, etc.)

// The SAME model...
let model = build_regression_model(&data);

// Can be executed different ways:
let (sample, _) = run(PriorHandler { /*...*/ }, model.clone());      // Forward sampling
let (_, scored_trace) = run(ScoreGivenTrace { /*...*/ }, model.clone()); // Compute likelihood  
let mcmc_chain = adaptive_mcmc_chain(rng, || model.clone(), 1000, 500); // MCMC inference

Key Takeaways

You now understand the foundational concepts:

✅ Separation of Concerns: Models describe computations, handlers execute them
✅ Addressing Strategy: Unique, stable addresses enable advanced inference
✅ Monadic Composition: pure, map, bind build complex models from simple parts
✅ Compositional Patterns: Collections, conditionals, and hierarchical structures
✅ Generative Thinking: Model the data generation process

What's Next?

You have the conceptual foundation to build sophisticated models! 🎉

Time: ~8 minutes • Next: Running Inference

Running Inference

You now know how to build probabilistic models. But models alone don't give you answers - you need inference to extract insights from them. Let's explore Fugue's inference algorithms!

What is Inference?

Inference is the process of learning about model parameters after seeing data. In Bayesian terms:

$$\text{Posterior}=\frac{\text{Prior}\times \text{Likelihood}}{\text{Evidence}}$$

graph LR
    subgraph "Before Data"
        P["Prior Beliefs<br/>p(theta)"]
    end

    subgraph "Observing Data"
        L["Likelihood<br/>p(y|theta)"]
        D["Data<br/>y₁, y₂, ..."]
    end

    subgraph "After Data"
        Post["Posterior Beliefs<br/>p(theta|y)"]
    end

    P --> Post
    L --> Post
    D --> Post

The Challenge

Most real models don't have analytical solutions. We need algorithms to approximate the posterior distribution.

Fugue's Inference Arsenal

1. MCMC (Markov Chain Monte Carlo) 🥇

Best for: Most general-purpose Bayesian inference

How it works: Generates samples that approximate the posterior distribution

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn coin_bias_model(heads: u64, total: u64) -> Model<f64> {
    sample(addr!("bias"), Beta::new(1.0, 1.0).unwrap())  // Prior
        .bind(move |bias| {
            observe(addr!("heads"), Binomial::new(total, bias).unwrap(), heads)  // Likelihood
                .map(move |_| bias)
        })
}

fn main() {
    let mut rng = StdRng::seed_from_u64(42);

    // Run adaptive MCMC
    let samples = inference::mh::adaptive_mcmc_chain(
        &mut rng,
        || coin_bias_model(7, 10),  // 7 heads out of 10 flips
        1000,  // number of samples
        500,   // warmup samples
    );

    // Extract bias estimates
    let bias_samples: Vec<f64> = samples.iter()
        .filter_map(|(_, trace)| trace.get_f64(&addr!("bias")))
        .collect();

    let mean_bias = bias_samples.iter().sum::<f64>() / bias_samples.len() as f64;
    println!("Estimated bias: {:.3}", mean_bias);

    // Compute 90% credible interval
    let mut sorted = bias_samples.clone();
    sorted.sort_by(|a, b| a.partial_cmp(b).unwrap());
    let lower = sorted[(0.05 * sorted.len() as f64) as usize];
    let upper = sorted[(0.95 * sorted.len() as f64) as usize];
    println!("90% credible interval: [{:.3}, {:.3}]", lower, upper);
}

When to use MCMC:

• ✅ Want exact posterior samples
• ✅ Moderate number of parameters (< 100)
• ✅ Can afford computation time
• ✅ Model evaluation is reasonably fast

2. SMC (Sequential Monte Carlo) 🎯

Best for: Sequential data and online learning

How it works: Uses particles to approximate the posterior, good for streaming data

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn main() {
    let mut rng = StdRng::seed_from_u64(42);

    // Generate particles from prior
    let particles = inference::smc::smc_prior_particles(
        &mut rng,
        1000,  // number of particles
        || coin_bias_model(7, 10),
    );

    println!("Generated {} particles", particles.len());

    // Compute weighted posterior mean
    let total_weight: f64 = particles.iter().map(|p| p.weight).sum();
    let weighted_mean: f64 = particles.iter()
        .filter_map(|p| {
            p.trace.get_f64(&addr!("bias"))
                .map(|bias| bias * p.weight)
        })
        .sum::<f64>() / total_weight;

    println!("Weighted posterior mean: {:.3}", weighted_mean);

    // Check effective sample size
    let weights: Vec<f64> = particles.iter().map(|p| p.weight).collect();
    let ess = 1.0 / weights.iter().map(|w| w * w).sum::<f64>();
    println!("Effective sample size: {:.1}", ess);
}

When to use SMC:

• ✅ Sequential/streaming data
• ✅ Online inference needed
• ✅ Many discrete latent variables
• ✅ Want to visualize inference process

3. Variational Inference (VI) ⚡

Best for: Fast approximate inference with many parameters

How it works: Finds the best approximation within a family of simple distributions

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn main() {
    let mut rng = StdRng::seed_from_u64(42);

    // Estimate ELBO (Evidence Lower BOund)
    let elbo = inference::vi::estimate_elbo(
        &mut rng,
        || coin_bias_model(7, 10),
        100,  // number of samples for estimation
    );

    println!("ELBO estimate: {:.3}", elbo);

    // For more sophisticated VI, you'd set up a variational guide
    // and optimize it (see the VI tutorial for details)
}

When to use VI:

• ✅ Need fast approximate inference
• ✅ Many parameters (> 100)
• ✅ Can accept approximation error
• ✅ Want predictable runtime

4. ABC (Approximate Bayesian Computation) 🎲

Best for: Models where likelihood is intractable or expensive

How it works: Simulation-based inference using distance between simulated and observed data

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn main() {
    let mut rng = StdRng::seed_from_u64(42);

    // ABC with summary statistics
    let observed_summary = 0.7; // 7/10 = 0.7 success rate
    let samples = inference::abc::abc_scalar_summary(
        &mut rng,
        || sample(addr!("bias"), Beta::new(1.0, 1.0).unwrap()), // Prior only
        |trace| trace.get_f64(&addr!("bias")).unwrap_or(0.0), // Extract bias
        observed_summary,  // Target summary statistic
        0.1,              // Tolerance
        1000,             // Max samples to try
    );

    println!("ABC accepted {} samples", samples.len());

    if !samples.is_empty() {
        let abc_estimates: Vec<f64> = samples.iter()
            .filter_map(|trace| trace.get_f64(&addr!("bias")))
            .collect();
        let abc_mean = abc_estimates.iter().sum::<f64>() / abc_estimates.len() as f64;
        println!("ABC estimated bias: {:.3}", abc_mean);
    }
}

When to use ABC:

• ✅ Likelihood is intractable or very expensive
• ✅ Can simulate from the model easily
• ✅ Have good summary statistics
• ✅ Can tolerate approximation error

Algorithm Comparison

Practical Inference Workflow

Here's a typical workflow for real inference:

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn inference_workflow() {
    let mut rng = StdRng::seed_from_u64(42);

    // 1. Define your model
    let model = || coin_bias_model(17, 25);  // 17 heads out of 25 flips

    // 2. Run inference (adaptive MCMC is often a good default)
    let samples = inference::mh::adaptive_mcmc_chain(
        &mut rng,
        model,
        2000,  // samples
        1000,  // warmup
    );

    // 3. Extract parameter values
    let bias_samples: Vec<f64> = samples.iter()
        .filter_map(|(_, trace)| trace.get_f64(&addr!("bias")))
        .collect();

    // 4. Compute summary statistics
    let mean = bias_samples.iter().sum::<f64>() / bias_samples.len() as f64;
    let variance = bias_samples.iter()
        .map(|&x| (x - mean).powi(2))
        .sum::<f64>() / (bias_samples.len() - 1) as f64;
    let std_dev = variance.sqrt();

    println!("Posterior Summary:");
    println!("  Mean: {:.3}", mean);
    println!("  Std Dev: {:.3}", std_dev);

    // 5. Check convergence (effective sample size)
    let ess = inference::diagnostics::effective_sample_size(&bias_samples);
    println!("  Effective Sample Size: {:.1}", ess);

    if ess > 100.0 {
        println!("  ✅ Good mixing!");
    } else {
        println!("  ⚠️ Poor mixing - consider more samples");
    }

    // 6. Make predictions
    println!("\nPredictions:");
    println!("  P(bias > 0.5) = {:.2}",
        bias_samples.iter().filter(|&&b| b > 0.5).count() as f64 / bias_samples.len() as f64);
}

Choosing the Right Algorithm

Decision Tree

graph TD
    A[Need inference?] -->|Yes| B[Real-time/online?]
    A -->|No| Z[Use PriorHandler<br/>for forward sampling]

    B -->|Yes| SMC[SMC]
    B -->|No| C[Likelihood tractable?]

    C -->|No| ABC[ABC]
    C -->|Yes| D[Many parameters?]

    D -->|Yes > 100| VI[Variational Inference]
    D -->|No < 100| E[Need exact samples?]

    E -->|Yes| MCMC[MCMC]
    E -->|No| VI2[VI for speed]

Rules of Thumb

1. Start with MCMC for most problems - it's the most general
2. Use SMC if you have sequential/streaming data
3. Use VI if you need speed and can accept approximation
4. Use ABC only when likelihood is truly intractable

Key Takeaways

You now know how to extract insights from your models:

✅ Inference Purpose: Learn parameters from data using Bayesian updating
✅ Algorithm Options: MCMC, SMC, VI, ABC each have their strengths
✅ Practical Workflow: Define model → Run inference → Extract parameters → Check diagnostics
✅ Algorithm Selection: Choose based on problem characteristics and requirements

What's Next?

You've completed Getting Started! 🎉

Time: ~5 minutes • Next: Complete Tutorials

How-To Guides

The How-To guides provide practical, task-oriented instructions for accomplishing specific goals with Fugue. Unlike tutorials that teach concepts step-by-step, these guides assume you understand the basics and want to solve particular problems efficiently.

Guide Overview

These guides are designed to be example-first and immediately actionable. Each guide includes comprehensive, executable code examples that serve as the canonical source of truth for the patterns they demonstrate.

📊 Working with Distributions

When to use: You need to understand Fugue's type-safe distribution system, parameter validation, or probability calculations.

What you'll learn:

• Type-safe distribution usage (bool, u64, usize, f64 return types)
• Parameter validation and error handling
• Continuous vs discrete distribution patterns
• Categorical distributions and safe indexing
• Distribution composition and practical modeling
• Probability calculations and testing strategies

Key patterns: Natural return types, parameter validation, distribution testing

🏗️ Building Complex Models

When to use: You want to compose sophisticated probabilistic models using Fugue's macro system and advanced patterns.

What you'll learn:

• prob! macro for do-notation style probabilistic programming
• plate! macro for vectorized operations and array processing
• scoped_addr! macro for hierarchical address management
• Model composition and sequential dependencies
• Hierarchical modeling patterns
• Bayesian linear regression and mixture models

Key patterns: Monadic composition, vectorization, hierarchical structure

⚡ Optimizing Performance

When to use: Your models need to run efficiently in production or handle large-scale inference workloads.

What you'll learn:

• Memory pooling with TracePool and PooledPriorHandler
• Numerical stability with log-space computations
• Efficient trace construction with TraceBuilder
• Copy-on-write traces for MCMC optimization
• Batch processing patterns
• Performance monitoring and measurement

Key patterns: Memory optimization, numerical stability, batch processing

🔍 Debugging Models

When to use: Your probabilistic models aren't behaving as expected, or you need to diagnose issues in inference.

What you'll learn:

• Comprehensive trace inspection and analysis
• Type-safe value access with proper error handling
• Model validation against analytical solutions
• Safe vs strict handler usage for error resilience
• MCMC diagnostics and convergence assessment
• Performance and memory debugging techniques

Key patterns: Trace analysis, validation testing, diagnostic metrics

🎛️ Custom Handlers

When to use: You need to extend Fugue's execution model with custom behavior, logging, or specialized inference algorithms.

What you'll learn:

• Complete Handler trait implementation
• Decorator pattern for cross-cutting concerns
• Stateful handlers for analytics and monitoring
• Conditional filtering and value modification
• Performance monitoring integration
• Custom inference algorithms (MCMC-like patterns)
• Handler composition and chaining

Key patterns: Algebraic effects, decorator composition, custom inference

🚀 Production Deployment

When to use: You're deploying probabilistic models to production environments and need reliability, monitoring, and operational excellence.

What you'll learn:

• Error handling and graceful degradation patterns
• Circuit breaker implementation for fault tolerance
• Configuration management for multiple environments
• Comprehensive metrics and Prometheus integration
• Automated health checks and system validation
• Input validation and security best practices
• Deployment strategies (blue-green, canary, rolling)

Key patterns: Fault tolerance, observability, operational readiness

How to Use These Guides

🎯 Task-Oriented Approach

Each guide focuses on solving specific problems:

• Need to understand distributions? → Start with "Working with Distributions"
• Building complex models? → "Building Complex Models" has the macros and patterns
• Performance issues? → "Optimizing Performance" covers memory and numerical techniques
• Models not working? → "Debugging Models" provides diagnostic approaches
• Need custom behavior? → "Custom Handlers" shows how to extend the system
• Going to production? → "Production Deployment" covers operational concerns

📚 Progressive Complexity

The guides are ordered by increasing complexity:

graph TD
    A["📊 Working with<br/>Distributions"] --> B["🏗️ Building Complex<br/>Models"]
    B --> C["⚡ Optimizing<br/>Performance"]
    A --> D["🚀 Production<br/>Deployment"]
    E["🔍 Debugging<br/>Models"] --> F["🎛️ Custom<br/>Handlers"]
    F --> D
    C --> D

• Start with distributions and model building
• Add performance optimization when needed
• Use debugging when things go wrong
• Extend with custom handlers for specialized needs
• Deploy with production patterns for real applications

🔄 Cross-References

Guides frequently reference each other:

• Performance optimization builds on complex models
• Debugging techniques apply to all model types
• Custom handlers can incorporate performance patterns
• Production deployment uses patterns from all previous guides

📝 Example-First Philosophy

Every guide follows the same structure:

1. Executable examples as the source of truth
2. Comprehensive code snippets with anchor tags
3. Practical explanations of when and how to use patterns
4. Testing strategies to verify correctness
5. Best practices learned from real-world usage

Code Examples

All code examples in these guides are:

• ✅ Executable: Run with cargo run --example <guide_name>
• ✅ Tested: Verified with cargo test --examples
• ✅ Documented: Included via {{#include}} from example files
• ✅ Comprehensive: Cover real-world usage patterns
• ✅ Type-Safe: Leverage Rust's type system throughout

Quick Reference

Integration with Other Documentation

🚀 Getting Started → How-To Guides

After completing the Getting Started tutorials, use these guides to solve specific problems in your own projects.

📖 How-To Guides → Tutorials

For deeper conceptual understanding, see the Tutorials section which provides comprehensive examples of complete applications.

🔧 How-To Guides → API Reference

For detailed API documentation, consult the API Reference section for specific functions and types mentioned in these guides.

Contributing

When adding new How-To guides:

1. Create executable examples first in examples/
2. Use anchor tags to mark code sections for inclusion
3. Write the guide using {{#include}} for all code snippets
4. Test thoroughly with both cargo test and mdbook test
5. Update this README with the new guide information

Each guide should solve specific, practical problems that users commonly encounter when working with Fugue in real applications.

Working with Distributions

Fugue's type-safe distribution system represents a principled approach to probabilistic programming, eliminating entire classes of runtime errors through rigorous type theory while preserving the full expressiveness of statistical modeling. This guide demonstrates the mathematical foundations and practical applications of Fugue's distribution architecture.

Type Safety in Practice

Traditional probabilistic programming libraries return f64 for everything, leading to casting overhead and runtime errors. Fugue distributions return their natural types:

    // Demonstrate natural return types
    let coin = Bernoulli::new(0.5).unwrap();
    let flip: bool = coin.sample(&mut rng); // Natural boolean

    if flip {
        println!("🪙 Coin flip result: Heads!");
    } else {
        println!("🪙 Coin flip result: Tails!");
    }

No casting, no comparisons with floating-point values—just natural boolean logic.

Continuous Distributions

Continuous distributions in Fugue model phenomena over uncountable domains $\mathcal{S}\subseteq {\mathbb{R}}^n$. The probability density function $f_X(x)$ satisfies the normalization condition:

$${\int }_{\mathcal{S}}{f}_X(x)dx=1$$

For computational stability, Fugue operates in log-space by default, computing $\log f_X(x)$ to avoid numerical underflow:

    // Working with continuous distributions
    let standard_normal = Normal::new(0.0, 1.0).unwrap();
    let sample: f64 = standard_normal.sample(&mut rng);
    println!("📊 Standard normal sample: {:.3}", sample);

    // Compute log-probability density
    let log_density = standard_normal.log_prob(&0.0); // Peak of standard normal
    println!(
        "📈 Log-density at x=0: {:.3} (peak of standard normal)",
        log_density
    );

    // Custom parameters
    let measurement_model = Normal::new(10.0, 0.5).unwrap();
    let measurement = measurement_model.sample(&mut rng);
    println!("🔬 Sensor measurement (μ=10.0, σ=0.5): {:.3}", measurement);

Key Points:

• sample() returns f64 for direct arithmetic
• log_prob() computes log-density (avoids numerical underflow)
• Parameter validation happens at construction time

Discrete Distributions

Discrete distributions operate over countable support sets $\mathcal{S}\subseteq {\mathbb{Z}}^n$ or finite sets. The probability mass function satisfies:

$$\sum _{x\in \mathcal{S}}P(X=x)=1$$

Fugue enforces this constraint at construction time and leverages natural integer types to eliminate precision loss from floating-point representation:

    // Working with discrete distributions

    // Count data
    let event_rate = Poisson::new(3.0).unwrap();
    let count: u64 = event_rate.sample(&mut rng);
    println!("📅 Event count (λ=3.0): {} events", count);

    // Log-probability mass
    let prob_3_events = event_rate.log_prob(&3);
    println!(
        "🎯 Log-probability of exactly 3 events: {:.3}",
        prob_3_events
    );

    // Use counts directly in calculations
    let total_cost = count * 50; // Direct arithmetic with u64
    println!("💰 Total cost ({} events × $50): ${}", count, total_cost);

Benefits:

• u64 counts support direct arithmetic without casting
• No precision loss from floating-point representations
• Natural integration with Rust's type system

Safe Categorical Sampling

Categorical distributions return usize for safe array indexing:

    // Safe categorical sampling
    let choices = vec![0.3, 0.5, 0.2]; // Three categories
    let categorical = Categorical::new(choices).unwrap();
    let selected: usize = categorical.sample(&mut rng);

    // Safe array indexing (no bounds checking needed)
    let options = ["Option A", "Option B", "Option C"];
    println!(
        "🎲 Categorical choice (weights: 0.3, 0.5, 0.2): {}",
        options[selected]
    );

    // Uniform categorical
    let uniform_choice = Categorical::uniform(5).unwrap();
    let idx: usize = uniform_choice.sample(&mut rng);
    println!("🎯 Uniform random index (0-4): {}", idx);

Parameter Validation

Fugue enforces mathematical constraints through compile-time and runtime validation. Each distribution family ${\mathcal{D}}_{\theta }$ has a parameter space $\Theta$ defining valid configurations:

graph TD
    A[Parameter Input θ] --> B{θ ∈ Θ?}
    B -->|Yes| C[Distribution Construction]
    B -->|No| D[ValidationError]
    C --> E[Type-Safe Sampling]
    D --> F[Early Failure Detection]

Constraint Examples:

• Normal Distribution: $\mathcal{N}(\mu ,{\sigma }^2)$ requires $\sigma >0$
• Beta Distribution: $\text{Beta}(\alpha ,\beta )$ requires $\alpha ,\beta >0$
• Categorical Distribution: ${\sum }_i{p}_i=1$ and $p_i\ge 0\forall i$

    // Distribution parameter validation

    // This will return an error
    match Normal::new(0.0, -1.0) {
        Ok(_) => println!("✅ Normal(μ=0.0, σ=-1.0) created successfully"),
        Err(e) => println!("❌ Normal(μ=0.0, σ=-1.0) failed: {:?}", e),
    }

    // Beta distribution parameters must be positive
    match Beta::new(0.0, 1.0) {
        Ok(_) => println!("✅ Beta(α=0.0, β=1.0) created successfully"),
        Err(e) => println!("❌ Beta(α=0.0, β=1.0) failed: {:?}", e),
    }

    // Poisson rate must be non-negative
    match Poisson::new(-1.0) {
        Ok(_) => println!("✅ Poisson(λ=-1.0) created successfully"),
        Err(e) => println!("❌ Poisson(λ=-1.0) failed: {:?}", e),
    }

Validation Rules:

• Normal: σ > 0
• Beta: α > 0, β > 0
• Poisson: λ ≥ 0
• Categorical: probabilities sum to 1, all non-negative

Storing Mixed Distributions

Use trait objects for collections of distributions with the same return type:

    // Storing different distributions together
    let continuous_dists: Vec<Box<dyn Distribution<f64>>> = vec![
        Normal::new(0.0, 1.0).unwrap().clone_box(),
        Beta::new(2.0, 5.0).unwrap().clone_box(),
        Uniform::new(-1.0, 1.0).unwrap().clone_box(),
    ];

    // Sample from each
    for (i, dist) in continuous_dists.iter().enumerate() {
        let sample = dist.sample(&mut rng);
        let dist_name = match i {
            0 => "Normal(0,1)",
            1 => "Beta(2,5)",
            2 => "Uniform(-1,1)",
            _ => "Unknown",
        };
        println!("📦 {} sample: {:.3}", dist_name, sample);
    }

This enables dynamic distribution selection and model composition patterns.

Practical Modeling Patterns

Common modeling scenarios demonstrate natural type usage:

    // Practical modeling examples

    // Model a sensor with noise
    let true_temperature = 20.5; // True value
    let sensor_noise = Normal::new(0.0, 0.2).unwrap(); // Measurement error
    let measured_temp = true_temperature + sensor_noise.sample(&mut rng);
    println!(
        "🌡️  True temperature: {:.2}°C → Measured: {:.2}°C",
        true_temperature, measured_temp
    );

    // Count model for arrivals
    let arrival_rate = Poisson::new(2.5).unwrap(); // 2.5 arrivals per hour
    let hourly_arrivals = arrival_rate.sample(&mut rng);
    println!(
        "🚪 Expected arrivals: 2.5/hour → Actual: {} arrivals",
        hourly_arrivals
    );

    // Decision model
    let decision_prob = 0.7;
    let decision = Bernoulli::new(decision_prob).unwrap();
    let will_buy = decision.sample(&mut rng);
    if will_buy {
        println!("🛒 Customer decision (p=0.7): Will make a purchase");
    } else {
        println!("🚶 Customer decision (p=0.7): Will not purchase");
    }

Each distribution serves its natural domain without artificial conversions.

Working with Log-Probabilities

Logarithmic probability computation is essential for numerical stability in probabilistic programming. Consider the log-sum-exp operation for computing:

$$\log \left(\sum _{i=1}^n{e}^{x_i}\right)=x_{\max }+\log \left(\sum _{i=1}^n{e}^{x_i-x_{\max }}\right)$$

where $x_{\max }={\max }_i{x}_i$. This formulation prevents overflow when $|x_i|$ is large:

    // Working with log-probabilities

    let normal = Normal::new(100.0, 15.0).unwrap();

    // Multiple observations
    let observations = vec![98.5, 102.1, 99.8, 101.5, 97.2];
    let mut total_log_prob = 0.0;

    println!(
        "📋 Evaluating {} observations under Normal(μ=100.0, σ=15.0):",
        observations.len()
    );
    for obs in &observations {
        let log_p = normal.log_prob(obs);
        total_log_prob += log_p;
        println!("   x={}: log P(x) = {:.3}", obs, log_p);
    }

    println!("🔢 Joint log-probability: {:.3}", total_log_prob);

    // Convert back to probability (be careful with underflow!)
    if total_log_prob > -700.0 {
        // Avoid underflow
        let probability = total_log_prob.exp();
        println!("📊 Joint probability: {:.2e}", probability);
    } else {
        println!("⚠️  Joint probability too small to represent as f64");
    }

Advanced Patterns

For complex modeling scenarios, see these patterns:

Hierarchical Models

    // Hierarchical prior structure
    let global_mean = Normal::new(0.0, 10.0).unwrap();
    let mu = global_mean.sample(&mut rng);

    let group_precision = Gamma::new(2.0, 0.5).unwrap();
    let tau = group_precision.sample(&mut rng);
    let sigma = (1.0 / tau).sqrt(); // Convert precision to std dev

    // Individual observations from hierarchical model
    let individual = Normal::new(mu, sigma).unwrap();
    let observation = individual.sample(&mut rng);

    println!("🌐 Global mean: {:.3}", mu);
    println!("📊 Group std dev: {:.3}", sigma);
    println!("👤 Individual observation: {:.3}", observation);

Mixture Components

    // Mixture model components
    let mixture_weights = vec![0.6, 0.3, 0.1];
    let component_selector = Categorical::new(mixture_weights).unwrap();
    let selected_component: usize = component_selector.sample(&mut rng);

    // Different components
    let components = [
        Normal::new(-2.0, 0.5).unwrap(),
        Normal::new(0.0, 1.0).unwrap(),
        Normal::new(3.0, 0.8).unwrap(),
    ];

    let sample = components[selected_component].sample(&mut rng);
    println!(
        "🎯 Selected component {}: sample = {:.3}",
        selected_component, sample
    );

Conjugate Priors

    // Beta-Bernoulli conjugacy
    let prior_alpha = 2.0;
    let prior_beta = 8.0;
    let prior = Beta::new(prior_alpha, prior_beta).unwrap();
    let p: f64 = prior.sample(&mut rng);

    // Simulate some trials
    let trials = 20;
    let mut successes = 0;
    let bernoulli = Bernoulli::new(p).unwrap();

    for _ in 0..trials {
        if bernoulli.sample(&mut rng) {
            successes += 1;
        }
    }

    // Posterior parameters (conjugate update)
    let posterior_alpha = prior_alpha + successes as f64;
    let posterior_beta = prior_beta + (trials - successes) as f64;
    let posterior = Beta::new(posterior_alpha, posterior_beta).unwrap();
    let updated_p = posterior.sample(&mut rng);

    println!("🎲 Prior p: {:.3}", p);
    println!("📈 Observed: {}/{} successes", successes, trials);
    println!("🔄 Posterior p: {:.3}", updated_p);

Testing Your Distributions

Always test distribution properties and parameter validation:

    #[test]
    fn test_distribution_properties() {
        let mut rng = thread_rng();

        // Test type safety
        let coin = Bernoulli::new(0.5).unwrap();
        let flip: bool = coin.sample(&mut rng);
        assert!(flip == true || flip == false); // Must be boolean

        // Test parameter validation
        assert!(Normal::new(0.0, -1.0).is_err());
        assert!(Beta::new(0.0, 1.0).is_err());
        assert!(Poisson::new(-1.0).is_err());

        // Test valid distributions
        let normal = Normal::new(0.0, 1.0).unwrap();
        let sample = normal.sample(&mut rng);
        let log_prob = normal.log_prob(&sample);
        assert!(log_prob.is_finite());
    }

Common Pitfalls

1. Underflow in probability space: Use log-probabilities for accumulation
2. Parameter validation: Check constructor errors, don't assume success
3. Precision with counts: Use u64 return types directly, avoid f64 conversion
4. Categorical indexing: Trust the usize return—it's guaranteed valid

Next Steps

• Complex Models: See Building Complex Models for compositional patterns
• Debugging: Check out Debugging Models for troubleshooting
• Custom Logic: Learn Custom Handlers for specialized inference

The type-safe distribution system eliminates entire classes of runtime errors while making statistical code more readable and maintainable.

Building Complex Models

Fugue's compositional architecture is grounded in category theory and monadic structures, enabling the systematic construction of sophisticated probabilistic models through principled composition operators. This guide explores the mathematical foundations and practical applications of Fugue's macro system for building complex probabilistic programs.

Do-Notation with prob!

The prob! macro implements monadic do-notation for probabilistic computations, providing a natural syntax for sequential dependence. Formally, it translates:

$$\begin{array}{rl}\text{prob!}\{& \\ & x\leftarrow {\mathcal{M}}_1\\ & y\leftarrow {\mathcal{M}}_2(x)\\ & \text{pure}(f(x,y))\end{array}\}$$

into the monadic composition ${\mathcal{M}}_1\gg =\lambda x.{\mathcal{M}}_2(x)\gg =\lambda y.\eta (f(x,y))$:

graph LR
    A[M₁] -->|bind| B[λx.M₂⁽ˣ⁾]
    B -->|bind| C[λy.η⁽f⁽ˣ'ʸ⁾⁾]
    C --> D[Result⁽ᶻ⁾]

    // Simple do-notation style probabilistic program
    let _simple_model = prob!(
        let x <- sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
        let y <- sample(addr!("y"), Normal::new(x, 0.5).unwrap());
        let sum = x + y;  // Regular variable assignment
        pure(sum)
    );
    println!("✅ Created simple model with prob! macro");

Key Features:

• <- for probabilistic binding (monadic bind)
• = for regular variable assignment
• pure() to lift deterministic values
• Natural control flow without callback nesting

Vectorized Operations with plate!

The plate! macro implements plate notation from graphical models, representing conditionally independent replications. Given $N$ independent observations, plate notation expresses:

$$P(\mathbf{x}\mid \theta )=\prod _{i=1}^{N}P(x_i\mid \theta )$$

The computational graph shows the independence structure:

graph TB
    subgraph "Plate: i ∈ {1..N}"
        A[θ] --> B1[x₁]
        A --> B2[x₂]
        A --> B3[...]
        A --> BN[xₙ]
    end

    // Independent samples using plate notation
    let _vector_model = plate!(i in 0..5 => {
        sample(addr!("sample", i), Normal::new(0.0, 1.0).unwrap())
    });
    println!("✅ Created vectorized model with {} samples", 5);

    // Plate with observations
    let observations = [1.2, -0.5, 2.1, 0.8, -1.0];
    let n_obs = observations.len();
    let _observed_model = plate!(i in 0..n_obs => {
        observe(addr!("obs", i), Normal::new(0.0, 1.0).unwrap(), observations[i])
    });
    println!("✅ Created observation model for {} data points", n_obs);

Benefits:

• Automatic address indexing prevents conflicts
• Natural iteration over data structures
• Vectorized likelihood computations
• Clear intent for independent operations

Hierarchical Address Management

Complex models require systematic parameter organization following a tree-structured address space. The address hierarchy $\mathcal{A}$ forms a prefix tree where each node represents a scope:

$$\mathcal{A}=\{a_1,a_1.a_2,a_1.a_2.a_3,\dots \}$$

This hierarchical structure prevents address collisions and enables efficient parameter lookup:

    // Hierarchical model using scoped addresses
    let _hierarchical_model = prob!(
        let global_mu <- sample(addr!("global_mu"), Normal::new(0.0, 10.0).unwrap());
        let group_mu <- sample(scoped_addr!("group", "mu", "{}", 0),
                              Normal::new(global_mu, 1.0).unwrap());
        pure((global_mu, group_mu))
    );
    println!("✅ Created hierarchical model with scoped addresses");

Address Strategy:

• scoped_addr! prevents parameter name collisions
• Hierarchical structure mirrors model dependencies
• Systematic naming aids debugging and introspection
• Indices enable parameter arrays

Sequential Dependencies

Sequential models exhibit temporal dependence where the state at time $t$ depends on previous states. This creates a Markov chain structure:

$$P(x_{1:T})=P(x_1)\prod _{t=2}^{T}P(x_t\mid x_{t-1})$$

The computational challenge lies in maintaining state consistency while enabling efficient inference:

    // Sequential model with dependencies
    let _sequential_model = prob! {
        let states <- plate!(t in 0..3 => {
            sample(addr!("x", t), Normal::new(0.0, 1.0).unwrap())
                .bind(move |x_t| {
                    observe(addr!("y", t), Normal::new(x_t, 0.5).unwrap(), 1.0 + t as f64)
                        .map(move |_| x_t)
                })
        });

        pure(states)
    };
    println!("✅ Created sequential model with observations");

Patterns:

• Explicit state threading through computations
• Observation conditioning at each time step
• Autoregressive dependencies
• Mixed probabilistic and deterministic updates

Composable Model Functions

Build reusable model components:

    // Helper function to create a component model
    fn create_normal_component(name: &str, mean: f64, std: f64) -> Model<f64> {
        sample(addr!(name), Normal::new(mean, std).unwrap())
    }

    // Compose multiple components
    let _composition_model = prob! {
        let param1 <- create_normal_component("param1", 0.0, 1.0);
        let param2 <- create_normal_component("param2", 2.0, 0.5);
        let combined = param1 * param2;
        pure(combined)
    };
    println!("✅ Created composed model with reusable components");

Design Principles:

• Functions return Model<T> for composability
• Pattern matching enables model selection
• Pure functions for deterministic transformations
• Higher-order functions for model templates

Advanced Address Patterns

For large-scale models like neural networks:

    // Complex addressing for large models
    let _neural_layer_model = plate!(layer in 0..3 => {
        let layer_size = match layer {
            0 => 4,
            1 => 8,
            2 => 1,
            _ => 1,
        };

        plate!(i in 0..layer_size => {
            sample(
                scoped_addr!("layer", "weight", "{}_{}", layer, i),
                Normal::new(0.0, 0.1).unwrap()
            )
        })
    });
    println!("✅ Created neural network parameter structure");

Scaling Strategies:

• Systematic parameter naming conventions
• Multi-level scoping for complex architectures
• Consistent indexing schemes
• Hierarchical parameter organization

Mixing Styles for Flexibility

Combine macros with traditional function composition:

    // Mixture model with component selection
    let _mixture_model = prob! {
        let component <- sample(addr!("component"), Bernoulli::new(0.3).unwrap());
        let mu = if component { -2.0 } else { 2.0 };
        let x <- sample(addr!("x"), Normal::new(mu, 1.0).unwrap());
        pure((component, x))
    };
    println!("✅ Created mixture model with 2 components");

Best Practices:

• Use functions for reusable components
• Use macros for readable composition
• Separate concerns (priors, likelihood, observations)
• Document parameter dependencies

Real-World Applications

Bayesian Linear Regression

Bayesian linear regression models the relationship $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ with uncertainty quantification:

$$\begin{array}{rl}\mathbf{\beta }& \sim \mathcal{N}({\mathbf{\mu }}_0,{\mathbf{\Sigma }}_0)\\ {\sigma }^2& \sim \text{InverseGamma}(\alpha ,\beta )\\ y_i\mid {\mathbf{x}}_i,\mathbf{\beta },{\sigma }^2& \sim \mathcal{N}({\mathbf{x}}_i^T\mathbf{\beta },{\sigma }^2)\end{array}$$

    // Complete Bayesian linear regression
    let x_data = [1.0, 2.0, 3.0, 4.0, 5.0];
    let y_data = [2.1, 3.9, 6.2, 8.1, 9.8];
    let n = x_data.len();

    let _regression_model = prob! {
        let intercept <- sample(addr!("intercept"), Normal::new(0.0, 10.0).unwrap());
        let slope <- sample(addr!("slope"), Normal::new(0.0, 10.0).unwrap());
        let precision <- sample(addr!("precision"), Gamma::new(1.0, 1.0).unwrap());
        let sigma = (1.0 / precision).sqrt();

        let _likelihood <- plate!(i in 0..n => {
            let predicted = intercept + slope * x_data[i];
            observe(addr!("y", i), Normal::new(predicted, sigma).unwrap(), y_data[i])
        });

        pure((intercept, slope, sigma))
    };
    println!("✅ Created Bayesian linear regression model");

Hierarchical Clustering

Hierarchical models implement partial pooling through multi-level parameter structures. The hierarchy enables information sharing across groups while maintaining group-specific effects:

$$\begin{array}{rl}{\mu }_{\text{pop}}& \sim \mathcal{N}(0,{\tau }_{\text{pop}}^2)\\ {\sigma }_{\text{group}}& \sim \text{HalfCauchy}({\sigma }_0)\\ {\mu }_j\mid {\mu }_{\text{pop}},{\sigma }_{\text{group}}& \sim \mathcal{N}({\mu }_{\text{pop}},{\sigma }_{\text{group}}^2)\\ y_{ij}\mid {\mu }_j,{\sigma }_j& \sim \mathcal{N}({\mu }_j,{\sigma }_j^2)\end{array}$$

    // Simplified hierarchy to avoid nested macro issues
    let _multilevel_model = prob!(
        let pop_mean <- sample(addr!("pop_mean"), Normal::new(0.0, 10.0).unwrap());
        let _pop_precision <- sample(addr!("pop_precision"), Gamma::new(2.0, 0.5).unwrap());
        let group_mean <- sample(scoped_addr!("group", "mean", "{}", 0),
                                Normal::new(pop_mean, 1.0).unwrap());
        pure((pop_mean, group_mean))
    );
    println!("✅ Created hierarchical model structure");

State Space Models

Sequential latent variable models:

    // Sequential model with dependencies
    let _sequential_model = prob! {
        let states <- plate!(t in 0..3 => {
            sample(addr!("x", t), Normal::new(0.0, 1.0).unwrap())
                .bind(move |x_t| {
                    observe(addr!("y", t), Normal::new(x_t, 0.5).unwrap(), 1.0 + t as f64)
                        .map(move |_| x_t)
                })
        });

        pure(states)
    };
    println!("✅ Created sequential model with observations");

Multi-Level Hierarchies

Population → Groups → Individuals structure:

    // Simplified hierarchy to avoid nested macro issues
    let _multilevel_model = prob!(
        let pop_mean <- sample(addr!("pop_mean"), Normal::new(0.0, 10.0).unwrap());
        let _pop_precision <- sample(addr!("pop_precision"), Gamma::new(2.0, 0.5).unwrap());
        let group_mean <- sample(scoped_addr!("group", "mean", "{}", 0),
                                Normal::new(pop_mean, 1.0).unwrap());
        pure((pop_mean, group_mean))
    );
    println!("✅ Created hierarchical model structure");

Key Features:

• Partial pooling across hierarchy levels
• Systematic parameter organization
• Natural shrinkage properties
• Scalable to large group structures

Configurable Model Factories

Dynamic model construction:

    // Helper function to create a component model
    fn create_normal_component(name: &str, mean: f64, std: f64) -> Model<f64> {
        sample(addr!(name), Normal::new(mean, std).unwrap())
    }

    // Compose multiple components
    let _composition_model = prob! {
        let param1 <- create_normal_component("param1", 0.0, 1.0);
        let param2 <- create_normal_component("param2", 2.0, 0.5);
        let combined = param1 * param2;
        pure(combined)
    };
    println!("✅ Created composed model with reusable components");

Flexibility Benefits:

• Runtime model configuration
• Conditional model components
• A/B testing different model structures
• Experiment management

Testing Complex Models

Model validation requires systematic testing across multiple dimensions: syntactic correctness, semantic validity, and statistical consistency:

graph TD
    A[Model M] --> B[Syntactic Tests]
    A --> C[Semantic Tests]
    A --> D[Statistical Tests]

    B --> E[Type Checking]
    B --> F[Address Uniqueness]

    C --> G[Trace Validity]
    C --> H[Parameter Bounds]

    D --> I[Prior Predictive]
    D --> J[Posterior Consistency]

    E --> K{All Pass?}
    F --> K
    G --> K
    H --> K
    I --> K
    J --> K

    K -->|Yes| L[Model Validated]
    K -->|No| M[Refinement Required]

Testing Hierarchy:

1. Unit Tests: Individual model components
2. Integration Tests: Model composition correctness
3. Statistical Tests: Distributional properties
4. Performance Tests: Scalability and efficiency

    #[test]
    fn test_model_composition() {
        // Test that models construct without errors
        let _simple = prob! {
            let x <- sample(addr!("test_x"), Normal::new(0.0, 1.0).unwrap());
            pure(x)
        };

        // Test plate notation
        let _plate_model = plate!(i in 0..3 => {
            sample(addr!("plate_test", i), Normal::new(0.0, 1.0).unwrap())
        });

        // Test scoped addresses
        let addr1 = scoped_addr!("test", "param");
        let addr2 = scoped_addr!("test", "param", "{}", 42);

        // Addresses should be different
        assert_ne!(addr1.0, addr2.0);
        assert!(addr2.0.contains("42"));

        // Test hierarchical model construction
        let _hierarchical = prob! {
            let global <- sample(addr!("global"), Normal::new(0.0, 1.0).unwrap());
            let locals <- plate!(i in 0..2 => {
                sample(scoped_addr!("local", "param", "{}", i),
                       Normal::new(global, 0.1).unwrap())
            });
            pure((global, locals))
        };

        // All models should construct successfully
        // (Actual execution would require handlers)
    }

Common Pitfalls

1. Address Conflicts: Use scoped_addr! for complex models
2. Memory Usage: Large plate operations can create big traces
3. Sequential Dependencies: Explicit state management required
4. Type Inference: Sometimes need explicit type annotations

Performance Considerations

• Plate Size: Very large plates may exceed memory limits
• Nesting Depth: Deep hierarchies increase trace size
• Address Complexity: Simple addresses are more efficient
• Function Composition: Pure functions are optimized away

Next Steps

• Optimization: See Optimizing Performance for efficiency techniques
• Debugging: Check Debugging Models for troubleshooting complex models
• Production: Learn Production Deployment for scaling

Complex models become tractable and maintainable through systematic composition, principled addressing, and mathematical abstraction. Fugue's macro system provides elegant syntactic sugar while preserving the underlying categorical structure that enables powerful inference algorithms and compositional reasoning about probabilistic programs.

Optimizing Performance

Performance optimization in probabilistic programming requires understanding both computational complexity and numerical analysis. This guide explores Fugue's systematic approach to memory optimization, numerical stability, and algorithmic efficiency for production-scale probabilistic workloads.

Memory-Optimized Inference

Memory allocation becomes the computational bottleneck in high-throughput scenarios due to garbage collection overhead. The allocation rate $R_{\text{alloc}}$ for naive inference scales as:

$$R_{\text{alloc}}=n\cdot |T|\cdot f_{\text{gc}}$$

where $n$ is the sample count, $|T|$ is the trace size, and $f_{\text{gc}}$ is the GC frequency. Fugue's object pooling reduces this to $\mathcal{O}(1)$ after warmup:

graph TD
    subgraph "Traditional Allocation"
        A1["Sample 1"] --> B1["Allocate Trace"]
        B1 --> C1["GC Pressure"]
        A2["Sample 2"] --> B2["Allocate Trace"]  
        B2 --> C2["GC Pressure"]
        A3["Sample n"] --> B3["Allocate Trace"]
        B3 --> C3["GC Pressure"]
    end
    
    subgraph "Pooled Allocation"  
        D1["Sample 1"] --> E1["Reuse from Pool"]
        D2["Sample 2"] --> E1
        D3["Sample n"] --> E1
        E1 --> F["Zero GC Pressure"]
    end

    // Create trace pool for zero-allocation inference
    let mut pool = TracePool::new(50); // Pool up to 50 traces
    let mut rng = thread_rng();

    // Define a model that would normally cause many allocations
    let make_model = || {
        prob!(
            let x <- sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
            let y <- sample(addr!("y"), Normal::new(x, 0.5).unwrap());
            observe(addr!("obs"), Normal::new(y, 0.1).unwrap(), 1.5);
            pure(x)
        )
    };

    // Time pooled vs non-pooled execution
    let start = Instant::now();
    for _iteration in 0..1000 {
        // Use pooled handler for efficient memory reuse
        let (_result, trace) =
            runtime::handler::run(PooledPriorHandler::new(&mut rng, &mut pool), make_model());
        // Return trace to pool for reuse
        pool.return_trace(trace);
    }
    let pooled_time = start.elapsed();

    let stats = pool.stats();
    println!("✅ Completed 1000 iterations with memory pooling");
    println!("   - Execution time: {:?}", pooled_time);
    println!("   - Hit ratio: {:.1}%", stats.hit_ratio());
    println!(
        "   - Pool stats - hits: {}, misses: {}",
        stats.hits, stats.misses
    );

Key Benefits:

• Zero-allocation execution after warm-up
• Configurable pool size for memory control
• Automatic trace recycling and cleanup
• Built-in performance monitoring with hit ratios

Numerical Stability

Numerical stability in probabilistic computing requires careful analysis of condition numbers and floating-point precision. The log-sum-exp operation is fundamental:

$$\text{LSE}(\mathbf{x})=\log \left(\sum _{i=1}^n{e}^{x_i}\right)=x_{\max }+\log \left(\sum _{i=1}^n{e}^{x_i-x_{\max }}\right)$$

Stability Analysis: Direct computation of $\sum e^{x_i}$ has condition number $\kappa \approx e^{x_{\max }-x_{\min }}$, which becomes ill-conditioned when $x_{\max }-x_{\min }\gg \log ({ϵ}_{\text{machine}})$.

    // Demonstrate stable log-probability computations
    let extreme_log_probs = vec![700.0, 701.0, 699.0, 698.0]; // Would overflow in linear space

    // Safe log-sum-exp prevents overflow
    let log_normalizer = log_sum_exp(&extreme_log_probs);
    let normalized_probs = normalize_log_probs(&extreme_log_probs);

    println!("✅ Stable computation with extreme log-probabilities");
    println!("   - Log normalizer: {:.2}", log_normalizer);
    println!(
        "   - Probabilities sum to: {:.10}",
        normalized_probs.iter().sum::<f64>()
    );

    // Weighted log-sum-exp for importance sampling
    let log_values = vec![-1.0, -2.0, -3.0, -4.0];
    let weights = vec![0.4, 0.3, 0.2, 0.1];
    let weighted_result = weighted_log_sum_exp(&log_values, &weights);

    println!("   - Weighted log-sum-exp: {:.4}", weighted_result);

    // Safe logarithm handling
    let safe_results: Vec<f64> = [1.0, 0.0, -1.0].iter().map(|&x| safe_ln(x)).collect();
    println!("   - Safe ln results: {:?}", safe_results);

Stability Features:

• log_sum_exp prevents overflow in mixture computations
• weighted_log_sum_exp for importance sampling
• safe_ln handles edge cases gracefully
• All operations maintain numerical precision across scales

Efficient Trace Construction

When building traces programmatically, use TraceBuilder for optimal performance:

    // Use TraceBuilder for efficient trace creation
    let mut builder = TraceBuilder::new();

    let start = Instant::now();
    for i in 0..100 {
        // Add choices efficiently without reallocations
        builder.add_sample(
            addr!("param", i),
            i as f64,
            0.0, // log_prob
        );
    }

    // Build final trace efficiently
    let constructed_trace = builder.build();
    let construction_time = start.elapsed();

    println!("✅ Efficient trace construction");
    println!(
        "   - Built trace with {} choices in {:?}",
        constructed_trace.choices.len(),
        construction_time
    );
    println!(
        "   - Total log weight: {:.2}",
        constructed_trace.total_log_weight()
    );

Construction Benefits:

• Pre-allocated data structures minimize reallocations
• Type-specific insertion methods avoid boxing overhead
• Batch operations for multiple choices
• Efficient conversion to immutable traces

Copy-on-Write for MCMC

MCMC algorithms exhibit temporal locality in parameter updates, modifying only $\mathcal{O}(\log d)$ parameters per iteration where $d$ is the total dimensionality. Copy-on-Write (COW) data structures exploit this pattern:

graph TD
    subgraph "MCMC Iteration Structure"
        A["Base Trace T₀"] --> B{"Proposal Step"}
        B --> C["Modified Parameters δ"]
        C --> D{"Small Changes?"}
        D -->|Yes| E["COW: Share + Δ"]
        D -->|No| F["Full Copy"]
        E --> G["O(1) Memory"]
        F --> H["O(d) Memory"]
    end

Complexity Analysis: Traditional MCMC requires $\mathcal{O}(d)$ space per sample. COW reduces this to $\mathcal{O}(\Delta +\log d)$ where $\Delta$ is the edit distance between traces.

    // Create base trace manually for MCMC
    let mut builder = TraceBuilder::new();
    builder.add_sample(addr!("mu"), 0.5, -0.5);
    builder.add_sample(addr!("sigma"), 1.0, -1.0);
    builder.add_sample_bool(addr!("component"), true, -0.69);
    let base_trace = builder.build();

    // Create COW trace for efficient copying
    let cow_base = CowTrace::from_trace(base_trace);

    let start = Instant::now();
    let mut mcmc_traces = Vec::new();

    for _proposal in 0..1000 {
        // Clone is O(1) until modification
        let mut proposal_trace = cow_base.clone();

        // Modify only one parameter (triggers COW)
        proposal_trace.insert_choice(
            addr!("mu"),
            Choice {
                addr: addr!("mu"),
                value: ChoiceValue::F64(0.6),
                logp: -0.4,
            },
        );

        mcmc_traces.push(proposal_trace);
    }
    let cow_time = start.elapsed();

    println!("✅ Copy-on-write MCMC proposals");
    println!("   - Created 1000 proposal traces in {:?}", cow_time);
    println!("   - Memory sharing until modification");

MCMC Optimizations:

• O(1) trace cloning until modification
• Shared memory for unchanged parameters
• Lazy copying only when traces diverge
• Perfect for Metropolis-Hastings and Gibbs sampling

Vectorized Model Patterns

Structure models for efficient batch processing:

    // Pre-allocate data structures for repeated use
    let observations: Vec<f64> = (0..100).map(|i| i as f64 * 0.1).collect();
    let n = observations.len();

    // Efficient vectorized model
    let vectorized_model = || {
        prob!(
            let mu <- sample(addr!("global_mu"), Normal::new(0.0, 10.0).unwrap());
            let precision <- sample(addr!("precision"), Gamma::new(2.0, 1.0).unwrap());
            let sigma = (1.0 / precision).sqrt();

            // Use plate for efficient vectorized operations
            let _likelihoods <- plate!(i in 0..n => {
                observe(addr!("obs", i), Normal::new(mu, sigma).unwrap(), observations[i])
            });

            pure((mu, sigma))
        )
    };

    let start = Instant::now();
    let (_result, _trace) = runtime::handler::run(
        PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        vectorized_model(),
    );
    let vectorized_time = start.elapsed();

    println!("✅ Optimized vectorized model");
    println!("   - Processed {} observations in {:?}", n, vectorized_time);

Vectorization Strategy:

• Pre-allocate data collections
• Use plate! for independent parallel operations
• Minimize dynamic allocations in hot paths
• Leverage compiler optimizations with static sizing

Performance Monitoring

Systematic performance monitoring requires tracking multiple performance metrics with their theoretical bounds:

$$\begin{array}{rl}\text{Throughput}& =\frac{\text{samples}}{\text{time}}\le \frac{1}{{\tau }_{\min }}\\ \text{Latency}& =\text{time per sample}\ge {\tau }_{\min }\\ \text{Memory Efficiency}& =\frac{\text{useful allocations}}{\text{total allocations}}\to 1\end{array}$$

where ${\tau }_{\min }$ is the theoretical minimum execution time per sample.

    // Monitor trace characteristics for optimization insights
    #[derive(Debug)]
    struct TraceMetrics {
        num_choices: usize,
        log_weight: f64,
        is_valid: bool,
        memory_size_estimate: usize,
    }

    impl TraceMetrics {
        fn from_trace(trace: &Trace) -> Self {
            let num_choices = trace.choices.len();
            let log_weight = trace.total_log_weight();
            let is_valid = log_weight.is_finite();

            // Rough memory estimate (actual implementation would be more precise)
            let memory_size_estimate = num_choices * 64; // Rough bytes per choice

            Self {
                num_choices,
                log_weight,
                is_valid,
                memory_size_estimate,
            }
        }
    }

    // Example: Monitor a complex model's performance
    let complex_model = || {
        prob!(
            let components <- plate!(c in 0..5 => {
                sample(addr!("weight", c), Gamma::new(1.0, 1.0).unwrap())
                    .bind(move |weight| {
                        sample(addr!("mu", c), Normal::new(0.0, 2.0).unwrap())
                            .map(move |mu| (weight, mu))
                    })
            });

            let selector <- sample(addr!("selector"),
                                  Categorical::new(vec![0.2, 0.2, 0.2, 0.2, 0.2]).unwrap());

            pure((components, selector))
        )
    };

    let (_result, trace) = runtime::handler::run(
        PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        complex_model(),
    );

    let metrics = TraceMetrics::from_trace(&trace);
    println!("✅ Performance monitoring active");
    println!("   - Trace choices: {}", metrics.num_choices);
    println!("   - Log weight: {:.2}", metrics.log_weight);
    println!("   - Valid: {}", metrics.is_valid);
    println!(
        "   - Memory estimate: {} bytes",
        metrics.memory_size_estimate
    );

Monitoring Approach:

• Collect trace characteristics for optimization insights
• Track memory usage patterns
• Validate numerical stability
• Profile execution bottlenecks

Batch Processing

Batch processing amortizes setup costs and exploits hardware parallelism. The optimal batch size $b^{\ast }$ balances memory usage and throughput:

$$b^{\ast }=\arg \underset{b}{\min }\left(\frac{C_{\text{setup}}}{b}+b\cdot C_{\text{memory}}+\frac{C_{\text{sync}}}{b}\right)$$

where:

• $C_{\text{setup}}$ is the per-batch initialization cost
• $C_{\text{memory}}$ is the per-sample memory cost
• $C_{\text{sync}}$ is the synchronization overhead

graph LR
    subgraph "Performance vs Batch Size"
        A["Small Batches<br/>b → 1"] --> B["High Setup<br/>Overhead"]
        C["Large Batches<br/>b → ∞"] --> D["Memory<br/>Pressure"]  
        E["Optimal Batch<br/>b*"] --> F["Balanced<br/>Performance"]
    end

    // Efficient batch inference using memory pooling
    let batch_size = 100;
    let mut batch_pool = TracePool::new(batch_size);

    let start = Instant::now();
    let mut batch_results = Vec::with_capacity(batch_size);

    for _batch in 0..batch_size {
        let (result, trace) = runtime::handler::run(
            PooledPriorHandler::new(&mut rng, &mut batch_pool),
            make_model(),
        );
        // Return trace to pool for reuse
        batch_pool.return_trace(trace);
        batch_results.push(result);
    }

    let batch_time = start.elapsed();
    let batch_stats = batch_pool.stats();

    println!("✅ Batch processing complete");
    println!("   - Processed {} samples in {:?}", batch_size, batch_time);
    println!(
        "   - Average time per sample: {:?}",
        batch_time / batch_size as u32
    );
    println!(
        "   - Memory efficiency: {:.1}% hit ratio",
        batch_stats.hit_ratio()
    );

Batch Benefits:

• Amortized setup costs across samples
• Memory pool reuse for consistent performance
• Scalable to large sample counts
• Predictable memory footprint

Numerical Precision Testing

Validate stability across different computational scales:

    // Test numerical stability across different scales
    let test_scales = vec![1e-10, 1e-5, 1.0, 1e5, 1e10];

    for &scale in &test_scales {
        let scale: f64 = scale;
        let log_vals = vec![scale.ln() + 1.0, scale.ln() + 2.0, scale.ln() + 0.5];

        let stable_sum = log_sum_exp(&log_vals);
        let log1p_result = log1p_exp(scale.ln());

        println!(
            "   Scale {:.0e}: log_sum_exp={:.4}, log1p_exp={:.4}",
            scale, stable_sum, log1p_result
        );
    }
    println!("✅ Numerical stability verified across scales");

Testing Strategy:

• Verify stability across extreme value ranges
• Test edge cases and boundary conditions
• Validate consistency of numerical operations
• Profile precision vs. performance trade-offs

Performance Testing

Implement systematic performance validation:

    #[test]
    fn test_memory_pool_efficiency() {
        let mut pool = TracePool::new(10);
        let mut rng = thread_rng();

        // Test pool reuse with PooledPriorHandler
        for _i in 0..20 {
            let (_, trace) = runtime::handler::run(
                PooledPriorHandler::new(&mut rng, &mut pool),
                sample(addr!("test"), Normal::new(0.0, 1.0).unwrap()),
            );
            // Return trace to pool for reuse
            pool.return_trace(trace);
        }

        let stats = pool.stats();
        assert!(
            stats.hit_ratio() > 50.0,
            "Pool should have good hit ratio, got {:.1}%",
            stats.hit_ratio()
        );
        assert!(stats.hits + stats.misses > 0, "Pool should have been used");
    }

    #[test]
    fn test_numerical_stability() {
        // Test log_sum_exp with extreme values
        let extreme_vals = vec![700.0, 701.0, 699.0];
        let result = log_sum_exp(&extreme_vals);
        assert!(
            result.is_finite(),
            "log_sum_exp should handle extreme values"
        );

        // Test normalization
        let normalized = normalize_log_probs(&extreme_vals);
        let sum: f64 = normalized.iter().sum();
        assert!(
            (sum - 1.0).abs() < 1e-10,
            "Normalized probabilities should sum to 1"
        );

        // Test weighted computation
        let weights = vec![0.5, 0.3, 0.2];
        let weighted_result = weighted_log_sum_exp(&extreme_vals, &weights);
        assert!(
            weighted_result.is_finite(),
            "Weighted log_sum_exp should be finite"
        );
    }

    #[test]
    fn test_trace_builder_efficiency() {
        let mut builder = TraceBuilder::new();

        // Add many choices efficiently
        for i in 0..100 {
            builder.add_sample(addr!("param", i), i as f64, -0.5);
        }

        let trace = builder.build();
        assert_eq!(trace.choices.len(), 100);
        assert!(trace.total_log_weight().is_finite());
    }

    #[test]
    fn test_cow_trace_sharing() {
        // Create base trace using builder
        let mut builder = TraceBuilder::new();
        builder.add_sample(addr!("x"), 1.0, -0.5);
        let base = builder.build();
        let cow_trace = CowTrace::from_trace(base);

        // Clone should be fast
        let clone1 = cow_trace.clone();
        let clone2 = cow_trace.clone();

        // Should share data until modification - convert to regular trace to test
        let trace1 = clone1.to_trace();
        let trace2 = clone2.to_trace();
        assert_eq!(trace1.get_f64(&addr!("x")), Some(1.0));
        assert_eq!(trace2.get_f64(&addr!("x")), Some(1.0));
    }

Testing Framework:

• Memory pool efficiency validation
• Numerical stability regression tests
• Trace construction benchmarking
• COW sharing verification

Production Deployment

Memory Configuration

• Size TracePool based on peak concurrent inference
• Monitor hit ratios to validate pool efficiency
• Use COW traces for MCMC workloads
• Pre-warm pools before production traffic

Numerical Strategies

• Always use log-space for probability computations
• Validate extreme value handling in testing
• Monitor for numerical instabilities in production
• Use stable algorithms for critical computations

Monitoring and Alerting

• Track inference latency and memory usage
• Monitor pool statistics and efficiency metrics
• Alert on numerical instabilities or performance degradation
• Profile hot paths for optimization opportunities

Common Performance Patterns

1. Pool First: Use TracePool for any repeated inference
2. Log Always: Work in log-space for numerical stability
3. Batch Everything: Amortize costs across multiple samples
4. Monitor Continuously: Track performance metrics in production
5. Test Extremes: Validate stability with extreme values

These optimization strategies enable Fugue to handle production-scale probabilistic programming workloads with consistent performance and numerical reliability.

Debugging Models

Debugging probabilistic models presents unique challenges due to their stochastic nature and high-dimensional parameter spaces. Unlike deterministic programs, probabilistic models require statistical validation, convergence analysis, and distributional testing. This guide establishes a systematic methodology for probabilistic model debugging using Fugue's comprehensive diagnostic framework.

Trace Inspection and Analysis

Execution traces form the foundation of probabilistic model debugging. Each trace $\mathcal{T}$ contains a complete record of the program's stochastic execution:

$$\mathcal{T}=\{(a_i,v_i,w_i){\}}_{i=1}^n$$

where $a_i$ is the address, $v_i$ is the sampled value, and $w_i$ is the log-weight contribution.

graph TD
    subgraph "Trace Analysis Workflow"
        A[Execution Trace T] --> B{Finite Log-Weight?}
        B -->|No| C[Numerical Instability]
        B -->|Yes| D[Choice Analysis]
        D --> E[Parameter Extraction]
        E --> F[Statistical Validation]
        F --> G{Passes Tests?}
        G -->|No| H[Model Refinement]
        G -->|Yes| I[Model Validated]
        C --> J[Debug Constraints]
        H --> A
        J --> A
    end

Mathematical Properties: A valid trace must satisfy the weight consistency condition: $$\log P(\mathcal{T})=\sum _{i=1}^n{w}_i<\infty$$

    // Execute a model and examine its trace structure
    let mut rng = thread_rng();

    let diagnostic_model = || {
        prob!(
            let mu <- sample(addr!("mu"), Normal::new(0.0, 2.0).unwrap());
            let sigma <- sample(addr!("sigma"), Gamma::new(2.0, 1.0).unwrap());
            observe(addr!("obs1"), Normal::new(mu, sigma).unwrap(), 1.5);
            observe(addr!("obs2"), Normal::new(mu, sigma).unwrap(), 1.2);
            factor(if mu.abs() < 3.0 { 0.0 } else { f64::NEG_INFINITY });
            pure((mu, sigma))
        )
    };

    let ((mu_val, sigma_val), trace) = runtime::handler::run(
        PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        diagnostic_model(),
    );

    println!("✅ Model execution complete");
    println!("   - Result: mu = {:.3}, sigma = {:.3}", mu_val, sigma_val);
    println!("   - Choices recorded: {}", trace.choices.len());
    println!("   - Prior log-weight: {:.6}", trace.log_prior);
    println!("   - Likelihood log-weight: {:.6}", trace.log_likelihood);
    println!("   - Factor log-weight: {:.6}", trace.log_factors);
    println!("   - Total log-weight: {:.6}", trace.total_log_weight());

    // Per-choice breakdown
    println!("   - Choice breakdown:");
    for (addr, choice) in &trace.choices {
        println!(
            "     {}: {:?} (logp: {:.6})",
            addr, choice.value, choice.logp
        );
    }

Key Debugging Insights:

• Choice count reveals model complexity and structure
• Log-weight decomposition identifies prior vs. likelihood vs. factor issues
• Per-choice analysis shows individual parameter contributions
• Finite log-weights indicate valid model execution

Type-Safe Value Access

Fugue provides robust access patterns that handle type mismatches gracefully:

    // Safe access patterns that handle type mismatches gracefully

    // Option-based access (returns None on mismatch)
    match trace.get_f64(&addr!("mu")) {
        Some(mu) => println!("✅ Retrieved mu = {:.3}", mu),
        None => println!("❌ Failed to get mu as f64"),
    }

    // Result-based access (returns detailed error info)
    match trace.get_f64_result(&addr!("sigma")) {
        Ok(sigma) => println!("✅ Retrieved sigma = {:.3}", sigma),
        Err(e) => println!("❌ Error getting sigma: {}", e),
    }

    // Handle missing addresses
    match trace.get_f64_result(&addr!("missing_param")) {
        Ok(_) => unreachable!(),
        Err(e) => println!("✅ Correctly caught missing address: {}", e),
    }

    // Handle type mismatches
    match trace.get_bool_result(&addr!("mu")) {
        Ok(_) => unreachable!(),
        Err(e) => println!("✅ Correctly caught type mismatch: {}", e),
    }

    // Iterate through all choices for debugging
    println!("   - All choices and their types:");
    for (addr, choice) in &trace.choices {
        let type_info = match &choice.value {
            ChoiceValue::F64(_) => "f64",
            ChoiceValue::Bool(_) => "bool",
            ChoiceValue::U64(_) => "u64",
            ChoiceValue::I64(_) => "i64",
            ChoiceValue::Usize(_) => "usize",
        };
        println!("     {} ({}): {:?}", addr, type_info, choice.value);
    }

Error Handling Strategies:

• Use get_*_result() for detailed error information
• Use get_*() for simple None-handling
• Always check for missing addresses before assuming success
• Iterate through all choices to understand model structure

Model Validation and Testing

Systematic validation ensures your model behaves as expected:

    // Test a simple conjugate model against analytical solution
    let conjugate_model = || {
        prob!(
            let theta <- sample(addr!("theta"), Beta::new(1.0, 1.0).unwrap());
            observe(addr!("successes"), Binomial::new(10, theta).unwrap(), 7u64);
            pure(theta)
        )
    };

    // Run a few samples to test basic functionality
    let mut theta_samples = Vec::new();
    for _ in 0..20 {
        let (theta, test_trace) = runtime::handler::run(
            PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            },
            conjugate_model(),
        );

        // Validate trace structure
        assert!(test_trace.choices.contains_key(&addr!("theta")));
        assert!(
            test_trace.total_log_weight().is_finite(),
            "Trace should have finite log-weight"
        );
        assert!(
            test_trace.log_likelihood.is_finite(),
            "Likelihood should be finite"
        );

        theta_samples.push(theta);
    }

    // Basic statistical checks
    let sample_mean = theta_samples.iter().sum::<f64>() / theta_samples.len() as f64;
    println!("✅ Validation tests passed");
    println!("   - Generated {} samples", theta_samples.len());
    println!(
        "   - Sample mean: {:.3} (expected ~0.7 for Beta-Binomial)",
        sample_mean
    );
    println!("   - All traces had finite log-weights");

Validation Best Practices:

• Test against known analytical solutions
• Verify all traces have finite log-weights
• Check basic statistical properties (means, variances)
• Test edge cases and boundary conditions

Safe vs Strict Error Handling

Fugue provides both strict (fail-fast) and safe (error-resilient) execution modes:

    // Create a trace with known structure for replay testing
    let mut base_trace = Trace::default();
    base_trace.insert_choice(addr!("param"), ChoiceValue::F64(1.5), -0.5);

    let test_model = || sample(addr!("param"), Normal::new(0.0, 1.0).unwrap());

    // Strict replay - will panic on mismatch (commented out for safety)
    // let strict_replay = ReplayHandler { base_trace: &base_trace };
    // let (strict_result, strict_trace) = runtime::handler::run(strict_replay, test_model());

    // Safe replay - handles errors gracefully
    let safe_replay = SafeReplayHandler {
        rng: &mut rng,
        base: base_trace.clone(),
        trace: Trace::default(),
        warn_on_mismatch: true,
    };
    let (safe_result, safe_trace) = runtime::handler::run(safe_replay, test_model());

    println!("✅ Safe replay succeeded");
    println!("   - Result: {:.3}", safe_result);
    println!(
        "   - Retrieved value: {:?}",
        safe_trace.get_f64(&addr!("param"))
    );

    // Test scoring with safe handler
    let safe_score = SafeScoreGivenTrace {
        base: base_trace,
        trace: Trace::default(),
        warn_on_error: false,
    };
    let (_, score_trace) = runtime::handler::run(safe_score, test_model());

    println!(
        "   - Score trace log-weight: {:.3}",
        score_trace.total_log_weight()
    );

When to Use Each:

• Strict handlers (ReplayHandler, ScoreGivenTrace): Development and testing
• Safe handlers (SafeReplayHandler, SafeScoreGivenTrace): Production systems
• Safe handlers log warnings instead of panicking on mismatches

MCMC Diagnostics

Markov Chain Monte Carlo convergence assessment requires statistical hypothesis testing and diagnostic metrics. The fundamental question is whether the chain has reached its stationary distribution $\pi (\theta )$.

Gelman-Rubin Diagnostic

The potential scale reduction factor $\hat{R}$ compares within-chain and between-chain variance:

$$\hat{R}=\sqrt{\frac{\hat{V}}{W}}$$

where:

• $W=\frac{1}{m}{\sum }_{j=1}^m{s}_j^2$ (within-chain variance)
• $B=\frac{n}{m-1}{\sum }_{j=1}^m({\bar{\theta }}_{j\cdot }-{\bar{\theta }}_{\cdot \cdot }{)}^2$ (between-chain variance)
• $\hat{V}=\frac{n-1}{n}W+\frac{1}{n}B$ (marginal posterior variance estimate)

Effective Sample Size

The effective sample size accounts for autocorrelation in MCMC samples:

$$\text{ESS}=\frac{mn}{1+2{\sum }_{t=1}^{\infty }{\rho }_t}$$

where ${\rho }_t$ is the lag-$t$ autocorrelation and $mn$ is the total number of samples.

    // Generate simple MCMC chains for diagnostic testing
    let mcmc_model = || {
        prob!(
            let mu <- sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap());
            observe(addr!("y"), Normal::new(mu, 0.5).unwrap(), 1.0);
            pure(mu)
        )
    };

    // Generate two short chains for R-hat calculation
    let n_samples = 50;
    let n_warmup = 10;

    let mut chain1_samples = Vec::new();
    let mut chain2_samples = Vec::new();

    // Chain 1
    let mut rng1 = rand::rngs::StdRng::seed_from_u64(42);
    let chain1 = adaptive_mcmc_chain(&mut rng1, mcmc_model, n_samples, n_warmup);
    for (_, trace) in &chain1 {
        if let Some(mu) = trace.get_f64(&addr!("mu")) {
            chain1_samples.push(mu);
        }
    }

    // Chain 2
    let mut rng2 = rand::rngs::StdRng::seed_from_u64(123);
    let chain2 = adaptive_mcmc_chain(&mut rng2, mcmc_model, n_samples, n_warmup);
    for (_, trace) in &chain2 {
        if let Some(mu) = trace.get_f64(&addr!("mu")) {
            chain2_samples.push(mu);
        }
    }

    // Compute diagnostics
    if !chain1_samples.is_empty() && !chain2_samples.is_empty() {
        // Extract traces for R-hat calculation
        let chain1_traces: Vec<Trace> = chain1.into_iter().map(|(_, trace)| trace).collect();
        let chain2_traces: Vec<Trace> = chain2.into_iter().map(|(_, trace)| trace).collect();
        let r_hat = r_hat_f64(&[chain1_traces, chain2_traces], &addr!("mu"));
        let ess1 = effective_sample_size_mcmc(&chain1_samples);
        let ess2 = effective_sample_size_mcmc(&chain2_samples);

        println!("✅ MCMC diagnostics computed");
        println!(
            "   - Chain 1: {} samples, ESS = {:.1}",
            chain1_samples.len(),
            ess1
        );
        println!(
            "   - Chain 2: {} samples, ESS = {:.1}",
            chain2_samples.len(),
            ess2
        );
        println!("   - R-hat: {:.4} (< 1.1 indicates convergence)", r_hat);

        if r_hat < 1.1 {
            println!("   - ✅ Chains appear to have converged");
        } else {
            println!("   - ⚠️  Chains may not have converged - run longer");
        }
    }

Convergence Indicators:

• R-hat < 1.1: Chains have converged
• High ESS: Efficient sampling without excessive correlation
• Multiple chains: Essential for reliable convergence assessment
• Visual inspection: Always examine trace plots when possible

Model Structure Analysis

Model structure analysis reveals the computational graph and parameter dependencies. This analysis is crucial for understanding model complexity and identifying potential issues:

graph TD
    subgraph "Model Structure Hierarchy"
        A[Model M] --> B[Parameter Groups]
        B --> C1[Hyperpriors θ₁]
        B --> C2[Primary Parameters θ₂] 
        B --> C3[Observations y]
        C1 --> D1[Constraint Analysis]
        C2 --> D2[Dependency Graph]
        C3 --> D3[Likelihood Terms]
        D1 --> E[Structure Validation]
        D2 --> E
        D3 --> E
    end

Structural Invariants to validate:

1. Address Uniqueness: $|\{a_i\}|=n$ (no collisions)
2. Parameter Hierarchy: $\forall i,j:a_i⪯a_j⟹\text{dependency}(i,j)$
3. Choice Count Consistency: Expected vs. actual parameter count
4. Type Safety: Each address maps to consistent value types

    // Create a complex model to demonstrate structure analysis
    let complex_model = || {
        prob!(
            // Hierarchical structure
            let global_scale <- sample(addr!("global_scale"), Gamma::new(2.0, 1.0).unwrap());

            let group_params <- plate!(g in 0..3 => {
                sample(addr!("group_mean", g), Normal::new(0.0, global_scale).unwrap())
                    .bind(move |mean| {
                        sample(addr!("group_precision", g), Gamma::new(2.0, 1.0).unwrap())
                            .map(move |prec| (mean, prec))
                    })
            });

            // Individual observations (simplified to avoid move issues)
            let observations = [1.2, 1.5, 0.8];
            let likelihoods <- plate!(i in 0..observations.len() => {
                // Use fixed parameters for demonstration
                observe(addr!("obs", i), Normal::new(0.0, 1.0).unwrap(), observations[i])
            });

            pure((global_scale, group_params, likelihoods))
        )
    };

    let (_result, complex_trace) = runtime::handler::run(
        PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        complex_model(),
    );

    // Analyze model structure
    let mut address_analysis = BTreeMap::new();
    for (addr, choice) in &complex_trace.choices {
        let addr_str = addr.0.clone();
        let category = if addr_str.contains("global") {
            "Global Parameters"
        } else if addr_str.contains("group") {
            "Group Parameters"
        } else if addr_str.contains("obs") {
            "Observations"
        } else {
            "Other"
        };

        address_analysis
            .entry(category)
            .or_insert(Vec::new())
            .push((addr_str, choice.logp));
    }

    println!("✅ Complex model structure analysis");
    println!("   - Total choices: {}", complex_trace.choices.len());
    println!("   - Address structure:");
    for (category, addresses) in address_analysis {
        println!("     {}: {} choices", category, addresses.len());
        for (addr, logp) in addresses.iter().take(3) {
            // Show first 3
            println!("       {} (logp: {:.3})", addr, logp);
        }
        if addresses.len() > 3 {
            println!("       ... and {} more", addresses.len() - 3);
        }
    }

Structure Analysis Benefits:

• Understand parameter organization and hierarchies
• Detect unexpected address patterns
• Verify choice counts match model expectations
• Identify bottlenecks in complex models

Performance Diagnostics

Monitor computational efficiency and identify bottlenecks:

    use std::time::Instant;

    // Benchmark model execution and trace construction
    let benchmark_model = || {
        prob!(
            let params <- plate!(i in 0..100 => {
                sample(addr!("param", i), Normal::new(0.0, 1.0).unwrap())
            });
            pure(params)
        )
    };

    let start = Instant::now();
    let (_, bench_trace) = runtime::handler::run(
        PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        benchmark_model(),
    );
    let execution_time = start.elapsed();

    // Analyze trace characteristics
    let choice_count = bench_trace.choices.len();
    let memory_estimate = choice_count * 64; // Rough estimate
    let log_weight_is_finite = bench_trace.total_log_weight().is_finite();

    println!("✅ Performance diagnostics");
    println!("   - Execution time: {:?}", execution_time);
    println!("   - Choices created: {}", choice_count);
    println!("   - Memory estimate: ~{} bytes", memory_estimate);
    println!("   - Log-weight valid: {}", log_weight_is_finite);

    // Check for potential issues
    if choice_count == 0 {
        println!("   - ⚠️  No choices recorded - possible model issue");
    }
    if !log_weight_is_finite {
        println!("   - ⚠️  Invalid log-weight - check factors and observations");
    }
    if execution_time.as_millis() > 100 {
        println!("   - ⚠️  Slow execution - consider optimization");
    }

Performance Warning Signs:

• Zero choices recorded (model execution failure)
• Infinite log-weights (constraint violations)
• Excessive execution time (optimization needed)
• Large memory footprint (consider streaming approaches)

Common Debugging Patterns

Systematic debugging follows a hierarchical validation strategy from basic correctness to statistical validity:

graph TD
    subgraph "Debugging Methodology"
        A[Model Implementation] --> B{Syntax Valid?}
        B -->|No| C[Fix Code Structure]
        B -->|Yes| D{Types Consistent?}
        D -->|No| E[Fix Type Errors]
        D -->|Yes| F{Finite Log-Weights?}
        F -->|No| G[Fix Constraints]
        F -->|Yes| H{Statistical Properties?}
        H -->|No| I[Validate Distributions]
        H -->|Yes| J{Convergence?}
        J -->|No| K[Tune Inference]
        J -->|Yes| L[Model Validated]
        
        C --> A
        E --> A
        G --> A
        I --> A
        K --> A
    end

Debug Level Hierarchy:

1. Syntactic: Code compiles and types check
2. Semantic: Model executes without runtime errors
3. Numerical: Computations remain stable and finite
4. Statistical: Results match theoretical expectations
5. Convergence: Inference algorithms reach stationarity

    // Pattern 1: Systematic model testing
    fn test_model_basic_properties<F, T>(
        model_fn: F,
        expected_choice_count: usize,
        description: &str,
    ) where
        F: Fn() -> Model<T>,
    {
        let mut rng = thread_rng();
        let (_, trace) = runtime::handler::run(
            PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            },
            model_fn(),
        );

        println!("Testing {}", description);

        // Basic trace validity
        assert!(
            trace.total_log_weight().is_finite(),
            "Log-weight should be finite"
        );
        assert_eq!(
            trace.choices.len(),
            expected_choice_count,
            "Choice count mismatch"
        );

        // Check for common issues
        if trace.log_prior.is_infinite() {
            println!("  - ⚠️  Infinite prior - check parameter ranges");
        }
        if trace.log_likelihood.is_infinite() {
            println!("  - ⚠️  Infinite likelihood - check observations");
        }
        if trace.log_factors.is_infinite() {
            println!("  - ⚠️  Infinite factors - check constraint satisfaction");
        }

        println!("  - ✅ {} passed basic tests", description);
    }

    // Test simple models
    test_model_basic_properties(
        || sample(addr!("x"), Normal::new(0.0, 1.0).unwrap()),
        1,
        "Simple normal sampling",
    );

    test_model_basic_properties(
        || {
            prob!(
                let x <- sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
                observe(addr!("y"), Normal::new(x, 0.5).unwrap(), 1.0);
                pure(x)
            )
        },
        1,
        "Normal model with observation",
    );

    // Pattern 2: Address collision detection
    fn check_address_collisions(trace: &Trace) -> Vec<String> {
        let mut collisions = Vec::new();
        let addresses: Vec<&str> = trace.choices.keys().map(|addr| addr.0.as_str()).collect();

        for (i, addr1) in addresses.iter().enumerate() {
            for addr2 in addresses.iter().skip(i + 1) {
                if addr1 == addr2 {
                    collisions.push(format!("Duplicate address: {}", addr1));
                }
            }
        }
        collisions
    }

    let test_trace = complex_trace; // Use trace from earlier
    let collisions = check_address_collisions(&test_trace);
    if collisions.is_empty() {
        println!("  - ✅ No address collisions detected");
    } else {
        for collision in collisions {
            println!("  - ⚠️  {}", collision);
        }
    }

    println!("✅ Debugging patterns demonstration complete");

Debugging Workflow:

1. Start Simple: Test individual components before complex composition
2. Validate Incrementally: Add complexity one piece at a time
3. Check Address Uniqueness: Prevent parameter collision bugs
4. Monitor Log-Weights: Track prior, likelihood, and factor contributions
5. Use Systematic Testing: Automated validation for all model components

Testing Framework Integration

Embed debugging checks in your test suite:

    #[test]
    fn test_trace_inspection_patterns() {
        let mut rng = thread_rng();

        let model = prob!(
            let x <- sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
            let y <- sample(addr!("y"), Beta::new(1.0, 1.0).unwrap());
            observe(addr!("obs"), Normal::new(x, 0.1).unwrap(), 1.5);
            pure((x, y))
        );

        let (_result, trace) = runtime::handler::run(
            PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            },
            model,
        );

        // Basic trace properties
        assert_eq!(trace.choices.len(), 2); // x and y samples
        assert!(trace.total_log_weight().is_finite());
        assert!(trace.log_likelihood.is_finite());

        // Type-safe access
        assert!(trace.get_f64(&addr!("x")).is_some());
        assert!(trace.get_f64(&addr!("y")).is_some());
        assert!(trace.get_bool(&addr!("x")).is_none()); // Type mismatch

        // Result access patterns
        assert!(trace.get_f64_result(&addr!("x")).is_ok());
        assert!(trace.get_f64_result(&addr!("missing")).is_err());
    }

    #[test]
    fn test_safe_vs_strict_handlers() {
        let mut rng = thread_rng();

        // Create base trace
        let mut base_trace = Trace::default();
        base_trace.insert_choice(addr!("param"), ChoiceValue::F64(2.5), -1.0);

        let model = sample(addr!("param"), Normal::new(0.0, 1.0).unwrap());

        // Safe replay should work
        let safe_handler = SafeReplayHandler {
            rng: &mut rng,
            base: base_trace,
            trace: Trace::default(),
            warn_on_mismatch: false,
        };
        let (result, trace) = runtime::handler::run(safe_handler, model);

        assert_eq!(result, 2.5);
        assert_eq!(trace.get_f64(&addr!("param")), Some(2.5));
    }

    #[test]
    fn test_model_structure_analysis() {
        let mut rng = thread_rng();

        let hierarchical_model = || {
            prob!(
                let global <- sample(addr!("global"), Normal::new(0.0, 1.0).unwrap());
                let locals <- plate!(i in 0..3 => {
                    sample(addr!("local", i), Normal::new(global, 0.1).unwrap())
                });
                pure((global, locals))
            )
        };

        let (_, trace) = runtime::handler::run(
            PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            },
            hierarchical_model(),
        );

        // Should have global + 3 local parameters
        assert_eq!(trace.choices.len(), 4);

        // Check address structure
        assert!(trace.choices.contains_key(&addr!("global")));
        assert!(trace.choices.contains_key(&addr!("local", 0)));
        assert!(trace.choices.contains_key(&addr!("local", 1)));
        assert!(trace.choices.contains_key(&addr!("local", 2)));
    }

    #[test]
    fn test_performance_diagnostics() {
        use std::time::Instant;
        let mut rng = thread_rng();

        let large_model = || {
            plate!(i in 0..50 => {
                sample(addr!("x", i), Normal::new(0.0, 1.0).unwrap())
            })
        };

        let start = Instant::now();
        let (_, trace) = runtime::handler::run(
            PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            },
            large_model(),
        );
        let duration = start.elapsed();

        assert_eq!(trace.choices.len(), 50);
        assert!(trace.total_log_weight().is_finite());

        // Performance should be reasonable
        assert!(duration.as_millis() < 1000, "Model execution too slow");
    }

Testing Strategy:

• Unit tests for individual model components
• Integration tests for complete workflows
• Performance regression tests
• Statistical validation against known results

Common Issues and Solutions

Issue: Infinite Log-Weights

Symptoms: trace.total_log_weight().is_infinite()

Causes:

• Factor statements with impossible constraints
• Parameters outside valid ranges
• Numerical overflow in likelihood computations

Solutions:

• Check factor conditions carefully
• Validate parameter ranges in constructors
• Use log-space computations for numerical stability

Issue: Missing or Wrong Parameter Values

Symptoms: get_*() returns None or wrong types

Causes:

• Address typos or inconsistencies
• Model structure doesn't match expectations
• Type mismatches in trace replay

Solutions:

• Use consistent address naming conventions
• Print all addresses for verification
• Use safe handlers for production resilience

Issue: Poor MCMC Convergence

Symptoms: High R-hat values, low ESS

Causes:

• Inappropriate step sizes
• Poor model parameterization
• Insufficient warm-up periods

Solutions:

• Increase warm-up iterations
• Reparameterize for better geometry
• Use adaptive algorithms with proper tuning

Issue: Slow Model Execution

Symptoms: High execution times, memory usage

Causes:

• Inefficient model structure
• Excessive address creation
• Large trace construction overhead

Solutions:

• Use plate! for vectorized operations
• Pre-allocate data structures when possible
• Profile with performance diagnostics

Best Practices Summary

1. Debug Incrementally: Start simple and add complexity systematically
2. Use All Tools: Combine trace inspection, validation, and diagnostics
3. Test Edge Cases: Verify behavior at parameter boundaries
4. Monitor Performance: Track execution time and memory usage
5. Validate Statistically: Compare against known theoretical results
6. Handle Errors Gracefully: Use safe handlers in production
7. Document Assumptions: Clear model specifications aid debugging

Effective debugging transforms probabilistic programming from guesswork into systematic model development. Fugue's comprehensive debugging toolkit enables confident deployment of complex probabilistic systems.

Custom Handlers

Fugue's handler system is grounded in algebraic effect theory, providing a principled approach to effect interpretation and computational extension. Custom handlers enable specialized execution strategies, monitoring systems, and novel inference algorithms through systematic effect handling and handler composition.

Understanding the Handler Trait

The Handler trait provides the algebraic signature for probabilistic effects. Each method represents an effect operation with its semantic interpretation:

graph TD
    subgraph "Handler Architecture"
        A[Effect E] --> B{Effect Type}
        B -->|sample| C[on_sample_T]
        B -->|observe| D[on_observe_T]  
        B -->|factor| E[on_factor]
        C --> F[Handler State H]
        D --> F
        E --> F
        F --> G[Updated State H']
        G --> H[Continue Execution]
    end

Effect Algebra: Each handler interprets the probabilistic effect signature:

$$\begin{array}{rl}\text{sample}& :\text{Dist}[T]\to T\\ \text{observe}& :\text{Dist}[T]\times T\to \text{Unit}\\ \text{factor}& :\mathbb{R}\to \text{Unit}\end{array}$$

where the carrier type varies by handler implementation.

/// Simple handler that just samples from priors (similar to PriorHandler)
struct BasicHandler<R: Rng> {
    rng: R,
    trace: Trace,
}

impl<R: Rng> Handler for BasicHandler<R> {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        let value = dist.sample(&mut self.rng);
        let log_prob = dist.log_prob(&value);

        // Store in trace
        self.trace.log_prior += log_prob;
        self.trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::F64(value),
                logp: log_prob,
            },
        );

        value
    }

    fn on_sample_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>) -> bool {
        let value = dist.sample(&mut self.rng);
        let log_prob = dist.log_prob(&value);

        self.trace.log_prior += log_prob;
        self.trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::Bool(value),
                logp: log_prob,
            },
        );

        value
    }

    fn on_sample_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>) -> u64 {
        let value = dist.sample(&mut self.rng);
        let log_prob = dist.log_prob(&value);

        self.trace.log_prior += log_prob;
        self.trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::U64(value),
                logp: log_prob,
            },
        );

        value
    }

    fn on_sample_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>) -> usize {
        let value = dist.sample(&mut self.rng);
        let log_prob = dist.log_prob(&value);

        self.trace.log_prior += log_prob;
        self.trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::Usize(value),
                logp: log_prob,
            },
        );

        value
    }

    fn on_observe_f64(&mut self, _addr: &Address, dist: &dyn Distribution<f64>, value: f64) {
        self.trace.log_likelihood += dist.log_prob(&value);
    }

    fn on_observe_bool(&mut self, _addr: &Address, dist: &dyn Distribution<bool>, value: bool) {
        self.trace.log_likelihood += dist.log_prob(&value);
    }

    fn on_observe_u64(&mut self, _addr: &Address, dist: &dyn Distribution<u64>, value: u64) {
        self.trace.log_likelihood += dist.log_prob(&value);
    }

    fn on_observe_usize(&mut self, _addr: &Address, dist: &dyn Distribution<usize>, value: usize) {
        self.trace.log_likelihood += dist.log_prob(&value);
    }

    fn on_factor(&mut self, log_weight: f64) {
        self.trace.log_factors += log_weight;
    }

    fn finish(self) -> Trace {
        self.trace
    }
}

Handler Responsibilities:

• Type-specific sampling: Handle f64, bool, u64, and usize distributions appropriately
• Observation handling: Process observed values and update likelihood components
• Factor management: Accumulate constraint and penalty terms
• Trace construction: Build execution traces with choices and log-weights
• Resource cleanup: Properly finalize and return traces

Decorator Pattern for Handler Composition

The decorator pattern implements handler composition through effect forwarding with computational augmentation. This pattern follows the mathematical principle of function composition:

$$(f\circ g)(x)=f(g(x))$$

Applied to handlers: $h_{\text{decorated}}=h_{\text{decorator}}\circ h_{\text{base}}$

graph LR
    subgraph "Handler Composition Chain"
        A[Effect] --> B[Decorator₁]
        B --> C[Decorator₂]
        C --> D[Base Handler]
        D --> E[Result]

        B -.->|"Log, Monitor"| F[Side Effects]
        C -.->|"Transform, Filter"| G[Modifications]
    end

Compositional Properties:

• Associativity: $(h_1\circ h_2)\circ h_3=h_1\circ (h_2\circ h_3)$
• Identity: $\text{id}\circ h=h\circ \text{id}=h$
• Effect Preservation: Core semantics remain unchanged

/// Handler decorator that logs all operations
struct LoggingHandler<H: Handler> {
    inner: H,
    log: Vec<String>,
    verbose: bool,
}

impl<H: Handler> LoggingHandler<H> {
    fn new(inner: H, verbose: bool) -> Self {
        Self {
            inner,
            log: Vec::new(),
            verbose,
        }
    }

    fn log_operation(&mut self, operation: String) {
        if self.verbose {
            println!("LOG: {}", operation);
        }
        self.log.push(operation);
    }
}

impl<H: Handler> Handler for LoggingHandler<H> {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        let value = self.inner.on_sample_f64(addr, dist);
        self.log_operation(format!("Sample f64 at {}: {:.3}", addr, value));
        value
    }

    fn on_sample_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>) -> bool {
        let value = self.inner.on_sample_bool(addr, dist);
        self.log_operation(format!("Sample bool at {}: {}", addr, value));
        value
    }

    fn on_sample_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>) -> u64 {
        let value = self.inner.on_sample_u64(addr, dist);
        self.log_operation(format!("Sample u64 at {}: {}", addr, value));
        value
    }

    fn on_sample_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>) -> usize {
        let value = self.inner.on_sample_usize(addr, dist);
        self.log_operation(format!("Sample usize at {}: {}", addr, value));
        value
    }

    fn on_observe_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>, value: f64) {
        self.log_operation(format!("Observe f64 at {}: {:.3}", addr, value));
        self.inner.on_observe_f64(addr, dist, value);
    }

    fn on_observe_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>, value: bool) {
        self.log_operation(format!("Observe bool at {}: {}", addr, value));
        self.inner.on_observe_bool(addr, dist, value);
    }

    fn on_observe_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>, value: u64) {
        self.log_operation(format!("Observe u64 at {}: {}", addr, value));
        self.inner.on_observe_u64(addr, dist, value);
    }

    fn on_observe_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>, value: usize) {
        self.log_operation(format!("Observe usize at {}: {}", addr, value));
        self.inner.on_observe_usize(addr, dist, value);
    }

    fn on_factor(&mut self, log_weight: f64) {
        self.log_operation(format!("Factor: {:.3}", log_weight));
        self.inner.on_factor(log_weight);
    }

    fn finish(self) -> Trace {
        let trace = self.inner.finish();
        println!("✅ Logged {} operations total", self.log.len());
        trace
    }
}

Decorator Benefits:

• Non-invasive functionality addition
• Composable and reusable components
• Separation of concerns between core logic and cross-cutting features
• Easy to enable/disable features dynamically

Stateful Handlers for Analytics

Handlers can maintain state to accumulate statistics and monitor model behavior:

/// Handler that accumulates statistics about model execution
#[derive(Debug)]
struct ExecutionStats {
    sample_counts: HashMap<String, u32>, // Type -> count
    observe_counts: HashMap<String, u32>,
    factor_count: u32,
    total_log_weight: f64,
    parameter_ranges: HashMap<String, (f64, f64)>, // Address -> (min, max) for f64 params
}

impl Default for ExecutionStats {
    fn default() -> Self {
        Self {
            sample_counts: HashMap::new(),
            observe_counts: HashMap::new(),
            factor_count: 0,
            total_log_weight: 0.0,
            parameter_ranges: HashMap::new(),
        }
    }
}

struct StatisticsHandler<H: Handler> {
    inner: H,
    stats: ExecutionStats,
}

impl<H: Handler> StatisticsHandler<H> {
    fn new(inner: H) -> Self {
        Self {
            inner,
            stats: ExecutionStats::default(),
        }
    }

    fn update_f64_range(&mut self, addr: &Address, value: f64) {
        let key = addr.0.clone();
        self.stats
            .parameter_ranges
            .entry(key)
            .and_modify(|(min, max)| {
                *min = min.min(value);
                *max = max.max(value);
            })
            .or_insert((value, value));
    }
}

impl<H: Handler> Handler for StatisticsHandler<H> {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        let value = self.inner.on_sample_f64(addr, dist);
        *self
            .stats
            .sample_counts
            .entry("f64".to_string())
            .or_insert(0) += 1;
        self.update_f64_range(addr, value);
        value
    }

    fn on_sample_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>) -> bool {
        let value = self.inner.on_sample_bool(addr, dist);
        *self
            .stats
            .sample_counts
            .entry("bool".to_string())
            .or_insert(0) += 1;
        value
    }

    fn on_sample_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>) -> u64 {
        let value = self.inner.on_sample_u64(addr, dist);
        *self
            .stats
            .sample_counts
            .entry("u64".to_string())
            .or_insert(0) += 1;
        value
    }

    fn on_sample_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>) -> usize {
        let value = self.inner.on_sample_usize(addr, dist);
        *self
            .stats
            .sample_counts
            .entry("usize".to_string())
            .or_insert(0) += 1;
        value
    }

    fn on_observe_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>, value: f64) {
        *self
            .stats
            .observe_counts
            .entry("f64".to_string())
            .or_insert(0) += 1;
        self.inner.on_observe_f64(addr, dist, value);
    }

    fn on_observe_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>, value: bool) {
        *self
            .stats
            .observe_counts
            .entry("bool".to_string())
            .or_insert(0) += 1;
        self.inner.on_observe_bool(addr, dist, value);
    }

    fn on_observe_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>, value: u64) {
        *self
            .stats
            .observe_counts
            .entry("u64".to_string())
            .or_insert(0) += 1;
        self.inner.on_observe_u64(addr, dist, value);
    }

    fn on_observe_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>, value: usize) {
        *self
            .stats
            .observe_counts
            .entry("usize".to_string())
            .or_insert(0) += 1;
        self.inner.on_observe_usize(addr, dist, value);
    }

    fn on_factor(&mut self, log_weight: f64) {
        self.stats.factor_count += 1;
        self.stats.total_log_weight += log_weight;
        self.inner.on_factor(log_weight);
    }

    fn finish(self) -> Trace {
        println!("✅ Execution Statistics:");
        println!("   - Samples by type: {:?}", self.stats.sample_counts);
        println!("   - Observations by type: {:?}", self.stats.observe_counts);
        println!("   - Factor operations: {}", self.stats.factor_count);
        println!("   - Parameter ranges:");
        for (addr, (min, max)) in &self.stats.parameter_ranges {
            println!("     {}: [{:.3}, {:.3}]", addr, min, max);
        }
        self.inner.finish()
    }
}

Analytics Applications:

• Model complexity analysis (parameter counts by type)
• Execution profiling and bottleneck identification
• Parameter range monitoring for numerical stability
• Distribution usage patterns for optimization

Conditional and Filtering Handlers

Implement business logic and constraints through conditional handling:

/// Handler that filters/modifies values based on conditions
struct FilteringHandler<H: Handler> {
    inner: H,
    f64_clamp_range: Option<(f64, f64)>,
    bool_flip_probability: f64,
    rng: rand::rngs::ThreadRng,
}

impl<H: Handler> FilteringHandler<H> {
    fn new(inner: H, f64_clamp_range: Option<(f64, f64)>, bool_flip_probability: f64) -> Self {
        Self {
            inner,
            f64_clamp_range,
            bool_flip_probability,
            rng: thread_rng(),
        }
    }
}

impl<H: Handler> Handler for FilteringHandler<H> {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        let mut value = self.inner.on_sample_f64(addr, dist);

        // Apply clamping if specified
        if let Some((min, max)) = self.f64_clamp_range {
            value = value.clamp(min, max);
        }

        value
    }

    fn on_sample_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>) -> bool {
        let mut value = self.inner.on_sample_bool(addr, dist);

        // Flip boolean with specified probability
        if self.rng.gen::<f64>() < self.bool_flip_probability {
            value = !value;
        }

        value
    }

    fn on_sample_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>) -> u64 {
        self.inner.on_sample_u64(addr, dist)
    }

    fn on_sample_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>) -> usize {
        self.inner.on_sample_usize(addr, dist)
    }

    fn on_observe_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>, value: f64) {
        self.inner.on_observe_f64(addr, dist, value);
    }

    fn on_observe_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>, value: bool) {
        self.inner.on_observe_bool(addr, dist, value);
    }

    fn on_observe_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>, value: u64) {
        self.inner.on_observe_u64(addr, dist, value);
    }

    fn on_observe_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>, value: usize) {
        self.inner.on_observe_usize(addr, dist, value);
    }

    fn on_factor(&mut self, log_weight: f64) {
        self.inner.on_factor(log_weight);
    }

    fn finish(self) -> Trace {
        self.inner.finish()
    }
}

Filtering Use Cases:

• Parameter clamping for numerical stability
• Outlier detection and handling
• Domain-specific constraints enforcement
• Robustness testing through perturbations

Performance Monitoring

Track and optimize computational characteristics with monitoring handlers:

use std::time::{Duration, Instant};

/// Handler that monitors performance characteristics
struct PerformanceHandler<H: Handler> {
    inner: H,
    start_time: Instant,
    operation_times: Vec<Duration>,
    sample_count: u32,
    observe_count: u32,
}

impl<H: Handler> PerformanceHandler<H> {
    fn new(inner: H) -> Self {
        Self {
            inner,
            start_time: Instant::now(),
            operation_times: Vec::new(),
            sample_count: 0,
            observe_count: 0,
        }
    }

    fn time_operation<F, R>(&mut self, operation: F) -> R
    where
        F: FnOnce(&mut H) -> R,
    {
        let start = Instant::now();
        let result = operation(&mut self.inner);
        let duration = start.elapsed();
        self.operation_times.push(duration);
        result
    }
}

impl<H: Handler> Handler for PerformanceHandler<H> {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        self.sample_count += 1;
        self.time_operation(|inner| inner.on_sample_f64(addr, dist))
    }

    fn on_sample_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>) -> bool {
        self.sample_count += 1;
        self.time_operation(|inner| inner.on_sample_bool(addr, dist))
    }

    fn on_sample_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>) -> u64 {
        self.sample_count += 1;
        self.time_operation(|inner| inner.on_sample_u64(addr, dist))
    }

    fn on_sample_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>) -> usize {
        self.sample_count += 1;
        self.time_operation(|inner| inner.on_sample_usize(addr, dist))
    }

    fn on_observe_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>, value: f64) {
        self.observe_count += 1;
        self.time_operation(|inner| inner.on_observe_f64(addr, dist, value))
    }

    fn on_observe_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>, value: bool) {
        self.observe_count += 1;
        self.time_operation(|inner| inner.on_observe_bool(addr, dist, value))
    }

    fn on_observe_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>, value: u64) {
        self.observe_count += 1;
        self.time_operation(|inner| inner.on_observe_u64(addr, dist, value))
    }

    fn on_observe_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>, value: usize) {
        self.observe_count += 1;
        self.time_operation(|inner| inner.on_observe_usize(addr, dist, value))
    }

    fn on_factor(&mut self, log_weight: f64) {
        self.time_operation(|inner| inner.on_factor(log_weight))
    }

    fn finish(self) -> Trace {
        let total_time = self.start_time.elapsed();
        let avg_op_time = if !self.operation_times.is_empty() {
            self.operation_times.iter().sum::<Duration>() / self.operation_times.len() as u32
        } else {
            Duration::ZERO
        };

        println!("✅ Performance Monitoring Results:");
        println!("   - Total execution time: {:?}", total_time);
        println!("   - Operations performed: {}", self.operation_times.len());
        println!("   - Sample operations: {}", self.sample_count);
        println!("   - Observe operations: {}", self.observe_count);
        println!("   - Average operation time: {:?}", avg_op_time);

        self.inner.finish()
    }
}

Performance Insights:

• Operation timing and bottleneck identification
• Memory allocation patterns
• Execution hotspots and optimization opportunities
• Scalability analysis for production deployment

Custom Inference Algorithms

Custom inference algorithms extend Fugue's effect interpretation to implement novel sampling strategies and approximate inference methods. Each algorithm provides a unique semantic mapping from probabilistic effects to computational actions:

graph TD
    subgraph "Inference Algorithm Architecture"
        A[Model M] --> B[Effect Sequence]
        B --> C{Handler Type}
        C -->|MCMC| D[Markov Chain<br/>Sampling]
        C -->|VI| E[Variational<br/>Approximation]
        C -->|SMC| F[Sequential<br/>Monte Carlo]
        C -->|ABC| G[Approximate<br/>Bayesian Computation]

        D --> H[Posterior Samples]
        E --> I[Approximate<br/>Distribution]
        F --> J[Weighted<br/>Particles]
        G --> K[Likelihood-Free<br/>Samples]
    end

Algorithm Design Principles:

1. Effect Consistency: $\forall e\in \mathcal{E}:h(e)$ preserves probabilistic semantics
2. Convergence Guarantees: Algorithm converges to target distribution under regularity conditions
3. Computational Tractability: Runtime complexity is polynomial in problem dimensions
4. Statistical Efficiency: Effective sample size scales appropriately with computational cost

Mathematical Framework: Each inference handler implements a stochastic operator $T:\mathcal{P}(\Theta )\to \mathcal{P}(\Theta )$ with fixed point $\pi$ such that $T\pi =\pi$.

/// Simple custom MCMC-like handler that perturbs existing values
struct SimpleMCMCHandler<R: Rng> {
    rng: R,
    base_trace: Trace,
    current_trace: Trace,
    perturbation_scale: f64,
}

impl<R: Rng> SimpleMCMCHandler<R> {
    fn new(rng: R, base_trace: Trace, perturbation_scale: f64) -> Self {
        Self {
            rng,
            base_trace,
            current_trace: Trace::default(),
            perturbation_scale,
        }
    }
}

impl<R: Rng> Handler for SimpleMCMCHandler<R> {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        let value = if let Some(base_value) = self.base_trace.get_f64(addr) {
            // Perturb existing value
            let perturbation = Normal::new(0.0, self.perturbation_scale).unwrap();
            base_value + perturbation.sample(&mut self.rng)
        } else {
            // Sample fresh if not in base trace
            dist.sample(&mut self.rng)
        };

        let log_prob = dist.log_prob(&value);
        self.current_trace.log_prior += log_prob;
        self.current_trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::F64(value),
                logp: log_prob,
            },
        );

        value
    }

    fn on_sample_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>) -> bool {
        let value = if let Some(base_value) = self.base_trace.get_bool(addr) {
            // Maybe flip the boolean with small probability
            if self.rng.gen::<f64>() < 0.1 {
                !base_value
            } else {
                base_value
            }
        } else {
            dist.sample(&mut self.rng)
        };

        let log_prob = dist.log_prob(&value);
        self.current_trace.log_prior += log_prob;
        self.current_trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::Bool(value),
                logp: log_prob,
            },
        );

        value
    }

    fn on_sample_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>) -> u64 {
        // For simplicity, just use base value or sample fresh
        let value = self
            .base_trace
            .get_u64(addr)
            .unwrap_or_else(|| dist.sample(&mut self.rng));

        let log_prob = dist.log_prob(&value);
        self.current_trace.log_prior += log_prob;
        self.current_trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::U64(value),
                logp: log_prob,
            },
        );

        value
    }

    fn on_sample_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>) -> usize {
        let value = self
            .base_trace
            .get_usize(addr)
            .unwrap_or_else(|| dist.sample(&mut self.rng));

        let log_prob = dist.log_prob(&value);
        self.current_trace.log_prior += log_prob;
        self.current_trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::Usize(value),
                logp: log_prob,
            },
        );

        value
    }

    fn on_observe_f64(&mut self, _addr: &Address, dist: &dyn Distribution<f64>, value: f64) {
        self.current_trace.log_likelihood += dist.log_prob(&value);
    }

    fn on_observe_bool(&mut self, _addr: &Address, dist: &dyn Distribution<bool>, value: bool) {
        self.current_trace.log_likelihood += dist.log_prob(&value);
    }

    fn on_observe_u64(&mut self, _addr: &Address, dist: &dyn Distribution<u64>, value: u64) {
        self.current_trace.log_likelihood += dist.log_prob(&value);
    }

    fn on_observe_usize(&mut self, _addr: &Address, dist: &dyn Distribution<usize>, value: usize) {
        self.current_trace.log_likelihood += dist.log_prob(&value);
    }

    fn on_factor(&mut self, log_weight: f64) {
        self.current_trace.log_factors += log_weight;
    }

    fn finish(self) -> Trace {
        self.current_trace
    }
}

Inference Handler Patterns:

• MCMC variants: Custom proposal mechanisms and acceptance criteria
• Variational methods: Gradient-based optimization with custom families
• Rejection sampling: Domain-specific acceptance/rejection logic
• Importance sampling: Custom proposal distributions and weight calculations

Handler Composition and Chaining

Handler chaining implements multi-stage effect processing through systematic composition operators. The composition forms a computational pipeline with well-defined data flow and effect propagation:

graph TD
    subgraph "Handler Composition Pipeline"
        A[Raw Effect E] --> B[Statistics Handler]
        B --> C[Logging Handler]
        C --> D[Performance Handler]
        D --> E[Base Handler]
        E --> F[Result + Trace]

        B -.->|Metrics| G[(Statistics DB)]
        C -.->|Events| H[(Log Stream)]
        D -.->|Timing| I[(Performance Monitor)]

        F --> J{Validation}
        J -->|Pass| K[Success]
        J -->|Fail| L[Error Recovery]
    end

Composition Laws:

1. Preservation: $h_n\circ \dots \circ h_1$ preserves effect semantics
2. Associativity: Composition order affects performance but not correctness
3. Commutativity: Decorators with disjoint side effects commute
4. Distributivity: $h\circ (g_1+g_2)=h\circ g_1+h\circ g_2$ for effect unions

Performance Analysis: Handler chain depth $d$ introduces overhead $\mathcal{O}(d\cdot c)$ where $c$ is the per-handler cost. Optimization strategies include handler fusion and effect batching.

fn main() {
    println!("=== Custom Handlers in Fugue ===\n");

    println!("1. Basic Custom Handler Implementation");
    println!("------------------------------------");

    // Test the basic handler
    let mut rng = thread_rng();
    let handler = BasicHandler {
        rng: &mut rng,
        trace: Trace::default(),
    };

    let test_model = || sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
    let (result, trace) = runtime::handler::run(handler, test_model());

    println!("✅ Basic handler executed");
    println!("   - Result: {:.3}", result);
    println!("   - Trace choices: {}", trace.choices.len());
    println!("   - Total log-weight: {:.3}", trace.total_log_weight());
    println!();

    println!("2. Logging Handler - Decorator Pattern");
    println!("-------------------------------------");

    // Test the logging handler
    let mut rng = thread_rng();
    let base_handler = PriorHandler {
        rng: &mut rng,
        trace: Trace::default(),
    };
    let logging_handler = LoggingHandler::new(base_handler, false); // Non-verbose

    let logged_model = || {
        prob!(
            let x <- sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap());
            observe(addr!("obs"), Normal::new(x, 0.5).unwrap(), 1.2);
            factor(-0.5);
            pure(x)
        )
    };

    let (result, _trace) = runtime::handler::run(logging_handler, logged_model());
    println!("   - Logged execution result: {:.3}", result);
    println!();

    println!("3. Statistics Accumulating Handler");
    println!("--------------------------------");

    // Test the statistics handler
    let mut rng = thread_rng();
    let base_handler = PriorHandler {
        rng: &mut rng,
        trace: Trace::default(),
    };
    let stats_handler = StatisticsHandler::new(base_handler);

    let complex_model = || {
        prob!(
            let mu <- sample(addr!("mu"), Normal::new(0.0, 2.0).unwrap());
            let is_outlier <- sample(addr!("outlier"), Bernoulli::new(0.1).unwrap());
            let count <- sample(addr!("count"), Poisson::new(3.0).unwrap());
            let category <- sample(addr!("category"), Categorical::new(vec![0.3, 0.4, 0.3]).unwrap());

            observe(addr!("y1"), Normal::new(mu, 1.0).unwrap(), 1.5);
            observe(addr!("y2"), Normal::new(mu, 1.0).unwrap(), 2.1);
            factor(if is_outlier { -2.0 } else { 0.0 });

            pure((mu, is_outlier, count, category))
        )
    };

    let (result, _trace) = runtime::handler::run(stats_handler, complex_model());
    println!(
        "   - Complex model result: {:?}",
        (result.0.round(), result.1, result.2, result.3)
    );
    println!();

    println!("4. Conditional Filtering Handler");
    println!("-------------------------------");

    // Test the filtering handler
    let mut rng = thread_rng();
    let base_handler = PriorHandler {
        rng: &mut rng,
        trace: Trace::default(),
    };
    let filtering_handler = FilteringHandler::new(
        base_handler,
        Some((-2.0, 2.0)), // Clamp f64 values to [-2, 2]
        0.1,               // 10% chance to flip booleans
    );

    let filter_test_model = || {
        prob!(
            let x <- sample(addr!("x"), Normal::new(0.0, 5.0).unwrap()); // Wide distribution
            let flag <- sample(addr!("flag"), Bernoulli::new(0.8).unwrap());
            pure((x, flag))
        )
    };

    let (result, _trace) = runtime::handler::run(filtering_handler, filter_test_model());
    println!("✅ Filtering handler executed");
    println!("   - Clamped value: {:.3} (should be in [-2, 2])", result.0);
    println!(
        "   - Boolean value: {} (may be flipped from original)",
        result.1
    );
    println!();

    println!("5. Performance Monitoring Handler");
    println!("--------------------------------");

    // Test the performance handler
    let mut rng = thread_rng();
    let base_handler = PriorHandler {
        rng: &mut rng,
        trace: Trace::default(),
    };
    let perf_handler = PerformanceHandler::new(base_handler);

    let perf_test_model = || {
        plate!(i in 0..10 => {
            sample(addr!("param", i), Normal::new(0.0, 1.0).unwrap())
        })
    };

    let (_result, _trace) = runtime::handler::run(perf_handler, perf_test_model());
    println!();

    println!("6. Custom Inference Handler");
    println!("---------------------------");

    // Test the custom inference handler
    let mut rng1 = thread_rng();
    let rng2 = thread_rng();

    // First get a base trace
    let base_handler = PriorHandler {
        rng: &mut rng1,
        trace: Trace::default(),
    };

    let inference_model = || {
        prob!(
            let mu <- sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap());
            observe(addr!("y"), Normal::new(mu, 0.5).unwrap(), 1.0);
            pure(mu)
        )
    };

    let (base_result, base_trace) = runtime::handler::run(base_handler, inference_model());

    // Now use custom MCMC handler to perturb it
    let base_log_weight = base_trace.total_log_weight();
    let mcmc_handler = SimpleMCMCHandler::new(rng2, base_trace, 0.1);
    let (mcmc_result, mcmc_trace) = runtime::handler::run(mcmc_handler, inference_model());

    println!("✅ Custom MCMC-like inference:");
    println!("   - Base result: {:.3}", base_result);
    println!("   - MCMC result: {:.3}", mcmc_result);
    println!("   - Base log-weight: {:.3}", base_log_weight);
    println!("   - MCMC log-weight: {:.3}", mcmc_trace.total_log_weight());
    println!();

    println!("7. Handler Composition and Chaining");
    println!("----------------------------------");

    // Demonstrate composing multiple handler decorators
    let mut rng = thread_rng();
    let base_handler = PriorHandler {
        rng: &mut rng,
        trace: Trace::default(),
    };

    // Chain multiple decorators: Statistics -> Logging -> Performance -> Base
    let stats_handler = StatisticsHandler::new(base_handler);
    let logging_handler = LoggingHandler::new(stats_handler, false);
    let performance_handler = PerformanceHandler::new(logging_handler);

    let composition_model = || {
        prob!(
            let x <- sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
            let y <- sample(addr!("y"), Bernoulli::new(0.7).unwrap());
            observe(addr!("obs"), Normal::new(x, 0.2).unwrap(), 0.5);
            factor(-0.3);
            pure((x, y))
        )
    };

    println!("✅ Handler composition example:");
    let (_result, _trace) = runtime::handler::run(performance_handler, composition_model());
    println!("   - Multiple handler layers executed successfully");
    println!();

    println!("=== Custom Handler Patterns Demonstrated! ===");
}

Composition Strategies:

• Layered approach: Statistics → Logging → Performance → Base
• Conditional activation: Enable decorators based on environment/configuration
• Feature flags: Runtime selection of handler combinations
• Pipeline optimization: Order decorators for minimal overhead

Advanced Handler Patterns

Caching Handler

struct CachingHandler<H: Handler> {
    inner: H,
    cache: HashMap<(Address, String), ChoiceValue>, // Address + dist info -> cached value
}

Distributed Handler

struct DistributedHandler<H: Handler> {
    inner: H,
    worker_id: usize,
    coordinator: Arc<Mutex<SharedState>>,
}

Fault-Tolerant Handler

struct FaultTolerantHandler<H: Handler> {
    inner: H,
    fallback_strategy: FallbackMode,
    error_count: u32,
    max_errors: u32,
}

Testing Custom Handlers

Systematic testing ensures handler correctness:

    #[test]
    fn test_basic_custom_handler() {
        let mut rng = thread_rng();
        let handler = BasicHandler {
            rng: &mut rng,
            trace: Trace::default(),
        };

        let model = sample(addr!("test"), Normal::new(0.0, 1.0).unwrap());
        let (result, trace) = runtime::handler::run(handler, model);

        assert!(trace.choices.contains_key(&addr!("test")));
        assert!(trace.total_log_weight().is_finite());
        assert!(result.is_finite());
    }

    #[test]
    fn test_logging_handler() {
        let mut rng = thread_rng();
        let base_handler = PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        };
        let logging_handler = LoggingHandler::new(base_handler, false);

        let model = prob!(
            let x <- sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
            observe(addr!("obs"), Normal::new(x, 0.1).unwrap(), 1.0);
            pure(x)
        );

        let (result, trace) = runtime::handler::run(logging_handler, model);

        assert!(trace.choices.contains_key(&addr!("x")));
        assert!(trace.log_likelihood.is_finite());
        assert!(result.is_finite());
    }

    #[test]
    fn test_statistics_handler() {
        let mut rng = thread_rng();
        let base_handler = PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        };
        let stats_handler = StatisticsHandler::new(base_handler);

        let model = prob!(
            let x <- sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
            let flag <- sample(addr!("flag"), Bernoulli::new(0.5).unwrap());
            pure((x, flag))
        );

        let (result, trace) = runtime::handler::run(stats_handler, model);

        assert_eq!(trace.choices.len(), 2);
        assert!(result.0.is_finite());
    }

    #[test]
    fn test_handler_composition() {
        let mut rng = thread_rng();
        let base_handler = PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        };

        // Compose multiple handlers
        let logged_handler = LoggingHandler::new(base_handler, false);
        let stats_handler = StatisticsHandler::new(logged_handler);

        let model = prob!(
            let x <- sample(addr!("param"), Normal::new(0.0, 1.0).unwrap());
            factor(-0.5);
            pure(x)
        );

        let (result, trace) = runtime::handler::run(stats_handler, model);

        assert!(trace.choices.contains_key(&addr!("param")));
        assert!(trace.log_factors.abs() > 0.0); // Factor was applied
        assert!(result.is_finite());
    }

Testing Strategy:

• Unit tests: Individual handler method behavior
• Integration tests: Handler with realistic models
• Property tests: Invariant verification across random inputs
• Composition tests: Multi-layer handler combinations

Production Considerations

Error Handling

impl Handler for ProductionHandler {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        match self.inner.on_sample_f64(addr, dist) {
            value if value.is_finite() => value,
            _ => {
                self.log_error(addr, "Non-finite sample");
                0.0 // Safe fallback
            }
        }
    }
}

Memory Management

struct MemoryEfficientHandler<H: Handler> {
    inner: H,
    choice_pool: Vec<Choice>, // Reusable allocations
    max_trace_size: usize,
}

Monitoring Integration

struct MetricsHandler<H: Handler> {
    inner: H,
    metrics_client: MetricsClient,
    model_name: String,
}

Common Patterns Summary

1. Decorator Pattern: Wrap handlers for additional functionality
2. State Accumulation: Track statistics and model behavior
3. Conditional Logic: Apply domain-specific constraints
4. Performance Monitoring: Identify bottlenecks and optimization opportunities
5. Custom Inference: Implement specialized algorithms
6. Composition: Chain multiple handlers for comprehensive capabilities
7. Error Handling: Graceful degradation and recovery
8. Resource Management: Efficient memory and computation usage

Best Practices

1. Single Responsibility: Each handler should have one clear purpose
2. Composability: Design handlers to work well in combination
3. Type Safety: Leverage Rust's type system for correctness
4. Performance: Minimize overhead in hot paths
5. Error Handling: Fail gracefully with meaningful diagnostics
6. Testing: Comprehensive unit and integration tests
7. Documentation: Clear API contracts and usage examples

Custom handlers transform Fugue from a probabilistic programming framework into a platform for building specialized inference systems, analytics tools, and production-ready probabilistic applications.

Production Deployment

Production deployment of probabilistic models requires reliability engineering, performance optimization, and operational excellence at scale. This guide establishes a mathematical framework for fault tolerance, service reliability, and system observability using Fugue's production-ready infrastructure patterns.

Error Handling and Graceful Degradation

Graceful degradation implements fault tolerance through systematic error recovery and service continuity. The mathematical foundation relies on Markov reliability models and circuit breaker theory:

stateDiagram-v2
    [*] --> Healthy
    Healthy --> Degraded : Error Rate > τ₁
    Degraded --> Failed : Error Rate > τ₂  
    Failed --> Recovery : Time > T_recovery
    Recovery --> Healthy : Success Rate > σ
    Degraded --> Healthy : Error Rate < τ₀
    
    note right of Healthy
        Error Rate: λ < τ₀
        SLA: 99.9%
        Full Functionality
    end note
    
    note right of Degraded  
        Error Rate: τ₀ < λ < τ₁
        SLA: 95%
        Limited Functionality
    end note
    
    note right of Failed
        Error Rate: λ > τ₂
        Circuit Open
        Fallback Mode
    end note

Circuit Breaker Mathematics: The failure rate follows a Poisson process with rate $\lambda (t)$. The circuit breaker transitions based on:

$$P(\text{trip})=1-e^{-{\int }_0^T\lambda (t)dt}$$

Error Budget Model: For SLA target $\alpha$, the error budget is: $$\text{Budget}(t)=(1-\alpha )\cdot t-{\int }_0^t{1}_{\text{error}}(s)ds$$

/// Production-ready handler that gracefully handles failures
struct RobustProductionHandler<H: Handler> {
    inner: H,
    error_count: u32,
    max_errors: u32,
    _fallback_values: HashMap<String, ChoiceValue>,
    circuit_breaker_open: bool,
}

impl<H: Handler> RobustProductionHandler<H> {
    fn new(inner: H, max_errors: u32) -> Self {
        let mut fallback_values = HashMap::new();
        fallback_values.insert("default_f64".to_string(), ChoiceValue::F64(0.0));
        fallback_values.insert("default_bool".to_string(), ChoiceValue::Bool(false));
        fallback_values.insert("default_u64".to_string(), ChoiceValue::U64(0));
        fallback_values.insert("default_usize".to_string(), ChoiceValue::Usize(0));

        Self {
            inner,
            error_count: 0,
            max_errors,
            _fallback_values: fallback_values,
            circuit_breaker_open: false,
        }
    }

    fn handle_error(&mut self, operation: &str, addr: &Address) -> bool {
        self.error_count += 1;
        eprintln!("PRODUCTION ERROR: {} failed at address {}", operation, addr);

        if self.error_count >= self.max_errors {
            self.circuit_breaker_open = true;
            eprintln!("CIRCUIT BREAKER: Too many errors, switching to fallback mode");
        }

        self.circuit_breaker_open
    }

    fn get_fallback_f64(&self, addr: &Address) -> f64 {
        // In production, this might come from a cache, configuration, or ML model
        match addr.0.as_str() {
            s if s.contains("temperature") => 20.0,
            s if s.contains("price") => 100.0,
            s if s.contains("probability") => 0.5,
            _ => 0.0,
        }
    }
}

impl<H: Handler> Handler for RobustProductionHandler<H> {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        if self.circuit_breaker_open {
            return self.get_fallback_f64(addr);
        }

        match std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
            self.inner.on_sample_f64(addr, dist)
        })) {
            Ok(value) if value.is_finite() => value,
            Ok(invalid_value) => {
                eprintln!("Invalid f64 sample: {} at {}", invalid_value, addr);
                self.handle_error("sample_f64", addr);
                self.get_fallback_f64(addr)
            }
            Err(_) => {
                self.handle_error("sample_f64_panic", addr);
                self.get_fallback_f64(addr)
            }
        }
    }

    fn on_sample_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>) -> bool {
        if self.circuit_breaker_open {
            return false; // Safe fallback
        }

        match std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
            self.inner.on_sample_bool(addr, dist)
        })) {
            Ok(value) => value,
            Err(_) => {
                self.handle_error("sample_bool_panic", addr);
                false
            }
        }
    }

    fn on_sample_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>) -> u64 {
        if self.circuit_breaker_open {
            return 1; // Safe default
        }

        match std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
            self.inner.on_sample_u64(addr, dist)
        })) {
            Ok(value) => value,
            Err(_) => {
                self.handle_error("sample_u64_panic", addr);
                1
            }
        }
    }

    fn on_sample_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>) -> usize {
        if self.circuit_breaker_open {
            return 0; // Safe array index
        }

        match std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
            self.inner.on_sample_usize(addr, dist)
        })) {
            Ok(value) => value,
            Err(_) => {
                self.handle_error("sample_usize_panic", addr);
                0
            }
        }
    }

    fn on_observe_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>, value: f64) {
        if !self.circuit_breaker_open
            && std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
                self.inner.on_observe_f64(addr, dist, value)
            }))
            .is_err()
        {
            self.handle_error("observe_f64_panic", addr);
        }
    }

    fn on_observe_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>, value: bool) {
        if !self.circuit_breaker_open
            && std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
                self.inner.on_observe_bool(addr, dist, value)
            }))
            .is_err()
        {
            self.handle_error("observe_bool_panic", addr);
        }
    }

    fn on_observe_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>, value: u64) {
        if !self.circuit_breaker_open
            && std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
                self.inner.on_observe_u64(addr, dist, value)
            }))
            .is_err()
        {
            self.handle_error("observe_u64_panic", addr);
        }
    }

    fn on_observe_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>, value: usize) {
        if !self.circuit_breaker_open
            && std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
                self.inner.on_observe_usize(addr, dist, value)
            }))
            .is_err()
        {
            self.handle_error("observe_usize_panic", addr);
        }
    }

    fn on_factor(&mut self, log_weight: f64) {
        if !log_weight.is_finite() {
            eprintln!("Invalid factor log-weight: {}", log_weight);
            self.error_count += 1;
            return; // Skip invalid factors
        }

        if !self.circuit_breaker_open {
            self.inner.on_factor(log_weight);
        }
    }

    fn finish(self) -> Trace {
        println!("✅ Production handler statistics:");
        println!(
            "   - Errors encountered: {}/{}",
            self.error_count, self.max_errors
        );
        println!(
            "   - Circuit breaker status: {}",
            if self.circuit_breaker_open {
                "OPEN (fallback mode)"
            } else {
                "CLOSED (normal)"
            }
        );

        if self.circuit_breaker_open {
            // Return minimal valid trace in fallback mode
            Trace::default()
        } else {
            self.inner.finish()
        }
    }
}

Robust Error Handling Features:

• Circuit Breaker Pattern: Prevents cascade failures by switching to fallback mode
• Panic Recovery: Catches panics and provides safe default values
• Input Validation: Ensures all inputs are finite and within expected ranges
• Fallback Values: Domain-specific defaults for different parameter types
• Error Counting: Tracks error rates to trigger circuit breaker activation

Configuration Management

Production models require flexible configuration for different environments:

#[derive(Debug, Clone)]
struct ModelConfig {
    // Model parameters
    temperature_prior_mean: f64,
    temperature_prior_std: f64,
    validity_probability: f64,
    sensor_noise_std: f64,

    // Runtime configuration
    max_inference_time_ms: u64,
    memory_pool_size: usize,
    enable_circuit_breaker: bool,
    error_threshold: u32,

    // Environment settings
    environment: String, // "development", "staging", "production"
    _log_level: String,
    enable_metrics: bool,
}

impl Default for ModelConfig {
    fn default() -> Self {
        Self {
            temperature_prior_mean: 20.0,
            temperature_prior_std: 5.0,
            validity_probability: 0.95,
            sensor_noise_std: 1.0,
            max_inference_time_ms: 1000,
            memory_pool_size: 100,
            enable_circuit_breaker: true,
            error_threshold: 10,
            environment: "production".to_string(),
            _log_level: "info".to_string(),
            enable_metrics: true,
        }
    }
}

struct ConfigurableModelRunner {
    config: ModelConfig,
    pool: TracePool,
    metrics: ProductionMetrics,
}

impl ConfigurableModelRunner {
    fn new(config: ModelConfig) -> Self {
        Self {
            pool: TracePool::new(config.memory_pool_size),
            metrics: ProductionMetrics::new(config.enable_metrics),
            config,
        }
    }

    fn create_model(&self) -> Model<(f64, bool)> {
        let config = self.config.clone();
        prob!(
            let temp <- sample(
                addr!("temperature"),
                Normal::new(config.temperature_prior_mean, config.temperature_prior_std).unwrap()
            );
            let valid <- sample(
                addr!("valid"),
                Bernoulli::new(config.validity_probability).unwrap()
            );
            // Simulate sensor reading with configured noise
            observe(
                addr!("sensor"),
                Normal::new(temp, config.sensor_noise_std).unwrap(),
                22.0
            );
            pure((temp, valid))
        )
    }

    fn run_inference(&mut self) -> Result<(f64, bool), String> {
        let start = Instant::now();

        // Configure handler based on environment
        let mut rng = thread_rng();
        let model = self.create_model(); // Create model before borrowing
        let result = if self.config.environment == "production" {
            // Use safe, fault-tolerant execution in production
            let base_handler = PooledPriorHandler::new(&mut rng, &mut self.pool);
            let robust_handler =
                RobustProductionHandler::new(base_handler, self.config.error_threshold);

            let (result, _trace) = runtime::handler::run(robust_handler, model);
            Ok(result)
        } else {
            // Use faster, less safe execution in development
            let handler = PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            };
            let (result, _trace) = runtime::handler::run(handler, model);
            Ok(result)
        };

        let duration = start.elapsed();
        self.metrics.record_inference_time(duration);

        // Check timeout
        if duration.as_millis() > self.config.max_inference_time_ms as u128 {
            self.metrics.increment_timeout_count();
            return Err(format!(
                "Inference timeout: {}ms > {}ms",
                duration.as_millis(),
                self.config.max_inference_time_ms
            ));
        }

        result
    }
}

Configuration Best Practices:

• Environment-Specific Settings: Different behavior for development/staging/production
• Model Parameter Configuration: Tunable priors, noise levels, and thresholds
• Runtime Configuration: Memory pool sizes, timeout limits, error thresholds
• Deployment Configuration: Circuit breaker settings, logging levels, metrics enablement
• Type-Safe Defaults: Sensible fallbacks for all configuration parameters

Production Metrics and Observability

Observability requires systematic metric collection with statistical analysis and anomaly detection. The metric taxonomy follows the USE method (Utilization, Saturation, Errors) and RED method (Rate, Errors, Duration):

graph TD
    subgraph "Observability Architecture"
        A[Model Execution] --> B[Metric Collection]
        B --> C{Metric Type}
        C -->|USE| D[Resource Metrics]
        C -->|RED| E[Service Metrics]  
        C -->|Business| F[Domain Metrics]
        
        D --> G[Utilization: ρ = λ/μ]
        D --> H[Saturation: Queue Length]
        D --> I[Error Rate: λₑ]
        
        E --> J[Request Rate: λᵣ]
        E --> K[Error Rate: εᵣ] 
        E --> L[Duration: T₉₉]
        
        F --> M[Inference Accuracy]
        F --> N[Model Drift]
        F --> O[Business KPIs]
        
        G --> P[(Time Series DB)]
        H --> P
        I --> P
        J --> P
        K --> P
        L --> P
        M --> Q[(Analytics DB)]
        N --> Q
        O --> Q
    end

Statistical Process Control: Metrics follow control chart theory with statistical control limits:

$$\begin{array}{rl}\text{UCL}& =\bar{X}+3\sigma /\sqrt{n}\\ \text{LCL}& =\bar{X}-3\sigma /\sqrt{n}\end{array}$$

Anomaly Detection: Using exponentially weighted moving averages: $${\text{EWMA}}_t=\alpha X_t+(1-\alpha ){\text{EWMA}}_{t-1}$$

#[derive(Debug, Clone)]
struct ProductionMetrics {
    enabled: bool,
    inference_count: u64,
    error_count: u64,
    timeout_count: u64,
    total_inference_time: Duration,
    start_time: SystemTime,
}

impl ProductionMetrics {
    fn new(enabled: bool) -> Self {
        Self {
            enabled,
            inference_count: 0,
            error_count: 0,
            timeout_count: 0,
            total_inference_time: Duration::ZERO,
            start_time: SystemTime::now(),
        }
    }

    fn record_inference_time(&mut self, duration: Duration) {
        if self.enabled {
            self.inference_count += 1;
            self.total_inference_time += duration;
        }
    }

    fn _increment_error_count(&mut self) {
        if self.enabled {
            self.error_count += 1;
        }
    }

    fn increment_timeout_count(&mut self) {
        if self.enabled {
            self.timeout_count += 1;
        }
    }

    fn get_stats(&self) -> HashMap<String, f64> {
        let mut stats = HashMap::new();
        if self.enabled {
            let uptime = self.start_time.elapsed().unwrap_or(Duration::ZERO);
            let avg_inference_time = if self.inference_count > 0 {
                self.total_inference_time.as_millis() as f64 / self.inference_count as f64
            } else {
                0.0
            };

            stats.insert("inference_count".to_string(), self.inference_count as f64);
            stats.insert("error_count".to_string(), self.error_count as f64);
            stats.insert("timeout_count".to_string(), self.timeout_count as f64);
            stats.insert(
                "error_rate".to_string(),
                if self.inference_count > 0 {
                    self.error_count as f64 / self.inference_count as f64
                } else {
                    0.0
                },
            );
            stats.insert("avg_inference_time_ms".to_string(), avg_inference_time);
            stats.insert("uptime_seconds".to_string(), uptime.as_secs() as f64);
            stats.insert(
                "throughput_per_second".to_string(),
                if uptime.as_secs() > 0 {
                    self.inference_count as f64 / uptime.as_secs() as f64
                } else {
                    0.0
                },
            );
        }
        stats
    }

    fn export_prometheus_metrics(&self) -> String {
        let mut metrics = String::new();
        let stats = self.get_stats();

        metrics.push_str("# HELP fugue_inference_total Total number of inference runs\n");
        metrics.push_str("# TYPE fugue_inference_total counter\n");
        metrics.push_str(&format!(
            "fugue_inference_total {}\n",
            stats.get("inference_count").unwrap_or(&0.0)
        ));

        metrics.push_str("# HELP fugue_errors_total Total number of errors\n");
        metrics.push_str("# TYPE fugue_errors_total counter\n");
        metrics.push_str(&format!(
            "fugue_errors_total {}\n",
            stats.get("error_count").unwrap_or(&0.0)
        ));

        metrics.push_str(
            "# HELP fugue_inference_duration_ms Average inference duration in milliseconds\n",
        );
        metrics.push_str("# TYPE fugue_inference_duration_ms gauge\n");
        metrics.push_str(&format!(
            "fugue_inference_duration_ms {}\n",
            stats.get("avg_inference_time_ms").unwrap_or(&0.0)
        ));

        metrics.push_str("# HELP fugue_error_rate Error rate (errors/total inferences)\n");
        metrics.push_str("# TYPE fugue_error_rate gauge\n");
        metrics.push_str(&format!(
            "fugue_error_rate {}\n",
            stats.get("error_rate").unwrap_or(&0.0)
        ));

        metrics
    }
}

/// Production monitoring handler that integrates with metrics systems
struct MetricsHandler<H: Handler> {
    inner: H,
    metrics: Arc<std::sync::Mutex<ProductionMetrics>>,
    _model_name: String,
}

impl<H: Handler> MetricsHandler<H> {
    fn new(
        inner: H,
        metrics: Arc<std::sync::Mutex<ProductionMetrics>>,
        model_name: String,
    ) -> Self {
        Self {
            inner,
            metrics,
            _model_name: model_name,
        }
    }
}

impl<H: Handler> Handler for MetricsHandler<H> {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        let start = Instant::now();
        let result = self.inner.on_sample_f64(addr, dist);

        if let Ok(mut metrics) = self.metrics.lock() {
            metrics.record_inference_time(start.elapsed());
        }

        result
    }

    fn on_sample_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>) -> bool {
        self.inner.on_sample_bool(addr, dist)
    }

    fn on_sample_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>) -> u64 {
        self.inner.on_sample_u64(addr, dist)
    }

    fn on_sample_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>) -> usize {
        self.inner.on_sample_usize(addr, dist)
    }

    fn on_observe_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>, value: f64) {
        self.inner.on_observe_f64(addr, dist, value);
    }

    fn on_observe_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>, value: bool) {
        self.inner.on_observe_bool(addr, dist, value);
    }

    fn on_observe_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>, value: u64) {
        self.inner.on_observe_u64(addr, dist, value);
    }

    fn on_observe_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>, value: usize) {
        self.inner.on_observe_usize(addr, dist, value);
    }

    fn on_factor(&mut self, log_weight: f64) {
        self.inner.on_factor(log_weight);
    }

    fn finish(self) -> Trace {
        self.inner.finish()
    }
}

Metrics Collection:

• Performance Metrics: Inference time, throughput, operation counts
• Error Tracking: Error rates, timeout counts, failure categorization
• System Health: Uptime, resource utilization, memory pool efficiency
• Prometheus Integration: Standard metrics format for monitoring systems
• Real-Time Dashboards: Live performance and health indicators

Health Checks and System Validation

Health monitoring implements continuous system validation through multi-level health checks with statistical thresholds and predictive alerting:

graph TD
    subgraph "Health Check Hierarchy"
        A[System Health H⁽ˢ⁾] --> B[Model Health H⁽ᵐ⁾]
        B --> C[Inference Health H⁽ⁱ⁾]
        C --> D[Resource Health H⁽ʳ⁾]

        A --> E{H⁽ˢ⁾ > θₛ?}
        B --> F{H⁽ᵐ⁾ > θₘ?}
        C --> G{H⁽ⁱ⁾ > θᵢ?}
        D --> H{H⁽ʳ⁾ > θʳ?}

        E -->|No| I[System Alert]
        F -->|No| J[Model Alert]
        G -->|No| K[Inference Alert]
        H -->|No| L[Resource Alert]

        E -->|Yes| M[System OK]
        F -->|Yes| M
        G -->|Yes| M
        H -->|Yes| M
    end

Health Score Calculation: Weighted combination of subsystem health: $$H^{(s)}=\sum _i{w}_i{H}^{(i)}$$

where $w_i$ are importance weights and $\sum w_i=1$.

Predictive Health Modeling: Using time series forecasting: $$H_{t+k}=\alpha H_t+\beta \frac{dH}{dt}{|}_t+\gamma \frac{d^{2}H}{d{t}^2}{|}_t$$

#[derive(Debug, Clone)]
struct HealthCheckResult {
    status: HealthStatus,
    message: String,
    details: HashMap<String, String>,
    timestamp: SystemTime,
}

#[derive(Debug, Clone, PartialEq)]
enum HealthStatus {
    Healthy,
    Degraded,
    Unhealthy,
}

struct ProductionHealthChecker {
    model_config: ModelConfig,
    metrics: Arc<std::sync::Mutex<ProductionMetrics>>,
}

impl ProductionHealthChecker {
    fn new(config: ModelConfig, metrics: Arc<std::sync::Mutex<ProductionMetrics>>) -> Self {
        Self {
            model_config: config,
            metrics,
        }
    }

    fn run_health_check(&self) -> HealthCheckResult {
        let mut details = HashMap::new();
        let mut overall_status = HealthStatus::Healthy;
        let mut messages = Vec::new();

        // Check 1: Model execution health
        match self.check_model_execution() {
            Ok(duration) => {
                details.insert("model_execution".to_string(), "healthy".to_string());
                details.insert(
                    "execution_time_ms".to_string(),
                    format!("{:.1}", duration.as_millis()),
                );
            }
            Err(e) => {
                overall_status = HealthStatus::Unhealthy;
                messages.push(format!("Model execution failed: {}", e));
                details.insert("model_execution".to_string(), "failed".to_string());
            }
        }

        // Check 2: Memory usage
        if let Some(pool_stats) = self.check_memory_health() {
            let hit_ratio = pool_stats.hit_ratio();
            details.insert("memory_hit_ratio".to_string(), format!("{:.2}%", hit_ratio));

            if hit_ratio < 50.0 {
                overall_status = HealthStatus::Degraded;
                messages.push("Low memory pool hit ratio".to_string());
            }
        }

        // Check 3: Error rates
        if let Ok(metrics) = self.metrics.lock() {
            let stats = metrics.get_stats();
            let error_rate = stats.get("error_rate").unwrap_or(&0.0) * 100.0;
            details.insert(
                "error_rate_percent".to_string(),
                format!("{:.2}%", error_rate),
            );

            if error_rate > 5.0 {
                overall_status = HealthStatus::Degraded;
                messages.push(format!("High error rate: {:.1}%", error_rate));
            } else if error_rate > 20.0 {
                overall_status = HealthStatus::Unhealthy;
                messages.push(format!("Critical error rate: {:.1}%", error_rate));
            }

            let avg_time = stats.get("avg_inference_time_ms").unwrap_or(&0.0);
            details.insert(
                "avg_inference_time_ms".to_string(),
                format!("{:.1}", avg_time),
            );

            if *avg_time > self.model_config.max_inference_time_ms as f64 * 0.8 {
                overall_status = HealthStatus::Degraded;
                messages.push("Inference time approaching timeout threshold".to_string());
            }
        }

        // Check 4: System resources
        details.insert(
            "memory_pool_size".to_string(),
            self.model_config.memory_pool_size.to_string(),
        );
        details.insert(
            "circuit_breaker".to_string(),
            if self.model_config.enable_circuit_breaker {
                "enabled".to_string()
            } else {
                "disabled".to_string()
            },
        );

        let message = if messages.is_empty() {
            "All systems healthy".to_string()
        } else {
            messages.join("; ")
        };

        HealthCheckResult {
            status: overall_status,
            message,
            details,
            timestamp: SystemTime::now(),
        }
    }

    fn check_model_execution(&self) -> Result<Duration, String> {
        let start = Instant::now();
        let mut rng = thread_rng();

        // Run a simplified version of the model for health checking
        let health_model = || {
            prob!(
                let value <- sample(addr!("health_check"), Normal::new(0.0, 1.0).unwrap());
                pure(value)
            )
        };

        let handler = PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        };

        match std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
            runtime::handler::run(handler, health_model())
        })) {
            Ok((result, trace)) => {
                if result.is_finite() && trace.total_log_weight().is_finite() {
                    Ok(start.elapsed())
                } else {
                    Err("Invalid model output".to_string())
                }
            }
            Err(_) => Err("Model execution panicked".to_string()),
        }
    }

    fn check_memory_health(&self) -> Option<fugue::runtime::memory::PoolStats> {
        // In a real implementation, this would check the actual memory pool
        // For demonstration, we'll create a temporary pool
        let pool = TracePool::new(10);
        Some(pool.stats().clone())
    }
}

Health Check Components:

• Model Execution Health: Verifies core functionality with simplified tests
• Memory Health: Monitors pool efficiency and memory usage patterns
• Error Rate Analysis: Tracks and categorizes different failure modes
• Performance Monitoring: Identifies degradation before it impacts users
• Multi-Level Status: Healthy/Degraded/Unhealthy with detailed diagnostics

Input Validation and Security

Robust input validation prevents security vulnerabilities and system failures:

/// Secure input validator for production model parameters
struct InputValidator;

impl InputValidator {
    fn validate_temperature(temp: f64) -> Result<f64, String> {
        match temp {
            t if !t.is_finite() => Err("Temperature must be finite".to_string()),
            t if t < -50.0 => Err("Temperature too low (< -50°C)".to_string()),
            t if t > 100.0 => Err("Temperature too high (> 100°C)".to_string()),
            t => Ok(t),
        }
    }

    fn validate_probability(p: f64) -> Result<f64, String> {
        match p {
            p if !p.is_finite() => Err("Probability must be finite".to_string()),
            p if p < 0.0 => Err("Probability must be non-negative".to_string()),
            p if p > 1.0 => Err("Probability must not exceed 1.0".to_string()),
            p => Ok(p),
        }
    }

    fn _validate_sensor_reading(reading: f64) -> Result<f64, String> {
        match reading {
            r if !r.is_finite() => Err("Sensor reading must be finite".to_string()),
            r if r.abs() > 1000.0 => Err("Sensor reading out of reasonable range".to_string()),
            r => Ok(r),
        }
    }

    fn sanitize_address_component(component: &str) -> Result<String, String> {
        // Prevent injection attacks in address components
        if component
            .chars()
            .any(|c| !(c.is_alphanumeric() || c == '_' || c == '-'))
        {
            return Err("Address component contains invalid characters".to_string());
        }

        if component.len() > 50 {
            Err("Address component too long".to_string())
        } else if component.is_empty() {
            Err("Address component cannot be empty".to_string())
        } else {
            Ok(component.to_string())
        }
    }
}

/// Production model with comprehensive input validation
fn create_validated_model(
    temperature_reading: f64,
    sensor_id: &str,
    prior_prob: f64,
) -> Result<Model<(f64, bool)>, String> {
    // Validate all inputs before model creation
    let validated_temp = InputValidator::validate_temperature(temperature_reading)?;
    let validated_prob = InputValidator::validate_probability(prior_prob)?;
    let sanitized_sensor_id = InputValidator::sanitize_address_component(sensor_id)?;

    // Additional business logic validation
    if sanitized_sensor_id.starts_with("test_") && validated_prob > 0.5 {
        return Err("Test sensors cannot have high prior probability".to_string());
    }

    Ok(prob!(
        let true_temp <- sample(
            addr!("temperature"),
            Normal::new(validated_temp, 2.0).unwrap()
        );
        let is_working <- sample(
            addr!("sensor_working", sanitized_sensor_id.clone()),
            Bernoulli::new(validated_prob).unwrap()
        );

        // Safe observation with validated input
        observe(
            addr!("reading", sanitized_sensor_id),
            Normal::new(true_temp, if is_working { 0.5 } else { 5.0 }).unwrap(),
            validated_temp
        );

        pure((true_temp, is_working))
    ))
}

Security Measures:

• Range Validation: Ensure parameters are within physically meaningful bounds
• Type Safety: Validate all inputs before model construction
• Sanitization: Clean address components to prevent injection attacks
• Business Rule Enforcement: Domain-specific validation logic
• Error Messages: Informative feedback without revealing system internals

Deployment Strategies and Patterns

Deployment strategies implement risk management through controlled rollout and statistical validation. Each strategy provides different risk-latency tradeoffs:

graph TD
    subgraph "Deployment Strategy Matrix"  
        A[New Model Version] --> B{Strategy Selection}
        B -->|Low Risk| C[Blue-Green]
        B -->|Medium Risk| D[Canary]
        B -->|High Risk| E[A/B Test]

        C --> F[Instant Switch<br/>Risk: High<br/>Rollback: Fast]
        D --> G[Gradual Rollout<br/>Risk: Medium<br/>Validation: Statistical]
        E --> H[Statistical Test<br/>Risk: Low<br/>Duration: Long]

        F --> I{Success?}
        G --> J{Performance > Baseline?}
        E --> K{Significance Test?}

        I -->|No| L[Instant Rollback]
        J -->|No| M[Gradual Rollback]
        K -->|No| N[Maintain Status Quo]

        I -->|Yes| O[Full Deployment]
        J -->|Yes| P[Continue Rollout]
        K -->|Yes| Q[Gradual Migration]
    end

Canary Analysis: Statistical significance testing for canary deployments:

$$\text{Z-score}=\frac{({\bar{X}}_{\text{canary}}-{\bar{X}}_{\text{control}})}{\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}}}$$

A/B Testing: Welch's t-test for unequal variances: $$t=\frac{{\bar{X}}_A-{\bar{X}}_B}{\sqrt{\frac{s_A^2}{n_A}+\frac{s_B^2}{n_B}}}$$

/// Production deployment manager with different strategies
#[derive(Debug, Clone)]
enum DeploymentStrategy {
    BlueGreen,
    CanaryRelease { percentage: f64 },
    RollingUpdate,
    _ImmediateSwitch,
}

struct ModelDeploymentManager {
    current_model_version: String,
    candidate_model_version: String,
    deployment_strategy: DeploymentStrategy,
    _rollback_threshold_error_rate: f64,
}

impl ModelDeploymentManager {
    fn new(strategy: DeploymentStrategy) -> Self {
        Self {
            current_model_version: "v1.0.0".to_string(),
            candidate_model_version: "v1.1.0".to_string(),
            deployment_strategy: strategy,
            _rollback_threshold_error_rate: 0.05, // 5% error rate triggers rollback
        }
    }

    fn should_use_candidate_model(&self, request_id: u64) -> bool {
        match &self.deployment_strategy {
            DeploymentStrategy::BlueGreen => {
                // In blue-green, we typically switch all traffic at once
                // For demo, we'll use request ID to simulate the switch
                request_id % 100 < 10 // 10% to candidate for testing
            }
            DeploymentStrategy::CanaryRelease { percentage } => {
                let hash = request_id % 100;
                (hash as f64) < (*percentage * 100.0)
            }
            DeploymentStrategy::RollingUpdate => {
                // Gradual rollout based on some criteria
                request_id % 10 < 3 // 30% rollout
            }
            DeploymentStrategy::_ImmediateSwitch => true,
        }
    }

    fn create_model(&self, use_candidate: bool) -> impl Fn() -> Model<f64> {
        let version = if use_candidate {
            self.candidate_model_version.clone()
        } else {
            self.current_model_version.clone()
        };

        move || {
            if version.starts_with("v1.1") {
                // Candidate model with improved parameters
                prob!(
                    let value <- sample(addr!("improved_param"), Normal::new(0.0, 0.8).unwrap());
                    factor(0.1); // Slight preference for this model
                    pure(value)
                )
            } else {
                // Current stable model
                prob!(
                    let value <- sample(addr!("stable_param"), Normal::new(0.0, 1.0).unwrap());
                    pure(value)
                )
            }
        }
    }

    fn process_request(&self, request_id: u64) -> Result<(f64, String), String> {
        let use_candidate = self.should_use_candidate_model(request_id);
        let version = if use_candidate {
            &self.candidate_model_version
        } else {
            &self.current_model_version
        };

        let model = self.create_model(use_candidate);

        // Execute with error handling
        let mut rng = thread_rng();
        let handler = PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        };

        match std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
            runtime::handler::run(handler, model())
        })) {
            Ok((result, trace)) => {
                if result.is_finite() && trace.total_log_weight().is_finite() {
                    Ok((result, version.clone()))
                } else {
                    Err(format!("Invalid result from model {}", version))
                }
            }
            Err(_) => Err(format!("Model {} panicked", version)),
        }
    }
}

Deployment Patterns:

• Blue-Green Deployment: Instant traffic switching between model versions
• Canary Releases: Gradual rollout to percentage of traffic for risk mitigation
• Rolling Updates: Progressive deployment across infrastructure
• A/B Testing: Compare model performance with statistical significance
• Rollback Capability: Quick reversion to previous version on issues

Performance Optimization Patterns

Memory Management

use fugue::runtime::memory::{TracePool, PooledPriorHandler};

// Production memory management
let mut pool = TracePool::new(1000);
let handler = PooledPriorHandler::new(&mut rng, &mut pool);

Batch Processing

// Process multiple inference requests efficiently
struct BatchProcessor {
    pool: TracePool,
    batch_size: usize,
}

impl BatchProcessor {
    fn process_batch(&mut self, requests: Vec<InferenceRequest>) -> Vec<InferenceResult> {
        requests.into_iter().map(|req| {
            let handler = PooledPriorHandler::new(&mut req.rng, &mut self.pool);
            self.run_single_inference(handler, req.model)
        }).collect()
    }
}

Connection Pooling

// Database connection management for model parameters
struct ModelParameterStore {
    connection_pool: Arc<ConnectionPool>,
    parameter_cache: LruCache<String, ModelParameters>,
}

Monitoring Integration

Prometheus Metrics

// Export metrics in Prometheus format
fn export_metrics(metrics: &ProductionMetrics) -> String {
    format!(
        "# HELP fugue_inference_total Total inference operations\n\
         TYPE fugue_inference_total counter\n\
         fugue_inference_total {}\n\
         HELP fugue_error_rate Current error rate\n\
         TYPE fugue_error_rate gauge\n\
         fugue_error_rate {}\n",
        metrics.inference_count,
        metrics.error_rate()
    )
}

Structured Logging

use serde_json::json;

// Structured logging for production debugging
fn log_inference_event(
    request_id: &str,
    model_version: &str,
    duration: Duration,
    result: &InferenceResult
) {
    let log_entry = json!({
        "event": "inference_completed",
        "request_id": request_id,
        "model_version": model_version,
        "duration_ms": duration.as_millis(),
        "success": result.is_success(),
        "timestamp": SystemTime::now(),
    });
    println!("{}", log_entry);
}

Alert Rules

// Define alerting thresholds
struct AlertRules {
    max_error_rate: f64,
    max_latency_ms: u64,
    min_throughput_per_sec: f64,
}

impl AlertRules {
    fn check_alerts(&self, metrics: &ProductionMetrics) -> Vec<Alert> {
        let mut alerts = Vec::new();

        if metrics.error_rate() > self.max_error_rate {
            alerts.push(Alert::HighErrorRate(metrics.error_rate()));
        }

        if metrics.avg_latency().as_millis() > self.max_latency_ms as u128 {
            alerts.push(Alert::HighLatency(metrics.avg_latency()));
        }

        alerts
    }
}

Testing in Production

Shadow Mode Testing

// Run new model versions in shadow mode
struct ShadowTester {
    primary_model: Box<dyn Fn() -> Model<f64>>,
    shadow_model: Box<dyn Fn() -> Model<f64>>,
    comparison_rate: f64,
}

impl ShadowTester {
    fn run_with_shadow(&mut self, input: &Input) -> (PrimaryResult, Option<ShadowResult>) {
        let primary = self.run_primary(input);

        let shadow = if rand::random::<f64>() < self.comparison_rate {
            Some(self.run_shadow(input))
        } else {
            None
        };

        (primary, shadow)
    }
}

Production Validation

// Continuous validation in production
fn validate_model_assumptions(trace: &Trace) -> ValidationResult {
    let mut issues = Vec::new();

    // Check log-weight stability
    if !trace.total_log_weight().is_finite() {
        issues.push("Non-finite log-weight detected".to_string());
    }

    // Check parameter ranges
    for (addr, choice) in &trace.choices {
        if let ChoiceValue::F64(value) = choice.value {
            if value.abs() > 1000.0 {
                issues.push(format!("Extreme value at {}: {}", addr, value));
            }
        }
    }

    ValidationResult { issues }
}

Operational Excellence

Infrastructure as Code

# Kubernetes deployment example
apiVersion: apps/v1
kind: Deployment
metadata:
  name: fugue-inference-service
spec:
  replicas: 3
  selector:
    matchLabels:
      app: fugue-inference
  template:
    metadata:
      labels:
        app: fugue-inference
    spec:
      containers:
      - name: inference-service
        image: fugue-inference:v1.2.0
        resources:
          requests:
            memory: "256Mi"
            cpu: "250m"
          limits:
            memory: "512Mi"
            cpu: "500m"
        livenessProbe:
          httpGet:
            path: /health
            port: 8080
          initialDelaySeconds: 30
          periodSeconds: 10

Service Level Objectives (SLOs)

// Define and monitor SLOs
struct ServiceLevelObjectives {
    availability_target: f64,    // 99.9%
    latency_p99_ms: u64,        // 100ms
    error_rate_threshold: f64,   // 0.1%
}

impl ServiceLevelObjectives {
    fn evaluate_slo_compliance(&self, metrics: &ProductionMetrics) -> SLOReport {
        SLOReport {
            availability: self.calculate_availability(metrics),
            latency_compliance: metrics.p99_latency() <= Duration::from_millis(self.latency_p99_ms),
            error_rate_compliance: metrics.error_rate() <= self.error_rate_threshold,
        }
    }
}

Capacity Planning

// Capacity planning and auto-scaling
struct CapacityPlanner {
    target_cpu_utilization: f64,
    target_memory_utilization: f64,
    scale_up_threshold: f64,
    scale_down_threshold: f64,
}

impl CapacityPlanner {
    fn recommend_scaling(&self, current_metrics: &SystemMetrics) -> ScalingRecommendation {
        if current_metrics.cpu_utilization > self.scale_up_threshold {
            ScalingRecommendation::ScaleUp(self.calculate_scale_factor(current_metrics))
        } else if current_metrics.cpu_utilization < self.scale_down_threshold {
            ScalingRecommendation::ScaleDown(0.5)
        } else {
            ScalingRecommendation::NoAction
        }
    }
}

Security Best Practices

Input Sanitization

Always validate and sanitize inputs before processing:

• Range checks for numerical parameters
• Character filtering for string inputs
• Business rule validation for domain constraints
• Rate limiting to prevent abuse
• Authentication and authorization for API access

Secret Management

// Secure configuration management
struct SecureConfig {
    database_url: SecretString,
    api_key: SecretString,
    model_parameters: ModelConfig,
}

impl SecureConfig {
    fn from_environment() -> Result<Self, ConfigError> {
        Ok(SecureConfig {
            database_url: env::var("DATABASE_URL")?.into(),
            api_key: env::var("API_KEY")?.into(),
            model_parameters: ModelConfig::from_file("model_config.toml")?,
        })
    }
}

Audit Logging

// Comprehensive audit trail
fn log_inference_request(
    user_id: &str,
    request: &InferenceRequest,
    response: &InferenceResponse,
) {
    let audit_log = AuditLogEntry {
        timestamp: SystemTime::now(),
        user_id: user_id.to_string(),
        action: "inference_request".to_string(),
        input_hash: hash_sensitive_data(&request.input),
        output_hash: hash_sensitive_data(&response.output),
        model_version: response.model_version.clone(),
        success: response.success,
    };

    audit_logger::log(audit_log);
}

Common Production Pitfalls

Memory Leaks

// Avoid: Creating new pools repeatedly
// for _ in 0..1000 {
//     let pool = TracePool::new(100); // Memory leak!
// }

// Do: Reuse pools across requests
let mut pool = TracePool::new(100);
for request in requests {
    let handler = PooledPriorHandler::new(&mut request.rng, &mut pool);
    process_request(handler, request);
}

Blocking Operations

// Avoid: Synchronous database calls in request handlers
// let result = database.query_sync(query); // Blocks event loop

// Do: Use async operations with proper timeouts
async fn process_request(request: Request) -> Result<Response, Error> {
    let timeout = Duration::from_millis(100);
    let result = tokio::time::timeout(timeout, database.query(query)).await??;
    Ok(result)
}

Error Propagation

// Avoid: Panicking on errors
// let value = risky_operation().unwrap(); // May crash service

// Do: Graceful error handling with fallbacks
let value = match risky_operation() {
    Ok(v) => v,
    Err(e) => {
        metrics.increment_error_count();
        log::warn!("Operation failed: {}, using fallback", e);
        fallback_value()
    }
};

Production deployment represents the culmination of probabilistic programming excellence, where theoretical foundations meet operational reality. Fugue's comprehensive tooling transforms academic probabilistic models into production-grade systems that deliver reliable, scalable, and maintainable probabilistic computing solutions.

Tutorials

Foundation Tutorials

Welcome to Fugue's Foundation Tutorials — your comprehensive introduction to the core concepts and unique features that make Fugue a revolutionary approach to probabilistic programming.

What You'll Learn

These tutorials build upon each other to give you a complete understanding of Fugue's foundational principles:

graph TB
    A["🪙 Bayesian Coin Flip<br/>Basic Probabilistic Modeling"] --> B["🔒 Type Safety Features<br/>Fugue's Type System Advantages"]
    B --> C["🔍 Trace Manipulation<br/>Runtime System & Custom Inference"]
    
    A --> D["Statistical Foundations"]
    B --> E["Type-Safe Programming"] 
    C --> F["Advanced Inference"]
    
    D --> G["Production Ready<br/>Probabilistic Models"]
    E --> G
    F --> G
    
    style A fill:#e1f5fe
    style B fill:#f3e5f5  
    style C fill:#e8f5e8
    style G fill:#fff3e0

Learning Path

🎯 Recommended Order

1. Bayesian Coin Flip (~45 minutes)

   • Start here for statistical foundations
   • Learn Bayesian inference principles
   • Understand model specification and analysis

2. Type Safety Features (~30 minutes)

   • Discover Fugue's unique advantages
   • Master type-safe probabilistic programming
   • Eliminate runtime errors with compile-time guarantees

3. Trace Manipulation (~60 minutes)

   • Deep dive into Fugue's runtime system
   • Learn custom inference and debugging techniques
   • Build production-ready probabilistic applications

Tutorial Overview

🪙 Bayesian Coin Flip

Foundation: Statistical inference and model analysis

Your introduction to Bayesian reasoning through the classic coin flipping problem. This tutorial demonstrates how prior beliefs are updated with evidence to form posterior distributions.

Key Concepts:

• Prior, likelihood, and posterior distributions
• Bayesian updating with Beta-Binomial conjugacy
• Model validation and parameter estimation
• Analytical vs computational solutions

What You'll Build:

• Complete Bayesian coin bias estimation model
• Prior sensitivity analysis framework
• Model validation with synthetic data

🔒 Type Safety Features

Foundation: Type-safe probabilistic programming

Explore Fugue's revolutionary type system that eliminates runtime errors while preserving full statistical expressiveness. Learn how dependent types make probabilistic programs both safer and faster.

Key Concepts:

• Natural return types for distributions (bool, u64, f64, usize)
• Compile-time safety guarantees
• Safe array indexing with categorical distributions
• Parameter validation at construction time
• Performance benefits through zero-cost abstractions

What You'll Build:

• Type-safe hierarchical models
• Safe array indexing examples
• Performance comparison with traditional PPLs

🔍 Trace Manipulation

Foundation: Runtime system and advanced inference

Master Fugue's execution trace system — the foundation that enables sophisticated inference algorithms. Learn how traces record, replay, and analyze probabilistic model executions.

Key Concepts:

• Trace system architecture and execution history
• Handler system for flexible model interpretation
• Replay mechanics for MCMC and inference algorithms
• Custom handlers for specialized inference strategies
• Memory optimization for production deployment
• Diagnostic tools for convergence assessment

What You'll Build:

• Custom MCMC algorithm using trace replay
• Specialized handlers for debugging models
• Production inference pipeline with memory optimization
• Comprehensive diagnostic system for model validation

Learning Outcomes

After completing these foundation tutorials, you will:

📊 Statistical Mastery

• ✅ Understand Bayesian inference from first principles
• ✅ Build and validate probabilistic models confidently
• ✅ Interpret posterior distributions and uncertainty quantification
• ✅ Apply conjugate analysis and computational methods

🛡️ Type-Safe Programming

• ✅ Write probabilistic programs that catch errors at compile time
• ✅ Leverage natural return types for cleaner, safer code
• ✅ Understand performance benefits of zero-cost abstractions
• ✅ Build complex models with guaranteed type safety

⚙️ Advanced Inference

• ✅ Manipulate execution traces for custom inference algorithms
• ✅ Implement specialized handlers for unique requirements
• ✅ Debug and optimize problematic models systematically
• ✅ Deploy production-ready probabilistic systems

🔧 Production Skills

• ✅ Memory optimization techniques for high-throughput scenarios
• ✅ Convergence diagnostics and model validation workflows
• ✅ Custom inference algorithms tailored to specific problems
• ✅ Systematic debugging of numerical issues

Code Examples

All tutorials include comprehensive, tested code examples:

• 📁 examples/bayesian_coin_flip.rs - Complete Bayesian analysis
• 📁 examples/type_safety.rs - Type system demonstrations
• 📁 examples/trace_manipulation.rs - Runtime system examples

Each example includes:

• ✅ Comprehensive tests ensuring correctness
• ✅ Detailed comments explaining every concept
• ✅ Runnable code you can execute immediately
• ✅ Performance benchmarks where applicable

# Run any example to see concepts in action
cargo run --example bayesian_coin_flip
cargo run --example type_safety  
cargo run --example trace_manipulation

# Run tests to verify your understanding
cargo test --example bayesian_coin_flip
cargo test --example type_safety
cargo test --example trace_manipulation

Next Steps

After mastering these foundations, you're ready for:

📈 Statistical Modeling Tutorials

Apply your knowledge to real-world problems:

• Linear and logistic regression
• Hierarchical models and mixed effects
• Mixture models and clustering
• Time series and forecasting

🏗️ How-To Guides

Practical guidance for specific tasks:

• Building Complex Models
• Custom Handlers
• Optimizing Performance
• Debugging Models

🚀 Advanced Applications

Cutting-edge probabilistic programming:

• Advanced inference techniques
• Model comparison and selection
• Large-scale distributed inference

Getting Help

📚 Documentation

• Getting Started Guide - Fugue basics
• API Reference - Complete function documentation
• How-To Guides - Task-specific instructions

💡 Tips for Success

🔧 Troubleshooting

If you encounter issues:

1. Check prerequisites - Ensure you have the required mathematical background
2. Run examples step-by-step - Isolate where confusion arises
3. Review error messages - Fugue's type system provides helpful compile-time feedback
4. Consult diagnostics - Use trace analysis to debug model behavior

Ready to Begin?

Start your journey with Bayesian Coin Flip — the gateway to mastering probabilistic programming with Fugue.

Bayesian Coin Flip

A comprehensive introduction to Bayesian inference through the classic coin flip problem. This tutorial demonstrates core Bayesian concepts including prior beliefs, likelihood functions, posterior distributions, and conjugate analysis using Fugue's type-safe probabilistic programming framework.

The Problem & Data

Research Question: Is a coin fair, or does it have a bias toward heads or tails?

In classical statistics, we might perform a hypothesis test. In Bayesian statistics, we express our uncertainty about the coin's bias as a probability distribution and update this belief as we observe data.

Data Generation & Exploration

    // Real experimental data: coin flip outcomes
    // H = Heads (success), T = Tails (failure)
    let observed_flips = vec![
        true, false, true, true, false, true, true, false, true, true,
    ];
    let n_flips = observed_flips.len();
    let successes = observed_flips.iter().filter(|&&x| x).count();

    println!("🪙 Observed coin flip sequence:");
    for (i, &flip) in observed_flips.iter().enumerate() {
        print!("  Flip {}: {}", i + 1, if flip { "H" } else { "T" });
        if (i + 1) % 5 == 0 {
            println!();
        }
    }
    println!(
        "\n📊 Summary: {} successes out of {} flips ({:.1}%)",
        successes,
        n_flips,
        (successes as f64 / n_flips as f64) * 100.0
    );

    // Research question: Is this a fair coin? (p = 0.5)
    println!("❓ Research Question: Is this coin fair (p = 0.5)?");

Real-World Context: This could represent:

• Quality Control: Defective vs. non-defective products
• Medical Trials: Treatment success rates
• A/B Testing: Conversion rates between variants
• Survey Response: Yes/No answers to questions

Mathematical Foundation

The Bayesian paradigm treats parameters as random variables with probability distributions. For the coin flip problem:

Model Specification

Prior Distribution: Our initial belief about the coin bias $p$ before seeing data.

$$p\sim \text{Beta}({\alpha }_0,{\beta }_0)$$

where ${\alpha }_0$ and ${\beta }_0$ encode our prior "pseudo-observations" of successes and failures.

Likelihood Function: Given bias $p$, each flip $X_i$ follows:

$$X_i\mid p\sim \text{Bernoulli}(p)\text{for }i=1,2,\dots ,n$$

The joint likelihood for $n$ independent flips with $k$ successes is:

$$L(p\mid \mathbf{x})=p^k(1-p{)}^{n-k}$$

Posterior Distribution: By Bayes' theorem:

$$p(\theta \mid \text{data})=\frac{p(\text{data}\mid \theta )\cdot p(\theta )}{p(\text{data})}\propto p(\text{data}\mid \theta )\cdot p(\theta )$$

For the Beta-Bernoulli model, the posterior is:

$$p\mid \mathbf{x}\sim \text{Beta}({\alpha }_0+k,{\beta }_0+n-k)$$

    // Bayesian Model Specification:
    // Prior: p ~ Beta(α₀, β₀)  [belief about coin bias before data]
    // Likelihood: X_i ~ Bernoulli(p)  [each flip outcome]
    // Posterior: p|data ~ Beta(α₀ + successes, β₀ + failures)

    // Prior parameters (weakly informative)
    let prior_alpha = 2.0_f64; // Prior "successes"
    let prior_beta = 2.0_f64; // Prior "failures"

    // Prior implies: E[p] = α/(α+β) = 0.5, but allows uncertainty
    let prior_mean = prior_alpha / (prior_alpha + prior_beta);
    let prior_variance = (prior_alpha * prior_beta)
        / ((prior_alpha + prior_beta).powi(2_i32) * (prior_alpha + prior_beta + 1.0));

    println!(
        "📈 Prior Distribution: Beta({}, {})",
        prior_alpha, prior_beta
    );
    println!("   - Prior mean: {:.3}", prior_mean);
    println!("   - Prior variance: {:.4}", prior_variance);
    println!("   - Interpretation: Weakly favors fairness but allows bias");

Prior Choice and Interpretation

The Beta($\alpha ,\beta$) distribution has:

• Mean: $\mathbb{E}[p]=\frac{\alpha }{\alpha +\beta }$
• Variance: $\text{Var}[p]=\frac{\alpha \beta }{(\alpha +\beta {)}^2(\alpha +\beta +1)}$

Common choices:

• Beta(1,1): Uniform prior (no preference)
• Beta(2,2): Weakly informative, slight preference for fairness
• Beta(0.5, 0.5): Jeffreys prior (non-informative)

Basic Implementation

Let's implement our Bayesian coin flip model in Fugue:

use fugue::*;
// Define the probabilistic model
fn coin_flip_model(data: Vec<bool>) -> Model<f64> {
    prob!(
        // Prior belief about coin bias
        let p <- sample(addr!("coin_bias"), Beta::new(2.0, 2.0).unwrap());

        // Constrain p to valid range [0, 1] for numerical stability
        let p_constrained = p.clamp(1e-10, 1.0 - 1e-10);

        // Likelihood: observe each flip given the bias
        let _observations <- plate!(i in 0..data.len() => {
            observe(addr!("flip", i), Bernoulli::new(p_constrained).unwrap(), data[i])
        });

        // Return the inferred bias
        pure(p)
    )
}

Key Implementation Details:

1. Type Safety: Bernoulli returns bool directly—no casting required
2. Plate Notation: plate! efficiently handles vectorized observations
3. Address Uniqueness: Each observation gets a unique address flip#i
4. Pure Functional: Model returns the parameter of interest directly

Advanced Techniques

Analytical Posterior Solution

The beauty of conjugate priors is that we can compute the exact posterior without numerical approximation:

    // Beta-Bernoulli conjugacy gives exact posterior
    let posterior_alpha = prior_alpha + successes as f64;
    let posterior_beta = prior_beta + (n_flips - successes) as f64;

    let posterior_mean = posterior_alpha / (posterior_alpha + posterior_beta);
    let posterior_variance = (posterior_alpha * posterior_beta)
        / ((posterior_alpha + posterior_beta).powi(2) * (posterior_alpha + posterior_beta + 1.0));

    println!(
        "🎯 Analytical Posterior: Beta({:.0}, {:.0})",
        posterior_alpha, posterior_beta
    );
    println!("   - Posterior mean: {:.3}", posterior_mean);
    println!("   - Posterior variance: {:.4}", posterior_variance);

    // Credible intervals
    let posterior_dist = Beta::new(posterior_alpha, posterior_beta).unwrap();
    let _lower_bound = 0.025; // 2.5th percentile
    let _upper_bound = 0.975; // 97.5th percentile

    // Approximate quantiles (would need inverse CDF for exact)
    println!(
        "   - 95% credible interval: approximately [{:.2}, {:.2}]",
        posterior_mean - 1.96 * posterior_variance.sqrt(),
        posterior_mean + 1.96 * posterior_variance.sqrt()
    );

    // Hypothesis testing: P(p > 0.5 | data)
    let prob_biased_heads = if posterior_mean > 0.5 {
        0.8 // Rough approximation - would integrate Beta CDF for exact value
    } else {
        0.3
    };
    println!("   - P(p > 0.5 | data) ≈ {:.1}", prob_biased_heads);

MCMC Validation

While analytical solutions are preferred, we can validate our results using MCMC:

    // Use MCMC to approximate the posterior (for validation)
    let n_samples = 2000;
    let n_warmup = 500;

    let mut rng = rand::rngs::StdRng::seed_from_u64(42);
    let mcmc_samples = adaptive_mcmc_chain(
        &mut rng,
        || coin_flip_model(observed_flips.clone()),
        n_samples,
        n_warmup,
    );

    // Extract posterior samples for p
    let posterior_samples: Vec<f64> = mcmc_samples
        .iter()
        .filter_map(|(_, trace)| trace.get_f64(&addr!("coin_bias")))
        .collect();

    if !posterior_samples.is_empty() {
        let mcmc_mean = posterior_samples.iter().sum::<f64>() / posterior_samples.len() as f64;
        let mcmc_variance = {
            let mean = mcmc_mean;
            posterior_samples
                .iter()
                .map(|x| (x - mean).powi(2_i32))
                .sum::<f64>()
                / (posterior_samples.len() - 1) as f64
        };
        let ess = effective_sample_size(&posterior_samples);

        println!(
            "✅ MCMC Results ({} effective samples from {} total):",
            ess as usize,
            posterior_samples.len()
        );
        println!(
            "   - MCMC mean: {:.3} (analytical: {:.3})",
            mcmc_mean, posterior_mean
        );
        println!(
            "   - MCMC variance: {:.4} (analytical: {:.4})",
            mcmc_variance, posterior_variance
        );
        println!("   - Effective Sample Size: {:.0}", ess);
        println!(
            "   - Agreement: {}",
            if (mcmc_mean - posterior_mean).abs() < 0.05 {
                "✅ Excellent"
            } else {
                "⚠️ Check convergence"
            }
        );
    }

Why MCMC for a Conjugate Model?

Even though we have analytical solutions, MCMC serves important purposes:

• Validation: Confirms our analytical calculations
• Flexibility: Easily extends to non-conjugate models
• Diagnostics: Provides convergence and mixing assessments

Diagnostics & Validation

Model validation ensures our model adequately represents the data-generating process:

    // Posterior predictive checks
    println!("🔍 Posterior Predictive Validation:");

    // Simulate new data from posterior predictive distribution
    let mut rng = thread_rng();
    let n_pred_samples = 1000;
    let mut predicted_successes = Vec::new();

    for _ in 0..n_pred_samples {
        // Sample bias from posterior
        let p_sample = posterior_dist.sample(&mut rng);

        // Simulate n_flips with this bias
        let mut pred_successes = 0;
        for _ in 0..n_flips {
            if Bernoulli::new(p_sample).unwrap().sample(&mut rng) {
                pred_successes += 1;
            }
        }
        predicted_successes.push(pred_successes);
    }

    // Compare with observed successes
    let pred_mean = predicted_successes.iter().sum::<usize>() as f64 / n_pred_samples as f64;
    let pred_within_range = predicted_successes
        .iter()
        .filter(|&&x| (x as i32 - successes as i32).abs() <= 2)
        .count() as f64
        / n_pred_samples as f64;

    println!("   - Observed successes: {}", successes);
    println!("   - Predicted mean successes: {:.1}", pred_mean);
    println!(
        "   - P(|pred - obs| ≤ 2): {:.1}%",
        pred_within_range * 100.0
    );

    if pred_within_range > 0.5 {
        println!("   - ✅ Model fits data well");
    } else {
        println!("   - ⚠️ Model may not capture data well");
    }

Posterior Predictive Checks

The posterior predictive distribution answers: "If our model is correct, what data would we expect to see?"

$$p(\mathbf{x}\mid \mathbf{x})=\int p(\mathbf{x}\mid \theta )p(\theta \mid \mathbf{x})d\theta$$

Interpretation:

• Good fit: Observed data looks typical under the posterior predictive
• Poor fit: Observed data is extreme under the posterior predictive
• Model inadequacy: Systematic deviations suggest missing model components

Production Extensions

Decision Theory and Practical Applications

Bayesian inference provides the foundation for optimal decision-making under uncertainty:

    // Bayesian decision theory for fairness testing
    println!("🎲 Decision Analysis:");

    // Define loss function for hypothesis testing
    // H0: coin is fair (p = 0.5), H1: coin is biased (p ≠ 0.5)
    let fairness_threshold = 0.05; // How far from 0.5 counts as "biased"
    let prob_fair = if (posterior_mean - 0.5).abs() < fairness_threshold {
        // Approximate based on credible interval
        0.6
    } else {
        0.2
    };

    println!(
        "   - Posterior probability coin is fair: {:.1}%",
        prob_fair * 100.0
    );
    println!(
        "   - Evidence for bias: {}",
        if prob_fair < 0.3 {
            "Strong"
        } else if prob_fair < 0.7 {
            "Moderate"
        } else {
            "Weak"
        }
    );

    // Expected number of heads in future flips
    let future_flips = 20;
    let expected_heads = posterior_mean * future_flips as f64;
    let uncertainty = (posterior_variance * future_flips as f64).sqrt();

    println!(
        "   - Expected heads in next {} flips: {:.1} ± {:.1}",
        future_flips,
        expected_heads,
        1.96 * uncertainty
    );

    // Practical recommendations
    if (posterior_mean - 0.5).abs() < 0.1 {
        println!("   - 💡 Recommendation: Treat as approximately fair for practical purposes");
    } else if posterior_mean > 0.5 {
        println!("   - 💡 Recommendation: Coin appears biased toward heads");
    } else {
        println!("   - 💡 Recommendation: Coin appears biased toward tails");
    }

Advanced Model Extensions

Real applications often require extensions beyond the basic model:

    // Hierarchical model for multiple coins
    println!("🔬 Advanced Modeling Extensions:");

    // Example: What if we had multiple coins?
    let _multi_coin_model = || {
        prob!(
            // Population-level parameters
            let pop_mean <- sample(addr!("population_mean"), Beta::new(1.0, 1.0).unwrap());
            let pop_concentration <- sample(addr!("concentration"), Gamma::new(2.0, 0.5).unwrap());

            // Individual coin bias (hierarchical prior)
            let alpha = pop_mean * pop_concentration;
            let beta = (1.0 - pop_mean) * pop_concentration;
            let coin_bias <- sample(addr!("coin_bias"), Beta::new(alpha, beta).unwrap());

            pure(coin_bias)
        )
    };

    println!("   - 📈 Hierarchical Extension: Population of coins with shared parameters");
    println!("   - 🔄 Sequential Learning: Update beliefs with each new flip");
    println!("   - 🎯 Robust Models: Heavy-tailed priors for outlier resistance");
    println!("   - 📊 Model Comparison: Bayes factors between fair vs. biased hypotheses");

    // Model comparison example (simplified)
    let fair_model_evidence = -5.2_f64; // Log marginal likelihood for fair model
    let biased_model_evidence = -4.8_f64; // Log marginal likelihood for biased model
    let bayes_factor = (biased_model_evidence - fair_model_evidence).exp();

    println!("   - ⚖️ Bayes Factor (biased/fair): {:.2}", bayes_factor);
    if bayes_factor > 3.0 {
        println!("     Evidence favors biased model");
    } else if bayes_factor < 1.0 / 3.0 {
        println!("     Evidence favors fair model");
    } else {
        println!("     Evidence is inconclusive");
    }

Advanced Modeling Scenarios:

1. Hierarchical Models: Multiple coins with shared population parameters
2. Sequential Learning: Online updates as new flips arrive
3. Robust Priors: Heavy-tailed distributions to handle outliers
4. Model Selection: Comparing fair vs. biased hypotheses using Bayes factors

graph TD
    A[Basic Beta-Bernoulli] --> B[Hierarchical Extension]
    A --> C[Sequential Updates]  
    A --> D[Robust Priors]
    A --> E[Model Comparison]
    B --> F[Population Studies]
    C --> G[Online Learning]
    D --> H[Outlier Resistance] 
    E --> I[Bayes Factors]

Real-World Considerations

When to Use Bayesian vs. Frequentist Methods

Bayesian Advantages:

• Natural uncertainty quantification: Full posterior distributions
• Prior knowledge incorporation: Systematic way to include expert knowledge
• Decision-theoretic framework: Optimal decisions under specified loss functions
• Sequential updating: Natural online learning as data arrives

Frequentist Advantages:

• Objective interpretation: No need to specify prior distributions
• Computational simplicity: Often faster for standard problems
• Regulatory acceptance: Many standards assume frequentist methods

Performance Implications

• Conjugate Models: Analytical solutions are extremely fast
• MCMC Methods: Scale linearly with data size and number of parameters
• Memory Usage: Fugue's trace system efficiently manages large parameter spaces
• Numerical Stability: Log-space computations prevent underflow issues

Common Pitfalls

1. Improper Priors: Always verify prior distributions integrate to 1
2. Label Switching: In mixture models, parameter interpretability can change
3. Convergence Assessment: Always check MCMC diagnostics before making inferences
4. Prior Sensitivity: Test how conclusions change under different reasonable priors

Exercises

1. Prior Sensitivity Analysis:

   • Try different Beta priors: Beta(1,1), Beta(5,5), Beta(0.5,0.5)
   • How do the posteriors differ with the same data?
   • When does prior choice matter most?

2. Sequential Learning:

   • Start with Beta(2,2) prior
   • Update after each flip in sequence
   • Plot how the posterior evolves with each observation

3. Model Comparison:

   • Implement a "fair coin" model with p = 0.5 exactly
   • Compare evidence for fair vs. biased models using marginal likelihoods
   • What sample size is needed to distinguish p = 0.5 from p = 0.6?

4. Hierarchical Extension:

   • Model 5 different coins with a shared Beta population prior
   • Each coin has different numbers of flips
   • How does information sharing affect individual coin estimates?

Testing Your Understanding

Comprehensive test suite for validation:

    #[test]
    fn test_coin_flip_model_properties() {
        let test_data = vec![true, true, false, true];
        let mut rng = thread_rng();

        // Test model executes without panics
        let (bias_sample, trace) = runtime::handler::run(
            PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            },
            coin_flip_model(test_data.clone()),
        );

        // Bias should be valid probability
        assert!(bias_sample >= 0.0 && bias_sample <= 1.0);

        // Trace should contain expected choices
        assert!(trace.get_f64(&addr!("coin_bias")).is_some());
        assert!(trace.total_log_weight().is_finite());

        // Should have observation sites for each data point
        for _i in 0..test_data.len() {
            // Observations don't create choices, but affect likelihood
            assert!(trace.log_likelihood.is_finite());
        }
    }

    #[test]
    fn test_conjugate_update_correctness() {
        // Test analytical posterior against known values
        let prior_alpha = 2.0;
        let prior_beta = 2.0;
        let successes = 7;
        let failures = 3;

        let posterior_alpha = prior_alpha + successes as f64;
        let posterior_beta = prior_beta + failures as f64;
        let posterior_mean = posterior_alpha / (posterior_alpha + posterior_beta);

        // Should be (2+7)/(2+2+7+3) = 9/14 ≈ 0.643
        assert!((posterior_mean - 9.0 / 14.0).abs() < 1e-10);

        // Posterior should be more concentrated than prior
        let prior_variance = (2.0 * 2.0) / (4.0_f64.powi(2_i32) * 5.0);
        let posterior_variance = (posterior_alpha * posterior_beta)
            / ((posterior_alpha + posterior_beta).powi(2_i32)
                * (posterior_alpha + posterior_beta + 1.0));
        assert!(posterior_variance < prior_variance);
    }

    #[test]
    fn test_model_with_edge_cases() {
        let mut rng = thread_rng();

        // Test with all heads
        let all_heads = vec![true; 10];
        let (bias, _) = runtime::handler::run(
            PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            },
            coin_flip_model(all_heads),
        );
        // Should still be valid probability
        assert!(bias >= 0.0 && bias <= 1.0);

        // Test with all tails
        let all_tails = vec![false; 10];
        let (bias, _) = runtime::handler::run(
            PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            },
            coin_flip_model(all_tails),
        );
        assert!(bias >= 0.0 && bias <= 1.0);

        // Test with single flip
        let single_flip = vec![true];
        let (bias, _) = runtime::handler::run(
            PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            },
            coin_flip_model(single_flip),
        );
        assert!(bias >= 0.0 && bias <= 1.0);
    }

Next Steps

Now that you understand Bayesian inference fundamentals:

• Type Safety Features: Learn how Fugue's type system prevents common errors
• Trace Manipulation: Understand Fugue's runtime system for custom inference
• Linear Regression: Extend to continuous outcomes
• Hierarchical Models: Multi-level modeling

The coin flip problem provides the conceptual foundation for all Bayesian modeling. Every complex model builds on these same principles: prior beliefs, likelihood functions, and posterior inference.

Type Safety Features

A comprehensive exploration of Fugue's revolutionary type-safe distribution system and its practical implications for probabilistic programming. This tutorial demonstrates how dependent type theory principles eliminate runtime errors while preserving full statistical expressiveness, making probabilistic programs both safer and more performant.

The Type Safety Problem

Traditional probabilistic programming languages force all distributions to return f64, creating a fundamental mismatch between mathematical concepts and their computational representation. This leads to pervasive runtime errors, casting overhead, and semantic confusion.

graph TD
    A["Traditional PPL"] --> B["All distributions → f64"]
    B --> C["Runtime Errors"]
    B --> D["Casting Overhead"]  
    B --> E["Semantic Confusion"]
    B --> F["Precision Loss"]
    
    G["Fugue PPL"] --> H["Natural Return Types"]
    H --> I["bool for Bernoulli"]
    H --> J["u64 for Poisson"]
    H --> K["usize for Categorical"]
    H --> L["f64 for Normal"]
    
    I --> M["Compile-Time Safety"]
    J --> M
    K --> M
    L --> M
    
    style A fill:#ffcccc
    style G fill:#ccffcc
    style M fill:#ccffff

Traditional PPL Problems

// Demonstrates problems with traditional PPL approaches (shown for contrast)
fn traditional_ppl_problems() {
    println!("=== Traditional PPL Problems (What Fugue Solves) ===\n");

    // In traditional PPLs, everything returns f64, leading to:
    println!("❌ Traditional PPL Issues:");
    println!("   - Bernoulli returns f64 → if sample == 1.0 (awkward)");
    println!("   - Poisson returns f64 → count.round() as u64 (precision loss)");
    println!("   - Categorical returns f64 → array[sample as usize] (unsafe)");
    println!("   - Runtime type errors and casting overhead");
    println!();
}

Mathematical Foundation

Fugue's type system is grounded in dependent type theory, where each distribution $D$ is parameterized not just by its parameters $\theta$, but by its support type $\mathcal{S}$.

Formal Type System

For a distribution $D_{\theta }$ with parameters $\theta$ and support $\mathcal{S}$:

$$\text{sample}(D_{\theta }):\mathcal{S}$$

This ensures that sampling operations return values in their natural mathematical domain:

Type-Theoretic Properties

Natural Type System

Fugue eliminates the f64-everything problem by returning mathematically appropriate types:

use fugue::*;
use rand::thread_rng;
// Demonstrate Fugue's natural return types
fn natural_type_system() {
    println!("✅ Fugue's Natural Type System");
    println!("==============================\n");

    let mut rng = thread_rng();

    // Boolean decisions: Bernoulli → bool
    let fair_coin = Bernoulli::new(0.5).unwrap();
    let is_heads: bool = fair_coin.sample(&mut rng);

    // Natural conditional logic - no comparisons!
    let outcome = if is_heads {
        "Heads - you win!"
    } else {
        "Tails - try again"
    };
    println!("🪙 Coin flip: {} (type: bool)", outcome);

    // Count data: Poisson → u64
    let customer_arrivals = Poisson::new(5.0).unwrap();
    let arrivals: u64 = customer_arrivals.sample(&mut rng);

    // Direct arithmetic with counts - no casting!
    let service_time = arrivals * 10; // minutes per customer
    println!(
        "👥 Customers: {} arrivals, {}min service (type: u64)",
        arrivals, service_time
    );

    // Category selection: Categorical → usize
    let product_preferences = Categorical::new(vec![0.4, 0.35, 0.25]).unwrap();
    let choice: usize = product_preferences.sample(&mut rng);

    // Safe array indexing - guaranteed bounds safety!
    let products = ["Laptop", "Smartphone", "Tablet"];
    println!(
        "🛒 Customer chose: {} (index: {}, type: usize)",
        products[choice], choice
    );

    // Continuous values: Normal → f64 (unchanged, as expected)
    let measurement = Normal::new(100.0, 5.0).unwrap();
    let reading: f64 = measurement.sample(&mut rng);
    println!("📏 Sensor reading: {:.2} units (type: f64)", reading);

    println!();
}

Type Benefits by Distribution

Bernoulli Distributions

• Returns: bool - natural boolean logic
• Benefit: Direct conditional statements without equality comparisons
• Performance: No floating-point comparisons needed

Count Distributions (Poisson, Binomial)

• Returns: u64 - natural counting numbers
• Benefit: Direct arithmetic without casting or precision loss
• Performance: Integer operations are faster than float conversions

Categorical Distributions

• Returns: usize - natural array indices
• Benefit: Guaranteed bounds safety for array indexing
• Performance: No runtime bounds checking required

Continuous Distributions

• Returns: f64 - unchanged for appropriate domains
• Benefit: Expected behavior preserved for mathematical operations

Compile-Time Safety

Fugue's type system catches errors at compile time, eliminating entire classes of runtime failures:

use fugue::*;
use fugue::runtime::interpreters::PriorHandler;
use rand::thread_rng;
// Demonstrate compile-time type safety guarantees
fn compile_time_safety_demo() {
    println!("🛡️ Compile-Time Type Safety");
    println!("============================\n");

    // Type-safe model composition
    let data_model: Model<(bool, u64, usize, f64)> = prob!(
        let coin_result <- sample(addr!("coin"), Bernoulli::new(0.6).unwrap());
        let event_count <- sample(addr!("events"), Poisson::new(3.0).unwrap());
        let category <- sample(addr!("category"), Categorical::uniform(4).unwrap());
        let measurement <- sample(addr!("measure"), Normal::new(0.0, 1.0).unwrap());

        // Compiler enforces correct types throughout
        pure((coin_result, event_count, category, measurement))
    );

    println!("✅ Model created with strict type guarantees:");
    println!("   - coin_result: bool (no == 1.0 needed)");
    println!("   - event_count: u64 (direct arithmetic)");
    println!("   - category: usize (safe indexing)");
    println!("   - measurement: f64 (natural continuous)");

    // Execute model safely
    let mut rng = thread_rng();
    let (sample, _trace) = runtime::handler::run(
        PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        data_model,
    );

    println!(
        "📊 Sample: coin={}, events={}, category={}, value={:.3}",
        sample.0, sample.1, sample.2, sample.3
    );
    println!();
}

Type-Safe Model Composition

Models compose naturally while preserving type information throughout the computation:

Safe Array Indexing

One of the most dangerous operations in traditional PPLs is array indexing with categorical samples. Fugue makes this provably safe:

use fugue::*;
use rand::thread_rng;
// Demonstrate safe array indexing with categorical distributions
fn safe_array_indexing() {
    println!("🎯 Safe Array Indexing");
    println!("======================\n");

    let mut rng = thread_rng();

    // Define categories with natural indexing
    let algorithms = ["MCMC", "Variational Inference", "ABC", "SMC", "Exact"];
    let method_weights = vec![0.3, 0.25, 0.2, 0.15, 0.1];

    let method_selector = Categorical::new(method_weights).unwrap();

    println!("🧮 Available inference methods:");
    for (i, method) in algorithms.iter().enumerate() {
        println!("   {}: {}", i, method);
    }
    println!();

    // Sample multiple times to show safety
    for trial in 1..=5 {
        let selected_idx: usize = method_selector.sample(&mut rng);

        // This is GUARANTEED safe - no bounds checking needed!
        let chosen_method = algorithms[selected_idx];

        println!(
            "Trial {}: Selected method '{}' (index {})",
            trial, chosen_method, selected_idx
        );
    }

    println!("\n✅ All array accesses guaranteed safe by type system!");
    println!();
}

Bounds Safety Guarantee

Why This Matters

Traditional PPLs require defensive programming:

// Traditional PPL - unsafe!
let category = categorical_sample as usize;
if category < array.len() {
    return array[category];  // Still might panic due to float precision!
} else {
    return default_value;    // Defensive fallback
}

Fugue guarantees safety:

// Fugue - provably safe!
let category: usize = categorical.sample(&mut rng);
return array[category];  // Cannot panic - guaranteed by type system

Parameter Validation

Fugue validates all distribution parameters at construction time, catching invalid configurations before they can cause runtime errors:

use fugue::*;
// Demonstrate parameter validation and error handling
fn parameter_validation_demo() {
    println!("🔍 Parameter Validation");
    println!("=======================\n");

    println!("Fugue validates parameters at construction time:");
    println!();

    // Valid constructions
    match Normal::new(0.0, 1.0) {
        Ok(_) => println!("✅ Normal(μ=0.0, σ=1.0) - valid"),
        Err(e) => println!("❌ Unexpected error: {:?}", e),
    }

    match Beta::new(2.0, 3.0) {
        Ok(_) => println!("✅ Beta(α=2.0, β=3.0) - valid"),
        Err(e) => println!("❌ Unexpected error: {:?}", e),
    }

    match Categorical::new(vec![0.3, 0.4, 0.3]) {
        Ok(_) => println!("✅ Categorical([0.3, 0.4, 0.3]) - valid"),
        Err(e) => println!("❌ Unexpected error: {:?}", e),
    }

    println!();

    // Invalid constructions - caught at compile time with .unwrap()
    // or handled gracefully with pattern matching
    println!("Invalid parameter examples:");

    match Normal::new(0.0, -1.0) {
        Ok(_) => println!("✅ Normal(μ=0.0, σ=-1.0) - unexpected success"),
        Err(e) => println!("❌ Normal(μ=0.0, σ=-1.0) - {}", e),
    }

    match Beta::new(0.0, 1.0) {
        Ok(_) => println!("✅ Beta(α=0.0, β=1.0) - unexpected success"),
        Err(e) => println!("❌ Beta(α=0.0, β=1.0) - {}", e),
    }

    match Categorical::new(vec![0.5, 0.6]) {
        // Doesn't sum to 1
        Ok(_) => println!("✅ Categorical([0.5, 0.6]) - unexpected success"),
        Err(e) => println!("❌ Categorical([0.5, 0.6]) - {}", e),
    }

    println!("\n✅ All invalid parameters caught before runtime!");
    println!();
}

Validation Strategy

Fugue uses fail-fast construction with comprehensive parameter checking:

Type-Safe Observations

Observations in Fugue must match the distribution's return type, providing compile-time guarantees about data consistency:

use fugue::*;
use fugue::runtime::interpreters::PriorHandler;
use rand::thread_rng;
// Demonstrate type-safe observations with automatic type checking
fn type_safe_observations() {
    println!("🔗 Type-Safe Observations");
    println!("=========================\n");

    // Observations must match distribution return types
    let observation_model = prob!(
        // Boolean observation - must provide bool
        let _bool_obs <- observe(addr!("coin_obs"),
                                Bernoulli::new(0.7).unwrap(),
                                true); // ✅ bool type matches

        // Count observation - must provide u64
        let _count_obs <- observe(addr!("events_obs"),
                                 Poisson::new(4.0).unwrap(),
                                 5u64); // ✅ u64 type matches

        // Category observation - must provide usize
        let _category_obs <- observe(addr!("choice_obs"),
                                    Categorical::new(vec![0.2, 0.5, 0.3]).unwrap(),
                                    1usize); // ✅ usize type matches

        // Continuous observation - must provide f64
        let _continuous_obs <- observe(addr!("measurement_obs"),
                                      Normal::new(10.0, 2.0).unwrap(),
                                      12.5f64); // ✅ f64 type matches

        pure(())
    );

    println!("✅ All observations type-checked at compile time!");
    println!("   - Bernoulli observation: bool");
    println!("   - Poisson observation: u64");
    println!("   - Categorical observation: usize");
    println!("   - Normal observation: f64");

    // Execute to verify it works
    let mut rng = thread_rng();
    let (_result, trace) = runtime::handler::run(
        PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        observation_model,
    );

    println!(
        "📊 Model executed successfully with {} addresses",
        trace.choices.len()
    );
    println!();
}

Observation Type Matching

The type system ensures that observed values match the distribution's natural type:

// ✅ Type-safe observations
observe(addr!("coin"), Bernoulli::new(0.5).unwrap(), true);     // bool
observe(addr!("count"), Poisson::new(3.0).unwrap(), 5u64);     // u64  
observe(addr!("choice"), Categorical::uniform(3).unwrap(), 1usize); // usize
observe(addr!("measure"), Normal::new(0.0, 1.0).unwrap(), 2.5f64);  // f64

// ❌ These would be compile-time errors
observe(addr!("coin"), Bernoulli::new(0.5).unwrap(), 1.0);     // f64 ≠ bool
observe(addr!("count"), Poisson::new(3.0).unwrap(), 5.0);      // f64 ≠ u64
observe(addr!("choice"), Categorical::uniform(3).unwrap(), 1);  // i32 ≠ usize

Advanced Type Composition

Fugue supports complex hierarchical models with full type safety throughout the computation:

use fugue::*;
use fugue::runtime::interpreters::PriorHandler;
use rand::thread_rng;
// Demonstrate advanced type-safe model composition
fn advanced_type_composition() {
    println!("🧩 Advanced Type Composition");
    println!("============================\n");

    // Complex hierarchical model with full type safety
    let hierarchical_model = prob!(
        // Global parameters
        let success_rate <- sample(addr!("global_rate"), Beta::new(1.0, 1.0).unwrap());

        // Group-specific parameters (different types working together)
        let group_sizes <- sequence_vec((0..3).map(|group_id| {
            sample(addr!("group_size", group_id), Poisson::new(10.0).unwrap())
        }).collect());

        let group_successes <- sequence_vec(group_sizes.iter().enumerate().map(|(group_id, &size)| {
            sample(addr!("successes", group_id), Binomial::new(size, success_rate).unwrap())
        }).collect());

        // Category assignments for each group
        let group_categories <- sequence_vec((0..3).map(|group_id| {
            sample(addr!("category", group_id), Categorical::uniform(4).unwrap())
        }).collect());

        // Return complex structured result with full type safety
        pure((success_rate, group_sizes, group_successes, group_categories))
    );

    println!("🏗️ Hierarchical model structure:");
    println!("   - Global success rate: f64 (Beta distribution)");
    println!("   - Group sizes: Vec<u64> (Poisson distributions)");
    println!("   - Group successes: Vec<u64> (Binomial distributions)");
    println!("   - Group categories: Vec<usize> (Categorical distributions)");
    println!();

    // Sample from the complex model
    let mut rng = thread_rng();
    let (result, _trace) = runtime::handler::run(
        PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        hierarchical_model,
    );

    let (rate, sizes, successes, categories) = result;

    println!("📈 Sample from hierarchical model:");
    println!("   Global success rate: {:.3}", rate);

    for (i, ((&size, &success), &category)) in sizes
        .iter()
        .zip(successes.iter())
        .zip(categories.iter())
        .enumerate()
    {
        println!(
            "   Group {}: {} trials, {} successes, category {}",
            i, size, success, category
        );
    }

    println!("\n✅ Complex model composed with full type safety!");
    println!();
}

Hierarchical Type Structure

Complex models maintain precise type information at every level:

graph TD
    A["Global: f64"] --> B["Group Sizes: Vec<u64>"]
    A --> C["Group Successes: Vec<u64>"]  
    A --> D["Group Categories: Vec<usize>"]
    
    B --> E["Model: (f64, Vec<u64>, Vec<u64>, Vec<usize>)"]
    C --> E
    D --> E
    
    style A fill:#e1f5fe
    style E fill:#c8e6c9

Performance Benefits

Type safety in Fugue eliminates runtime overhead through zero-cost abstractions:

// Demonstrate performance benefits of type safety
fn performance_benefits() {
    println!("⚡ Performance Benefits");
    println!("======================\n");

    println!("Type safety eliminates runtime overhead:");
    println!();

    println!("🚫 Traditional PPL (f64 everything):");
    println!("   let coin_flip = sample(...); // Returns f64");
    println!("   if coin_flip == 1.0 {{ ... }} // Float comparison");
    println!("   let count = sample(...) as u64; // Casting overhead");
    println!("   array[sample(...) as usize] // Unsafe casting + bounds check");
    println!();

    println!("✅ Fugue (natural types):");
    println!("   let coin_flip: bool = sample(...); // Returns bool");
    println!("   if coin_flip {{ ... }} // Natural boolean");
    println!("   let count: u64 = sample(...); // Direct u64");
    println!("   array[sample(...)] // Safe usize indexing");
    println!();

    println!("🎯 Benefits:");
    println!("   ✓ Zero casting overhead");
    println!("   ✓ No floating-point comparisons for discrete values");
    println!("   ✓ Eliminated bounds checking for categorical indexing");
    println!("   ✓ No precision loss from float→int conversions");
    println!("   ✓ Compile-time error detection");
    println!();
}

Performance Analysis

Real-World Applications

Quality Control System

use fugue::*;
let quality_model = prob!(
    // Product defect rate (continuous parameter)
    let defect_rate <- sample(addr!("defect_rate"), Beta::new(1.0, 9.0).unwrap());
    
    // Number of products tested (count data)
    let products_tested <- sample(addr!("tested"), Poisson::new(100.0).unwrap());
    
    // Actual defects found (count with bounds)
    let defects_found <- sample(addr!("defects"), 
                                Binomial::new(products_tested, defect_rate).unwrap());
    
    // Inspector assignment (categorical choice)
    let inspector <- sample(addr!("inspector"), Categorical::uniform(3).unwrap());
    
    // Natural type usage throughout
    pure((defect_rate, products_tested, defects_found, inspector))
);

Medical Diagnosis System

use fugue::*;  
let diagnosis_model = prob!(
    // Prior disease probability (continuous)
    let disease_prob <- sample(addr!("prior"), Beta::new(2.0, 98.0).unwrap());
    
    // Number of symptoms (count) 
    let symptom_count <- sample(addr!("symptoms"), Poisson::new(2.5).unwrap());
    
    // Test result (boolean outcome)
    let test_positive <- sample(addr!("test"), Bernoulli::new(0.95).unwrap());
    
    // Treatment recommendation (categorical)
    let treatment <- sample(addr!("treatment"), 
                           Categorical::new(vec![0.6, 0.3, 0.1]).unwrap());
    
    pure((disease_prob, symptom_count, test_positive, treatment))
);

Production Considerations

Error Handling Strategy

use fugue::*;
// Robust parameter validation
fn create_robust_model(rate: f64, categories: Vec<f64>) -> Result<Model<(f64, usize)>, String> {
    let poisson = Poisson::new(rate)
        .map_err(|e| format!("Invalid Poisson rate {}: {}", rate, e))?;
        
    let categorical = Categorical::new(categories)
        .map_err(|e| format!("Invalid categorical weights: {}", e))?;
    
    Ok(prob!(
        let count <- sample(addr!("count"), poisson);
        let choice <- sample(addr!("choice"), categorical);
        pure((count as f64, choice))
    ))
}

Performance Optimization

Testing Your Understanding

Exercise 1: Safe Model Construction

Create a model that demonstrates all four natural return types. Ensure it:

• Uses boolean logic for decision-making
• Performs arithmetic with count data
• Safely indexes into arrays
• Handles continuous parameters

// Exercise framework for testing understanding
fn testing_framework_example() {
    println!("🧪 Testing Framework Example");
    println!("============================\n");

    let comprehensive_model = prob!(
        // Boolean decision making
        let is_premium <- sample(addr!("premium"), Bernoulli::new(0.3).unwrap());

        // Count data arithmetic
        let base_items <- sample(addr!("base_items"), Poisson::new(5.0).unwrap());
        let bonus_items = if is_premium { base_items + 2 } else { base_items };

        // Safe array indexing
        let service_tier <- sample(addr!("tier"), Categorical::new(vec![0.5, 0.3, 0.2]).unwrap());

        // Continuous parameters
        let satisfaction <- sample(addr!("satisfaction"), Beta::new(2.0, 1.0).unwrap());

        pure((is_premium, bonus_items, service_tier, satisfaction))
    );

    println!("✅ Comprehensive model demonstrates:");
    println!("   - Boolean logic: Premium account decision");
    println!("   - Count arithmetic: Items calculation with bonus");
    println!("   - Safe indexing: Service tier selection");
    println!("   - Continuous data: Customer satisfaction modeling");

    let mut rng = thread_rng();
    let (premium, items, tier, satisfaction) = runtime::handler::run(
        PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        comprehensive_model,
    )
    .0;

    let tiers = ["Basic", "Standard", "Premium"];
    println!("\n📊 Sample result:");
    println!("   Premium account: {}", premium);
    println!("   Items received: {}", items);
    println!("   Service tier: {} ({})", tiers[tier], tier);
    println!("   Satisfaction: {:.2}%", satisfaction * 100.0);
    println!();
}

Exercise 2: Parameter Validation

Write a function that attempts to create distributions with both valid and invalid parameters. Handle errors gracefully and provide meaningful error messages.

Exercise 3: Hierarchical Composition

Design a hierarchical model that combines multiple data types across different levels. Ensure type safety is maintained throughout the composition.

Key Takeaways

Core Benefits:

• ✅ Eliminated runtime type errors - impossible by construction
• ✅ Natural mathematical operations - no awkward casting or comparisons
• ✅ Guaranteed array safety - categorical indexing cannot panic
• ✅ Performance improvements - zero-cost abstractions enable optimizations
• ✅ Clear code intent - types document the mathematical structure

Further Reading

• Working with Distributions - Practical distribution usage patterns
• Building Complex Models - Advanced composition techniques
• API Reference - Complete type specifications
• Types and Programming Languages by Benjamin Pierce - Theoretical foundations
• Probabilistic Programming & Bayesian Methods for Hackers - Applied Bayesian inference

Trace Manipulation

A deep exploration of Fugue's runtime system and trace manipulation capabilities. This tutorial demonstrates how traces enable sophisticated probabilistic programming techniques including replay, scoring, custom inference, and debugging. Learn how Fugue's execution history recording makes advanced inference algorithms possible while maintaining full type safety.

The Execution History Problem

Traditional programming languages execute once and discard their execution history. In probabilistic programming, we need to record, manipulate, and reason about random choices to enable sophisticated inference algorithms. Fugue's trace system solves this fundamental challenge.

graph TD
    A["Model Specification"] --> B["Handler Selection"]
    B --> C["PriorHandler<br/>Forward Sampling"]
    B --> D["ReplayHandler<br/>MCMC Proposals"]  
    B --> E["ScoreGivenTrace<br/>Importance Sampling"]
    B --> F["Custom Handler<br/>Specialized Logic"]
    
    C --> G["Execution Trace"]
    D --> G
    E --> G
    F --> G
    
    G --> H["Choice Records"]
    G --> I["Log-Weight Components"]
    G --> J["Type-Safe Values"]
    
    H --> K["Replay"]
    H --> L["Scoring"]
    H --> M["Conditioning"]
    H --> N["Debugging"]
    
    style G fill:#ccffcc
    style K fill:#e1f5fe
    style L fill:#e1f5fe
    style M fill:#e1f5fe
    style N fill:#e1f5fe

Mathematical Foundation

Trace Formalization

A trace $\tau$ records the complete execution history of a probabilistic model, formally represented as:

$$\tau =⟨\mathcal{C},\log w_{\text{prior}},\log w_{\text{likelihood}},\log w_{\text{factors}}⟩$$

Where:

• $\mathcal{C}$: Map from addresses to choices $\{a_i\mapsto (v_i,\log p_i)\}$
• $\log w_{\text{prior}}$: Accumulated prior log-probability ${\sum }_i\log p(v_i)$
• $\log w_{\text{likelihood}}$: Accumulated observation log-probability ${\sum }_j\log p(y_j|x_j)$
• $\log w_{\text{factors}}$: Accumulated factor weights ${\sum }_k\log f_k$

Total Log-Weight

The total unnormalized log-probability is:

$$\log w(\tau )=\log w_{\text{prior}}+\log w_{\text{likelihood}}+\log w_{\text{factors}}$$

This decomposition enables sophisticated inference algorithms to reason about different sources of probability mass.

Basic Trace Inspection

Let's start by understanding how Fugue records execution history:

use fugue::*;
use fugue::runtime::interpreters::PriorHandler;
use rand::{SeedableRng, rngs::StdRng};
// Demonstrate basic trace inspection and manipulation
fn basic_trace_inspection() {
    println!("=== Basic Trace Inspection ===\n");

    // Define a model with multiple types of choices
    let model = prob!(
        let coin <- sample(addr!("coin"), Bernoulli::new(0.7).unwrap());
        let count <- sample(addr!("count"), Poisson::new(3.0).unwrap());
        let category <- sample(addr!("category"), Categorical::uniform(3).unwrap());
        let measurement <- sample(addr!("measurement"), Normal::new(0.0, 1.0).unwrap());

        // Observation adds to likelihood
        let _obs <- observe(addr!("obs"), Normal::new(measurement, 0.1).unwrap(), 0.5);

        pure((coin, count, category, measurement))
    );

    // Execute and inspect the trace
    let mut rng = StdRng::seed_from_u64(12345);
    let (result, trace) = runtime::handler::run(
        PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        model,
    );

    println!("🔍 Trace Inspection:");
    println!("   - Total choices: {}", trace.choices.len());
    println!("   - Prior log-weight: {:.4}", trace.log_prior);
    println!("   - Likelihood log-weight: {:.4}", trace.log_likelihood);
    println!("   - Factor log-weight: {:.4}", trace.log_factors);
    println!("   - Total log-weight: {:.4}", trace.total_log_weight());
    println!();

    println!("📊 Individual Choices:");
    for (addr, choice) in &trace.choices {
        println!(
            "   - {}: {:?} (logp: {:.4})",
            addr, choice.value, choice.logp
        );
    }
    println!();

    println!("🎯 Type-Safe Value Access:");
    println!("   - Coin (bool): {:?}", trace.get_bool(&addr!("coin")));
    println!("   - Count (u64): {:?}", trace.get_u64(&addr!("count")));
    println!(
        "   - Category (usize): {:?}",
        trace.get_usize(&addr!("category"))
    );
    println!(
        "   - Measurement (f64): {:?}",
        trace.get_f64(&addr!("measurement"))
    );

    let (coin, count, category, measurement) = result;
    println!(
        "   - Result: coin={}, count={}, category={}, measurement={:.3}",
        coin, count, category, measurement
    );
    println!();
}

Trace Structure Analysis

Every trace contains three critical components:

1. Choices Map: Records every random decision with its address, value, and log-probability
2. Weight Decomposition: Separates prior, likelihood, and factor contributions
3. Type-Safe Values: Maintains natural types throughout execution

Replay Mechanics

The replay system is the foundation of MCMC algorithms. It allows deterministic re-execution with modified random choices:

use fugue::*;
use fugue::runtime::interpreters::*;
use fugue::runtime::trace::*;
use rand::{SeedableRng, rngs::StdRng};
// Demonstrate trace replay mechanics for MCMC
fn replay_mechanics() {
    println!("=== Trace Replay Mechanics ===\n");

    // Define a simple model
    let make_model = || {
        prob!(
            let mu <- sample(addr!("mu"), Normal::new(0.0, 2.0).unwrap());
            let _obs <- observe(addr!("y"), Normal::new(mu, 1.0).unwrap(), 1.5);
            pure(mu)
        )
    };

    // 1. Generate initial trace
    let mut rng = StdRng::seed_from_u64(42);
    let (mu1, trace1) = runtime::handler::run(
        PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        make_model(),
    );

    println!("🎲 Original Execution:");
    println!("   - mu = {:.3}", mu1);
    println!("   - Prior logp: {:.3}", trace1.log_prior);
    println!("   - Likelihood logp: {:.3}", trace1.log_likelihood);
    println!("   - Total logp: {:.3}", trace1.total_log_weight());
    println!();

    // 2. Replay with exact same trace
    let mut rng2 = StdRng::seed_from_u64(42); // New RNG for replay
    let (mu2, trace2) = runtime::handler::run(
        ReplayHandler {
            rng: &mut rng2,
            base: trace1.clone(),
            trace: Trace::default(),
        },
        make_model(),
    );

    println!("🔄 Exact Replay:");
    println!("   - mu = {:.3} (should match original)", mu2);
    println!("   - Values match: {}", mu1 == mu2);
    println!(
        "   - Traces match: {}",
        trace1.total_log_weight() == trace2.total_log_weight()
    );
    println!();

    // 3. Modify trace for proposal
    let mut modified_trace = trace1.clone();
    // Modify the mu value (MCMC proposal)
    if let Some(choice) = modified_trace.choices.get_mut(&addr!("mu")) {
        let old_value = choice.value.as_f64().unwrap();
        let new_value = old_value + 0.1; // Small proposal step
        choice.value = ChoiceValue::F64(new_value);

        // Recompute log-probability under the distribution
        let normal_dist = Normal::new(0.0, 2.0).unwrap();
        choice.logp = normal_dist.log_prob(&new_value);

        println!("🔧 Modified Trace (Proposal):");
        println!("   - Old mu: {:.3}", old_value);
        println!("   - New mu: {:.3}", new_value);
        println!(
            "   - Old logp: {:.3}",
            trace1
                .get_f64(&addr!("mu"))
                .map(|v| Normal::new(0.0, 2.0).unwrap().log_prob(&v))
                .unwrap_or(0.0)
        );
        println!("   - New logp: {:.3}", choice.logp);
    }

    // 4. Score the modified trace
    let mut rng3 = StdRng::seed_from_u64(42); // New RNG for proposal
    let (mu3, trace3) = runtime::handler::run(
        ReplayHandler {
            rng: &mut rng3,
            base: modified_trace,
            trace: Trace::default(),
        },
        make_model(),
    );

    println!("   - Proposal result: mu = {:.3}", mu3);
    println!("   - Proposal total logp: {:.3}", trace3.total_log_weight());
    println!(
        "   - Accept/Reject ratio: {:.3}",
        (trace3.total_log_weight() - trace1.total_log_weight()).exp()
    );
    println!();
}

MCMC Proposal Mechanism

sequenceDiagram
    participant M as Model
    participant T1 as Current Trace
    participant T2 as Proposal Trace
    participant A as Accept/Reject
    
    M->>T1: Execute with PriorHandler
    Note over T1: Record current state
    
    T1->>T2: Modify choice values
    Note over T2: Create proposal
    
    M->>T2: Execute with ReplayHandler  
    Note over T2: Score proposal
    
    T2->>A: Compare log-weights
    Note over A: Accept if log(u) < Δlog(w)
    
    A->>T1: Keep current (if rejected)
    A->>T2: Accept proposal (if accepted)

The key insight: the same model specification can be executed with different random choices by manipulating the trace and using replay.

Custom Handlers

Handlers define how probabilistic effects are interpreted. You can create custom handlers for specialized inference algorithms:

use fugue::*;
use fugue::runtime::{handler::Handler, trace::*};
use rand::{SeedableRng, rngs::StdRng};
// Demonstrate custom handler implementation
struct DebugHandler<R: rand::Rng> {
    rng: R,
    trace: Trace,
    debug_info: Vec<String>,
}

impl<R: rand::Rng> DebugHandler<R> {
    fn new(rng: R) -> Self {
        Self {
            rng,
            trace: Trace::default(),
            debug_info: Vec::new(),
        }
    }
}

impl<R: rand::Rng> Handler for DebugHandler<R> {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        let value = dist.sample(&mut self.rng);
        let logp = dist.log_prob(&value);

        self.debug_info.push(format!(
            "SAMPLE f64 at {}: {} (logp: {:.3})",
            addr, value, logp
        ));

        self.trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::F64(value),
                logp,
            },
        );
        self.trace.log_prior += logp;

        value
    }

    fn on_sample_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>) -> bool {
        let value = dist.sample(&mut self.rng);
        let logp = dist.log_prob(&value);

        self.debug_info.push(format!(
            "SAMPLE bool at {}: {} (logp: {:.3})",
            addr, value, logp
        ));

        self.trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::Bool(value),
                logp,
            },
        );
        self.trace.log_prior += logp;

        value
    }

    fn on_sample_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>) -> u64 {
        let value = dist.sample(&mut self.rng);
        let logp = dist.log_prob(&value);

        self.debug_info.push(format!(
            "SAMPLE u64 at {}: {} (logp: {:.3})",
            addr, value, logp
        ));

        self.trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::U64(value),
                logp,
            },
        );
        self.trace.log_prior += logp;

        value
    }

    fn on_sample_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>) -> usize {
        let value = dist.sample(&mut self.rng);
        let logp = dist.log_prob(&value);

        self.debug_info.push(format!(
            "SAMPLE usize at {}: {} (logp: {:.3})",
            addr, value, logp
        ));

        self.trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::Usize(value),
                logp,
            },
        );
        self.trace.log_prior += logp;

        value
    }

    fn on_observe_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>, value: f64) {
        let logp = dist.log_prob(&value);

        self.debug_info.push(format!(
            "OBSERVE f64 at {}: {} (logp: {:.3})",
            addr, value, logp
        ));

        self.trace.log_likelihood += logp;
    }

    fn on_observe_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>, value: bool) {
        let logp = dist.log_prob(&value);

        self.debug_info.push(format!(
            "OBSERVE bool at {}: {} (logp: {:.3})",
            addr, value, logp
        ));

        self.trace.log_likelihood += logp;
    }

    fn on_observe_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>, value: u64) {
        let logp = dist.log_prob(&value);

        self.debug_info.push(format!(
            "OBSERVE u64 at {}: {} (logp: {:.3})",
            addr, value, logp
        ));

        self.trace.log_likelihood += logp;
    }

    fn on_observe_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>, value: usize) {
        let logp = dist.log_prob(&value);

        self.debug_info.push(format!(
            "OBSERVE usize at {}: {} (logp: {:.3})",
            addr, value, logp
        ));

        self.trace.log_likelihood += logp;
    }

    fn on_factor(&mut self, logw: f64) {
        self.debug_info.push(format!("FACTOR: {:.3}", logw));
        self.trace.log_factors += logw;
    }

    fn finish(self) -> Trace {
        self.trace
    }
}

fn custom_handler_demo() {
    println!("=== Custom Handler Demo ===\n");

    let model = prob!(
        let prior_mu <- sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap());
        let success <- sample(addr!("success"), Bernoulli::new(0.6).unwrap());

        let _obs1 <- observe(addr!("data1"), Normal::new(prior_mu, 0.5).unwrap(), 1.2);
        let _obs2 <- observe(addr!("data2"), Bernoulli::new(if success { 0.8 } else { 0.2 }).unwrap(), true);

        // Add a factor for soft constraints
        let _factor_result <- factor(if prior_mu > 0.0 { 0.1 } else { -0.1 });

        pure((prior_mu, success))
    );

    let rng = StdRng::seed_from_u64(67890);
    let debug_handler = DebugHandler::new(rng);

    let (result, final_trace) = runtime::handler::run(debug_handler, model);

    println!("🔍 Debug Handler Output:");
    println!("   - Result: {:?}", result);
    println!(
        "   - Total log-weight: {:.4}",
        final_trace.total_log_weight()
    );
    println!();

    println!("📝 Execution Log:");
    println!("   - {} operations recorded", final_trace.choices.len());

    println!();
}

Handler Architecture

graph TD
    A["Model Effects"] --> B["Handler Dispatch"]
    B --> C["on_sample_f64()"]
    B --> D["on_sample_bool()"] 
    B --> E["on_sample_u64()"]
    B --> F["on_sample_usize()"]
    B --> G["on_observe_*()"]
    B --> H["on_factor()"]
    
    C --> I["Custom Logic"]
    D --> I
    E --> I
    F --> I
    G --> I
    H --> I
    
    I --> J["Trace Update"]
    I --> K["Side Effects"]
    I --> L["Logging/Debug"]
    
    style I fill:#ccffcc

Built-in Handler Types

Trace Scoring

Scoring computes the log-probability of a specific execution path, essential for importance sampling and model comparison:

use fugue::*;
use fugue::runtime::interpreters::*;
use fugue::runtime::trace::*;
use rand::{SeedableRng, rngs::StdRng};
// Demonstrate trace scoring for importance sampling
fn trace_scoring_demo() {
    println!("=== Trace Scoring Demo ===\n");

    let make_model = || {
        prob!(
            let theta <- sample(addr!("theta"), Beta::new(2.0, 2.0).unwrap());

            // Multiple observations
            let _obs1 <- observe(addr!("y1"), Bernoulli::new(theta).unwrap(), true);
            let _obs2 <- observe(addr!("y2"), Bernoulli::new(theta).unwrap(), true);
            let _obs3 <- observe(addr!("y3"), Bernoulli::new(theta).unwrap(), false);

            pure(theta)
        )
    };

    // Generate a trace from the prior
    let mut rng = StdRng::seed_from_u64(111);
    let (theta_val, prior_trace) = runtime::handler::run(
        PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        make_model(),
    );

    println!("🎲 Prior Sample:");
    println!("   - theta = {:.3}", theta_val);
    println!("   - Prior logp: {:.3}", prior_trace.log_prior);
    println!("   - Likelihood logp: {:.3}", prior_trace.log_likelihood);
    println!("   - Total logp: {:.3}", prior_trace.total_log_weight());
    println!();

    // Now score this trace under the model (should get same result)
    let (theta_scored, scored_trace) = runtime::handler::run(
        ScoreGivenTrace {
            base: prior_trace.clone(),
            trace: Trace::default(),
        },
        make_model(),
    );

    println!("📊 Scoring Same Trace:");
    println!("   - theta = {:.3} (should match)", theta_scored);
    println!("   - Prior logp: {:.3}", scored_trace.log_prior);
    println!("   - Likelihood logp: {:.3}", scored_trace.log_likelihood);
    println!("   - Total logp: {:.3}", scored_trace.total_log_weight());
    println!(
        "   - Weights match: {}",
        (prior_trace.total_log_weight() - scored_trace.total_log_weight()).abs() < 1e-10
    );
    println!();

    // Create a modified trace for importance sampling
    let mut importance_trace = prior_trace.clone();

    // Change theta to a different value
    if let Some(choice) = importance_trace.choices.get_mut(&addr!("theta")) {
        let new_theta = 0.8; // High success probability
        choice.value = ChoiceValue::F64(new_theta);
        choice.logp = Beta::new(2.0, 2.0).unwrap().log_prob(&new_theta);

        println!("🎯 Importance Sample:");
        println!("   - Modified theta to: {:.3}", new_theta);
        println!("   - New prior logp: {:.3}", choice.logp);
    }

    // Score under original model
    let (theta_is, is_trace) = runtime::handler::run(
        ScoreGivenTrace {
            base: importance_trace,
            trace: Trace::default(),
        },
        make_model(),
    );

    println!("   - IS result: theta = {:.3}", theta_is);
    println!("   - IS total logp: {:.3}", is_trace.total_log_weight());
    println!(
        "   - Importance weight: {:.3}",
        is_trace.total_log_weight() - prior_trace.total_log_weight()
    );
    println!();
}

Importance Sampling Theory

Given proposal trace ${\tau }_q$ and target model $p$:

$$w({\tau }_q)=\frac{p({\tau }_q)}{q({\tau }_q)}$$

Where the importance weight is: $$\log w=\log p({\tau }_q)-\log q({\tau }_q)$$

Fugue's scoring system automatically computes $\log p({\tau }_q)$ for any trace under any model.

Memory Optimization

For production workloads, efficient memory management is crucial:

use fugue::*;
use fugue::runtime::{interpreters::PriorHandler, memory::*};
use rand::{SeedableRng, rngs::StdRng};
// Demonstrate memory-optimized trace handling
fn memory_optimization_demo() {
    println!("=== Memory Optimization Demo ===\n");

    // Simple model for batch processing
    let make_model = |obs_val: f64| {
        prob!(
            let mu <- sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap());
            let sigma <- sample(addr!("sigma"), Gamma::new(2.0, 0.5).unwrap());
            let _obs <- observe(addr!("y"), Normal::new(mu, sigma).unwrap(), obs_val);
            pure((mu, sigma))
        )
    };

    println!("🏭 Batch Processing with Memory Pool:");

    // Simulate batch inference with trace reuse
    let observations = [1.0, 1.2, 0.8, 1.5, 0.9];
    let mut results = Vec::new();

    // Use copy-on-write traces for efficiency
    let base_trace = CowTrace::new();

    for (i, &obs) in observations.iter().enumerate() {
        let mut rng = StdRng::seed_from_u64(200 + i as u64);
        let handler = PriorHandler {
            rng: &mut rng,
            trace: base_trace.to_trace(), // Convert to regular trace
        };

        let (result, trace) = runtime::handler::run(handler, make_model(obs));
        results.push((result, trace));

        println!(
            "   Sample {}: mu={:.3}, sigma={:.3}, obs={:.1}, logp={:.3}",
            i + 1,
            result.0,
            result.1,
            obs,
            results[i].1.total_log_weight()
        );
    }

    println!();
    println!("📊 Batch Statistics:");
    let mu_mean = results.iter().map(|((mu, _), _)| mu).sum::<f64>() / results.len() as f64;
    let sigma_mean =
        results.iter().map(|((_, sigma), _)| sigma).sum::<f64>() / results.len() as f64;
    let logp_mean = results
        .iter()
        .map(|(_, trace)| trace.total_log_weight())
        .sum::<f64>()
        / results.len() as f64;

    println!("   - Average mu: {:.3}", mu_mean);
    println!("   - Average sigma: {:.3}", sigma_mean);
    println!("   - Average log-probability: {:.3}", logp_mean);
    println!();

    println!("🔧 Trace Builder Demo:");

    // Demonstrate efficient trace building
    let _builder = TraceBuilder::new();
    // Note: TraceBuilder API may not have reserve_choices method
    // This is a conceptual example of memory pre-allocation

    // Manually construct a trace (rarely needed, but shows internals)
    let demo_trace = Trace {
        choices: [
            (
                addr!("param1"),
                Choice {
                    addr: addr!("param1"),
                    value: ChoiceValue::F64(0.5),
                    logp: -1.4,
                },
            ),
            (
                addr!("param2"),
                Choice {
                    addr: addr!("param2"),
                    value: ChoiceValue::Bool(true),
                    logp: -0.7,
                },
            ),
        ]
        .iter()
        .cloned()
        .collect(),
        log_prior: -2.1,
        log_likelihood: -0.5,
        log_factors: 0.0,
    };

    println!(
        "   - Manual trace: {} choices, total logp: {:.3}",
        demo_trace.choices.len(),
        demo_trace.total_log_weight()
    );
    println!();
}

Production Memory Strategies

1. Copy-on-Write Traces: Share read-only data, copy only when modified
2. Trace Pooling: Reuse allocated memory across multiple inferences
3. Pre-sized Allocation: Reserve space for expected number of choices
4. Batch Processing: Amortize allocation costs across many executions

Diagnostic Tools

Fugue provides comprehensive tools for analyzing trace quality and convergence:

use fugue::*;
use fugue::runtime::interpreters::PriorHandler;
use fugue::inference::diagnostics::*;
use rand::{SeedableRng, rngs::StdRng};
// Demonstrate diagnostic tools for trace analysis
fn diagnostic_tools_demo() {
    println!("=== Diagnostic Tools Demo ===\n");

    // Generate multiple MCMC-like traces for diagnostics
    let make_model = || {
        prob!(
            let mu <- sample(addr!("mu"), Normal::new(0.0, 2.0).unwrap());
            let precision <- sample(addr!("precision"), Gamma::new(2.0, 1.0).unwrap());

            // Multiple observations
            let _obs1 <- observe(addr!("y1"), Normal::new(mu, 1.0/precision.sqrt()).unwrap(), 1.0);
            let _obs2 <- observe(addr!("y2"), Normal::new(mu, 1.0/precision.sqrt()).unwrap(), 1.2);
            let _obs3 <- observe(addr!("y3"), Normal::new(mu, 1.0/precision.sqrt()).unwrap(), 0.8);

            pure((mu, precision))
        )
    };

    // Simulate two chains
    println!("🔗 Generating MCMC-like traces:");
    let mut chain1 = Vec::new();
    let mut chain2 = Vec::new();

    // Chain 1
    for i in 0..20 {
        let mut rng = StdRng::seed_from_u64(300 + i);
        let (_, trace) = runtime::handler::run(
            PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            },
            make_model(),
        );
        chain1.push(trace);
    }

    // Chain 2 (different seed)
    for i in 0..20 {
        let mut rng = StdRng::seed_from_u64(400 + i);
        let (_, trace) = runtime::handler::run(
            PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            },
            make_model(),
        );
        chain2.push(trace);
    }

    println!("   - Chain 1: {} samples", chain1.len());
    println!("   - Chain 2: {} samples", chain2.len());
    println!();

    // Extract parameter values
    let mu_values1 = extract_f64_values(&chain1, &addr!("mu"));
    let mu_values2 = extract_f64_values(&chain2, &addr!("mu"));
    let precision_values1 = extract_f64_values(&chain1, &addr!("precision"));
    let _precision_values2 = extract_f64_values(&chain2, &addr!("precision"));

    println!("📈 Parameter Summaries:");
    println!(
        "   - mu chain1: mean={:.3}, min={:.3}, max={:.3}",
        mu_values1.iter().sum::<f64>() / mu_values1.len() as f64,
        mu_values1.iter().fold(f64::INFINITY, |a, &b| a.min(b)),
        mu_values1.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b))
    );

    println!(
        "   - mu chain2: mean={:.3}, min={:.3}, max={:.3}",
        mu_values2.iter().sum::<f64>() / mu_values2.len() as f64,
        mu_values2.iter().fold(f64::INFINITY, |a, &b| a.min(b)),
        mu_values2.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b))
    );

    println!(
        "   - precision chain1: mean={:.3}, min={:.3}, max={:.3}",
        precision_values1.iter().sum::<f64>() / precision_values1.len() as f64,
        precision_values1
            .iter()
            .fold(f64::INFINITY, |a, &b| a.min(b)),
        precision_values1
            .iter()
            .fold(f64::NEG_INFINITY, |a, &b| a.max(b))
    );

    // Compute R-hat (simplified version)
    let chains_mu = vec![chain1.clone(), chain2.clone()];
    let r_hat_mu = r_hat_f64(&chains_mu, &addr!("mu"));
    let r_hat_precision = r_hat_f64(&chains_mu, &addr!("precision"));

    println!();
    println!("🎯 Convergence Diagnostics:");
    println!("   - R-hat for mu: {:.4} (< 1.1 is good)", r_hat_mu);
    println!(
        "   - R-hat for precision: {:.4} (< 1.1 is good)",
        r_hat_precision
    );

    // Parameter summary
    let mu_summary = summarize_f64_parameter(&chains_mu, &addr!("mu"));
    let q5 = mu_summary.quantiles.get("2.5%").unwrap_or(&f64::NAN);
    let q95 = mu_summary.quantiles.get("97.5%").unwrap_or(&f64::NAN);
    println!(
        "   - mu summary: mean={:.3}, std={:.3}, q2.5={:.3}, q97.5={:.3}",
        mu_summary.mean, mu_summary.std, q5, q95
    );

    println!();
    println!("📊 Trace Quality Assessment:");
    let total_logp_chain1: f64 = chain1.iter().map(|t| t.total_log_weight()).sum();
    let total_logp_chain2: f64 = chain2.iter().map(|t| t.total_log_weight()).sum();
    let avg_logp1 = total_logp_chain1 / chain1.len() as f64;
    let avg_logp2 = total_logp_chain2 / chain2.len() as f64;

    println!("   - Chain 1 avg log-probability: {:.3}", avg_logp1);
    println!("   - Chain 2 avg log-probability: {:.3}", avg_logp2);
    println!(
        "   - Chains similar quality: {}",
        (avg_logp1 - avg_logp2).abs() < 0.5
    );
    println!();
}

Convergence Assessment

R-hat Statistic: Compares between-chain variance to within-chain variance

$$\hat{R}=\sqrt{\frac{\hat{V}}{W}}$$

Where:

• $\hat{V}$: Estimated marginal posterior variance
• $W$: Within-chain variance

Interpretation:

• $\hat{R}\approx 1.0$: Good convergence
• $\hat{R}>1.1$: Chains haven't mixed well, need more samples
• $\hat{R}>1.2$: Poor convergence, investigate model or algorithm

Parameter Summaries

For each parameter, compute:

• Mean and Standard Deviation: Central tendency and spread
• Quantiles: 5%, 25%, 50%, 75%, 95% for uncertainty intervals
• Effective Sample Size: Accounting for autocorrelation

Advanced Debugging

When models behave unexpectedly, trace analysis reveals the root causes:

use fugue::*;
use fugue::runtime::interpreters::*;
use fugue::runtime::trace::*;
use rand::{SeedableRng, rngs::StdRng};
// Demonstrate advanced debugging techniques
fn advanced_debugging_demo() {
    println!("=== Advanced Debugging Techniques ===\n");

    // Model with potential numerical issues
    let _problematic_model = prob!(
        let scale <- sample(addr!("scale"), Exponential::new(1.0).unwrap());

        // This could cause numerical issues if scale is very small
        let precision <- sample(addr!("precision"), Gamma::new(1.0, scale).unwrap());

        let mu <- sample(addr!("mu"), Normal::new(0.0, 1.0 / precision.sqrt()).unwrap());

        // Observation that might conflict
        let _obs <- observe(addr!("y"), Normal::new(mu, 0.01).unwrap(), 10.0);

        pure((scale, precision, mu))
    );

    println!("🚨 Debugging Problematic Model:");

    // Try multiple executions to find issues
    for attempt in 1..=5 {
        let mut rng = StdRng::seed_from_u64(500 + attempt);
        let problematic_model_copy = prob!(
            let scale <- sample(addr!("scale"), Exponential::new(1.0).unwrap());

            // This could cause numerical issues if scale is very small
            let precision <- sample(addr!("precision"), Gamma::new(1.0, scale).unwrap());

            let mu <- sample(addr!("mu"), Normal::new(0.0, 1.0 / precision.sqrt()).unwrap());

            // Observation that might conflict
            let _obs <- observe(addr!("y"), Normal::new(mu, 0.01).unwrap(), 10.0);

            pure((scale, precision, mu))
        );

        let (result, trace) = runtime::handler::run(
            PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            },
            problematic_model_copy,
        );

        let (scale, precision, mu) = result;
        let total_logp = trace.total_log_weight();

        println!(
            "   Attempt {}: scale={:.6}, precision={:.6}, mu={:.3}, logp={:.3}",
            attempt, scale, precision, mu, total_logp
        );

        // Check for numerical issues
        if !total_logp.is_finite() {
            println!("     ⚠️ Non-finite log-probability detected!");
        }

        if precision < 1e-6 {
            println!("     ⚠️ Very small precision: {:.8}", precision);
        }

        if mu.abs() > 5.0 {
            println!("     ⚠️ Extreme mu value: {:.3}", mu);
        }

        // Examine individual components
        println!(
            "     Components: prior={:.3}, likelihood={:.3}, factors={:.3}",
            trace.log_prior, trace.log_likelihood, trace.log_factors
        );
    }

    println!();
    println!("🔍 Trace Validation:");

    // Create a trace with known good values for validation
    let validation_trace = Trace {
        choices: [
            (
                addr!("scale"),
                Choice {
                    addr: addr!("scale"),
                    value: ChoiceValue::F64(1.0),
                    logp: Exponential::new(1.0).unwrap().log_prob(&1.0),
                },
            ),
            (
                addr!("precision"),
                Choice {
                    addr: addr!("precision"),
                    value: ChoiceValue::F64(2.0),
                    logp: Gamma::new(1.0, 1.0).unwrap().log_prob(&2.0),
                },
            ),
            (
                addr!("mu"),
                Choice {
                    addr: addr!("mu"),
                    value: ChoiceValue::F64(0.5),
                    logp: Normal::new(0.0, 1.0 / (2.0_f64).sqrt())
                        .unwrap()
                        .log_prob(&0.5),
                },
            ),
        ]
        .iter()
        .cloned()
        .collect(),
        log_prior: 0.0,
        log_likelihood: 0.0,
        log_factors: 0.0,
    };

    // Score this validation trace
    let validation_model = prob!(
        let scale <- sample(addr!("scale"), Exponential::new(1.0).unwrap());
        let precision <- sample(addr!("precision"), Gamma::new(1.0, scale).unwrap());
        let mu <- sample(addr!("mu"), Normal::new(0.0, 1.0 / precision.sqrt()).unwrap());
        let _obs <- observe(addr!("y"), Normal::new(mu, 0.01).unwrap(), 10.0);
        pure((scale, precision, mu))
    );

    let (val_result, val_trace) = runtime::handler::run(
        ScoreGivenTrace {
            base: validation_trace,
            trace: Trace::default(),
        },
        validation_model,
    );

    println!("   - Validation result: {:?}", val_result);
    println!("   - Validation logp: {:.3}", val_trace.total_log_weight());
    println!(
        "   - Validation finite: {}",
        val_trace.total_log_weight().is_finite()
    );
    println!();
}

Debugging Strategy

graph TD
    A["Model Issues"] --> B["Check Log-Weights"]
    B --> C["Infinite/NaN Values?"]
    B --> D["Extreme Parameter Values?"]
    B --> E["Prior-Likelihood Conflict?"]
    
    C --> F["Examine Individual Choices"]
    D --> G["Check Parameter Ranges"]
    E --> H["Validate Observations"]
    
    F --> I["Fix Distribution Parameters"]
    G --> J["Add Constraints/Priors"] 
    H --> K["Verify Data Consistency"]
    
    I --> L["Test with Validation Trace"]
    J --> L
    K --> L
    
    style C fill:#ffcccc
    style D fill:#ffcccc
    style E fill:#ffcccc

Common Issues and Solutions

Real-World Applications

Custom MCMC Algorithm

use fugue::*;
use fugue::runtime::{interpreters::*, trace::*};

struct CustomMCMC<R: rand::Rng> {
    rng: R,
    current_trace: Trace,
    step_size: f64,
}

impl<R: rand::Rng> CustomMCMC<R> {
    fn step<F>(&mut self, model_fn: F) -> bool 
    where F: Fn() -> Model<f64>
    {
        // Create proposal by modifying current trace
        let mut proposal_trace = self.current_trace.clone();
        
        // Modify a random choice (simplified)
        if let Some((addr, choice)) = proposal_trace.choices.iter_mut().next() {
            if let Some(current_val) = choice.value.as_f64() {
                let proposal_val = current_val + self.step_size * 
                    Normal::new(0.0, 1.0).unwrap().sample(&mut self.rng);
                choice.value = ChoiceValue::F64(proposal_val);
            }
        }
        
        // Score proposal
        let (_, scored_trace) = runtime::handler::run(
            ScoreGivenTrace::new(proposal_trace),
            model_fn()
        );
        
        // Accept/reject based on Metropolis criterion
        let log_alpha = scored_trace.total_log_weight() - 
                       self.current_trace.total_log_weight();
        
        if log_alpha > 0.0 || 
           self.rng.gen::<f64>().ln() < log_alpha {
            self.current_trace = scored_trace;
            true // Accepted
        } else {
            false // Rejected  
        }
    }
}

Production Inference Pipeline

use fugue::*;
use fugue::runtime::memory::TracePool;

struct InferencePipeline {
    pool: TracePool,
    diagnostics: Vec<f64>,
}

impl InferencePipeline {
    fn run_batch<F>(&mut self, 
                   model_fn: F, 
                   n_samples: usize) -> Vec<(f64, Trace)>
    where F: Fn() -> Model<f64> + Copy
    {
        let mut results = Vec::with_capacity(n_samples);
        
        for _ in 0..n_samples {
            // Get pooled trace to avoid allocation
            let pooled_trace = self.pool.get_trace();
            
            let mut rng = rand::thread_rng();
            let handler = PriorHandler { 
                rng: &mut rng, 
                trace: pooled_trace 
            };
            
            let (result, trace) = runtime::handler::run(handler, model_fn());
            
            // Record diagnostics
            self.diagnostics.push(trace.total_log_weight());
            
            results.push((result, trace));
        }
        
        results
    }
    
    fn convergence_summary(&self) -> (f64, f64) {
        let mean = self.diagnostics.iter().sum::<f64>() / self.diagnostics.len() as f64;
        let var = self.diagnostics.iter()
            .map(|x| (x - mean).powi(2))
            .sum::<f64>() / (self.diagnostics.len() - 1) as f64;
        (mean, var.sqrt())
    }
}

Best Practices

Handler Development

Memory Management

Debugging Workflow

Testing Your Understanding

Exercise 1: Custom Proposal Mechanism

Implement a custom handler that uses adaptive proposals based on the acceptance rate history:

use fugue::*;
use fugue::runtime::{handler::Handler, trace::*};

struct AdaptiveMCMCHandler<R: rand::Rng> {
    rng: R,
    current_trace: Trace,
    proposal_scale: f64,
    acceptance_history: Vec<bool>,
    adaptation_interval: usize,
}

// TODO: Implement Handler trait with adaptive step size

Exercise 2: Multi-Chain Diagnostics

Create a system that runs multiple MCMC chains in parallel and automatically assesses convergence:

use fugue::inference::diagnostics::*;

fn multi_chain_inference<F>(
    model_fn: F,
    n_chains: usize, 
    n_samples: usize
) -> (Vec<Vec<Trace>>, bool)
where F: Fn() -> Model<f64> + Copy
{
    // TODO: Run multiple chains and check R-hat convergence
    unimplemented!()
}

Exercise 3: Memory-Optimized Batch Processing

Design a system for processing thousands of similar models efficiently:

use fugue::runtime::memory::*;

struct BatchProcessor {
    pool: TracePool,
    // TODO: Add fields for efficient batch processing
}

impl BatchProcessor {
    fn process_batch<F>(&mut self, 
                       models: Vec<F>) -> Vec<(f64, Trace)>
    where F: Fn() -> Model<f64>
    {
        // TODO: Implement memory-efficient batch processing
        unimplemented!()
    }
}

Key Takeaways

Core Capabilities:

• ✅ Complete execution recording with type safety and weight decomposition
• ✅ Flexible interpretation through the handler system
• ✅ MCMC foundation via deterministic replay mechanics
• ✅ Custom inference algorithms through handler extensibility
• ✅ Production optimization with memory pooling and efficient allocation
• ✅ Comprehensive diagnostics for convergence assessment and debugging

Further Reading

• Custom Handlers Guide - Building specialized interpreters
• Optimizing Performance - Production deployment strategies
• Debugging Models - Troubleshooting problematic models
• API Reference - Complete runtime system specification
• The Elements of Statistical Learning - Theoretical foundations of inference algorithms
• Monte Carlo Statistical Methods - MCMC theory and practice

Statistical Modeling