Vector Function Tutorial

One of the goals of ASSET is to provide users the ability to construct functions dynamically within Python that are able to be used by ASSET. By doing this we can simplify a user’s work-flow, where the benefits of high speed C++ code can be combined with the ease of use Python provides. In this section we will give an depth overview of ASSET’s VectorFunction type. At a high-level, this is simple functional (in the programming sense) domain specific language for defining mathematical Vector Functions that take a fixed number of inputs and produce a fixed number outputs, with both inputs and outputs assumed to be column vectors.

Arguments

Let us begin by importing ASSET and the VectorFunctions module which contains all types and functions defining the language. From this module, we will then import the Arguments type and give it a shorthand name. The Arguments type is the base expression in the VectorFunction system and represents a function that simply takes some number of input arguments and returns them as outputs. It always serves as the starting point for defining more complicated functional expressions involving some or all of its outputs.

vf        = ast.VectorFunctions
Args      = vf.Arguments

We can construct the object by simply specifying the number of arguments, which in this case is 6. This instance X is now a first class function thats takes any vector of size 6 and returns that vector. Since it is a VectorFunction, we can compute its output value using the () operator, first derivative using the .jacobian method, and second derivative using the .adjointhessian method. To do this, we provide either a numpy vector or python list of real valued inputs, and additionally for the second derivative, a vector list of lagrange multipliers with same dimensions as the output of the function. One important note, ASSET does not compute the full 3D tensor second derivative of vector valued functions, instead it computes the second derivative dotted with a vector lagrange multipliers, resulting in a 2D symmetric matrix with rows and columns equal to the number of inputs. We refer to this as the adjointhessian, and in the case of a function with a single output it is equivalent to the normal hessian. Since X here is a simple linear function, the first derivative is simply the identity matrix and the adjointhessian is zero. This is a rather trivial example, but the same methods can be applied to any ASSET VectorFunction that we can construct. We should also note that while these methods are available for all VectorFunctions, for most applications and examples you won’t ever actually need to explicitly call the function or its derivatives at real arguments, as that will be handled for you by some other interface such as an integrator or optimal control problem.

xvals = np.array([0,1,2,3,4,5])
lvals = np.ones((6))

X = Args(6)


print( X(xvals) ) #prints [0,1,2,3,4,5]

print( X.jacobian(xvals) ) #prints Identiy matrix of size 6

print( X.adjointhessian(xvals,lvals) ) #prints zero matrix of size 6

As you can see, Arguments itself does not do anything very interesting, but what it does do is serve as a starting point for defining functions of elements and sub-vectors. For example, we may make a new object referencing one of its elements by using the bracket operator. This will return an object of another fundamental type, Element, which is itself a function that takes all input arguments and returns the specified element as a scalar output. Attempting to address an element out of bounds of the Arguments will immediately throw an error at the offending operation.

xvals = np.array([0,1,2,3,4,5])

X = Args(6)

x0 = X[0]
x5 = X[5]

#x42 =X[42]  #throws an error

print(x0(xvals))  # prints [0.0]
print(x5(xvals))  # prints [5.0]

Often times you will write an expression where the Arguments represent many separate distinct scalar elements that will be combined to construct a more complicated expression. In this case, we can bypass the many lines necessary to index them individually, by using the .tolist() method of Arguments. This method will return all the individual elements concatenated in order inside of a single python list, which we can then unpack into individual named elements in a single line. You may use whatever method you wish, but we personally prefer the .tolist() method in this case.

xvals = np.array([0,1,2,3,4,5])

X = Args(6)

x0 = X[0]
x1 = X[1]
x2 = X[2]
x3 = X[3]
x4 = X[4]
x5 = X[5]

## Equivalent to

x0,x1,x2,x3,x4,x5 = X.tolist()

In addition to scalar Elements, one may also address contiguous sub-vectors in a set of arguments using the .head(), .tail(), and .segment() methods of Arguments, or standard python (contiguous) list indexing. For example, if we want to treat the first three arguments of the Arguments below as a single vector R, we can use the .head(n) method. The .head(n) method returns a sub-vector of size n starting at the first element. This syntax mirrors the Eigen C++ library, which we find to be quite nice, but you may also use standard python list indexing to accomplish the same goal. Similarly, if we want to address the last three arguments as a single vector V , we can use the .tail(n) method which returns the last n elements of some arguments. Finally, we can address vectors of length \(n\) starting at index i using the .segment(i,n) method. The return type of all of these methods is the fundamental Segment type, which is a function that returns as its output the specified sub-vector of the arguments.

xvals = np.array([1,2,3,4,5,6])

X = Args(6)

R = X.head(3)
R = X[0:3]    # Same as above

print(R(xvals)) #prints [1,2,3]


V = X.tail(3)
V = X[3:6]    # same as above

print(V(xvals))  #prints [4,5,6]


R = X.segment(0,3) # same as R above
V = X.segment(3,3) # same as V above


N = X.segment(1,4) # first argument is starting index, second is size
N = X[1:5]         #same as above but python style

print(N(xvals))    # prints [2,3,4,5]

Paralleling what we did before with elements, we can also partition an input argument list list into segments and elements using the .tolist([(start,size), ..]) method. In this case, we should pass a python list of tuples, where the first element of each tuple is the starting index of the sub-vector and the second is the size; sub-vectors of size one are returned as elements. Note that this method does not require you to partition all of the argument set, though this example does. Furthermore, it is not required that the sub-vectors specified be in any particular order, though we highly recommend you sort them according to starting index.

xvals = np.array([1,2,3,4,5,6,7,8])

X = Args(8)

R = X.head(3)
V = X.segment(3,3)
t = X[6]
u = X[7]

## Equivalent to the Above
R,V,t,u = X.tolist([ (0,3), (3,3), (6,1),(7,1) ])



print(R(xvals)) #prints [1,2,3]

print(V(xvals))  #prints [4,5,6]

print(t(xvals))  #prints [7]

print(u(xvals))  #prints [8]

Finally, all of the above indexing methods behave exactly the same when applied to segments rather than arguments, and we can address their individual components as elements, and split them in smaller segments. For example, we may split R into its scalar components using .tolist, address a single component using brackets, or a subsegment using head, tail, segment, etc..

xvals = np.array([1,2,3,4,5,6])

X = Args(6)

R = X.head(3)
V = X.tail(3)

r0,r1,r2 = R.tolist()

print(r0(xvals))  #prints([1])

v0 = V[0]

print(v0(xvals))  #prints([4])

V12 = V.tail(2)

print(V12(xvals))  #prints([5,6])

Standard Math Operations

Having covered most everything related to constructing arguments, and their elements and sub-vectors, we can move on the to combining them together into meaningful mathematical functions. We should note that the result of any mathematical non-indexing operation will have the generic type VectorFunction (more than one output) or ScalarFunction (one output), which themselves may be operated on and combined with the three fundamental types using the same rules. In general, types will be converted automatically, and users should not concern themselves with the types of resulting expressions and should only make sure that their expressions are mathematically consistent. We may add, subtract, multiply, and divide functions by other functions and numerical constants using the standard rules of vector math. For example, we may add or subtract two functions of the same output size to together, add or subtract vectors of constants or constant scalars, multiply functions by constant scalars, multiply functions by Scalar functions, etc.

xvals = np.array([1,2,3,4,5,6])

X = Args(6)

R = X.head(3)
V = X.tail(3)

S = R[0]*V[0]*V[1]*5.0

RpV = R + V

RmC = R - np.array([1.0,1.0,1.0])

Rtv0 = R*V[0]

RtC   = R*2

RdC   = R/2

Vdr0 = V/R[0]

N = Rtv0 + RdC

v1pv0 = (V[1]+V[0] + 9.0)*2.0

inv0 = 1.0/v0

As this is a vector math language, certain operations involving vectors are not allowed via standard multiply and divide operator overloads. For example, one may not multiply two VectorFunctions together using the * operator as is possible with two arrays in numpy. This is an explicit choice because in our opinion, for the types of expressions written using ASSET, allowing element-wise vector multiplication creates more problems in terms of incorrect problem formulation than it solves. However, these operations can be accomplished using methods we describe later. Note, this does not apply to ScalarFunctions such as Element or ScalarFunction, which may be multiplied together with no issue, and may also scale a VectorFunction.

## RmV = R*V  # Throws and Error
## RdV = R/V  # Throws and Error

Scalar Math Operations

Next, we will move on to describe the standard mathematical functions that can be applied to scalar valued functions. These encompass most of the standard functions that can be found in python or C math libraries, such as sin, cos, tan, etc. All of these functions are stored inside the VectorFunctions module (which we have imported as vf), and can be called as shown below. A complete list of functions is given in the table below.

X = Args(6)

a = vf.sin(X[0])
b = vf.cos(X[1])
c = vf.tan(X[1])

d = vf.cosh((X[1]+X[0])*X[1])

e = vf.arctan2(X[0],X[1]/3.14)

f = X[0]**2  # power operator

g = vf.abs(X[0])

h  = vf.sign(-X[1])
List of Scalar Math Functions

Function

Description

vf.sin(f)

Returns the sine of an input Element or ScalarFunction

vf.cos(f)

Returns the cosine of an input Element or ScalarFunction

vf.tan(f)

Returns the tangent of an input Element or ScalarFunction

vf.arcsin(f)

Returns the inverse sine of an input Element or ScalarFunction

vf.arccos(f)

Returns the inverse cosine of an input Element or ScalarFunction

vf.arctan(f)

Returns the inverse tangent of an input Element or ScalarFunction

vf.sinh(f)

Returns the hyperbolic sine of an input Element or ScalarFunction

vf.cosh(f)

Returns the hyperbolic cosine of an input Element or ScalarFunction

vf.tanh(f)

Returns the hyperbolic tangent of an input Element or ScalarFunction

vf.arcsinh(f)

Returns the inverse hyperbolic sine of an input Element or ScalarFunction

vf.arccosh(f)

Returns the inverse hyperbolic cosine of an input Element or ScalarFunction

vf.arctanh(f)

Returns the inverse hyperbolic tangent of an input Element or ScalarFunction

vf.log(f)

Returns the natural logarithm of an input Element or ScalarFunction

vf.exp(f)

Returns the exponential function of an input Element or ScalarFunction

vf.sqrt(f)

Returns the square root of an input Element or ScalarFunction

vf.sign(f)

Returns the sign(+1.0,-1.0) of an input Element or ScalarFunction

vf.abs(f)

Returns the absolute value an input Element or ScalarFunction

Vector Norms and Normalizations

For Vector valued functions, we also provide member functions that will compute various useful norms and transformations on vectors. While most of these could be computed using the math operations we have already covered, users should always use one of these methods if applicable, as the resulting expressions will be much faster when evaluated. A few examples are illustrated here, and a complete list of such functions is given in the table below.

X = Args(6)

R = X.head(3)
V = X.tail(3)

r   = R.norm()
r   = vf.sqrt(R[0]**2 + R[1]**2 + R[2]**2)  # Same as above but slower

v2 =  V.squared_norm()
v2 = V[0]**2 + V[1]**2 + V[2]**2 # Same as above but slower


Vhat = V.normalized()
Vhat = V/V.norm()        # Same as above but slower



r3 =  R.cubed_norm()

Grav = - R.normalized_power3()  # R/|R|^3
Grav2 = - R/r3         # Same as above but slower
List of Vector Norms and Normalizations

Function

Math Form

Description

F.norm()

\(|\vec{F}|\)

Returns the euclidean norm of VectorFunction or Segment F

F.squared_norm()

\(|\vec{F}|^2\)

Returns the square of the euclidean norm of VectorFunction or Segment F

F.cubed_norm()

\(|\vec{F}|^3\)

Returns the cube of the euclidean norm of VectorFunction or Segment F

F.inverse_norm()

\(1/|\vec{F}|\)

Returns the inverse of the euclidean norm of VectorFunction or Segment F

F.inverse_squared_norm()

\(1/|\vec{F}|^2\)

Returns the inverse square of the euclidean norm of VectorFunction or Segment F

F.inverse_cubed_norm()

\(1/|\vec{F}|^3\)

Returns the inverse cube of the euclidean norm of VectorFunction or Segment F

F.normalized()

\(\frac{\vec{F}}{|\vec{F}|}\)

Returns the normalized output of VectorFunction or Segment F

F.normalized_power2()

\(\frac{\vec{F}}{|\vec{F}|^2}\)

Returns the output of VectorFunction or Segment F divided by its euclidean norm squared.

F.normalized_power3()

\(\frac{\vec{F}}{|\vec{F}|^3}\)

Returns the output of VectorFunction or Segment F divided by its euclidean norm cubed.

F.normalized_power4()

\(\frac{\vec{F}}{|\vec{F}|^4}\)

Returns the output of VectorFunction or Segment F divided by its euclidean norm to the fourth power.

F.normalized_power5()

\(\frac{\vec{F}}{|\vec{F}|^5}\)

Returns the output of VectorFunction or Segment F divided by its euclidean norm to the fifth power.

Vector Products

In addition to the standard binary math operations supported via operator overloads, we also provide member functions and free functions for performing various common vector operations. The most commonly used are the dot, cross, quaternion, and coefficient-wise products. A few examples of how these can be used are shown below. All functions appearing in these expressions must have the correct output size, otherwise an error will be immediately thrown. You may also mix and match constant numpy arrays and VectorFunctions as needed to define your function. It should be noted that our quaternion products assume that the vector part of the quaternion is the first three components of the output while the real part is the 4th element(ie: q =[qv,q4]).

R,V,N,K = Args(14).tolist([(0,3),(3,3),(6,4),(10,4)])

C2 = np.array([1.0,1.0])
C3 = np.array([1.0,1.0,2.0])
C4 = np.array([1.0,1.0,2.0,3.0])


dRV = R.dot(V)
dRV = vf.dot(R,V)

dRC = R.dot(C3)     # use .dot with a constant vector of size 3
dRC = vf.dot(C3,R)  # Or do it with a free function

#dRC = R.dot(C4)  # throws ERROR because vector is incorrect size


RcrossV = R.cross(V)
RcrossV = vf.cross(R,V)
RcrossC3 = vf.cross(R,C3)

RcVcNdC3 = (R.cross(V)).cross(N.head(3)).dot(C3)

#RcrossC4 = vf.cross(R,C4)  # throws an error

KqpN = vf.quatProduct(K,N) # Hamiltonian quaternion product
Krn  = vf.quatRotate(K,V)  ## Rotates 3x1 vector V using quaternion K


KpN  = K.cwiseProduct(N)
NpC4 = N.cwiseProduct(C4)
KdN = K.cwiseQuotient(N)

Stacking Outputs

Up to this point, we have looked at partitioning and operating on the outputs of other functions, and have not addressed how the outputs of functions may be combined together into a larger single function. This can be accomplished using the VERY IMPORTANT vf.stack() method. In general stack takes a list of ASSET function types and produces another function whose output is the concatenation of all the outputs. There are two signatures for stack, The first one (vf.stack([f1,f2,...])) takes a python list containing only explicit ASSET function types (ie: Element, ScalarFunction , VectorFunction, Segment etc..). This version does not allow one to mix in floats or numpy vectors. The second signature (vf.stack(f1,f2,...)) does the same thing as the first but does not enclose the objects to be stacked inside of a list. Additionally, for this second signature, you may mix in arbitrary floats and numpy vectors that will be included in the output.

xvals = np.array([1,0,0,
          0,1,0])

R,V = Args(6).tolist([(0,3),(3,3)])

Rhat = R.normalized()
Nhat = R.cross(V).normalized()
That = Nhat.cross(Rhat).normalized()

RTN = vf.stack([Rhat,That,Nhat])
print(RTN(xvals))  #prints [1. 0. 0. 0. 1. 0. 0. 0. 1.]

#Err = vf.stack([Rhat,That,np.array([1.0,1.0])]) # Throws Error, numpy array not allowed

RTN = vf.stack(Rhat,That,Nhat)  # Same as above

Stuff = vf.stack(7.0, Rhat,42.0,That,Nhat, np.array([2.71,3.14]) )

print(Stuff(xvals))  #prints [ 7., 1., 0.,  0., 42., 0., 1., 0. ,0., 0. ,1. ,2.71,3.14]

Matrix Operations

While ASSET is and always will be a language for defining functions with vector valued inputs and outputs, we do have limited but growing support for interpreting VectorFunctions as matrices inside of expressions. This is supported through the vf.ColMatrix and vf.RowMatrix types. These are types constructed from some VectorFunction and interprets the outputs as \(n \times m\) matrix. A ColMatrix will interpret the coefficients of the output as a column major matrix, whereas RowMatrix interprets them as a row major matrix. Once constructed you may multiply matrices by any other appropriately sized Row/ColMatrix functions in any order, or multiply them on the right by appropriately sized VectorFunctions. The result of all matrix on matrix operations are assumed to be ColMatrix type. The result of Matrix*vector operations is VectorFunction. You may also return the transpose of an existing Row/ColMatrix. Furthermore, square matrices may be inverted resulting in a Matrix type with same row/col type. For now, one may only add matrices together if they have the same Row/Col type, though we will support adding different Row/Col types in the future.

R,V,U = Args(9).tolist([(0,3),(3,3),(6,3)])

## Three orthonormal basis vectors
Rhat = R.normalized()
Nhat = R.cross(V).normalized()
That = Nhat.cross(Rhat).normalized()

RTNcoeffs = vf.stack([Rhat,That,Nhat])

RTNmatC = vf.ColMatrix(RTNcoeffs,3,3)  # Interpret as col major 3x3 Rotation matrix
RTNmatR = vf.RowMatrix(RTNcoeffs,3,3)  # Interpret as row major 3x3 Rotation matrix

M2 = RTNmatC*RTNmatR # Multiply matrices together result is column major

U1 = RTNmatC*U       # Multiply on the right by a VectorFunction of size (3x1)
U2 = RTNmatR*U
U3 = M2*U

Zero = RTNmatR.inverse()*U -RTNmatC*U
Identity = RTNmatC*RTNmatC.transpose()

RTNmatC +RTNmatC

Conditional Statement/Operations

ASSET’s intended use case is for defining constraints, objectives, and dynamical models that will eventually be put to use inside of a second derivative optimizer. As a general rule of thumb, it is a bad idea for such functions to contain conditional statements, as this could potentially result in non-smooth derivatives. In these cases we always recommend considering whether what you were trying to accomplish with the conditional statement can be reformulated in another way. However, if this is not possible, or you are writing a function that will not see the inside of an optimizer, we do offer support for simple conditional statements and boolean operations with VectorFunction expressions. To be precise, we support constructing boolean statements involving the outputs of scalar valued functions, and then using those as conditional statements to control the output of another expression. Conditional statements are constructed by applying the comparison operators (>,<,<=,>=) to the outputs of ScalarFunctions. This can be used to dispatch one of two functions using the vf.ifelse() function as shown below. Note that the output sizes of both the true and false functions MUST be the same. Conditional statements may also be combined together using the bitwise or/and operators (|,&).

x0,x1,x2 = Args(3).tolist()

condition = x0<1.0

output_if_true = x1*2
output_if_false = x1+x2

func = vf.ifelse(condition,output_if_true,output_if_false)


print(func([0,  2,3]))  # prints [4.0]
print(func([1.5,2,3]))  # prints [5.0]


Fine = vf.ifelse(condition,vf.stack(x1,x2),vf.stack(x2,x1))
#Error = vf.ifelse(condition,vf.stack(x1,x2),output_if_false)


combo_condition = (x0<1.0)|(x0>x1)

func = vf.ifelse(combo_condition,output_if_true,output_if_false)


print(func([0,  2,3]))  # prints [4.0]
print(func([1.5,2,3]))  # prints [5.0]
print(func([2.5,2,3]))  # prints [4.0]

Some Notes on Input Arguments

Before moving on any further, we need to make one very important note about how the VectorFunction type system works. In all of our previous examples, we have created and partitioned one set of arguments of a certain size, from which we constructed other functions. You might ask, what happens if we try to mix expressions formulated out of arguments of different sizes? This is strictly not allowed, as our entire type system is predicated on the fact that expressions can only be combined if they have the same sized input arguments. For example, the following code will throw an error to alert you that you have made a mistake. However, we should also note as shown below, that there is nothing unique about any two sets of arguments of the same size. Thus, you may (though it is pointless) combine expressions derived from two arguments objects of the same size.

X1 = Args(9)
X2 = Args(12)
X3 = Args(12)

R1,V1,U1 = X1.tolist([(0,3),(3,3),(6,3)])
R2,V2,U2 = X2.tolist([(0,3),(3,3),(6,3)])
R3,V3,U3 = X3.tolist([(0,3),(3,3),(6,3)])

#Error = R1 + R2
#Error = R1.dot(V2)

## These two functions do identical things
Fine = R2.dot(V3)
Fine = R3.dot(V2)

Suggested Style and Organization

At this point we have covered most all of the operations one can and can’t perform with ASSET VectorFunctions, with the important exception of function composition (which we will cover in the next section). As you might have noticed, in all of our scratch pad examples, we simply created a single set of arguments and operated on them in the same scope. Every one of these functions is a fully formed ASSET type and can be immediately passed off to other parts of the library to be used as constraints/ODEs/controllers etc. However, obviously it is not a recipe for longterm success to simply write expressions inline wherever they are needed. How you package or encapsulate the construction of ASSET VectorFunctions is up to you, but we suggest one of the following two methods.

Method one involves simply writing a standard python function that takes as arguments any meta data or constants, needed to define the function, then writing and returning your ASSET VectorFunction. A trivial example of this is shown below, and you can find many others throughout our problem specific examples contained in other sections.

def FuncMethod(a,b,c):
        x0,x1,x2 = Args(3).tolist()
        eq1 = x0 +a - x1
        eq2 = x2*b + x1*c
        return vf.stack(eq1,eq2)

func = FuncMethod(1,2,3)

print(func([1,1,1]))  # prints [1,5]

Method two involves defining a new class that inherits from the appropriate ASSET type (vf.VectorFunction if output size is >1, vf.ScalarFunction if output size =1) and then defining and initializing the expression in the constructor. This method should only be preferred if you need to store the meta-data as part of the class or add additional methods to the object. Otherwise, this method is functionally identical to the one above.

class FuncClass(vf.VectorFunction):
        def __init__ (self,a,b,c):
                self.a =a
                self.b =b
                self.c =c

                x0,x1,x2 = Args(3).tolist()
                eq1 = x0 +a - x1
                eq2 = x2*b + x1*c

                super().__init__(vf.stack(eq1,eq2)) #Do not forget to call CTOR of Base!!

        def get_a(self):return self.a

func = FuncClass(1,2,3)

print(func([1,1,1]))  # prints [1,5]
print(func.get_a())   # prints 1

Function Composition

Now that we have a good understanding of the rules and style for defining single VectorFunctions, we can cover how to call them inside of other functions. For this final example let us tackle a concrete problem that occurs in astrodynamics: frame conversions. Specifically, we wish to write a function that takes the position and velocity of some object in Cartesian coordinates, as well as some other vector, and then transforms that vector into the RTN frame. The RTN basis vectors can be computed purely as a function of position and velocity, so let us first write a function that does just that.

def RTNBasis():

        R,V = Args(6).tolist([(0,3),(3,3)])

        Rhat = R.normalized()
        Nhat = R.cross(V).normalized()
        That = Nhat.cross(R).normalized()

        return vf.stack(Rhat,That,Nhat)

We can then write another function that takes position and velocity as well as the vector to be transformed. We then instantiate our previously defined function that computes basis vectors and then “call” it with the position and velocity arguments defined inside our new function. Calling the already instantiated function can be accomplished by passing other VectorFunctions, to the () call operator the same way we do for real number arguments. In this case, providing the contiguous segment of size 6 RV, is the most efficient way to the define the expression. However, if this were not the case, we could also use the other call signatures shown. We can provide two separate functions, in this case R and V, either as individual arguments or grouped together in a python list. These will be implicitly stacked using the same rules governing vf.stack and then forwarded to the function.

def RTNTransform():


        X = Args(9)

        RV,U = X.tolist([(0,6),(6,3)])

        R,V = X.tolist([(0,3),(3,3)])

        RTNBasisFunc = RTNBasis() # Instantiate function object


        RTNcoeffs = RTNBasisFunc(RV)  ### Call Function at new vectorfunction arguments

        RTNcoeffs = RTNBasisFunc(R,V) # Same effect as original but slower
        RTNcoeffs = RTNBasisFunc(vf.stack(R,V)) # Does Exactly the same thing as above calls stack on R,V explicitly

        RTNcoeffs = RTNBasisFunc([R,V]) # Same effect as original but slower
        RTNcoeffs = RTNBasisFunc(vf.stack([R,V])) # Does Exactly the same thing as above calls stack on [R,V] explicitly



        RTNmat = vf.RowMatrix(RTNcoeffs,3,3)

        U_RTN = RTNmat*U

        return U_RTN

Repeated Sub Expressions

Being a functional programming language, it is important to note that an ASSET expression is evaluated everywhere it appears in a statement. There is no notion of assigning it to a temporary variable and then reusing it later without recalculating it. For example, in the following code, just because we bind the complicated expression to the name expensive, the function answer will still require actually evaluating expensive three times.

R,V = X.tolist([(0,3),(3,3)])

expensive = 1.0/(R.normalized().cross(V.normalized_power3()).dot(R+V.cross(R).normalized()))**3.14

answer = R+ vf.stack(expensive,expensive+1,expensive)

In the vast majority of cases you should not worry about the cost of reevaluating subexpressions, as the run time hit is marginal. There is, however, one way to explicitly reduce the cost of expensive repeated sub expressions, should you need to. You can do this by writing a second function where the subexpression appears linearly as additional arguments or segments and then using the call operator to compose this new function and the original arguments and subexpression together. For example, the following code will produce the same output as above while only ever evaluating expensive once.

R,V = X.tolist([(0,3),(3,3)])

expensive = 1.0/(R.normalized().cross(V.normalized_power3()).dot(R+V.cross(R).normalized()))**3.14


## New args for defining function of only R and expensive
R_temp, expensive_tmp = Args(4).tolist([(0,3),(3,1)])

answer_tmp = R_temp+ vf.stack(expensive_tmp,expensive_tmp+1,expensive_tmp)

answer = answer_tmp([R,expensive])

Tabular Data and Interpolation

We also have support for interpreting tabular data as a VectorFunction using differentiable interpolation table objects.

1-D Interpolation

Interpolation of vector or scalar data with one input dimension may be accomplished using vf.InterpTable1D. This class is designed to be constructed and behave similarly to scipy’s interp1D class. To construct a table for interpolating vector data, we pass a list of sorted coordinate values along with an array whose rows or columns are the vectors of values at each coordinate. You may also pass in the data as a list of numpy arrays or lists which include the coordinate element. You have the option to choose between linear or cubic interpolation using the kind parameter in the constructor. Cubic interpolation is recommended to ensure that the function is twice differentiable.

ts = np.linspace(0,2*np.pi,1000)

VecDat = np.array([ [np.sin(t),np.cos(t)] for t in ts])

kind = 'cubic' # or 'linear'

Tab = vf.InterpTable1D(ts,VecDat,axis=0,kind=kind)
print(Tab(np.pi/2.0)) #prints [1,.0]

# Or if data is transposed
Tab = vf.InterpTable1D(ts,VecDat.T,axis=1,kind=kind)
print(Tab(np.pi/2.0)) #prints [1,.0]

# Or if data is a list of arrays or lists with time included as one the elements
VecList = [ [np.sin(t), np.cos(t), t] for t in ts]

Tab = vf.InterpTable1D(VecList,tvar=2,kind=kind)
print(Tab(np.pi/2.0)) #prints [1,.0]

To construct a table for interpolating scalar data, you may just pass in the list of coordinates along with a 1-D numpy array or python list of the values of the function at each point.

ScalDat = [np.sin(t) for t in ts]
STab =vf.InterpTable1D(ts,ScalDat,kind=kind)
print(STab(np.pi/2.0)) # prints [1.0]

The outputs of vf.InterpTable1D are only well defined within the domain of the supplied coordinates. By default attempts to interpolate outside of the domain will result in inaccurate extrapolation and print a warning the screen. This can be disabled using the .WarnOutOfBounds field of the object. Additionally, you may specify that you would like an exception to be thrown when this occurs by setting .ThrowOutOfBounds to True.

Tab.WarnOutOfBounds=True   # By default
print(Tab(-.00001))        # prints [-1.0e-5,1] and a warning
Tab.ThrowOutOfBounds=True
#print(Tab(-.00001))       # throws an exception

Once you have constructed a table object, it can be composed with other ASSET VectorFunctions or ScalarFunctions by providing a ScalarFunction argument to the call operator.

x,V,t = Args(4).tolist([(0,1),(1,2),(3,1)])

f1 = STab(t) + x  # STab(t) is an asset scalar function
f2 = Tab(t) + V   # Tab(t) is an asset vector function

2-D Interpolation

Similarly, you can also interpret scalar data defined on a 2-D rectilinear grid of coordinates as an ASSET ScalarFunction using the vf.InterpTable2D class. The class may be constructed by supplying the grid coordinates as either python lists or numpy arrays along with the function values formatted like a numpy meshgrid. The function values MUST be formatted like an xy indexed meshgrid, which is the default behavior for np.meshgrid in 2 dimensions

Note

The function values MUST be formatted like an xy indexed meshgrid.If you are using np.meshgrid to generate the function values, take note that it assumes xy indexing by default. If you index the values wrong and your coordinate dimensions are the same size, we will not be able to detect an error.

nx =500
ny =800
lim = 2*np.pi

xs = np.linspace(-lim,lim,nx)
ys = np.linspace(-lim,lim,ny)

def f(x,y):return np.sin(x)*np.cos(y)
X, Y = np.meshgrid(xs, ys)
Z    = f(X,Y)             #Scalar data defined on 2-D meshgrid

kind = 'cubic' # or 'linear'

Tab2D = vf.InterpTable2D(xs,ys,Z,kind=kind)

print(Tab2D(np.pi/2,0))  #prints 1.0


Tab2D.WarnOutOfBounds=True   # By default
print(Tab2D(-6.3,0))        # prints a warning
Tab2D.ThrowOutOfBounds=True
#print(Tab2D(-6.3,0))       # throws exception

Once constructed, vf.InterpTable2D can be converted into an ASSET ScalarFunction by supplying the x and y coordinates to the table’s call operator as a singe VectorFunction or two separate ScalarFunctions .

xy,c= Args(3).tolist([(0,2),(2,1)])
x,y = xy.tolist()

# Use it as scalar function inside a statement
Tab2sf = Tab2D(xy)
Tab2sf = Tab2D(x,y)             # Or
Tab2sf = Tab2D(vf.stack([x,y])) # Or

Func = Tab2sf + c   # Use it as a normal scalar function

print(Func([np.pi/2,0,1.0]))  # prints [2.0]

3-D Interpolation

As of version 0.1.0, we also support interpolating scalar data defined on 3 dimensional rectilinear grids with the vf.InterpTable3D class. To construct it, we supply the coordinate values as three python lists/numpy vectors along with the function values.

Note

The function values MUST be formatted like an ij indexed meshgrid. If you are using np.meshgrid to generate the function values, take note that it assumes xy indexing by default, and this should be changed to ij.If you index the values wrong and your coordinate dimensions are the same size, we will not be able to detect an error.

As with the 1D and 2D tables, you can also select between linear or cubic interpolation. Cubic interpolation should be preferred for anything that will end up in an optimizer. Additionally, for cubic interpolation you can specify that you want to pre-calculate and cache all possible values of the interpolation coefficients in each voxel of the domain. This requires the up front calculation of (nx-1)*(ny-1)*(nz-1) matrix vector products of size (64x64)x(64x1) and 8 times the storage of the non-cached algorithm. However, interpolation will be 10-20 times faster. So use this if the dimensions are small and the interpolator will be called a lot.

def f(x,y,z):return np.cos(x)*np.cos(y)*np.cos(z)

nx = 100
ny = 100
nz = 100

xlim = np.pi
ylim = np.pi
zlim = np.pi

xs = np.linspace(-xlim,xlim,nx)
ys = np.linspace(-ylim,ylim,ny)
zs = np.linspace(-zlim,zlim,nz)

X,Y,Z = np.meshgrid(xs, ys,zs,indexing = 'ij')
Fs    = f(X,Y,Z)    #Scalar data defined on 3-D meshgrid in ij format!!!

kind = 'cubic' # or 'linear', defaults to 'cubic'
cache = False # defaults to False
#cache = True # will precalculate and cache all interpolation coeffs

Tab3D = vf.InterpTable3D(xs,ys,zs,Fs,kind=kind,cache=cache)

print(Tab3D(0,0,0))  #prints 1.0

Tab3D.WarnOutOfBounds=True   # By default
print(Tab3D(-10,0,0))        # prints a warning
print(Tab3D(0,-10,0))        # prints a warning
print(Tab3D(0,0,-10))        # prints a warning

Tab3D.ThrowOutOfBounds=True
#print(Tab3D(-10,0,0))       # throws exception

4-D Interpolation

As of version 0.4.0, we also support interpolating scalar data defined on 4 dimensional rectilinear grids with the vf.InterpTable4D class. To construct it, we supply the coordinate values as four python lists/numpy vectors along with the function values.

Note

The function values MUST be formatted like an ij indexed meshgrid. If you are using np.meshgrid to generate the function values, take note that it assumes xy indexing by default, and this should be changed to ij.If you index the values wrong and your coordinate dimensions are the same size, we will not be able to detect an error.

As with the previous tables, you can select between linear or cubic interpolation. Cubic interpolation should be preferred for anything that will end up in an optimizer. Additionally, for cubic interpolation you can specify that you want to pre-calculate and cache all possible values of the interpolation coefficients in each voxel of the domain. This requires the up front calculation of (nx-1)*(ny-1)*(nz-1)*(nw-1) matrix vector products of size (256x256)x(256x1) and 16 times the storage of the non-cached algorithm. However, interpolation will be 5-6 times faster. So use this if the dimensions are small and the interpolator will be called a lot.

def f(x,y,z,w):return np.cos(x)*np.cos(y)*np.cos(z)*np.cos(w)

nx = 50
ny = 50
nz = 50
nw = 50

xlim = np.pi
ylim = np.pi
zlim = np.pi
wlim = np.pi

xs = np.linspace(-xlim,xlim,nx)
ys = np.linspace(-ylim,ylim,ny)
zs = np.linspace(-zlim,zlim,nz)
ws = np.linspace(-zlim,zlim,nw)


X,Y,Z,W = np.meshgrid(xs,ys,zs,ws,indexing = 'ij')
Fs    = f(X,Y,Z,W)    #Scalar data defined on 4-D meshgrid in ij format!!!

kind = 'cubic' # or 'linear', defaults to 'cubic'
cache = False # defaults to False
#cache = True # will precalculate and cache all interpolation coeffs

Tab4D = vf.InterpTable4D(xs,ys,zs,ws,Fs,kind=kind,cache=cache)

print(Tab4D(0,0,0,0))  #prints 1.0

Tab4D.WarnOutOfBounds=True   # By default
print(Tab4D(-10,0,0,0))        # prints a warning
print(Tab4D(0,-10,0,0))        # prints a warning
print(Tab4D(0,0,-10,0))        # prints a warning
print(Tab4D(0,0,0,-10))        # prints a warning

Tab4D.ThrowOutOfBounds=True
#print(Tab4D(-10,0,0,0))       # throws exception

Once constructed, vf.InterpTable4D can be converted into an ASSET ScalarFunction by supplying the x,y,z, and w coordinates to the table’s call operator as a singe VectorFunction or four separate ScalarFunctions.

xyzw,c= Args(5).tolist([(0,4),(4,1)])
x,y,z,w = xyzw.tolist()

# Use it as scalar function inside a statement
Tab4sf = Tab4D(xyzw)
Tab4sf = Tab4D(x,y,z,w)             # Or
Tab4sf = Tab4D(vf.stack([x,y,z,w])) # Or

Func = Tab4sf + c

print(Func([0,0,0,0,1]))  # prints [2.0]

Note on Size of VectorFunctions

The VectorFunction type system has been designed to have good performance for evaluating the value and derivatives of dense VectorFunctions with a small number of arguments (<50). It will work for larger expressions, but performance will begin to degrade considerably. This may seem strange since it is ostensibly designed to be used to define constraints and objective inside of large sparse non-linear programs. However, in our experience these problems are almost never composed of single monolithic functions, and can generally be decomposed into smaller dense functions that only take partial subsets of the problem variables. In that case, we can define our functions in terms of only the arguments they take, and then under the hood, ASSET will ensure that the inputs and outputs are gathered and scattered to the correct locations inside the larger problem. The specifics of how this works will be discussed in later sections.

X = Args(1000) # Legal but Bad

Binding Raw Python Functions (DON’T DO THIS)

You also have the option, should you need to, to bind raw python functions as ASSET VectorFunctions and ScalarFunctions. This can be accomplished using the vf.PyVectorFunction and vf.PyScalarFunction types as shown below. The function must have a signature accepting as the first argument a 1 dimensional numpy array of input arguments (named X in this case) and returning a numpy array. Additional parameters on which the implementation depends (these are not mathematical input variables) may be included as additional arguments. You must also explicitly state the input and output (if not scalar) sizes of the function. The function jacobian and hessian will be computed with finite differences using the user specified jacobian and hessian step sizes.

# A vector function
def VFunc(X,a,b):
        return np.array([a*X[0]**2,X[1]*b])

InputSize = 2
OutputSize =2

PyVfunc =vf.PyVectorFunction(InputSize,OutputSize,VFunc,Jstepsize=1.0e-6,Hstepsize=1.0e-4,args = (3,7)) ## a and b will be 2

print(PyVfunc([2,2]))  # prints [12,14]


# A scalar function
def SFunc(X,a,b,c):
        return np.array([a*X[0]**2 + X[1]*b + c]) # output is 1x1 array

InputSize = 2

PySfunc =vf.PyScalarFunction(InputSize,SFunc,Jstepsize=1.0e-6,Hstepsize=1.0e-5,args = (1,2,3))

You should be warned that extensive use of these objects inside of the optimizer or an ODE will result in VERY slow and non-parrallelizable code with inexact derivatives. If you find yourself in a situation where you don’t think you can write an expression without using vf.PyVectorFunction or vf.PyScalarFunction, please submit an issue on GitHub. We will happily give suggestions on how you might be able to accomplish your task with the standard VectorFunctions. If it’s truly not possible, we will consider adding the missing expression to the core library in a future release.