SICP Exercise 2.16 notes

sicp-216.rkt

#lang sicp

#| 

SICP's exercise 2.16 is marked with a warning that it's very difficult. There
are a handful of one-paragraph or even one-sentence attempts at answering it
floating around the internet, but nothing that satisfied me, so I'm trying to
build on what I found, and what I determined for myself, and write something a
little more complete.

1. Linear expressions only, please!

As defined in the exercise, interval arithemtic is treated sort of like a linear
optimization problem -- each input interval is a constraint, and we check all
the 'corners' of the constrained space -- each combination where a value is
exactly its lower-bound or its upper-bound.

Squaring, in general, offers a clear picture of how wrong this can get. Consider
(x - y)^2, where x = [1 5] and y = [2 3]. The actual minimum of the whole
expression is 0; it occurs when x and y are equal, which doesn't occur at any of
the 'corners' -- it only occurs in the middle of the x interval. Multiplication,
as defined in 2.14, will never return the right answer for this kind of
expression (and also fails on an infinite number of other non-linear
expressions).

2. Can we make a meaningful algebrea of intervals?

We want to be able to add, subtract, multiply and divide intervals in a way that
mimics the way we add, subtract, multiply and divide real numbers, with the
freedom to combine partial results any way we like.

First, we need to be specific about what it means to "mimick the way we ... real
numbers". Luckily, mathematics has already formalized that idea; a set of
operators and operands that behaves like the real numbers is called a "field".
To be a field, our interval algebra must satisfy these constraints:

- Associativity of addition and multiplication: 
    a + (b + c) = (a + b) + c and a · (b · c) = (a · b) · c.

- Commutativity of addition and multiplication: 
    a + b = b + a and a · b = b · a.

- Additive and multiplicative identity: there exist two different elements 0 and 1 in F
    such that a + 0 = a and a · 1 = a.

- Additive inverses: for every a in F, there exists an element in F, denoted -a, 
    called the additive inverse of a, such that a + (-a) = 0.

- Multiplicative inverses: for every a ≠ 0 in F, there exists an element in F, 
    denoted by a^-1, 1/a, called the multiplicative inverse of a, such that a · a^-1 = 1.

- Distributivity of multiplication over addition: 
    a · (b + c) = (a · b) + (a · c).

The solution SICP proposes has some of these properties, but fails to have a 
multiplicative inverse: a · b / b != a for every a and b != 0.

I *believe* no appropriate field for intervals exists:

Each interval involces two numbers, a lower-bound and an upper-bound.

Adding intervals treats upper-bounds and lower-bounds separately; there's an
operation for lower bounds, and an operation for upper bounds, each of which
satisfies the rules for addition in a field (because it's literally addition!)

Multiplying intervals, however, mixes the streams, especially when they cross
zero, and loses information (sometimes an incoming bound gets completely
ignored); SICP even notes that division is undefined for all intervals that
cross zero. We're allowed to ignore the additive identity ([0 0]) but interval
division must ignore an infinite set, not just the identity.

3. Can we define a method for calculating arbitrary expressions, even if the
results of those expressions can't be combined later?

A lot of solutions people offer to 2.16 involve doing some manipulation of an
entire expression before doing the calculation. Some people want to try and
reduce the expression so that it only includes each incoming interval once --
this idea is hinted at by the problem. Unfortunately, that can only be done for
linear combinations of intervals; even expressions as simple as x^2, or
(x+y)/(x+z) inevitably involve repetition. See 1, above :). 

(One could go further, and symbolically differentiate the expression, find
zeros, and include them in the 'corners' to be explored, but still it would fail
for non-differentiatable expressions)

4. Can we produce *approximate* solutions for arbitrary expressions, even if the
results of those expressions can't be combined later?

A more compelling option is monte-carlo simulation. MC simulation would involve
repeatedly picking a random value within the range of each interval, and
calculating the expression using those values. Do this thousands of times, and
get thousands of possible outputs. The range of these potential outputs is an
approximation, but always narrower-or-equal to the true total range of the
expression. 

Monte carlo is NOT efficient, but it DOES have advantages -- not the least of
which is that it allows the use of more realistic distributions of values within
an interval, for those cases when an error range is normally distributed (a bell
curve rather than a flat line).

Simulation is not a satisfyingly pure algebraic solution, but it *is*
functional, and would be reasonably easy to write if we'd learned about
side-effects yet, or had macros to play with. Technically, we don't, so here's
an interval math system that uses monte-carlo simulation, built using techniques
we've (barely) learned so far, plus three tiny new ideas:

- (random 1.0) returns a random number between 0 and 1.
- (par? x) returns true if x is a cons pair, and false otherwise
- '() is nil, a default, falsy, non-pair-y value

The general approach taken below is to represent each interval as a looong list
of (randomly selected) potential values within that interval. Doing arithmetic
on two intervals means doing that arithemetic on the first value of each list,
then the second value of each list, then the third... building up an equally
long list of potential results. This means we're treating every operation the
way we had been treating addition (we never 'cross the streams') so our values
and operations form a field.

The lower-bound of an interval is the lowest value in the list; and the
upper-bound is the highest.

|#

(define simulations 5000) ; higher means slower but more accurate

(define (nth n interval)
  ; returns the nth value in a cons list
  (if (= n 0)
      (car interval)
      (nth (- n 1) (cdr interval))))

(define (random-value-in x y)
  ; returns a random value between x and y
  (+ (* (random 1.0) (- y x)) x))

(define (n-random-values-in n x y)
  ; returns a const list of n random values between x and y
  (if (> n 0)
      (cons (random-value-in x y) (n-random-values-in (- n 1) x y))
      '()))

(define (make-interval x y)
  ; An interval is represented as a cons list of potential values
  (n-random-values-in simulations x y))

(define (make-interval-arithmetic fn)
  #|
  An interval-math function performs the non-interval function
  on each of the simulated values of its arguments in turn.
  This function turns the non-interval function fn into its interval
  counterpart.
  |#
  (define (res x y)
    (define (simulation n)
      (if (< n simulations)
          (cons
           (fn (nth n x) (nth n y))
           (simulation (+ n 1)))
          '()))
    (simulation 0))
  res)

(define add-interval (make-interval-arithmetic +))
(define sub-interval (make-interval-arithmetic -))
(define mul-interval (make-interval-arithmetic *))
(define div-interval (make-interval-arithmetic /))

(define (upper-bound interval)
  ; The upper-bound of an interval is the largest of its possible values
  (define (find-max max interval)
    (if (pair? interval)
        (if (> max (car interval))
            (find-max max (cdr interval))
            (find-max (car interval) (cdr interval)))
        max))
  (find-max (car interval) interval))

(define (lower-bound interval)
  ; The lower-bound of an interval is the smallest of its possible values
  (define (find-min min interval)
    (if (pair? interval)
        (if (< min (car interval))
            (find-min min (cdr interval))
            (find-min (car interval) (cdr interval)))
        min))
  (find-min (car interval) interval))

(define (print-interval i)
  (newline)
  (display "[")
  (display (lower-bound i))
  (display " - ")
  (display (upper-bound i))
  (display "]"))

#| USAGE: Same as in 2.1.4, just slower. |#


; The example from the text which fails for their approach
(define (par1 r1 r2)
  (div-interval (mul-interval r1 r2)
                (add-interval r1 r2)))
(define (par2 r1 r2)
  (let ((one (make-interval 1 1))) 
    (div-interval one
                  (add-interval (div-interval one r1)
                                (div-interval one r2)))))

(define r1 (make-interval 9 11))
(define r2 (make-interval 19 21))
(print-interval (par1 r1 r2))
(print-interval (par2 r1 r2))

; Both methods produce the same result! But that result is a narrower
; interval than the algebraically-correct result.

; My example from above of failing for non-linear expressions:

(define (square-diff x y)
  (let ((diff (sub-interval x y)))
    (mul-interval diff diff)))

(define i1 (make-interval 1 5))
(define i2 (make-interval 2 3))
(print-interval (square-diff i1 i2))

; Note how poor the accuracy is for some conditions. The lower-bound is very
; close to 0, because there are many combinations of values that give that result;
; but the upper bound is several tenths away from 9, because only one specific
; combination of potential values will yeild that result.