```#!/usr/bin/env sh
# -*- wisp -*-
guile -L \$(dirname \$(dirname \$(realpath "\$0"))) -c '(import (language wisp spec))'
exec guile -L \$(dirname \$(dirname \$(realpath "\$0"))) --language=wisp -e '(@@ (examples cholesky) main)' -s "\$0" "\$@"
; !#

;; Cholesky decomposition, following https://de.wikipedia.org/wiki/Cholesky-Zerlegung#Pseudocode

define-module : examples cholesky
. #:export : cholesky! matrix-ref matrix-set! matrix-transpose matrix-multiply

use-modules : srfi srfi-42 ; list-comprehension
srfi srfi-11 ; let-values

define : ->exact-matrix list-of-lists
. "Turn a list of lists into a matrix"
map
λ : x
apply list : map inexact->exact x
. list-of-lists

define : ->inexact-matrix list-of-lists
. "Turn a list of lists into a matrix"
map
λ : x
apply list : map exact->inexact x
. list-of-lists

define : matrix-ref X row col
list-ref (list-ref X row) col

define : matrix-set! X row col val
let : : sublist : list-ref X row
list-set! sublist col val
list-set! X row sublist

define : matrix-transpose X
. "Swap columns and rows of a matrix"
list-ec (: outer (length X)) ; outer
list-ec (: inner (length (list-ref X outer))) ; inner
matrix-ref X inner outer

define : matrix-multiply X Y
. "Calculate the matrix product of X and Y"
list-ec (: row (length Y))
list-ec (: col (length X))
sum-ec (: inner (length (list-ref Y row)))
* : matrix-ref Y inner col
matrix-ref X row inner

define : mostly-exact-sqrt n
. "Calculate an exact sqrt if possible, else use a good approximation"
inexact->exact : sqrt : inexact->exact n

define : cholesky! a
. "Modifies the square matrix a to contain its cholesky decomposition.

sets a to g with a = ggT,

a is represented as list of lists."
let : : n : length a
do-ec (: i n)
do-ec (: j (+ 1 i))
let : : sum : matrix-ref a i j
; format #t "n: ~A i: ~A j: ~A\n" n i j
when : >= j 1
do-ec (: k j)
set! sum : - sum : * (matrix-ref a i k) (matrix-ref a j k)
cond
: > i j ; lower triangle
matrix-set! a i j
/ sum
matrix-ref a j j
: > sum 0 ; diagonal element
matrix-set! a i i
mostly-exact-sqrt sum ; preserves the exactness, since the result is an exact number, though not an exact result
else
throw 'matrix-numerically-not-symmetric-positive-definite a
do-ec (: i n)
do-ec (: j (+ 1 i) n)
matrix-set! a i j 0
. a

define : main args
let
:  X : ->exact-matrix '(( 1  -1   1)
(-1   3 -.5)
( 1 -.5   4))
L : ->exact-matrix '(( 1 0          0)
(-1 1.41421356 0)
( 1 0.35355339 1.6955825))
format #t "X\n"
display : ->inexact-matrix X
newline
format #t "cholesky\n"
display : ->inexact-matrix : cholesky! X
newline
format #t "L\n"
display : ->inexact-matrix L
newline
format #t "L·Lt\n"
display
->inexact-matrix
matrix-multiply L : matrix-transpose L
newline
format #t "X·Xt\n"
display
->inexact-matrix
matrix-multiply X : matrix-transpose X
newline
```