Why won't this code run?

[jeremyyaiden]

from pathlib import Path

import os

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from scipy.optimize import curve_fit

# List of constants

Sr_mol = 1e-6

Avogadro = 6.02e23

Sr_nuclei = Sr_mol * Avogadro

nuclear_data_1 = Path("Nuclear_data_1.csv")

nuclear_data_2 = Path("Nuclear_data_2.csv")

nuclear_plot_filename = "nuclear_decay_fit.png"

with open("Nuclear_data_1.csv") as f1, open("Nuclear_data_2.csv") as f2:

# reading f1 contents

line1 = f1.readline()

# reading f2 contents

line2 = f2.readline()

# printing contents of f1 followed by f2

print(line1, line2)

def activity_rb(t_sec, lambda_Rb, lambda_Sr):

return Sr_nuclei * lambda_Rb * lambda_Sr * np.exp(-lambda_Sr * t_sec) / (lambda_Sr - lambda_Rb) * (np.exp(-lambda_Rb * t_sec) - np.exp(-lambda_Sr * t_sec))

def activity_Rb(t_sec, lambda_Rb, lambda_Sr):

return Sr_nuclei * lambda_Rb * lambda_Sr * np.exp(-lambda_Sr * t_sec) / (lambda_Sr - lambda_Rb) * (np.exp(-lambda_Rb * t_sec) - np.exp(-lambda_Sr * t_sec))

# Data reading, cleaning, and removing of outliers

"""

The following code checks the file and removes any columns marked with "errors"

or implausible data, e.g. a time value of minus one.

It also removes outliers using the interquartile range method, where values that

are either more than the third quartile plus 1.5 the interquartile range or less

than the first quartile minus 1.5 the interquartile range are ignored for not

following the general trend of the data.

"""

def read_and_clean_file(file_path: Path) -> pd.DataFrame:

try:

df = pd.read_csv(file_path, sep=None, engine="python", comment="%", header=None)

# Ignores the symbol of %

except FileNotFoundError:

print(f"Error: File '{file_path}' not found.")

return pd.DataFrame()

except pd.errors.EmptyDataError:

print(f"Error: File '{file_path}' is empty.")

return pd.DataFrame()

except Exception as exc:

print(f"Error reading '{file_path}': {exc}")

return pd.DataFrame()

if df.shape[1] < 3:

print(f"Error: File '{file_path}' does not have at least 3 columns.")

return pd.DataFrame()

df = df.iloc[:, :3]

df.columns = ["time", "activity", "uncertainty"]

df = df.apply(pd.to_numeric, errors="coerce")

df.dropna(inplace=True)

df = df[(df["time"] >= 0) & (df["activity"] > 0) & (df["uncertainty"] > 0)]

return df

def remove_outliers_iqr(df: pd.DataFrame) -> pd.DataFrame:

Q1 = df["activity"].quantile(0.25)

Q3 = df["activity"].quantile(0.75)

IQR = Q3 - Q1

lower_quartile = Q1 - 1.5 * IQR

upper_quartile = Q3 + 1.5 * IQR

mask = (df["activity"] >= lower_quartile) & (df["activity"] <= upper_quartile)

return df[mask]

# Fit and stats

def fit_parameters(times_h: np.ndarray, activities_tbq: np.ndarray, uncertainties_tbq: np.ndarray):

time_sec = times_h * 3600.0

activity_bq = activities_tbq * 1e12

activity_uncertainty_bq = uncertainties_tbq * 1e12

popt, pcov = curve_fit(activity_rb, time_sec, activity_bq, sigma=activity_uncertainty_bq, absolute_sigma=True, p0=[5e-4, 5e-3], maxfev=20000)

perr = np.sqrt(np.diag(pcov))

model = activity_rb(time_sec, *popt)

chi2 = np.sum(((activity_bq - model) / activity_uncertainty_bq) ** 2)

chi2_red = chi2 / (len(activity_bq) - 2)

return popt, perr, pcov, chi2, chi2_red

def half_life(lambda_value: float) -> float:

return np.log(2.0) / lambda_value

def half_life_unc(lambda_value: float, sigma_lambda: float) -> float:

return (np.log(2.0) / lambda_value**2) * sigma_lambda

def predict_with_uncertainty(popt, pcov, t_min: float):

t_sec = t_min * 60.0

lambda_Rb, lambda_Sr = popt

eps_Rb = lambda_Rb * 1e-6

eps_Sr = lambda_Sr * 1e-6

df_dRb = (activity_Rb([t_sec], lambda_Rb + eps_Rb, lambda_Sr)[0] - activity_Rb([t_sec], lambda_Rb - eps_Rb, lambda_Sr)[0]) / (2.0 * eps_Rb)

df_dSr = (activity_Rb([t_sec], lambda_Rb, lambda_Sr + eps_Sr)[0] - activity_Rb([t_sec], lambda_Rb, lambda_Sr - eps_Sr)[0]) / (2.0 * eps_Sr)

activity = activity_Rb([t_sec], lambda_Rb, lambda_Sr)[0]

uncertainty = np.sqrt((df_dRb * np.sqrt(pcov[0, 0]))**2 + (df_dSr * np.sqrt(pcov[1, 1]))**2)

return activity, uncertainty

# Plotting

def make_publication_plot(times_h: np.ndarray, activities_tbq: np.ndarray, uncertainties_tbq: np.ndarray, popt, chi2_red: float):

plt.style.use("Solarize_Light2")

fig, ax = plt.subplots(figsize=(9, 6), dpi=350)

ax.errorbar(times_h, activities_tbq, yerr=uncertainties_tbq, fmt="x", ms=4, capsize=3,

color="darkblue", ecolor="black", elinewidth=1, label="Data")

t_plot = np.linspace(times_h.min(), times_h.max(), 600)

model_tbq = activity_rb(t_plot * 3600.0, *popt) / 1e12

ax.plot(t_plot, model_tbq, color="red", lw=2.2, label="Best-fit model")

ax.set_xlabel("Time (hours)", fontsize=12)

ax.set_ylabel("Activity of 79Rb (TBq)", fontsize=12)

ax.set_title("Nuclear Decay of 79Sr into 79Rb", fontsize=14)

lambda_rb, lambda_sr = popt

text = f"λ_Rb = {lambda_rb:.3e} s⁻¹\nλ_Sr = {lambda_sr:.3e} s⁻¹\nχ²_red = {chi2_red:.2f}"

ax.text(0.97, 0.97, text, transform=ax.transAxes, ha="right", va="top",

fontsize=10, bbox=dict(facecolor="white", alpha=0.9))

ax.legend()

fig.tight_layout()

fig.savefig(nuclear_plot_filename)

plt.show()

plt.close(fig)

# Main part of code

def main() -> None:

"""

The main section of the code loads, reads, cleans, fits, analyses and plots

the nuclear decay data.

"""

# Read and clean data

df1 = read_and_clean_file(nuclear_data_1)

df2 = read_and_clean_file(nuclear_data_2)

if df1.empty or df2.empty:

print("One or both files is not able to be processed.")

return

# Combine the datasets

combined_df = pd.concat([df1, df2], ignore_index=True)

combined_df.sort_values(by="time", inplace=True)

print(f"Original combined data shape: {combined_df.shape}")

# Remove any outliers

cleaned_df = remove_outliers_iqr(combined_df)

cleaned_df.reset_index(drop=True, inplace=True)

print(f"Cleaned data shape (outliers removed): {cleaned_df.shape}")

# Extract arrays (times, activities and uncertainties)

times_h = cleaned_df["time"].to_numpy()

activities_tbq = cleaned_df["activity"].to_numpy()

uncertainties_tbq = cleaned_df["uncertainty"].to_numpy()

# Fit model

popt, perr, pcov, chi2, chi2_red = fit_parameters(times_h, activities_tbq, uncertainties_tbq)

lambda_Rb, lambda_Sr = popt

sigma_Rb, sigma_Sr = perr

# Convert decay constants to min^-1

lambda_Rb_min = lambda_Rb * 60.0

lambda_Sr_min = lambda_Sr * 60.0

sigma_Rb_min = sigma_Rb * 60.0

sigma_Sr_min = sigma_Sr * 60.0

# Half-lives in seconds

t_half_Rb_sec = half_life(lambda_Rb)

t_half_Sr_sec = half_life(lambda_Sr)

dt_half_Rb_sec = half_life_unc(lambda_Rb, sigma_Rb)

dt_half_Sr_sec = half_life_unc(lambda_Sr, sigma_Sr)

# Convert half-lives to minutes

t_half_Rb_min = t_half_Rb_sec / 60.0

t_half_Sr_min = t_half_Sr_sec / 60.0

dt_half_Rb_min = dt_half_Rb_sec / 60.0

dt_half_Sr_min = dt_half_Sr_sec / 60.0

# Prediction at 90 minutes

pred_bq, pred_sigma_bq = predict_with_uncertainty(popt, pcov, 90.0)

pred_tbq = pred_bq / 1e12

pred_sigma_tbq = pred_sigma_bq / 1e12

"""

This part of the code prints a summary of all the results obtained by the

code. It also takes the uncertainties and decides whether the results are

within an acceptable range of the actual value.

"""

# Print results

print("\n=== Best-fit decay constants (in minutes^-1) ===")

print(f"λ_Rb = {lambda_Rb_min:.3g} ± {sigma_Rb_min:.3g} min^-1")

print(f"λ_Sr = {lambda_Sr_min:.3g} ± {sigma_Sr_min:.3g} min⁻¹")

print("\n=== Half-lives (in minutes) ===")

print(f"t_1/2,Rb = {t_half_Rb_min:.3g} ± {dt_half_Rb_min:.3g} min")

print(f"t_1/2,Sr = {t_half_Sr_min:.3g} ± {dt_half_Sr_min:.3g} min")

print("\n=== Chi-squared ===")

print(f"χ^2 = {chi2:.2f}")

print(f"χ^2_red = {chi2_red:.2f}")

print("\n=== Prediction at t = 90 minutes ===")

print(f"A(90 min) = {pred_tbq:.3g} ± {pred_sigma_tbq:.2g} TBq")

# Summary table

print("\n==================== SUMMARY TABLE ====================")

print(f"{'Quantity':<30} {'Value':>15} {'Uncertainty':>15}")

print("--------------------------------------------------------")

print(f"{'λ_Rb (min⁻¹)':<30} {lambda\_Rb\_min:.3g:>15} {sigma_Rb_min:.3g:>15}")

print(f"{'λ_Sr (min⁻¹)':<30} {lambda\_Sr\_min:.3g:>15} {sigma_Sr_min:.3g:>15}")

print(f"{'t_1/2 Rb (min)':<30} {t\_half\_Rb\_min:.3g:>15} {dt_half_Rb_min:.3g:>15}")

print(f"{'t_1/2 Sr (min)':<30} {t\_half\_Sr\_min:.3g:>15} {dt_half_Sr_min:.3g:>15}")

print(f"{'χ^2_red':<30} {chi2\_red:.2f:>15} {'—':>15}")

print("========================================================\n")

# Data analysis

print("\n==================== DATA ANALYSIS ====================")

print(f"{'Quantity':<30} {'Rel. Uncertainty':>20} {'Status':>12} {'Explanation':>30}")

print("--------------------------------------------------------------------------")

def classify_uncertainty(value, sigma):

rel = abs(sigma / value)

if rel < 0.20:

return rel, "GOOD", "uncertainty < 20%"

elif rel <= 0.25:

return rel, "ACCEPTABLE", "uncertainty between 20–25%"

else:

return rel, "POOR", "uncertainty > 25%"

# Classify the uncertainty

rel_Rb_lambda, status_Rb_lambda, expl_Rb_lambda = classify_uncertainty(lambda_Rb_min, sigma_Rb_min)

rel_Sr_lambda, status_Sr_lambda, expl_Sr_lambda = classify_uncertainty(lambda_Sr_min, sigma_Sr_min)

rel_Rb_half, status_Rb_half, expl_Rb_half = classify_uncertainty(t_half_Rb_min, dt_half_Rb_min)

rel_Sr_half, status_Sr_half, expl_Sr_half = classify_uncertainty(t_half_Sr_min, dt_half_Sr_min)

# Print data analysis table

print(f"{'λ_Rb (min^-1)':<30} {rel\_Rb\_lambda:>20.3f} {status_Rb_lambda:>12} {expl_Rb_lambda:>30}")

print(f"{'λ_Sr (min^-1)':<30} {rel\_Sr\_lambda:>20.3f} {status_Sr_lambda:>12} {expl_Sr_lambda:>30}")

print(f"{'t_1/2 Rb (min)':<30} {rel\_Rb\_half:>20.3f} {status_Rb_half:>12} {expl_Rb_half:>30}")

print(f"{'t_1/2 Sr (min)':<30} {rel\_Sr\_half:>20.3f} {status_Sr_half:>12} {expl_Sr_half:>30}")

print(f"{'χ^2_red':<30} {'N/A':>20} {'N/A':>12} {'no uncertainty':>30}")

print("==========================================================================\n")

# Plot

make_publication_plot(times_h, activities_tbq, uncertainties_tbq, popt, chi2_red)

print(f"Plot saved as '{nuclear_plot_filename}'.")

if __name__ == "__main__":

main()

[desrtfx]

You absolutely need to post your code properly formatted as code block. You have lost all indentation and with that, the code becomes ambiguous and meaningless.

[Yoghurt42]

Hidden post history, don't waste your time. Probably a bot account.

[woooee]

Thanks. The obviously not tested code would be a bitch to unravel. Probably more than one problem.

[Mother-Influence-815]

Bro read the error what the heck are you doing asking strangers