Write configuration index in model[vec_types] (dba51b49) · Commits · Giacomo Mulas / NP_TMcode

src/scripts/guess_resources.py

+276 −11

Original line number	Diff line number	Diff line
		@@ -24,6 +24,7 @@
		# The script requires python3.

		import math
		import numpy as np
		import pdb
		import re
		from sys import argv
		@@ -58,6 +59,7 @@ def main():
		model = scan_model(config['edfb_file'], config['geom_file'])
		if (model is not None):
		print_model_summary(model)
		print_resource_summary(model)
		else:
		result = 1
		return result
		@@ -72,6 +74,76 @@ def iter_to_array(my_iter):
		array.append(item)
		return array

		## \brief Compute the particle's filling factor.
		#
		# In this application the filling factor is defined as the ratio of the
		# sum of the volumes of the spheres that compose the particle with respect
		# to the volume of a homogeneous ellipsoid that has the same rotational
		# properties as the particle. The result provided by this script is an
		# approximation that assumes the same specific weight for all monomers.
		# More accurate estimates (including the actual mass of the particle)
		# should be computed via the `inertia.py` script.
		#
		# \param model: `dict` Model description dictionary.
		# \return ff: `float` The particle's filling factor.
		def filling_factor(model):
		ff = 0.0
		N = model['nsph']
		Itot = np.zeros((3, 3))
		#vtot = get_total_volume(stypes, config)
		vtot = 0.0
		volumes = np.zeros((N))
		gf = 4.0 * math.pi / 3.0
		for i in range(model['nsph']):
		sph_type_index = model['vec_types'][i] - 1
		ros = model['ros'][sph_type_index]
		ros3 = ros * ros * ros
		volumes[i] = gf * ros3
		vtot += volumes[i]
		# for i ends here
		masses = volumes * 1000.0
		mtot = np.sum(masses, axis=0)

		positions = np.array([
		[model['vec_x'][i], model['vec_y'][i], model['vec_z'][i]] for i in range(N)
		])
		center_of_mass = np.sum(masses[:, np.newaxis] * positions, axis=0) / mtot
		for i in range(N):
		M = masses[i]
		typeID = model['vec_types'][i] - 1
		R = model['ros'][typeID]
		x = model['vec_x'][i] - center_of_mass[0]
		y = model['vec_y'][i] - center_of_mass[1]
		z = model['vec_z'][i] - center_of_mass[2]

		# 1. Get the inertia tensor of each sphere
		I_sph_val = (2/5) * M * R**2
		I_sph = np.eye(3) * I_sph_val

		# 2. Get the transfer tensor through the parallel axes theorem.
		# Diagonals: M * (sum of squared other coordinates)
		# Off-diagonals: -M * (product of coordinates)
		I_transfer = M * np.array([
		[y2 + z2, -x * y, -x * z],
		[-y * x, x2 + z2, -y * z],
		[-z * x, -z * y, x2 + y2]
		])
		Itot += (I_sph + I_transfer)
		# i loop ends here

		eigenval, eigenvec = np.linalg.eigh(Itot)
		print("INFO: ellipsoid directions:")
		for vi in range(3): print(" "*5, eigenvec[vi])
		gf = 5.0 / (2.0 * mtot)
		I1, I2, I3 = eigenval
		a = np.sqrt(max(0.0, gf * (I2 + I3 - I1)))
		b = np.sqrt(max(0.0, gf * (I1 + I3 - I2)))
		c = np.sqrt(max(0.0, gf * (I1 + I2 - I3)))
		# Volume of the ellipsoid that has the same moments of inertia as the particle
		vell = 4.0 * math.pi * a * b * c / 3.0
		ff = vtot / vell
		return ff

		## \brief Parse the command line arguments.
		#
		# The script behaviour can be modified through a set of optional arguments.
		@@ -139,20 +211,29 @@ def print_model_summary(model):
		avgX = 0.0
		avgY = 0.0
		avgZ = 0.0
		avgR = 0.0
		cmX = 0.0
		cmY = 0.0
		cmZ = 0.0
		Rmin = 0.0
		Rmax = 0.0
		Reqm = 0.0
		R3tot = 0.0
		Rcirc = 0.0
		square_farthest = 0.0
		rad_farthest = 0.0
		# breakpoint()
		for i in range(model['nsph']):
		sph_type_index = model['vec_types'][i] - 1
		ros = model['ros'][sph_type_index]
		ros3 = ros * ros * ros
		avgX += model['vec_x'][i]
		avgY += model['vec_y'][i]
		avgZ += model['vec_z'][i]
		R3tot += math.pow(ros, 3.0)
		avgR += ros
		cmX += (ros3 * model['vec_x'][i])
		cmY += (ros3 * model['vec_y'][i])
		cmZ += (ros3 * model['vec_z'][i])
		R3tot += ros3
		if (ros > Rmax):
		Rmax = ros
		if (Rmin == 0.0 or ros < Rmin):
		@@ -161,6 +242,11 @@ def print_model_summary(model):
		avgX /= model['nsph']
		avgY /= model['nsph']
		avgZ /= model['nsph']
		avgR /= model['nsph']
		avgR3 = avgR * avgR * avgR
		cmX /= (avgR3 * model['nsph'])
		cmY /= (avgR3 * model['nsph'])
		cmZ /= (avgR3 * model['nsph'])
		for i in range(model['nsph']):
		sph_type_index = model['vec_types'][i] - 1
		ros = model['ros'][sph_type_index]
		@@ -170,8 +256,12 @@ def print_model_summary(model):
		square_range = dX * dX + dY * dY + dZ * dZ + ros * ros
		if (square_range > square_farthest):
		square_farthest = square_range
		rad_farthest = ros
		Rcirc = math.sqrt(square_farthest) + rad_farthest
		Rcirc = math.sqrt(square_farthest)
		if (model['app'] == "SPHERE"):
		print("INFO: maximum expansion order LM = %d."%model['li'])
		else:
		print("INFO: maximum internal expansion order LI = %d."%model['li'])
		print("INFO: maximum external expansion order LE = %d."%model['le'])
		if (Rmax == Rmin):
		print("INFO: monomer radius is Rsph = %.5em"%Rmin)
		else:
		@@ -179,6 +269,155 @@ def print_model_summary(model):
		print("INFO: largest monomer radius is Rmax = %.5em"%Rmax)
		print("INFO: equivalent volume radius is Reqv = %.5em"%Reqm)
		print("INFO: minimum encircling radius is Rcirc = %.5em"%Rcirc)
		print(
		"INFO: geometric center at [{0:.5e}, {1:.5e}, {2:.5e}]".format(
		avgX, avgY, avgZ
		)
		)
		print(
		"INFO: center of filled volume at [{0:.5e}, {1:.5e}, {2:.5e}]".format(
		cmX, cmY, cmZ
		)
		)
		fill_factor = filling_factor(model)
		print("INFO: the particle's filling factor is %.5g"%fill_factor)

		## \brief Print a summary of resources needed per iteration.
		#
		# \param[in] model: `dict` A model description dictionary.
		def print_resource_summary(model):
		host_mem_gb = 0.0
		gpu_mem_gb = 0.0
		half_layers = 1 + int(model['max_layers'] / 2)
		nhspo = max(model['npnt'], model['npntts']) * 2 - 1
		if (model['app'] == 'SPHERE'):
		nlmmt = 2 * model['li'] * (model['li'] + 2)
		host_mem_gb += 16 * 2 * model['li'] # RMI + REI
		host_mem_gb += 16 * 4 * nlmmt # W
		host_mem_gb += 16 * 21 # FSAS + SAS + VINTS
		host_mem_gb += 8 * 7 # SSCS + SEXS + SABS + SQSCS + SQEXS + SQABS + GCSV
		host_mem_gb += 8 * (1 + model['max_layers']) # ROS + RC
		host_mem_gb += 8 # IOG + NSHL
		host_mem_gb += 16 * 2 * nhspo # RIS + DLRI
		host_mem_gb += 16 * half_layers # DC0
		host_mem_gb += 16 + 8 # VKT + VSZ
		host_mem_gb += 8 * model['max_layers'] # RCF
		host_mem_gb += 8 * (16 + 16) # CMULLR + CMUL
		host_mem_gb += 8 * 3 * 11 # UNIT VECTORS
		host_mem_gb += 8 + 8 # ARGI + ARGS
		host_mem_gb += 16 * 16 # VINT
		host_mem_gb += 8 * (model['li'] * 12 + 2) # ZPV + GAPS
		host_mem_gb += 8 * (2 * 4 + 4 * 4) # TQSE + TQSPE + TQSS + TQSPS
		host_mem_gb += 16 * half_layers # DC0M
		host_mem_gb += 8 * model['nxi'] # XIV
		host_mem_gb += 16 * 3 # ARG + S0 + TFSAS
		elif (model['app'] == 'CLUSTER'):
		nsph = model['nsph']
		ncou = nsph * (nsph - 1)
		litpo = 2 * model['li'] + 1
		litpos = litpo * litpo
		lmtpo = 2 * model['li'] * model['le'] + 1
		lmtpos = lmtpo * lmtpo
		lm = max(model['li'], model['le'])
		lmpo = lm + 1
		ndi = nsph * model['li'] * (model['li'] + 2)
		ndit = ndi + ndi
		nlem = model['le'] * (model['le'] + 2)
		nlemt = nlem + nlem
		nv3j = lm * lmpo * (2 * lm + 7) / 6
		host_mem_gb += 16 * ncou * litpo # VH
		host_mem_gb += 16 * nsph * lmtpo # VJ0
		host_mem_gb += 16 * ncou * litpos # VYHJ
		host_mem_gb += 16 * nsph * lmtpos # VYJ0
		host_mem_gb += 16 * 2 * nsph * model['li'] # RMI + REI
		host_mem_gb += 16 * nlemt * 4 # W
		host_mem_gb += 16 * nlemt * nlemt # AM0M
		host_mem_gb += 16 * nsph * 5 # FSAS + SAS
		host_mem_gb += 16 * nsph * 16 # VINTS
		host_mem_gb += 8 * nv3j # V3J0
		host_mem_gb += 8 * nsph * 11 # sphere CSs and Qs + coords and radii
		host_mem_gb += 8 * nsph * model['max_layers'] # RC
		host_mem_gb += 4 * lmpo * lm # IND3J
		host_mem_gb += 4 * 2 * nsph # IOG + NSHL
		host_mem_gb += 16 * (16 + 16 + 8) # VINT + VINTM + SCSCP + ECSCP + SCSCPM + ECSCPM
		host_mem_gb += 16 * nhspo * 2 # RIS + DLRI
		host_mem_gb += 16 * half_layers # DC0
		host_mem_gb += (16 + 8) * nsph # VKT + VSZ
		host_mem_gb += 16 * 5 # TFSAS + TSAS
		host_mem_gb += 8 * 4 # GCS + SCS + ECS + ACS
		host_mem_gb += 8 * lmtpo # RAC3J
		host_mem_gb += 16 * 2 * ndi * nlem # GIS + GLS
		host_mem_gb += 16 * ndit * nlemt # SAM
		host_mem_gb += 16 * ndit * ndit # AM
		gpu_mem_gb += 16 * ndit * ndit # AM
		host_mem_gb += 8 * nsph * model['max_layers'] # RCF
		host_mem_gb += 8 * (16 + 16 + 16 + 16) # CMULLR + CMUL + CEXTLR + CEXT
		host_mem_gb += 8 * 3 * 11 # UNIT VECTORS
		host_mem_gb += 8 + 8 # ARGI + ARGS
		host_mem_gb += 8 * (lm * 12 + nsph) # ZPV + GAPS
		host_mem_gb += 8 * (3 + 6 + 6) # GAP + GAPV + GAPM
		host_mem_gb += 16 * (6 + 6) # GAPP + GAPPM
		host_mem_gb += 8 * (2 * nsph + 4 * nsph) # TQSE + TQSPE + TQSS + TQSPS
		host_mem_gb += 8 * (2 * 4 + 4 * 4) # TQCE + TQCPE + TQCS + TQCPS
		host_mem_gb += 8 * (3 + 3) # TQEV + TQSV
		host_mem_gb += 16 * nsph * half_layers # DC0M
		host_mem_gb += 8 * model['nxi'] # XIV
		host_mem_gb += 16 * 4 # ARG + S0 + S0M + CCSAM
		elif (model['app'] == "INCLUSION"):
		nsph = model['nsph']
		ncou = nsph * (nsph - 1)
		litpo = 2 * model['li'] + 1
		litpos = litpo * litpo
		lmtpo = 2 * model['li'] * model['le'] + 1
		lmtpos = lmtpo * lmtpo
		lm = max(model['li'], model['le'])
		lmpo = lm + 1
		ndi = nsph * model['li'] * (model['li'] + 2)
		ndit = ndi + ndi
		nlem = model['le'] * (model['le'] + 2)
		nlemt = nlem + nlem
		nv3j = lm * lmpo * (2 * lm + 7) / 6
		host_mem_gb += 16 * ncou * litpo # VH
		host_mem_gb += 16 * nsph * lmtpo # VJ0
		host_mem_gb += 16 * ncou * litpos # VYHJ
		host_mem_gb += 16 * nsph * lmtpos # VYJ0
		host_mem_gb += 16 * 2 * nsph * model['li'] # RMI + REI
		host_mem_gb += 16 * nlemt * 4 # W
		host_mem_gb += 16 * (8 * model['le']) # RM0 + RE0 + RMW + REW + TM + TE + TM0 + TE0
		host_mem_gb += 16 * nlemt * nlemt # AM0M
		host_mem_gb += 8 * nv3j # V3J0
		host_mem_gb += 16 * 3 * 4 # FSAC + SAC + FSACM
		host_mem_gb += 16 * (16 + 16 + 8) # VINT + VINTM + SCSCP + ECSCP + SCSCPM + ECSCPM
		host_mem_gb += 8 * nsph * 4 # sphere coords and radii
		host_mem_gb += 8 * nsph * model['max_layers'] # RC
		host_mem_gb += 4 * lmpo * lm # IND3J
		host_mem_gb += 4 * 2 * nsph # IOG + NSHL
		host_mem_gb += 16 * nhspo * 2 # RIS + DLRI
		host_mem_gb += 16 * (half_layers + 1) # DC0
		host_mem_gb += (16 + 8) * nsph # VKT + VSZ
		host_mem_gb += 8 * lmtpo # RAC3J
		host_mem_gb += 16 * nlemt * ndi # AT
		host_mem_gb += 16 * ndit * ndit # AM
		gpu_mem_gb += 16 * ndit * ndit # AM
		host_mem_gb += 8 * nsph * model['max_layers'] # RCF
		host_mem_gb += 8 * (16 + 16 + 16 + 16) # CMULLR + CMUL + CEXTLR + CEXT
		host_mem_gb += 8 * 3 * 11 # UNIT VECTORS
		host_mem_gb += 8 + 8 # ARGI + ARGS
		host_mem_gb += 8 * (lm * 12 + nsph) # ZPV + GAPS
		host_mem_gb += 8 * (3 + 6 + 6) # GAP + GAPV + GAPM
		host_mem_gb += 16 * (6 + 6) # GAPP + GAPPM
		host_mem_gb += 8 * (2 * nsph + 4 * nsph) # TQSE + TQSPE + TQSS + TQSPS
		host_mem_gb += 8 * (2 * 4 + 4 * 4) # TQCE + TQCPE + TQCS + TQCPS
		host_mem_gb += 8 * (3 + 3) # TQEV + TQSV
		host_mem_gb += 16 * nsph * (half_layers + 1) # DC0M
		host_mem_gb += 8 * model['nxi'] # XIV
		host_mem_gb += 16 * 5 # ENT + ENTN + ARG + S0 + S0M

		host_mem_gb /= (1024.0 * 1024.0 * 1024.0)
		gpu_mem_gb /= (1024.0 * 1024.0 * 1024.0)
		print("INFO: 1 wavelength iteration uses %.5g Gb of host memory."%host_mem_gb)
		print("INFO: 1 wavelength iteration uses %.5g Gb of GPU memory"%gpu_mem_gb)
		print(" (or %.5g Gb, if iterative refinement is on)."%(3.0 * gpu_mem_gb))

		## \brief Scan the calculation model.
		#
		@@ -195,8 +434,13 @@ def scan_model(edfb_name, geom_name):
		'ies': 0,
		'li': 0,
		'le': 0,
		'npnt': 0,
		'npntts': 0,
		'nxi': 0,
		'configurations': 0,
		'max_layers': 0,
		'ros': [],
		'vec_nshl': [],
		'vec_types': [],
		'vec_x': [],
		'vec_y': [],
		@@ -204,6 +448,7 @@ def scan_model(edfb_name, geom_name):
		}
		# PARSING OF SCATTERER CONFIGURATION
		read_lines = 0
		max_layers = 0
		edfb_file = open(edfb_name, "r")
		file_line = edfb_file.readline()
		read_lines += 1
		@@ -243,6 +488,7 @@ def scan_model(edfb_name, geom_name):
		nxi = int(array_values[1].group())
		instpc = int(array_values[2].group())
		insn = int(array_values[3].group())
		model['nxi'] = nxi
		if (instpc != 0):
		file_line.readline()
		file_line.readline()
		@@ -255,6 +501,7 @@ def scan_model(edfb_name, geom_name):
		found_spheres = 0
		configurations = 0
		last_configuration = 0
		last_type = 0
		while (found_spheres < model['nsph']):
		file_line = edfb_file.readline()
		read_lines += 1
		@@ -269,10 +516,11 @@ def scan_model(edfb_name, geom_name):
		found_spheres += len(array_values)
		for ci in range(len(array_values)):
		type_id = int(array_values[ci].group())
		model['vec_types'].append(type_id)
		if (type_id > last_configuration):
		last_configuration = type_id
		if (type_id > last_type):
		last_type = type_id
		configurations += 1
		last_configuration += 1
		model['vec_types'].append(last_configuration)
		# end for ci block
		# end while(found_spheres) block
		model['configurations'] = configurations
		@@ -288,7 +536,9 @@ def scan_model(edfb_name, geom_name):
		edfb_file.close()
		return None
		nsh = int(array_values[0].group())
		model['vec_nshl'].append(nsh)
		if (ci == 0): nsh += 1
		if (nsh > max_layers): max_layers = nsh
		str_line = file_line[array_values[0].end():].replace('D', 'E').replace('d', 'E')
		iter_values = number_reg.finditer(str_line)
		array_values = iter_to_array(iter_values)
		@@ -317,6 +567,7 @@ def scan_model(edfb_name, geom_name):
		# end of if (ci == 0) block
		model['ros'].append(radius)
		# end of for ci block
		model['max_layers'] = max_layers
		else:
		# ies == 0
		if model['nsph'] == 1:
		@@ -344,6 +595,7 @@ def scan_model(edfb_name, geom_name):
		nxi = int(array_values[1].group())
		instpc = int(array_values[2].group())
		insn = int(array_values[3].group())
		model['nxi'] = nxi
		if (instpc != 0):
		file_line.readline()
		file_line.readline()
		@@ -356,7 +608,7 @@ def scan_model(edfb_name, geom_name):
		file_line = edfb_file.readline()
		file_line = edfb_file.readline()
		read_lines += 2
		model['vec_types'].append(1)
		model['vec_types'].append(0)
		iter_values = int_reg.finditer(file_line)
		array_values = iter_to_array(iter_values)
		if (len(array_values) < 1):
		@@ -365,6 +617,8 @@ def scan_model(edfb_name, geom_name):
		print(" at line %d in %s."%(read_lines, edfb_name))
		edfb_file.close()
		return None
		max_layers = int(array_values[0].group())
		model['vec_nshl'].append(max_layers)
		str_line = file_line[array_values[0].end():].replace('D', 'E').replace('d', 'E')
		iter_values = number_reg.finditer(str_line)
		array_values = iter_to_array(iter_values)
		@@ -376,6 +630,7 @@ def scan_model(edfb_name, geom_name):
		return None
		radius = float(array_values[0].group())
		model['ros'].append(radius)
		model['max_layers'] = max_layers
		else:
		model['app'] = "CLUSTER"
		file_line = edfb_file.readline()
		@@ -402,6 +657,7 @@ def scan_model(edfb_name, geom_name):
		nxi = int(array_values[1].group())
		instpc = int(array_values[2].group())
		insn = int(array_values[3].group())
		model['nxi'] = nxi
		if (instpc != 0):
		file_line.readline()
		file_line.readline()
		@@ -414,6 +670,7 @@ def scan_model(edfb_name, geom_name):
		found_spheres = 0
		configurations = 0
		last_configuration = 0
		last_type = 0
		while (found_spheres < model['nsph']):
		file_line = edfb_file.readline()
		read_lines += 1
		@@ -428,10 +685,11 @@ def scan_model(edfb_name, geom_name):
		found_spheres += len(array_values)
		for ci in range(len(array_values)):
		type_id = int(array_values[ci].group())
		model['vec_types'].append(type_id)
		if (type_id > last_configuration):
		last_configuration = type_id
		if (type_id > last_type):
		last_type = type_id
		configurations += 1
		last_configuration += 1
		model['vec_types'].append(last_configuration)
		# end for ci block
		# end while(found_spheres) block
		model['configurations'] = configurations
		@@ -447,6 +705,8 @@ def scan_model(edfb_name, geom_name):
		edfb_file.close()
		return None
		nsh = int(array_values[0].group())
		model['vec_nshl'].append(nsh)
		if (nsh > max_layers): max_layers = nsh
		str_line = file_line[array_values[0].end():].replace('D', 'E').replace('d', 'E')
		iter_values = number_reg.finditer(str_line)
		array_values = iter_to_array(iter_values)
		@@ -462,6 +722,7 @@ def scan_model(edfb_name, geom_name):
		read_lines += 1
		model['ros'].append(radius)
		# end of for ci block
		model['max_layers'] = max_layers
		# end if model['nsph'] block
		# end if model['ies'] block
		edfb_file.close()
		@@ -487,6 +748,8 @@ def scan_model(edfb_name, geom_name):
		return None
		model['li'] = int(array_values[1].group())
		model['le'] = int(array_values[2].group())
		model['npnt'] = int(array_values[5].group())
		model['npntts'] = int(array_values[6].group())
		for si in range(nsph):
		file_line = geom_file.readline()
		read_lines += 1
		@@ -519,6 +782,8 @@ def scan_model(edfb_name, geom_name):
		geom_file.close()
		return None
		model['li'] = int(array_values[1].group())
		model['npnt'] = int(array_values[3].group())
		model['npntts'] = int(array_values[4].group())
		model['vec_x'].append(0.0)
		model['vec_y'].append(0.0)
		model['vec_z'].append(0.0)