<jupyter_start><jupyter_text>Sample calculations with stokes-numerics Note: This notebook needs to be able to import modules from the stokes-numerics repository.  If the notebook file is moved to a location other than the root directory of the repository, it will be necessary to adjust `sys.path` so that the imports can find stokes-numerics modules.  Setup First we set up logging.<jupyter_code>import logging
import logconfig
logconfig.logconfig(filename=None)
logconfig.loglevel(logging.WARN)<jupyter_output><empty_output><jupyter_text>Next we import the main modules.<jupyter_code>import integralequations
import framedata
import comparisons
import hkmetric<jupyter_output><empty_output><jupyter_text>Some quick IEQ computations Turn on more verbose output in this section for reassurance, since some of the commands take a little while to run.<jupyter_code>logconfig.loglevel(logging.INFO)<jupyter_output><empty_output><jupyter_text>Compute the X-functions in the (A1,A2) example, for the Hitchin section, at R=1.<jupyter_code>xar = integralequations.computeXar("A1A2", R = 1)<jupyter_output>[06-29 13:08] [MainProcess] [integralequations] [INFO] Starting to compute xars for theory A1A2, R = 1.000E+00, tolerance = 2.000E-15, steps = 131072, maxiter = 1000, L = 200, damping = 0.3
[06-29 13:08] [MainProcess] [integralequations] [INFO] IEQ: iteration 1, maxdelta = 6.928E-04 (0.5 s)
[06-29 13:08] [MainProcess] [integralequations] [INFO] IEQ: iteration 2, maxdelta = 2.074E-04 (0.5 s)
[06-29 13:08] [MainProcess] [integralequations] [INFO] IEQ: iteration 3, maxdelta = 6.208E-05 (0.5 s)
[06-29 13:08] [MainProcess] [integralequations] [INFO] IEQ: iteration 4, maxdelta = 1.858E-05 (0.5 s)
[06-29 13:08] [MainProcess] [integralequations] [INFO] IEQ: iteration 5, maxdelta = 5.563E-06 (0.5 s)
[06-29 13:08] [MainProcess] [integralequations] [INFO] IEQ: iteration 6, maxdelta = 1.665E-06 (0.5 s)
[06-29 13:08] [MainProcess] [integralequations] [INFO] IEQ: iteration 7, maxdelta = 4.985E-07 (0.5 s)
[06-29 13:08] [MainProcess] [integralequations] [INFO] IEQ: iteration 8, maxdelta = 1.492E-07 (0[...]<jupyter_text>Use the X-functions to get the spectral coordinates X1 and X2, at zeta=1.<jupyter_code>xar.computeX(charge = [1,0], zeta=1)
xar.computeX(charge = [0,1], zeta=1)<jupyter_output><empty_output><jupyter_text>Compute the X-functions at other values of zeta.<jupyter_code>xar.computeX(charge = [1,0], zeta = 3 + 5j)<jupyter_output><empty_output><jupyter_text>Compute the cluster: this is a packaging of the spectral coordinates at zeta=e^(i theta), convenient for comparing with the PDE results.<jupyter_code>xar.getCluster(theta = 0)<jupyter_output><empty_output><jupyter_text>A tricky point is that `getCluster()` is trying to compute some specific determinantal invariants, and this requires that it should return different combinations of spectral coordinates depending on the value of theta. The particular combinations of spectral coordinates which we compute are specified in the `theorydata` module. If we pick a `theta` too far from zero, then `getCluster()` usually doesn't know what combination of spectral coordinates it should compute; in that case it will return an error. In the (A1,A2) example we have 6 nonzero BPS invariants and correspondingly 6 distinguished rays in the zeta-plane, along which the integrals in the IEQ method run.<jupyter_code>xar.theory.activerays
xar.theory.activerays[0].phase
xar.theory.activerays[0].id<jupyter_output><empty_output><jupyter_text>The code stores information about X-functions on only 3 of these rays, since the other 3 are related by a symmetry.<jupyter_code>xar.theory.rayequivclasses<jupyter_output><empty_output><jupyter_text>Each of these 3 rays has a list of charges, which contains just 1 element, namely the charge gamma associated with this ray, and also has a "rayid" used to identify the ray.<jupyter_code>xar.raydata[0].charges
xar.raydata[0].rayid<jupyter_output><empty_output><jupyter_text>Choosing one of the rays, we can plot the associated X-function along the ray.<jupyter_code>import integralequationsplotting
integralequationsplotting.xarrayrayplots(xar, part="real", rayid=0, tcutoff=15)<jupyter_output><empty_output><jupyter_text>To make similar computations for opers we just add the keyword `oper=True` in various places. Here is a computation in the (A1,A2) example again. In this case we do not have to specify a parameter R, nor do we have to specify h at the time we run the iterative computation of the X-functions; we only have to specify h later, when we want to evaluate.<jupyter_code>xaroper = integralequations.computeXar("A1A2", oper = True)
xaroper.computeX(charge = [1,0], h=1)
xaroper.computeX(charge = [0,1], h=1)
integralequationsplotting.xarrayrayplots(xaroper, part="real", rayid=0, tcutoff=15)<jupyter_output><empty_output><jupyter_text>A few computations in the (A1,A3) example. For fun we choose a smaller R this time; one consequence is that we will see a plateau-like structure in the plotted X function.<jupyter_code>xar2 = integralequations.computeXar("A1A3", R = 1e-5)<jupyter_output>[06-29 13:08] [MainProcess] [integralequations] [INFO] Starting to compute xars for theory A1A3, R = 1.000E-05, tolerance = 2.000E-15, steps = 131072, maxiter = 1000, L = 200, damping = 0.3
[06-29 13:08] [MainProcess] [integralequations] [INFO] IEQ: iteration 1, maxdelta = 1.455E+00 (2.3 s)
[06-29 13:08] [MainProcess] [integralequations] [INFO] IEQ: iteration 2, maxdelta = 4.218E-01 (2.0 s)
[06-29 13:08] [MainProcess] [integralequations] [INFO] IEQ: iteration 3, maxdelta = 7.130E-02 (2.1 s)
[06-29 13:09] [MainProcess] [integralequations] [INFO] IEQ: iteration 4, maxdelta = 1.003E-02 (2.0 s)
[06-29 13:09] [MainProcess] [integralequations] [INFO] IEQ: iteration 5, maxdelta = 3.364E-03 (2.2 s)
[06-29 13:09] [MainProcess] [integralequations] [INFO] IEQ: iteration 6, maxdelta = 5.102E-04 (2.2 s)
[06-29 13:09] [MainProcess] [integralequations] [INFO] IEQ: iteration 7, maxdelta = 2.305E-04 (2.2 s)
[06-29 13:09] [MainProcess] [integralequations] [INFO] IEQ: iteration 8, maxdelta = 4.750E-05 (2[...]<jupyter_text>The plateau value in this case for X1 is just 3, as opposed to the golden ratio which we get in the (A1,A2) case (see the arXiv paper for discussion of that case).<jupyter_code>xar2.computeX(charge = [1,0,0], zeta = 1)
xar2.raydata[0].charges
integralequationsplotting.xarrayrayplots(xar2, part="real", rayid=0, tcutoff=25)<jupyter_output><empty_output><jupyter_text>Some quick ODE/PDE computations<jupyter_code>logconfig.loglevel(logging.INFO)<jupyter_output><empty_output><jupyter_text>DE computation of the subdominant sections for the (A1,A2) example at R=1.<jupyter_code>fd = framedata.computeFrames("A1A2", R = 1)<jupyter_output>[06-29 13:10] [MainProcess] [framedata] [INFO] Starting to compute frames for theory A1A2, R = 1.000E+00, theta = 0.000E+00, pde_nmesh = 511
[06-29 13:10] [MainProcess] [framedata] [INFO] Raw coefs: [[-1, 0, 0, 1]]
[06-29 13:10] [MainProcess] [framedata] [INFO] Rescaled coefs: [[(-1+0j), 0j, 0j, (1+0j)]]
[06-29 13:10] [MainProcess] [framedata] [INFO] Monicized coefs: [[(-1+0j), 0j, 0j, (1+0j)]]
[06-29 13:10] [MainProcess] [framedata] [INFO] Using rmax: 2.535
[06-29 13:10] [MainProcess] [framedata] [INFO] Error model prediction: 2.57885e-05
[06-29 13:10] [MainProcess] [nlp_fourier] [INFO] Computing Laplacian of initial guess and its transform
[06-29 13:10] [MainProcess] [nlp_fourier] [INFO] Solving PDE: 511x511 grid, thresh=1e-09
[06-29 13:10] [MainProcess] [nlp_fourier] [INFO] PDE: iter=1 L2error=0.0761667 L2delta=0	(0.08s)
[06-29 13:10] [MainProcess] [nlp_fourier] [INFO] PDE: iter=2 L2error=0.0601749 L2delta=0.000165911	(0.12s)
[06-29 13:10] [MainProcess] [nlp_fourier] [INFO] PDE: ite[...]<jupyter_text>Compute determinantal invariants from the subdominant sections.<jupyter_code>fd.getCluster()<jupyter_output><empty_output><jupyter_text>DE computation of the subdominant sections for the (A1,A2) example in the oper case, with hbar = 1.<jupyter_code>fdoper = framedata.computeFrames("A1A2", oper = True, absh = 1)<jupyter_output>[06-29 13:11] [MainProcess] [framedata] [INFO] Starting to compute frames for theory A1A2, oper, abs(h) = 1.00000000E+00, theta = 0.000E+00
[06-29 13:11] [MainProcess] [framedata] [INFO] Raw coefs: [[-1, 0, 0, 1]]
[06-29 13:11] [MainProcess] [framedata] [INFO] Rescaled coefs: [[(-1+0j), 0j, 0j, (1+0j)]]
[06-29 13:11] [MainProcess] [framedata] [INFO] Monicized coefs: [[(-1+0j), 0j, 0j, (1+0j)]]
[06-29 13:11] [MainProcess] [framedata] [INFO] Root enclosing radius: 1
[06-29 13:11] [MainProcess] [framedata] [INFO] Using rmax: 8
[06-29 13:11] [MainProcess] [framedata] [INFO] Solving ODE with initial angle = 0.2000*pi, r = 8.00
[06-29 13:11] [MainProcess] [planarode] [INFO] Calling integrator at angle 0.2000*pi, r = 8.00
[06-29 13:11] [MainProcess] [planarode] [INFO] Calling integrator at angle 0.6000*pi, r = 8.00
[06-29 13:11] [MainProcess] [planarode] [INFO] Calling integrator at angle 1.0000*pi, r = 8.00
[06-29 13:11] [MainProcess] [planarode] [INFO] Calling integrator at angle 1.4000*pi,[...]<jupyter_text>Again compute determinantal invariants from the subdominant sections.<jupyter_code>fdoper.getCluster()<jupyter_output><empty_output><jupyter_text>Comparing IEQ and DE computations Perform both the IEQ and DE computations in the (A1,A2) example, with the parameter Lambda = 0.2, at R = 0.07, and compare the results. (The parameter `scratch=True` tells `compareClusters()` not to look for saved data files but rather to compute both sides directly.)<jupyter_code>comparisons.compareClusters("A1A2_Lambda=0.2", R = 0.07, scratch = True, pde_nmesh = 511)<jupyter_output>[06-29 13:11] [MainProcess] [theorydata] [INFO] Computed turning points: [(1.0665744498072833+0j), (-0.533287224903211+0.8081991046377994j), (-0.5332872249032111-0.8081991046377994j)]
[06-29 13:11] [MainProcess] [theorydata] [INFO] Computed base sheets: [(1+0j), (-1-0j)]
[06-29 13:11] [MainProcess] [theorydata] [INFO] Computed central charges: [(-2.782474924933036+1.3077381900528278j), (-8.326672684688674e-17-2.6154763801050676j)]
[06-29 13:11] [MainProcess] [integralequations] [INFO] Starting to compute xars for theory A1A2_Lambda=0.2, R = 7.000E-02, tolerance = 2.000E-15, steps = 131072, maxiter = 1000, L = 200, damping = 0.3
[06-29 13:11] [MainProcess] [integralequations] [INFO] IEQ: iteration 1, maxdelta = 2.761E-01 (0.5 s)
[06-29 13:11] [MainProcess] [integralequations] [INFO] IEQ: iteration 2, maxdelta = 3.354E-02 (0.5 s)
[06-29 13:11] [MainProcess] [integralequations] [INFO] IEQ: iteration 3, maxdelta = 4.614E-03 (0.5 s)
[06-29 13:11] [MainProcess] [integralequations] [INFO] IEQ[...]<jupyter_text>In the output above, `xarcluster` is the IEQ result, `fdcluster` is the DE result, `reldiff` is the relative difference. The relative error in X1 is about 7e-6, which we believe is mostly due to PDE discretization error; in particular, if we pass a larger `pde_nmesh` to compareClusters, the error would likely decrease. For sufficiently large Lambda the agreement will break down; in particular this would be expected to happen if one of the periods crosses the real axis (or more generally if the phase of one the periods crosses theta), because then the formulas for the spectral coordinates as determinantal invariants would change. The IEQ side of the above computation uses period integrals which are computed numerically. The commands below compute just the periods.<jupyter_code>import theory
t = theory.theory("A1A2_Lambda=0.2")
t.Z(charge = [1,0])
t.Z(charge = [0,1])<jupyter_output><empty_output><jupyter_text>Perform both the IEQ and DE computations in the (A2,A1) example, oper case, with the parameter c = 0.3 + 0.2i, at hbar = 1. (This example is interesting because it involves both the quadratic and cubic differentials nonzero. Again, we expect it to work only for small enough c.)<jupyter_code>comparisons.compareClusters("A2A1_c=0.3+0.2j", oper = True, scratch = True)
comparisons.compareClusters("A1A3", scratch = True)<jupyter_output>[06-29 13:14] [MainProcess] [integralequations] [INFO] Starting to compute xars for theory A1A3, R = 1.000E+00, tolerance = 2.000E-15, steps = 131072, maxiter = 1000, L = 200, damping = 0.3
[06-29 13:14] [MainProcess] [integralequations] [INFO] IEQ: iteration 1, maxdelta = 4.282E-03 (2.1 s)
[06-29 13:14] [MainProcess] [integralequations] [INFO] IEQ: iteration 2, maxdelta = 1.277E-03 (2.1 s)
[06-29 13:14] [MainProcess] [integralequations] [INFO] IEQ: iteration 3, maxdelta = 3.807E-04 (1.9 s)
[06-29 13:14] [MainProcess] [integralequations] [INFO] IEQ: iteration 4, maxdelta = 1.135E-04 (2.1 s)
[06-29 13:14] [MainProcess] [integralequations] [INFO] IEQ: iteration 5, maxdelta = 3.384E-05 (2.1 s)
[06-29 13:14] [MainProcess] [integralequations] [INFO] IEQ: iteration 6, maxdelta = 1.009E-05 (2.1 s)
[06-29 13:15] [MainProcess] [integralequations] [INFO] IEQ: iteration 7, maxdelta = 3.009E-06 (2.1 s)
[06-29 13:15] [MainProcess] [integralequations] [INFO] IEQ: iteration 8, maxdelta = 8.972E-07 (2[...]<jupyter_text>Some quick metric computations<jupyter_code>hkmetric.comparemetrics(c = 0.5, Lambda = 0, pde_nmesh = 500)
hkmetric.comparemetrics(c = 0.5, Lambda = 0.5, pde_nmesh = 500)<jupyter_output>[06-29 13:20] [MainProcess] [hkmetric] [INFO] Starting PDE computation of hyperkahler metric element
[06-29 13:20] [MainProcess] [hkmetric] [INFO] coefs: [-0.5, -0.5, 0, 1]
[06-29 13:20] [MainProcess] [hkmetric] [INFO] Monicized coefs: [-0.5, -0.5, 0.0, 1.0]
[06-29 13:20] [MainProcess] [discdiff] [INFO] Loading discrete laplacian from grid500.v2.disclaplacian.cache.npz
[06-29 13:20] [MainProcess] [nlp_euler] [INFO] Solving PDE
[06-29 13:20] [MainProcess] [nlp_euler] [INFO] PDE: GOAL=5e-11
[06-29 13:21] [MainProcess] [nlp_euler] [INFO] PDE: iter=1 error=17.9155 delta=0.247435	(20.86s)
[06-29 13:21] [MainProcess] [nlp_euler] [INFO] PDE: iter=2 error=0.361531 delta=0.022183	(22.35s)
[06-29 13:22] [MainProcess] [nlp_euler] [INFO] PDE: iter=3 error=0.00162679 delta=9.20838e-05	(22.45s)
[06-29 13:22] [MainProcess] [nlp_euler] [INFO] PDE: iter=4 error=2.75986e-08 delta=1.51395e-09	(23.55s)
[06-29 13:22] [MainProcess] [nlp_euler] [INFO] PDE: iter=5 error=1.31877e-11 delta=6.36145e-16	(21.93s)
[...]<jupyter_text>Looking at stored data<jupyter_code>logconfig.loglevel(logging.WARN)<jupyter_output><empty_output><jupyter_text>The commands in this section use the stored data files; these files are available as the dataset at [https://doi.org/10.7910/DVN/W0V4D9](https://doi.org/10.7910/DVN/W0V4D9). The code expects them to be installed at paths specified in `paths.conf`. The command below interrogates the code about what paths it is looking in.<jupyter_code>(integralequations.fileinfo.XARPATH, integralequations.fileinfo.FRAMEPATH, integralequations.fileinfo.HKMETRICPATH, integralequations.fileinfo.CACHEPATH)<jupyter_output><empty_output><jupyter_text>A comparison between the two stored clusters in the (A1,A2) example, at R=1 and theta=0.1.`xarcluster` is the IEQ result, `fdcluster` is the DE result, `reldiff` is the relative difference.<jupyter_code>comparisons.compareClusters("A1A2", oper=False, R=1, pde_nmesh=8191, theta = 0.1)<jupyter_output><empty_output><jupyter_text>All the clusters computed by IEQ method in the (A1,A2) example. The list of R values is drawn from the `theorydata` module; if the data file for some R were missing, this command would generate an error.<jupyter_code>comparisons.xarclustertable("A1A2", theta = 0.1)<jupyter_output><empty_output><jupyter_text>All the clusters computed by DE method in the (A1,A2) example. Since data is stored for several values of `pde_nmesh` we have to specify which we want. (Otherwise the default will be used, which is stored in framedata.PDE_NMESH; when the framedata module is loaded, this is set to 511; we actually don't have data stored for that value, so it would generate an error here.)<jupyter_code>comparisons.fdclustertable("A1A2", theta = 0.1, pde_nmesh = 8191)<jupyter_output><empty_output><jupyter_text>A quick and dirty comparison table.<jupyter_code>print(comparisons.plaintables("A1A2", theta = 0.1, pde_nmesh = 8191, precision=5)["data"])<jupyter_output>R              X1f             X1i          X1logdiff          X2f             X2i          X2logdiff    
 1.60523E-09    -6.18034E-01    -6.18034E-01     9.88098E-15    -1.00000E+00    -1.00000E+00    -1.33227E-14    
 2.33559E-09    -6.18034E-01    -6.18034E-01    -2.33147E-15    -1.00000E+00    -1.00000E+00    -5.32907E-15    
 3.39827E-09    -6.18034E-01    -6.18034E-01    -7.27751E-14    -1.00000E+00    -1.00000E+00    -1.43219E-14    
 4.94445E-09    -6.18034E-01    -6.18034E-01    -8.77076E-15    -1.00000E+00    -1.00000E+00     6.88338E-15    
 7.19413E-09    -6.18034E-01    -6.18034E-01    -9.10938E-14    -1.00000E+00    -1.00000E+00    -1.33227E-15    
 1.04674E-08    -6.18034E-01    -6.18034E-01     3.35842E-14    -1.00000E+00    -1.00000E+00    -3.07532E-14    
 1.52300E-08    -6.18034E-01    -6.18034E-01    -2.37033E-14    -1.00000E+00    -1.00000E+00    -1.40998E-14    
 2.21595E-08    -6.18034E-01    -6.18034E-01    -8.51541E-14    -1.00000E+00    -1.00000E+00    [...]
{'fraction_non_alphanumeric': 0.1292177246065393, 'fraction_numerical': 0.11568616063650619, 'mean_word_length': 5.209880995311937, 'pattern_counts': {'":': 0, '<': 109, '<?xml version=': 0, '>': 109, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '29', 'dir': 'jupyter-structured-clean-dedup', 'id': '513116', 'license': 'MIT', 'n_tokens_mistral': 7638, 'n_tokens_neox': 6392, 'n_words': 1846, 'path': 'sample-computations.ipynb', 'repo_name': 'neitzke/stokes-numerics', 'size': '17219'}