GEE/Constraint networks

A constraint network is a way to perform a computation represented by an acyclic graph. When the base library of GEE is installed: it makes available the (gee math constraint) module, which implements a constraint network evaluator built upon the GOOPS extension to Guile.

The module is inspired by/derived from the chapter Propagation of constraints in the book:

 Structure and Interpretation of Computer Programs (sicp)
 Harold Abelson and Gerald Sussman
 Copyright (C) 1996 Massachussets Institute of Technology


The nodes of the acyclic graph are named constraints and the links are named connectors. Constraints hold functions, connectors hold values. The functions accept connector values as arguments and yield connector values as result.

A constraint is satisfied when all its output connectors hold values that can be computed, by the constraint's function, using its input connectors' values. The purpose of a network of constraints is to fill each connector with a value, so that all the constraints are satisfied.


It is possible to "solve" a network by evaluating one of the constraints or by setting/getting the value of a connector. The evaluation propagation rules can be summarised with:


Connectors hold their values until they are reset. If we want to change the input values of the problem and reevaluate the network we have to reset it. The reset propagation rules can be summarised with:

Delayed connectors

Special connectors follow different rules: values can be set and got from them without triggering the evaluation of the network; they propagate the reset request to adjacent object, but do not reset themselves. To reset a special connector we have to explicitly set its value to #f.

One connector, multiple inputs

In this case more than one constraint sets its output value into the same connector. If a connector already holds a value: trying to set it to a different value causes an error to be raised. A predicate can be registered to determine when two values are equal; the default works for real numbers that are approximately equal.

Loading the module

To load the module we do:

 ;; beginning of file 
 (define-module (this) 
   #:use-module (oop goops) 
   #:use-module (gee math constraint) 
   #:duplicates merge-generics) 
 ;; your code goes here 
 ;; end of file 

There are no object maker functions so we need to include the {{{oop goops}}} module.

Simple example

A quick example to get the idea of the interface: let's say that we want to compute the Ohm's law: V = R * I - E = 10 * 5 - 7.5 = 42.5, we code the network like this:

 (let* ((A       (make <constraint> 
                   #:constraint (lambda (R I) 
                                  (list (* R I))))) 
        (B       (make <constraint> 
                   #:constraint (lambda (V_R E) 
                                  (list (- V_R E))))) 
        (result  (begin 
                   (make <connector> 
                     #:name 'R 
                     #:value 10 
                     #:outputs (list A)) 
                   (make <connector> 
                     #:name 'I 
                     #:value 5 
                     #:outputs (list A)) 
                   (make <connector> 
                     #:name 'V_R 
                     #:inputs (list A) 
                     #:outputs (list B)) 
                   (make <connector> 
                     #:name 'E 
                     #:value 7.5 
                     #:outputs (list B)) 
                   (make <connector> 
                     #:name 'V 
                     #:inputs (list B))))) 
   (evaluate! A) 
   (format #t "V = ~A~%" (value result))) 

we notice that we need to put the <constraint> objects into variables, because we need them in inputs and outputs lists, but only the result's connector needs to be put into a variable.


For the full documentation see the Texinfo manual in the distribution. Basically there are three GOOPS classes: <constraint>, <connector> and <delayed-connector>.

To instantiate a <constraint> we do:

 (define A (make <constraint> 
             #:constraint (lambda ...))) 

the constraint functions must accept the input values as arguments and return a list of output values. When making a constraint we do not select its connectors.

To instantiate a <connector> we do:

 (define c (make <connector> 
             #:inputs (list A ...)  
             #:outputs (list B ...))) 

the argument of the #:inputs keyword is the list of input constraints, the arguments of the #:outputs keyword is the list of output constraints.

We make the constraints first, then the connectors. The order of making of the connectors does matter: each connector registers itself in its constraints' states by appending itself to lists. The order in which the connectors appear in these lists is the order of arguments to the constraint's function and the order of values in the function's return value.

Example: to implement:

 connector 1 |              | connector 3
------>------|              |----->-------
             | constraint A |
 connector 2 |              | connector 4
------>------|              |----->-------

we do:

 (let* ((A       (make <constraint> 
                   #:constraint (lambda (one two) 
                                  (list (+ one two) (* one two))))) 
        (one     (make <connector> 
                   #:name 'one 
                   #:value 3 
                   #:outputs (list A))) 
        (two     (make <connector> 
                   #:name 'two 
                   #:value 2 
                   #:outputs (list A))) 
        (sum     (make <connector> 
                   #:name 'sum 
                   #:inputs (list A))) 
        (mul     (make <connector> 
                   #:name 'mul 
                   #:inputs (list A)))) 
  (format #t "sum = ~A, mul = ~A~%" (value sum) (value mul))) 

By making first one and then two: we state that the constraint function of A will be invoked with the value of one first and the value of two next.

By making first sum and then mul: we state that the constraint function of A must have a list of two elements as return value, the first be the value of sum and the second will be the value of mul.

If we make a mess with the network and we are not sure about the order: we can give a name to each connector (with the #:name keyword), and the just write it to the terminal:

 (write A)(newline) 

the string representation shows the lists of names.

Avoiding infinite loops

A constraint network could enter an infinite loop when the values of the inputs depend on the values of the outputs. The second most simple case is:

           connector 1
       |                |
 --------------   --------------
| constraint A | | constraint B |
 --------------   --------------
       |                |
           connector 2

a request to evaluate constraint A causes a request for the value of connector 2, which causes constraint B to be evaluated, which causes a request for the value of constraint 1, which must request an evaluation of constraint A.

The most simple case is:

       |         |

 --------------  |
| constraint A | v connector 1
 --------------  |
       |         |

Such networks can be solved only if we set explicitly the correct value for one of the constraints, that is: if we break the loops. It is mandatory for a constraint network evaluator to detect such cases and output an error; gee math constraint does this at evaluation time: for the first example network:

 (let* ((A       (make <constraint> 
                   #:constraint (lambda (two) (list two)))) 
        (B       (make <constraint> 
                   #:constraint (lambda (one) (list one))))) 
   (make <connector> 
     #:name 'one 
     #:inputs (list A) 
     #:outputs (list B)) 
   (make <connector> 
     #:name 'two 
     #:inputs (list B) 
     #:outputs (list A)) 
   (catch #t 
          (lambda () 
            (evaluate! A)) 
          (lambda (key . args) 
            (apply format #t (cadr args) (caddr args)) 

the output is the message unsolvable connector 'one', and for the second example network:

 (let* ((A       (make <constraint> 
                   #:constraint (lambda (one) (list one))))) 
   (make <connector> 
     #:name 'one 
     #:inputs (list A) 
     #:outputs (list A)) 
   (catch #t 
          (lambda () 
            (evaluate! A)) 
          (lambda (key . args) 
            (apply format #t (cadr args) (caddr args)) 

it is the same unsolvable connector 'one'. But if we set the correct value in the first network (any value in this case, because the constraints are transparent):

 (let* ((A       (make <constraint> 
                   #:constraint (lambda (two) (list two)))) 
        (B       (make <constraint> 
                   #:constraint (lambda (one) (list one)))) 
        (result  (begin 
                   (make <connector> 
                     #:name 'one 
                     #:inputs (list A) 
                     #:outputs (list B)) 
                   (make <connector> 
                     #:name 'two 
                     #:value 10 
                     #:inputs (list B) 
                     #:outputs (list A))))) 
   (evaluate! A) 
   (format #t "result = ~A~%" (value result))) 

the result is 10 and no error is raised.


Many problems require the iteration of the computation's acyclic graph: output values are recycled as input values for the next iterations. (gee math constraint) supports this with delayed connectors.

Let's say that we want to solve the system:

y_1 = -4 * x_1 + 2
y_2 =  5 * x_2 + 3

that is: to find the abscissa x at which y_1 = y_2; to do it we write:

x_2 = (y_2 - 3) / 5

choose an initial guess for x, for example zero, and then iterate with the network:

   x    ---   y   ---  
 -->---| A |-->--| B |--
|       ---       ---   |

where the constraint's function of A is y = -4 * x + 2 and the constraint function of B is x = (y - 3) / 5; to close the loop we use a delayed connector:

 (let* ((A       (make <constraint> 
                   #:constraint (lambda (x) 
                                  (list (+ (* -4 x) 2))))) 
        (B       (make <constraint> 
                   #:constraint (lambda (y) 
                                  (list (/ (- y 3) 5))))) 
        (result  (make <delayed-connector> 
                   #:value 0 
                   #:inputs (list B) 
                   #:outputs (list A))) 
        (tmp     (make <connector> 
                   #:inputs (list A) 
                   #:outputs (list B)))) 
   (define (iterate num) 
     (do ((i 0 (1+ i))) 
         ((> i num)) 
       (evaluate! A) 
       (reset! A))) 
   (define (show) 
     (let ((X    (exact->inexact (value result)))) 
       (format #t "X = ~A, y_1 = ~A, y_2 = ~A~%" X 
               (+ (* -4 X) 2) 
               (+ (* 5 X) 3)))) 
   (iterate 10)(show) 
   (iterate 50)(show) 
   (iterate 100)(show)) 

Guile-GSL ODE solver

The original purpose of existence for (gee math constraint) is to organise the computation of evolution of systems of ordinary differential equations (ODE) using a Guile interface to the GNU Scientific Library (GSL). The interface is part of GEE.

It is always possible to write a The One/total/global/all inclusive system of equations and integrate it using a single ODE solver, but is it often impractical to handle it. Once we have decided to use more than one solver, and to make them evolve in parallel, it is useful to have an infrastructure to share values.

Let's say that we have the following systems (lower case things are vectors, upper case things are matrices):

 d y1
 ---- = A1 y1 + B1 u + C11 y2 + C12 y3
 d t

 d y2
 ---- = A2 y2 + B2 u + C21 y1 + C22 y3
 d t

 d y3
 ---- = A3 y3 + B3 u + C31 y1 + C32 y2
 d t

along with the initial conditions y1(0), y2(0) and y3(0).

We make 3 ODE solvers and we must share the 3 state vectors and the vector of inputs; this is just an example, so the matrices of coefficients are taken from a single 6x6 matrix built to have negative real valued eigenvalues:

 (define A #,(gmr (6 . 6) 
    -0.81407    0.29826    0.37469     0.16042   -0.067375 -0.10942 
     0.29826   -0.74762    0.055107   -0.40254    0.3763   -0.017351 
     0.37469    0.055107  -0.57251     0.0008632 -0.18595  -0.021125 
     0.16042   -0.40254    0.0008632  -1.1071    -0.22968  -0.15622 
    -0.067375   0.3763    -0.18595    -0.22968   -0.7432   -0.074531 
    -0.10942   -0.017351  -0.021125   -0.15622   -0.074531 -0.65386)) 
 ;; 0.0 0.1  0.2 0.3  0.4 0.5 
 ;; 1.0 1.1  1.2 1.3  1.4 1.5 
 ;; 2.0 2.1  2.2 2.3  2.4 2.5 
 ;; 3.0 3.1  3.2 3.3  3.4 3.5 
 ;; 4.0 4.1  4.2 4.3  4.4 4.5 
 ;; 5.0 5.1  5.2 5.3  5.4 5.5 
 (define A1 (elm A '((0 . 0) . (1 . 1)))) 
 (define A2 (elm A '((2 . 2) . (3 . 3)))) 
 (define A3 (elm A '((4 . 4) . (5 . 5)))) 
 (define B1 #,(gmr (2 . 2) 4.1 -4.2 4.3 4.4)) 
 (define B2 #,(gmr (2 . 2) 5.1 -5.2 5.3 5.4)) 
 (define B3 #,(gmr (2 . 2) 6.1 -6.2 6.3 6.4)) 
 (define C11 (elm A '((0 . 2) . (1 . 3)))) 
 (define C12 (elm A '((0 . 4) . (1 . 5)))) 
 (define C21 (elm A '((2 . 0) . (3 . 1)))) 
 (define C22 (elm A '((2 . 4) . (3 . 5)))) 
 (define C31 (elm A '((4 . 0) . (5 . 1)))) 
 (define C32 (elm A '((4 . 2) . (5 . 3)))) 
 (define omega GSL_CONST_NUM_TWICE_PI) 
 (define (u t) 
   (read-gsl-vector-real 2 (cos (* omega t)) (sin (* 5 omega t)))) 
 (define (dy/dt-1 t y1 y2 y3) 
   (+ (dot A1 y1) (dot B1 (u t)) (dot C11 y2)  (dot C12 y3))) 
 (define (dy/dt-2 t y2 y1 y3) 
   (+ (dot A2 y2) (dot B2 (u t)) (dot C21 y1)  (dot C22 y3))) 
 (define (dy/dt-3 t y3 y1 y2) 
   (+ (dot A3 y3) (dot B3 (u t)) (dot C31 y1)  (dot C32 y2))) 
 (define y1-initial #,(gvr 2  -0.5 0.5)) 
 (define y2-initial #,(gvr 2 -1 1)) 
 (define y3-initial #,(gvr 2 -2 2)) 
 (define (make-state-matrix initial points) 
   (let ((m (make-gsl-matrix-real (cons 2 points)))) 
     (set! (column m 0) initial) 
 (define (make-ode solver dy/dt initial time) 
   (make-ode-matrix-evolution (make <ode-system> #:dy/dt dy/dt) solver 
                              #:state-matrix (make-state-matrix initial (dim time)) 
                              #:indep-vector time)) 
 (let* ((solver          (make <ode-solver>)) 
        (time            (linspace 0 4 400)) 
        (ode1            (make-ode solver dy/dt-1 y1-initial time)) 
        (ode2            (make-ode solver dy/dt-2 y2-initial time)) 
        (ode3            (make-ode solver dy/dt-3 y3-initial time)) 
        (cntr1           (make <constraint> 
                           #:constraint (lambda (y2 y3) 
                                          (evolve ode1 y2 y3) 
                                          (list (evolution ode1))))) 
        (cntr2           (make <constraint> 
                           #:constraint (lambda (y1 y3) 
                                          (evolve ode2 y1 y3) 
                                          (list (evolution ode2))))) 
        (cntr3           (make <constraint> 
                           #:constraint (lambda (y1 y2) 
                                          (evolve ode3 y1 y2) 
                                          (list (evolution ode3))))) 
   (make <connector> 
     #:name 'one 
     #:inputs (list cntr1) 
     #:outputs (list cntr2 cntr3)) 
   (make <delayed-connector> 
     #:name 'two 
     #:value y2-initial 
     #:inputs (list cntr2) 
     #:outputs (list cntr1 cntr3)) 
   (make <delayed-connector> 
     #:name 'three 
     #:value y3-initial 
     #:inputs (list cntr3) 
     #:outputs (list cntr1 cntr2)) 
   (evaluate! cntr1) 
   (while (running? ode1) 
     (reset! cntr1) 
     (evaluate! cntr1)) 
   ;; now to access the resulting states: 
   (... (evolution ode1)) 
   (... (evolution ode2)) 
   (... (evolution ode3)) 

(try it with $imulink (muah! ha! ha! ha! ha! ha! haaah!!)).