Update model maker to estimate system requirements and use exception handling+ù (96cc3945) · Commits · Giacomo Mulas / NP_TMcode

src/scripts/model_maker.py

+199 −173

Original line number	Diff line number	Diff line
		@@ -31,6 +31,7 @@ import pdb
		import random
		import yaml

		## \brief 3D software generation capability flag.
		allow_3d = True
		try:
		import pyvista as pv
		@@ -46,7 +47,7 @@ from sys import argv
		# `main()` is the function that handles the creation of the code configuration.
		# It returns an integer value as exit code, using 0 to signal successful execution.
		#
		# \returns result: `int` Number of detected error-level inconsistencies.
		# \returns result: `int` Exit code (0 = SUCCESS).
		def main():
		result = 0
		config = parse_arguments()
		@@ -68,11 +69,14 @@ def main():
		# \return result: `int` An exit code (0 if successful).
		def interpolate_constants(sconf):
		result = 0
		err_arg = ""
		try:
		for i in range(sconf['configurations']):
		for j in range(sconf['nshl'][i]):
		file_idx = sconf['dielec_id'][i][j]
		dielec_path = Path(sconf['dielec_path'], sconf['dielec_file'][int(file_idx) - 1])
		file_name = str(dielec_path)
		err_arg = str(dielec_path)
		file_name = err_arg
		dielec_file = open(file_name, 'r')
		wavelengths = []
		rpart = []
		@@ -125,6 +129,9 @@ def interpolate_constants(sconf):
		else:
		print("ERROR: file %s does not cover requested wavelengths!"%file_name)
		return 2
		except FileNotFoundError as ex:
		print("ERROR: file not found %s!"%err_arg)
		return 3
		return result

		## \brief Create tha calculation configuration structure from YAML input.
		@@ -261,12 +268,16 @@ def load_model(model_file):
		[0.0 for k in range(sconf['nxi'])] for j in range(sconf['configurations'])
		] for i in range(max_layers)
		]
		interpolate_constants(sconf)
		check = interpolate_constants(sconf)
		if (check != 0):
		return (None, None)
		else: # sconf[idfc] != 0 and scaling on wavelength
		print("ERROR: for wavelength scaling, optical constants must be tabulated!")
		return (None, None)
		elif (model['material_settings']['match_mode'] == "GRID"):
		match_grid(sconf)
		check = match_grid(sconf)
		if (check != 0):
		return(None, None)
		else:
		print("ERROR: %s is not a recognized match mode!"%(model['material_settings']['match_mode']))
		return (None, None)
		@@ -384,17 +395,25 @@ def load_model(model_file):
		rnd_engine = "COMPACT"
		if (rnd_engine == "COMPACT"):
		check = random_compact(sconf, gconf, rnd_seed, max_rad)
		if (check == 1):
		if (check == -1):
		print("ERROR: compact random generator works only when all sphere types have the same radius.")
		return (None, None)
		elif (check == -2):
		print("ERROR: sub-particle radius larger than particle radius.")
		return (None, None)
		elif (check == -3):
		print("ERROR: requested number of spheres cannot fit in allowed volume.")
		return (None, None)
		elif (rnd_engine == "LOOSE"):
		# random_aggregate() checks internally whether application is INCLUSION
		check = random_aggregate(sconf, gconf, rnd_seed, max_rad)
		else:
		print("ERROR: unrecognized random generator engine.")
		return (None, None)
		if (check != 0):
		print("WARNING: %d sphere(s) could not be placed."%check)
		if (check != sconf['nsph']):
		print("WARNING: placed only %d out of %d requested spheres."%(check, sconf['nsph']))
		sconf['nsph'] = check
		gconf['nsph'] = check
		else:
		if (len(model['geometry_settings']['x_coords']) != gconf['nsph']):
		print("ERROR: coordinate vectors do not match the number of spheres!")
		@@ -409,13 +428,16 @@ def load_model(model_file):
		write_obj(sconf, gconf, max_rad)
		try:
		max_gpu_ram = int(model['system_settings']['max_gpu_ram'])
		if (max_gpu_ram > 0):
		max_gpu_ram_bytes = max_gpu_ram * 1024 * 1024 * 1024
		matrix_dim = 2 * gconf['nsph'] * gconf['li'] * (gconf['li'] + 2)
		matrix_size_bytes = 16 * matrix_dim * matrix_dim
		matrix_size_Gb = float(matrix_size_bytes) / 1024.0 / 1024.0 / 1024.0
		print("INFO: estimated matrix size is {0:.3g} Gb.".format(matrix_size_Gb))
		if (max_gpu_ram > 0):
		max_gpu_ram_bytes = max_gpu_ram * 1024 * 1024 * 1024
		if (matrix_size_bytes < max_gpu_ram_bytes):
		max_gpu_processes = int(max_gpu_ram_bytes / matrix_size_bytes)
		print("INFO: system supports up to %d simultaneous processes on GPU."%max_gpu_processes)
		print("INFO: only %d GPU processes allowed, if using refinement."%(max_gpu_processes / 3))
		else:
		print("WARNING: estimated matrix size is larger than available GPU memory!")
		else:
		@@ -454,6 +476,8 @@ def match_grid(sconf):
		max_layers = 0
		nxi = 0
		sconf['vec_xi'] = []
		err_arg = ""
		try:
		for i in range(sconf['configurations']):
		layers = sconf['nshl'][i]
		if (sconf['application'] == "INCLUSION" and i == 0):
		@@ -461,7 +485,8 @@ def match_grid(sconf):
		for j in range(layers):
		file_idx = sconf['dielec_id'][i][j]
		dielec_path = Path(sconf['dielec_path'], sconf['dielec_file'][int(file_idx) - 1])
		file_name = str(dielec_path)
		err_arg = str(dielec_path)
		file_name = err_arg
		dielec_file = open(file_name, 'r')
		wavelengths = []
		rpart = []
		@@ -520,6 +545,9 @@ def match_grid(sconf):
		wi += 1
		sconf['rdc0'][j][i][dci] = ry
		sconf['idc0'][j][i][dci] = iy
		except FileNotFoundError as ex:
		print("ERROR: file not found %s!"%err_arg)
		return 3
		return result

		## \brief Parse the command line arguments.
		@@ -587,15 +615,19 @@ def random_aggregate(scatterer, geometry, seed, max_rad, max_attempts=100):
		vec_thetas = [0.0 for i in range(nsph)]
		vec_phis = [0.0 for i in range(nsph)]
		vec_rads = [0.0 for i in range(nsph)]
		vec_types = []
		n_types = scatterer['configurations']
		if (0 in scatterer['vec_types']):
		tincrement = 1 if scatterer['application'] != "INCLUSION" else 2
		for ti in range(nsph):
		itype = tincrement + int(n_types * random.random())
		scatterer['vec_types'][ti] = itype
		if (scatterer['application'] == "INCLUSION"):
		scatterer['vec_types'][0] = 1
		sph_type_index = scatterer['vec_types'][0] - 1
		vec_spheres = [{'itype': sph_type_index + 1, 'x': 0.0, 'y': 0.0, 'z': 0.0}]
		vec_rads[0] = scatterer['ros'][sph_type_index]
		vec_types.append(sph_type_index + 1)
		placed_spheres = 1
		attempts = 0
		for i in range(1, nsph):
		@@ -669,7 +701,9 @@ def random_aggregate(scatterer, geometry, seed, max_rad, max_attempts=100):
		})
		is_placed = True
		placed_spheres += 1
		vec_types.append(sph_type_index + 1)
		attempts = 0
		scatterer['vec_types'] = vec_types
		sph_index = 0
		for sphere in sorted(vec_spheres, key=lambda item: item['itype']):
		scatterer['vec_types'][sph_index] = sphere['itype']
		@@ -677,6 +711,7 @@ def random_aggregate(scatterer, geometry, seed, max_rad, max_attempts=100):
		geometry['vec_sph_y'][sph_index] = sphere['y']
		geometry['vec_sph_z'][sph_index] = sphere['z']
		sph_index += 1
		result = placed_spheres
		return result

		## \brief Generate a random compact cluster from YAML configuration options.
		@@ -699,45 +734,30 @@ def random_compact(scatterer, geometry, seed, max_rad):
		random.seed(seed)
		nsph = scatterer['nsph']
		n_types = scatterer['configurations']
		if (0 in scatterer['vec_types']):
		tincrement = 1 if scatterer['application'] != "INCLUSION" else 2
		for ti in range(nsph):
		itype = tincrement + int(n_types * random.random())
		scatterer['vec_types'][ti] = itype
		radius = scatterer['ros'][0]
		# Return an error code if types have different radii
		if (max(scatterer['ros']) != min(scatterer['ros'])):
		result = 1
		result = -1
		elif (radius > max_rad):
		# Requested spheres are larger than the maximum allowed volume.
		# End function with error code -2.
		result = -2
		else:
		radius = scatterer['ros'][0]
		x_centers = np.arange(-1.0 * max_rad + radius, max_rad, 2.0 * radius)
		y_centers = np.arange(-1.0 * max_rad + radius, max_rad, math.sqrt(3.0) * radius)
		z_centers = np.arange(-1.0 * max_rad + radius, max_rad, math.sqrt(3.0) * radius)
		x_offset = radius
		y_offset = radius
		x_layer_offset = radius
		y_layer_offset = radius / math.sqrt(3.0)
		x_centers = np.arange(-1.0 * max_rad + 2.0 * radius, max_rad, 2.0 * radius)
		x_size = len(x_centers)
		y_size = int(2.0 * max_rad / ((1.0 + math.sqrt(3.0) / 3.0) * radius))
		z_size = int(2.0 * max_rad / ((1.0 + 2.0 * math.sqrt(6.0) / 3.0) * radius))
		tmp_spheres = []
		n_cells = len(x_centers) * len(y_centers) * len(z_centers)
		n_cells = x_size * y_size * z_size
		print("INFO: the cubic space would contain %d spheres."%n_cells)
		n_max_spheres = int((max_rad / radius) * (max_rad / radius) * (max_rad / radius) * 0.74)
		print("INFO: the maximum radius allows for %d spheres."%n_max_spheres)
		for zi in range(len(z_centers)):
		if (x_layer_offset == 0.0):
		x_layer_offset = radius
		else:
		x_layer_offset = 0.0
		if (y_offset == 0.0):
		y_offset = radius
		else:
		y_offset = 0.0
		for yi in range(len(y_centers)):
		if (x_offset == 0.0):
		x_offset = radius
		else:
		x_offset = 0.0
		for xi in range(len(x_centers)):
		x = x_centers[xi] + x_offset + x_layer_offset
		y = y_centers[yi] + y_offset
		z = z_centers[zi]
		k = 0
		z = -max_rad + radius
		while (z < max_rad - radius):
		j = 0
		y = -max_rad + radius
		while (y < max_rad - radius):
		for i in range(len(x_centers)):
		x = (2 * (i + 1) + (j + k) % 2) * radius - max_rad
		extent = radius + math.sqrt(x * x + y * y + z * z)
		if (extent < max_rad):
		tmp_spheres.append({
		@@ -746,6 +766,11 @@ def random_compact(scatterer, geometry, seed, max_rad):
		'y': y,
		'z': z
		})
		#
		j += 1
		y = math.sqrt(3.0) * (j + (k % 2) / 3.0) * radius - max_rad + radius
		k += 1
		z = 2.0 / 3.0 * math.sqrt(6.0) * k * radius - max_rad + radius
		#tmp_spheres = [{'itype': 1, 'x': 0.0, 'y': 0.0, 'z': 0.0}]
		current_n = len(tmp_spheres)
		print("INFO: before erosion there are %d spheres in use."%current_n)
		@@ -777,20 +802,20 @@ def random_compact(scatterer, geometry, seed, max_rad):
		current_n -= 1
		vec_spheres = []
		sph_index = 0
		# Generate a vector of types if none is given
		if (0 in scatterer['vec_types']):
		tincrement = 1 if scatterer['application'] != "INCLUSION" else 2
		for ti in range(current_n):
		itype = tincrement + int(n_types * random.random())
		scatterer['vec_types'][ti] = itype
		if (scatterer['application'] == "INCLUSION"):
		scatterer['vec_types'][0] = 1
		for ti in range(len(tmp_spheres)):
		sphere = tmp_spheres[ti]
		if (sphere['x'] < max_rad):
		sphere['itype'] = scatterer['vec_types'][sph_index]
		sph_index += 1
		vec_spheres.append(sphere)
		#pl = pv.Plotter()
		#for si in range(len(vec_spheres)):
		# x = vec_spheres[si]['x'] / max_rad
		# y = vec_spheres[si]['y'] / max_rad
		# z = vec_spheres[si]['z'] / max_rad
		# mesh = pv.Sphere(radius / max_rad, (x, y, z))
		# pl.add_mesh(mesh)
		#pl.export_obj("scene.obj")
		sph_index = 0
		for sphere in sorted(vec_spheres, key=lambda item: item['itype']):
		scatterer['vec_types'][sph_index] = sphere['itype']
		@@ -798,7 +823,7 @@ def random_compact(scatterer, geometry, seed, max_rad):
		geometry['vec_sph_y'][sph_index] = sphere['y']
		geometry['vec_sph_z'][sph_index] = sphere['z']
		sph_index += 1
		return result
		return current_n

		## \brief Write the geometry configuration dictionary to legacy format.
		#
		@@ -952,9 +977,10 @@ def write_legacy_sconf(conf):

		## \brief Export the model to a set of OBJ files for 3D visualization.
		#
		# This function exports the model as a set of OBJ files (one for every
		# spherical unit, plus a single scene file) to allow for model visualization
		# with 3D software tools.
		# This function exports the model as a single OBJ file, containing the
		# information to visualize the particle with 3D software tools. The model
		# file is associated with a MTL material libray file, used to assign colors
		# to spheres of different type.
		#
		# \param scatterer: `dict` Scatterer configuration dictionary (gets modified)
		# \param geometry: `dict` Geometry configuration dictionary (gets modified)