diff --git a/Geo_Manim.py b/Geo_Manim.py deleted file mode 100644 index 3d8bae6..0000000 --- a/Geo_Manim.py +++ /dev/null @@ -1,84 +0,0 @@ -#Librerías necesarias -import matplotlib.pyplot as plt -import numpy as np -from shapely.geometry import LineString -from sympy.parsing.sympy_parser import parse_expr, standard_transformations, convert_equals_signs -from geometric import geometric_aproach -from utils import * -from manim import * -# from manimlib import * -config.frame_width =50 -config.frame_height =50 -w = config.frame_width/7 -h = config.frame_height/7 -x1 = np.arange(-100, 150, 10) -y1 = np.arange(-100, 150, 10) -p1 = np.arange(-100, 150, 10) - -with open('geometric_aproach.json') as settings: - data = json.load(settings) -function_ = data["func"] -constraints_ = data["constraints"] -(pairs, _lines, m, n, intersections, intersects_evals, Xmax, Ymax, Zmax, Xmin, Ymin, Zmin) = geometric_aproach(constraints_, function_, x1, y1) - - - -dots = [] -lines = [] -intercepts_dots = [] -intercepts = [] -if Xmax is not None: - for pair in pairs: - x_pair, y_pair = pair - # for i in range(len(x_pair)): - try: - x_i=x_pair[0].p/x_pair[0].q - x_e=x_pair[-1].p/x_pair[0].q - except: - x_i = x_pair[0].item() - x_e = x_pair[-1].item() - try: - y_i= y_pair[0].p/y_pair[0].q - y_e= y_pair[-1].p/y_pair[0].q - except: - y_i = y_pair[0].item() - y_e = y_pair[-1].item() - init_dot = Dot([x_i/w, y_i/h, 0]).set_color(ORANGE) - end_dot = Dot([x_e/w, y_e/h, 0]) - # self.add(Dot([6, -10, 0])) - dots.append(init_dot) - dots.append(end_dot) - # dots.append((init_dot, end_dot)) - lines.append(Line([x_i/w, y_i/h, 0], [x_e/w, y_e/h, 0]).set_color(ORANGE)) - - for (x, y) in intersections: - intercepts_dots.append(Dot([x/w, y/h, 0], stroke_width=10)) - intercepts.append([x/w, y/h, 0]) - a = 0 - - max_dot = Dot([Xmax/w, Ymax/h, 0], stroke_width=25).set_color(RED) - max_text = Text(f'Max: {Zmax} at ({Xmax}, {Ymax}) with red color') - min_dot = Dot([Xmin/w, Ymin/h, 0], stroke_width=25).set_color(BLUE) - min_text = Text(f'Min: {Zmin} at ({Xmin}, {Ymin}) with blue color') - - -class Geo_Manim(Scene): - def construct(self): - - for line in lines: - self.play(Create(line)) - for d in intercepts_dots: - self.play(Create(d)) - if Xmax is not None: - pol = Polygon(*intercepts).set_fill(RED,opacity=0.5) - self.play(Create(pol)) - self.play(Create(max_dot)) - self.play(Create(min_dot)) - self.play(Write(max_text)) - self.play(max_text.animate.move_to(DOWN*5)) - self.play(Write(min_text)) - self.play(min_text.animate.move_to(DOWN*3)) - self.wait(3) - else: - self.play(Write(Text("There is no Ansewer"))) - \ No newline at end of file diff --git a/Informe para el usuario.md b/Informe para el usuario.md deleted file mode 100644 index e3ea2f9..0000000 --- a/Informe para el usuario.md +++ /dev/null @@ -1,41 +0,0 @@ -### Informe para el usuario: - -Este sería un ejemplo de como ejecutar el programa, el flag __-ql__ es para la calidad, la __l__ significa low, existen otros como __m__, __h__, __k__ los cuales corresponden a medium, high y 4k respectivamente. Se debe tener cuidado al usar estos por el consumo de CPU. También podemos añadirle __-p__ al flag anterior quedando __-pql__ para que al terminar de ejecutar el programa muestre el resultado en un video. Los demás parametros deben ser el nombre del .py a ejecutar y el nombre de la clase que implementa las clases de Manim. El flag --format=gif es usado para que el archivo que retorne sea en formato gif y no un video. - -![drag-img](/home/regnod/Pictures/Screenshot from 2022-02-18 18-04-38.png) - -Un ejemplo del proceso del programa. - -![drag-img](/home/regnod/Pictures/Screenshot from 2022-02-18 17-37-20.png) - - - -Cuando llegue aqui ya terminó el proceso y muestra el path donde se encuentra el archivo. - -![drag-img](/home/regnod/Pictures/Screenshot from 2022-02-18 17-38-00.png) - - - -Esta sería la configuración del __penalty_settings.json__ la cual consiste en 2 objetos, el primer string corresponde a la funcion objetivo y el segundo es una lista con las restricciones descritas a través de strings. Los otros parametros son más especificos del método de penalización, la cantidad de iteraciones como caso de parada, el coeficiente de penalización, el valor para aumentar este, los rangos de x y las y. - -![drag-img](/home/regnod/Pictures/Screenshot from 2022-02-18 17-47-22.png) - - - -Esta sería la configuración del __geometric_aproach.json__ la cual consiste en 2 objetos, el primer string corresponde a la funcion objetivo y el segundo es una lista con las restricciones descritas a través de strings. - -![drag-img](/home/regnod/Pictures/Screenshot from 2022-02-18 17-47-32.png) - -A continuación tenemos ejemplos de lo realizado. - -![drag-gif](/home/regnod/Pictures/Screenshot from 2022-02-18 18-15-49.png) - -![drag-gif](/home/regnod/Pictures/Screenshot from 2022-02-18 18-15-17.png) - - - -Aquí podemos observar como sería el resultado de graficar una función con esta erramienta en 3d. - -![drag-gif](/home/regnod/Documents/4to/modelos de optimización 2/media/videos/Penalty_Manim/720p30/Penalty_ManimCE.gif) - -![drag-gif](/home/regnod/Pictures/Screenshot from 2022-02-18 18-15-38.png) \ No newline at end of file diff --git a/Penalty_Manim.py b/Penalty_Manim.py deleted file mode 100644 index 79a24ad..0000000 --- a/Penalty_Manim.py +++ /dev/null @@ -1,89 +0,0 @@ -from manim import * -from utils import * -import numpy as np -from penalty_newton import penalty_newton - -with open('penalty_settings.json') as settings: - data = json.load(settings) - -_dict = penalty_newton() -variables, function, constraints = read_json('penalty_settings.json') -lam = sym.Lambda(convert_list_to_tuples(variables), function) -restrictions_lambdas = [sym.Lambda(convert_list_to_tuples(variables), con) for con in constraints] -path_dots = read_dots_from_json(_dict["points"]) -p_dots = [] -for (x_, y_) in path_dots: - p_dots.append([x_, y_, lam(x_, y_)]) -# p_dots.append(Dot([*path_dots[-1], lam(*path_dots[-1])], stroke_width=35).set_color(RED_D)) - - -a = 0 -# x_ = np.arange(-10,10,1) - -# y_ = np.arange(-10,10,1) -# for i in x_: -# print(str([i, i, lam(i, i)])) -# config.frame_width =50 -# config.frame_height =50 -# w = config.frame_width/7 -# h = config.frame_height/7 -x1 = np.arange(-100, 150, 10) -y1 = np.arange(-100, 150, 10) - -try: - x_range = data['Penalty_x_range'] - y_range = data['Penalty_y_range'] -except: - raise Exception("You must specify the range of coordinates, i.e x_range = [-4,4], y_range = [0,10]") -class Penalty(ThreeDScene): - def construct(self): - axes = ThreeDAxes(tips=False) - def shape(u,v): - z = lam(u, v) - return np.array([u, v, z]) - graph = Surface( - shape, v_range=[-2,2], u_range=[-2,2], - checkerboard_colors=[RED_D, RED_E], resolution=(50, 50), fill_opacity=0.5 - ) - - self.set_camera_orientation(phi=75 * DEGREES, theta=-30 * DEGREES) - # self.renderer.camera.light_source.move_to(3*IN) # changes the source of the light - - # self.play(Rotate(graph, 0*DEGREES)) - # self.play(Rotate(axes, 0*DEGREES)) - - itt_points = [Dot3D(axes.coords_to_point(*i), color=BLUE, radius=0.07) for i in p_dots[1:-1]] # the middle steps points - ppoints = [Dot3D(axes.coords_to_point(*p_dots[0]), color=RED)] - ppoints.extend(itt_points) - ppoints.append(Dot3D(axes.coords_to_point(*p_dots[-1]), color=GREEN)) - - # self.add(graph, axes) - self.play(Create(axes)) - self.play(Create(graph)) - method_text = Text( - f"The point reached was at \n({p_dots[-1][0]}, {p_dots[-1][1]}) \nwith a value of {p_dots[-1][2]}", font_size=20, - - ) - self.add_fixed_in_frame_mobjects(method_text) - method_text.to_corner(UL) - # self.add(method_text) - - # self.play(Write(method_text.to_corner(UL))) - - for p in ppoints: - # self.add(p) - self.play(Create(p)) - self.begin_ambient_camera_rotation(rate=0.3) - self.wait() - self.stop_ambient_camera_rotation() - - self.wait() - # self.add(axes, graph) - # self.begin_ambient_camera_rotation(rate=0.7) - # self.wait(4) - # self.stop_ambient_camera_rotation() - # self.move_camera(phi=75 * DEGREES, theta=30 * DEGREES) - # self.wait() - - - \ No newline at end of file diff --git a/Richard.md b/Richard.md deleted file mode 100644 index 6b70f7c..0000000 --- a/Richard.md +++ /dev/null @@ -1,81 +0,0 @@ -### Solución Geométrica para un problema de optimización lineal: - -La solución geométrica asociada a un problema de optimización lineal se realizó de la siguiente forma: - -Paso 1: Se usan todas las restricciones del problema y se grafican en el plano como rectas. - -Paso 2: Se hallan las intercepciones entre estas rectas, las cuales constituyen puntos. - -Paso 3: Se elije entre todos los puntos de intercepción solo aquellos que cumplan todas las restricciones del problema, los demás se descartan. - -Paso 4: Luego se usa la función objetivo para determinar cuáles son las evaluaciones de estos puntos en la misma y ordenarlos en cuanto a su valor para encontrar el mínimo y el máximo. - -Paso 5: Los resultados de la evaluación se almacenan en una estructura diccionario donde el identificador es un entero entre 0 y la cantidad de puntos de intersección encontrados, después de haber sido ordenados. - -Paso 6: Se crean dos listas __m__, __n__ que contienen las los valores de __x__, __y__ respectivamente de los puntos ordenados. - -Paso 7: Para encontrar el mínimo se escoge de la estructura diccionario anteriormente descrita el identificador que pertenece al menor valor almacenado en este y se utiliza como posición para indexar en las listas __m__, __n__ y de esta forma obtenemos el mínimo. De igual forma obtenemos el máximo buscando el identificador del mayor valor en el diccionario. - -Luego para graficar en Manim, lo que hacemos es ir siguiendo los pasos hasta obtener todos los datos necesarios tales como las rectas, las intersecciones de estas rectas que son válidas para el problema así como el punto máximo y mínimo. Hacemos aparecer las líneas, luego los puntos intersección, luego el poligono asociado a estos puntos y por último los puntos de máximo y mínimo en caso de existir. - -##### Detalles del código: - -Todo el proceso para hayar las intersecciones así como mínimo y máximo están en el archivo __geometric.py__ en un método llamado __geometric_aproach__, el cual recibe 4 parámetros, las restricciones en forma de listas de strings, la ecuación objetivo en forma de string y dos listas de valores __x__, __y__ con los cuales se formaran las rectas. - -Luego el proceso para graficar com Manim se encuentra en el archivo __Geo_Manim.py__ donde primero se realiza el llamado a __geometric_aproach__ para obtener todos los datos necesarios y luego se procede a graficar en orden las líneas, puntos de intersección, el polígono asociado a estos puntos y por último los puntos de máximo y mínimo. - -Por último, los detalles como la función y las restricciones se almacenan en un archivo __geometric_aproach.json__ el cual se utiliza como configuración del problema a resolver. - -Se utilizó como nombre de variables "x" y "y". - - - -### Penalty: - -Método de Penalización: - -Este método reemplaza la función objetivo con restricciones por una función sin restricciones, la cual es formada por: - --La función objetivo - --Un término adicional por cada constraint, que es positivo si el punto actual __x__ viola las restricciones o 0 de otra forma. - -La mayoría de las formas de este método usan un coeficinetepositivo el cual se utiliza para multiplicar la parte de la funcion de las restricciones. Haciendo que este coeficiente sea más grande se penaliza más fuerte cada ves la violación de las restricciones para llegar más cerca a la región factible para el problema restringido. - -Estos son considerados "exterior penalty method", dado que el término de penalización en la función para cada restricción no es 0 cuando x es no factible con respecto a esa restricción. - -El método cuadrático de penalización es el que cada término de penalización es elevado al cuadrado. Además las restricciones de igualdad y desigualdad tienen distintas formas: - --igualdad: h(x) => $(h(x))^2$ - --desigualdad: g(x) => $max(0, g(x))^2$ - -$\Phi(x)$ - -##### Detalles del código: - -Todo el proceso para hayar las intersecciones así como mínimo y máximo están en el archivo __penalty_newton.py__ en un método llamado __penalty_newton__. - -El método se puede describir en pseudo código de la siguiente forma: - -dados el coeficiente de penalización ($\mu_0$), una tolerancia ($\tau_0>0$) y un punto inicial ($x_0$) - -for $k$ = 0, 1,2, . . . - -​ encontrar un minimizador aproximado x_k para la función y termina cuando $||\nabla \Phi(x_k)||<= \tau_k$ - -​ if se satisface el test de convergencia final parar con la aproximacion $x_k$. - -​ else aumentar el coeficiente de penalización y coger el punto hayado como nuevo start_point, repetir el proceso. - -end ( for) - -En el proceso se va modificando la función de penalización dado que el coeficiente por el que se multiplican los términos de penalización varía. - -Al final se devuelve un diccionario con la cantidad de iteraciones realizadas, los puntos intermedios y el una lista con los valores de __x__, __y__ y el resultado de evaluar la funcion en estos (__z__). - -Luego el proceso para graficar com Manim se encuentra en el archivo __Penalty_Manim.py__ donde primero se realiza el llamado a __penalty_newton__ para obtener todos los datos necesarios y luego se procede a graficar en orden el eje de coordenadas, la función y los puntos resultantes. - -Por último, los detalles como la función y las restricciones se almacenan en un archivo __penalty_settings.json__ el cual se utiliza como configuración del problema a resolver. - -Se utilizó como nombre de variables "x" y "y". \ No newline at end of file diff --git a/Richard1.pdf b/Richard1.pdf deleted file mode 100644 index cd94d91..0000000 Binary files a/Richard1.pdf and /dev/null differ diff --git a/asd.py b/asd.py deleted file mode 100644 index c0ec3ac..0000000 --- a/asd.py +++ /dev/null @@ -1,189 +0,0 @@ - - -# from scipy.optimize import minimize, line_search -# import sympy as sym -# from sympy.parsing.sympy_parser import parse_expr, standard_transformations, convert_equals_signs -# from utils import func_eval, get_lambdas, get_constraints_cleared, get_penalty_func, get_penalty_ineq_constraints_minimize, get_eq, get_penalty_eq_constraints -# import numpy as np -# from gradient import gradient -# # rz = lambda x: (1-x[0])**2 + 100*(x[1] - x[0]**2)**2 -# # h_1 = lambda x: (x[0] - 2 * x[1] + 2) -# # h_2 = lambda x: (-x[0] - 2 * x[1] + 6) -# # h_3 = lambda x: (-x[0] + 2 * x[1] + 2) - -# # x0 = [2.3, 5] -# # cons = ({'type': 'ineq', 'fun': h_1}, -# # {'type': 'ineq', 'fun': h_2}, -# # {'type': 'ineq', 'fun': h_3}) -# # minimize(rz, x0, constraints=cons) - -# # x = sym.Symbol('x') -# # y = sym.Symbol('y') - -# ineqs_ = [] -# eqs_ = [] -# # ineqs_.append('20*x+50*y <= 3000') -# # ineqs_.append('x - 2 * y + 2 >= 0') -# # ineqs_.append('-x - 2 * y + 6 >= 0') -# # ineqs_.append('-x + 2 * y + 2 >= 0') -# # # eqs_.append('2 = x + y')#, transformations = standard_transformations+(convert_equals_signs,))) -# # # eqs_.append('2 = x + y')#, transformations = standard_transformations+(convert_equals_signs,))) -# # func = '(1-x)**2+100*(y-x**2)**2' -# ineqs_.append('(x-5)**2+y**2-26>=0') -# func = '(x**2+y-11)**2+(x+y**2-7)**2' -# # a = parse_expr(func) -# # x = ['x'] -# # la = lambda x: func_eval(['x', 'y'], [x[0], x[1]], a).p/func_eval(['x', 'y'], [x[0], x[1]], a).q -# # b = sym.Lambda(x, a) -# # c = la([2, 3]) -# # d = 0 -# # la funcion get_penalty_func recibe en strings: -# # la funcion objetivo, las restricciones, una tupla con las variables, -# # el parametro de penalización y 0 o 1 para minimizar y maximizar respectivamente -# # Devuelve: -# # las inecuaciones -# # la funcion en formato sympy, la funcion de penalizacion en formato de sympy, -# # la funcion de penalización en formato lambda, las inecuaciones y las ecuaciones - - -# def penalty(x_, iters, epsilon_, mm = 0): -# global func, ineqs_, eqs_ -# x_c = initial_point = [0.11, 0.1] - -# c = 1 -# while c < 1000: -# ineqs, func_, penalty_func, penalty_lambda, eqs, parsed_ineqs = get_penalty_func(func, ineqs_+eqs_, x_, c, mm)#('x', 'y'), c, mm) -# fp = lambda x: penalty_lambda(x[0], x[1])#float(func_eval(x_, [x[0], x[1]], penalty_func))#.p/func_eval(x_, [x[0], x[1]], penalty_func).q - -# m = gradient(x_, str(penalty_func), x_c) -# min_ = m["min"] -# min_values = m["min_value"] -# # x_1 = minimize(fp, x_c).x -# x_1 = min_ -# print(x_1) -# track = np.linalg.norm(x_1)/np.linalg.norm(x_c) -# x_c = [x_1[0], x_1[1]] -# if track <= 0.001: -# break -# c*= 2 -# return x_c - -# x_c_ = penalty(['x', 'y'],0,0,0) - -# #3.00000001176746 -# #1.9999999751723971 -# #0.83187365686818, 2.932764818735715 -# a = 0 -# # ineqs = get_penalty_eq_constraints(eqs) -# # converted = get_constraints_cleared(ineqs) -# # eq = get_eq(converted) -# # total_ineq = eq[0] -# # for e in eq[1:]: -# # total_ineq += e -# # total_ineq += 0 -# # x_constraints = [] -# # y_constraints = [] -# # for c in converted: -# # l = c.lhs -# # r = c.rhs -# # op = c.rel_op -# # if len (l.free_symbols) == 1: -# # if l.free_symbols.__contains__(x): -# # x_constraints.append([c, op]) -# # else: -# # y_constraints.append([c, op]) -# # elif len (r.free_symbols) == 1: -# # op = '<=' if op.__contains__('>') else '>=' -# # if r.free_symbols.__contains__(x): -# # x_constraints.append([c, op]) -# # else: -# # y_constraints.append([c, op]) - -# # lambdas_x = get_lambdas([i[0] for i in x_constraints]) -# # lambdas_y = get_lambdas([i[0] for i in y_constraints]) - -# # hasta aqui todo en talla - - -# # x_c = [2.3, 3] -# # i = 1 -# # while i < 1000: -# # curr_func = lambda x: rz(x) + i*(h_1(x)**2 + h_2(x)**2 + h_3(x)**2) -# # x_c = minimize(curr_func, x_c).x -# # i *= 1.2 -# # print(x_c.x) - -# # def eval_eq_constraints(eq, c, x_, vars_): -# # return c * (func_eval(vars_, x_, eq))**2 - -# # def eval_ineq_constraints(ineq, c, x_, vars_): -# # return c * (max(0, func_eval(vars_, x_, ineq)))**2 - -# # def eval_(equation, eq_cons, ineq_cons, x_, vars_, c): -# # valor = 0 -# # for con in eq_cons: -# # valor+= eval_eq_constraints(con, c, x_, vars_) - -# # for con in ineq_cons: -# # valor+= eval_ineq_constraints(con, c, x_, vars_) - -# # valor+= func_eval(vars_, x_, equation) -# # return valor - - -# if __name__ == "__main__": -# func_key_name = get_key_name(function, constraints) -# # x_ip = np.array([0.11, 0.1]) -# x_ip = np.array([2.3, 3]) -# M = 2 ## Specifies the total dimension we are working with.. //Global Var. -# c = 1.55 ##factor for updating r. -# numseq = 20 ## number of sequences for penalty function method. -# r = 0.1 -# grad = np.zeros(M) -# sol = np.zeros(M) -# sol_ = [] -# x_2 = np.zeros(M) -# ## has the principal values in the key name. -# func_key_name += ' ('+str(x_ip)[1:len(str(x_ip))-1]+') '+str(numseq)+' '+str(r)+' '+str(c) -# dic = load_saved_data(func_key_name) -# if dic is None: - -# ## BRACKET OPERATOR PENALTY METHOD..1.1,1.1 -# for k in range(numseq): -# print('\n') -# print('sequence number ', k) -# print('current function value is ', multi_f(x_ip)) -# x_1 = np.copy(x_ip) -# m = newton(_variables, str(penalty_func), x_ip) -# sol = m["min"] -# min_values = m["min_value"] -# # sol = DFP(x_ip) -# x_ip = np.copy(x_1) - -# for j in range(M): -# x_2[j] = sol[j] - x_ip[j] -# track = np.linalg.norm(x_2)/np.linalg.norm(x_ip) -# sol_ = [] -# for i in range(M): -# print(sol[i]) -# sol_.append(sol[i]) -# x_ip[i] = sol[i] -# sol_points.append(sol_) -# r*= c -# if track <= 0.001: -# break -# sol_.append(func_eval(variables, sol_, function)) -# dic = { -# func_key_name: { -# 'iterations': k, -# 'points': str(sol_points), -# 'min': str(sol_), -# } -# } -# save_process(dic) -# else: -# dic = load_saved_data(func_key_name) -# print('\n') -# print('sequence number ', dic['iterations']) -# print('the sequence of points was ', '\n['.join(dic['points'][1:len(dic['points'])-1].split('['))+'\n') -# print('the final solution was ', dic['min']) \ No newline at end of file diff --git a/data.json b/data.json deleted file mode 100644 index 2269c29..0000000 --- a/data.json +++ /dev/null @@ -1 +0,0 @@ -{"(x**2+y-11)**2+(x+y**2-7)**2 (x-5)**2+y**2-26>=0 (2.3, 3) 20 0.1 0.155": {"iterations": 4, "points": "[[2.62787032893929, 2.4747184715759922], [2.9431385041327696, 2.093940561116741], [2.9911746917593227, 2.0152361277836834], [2.9986316182431034, 2.002378938130111], [2.9997878896020325, 2.0003691566158373, 2.41558692916641e-6]]", "min": "[2.9997878896020325, 2.0003691566158373, 2.41558692916641e-6]"}, "x**2+y**2 x>=2 y<=10 (2.3, 3) 20 0.1 0.155": {"iterations": 10, "points": "[[0.18181818159756785, -1.3082008362719662e-10], [0.030526834239981454, -1.9967616362109277e-12], [0.004793483665491628, -4.785722905642079e-15], [0.0007444977473717645, -1.7814670272994916e-18], [0.00011543346197413372, -1.0281988922348497e-22], [1.7893059360695937e-05, -9.199527759859884e-28], [2.773445117832538e-06, -1.2733288043608067e-33], [4.29884499419984e-07, -2.7170138035233225e-40], [6.66321095741996e-08, -8.673576338753747e-48], [1.0327977240182738e-08, -2.7383396660629417e-56], [1.0327977240182738e-08, -2.7383396660629417e-56, 1.06667113873733e-16]]", "min": "[1.0327977240182738e-08, -2.7383396660629417e-56, 1.06667113873733e-16]"}, "(x**2 + y - 11)+(x + y**2 - 7)**2 x>=0 y>=0 (x-5)**2 + y**2 - 26 >= 0 (2.3, 3) 20 0.1 0.155": {"iterations": 6, "points": "[[0.5141717166435834, 2.508878164932901], [0.38394013087724993, 2.54604942945742], [0.19409650478351784, 2.588241173510102], [0.1147650408041444, 2.6052554988113568], [0.09886771536559803, 2.608637458675471], [0.09628833672392983, 2.6091853922232375], [0.09588558158951764, 2.609270930384672, -8.37235362331541]]", "min": "[0.09588558158951764, 2.609270930384672, -8.37235362331541]"}, "(x**2 + y - 11)+(x + y**2 - 7)**2 x>=0 y>=0 (x-5)**2 + y**2 - 26 >= 0 (0.11, 0.1) 20 0.1 0.155": {"iterations": 6, "points": "[[0.5141717155708564, 2.5088781652906116], [0.38394012987086157, 2.5460494294077347], [0.19409650457144298, 2.588241178554411], [0.11476504331504782, 2.605255498044831], [0.09886773223190584, 2.6086374555129255], [0.09628830754776765, 2.609185398013029], [0.09588558397051618, 2.609270928921136, -8.37235362331498]]", "min": "[0.09588558397051618, 2.609270928921136, -8.37235362331498]"}, "(x**2 + y - 11)+(x + y**2 - 7)**2 x>=0 y>=0 (x-5)**2 + y**2 - 26 >= 0 (0.11, 0.11) 20 0.1 0.155": {"iterations": 1, "points": "[[0.09581163371728686, 2.609286627357837], [0.09581163338821046, 2.60928662833568, -8.37235362880134]]", "min": "[0.09581163338821046, 2.60928662833568, -8.37235362880134]"}} \ No newline at end of file diff --git a/geometric.py b/geometric.py deleted file mode 100644 index ffa1984..0000000 --- a/geometric.py +++ /dev/null @@ -1,114 +0,0 @@ -import numpy as np -from shapely.geometry import LineString -from sympy.parsing.sympy_parser import parse_expr -import sympy as sym -from utils import * - -def geometric_aproach(ineqs, equation, x1, y1): - equation = parse_expr(equation) - equation = sym.Lambda(('x', 'y'), equation) - ineqs_, eqs_ = get_geometric_ineqs_and_eqs(ineqs) - converted = get_constraints_cleared(ineqs_) - eq_converted = [get_eq_cleared_constraint(eq_)for eq_ in eqs_] - x_constraints = [] - y_constraints = [] - for c in converted: - l = c.lhs - r = c.rhs - op = c.rel_op - if len (l.free_symbols) == 1: - if l.free_symbols.__contains__(x): - x_constraints.append([c, op]) - else: - y_constraints.append([c, op]) - elif len (r.free_symbols) == 1: - op = '<=' if op.__contains__('>') else '>=' - if r.free_symbols.__contains__(x): - x_constraints.append([c, op]) - else: - y_constraints.append([c, op]) - - lambdas_x = get_lambdas([i[0] for i in x_constraints]) - lambdas_y = get_lambdas([i[0] for i in y_constraints]) - lambdas_x_eq = [] - lambdas_y_eq = [] - for (eq_, symbol) in eq_converted: - if symbol == sym.Symbol('x'): - lambdas_x_eq.append(sym.Lambda(eq_.free_symbols, eq_)) - else: - lambdas_y_eq.append(sym.Lambda(eq_.free_symbols, eq_)) - - pairs = [] - for lam in lambdas_y: - pairs.append((x1, [lam(i) for i in x1])) - - for lam in lambdas_x: - pairs.append(([lam(i)for i in y1], y1)) - - for lam in lambdas_y_eq: - pairs.append((x1, [lam(i) for i in x1])) - - for lam in lambdas_x_eq: - pairs.append(([lam(i)for i in y1], y1)) - - - lines = [] - for pair in pairs: - lines.append(LineString(np.column_stack(pair))) - - # aqui se plotearian las líneas en manim - # for l, (p_x, p_y) in zip(lines, pairs): - # plt.plot(p_x, p_y, '-', linewidth=2, color='k') - - - intersections = [] - for i in range(len(lines)): - j = i+1 - while j < len(lines): - (i_x_, i_y_) = (lines[i].intersection(lines[j])).xy - if len(i_x_) != 0 and check_constraint_point(i_x_[0],i_y_[0], y_constraints, lambdas_y) and check_constraint_point(i_y_[0], i_x_[0], x_constraints, lambdas_x): - if check_eq_constraint_point(i_x_[0],i_y_[0], lambdas_y_eq) and check_eq_constraint_point(i_y_[0],i_x_[0], lambdas_x_eq): - intersections.append((i_x_[0],i_y_[0])) - j+=1 - intersections.sort() - - # aqui se pintarian los puntos intercepción - # for intersection in intersections: - # plt.plot(*intersection, 'o', color='r') - - intersects_evals = [] - for intersection in intersections: - (int_x, int_y) = intersection - intersects_evals.append(round(equation(int_x, int_y), 2)) - intersects_evals.sort() - - dict1 = {} - for i in range(len(intersects_evals)): - dict1[i] = intersects_evals[i] - if len(intersections) >= 1: - #varia según min o max - # Z = intersects_evals[-1] if min_max == 1 else intersects_evals[0] - Zmax = intersects_evals[-1] - Zmin = intersects_evals[0] - - m = [xi[0] for xi in intersections] - n = [yi[1] for yi in intersections] - - # rellenar el polígono - - # posicion = max(dict1, key=dict1.get) if min_max == 1 else min(dict1, key=dict1.get) - posicionMax = max(dict1, key=dict1.get) - posicionMin = min(dict1, key=dict1.get) - - Xmax = m[posicionMax] - Xmin = m[posicionMin] - Ymax = n[posicionMax] - Ymin = n[posicionMin] - - intersects_evals = [round(i.num, 2) for i in intersects_evals] - Zmax = round(Zmax.num, 2) - Zmin = round(Zmin.num, 2) - - return (pairs, lines, m, n, intersections, intersects_evals, Xmax, Ymax, Zmax, Xmin, Ymin, Zmin) - else: - return (pairs, lines, [], [], intersections, [], None, None, None, None, None, None) \ No newline at end of file diff --git a/geometric_aproach.json b/geometric_aproach.json deleted file mode 100644 index 074ef5a..0000000 --- a/geometric_aproach.json +++ /dev/null @@ -1,4 +0,0 @@ -{ - "func": "10000*x + 6000*y", - "constraints": ["20*x+50*y <= 3000", "x+y <= 90", "y >= 10", "y >= 0", "x >= 0"] -} \ No newline at end of file diff --git a/geometric_matplotlib.py b/geometric_matplotlib.py deleted file mode 100644 index 5dc8bfc..0000000 --- a/geometric_matplotlib.py +++ /dev/null @@ -1,114 +0,0 @@ -#Librerías necesarias -import matplotlib.pyplot as plt -import numpy as np -from shapely.geometry import LineString - -#Ecuaciones e intervalos (Para tabular) -x = np.arange(-100, 150, 1) -y = np.arange(-100, 150, 1) -y1 = 3 + (0 * x) -x1 = 20 + (0 * y) -y2 = 20 - (2 * x) -y3 = (7 - (3 * x))/ -2 -z = (-50 * x) / 20 - -#Identificadores para las líneas -primera_linea = LineString(np.column_stack((x, y1))) -segunda_linea = LineString(np.column_stack((x, y2))) -tercera_linea = LineString(np.column_stack((x, y3))) -cuarta_linea = LineString(np.column_stack((x1, y))) -quinta_linea = LineString(np.column_stack((x, z))) - -#Graficando las líneas -plt.plot(x, y1, '-', linewidth=2, color='b') -plt.plot(x, y2, '-', linewidth=2, color='g') -plt.plot(x, y3, '-', linewidth=2, color='r') -plt.plot(x1, y, '-', linewidth=2, color='y') -plt.plot(x, z, ':', linewidth=1, color='k') - -#Generando las intersecciones (vértices) -primera_interseccion = cuarta_linea.intersection(primera_linea) -segunda_interseccion = primera_linea.intersection(segunda_linea) -tercera_interseccion = segunda_linea.intersection(tercera_linea) -cuarta_interseccion = tercera_linea.intersection(cuarta_linea) - -#Graficando los vértices -plt.plot(*primera_interseccion.xy, 'o') -plt.plot(*segunda_interseccion.xy, 'o') -plt.plot(*tercera_interseccion.xy, 'o') -plt.plot(*cuarta_interseccion.xy, 'o') - -#Imprimiendo las coordenadas de los vértices en la consola -print('\n COORDENADAS DE LAS INTERSECCIONES') -print('Coordenadas de la primera intersección: {} '.format(primera_interseccion)) -print('Coordenadas de la segunda intersección: {} '.format(segunda_interseccion)) -print('Coordenadas de la tercera intersección: {} '.format(tercera_interseccion)) -print('Coordenadas de la cuarta intersección: {} '.format(cuarta_interseccion)) - -#Identificando los valores de las coordenadas x y y de cada vértice -xi1m, yi1m = primera_interseccion.xy -xi2m, yi2m = segunda_interseccion.xy -xi3m, yi3m = tercera_interseccion.xy -xi4m, yi4m = cuarta_interseccion.xy - -#Cambiamos el formato de matriz a float -xi1 = np.float64(np.array(xi1m)) -xi2 = np.float64(np.array(xi2m)) -xi3 = np.float64(np.array(xi3m)) -xi4 = np.float64(np.array(xi4m)) -yi1 = np.float64(np.array(yi1m)) -yi2 = np.float64(np.array(yi2m)) -yi3 = np.float64(np.array(yi3m)) -yi4 = np.float64(np.array(yi4m)) - -#Evaluando la función objetivo en cada vértice -FOi1 = (xi1 * 50) + (yi1 * 20) -FOi2 = (xi2 * 50) + (yi2 * 20) -FOi3 = (xi3 * 50) + (yi3 * 20) -FOi4 = (xi4 * 50) + (yi4 * 20) - -#Imprimiendo las evaluaciones de la FO en cada vértice (Consola) -print('\n EVALUACIÓN DE LA FO EN LOS VÉRTICES') -print('Función objetivo en la intersección 1: {} pesos'.format(FOi1)) -print('Función objetivo en la intersección 2: {} pesos'.format(FOi2)) -print('Función objetivo en la intersección 3: {} pesos'.format(FOi3)) -print('Función objetivo en la intersección 4: {} pesos'.format(FOi4)) - -#Calculando el mejor resultado (Maximizar) -ZMAX = max(FOi1, FOi2, FOi3, FOi4) - -#Imprimiendo la solución óptima en la consola -print('\n SOLUCIÓN ÓPTIMA') -print('Solución óptima: {} pesos'.format(ZMAX)) - -#Ordenando las coordenadas de los vértices (Las coordenadas x en m y las coordenadas y en n) -m = [xi1, xi2, xi3, xi4] -n = [yi1, yi2, yi3, yi4] - -#Graficando el polígono solución a partir de las coordenadas de los vértices -plt.fill(m, n, color='silver') - -#Identificando el índice del vértice de la mejor solución -dict1 = {0:FOi1, 1:FOi2, 2:FOi3, 3:FOi4} -posicion = max(dict1, key=dict1.get) - -#Obteniendo las coordenadas del vértice de la mejor solución de acuerdo al índice -XMAX = m[posicion] -YMAX = n[posicion] - -#Imprimiendo las coordenadas del vértice de la mejor solución (variables de decisión) -print('\n VARIABLES DE DECISIÓN') -print('Cantidad de asientos a reservar para fumadores: {} '.format(XMAX)) -print('Cantidad de asientos a reservar para no fumadores: {} '.format(YMAX)) - -#Generando las anotaciones de las coordenadas y solución óptima en el gráfico -plt.annotate(' X: {0} / Y: {1} (Coordenadas)'.format(XMAX, YMAX), (XMAX, YMAX)) -plt.annotate(' Solución óptima: {}'.format(ZMAX), (XMAX, YMAX+3)) - - -#Configuraciones adicionales del gráfico -plt.grid() -plt.xlabel('Asientos para fumadores') -plt.ylabel('Asientos para no fumadores') - -plt.show() diff --git a/gradient.py b/gradient.py deleted file mode 100644 index fc67ba2..0000000 --- a/gradient.py +++ /dev/null @@ -1,103 +0,0 @@ -from functools import reduce -import sympy as sp -from scipy.optimize import line_search -import numpy as np -from pprint import pprint -import json - - -EPSILON = 0.1 - - -def convert_list_to_tuples(list): - return reduce(lambda a, b: a + (b,), list, ()) - - -def gradient(vars, func, initial_point, cycles=100, verbose=False): - vars = convert_list_to_tuples([sp.symbols(v) for v in vars]) - - func = sp.parse_expr(func) - func_lambda = sp.Lambda(vars, func) - func_evaluated = lambda x: np.array([func_lambda(*x)], dtype=float) - - gradient_lambda = np.array([sp.Lambda(vars, sp.diff(func, var)) for var in vars]) - gradient_evaluated = lambda x: np.array( - [g(*x) for g in gradient_lambda], dtype=float - ) - - initial_point = np.array(initial_point, dtype=float) - - points = [ - { - "point": initial_point, - "value": func_evaluated(initial_point), - "iteration": 1, - "gradient": gradient_evaluated(initial_point), - "step_arrived": 0, - } - ] - - if verbose: - pprint(points[-1]) - - # start_point = points[-1]["point"] - for _ in range(cycles): - start_point = points[-1]["point"] - step = points[-1]["iteration"] - gradient = gradient_evaluated(start_point) - search_gradient = -1 * gradient - - try: - res = line_search( - func_evaluated, gradient_evaluated, start_point, search_gradient - ) - print(res[0]) - r_point = start_point + res[0] * search_gradient - # search_gradient = res[5] - # x_2 = [0 for _ in range(len(vars))] - # for j in range(len(vars)): - # x_2[j] = r_point[j] - start_point[j] - - # if((np.linalg.norm(gradient_evaluated(r_point)) <= 0.001) or (np.linalg.norm(x_2)/np.linalg.norm(start_point)) <= 0.001): - # break - - if points[-1]["value"] - func_evaluated(r_point) < EPSILON: - break - - points.append( - { - "point": r_point, - "value": func_evaluated(r_point), - "iteration": step + 1, - "gradient": gradient_evaluated(r_point), - "step_arrived": res[0], - } - ) - - if verbose: - pprint(points[-1]) - - except Exception as e: - break - - return { - "min": points[-1]["point"], - "min_value": points[-1]["value"], - "points": points, - } - - -def main(): - f = open("gradient.json") - data = json.load(f) - - m = gradient(**data) - min = m["min"] - min_values = m["min_value"] - - print(min, min_values) - pprint(m["points"]) - - -if __name__ == "__main__": - main() \ No newline at end of file diff --git a/myPenaltyMethod.py b/myPenaltyMethod.py deleted file mode 100644 index f844ac8..0000000 --- a/myPenaltyMethod.py +++ /dev/null @@ -1,289 +0,0 @@ -## Penalty Function method. -## Loading the libraries. -import numpy as np -from utils import * - -sol_points = [] -eps = 0.000001 #Global Var. -delta = 0.001 #Global Var. -nc = 1 ## number of constraints. -g = np.zeros(nc) - -# Devuelve el valor de la función penalizada Φ(x, u) en un punto -variables, function, constraints = read_json('settings.json') - -def multi_f(x_vector): - global variables, function, constraints - sum_ = func_eval(variables, x_vector, function) - g = get_evaluated_constraints(variables, x_vector, constraints) - - for i in range(nc): - if(g[i] < 0.0): ## meaning that the constraint is violated. - sum_ = sum_ + r*g[i]*g[i] - - return(sum_) - -## This function will give us the gradient vector.. -def grad_multi_f(grad, x_ip): - d1 = np.zeros(M) - d2 = np.zeros(M) - delta_ = 0.001 - for i in range(M): - for j in range(M): - d1[j] = x_ip[j] - d2[j] = x_ip[j] - d1[i] = d1[i] + delta_ - d2[i] = d2[i] - delta_ - - fdiff = multi_f(d1) - bdiff = multi_f(d2) - grad[i] = (fdiff-bdiff)/(2*delta_) - return grad - - -## This function finds the alphastar <<-- UNIDIRECTIONAL SEARCH,, -def uni_search(guess, x_0, s_0, a1, b1): - xd = np.zeros(M) - xw = a1 + (b1-a1)*guess - for i in range(M): - xd[i] = x_0[i] + xw*s_0[i] - return multi_f(xd) - -## BOUNDING PHASE METHOD,, TO BRACKET THE ALPHASTAR. -def bracketing_(x_0, s_0): - k=0 - - ## Step__1 - ## "choose an initial guess and increment(delta) for the Bounding phase method". - del_=0.0000001; w_0=0.5;#initial guess for bounding phase. - - while(1): - f0 = uni_search(w_0, x_0, s_0, 0.0, 1.0) - fp = uni_search(w_0+del_, x_0, s_0, 0.0, 1.0) - fn = uni_search(w_0-del_, x_0, s_0, 0.0, 1.0) - - if(fn >= f0): - if(f0 >= fp): - del_ = abs(del_) ## step__2 - break - - else: - a=w_0-del_ - b=w_0+del_ - - elif((fn <= f0) & (f0 <= fp)): - del_ = -1*abs(del_) - break - - else: - del_ /= 2.0 - - wn = w_0 - del_## wn = lower limit and w_1 = upper limit. - - ## Step__3 - w_1=w_0 + del_*pow(2,k) ## << exponential purterbation - f1=uni_search(w_1, x_0, s_0, 0.0, 1.0) - - ## Step__4 - while(f1 < f0): - k+=1 - wn=w_0 - fn=f0 - w_0=w_1 - f0=f1 - w_1=w_0 + del_*pow(2,k) - f1=uni_search(w_1, x_0, s_0, 0.0, 1.0) - - a=wn - b=w_1 - ## if del_ is + we will bracket the minima in (a, b) - ## If del_ was -, - if(b < a): - temp=a;a=b;b=temp;#swapp - return(a, b) - -# In order to apply Secant Method(In order to compute z). We need to also calculate the derivative of the F(alpha) function, -def f_dash(x_0, s_0, xm): - xd = np.zeros(M) - - # using central difference. - for i in range(M): - xd[i] = x_0[i] + (xm+delta)*s_0[i] - # fdif = multi_f(xd) - fdif = multi_f(xd) - - for i in range(M): - xd[i] = x_0[i] + (xm-delta)*s_0[i] - # bdif = multi_f(xd) - bdif = multi_f(xd) - - f_ = (fdif - bdif)/(2*delta) - return(f_) - -## This is a function to compute the "z" used in the formula of SECANT Method.. -def compute_z(x_0, s_0, x1, x2): - z_ = x2 - ((x2-x1)*f_dash(x_0, s_0, x2))/(f_dash(x_0, s_0, x2)-f_dash(x_0, s_0,x1)) - return(z_) - -def secant_minima(a, b, x_0, s_0): - ## Step__1 - x1=a - x2=b - - ##Step__2 --> Compute New point. - z = compute_z(x_0, s_0, x1, x2) - - ##Step__3 - while(abs(f_dash(x_0, s_0, z)) > eps): - if(f_dash(x_0, s_0, z) >= 0): - x2=z # i.e. eliminate(z, x2) - else: - x1=z # i.e. eliminate(x1, z) - z = compute_z(x_0, s_0, x1, x2) #Step__2 - - return(z) - -def DFP(x_ip): - eps1 = 0.0001 - eps2 = 0.0001 - - grad0=np.zeros(M) - grad1=np.zeros(M) - xdiff=np.zeros(M) - Ae=np.zeros(M) - graddiff=np.zeros(M) - x_1=np.zeros(M) - x_2=np.zeros(M) - - ## step_1 - k=0 - x_0=x_ip #This is the initial Guess_. - - ## step_2 - # s_0=np.zeros(M) ##Search Direction. - # grad0 = grad_multi_f(grad0, x_ip) - # s_0 = -grad0 - s_0 = -grad_multi_f(grad0, x_ip) - - ## step_3 --> unidirectional search along s_0 from x_0. Aim -- to find alphastar and eventually x(1). - ## find alphastar such that f(x + alphastar*s) is minimum. - a,b = bracketing_(x_0, s_0) - alphastar = secant_minima(a, b, x_0, s_0) - - for i in range(M): - x_1[i] = x_0[i] + alphastar*s_0[i] - - grad1 = grad_multi_f(grad1, x_1) # Now computing grad(x(1)).. - - ## step_4 - ## Lets initialize the matrix A, - A = np.eye(M) #Identity Matrix. - dA = np.zeros((M,M)) - A1 = np.zeros((M,M)) - - for j in range(50): # 50 is the total iterations we will perform. - k+=1 - #print("Iteration Number -- ", k) - for i in range(M): - xdiff[i]=x_1[i]-x_0[i] - graddiff[i]= grad1[i]-grad0[i] - - # Now, Applying that GIANT FORMULA.. - dr=np.inner(xdiff, graddiff) - for i in range(M): - Ae[i]=0 - for j in range(M): - dA[i][j] = (xdiff[i]*xdiff[j])/dr - Ae[i] += (A[i][j]*graddiff[j]) - - dr=np.inner(graddiff, Ae) - for i in range(M): - for j in range(M): - A1[i][j] = A[i][j] + dA[i][j] - (Ae[i]*Ae[j])/dr - A[i][j] =A1[i][j] - - ## Finding the New Search Direction s_1.. -A_1*grad(x_1) - s_1 =np.zeros(M) - for i in range(M): - s_1[i]=0 - for j in range(M): - s_1[i] += (-A1[i][j]*grad1[j]) - - ## Unit vector along s_1; - unitv = np.linalg.norm(s_1) - for j in range(M): - s_1[j] /= unitv - - ## Step__5 - # Performing Unidirectional Search along s_1 from x_1; to find x_2;) - a,b=bracketing_(x_1, s_1) - alphastar = secant_minima(a, b, x_1, s_1) - alphastar = (int)(alphastar * 1000.0)/1000.0 #Rounding alphastar to three decimal places. - for i in range(M): - x_0[i] = x_1[i] - x_1[i] = x_1[i] + alphastar*s_1[i] - - ## Step__6 --->> TERMINATION CONDITION,, - grad1 = grad_multi_f(grad1, x_1) - grad0 = grad_multi_f(grad0, x_0) - - for j in range(M): - x_2[j] = x_1[j] - x_0[j] - - if((np.linalg.norm(grad1) <= eps1) or (np.linalg.norm(x_2)/np.linalg.norm(x_0)) <= eps2): - break - - return(x_1) - -if __name__ == "__main__": - func_key_name = get_key_name(function, constraints) - x_ip = np.array([0.11, 0.1]) - M = 2 ## Specifies the total dimension we are working with.. //Global Var. - c = 1.55 ##factor for updating r. - numseq = 20 ## number of sequences for penalty function method. - r = 0.1 - grad = np.zeros(M) - sol = np.zeros(M) - sol_ = [] - x_2 = np.zeros(M) - ## has the principal values in the key name. - func_key_name += ' ('+str(x_ip)[1:len(str(x_ip))-1]+') '+str(numseq)+' '+str(r)+' '+str(c) - dic = load_saved_data(func_key_name) - if dic is None: - - ## BRACKET OPERATOR PENALTY METHOD..1.1,1.1 - for k in range(numseq): - print('\n') - print('sequence number ', k) - print('current function value is ', multi_f(x_ip)) - x_1 = np.copy(x_ip) - sol = DFP(x_ip) - x_ip = np.copy(x_1) - - for j in range(M): - x_2[j] = sol[j] - x_ip[j] - track = np.linalg.norm(x_2)/np.linalg.norm(x_ip) - sol_ = [] - for i in range(M): - print(sol[i]) - sol_.append(sol[i]) - x_ip[i] = sol[i] - sol_points.append(sol_) - r*= c - if track <= 0.001: - break - sol_.append(func_eval(variables, sol_, function)) - dic = { - func_key_name: { - 'iterations': k, - 'points': str(sol_points), - 'min': str(sol_), - } - } - save_process(dic) - else: - dic = load_saved_data(func_key_name) - print('\n') - print('sequence number ', dic['iterations']) - print('the sequence of points was ', '\n['.join(dic['points'][1:len(dic['points'])-1].split('['))+'\n') - print('the final solution was ', dic['min']) \ No newline at end of file diff --git a/newton.py b/newton.py deleted file mode 100644 index 2f4890c..0000000 --- a/newton.py +++ /dev/null @@ -1,55 +0,0 @@ -from functools import reduce -import sympy as sp -from scipy.optimize import minimize -import numpy as np -# from pprint import pprint -# import json - - -def convert_list_to_tuples(list): - return reduce(lambda a, b: a + (b,), list, ()) - - -def newton(vars, func, initial_point, cycles=100, verbose=False): - vars = convert_list_to_tuples([sp.symbols(v) for v in vars]) - - func = sp.parse_expr(func) - func_lambda = sp.Lambda(vars, func) - func_evaluated = lambda x: np.array([func_lambda(*x)], dtype=float) - - gradient_lambda = np.array([sp.Lambda(vars, sp.diff(func, var)) for var in vars]) - gradient_evaluated = lambda x: np.array( - [g(*x) for g in gradient_lambda], dtype=float - ) - - initial_point = np.array(initial_point, dtype=float) - - m = minimize( - func_evaluated, - initial_point, - method="Newton-CG", - jac=gradient_evaluated, - options={"return_all": True, "maxiter": cycles}, - ) - - return { - "min": m["x"], - "min_value": m["fun"], - "points": m["allvecs"], - } - - -# def main(): -# f = open("newton.json") -# data = json.load(f) - -# m = newton(**data) -# min = m["min"] -# min_values = m["min_value"] - -# print(min, min_values) -# pprint(m["points"]) - - -# if __name__ == "__main__": -# main() \ No newline at end of file diff --git a/penalty.py b/penalty.py deleted file mode 100644 index 8ee042e..0000000 --- a/penalty.py +++ /dev/null @@ -1,262 +0,0 @@ - -## Penalty Function method. -## Loading the libraries. -import numpy as np - -eps = 0.000001 #Global Var. -delta = 0.001 #Global Var. -nc = 1 ## number of constraints. -g = np.zeros(nc) - -## This function will give the functional value at a point, -def multi_f(x): - ## Himmel Blau function! - sum_ = pow((pow(x[0],2) + x[1] - 11),2) + pow((pow(x[1],2) + x[0] - 7),2) - g[0] = -26.0 + pow((x[0]-5.0), 2) + pow(x[1],2);#constraints. - - for i in range(nc): - if(g[i] < 0.0): ## meaning that the constraint is violated. - sum_ = sum_ + r*g[i]*g[i]; - - return(sum_) - -## This function will give us the gradient vector.. -def grad_multi_f(grad, x_ip): - d1 = np.zeros(M); - d2 = np.zeros(M); - delta_=0.001; - - for i in range(M): - for j in range(M): - d1[j]=x_ip[j]; d2[j]=x_ip[j]; - - d1[i] = d1[i] + delta_; - d2[i] = d2[i] - delta_; - - fdiff = multi_f(d1); bdiff = multi_f(d2); - grad[i] = (fdiff - bdiff)/(2*delta_); - - return(grad) - -## This function finds the alphastar <<-- UNIDIRECTIONAL SEARCH,, -def uni_search(guess, x_0, s_0, a1, b1): - xd = np.zeros(M); - xw = a1 + (b1-a1)*guess; - for i in range(M): - xd[i] = x_0[i] + xw*s_0[i] - return multi_f(xd) - -## BOUNDING PHASE METHOD,, TO BRACKET THE ALPHASTAR. -def bracketing_(x_0, s_0): - k=0; - - ## Step__1 - ## "choose an initial guess and increment(delta) for the Bounding phase method". - del_=0.0000001; w_0=0.5;#initial guess for bounding phase. - - while(1): - f0 = uni_search(w_0, x_0, s_0, 0.0, 1.0); - fp = uni_search(w_0+del_, x_0, s_0, 0.0, 1.0); - fn = uni_search(w_0-del_, x_0, s_0, 0.0, 1.0); - - if(fn >= f0): - if(f0 >= fp): - del_ = abs(del_); ## step__2 - break; - - else: - a=w_0-del_; - b=w_0+del_; - - elif((fn <= f0) & (f0 <= fp)): - del_ = -1*abs(del_); - break; - - else: - del_ /= 2.0; - - wn = w_0 - del_;## wn = lower limit and w_1 = upper limit. - - ## Step__3 - w_1=w_0 + del_*pow(2,k); ## << exponential purterbation - f1=uni_search(w_1, x_0, s_0, 0.0, 1.0); - - ## Step__4 - while(f1 < f0): - k+=1; - wn=w_0; - fn=f0; - w_0=w_1; - f0=f1; - w_1=w_0 + del_*pow(2,k); - f1=uni_search(w_1, x_0, s_0, 0.0, 1.0); - - a=wn; - b=w_1; - ## if del_ is + we will bracket the minima in (a, b); - ## If del_ was -, - if(b < a): - temp=a;a=b;b=temp;#swapp - return(a, b) - -# In order to apply Secant Method(In order to compute z). We need to also calculate the derivative of the F(alpha) function, -def f_dash(x_0, s_0, xm): - xd = np.zeros(M); - - # using central difference. - for i in range(M): - xd[i] = x_0[i] + (xm+delta)*s_0[i]; - fdif = multi_f(xd); - - for i in range(M): - xd[i] = x_0[i] + (xm-delta)*s_0[i]; - bdif = multi_f(xd); - - f_ = (fdif - bdif)/(2*delta); - return(f_) - -## This is a function to compute the "z" used in the formula of SECANT Method.. -def compute_z(x_0, s_0, x1, x2): - z_ = x2 - ((x2-x1)*f_dash(x_0, s_0, x2))/(f_dash(x_0, s_0, x2)-f_dash(x_0, s_0,x1)); - return(z_) - -def secant_minima(a, b, x_0, s_0): - ## Step__1 - x1=a; - x2=b; - - ##Step__2 --> Compute New point. - z = compute_z(x_0, s_0, x1, x2); - - ##Step__3 - while(abs(f_dash(x_0, s_0, z)) > eps): - if(f_dash(x_0, s_0, z) >= 0): - x2=z; # i.e. eliminate(z, x2) - else: - x1=z; # i.e. eliminate(x1, z) - z = compute_z(x_0, s_0, x1, x2); #Step__2 - - return(z) - -def DFP(x_ip): - eps1 = 0.0001; - eps2 = 0.0001; - - grad0=np.zeros(M); - grad1=np.zeros(M); - xdiff=np.zeros(M); - Ae=np.zeros(M); - graddiff=np.zeros(M); - x_1=np.zeros(M); - x_2=np.zeros(M); - - ## step_1 - k=0; - x_0=x_ip #This is the initial Guess_. - - ## step_2 - s_0=np.zeros(M); ##Search Direction. - grad0 = grad_multi_f(grad0, x_ip); - s_0 = -grad0 - - ## step_3 --> unidirectional search along s_0 from x_0. Aim -- to find alphastar and eventually x(1). - ## find alphastar such that f(x + alphastar*s) is minimum. - a,b = bracketing_(x_0, s_0); - alphastar = secant_minima(a, b, x_0, s_0); - - for i in range(M): - x_1[i] = x_0[i] + alphastar*s_0[i]; - - grad1 = grad_multi_f(grad1, x_1); # Now computing grad(x(1)).. - - ## step_4 - ## Lets initialize the matrix A, - A = np.eye(M) #Identity Matrix. - dA = np.zeros((M,M)) - A1 = np.zeros((M,M)) - - for j in range(50): # 50 is the total iterations we will perform. - k+=1; - #print("Iteration Number -- ", k) - for i in range(M): - xdiff[i]=x_1[i]-x_0[i]; - graddiff[i]= grad1[i]-grad0[i]; - - # Now, Applying that GIANT FORMULA.. - dr=np.inner(xdiff, graddiff); - for i in range(M): - Ae[i]=0; - for j in range(M): - dA[i][j] = (xdiff[i]*xdiff[j])/dr; - Ae[i] += (A[i][j]*graddiff[j]); - - dr=np.inner(graddiff, Ae) - for i in range(M): - for j in range(M): - A1[i][j] = A[i][j] + dA[i][j] - (Ae[i]*Ae[j])/dr; - A[i][j] =A1[i][j]; - - ## Finding the New Search Direction s_1.. -A_1*grad(x_1) - s_1 =np.zeros(M); - for i in range(M): - s_1[i]=0; - for j in range(M): - s_1[i] += (-A1[i][j]*grad1[j]); - - ## Unit vector along s_1; - unitv = np.linalg.norm(s_1); - for j in range(M): - s_1[j] /= unitv; - - ## Step__5 - # Performing Unidirectional Search along s_1 from x_1; to find x_2;) - a,b=bracketing_(x_1, s_1); - alphastar = secant_minima(a, b, x_1, s_1); - alphastar = (int)(alphastar * 1000.0)/1000.0; #Rounding alphastar to three decimal places. - for i in range(M): - x_0[i] = x_1[i]; - x_1[i] = x_1[i] + alphastar*s_1[i]; - - ## Step__6 --->> TERMINATION CONDITION,, - grad1 = grad_multi_f(grad1, x_1); - grad0 = grad_multi_f(grad0, x_0); - - for j in range(M): - x_2[j] = x_1[j] - x_0[j]; - - if((np.linalg.norm(grad1) <= eps1) or (np.linalg.norm(x_2)/np.linalg.norm(x_0)) <= eps2): - break; - - return(x_1) - -if __name__ == "__main__": - M = 2 ## Specifies the total dimension we are working with.. //Global Var. - c=1.55; ## factor for updating r. - numseq = 20; ## number of sequences for penalty function method. - r = 0.1 - - x_ip = np.array([0.11,0.1])## Our initial guess, - grad = np.zeros(M) - sol = np.zeros(M); - x_2 = np.zeros(M); - - ## BRACKET OPERATOR PENALTY METHOD..1.1,1.1 - for k in range(numseq): - print("\n") - print("sequence number ", k) - print("current function value is ", multi_f(x_ip)) - x_1 = np.copy(x_ip); - sol = DFP(x_ip); - x_ip = np.copy(x_1); - - for j in range(M): - x_2[j] = sol[j] - x_ip[j]; - - track = np.linalg.norm(x_2)/np.linalg.norm(x_ip); - for i in range(M): - print(sol[i]); - x_ip[i] = sol[i] - - r *= c; - if(track <= 0.001):# TERMINATION CONDITION. - break; \ No newline at end of file diff --git a/penalty_gradient.py b/penalty_gradient.py deleted file mode 100644 index 04a985b..0000000 --- a/penalty_gradient.py +++ /dev/null @@ -1,88 +0,0 @@ -## Penalty Function method. -## Loading the libraries. -from scipy.optimize import line_search -from gradient import gradient -from utils import * -import numpy as np - - -# Devuelve el valor de la función penalizada Φ(x, u) en un punto - -# minimize = 0, maximize = 1 -def penalty_gradient(): - with open('penalty_settings.json') as settings: - data = json.load(settings) - numseq = data['Penalty_number_of_sequence'] - r = data['Penalty_penalty_factor'] - c = data['Penalty_update_factor'] - constraints = data['Penalty_constraints'] - min_or_max = data['Penalty_max_or_min'] - x_ip = data['Penalty_init_point'] - function = data['Penalty_func'] - variables = data['Penalty_vars'] - - variables = convert_list_to_tuples(variables) - M = len(variables) ## Specifies the total dimension we are working with.. //Global Var. - sol = np.zeros(M) - x_2 = np.zeros(M) - penalty_func = '' - sol_ = [] - sol_points = [] - - if min_or_max == 1: - function = '-('+function+')' - - func_key_name = get_key_name(function, constraints) - ## has the principal values in the key name. - func_key_name += ' ('+str(x_ip)[1:len(str(x_ip))-1]+') '+str(numseq)+' '+str(r)+' '+str(c) - dic = load_saved_data(func_key_name) - - if dic is None: - for k in range(numseq): - print('\n') - print('sequence number ', k) - _, func, penalty_func, penalty_lambda, _, _ = get_penalty_func(function, constraints, variables, r, min_or_max) - - sum_ = penalty_lambda(*x_ip) - print('current function value is ', sum_) - x_1 = np.copy(x_ip) - m = gradient(variables, str(penalty_func), x_ip) - sol = m["min"] - x_ip = np.copy(x_1) - - for j in range(M): - x_2[j] = sol[j] - x_ip[j] - - track = np.linalg.norm(x_2)/np.linalg.norm(x_ip) - sol_ = [] - - for i in range(M): - print(sol[i]) - sol_.append(sol[i]) - x_ip[i] = sol[i] - - sol_points.append(sol_) - r*= c - - if track <= 0.001: - break - - sol_.append(func_eval(variables, sol_, func)) - dic = { - func_key_name: { - 'iterations': k, - 'points': str(sol_points), - 'min': str(sol_), - } - } - save_process(dic) - return dic - else: - dic = load_saved_data(func_key_name) - print('\n') - print('sequence number ', dic['iterations']) - print('the sequence of points was ', '\n['.join(dic['points'][1:len(dic['points'])-1].split('['))+'\n') - print('the final solution was ', dic['min']) - return dic - -# dict_ = penalty_gradient() \ No newline at end of file diff --git a/penalty_newton.py b/penalty_newton.py deleted file mode 100644 index aaee8a4..0000000 --- a/penalty_newton.py +++ /dev/null @@ -1,88 +0,0 @@ -## Penalty Function method. -## Loading the libraries. -from scipy.optimize import line_search -from newton import newton -from utils import * -import numpy as np - - -# Devuelve el valor de la función penalizada Φ(x, u) en un punto - -# minimize = 0, maximize = 1 -def penalty_newton(): - with open('penalty_settings.json') as settings: - data = json.load(settings) - numseq = data['Penalty_number_of_sequence'] - r = data['Penalty_penalty_factor'] - c = data['Penalty_update_factor'] - constraints = data['Penalty_constraints'] - min_or_max = data['Penalty_max_or_min'] - x_ip = data['Penalty_init_point'] - function = data['Penalty_func'] - variables = data['Penalty_vars'] - - variables = convert_list_to_tuples(variables) - M = len(variables) ## Specifies the total dimension we are working with.. //Global Var. - sol = np.zeros(M) - x_2 = np.zeros(M) - penalty_func = '' - sol_ = [] - sol_points = [] - - if min_or_max == 1: - function = '-('+function+')' - - func_key_name = get_key_name(function, constraints) - ## has the principal values in the key name. - func_key_name += ' ('+str(x_ip)[1:len(str(x_ip))-1]+') '+str(numseq)+' '+str(r)+' '+str(c) - dic = load_saved_data(func_key_name) - - if dic is None: - for k in range(numseq): - print('\n') - print('sequence number ', k) - _, func, penalty_func, penalty_lambda, _, _ = get_penalty_func(function, constraints, variables, r)#min_or_max) - - sum_ = penalty_lambda(*x_ip) - print('current function value is ', sum_) - x_1 = np.copy(x_ip) - m = newton(variables, str(penalty_func), x_ip) - sol = m["min"] - x_ip = np.copy(x_1) - - for j in range(M): - x_2[j] = sol[j] - x_ip[j] - - track = np.linalg.norm(x_2)/np.linalg.norm(x_ip) - sol_ = [] - - for i in range(M): - print(sol[i]) - sol_.append(sol[i]) - x_ip[i] = sol[i] - - sol_points.append(sol_) - r*= c - - if track <= 0.001: - break - - sol_.append(func_eval(variables, sol_, func)) - dic = { - func_key_name: { - 'iterations': k, - 'points': str(sol_points), - 'min': str(sol_), - } - } - save_process(dic) - return dic - else: - dic = load_saved_data(func_key_name) - print('\n') - print('sequence number ', dic['iterations']) - print('the sequence of points was ', '\n['.join(dic['points'][1:len(dic['points'])-1].split('['))+'\n') - print('the final solution was ', dic['min']) - return dic - -# dict_ = penalty_newton() \ No newline at end of file diff --git a/penalty_settings.json b/penalty_settings.json deleted file mode 100644 index e1a077b..0000000 --- a/penalty_settings.json +++ /dev/null @@ -1,12 +0,0 @@ -{ - "Penalty_number_of_sequence": 20, - "Penalty_penalty_factor": 0.1, - "Penalty_update_factor": 0.155, - "Penalty_constraints" : ["x>=0", "y>=0", "(x-5)**2 + y**2 - 26 >= 0"], - "Penalty_max_or_min": 0, - "Penalty_init_point": [0.11, 0.11], - "Penalty_x_range": [-20, 20, 1], - "Penalty_y_range": [-20, 20, 1], - "Penalty_func" : "(x**2 + y - 11)+(x + y**2 - 7)**2", - "Penalty_vars" : ["x", "y"] -} \ No newline at end of file diff --git a/settings.json b/settings.json deleted file mode 100644 index 19633c6..0000000 --- a/settings.json +++ /dev/null @@ -1,12 +0,0 @@ -{ - "Penalty_number_of_sequence": 20, - "Penalty_penalty_factor": 0.1, - "Penalty_update_factor": 0.155, - "Penalty_constraints" : ["x>=2", "y<=10"], - "Penalty_max_or_min": 0, - "Penalty_init_point": [2.3, 3], - "Penalty_x_range": [-20, 20, 1], - "Penalty_y_range": [-20, 20, 1], - "Penalty_func" : "x**2+y**2", - "Penalty_vars" : ["x", "y"] -} \ No newline at end of file diff --git a/testnp.py b/testnp.py deleted file mode 100644 index 96c9109..0000000 --- a/testnp.py +++ /dev/null @@ -1,75 +0,0 @@ -import json -import sympy as sym -from sympy.parsing.sympy_parser import parse_expr - -g = list() -# x = list() -# x.append("x0") -# for i in range(1,20): -# q="x"+str(i) -# x.append(sym.Symbol(q)) -# print(x) - -# f='x**2+2y+1'#input("f(x)=x**2+x") -# f=parse_expr(f) -# print(f.subs('x', 2)) -def move_inequality_constants(ineq): - l = ineq.lhs - r = ineq.rhs - op = ineq.rel_op - - if op.__contains__('<'): - return l - r - else: - return r - l -def func_eval(x_vector, value_vector, func): - f = func - for x_i, v_i in zip(x_vector, value_vector): - f = f.subs(x_i, v_i) - return f -# f1 = lambda x,y: x**2+y+1 -# print(f1(1,2)) -with open('settings.json') as settings: - data = json.load(settings) - -x = data['vars'] -f = data['func'] -constraints = data['constraints'] - -c = parse_expr(constraints[0]) -c = move_inequality_constants(c) -print(c) - -f = parse_expr(f) -values = [2,2] - -print(func_eval(x, values, f)) - -gfg = sym.Lambda(x, f) -print(gfg(2,2)) - -def get_constraints(constraints): - g = list() - for c in constraints: - g.append((move_inequality_constants(parse_expr(c)))) - - -# def multi_f(x_vector): -# with open('settings.json') as settings: -# data = json.load(settings) - -# x = data['vars'] -# f = data['func'] -# constraints = data['constraints'] -# c = parse_expr(constraints[0]) -# c = move_inequality_constants(c) -# function = parse_expr(f) -# sum_ = func_eval(x_vector, values, function) -# # sum_ = pow((pow(x[0],2) + x[1] - 11),2) + pow((pow(x[1],2) + x[0] - 7),2) -# g[0] = -26.0 + pow((x[0]-5.0), 2) + pow(x[1],2)#constraints. - -# for i in range(nc): -# if(g[i] < 0.0): ## meaning that the constraint is violated. -# sum_ = sum_ + r*g[i]*g[i] - -# return(sum_) \ No newline at end of file diff --git a/utils.py b/utils.py deleted file mode 100644 index 7813291..0000000 --- a/utils.py +++ /dev/null @@ -1,396 +0,0 @@ -import json -from functools import reduce -import sympy as sym -from sympy.parsing.sympy_parser import parse_expr, standard_transformations, convert_equals_signs -import math -def move_inequality_constants(ineq): - l = ineq.lhs - r = ineq.rhs - op = ineq.rel_op - - #revisar como es para minimo - # if op.__contains__('<'): - if op.__contains__('>'): - return l - r - else: - return r - l - -def func_eval(x_vector, value_vector, func): - f = func - for x_i, v_i in zip(x_vector, value_vector): - f = f.subs(x_i, v_i) - return f - -def get_lambda(x_vector, f): - return sym.Lambda(('x', 'y'), f) - # return sym.Lambda(x_vector, f) - -def get_constraints(constraints): - result = list() - for c in constraints: - result.append((move_inequality_constants(parse_expr(c)))) - return result - -def get_evaluated_constraints(variables, values, constraints): - result = list() - for con in constraints: - result.append(func_eval(variables, values, con)) - return result - -def read_json(path): - with open(path) as settings: - data = json.load(settings) - - f = data['Penalty_func'] - variables = data['Penalty_vars'] - c = data['Penalty_constraints'] - - constraints = get_constraints(c) - function = parse_expr(f) - - return variables, function, constraints - -def convert_list_to_tuples(list): - return reduce(lambda a, b: a + (b,), list, ()) - -def save_process(my_dict): - with open('data.json', 'r') as fp1: - data = json.load(fp1) - for key in my_dict.keys(): - data[key] = my_dict[key] - with open('data.json', 'w') as fp: - json.dump(data, fp) - -def load_saved_data(key): - with open('data.json', 'r') as fp: - data = json.load(fp) - try: - return data[key] - except KeyError: - return None - -def get_key_name(func, constraints): - f = str(func) + ' ' - for i, c in enumerate(constraints): - f+= str(c) + ' ' if i != len(constraints) else str(c) - return f - - -x = sym.Symbol('x') -y = sym.Symbol('y') -# def get_min_cleared_constraint(ineq): -# l = ineq.lhs -# r = ineq.rhs -# if len(l.free_symbols) == 2: -# # estan ambas dos en la izquierda aún -# newr = r -# for arg in l.args: -# if arg.free_symbols.__contains__(y): -# newl = arg -# continue -# if l.is_Add: -# newr -= arg -# elif l.is_Mul: -# newr /= arg -# return get_min_cleared_constraint(newl >= newr) -# if len(l.free_symbols) == 1: -# if l.free_symbols.__contains__(y): -# if len(l.args) == 0: -# return ineq -# else: -# newr = r -# for arg in l.args: -# if arg.free_symbols.__contains__(y): -# newl = arg -# continue -# if l.is_Add: -# newr -= arg -# elif l.is_Mul: -# newr /= arg -# return get_min_cleared_constraint(newl >= newr) -# elif l.free_symbols.__contains__(x): -# if len(l.args) == 0: -# return ineq -# else: -# newr = r -# for arg in l.args: -# if arg.free_symbols.__contains__(x): -# newl = arg -# continue -# if l.is_Add: -# newr -= arg -# elif l.is_Mul: -# newr /= arg -# return get_min_cleared_constraint(newl >= newr) - # if ineq.rel_op.__contains__('>'): - # return newl >= newr - # else: - # return newl <= newr - -# x**2 = 4 => x = log4_2 - -def get_cleared_constraint(ineq): - if len(ineq.free_symbols) >= 2: - if ineq.free_symbols.__contains__(y): - a = sym.solve(ineq, y) - return a.args[1] - else: - a = sym.solve(ineq, x) - return a.args[1] - else: - if ineq.free_symbols.__contains__(y): - a = sym.solve(ineq, y) - return a.args[0] - else: - a = sym.solve(ineq, x) - return a.args[0] -def get_eq_cleared_constraint(ineq): - if len(ineq.free_symbols) >= 2: - if ineq.free_symbols.__contains__(y): - a = sym.solve(ineq, y) - return a[0], y - else: - a = sym.solve(ineq, x) - return a[0], x - else: - if ineq.free_symbols.__contains__(y): - a = sym.solve(ineq, y) - return a[0], y - else: - a = sym.solve(ineq, x) - return a[0], x - -def move_inequality_to_min(ineq): - l = ineq.lhs - r = ineq.rhs - op = ineq.rel_op - - if op.__contains__('>'): - return ineq - else: - return -1*l >= -1*r -def move_inequality_to_max(ineq): - l = ineq.lhs - r = ineq.rhs - op = ineq.rel_op - - if op.__contains__('<'): - return ineq - else: - return -1*l <= -1*r -#max_or_min = 0 means doing minimun -def get_constraints_cleared(ineqs_, max_or_min=0): - ineqs = [get_cleared_constraint(ineq) for ineq in ineqs_] - # ineqs = [move_inequality_to_min(i) if max_or_min == 0 else move_inequality_to_max(i) for i in ineqs] - return ineqs -def get_lambdas(ineqs): - lam = [] - for i in ineqs: - l = i.lhs - r = i.rhs - if len (l.free_symbols) == 1: - syms = [s for s in r.free_symbols] - if len(syms) == 0: - syms = [x] if r.free_symbols.__contains__(y) else [y] - lam.append(sym.Lambda(syms, r)) - elif len (l.free_symbols) == 0: - syms = [x] if r.free_symbols.__contains__(y) else [y] - lam.append(sym.Lambda(syms, l)) - else: - syms = [s for s in l.free_symbols] - lam.append(sym.Lambda(syms, l)) - return lam -def get_eq(ineqs): - lam = [] - for i in ineqs: - l = i.lhs - r = i.rhs - if len (l.free_symbols) == 1: - lam.append(r) - else: - lam.append(l) - return lam - -def check_constraint_point(x_, y_, constraints_, lambdas_): - for c, l in zip(constraints_,lambdas_): - op = c[1] - _y = l(x_) - if op.__contains__('>') and y_ < _y: - return False - elif op.__contains__('<') and y_ > _y: - return False - return True -def check_eq_constraint_point(x_, y_, lambdas_): - for l in lambdas_: - _y = l(x_) - if _y != y_: - return False - return True -# ineqs = [] -# ineqs.append(parse_expr('20*x+50*y <= 3000')) -# ineqs.append(parse_expr('x+y <= 90')) -# ineqs.append(parse_expr('y >= 10')) -# ineqs.append(parse_expr('y >= 0')) -# ineqs.append(parse_expr('x >= 0')) -# ineqs.append(parse_expr('10000*x+6000*y >= z')) -# converted = get_min_constraints_cleared(ineqs) - -# a = 0 - -# def move_inequality_constants(ineq): -# l = ineq.lhs -# r = ineq.rhs -# op = ineq.rel_op - -# if op == "<=": -# return sym.GreaterThan(r - l, 0) -# elif op == "<": -# return sym.StrictGreaterThan(r - l, 0) -# elif op == ">=": -# return sym.GreaterThan(l - r, 0) -# elif op == ">": -# return sym.StrictGreaterThan(l - r, 0) -# elif op == "=": -# return sym.Eq(l - r, 0) -################## geometric ################### -def get_geometric_ineqs_and_eqs(constraints): - ineqs = [] - eqs = [] - for c in constraints: - con = parse_expr(c, transformations = standard_transformations+(convert_equals_signs,)) - l = con.lhs - r = con.rhs - op = con.rel_op - - if op == "<=": - ineqs.append(con) - elif op == "<": - ineqs.append(con) - elif op == ">=": - ineqs.append(con) - elif op == ">": - ineqs.append(con) - elif op == "==": - eqs.append(con) - return ineqs, eqs - -################### penalty ###################### -def get_penalty_ineqs_and_eqs(constraints): - ineqs = [] - eqs = [] - for c in constraints: - con = parse_expr(c, transformations = standard_transformations+(convert_equals_signs,)) - l = con.lhs - r = con.rhs - op = con.rel_op - - if op == "<=": - ineqs.append(con) - elif op == "<": - ineqs.append(con) - elif op == ">=": - ineqs.append(con) - elif op == ">": - ineqs.append(con) - elif op == "==": - eqs.append(con) - return ineqs, eqs - -def get_penalty_eq_constraints(eqs): - if len(eqs) == 0: - return [], None - result = list() - total_eq = eqs[0] - result.append(total_eq.lhs-total_eq.rhs) - total_eq = (result[0])**2 - for eq in eqs[1:]: - parsed_eq = eq - result.append(parsed_eq.lhs-parsed_eq.rhs) - total_eq += (parsed_eq.lhs-parsed_eq.rhs)**2 - return result, total_eq - -################### penalty for minimize -def move_inequality_constants_minimize(ineq): - l = ineq.lhs - r = ineq.rhs - op = ineq.rel_op - - if op.__contains__('>'): - return l - r - else: - return r - l - - -def get_penalty_ineq_constraints_minimize(ineqs): - if len(ineqs) == 0: - return [], None - result = list() - total_ineq = move_inequality_constants_minimize(ineqs[0]) - result.append(total_ineq) - total_ineq = (sym.Max(0, total_ineq))**2 - for c in ineqs[1:]: - ineq = move_inequality_constants_minimize(c) - result.append(ineq) - total_ineq += (sym.Max(0, ineq))**2 - return result, total_ineq - -def get_penalty_func_minimize(func, ineqs, eqs, variables, step): - penalty_func = func = parse_expr(func) - parsed_eqs, eqs_func = get_penalty_eq_constraints(eqs) - parsed_ineqs, ineqs_func = get_penalty_ineq_constraints_minimize(ineqs) - penalty_func = penalty_func + step * eqs_func if eqs_func is not None else penalty_func - penalty_func = penalty_func + step * ineqs_func if ineqs_func is not None else penalty_func - penalty_lambda = sym.Lambda((variables[0], variables[1]), penalty_func) - return func, penalty_func, penalty_lambda, parsed_eqs, parsed_ineqs - -################### penalty for maximize -def move_inequality_constants_maximize(ineq): - l = ineq.lhs - r = ineq.rhs - op = ineq.rel_op - - if op.__contains__('<'): - return l - r - else: - return r - l - -def get_penalty_ineq_constraints_maximize(ineqs): - if len(ineqs) == 0: - return [], None - result = list() - total_ineq = move_inequality_constants_maximize(ineqs[0]) - result.append(total_ineq) - total_ineq = sym.Max(0, total_ineq) - for c in ineqs[1:]: - result.append((move_inequality_constants_maximize(c))) - total_ineq += sym.Max(0, move_inequality_constants_maximize(c)) - return result, total_ineq - -def get_penalty_func_maximize(func, ineqs, eqs, variables, step): - penalty_func = func = parse_expr(func) - parsed_eqs, eqs_func = get_penalty_eq_constraints(eqs) - parsed_ineqs, ineqs_func = get_penalty_ineq_constraints_maximize(ineqs) - penalty_func = penalty_func + step * eqs_func if eqs_func is not None else penalty_func - penalty_func = penalty_func + step * ineqs_func if ineqs_func is not None else penalty_func - penalty_lambda = sym.Lambda((variables[0], variables[1]), penalty_func) - return func, penalty_func, penalty_lambda, parsed_eqs, parsed_ineqs - -################################################################# -# 0 for minimize, 1 for maximize -def get_penalty_func(func, constraints, variables, step, max_or_min = 0): - ineqs, eqs = get_penalty_ineqs_and_eqs(constraints) - if max_or_min == 0: # minimize - return ineqs, *get_penalty_func_minimize(func, ineqs, eqs, variables, step) - else: # maximize - return ineqs, *get_penalty_func_maximize(func, ineqs, eqs, variables, step) - -def read_dots_from_json(string_): - dots = [] - splited_text = (string_[2:len(string_)-2].split('[')) - for splited in splited_text[:len(splited_text)-1]: - new_splited = splited.split(', ') - dots.append((float(new_splited[0]), float(new_splited[1][:len(new_splited[1])-1]))) - a = 0 - last_row = splited_text[-1].split(', ') - dots.append((float(last_row[0]), float(last_row[1]))) - return dots \ No newline at end of file