<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Figure 1
Start by loading some boiler plate
Step1:  And some more specialized dependencies
Step2:  Configuration for this figure.
Step3:  Open a chest located on a remote globus endpoint and load a remote json configuration file.
Step4:  We want to plot the spike depth, which is the 'H' field in the chest.
Chests can prefetch lists of keys more quickly than individual ones, so we'll prefetch the keys we want.
Step5:  Plot the bubble height, the 'H' keys, vs. time.
Use a spline to compute the derivative.
Step6:  What if we want to dig deeper?
The spike depth, 'H', is the smallest depth for which the span-wise average of the density is 99% of its asymptotic value.
What do the density profiles themselves look like?
The density profiles are stored as 't_proj_z'
Step7:  Now let's plot it.
Step8:  That doesn't look great.  Let's zoom in.  The data is cached, so this is fast!
Step9:  Now where's that 99% threshold?
Step10:  Let's re-do the analysis for a 95% threshold to make sure nothing changes.
Step11:  That didn't change much.  What about 90%?
Step12:  That introduces an offset, but the spike velocity (slope) looks the same for late time.
I've convinced myself that the results in the original figure are robust to changes in parameters.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib
matplotlib.rcParams['figure.figsize'] = (10.0, 8.0)
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d, InterpolatedUnivariateSpline
from scipy.optimize import bisect
import json
from functools import partial
class Foo: pass
Explanation: Figure 1
Start by loading some boiler plate: matplotlib, numpy, scipy, json, functools, and a convenience class.
End of explanation
from chest import Chest
from slict import CachedSlict
from glopen import glopen, glopen_many
Explanation: And some more specialized dependencies:
 1. Slict provides a convenient slice-able dictionary interface
 2. Chest is an out-of-core dictionary that we'll hook directly to a globus remote using...
 3. glopen is an open-like context manager for remote globus files
End of explanation
config = Foo()
config.name     = "HighAspect/HA_conductivity_1.0E-5/HA_conductivity_1.0E-5"
#config.arch_end = "maxhutch#alpha-admin/~/pub/"
config.arch_end = "alcf#dtn_mira/projects/alpha-nek/experiments"
Explanation: Configuration for this figure.
End of explanation
c = Chest(path      = "{:s}-results".format(config.name),
          open      = partial(glopen,      endpoint=config.arch_end),
          open_many = partial(glopen_many, endpoint=config.arch_end),
          available_memory=1e12)
sc = CachedSlict(c)
with glopen(
            "{:s}.json".format(config.name), mode='r',
            endpoint = config.arch_end,
            ) as f:
    params = json.load(f)
Explanation: Open a chest located on a remote globus endpoint and load a remote json configuration file.
End of explanation
c.prefetch(sc[:,'H'].full_keys())
c.prefetch(sc[:,'h'].full_keys())
c.prefetch(sc[:,'flux_proj_z'].full_keys())
c.prefetch(sc[:,'frame'].full_keys())
c.prefetch(sc[:,'t_yz'].full_keys())
c.prefetch(sc[:,'w_yz'].full_keys())
Explanation: We want to plot the spike depth, which is the 'H' field in the chest.
Chests can prefetch lists of keys more quickly than individual ones, so we'll prefetch the keys we want.
End of explanation
spl = InterpolatedUnivariateSpline(sc[:,'H'].keys(), sc[:,'h'].values(), k=3)
dHdt = spl.derivative()
Ts = np.linspace(sc[:,'H'].keys()[0], sc[:,'H'].keys()[-1], 1000)
V = -dHdt(Ts)
theory0 = np.sqrt(params["atwood"] * params["g"] * 2 * params["extent_mesh"][0] / np.pi)
                
theory1 = np.sqrt(
                  params["atwood"] * params["g"] * params["extent_mesh"][0] / np.pi 
                + (2.*np.pi*params["viscosity"] / params["extent_mesh"][0])**2
                ) - (2.*np.pi*params["viscosity"] / params["extent_mesh"][0]) 
theory2 = params["atwood"] * params["g"] * params["extent_mesh"][0]**2 / (
          113.816*params["viscosity"])
theory3 = 0.008440 * params["atwood"] * params["g"] * params["extent_mesh"][0]**2 / params["viscosity"]
flux_midplane = np.array([sc[t,'flux_proj_z'][1892] for t in sc[:,'flux_proj_z'].keys()])
flux_total = np.array([np.sum(sc[t,'flux_proj_z']) for t in sc[:,'flux_proj_z'].keys()])
fig, axs = plt.subplots(2,1, sharex=True)
axs[0].plot(Ts, spl(Ts)/128);
axs[0].plot(sc[:,'H'].keys(), -np.array(sc[:,'H'].values()), 'yo');
axs[0].set_ylabel('Depth (m)')
#axs[1].plot(sc[:,'frame'].values(), -V/128);
#axs[1].plot([Ts[0],Ts[-1]], [theory0, theory0], 'k--');
#axs[1].plot([Ts[0],Ts[-1]], [theory1, theory1], 'k-.');
#axs[1].plot([Ts[0],Ts[-1]], [theory2, theory2], 'k--');
#axs[1].plot([Ts[0],Ts[-1]], [theory3, theory3], 'k-.');
#axs[1].plot(sc[:,'frame'].values(), 16.*flux_midplane)
axs[1].plot(sc[:,'frame'].values(), sc[:,'Dissipated'].values())
#axs[1].plot(sc[:,'frame'].values(), flux_total)
axs[1].set_ylabel('Velocity (m/s)')
axs[1].set_xlabel('Time (s)');
#plt.xlim([20,40])
plt.savefig('Figure1.eps')
diss = np.array([np.sum(sc[T,'d_yz']) for T in sc[:,'d_yz'].keys()]);
plt.plot(diss)
flux = []
vflux = []
for T in sc[:,'w_yz'].keys():
    tmp = sc[T,'w_yz'] * sc[T,'t_yz']
    tmp[sc[T,'t_yz'] < 0] = 0
    tmp[tmp < 0] = 0.
    flux.append(np.sum(tmp) / sc[T,'h'])
    tmp = sc[T,'w_yz'].copy()
    tmp[tmp < 0] = 0.
    vflux.append(np.sum(tmp) / sc[T,'h'])
plt.plot(flux)
plt.plot(vflux)
Explanation: Plot the bubble height, the 'H' keys, vs. time.
Use a spline to compute the derivative.
End of explanation
c.prefetch(sc[:,'t_proj_z'].full_keys())
Explanation: What if we want to dig deeper?
The spike depth, 'H', is the smallest depth for which the span-wise average of the density is 99% of its asymptotic value.
What do the density profiles themselves look like?
The density profiles are stored as 't_proj_z': the temperature projected onto the z axis.
    Let's prefech it.
End of explanation
Ts = sc[:,'t_proj_z'].keys()
Zs = np.linspace(params["root_mesh"][2], params["extent_mesh"][2], sc[Ts[0], 't_proj_z'].shape[0])
for t in Ts:
    plt.plot(Zs,sc[t,'t_proj_z'], 'k-')
plt.xlim(xmin=params["root_mesh"][2],xmax=params["extent_mesh"][2]);
plt.ylim(ymin=-1.1, ymax=1.1)
plt.xlabel("Height (m)");
plt.ylabel(r"$\Delta \rho / \bar{\rho}$");
plt.title("Density profiles");
Explanation: Now let's plot it.
End of explanation
Ts = sc[:,'t_proj_z'].keys()
Zs = np.linspace(params["root_mesh"][2], params["extent_mesh"][2], sc[Ts[0], 't_proj_z'].shape[0])
for t in Ts:
    plt.plot(Zs,sc[t,'t_proj_z'], 'k-')
plt.xlim(xmin=-1.5,xmax=1.5);
plt.ylim(ymin=-1.1, ymax=1.1);
plt.xlabel("Height (m)");
plt.ylabel(r"$\Delta \rho / \bar{\rho}$");
plt.title("Density profiles");
Explanation: That doesn't look great.  Let's zoom in.  The data is cached, so this is fast!
End of explanation
Ts = sc[:,'t_proj_z'].keys()
Zs = np.linspace(params["root_mesh"][2], params["extent_mesh"][2], sc[Ts[0], 't_proj_z'].shape[0])
for t in Ts:
    plt.plot(Zs,sc[t,'t_proj_z'], 'k-')
plt.plot([-1.5, 0.], np.array([.98, .98]));
plt.xlim(xmin=-1.5,xmax=0);
plt.xlabel("Height (m)");
plt.ylabel(r"$\Delta \rho / \bar{\rho}$");
plt.title("Density profiles");   
plt.ylim(ymin=0.9,ymax=1.1);
Explanation: Now where's that 99% threshold?
End of explanation
Ts = sc[:,'t_proj_z'].keys()
Zs = np.linspace(params["root_mesh"][2], params["extent_mesh"][2], sc[Ts[0], 't_proj_z'].shape[0])
my_H = []
for t in Ts:
    f = interp1d(Zs, sc[t, 't_proj_z']-(.90))
    my_H.append(bisect(f,params["root_mesh"][2], 0))
my_spl = InterpolatedUnivariateSpline(sc[:,'H'].keys(), my_H, k=4)
my_dHdt = my_spl.derivative()
Ts = np.linspace(sc[:,'H'].keys()[0], sc[:,'H'].keys()[-1], 1000)
my_V = -my_dHdt(Ts) 
theory = np.sqrt(
                  params["atwood"] * params["g"] * params["extent_mesh"][0] / np.pi 
                + (2.*np.pi*params["viscosity"] / params["extent_mesh"][0])**2
                ) - (2.*np.pi*params["viscosity"] / params["extent_mesh"][0]) 
fig, axs = plt.subplots(2,1, sharex=True)
axs[0].plot(Ts, -spl(Ts));
axs[0].plot(Ts, -my_spl(Ts), 'k--');
axs[0].plot(sc[:,'H'].keys(), -np.array(my_H), 'yo');
axs[0].set_ylabel('Depth (m)')
axs[1].plot(Ts, V);
axs[1].plot(Ts, my_V, 'k');
axs[1].plot([Ts[0],Ts[-1]], [theory, theory], 'k--');
axs[1].set_ylabel('Velocity (m/s)')
axs[1].set_xlabel('Time (s)');
Explanation: Let's re-do the analysis for a 95% threshold to make sure nothing changes.
End of explanation
Ts = sc[:,'t_proj_z'].keys()
Zs = np.linspace(params["root_mesh"][2], params["extent_mesh"][2], sc[Ts[0], 't_proj_z'].shape[0])
my_H = []
for t in Ts:
    f = interp1d(Zs, sc[t, 't_proj_z']-(.80))
    my_H.append(bisect(f,params["root_mesh"][2], 0))
my_spl = InterpolatedUnivariateSpline(sc[:,'H'].keys(), my_H, k=4)
my_dHdt = my_spl.derivative()
Ts = np.linspace(sc[:,'H'].keys()[0], sc[:,'H'].keys()[-1], 1000)
my_V = -my_dHdt(Ts) 
theory = np.sqrt(
                  params["atwood"] * params["g"] * params["extent_mesh"][0] / np.pi 
                + (2.*np.pi*params["viscosity"] / params["extent_mesh"][0])**2
                ) - (2.*np.pi*params["viscosity"] / params["extent_mesh"][0]) 
fig, axs = plt.subplots(2,1, sharex=True)
axs[0].plot(Ts, -spl(Ts));
axs[0].plot(Ts, -my_spl(Ts), 'k--');
axs[0].plot(sc[:,'H'].keys(), -np.array(my_H), 'yo');
axs[0].set_ylabel('Depth (m)')
axs[1].plot(Ts, V);
axs[1].plot(Ts, my_V, 'k');
axs[1].plot([Ts[0],Ts[-1]], [theory, theory], 'k--');
axs[1].set_ylabel('Velocity (m/s)')
axs[1].set_xlabel('Time (s)');
Explanation: That didn't change much.  What about 90%?
End of explanation
%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py 
%load_ext version_information 
%version_information numpy, matplotlib, slict, chest, glopen, globussh
Explanation: That introduces an offset, but the spike velocity (slope) looks the same for late time.
I've convinced myself that the results in the original figure are robust to changes in parameters.
End of explanation
<END_TASK>