gene.hs

import Control.Monad
import Control.Monad.Random

import Data.Bits
import Data.List

type Nucleotido = Bool

type Individuo = [Nucleotido]

type Poblacion = [Individuo]

type Probabilidad = Double

type TasaMutacion = Double

type ModeloAptitud = Individuo -> Double

type FuncionAjuste = Individuo -> Individuo

seleccion :: ModeloAptitud -> [(Individuo, Probabilidad)] -> [Individuo]
seleccion aptitud pob' = map (muestra . snd) pob'
  where
    pob = map fst pob'
    aptitudes = map aptitud pob
    totalAptitud = sum aptitudes
    prob = map (/ totalAptitud) aptitudes
    probAcum = tail (scanl (+) 0 prob)
    probAcum' = map (1-) probAcum
    elementos = zip probAcum' pob
    muestra u = snd . head . dropWhile (\x -> fst x > u) $ elementos

mutacion :: TasaMutacion -> [(Nucleotido, Probabilidad)] -> Individuo
mutacion pmut crom' = map mutar crom'
  where
    crom = map fst crom'
    mutar (b, u)
      | u < pmut  = not b
      | otherwise = b

cruce :: Probabilidad -> (Individuo, Individuo) -> (Individuo, Individuo)
cruce u (ind1, ind2) = (ind1inf ++ ind2sup, ind2inf ++ ind1sup)
  where
    long = fromIntegral (length ind1)
    puntoCruz = fromIntegral $ floor (u * long)
    (ind1inf, ind1sup) = splitAt puntoCruz ind1
    (ind2inf, ind2sup) = splitAt puntoCruz ind2

data ProblemaGenetico = PG  { pgAptitud :: ModeloAptitud
                            , pgAjuste  :: FuncionAjuste
                            , pgTasaMut :: TasaMutacion
                            , pgTamPob  :: Int
                            , pgNumGen  :: Int
                            , pgTamCrom :: Int
                            }

generacion :: (RandomGen g) => ProblemaGenetico -> Poblacion -> Rand g Poblacion
generacion pg pob = (ajustar >=> seleccionar) pob
  where
    ajustar pob = do
      return $ map (pgAjuste pg) pob
    
    seleccionar pob = do
      psel <- replicateM (pgTamPob pg) getRandom
      return $ seleccion (pgAptitud pg) (zip pob psel)
    
    cruzar pob = return pob
    
    mutar pob = forM pob $ \ind -> do
      pmut <- replicateM (pgTamPob pg) getRandom
      return $ mutacion (pgTasaMut pg) (zip ind pmut)

simulacion :: (RandomGen g) => ProblemaGenetico -> Rand g Poblacion
simulacion pg = do
  inicial <- replicateM (pgTamPob pg) genIndividuo
  foldM paso inicial [1..pgNumGen pg]
  where
    paso pob _ = generacion pg pob
    genIndividuo = replicateM (pgTamCrom pg) getRandom

-------------------------

logBase2 x = finiteBitSize x - 1 - countLeadingZeros x

data Planta = Planta Int Int Double Double Double

plantaLimite (Planta _ m _ _ _) = m

plantaLong (Planta n _ _ _ _) = n

plantaCosto (Planta _ _ a b c) x = a * x * x + b * x + c

problemaPlantas :: ProblemaGenetico
problemaPlantas = PG aptitud ajuste tasaMut tamPob numGen tamCrom

codificar :: [Int] -> Individuo
codificar sols = concat (zipWith cod plantas sols)
  where
    cod planta sol = map (testBit sol) [0..(plantaLong planta - 1)]

decodificar :: Individuo -> [Int]
decodificar ind = map dec partes
  where
    longitudes = map plantaLong plantas
    indices = init $ scanl (+) 0 longitudes
    partes = zipWith (\i n -> take n $ drop i ind) indices longitudes
    dec bs = foldl set zeroBits (zip bs [0..])
      where
        set x (b, i) = if b then x `setBit` i else x

costo :: Individuo -> Double
costo ind = sum $ zipWith plantaCosto plantas (map fromIntegral $ decodificar ind)

aptitud :: ModeloAptitud
aptitud ind = if falla == maxEnergia then 0 else 1 / (costo ind * (falla + 1))
  where
    sols = map fromIntegral $ decodificar ind
    falla = abs (maxEnergia - sum sols)

ajuste :: FuncionAjuste
ajuste ind = codificar sols'
  where
    sols  = decodificar ind
    sols' = zipWith ajustar plantas sols
    ajustar planta sol = sol `rem` (1 + plantaLimite planta)

tasaMut = 0.01
tamPob  = 1000
numGen  = 20
tamCrom = sum (map plantaLong plantas)

maxEnergia = 180

plantas = [ planta 100 0.3    (-18.0)   738.0
          , planta 100 0.166  (-6.66)   616.0
          , planta 100 0.208  (-12.08)  733.3
          , planta 100 0.255  (-15.4)   685.0 ]
  where
    planta m = Planta (1 + logBase2 m) m

main = do
  res <- evalRandIO (simulacion problemaPlantas)
  let mejor = last $ sortOn (pgAptitud problemaPlantas) res
  
  print (decodificar mejor)
  print (costo mejor)

gene.py

from random import random

long_cromosoma = 28
tam_poblacion = 300
num_generaciones = 10
tasa_cruzamiento = 0.85
tasa_mutacion = 0.02

num_plantas = 4

class blank(object):
    pass

def mk_planta(lim_superior, longitud, a, b, c):
    p = blank()
    p.lim_superior = lim_superior
    p.longitud = longitud
    p.a = a
    p.b = b
    p.c = c
    return p

def mk_individuo():
    c = blank()
    c.cromosoma = [False for _ in range(long_cromosoma)]
    c.aptitud = 0
    c.soluciones = [0 for _ in range(num_plantas)]
    c.falla = 0
    c.costo = 0
    return c

pob_anterior = [mk_individuo() for _ in range(tam_poblacion)]
pob_actual = [mk_individuo() for _ in range(tam_poblacion)]
aptitud_pob = 0

plantas = [
    mk_planta(100, 7, 0.3, -18, 738),
    mk_planta(100, 7, 0.166, -6.66, 616),
    mk_planta(100, 7, 0.208, -12.08, 733.3),
    mk_planta(100, 7, 0.255, -15.4, 685)
]
capacidad = 180

def crear_poblacion_inicial():
    global pob_anterior
    for i in range(tam_poblacion):
        for j in range(long_cromosoma):
            pob_anterior[i].cromosoma[j] = random() < 0.5

def decodificar_cromosoma(cromosoma, soluciones):
    posicion = 0
    for i in range(num_plantas):
        acumulado = 0
        potencia = 1
        for j in range(plantas[i].longitud):
            if cromosoma[posicion]:
                acumulado += potencia
            potencia *= 2
            posicion += 1
        soluciones[i] = ajustar_solucion(acumulado, plantas[i].lim_superior)

def ajustar_solucion(valor, limite):
    return valor % (limite + 1)

def evaluar_poblacion(poblacion):
    global aptitud_pob
    aptitud_pob = 0
    for i in range(tam_poblacion):
        decodificar_cromosoma(poblacion[i].cromosoma, poblacion[i].soluciones)
        
        evaluar_individuo(poblacion[i])
        
        #evaluar_cromosoma(poblacion[i].soluciones, poblacion[i].aptitud, poblacion[i].falla, poblacion[i].costo)
        
        aptitud_pob += poblacion[i].aptitud

def evaluar_individuo(individuo):
    falla = capacidad
    costo = 0
    for i in range(num_plantas):
        x = individuo.soluciones[i]
        costo += plantas[i].a * x * x + plantas[i].b * x + plantas[i].c
        falla -= x
    falla = abs(falla)
    if falla == capacidad:
        aptitud = 0
    else:
        aptitud = 1 / (costo * (falla + 1))
    individuo.costo = costo
    individuo.aptitud = aptitud
    individuo.falla

def evaluar_cromosoma(soluciones, aptitud, falla, costo):
    falla = capacidad
    costo = 0
    for i in range(num_plantas):
        x = soluciones[i]
        costo += plantas[i].a * x * x + plantas[i].b * x + plantas[i].c
        falla -= x
    falla = abs(falla)
    if falla == capacidad:
        aptitud = 0
    else:
        aptitud = 1 / (costo * (falla + 1))

def seleccion():
    global pob_anterior, pob_actual, aptitud_pob
    for i in range(tam_poblacion):
        randomico = random() * aptitud_pob
        acumulado = 0
        for j in range(tam_poblacion):
            acumulado += pob_anterior[j].aptitud
            if randomico <= acumulado:
                break
        pob_actual[i].cromosoma = pob_anterior[j].cromosoma

def cruzamiento():
    global pob_actual
    vector_rand = [random() for _ in range(tam_poblacion)]
    for i in range(tam_poblacion):
        if vector_rand[i] >= tasa_cruzamiento:
            continue
        for j in range(i+1, tam_poblacion):
            if vector_rand[j] >= tasa_cruzamiento:
                continue
            puntocruz = int(long_cromosoma * random() + 1)
            if puntocruz > 0 and puntocruz < long_cromosoma:
                for k in range(puntocruz+1, long_cromosoma):
                    aux = pob_actual[i].cromosoma[k]
                    pob_actual[i].cromosoma[k] = pob_actual[j].cromosoma[k]
                    pob_actual[j].cromosoma[k] = aux
        # i = j

def mutacion():
    global pob_actual
    for i in range(tam_poblacion):
        for j in range(long_cromosoma):
            rnd = random()
            if rnd < tasa_mutacion:
                pob_actual[i].cromosoma[j] = not pob_actual[i].cromosoma[j]

def resultados_poblacion(poblacion):
    aptitud = 0
    j = 0
    for i in range(tam_poblacion):
        if poblacion[i].aptitud > aptitud:
            aptitud = poblacion[i].aptitud
            j = i
    val_cromosoma = ''
    for i in range(long_cromosoma):
        val_cromosoma += str(int(poblacion[j].cromosoma[i]))
    val_solucion = ''
    for i in range(num_plantas):
        val_solucion += str(poblacion[j].soluciones[i]) + ' '
    print(val_solucion)
    print(poblacion[j].costo)

def simular():
    global pob_anterior, pob_actual
    generacion = 0
    crear_poblacion_inicial()
    evaluar_poblacion(pob_anterior)
    resultados_poblacion(pob_anterior)
    for generacion in range(num_generaciones):
        seleccion()
        cruzamiento()
        mutacion()
        evaluar_poblacion(pob_actual)
        resultados_poblacion(pob_actual)
        pob_anterior = pob_actual

simular()