# -*- coding:latin-1 -*- pour les caractères accentués avec certains éditeurs

# Méthodes d'estimation
# SIE 3e semestre, GC 5e semestre

print('\n','---------------------------------------------------------')
print(' Gaz parfait mesuré dans 5 états')
print(' compensation conditionnelle avec paramètres = cas combiné')
print(' ---------------------------------------------------------')

"""
On peut poser les conditions de diverses manières.
Ici on choisit l'expression simple : wi = Pi*Vi - cst*Ti,
où le paramètre "cst" correspond à la constante du gaz n*R.

Les données incluent la correction de la faute détectée lors de l'exercice 5.

TOPO-ESO / 12.12.23
"""

import numpy as np
np.set_printoptions(precision=2)

"""
-------------------------------------------------------------------------
Gaz parfait - données 2022

1. Mesures

[mb]        [ml]        [K]
"""

P1 = 1002 ; V1 =  993 ; T1 = 298 ;
P2 = 1994 ; V2 = 1504 ; T2 = 901 ;
# P3 = 2210 ; V3 =  703 ; T3 = 450 ; # avant correction
P3 = 2120 ; V3 =  703 ; T3 = 450 ; # après correction
P4 = 1503 ; V4 = 1482 ; T4 = 665 ;
P5 =  751 ; V5 = 2670 ; T5 = 602 ;



# 2. Ecarts-types (a priori)

print('\n','modèle pour les écarts-types')
print("1 = écart-type (absolu) constant pour chaque type d'observation")
print("2 = écart-type relatif constant pour chaque type d'observation")
cas = int(input(' indice du modèle stochastique choisi: '))
#cas =1 

if cas==1:
    sigmaP = 5.0   # [mb]
    sigmaV = 3.0   # [ml]
    sigmaT = 2.0   # [K]
else:
    sigrelP = 0.4   # erreur relative en % [=mb/100mb]
    sigrelV = 0.3   # erreur relative en % [=ml/100ml]
    sigrelT = 0.5   # erreur relative en % [=K/100K]
    
"""
-------------------------------------------------------------------------
formation du vecteur des observations

d'éliminer ou d'ajouter un état.
"""

obs = np.array([P1, V1, T1, P2, V2, T2, P3, V3, T3, P4, V4, T4, P5, V5, T5])
obs = obs.reshape((-1, 1))

n = len(obs)
m = int(n/3) # nombre de conditions (1 pour chaque état)

# nombre de paramètres, valeurs approchées
u = 1
cstbul = obs[0,0]*obs[1,0]/obs[2,0]

# SURDETERMINATION
r = m-u

"""
-------------------------------------------------------------------------
MODELE FONCTIONNEL

conditions et écarts de fermeture w
construction des matrices des dérivées partielles A et B
linéarisation des conditions en fonction de a
linéarisation des conditions en fonction de P1, V1, T1, .... , Pm, Vm, Tm
"""

# réservation d'espaces contigus en mémoire
w = np.zeros((m,1))
A = np.zeros((m,u))
B = np.zeros((m,n))


# i = indice de la condition = ligne de w, de A et de B
# j = indice de l'observation = colonne de B

for i in range(0,m):

	j = 3*i # indice de la pression de l'état i (= condition)

	w[i,0] = obs[j,0]*obs[j+1,0] - cstbul*obs[j+2,0]
	A[i,0] = -obs[j+2,0]
	B[i,j:j+3] = [obs[j+1,0], obs[j,0], -cstbul]

print('\n','w =\n',w)
print('\n','A =\n',A)
print('\n','B =\n',B)

"""
-------------------------------------------------------------------------
MODELE STOCHASTIQUE

comme pour la compensation conditionnelle ou paramétrique
"""

if cas==1:

    # matrice de covariance des observations pour un état

    sigma = np.array([sigmaP, sigmaV, sigmaT])
    K33 = np.diag(sigma*sigma)

    # Pour plusieurs états non corrélés

    Kll = np.zeros((m*3,m*3))

    for i in range(0,m):
        j = 3*i
        Kll[j:j+3,j:j+3]=K33

else:

    # vecteur des écarts-types pour tous les états / toutes les observations

    sigrel = np.array([sigrelP, sigrelV, sigrelT])
    sigrel = sigrel.reshape((-1, 1))
    sigma = np.zeros((3*m,1))

    for i in range(0,m):
        j = 3*i
        sigma[j:j+3] = sigrel*obs[j:j+3]/100

    Kll = np.diagflat(sigma*sigma)

print('\n','Kll =\n',Kll)

"""
-------------------------------------------------------------------------
COMPENSATION
 
Selon le énoncé, on peut calculer dx comme suit.
dx = -inv(A'*(inv(B*Kll*B')*A) * A'inv(B*Kll*B')*w
Toutefois l'analyse requiert d'autres résultats: résidus, covariances, ...

Par propagation on établit la covariance des paramètres.
Kxx = inv(A'*(inv(B*Kll*B')*A)
Par égard pour le processeur, on garde certains résultats intermédiaires.
Par exemple, et selon l'algorithme dans l'énoncé on peut procéder comme suit.

Pstar = np.linalg.inv(B@Kll@B.T)
Kxx = np.linalg.inv(A.T@Pstar@A)
dx = -Kxx @ A@*Pstar@w

Mieux: on peut décomposer en 2 compensations successives.
D'abord on définit des matrices auxiliaires.
"""

Kstar = B@Kll@B.T
S = Kll@B.T@np.linalg.inv(Kstar)

# Ensuite on définit un résidu provisoire.

vstar = S@w

# C'est la solution d'une compensation conditionnelle classique (formule 3.10).
# vstar = Kll*@B.T@np.linalg.inv(B@Kll@B.T)@w

# On définit encore une matrice auxiliaire.

Astar = -S@A

"""
La formule (6.2) du polycopié devient celle d'une compensation
paramétrique classique (4.13), avec vbulle = vstar et A = Astar.
dx = inv(Astar'*P*Astar) * Astar'*P*vstar
En vue de la suite, on préfère calculer comme suit.
"""

P = np.linalg.inv(Kll)
Kxx = np.linalg.inv(Astar.T@P@Astar)
dx = Kxx @ Astar.T@P@vstar
print('\n','dx =',dx)

# Il faut envisager une itération pour assurer la convergence.

cstcomp = cstbul + dx
print('\n','constante du gaz =',cstcomp)

# On extrait aussi l'écart-type de la constante.

sigmacst = np.sqrt(Kxx)
print('\n','sigma de la constante =',sigmacst)
"""
-------------------------------------------------------------------------
Autres résultats: résidus compensés et observations compensées

Dans la littérature, on trouve une formule très indigeste pour vcomp en
fonction des éléments originaux Kll, A, B et w. Grâce à la décomposition
opérée ci-dessus, les résidus résultent de la 2e étape, soit la compensation
paramétrique avec vbulle = vstar et A = Astar (formule 4.6). Avec les
substitutions nécessaires, on peut appliquer la formule 4.16.
"""

vcomp = vstar - Astar@dx
print('\n','vcomp =\n',vcomp)

# Le contrôle est facile: on doit obtenir les mêmes résidus compensés que
# lors de l'exercice 6 (compensation conditionnelle).

# L'expression des observations compensées est inchangée.

lcomp = obs - vcomp
print('\n','lcomp =\n',lcomp)

"""
-------------------------------------------------------------------------
Contrôle
Les écarts de fermeture calculés avec les valeurs compensées devraient
être nuls. Généralement, ils ne le sont pas pour des raisons numériques,
mais leur valeur doit être négligeable compte tenu de la nature du
problème.
"""

wcomp = np.zeros((m,1))
for i in range(0,m):

    j = 3*i
    wcomp[i,0] = lcomp[j]*lcomp[j+1] - cstcomp*lcomp[j+2]

print('\n','wcomp =\n',wcomp)

"""
-------------------------------------------------------------------------
Description et interprétation des résultats

matrices de covariance des résidus et des observations compensés
Les matrices auxiliaires d'avant sont bien pratiques!

Par rapport à l'exercice 6, l'expression de la matrice de covariance des
résidus change. On l'obtient par propagation de variance de vcomp.
La formule suivante, donnée ici sans démonstration, fait ressortir les
contributions des observations et des paramètres.
S@Kstar@S.T correspond à une compensation conditionnelle classique.
En ajoutant de la flexibilité dans la compensation, les paramètres
du modèle combiné réduisent cette covariance de S@A@Kxx@A.T@S.T.
"""

Kvcomp = S@(Kstar - A@Kxx@A.T)@S.T

# L'expression de la matrice de covariance des observations compensées est
# inchangée.

Klcomp = Kll - Kvcomp

# Telles quelles, l'interprétation de ces matrices est difficile,  
# mais comme toujours pour des matrices de covariance, on peut calculer 
# les écarts-types des valeurs compensées, ainsi que leurs corrélations.

from covmat2cormat import covmat2cormat

[sigmavcomp,Rvcomp] = covmat2cormat(Kvcomp)
sigmavcomp = sigmavcomp.reshape((-1, 1))
print('\n','sigmavcomp =\n',sigmavcomp)
print('\n','Rvcomp =\n',Rvcomp)

# Notamment, on remarque que les résidus compensés sont entièrement
# corrélées au sein d'un même état. Par exemple pour le 2e état
# (lignes et colonnes 3 à 6 de Rvcomp),
# la sous-matrice vaut [+1 +1 -1; +1 +1 -1; -1 -1 +1].
# Si P2*V2/T2 doit augmenter pour s'adapter aux conditions, le changement
# se répartit: P2 et V2 augmentent et T2 diminue dans la même proportion.

[sigmalcomp,Rlcomp] = covmat2cormat(Klcomp)
sigmalcomp = sigmalcomp.reshape((-1, 1))
print('\n','sigmalcomp =\n',sigmalcomp)

# Les écarts-types des observations compensées sont plus petits que les 
# valeurs a priori attribuées aux observations.

print('\n','Rlcomp =\n',Rlcomp)

# Les corrélations ne sont pas triviales. Tous les états sont liés par les
# conditions et toutes les observations compensées sont corrélées.

"""
-------------------------------------------------------------------------
Estimation de la précision
Vu que les observations sont de nature différente, on a renoncé à
introduire des cofacteurs.
Implicitement, on a posé sigma0=1 [sans dimension], donc Kll=Qll et
P=np.linalg.inv(Kll). Sous cette forme, le calcul de l'écart-type a posteriori
donne un nombre [sans dimension] qui exprime le rapport global des
écarts-types a posteriori / a priori pour tous les types de mesures.
"""

P = np.linalg.inv(Kll)
quoglo = np.sqrt((np.transpose(vcomp) @ P @ vcomp)/r)
print('\n','Quotient Global =',quoglo)

"""
-------------------------------------------------------------------------
Analyse complémentaire
Le code suivant fait appel à la théorie de la fiabilité (chapitre 5),
mais on aurait pu l'inclure tel quel pour le modèle conditionnel
et pour le modèle paramétrique.

On peut analyser les résidus de manière plus fine en estimant la
précision de chaque type de mesure séparément.

La part de redondance de chaque observation indique à quel point elle est
contrôlée par les autres observations. Elle est comprise entre 0 (sans
contrôle, donc indispensable) et 1 (entièrement contrôlée, donc inutile).
"""

z = np.diagonal(Kvcomp @ P)
z = z.reshape((-1, 1))
print('\n','z =\n',z)

# La somme des parts de redondance correspond à la redondance globale
# (surdétermination). Autrement dit, on doit avoir: np.sum(z)=r.

print('\n','sum(z) =',np.sum(z))

# On sépare les résidus et les parts de redondance selon le type
#d'observation. On forme des vecteurs distincts pour P, V et T.

vP = np.zeros((m,1))
vV = np.zeros((m,1))
vT = np.zeros((m,1))

KllP = np.zeros((m,m))
KllV = np.zeros((m,m))
KllT = np.zeros((m,m))

zP = np.zeros((m,1))
zV = np.zeros((m,1))
zT = np.zeros((m,1))

for i in range(0,m):
    j = i*3
     
    vP[i,0] = vcomp[j,0]
    vV[i,0] = vcomp[j+1,0]
    vT[i,0] = vcomp[j+2,0]

    KllP[i,i] = Kll[j,j]
    KllV[i,i] = Kll[j+1,j+1]
    KllT[i,i] = Kll[j+2,j+2]
    
    zP[i,0] = z[j]
    zV[i,0] = z[j+1]
    zT[i,0] = z[j+2]
 
print('\n','vP =\n',vP)
print(' zP =\n',zP)

print('\n','vV =\n',vV)
print(' zV =\n',zV)

print('\n','vT =\n',vT)
print(' zT =\n',zT)

PP = np.linalg.inv(KllP)
rP = np.sum(zP)
quoP = np.sqrt((vP.T@PP@vP) / rP)

PV = np.linalg.inv(KllV)
rV = np.sum(zV)
quoV = np.sqrt((vV.T@PV@vV) / rV)

PT = np.linalg.inv(KllT)
rT = np.sum(zT)
quoT = np.sqrt((vT.T@PT@vT) / rT)

print('\n','Quotient P = ',quoP)
print(' Quotient V = ',quoV)
print(' Quotient T = ',quoT)

"""
On aime bien que les différents quotients prennent des valeurs proches.
Sinon, on peut envisager de modifier les écarts-types a priori, car la
précision d'un type d'observation peut être mal estimée par rapport aux
autres.
"""