<jupyter_start><jupyter_text>CS207 Systems Development Final Project:  Automatic Differentiation package: genericdiff TopCoderKitty-ML**Collaborators**: <NAME>, <NAME> and <NAME>   Introduction  OverviewDerivatives play an integral role in computational science, ranging from its use in gradient descent, Newton's method, to finding the posteriors of Bayesian models. We discuss numerical differentiation, symbolic differentiation, how both demonstrate limitations, and automatic differentiation, the focus of our software. We acknowledge its effectiveness in both its accuracy and efficiency when evaluating derivatives, provide real life applications in the biological context, as well as a guide to how to use our software.  Motivation for Automatic DifferentiationBecause functions are often too complex to solve analytically, instead, we look to alternative methods that automatically calculate derivatives. There are three main ways to approach this issue: numerical differentiation from finding finite difference approximations, symbolic differentiation through expression manipulation, and automatic differentiation (AD or algorithmic differentiation). While numerical differentiation is easy to code, it is also subject to floating point errors; symbolic differentiation gives exact and accurate results, but is too computationally expensive. Thus, automatic differentiation proves to be the most effective method as it works to resolve both of these issues; AD is both exact/numerically stable and computationally efficient.   Background What is AD?Conceptually straightforward, automatic differentiation can be defined as a family of techniques that evaluate the derivative through the use of elementary arithmetic operations (ie. addition, subtraction, etc.), elementary arithmetic functions (ie. exp, log, sin), and the chain rule. AD is a recursive process that involves repeatedly taking the chain rule to the elementary operations at hand, and allows us to calculate the individual components of the gradient (a list of partial derivatives in terms of each of the inputs) evaluations, to produce results that are automatic and precise. Because AD involves a specific family of techniques that compute derivatives through accumulation of values during code execution to generate numerical derivative evaluations rather than derivative expressions, it can attain machine precision. There are two modes in AD: the forward mode and reverse mode.  Forward ModeThe forward mode begins at the innermost portion of the function and repeatedly, or recursively, applies the chain rule while traversing out. Thus, the forward pass creates the evaluation trace, which is a composition of the finite set of elementary operations for which derivatives are known. This is then combined to evaluate the derivative of the overall composition. Notably, the derivative at subsequent steps are calculated based on the derivatives calculated in preceding steps. It is also important to note that the forward pass also finds the derivatives and the values of the partial derivatives at each step, thus requiring the setting of the seed vector, which indicates which input variable to the take the partial derivative in terms of. Taking the derivative of the m dependent output variables in terms of a single independent input variable make up one column of the jacobian matrix. Thus, the full jacobian matrix can be defined as the partial derivative of the m output variables in terms of the n input variables, or applying the forward pass across n evaluations. These recursive steps can be documented in a table, and visually represented through a computational graph.The forward mode can be simplified by utilizing another important component of automatic differentiation: dual numbers. Dual numbers are a type of number that uses  and allows for simultaneously automatically differentiating a function while also evaluating the value of the function. The forward mode is efficient and straight- forward because it is able to compute all the derivatives in terms of one input with just one forward pass. Reverse modeIn the backward mode, a forward pass creates the evaluation trace and indicates the partial derivatives at each step, but does not find the values of the partial derivatives. At the end of this process, the final node’s derivative is evaluated by using an arbitrary seed. Then, the values of the partial derivatives that constitute the end node’s derivative are found by performing a backward pass through the tree to get all the values. During both the forward and backward modes, all intermediate variables are evaluated, and their values are stored; these steps can be represented in a table, and further visualized in a computational graph. The graph (and table) essentially outlines this repeated chain rule process; it also serves as the basis of the logic behind our automatic differentiation software library. Jacobian productAuto differentiation can take in a vector of input functions with multiple variables. In this case, the autodifferentiation algorithm returns the Jacobian Product matrix which is just a matrix where each row represents the partial derivatives of the variables of a function within the vector. If the vector has m functions then the Jacobian product will have m rows. If the functions contain n variables, the Jacobian will contain n columns. The Jacobian product matrix is handy in first order derivative applications in the sciences when you are dealing with lots of different systems (functions) and unknowns in the systems (variables). We will see an application of this in the RNA velocity package we will be building.  Application of ADAD can be applied to many branches of computational science, ranging from areas in biology to politics. We chose to focus our application in a biological context, namely RNA velocity. This involves a system of differential equations with pretty complicated solutions, that is able to predict the future state of individual cells through single cell RNA gene expression. Estimating this is important in aiding the analysis of developmental lineages and cellular dynamics, particularly in humans.  Software organizationOur implentation philosophy was to create an object such that all the elemental functions and constants could be arbritrarily combined and nested for differentiation.We created a ```genericdiff``` library that contains the following modules:- ```__init__.py``` : Initializes the genericdiff package - ```generic_diff.py``` : Contains GenericDiff class where any specialized classes such as sine or cos class were reduced to just value and derivative attributes such that we could use our overloaded operators. +, -, /, * , **, and negation were defined according to to the following differentiation rules: chain rule, generalized power rule, quotient rule and product rule. It also overloaded comparison operators, which include , =, ==, and !=.- ```elemental_functions.py```: Specialized classes in the following modules that inherited from GenericDiff- exponential, sin, cos, tan, sinh, cosh, tanh, acos, asin, atan, exp, log, logit, sqrt classes were defined here. - ```vector_jacobian.py```: A class to handle taking in a vector of functions, multiple inputs and producing a jacobian product matrixWe created a test package to test our generic_diff, elemental_functions, and vector_jacobian modules for all methods and error handling:- ```test_elemental_functions.py```- ```test_generic_diff.py```- ```test_generic_diff_comparisons.py```- ```test_vector_jacobian.py```The test suite is run through pytest and called through Travis CI using the travis yml file - it sends a report to codecov for code coverage.The driver scripts are located in /driverscripts is:- ```driver_univariate.py``` This solves for roots using newton's method - the univariate case is using GenericDiff class- ```driver_multiivariate.py``` This solves for roots using newton's method - the multivariate case is using JacobianProduct class- ```driver_jacobian_product.py``` Finds a jacobian product matrix - uses JacobianProduct class  ImplementationWe wanted to keep our implementation lightweight. Therefore, we do not use any imports except python's standard math library for our single function, single input use case. In our jacobian vector product use case, we will use numpy.The core data structure is our GenericDiff class which has at its base, a value attribute and a derivative attribute. This can be combined or taken in by a more specialized class that differentiates specialized mathematical functions like sine, cosine, exp, etc. We have explained how this works based on the software organization above. Given this data structure, we can arbritrarily combine objects to represent extremely complex mathematical functions.In our jacobian product matrix we will utilize arrays to take care of multiple inputs.Our core classes are:- JacobianProduct: This class takes in function vectors and allows us to find partial derivatives and jacobian products given a set of input values. This is the user-facing class that takes care of multiple function, multiple input use cases.- GenericDiff: This is the underlying class that powers JacobianProduct. It takes care of uni-variate differentiation, all overloaded operator functions, variable and constant instantiation rules    - Var: this is a variable class that inherits from GenericDiff. This identifies a variable in a function expression.    - Constant: This is a constant class that inherits from GenericDiff. This converts any floats or integers in function expressions into a constant to assign constants with a derivative value of 0- Elemental functions: These classes are taken care of in the elemental_functions.py moduleThe important attributes in our GenericDiff class are:- .val = value- .der = derivative- NOTE the end user does not need to interact directly with these, as the user can just use the JacobianProduct class and call methods in that class to find the values and derivatives.The external dependencies are:- math- numpyElementary functions covered: - sin- cos- tan- sinh- cosh- tanh- acos- asin- atan- log- logit- sqrt- exp- powers- multiplication- division- addition- subtraction- unary negation  How to use```genericdiff``` differentiates by using the forward mode of automatic differentiation.To do this, one just has to instantiate the ```JacobianProduct``` object, which takes in vectors of functions with multiple inputs to return the jacobian product matrices. ```JacobianProduct``` classes can take in single function, single inputs as well.\* Developers looking to customize further, can use the single variable, single input ```GenericDiff``` class by instantiating a Var object. One can then combine Var objects with constants using various operations to arrive at derivative evaluations.   Getting setupInstalling the library is straight- forward. 1. Install the following from pip using the command line.<jupyter_code>python -m pip install rnavelocity-genericdiff<jupyter_output><empty_output><jupyter_text>2. Next, import the package:<jupyter_code>import sys
sys.path.append("../")
from genericdiff import *<jupyter_output><empty_output><jupyter_text>Class and method examples  JacobianProductThe JacobianProduct class is the user-facing tool that a layperson should use. This takes care of all single, multi input and multi function use cases. We will highlight some code showing how our underlying GenericDiff class (used for the single variate case and powers JacobianProduct) works later in this documentation but it is not something a user needs to call directly.The way in which this package can differentiate automatically is through the instantiation of a JacobianProduct class, which contains a vector of functions as its input. The partial method calculates the partial derivative with respect to a given variable, and the jacobian_product method gets the jacobian product matrix with respect to all variables. Here is a demo of the partial method being called with respect to x, as well as the jacobian_product method:<jupyter_code>f = lambda x, y: x**2 - y**3
h = lambda x, y: x**3 + y**3

function_vector = [f, h]

jp_object = JacobianProduct(function_vector)
# getting partial (derivative) with respect to x (index 0 in lambdas)
inputs = [[1, 2, 3], 0]
partial_der_wrt = jp_object.partial_ders(wrt=0, inputs=inputs)
print(partial_der_wrt)
# getting partial (value) with respect to x (index 0 in lambdas)
inputs = [[1, 2, 3], 0]
partial_val_wrt = jp_object.partial_vals(wrt=0, inputs=inputs)
print(partial_val_wrt)
f = lambda x, y: x**2 - y**3
h = lambda x, y: x**3 + y**3

function_vector = [f, h]
jp_object = JacobianProduct(function_vector)

# inputs are x = {1, 2, 3} and y = {1, 2, 3}
# this means we will calculate derivatives for
# (1, 1), (2, 2), and (3, 3)

inputs = [[1, 2, 3], [1, 2, 3]]

# getting jp matrix with respect to all variables
jp_matrix = jp_object.jacobian_product(inputs)

print(jp_matrix)<jupyter_output>[array([[ 2., -3.],
       [ 3.,  3.]]), array([[  4., -12.],
       [ 12.,  12.]]), array([[  6., -27.],
       [ 27.,  27.]])]<jupyter_text>GenericDiffNOTE this class powers JacobianProduct and might be useful for developers, but all functionality of this class is captured in JacobianProduct in a more user friendly way.We can instantiate our GenericDiff objects for the uni-variate use case using the ```Var``` class to instantiate a variable. We can retrieve the .val and .der attributes from the GenericDiff object.Let's demonstrate this using our operator functions. We can combine two objects together (variable or constant) and get their values and derivatives using the mathematical operators:1. subtraction -2. addition +3. division /4. multiplication *5. power **6. negation -( )Here we demo the + operator for the single function single input use case:<jupyter_code>h = Var(4) + 2
print("The value of x + 2 is:", h.val)
print("The derivative of x + 2 is:", h.der)<jupyter_output>The value of x + 2 is: 6
The derivative of x + 2 is: 1<jupyter_text>Here we demo the power operator:<jupyter_code>g = Var(4) ** 2

print("The value of x ** 2 is:", g.val)
print("The derivative of x ** 2 is:", g.der)<jupyter_output>The value of x ** 2 is: 16
The derivative of x ** 2 is: 8.0<jupyter_text>We can also apply trigonometric and exp functions, which inherit from the GenericDiff class.The available functions are:We can also apply trigonometric and exp functions, which inherit from the GenericDiff class.The available functions are:1. sine    ```sin(x)```2. cosine  ```cos(x)```3. tangent ```tan(x)```4. $e^x$   ```exponential(x)```5. hyperbolic sine ```sinh(x)```6. hyperbolic cosine ```cosh(x)```7. hyperbolic tangent ```tanh(x)```8. arc sine  ```asin(x)```9. arc cosine ```acos(x)```10. arc tangent ```atan(x)```11. log ```log(x, base=e)``` numpy default base is e, but user can set any base they want12. logit ```logit(x)``` $\frac{e^x}{1+e^x}$13. square root ```sqrt(x)``` $\sqrt{x}$Here we demo the sinh function by applying it to the x + 2 function we created above:<jupyter_code>j = sinh(Var(4) + 2)

print("The value of sinh(x+2) is:", j.val)
print("The derivative of sinh(x+2) is:", j.der)<jupyter_output>The value of sinh(x+2) is: 201.71315737027922
The derivative of sinh(x+2) is: 201.7156361224559<jupyter_text>Note that the derivative in this case is ```j.der```. It is a scalar value since we are only doing single input single functions. We can further complicate this by exponentiating:<jupyter_code>g = exp(j)

print("The value of exp(sinh(x+2)) is:", g.val)
print("The derivative of exp(sinh(x+2)) is:", g.der)<jupyter_output>The value of exp(sinh(x+2)) is: 4.007847362293387e+87
The derivative of exp(sinh(x+2)) is: 8.084454801667175e+89<jupyter_text>We can use the comparison operators to compare the derivatives of the objects or compmare their derivatives to a value<jupyter_code>print(j == g)
print(j != g)
print(j <= g)
print(j >= g)
print(j > g)
print(j < g)
print(j > 201)<jupyter_output>False
True
True
False
False
True
True<jupyter_text>Advanced featuresWe have applied our package to creating a biology package focused on ordering cells given gene expression values. Please read the documentation in the ```/rnavelocity folder``` for more information. Here is a brief primer. rnavelocity: derivative of single-cell gene expression stateThis package calculates the RNA velocities for a given set of genes and their cells using their expression values. We can think of RNA velocities as first order derivatives specified from a functional form. The applications of RNA velocities are:1. To serve as important and powerful indicator of future state of individual cells 2. To distinguish between unspliced and spliced mRNAs in common single-cell RNA sequencing protocols3. To predict future state of individual cells on a timescale of hours4. To aid analysis of developmental lineages and cellular dynamics, particularly in humansIn the rnavelocity package we specifically attack application 4 where we will optimize for parameters:1. alpha $\alpha$ (transcription rate)2. gamma $\gamma$ (decay rate of the product)We can then plug these optimized parameters into a time sorting algorithm to find the correct ordering of the lineage of cells given their genes. In simpler terms, we want to answer the chicken or egg question: "which cells were reactants that allowed other cells to develop down the line?"  User Documentation  JacobianProduct Class The JacobianProduct class contains the user facing methods for single variate, multi variate, and multi function use cases.  Methods in JacobianProduct Class:```JacobianProduct(function_vector)```This instantiates the JacobianProduct class with the user defined functions.Parameters:        function_vector: list        a list of functions on which to differentiate. For example, function_vector = [lambda x, y: x, lambda x, y: sin(x) + cos(y)]```JacobianProduct.partial_ders(self, wrt, inputs, fun_idx=-1)```Returns the partial derivative with respect to a chosen variableParameters:        wrt: int         index of variable partial derivative is taken with respect to. For example if x is the first variable in the function declaration, setting wrt=0 will differentiate with respect to x.        inputs: list        chosen input values for variables. Note that the inputs list is sensitive to the order in which you declared your variables in the function vector. The variable with which you are differentiating can take on multiple values. For example, if there are three variables x, y, z and wrt = 0, then you can specify multiple values as such [[1, 2, 3], 0, 0]. Note that variables held constant can only assume ONE value. In this case, y and z are held constant. If non-differentiating variables (in this case, y and z) are inputted with multiple values, this method will throw an error.            fun_idx: int        index of the functions in which partial derivative is applied. default is set to -1 differentiate for all functions in the function vector. For example, if the fun_idx = 0 then it will only differentiate against the first function in the function vectorOutput:        returns a list of np arrays with partial DERIVATIVES for each function. For example, if we have two functions f and g with three x points: [[df/dx1, df/dx2, df/dx3], [dh/dx1, dh/dx2, dh/dx3]]```JacobianProduct.partial_vals(self, wrt, inputs, fun_idx=-1)```Returns the partial derivative values with respect to a chosen variable. This is a convenience function for those looking to get the function values in the same format as the partial derivative lists.Parameters:        wrt: int         index of variable partial derivative is taken with respect to. For example if x is the first variable in the function declaration, setting wrt=0 will differentiate with respect to x.        inputs: list        chosen input values for variables. Note that the inputs list is sensitive to the order in which you declared your variables in the function vector. The variable with which you are differentiating can take on multiple values. For example, if there are three variables x, y, z and wrt = 0, then you can specify multiple values as such [[1, 2, 3], 0, 0]. Note that variables held constant can only assume ONE value. In this case, y and z are held constant. If non-differentiating variables (in this case, y and z) are inputted with multiple values, this method will throw an error.             fun_idx: int        index of the functions in which partial derivative is applied. default is set to -1 differentiate for all functions in the function vector. For example, if the fun_idx = 0 then it will only differentiate against the first function in the function vector    Output:        returns a list of np arrays with partial VALUES for each function. For example, if we have two functions f and g with three x points: [[f(x1), f(x2), f(x3)], [h(x1), h(x2), h(x3)]]        ```JacobianProduct.jacobian_product(self, inputs, fun_idx=-1)```Returns list of jacobian productsParameters:    inputs: list        chosen inputs for variables. For exammple, if we have 3 inputs x, y and z. Single inputs can take the form [1, 2, 3] or [[1], [2], [3]]. Multiple inputs can take the form [[1, 2], [1, 2], [1, 2]]. If multiple inputs per variable, the length of inputs for each variable MUST be the same. If not, the method returns an error. The jacobian product will be found for all element-wise sets. In this example, a jacobian product is found for (1, 1, 1) and (2, 2, 2).            fun_idx: int         index of the functions in which jacobian product is applied. Default is set to -1 differentiate for all functions in the function vector. For example, if the fun_idx = 0 then it will only differentiate against the first function in the function vectorOutput:    returns a list of np arrays which are jacobian products for all input sets. For example, if we have the input [[1, 2], [1, 2], [1, 2]], the output will be a list of two jacobian product matrices - one for set (1, 1, 1) and one for set (2, 2, 2).  ExampleFor more examples please see the vignette section<jupyter_code>##### Partial Derivative

f = lambda x, y: x**2 - cos(y)**3
h = lambda x, y: x**3 + cos(y)**3

function_vector = [f, h]

jp_object = JacobianProduct(function_vector)
# getting partial (derivative) with respect to x (index 0 in lambdas)
inputs = [[1, 2, 3], 0]
partial_der_wrt = jp_object.partial_ders(wrt=0, inputs=inputs)

print("======Partials wrt x======")
print(partial_der_wrt)

####### Jacobian Product

# inputs are x = {1, 2, 3} and y = {1, 2, 3}
# this means we will calculate jacobian products for
# (1, 1), (2, 2), and (3, 3)

inputs = [[1, 2, 3], [1, 2, 3]]

# getting jp matrix with respect to all variables
jp_matrix = jp_object.jacobian_product(inputs)

print("======Jacobian Product=======")
print(jp_matrix)<jupyter_output>======Partials wrt x======
[[2.0, 4.0, 6.0], [3.0, 12.0, 27.0]]
======Jacobian Product=======
[array([[ 2.        ,  0.73694324],
       [ 3.        , -0.73694324]]), array([[ 4.        ,  0.47241145],
       [12.        , -0.47241145]]), array([[ 6.        ,  0.41492887],
       [27.        , -0.41492887]])]<jupyter_text>We've overloaded elemental functions such that a user can use these in their function definitions and they will be automatically converted to GenericDiff objects that are then differentiated in the JacobianProduct class. For examples, see genericdiff_vignette.ipynb. The following are available to use:  Elemental Functions1. sine    ```sin(x)```2. cosine  ```cos(x)```3. tangent ```tan(x)```4. $e^x$   ```exponential(x)```5. hyperbolic sine ```sinh(x)```6. hyperbolic cosine ```cosh(x)```7. hyperbolic tangent ```tanh(x)```8. arc sine  ```asin(x)```9. arc cosine ```acos(x)```10. arc tangent ```atan(x)```11. log ```log(x, base=e)``` numpy default base is e, but user can set any base they want12. logit ```logit(x)``` $\frac{e^x}{1+e^x}$13. square root ```sqrt(x)``` $\sqrt{x}$For the JacobianProduct class function vector, the user can define a vector using a lambda function using the elemental functions in the following way:<jupyter_code>f = lambda x: sin(x)
g = lambda x: cos(x)
function_vector = [f, g]
jp_object = JacobianProduct(function_vector)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07169873376089458, 'fraction_numerical': 0.015704653839828976, 'mean_word_length': 4.816594930655189, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 40, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '65512', 'license': 'MIT', 'n_tokens_mistral': 6480, 'n_tokens_neox': 5991, 'n_words': 3411, 'path': 'genericdiff/documentation.ipynb', 'repo_name': 'Topcoder-Kitty-ML/cs207-FinalProject', 'size': '24339'}