OiO.lk Community platform!

Oio.lk is an excellent forum for developers, providing a wide range of resources, discussions, and support for those in the developer community. Join oio.lk today to connect with like-minded professionals, share insights, and stay updated on the latest trends and technologies in the development field.
  You need to log in or register to access the solved answers to this problem.
  • You have reached the maximum number of guest views allowed
  • Please register below to remove this limitation

Calculating the locations of the maximum and mean, and the standard deviation between them, with respect to local and absolute maxima, within dataset

  • Thread starter Thread starter PuppyLord
  • Start date Start date
P

PuppyLord

Guest
I have multiple datasets of data, for which a compound sigmoid fitting curve is plotted to each. The derivatives of each fit curve are taken, resulting in 6 different theoretical derivative curves. Each derivative curve has two maxima, an absolute maximum at lower temperatures (x-axis) and a local maxima at higher ranges.

Per each of the 6 datasets, I need:

  1. the location of the local and absolute maximum (x-axis position) calculated and printed
  2. the location of their means, with respect to each peak/distribution calculated
  3. the location of the medians with respect to each peak/distribution calculated
  4. the average standard deviation for each peak/distribution
  5. Finally, the distance between the means and their respective maxima, converted into a standard deviation value.

We would therefore have printed in terminal:

Dataset 1: Location of absolute maximum: Mean of absolute maximum: Median of absolute maximum: Error of mean from absolute maximum (not distance; the standard deviation equiv. max. - mean location): Average standard deviation of that maximum:

Location of local maximum: Mean of local maximum: Median of local maximum: Error of mean from local maximum (not distance; the standard deviation equiv. max. - mean location): Average standard deviation of that maximum:

I have tried extremely hard and ended up creating a Frankenstein code from other types of specimens I analyzed. So, the terminal plots multiple values.

Code:
Often, the standard deviations are too big, as can be visually seen in the PLOTS that are generated. The code can't just print the absolute maximum and local maximum-specific statistics. And can't accurately convert the distance between the means and their respective maxima into standard deviation values (Standard Error). The plots are perfect, please use this as reference. 

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.signal import find_peaks
from sklearn.metrics import r2_score

#1. User datasets: 
x_data = np.array([25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115,       
120, 125, 130])
y_data50 = np.array([0, 0.50031, 4.00248, 9.50589, 13.50837, 15.0093, 16.00992, 17.01054, 18.01116,   
19.01178, 20.51271, 21.51333, 23.01426, 25.0155, 27.51705, 30.51891, 32.52015, 33.270615, 33.770925, 
33.770925, 34.02108, 34.271235]) 
y_data100 = np.array([0, 4.95405, 12.88053, 20.80701, 28.73349, 36.65997, 39.6324, 42.109425,  
44.58645, 47.55888, 50.035905, 52.51293, 55.48536, 59.4486, 63.41184, 68.36589, 72.32913, 73.31994, 
74.31075, 75.30156, 75.30156, 75.30156])
y_data150 = np.array([0, 2.96262, 8.88786, 35.55144, 51.84585, 57.77109, 62.21502, 65.17764, 
68.14026, 72.58419, 75.54681, 79.250085, 83.694015, 88.8786, 94.063185, 99.24777, 107.394975, 
115.54218, 119.98611, 122.94873, 124.43004, 124.43004])
y_data200 = np.array([0, 5.91543, 29.57715, 70.98516, 88.73145, 92.67507, 98.5905, 102.53412, 
108.44955, 112.39317, 118.3086, 126.19584, 132.11127, 139.99851, 151.82937, 161.68842, 173.51928,   
175.49109, 177.4629, 179.43471, 179.43471, 179.43471])
y_data250 = np.array([0, 12.31155, 39.39696, 83.71854, 105.87933, 115.72857, 120.65319, 128.04012, 
135.42705, 142.81398, 152.66322, 162.51246, 169.89939, 182.21094, 194.52249, 206.83404, 219.14559, 
221.6079, 224.07021, 226.53252, 226.53252, 226.53252])
y_data300 = np.array([0, 11.81124, 38.38653, 76.77306, 121.06521, 132.87645, 138.78207, 147.6405, 
153.54612, 159.45174, 168.31017, 180.12141, 188.97984, 203.74389, 218.50794, 236.2248, 248.03604,   
253.94166, 259.84728, 261.323685, 262.80009, 265.7529])

#2a. Sigmoid function definition: 
def compound_sigmoid(x, U1, x01, k1, U2, x02, k2):
    return U1 / (1 + np.exp(-k1 * (x - x01))) + U2 / (1 + np.exp(-k2 * (x - x02)))

#2b. Derivative of sigmoid function:
def compound_sigmoid_derivative(x, U1, x01, k1, U2, x02, k2):
    sig1 = (U1 * k1 * np.exp(-k1 * (x - x01))) / ((1 + np.exp(-k1 * (x - x01))) ** 2)
    sig2 = (U2 * k2 * np.exp(-k2 * (x - x02))) / ((1 + np.exp(-k2 * (x - x02))) ** 2)
    return sig1 + sig2

#3. Initial guesses for calculating fitting parameters: USER SELECTION REQUIRED
#(1-A.) [PERSONAL/EACH] >WORK DENSITY< Initial Guesses for [STRETCHED]: 
p0_50 = [max(y_data50)/2, 38, -0.1, max(y_data50), 65, -0.1] 
p0_100 = [max(y_data100)/2, 38, -0.1, max(y_data100), 65, -0.1] 
p0_150 = [max(y_data150)/2, 38, -0.1, max(y_data150), 65, -0.1] 
p0_200 = [max(y_data200)/2, 38, -0.1, max(y_data200), 65, -0.1] 
p0_250 = [max(y_data250)/2, 38, -0.1, max(y_data250), 65, -0.1] 
p0_300 = [max(y_data300)/2, 38, -0.1, max(y_data300), 65, -0.1] 

#4a. x-value plot range fitting:
x_fit = np.linspace(25, 130, 100)

#4b. Fitting operation for calculating theoretical fit parameters: 
def fit_and_evaluate(x_data, y_data, p0):
    popt, _ = curve_fit(compound_sigmoid, x_data, y_data, p0, maxfev=10000)
    y_fit = compound_sigmoid(x_fit, *popt)
    dy_dx_fit = compound_sigmoid_derivative(x_fit, *popt)
    r_squared = r2_score(y_data, compound_sigmoid(x_data, *popt))
    return popt, y_fit, dy_dx_fit, r_squared

popt50, y_fit50, dy_dx_fit50, r_squared50 = fit_and_evaluate(x_data, y_data50, p0_50)
popt100, y_fit100, dy_dx_fit100, r_squared100 = fit_and_evaluate(x_data, y_data100, p0_100)
popt150, y_fit150, dy_dx_fit150, r_squared150 = fit_and_evaluate(x_data, y_data150, p0_150)
popt200, y_fit200, dy_dx_fit200, r_squared200 = fit_and_evaluate(x_data, y_data200, p0_200)
popt250, y_fit250, dy_dx_fit250, r_squared250 = fit_and_evaluate(x_data, y_data250, p0_250)
popt300, y_fit300, dy_dx_fit300, r_squared300 = fit_and_evaluate(x_data, y_data300, p0_300) 

#4c. Print the R^2 Values for each compound sigmoidal fit: 
print(f"R^2 Value: {r_squared50:.3f}")
print(f"R^2 Value: {r_squared100:.3f}")
print(f"R^2 Value: {r_squared150:.3f}")
print(f"R^2 Value: {r_squared200:.3f}")
print(f"R^2 Value: {r_squared250:.3f}")
print(f"R^2 Value: {r_squared300:.3f}")

#5c. Calculation of Sigmoid Fit y-values & derivative y-values using the estimated parameters per   
dataset:
#(i.) [m(w) = 50m(f)]:
y_fit50 = compound_sigmoid(x_fit, *popt50)
dy_dx_fit50 = compound_sigmoid_derivative(x_fit, *popt50)
#(ii.) [m(w) = 100m(f)]:
y_fit100 = compound_sigmoid(x_fit, *popt100)
dy_dx_fit100 = compound_sigmoid_derivative(x_fit, *popt100)
#(iii.) [m(w) = 150m(f)]:
y_fit150 = compound_sigmoid(x_fit, *popt150)
dy_dx_fit150 = compound_sigmoid_derivative(x_fit, *popt150)
#(iv.) [m(w) = 200m(f)]:
y_fit200 = compound_sigmoid(x_fit, *popt200)
dy_dx_fit200 = compound_sigmoid_derivative(x_fit, *popt200)
#(v.) [m(w) = 250m(f)]:
y_fit250 = compound_sigmoid(x_fit, *popt250)
dy_dx_fit250 = compound_sigmoid_derivative(x_fit, *popt250)
#(vi.) [m(w) = 300m(f)]:
y_fit300 = compound_sigmoid(x_fit, *popt300)
dy_dx_fit300 = compound_sigmoid_derivative(x_fit, *popt300)

#6. Calculation of the mean mean, median, and maxima locations for the derivative curves: 
#6a. Input method for weighing the distribution for accurately computed statistical parameters   
considering skew: 
def calculate_mean_location(x, y):
    weighted_sum = np.sum(x * y)
    total_sum = np.sum(y)
    return weighted_sum / total_sum

#6b. Median locations: 
def calculate_median_location(x, y):
    cumulative = np.cumsum(y)
    half_max = cumulative[-1] / 2
    median_index = np.where(cumulative >= half_max)[0][0]
    return x[median_index]

median_locations = [
    calculate_median_location(x_fit, dy_dx_fit50),
    calculate_median_location(x_fit, dy_dx_fit100),
    calculate_median_location(x_fit, dy_dx_fit150),
    calculate_median_location(x_fit, dy_dx_fit200),
    calculate_median_location(x_fit, dy_dx_fit250),
    calculate_median_location(x_fit, dy_dx_fit300)
]

#6c. Mean locations: 
mean_locations = [
    calculate_mean_location(x_fit, dy_dx_fit50),
    calculate_mean_location(x_fit, dy_dx_fit100),
    calculate_mean_location(x_fit, dy_dx_fit150),
    calculate_mean_location(x_fit, dy_dx_fit200),
    calculate_mean_location(x_fit, dy_dx_fit250),
    calculate_mean_location(x_fit, dy_dx_fit300)
]

#6d. Maxima locations: 
max_locations = [
    x_fit[np.argmax(dy_dx_fit50)],
    x_fit[np.argmax(dy_dx_fit100)],
    x_fit[np.argmax(dy_dx_fit150)],
    x_fit[np.argmax(dy_dx_fit200)],
    x_fit[np.argmax(dy_dx_fit250)],
    x_fit[np.argmax(dy_dx_fit300)]
]

#7a. Specifications for computing standard deviation parameters: Local peak isolation for FWHM-ROI  
calculation accuracy. 
def calculate_statistics(x, y):
    std = np.std(y)
    mean = np.mean(y)
    mean_location = np.sum(x * y / np.sum(y))
    median_location = np.median(x[np.argsort(y)[len(y) // 2]])
    return std, mean, mean_location, median_location

def calculate_local_statistics(x, y, peak_indices, window_size=5):
    statistics = []
    for peak in peak_indices:
        start = max(0, peak - window_size)
        end = min(len(x), peak + window_size + 1)
        x_local = x[start:end]
        y_local = y[start:end]
        std, mean, mean_location, median_location = calculate_statistics(x_local, y_local)
        statistics.append((std, mean, mean_location, median_location))
    return statistics

#7b. Finding the local maxima for then calculating: Standard deviation as 'average deviation' of the   
distribution:
#Print the mean, median, and maximum locations along with the specific standard error for the mean:
dy_dx_fits = [dy_dx_fit50, dy_dx_fit100, dy_dx_fit150, dy_dx_fit200, dy_dx_fit250, dy_dx_fit300]

def find_maxima(x, y):
    peaks, _ = find_peaks(y)
    local_maximum = [(x[p], y[p]) for p in peaks] 
    absolute_maximum = (x[np.argmax(y)], np.max(y))
    return local_maximum, absolute_maximum

def calculate_average_deviation_max(x, y, peak_indices, window_size=5):
    statistics = []
    for peak in peak_indices:
        avg_std_max = np.std(y)
        start = max(0, peak - window_size)
        end = min(len(x), peak + window_size + 1)
        x_local = x[start:end]
        y_local = y[start:end]
        average_dev = calculate_average_deviation_max(x_local, y_local)
    return statistics

def calculate_average_deviation_loc(x, y, peak_indices, window_size=5):
    statistics = []
    for peak in peak_indices:
        avg_std_loc = np.std(y)
        start = max(0, peak - window_size)
        end = min(len(x), peak + window_size + 1)
        x_local = x[start:end]
        y_local = y[start:end]
        average_dev = calculate_average_deviation_loc(x_local, y_local)
    return statistics

datasets = [
    ('50m(f)', dy_dx_fit50),
    ('100m(f)', dy_dx_fit100),
    ('150m(f)', dy_dx_fit150),
    ('200m(f)', dy_dx_fit200),
    ('250m(f)', dy_dx_fit250),
    ('300m(f)', dy_dx_fit300)
]

#7c. Calculation of the standard deviation of the mean from its local maximum: 
def calculate_std(x, y):
    mean_loc = np.sum(x * y) / np.sum(y)
    variance = np.sum(y * (x - mean_loc)**2) / np.sum(y)
    return np.sqrt(variance)

def calculate_std_distance(x, y, mean_loc, max_loc):
    std_dev = calculate_std(x, y)
    distance = abs(mean_loc - max_loc)
    std_distance = distance / std_dev
    return std_distance

#8a. Printing statistical parameters and the mean-to-max deviation ('Standard Error'): 
for i, (mean_loc, median_loc, max_loc) in enumerate(zip(mean_locations, median_locations,  
max_locations)):
    std_distance = calculate_std_distance(x_fit, dy_dx_fits[i], mean_loc, max_loc)
    average_dev = calculate_std(x_fit, dy_dx_fits[i])
    print(f"Dataset {i+1}:")
    print(f"Mean Location: {mean_loc:.2f}")
    print(f"Median Location: {median_loc:.2f}")
    print(f"Maximum Location: {max_loc:.2f}") 
    print(f"Standard Error [Mean to Max]: {std_distance:.4f}\n") 

for label, dy_dx_fit in datasets:
    local_maximum, absolute_maximum = find_maxima(x_fit, dy_dx_fit)
    peak_indices = [np.where(x_fit == m[0])[0][0] for m in local_maximum]
    local_stats = calculate_local_statistics(x_fit, dy_dx_fit, peak_indices)

    print(f"\nStatistics for {label}:")


    for i, ((x_peak, y_peak), (avg_std_loc, mean, mean_location, median_location)) in 
enumerate(zip(local_maximum, local_stats)):
        print(f"  Local Maximum {i+1} - X: {x_peak}, Y: {y_peak}, Avg. Std: {avg_std_loc}, Mean:  
{mean}, Mean Location: {mean_location}, Median Location: {median_location}")

    avg_std_max, mean, mean_location, median_location = calculate_statistics(x_fit, dy_dx_fit)
    x_abs_max, y_abs_max = absolute_maximum
    print(f"  Absolute Maximum - X: {x_abs_max}, Y: {y_abs_max}, Avg. Std: {avg_std_max}, Mean:  
{mean}, Mean Location: {mean_location}, Median Location: {median_location}")


#8b. Printing the average deviation for the entire curve: 
for label, dy_dx_fit in datasets:
    local_maximum, absolute_maximum = find_maxima(x_fit, dy_dx_fit)
    peak_indices = [np.where(x_fit == m[0])[0][0] for m in local_maximum]
    local_stats = calculate_local_statistics(x_fit, dy_dx_fit, peak_indices)

    print(f"\nAverage Deviation for {label}:")
    for i, ((x_peak, y_peak), (std, mean, mean_location, median_location)) in 
enumerate(zip(local_maximum, local_stats)):
        print(f" Std: {std}")

#6. Plot data, logistic sigmoid function fit & the derivative: 
plt.figure(figsize = (9, 10))
plt.subplots_adjust(hspace=0.275)

#6a. Datas' & Sigmoid Fits' subplot: 
plt.subplot(2, 1, 1) 

plt.plot(x_data, y_data300, 'o', color = 'black', label = f'$m_{"w"}$ = $300m_{"f"}$' + ' Empirical') 
plt.plot(x_fit, y_fit300, '-', color = 'black', label=f'$m_{"w"}$ = $300m_{"f"}$') # | $R^{"2"}$ 
Value = {r_squared300:.3f}') 
plt.plot(x_fit, y_fit250, '-', color = 'purple', label=f'$m_{"w"}$ = $250m_{"f"}$') # | $R^{"2"}$  
Value = {r_squared300:.3f}') 
plt.plot(x_fit, y_fit200, '-', color = 'red', label=f'$m_{"w"}$ = $200m_{"f"}$') # | $R^{"2"}$ Value 
= {r_squared300:.3f}') 
plt.plot(x_fit, y_fit150, '-', color = 'darkorange', label=f'$m_{"w"}$ = $150m_{"f"}$') # | $R^{"2"}$ 
Value = {r_squared300:.3f}') 
plt.plot(x_fit, y_fit100, '-', color = 'gold', label=f'$m_{"w"}$ = $100m_{"f"}$') # | $R^{"2"}$ Value 
= {r_squared300:.3f}') 
plt.plot(x_fit, y_fit50, '-', color = 'green', label=f'$m_{"w"}$ = $50m_{"f"}$') # | $R^{"2"}$ Value  
= {r_squared300:.3f}') 
plt.plot(x_data, y_data50, 'o', color = 'green', label = f'$m_{"w"}$ = $50m_{"f"}$' + ' Empirical') 

#6b. Work Density Title [Stretched]: USER SELECTION REQUIRED 
plt.xlabel('Temperature (°C)', size = 11, y = 0.97) 
plt.xlim((25, 130)) 
plt.ylim((0, 275)) #- Suitable for presenting Stretched. 
plt.xticks(np.arange(25, 130, 15.0))
plt.ylabel('Work Density ' + '$\it{WD}$ ' + '(Nm/kg)', size = 11, x = 1.02)
plt.title('$\mathregular{[l_{i}}$ = $\mathregular{(1.5)L_{0}]}$ '+ 'Logistic Sigmoidal Fit to Work   
Density Data', size = 13, y = 1.02) 
plt.legend(prop = {"size": 8}) 
plt.grid(True) 

#6b. Strain Title [Stretched]: USER SELECTION REQUIRED 
#plt.xlabel('Temperature (°C)', size = 11, y = 0.97) 
#plt.xlim((25, 130)) 
#plt.ylim((0, 0.8)) #- Suitable for presenting Stretched. 
#plt.xticks(np.arange(25, 130, 15.0))
#plt.ylabel('Strain ' + '$\it{ε}$ ' + '(Nm/kg)', size = 11, x = 1.02)
#plt.title('$\mathregular{[l_{i}}$ = $\mathregular{(1.5)L_{0}]}$ '+ 'Logistic Sigmoidal Fit to Strain  
Data', size = 13, y = 1.02) 
#plt.legend(prop = {"size": 8}) 
#plt.grid(True) 

#6b. Work Density Title [Stretched-Annealed]: USER SELECTION REQUIRED
#plt.xlabel('Temperature (°C)', size = 11, y = 0.97) 
#plt.xlim((25, 130)) 
#plt.ylim((0, 220)) #- Suitable for presenting Stretched-Annealed. 
#plt.xticks(np.arange(25, 130, 15.0))
#plt.ylabel('Work Density ' + '$\it{WD}$ ' + '(Nm/kg)', size = 11, x = 1.02)
#plt.title('$\mathregular{[l_{i}}$ = $\mathregular{(1.5)L_{0}]}$ '+ '$\mathregular{[T_{ann}}$ = 
$\mathregular{(5.2)T_{0}]}$ ' + ' Logistic Sigmoidal Fit to Work Density Data', size = 13, y = 1.02) 
#plt.legend(prop = {"size": 8}) 
#plt.grid(True) 

#6b. Strain Title [Stretched-Annealed]: USER SELECTION REQUIRED
#plt.xlabel('Temperature (°C)', size = 11, y = 0.97) 
#plt.xlim((25, 130)) 
#plt.ylim((0, 0.75)) #- Suitable for presenting Stretched-Annealed. 
#plt.xticks(np.arange(25, 130, 15.0))
#plt.ylabel('Strain ' + '$\it{ε}$ ' + '(Nm/kg)', size = 11, x = 1.02)
#plt.title('$\mathregular{[l_{i}}$ = $\mathregular{(1.5)L_{0}]}$ '+ '$\mathregular{[T_{ann}}$ = 
$\mathregular{(5.2)T_{0}]}$ ' + ' Logistic Sigmoidal Fit to Strain Data', size = 13, y = 1.02) 
#plt.legend(prop = {"size": 8}) 
#plt.grid(True) 

#7. Derivatives of the sigmoids' subplot: 
plt.subplot(2, 1, 2) 

plt.plot(x_fit, dy_dx_fit300, '--', color = 'black', label = f'$m_{"w"}$ = $300m_{"f"}$') 
plt.plot(x_fit, dy_dx_fit250, '--', color = 'purple', label = f'$m_{"w"}$ = $250m_{"f"}$') 
plt.plot(x_fit, dy_dx_fit200, '--', color = 'red', label = f'$m_{"w"}$ = $200m_{"f"}$') 
plt.plot(x_fit, dy_dx_fit150, '--', color = 'darkorange', label = f'$m_{"w"}$ = $150m_{"f"}$') 
plt.plot(x_fit, dy_dx_fit100, '--', color = 'gold', label = f'$m_{"w"}$ = $100m_{"f"}$')
plt.plot(x_fit, dy_dx_fit50, '--', color = 'green', label = f'$m_{"w"}$ = $50m_{"f"}$') 

#7b. Derivative Subplot Work Density Title: USER SELECTION REQUIRED 
plt.xlabel('Temperature (°C)', size = 11, y = 0.97) 
plt.xlim((25, 130)) 
plt.ylim((0, 10)) #- Suitable for presenting Stretched: USER SELECTION REQUIRED. 
#plt.ylim((0, 15)) #- Suitable for presenting Stretched-Annealed: USER SELECTION REQUIRED. 
plt.xticks(np.arange(25, 130, 15.0))
plt.ylabel('Rate of Change ' + '$\it{ΔWD/ΔT}$ ' + '(Nm/kg°C)', size = 11, x = 1.02) 
plt.legend(prop = {"size": 10}) 
plt.title('Work Density Rates (Derivatives of Characteristic Sigmoids)', size = 13, y = 1.02) 
plt.grid(True)

plt.show()
<p>I have multiple datasets of data, for which a compound sigmoid fitting curve is plotted to each. The derivatives of each fit curve are taken, resulting in 6 different theoretical derivative curves. Each derivative curve has two maxima, an absolute maximum at lower temperatures (x-axis) and a local maxima at higher ranges.</p>
<p>Per each of the 6 datasets, I need:</p>
<ol>
<li>the location of the local and absolute maximum (x-axis position) calculated and printed</li>
<li>the location of their means, with respect to each peak/distribution calculated</li>
<li>the location of the medians with respect to each peak/distribution calculated</li>
<li>the average standard deviation for each peak/distribution</li>
<li>Finally, the distance between the means and their respective maxima, converted into a standard deviation value.</li>
</ol>
<p>We would therefore have printed in terminal:</p>
<p>Dataset 1:
Location of absolute maximum:
Mean of absolute maximum:
Median of absolute maximum:
Error of mean from absolute maximum (not distance; the standard deviation equiv. max. - mean location):
Average standard deviation of that maximum:</p>
<p>Location of local maximum:
Mean of local maximum:
Median of local maximum:
Error of mean from local maximum (not distance; the standard deviation equiv. max. - mean location):
Average standard deviation of that maximum:</p>
<p>I have tried extremely hard and ended up creating a Frankenstein code from other types of
specimens I analyzed. So, the terminal plots multiple values.</p>
<pre><code>Often, the standard deviations are too big, as can be visually seen in the PLOTS that are generated. The code can't just print the absolute maximum and local maximum-specific statistics. And can't accurately convert the distance between the means and their respective maxima into standard deviation values (Standard Error). The plots are perfect, please use this as reference.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.signal import find_peaks
from sklearn.metrics import r2_score

#1. User datasets:
x_data = np.array([25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115,
120, 125, 130])
y_data50 = np.array([0, 0.50031, 4.00248, 9.50589, 13.50837, 15.0093, 16.00992, 17.01054, 18.01116,
19.01178, 20.51271, 21.51333, 23.01426, 25.0155, 27.51705, 30.51891, 32.52015, 33.270615, 33.770925,
33.770925, 34.02108, 34.271235])
y_data100 = np.array([0, 4.95405, 12.88053, 20.80701, 28.73349, 36.65997, 39.6324, 42.109425,
44.58645, 47.55888, 50.035905, 52.51293, 55.48536, 59.4486, 63.41184, 68.36589, 72.32913, 73.31994,
74.31075, 75.30156, 75.30156, 75.30156])
y_data150 = np.array([0, 2.96262, 8.88786, 35.55144, 51.84585, 57.77109, 62.21502, 65.17764,
68.14026, 72.58419, 75.54681, 79.250085, 83.694015, 88.8786, 94.063185, 99.24777, 107.394975,
115.54218, 119.98611, 122.94873, 124.43004, 124.43004])
y_data200 = np.array([0, 5.91543, 29.57715, 70.98516, 88.73145, 92.67507, 98.5905, 102.53412,
108.44955, 112.39317, 118.3086, 126.19584, 132.11127, 139.99851, 151.82937, 161.68842, 173.51928,
175.49109, 177.4629, 179.43471, 179.43471, 179.43471])
y_data250 = np.array([0, 12.31155, 39.39696, 83.71854, 105.87933, 115.72857, 120.65319, 128.04012,
135.42705, 142.81398, 152.66322, 162.51246, 169.89939, 182.21094, 194.52249, 206.83404, 219.14559,
221.6079, 224.07021, 226.53252, 226.53252, 226.53252])
y_data300 = np.array([0, 11.81124, 38.38653, 76.77306, 121.06521, 132.87645, 138.78207, 147.6405,
153.54612, 159.45174, 168.31017, 180.12141, 188.97984, 203.74389, 218.50794, 236.2248, 248.03604,
253.94166, 259.84728, 261.323685, 262.80009, 265.7529])

#2a. Sigmoid function definition:
def compound_sigmoid(x, U1, x01, k1, U2, x02, k2):
return U1 / (1 + np.exp(-k1 * (x - x01))) + U2 / (1 + np.exp(-k2 * (x - x02)))

#2b. Derivative of sigmoid function:
def compound_sigmoid_derivative(x, U1, x01, k1, U2, x02, k2):
sig1 = (U1 * k1 * np.exp(-k1 * (x - x01))) / ((1 + np.exp(-k1 * (x - x01))) ** 2)
sig2 = (U2 * k2 * np.exp(-k2 * (x - x02))) / ((1 + np.exp(-k2 * (x - x02))) ** 2)
return sig1 + sig2

#3. Initial guesses for calculating fitting parameters: USER SELECTION REQUIRED
#(1-A.) [PERSONAL/EACH] >WORK DENSITY< Initial Guesses for [STRETCHED]:
p0_50 = [max(y_data50)/2, 38, -0.1, max(y_data50), 65, -0.1]
p0_100 = [max(y_data100)/2, 38, -0.1, max(y_data100), 65, -0.1]
p0_150 = [max(y_data150)/2, 38, -0.1, max(y_data150), 65, -0.1]
p0_200 = [max(y_data200)/2, 38, -0.1, max(y_data200), 65, -0.1]
p0_250 = [max(y_data250)/2, 38, -0.1, max(y_data250), 65, -0.1]
p0_300 = [max(y_data300)/2, 38, -0.1, max(y_data300), 65, -0.1]

#4a. x-value plot range fitting:
x_fit = np.linspace(25, 130, 100)

#4b. Fitting operation for calculating theoretical fit parameters:
def fit_and_evaluate(x_data, y_data, p0):
popt, _ = curve_fit(compound_sigmoid, x_data, y_data, p0, maxfev=10000)
y_fit = compound_sigmoid(x_fit, *popt)
dy_dx_fit = compound_sigmoid_derivative(x_fit, *popt)
r_squared = r2_score(y_data, compound_sigmoid(x_data, *popt))
return popt, y_fit, dy_dx_fit, r_squared

popt50, y_fit50, dy_dx_fit50, r_squared50 = fit_and_evaluate(x_data, y_data50, p0_50)
popt100, y_fit100, dy_dx_fit100, r_squared100 = fit_and_evaluate(x_data, y_data100, p0_100)
popt150, y_fit150, dy_dx_fit150, r_squared150 = fit_and_evaluate(x_data, y_data150, p0_150)
popt200, y_fit200, dy_dx_fit200, r_squared200 = fit_and_evaluate(x_data, y_data200, p0_200)
popt250, y_fit250, dy_dx_fit250, r_squared250 = fit_and_evaluate(x_data, y_data250, p0_250)
popt300, y_fit300, dy_dx_fit300, r_squared300 = fit_and_evaluate(x_data, y_data300, p0_300)

#4c. Print the R^2 Values for each compound sigmoidal fit:
print(f"R^2 Value: {r_squared50:.3f}")
print(f"R^2 Value: {r_squared100:.3f}")
print(f"R^2 Value: {r_squared150:.3f}")
print(f"R^2 Value: {r_squared200:.3f}")
print(f"R^2 Value: {r_squared250:.3f}")
print(f"R^2 Value: {r_squared300:.3f}")

#5c. Calculation of Sigmoid Fit y-values & derivative y-values using the estimated parameters per
dataset:
#(i.) [m(w) = 50m(f)]:
y_fit50 = compound_sigmoid(x_fit, *popt50)
dy_dx_fit50 = compound_sigmoid_derivative(x_fit, *popt50)
#(ii.) [m(w) = 100m(f)]:
y_fit100 = compound_sigmoid(x_fit, *popt100)
dy_dx_fit100 = compound_sigmoid_derivative(x_fit, *popt100)
#(iii.) [m(w) = 150m(f)]:
y_fit150 = compound_sigmoid(x_fit, *popt150)
dy_dx_fit150 = compound_sigmoid_derivative(x_fit, *popt150)
#(iv.) [m(w) = 200m(f)]:
y_fit200 = compound_sigmoid(x_fit, *popt200)
dy_dx_fit200 = compound_sigmoid_derivative(x_fit, *popt200)
#(v.) [m(w) = 250m(f)]:
y_fit250 = compound_sigmoid(x_fit, *popt250)
dy_dx_fit250 = compound_sigmoid_derivative(x_fit, *popt250)
#(vi.) [m(w) = 300m(f)]:
y_fit300 = compound_sigmoid(x_fit, *popt300)
dy_dx_fit300 = compound_sigmoid_derivative(x_fit, *popt300)

#6. Calculation of the mean mean, median, and maxima locations for the derivative curves:
#6a. Input method for weighing the distribution for accurately computed statistical parameters
considering skew:
def calculate_mean_location(x, y):
weighted_sum = np.sum(x * y)
total_sum = np.sum(y)
return weighted_sum / total_sum

#6b. Median locations:
def calculate_median_location(x, y):
cumulative = np.cumsum(y)
half_max = cumulative[-1] / 2
median_index = np.where(cumulative >= half_max)[0][0]
return x[median_index]

median_locations = [
calculate_median_location(x_fit, dy_dx_fit50),
calculate_median_location(x_fit, dy_dx_fit100),
calculate_median_location(x_fit, dy_dx_fit150),
calculate_median_location(x_fit, dy_dx_fit200),
calculate_median_location(x_fit, dy_dx_fit250),
calculate_median_location(x_fit, dy_dx_fit300)
]

#6c. Mean locations:
mean_locations = [
calculate_mean_location(x_fit, dy_dx_fit50),
calculate_mean_location(x_fit, dy_dx_fit100),
calculate_mean_location(x_fit, dy_dx_fit150),
calculate_mean_location(x_fit, dy_dx_fit200),
calculate_mean_location(x_fit, dy_dx_fit250),
calculate_mean_location(x_fit, dy_dx_fit300)
]

#6d. Maxima locations:
max_locations = [
x_fit[np.argmax(dy_dx_fit50)],
x_fit[np.argmax(dy_dx_fit100)],
x_fit[np.argmax(dy_dx_fit150)],
x_fit[np.argmax(dy_dx_fit200)],
x_fit[np.argmax(dy_dx_fit250)],
x_fit[np.argmax(dy_dx_fit300)]
]

#7a. Specifications for computing standard deviation parameters: Local peak isolation for FWHM-ROI
calculation accuracy.
def calculate_statistics(x, y):
std = np.std(y)
mean = np.mean(y)
mean_location = np.sum(x * y / np.sum(y))
median_location = np.median(x[np.argsort(y)[len(y) // 2]])
return std, mean, mean_location, median_location

def calculate_local_statistics(x, y, peak_indices, window_size=5):
statistics = []
for peak in peak_indices:
start = max(0, peak - window_size)
end = min(len(x), peak + window_size + 1)
x_local = x[start:end]
y_local = y[start:end]
std, mean, mean_location, median_location = calculate_statistics(x_local, y_local)
statistics.append((std, mean, mean_location, median_location))
return statistics

#7b. Finding the local maxima for then calculating: Standard deviation as 'average deviation' of the
distribution:
#Print the mean, median, and maximum locations along with the specific standard error for the mean:
dy_dx_fits = [dy_dx_fit50, dy_dx_fit100, dy_dx_fit150, dy_dx_fit200, dy_dx_fit250, dy_dx_fit300]

def find_maxima(x, y):
peaks, _ = find_peaks(y)
local_maximum = [(x[p], y[p]) for p in peaks]
absolute_maximum = (x[np.argmax(y)], np.max(y))
return local_maximum, absolute_maximum

def calculate_average_deviation_max(x, y, peak_indices, window_size=5):
statistics = []
for peak in peak_indices:
avg_std_max = np.std(y)
start = max(0, peak - window_size)
end = min(len(x), peak + window_size + 1)
x_local = x[start:end]
y_local = y[start:end]
average_dev = calculate_average_deviation_max(x_local, y_local)
return statistics

def calculate_average_deviation_loc(x, y, peak_indices, window_size=5):
statistics = []
for peak in peak_indices:
avg_std_loc = np.std(y)
start = max(0, peak - window_size)
end = min(len(x), peak + window_size + 1)
x_local = x[start:end]
y_local = y[start:end]
average_dev = calculate_average_deviation_loc(x_local, y_local)
return statistics

datasets = [
('50m(f)', dy_dx_fit50),
('100m(f)', dy_dx_fit100),
('150m(f)', dy_dx_fit150),
('200m(f)', dy_dx_fit200),
('250m(f)', dy_dx_fit250),
('300m(f)', dy_dx_fit300)
]

#7c. Calculation of the standard deviation of the mean from its local maximum:
def calculate_std(x, y):
mean_loc = np.sum(x * y) / np.sum(y)
variance = np.sum(y * (x - mean_loc)**2) / np.sum(y)
return np.sqrt(variance)

def calculate_std_distance(x, y, mean_loc, max_loc):
std_dev = calculate_std(x, y)
distance = abs(mean_loc - max_loc)
std_distance = distance / std_dev
return std_distance

#8a. Printing statistical parameters and the mean-to-max deviation ('Standard Error'):
for i, (mean_loc, median_loc, max_loc) in enumerate(zip(mean_locations, median_locations,
max_locations)):
std_distance = calculate_std_distance(x_fit, dy_dx_fits, mean_loc, max_loc)
average_dev = calculate_std(x_fit, dy_dx_fits)
print(f"Dataset {i+1}:")
print(f"Mean Location: {mean_loc:.2f}")
print(f"Median Location: {median_loc:.2f}")
print(f"Maximum Location: {max_loc:.2f}")
print(f"Standard Error [Mean to Max]: {std_distance:.4f}\n")

for label, dy_dx_fit in datasets:
local_maximum, absolute_maximum = find_maxima(x_fit, dy_dx_fit)
peak_indices = [np.where(x_fit == m[0])[0][0] for m in local_maximum]
local_stats = calculate_local_statistics(x_fit, dy_dx_fit, peak_indices)

print(f"\nStatistics for {label}:")


for i, ((x_peak, y_peak), (avg_std_loc, mean, mean_location, median_location)) in
enumerate(zip(local_maximum, local_stats)):
print(f" Local Maximum {i+1} - X: {x_peak}, Y: {y_peak}, Avg. Std: {avg_std_loc}, Mean:
{mean}, Mean Location: {mean_location}, Median Location: {median_location}")

avg_std_max, mean, mean_location, median_location = calculate_statistics(x_fit, dy_dx_fit)
x_abs_max, y_abs_max = absolute_maximum
print(f" Absolute Maximum - X: {x_abs_max}, Y: {y_abs_max}, Avg. Std: {avg_std_max}, Mean:
{mean}, Mean Location: {mean_location}, Median Location: {median_location}")


#8b. Printing the average deviation for the entire curve:
for label, dy_dx_fit in datasets:
local_maximum, absolute_maximum = find_maxima(x_fit, dy_dx_fit)
peak_indices = [np.where(x_fit == m[0])[0][0] for m in local_maximum]
local_stats = calculate_local_statistics(x_fit, dy_dx_fit, peak_indices)

print(f"\nAverage Deviation for {label}:")
for i, ((x_peak, y_peak), (std, mean, mean_location, median_location)) in
enumerate(zip(local_maximum, local_stats)):
print(f" Std: {std}")

#6. Plot data, logistic sigmoid function fit & the derivative:
plt.figure(figsize = (9, 10))
plt.subplots_adjust(hspace=0.275)

#6a. Datas' & Sigmoid Fits' subplot:
plt.subplot(2, 1, 1)

plt.plot(x_data, y_data300, 'o', color = 'black', label = f'$m_{"w"}$ = $300m_{"f"}$' + ' Empirical')
plt.plot(x_fit, y_fit300, '-', color = 'black', label=f'$m_{"w"}$ = $300m_{"f"}$') # | $R^{"2"}$
Value = {r_squared300:.3f}')
plt.plot(x_fit, y_fit250, '-', color = 'purple', label=f'$m_{"w"}$ = $250m_{"f"}$') # | $R^{"2"}$
Value = {r_squared300:.3f}')
plt.plot(x_fit, y_fit200, '-', color = 'red', label=f'$m_{"w"}$ = $200m_{"f"}$') # | $R^{"2"}$ Value
= {r_squared300:.3f}')
plt.plot(x_fit, y_fit150, '-', color = 'darkorange', label=f'$m_{"w"}$ = $150m_{"f"}$') # | $R^{"2"}$
Value = {r_squared300:.3f}')
plt.plot(x_fit, y_fit100, '-', color = 'gold', label=f'$m_{"w"}$ = $100m_{"f"}$') # | $R^{"2"}$ Value
= {r_squared300:.3f}')
plt.plot(x_fit, y_fit50, '-', color = 'green', label=f'$m_{"w"}$ = $50m_{"f"}$') # | $R^{"2"}$ Value
= {r_squared300:.3f}')
plt.plot(x_data, y_data50, 'o', color = 'green', label = f'$m_{"w"}$ = $50m_{"f"}$' + ' Empirical')

#6b. Work Density Title [Stretched]: USER SELECTION REQUIRED
plt.xlabel('Temperature (°C)', size = 11, y = 0.97)
plt.xlim((25, 130))
plt.ylim((0, 275)) #- Suitable for presenting Stretched.
plt.xticks(np.arange(25, 130, 15.0))
plt.ylabel('Work Density ' + '$\it{WD}$ ' + '(Nm/kg)', size = 11, x = 1.02)
plt.title('$\mathregular{[l_{i}}$ = $\mathregular{(1.5)L_{0}]}$ '+ 'Logistic Sigmoidal Fit to Work
Density Data', size = 13, y = 1.02)
plt.legend(prop = {"size": 8})
plt.grid(True)

#6b. Strain Title [Stretched]: USER SELECTION REQUIRED
#plt.xlabel('Temperature (°C)', size = 11, y = 0.97)
#plt.xlim((25, 130))
#plt.ylim((0, 0.8)) #- Suitable for presenting Stretched.
#plt.xticks(np.arange(25, 130, 15.0))
#plt.ylabel('Strain ' + '$\it{ε}$ ' + '(Nm/kg)', size = 11, x = 1.02)
#plt.title('$\mathregular{[l_{i}}$ = $\mathregular{(1.5)L_{0}]}$ '+ 'Logistic Sigmoidal Fit to Strain
Data', size = 13, y = 1.02)
#plt.legend(prop = {"size": 8})
#plt.grid(True)

#6b. Work Density Title [Stretched-Annealed]: USER SELECTION REQUIRED
#plt.xlabel('Temperature (°C)', size = 11, y = 0.97)
#plt.xlim((25, 130))
#plt.ylim((0, 220)) #- Suitable for presenting Stretched-Annealed.
#plt.xticks(np.arange(25, 130, 15.0))
#plt.ylabel('Work Density ' + '$\it{WD}$ ' + '(Nm/kg)', size = 11, x = 1.02)
#plt.title('$\mathregular{[l_{i}}$ = $\mathregular{(1.5)L_{0}]}$ '+ '$\mathregular{[T_{ann}}$ =
$\mathregular{(5.2)T_{0}]}$ ' + ' Logistic Sigmoidal Fit to Work Density Data', size = 13, y = 1.02)
#plt.legend(prop = {"size": 8})
#plt.grid(True)

#6b. Strain Title [Stretched-Annealed]: USER SELECTION REQUIRED
#plt.xlabel('Temperature (°C)', size = 11, y = 0.97)
#plt.xlim((25, 130))
#plt.ylim((0, 0.75)) #- Suitable for presenting Stretched-Annealed.
#plt.xticks(np.arange(25, 130, 15.0))
#plt.ylabel('Strain ' + '$\it{ε}$ ' + '(Nm/kg)', size = 11, x = 1.02)
#plt.title('$\mathregular{[l_{i}}$ = $\mathregular{(1.5)L_{0}]}$ '+ '$\mathregular{[T_{ann}}$ =
$\mathregular{(5.2)T_{0}]}$ ' + ' Logistic Sigmoidal Fit to Strain Data', size = 13, y = 1.02)
#plt.legend(prop = {"size": 8})
#plt.grid(True)

#7. Derivatives of the sigmoids' subplot:
plt.subplot(2, 1, 2)

plt.plot(x_fit, dy_dx_fit300, '--', color = 'black', label = f'$m_{"w"}$ = $300m_{"f"}$')
plt.plot(x_fit, dy_dx_fit250, '--', color = 'purple', label = f'$m_{"w"}$ = $250m_{"f"}$')
plt.plot(x_fit, dy_dx_fit200, '--', color = 'red', label = f'$m_{"w"}$ = $200m_{"f"}$')
plt.plot(x_fit, dy_dx_fit150, '--', color = 'darkorange', label = f'$m_{"w"}$ = $150m_{"f"}$')
plt.plot(x_fit, dy_dx_fit100, '--', color = 'gold', label = f'$m_{"w"}$ = $100m_{"f"}$')
plt.plot(x_fit, dy_dx_fit50, '--', color = 'green', label = f'$m_{"w"}$ = $50m_{"f"}$')

#7b. Derivative Subplot Work Density Title: USER SELECTION REQUIRED
plt.xlabel('Temperature (°C)', size = 11, y = 0.97)
plt.xlim((25, 130))
plt.ylim((0, 10)) #- Suitable for presenting Stretched: USER SELECTION REQUIRED.
#plt.ylim((0, 15)) #- Suitable for presenting Stretched-Annealed: USER SELECTION REQUIRED.
plt.xticks(np.arange(25, 130, 15.0))
plt.ylabel('Rate of Change ' + '$\it{ΔWD/ΔT}$ ' + '(Nm/kg°C)', size = 11, x = 1.02)
plt.legend(prop = {"size": 10})
plt.title('Work Density Rates (Derivatives of Characteristic Sigmoids)', size = 13, y = 1.02)
plt.grid(True)

plt.show()
</code></pre>
 
Top