From 69f31bc7c90a1182392ccb22b549350c36cbc1f5 Mon Sep 17 00:00:00 2001 From: Mark Dewing Date: Fri, 1 Feb 2019 14:27:29 -0600 Subject: [PATCH 1/2] Add unit tests for cubic splines Solve for the coefficients using a straightforward (but inefficient) derivation from the defining equations. The functions in eqn_manip.py perform the equation manipulations one might use when working by hand to put equations in a particular form - eg, moving all of one type of term to one side, in order to set up a matrix equation. Test the coefficient solver and the evaluation class (OneDimCubicSpline). --- src/Numerics/tests/CMakeLists.txt | 4 +- src/Numerics/tests/eqn_manip.py | 47 +++ src/Numerics/tests/gen_cubic_spline.py | 337 ++++++++++++++++++ src/Numerics/tests/test_nr_spline.cpp | 39 -- .../tests/test_one_dim_cubic_spline.cpp | 178 +++++++++ 5 files changed, 564 insertions(+), 41 deletions(-) create mode 100644 src/Numerics/tests/eqn_manip.py create mode 100644 src/Numerics/tests/gen_cubic_spline.py delete mode 100644 src/Numerics/tests/test_nr_spline.cpp create mode 100644 src/Numerics/tests/test_one_dim_cubic_spline.cpp diff --git a/src/Numerics/tests/CMakeLists.txt b/src/Numerics/tests/CMakeLists.txt index e6d8dcb54..f2c433f9b 100644 --- a/src/Numerics/tests/CMakeLists.txt +++ b/src/Numerics/tests/CMakeLists.txt @@ -19,10 +19,10 @@ SET(UTEST_NAME deterministic-unit_test_${SRC_DIR}) -SET(UTEST_SRCS test_grid_functor.cpp test_nr_spline.cpp test_stdlib.cpp test_bessel.cpp +SET(UTEST_SRCS test_grid_functor.cpp test_stdlib.cpp test_bessel.cpp test_ylm.cpp test_e2iphi.cpp test_aligned_allocator.cpp test_gaussian_basis.cpp test_cartesian_tensor.cpp test_soa_cartesian_tensor.cpp - test_transform.cpp test_min_oned.cpp) + test_transform.cpp test_min_oned.cpp test_one_dim_cubic_spline.cpp) # Run gen_gto.py to create these files. They may take a long time to compile. #SET(UTEST_SRCS ${UTEST_SRCS} test_full_cartesian_tensor.cpp test_full_soa_cartesian_tensor.cpp) diff --git a/src/Numerics/tests/eqn_manip.py b/src/Numerics/tests/eqn_manip.py new file mode 100644 index 000000000..e8e3bc3d6 --- /dev/null +++ b/src/Numerics/tests/eqn_manip.py @@ -0,0 +1,47 @@ +from __future__ import print_function + +from sympy import Eq + +# Move symbols in sym_list from left hand side of equation to right hand side +def move_terms(eqn, sym_list): + new_lhs = eqn.lhs + new_rhs = eqn.rhs + for sym in sym_list: + c = eqn.lhs.coeff(sym) + new_lhs = new_lhs - c*sym + new_rhs = new_rhs - c*sym + return Eq(new_lhs, new_rhs) + +# Move symbols in sym_list from right hand side of equation to left hand side +def move_terms_left(eqn, sym_list): + new_lhs = eqn.lhs + new_rhs = eqn.rhs + for sym in sym_list: + c = eqn.rhs.coeff(sym) + new_lhs = new_lhs - c*sym + new_rhs = new_rhs - c*sym + return Eq(new_lhs, new_rhs) + +def divide_terms(eqn, sym_list_left, sym_list_right): + #print 'start',eqn + eqn1 = move_terms(eqn, sym_list_right) + #print 'middle ',eqn1 + eqn2 = move_terms_left(eqn1, sym_list_left) + return eqn2 + +# Multiply equation by term +def mult_eqn(eqn, e): + return Eq(eqn.lhs*e, eqn.rhs*e) + +# for all values other than the target, set to zero +def get_coeff_for(expr, sym, symlist): + # expression + # symbol to get coefficient for + # symlist - total list of symbols + subslist = {} + subslist[sym] = 1 + for s in symlist: + if s != sym: + subslist[s] = 0 + coeff = expr.subs(subslist) + return coeff diff --git a/src/Numerics/tests/gen_cubic_spline.py b/src/Numerics/tests/gen_cubic_spline.py new file mode 100644 index 000000000..5516dd30f --- /dev/null +++ b/src/Numerics/tests/gen_cubic_spline.py @@ -0,0 +1,337 @@ + +from __future__ import print_function +from sympy import * +from eqn_manip import * + +# See for more explanation, see +# https://github.com/QMCPACK/qmc_algorithms/blob/master/Wavefunctions/Cubic%20Splines%20Basic.ipynb + +# Distance to knot, for non-uniform knots +t = IndexedBase('t') + +y = IndexedBase('y') +x = IndexedBase('x') + + +# Create the equations that define cubic splines and solve by brute force +def create_solution_for_val(nknots, naturalBC=(True, True), firstDeriv=(0.0, 0.0)): + + n = Symbol('n', integer=True) + i = Symbol('i', integer=True) + + a,b,c,d = [IndexedBase(s) for s in 'a b c d'.split()] + # Non-uniform knots + si = a[i] + b[i]*t[i] + c[i]*t[i]*t[i] + d[i]*t[i]**3 + + # Value at knots (t=0) + eq1 = Eq(si.subs(t[i],0), y[i]) + + # Value at knots (t=1) + eq2 = Eq(si.subs(t[i],x[i+1]-x[i]), y[i+1]) + + # Continuity of first derivatives at knots + dsi = diff(si, t[i]) + eq3 = Eq(dsi.subs(t[i],x[i+1]-x[i]), dsi.subs(i, i+1).subs(t[i+1],0)) + + # Continuity of second derivatives at knots + d2si = diff(si, t[i], 2) + eq4 = Eq(d2si.subs(t[i],x[i+1]-x[i]).subs(t[0],x[0]), d2si.subs(i, i+1).subs(t[i+1],0)) + + + sym_lhs = [a[i],b[i],c[i],d[i],a[i+1],b[i+1],c[i+1],d[i+1]] + sym_rhs = [y[i]] + eq3b = divide_terms(eq3, sym_lhs, sym_rhs) + eq4b = divide_terms(eq4, sym_lhs, sym_rhs) + + if naturalBC[0]: + # Natural BC (second deriv is zero at both ends) + eq5 = Eq(d2si.subs(i,0).subs(t[0],0), 0) + else: + # Specified first derivatives at boundaries + eq5 = Eq(dsi.subs(i,0).subs(t[0],0), firstDeriv[0]) + + + if naturalBC[1]: + eq6 = Eq(d2si.subs(t[i],x[i+1]-x[i]).subs(i,n-1), 0) + else: + eq6 = Eq(dsi.subs(t[i],x[i+1]-x[i]).subs(i,n-1), firstDeriv[1]) + + + + nval = nknots - 1 + neqn = 4*nval + m = Matrix.eye(neqn, neqn) + rhs = Matrix.eye(neqn,1) + + symlist = list() + for idx in range(nval): + symlist.append(a[idx]) + symlist.append(b[idx]) + symlist.append(c[idx]) + symlist.append(d[idx]) + #print(symlist) + + for idx,sym in enumerate(symlist): + m[0,idx] = get_coeff_for(eq5.lhs.subs(i,0), sym, symlist) + m[1,idx] = get_coeff_for(eq1.lhs.subs(i,0), sym, symlist) + m[2,idx] = get_coeff_for(eq2.lhs.subs(i,0), sym, symlist) + rhs[0] = eq5.rhs.subs(i,0) + rhs[1] = eq1.rhs.subs(i,0) + rhs[2] = eq2.rhs.subs(i,0) + + jdx = 3 + for nv in range(1,nval): + for idx,sym in enumerate(symlist): + m[jdx+0,idx] = get_coeff_for(eq1.lhs.subs(i,nv), sym, symlist) + m[jdx+1,idx] = get_coeff_for(eq2.lhs.subs(i,nv), sym, symlist) + m[jdx+2,idx] = get_coeff_for(eq3b.lhs.subs(i,nv-1), sym, symlist) + m[jdx+3,idx] = get_coeff_for(eq4b.lhs.subs(i,nv-1), sym, symlist) + rhs[jdx+0] = eq1.rhs.subs(i,nv) + rhs[jdx+1] = eq2.rhs.subs(i,nv) + rhs[jdx+2] = eq3b.rhs.subs(i,nv) + rhs[jdx+3] = eq4b.rhs.subs(i,nv) + jdx += 4 + + for idx,sym in enumerate(symlist): + m[jdx,idx] = get_coeff_for(eq6.lhs.subs(n,nval), sym , symlist) + rhs[jdx] = eq6.rhs.subs(n,nval) + + #for idx in range(4*nval): + # print('m = ',m.row(idx)) + #print('m = ',m) + + sln=m.LUsolve(rhs) + + subslist={} + for idx in range(nval): + subslist[a[idx]] = sln[4*idx + 0] + subslist[b[idx]] = sln[4*idx + 1] + subslist[c[idx]] = sln[4*idx + 2] + subslist[d[idx]] = sln[4*idx + 3] + + spline = dict() + for idx in range(nval): + spline[idx] = si.subs(i,idx).subs(subslist) + return spline + + +# Create C++ brace initializer list from a python list +def convert_to_brace_list(xvals, vertical_layout=False, constructor=None): + joiner = "," + if vertical_layout: + joiner = ",\n" + s = "{" + if constructor is None: + s += joiner.join([str(x) for x in xvals]) + else: + # This assumes the representation (str(x)) already encloses the value in parentheses. + # This is true if the entries in xvals are tuples. + s += joiner.join([constructor + str(x) for x in xvals]) + s+= "}" + return s; + +# -------------------------------------------------------------- +# Tests of the coefficients from the cubic spline solver routine +# -------------------------------------------------------------- + +spline_solver_template =""" +// Generated from gen_cubic_spline.py +TEST_CASE("spline_function_{test_case}", "[numerics]") +{{ + const int n = {n}; + double x[n] = {xvals}; + double y[n] = {yvals}; + double y2[n]; + + NRCubicSpline(x, y, n, {leftDeriv}, {rightDeriv}, y2); + +{y2_checks} + +}} +""" + + +# Test the coefficient values produced by the cubic spline solver +def generate_cubic_spline_coeff_test_case(n, xvals, yvals, naturalBC=(True,True), + firstDeriv=(0.0, 0.0), test_idx=0): + sln = create_solution_for_val(n, naturalBC=naturalBC, firstDeriv=firstDeriv) + + x_init_expr = convert_to_brace_list(xvals) + y_init_expr = convert_to_brace_list(yvals) + + ysubs = {y[i]:yv for i,yv in enumerate(yvals)} + xsubs = {x[i]:xv for i,xv in enumerate(xvals)} + + nd = 1e33 # signal for natural derivative + + checks_str = "" + for i in range(n-1): + si2 = diff(sln[i], t[i], 2) + sval = si2.subs(xsubs).subs(ysubs).subs(t[i], 0.0) + #print('sval = ',sval) + check = " REQUIRE(y2[%d] == Approx(%g));"%(i, sval) + checks_str += check + "\n" + + # Last second derivative must be obtained from the end of the last interval + si2 = diff(sln[n-2], t[n-2], 2) + sval = si2.subs(xsubs).subs(ysubs).subs(t[n-2], xvals[n-1]-xvals[n-2]) + #print('last sval = ',sval) + check = " REQUIRE(y2[%d] == Approx(%g));"%(n-1, sval) + checks_str += check + "\n" + + leftDeriv = nd if naturalBC[0] else firstDeriv[0] + rightDeriv = nd if naturalBC[1] else firstDeriv[1] + out = spline_solver_template.format(test_case=test_idx, n=n, xvals=x_init_expr, yvals=y_init_expr, + leftDeriv=leftDeriv, rightDeriv=rightDeriv,y2_checks=checks_str) + print(out) + + +# Create multiple test cases with different number of knots and different boundary conditions +def generate_cubic_spline_coeff_test_cases(): + xvals = [0.0, 1.0, 2.0] + yvals = [1.0, 2.0, 1.5] + generate_cubic_spline_coeff_test_case(len(xvals), xvals, yvals, test_idx=1) + generate_cubic_spline_coeff_test_case(len(xvals), xvals, yvals, naturalBC=(True,False), firstDeriv=(0.0, 1.0), test_idx=2) + generate_cubic_spline_coeff_test_case(len(xvals), xvals, yvals, naturalBC=(False,False), firstDeriv=(1.0, 2.0), test_idx=3) + + xvals = [0.0, 1.2, 2.4, 3.0] + yvals = [1.0, 2.0, 1.5, 1.8] + generate_cubic_spline_coeff_test_case(len(xvals), xvals, yvals, test_idx=4) + generate_cubic_spline_coeff_test_case(len(xvals), xvals, yvals, naturalBC=(False,False), firstDeriv=(0.0, 1.0), test_idx=5) + + + +# --------------------------------------------- +# Tests of the cubic spline evaluation routines +# --------------------------------------------- + +spline_evaluate_data_structure = """ +// Structure for holding value, first and second derivatives +struct D2U +{ + D2U(double val_, double du_, double d2u_) : val(val_), du(du_), d2u(d2u_) {} + double val; + double du; + double d2u; +}; +""" + + +spline_evaluate_template =""" +// Generated from gen_cubic_spline.py +TEST_CASE("one_dim_cubic_spline_{test_case}", "[numerics]") +{{ + + const int n = {n}; + std::vector yvals = {yvals}; + + LinearGrid grid; + grid.set(0.0, 2.0, n); + + OneDimCubicSpline cubic_spline(&grid, yvals); + + int imin = 0; + int imax = {n}-1; + + double yp0 = {leftDeriv}; + double ypn = {rightDeriv}; + + cubic_spline.spline(imin, yp0, imax, ypn); + + std::vector check_xvals = {check_xvals}; + std::vector check_yvals; + std::vector check_yvals_d2u; + + for (int i = 0; i < check_xvals.size(); i++) {{ + double r = check_xvals[i]; + double val = cubic_spline.splint(r); + check_yvals.push_back(val); + + double du, d2u; + double val2 = cubic_spline.splint(r, du, d2u); + check_yvals_d2u.push_back(D2U(val2, du, d2u)); + + //std::cout << i << " r = " << r << " val = " << val << " " << check_yvals[i] << std::endl; + }} + +{output_check} + +}} +""" + +def generate_cubic_spline_evaluation_tests(): + xvals = [0.0, 1.0, 2.0] + yvals = [1.0, 2.0, 1.5] + leftDeriv = 1.0 + rightDeriv = 2.0 + + n = len(xvals) + + sln = create_solution_for_val(n, naturalBC=(False, False), firstDeriv=(leftDeriv, rightDeriv)) + + ysubs = {y[i]:yv for i,yv in enumerate(yvals)} + xsubs = {x[i]:xv for i,xv in enumerate(xvals)} + + checks_str = "" + for i in range(n-1): + si2 = diff(sln[i], t[i], 2) + sval = si2.subs(xsubs).subs(ysubs).subs(t[i], 0.0) + #print('2nd deriv %i %g'%(i,sval)) + + # Shrink the range slightly so the last point doesn't fall exactly on the boundary + grid_range = xvals[-1] - xvals[0] - 0.0000000001 + ncheck = 6 + delta = grid_range/(ncheck-1) + check_xvals = [] + check_yvals = [] + check_y2vals = [] + check_y3vals = [] + for i in range(ncheck): + rval = i*delta + xvals[0] + check_xvals.append(rval) + # Find the correct spline interval + tval = 0.0 + for idx in range(n-1): + tval = rval - xvals[idx] + if rval >= xvals[idx] and rval < xvals[idx+1]: + break + + # Compute spline value and derivatives + val = sln[idx].subs(xsubs).subs(ysubs).subs(t[idx], tval) + dval = diff(sln[idx],t[idx]).subs(xsubs).subs(ysubs).subs(t[idx], tval) + d2val = diff(sln[idx],t[idx],2).subs(xsubs).subs(ysubs).subs(t[idx], tval) + + check_yvals.append(val) + check_y2vals.append((val,dval,d2val)) + #print(rval,val) + + y_expr = convert_to_brace_list(yvals) + x_check_expr = convert_to_brace_list(check_xvals, vertical_layout=True) + #y_check_expr = convert_to_brace_list(check_yvals, vertical_layout=True) + #y2_check_expr = convert_to_brace_list(check_y2vals, vertical_layout=True, constructor='D2U') + + output_check = '' + for i in range(ncheck): + output_check += " REQUIRE(check_yvals[%d] == Approx(%12.8g));\n"%(i,check_yvals[i]) + + output_check += '\n' + + for i in range(ncheck): + output_check += " REQUIRE(check_yvals_d2u[%d].val == Approx(%12.8g));\n"%(i,check_y2vals[i][0]) + output_check += " REQUIRE(check_yvals_d2u[%d].du == Approx(%12.8g));\n"%(i,check_y2vals[i][1]) + output_check += " REQUIRE(check_yvals_d2u[%d].d2u == Approx(%12.8g));\n"%(i,check_y2vals[i][2]) + output_check += '\n' + + test_idx = 1 + out = spline_evaluate_template.format(test_case=test_idx, n=n, yvals=y_expr, + check_xvals=x_check_expr, + leftDeriv=leftDeriv, rightDeriv=rightDeriv, output_check=output_check) + + print(spline_evaluate_data_structure) + print("\n") + print(out) + +if __name__ == '__main__': + # Put these in test_one_dim_cubic_spline.cpp, and then run through clang-format + generate_cubic_spline_coeff_test_cases() + generate_cubic_spline_evaluation_tests() diff --git a/src/Numerics/tests/test_nr_spline.cpp b/src/Numerics/tests/test_nr_spline.cpp deleted file mode 100644 index 61e8cc141..000000000 --- a/src/Numerics/tests/test_nr_spline.cpp +++ /dev/null @@ -1,39 +0,0 @@ -////////////////////////////////////////////////////////////////////////////////////// -// This file is distributed under the University of Illinois/NCSA Open Source License. -// See LICENSE file in top directory for details. -// -// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers. -// -// File developed by: Mark Dewing, markdewing@gmail.com, University of Illinois at Urbana-Champaign -// -// File created by: Mark Dewing, markdewing@gmail.com, University of Illinois at Urbana-Champaign -////////////////////////////////////////////////////////////////////////////////////// - - -#define CATCH_CONFIG_MAIN -#include "catch.hpp" -#include "Numerics/NRSplineFunctions.h" - -#include -#include - -using std::string; - -namespace qmcplusplus -{ - -TEST_CASE("NR_spline_functions", "[numerics]") -{ - - const int n = 2; - double x[n] = {0.0, 1.0}; - double y[n] = {1.0, 2.0}; - double y2[n]; - - NRCubicSpline(x, y, n, 1.0, 2.0, y2); - - REQUIRE(y2[0] != 0.0); -} - -} - diff --git a/src/Numerics/tests/test_one_dim_cubic_spline.cpp b/src/Numerics/tests/test_one_dim_cubic_spline.cpp new file mode 100644 index 000000000..eed74bf42 --- /dev/null +++ b/src/Numerics/tests/test_one_dim_cubic_spline.cpp @@ -0,0 +1,178 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2019 Jeongnim Kim and QMCPACK developers. +// +// File developed by: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory +// +// File created by: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory +////////////////////////////////////////////////////////////////////////////////////// + +#define CATCH_CONFIG_MAIN +#include "catch.hpp" +#include "Numerics/OneDimCubicSpline.h" + +#include +#include + +using std::string; + +namespace qmcplusplus +{ +// Generated from gen_cubic_spline.py +TEST_CASE("spline_function_1", "[numerics]") +{ + const int n = 3; + double x[n] = {0.0, 1.0, 2.0}; + double y[n] = {1.0, 2.0, 1.5}; + double y2[n]; + + NRCubicSpline(x, y, n, 1e+33, 1e+33, y2); + + REQUIRE(y2[0] == Approx(0)); + REQUIRE(y2[1] == Approx(-2.25)); + REQUIRE(y2[2] == Approx(0)); +} + + +// Generated from gen_cubic_spline.py +TEST_CASE("spline_function_2", "[numerics]") +{ + const int n = 3; + double x[n] = {0.0, 1.0, 2.0}; + double y[n] = {1.0, 2.0, 1.5}; + double y2[n]; + + NRCubicSpline(x, y, n, 1e+33, 1.0, y2); + + REQUIRE(y2[0] == Approx(0)); + REQUIRE(y2[1] == Approx(-3.85714)); + REQUIRE(y2[2] == Approx(6.42857)); +} + + +// Generated from gen_cubic_spline.py +TEST_CASE("spline_function_3", "[numerics]") +{ + const int n = 3; + double x[n] = {0.0, 1.0, 2.0}; + double y[n] = {1.0, 2.0, 1.5}; + double y2[n]; + + NRCubicSpline(x, y, n, 1.0, 2.0, y2); + + REQUIRE(y2[0] == Approx(2.75)); + REQUIRE(y2[1] == Approx(-5.5)); + REQUIRE(y2[2] == Approx(10.25)); +} + + +// Generated from gen_cubic_spline.py +TEST_CASE("spline_function_4", "[numerics]") +{ + const int n = 4; + double x[n] = {0.0, 1.2, 2.4, 3.0}; + double y[n] = {1.0, 2.0, 1.5, 1.8}; + double y2[n]; + + NRCubicSpline(x, y, n, 1e+33, 1e+33, y2); + + REQUIRE(y2[0] == Approx(0)); + REQUIRE(y2[1] == Approx(-2.12121)); + REQUIRE(y2[2] == Approx(2.23485)); + REQUIRE(y2[3] == Approx(0)); +} + + +// Generated from gen_cubic_spline.py +TEST_CASE("spline_function_5", "[numerics]") +{ + const int n = 4; + double x[n] = {0.0, 1.2, 2.4, 3.0}; + double y[n] = {1.0, 2.0, 1.5, 1.8}; + double y2[n]; + + NRCubicSpline(x, y, n, 0.0, 1.0, y2); + + REQUIRE(y2[0] == Approx(3.60507)); + REQUIRE(y2[1] == Approx(-3.04348)); + REQUIRE(y2[2] == Approx(2.31884)); + REQUIRE(y2[3] == Approx(1.34058)); +} + + +// Structure for holding value, first and second derivatives +struct D2U +{ + D2U(double val_, double du_, double d2u_) : val(val_), du(du_), d2u(d2u_) {} + double val; + double du; + double d2u; +}; + + +// Generated from gen_cubic_spline.py +TEST_CASE("one_dim_cubic_spline_1", "[numerics]") +{ + const int n = 3; + std::vector yvals = {1.0, 2.0, 1.5}; + + LinearGrid grid; + grid.set(0.0, 2.0, n); + + OneDimCubicSpline cubic_spline(&grid, yvals); + + int imin = 0; + int imax = 3 - 1; + + double yp0 = 1.0; + double ypn = 2.0; + + cubic_spline.spline(imin, yp0, imax, ypn); + + std::vector check_xvals = {0.0, 0.39999999998, 0.79999999996, 1.19999999994, 1.59999999992, 1.9999999999}; + std::vector check_yvals; + std::vector check_yvals_d2u; + + for (int i = 0; i < check_xvals.size(); i++) + { + double r = check_xvals[i]; + double val = cubic_spline.splint(r); + check_yvals.push_back(val); + + double du, d2u; + double val2 = cubic_spline.splint(r, du, d2u); + check_yvals_d2u.push_back(D2U(val2, du, d2u)); + + //std::cout << i << " r = " << r << " val = " << val << " " << check_yvals[i] << std::endl; + } + + REQUIRE(check_yvals[0] == Approx(1)); + REQUIRE(check_yvals[1] == Approx(1.532)); + REQUIRE(check_yvals[2] == Approx(1.976)); + REQUIRE(check_yvals[3] == Approx(1.836)); + REQUIRE(check_yvals[4] == Approx(1.352)); + REQUIRE(check_yvals[5] == Approx(1.5)); + + REQUIRE(check_yvals_d2u[0].val == Approx(1)); + REQUIRE(check_yvals_d2u[0].du == Approx(1)); + REQUIRE(check_yvals_d2u[0].d2u == Approx(2.75)); + REQUIRE(check_yvals_d2u[1].val == Approx(1.532)); + REQUIRE(check_yvals_d2u[1].du == Approx(1.44)); + REQUIRE(check_yvals_d2u[1].d2u == Approx(-0.55)); + REQUIRE(check_yvals_d2u[2].val == Approx(1.976)); + REQUIRE(check_yvals_d2u[2].du == Approx(0.56)); + REQUIRE(check_yvals_d2u[2].d2u == Approx(-3.85)); + REQUIRE(check_yvals_d2u[3].val == Approx(1.836)); + REQUIRE(check_yvals_d2u[3].du == Approx(-1.16)); + REQUIRE(check_yvals_d2u[3].d2u == Approx(-2.35)); + REQUIRE(check_yvals_d2u[4].val == Approx(1.352)); + REQUIRE(check_yvals_d2u[4].du == Approx(-0.84)); + REQUIRE(check_yvals_d2u[4].d2u == Approx(3.95)); + REQUIRE(check_yvals_d2u[5].val == Approx(1.5)); + REQUIRE(check_yvals_d2u[5].du == Approx(2)); + REQUIRE(check_yvals_d2u[5].d2u == Approx(10.25)); +} + +} // namespace qmcplusplus From 1bf145f87d7b3fac42bf4b3a738e9a9c80859db5 Mon Sep 17 00:00:00 2001 From: Mark Dewing Date: Sat, 2 Feb 2019 11:27:47 -0600 Subject: [PATCH 2/2] Expand comments in cubic spline test scripts Explain more about the purpose and usage of the cubic spline unit test generation script. --- src/Numerics/tests/eqn_manip.py | 12 +++++++++++- src/Numerics/tests/gen_cubic_spline.py | 16 +++++++++++++++- 2 files changed, 26 insertions(+), 2 deletions(-) diff --git a/src/Numerics/tests/eqn_manip.py b/src/Numerics/tests/eqn_manip.py index e8e3bc3d6..dec47e26c 100644 --- a/src/Numerics/tests/eqn_manip.py +++ b/src/Numerics/tests/eqn_manip.py @@ -2,6 +2,14 @@ from __future__ import print_function from sympy import Eq +# Perform equation manipulations as one might use when working by hand to put +# equations into a particular form - e.g. moving all the terms of one type to +# one side. +# There is also a function to extract coefficients for a particular symbol, +# used for creating the matrix of coefficients from a set of equations. + +# Used by scripts that symbolically derive equations (cubic splines, etc) + # Move symbols in sym_list from left hand side of equation to right hand side def move_terms(eqn, sym_list): new_lhs = eqn.lhs @@ -22,6 +30,7 @@ def move_terms_left(eqn, sym_list): new_rhs = new_rhs - c*sym return Eq(new_lhs, new_rhs) +# Move all there terms in the symbol lists to the respective side of the equation def divide_terms(eqn, sym_list_left, sym_list_right): #print 'start',eqn eqn1 = move_terms(eqn, sym_list_right) @@ -33,7 +42,8 @@ def divide_terms(eqn, sym_list_left, sym_list_right): def mult_eqn(eqn, e): return Eq(eqn.lhs*e, eqn.rhs*e) -# for all values other than the target, set to zero +# Extract coefficient from an expression for the symbol 'sym'. +# Works by setting all values in symlist to zero, except the target in 'sym'. def get_coeff_for(expr, sym, symlist): # expression # symbol to get coefficient for diff --git a/src/Numerics/tests/gen_cubic_spline.py b/src/Numerics/tests/gen_cubic_spline.py index 5516dd30f..d33b662ad 100644 --- a/src/Numerics/tests/gen_cubic_spline.py +++ b/src/Numerics/tests/gen_cubic_spline.py @@ -3,9 +3,23 @@ from __future__ import print_function from sympy import * from eqn_manip import * -# See for more explanation, see +# Solve for cubic spline coefficients using a straightforward (but inefficient) derivation +# from the defining equations. + +# Generates unit tests for the spline coefficient solver (NRCubicSpline in NRSplineFunctions.h) +# and the evaluation routines (OneDimCubicSpline in OneDimCubicSpline.h) +# Run this script when these unit tests need updating or expanding. +# (The computation for the 4-knot case may take a few minutes) +# Put the output of this script into test_one_dim_cubic_spline.cpp, and then run through +# clang-format to clean it up. + +# See for more explanation of the derivation, see # https://github.com/QMCPACK/qmc_algorithms/blob/master/Wavefunctions/Cubic%20Splines%20Basic.ipynb + + +# Symbols useful enough to be global + # Distance to knot, for non-uniform knots t = IndexedBase('t')