# Difference between revisions of "ALPS 2 Tutorial:MC-09 Directed Worm Algorithm:Intensive Example 2"

Jump to: navigation, search

# 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]);
```

### 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 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.1", params);
```

Run them:

```for h5_infile in h5_infiles:
parallel_automated_processing(h5_infile);
```

## Condensate fraction

### Turn on respective measurements

Gather the files:

```h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.V0.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 jobV0 -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 jobV0 -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 jobV0 -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 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');
h5_infiles = [];
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);
if not thermalized:
h5_infiles.append(resultFile);
print(resultFile + ' , ratio_wlenx : ' + str(ratio_wlenx) + ' , mu : ' + str(mu) + ' , N : ' + str(N) + ' , thermalized : ' + str(thermalized));
h5_infiles = [h5_infile.replace('out.h5','in.h5') for h5_infile in h5_infiles]
```

### 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

First the headers again:

```import numpy;
import pyalps;
import pyalps.dwa;
```

Easily done like this:

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

Prepare new input files:

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

Run them:

```for h5_infile in h5_infiles:
parallel_automated_processing(h5_infile);
```

## Condensate fraction

### Turn on respective measurements

Gather the files:

```h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.wlen.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 jobwlen -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 jobwlen -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 jobwlen -n1 -W60' ,
loc = locals()
);
```
```for h5_infile in h5_infiles:
parallel_automatic_processing(h5_infile);
```

### Extracting final data(plot)

This section is under construction.