{ "cells": [ { "cell_type": "markdown", "id": "7516a4a2", "metadata": {}, "source": [ "# 800100715151 Astronomide Veritabanları #\n", "\n", "## Ders - 04c Katalogların k-d Ağaçlarıyla Çapraz Eşleştirilmesi##" ] }, { "cell_type": "markdown", "id": "6d2179b1", "metadata": {}, "source": [ "Önemli ölçüde kısıtlanmış kataloglarda bile çapraz eşleştirmenin uzun zaman aldığı açıktır. Bunun nedeni hesap etkinliği bakımından yetersiz bir algoritma olmasıdır. Çapraz eşleştiricinin [uygulanma şekli](Ders4b_Kataloglarin_Capraz_Eslestirilmesi.ipynb), BSS'deki her nesne için SuperCOSMOS'taki her nesneye olan mesafenin hesaplanmasını gerektirir. Bu kısıtlanmış kataloglar arası çapraz eşleştirme bile 160 × 500 = 80.000 mesafe hesabını gerektirmektedir. \n", "\n", "Her bir uzaklık hesabı birkaç mikrosaniye alırken, bu saniyeler veya dakikalar hızla birikir. Saniyeler başta önemsiz görünebilir, ancak kısıtlanmamış tam SuperCOSMOS kataloğunun bir önceki örnekte çalışılandan 250000 kat daha büyük olduğu ve 126 milyon nesne içerdiği düşünülürse hesap sayısı ve dolayısıyla süresinin ne kadar uzayacağını tahmin etmek zor değildir. SuperCOSMOS ile karşılaştırılabilir bir boyuta sahip AT20G BSS'den farklı bir katalog bununla çapraz eşleştirilmeye çalışıldığında bu yük içinden çıkılmaz hale gelebilir. Bunun gibi bir çapraz eşleştirme işlemi kişisel bir bilgisayarda günler hatta aylar alabilir.\n", "\n", "Daha \"akıllı\" bir algoritmaya ihtiyacımız olacağı açıktır. Bir önceki uygulamada geliştirilen [çapraz eşleştiricide](Ders4b_Kataloglarin_Capraz_Eslestirilmesi.ipynb), trigonometrik fonksiyonların doğru çalışabilmesi için $RA$ ve $DEC$'i radyana çevirmek gerekiyordu. Bu dönüşüm her seferinde cisimler arası uzaklığı hesaplayan fonksiyonda yapılırsa, çapraz eşleştirme işlemi sırasında aynı koordinatlar birçok kez dönüştürülmüş olur.\n", "\n", "Bir sonraki problemde, herhangi bir uzaklık hesabından önce dönüşümün yalnızca bir kez gerçekleşmesi için çapraz eşleştirici algoritmasının değiştirilmesi istenecektir. Zamandan yapılan bu küçük tasarruf bile önemlidir ve birikimli olarak toplam kod çalışma süresini etkiler.\n", "\n", "Ayrıca geliştirilmesi istenen algoritma SuperCOSMOS vee AT20G BSS katalogları dışındaki kataloglara da uygulanabilir olmalıdır. Bu daha genel bir çözümdür ve istenen iki kataloğun birleştirilmesi için kullanılabilir.\n", "\n", "Not: Bu uygulama coursera.org 'da University of Sydney tarafından verilen Data Driven Astronomy dersinden adapte edilmiştir." ] }, { "cell_type": "markdown", "id": "18191bc1", "metadata": {}, "source": [ "# Katalogların Yüklenmesi ve Yardımcı Fonksiyonlar:\n", "\n", "Öncelikle çapraz eşleştirme uygulamasının \"naif\" bir [versiyonu](Ders4b_Kataloglarin_Capraz_Eslestirilmesi.ipynb) için yazılan yardımcı fonksiyonlar ve kullanılan katalogların karşılaştırmalar için uygun formatta `Pandas` veriçerçevelerine alınmasını sağlayan `import` fonksiyonlarını alıntılamaya ihtiyaç olacaktır." ] }, { "cell_type": "code", "execution_count": null, "id": "4be05930", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "def import_bss(bss_dosya):\n", " bss = pd.read_fwf(bss_dosya, header=None, usecols=range(1, 7))\n", " bss['RA'] = np.zeros(len(bss[1]))\n", " bss['DEC'] = np.zeros(len(bss[1]))\n", " for i in range(len(bss[1])):\n", " bss.at[i,'RA'] = hms2dec(bss.loc[i,1],bss.loc[i,2],bss.loc[i,3])\n", " bss.at[i,'DEC'] = dms2dec(bss.loc[i,4],bss.loc[i,5],bss.loc[i,6])\n", " bss.drop(range(1,7), axis=1, inplace=True)\n", " return bss\n", "\n", "def import_super(superCOSMOS_dosya):\n", " sc = pd.read_csv(superCOSMOS_dosya, usecols=[0,1])\n", " sc.rename(columns = {'Dec': 'DEC'}, inplace=True)\n", " return sc\n", "\n", "\n", "def hms2dec(h,m,s):\n", " # HH:MM:SS yapisindaki sagacikligi\n", " # derece biriminde decimal formata donusturen fonksiyon\n", " if h > 24 or h < 0:\n", " h %= 24\n", " return 15*(h + m/60 + s/3600)\n", "\n", "def dms2dec(d,m,s):\n", " # DD:MM:SS yapisindaki sagacikligi\n", " # derece biriminde decimal formata donusturen fonksiyon\n", " if d > 90 or d < -90:\n", " d %= 90\n", " if d < 0:\n", " return d - m/60 - s/3600\n", " else:\n", " return d + m/60 + s/3600\n", "\n", "def angular_dist(ra1, dec1, ra2, dec2):\n", " # Koordinatlari radyan biriminde verilen\n", " # İki nokta arasi uzakligi Haversine formuluyle\n", " # hesaplayan fonksiyon\n", " a = np.sin(abs(dec1 - dec2)/2)**2\n", " b = np.cos(dec1)*np.cos(dec2)*np.sin(abs(ra1 - ra2)/2)**2\n", " d = 2*np.arcsin(np.sqrt(a + b))\n", " return d" ] }, { "cell_type": "markdown", "id": "59d03c32", "metadata": {}, "source": [ "## Soru-1:\n", "\n", "### Mikro optimizasyon\n", "\n", "Derece cinsinden $RA$ ve $DEC$ 2-boyutlu dizilerini eşleştirmek üzere `crossmatch` isimli bir fonksiyon yazınız. \n", "\n", "Fonksiyonunuz, çapraz eşleştirmeye başlamadan önce tüm koordinatları radyana çevirmelidir ve aşağıdaki 3 değeri döndürmelidir:\n", "\n", "1. Eşleştirilen cisimlerin ID'leri ve aralarındaki mesafeden oluşan demetlerin bir listesi;\n", "2. Birinci katalogda olup ikinciyle eşleşmeyen cisimlerin ID'lerinden oluşan bir liste;\n", "3. Fonksiyonun çalışmasının kaç saniye sürdüğünü gösteren bir kayan noktalı sayı.\n", "\n", "Her iki katalog da 2-boyutlu `NumPy` dizileri olarak fonksiyona gönderilir. Her satır, tek bir nesnenin koordinatlarını içerir. Verinin ilk iki sütunu $RA$ ve $DEC$ iken cisim ID'leri 0'dan başlayarak sayılan satır numarasıdır. \n", "\n", "`time` modülü fonksiyonlarından (`perf_counter()`) yararlanarak tüm işlemin kaç saniyede yapıldığını hesaplayan bir yapıyı fonksiyonunuza ekleyiniz ve bu hesap süresini de fonksiyondan döndürünüz." ] }, { "cell_type": "code", "execution_count": null, "id": "153afa3f", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Cozum\n", "import numpy as np\n", "import time\n", "def crossmatch(cat1, cat2, tolerans):\n", " start = time.perf_counter()\n", " matched = []\n", " unmatched = []\n", " cat1 = np.radians(cat1)\n", " cat2 = np.radians(cat2)\n", " tolerans = np.radians(tolerans)\n", " for i in range(len(cat1['RA'])):\n", " r = cat1.loc[i,'RA']\n", " d = cat1.loc[i,'DEC']\n", " minsep = angular_dist(r,d,cat2.loc[0,'RA'],cat2.loc[0,'DEC'])\n", " idmin = 0\n", " for j in range(len(cat2['RA'])):\n", " angsep = angular_dist(r,d,cat2.loc[j,'RA'],cat2.loc[j,'DEC'])\n", " if angsep < minsep:\n", " minsep = angsep\n", " idmin = j\n", " if minsep < tolerans:\n", " matched.append((i,idmin,np.degrees(minsep)))\n", " else:\n", " unmatched.append(i)\n", " tot_time = time.perf_counter() - start\n", " return matched,unmatched,tot_time\n", "\n", "# Test fonksiyonu\n", "if __name__ == '__main__':\n", " cat1 = import_bss(\"veri/bss.dat\")\n", " cat2 = import_super(\"veri/superCOSMOS.csv\")\n", " tol = 40 / 3600. # 40 yaysaniyesi\n", " matches, no_matches, time_taken = crossmatch(cat1, cat2, tol)\n", " print('Eslesenler (Ilk 5) :', matches[:5])\n", " print(\"Toplam eslesen koordinat sayisi: \", len(matches))\n", " print('Eslesmeyenler (Ilk 5) :', no_matches[:5])\n", " print(\"Toplam eslesmeyen koordinat sayisi: \", len(no_matches))\n", " print('Toplam sure:', time_taken)" ] }, { "cell_type": "markdown", "id": "c0058030", "metadata": {}, "source": [ "$NumPy$ bu hesapları C ve Fortran dillerinde geliştirilmiş nümerik yapılardan faydalanarak yaptığı için Python listeleri üzerinde çalışmakta olan $for$ döngülerinden çok daha hızlıdır.\n", "\n", "Listelerin elemanlarına tek tek uygulamak yerine tüm bir diziye uygulanan `np.sqrt`, `np.sin`, `np.cos` ve `np.arcsin)` fonksiyonları kullanılarak çapraz eşleştirmenin ikinci katalog içindeki (aşağıda) tarama bölümü \n", "\n", "```\n", "for id2, (ra2, dec2) in enumerate(cat2):\n", " dist = angular_dist_rad(ra1, dec1, ra2, dec2)\n", " if dist < min_dist:\n", " min_dist = dist\n", "```\n", "\n", "aşağıdaki şekilde değiştirilecek olursa:\n", "\n", "```\n", "ra2s = cat2[:, 0]\n", "dec2s = cat2[:, 1]\n", "dists = angular_distance_rad(ra1, dec1, ra2s, dec2s)\n", "min_dist = np.min(dists)\n", "```\n", "\n", "Artık cisimler arası uzaklıkların hesaplanması, ikinci kataloğun tüm $RA$'ları ve $DEC$'leri üzerinde aynı anda gerçekleşir. " ] }, { "cell_type": "markdown", "id": "c0236fc7", "metadata": {}, "source": [ "## Soru-2:\n", "\n", "### Vektörleştirme\n", "\n", "Daha önce yazdığınız `angular_dist` ve `crossmatch` `NumPy` dizileriyle çalışacak şekilde düzenleyiniz (vektörleştiriniz)." ] }, { "cell_type": "code", "execution_count": null, "id": "47f9fce6", "metadata": {}, "outputs": [], "source": [ "# Solution\n", "# Write your crossmatch function here.\n", "import numpy as np\n", "import time\n", "def crossmatch(cat1, cat2, tolerans):\n", " start = time.perf_counter()\n", " matched = []\n", " unmatched = []\n", " cat1 = np.radians(cat1)\n", " cat2 = np.radians(cat2)\n", " tolerans = np.radians(tolerans)\n", " for i in range(len(cat1['RA'])):\n", " r = cat1.loc[i,'RA']\n", " d = cat1.loc[i,'DEC']\n", " dists = angular_dist(r,d,cat2['RA'],cat2['DEC'])\n", " minsep = np.min(dists)\n", " idmin = np.argmin(dists)\n", " if minsep < tolerans:\n", " matched.append((i,idmin,np.degrees(minsep)))\n", " else:\n", " unmatched.append(i)\n", " tot_time = time.perf_counter() - start\n", " return matched,unmatched,tot_time\n", " \n", "# Test fonksiyonu\n", "if __name__ == '__main__':\n", " cat1 = import_bss(\"veri/bss.dat\")\n", " cat2 = import_super(\"veri/superCOSMOS.csv\")\n", " tol = 40 / 3600. # 40 yaysaniyesi\n", " matches, no_matches, time_taken = crossmatch(cat1, cat2, tol)\n", " print('Eslesenler (Ilk 5) :', matches[:5])\n", " print(\"Toplam eslesen koordinat sayisi: \", len(matches))\n", " print('Eslesmeyenler (Ilk 5) :', no_matches[:5])\n", " print(\"Toplam eslesmeyen koordinat sayisi: \", len(no_matches))\n", " print('Toplam sure:', time_taken)" ] }, { "cell_type": "markdown", "id": "40f97946", "metadata": {}, "source": [ "Not: Hantal `for` döngüleri yerine `NumPy` dizileri ile çalışılması toplam çalışma süresini neredeyse 1 / 10'una düşürmüştür!" ] }, { "cell_type": "markdown", "id": "98cc8433", "metadata": {}, "source": [ "## Tekrarlardan kaçınmak\n", "\n", "Yapılabilecek bir diğer optimizasyon, ikinci katalogdaki cisimleri, o anda eşleştirilmekte olan ilk katalog cisminden uzak bir DEC değerine sahipse gözardı etmektir. Bunu yapmak için\n", "\n", "* Diziyi ID yerine DEC'e göre sıralayıp ikinci katalog cisimlerini taramak,\n", "\n", "* Ulaşılan DEC değerleri arananı maksimum toleranstan fazla geçtiğinde aramayı durdurmak\n", "\n", "iyi bir çözümdür. \n", "\n", "Örneğin, $\\delta$ dikaçıklıklı bir ilk katalog cismi (hedef) olduğunu ve maksimum eşleşme yarıçapının da $r$ olarak belirlendiğini varsayalım. Yeni algoritma, döngüden çıkmadan önce yalnızca $-90$ (başlangıç) ile $\\delta + r$ derece arasındaki deklinasyona sahip ikinci katalog cisimleri üzerinde çalışacaktır. Bu da ortalamada kataloğun neredeyse yarısının hiç taranmayacağı anlamına gelir ki önemli bir kazançtır. Bu tekrar döngü yapısını kullanmayı gerektirse de performans açısından fayda sağlayıp sağlamayacağını denetlemekte fayda vardır." ] }, { "cell_type": "markdown", "id": "415615ff", "metadata": {}, "source": [ "## Soru-3:\n", "\n", "### Aramayı daraltma\n", "\n", "`crossmatch` fonksiyonunu 2. kataloğu deklinasyonda sıralayıp sadece tolerans limitleri dahilinde arama yapılacak şekilde düzenleyiniz. Sıralama işini `Pandas` veriçerçevelerine ve serilerine uygulanabilen `sort_values` fonksiyonuyla yapmak çok daha kolaydır; çünkü sıralama sonrası satırların (cisim ID'lerinin) orjinal indeksleri korunur! Ancak ikinci katalogdaki $RA$ ve $DEC$ değerlerini karşılaştırmalar ve tarama için `NumPy` dizilerine dönüştürürken bu indeksleri de bir `NumPy` dizisine almak gerekecektir." ] }, { "cell_type": "code", "execution_count": null, "id": "e37c046b", "metadata": {}, "outputs": [], "source": [ "# Cozum\n", "import numpy as np\n", "import time\n", "def crossmatch(cat1, cat2, tolerans):\n", " start = time.perf_counter()\n", " matched = []\n", " unmatched = []\n", " cat1 = np.radians(cat1)\n", " cat2 = np.radians(cat2)\n", " tolerans = np.radians(tolerans)\n", " cat2_sorted = cat2.sort_values(by='DEC')\n", " # tarama ve karsilastirma icin numpy dizileri\n", " ra2 = np.array(cat2_sorted['RA'])\n", " dec2 = np.array(cat2_sorted['DEC'])\n", " # orjinal indeksler\n", " ind2 = np.array(cat2_sorted.index)\n", " for i in range(len(cat1['RA'])):\n", " r = cat1.loc[i,'RA']\n", " d = cat1.loc[i,'DEC']\n", " minsep = angular_dist(r,d,ra2[0],dec2[0])\n", " idmin = ind2[0]\n", " for j in range(len(ra2)):\n", " angsep = angular_dist(r,d,ra2[j],dec2[j])\n", " if angsep < minsep:\n", " minsep = angsep\n", " idmin = ind2[j]\n", " if dec2[j] > d + tolerans:\n", " break\n", " if minsep < tolerans:\n", " matched.append((i,idmin,np.degrees(minsep)))\n", " else:\n", " unmatched.append(i)\n", " tot_time = time.perf_counter() - start\n", " return matched,unmatched,tot_time\n", " \n", "\n", "# Test fonksiyonu\n", "if __name__ == '__main__':\n", " cat1 = import_bss(\"veri/bss.dat\")\n", " cat2 = import_super(\"veri/superCOSMOS.csv\")\n", " tol = 40 / 3600. # 40 yaysaniyesi\n", " matches, no_matches, time_taken = crossmatch(cat1, cat2, tol)\n", " print('Eslesenler (Ilk 5) :', matches[:5])\n", " print(\"Toplam eslesen koordinat sayisi: \", len(matches))\n", " print('Eslesmeyenler (Ilk 5) :', no_matches[:5])\n", " print(\"Toplam eslesen koordinat sayisi: \", len(no_matches))\n", " print('Toplam sure:', time_taken)" ] }, { "cell_type": "markdown", "id": "be757c8d", "metadata": {}, "source": [ "Görüldüğü gibi her ne kadar full taramaya göre süre açısından önemli bir avantaj sağladıysa da bu çözüm, `for` döngüsü kullanmaksızın `NumPy` dizilerine doğrudan dayanan çözüme göre bir miktar daha uzun sürmüştür. Bu durum aranan değerin daha erken bulunabileceği durumlar için daha kısa çalışma süreleri getirerek uzun vadede pek çok başka katalog eşleştirmeleri için daha avantajlı dahi olabilir. Ayrıca `NumPy` dizilerine dönüşümler gibi bazı ek işlemlerin de koda eklenmiş olduğu dikkate alınmalıdır." ] }, { "cell_type": "markdown", "id": "c1772545", "metadata": {}, "source": [ "## Binary search (İkili Arama)\n", "\n", "Aranan deklinasyon değerini tolerans kadar geçtikten sonra aramayı durdurmakla kalmayıp, aramayı aradığımız cisme mümkün olduğunca yakın başlatarak önceki optimizasyonu daha da geliştirebiliriz. Bunun için\n", "\n", "* İkinci katalog cisimlerini deklinasyon sırasına göre sırala;\n", "* Aramayı ikinci katalogda $\\delta - r$ değerinden büyük deklinasyondan başlat;\n", "* İkinci katalogda aramayı $\\delta + r$ değerinde sonlandır.\n", "\n", "Aramayı bu şekilde daraltmak performansı önemli ölçüde arttıracaktır.\n", "\n", "$[\\delta - r, \\delta + r]$ sınırlarına en yakın ikinci katalog nesnelerini bulmanın hızlı bir yolunu bulmamız gerekiyor, böylece aramaya nereden başlayıp nerede bitirileceği belirlenmiş olur. Aksi takdirde yine tarayarak bu değerlere erişilmiş olunur ki tarama yapılmış olacağı için istenen bu değildir ve performans artışı da getirmeyecektir.\n", "\n", "Bunu yapmanın çok etkin bir yolu \"ikili arama\" (ing. binary search) algoritmasından faydalanmaktır.\n", "\n", "[İkili arama](https://tr.wikipedia.org/wiki/%C4%B0kili_arama_algoritmas%C4%B1) algoritması sıralı bir listede aranan elemanın indeks değerini bulmak için listedeki tüm değerlerle tek tek karşılaştırmadan çok daha etkilidir. İkili arama, listeyi art arda ikiye bölerek sadece aranan nesneyi içerebilecek yarıya ulaşılanan kadar aramayı sürdürür.\n", "\n", "\n", "Aşağıdaki listede 15 sayısının bulunması örnek olarak verilebilir:\n", "\n", "```\n", "# 0 1 2 3 4 5 6 7 8 9\n", "s = [10, 11, 12, 13, 14, 15, 16, 17, 18, 19]\n", "```\n", "\n", "1. Bu listenin ortasındaki değer $s[4] = 14$ 'tür.\n", "2. 14 değeri 15'ten küçük olduğundan, 15 $s[5:10]$ diliminde olmalıdır.\n", "3. $s[5:10]$'un ortası $s[7] = 17$'dir.\n", "4. 17 değeri 15'ten büyüktür, öyleyse 15 $s[5:7]$ diliminde olmalıdır.\n", "5. $s[5:7]$ 'nin orta değeri $s[5] = 15$ 'tir.\n", "6. 15 aradığımız değer olup indeksi 5'tir.\n", "\n", "Bu, 10'dan 19'a kadar olan bir listede 15'i bulmanın dolambaçlı bir yolu gibi görünüyor, ancak yalnızca 3 karşılaştırmanın yapıldığını, ancak tüm liste aransa 6 karşılaştırmanın yapılacağı unutulmamalıdır. Büyük dizilerde bu tür küçük tasarruflar önemli fark yaratabilir. Doğrudan arama ile 1000 uzunluğundaki bir listedeki bir öğeyi bulmak için ortalama olarak 500 karşılaştırma gerekliyken, ikili aramada yalnızca 10 karşılaştırma gereklidir.\n", "\n", "\n", "## Numpy ve İkili Arama Algoritması\n", "\n", "`NumPy` modülünün `searchsorted` fonksiyonu (aynı fonksiyon `Pandas` modülünde de bulunmaktadır) bu algoritmayla arama yapar. Fonksiyonun nasıl kullanıldığı aşağıda örneklenmiştir." ] }, { "cell_type": "code", "execution_count": null, "id": "d2d882ef", "metadata": {}, "outputs": [], "source": [ "######0 1 2 3 4 5 6 7 8 9\n", "s = [10, 11, 12, 13, 14, 15, 16, 17, 18, 19]\n", "index = np.searchsorted(s, 15, side='left')\n", "print(index)" ] }, { "cell_type": "code", "execution_count": null, "id": "4d149635", "metadata": {}, "outputs": [], "source": [ "################0 1 2 3 4 5 6 7 8 9\n", "s = pd.Series([10, 11, 12, 13, 14, 15, 16, 17, 18, 19])\n", "index = s.searchsorted(15, side='left')\n", "print(index)" ] }, { "cell_type": "markdown", "id": "ceb88547", "metadata": {}, "source": [ "Aynı dizide aranan değerden birden fazla varsa (örneğin 15 değerinden) `side='left'` ayarı yapıldığında ilk bulunan değerin (15)'in, `side='right'` yapılırsa son bulunan değerin (15) indeksi döndürülür." ] }, { "cell_type": "markdown", "id": "44d73d63", "metadata": {}, "source": [ "## Soru-4:\n", "\n", "### Aramayı Daraltma\n", "\n", "`crossmatch` fonksiyonunu ikinci katalogda sadece $dec1 - tolerans$ ile $dec1 + tolerans$ arasındaki değerleri ikili arama yapacak şekilde düzenleyiniz. `searchsorted` fonksiyonu hem aramanın yapılacağı $dec1 - tolerans$, hem de döngüden çıkılacak $dec1 + tolerans$ değerlerini bulmak için kullanılmalıdır." ] }, { "cell_type": "code", "execution_count": null, "id": "a4980279", "metadata": {}, "outputs": [], "source": [ "# Cozum\n", "import numpy as np\n", "import time\n", "def crossmatch(cat1, cat2, tolerans):\n", " start = time.perf_counter()\n", " matched = []\n", " unmatched = []\n", " cat1 = np.radians(cat1)\n", " cat2 = np.radians(cat2)\n", " tolerans = np.radians(tolerans)\n", " cat2_sorted = cat2.sort_values(by='DEC')\n", " ra2 = np.array(cat2_sorted['RA'])\n", " dec2 = np.array(cat2_sorted['DEC'])\n", " ind2 = np.array(cat2_sorted.index)\n", " for i in range(len(cat1['RA'])):\n", " r = cat1.loc[i,'RA']\n", " d = cat1.loc[i,'DEC']\n", " index_start = np.searchsorted(dec2, d - tolerans,side='left')\n", " index_stop = np.searchsorted(dec2, d + tolerans, side=\"left\")\n", " minsep = angular_dist(r,d,ra2[index_start],dec2[index_start])\n", " idmin = ind2[index_start]\n", " for j in range(index_start,index_stop):\n", " angsep = angular_dist(r,d,ra2[j],dec2[j])\n", " if angsep < minsep:\n", " minsep = angsep\n", " idmin = ind2[j]\n", " if minsep < tolerans:\n", " matched.append((i,idmin,np.degrees(minsep)))\n", " else:\n", " unmatched.append(i)\n", " tot_time = time.perf_counter() - start\n", " return matched,unmatched,tot_time\n", " \n", "# Test fonksiyonu\n", "if __name__ == '__main__':\n", " cat1 = import_bss(\"veri/bss.dat\")\n", " cat2 = import_super(\"veri/superCOSMOS.csv\")\n", " tol = 40 / 3600. # 40 yaysaniyesi\n", " matches, no_matches, time_taken = crossmatch(cat1, cat2, tol)\n", " print('Eslesenler (Ilk 5) :', matches[:5])\n", " print(\"Toplam eslesen koordinat sayisi: \", len(matches))\n", " print('Eslesmeyenler (Ilk 5) :', no_matches[:5])\n", " print(\"Toplam eslesmeyen koordinat sayisi: \", len(no_matches))\n", " print('Toplam sure:', time_taken)" ] }, { "cell_type": "markdown", "id": "3a8eb0c5", "metadata": {}, "source": [ "Çalışma süresi bu kez önemli miktarda azalmıştır (1/25 civarında). Büyük boyutlu katalogların karşılaştırılmasında bu algoritmanın ne kadar büyük bir avantaj sağlayacağı açıktır." ] }, { "cell_type": "markdown", "id": "25cfd310", "metadata": {}, "source": [ "# k-d Ağaçları\n", "\n", "Çapraz eşleştirme, astronomide sıklıkla yapıaln bir işlem olduğundan optimize edilmiş uygulamalarının zaten kodlanmış olması doğaldır. `Astropy` modülü [k-d ağaçlarına](https://en.wikipedia.org/wiki/K-d_tree) dayalı bir çapraz korelasyon fonksiyonunu bu amaçla sunmaktadadır.\n", "\n", "Astropy, içinde arama yapılan ikinci katalogdan bir k-d ağacı oluşturarak, ilk katalogdaki her cisim için ikili aramaya benzer bir şekilde bir eşleşme arar. Bu algoritmada k-boyutlu uzay, her seferinde bölümün bir tarafı yalnızca tek bir cisim içerene kadar özyinelemeli (recursive) olarak iki parçaya bölünür.\n", "\n", "Katalog eşleştirme uygulamasında bu algoritma;\n", "\n", "1. Sağaçıklığın ortanca (medyan) değerini bul, kataloğu bunun sağı ve solundan ikiye böl.\n", "\n", "2. Her bölümde deklinasyonun ortanca değerini bul, bölümü bunun aşağı ve yukarısından ikiye böl.\n", "\n", "3. Bu bölümlerin her birinde sağaçıklığın ortanca (medyan) değerini bul, bu bölümü bunun sağı ve solundan ikiye böl.\n", "\n", "4. 2 ve 3. adımları her bölümde tek bir nesne kalıncaya kadar sürdür.\n", "\n", "Bu algoritma içinde arama yapılacak aşağıdaki şekilde bir ikili ağaç oluşturur.\n", "\n", "
\n", " \n", "
\n", "\n", "Yukarıdaki şekilde 10 cisimden 6. sı olan $A$ sağ açıklıktaki medyan (ortanca) olarak seçilebilir (5 de seçilebilirdi ama sağdaki (daha büyük sağ açıklık değerine sahip olan tercih edilir.). Buradan yapılacak bir bölme sonrası, bu değerin sol ya da sağ tarafına geçilerek ilerlenebilir. Sol tarafa geçilecek olunursa oradaki 5 koordinatın deklinasyonda ortancası $E$'dir. Buradan bölündükten sonra üst ya da alt tarafa geçilerek devam edilebilir. Bu bölmenin üst tarafına geçilecek olursa iki koordinat kaldığı görülür ($D$,$B$); medyan olarak bunlardan sağ açıklığı büyük olan $B$ koordinatı seçilir. Burada deklinasyon ekseninde ilerlenebilecek bir düğüm kalmamıştır; zira $D$ tek başınadır ve alt bölmeye geçilir. Alt bölmenin sağ açıklıktaki medyanı sağ açıklığı büyük olan $F$ olarak belirlenir. Daha sonra $A$'nın sağ tarafına geçilir. Bu bölmede bulunan 4 koordinattan deklinasyonda daha büyük olan $G$ medyan olarak belirlenir. Bunun üzerinde tek bir nokta kalmış olduğundan alt bölmeye inilir ve sağ açıklıkta daha büyük olan $H$ medyan olarak belirlenir. Bunun sol tarafında sadece bir koordinat kaldığından ağaç yapısı tamamlanmış olur.\n", "\n", "## k-d Ağaçları İçinde Arama\n", "\n", "k-d ağacı oluşturulduktan sonra arama aşağıdaki şekilde yapılır.\n", "\n", "1. Aranan cisimle en yüksek seviyedeki düğüm (kök düğüm) arasındaki uzaklığı hesapla, ardından bu düğüme sağ açıklıkta en yakın alt düğüme git,\n", "\n", "2. Aranan cisimle bu düğümdeki (child) cisim arasındaki uzaklığı hesapla, ardından bu alt düğüme dik açıklıkta en yakın alt düğüme git,\n", "\n", "3. Arnanan cisimle bu düğümdeki (child) cisim arasındaki uzaklığı hesapla, ardından alt düğüme sağ açıklıkta en yakın alt düğüme git,\n", "\n", "4. 2-3 arasındaki adımları artık alt düğüm (child) kalmayıncaya kadar tekrarla ve en alt düğüme (leaf node) ulaş,\n", "\n", "5. Hesaplanan tüm uzaklıkların en kısasını bul, bu uzaklık aranan cisme ikinci katalogdaki en yakın cismin uzaklığıdır.\n", "\n", "6. Her bir düğüm iki \"çocuk\" düğüme ayrıldığı için\n", "\n", "$N$ cisim için \"kökten\" (root) \"yaprağa\" (leaf) $log_2(N) $ düğümde inilir. Bu şekilde SuperCOSMOS kataloğu gibi 250 million cisim içeren bir kataloda yalnız 28 uzaklık hesabı yeterli olmaktadır! \n", "\n", "Aşağıdaki şekilde koordinatları bilinen $T$ cismine k-d ağacına dönüştürülen bir katalogdaki en yakın cisim aranmaktadır.\n", "\n", "
\n", " \n", "
\n", "\n", "$T$ cismi ile öncelikle kök düğümdeki cisim ($A$) arasındaki uzaklık hesaplanır. Daha sonra $A$'ya sağaçıklıkta en yakın alt düğümdeki cisme ($E$) gidilir ve bununla $T$ cismi arasındaki mesafe hesaplanır. Bu uzaklık $A$'ya uzaklığa göre küçük olduğundan $E$'de kalınır ve $E$'ye deklinasyonda en yakın cismin bulunduğu alt düğüme geçilir. Bu cisim $B$'dir. $B$ ile $T$ arasındaki uzaklık $E$ ile olandan küçüktür ve devam edilir. Bu kez $B$'ye sağ açıklıkça en yakın alt düğüme gidilir. Bu alt düğümde $D$ bulunmaktadır. $D$ ile $T$ arasındaki mesafe $B$ ile olana göre büyüktür. $T$ cismine k-d ağacına dönüştürülmüş katalogdaki en yakın cisim $B$ olarak 4 karşılaştırma 5 uzaklık ölçümü sonrası bulunmuş olur. Eğer bu cisim istenen bir yarıçap içindeyse $T$ cismi ile eşleşen bir cisim olarak değerlendirilebilir.\n", "\n", "`Astropy` ile bu şekilde yapılan bir arama aşağıda örnek olarak verilmiştir." ] }, { "cell_type": "code", "execution_count": null, "id": "97fb5785", "metadata": {}, "outputs": [], "source": [ "from astropy.coordinates import SkyCoord\n", "from astropy import units as u\n", "coords1 = [[270, -30], [185, 15], [120, -10], [50, 25], [300,5]]\n", "coords2 = [[185, 20], [280, -30], [45, 20], [295, 0], [100, -15]]\n", "sky_cat1 = SkyCoord(coords1*u.degree, frame='icrs')\n", "sky_cat2 = SkyCoord(coords2*u.degree, frame='icrs')\n", "closest_ids, closest_dists, closest_dists3d = sky_cat1.match_to_catalog_sky(sky_cat2)\n", "print(closest_ids)\n", "print(closest_dists)" ] }, { "cell_type": "markdown", "id": "6650aa76", "metadata": {}, "source": [ "`SkyCoord` nesnesi, `Astropy`'ın gökyüzü kataloğu nesnesidir. Birimleri (burada derece cinsinden `u.degree` ile verilmiştir) ve bir referans çerçevesi (burada [ICRS](https://en.wikipedia.org/wiki/)) belirtildiği sürece her koordinat dizisiyle çalışırlar. International Celestial Reference System esasen ekvatoral koordinatlardan oluşan bir katalogdur.`closest_id` ve `closest_dists` çıktıları, `sky_cat1`'de olup `sky_cat2` 'deki bir cisim eşleşen cismin `sky_cat2` içindeki satır indeksini ve ona olan uzaklığını verir.`closest_dists` açısal uzaklık `closest_dists3d` ise 3-boyuttaki uzaklıktır.\n", "\n", "#### Not:\n", "\n", "`Astropy` uzaklık değerlerini `Quantity` nesnesi olarak döndürür. Bu nesne türü bir `NumPy` dizisine çevrilirken bu nesnenin değerini `value` metodu ile almak yeterlidir; birimi de akılda tutmakta fayda vardır.\n", "\n", "```\n", "closest_dists_array = closest_dists.value\n", "```" ] }, { "cell_type": "markdown", "id": "7860f715", "metadata": {}, "source": [ "## Soru-5\n", "\n", "`crossmatch` fonksiyonunu `Astropy` fonksiyonalitesiyle arama yapacak şekilde düzenleyiniz." ] }, { "cell_type": "code", "execution_count": null, "id": "8b8060f7", "metadata": {}, "outputs": [], "source": [ "# Cozum\n", "import numpy as np\n", "import time\n", "from astropy.coordinates import SkyCoord\n", "from astropy import units as u\n", "def crossmatch(cat1, cat2, tolerans):\n", " start = time.perf_counter()\n", " matched = []\n", " unmatched = []\n", " cat1 = SkyCoord(ra=cat1['RA']*u.degree, dec=cat1['DEC']*u.degree, frame=\"icrs\")\n", " cat2 = SkyCoord(ra=cat2['RA']*u.degree, dec=cat2['DEC']*u.degree, frame=\"icrs\")\n", " closest_ids, closest_dists, _ = cat1.match_to_catalog_sky(cat2)\n", " for i,j in enumerate(closest_ids):\n", " angsep = closest_dists[i].value\n", " if angsep < tolerans:\n", " matched.append((i, j, angsep))\n", " else:\n", " unmatched.append(i)\n", " tot_time = time.perf_counter() - start\n", " return matched,unmatched,tot_time\n", "\n", "# Test fonksiyonu\n", "if __name__ == '__main__':\n", " cat1 = import_bss(\"veri/bss.dat\")\n", " cat2 = import_super(\"veri/superCOSMOS.csv\")\n", " tol = 40 / 3600. # 40 yaysaniyesi\n", " matches, no_matches, time_taken = crossmatch(cat1, cat2, tol)\n", " print('Eslesenler (Ilk 5) :', matches[:5])\n", " print(\"Toplam eslesen koordinat sayisi: \", len(matches))\n", " print('Eslesmeyenler (Ilk 5) :', no_matches[:5])\n", " print(\"Toplam eslesen koordinat sayisi: \", len(no_matches))\n", " print('Toplam sure:', time_taken)" ] }, { "cell_type": "markdown", "id": "4624efcd", "metadata": {}, "source": [ "Sonuç olarak hızlı olmakla birlikte ikili aramaya dayanan bir önceki çözümden neredeyse 5 kat yavaş sürdü ancak daha uzun boyutlu katalogların karşılaştırılması için kullanıldığında bir performans avantajı sağlayacaktır. Kataloglar $RA$ ya da $DEC$'e göre sıraladığınızda daha hızlı sonuç elde edildiğini deneyerek görebilirsiniz!" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.12" } }, "nbformat": 4, "nbformat_minor": 5 }