bench_nemo_naïf.py

# -*- coding: utf-8 -*-
"""
Created on Thu Oct 31 15:22:51 2019

@author: David N.
"""

import random as rng
import numpy as np
import time
import matplotlib.pyplot as plt

def rand_in_range(mini, maxi):
    return mini + rng.random() * (maxi - mini)

def create_fog(n):
    return [np.array((rand_in_range(-10,10), rand_in_range(-10,10))) for _ in range(n)]

def center_circle(A,B,C):    
    I = (A+B)/2
    J = (B+C)/2
    AB = B-A
    BC = C-B
    a = -AB[0]/AB[1]
    b = -a*I[0]+I[1]
    ap= -BC[0]/BC[1]
    bp= -ap*J[0]+J[1]
    return np.linalg.solve([[a,-1.],[ap,-1]],[-b,-bp])

def dist_sq(A,B):
    return sum((A-B)**2)

def dist(A,B):
    return np.sqrt(dist_sq(A,B))

n = 2
delta_t = 0.
N = create_fog(n)

list_n = []
list_delta_t = []

while delta_t < 10.:
    n += 1
    N.append(np.array((rand_in_range(-10,10), rand_in_range(-10,10))))
    
    start = time.time()
    
    max_A,max_B,max_C = (),(),()
    max_O = ()
    max_r = 0.
    
    for i in range(n-2):
        for j in range(i+1, n-1):
            for k in range(j+1, n):
                A,B,C = N[i],N[j],N[k]
                valid = True
                
                O = center_circle(A,B,C)
                r_sq = dist_sq(O,A)
                
                for P in N:
                    if (P == A).all() or (P == B).all() or (P == C).all():
                        continue
                    valid = dist_sq(O,P) >= r_sq
                    if not valid:
                        break
                
                r = np.sqrt(r_sq)
                
                if valid and r > max_r:
                    max_A,max_B,max_C = A,B,C
                    max_O = O
                    max_r = r
                    
    delta_t = time.time() - start
    
    print("Nuage de taille ", n, " résolu en ", delta_t, "s")
    
    list_n.append(n)
    list_delta_t.append(delta_t)
    
plt.plot(list_n, list_delta_t, label='mesure')

coef = list_delta_t[-1] / (list_n[-1]**2)
plt.plot(list_n, [coef*(n**2) for n in list_n], label='n^2')

coef = list_delta_t[-1] / (list_n[-1]**3)
plt.plot(list_n, [coef*(n**3) for n in list_n], label='n^3')

coef = list_delta_t[-1] / (list_n[-1]**4)
plt.plot(list_n, [coef*(n**4) for n in list_n], label='n^4')

plt.show()