View source for ALPS 2 Tutorial:MC-09 Directed Worm Algorithm:Intensive Example 2

= Description =

This example is about error budgeting.

'''This section is under construction.'''

= Topic 1: Error budgeting in total particle number =

We want to investigate how the condensate fraction and visibility varies with up to 5% fluctuation in the total particle number.

== A first reference to start the business ==

The usual headers:

 import numpy;
 import pyalps;
 import pyalps.dwa;

Parameters:

 tof_phase = pyalps.dwa.tofPhase(time_of_flight=15.5, wavelength=[765,843,843], mass=86.99)

 params=[]

 params.append(
   { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
   , 'MODEL'           : 'boson Hubbard'                        
   , 'L'               : 160                             
   , 'Nmax'            : 20                 
   , 't'               : 1.                 
   , 'U'               : 8.11               
   , 'T'               : 1.                 
   , 'mu_homogeneous'  : 4.05                
   , 'mu'              : 'mu_homogeneous - (0.0073752*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.0036849*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.0039068155*(z-(L-1)/2.)*(z-(L-1)/2.))'
   , 'tof_phase'       : str(tof_phase) 
   , 'SWEEPS'          : 100000             
   , 'SKIP'            : 100                 
   }
 )

Write:

 h5_infiles = pyalps.writeInputH5Files("parm.errorbudget.N",params);

Or read:

 h5_infiles = pyalps.getInputH5Files(prefix = "parm.errorbudget.N");

A trial run of course:

 pyalps.runApplication('dwa', h5_infiles[0]);

== Pinning chemical potential ==

=== 5% fluctuation in chemical potential? ===

The usual headers:

 import numpy;
 import pyalps;
 import pyalps.dwa;

Write from scratch:

 tof_phase = pyalps.dwa.tofPhase(time_of_flight=15.5, wavelength=[765,843,843], mass=86.99)
 
 params=[]
 for ratio in [0.95, 0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04, 1.05]:
   params.append(
     { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
     , 'MODEL'           : 'boson Hubbard'                        
     , 'L'               : 160                             
     , 'Nmax'            : 20               
     , 't'               : 1.              
     , 'U'               : 8.11               
     , 'T'               : 1.                 
     , 'mu_homogeneous'  : ratio * 4.05                
     , 'mu'              : 'mu_homogeneous - (0.0073752*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.0036849*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.0039068155*(z-(L-1)/2.)*(z-(L-1)/2.))'
     , 'tof_phase'       : str(tof_phase) 
     , 'SWEEPS'          : 100000             
     , 'SKIP'            : 100                 
     }
   )
 
 h5_infiles = pyalps.writeInputH5Files('parm.errorbudget.N',params);

Read from archive:

 h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.N');

Check for yourself, for instance:

 pyalps.getParameters(h5_infiles[0]);

=== <div id="ex2_automation">Automated processing</div> ===

 def parallel_automated_processing(taskfile):
   pyalps.dwa.recursiveRun( 
     "pyalps.runApplication('dwa', taskfile, T=2700)" ,  cmd_lang = 'python' ,
     follow_up_script = "pyalps.dwa.extract_worldlines(infile=pyalps.input2output(taskfile), outfile=taskfile)" ,
     break_if = "pyalps.dwa.thermalized(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
     write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Thermalizing')" ,
     loc = locals() ,
     batch_submit = True ,
     batch_cmd_prefix = 'bsub -J jobN -n1 -W60' ,
     batch_run_script = taskfile + '.thermalize.script' ,
     batch_next_run_script = taskfile + '.converge.script' ,
     batch_noRun = True
   );
   
   pyalps.dwa.recursiveRun( 
     "pyalps.runApplication('dwa', pyalps.input2output(taskfile), T=2700)" ,  cmd_lang = 'python' ,
     break_if = "pyalps.dwa.converged(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
     write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Converging')" ,
     loc = locals() ,
     batch_submit = True ,
     batch_cmd_prefix = 'bsub -J jobN -n1 -W60' ,
     batch_run_script = taskfile + '.converge.script' ,
     batch_noRun = True
   );
   
   pyalps.dwa.startRunScript(
     batch_run_script= taskfile + '.thermalize.script' ,
     batch_cmd_prefix = 'bsub -J jobN -n1 -W60' ,
     loc = locals()
   );

 for h5_infile in h5_infiles:
   parallel_automated_processing(h5_infile);

=== Observing on your own ===

First the headers again:

 import matplotlib;
 matplotlib.use('Agg');

 import pyalps;
 import pyalps.plot;
 import pyalps.dwa;

We collect useful data:

 resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.N');
 data = pyalps.loadMeasurements(resultFiles, 'Total Particle Number');
 N = pyalps.collectXY(data, x='mu_homogeneous', y='Total Particle Number')

then plot it:

 fig = matplotlib.pyplot.figure();
 pyalps.plot.plot(N);
 fig.savefig('parm.errorbudget.N.N.jpg');
 fig.show();

or send it to yourself:

 pyalps.sendmail( 'pingnang@phys.ethz.ch'    # email address of recipients 
                , attachment='parm.errorbudget.N.N.jpg'
                , subject='Picture: parm.errorbudget.N.N.jpg'
                );

=== Correcting chemical potential for fixed N ===

First the headers again:

 import numpy;
 import pyalps;
 import pyalps.dwa;

Easily done like this:

 resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.N');
 data = pyalps.loadMeasurements(resultFiles, 'Total Particle Number');
 params = pyalps.paramsAtFixedY(data, x='mu_homogeneous', y='Total Particle Number', fixedY=280000*numpy.array([0.96, 0.97, 0.98, 0.99, 1., 1.01, 1.02, 1.03, 1.04]));

Prepare new input files:

 h5_infiles = pyalps.writeInputH5Files("parm.errorbudget.N.1", params);

Run them:

 for h5_infile in h5_infiles:
   parallel_automated_processing(h5_infile);

To check whether they are thermalized:

 resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.N.1');
 for resultFile in resultFiles:
   mu=pyalps.getParameters(resultFile)[0]['mu_homogeneous'];
   N=pyalps.getMeasurements(resultFile,observable='Total Particle Number')[0]['mean']['value'];
   thermalized=pyalps.dwa.thermalized(resultFile, 'Total Particle Number', simplified=True);
   print('mu : ' + str(mu) + ' , N : ' + str(N) + ' , thermalized : ' + str(thermalized));

== Condensate fraction ==

=== Turn on respective measurements ===

Gather the files:

 h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.N.1')

Switch on local density measurements:

 for h5_infile in h5_infiles:
   pyalps.dwa.switchParameter(h5_infile, 'MEASURE[Local Density]', 1);
   pyalps.dwa.switchParameter(pyalps.input2output(h5_infile), 'MEASURE[Local Density]', 1);

Switch on green function measurements:

 for h5_infile in h5_infiles:
   pyalps.dwa.switchParameter(h5_infile, 'MEASURE[Green Function]', 1);
   pyalps.dwa.switchParameter(pyalps.input2output(h5_infile), 'MEASURE[Green Function]', 1);

=== Automated processing ===

 def parallel_automatic_processing(taskfile):
   pyalps.dwa.recursiveRun( 
     "pyalps.runApplication('dwa', taskfile, T=2700)" ,  cmd_lang = 'python' ,
     write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Measuring')" ,
     loc = locals() ,
     batch_submit = True ,
     batch_cmd_prefix = 'bsub -J jobN -n1 -W60' ,
     batch_run_script = taskfile + '.measure1.script' ,
     batch_next_run_script = taskfile + '.measure.script' ,
     batch_noRun = True
   ); 
   
   pyalps.dwa.recursiveRun( 
     "pyalps.runApplication('dwa', pyalps.input2output(taskfile), T=2700)" ,  cmd_lang = 'python' ,
     n = 20 ,
     write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Measuring')" ,
     loc = locals() ,
     batch_submit = True ,
     batch_cmd_prefix = 'bsub -J jobN -n1 -W60' ,
     batch_run_script = taskfile + '.measure.script' ,
     batch_next_run_script = taskfile + '.sendmail.script' ,
     batch_noRun = True
   ); 
   
   pyalps.dwa.recursiveRun(
      "pyalps.sendmail('pingnang@phys.ethz.ch', message= taskfile + ' has finished measurements for 20 iterations.')" ,
      cmd_lang = 'python' ,
      loc = locals() ,
      batch_submit = True ,
      batch_run_script = taskfile + '.sendmail.script' ,
      batch_noRun = True
   ); 
   
   pyalps.dwa.startRunScript(
     batch_run_script= taskfile + '.measure1.script' ,
     batch_cmd_prefix = 'bsub -J jobN -n1 -W60' ,
     loc = locals()
   );

 for h5_infile in h5_infiles:
   parallel_automatic_processing(h5_infile);

=== Extracting final data(plot) ===

First, the headers again:

 import matplotlib; 
 matplotlib.use('Agg');

 import matplotlib.pyplot;
 from mpl_toolkits.mplot3d.axes3d import Axes3D;
 import numpy;
 import numpy.fft;
 import pyalps;
 import pyalps.dwa;

Then, evaluate visibility:

 resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.N.1');

 results = []; 
 for resultFile in resultFiles:
   L = int(pyalps.getParameters(resultFile)[0]['L']);
   green_tof = pyalps.getMeasurements(resultFile, observable='Green Function:TOF')[0]['mean']['value'].reshape([L,L,L]);
   momentum_density = numpy.fft.fftn(green_tof).real; 
   
   V0   = numpy.array([8.805, 8. ,  8. ]);
   wlen = numpy.array([765., 843., 843.]);
   a    = 101;
   m    = 86.99;
   
   band = pyalps.dwa.bandstructure(V0, wlen, a, m, L);
   
   q_z   = numpy.array(band.q(2));
   wk2_z = numpy.array(band.wk2(2));
   wk2_z = numpy.array([q_z, wk2_z]).transpose();
   wk2_z0= wk2_z[wk2_z[:,0] == 0.][0][1];
   wk2_z = numpy.transpose(wk2_z[wk2_z[:,1]/wk2_z0 > 1e-4]);
   q_z   = wk2_z[0]
   wk2_z = wk2_z[1]
   
   momentum_density = numpy.tile(momentum_density, reps=(1,1,2*((q_z[q_z[:] >= 0].size / L)+1)));
   
   dummy = numpy.zeros(momentum_density.shape[2]);
   dummy[dummy.size/2:dummy.size/2 + q_z[q_z[:] >= 0].size] = wk2_z[q_z[:] >= 0]
   dummy[dummy.size/2-q_z[q_z[:] < 0].size:dummy.size/2]  = wk2_z[q_z[:] < 0]    
   
   momentum_density = numpy.tensordot(momentum_density, dummy, axes=([2],[0]))
   momentum_density = numpy.tile(momentum_density, reps=(4,4));
   
   q_x   = numpy.array(band.q(0));
   wk2_x = numpy.array(band.wk2(0));
   wk2_x = wk2_x[(q_x[:] >= -2.)*(q_x[:] < 2.)];
   q_x   = q_x[(q_x[:] >= -2.)*(q_x[:] < 2.)] 
   
   q_y   = numpy.array(band.q(1));
   wk2_y = numpy.array(band.wk2(1));
   wk2_y = wk2_y[(q_y[:] >= -2.)*(q_y[:] < 2.)];
   q_y   = q_y[(q_y[:] >= -2.)*(q_y[:] < 2.)]
   
   wk2 = numpy.outer(wk2_x, wk2_y);
   
   q_x = numpy.array([q_x] * wk2.shape[1], float);
   q_y = numpy.array([q_y] * wk2.shape[0], float).transpose();
   tof_image = wk2 * momentum_density;
   
   selection_max = numpy.outer([(q_x[0,:] >= -0.05) * (q_x[0,:] <= 0.05)] , [(q_y[:,0] >= -0.05) * (q_y[:,0] <= 0.05)])
   selection_min = numpy.outer([(q_x[0,:] >= 0.45) * (q_x[0,:] <= 0.55)] , [(q_y[:,0] >= 0.45) * (q_y[:,0] <= 0.55)])                   
   
   nk_max = sum(tof_image[selection_max])
   nk_min = sum(tof_image[selection_min])
   
   condensate_fraction = nk_max;
   visibility = (nk_max - nk_min) / (nk_max + nk_min);
     
   N = pyalps.getMeasurements(resultFile, observable='Total Particle Number')[0]['mean']['value'];
   
   results.append([N, condensate_fraction]);

'''This section is under construction.'''

= Topic 2: Error budgeting in scattering length =

We want to investigate how the condensate fraction and visibility varies with up to 5% fluctuation in the scattering length.

== Varying scattering length at fixed chemical potential ==

Let's first fix the chemical potential while varying scattering length, and see what happens out of it.

=== Parameters ===

The usual headers:

 import numpy;
 import pyalps;
 import pyalps.dwa;

Write from scratch:

 tof_phase = pyalps.dwa.tofPhase(time_of_flight=15.5, wavelength=[765,843,843], mass=86.99)
 
 ps = [];
 for ratio in [0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04]:
   V0   = numpy.array([8.805, 8. ,  8. ]);
   wlen = numpy.array([765., 843., 843.]);
   a    = ratio*101;
   m    = 86.99;
   L    = 160;
  
   band = pyalps.dwa.bandstructure(V0, wlen, a, m, L);
   p = {};
   p.update({'ratio_a' : ratio , 'U' : band.Ut()[2]})
   ps.append(p);
 
 params=[]
 for p in ps:
   params.append(
     { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
     , 'MODEL'           : 'boson Hubbard'                        
     , 'L'               : 160                             
     , 'Nmax'            : 20
     , 'ratio_a'         : p['ratio_a']                 
     , 't'               : 1.          
     , 'U'               : p['U']               
     , 'T'               : 1.                 
     , 'mu_homogeneous'  : 4.0265                
     , 'mu'              : 'mu_homogeneous - (0.0073752*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.0036849*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.0039068155*(z-(L-1)/2.)*(z-(L-1)/2.))'
     , 'tof_phase'       : str(tof_phase) 
     , 'SWEEPS'          : 100000             
     , 'SKIP'            : 100                 
     }
   )
 
 h5_infiles = pyalps.writeInputH5Files("parm.errorbudget.a",params);

Read from archive:

 h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.a');

Check for yourself, for instance:

 pyalps.getParameters(h5_infiles[4]);

Have a first taste of simulation:

 pyalps.runApplication('dwa', h5_infiles[4], T=2700);

=== Automated processing ===

 def parallel_automated_processing(taskfile):
   pyalps.dwa.recursiveRun( 
     "pyalps.runApplication('dwa', taskfile, T=2700)" ,  cmd_lang = 'python' ,
     follow_up_script = "pyalps.dwa.extract_worldlines(infile=pyalps.input2output(taskfile), outfile=taskfile)" ,
     break_if = "pyalps.dwa.thermalized(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
     write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Thermalizing')" ,
     loc = locals() ,
     batch_submit = True ,
     batch_cmd_prefix = 'bsub -J joba -n1 -W60' ,
     batch_run_script = taskfile + '.thermalize.script' ,
     batch_next_run_script = taskfile + '.converge.script' ,
     batch_noRun = True
   );
   
   pyalps.dwa.recursiveRun( 
     "pyalps.runApplication('dwa', pyalps.input2output(taskfile), T=2700)" ,  cmd_lang = 'python' ,
     break_if = "pyalps.dwa.converged(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
     write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Converging')" ,
     loc = locals() ,
     batch_submit = True ,
     batch_cmd_prefix = 'bsub -J joba -n1 -W60' ,
     batch_run_script = taskfile + '.converge.script' ,
     batch_noRun = True
   );
   
   pyalps.dwa.startRunScript(
     batch_run_script= taskfile + '.thermalize.script' ,
     batch_cmd_prefix = 'bsub -J joba -n1 -W60' ,
     loc = locals()
   );

 for h5_infile in h5_infiles:
   parallel_automated_processing(h5_infile);

=== Observing on your own ===

First the headers again:

 import matplotlib;
 matplotlib.use('Agg');

 import pyalps;
 import pyalps.plot;
 import pyalps.dwa;

We collect useful data:

 resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.a');
 data = pyalps.loadMeasurements(resultFiles, 'Total Particle Number');
 N = pyalps.collectXY(data, x='ratio_a', y='Total Particle Number')

then plot it:

 fig = matplotlib.pyplot.figure();
 pyalps.plot.plot(N);
 fig.savefig('parm.errorbudget.a.N.jpg');
 fig.show();

or send it to yourself:

 pyalps.sendmail( 'pingnang@phys.ethz.ch'    # email address of recipients 
                , attachment='parm.errorbudget.a.N.jpg'
                , subject='Picture: parm.errorbudget.a.N.jpg'
                );

== Varying scattering length at fixed N ==

=== Parameters ===

The usual headers:

 import numpy;
 import pyalps;
 import pyalps.dwa;

Write from scratch:

 tof_phase = pyalps.dwa.tofPhase(time_of_flight=15.5, wavelength=[765,843,843], mass=86.99)
 
 ps = [];
 for ratio in [0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04]:
   V0   = numpy.array([8.805, 8. ,  8. ]);
   wlen = numpy.array([765., 843., 843.]);
   a    = ratio*101;
   m    = 86.99;
   L    = 160;
  
   band = pyalps.dwa.bandstructure(V0, wlen, a, m, L);
   p = {};
   p.update({'ratio_a' : ratio , 'U' : band.Ut()[2]})
   ps.append(p);
 
 params=[]
 for p in ps:
   for ratio in [0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04]:
     params.append(
       { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
       , 'MODEL'           : 'boson Hubbard'                        
       , 'L'               : 160                             
       , 'Nmax'            : 20
       , 'ratio_a'         : p['ratio_a']                 
       , 't'               : 1.          
       , 'U'               : p['U']               
       , 'T'               : 1.                 
       , 'mu_homogeneous'  : ratio*4.0265                
       , 'mu'              : 'mu_homogeneous - (0.0073752*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.0036849*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.0039068155*(z-(L-1)/2.)*(z-(L-1)/2.))'
       , 'tof_phase'       : str(tof_phase) 
       , 'SWEEPS'          : 100000             
       , 'SKIP'            : 100                 
       }
     )
 
 h5_infiles = pyalps.writeInputH5Files("parm.errorbudget.a",params);

Read from archive:

 h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.a');

=== Observing on your own ===

First the headers again:

 import matplotlib;
 matplotlib.use('Agg');

 import pyalps;
 import pyalps.plot;
 import pyalps.dwa;

We collect useful data:

 resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.a');
 data = pyalps.loadMeasurements(resultFiles, 'Total Particle Number');
 N = pyalps.collectXY(data, x='mu_homogeneous', y='Total Particle Number', foreach=['ratio_a'])

then plot it:

 fig = matplotlib.pyplot.figure();
 pyalps.plot.plot(N);
 fig.savefig('parm.errorbudget.a.N.jpg');
 fig.show();

or send it to yourself:

 pyalps.sendmail( 'pingnang@phys.ethz.ch'    # email address of recipients 
                , attachment='parm.errorbudget.a.N.jpg'
                , subject='Picture: parm.errorbudget.a.N.jpg'
                );

=== Fine tuning chemical potential ===

First the headers again:

 import numpy;
 import pyalps;
 import pyalps.dwa;

Easily done like this:

 resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.a');
 data = pyalps.loadMeasurements(resultFiles, 'Total Particle Number');
 params = pyalps.paramsAtFixedY(data, x='mu_homogeneous', y='Total Particle Number', foreach=['ratio_a'], fixedY=280000);

Prepare new input files:

 h5_infiles = pyalps.writeInputH5Files("parm.errorbudget.a.1", params);

Run them:

 for h5_infile in h5_infiles:
   parallel_automated_processing(h5_infile);

To check whether they are thermalized:

 resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.a.1');
 for resultFile in resultFiles:
   ratio_a = pyalps.getParameters(resultFile)[0]['ratio_a'];
   mu=pyalps.getParameters(resultFile)[0]['mu_homogeneous'];
   N=pyalps.getMeasurements(resultFile,observable='Total Particle Number')[0]['mean']['value'];
   thermalized=pyalps.dwa.thermalized(resultFile, 'Total Particle Number', simplified=True);
   print('ratio_a : ' + str(ratio_a) + ' , mu : ' + str(mu) + ' , N : ' + str(N) + ' , thermalized : ' + str(thermalized));

== Condensate fraction ==

=== Turn on respective measurements ===

Gather the files:

 h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.a.1')

Switch on local density measurements:

 for h5_infile in h5_infiles:
   pyalps.dwa.switchParameter(h5_infile, 'MEASURE[Local Density]', 1);
   pyalps.dwa.switchParameter(pyalps.input2output(h5_infile), 'MEASURE[Local Density]', 1);

Switch on green function measurements:

 for h5_infile in h5_infiles:
   pyalps.dwa.switchParameter(h5_infile, 'MEASURE[Green Function]', 1);
   pyalps.dwa.switchParameter(pyalps.input2output(h5_infile), 'MEASURE[Green Function]', 1);

=== Automated processing ===

 def parallel_automatic_processing(taskfile):
   pyalps.dwa.recursiveRun( 
     "pyalps.runApplication('dwa', taskfile, T=2700)" ,  cmd_lang = 'python' ,
     write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Measuring')" ,
     loc = locals() ,
     batch_submit = True ,
     batch_cmd_prefix = 'bsub -J joba -n1 -W60' ,
     batch_run_script = taskfile + '.measure1.script' ,
     batch_next_run_script = taskfile + '.measure.script' ,
     batch_noRun = True
   ); 
   
   pyalps.dwa.recursiveRun( 
     "pyalps.runApplication('dwa', pyalps.input2output(taskfile), T=2700)" ,  cmd_lang = 'python' ,
     n = 20 ,
     write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Measuring')" ,
     loc = locals() ,
     batch_submit = True ,
     batch_cmd_prefix = 'bsub -J joba -n1 -W60' ,
     batch_run_script = taskfile + '.measure.script' ,
     batch_next_run_script = taskfile + '.sendmail.script' ,
     batch_noRun = True
   ); 
   
   pyalps.dwa.recursiveRun(
      "pyalps.sendmail('pingnang@phys.ethz.ch', message= taskfile + ' has finished measurements for 20 iterations.')" ,
      cmd_lang = 'python' ,
      loc = locals() ,
      batch_submit = True ,
      batch_run_script = taskfile + '.sendmail.script' ,
      batch_noRun = True
   ); 
   
   pyalps.dwa.startRunScript(
     batch_run_script= taskfile + '.measure1.script' ,
     batch_cmd_prefix = 'bsub -J joba -n1 -W60' ,
     loc = locals()
   );

 for h5_infile in h5_infiles:
   parallel_automatic_processing(h5_infile);

=== Extracting final data(plot) ===

'''This section is under construction.'''

= Topic 3: Error budgeting in lattice potential =

We want to investigate how the condensate fraction and visibility varies with up to 5% fluctuation in the lattice potential along the x-direction.

=== Parameters ===

The usual headers:

 import numpy;
 import pyalps;
 import pyalps.dwa;

Write from scratch:

 tof_phase = pyalps.dwa.tofPhase(time_of_flight=15.5, wavelength=[765,843,843], mass=86.99)
 
 ps = [];
 for ratio in [0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04]:
   V0   = numpy.array([ratio*8.805, 8. ,  8. ]);
   wlen = numpy.array([765., 843., 843.]);
   a    = 101;
   m    = 86.99;
   L    = 160;
  
   band = pyalps.dwa.bandstructure(V0, wlen, a, m, L);
   p = {};
   p.update({'ratio_V0x' : ratio , 'U' : band.Ut()[2] , 'tx_t' : band.t()[0]/band.t()[2]})
   ps.append(p);
 
 params=[]
 for p in ps:
   for ratio in [0.8, 0.82, 0.84, 0.88, 0.9, 0.92, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.08, 1.10, 1.12, 1.14, 1.16, 1.18, 1.2]:
     params.append(
       { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
       , 'MODEL'           : 'boson Hubbard'                        
       , 'L'               : 160                             
       , 'Nmax'            : 20
       , 'ratio_V0x'       : p['ratio_V0x']                 
       , 't'               : 1.
       , 'tx_t'            : p['tx_t']                 
       , 'U'               : p['U']               
       , 'T'               : 1.                 
       , 'mu_homogeneous'  : ratio*4.0265                
       , 'mu'              : 'mu_homogeneous - (0.0073752*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.0036849*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.0039068155*(z-(L-1)/2.)*(z-(L-1)/2.))'
       , 'tof_phase'       : str(tof_phase) 
       , 'SWEEPS'          : 100000             
       , 'SKIP'            : 100                 
       }
     )
 
 h5_infiles = pyalps.writeInputH5Files("parm.errorbudget.V0",params);

Read from archive:

 h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.V0');

=== Automated processing ===

 def parallel_automated_processing(taskfile):
   pyalps.dwa.recursiveRun( 
     "pyalps.runApplication('dwa', taskfile, T=2700)" ,  cmd_lang = 'python' ,
     follow_up_script = "pyalps.dwa.extract_worldlines(infile=pyalps.input2output(taskfile), outfile=taskfile)" ,
     break_if = "pyalps.dwa.thermalized(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
     write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Thermalizing')" ,
     loc = locals() ,
     batch_submit = True ,
     batch_cmd_prefix = 'bsub -J jobV0 -n1 -W60' ,
     batch_run_script = taskfile + '.thermalize.script' ,
     batch_next_run_script = taskfile + '.converge.script' ,
     batch_noRun = True
   );
   
   pyalps.dwa.recursiveRun( 
     "pyalps.runApplication('dwa', pyalps.input2output(taskfile), T=2700)" ,  cmd_lang = 'python' ,
     break_if = "pyalps.dwa.converged(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
     write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Converging')" ,
     loc = locals() ,
     batch_submit = True ,
     batch_cmd_prefix = 'bsub -J jobV0 -n1 -W60' ,
     batch_run_script = taskfile + '.converge.script' ,
     batch_noRun = True
   );
   
   pyalps.dwa.startRunScript(
     batch_run_script= taskfile + '.thermalize.script' ,
     batch_cmd_prefix = 'bsub -J jobV0 -n1 -W60' ,
     loc = locals()
   );

 for h5_infile in h5_infiles:
   parallel_automated_processing(h5_infile);

=== Unexpected norun ===

Sometimes, due to certain node failure, jobs may not even start. The following is the way to fix this:

 import pyalps;
 import pyalps.dwa;

 all_h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.V0');
 h5_infiles = [];
 for h5_infile in all_h5_infiles:
    ar = pyalps.hdf5.h5ar(h5_infile);
    if not ar.is_group('/simulation/worldlines'):
       h5_infiles.append(h5_infile);

 for h5_infile in h5_infiles:
   parallel_automated_processing(h5_infile);

=== Unexpected worldlines configuration corruption ===

You may check worldlines configuration like the following:

 import pyalps;
 import pyalps.dwa;

Checking infiles:

 all_h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.V0');
 h5_infiles = [];
 for h5_infile in all_h5_infiles:
    ar = pyalps.hdf5.h5ar(h5_infile);
    if ar.is_group('/simulation/worldlines'):
      wl = pyalps.dwa.extract_worldlines(h5_infile);
      valid = wl.is_valid(20);
      if not valid:
        h5_infiles.append(h5_infile);
        print('h5_infile : ' + h5_infile + ' , valid : ' + str(valid));

Checking outfiles:

 all_h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.V0');
 h5_infiles = [];
 for h5_infile in all_h5_infiles:
    ar = pyalps.hdf5.h5ar(h5_infile);
    if ar.is_group('/simulation/worldlines'):
      wl = pyalps.dwa.extract_worldlines(pyalps.input2output(h5_infile));
      valid = wl.is_valid(20);
      if not valid:
        h5_infiles.append(h5_infile);
        print('h5_infile : ' + h5_infile + ' , valid : ' + str(valid));

=== Checking thermalization ===

 resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.V0');
 for resultFile in resultFiles:
   ratio_V0x = pyalps.getParameters(resultFile)[0]['ratio_V0x'];
   mu=pyalps.getParameters(resultFile)[0]['mu_homogeneous'];
   N=pyalps.getMeasurements(resultFile,observable='Total Particle Number')[0]['mean']['value'];
   thermalized=pyalps.dwa.thermalized(resultFile, 'Total Particle Number', simplified=True);
   print(resultFile + ' , ratio_V0x : ' + str(ratio_V0x) + ' , mu : ' + str(mu) + ' , N : ' + str(N) + ' , thermalized : ' + str(thermalized));

=== Observing on your own ===

First the headers again:

 import matplotlib;
 matplotlib.use('Agg');

 import pyalps;
 import pyalps.plot;
 import pyalps.dwa;

We collect useful data:

 resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.V0');
 data = pyalps.loadMeasurements(resultFiles, 'Total Particle Number');
 N = pyalps.collectXY(data, x='mu_homogeneous', y='Total Particle Number', foreach=['ratio_V0x'])

then plot it:

 fig = matplotlib.pyplot.figure();
 pyalps.plot.plot(N);
 fig.savefig('parm.errorbudget.V0.N.jpg');
 fig.show();

or send it to yourself:

 pyalps.sendmail( 'pingnang@phys.ethz.ch'    # email address of recipients 
                , attachment='parm.errorbudget.V0.N.jpg'
                , subject='Picture: parm.errorbudget.V0.N.jpg'
                );

=== Fine tuning the chemical potential ===

First the headers again:

 import numpy;
 import pyalps;
 import pyalps.dwa;

Easily done like this:

 resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.V0');
 data = pyalps.loadMeasurements(resultFiles, 'Total Particle Number');
 params = pyalps.paramsAtFixedY(data, x='mu_homogeneous', y='Total Particle Number', foreach=['ratio_V0x'], fixedY=280000);

Prepare new input files:

 h5_infiles = pyalps.writeInputH5Files("parm.errorbudget.V0.2", params);

Run them:

 for h5_infile in h5_infiles:
   parallel_automated_processing(h5_infile);

= Topic 4: Error budgeting in laser wavelength =

We want to investigate how the condensate fraction and visibility varies with up to 5% fluctuation in the laser wavelength along the x-direction.

=== Parameters ===

The usual headers:

 import numpy;
 import pyalps;
 import pyalps.dwa;

Write from scratch:

 ps = [];
 for ratio in [0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04]:
   V0   = numpy.array([ratio*ratio*8.805, 8. ,  8. ]);
   wlen = numpy.array([ratio*765., 843., 843.]);
   a    = 101;
   m    = 86.99;
   L    = 160;
  
   band = pyalps.dwa.bandstructure(V0, wlen, a, m, L);
   p = {};
   p.update({'ratio_wlenx' : ratio , 'U' : band.Ut()[2] , 'tx_t' : band.t()[0]/band.t()[2] , 'tof_phase' : pyalps.dwa.tofPhase(time_of_flight=15.5, wavelength=[ratio*765,843,843], mass=86.99) })
   ps.append(p);
 
 params=[]
 for p in ps:
   for ratio in [0.8, 0.82, 0.84, 0.88, 0.9, 0.92, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.08, 1.10, 1.12, 1.14, 1.16, 1.18, 1.2]:
     params.append(
       { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
       , 'MODEL'           : 'boson Hubbard'                        
       , 'L'               : 160                             
       , 'Nmax'            : 20
       , 'ratio_wlenx'     : p['ratio_wlenx']                 
       , 't'               : 1.
       , 'tx_t'            : p['tx_t']                 
       , 'U'               : p['U']               
       , 'T'               : 1.                 
       , 'mu_homogeneous'  : ratio*4.0265                
       , 'mu'              : 'mu_homogeneous - (0.0073752*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.0036849*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.0039068155*(z-(L-1)/2.)*(z-(L-1)/2.))'
       , 'tof_phase'       : str(p['tof_phase']) 
       , 'SWEEPS'          : 100000             
       , 'SKIP'            : 100                 
       }
     )
 
 h5_infiles = pyalps.writeInputH5Files("parm.errorbudget.wlen",params);

Read from archive:

 h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.wlen');

=== Automated processing ===

 def parallel_automated_processing(taskfile):
   pyalps.dwa.recursiveRun( 
     "pyalps.runApplication('dwa', taskfile, T=2700)" ,  cmd_lang = 'python' ,
     follow_up_script = "pyalps.dwa.extract_worldlines(infile=pyalps.input2output(taskfile), outfile=taskfile)" ,
     break_if = "pyalps.dwa.thermalized(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
     write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Thermalizing')" ,
     loc = locals() ,
     batch_submit = True ,
     batch_cmd_prefix = 'bsub -J jobwlen -n1 -W60' ,
     batch_run_script = taskfile + '.thermalize.script' ,
     batch_next_run_script = taskfile + '.converge.script' ,
     batch_noRun = True
   );
   
   pyalps.dwa.recursiveRun( 
     "pyalps.runApplication('dwa', pyalps.input2output(taskfile), T=2700)" ,  cmd_lang = 'python' ,
     break_if = "pyalps.dwa.converged(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
     write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Converging')" ,
     loc = locals() ,
     batch_submit = True ,
     batch_cmd_prefix = 'bsub -J jobwlen -n1 -W60' ,
     batch_run_script = taskfile + '.converge.script' ,
     batch_noRun = True
   );
   
   pyalps.dwa.startRunScript(
     batch_run_script= taskfile + '.thermalize.script' ,
     batch_cmd_prefix = 'bsub -J jobwlen -n1 -W60' ,
     loc = locals()
   );

 for h5_infile in h5_infiles:
   parallel_automated_processing(h5_infile);

=== Unexpected norun ===

Sometimes, due to certain node failure, jobs may not even start. The following is the way to fix this:

 import pyalps;
 import pyalps.dwa;

 all_h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.wlen');
 h5_infiles = [];
 for h5_infile in all_h5_infiles:
    ar = pyalps.hdf5.h5ar(h5_infile);
    if not ar.is_group('/simulation/worldlines'):
       h5_infiles.append(h5_infile);

 for h5_infile in h5_infiles:
   parallel_automated_processing(h5_infile);

=== Unexpected worldlines configuration corruption ===

You may check worldlines configuration like the following:

 import pyalps;
 import pyalps.dwa;

Checking infiles:

 all_h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.wlen');
 h5_infiles = [];
 for h5_infile in all_h5_infiles:
    ar = pyalps.hdf5.h5ar(h5_infile);
    if ar.is_group('/simulation/worldlines'):
      wl = pyalps.dwa.extract_worldlines(h5_infile);
      valid = wl.is_valid(20);
      if not valid:
        h5_infiles.append(h5_infile);
        print('h5_infile : ' + h5_infile + ' , valid : ' + str(valid));

Checking outfiles:

 all_h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.wlen');
 h5_infiles = [];
 for h5_infile in all_h5_infiles:
    ar = pyalps.hdf5.h5ar(h5_infile);
    if ar.is_group('/simulation/worldlines'):
      wl = pyalps.dwa.extract_worldlines(pyalps.input2output(h5_infile));
      valid = wl.is_valid(20);
      if not valid:
        h5_infiles.append(h5_infile);
        print('h5_infile : ' + h5_infile + ' , valid : ' + str(valid));

=== Checking thermalization ===

 resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.wlen');
 for resultFile in resultFiles:
   ratio_wlenx = pyalps.getParameters(resultFile)[0]['ratio_wlenx'];
   mu=pyalps.getParameters(resultFile)[0]['mu_homogeneous'];
   N=pyalps.getMeasurements(resultFile,observable='Total Particle Number')[0]['mean']['value'];
   thermalized=pyalps.dwa.thermalized(resultFile, 'Total Particle Number', simplified=True);
   print('ratio_wlenx : ' + str(ratio_wlenx) + ' , mu : ' + str(mu) + ' , N : ' + str(N) + ' , thermalized : ' + str(thermalized));

=== Observing on your own ===

First the headers again:

 import matplotlib;
 matplotlib.use('Agg');

 import pyalps;
 import pyalps.plot;
 import pyalps.dwa;

We collect useful data:

 resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.wlen');
 data = pyalps.loadMeasurements(resultFiles, 'Total Particle Number');
 N = pyalps.collectXY(data, x='mu_homogeneous', y='Total Particle Number', foreach=['ratio_wlenx'])

then plot it:

 fig = matplotlib.pyplot.figure();
 pyalps.plot.plot(N);
 fig.savefig('parm.errorbudget.wlen.N.jpg');
 fig.show();

or send it to yourself:

 pyalps.sendmail( 'pingnang@phys.ethz.ch'    # email address of recipients 
                , attachment='parm.errorbudget.wlen.N.jpg'
                , subject='Picture: parm.errorbudget.wlen.N.jpg'
                );

=== Fine tuning the chemical potential ===

'''This section is under construction'''