<jupyter_start><jupyter_text>Identifikace systému, odhad parametrů modeluTento notebook je interaktivní, každou buňku můžete vykonat stisknutím `Shift-Enter`, můžete změnit její obsah a opětovně ji vykonat.V kódu jsou tzv. elipsy, buď jako tři tečky `...` nebo podtržítko `_` v místě, kdě je potřeba doplnit nebo naprogramovat část dle zadání úkolu.  Příprava v Modelice- na fyzické mašině, naklonujte si `Physiolibrary-models` nejlépe do stejného adresáře jako je Bodylight-notebooks nebo Bodylight-VirtualMachine `git clone https://github.com/creative-connections/Physiolibrary-models`- ve virtuální mašině v OMEdit (nebo Dymola) otevřete `/vagrant_data/Physiolibrary_models/Metabolism/package.mo`- vytvořte GlucoseToleranceTest1, který v čase 50h od začátku simulace bude pumpovat glukózu rychlostí 100 g za hodinu po dobu 20 minut, přepočtěte na mg/s a spojte je s GIRegulationBlock, simulujte 3 dny, zobrazte koncetraci glukózy a inzulínu (out1 a out2)- vytvořte GlucoseToleranceTest2, který v čase 50h od začátku simulace bude pumpovat glukózu rychlostí 100 g za hodinu po dobu 20 minut, přepočtěte na kg/s a spojte je s GIRegulationComponent, simulujte 3 dny- zkuste měnit parametry beta a nu (na 10% původní hodnoty), simulujte 3 dny- exportujte GlucoseToleranceTest2 od FMU (mód ME), soubor dejte do adresáře `Seminar8FmiIdentification/fmus/`   1 ProblémBudeme zpracovávat data měření tzv. glukózového tolerančního testu. Máme informace, že pacient dostával intravenózně glukózu do krve rychlostí 28 mg/s. po dobu 12 minut. Pak byly v 12 minutových intervalech prováděny testy koncentrace glukózy a inzulínu hodnoty jsou v souborech `data/PatientX_y.csv` kde `X` je označení pacienta a `y` je buď `g`lukóza nebo `i`nzulín.<jupyter_code># nacteme soubor s daty glukozy pro pacienta A 
import numpy
data = numpy.genfromtxt('data/PatientA_g.csv',delimiter=',')
data
# casovou osu vygenerujeme, dle zadani to jde po 12 minutach (tj. po 720 sekundach)
time = numpy.arange(0,720*data.size,720)
time
# definujeme rutinu pro kresleni grafu
def plot(x,y,z=None,labely=None,labelz=None):
    %matplotlib inline
    import matplotlib.pyplot as plt
    fig = plt.figure()
    plt.plot(x,y,label=labely,color='red')
    if z is not None:
        plt.plot(x,z, label=labelz,color='blue')
    plt.legend()
    plt.grid(True)
    plt.show()

plot(time,data,labely="data patient A glucose")<jupyter_output><empty_output><jupyter_text>2. FmPy nebo PyFMI knihovny a simulace modelu v Pythonu Model z Modeliky exportovat jako FMU (nejlépe verze 2.0 v módu ModelExchange nebo Co-Simulation). V OpenModelice 1.16.x doporučujeme exportovat v módu ModelExchange. V Dymole doporučujeme exportovat v módu Co-Simulation (solver DASSL je robustní).Knihovny FmPy a PyFMI zjednodušují simulace FMU v Pythonu. FmPy implementuje solvery v Pythonu. Knihovna PyFMI používá implementace z jiných knihoven v C++, které jsou obecně rychlejší.V dalším textu budeme používat FmPy.<jupyter_code>import fmpy
model1 = 'fmus/Metabolism.GlucoseInsulin.Test.GlucoseToleranceTest.fmu' 
result1 = result = fmpy.simulate_fmu(
    model1,
    start_time=0,
    stop_time=259200,
    step_size=720,
    record_events=False,
    output=['glucoseInsulinRegulation.glucoseConc','glucoseInsulinRegulation.insulinConc'],
    start_values={'glucoseInsulinRegulation.tissueUtilizationInsulinDependent.Nu_permu': 139000,'glucoseInsulinRegulation.insulinProduction.beta_mu': 1430},
    solver='CVode')
fmpy.util.plot_result(result1)<jupyter_output><empty_output><jupyter_text>Jestli se objeví chyba, je třeba nainstalovat fmpy, např.: `conda install -c conda-forge fmpy`.  3. Kalibrace dat a simulace pro odhad parametrů<jupyter_code>plot(result1['time'],result1['glucoseInsulinRegulation.glucoseConc'],labely='model')
plot(time,data,labely="data patient A glucose")<jupyter_output><empty_output><jupyter_text>3.1 data a simulace modelu jsou vzorkována stejně? Ano - interval 12 minut (720s) 3.2 Odpovídají vzorky simulace modelu odpovídají vzorkům dat podle simulovaného experimentu?<jupyter_code>print('data size',time.size,' model simulation size:',result1['time'].size)<jupyter_output>data size 10  model simulation size: 361<jupyter_text>Musíme vybrat začátek, kde simulace modelu začne odpovídat tomu co se děje v reálných datech<jupyter_code>model = result1['glucoseInsulinRegulation.glucoseConc'][249:259]
plot(time,data,model,labely='real data',labelz='model')
cindex = 249 # calibrating index


# Porovnani modelu a simulace jako funkce
def simulatemodelbeta(relativebeta=1,modelfile=model1):
    # odsimulovat s parametrem beta*relativebeta
    modelresult = fmpy.simulate_fmu(
        modelfile,
        start_time=0,
        stop_time=259200,
        step_size=720,
        record_events=False,
        output=['glucoseInsulinRegulation.glucoseConc','glucoseInsulinRegulation.insulinConc'],
        start_values={
            'glucoseInsulinRegulation.tissueUtilizationInsulinDependent.Nu_permu': 139000,
            'glucoseInsulinRegulation.insulinProduction.beta_mu': 1430*relativebeta
        },
        solver='CVode')
    return modelresult;

def comparemodeldatabeta(relativebeta=1,realdata=data,realtime=time,modelfile=model1):
    modelresult = simulatemodelbeta(relativebeta,modelfile)
    # vykreslit do grafu
    plot(realtime,realdata,modelresult['glucoseInsulinRegulation.glucoseConc'][cindex:cindex+10],labely='real data',labelz='model')
    


comparemodeldatabeta(0.1)<jupyter_output><empty_output><jupyter_text>Zkusíme najit hrubou silou parametr `beta` tak, aby model odpovídal co nejvíc datům - rozdíl mezi modelem a daty byl minimální.Nejprve definujeme metriku rozdílu mezi simulací modelu a daty, tj. účelová funkce (objective function) $ y=\sum_{i=1}^n (s_i - d_i)^2 $ kde $s_i$ je hodnota sledované veličiny simulace v bodu $i$ a $d_i$ je hodnota sledované veličiny naměřených dat.<jupyter_code>def mydiff(s,d):
    sum = 0
    for i in range(0,s.size):
        sum+=(s[i]-d[i])**2
    return sum<jupyter_output><empty_output><jupyter_text>Nyní definujeme cyklus, který odsimuluje model s parametry 1 promile až 1000 promile normální hodnoty beta.<jupyter_code># hrubou silou odsimulujeme vsech 1000 moznych hodnot
#diffs = [] # pole s hodnotou ucelove funkci a parametrem
#for i in range(1,1000):
#    modelresult = simulatemodelbeta(i/1000)
#    diff= mydiff(modelresult['glucoseInsulinRegulation.glucoseConc'][cindex:cindex+10],data)
#    diffs.append([i/1000,diff])
# v poli diffs mam nyni parametr a hodnotu ucelove funkce (vzdalenost mezi simulaci a realnymi daty)
# varianta, pokud predchozi vypocet spadl, hrubou silou zkusime procenta
diffs = [] # pole s hodnotou ucelove funkci a parametrem
for i in range(1,100):
    modelresult = simulatemodelbeta(i/100)
    diff= mydiff(modelresult['glucoseInsulinRegulation.glucoseConc'][cindex:cindex+10],data)
    diffs.append([i/100,diff])

# muzeme zobrazit jak procento beta ovlivnuje rozdil mezi daty a modelem
ndiffs = numpy.array(diffs)
plot(ndiffs[:,0],ndiffs[:,1],labely="diff between data and model")<jupyter_output><empty_output><jupyter_text>Hodnota parametru s nejnizsi hodnotou ucelove funkce je hledany odhad parametru, ktery odpovida realnym datum.<jupyter_code>diffs.sort(key=lambda x:x[1])<jupyter_output><empty_output><jupyter_text>Nejlepsi shoda modelu s daty je prvni (v Pythonu nulty prvek setrideneho pole).<jupyter_code>modelbeta = diffs[0][0]
print('Nejlepsi shoda modelu s daty je pro parametr beta (v podilu k normalni hodnote)= ',modelbeta)
# jak vypada srovnani dat a modelu (hodnota parametru je v procentech - musime vydelit 100)
comparemodeldatabeta(diffs[0][0])
# jak vypada srovnani dat a modelu s 2 nejlepsi hodnotou parametru?
print(diffs[1][0])
comparemodeldatabeta(diffs[1][0])
# Jak je to s parametrem ny (nu)
# Porovnani modelu a simulace jako funkce
def simulatemodelnu(relativenu=1,modelfile=model1):
    # odsimulovat s parametrem beta*relativebeta
    modelresult = fmpy.simulate_fmu(
        modelfile,
        start_time=0,
        stop_time=259200,
        step_size=720,
        record_events=False,
        output=['glucoseInsulinRegulation.glucoseConc','glucoseInsulinRegulation.insulinConc'],
        start_values={
            'glucoseInsulinRegulation.tissueUtilizationInsulinDependent.Nu_permu': 139000*relativenu,
            'glucoseInsulinRegulation.insulinProduction.beta_mu': 1430
        },
        solver='CVode')
    return modelresult;
def comparemodeldatanu(relativenu=1,realdata=data,realtime=time,modelfile=model1):
    modelresult = simulatemodelnu(relativenu,modelfile)
    # vykreslit do grafu
    plot(realtime,realdata,modelresult['glucoseInsulinRegulation.glucoseConc'][cindex:cindex+10],labely='real data',labelz='model')
diffs2 = [] # pole s hodnotou ucelove funkci a parametrem
for i in range(1,100):
    modelresult = simulatemodelnu(i/100)
    diff2= mydiff(modelresult['glucoseInsulinRegulation.glucoseConc'][cindex:cindex+10],data)
    diffs2.append([i/100,diff2])
# muzeme zobrazit jak procento beta ovlivnuje rozdil mezi daty a modelem
ndiffs2 = numpy.array(diffs2)
plot(ndiffs2[:,0],ndiffs2[:,1],labely="diff between data and model(nu)")    
diffs2.sort(key=lambda x:x[1])
print('Nejlepsi shoda dat a modelu s parametrem nu v procentech s normou= ',diffs2[0][0])<jupyter_output><empty_output><jupyter_text>InterpretaceSnížená hodnota (10% normy) parametru beta v modelu může svědčit pro diagnózu diabetu 1.typu (snížená produkce inzulínu beta buněk).Snížená hodnota (10% normy) parametru nu v modelu může svědčit pro diagnózu diabetu 2.typu (snížená citlivost tkání na inzulin)Musíme vzít v potaz i naměřená data hladiny inzulinu během glukózového tolerančního testu.<jupyter_code>datai = numpy.genfromtxt('data/PatientA_i.csv',delimiter=',')
plot(time,datai,labely='insulin data')
# zobrazime hladinu insulinu modelu, v normalnim fyziologickem stavu
plot(result1['time'],result1['glucoseInsulinRegulation.insulinConc'],labely='model')
result2 = simulatemodelnu(0.12)
plot(result2['time'],result2['glucoseInsulinRegulation.insulinConc'],labely='model')
plot(time,datai,result2['glucoseInsulinRegulation.insulinConc'][cindex:cindex+10],labely='real data',labelz='model')
result2['glucoseInsulinRegulation.insulinConc'][cindex:cindex+10]<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup