Probability long term strength model

import math
import operator
import numpy as np
from sklearn import linear_model
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from numpy import arange,array,ones,polyfit,polyval#,random,linalg
from pylab import plot,show,figure,legend
from scipy import stats
import math
#
#
# stress = [70, 70, 70, 70, 70, 45, 45, 45, 300, 300, 300, 300,\
#           300, 300, 300, 300, 300, 300, 200, 200, 200, 200, 200]
# temperature = [523, 523, 523, 523, 523, 573, 573, 573, 423,\
#                423, 423, 423, 423, 423, 423, 423, 423, 423, 473, 473, 473, 473, 473]
# time = [254, 253, 255, 256, 252, 160, 162, 158, 165, 155, 157,\
#         163, 164, 156, 153, 167, 155, 165, 17.3, 17, 16.7, 17.5, 16.5]
import pandas as pd

parameters = ['Stress', 'Temperature', 'MinimumCreepRate', 'Priorplasticdeformation', 'Durationofprimarycreep',\
              'Timetothirdcreepperiod', 'Timetorupture','Deformation', 'Time', 'MinimumCreepRate2', 'Time3',    'MinimumCreepRate3']

path = "C:\\PostGraduate\\Experiments_new.csv"

def read_csv():
    data = pd.read_csv(path, sep=',',
                          names=parameters, skiprows=1).astype('float')
    return data

def process_calc():
    data = read_csv()
    result = {}
    n, h, a0 = get_first_params(data, result)
    l = get_l(data, result)
    m, p, b0 = get_second_params(data, result)
    data['a'] = (data['MinimumCreepRate2'] / (data['Stress'].pow(n) * np.exp(-h/data['Temperature'])))
    data['b'] = (1/((l+1)*data['Timetorupture']*np.exp(-p/data['Temperature'])*data['Stress'].pow(m)))
    data['def_first'] = data['MinimumCreepRate2'] * data['Durationofprimarycreep']
    # data['teta'] = data['a'] * np.exp(-h/data['Temperature']) * data['Stress'].pow(n) * data['Durationofprimarycreep']

    temp = 1-(l+1) * data['b'] * np.exp(-p/data['Temperature']) * data['Stress'].pow(m) * data['Durationofprimarycreep']
    data['temp'] = temp
    data['teta'] = data['a'] * np.exp((p-h)/data['Temperature']) * data['Stress'].pow(n-m) * \
                   ( data['temp'].pow((l - n +1) / (l + 1)) - 1) / (data['b'] * (n - l - 1))

    return result, data

def get_first_params(data, result):
    x1 = np.log(data['Stress'])
    x2 = - data['Temperature'].pow(-1)
    y = np.log(data['MinimumCreepRate2'])
    # linregress doesn't do multi regression, so we use sklearn
    x = np.vstack([x1, x2]).T # don't need ones for
    regr = linear_model.LinearRegression()
    regr.fit( x, y )
    coef = regr.coef_
    intercept = regr.intercept_
    result['a0'] = np.exp(intercept)
    result['n'] = coef[0]
    result['h'] = coef[1]
    return coef[0], coef[1], np.exp(intercept)

def get_l(data, result):
    n = result['n']
    a0 = result['a0']
    h = result['h']
    data['Minus_Temperature_reverse'] = - data['Temperature'].pow(-1)
    data['l'] = n - 1 + n/((data['Deformation']/\
                            (a0*np.exp(h*data['Minus_Temperature_reverse'])*data['Stress'].pow(n)*data['Timetorupture'])) - 1)
    l = data['l'].mean()
    result['l'] = l
    return l

def get_second_params(data, result):
    x1 = - np.log(data['Stress'])
    x2 = data['Temperature'].pow(-1)
    y = np.log(data['Timetorupture'])
    l = result['l']
    # linregress doesn't do multi regression, so we use sklearn
    x = np.vstack([x1, x2]).T # don't need ones for
    regr = linear_model.LinearRegression()
    regr.fit( x, y )
    coef = regr.coef_
    intercept = regr.intercept_
    result['b0'] = np.exp(-intercept)/(l+1)
    result['m'] = coef[0]
    result['p'] = coef[1]
    return coef[0], coef[1], np.exp(-intercept)/(l+1)


def get_characteristics(data):

    means = []
    stds = []
    stresses = []
    temperatures = []
    for k, group in data.groupby(['Stress', 'Temperature']):
        stress = k[0]
        temperature = k[1]
        data = group['b']
        # mean = data.mean()
        mean, std = norm.fit(data)
        # var = data.var()
        print('stress = ', stress, 'temperature = ', temperature, \
              'mean = ', mean, 'std = ', std)
        stresses.append(stress)
        temperatures.append(temperature)
        means.append(mean)
        stds.append(std)
    return stresses, temperatures, means, stds

def get_params_multiple_regressions(depend_var1, depend_var2, independ_var):
    x = np.vstack([depend_var1, depend_var2]).T # don't need ones for
    regr = linear_model.LinearRegression()
    regr.fit(x, independ_var)
    return regr.coef_, regr.intercept_

def get_params_multiple_regressions_many(depend_var1, depend_var2, independ_var):
    # depend_var1 - stress
    # depend_var2 - temperature
    depend_var1_sq = [dep*dep for dep in depend_var1]
    depend_var2_sq = [dep*dep for dep in depend_var2]
    depend_var1_var2_sq = [depend_var1[index]*depend_var2[index] for index in range(len(depend_var1))]
    x = np.vstack([depend_var1_sq, depend_var2_sq, depend_var1_var2_sq, depend_var1, depend_var2]).T # don't need ones for
    regr = linear_model.LinearRegression()
    regr.fit(x, independ_var)
    return regr.coef_, regr.intercept_

def get_value_by_params(coefficient, intercept, params):
    summa = 0
    for i in range(len(coefficient)):
        summa+=params[i]*coefficient[i]
    summa+=intercept
    return summa

def calculate_time(mean_b, std_b, result, stress, temperature):
    from numpy import random
    r = random.normal(mean_b, std_b, size = 1000)
    return 1/((result['l']+1)*r*np.exp(-result['p']/temperature)*math.pow(stress, result['m']))

def plot_graphs(time):
    plot_pdf(time)
    plot_cdf(time)

def plot_pdf(time):
    mu, std = norm.fit(time)

    # Plot the histogram.
    plt.hist(time, bins=25, normed=True, alpha=0.6, color='g')

    # Plot the PDF.
    xmin, xmax = plt.xlim()
    x = np.linspace(xmin, xmax, 100)
    p = norm.pdf(x, mu, std)
    plt.plot(x, p, 'k', linewidth=2)
    title = "Fit results: mu = %.2f,  std = %.2f" % (mu, std)
    plt.title(title)

    plt.show()

def plot_cdf(time):
    mu, std = norm.fit(time)

    # Plot the histogram.
    plt.hist(time, bins=25, normed=True, alpha=0.6, color='g')

    # Plot the CDF.
    xmin, xmax = plt.xlim()
    x = np.linspace(xmin, xmax, 100)
    c = norm.cdf(x, mu, std)
    plt.plot(x, c, 'k', linewidth=2)
    title = "Fit results: mu = %.2f,  std = %.2f" % (mu, std)
    plt.title(title)

    plt.show()

def get_selected_stresses(result, data):
    temperature_2_selected_stress = {}
    for temperature, data_temp in data.groupby(['Temperature']):
        # data_filtered = data[data['Temperature']==temperature]
        median_deform = data_temp['def_first'].median()
        selected_stress = get_selected_stress(data_temp, median_deform)
        temperature_2_selected_stress[temperature] = selected_stress
    return list(set(temperature_2_selected_stress.values()))

def get_params_c_k(result, data):
    c_list = []
    k_list = []
    selected_stresses = get_selected_stresses(result, data)
    data_filtered = data[data['Stress'].isin(selected_stresses)]
    for stress, group in data_filtered.groupby(['Stress']):
        def_first = group['def_first'].values
        teta = group['teta'].values
        c, k = get_estimations_params(def_first, teta)
        c_list.append(c)
        k_list.append(k)
    filtered_c_list = [c for c in c_list if c != 0]
    filtered_k_list = [k for k in k_list if k != 0]
    return np.mean(filtered_c_list), np.mean(filtered_k_list)

def get_charact_c_k_2(result, data):
    stresses = []
    temperatures = []
    c_list = []
    k_list = []
    for k, group in data.groupby(['Stress', 'Temperature']):
        stress = k[0]
        temperature = k[1]
        def_first = group['def_first'].values
        teta = group['teta'].values
        c, k = get_estimations_params(def_first, teta)
        if c != 0 or k != 0:
            stresses.append(stress)
            temperatures.append(temperature)
            c_list.append(c)
            k_list.append(k)
    return stresses, temperatures, c_list, k_list

def get_params_c_k_2(stress, temperature, result, data):
    stresses, temperatures, c_list, k_list = get_charact_c_k_2(result, data)
    params = [stress * stress, temperature * temperature, stress * temperature, stress, temperature]

    # estimate c
    mean_coef_c, mean_intercept_c = get_params_multiple_regressions_many(stresses, temperatures, c_list)
    mean_c = get_value_by_params(mean_coef_c, mean_intercept_c, params)

    # estimate k
    mean_coef_k, mean_intercept_k = get_params_multiple_regressions_many(stresses, temperatures, k_list)
    mean_k = get_value_by_params(mean_coef_k, mean_intercept_k, params)
    return mean_c, mean_k

def estimation_function(x, c, k):
    return k*np.log((1+c)*np.exp(x/k) - c)

def func(x, p1,p2):
  return p1*np.cos(p2*x) + p2*np.sin(p1*x)

def get_estimations_params(def_first, teta):
    try:
        popt, pcov = curve_fit(estimation_function, teta, def_first)
        c = popt[0]
        k = popt[1]
        data_checked = estimation_function(teta, c, k)
        residuals = def_first - data_checked
        error = np.sqrt(sum(residuals**2))
        return c, k
    except:
        return 0, 0

def get_selected_stress(data_filtered, median_deform):
    stress_2_distance = {}
    for stress, data_stress in data_filtered.groupby(['Stress']):
        # deform_first = data_filtered[data_filtered['Stress'] == stress]
        median = data_stress['def_first'].median()
        stress_2_distance[stress] = math.fabs(median_deform - median)
    stress_res = sorted(stress_2_distance.items(), key=operator.itemgetter(1))
    return stress_res[0][0]

def calculate_time_to_failure(data, result, stress, temperature):
    stresses, temperatures, means, stds = get_characteristics(data)
    # mean_coef, mean_intercept = get_params_multiple_regressions(stresses, temperatures, means)
    # std_coef, std_intercept = get_params_multiple_regressions(stresses, temperatures, stds)
    # params = [stress, temperature]
    mean_coef, mean_intercept = get_params_multiple_regressions_many(stresses, temperatures, means)
    std_coef, std_intercept = get_params_multiple_regressions_many(stresses, temperatures, stds)
    params = [stress * stress, temperature * temperature, stress * temperature, stress, temperature]
    mean_b = get_value_by_params(mean_coef, mean_intercept, params)
    std_b = get_value_by_params(std_coef, std_intercept, params)
    time = calculate_time(mean_b, std_b, result, stress, temperature)
    return time

def mean_calculation(stress, temperature):
    result, data = process_calc()
    time = calculate_time_to_failure(data, result, stress, temperature)
    # plot_graphs(time)
    # print('mean_Time = ', np.mean(time), 'std_Time = ', np.std(time))
    # c, k = get_params_c_k(result, data)
    c, k = get_params_c_k_2(stress, temperature, result, data)
    result['c'] = c
    result['k'] = k
    print('res')


if __name__ == '__main__':
    # stress = 45
    # temperature = 573
    stresses_test = [70, 45, 300, 200]
    temperatures_test = [523, 573, 423, 473]
    for i in range(len(stresses_test)):
        mean_calculation(stresses_test[i], temperatures_test[i])

# mean_Time =  202.977268144 std_Time =  35.3658215247
# mean_Time =  229.291468263 std_Time =  23.530169309
# mean_Time =  36.3147799399 std_Time =  21.4642344076
# mean_Time =  362.608885147 std_Time =  48.3941602839