precipitation_fitter.py

# -*- coding: utf-8 -*-
"""
Created on Tue Oct 26 15:33:11 2021

@author: pgiol
"""
import numpy as np

import itertools
import scipy.optimize as optimize
from sklearn.metrics import r2_score

import matplotlib.pyplot as plt
import matplotlib.cm as cm

from decimal import *

import scipy.sparse
import scipy.sparse.linalg
from scipy import sparse
from scipy.integrate import odeint



getcontext().prec = 256



TitraM = np.array([
    [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0] #%WaterContent
    ,[0.017202142,0.052222222,0.59832664,0.691097724,0.921659973,0.974431058,1.147215529,1.044431058,0.901659973,0.958875502]
    ,[0.005612,0.016823,0.08972,0.120561,0.072897,0.075708,0.145802,0.131775,0.120561,0.092523]
    ,[0.003262,0.215883,0.871488,1.028035,1.013547,1.032717,1.037855,1.014956,0.997656,1.042054] #3A
    ,[0.011212,0.056068,0.039252,0.056075,0.028038,0.028037,0.016822,0.025234,0.014018,0.019626]
    ,[0.014469,0.580369,0.882703,1.061687,1.033173,0.976636,0.987381,0.998126,0.989252,0.941112] #4A
    ,[0.011215,0.098131,0.047663,0.044867,0.019626,0.042057,0,0.036441,0.019626,0.019626]
    ,[0.008876,0.642051,1.011675,0.927108,1.094848,1.021481,0.964958,0.972892,1.025701,1.005605] #4.5A
    ,[0.005607,0.064486,0.095327,0.014019,0.058871,0.042049,0.019626,0.011215,0.030841,0.056075]
    ,[0.01168,0.54392,0.9556,0.991586,1.047192,1.052336,1.00421,0.992519,1.042523,1.033642]       #7E
    ,[0.011215,0.064485,0.061682,0.022429,0.011215,0.044866,0.030841,0.100935,0.070093,0.036448]])

    # [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0] #%WaterContent
    # ,[0.017202142,0.052222222,0.59832664,0.691097724,0.921659973,0.974431058,1.147215529,1.044431058,0.901659973,0.958875502] #2A
    # ,[0.005612,0.016823,0.08972,0.120561,0.072897,0.075708,0.145802,0.131775,0.120561,0.092523] 
    # ,[1.78E-15,0.216216216,0.872972973,1.02972973,1.016216216,1.040540541,1.040540541,1.013513514,1.013513514,1.045945946]
    # ,[0.021621622,0.056756757,0.067567568,0.054054054,0.089189189,0.021621622,0.018918919,0.027027027,0.035135135,0.02972973]
    # ,[2.70E-03,0.575675676,0.878378378,1.062162162,1.02972973,0.986486486,0.986486486,1.002702703,0.997297297,0.945945946]
    # ,[0,0,0,0,0,0,0,0,0,0]
    # ,[0.008876,0.642051,1.011675,0.927108,1.094848,1.021481,0.964958,0.972892,1.025701,1.005605] #4.5A
    # ,[0.005607,0.064486,0.095327,0.014019,0.058871,0.042049,0.019626,0.011215,0.030841,0.056075]
    # ,[0.01168,0.54392,0.9556,0.991586,1.047192,1.052336,1.00421,0.992519,1.042523,1.033642]       #7E
    # ,[0.011215,0.064485,0.061682,0.022429,0.011215,0.044866,0.030841,0.100935,0.070093,0.036448]])



def Fit_model(inputs, a, b, c, d):
    a = 1
    b = 0
    x = (inputs[:])
    return  a + (b-a)/(a + np.power(x/c,d))

Ti = 5

newp, pcov_new = optimize.curve_fit( Fit_model, (TitraM[0,:]*0.6), TitraM[Ti,:], ftol=1e-15, xtol=1e-15, maxfev=800000)

xi = np.linspace(min(TitraM[0,:]*0.6)*0.25,max(TitraM[0,:]*0.6)*1.25,100)#len(TitraM[0,:]))
Fnew = Fit_model((TitraM[0,:]*0.6), *newp)
rsq_new = r2_score(TitraM[1,:], Fnew)
Fit = Fit_model((xi),*newp)

print(newp)

def PlotFit(): #Plots power law fit over data
    fig = plt.figure(15)
    plt.plot(TitraM[0,:]*0.6,TitraM[Ti,:],'or',label='Experimental Data')
    plt.plot(xi,Fit,'red',label='P=%5.2f $\delta^{%5.2f}, R^2$ = %1.5f' %(newp[0],newp[1],rsq_new))
    plt.legend(loc='best')
    plt.xlabel('Disp.')
    plt.ylabel('Load')
    plt.show()


PlotFit()