226 lines
7.1 KiB
Python
226 lines
7.1 KiB
Python
import numpy as np
|
|
import os
|
|
import pickle
|
|
from scipy.spatial.distance import euclidean
|
|
from itertools import combinations, product
|
|
import matplotlib.pyplot as plt
|
|
import sys
|
|
sys.path.append("/home/ritchie46/code/python/anaStruct")
|
|
from anastruct.fem.system import SystemElements
|
|
|
|
|
|
class DNA:
|
|
def __init__(self, length, height, pop_size=600, cross_rate=0.8, mutation_rate=0.0001):
|
|
self.length = length
|
|
self.height = height
|
|
self.mirror_line = length // 2
|
|
self.pop_size = pop_size
|
|
self.cross_rate = cross_rate
|
|
self.mutation_rate = mutation_rate
|
|
# Assumed that length > height
|
|
# product: permutations with replacement.
|
|
self.loc = np.array(list(filter(lambda x: x[1] <= height, product(range(self.mirror_line + 1), repeat=2))))
|
|
|
|
# Index tuples of possible connections
|
|
# filters all the vector combinations with an euclidean distance < 1.5.
|
|
# dna
|
|
self.comb = np.array(list(filter(lambda x: euclidean(self.loc[x[1]], self.loc[x[0]]) < 1.5,
|
|
combinations(range(len(self.loc)), 2))))
|
|
|
|
self.pop = np.random.randint(0, 2, size=(pop_size, len(self.comb)))
|
|
|
|
self.builds = None
|
|
|
|
def build(self):
|
|
builds = np.zeros(self.pop_size, dtype=object)
|
|
middle_node = np.zeros(self.pop_size, dtype=int)
|
|
all_lengths = np.zeros(self.pop_size, dtype=int)
|
|
n_elements = np.zeros(self.pop_size, dtype=int)
|
|
|
|
for i in range(self.pop.shape[0]):
|
|
ss = SystemElements()
|
|
on = np.argwhere(self.pop[i] == 1)
|
|
|
|
for j in on.flatten():
|
|
n1, n2 = self.comb[j]
|
|
l1 = self.loc[n1]
|
|
l2 = self.loc[n2]
|
|
|
|
ss.add_element([l1, l2])
|
|
# add mirror
|
|
ss.add_element([mirror(l1, self.mirror_line), mirror(l2, self.mirror_line)])
|
|
|
|
# Placing the supports on the outer nodes, and the point load on the middle node.
|
|
x_range = ss.nodes_range('x')
|
|
if len(x_range) <= 2:
|
|
builds[i] = None
|
|
all_lengths[i] = 0
|
|
n_elements[i] = 0
|
|
else:
|
|
length = max(x_range)
|
|
start = min(x_range)
|
|
ids = list(ss.node_map.keys())
|
|
|
|
max_node_id = ids[np.argmax(x_range)]
|
|
|
|
for j in range(self.height):
|
|
middle_node_id = ss.nearest_node("both", np.array([(length + start) / 2, self.height - j]))
|
|
if middle_node_id:
|
|
break
|
|
|
|
if middle_node_id is None:
|
|
middle_node_id = ids[np.argmin(np.abs(np.array(x_range) - (length + start) / 2))]
|
|
|
|
ss.add_support_hinged(1)
|
|
ss.add_support_hinged(max_node_id)
|
|
ss.point_load(middle_node_id, Fz=-100)
|
|
|
|
builds[i] = ss
|
|
middle_node[i] = middle_node_id
|
|
all_lengths[i] = length
|
|
n_elements[i] = on.size
|
|
|
|
self.builds = builds
|
|
return builds, middle_node, all_lengths, n_elements
|
|
|
|
def get_fitness(self):
|
|
builds, middle_node, fitness_l, fitness_n = self.build()
|
|
fitness_w = np.zeros(self.pop_size)
|
|
|
|
for i in range(builds.shape[0]):
|
|
if validate_calc(builds[i]):
|
|
|
|
w = np.abs(builds[i].get_node_displacements(middle_node[i])["uy"])
|
|
|
|
x_range = builds[i].nodes_range('x')
|
|
length = max(x_range) - min(x_range)
|
|
fitness_w[i] = 1.0 / (w / ((100 * length**3) / (48 * builds[i].EI)))
|
|
|
|
fitness_n = (400 / fitness_n)**2
|
|
|
|
return fitness_l**2 + fitness_n + fitness_w, fitness_w, fitness_n
|
|
|
|
def crossover(self, parent, pop, fitness):
|
|
if np.random.rand() < self.cross_rate:
|
|
# i = np.random.choice(np.arange(self.pop_size), size=1, p=fitness / np.sum(fitness))
|
|
i = np.random.randint(0, self.pop_size, size=1)
|
|
cross_index = np.random.randint(0, 2, size=self.comb.shape[0]).astype(np.bool)
|
|
parent[cross_index] = pop[i, cross_index]
|
|
|
|
return parent
|
|
|
|
def mutate(self, child):
|
|
i = np.where(np.random.random(self.comb.shape[0]) < self.mutation_rate)[0]
|
|
child[i] = np.random.randint(0, 2, size=i.shape)
|
|
return child
|
|
|
|
def evolve(self, fitness):
|
|
pop = rank_selection(self.pop, fitness)
|
|
pop_copy = pop.copy()
|
|
|
|
for i in range(pop.shape[0]):
|
|
parent = pop[i]
|
|
child = self.crossover(parent, pop_copy, fitness)
|
|
child = self.mutate(child)
|
|
parent[:] = child
|
|
|
|
self.pop = pop
|
|
|
|
|
|
def rank_selection(pop, fitness):
|
|
order = np.argsort(fitness)[::-1]
|
|
pop = pop[order]
|
|
|
|
rank_p = 1 / np.arange(1, pop.shape[0] + 1)
|
|
idx = np.random.choice(np.arange(pop.shape[0]), size=pop.shape[0], replace=True, p=rank_p / np.sum(rank_p))
|
|
return pop[idx]
|
|
|
|
|
|
def validate_calc(ss):
|
|
try:
|
|
a = ss.validate()
|
|
displacement_matrix = ss.solve()
|
|
return not np.any(np.abs(displacement_matrix) > 1e9) and a
|
|
except (np.linalg.LinAlgError, AttributeError):
|
|
return False
|
|
|
|
|
|
def normalize(x):
|
|
if np.allclose(x, x[0]):
|
|
return np.ones(x.shape)*0.1
|
|
return (x - np.min(x)) / (np.max(x) - np.min(x))
|
|
|
|
|
|
def choose_fit_parent(pop):
|
|
"""
|
|
https://www.electricmonk.nl/log/2011/09/28/evolutionary-algorithm-evolving-hello-world/
|
|
|
|
:param pop: population sorted by fitness
|
|
:return:
|
|
"""
|
|
# product uniform distribution
|
|
i = int(np.random.random() * np.random.random() * (pop.shape[1] - 1))
|
|
return pop[i]
|
|
|
|
|
|
def mirror(v, m_x):
|
|
"""
|
|
|
|
:param v: (array) vertex
|
|
:param m_x: (int) mirror x value
|
|
:return: (array) vertex
|
|
"""
|
|
|
|
return np.array([m_x + m_x - v[0], v[1]])
|
|
|
|
|
|
a = DNA(10, 6, 200, cross_rate=0.8, mutation_rate=0.05)
|
|
# plt.ion()
|
|
|
|
|
|
base_dir = "/home/ritchie46/code/machine_learning/vanilla-machine-learning/genetic_algorithms/img/"
|
|
name = "n3"
|
|
os.makedirs(os.path.join(base_dir, f"best_{name}"), exist_ok=1)
|
|
|
|
with open(os.path.join(base_dir, f"best_{name}", "save.pkl"), "rb") as f:
|
|
a = pickle.load(f)
|
|
# a.mutation_rate = 0.1
|
|
# a.cross_rate= 0.8
|
|
f, w, n = a.get_fitness()
|
|
f[np.argwhere(w == 0)] = 0
|
|
idx = np.argmax(f)
|
|
print(w[idx], n[idx])
|
|
a.builds[idx].show_bending_moment()
|
|
|
|
last_fitness = 0
|
|
|
|
for i in range(100, 150):
|
|
fitness, w, n = a.get_fitness()
|
|
|
|
fitness[np.argwhere(w == 0)] = 0
|
|
|
|
a.evolve(fitness)
|
|
|
|
max_idx = np.argmax(fitness)
|
|
print("gen", i, "max fitness", fitness[max_idx], "w", w[max_idx], "n", n[max_idx])
|
|
|
|
if i % 1 == 0:
|
|
|
|
plt.cla()
|
|
|
|
if last_fitness != fitness[max_idx]:
|
|
try:
|
|
fig = a.builds[max_idx].show_structure(show=False, verbosity=1)
|
|
plt.title(f"fitness = {round(fitness[max_idx], 3)}")
|
|
fig.savefig(os.path.join(base_dir, f"best_{name}", f"ga{i}.png"))
|
|
except AttributeError:
|
|
pass
|
|
|
|
last_fitness = fitness[max_idx]
|
|
# plt.pause(0.5)
|
|
|
|
if i % 1 == 0:
|
|
with open(os.path.join(base_dir, f"best_{name}", "save.pkl"), "wb") as f:
|
|
pickle.dump(a, f)
|
|
|