dimension-reduction/dimensiereductie_AP_veelpunten.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

# AP model (phase_1000).
# Dit bouwt verder op dimensiereductie_AP_verbeterd.py
# door heel veel meer punten te nemen, en maakt gebruik
# van de bevinding dat door epsilon handmatig in te stellen,
# er geen convergentieproblemen ontstaan bij ARPACK.

# Laadt eerst de data in
data = ith.SimData.from_stem("phase_1000/phase_1000")

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



def als_kolomvector(waardes):
    return waardes.reshape((waardes.shape[0], 1))

def onderzoek(aantal_punten):
    begintijd = 2
    var_u = snij_randen_weg(data, 'u', begintijd).flatten()
    var_v = snij_randen_weg(data, 'v', begintijd).flatten()
    # Verlaag aantal punten om geheugengebruik en tijdsduur te beperken.
    # replace=True strikt genomen incorrect maar wel een stuk sneller
    np.random.seed(0)
    keuzes = np.random.choice(var_u.shape[0], aantal_punten, replace=True)
    var_u = var_u[keuzes]
    var_v = var_v[keuzes]
    # Plot het!
    #    print("scatter")
    #    plt.scatter(var_u,var_v)
    #    plt.title("phaseplot of steekproef of AP model (%s points)" % (aantal_punten,)) # todo engels
    #    plt.xlabel("u")
    #    plt.ylabel("v")
    #    plt.savefig("AP %s punten 0.04 eps 2 eigvectoren (scatter).png" % (aantal_punten,))
    #    plt.savefig("AP %s punten 0.04 eps 2 eigvectoren (scatter).pdf" % (aantal_punten,))
    #    plt.close()
    # Maak nu diffusion maps
    faseruimte = np.column_stack((var_u,var_v))
    np.random.seed(100)
    dmap = DiffusionMap.from_sklearn(epsilon=0.04, n_evecs=2) #0.01)
    np.random.seed(200)
    print("fitting")
    dmap.fit(faseruimte)
    print("data")
    fig = diff_visualization.data_plot(dmap,show=False)
    fig.gca().set_title("Data coloured with first DC (AP model, %s points)" % (aantal_punten,)) # todo engels
    plt.xlabel("u")
    plt.ylabel("v")
    fig.savefig("AP %s punten 0.04 eps 2 eigvectoren (data).png" % (aantal_punten,))
    fig.savefig("AP %s punten 0.04 eps 2 eigvectoren (data).pdf" % (aantal_punten,))
    plt.close()
    print("embedding")
    fig = diff_visualization.embedding_plot(dmap,show=False)
    fig.gca().set_title("Embedding given by first two DCs (AP model, %s points)" % (aantal_punten,)) # todo engels
    fig.savefig("AP %s punten 0.04 eps 2 eigvectoren (embedding).png" % (aantal_punten,))
    fig.savefig("AP %s punten 0.04 eps 2 eigvectoren (embedding).pdf" % (aantal_punten,))
    plt.close()

# Voor 7000 moet ik al een tijdje wachten
for aantal_punten in [100, 200, 500, 1000, 2000, 3000, 4000, 5000, 6000]:
    onderzoek(aantal_punten)

# AP 6000 punten 0.04 eps 2 eigvectoren.png: twee deuken in een cirkel, een paar uitschieters
# AP 3000 punten 0.04 eps 2 eigvectoren.png: twee deukjes, een aantal uitschieters
# AP 2000 punten 0.04 eps 2 eigvectoren.png: een appel met twee ogen
# AP 1000 punten 0.04 eps 2 eigvectoren.png: een neus met een paar stippen