some development to refactor the entire code to use funcanimation and…

… a more oop approach to building the stratigraphy

LICENSE.txt

@@ -2,20 +2,8 @@ MIT License
 
 Copyright (c) 2018 Andrew Moodie
 
-Permission is hereby granted, free of charge, to any person obtaining a copy
-of this software and associated documentation files (the "Software"), to deal
-in the Software without restriction, including without limitation the rights
-to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-copies of the Software, and to permit persons to whom the Software is
-furnished to do so, subject to the following conditions:
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 
-The above copyright notice and this permission notice shall be included in all
-copies or substantial portions of the Software.
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
 
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-SOFTWARE.
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

src/bayes.py

@@ -0,0 +1,56 @@
+import math
+
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.animation import FuncAnimation
+
+
+def beta_pdf(x, a, b):
+    return (x**(a-1) * (1-x)**(b-1) * math.gamma(a + b)
+            / (math.gamma(a) * math.gamma(b)))
+
+
+class UpdateDist(object):
+    def __init__(self, ax, prob=0.5):
+        self.success = 0
+        self.prob = prob
+        self.line, = ax.plot([], [], 'k-')
+        self.x = np.linspace(0, 1, 200)
+        self.ax = ax
+
+        # Set up plot parameters
+        self.ax.set_xlim(0, 1)
+        self.ax.set_ylim(0, 15)
+        self.ax.grid(True)
+
+        # This vertical line represents the theoretical value, to
+        # which the plotted distribution should converge.
+        self.ax.axvline(prob, linestyle='--', color='black')
+
+    def init(self):
+        self.success = 0
+        self.line.set_data([], [])
+        return self.line,
+
+    def __call__(self, i):
+        # This way the plot can continuously run and we just keep
+        # watching new realizations of the process
+        if i == 0:
+            return self.init()
+
+        # Choose success based on exceed a threshold with a uniform pick
+        if np.random.rand(1,) < self.prob:
+            self.success += 1
+        y = beta_pdf(self.x, self.success + 1, (i - self.success) + 1)
+        self.line.set_data(self.x, y)
+        return self.line,
+
+# Fixing random state for reproducibility
+np.random.seed(19680801)
+
+
+fig, ax = plt.subplots()
+ud = UpdateDist(ax, prob=0.7)
+anim = FuncAnimation(fig, ud, frames=np.arange(100), init_func=ud.init,
+                     interval=100, blit=True)
+plt.show()

src/rain.py

@@ -0,0 +1,60 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.animation import FuncAnimation
+
+# Fixing random state for reproducibility
+np.random.seed(19680801)
+
+
+# Create new Figure and an Axes which fills it.
+fig = plt.figure(figsize=(7, 7))
+ax = fig.add_axes([0, 0, 1, 1], frameon=False)
+ax.set_xlim(0, 1), ax.set_xticks([])
+ax.set_ylim(0, 1), ax.set_yticks([])
+
+# Create rain data
+n_drops = 50
+rain_drops = np.zeros(n_drops, dtype=[('position', float, 2),
+                                      ('size',     float, 1),
+                                      ('growth',   float, 1),
+                                      ('color',    float, 4)])
+
+# Initialize the raindrops in random positions and with
+# random growth rates.
+rain_drops['position'] = np.random.uniform(0, 1, (n_drops, 2))
+rain_drops['growth'] = np.random.uniform(50, 200, n_drops)
+
+# Construct the scatter which we will update during animation
+# as the raindrops develop.
+scat = ax.scatter(rain_drops['position'][:, 0], rain_drops['position'][:, 1],
+                  s=rain_drops['size'], lw=0.5, edgecolors=rain_drops['color'],
+                  facecolors='none')
+
+
+def update(frame_number):
+    # Get an index which we can use to re-spawn the oldest raindrop.
+    current_index = frame_number % n_drops
+
+    # Make all colors more transparent as time progresses.
+    rain_drops['color'][:, 3] -= 1.0/len(rain_drops)
+    rain_drops['color'][:, 3] = np.clip(rain_drops['color'][:, 3], 0, 1)
+
+    # Make all circles bigger.
+    rain_drops['size'] += rain_drops['growth']
+
+    # Pick a new position for oldest rain drop, resetting its size,
+    # color and growth factor.
+    rain_drops['position'][current_index] = np.random.uniform(0, 1, 2)
+    rain_drops['size'][current_index] = 5
+    rain_drops['color'][current_index] = (0, 0, 0, 1)
+    rain_drops['growth'][current_index] = np.random.uniform(50, 200)
+
+    # Update the scatter collection, with the new colors, sizes and positions.
+    scat.set_edgecolors(rain_drops['color'])
+    scat.set_sizes(rain_drops['size'])
+    scat.set_offsets(rain_drops['position'])
+
+
+# Construct the animation, using the update function as the animation director.
+animation = FuncAnimation(fig, update, interval=10)
+plt.show()