diff options
Diffstat (limited to 'services/nnls.rkt')
| -rw-r--r-- | services/nnls.rkt | 138 |
1 files changed, 57 insertions, 81 deletions
diff --git a/services/nnls.rkt b/services/nnls.rkt index 649d416..96d37ec 100644 --- a/services/nnls.rkt +++ b/services/nnls.rkt @@ -20,23 +20,16 @@ (define nutrients (get-nutrients)) (define fertilizer-product-matrix (get-fertilizer-product-matrix nutrients fertilizers)) (define deficits - (->col-matrix - (for/list ([n nutrients]) - (define latest-measurement (get-latest-nutrient-measurement-value n)) - (define latest-target (get-latest-nutrient-target-value n)) - (define deficit - (cond - [(false? latest-target) - 0] - [(or (false? latest-measurement) - (zero? latest-measurement)) - latest-target] - [(and (number? latest-measurement) - (number? latest-target)) - (* 100 - (/ (- latest-target latest-measurement) - latest-measurement))])) - deficit))) + (->col-matrix (for/list ([n nutrients]) + (define latest-measurement (get-latest-nutrient-measurement-value n)) + (define latest-target (get-latest-nutrient-target-value n)) + (define deficit + (cond + [(false? latest-target) 0] + [(or (false? latest-measurement) (zero? latest-measurement)) latest-target] + [(and (number? latest-measurement) (number? latest-target)) + (* 100 (/ (- latest-target latest-measurement) latest-measurement))])) + deficit))) (define error-threshold 10e-4) (lawson-hanson-1974 fertilizer-product-matrix deficits error-threshold)) @@ -46,16 +39,14 @@ ;; Inputs ;;;;;;;;; - (-> matrix? ; Real-valued matrix A of dimension m × n + (-> matrix? ; Real-valued matrix A of dimension m × n col-matrix? ; Real-valued column matrix (vector) y of dimension m - real? ; Real-value ɛ, tolerance for the stopping criterion - col-matrix?) ; Real-valued solution column matrix x - (define-values - (m ; Number of nutrients - n) ; Number of fertilizer products + real? ; Real-value ɛ, tolerance for the stopping criterion + col-matrix?) ; Real-valued solution column matrix x + (define-values (m ; Number of nutrients + n) ; Number of fertilizer products (matrix-shape A)) - ;;;;;;;;;;;;; ;; Initialize ;;;;;;;;;;;;; @@ -70,12 +61,12 @@ ;; Gradient-like vector for residual error. (define (compute-w x) - (matrix* (matrix-transpose A) - (matrix- y (matrix* A x)))) + (matrix* (matrix-transpose A) (matrix- y (matrix* A x)))) ;; max over j in R of w_j, returning (values max-val j*) (define (max-w-in-R w) - (for/fold ([max-val -inf.0] [max-j #f]) + (for/fold ([max-val -inf.0] + [max-j #f]) ([j (in-set R)]) (define v (colv-ref w j)) (if (> v max-val) @@ -88,24 +79,25 @@ (if (set-empty? P) (make-matrix n 1 0) (let* ([idxs (sort (set->list P) <)] - [AP (submatrix A (::) idxs)] - [sP (matrix* (matrix-inverse - (matrix* (matrix-transpose AP) AP)) - (matrix-transpose AP) - y)]) + [AP (submatrix A (::) idxs)] + [sP (matrix* (matrix-inverse (matrix* (matrix-transpose AP) AP)) + (matrix-transpose AP) + y)]) ;; map: column index j in P -> corresponding sP entry (define mapping - (for/list ([j idxs] [k (in-naturals)]) + (for/list ([j idxs] + [k (in-naturals)]) (cons j (colv-ref sP k)))) (define (s-at i) (define p (assoc i mapping)) - (if p (cdr p) 0)) + (if p + (cdr p) + 0)) (build-matrix n 1 (λ (i j) (s-at i)))))) ;; The "first try" x represents no addition of any fertilizer. (define x (make-matrix n 1 0)) - ;;;;;;;;;;;;; ;; Outer loop ;;;;;;;;;;;;; @@ -115,16 +107,14 @@ (cond ;; If no remaining candidates in R, we're done. - [(set-empty? R) - x] + [(set-empty? R) x] [else (define-values (max-val j*) (max-w-in-R w)) ;; Stopping criterion: max(w_R) <= ε (cond - [(or (not j*) (<= max-val ε)) - x] + [(or (not j*) (<= max-val ε)) x] [else ;; Add j* to P, remove from R @@ -139,8 +129,7 @@ (define min-sP (if (set-empty? P) +inf.0 - (for/fold ([mn +inf.0]) - ([j (in-set P)]) + (for/fold ([mn +inf.0]) ([j (in-set P)]) (min mn (colv-ref s j))))) (cond @@ -152,11 +141,10 @@ [else ;; Compute α = min_{i in P, s_i <= 0} x_i / (x_i - s_i) (define α - (for/fold ([a +inf.0]) - ([j (in-set P)]) + (for/fold ([a +inf.0]) ([j (in-set P)]) (define sj (colv-ref s j)) (if (<= sj 0) - (let* ([xj (colv-ref x j)] + (let* ([xj (colv-ref x j)] [den (- xj sj)]) (if (> den 0) (min a (/ xj den)) @@ -167,8 +155,7 @@ (error 'lawson-hanson-1974 "no valid α in inner loop")) ;; x ← x + α (s − x) - (define new-x - (matrix+ x (matrix* α (matrix- s x)))) + (define new-x (matrix+ x (matrix* α (matrix- s x)))) ;; Move to R all indices j in P with x_j <= 0 (define to-remove '()) @@ -189,9 +176,10 @@ (λ (i j) (define selected-nutrient (list-ref nutrients i)) (define product (list-ref fertilizers j)) - (define pair (assoc selected-nutrient - (fertilizer-product-values product))) - (if pair (cdr pair) 0)))) + (define pair (assoc selected-nutrient (fertilizer-product-values product))) + (if pair + (cdr pair) + 0)))) (module+ test (require rackunit @@ -202,35 +190,23 @@ (define test-date "2025-01-01") (run-tests - (test-suite - "Nutrient measurement model" - #:before (λ () - (connect! #:path 'memory) - ;; (connect! #:path "test.sqlite3") - (migrate-all!) - - (define nitrogen (create-nutrient! "Nitrogen" "N")) - (define phosphorus (create-nutrient! "Phosphorus" "P")) - - (create-nutrient-measurement! test-date - `((,nitrogen . 0) - (,phosphorus . 0))) - (create-nutrient-target! test-date - `((,nitrogen . 100) - (,phosphorus . 50))) - - (create-fertilizer-product! "King Nitrogen" - `((,nitrogen . 100))) - (create-fertilizer-product! "Phosphorescent Baboon" - `((,nitrogen . 10) - (,phosphorus . 100))) - (create-fertilizer-product! "John's Phosphorus" - `((,nitrogen . 3) - (,phosphorus . 30)))) - #:after (λ () - (disconnect!)) - - (test-case "Solve for NNLS" - (displayln - (format "Final solution for fertilizers is combination ~a" - (find-ferti-recipe))))))) + (test-suite "Nutrient measurement model" + #:before + (λ () + (connect! #:path 'memory) + ;; (connect! #:path "test.sqlite3") + (migrate-all!) + + (define nitrogen (create-nutrient! "Nitrogen" "N")) + (define phosphorus (create-nutrient! "Phosphorus" "P")) + + (create-nutrient-measurement! test-date `((,nitrogen . 0) (,phosphorus . 0))) + (create-nutrient-target! test-date `((,nitrogen . 100) (,phosphorus . 50))) + + (create-fertilizer-product! "King Nitrogen" `((,nitrogen . 100))) + (create-fertilizer-product! "Phosphorescent Baboon" `((,nitrogen . 10) (,phosphorus . 100))) + (create-fertilizer-product! "John's Phosphorus" `((,nitrogen . 3) (,phosphorus . 30)))) + #:after (λ () (disconnect!)) + + (test-case "Solve for NNLS" + (displayln (format "Final solution for fertilizers is combination ~a" (find-ferti-recipe))))))) |