1 #lang racket 2 3 (provide find-ferti-recipe) 4 5 (require math/array 6 math/matrix 7 "../models/nutrient.rkt" 8 "../models/fertilizer-product.rkt" 9 "../models/nutrient-measurement.rkt" 10 "../models/crop-rotation.rkt" 11 "../models/crop-requirement.rkt") 12 13 (define (find-ferti-recipe date-string) 14 (define fertilizers (get-fertilizer-products)) 15 (define measurement-values 16 (nutrient-measurement-values 17 (or (get-nutrient-measurement #:date date-string) 18 (error 'nnls 19 "Can't compute the Ferti recipe (missing nutrient measurement for ~a)" 20 date-string)))) 21 (define rotation-values 22 (average-crop-requirement-nutrient-values 23 (crop-rotation-requirements 24 (or 25 (get-crop-rotation #:date date-string) 26 (error 'nnls "Can't compute the Ferti recipe (missing crop rotation for ~a)" date-string))))) 27 (define solution-array (solve-nnls fertilizers measurement-values rotation-values)) 28 (for/hash ([fertilizer (in-list fertilizers)] 29 [quantity (in-array solution-array)]) 30 (values fertilizer quantity))) 31 32 (define (solve-nnls fertilizers measurement-values rotation-values) 33 (define nutrients (get-nutrients)) 34 (define fertilizer-product-matrix (get-fertilizer-product-matrix nutrients fertilizers)) 35 (define deficits 36 (->col-matrix (for/list ([n nutrients]) 37 (define measured (hash-ref measurement-values n 0)) 38 (define required (hash-ref rotation-values n 0)) 39 (define deficit 40 (if (zero? required) 41 0 42 (* 100 (/ (- required measured) required)))) 43 deficit))) 44 (define error-threshold 10e-4) 45 (lawson-hanson-1974 fertilizer-product-matrix deficits error-threshold)) 46 47 ;; Algorithm lifted from the Wikipedia article on NNLS 48 (define (lawson-hanson-1974 A y ε) 49 ;;;;;;;;; 50 ;; Inputs 51 ;;;;;;;;; 52 53 (-> matrix? ; Real-valued matrix A of dimension m × n 54 col-matrix? ; Real-valued column matrix (vector) y of dimension m 55 real? ; Real-value ɛ, tolerance for the stopping criterion 56 col-matrix?) ; Real-valued solution column matrix x 57 (define-values (m ; Number of nutrients 58 n) ; Number of fertilizer products 59 (matrix-shape A)) 60 61 ;;;;;;;;;;;;; 62 ;; Initialize 63 ;;;;;;;;;;;;; 64 65 ;; The passive set P is initially empty. 66 (define P (mutable-set)) 67 ;; The active set R initially contains the indexes to the nutrients allowed to be... 68 (define R (list->mutable-set (range n))) 69 70 (define (colv-ref v i) 71 (matrix-ref v i 0)) 72 73 ;; Gradient-like vector for residual error. 74 (define (compute-w x) 75 (matrix* (matrix-transpose A) (matrix- y (matrix* A x)))) 76 77 ;; max over j in R of w_j. 78 (define (max-w-in-R w) 79 (if (set-empty? R) 80 (values -inf.0 #f) 81 (let ([max-j (argmax (λ (j) (colv-ref w j)) (set->list R))]) 82 (values (colv-ref w max-j) max-j)))) 83 84 ;; Build full candidate vector s from current P: 85 ;; s_P = (A_Pᵀ A_P)⁻¹ A_Pᵀ y, s_R = 0 86 (define (make-s-from-P) 87 (if (set-empty? P) 88 (make-matrix n 1 0) 89 (let* ([idxs (sort (set->list P) <)] 90 [AP (submatrix A (::) idxs)] 91 [sP (matrix* (matrix-inverse (matrix* (matrix-transpose AP) AP)) 92 (matrix-transpose AP) 93 y)]) 94 ;; map: column index j in P -> corresponding sP entry 95 (define mapping 96 (for/list ([j idxs] 97 [k (in-naturals)]) 98 (cons j (colv-ref sP k)))) 99 (define (s-at i) 100 (define p (assoc i mapping)) 101 (if p 102 (cdr p) 103 0)) 104 (build-matrix n 1 (λ (i j) (s-at i)))))) 105 106 ;; The "first try" x represents no addition of any fertilizer. 107 (define x (make-matrix n 1 0)) 108 109 ;;;;;;;;;;;;; 110 ;; Outer loop 111 ;;;;;;;;;;;;; 112 113 (let outer-loop () 114 (define w (compute-w x)) 115 116 (cond 117 ;; If no remaining candidates in R, we're done. 118 [(set-empty? R) x] 119 120 [else 121 (define-values (max-val j*) (max-w-in-R w)) 122 123 ;; Stopping criterion: max(w_R) <= ε 124 (cond 125 [(or (not j*) (<= max-val ε)) x] 126 127 [else 128 ;; Add j* to P, remove from R 129 (set-remove! R j*) 130 (set-add! P j*) 131 132 ;; Inner loop: adjust until s_P > 0 133 (let inner-loop () 134 (define s (make-s-from-P)) 135 (define min-sP 136 (if (set-empty? P) 137 +inf.0 138 (for/fold ([mn +inf.0]) ([j (in-set P)]) 139 (min mn (colv-ref s j))))) 140 (cond 141 ;; If all s_P > 0 (or P empty), accept s as new x and go back to outer loop 142 [(or (set-empty? P) (> min-sP 0)) 143 (set! x s) 144 (outer-loop)] 145 146 [else 147 ;; Compute α = min_{i in P, s_i <= 0} x_i / (x_i - s_i) 148 (define α 149 (for/fold ([a +inf.0]) ([j (in-set P)]) 150 (define sj (colv-ref s j)) 151 (if (<= sj 0) 152 (let* ([xj (colv-ref x j)] 153 [den (- xj sj)]) 154 (if (> den 0) 155 (min a (/ xj den)) 156 a)) 157 a))) 158 159 (when (or (equal? α +inf.0) (<= α 0)) 160 (error 'lawson-hanson-1974 "no valid α in inner loop")) 161 162 ;; x ← x + α (s − x) 163 (define new-x (matrix+ x (matrix-scale (matrix- s x) α))) 164 165 ;; Move to R all indices j in P with x_j <= 0 166 (define to-remove '()) 167 (for ([j (in-set P)]) 168 (when (<= (colv-ref new-x j) 0) 169 (set! to-remove (cons j to-remove)))) 170 (for ([j to-remove]) 171 (set-remove! P j) 172 (set-add! R j)) 173 174 (set! x new-x) 175 (inner-loop)]))])]))) 176 177 (define (get-fertilizer-product-matrix nutrients fertilizers) 178 ;; Lines are nutrients, columns are fertilizers 179 (build-matrix (length nutrients) 180 (length fertilizers) 181 (λ (i j) 182 (define selected-nutrient (list-ref nutrients i)) 183 (define product (list-ref fertilizers j)) 184 (hash-ref (fertilizer-product-values product) selected-nutrient 0)))) 185 186 (module+ test 187 (require rackunit 188 rackunit/text-ui 189 "../db/conn.rkt" 190 "../db/migrations.rkt") 191 192 (run-tests 193 (test-suite "NNLS" 194 #:before (λ () 195 (connect! #:path 'memory) 196 (migrate-all!) 197 (create-nutrient! "Nitrogen" "Azote" "N") 198 (create-nutrient! "Phosphorus" "Phosphore" "P") 199 (create-nutrient! "Potassium" "Potassium" "K")) 200 #:after (λ () (disconnect!)) 201 202 (test-case "Build fertilizer product matrix" 203 (define n (get-nutrient #:name "Nitrogen")) 204 (define p (get-nutrient #:name "Phosphorus")) 205 206 (define nutrients (list n p)) 207 208 (define f1 (create-fertilizer-product! "F1" "F1" (hash n 10 p 20))) 209 (define f2 (create-fertilizer-product! "F2" "F2" (hash n 30 p 5))) 210 (define fertilizers (list f1 f2)) 211 212 (define matrix (get-fertilizer-product-matrix nutrients fertilizers)) 213 214 (check-= (matrix-ref matrix 0 0) 10 0 "N1 in F1") 215 (check-= (matrix-ref matrix 0 1) 30 0 "N1 in F2") 216 (check-= (matrix-ref matrix 1 0) 20 0 "N2 in F1") 217 (check-= (matrix-ref matrix 1 1) 5 0 "N2 in F2")) 218 219 (test-case "Single nutrient, single fertilizer" 220 (define A (matrix [[2]])) 221 (define y (col-matrix [10])) 222 (define ε 1e-6) 223 224 (define result (lawson-hanson-1974 A y ε)) 225 226 (check-= (matrix-ref result 0 0) 5.0 ε "Should give x = 5 since 2*5 = 10")) 227 228 (test-case "Two variables, known solution" 229 (define A (matrix [[1 0] [0 1]])) 230 (define y (col-matrix [3 4])) 231 (define ε 1e-6) 232 233 (define result (lawson-hanson-1974 A y ε)) 234 235 (check-= (matrix-ref result 0 0) 3.0 ε "x1 should be 3") 236 (check-= (matrix-ref result 1 0) 4.0 ε "x2 should be 4")) 237 238 (test-case "Overdetermined system" 239 (define A (matrix [[1 1] [2 1] [1 2]])) 240 (define y (col-matrix [3 5 5])) 241 (define ε 1e-4) 242 243 (define result (lawson-hanson-1974 A y ε)) 244 245 ;; Solution should be approximately [1.636, 1.636] (least squares fit) 246 (check-= (matrix-ref result 0 0) 1.636 0.01 "x1 approximately 1.636") 247 (check-= (matrix-ref result 1 0) 1.636 0.01 "x2 approximately 1.636")) 248 249 (test-case "Non-negativity enforcement" 250 (define A (matrix [[1 -1] [1 1]])) 251 (define y (col-matrix [1 3])) 252 (define ε 1e-6) 253 254 (define result (lawson-hanson-1974 A y ε)) 255 256 ;; All results should be non-negative 257 (check-true (>= (matrix-ref result 0 0) 0) "x1 should be non-negative") 258 (check-true (>= (matrix-ref result 1 0) 0) "x2 should be non-negative")) 259 260 (test-case "Zero target" 261 (define A (matrix [[1 2] [3 4]])) 262 (define y (col-matrix [0 0])) 263 (define ε 1e-6) 264 265 (define result (lawson-hanson-1974 A y ε)) 266 267 (check-= (matrix-ref result 0 0) 0.0 ε "x1 should be 0") 268 (check-= (matrix-ref result 1 0) 0.0 ε "x2 should be 0")) 269 270 (test-case "Ferti recipe" 271 (define test-date "2025-01-01") 272 (define n (get-nutrient #:name "Nitrogen")) 273 (define p (get-nutrient #:name "Phosphorus")) 274 (define k (get-nutrient #:name "Potassium")) 275 276 (create-nutrient-measurement! (nutrient-measurement #f test-date (hash n 0 p 0))) 277 (define test-crop-requirement (create-crop-requirement! "Test requirement" (hash n 100 p 50))) 278 (create-crop-rotation! test-date (hash test-crop-requirement 100)) 279 280 (create-fertilizer-product! "Nitrogen" "King Nitrogen" (hash n 100)) 281 (create-fertilizer-product! "Phosphorus" "Phosphorescent Baboon" (hash n 10 p 100)) 282 (create-fertilizer-product! "Diluted phosphorus" "John's Phosphorus" (hash n 3 p 30)) 283 284 (define recipe (find-ferti-recipe test-date)) 285 286 (check-equal? (hash-count recipe) 5 "Should have 5 fertilizer products") 287 288 (for ([(fertilizer quantity) (in-hash recipe)]) 289 (check-true (>= quantity 0) "Fertilizer quantity should be non-negative"))) 290 291 (test-case "Recipe calculation with real-world scenario" 292 (define test-date "2025-01-02") 293 (define n (get-nutrient #:name "Nitrogen")) 294 (define p (get-nutrient #:name "Phosphorus")) 295 (define k (get-nutrient #:name "Potassium")) 296 297 ;; Current levels 298 (create-nutrient-measurement! (nutrient-measurement #f test-date (hash n 50 p 10 k 100))) 299 300 ;; Target levels 301 (define test-crop-requirement 302 (create-crop-requirement! "Test requirement 2" (hash n 150 p 30 k 200))) 303 (create-crop-rotation! test-date (hash test-crop-requirement 100)) 304 305 ;; Available fertilizer products 306 (create-fertilizer-product! "" "Balanced" (hash n 100 p 100 k 100)) 307 (create-fertilizer-product! "Nitrogen blend" "High-N" (hash n 200 p 50 k 50)) 308 (create-fertilizer-product! "Phosphorus blend" "High-P" (hash n 50 p 200 k 50)) 309 (create-fertilizer-product! "Potassium blend" "High-K" (hash n 50 p 50 k 200)) 310 311 (define recipe (find-ferti-recipe test-date)) 312 313 (check-equal? (hash-count recipe) 9 "Recipe should have 9 fertilizers") 314 315 ;; Verify solution is non-negative 316 (for ([(fertilizer quantity) (in-hash recipe)]) 317 (check-true (>= quantity 0) 318 (format "~a quantity must be non-negative" 319 (fertilizer-product-name fertilizer)))))))) 320