import timeit
import numpy as np
import cosmotool as ct

def myfun(N):
  f=0.10
  d=np.random.normal(size=(N,)*3)
  rho = d + f *(d*d - np.average(d*d))
  
  B = ct.bispectrum(rho, 1, N, fourier=False)
  P = ct.powerspectrum(rho, 1, N, fourier=False)
  PP = P[1]/P[0]/N**3
  
  x = PP[:,None,None] * PP[None,:,None] + PP[:,None,None]*PP[None,None,:] + PP[None,:,None]*PP[None,None,:]
  
  BB = B[1]/B[0]/N**6
  
  y = BB/x

  np.savez("bispec_%d.npy" % N, y=y, B=B, P=P, rho=rho);
  
  
print( timeit.timeit('from __main__ import myfun; myfun(16)', number=1) )
print( timeit.timeit('from __main__ import myfun; myfun(32)', number=1) )
print( timeit.timeit('from __main__ import myfun; myfun(64)', number=1) )