# Use a better version of the octave interface
from octave import octave
%default_mode octave

a = [1 2 3]

randn(10)

randn(10)

x + 1 = 2

x = 2 - 1

#Find a good increment to use to explore points near tails
endval=0.999;
incval=2*endval/20;
xuniform=[-endval:incval:endval]'; %For the simple test
yuniform=[-endval:incval:endval]'; %For the simple test

#For a Gaussian 2D density with zero covariance the density is factorable permitting...
#Use inverse error function to generate a Gaussian marginal distribution of trajectories

xvalues=erfinv(xuniform);
yvalues=erfinv(yuniform);
nxv=size(xvalues,1);
nyv=size(yvalues,1);
N=nxv*nyv; %Number of worlds

#We need to build the 2D coordinates obtained from Cartesian product of x and y coordinates
[Xc, Yc] = meshgrid (xvalues, yvalues);

#X is the state vector, a single column with every x and y coordinate as an entry followed
#by every time derivative of these coordinates 
#For the 2D case our stacking will be of all x coordinates first, then all y, then derivatives
Xcol=reshape(Xc,N,1);
Ycol=reshape(Yc,N,1);

#The function voronoin will generate two lists:
#C contains the vertices of the Voronoi polygons. 
#F contains, for each original point, the indices of the Voronoi polygon vertices. 
#while polyarea will create the areas of each voronoi polygon for us
[C, F] = voronoin ([Xcol, Ycol]);
af = zeros (size (F));
for i = 1 : length (F)
  af(i) = polyarea (C (F {i, :}, 1), C (F {i, :}, 2));
endfor
#display(af) #display areas

#So the probability density associated with each trajectory point is
Pin=(N*af).^(-1);

%Plot the input trajectory density
#Next 2 lines only needed to avoid plot data error
Pf=Pin;
Pf(isnan(Pin))=0;

%sage
m=matrix(RR,octave("[Xcol(:), Ycol(:), Pf(:)]"))
data=[(x,y,z) for x,y,z in m.rows()]
surface=list_plot3d(data,aspect_ratio=[1,1,8],texture={'color':'green','opacity':0.85,'ambient':0.9})
pts=points(data,size=2,color='red')
show(surface+pts,spin=10,background='pink',height=int(600),width=int(1000),renderer='webgl',frame={'thickness':4.0, 'color ':'red','labels' : True,'fontface' : 'roman'} )