Coverage for icet/tools/convex_hull.py: 97%
96 statements
« prev ^ index » next coverage.py v7.10.1, created at 2025-09-14 04:08 +0000
1"""
2This module provides the ConvexHull class.
3"""
5import itertools
6from typing import Union
7import numpy as np
8from typing import Sized
9from scipy.interpolate import griddata
10from scipy.spatial import ConvexHull as ConvexHullSciPy
11from scipy.spatial import QhullError
14class ConvexHull:
15 """This class provides functionality for extracting the convex hull
16 of the (free) energy of mixing. It is based on the `convex hull
17 calculator in SciPy
18 <http://docs.scipy.org/doc/scipy-dev/reference/\
19generated/scipy.spatial.ConvexHull.html>`_.
21 Parameters
22 ----------
23 concentrations
24 Concentrations for each structure listed as ``[[c1, c2], [c1, c2],
25 ...]``; for binaries, in which case there is only one independent
26 concentration, the format ``[c1, c2, c3, ...]`` works as well.
27 energies
28 Energy (or energy of mixing) for each structure.
30 Attributes
31 ----------
32 concentrations : np.ndarray
33 Concentrations of the structures on the convex hull.
34 energies : np.ndarray
35 Energies of the structures on the convex hull.
36 dimensions : int
37 Number of independent concentrations needed to specify a point in
38 concentration space (1 for binaries, 2 for ternaries and so on).
39 structures : list[int]
40 Indices of structures that constitute the convex hull (indices are
41 defined by the order of their concentrations and energies are fed when
42 initializing the :class:`ConvexHull` object).
44 Examples
45 --------
46 A :class:`ConvexHull` object is easily initialized by providing lists of
47 concentrations and energies::
49 >>> data = {'concentration': [0, 0.2, 0.2, 0.3, 0.4, 0.5, 0.8, 1.0],
50 ... 'mixing_energy': [0.1, -0.2, -0.1, -0.2, 0.2, -0.4, -0.2, -0.1]}
51 >>> hull = ConvexHull(data['concentration'], data['mixing_energy'])
53 Now one can for example access the points along the convex hull directly::
55 >>> for c, e in zip(hull.concentrations, hull.energies):
56 ... print(c, e)
57 0.0 0.1
58 0.2 -0.2
59 0.5 -0.4
60 1.0 -0.1
62 or plot the convex hull along with the original data using e.g., matplotlib::
64 >>> import matplotlib.pyplot as plt
65 >>> plt.scatter(data['concentration'], data['mixing_energy'], color='darkred')
66 >>> plt.plot(hull.concentrations, hull.energies)
67 >>> plt.show(block=False)
69 It is also possible to extract structures at or close to the convex hull::
71 >>> low_energy_structures = hull.extract_low_energy_structures(
72 ... data['concentration'], data['mixing_energy'],
73 ... energy_tolerance=0.005)
75 A complete example can be found in the :ref:`basic tutorial
76 <tutorial_enumerate_structures>`.
77 """
79 def __init__(self,
80 concentrations: Union[list[float], list[list[float]]],
81 energies: list[float]) -> None:
82 assert len(concentrations) == len(energies)
83 # Prepare data in format suitable for SciPy-ConvexHull
84 concentrations = np.array(concentrations)
85 energies = np.array(energies)
86 points = np.column_stack((concentrations, energies))
87 self.dimensions = len(points[0]) - 1
89 # Construct convex hull
90 hull = ConvexHullSciPy(points)
92 # Collect convex hull points in handy arrays
93 concentrations = []
94 energies = []
95 for vertex in hull.vertices:
96 if self.dimensions == 1:
97 concentrations.append(points[vertex][0])
98 else:
99 concentrations.append(points[vertex][0:-1])
100 energies.append(points[vertex][-1])
101 concentrations = np.array(concentrations)
102 energies = np.array(energies)
104 structures = hull.vertices
105 # If there is just one independent concentration, we'd better sort
106 # according to it
107 if self.dimensions == 1:
108 ces = list(zip(*sorted(zip(concentrations, energies, structures))))
109 self.concentrations = np.array(ces[0])
110 self.energies = np.array(ces[1])
111 self.structures = np.array(ces[2])
112 else:
113 self.concentrations = concentrations
114 self.energies = energies
115 self.structures = structures
117 # Remove points that are above the "pure components plane"
118 self._remove_points_above_tie_plane()
120 def __str__(self) -> str:
121 width = 40
122 s = []
123 s += ['{s:=^{n}}'.format(s=' Convex Hull ', n=width)]
124 s += [' {:24} : {}'.format('dimensions', self.dimensions)]
125 s += [' {:24} : {}'.format('number of points', len(self.concentrations))]
126 s += [' {:24} : {}'.format('smallest concentration', np.min(self.concentrations))]
127 s += [' {:24} : {}'.format('largest concentration', np.max(self.concentrations))]
128 s += [''.center(width, '=')]
129 return '\n'.join(s)
131 def _repr_html_(self) -> str:
132 """ HTML representation. Used, e.g., in jupyter notebooks. """
133 s = ['<h4>Convex Hull</h4>']
134 s += ['<table border="1" class="dataframe">']
135 s += ['<tbody>']
136 s += [f'<tr><td style="text-align: left;">Dimensions</td><td>{self.dimensions}</td></tr>']
137 s += ['<tr><td style="text-align: left;">Number of points</td>'
138 f'<td>{len(self.concentrations)}</td></tr>']
139 s += ['<tr><td style="text-align: left;">Smallest concentration</td>'
140 f'<td>{np.min(self.concentrations)}</td></tr>']
141 s += ['<tr><td style="text-align: left;">Largest concentration</td>'
142 f'<td>{np.max(self.concentrations)}</td></tr>']
143 s += ['</tbody>']
144 s += ['</table>']
145 return ''.join(s)
147 def _remove_points_above_tie_plane(self, tol: float = 1e-6) -> None:
148 """
149 Remove all points on the convex hull that correspond to maximum rather
150 than minimum energy.
152 Parameters
153 ----------
154 tol
155 Tolerance for what energy constitutes a lower one.
156 """
158 # Identify the "complex concentration hull", i.e. the extremal
159 # concentrations. In the simplest case, these should simply be the
160 # pure components.
161 if self.dimensions == 1:
162 # Then the ConvexHullScipy function doesn't work, so we just pick
163 # the indices of the lowest and highest concentrations.
164 vertices = []
165 vertices.append(np.argmin(self.concentrations))
166 vertices.append(np.argmax(self.concentrations))
167 vertices = np.array(vertices)
168 else:
169 concentration_hull = ConvexHullSciPy(self.concentrations)
170 vertices = concentration_hull.vertices
172 # Remove all points of the convex energy hull that have an energy that
173 # is higher than what would be gotten with pure components at the same
174 # concentration. These points are mathematically on the convex hull,
175 # but in the physically uninteresting upper part, i.e. they maximize
176 # rather than minimize energy.
177 to_delete = []
178 for i, concentration in enumerate(self.concentrations):
179 # The points on the convex concentration hull should always be
180 # included, so skip them.
181 if i in vertices:
182 continue
184 # The energy obtained as a linear combination of concentrations on
185 # the convex hull is the "z coordinate" of the position on a
186 # (hyper)plane in the (number of independent concentrations +
187 # 1)-dimensional (N-D) space. This plane is spanned by N points.
188 # If there are more vertices on the convex hull, we need to loop
189 # over all combinations of N vertices.
190 for plane in itertools.combinations(vertices,
191 min(len(vertices),
192 self.dimensions + 1)):
193 # Calculate energy that would be gotten with pure components
194 # with ascribed concentration.
195 energy_pure = griddata(self.concentrations[np.array(plane)],
196 self.energies[np.array(plane)],
197 concentration,
198 method='linear')
200 # Prepare to delete if the energy was lowered. `griddata` gives
201 # NaN if the concentration is outside the triangle formed by
202 # the three vertices. The result of the below comparison is
203 # then False, which is what we want.
204 if energy_pure < self.energies[i] - tol:
205 to_delete.append(i)
206 break
208 # Finally remove all points
209 self.concentrations = np.delete(self.concentrations, to_delete, 0)
210 self.energies = np.delete(self.energies, to_delete, 0)
211 self.structures = list(np.delete(self.structures, to_delete, 0))
213 def get_energy_at_convex_hull(self, target_concentrations:
214 Union[list[float],
215 list[list[float]]]) -> np.ndarray:
216 """Returns the energy of the convex hull at specified concentrations.
217 If any concentration is outside the allowed range, NaN is
218 returned.
220 Parameters
221 ----------
222 target_concentrations
223 Concentrations at target points.
225 If there is one independent concentration, a list of
226 floats is sufficient. Otherwise, the concentrations ought
227 to be provided as a list of lists, such as ``[[0.1, 0.2],
228 [0.3, 0.1], ...]``.
229 """
230 if self.dimensions > 1 and isinstance(target_concentrations[0], Sized):
231 assert len(target_concentrations[0]) == self.dimensions
233 # Loop over all complexes of N+1 points to make sure that the lowest
234 # energy plane is used in the end. This is needed in two dimensions
235 # but in higher.
236 hull_candidate_energies = []
237 for plane in itertools.combinations(range(len(self.energies)),
238 min(len(self.energies),
239 self.dimensions + 1)):
240 try:
241 plane_energies = griddata(self.concentrations[list(plane)],
242 self.energies[list(plane)],
243 np.array(target_concentrations),
244 method='linear')
245 except QhullError:
246 # If the points lie on a line, the convex hull will fail, but
247 # we do not need to care about these "planes" anyway
248 continue
249 hull_candidate_energies.append(plane_energies)
251 # Pick out the lowest energies found
252 hull_energies = np.nanmin(hull_candidate_energies, axis=0)
253 return hull_energies
255 def extract_low_energy_structures(self, concentrations:
256 Union[list[float],
257 list[list[float]]],
258 energies: list[float],
259 energy_tolerance: float) -> list[int]:
260 """Returns the indices of energies that lie within a certain
261 tolerance of the convex hull.
263 Parameters
264 ----------
265 concentrations
266 Concentrations of candidate structures.
268 If there is one independent concentration, a list of
269 floats is sufficient. Otherwise, the concentrations must
270 be provided as a list of lists, such as ``[[0.1, 0.2],
271 [0.3, 0.1], ...]``.
272 energies
273 Energies of candidate structures.
274 energy_tolerance
275 Include structures with an energy that is at most this far
276 from the convex hull.
277 """
278 # Convert to numpy arrays, can be necessary if, for example,
279 # they are Pandas Series with "gaps"
280 concentrations = np.array(concentrations)
281 energies = np.array(energies)
283 n_points = len(concentrations)
284 if len(energies) != n_points: 284 ↛ 285line 284 didn't jump to line 285 because the condition on line 284 was never true
285 raise ValueError('concentrations and energies must have '
286 'the same length')
288 # Calculate energy at convex hull for specified concentrations
289 hull_energies = self.get_energy_at_convex_hull(concentrations)
291 # Extract those that are close enough
292 close_to_hull = [i for i in range(n_points)
293 if energies[i] <= hull_energies[i] + energy_tolerance]
295 return close_to_hull