Fonctions Floues

 avatar
unknown
plain_text
2 years ago
42 kB
9
Indexable
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


def xy_pred(x0, y0, u_large, v_large, n_grad_x, n_grad_y, r_TC, t_1s, dts, K,SF3_ouput):
    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 = []
        n_grad_y_r1 = []
        n_grad_x_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)))
            n_grad_y_r1.append(np.nanmean(array_r1_compute(r, r_map, angle_map, n_grad_y)))
            n_grad_x_r1.append(np.nanmean(array_r1_compute(r, r_map, angle_map, n_grad_x)))

        # FLOU
        # Ces lignes là seront à changer pour le flou
        u = np.mean(u_large_r1)
        v = np.mean(v_large_r1)
        ny = np.mean(n_grad_y)
        nx = np.mean(n_grad_x)
        #
        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)

    u_TC = np.gradient(psi_TC, lat, axis=0)
    v_TC = -np.flip(np.gradient(psi_TC, lon, 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)

    n_grad_y = np.array(np.gradient(n, lat, axis=0))
    n_grad_x = np.array(np.gradient(n, lon, axis=1))

    # Il vont être utilisés via une soustraction
    n_grad_y = np.flip(n_grad_y, axis=2)

    # On lisse les données
    # smooth = 5

    # n_grad_y = smooth_array(n_grad_y, lon, smooth, axis=1)
    # n_grad_x = smooth_array(n_grad_x, lat, smooth, axis=0)

    # 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],
                                                           n_grad_x[:, :, j], n_grad_y[:, :, 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)))
    print('n_grad_x.mean = ', np.mean(n_grad_x, axis=(0, 1)))
    print('n_grad_y.mean = ', np.mean(np.flip(n_grad_y, axis=2), 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)
        mask = np.invert(xy_grid == 0)
        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
Editor is loading...