Coverage for icet/tools/constraints.py: 100%

1import itertools 

2import numpy as np 

3from scipy.linalg import null_space 

4 

5 

6class Constraints: 

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

8 

9 Parameters 

10 ---------- 

11 n_params 

12 Number of parameters in model. 

13 

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:: 

18 

19 >>> import numpy as np 

20 >>> from icet.tools import Constraints 

21 >>> from trainstation import Optimizer 

22 

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

24 >>> n_params = 10 

25 >>> n_energies = 20 

26 >>> A = np.random.random((n_energies, n_params)) 

27 >>> y = np.random.random(n_energies) 

28 

29 >>> # Define constraints 

30 >>> c = Constraints(n_params=n_params) 

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

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

33 >>> c.add_constraint(M) 

34 

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

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

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

38 >>> opt.train() 

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

40 

41 """ 

42 

43 def __init__(self, n_params: int): 

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

45 self.constraint_vectors = np.eye(n_params) 

46 

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

48 """ Transform array to constrained parameter space. 

49 

50 Parameters 

51 ---------- 

52 A 

53 Array to be transformed. 

54 """ 

55 return A.dot(self.constraint_vectors) 

56 

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

58 """ Inverse transform array from constrained parameter space 

59 to unconstrained space. 

60 

61 Parameters 

62 ---------- 

63 A 

64 Array to be inverse transformed. 

65 """ 

66 return self.constraint_vectors.dot(A) 

67 

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

69 """ Add a constraint matrix and resolve for the constraint space. 

70 

71 Parameters 

72 ---------- 

73 M 

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

75 cluster vectors. 

76 """ 

77 M = np.array(M) 

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

79 self.constraint_vectors = null_space(self.M) 

80 

81 

82def get_mixing_energy_constraints(cluster_space) -> Constraints: 

83 """ 

84 A cluster expansion of the *mixing energy* should ideally predict zero energy 

85 for concentrations 0 and 1. This function constructs a :class:`Constraints` 

86 object that enforces that condition during training. 

87 

88 Parameters 

89 ---------- 

90 cluster_space 

91 Cluster space corresponding to cluster expansion for which constraints 

92 should be imposed. 

93 

94 Example 

95 ------- 

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

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

98 random target energies:: 

99 

100 >>> import numpy as np 

101 >>> from ase.build import bulk 

102 >>> from icet import ClusterSpace 

103 >>> from icet.tools import get_mixing_energy_constraints 

104 >>> from trainstation import Optimizer 

105 

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

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

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

109 >>> n_params = len(cs) 

110 >>> n_energies = 20 

111 >>> A = np.random.random((n_energies, n_params)) 

112 >>> y = np.random.random(n_energies) 

113 

114 >>> # Define constraints 

115 >>> c = get_mixing_energy_constraints(cs) 

116 

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

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

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

120 >>> opt.train() 

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

122 

123 Warning 

124 ------- 

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

126 fit quality of the others. Always expect that your cross-validation scores will 

127 increase somewhat when using this function. 

128 """ 

129 M = [] 

130 

131 prim = cluster_space.primitive_structure.copy() 

132 sublattices = cluster_space.get_sublattices(prim) 

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

134 

135 # Loop over all combinations of pure phases 

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

137 structure = prim.copy() 

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

139 for atom_index in subl.indices: 

140 structure[atom_index].symbol = symbol 

141 

142 # Add constraint for this pure phase 

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

144 

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

146 c.add_constraint(M) 

147 return c