<jupyter_start><jupyter_text>Ball Drop with Varying RadiiConsider dropping balls of various radii $R$ from a tower at heights 5, 10, 15, and 20 meters. Each dropped ball experiment produces vector output in the form of a height-time curve, i.e., a curve showing time as a function of the current height $t=t(h)$.Rather than performing these experiments, we generate synthetic experimental data from:\begin{equation}\frac{d^2h}{dt^2} = g - \frac{C}{R}\bigg(\frac{dh}{dt}\bigg)^2\end{equation}where $g$ is the acceleration due to gravity, $C$ is the coefficient of drag, and R is the radius of the ball. We will refer to data generated from equation (1) as system observations. To illustrate the modeling of discrepency between simulator and reality, we will generate simulated data from:\begin{equation}\frac{d^2h}{dt^2} = g - \frac{C}{R}\frac{dh}{dt}\end{equation}Note the simulated data will be systematically biased from the experimental data as it is generated with a linear velocity term $(dh/dt)$ rather than a squared velocity term. We will explicitly model this discrepecy.These differential equations can be solved for height $h$ as a function of time $t$, ball radius $R$, and drag $C$. Observations will be the inverse: time recorded as a function of height, radius, and drag which are calculated by an optimizer from the scipy library in Python.Acceleration due to gravity $g$ is assumed to be known, but the coefficient of drag $C$ is not, and is the target of inference. That is, the goal is to find the posterior distribution of $C$ that corresponds to the Sepia statistical model fit to all available data.The first step is to compute system observations from (1) and simulated data from (2). As is standard practice in the study of computer codes a (space-filling) Latin hypercube (LH) design is used to efficiently collect simulator data over the domain of model parameters, in this case $(R,C)$.Suppose we have data from $n=3$ experiments, one each for balls of radius $R={0.1,0.2,0.4}$ meters. Each experiment produces a curve of drop times made of four height-time pairs. The experimental heights are $h_e={5,10,15,20}$ meters. We assume a zero mean, normally distributed observation error $\epsilon \sim N(0,et)$ (IS THIS TRUE?). We fix $C=0.1*4\pi/3$, the coefficient of drag for a sphere, and $g=9.8 \; m/s^2$.<jupyter_code>et = 0.01                      # observation error
data_dict = gen_data(et)

# field data
R = data_dict['R']             # radii of balls   .1,.2,.4 (m)
h_field = data_dict['h_field'] # observed heights 5,10,15,20 (m)
y_field = data_dict['y_field'] # observed times<jupyter_output>generating data with C =  0.023873241463784306<jupyter_text>For simulation design points, we use a scaled latin hypercube to select $m=25\;(R,C)$ pairs. We then scale the design so $R\in[0.05,0.45]$ to better match experimental data and $C\in[0.05,0.25]$ is chosen with some hindsight. The simulator will be run at these 25 design points and will return a time curve at 16 evenly spaced heights between 1.5 and 25 meters.<jupyter_code># sim data
sim_design = data_dict['sim_design']

# to generate a different design
#import pyDOE
#sim_design = pyDOE.lhs(2,data_dict['m'])*np.array([.4,.2])+.05
#data_dict['sim_design'] = sim_design

R_sim = sim_design[:,0]
C_sim = sim_design[:,1]
h_sim = data_dict['h_sim']
y_sim = data_dict['y_sim']<jupyter_output><empty_output><jupyter_text>See invertH.py for details on solving and inverting equations (1) and (2) to obtain experimental and simulated data. Below we show, for each Radius R, the physical reality, defined by equation (1), the noisy observations, and all 25 simulation realizations. For each R, the 3 simulations with design points nearest to R are shown as red, green, and blue lines.<jupyter_code># see GenDataBallDrop1.py for ploting details
plot_data(data_dict)<jupyter_output><empty_output><jupyter_text>Preparing the data for SepiaTo use Sepia, we must package our data into a SepiaData object. In this example, the known inputs to the simulator are simply the vector of radii R_sim. We pass this into x_sim as a column vector. Similarly, C_sim is passed into t_sim as a column vector and is the parameter to be calibrated. We also pass in y_sim, the simulated time-height curves, and h_sim the heights associated with those times in y_sim. For the observed data, x_obs get the experimental radii R, and y_obs gets the experimental time-height curves generated from (1). Finally we pass in the heights at which the experimental time observations were taken, y_ind_obs = h_field.<jupyter_code>data = SepiaData(x_sim = np.reshape(R_sim,(len(R_sim),1)),
                 t_sim = np.reshape(C_sim,(len(C_sim),1)), 
                 y_sim = y_sim, y_ind_sim = h_sim,
                 x_obs = np.reshape(R,(len(R),1)), y_obs = y_field, y_ind_obs=h_field)<jupyter_output><empty_output><jupyter_text>Transforming x, t, and ySepia required that the inputs $x,t$ lie in the interval $[0,1]^{p+q}$, and the responses $y_{sim},y_{obs}$ be $N(0,1)$.<jupyter_code>data.transform_xt()
data.standardize_y()<jupyter_output><empty_output><jupyter_text>Generate K basisSepia models multivariate observations and responses using a linear basis. These *principal components*, or scaled eigenvectors, are computed by the singular value decomposition.<jupyter_code>data.create_K_basis(2)
fig_sim, fig_obs = SepiaPlot.plot_K_basis(data)<jupyter_output><empty_output><jupyter_text>We can visually check that our principal component weights are reasonable. In an ideal setting, weights would follow a standard normal distribution.<jupyter_code>fig_uw = SepiaPlot.plot_K_weights(data)
fig_uw_pairs = SepiaPlot.plot_u_w_pairs(data)<jupyter_output><empty_output><jupyter_text>Model initialization and MCMCNow that our data is appropiately transformed, and out bases set, the model can be initialized and the posterior distributions of parameters sampled via Markov chain Monte Carlo (MCMC). Before performing MCMC, we call tune_step_sizes() to optimize proposal widths for each parameter.<jupyter_code>model = SepiaModel(data)
nsamp = 1000
if os.path.isfile(model_file+'.pkl') and load_model:
    model.restore_model_info(cachefile_name)
else:
    model.tune_step_sizes(50, 20)
    model.do_mcmc(1000)
    if save_model: model.save_model_info(cachefile_name)
model = SepiaModel(data)
nsamp = 1000
if os.path.isfile(model_file+'.pkl') and load_model:
    model.restore_model_info(cachefile_name)
else:
    model.tune_step_sizes(50, 20)
    model.do_mcmc(1000)
    if save_model: model.save_model_info(cachefile_name)
samples_dict = model.get_samples()
theta_names = ['C']
theta_pairs = SepiaPlot.theta_pairs(samples_dict,theta_names)
mcmc_trace = SepiaPlot.mcmc_trace(samples_dict,theta_names)
p_stats = SepiaPlot.param_stats(samples_dict,theta_names=theta_names,q1=.05,q2=.95,digits=4)
print(p_stats)
acf = SepiaPlot.plot_acf(model,50)<jupyter_output><empty_output><jupyter_text>Assessing emulator adequacyIt is important to examine diagnostics to understand the quality of the emulator in reproducing the simulator output at new, untried settings. The following are usful ways to asses the quality of the emulator fit.  Principal components  Parameters in the Gaussian Process FitUsing the MCMC draws of the spacial dependence parameters $\beta$, $\rho=exp\{-\beta/4\}$ can be calculated to show the values on the more interpretable $[0,1]$ scale. The values of $\rho$ give us information about the dependence of the simulation output on each input parameter $x$ and $\theta$. We show boxplots of the posterior draws for $\rho$ for each $x$ and $\theta$. and for each principal component.  When $\rho$ is near 1, it suggests that particular component of the simulator output in linear in that dimention. As $\rho$ goes smaller, nonlinear activity is associated with that input. In other words larger values of $\rho$ indicate more smoothness. Additionally, smaller values of $\rho$ indicate increasing flexibility in the model response. Very small values could indicate overfitting. Thus, if any of the boxplots show values very near to zero, the model is suspect. The boxplots below are for correlation parameters associated with the K-basis. Discrepency basis correlation parameters could be similarly assesed.<jupyter_code>rho_box = SepiaPlot.rho_box_plots(model,labels=[r'$R_{ball}$','C'])<jupyter_output><empty_output><jupyter_text>---- Predictions of model output from the full modelFirst, predictions showing just the emulator response, usingSepiaEmulatorPredict<jupyter_code>model.verbose=True
x_pred=data.obs_data.x
n_pred=x_pred.shape[0]
pred_samples=model.get_samples(nburn=int(.1*nsamp),effectivesamples=True)
pred=SepiaEmulatorPrediction(x_pred=x_pred, samples=pred_samples, model=model)#, t_pred=t_pred)
predw=pred.get_w()
plt.figure()
plt.plot(np.mean(predw,0).T,'--o')
plt.xlabel('PCA weights w')
plt.ylabel('w value for each predicted x \n mean over samples')
plt.show()
predystd=pred.get_y(std=True)
plt.figure()
plt.plot(model.data.sim_data.y_ind,np.mean(predystd,0).T)
plt.xlabel('native data space y')
plt.ylabel('standardized predicted y for each predicted x \n mean over samples')
plt.show()
predy=pred.get_y()
plt.figure()
plt.plot(model.data.sim_data.y_ind,np.mean(predy,0).T)
plt.xlabel('native data space y ')
plt.ylabel('native (original) predicted value for each x_pred \n mean over samples')
plt.show()
plt.figure()
plt.plot(model.data.sim_data.y_ind,np.mean(predy[:,0,:],0).T,label='R=.1')
plt.plot(model.data.sim_data.y_ind,np.mean(predy[:,1,:],0).T,label='R=.2')
plt.plot(model.data.sim_data.y_ind,np.mean(predy[:,2,:],0).T,label='R=.4')
plt.plot(model.data.sim_data.y_ind,np.quantile(predy[:,0,:],[0.05, 0.95],axis=0).T,'k:')
plt.plot(model.data.sim_data.y_ind,np.quantile(predy[:,1,:],[0.05, 0.95],axis=0).T,'k:')
plt.plot(model.data.sim_data.y_ind,np.quantile(predy[:,2,:],[0.05, 0.95],axis=0).T,'k:')
plt.xlabel('native data space y ')
plt.ylabel('native (original) mean and (0.05,0.95) quantiles')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Summary plot of predicted uncertainty for an observation.The 5th prediction done corresponds to the observation in the model, with xpred=0.5.<jupyter_code>predysim=pred.get_y()
colors = ('r','g','b')
for R_idx in range(3):
    plt.figure()
    plt.plot(model.data.obs_data.y_ind,model.data.obs_data.y[R_idx,:],marker='o',color=colors[R_idx])
    plt.plot(model.data.sim_data.y_ind,np.mean(predysim[:,R_idx,:],0).T,'--',c='lightgreen')
    plt.plot(model.data.sim_data.y_ind,np.quantile(predysim[:,R_idx,:],.05,0),'k:')
    plt.plot(model.data.sim_data.y_ind,np.quantile(predysim[:,R_idx,:],.95,0),'k:')
    plt.xlabel('native scale (obs grid)')
    plt.ylabel('native scale response')
    plt.title('Observation and model prediction')
    plt.legend(['obs R=.1','pred mean','pred q5/q95'])
    plt.ylim(0,4.5)
    plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11230191444536045, 'fraction_numerical': 0.019510854630392965, 'mean_word_length': 5.573148705599037, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 44, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '620857', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3676, 'n_tokens_neox': 3471, 'n_words': 1153, 'path': 'examples/Ball_Drop/ball_drop_1_noD.ipynb', 'repo_name': 'lanl/SEPIA', 'size': '10917'}