 r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import numpy as np

2from scipy.linalg import null_space

3import itertools

6class Constraints:

7 """ Class for handling linear constraints with right hand side equal to zero.

9 Parameters

10 ----------

11 n_params

12 number of parameters in model

14 Example

15 -------

16 The following example demonstrates fitting of a cluster expansion under the

17 constraint that parameter 2 and parameter 4 should be equal::

19 >>> from icet.tools import Constraints

20 >>> from icet.fitting import Optimizer

21 >>> import numpy as np

23 >>> # Set up random sensing matrix and target "energies"

24 >>> n_params = 10

25 >>> A = np.random.random((10, n_params))

26 >>> y = np.random.random(10)

28 >>> # Define constraints

29 >>> c = Constraints(n_params=n_params)

30 >>> M = np.zeros((1, n_params))

31 >>> M[0, [2, 4]] = 1

34 >>> # Do the actual fit and finally extract parameters

35 >>> A_constrained = c.transform(A)

36 >>> opt = Optimizer((A_constrained, y), fit_method='ridge')

37 >>> opt.train()

38 >>> parameters = c.inverse_transform(opt.parameters)

40 """

42 def __init__(self, n_params: int):

43 self.M = np.empty((0, n_params))

44 self.constraint_vectors = np.eye(n_params)

46 def transform(self, A: np.ndarray) -> np.ndarray:

47 """ Transform array to constrained parameter space

49 Parameters

50 ----------

51 A

52 array to be transformed

53 """

54 return A.dot(self.constraint_vectors)

56 def inverse_transform(self, A: np.ndarray) -> np.ndarray:

57 """ Inverse transform array from constrained parameter space

58 to unconstrained space

60 Parameters

61 ----------

62 A

63 array to be inversed transformed

64 """

65 return self.constraint_vectors.dot(A)

67 def add_constraint(self, M: np.ndarray) -> None:

68 """ Add a constraint matrix and resolve for the constraint space

70 Parameters

71 ----------

72 M

73 Constraint matrix with each constraint as a row. Can (but need not be)

74 cluster vectors.

75 """

76 M = np.array(M)

77 self.M = np.vstack((self.M, M))

78 self.constraint_vectors = null_space(self.M)

81def get_mixing_energy_constraints(cluster_space) -> Constraints:

82 """

83 A cluster expansion of *mixing energy* should ideally predict zero energy

84 for concentration 0 and 1. This function constructs a :class:`Constraints`

85 object that enforces that condition during fitting.

87 Parameters

88 ----------

89 cluster_space : ClusterSpace

90 Cluster space corresponding to cluster expansion for which constraints

91 should be imposed

93 Example

94 -------

95 This example demonstrates how to constrain the mixing energy to zero

96 at the pure phases in a toy example with random cluster vectors and

97 random target energies::

99 >>> from icet.tools import get_mixing_energy_constraints

100 >>> from icet.fitting import Optimizer

101 >>> from icet import ClusterSpace

102 >>> from ase.build import bulk

103 >>> import numpy as np

105 >>> # Set up cluster space along with random sensing matrix and target "energies"

106 >>> prim = bulk('Au')

107 >>> cs = ClusterSpace(prim, cutoffs=[6.0, 5.0], chemical_symbols=['Au', 'Ag'])

108 >>> n_params = len(cs)

109 >>> A = np.random.random((10, len(cs)))

110 >>> y = np.random.random(10)

112 >>> # Define constraints

113 >>> c = get_mixing_energy_constraints(cs)

115 >>> # Do the actual fit and finally extract parameters

116 >>> A_constrained = c.transform(A)

117 >>> opt = Optimizer((A_constrained, y), fit_method='ridge')

118 >>> opt.train()

119 >>> parameters = c.inverse_transform(opt.parameters)

121 Warning

122 -------

123 Constraining the energy of one structure is always done at the expense of the

124 fit quality of the others. Always expect that your :term:`CV` scores will increase

125 somewhat when using this function.

126 """

127 M = []

129 prim = cluster_space.primitive_structure.copy()

130 sublattices = cluster_space.get_sublattices(prim)

131 chemical_symbols = [subl.chemical_symbols for subl in sublattices]

133 # Loop over all combinations of pure phases

134 for symbols in itertools.product(*chemical_symbols):

135 structure = prim.copy()

136 for subl, symbol in zip(sublattices, symbols):

137 for atom_index in subl.indices:

138 structure[atom_index].symbol = symbol

140 # Add constraint for this pure phase

141 M.append(cluster_space.get_cluster_vector(structure))

143 c = Constraints(n_params=len(cluster_space))