Biljana

import numpy as np
from math import cos, sin, sqrt
from random import randint

def f_squared(t, u, v, p1, p2, p3):
    # Computes the squared area of the triangle
    # defined by p1(t), p2(u), and p3(v)
    # with circle radii r1, r2, and r3
    a = np.linalg.norm(np.cross(p2(u), p3(v)))
    b = np.linalg.norm(np.cross(p3(v), p1(t)))
    c = np.linalg.norm(np.cross(p1(t), p2(u)))
    s = (a + b + c) / 2
    area_squared = s * (s - a) * (s - b) * (s - c)
    return area_squared

def grad_f_squared(t, u, v, p1, p2, p3):
    # Computes the gradient of -f^2 at (t, u, v)
    # with circle radii r1, r2, and r3
    grad = np.zeros(3)
    grad[0] = 0.5 * np.dot(p1(t) - p2(u), np.cross(p2(u), p3(v)))[0] / f_squared(t, u, v, p1, p2, p3)
    grad[1] = 0.5 * np.dot(p2(u) - p3(v), np.cross(p3(v), p1(t)))[0] / f_squared(t, u, v, p1, p2, p3)
    grad[2] = 0.5 * np.dot(p3(v) - p1(t), np.cross(p1(t), p2(u)))[0] / f_squared(t, u, v, p1, p2, p3)
    return -2 * grad

def steepest_descent(p1, p2, p3, epsilon=1e-2, max_iter=100, step_size=0.05):
    # Performs steepest descent to find the minimum value of -f^2
    # with circle radii r1, r2, and r3
    # Returns the optimal (t, u, v) values and the corresponding minimum value
    t, u, v = 0.0, 0.0, 0.0
    t, u, v = randint(0, 628) / 100, randint(0, 628) / 100, randint(0, 628) / 100
    for i in range(max_iter):
        grad = grad_f_squared(t, u, v, p1, p2, p3)
        if np.linalg.norm(grad) < epsilon:
            break
        t -= step_size * grad[0]
        u -= step_size * grad[1]
        v -= step_size * grad[2]
    return (t, u, v), abs(-f_squared(t, u, v, p1, p2, p3))

K = [[3, 4, 2], [0, 0, 1], [-3, 3, 1.5]]

def p1(t):
    return np.array([K[0][0] + K[0][2] * cos(t), K[0][1] + K[0][2] * sin(t)])
def p2(t):
    return np.array([K[1][0] + K[1][2] * cos(t), K[1][1] + K[1][2] * sin(t)])
def p3(t):
    return np.array([K[2][0] + K[2][2] * cos(t), K[2][1] + K[2][2] * sin(t)])

best_sol = None
best_area = 0
for i in range(100):
    sol, area = steepest_descent(p1, p2, p3)
    print("Area: ", sqrt(area))
    print(p1(sol[0]), p2(sol[1]), p3(sol[2]))
    if area > best_area:
        best_area = area
        best_sol = sol

print("Best Area: ", sqrt(best_area))
print(p1(best_sol[0]), p2(best_sol[1]), p3(best_sol[2]))