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] = 1 

33 >>> M[0, 4] = -1 

34 >>> c.add_constraint(M) 

35 

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

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

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

39 >>> opt.train() 

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

41 

42 """ 

43 

44 def __init__(self, n_params: int): 

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

46 self.constraint_vectors = np.eye(n_params) 

47 

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

49 """ Transform array to constrained parameter space. 

50 

51 Parameters 

52 ---------- 

53 A 

54 Array to be transformed. 

55 """ 

56 return A.dot(self.constraint_vectors) 

57 

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

59 """ Inverse transform array from constrained parameter space 

60 to unconstrained space. 

61 

62 Parameters 

63 ---------- 

64 A 

65 Array to be inverse transformed. 

66 """ 

67 return self.constraint_vectors.dot(A) 

68 

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

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

71 

72 Parameters 

73 ---------- 

74 M 

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

76 cluster vectors. 

77 """ 

78 M = np.array(M) 

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

80 self.constraint_vectors = null_space(self.M) 

81 

82 

83def get_mixing_energy_constraints(cluster_space) -> Constraints: 

84 """ 

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

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

87 object that enforces that condition during training. 

88 

89 Parameters 

90 ---------- 

91 cluster_space 

92 Cluster space corresponding to cluster expansion for which constraints 

93 should be imposed. 

94 

95 Example 

96 ------- 

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

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

99 random target energies:: 

100 

101 >>> import numpy as np 

102 >>> from ase.build import bulk 

103 >>> from icet import ClusterSpace 

104 >>> from icet.tools import get_mixing_energy_constraints 

105 >>> from trainstation import Optimizer 

106 

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

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

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

110 >>> n_params = len(cs) 

111 >>> n_energies = 20 

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

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

114 

115 >>> # Define constraints 

116 >>> c = get_mixing_energy_constraints(cs) 

117 

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

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

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

121 >>> opt.train() 

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

123 

124 Warning 

125 ------- 

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

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

128 increase somewhat when using this function. 

129 """ 

130 M = [] 

131 

132 prim = cluster_space.primitive_structure.copy() 

133 sublattices = cluster_space.get_sublattices(prim) 

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

135 

136 # Loop over all combinations of pure phases 

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

138 structure = prim.copy() 

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

140 for atom_index in subl.indices: 

141 structure[atom_index].symbol = symbol 

142 

143 # Add constraint for this pure phase 

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

145 

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

147 c.add_constraint(M) 

148 return c