Вопрос по оптимизации скрипта [R]

[Alibek24]

Доброго времени суток!

В данный момент пишу диссертацию по геостатистике. Возникла проблема расчета вариограммы. Напомню что вариаграмма считается по формуле:
G=1/2N(h) * Sum[(x1-x2)^2]
где N - это число пар

У меня есть база данных состоящая из 5 колонок, 1ая это имя скважины, 2-3-4 это координаты x, y, z, и 5ая колонка это значение ллюбое (допустим кажущееся сопротивление). Вот как она выглядит:

Имя X Y Z KS
212 25 13 -5 0.5
.........................


И таких строк 28 тысяч.

Я написал скрипт на языке R:

pairs=function(well, h){
k=0
for(i in 1:nrow(well)-1){
for(j in i:nrow(well)){
if(well[i, "Имя"]!=well[j, "Имя"] $ well[i, "Z"]-well[j, "z"]<0.05 &
sqrt((well[i, "Z"]-well[i, "Z"])**2+(well[i, "Z"]-well[i, "Z"])^2+(well[i, "Z"]-well[i, "Z"])^2)<h){
k=k+1
}}}return(k)}

diff=function(well, h){
g=0
for(i in 1:nrow(well)-1){
for(j in i:nrow(well)){
if(well[i, "Имя"]!=well[j, "Имя"] $ well[i, "Z"]-well[j, "z"]<0.05 &
sqrt((well[i, "Z"]-well[i, "Z"])**2+(well[i, "Z"]-well[i, "Z"])^2+(well[i, "Z"]-well[i, "Z"])^2)<h){
g=g+(well[i, "KS"]- well[j, "KS"])^2
}}}return(g)}

c=seq(5,101, by=5)

for(i in c){
var=1/2*pairs(well, i) * diff(well, i)
print(var)}

Условия if там для того чтобы не брать одну и ту же скважину, примерно на одном уровне Z+-5cm, и расстояние было меньше h. Все работает но очень долго. Если кто сможет помочь улучшить алгоритм буду очень признателен. Спасибо!

[Pavia]

Язык R не знаю. Но общие рекомендации таковы.
1. Отказаться в цикле от именованных массивов в пользу индексных. К примеру завести временные массивы.
2. Вынести условия из цикла.
3. Степень заменить умножением.
4. Алгоритм сам по себе тяжёлый. Квадратичный. Тут нечем не поделаешь.

Цитата:

sqrt((well[i, "Z"]-well[i, "Z"])**2+(well[i, "Z"]-well[i, "Z"])^2+(well[i, "Z"]-well[i, "Z"])^2)<h)

Обычный трюк возвести в квадрат обе части неравенства. Тем самым избавляемся от корня. А h*h вычисляется до цикла. Но возможно это условие стоит выкинуть?
Или использовать разбить всё множество точек решёткой и искать пары только в квадратах это ускорит работу.

Надеюсь R работает с числами двойной точности(double float)

[iifat]

Как минимум — слить функции в один цикл, считать сразу количество и сумму. Раза в два быстрее получится.

[arseniiv]

Pavia в сообщении #751401 писал(а):

1. Отказаться в цикле от именованных массивов в пользу индексных. К примеру завести временные массивы.

Это может никак не действовать по многим причинам, так что сгодилось бы, возможно, лишь как последний пункт. («Именованных массивов»? Я думал, это называется нецелочисленные индексы.)

[Sphinx Pinastri]

Можно переписать на C.
А еще можно просто подождать. За выходные, наверно, сосчитает.
Чтобы было не скучно, можно выводить счетчик внешнего цикла.
Железо и электричество нынче дешевы, а труд дорог.

[Alibek24]

Pavia в сообщении #751401 писал(а):

Язык R не знаю. Но общие рекомендации таковы.
1. Отказаться в цикле от именованных массивов в пользу индексных. К примеру завести временные массивы.
2. Вынести условия из цикла.
3. Степень заменить умножением.
4. Алгоритм сам по себе тяжёлый. Квадратичный. Тут нечем не поделаешь.

Цитата:

sqrt((well[i, "Z"]-well[i, "Z"])**2+(well[i, "Z"]-well[i, "Z"])^2+(well[i, "Z"]-well[i, "Z"])^2)<h)

Обычный трюк возвести в квадрат обе части неравенства. Тем самым избавляемся от корня. А h*h вычисляется до цикла. Но возможно это условие стоит выкинуть?
Или использовать разбить всё множество точек решёткой и искать пары только в квадратах это ускорит работу.

Надеюсь R работает с числами двойной точности(double float)

Я сам в прораммировании не ас, но разве я не с индексами работаю? Условия из цикла не получится вынести, они должны же на каждой линии проверятся. А так по остальным советам все сделал, прирост не большой.

-- Сб авг 03, 2013 10:33:14 --

Sphinx Pinastri в сообщении #751458 писал(а):

Можно переписать на C.
А еще можно просто подождать. За выходные, наверно, сосчитает.
Чтобы было не скучно, можно выводить счетчик внешнего цикла.
Железо и электричество нынче дешевы, а труд дорог.

C не знаю, лучше знаю питон, но думаю питон будет по слабее с таким объемом данных. R как никак создан для дата анализа. А так на си можно было бы попробовать, но не подскажете как считать данные с формата .csv?

-- Сб авг 03, 2013 10:34:28 --

arseniiv в сообщении #751456 писал(а):

Pavia в сообщении #751401 писал(а):

1. Отказаться в цикле от именованных массивов в пользу индексных. К примеру завести временные массивы.

Это может никак не действовать по многим причинам, так что сгодилось бы, возможно, лишь как последний пункт. («Именованных массивов»? Я думал, это называется нецелочисленные индексы.)

Можно поподробнее про индексные массивы?

-- Сб авг 03, 2013 11:20:44 --

Еще один вопрос. На си можно было в цикле создавать списки с именами по порядку. Например если цикл ну одного до пяти то он создаст пять массивов с именами в1 в2 в3 в4 в5. Как этот прием назывался? Не знаю как запрос в гугле сделать

[Pavia]

arseniiv в сообщении #751456 писал(а):

«Именованных массивов»? Я думал, это называется нецелочисленные индексы.

Разные авторы по разному, называют. Также это называют ассоциативным массивом.
Наверно правильнее не целочисленные. Просто вчера не смог вспомнить как правильнее.

Alibek24 в сообщении #751460 писал(а):

Можно поподробнее про индексные массивы?

Я считаю что основные тормоза за обращения к массиву well[j, "Имя"]. Да теоретически оптимизатор мог справиться, но мой опыт подсказывает что обычно нет. Так что тут скорее всего основной тормоз прирост должен быть больным.

Массив в котором в качестве индекса используют целые числа работает быстрее нежели чем массив в котором в качестве индекса используют строки.
well[j, "Имя"] ; well[j, "X"]; well[j, "Y"]; well[j, "Z"] надо заменить на 4 массива.
name[j]; X[j]; Y[j]; Z[j]

Alibek24 в сообщении #751460 писал(а):

как считать данные с формата .csv?

Это простой формат. Он является текстовым. Столбцы разделены либо точкой либо точкой с запятой.
Строки разделены стандартным текстовым кодом новой строки - CRLF.

[Sphinx Pinastri]

Вы говорите, что программа работает, но в это верится с трудом.
Почему одно 'z' на нижнем регистре.
Зачем вычитать одинаковые величины. Ясно, что получится 0.
Зачем считать одно и то же много раз.

if(well[i, "Имя"] != well[j, "Имя"] $
      well[i, "Z"]-well[j, "z"] < 0.05 &
sqrt((well[i, "Z"]-well[i, "Z"])**2 + 
     (well[i, "Z"]-well[i, "Z"])^2  +
     (well[i, "Z"]-well[i, "Z"])^2) < h) { ... }

Наверно, этот фрагмент должен быть таким.

if(well[i, "Имя"] != well[j, "Имя"] $
      well[i, "Z"]-well[j, "Z"] < 0.05 &
sqrt((well[i, "Z"]-well[j, "Z"])**2 + 
     (well[i, "X"]-well[j, "X"])^2  +
     (well[i, "Y"]-well[j, "Y"])^2) < h) { ... }

Если можно, опубликуйте данные, а программистов здесь много.

[Alibek24]

Может быть где то и ошибся при написании кода, так как она была написана по памяти. Z это координата. для одной скважины таких координат может быть 500 и больше. Мне же надо сравнить данные на одном уровне для всех скважин, в радиусе h. Я проверил для более мелких циклов, работала

[Sphinx Pinastri]

Пожалуйста приведите программу, которая работает.
Поиск опечаток силами форума в том, что написано -- неэффективное использование
кредита доверия к автору.

Вы не ответили на вопросы в моем предыдущем сообщении, т.е. не видите бессмысленность
выполняемых операций даже когда на них указано явно.
Списывать сейчас не сезон, но другое объяснение придумать трудно.

Вам надо отсортировать результаты по Z и пробовать только точки с близкими Z.

[Alibek24]

Sphinx Pinastri в сообщении #751554 писал(а):

Пожалуйста приведите программу, которая работает.
Поиск опечаток силами форума в том, что написано -- неэффективное использование
кредита доверия к автору.

Вы не ответили на вопросы в моем предыдущем сообщении, т.е. не видите бессмысленность
выполняемых операций даже когда на них указано явно.
Списывать сейчас не сезон, но другое объяснение придумать трудно.

Вам надо отсортировать результаты по Z и пробовать только точки с близкими Z.

Ok завтра приведу код. После улучшении показанных здесь переписал код на питон, и по ощущениям работает быстрее. За 3 часа нашел вариаграммы на 5 и 10 метров.

Сортировка по Z даст прирост в одном, но в тоже время по именам скважин будет бардак, не? Так как именно название скважин должно быть разным.

Может быть ход моих мыслей неверен, и есть какой нибудь алгоритм типо Эратосфена, который даст прирост не в 2 или 3 раза, а на порядок?

[Sphinx Pinastri]

Значений Z много, 500 и больше, сортировка позволит сразу отбросить
неподходящие точки. Это должно ускорить программу на 2 порядка.
Имен тоже много. Значит, вклад от совпадающих имен несущественный.

Далее, современные процессоры имеют 6-8 потоков комманд. Значит можно
пускать 8 процессов сразу. Если своих компьютеров не хватает, то стоит
попроситься к кому-нибудь из приятелей.

[Alibek24]

Вот код на питоне

from numpy import recfromcsv
GT=recfromcsv('GT.csv', delimiter = ',')
def variogram(name, h):
gamma=0.0
k=0
dist=float(h*h)
for i in range(1, 27274):
line=name[i][0]
z=name[i][3]
x=name[i][1]
y=name[i][2]
gt=name[i][4]

for j in range(i, 27275):
if (line!=name[j][0]) and ((z-name[j][3])<0.1) and (((x-name[j][1])*(x-name[j][1])+(y-name[j][2])*(y-name[j][2])+(z-name[j][3])*(z-name[j][3]))<=dist):
gamma=gamma+(gt-name[j][4])*(gt-name[j][4])
k+=1
return gamma/(2*k)
a=[]
for i in range(5, 101, 5):
a.append(variogram(GT, i))
print a

[Sphinx Pinastri]

Программу поместите в блок code а то опять гадать надо что к чему относится.

За одно проверьте, не нужен ли модуль в

(z-name[j][3])<0.1

[Alibek24]

Sphinx Pinastri в сообщении #751657 писал(а):

Программу поместите в блок code а то опять гадать надо что к чему относится.

За одно проверьте, не нужен ли модуль в

(z-name[j][3])<0.1

про модуль очень хорошо подсказали, совсем забыл про это))

from numpy import recfromcsv
GT=recfromcsv('GT.csv', delimiter = ',')
def variogram(name, h):
     gamma=0.0
     k=0
     dist=float(h*h)
     for i in range(1, 27274):
          line=name[i][0]
          z=name[i][3]
          x=name[i][1]
          y=name[i][2]
          gt=name[i][4]

          for j in range(i, 27275):
               if (line!=name[j][0]) and ((z-name[j][3])<0.1) and (((x-name[j][1])*(x-name[j][1])+(y-name[j][2])*(y-name[j][2])+(z-name[j][3])*(z-name[j][3]))<=dist):
                   gamma=gamma+(gt-name[j][4])*(gt-name[j][4])
                   k+=1
     return gamma/(2*k)
a=[]
for i in range(5, 101, 5):
     a.append(variogram(GT, i))
     print a