#!/usr/bin/env python
# coding: utf-8

# # Interpolacja

# Do interpolacji użyję funkcji `interp1d` (interpolacja jednowymiarowa) z pakietu [scipy](https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html). Żeby móc z niej korzystać trzeba ją „zaimportować” (włączyć do projektu z biblioteki).
# 
# Funkcja `interp1d` jako argumenty przyjmuje dwa dektory danych `x` i `y`. Można jej również określić rodzaj interpolacji. Może to być: `'linear'`, `'nearest'`, `'zero'`, `'slinear'`, `'quadratic'`, `'cubic'` lub stopień wielomianu interpolacyjnego.
# 
# Biblioteka [numpy](https://docs.scipy.org/doc/numpy-1.13.0/reference/) to wygodna biblioteka do operacji na tablicach z danymi.
# 
# [matplotlib.pyplot](https://matplotlib.org/users/pyplot_tutorial.html) to biblioteka oferująca podstawowe funkcje rysowania wykresów (dosyć podobna do tej z Matlaba).

# In[1]:


import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (12,9) # Tu definiujemy rozmiary wykresu. Standardowo wydaje się być zbyt mały.


# In[2]:


import numpy as np
from scipy.interpolate import interp1d
from scipy.interpolate import lagrange
from scipy.interpolate import BarycentricInterpolator


# Nasze done to wartości losowe z zakresu od 0 do 1

# In[3]:


N = 21 # liczba losowanych wartości
M = 5 * N # liczba wartości do wyliczenia funkcji interpolowanej
p = 3 # stopień wielomianu interpolacyjnego


# In[4]:


x0 = np.linspace(-1, 1, N) # Oś x zawiera się między -1 a 1
y0 = np.random.random(N)


# Przygotowujemy wykres

# In[5]:


plt.plot(x0, y0, 'o', label='Data');


# In[6]:


x = np.linspace(-1, 1, M)


# tablica `options` zawiara możliwe typy interpoacji. W pętli będziemy je przebiegać.

# In[7]:


options = ('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic', p)


# In[8]:


for o in options:
    f = interp1d(x0, y0, kind=o)    # interpolacja
    plt.plot(x, f(x), label=o)      # wykres funkcji interpolującej


# In[9]:


plt.legend() # opis wykresu
plt.show()   # generacja wykresu


# Jak widać w każdym z tych przypadków funkcja interpolująca przechodzi przez wygenerowane punkty. między nimi, zaś, może zachowywać się dosyć dziwnie. I odczytywanie danych pomiędzy węzłami interpolacji traktowane jest — na ogół — jako czary albo wróżenie z fusów.

# ## Teraz dane z pomiarów

# In[10]:


data = np.loadtxt("AVERAGE300.dat")


# In[11]:


plt.plot(data[:,0],data[:,1],label='.')


# In[12]:


plt.show()


# ### Fltracja danych
# 
# W tym wypadku filtracja będzie oznaczać odrzucenie danych złych.
# 
# Jak widać dane zawierają błędy. Prawdopodobnie wszystkie wartości mniejsze ok -25 są błędne. Trzeba je usunąć. Funkcja `min` podaje minimalną wartość z tablicy. Funkcja `np.argmin` numer indeksu na którym ta wartość się realizuje.

# In[13]:


print('min =', min(data[:,1]))
print('argmin =', np.argmin(data[:,1]))


# Pojedyńczą wartość można usunąć za pomocą funkcji `np.delete`, ale nie jest to najwygodniejsza metoda gdy złych wartości jest więcej.

# In[14]:


temperatura = np.delete(data, (178), axis=0);
min(temperatura[:,1])
print(data.shape[0]) # to będzie liczba wierszy
print(temperatura.shape[0])


# ### Masked array
# 
# Chyba lepszym pomysłem będzie stworzenie macierzy „z maską” ([*masked array*](https://docs.scipy.org/doc/numpy-1.13.0/reference/maskedarray.generic.html#rationale)). Warto nauczyć się tej techniki.

# In[15]:


import numpy.ma as ma 


# In[16]:


temperatura = ma.masked_less(data[:,1], -25)


# Z pierwszej kolumny wybieramy dane które są mniejsze niż -25 (stopni celsjusza). Powstaje obiekt (`temperatura`) składający się z danych oraz tablicy logicznej zawierających prawdę wszędzie tam gdzie warunke jest spełniony.
# 
# `ma.getmask(temperatura)` (albo `temperatura.mask`) zwraca tablicę logiczną.
# 
# `temperatura.data` (albo `ma.getdata(temperatura)` zwraca „dobre” wartości. 

# In[17]:


ma.getmask(temperatura)


# In[18]:


np.shape(temperatura[~temperatura.mask].data)[0]


# Spróbujmy obejrzeć „zamaskowane” wartości:

# In[19]:


temperatura[temperatura.mask].data


# Teraz musimy wybrać współrzędne czasowe poprawnych wartości z pierwszej kolumny macierzy data. Stworzymy kolejną tablicę z maską (używając maski uzyskanej z temperatury). Funkcja `ma.array` służy do tworzenia tablicy z maską:

# In[20]:


czas=ma.array(data[:,0], mask = temperatura.mask)


# Skoro dostęp do zamaskowanych wartości daje `temperatura.mask` to dostęp do warości poprawnych da negacja, czyli `~temparatura.mask`.

# In[21]:


CZAS = czas[~czas.mask].data # tylko dane
CZAS = CZAS - min(CZAS) # dzięki temu współrzędne x są liczone względem początku doby
TEMPERATURA = temperatura[~temperatura.mask].data # tylko dane


# In[22]:


plt.plot(CZAS[0:30],TEMPERATURA[0:30],label='.')


# In[23]:


f = lagrange(CZAS[0:30], TEMPERATURA[0:30])
x = np.linspace(0,9000, 200)
plt.plot(x, f(x))


# In[24]:


plt.show()


# In[24]:


CZAS[0:30]


# Jak widać ten rodzaj interpolacji, w tym przypadku nie sprawdza się (choć wyniki dla 10 węzłów wyglądają interesująco)

# In[25]:


f = BarycentricInterpolator(CZAS[0:30], TEMPERATURA[0:30])
x = np.linspace(1000,7800, 200)
plt.plot(CZAS[5:30],TEMPERATURA[5:30],label='.')
plt.plot(x, f(x))
plt.show()


# ## Teraz funkcja trygonometryczna (sin)

# In[26]:


from math import pi
from math import sqrt
from math import sin


# In[27]:


X=np.array([0,  pi/6,   pi/4,   pi/3,   pi/2,   2*pi/3, 3*pi/4, 5*pi/6, 
            pi, 7*pi/6, 5*pi/4, 4*pi/3, 3*pi/2, 5*pi/3, 7*pi/4, 11*pi/6, 
            2*pi])


# In[28]:


Y=np.array([0,  1/2,  sqrt(2)/2,  sqrt(3)/2,   1,  sqrt(3)/2,  sqrt(2)/2,  1/2, 
            0, -1/2, -sqrt(2)/2, -sqrt(3)/2,  -1, -sqrt(3)/2, -sqrt(2)/2, -1/2, 
            0])


# In[29]:


my_sin=BarycentricInterpolator(X,Y)


# In[30]:


x=np.linspace(0, 2*pi, 100)


# In[31]:


plt.plot(X,Y)
plt.plot(x,my_sin(x))
plt.show()


# In[32]:


s = np.sin(x)


# In[33]:


plt.plot(x, s-my_sin(x))


# In[35]:


plt.show()


# Zwracam uwagę, że interpolaca jest bardzo dobra! „Przeregulowania” na początku i końcu przedziału co do wartości bezwzględnej są mniejsze niż $1^{-9}$.

# In[ ]: