import monteCarloPCP as mc
import numpy as np
import scipy.linalg as la
import sys
import os
import subprocess
import time
import glob
import pdb


gold = 0.61803398875

def pbcDist(pos1,pos2):
	delta = pos1 - pos2
	delta -= np.round(delta)
	distance = np.sqrt(np.sum(delta**2))
	return distance

def recalcNeighbors(neighborCutDistance,lastPositions):
	# if (p.getSumOfMaxDisplacements()/np.min(p.getRadii()) >= neighborCutDistance+1):
	if (p.getSumOfMaxDisplacements() >= 2*neighborCutDistance*np.min(p.getRadii())):
		p.setPositions(lastPositions)
		p.calcNeighborsCut(neighborCutDistance)
		p.resetDisplacement()

def getGaps(neighborCut = 1.0):
	# print("blah 1")
	p.setPotentialType(mc.potentialEnum.log)
	# print("blah 2")
	# print("blah 3")
	# print(p.getPositions())
	p.calcNeighborsCut(neighborCut)
	# print("blah 3")
	gaps = p.getSmallGaps(neighborCut)
	# while len(gaps.data) == 0:
	#       print("no gaps, recalculating, neighborCut=" +str(neighborCut))
	#       neighborCut = 10*neighborCut
	#       p.calcNeighborsCut(neighborCut)
	#       gaps = p.getSmallGaps(neighborCut)
	# print("blah 4 " + str(neighborCut) + " " + str(len(p.getNeighbors().data)))
	minGap = p.getMinGap(10)
	# print("blah 5")
	return minGap, gaps.data

def findJammingGap(precision = 1e-2,neighborCut = 1.0, hi = 1.0, lo = 0.0):
	p.setPotentialType(mc.potentialEnum.contactPower)
	recalcNeighbors(neighborCut,p.getPositions())
	p.calcForceEnergy()
	resetPhi = p.getPhi()
	check = (hi-lo)*gold
	numRattlers = 0
	while (hi-lo)/hi > precision:
		check = (hi-lo)*gold + lo
		p.setPhi(resetPhi + check)
		# recalcNeighbors(neighborCut)
		# p.calcForceEnergy()
		stableAndExcess = p.getStableListAndExcess()
		if stableAndExcess[1] is np.nan:
			lo = check
		else:
			numRattlers = nParticles - np.sum(stableAndExcess[0])
			hi = check
		print("hi is " + str(hi) + " , lo is " + str(lo) + " , and check is " + str(check))
	p.setPotentialType(mc.potentialEnum.log)
	print("estimated phi error is  " + str(check) + " and estimated phi is " + str(check + resetPhi))
	print("Rattler fraction is " + str(numRattlers/nParticles))
	p.setPhi(resetPhi)
	recalcNeighbors(neighborCut,p.getPositions())
	p.calcForceEnergy()
	return hi, lo, check

def findJammingGapSoft(precision = 1e-2,neighborCut = 1.0, hi = 1.0, lo = 0.0):
	p.setPotentialType(mc.potentialEnum.contactPower)
	p.calcForceEnergy()
	resetPhi = p.getPhi()
	check = (hi-lo)*gold
	numRattlers = 0
	beforePos = p.getPositions()
	preNeigh = p.getNeighbors()
	while (hi-lo)/hi > precision:
		check = (hi-lo)*gold + lo
		p.setPositions(beforePos)
		p.setPhi(resetPhi + check)
		recalcNeighbors(neighborCut,beforePos)
		p.calcForceEnergy()
		p.minimizeFIRE(criticalForce = 1e-15)
		stableAndExcess = p.getStableListAndExcess()
		if stableAndExcess[1] >= 0:
			numRattlers = nParticles - np.sum(stableAndExcess[0])
			hi = check
		else:
			lo = check
		print("hi is " + str(hi) + " , lo is " + str(lo) + " , and check is " + str(check))
	# p.setPotentialType(mc.potentialEnum.log)
	print("estimated phi error is  " + str(check) + " and estimated phi is " + str(check + resetPhi))
	print("Rattler fraction is " + str(numRattlers/nParticles))
	p.setPhi(resetPhi)
	p.setPositions(beforePos)
	recalcNeighbors(neighborCut,beforePos)
	p.calcForceEnergy()
	return hi, lo, check

def findNewExpansion(scaleFactor):
	minGap, normGaps = getGaps(p.getLogZero()*5.0)
	print("minGap is " + str(minGap))
	# print("maxR is " + str(maxR))
	#expansionFactor = scaleFactor*dist/radSum + np.quad('1') - scaleFactor
	#expansionFactor = np.quad('1')/((np.quad('1')-scaleFactor) + scaleFactor*radSum/dist)
	expansionFactor = 1.0/((1.0-scaleFactor) + scaleFactor/(1.0 + minGap))
	if (expansionFactor)**p.getNDim() > 1:
		return (expansionFactor)**p.getNDim(), minGap, normGaps
	else:
		return 1, minGap, normGaps

def setNewPhi(scaleFactor = 0.9, newPhi = None, gap = 1e-5,neighborCut=1.0):
	p.setPotentialType(mc.potentialEnum.contactPower)
	recalcNeighbors(neighborCut,p.getPositions())
	oldPhi = p.getPhi()
	oldPositions = p.getPositions()
	# print("Old phi is " + str(p.getPhi()))
	if newPhi is None:
		expansionFactor, minGap, normGaps = findNewExpansion(scaleFactor)
		p.setPhi(expansionFactor*oldPhi)
	else:
		p.setPhi(newPhi)
	recalcNeighbors(neighborCut,oldPositions)
	print("New phi is " + str(p.getPhi()))
	p.setPotentialType(mc.potentialEnum.log)
	return minGap, normGaps

def findLogCutoff(zeroFactor,neighborCut=5.0):
	# print("Shit we here 1 ")
	minGap, normGaps = getGaps(neighborCut = neighborCut)
	# print("Shit we here 2")
	sortedGaps = np.sort(normGaps)
	# print("Shit we here 3")
	return sortedGaps[np.min([int(zeroFactor*nParticles*dim + 1),len(sortedGaps)-1])], minGap, normGaps, sortedGaps

def simpleLogMin(oldPhi, criticalForce=1e-15,dt_max=0.1,dt=0.001,dtMaxFactor = 10.0,neighborCut=1.0,scaleFactor = 0.5,maxIterations = 100000,zeroFactor = 2.5,saveBeforeExpansion = False, phiJError = None,gardnerFound = False):
	p.setPotentialType(mc.potentialEnum.log)
	# print("here 1")
	logZero, minGap, normGaps, sortedGaps = findLogCutoff(zeroFactor,neighborCut)
	neighborCut = logZero*5.0
	# print("here 2")
	dtMax = np.sqrt(minGap)
	# logZero = np.min([logZero,0.1])
	p.setLogZero(logZero)
	recalcNeighbors(neighborCut,p.getPositions())
	# print("here 3")
	p.calcForceEnergy()
	if saveBeforeExpansion:
		p.save(saveDir + '/beforeMinimization/' + str(p.getPhi()))
	dtMax1 = dtMax
	# print("here 4")
	runningTotal = 0
	fireIterations = 0
	thisLastPositions = p.getPositions()
	while fireIterations < 1:
		#Only pass in the gaps that are used
		# print("here 5")
		fireIterations, minGap1 = minimizeFIRESave(lastPositions = thisLastPositions, oldPhi = oldPhi, criticalForce=criticalForce,neighborCut = neighborCut,dt_max = dtMax1,dt = dtMax1/50.0,recalcNeighbors=True,scaleFreq=1600, maxIterations = maxIterations,saveDuringFIRE = saveBeforeExpansion,scaleFactor=scaleFactor, sortedGaps = sortedGaps[sortedGaps < logZero],gardnerFound = gardnerFound)
		dtMax1 /= 1.1
		runningTotal += fireIterations
		if fireIterations < 1:
			neighborCut *= 1.1
			p.calcNeighborsCut(neighborCut)
			p.calcForceEnergy()
		# print("here 6")
			# print("minGap is " + str(minGap) + " and minGap1 is " + str(minGap1))
		# if dtMax1 > np.quad('1e-3') or minGap1 > np.quad('10')*minGap:
		# else:
		#       print("phiJError is " + str(phiJError))
		#       _,phiJError,_ = findJammingGap(hi = np.quad('2')*phiJError, lo = phiJError/np.quad('10'))
	if saveBeforeExpansion:
		p.save(saveDir + '/beforeExpansion/' + str(p.getPhi()))
	minGap, normGaps = setNewPhi(scaleFactor = scaleFactor,neighborCut=neighborCut)
	if saveBeforeExpansion:
		p.save(saveDir + '/afterExpansion/' + str(p.getPhi()))
	return dtMax1, minGap1, normGaps, runningTotal

def minimizeFIRESave(lastPositions, sortedGaps, oldPhi, saveDuringFIRE = False, neighborCut = 1.0, scaleFactor=0.9, dt_max = 0.1, dt = 0.001, scaleFreq = 16,maxIterations = -1,gardnerFound = False,*args, **kwargs):
	"""
	Runs a single FIRE minimization. See FIREMinimizer for further
	documentation
	"""
	iteration = 0
	minGap = p.getMinGap(100)
	gapCutoff = (1.0/(1.0 - scaleFactor))*minGap
	# gapCutoff = np.min([(1.0/(1.0 - scaleFactor))*minGap,0.1])
	lastMinGap = minGap
	# p.save(saveDir + '/' + str(p.getPhi()))
	# mc.fileIO.saveScalar(saveDir + '/' + str(p.getPhi()),'logZero',p.getLogZero())
	if p.getSmallOverlaps().size > 0:
		p.setPositions(lastPositions)
		p.calcNeighborsCut(neighborCut)
		p.resetDisplacement()
	lastPos = p.getPositions()
	totalIterations = 0
	# print("minGap " + str(minGap) + ", gap cutoff " + str(gapCutoff) + " " + str(p.getMaxUnbalancedForce()))
	startingMinGap = minGap
	if gardnerFound:
		gapCutoff = 1000000*minGap
	neighborCut = 5.0*p.getLogZero()
	while ((minGap < gapCutoff and totalIterations < maxIterations)):
		# print("we here")
		# print("max un force = " + str(p.getMaxUnbalancedForce()))
		# print("minGap " + str(minGap) + ", gap cutoff " + str(gapCutoff) + " " + str(p.getMaxUnbalancedForce()))
		# print("total it " + str(totalIterations))
		minGap = p.getMinGap(100)
		# print("done gotem")
		dt_max = np.sqrt(minGap)
		dt = np.min([dt,dt_max/10.0])
		minForce = 0.1
		# print("now here")
		if minGap > 0.5:
			dt = np.sqrt(minGap)/10
			dt_max = np.sqrt(minGap)/2
		for iteration in p.FIREMinimizer(dt_max = dt_max, dt = dt,maxIterations = maxIterations-totalIterations,scaleFreq=1600, cutDistance = neighborCut,minCut = neighborCut,initialCut=neighborCut,*args, **kwargs):
			smallOverlaps = p.getSmallOverlaps().data
			# print("max force: " + str(p.getMaxUnbalancedForce()))
			# print("log zero: " + str(p.getLogZero()))
			# print(iteration)
			# print(minForce*dt**2)
			# print(minGap)
			if p.getMaxUnbalancedForce() == 0.0 or len(smallOverlaps):
				if np.any(np.isnan(p.getPositions())) or len(smallOverlaps):
					# print("max force: " + str(p.getMaxUnbalancedForce()))
					# print("log zero: " + str(p.getLogZero()))
					p.setPositions(lastPos)
					# print("No here")
					# print(lastPos)
					iteration -= 1
					# print(minForce)
				p.zeroVelocities()
				minForce = minForce/2
				# print("done gotta do this part 1")
				neighborCut = 1.1*neighborCut
				# print("done gotta do this part2")
				p.calcNeighborsCut(neighborCut)
				# print("done gotta do this part3")
				# recalcNeighbors(neighborCut)
				p.calcForceEnergy()
				p.setMaxForceRescale(minForce)
				p.rescaleForces(minForce)
				if p.getMaxUnbalancedForce() == 0.0:
					logZero,minGap,_,_ = findLogCutoff(zeroFactor,neighborCut)
					# print("done gotta do this part4")
					# print(logZero)
					# print("done gotta do this part5")
					# neighborCut = 2*logZero
					# print("done gotta do this part6")
					# logZero = np.min([logZero,0.1])
					p.setLogZero(logZero)
					p.calcForceEnergy()
				# p.rescaleForces(0.1)
				# print("done gotta do this part8")
				minGap = p.getMinGap(minGap*10)
			else:
				minForce = 0.1
				p.setMaxForceRescale(minForce)
				# print("Here")
				# print("y aqui " + str(iteration))
				# print("gap is " + str(len(p.getSmallGaps().data)))
				# print("y aqui 2 " + str(iteration))
				minGap = p.getMinGap(minGap*10)
				# print("y aqui 1 " + str(iteration))
				dt_max = np.sqrt(minGap)
				if minGap > 0.5:
					dt = np.sqrt(minGap)/10
					dt_max = np.sqrt(minGap)/2
					neighborCut = 10.0*minGap
					p.calcNeighborsCut(neighborCut)
					# print("recalc")
				if neighborCut < 3.0*minGap:
					neighborCut = 10.0*minGap
					p.calcNeighborsCut(neighborCut)
				# neighborCut = np.max([neighborCut,3.0*minGap])
				dt = np.min([dt,dt_max])
				# print("dt is " + str(p.getFIRETimeStep()))
				# if p.getSmallOverlaps().size > 0:
				# 	p.setPositions(lastPos)
				# 	recalcNeighbors(neighborCut)
				# 	p.calcForceEnergy()
				# 	p.rescaleForces(0.1)
				# 	print("Had to go back")
				# 	p.setFIRETimeStep(p.getFIRETimeStep()/100.0)
				# 	print("This was iteration " + str(iteration))
				# 	minGap = p.getMinGap(100)
				# 	iteration = iteration - 1
				# 	# print("breaking here")
				# 	break
				# 	return iteration - 1, lastMinGap
				# else:
				lastPos = p.getPositions()
				lastMinGap = minGap
				if saveDuringFIRE:
					p.save(saveDir + '/duringFIRE/' + str(p.getPhi()) + '/' + str(iteration) + '/')
				if minGap > gapCutoff:
					# print("breaking here 1, totalIterations is  " + str(totalIterations) + " , iterations is " + str(iteration))
					break
				# print("iterations is " + str(iteration))
				# print("totalIterations is " + str(totalIterations))
				# print("minGap is " + str(minGap))
				# print("gapCutoff is " + str(gapCutoff))
				# print("max unbalanced force is " + str(p.getMaxUnbalancedForce()))
				# print("num overlaps " + str(len(p.getSmallOverlaps().data)))
				# print("num neighbors is " + str(len(p.getNeighbors().data)))
				# print("neighbor cutoff is " + str(neighborCut))
			pass
		totalIterations += iteration
		# print("now we up in this b")
		# if p.getMaxUnbalancedForce() == 0:
		#       neighborCut = 10.0*minGap
		#       p.calcNeighborsCut(neighborCut)
		#       logZero,minGap,_,_ = findLogCutoff(zeroFactor,neighborCut)
		#       p.setLogZero(zeroFactor)
		# print("iterations is " + str(iteration))
		# print("totalIterations is " + str(totalIterations))
		# print("minGap is " + str(minGap))
		# print("gapCutoff is " + str(gapCutoff))
		# print("max unbalanced force is " + str(p.getMaxUnbalancedForce()))
		# print("num overlaps " + str(len(p.getSmallOverlaps().data)))
		# print("num neighbors is " + str(len(p.getNeighbors().data)))
	return totalIterations, minGap

def isoSearch(beforePos, phiInitial, phiC, criticalForce = 1e-10):
	p.setPositions(beforePos)
	p.setPhi(phiInitial)
	recalcNeighbors(1.0,p.getPositions())
	p.calcForceEnergy()
	for iteration in p.isostaticitySearch(phiC = phiC, phiInitial = phiInitial, criticalForce = criticalForce):
		pass

def file_len(fname):
	with open(fname) as f:
		for i, l in enumerate(f):
			pass
	return i + 1

# seed = 35
# dim = 3
# nParticles = 2048
# phiString = '0.54'
# maxIterationFactorString = str(2*dim)
# scaleFactorString = '0.9'
# zeroFactorString = '2'
# biDisperse = 1
# saveMiddleSteps = 0
# # basedir = '/work/pkm8/Data/logMCMinGapReworkGardnerCopies'
# basedir = '/work/pkm8/Data/logMCMinGapRework'
# basedir = '/pylon5/mr4s81p/pkmorse/Data/logMCMinGapRework'

seed = int(sys.argv[1])
dim = int(sys.argv[2])
nParticles = int(sys.argv[3])
phiString = sys.argv[4]
maxIterationFactorString = sys.argv[5]
scaleFactorString = sys.argv[6]
zeroFactorString = sys.argv[7]
biDisperse = int(sys.argv[8])
saveMiddleSteps = int(sys.argv[9])
basedir = sys.argv[10]

start = time.time()


maxIterationFactor = np.float64(maxIterationFactorString)
maxIterations = np.floor(maxIterationFactor * nParticles)
criticalForce = 1e-10
zeroFactor = np.float64(zeroFactorString)
scaleFactor = np.float64(scaleFactorString)
expCutoff = 0.367879441171442

saveDir = basedir + '/' + str(dim) + '/' + str(nParticles) + '/' + phiString + '/' + str(seed) + '/' + zeroFactorString + '/' + maxIterationFactorString + '/' + scaleFactorString + '/'
startingFolder = basedir + '/' + str(dim) + '/' + str(nParticles) + '/' + phiString + '/' + str(seed) + '/beforeMin/'
np.random.seed(seed)

p = mc.MCPacking(dim,nParticles,np.float64(phiString))
p.setPotentialPower(2)
p.setRandomSeed(seed)
if dim > 6:
	p.setGeometryType(mc.geometryEnum.Dn)

if dim == 8:
	p.setGeometryType(mc.geometryEnum.E8)

p.setMonoRadii()
if biDisperse == 1:
	p.setBiRadii()

if biDisperse == 0 and dim == 3:
	saveDir = basedir + '/' + str(dim) + '-mono/' + str(nParticles) + '/' + phiString + '/' + str(seed) + '/' + zeroFactorString + '/' + maxIterationFactorString + '/' + scaleFactorString + '/'
	startingFolder = basedir + '/' + str(dim) + '-mono/' + str(nParticles) + '/' + phiString + '/' + str(seed) + '/beforeMin/'


p.setRandomPositions()

beforeFolder = saveDir + '/beforeMin/' + phiString
finishFile = saveDir + '/finished.dat'
tempFolder = saveDir + '/temp/'
gardnerFolder = saveDir + '/gardner/'
#Check to see if the directory exists. If so, read the last step and go on,
#if not, do the monte carlo

toRestart = False

neighborCut = 1.0
if not os.path.exists(finishFile):
	if os.path.exists(tempFolder):
		p.load(tempFolder)
		p.calcNeighborsCut(1)
		toRestart = p.getSmallOverlaps().getnnz() != 0
	elif os.path.exists(startingFolder):
		p.load(startingFolder)
		p.calcNeighborsCut(1)
		toRestart = p.getSmallOverlaps().getnnz() != 0
	if (not os.path.exists(startingFolder) or toRestart):
		criticalForceBegin = 1e-12
		p.calcNeighborsCut(neighborCut)
		p.calcForceEnergy()
		print('about to minimize')
		p.minimizeFIRE(criticalForce=criticalForceBegin)
		initialPosition = p.getPositions()
		numSteps, tauAlpha, neighborCutDistance, maxMCDisp = p.thermalizePacking(verbose=True,numTrials=100,numCosts=25,numTau = 3)
		p.save(beforeFolder)
		p.save(startingFolder)
		with open(saveDir + '/tauAlpha.dat', 'w') as file:
			_ = file.write(str(tauAlpha))
		with open(startingFolder + '/tauAlpha.dat', 'w') as file:
			_ = file.write(str(tauAlpha))
		with open(startingFolder + '/neighborCut.dat', 'w') as file:
			_ = file.write(str(neighborCut))
		with open(startingFolder + '/maxMCDisp.dat', 'w') as file:
			_ = file.write(str(maxMCDisp))
		count = 0
	else:
		if not os.path.exists(tempFolder):
			print("loading starting folder")
			p.load(startingFolder)
			p.save(beforeFolder)
			count = 0
			with open(startingFolder + '/tauAlpha.dat', 'r') as file:
				tauAlpha = int(file.read())
			with open(saveDir + '/tauAlpha.dat', 'w') as file:
				_ = file.write(str(tauAlpha))
			with open(startingFolder + '/tauAlpha.dat', 'w') as file:
				_ = file.write(str(tauAlpha))
			count = 0
		else:
			print("loading previous data")
			p.load(tempFolder)
			count = file_len(saveDir + '/phiList.dat')
	if len(glob.glob(gardnerFolder + '/*/positions.dat')) == 0:
		gardnerFound = False
	else:
		gardnerFound = True
	recalcNeighbors(neighborCut,p.getPositions())
	logZero,minGap,_,_ = findLogCutoff(zeroFactor,neighborCut)
	print(logZero)
	neighborCut = 5.0*logZero
	# logZero = np.min([logZero,0.1])
	p.setLogZero(logZero)
	p.calcForceEnergy()
	dtMax = np.sqrt(minGap)
	lastPos = p.getPositions()
	lastPhi = p.getPhi()
	phiList = np.array((),dtype=np.float64)
	criticalForceList = np.array((),dtype=np.float64)
	stdGapList = np.array((),dtype=np.float64)
	minGapList = np.array((),dtype=np.float64)
	continueRunning = True
	recalcNeighbors(neighborCut,p.getPositions())
	p.calcForceEnergy()
	phiJError = 1.0
	lastRecalc = 0
	thisMaxIterations = maxIterations
	count = 0
	# _,phiJError,_ = findJammingGap(hi = 1.0, lo = 0)
	while ((phiJError/p.getPhi() > 2e-5 and count < 10000) or (not gardnerFound)):
		count += 1
		lastRecalc += 1
		oldPhi = p.getPhi()
		dtMax, minGap, normGaps, fireIterations = simpleLogMin(oldPhi = oldPhi, scaleFactor=scaleFactor,zeroFactor=zeroFactor,criticalForce=-1,maxIterations=thisMaxIterations,dt_max=dtMax,dt=dtMax/10.0,saveBeforeExpansion = False, phiJError=phiJError,neighborCut=neighborCut,gardnerFound=gardnerFound)
		sortedGaps = np.sort(normGaps)
		allGaps = sortedGaps[0:int(zeroFactor*nParticles*dim + 1)]
		# if (np.abs(oldPhi - p.getPhi())/oldPhi < 1e-10) or (p.getMaxUnbalancedForce() == 0.0):
		if (np.abs(oldPhi - p.getPhi())/oldPhi < 1e-10):
			break
		with open(saveDir + '/FIREIterations.dat', 'a') as file:
			_ = file.write(str(fireIterations) + '\n')
		with open(saveDir + '/criticalForceList.dat', 'a') as file:
			_ = file.write(str(p.getMaxUnbalancedForce()) + '\n')
		with open(saveDir + '/phiList.dat', 'a') as file:
			_ = file.write(str(p.getPhi()) + '\n')
		with open(saveDir + '/minGapList.dat', 'a') as file:
			_ = file.write(str(minGap) + '\n')
		with open(saveDir + '/meanGap.dat', 'a') as file:
			_ = file.write(str(np.mean(allGaps)) + '\n')
		p.save(tempFolder)
		if fireIterations >= maxIterations and not gardnerFound:
			sortedGaps = np.sort(normGaps)
			allGaps = sortedGaps[0:int(zeroFactor*nParticles*dim + 1)]
			p.save(gardnerFolder + '/' + str(p.getPhi()))
			mc.fileIO.saveScalar(gardnerFolder + '/' + str(p.getPhi()),'zeroFactor',zeroFactor)
			mc.fileIO.saveScalar(gardnerFolder + '/' + str(p.getPhi()),'maxIterations',maxIterations)
			mc.fileIO.saveScalar(gardnerFolder + '/' + str(p.getPhi()),'scaleFactor',scaleFactorString)
			mc.fileIO.saveArray(gardnerFolder + '/allGaps.dat',allGaps)
			with open(gardnerFolder + '/meanGap.dat', 'w') as file:
				_ = file.write(str(np.mean(allGaps)))
			gardnerFound = True
		if saveMiddleSteps:
			p.save(saveDir + '/' + str(p.getPhi()))
		print("FIREIterations is " + str(fireIterations))
		if (gardnerFound and lastRecalc%10 == 0):
			_,phiJError,_ = findJammingGap(hi = 1.0, lo = 0)
		with open(saveDir + '/dtMax.dat', 'w') as file:
			_ = file.write(str(dtMax))
	_,phiJError,_ = findJammingGap(hi = 1.0, lo = 0)
	with open(saveDir + '/phiJError.dat', 'w') as file:
		_ = file.write(str(phiJError))
	if phiJError/p.getPhi() <= np.max([0.5/nParticles, 2e-5]):
		p.save(saveDir + '/' + str(p.getPhi()))
		mc.fileIO.saveScalar(saveDir + '/' + str(p.getPhi()),'zeroFactor',zeroFactor)
		mc.fileIO.saveScalar(saveDir + '/' + str(p.getPhi()),'maxIterations',maxIterations)
		mc.fileIO.saveScalar(saveDir + '/' + str(p.getPhi()),'scaleFactor',scaleFactorString)
		with open(finishFile, 'w') as file:
			_ = file.write(str(1))

