Category: наука

Category was added automatically. Read all entries about "наука".

slow loris

Big Data, чтоб их ... (4)

Наконец-то стряхнул пыль с компьютерной модели сближения, добавил в неё код, чтобы мы могли определить интересующие нас точки, и выписать ковариационные матрицы Qk, найденные там, чтобы потом преобразовать в весовые матрицы по формуле



где



Кстати, если сложить все Ak, мы должны получить единичную матрицу, что элементарно доказывается - матрица R вынесется за скобки, и тогда получится R-1R=E. Так что несмещённость результатов обеспечена.

Но когда я попробовал всё это посчитать, поначалу вышел полный бред - не выполнялось это условие, и всё тут, какая там единица на диагонали - нескольких порядков не хватало! Всё дело было в том, что я с какого-то перепуга подумал, что матрицы Ak получаются симметричными! Ну, как бы все матрицы Qk симметричны (для ковариационных это обязательно), и обратные к ним симметричны, и потом их сумма симметрична, и обратная к ней, т.е матрица R, также симметрична. А вот произведение двух симметричных - УЖЕ НЕТ.

Простейший пример:


А когда-то я об этом точно знал и написал процедуру перемножения симметричных матриц, дающую на выходе матрицу "общего вида". А сейчас у меня лежали эти матрицы в массиве 10х21 (десять симметричных матриц 6х6, для которых нужно 21 элемент), я попытался одну за другой умножать "на месте", а результат умножения (36 элементов) затирал соседнюю матрицу, вот всё и посыпалось как костяшки домино...

Кстати, из-за этого и выражение, как данные усреднять, я в прошлый раз неправильное привёл. Правильное вот:



(отличие в транспонировании матриц! Изначально они транспонированные и получались, но я подумал - они же симметричные, можно знак убрать, так оно красивее...)

Когда это дело подправил, получил довольно интересные результаты! Я боялся, что вообще получу практически единичные матрицы, или, по крайней мере диагональные (см. Big Data, чтоб их... (3)), но не тут-то было! Но чтобы понять, что это мне даёт, нужно ещё чуть подумать - вывести итоговую ковариационную матрицу результатов, если мы воспользуемся этими весами!

Collapse )

Приходим к удивительному результату: матрица R, которую мы посчитали ранее - это и есть результирующая ковариационная матрица, выражающая, как шумит наше усреднённое "по науке" значение.

Можно ли было это понять методом "пристального взгляда", обязательно ли было эти выкладки проводить - пока не знаю.

Завтра, наконец, узнаем, насколько будут отличаться результаты разных усреднений в условиях, "максимально приближенных к боевым"...
doomguy

Ликбез по кватернионам, Мортари-Маркли берут реванш!

Часть 1 - история вопроса
Часть 2 - основные операции
Часть 3 - запись вращения через кватернионы
Часть 4 - кватернионы и спиноры; порядок перемножения
Часть 5 - практическая реализация поворота
Часть 5 1/2 - введение метрики, "расстояния" между поворотами
Часть 6 - поворот по кратчайшему пути
Часть 6 1/4 - кратчайший поворот в общем случае
Часть 6 2/4 - поворот, совмещающий два направления
Часть 6 3/4 - кватернион из синуса и косинуса угла
Часть 7 - интегрирование угловых скоростей, углы Эйлера-Крылова
Часть 8 - интегрирование угловых скоростей, матрицы поворота
Часть 8 1/2 - ортонормирование матрицы и уравнения Пуассона
Часть 9 - интегрирование угловых скоростей с помощью кватернионов
Часть 10 - интегрирование угловых скоростей, методы 2-го порядка
Часть 10 1/2 - интегрирование с поддержанием нормы
Часть 11 - интегрирование угловых скоростей, методы высших порядков (в разработке)
Часть 12 - навигационная задача
Часть 13 - Дэвенпорт берёт след!
Часть 14 - линейный метод Мортари-Маркли
Часть 15 - среднее от двух кватернионов
Часть 15 1/2 - проверка и усреднение кватернионов
Часть 16 - разложение кватерниона на повороты
Часть 17 - лидарная задача
Задачка к части 17
Дэвенпорт VS Мортари-Маркли
Мортари-Маркли берут реванш!

В пятницу мы попробовали найти поворот из зашумлённых векторов, используя два разных метода - более классический Дэвенпорта и относительно новый (2007 года) Мортари-Маркли. Первый, будучи куда более сложным в реализации (особенно на очень простом "космическом" железе, умеющем складывать и умножать в целых числах, и ничего более), вышел на 11% точнее.

Но это была всего лишь одна выборка случайного шума, по одной выборке делать результаты - совершенно недопустимо! Поэтому решил всё-таки собрать статистику чуточку побольше. Но тут уже нельзя опираться на Wolfram Alpha, который найдёт все собственные значения и собственные вектора - надо научиться делать это самостоятельно, "малой кровью". К счастью, для данной конкретной задачи есть несколько интересных наблюдений, существенно упрощающих нам жизнь!
Collapse )

Что ж, метод Мортари-Маркли частично оправдан: в нашем примере он в среднем даёт ошибку на 6,5% больше при точности измерений в 1 мм, и на 3% больше при точности измерений в 0,1 мм. Хотя статистики до сих пор маловато: мы рассмотрели только один поворот, на 90°, и ровно одну конкретную мишень.

Что интересно, и метод Дэвенпорта мы частично оправдали: оказалось, что и по нему можно находить кватернион "в пол-пинка", поскольку нужное нам собственное значение мы можем узнать по простой формуле через входные данные, с большой точностью. Собственно, так и начинается алгоритм QUEST (QUaternion ESTimator).

Наверное, пора немножко отвлечься от этой математики (есть не один, а целых два алгоритма, которые меня вполне устраивают!) и вернуться к макету. Скоро ехать к заказчику, надо как-то немножко лоску навести на всё это безобразие. Не знаю, правда, за что браться.
QuatCore

Алгоритм Ремеза в экселе

Вот и до него руки дошли, причина станет ясна в следующем посте. Изучать чужие библиотеки было лениво (в том же BOOSTе сам чёрт ногу сломит), писать самому - тоже, а вот "в полуавтоматическом режиме" в Excel - вышло в пол-пинка.

Для примера возьмём приближение


при условии x2+y2=1.

Чтобы подобрать коэффициенты, делаем "замену" x=cos(z), y=sin(z), и по сути:


Используя метод наименьших квадратов и "вручную" подбирая диапазон оптимизации, мы смогли подобрать α, β чтобы арктангенс считался с точностью не хуже 2,5 угловых минут. Позже, сделав интерполяцию по узлам Чебышёва, смог получить результат чуть лучше: 2,36 угловых минут. Интересно, можно ли улучшить этот результат, и насколько сильно.

Нам понадобится график ошибки, который приводили в прошлый раз:
ChebError.png

Collapse )

Как я и ожидал, удалось добиться улучшения, но не очень большого: с 2,36 угловых минут до 2,28 угловых минут (ещё минус 3,3%) и это, похоже, предел для такой формулы.

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

Таблица коэффициентов для atan1

Пора уже добить этот несчастный арктангенс, что-то слишком долго его ковыряю. Мы упоминали, что в оперативной памяти будет лежать таблица в 15 слов: 5 углов и 5 комплексных чисел, им соответствующих. Вопрос был, в каком порядке всё это дело расположить, ну и сами константы "придумать". Пожалуй, вот так:

;таблица поворотов для процедуры atan1
AtanTable:
;завершающая итерация
;пока длина вектора составляет 1,0000032 от исходной (в теории, т.е если коэффициенты заданы ровно эти, но вычисления с абсолютной точностью)
Atan4a	Int16 1608		;pi/32 в формате Q2.14
;если cos(pi/32) и sin(pi/32), это было бы 32610 и 3212.
;если cos(1608/16384) и sin(1608/16384), то 32610 и 3211
;но мы ещё должны поделить на 2, т.к формат 2.14, а вектор был в 1.15
;(и это деление проведём, заменив ADD на DIV2A в конце процедуры)
;а ещё умножить на k=1,00120570556/1,0000032
;в итоге возьмём cos(1608/16384)*k и sin(1608/16384)*k
Atan4x	Int16	32650
Atan4y	Int16	3215


;к этому моменту длина вектора составляет 1,000013 от исходной
Atan3a	Int16 3217		;pi/16 в формате Q2.14
Atan3x	Int16 32138		;exp(i*pi/16) = 0,9808 + 0,1951i
Atan3y	Int16 6393

;к этому моменту длина вектора составляет 1,00000107 от исходной (в теории, это sqrt(23170^2+23171^2)/32768)
Atan2a	Int16 6434		;pi/8 в формате Q2.14
Atan2x	Int16 30274		;exp(i*pi/8) = 0,9239 + 0,3827i
Atan2y	Int16 12540

;к этому моменту длина вектора не менялась
Atan1a	Int16 12868		;pi/4, выраженные в формате Q2.14 (в который влезают значения от -4 до 4)
Atan1x	Int16 23170		;exp(i*pi/4) = 0,7071 + 0,7071i
Atan1y	Int16 23171		;а точнее, exp(i * 12868 / 16384), чтобы осуществить именно тот поворот, что указан в Atan1a

Atan0a	Int16 0
Atan0x	Int16 0		;первый поворот, на +- 90 градусов (0+1i)
Atan0y	Int16 -32768	;взяли знак "-" потому что первый раз мы берём обратный знак (из нуля вычитаем "x")
;длина вектора остаётся прежней


Как всегда, здесь всё задом наперёд!
Collapse )

Наверное, в очередной раз заморочался куда сильнее, чем следовало бы. Очень странная форма перфекционизма - придумать приближённое решение, но потом попытаться выжать из него все-все соки.
QuatCore

atan(y/x) на QuatCore

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

Но всё же, чтобы тупо выполнить условия ТЗ и выдать УГЛЫ, которые у нас просят (курс, тангаж, крен), нужно под самый конец преобразовать кватернион в эти углы.

Я так прикинул, для этого нужна ровно одна обратная тригонометрическая функция, которую условно назвал atan(y/x), не путать с atan2(y,x)!

Последняя выдаёт угол между осью OX и направлением на точку (x,y), с учётом квадранта. Как правило, ответ лежит в диапазоне от -pi до +pi, т.е "полный оборот". Но у нас сплошь используются половинные углы в кватернионах, поэтому результат надо ещё умножить на 2, и как раз такая "разборчивость к квадранту" оказывается излишней! Как раз-таки самый простой atan(y/x), выдающий значения от -pi/2 до pi/2, оказывается предпочтительнее, только вот "в явном виде" выполнять деление - наверное не лучшее решение. Когда мы приближаемся к плюс-минус pi/2, x стремится к нулю, аргумент арктангенса стремится к бесконечности и "достигает" её на краях.

Уж лучше этой функции дать два аргумента: x,y, т.е реализовать некий atan(y,x), в котором на самом деле мы пойдём "окольным путём".

Я хочу реализацию в целых 16-битных числах. Каким считать масштаб y,x - не так уж важно, главное, он один и тот же. А результат будет иметь формат Q3.13, т.е 3 бита перед запятой, 13 после запятой, со знаком, т.е диапазон от -4 до 3,9999 радиан. Цена младшего разряда около 25 угловых секунд, значит, точность порядка 12 угловых секунд, что соответствует точности задания угла кватернионом с 16-битными компонентами.

Хочу точность как минимум 0,1°, а лучше точнее. При этом максимально компактный код, т.к память "поджимает". А вот время выполнения меня не сильно заботит: нам отводится 200 мс на обработку измерений, а пока мы укладываемся в 200..300 мкс.

Collapse )

Похоже, нащупали неплохой способ, осталось его реализовать на ассемблере.

PS. Хотя забавно, как этот способ допускает максимальную ошибку при вычислении atan(0). У него получается не ноль, а 8 угловых секунд, которые, надеюсь, округлятся-таки до нуля при выдаче 16-битного результата. Хотя казалось бы, уж в нуле ошибиться никак нельзя, тут практически любой сходу даст правильный ответ!
QuatCore

Аффинный алгоритм+информационный обмен, минус два бага

Посидел глубоко в отладке, вроде разобрался, что к чему. Всё неправильное поведение свелось к двум багам, один в модуле RTC (Real Time Clock), второй - в программе, я в неё влез по дурости, показалось, что можно получше сделать, в итоге выкинув одну строчку, запорол одновременно и нормировку кватерниона, и компоненты Y,Z вектора параллельного переноса. Это ещё умудриться надо было...

Collapse )

Завтра проверю, что всё правильно - и придётся в очередной раз резко сменить род деятельности. Сегодня мне в РКК полтора часа мозги компостировали насчёт юстировки прибора. У нас для этого в теории существуют КПАшники (те, кто занимаются Контрольно-Проверочной Аппаратурой), целый отдел, но они самоустранились, сославшись на Стандарт Предприятия (СТП), по которому они делают только ту часть, что отправляется на входной контроль, это абсолютный минимум, способный лишь показать "ничего не повредили пока везли", а вот рабочее место, на котором проверяется львиная доля требований, а также всё, что касается работы заказчиков с прибором, включая измерение положения прибора на аппарате - должен придумать разработчик прибора, то есть я. Да, идейки были, но пора их оформить...
Doc

Самодельный "холодный потолок" - дополнение

О том, насколько полипропиленовые трубы по теплоотдаче проигрывают металлическим, сколько сейчас мощности уходит на прокачку воды, и какое гидравлическое сопротивление этой системы должно быть "в теории", т.е сильно ли мне навредили эти сужения?

Collapse )

В общем, ясно, что ничего не ясно :)

Если в эти выходные снова будет жарко, возможно, продолжу эксперимент.
QuatCore

Алг. ближ. дист, вычисление масштаба

Вот знал, что имеющийся код, который предполагался универсальным и для ближней дистанции, и для дальней, всё же придётся причесать. Вон, с креном обнаружилось переполнение при работе с полметра (минимально возможной дистанции). Это мы поправили. Следующая на очереди - процедура "вычленения" масштаба из матрицы масштаба и ракурсов.

Старый код:
	;Состояние регистров к этому моменту:
	;X=AfTransf,
	;Y=QuatY,
	;Z=Matrix,
	;i=j=k=Inv=0,
	;C=Acc=Txx
	FindScale	proc
	;1. Вычисляем величины |X-Z| (т.е |Txx - Tyy| ) и |2Y| (например, |Txy+Tyx|, хотя можно 2|Txy| или 2 |Tyx| )
	;важно, что в данный момент Inv=0
				[SP]		0	;манхэттенская метрика, сумма модулей (метрика городских кварталов)
				[SP+1]		0	;чебышевская, максимум от модулей 
				;текущее значение будем хранить в регистре C - дешево и сердито!
				i		1
				j		3					
		@@loop:		Acc		[X+i]	;+Tyx (i=1), -Tyy (i=0)
				PM		[X+i^j]	;Txy (i=1), Txx (i=0)   ;поменяли местами чтобы проверить модуль
				Inv		1		;уже можно! Так что на следующей итерации будет "-"
				[Y+i]		0		;очистим две компоненты кватерниона
				ABS		Acc		;|Txy+Tyx| (i=1), |Txx-Tyy| (i=0)
				C		Acc		;сохранили это значение
				ADD		[SP]	
				[SP]		Acc	;обновили манхэттенскую
				Acc		C
				SUB		[SP+1]
				JL		@@skipCheb ;лучше подошел бы JGE и поменять местами C и [SP+1]. Разница при нуле-на 1 больше операций.
				[SP+1]		C	;обновили Чебышевскую					
		@@skipCheb:	iLOOP		@@loop
	;теперь в [SP] лежит манхэттенская метрика, а в [SP+1] - чебышевская
	;осталось посчитать выражение целиком - сумму Txx+Tyy с весом 1/2, а эти - с весами 0,332 и 0,168
	;или сначала всё сложить - а потом поделить пополам...
	;можно, если ввести команду UDIV2 - беззнаковое деление
	;по сути, SHR1 
	
	;кстати, у нас до сих пор j=3 :)
	;i=k=0
				Acc		[X+i]	;добавили Txx
				ADD		[X+i^j]	;и Tyy
	;теперь наши метрики берём с весами	
				C		[SP]
				FMA		11010	;ManhW
				C		[SP+1]
				FMA		21758	;ChebW
				
	;ага, сделали
	;продемонстрируем результат
				JO		@@endLoop
		@@shl1:		ADD		Acc
				j++		0
				JNO		@@shl1		
	;ага, сделали
		@@endLoop:	DIV2		UAC ;получается отрицательное значение от -16384 до -1 (число -1 при "делении на 2" так и останется -1)
				SUB		32767 ;вычитает 32767, что даст отрицательное значение от -49151 до -32766. Но в UAC отобразится от 16385 до 32768.
				;если теперь найти обратную величину, это должно получиться в диапазоне 32768..65532
				;а если мы вычтем ещё половинку, может выйти от 16384 до 32768, и диапазон 32768..65536.
				;так что ну её нахрен, эту половинку?
				;или можем уже на этапе Ньютона этот случай проверить.
				DIV2S		1 ;это для нужд округления. 					
				[SP]		UAC							
	FindScale	endp


Содержание регистров к началу процедуры претерпели изменения, но уже не столь серьёзные. А ещё мы решили, что без команды ijk некуда, поэтому можно её здесь применить и немножко сэкономить.

И ещё одна вещь: в [SP+1] у нас до сих пор лежит число, позволяющее отличить "ближнюю дистанцию" от "дальней". В первом случае там лежит пятёрка (т.е точки от 0 до 5), во втором - тройка (т.е от 0 до 3). Не затереть бы его ненароком!

Collapse )

Компилируется всё это дело (алг. ближ дистанции) в 191 слово кода (382 байта). Возможно, здесь я чуть перемудрил с нахождением масштаба, можно львиную долю этого кода выполнить с незначительным снижением точности. Это когда появится алгоритм сопровождения, надо будет основательную отработку устроить, смотреть, удаётся ли "втянуться" в сопровождение с различных ракурсов и различных масштабов, или происходит срыв.

А сначала хоть этот кусок надо потестировать...
QuatCore neuro CPU

Нейросетевой процессор на Воронежской ПЛИС!

Всё я какие-то устаревшие решения пытаюсь применять для этого ВидеоИзмерителя Параметров Сближения. 16-битный одноядерный процессор, ассемблер, связанные списки (обычные, не блокчейн) - никто сейчас так не делает! То, что неделю назад поразвлекался с теорией оптимального обнаружения, БИХ-фильтрами и иже с ними - не сильно лучше, этим методам лет 60 как минимум :)

Пора уже вспомнить, что на дворе 2021 год - и сделать обнаружение точек на основе нейросети!
IMG20210401204530.jpg

Причём взять не какой-то там дряхлый персептрон, а самую нынче модную технологию - свёрточные нейросети (Convolutional Neural Networks, CNN), а если точнее - трёхслойную структуру, где первый слой отвечает за распознавание горизонтальных признаков, второй слой - за вертикальные признаки, а третий слой непосредственно обнаруживает пятна, т.е когда один из нейронов "вспыхивает", это означает: в этом месте у нас пятно, в смысле, мишень!

Если делать это "в лоб" на компьютере, то даже довольно мощные персоналки могут задуматься всерьёз и надолго. Но ПЛИС - совсем другое дело, и такая сеть может быть выполнена на Воронежской ПЛИС 5576ХС4Т, и работать в реальном времени (25 кадров в секунду 1024х720) при тактовой частоте 25 МГц :)

Collapse )
Doc

Более тяжёлые тела падают быстрее!

Увидел не так давно видео от Flammable Maths с таким заголовком, и подумал поначалу - он опять нас троллит. Это немецкий препод математики (насколько я знаю), и чувство юмора у него очень специфическое, особенно любит над инженерами издеваться, дескать e=π=3, cos(x) = 1, sin(x) = x.

Но нет, всё корректно: в безвоздушном пространстве, согласно законам Ньютона, тяжёлые тела действительно будут падать быстрее!

Быстренько изложу своими словами. Есть объект массой m и Земля массой M. На объект со стороны Земли действует сила



И согласно 2-му закону Ньютона она придаёт ему ускорение:



Тут можно вставить несколько страниц дискуссии, почему "инерционная масса" (слева) оказалась равна "гравитационной массе" (справа) и причём тут Эйнштейн, но сейчас мы о другом. Сокращаем их с чистой совестью, и получаем:



Ускорение зависит лишь от массы Земли и расстояния до неё, и не зависит от массы самого объекта, что как бы говорит все тела падают одинаково.

Вот только мы кое-чего забыли!

Collapse )