sejeteralo/eau.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Jun 18 2022
Backend for water price computation
@author: alex delga
"""

import xlrd
import xlwt
import numpy as np
import matplotlib.pyplot as plt # main ploting module
from matplotlib.widgets import Slider
import pandas as pd
import seaborn as sns
sns.set_theme()

class NewModel():
    def __init__(self,pmax=20,vmax=2100): # Set the static parameters
        self.pmax=pmax #maximum price per m3 acceptable
        self.vmax=vmax #consumed volume in m3 for which such price is reached. All inhabitants must be under. 
        
        self.nbpts=1023
        self.tt=np.linspace(0,1-1e-6,self.nbpts) # the -1e-6 is to avoid pb at vv=vmax later
        self.vv=np.zeros((2*self.nbpts,))#must be reset parametrically using tt everytime a,b,c,d,e,vinf are changed. 
        
        self.abopArray = np.linspace(0,200,201) # potential values for abop
        self.abosArray = np.linspace(0,200,201) # potential values for abos
        self.recettesArray = np.linspace (50000,100000,51) # potential values for recettes
        self.vinfArray = np.linspace(0,vmax,vmax+1) # potential values for volume at inflexion point 
        self.aArray = np.linspace(0,1,51) # potential values for a
        self.bArray = np.linspace(0,1,51) # potential values for b
        self.cArray = np.linspace(0,1,51) # potential values for c
        self.dArray = np.linspace(0,1,51) # potential values for d
        self.eArray = np.linspace(0,1,51) # potential values for e
        
        #loading the inhabitant data
        book = xlrd.open_workbook('Eau2018.xls')
        sheet = book.sheet_by_name('CALCULS')
        self.nbHab=363
        self.volume = np.array([sheet.cell_value(r,4) for r in range(1,self.nbHab+1)])#get the volume used in 2018 in m3
        self.vTot=np.sum(self.volume)
        self.volume[self.volume==0]=1e-5 #little trick to avoid divide by 0
        self.status = np.array([sheet.cell_value(r,3) for r in range(1,self.nbHab+1)])#get the RS, RP, PRO status
        self.prix2018 = np.array([sheet.cell_value(r,33) for r in range(1,self.nbHab+1)])#get the price paid in 2018 in EUR
        self.prix2018_m3 = self.prix2018 / self.volume
        

    def mainFunction(self):
        '''
        main function, the one called in the notebook, and that update computations depending on the slider values.
        '''
        # initial values
        self.abop= self.abopArray[100]
        self.abos= self.abosArray[100]
        self.recettes= self.recettesArray[25]
        self.vinf=self.vinfArray[int(self.vmax/2)]
        self.a=self.aArray[25] 
        self.b=self.bArray[25] 
        self.c=self.cArray[25] 
        self.d=self.dArray[25] 
        self.e=self.eArray[25] 
        
        # initial plot
        self.updateComputation()
        
        self.fig = plt.figure(figsize=(9,6))
        self.ax1  = plt.subplot(211)
        self.ax2  = plt.subplot(212,sharex=self.ax1)
        #self.ax1.set_xlabel(r'volume ($m^3$) ')
        self.ax2.set_xlabel(r'volume ($m^3$) ')
        self.ax1.set_ylabel('Facture totale (€)')
        self.ax2.set_ylabel('Prix au m3 (€)')
        
        #default values of ax limits
        self.axx=(min(self.vv),max(self.vv))
        self.ax1y=(0,max(self.prixp))
        self.ax2y=(0,self.pmax)
        
        self.updatePlot()
    
        # SLIDERS
        axcolor   = 'lightgoldenrodyellow'
        abopAxis = plt.axes([0.7, 0.70, 0.2, 0.03], facecolor=axcolor)
        self.abopSlider = Slider(abopAxis, r'Abo. P/PRO (€)', np.min(self.abopArray), np.max(self.abopArray), self.abop)
        self.abopSlider.on_changed(self.updateParam)
        
        abosAxis = plt.axes([0.7, 0.65, 0.2, 0.03], facecolor=axcolor)
        self.abosSlider = Slider(abosAxis, r'Abo. S (€)', np.min(self.abosArray), np.max(self.abosArray), self.abos)
        self.abosSlider.on_changed(self.updateParam)        
        
        recettesAxis = plt.axes([0.7, 0.60, 0.2, 0.03], facecolor=axcolor)
        self.recettesSlider = Slider(recettesAxis, r'Recettes (€)', np.min(self.recettesArray), np.max(self.recettesArray), self.recettes)
        self.recettesSlider.on_changed(self.updateParam)
        
        vinfAxis = plt.axes([0.7, 0.55, 0.2, 0.03], facecolor=axcolor)
        self.vinfSlider = Slider(vinfAxis, r'$v_{0}$', np.min(self.vinfArray), np.max(self.vinfArray), self.vinf)
        self.vinfSlider.on_changed(self.updateParam) 
        
        aAxis = plt.axes([0.7, 0.5, 0.2, 0.03], facecolor=axcolor)
        self.aSlider = Slider(aAxis, r'a', np.min(self.aArray), np.max(self.aArray), self.a)
        self.aSlider.on_changed(self.updateParam)
        
        bAxis = plt.axes([0.7, 0.45, 0.2, 0.03], facecolor=axcolor)
        self.bSlider = Slider(bAxis, r'b', np.min(self.bArray), np.max(self.bArray), self.b)
        self.bSlider.on_changed(self.updateParam)
        
        cAxis = plt.axes([0.7, 0.4, 0.2, 0.03], facecolor=axcolor)
        self.cSlider = Slider(cAxis, r'c', np.min(self.cArray), np.max(self.cArray), self.c)
        self.cSlider.on_changed(self.updateParam)
        
        dAxis = plt.axes([0.7, 0.35, 0.2, 0.03], facecolor=axcolor)
        self.dSlider = Slider(dAxis, r'd', np.min(self.dArray), np.max(self.dArray), self.dArray[25])
        self.dSlider.on_changed(self.updateParam)
        
        eAxis = plt.axes([0.7, 0.3, 0.2, 0.03], facecolor=axcolor)
        self.eSlider = Slider(eAxis, r'e', np.min(self.eArray), np.max(self.eArray), self.eArray[25])
        self.eSlider.on_changed(self.updateParam)
        
    def updateComputation(self): # The computation, parameter dependant
        #first deal with abonnements
        mask= self.status=='RS'
        self.abo= self.abop*np.ones((self.nbHab,))
        self.abo[self.status=='RS']=self.abos
        
        #second deal with consumption and total prices
        self.prixConso=np.zeros(self.vv.shape)
        self.prixp=np.zeros(self.vv.shape)
        self.prixs=np.zeros(self.vv.shape)
        if np.sum(self.abo)>=self.recettes:
            self.p0=0
            self.prixp = self.abop
            self.prixs = self.abos
        else: 
            #find p0 to balance income. 
            alpha1=np.zeros((self.nbHab,))# préfacteur de p_0 dans la tranche 1
            alpha2=np.zeros((self.nbHab,))# préfacteur de p_0 dans la tranche 2
            beta2=np.zeros((self.nbHab,))# constante dans la tranche 2

            #compute consumption integral on all inhabitants
            for ih in range(self.nbHab):
                alpha1[ih],alpha2[ih],beta2[ih]=self.computeIntegrals(self.volume[ih])
        
            self.p0=(self.recettes-np.sum(beta2+self.abo))/np.sum(alpha1+alpha2) # prix d'inflexion
            self.prixProj = (alpha1+alpha2)*self.p0+beta2 #Calcul du prix payé par les habitantes dans le nouveau scénario
            
            #compute consumption price per m3
            #tranche 1
            self.vv[0:self.nbpts]=self.vinf*((1-3*self.b)*(self.tt**3)+3*self.b*(self.tt**2))
            self.prixConso[0:self.nbpts]=self.p0*((3*self.a-2)*(self.tt**3)+(-6*self.a+3)*(self.tt**2)+3*self.a*self.tt) 
            #tranche 2
            self.vv[self.nbpts:]=self.vinf+(self.vmax-self.vinf)*((3*(self.c+self.d-self.c*self.d)-2)*(self.tt**3)+\
                                                             3*(1-2*self.c-self.d+self.c*self.d)*(self.tt**2)+3*self.c*self.tt)
            self.prixConso[self.nbpts:]=self.p0+(self.pmax-self.p0)*((1-3*self.e)*(self.tt**3)+3*self.e*(self.tt**2)) 
            
            #compute full consumption price (integral of price per m3)
            alpha1=np.zeros(self.vv.shape)# préfacteur de p_0 dans la tranche 1
            alpha2=np.zeros(self.vv.shape)# préfacteur de p_0 dans la tranche 2
            beta2=np.zeros(self.vv.shape)# constante dans la tranche 2
            for iv,vv in enumerate(self.vv): 
                alpha1[iv],alpha2[iv],beta2[iv]=self.computeIntegrals(vv)
            self.prixp = self.abop+(alpha1+alpha2)*self.p0+beta2
            self.prixs = self.abos+(alpha1+alpha2)*self.p0+beta2
            
        #now compute prices
        self.prixp_m3 = self.prixp/self.vv
        self.prixs_m3 = self.prixs/self.vv
            
    def computeIntegrals(self,vv):#return alpha1, alpha2 and beta for an input volume
        a=self.a
        b=self.b
        c=self.c
        d=self.d
        e=self.e
        if vv<=self.vinf:# on ne travaille que dans la tranche 1.
             #Find the value of T by solving the equation V=v(T)
            if vv==0: #remove extrema points where root finding can be problematic
                T=0.0
            elif vv==self.vinf:
                T=1.0
            else: 
                p = [1-3*b, 3*b, 0, -vv/self.vinf]# polynomial representation of the equation
                roots = np.roots(p)
                roots = np.unique(roots)#remove duplicates
                roots2 = np.real(roots[np.isreal(roots)]) # take only real roots 
                mask = (roots2<=1.) & (roots2>=0.)#check that T is real and between 0 and 1
                T=float(roots2[mask])
            alpha1=3*self.vinf*(T**6/6*(-9*a*b+3*a+6*b-2)+T**5/5*(24*a*b-6*a-13*b+3)+3*T**4/4*(-7*a*b+a+2*b)+T**3/3*6*a*b) 
            return alpha1,0,0
        else:
            alpha1=3*self.vinf*(1/6*(-9*a*b+3*a+6*b-2)+1/5*(24*a*b-6*a-13*b+3)+3/4*(-7*a*b+a+2*b)+1/3*6*a*b)# tranche 1 avec T=1
            #tranche 2
            wmax=self.vmax-self.vinf
            #Find the value of T by solving the equation V-vinf=w(T)
            if vv==self.vinf: #remove extrema points where root finding can be problematic
                T=0.0
            elif vv==self.vmax:
                T=1.0
            else:
                p = [3*(c+d-c*d)-2, 3*(1-2*c-d+c*d), 3*c, -(vv-self.vinf)/wmax]
                roots = np.roots(p)
                roots = np.unique(roots)#remove duplicates
                roots2 = np.real(roots[np.isreal(roots)]) # take only real roots 
                mask = (roots2<=1.) & (roots2>=0.)#check that T is real and between 0 and 1#
                if not mask.any(): # for a weird reason, if vv=vmax, we have roots2=1.0 but it does not pass roots2<=1.0 ???? 
                    print(vv,roots,roots2,mask)
                T=float(np.real(roots2[mask]))
            uu=(-3*c*d+9*e*c*d+3*c-9*e*c+3*d-9*e*d+6*e -2)*T**6/6+(2*c*d-15*e*c*d-4*c+21*e*c-2*d+15*e*d-12*e+2)*T**5/5 +(6*e*c*d+c-15*e*c-6*e*d+6*e)*T**4/4+(3*e*c)*T**3/3
            alpha2=vv-self.vinf-3*uu*wmax
            beta2=3*self.pmax*wmax*uu
            return alpha1,alpha2,beta2

    def updatePlot(self): # The plot update
        self.ax1.plot(self.vv,self.prixp,label='RP/PRO') 
        self.ax1.plot(self.vv,self.prixs,label='RS') 
        self.ax1.scatter(self.volume,self.prix2018,s=3,color='g',label='2018')
        self.ax2.plot(self.vv,self.prixp_m3,label='RP/PRO') 
        self.ax2.plot(self.vv,self.prixs_m3,label='RS') 
        self.ax2.plot(self.vv,self.prixConso,'k-',label='Conso') 
        self.ax2.plot([self.vinf,self.vinf],[0,self.p0],color = '0.75')
        self.ax2.plot([0,self.vinf],[self.p0,self.p0],color = '0.75')
        self.ax2.scatter(self.volume,self.prix2018_m3,s=3,color='g',label='2018')
        
        self.ax1.set_xlim(self.axx)
        self.ax1.set_ylim(self.ax1y)
        self.ax2.set_ylim(self.ax2y)
        #self.ax2.text(0.9, 0.2, '$p_0=${:.2f} €'.format(self.p0),ha='right', va='bottom', transform=self.ax2.transAxes)
        self.ax2.legend()
        self.ax2.set_xlabel(r'volume ($m^3$) ')
        self.ax1.set_ylabel('Facture totale (€)')
        self.ax2.set_ylabel('Prix au m3 (€)')
        plt.tight_layout()
        plt.subplots_adjust(top = 0.9,right = 0.5) 

    def updateParam(self,val): # Update the parameter using the slider values
        #get current values of ax limits
        self.axx=self.ax1.get_xlim()
        self.ax1y=self.ax1.get_ylim()
        self.ax2y=self.ax2.get_ylim()
        self.ax1.clear()
        self.ax2.clear()
        self.abop= self.abopSlider.val
        self.abos= self.abosSlider.val
        self.recettes= self.recettesSlider.val
        self.vinf=self.vinfSlider.val
        self.a=self.aSlider.val
        self.b=self.bSlider.val
        self.c=self.cSlider.val
        self.d=self.dSlider.val
        self.e=self.eSlider.val
        
        self.updateComputation()
        self.updatePlot()  # 

        
class CurrentModel():
    def __init__(self): # Set the static parameters
        
        self.vv = np.linspace(0,2100,2101) # volume values
        self.vv[0]=1e-5 #little trick to avoid divide by 0
        self.abopArray = np.linspace(0,200,201) # potential values for abop
        self.abosArray = np.linspace(0,200,201) # potential values for abos
        self.recettesArray = np.linspace (50000,100000,51) # potential values for recettes
        
        #loading the inhabitant data
        book = xlrd.open_workbook('Eau2018.xls')
        sheet = book.sheet_by_name('CALCULS')
        self.nbHab=363
        self.volume = np.array([sheet.cell_value(r,4) for r in range(1,self.nbHab+1)])#get the volume used in 2018 in m3
        self.vTot=np.sum(self.volume)
        self.volume[self.volume==0]=1e-5 #little trick to avoid divide by 0
        self.status = np.array([sheet.cell_value(r,3) for r in range(1,self.nbHab+1)])#get the RS, RP, PRO status
        self.prix2018 = np.array([sheet.cell_value(r,33) for r in range(1,self.nbHab+1)])#get the price paid in 2018 in EUR
        self.prix2018_m3 = self.prix2018 / self.volume
        

    def mainFunction(self):
        '''
        main function, the one called in the notebook, and that update computations depending on the slider values.
        '''
        self.fig = plt.figure(figsize=(8,5))
        self.ax1  = plt.subplot(121)
        self.ax2  = plt.subplot(122,sharex=self.ax1)
        self.ax1.set_xlabel(r'volume ($m^3$) ')
        self.ax2.set_xlabel(r'volume ($m^3$) ')
        self.ax1.set_ylabel('Facture totale (€)')
        self.ax2.set_ylabel('Prix au m3 (€)')
        
        # initial values
        self.abop= self.abopArray[100]
        self.abos= self.abosArray[100]
        self.recettes= self.recettesArray[25]
        
        # initial plot
        self.updateComputation()
        self.updatePlot()
    
        # SLIDERS
        axcolor   = 'lightgoldenrodyellow'
        abopAxis = plt.axes([0.2, 0.15, 0.65, 0.05], facecolor=axcolor)
        self.abopSlider = Slider(abopAxis, r'Abo. P/PRO (EUR)', np.min(self.abopArray), np.max(self.abopArray), self.abopArray[100])
        self.abopSlider.on_changed(self.updateParam)
        abosAxis = plt.axes([0.2, 0.1, 0.65, 0.05], facecolor=axcolor)
        self.abosSlider = Slider(abosAxis, r'Abo. S (EUR)', np.min(self.abosArray), np.max(self.abosArray), self.abosArray[100])
        self.abosSlider.on_changed(self.updateParam)        
        recettesAxis = plt.axes([0.2, 0.05, 0.65, 0.05], facecolor=axcolor)
        self.recettesSlider = Slider(recettesAxis, r'Recettes (EUR)', np.min(self.recettesArray), np.max(self.recettesArray), self.recettesArray[25])
        self.recettesSlider.on_changed(self.updateParam)
        
    def updateComputation(self): # The computation, parameter dependant
        #first deal with abonnements
        mask= self.status=='RS'
        self.abo= self.abop*np.ones((self.nbHab,))
        self.abo[self.status=='RS']=self.abos
        
        #second deal with consumption
        if np.sum(self.abo)>=self.recettes:
            self.p0=0
        else: 
            self.p0= (self.recettes-np.sum(self.abo))/self.vTot
        
        #now deal with prices
        self.prixp = self.abop+self.p0*self.vv
        self.prixs = self.abos+self.p0*self.vv
        self.prixp_m3 = self.abop/self.vv+self.p0
        self.prixs_m3 = self.abos/self.vv+self.p0

    def updatePlot(self): # The plot update
        self.ax1.plot(self.vv,self.prixp,label='RP/PRO') 
        self.ax1.plot(self.vv,self.prixs,label='RS') 
        self.ax1.scatter(self.volume,self.prix2018,s=3,color='g',label='2018')

        self.ax2.plot(self.vv,self.prixp_m3,label='RP/PRO') 
        self.ax2.plot(self.vv,self.prixs_m3,label='RS') 
        self.ax2.scatter(self.volume,self.prix2018_m3,s=3,color='g',label='2018')
        self.ax2.set_ylim(0,20)
        self.ax2.text(0.9, 0.2, '$p_0=${:.2f} €'.format(self.p0),ha='right', va='bottom', transform=self.ax2.transAxes)
        self.ax2.legend()
        
        self.ax1.set_xlabel(r'volume ($m^3$) ')
        self.ax2.set_xlabel(r'volume ($m^3$) ')
        self.ax1.set_ylabel('Facture totale (€)')
        self.ax2.set_ylabel('Prix au m3 (€)')
        plt.tight_layout()
        plt.subplots_adjust(top = 0.9,bottom = 0.35) 
        
    def updateParam(self,val): # Update the parameter using the slider values
        self.ax1.clear()
        self.ax2.clear()
        self.abop= self.abopSlider.val
        self.abos= self.abosSlider.val
        self.recettes= self.recettesSlider.val
        self.updateComputation()
        self.updatePlot()  #