[Racket] Ferti hydroponic nutrient solver, redux.
First jab at the Lawson-Hanson NNLS algorithm.
services/nnls.rkt
@@ -0,0 +1,237 @@
1
Added:
#lang racket
2
Added:
3
Added:
(provide find-ferti-recipe)
4
Added:
5
Added:
(require math/array
6
Added:
math/matrix
7
Added:
"../models/nutrient.rkt"
8
Added:
"../models/nutrient-measurement.rkt"
9
Added:
"../models/nutrient-target.rkt"
10
Added:
"../models/fertilizer-product.rkt")
11
Added:
12
Added:
(define (find-ferti-recipe)
13
Added:
(define fertilizers (get-fertilizer-products))
14
Added:
(define solution-array (find-nnls fertilizers))
15
Added:
(for/list ([i (length (get-fertilizer-products))])
16
Added:
(cons (list-ref fertilizers i)
17
Added:
(array-ref solution-array (vector i 0)))))
18
Added:
19
Added:
(define (find-nnls fertilizers)
20
Added:
(define nutrients (get-nutrients))
21
Added:
(define fertilizer-product-matrix (get-fertilizer-product-matrix nutrients fertilizers))
22
Added:
(define deficits
23
Added:
(->col-matrix
24
Added:
(for/list ([n nutrients])
25
Added:
(define latest-measurement (get-latest-nutrient-measurement-value n))
26
Added:
(define latest-target (get-latest-nutrient-target-value n))
27
Added:
(define deficit
28
Added:
(cond
29
Added:
[(or (false? latest-measurement)
30
Added:
(zero? latest-measurement))
31
Added:
latest-target]
32
Added:
[(false? latest-target)
33
Added:
0]
34
Added:
[(and (number? latest-measurement)
35
Added:
(number? latest-target))
36
Added:
(* 100
37
Added:
(/ (- latest-target latest-measurement)
38
Added:
latest-measurement))]
39
Added:
[else (error "either the target or measurement are not numbers")]))
40
Added:
deficit)))
41
Added:
(define error-threshold 10e-4)
42
Added:
(lawson-hanson-1974 fertilizer-product-matrix deficits error-threshold))
43
Added:
44
Added:
;; Algorithm lifted from the Wikipedia article on NNLS
45
Added:
(define/contract (lawson-hanson-1974 A y ε)
46
Added:
;;;;;;;;;
47
Added:
;; Inputs
48
Added:
;;;;;;;;;
49
Added:
50
Added:
(-> matrix? ; Real-valued matrix A of dimension m × n
51
Added:
col-matrix? ; Real-valued column matrix (vector) y of dimension m
52
Added:
real? ; Real-value ɛ, tolerance for the stopping criterion
53
Added:
col-matrix?) ; Real-valued solution column matrix x
54
Added:
(define-values
55
Added:
(m ; Number of nutrients
56
Added:
n) ; Number of fertilizer products
57
Added:
(matrix-shape A))
58
Added:
59
Added:
60
Added:
;;;;;;;;;;;;;
61
Added:
;; Initialize
62
Added:
;;;;;;;;;;;;;
63
Added:
64
Added:
;; The passive set P is initially empty.
65
Added:
(define P (mutable-set))
66
Added:
;; The active set R initially contains the indexes to the nutrients allowed to be...
67
Added:
(define R (list->mutable-set (range n)))
68
Added:
69
Added:
(define (colv-ref v i)
70
Added:
(matrix-ref v i 0))
71
Added:
72
Added:
;; Gradient-like vector for residual error.
73
Added:
(define (compute-w x)
74
Added:
(matrix* (matrix-transpose A)
75
Added:
(matrix- y (matrix* A x))))
76
Added:
77
Added:
;; max over j in R of w_j, returning (values max-val j*)
78
Added:
(define (max-w-in-R w)
79
Added:
(for/fold ([max-val -inf.0] [max-j #f])
80
Added:
([j (in-set R)])
81
Added:
(define v (colv-ref w j))
82
Added:
(if (> v max-val)
83
Added:
(values v j)
84
Added:
(values max-val max-j))))
85
Added:
86
Added:
;; Build full candidate vector s from current P:
87
Added:
;; s_P = (A_Pᵀ A_P)⁻¹ A_Pᵀ y, s_R = 0
88
Added:
(define (make-s-from-P)
89
Added:
(if (set-empty? P)
90
Added:
(make-matrix n 1 0)
91
Added:
(let* ([idxs (sort (set->list P) <)]
92
Added:
[AP (submatrix A (::) idxs)]
93
Added:
[sP (matrix* (matrix-inverse
94
Added:
(matrix* (matrix-transpose AP) AP))
95
Added:
(matrix-transpose AP)
96
Added:
y)])
97
Added:
;; map: column index j in P -> corresponding sP entry
98
Added:
(define mapping
99
Added:
(for/list ([j idxs] [k (in-naturals)])
100
Added:
(cons j (colv-ref sP k))))
101
Added:
(define (s-at i)
102
Added:
(define p (assoc i mapping))
103
Added:
(if p (cdr p) 0))
104
Added:
(build-matrix n 1 (λ (i j) (s-at i))))))
105
Added:
106
Added:
;; The "first try" x represents no addition of any fertilizer.
107
Added:
(define x (make-matrix n 1 0))
108
Added:
109
Added:
110
Added:
;;;;;;;;;;;;;
111
Added:
;; Outer loop
112
Added:
;;;;;;;;;;;;;
113
Added:
114
Added:
(let outer-loop ()
115
Added:
(define w (compute-w x))
116
Added:
117
Added:
(cond
118
Added:
;; If no remaining candidates in R, we're done.
119
Added:
[(set-empty? R)
120
Added:
x]
121
Added:
122
Added:
[else
123
Added:
(define-values (max-val j*) (max-w-in-R w))
124
Added:
125
Added:
;; Stopping criterion: max(w_R) <= ε
126
Added:
(cond
127
Added:
[(or (not j*) (<= max-val ε))
128
Added:
x]
129
Added:
130
Added:
[else
131
Added:
;; Add j* to P, remove from R
132
Added:
(set-remove! R j*)
133
Added:
(set-add! P j*)
134
Added:
135
Added:
;; Inner loop: adjust until s_P > 0
136
Added:
(let inner-loop ()
137
Added:
(define s (make-s-from-P))
138
Added:
139
Added:
;; min(s_P)
140
Added:
(define min-sP
141
Added:
(if (set-empty? P)
142
Added:
+inf.0
143
Added:
(for/fold ([mn +inf.0])
144
Added:
([j (in-set P)])
145
Added:
(min mn (colv-ref s j)))))
146
Added:
147
Added:
(cond
148
Added:
;; If all s_P > 0 (or P empty), accept s as new x and go back to outer loop
149
Added:
[(or (set-empty? P) (> min-sP 0))
150
Added:
(set! x s)
151
Added:
(outer-loop)]
152
Added:
153
Added:
[else
154
Added:
;; Compute α = min_{i in P, s_i <= 0} x_i / (x_i - s_i)
155
Added:
(define α
156
Added:
(for/fold ([a +inf.0])
157
Added:
([j (in-set P)])
158
Added:
(define sj (colv-ref s j))
159
Added:
(if (<= sj 0)
160
Added:
(let* ([xj (colv-ref x j)]
161
Added:
[den (- xj sj)])
162
Added:
(if (> den 0)
163
Added:
(min a (/ xj den))
164
Added:
a))
165
Added:
a)))
166
Added:
167
Added:
(when (or (equal? α +inf.0) (<= α 0))
168
Added:
(error 'lawson-hanson-1974 "no valid α in inner loop"))
169
Added:
170
Added:
;; x ← x + α (s − x)
171
Added:
(define new-x
172
Added:
(matrix+ x (matrix* α (matrix- s x))))
173
Added:
174
Added:
;; Move to R all indices j in P with x_j <= 0
175
Added:
(define to-remove '())
176
Added:
(for ([j (in-set P)])
177
Added:
(when (<= (colv-ref new-x j) 0)
178
Added:
(set! to-remove (cons j to-remove))))
179
Added:
(for ([j to-remove])
180
Added:
(set-remove! P j)
181
Added:
(set-add! R j))
182
Added:
183
Added:
(set! x new-x)
184
Added:
(inner-loop)]))])])))
185
Added:
186
Added:
(define (get-fertilizer-product-matrix nutrients fertilizers)
187
Added:
;; Lines are nutrients, columns are fertilizers
188
Added:
(build-matrix (length nutrients)
189
Added:
(length fertilizers)
190
Added:
(λ (i j)
191
Added:
(define selected-nutrient (list-ref nutrients i))
192
Added:
(define product (list-ref fertilizers j))
193
Added:
(define pair (assoc selected-nutrient
194
Added:
(fertilizer-product-values product)))
195
Added:
(if pair (cdr pair) 0))))
196
Added:
197
Added:
(module+ test
198
Added:
(require rackunit
199
Added:
rackunit/text-ui
200
Added:
"../db/conn.rkt"
201
Added:
"../db/migrations.rkt")
202
Added:
203
Added:
(define test-date "2025-01-01")
204
Added:
205
Added:
(run-tests
206
Added:
(test-suite
207
Added:
"Nutrient measurement model"
208
Added:
#:before (λ ()
209
Added:
(connect! #:path 'memory)
210
Added:
;; (connect! #:path "test.sqlite3")
211
Added:
(migrate-all!)
212
Added:
213
Added:
(define nitrogen (create-nutrient! "Nitrogen" "N"))
214
Added:
(define phosphorus (create-nutrient! "Phosphorus" "P"))
215
Added:
216
Added:
(create-nutrient-measurement! test-date
217
Added:
`((,nitrogen . 0)
218
Added:
(,phosphorus . 0)))
219
Added:
(create-nutrient-target! test-date
220
Added:
`((,nitrogen . 100)
221
Added:
(,phosphorus . 50)))
222
Added:
223
Added:
(create-fertilizer-product! "King Nitrogen"
224
Added:
`((,nitrogen . 100)))
225
Added:
(create-fertilizer-product! "Phosphorescent Baboon"
226
Added:
`((,nitrogen . 10)
227
Added:
(,phosphorus . 100)))
228
Added:
(create-fertilizer-product! "John's Phosphorus"
229
Added:
`((,nitrogen . 3)
230
Added:
(,phosphorus . 30))))
231
Added:
#:after (λ ()
232
Added:
(disconnect!))
233
Added:
234
Added:
(test-case "Solve for NNLS"
235
Added:
(displayln
236
Added:
(format "Final solution for fertilizers is combination ~a"
237
Added:
(find-ferti-recipe)))))))