diff --git a/ciderpress/dft/baselines.py b/ciderpress/dft/baselines.py index 6c9349c..95eef1c 100644 --- a/ciderpress/dft/baselines.py +++ b/ciderpress/dft/baselines.py @@ -1,6 +1,6 @@ import numpy as np -from ciderpress.density import LDA_FACTOR +from ciderpress.dft.settings import LDA_FACTOR def nsp_rho_basline(X0T, e=None, dedx=None): diff --git a/ciderpress/dft/settings.py b/ciderpress/dft/settings.py index e07bf0c..a4100d4 100644 --- a/ciderpress/dft/settings.py +++ b/ciderpress/dft/settings.py @@ -78,7 +78,7 @@ def get_cider_exponent( dadtau = C / rho if grad_mul != 0: grad_fac = grad_mul * 1.2 * (6 * np.pi**2) ** (2.0 / 3) / np.pi - grad_fac *= np.pi / 2 ** (2.0 / 3) * 0.125 + grad_fac *= np.pi / 2 ** (2.0 / 3) * 0.125 / CFC ascale += grad_fac * sigma / (rho * rho) dadsigma = grad_fac / (rho * rho) dadrho -= 2 * grad_fac * sigma / (rho * rho * rho) diff --git a/ciderpress/dft/tests/__init__.py b/ciderpress/dft/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/ciderpress/dft/tests/debug_normalizers.py b/ciderpress/dft/tests/debug_normalizers.py new file mode 100644 index 0000000..c3dc6c9 --- /dev/null +++ b/ciderpress/dft/tests/debug_normalizers.py @@ -0,0 +1,179 @@ +""" +This class provides tools for normalizing features to have more +desirable behavior under uniform scaling. +""" +from abc import ABC, abstractmethod + +import numpy as np + +from ciderpress.dft.settings import get_cider_exponent + + +def regularize_cider_exponent_(a, dadrho, dadsigma, dadtau, amin): + cond = a < amin + a[cond] = amin * np.exp(a[cond] / amin - 1) + grad_damp = a[cond] * np.exp(a[cond] / amin - 1) + dadrho[cond] *= grad_damp + dadsigma[cond] *= grad_damp + dadtau[cond] *= grad_damp + + +class FeatNormalizer(ABC): + @abstractmethod + def evaluate(self, rho, sigma, tau): + pass + + @abstractmethod + def evaluate_with_grad(self, rho, sigma, tau): + pass + + def apply_norm_fwd_(self, feat, rho, sigma, tau): + feat[:] *= self.evaluate(rho, sigma, tau) + + def apply_norm_bwd_(self, vfeat, feat, rho, sigma, tau, vrho, vsigma, vtau): + c, d1, d2, d3 = self.evaluate_with_grad(rho, sigma, tau) + vfeat[:] *= c + vrho[:] += d1 * vfeat * feat + vsigma[:] += d2 * vfeat * feat + vtau[:] += d3 * vfeat * feat + + @abstractmethod + def get_usp(self): + pass + + +class ConstantNormalizer(FeatNormalizer): + def __init__(self, const): + self.const = const + + def evaluate(self, rho, sigma, tau): + res = np.ones_like(rho) + res[:] *= self.const + return res + + def evaluate_with_grad(self, rho, sigma, tau): + res = self.evaluate(rho, sigma, tau) + dres = np.zeros_like(rho) + return res, dres, dres, dres + + def get_usp(self): + return 0 + + +class DensityPowerNormalizer(FeatNormalizer): + def __init__(self, const, power, cutoff=1e-10): + self.const = const + self.power = power + self.cutoff = cutoff + + def evaluate(self, rho, sigma, tau): + rho = np.maximum(self.cutoff, rho) + return self.const * rho**self.power + + def evaluate_with_grad(self, rho, sigma, tau): + cond = rho < self.cutoff + rho = np.maximum(self.cutoff, rho) + res = self.const * rho**self.power + dres = self.power * res / rho + dres[cond] = 0 + zeros = np.zeros_like(res) + return res, dres, zeros, zeros + + def get_usp(self): + return 3 * self.power + + +class CiderExponentNormalizer(FeatNormalizer): + def __init__(self, const1, const2, const3, power, cutoff=1e-10): + self.const1 = const1 + self.const2 = const2 + self.const3 = const3 + self.power = power + self.cutoff = cutoff + + def evaluate(self, rho, sigma, tau): + rho = np.maximum(self.cutoff, rho) + return ( + get_cider_exponent( + rho, + sigma, + tau, + a0=self.const1, + grad_mul=self.const2, + tau_mul=self.const3, + rhocut=self.cutoff, + nspin=1, + )[0] + ** self.power + ) + + def evaluate_with_grad(self, rho, sigma, tau): + cond = rho < self.cutoff + res0, d1, d2, d3 = get_cider_exponent( + rho, + sigma, + tau, + a0=self.const1, + grad_mul=self.const2, + tau_mul=self.const3, + rhocut=self.cutoff, + nspin=1, + ) + res = res0**self.power + tmp = self.power * res0 ** (self.power - 1) + d1 *= tmp + d2 *= tmp + d3 *= tmp + d1[cond] = 0 + d2[cond] = 0 + d3[cond] = 0 + return res, d1, d2, d3 + + def get_usp(self): + return 2 * self.power + + +class GeneralNormalizer(FeatNormalizer): + def __init__(self, const1, const2, const3, power1, power2, cutoff=1e-10): + self.norm1 = CiderExponentNormalizer(const1, const2, const3, power1, cutoff) + self.norm2 = DensityPowerNormalizer(1.0, power2, cutoff) + + def evaluate(self, rho, sigma, tau): + return self.norm1.evaluate(rho, sigma, tau) * self.norm2.evaluate( + rho, sigma, tau + ) + + def evaluate_with_grad(self, rho, sigma, tau): + n1, d1a, d1b, d1c = self.norm1.evaluate_with_grad(rho, sigma, tau) + n2, d2a, d2b, d2c = self.norm2.evaluate_with_grad(rho, sigma, tau) + return ( + n1 * n2, + n1 * d2a + n2 * d1a, + n1 * d2b + n2 * d1b, + n1 * d2c + n2 * d1c, + ) + + def get_usp(self): + return self.norm1.get_usp() + self.norm2.get_usp() + + +class FeatNormalizerList: + def __init__(self, normalizers): + """ + Args: + normalizers (list[FeatNormalizer]): list of feature normalizers + """ + self._normalizers = normalizers + + @property + def nfeat(self): + return len(self._normalizers) + + def get_usps(self): + return [n.get_usp() for n in self._normalizers] + + def __getitem__(self, item): + return self._normalizers[item] + + def __setitem__(self, key, value): + raise RuntimeError("Cannot set element of normalizer list") diff --git a/ciderpress/dft/tests/test_feat_normalizer.py b/ciderpress/dft/tests/test_feat_normalizer.py new file mode 100644 index 0000000..4a6c8fc --- /dev/null +++ b/ciderpress/dft/tests/test_feat_normalizer.py @@ -0,0 +1,154 @@ +import unittest + +import numpy as np +from numpy.testing import assert_almost_equal + +from ciderpress.dft.feat_normalizer import ( + CFC, + ConstantNormalizer, + DensityNormalizer, + FeatNormalizerList, + get_normalizer_from_exponent_params, +) +from ciderpress.dft.settings import get_alpha, get_s2 +from ciderpress.dft.tests import debug_normalizers as rho_normalizers + +DELTA = 1e-7 + + +def eval_fd(func, inps, i): + inps = [v.copy() for v in inps] + inps[i] += 0.5 * DELTA + fp = func(*inps) + inps[i] -= DELTA + fm = func(*inps) + return (fp - fm) / DELTA + + +class TestFeatNormalizer(unittest.TestCase): + + _x = None + norm_list: FeatNormalizerList = None + ref_norm_list: rho_normalizers.FeatNormalizerList = None + + @classmethod + def setUpClass(cls): + N = 50 + cls._x = np.linspace(0, 1, N) + np.random.seed(53) + directions = 0.25 * np.random.normal(size=(3, N)) + np.random.seed(53) + alpha = np.random.normal(size=N) ** 2 + rho = 1 + 0.5 * np.sin(np.pi * cls._x) + drho = directions * rho + tauw = np.einsum("xg,xg->g", drho, drho) / (8 * rho) + tau = tauw + alpha * CFC * rho ** (5.0 / 3) + rho_data = np.concatenate([[rho], drho, [tau]], axis=0) + assert rho_data.shape[0] == 5 + cls.rho = rho_data + sigma = np.einsum("xg,xg->g", drho, drho) + cls.norm_list = FeatNormalizerList( + [ + None, + None, + None, + ConstantNormalizer(3.14), + DensityNormalizer(1.2, 0.75), + get_normalizer_from_exponent_params(0.0, 1.5, 1.0, 0.03125), + get_normalizer_from_exponent_params( + 0.75, 1.5, 1.0 * 1.2, 0.03125 * 1.2 + ), + ] + ) + np.random.seed(72) + feat = np.random.normal(size=(cls.norm_list.nfeat - 3, N)) + dfeat = np.random.normal(size=(cls.norm_list.nfeat, N)) + cls.x = np.stack( + [ + rho, + get_s2(rho, sigma), + get_alpha(rho, sigma, tau), + ] + + [feat[i] for i in range(feat.shape[0])] + )[None, :, :] + cls.dx = dfeat + cls.ref_norm_list = rho_normalizers.FeatNormalizerList( + [ + rho_normalizers.ConstantNormalizer(3.14), + rho_normalizers.DensityPowerNormalizer(1.2, 0.75), + rho_normalizers.CiderExponentNormalizer(1.0, 0.0, 0.03125, 1.5), + rho_normalizers.GeneralNormalizer( + 1.0 * 1.2, 0.0, 0.03125 * 1.2, 1.5, 0.75 + ), + ] + ) + + def test_evaluate(self): + rho_data = self.rho.copy() + rho = rho_data[0] + sigma = np.einsum("xg,xg->g", rho_data[1:4], rho_data[1:4]) + tau = rho_data[4] + initial_vals = [] + + feat = self.norm_list.get_normalized_feature_vector(self.x) + dfeat = self.norm_list.get_derivative_of_normed_features(self.x[0], self.dx) + featp = self.norm_list.get_normalized_feature_vector(self.x + DELTA * self.dx) + featm = self.norm_list.get_normalized_feature_vector(self.x - DELTA * self.dx) + dfeat_ref = (featp - featm) / (2 * DELTA) + for i in range(self.norm_list.nfeat): + assert_almost_equal(dfeat[i], dfeat_ref[0, i]) + v = 2 * feat + vfeat = self.norm_list.get_derivative_wrt_unnormed_features(self.x, v) + for i in range(self.norm_list.nfeat): + x = self.x.copy() + x[:, i] += 0.5 * DELTA + featp = self.norm_list.get_normalized_feature_vector(x) + ep = (featp**2).sum(axis=1) + x[:, i] -= DELTA + featm = self.norm_list.get_normalized_feature_vector(x) + em = (featm**2).sum(axis=1) + assert_almost_equal(vfeat[:, i], (ep - em) / DELTA, 5) + + for i in range(self.ref_norm_list.nfeat): + res0 = self.ref_norm_list[i].evaluate(rho, sigma, tau) + res, drdn, drds, drdt = self.ref_norm_list[i].evaluate_with_grad( + rho, sigma, tau + ) + inps = [rho, sigma, tau] + drdn_ref = eval_fd(self.ref_norm_list[i].evaluate, inps, 0) + drds_ref = eval_fd(self.ref_norm_list[i].evaluate, inps, 1) + drdt_ref = eval_fd(self.ref_norm_list[i].evaluate, inps, 2) + assert_almost_equal(res, res0) + assert_almost_equal(feat[0, i + 3], res0 * self.x[0, 3 + i]) + assert_almost_equal(drdn, drdn_ref) + assert_almost_equal(drds, drds_ref) + assert_almost_equal(drdt, drdt_ref) + initial_vals.append((res, drdn, drds, drdt)) + + lambd = 1.7 + rho = lambd**3 * rho_data[0] + sigma = lambd**8 * np.einsum("xg,xg->g", rho_data[1:4], rho_data[1:4]) + tau = lambd**5 * rho_data[4] + for i in range(self.ref_norm_list.nfeat): + res, drdn, drds, drdt = self.ref_norm_list[i].evaluate_with_grad( + rho, sigma, tau + ) + res_ref, drdn_ref, drds_ref, drdt_ref = initial_vals[i] + usp = self.ref_norm_list[i].get_usp() + assert_almost_equal(res, lambd**usp * res_ref) + assert_almost_equal(drdn, lambd ** (usp - 3) * drdn_ref) + assert_almost_equal(drds, lambd ** (usp - 8) * drds_ref) + assert_almost_equal(drdt, lambd ** (usp - 5) * drdt_ref) + + x = self.x.copy() + x[:, 0] *= lambd**3 + scaled_feat = self.norm_list.get_normalized_feature_vector(x) + usps = np.array(self.norm_list.get_usps()) + usps[0] = 3 + scaled_ref = lambd ** (usps[None, :, None]) * feat + for i in range(scaled_ref.shape[1]): + assert_almost_equal(scaled_feat[:, i], scaled_ref[:, i]) + + +if __name__ == "__main__": + unittest.main() diff --git a/ciderpress/dft/tests/test_plans.py b/ciderpress/dft/tests/test_plans.py new file mode 100644 index 0000000..1ce84ed --- /dev/null +++ b/ciderpress/dft/tests/test_plans.py @@ -0,0 +1,377 @@ +import unittest + +import numpy as np +from numpy.testing import assert_almost_equal, assert_equal + +from ciderpress.density import CFC +from ciderpress.dft.plans import NLDFGaussianPlan, NLDFSplinePlan +from ciderpress.dft.settings import ( + NLDFSettingsVI, + NLDFSettingsVIJ, + NLDFSettingsVJ, + NLDFSettingsVK, +) + +DELTA = 1e-7 + + +def ovlp_reference(a, b, featid, extra_args=None): + if featid == "se": + return (np.pi / (a + b)) ** 1.5 + elif featid == "se_ar2": + return a * 1.5 * np.pi**1.5 / (a + b) ** 2.5 + elif featid == "se_a2r4": + return a * a * 15.0 / 4 * np.pi**1.5 / (a + b) ** 3.5 + elif featid == "se_erf_rinv": + assert extra_args is not None + c = extra_args[0] * a + return (np.pi / (a + b)) ** 1.5 / np.sqrt(1 + c / (a + b)) + else: + raise ValueError + + +class TestNLDFGaussianPlan(unittest.TestCase): + + _x = None + _a = None + _b = None + _c1 = None + _c2 = None + alpha0 = None + lambd = None + nalpha = None + PlanClass = NLDFGaussianPlan + + @classmethod + def setUpClass(cls): + N = 50 + cls._x = np.linspace(0, 1, N) + + # construct a reasonable set of densities + np.random.seed(53) + directions = 0.25 * np.random.normal(size=(3, N)) + np.random.seed(53) + alpha = np.random.normal(size=N) ** 2 + rho = 1 + 0.5 * np.sin(np.pi * cls._x) + drho = directions * rho + tauw = np.einsum("xg,xg->g", drho, drho) / (8 * rho) + tau = tauw + alpha * CFC * rho ** (5.0 / 3) + rho_data = np.concatenate([[rho], drho, [tau]], axis=0) + assert rho_data.shape[0] == 5 + cls.rho = rho_data + + # set up parameters for testing nonlocal parts + cls._a = 0.07 + cls._b = 5.3 + # test function f(x, r) = c1(x) * exp(-a*r^2) + c2(x) * exp(-b*r^2) + cls._c1 = 1.3 * cls._x**2 + # cls._c1 = 0 * cls._x ** 2 + cls._c2 = 0.2 - 0.5 * cls._x + # cls._c2 = 0.0 * cls._x + cls.alpha0 = 0.003 + cls.lambd = 1.6 + cls.nalpha = 30 + thetap = [2.0, 0.5, 0.0625] + feat_params = [ + [1.0, 0.25, 0.03125], + [1.0, 0.00, 0.03125], + [4.0, 0.00, 0.06250, 1.0], + [1.0, 0.25, 0.03125], + [1.0, 0.00, 0.03125], + ] + feat_specs = ["se", "se_a2r4", "se_erf_rinv", "se_ar2", "se"] + i0specs = ["se_lapl", "se_ap2r2", "se_apr2", "se", "se_ap"] + i1specs = ["se_rvec", "se_grad"] + i1dots = [(-1, 1), (-1, 0), (1, 0)] + cls.example_vi_settings = NLDFSettingsVI( + thetap, + "one", + i0specs, + i1specs, + i1dots, + ) + cls.example_vj_settings = NLDFSettingsVJ(thetap, "one", feat_specs, feat_params) + cls.example_vij_settings = NLDFSettingsVIJ( + thetap, + "one", + i0specs, + i1specs, + i1dots, + feat_specs, + feat_params, + ) + cls.example_vk_settings = NLDFSettingsVK( + thetap, + "one", + feat_params[:2], + rho_damp="exponential", + ) + + def _get_plan(self, settings, nspin, **kwargs): + return self.PlanClass( + settings, nspin, self.alpha0, self.lambd, self.nalpha, **kwargs + ) + + def test_eval_feat_exp(self): + plan = self._get_plan(self.example_vj_settings, 1) + plan2 = self._get_plan(self.example_vj_settings, 2) + sigma = np.einsum("xg,xg->g", self.rho[1:4], self.rho[1:4]) + results = [] + for i in range(5): + results.append(plan.eval_feat_exp(self.rho[0], sigma, self.rho[4], i=i)) + results.append(plan.eval_feat_exp(self.rho[0], sigma, self.rho[4], i=-1)) + for j in range(4): + assert_equal(results[0][j], results[3][j]) + assert_equal(results[1][j], results[4][j]) + assert_almost_equal(results[-1][j], 2 * results[0][j]) + assert_almost_equal(results[-1][j], 2 * results[3][j]) + for i in [0, 1, 2, 3, 4, -1]: + rho_tmp = self.rho[0].copy() + sigma_tmp = sigma.copy() + tau_tmp = self.rho[4].copy() + a, dadn, dadsigma, dadtau = plan.eval_feat_exp( + rho_tmp, sigma_tmp, tau_tmp, i=i + ) + a2, da2dn, da2dsigma, da2dtau = plan2.eval_feat_exp( + 0.5 * rho_tmp, 0.25 * sigma_tmp, 0.5 * tau_tmp, i=i + ) + assert_almost_equal(a2, a, 12) + assert_almost_equal(da2dn, 2 * dadn, 12) + assert_almost_equal(da2dsigma, 4 * dadsigma, 12) + assert_almost_equal(da2dtau, 2 * dadtau, 12) + inputs = [rho_tmp, sigma_tmp, tau_tmp] + derivs = [dadn, dadsigma, dadtau] + for input, deriv in zip(inputs, derivs): + input[:] += 0.5 * DELTA + ap = plan.eval_feat_exp(rho_tmp, sigma_tmp, tau_tmp, i=i)[0] + input[:] -= DELTA + am = plan.eval_feat_exp(rho_tmp, sigma_tmp, tau_tmp, i=i)[0] + input[:] += 0.5 * DELTA + assert_almost_equal(deriv, (ap - am) / DELTA, 7) + rho_tiny = 5e-11 * np.ones_like(rho_tmp) + sigma_tiny = sigma_tmp * rho_tiny**2 / rho_tmp**2 + tau_tiny = tau_tmp * rho_tiny / rho_tmp + a, dadn, dadsigma, dadtau = plan.eval_feat_exp( + rho_tiny, sigma_tiny, tau_tiny, i=i + ) + assert (dadn == 0).all() + assert (dadsigma == 0).all() + assert (dadtau == 0).all() + + def test_get_interpolation_arguments(self): + plan = self.PlanClass( + self.example_vij_settings, + 1, + self.alpha0, + self.lambd, + self.nalpha, + ) + rho = self.rho[0] + tau = self.rho[4] + sigma = np.einsum("xg,xg->g", self.rho[1:4], self.rho[1:4]) + for i in range(-1, 5): + rho_tmp = rho.copy() + sigma_tmp = sigma.copy() + tau_tmp = tau.copy() + a, dadn, dadsigma, dadtau = plan.get_interpolation_arguments( + rho_tmp, sigma_tmp, tau_tmp, i=i + ) + inputs = [rho_tmp, sigma_tmp, tau_tmp] + derivs = [dadn, dadsigma, dadtau] + for input, deriv in zip(inputs, derivs): + input[:] += 0.5 * DELTA + ap = plan.get_interpolation_arguments(rho_tmp, sigma_tmp, tau_tmp, i=i)[ + 0 + ] + input[:] -= DELTA + am = plan.get_interpolation_arguments(rho_tmp, sigma_tmp, tau_tmp, i=i)[ + 0 + ] + input[:] += 0.5 * DELTA + assert_almost_equal(deriv, (ap - am) / DELTA, 7) + rho_tiny = 5e-11 * np.ones_like(rho_tmp) + sigma_tiny = sigma_tmp * rho_tiny**2 / rho_tmp**2 + tau_tiny = tau_tmp * rho_tiny / rho_tmp + a, dadn, dadsigma, dadtau = plan.eval_feat_exp( + rho_tiny, sigma_tiny, tau_tiny, i=i + ) + assert (dadn == 0).all() + assert (dadsigma == 0).all() + assert (dadtau == 0).all() + try: + plan.get_interpolation_arguments(rho, sigma, tau, i=5) + raise AssertionError + except ValueError: + pass + + def check_get_interpolation_coefficients(self, **kwargs): + plan = self._get_plan(self.example_vij_settings, 1, **kwargs) + rho = self.rho[0] + tau = self.rho[4] + sigma = np.einsum("xg,xg->g", self.rho[1:4], self.rho[1:4]) + feat_parts = self._get_feat(plan.alphas[:, None], "se", None) + feat_parts *= plan.alpha_norms[:, None] + for i in range(-1, 5): + exp_g = plan.eval_feat_exp(rho, sigma, tau, i=i)[0] + arg_g = plan.get_interpolation_arguments(rho, sigma, tau, i=i)[0] + p, dp = plan.get_interpolation_coefficients(arg_g, i=i) + p2, dp2 = np.empty_like(p), np.empty_like(p) + plan.get_interpolation_coefficients(arg_g, i=i, vbuf=p2, dbuf=dp2) + assert_equal(p, p2) + assert_equal(dp, dp2) + arg_g[:] += 0.5 * DELTA + pp, _ = plan.get_interpolation_coefficients(arg_g, i=i) + arg_g[:] -= DELTA + pm, _ = plan.get_interpolation_coefficients(arg_g, i=i) + arg_g[:] += 0.5 * DELTA + assert_almost_equal(dp, (pp - pm) / DELTA, 4) + p = plan.get_transformed_interpolation_terms( + p, i=i, fwd=True, inplace=False + ) + if i == -1: + feat_id = "se" + else: + feat_id = plan.nldf_settings.feat_spec_list[i] + extra_args = plan._get_extra_args(i) if i > -1 else None + ref_feat = self._get_feat(exp_g, feat_id, extra_args) + test_feat = np.einsum("qg,{}->g".format(plan.coef_order), feat_parts, p) + assert_almost_equal(test_feat, ref_feat, 3) + + def test_get_interpolation_coefficients(self): + # self.check_get_interpolation_coefficients(alpha_formula='zexp') + for alpha_formula in ["etb", "zexp"]: + for coef_order in ["gq", "qg"]: + self.check_get_interpolation_coefficients( + coef_order=coef_order, + alpha_formula=alpha_formula, + ) + + def check_eval_rho_full(self, settings, coef_order="qg"): + plan = self._get_plan(settings, 1, coef_order=coef_order) + rho = self.rho[0] + drho = self.rho[1:4] + tau = self.rho[4] + sigma = np.einsum("xg,xg->g", drho, drho) + feat_parts = self._get_feat(plan.alphas[:, None], "se", None) + feat_parts *= plan.alpha_norms[:, None] + plan.get_transformed_interpolation_terms( + feat_parts.T if coef_order == "gq" else feat_parts, + i=0, + fwd=False, + inplace=True, + ) + f = plan.zero_coefs_full( + feat_parts.shape[1], + buffer_nrows=len(plan.nldf_settings.l1_feat_specs), + ) + vf = np.zeros_like(f) + if coef_order == "gq": + if isinstance(settings, (NLDFSettingsVJ, NLDFSettingsVIJ, NLDFSettingsVK)): + f[:, : plan.nalpha] = feat_parts.T + start = plan.nalpha + else: + start = 0 + if isinstance(settings, (NLDFSettingsVI, NLDFSettingsVIJ)): + f[:, start:] = np.linspace(0.5, 1.5, f.shape[0])[:, None] + else: + if isinstance(settings, (NLDFSettingsVJ, NLDFSettingsVIJ, NLDFSettingsVK)): + f[: plan.nalpha] = feat_parts + start = plan.nalpha + else: + start = 0 + if isinstance(settings, (NLDFSettingsVI, NLDFSettingsVIJ)): + f[start:] = np.linspace(0.5, 1.5, f.shape[1]) + feat, dfeat = plan.eval_rho_full(f, rho, drho, tau) + assert feat.shape == (plan.nldf_settings.nfeat, rho.size) + assert dfeat.shape == (plan.nldf_settings.num_feat_param_sets, rho.size) + for i in range(plan.nldf_settings.num_feat_param_sets): + exp_g = plan.eval_feat_exp(rho, sigma, tau, i=i)[0] + feat_id = plan.nldf_settings.feat_spec_list[i] + extra_args = plan._get_extra_args(i) if i > -1 else None + ref_feat = self._get_feat(exp_g, feat_id, extra_args) + assert_almost_equal(feat[i], ref_feat, 5) + + def get_e_and_vfeat(feat): + energy = (1 + (feat**2).sum(axis=0)) ** 0.5 + return (energy, feat / energy) + + energy, vfeat = get_e_and_vfeat(feat) + vrho = np.zeros_like(rho) + vdrho = np.zeros_like(drho) + vtau = np.zeros_like(tau) + plan.eval_vxc_full(vfeat, vrho, vdrho, vtau, dfeat, rho, drho, tau, vf=vf) + inputs = [rho, drho[0], drho[1], drho[2], tau] + derivs = [vrho, vdrho[0], vdrho[1], vdrho[2], vtau] + if coef_order == "gq": + inputs += [f.T[q] for q in range(f.shape[1])] + derivs += [vf.T[q] for q in range(f.shape[1])] + else: + inputs += [f[q] for q in range(f.shape[0])] + derivs += [vf[q] for q in range(f.shape[0])] + i = 0 + for target, deriv in zip(inputs, derivs): + i += 1 + target[:] += 0.5 * DELTA + feat_tmp, _ = plan.eval_rho_full(f, rho, drho, tau) + ep, _ = get_e_and_vfeat(feat_tmp) + target[:] -= DELTA + feat_tmp, _ = plan.eval_rho_full(f, rho, drho, tau) + em, _ = get_e_and_vfeat(feat_tmp) + target[:] += 0.5 * DELTA + assert_almost_equal(deriv, (ep - em) / DELTA, 6) + + def test_eval_rho_and_vxc_full(self): + self.check_eval_rho_full(self.example_vij_settings, coef_order="qg") + self.check_eval_rho_full(self.example_vi_settings, coef_order="qg") + self.check_eval_rho_full(self.example_vj_settings, coef_order="qg") + self.check_eval_rho_full(self.example_vk_settings, coef_order="qg") + self.check_eval_rho_full(self.example_vij_settings, coef_order="gq") + self.check_eval_rho_full(self.example_vi_settings, coef_order="gq") + self.check_eval_rho_full(self.example_vj_settings, coef_order="gq") + self.check_eval_rho_full(self.example_vk_settings, coef_order="gq") + + def check_version_k(self, coef_order="gq"): + plan = self._get_plan(self.example_vk_settings, 1, coef_order=coef_order) + rho = self.rho[0] + sigma = np.einsum("xg,xg->g", self.rho[1:4], self.rho[1:4]) + tau = self.rho[4] + exp_g = plan.eval_feat_exp(rho, sigma, tau, i=0)[0] + # make some other exponent for testing + ref_exp_g = exp_g / (0.4 + exp_g) + 0.5 * exp_g + # arg1_g = plan.get_interpolation_arguments(rho, sigma, tau, i=-1)[0] + arg2_g = plan.get_interpolation_arguments(rho, sigma, tau, i=0)[0] + p1, dp1 = plan.get_interpolation_coefficients(ref_exp_g, i=-1) + p2, dp2 = plan.get_interpolation_coefficients(arg2_g, i=0) + plan.get_transformed_interpolation_terms(p2, i=0, fwd=True, inplace=True) + refval = np.exp(-1.5 * ref_exp_g / exp_g) + if coef_order == "gq": + p2 *= plan.alpha_norms + else: + p2 *= plan.alpha_norms[:, None] + testval = np.einsum("{},{}->g".format(coef_order, coef_order), p1, p2) + assert_almost_equal(testval, refval, 3) + + def test_version_k(self): + self.check_version_k("gq") + self.check_version_k("qg") + + def _get_feat(self, arg_g, feat_id, extra_args): + res = self._c1 * ovlp_reference(arg_g, self._a, feat_id, extra_args=extra_args) + res += self._c2 * ovlp_reference(arg_g, self._b, feat_id, extra_args=extra_args) + return res + + +class TestNLDFSplinePlan(TestNLDFGaussianPlan): + + PlanClass = NLDFSplinePlan + + def _get_plan(self, settings, nspin, **kwargs): + spline_size = kwargs.pop("spline_size", 160) + kwargs["spline_size"] = spline_size + return self.PlanClass( + settings, nspin, self.alpha0, self.lambd, self.nalpha, **kwargs + ) + + +if __name__ == "__main__": + unittest.main() diff --git a/ciderpress/dft/tests/test_transform_data.py b/ciderpress/dft/tests/test_transform_data.py new file mode 100644 index 0000000..b52a898 --- /dev/null +++ b/ciderpress/dft/tests/test_transform_data.py @@ -0,0 +1,116 @@ +import unittest + +import numpy as np +from numpy.testing import assert_almost_equal + +from ciderpress.dft.transform_data import * + +TMP_TEST = "test_files/tmp" + +TEST_VEC = np.array( + [ + [0.0, 0.4, 8.0], + [0.0, 0.8, 5.0], + [0.0, 2.1, 3.0], + [0.0, 5.0, 20.0], + [-12.0, 0.0, 5.0], + [0.0, 5.0, 20.0], + ] +) + + +def get_grad_fd(vec, feat_list, delta=1e-8): + deriv = np.zeros(vec.shape) + for i in range(vec.shape[0]): + dvec = vec.copy() + dvec[i] += delta + deriv[i] += np.sum((feat_list(dvec.T).T - feat_list(vec.T).T) / delta, axis=0) + return deriv + + +def get_grad_a(vec, feat_list): + deriv = np.zeros(vec.shape) + feat = np.zeros((feat_list.nfeat, vec.shape[1])) + feat_list.fill_vals_(feat, vec) + fderiv = np.ones(feat.shape) + feat_list.fill_derivs_(deriv, fderiv, vec) + return deriv + + +class TestFeatureNormalizer(unittest.TestCase): + def test_lmap(self): + flist = FeatureList([LMap(0), LMap(3)]) + gradfd = get_grad_fd(TEST_VEC, flist) + grada = get_grad_a(TEST_VEC, flist) + assert_almost_equal(gradfd, grada) + + def test_tmap(self): + flist = FeatureList([TMap(0, 0, 1), TMap(1, 3, 2)]) + gradfd = get_grad_fd(TEST_VEC, flist) + grada = get_grad_a(TEST_VEC, flist) + assert_almost_equal(gradfd, grada) + + def test_umap(self): + flist = FeatureList([UMap(0, 0, 0.25), UMap(1, 3, 0.25)]) + gradfd = get_grad_fd(TEST_VEC, flist) + grada = get_grad_a(TEST_VEC, flist) + assert_almost_equal(gradfd, grada) + + def test_signed_umap(self): + flist = FeatureList([SignedUMap(0, 0.25), SignedUMap(3, 0.25)]) + gradfd = get_grad_fd(TEST_VEC, flist) + grada = get_grad_a(TEST_VEC, flist) + assert_almost_equal(gradfd, grada) + + def test_vmap(self): + flist = FeatureList( + [VMap(1, 0.25, scale=2.0, center=1.0), VMap(5, 0.5, scale=2.0, center=1.0)] + ) + gradfd = get_grad_fd(TEST_VEC, flist) + grada = get_grad_a(TEST_VEC, flist) + assert_almost_equal(gradfd, grada) + + def test_wmap(self): + flist = FeatureList([WMap(0, 1, 4, 0.5, 0.5)]) + gradfd = get_grad_fd(TEST_VEC, flist) + grada = get_grad_a(TEST_VEC, flist) + assert_almost_equal(gradfd, grada) + + def test_xmap(self): + flist = FeatureList([XMap(0, 1, 4, 0.5, 0.5)]) + gradfd = get_grad_fd(TEST_VEC, flist) + grada = get_grad_a(TEST_VEC, flist) + assert_almost_equal(gradfd, grada) + + def test_ymap(self): + flist = FeatureList([YMap(0, 1, 3, 4, 0.5, 0.3, 0.2)]) + gradfd = get_grad_fd(TEST_VEC, flist) + grada = get_grad_a(TEST_VEC, flist) + assert_almost_equal(gradfd, grada) + + def test_zmap(self): + flist = FeatureList( + [ZMap(1, 0.25, scale=2.0, center=1.0), ZMap(5, 0.5, scale=2.0, center=1.0)] + ) + gradfd = get_grad_fd(TEST_VEC, flist) + grada = get_grad_a(TEST_VEC, flist) + assert_almost_equal(gradfd, grada) + + def test_integration(self): + flist = FeatureList( + [ + UMap(0, 0.25), + UMap(1, 3, 0.25), + VMap(1, 0.25, scale=2.0, center=1.0), + WMap(0, 1, 4, 0.5, 0.5), + XMap(0, 1, 4, 0.5, 0.5), + YMap(0, 1, 3, 4, 0.5, 0.3, 0.2), + ] + ) + gradfd = get_grad_fd(TEST_VEC, flist) + grada = get_grad_a(TEST_VEC, flist) + assert_almost_equal(gradfd, grada, 6) + + +if __name__ == "__main__": + unittest.main() diff --git a/ciderpress/dft/tests/test_ueg.py b/ciderpress/dft/tests/test_ueg.py new file mode 100644 index 0000000..2e85d44 --- /dev/null +++ b/ciderpress/dft/tests/test_ueg.py @@ -0,0 +1,148 @@ +import unittest + +import numpy as np +from numpy.testing import assert_almost_equal + +from ciderpress.dft.debug_numint import get_get_exponent, get_nonlocal_features +from ciderpress.dft.settings import ( + ALLOWED_I_SPECS_L0, + ALLOWED_I_SPECS_L1, + ALLOWED_J_SPECS, + NLDFSettingsVI, + NLDFSettingsVIJ, + NLDFSettingsVJ, + NLDFSettingsVK, +) + +# constants for uniform grid +M = 100 +XMAX = 10.0 +N = 2 * M + 1 +H = (XMAX / M) ** 3 +TOL = 10 + + +class TestUEGVector(unittest.TestCase): + def _run_ueg_test(self, rho_const, gg_kwargs, vv_gg_kwargs): + tau_const = 0.3 * (3 * np.pi**2) ** (2.0 / 3) * rho_const ** (5.0 / 3) + + x = np.linspace(-XMAX, XMAX, N) + y = np.linspace(-XMAX, XMAX, N) + z = np.linspace(-XMAX, XMAX, N) + + xg, yg, zg = np.meshgrid(x, y, z) + coords = np.stack([xg.flatten(), yg.flatten(), zg.flatten()]).T + coords = np.ascontiguousarray(coords) + + G = coords.shape[0] + + vvrho = np.zeros((5, G)) + vvrho[0] += rho_const + vvrho[4] += tau_const + vvweight = H * np.ones(G) + rho = vvrho[:, :1].copy() + + get_exponent_r = get_get_exponent(gg_kwargs) + get_exponent_rp = get_get_exponent(vv_gg_kwargs) + theta_params = [vv_gg_kwargs["a0"], 0.0, vv_gg_kwargs["fac_mul"]] + feat_params = [] + feat_params_k = [] + for i in range(4): + feat_params.append([gg_kwargs["a0"], 0.0, gg_kwargs["fac_mul"]]) + feat_params_k.append( + [ + gg_kwargs["a0"] * 2 ** (i - 1), + 0.0, + gg_kwargs["fac_mul"] * 2 ** (i - 1), + ] + ) + feat_params[3].append(2.0) + + isettings = NLDFSettingsVI( + theta_params, + "one", + l0_feat_specs=ALLOWED_I_SPECS_L0, + l1_feat_specs=ALLOWED_I_SPECS_L1, + l1_feat_dots=[(-1, 0), (-1, 1), (0, 1)], + ) + jsettings = NLDFSettingsVJ( + theta_params, + "one", + feat_specs=ALLOWED_J_SPECS, + feat_params=feat_params, + ) + ijsettings = NLDFSettingsVIJ( + theta_params, + "one", + l0_feat_specs_i=ALLOWED_I_SPECS_L0, + l1_feat_specs_i=ALLOWED_I_SPECS_L1, + l1_feat_dots_i=[(-1, 0), (-1, 1), (0, 1)], + feat_specs_j=ALLOWED_J_SPECS, + feat_params_j=feat_params, + ) + ksettings = NLDFSettingsVK( + theta_params, "one", feat_params=feat_params_k, rho_damp="exponential" + ) + + zero = np.zeros((1, 3)) + ifeat = get_nonlocal_features( + rho, + zero, + vvrho, + vvweight, + coords, + get_exponent_r, + get_exponent_rp, + version="i", + )[0] + ifeat0 = ifeat[:6] + ifeat1 = np.zeros(3) + ifeat = np.append(ifeat0, ifeat1) + jfeat = get_nonlocal_features( + rho, + zero, + vvrho, + vvweight, + coords, + get_exponent_r, + get_exponent_rp, + version="j", + )[0] + kfeat = get_nonlocal_features( + rho, + zero, + vvrho, + vvweight, + coords, + get_exponent_r, + get_exponent_rp, + version="k", + )[0] + ijfeat = np.append(jfeat, ifeat) + + ifeat_ana = isettings.ueg_vector(rho_const) + jfeat_ana = jsettings.ueg_vector(rho_const) + ijfeat_ana = ijsettings.ueg_vector(rho_const) + kfeat_ana = ksettings.ueg_vector(rho_const) + assert_almost_equal(ijfeat_ana, np.append(jfeat_ana, ifeat_ana), TOL) + + assert_almost_equal(ifeat_ana, ifeat, TOL) + assert_almost_equal(jfeat_ana, jfeat, TOL) + assert_almost_equal(ijfeat_ana, ijfeat, TOL) + assert_almost_equal(kfeat_ana, kfeat, TOL) + + def test_ueg(self): + num = 1.0 + vv_gg_kwargs = {"a0": 1.0 * num, "fac_mul": 0.03125 * num, "amin": 0.0625 * num} + num *= 1 + gg_kwargs = {"a0": 1.0 * num, "fac_mul": 0.03125 * num, "amin": 0.0625 * num} + self._run_ueg_test(1.0, gg_kwargs, vv_gg_kwargs) + num = 2.3 + vv_gg_kwargs = {"a0": 1.0 * num, "fac_mul": 0.03125 * num, "amin": 0.0625 * num} + num *= 0.7 + gg_kwargs = {"a0": 1.0 * num, "fac_mul": 0.03125 * num, "amin": 0.0625 * num} + self._run_ueg_test(0.8, gg_kwargs, vv_gg_kwargs) + + +if __name__ == "__main__": + unittest.main()