diff --git a/setup.py b/setup.py index 1d86313..554425d 100644 --- a/setup.py +++ b/setup.py @@ -14,6 +14,7 @@ url='https://github.com/hpparvi/PyDE', package_dir={'pyde':'src'}, packages=['pyde'], + ext_modules=[Extension('pyde.de_f', ['src/de.f90'], libraries=['gomp','m'])], install_requires=["numpy"], license='GPLv2', classifiers=[ diff --git a/src/de.f90 b/src/de.f90 new file mode 100644 index 0000000..64c0529 --- /dev/null +++ b/src/de.f90 @@ -0,0 +1,62 @@ +module de_f +use omp_lib +implicit none + +contains + subroutine evolve_population(pop, f, c, ndim, npop, pout) + implicit none + integer, intent(in) :: ndim, npop + real(8), intent(in) :: f, c + real(8), intent(in), dimension(npop,ndim) :: pop + real(8), intent(out), dimension(npop,ndim) :: pout + integer :: i + + ! $omp parallel do shared(pop,f,c,ndim,npop,pout) private(i) + do i=0,npop-1 + call evolve_pv(i, pop, f, c, ndim, npop, pout(i+1,:)) + end do + ! $omp end parallel do + end subroutine evolve_population + + + subroutine evolve_pv(i, pop, f, c, ndim, npop, vout) + implicit none + integer, intent(in) :: i, ndim, npop + real(8), intent(in) :: f, c + real(8), intent(in), dimension(npop,ndim) :: pop + real(8), intent(out), dimension(ndim) :: vout + integer :: j, ids(3) + real(8), dimension(ndim) :: co, v + real(8) :: x + + ids = i+1 + do while (ids(1) == i+1) + call random_number(x) + ids(1) = floor(x*npop) + 1 + end do + do while ((ids(2) == i+1) .or. (ids(2) == ids(1))) + call random_number(x) + ids(2) = floor(x*npop) + 1 + end do + do while ((ids(3) == i+1) .or. (ids(3) == ids(2)) .or. (ids(3) == ids(1))) + call random_number(x) + ids(3) = floor(x*npop) + 1 + end do + + v = pop(ids(1),:) + f * (pop(ids(2),:) - pop(ids(3),:)) + + !! --- CROSS OVER --- + call random_number(co) + where (co <= c) + vout = v + elsewhere + vout = pop(i+1,:) + end where + + !! --- FORCED CROSSING -- - + call random_number(x) + j = floor(x*ndim) + 1 + vout(j) = v(j) + + end subroutine evolve_pv +end module de_f diff --git a/src/de.py b/src/de.py index fda1271..c8e077f 100644 --- a/src/de.py +++ b/src/de.py @@ -8,6 +8,7 @@ import numpy as np from numpy.random import random, randint +from de_f import de_f class DiffEvol(object): """ @@ -36,8 +37,9 @@ class DiffEvol(object): :param maximize: (optional) Switch setting whether to maximize or minimize the function. Defaults to minimization. """ - def __init__(self, fun, bounds, npop, F=0.5, C=0.5, seed=0, maximize=False): - np.random.seed(seed) + def __init__(self, fun, bounds, npop, F=0.5, C=0.5, seed=None, maximize=False, vfun=False, cbounds=[0.25, 1]): + if seed: + np.random.seed(seed) self.minfun = fun self.bounds = np.asarray(bounds) @@ -47,16 +49,27 @@ def __init__(self, fun, bounds, npop, F=0.5, C=0.5, seed=0, maximize=False): self.bw = np.tile(self.bounds[:,1]-self.bounds[:,0],[npop,1]) self.m = -1 if maximize else 1 + self.cmin = cbounds[0] + self.cmax = cbounds[1] + self.seed = seed self.F = F self.C = C - self._population = self.bl + random([self.n_pop, self.n_par]) * self.bw + self._population = np.asarray(self.bl + random([self.n_pop, self.n_par]) * self.bw) self._fitness = np.zeros(npop) self._minidx = None + self._trial_pop = np.zeros_like(self._population) + self._trial_fit = np.zeros_like(self._fitness) + + if vfun: + self._eval = self._eval_vfun + else: + self._eval = self._eval_sfun + @property - def population(self): + def population(self): """The parameter vector population""" return self._population @@ -82,36 +95,49 @@ def optimize(self, ngen): return res def __call__(self, ngen=1): - t = np.zeros(3, np.int) + return self._eval(ngen) + + def _eval_sfun(self, ngen=1): + """Run DE for a function that takes a single pv as an input and retuns a single value.""" + popc, fitc = self._population, self._fitness + popt, fitt = self._trial_pop, self._trial_fit - for i in xrange(self.n_pop): - self._fitness[i] = self.m * self.minfun(self._population[i,:]) - - for j in xrange(ngen): - for i in xrange(self.n_pop): - t[:] = i - while t[0] == i: - t[0] = randint(self.n_pop) - while t[1] == i or t[1] == t[0]: - t[1] = randint(self.n_pop) - while t[2] == i or t[2] == t[0] or t[2] == t[1]: - t[2] = randint(self.n_pop) + for ipop in xrange(self.n_pop): + fitc[ipop] = self.m * self.minfun(popc[ipop,:]) + + for igen in xrange(ngen): + popt[:,:] = de_f.evolve_population(popc, self.F, self.C) + + for ipop in xrange(self.n_pop): + fitt[ipop] = self.m * self.minfun(popt[ipop,:]) - v = self._population[t[0],:] + self.F * (self._population[t[1],:] - self._population[t[2],:]) + msk = fitt < fitc + popc[msk,:] = popt[msk,:] + fitc[msk] = fitt[msk] - ## --- CROSS OVER --- - crossover = random(self.n_par) <= self.C - u = np.where(crossover, v, self._population[i,:]) + self._minidx = np.argmin(fitc) + yield popc[self._minidx,:], fitc[self._minidx] - ## --- FORCED CROSSING --- - ri = randint(self.n_par) - u[ri] = v[ri].copy() - ufit = self.m * self.minfun(u) - - if ufit < self._fitness[i]: - self._population[i,:] = u[:].copy() - self._fitness[i] = ufit + def _eval_vfun(self, ngen=1): + """Run DE for a function that takes the whole population as an input and retuns a value for each pv.""" + popc, fitc = self._population, self._fitness + popt, fitt = self._trial_pop, self._trial_fit + + fitc[:] = self.m * self.minfun(self._population) + + for igen in xrange(ngen): + x = float(ngen-igen)/float(ngen) + + self.F = np.random.uniform(0.25,0.75) + self.C = x*self.cmax + (1-x)*self.cmin + + popt[:,:] = de_f.evolve_population(popc, self.F, self.C) + fitt[:] = self.m * self.minfun(popt) + msk = fitt < fitc + popc[msk,:] = popt[msk,:] + fitc[msk] = fitt[msk] + + self._minidx = np.argmin(fitc) + yield popc[self._minidx,:], fitc[self._minidx] - self._minidx = np.argmin(self._fitness) - yield self.population[self._minidx,:], self._fitness[self._minidx]