Созвездия

Исходные данные

Сначала загрузим необходимые пакеты.

library(maditr)
library(ggplot2)
library(cluster) # для k-medoid
# library(kmed) # для k-medoid, но pam из пакета cluster работает лучше
library(akmeans) # для k-means с косинусной метрикой

Скачаем данные из репозитория Stellarium:

CONSTELLATIONS = "data/constellations.txt"
CONSTELLATION_NAMES = "data/constellation_names.txt"

# файл с координатами звезд
if(!file.exists(CONSTELLATIONS)){
    download.file("https://raw.githubusercontent.com/Stellarium/stellarium/master/skycultures/western_SnT/SnT_constellations.txt", CONSTELLATIONS)
}
# файл с расшифровкой названий созвездий
if(!file.exists(CONSTELLATION_NAMES)){
    download.file("https://raw.githubusercontent.com/Stellarium/stellarium/master/skycultures/western_SnT/constellation_names.eng.fab", CONSTELLATION_NAMES)
}

Загрузка полных названий созвездий:

constellations_names = fread(CONSTELLATION_NAMES, sep = "", header = FALSE) %>% 
    let(
        V1 = trimws(gsub("_.+$", "", V1)),
        V1 = gsub('"', "", V1)
    )

constellations_names = strsplit(constellations_names$V1, "\t") %>% 
    to_df() %>% 
    setnames(c("abbr", "constellation"))

head(constellations_names)

##    abbr constellation
## 1:  And     Andromeda
## 2:  Ant        Antlia
## 3:  Aps          Apus
## 4:  Aqr      Aquarius
## 5:  Aql        Aquila
## 6:  Ara           Ara

Данные в файле с координатами содержатся в позиционном формате, поэтому загрузка данных не совсем прямолинейна. Описание:

1-5   visual magnitude of star
6     <blank>
7-14  RA (hours and decimals, J2000.0)*
15    <blank>
16-23 North polar distance (deg. and dec. J2000.0)*
24    <blank>
25-28 Bayer Greek letter, maybe with superscript
29    code for line weight*
30-32 constellation*

constellations = fread(CONSTELLATIONS, sep = "", header = FALSE, strip.white = FALSE) %>% 
    take_if(!startsWith(V1, "#"),
            magnitude = as.numeric(substr(V1, 1, 5)), 
            ra = as.numeric(substr(V1, 7, 14)), 
            npd = as.numeric(substr(V1, 16, 23)), 
            const_abbr = substr(V1, 30, 32) 
    ) %>% 
    dt_left_join(constellations_names, by = c(const_abbr = "abbr")) %>% 
    na.omit() # у нас несколько пропусков, удаляем их


head(constellations)

##    const_abbr magnitude       ra     npd constellation
## 1:        And      2.06  0.13981 60.9094     Andromeda
## 2:        And      3.28  0.65547 59.1392     Andromeda
## 3:        And      2.06  1.16219 54.3794     Andromeda
## 4:        And      2.10  2.06500 47.6703     Andromeda
## 5:        And      3.62 23.03203 47.6742     Andromeda
## 6:        And      4.30 23.63561 46.7319     Andromeda

summary(constellations)

##   const_abbr          magnitude            ra               npd          
##  Length:969         Min.   :-1.470   Min.   : 0.0637   Min.   :  0.7358  
##  Class :character   1st Qu.: 2.910   1st Qu.: 6.3827   1st Qu.: 60.9094  
##  Mode  :character   Median : 3.650   Median :12.6197   Median : 90.4928  
##                     Mean   : 3.467   Mean   :12.4378   Mean   : 92.2207  
##                     3rd Qu.: 4.200   3rd Qu.:18.1257   3rd Qu.:124.0742  
##                     Max.   : 6.720   Max.   :23.9986   Max.   :173.6678  
##  constellation     
##  Length:969        
##  Class :character  
##  Mode  :character  
##                    
##                    
##

Разведочный анализ

Список всех созвездий:

dt_count(constellations, constellation, sort = TRUE)

##        constellation  n
##  1:         Hercules 32
##  2:            Orion 32
##  3:        Centaurus 31
##  4:       Ursa Major 29
##  5:        Ophiuchus 28
##  6:          Perseus 27
##  7:           Gemini 25
##  8:      Sagittarius 24
##  9:        Andromeda 23
## 10:         Aquarius 23
## 11:            Hydra 22
## 12:              Leo 21
## 13:         Scorpius 21
## 14:            Virgo 21
## 15:           Carina 19
## 16:            Draco 19
## 17:           Aquila 18
## 18:           Bootes 17
## 19:             Pavo 17
## 20:           Pisces 17
## 21:           Cygnus 16
## 22:          Pegasus 16
## 23:          Serpens 16
## 24:      Canis Major 15
## 25:            Cetus 15
## 26:          Lacerta 15
## 27:            Lepus 15
## 28:          Cepheus 14
## 29:        Monoceros 14
## 30:           Taurus 14
## 31:           Auriga 13
## 32:           Crater 13
## 33:            Lupus 13
## 34:           Puppis 13
## 35:   Camelopardalis 10
## 36:      Capricornus 10
## 37:         Eridanus 10
## 38:            Libra 10
## 39: Piscis Austrinus 10
## 40:              Ara  9
## 41:           Dorado  9
## 42:             Grus  9
## 43:             Lyra  9
## 44:          Phoenix  9
## 45:             Vela  9
## 46:          Columba  8
## 47:  Corona Borealis  8
## 48:             Lynx  8
## 49:       Ursa Minor  8
## 50:           Corvus  7
## 51:            Musca  7
## 52:           Tucana  7
## 53:           Volans  7
## 54:       Chamaeleon  6
## 55:           Cancer  6
## 56:             Crux  6
## 57:        Delphinus  6
## 58:       Horologium  6
## 59:            Indus  6
## 60:        Leo Minor  6
## 61:             Apus  5
## 62:            Aries  5
## 63:       Cassiopeia  5
## 64: Corona Australis  5
## 65:           Hydrus  5
## 66:            Norma  5
## 67:        Reticulum  5
## 68:           Scutum  5
## 69:          Sagitta  5
## 70:           Caelum  4
## 71:   Coma Berenices  4
## 72:           Octans  4
## 73:            Pyxis  4
## 74:         Sculptor  4
## 75:          Sextans  4
## 76:       Triangulum  4
## 77:           Antlia  3
## 78:         Circinus  3
## 79:         Equuleus  3
## 80:           Fornax  3
## 81:           Pictor  3
## 82:      Telescopium  3
## 83:      Canis Minor  2
## 84:        Vulpecula  2
##        constellation  n

# количество созвездий
n_const = uniqueN(constellations$constellation)

Всего у нас 84 созвездия. Это число нам понадобится, когда мы будем делать кластерный анализ. Нарисуем импровизированную карту звездного неба. Каждый цвет обозначает свое созвездие.

ggplot(constellations) + 
    geom_point(aes(x = ra, y = npd, alpha = magnitude + 2), color = "orange") +
    geom_point(aes(x = ra, y = npd, color = constellation), size = 4, alpha = .3) +
    geom_path(aes(x = ra, y = npd, group = constellation)) +
    theme_minimal() +
    theme(legend.position="none")

Кластеризация

Так как кластерный анализ не знает реальных названий созвездий, то нам нужна функция, которая будет сопоставлять номер кластера с созвездием. Она должна делать это таким образом, чтобы максимизировать итоговую accuracy. Мы считаем, что сопоставление должно быть один к одному. Алгоритм такой:

Подсчитываем, сколько звезд попало в данное пересечение “созвездие - кластер” для каждой пары
Сортируем пары по убыванию
Идем сверху вниз и берем соотвествие с максимальным количеством звезд. Повторно использовать созвездие или номер кластера нельзя.

Надо обратить внимание, что в таком подходе некоторые кластеры совсем не попадут в цель.

match_clusters = function(target, predicted){
    # функция возвращает вектор с проставленными предсказанными созвездиями
    # там, где созвездие предсказнно неверно, будет NA
    # данный вектор должен максимизировать accuracy
    dt = data.table(target, predicted)
    res =  dt %>% 
        take(n = .N, by = .(target, predicted)) %>% 
        sort_by(-n)  

    all_targets = unique(res$target)
    used_clusters = c()
    for(each_target in all_targets){
        curr_target = query(res, target == each_target & !(predicted %in% used_clusters))
        curr_pred = which.max(res$n[curr_target])[1]
        curr_match = res$predicted[curr_target][curr_pred]
        used_clusters = c(used_clusters, curr_match)
    }
    to_join = na.omit(data.table(predicted = used_clusters, matches = all_targets))
    (anyDuplicated(to_join$predicted) || anyDuplicated(to_join$matches)) && stop("Что-то пошло не так")
    
    dt_left_join(dt,
                 to_join,
                 by = "predicted"
                 ) %>% 
        query(matches)
}

is_correct = function(target, matched){
    # возвращает TRUE для правильно восстановленного target
    # FALSE для прочих. 
    !is.na(matched) & matched == target
}

Делаем три кластерных анализа с разным типом дистанции.

set.seed(123)
constellations = constellations %>% 
    let(
        random_baseline = sample(constellation), # будем сравнивать со случайным угадыванием
        res_random_baseline = is_correct(constellation, random_baseline),
        constant_baseline = "Hercules", # самое большое созвездие
        res_constant_baseline = is_correct(constellation, constant_baseline),
        # евклидово расстояние
        k_euclid = kmeans(cbind(ra, npd), centers = n_const)$cluster,
        const_euclid = match_clusters(constellation, k_euclid),
        res_euclid = is_correct(constellation, const_euclid),
        # k-medoids
        k_medoid = pam(cbind(ra, npd), k = n_const, cluster.only=TRUE),
        # k_medoid = fastkmed(dist(cbind(ra, npd)), ncluster = n_const)$cluster,
        const_medoid = match_clusters(constellation, k_medoid),
        res_medoid = is_correct(constellation, const_medoid),
        # k-means с косинусом
        k_cosine = norm.sim.ksc(cbind(ra, npd), k = n_const)$cluster,
        const_cosine = match_clusters(constellation, k_cosine),
        res_cosine = is_correct(constellation, const_cosine),
        # все то же самое, но со стандартизацией
        # евклидово расстояние
        k_scaled_euclid = kmeans(scale(cbind(ra, npd)), centers = n_const)$cluster,
        const_scaled_euclid = match_clusters(constellation, k_scaled_euclid),
        res_scaled_euclid = is_correct(constellation, const_scaled_euclid),
        # k-medoids
        k_scaled_medoid = pam(scale(cbind(ra, npd)), k = n_const, cluster.only=TRUE),
        # k_scaled_medoid = fastkmed(dist(cbind(ra, npd)), ncluster = n_const)$cluster,
        const_scaled_medoid = match_clusters(constellation, k_scaled_medoid),
        res_scaled_medoid = is_correct(constellation, const_scaled_medoid),
        # k-means с косинусом
        k_scaled_cosine = norm.sim.ksc(scale(cbind(ra, npd)), k = n_const)$cluster,
        const_scaled_cosine = match_clusters(constellation, k_scaled_cosine),
        res_scaled_cosine = is_correct(constellation, const_scaled_cosine)
    )

Подводим итоги. Выведем таблицу с accuracy в процентах:

result = take_all(constellations, if(startsWith(.name, "res_")) round(mean(.x)*100,1)) 
knitr::kable(t(result), col.names = c( "Accuracy, %"))

	Accuracy, %
res_random_baseline	2.2
res_constant_baseline	3.3
res_euclid	39.8
res_medoid	40.2
res_cosine	6.4
res_scaled_euclid	57.2
res_scaled_medoid	60.8
res_scaled_cosine	24.4

Отрисуем картинку - черным закрашены звезды с правильно восстановленным названием созвездия:

ggplot(constellations) + 
    geom_point(aes(x = ra, y = npd, color = ifelse(res_scaled_medoid, "correct", "error")), size = 4, alpha = .5) +
    scale_color_manual(values = c("correct" = "black", "error" = "orange")) +
    geom_path(aes(x = ra, y = npd, group = constellation)) +
    theme_minimal() +
    theme(legend.position="none")

Выводы

Наилучшую точность обеспечивает метод k-medoid по стандартизированным данным: 60.8 процентов. Практически такую же точность дает евклидово расстояние. Косинусная мера сходства сильно отстает. Однако все они значительно превосходят случайное угадывание.

Все это не очень помогает ответить на вопрос от booking.com. Можно использовать любой из трех методов кластеризации, они все дают результат, далекий от случайного. Но этот результат не идеально точный, поэтому возожно, что правильный ответ “Ничего из перечисленного”.

В дальнейших изысканиях можно попытаться получить угловые расстояния между звездами для наблюдателя с Земли и использовать их для кластеризации. Так же могут хорошо сработать методы кластеризации на графах. Однако все это пусть останется для будущих исследователей.

Github с исходным кодом и данными: Constellations

Созвездия

Григорий Демин

2020 06 13

Постановка задачи

Исходные данные

Разведочный анализ

Кластеризация

Выводы