# Statistics. Identification of 2D shapes  in data
# Performs calculations of eccentricities 
# authors:  C.Levy and S.Chekanov. Phys.Rev.D82:094029,2010
# based on https://atlaswww.hep.anl.gov/asc/statshape/

from jhplot  import *
from java.util import Random
from java.awt import Color, Font
from jhplot.stat  import StatShape

c1 = HPlot("Canvas",500,450)
fo1=Font("Arial", Font.BOLD, 24)
c1.setMarginLeft(90)
c1.setNameX("η", fo1)
c1.setNameY("φ", fo1)
c1.setRange(-5.5,5.5,-5.5,5.5)
c1.visible()
c1.setLegend(0)

g1=HLabel( "(1)",0.8,0.25,"NDC")
g1.setFont(fo1)
c1.add(g1)

r=Random(22);
x,y,w=[],[],[]
for i in range(200):
          x.append(r.nextGaussian())
          y.append(r.nextGaussian())
          w.append(1)
                    
for i in range(200):
          x.append(r.nextGaussian()) 
          y.append(r.nextGaussian())
          w.append(1)
 

ss=StatShape(x,y,w,len(w))

# use unweighted linear regression and weighted means in quadrants
ss.process(1)
data=ss.getData()
data.setSymbolSize(4)
c1.draw(data)

sum=ss.getSummary()
# draw major axis
f=ss.getFitFunction()
f.setPenWidth(4)
f.setColor(Color.red)
c1.draw(f)
means = ss.getMeans()
print "weighted meanX=",means[0]
print "weighted meanY=",means[1]

# draw minor axis
ff=ss.getFitFunctionPerp() 
ff.setPenWidth(2)
c1.draw(ff)
 
print "Major=" ,sum[0]
print "Minor=", sum[1]
# eccentricity is define as 1- Minor/Major
print "Eccentricity=", sum[2]
print "MajorEccentricity=",sum[20]
print "MinorEccentricity=", sum[21]
print "MajorLength=",sum[14], " Global Major Length=",sum[9]
print "MinorLength=", sum[15], " Global Minor Length=", sum[10]
print "Fmax (fraction of weight of a point with max weight)=", sum[22]

k1 = ss.getCenters(0)
k2 = ss.getCenters(1)
k3 = ss.getCenters(2)
k4 = ss.getCenters(3)


type=0
p1=P1D("mean1",Color.red)
p1.setSymbol(type)
p1.add(k1[0],k1[1])

p2=P1D("mean2",Color.red)
p2.setSymbol(type)
p2.add(k2[0],k2[1])

p3=P1D("mean3",Color.red)
p3.setSymbol(type)
p3.add(k3[0],k3[1])

p4=P1D("mean4",Color.red)
p4.setSymbol(type)
p4.add(k4[0],k4[1])

size=16
p1.setSymbolSize(size)
p2.setSymbolSize(size)
p3.setSymbolSize(size)
p4.setSymbolSize(size)
c1.draw(p1)
c1.draw(p2)
c1.draw(p3)
c1.draw(p4)

s1="ECC= %.2f" % sum[16]
s2="ECC_{MJ}= %.2f" % sum[20]
s3="ECC_{MI}= %.2f" % sum[21]
xpos=0.19
font=Font("Arial", Font.BOLD, 14)
lab0=HLabel( "NQM:",xpos,0.85,"NDC")
lab1=HLabel( s1,xpos,0.81,"NDC")
lab2=HLabel( s2,xpos,0.77,"NDC")
lab3=HLabel( s3,xpos,0.73,"NDC")
lab0.setFont(Font("Arial", Font.BOLD, 18))
lab1.setFont(font)
lab2.setFont(font)
lab3.setFont(font)
c1.add(lab0)
c1.add(lab1)
c1.add(lab2)
c1.add(lab3)

s1="ECC= %.2f" % sum[2]
s2="ECC_{MJ} = %.2f" % sum[5]
s3="ECC_{MI}= %.2f" % sum[8]
xpos=0.69
font=Font("Arial", Font.BOLD, 14)
lab0=HLabel( "QM:",xpos,0.85,"NDC")
lab1=HLabel( s1,xpos,0.81,"NDC")
lab2=HLabel( s2,xpos,0.77,"NDC")
lab3=HLabel( s3,xpos,0.73,"NDC")
lab0.setFont(Font("Arial", Font.BOLD, 18))
lab1.setFont(font)
lab2.setFont(font)
lab3.setFont(font)
c1.add(lab0)
c1.add(lab1)
c1.add(lab2)
c1.add(lab3)

c1.update()
c1.export("stat1.eps")