Week 5 - Curve fitting

# importing modules
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit


# Linear Curve fit function
def function_linear_curve(t, a, b):
    return a * t + b


# Cubic Curve fit function
def function_cubic_curve(t, a, b, c, d):
    return a * pow(t, 3) + b * pow(t, 2) + c * t + d


# Reading thermodynamic data file
def read_file():
    temperature = []
    cp = []
    for line in open('data', 'r'):
        values = line.split(',')
        temperature.append(float(values[0]))
        cp.append(float(values[1]))

    return [temperature, cp]

# read file data by calling function
temperature, cp = read_file()

print("temperature",temperature)
print("cp",cp)
# calculation of coefficient for linear curve
popt, pcov = curve_fit(function_linear_curve, temperature, cp)
fit_cp_l = function_linear_curve(np.array(temperature), *popt)
# print("fit_cp_c",fit_cp_l)

# plotting data for linear curve
plt.figure(1)
plt.plot(temperature, cp, 'b--')
plt.plot(temperature, fit_cp_l, color='red')
plt.legend(['Actual data', 'Curve fit'])
plt.xlabel('Temperature (K)')
plt.ylabel('Cp')
plt.title('Linear Curve fitting')
plt.show()

# calculation of coefficient for cubic curve
popt, pcov = curve_fit(function_cubic_curve, temperature, cp)
fit_cp_c = function_cubic_curve(np.array(temperature), *popt)

# plotting data for cubic curve
plt.figure(2)
plt.plot(temperature, cp, 'b--')
plt.plot(temperature, fit_cp_c, color='red')
plt.legend(['Actual data', 'Curve fit'])
plt.xlabel('Temperature (K)')
plt.ylabel('Cp')
plt.title('Cubic Curve fitting')
plt.show()


#Measuring the fitness characetristics of the curves
s = np.sum(cp)
l = np.size(cp)

# Finding mean of Cp data
m = s / l 
print('Mean of all cp :', m,'\n')

# For Linear curve fit
linear_sse = 0
linear_ssr = 0
for i in range(l):
    linear_error = abs((np.sum((cp[i] - fit_cp_l[i]))))
   
    linear_sse = linear_sse + pow(linear_error, 2)
    linear_ssr = linear_ssr + np.sum(pow((fit_cp_l[i] - m), 2))

print('For Linear Curve Fit')
linear_sst = linear_sse + linear_ssr
print('Linear SST :', linear_sst)
linear_R2 = linear_ssr / linear_sst
print('Linear R2 :', linear_R2)
linear_rmse = pow((linear_sse / l), 0.5)
print('Linear RMSE :', linear_rmse ,'\n')

# For Cubic curve fit
cubic_sse = 0
cubic_ssr = 0
for j in range(l):
    cubic_error = abs(np.sum((cp[i] - fit_cp_c[i])))
    cubic_sse = cubic_sse + (pow(cubic_error, 2))
    cubic_ssr = cubic_ssr + np.sum(pow((fit_cp_c[i] - m), 2))

print('For Cubic Curve Fit')
cubic_sst = cubic_sse + cubic_ssr
print('Cubic SST :', cubic_sst)
cubic_R2 = cubic_ssr / cubic_sst
print('Cubic R2 :', cubic_R2)
cubic_rmse = pow((cubic_sse / l), 0.5)
print('Cubic RMSE :', cubic_rmse)