epsi
/
msspec_python3
3
0
Fork 0
Code Issues 6 Pull Requests Releases Wiki Activity

Merge branch 'alter_devel' into devel

This commit is contained in:
Sylvain Tricot 2025-02-27 16:33:09 +01:00
parent 4aec48b4f6 5c5c57be2b
commit 8fe37b1509
4 changed files with 144 additions and 50 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
Makefile
View File

@ -38,10 +38,9 @@ light: venv
@$(INSIDE_VENV) pip install src/
nogui: VENV_PATH = ./_venv
nogui: venv
@$(INSIDE_VENV) pip install --no-cache-dir --upgrade -r src/pip.freeze
nogui: venv pybinding
@$(INSIDE_VENV) pip install -e src/
@+$(INSIDE_VENV) $(MAKE) -C src pybinding
_attrdict:
# Check if virtualenv python version > 3.3.0

2
src/msspec/iodata.py
View File

@ -17,7 +17,7 @@
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/iodata.py
# Last modified: Tue, 22 Oct 2024 12:39:54 +0200
# Last modified: Thu, 27 Feb 2025 16:33:09 +0100
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>

85
src/msspec/looper.py
View File

@ -17,8 +17,8 @@
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/looper.py
# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
# Last modified: Thu, 27 Feb 2025 16:33:09 +0100
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>
from collections import OrderedDict
from functools import partial
@ -92,9 +92,8 @@ class Sweep:
class SweepRange:
def __init__(self, *sweeps, passindex=False):
def __init__(self, *sweeps):
self.sweeps = sweeps
self.passindex = passindex
self.index = 0
# First check that sweeps that are linked to another on are all included
@ -158,17 +157,15 @@ class SweepRange:
for s in [sweep,] + children:
key, value = s[idx]
row[key] = value
if self.passindex:
row['sweep_index'] = i
row['sweep_index'] = i
return row
else:
raise StopIteration
@property
def columns(self):
cols = [sweep.key for sweep in self.sweeps]
if self.passindex:
cols.append('sweep_index')
cols = ['sweep_index']
cols += [sweep.key for sweep in self.sweeps]
return cols
@property
@ -202,31 +199,27 @@ class Looper:
logger.debug("Pipeline called with {}".format(x))
return self.pipeline(**x)
def run(self, *sweeps, ncpu=1, passindex=False):
def run(self, *sweeps, ncpu=1, **kwargs):
logger.info("Loop starts...")
# prepare the list of inputs
sr = SweepRange(*sweeps, passindex=passindex)
sr = SweepRange(*sweeps)
items = sr.items
data = []
t0 = time.time()
if ncpu == 1:
# serial processing...
logger.info("serial processing...")
t0 = time.time()
for i, values in enumerate(items):
values.update(kwargs)
result = self._wrapper(values)
data.append(result)
t1 = time.time()
dt = t1 - t0
logger.info("Processed {:d} sets of inputs in {:.3f} seconds".format(
len(sr), dt))
else:
# Parallel processing...
chunksize = 1 #int(nsets/ncpu)
[values.update(kwargs) for values in items]
logger.info(("Parallel processing over {:d} cpu (chunksize={:d})..."
"").format(ncpu, chunksize))
t0 = time.time()
@ -236,21 +229,23 @@ class Looper:
pool.close()
pool.join()
t1 = time.time()
dt = t1 - t0
logger.info(("Processed {:d} sets of inputs in {:.3f} seconds"
"").format(len(sr), dt))
t1 = time.time()
dt = t1 - t0
logger.info(("Processed {:d} sets of inputs in {:.3f} seconds"
"").format(len(sr), dt))
# Create the DataFrame
dfdata = []
columns = sr.columns + ['output',]
columns = sr.columns + list(kwargs.keys()) + ['output',]
for i in range(len(sr)):
row = list(items[i].values())
row.append(data[i])
dfdata.append(row)
df = pd.DataFrame(dfdata, columns=columns)
df = df.drop(columns=['sweep_index'])
self.data = df
@ -259,14 +254,14 @@ class Looper:
# of corresponding dict of parameters {'keyA': [val0,...valn],
# 'keyB': [val0,...valn], ...}
all_xy = []
for irow, row in df.iterrows():
all_xy.append(row.output[0])
all_xy.append(row.output[1])
parameters = df.to_dict()
parameters.pop('output')
# all_xy = []
# for irow, row in df.iterrows():
# all_xy.append(row.output[0])
# all_xy.append(row.output[1])
# parameters = df.to_dict()
# parameters.pop('output')
return all_xy, parameters
return self.data #all_xy, parameters
@ -276,17 +271,16 @@ class Looper:
if __name__ == "__main__":
import numpy as np
import time
import logging
logging.basicConfig(level=logging.DEBUG)
logger.setLevel(logging.DEBUG)
def bar(**kwargs):
return 0
def post_process(data):
x = data.x.unique()
y = data.y.unique()
i = kwargs.get('sweep_index')
return np.linspace(0,i,10)
theta = Sweep(key='theta', comments="The polar angle",
start=-70, stop=70, num=3,
@ -314,7 +308,16 @@ if __name__ == "__main__":
looper = Looper()
looper.pipeline = bar
data = looper.run(emitter, emitter_plane, uij, theta, levels, ncpu=4,
passindex=True)
other_kws = {'un':1, 'deux':2}
data = looper.run(emitter, emitter_plane, uij, theta, levels, ncpu=4, **other_kws)
# Print the dataframe
print(data)
#print(data[data.emitter_plane.eq(0)].theta.unique())
# Accessing the parameters and ouput values for a given sweep (e.g the last one)
print(looper.data.iloc[-1])
# Post-process the output values. For example here, the output is a 1D-array,
# make the sum of sweeps with 'Sr' emitter
X = np.array([ x for x in data[data.emitter == 'Sr'].output]).sum(axis=0)
print(X)

102
src/msspec/utils.py
View File

@ -19,8 +19,8 @@
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/utils.py
# Last modified: Thu, 06 Oct 2022 18:27:24 +0200
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr> 1665073644 +0200
# Last modified: Thu, 27 Feb 2025 16:33:09 +0100
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>
"""
@ -468,8 +468,12 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0,
# the symbol of your emitter
symbol = cluster[np.where(cluster.get_tags() == emitter_tag)[0][0]].symbol
assert (diameter != 0 or planes != 0), \
"At least one of diameter or planes parameter must be use."
if shape.lower() in ('spherical'):
assert (diameter != 0 or planes != 0), \
"At least one of diameter or planes parameter must be use."
elif shape.lower() in ('cylindrical'):
assert (diameter != 0 and planes != 0), \
"Diameter and planes parameters must be defined for cylindrical shape."
if diameter == 0:
# calculate the minimal diameter according to the number of planes
@ -479,6 +483,7 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0,
# number of repetition in each direction
rep = int(3*min_diameter/min(a, c))
#print("rep = ", rep)
# repeat the cluster
cluster = cluster.repeat((rep, rep, rep))
@ -542,7 +547,7 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0,
xplan, yplan = get_xypos(cluster, zplan)
radius = np.sqrt(xplan**2 + yplan**2 + zplan**2)
if diameter != 0:
if diameter != 0 and shape in ('spherical'):
assert (radius <= diameter/2), ("The number of planes is too high "
"compared to the diameter.")
radius = max(radius, diameter/2)
@ -575,3 +580,90 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0,
Atoms.translate(cluster, [0, 0, -ze]) # put the emitter in (0,0,0)
return cluster
def shape_cluster(primitive, emitter_tag=0, emitter_plane=0, diameter=0,
planes=0, shape='spherical'):
"""Creates and returns a cluster based on an Atoms object and some
parameters.
:param cluster: the Atoms object used to create the cluster
:type cluster: Atoms object
:param emitter_tag: the tag of your emitter
:type emitter_tag: integer
:param diameter: the diameter of your cluster in Angströms
:type diameter: float
:param planes: the number of planes of your cluster
:type planes: integer
:param emitter_plane: the plane where your emitter will be starting by 0
for the first plane
:type emitter_plane: integer
See :ref:`hemispherical_cluster_faq` for more informations.
"""
# We need the radius of the cluster and the number of planes
if shape.lower() in ('ispherical', 'cylindrical'):
assert (nplanes != 0 and diameter != 0), "nplanes and diameter cannot be zero for '{}' shape".format(shape)
elif shape.lower() in ('spherical'):
if diameter <= 0:
# find the diameter based on the number of planes
assert planes != 0, "planes should be > 0"
n = 3
natoms = 0
while True:
n += 2
cluster = primitive.copy()
# Repeat the primitive cell
cluster = cluster.repeat((n, n, n))
center_cluster(cluster)
# Find the emitter closest to the origin
all_tags = cluster.get_tags()
are_emitters = all_tags == emitter_tag
_ie = np.linalg.norm(cluster[are_emitters].positions, axis=1).argmin()
ie = np.nonzero(are_emitters)[0][_ie]
# Translate the cluster to this emitter position
cluster.translate(-cluster[ie].position)
# cut plane at surface and at bottom
all_z = np.unique(cluster.positions[:,2])
try:
zsurf = all_z[all_z >= 0][emitter_plane]
except IndexError:
# There are not enough planes above the emitter
zsurf = all_z.max()
try:
zbottom = all_z[all_z <= 0][::-1][planes - (emitter_plane+1)]
except IndexError:
# There are not enough planes below the emitter
zbottom = all_z.min()
cluster = cut_plane(cluster, z=[zbottom,zsurf])
# spherical shape
if shape.lower() in ('spherical'):
cluster = cut_sphere(cluster, radius=diameter/2, center=(0,0,zsurf))
if shape.lower() in ('ispherical'):
cluster = cut_sphere(cluster, radius=diameter/2, center=(0,0,0))
elif shape.lower() in ('cylindrical'):
cluster = cut_cylinder(cluster, radius=diameter/2)
else:
raise NameError("Unknown shape")
cluster.set_cell(primitive.cell)
if len(cluster) <= natoms:
break
else:
natoms = len(cluster)
return cluster