# Author: Damodar Rajbhandari (2023)

%reset -f

%matplotlib notebook

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc
from matplotlib.animation import FuncAnimation
from IPython.display import clear_output
from IPython.display import Image

plt.rcParams.update({
    "text.usetex": True,
    "font.family": "Computer Modern Roman",
    "font.size": 24,
})

n_dots = 1_000_000  # total number of dots
k = np.linspace(0.0, 4.0, n_dots)
iterations = 100

fig = plt.figure(figsize=(20, 15))
ax = plt.axes()
ax.set_xlim(0.0, 4.0)
ax.set_ylim(0.0, 1.0)

rc('animation', html='html5')

w = 1.5  # \omega bias setting to 1.5

def animate(i):
  plt.cla()
  ax.set_title(r"Feigenbaum Tree with bias $\omega$; $x^{(k)}_{n + 1} = \frac{k x^{(k)}_{n}(1 - x^{(k)}_{n}) + \omega x^{(k)}_{n}}{1 + \omega}$ for $\omega = 1.5$, $n= $"+" {}".format(i +1), fontsize=24)
  x = np.random.uniform(0.0, 1.0, n_dots)
  for _ in range(i + 1): 
    x = (k * x * (1 - x) + w * x)/(1 + w)
  
  ax.plot(k, x, '.k', markersize=.04)
  ax.set(xlabel=r'$k$', ylabel=r'$x^{(k)}_{n}$')  

  print("Done iteration = {}".format(i))
  clear_output(wait=True)

fig = FuncAnimation(fig, animate, frames=iterations, interval=300, blit=True)
fig.save('feigenbaum_tree_with_bias.gif', writer='imagemagick', fps=10)

with open('feigenbaum_tree_with_bias.gif','rb') as f:
    display(Image(data=f.read(), format='png'))