Understanding and Exploring Models

This tutorial covers how to work with PyHS3 models - understanding their structure, exploring their contents, and evaluating them.

What is a Model?

A Model is the computational representation of your statistical model created from a workspace. It contains:

  • Parameters: Symbolic tensor variables representing your model parameters

  • Distributions: Compiled probability distribution functions

  • Functions: Compiled mathematical functions

  • Computational Graph: The dependency structure between all components

Models are created from workspaces and provide the interface for evaluating PDFs, generating samples, and performing statistical analysis.

Creating Models

Models are created from workspaces using the .model() method. The following example uses the :hs3:label:`gaussian_dist <hs3.gaussian-normal-distribution>` distribution type:

Exploring Model Structure

Once you have a model, you can explore its structure:

Understanding the Computational Graph

PyHS3 models are built as computational graphs where:

  • Parameters are leaf nodes (input variables)

  • Functions transform parameters into intermediate values

  • Distributions depend on parameters and/or function outputs

  • Dependencies define the evaluation order

You can visualize the computational graph:

Parameter Discovery and Bounds

PyHS3 automatically discovers parameters from your distributions and functions. Parameters are created with domain bounds applied:

Evaluating Models

The primary use of models is to evaluate probability density functions. PyHS3 provides both type-safe and convenience methods for PDF evaluation:

Convenience API (For Testing and Interactive Use)

The pdf_unsafe() and logpdf_unsafe() methods automatically convert Python floats and lists to numpy arrays. Use these for convenience in testing or interactive sessions:

# Floats are automatically converted to numpy arrays
pdf_value = model.pdf_unsafe("gaussian_model", observable=0.0, mu=0.0, sigma=1.0)
print(f"PDF(0.0) = {pdf_value}")

# Lists are also supported (converted to 1-d arrays)
pdf_values = model.pdf_unsafe("gaussian_model", observable=1.5, mu=0.0, sigma=1.0)

# Log PDF with automatic conversion
logpdf_value = model.logpdf_unsafe("gaussian_model", observable=0.0, mu=0.0, sigma=1.0)

Warning

The *_unsafe methods have a small performance overhead due to type conversion. For production code or performance-critical loops, prefer pdf() and logpdf() with pre-converted numpy arrays.

Vectorized Evaluation

For evaluating PDFs across multiple points, use numpy arrays:

import numpy as np

# Evaluate at multiple points
x_values = np.linspace(-3, 3, 100)
pdf_values = [
    model.pdf(
        "gaussian_model", observable=np.array(x), mu=np.array(0.0), sigma=np.array(1.0)
    )
    for x in x_values
]

Model Compilation and Performance

Models use PyTensor for fast compilation and evaluation:

# Models support different compilation modes
fast_model = ws.model(mode="FAST_RUN")  # Maximum optimization
debug_model = ws.model(mode="FAST_COMPILE")  # Faster compilation

# Check compilation status
print(f"Model mode: {model.mode}")
summary = model.graph_summary("gaussian_model")
print("Compiled:" in summary)  # Shows if function is compiled

Working with Complex Models

For models with multiple distributions and functions, including :hs3:label:`generic_dist <hs3.sec:generic_expression>` with mathematical expressions and :hs3:label:`sum functions <hs3.sum>`:

complex_model = {
    "metadata": {"hs3_version": "0.2"},
    "distributions": [
        {
            "name": "signal",
            "type": "gaussian_dist",
            "x": "mass",
            "mean": "signal_mean",
            "sigma": "resolution",
        },
        {
            "name": "background",
            "type": "generic_dist",
            "x": "mass",
            "expression": "exp(-mass/slope)",
        },
    ],
    "functions": [
        {
            "name": "total_yield",
            "type": "sum",
            "summands": ["signal_events", "background_events"],
        },
        {
            "name": "signal_fraction",
            "type": "generic_function",
            "expression": "signal_events / total_yield",
        },
    ],
    "parameter_points": [
        {
            "name": "physics_point",
            "parameters": [
                {"name": "signal_mean", "value": 125.0},
                {"name": "resolution", "value": 2.5},
                {"name": "signal_events", "value": 100.0},
                {"name": "background_events", "value": 1000.0},
                {"name": "slope", "value": 50.0},
            ],
        }
    ],
}

complex_ws = pyhs3.Workspace(**complex_model)
complex_model = complex_ws.model()

# Evaluate individual components
signal_pdf = complex_model.pdf("signal", mass=125.0, signal_mean=125.0, resolution=2.5)
background_pdf = complex_model.pdf("background", mass=125.0, slope=50.0)

# Evaluate functions
total = complex_model.pdf("total_yield", signal_events=100.0, background_events=1000.0)
fraction = complex_model.pdf(
    "signal_fraction", signal_events=100.0, background_events=1000.0
)

print(f"Signal PDF: {signal_pdf}")
print(f"Background PDF: {background_pdf}")
print(f"Total yield: {total}")
print(f"Signal fraction: {fraction}")

Debugging and Troubleshooting

When working with models, you can debug issues using:

# 1. Check model structure
print(model)

# 2. Examine computational graph
summary = model.graph_summary("distribution_name")
print(summary)

# 3. Use debug compilation mode
debug_model = ws.model(mode="DebugMode")

# 4. Visualize dependencies
try:
    model.visualize_graph("distribution_name")
except ImportError:
    print("Install pydot for graph visualization")

# 5. Check parameter discovery
print("Available parameters:", list(model.parameters.keys()))
print("Available distributions:", list(model.distributions.keys()))
print("Available functions:", list(model.functions.keys()))

Advanced Topics

Tensor Types

Parameters can have different tensor types based on their intended use:

# In parameter_points, you can specify tensor kinds:
vector_params = {
    "name": "vector_params",
    "parameters": [
        {"name": "scalar_param", "value": 1.0, "kind": "scalar"},  # Default
        {"name": "vector_param", "value": [1.0, 2.0], "kind": "vector"},
    ],
}

Custom Functions

You can define custom mathematical expressions:

custom_function = {
    "name": "custom_calc",
    "type": "generic_function",
    "expression": "sqrt(x**2 + y**2)",  # Uses SymPy syntax
}

Performance Optimization

For better performance:

  • Use mode="FAST_RUN" for production models

  • Avoid repeated model creation

  • Cache compiled functions when possible

  • Use appropriate tensor types for your data

# Good: reuse model
model = ws.model(mode="FAST_RUN")
results = []
for x in data_points:
    result = model.pdf("my_dist", observable=x, mu=0.0, sigma=1.0)
    results.append(result)

# Less efficient: recreate model each time
# for x in data_points:
#     model = ws.model()  # Don't do this
#     result = model.pdf("my_dist", observable=x, mu=0.0, sigma=1.0)