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:

- the evaluation of a constraint gets the values from all its input connectors and sets the values of all its output connectors;
- getting the value from an empty connector causes all its input constraints to be evaluated;

- setting the value of a connector causes all its output constraints to be evaluated.

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:

- resetting a constraint resets all its output connectors and all its input connectors;

- resetting a connector resets all its output constraints and all its input constraints.

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`.

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.

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.

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.

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))
(newline))))
```

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))
(newline))))
```

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))
```

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) m)) (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!!)).