# Fonctions Floues

unknown
plain_text
a year ago
42 kB
5
Indexable
Never
```from fonctions_autres import *

import numpy as np

import matplotlib.pyplot as plt
import scipy

import skfuzzy as fuzz
from skfuzzy import control as ctrl

from tqdm import tqdm

def input_fuzz(variable, input_to_fuzz):
fuzz_input = []
for vari in variable.terms.keys():
fuzz_input.append(fuzz.interp_membership(variable.universe, variable.terms[vari].mf, input_to_fuzz))
return fuzz_input

def IRR_2var(SF_rules, input_fuzz1, input_fuzz2):
IRR = np.empty([len(SF_rules), 2])
k = -1
for j in range(0, len(SF_rules)):
if (j % len(input_fuzz2)) == 0:
k += 1
IRR[j, 0] = input_fuzz1[k]
IRR[j, 1] = input_fuzz2[j % len(input_fuzz2)]

return IRR

def IRR_3var(SF_rules, input_fuzz1, input_fuzz2, input_fuzz3):
IRR = np.empty([len(SF_rules), 3])

k = -1
m = -1
for j in range(0, len(SF_rules)):
if (j % len(input_fuzz3)) == 0:
k += 1
if (k % len(input_fuzz2)) == 0:
m += 1
IRR[j, 0] = input_fuzz1[m]
IRR[j, 1] = input_fuzz2[j % len(input_fuzz2)]
IRR[j, 2] = input_fuzz3[j % len(input_fuzz3)]
return IRR

def output_fuzz(csq, output_var):
first_rule = True
output = np.zeros(1)
for membi in output_var.terms:
csq_i = '[' + str(output_var.terms[membi]) + ']'
if csq[csq_i] != 0:
x_csq = fuzz.fuzzymath.interp_universe(output_var.universe, output_var.terms[membi].mf, 1)
if first_rule:
output = np.array([[np.max(x_csq), csq[csq_i]]])
first_rule = False
else:
output = np.append(output, [[(np.min(x_csq)), csq[csq_i]]], axis=0)
return output

def SF1_compute(input_temperature, input_d_speed, input_time, want_to_plot):
pas = 0.01
temperature = ctrl.Antecedent(np.arange(263, 323, pas), 'Temperature (K)')
d_speed = ctrl.Antecedent(np.arange(-100, 100, pas), 'Variation of wind speed (m/s)')
timeTC = ctrl.Antecedent(np.arange(0, 15, pas), 'Time since TC formation (days)')
lifespan = ctrl.Consequent(np.arange(0, 15, pas), 'Lifespan')

temperature['Low'] = fuzz.trapmf(temperature.universe, [263, 263, 293, 303])
temperature['Ok'] = fuzz.trapmf(temperature.universe, [293, 303, 323, 323])

d_speed['--'] = fuzz.trapmf(d_speed.universe, [-100, -100, -60, -20])
d_speed['-'] = fuzz.trimf(d_speed.universe, [-60, -30, 0])
d_speed['St'] = fuzz.trimf(d_speed.universe, [-30, 0, 20])
d_speed['+'] = fuzz.trimf(d_speed.universe, [0, 30, 60])
d_speed['++'] = fuzz.trapmf(d_speed.universe, [20, 60, 100, 100])

timeTC['Young'] = fuzz.zmf(timeTC.universe, 3, 6)
timeTC['Middle'] = fuzz.gaussmf(timeTC.universe, 6, 1)
timeTC['Old'] = fuzz.smf(timeTC.universe, 6, 9)

lifespan['Dying+'] = fuzz.trapmf(lifespan.universe, [0, 0, 3, 5])
lifespan['Dying'] = fuzz.trimf(lifespan.universe, [3, 5, 7])
lifespan['St'] = fuzz.trimf(lifespan.universe, [5, 7, 10])
lifespan['Living'] = fuzz.trimf(lifespan.universe, [7, 10, 13])
lifespan['Living+'] = fuzz.trapmf(lifespan.universe, [10, 13, 15, 15])

if want_to_plot:
lifespan.view()
temperature.view()
d_speed.view()
timeTC.view()

SF1_rules = [ctrl.Rule(temperature['Low'] & d_speed['--'] & timeTC['Young'], lifespan['Dying+']),
ctrl.Rule(temperature['Low'] & d_speed['--'] & timeTC['Middle'], lifespan['Dying+']),
ctrl.Rule(temperature['Low'] & d_speed['--'] & timeTC['Old'], lifespan['Dying+']),
ctrl.Rule(temperature['Low'] & d_speed['-'] & timeTC['Young'], lifespan['Dying']),
ctrl.Rule(temperature['Low'] & d_speed['-'] & timeTC['Middle'], lifespan['Dying+']),
ctrl.Rule(temperature['Low'] & d_speed['-'] & timeTC['Old'], lifespan['Dying+']),
ctrl.Rule(temperature['Low'] & d_speed['St'] & timeTC['Young'], lifespan['Dying']),
ctrl.Rule(temperature['Low'] & d_speed['St'] & timeTC['Middle'], lifespan['Dying']),
ctrl.Rule(temperature['Low'] & d_speed['St'] & timeTC['Old'], lifespan['Dying']),
ctrl.Rule(temperature['Low'] & d_speed['+'] & timeTC['Young'], lifespan['St']),
ctrl.Rule(temperature['Low'] & d_speed['+'] & timeTC['Middle'], lifespan['St']),
ctrl.Rule(temperature['Low'] & d_speed['+'] & timeTC['Old'], lifespan['St']),
ctrl.Rule(temperature['Low'] & d_speed['++'] & timeTC['Young'], lifespan['Living']),
ctrl.Rule(temperature['Low'] & d_speed['++'] & timeTC['Middle'], lifespan['Living']),
ctrl.Rule(temperature['Low'] & d_speed['++'] & timeTC['Old'], lifespan['Living']),
ctrl.Rule(temperature['Ok'] & d_speed['--'] & timeTC['Young'], lifespan['Dying']),
ctrl.Rule(temperature['Ok'] & d_speed['--'] & timeTC['Middle'], lifespan['Dying']),
ctrl.Rule(temperature['Ok'] & d_speed['--'] & timeTC['Old'], lifespan['Dying']),
ctrl.Rule(temperature['Ok'] & d_speed['-'] & timeTC['Young'], lifespan['Dying']),
ctrl.Rule(temperature['Ok'] & d_speed['-'] & timeTC['Middle'], lifespan['Dying']),
ctrl.Rule(temperature['Ok'] & d_speed['-'] & timeTC['Old'], lifespan['Dying']),
ctrl.Rule(temperature['Ok'] & d_speed['St'] & timeTC['Young'], lifespan['Living']),
ctrl.Rule(temperature['Ok'] & d_speed['St'] & timeTC['Middle'], lifespan['St']),
ctrl.Rule(temperature['Ok'] & d_speed['St'] & timeTC['Old'], lifespan['Dying']),
ctrl.Rule(temperature['Ok'] & d_speed['+'] & timeTC['Young'], lifespan['Living']),
ctrl.Rule(temperature['Ok'] & d_speed['+'] & timeTC['Middle'], lifespan['Living']),
ctrl.Rule(temperature['Ok'] & d_speed['+'] & timeTC['Old'], lifespan['Living']),
ctrl.Rule(temperature['Ok'] & d_speed['++'] & timeTC['Young'], lifespan['Living+']),
ctrl.Rule(temperature['Ok'] & d_speed['++'] & timeTC['Middle'], lifespan['Living+']),
ctrl.Rule(temperature['Ok'] & d_speed['++'] & timeTC['Old'], lifespan['Living+']),
]
SF1_ctrl = ctrl.ControlSystem(SF1_rules)
SF1 = ctrl.ControlSystemSimulation(SF1_ctrl)

input_temperature_fuzz = input_fuzz(temperature, input_temperature)
input_d_speed_fuzz = input_fuzz(d_speed, input_d_speed)
input_time_fuzz = input_fuzz(timeTC, input_time)

IRR = IRR_3var(SF1_rules, input_temperature_fuzz, input_d_speed_fuzz, input_time_fuzz)

declenchement = IRR.min(axis=1)

csq_SF1 = {}  # variable de type dictionnaire
for rule_i in range(0, len(SF1_rules)):  # on initialise à 0 pour chaque règle
csq_SF1[str(SF1_rules[rule_i].consequent)] = 0
for rule_i in range(0, len(SF1_rules)):
csq_i = str(SF1_rules[rule_i].consequent)
csq_SF1[csq_i] = max(csq_SF1[csq_i], declenchement[rule_i])  # Comme dans le TP6

return output_fuzz(csq_SF1, lifespan)

def SF2_compute(input_humidity, input_TC_size, want_to_plot):
pas = 0.01
humidity = ctrl.Antecedent(np.arange(0, 100, pas), 'humidity')  # Humidity (%)'
TC_size = ctrl.Antecedent(np.arange(0, 350, pas), 'TC_size')  # TC size (km)
TC_size_variation = ctrl.Consequent(np.arange(0, 60, pas), 'TC_var')  # Variation of TC size for next day (km)

humidity['--'] = fuzz.trapmf(humidity.universe, [0, 0, 20, 40])  # 20%
humidity['-'] = fuzz.trimf(humidity.universe, [20, 40, 60])  # 40%
humidity['+'] = fuzz.trimf(humidity.universe, [40, 60, 80])  # 60%
humidity['++'] = fuzz.trapmf(humidity.universe, [60, 80, 100, 100])  # 80%

TC_size['Little+'] = fuzz.trapmf(TC_size.universe, [0, 0, 75, 120])
TC_size['Little'] = fuzz.trimf(TC_size.universe, [75, 150, 200])
# TC_size['Normal'] = fuzz.trimf(TC_size.universe, [120, 150, 200])
TC_size['Huge'] = fuzz.trimf(TC_size.universe, [150, 200, 250])
TC_size['Huge+'] = fuzz.trapmf(TC_size.universe, [200, 250, 350, 350])

# 10 13 15 20 25 50

TC_size_variation['Very Low'] = fuzz.trapmf(TC_size_variation.universe, [0, 0, 10, 13])  # 10
TC_size_variation['Low'] = fuzz.trimf(TC_size_variation.universe, [10, 13, 15])  # 13
TC_size_variation['Quite Low'] = fuzz.trimf(TC_size_variation.universe, [13, 15, 20])  # 15
TC_size_variation['Quite High'] = fuzz.trimf(TC_size_variation.universe, [15, 20, 25])  # 20
TC_size_variation['High'] = fuzz.trimf(TC_size_variation.universe, [20, 25, 50])  # 25
TC_size_variation['Very High'] = fuzz.trapmf(TC_size_variation.universe, [25, 50, 50, 50])  # 50

if want_to_plot:
humidity.view()
TC_size.view()
TC_size_variation.view()
SF2_rules = [ctrl.Rule(humidity['--'] & TC_size['Little+'], TC_size_variation['Very Low']),
ctrl.Rule(humidity['--'] & TC_size['Little'], TC_size_variation['Quite Low']),
ctrl.Rule(humidity['--'] & TC_size['Huge'], TC_size_variation['Quite Low']),
ctrl.Rule(humidity['--'] & TC_size['Huge+'], TC_size_variation['Quite Low']),
ctrl.Rule(humidity['-'] & TC_size['Little+'], TC_size_variation['Low']),
ctrl.Rule(humidity['-'] & TC_size['Little'], TC_size_variation['Quite Low']),
ctrl.Rule(humidity['-'] & TC_size['Huge'], TC_size_variation['Quite Low']),
ctrl.Rule(humidity['-'] & TC_size['Huge+'], TC_size_variation['Quite Low']),
ctrl.Rule(humidity['+'] & TC_size['Little+'], TC_size_variation['Quite High']),
ctrl.Rule(humidity['+'] & TC_size['Little'], TC_size_variation['Quite High']),
ctrl.Rule(humidity['+'] & TC_size['Huge'], TC_size_variation['Quite High']),
ctrl.Rule(humidity['+'] & TC_size['Huge+'], TC_size_variation['Quite High']),
ctrl.Rule(humidity['++'] & TC_size['Little+'], TC_size_variation['High']),
ctrl.Rule(humidity['++'] & TC_size['Little'], TC_size_variation['Very High']),
ctrl.Rule(humidity['++'] & TC_size['Huge'], TC_size_variation['Very High']),
ctrl.Rule(humidity['++'] & TC_size['Huge+'], TC_size_variation['Very Low'])]

SF2_ctrl = ctrl.ControlSystem(SF2_rules)
SF2 = ctrl.ControlSystemSimulation(SF2_ctrl)

input_TC_size_fuzz = input_fuzz(TC_size, input_TC_size)
input_humidity_fuzz = input_fuzz(humidity, input_humidity)

IRR = IRR_2var(SF2_rules, input_humidity_fuzz, input_TC_size_fuzz)

declenchement = IRR.min(axis=1)

csq_SF2 = {}  # variable de type dictionnaire
for rule_i in range(0, len(SF2_rules)):  # on initialise à 0 pour chaque règle
csq_SF2[str(SF2_rules[rule_i].consequent)] = 0
for rule_i in range(0, len(SF2_rules)):
csq_i = str(SF2_rules[rule_i].consequent)
csq_SF2[csq_i] = max(csq_SF2[csq_i], declenchement[rule_i])  # Comme dans le TP6

return output_fuzz(csq_SF2, TC_size_variation)  # Dictionnaire des conséquences, et variable de sortie

def fuzz_CAF1(uwind, vwind, x_TC, y_TC):
# Fuzzification des matrices, avec des NFT
# Soit x le nombre net, alors cela revient à l'ajout d'une dimension à la matrice, de valeur : [x+a,x+b,x+c,x+d]
a, b, c, d = [-2, -1, 1, 2] # Augmenter l'envergure augmente la zone des cases dans le graph final

uwind = uwind[:, :, np.newaxis]  # Extension de la matrice en 3D

# Fuzzification uwind vwind
# Créer une dimension pour les coordonnées des nombres trapézoidaux, de forme [a,b,c,d]

uwind = np.dstack((uwind + a, uwind + b, uwind + c, uwind + d))  # La matrice est maintenant fuzzifiée

a, b, c, d = [-2, -1, 1,2]

vwind = vwind[:, :, np.newaxis]

vwind = np.dstack((vwind + a, vwind + b, vwind + c, vwind + d))  # La matrice est maintenant fuzzifiée

dNFT = [-10000, -5000, 5000, 10000] #Faut pas trop augmenter sinon ça deviens incohérent
x_TC = np.array(x_TC + dNFT)
y_TC = np.array(y_TC + dNFT)
return uwind, vwind, x_TC, y_TC

def CAF1_compute(uwind, vwind, x_TC, y_TC):
# Calcul PSI, PSI_TC
# ---------
# PSI : Ce dernier est calculé à partir d'une intégrale, donc d'une somme et d'une division par un nombre net
# U et V sont des IFT de la forme (a,b,c,d)
# Nous pouvons calculer l'intégrale pour a,b,c et d, 1 à 1

# Intégrales de U et V pour déterminer psi, à l'instant t
# Cf.Formule
intx = scipy.integrate.cumtrapz(vwind, lon, axis=1, initial=0)
inty = np.flip(scipy.integrate.cumtrapz(uwind, lat, axis=0, initial=0),axis=2)
#############
# POURQUOI FAUT FLIP INTY??????
#############
# La soustraction nécessite un np.flip
psi = (inty - np.flip(intx,axis=2)) / 2

# ----------
# PSI_TC : Comme au dessus, PSI_TC se base sur une moyenne, donc par le calcul d'additions et d'une division par un nombre net
# Nous pouvons calculer directement PSI_TC pour les 4 termes qui composent notre IFT
# Nous allons le déterminer avec les coordonnées polaires
# L'idée est de longer selon l'axe x, et de calculer PSI_TC sur le cercle de rayon x

res = 3000
psi_TC_xy = np.zeros(np.shape(psi))
psi_TC_r = np.zeros((res, 4))

K = np.zeros(4)
for i in range(0, 4):
univers = np.linspace(x_TC[i], max(lon) + 5000000, res)
# On va au delà de lon, nécessaire car on travaille avec des cercles, faut sortir du carré

# L'Univers est pour définir l'espace de travail de x
# On veut PSI_TC selon R (une liste), et PSI_TC selon (x,y) (un champ 2D)
psi_TC_r[:, i], psi_TC_xy[:, :, i] = psi_TC_compute(x_TC[i], y_TC[i], univers, psi[:, :, i], is_fuzzy=True)

# Soustraction par psi_inf.
psi_inf = np.mean(np.flip(psi, axis=2)[:, :, i])

psi_TC_r[:, i] = [x - psi_inf for x in psi_TC_r[:, i]]
psi_TC_xy[:, :, i] = psi_TC_xy[:, :, i] - psi_inf
K[i] = K_compute(psi_TC_r[:, i], univers)
# Le paramètre K nous sera aussi nécessaire

psi_large = (psi - np.flip(psi_TC_xy, axis=2))

return psi, psi_large, psi_TC_xy, K

x_TC_out = x0  # On prépare la sortie
y_TC_out = y0
# A partir des champs et de la position du TC, nous pouvons prédire sa trajectoire pour dans un temps t_1s
# La vitesse varie dans le champ, nous ferons donc des "petits pas" de temps dts, pour arriver à l'instant final
# Retourne la position finale (x;y)
# Cependant il n'est pas intéressant de prendre les valeurs de vitesse en 1 point : on s'intéresse à une zone
# de r mètres autour du cyclone. On prendra la vitesse moyenne dans cette zone
# _____________
# Divers paramètres pour le calcul :
r_min = 20000
C = 0.7156
nm = 1.2830
a = 1.5

for k in range(0, int(np.round(t_1s)), int(np.round(dts))):

r_map, angle_map = cart_to_polar(x_TC_out, y_TC_out, lon, lat)

u_large_r1 = []
v_large_r1 = []

r_univers = np.arange(r_min, r_TC, 10000)

for r in r_univers:
u_large_r1.append(np.nanmean(array_r1_compute(r, r_map, angle_map, u_large)))
v_large_r1.append(np.nanmean(array_r1_compute(r, r_map, angle_map, v_large)))

# FLOU
# Ces lignes là seront à changer pour le flou
u = np.mean(u_large_r1)
v = np.mean(v_large_r1)
#
TC_lat = np.radians(met_to_deg(y_TC_out))  # Détermination de la latitude du cyclone, en rad
# Calcul du paramètre m. Celui-ci corrige 'l'erreur' lié au fait que nous travaillons sur une sphère, et donc
# que les distances sont courbées comparées à une étude sur un plan 2D
m = (C * (np.cos(TC_lat)) ** (-1)) * (np.tan(0.25 * (np.pi - 2 * TC_lat))) ** nm

# idx_x, idx_y = xy_to_indx_lonlat(x_TC_out, y_TC_out, lon, lat)  # On récupère les index selon [lon,lat] du TC

Cx0 = a * u # - K * ny  # Vélocité selon x au centre du cyclone, en m/s
Cy0 = a * v # + K * nx  # Vélocité selon y au centre du cyclone, en m/s

# On actualise la position
x_TC_out = x_TC_out + Cx0 * dts * m  # m > 0
y_TC_out = y_TC_out + Cy0 * dts * m

return x_TC_out, y_TC_out

def CAF2_compute(uwind, vwind, psi, psi_TC, psi_large, x_TC, y_TC, r_TC, temps_etude, K,SF3_output):
f = np.expand_dims(lat, axis=1) * np.ones(np.shape(psi[:, :, 0]))
f = parametre_coriolis(f)
f = np.dstack(([f, f, f, f]))

delta = lat[1] - lat[0]  # Résolution de nos données
n = scipy.ndimage.laplace(psi_large) / (delta ** 2) + 1 * f

# u_TC = fuzzy_gradient(psi_TC, delta=lat[1] - lat[0], axis=0)
# v_TC = -np.flip(fuzzy_gradient(psi_TC, delta=lon[1] - lon[0], axis=1), axis=2)

v_large = vwind[:, :, :] - np.flip(v_TC, axis=2)  # La négation du nombre flou demande le np.flip()
u_large = uwind[:, :, :] - np.flip(u_TC, axis=2)

# Il vont être utilisés via une soustraction

# On lisse les données
# smooth = 5

# On a tous les paramètres nécessaires au calcul de la vélocité, et donc de la position future du cyclone
# La vélocité dépend des champs n et PSI. Ces derniers n'étant pas uniformes, la vélocité n'est pas la même partout
# Il est donc nécessaire d'étudier la vélocité à petit pas, avec un pas de temps dt
# Soit t_1 (heures) l'instant futur du cyclone
# Nous travaillons en unités SI

t_1s = temps_etude * 60 * 60
dts = (temps_etude / 3) * 60 * 60

# Idée pour le flou : Ici r_TC c'est la taille max ; Mais ça sera un nb flou.
# On trace la moyenne de U,V,N en fonction du rayon
# Pour chaque valeur du rayon on multiplie par le degré de vérité de r_TC
# Pour chaque pas de temps dt, allant de l'instant t0 (maintenant) à l'instant t_1

# Désormais il sera obligatoire de différencier les x en fonction de leur hauteur
h = 1
x_TC = np.append(x_TC, h)
y_TC = np.append(y_TC, h)
# Constitué de plein de IFT. 1ère dim : chaque point.
x_field = np.zeros((1, 5))
y_field = np.zeros((1, 5))
x_field[0] = x_TC
y_field[0] = y_TC
# Au départ x_field n'est composé que d'un point

for t in range(0, t_1s, int(dts)):
# Pour chaque point i flou ; donc pour chaque ligne
x_field_1 = np.zeros((1, 5))
y_field_1 = np.zeros((1, 5))
print("Calcul...", t / dts)
for i in tqdm(range(0, np.size(x_field[:, 0]))):

# On récupère les coordonnées
xi_fuzz = np.transpose(np.expand_dims(x_field[i, :-1], axis=0))  # Vecteur colonne
yi_fuzz = np.transpose(np.expand_dims(y_field[i, :-1], axis=0))

# On met le nombre (singleton) sous l'écriture de Kaufman.
yi_fuzz = yi_fuzz * np.ones((4, 5))
xi_fuzz = xi_fuzz * np.ones((4, 5))

# Pour la hauteur :

xi_fuzz[:, 4] = x_field[i, 4]
yi_fuzz[:, 4] = y_field[i, 4]

# A l'extrémité des IFT, on a forcément h=0
xi_fuzz[3::4, 4] = 0
yi_fuzz[3::4, 4] = 0
xi_fuzz[::4, 4] = 0
yi_fuzz[::4, 4] = 0

# Ajout des alpha_coupes
# Alpha = 0.5 ; La hauteur est de moitié
to_insert = (xi_fuzz[0, :] + xi_fuzz[1, :]) / 2
to_insert1 = (xi_fuzz[2, :] + xi_fuzz[3, :]) / 2

xi_fuzz = np.insert(xi_fuzz, [1], to_insert, axis=0)
xi_fuzz = np.insert(xi_fuzz, [-1], to_insert1, axis=0)

to_insert = (yi_fuzz[0, :] + yi_fuzz[1, :]) / 2
to_insert1 = (yi_fuzz[2, :] + yi_fuzz[3, :]) / 2

yi_fuzz = np.insert(yi_fuzz, [1], to_insert, axis=0)
yi_fuzz = np.insert(yi_fuzz, [-1], to_insert1, axis=0)

# Nous avons xy_fuzz (net), en 1 point, nous calculons son déplacement (flou)
# Comme il s'agit plus précisément d'une addition, alors on procède élément par élément
for m in range(0, np.size(xi_fuzz[:, 0])):
for j in range(0, 4):  # Il ne va pas à 4
xi_fuzz[m, j], yi_fuzz[m, j] = xy_pred(xi_fuzz[m, j], yi_fuzz[m, j], u_large[:, :, j],
v_large[:, :, j],
r_TC, dts, dts/2, K[j],SF3_output)

# Nous avons le déplacement flou
x_field_1 = np.append(x_field_1, xi_fuzz, axis=0)
y_field_1 = np.append(y_field_1, yi_fuzz, axis=0)

# On enlève ceux avec une possibilité trop faible
x_field = x_field_1[x_field_1[:, 4] >= 0.1]
y_field = y_field_1[y_field_1[:, 4] >= 0.1]

print('u_large.mean = ', np.mean(u_large, axis=(0, 1)))
print('v_large.mean = ', np.mean(v_large, axis=(0, 1)))

return x_field, y_field

def IFT_appartenance(a, univers):
a1, a2, a3, a4, h = a

mu = float("nan") * univers
mu[(univers < a1)] = 0
mu[(univers > a1) & (univers < a2)] = h * (univers[(univers > a1) & (univers < a2)] - a1) / (a2 - a1)
mu[(univers > a2) & (univers < a3)] = h
mu[(univers > a3) & (univers < a4)] = h * (a4 - univers[(univers > a3) & (univers < a4)]) / (a4 - a3)
mu[(univers > a4)] = 0

return mu

def xy_fuzz_plot(x, y, lon_plot, lat_plot):
# Soit un point de coordonnées X et Y IFT notation de Kaufmann
# Ressort un champ /meshgrid avec les degrés de possibilités
# Unités SI (mètres)

grid = np.meshgrid(lon_plot, lat_plot)
plot_grid = 0 * grid[0]
print("___________AFFICHAGE_EN_COURS_________")
for i in tqdm(range(0, np.size(x[:, 0]))):  # Pour chaque point flou

mu_x = IFT_appartenance(x[i, :], grid[0])
mu_y = IFT_appartenance(y[i, :], grid[1])
xy_grid = mu_x * mu_y  # FAUT DEFINIR UNE T NORME

# là où il n'y a pas de 0 dans xy_grid on rajoute xy-grid, sinon on met plot_grid (= on change rien)
plot_grid = np.where(mask, np.maximum(plot_grid, xy_grid), plot_grid)

plot_array(plot_grid, lon_plot, lat_plot)

def SF3_compute(input_TC_velocity, input_surface_pressure, input_dsurface_pressure, want_to_plot):
velocity = ctrl.Antecedent(np.arange(0, 100, 0.01), 'Velocity')  # Velocity (m/s)
surface_pressure = ctrl.Antecedent(np.arange(92000, 103000, 10), 'Surface_Pressure')  # Surface pressure (Pa)
dsurface_pressure = ctrl.Antecedent(np.arange(-2000, 2000, 10),
'Surface_Pressure_Variation')  # Surface pressure variation (Pa/h) Surface variation per hour or per day?

velocity_variation = ctrl.Consequent(np.arange(-25, 25, 0.1),
'Velocity_Variation')  # Variation of the velocity next hour (m/s). Next hour / next day  ?

velocity['Little'] = fuzz.trapmf(velocity.universe, [0, 0, 20, 40])
velocity['Intermediate'] = fuzz.trapmf(velocity.universe, [20, 30, 50, 60])
velocity['High'] = fuzz.trapmf(velocity.universe, [40, 60, 100, 100])

surface_pressure['Little'] = fuzz.trapmf(surface_pressure.universe, [92000, 92000, 99000, 101000])
surface_pressure['High'] = fuzz.trapmf(surface_pressure.universe, [99000, 101000, 103000, 103000])

dsurface_pressure['--'] = fuzz.trapmf(dsurface_pressure.universe, [-2000, -2000, -1500, -500])
dsurface_pressure['-'] = fuzz.trimf(dsurface_pressure.universe, [-1500, -500, 0])
dsurface_pressure['0'] = fuzz.trimf(dsurface_pressure.universe, [-500, 0, 500])
dsurface_pressure['+'] = fuzz.trimf(dsurface_pressure.universe, [0, 500, 1500])
dsurface_pressure['++'] = fuzz.trapmf(dsurface_pressure.universe, [500, 1500, 2000, 2000])

velocity_variation['--'] = fuzz.trimf(velocity_variation.universe, [-25, -25, -10])
velocity_variation['-'] = fuzz.trapmf(velocity_variation.universe, [-25, -10, -5, 0])
velocity_variation['St'] = fuzz.trimf(velocity_variation.universe, [-5, 0, 5])
velocity_variation['+'] = fuzz.trapmf(velocity_variation.universe, [0, 5, 10, 25])
velocity_variation['++'] = fuzz.trimf(velocity_variation.universe, [10, 25, 25])

if want_to_plot:
velocity.view()
surface_pressure.view()
dsurface_pressure.view()
velocity_variation.view()

SF3_rules = [
ctrl.Rule(surface_pressure['Little'] & velocity['Little'] & dsurface_pressure['--'], velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['Intermediate'] & dsurface_pressure['-'],
velocity_variation['+']),
ctrl.Rule(surface_pressure['Little'] & velocity['High'] & dsurface_pressure['0'], velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['Little'] & dsurface_pressure['+'], velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['Intermediate'] & dsurface_pressure['++'],
velocity_variation['-']),
ctrl.Rule(surface_pressure['Little'] & velocity['High'] & dsurface_pressure['--'], velocity_variation['++']),
ctrl.Rule(surface_pressure['Little'] & velocity['Little'] & dsurface_pressure['-'], velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['Intermediate'] & dsurface_pressure['0'],
velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['High'] & dsurface_pressure['+'], velocity_variation['-']),
ctrl.Rule(surface_pressure['Little'] & velocity['Little'] & dsurface_pressure['++'], velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['Intermediate'] & dsurface_pressure['--'],
velocity_variation['+']),
ctrl.Rule(surface_pressure['Little'] & velocity['High'] & dsurface_pressure['-'], velocity_variation['++']),
ctrl.Rule(surface_pressure['Little'] & velocity['Little'] & dsurface_pressure['0'], velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['Intermediate'] & dsurface_pressure['+'],
velocity_variation['-']),
ctrl.Rule(surface_pressure['Little'] & velocity['High'] & dsurface_pressure['++'], velocity_variation['--']),
ctrl.Rule(surface_pressure['High'] & velocity['Little'] & dsurface_pressure['--'], velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['Intermediate'] & dsurface_pressure['-'],
velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['High'] & dsurface_pressure['0'], velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['Little'] & dsurface_pressure['+'], velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['Intermediate'] & dsurface_pressure['++'],
velocity_variation['-']),
ctrl.Rule(surface_pressure['High'] & velocity['High'] & dsurface_pressure['--'], velocity_variation['+']),
ctrl.Rule(surface_pressure['High'] & velocity['Little'] & dsurface_pressure['-'], velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['Intermediate'] & dsurface_pressure['0'],
velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['High'] & dsurface_pressure['+'], velocity_variation['-']),
ctrl.Rule(surface_pressure['High'] & velocity['Little'] & dsurface_pressure['++'], velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['Intermediate'] & dsurface_pressure['--'],
velocity_variation['+']),
ctrl.Rule(surface_pressure['High'] & velocity['High'] & dsurface_pressure['-'], velocity_variation['+']),
ctrl.Rule(surface_pressure['High'] & velocity['Little'] & dsurface_pressure['0'], velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['Intermediate'] & dsurface_pressure['+'],
velocity_variation['-']),
ctrl.Rule(surface_pressure['High'] & velocity['High'] & dsurface_pressure['++'], velocity_variation['-'])]

SF3_ctrl = ctrl.ControlSystem(SF3_rules)
SF3 = ctrl.ControlSystemSimulation(SF3_ctrl)

input_velocity_fuzz = input_fuzz(velocity, input_TC_velocity)
input_surface_pressure_fuzz = input_fuzz(surface_pressure, input_surface_pressure)
input_dsurface_pressure_fuzz = input_fuzz(dsurface_pressure, input_dsurface_pressure)

IRR = IRR_3var(SF3_rules, input_surface_pressure_fuzz, input_velocity_fuzz, input_dsurface_pressure_fuzz)

declenchement = IRR.min(axis=1)

csq_SF3 = {}  # variable de type dictionnaire
for i in range(0, len(SF3_rules)):  # on initialise à 0 pour chaque règle
csq_SF3[str(SF3_rules[i].consequent)] = 0
for i in range(0, len(SF3_rules)):
csq_i = str(SF3_rules[i].consequent)
csq_SF3[csq_i] = max(csq_SF3[csq_i], declenchement[i])  # Comme dans le TP6

return output_fuzz(csq_SF3, velocity_variation)

def estimation2(t0, t1, dt, uwind, vwind, study_grid):
# Estimation flou de la trajectoire (sortie CAF2), de t0 à t1

print("_________________________SF2_______________________")

SF2_humidity = np.where(study_grid, humidity[t0, :, :], float('nan'))
SF2_humidity = np.nanmean(SF2_humidity)

SF2_taille_TC0 = taille_TC(t0, study_grid)  # Rayon, en km

dtaille = SF2_compute(SF2_humidity, SF2_taille_TC0,
want_to_plot=False)  # Prediction pour dans 24 heures, en pourcentage

dtaille[:, 0] = (t1 - t0) * (dtaille[:, 0] / 24)  # On le ramène à une prévision au temps t1
dtaille[:, 0] = (dtaille[:, 0] * 0.01)

dtaille[:, 0] = (dt) * dtaille[:, 0] / 24  # Pour les calculs de vitesse nous voudrons l'évolution pour chaque dt
SF2_taille_TC1 = dtaille
SF2_taille_TC1[:, 0] = SF2_taille_TC0 * (1+dtaille[:, 0])

print("Taille Actuelle du cyclone : ", np.rint(SF2_taille_TC0), "Humidité :", np.rint(SF2_humidity))
print("_________________________SORTIE SF2 - Future Taille du Cyclone à dt___")
print("Prédiction :", SF2_taille_TC1, 'Réel :', np.rint(taille_TC(t0 + dt, study_grid)))

print("_____________________________CAF1____________")
uwind_fuzz, vwind_fuzz, x_TC_fuzz, y_TC_fuzz = fuzz_CAF1(uwind[t0, :, :], vwind[t0, :, :], x_TC[t0], y_TC[t0])
print('X =', x_TC_fuzz, 'Y = ', y_TC_fuzz)
print('uwind.mean = ', np.rint(np.mean(uwind_fuzz, axis=(0, 1))))
print('vwind.mean = ',np.rint(np.mean(vwind_fuzz, axis=(0, 1))))
print('____')
psi, psi_large, psi_TC, K = CAF1_compute(uwind_fuzz, vwind_fuzz, x_TC_fuzz, y_TC_fuzz)
print('Psi.mean =', np.rint(np.mean(psi, axis=(0, 1))))
print('Psi_large.mean = ', np.rint(np.mean(psi_large, axis=(0, 1))))
print('Psi_TC.mean = ', np.rint(np.mean(psi_TC, axis=(0, 1))))
print('K', np.rint(K))
print('_______________CAF1 - DONE___________________')

print("_______________SF3___________________")
#SF3_speed_TC = vitesse_TC(0,t0,calcul_pth())
#SF3_pressure_TC = pressure[t0,idx_y_TC[t0],idx_x_TC[t0]]
#SF3_dpressure_TC = pressure[t0,idx_y_TC[t0],idx_x_TC[t0]]
SF3_output = 0#SF3_compute(SF3_speedTC,SF3_pressure_TC,SF3_dpressure_TC)
print("_____________________SF3 - DONE______________")

x, y = CAF2_compute(uwind_fuzz, vwind_fuzz, psi, psi_TC, psi_large, x_TC_fuzz, y_TC_fuzz, 200000,
t1 - t0,
K,SF3_output)

print('Il y a ', np.size(x[:, 0]), "points !")
print('________________CAF2 - DONE________________________')

# plot_array(psi[:, :, 0], lon, lat)

# for i in range(0, np.size(x[:, 0])):
#    if x[i, 4] == 1:
# plt.plot(met_to_deg(x[i, 0:4]), met_to_deg(y[i, 0:4]), ms=1, marker='x', markeredgecolor='red')
plt.plot(met_to_deg(x_TC[t0]), met_to_deg(y_TC[t0]), ms=10, marker='o', markeredgecolor='red')
plt.plot(met_to_deg(x_TC[t1]), met_to_deg(y_TC[t1]), ms=10, marker='x', markeredgecolor='red')
# plt.show()

# On augmente la résolution pour l'affichage des degrés de possibilités, sinon on n'est pas assez précis
lon_fuzz_plot = np.arange(np.min(lon), np.max(lon), ((lon[1]) - (lon[0])) / 6)
lat_fuzz_plot = np.arange(np.max(lat), np.min(lat), ((lat[1]) - (lat[0])) / 6)
xy_fuzz_plot(x, y, lon_fuzz_plot, lat_fuzz_plot)

plt.show()

def SF4_compute(input_actu, input_ant, want_to_plot):
tendance_actuelle = ctrl.Antecedent(np.arange(-90, 90, 0.1), 'Angle calcul')
tendance_antecedents = ctrl.Antecedent(np.arange(-90, 90, 0.1), 'Angle évalué par la forme')

angle = ctrl.Consequent(np.arange(-90, 90, 0.1), 'Prédiction courbature')

tendance_actuelle['Courbé Gauche Fort'] = fuzz.trapmf(tendance_actuelle.universe, [-90, -90, -60, -45])
tendance_actuelle['Courbé Gauche'] = fuzz.trapmf(tendance_actuelle.universe, [-60, -45, -10, 0])
tendance_actuelle['Linéaire'] = fuzz.trapmf(tendance_actuelle.universe, [-10, 0, 0, 10])
tendance_actuelle['Courbé Droite'] = fuzz.trapmf(tendance_actuelle.universe, [0, 10, 45, 60])
tendance_actuelle['Courbé Droite Fort'] = fuzz.trapmf(tendance_actuelle.universe, [45, 60, 90, 90])

tendance_antecedents['Courbé Gauche Fort'] = fuzz.trapmf(tendance_antecedents.universe, [-90, -90, -60, -45])
tendance_antecedents['Courbé Gauche'] = fuzz.trapmf(tendance_antecedents.universe, [-60, -45, -20, 0])
tendance_antecedents['Linéaire'] = fuzz.trapmf(tendance_antecedents.universe, [-20, 0, 5, 20])
tendance_antecedents['Courbé Droite'] = fuzz.trapmf(tendance_antecedents.universe, [5, 20, 45, 60])
tendance_antecedents['Courbé Droite Fort'] = fuzz.trapmf(tendance_antecedents.universe, [45, 60, 90, 90])

angle['Courbé Gauche Fort'] = fuzz.trapmf(angle.universe, [-90, -90, -60, -45])
angle['Courbé Gauche'] = fuzz.trapmf(angle.universe, [-60, -45, -20, 0])
angle['Linéaire'] = fuzz.trapmf(angle.universe, [-20, -5, 5, 20])
angle['Courbé Droite'] = fuzz.trapmf(angle.universe, [5, 20, 45, 60])
angle['Courbé Droite Fort'] = fuzz.trapmf(angle.universe, [45, 60, 90, 90])

if want_to_plot:
tendance_actuelle.view()
tendance_antecedents.view()
angle.view()

SF4_rules = [ctrl.Rule(tendance_actuelle['Courbé Gauche Fort'] & tendance_antecedents['Courbé Gauche Fort'],
angle['Courbé Gauche Fort']),
ctrl.Rule(tendance_actuelle['Courbé Gauche Fort'] & tendance_antecedents['Courbé Gauche'],
angle['Courbé Gauche Fort']),
ctrl.Rule(tendance_actuelle['Courbé Gauche Fort'] & tendance_antecedents['Linéaire'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Courbé Gauche Fort'] & tendance_antecedents['Courbé Droite'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Courbé Gauche Fort'] & tendance_antecedents['Courbé Droite Fort'],
angle['Linéaire']),
ctrl.Rule(tendance_actuelle['Courbé Gauche'] & tendance_antecedents['Courbé Gauche Fort'],
angle['Courbé Gauche Fort']),
ctrl.Rule(tendance_actuelle['Courbé Gauche'] & tendance_antecedents['Courbé Gauche'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Courbé Gauche'] & tendance_antecedents['Linéaire'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Courbé Gauche'] & tendance_antecedents['Courbé Droite'],
angle['Linéaire']),
ctrl.Rule(tendance_actuelle['Courbé Gauche'] & tendance_antecedents['Courbé Droite Fort'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Linéaire'] & tendance_antecedents['Courbé Gauche Fort'],
angle['Courbé Gauche Fort']),
ctrl.Rule(tendance_actuelle['Linéaire'] & tendance_antecedents['Courbé Gauche'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Linéaire'] & tendance_antecedents['Linéaire'], angle['Linéaire']),
ctrl.Rule(tendance_actuelle['Linéaire'] & tendance_antecedents['Courbé Droite'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Linéaire'] & tendance_antecedents['Courbé Droite Fort'],
angle['Courbé Droite Fort']),
ctrl.Rule(tendance_actuelle['Courbé Droite'] & tendance_antecedents['Courbé Gauche Fort'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Courbé Droite'] & tendance_antecedents['Courbé Gauche'],
angle['Linéaire']),
ctrl.Rule(tendance_actuelle['Courbé Droite'] & tendance_antecedents['Linéaire'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Courbé Droite'] & tendance_antecedents['Courbé Droite'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Courbé Droite'] & tendance_antecedents['Courbé Droite Fort'],
angle['Courbé Droite Fort']),
ctrl.Rule(tendance_actuelle['Courbé Droite Fort'] & tendance_antecedents['Courbé Gauche Fort'],
angle['Linéaire']),
ctrl.Rule(tendance_actuelle['Courbé Droite Fort'] & tendance_antecedents['Courbé Gauche'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Courbé Droite Fort'] & tendance_antecedents['Linéaire'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Courbé Droite Fort'] & tendance_antecedents['Courbé Droite'],
angle['Courbé Droite Fort']),
ctrl.Rule(tendance_actuelle['Courbé Droite Fort'] & tendance_antecedents['Courbé Droite Fort'],
angle['Courbé Droite Fort'])]

SF4_ctrl = ctrl.ControlSystem(SF4_rules)
SF4 = ctrl.ControlSystemSimulation(SF4_ctrl)

input_tendance_actuelle_fuzz = input_fuzz(tendance_actuelle, input_actu)
input_tendance_antecedents_fuzz = input_fuzz(tendance_antecedents, input_ant)

#print('Tendance calculée fuzz')
#print(input_tendance_actuelle_fuzz)
#print('Tendance évaluée fuzz')
#print(input_tendance_antecedents_fuzz)

IRR = IRR_2var(SF4_rules, input_tendance_actuelle_fuzz, input_tendance_antecedents_fuzz)
#print(IRR)

declenchement = IRR.min(axis=1)
#print(declenchement)
csq_SF4 = {}  # variable de type dictionnaire
for i in range(0, len(SF4_rules)):  # on initialise à 0 pour chaque règle
csq_SF4[str(SF4_rules[i].consequent)] = 0
for i in range(0, len(SF4_rules)):
csq_i = str(SF4_rules[i].consequent)
csq_SF4[csq_i] = max(csq_SF4[csq_i], declenchement[i])  # Comme dans le TP6

# print("Conséquences SF4 :", csq_SF4)
return (output_fuzz(csq_SF4, angle))

# Différence d'angle nécessaire au calcul de SF4, entre 2 points

# ________

def SF4_angle(t,dt):
return get_angle(x_TC[t],y_TC[t],x_TC[t+dt],y_TC[t+dt]) - get_angle(x_TC[t+dt],y_TC[t+dt],x_TC[t+2*dt],y_TC[t+2*dt])

def SF4_conclusion(angle_TC,angle_tendance,t):
# Prediction pour 12 heures

prediction = np.zeros(np.shape(uwind[0, :, :]))

SF4_conc = SF4_compute(angle_TC, angle_tendance,want_to_plot=False)

for lignes in range(np.shape(SF4_conc)[0]):

a = get_angle(x_TC[t],y_TC[t],x_TC[t-12],y_TC[t-12]) + SF4_conc[lignes][0]

for r in np.arange(0, 1000000, 20000):
xr, yr = polar_to_cart(r, a)

x = xr + x_TC[t]
y = - yr + y_TC[t]

idx_lon, idx_lat = xy_to_indx_lonlat(x,y, lon, lat)

prediction[idx_lat, idx_lon] = SF4_conc[lignes][1]

plot_array(prediction[:, :], lon, lat)
plt.plot(met_to_deg(x_TC[t]),met_to_deg(y_TC[t]),ms=5, marker='o', markeredgecolor='red')
plt.show()

return prediction,SF4_conc
```