nested for loops and suggested alternatives

[travis.h...@gmail.com]

Hi All,

I'm an R programmer that has recently started learning Racket. I decided to start by trying to create a simple age-structured population model. In R, I would initialize a matrix and use nested for loops to move through the elements of the matrix and propagate the population forward through time. For my first attempt in Racket (https://gist.github.com/hinkelman/3ee6115cdd7f0a4c8f1672b7d8df5c27), I used for* to loop through a vector of vectors. The code in that gist doesn't quite work. There is apparently something wrong with how I'm using vector-set! such that "rows" in my vector of vectors are being updated prematurely. I would greatly appreciate it if someone could point out what I'm doing wrong there. I'm also interested in suggestions for alternative approaches because it seems unlikely that I have written this code in idiomatic Racket.

Thanks,

Travis

P.S. If it is helpful, here is a gist (https://gist.github.com/hinkelman/d5b8414b0c6383057d7846509a724bbf) with the R code that I was trying to write in Racket.

[Alex Harsanyi]

This line looks suspicious:

     (define results (make-vector years (make-vector (vector-length fecundity) 0)))

The "(make-vector (vector-length fecundity) 0)" expression will create a single vector, than it creates the outer vector will all elements pointing to it.  It is not a matrix, but a "column" vector where each element is referencing the same row vector.  This means that if you update an element in one of the rows, the same value will "appear" in all other rows. The only row that is different is the first one which you initialize in the line:

    (vector-set! results 0 (make-vector (vector-length fecundity) 10))

What you probably want is a vector of vectors, which can be built like this

    (define results (for/vector ([index (in-range years)]) (make-vector (vector-length fecundity) 0)))

Alex.

[Philip McGrath]

You may also want to look into the math/array and math/matrix libraries.

-Philip

--
You received this message because you are subscribed to the Google Groups "Racket Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to racket-users...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Alex Harsanyi]

On Sunday, February 10, 2019 at 4:47:52 PM UTC+8, Philip McGrath wrote:

You may also want to look into the math/array and math/matrix libraries.

looking at the code, it seems to me that each iteration only the "current" and "next" generation are needed.  The best structure might be a list of vectors which is constructed step by step -- every iteration a new vector is added to this list.   This assumes that the entire "evolution" is needed at the end.  If all is needed at the end is the "current" generation, it might be more efficient to just use two vectors, the current and next one and swap them each iteration.

But before we get into improving the code, let's get it working first :-)

Alex.

[travis.h...@gmail.com]

Yes, this was the problem. I now have results that match the output from R. I've updated the gist with your line for the correct way to create a vector of vectors. I will have to spend some more time to understand the make-vector behavior. Perhaps my thinking is too constrained by my R experience where many functions are "vectorized"? I guess the part that was confusing me (and still is) is why make-vector worked as expected (by me) for my inner vector but not my outer vector. Thanks!

[Daniel Prager]

Hi Travis

Would you mind sharing your results from R (and now Racket)?

I'm having a play with the Racket code, and would like to check I'm not breaking anything as I refactor it.

Dan

[travis.h...@gmail.com]

I played around a bit with the math/matrix library. Actually, my first idea for a Racket learning project was to rewrite the code in the R popbio package (https://cran.r-project.org/web/packages/popbio/popbio.pdf) in Racket. But I bailed on that when I saw that eigendecomposition was still on the list of unimplemented but useful algorithms for math/matrix.

Performance is not important for this example but this math/array warning gave me the impression that perhaps it was not the best place to start for me. 

Performance Warning: Indexing the elements of arrays created in untyped Racket is currently 25-50 times slower than doing the same in Typed Racket, due to the overhead of checking higher-order contracts. We are working on it.

I also thought it would be good to stick to racket/base for my early forays with Racket but I would not give that same advice to someone learning R (i.e., some R packages greatly ease the introduction to R).

[travis.h...@gmail.com]

I've attached a screenshot of the resulting matrix from R. Hopefully, that is sufficient.

[Jens Axel Søgaard]

Performance Warning: Indexing the elements of arrays created in untyped Racket is currently 25-50 times slower than doing the same in Typed Racket, due to the overhead of checking higher-order contracts. We are working on it.

I have no idea how to improve that (the timings are old - perhaps it has become better already). 

/Jens Axel

--

--
Jens Axel Søgaard

[Daniel Prager]

Thanks for the screenshot Travis.

Just for fun, here's a version in Racket that eschews assignment, vectors and for loops, in favour of recursion and lists ...

#lang racket

(define years 30)

(define prop-female 0.5)

(define egg-surv 0.6)

(define fecundity '(0 0 200 400 800))

(define survival '(0.2 0.4 0.6 0.8 0))

(define capacity '(1e6 1e5 1e4 1e3 1e2 -9999))

(define cap0 (first capacity))

(define (beverton-holt N p c) (/ N (+ (/ 1 p) (/ N c))))

(define (evolve N [f fecundity] [s survival] [cap (rest capacity)] [Nt0 0] [Nt null])

    (if (null? f)

        (cons Nt0 (reverse Nt))

        (evolve (rest N) (rest f) (rest s) (rest cap)

                (+ Nt0 (if (= (first f) 0)

                           0

                           (beverton-holt (* (first N) prop-female)

                                          (* (first f) egg-surv)

                                          (- cap0 Nt0))))

                (if (= (first s) 0)

                    Nt

                    (cons (beverton-holt (first N) (first s) (first cap)) Nt)))))

(define (iterate N n [i 1])

  (displayln (list i N))

  (unless (= i n) (iterate (evolve N) n (+ i 1))))

(iterate (make-list (length fecundity) 10) years)

[Gustavo Massaccesi]

To understand the problem, it will be useful to understand the difference between

-- (make-vector 5 (make-vector 5))
-- (build-vector 5 (lambda (_) (build-vector 5 (lambda (_) 0))))
-- (for/vector ([_ 5]) (for/vector ([_ 5]) 0))

Gustavo

--

[Laurent]

perhaps more importantly:

-- (build-vector 5 (lambda (_) (make-vector 5 0)))

-- (make-vector 5 (build-vector 5 (lambda (_) 0)))

[travis.h...@gmail.com]

Thanks, Daniel, this is helpful. I think that I understand your code, but it is a still a foreign way of thinking for me. Of course, that is a big part of why I'm learning Racket, i.e., to make programming with lists and recursion more natural. One key gap for me is how to build up data structures using lists and recursion. Perhaps the answer is that you don't build up the data structure in memory but by writing to disk. In this particular example, it is useful to have the values at every iteration for plotting the population trajectory (by age class or summed across age classes). 

Below is a little example I was exploring that involves matrix multiplication (col-matrix is abundance in different age classes; square matrix is transition probabilities among age classes). It was obvious to me how to use recursion to get the final col-matrix of abundances but not how to build up a data structure that included the abundances at every iteration.

(require math/matrix)

(define A (matrix [[0 0 5.905]

                   [0.368 0.639 0.025]

                   [0.001 0.152 0.051]]))

(define n (col-matrix [5 5 5]))

(define (pop-projection A n iter)

  (if (zero? iter) n

      (pop-projection A (matrix* A n) (- iter 1))))

(pop-projection A n 25)

[Alex Harsanyi]

On Monday, February 11, 2019 at 8:01:47 AM UTC+8, travis.h...@gmail.com wrote:

Thanks, Daniel, this is helpful. I think that I understand your code, but it is a still a foreign way of thinking for me. Of course, that is a big part of why I'm learning Racket, i.e., to make programming with lists and recursion more natural. One key gap for me is how to build up data structures using lists and recursion. Perhaps the answer is that you don't build up the data structure in memory but by writing to disk. In this particular example, it is useful to have the values at every iteration for plotting the population trajectory (by age class or summed across age classes). 

Below is a little example I was exploring that involves matrix multiplication (col-matrix is abundance in different age classes; square matrix is transition probabilities among age classes). It was obvious to me how to use recursion to get the final col-matrix of abundances but not how to build up a data structure that included the abundances at every iteration.

One way to do this is for `pop-abundances` to have an extra parameter, the list of previous abundances, and whenever the function is called recursively, it adds the current abundance to this list and passes it on to the next call.  The final call will than return this result instead of the last abundance.  In the example below, "cons" adds to the front of the list, so "result" contains the most recent values first and  this list is reversed before being returned to the user.  Also, when `pop-abundances` is invoked by the user, there are no "previous abundances" , so it needs to be invoked with an empty list -- this is handled by a default parameter for 'result':

#lang racket

(require math/matrix)

(define A (matrix [[0 0 5.905]

                   [0.368 0.639 0.025]

                   [0.001 0.152 0.051]]))

(define n (col-matrix [5 5 5]))

(define (pop-projection A n iter [result '()])

  (if (zero? iter)

      (reverse (cons n result))

      (pop-projection A (matrix* A n) (- iter 1) (cons n result))))

(pop-projection A n 25)

Alex. 

[Matthias Felleisen]

On Feb 10, 2019, at 7:26 PM, Alex Harsanyi <alexha...@gmail.com> wrote:

One way to do this is for `pop-abundances` to have an extra parameter, the list of previous abundances, and whenever the function is called recursively, it adds the current abundance to this list and passes it on to the next call.  The final call will than return this result instead of the last abundance.  In the example below, "cons" adds to the front of the list, so "result" contains the most recent values first and  this list is reversed before being returned to the user.  Also, when `pop-abundances` is invoked by the user, there are no "previous abundances" , so it needs to be invoked with an empty list -- this is handled by a default parameter for 'result':

#lang racket

(require math/matrix)

(define A (matrix [[0 0 5.905]

                   [0.368 0.639 0.025]

                   [0.001 0.152 0.051]]))

(define n (col-matrix [5 5 5]))

(define (pop-projection A n iter [result '()])

  (if (zero? iter)

      (reverse (cons n result))

      (pop-projection A (matrix* A n) (- iter 1) (cons n result))))

(pop-projection A n 25)

Don’t use accumulators if the function has all the information: 

#lang racket

(require math/matrix)

(define A

  (matrix

   [[0     0     5.905]

    [0.368 0.639 0.025]

    [0.001 0.152 0.051]]))

(define n (col-matrix [5 5 5]))

(define (pop-projection A n iter [result '()])

  (if (zero? iter)

      (reverse (cons n result))

      (pop-projection A (matrix* A n) (- iter 1) (cons n result))))

(define (pop-projection.v1 A n iter)

  (if (zero? iter)

      (list n)

      (cons n (pop-projection A (matrix* A n) (- iter 1)))))

(equal? (pop-projection A n 25) (pop-projection.v1 A n 25))

[travis.h...@gmail.com]

For the benefit of other beginners, I think there was a small typo in prop-projection.v1 where the intention was to call that function recursively rather than calling pop-projection in the body of pop-projection.v1. The typo is corrected below:

(define (pop-projection.v1 A n iter)

  (if (zero? iter)

      (list n)

      (cons n (pop-projection.v1 A (matrix* A n) (- iter 1)))))

Thanks to everyone for the feedback. It has been very illuminating!

[Daniel Prager]

Hi Travis

Glad you found it instructive.

In mathematical terms these kinds of systems map from the current state at t=n to the next time-step at t=n+1. 

Given a transition function I tend to just use for/list (or for/vector) to build up a history.

E.g.

#lang racket

(require math/matrix)

(define A (matrix [[0 0 5.905]

                   [0.368 0.639 0.025]

                   [0.001 0.152 0.051]]))

(define n (col-matrix [5 5 5]))

(define (iterate A n iter)

  (for/list ([i iter])

    (define temp n)

    (set! n (matrix* A n))

    (list i temp)))

(iterate A n 25)

You could use for/fold, which avoids the use of set!, but the syntax is a little heavier.

Dan