DjangoBB LoFi version

Полная версия: Дистанция хи-квадрат для scipy.spatial.distance.pdist

Начало » Центр помощи » Дистанция хи-квадрат для scipy.spatial.distance.pdist

oldbay

Апрель 2, 2015 21:10:41

Добрый день
Делаю алгоритм горной кластеризации для анализа многослойных сетов растровых геоданных.
Для вычислении растояния в алгоритме использую scipy.spatial.distance pdist и cdist. К сожалению, набор измерений дистанции по умолчанию у scipy не даёт нужного эффета - нужен метод определения подобия кривых хи-квадрат описан тут.
pdist и cdist позволяют производить обсчёт дистанций методами отличными от стандартных документация по scipy.spatial.distance.pdist пункт 23
В результате сделал следующий тестовый скрипт по хи-квадрат, показавший хорошие результаты:

import numpy as np
from scipy.spatial.distance import pdist, squareform
a = np.array(
    [                                                                                          
        [1,1,1,2],                                                                             
        [3,3,3,5],                                                                             
        [1,2,3,4],                                                                             
        [2,4,6,8],
        [1,2,3,2],
        [3,4,5,4],
    ]
)
X2_1=lambda u, v: np.sqrt(
        (
            (u - np.mean(u)*(u+v)/(np.mean(u)+np.mean(v)))**2/(
                np.mean(u)*(u+v)/(np.mean(u)+np.mean(v)))
        ).sum() + (
            (v - np.mean(v)*(u+v)/(np.mean(u)+np.mean(v)))**2/(
                np.mean(v)*(u+v)/(np.mean(u)+np.mean(v)))
        ).sum()
    )
b = pdist(
    a,
    X2_1
)
c = pdist(
    a,
    "correlation"
)
print a
print squareform(b)
print squareform(c)

Для сравнения рассчитал расстояния для того же массива с применением более-менее близкого метода correlation
массив

[[1 1 1 2]
 [3 3 3 5]
 [1 2 3 4]
 [2 4 6 8]
 [1 2 3 2]
 [3 4 5 4]]

расстояния хи-квадрат

[[ 0.          0.17053338  0.61237244  0.69006556  0.82285074  0.7154544 ]
 [ 0.17053338  0.          0.81416039  1.00385442  0.95791547  0.78956039]
 [ 0.61237244  0.81416039  0.          0.          0.67082039  0.90886834]
 [ 0.69006556  1.00385442  0.          0.          0.74833148  1.09295263]
 [ 0.82285074  0.95791547  0.67082039  0.74833148  0.          0.4330127 ]
 [ 0.7154544   0.78956039  0.90886834  1.09295263  0.4330127   0.        ]]

расстояния correlation

[[  0.00000000e+00  -2.22044605e-16   2.25403331e-01   2.25403331e-01
    1.00000000e+00  1.00000000e+00]
 [ -2.22044605e-16   0.00000000e+00   2.25403331e-01   2.25403331e-01
    1.00000000e+00   1.00000000e+00]
 [  2.25403331e-01   2.25403331e-01   0.00000000e+00   2.22044605e-16
    3.67544468e-01   3.67544468e-01]
 [  2.25403331e-01   2.25403331e-01   2.22044605e-16   0.00000000e+00
    3.67544468e-01   3.67544468e-01]
 [  1.00000000e+00   1.00000000e+00   3.67544468e-01   3.67544468e-01
    0.00000000e+00   2.22044605e-16]
 [  1.00000000e+00   1.00000000e+00   3.67544468e-01   3.67544468e-01
    2.22044605e-16   0.00000000e+00]]

Вроде бы всё хорошо, но есть большое НО:
При расчёте реальных растров процесс вычисления дистанций превращается для хи-квадрат в многочасовой. Понятно что в методе с пользовательской функцией начинают крутится питоновские циклы вместо бинарных стандартных методов написанных на C или C++

Могу переписать функцию для формулы в C++ но тогда встаёт вопрос как потом её интегрировать в качетве пользовательской функции для scipy.spatial.distance.pdist и sdist?

п.с:
Пробовал реализовать пример последний с instant Но установленная версия instant-1.5.0 не имеет функции Instant() c методом create_extension и вообще по данному пакету информация очень фрагментарная … не исключаю что плохо искал.

terabayt

Апрель 3, 2015 22:15:52

а так правильно выдает результат?
время сократил в 8,5 раз

import numpy as np
from scipy.spatial.distance import pdist, squareform
a = np.array(
    [                                                                                          
        [1,1,1,2],                                                                             
        [3,3,3,5],                                                                             
        [1,2,3,4],                                                                             
        [2,4,6,8],
        [1,2,3,2],
        [3,4,5,4],
    ]
)
def X2_1(u, v):
    a = np.mean(u)*(u+v)
    b = np.mean(v)*(u+v)
    return np.sqrt(
        (
            (u - a/(np.mean(u)+np.mean(v)))**2/(
                a/(np.mean(u)+np.mean(v)))
        ).sum() + (
            (v - b/(np.mean(u)+np.mean(v)))**2/(
                b/(np.mean(u)+np.mean(v)))
        ).sum()
    )
b = pdist(
    a,
    X2_1
)
c = pdist(
    a,
    "correlation"
)
print a
print squareform(b)
print squareform(c)

kamisama

Апрель 3, 2015 22:43:48

Тут нет питоноциклов. Переписывание на си, если ты не байто-джедай ощутимой прибавки не даст.
В X2_1 много повторяющихся вычислений. К коду терабайта можно еще предвычисления средних добавить.

terabayt

Апрель 3, 2015 22:47:48

kamisama
К коду терабайта можно еще предвычисления средних добавить.

в алгоритм я не вникал, но идею дал как можно улучшить

oldbay

Апрель 4, 2015 02:27:46

terabayt
в алгоритм я не вникал, но идею дал как можно улучшить

Спасибо - идея дельная: буду использовать.

Но вот дальнейшие раскопки показали что в тормозах виновата как функция так и pdist. который принимая стороннию функцию: выполняется не уровне python-а и банально крутит его циклы(это в принципе указано и в документации к pdist-у)

Потому всё таки отимизировал скорость выполнения функции хи-квадрат через instant(переписал с C):

import numpy as np
from instant import inline_with_numpy
from scipy.spatial.distance import pdist, squareform
a = np.array(
    [
        [0.001,0.001,0.001,0.002,0.001,0.002],
        [0.003,0.003,0.003,0.005,0.001,0.002],
        [0.001,0.002,0.003,0.004,0.001,0.002],
        [0.002,0.004,0.006,0.008,0.001,0.002],
        [0.001,0.002,0.003,0.002,0.001,0.002],
        [0.003,0.004,0.005,0.004,0.001,0.002],
    ]
)
print a
hi2_code = """
double hi2_c(int u, double* uarr, int v, double* varr)
    {
        double usum = 0.0;
        for (int i=0; i<u; i++) {
            usum += uarr[i];
        }
        double umean = usum/u;
        double vsum = 0.0;
        for (int i=0; i<v; i++) {
            vsum += varr[i];
        }
        double vmean = vsum/v;
        double uhi = 0.0;
        for (int i=0; i<u; i++) {
            double Eu = umean * (uarr[i] + varr[i]) / (umean + vmean);
            uhi += pow((uarr[i] - Eu),2) / Eu;
        }
        double vhi = 0.0;
        for (int i=0; i<v; i++) {
            double Ev = vmean * (uarr[i] + varr[i]) / (umean + vmean);
            vhi += pow((varr[i] - Ev),2) / Ev;
        }
        return sqrt(uhi + vhi);
    }
"""
Hi2_c = inline_with_numpy(
    hi2_code,
    arrays = [
        ['u', 'uarr', 'in'],
        ['v', 'varr', 'in'],
    ]
)
b1 = pdist(
    a,
    Hi2_c
)
print squareform(b1)

Но этого тоже пока недостаточно - теперь тормозит pdist = придётся и его функционал в C перевести.