diff --git a/cocofest/examples/dynamics/reaching_task/reaching_task_pulse_duration_optimization_testing.py b/cocofest/examples/dynamics/reaching_task/reaching_task_pulse_duration_optimization_testing.py deleted file mode 100644 index 4eacc32..0000000 --- a/cocofest/examples/dynamics/reaching_task/reaching_task_pulse_duration_optimization_testing.py +++ /dev/null @@ -1,143 +0,0 @@ -""" -This example will do a pulse duration optimization to either minimize overall muscle force or muscle fatigue -for a reaching task. Those ocp were build to move from starting position (arm: 0°, elbow: 5°) to a target position -defined in the bioMod file. At the end of the simulation 2 files will be created, one for each optimization. -The files will contain the time, states, controls and parameters of the ocp. -""" - -import pickle - -from bioptim import ( - Axis, - ConstraintFcn, - ConstraintList, - Node, - Solver, -) - -from cocofest import DingModelPulseDurationFrequencyWithFatigue, OcpFesMsk, SolutionToPickle - -# Scaling alpha_a and a_scale parameters for each muscle proportionally to the muscle PCSA and fiber type 2 proportion -# Fiber type proportion from [1] -biceps_fiber_type_2_proportion = 0.607 -triceps_fiber_type_2_proportion = 0.465 -brachioradialis_fiber_type_2_proportion = 0.457 -alpha_a_proportion_list = [ - biceps_fiber_type_2_proportion, - biceps_fiber_type_2_proportion, - triceps_fiber_type_2_proportion, - triceps_fiber_type_2_proportion, - triceps_fiber_type_2_proportion, - brachioradialis_fiber_type_2_proportion, -] - -# PCSA (cm²) from [2] -triceps_pcsa = 28.3 -biceps_pcsa = 12.7 -brachioradialis_pcsa = 11.6 -triceps_a_scale_proportion = 1 -biceps_a_scale_proportion = biceps_pcsa / triceps_pcsa -brachioradialis_a_scale_proportion = brachioradialis_pcsa / triceps_pcsa -a_scale_proportion_list = [ - biceps_a_scale_proportion, - biceps_a_scale_proportion, - triceps_a_scale_proportion, - triceps_a_scale_proportion, - triceps_a_scale_proportion, - brachioradialis_a_scale_proportion, -] - -# Build the functional electrical stimulation models according -# to number and name of muscle in the musculoskeletal model used -fes_muscle_models = [ - DingModelPulseDurationFrequencyWithFatigue(muscle_name="BIClong"), - DingModelPulseDurationFrequencyWithFatigue(muscle_name="BICshort"), - DingModelPulseDurationFrequencyWithFatigue(muscle_name="TRIlong"), - DingModelPulseDurationFrequencyWithFatigue(muscle_name="TRIlat"), - DingModelPulseDurationFrequencyWithFatigue(muscle_name="TRImed"), - DingModelPulseDurationFrequencyWithFatigue(muscle_name="BRA"), -] - -# Applying the scaling -for i in range(len(fes_muscle_models)): - fes_muscle_models[i].alpha_a = fes_muscle_models[i].alpha_a * alpha_a_proportion_list[i] - fes_muscle_models[i].a_scale = fes_muscle_models[i].a_scale * a_scale_proportion_list[i] - -minimum_pulse_duration = DingModelPulseDurationFrequencyWithFatigue().pd0 -# pickle_file_list = ["minimize_muscle_fatigue", "minimize_muscle_force"] -pickle_file_list = ["minimize_muscle_force"] -n_stim = 60 -n_shooting = 25 -# Step time of 1ms -> 1sec / (40Hz * 25) = 0.001s - -constraint = ConstraintList() -constraint.add( - ConstraintFcn.SUPERIMPOSE_MARKERS, - first_marker="COM_hand", - second_marker="reaching_target", - phase=39, - node=Node.END, - axes=[Axis.X, Axis.Y], -) - -counter = 0 -for i in range(len(pickle_file_list)): - while counter == 0 or sol.status == 0: - if counter == 0: - ocp = OcpFesMsk.prepare_ocp( - biorbd_model_path="../../msk_models/arm26.bioMod", - bound_type="start_end", - bound_data=[[0, 5], [0, 5]], - fes_muscle_models=fes_muscle_models, - n_stim=n_stim, - n_shooting=n_shooting, - final_time=1.5, - pulse_duration={ - "min": minimum_pulse_duration, - "max": 0.0006, - "bimapping": False, - }, - with_residual_torque=False, - custom_constraint=constraint, - activate_force_length_relationship=True, - activate_force_velocity_relationship=True, - minimize_muscle_fatigue=True if pickle_file_list[i] == "minimize_muscle_fatigue" else False, - minimize_muscle_force=True if pickle_file_list[i] == "minimize_muscle_force" else False, - use_sx=False, - ) - else: - with open( - f"result_file/" + "pulse_duration_" + pickle_file_list[i] + "_" + str(counter - 1) + ".pkl", "rb" - ) as file: - data_from_previous_motion = pickle.load(file) - q_init = data_from_previous_motion["states"]["q"][:, :-1] - qdot_init = data_from_previous_motion["states"]["qdot"][:, :-1] - - for j in range(len(ocp.nlp)): - ocp.nlp[j].x_init["q"].init[0] = q_init[0][i * 2] - ocp.nlp[j].x_init["q"].init[1] = q_init[1][i * 2] - ocp.nlp[j].x_init["qdot"].init[0] = qdot_init[0][i * 2] - ocp.nlp[j].x_init["qdot"].init[1] = qdot_init[1][i * 2] - for key in ocp.nlp[i].x_bounds.keys(): - if key == "q" or key == "qdot": - ocp.nlp[0].x_bounds[key].max[0][0] = data_from_previous_motion["states"][key][0][-1] - ocp.nlp[0].x_bounds[key].max[1][0] = data_from_previous_motion["states"][key][1][-1] - ocp.nlp[0].x_bounds[key].min[0][0] = data_from_previous_motion["states"][key][0][-1] - ocp.nlp[0].x_bounds[key].min[1][0] = data_from_previous_motion["states"][key][1][-1] - else: - ocp.nlp[0].x_bounds[key].max[0][0] = data_from_previous_motion["states"][key][-1] - ocp.nlp[0].x_bounds[key].min[0][0] = data_from_previous_motion["states"][key][-1] - - sol = ocp.solve(Solver.IPOPT(_max_iter=10000)) - SolutionToPickle( - sol, "pulse_duration_" + pickle_file_list[i] + "_" + str(counter) + ".pkl", "result_file/" - ).pickle() - counter += 1 - -# [1] Dahmane, R., Djordjevič, S., Šimunič, B., & Valenčič, V. (2005). -# Spatial fiber type distribution in normal human muscle: histochemical and tensiomyographical evaluation. -# Journal of biomechanics, 38(12), 2451-2459. - -# [2] Klein, C. S., Allman, B. L., Marsh, G. D., & Rice, C. L. (2002). -# Muscle size, strength, and bone geometry in the upper limbs of young and old men. -# The Journals of Gerontology Series A: Biological Sciences and Medical Sciences, 57(7), M455-M459. diff --git a/cocofest/examples/dynamics/reaching_task/visualization/pulse_duration_subplot.py b/cocofest/examples/dynamics/reaching_task/visualization/pulse_duration_subplot.py index c0546ad..b02df00 100644 --- a/cocofest/examples/dynamics/reaching_task/visualization/pulse_duration_subplot.py +++ b/cocofest/examples/dynamics/reaching_task/visualization/pulse_duration_subplot.py @@ -25,13 +25,11 @@ nb_stim = len(data_minimize_fatigue["parameters"][pulse_duration_keys[0]]) width = round(data_minimize_fatigue["time"][-1], 2) / nb_stim -pw_for_force_optim = [ - data_minimize_force["parameters"][pulse_duration_keys[i]] * 1000000 for i in range(len(pulse_duration_keys)) -] -pw_for_fatigue_optim = [ - data_minimize_fatigue["parameters"][pulse_duration_keys[i]] * 1000000 for i in range(len(pulse_duration_keys)) -] -pw_list = [pw_for_force_optim, pw_for_fatigue_optim] +pw_data_list = [data_minimize_force["parameters"], data_minimize_fatigue["parameters"]] + +pw_list = [] +for j in range(2): + pw_list.append([pw_data_list[j][pulse_duration_keys[i]] * 1000000 for i in range(len(pulse_duration_keys))]) plasma = cm = plt.get_cmap("plasma") cNorm = colors.Normalize(vmin=100, vmax=600) @@ -57,7 +55,7 @@ def plot_graph(datas): axs[i].set_xlabel("Time (s)") for j in range(nb_stim): - value = pw_for_force_optim[i][j] + value = datas[i][j] color = scalarMap.to_rgba(value) axs[i].barh(muscle_names[i], width, left=j * width, height=0.5, color=color)