LINUX.ORG.RU

Создание алгоритма методом «научного тыка»?


0

4

Задача из геометрической оптики (акустики). Есть треугольная лучевая трубка, образованная тремя лучами. Лучи движутся в среде с переменным показателем преломления, показатель меняется непрерывно. Местами лучевую трубку «выворачивает» (лучи перекрещиваются), нужно посчитать сколько раз ее «вывернуло».

Лучи считаются численно, методом Рунге-Кутты. Есть координаты лучей в начале шага и в конце шага, есть нормали лучей (единичные вектора, описывающие куда с-но лучи направлены). Казалось бы чего проще - считаем смешанные произведения чего нить с чем нить (вариантов много), если знак сменился значит трубку «вывернуло». Поскольку из общих соображений выражения должны быть симметричны (лучи равноправны), то следим за суммой знаков смешанных произведений (или за произведением знаков).

Но «дьявол кроется в деталях». Каждый из вариантов более-менее пристоен, но иногда дает ложные срабатывания, а иногда наоборот не срабатывает. Т.е. я могу посмотреть глазами на трубку и понять сколько раз ее вывернуло, но не могу научить машину определять это безошибочно. Число ошибок невелико, но они в итоге неприемлемо портят результат. Число трубок в боевой версии кода будет огромно, т.е. какой то ручной контроль будет исключен.

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

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

★★★★★

Ответ на: комментарий от anonymous

Фактически, нам надо определять моменты, когда меняется сонаправленность векторного произведени двух сторон и нормали.

Ну а нормаль ты откуда брать собрался?

ebantrop
()
Ответ на: комментарий от ebantrop

Я же выше написал. Вначале выбираем направление произвольным образом, а потом на каждом шаге исходим из того что за один шаг больше чем на 180 градусов треугольник не повернется.

anonymous
()
Ответ на: комментарий от anonymous

Я же выше написал. Вначале выбираем направление произвольным образом, а потом на каждом шаге исходим из того что за один шаг больше чем на 180 градусов треугольник не повернется.

Ну смотри, представь у тебя есть три отрезка:

A (-1,0,0)-(-1,0,1)

B (1,0,0)-(1,0,1)

и

C (0,-1,0)-(0,1,1)

То есть A и B лежат в плоскости xz, а C протыкает ее в точке (0,0,1/2). При этом может быть так что в начале шаги по отрезкам A и B будут малы, а C проткнет xz. Потом A B догонят, при этом твой критерий будет выполнен всегда, хотя ежу понятно что есть вывех.

ebantrop
()
Ответ на: комментарий от ebantrop

Еще раз, откуда ты берешь отрезки? В данных нету отрезков, там только точки. По три точки на шаг. Из этих точек составляется треугольник.

То есть A и B лежат в плоскости xz, а C протыкает ее в точке (0,0,1/2). При этом может быть так что в начале шаги по отрезкам A и B будут малы, а C проткнет xz. Потом A B догонят, при этом твой критерий будет выполнен всегда, хотя ежу понятно что есть вывех.

Да не важно кто кого и как догоняет. Нам нужно только два шага. На первом шаге у треугольника будет левая ориентация, на втором - правая. Это значит что между шагами произошло выворачивание. Если между шагами ориентация не меняется, то значит выворачивания не было. Ну или было дважды, но это уже не отследить.

anonymous
()
Ответ на: комментарий от ebantrop

А, я понял, у тебя первый треугольник - это (-1,0,0)-(1,0,0)-(0,-1,0), а второй - (-1,0,1)-(1,0,1)-(0,1,1)?

Соответственно, за шаг треугольник отражается. При этом нормаль меняется на 180 градусов - то есть произошло выворачивание, все ок.

anonymous
()
Ответ на: комментарий от ebantrop

Ха, прогнал сейчас с данными, которые давал ОП - там полный пиздец. На первом «вывороте» (шаги 163-164) треугольник практически вырождается в отрезок, и на глаз то кажется, что там трубка «вывернулась», но это на глаз и именно _кажется_ - на деле то хрен его знает действительно она там вывернулась, или просто треугольник немного провернулся, картинка-то будет такой же, а провернуться ему надо было совсем немного. То есть вопрос в том, прошла ли точка по прямой или отклонилась немного и ушла в дугу, но это никак не выяснить. А могло быть и так и так. И выворот или был или не был - и не важно, чего там на глаз «показалось».

Если там все случаи такие то на этих данных каши не сваришь.

anonymous
()
Ответ на: комментарий от anonymous

http://pastebin.com/3pnM8A6b

Твое поделие показывает 3 вывиха, а там он один.

UPD: рисует тоже какойто бред, должно быть так:

http://picpaste.com/untitled-T4wEx4le.png

ebantrop
()
Последнее исправление: ebantrop (всего исправлений: 1)
Ответ на: комментарий от anonymous

мне для визуализации тех данных помогло следующее:

$ awk '
BEGIN{
  i=0; X=0;Y=1;Z=2;
  a[0]=1;a[1]=2;a[ 2]=3;a[ 3]=1;
  a[4]=4;a[5]=5;a[ 6]=2;a[ 7]=5;
  a[8]=6;a[9]=3;a[10]=6;a[11]=4;
}
/^#/{next}
{
  if(NF!=18) next;
  if(i!=0) {
    s[4,X]=$10; s[4,Y]=$11; s[4,Z]=$12;
    s[5,X]=$13; s[5,Y]=$14; s[5,Z]=$15;
    s[6,X]=$16; s[6,Y]=$17; s[6,Z]=$18;
    for(j=0;j<12;j++)
      print s[a[j],X],s[a[j],Y],s[a[j],Z] >sprintf("%i.dat",i)
  }
  s[1,X]=$10; s[1,Y]=$11; s[1,Z]=$12;
  s[2,X]=$13; s[2,Y]=$14; s[2,Z]=$15;
  s[3,X]=$16; s[3,Y]=$17; s[3,Z]=$18;
  i++
}
' "tube.dat"
$ gnuplot
gnuplot> splot "1.dat" using 1:2:3 with line, "2.dat" using 1:2:3 with line, "3.dat" using 1:2:3 with line, "4.dat" using 1:2:3 with line, "5.dat" using 1:2:3 with line, "6.dat" using 1:2:3 with line, "7.dat" using 1:2:3 with line, "8.dat" using 1:2:3 with line

anonymous
()
Ответ на: комментарий от anonymous

Но прокол же не является выворотом, выворот происходит только тогда, когда луч пересекает сторону _текущего_ треугольника, а это далеко не любой прокол. То есть в этом случае в какой-то точке треугольник превращается в прямую. Понятно, что эта точка обычно будет где-то между шагов, но тут ее очевидно нет, нигде между шагов треугольник в прямую не превращался.

А если любой прокол считать, то у тебя выворотом будет банальный разворот треугольника «на месте» (ну или просто с низкой скоростью. Хотя тут надо, конечно, чтобы ОП пришел и дал точное определение, что считать выворотом, а что нет.)

anonymous
()
Ответ на: комментарий от anonymous

Хотя тут надо, конечно, чтобы ОП пришел и дал точное определение, что считать выворотом, а что нет.)

Походу ТС забил на топик. Как я понял выворот это когда один из лучей (луч в смысле поста) пересекает поверхность, образованную двумя другими. Поэтому имеет смысл говорить об отрезках, которые описывают эти три луча. То есть разворотом треугольника тут не обойдешься, либо надо строить все возможные треугольники, что в лоб есть O(n^3).

ebantrop
()
Ответ на: комментарий от ebantrop

а чем та поделка не устраивает? для какого случая она не работает? или там сложность тоже O(n^3)?

anonymous
()
Ответ на: комментарий от ebantrop

Я не забил, я был в дороге. В наших РЖД-шных поездах интернетов нету;-(

Я сейчас смотрю отдельные призмы, соответствующие «вывертам». Чудовищно муторное занятие - текущая версия алгоритма почти всегда дает правильный ответ (в т.ч. и в тех случаях где я сам сходу сказать не могу что к чему), но иногда лажает.

Вот чекер http://a-iv.ru/trash/check-tubes.py http://a-iv.ru/trash/check-tubes.gp - надо запускать питоний файл, указывая интересующий трубки в качестве аргумента. si - это знаки смешанных произведений, I - индекс (число вывертов).

Трубки вот тут http://a-iv.ru/trash/tubes-G1.tgz какие из них плохие сказать все еще не могу;-( в имени файла I указывает на число вывертов определенное текущей версией алгоритма, остальное неважно.

Насчет смешанных произведений анонимус прав, только там много вариантов смешанных произведений, иногда реагирует один вариант, иногда другой.

psv1967

Анализ принципиальных компонент. Идея в том, что траектория будет развернута этим методом вдоль основных осей.

А что в данном случае является основными осями? Навскидку это кривая проходящая через центры сечений, и еще какой то угол вращения вокруг это кривой?

AIv ★★★★★
() автор топика
Ответ на: комментарий от AIv

В экстремальном случае как я показал будет срабатывать?

http://pastebin.com/3pnM8A6b

Как мне кажется если работать с треугольниками или призмами, выверты надо смотреть еще +/- несколько точек вокруг по каждому лучу. Поскольку выверт не обязательно произойдет на одном шаге. Если на одном, то таки да, анон абсолютно прав.

ebantrop
()
Ответ на: комментарий от ebantrop

Походу ТС забил на топик. Как я понял выворот это когда один из лучей (луч в смысле поста) пересекает поверхность, образованную двумя другими.

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

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

anonymous
()
Ответ на: комментарий от ebantrop

надо смотреть еще +/- несколько точек вокруг по каждому лучу. Поскольку выверт не обязательно произойдет на одном шаге. Если на одном, то таки да, анон абсолютно прав.

Не несколько, а вообще все. Если убрать условие «одношаговости», то прокол траектории может происходить в любой момент, после того, как эта траектория появилась. Хоть через сто шагов, хоть через миллион.

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

Опять же, если треугольник двигался вниз, а потом уменьшился и начал двигаться вверх (оказавшись внутри старой траектории) - это считается за выворачивание? А если на каком-то шаге он потом опять расширился (сделав прокол?

anonymous
()
Ответ на: комментарий от ebantrop

В экстремальном случае как я показал будет срабатывать?

ЯННП (что именно Вы показали). Это пример трубки, алгоритм или что?

AIv ★★★★★
() автор топика
Ответ на: комментарий от anonymous

Не несколько, а вообще все. Если убрать условие «одношаговости», то прокол траектории может происходить в любой момент, после того, как эта траектория появилась

Ну да, тогда сложность О(n^3). Луча то три.

ebantrop
()
Ответ на: комментарий от ebantrop

Такие случаи бывают (обычно запаздывание на один шаг), и текущая версия алгоритма их обрабатывает неправильно.

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

AIv ★★★★★
() автор топика
Ответ на: комментарий от AIv

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

наверное не должно быть пересечений этих трех «лучей» в проекциях ориентированных по «макроизгибам»?

если добавить другие характеристики (локальные для конкретного сечения трубки), то надо строить окно из сечения вперед, текущего и сечения назад и смотреть или просто PCA картину на предмет кластеров естественных. или учить какой нибудь рандом форест по образцу.

psv1967 ★★★★★
()
Ответ на: комментарий от AIv

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

код для файла «pos0x17-I4-rank4-13.dat» формирует результат вида:

104 2.1 1
106 2.14 2
217 4.36 3
count: 3
последний столбец - вид пересечения ребер:

  • 1 - одно ребро пересекает плоскость образованную двумя другими ребрами,
  • 2 - два ребра скрещиваются в одной точке
  • 3 - все 3 ребра скрещиваются в одной точке
anonymous
()
Ответ на: комментарий от anonymous

«просмотр фрагментов» первым параметром принимает номер точки начала просмотра и от неё рисует 8 следующих призм.

anonymous
()
Ответ на: комментарий от AIv

вот собственно три соседних выверта

l <- system2(«ls», "./tubes", stdout = T)

schift.3.data <- function(data) {
  cbind(data[1:(nrow(data)-2),],
        data[2:(nrow(data)-1),],
        data[3:nrow(data),])
}


data <- do.call(rbind, lapply(1:length(l), function(i) schift.3.data(read.table(paste0("./tubes/",l))[,2:11]) )) 

plot(prcomp(data)$x[,1:2])

получается вся изменчивость вывертов практически описывается плоскостью.

на плоскости симметричная картинка, многие ситуации накладываются друг на друга. есть редкие ситуации.

нужны сведения какие из файлов описывают ошибки. или на какой строке файла наступает событие (в принципе подойдет по методу исключения и указание какой файл совсем без ошибок)

psv1967 ★★★★★
()
Последнее исправление: psv1967 (всего исправлений: 1)
Ответ на: комментарий от anonymous

у него «граф переходов» в «вывертах» выглядит вот так


which(l==«pos0x21-I3-rank4-6.dat»)

[1] 132

result<-prcomp(data)

plot(predict(result, schift.3.data(read.table(paste0("./tubes/",l[132]))[,2:11]))[,1:2], type=«b»)

а pos0x21-I4-rank4-11.dat имеет пересечения? или pos0x22-I3-rank4-3.dat ?

psv1967 ★★★★★
()
Последнее исправление: psv1967 (всего исправлений: 1)
Ответ на: комментарий от psv1967

мой код на файл «pos0x21-I3-rank4-6.dat» дает следующий результат:

147 2.96 2
148 2.98 1
149 3 3
302 6.06 1
count: 4
неверно трактуется призмы 147-149, в них находится 2,1,3 переходы. на глаз там видится в 147-148 призмах двойной выверт, и 149 призма одиночный выверт. кроме них ещё в 302 призме выверт, вполне обычный. если я правильно трактую слово «выверт».

на файл «pos0x21-I4-rank4-11.dat» результат:

149 3 1
155 3.12 1
290 5.82 1
count: 3
код не нашел 4-тый «выверт»

на файл «pos0x22-I3-rank4-3.dat» результат:

159 3.2 2
160 3.22 1
161 3.24 1
count: 3
сомнения на счет 2-ного «выверта» 159-той призмы, в остальных кажется есть «выверты»

anonymous
()
Ответ на: комментарий от anonymous

первый столбец - номер второго основания, второй столбец - это время второго основания призмы (время текущего треугольника), третий столбец - вид «пересечения»

anonymous
()
Ответ на: комментарий от psv1967

ЯННП. Какой ЯП то хоть?;-)

нужны сведения какие из файлов описывают ошибки

Все *-I3-* гарантированно лажовые, I=3 в той постановке быть не могло, т.е. текущий алгоритм сработал лишний раз.

Мне как то удалось эту хрень классифицировать, но пока не работает;-(

AIv ★★★★★
() автор топика
Ответ на: комментарий от anonymous

вот новый код,

Круто! Жалко я баша толком не знаю, но постараюсь раскурить, спасибо.

AIv ★★★★★
() автор топика
Ответ на: комментарий от AIv

там awk главную задачу выполняет:

...
awk -v DST=$DST -v DBG=$DBG
смотри вот этот код, остальное графики
' $SRC
...

anonymous
()
Ответ на: комментарий от AIv

R на сколько я понимаю, только у меня к сожалению не получается ни один из его примеров повторить, не хватает знаний R. (если что пробовал на RStudio v0.97.551, R version 2.15.2 (2012-10-26))

ошибки следующие, первая:

> l <- system2("ls", "./tubes", stdout = T)
> schift.3.data <- function(data) {cbind(data[1:(nrow(data)-2),], data[2:(nrow(data)-1),], data[3:nrow(data),])}
> data <- do.call(rbind, lapply(1:length(l), function(i) schift.3.data(read.table(paste0("./tubes/",l))[,2:11]) ))
Ошибка в file(file, "rt") : неправильный аргумент 'description'
вторая:
> l <- system2("ls", "./tubes", stdout = T)
> schift.3.data <- function(data) {cbind(data[1:(nrow(data)-2),], data[2:(nrow(data)-1),], data[3:nrow(data),])}
> result<-prcomp(data)
Ошибка в as.vector(x, mode) : 
  cannot coerce type 'closure' to vector of type 'any'

anonymous
()
Ответ на: комментарий от anonymous

> plot(result$x[,1:2])
> lines(unique(predict(result, schift.3.data(read.table(paste0("./tubes/",l[132]))[,2:11]))[,1:2]), col="red")
> which(l=="pos0x21-I4-rank4-11.dat")
[1] 134
> lines(unique(predict(result, schift.3.data(read.table(paste0("./tubes/",l[134]))[,2:11]))[,1:2]), col="blue")
> which(l=="pos0x22-I3-rank4-3.dat")
[1] 151
> lines(unique(predict(result, schift.3.data(read.table(paste0("./tubes/",l[151]))[,2:11]))[,1:2]), col="green")


довольно характерный рисунок пишет. причем симметричный, сейчас попробую все траектории «вывертов» вывести, может они еще чем то отличаются...

psv1967 ★★★★★
()
Последнее исправление: psv1967 (всего исправлений: 1)
Ответ на: комментарий от anonymous

попробуйте просто в директории где распаковали архив tubes (он распаковывается в ./tubes) запустить R

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

psv1967 ★★★★★
()
Ответ на: комментарий от psv1967

да так и делаю, не знаю. и если просто путь указывать то строит все нормально. правда не понимаю что вижу. то ли это, что должен видеть.

anonymous
()
Ответ на: комментарий от AIv

функция solve: первые 3 параметра - первый треугольник в пространстве, вторые 3 параметра - второй треугольник. каждый из треугольников образует плоскость (S1 и S2). 7, 8 параметры, проверяемые точки (p1 и p2). остальные параметры для отладки.

проверка заключается в определении расположения точек p1 и p2 относительно плоскостей S1 и S2. для этого используется функция getLength, определение расстояния от точки до поверхности с учетом знака («сверху» или «снизу») относительно поверхности.

пронумеруем точки 1,2,3 - точки первого треугольника, 4,5,6 - точки второго треугольника. если выбрать произвольную грань призмы 1-4, то точки 5 и 6 (p1 и p2) должны быть расположены с одной стороны от двух плоскостей: первой S1 - проведенной через основание (1,2,3) и S2 - через новый треугольник (4,2,3). если произведение расстояний для каждой из точек имеет отрицательное значение, то выбранная грань «вероятно» проходит через «плоскость» образованную двумя оставшимися гранями. анализ по всем граням начиная со старого и с нового треугольника дает какое-то понимание проблемы, правда иногда не правильное.

и код показывающий идею:

...
function solve(s1,s2,s3,s4,s5,s6, p1,p2, l,i) { 
  ...
  l0=getLength(S1,p1); l1=getLength(S2,p1);
  l2=getLength(S1,p2); l3=getLength(S2,p2);
  l[i+0]=l0*l1;
  l[i+1]=l2*l3;
  if(l[i+0]<-1E-20) cnt++;
  if(l[i+1]<-1E-20) cnt++;
  if(cnt==2) return 1;
  ...
}
...
  res1=solve(s1,s2,s3,s4,s2,s3, s5,s6, l,0);
...
остальной код - на различные тонкие случаи.

anonymous
()
Ответ на: комментарий от psv1967

кажется нашел ошибку в этой функции:

data <- do.call(rbind, lapply(1:length(l), function(i) schift.3.data(read.table(paste0("./tubes/",l))[,2:11]) ))

нужно было с написать l[i], а я просто повторил. (лучше код писать в тегах [code], а то [i] съедается)

data <- do.call(rbind, lapply(1:length(l), function(i)
schift.3.data(read.table(paste0("./tubes/",l[i]))[,2:11]) )) 

anonymous
()
Ответ на: комментарий от anonymous

да, ещё, в том коде не проверяется случай построения плоскости через точки расположенные на одной прямой (или близкие к одной прямой), может оттуда ошибки оставшиеся берутся.

anonymous
()
Ответ на: комментарий от anonymous

вот и следующая ошибка:

> l <- system2("ls", "./tubes", stdout = T)
> schift.3.data <- function(data) {cbind(data[1:(nrow(data)-2),], data[2:(nrow(data)-1),], data[3:nrow(data),])}
> data <- do.call(rbind, lapply(1:length(l), function(i) schift.3.data(read.table(paste0("./tubes/",l[i]))[,2:11]) ))
> result<-prcomp(data)
> plot(result$x[,1:2])
Ошибка в plot.new() : края рисунка слишком велики
решилась совершенно глупо, развертыванием RStudio на весь экран

anonymous
()
Ответ на: комментарий от psv1967

построить получилось, но я не понимаю, что я вижу на этом рисунке. все равно не хватает знаний о смысле функций: prcomp и predict.

anonymous
()
Ответ на: комментарий от anonymous

там 12 основных состояний получается (для этого надо смотреть как накапливаются точки в пространстве первых двух компонент, остальные компоненты по графику нагрузок plot(prcomp()) похоже не существенны).

между этими основными состояниями небольшая кучка переходных.

процесс между этими 12 состояниями «переключается», очевидно некоторые из сочетаний переходов «запретны».

надо обрабатывать последовательность событий получается :)

тут или марковские цепочки (там есть возможность прямо первые две компоненты в плоскости которых всё происходит обработать), или из кернеллаба вариант представляющий последовательность (только надо четко свести к состояниям) в виде сочетания подпоследовательностей.

PS эти 12 состояний симметричны по первой выделенной компоненте 6 на 6.

psv1967 ★★★★★
()
Последнее исправление: psv1967 (всего исправлений: 1)
Ответ на: комментарий от AIv

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

При проколе ничего не меняется-то. Он вообще не является топологической особенностью :(

anonymous
()
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.