# code for investigating samples and populations
# Anthony Smith 2014

import numpy
import matplotlib.pyplot as plt

inputtype=1

if inputtype==1:
	samplefraction=numpy.double(raw_input("Total Fraction of the Population Sampled: "))
	
	finaldensity=numpy.double(raw_input("Density of the Population: "))
	L=numpy.double(raw_input("Structure Parameter: "))

if inputtype==2:
	samplefraction=0.5
	
	finaldensity=1000
	L=0 # L=0 for an unstructured random set

density=0

scalefact=2

areaclusteredxs=[]
areaclusteredys=[]

while density<finaldensity:
	l=0;
	while l<=L:
		if l==0:
			seedxs=numpy.random.rand(1,1)
			seedys=numpy.random.rand(1,1)
			
			clusteredxs=seedxs;
			clusteredys=seedys;
			
			l=1;
			originalseparation=((numpy.random.rand(1))*0.5+0.5);
		numseeds=seedxs.shape[0]
		thetas=(numpy.random.rand(numseeds,1))*2*3.141592653589793;
		separation=originalseparation/(numpy.power(scalefact,(l-1)));
		
		halfseparation=0.5*separation;
		
		newseedxdiffs=halfseparation*(numpy.cos(thetas));
		newseedydiffs=halfseparation*(numpy.sin(thetas));
		positivenewseedxs=seedxs+newseedxdiffs #*numpy.ones((seedxs.shape[0],seedxs.shape[1]));
		positivenewseedys=seedys+newseedydiffs;
		negativenewseedxs=seedxs-newseedxdiffs;
		negativenewseedys=seedys-newseedydiffs;
		
		seedxs=numpy.concatenate((positivenewseedxs,negativenewseedxs),axis=0)
		seedys=numpy.concatenate((positivenewseedys,negativenewseedys),axis=0)
		clusteredxs=numpy.concatenate((clusteredxs,numpy.reshape(positivenewseedxs,(positivenewseedxs.shape[0]*positivenewseedxs.shape[1],1)),numpy.reshape((negativenewseedxs),(negativenewseedxs.shape[0]*negativenewseedxs.shape[1],1))),axis=0)
		clusteredys=numpy.concatenate((clusteredys,numpy.reshape(positivenewseedys,(positivenewseedys.shape[0]*positivenewseedys.shape[1],1)),numpy.reshape((negativenewseedys),(negativenewseedys.shape[0]*negativenewseedys.shape[1],1))),axis=0)
		
		l=l+1;
	areaclusteredxs=numpy.concatenate((areaclusteredxs,clusteredxs[numpy.nonzero((0<clusteredxs)*(clusteredxs<1)*(0<clusteredys)*(clusteredys<1)*(clusteredxs))]),axis=0)
	areaclusteredys=numpy.concatenate((areaclusteredys,clusteredys[numpy.nonzero((0<clusteredxs)*(clusteredxs<1)*(0<clusteredys)*(clusteredys<1)*(clusteredys))]),axis=0)
	density=areaclusteredxs.shape[0]

onesamplearea=samplefraction/4
onesampleside=numpy.power(onesamplearea,0.5)

#sample1count=sum(((0)<areaclusteredxs)*(areaclusteredxs<(onesampleside))*((0)<areaclusteredys)*(areaclusteredys<(onesampleside))) # lower left
#sample2count=sum(((0)<areaclusteredxs)*(areaclusteredxs<(onesampleside))*((1-onesampleside)<areaclusteredys)*(areaclusteredys<(1))) # upper left
#sample3count=sum(((1-onesampleside)<areaclusteredxs)*(areaclusteredxs<(1))*((0)<areaclusteredys)*(areaclusteredys<(onesampleside))) # lower right
#sample4count=sum(((1-onesampleside)<areaclusteredxs)*(areaclusteredxs<(1))*((1-onesampleside)<areaclusteredys)*(areaclusteredys<(1))) # upper right

sample1count=sum(((0.25-0.5*onesampleside)<areaclusteredxs)*(areaclusteredxs<(0.25+0.5*onesampleside))*((0.25-0.5*onesampleside)<areaclusteredys)*(areaclusteredys<(0.25+0.5*onesampleside))) # lower left
sample2count=sum(((0.25-0.5*onesampleside)<areaclusteredxs)*(areaclusteredxs<(0.25+0.5*onesampleside))*((0.75-0.5*onesampleside)<areaclusteredys)*(areaclusteredys<(0.75+0.5*onesampleside))) # upper left
sample3count=sum(((0.75-0.5*onesampleside)<areaclusteredxs)*(areaclusteredxs<(0.75+0.5*onesampleside))*((0.25-0.5*onesampleside)<areaclusteredys)*(areaclusteredys<(0.25+0.5*onesampleside))) # lower right
sample4count=sum(((0.75-0.5*onesampleside)<areaclusteredxs)*(areaclusteredxs<(0.75+0.5*onesampleside))*((0.75-0.5*onesampleside)<areaclusteredys)*(areaclusteredys<(0.75+0.5*onesampleside))) # upper right

sample1density=sample1count/onesamplearea
sample2density=sample2count/onesamplearea
sample3density=sample3count/onesamplearea
sample4density=sample4count/onesamplearea

meandensities=numpy.mean([sample1density,sample2density,sample3density,sample4density])
stddensities=numpy.std([sample1density,sample2density,sample3density,sample4density])

plt.figure(num=None,figsize=(10,10))

plt.subplot(2,1,1)

plt.hold(True)

#numpy.random.shuffle(areaclusteredxs)
#numpy.random.shuffle(areaclusteredys)

plt.plot(areaclusteredxs,areaclusteredys,'k.')

#plt.plot(numpy.array([0,onesampleside,onesampleside]),numpy.array([onesampleside,onesampleside,0]),'r-') # lower left
#plt.plot(numpy.array([0,onesampleside,onesampleside]),numpy.array([1-onesampleside,1-onesampleside,1]),'r-') # upper left
#plt.plot(numpy.array([1-onesampleside,1-onesampleside,1]),numpy.array([0,onesampleside,onesampleside]),'r-') # lower right
#plt.plot(numpy.array([1-onesampleside,1-onesampleside,1]),numpy.array([1,1-onesampleside,1-onesampleside]),'r-') # upper right

plt.plot(numpy.array([0.25-0.5*onesampleside,0.25+0.5*onesampleside,0.25+0.5*onesampleside,0.25-0.5*onesampleside,0.25-0.5*onesampleside]),numpy.array([0.25+0.5*onesampleside,0.25+0.5*onesampleside,0.25-0.5*onesampleside,0.25-0.5*onesampleside,0.25+0.5*onesampleside]),'r-') # lower left
plt.plot(numpy.array([0.25-0.5*onesampleside,0.25+0.5*onesampleside,0.25+0.5*onesampleside,0.25-0.5*onesampleside,0.25-0.5*onesampleside]),numpy.array([0.75-0.5*onesampleside,0.75-0.5*onesampleside,0.75+0.5*onesampleside,0.75+0.5*onesampleside,0.75-0.5*onesampleside]),'r-') # upper left
plt.plot(numpy.array([0.75-0.5*onesampleside,0.75-0.5*onesampleside,0.75+0.5*onesampleside,0.75+0.5*onesampleside,0.75-0.5*onesampleside]),numpy.array([0.25-0.5*onesampleside,0.25+0.5*onesampleside,0.25+0.5*onesampleside,0.25-0.5*onesampleside,0.25-0.5*onesampleside]),'r-') # lower right
plt.plot(numpy.array([0.75-0.5*onesampleside,0.75-0.5*onesampleside,0.75+0.5*onesampleside,0.75+0.5*onesampleside,0.75-0.5*onesampleside]),numpy.array([0.75+0.5*onesampleside,0.75-0.5*onesampleside,0.75-0.5*onesampleside,0.75+0.5*onesampleside,0.75+0.5*onesampleside]),'r-') # upper right

plt.xticks([])
plt.yticks([])

plt.subplot(2,1,2)

plt.plot(numpy.array([0.5,0.5,0.5,0.5]),numpy.array([sample1density,sample2density,sample3density,sample4density]),'k.',label='Individual Samples')
plt.plot(0.5,density,'b.',label='Population')
plt.plot(numpy.array([0.4,0.6]),numpy.array([meandensities+stddensities,meandensities+stddensities]),'r-',label='Standard Deviation from Samples')
plt.plot(numpy.array([0.4,0.6]),numpy.array([meandensities-stddensities,meandensities-stddensities]),'r-')

plt.ylabel('Density',size=20)

plt.xlim([0,2])
plt.xticks([])

plt.legend()

plt.show()