API tools faq
paste
Advertisement
Guest User

PSO Layara

a guest
Nov 14th, 2019
91
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 10.67 KB | None | 0 0
raw download clone embed print report
  1. import numpy as np
  2. import time  as t
  3.  
  4. from   numpy.random import random as np_rand
  5. from   matplotlib   import pyplot as plt
  6.  
  7. class GBest:
  8.     """
  9.    Serve para contar num_fx.
  10.    Serve para definir qual o g_best atual.
  11.    Serve para atualizar os grupos de partículas (grupo = partícula central + orbitais).
  12.    Atualizar significa que se tiver alguma melhor orbital melhor que a central do grupo dela, a orbital que é boa será a nova central."""
  13.  
  14.     def __init__(self, num_fx, max_fx):
  15.         self.g_best_coords  = None
  16.         self.g_best_fitness = +float('inf')
  17.  
  18.         self.num_fx = num_fx
  19.         self.max_fx = max_fx
  20.    
  21.     def updateGBest(self, coords, fitness):
  22.         "Atualiza g best."
  23.         if fitness < self.g_best_fitness:
  24.             self.g_best_coords  = coords.copy()
  25.             self.g_best_fitness = fitness.copy()
  26.    
  27.     def updateGroups(self, particulas):
  28.         "Vai percorrer todas as partículas e atualizar todos os grupos de centrais + orbitais."
  29.         global num_orbitais
  30.         global alpha_orbitais, alpha_centrais
  31.  
  32.         old_centrais = [particula for particula in particulas if type(particula) == ParticulaPrincipal]
  33.         index_old_centrais = [particulas.index(particula) for particula in old_centrais]
  34.  
  35.         nova_pop = particulas.copy()
  36.  
  37.         for index_central, central in zip(index_old_centrais, old_centrais):
  38.             melhor = 0
  39.             for particula in range(len(central.orbita)):
  40.                 if central.orbita[particula].fitness < central.orbita[melhor].fitness:
  41.                     melhor = particula
  42.        
  43.             if central.orbita[melhor].fitness < central.fitness:
  44.                 # substituo a antiga central com a nova central.
  45.                 nova_pop[index_central] = ParticulaPrincipal(central.orbita[melhor], num_orbitais, alpha_centrais)
  46.                
  47.                 # Transformo as partículas não melhores em orbitais para a nova central.
  48.                 for index, orbital in enumerate(central.orbita):
  49.                     if index != melhor: # não posso transformar a melhor em orbital também!
  50.                         nova_pop[nova_pop.index(central.orbita[index])] = ParticulaOrbital(orbital, nova_pop, alpha_orbitais)
  51.  
  52.                 # Agora a que era central vou transformar em uma orbital. Em nova_pop vou substituir a que era melhor orbital com a antiga central (que agora é orbital)
  53.                 nova_pop[nova_pop.index(central.orbita[melhor])] = ParticulaOrbital(central, nova_pop, alpha_orbitais)
  54.             else:
  55.                 pass # Não é necessário atualizar grupo.
  56.        
  57.         return nova_pop
  58.  
  59. class Particula:
  60.     def __init__(self, dim, xmin, xmax, g_best):
  61.         global alpha_orbitais
  62.         self.xmin = xmin
  63.         self.xmax = xmax
  64.         self.dim  = dim
  65.  
  66.         self.coords   = xmin + np_rand(dim) * (xmax - xmin)
  67.         self.fitness  = rastrigin([self.coords])[0]
  68.  
  69.         self.velocity = 0.5*np_rand(dim)*(xmax - xmin)
  70.         self.max_velo = 0.2*(xmax - xmin)
  71.  
  72.         self.p_best_coords  = self.coords
  73.         self.p_best_fitness = self.fitness
  74.  
  75.         self.coef_velo   = alpha_orbitais[0]
  76.         self.coef_p_best = alpha_orbitais[1]
  77.         self.coef_g_best = alpha_orbitais[2]
  78.  
  79.         self.g_best = g_best
  80.  
  81.         self.g_best.updateGBest(self.coords, self.fitness)
  82.    
  83.     def updateFitness(self):
  84.         self.fitness = rastrigin([self.coords])[0]
  85.         self.g_best.num_fx += 1
  86.  
  87.         if self.fitness < self.p_best_fitness:
  88.             self.p_best_coords  = self.coords.copy()
  89.             self.p_best_fitness = self.fitness.copy()
  90.  
  91.             self.g_best.updateGBest(self.coords, self.fitness)
  92.  
  93.     def updateVelocity(self, update_coef_velo=True):
  94.         "Atualiza a velocidade da partícula."
  95.         if update_coef_velo:
  96.             self.coef_velo = 1 - self.g_best.num_fx / self.g_best.max_fx
  97.  
  98.         self.velocity = self.coef_velo*self.velocity + self.coef_p_best*(self.p_best_coords - self.coords) + self.coef_g_best*(self.g_best.g_best_coords - self.coords)
  99.  
  100.         for d in range(self.dim):
  101.             self.velocity[d] = max(self.velocity[d], -self.max_velo)
  102.             self.velocity[d] = min(self.velocity[d], +self.max_velo)
  103.  
  104.     def move(self):
  105.         "Movimenta a partícula já atualizando sua velocidade."
  106.  
  107.         self.updateVelocity()
  108.  
  109.         self.coords += self.velocity
  110.  
  111.         for d in range(dim):
  112.             self.coords[d] = max(self.coords[d], self.xmin)
  113.             self.coords[d] = min(self.coords[d], self.xmax)
  114.  
  115. class ParticulaPrincipal(Particula):
  116.     def __init__(self, particula, max_orbitais, alpha):
  117.         "Transforma uma particula orbital ou sem classe em uma principal."
  118.  
  119.         self.xmin = particula.xmin
  120.         self.xmax = particula.xmax
  121.         self.dim  = particula.dim
  122.  
  123.         self.coords, self.fitness = particula.coords, particula.fitness
  124.  
  125.         self.velocity = particula.velocity
  126.         self.max_velo = particula.max_velo
  127.  
  128.         self.p_best_coords  = particula.p_best_coords
  129.         self.p_best_fitness = particula.p_best_fitness
  130.  
  131.         self.coef_velo   = alpha[0]
  132.         self.coef_p_best = alpha[1]
  133.         self.coef_g_best = alpha[2]
  134.  
  135.         self.g_best = particula.g_best
  136.  
  137.         self.num_orbitais = 0
  138.         self.max_orbitais = max_orbitais
  139.        
  140.         self.orbita = [0]*self.max_orbitais
  141.    
  142.     def addParticula(self, particula):
  143.         "Assume que já foi verificado se pode adicionar na órbita dessa."
  144.         # print("@", particula)
  145.         self.orbita[self.num_orbitais] = particula
  146.         self.num_orbitais += 1
  147.    
  148.     def removeParticula(self, particula=None, index=None):
  149.         if not particula == None:
  150.             index = self.orbita.index(particula)
  151.         else:
  152.             print("PP:remove:Argumento inválido.")
  153.             return
  154.        
  155.         self.orbita.pop(index)
  156.         self.num_orbitais -= 1
  157.  
  158. class ParticulaOrbital(Particula):
  159.     """Ao transformar uma partícula em orbital, a coordenada da nova partícula orbital será a coordenada da central mais de -2 a 2 vezes a velocidade da partícula.
  160.    Deve receber também a lista de partículas, para assim escolher uma principal para orbitar."""
  161.     def __init__(self, particula, particulas, alpha):
  162.  
  163.         particulas_principais_com_orbitas_vazias = [p for p in particulas if type(p) == ParticulaPrincipal and p.num_orbitais < p.max_orbitais]
  164.         index_p_principal_escolhida = np.random.randint(0, len(particulas_principais_com_orbitas_vazias))
  165.         self.orbiting = particulas_principais_com_orbitas_vazias[index_p_principal_escolhida]
  166.         self.orbiting.addParticula(self)
  167.  
  168.         self.xmin = self.orbiting.xmin
  169.         self.xmax = self.orbiting.xmax
  170.         self.dim  = self.orbiting.dim
  171.  
  172.         self.velocity = particula.velocity
  173.         self.max_velo = particula.max_velo
  174.  
  175.         self.p_best_coords  = particula.p_best_coords
  176.         self.p_best_fitness = particula.p_best_fitness
  177.  
  178.         self.coef_velo   = alpha[0]
  179.         self.coef_p_best = alpha[1]
  180.         self.coef_g_best = alpha[2]
  181.  
  182.         self.g_best = particula.g_best
  183.  
  184.         self.coords = self.orbiting.coords + np.random.randint(-2, 2 + 1)*particula.velocity
  185.         self.updateFitness()
  186.  
  187.         # self.coords, self.fitness = particula.coords, particula.fitness
  188.  
  189.     def updateVelocity(self):
  190.         Particula.updateVelocity(self)
  191.  
  192.         # Vai reescrever updateVelocity de forma que com ela seja possível implementar o movimento de órbita.
  193.         # E se eu mudar apenas isso, script não precisa de grandes modificações, dado que self.move() vai realizar o movimento.
  194.         #self.coef_velo = 1 - self.g_best.num_fx / self.g_best.max_fx
  195.  
  196.         # while abs((sinal := np.random.randint(-1, 1 + 1))*np.sqrt(sum((mov := np.random.rand(self.dim))**2))) >= 1 or sinal == 0:
  197.             # pass
  198.         # else:
  199.             # self.coords = ( self.orbiting.coords
  200.                             # + mov
  201.                             # # + self.coef_velo*mov
  202.                             # #
  203.                             # # + self.coef_p_best*(self.p_best_coords - self.coords)
  204.                             # #
  205.                             # # + self.coef_g_best*(self.g_best.g_best_coords - self.coords)
  206.                         # )
  207.             # # self.coef_velo*self.velocity + self.coef_p_best*(self.p_best_coords - self.coords) + self.coef_g_best*(self.g_best.g_best_coords - self.coords)
  208.  
  209.  
  210.  
  211. def rastrigin(pop):
  212.     return np.array([10*len(x) + sum([(_x**2 - 10 * np.cos(2 * np.pi * _x)) for _x in x]) for x in pop])
  213.  
  214. g_best, pop     = GBest(num_fx := 0, max_fx := 1e3), []
  215. xmin, xmax, dim = -5.12, +5.12, 5
  216.  
  217. global alpha_orbitais, alpha_centrais
  218. alpha_orbitais = np.array([1., 1., 1.])
  219. alpha_centrais = np.array([1., .5, 1.])
  220.  
  221. plot  = True
  222.  
  223. num_pop, qtd_centrais = 20, 5
  224. global num_orbitais
  225. num_orbitais = num_pop//qtd_centrais - 1
  226.  
  227. for _ in range(0, num_pop):
  228.     pop.append(Particula(dim, xmin, xmax, g_best)) # cria todas as partículas
  229.  
  230. for i in range(qtd_centrais):
  231.     index = int(np.random.randint(0, num_pop))
  232.     while type(pop[index]) == ParticulaPrincipal:
  233.         index = int(np.random.randint(0, num_pop))
  234.     else:
  235.         pop[index] = ParticulaPrincipal(pop[index], max_orbitais=num_orbitais, alpha=alpha_centrais) # estou assumindo que sempre é 3.
  236.  
  237. for i in range(num_pop):
  238.     if type(pop[i]) == Particula:
  239.         pop[i] = ParticulaOrbital(pop[i], pop, alpha=alpha_orbitais) # transforma todas as que não são centrais em orbitais.
  240.  
  241. if plot:
  242.     plt.ion()
  243.  
  244. tempo = 0
  245.  
  246. while g_best.num_fx < g_best.max_fx:
  247.     begin = t.perf_counter()
  248.  
  249.     particles = []
  250.     p_best    = []
  251.  
  252.     for i in range(num_pop):
  253.         pop[i].move()
  254.         pop[i].updateFitness()
  255.  
  256.         particles.append(pop[i].coords)
  257.         p_best.append(pop[i].p_best_coords)
  258.    
  259.     pop = g_best.updateGroups(pop)
  260.    
  261.     tempo += t.perf_counter() - begin
  262.    
  263.     if plot:
  264.         plt.clf()
  265.         axes = plt.gca()
  266.         axes.set_xlim([xmin, xmax])
  267.  
  268.         plt.plot([x[0] for x in particles], [x[1] for x in particles], 'bo')
  269.         plt.plot([x[0] for x in p_best],    [x[1] for x in p_best],    'ro')
  270.         plt.plot(g_best.g_best_coords[0],   g_best.g_best_coords[1],   'gx')
  271.  
  272.         plt.title("numfx %d " %g_best.num_fx + str(g_best.g_best_fitness) + " : " + str(round(tempo, 4)) + " segundos.")
  273.    
  274.         plt.pause(1e-323)
  275.         plt.show()
  276. else:
  277.     print("Parou! Tempo gasto pelo algoritmo:", tempo)
  278.  
  279.     print("G BEST:")
  280.     print(g_best.g_best_coords, g_best.g_best_fitness)
  281.  
  282.     for p in pop:
  283.         print(p.coords)
  284.    
  285.     if plot:
  286.         try:
  287.             plt.pause(999999999999)
  288.         except Exception:
  289.             pass
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!