Doç. Dr. Özgür Baştürk
Ankara Üniversitesi, Astronomi ve Uzay Bilimleri Bölümü
obasturk at ankara.edu.tr
http://ozgur.astrotux.org
Bu derste neler öğreneceksiniz?¶
İnterpolasyon¶
İnterpolasyon¶
Sayısal çözümlemede interpolasyon, bir süreksiz veri setinde bilinen değerlerden yola çıkarak bilinmeyen ara değerlerin hesaplanmasına dayanan her türlü yöntem olarak tanımlanır. Astronomi, interpolasyonun en sık ihtiyaç duyulduğu alanlardan biridir. Zira astronomide sürekli veri almak ancak uzay teleskoplarıyla mümkün olabilmekte, onlarda dahi sürekli veri alınmasına engel kısıtlar bulunmaktadır. Gece / gündüz, kapalı / açık hava, aletsel problemler, kullanılan dedektörlerin veri alabilme sıklığı, çok bant gözlem yapma ihtiyacı, uzay teleskoplarıın veriyi Dünya'ya iletmek için veri alımına ara vermesi gibi nedenlerle astronomide kullanılan verinin içerisinde büyük boşluklar oluşabilmektedir. Bu boşlukları, var olan veriden yola çıkarak doldurmak önemli bir ihtiyaç olarak karşımıza çıkmaktadır.
İnterpolasyonda izlenen geleneksel yöntemler, örneklerle birlikte İnterpolasyon Yöntemleri kaynağında verilmiştir. Ayrıca astronoimde en sık başvurulan interpolasyon yöntemi olan </i>Spline İnterpolasyonu</i> dokümanında ayrıntılı olarak anlatılmş ve örneklendirilmiştir. Bu derste Python programlama dilinin ve ek modüllerin interpolasyona ilişkin sağladığı olanaklar örnek kodlarla işlenecektir. Dersten önce bu iki dokümanın incelenmesinde fayda vardır.
scipy.Interpolate Modülü ile İnterpolasyon¶
1 Boyutlu İnterpolasyon: interp1d¶
scipy.interpolate modülü fonksiyonlarından interp1d verilen verideki birbirini takip eden her bir $(x_i,y_i)$ ikilisinin arasına bir doğru uyumlanmasına dayanır. Veri noktası bulunmayan $x_k$ ara değerleri değerleri için bu doğruların karşılık geldiği $y_k$ değerleri hesaplanarak aradeğer hesabı (interpolasyon) yapılır.
Aşağıdaki örnekte ilk olarak $y = e^{\frac{-x}{3}}$ fonksiyonu $x = [0,10)$ aralığında 50 nokta için oluşturulmaktadır. Daha sonra bu $x_i,y_i$ ikilileri interpolate.interp1d fonksiyonuna sağlanarak $f$ doğrusal interpolasyon fonksiyonu oluşturulmuştur. $x_{yeni} = [1,8]$ aralığında oluşturulan 5 yeni $x_k$ değeri bu fonksiyona sağlanmış, bu noktalar için $y_k$ aradeğerleri elde edilmiştir. Bu ara değerler ve $(x_i,y_i)$ ikilileri arasında interpolasyonun üzerinden yapıldığı doğrular, veriyi oluşturmak için kullanılan fonksiyonla birlikte çizdirilmiştir.
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate
x = np.arange(0,10)
y = np.exp(-x/3.0) # y = e^(-x / 3)
f = interpolate.interp1d(x,y)
xyeni = np.linspace(1,8,5)
yyeni = f(xyeni)
plt.plot(x,y,'b-')
plt.plot(xyeni,yyeni,'ro')
plt.plot(xyeni,yyeni,'r-')
plt.xlabel("x")
plt.ylabel("y")
plt.show()
Aksi söyllenmediği takdirde fonksiyon, veri noktaları arasında birer doğru uyumlayacak şekilde oluşturulur. Eğer kind parametresi cubic şekilde belirlenmişse bu kez noktaların arasına 3. dereceden birer polinom uyarlanır. Bu yapısıyla fonksiyon gerçekte lineer ya da kübik spline interpolasyonu uygulamaktadır.
Aşağıdaki örnekte bu kez önce $y = cos(-x^2 / 9.0)$ fonksiyonu $x = [0,10]$ aralığındaki 11 nokta için oluşturulmakta; daha sonra, ardışık her $(x,y)$ ikilisine sırasıyla doğrusal ve 3. dereceden polinomlar aynı aralıktaki 41 nokta için uyarlanmaktadır. Burada dikkat edilmesi gereken husus, veri setninin tamamına doğrusal ya da 3. dereceden bir polinomun uyarlanmadığı, bu fonksiyonların veri setindeki ardışık ikililerin arasına uyarlandığıdır.
from matplotlib import pyplot as plt
import numpy as np
from scipy.interpolate import interp1d
x = np.linspace(0, 10, 11)
y = np.cos(-x**2/9.0)
f1 = interp1d(x,y)
f2 = interp1d(x,y,kind='cubic')
xyeni = np.linspace(0, 10, 41)
y1 = f1(xyeni)
y2 = f2(xyeni)
plt.plot(x, y, 'o', xyeni, y1, '-', xyeni, y2, '-.')
plt.legend(['veri', 'dogrusal', 'kubik'], loc = 'best')
plt.plot
plt.show()
Çok Boyutlu İnterpolasyon: griddata¶
scipy.interpolated sadece bir boyutta tanımlı $y = f(x)$ yapısındaki veriler için değil çok boyutlu $z = f(x,y,\omega,...)$ şeklindeki veriler için de interpolasyon yapmaya yönelik olanaklar sağlamaktadır. Öncelikle veriyi oluşturmak üzere
fonksiyonunu tanımlayalım
import numpy as np
def fonk(x, y):
return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2
$z$ fonksiyonunun $[0, 1] x [0, 1]$ gridindeki çok sayıda nokta için interpolasyonunu hesaplamak üzere değerini bu griddeki rastgele 1000 nokta için bildiğimizi varsayalım. Aşağıda np.random.rand fonksiyonu ile ortalama değeri $\mu = 0$, standart sapması $\sigma = 1$ olan birer normal dağılımdan oluşturulan $noktalar$ dizisi iki boyutludur ve fonksiyona gönderilen hem $x$, hem de $y$ değerlerini tutmaktadır. Fonksiyon bu ikililer için hesaplanan $z$ değerlerini döndürür ve böylece arasında interpolasyon yapmak istediğimiz $z$ noktaları oluşturulmuş olur.
noktalar = np.random.rand(1000, 2)
degerler = fonk(noktalar[:,0], noktalar[:,1])
Bu noktalardan hareketle istediğimiz $(x,y)$ ikilileri için bu kez $z_{interp.}$ interpolasyon değerlerini hesaplayalım. Bunun için interpolasyon yapmak istedğimiz $(x,y)$ ikililerini np.mgrid fonksiyonuyla $[0,1]$ arasında $x$ ekseninde 100, $y$ ekseninde 200 nokta içerecek şekilde oluşturalım. ve interpolate.griddata fonksiyonuna orjinal verimizi içeren $noktalar$ ve $degerler$ dizilerinin yanı sıra, interpolasyon yapılmasını istedğimiz bu $(x,y)$ ikililerini ve interpolasyon yöntemimizi sağlayalım.
from scipy.interpolate import griddata
grid_x, grid_y = np.mgrid[0:1:100j, 0:1:200j]
grid_z0 = griddata(noktalar, degerler, (grid_x, grid_y), method='nearest')
grid_z1 = griddata(noktalar, degerler, (grid_x, grid_y), method='linear')
grid_z2 = griddata(noktalar, degerler, (grid_x, grid_y), method='cubic')
Son olarak orjinal fonksiyon ve interpolasyonla elde ettiğimiz noktaları bir grafikle görelim.
import matplotlib.pyplot as plt
%matplotlib inline
plt.subplot(221)
plt.imshow(fonk(grid_x, grid_y).T, extent=(0,1,0,1), origin='lower')
plt.plot(noktalar[:,0], noktalar[:,1], 'k.', ms=1)
plt.title('Veri')
plt.subplot(222)
plt.imshow(grid_z0.T, extent=(0,1,0,1), origin='lower')
plt.title('En Yakin')
plt.subplot(223)
plt.imshow(grid_z1.T, extent=(0,1,0,1), origin='lower')
plt.title('Dogrusal')
plt.subplot(224)
plt.imshow(grid_z2.T, extent=(0,1,0,1), origin='lower')
plt.title('Kubik')
plt.gcf().set_size_inches(6, 6)
plt.show()
Yapısal Spline İnterpolasyonu¶
Yapısal (ing. structural) spline interpolasyonu iki adımda gerçekleştirilir: i) verinin spline interpolasyonu ile temsili oluşturulur, ii) istenen noktalarda spline interpolasyonu hesabı yapılır. Spline temsili için gereken katsayılar doğrudan ya da parametrik olarak elde edilebilir. Verinin doğrudan spline temsili (katsayıların hesabı) splrep fonksiyonu ile gerçekleştirilir. Fonksiyona verinin $(x,y)$ noktaları geçirilir; fonksiyon $t$: düğüm noktalarını, $c$: katsayıları, $k$: spline derecesini göstermek üzere $(t, c, k)$ demet değişkenini döndürür. Varsayılan spline derecesi kübiktir ancak istenirse $k$ parametresi ile değişitirilebilir.
N-boyutlu uzayda veri $splprep$ fonksiyonu ile parametrik olarak temsil edilir. Bu fonksiyon tek bir değişken kabul ettiği için N-boyutlu veri N-boyutlu bir dizi ile tanımlanır. Dizilerin uzunluğu verideki nokta sayısını, dizilerin sayısı ise verinin kaç boyut içerdiğini gösterir. Parametre değişkeni $u$ parametresi ile gösterilir ve $[0-1]$ aralığında monotonik artan bir dizidir. Fonksiyon yine spline parametrelerini tutan $(t, c, k)$ demetini ve $u$ parametresini döndürür.
$s$ interpolasyon sırasındaki düzleştirmenin (ing. smoothing) derecesini belirlemek üzere kulanılabilecek bir parametredir. Varsayılan değeri $m$: veri sayısını göstermek üzere $s = m - \sqrt{2m}$ ile verilir. Düzleştirme istenmiyorsa $s = 0$ olduğu belirtilmelidir.
Spline temsili belirlendikten sonra istenen noktalarda spline interpolasyon (ve türevi ile integrailinin) hesabını yapan fonksiyonlar $splev$, $splde$, $splint$ fonksiyonlarıdır. 8 ya da daha fazla düğüm içeren kübik spline fonksiyonları için kök hesabı yapan bir de $splroot$ fonksiyonu bulunmaktadır.
Bir Boyutlu Yapısal Spline İnterpolasyonu¶
Aşağıda $y = sin(x)$ fonksiyonuyla oluşturulan veri üzerinden bu yöntem kullanılarak yapılan bir kübik interpolasyon örneği verilmiştir. Görüldüğü gibi orjinal verideki $(x,y)$ ikilileri scipy.interpolate.splev fonksiyonuna sağlanarak spline temsili oluşturulmakta, bu temsil interpolasyon yapılmak istenen $x_{yeni}$ noktaları ile scipy.interpolate.splev fonksiyonuna sağlanmakta ve bu fonksiyon verilen parametrelerle bu interpolasyon işlemini $x_{yeni}$ noktalarında gerçekleştirerek $y_{yeni}$ değerlerini hesaplamaktadır. $[0, 2\pi]$ aralığında sadece $\pi / 4$ aralıkla 6 nokta kullanılarak oluşturulan veri üzerinde dahi kübik interpolasyonun $sin(x)$ fonksiyonunu başarıyla oluşturduğu görülmektedir.
import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
np.random.seed(1234)
x = np.arange(0, 2*np.pi+np.pi/4, 2*np.pi/8)
y = np.sin(x) + np.random.randn(len(x))*0.1
tck = interpolate.splrep(x, y, s=0)
xyeni = np.arange(0, 2*np.pi, np.pi/50)
aradeger = interpolate.splev(xyeni, tck, der=0)
plt.figure()
plt.plot(x,y,'go',label="veri")
plt.plot(xyeni,aradeger,'b-',label="Kubik spline")
plt.plot(xyeni,np.sin(xyeni), 'r-', label="sin(x)")
plt.legend(loc="upper right")
plt.axis([-0.05, 6.33, -1.25, 1.25])
plt.title('Kubik Spline Interpolasyonu')
plt.show()
İki Boyutlu Yapısal Spline İnterpolasyonu¶
İki boyutlu bir yüzey üzerindeki veri için ara değer hesabı yapmak (interpolasyon) üzere interpolate.bsplrep fonksiyonu kullanılır. Bu fonksiyona tıpkı 1-boyutlu yapısal (procedural) spline interpolasyonunda olduğu gibi $z = f(x,y)$ yüzeyini oluşturan x, y ve z dizileri gönderilir. Fonksiyon $tx, ty, tz$: düğüm noktalarını, $c$: katsayıları, $kx$, $ky$: spline derecesini göstermek üzere $(tx, ty, tz, c, kx, ky)$ demet değişkenini döndürür. Bu demet değişken interpolate.splev fonksiyonuna geçirilir ve istenen $(x, y)$ noktaları için interpolasyon (ara değer) hesabı yapılır.
Tıpkı 1-boyutlu spline interpolasyonunda olduğu gibi 2-boyutlu spline interpolasyounda da $s$ parametresi düzleştirme (ing. smoothing) için kullanılır. $s = 0$ düzleştirme yapılmak istenmediği anlamına gelir. $m$: nokta sayısını göstermek üzere varsayılan düzleştirme değeri $s = m - \sqrt{2m}$ 'dir.
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
# 20x20 'lik bir grid uzerinde bir fonksiyon tanimlayalim
x, y = np.mgrid[-1:1:20j, -1:1:20j]
z = (x+y) * np.exp(-6.0*(x*x+y*y))
# grafik
plt.figure()
plt.pcolor(x, y, z)
plt.colorbar()
plt.title("Seyrek Orneklenmis Fonksiyon")
plt.show()
# 70x70'lik daha genis bir grid uzerinden interpolasyonla ornekleme
xyeni, yyeni = np.mgrid[-1:1:70j, -1:1:70j]
tck = interpolate.bisplrep(x, y, z, s=0, kx=3, ky=3)
zyeni = interpolate.bisplev(xyeni[:,0], yyeni[0,:], tck)
# grafik
plt.figure()
plt.pcolor(xyeni, yyeni, zyeni)
plt.colorbar()
plt.title("Spline Interpolasyonu")
plt.show()
UnivariateSpline İle Nesne Yönelimli Spline İnterpolasyonu¶
Python'da spline interpolasyonu için nesne yönelimli (ing. object oriented) bir alternatif de bulunmaktadır. Bu alternatifte başlatıcı (\__init__) metoduna $(x,y)$ verisi geçirilerek bir $spline$ nesnesi oluşturulur. Sınıfın \__call__ metodu gönderilen $x$ değerleri için spline interpolasyonu ile y değerlerini hesaplar ve döndürür. InterpolatedUnivariateSpline alt sınıfı bu işlemi yapmasının yanı sıra ayrıca integral, türev ve kök bulma için de metodlara sahiptir.
Bunun yanı sıra düzleştirme istendiğinde tıpkı yapısal alternatifte olduğu gibi $s$ parametresi kullanılır. Bu durumda verideki düğüm sayısından daha az (düzleştirme miktarı kadar az) düğüm içeren bir spline fonksiyonu oluşturulur. Bu şekilde düzleştirme kullanılarak oluşturulan spline fonksiyonuna “interpole spline” yerine “düzleştirilmiş spline” adı verilir.
UnivariateSpline sınıfının bir başka alt sınıfı olan LSQUnivariateSpline $t$ parametresi ile düğümleri dışarıdan yerleri ile birlikte sağlamak için kullanılabilir. Aralarında eşit uzaklıklar olmayan ve her noktadan geçmeyen; bazı aralıklarda düzleştirme içeren bazılarında içermeyen özel nitelikli spline fonksiyonları yaratmak için kullanışlıdır.
Şimdi bir önceki örnekte yapısal yöntemle gerçekleştirdiğimizi interpolasyonu bir de nesne yönelimli yöntemle gerçekleştirelim.
import numpy as np
from scipy import interpolate
from matplotlib import pyplot as plt
np.random.seed(1234)
x = np.arange(0, 2*np.pi+np.pi/4, 2*np.pi/8)
y = np.sin(x) + np.random.randn(len(x))*0.1
s = interpolate.InterpolatedUnivariateSpline(x, y, k=3)
xyeni = np.arange(0, 2*np.pi, np.pi/50)
yyeni = s(xyeni)
plt.plot(x,y,'go',label="veri")
plt.plot(xyeni,yyeni,'b+',label="Kubik spline")
plt.legend(loc="upper right")
plt.axis([-0.05, 6.33, -1.25, 1.25])
plt.title('Kubik Spline (Nesne Yonelimli) Interpolasyonu')
plt.show()
Tek bir nokta için de (sözgelimi $x = 5$ için) spline ile aradeğer bulunabilir.
x = 5
print(s(x))
Örnek: Spline Fonksiyonları ile İnterpolasyon¶
Spline fonksiyonlarında türev, integral ve kök bulmaya yönelik genel bir örnek aşağıdaki kodla verilmektedir. Koda ilişkin açıklamalar iç dokümatasyonda "#" ile sınırlanan yorum ifadeleri ile kod içinden sağlanmıştır.
Spline interpolasyonu sırasında nokta ikilileri arasına yapılan 3. dereceden polinom uyumlamalarının (fonksiyonlarının) türev ve integrallerinin verinin oluşturulduğu sinüs foksiyonunun türev ve integraliyle uyumunu dikkatle inceleyiniz. Her iki eğri de genel yapı olarak benzer şekilde ilerlerken, bir nokta aralığından diğerine değişen 3. derece polinomlarda genel yapı bir kenara bırakılırsa bir aralıktan diğerine hem türev, hem de integral de genel yapıdan sapmalar olduğuna dikkat ediniz.
import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
np.random.seed(1234)
# Kubik Spline
x = np.arange(0, 2*np.pi+np.pi/4, 2*np.pi/8)
y = np.sin(x)+np.random.randn(len(x))*0.1
tck = interpolate.splrep(x, y, s=0, k=3)
xyeni = np.arange(0, 2*np.pi, np.pi/50)
yyeni = interpolate.splev(xyeni, tck, der=0)
# Grafik
plt.figure()
plt.plot(x, y, 'x', xyeni, yyeni, xyeni, np.sin(xyeni))
plt.legend(['Veri', 'Kubik Spline', 'sin(x)'])
plt.axis([-0.05, 6.33, -1.25, 1.25])
plt.title('Kubik Spline Interpolasyonu')
plt.show()
# Turevi
yturev = interpolate.splev(xyeni, tck, der=1)
plt.figure()
plt.plot(xyeni, yturev, xyeni, np.cos(xyeni),'--')
plt.legend(['Kubik Spline', 'cos(x)'])
plt.axis([-0.05, 6.33, -1.25, 1.25])
plt.title('Spline Ile Hesaplanan Turev')
plt.show()
# Integrali
def integ(x, tck, constant=-1):
x = np.atleast_1d(x)
sonuc = np.zeros(x.shape, dtype=x.dtype)
for n in range(len(sonuc)):
sonuc[n] = interpolate.splint(0, x[n], tck)
sonuc += constant
return sonuc
yint = integ(xyeni, tck)
plt.figure()
plt.plot(xyeni, yint, xyeni, -np.cos(xyeni), '--')
plt.legend(['Kubik Spline', 'Veri'])
plt.axis([-0.05, 6.33, -1.25, 1.25])
plt.title('Spline Ile Hesaplanan Integral')
plt.show()
# Kokleri
interpolate.sproot(tck)
# array([3.1416])
# sproot fonksiyonunun x = 0 kokunu bulamadigina dikkat ediniz!
# bu koku bulabilmek icin spline'i daha genis bir aralik uzerinden
# hesaplamaniz gerekir!
x = np.linspace(-np.pi/4, 2.*np.pi + np.pi/4, 21)
y = np.sin(x)
tck = interpolate.splrep(x, y, s=0)
interpolate.sproot(tck)
Kaynaklar¶
Data Reduction and Error Analysis for the Physical Sciences, 3rd ed., Philip R.Bevington & D. Keith Robinson, 2003, McGraw Hill Higher Education
Tanımlar için: wikipedia Interpolasyon başlığı