5. Ordinary Differential Equations on GraphsGraphs
5.4. Rate Assignment Operator
Listing 5.4: Main function for diffusion in a tree using advanced integration methods.
1 // c a l c u l a t e l e n g t h o f s t a t e v e c t o r
We have seen in the previous section 5.3.2 that even in a simple model the creation of the mapping between node attributes and the elements of the rate/state vector can become complex swiftly. In a more general setting, this can be even quite challenging, especially if class hierarchies are involved. Therefore, the system should assist the user by creating this mapping automatically.
Consider the class diagram shown in figure 5.6a. It is assumed that all attributes Memory Allocation except x participate in integration. So for all of these attributes a location in the rate
and state vector must be assigned and corresponding memory must be allocated. Let nA, nB, nC, nD designate the number of instances of the respective types A,B,C, and D.
Let further νA = 1, νB = 0, νC = 2, νD = 1 indicate the number of attributes of these
Listing 5.5: Rate function for diffusion in a tree using advanced integration methods.
1 void g e t R a t e (double[ ] r a t e , double time , double[ ] s t a t e ) {
2 // z e r o o u t p u t a r r a y
3 j a v a . u t i l . A r r a y s . f i l l ( r a t e , 0 ) ;
4
5 // copy s t a t e y t o g r a p h
6 [ c : C : :> c [ c a r b o n ] = y [ c [ i n d e x ] + 0 ] ; ]
7 [ c : C : :> c [ l e n g t h ] = y [ c [ i n d e x ] + 1 ] ; ]
8
9 // c a l c u l a t e r a t e v e c t o r
10 [
11 // a p p l y p r o d u c t i o n t o A n o d e s
12 c :A : :> r a t e [ c [ i n d e x ] ] += PROD;
13
14 // a p p l y c o n s u m p t i o n and c o n v e r t t o g r o w t h
15 c : B : :> {
16 double r = CONS ∗ c [ c a r b o n ] ;
17 r a t e [ c [ i n d e x ] + 0 ] −= r ;
18 r a t e [ c [ i n d e x ] + 1 ] += γ ∗ r ;
19 }
20
21 // p e r f o r m d i f f u s i o n b e t w e e n n o d e s
22 ca : C (−−>)+ : ( cb : C) : :> {
23 double r = D ∗ ( ca [ c a r b o n ] − cb [ c a r b o n ] ) ;
24 r a t e [ ca [ i n d e x ] ] −= r ;
25 r a t e [ cb [ i n d e x ] ] += r ;
26 }
27 ]
28 }
types. Then the total lengthnof the state vector can be computed as:
n=nA·νA+nB·(νA+νB) +nC·(νA+νB+νC) +nD·(νA+νD).
It follows that for each type the number of attributes it declares together with those inherited from its supertypes must be counted. With respect to integration, a reduced class hierarchy is obtained by omitting attributes and classes that do not participate in integration (see figure 5.6b).
An alternative way to compute n can be realized if tA = nA+nB+nC +nD, tB = nB+nC, tC =nC, tD =nD denote the total number of instances that can be cast to a certain type. Then the size of the rate and state vector is evaluated as:
n=tA·νA+tB·νB+tC·νC +tD·νD. (5.2) This latter variant is beneficial as the values forν∗ can be precomputed. Also we will see later on that the values fort∗are directly provided by the GroIMP implementation at runtime. This makes calculation of the size of the rate and state vector straightforward.
5.4. Rate Assignment Operator
Figure 5.6.: Class diagram of module hierarchy used in an XL model.
Once memory for rate and state vector has been allocated, for each attribute that Mapping participates in integration an index into these vectors must be assigned. Aseparatedview
on the class hierarchy is created (see figure 5.6c), where attributes have been merged for each class in the order they appear in the derivation chain. The attributes in each separated class can then be assigned offsets starting from zero, and the size of each class can be determined (in the example sizes are σA= 1, σC = 3, σD = 2). For each node, a base index is assigned, so that consecutive elements in rate and state vector belong to this node. For a node nwith base index nbase and a propertyp of this node with offset poffset, the index can be calculated as:
index=nbase+poffset.
What remains to do is to provide the user with a tool to indicate which attributes Rate Assignment Operator
participate in integration and to assign rates to these attributes. Therefore, a new operator :’=was introduced into the language XL, called the deferred rate assignment operator. For a nodenwith a property pa rate rcan be assigned by calling
n [ p ] : ’= r ;
The operator does not modify the attribute itself, but rather writes to a hidden rate vector. Multiple calls to this operator result in the rates being accumulated. There must be no step sizeh anymore, as managing this is up to the integrator.
When working in RGG mode (see [Kni08, appendix B.4]), the user benefits from the fact that the model class then implicitly extends a class RGG, which in turn is derived fromNode. This allows to write rate assignments as
p : ’= r ;
for a (global) propertypof the model and the current instance of RGG is used as node.
At compile time, rate assignments are written to a table as pairs of node type and property. At run time, these data are analyzed to calculate offsets for the attributes and sizes of the node types, as was described above. Then the system performs the steps that were described in section 5.3.2. These are now automatic and don’t need any user assistance, as all necessary data is already available.
For the integrator to execute the rate assignments, the user is responsible to group them into a function getRate. The user must be aware that the proposed state does not necessarily belong to the solution curve, but is solely provided to compute rates.
Enforcing to group all rate assignments into a single function helps the user to stay alert for this.
5.4.1. Circular Transport with Inhibition
The rule for circular transport from sections 5.2.1 and 5.3.1 can now be written as:
x : S −EDGE 0−> y : S −EDGE 0−> z : S : :> { double r a t e = x [ c ] > T ? 0 : µ ∗ y [ c ] ; y [ c ] : ’= −r a t e ;
z [ c ] : ’= +r a t e ; }
It attains the simplicity of the Euler version, while it maintains support for advanced numerical methods. By using a rule, the differential equation can be directly applied to a graph.
5.4.2. Transport in a Tree
Similarly, we can reformulate the tree model in terms of rate assignments (see listing 5.6).
Again, the simplicity of the Euler version from section 5.2.2 is attained (listing 5.1), while keeping the support for better integration methods.
One modification that is necessary is the introduction of a new propertylenfor nodes of typeB. This is due to the fact that the rate assignment operator is only defined on properties of typedouble, which reflects the observation that all investigated libraries for numerical integration work ondoubleas well. However, the inherited attributelength of a cylinder is defined as single precisionfloat.
Technically it would be possible to also define the rate assignment operator for at-tributes of typefloat. When passing values to the integrator, these would have to be converted implicitly to typedouble, and when such attributes are copied into the graph they would be converted back to float. However, this rounding to a lower precision might cause numerical issues and instability in integration, which may result in subtle