731329f77f8b — Leonard Ritter 10 months ago
`* lebesgue integrals, varignon tests`
```4 files changed, 150 insertions(+), 120 deletions(-)

M lib/tukan/math.sc
A => testing/lebesgue_integral.sc
M testing/randproj_qt.sc
A => testing/varignon.sc
```
`M lib/tukan/math.sc +3 -1`
```@@ 35,9 35,11 @@ fn gaussian (x sigma)
a * (exp h)

fn wendland (x)
-    # produces a truncated gaussian similar to sigma 0.378
+    # approximates the right half of a truncated gaussian
# x must be in the range 0..1
+    # first derivative is zero at x = 0
# first three derivatives are zero at x = 1
+    # the integral of the interval 0..1 is 1/3
h := (1.0 - x)
h2 := h * h
h2 * h2 * (4.0 * x + 1.0)

```
`A => testing/lebesgue_integral.sc +46 -0`
```@@ 0,0 1,46 @@
+
+
+# do it in 100000 partitions
+let S = 100000
+
+# integral over begin..end
+let begin = 0.3
+let end = 1.1
+
+# f(x) = sqrt(x), f^-1(x) = x²
+fn f (x)
+    log x
+
+fn f^-1 (x)
+    exp x
+
+# integral of sqrt(x) over 0..1
+define integral
+    # our interval is begin..end
+    local k = 0.0
+    w := ((end - begin) / S)
+    for i in (range S)
+        x := begin + i as f32 * w
+        k += (f x) * w
+    k
+
+# lebesgue integral of sqrt(x)
+define lebesgue_integral
+    # our interval is f(begin)..f(end)
+    local k = 0.0
+    f_begin := (f begin)
+    f_end := (f end)
+    w := ((f_end - f_begin) / S)
+    # add constant part in front of interval
+    k += (end - begin) * f_begin
+    for i in (range S)
+        y := f_begin + i as f32 * w
+        k += (end - (f^-1 y)) * w
+    k
+
+
+print integral lebesgue_integral
+# > what is the difference between these two pictures?
+# they're the same!
+assert ((abs (integral - lebesgue_integral)) < 1e-5)
+

```
`M testing/randproj_qt.sc +42 -119`
```@@ 11,10 11,26 @@ using import tukan.bitmap
local rng : (Random)
'seed rng 712

-LEVELS := 5
+fn f (w c)
+    #(1.0 - (pow w 2.0)) ** c
+    #(1.0 - (pow w 0.125)) ** c
+    q := (pow (w * 4.0) 2.0)
+    pow 2.0 (- (q * c))
+
+#print
+    f (1.0 / 16.0) 4.0
+#print
+    f (2.0 / 16.0) 2.0
+#print
+    f (4.0 / 16.0) 1.0
+
+#if true
+    exit 0
+
+LEVELS := 6
Q := (1 << LEVELS)
N := Q * Q
-T := 3
+T := N // 2

fn weights (u v)
# produce the nine weights for a given offset of a 2x2 sized sample

@@ 51,8 67,8 @@ fn pick1 (levels level x y)

fn place1 (levels level x y v)
N := (2 << level)
-    if ((x < 0) | (y < 0) | (x >= N) | (y >= N))
-        return;
+    x := (x + N) % N
+    y := (y + N) % N
m := levels @ level
m @ (y * N + x) *= v
#print level "@" x y "+=" v

@@ 64,17 80,17 @@ fn place (levels x y)
(x as f32 + 0.5) / Qf; (y as f32 + 0.5) / Qf
for i in (rrange 0 LEVELS)
N := (2 << i) as f32
-        #w := (pow (N / Q) 0.25)
-        w := 1.0
+        w := (pow (N / Q) 2.0)
+        #w := 1.0
xf := u * N
yf := v * N
su := xf - (floor xf)
sv := yf - (floor yf)
x := (xf as i32)
y := (yf as i32)
-        print N x y su sv w
+        #print N x y su sv w
inline tf (x)
-            1.0 - x * 2.0 * w
+            1.0 - 2.0 * x
let p00 p10 p20 p01 p11 p21 p02 p12 p22 =
weights su sv
place1 levels i (x - 1) (y - 1) (tf p00)

@@ 90,13 106,16 @@ fn place (levels x y)
fn pick (levels rng)
let r = ('random rng)
fold (r x y = r 0 0) for i in (range 0 LEVELS)
+        let r = ('random rng)
x := x << 1
y := y << 1
-        v0 := (pick1 levels i x y)
-        v1 := (pick1 levels i x (y + 1))
-        v2 := (pick1 levels i (x + 1) y)
-        v3 := (pick1 levels i (x + 1) (y + 1))
-        #sum := 0.0
+        inline tf (x)
+            #x * x
+            x
+        v0 := (tf (pick1 levels i x y))
+        v1 := (tf (pick1 levels i x (y + 1)))
+        v2 := (tf (pick1 levels i (x + 1) y))
+        v3 := (tf (pick1 levels i (x + 1) (y + 1)))
sum := (+ v0 v1 v2 v3)
let v0 v1 v2 v3 sum =
if (sum == 0)

@@ 123,12 142,12 @@ fn pick (levels rng)
local level0 : (Array i32)
'resize level0 N

-#place levels 11 7
+place levels 32 32
for i in (range T)
let r x y = (pick levels rng)
-    print r x y
-    place levels x y
-    level0 @ (y * Q + x) = i
+    #print r x y
+    #place levels x y
+    level0 @ (y * Q + x) = (i + 1)

fn... ansicolor (col : vec3)
let r g b = (unpack (ivec3 (* 255.0 (clamp col (vec3 0) (vec3 1)))))

@@ 137,7 156,7 @@ fn... ansicolor (col : vec3)
b := (tostring b)
.. "\x1b[38;2;" r ";" g ";" b "m"

-do
+if true
for y in (range Q)
fy := y / Q
for x in (range Q)

@@ 151,22 170,22 @@ do
print (ansicolor (vec3 1)) "\n"

do
-    L := 4
-    LQ := 2 << L
-    T := 1.0
-    for y in (range LQ)
-        for x in (range LQ)
+    for y in (range Q)
+        for x in (range Q)
vvv bind w
fold (w = 1.0) for i in (range LEVELS)
+                #w := 1.0; i := 0
W := 2 << i
level := levels @ i
x := x >> (LEVELS - i - 1)
y := y >> (LEVELS - i - 1)
n := y * W + x
#w * (1.0 - (level @ n) * 2.0)
+                #w * (pow 2.0 (- (level @ n)))
w * (level @ n)
#w := (viridis (w / ((1 * LEVELS) as f32)))
#w := (clamp w 0.0 1.0)
+            #w := (viridis (w * T * T * T * T))
w := (viridis w)
io-write!
ansicolor w

@@ 178,105 197,9 @@ do
if true
exit 0

-local weights : (Array f32)
-'resize weights N 1.0
-
-fn coord (n)
-    _ (n % Q) (n // Q)
fn index (x y)
((y + Q) % Q) * Q + ((x + Q) % Q)

-fn moddist (a b q)
-    let d = (abs (a - b))
-    min d (q - d)
-
-local totalcount = 0.0
-phi := 2.0 / ((sqrt 5.0) + 1.0)
-for i in (range T)
-    if ((i % Q) == 0)
-        print "progress:" (i * 100 / T) "%"
-    let slots = ((countof level0) as i32)
-    #let k = ('random rng)
-    #let k = ((0.5 + phi * (i as f32)) % 1.0)
-    let k =
-        #'random rng
-        do
-            let g = 1.32471795724474602596
-            let a1 = (1.0 / g)
-            let a2 = (1.0 / (g * g))
-            /
-                index
-                    (((0.5 + a1 * (i as f32)) % 1.0) * Q) as i32
-                    (((0.5 + a2 * (i as f32)) % 1.0) * Q) as i32
-                N
-    local totalweight = 0.0
-    local maxweight = 0.0
-    for w in weights
-        totalweight += w
-        maxweight = (max maxweight w)
-    #if (i == (T - 1))
-        for y in (range Q)
-            for x in (range Q)
-                n := y * Q + x
-                let w = ((weights @ n) / maxweight)
-                io-write!
-                    ansicolor (vec3 w)
-                io-write! "██"
-            io-write! "\n"
-        if true
-            exit 0
-
-    fold (w = 0.0) for n in (range slots)
-        let inthis = (weights @ n)
-        let prob = (inthis / totalweight)
-        #let prob = (1.0 / N)
-        let w = (w + prob)
-        if (k < w)
-            level0 @ n = i
-            let cx cy = (coord n)
-
-            let maxlen = (length (vec2 Q Q))
-
-            for i u in (enumerate weights)
-                let x y = (coord i)
-                let dx = ((moddist x cx Q) / Q)
-                let dy = ((moddist y cy Q) / Q)
-                let d = (length (vec2 dx dy))
-                #let phi = (sqrt 0.04)
-                #let phi = (sqrt 0.001)
-                #let phi = (sqrt 0.1)
-                #let phi = (sqrt 0.01)
-                let g =
-                    /
-                        exp
-                            /
-                                - (pow d 2.0)
-                                2.0 * phi * phi
-                        1.0 #phi * (sqrt tau)
-                let w = (1.0 - g)
-                u = (u / maxweight) * w
-                #u *= w
-            break w
-        w
-    totalcount += 1.0
-#for y in (range Q)
-    fy := y / Q
-    for x in (range Q)
-        n := y * Q + x
-        let w = ((level0 @ n) / T)
-        io-write!
-            ansicolor w
-            #ansicolor (step fy w)
-        io-write! "██"
-    io-write! "\n"
-#print (ansicolor (vec3 1)) "\n"
-#fn... ansicolor (col : vec3)
-    let r g b = (unpack (ivec3 (* 255.0 (clamp col (vec3 0) (vec3 1)))))
-    r := (tostring r)
-    g := (tostring g)
-    b := (tostring b)
-    .. "\x1b[38;2;" r ";" g ";" b "m"
-
let outimage =
Bitmap4 (ivec2 Q Q)
for x y in (dim Q Q)

```
`A => testing/varignon.sc +59 -0`
```@@ 0,0 1,59 @@
+
+# experiments with varignon's theorem
+
+using import glm
+
+#embed
+    p0 := (vec3 0 0 0)
+    p1 := (vec3 2 0 0)
+    p2 := (vec3 0 2 0)
+    p3 := (vec3 0 0 2)
+
+embed
+    p0 := (vec3 0 0 0) + 1
+    p1 := (vec3 1 3 0) + 1
+    p2 := (vec3 0 6 0) + 1
+    p3 := (vec3 -1 3 0) + 1
+
+p0123 := (p0 + p1 + p2 + p3) * 0.25
+
+p01 := (p0 + p1) * 0.5
+p02 := (p0 + p2) * 0.5
+p03 := (p0 + p3) * 0.5
+p12 := (p1 + p2) * 0.5
+p13 := (p1 + p3) * 0.5
+p23 := (p2 + p3) * 0.5
+p30 := (p3 + p0) * 0.5
+
+gx := (p23 - p01)
+gy := (p13 - p02)
+gz := (p12 - p03)
+
+print gx gy gz
+
+p012 := p12 - p01
+p123 := p23 - p12
+p032 := p23 - p30
+p103 := p30 - p01
+
+print p012
+print p123
+print p032
+print p103
+
+u := p2 - p0
+v := p3 - p1
+
+print u v
+print (degrees (acos (dot (normalize u) (normalize v))))
+
+
+n := (cross u v)
+q := (length n)
+n := (n / q)
+
+print (dot n p0123)
+print (dot n p01)
+print (dot n p12)
+print (dot n p23)
+print (dot n p30)
No newline at end of file

```