Imposed SSCs with optimal MTC shortening/lengthening velocity

The analysis of this page corresponds to the sections ‘Imposed stretch-shortening cycles — optimal MTC shortening/lengthening velocity’.

Because the analysis is relatively extensive, it is split into separate scripts included on this page:

1. Predict the maximally attainable AMPO for one imposed SSC parameter (i.e. cycle frequency, FTS or MTC length excursion), while leaving CE length, MTC length and stimulation over time otherwise unconstrained.
2. Check convergence of the optimal control solutions by comparing repeated optimisation runs for each imposed parameter value and specimen.
3. Compare the maximally attainable AMPO of the unconstrained optimal-control SSCs with the constant MTC shortening/lengthening velocity SSCs.

Custom functions used:

• run_ssc_ocp(muspar, cf=None, fts=None, mle=None, N=200, d=3, do_plot=False, do_print=True)

  Solve the optimal control problem for one specimen and one imposed SSC parameter, returning MTC length, stimulation, SEE force, active state and CE length over time.

• ca_func.create_funcs(muspar)

  Create CasADi functions for the Hill-type MTC dynamics, mechanical power objective and helper variables used in the optimal control problem.

• ca_func.get_sim_guess(f_dyn, cf0, fts0, mle0, N)

  Generate an initial guess for the optimal control problem from a constant-velocity SSC simulation.

• ca_func.ode(x, u, muspar, b=1e3)

  Evaluate excitation dynamics, CE force, CE velocity, MTC length, active state and SEE force for the Hill-type MTC model.

• helpers.get_ampo(filepaths)

  Load simulation files and compute AMPO values for comparison between constant-velocity and optimal-control SSCs.

Compute AMPO for imposed SSC parameters

"""
This script predicts the maximally attainable AMPO when one SSC parameter is
imposed and MTC length, CE length and stimulation over time are otherwise free
to vary. The optimisation is formulated as a direct-collocation optimal control
problem using the Hill-type MTC model.

Specifically, the following steps were taken:

-   Define the optimal control problem for a periodic SSC.
-   Impose either cycle frequency, FTS, or MTC length excursion.
-   Optimise CE velocity and stimulation over time to maximise AMPO.
-   Save repeated optimisation results for later convergence checks.

Custom functions used:

-   `run_ssc_ocp(muspar, cf=None, fts=None, mle=None, N=200, d=3, do_plot=False, do_print=True)`
    :   Solve the optimal control problem for one specimen and one imposed SSC
        parameter.
-   `ca_func.create_funcs(muspar)`
    :   Create CasADi functions for the Hill-type MTC dynamics, mechanical
        power objective and helper variables.
-   `ca_func.get_sim_guess(f_dyn, cf0, fts0, mle0, N)`
    :   Generate an initial guess from a constant-velocity SSC simulation.
-   `ca_func.ode(x, u, muspar, b=1e3)`
    :   Evaluate the Hill-type MTC model states and outputs.
"""

#%% Load CasADi & set directories
import os, sys, pickle
import casadi as ca
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import integrate
from pathlib import Path

# Set directories
cwd = Path.cwd()
baseDir = cwd.parent.parent
dataDir = baseDir / 'data'
funcDir = baseDir / 'analysis' / 'functions'
sys.path.append(str(funcDir))

import ca_func

plt.close('all')

#%% Create optimisation function
def run_ssc_ocp(muspar, cf=None, fts=None, mle=None, N=200, d=3, do_plot=False, do_print=True):
    """
    Solve one optimal-control SSC problem.

    One of the SSC parameters can be imposed by passing `cf`, `fts` or `mle`.
    Parameters that are passed as `None` become optimisation variables. CE
    velocity and stimulation are optimised over a periodic cycle, subject to
    excitation dynamics, force-length and force-velocity properties, and basic
    path constraints. Successful solutions are returned as a table containing
    time, MTC length, stimulation, SEE force, active state and relative CE
    length.
    """
    f_dyn, f_cost, f_var = ca_func.create_funcs(muspar)
    tau = ca.collocation_points(d, 'legendre')
    [C, D, B] = ca.collocation_coeff(tau)
    
    N2 = N//2
    N1 = N-N2
    
    # Set up optimization problem
    opti = ca.Opti()
    
    # Set cf, fts, mle (and their initial guess)
    if fts is None:
        fts = opti.variable()
        opti.subject_to(opti.bounded(0.01,fts,0.99))
        fts0 = np.random.uniform(low=0.5, high=0.8)
        opti.set_initial(fts, fts0)
    else:
        fts0 = fts
        
    if cf is None:
        cf = opti.variable()
        opti.subject_to(opti.bounded(0.2,cf,12))
        cf0 = np.random.uniform(low=1, high=fts0/0.05) # minimal 50ms Tshort
        opti.set_initial(cf, cf0)
    else:
        cf0 = cf 
    
    if mle is None:
        mle0 = np.random.uniform(low=0.05, high=0.2)
    else:
        mle0 = mle
    
    # Decision variables
    x_k         = opti.variable(2,N+1)
    u_k         = opti.variable(2,N)
    p_k         = opti.variable(1)
    
    # Simulate the system first (to obtain initial guess)
    x0,u0,t0 = ca_func.get_sim_guess(f_dyn,cf0,fts0,mle0/muspar['lce_opt'],N)
    
    # Initial guesses 
    opti.set_initial(x_k, x0)
    opti.set_initial(u_k, u0)
    p_k = 1e3
        
    # Dynamics constraints
    J = 0
    t_k = [0] # [s] time-axis
    lmtc_k = [] 
    for k in range(N):
        if k<N1:
            dt = fts/cf/N1
        else:
            dt = (1-fts)/cf/N2
                
        # Collect states at k and k+1
        Xk = x_k[:, k]
        Xk_next = x_k[:, k + 1]
    
        # Collect controls at k
        Uk = u_k[:, k]
    
        # Collect helper variables at each collocation point
        Xc = opti.variable(2,d)
        # opti.set_initial(Xc, np.repeat(np.array([[1.0, 0.5]]),d,axis=0).T)
        opti.set_initial(Xc, np.tile(x0[:,k], (d, 1)).T)
        
        # Compute cost function and mechanical work
        Jp = -f_cost(Xc, Uk, p_k)
        J = J + ca.mtimes(Jp,B) * dt # cost
    
        # Get interpolating points of collocation polynomial
        Z = ca.horzcat(Xk, Xc)
        # Get slope of interpolating polynomial (normalized)
        Pidot = ca.mtimes(Z, C) / dt
        # State at end of collocation interval
        Xk_end = ca.mtimes(Z, D)
        opti.subject_to(Xk_end == Xk_next)
    
        # Explicit dynamics
        ode = f_dyn(Xc, Uk, p_k)
        opti.subject_to(Pidot == ode)
        
        # Constrain on CE force
        fce_c = f_var(Xc,Uk, p_k)[0]
        opti.subject_to(fce_c >= 0)
        
        # Time-axis
        t_k.append(t_k[-1]+dt)
        
        # MTC length
        lmtc_k.append(f_var(Xk,Uk, p_k)[1])
        
    t_k =  ca.vertcat(*t_k) 
    lmtc_k =  ca.vertcat(*lmtc_k)
    
    # Periodic constraints
    opti.subject_to(x_k[:,0]    == x_k[:,-1])
    
    # Path constraints 
    gamma_k     = x_k[0,:]
    lcerel_k    = x_k[1,:];     lce_k = lcerel_k*muspar['lce_opt']
    
    stim_k      = u_k[0,:]
    vcerel_k    = u_k[1,:];     vcerel_k1 = vcerel_k[:,:N1+1];  vcerel_k2 = vcerel_k[:,N1:]
    
    if mle is not None:
        opti.subject_to(opti.bounded(0.8*mle,lce_k[0]-lce_k[N1],mle*1.2))
        opti.subject_to(lmtc_k[0]-lmtc_k[N1] == mle)
        opti.subject_to(vcerel_k1[0] >= 2*vcerel_k1[1:])
        # opti.subject_to(lce_k[0]-lce_k[N1] == mle)
    opti.subject_to(opti.bounded(muspar['gamma_0'],gamma_k,1))
    opti.subject_to(opti.bounded(1-muspar['w'],lcerel_k,1+muspar['w']))
    opti.subject_to(opti.bounded(0,stim_k,1))
    opti.subject_to(vcerel_k1 <= 0)
    opti.subject_to(vcerel_k2 >= 0)
    
    # Optimize    
    opti.minimize(J*cf*1e1) # Objective
    
    # Set solver and solve the optimization problem
    try:
        opti.solver('ipopt',{'ipopt.max_iter': 500, 
                              'ipopt.mu_strategy': 'adaptive', # monotone = standard, other option: adaptive
                              'ipopt.hessian_approximation': 'limited-memory', # exact vs. limited-memory
                              'ipopt.tol': 5e-3
                              })
        sol = opti.solve()
    except:
        sol = opti.debug
    if mle is None:
        hessian_approx = 'limited-memory' # apparantly works  better for imposed MLE..
    else:
        hessian_approx = 'exact'
    try:        
        opti.solver('ipopt', {
            'ipopt.max_iter': 1000,
            'ipopt.mu_strategy': 'monotone',       # Good for warm start
            'ipopt.hessian_approximation': hessian_approx,
            'ipopt.warm_start_init_point': 'yes',
            'ipopt.warm_start_bound_push': 1e-4,
            'ipopt.warm_start_mult_bound_push': 1e-4,   # Also warm-start dual feasibility
            'ipopt.bound_push': 1e-4,
            'ipopt.bound_frac': 1e-4,
            'ipopt.accept_every_trial_step': 'no'  # Make Ipopt more conservative accepting steps
        })
        opti.set_initial(opti.x, sol.value(opti.x)) # Set primal warm-start
        opti.set_initial(opti.lam_g, sol.value(opti.lam_g)) #  Set constraint duals
        sol = opti.solve();
    except:
        sol = opti.debug

    # Extract and display the solution
    if sol.stats()['success'] is True: # only if succesfully optimised!
        if not isinstance(cf, float): cf = sol.value(cf)
        if not isinstance(fts, float): fts = sol.value(fts)
        
        t   = sol.value(t_k)
        x   = sol.value(x_k)
        u   = sol.value(u_k)
        u = np.hstack((u,u[:,0:1]))
        _,fce,vce,lmtc,q,fsee = ca_func.ode(x,u,muspar)[0:6]
        fce = np.squeeze(fce)
        lmtc = np.squeeze(lmtc)
        fsee = np.squeeze(fsee)
        
        gamma, lcerel = x
        stim, vcerel = u
        
        data = np.vstack((t,lmtc,stim,fsee,gamma,lcerel)).T
        df = pd.DataFrame(data)
        
        if do_plot == True:
            plt.figure()
            plt.subplot(311)
            plt.plot(t, x.T)
            
            plt.subplot(312)
            plt.plot(t, u.T)
        
        if do_print == True:
            # Calculate AMPO
            Wmech = -integrate.cumulative_trapezoid(fce,lcerel*muspar['lce_opt'])
            AMPO = Wmech[-1]*cf # [mW]
            print("AMPO = %1.2f mW" % (AMPO*1e3))
            
            # Check SSC parameters
            cf = 1/t[-1]
            iMin = np.argmin(lcerel)
            fts = t[iMin]*cf
            mle = lmtc[0]-lmtc[N1]
            print("CF = %1.2f Hz" % cf)
            print("FTS = %1.2f" % fts)
            print("MLE = %1.2f mm" % (mle*1e3))
    else:
        df = None
    return df

#%% Load muscle parameters
mus = 'GMe3'
parFile = os.path.join(dataDir, mus, mus + '_IM.pkl')
muspar = pickle.load(open(parFile, 'rb'))[0]
# Now we change lsee0 to 3mm
muspar_lsee0_new = 3e-3
muspar['ksee'] = muspar['ksee']*(muspar['lsee0']/muspar_lsee0_new)**2
muspar['lsee0'] = muspar_lsee0_new
dataDirSim = os.path.join(dataDir,mus,'simsOC','it','')

#%% Perform OC optimisation
for cf in np.arange(0.5,8.1,0.5):
    for i in range(1,6):
        df = run_ssc_ocp(muspar,cf,None,None)
        fileName = mus+f'_cf{cf:0.1f}Hz_ftsOpt_mleOpt_it{i:02d}'
        if df is not None:
            df.to_csv(dataDirSim+fileName+'.csv',index=False,header=['Time [s]','Lmtc [m]','STIM [ ]', 'Fsee [N]', 'Gamma [ ]','Lcerel [ ]'])

for fts in np.arange(0.05,0.96,0.05):
    for i in range(1,6):
        df = run_ssc_ocp(muspar,None,fts,None)
        fileName = mus+f'_cfOpt_fts{fts:0.2f}_mleOpt_it{i:02d}'
        if df is not None:
            df.to_csv(dataDirSim+fileName+'.csv',index=False,header=['Time [s]','Lmtc [m]','STIM [ ]', 'Fsee [N]', 'Gamma [ ]','Lcerel [ ]'])

for mle in np.arange(1,11.1,1):
    for i in range(1,6):
        df = run_ssc_ocp(muspar,None,None,mle/1e3)
        fileName = mus+f'_cfOpt_ftsOpt_mle{mle:04.1f}mm_it{i:02d}'
        if df is not None:
            df.to_csv(dataDirSim+fileName+'.csv',index=False,header=['Time [s]','Lmtc [m]','STIM [ ]', 'Fsee [N]', 'Gamma [ ]','Lcerel [ ]'])

Check convergence of Optimal Control solutions

"""
This script checks whether the optimal-control predictions are robust to the
random initial guesses used during optimisation.

For every specimen and imposed SSC parameter value, five optimisation attempts
were run by `ssc-oc-run-predictions.py`. This script reloads those attempts,
computes AMPO for each successful solution and summarises the convergence
quality using:

-   the number of failed or missing optimisation runs;
-   the percentage spread between repeated successful solutions;
-   the overall success rate across cycle frequency, FTS and MTC length
    excursion sweeps.
"""

#%% Load packages & set directories
import os, sys, pickle
import numpy as np
import pandas as pd
from scipy import integrate
from pathlib import Path

# Set directories
cwd = Path.cwd()
baseDir = cwd.parent.parent
dataDir = baseDir / 'data'
funcDir = baseDir / 'analysis' / 'functions'
sys.path.append(str(funcDir))

#%% Check convergence for imposed cycle frequency
iters = range(1, 6)
muscles = ['GMe1', 'GMe2', 'GMe3']

iTotal = 0
iFail = 0
pct_diffs = []

#%% Imposed CF
# Define CFs, iterations, and mus
cf_vals = np.arange(1, 6.1, 0.5)
cols = [f"{cf:.1f}" for cf in cf_vals]

# Initialize DataFrame (3D by adding 'mus' as another index)
df_cf = pd.DataFrame(index=iters, columns=pd.MultiIndex.from_product([muscles, cols], names=['mus', 'cf']), dtype=float)

# Loop over mus, cf, and iterations
for mus in muscles:
    parFile = os.path.join(dataDir, mus, mus + '_IM.pkl')
    muspar = pickle.load(open(parFile, 'rb'))[0]
    dataDirSim = os.path.join(dataDir, mus, 'simsOC', 'it', '')
    for cf_str in cols:
        cf = float(cf_str)
        AMPOset = []
        for i in iters:
            filepath = os.path.join(dataDirSim, f"{mus}_cf{cf_str}Hz_ftsOpt_mleOpt_it{i:02d}.csv")
            try:
                df = pd.read_csv(filepath)
                _, _, _, fsee, _, lcerel = df.to_numpy().T
                fce = fsee  # TEMP
                Wmech = -integrate.cumulative_trapezoid(fce, lcerel * muspar['lce_opt'])
                df_cf.loc[i, (mus, cf_str)] = Wmech[-1] * 1e3 * cf
            except:
                continue  # leave as NaN
        # Compute stats if possible
        if len(AMPOset) >= 2:
            AMPOset = -np.sort(-np.array(AMPOset))
            max_val = max(AMPOset)
            pct_diff = [(max_val - val) / max_val * 100 for val in AMPOset][1:]
            pct_diffs.append(pct_diff)   
        
# Compute stats (here stats are computed across the 2D slice per mus)
stats = pd.DataFrame({
    'percent_diff': ((df_cf.max(axis=0) - df_cf.min(axis=0)) / df_cf.max(axis=0) * 100).round(2),
    'nan_count': df_cf.isna().sum(axis=0)
}).T

# Append stats to each level of 'mus'
df_cf = pd.concat([df_cf, stats])

# Calculate total nan_count and highest percent_diff
total_nan_count = stats.loc['nan_count'].sum()
highest_percent_diff = stats.loc['percent_diff'].max()
n_opt = len(iters)*len(cf_vals)*len(muscles)

# Print total nan_count and highest percent_diff
print("\nTotal nan_count across all mus and cf combinations:", total_nan_count)
print(f"Succes rate: {(1-total_nan_count/n_opt)*100:0.2f}%")
print("Highest percent_diff across all mus and cf combinations:", highest_percent_diff)

iTotal += n_opt
iFail += total_nan_count

#%% Check convergence for imposed FTS
# Define FTSs, iterations, and mus
fts_vals = np.arange(0.05, 0.96, 0.05)
cols = [f"{fts:.2f}" for fts in fts_vals]

# Initialize DataFrame (3D by adding 'mus' as another index)
df_fts = pd.DataFrame(index=iters, columns=pd.MultiIndex.from_product([muscles, cols], names=['mus', 'cf']), dtype=float)

# Loop over mus, cf, and iterations
for mus in muscles:
    parFile = os.path.join(dataDir, mus, mus + '_IM.pkl')
    muspar = pickle.load(open(parFile, 'rb'))[0]
    dataDirSim = os.path.join(dataDir, mus, 'simsOC', 'it', '')
    for fts_str in cols:
        fts = float(fts_str)
        AMPOset = []
        for i in iters:
            filepath = os.path.join(dataDirSim, f"{mus}_cfOpt_fts{fts:0.2f}_mleOpt_it{i:02d}.csv")
            try:
                df = pd.read_csv(filepath)
                time, _, _, fsee, _, lcerel = df.to_numpy().T
                fce = fsee  # TEMP
                Wmech = -integrate.cumulative_trapezoid(fce, lcerel * muspar['lce_opt'])
                df_fts.loc[i, (mus, fts_str)] = Wmech[-1] * 1e3 / time[-1]
            except:
                continue  # leave as NaN
        # Compute stats if possible
        if len(AMPOset) >= 2:
            AMPOset = -np.sort(-np.array(AMPOset))
            max_val = max(AMPOset)
            pct_diff = [(max_val - val) / max_val * 100 for val in AMPOset][1:]
            pct_diffs.append(pct_diff)        

# Compute stats (here stats are computed across the 2D slice per mus)
stats = pd.DataFrame({
    'percent_diff': ((df_fts.max(axis=0) - df_fts.min(axis=0)) / df_fts.max(axis=0) * 100).round(2),
    'nan_count': df_fts.isna().sum(axis=0)
}).T

# Append stats to each level of 'mus'
df_fts = pd.concat([df_fts, stats])

# Calculate total nan_count and highest percent_diff
total_nan_count = stats.loc['nan_count'].sum()
highest_percent_diff = stats.loc['percent_diff'].max()
n_opt = len(iters)*len(fts_vals)*len(muscles)

# Print total nan_count and highest percent_diff
print("\nTotal nan_count across all mus and fts combinations:", total_nan_count)
print(f"Succes rate: {(1-total_nan_count/n_opt)*100:0.2f}%")
print("Highest percent_diff across all mus and fts combinations:", highest_percent_diff)

iTotal += n_opt
iFail += total_nan_count

#%% Check convergence for imposed MTC length excursion
mle_vals = np.arange(1, 11.1, 1)
cols = [f'{mle:04.1f}' for mle in mle_vals]

# Initialize DataFrame (3D by adding 'mus' as another index)
df_mle = pd.DataFrame(index=iters, columns=pd.MultiIndex.from_product([muscles, cols], names=['mus', 'mle']), dtype=float)

# Main loop
for mus in muscles:
    parFile = os.path.join(dataDir, mus, mus + '_IM.pkl')
    muspar = pickle.load(open(parFile, 'rb'))[0]
    dataDirSim = os.path.join(dataDir, mus, 'simsOC', 'it', '')
    
    for mle_str in cols:
        mle = float(mle_str)
        AMPOset = []
        for i in iters:
            filepath = os.path.join(dataDirSim, f"{mus}_cfOpt_ftsOpt_mle{mle:04.1f}mm_it{i:02d}.csv")
            try:
                df = pd.read_csv(filepath)
                time, _, _, fsee, _, lcerel = df.to_numpy().T
                fce = fsee  # TEMP
                AMPO = -integrate.trapezoid(fce, lcerel * muspar['lce_opt']) * 1e3 / time[-1]
                AMPOset.append(AMPO)
                df_mle.loc[i, (mus, mle_str)] = AMPO
            except:
                continue  # leave as NaN
        # Compute stats if possible
        if len(AMPOset) >= 2:
            AMPOset = -np.sort(-np.array(AMPOset))
            max_val = max(AMPOset)
            pct_diff = [(max_val - val) / max_val * 100 for val in AMPOset][1:]
            pct_diffs.append(pct_diff)        

# Compute stats (here stats are computed across the 2D slice per mus)
stats = pd.DataFrame({
    'percent_diff': ((df_mle.max(axis=0) - df_mle.min(axis=0)) / df_mle.max(axis=0) * 100).round(2),
    'nan_count': df_mle.isna().sum(axis=0),
}).T

# Append stats to each level of 'mus'
df_mle = pd.concat([df_mle, stats])

# Calculate total nan_count and highest percent_diff
total_nan_count = stats.loc['nan_count'].sum()
highest_percent_diff = stats.loc['percent_diff'].max()
n_opt = len(iters)*len(mle_vals)*len(muscles)

# Print total nan_count and highest percent_diff
print("\nTotal nan_count across all mus and fts combinations:", total_nan_count)
print(f"Succes rate: {(1-total_nan_count/n_opt)*100:0.2f}%")
print("Highest percent_diff across all mus and fts combinations:", highest_percent_diff)

iTotal += n_opt
iFail += total_nan_count

#%% Summarise convergence across all optimisation runs
pct_diffs = [item for sublist in pct_diffs for item in sublist]
print(np.mean(pct_diffs))
print(np.std(pct_diffs))
print(iFail/iTotal*100) # thus 98.2% converged

# iTotal = 615, iFail = 11. Thus 604 succesfull..

Compare optimal-control and constant-velocity SSCs

"""
This script quantifies how much AMPO increases when MTC length over time is
optimised freely, compared with SSCs that impose constant MTC shortening and
lengthening velocities.

For each specimen, AMPO is loaded for matching constant-velocity and
optimal-control simulations. Ratios are computed separately for sweeps with
imposed cycle frequency, imposed FTS and imposed MTC length excursion. The
printed values report the mean percentage increase of the optimal-control
solutions over the constant-velocity solutions.
"""

#%% Load packages & set directories
import os, sys
import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
from pathlib import Path

# Set directories
cwd = Path.cwd()
baseDir = cwd.parent.parent
dataDir = baseDir / 'data'
funcDir = baseDir / 'analysis' / 'functions'
sys.path.append(str(funcDir))

import helpers

plt.close('all')
        
#%% Compare simulations with imposed cycle frequency
cfSet = np.arange(1.0,6.1,0.5) # Hz

d = []
for mus in ['GMe1', 'GMe2', 'GMe3']:
    dataFolder = os.path.join(dataDir,mus,'simsCV','')
    filepaths_cf    = [os.path.join(dataFolder, f'{mus}_cf{cf:0.1f}Hz_ftsOpt_mleOpt.csv') for cf in cfSet]
    AMPOssc = helpers.get_ampo(filepaths_cf)
    
    dataFolder = os.path.join(dataDir,mus,'simsOC','')
    filepaths_cf    = [os.path.join(dataFolder, f'{mus}_cf{cf:0.1f}Hz_ftsOpt_mleOpt.csv') for cf in cfSet]
    AMPOoc = helpers.get_ampo(filepaths_cf)

    d.append(AMPOoc/AMPOssc)
avg,std = (np.nanmean(d)-1)*100, np.nanstd(d)*100
print(f'For imposed cycle frequency OC is {avg:0.1f} +- {std:0.1f}% higher')

#%% Compare simulations with imposed FTS
ftsSet = np.arange(0.25,0.96,0.05)

d = []
for mus in ['GMe1', 'GMe2', 'GMe3']:
    dataFolder = os.path.join(dataDir,mus,'simsCV','')
    filepaths_fts   = [os.path.join(dataFolder, f'{mus}_cfOpt_fts{fts:0.2f}_mleOpt.csv') for fts in ftsSet]
    AMPOssc = helpers.get_ampo(filepaths_fts)
    
    dataFolder = os.path.join(dataDir,mus,'simsOC','')
    filepaths_fts   = [os.path.join(dataFolder, f'{mus}_cfOpt_fts{fts:0.2f}_mleOpt.csv') for fts in ftsSet]
    AMPOoc = helpers.get_ampo(filepaths_fts)

    d.append(AMPOoc/AMPOssc)
avg,std = (np.nanmean(d)-1)*100, np.nanstd(d)*100
print(f'For imposed FTS OC is {avg:0.1f} +- {std:0.1f}% higher')

#%% Compare simulations with imposed MTC length excursion
mleSet = np.arange(1,12.1,1) # mm

d = []
for mus in ['GMe1', 'GMe2', 'GMe3']:
    dataFolder = os.path.join(dataDir,mus,'simsCV','')
    filepaths_mle   = [os.path.join(dataFolder, f'{mus}_cfOpt_ftsOpt_mle{mle:0.1f}mm.csv') for mle in mleSet]
    AMPOssc = helpers.get_ampo(filepaths_mle)
    
    dataFolder = os.path.join(dataDir,mus,'simsOC','')
    filepaths_mle   = [os.path.join(dataFolder, f'{mus}_cfOpt_ftsOpt_mle{mle:04.1f}mm.csv') for mle in mleSet]
    AMPOoc = helpers.get_ampo(filepaths_mle)

    d.append(AMPOoc/AMPOssc)
avg,std = (np.nanmean(d)-1)*100, np.nanstd(d)*100
print(f'For imposed MTC length excursion OC is {avg:0.1f} +- {std:0.1f}% higher')