import numpy as np import matplotlib.pyplot as plt import scipy

import tempfile import os import imageio

def anscombe_transform(samples, m):

   return 2 * np.sqrt(samples + 3/8) - (2*np.sqrt(m+3/8) - 1/(4*np.sqrt(m)))

def plot_anscombe(m=10, n_samples=1000000):

   fig, axes = plt.subplot_mosaic("A", figsize=(8, 4))
   ax1 = axes["A"]
   samples = anscombe_transform(np.random.poisson(m, n_samples), m)
   mean_diff = np.mean(samples)
   bins = sorted(list(set(samples)))
   # Plot the histogram of the samples
   ax1.hist(samples, bins=bins, align='right', rwidth=2, density=True)
   xs = np.linspace(-3.5, 3.5, 1000)
   ax1.plot(xs, scipy.stats.norm.pdf(xs))
   ax1.vlines([mean_diff], 0,0.4, color='k')
   # Set the x-axis label and title
   ax1.set_xlabel('Number of Events')
   ax1.set_xlim(-4,+4)
   ax1.set_ylim(0, 0.44)
   ax1.set_title('Anscombe transform of Poisson(m)')
   
   text_lines = [r'$m =$' + f'{m}',
                 r'$m^{3/2}\mu =$' + f'{m**1.5 * mean_diff:.2f}, ', 
                 r'$m^{2}(\sigma-1) =$' + f'{m**2 * (np.std(samples)-1):.2f}',]
   text_x = 0.03
   text_y = 0.9
   text_color = 'black'
   text_size = 12
   for i, line in enumerate(text_lines):
       ax1.text(text_x, text_y-(i*0.08), line, 
                color=text_color, fontsize=text_size, 
                ha='left', va='bottom', transform=ax1.transAxes)
   fig.tight_layout()
   return fig

def interpolate_counts(counts, frames_per_step):

   interpolated_counts = [counts[0]]
   for i in range(1,len(counts)):
       interval = (counts[i] - counts[i-1]) // i
       interpolated_counts += list(range(counts[i-1], counts[i], interval))
   return interpolated_counts + [counts[-1]]

with tempfile.TemporaryDirectory() as temp_dir:

   n_steps = 10
   frames_per_step = 10
   ms = interpolate_counts([2**n for n in range(n_steps)], frames_per_step)
   n_frames = len(ms)-1
   
   for i in range(n_frames):
       fig = plot_anscombe(m=ms[i], n_samples=10000000)
       filename = os.path.join(temp_dir, f"plot_{i:03d}.png")
       fig.savefig(filename)
       plt.close(fig)

   # Compile images into GIF
   fps = 12
   images = []
   for i in range(n_frames):
       filename = os.path.join(temp_dir, f"plot_{i:03d}.png")
       images.append(imageio.imread(filename))
   imageio.mimsave(f"Anscombe transform.gif", images, duration=1/fps)