function [bandwidth,density,X,Y]=kde2d(data,n,MIN_XY,MAX_XY) 
global N A2 I 
if nargin<2 
    n=2^8; 
end 
n=2^ceil(log2(n)); 
N=size(data,1); 
if nargin<3 
    MAX=max(data,[],1); MIN=min(data,[],1); Range=MAX-MIN; 
    MAX_XY=MAX+Range/4; MIN_XY=MIN-Range/4; 
end 
scaling=MAX_XY-MIN_XY; 
if N<=size(data,2) 
    error('data has to be an N by 2 array where each row represents a two dimensional observation') 
end 
transformed_data=(data-repmat(MIN_XY,N,1))./repmat(scaling,N,1); 

initial_data=ndhist(transformed_data,n); 

a= dct2d(initial_data); 

  I=(0:n-1).^2; A2=a.^2; 
 
 t_star=fzero( @(t)(t-evolve(t)),[0,0.1]); 
 
p_02=func([0,2],t_star);p_20=func([2,0],t_star); p_11=func([1,1],t_star); 
t_y=(p_02^(3/4)/(4*pi*N*p_20^(3/4)*(p_11+sqrt(p_20*p_02))))^(1/3); 
t_x=(p_20^(3/4)/(4*pi*N*p_02^(3/4)*(p_11+sqrt(p_20*p_02))))^(1/3); 

a_t=exp(-(0:n-1)'.^2*pi^2*t_x/2)*exp(-(0:n-1).^2*pi^2*t_y/2).*a;  

if nargout>1 
    density=idct2d(a_t)*(numel(a_t)/prod(scaling)); 
    [X,Y]=meshgrid(MIN_XY(1):scaling(1)/(n-1):MAX_XY(1),MIN_XY(2):scaling(2)/(n-1):MAX_XY(2)); 
end 
bandwidth=sqrt([t_x,t_y]).*scaling;  
end 
%####################################### 
function  [out,time]=evolve(t) 
global N 
Sum_func = func([0,2],t) + func([2,0],t) + 2*func([1,1],t); 
time=(2*pi*N*Sum_func)^(-1/3); 
out=(t-time)/time; 
end 
%####################################### 
function out=func(s,t) 
global N 
if sum(s)<=4 
    Sum_func=func([s(1)+1,s(2)],t)+func([s(1),s(2)+1],t); const=(1+1/2^(sum(s)+1))/3; 
    time=(-2*const*K(s(1))*K(s(2))/N/Sum_func)^(1/(2+sum(s))); 
    out=psi(s,time); 
else 
    out=psi(s,t); 
end 
 
end 
%####################################### 
function out=psi(s,Time) 
global I A2 
% s is a vector 
w=exp(-I*pi^2*Time).*[1,.5*ones(1,length(I)-1)]; 
wx=w.*(I.^s(1)); 
wy=w.*(I.^s(2)); 
out=(-1)^sum(s)*(wy*A2*wx')*pi^(2*sum(s)); 
end 
%####################################### 
function out=K(s) 
out=(-1)^s*prod((1:2:2*s-1))/sqrt(2*pi); 
end 
%####################################### 
function data=dct2d(data) 
% computes the 2 dimensional discrete cosine transform of data 
% data is an nd cube 
[nrows,ncols]= size(data); 
if nrows~=ncols 
    error('data is not a square array!') 
end 
% Compute weights to multiply DFT coefficients 
w = [1;2*(exp(-i*(1:nrows-1)*pi/(2*nrows))).']; 
weight=w(:,ones(1,ncols)); 
data=dct1d(dct1d(data)')'; 
    function transform1d=dct1d(x) 
 
        % Re-order the elements of the columns of x 
        x = [ x(1:2:end,:); x(end:-2:2,:) ]; 
 
        % Multiply FFT by weights: 
        transform1d = real(weight.* fft(x)); 
    end 
end 
%####################################### 
function data = idct2d(data) 
% computes the 2 dimensional inverse discrete cosine transform 
[nrows,ncols]=size(data); 
% Compute wieghts 
w = exp(i*(0:nrows-1)*pi/(2*nrows)).'; 
weights=w(:,ones(1,ncols)); 
data=idct1d(idct1d(data)'); 
    function out=idct1d(x) 
        y = real(ifft(weights.*x)); 
        out = zeros(nrows,ncols); 
        out(1:2:nrows,:) = y(1:nrows/2,:); 
        out(2:2:nrows,:) = y(nrows:-1:nrows/2+1,:); 
    end 
end 
%####################################### 
function binned_data=ndhist(data,M) 
% this function computes the histogram 
% of an n-dimensional data set; 
% 'data' is nrows by n columns 
% M is the number of bins used in each dimension 
% so that 'binned_data' is a hypercube with 
% size length equal to M; 
[nrows,ncols]=size(data); 
bins=zeros(nrows,ncols); 
for i=1:ncols 
    [dum,bins(:,i)] = histc(data(:,i),[0:1/M:1],1); 
    bins(:,i) = min(bins(:,i),M); 
end 
% Combine the  vectors of 1D bin counts into a grid of nD bin 
% counts. 
binned_data = accumarray(bins(all(bins>0,2),:),1/nrows,M(ones(1,ncols))); 
end 