#!/usr/bin/env python
# coding: utf-8

# ###Note: this tutorial needs updating and has not been recently tested for basic functionality: use at own risk

# In[1]:


import os
import sys
import time
import json

import pygsti

#The following two lines should be changed to reflect your own paths to GST/src and GSTData/data/Ion_Trap_SNL/2015_03_30/GST_BB1_XX_512_condensed.txt 
#sys.path.append('/Users/kmrudin/GST/src/')
#dataFile = '/Users/kmrudin/GSTData/data/Ion_Trap_SNL/2015_03_30/GST_BB1_XX_512_condensed.txt'
#import Bootstrap
get_ipython().run_line_magic('pylab', 'inline')


# In[2]:


#Load example quantities from files
gs_target = pygsti.io.load_gateset("tutorial_files/Example_Gateset.txt")
gs_mc2gst = pygsti.io.load_gateset("tutorial_files/Example_MC2GST_Gateset.txt")

ds = pygsti.io.load_dataset("tutorial_files/Example_Dataset.txt", cache=True)

fiducials = pygsti.io.load_gatestring_list("tutorial_files/Example_FiducialList.txt")
germs = pygsti.io.load_gatestring_list("tutorial_files/Example_GermsList.txt")
maxLengths = json.load(open("tutorial_files/Example_maxLengths.json","r"))

specs = pygsti.construction.build_spam_specs(fiducials)


# ##Parametric Bootstrapping
# Here we do parametric bootstrapping, as indicated by the 'parametric' argument below.
# The output is eventually stored in the "mean" and "std" GateSets, which hold the mean and standard deviation values of the set of bootstrapped gatesets (after gauge optimization).  It is this latter "standard deviation Gateset"
# which holds the collection of error bars.  Note: due to print setting issues, the outputs that are printed here will not necessarily reflect the true accuracy of the estimates made.
# 

# In[3]:


#The number of simulated datasets & gatesets made for bootstrapping purposes.  
# For good statistics, should probably be greater than 10.
numGatesets=10

param_boot_gatesets = pygsti.drivers.make_bootstrap_gatesets(
                        numGatesets, ds, 'parametric', fiducials, fiducials, germs, maxLengths,
                        inputGateSet=gs_mc2gst, startSeed=0, constrainToTP=True, returnData=False,
                        verbosity=2)


# In[4]:


gauge_opt_pboot_gatesets = pygsti.drivers.gauge_optimize_gs_list(param_boot_gatesets, gs_mc2gst,
                                                                 constrainToTP=True, plot=True)


# In[6]:


pboot_mean = pygsti.drivers.to_mean_gateset(gauge_opt_pboot_gatesets, gs_mc2gst)
pboot_std  = pygsti.drivers.to_std_gateset(gauge_opt_pboot_gatesets, gs_mc2gst)

#Summary of the error bars
print "Parametric bootstrapped error bars, with",numGatesets,"resamples\n"
print "Error in rho vec:" 
print pboot_std.rhoVecs[0]
print
print "Error in E vec:"
print pboot_std.EVecs[0]
print
print "Error in Gi:"
print pboot_std['Gi']
print
print "Error in Gx:"
print pboot_std['Gx']
print
print "Error in Gy:"
print pboot_std['Gy']


# ##Non-parametric Bootstrapping
# Here we do non-parametric bootstrapping, as indicated by the 'nonparametric' argument below.
# The output is again eventually stored in the "mean" and "std" GateSets, which hold the mean and standard deviation values of the set of bootstrapped gatesets (after gauge optimization).  It is this latter "standard deviation Gateset"
# which holds the collection of error bars.  Note: due to print setting issues, the outputs that are printed here will not necessarily reflect the true accuracy of the estimates made.
# 
# (Technical note: ddof = 1 is by default used when computing the standard deviation -- see numpy.std -- meaning that we are computing a standard deviation of the sample, not of the population.)

# In[7]:


#The number of simulated datasets & gatesets made for bootstrapping purposes.  
# For good statistics, should probably be greater than 10.
numGatesets=10

nonparam_boot_gatesets = pygsti.drivers.make_bootstrap_gatesets(
                          numGatesets, ds, 'nonparametric', fiducials, fiducials, germs, maxLengths,
                          targetGateSet=gs_mc2gst, startSeed=0, constrainToTP=True, returnData=False,
                          verbosity=2)


# In[8]:


gauge_opt_npboot_gatesets = pygsti.drivers.gauge_optimize_gs_list(nonparam_boot_gatesets, gs_mc2gst,
                                                                 constrainToTP=True, plot=True)


# In[9]:


npboot_mean = pygsti.drivers.to_mean_gateset(gauge_opt_npboot_gatesets, gs_mc2gst)
npboot_std  = pygsti.drivers.to_std_gateset(gauge_opt_npboot_gatesets, gs_mc2gst)

#Summary of the error bars
print "Non-parametric bootstrapped error bars, with",numGatesets,"resamples\n"
print "Error in rho vec:" 
print npboot_std.rhoVecs[0]
print
print "Error in E vec:"
print npboot_std.EVecs[0]
print
print "Error in Gi:"
print npboot_std['Gi']
print
print "Error in Gx:"
print npboot_std['Gx']
print
print "Error in Gy:"
print npboot_std['Gy']


# In[15]:


loglog(npboot_std.to_vector(),pboot_std.to_vector(),'.')
loglog(np.logspace(-4,-2,10),np.logspace(-4,-2,10),'--')
xlabel('Non-parametric')
ylabel('Parametric')
xlim((1e-4,1e-2)); ylim((1e-4,1e-2))
title('Scatter plot comparing param vs. non-param bootstrapping error bars.')


# In[ ]: