Spaces:
Runtime error
Runtime error
| from fmpy import * | |
| from fmpy import read_model_description, extract | |
| from fmpy.fmi2 import FMU2Slave | |
| import numpy as np | |
| import shutil | |
| import pandas as pd | |
| import random | |
| import plotly.graph_objects as go | |
| def add_noise(min, max, std, mean): | |
| return np.clip(np.random.normal(mean, std), min, max) | |
| def gen_pulsive_noise(count, min_val, max_val): | |
| pulse = 0 | |
| n = count // 100 # Use integer division to get the floor of count / 100 | |
| start = 100 * n | |
| end = start + 50 | |
| if (count > start and count < end and n % 15 == 0): | |
| pulse = random.randint(min_val, max_val) | |
| return pulse | |
| def simulation(input_noise, output_noise, min, max, std, mean, min_pul, max_pul, kp, ki, kd, show_original): | |
| vrs = {} | |
| fmu = 'Pharmacokinetics_4_comportmental_model_PI_ref_FMU_base4_OAAS_lnx.fmu' | |
| model_description = read_model_description(fmu) | |
| for variable in model_description.modelVariables: | |
| vrs[variable.name] = variable.valueReference | |
| start_time = 0.0 | |
| stop_time = 7000 | |
| step_size = 1 | |
| unzipdir = extract(fmu) | |
| fmu = FMU2Slave(guid=model_description.guid, | |
| unzipDirectory=unzipdir, | |
| modelIdentifier=model_description.coSimulation.modelIdentifier, | |
| instanceName='instance1') | |
| # initialize | |
| fmu.instantiate() | |
| fmu.setupExperiment(startTime=start_time) | |
| fmu.enterInitializationMode() | |
| fmu.exitInitializationMode() | |
| fmu.setReal([vrs["amesim_interface.Age_year"]], [60]) | |
| fmu.setReal([vrs["amesim_interface.BIS0"]], [95.6]) | |
| fmu.setReal([vrs["amesim_interface.BISmin"]], [8.9]) | |
| fmu.setReal([vrs["amesim_interface.Drug_concentration_mgmL"]], [20]) | |
| fmu.setReal([vrs["amesim_interface.EC50"]], [2.23]) | |
| fmu.setReal([vrs["amesim_interface.Gamma"]], [1.58]) | |
| fmu.setReal([vrs["amesim_interface.Gender_0male_1female"]], [1]) | |
| fmu.setReal([vrs["amesim_interface.Height_cm"]], [168]) | |
| fmu.setReal([vrs["amesim_interface.Infusion_rate_mLh"]], [200]) | |
| fmu.setReal([vrs["amesim_interface.Weight_kg"]], [75]) | |
| vr_input = vrs["amesim_interface.Infusion_rate_mLh"] | |
| vr_output = vrs["amesim_interface.BIS_Index"] | |
| rows = [] # list to record the results | |
| time = start_time | |
| infusion_rate = 200 | |
| i = 0 | |
| target = 30 | |
| last_error = 0 | |
| # simulation loop | |
| while time < stop_time: | |
| if time >= 2.4e3 and time < 4.5e3: | |
| target = 70 | |
| p = 0 | |
| i = 0 | |
| if time >= 4.5e3: | |
| target = 30 | |
| p = 0 | |
| i = 0 | |
| bis = fmu.getReal([int(vr_output)])[0] if time > step_size else 95.6 | |
| if input_noise: | |
| bis += add_noise(min, max, std, mean) | |
| p = bis - target | |
| i = i + p | |
| d = p - last_error | |
| last_error = p | |
| infusion_rate = np.clip(kp*p + ki*i + kd*d, 0, 200) | |
| if output_noise: | |
| infusion_rate += gen_pulsive_noise(time, min_pul, max_pul) | |
| fmu.setReal([vr_input], [int(infusion_rate)]) | |
| # perform one step | |
| fmu.doStep(currentCommunicationPoint=time, communicationStepSize=step_size) | |
| # advance the time | |
| time += step_size | |
| # get the values for 'inputs' and 'outputs[4]' | |
| inputs, outputs = fmu.getReal([int(vr_input), int(vr_output)]) | |
| # append the results | |
| rows.append((time, bis, inputs)) | |
| fmu.terminate() | |
| fmu.freeInstance() | |
| shutil.rmtree(unzipdir, ignore_errors=True) | |
| result = np.array(rows, dtype=np.dtype([('time', np.float64), ('BIS', np.float64), ('Infusion', np.float64)])) | |
| result_baseline = np.load("result.npy") | |
| df = pd.DataFrame(result) | |
| df_original = pd.DataFrame(result_baseline) | |
| trace1 = go.Scatter(x=df.index, y=df['BIS'], mode='lines', name='BIS') | |
| fig1 = go.Figure(data=trace1) | |
| fig1.update_layout(height=400, width=1200, title_text="BIS evolution") | |
| # Add a line trace for column_2 in the second subplot | |
| trace2 = go.Scatter(x=df.index, y=df['Infusion'], mode='lines', name='Infusion') | |
| fig2 = go.Figure(data=trace2) | |
| fig2.update_layout(height=400, width=1200, title_text="Infusion evolution") | |
| if show_original: | |
| fig1.add_trace(go.Scatter(x=df_original.index, y=df_original['BIS'], mode='lines', name='BIS original', line=dict(color="red"), opacity=0.5)) | |
| fig2.add_trace(go.Scatter(x=df_original.index, y=df_original['Infusion'], mode='lines', name='Infusion original', line=dict(color="red"), opacity=0.5)) | |
| return fig1, fig2 | |
| import gradio as gr | |
| with gr.Blocks() as demo: | |
| with gr.Row(): | |
| with gr.Column(scale=1): | |
| gr.Markdown("**BIS noise**") | |
| input_noise = gr.inputs.Checkbox(label="Add Gaussian noise") | |
| with gr.Accordion("noise range"): | |
| min_gaussian = gr.inputs.Slider(minimum=0, maximum=20, step=1, default=5, label="noise min") | |
| max_gaussian = gr.inputs.Slider(minimum=0, maximum=20, step=1, default=10, label="noise max") | |
| std_gaussian = gr.inputs.Slider(minimum=0, maximum=10, step=1, default=2, label="noise standard deviation") | |
| mean_gaussian = gr.inputs.Slider(minimum=0, maximum=20, step=1, default=7, label="noise mean") | |
| gr.Markdown("**Infusion noise**") | |
| with gr.Blocks(): | |
| output_noise = gr.inputs.Checkbox(label="Add Pulsive noise") | |
| with gr.Accordion("noise range"): | |
| min_pul = gr.inputs.Slider(minimum=0, maximum=50, step=1, default=50, label="noise min") | |
| max_pul = gr.inputs.Slider(minimum=0, maximum=150, step=1, default=150, label="noise max") | |
| gr.Markdown("**PID controller paramters**") | |
| with gr.Blocks(): | |
| kp_slider = gr.inputs.Slider(minimum=0, maximum=20, default=4, label="kp") | |
| ki_slider = gr.inputs.Slider(minimum=0, maximum=1, default=0.01, label="ki") | |
| kd_slider = gr.inputs.Slider(minimum=0, maximum=200, default=120, label="kd") | |
| show_original = gr.inputs.Checkbox(label="Show original") | |
| button = gr.Button("Simulate") | |
| with gr.Column(scale=5): | |
| plot1 = gr.Plot(label="BIS evolution") | |
| plot2 = gr.Plot(label="Infusion evolution") | |
| button.click(simulation, inputs=[input_noise, output_noise, min_gaussian, max_gaussian, std_gaussian, mean_gaussian, min_pul, max_pul, kp_slider, ki_slider, kd_slider, show_original], outputs=[plot1, plot2]) | |
| demo.launch() |