Optimization on MSE Not Working Correctly

Hello. I am trying to optimize the MSE (mean squared error) between experimental data and theoretical data. I want to constrain my parameter values within their bounds and to the constraint function, calculated the optimal delay time that gives the minimal MSE. I tried an initial guess that matches another block of code I have (which is correct, the one following the #), but when trying the initial guess I have written out (which I can verify gives larger theoretical data and thus larger/worse MSE), it returns the actual theoretical value that I have from another block of code and just the initial guess parameters as the optimal ones. Any advice? Could the blocks of code be interfering with one another? Should I open a new notebook? Or is something fundamentally wrong?

import numpy as np
import scipy.optimize
from scipy.optimize import LinearConstraint
from scipy.optimize import Bounds
from scipy.optimize import minimize

lab_delayvalues_glu = {2: 2.58, 
                       4: 3.84, 
                       5: 5.68, 
                       6: 4.43, 
                       8: 3.60, 
                       10: 3.10}

def function(z, N):
    z = k_off, P_in, P_out
    
    # Set initial conditions
    y0 = [0, r_min+(lambda_0/kappa_t), 0, sigma_i]

    # Define the time range
    tMax=(N+1)+4+(30/lambda_0)
    time=np.linspace(0,tMax,int(10*tMax))

    growth_rates_list = []

    # Solve the ODEs for each value of N
    solution = solve_ivp(model_equations, [time[0], time[-1]], y0, method='LSODA', t_eval=time)

    gr1_values = np.zeros_like(time)
    gr2_values = np.zeros_like(time)

    # Calculate gr1 and gr2 for all time points
    for i, t in enumerate(time):
        if 0 < t < N + 1:
            gr1_values[i] = ((solution.y[1, i] - r_min) * kappa_t)
        else:
            gr2_values[i] = ((r_max - solution.y[1, i] - solution.y[2, i]) * kappa_n)

    # Calculate the total growth rate by element-wise addition
    growth_rate = gr1_values + gr2_values

    # Append growth rate to the list
    growth_rates_list.append(growth_rate)
    
    ODs_and_x_intercepts = []
    data_list = []

    # Calculate OD and x-intercept data for each N
    OD = np.cumsum(growth_rate) * (time[1] - time[0])

    start_x = time[-1]
    start_index = np.where(time == start_x)[0][0]
    start_y = OD[start_index]
    
    if abs(1-(OD[start_index]-OD[start_index - 1])/(time[1]-time[0])/lambda_0) > 0.01:
        print(N, lambda_0)

    desired_slope = lambda_0

    x_intercept = start_x - (start_y / desired_slope)

    return x_intercept

# Define the Mean Squared Error (MSE) objective function for a specific N
def mse_objective(params, N):
    # Get the experimental delay time for the current N
    exp_delay_time = lab_delayvalues_glu[N]
    
    # Calculate delay time using current parameters
    calc_delay_time = function(params, N)
    
    # Calculate squared error and return MSE
    mse = (exp_delay_time - calc_delay_time) ** 2
    return mse

# Define the constraint function
def constraint_function(params):
    # Calculate the value of the constraint (IC_50 - 1.75)
    k_off, P_in, P_out = params
    calc_IC_50 = (k_off/k_on)*(P_out/P_in)*(lambda_max/lambda_0)
    return abs(calc_IC_50 - 1.75)  # Ensure IC_50 is greater than or equal to 1.75

# Initial guess for the parameters
initial_guess = [513.0, 1360.0, 80.0] #[4000.0, 1540.8822857142854, 100.0]

# Define bounds for the parameters
param_bounds = [(513, 20520), (1360, 2040), (80, 120)]  # Example bounds for k_off, P_in, P_out

# Define the constraint dictionary
constraints = {'type': 'ineq', 'fun': constraint_function}

# Dictionary to store experimental and optimized delay times for plotting
experimental_delay_times = []
optimized_delay_times = []

# Iterate over all N values in lab_delayvalues_glu and minimize MSE for each
for N in lab_delayvalues_glu:
    # Experimental delay time for current N
    experimental_delay_time = lab_delayvalues_glu[N]
    experimental_delay_times.append(experimental_delay_time)
    
    # Minimize MSE for current N
    result = minimize(mse_objective, initial_guess, args=(N,), method='trust-constr', bounds=param_bounds, constraints=constraints)
    
    # Optimized delay time for current N
    optimized_delay_time = function(result.x, N)
    optimized_delay_times.append(optimized_delay_time)
    
    # Print the results for current N
    print(f"N={N}:")
    print("Experimental Delay Time:", experimental_delay_time)
    print("Optimized Delay Time:", optimized_delay_time)
    print("Optimized Parameters [k_off, P_in, P_out]:", result.x)
    print("Mean Squared Error (MSE):", result.fun)
    print()

# Plot experimental and calculated delay times for each N
plt.figure(figsize=(8, 6))
plt.plot(list(lab_delayvalues_glu.keys()), experimental_delay_times, marker='o', color='blue', label='Experimental')
plt.plot(list(lab_delayvalues_glu.keys()), optimized_delay_times, marker='o', color='red', label='Calculated')

plt.title('Experimental vs Calculated Delay Times')
plt.xlabel('N Value')
plt.ylabel('Delay Time')
plt.xticks(list(lab_delayvalues_glu.keys()))  # Set x-axis ticks to N values
plt.legend()
plt.grid(True)
plt.show()

Hi,

there are a few errors in your code; I won’t name them all. Here, I will only highlight some issues with the first function in your code. The same principles apply to the rest of the code.

In your first function (def function(z, N):), you have z as in input parameter, yet it it is not used as such in the body. It is used as a dependent variable of the function and assigned to values of three variables which are not defined / assigned anywhere in your code.

z = k_off, P_in, P_out 

The following variables are also not defined anywhere in your code:

r_min, lambda_0, kappa_t, sigma_i, tMax, 

All variables must be assigned (defined) before they may be used.

I think that you need to review the fundamentals of how functions work, and how they are defined.

Where is the function solve_ivp defined? It is not imported. I take it is to solve ordinary differential equations?

1 Like

The error that @onePythonUser pointed out is a really fundamental error. In a Jupyter Notebook this error will be more pernicious than usual, since you might have been using/definining global variables called k_off, P_in, P_out, so the code would not throw an exception but just give invalid results.

You can also make your code more readable (which implies easier to debug/maintain/understand) by choosing more descriptive (and non-redundant) names for functions. For instance function is about the worst function name to choose - since it’s both non-descriptive and redundant. Other example: constraint_function is somewhat decriptive but is ambiguous and contains a redundant suffix _function. It’s ambiguous because it’s unclear from the name whether it is intended to apply a constraint or to compute one. So, it seems better to leave out the suffix, and to change the name to sth more specific like calculate_constraint.

1 Like

Hello, thank you for the responses. I have since rectified my code by using Monte Carlo simulations, however the run time is incredibly long. I’d like to put a stronger constraint on my parameters and have something like the following. However, it is not returning any optimal parameters which I find incorrect. Any advice?

import numpy as np
import scipy.optimize
from scipy.optimize import LinearConstraint
from scipy.optimize import Bounds
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import math
from timeit import default_timer as timer

start = timer()

# Define parameter values
r_min = 5.8
r_max = 52.3
P_in = 1485.10542448074
P_out = 100.0
k_on = 1026.0
k_off = 2000.0
K_D = k_off/k_on
gamma = 11.0
lambda_0 = 0.40
IC_50 = 1.45
delta_r = 46.5 #r_max - r_min
rho = 116.25 #delta_r/delta_phi_r uM
kappa_t = gamma/rho
kappa_n = 1/((delta_r * ((1 / lambda_0) - (1 / (kappa_t * delta_r)))))
phi_r_min = 0.05 #r_min/rho
phi_r_max = 0.45 #r_max/rho
nu = rho*kappa_n #h^-1
chi_r_min = phi_r_max-(lambda_0/nu)
sigma_i = lambda_0/chi_r_min
lambda_0_star = 2*math.sqrt(P_out*kappa_t*K_D)
IC_50_star = (delta_r*lambda_0_star)/(2*P_in)
lambda_max = kappa_t*delta_r
N_values = [4, 6, 8, 10]

exp_delay_times = {4: 2.22, 6: 2.71, 8: 2.86, 10: 3.03}

def model_equations(t, y, N, z):
    a = y[0] # Antibiotic concentration
    r_u = y[1] # Unbound ribosomes
    r_b = y[2] # Bound ribosomes
    sigma = y[3] # Translational activity

    # Iterate step-pulse addition of antibiotics
    if 0<=t<=1:
        a_ex = 0
        chi_r = chi_r_min
        x = (r_u-r_min)*kappa_t
    elif 1<=t<=N+1:
        a_ex = (4*IC_50)/N
        lam_f = lambda_0/(1+(4/N))
        chi_r_max = phi_r_max-(lam_f/nu)
        chi_r = chi_r_max
        x = (r_u-r_min)*kappa_t 
    else:
        a_ex = 0
        chi_r = chi_r_min
        x = (r_max-r_u-r_b)*kappa_n
    
    # Calculate synthesis rate
    #s = x * (r_max - ((x * delta_r) * ((1 / lambda_0) - (1 / (kappa_t*delta_r)))))
    #s = x*rho*(x/sigma)
    s = x*rho*(phi_r_max/(1+(sigma/nu)))
    
    # Calculate F - binding/unbinding
    F = (k_on * a * (r_u - r_min)) - (z[0] * r_b)

    # Calculate the rate of change of a, r_u, r_b, and sigma
    da_dt = -F - (x * a) + (z[1] * a_ex) - (z[2] * a)
    dr_u_dt = -F - (x * r_u) + s
    dr_b_dt = F - (x * r_b)
    dsigma_dt = sigma*(x-(sigma*chi_r))

    # Return the rate of change of a, r_u, r_b, and sigma
    return[da_dt, dr_u_dt, dr_b_dt, dsigma_dt]
    pass

#Model equations - check

#Now need to have the solution to the ODEs be a function that will accept the parameter changes for
#k_off, P_in, P_out

def delaytime(z):
    x_intercepts = []
    
    # Set initial conditions
    y0 = [0, r_min + (lambda_0 / kappa_t), 0, sigma_i]

    # Define the time range
    for N in N_values:
        tMax=(N+1)+4+(30/lambda_0)
        time=np.linspace(0,tMax,int(10*tMax))

    growth_rates_list = []

    # Solve the ODEs for each value of N
    for N in N_values:
        solution = solve_ivp(lambda t, y: model_equations(t, y, N, z), [time[0], time[-1]], y0, method='LSODA', t_eval=time)

        gr1_values = np.zeros_like(time)
        gr2_values = np.zeros_like(time)

        # Calculate gr1 and gr2 for all time points
        for i, t in enumerate(time):
            if 0 < t < N + 1:
                gr1_values[i] = ((solution.y[1, i] - r_min) * kappa_t)
            else:
                gr2_values[i] = ((r_max - solution.y[1, i] - solution.y[2, i]) * kappa_n)

        # Calculate the total growth rate by element-wise addition
        growth_rate = gr1_values + gr2_values

        # Append growth rate to the list
        growth_rates_list.append(growth_rate)

    # Calculate OD and x-intercept data
    for N, growth_rate in zip(N_values, growth_rates_list):
        OD = np.cumsum(growth_rate) * (time[1] - time[0])

        start_x = time[-1]
        start_index = np.where(time == start_x)[0][0]
        start_y = OD[start_index]

        if abs(1 - (OD[start_index] - OD[start_index - 1]) / (time[1] - time[0]) / lambda_0) > 0.01:
            print(N, lambda_0)

        desired_slope = lambda_0

        x_intercept = start_x - (start_y / desired_slope)
        x_intercepts.append(x_intercept)
        #print(f"x_intercept for N={N}: {x_intercept}")

    return x_intercepts

def objective_function(parameters):
    k_off, P_in, P_out = parameters
    
    # Initialize the smallest MSE and corresponding parameter values
    smallest_mse = float('inf')
    best_params = None
    best_delay_times = None
    
    # Number of random parameter sets to generate
    num_samples = 10
    
    # Iterate over random parameter sets
    for _ in range(num_samples):
        # Randomly sample k_off, P_in, and P_out within the specified bounds
        k_off_sampled = np.random.uniform(*param_bounds[0])
        P_in_sampled = np.random.uniform(*param_bounds[1])
        P_out_sampled = np.random.uniform(*param_bounds[2])
        
        current_z = [k_off_sampled, P_in_sampled, P_out_sampled]
        
        # Calculate x_intercepts for the current parameter set
        x_intercepts = delaytime(current_z)
        
        # Compute lambda_0_star and IC_50_star based on the current parameter set
        lambda_0_star = 2 * math.sqrt(current_z[2] * kappa_t * (current_z[0] / k_on))
        IC_50_star = (delta_r * lambda_0_star) / (2 * current_z[1])
        
        # Initialize a list to store delay times for this parameter set
        delay_times = []  
        
        # Calculate delay times for each N value
        for i, N in enumerate(N_values):
            calculated_delay_time = x_intercepts[i]
            exp_delay_time = exp_delay_times[N]
            expected_P_out = (IC_50*lambda_0*k_on*current_z[1])/(lambda_max*current_z[0])

            # Store delay time for N if condition is met
            if abs(expected_P_out - current_z[2]) < 0.00000001:
                delay_times.append((N, calculated_delay_time)) 
        
        # If delay times are found for the current parameter set
        if delay_times:  
            # Compute MSE for the parameter set
            mse = 0
            num_delay_times = len(delay_times)
            for N, calculated_delay_time in delay_times:
                exp_delay_time = exp_delay_times[N]
                mse += (exp_delay_time - calculated_delay_time) ** 2
            if num_delay_times > 0:
                mse /= num_delay_times  # Divide by the number of delay times to get the average MSE
                
                # Check if this MSE is the smallest encountered so far
                if mse < smallest_mse:
                    smallest_mse = mse
                    best_params = (k_off_sampled, P_in_sampled, P_out_sampled)
                    best_delay_times = delay_times
    
    # Return the smallest MSE and corresponding parameter values
    return smallest_mse, best_params, best_delay_times

# Bounds for the parameters
param_bounds = [(513, 20520), (1360, 2040), (80, 120)]

Y = objective_function(param_bounds)
print("Smallest MSE:", Y[0])
print("Optimal Parameters (k_off, P_in, P_out):", Y[1])

for N, delay_time in Y[2]:
    print(f"N={N}: Theoretical Delay Time = {delay_time}")

for N, exp_delay_time in exp_delay_times.items():
    print(f"N={N}: Experimental Delay Time = {exp_delay_time}")
    
end = timer()
print((end - start)/60, "minutes.")

I’d rather something along the lines of "if current_z[2] == (IC_50 x lambda_0 x k_on x current_z[1])/(lambda_max x current_z[0]), then delay_times.append((N_ calculated_delay_time)).

The following:

Can be simplified to:

from scipy.optimize import LinearConstraint, Bounds, minimize

Here, assuming that you do not need scipy.optimize anywhere else.

Please try to zoom in on a particular part of the code that you are having an issue with. What you are requesting is rather open ended (what in particular is not working, etc).

When using constants in your code, it is generally good practice to have them in all capitals:

to:

R_MIN = 5.8
R_MAX = 52.3
P_IN= 1485.10542448074
P_OUT = 100.0
k_ON = 1026.0
k_OFF = 2000.0
K_D = K_OFF/K_ON
.
.
.

I don’t use any of the optimization, so I’ve removed them.

It would be this portion of the code, namely the condition after “#Store delay time for N if condition is met” to be replaced with something more strict like “if current_z[2] == (IC_50 x lambda_0 x k_on x current_z[1])/(lambda_max x current_z[0]), then delay_times.append((N_ calculated_delay_time))”.

def objective_function(parameters):
    k_off, P_in, P_out = parameters
    
    # Initialize the smallest MSE and corresponding parameter values
    smallest_mse = float('inf')
    best_params = None
    best_delay_times = None
    
    # Number of random parameter sets to generate
    num_samples = 10
    
    # Iterate over random parameter sets
    for _ in range(num_samples):
        # Randomly sample k_off, P_in, and P_out within the specified bounds
        k_off_sampled = np.random.uniform(*param_bounds[0])
        P_in_sampled = np.random.uniform(*param_bounds[1])
        P_out_sampled = np.random.uniform(*param_bounds[2])
        
        current_z = [k_off_sampled, P_in_sampled, P_out_sampled]
        
        # Calculate x_intercepts for the current parameter set
        x_intercepts = delaytime(current_z)
        
        # Compute lambda_0_star and IC_50_star based on the current parameter set
        lambda_0_star = 2 * math.sqrt(current_z[2] * kappa_t * (current_z[0] / k_on))
        IC_50_star = (delta_r * lambda_0_star) / (2 * current_z[1])
        
        # Initialize a list to store delay times for this parameter set
        delay_times = []  
        
        # Calculate delay times for each N value
        for i, N in enumerate(N_values):
            calculated_delay_time = x_intercepts[i]
            exp_delay_time = exp_delay_times[N]
            expected_P_out = (IC_50*lambda_0*k_on*current_z[1])/(lambda_max*current_z[0])

            # Store delay time for N if condition is met
            if abs(expected_P_out - current_z[2]) < 0.00000001:
                delay_times.append((N, calculated_delay_time)) 
        
        # If delay times are found for the current parameter set
        if delay_times:  
            # Compute MSE for the parameter set
            mse = 0
            num_delay_times = len(delay_times)
            for N, calculated_delay_time in delay_times:
                exp_delay_time = exp_delay_times[N]
                mse += (exp_delay_time - calculated_delay_time) ** 2
            if num_delay_times > 0:
                mse /= num_delay_times  # Divide by the number of delay times to get the average MSE
                
                # Check if this MSE is the smallest encountered so far
                if mse < smallest_mse:
                    smallest_mse = mse
                    best_params = (k_off_sampled, P_in_sampled, P_out_sampled)
                    best_delay_times = delay_times
    
    # Return the smallest MSE and corresponding parameter values
    return smallest_mse, best_params, best_delay_times

The name of the function is: objective_function. This means nothing to someone who is reviewing your code. Give it a name that has meaning - relevant to what it is actually doing.

You’re converting a string to a float? Why?

Try to give meaning to what something is doing. An underscore is vague.

Change it to:

for sample in range(num_samples):`

In the following line(s), what is the purpose of the askerisk *? Is this required syntax by the numpy object?

`k_off_sampled = np.random.uniform(*param_bounds[0])`

Your question is still too vague. What is it not doing that you expect it to?

Example: The output is 100. But based on my calculations, I should be getting 345. etc.

  • By the way, the following comparison, doesn’t make sense:

The reason is that expected_P_out has units of watts and current_z[2] has units of amps.

Thus, they can’t be subtracted from each other. The set up of the code is wrong.

Thank you for the reply.

  1. Changed to a more understandable function name;
  2. This is just a placeholder value before smallest_mse is calculated later;
  3. Changed from vague underscore;
  4. If the form *identifier is present, it is initialized to a tuple receiving any excess positional parameters, defaulting to the empty tuple.

Now with the formalities out of the way, what I am trying to accomplish is a stronger constraint on a set of parameter values. I am using a Monte Carlo simulation to take the smallest MSE taken from uniform random parameter sets within specified parameter bounds, but each value needs to satisfy a very strict constraint (due to biological relevance). It computes the MSE via the difference between my experimental data and theoretical (Python calculated) data. I want my code to only calculate the delay time if the following condition is met:

if current_z[2] == (IC_50*lambda_0*k_on*current_z[1])/(lambda_max x current_z[0]), then delay_times.append((N_ calculated_delay_time)), i.e. only calculate the delay time with the parameter set z that satisfies that condition. A tolerance-based constraint does not make sense in this setting. I tried earlier an == condition, as well as a != condition, both of which returned no MSE and no parameter sets which should not happen. I attach the code again with more thorough annotation:

def minimize_MSE(parameters):
    k_off, P_in, P_out = parameters
    
    # Initialize the smallest MSE and corresponding parameter values
    smallest_mse = float('inf') #placeholder to initialize smallest MSE
    best_params = None
    best_delay_times = None
    
    # Number of random parameter sets to generate
    num_samples = 10
    
    # Iterate over random parameter sets
    for samples in range(num_samples):
        # Randomly sample k_off, P_in, and P_out within the specified bounds
        k_off_sampled = np.random.uniform(*param_bounds[0])
        P_in_sampled = np.random.uniform(*param_bounds[1])
        P_out_sampled = np.random.uniform(*param_bounds[2])
        
        current_z = [k_off_sampled, P_in_sampled, P_out_sampled]
        
        # Calculate x_intercepts for the current parameter set
        x_intercepts = delaytime(current_z)
        
        # Compute lambda_0_star and IC_50_star based on the current parameter set
        lambda_0_star = 2 * math.sqrt(current_z[2] * kappa_t * (current_z[0] / k_on))
        IC_50_star = (delta_r * lambda_0_star) / (2 * current_z[1])
        
        # Initialize a list to store delay times for this parameter set
        delay_times = []  
        
        # Calculate delay times for each N value
        for i, N in enumerate(N_values):
            calculated_delay_time = x_intercepts[i]
            exp_delay_time = exp_delay_times[N]
            expected_P_out = (IC_50*lambda_0*k_on*current_z[1])/(lambda_max*current_z[0])

            # Store delay time for N if condition is met
            # This condition is incorrect and not what I am looking to accomplish, this needs to be a more strict "equals to" for expected_P_out. The difference between the P_out's should not be considered
            if abs(expected_P_out - current_z[2]) < 0.00000001:
                delay_times.append((N, calculated_delay_time)) 
        
        # If delay times are found for the current parameter set
        if delay_times:  
            # Compute MSE for the parameter set
            mse = 0
            num_delay_times = len(delay_times)
            for N, calculated_delay_time in delay_times:
                exp_delay_time = exp_delay_times[N]
                mse += (exp_delay_time - calculated_delay_time) ** 2
            if num_delay_times > 0:
                mse /= num_delay_times  # Divide by the number of delay times to get the average MSE
                
                # Check if this MSE is the smallest encountered so far
                if mse < smallest_mse:
                    smallest_mse = mse
                    best_params = (k_off_sampled, P_in_sampled, P_out_sampled)
                    best_delay_times = delay_times
    
    # Return the smallest MSE and corresponding parameter values
    return smallest_mse, best_params, best_delay_times

The following conditional statement still needs to be modified to something that makes sense:

You can’t have power - current. They have different units. So it has no meaning.

The following conditional statement, also doesn’t make sense:

if current_z[2] == (IC_50*lambda_0*k_on*current_z[1])/(lambda_max x current_z[0])

This is an analog value (not digital). It makes more sense to have >= or <= instead of ==. The chances that the values are EXACTLY equal are almost zero. What if the value on the left side is 25.012244 and the value on the right-side is 25.011343. They wont be equal even if they’re off by ~0.001 (1 mA).

Sorry, I should clarify - z[2] and expected_P_out have the same units, z[2]=P_out which is an outward cellular permeability rate. If they cannot be exactly equal, do you have a recommendation for the constraint? Essentially the original constraint I had was where IC_50=1.75, calc_IC_50=((kappa_t*delta_r)/lambda_0)*(k_off/k_on)*(P_out/P_in)and was testing if abs(calc_IC_50-IC_50)<0.01, then use the parameters in z to calculate delay time. However, iterating over an O(n^3) took an extremely long time for a large number of sample. Instead, I want to isolate for one of the parameters (in this case, P_out) and reduce the problem to an O(n^2) problem. But I don’t have a condition like abs(X-Y) for this constraint, as I just want P_out=RHS.

One other observation of which I will start with a small example:

for number in range(1, 6):
    print(number, end= ' ')

Notice that the number variable is in the body to be used.

In your iterable loop statement:

for samples in range(num_samples):

where is samples used in your code? As far as I can tell, it is not being used since it is not in the body. Therefore, the body results will always be static - no changes will be generated.

Please don’t include code inside text paragraphs. It makes it a bit more difficult to read. Instead, give them their own lines separately:


IC_50 = 1.75
calc_IC_50 = ((kappa_t*delta_r)/lambda_0) * (k_off/k_on) * (P_out/P_in)

if abs(calc_IC_50 - IC_50) < 0.01

The recommendation that I provided in my previous post is still valid. Don’t use ==.

Sure, I understand why == won’t work. Then why is it when I use negation as opposed to satisfaction (i.e.:

if current_z[2] != (IC_50*lambda_0*k_on*current_z[1])/(lambda_max x current_z[0]):
continue

it still returns no parameter set z nor any calculated MSE? I cannot use <= or >=, as it must be strictly equal to (though I understand based on the arithmetic it won’t be “equal”).

Notice the following:

(lambda_max x current_z[0]) # The 'x' has no meaning since it is not multiplication

You need to use an asterisk *.

Sorry, that’s a typo by me. In my Python console, it’s correct. That still doesn’t answer my question all the same.

From my previous post, samples is not being used anywhere in your code so the outputs will remain the same for every iteration.

Right, this is why I had the _ in the beginning. I don’t care about the variable, I care about the iteration over the bounds of the three parameters.

How else would you test and verify that your code is working? Every system has inputs and outputs. Your code should be able to provide expected outputs provided a set of inputs.

The underscore has nothing to do with this.

My recommendation is for you to test smaller parts of your code in isolation and ensure that they are providing the expected results based on test inputs. Use your own test values. Are you getting the expected outputs? And so on before creating an application.

Not a problem. I’ll do some more debugging, fine-tuning and will figure it out. Thank your for your input!

Just as an aside - I only included the code relevant to the problem at hand here. I have all “fundamentals” like function definitions and variable definitions in my code and written correctly, as well as having imported all relevant packages.