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))