From d3f20a07516069f5bb7921650299026af58df52d Mon Sep 17 00:00:00 2001
From: Wesley Kuhron Jones <jazzzman627@yahoo.com>
Date: Mon, 24 Oct 2022 00:01:12 +0000
Subject: [PATCH] playing with a chaotic map

---
 Chaos/ShearMap.py | 67 +++++++++++++++++++++++++++++++++++++++++++++++
 PltContentOnly.py |  4 +--
 2 files changed, 69 insertions(+), 2 deletions(-)
 create mode 100644 Chaos/ShearMap.py

diff --git a/Chaos/ShearMap.py b/Chaos/ShearMap.py
new file mode 100644
index 00000000..da49d774
--- /dev/null
+++ b/Chaos/ShearMap.py
@@ -0,0 +1,67 @@
+# map square [0,1]^2 onto itself with shear and cut function
+
+import random
+import numpy as np
+import matplotlib.pyplot as plt
+from datetime import datetime
+
+import sys
+sys.path.insert(0, "/home/wesley/programming/")
+from PltContentOnly import imshow_content_only
+
+
+get_signed_pareto = lambda: random.choice([-1, 1]) * (np.random.pareto(1))
+get_normal = lambda: np.random.normal()
+
+
+def get_polynomial_func(degree, mod=True):
+    coeffs = [get_signed_pareto() for i in range(degree+1)]
+    print(coeffs)
+    f = lambda x, coeffs=coeffs: sum(coeffs[i] * x**i for i in range(degree+1))
+    if mod:
+        return lambda x: f(x) % 1
+    else:
+        return f
+
+
+def get_bivariate_polynomial_func(degree):
+    fxx = get_polynomial_func(degree, mod=False)
+    fxy = get_polynomial_func(degree, mod=False)
+    fyx = get_polynomial_func(degree, mod=False)
+    fyy = get_polynomial_func(degree, mod=False)
+    f = lambda x, y, fxx=fxx, fxy=fxy, fyx=fyx, fyy=fyy: ((fxx(x) + fxy(y)) % 1, (fyx(x) + fyy(y)) % 1)
+    return np.vectorize(f)
+
+
+def run():
+    f = get_bivariate_polynomial_func(2)
+
+    xs = np.linspace(0, 1, 100)
+    ys = np.linspace(0, 1, 100)
+    X, Y = np.meshgrid(xs, ys)
+    X = X.flatten()
+    Y = Y.flatten()
+
+    heatmap_resolution = 250
+    bin_counts_2d = np.zeros((heatmap_resolution, heatmap_resolution))
+
+    n_steps = 1000
+    for i in range(n_steps):
+        if i % 10 == 0:
+            print(f"step {i}/{n_steps}")
+        X, Y = f(X, Y)
+        for x, y in zip(X, Y):
+            # find what bin this point is in
+            x_bin_i = int(x / (1 / heatmap_resolution))
+            y_bin_i = int(y / (1 / heatmap_resolution))
+            bin_counts_2d[x_bin_i, y_bin_i] += 1
+
+    now_str = datetime.utcnow().strftime("%Y-%m-%d-%H:%M:%S")
+    save_fp = f"Images/ShearMap/{now_str}.png"
+    imshow_content_only(bin_counts_2d.T, size_inches=(10,10), save_fp=save_fp, origin="lower")
+    plt.gcf().clear()
+
+
+if __name__ == "__main__":
+    while True:
+        run()
diff --git a/PltContentOnly.py b/PltContentOnly.py
index 42690e6b..c80c578e 100644
--- a/PltContentOnly.py
+++ b/PltContentOnly.py
@@ -13,10 +13,10 @@ def get_fig_ax(size_inches):
     return fig, ax
 
 
-def imshow_content_only(arr, size_inches=(6,6), save_fp=None):
+def imshow_content_only(arr, size_inches=(6,6), save_fp=None, **kwargs):
     # https://stackoverflow.com/questions/8218608/scipy-savefig-without-frames-axes-only-content
     fig, ax = get_fig_ax(size_inches)
-    ax.imshow(arr, aspect="auto")
+    ax.imshow(arr, aspect="auto", **kwargs)
     if save_fp is not None:
         fig.savefig(save_fp)