Python optimize.LBFGS_CPP类代码示例

 def test_run_niter2(self):
     lbfgs1 = LBFGS_CPP(_xrand, _EG())
     res1 = lbfgs1.run()
     lbfgs2 = LBFGS_CPP(_xrand, _EG())
     res2 = lbfgs2.run(res1.nsteps / 2)
     res2 = lbfgs2.run()
     self.assert_same(res1, res2)

 def test_run_niter3(self):
     lbfgs1 = LBFGS_CPP(_xrand, _EG())
     res1 = lbfgs1.run(10)
     lbfgs2 = LBFGS_CPP(_xrand, _EG())
     res2 = lbfgs2.run(5)
     res2 = lbfgs2.run(5)
     self.assert_same(res1, res2)

 def optimize(self, nr_samples = 1):
     self.optimizer = ModifiedFireCPP(self.x_initial, self.potential, tol = self.tol, maxstep = self.maxstep)
     self.optimizer_ = LBFGS_CPP(self.x_initial, self.potential)
     self.optimizer_cells = ModifiedFireCPP(self.x_initial, self.potential_cells, tol = self.tol, maxstep = self.maxstep)
     self.optimizer_cells_ = LBFGS_CPP(self.x_initial, self.potential_cells)
     t0 = time.time()
     self.optimizer.run(self.nstepsmax)
     t1 = time.time()
     self.optimizer_cells.run(self.nstepsmax)
     t2 = time.time()
     self.optimizer_.run(self.nstepsmax)
     t3 = time.time()
     self.optimizer_cells_.run(self.nstepsmax)
     t4 = time.time()
     res_x_final = self.optimizer.get_result()
     res_x_final_cells = self.optimizer_cells.get_result()
     self.x_final = res_x_final.coords
     self.x_final_cells = res_x_final_cells.coords
     print "number of particles:", self.N
     print "time no cell lists:", t1 - t0, "sec"
     print "time cell lists:", t2 - t1, "sec"
     print "ratio:", (t1 - t0) / (t2 - t1)
     assert(res_x_final.success)
     assert(res_x_final_cells.success)
     for (xci, xi) in zip(self.x_final_cells, self.x_final):
         passed = (np.abs(xci - xi) < 1e-10 * self.N)
         if (passed is False):
             print "xci", xci
             print "xi", xi
             assert(passed)
     print "energy no cell lists:", res_x_final.energy
     print "energy cell lists:", res_x_final_cells.energy
     self.t_ratio = (t1 - t0) / (t2 - t1)
     self.t_ratio_lbfgs = (t3 - t2) / (t4 - t3)

 def do_check(self, pot, **kwargs):
     lbfgs = LBFGS_CPP(_xrand, pot, **kwargs)
     res = lbfgs.run()
     self.assertAlmostEqual(res.energy, _emin, 4)
     self.assertTrue(res.success)
     self.assertLess(np.max(np.abs(res.coords - _xmin)), 1e-2)
     self.assertGreater(res.nfev, 0)

def lbfgs_cpp(coords, pot, **kwargs):
    if not hasattr(pot, "getEnergyGradient"):
        # for compatibility with old quenchers.
        # assume pot is a getEnergyGradient function
        pot = _getEnergyGradientWrapper(pot)
    lbfgs = LBFGS_CPP(coords, pot, **kwargs)    
    return lbfgs.run()

class Config2D(object):
    def __init__(self, nparticles_x, amplitude):
        self.LX = nparticles_x
        self.LY = self.LX
        self.nparticles_x = nparticles_x
        self.N = self.nparticles_x ** 2
        self.amplitude = amplitude
        self.x = np.zeros(2 * self.N)
        for particle in xrange(self.N):
            pid = 2 * particle
            self.x[pid] = particle % self.LX
            self.x[pid + 1] = int(particle / self.LX)
        self.x_initial = [xi + np.random.uniform(- self.amplitude, self.amplitude) for xi in self.x]
        self.radius = 0.25
        self.sca = 1.1
        self.radii = np.ones(self.N) * self.radius
        self.eps = 1
        self.boxvec = np.array([self.LX, self.LY])
        self.potential = HS_WCA(self.eps, self.sca, self.radii, ndim=2, boxvec=self.boxvec)
        self.rcut = 2 * (1 + self.sca) * self.radius
        self.ncellx_scale = 1
        self.potential_cells = HS_WCAPeriodicCellLists(self.eps, self.sca, self.radii, self.boxvec, self.x_initial, self.rcut, ndim = 2, ncellx_scale = self.ncellx_scale)
        self.tol = 1e-7
        self.maxstep = 1
        self.nstepsmax = 1e5
    def optimize(self, nr_samples = 1):
        self.optimizer = ModifiedFireCPP(self.x_initial, self.potential, tol = self.tol, maxstep = self.maxstep)
        self.optimizer_ = LBFGS_CPP(self.x_initial, self.potential)
        self.optimizer_cells = ModifiedFireCPP(self.x_initial, self.potential_cells, tol = self.tol, maxstep = self.maxstep)
        self.optimizer_cells_ = LBFGS_CPP(self.x_initial, self.potential_cells)
        t0 = time.time()
        self.optimizer.run(self.nstepsmax)
        t1 = time.time()
        self.optimizer_cells.run(self.nstepsmax)
        t2 = time.time()
        self.optimizer_.run(self.nstepsmax)
        t3 = time.time()
        self.optimizer_cells_.run(self.nstepsmax)
        t4 = time.time()
        res_x_final = self.optimizer.get_result()
        res_x_final_cells = self.optimizer_cells.get_result()
        self.x_final = res_x_final.coords
        self.x_final_cells = res_x_final_cells.coords
        print "number of particles:", self.N
        print "time no cell lists:", t1 - t0, "sec"
        print "time cell lists:", t2 - t1, "sec"
        print "ratio:", (t1 - t0) / (t2 - t1)
        assert(res_x_final.success)
        assert(res_x_final_cells.success)
        for (xci, xi) in zip(self.x_final_cells, self.x_final):
            passed = (np.abs(xci - xi) < 1e-10 * self.N)
            if (passed is False):
                print "xci", xci
                print "xi", xi
                assert(passed)
        print "energy no cell lists:", res_x_final.energy
        print "energy cell lists:", res_x_final_cells.energy
        self.t_ratio = (t1 - t0) / (t2 - t1)
        self.t_ratio_lbfgs = (t3 - t2) / (t4 - t3)

    def test_event_raise(self):
        class EventException(BaseException): pass

        def myevent(*args, **kwargs): raise EventException

        with self.assertRaises(EventException):
            lbfgs = LBFGS_CPP(_xrand, _EG(), events=[myevent])
            lbfgs.run()

 def test_result(self):
     lbfgs = LBFGS_CPP(_xrand, _EG())
     res = lbfgs.one_iteration()
     self.assertIn("H0", res)
     self.assertIn("energy", res)
     self.assertIn("grad", res)
     self.assertIn("success", res)
     self.assertIn("coords", res)
     self.assertIn("rms", res)
     self.assertIn("nsteps", res)
     self.assertIn("nfev", res)

def mytest():
    system = test_functions.BealeSystem()
    print "do pot"
    pot = system.get_potential()
    print "done pot"
    x = pot.target_coords.copy() 
    x += np.random.uniform(-0.2, 0.2, x.shape)
    from pele.optimize import LBFGS_CPP
    lbfgs = LBFGS_CPP(x, pot, verbosity=100)
    print "done setting up"
    res = lbfgs.run()
    res = _quench.lbfgs_cpp(x, pot, verbosity=100)
#    print res
    print res

 def __init__(self, boxdim=2, nr_particles=100, hard_phi=0.4,
              nr_steps=1e6, epsilon=1, alpha=0.1, verbose=False):
     # Settings.
     np.random.seed(42)
     # Input parameters.
     self.boxdim = boxdim
     self.nr_particles = nr_particles
     self.hard_phi = hard_phi
     self.nr_steps = nr_steps
     self.epsilon = epsilon
     self.alpha = alpha
     self.verbose = verbose
     # Derived quantities.
     self.hard_radii = np.ones(self.nr_particles)
     def volume_nball(radius, n):
         return np.power(np.pi, n / 2) * np.power(radius, n) / gamma(n / 2 + 1)
     self.box_length = np.power(np.sum(np.asarray([volume_nball(r, self.boxdim) for r in self.hard_radii])) / self.hard_phi, 1 / self.boxdim)
     self.box_vector = np.ones(self.boxdim) * self.box_length
     # HS-WCA potential.
     self.potential = HS_WCA(use_periodic=True, use_cell_lists=True,
                             ndim=self.boxdim, eps=self.epsilon,
                             sca=self.alpha, radii=self.hard_radii,
                             boxvec=self.box_vector)
     # Initial configuration by minimization.
     self.nr_dof = self.boxdim * self.nr_particles
     self.x = np.random.uniform(-0.5 * self.box_length, 0.5 * self.box_length, self.nr_dof)
     optimizer = LBFGS_CPP(self.x, self.potential)
     optimizer.run()
     if not optimizer.get_result().success:
         print ("warning: minimization has not converged")
     self.x = optimizer.get_result().coords.copy()
     # Potential and MC rules.
     self.temperature = 1
     self.mc = MC(self.potential, self.x, self.temperature, self.nr_steps)
     self.step = RandomCoordsDisplacement(42, 1, single=True, nparticles=self.nr_particles, bdim=self.boxdim)
     if self.verbose:
         print ("initial MC stepsize")
         print self.step.get_stepsize()
     self.mc.set_takestep(self.step)
     self.eq_steps = self.nr_steps / 2
     self.mc.set_report_steps(self.eq_steps)
     self.gr_quench = RecordPairDistHistogram(self.box_vector, 50, self.eq_steps, self.nr_particles, optimizer=optimizer)
     self.gr = RecordPairDistHistogram(self.box_vector, 50, self.eq_steps, self.nr_particles)
     self.mc.add_action(self.gr_quench)
     self.mc.add_action(self.gr)
     self.test = MetropolisTest(44)
     self.mc.add_accept_test(self.test)

class MinimizeUniformHardsSpheres(object):
    def __init__(self, nr_particles=42, hard_volume_fraction=0.5, epsilon=1, alpha=0.2):
        np.random.seed(42)
        self.nr_particles = nr_particles
        self.hard_volume_fraction = hard_volume_fraction
        self.epsilon = epsilon
        self.alpha = alpha
        self.hard_radii = np.random.normal(loc=1, scale=0.1, size=self.nr_particles)
        self.box_length = np.power(np.sum(np.asarray([4 * np.pi * r**3 / 3 for r in self.hard_radii])) / self.hard_volume_fraction, 1/3)
        self.nr_dof = 3 * self.nr_particles
        self.x = np.random.uniform(-0.5 * self.box_length, 0.5 * self.box_length, self.nr_dof)
        self.box_vector = np.ones(3) * self.box_length
        self.rcut = 2 * (1 + alpha) * np.amax(self.hard_radii)
        self.potential = HS_WCA(use_periodic=True, use_cell_lists=False, eps=self.epsilon, sca=self.alpha, radii=self.hard_radii, boxvec=self.box_vector, rcut=self.rcut)
        self.optimizer = LBFGS_CPP(self.x, self.potential)
        print "energy before:", self.potential.getEnergy(self.x)
        self.optimizer.run()
        print "minimization converged", self.optimizer.get_result().success
        print "energy after:", self.potential.getEnergy(self.optimizer.get_result().coords)

 def test_normal_estimates_OK(self):
     observed = np.sqrt(self.ss) * randn(self.nr_points) + self.mu
     pot = MLCost(observed, probf=gauss)
     potl = MLCost(observed, log_probf=log_gauss)
     parameters = [self.mu + 1, self.ss - 2]
     optimizer = LBFGS_CPP(parameters, pot)
     optimizerl = LBFGS_CPP(parameters, potl)
     result = optimizer.run()
     resultl = optimizerl.run()
     opt_parameters = result.coords
     opt_parametersl = resultl.coords
     opt_mu = opt_parameters[0]
     opt_mul = opt_parametersl[0]
     opt_ss = opt_parameters[1]
     opt_ssl = opt_parametersl[1]
     self.assertAlmostEqual(opt_mu, self.mu, delta=2*self.mu/np.sqrt(self.nr_points))
     self.assertAlmostEqual(opt_ss, self.ss, delta=2*self.ss/np.sqrt(self.nr_points))
     self.assertAlmostEqual(opt_mu, opt_mul, delta = 1e-4)
     self.assertAlmostEqual(opt_ss, opt_ssl, delta = 1e-4)
     confidence_intervalsl = potl.get_error_estimate(opt_parametersl)
     for i, par in enumerate(opt_parameters):
         self.assertLessEqual(confidence_intervalsl[i][0], par)
         self.assertLessEqual(par, confidence_intervalsl[i][1])

 def test_raises(self):
     with self.assertRaises(NotImplementedError):
         lbfgs = LBFGS_CPP(_xrand, _Raise())
         lbfgs.run()

 def optimize(self, nr_samples=1):
     self.x_initial_red = reduce_coordinates(self.x_initial,self.frozen_atoms,self.ndim)
     self.optimizer = ModifiedFireCPP(self.x_initial_red.copy(), self.potential, tol = self.tol, maxstep = self.maxstep)
     self.optimizer_ = LBFGS_CPP(self.x_initial_red.copy(), self.potential)
     self.optimizer_cells = ModifiedFireCPP(self.x_initial_red.copy(), self.potential_cells, tol = self.tol, maxstep = self.maxstep)
     self.optimizer_cells_ = LBFGS_CPP(self.x_initial_red.copy(), self.potential_cells)
     t0 = time.time()
     print "self.optimizer.run(self.nstepsmax)", self.nstepsmax
     self.optimizer.run(self.nstepsmax)
     self.res_x_final = self.optimizer.get_result()
     t1 = time.time()
     self.optimizer_cells.run(self.nstepsmax)
     self.res_x_final_cells = self.optimizer_cells.get_result()
     t2 = time.time()
             
     self.x_final = self.res_x_final.coords
     self.x_final_cells = self.res_x_final_cells.coords
     print "fire final E, x:", self.optimizer.get_result().energy
     print "fire final E, x_cells:", self.optimizer_cells.get_result().energy
     print "fire final E, plain: ", self.potential.getEnergy(self.x_final)
     print "fire final E, cell: ", self.potential_cells.getEnergy(self.x_final_cells)
     print "fire number of particles:", self.N
     print "fire time no cell lists:", t1 - t0, "sec"
     print "fire time cell lists:", t2 - t1, "sec"
     print "fire ratio:", (t1 - t0) / (t2 - t1)
     
     if not self.res_x_final.success or not self.res_x_final_cells.success:
         print "-------------"
         print "res_x_final.rms:", self.res_x_final.rms
         print "res_x_final.nfev:", self.res_x_final.nfev
         print "res_x_final_cells.rms:", self.res_x_final_cells.rms
         print "res_x_final_cells.nfev:", self.res_x_final_cells.nfev
         print "self.res_x_final.success", self.res_x_final.success
         print "self.res_x_final_cells.success", self.res_x_final_cells.success
         print "-------------"
         plot_disks(self.x_initial, self.radii, self.boxvec, sca=self.sca)
         plot_disks(self.x_final, self.radii, self.boxvec, sca=self.sca)
         plot_disks(self.x_final_cells, self.radii, self.boxvec, sca=self.sca)
             
     self.optimizer_.run(self.nstepsmax)
     self.res_x_final_ = self.optimizer_.get_result()
     t3 = time.time()
     self.optimizer_cells_.run(self.nstepsmax)
     self.res_x_final_cells_ = self.optimizer_cells_.get_result()
     t4 = time.time()
     
     self.x_final_ = self.res_x_final_.coords
     self.x_final_cells_ = self.res_x_final_cells_.coords
     print "lbfgs final E, x:", self.optimizer_.get_result().energy
     print "lbfgs final E, x_cells:", self.optimizer_cells_.get_result().energy
     print "lbfgs final E, plain: ", self.potential.getEnergy(self.x_final_)
     print "lbfgs final E, cell: ", self.potential_cells.getEnergy(self.x_final_cells_)
     print "lbfgs number of particles:", self.N
     print "lbfgs time no cell lists:", t3 - t2, "sec"
     print "lbfgs time cell lists:", t4 - t3, "sec"
     print "lbfgs ratio:", (t3 - t2) / (t4 - t3)
     
     if not self.res_x_final_.success or not self.res_x_final_cells_.success or not self.res_x_final.success or not self.res_x_final_cells.success:
         print "-------------"
         print "res_x_final_.rms:", self.res_x_final_.rms
         print "res_x_final_.nfev:", self.res_x_final_.nfev
         print "res_x_final_cells_.rms:", self.res_x_final_cells_.rms
         print "res_x_final_cells_.nfev:", self.res_x_final_cells_.nfev
         print "self.res_x_final_.success", self.res_x_final.success
         print "self.res_x_final_cells_.success", self.res_x_final_cells.success
         print "-------------"
         plot_disks(self.x_initial, self.radii, self.boxvec, sca=self.sca)
         plot_disks(self.x_final_, self.radii, self.boxvec, sca=self.sca)
         plot_disks(self.x_final_cells_, self.radii, self.boxvec, sca=self.sca)
     
     assert(self.res_x_final.success)
     assert(self.res_x_final_cells.success)
     assert(self.res_x_final_.success)
     assert(self.res_x_final_cells_.success)
     for (xci, xi) in zip(self.x_final_cells, self.x_final):
         passed = (np.abs(xci - xi) < 1e-10)
         if (passed is False):
             print "xci", xci
             print "xi", xi
             assert(passed)
     print "energy no cell lists:", self.res_x_final.energy
     print "energy cell lists:", self.res_x_final_cells.energy
     self.t_ratio = (t1 - t0) / (t2 - t1)
     self.t_ratio_lbfgs = (t3 - t2) / (t4 - t3)

class Config2DFrozenBoundary(object):
    def __init__(self, nparticles_x, amplitude):
        self.ndim = 2
        self.LX = nparticles_x
        self.LY = self.LX
        self.nparticles_x = nparticles_x
        self.N = self.nparticles_x ** self.ndim
        self.amplitude = amplitude
        self.dof = self.ndim * self.N
        self.x = np.zeros(self.dof)
        self.frozen_atoms = []
        for particle in xrange(self.N):
            pid = self.ndim * particle
            xcoor = particle % self.LX
            ycoor = int(particle / self.LX)
            self.x[pid] = xcoor
            self.x[pid + 1] = ycoor
            if xcoor == 0 or xcoor == self.LX - 1 or ycoor == 0 or ycoor == self.LY - 1:
                self.frozen_atoms.append(particle)
        self.x_initial = copy.copy(self.x)
        for particle in xrange(self.N):
            if particle not in self.frozen_atoms:
                pid = self.ndim * particle
                self.x_initial[pid] += np.random.uniform(- self.amplitude, self.amplitude)
                self.x_initial[pid + 1] += np.random.uniform(- self.amplitude, self.amplitude)
        self.x_initial = np.reshape(self.x_initial, (self.N,2))
        self.x_initial[:,0] -= np.mean(self.x_initial[:,0])
        self.x_initial[:,1] -= np.mean(self.x_initial[:,1])
        self.x_initial = self.x_initial.flatten()
        min_x = np.amin(self.x_initial)
        if min_x < 0:
            self.x_initial -= min_x
        #self.radius = 0.3
        #self.sca = 1.5
        self.radius = 0.25
        self.sca = 1.8
        self.radii = np.ones(self.N) * self.radius
        self.eps = 1.0
        max_edge = np.amax([np.amax(self.x_initial), np.abs(np.amin(self.x_initial))]) + 2 * self.amplitude + (1 + self.sca) * self.radius
        self.boxvec = np.array([max_edge, max_edge])
        self.frozen_atoms1 = np.array(self.frozen_atoms)
        self.frozen_atoms2 = np.array(self.frozen_atoms)
        print "self.frozen_atoms1", self.frozen_atoms1
        self.potential = HS_WCA(use_frozen=True, use_periodic=use_periodic_frozen,
                         reference_coords=self.x_initial,
                         frozen_atoms=self.frozen_atoms1,
                         eps=self.eps, sca=self.sca, radii=self.radii,
                         ndim=self.ndim, boxvec=self.boxvec)
        self.rcut =  2 * (1 + self.sca) * self.radius
        self.ncellx_scale = 1.0
        self.potential_cells = HS_WCA(use_frozen=True,
                               use_periodic=use_periodic_frozen, use_cell_lists=True,
                               eps=self.eps, sca=self.sca,
                               radii=self.radii, boxvec=self.boxvec,
                               reference_coords=self.x_initial,
                               rcut=self.rcut, ndim=self.ndim,
                               ncellx_scale=self.ncellx_scale,
                               frozen_atoms=self.frozen_atoms2)
        self.tol = 1e-7
        self.maxstep = np.amax(self.radii)
        self.nstepsmax = int(1e6)
        assert(self.boxvec[0]==self.boxvec[1])
        self.x_initial_red = reduce_coordinates(self.x_initial,self.frozen_atoms,self.ndim)
        print "x_initial energy:", self.potential.getEnergy(self.x_initial_red)
        print "x_initial cells energy:", self.potential_cells.getEnergy(self.x_initial_red)
        #assert abs(self.potential.getEnergy(self.x_initial_red) - self.potential_cells.getEnergy(self.x_initial_red)) < 1e-10
        assert np.allclose(self.potential.getEnergy(self.x_initial_red), self.potential_cells.getEnergy(self.x_initial_red), rtol=1e-10)
        print self.boxvec
    
    def optimize(self, nr_samples=1):
        self.x_initial_red = reduce_coordinates(self.x_initial,self.frozen_atoms,self.ndim)
        self.optimizer = ModifiedFireCPP(self.x_initial_red.copy(), self.potential, tol = self.tol, maxstep = self.maxstep)
        self.optimizer_ = LBFGS_CPP(self.x_initial_red.copy(), self.potential)
        self.optimizer_cells = ModifiedFireCPP(self.x_initial_red.copy(), self.potential_cells, tol = self.tol, maxstep = self.maxstep)
        self.optimizer_cells_ = LBFGS_CPP(self.x_initial_red.copy(), self.potential_cells)
        t0 = time.time()
        print "self.optimizer.run(self.nstepsmax)", self.nstepsmax
        self.optimizer.run(self.nstepsmax)
        self.res_x_final = self.optimizer.get_result()
        t1 = time.time()
        self.optimizer_cells.run(self.nstepsmax)
        self.res_x_final_cells = self.optimizer_cells.get_result()
        t2 = time.time()
                
        self.x_final = self.res_x_final.coords
        self.x_final_cells = self.res_x_final_cells.coords
        print "fire final E, x:", self.optimizer.get_result().energy
        print "fire final E, x_cells:", self.optimizer_cells.get_result().energy
        print "fire final E, plain: ", self.potential.getEnergy(self.x_final)
        print "fire final E, cell: ", self.potential_cells.getEnergy(self.x_final_cells)
        print "fire number of particles:", self.N
        print "fire time no cell lists:", t1 - t0, "sec"
        print "fire time cell lists:", t2 - t1, "sec"
        print "fire ratio:", (t1 - t0) / (t2 - t1)
        
        if not self.res_x_final.success or not self.res_x_final_cells.success:
            print "-------------"
            print "res_x_final.rms:", self.res_x_final.rms
            print "res_x_final.nfev:", self.res_x_final.nfev
            print "res_x_final_cells.rms:", self.res_x_final_cells.rms
#.........这里部分代码省略.........