dimension-reduction/dimensiereductie_courtemanche.py

#!/bin/env python3
import ithildin as ith
import numpy as np
import matplotlib.pyplot as plt
import sys
from typing import List
from scipy.interpolate import RectBivariateSpline
from pydiffmap.diffusion_map import DiffusionMap
from pydiffmap import visualization as diff_visualization
import matplotlib

# Courtemanche1998 model, 1 ruimtedimensie, 20 variabelen

data = ith.SimData.from_stem("myokit12/myokit_1")

# Vermijd randeffecten en initialiatie
def snij_randen_weg(variabele,begintijd, eindtijd):
    def begin(index):
        if index == 0:
            return begintijd
        return variabele.shape[index]//8
    def eind(index):
        if index == 0:
            return eindtijd
        if variabele.shape[index] < 20:
            return variabele.shape[index]
        return (7 * variabele.shape[index])//8
    return variabele[begin(0):eind(0),begin(1):eind(1),begin(2):eind(2),begin(3):eind(3)]

def snij_randen_weg2(variables,begintijd,eindtijd):
    newdata = dict()
    for key in variables.keys():
        newdata[key] = snij_randen_weg(variables[key], begintijd,eindtijd)
    return newdata

def onderzoek(t0, t1, aantal_punten, eps):
    vars = snij_randen_weg2(data.vars, t0, t1)
    # We bekijken de faseruimte, tijds- en ruimtecoordinaten zijn daarin niet van belang.
    for key in vars.keys():
        vars[key] = vars[key].ravel()
    # Verlaag aantal punten om geheugengebruik en tijdsduur te beperken.
    # replace=True strikt genomen incorrect maar verlaagt geheugengebruik
    # en wegens de hoge hoeveelheid punten weinig belangrijk
    np.random.seed(1)
    keuzes = np.random.choice(vars['u'].shape[0], aantal_punten,replace=True)
    for key in data.vars.keys():
        vars[key] = vars[key][keuzes]
    # Nu de randen verwijderd zijn en we punten in een faseruimte hebben, kunnen we proberen
    # diffusion maps te gebruiken.
    dmap = DiffusionMap.from_sklearn(epsilon=eps,n_evecs=4)
    dmap.fit(np.column_stack((vars['gateui'], vars['gatexs'], vars['gated'], vars['Ca'], vars['gatew'], vars['gateu'], vars['Na'], vars['gateoa'], vars['Caup'], vars['gatef'], vars['gateua'], vars['gateh'], vars['gatexr'], vars['K'], vars['gatefCa'], vars['u'], vars['gatem'], vars['gatev'], vars['gatej'], vars['gateoi'], vars['Carel'])))
    diff_visualization.embedding_plot(dmap,dim=3)
    return dmap


# Zoals ‘Embedding given by first three DCs, coloured by fourth (BOFC model, %s points)’,
# maar plot ook de punten die niet gebruikt worden om de diffusion map op te stellen.
# Hopelijk maakt dat de figuur wat duidelijker ...
def plot_veel_punten(t0, t1, aantal_punten, dmap):
    vars = snij_randen_weg2(data.vars, t0, t1)
    # copied from 'onderzoek'
    for key in vars.keys():
        vars[key] = vars[key].ravel()
    print("raveled")
    np.random.seed(1) # determinism
    keuzes = np.random.choice(vars['u'].shape[0], aantal_punten, replace=True)
    # calculate the diffusion coordinates # TODO waarom vind Emacs+Python de o-trema niet goed? 
    var_phi1234 = dmap.transform(np.column_stack((vars['gateui'], vars['gatexs'], vars['gated'], vars['Ca'], vars['gatew'], vars['gateu'], vars['Na'], vars['gateoa'], vars['Caup'], vars['gatef'], vars['gateua'], vars['gateh'], vars['gatexr'], vars['K'], vars['gatefCa'], vars['u'], vars['gatem'], vars['gatev'], vars['gatej'], vars['gateoi'], vars['Carel'])))
    print("transformed")
    # plot it!
    fig = plt.figure(figsize=(6,6))
    ax = fig.add_subplot(111,projection='3d')
    ax.scatter(var_phi1234[:,0],var_phi1234[:,1],var_phi1234[:,2],var_phi1234[:,3],c=var_phi1234[:,3],cmap='viridis')
    ax.set_title("Embedding given by first three DCs, coloured by fourth (Courtemanche model, 1D)") # todo engels
    ax.set_xlabel(r'$\psi_1$')
    ax.set_ylabel(r'$\psi_2$')
    ax.set_zlabel(r'$\psi_3$')
    plt.axis('tight') # ? copied from pydiffmap
    fig.savefig("Courtemanche1D 3+1 eigenvectors, extended embedding (uvw, s).png")
    #fig.savefig("BOFC 3+1 eigenvectors, extended embedding (uvw, s).pdf")
    plt.show()
    plt.close()

# vanaf t < 50 komen er losse puntjes (misschien te ...)
matplotlib.use("TkAgg")
dmap = onderzoek(50, 151, 7000, 0.2)
plot_veel_punten(50, 151, 9, dmap)
#del var1, var2, var3, var4 # geheugengebruik beperken

#dmap = DiffusionMap.from_sklearn(n_evecs = 4)
#map.fit(faseruimte)