[MMTK] Problems adding configuration to Universe Copy
vanitha@cs.wisc.edu
vanitha at cs.wisc.edu
Wed Feb 8 19:33:10 CET 2006
Sorry for the copy and paste error in the code snippets. Some of the
sections of the code are repeated twice. I wouldn't want to resend this
again to the mailing list but the order is as follows:
1) Termination criteria
2) Choose initial approximation to inverse Hessian
3) Get the descent direction
4) Get the Step size
5) Update the universe
6) Update/Discard stored vectors
7) Update initial energy, gradient and configuration to current one.
Hope this is clear!
Thanks,
- Vanitha
> 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
>
>
>
> _______________________________________________
> mmtk maillist - mmtk at starship.python.net
> http://starship.python.net/mailman/listinfo/mmtk
>
More information about the mmtk
mailing list