# HG changeset patch # User Leonard Ritter # Date 1635435568 -7200 # Thu Oct 28 17:39:28 2021 +0200 # Node ID 731329f77f8b625f4facee2e1f14933fcd2450fd # Parent 1f76611de99d62ba6637d9dcb75bf25da0479bc0 * lebesgue integrals, varignon tests diff --git a/lib/tukan/math.sc b/lib/tukan/math.sc --- a/lib/tukan/math.sc +++ b/lib/tukan/math.sc @@ -35,9 +35,11 @@ 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) diff --git a/testing/lebesgue_integral.sc b/testing/lebesgue_integral.sc new file mode 100644 --- /dev/null +++ b/testing/lebesgue_integral.sc @@ -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) + diff --git a/testing/randproj_qt.sc b/testing/randproj_qt.sc --- a/testing/randproj_qt.sc +++ b/testing/randproj_qt.sc @@ -11,10 +11,26 @@ 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 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 @@ (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 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 @@ 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 @@ 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 @@ 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 @@ 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) diff --git a/testing/varignon.sc b/testing/varignon.sc new file mode 100644 --- /dev/null +++ b/testing/varignon.sc @@ -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