{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
"<script>\n",
" window.dataLayer = window.dataLayer || [];\n",
" function gtag(){dataLayer.push(arguments);}\n",
" gtag('js', new Date());\n",
"\n",
" gtag('config', 'UA-59152712-8');\n",
"</script>\n",
"\n",
"# Symbolic Partial Derivative Routine\n",
"\n",
"## Authors: Zach Etienne & Tyler Knowles\n",
"\n",
"## This module contains a routine for computing partial derivatives of a mathematical expression that is written as several subexpressions. Leveraging SymPy syntax, it offers a robust platform for creating derivative code by performing common subexpression elimination (CSE) on Hamiltonian expressions.\n",
"\n",
"**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
"\n",
"**Validation Notes:** This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented [below](#code_validation). Additionally, this notebook has been validated by checking that results are consistent with exact derivative expressions used in the SEOBNRv3_opt appointment of [LALSuite](https://git.ligo.org/lscsoft/lalsuite).\n",
"\n",
"### NRPy+ Source Code for this module: [SEOBNR_Derivative_Routine.py](../edit/SEOBNR/SEOBNR_Derivative_Routine.py)\n",
"\n",
"## Introduction\n",
"$$\\label{intro}$$\n",
"\n",
"This notebook documents the symbolic partial derivative routine used to generate analytic derivatives of the [SEOBNRv3](https://git.ligo.org/lscsoft/lalsuite) Hamiltonian (documented [here](../Tutorial-SEOBNR_v3_Hamiltonian.ipynb)) and described in [this article](https://arxiv.org/abs/1803.06346). In general, this notebook takes as input a file of inter-dependent mathematical expressions (in SymPy syntax), a file listing the names of values within those expressions, and a file listing all variables with which to take partial derivatives of each expression. The output is a text file containing the original expression and those for each partial derivative computation. The intention is to perform CSE on these expressions to create efficient partial derivative code!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='toc'></a>\n",
"\n",
"# Table of Contents\n",
"$$\\label{toc}$$\n",
"\n",
"This notebook is organized as follows\n",
"\n",
"1. [Step 1](#initializenrpy): Initialize core Python/NRPy+ modules\n",
"1. [Step 2:](#read_expressions) Read in Hamiltonian expressions from `Hamstring.txt`\n",
"1. [Step 3:](#list_constants) Specify constants and variables in Hamiltonian expression\n",
"1. [Step 4:](#list_free_symbols) Extract free symbols\n",
"1. [Step 5:](#convert_to_func) Convert variables to function notation; e.g., `var` goes to `var(xx)`\n",
"1. [Step 6:](#differentiate) Differentiate with respect to `xx`\n",
"1. [Step 7:](#remove_zeros) Remove derivatives (of constants) that evaluate to zero, simplifying derivative expressions\n",
"1. [Step 8:](#partial_derivative) Simplify derivatives with respect to a specific variable\n",
"1. [Step 9:](#store_results) Store partial derivatives to SymPy notebook `partial_derivatives.txt-VALIDATION.txt`\n",
"1. [Step 10:](#code_validation) Validate against LALSuite and trusted `SEOBNR_Derivative_Routine` NRPy+ module\n",
"1. [Step 11:](#latex_pdf_output) Output this notebook to $\\LaTeX$-formatted PDF file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='initializenrpy'></a>\n",
"\n",
"# Step 1: Initialize core Python/NRPy+ modules \\[Back to [top](#toc)\\]\n",
"$$\\label{initializenrpy}$$\n",
"\n",
"Let's start by importing all the needed modules from Python/NRPy+ and creating the output directory (if it does not already exist):"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:22:48.882278Z",
"iopub.status.busy": "2021-03-07T17:22:48.881327Z",
"iopub.status.idle": "2021-03-07T17:22:49.205019Z",
"shell.execute_reply": "2021-03-07T17:22:49.204239Z"
}
},
"outputs": [],
"source": [
"# Step 1.a: import all needed modules from Python/NRPy+:\n",
"import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
"import sys, os # Standard Python modules for multiplatform OS-level functions\n",
"from outputC import superfast_uniq, lhrh # Remove duplicate entries from a Python array; store left- and right-\n",
" # hand sides of mathematical expressions\n",
"\n",
"# Step 1.b: Check for a sufficiently new version of SymPy (for validation)\n",
"# Ignore the rc's and b's for release candidates & betas.\n",
"sympy_version = sp.__version__.replace('rc', '...').replace('b', '...')\n",
"sympy_version_decimal = float(sympy_version.split(\".\")[0]) + float(sympy_version.split(\".\")[1])/10.0\n",
"if sympy_version_decimal < 1.2:\n",
" print('Error: NRPy+ does not support SymPy < 1.2')\n",
" sys.exit(1)\n",
"\n",
"# Step 1.c: Name of the directory containing the input file\n",
"inputdir = \"SEOBNR\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='read_expressions'></a>\n",
"\n",
"# Step 2: Read in Hamiltonian expressions from `Hamstring.txt` \\[Back to [top](#toc)\\]\n",
"$$\\label{read_expressions}$$\n",
"\n",
"We read in the expressions of which we will compute partial derivatives in a single large string before splitting the string by line (carriage return) and by \"=\". Doing so allows us to manipulate the right- and left-hand sides of the expressions appropriately. We store the left- and right-hand sides in the array `lr`, which consists of `lhrh` arrays with left-hand sides `lhs` and right-hand sides `rhs`. Note that `Lambda` is a protected keyword in Python, so the variable $\\Lambda$ in the Hamiltonian is renamed `Lamb`."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:22:49.215900Z",
"iopub.status.busy": "2021-03-07T17:22:49.212777Z",
"iopub.status.idle": "2021-03-07T17:22:49.777402Z",
"shell.execute_reply": "2021-03-07T17:22:49.776841Z"
}
},
"outputs": [],
"source": [
"# Step 2.a: Read in expressions as a (single) string\n",
"with open(os.path.join(inputdir,'Hamstring.txt'), 'r') as file:\n",
" expressions_as_lines = file.readlines()\n",
"\n",
"# Step 2.b: Create and populate the \"lr\" array, which separates each line into left- and right-hand sides\n",
"# Each entry is a string of the form lhrh(lhs='',rhs='')\n",
"lr = []\n",
"\n",
"for i in range(len(expressions_as_lines)):\n",
" # Ignore lines with 2 or fewer characters and those starting with #\n",
" if len(expressions_as_lines[i]) > 2 and expressions_as_lines[i][0] != \"#\":\n",
" # Split each line by its equals sign\n",
" split_line = expressions_as_lines[i].split(\"=\")\n",
" # Append the line to \"lr\", removing spaces, \"sp.\" prefixes, and replacing Lambda->Lamb\n",
" # (Lambda is a protected keyword):\n",
" lr.append(lhrh(lhs=split_line[0].replace(\" \",\"\").replace(\"Lambda\",\"Lamb\"),\n",
" rhs=split_line[1].replace(\" \",\"\").replace(\"sp.\",\"\").replace(\"Lambda\",\"Lamb\")))\n",
"\n",
"# Step 2.c: Separate and sympify right- and left-hand sides into separate arrays\n",
"lhss = []\n",
"rhss = []\n",
"\n",
"# As of April 2021, \"sp.sympify(\"Q+1\")\" fails because Q is a reserved keyword.\n",
"# This is the workaround, courtesy Ken Sible.\n",
"custom_global_dict = {}\n",
"exec('from sympy import *', custom_global_dict)\n",
"del custom_global_dict['Q']\n",
"if sympy_version_decimal >= 1.6:\n",
" custom_parse_expr = lambda expr: sp.parse_expr(expr, global_dict=custom_global_dict)\n",
"else:\n",
" custom_parse_expr = lambda expr: sp.sympify(expr)\n",
"\n",
"for i in range(len(lr)):\n",
" lhss.append(custom_parse_expr(lr[i].lhs))\n",
" rhss.append(custom_parse_expr(lr[i].rhs))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='list_constants'></a>\n",
"\n",
"# Step 3: Specify constants and variables in Hamiltonian expression \\[Back to [top](#toc)\\]\n",
"$$\\label{list_constants}$$\n",
"\n",
"We read in and declare as SymPy symbols the constant values; derivatives with respect to these variables will be set to zero. We then read in the variables with respect to which we want to take derivatives and declare those as SymPy variables as well."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:22:49.784164Z",
"iopub.status.busy": "2021-03-07T17:22:49.783560Z",
"iopub.status.idle": "2021-03-07T17:22:49.785792Z",
"shell.execute_reply": "2021-03-07T17:22:49.786259Z"
}
},
"outputs": [],
"source": [
"# Step 3.a: Create `input_constants` array and populate with SymPy symbols\n",
"m1,m2,tortoise,eta,KK,k0,k1,EMgamma,d1v2,dheffSSv2 = sp.symbols('m1 m2 tortoise eta KK k0 k1 EMgamma d1v2 dheffSSv2',\n",
" real=True)\n",
"input_constants = [m1,m2,tortoise,eta,KK,k0,k1,EMgamma,d1v2,dheffSSv2]\n",
"\n",
"# Step 3.b: Create `dynamic_variables` array and populate with SymPy symbols\n",
"x,y,z,px,py,pz,s1x,s1y,s1z,s2x,s2y,s2z = sp.symbols('x y z px py pz s1x s1y s1z s2x s2y s2z', real=True)\n",
"dynamic_variables = [x,y,z,px,py,pz,s1x,s1y,s1z,s2x,s2y,s2z]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='list_free_symbols'></a>\n",
"\n",
"# Step 4: Extract free symbols \\[Back to [top](#toc)\\]\n",
"$$\\label{list_free_symbols}$$\n",
"\n",
"By ''free symbols'' we mean the variables in the right-hand sides. We first create a list of all such terms (using SymPy's built-in free_symbol attribute), including duplicates, and then strip the duplicates. We then remove input constants from the symbol list."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:22:49.848133Z",
"iopub.status.busy": "2021-03-07T17:22:49.825046Z",
"iopub.status.idle": "2021-03-07T17:22:49.850966Z",
"shell.execute_reply": "2021-03-07T17:22:49.850273Z"
}
},
"outputs": [],
"source": [
"# Step 4.a: Prepare array of \"free symbols\" in the right-hand side expressions\n",
"full_symbol_list_with_dups = []\n",
"for i in range(len(lr)):\n",
" for variable in rhss[i].free_symbols:\n",
" full_symbol_list_with_dups.append(variable)\n",
"\n",
"# Step 4.b: Remove duplicate free symbols\n",
"full_symbol_list = superfast_uniq(full_symbol_list_with_dups)\n",
"\n",
"# Step 4.c: Remove input constants from symbol list\n",
"for inputconst in input_constants:\n",
" for symbol in full_symbol_list:\n",
" if str(symbol) == str(inputconst):\n",
" full_symbol_list.remove(symbol)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='convert_to_func'></a>\n",
"\n",
"# Step 5: Convert variables to function notation; e.g., `var` goes to `var(xx)` \\[Back to [top](#toc)\\]\n",
"$$\\label{convert_to_func}$$\n",
"\n",
"In order to compute the partial derivative of each right-hand side, we mark each variable (left-hand side) and each free symbol (in right-hand sides) as a function with argument $\\texttt{xx}$."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:22:49.925636Z",
"iopub.status.busy": "2021-03-07T17:22:49.889818Z",
"iopub.status.idle": "2021-03-07T17:22:53.604953Z",
"shell.execute_reply": "2021-03-07T17:22:53.605440Z"
}
},
"outputs": [],
"source": [
"# Step 5.a: Convert each left-hand side to function notation\n",
"# while separating and simplifying left- and right-hand sides\n",
"xx = sp.Symbol('xx',real=True)\n",
"func = []\n",
"for i in range(len(lr)):\n",
" func.append(sp.sympify(sp.Function(lr[i].lhs,real=True)(xx)))\n",
"\n",
"# Step 5.b: Mark each free variable as a function with argument xx\n",
"full_function_list = []\n",
"for symb in full_symbol_list:\n",
" func = sp.sympify(sp.Function(str(symb),real=True)(xx))\n",
" full_function_list.append(func)\n",
" for i in range(len(rhss)):\n",
" for var in rhss[i].free_symbols:\n",
" if str(var) == str(symb):\n",
" rhss[i] = rhss[i].subs(var,func)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='differentiate'></a>\n",
"\n",
"# Step 6: Differentiate with respect to `xx` \\[Back to [top](#toc)\\]\n",
"$$\\label{differentiate}$$\n",
"\n",
"Now we differentiate the right-hand expressions with respect to `xx`. We use the SymPy $\\texttt{diff}$ command, differentiating with respect to $\\texttt{xx}$. After so doing, we remove $\\texttt{(xx)}$ and \"Derivative\" (which is output by $\\texttt{diff}$), and use \"prm\" suffix to denote the derivative with respect to $\\texttt{xx}$."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:22:53.680022Z",
"iopub.status.busy": "2021-03-07T17:22:53.644184Z",
"iopub.status.idle": "2021-03-07T17:22:56.040496Z",
"shell.execute_reply": "2021-03-07T17:22:56.041098Z"
}
},
"outputs": [],
"source": [
"# Step 6: Use SymPy's diff function to differentiate right-hand sides with respect to xx\n",
"# and append \"prm\" notation to left-hand sides\n",
"lhss_deriv = []\n",
"rhss_deriv = []\n",
"for i in range(len(rhss)):\n",
" lhss_deriv.append(custom_parse_expr(str(lhss[i])+\"prm\"))\n",
" newrhs = custom_parse_expr(str(sp.diff(rhss[i],xx)).replace(\"(xx)\",\"\").replace(\", xx\",\"prm\").replace(\"Derivative\",\"\"))\n",
" rhss_deriv.append(newrhs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='remove_zeros'></a>\n",
"\n",
"# Step 7: Remove derivatives (of constants) that evaluate to zero, simplifying derivative expressions \\[Back to [top](#toc)\\]\n",
"$$\\label{remove_zeros}$$\n",
"\n",
"We declare a function to simply the derivative expressions. In particular, we want to remove terms equal to zero."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:22:56.115471Z",
"iopub.status.busy": "2021-03-07T17:22:56.079868Z",
"iopub.status.idle": "2021-03-07T17:22:56.420571Z",
"shell.execute_reply": "2021-03-07T17:22:56.420043Z"
}
},
"outputs": [],
"source": [
"# Step 7.a: Define derivative simplification function\n",
"def simplify_deriv(lhss_deriv,rhss_deriv):\n",
" # Copy expressions into another array\n",
" lhss_deriv_simp = []\n",
" rhss_deriv_simp = []\n",
" for i in range(len(rhss_deriv)):\n",
" lhss_deriv_simp.append(lhss_deriv[i])\n",
" rhss_deriv_simp.append(rhss_deriv[i])\n",
" # If a right-hand side is 0, substitute value 0 for the corresponding left-hand side in later terms\n",
" for i in range(len(rhss_deriv_simp)):\n",
" if rhss_deriv_simp[i] == 0:\n",
" for j in range(i+1,len(rhss_deriv_simp)):\n",
" for var in rhss_deriv_simp[j].free_symbols:\n",
" if str(var) == str(lhss_deriv_simp[i]):\n",
" rhss_deriv_simp[j] = rhss_deriv_simp[j].subs(var,0)\n",
" zero_elements_to_remove = []\n",
" # Create array of indices for expressions that are zero\n",
" for i in range(len(rhss_deriv_simp)):\n",
" if rhss_deriv_simp[i] == sp.sympify(0):\n",
" zero_elements_to_remove.append(i)\n",
"\n",
" # When removing terms that are zero, we need to take into account their new index (after each removal)\n",
" count = 0\n",
" for i in range(len(zero_elements_to_remove)):\n",
" del lhss_deriv_simp[zero_elements_to_remove[i]+count]\n",
" del rhss_deriv_simp[zero_elements_to_remove[i]+count]\n",
" count -= 1\n",
" return lhss_deriv_simp,rhss_deriv_simp\n",
"\n",
"# Step 7.b: Call the simplication function and then copy results\n",
"lhss_deriv_simp,rhss_deriv_simp = simplify_deriv(lhss_deriv,rhss_deriv)\n",
"lhss_deriv = lhss_deriv_simp\n",
"rhss_deriv = rhss_deriv_simp"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='partial_derivative'></a>\n",
"\n",
"# Step 8: Simplify derivatives with respect to a specific variable \\[Back to [top](#toc)\\]\n",
"$$\\label{partial_derivative}$$\n",
"\n",
"In [Step 6](#differentiate) we took a generic derivative of each expression, assuming all variables were functions of `xx`. We now define a function that will select a specific dynamic variable (element of `dynamic_variables`) and set the derivative of the variable to 1 and all others to 0."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:22:56.534205Z",
"iopub.status.busy": "2021-03-07T17:22:56.498465Z",
"iopub.status.idle": "2021-03-07T17:23:22.940124Z",
"shell.execute_reply": "2021-03-07T17:23:22.939324Z"
}
},
"outputs": [],
"source": [
"# Step 8.a: Define onevar derivative function\n",
"def deriv_onevar(lhss_deriv,rhss_deriv,variable_list,index):\n",
" # Denote each variable with prm\n",
" variableprm_list = []\n",
" for variable in variable_list:\n",
" variableprm_list.append(str(variable)+\"prm\")\n",
"\n",
" # Copy expressions into another array\n",
" lhss_deriv_new = []\n",
" rhss_deriv_new = []\n",
" for i in range(len(rhss_deriv)):\n",
" lhss_deriv_new.append(lhss_deriv[i])\n",
" rhss_deriv_new.append(rhss_deriv[i])\n",
" # For each free symbol's derivative, replace it with:\n",
" # 1, if we are differentiating with respect to the variable, or\n",
" # 0, if we are note differentiating with respect to that variable\n",
" for i in range(len(rhss_deriv_new)):\n",
" for var in variableprm_list:\n",
" if variableprm_list.index(str(var))==index:\n",
" rhss_deriv_new[i] = rhss_deriv_new[i].subs(var,1)\n",
" else:\n",
" rhss_deriv_new[i] = rhss_deriv_new[i].subs(var,0)\n",
" # Simplify derivative expressions again\n",
" lhss_deriv_simp,rhss_deriv_simp = simplify_deriv(lhss_deriv_new,rhss_deriv_new)\n",
" return lhss_deriv_simp,rhss_deriv_simp\n",
"\n",
"# Step 8.b: Call the derivative function and populate dictionaries with the result\n",
"lhss_derivative = {}\n",
"rhss_derivative = {}\n",
"for index in range(len(dynamic_variables)):\n",
" lhss_temp,rhss_temp = deriv_onevar(lhss_deriv,rhss_deriv,dynamic_variables,index)\n",
" lhss_derivative[dynamic_variables[index]] = lhss_temp\n",
" rhss_derivative[dynamic_variables[index]] = rhss_temp"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='store_results'></a>\n",
"\n",
"# Step 9: Store partial derivatives to SymPy notebook `partial_derivatives.txt-VALIDATION.txt` \\[Back to [top](#toc)\\]\n",
"$$\\label{store_results}$$\n",
"\n",
"We write the resulting derivatives in SymPy syntax. Each partial derivative is output in its own file, in a similar format to the input expressions."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:23:22.950621Z",
"iopub.status.busy": "2021-03-07T17:23:22.949063Z",
"iopub.status.idle": "2021-03-07T17:23:24.016404Z",
"shell.execute_reply": "2021-03-07T17:23:24.015807Z"
}
},
"outputs": [],
"source": [
"# Step 9: Output original expression and each partial derivative expression in SymPy snytax\n",
"with open(os.path.join(inputdir,'partial_derivatives.txt-VALIDATION'), 'w') as output:\n",
" for i in range(len(lr)):\n",
" right_side = lr[i].rhs\n",
" right_side_in_sp = right_side.replace(\"sqrt(\",\"sp.sqrt(\").replace(\"log(\",\"sp.log(\").replace(\"pi\",\n",
" \"sp.pi\").replace(\"sign(\",\"sp.sign(\").replace(\"Abs(\",\n",
" \"sp.Abs(\").replace(\"Rational(\",\"sp.Rational(\")\n",
" output.write(str(lr[i].lhs)+\" = \"+right_side_in_sp)\n",
" for var in dynamic_variables:\n",
" for i in range(len(lhss_derivative[var])):\n",
" right_side = str(rhss_derivative[var][i])\n",
" right_side_in_sp = right_side.replace(\"sqrt(\",\"sp.sqrt(\").replace(\"log(\",\"sp.log(\").replace(\"pi\",\n",
" \"sp.pi\").replace(\"sign(\",\"sp.sign(\").replace(\"Abs(\",\n",
" \"sp.Abs(\").replace(\"Rational(\",\"sp.Rational(\").replace(\"prm\",\n",
" \"prm_\"+str(var))\n",
" output.write(str(lhss_derivative[var][i]).replace(\"prm\",\"prm_\"+str(var))+\" = \"+right_side_in_sp+\"\\n\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='code_validation'></a>\n",
"\n",
"# Step 10: Validate against LALSuite and trusted `SEOBNR_Derivative_Routine` NRPy+ module \\[Back to [top](#toc)\\]\n",
"$$\\label{code_validation}$$\n",
"\n",
"We validate the output of this notebook against known LALSuite values of the Hamiltonian partial derivatives and the output of the `SEOBNR_Derivative_Routine` NRPy+ module. We note that due to cancellations in the derivative terms, various versions of SymPy may result in relative errors that differ as much as an order of magnitude. Furthermore, even changing the set of input parameters can affect the relative error by as many as two orders of magnitude. Therefore we look for agreement with LALSuite to at least 10 significant digits.\n",
"\n",
"When comparing the notebook output to that of the NRPy+ module, we compare term-by-term using SymPy to check that each right-hand side is equivalent."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:23:24.085946Z",
"iopub.status.busy": "2021-03-07T17:23:24.051401Z",
"iopub.status.idle": "2021-03-07T17:24:15.454913Z",
"shell.execute_reply": "2021-03-07T17:24:15.455364Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Computing the difference between notebook output and trusted LALSuite derivative values...\n",
"Test passed: the notebook agrees with LALSuite to at least 10 significant digits!\n",
"Printing difference between notebook output and trusted NRPy+ module output...\n",
"Checking file partial_derivatives.txt-VALIDATION\n",
"No difference. TEST PASSED!\n"
]
}
],
"source": [
"# Define a function to return a set of reasonable input parameters\n",
"# This function contains three distinct sets of input parameters, and index differentiates between them\n",
"def reset_values(tort_value, index):\n",
" # Check that a reasonable tortoise value has been passed\n",
" if tort_value!=1 and tort_value!=2:\n",
" print(\"Error: tortoise must be 1 or 2.\")\n",
" sys.exit(1)\n",
" # Each index corresponds to a specific set of input parameters\n",
" if index==0:#-M 13 -m 11 -X 0.1 -Y -0.2 -Z 0.3 -x -0.3 -y 0.2 -z -0.1\n",
" values = {'m1': 1.300000000000000e+01, 'm2': 1.100000000000000e+01, 'eta': 2.482638888888889e-01,\n",
" 'x': 1.658426645098320e+01, 'y': 3.975021008701605e-02, 'z': -1.820682538442627e-07,\n",
" 's1x': 2.934751675254397e-02, 's1y': -5.867672205485316e-02, 's1z': 8.802097562761332e-02,\n",
" 's2x': -6.302678133897792e-02, 's2y': 4.200490780215727e-02, 's2z': -2.100705983874398e-02,\n",
" 'KK': 3.913980338468737e-01, 'k0': -7.447639215330089e-01, 'k1': -6.380586501824999e-01,\n",
" 'd1v2': -7.476323019145448e+01,'dheffSSv2':2.105103187692902e+01,\n",
" 'EMgamma': 0.577215664901532860606512090082402431}\n",
" # Note that we transform the momentum based on the tortoise values\n",
" if tort_value==1:\n",
" values.update({'px': -1.517631642228534e-03, 'py': 2.693180445886167e-01, 'pz': -1.320499830947482e-04,\n",
" 'tortoise': 1})\n",
" else:\n",
" values.update({'px': -1.633028076483384e-03, 'py': 2.693177679992048e-01, 'pz': -1.320499918278832e-04,\n",
" 'tortoise': 2})\n",
" elif index==1:#-M 25 -m 10 -X 0.1 -Y -0.0 -Z 0.1 -x -0.2 -y 0.0 -z -0.2\n",
" values = {'m1': 2.500000000000000e+01, 'm2': 1.000000000000000e+01, 'eta': 2.040816326530612e-01,\n",
" 'x': 1.289689003662444e+01, 'y': 5.495441315063273e-03, 'z': -1.717482806041791e-11,\n",
" 's1x': 5.102040816179230e-02, 's1y': 9.846215537206260e-07, 's1z': 5.102040816473832e-02,\n",
" 's2x': -1.632653061189792e-02, 's2y': -6.762952223804450e-07, 's2z': -1.632653061259188e-02,\n",
" 'KK': 5.642540639599580e-01, 'k0': -1.063532077165767e+00, 'k1': -8.835684149841774e-01,\n",
" 'd1v2': -8.041179092044979e+01,'dheffSSv2':1.125986130778842e+01,\n",
" 'EMgamma': 0.577215664901532860606512090082402431}\n",
" if tort_value==1:\n",
" values.update({'px': -1.898773926867491e-03, 'py': 3.160984442121970e-01, 'pz': 1.171602901570564e-07,\n",
" 'tortoise': 1})\n",
" else:\n",
" values.update({'px': -2.209215477700561e-03, 'py': 3.160983119312114e-01, 'pz': 1.171602905704723e-07,\n",
" 'tortoise': 2})\n",
" elif index==2:#-M 7 -m 5 -X 0.01 -Y -0.5 -Z 0.03 -x -0.04 -y 0.05 -z -0.06\n",
" values = {'m1': 7.000000000000000e+00, 'm2': 5.000000000000000e+00, 'eta': 2.430555555555556e-01,\n",
" 'x': 2.633506161699224e+01, 'y': 7.574563213724998e-02, 'z': -2.789625823248071e-08,\n",
" 's1x': 3.417297286269225e-03, 's1y': -1.701385963191495e-01, 's1z': 1.020835932957879e-02,\n",
" 's2x': -6.945454346305877e-03, 's2y': 8.679766617922793e-03, 's2z': -1.041665076794264e-02,\n",
" 'KK': 4.052853693162246e-01, 'k0': -7.706473492549312e-01, 'k1': -6.587426366263742e-01,\n",
" 'd1v2': -7.555647472993827e+01,'dheffSSv2':1.972817669753086e+01,\n",
" 'EMgamma': 0.577215664901532860606512090082402431}\n",
" if tort_value==1:\n",
" values.update({'px': -7.883793607066706e-04, 'py': 2.068742709904638e-01, 'pz': -7.338789145500886e-04,\n",
" 'tortoise': 1})\n",
" else:\n",
" values.update({'px': -8.039726989861640e-04, 'py': 2.068742261404732e-01, 'pz': -7.338789145335709e-04,\n",
" 'tortoise': 2})\n",
" else:\n",
" # If an improper index is passed, exit\n",
" print(\"Error: invalid index (only three sets of input parameters available).\")\n",
" sys.exit(1)\n",
" # Return the input values\n",
" return values\n",
"\n",
"# Numerically evaluate right-hand sides using input values\n",
"def evaluate_expression(left_sides,right_sides,input_values):\n",
" new_right_sides = []\n",
" for i in range(len(right_sides)):\n",
" term = custom_parse_expr(str(right_sides[i]).replace(\"(xx)\",\"\"))\n",
" # Only look for the free variables in each expression to reduce computation time\n",
" free_vars = term.free_symbols\n",
" for variable in free_vars:\n",
" term = term.subs(variable, input_values[str(variable)])\n",
" # Evaluate each term to reduce computation time\n",
" new_right_sides.append(sp.sympify(term.evalf()))\n",
" # Store each subexpression in values numerically\n",
" input_values[str(left_sides[i])] = new_right_sides[i]\n",
" # Return the input values dictionary with all numerical right-hand added\n",
" return input_values\n",
"\n",
"# Create array of trusted LALSuite derivative values\n",
"# Note that position in the array corresponds to the index of the corresponding input values\n",
"LALSuite_validated_values = []\n",
"#-M 13 -m 11 -X 0.1 -Y -0.2 -Z 0.3 -x -0.3 -y 0.2 -z -0.1\n",
"LALSuite_validated_values.append({'Hreal': 9.928923110195770e-01,'dHreal_dx': 9.932484846748471e-04,\n",
" 'dHreal_dy': 2.813294366789505e-06, 'dHreal_dz': 1.926378549762488e-06,\n",
" 'dHreal_dpx': -3.710666135737856e-04, 'dHreal_dpy': 6.116199124763537e-02,\n",
" 'dHreal_dpz': -5.600910364542288e-07, 'dHreal_ds1x': -1.438467658934620e-05,\n",
" 'dHreal_ds1y': -1.319462868057848e-06, 'dHreal_ds1z': 7.665413183773232e-04,\n",
" 'dHreal_ds2x': -2.075691477548065e-05,'dHreal_ds2y': 2.456427688083135e-06,\n",
" 'dHreal_ds2z': 8.762835349889455e-04})\n",
"#-M 25 -m 10 -X 0.1 -Y -0.0 -Z 0.1 -x -0.2 -y 0.0 -z -0.2\n",
"LALSuite_validated_values.append({'Hreal': 9.926852598351464e-01, 'dHreal_dx': 1.397519118422771e-03,\n",
" 'dHreal_dy': 1.928133240540033e-06, 'dHreal_dz': -1.215449398950413e-06,\n",
" 'dHreal_dpx': -4.004159849919695e-04, 'dHreal_dpy': 5.701850933742150e-02,\n",
" 'dHreal_dpz': 4.329487960716782e-08, 'dHreal_ds1x': 2.259457049322466e-06,\n",
" 'dHreal_ds1y': -2.544122765762015e-09, 'dHreal_ds1z': 9.834156257814124e-04,\n",
" 'dHreal_ds2x': 5.185557993931246e-06,'dHreal_ds2y': 2.437768415468806e-10,\n",
" 'dHreal_ds2z': 2.111169766641698e-03})\n",
"#-M 7 -m 5 -X 0.01 -Y -0.5 -Z 0.03 -x -0.04 -y 0.05 -z -0.06\n",
"LALSuite_validated_values.append({'Hreal': 9.955293642650920e-01, 'dHreal_dx': 3.734697245297603e-04,\n",
" 'dHreal_dy': 1.105998063449349e-06, 'dHreal_dz': 5.367207414282669e-08,\n",
" 'dHreal_dpx': -1.848412708548443e-04, 'dHreal_dpy': 4.754239153769983e-02,\n",
" 'dHreal_dpz': -3.549083643069269e-08, 'dHreal_ds1x': -4.819261725948465e-07,\n",
" 'dHreal_ds1y': 3.333280059627902e-06, 'dHreal_ds1z': 2.201786563823208e-04,\n",
" 'dHreal_ds2x': -7.576810957551029e-07,'dHreal_ds2y': 6.818093508597533e-06,\n",
" 'dHreal_ds2z': 2.922663340179887e-04})\n",
"\n",
"# Sort variables by which tortoise value we use to compute the derivatives\n",
"variables_tort2 = [x,y,z]\n",
"variables_tort1 = [px,py,pz,s1x,s1y,s1z,s2x,s2y,s2z]\n",
"\n",
"# Call evaluation function\n",
"print(\"Computing the difference between notebook output and trusted LALSuite derivative values...\")\n",
"for index in range(3):\n",
" values = reset_values(2,index)\n",
" values = evaluate_expression(lhss,rhss,values)\n",
" Hreal = values['Hreal']\n",
" trusted_Hreal = LALSuite_validated_values[index]['Hreal']\n",
" relative_difference = (trusted_Hreal - Hreal)/trusted_Hreal\n",
" if abs(relative_difference) > 1e-9:\n",
" print(\"The relative difference in Hreal is too large: %.15e\" % relative_difference)\n",
" sys.exit(1)\n",
"\n",
" for var in variables_tort2:\n",
" Hrealprm = evaluate_expression(lhss_derivative[var],rhss_derivative[var],values)['Hrealprm']\n",
" trusted_Hrealprm = LALSuite_validated_values[index]['dHreal_d'+str(var)]\n",
" relative_difference = (trusted_Hrealprm - Hrealprm)/trusted_Hrealprm\n",
" if abs(relative_difference) > 1e-9:\n",
" print(\"The relative difference in dHreal_d%s is too large: %.15e\" % (var,relative_difference))\n",
" sys.exit(1)\n",
"\n",
" values = reset_values(1,index)\n",
" values = evaluate_expression(lhss,rhss,values)\n",
"\n",
" for var in variables_tort1:\n",
" Hrealprm = evaluate_expression(lhss_derivative[var],rhss_derivative[var],values)['Hrealprm']\n",
" trusted_Hrealprm = LALSuite_validated_values[index]['dHreal_d'+str(var)]\n",
" relative_difference = (trusted_Hrealprm - Hrealprm)/trusted_Hrealprm\n",
" if abs(relative_difference) > 1e-9:\n",
" print(\"The relative difference in dHreal_d%s is too large: %.15e\" % (var,relative_difference))\n",
" sys.exit(1)\n",
"print(\"Test passed: the notebook agrees with LALSuite to at least 10 significant digits!\")\n",
"\n",
"print(\"Printing difference between notebook output and trusted NRPy+ module output...\")\n",
"# Open the files to compare\n",
"file = 'partial_derivatives.txt'\n",
"outfile = 'partial_derivatives.txt-VALIDATION'\n",
"\n",
"print(\"Checking file \" + outfile)\n",
"with open(os.path.join(inputdir,file), \"r\") as file1, open(os.path.join(inputdir,outfile), \"r\") as file2:\n",
" # Read the lines of each file\n",
" file1_lines = file1.readlines()\n",
" file2_lines = file2.readlines()\n",
" # Compare right-hand sides of the expressions by computing the difference between them\n",
" num_diffs = 0\n",
" for i in range(len(file1_lines)):\n",
" expr_new = custom_parse_expr(file1_lines[i].split(\"=\")[1].replace(\"sp.\",\"\"))\n",
" expr_validated = custom_parse_expr(file2_lines[i].split(\"=\")[1].replace(\"sp.\",\"\"))\n",
" difference = sp.simplify(expr_new - expr_validated)\n",
" if difference != 0:\n",
" num_diffs += 1\n",
" print(difference)\n",
" if num_diffs == 0:\n",
" print(\"No difference. TEST PASSED!\")\n",
" else:\n",
" print(\"WARNING: Possible disagreement found with the trusted file. See differences above.\")\n",
"# sys.exit(1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='latex_pdf_output'></a>\n",
"\n",
"# Step 11: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
"$$\\label{latex_pdf_output}$$\n",
"\n",
"The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n",
"[Tutorial-SEOBNR_Derivative_Routine.pdf](Tutorial-SEOBNR_Derivative_Routine.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:24:15.459650Z",
"iopub.status.busy": "2021-03-07T17:24:15.459038Z",
"iopub.status.idle": "2021-03-07T17:24:18.942926Z",
"shell.execute_reply": "2021-03-07T17:24:18.943661Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Created Tutorial-SEOBNR_Derivative_Routine.tex, and compiled LaTeX file to\n",
" PDF file Tutorial-SEOBNR_Derivative_Routine.pdf\n"
]
}
],
"source": [
"import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n",
"cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-SEOBNR_Derivative_Routine\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.1"
}
},
"nbformat": 4,
"nbformat_minor": 4
}