09: Fix fix fix

Exercise sheet 9/latticescalar.py

@@ -5,8 +5,6 @@ import argparse
 import time
 import h5py
 
-starttime = time.asctime()
-
 rng = np.random.default_rng()  
 
 def potential_v(x,lamb):
@@ -68,13 +66,13 @@ def main():
 	
 	# perform sanity checks on the arguments
 	if args.w is None or args.w < 1:
-	    parser.error("Please specify a positive lattice size!")
+		parser.error("Please specify a positive lattice size!")
 	if args.k is None or args.k <= 0.0:
-	    parser.error("Please specify a positive kappa!")
+		parser.error("Please specify a positive kappa!")
 	if args.e < 10:
-	    parser.error("Need at least 10 equilibration sweeps to determine the average cluster size")
+		parser.error("Need at least 10 equilibration sweeps to determine the average cluster size")
 	if args.d < 1:
-	    parser.error("Delta should be >= 1.")
+		parser.error("Delta should be >= 1.")
 	
 	# fix parameters
 	lamb = args.l
@@ -84,17 +82,27 @@ def main():
 	delta = args.d
 	equil_sweeps = args.e
 	measure_sweeps = args.m
-	measurements = args.n
+	measurements = args.n or 0
 	if args.f is None:
-	    # construct a filename from the parameters plus a timestamp (to avoid overwriting)
+		 # construct a filename from the parameters plus a timestamp (to avoid overwriting)
-	    output_filename = "data_l{}_k{}_w{}_d{}_{}.hdf5".format(lamb,kappa,width,delta,time.strftime("%Y%m%d%H%M%S"))
+		 output_filename = "data_l{}_k{}_w{}_d{}_{}.hdf5".format(lamb,kappa,width,delta,time.strftime("%Y%m%d%H%M%S"))
 	else:
-	    output_filename = args.f
+		 output_filename = args.f
 
+	# Equilibration
+	phi_state = np.zeros((width,width,width,width))
+	run_scalar_MH(phi_state,lamb,kappa,delta,equil_sweeps * num_sites)
+
+	# Last things before we go-go
+	measurements_done = 0
+	starttime = time.time()
+	last_output_time = time.time()
+
+	# Prepare file
 	last_output_time = time.time()	
 	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)
+			dataset = f.create_dataset("magnetizations", (0,), maxshape=(None,), chunks=True)
 			# store some information as metadata for the data set
 			dataset.attrs["lamb"] = lamb
 			dataset.attrs["kappa"] = kappa
@@ -105,19 +113,31 @@ def main():
 			dataset.attrs["measure_sweeps"] = measure_sweeps
 			dataset.attrs["measurements"] = measurements
 			dataset.attrs["start time"] = starttime
-			dataset.attrs["stop time"] = time.asctime()
+
 
 	# Measure
 	# TODO: Does mean_magn need to be a list?
-	mean_magn = []
+	magnetizations = []
-	phi_state = np.zeros((width,width,width,width))
+	while True:
-	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)
+		magnetizations.append(np.mean(phi_state))
-	mean, err = batch_estimate(np.abs(magnetizations),lambda x:np.mean(x),10)
+		measurements_done += 1
-	mean_magn.append([mean,err])
+		if measurements == 0 or measurements_done == measurements or time.time() - last_output_time > args.o:
+			# time to output data again
+			with h5py.File(output_filename,'a') as f:
+				# enlarge the data set
+				f["magnetizations"].resize(measurements, axis=0)
+				# copy the data to the new space at the end of the dataset
+				f["magnetizations"][-len(magnetizations):] = magnetizations
+				f["magnetizations"].attrs["current_time"] = time.asctime()
+				magnetizations.clear()
+			if measurements_done == measurements:
+				break
+			else:
+				last_output_time = time.time()
+	# TODO: Store stoptime if any.
+	#dataset.attrs["stop time"] = time.asctime()
+	# TODO: Save if n is unset, m is set.
 
 if __name__ == "__main__":
 	main()