Add random number to the DeterminantRatio.template_params

Since random numbers are introduced, the template_params is not
determininstic and causes some headache with testing on different
systems (see Issue #59 and PR #63). Even still, random numbers are
introduced to easily obtain "unbaised" denominator that gives
determinant that are approximately 1 (at the very least, nonsingular)

wfns/wfn/quasiparticle/det_ratio.py

@@ -261,20 +261,30 @@ def template_params(self):
 
  Notes
  -----
- Depends on attribute `dimension`.
+ Depends on attribute `numerator_mask`.
 
  """
+ # NOTE: assume that same columns are selected for all matrices
+ # NOTE: all of the rows are assumed to be selected when the columns are selected.
+
  # get ground slater determinant
  ground_sd = slater.ground(*self.matrix_shape)
+ columns = self.get_columns(ground_sd, 0)
 
- # make matrices
+ # make matrix
  matrix = np.zeros(self.matrix_shape, dtype=self.dtype)
- # NOTE: assume that same columns are selected for all matrices
- # NOTE: all of the rows are assumed to be selected when the columns are selected.
- matrix[np.arange(self.matrix_shape[0]), self.get_columns(ground_sd, 0)] = 1
+ matrix[np.arange(self.matrix_shape[0]), columns] = 1
+
+ # make denominator
+ denominator = np.copy(matrix)
+ indices = [i for i in range(self.matrix_shape[1]) if i not in columns]
+ denominator[:, indices] = np.random.rand(self.matrix_shape[0], len(indices))
+ denominator /= np.linalg.norm(denominator, axis=0)
 
  # flatten and join together
- return np.hstack([matrix.flatten()] * self.num_matrices)
+ matrices = np.array([matrix] * self.num_matrices)
+ matrices[np.logical_not(self.numerator_mask)] = denominator
+ return matrices.flatten()
 
  # TODO: each matrix is assumed to have the same shape (nelec, nspin). It may need to be made
  # more flexible (let them have different shape or let each matrix have different shapes

wfns/wfn/quasiparticle/test/test_det_ratio.py

@@ -91,7 +91,17 @@ def test_template_params():
  test.nelec = 2
  test.nspin = 4
  test.numerator_mask = np.array([True, False, False])
- assert np.allclose(test.template_params, np.array([1, 0, 0, 0, 0, 0, 1, 0] * 3))
+
+ matrix = np.array([[1, 0, 0, 0], [0, 0, 1, 0]])
+ indices = np.array([True, False, True, False])
+ assert np.allclose(test.template_params[:8], matrix.flat)
+
+ template_params = test.template_params
+ template_matrix = template_params[8: 16].reshape(2, 4)
+ assert np.allclose(template_matrix[:, indices], matrix[:, indices])
+ assert np.allclose(np.linalg.norm(template_matrix[:, np.logical_not(indices)], axis=0), 1)
+
+ assert np.allclose(template_params[8: 16], template_params[16:])
 
 
 def test_assign_params():
@@ -102,7 +112,15 @@ def test_assign_params():
  test.nspin = 4
  test.numerator_mask = np.array([True, False, False])
  test.assign_params(params=None, add_noise=False)
- assert np.allclose(test.params, np.array([1, 0, 0, 0, 0, 0, 1, 0] * 3))
+
+ # check if template is correct (cannot check directly because template changes with every call)
+ matrix = np.array([[1, 0, 0, 0], [0, 0, 1, 0]])
+ indices = np.array([True, False, True, False])
+ assert np.allclose(test.params[:8], matrix.flat)
+ template_matrix = test.params[8: 16].reshape(2, 4)
+ assert np.allclose(template_matrix[:, indices], matrix[:, indices])
+ assert np.allclose(np.linalg.norm(template_matrix[:, np.logical_not(indices)], axis=0), 1)
+ assert np.allclose(test.params[16:], test.params[8: 16])
 
  params = np.array([1, 2, 3, 4, 5, 6, 7, 8] * 3, dtype=float)
  test.assign_params(params=params, add_noise=False)