diff --git a/Exercise sheet 9/latticescalar.py b/Exercise sheet 9/latticescalar.py index 90da42a..6a5be34 100644 --- a/Exercise sheet 9/latticescalar.py +++ b/Exercise sheet 9/latticescalar.py @@ -48,38 +48,42 @@ def batch_estimate(data,observable,k): values = np.apply_along_axis(observable, 1, batches) return np.mean(values), np.std(values)/np.sqrt(k-1) -lamb = 1.5 -kappas = np.linspace(0.08,0.18,11) -width = 3 -num_sites = width**4 -delta = 1.5 # chosen to have ~ 50% acceptance -equil_sweeps = 10 -measure_sweeps = 2 -measurements = 20 +def main(): + lamb = 1.5 + kappas = np.linspace(0.08,0.18,11) + width = 3 + num_sites = width**4 + delta = 1.5 # chosen to have ~ 50% acceptance + equil_sweeps = 10 + measure_sweeps = 2 + measurements = 20 + + mean_magn = [] + for kappa in kappas: + phi_state = np.zeros((width,width,width,width)) + run_scalar_MH(phi_state,lamb,kappa,delta,equil_sweeps * num_sites) + magnetizations = np.empty(measurements) + for i in range(measurements): + run_scalar_MH(phi_state,lamb,kappa,delta,measure_sweeps * num_sites) + magnetizations[i] = np.mean(phi_state) + mean, err = batch_estimate(np.abs(magnetizations),lambda x:np.mean(x),10) + mean_magn.append([mean,err]) + + output_filename = 'preliminary_simulation.h5' + with h5py.File(output_filename,'a') as f: + if not "mean-magn" in f: + dataset = f.create_dataset("mean-magn", chunks=True, data=mean_magn) + # store some information as metadata for the data set + dataset.attrs["lamb"] = lamb + dataset.attrs["kappas"] = kappas + dataset.attrs["width"] = width + dataset.attrs["num_sites"] = num_sites + dataset.attrs["delta"] = delta + dataset.attrs["equil_sweeps"] = equil_sweeps + dataset.attrs["measure_sweeps"] = measure_sweeps + dataset.attrs["measurements"] = measurements + dataset.attrs["start time"] = starttime + dataset.attrs["stop time"] = time.asctime() -mean_magn = [] -for kappa in kappas: - phi_state = np.zeros((width,width,width,width)) - run_scalar_MH(phi_state,lamb,kappa,delta,equil_sweeps * num_sites) - magnetizations = np.empty(measurements) - for i in range(measurements): - run_scalar_MH(phi_state,lamb,kappa,delta,measure_sweeps * num_sites) - magnetizations[i] = np.mean(phi_state) - mean, err = batch_estimate(np.abs(magnetizations),lambda x:np.mean(x),10) - mean_magn.append([mean,err]) - -output_filename = 'preliminary_simulation.h5' -with h5py.File(output_filename,'a') as f: - if not "mean-magn" in f: - dataset = f.create_dataset("mean-magn", chunks=True, data=mean_magn) - # store some information as metadata for the data set - dataset.attrs["lamb"] = lamb - dataset.attrs["kappas"] = kappas - dataset.attrs["width"] = width - dataset.attrs["num_sites"] = num_sites - dataset.attrs["delta"] = delta - dataset.attrs["equil_sweeps"] = equil_sweeps - dataset.attrs["measure_sweeps"] = measure_sweeps - dataset.attrs["measurements"] = measurements - dataset.attrs["start time"] = starttime - dataset.attrs["stop time"] = time.asctime() +if __name__ == "__main__": + main()