[MMTK] Problems adding configuration to Universe Copy
vanitha@cs.wisc.edu
vanitha at cs.wisc.edu
Wed Feb 8 19:18:47 CET 2006
I've written an implementation of the Limited Memory BFGS minimization
algorithm in Python within the MMTK framework. I am running into a problem
when I create a universe, figure out a "search direction" and then I add
this to a copy of my universe. I get the following problem when I try to
do this. (see traceback and code below). How can I resolve this? I can
repeat code all over again, but that really won't help. I'm trying to see
how I can avoid computing the search direction twice, once with the
original universe and then with a copy of the universe.
Any help will be greatly appreciated!
Thanks,
- Vanitha
Traceback (most recent call last):
File "testlbfgs.py", line 13, in ?
Minimization.LBFGSMinimize(universe, 7, 1000, 0.0001)
File
"/afs/cs.wisc.edu/u/v/a/vanitha/private/MMTK-2.5.11/Examples/Minimization/Test3/Minimization.py",
line 105, in LBFGSMinimize
alfa = backtrackingLineSearch( universeCopy, initialAlfa, direction,
maxIter)
File
"/afs/cs.wisc.edu/u/v/a/vanitha/private/MMTK-2.5.11/Examples/Minimization/Test3/Minimization.py",
line 147, in backtrackingLineSearch
Dphi0 = gradient.dotProduct(direction)
File
"/u/v/a/vanitha/lib/python2.4/site-packages/MMTK/ParticleProperties.py",
line 286, in dotProduct
raise ValueError('Variables are for different universes')
ValueError: Variables are for different universes
def LBFGSMinimize( universe, m, maxIter, tolerance):
# Iteration counter
k = 0
# The current configuration of the universe is the initial guess
energy, gradient = universe.energyAndGradients()
configuration = universe.copyConfiguration()
while 1:
print energy, gradient.norm()
###############################################################
# Termination Criteria
###############################################################
if( k >= maxIter):
break
if( gradient.norm() <= tolerance):
break
################################################################
# Choose the initial approximation to the inverse of the
Hessian
# This is given by Hk0 = gamma*I
# The following steps compute the value of gamma
################################################################
if( k == 0):
# Initialize the rho, s and y vectors
rho = []
s = []
y = []
gamma = 1
else:
skMinusm = s[len(s)-1]
ykMinusm = y[len(y)-1]
# Skip the update if negative curvature
# was detected in the previous step.
#if( skMinusm.dotProduct(ykMinusm) > 0):
gammaTop = skMinusm.dotProduct(ykMinusm)
gammaBottom = ykMinusm.dotProduct(ykMinusm)
gamma = gammaTop/gammaBottom
####################################################
# Get the descent direction
####################################################
universeCopy = copy(universe)
energy, gradient = universeCopy.energyAndGradients()
q = gradient
alpha = []
# Start from len(s)-1 in steps of -1 until 0
for i in range(len(s)-1, -1, -1):
alpha.append(rho[i]*(s[i].dotProduct(q)))
q = q - alpha[len(s)-1-i]*y[i];
alpha.reverse()
r = gamma*q
while 1:
print energy, gradient.norm()
###############################################################
# Termination Criteria
###############################################################
if( k >= maxIter):
break
if( gradient.norm() <= tolerance):
break
################################################################
# Choose the initial approximation to the inverse of the
Hessian
# This is given by Hk0 = gamma*I
# The following steps compute the value of gamma
################################################################
if( k == 0):
# Initialize the rho, s and y vectors
rho = []
s = []
y = []
gamma = 1
else:
skMinusm = s[len(s)-1]
ykMinusm = y[len(y)-1]
# Skip the update if negative curvature
# was detected in the previous step.
#if( skMinusm.dotProduct(ykMinusm) > 0):
gammaTop = skMinusm.dotProduct(ykMinusm)
gammaBottom = ykMinusm.dotProduct(ykMinusm)
gamma = gammaTop/gammaBottom
####################################################
# Get the descent direction
####################################################
universeCopy = copy(universe)
energy, gradient = universeCopy.energyAndGradients()
q = gradient
alpha = []
# Start from len(s)-1 in steps of -1 until 0
for i in range(len(s)-1, -1, -1):
alpha.append(rho[i]*(s[i].dotProduct(q)))
#if( skMinusm.dotProduct(ykMinusm) > 0):
gammaTop = skMinusm.dotProduct(ykMinusm)
gammaBottom = ykMinusm.dotProduct(ykMinusm)
gamma = gammaTop/gammaBottom
####################################################
# Get the descent direction
####################################################
universeCopy = copy(universe)
energy, gradient = universeCopy.energyAndGradients()
q = gradient
alpha = []
# Start from len(s)-1 in steps of -1 until 0
for i in range(len(s)-1, -1, -1):
alpha.append(rho[i]*(s[i].dotProduct(q)))
q = q - alpha[len(s)-1-i]*y[i];
alpha.reverse()
r = gamma*q
# Start from 0 in steps of 1 until len(s)-1
for i in range(0, len(s)):
beta = rho[i]*(y[i].dotProduct(r))
r = r + s[i]*(alpha[i] - beta)
####################################################
# Get the step size
####################################################
initialAlfa = 1
alfa = initialAlfa
direction = -r
alfa = backtrackingLineSearch( universeCopy, initialAlfa,
direction, maxIter)
####################################################
# Update the universe
####################################################
universe.addToConfiguration(alfa*direction)
newEnergy, newGradient = universe.energyAndGradients()
newConfiguration = universe.copyConfiguration()
####################################################
# Update/discard the stored vectors
####################################################
# Discard the first vector pair from storage
if( k >= m):
del s[0]
del y[0]
del rho[0]
s.append(newConfiguration - configuration)
y.append(newGradient - gradient)
totalStoredVectors = len(y)
val =
y[totalStoredVectors-1].dotProduct(s[totalStoredVectors-1])
rho.append(1/val)
#######################################################
# Update the initial gradient, energy and configuration
# to the current one!
#######################################################
gradient = newGradient
energy = newEnergy
configuration = newConfiguration
k=k+1
def backtrackingLineSearch( universeCopy, initialAlfa, direction, maxIter):
# Initialization
rho = 0.05
c = 0.9
maxit = 100
energy, gradient = universeCopy.energyAndGradients()
Dphi0 = gradient.dotProduct(direction)
phi0 = energy
alfa = initialAlfa
if ((alfa <= 0) or (Dphi0 > 0)):
print "Initialization of step incorrect!"
for i in range(maxit):
universeCopy.addToConfiguration(alfa*direction)
new_energy, new_gradient = universeCopy.energyAndGradients()
#print "New Energy=", new_energy
#print "New gradient=", new_gradient.norm()
if(new_energy >= phi0 + alfa*c*Dphi0):
universeCopy.addToConfiguration(-alfa*direction)
alfa = rho*alfa
else:
break
# Reset the value of alfa back to its original value
# in case the max iterations have been reached and if
# no decrease in function value could be obtained with alfa.
if( (new_energy >= phi0 + alfa*c*Dphi0) & (i >= maxit) ):
alfa = initialAlfa
#print "alfa =", alfa
return alfa
More information about the mmtk
mailing list