import stats;
import graph;
import palette;
import contour;

size(20cm);

scale(false);

pair[] data=new pair[50000];
for(int i=0; i < data.length; ++i)
  data[i]=Gaussrandpair();

// Histogram limits and number of bins
pair datamin=(-0.15,-0.15);
pair datamax=(0.15,0.15);
int Nx=30;
int Ny=30;

int[][] bins=frequency(data,datamin,datamax,Nx,Ny);

real[] values=new real[Nx*Ny];
pair[] points=new pair[Nx*Ny];
int k=0;
real dx=(datamax.x-datamin.x)/Nx;
real dy=(datamax.y-datamin.y)/Ny;
for(int i=0; i < Nx; ++i) {
  for(int j=0; j < Ny; ++j) {
    values[k]=bins[i][j];
    points[k]=(datamin.x+(i+0.5)*dx,datamin.y+(j+0.5)*dy);
    ++k;
  }
}

// Create a color palette
pen[] InvGrayscale(int NColors=256) {
  real ninv=1.0/(NColors-1.0);
  return sequence(new pen(int i) {return gray(1-17*i*ninv);},NColors);
}

// Draw the histogram, with axes
bounds range=image(points,values,Range(0,40),InvGrayscale());
draw(contour(points,values,new real[] {1,2,3,4,8,12,16,20,24,28,32,36,40},
             operator--),blue);
xaxis("$x$",BottomTop,LeftTicks,above=true);
yaxis("$y$",LeftRight,RightTicks,above=true);