Vector-valued functions¶

In the case of a vector-valued function (i.e. functions which return an iterable object), the derive class will automatically return a list of jet objects upon evaluation, one for each component. Although derive can handle any composition of (vector-valued) functions in this way, it may be of use to combine multi-dimensional jet output in case they were produced in different processes. Support for the handling of such cases is given in the dedicated njet.extras module. In particular, there exist a general_fa_di_bruno routine and the cderive class which we shall describe in the following.

Of course, everything what will be discussed here can also be applied to the special 1D case. However, some of the routines in njet.extras (in particular those of the cderive class) will take advantage of jets carrying NumPy arrays and so this should be taken into account. A jet representing a single-valued function may be converted to a jet carrying a 1D NumPy array by means of the njet.extras.tile routine.

Faa di Bruno’s formula¶

Faa di Bruno’s formula describes the relation between the higher-order derivatives of a composition function and those of its constituents. In the one-dimensional case this relation reads: 1

\[\frac{d^n}{dt^n} f(g(t)) = \sum_{k = 1}^n f^{(k)}(g(t)) \cdot B_{n, k} (g'(t), g''(t), ..., g^{(n - k + 1)}(t)) ,\]

where \(B_{n, k}\) denotes the exponential incomplete Bell-polynomials, which can be defined with a generating function:

\[\exp \left(u \sum_{j = 1}^\infty x_j \frac{t^j}{j!} \right) =: \sum_{n \geq k \geq 0} B_{n, k} (x_1, ..., x_{n - k + 1}) \frac{t^n}{n!} u^k .\]

These polynomials obey a recurrence relation by which internally njet computes the composition of two jets (by means of the @-operator):

\[\begin{split}B_{n, k} = \sum_{i = 1}^{n - k + 1} \left(\begin{array}{c} n - 1 \\ i - 1 \end{array}\right) x_i B_{n - i, k - 1} .\end{split}\]

A similar – but more complicated – expression can be derived in the multi-dimensional case, i.e. where \(f\) takes more than one argument and \(g\) is vector-valued. 2 Without loss of generality let \(f\) be single-valued, so let \(f \colon \mathbb{K}^m \to \mathbb{K}\) and \(g \colon \mathbb{K}^l \to \mathbb{K}^m\) be two (sufficiently often) differentiable functions. Since the higher-order partial derivatives of the map \(f \circ g \; \colon \mathbb{K}^l \to \mathbb{K}\) commute, the respective symmetric multilinear map \((f \circ g)^{(n)} \colon (\mathbb{K}^l)^{\otimes n} \to \mathbb{K}\) is characterized by its action on the ‘diagonal’ \(z^{\otimes n}\), \(z \in \mathbb{K}^l\) as

\[\begin{split}\frac{(f \circ g)^{(n)}(z^{\otimes n})}{n!} = \sum_{\substack{n_1 + ... + n_r = n\\r \geq 1, n_j \geq 1}} \frac{(f^{(r)} \circ g)}{r!} \left(\frac{g^{(n_1)} (z^{\otimes n_1})}{n_1!} \otimes ... \otimes \frac{g^{(n_r)} (z^{\otimes n_r})}{n_r!} \right) ,\end{split}\]

by means of the polarization formula 3

\[z_1 \otimes z_2 \otimes ... \otimes z_k = \frac{1}{2^k k!} \sum_{\epsilon_j = \pm 1} \epsilon_1 \cdot ... \cdot \epsilon_k (\epsilon_1 z_1 + ... + \epsilon_k z_k)^{\otimes k} .\]

Regarding the polynomials in the njet module, this generalized Faa di Bruno formula can therefore be implemented by summing up those jet-components which correspond to a partition of the integer \(n\) according to the above equation.

Implementation¶

The generalized Faa di Bruno formula can be imported with

from njet.extras import general_faa_di_bruno

Define two vector-valued functions and compute their jet evaluation at a specific point z:

f = lambda *x: [x[0]**2 - (x[1] + 10 - 0.5*1j)**(-1) + x[2]*x[0], 1 + x[0]**2*x[1] + x[1]**3, x[2]]
g = lambda *x: [(x[0] + x[1]*x[2]*(0.9*1j + 0.56) + 1)**(-3), 2*x[1]*4j + 5, x[0]*x[1]**6]

from njet import derive

z = [0.4, 1.67 + 0.01*1j, -3.5 + 2.1*1j]

dg = derive(g, order=2, n_args=3)
df = derive(f, order=2, n_args=3)

evg = dg.eval(*z)
evf = df.eval(*g(*z))

The two lists evg and evf contain the jet-evaluations of the two functions at the points z and g(*z) for every component, and hence their Taylor-coefficients. For example, we can get the Taylor-coefficients of f at g(*z) for the first component (the n_args argument here is just to tell njet how to present the 0th-order and not immediately necessary):

  evf[0].taylor_coefficients(n_args=3)
> {(0, 0, 0): (-0.032316721543419864+0.07248263230025642j),
   (0, 0, 1): (0.0008719493399318281+0.0045037973966777466j),
   (1, 0, 0): (8.673861926635063+0.3207111009469555j),
   (0, 1, 0): (0.00038016081672203066-0.002549222116121679j),
   (1, 0, 1): (1+0j),
   (2, 0, 0): (2+0j),
   (0, 2, 0): (0.00013975244997413352+0.00022126191190187668j)}

Of course, we could have obtained this result directly by calling df at g(*z):

  df(*g(*z))[0]
> {(0, 0, 0): (-0.032316721543419864+0.07248263230025642j),
   (0, 0, 1): (0.0008719493399318281+0.0045037973966777466j),
   (1, 0, 0): (8.673861926635063+0.3207111009469555j),
   (0, 1, 0): (0.00038016081672203066-0.002549222116121679j),
   (1, 0, 1): (1+0j),
   (2, 0, 0): (2+0j),
   (0, 2, 0): (0.00013975244997413352+0.00022126191190187668j)}

Here we are interested in the Taylor-coefficients of the composition function \(f \circ g\). In the conventional approach we would have to derive the composition function:

  dfg = derive(lambda *x: f(*g(*x)), order=2, n_args=3)
  ref = dfg(*z)
  ref[0]
> {(0, 0, 0): (-0.032316721543419864+0.07248263230025642j),
   (0, 0, 1): (-0.009661433404866623+0.03378161852355409j),
   (0, 1, 0): (0.020635651416517554+0.06139363375476808j),
   (1, 0, 0): (0.028801731735594405+0.11295869722633461j),
   (1, 0, 1): (-0.01697989169984845+0.10671633370097759j),
   (1, 1, 0): (0.0058767806547268325+0.16073332310475294j),
   (0, 0, 2): (-0.02701945022246847+0.031297336355818564j),
   (0, 2, 0): (-0.01332575692398203+0.04362074226492914j),
   (0, 1, 1): (-0.022839564957718664+0.04217330624559906j),
   (2, 0, 0): (0.07990334788570935+0.07626091978109452j)}

However, making use of the general Faa di Bruno formula, we can deduce the same result by combining the previously computed multi-dimensional jet-evaluations evg and evf:

  gfb = general_faa_di_bruno(evf, evg)
  gfb[0].taylor_coefficients(n_args=3)
> {(0, 0, 0): (-0.032316721543419864+0.07248263230025642j),
   (0, 0, 1): (-0.009661433404866623+0.033781618523554095j),
   (1, 0, 0): (0.02880173173559441+0.11295869722633461j),
   (0, 1, 0): (0.02063565141651755+0.06139363375476808j),
   (1, 0, 1): (-0.016979891699848447+0.1067163337009776j),
   (1, 1, 0): (0.005876780654726825+0.16073332310475294j),
   (0, 2, 0): (-0.013325756923981996+0.04362074226492907j),
   (0, 0, 2): (-0.02701945022246847+0.031297336355818564j),
   (0, 1, 1): (-0.02283956495771866+0.04217330624559905j),
   (2, 0, 0): (0.07990334788570935+0.07626091978109452j)}

In this way it is possible to calculate and combine intermediate steps of a chain, without the need of passing the jets through all of the members of the chain in a single large calculation. This can become effective in particular if there are repetitions of functions in the chain.

Function chains¶

In the case that a chain of functions of the form \(f_N \circ f_{N - 1} ... \circ f_1\) needs to be differentiated, and there are repetitions of \(f_k\)’s in the chain, the generalized Faa di Bruno formula may help in reducing the amount of calculations required.

For example, if we want to compute a chain of \(2^N\) repetitions of the same function \(f\), then we can first compute the composition \(f_2 := f \circ f\), then the compositon \(f_4 := f_2 \circ f_2\) and so on, until we are finished after \(N\) steps.

Moreover, if a chain indeed admits repetitions, a single point will pass through a specific function several times while it is transported from \(f_1\) to \(f_N\). We can take advantage of this by calculating the derivatives at these ‘intersecting’ points for each unique element of the chain in parallel, using NumPy.

To conveniently manage chains of functions, the dedicated class cderive (‘chain-derive’) has been implemented. This class takes a series of functions, defining the unique functions in the chain, an order parameter, defining the number of derivatives to be computed and an optional ordering parameter, which is a list of integers and defines the chain itself. It can be imported with the following line:

from njet.extras import cderive

As an example, consider an alternation of a rotation and some polynomial perturbation \(M = 15\) times:

from njet.functions import sin, cos

g = 1
per = lambda *z, **kwargs: [z[0], z[1] - g*z[1]**2]
rot = lambda *z, alpha=0: [cos(alpha)*z[0] - sin(alpha)*z[1], sin(alpha)*z[0] + cos(alpha)*z[1]]

M = 15
ordering = [0, 1]*M
drp = cderive(rot, per, order=2, ordering=ordering, n_args=2)

To compute the derivatives of this chain, we simply call the cderive object at the point of interest. Hereby we can pass any additional keyworded arguments to the respective functions (currently every function in the chain will be called with the same set of keyworded parameters, that’s why we have provided per with an **kwargs parameter). For example:

  z0 = [0.2, 0.1]
  alpha = 1.22
  drp(*z0, alpha=alpha)
> ({(0, 0): (0.18073767467258015+0j),
    (1, 0): (0.3704449136384162+0j),
    (0, 1): (0.3778714745973873+0j),
    (0, 2): (-1.2481054532394522+0j),
    (2, 0): (-2.5881230191830347+0j),
    (1, 1): (-0.74624063677712+0j)},
   {(0, 0): (-0.012836565325653162+0j),
    (1, 0): (-0.4426107548482863+0j),
    (0, 1): (0.8026126813486407+0j),
    (0, 2): (-0.49946504768391176+0j),
    (2, 0): (0.7197445629662029+0j),
    (1, 1): (-0.3181227776491683+0j)})

The cderive object drp in the above example contains two unique elements, which are instances of the derive class. These unique elements can be accessed in the .dfunctions field:

  drp.dfunctions
> [<njet.ad.derive at 0x7f6a4c09db40>, <njet.ad.derive at 0x7f6a4c09f250>]

If one is interested in obtaining the function at position k in the chain, then one can type drp[k], e.g.

  drp[13]
> <njet.ad.derive at 0x7f6a4c09f250>

In fact, one can iterate over the drp object as with ordinary lists. The ordering of the chain is stored in drp.ordering. If one wishes to change the ordering on the same object, there is a dedicated routine for it, drp.set_ordering, because with a change in the ordering also other internal changes will become necessary. In particular, any previously computed jet evaluations may become invalid, since points passing through the newly ordered chain will ‘hit’ the unique functions at different places. Access to these jet-evaluations (the derivatives) are conveniently done by the .jev function. For example, the Taylor-coefficients of the first component in the current jet evaluation at z0 in the above chain at position 4 can be obtained as follows:

  drp.jev(4)[0].taylor_coefficients(n_args=2)
> {(0, 0): (-0.0911100230191827+0j),
   (1, 0): 0.34364574631604705,
   (0, 1): -0.9390993563190676}

In some circumstances, however, previously computed jet-evaluations remain valid. For example, this is the case if we want to extract a subchain of a given chain, defined by a slice of neighbouring indices. We can extract subchains from the original chain in the same manner as it can be done with lists. For example:

  mysubchain = drp[1:13]
  mysubchain.jev(3)[0].taylor_coefficients(n_args=2)
> {(0, 0): (-0.0911100230191827+0j),
   (1, 0): 0.34364574631604705,
   (0, 1): -0.9390993563190676}

Notice that the index of the same data is now at 3, because the sub-chain starts at index 1 of the original chain.

Another scenario where jet-evaluation data remains valid is by merging patterns occuring in the ordering by means of the general Faa di Bruno formula. A pattern of functions in a chain can be merged by the .merge command, which will return a new cderive object in which the pattern has been merged. For example, if we want to merge the first occurence of the pattern (1, 0, 1) in the chain drp we could do the following:

drpm = drp.merge(pattern=(1, 0, 1), positions=[1])

The new chain drpm now has len(drpm) = 28, since we have merged 3 elements of the original chain of length 30 and added 1 new ‘merged’ element. Moreover, the new chain still maintains any previously computed jet-evaluations (although the ones of the merged functions will vanish internally).

The .merge command requires that jet-evaluation data has been computed in advance. This can be done with the .eval command (in fact, .eval is always called whenever the cderive class is taken at a specific point, but .eval omits the calculation of the Taylor-coefficients, which is not always necessary). The .eval command can preferably be invoked with the compose=False option, to prevent the automatic composition calculation. After merging, the jet-evaluation of the entire chain can be combined by means of successive ‘Faa di Bruno’ operations, using the .compose routine. In scenarios where chains admit repeated function patterns, these three steps may drastically improve performance in comparison to the conventional derive approach.

Coming back to our example, we confirm that the results before and after merging are the same:

  c1 = drp.compose()
  c2 = drpm.compose()

  from njet import taylor_coefficients

  taylor_coefficients(c1, n_args=2)
> ({(0, 0): (0.18073767467258015+0j),
    (1, 0): (0.3704449136384162+0j),
    (0, 1): (0.3778714745973873+0j),
    (0, 2): (-1.2481054532394522+0j),
    (2, 0): (-2.5881230191830347+0j),
    (1, 1): (-0.74624063677712+0j)},
   {(0, 0): (-0.012836565325653162+0j),
    (1, 0): (-0.4426107548482863+0j),
    (0, 1): (0.8026126813486407+0j),
    (0, 2): (-0.49946504768391176+0j),
    (2, 0): (0.7197445629662029+0j),
    (1, 1): (-0.3181227776491683+0j)})

  taylor_coefficients(c2, n_args=2)
> ({(0, 0): (0.18073767467258015+0j),
    (1, 0): (0.3704449136384162+0j),
    (0, 1): (0.3778714745973873+0j),
    (0, 2): (-1.2481054532394522+0j),
    (2, 0): (-2.5881230191830347+0j),
    (1, 1): (-0.7462406367771199+0j)},
   {(0, 0): (-0.012836565325653162+0j),
    (1, 0): (-0.4426107548482863+0j),
    (0, 1): (0.8026126813486407+0j),
    (0, 2): (-0.49946504768391176+0j),
    (2, 0): (0.7197445629662029+0j),
    (1, 1): (-0.3181227776491683+0j)})

In case one has a periodic system, where the chain is traversed again and again, one might be interested in the derivatives along the chain for every possible cycle. A cycle of start index \(k\) is hereby understood as a chain of \(N\) functions of the form \(f_{k - 1} \circ ... \circ f_1 \circ f_N \circ f_{N - 1} \circ ... \circ f_k\). Such a calculation would have to be done for every \(k\), but it can also be performed in parallel using jets carrying NumPy entries. For this purpose there exist a dedicated routine: .cycle. This routine is called with the same input parameters as the ordinary cderive function. Returning to our example this would read:

cyc = drp.cycle(*z0, alpha=alpha)

The object cyc is a list of length len(drp) which contains the jet-evaluation results for every cycle. This includes our previous result:

  taylor_coefficients(cyc[0], n_args=2)
> ({(0, 0): (0.18073767467258015+0j),
    (1, 0): (0.3704449136384162+0j),
    (0, 1): (0.3778714745973873+0j),
    (0, 2): (-1.2481054532394522+0j),
    (2, 0): (-2.5881230191830347+0j),
    (1, 1): (-0.74624063677712+0j)},
   {(0, 0): (-0.012836565325653162+0j),
    (1, 0): (-0.4426107548482863+0j),
    (0, 1): (0.8026126813486407+0j),
    (0, 2): (-0.49946504768391176+0j),
    (2, 0): (0.7197445629662029+0j),
    (1, 1): (-0.3181227776491683+0j)})

as well as any other result around the cycle. For example, at the next position we check:

  ordering2 = [1] + [0, 1]*(M - 1) + [0]
  drp2 = cderive(rot, per, order=2, ordering=ordering2, n_args=2)

  taylor_coefficients(cyc[1], n_args=2) == drp2(*rot(*z0, alpha=alpha), alpha=alpha)
> True

In case that the underlying functions take parameters, it is important to remember that these parameters must be included in the call of the .cycle routine. In case some results look weird, there is a warn switch which might provide some clues.

Similar to other routines, also the .cycle routine can handle multi-dimensional NumPy arrays:

  z1 = np.array([[0.02, -0.056, z0[0]], [0.0031, 0.0118, z0[1]]])
  cyc = drp.cycle(*z1, alpha=alpha)
  taylor_coefficients(cyc[0], n_args=2)
> ({(0, 0): array([ 0.01863169+0.j, -0.04082189+0.j,  0.18073767+0.j]),
    (1, 0): array([0.84531442+0.j, 0.80706613+0.j, 0.37044491+0.j]),
    (0, 1): array([0.5277127 +0.j, 0.51823247+0.j, 0.37787147+0.j]),
    (0, 2): array([ 0.57344527+0.j,  1.20526504+0.j, -1.24810545+0.j]),
    (2, 0): array([-0.72251201+0.j,  1.81211122+0.j, -2.58812302+0.j]),
    (1, 1): array([ 0.10110713+0.j,  0.33841931+0.j, -0.74624064+0.j])},
   {(0, 0): array([-0.00775766+0.j,  0.0386481 +0.j, -0.01283657+0.j]),
    (1, 0): array([-0.51833409+0.j, -0.49550034+0.j, -0.44261075+0.j]),
    (0, 1): array([0.85906428+0.j, 0.81893139+0.j, 0.80261268+0.j]),
    (0, 2): array([ 0.55558523+0.j, -0.11809903+0.j, -0.49946505+0.j]),
    (2, 0): array([ 0.25984599+0.j, -0.73940164+0.j,  0.71974456+0.j]),
    (1, 1): array([ 0.13269792+0.j,  0.94135934+0.j, -0.31812278+0.j])})

We can see in the above example that the last entry in each coefficient agrees with our previous result, as expected.

The .cycle routine has the option to return a cderive object instead of a list of jet-evaluations. This cderive object hereby contains the (multi-dimensional) NumPy arrays, representing the different traces which are to be calculated in parallel (see the docstring in .cycle for an illustration of these traces), and so it will have a different length of 2*len(drp) - 1 with ordering (drp.ordering*2)[:-1]. If the traces are combined, they will yield a jet-evaluation result of NumPy arrays, where every entry corresponds to the list above. To demonstrate this we consider position 5 in our example-chain with three different (two-dimensional) input points z1. In the 0-th component we have:

  pos = 5
  cyc[pos][0].taylor_coefficients(n_args=2)
> {(0, 0): array([ 0.01178258+0.j, -0.01398243+0.j,  0.12224696+0.j]),
   (1, 0): array([0.84452107+0.j, 0.83519304+0.j, 0.49549387+0.j]),
   (0, 1): array([0.52872394+0.j, 0.50442847+0.j, 0.50766504+0.j]),
   (0, 2): array([-0.34437699+0.j, -0.4025439 +0.j, -1.51963413+0.j]),
   (2, 0): array([-0.62792391+0.j,  1.18315592+0.j, -2.73434645+0.j]),
   (1, 1): array([ 0.21824765+0.j,  0.0447755 +0.j, -0.48232797+0.j])}

  cc = drp.cycle(*z1, alpha=alpha, outf=0)
  cc_m = cc.compose()
  cc_m[0][pos].taylor_coefficients(n_args=2)
> {(0, 0): array([ 0.01178258+0.j, -0.01398243+0.j,  0.12224696+0.j]),
   (1, 0): array([0.84452107+0.j, 0.83519304+0.j, 0.49549387+0.j]),
   (0, 1): array([0.52872394+0.j, 0.50442847+0.j, 0.50766504+0.j]),
   (0, 2): array([-0.34437699+0.j, -0.4025439 +0.j, -1.51963413+0.j]),
   (2, 0): array([-0.62792391+0.j,  1.18315592+0.j, -2.73434645+0.j]),
   (1, 1): array([ 0.21824765+0.j,  0.0447755 +0.j, -0.48232797+0.j])}

The option to return the cderive object cc has the feature that, before combining, we may take advantage of any possible pattern repetitions (of the cc chain) with merge command(s), e.g.:

  cc1 = cc.merge(pattern=(1, 0, 1), positions=[1])
  cc1_m[0][pos].taylor_coefficients(n_args=2)
> {(0, 0): array([ 0.01178258+0.j, -0.01398243+0.j,  0.12224696+0.j]),
   (1, 0): array([0.84452107+0.j, 0.83519304+0.j, 0.49549387+0.j]),
   (0, 1): array([0.52872394+0.j, 0.50442847+0.j, 0.50766504+0.j]),
   (0, 2): array([-0.34437699+0.j, -0.4025439 +0.j, -1.51963413+0.j]),
   (2, 0): array([-0.62792391+0.j,  1.18315592+0.j, -2.73434645+0.j]),
   (1, 1): array([ 0.21824765+0.j,  0.0447755 +0.j, -0.48232797+0.j])}

Notice that the pos variable can be used just fine even after merging (which will produce the chain cc1, having a different length than cc), because now the pos variable represents a component in the multi-dimensional NumPy array which is tracked through the chain, and not a position in the cc or cc1 chain.

Note

In case of larger chains it is recommended to delete any previously computed result of the .cycle module explicitly, to free up memory before .cycle is called again.

Warning

The .extras module is more experimental. Therefore, please check your results carefully and let me know if you encounter any problems or bugs.

Module synopsis¶

class njet.extras.cderive(*functions, order: int = 1, ordering=None, run_params={}, reset=True, **kwargs)¶

Class to handle the derivatives of a chain of vector-valued functions with repetitions. The given functions should be unique, while their repetition in the chain will be given by an optional ordering.

Parameters

functions (list) –
A list of callable(s) or ‘derive’ classes related to the unique vector-valued*) functions in the chain to be derived.

Alternatively, one can pass a list of ‘derive’ objects, containing vector-valued jet evaluation output. This second case can be used to avoid re-calculating the derivatives in circumstances where it is not required.

Attention: By default, the first function is assumed to be executed first, in contrast to the mathematical notation. This can be changed by passing an ‘ordering’ argument (see below).

*) This means that the functions must return iterables in any case.
order (int, optional) – The maximal order of the derivatives to be computed. If nothing specified, the first derivative(s) will be computed.
ordering (list, optional) – The order defining how the unique functions are arranged in the chain. Hereby the index j must refer to the function at position j in the sequence of functions.
run_params (dict, optional) – A dictionary containing the output of njet.extras._make_run_params. These parameters can be send to the routine to avoid internal re-calculation when the routine it is called repeatedly.
reset (boolean, optional) – Parameter given to self.set_ordering. If ‘True’ (default), erase any jet-evaluation data coming along with the input (the ._evaluation fields). To prevent this behavior, set this parameter to ‘False’. However, in this case one should ensure that correct data has been stored – in particular concerning the requested ordering of the chain.
**kwargs – Optional keyworded arguments passed to njet.ad.derive init.

build_tensor(*args, **kwargs)¶

compose(**kwargs)¶

Compose the current jet-evaluations in the chain (requires function evaluations in advance).

See njet.extras.compose for details; its output is stored in self._evaluation._chain, while the last element is returned.

cycle(*args, outf='default', **kwargs)¶

Cycle through the given chain: Compute the derivatives at each point, assuming a periodic structure of the entire chain.

The basic idea goes as follows: Let c1 –> c2 –> c3 –> c4 –> c5 denote the current chain (in this example consisting of 5 elements). Then this routine will extend the internal jet-evaluation data to numpy arrays of length 5 by suitable clones and identity operators, so that we can compute (compose) in parallel the following chain of length 2*5 - 1:

c1 —> c2 —> c3 —> c4 —> c5*

c2 —> c3 —> c4 —> c5 —> c1*

c3 —> c4 —> c5 —> c1 —> c2*

c4 —> c5 —> c1 —> c2 —> c3*

c5 —> c1 —> c2 —> c3 —> c4*

The final results will then be returned as the starred values in the above diagram.

Attention: For larger cycles it may become necessary to delete any previously computed output of this routine in advance, in order to free up memory.

Parameters

*point (float or iterable, optional) – The point at which to evaluate the derivative(s). It has to be provided if the chain has not yet stored any jet-evaluation data.
periodic (boolean, optional) – If true, assume that the final point equals the given point. In particular, the data taken in the calculation will be taken from the current jet-evaluations (if it exists).
outf (str, optional) – Output format; if ‘default’, compose the extended jet-evaluations and return a list of jet-evaluations, representing the derivatives along the chain. Otherwise, return a cderive object of length len(self)*2 - 1 for further processing. The object can be composed to yield the jet-evaluations along the current chain of len(self) in form of numpy arrays. This second option is intended to be used for performance improvements.
**kwargs – Optional keyworded arguments passed to a (possible) chain evaluation.

Returns

Depending on the ‘outf’ parameter the output is either 1) A list of jet-evaluations, where the k-th entry corresponds to the derivative of the chain from position k to position k + L, where L denotes the length of the chain (using a periodic chain structure). 2) A cderive class with the property mentioned above.

Return type

list or cderive

eval(*point, compose=True, **kwargs)¶

Evaluate the functions in the chain at the requested point.

Parameters

*point (single value or array-like) – The point at which the evaluation should be performed.
compose (boolean, optional) – If true (default), compose the evaluations after they have been computed and return the result, resembling the behaviour of the original njet.derive.eval routine. If false, only generate the jet-evaluation data (e.g. in preparation for a succeeding merge command).
**kwargs (dict, optional) –
Keyworded arguments passed to the underlying jet-functions.

One can pass the boolean parameter disable_tqdm to enable (if False) a progress bar.

Returns

The outcome of self.compose routine (containing the jet-evaluations for each component).

Return type

list

grad(*args, **kwargs)¶

hess(*args, **kwargs)¶

index(value)¶: Return the element number at a given position in the beamline chain (similar as .index for lists)

jetfunc(*point, **kwargs)¶

Pass a point through the chain.

Parameters

*point (single value or array-like) – The point at which the chain should be evaluated.
**kwargs (dict, optional) – Keyworded arguments passed to the underlying jet-functions.

Returns

The final value after the chain of functions has been traversed.

Return type

single value or array-like

jev(pos: int)¶

Convenience function to obtain jet-evaluation data at a specific position (such data can be produced by the ‘eval’ command).

Parameters: pos (integer) – The position within the chain of functions.
Returns: A list of jet-evaluations.
Return type: list

merge(pattern=(), positions=None, **kwargs)¶

Merge one or more sections in the current chain simultaneously.

If a pattern may occur on several places, an additional parameter ‘positions’ has to be provided, so that all patterns are non-overlapping.

Parameters

pattern (tuple or list, optional) – Tuple of integers which defines a subsequence in self.ordering. If nothings specified, the entire sequence will be used, and so this routine becomes very similar to self.compose (with the difference that a cderive object will be returned here).
positions (list, optional) –
List of integers which defines the start indices of the above pattern in self.ordering. If nothing specified, every occurence of ‘pattern’ in self.ordering will be used.

In any case, only the members of ‘pattern’ are merged internally, so the number of occurences of that pattern should not affect the computational cost.
**kwargs – Optional keyworded arguments passed to the individual derive class(es).

Returns

A cderive object which contains a sequence of ‘derive’ classes, representing a new chain of functions in which the selected pattern(s) have been merged.

Return type

cderive

reset()¶: Remove all jet evaluations (if they exist).

set_ordering(ordering=None, reset=True)¶

Set the ordering of the current chain. Also compute the number of passages expected through each element according to the ordering.

Parameters

ordering (list, optional) – The order defining how the unique functions are arranged in the chain. Hereby the index j must refer to the function at position j in the sequence of functions.
reset (boolean, optional) – Reset the jet-evaluations as well. This may become necessary because these evaluations depend strongly on the current ordering.

njet.extras.compose(*evals, run_params={}, **kwargs)¶

Compose derivatives of a chain of vector-valued jet-evaluations.

Parameters

evals – A series of (multi-dimensional) jet-evaluations. Hereby every jet-evaluation is represented by a list of n-jets, one for each individual component.
**kwargs (dict, optional) – One can pass the boolean parameter disable_tqdm to enable (if False) a progress bar.

Returns

A list of jet evaluations for each vector component of the chain, representing the Taylor-expansion of the chain.

Return type

list

njet.extras.general_faa_di_bruno(f, g, run_params={})¶

Faa di Bruno for vector-valued functions.

Let G: K^l -> K^m and F: K^m -> K^w be two vector-valued functions. The goal is to compute the higher-order derivatives of the composition F o G.

Assume that f = [f1, …, fw] and g = [g1, …, gm] represent the n-jet-collections of F and G, i.e. fk = n-jet(fk_1, …, fk_n) where fk_j represents the j-th derivative of the k-th component of F etc.

Note that this routine requires that fk (and gl etc.) are represented by jetpoly objects in their higher-order array entries (they store the various partial derivatives for the given orders). In particular, the ordinary derive.eval routine will produce such objects. This is the major difference to the one-dimensional Faa di Bruno formula in njet.jet.faa_di_bruno, which has no such requirement.

Parameters

f (list) – A list of n-jet objects, representing the derivatives of the function F at position G(z).
g (list) – A list of n-jet objects, representing the derivatives of the function G at a given position z.
run_params (dict, optional) – A dictionary containing the output of njet.extras._make_run_params. These parameters can be send to the routine to avoid internal re-calculation when the routine it is called repeatedly.

Returns

A list of jet objects, representing the higher-orders of the components of the compositional map F o G.

Return type

list

njet.extras.tile(jev, ncopies: int)¶

Construct a new jet in which its entries contain copies of itself.

Intended to work if the jet entries are numpy arrays. A jet with scalar values can be converted to a jet containing numpy arrays of shape (1,) by using ncopies=1.

Parameters

jev (jet) –
ncopies (int) – Number of desired copies.

Returns

A jet having arrays of length ncopies + 1

Return type

jet

1: https://en.wikipedia.org/wiki/Fa%C3%A0_di_Bruno%27s_formula
2: https://mathoverflow.net/questions/106323/faa-di-brunos-formula-for-vector-valued-functions
3: Note that the operator ‘\(\otimes\)’ can be considered as being commutative when used in an argument of a symmetric tensor.