### Install gpcam

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_NonEuclideanInputSpaces.ipynb

Installs the latest version of the gpcam library. This is a prerequisite for running the example.

```python
# Install the newest version of gpcam
#!pip install gpcam==8.4.0
```

--------------------------------

### Install and Test gpCAM from Source

Source: https://github.com/lbl-camera/gpcam/blob/master/CLAUDE.md

Install the package in editable mode with test dependencies. Run the full test suite or a specific test. Build the package using hatch.

```bash
pip install -e .[tests]
pytest tests
pytest tests --cov=./ --cov-report=xml
pytest tests/test_gpCAM.py::test_basic_1task
hatch build
```

--------------------------------

### Install gpCAM and Plotting Tools

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_MultiTaskTest.ipynb

Installs the specified version of gpCAM and necessary plotting libraries. Run this before proceeding with other examples.

```python
## First, install the newest version of gpcam and some plotting tools
#!pip install gpCAM==8.4.0
#!pip install matplotlib
#!pip install plotly
```

--------------------------------

### Install gpcam and matplotlib

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_Minimal.ipynb

Installs specific versions of gpcam and matplotlib required for the example. Run these commands in your environment.

```python
#!pip install gpcam==8.4.0
#!pip install matplotlib
```

--------------------------------

### Basic Setup with Default Kernel

Source: https://github.com/lbl-camera/gpcam/blob/master/skills/multi-task-advanced/SKILL.md

Initializes an fvGPOptimizer with default settings, suitable for multi-task GPs. It uses a built-in deep kernel and does not require explicit hyperparameter bounds.

```python
import numpy as np
from gpcam import fvGPOptimizer

# 100 input points, 5 output channels
x_data = np.random.uniform(0, 1, (100, 2))  # 2D input
y_data = np.random.randn(100, 5)             # 5 outputs per point

# Default path — uses a built-in deep kernel, no hyperparameter bounds required:
gpo = fvGPOptimizer(x_data, y_data)
gpo.train(max_iter=20)
```

--------------------------------

### Initial Suggestion Example

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/1dSingleTaskAcqFuncTest.ipynb

This output shows an initial suggestion from the ask() function, likely before any optimization has occurred or with default parameters.

```text
Output:
led to new suggestion: 
 {'x': array([[0.]]), 'f_a(x)': array([-31.19240075]), 'opt_obj': <hgdl.hgdl.HGDL object at 0x7f325e629010>}
```

--------------------------------

### Install gpCAM Plugin Marketplace

Source: https://github.com/lbl-camera/gpcam/blob/master/docs/source/claude-skills.md

Registers the gpCAM repository as a plugin marketplace and installs the gpcam plugin. This allows Claude Code to access gpCAM skills without local cloning.

```text
/plugin marketplace add lbl-camera/gpCAM
/plugin install gpcam@gpcam
```

--------------------------------

### HPC Setup with SLURM and Dask

Source: https://github.com/lbl-camera/gpcam/blob/master/skills/gp2scale-advanced/SKILL.md

Configures a Dask SLURM cluster for high-performance computing environments. This example sets up cluster parameters like cores, memory, and walltime, then scales the jobs and creates a Dask client.

```python
from dask_jobqueue import SLURMCluster
from distributed import Client

cluster = SLURMCluster(
    cores=32,
    memory="64GB",
    walltime="01:00:00",
)
cluster.scale(jobs=4)  # 4 nodes × 32 cores
client = Client(cluster)
```

--------------------------------

### Basic gp2Scale Setup with Dask

Source: https://github.com/lbl-camera/gpcam/blob/master/skills/gp2scale-advanced/SKILL.md

Initializes a GPOptimizer for large-scale experiments using gp2Scale mode and a Dask distributed client. It's recommended to wait for workers before constructing the optimizer.

```python
from distributed import Client
from gpcam import GPOptimizer

# Start a local Dask cluster
client = Client()                  # uses all available cores
client.wait_for_workers(4)         # good practice: wait for workers before constructing

gpo = GPOptimizer(
    x_data=x_data,
    y_data=y_data,
    gp2Scale=True,
    dask_client=client,
    gp2Scale_batch_size=500,       # typical; tune up for large clusters
    init_hyperparameters=np.array([0.73, 0.0014]),  # signal var, length scale
)

gpo.train(hyperparameter_bounds=hps_bounds, max_iter=25, info=True)
```

--------------------------------

### Install Required Libraries

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_gp2ScaleTest.ipynb

Installs the necessary libraries: gpcam, matplotlib, and imate. Ensure you have gpcam version 8.4.0 or later.

```python
## First, install the newest version of gpcam, matplotlib, and imate for random linear algebra
#!pip install gpcam==8.4.0
#!pip install matplotlib
#!pip install imate
```

--------------------------------

### Example Hyperparameter Layout and Bounds

Source: https://github.com/lbl-camera/gpcam/blob/master/skills/prior-mean-functions/SKILL.md

Illustrates a common layout for coordinating hyperparameters between kernel, mean, and noise functions. It shows which indices are used for each component and provides example bounds for each hyperparameter.

```python
# Example layout:
# hps[0]     = signal variance     (kernel)
# hps[1:3]   = length scales       (kernel, 2D input)
# hps[3]     = intercept           (mean function)
# hps[4:6]   = slopes              (mean function)
# hps[6]     = noise amplitude     (noise function)
#
# Total: 7 hyperparameters

hp_bounds = np.array([
    [0.001, 100.0],   # signal variance
    [0.01, 50.0],     # length scale dim 0
    [0.01, 50.0],     # length scale dim 1
    [-10.0, 10.0],    # intercept
    [-5.0, 5.0],      # slope dim 0
    [-5.0, 5.0],      # slope dim 1
    [0.001, 10.0],    # noise amplitude
])
```

--------------------------------

### Install gpcam and matplotlib

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_Optimization.ipynb

Installs the necessary libraries for using gpCAM and plotting results. Ensure you have the latest versions.

```python
## First, install the latest version of gpCAM and matplotlib
#!pip install gpcam==8.4.0
#!pip install matplotlib
```

--------------------------------

### Asynchronous Point Suggestion with GPOptimizer

Source: https://github.com/lbl-camera/gpcam/blob/master/skills/experiment-designer/SKILL.md

Start a background search for the next optimal point using `method="hgdlAsync"`. The returned `opt_obj` can be used with `kill_client()` once the suggestion is utilized.

```python
opt_obj = gpo.ask(..., method="hgdlAsync")
# ... use suggestion ...
gpo.kill_client(opt_obj)
```

--------------------------------

### Custom Kernel Setup for fvGPOptimizer

Source: https://github.com/lbl-camera/gpcam/blob/master/skills/multi-task-advanced/SKILL.md

Configures an fvGPOptimizer with a custom kernel function and explicitly defines hyperparameter initialization and bounds. This allows for fine-grained control over the GP model.

```python
gpo = fvGPOptimizer(
    x_data=x_data,
    y_data=y_data,
    init_hyperparameters=np.ones(4) / 10.0,
    kernel_function=my_multi_task_kernel,
)
gpo.train(hyperparameter_bounds=np.array([
    [0.01, 10.0],  # signal variance
    [0.01, 10.0],  # length scale dim 0
    [0.01, 10.0],  # length scale dim 1
    [0.01, 10.0],  # length scale for task dimension
]))
```

--------------------------------

### Hyperparameter Layout Example

Source: https://github.com/lbl-camera/gpcam/blob/master/skills/noise-functions/SKILL.md

Illustrates a typical layout for hyperparameters, showing how noise function hyperparameters fit within the overall vector shared with kernel and mean parameters.

```python
# Example layout:
# hps[0]     = signal variance     (kernel)
# hps[1:3]   = length scales       (kernel)
# hps[3]     = noise std dev       (noise function) ← YOUR NOISE HP
#
# Total: 4 hyperparameters

hp_bounds = np.array([
    [0.001, 100.0],  # signal variance
    [0.01, 50.0],    # length scale dim 0
    [0.01, 50.0],    # length scale dim 1
    [0.001, 10.0],   # noise std dev
])

```

--------------------------------

### Import necessary libraries

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/1dSingleTaskAcqFuncTest.ipynb

Imports core libraries including numpy, matplotlib, GPOptimizer, time, loguru, and distributed client. This setup is required for most operations.

```python
import numpy as np
import matplotlib.pyplot as plt
from gpcam import GPOptimizer
import time
from loguru import logger
from distributed import Client
client = Client()


%load_ext autoreload
%autoreload 2
```

--------------------------------

### Initialize and Train GPOptimizer with Multiple GPs

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_SingleTaskTest.ipynb

Initializes the GPOptimizer with duplicated data to simulate running multiple Gaussian Processes in parallel. It then trains the hyperparameters using MCMC. Ensure 'st' (start time) is defined before this snippet.

```python
#duplicate data: in practice, this would be different data in every column
y_data = np.broadcast_to(y_data[:, None], (y_data.size, 10))

my_gpo = GPOptimizer(x_data,y_data,
            init_hyperparameters = np.ones((2))/10.,  # we need enough of those for kernel, noise, and prior mean functions
            noise_variances=np.ones(y_data.shape[0]) * 0.1, # providing noise variances and a noise function will raise a warning 
            compute_device='cpu',
            )


hps_bounds = np.array([[0.01,10.], #signal variance for the kernel
                       [0.01,10.], #length scale for the kernel
                      ])

print("Standard Training (MCMC)")
hps = my_gpo.train(hyperparameter_bounds=hps_bounds, info = True, max_iter = 100)
print("Result=", hps, "after ", time.time() - st, " seconds")
print("")
```

--------------------------------

### MCMC Setup for Hyperparameter Tuning

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_gp2ScaleTest.ipynb

Sets up a Markov Chain Monte Carlo (MCMC) sampler using `gpMCMC` to tune hyperparameters. It defines proposal distributions and a prior function to constrain hyperparameters within bounds.

```python
from gpcam import gpMCMC, ProposalDistribution

def in_bounds(v,bounds):
    if any(v<bounds[:,0]) or any(v>bounds[:,1]): return False
    return True

def prior_function(theta, bounds, args):
    if in_bounds(theta, bounds): return 0.
    else: return -np.inf

def func(hps,args):
  result = my_gp2S.log_likelihood(hyperparameters=hps)
  return result


pd1 = ProposalDistribution([0,1], init_prop_Sigma=np.identity(2) * 0.01)
my_mcmc = gpMCMC(func, bounds = hps_bounds, prior_function=prior_function, proposal_distributions= [pd1])

mcmc_result = my_mcmc.run_mcmc(x0=np.array([1.,0.01]), info=True, n_updates=2) ##usually you would do >200 updates
```

--------------------------------

### Basic gpCAM API Usage

Source: https://github.com/lbl-camera/gpcam/blob/master/README.md

Demonstrates initializing a GPOptimizer, training it, and iteratively asking for new points and telling the results. Includes optional training at specific intervals.

```python
!pip install gpcam

from gpCAM import GPOptimizer

my_gp = GPOptimizer(x_data,y_data,)
my_gp.train()

train_at = [10,20,30] #optional
for i in range(100):
    new = my_gp.ask(np.array([[0.,1.]]))["x"]
    my_gp.tell(new, f1(new).reshape(len(new)))
    if i in train_at: my_gp.train()

```

--------------------------------

### Testing 'ask()' with Different Optimization Methods

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/1dSingleTaskAcqFuncTest.ipynb

Demonstrates how to use the 'ask()' method to find new suggestions using various acquisition functions and optimization methods ('global', 'local', 'hgdl'). Requires initialized 'my_gpo', 'client', and defined 'bounds'.

```python
#let's test the asks:
bounds = np.array([[0.0,1.0]])
for acq_func in acquisition_functions:
    for method in ["global","local","hgdl"]:
        print("Acquisition function ", acq_func," and method ",method)
        new_suggestion = my_gpo.ask(bounds, acquisition_function=acq_func, 
                                    method=method, max_iter = 2, dask_client=client)
        print("led to new suggestion: 
", new_suggestion)
        print("")
```

--------------------------------

### Optimizer Initialization with Fixed Noise

Source: https://github.com/lbl-camera/gpcam/blob/master/skills/noise-functions/SKILL.md

Demonstrates initializing `GPOptimizer` with fixed, known noise variances using the `noise_variances` argument.

```python
# Option A: fixed known noise
gpo = GPOptimizer(x_data, y_data, noise_variances=np.full(N, 0.01))

```

--------------------------------

### Initialize GPOptimizer and Train with Different Methods

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_SingleTaskTest.ipynb

Initializes the GPOptimizer with specified data and hyperparameters, then trains it using MCMC, ADAM, global, local, and HGDL methods. Ensure necessary imports and data (x_data, y_data, skernel, meanf, f1, client) are available.

```python
my_gpo = GPOptimizer(x_data,y_data,
            init_hyperparameters = np.ones((4))/10.,  # We need enough of those for kernel, noise, and prior mean functions
            noise_variances=np.ones(y_data.shape) * 0.01, # providing noise variances and a noise function will raise a warning 
            compute_device='cpu', 
            kernel_function=skernel, 
            kernel_function_grad=None, 
            prior_mean_function=meanf, 
            prior_mean_function_grad=None,
            gp2Scale = False, 
            ram_economy=False, 
            args=None,
            )

hps_bounds = np.array([[0.01,10.], #signal variance for the kernel
                       [0.01,10.], #length scale for the kernel
                       [0.001,1.],  #noise
                       [0.01,1.]  #mean
                      ])

x_update = np.array([0.1,0.2,0.5]).reshape(3,1)
y_update = f1(x_update[:,0]) + (np.random.rand(len(x_update))-0.5) * 0.5
my_gpo.tell(x_update, 
            y_update, 
            noise_variances=np.ones(y_update.shape) * 0.05,
            append=True, rank_n_update=True)


st = time.time()
print("Standard Training (MCMC)")
hps = my_gpo.train(hyperparameter_bounds=hps_bounds, info = False)
print("Result=", hps, "after ", time.time() - st, " seconds")
print("ML: ",my_gpo.log_likelihood())
print("")

print("ADAM")
hps = my_gpo.train(hyperparameter_bounds=hps_bounds, info = True, max_iter = 100, method="adam")
print("Result=", hps, "after ", time.time() - st, " seconds")
print("ML: ",my_gpo.log_likelihood())
print("")

print("Global Training")
hps = my_gpo.train(hyperparameter_bounds=hps_bounds, method='global', max_iter = 20)
print("Result=", hps, "after ", time.time() - st, " seconds")
print("ML: ",my_gpo.log_likelihood())
print("")

print("Local Training")
hps = my_gpo.train(hyperparameter_bounds=hps_bounds, method='local')
print("Result=", hps, "after ", time.time() - st, " seconds")
print("ML: ",my_gpo.log_likelihood())
print("")

print("HGDL Training")
hps = my_gpo.train(hyperparameter_bounds=hps_bounds, method='hgdl', max_iter=2, dask_client=client)
print("Result=", hps, "after ", time.time() - st, " seconds")
print("ML: ",my_gpo.log_likelihood())
print("")

```

--------------------------------

### Initialize, train, and predict with GP

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_NonEuclideanInputSpaces.ipynb

Initializes a GPOptimizer with custom data and kernel, sets hyperparameter bounds, trains the model, and prints the learned hyperparameters, prediction, and uncertainty.

```python
my_gp = GPOptimizer(x_data,y_data, init_hyperparameters=np.ones((2)),
                    kernel_function=kernel)
bounds = np.array([[0.001,100.],[0.001,100]])
my_gp.train(hyperparameter_bounds=bounds)

print("hyperparameters: ", my_gp.hyperparameters)
print("prediction : ",my_gp.posterior_mean(['full'])["m(x)"])
print("uncertainty: ",np.sqrt(my_gp.posterior_covariance(['full'])["v(x)"])))
```

--------------------------------

### Define a custom mean function

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/1dSingleTaskAcqFuncTest.ipynb

Defines a custom mean function for the Gaussian Process. This example implements a function based on sine and a hyperparameter.

```python
def meanf(x, hps):
    #This is a simple mean function but it can be arbitrarily complex using many hyperparameters.
    return 1.-np.sin(hps[3] * x[:,0])
```

--------------------------------

### Define a custom noise function

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/1dSingleTaskAcqFuncTest.ipynb

Defines a custom noise function for the Gaussian Process. This example returns a constant noise level based on a hyperparameter.

```python
def my_noise(x,hps):
    #This is a simple noise function but can be made arbitrarily complex using many hyperparameters.
    #The noise function can return a matrix or a vector 
    return np.zeros((len(x))) + hps[2]
```

--------------------------------

### Initialize Dask Client and GPOptimizer

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_gp2ScaleTest.ipynb

Sets up a Dask distributed client and initializes the GPOptimizer with sample data. It's crucial to wait for workers to be ready before proceeding.

```python
import numpy as np
import matplotlib.pyplot as plt
from gpcam import GPOptimizer
from dask.distributed import Client
%load_ext autoreload
%autoreload 2

client = Client() ##this is the client you can make locally like this or 
#your HPC team can provide a script to get it. We included an example to get gp2Scale going
#on NERSC's Perlmutter


#It's good practice to make sure to wait for all the workers to be ready
client.wait_for_workers(4)
```

--------------------------------

### Define a custom stationary kernel function

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/1dSingleTaskAcqFuncTest.ipynb

Defines a custom stationary kernel function using `gpcam.kernels`. This example uses the matern_kernel_diff1 and allows for custom hyperparameters.

```python
#stationary
from gpcam.kernels import *
def skernel(x1,x2,hps):
    #The kernel follows the mathematical definition of a kernel. This
    #means there is no limit to the variety of kernels you can define.
    d = get_distance_matrix(x1,x2)
    return hps[0] * matern_kernel_diff1(d,hps[1])
```

--------------------------------

### Initialize and Train GPOptimizer with Custom Kernel

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_MultiTaskTest.ipynb

Initializes a GPOptimizer with prepared data and a custom kernel function. It then trains the model for a few iterations, demonstrating the integration of custom kernels.

```python
my_gp2 = fvGPOptimizer(x_data3,y_data3,)
print("Global Training in progress")
my_gp2.train(max_iter = 2)
```

--------------------------------

### Get Raw Posterior Samples

Source: https://github.com/lbl-camera/gpcam/blob/master/skills/transformed-optimizers-advanced/SKILL.md

Pass return_samples=True to obtain an array of original-scale posterior draws. This is useful for computing expectations of arbitrary functions or generating histograms.

```python
post = gpo.evaluate_posterior(x_query, return_samples=True, n_samples=8000)
samples = post["samples"]            # shape (len(x_query), 8000)

# Histogram at the first query point:
import matplotlib.pyplot as plt
plt.hist(samples[0], bins=50, density=True)
```

--------------------------------

### GPOptimizer Initialization and Data Update

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/1dSingleTaskAcqFuncTest.ipynb

Initializes GPOptimizer with custom functions and hyperparameters. Demonstrates appending and overwriting data using the 'tell' method.

```python
my_gpo = GPOptimizer(x_data,y_data,
            init_hyperparameters = np.ones((4))/10.,  # We need enough of those for kernel, noise, and prior mean functions 
            compute_device='cpu', 
            kernel_function=skernel, 
            kernel_function_grad=None, 
            prior_mean_function=meanf, 
            prior_mean_function_grad=None,
            noise_function=my_noise,
            #noise_variances=np.zeros(y_data.shape) + 0.1,
            gp2Scale = False, 
            ram_economy=False, 
            args={'a': 1.5, 'b':2.},
            )

hps_bounds = np.array([[0.01,10.], #signal variance for the kernel
                       [0.01,10.], #length scale for the kernel
                       [0.00001,0.1],  #noise
                       [0.00001,1.]  #mean
                      ])

#the following is not needed, this is just to show how data is replced or appended
x_update = np.array([0.1,0.2,0.5]).reshape(3,1)
y_update = f1(x_update[:,0]) + (np.random.rand(len(x_update))-0.5) * 0.5
my_gpo.tell(x_update, y_update, append=True) ##append to the data
my_gpo.tell(x_data, y_data, append=False) ## back to normal overwriting the updated data
```

--------------------------------

### Custom Constrained Acquisition Function (Avoid Regions)

Source: https://github.com/lbl-camera/gpcam/blob/master/skills/acquisition-functions/SKILL.md

Implement a constrained acquisition function that penalizes exploration in specific forbidden zones. This example avoids a circular region.

```python
def constrained_variance(x, gpo):
    """Explore but avoid a circular forbidden zone."""
    var = gpo.posterior_covariance(x, variance_only=True)("v(x)")
    
    # Forbidden zone: circle at (5, 5) with radius 1
    center = np.array([5.0, 5.0])
    dist = np.linalg.norm(x - center, axis=1)
    penalty = np.where(dist < 1.0, -1e6, 0.0)
    
    return var + penalty
```

--------------------------------

### Initialize and Train GPOptimizer

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_MultiTaskTest.ipynb

Initializes the GPOptimizer with custom kernel and data, then performs global training using MCMC. Hyperparameter bounds must be defined before training.

```python
my_gp2 = fvGPOptimizer(x_data3,y_data3,
              init_hyperparameters=np.ones((n.number_of_hps+2))*0.1, kernel_function=deep_multi_task_kernel
             )
print("Global Training in progress")


bounds = np.zeros((n.number_of_hps+2,2))
bounds[0] = np.array([0.01,1.])
bounds[1] = np.array([0.1,1.])
bounds[2:] = np.array([-1,1])
my_gp2.train(hyperparameter_bounds=bounds,max_iter = 1000, method = "mcmc")
```

--------------------------------

### Acquiring Next Best Points with fvGPOptimizer

Source: https://github.com/lbl-camera/gpcam/blob/master/skills/multi-task-advanced/SKILL.md

Shows how to use the `ask` method of fvGPOptimizer to find the next best points for exploration. This can be done for a single task, multiple tasks, or in batches using different acquisition functions.

```python
# Ask for the next best input point across all tasks:
gpo.ask(parameter_bounds, x_out=np.array([0, 1, 2, 3, 4]), n=1)

# Ask for a batch of 4 points using a batch-aware acquisition:
gpo.ask(parameter_bounds, x_out=np.array([0, 1]), n=4,
        acquisition_function="relative information entropy set", vectorized=True)
```

--------------------------------

### Initialize and Train LogitGPOptimizer

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_LogAndLogit.ipynb

Initializes both a standard GPOptimizer and a LogitGPOptimizer with the simulated data. It then trains both models and captures any warnings issued by LogitGPOptimizer, specifically looking for messages related to clipping boundary observations.

```python
plain_gp2 = GPOptimizer(x_train2, y_train2)
plain_gp2.train(hyperparameter_bounds=hp_bounds_1d)

with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    logit_gp = LogitGPOptimizer(x_train2, y_train2)
    logit_gp.train(hyperparameter_bounds=hp_bounds_1d)

for w in caught:
    if "clipped" in str(w.message):
        print(f"  {w.message}")
```

--------------------------------

### Get GPOptimizer Data Keys

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_SingleTaskTest.ipynb

Retrieves the keys of the data dictionary stored within the GPOptimizer object. This is useful for understanding what information is available after training or data updates.

```python
my_gpo.get_data().keys()
```

--------------------------------

### Initialize GPOptimizer with a Cost Function

Source: https://github.com/lbl-camera/gpcam/blob/master/skills/cost-functions/SKILL.md

Instantiate GPOptimizer with data, a cost function, and optional arguments. The cost function is automatically applied during the 'ask()' method.

```python
gpo = GPOptimizer(
    x_data=x_data,
    y_data=y_data,
    cost_function=l2_cost,
    # args={\"speeds\": np.array([1.0, 2.0])}  # for anisotropic_cost
)
```

--------------------------------

### Optimizer Initialization with Learnable Noise Function

Source: https://github.com/lbl-camera/gpcam/blob/master/skills/noise-functions/SKILL.md

Shows how to initialize `GPOptimizer` with a custom learnable noise function using the `noise_function` argument.

```python
# Option B: learnable noise function
gpo = GPOptimizer(x_data, y_data, noise_function=learnable_noise)

```

--------------------------------

### Using the `ask` Method with Acquisition Functions

Source: https://github.com/lbl-camera/gpcam/blob/master/skills/acquisition-functions/SKILL.md

Demonstrates how to call the `ask` method, accepting either a built-in string name or a callable for the acquisition function. This is the primary way to query for new points during an experiment.

```python
# Built-in string or a callable are both accepted:
result = gpo.ask(
    input_set=parameter_bounds,
    acquisition_function=ucb,   # or "ucb", "expected improvement", ...
)
```

--------------------------------

### Quadratic Polynomial Prior Mean Function

Source: https://github.com/lbl-camera/gpcam/blob/master/skills/prior-mean-functions/SKILL.md

Implements a quadratic polynomial background as a prior mean for 1D input. The coefficients `a`, `b`, and `c` are read from hyperparameters starting at index `K`.

```python
def quadratic_mean(x, hps):
    """Quadratic background: a + b*x + c*x^2 (1D)."""
    K = 2  # adjust
    return hps[K] + hps[K+1] * x[:, 0] + hps[K+2] * x[:, 0]**2
```

--------------------------------

### Define a Simple Custom Kernel with Imported Functions

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_MultiTaskTest.ipynb

Defines a simple, stationary Matern kernel using functions imported from gpcam.kernels. This kernel is intended as a basic example and might not yield optimal results.

```python
#A simple kernel, that won't lead to good performance because it's stationary
from gpcam.kernels import *
def mkernel(x1,x2,hps):
    d = get_distance_matrix(x1,x2)
    return hps[0] * matern_kernel_diff1(d,hps[1])
```

--------------------------------

### Initialize and train GPOptimizer

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_Minimal.ipynb

Initializes the GPOptimizer with the generated data and trains the model. It then iteratively asks for new points, tells the model the function value, and retrains periodically.

```python
my_gp = GPOptimizer(x_data,y_data,)
my_gp.train()

train_at = [10,20,30] #optional
for i in range(100):
    new = my_gp.ask(np.array([[0.,1.]]))["x"]
    my_gp.tell(new, f1(new).reshape(len(new)))
    if i in train_at: my_gp.train()
```

--------------------------------

### Asynchronous Adam Training

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/1dSingleTaskAcqFuncTest.ipynb

Starts asynchronous training with the Adam optimizer. This is useful for distributed or local training tasks. A Dask client must be provided. The training loop updates hyperparameters iteratively and can be monitored.

```python
my_gpo.set_hyperparameters(np.array([1.,1.,1.,1.]))
print(my_gpo.hyperparameters)
opt_obj = my_gpo.train(hyperparameter_bounds=hps_bounds, dask_client=client, asynchronous=True, method='adam')
# The result won't change much (or at all) since this is such a simple optimization
for i in range(20):
    my_gpo.update_hyperparameters(opt_obj)
    print("iteration ", i, " : ",my_gpo.hyperparameters)
    time.sleep(0.1)
my_gpo.stop_training(opt_obj) ##this leaves the dask client alive, kill_client() will shut it down.
```

--------------------------------

### Optimize using the 'ask' method with string inputs

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_NonEuclideanInputSpaces.ipynb

Demonstrates the use of the 'ask' method for optimization with string inputs. It returns potential next points to query based on the current model, along with their associated objective function values.

```python
my_gp2.ask([('who'),('could'),("it"),("be")], n = 4, x_out=np.array([0,1,2,3]), vectorized=True)
```

--------------------------------

### Define a Stationary Matern Kernel Function

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_MultiTaskTest.ipynb

Defines a custom stationary Matern kernel function for Gaussian Processes. This kernel may not lead to optimal performance due to its stationary nature but serves as a basic example.

```python
#As imple kernel, that won't lead to good performance because its stationary
def mkernel(x1,x2,hps,obj):
    d = obj.get_distance_matrix(x1,x2)
    return hps[0] * obj.matern_kernel_diff1(d,hps[1])
```

--------------------------------

### Ask for Next Point

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_SingleTaskTest.ipynb

Queries the GPOptimizer for the next point to evaluate. Requires the optimizer to be initialized and trained.

```python
next_point = my_gpo.ask(np.array([[0.,1.]]))
print(next_point)
```

--------------------------------

### Evaluate Posterior on Original Scale

Source: https://github.com/lbl-camera/gpcam/blob/master/skills/transformed-optimizers-advanced/SKILL.md

Use evaluate_posterior to get the posterior distribution on the original scale. The median and credible interval transform exactly due to monotone links. Mean/std are exact for LogGPOptimizer and Monte-Carlo for LogitGPOptimizer.

```python
post = gpo.evaluate_posterior(x, level=0.95)
# returns: {"median", "mean", "std", "lower", "upper", "level"}
```

--------------------------------

### List Available Acquisition Functions

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/1dSingleTaskAcqFuncTest.ipynb

Displays a list of acquisition functions available for single-task Gaussian Process Optimization.

```python
#available acquisition function for the single-task case:
acquisition_functions = ["variance","relative information entropy","relative information entropy set",
                        "ucb","lcb","maximum","minimum","gradient","expected improvement",
                         "probability of improvement", "target probability", "total correlation"]
```

--------------------------------

### Linear Trend Prior Mean Function

Source: https://github.com/lbl-camera/gpcam/blob/master/skills/prior-mean-functions/SKILL.md

Defines a linear prior mean function of the form `m(x) = a + b0*x0 + b1*x1 + ...`. The intercept and slopes are read from the hyperparameter vector starting at index `K`.

```python
def linear_mean(x, hps):
    """
    Linear prior mean: m(x) = a + b0*x0 + b1*x1 + ...
    Uses hps[K], hps[K+1], ..., hps[K+D] where K is where mean hps start.
    """
    K = 3  # INDEX WHERE MEAN HYPERPARAMETERS START — adjust for your kernel
    D = x.shape[1]
    intercept = hps[K]
    slopes = hps[K+1:K+1+D]
    return intercept + x @ slopes
```

--------------------------------

### Initialize and Train GPOptimizer Model

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_MultiTaskTest.ipynb

Initializes the fvGPOptimizer with multi-task data and trains the model using different optimization methods: hybrid, MCMC, local, and global. Requires the Dask client for the 'hgdl' method.

```python
my_gp2 = fvGPOptimizer(x_data.reshape(len(x_data),1), np.column_stack([y_data1, y_data2]))
print("Hybrid Training in progress")
my_gp2.train(max_iter = 20, method = "hgdl", dask_client=client)

print("MCMC Training in progress")
my_gp2.train(max_iter = 20, method = "mcmc")

print("Local Training in progress")
my_gp2.train(max_iter = 20, method = "local")

print("Global Training in progress")
my_gp2.train(max_iter = 20, method="global")
```

--------------------------------

### Prepare string and numerical data

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_NonEuclideanInputSpaces.ipynb

Initializes sample string data (x_data) and corresponding numerical data (y_data) for training the Gaussian Process model. It prints the dimensions of the data for verification.

```python
x_data = ['frf','ferfe','ferf','febhn']
y_data = np.zeros((len(x_data),5))


y_data[:,0] = np.random.rand(len(x_data))
y_data[:,1] = np.random.rand(len(x_data))
y_data[:,2] = np.random.rand(len(x_data))
y_data[:,3] = np.random.rand(len(x_data))
y_data[:,4] = np.random.rand(len(x_data))

#it is good practice to check the format of the data
print(len(x_data))
print(y_data.shape)
```

--------------------------------

### Initialize and Train GPOptimizer with Missing Data

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_MultiTaskTest.ipynb

Initializes a GPOptimizer with prepared multi-task data and trains it for a specified number of iterations. This snippet shows the basic workflow for handling missing data points during training.

```python
my_gp2 = fvGPOptimizer(x_data.reshape(len(x_data),1), y_data, noise_variances=noise_variances)
print("Global Training in progress")
my_gp2.train(max_iter = 20)
```

--------------------------------

### Canonical Ask/Tell/Train Loop Usage

Source: https://github.com/lbl-camera/gpcam/blob/master/CLAUDE.md

Demonstrates the standard usage of GPOptimizer for Bayesian optimization, including training and iterative data acquisition.

```python
gp = GPOptimizer(x_data, y_data)
gp.train()                                   # inherited from fvgp.GP
for i in range(N):
    new = gp.ask(bounds, acquisition_function="variance")["x"]
    gp.tell(new, measure(new))               # appends + rank-n updates the GP
    if i in train_at: gp.train()
```

--------------------------------

### Initialize and train fvGPOptimizer

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_NonEuclideanInputSpaces.ipynb

Initializes the fvGPOptimizer with string data and the custom kernel function. It then trains the model using the specified hyperparameter bounds and maximum iterations.

```python
my_gp2 = fvGPOptimizer(x_data,y_data,init_hyperparameters=np.ones((2)),
              kernel_function=kernel
             )
print("MCMC Training in progress")
#use the next two lines if kernel `mkernel` is used
#if not a default deep kernel will be used that will set initi hyperparameters and bounds
#hps_bounds = np.array([[0.001,10000.],[1.,1000.]])
#my_gp2.train(hyperparameter_bounds = hps_bounds, max_iter = 2)

#use this next line if the default (deep) kernel is used (no bounds required)
my_gp2.train(hyperparameter_bounds=bounds, max_iter = 20)
```

--------------------------------

### Initialize and optimize vector-valued function

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_Optimization.ipynb

Initializes fvGPOptimizer for a vector-valued function and runs the optimization. The callback function visualizes the progress for multiple tasks simultaneously.

```python
from IPython.display import clear_output
my_gp2 = fvGPOptimizer()
def callb(x,y):  #OPTIONAL FUNCTION FOR PLOTTING
    task_ind0 = np.where(x[:,1]==0.)[0]
    task_ind1 = np.where(x[:,1]==1.)[0]
    clear_output(wait=True)
    y1 = y[task_ind0]
    y2 = y[task_ind1]
    x = x[task_ind0,0:1]
    plt.scatter(x[:,0], y1, color = 'grey', label = "data points t1")
    plt.scatter(x[:,0], y2, color = 'black', label = "data points t2")
    plt.scatter(x[-1,0], y1[-1], color = 'red', label = "last point")
    plt.scatter(x[-1,0], y2[-1], color = 'red', label = "last point")
    m = my_gp2.posterior_mean(x_pred1D, x_out = np.array([0,1]))["m(x)"]
    s = np.sqrt(my_gp2.posterior_covariance(x_pred1D, x_out = np.array([0,1]))["v(x)"])
    m1 = m[:,0]
    m2 = m[:,1]
    s1 = s[:,0]
    s2 = s[:,1]
    plt.plot(x_pred1D, m1, label = "post. mean 1")
    plt.plot(x_pred1D, m2, label = "post. mean 2")
    #plt.plot(x_pred1D, -(m-3.*s), label =  "acq func")
    plt.fill_between(x_pred1D.flatten(), m1-3.*s1, m1+3.*s1, label = "uncertainty1", color = "grey", alpha = .5)
    plt.fill_between(x_pred1D.flatten(), m2-3.*s2, m2+3.*s2, label = "uncertainty2", color = "grey", alpha = .5)
    plt.xlim([0,1])
    plt.ylim(-3,3)
    plt.legend(loc = 'lower left')
    plt.show()
    
result = my_gp2.optimize(func = f2, search_space =  np.array([[0,1]]), callback = callb)
```

--------------------------------

### GPOptimizer Training with Different Methods

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/1dSingleTaskAcqFuncTest.ipynb

Trains the GPOptimizer using various methods including MCMC, ADAM, global, local, and HGDL. Shows how to set hyperparameter bounds and iteration limits.

```python
st = time.time()
print("Standard Training (MCMC)")
hps = my_gpo.train(hyperparameter_bounds=hps_bounds, info = True, max_iter = 100)
print("Result=", hps, "after ", time.time() - st, " seconds")
print("")

print("ADAM")
hps = my_gpo.train(hyperparameter_bounds=hps_bounds, info = True, max_iter = 100, method="adam")
print("Result=", hps, "after ", time.time() - st, " seconds")
print("")

print("Global Training")
my_gpo.train(hyperparameter_bounds=hps_bounds, method='global', max_iter = 20)
print("Result=", hps, "after ", time.time() - st, " seconds")
print("")

print("Local Training")
my_gpo.train(hyperparameter_bounds=hps_bounds, method='local')
print("Result=", hps, "after ", time.time() - st, " seconds")
print("")

print("HGDL Training")
my_gpo.train(hyperparameter_bounds=hps_bounds, method='hgdl', max_iter=2, dask_client=client)
print("Result=", hps, "after ", time.time() - st, " seconds")
print("")
```

--------------------------------

### Configure and Train GPOptimizer with gp2Scale

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_gp2ScaleTest.ipynb

Initializes the GPOptimizer with `gp2Scale=True`, specifying batch size and Dask client. The model is then trained using specified hyperparameter bounds and maximum iterations.

```python
hps_n = 2

hps_bounds = np.array([[0.01,1.],      ##signal var of Wendland kernel
                       [0.001,0.02]])  ##length scale for Wendland kernel

init_hps = np.random.uniform(size = len(hps_bounds), low = hps_bounds[:,0], high = hps_bounds[:,1])

my_gp2S = GPOptimizer(x_data,y_data,init_hyperparameters=np.array([0.73118673, 0.0013813191]), 
            gp2Scale = True, gp2Scale_batch_size= 500, dask_client = client
            )

my_gp2S.train(hyperparameter_bounds=hps_bounds, max_iter = 25, info = True) ##usually you would do >200 updates
```

--------------------------------

### Initialize and optimize scalar-valued function

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_Optimization.ipynb

Initializes GPOptimizer for a scalar-valued function and runs the optimization process. Includes an optional callback function for real-time plotting of optimization progress.

```python
from IPython.display import clear_output
my_gp1 = GPOptimizer()
def callb(x,y): #OPTIONAL FUNCTION FOR PLOTTING
    clear_output(wait=True)
    plt.scatter(x, y, color = 'black', label = "data points")
    plt.scatter(x[-1], y[-1], color = 'red', label = "last point")
    m = my_gp1.posterior_mean(x_pred1D)['m(x)']
    s = np.sqrt(my_gp1.posterior_covariance(x_pred1D)['v(x)'])
    plt.plot(x_pred1D, m, label = "post. mean")
    plt.plot(x_pred1D, -(m-3.*s), label =  "acq func")
    plt.fill_between(x_pred1D.flatten(), m-3.*s, m+3.*s, label = "uncertainty", color = "grey", alpha = .5)
    plt.xlim([0,1])
    plt.ylim(-3,3)
    plt.legend(loc = 'lower left')
    plt.show()
    
result = my_gp1.optimize(func = f1, search_space =  np.array([[0,1]]), callback=callb)
```

--------------------------------

### Suggest next measurement point

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_NonEuclideanInputSpaces.ipynb

Uses the 'ask' method of the GPOptimizer to suggest the next data points to measure based on the current model. This is useful for active learning or experimental design.

```python
##which one should I measure next?
my_gp.ask([('who'),('could'),("it"),("be")], n = 4)
```

--------------------------------

### Multi-Task Optimization with Various Acquisition Functions

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_MultiTaskTest.ipynb

Demonstrates calling the `ask` method with different acquisition functions for both n=1 and n>1 scenarios. Note that 'total correlation' and 'relative information entropy' may automatically set `vectorized=False`.

```python
print("n=1")
my_gp2.ask(np.array([[0,1],[0,1]]), n = 1, max_iter=2,pop_size=2, info = True,
           acquisition_function = 'relative information entropy set', x_out = np.array([0.,1.]), vectorized = True)
my_gp2.ask(np.array([[0,1],[0,1]]), n = 1, max_iter=2,pop_size=2, info = True,
           acquisition_function = 'relative information entropy', x_out = np.array([0.,1.]), vectorized = True)
my_gp2.ask(np.array([[0,1],[0,1]]), n = 1, max_iter=2,pop_size=2,info = True,
           acquisition_function = 'variance', x_out = np.array([0.,1.]), vectorized = True)
my_gp2.ask(np.array([[0,1],[0,1]]), n = 1, max_iter=2,pop_size=2,info = True,
           acquisition_function = 'total correlation', x_out = np.array([0.,1.]), vectorized = True)

print("n>1")
my_gp2.ask(np.array([[0,1],[0,1]]), n = 4, max_iter=2,pop_size=2,info = True,
           acquisition_function = 'relative information entropy set', x_out = np.array([0.,1.]), vectorized = True)
my_gp2.ask(np.array([[0,1],[0,1]]), n = 5, max_iter=2,pop_size=2,info = True,
           acquisition_function = 'relative information entropy', x_out = np.array([0.]), vectorized = True)
my_gp2.ask(np.array([[0,1],[0,1]]), n = 3, max_iter=2,pop_size=2,info = True,
           acquisition_function = 'variance', x_out = np.array([0.,1.]), vectorized = True)
my_gp2.ask(np.array([[0,1],[0,1]]), n = 2, max_iter=2,pop_size=2,info = True,
           acquisition_function = 'total correlation', vectorized = True)

```

--------------------------------

### Train GPOptimizer with Custom Kernel and Hyperparameter Bounds

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_MultiTaskTest.ipynb

Initializes a GPOptimizer with specific data, a custom kernel function, and initial hyperparameters. It then trains the model with defined hyperparameter bounds for a limited number of iterations.

```python
my_gp2 = fvGPOptimizer(x_data3,y_data3,
              init_hyperparameters=np.ones((2)), kernel_function=mkernel
             )
print("Global Training in progress")


bounds = np.array([[0.01,1.],[0.01,1.]])
my_gp2.train(hyperparameter_bounds=bounds,max_iter = 2)
```

--------------------------------

### Import Libraries and Initialize Dask Client

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/GPOptimizer_MultiTaskTest.ipynb

Imports essential libraries for GPOptimizer, plotting, and Dask distributed computing. Initializes a Dask client for parallel processing.

```python
import numpy as np
import matplotlib.pyplot as plt
from gpcam import fvGPOptimizer
import plotly.graph_objects as go
import plotly.io as pio
pio.renderers.default = "png"
from itertools import product
from distributed import Client
client = Client()

%load_ext autoreload
%autoreload 2
```

--------------------------------

### Parallel and Vectorized ask() Evaluations

Source: https://github.com/lbl-camera/gpcam/blob/master/examples/1dSingleTaskAcqFuncTest.ipynb

This snippet demonstrates different strategies for evaluating the acquisition function on a list of candidates using the ask() function. It compares sequential asking, parallel asking on Dask workers with sequential evaluation on each worker, parallel asking on Dask workers with vectorized evaluation on each worker, and purely vectorized asking. The output shows that all methods yield the same suggestions.

```python
#ask sequentially
print("suggestions=", my_gpo.ask(candidate_list, n = 30, acquisition_function="variance", vectorized=False)[ "x"] [0])
#ask in parallel on DASK workers, but sequentially on each worker:
print("suggestions=", my_gpo.ask(candidate_list, n = 30, acquisition_function="variance", vectorized=False, batch_size = 10, dask_client=client)[ "x"] [0])
#ask in parallel on DASK workers, and vectorized (if possible) on each worker:
print("suggestions=", my_gpo.ask(candidate_list, n = 30, acquisition_function="variance", vectorized=True, batch_size = 10, dask_client=client)[ "x"] [0])
#ask vectorized (if possible):
print("suggestions=", my_gpo.ask(candidate_list, n = 30, acquisition_function="variance", vectorized=True)[ "x"] [0])
print("They should be the same!")
```