private import graph;

real legendmarkersize=2mm;

real mean(real A[])
{
  return sum(A)/A.length;
}

// unbiased estimate
real variance(real A[])
{
  return sum((A-mean(A))^2)/(A.length-1);
}

real variancebiased(real A[])
{
  return sum((A-mean(A))^2)/A.length;
}

// unbiased estimate
real stdev(real A[])
{
  return sqrt(variance(A));
}

real rms(real A[])
{
  return sqrt(sum(A^2)/A.length);
}

real skewness(real A[])
{
  real[] diff=A-mean(A);
  return sum(diff^3)/sqrt(sum(diff^2)^3/A.length);
}

real kurtosis(real A[])
{
  real[] diff=A-mean(A);
  return sum(diff^4)/sum(diff^2)^2*A.length;
}

real kurtosisexcess(real A[])
{
  return kurtosis(A)-3;
}

real Gaussian(real x, real sigma)
{
  static real sqrt2pi=sqrt(2pi);
  return exp(-0.5*(x/sigma)^2)/(sigma*sqrt2pi);
}

real Gaussian(real x)
{
  static real invsqrt2pi=1/sqrt(2pi);
  return exp(-0.5*x^2)*invsqrt2pi;
}

// Return frequency count of data in [bins[i],bins[i+1]) for i=0,...,n-1.
int[] frequency(real[] data, real[] bins)
{
  int n=bins.length-1;
  int[] freq=new int[n];
  for(int i=0; i < n; ++i)
    freq[i]=sum(bins[i] <= data & data < bins[i+1]);
  return freq;
}

// Return frequency count in n uniform bins from a to b
// (faster than the above more general algorithm).
int[] frequency(real[] data, real a, real b, int n)
{
  int[] freq=sequence(new int(int x) {return 0;},n);
  real h=n/(b-a);
  for(int i=0; i < data.length; ++i) {
    int I=Floor((data[i]-a)*h);
    if(I >= 0 && I < n)
      ++freq[I];
  }
  return freq;
}

// Return frequency count in [xbins[i],xbins[i+1]) and [ybins[j],ybins[j+1]).
int[][] frequency(real[] x, real[] y, real[] xbins, real[] ybins)
{
  int n=xbins.length-1;
  int m=ybins.length-1;
  int[][] freq=new int[n][m];
  bool[][] inybin=new bool[m][y.length];
  for(int j=0; j < m; ++j)
    inybin[j]=ybins[j] <= y & y < ybins[j+1];
  for(int i=0; i < n; ++i) {
    bool[] inxbini=xbins[i] <= x & x < xbins[i+1];
    int[] freqi=freq[i];
    for(int j=0; j < m; ++j)
      freqi[j]=sum(inxbini & inybin[j]);
  }
  return freq;
}

// Return frequency count in nx by ny uniform bins in box(a,b).
int[][] frequency(real[] x, real[] y, pair a, pair b, int nx, int ny=nx)
{
  int[][] freq=new int[nx][];
  for(int i=0; i < nx; ++i)
    freq[i]=sequence(new int(int x) {return 0;},ny);
  real hx=nx/(b.x-a.x);
  real hy=ny/(b.y-a.y);
  real ax=a.x;
  real ay=a.y;
  for(int i=0; i < x.length; ++i) {
    int I=Floor((x[i]-ax)*hx);
    int J=Floor((y[i]-ay)*hy);
    if(I >= 0 && I <= nx && J >= 0 && J <= ny)
      ++freq[I][J];
  }
  return freq;
}

int[][] frequency(pair[] z, pair a, pair b, int nx, int ny=nx)
{
  int[][] freq=new int[nx][];
  for(int i=0; i < nx; ++i)
    freq[i]=sequence(new int(int x) {return 0;},ny);
  real hx=nx/(b.x-a.x);
  real hy=ny/(b.y-a.y);
  real ax=a.x;
  real ay=a.y;
  for(int i=0; i < z.length; ++i) {
    int I=Floor((z[i].x-ax)*hx);
    int J=Floor((z[i].y-ay)*hy);
    if(I >= 0 && I < nx && J >= 0 && J < ny)
      ++freq[I][J];
  }
  return freq;
}

path halfbox(pair a, pair b)
{
  return a--(a.x,b.y)--b;
}

path topbox(pair a, pair b)
{
  return a--(a.x,b.y)--b--(b.x,a.y);
}

// Draw a histogram for bin boundaries bin[n+1] of frequency data in count[n].
void histogram(picture pic=currentpicture, real[] bins, real[] count,
               real low=-infinity,
	       pen fillpen=nullpen, pen drawpen=nullpen, bool bars=false,
	       Label legend="", real markersize=legendmarkersize)
{
  if((fillpen == nullpen || bars == true) && drawpen == nullpen)
    drawpen=currentpen;
  bool[] valid=count > 0;
  real m=min(valid ? count : null);
  real M=max(valid ? count : null);
  bounds my=autoscale(pic.scale.y.scale.T(m),pic.scale.y.T(M),
                      pic.scale.y.scale);
  if(low == -infinity) low=pic.scale.y.scale.Tinv(my.min);
  real last=low;
  int n=count.length;
  begingroup(pic);
  for(int i=0; i < n; ++i) {
    if(valid[i]) {
      real c=count[i];
      pair b=Scale(pic,(bins[i+1],c));
      pair a=Scale(pic,(bins[i],low));
      if(fillpen != nullpen) {
	fill(pic,box(a,b),fillpen);
	if(!bars) draw(pic,b--(b.x,a.y),fillpen);
      }
      if(!bars)
	draw(pic,halfbox(Scale(pic,(bins[i],last)),b),drawpen);
      else draw(pic,topbox(a,b),drawpen);
      last=c;
    } else {
      if(!bars && last != low) {
	draw(pic,Scale(pic,(bins[i],last))--Scale(pic,(bins[i],low)),drawpen);
        last=low;
      }
    }
  }
  if(!bars && last != low)
    draw(pic,Scale(pic,(bins[n],last))--Scale(pic,(bins[n],low)),drawpen);
  endgroup(pic);

  if(legend.s != "") {
    marker m=marker(scale(markersize)*shift((-0.5,-0.5))*unitsquare,
		    drawpen,fillpen == nullpen ? Draw :
		    (drawpen == nullpen ? Fill(fillpen) : FillDraw(fillpen)));
    legend.p(drawpen);
    pic.legend.push(Legend(legend.s,legend.p,invisible,m.f));
  }
}

// Draw a histogram for data in n uniform bins between a and b
// (optionally normalized).
void histogram(picture pic=currentpicture, real[] data, real a, real b, int n,
               bool normalize=false, real low=-infinity,
	       pen fillpen=nullpen, pen drawpen=nullpen, bool bars=false,
	       Label legend="", real markersize=legendmarkersize)
{
  real dx=(b-a)/n;
  real[] freq=frequency(data,a,b,n);
  if(normalize) freq /= dx*sum(freq);
  histogram(pic,a+sequence(n+1)*dx,freq,low,fillpen,drawpen,bars,legend,
	    markersize);
}

// Method of Shimazaki and Shinomoto for selecting the optimal number of bins.
// Shimazaki H. and Shinomoto S., A method for selecting the bin size of a
// time histogram, Neural Computation (2007), Vol. 19(6), 1503-1527.
// cf. http://www.ton.scphys.kyoto-u.ac.jp/~hideaki/res/histogram.html
int bins(real[] data, int max=100)
{
  real m=min(data);
  real M=max(data)*(1+epsilon);
  real n=data.length;
  int bins=1;
  real minC=2n-n^2; // Cost function for N=1.
  for(int N=2; N <= max; ++N) {
    real C=N*(2n-sum(frequency(data,m,M,N)^2));
    if(C < minC) {
      minC=C;
      bins=N;
    }
  }

  return bins;
}

// return a pair of central Gaussian random numbers with unit variance
pair Gaussrandpair()
{
  real r2,v1,v2;
  do {
    v1=2.0*unitrand()-1.0;
    v2=2.0*unitrand()-1.0;
    r2=v1*v1+v2*v2;
  } while(r2 >= 1.0 || r2 == 0.0);
  return (v1,v2)*sqrt(-log(r2)/r2);
}

// return a central Gaussian random number with unit variance
real Gaussrand()
{
  static real sqrt2=sqrt(2.0);
  static pair z;
  static bool cached=true;
  cached=!cached;
  if(cached) return sqrt2*z.y;
  z=Gaussrandpair();
  return sqrt2*z.x;
}

struct linefit {
  real m,b;     // slope, intercept
  real dm,db;   // standard error in slope, intercept
  real r;       // correlation coefficient
  real fit(real x) {
    return m*x+b;
  }
}

// Do a least-squares fit of data in real arrays x and y to the line y=m*x+b
linefit leastsquares(real[] x, real[] y)
{
  linefit L;
  int n=x.length;
  if(n == 1) abort("Least squares fit requires at least 2 data points");
  real sx=sum(x);
  real sy=sum(y);
  real sxx=n*sum(x^2)-sx^2;
  real sxy=n*sum(x*y)-sx*sy;
  L.m=sxy/sxx;
  L.b=(sy-L.m*sx)/n;
  if(n > 2) {
    real syy=n*sum(y^2)-sy^2;
    if(sxx == 0 || syy == 0) return L;
    L.r=sxy/sqrt(sxx*syy);
    real arg=syy-sxy^2/sxx;
    if(arg <= 0) return L;
    real s=sqrt(arg/(n-2));
    L.dm=s*sqrt(1/sxx);
    L.db=s*sqrt(1+sx^2/sxx)/n;
  }
  return L;
}

// Do a least-squares fit of data in real arrays x and y weighted by w
// to the line y=m*x+b, by minimizing sum(w*(y-m*x-b)^2).
linefit leastsquares(real[] x, real[] y, real[] w)
{
  linefit L;
  int n=x.length;
  if(n == 1) abort("Least squares fit requires at least 2 data points");
  real sx=sum(w*x);
  real sy=sum(w*y);
  real W=sum(w);
  real sxx=W*sum(w*x^2)-sx^2;
  real sxy=W*sum(w*x*y)-sx*sy;
  L.m=sxy/sxx;
  L.b=(sy-L.m*sx)/W;
  if(n > 2) {
    real syy=W*sum(w*y^2)-sy^2;
    if(sxx == 0 || syy == 0) return L;
    L.r=sxy/sqrt(sxx*syy);
    real arg=syy-sxy^2/sxx;
    if(arg <= 0) return L;
    real s=sqrt(arg/(n-2));
    L.dm=s*sqrt(1/sxx);
    L.db=s*sqrt(1+sx^2/sxx)/W;
  }
  return L;
}