Случайны выбор дневника Раскрыть/свернуть полный список возможностей


Найдено 4540 сообщений
Cообщения с меткой

алгоритмы - Самое интересное в блогах

Следующие 30  »
rss_rss_hh_new

[Из песочницы] Алгоритм Левенберга — Марквардта для нелинейного метода наименьших квадратов и его реализация на Python

Пятница, 26 Августа 2016 г. 19:16 (ссылка)





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



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





Алгоритм Левенберга — Марквардта используется для решения нелинейного метода наименьших квадратов. Статья содержит:




  • объяснение алгоритма

  • объяснение методов: наискорейшего спуска, Ньтона, Гаусса-Ньютона

  • приведена реализация на Python с исходниками на github

  • сравнение методов



В коде использованы дополнительные библиотеки, такие как numpy, matplotlib. Если у вас их нет — очень рекомендую установить их из пакета Anaconda for Python



Зависимости



Алгоритм Левенберга — Марквардта опирается на методы, приведённые в блок-схеме







Поэтому, сначала, необходимо изучить их. Этим и займёмся



Определения






  • inline_formula — наша целевая функция. Мы будем минимизировать inline_formula. В этом случае, inline_formula является функцией потерь

  • inline_formula — градиент функции inline_formula в точке inline_formula

  • inline_formulainline_formula, при котором inline_formula является локальным минимумом, т.е. если существует проколотая окрестность inline_formula, такая что inline_formula

  • inline_formula — глобальный минимум, если inline_formula, т.е. inline_formula не имеет значений меньших, чем inline_formula

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

  • inline_formulaматрица Гессе (матрица вторых производных). Необходима для квадратичной аппроксимации



Выбор функции





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





Я принял inline_formula.Добавлен множитель 0.5 перед первой частью для удобства. Т.е. функция имеет вид:





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





Эта функция определена неотрицательно, имеет минимум inline_formula в точке inline_formula



В коде проще инкапсулировать все данные о функции в один класс и брать класс той функции, которая потребуется. Результат зависит от начальной точки оптимизации. Выберем её как inline_formula. Как видно из графика, в этой точке функция принимает наибольшее значение на интервале.



functions.py



class Rosenbrock:
initialPoint = (-2, -2)
camera = (41, 75)
interval = [(-2, 2), (-2, 2)]

"""
Целевая функция
"""
@staticmethod
def function(x):
return 0.5*(1-x[0])**2 + 0.5*(x[1]-x[0]**2)**2

"""
Для нелинейного МНК - возвращает вектор-функцию r
"""
@staticmethod
def function_array(x):
return np.array([1 - x[0] , x[1] - x[0] ** 2]).reshape((2,1))

@staticmethod
def gradient(x):
return np.array([-(1-x[0]) - (x[1]-x[0]**2)*2*x[0], (x[1] - x[0]**2)])

@staticmethod
def hesse(x):
return np.array(((1 -2*x[1] + 6*x[0]**2, -2*x[0]), (-2 * x[0], 1)))

@staticmethod
def jacobi(x):
return np.array([ [-1, 0], [-2*x[0], 1]])

"""
Векторизация для отрисовки поверхности
Детали: http://www.mathworks.com/help/matlab/matlab_prog/vectorization.html
"""
@staticmethod
def getZMeshGrid(X, Y):
return 0.5*(1-X)**2 + 0.5*(Y - X**2)**2


Метод наискорейшего спуска (Steepest Descent)



Сам метод крайне прост. Принимаем inline_formula, т.е. целевая функция совпадает с заданной.



Нужно найти inline_formula — направление наискорейшего спуска функции inline_formula в точке inline_formula.



inline_formula может быть линейно аппроксимирована в точке inline_formula:







где inline_formula — угол между вектором inline_formula. inline_formula следует из скалярного произведения



Так как мы минимизируем inline_formula, то чем больше разница в inline_formula, тем лучше. При inline_formula выражение будет максимально( inline_formula, норма вектора всегда неотрицательна), а inline_formula будет только если вектора inline_formula будут противоположны, поэтому





Направление у нас верное, но делая шаг длиной inline_formula можно уйти не туда. Делаем шаг меньше:





Теоретически, чем меньше шаг, тем лучше. Но тогда пострадает скорость сходимости. Рекомендуемое значение inline_formula



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



class Optimizer:
def _init_(self, function, initialPoint, gradient=None, jacobi=None, hesse=None,
interval=None, epsilon=1e-7, function_array=None, metaclass=ABCMeta):
self.function_array = function_array
self.epsilon = epsilon
self.interval = interval
self.function = function
self.gradient = gradient
self.hesse = hesse
self.jacobi = jacobi
self.name = self.class.name.replace('Optimizer', '')
self.x = initialPoint
self.y = self.function(initialPoint)

"Возвращает следующую координату по ходу оптимизационного процесса"
@abstractmethod
def next_point(self):
pass

"""
Движемся к следующей точке
"""
def move_next(self, nextX):
nextY = self.function(nextX)
self.y = nextY
self.x = nextX
return self.x, self.y


Код самого оптимизатора:



class SteepestDescentOptimizer(Optimizer):
...
def next_point(self):
nextX = self.x - self.learningRate * self.gradient(self.x)
return self.move_next(nextX)


Результат оптимизации:






























Итерация X Y Z
25 0.383 -0.409 0.334
75 0.693 0.32 0.058
532 0.996 0.990 inline_formula


Бросается в глаза: как быстро шла оптимизация в 0-25 итерациях, в 25-75 уже медленне, а в конце потребовалось 457 итераций, чтобы приблизиться к нулю вплотную. Такое поведение очень свойственно для МНС: очень хорошая скорость сходимости вначале, плохая в конце.



Метод Ньютона



Сам Метод Ньютона ищет корень уравнения, т.е. такой inline_formula, что inline_formula. Это не совсем то, что нам нужно, т.к. функция может иметь экстремум не обязательно в нуле.



А есть ещё Метод Ньютона для оптимизации. Когда говорят о МН в контексте оптимизации — имеют в виду его. Я сам, учась в институте, спутал по глупости эти методы и не мог понять фразу «Метод Ньютона имеет недостаток — необходимость считать вторые производные».



Рассмотрим для inline_formula



Принимаем inline_formula, т.е. целевая функция совпадает с заданной.



Разлагаем inline_formula в ряд Тейлора, только в отличии от МНС нам нужно квадратичное приближение:





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



Продифференцируем обе части по inline_formula. Наша цель, чтобы inline_formula, поэтому решаем уравнение:





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



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





Аналогично одномерному случаю — нужно проверить, правильно ли мы идём? Если матрица Гессе положительно определена, значит направление верное, иначе используем МНС.



В коде:



def is_pos_def(x):
return np.all(np.linalg.eigvals(x) > 0)

class NewtonOptimizer(Optimizer):
def next_point(self):
hesse = self.hesse(self.x)
# Если H - положительно определённая - Ок, иначе мы идём не в том направлении, используем МНС
if is_pos_def(hesse):
hesseInverse = np.linalg.inv(hesse)
nextX = self.x - self.learningRate * np.dot(hesseInverse, self.gradient(self.x))
else:
nextX = self.x - self.learningRate * self.gradient(self.x)

return self.move_next(nextX)


Результат:




























Итерация X Y Z
25 -1.49 0.63 4.36
75 0.31 -0.04 0.244
179 0.995 -0.991 inline_formula


Сравните с МНС. Там был очень сильный спуск до 25 итерации (практически упали с горы), но потом сходимость сильно замедлилась. В МН, напротив, мы сначала медленно спускаемся с горы, но затем движемся быстрее. У МНС ушло с 25 по 532 итерацию, чтобы дойти до нуля с inline_formula. МН же оптимизировал inline_formula за 154 последних итераций.



Это частое поведение: МН обладает квадратичной скоростью сходимости, если начинать с точки, близкой к локальному экстремуму. МНС же хорошо работает далеко от экстремума.



МН использует информацию кривизны, что было видно на рисунке выше (плавный спуск с горки). Ещё пример, демонстрирующий эту идею: на рисунке ниже красный вектор — это направление МНС, а зелёный — МН







[Нелинейный vs линейный] метод наименьших квадратов



В МНК у нас есть модель inline_formula, имеющая inline_formula параметров, которые настраиваются так, чтобы минимизировать inline_formula, где inline_formulainline_formula-е наблюдение.



В линейном МНК у нас есть inline_formula уравнений, каждое из которых мы можем представить как линейное уравнение





Для линейного МНК решение единственно. Существуют мощные методы, такие как QR декомпозиция, SVD декомпозиция, способные найти решение для линейного МНК за 1 приближённое решение матричного уравнения inline_formula.



В нелинейном МНК параметр inline_formula может сам быть представлен функцией, например inline_formula. Так же, может быть произведение параметров, например inline_formula и т.д.



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



Методы ниже имеют дело как раз с нелинейным случаем. Но, сперва, рассмотрим нелиненый МНК в контексте нашей задачи — минимизации функции





Ничего не напоминает? Это как раз форма МНК! Введём вектор-функцию inline_formula





и будем подбирать inline_formula так, чтобы решить систему уравнений (хотя бы приближённо):





Тогда нам нужна мера — насколько хороша наша аппроксимация. Вот она:





Я применил обратную операцию: подстроил вектор-функцию inline_formula под целевую inline_formula. Но можно и наоборот: если дана вектор-функция inline_formula, строим inline_formula из (5). Например:





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



Метод Гаусса-Ньютона



Метод основан на всё той же линейной аппроксимации, только теперь имеем дело с двумя функциями:





Далее делаем то же, что и в методе Ньютона — решаем уравнение (только для inline_formula):





Несложно показать, что вблизи inline_formula:





Код оптимизатора:



class NewtonGaussOptimizer(Optimizer):
def next_point(self):
# Решаем (J_t * J)d_ng = -J*f
jacobi = self.jacobi(self.x)
jacobisLeft = np.dot(jacobi.T, jacobi)
jacobiLeftInverse = np.linalg.inv(jacobisLeft)
jjj = np.dot(jacobiLeftInverse, jacobi.T) # (J_t * J)^-1 * J_t
nextX = self.x - self.learningRate * np.dot(jjj, self.function_array(self.x)).reshape((-1))
return self.move_next(nextX)


Результат превысил мои ожидания. Всего за 3 итерации мы пришли в точку inline_formula. Чтобы продемонстрировать траекторию движения я уменьшил learningrate до 0.2







Алгоритм Левенберга — Марквардта



Он основан на одной из версий Методе Гаусса-Ньютона ("damped version" ):





Для больших inline_formula получается метод наискорейшего спуска, для маленьких — метод Ньютона.

Сам алгоритм в процессе оптимизации подбирает нужный inline_formula на основе gain ratio, определяющийся как:





Если inline_formula, то inline_formula — хорошая аппроксимация для inline_formula, иначе — нужно увеличить inline_formula.



Начальное значение inline_formula задаётся как inline_formula, где inline_formula — элементы матрицы inline_formula.



inline_formula рекомендовано назначать за inline_formula. Критерием остановки является достижение глобального минимуму, т.е. inline_formula





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



class LevenbergMarquardtOptimizer(Optimizer):
def __init__(self, function, initialPoint, gradient=None, jacobi=None, hesse=None,
interval=None, function_array=None, learningRate=1):
self.learningRate = learningRate
functionNew = lambda x: np.array([function(x)])
super().__init__(functionNew, initialPoint, gradient, jacobi, hesse, interval, function_array=function_array)
self.v = 2
self.alpha = 1e-3
self.m = self.alpha * np.max(self.getA(jacobi(initialPoint)))

def getA(self, jacobi):
return np.dot(jacobi.T, jacobi)

def getF(self, d):
function = self.function_array(d)
return 0.5 * np.dot(function.T, function)

def next_point(self):
if self.y==0: # Закончено. Y не может быть меньше нуля
return self.x, self.y

jacobi = self.jacobi(self.x)
A = self.getA(jacobi)
g = np.dot(jacobi.T, self.function_array(self.x)).reshape((-1, 1))
leftPartInverse = np.linalg.inv(A + self.m * np.eye(A.shape[0], A.shape[1]))
d_lm = - np.dot(leftPartInverse, g) # moving direction
x_new = self.x + self.learningRate * d_lm.reshape((-1)) # линейный поиск
grain_numerator = (self.getF(self.x) - self.getF(x_new))
gain_divisor = 0.5* np.dot(d_lm.T, self.m*d_lm-g) + 1e-10
gain = grain_numerator / gain_divisor
if gain > 0: # Ок, хорошая оптимизация
self.move_next(x_new) # ok, шаг принят
self.m = self.m * max(1 / 3, 1 - (2 * gain - 1) ** 3)
self.v = 2
else:
self.m *= self.v
self.v *= 2

return self.x, self.y


Результат получился тоже хороший:


























Итерация X Y Z
0 -2 -2 22.5
4 0.999 0.998 inline_formula
11 1 1 0


При learningrate =0.2:



Сравнение методов







































Название метода Целевая функция Достоинства Недостатки Сходимость
Метод наискорейший спуск дифференцируемая -широкий круг применения

-простая реализация



-низкая цена одной итерации

-глобальный минимум ищется хуже, чем в остальных методах



-низкая скорость сходимости вблизи экстремума

локальная
Метод Нютона дважды дифференцируемая -высокая скорость сходимости вблизи экстремума



-использует информацию о кривизне

-функция должны быть дважды дифференцируема



-вернёт ошибку, если матрица Гессе вырождена (не имеет обратной)



-есть шанс уйти не туда, если находится далеко от экстремума

локальная
Метод Гаусса-Нютона нелинейный МНК -очень высокая скорость сходимости



-хорошо работает с задачей curve-fitting

-колонки матрицы J должны быть линейно-независимы



-налагает ограничения на вид целевой функции



локальная
Алгоритм Левенберга — Марквардта нелинейный МНК -наибольная устойчивость среди рассмотренных методов



-наибольшие шансы найти глобальный экстремум



-очень высокая скорость сходимости (адаптивная)



-хорошо работает с задачей curve-fitting

-колонки матрицы J должны быть линейно-независимы



-налагает ограничения на вид целевой функции



-сложность в реализации

локальная


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



Совмещённый результат (специально понижена скорость последних двух методов):







» Исходники можно скачать с github



» Источники




  1. K. Madsen, H.B. Nielsen, O. Tingleff (2004): Methods for non-linear least square

  2. Florent Brunet (2011): Basics on Continuous Optimization

  3. Least Squares Problems


Original source: habrahabr.ru.

https://habrahabr.ru/post/308626/?utm_source=habrahabr&utm_medium=rss&utm_campaign=best

Метки:   Комментарии (0)КомментироватьВ цитатник или сообщество
rss_rss_hh_new

Навигатор 2ГИС: Экстраполяция позиции автомобиля

Пятница, 26 Августа 2016 г. 10:26 (ссылка)





В приложении 2ГИС теперь есть навигатор. Мы научились «ехать» по треку, озвучивать манёвры, автоматически перестраивать маршрут, рассчитывать время в пути, доводить пользователя до входа в здание или организацию, учитывая заборы и шлагбаумы, — и всё это в честном офлайне. Пробки (вот разве что для них нужен интернет), разведённые мосты и перекрытые улицы учитываем давно. Пока в нашем навигаторе — необходимый минимум. Чуть позже научим его предупреждать о слишком высокой скорости, лежачих полицейских и камерах ГИБДД, настроим ночной режим, сделаем маршруты по платным и грунтовым дорогам опциональными. Чтобы воспользоваться им, нужно обновить 2ГИС в своем смартфоне или скачать в AppStore или Windows Store. Для Android обновление выходит постепенно, начиная с 22 августа (будет доступно на всю аудиторию к сентябрю).



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



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



Маркер GPS и маршрут



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



В самом простом случае можно отображать позицию устройства на местности, считывая координаты с датчика GPS и размещая маркер в соответствующее место на карте. Уже здесь мы сталкиваемся с первой проблемой — измерительной погрешностью, которая даже в условиях неплохого сигнала вполне может достигать 20–30 метров.



Для ответа на обычный вопрос «Где я нахожусь?» такого способа отображения будет вполне достаточно, особенно если вокруг маркера нарисовать ещё и круг точности с радиусом, равным оценке погрешности. Однако для навигации нужно придумать что-то получше, ведь водителя, движущегося по городской улице, вряд ли устроит маркер GPS, расположенный внутри соседнего дома или, того хуже, на каком-нибудь внутриквартальном проезде.



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



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



Отображение геопозиции во времени



Дискретный характер поступления данных от датчика GPS — ещё одна проблема при реализации ведения пользователя по маршруту. В идеальном случае координаты обновляются один раз в секунду. Рассмотрим несколько вариантов отображения геопозиции во времени и выберем наиболее подходящий для наших задач.



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



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



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



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



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



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



3. С учётом появившегося требования к точности позиционирования стоит заметить, что теперь от нас требуется незадолго до прихода нового отсчёта GPS расположить маркер в точке, максимально приближенной к этому новому отсчёту. То есть, по сути, заглянуть будущее, пусть и ненадолго. Хотя с изобретением машины времени у человечества пока дела обстоят из рук вон плохо, для нас спасение всё же есть. Движение автомобиля инертно, поэтому скорость и направление его движения не могут меняться мгновенно, а раз так, мы можем попытаться с некоторой точностью спрогнозировать, где пользователь будет находиться в интервале между последним отсчётом позиции и будущим. Если нам удастся добиться того, что ошибка прогнозирования в большинстве случаев будет меньше, чем погрешность второго способа, то мы здорово облегчим жизнь пользователей нашего навигатора.



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



Ведение по маршруту с экстраполяцией позиции



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



Вспомним поступающие к нам данные и введём для них обозначения:

s_i, i\in \mathbb{N} — реальные отсчёты смещения, получаемые притяжкой позиции GPS к линии маршрута;

t_i, i\in \mathbb{N} — время прихода соответствующих отсчётов смещения.

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



В конечном итоге нам необходимо построить функцию экстраполяции смещения s=s(t),t \in [t_0; +\infty), которая будет приближена к реальной динамике автомобиля и при этом обеспечит плавность движения маркера GPS по всему нашему маршруту (его длина ни на что не повлияет, так как завершение маршрута обрабатывается отдельно, поэтому условно будем считать маршрут бесконечным). Для обеспечения хорошей визуальной плавности достаточно будет условия гладкости s(t), то есть ни позиция, ни скорость маркера не должны меняться скачком. Другими словами, функция s(t) обязана быть непрерывной вместе со своей первой производной (здесь и далее — по времени) на всей области определения.



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



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



Построим такую функцию экстраполяции смещения s_i (t),t\in[t_i;+\infty), чтобы после i-го отсчёта её значения предсказывали реальное расположение пользователя в течение достаточного времени до прихода (i+1)-го отсчёта. Все полезные данные, которыми мы обладаем, — последовательность отсчётов до i-го включительно вместе со временем получения каждого из них.



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



v_i\approx \frac{s_i-s_{i-1}}{t_i-t_{i-1}}=s_i{^'} (t_i)


, где v_i — оценка скорости по отсчётам, а s_i{^'} (t_i) — производная экстраполяционной функции s_i (t), которую мы пытаемся построить.



Аналогично для производных более высокого порядка — ускорения, рывка и т.д.:



\left\{\begin{aligned}<br />
a_i\approx \frac{v_i-v_{i-1}} {t_i-t_{i-1}} = s_i^{''} \\<br />
j_i\approx \frac{a_i-a_{i-1}} {t_i-t_{i-1}} = s_i^{'''} \\<br />
\end{aligned}<br />
\right.


Как видно из этих формул, для получения оценки всё более старших производных смещения требуется учитывать всё больше отсчётов, предшествующих текущему: для определения скорости нужны два отсчёта, для ускорения — три, для рывка — четыре и т.д. С одной стороны, чем больше динамических характеристик движения мы будем учитывать в своём прогнозе, тем большую моделирующую способность получим; с другой — полезная информация, содержащаяся во всё более «старых» отсчётах, драматически теряет актуальность. Например, тот факт, что мы ехали со скоростью 30 км/ч минуту назад ничем не поможет нам в текущий момент времени: с тех пор можно было несколько раз разогнаться, затормозить или вообще остановиться. По этой причине оценки всё более старших производных смещения становятся всё дальше от реальности; кроме того, вклад погрешности вычисления некоторой производной в общую аналитическую модель смещения тоже растёт с увеличением порядка этой производной. Раз так, то, начиная с какого-то порядка, динамические характеристики, оценённые при помощи конечных разностей, вместо уточнения будут только портить нашу модель.



По результатам проверок на реальных данных выяснилось, что оценка рывка j_i, особенно в случаях «среднего» качества сигнала GPS, уже достаточно плоха, чтобы от неё было больше вреда, чем пользы. С другой стороны, к счастью, наиболее частые сценарии динамики автомобиля — это покой, равномерное и равнопеременное движение, описываемые полиномиальными уравнениями 0-й, 1-й и 2-й степени от времени соответственно.



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



s_i (t)=s_i+v_i t+\frac12 a_i t^2


Осталось сделать всего один шаг: область определения s_i (t) начинается с момента времени t_i, поэтому время в вычислениях удобнее считать с этого же момента.



В итоге функция примет вид:



s_i (t)=s_i+v_i (t-t_i )+\frac12 a_i (t-t_i )^2,t

https://habrahabr.ru/post/308500/?utm_source=habrahabr&utm_medium=rss&utm_campaign=best

Метки:   Комментарии (0)КомментироватьВ цитатник или сообщество
rss_rss_hh_new

[Из песочницы] Reverse engineering тестового crackme от Лаборатории Касперского

Четверг, 25 Августа 2016 г. 18:26 (ссылка)

Приветствую сообщество! Давным давно, в 2013 году на Хабре был опубликован пост «Reverse engineering на собеседовании: как мы нанимаем на работу». В нём был предложен тестовый crackme для претендентов на позицию вирусного аналитика. Убедившись, что полного разбора тестового файла в интернете нет, я решил написать свой разбор. И так, приступим. Crackme 64-разрядный. Запустим его в IDA Pro.



image




Видим слева в списке функций три функции: start — функция, с которой начинается выполнение программы, DialogFunc — эта функция общается с нами и некоторая функция sub_140001000. Рассмотрим диалоговую функцию. Декомпилируем её Hex Rays-ом.



image


В глаза бросается ветвление условий, согласно которому, если некоторая функция sub_140001000 вернет TRUE, то появится сообщение, информирующее нас о отлично проделанной работе, иначе неверно. Разберем нашу заветную функцию sub_140001000. Если мы пропустим её через декомпилятор, то увидим, что в качестве аргумента передается указатель всего на одно значение. Вероятно, это значение берется из диалогового окна и является вводимым ключом. Теперь рассмотрим ассемблерный листинг. Имеется первая проверка условия верности введенных данных. Если условие выполняется, то программа исполняется дальше, если не выполняется, то идет возврат из подпрограммы.



image


Запустим наш crackme под отладчиком. Воспользуемся x64dbg. Поставим breakpoint на нашей первой проверке. В качестве вводимого ключа используем набор цифр 1234567.



image


Как видно, происходит проверка значения регистра edx и числа 13h (в десятичной системе исчисления это 19). Вероятно, это проверка на количество введенных знаков ключа (у нас их 7 и в регистре edx число 7). Попробуем ввести другое количество символов. Запустим отладчик заново. Введем 9 цифр 123456789.



image


Похоже, что так оно и есть. Значит наш ключ должен содержать 19 символов. Вводим 19 символов 1234567890123456789 и переходим к следующему этапу проверки.



image


На этом шаге осуществляется проверка каждого пятого символа ключа на равенство значению 2Dh. Дело в том, что число 2Dh — это шестнадцатиричный код символа "-". Т.е. наш ключ должен иметь вид xxxx-xxxx-xxxx-xxxx. Используем в качестве ключа 1234-5678-9012-3456 и переходим к следующему шагу.



image


А на следующем шаге происходит проверка символов на числовую принадлежность. Порядок проверки такой: берется символ из ключа (в счет не идут каждый пятый символ ключа) и к его шестнадцатиричному коду прибавляется число -30 и полученный результат сравнивается с числом 9. Если меньше, то на проверку берется следующий символ ключа, если больше, то выводится сообщение, что ключ неверный. Идем дальше.



image


На этом шаге осуществляется проверка того, чтобы суммы чисел в блоках были равны. На рисунке выше выделен блок кода, который производит подсчет суммы чисел и область стека, куда эти суммы заносятся. Параллельно суммы блоков суммируются между собой и заносятся в регистр r10. Далее идет деление результата в регистре r10 на 4 (shr r10d,2 — сдвиг на 2 байта равносилен делению на 4) и сравнение значения из регистра r10 с ранее занесенными значениями в стеке. Отлично. Делаем так, чтобы суммы цифр каждого блока ключа были равны (например 1122-0123-2112-0006) и двигаемся дальше на следующий шаг проверки.



image


Участок кода, выделенный на рисунке выше, проверяет, чтобы расположение символов в каждом последующем блоке ключа не совпадало с предыдущим. В итоге наш ключ имеет вид 1478-7814-1478-7814. Проверяем.



image


image


Отличная работа!
Original source: habrahabr.ru.

https://habrahabr.ru/post/308548/?utm_source=habrahabr&utm_medium=rss&utm_campaign=best

Метки:   Комментарии (0)КомментироватьВ цитатник или сообщество
rss_rss_hh_new

[Из песочницы] Нейронные сети для любопытных программистов (с примером на c#)

Среда, 24 Августа 2016 г. 17:57 (ссылка)

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



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



пример на java и полезные ссылки

наглядная реализацыя с исползованием ООП



Поскольку теории очень много по этой теме хотелось бы приступить к реализации.



Реализация


using UnityEngine;
using System.Collections;
using System.Xml.Serialization;

public class Neuron {

[XmlAttribute("weight")]
public string data;

[XmlIgnore]
public int[,] weight; // веса нейронов

[XmlIgnore]
public int minimum = 50; // порог

[XmlIgnore]
public int row = 64,column = 64;

/**
* Конструктор нейрона, создает веси и устанавливает случайные значения
*/
public Neuron()
{
weight = new int[row,column];
randomizeWeights();
}

/**
* ответы нейронов, жесткая пороговая
* @param input - входной вектор
* @return ответ 0 или 1
*/
public int transferHard(int[,] input)
{
int Power = 0;
for(int r = 0; r < row;r++)
{
for(int c = 0; c < column;c++)
{
Power += weight[r,c]*input[r,c];
}
}
//Debug.Log("Power: " + Power);
return Power >= minimum ? 1 : 0;
}

/**
* ответы нейронов с вероятностями
* @param input - входной вектор
* @return n вероятность
*/
public int transfer(int[,] input)
{
int Power = 0;
for(int r = 0; r < row;r++)
for(int c = 0; c < column;c++)
Power += weight[r,c]*input[r,c];

//Debug.Log("Power: " + Power);
return Power;
}

/**
* устанавливает начальные произвольные значения весам
*/
void randomizeWeights()
{
for(int r = 0; r < row;r++)
for(int c = 0; c < column;c++)
weight[r,c] = Random.Range(0,10);
}

/**
* изменяет веса нейронов
* @param input - входной вектор
* @param d - разница между выходом нейрона и нужным выходом
*/
public void changeWeights(int[,] input,int d)
{
for(int r = 0; r < row;r++)
for(int c = 0; c < column;c++)
weight[r,c] += d*input[r,c];
}

public void prepareForSerialization()
{
data = "";
for(int r = 0; r < row;r++)
{
for(int c = 0; c < column;c++)
{
data += weight[r,c] + " ";
}
data += "\n";
}
}

public void onDeserialize()
{
weight = new int[row,column];

string[] rows = data.Split(new char[]{'\n'});
for(int r = 0; r < row;r++)
{
string[] columns = rows[r].Split(new char[]{' '});
for(int c = 0; c < column;c++)
{
weight[r,c] = int.Parse(columns[c]);
}
}
}
}





Класс нейронов который содержит weight — двоичный массив весов, minimum — порог функции. Функция transferHard возвращает ответ на входной вектор. Поськольку ответ функции жесткий, я использоваю его для обучения. На мой взгляд это более эффективно обучает нейроны. Я буду очень блогодарен если будут отзывы по этому поводу. Функция transfer возвращает ответ на входной вектор но с вероятностю, сумма может быть ближе к нулю или отрицательной если нейрон обучен для другово символа.



using UnityEngine;
using System.Collections;
using System.Xml.Serialization;
using System.Xml;
using System.IO;

public class NeuralNetwork {

[XmlArray("Neurons")]
public Neuron[] neurons;

/**
* Конструктор сети создает нейроны
*/
public NeuralNetwork()
{
neurons = new Neuron[10];

for(int i = 0;i**
* функция распознавания символа, используется для обучения
* @param input - входной вектор
* @return массив из нуллей и единиц, ответы нейронов
*/
int[] handleHard(int[,] input)
{
int[] output = new int[neurons.Length];
for(int i = 0;i**
* функция распознавания символа, используется для конечново ответа
* @param input - входной вектор
* @return массив из вероятностей, ответы нейронов
*/
int[] handle(int[,] input)
{
int[] output = new int[neurons.Length];
for(int i = 0;i**
* ответ сети
* @param input - входной вектор
* @return индекс нейронов предназначенный для конкретного символа
*/
public int getAnswer(int[,] input)
{
int[] output = handle(input);
int maxIndex = 0;
for(int i = 1; i < output.Length;i++)
if(output[i] > output[maxIndex])
maxIndex = i;

return maxIndex;
}

/**
* функция обучения
* @param input - входной вектор
* @param correctAnswer - правильный ответ
*/
public void study(int[,] input,int correctAnswer)
{
int[] correctOutput = new int[neurons.Length];
correctOutput[correctAnswer] = 1;

int[] output = handleHard(input);
while(!compareArrays(correctOutput,output))
{
for(int i = 0; i < neurons.Length;i++)
{
int dif = correctOutput[i]-output[i];
neurons[i].changeWeights(input,dif);
}
output = handleHard(input);
}
}

/**
* сравнение двух вектор
* @param true - если массивы одинаковые, false - если нет
*/
bool compareArrays(int[] a,int[] b)
{
if(a.Length != b.Length)
return false;

for(int i = 0;iNeuralNetwork.txt", FileMode.Create);
XmlWriter writer = new XmlTextWriter(stream, new System.Text.ASCIIEncoding());
using(writer)
{
serializer.Serialize(writer, this);
}
}

public static NeuralNetwork fromXml()
{
string xml = "";
FileStream fStream = new FileStream(Application.dataPath + "/NeuralNetwork.txt",
FileMode.OpenOrCreate);
if(fStream.Length > 0)
{
byte[] tempData = new byte[fStream.Length];
fStream.Read(tempData, 0, tempData.Length);

xml = System.Text.Encoding.ASCII.GetString(tempData);
}
fStream.Close();

if(string.IsNullOrEmpty(xml))
return new NeuralNetwork();


NeuralNetwork data;

XmlSerializer serializer = new XmlSerializer(typeof(NeuralNetwork));
using(TextReader reader = new StringReader(xml))
{
data = serializer.Deserialize(reader) as NeuralNetwork;
}

data.onDeserialize();

return data;
}
}





Класс NeuralNetwork содержит массив нейронов, каждый из них предназначен для конкретного символа. В конструкторе создается массив из десяти элементов, потому что пример сделан для разпознования цифр(0-9). Если вы хотите использовать сеть для распознавания букв то поменяйте размер массива соответствующим образом. Функция handleHard вызывается для обучение нейронов, возвращает массив из нуллей и единиц. Функция getAnswer ответ сети для входного вектора, использует функцию handle для получения массива ответов, каждый элемент массива содержит ответ нейрона с вероятностью. Функция getAnswer выбирает индекс элемента который содержит наибольшую вероятность и возвращает эго как ответ.



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



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



Заключение


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



Хотелось увидеть различные отзывы на тему и на пост.
Original source: habrahabr.ru (comments, light).

https://habrahabr.ru/post/308448/

Метки:   Комментарии (0)КомментироватьВ цитатник или сообщество
rss_rss_hh_new

[Перевод] Музыка, Mathematica и вычислительная вселенная: автоматическое создание музыки на основе клеточных автоматов

Среда, 24 Августа 2016 г. 17:03 (ссылка)



Перевод поста Стивена Вольфрама (Stephen Wolfram) "Music, Mathematica, and the Computational Universe" о замечательном ресурсе WolframTones, работа которого была недавно возобновлена на новой площадке Wolfram Cloud (сайт, созданный в 2005 г., был недоступен пару лет, так как использовал не поддерживаемые современными браузерами решения).

Выражаю огромную благодарность Кириллу Гузенко за помощь в переводе.



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



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



Но что есть творчество? Это то, что было необходимо в течение всей биологической и культурной эволюции? И может ли оно также существовать в системах, которые не имеют ничего общего с людьми?



В своей работе над книгой Новый вид науки (A New Kind of Science) я исследовал вычислительную вселенную возможных программ и обнаружил, что даже очень простые программы могут показывать поразительно богатый и сложный характер, наравне, например, с тем, что можно встретить в природе. И, опираясь на разработанный принцип вычислительной эквивалентности, я пришел к убеждению, что не может быть ничего, что принципиально отличает наши человеческие способности от любых процессов, которые происходят в природе, или даже в очень простых программах.



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



И мне стало любопытно: действительно ли музыка есть что-то особенное и исключительно человеческое? Или всё таки её можно прекрасно создавать автоматически, с помощью вычислений?



В 2003 году, после десятка лет моей затворнической работы над A New Kind of Science, мне постепенно переставали быть чужды мирские проблемы, и одна из них заключалась в том, что рингтон на моём сотовом был как у всех. Так что я подумал: если какая-то оригинальная индивидуализированная музыка на самом деле может создаваться автоматически, то можно было бы просто поменять всем мелодии на мобильном, и каждый бы имел свою собственную.



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



Получилась долгая история попыток создания музыки по правилам. Большинство полученного производило впечатление слишком "роботского" или случайного. Но то, к чему я пришёл в A New Kind of Science, казалось, давало нам новые возможности, ведь там было убедительно показано, что даже с правилами простых программ можно создавать такие сложные и живые вещи, которые мы наблюдаем в природе и которыми мы восхищаемся.



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



Факт в том, что за логикой в музыке всегда стоит какая-то простая программа. Но ключевым моментом A New Kind of Science является то, что основная программа может быть простой, однако производить богатую и сложную картину.



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



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



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



Проходящие мимо нашего офиса люди, слыша воспроизводящуюся музыку, часто останавливались и спрашивали: "Что это за песню вы слушаете?" Мы делали музыку, которая была достаточно хороша для того, чтобы люди считали её сделанной людьми; можно сказать, что мы успешно прошли аналог теста Тьюринга для музыки.



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







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



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



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



image

(Прослушать аудио, соответствующее данному клеточному автомату)



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







Сайт WolframTones был запущен 16 сентября 2005 года. И вот, с тех пор он есть в интернете и с помощью Mathematica создаёт музыку. Должен признаться, что я некоторое время не смотрел его логи. Однако заглянув туда сейчас, я обнаружил, что он использовался десятки миллионов раз для создания десятков миллионов музыкальных композиций.



image

(Прослушать аудио, соответствующее данному клеточному автомату)



По меркам, к примеру, использования Wolfram|Alpha — это ничто. Но для музыкальных композиций это огромная цифра. Itunes в настоящее время содержит около 14 миллионов музыкальных композиций и представляет собой большинство того, что когда-либо было выпущено.  Но всего за несколько лет WolframTones создал большее количество композиций. Используя лишь вычисления, он в некотором смысле превзошёл наш вид в создании музыкальной продукции, в одиночку создав больше оригинальной музыки, чем за всё время до него.



Для того, чтобы результат выдавался мгновенно, сайт кодирует музыку в MIDI (то, что Mathematica теперь поддерживает напрямую в символьном виде). Было создано множество композиций WolframTones в MP3 формате. Полным перераспределением ролей можно назвать случай, когда я несколько лет назад ходил на концерт, где люди играли на скрипках фрагмент, полностью созданный с помощью WolframTones.



Могут ли простые программы полностью создать целую симфонию? Композиции в WolframTones ведают нам некоторые истории из вычислительного мира, которые имеют продолжительность, пожалуй, не более минуты. Судя по моему опыту, для того, чтобы создать больший фрагмент — тот, что поведает нам какую-то более длинную историю — нам потребуются структуры более высокого уровня. Однако у нас нет каких-то препятствий в создании простой программы, которая бы обеспечила подобную структуру. Главенствующая история может быть абсолютно полной, как и многие другие истории, что звучат в нашем мире музыкой законов природы.



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



Очень легко начать создавать программы в Mathematica:



Sound[
Map[
Function[
SoundNote[DeleteCases[3 * Range[21] * Reverse[#], 0] - 24,
0.1
]
],
Transpose @ CellularAutomaton[90, {{1}, 0}, 20]
]
]






(Прослушать аудио, соответствующее данному клеточному автомату)



Sound[
Map[Function[SoundNote[#, 1 / 6, "Warm"]],
Map[Pick[{0, 5, 9, 12, 16, 21}, #, 1] &,
CellularAutomaton[30,
{{1, 0, 0, 0, 0, 0}, 0},
13,
{13, 5}
]
]
]
]






(Прослушать аудио, соответствующее данному клеточному автомату)



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



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



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



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



Когда я смотрю на Wolfram|Alpha, меня захватывает радость от осознания того, как много человеческих знаний можно запечатлеть и сделать их вычислимыми. Новый рубеж — запечатлеть не только знания, но и творчество. Чтобы иметь возможность, например, добиваясь некоторой цели, понять творческий путь её достижения. Музыка есть суть плод чистого творчества, и то, что мы узнали, как предполагает принцип вычислительной эквивалентности, что даже в этой области, такие вещи, как WolframTones, удивительно хороши в получении творчески насыщенного результата.



Мы хотим воплотить более высокий уровень автоматизации, тем самым резко расширив творческие возможности, что, без сомнения, добавит множество новых увлекательных возможностей.
Original source: habrahabr.ru.

https://habrahabr.ru/post/308428/?utm_source=habrahabr&utm_medium=rss&utm_campaign=best

Комментарии (0)КомментироватьВ цитатник или сообщество
rss_rss_hh_new

А вот про нейронные сети, ИИ и т.д. [Опрос]

Вторник, 23 Августа 2016 г. 20:43 (ссылка)

На просторах интернета часто доводилось видеть мнения, что де "нейросеть — панацея от всего и вся", т.е. например "натравите нейросеть — и все, профит" или еще брутальней "скоро создадут ИИ на базе нейронной сети, которая сможет заменить даже программистов / администраторов / аналитиков и т.д.".

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



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

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



Написать статью (и опрос) хотел уже довольно давно, но все как-то руки не доходили. А после очередного вопроса-предложения по е-майлу "натравите же нейронную сеть" на проблему из моей прошлой статьи "Мониторинг лог-журналов: Такой уязвимый лог...", все-таки понял — нет — надо писать.



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

Т.е. чтобы можно было ткнуть носом в статью очередного такого "все-лучше-всех-знающего" от интернета...



Хотя вдруг это я таки нещадно отстал от жизни...



Итак, для начала озвучим про или скорее контра использования нейронных сетей, что касается ее "интеллекта":




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

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

  • данные из обучающей выборки, на которых собственно сеть обучается, должны быть:


    • репрезентативны (должны иллюстрировать истинное положение вещей в предметной области);

    • непротиворечивы (противоречивые данные приведут к плохому качеству обучения сети);

    • преобразованы к виду который можно отправить входным нейронам (т.е. должны быть как минимум первично обработаны другими алгоритмами)


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

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



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

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



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



Кроме, того нередко результат работы такого "алгоритма" есть вероятностная оценка (или группа их).

Что я хочу этим сказать. Ну возьмем к примеру лог-строчку из вышеупомянутого поста:



Aug 18 08:04:51 srv sshd[2131]: Failed password for invalid user test from 1.2.3.4 port 46589 ssh2 from 4.3.2.1 port 58946 ssh2


Ну оценит она это как "неавторизованная попытка с хоста 1.2.3.4 с вероятностью 99%", если она обучалась на миллионах строк, где после первого " from " всегда стоит плохой адрес (и страшно ошибется). Или в лучшем случае — пропустит ее как "мусорную" строку или найдет там оба адреса (что как минимум должно быть предусмотрено в самом flow ее создателем).



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



Как заставить ее "думать" как человек (т.е. начать как минимум сомневаться в том что видим) и главное принять в итоге правильное решение — это как раз очень сложная задача, которую пока решает разработчик (и нередко алгоритмически, вне или в обвязке нейросети) в каждом конкретном случае.

И это, как мне думается, еще очень и очень надолго.



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



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



Если же вы считаете, что создали уже что-то такое мега-умное, у меня для вас есть новости:




  • вы уже безработный

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

  • вас ненавидит лучшая часть человечества (другими словами вы выбрали себе врагами миллионы умнейших и образованнейших людей планеты)

  • скайнет уже близко (и скоро нас всех поработит)

  • за вами вероятно уже выехали (поэтому продумайте таки вариант с Бора-Бора)



Теперь собственно к опросам.



Чтобы исключить "шум" в опросах, они поделены на три группы (две для тех кто знает, и одна — для предположений), по два опроса для каждой целевой аудитории:




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

  • для теоретиков, т.е. людей, разбирающихся в теории, как это работает (но ни разу не использовавших это в бою)

  • предположения: для всех остальных (ну дайте уже кликнуть где-нибудь).



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



Итак, поехали!





[Для профессионалов] Насколько «умны» современные нейронные сети и системы на их базе (что они для вас)


































































































Проголосовало 48 человек. Воздержалось 106 человек.





[Для профессионалов] Смогут ли ИИ в обозримом будующем заменить разработчика


























































Проголосовало 58 человек. Воздержалось 96 человек.





[Для теоретиков] Насколько «умны» современные нейронные сети и системы на их базе (что они для вас)


































































































Проголосовал 41 человек. Воздержалось 89 человек.





[Для теоретиков] Смогут ли ИИ в обозримом будующем заменить разработчика


























































Проголосовало 48 человек. Воздержалось 90 человек.





[Как вы предполагаете] Насколько «умны» современные нейронные сети и системы на их базе (что они для вас)


































































































Проголосовало 57 человек. Воздержалось 76 человек.





[Как вы предполагаете] Смогут ли ИИ в обозримом будующем заменить разработчика


























































Проголосовало 67 человек. Воздержалось 73 человека.





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


Original source: habrahabr.ru.

https://habrahabr.ru/post/308356/?utm_source=habrahabr&utm_medium=rss&utm_campaign=best

Метки:   Комментарии (0)КомментироватьВ цитатник или сообщество
rss_rss_hh_new

Быстрее быстрого или глубокая оптимизация Медианной фильтрации для GPU Nvidia

Вторник, 23 Августа 2016 г. 10:57 (ссылка)

Введение



В предыдущем посте я постарался описать, как легко можно воспользоваться преимуществом GPU для обработки изображений. Судьба сложилась так, что мне подвернулась возможность попробовать улучшить медианную фильтрацию для GPU. В данном посте я постараюсь рассказать каким образом можно получить еще больше производительности от GPU в обработке изображений, в частности, на примере медианной фильтрации. Сравнивать будем GPU GTX 780 ti с оптимизированным кодом, запущенном на современном процессоре Intel Core i7 Skylake 4.0 GHz с набором векторных регистров AVX2. Достигнутая скорость фильтрации квадратом 3х3 в 51 GPixels/sec для GPU GTX 780Ti и удельная скорость фильтрации квадратом 3х3 в 10.2 GPixels/sec на 1 TFlops для одинарной точности на данное время являются самыми высокими из всех известных в мире.



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




  1. Для каждой точки исходного изображения берется некоторая окрестность (в нашем случае 3x3 / 5х5).

  2. Точки данной окрестности сортируются по возрастанию яркости.

  3. Средняя точка отсортированной окрестности записывается в итоговое изображение.



Медианный фильтр квадратом 3х3



Один из способов оптимизации кода по обработке изображений на GPU — найти и выделить повторные операции и сделать их только один раз внутри одного warp'а. В некоторых случаях потребуется немного изменить алгоритм. Рассматривая фильтрацию, можно заметить, что, загружая квадрат 3х3 вокруг одного пикселя, соседние нити делают одинаковые логические операции и одинаковые загрузки из памяти:



__global__ __launch_bounds__(256, 8) 
void mFilter_3х3(const int H, const int W, const unsigned char * __restrict in, unsigned char *out)
{
int idx = threadIdx.x + blockIdx.x * blockDim.x + 1;
int idy = threadIdx.y + blockIdx.y * blockDim.y * 4 + 1;
unsigned a[9] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 };

if (idx < W - 1)
{
for (int z3 = 0; z3 < 4; ++z3)
{
const int shift = 8 * (3 - z3);
for (int z2 = 0; z2 < 3; ++z2)
for (int z1 = 0; z1 < 3; ++z1)
a[3 * z2 + z1] |= in[(idy - 1 + z2) * W + idx - 1 + z1] << shift; // <----- Здесь каждый warp делает лишние
// 6 операций загрузки и 12 логических
// операций
idy += BL_Y;
}

/* Остальной код далее */
}
}


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

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

Старая сортировка:





Новая сортировка:





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





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



__global__ BOUNDS(BL_X, 2048 / BL_X)
void mFilter_3x3(const int H, const int W, const unsigned char * __restrict in, unsigned char *out)
{
// индекс по Оси Х
int idx = threadIdx.x + blockIdx.x * blockDim.x;
// смещаемся немного назад, так как в каждом warp'е первая и последняя нить не участвует в итоговом вычислении
idx = idx - 2 * (idx / 32);

// умножаем на 4, так как мы обрабатываем по 4 точки сразу
int idy = threadIdx.y + blockIdx.y * blockDim.y * 4 + 1;
unsigned a[9];
unsigned RE[3];

// первая и последняя нить не участвует в итоговом вычислении
bool bound = (idxW > 0 && idxW < 31 && idx < W - 1);

// считываем и сортируем столбцы
RE[0] = in[(idy - 1) * W + idx] | (in[(idy + 0) * W + idx] << 8) | (in[(idy + 1) * W + idx] << 16) | (in[(idy + 2) * W + idx] << 24);
RE[1] = in[(idy + 0) * W + idx] | (in[(idy + 1) * W + idx] << 8) | (in[(idy + 2) * W + idx] << 16) | (in[(idy + 3) * W + idx] << 24);
RE[2] = in[(idy + 1) * W + idx] | (in[(idy + 2) * W + idx] << 8) | (in[(idy + 3) * W + idx] << 16) | (in[(idy + 4) * W + idx] << 24);

Sort(RE[0], RE[1]);
Sort(RE[1], RE[2]);
Sort(RE[0], RE[1]);

// делаем обмены
a[0] = __shfl_down(RE[0], 1);
a[1] = RE[0];
a[2] = __shfl_up(RE[0], 1);

a[3] = __shfl_down(RE[1], 1);
a[4] = RE[1];
a[5] = __shfl_up(RE[1], 1);

a[6] = __shfl_down(RE[2], 1);
a[7] = RE[2];
a[8] = __shfl_up(RE[2], 1);

// ищем максимум
a[2] = __vmaxu4(a[0], __vmaxu4(a[1], a[2]));

// ищем медиану
Sort(a[3], a[4]);
a[4] = __vmaxu4(__vminu4(a[4], a[5]), a[3]);

// ищем минимум
a[6] = __vminu4(a[6], __vminu4(a[7], a[8]));

// ищем финально медиану
Sort(a[2], a[4]);
a[4] = __vmaxu4(__vminu4(a[4], a[6]), a[2]);

// записываем результат, если не вышли за границу картинки
if (idy + 0 < H - 1 && bound)
out[(idy) * W + idx] = a[4] & 255;
if (idy + 1 < H - 1 && bound)
out[(idy + 1) * W + idx] = (a[4] >> 8) & 255;
if (idy + 2 < H - 1 && bound)
out[(idy + 2) * W + idx] = (a[4] >> 16) & 255;
if (idy + 3 < H - 1 && bound)
out[(idy + 3) * W + idx] = (a[4] >> 24) & 255;
}


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



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




  1. Обычно размеры картинок имеют фиксированные значения, например, популярные Full HD, Ultra HD, 720p. Можно просто иметь набор предкомпилированных ядер. Данная оптимизация дает порядка 10-15% к производительности.

  2. С версии Cuda Toolkit 7.5 появилась возможность выполнять динамическую компиляцию, которая позволяет выполнить подстановку значений во время выполнения.



Полностью оптимизированный код будет выглядеть так. Варьируя количество точек, можно получить максимальную производительность. В моем случае максимальная производительность была достигнута при numP_char = 3, то есть три набора по 4 точки или 12 точек на одну нить.



__global__ BOUNDS(BL_X, 2048 / BL_X)
void mFilter_3x3(const unsigned char * __restrict in, unsigned char *out)
{
int idx = threadIdx.x + blockIdx.x * blockDim.x;
idx = idx - 2 * (idx / 32);
int idxW = threadIdx.x % 32;

int idy = threadIdx.y + blockIdx.y * blockDim.y * (4 * numP_char) + 1;
unsigned a[numP_char][9];
unsigned RE[numP_char][3];

bool bound = (idxW > 0 && idxW < 31 && idx < W - 1);

#pragma unroll
for (int z = 0; z < numP_char; ++z)
{
RE[z][0] = in[(idy - 1 + z * 4) * W + idx] | (in[(idy - 1 + 1 + z * 4) * W + idx] << 8) | (in[(idy - 1 + 2 + z * 4) * W + idx] << 16) | (in[(idy - 1 + 3 + z * 4) * W + idx] << 24);
RE[z][1] = in[(idy + z * 4) * W + idx] | (in[(idy + 1 + z * 4) * W + idx] << 8) | (in[(idy + 2 + z * 4) * W + idx] << 16) | (in[(idy + 3 + z * 4) * W + idx] << 24);
RE[z][2] = in[(idy + 1 + z * 4) * W + idx] | (in[(idy + 1 + 1 + z * 4) * W + idx] << 8) | (in[(idy + 1 + 2 + z * 4) * W + idx] << 16) | (in[(idy + 1 + 3 + z * 4) * W + idx] << 24);

Sort(RE[z][0], RE[z][1]);
Sort(RE[z][1], RE[z][2]);
Sort(RE[z][0], RE[z][1]);

a[z][0] = __shfl_down(RE[z][0], 1);
a[z][1] = RE[z][0];
a[z][2] = __shfl_up(RE[z][0], 1);

a[z][3] = __shfl_down(RE[z][1], 1);
a[z][4] = RE[z][1];
a[z][5] = __shfl_up(RE[z][1], 1);

a[z][6] = __shfl_down(RE[z][2], 1);
a[z][7] = RE[z][2];
a[z][8] = __shfl_up(RE[z][2], 1);

a[z][2] = __vmaxu4(a[z][0], __vmaxu4(a[z][1], a[z][2]));
Sort(a[z][3], a[z][4]);
a[z][4] = __vmaxu4(__vminu4(a[z][4], a[z][5]), a[z][3]);
a[z][6] = __vminu4(a[z][6], __vminu4(a[z][7], a[z][8]));
Sort(a[z][2], a[z][4]);
a[z][4] = __vmaxu4(__vminu4(a[z][4], a[z][6]), a[z][2]);
}

#pragma unroll
for (int z = 0; z < numP_char; ++z)
{
if (idy + z * 4 < H - 1 && bound)
out[(idy + z * 4) * W + idx] = a[z][4] & 255;
if (idy + 1 + z * 4 < H - 1 && bound)
out[(idy + 1 + z * 4) * W + idx] = (a[z][4] >> 8) & 255;
if (idy + 2 + z * 4 < H - 1 && bound)
out[(idy + 2 + z * 4) * W + idx] = (a[z][4] >> 16) & 255;
if (idy + 3 + z * 4 < H - 1 && bound)
out[(idy + 3 + z * 4) * W + idx] = (a[z][4] >> 24) & 255;
}
}


Медианный фильтр квадратом 5х5



К сожалению, придумать какую то другую сортировку для квадрата 5х5 мне так и не удалось. Единственное, на чем можно сэкономить — на загрузке и объединении 4х точек в unsigned int. Приводить еще более длинный код я не вижу смысла, так как все преобразования можно проделать по аналогии с квадратом 3х3.

В данной статье описаны некоторые идеи по совмещению операций для двух квадратов с наложением в 20 элементов. Но предложенный авторами метод forgetful selection sort делает все равно больше операций, чем неполная сеть сортировки Бэтчера для 25 элементов, даже при объединении двух рядом расположенных квадратов 5х5.



Производительность



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

























































Время выполнения, ms Ускорение
Скалярный ЦПУ AVX2 ЦПУ GPU Скалярный / GPU AVX2 / GPU
3x3 1920x1080 22.9 0.255 0.044 (47.1 GP/s) 520 5.8
3x3 4096x2160 97.9 1.33 0.172 (51.4 GP/s) 569 7.73
5x5 1920x1080 134.3 1.47 0.260 (7.9 GP/s) 516 5.6
5x5 4096x2160 569.2 6.69 1.000 (8.8 GP/s) 569 6.69


Заключение



В заключении хочется отметить, что среди всех найденных статей и реализаций медианной фильтрации для GPU, интерес представляет уже упомянутая статья "Fine-tuned high-speed implementation of a GPU-based median filter", опубликованная в 2013 году. В данной статье авторы предложили совершенно другой подход к сортировке, а именно — метод forgetful selection. Суть данного метода заключается в том, что мы берем roundUp(N / 2) + 1 элементов, проходим слева направо и обратно, получая тем самым минимальный и максимальный элементы по краям. Далее забываем эти элементы, добавляем один из оставшихся элементов в конец и повторяем процесс. Когда добавлять уже будет нечего, мы получим массив из 3х элементов, из которых медиану выбрать будет уже просто. Один из плюсов данного подхода — уменьшенное количество используемых регистров, по сравнению с известными сортировками.



В статье указано, что авторы получили результат в 1.8 GPixels / sec на Tesla C2050. Мощность данной карты в одинарной точности оценивается в 1 TFlops. Мощность участвующей в тестировании 780Ti оценивается в 5 TFlops. Тем самым, удельная скорость вычислений на 1 TFlops предложенного мной алгоритма примерно в 5.5 раз больше для квадрата 3х3 и в 2 раза больше для квадрата 5х5, чем у предложенного авторами статьи. Данное сравнение является не совсем корректным, но более приближенным к истине. Также в данной статье было упомянуто, что на тот момент их реализация была самой быстрой из всех известных.



Достигнутое ускорение по сравнению с AVX2 версией составило в среднем в 6 раз. Если использовать новые карты на базе архитектуры Pascal, данное ускорение может увеличиться как минимум в 2 раза, что составит примерно 100 GPixels / sec.


Original source: habrahabr.ru (comments, light).

https://habrahabr.ru/post/308214/

Метки:   Комментарии (0)КомментироватьВ цитатник или сообщество
rss_rss_hh_new

Не Персеидами едиными или Моделируем вспышки спутников своими руками

Понедельник, 22 Августа 2016 г. 14:08 (ссылка)

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





Вспышка Иридиума, первое фото своими руками – навелся не туда, затвор открыл поздно, горизонт завалил :)



С чего все началось



Весной на Гиктаймсе появился пост про спутник «Маяк». Это кубсат размером 30х10х10 см, основная задача которого после выхода на орбиту – раскрыть зеркальную пленку, которая будет отражать солнечный свет в сторону Земли. По замыслу авторов, это должно сделать Маяк одним из ярчайших спутников в ночном небе. Пикантности добавляет тот факт, что Маяк – первый в мире спутник, средства на который собирали краудфандингом, причем весьма успешно.





(картинка отсюда)



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





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



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



Алгоритм



Основная задача модели – рассчитать максимальную яркость вспышки для заданных параметрах спутника и заданного положения наблюдателя на Земле. При этом предсказание вспышек на несколько дней уходит на второй план – для этого есть, например, heavens-above.com. Логично будет взять один пролет спутника над наблюдателем, и уже для него вычислить яркость отражения с небольшим шагом по времени. После этого, варьируя начальные параметры, можно будет вычислить максимально возможную яркость вспышки.



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




  • Первая функция считает координаты спутника и наблюдателя в течение одного витка с большим шагом по времени (1-5 секунд) и проверяет, в какие моменты спутник виден над горизонтом и не находится в тени. Если спутник появляется над горизонтом, функция рассчитывает интервалы времени, когда он виден, и передает их второй функции.

  • Вторая функция работает с маленьким шагом по времени (0.001 – 1 с) и на каждом шаге рассчитывает яркость вспышки, которую видит наблюдатель.

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



Физическая модель



Моделировать задачу будем в системе отсчета ECI (Earth-centered inertial) – ее начало совпадает с центром Земли, ось Z направлена на северный полюс, а Солнце находится в плоскости X-Z.





Почему именно ECI? Потому, что Солнце в ней неподвижно, наблюдатель вращается вместе с Землей вокруг оси Z, а орбита спутника не меняет своего положения в пространстве. Это значит, что координаты и спутника, и наблюдателя легко параметризуются через углы. Подробности – под спойлером.



Больше координат богу координат!

Координаты наблюдателя задаются широтой fiobs и текущей долготой thetaobs.





По сути thetaobs является местным временем: полдень соответствует thetaobs = 0, а полночь – thetaobs = 180. Еще нужно знать местное время в начальный момент времени thetaobs0 и скорость вращения Земли вокруг своей оси. Тогда координаты наблюдателя в ECI выражаются так:





Орбиту спутника для облегчения расчетов можно считать круговой. В этом случае ее можно задать тремя параметрами: радиусом Rorb, наклонением fi и долготой восходящего узла thetaAN. Восходящий узел (ascending node) – это точка, в которой орбита пересекает экватор по направлению к северу.





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





Теперь, когда мы знаем все это, положение спутника на орбите легко параметризовать через угол alpha:





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







Солнце. Фишка системы отсчета ECI в том, что Солнце всегда находится в плоскости X-Z. В зависимости от времени года оно меняет положение: когда в северном полушарии лето, Солнце находится выше оси Z, а зимой – наоборот. Нас интересует не само положение Солнца, а направление на него. Вернее, от него, то есть вдоль направления движения солнечных лучей. На картинках он называется ksun.



Посчитали координаты? Теперь надо проверить, виден ли спутник над горизонтом – иначе считать что-либо бессмысленно. Это определяется сонаправленностью векторов от наблюдателя на спутник -dirsat_obs и в зенит dirobs: если их скалярное произведение положительно, то спутник где-то над нами.



Если спутник над горизонтом, то нужно узнать его координаты на небе. Для этого направление на спутник надо спроецировать на плоскость наблюдателя. Плоскость наблюдателя задается векторами на север dirnord и запад dirwest. Первый находится через направление на зенит и северный полюс Земли dNP, а уже по нему и направлению в зенит вычисляется вектор на запад.





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




  • угол между ksun и dirsat (вектор на спутник) острый;

  • d < R





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



Наконец, если спутник освещен Солнцем и находится где-то над нами, остается посчитать яркость вспышки. Точнее, ее видимую звездную величину, которая показывает яркость по сравнению с другими звездами. Шкала звездных величин логарифмическая: она определяется как





то есть изменение яркости в 100 раз соответствует изменению звездной величины на 5 единиц. А еще отсюда видно, что шкала звездных величин инвертирована: более ярким объектам соответствуют меньшие значения. Например, звездная величина Арктура и Веги – ярчайших звезд Северного полушария – составляет +0.0, Сириуса – -1.7, а вот Венера достигает -4.6 звездной величины. Глаз различает звезды вплоть до +6 звездной величины, хотя в крупных городах засветка не дает увидеть что-то ярче +4.



Еще немного фотометрии

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



Как посчитать освещенность, которую создает отражение от спутника? Прежде всего, мощность отраженного света (то есть сколько фотонов за секунду отражаются от спутника) находится перемножением площади зеркала на освещенность, создаваемую Солнцем. Нужно не забыть, что зеркало может быть под углом к солнечным лучам. А еще про то, что зеркало неоднородно и отражает свет в большой телесный угол – «зайчик» на поверхности Земли будет протяженным, и мы можем быть вдали от его центра. Наконец, освещенность квадратично зависит от расстояния до спутника – чем дальше спутник, тем меньше света достигнет наблюдателя. Конечное выражение для освещенности выглядит так:





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





Сравнивая это число освещенностью, соответствующей нулевой звездной величине (2,5•10-8 Вт/м2), находим звездную величину спутника.



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





Легко и просто? Отнюдь! На практике быстрее считать, как должно располагаться зеркало, чтобы при заданном положении спутника и наблюдателя отражение падало бы прямо на наблюдателя. И затем сравнивать с текущим положением зеркала. Преимущество такого подхода в том, что он позволяет изменять шаг по времени. Пока зеркало далеко от оптимального положения, шаг можно увеличить – вспышки в ближайшее время все равно не будет. Если же вспышка уже началась, то шаг по времени нужно уменьшить, чтобы посчитать ее с максимальным разрешением.



Наконец, спутник вращается, а вместе с ним вращается и зеркало. Поэтому нормаль к зеркалу nmir удобно разложить на две компоненты – параллельную и перпендикулярную оси вращения спутника dmir. Первая остается неподвижной, вторая прецессирует – математика этого точь-в-точь такая же, как и движение спутника по орбите.



Код



Код написан на матлабе. Как я рассказал выше, он разбит на три функции, которые по очереди вызываются основным скриптом calc_one. В основном скрипте задаются константы, начальные условия и параметры симуляции, которые передаются функциям через структуру param следующим образом:



% satellite orbit
param.incl_sat = 97.75; % inclination, deg (ISS: 51.6, SSO for 600 km: 97.75)
param.theta_asc = -7.5; % longitude of the ascending node, deg
param.alpha_sat_0 = 0; % angle between AN and the sat in the plane of orbit at t=0, deg
param.h_orbit = 600 * 1E3; % orbit height, m


Еще больше параметров!
%%%%%%%%%%%%%%%%%%     SIMULATION PARAMETERS     %%%%%%%%%%%%%%%%%%%%%%%%%%

% world
param.rad_earth = 6400 * 1E3; % Earth radius, m
param.GM = 6.67*1E-11 * 5.972*1E24; % grav. constant times Earth mass, m^3/s^2
param.omega_obs = 1/(60*60*24); % Earth angular velocity, turns/s
param.k_sun = [-0.99, 0, 0.16]; % sunlight wavevector (22 Jun: [-0.92,0,-0.4]; 21 Dec: [-0.92,0,0.4], 15 Oct: [-0.99,0,0.16])
param.atmosphere_transmission = 0.7; % atmosphere transmission
param.flux_sun = 1367; % Sun irradiant flux on the Earth orbit, W/m^2
param.flux_m0 = 2.5*1E-8; % flux corresponding to the apparent magnitude = 0, W/m^2
param.min_m = 6; % minimal visible apparent magnitude

% observer
param.phi_obs = 55.75; % latitude, deg (Moscow: 55.75)
param.theta_obs_0 = 170; % longitude, deg (noon: 0, midnight: 180)

% satellite orbit
param.incl_sat = 97.75; % inclination, deg
param.theta_asc = -7.5; % longitude of the ascending node, deg
param.alpha_sat_0 = 0; % angle between AN and the sat in the plane of orbit at t=0, deg
param.h_orbit = 600 * 1E3; % orbit height, m

% satellite mirror
param.omega_rot = 1; % rotation speed of the satellite, turns/s
param.d_rot = [0 0.7 1]; % rotation axis of the satellite
param.n_mirror_0 = [1 0 0]; % normal to the mirror at t=0
param.gauss_w = 17.2; % divergence angle of a normal distrubution, deg
param.sq_mirror = 3.8; % mirror square, m^2

% simulation parameters
param.dt1 = 1; % time step if the satellite is below horizon or in the shadow, s
param.dt2 = 1; % same if satellite is above horizon, but reflection is far away
param.dt3 = 1 * 1E-2; % same if satellite is above horizon, reflection is close
param.dt4 = 5 * 1E-3; % same if reflection is seen
param.frame_duration = 0.03; % exposition time of a camera/eye, s


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



    % calculate some aux values
rad_orbit = rad_earth + h_orbit; % radius of the orbit from the center of Earth
r_np = [0 0 rad_earth]; % coordinates of the North pole
dir_an = [cosd(theta_asc), sind(theta_asc), 0]; % direction on the acsending node
dir_b = [-sind(theta_asc)*cosd(incl_sat),...
cosd(theta_asc)*cosd(incl_sat), sind(incl_sat)]; % direction on the highest point of the sat
omega_sat = sqrt(2*GM/rad_orbit)/(2*pi*rad_orbit); % satellite angular velocity on the orbit, s^-1
dir_obs_z = sind(phi_obs); % z-coordinate of the observer, does not change


После этого запускается while-цикл по времени с постоянным шагом, равным 1 секунде. Условие выхода из цикла – уход спутника за горизонт (или же два оборота спутника без единого появления над горизонтом, что нереально на низких орбитах). В цикле последовательно рассчитываются координаты наблюдателя и спутника:



        % coordinates of the observer
theta_obs = mod(theta_obs_0 + 360*omega_obs*current_time, 360);
dir_obs = [cosd(phi_obs)*cosd(theta_obs), cosd(phi_obs)*sind(theta_obs), dir_obs_z];
r_obs = rad_earth * dir_obs;

% coordinates of the satellite
theta_sat = mod(alpha_sat_0 + 360*omega_sat*current_time, 360);
dir_sat = cosd(theta_sat)*dir_an + sind(theta_sat)*dir_b;
r_sat = rad_orbit*dir_sat;

% vector from the satellite to the observer
r_sat_obs = r_obs-r_sat;
dir_sat_obs = fn_norm(r_sat_obs);


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



        % CHECK1: check if the satellite is above the horizon
if fn_dot(dir_obs, dir_sat_obs) < 0


            % basis on the skychart 
r_obs_np = r_np - r_obs; % vector from the observer to the North pole
dist_obs_np = sqrt(sum(r_obs_np.^2));
cos_angle_obs_np = fn_dot(r_obs_np, -dir_obs)/dist_obs_np;
dir_nord = r_obs_np + cos_angle_obs_np*dist_obs_np*dir_obs;
dir_nord = fn_norm(dir_nord); % vector to the nord
dir_west = fn_cross(dir_obs, dir_nord);
dir_west = fn_norm(dir_west); % vector to the west

% coordinates of the satellite on the sky chart
pass_path.alt(pass_step,1) = 90 - acosd(fn_dot(-dir_sat_obs, dir_obs));
map_proj_nord = fn_dot(-dir_sat_obs, dir_nord);
map_proj_west = fn_dot(-dir_sat_obs, dir_west);
pass_path.azimuth(pass_step,1) = mod(-atan2d(map_proj_west, map_proj_nord), 360);


и, наконец, проверяется, находится ли спутник в тени:



            % CHECK 2: check that the satellite is not in the shadow
cos_angle_sat_ksun = fn_dot(dir_sat,k_sun);
dist_sat_ksun = rad_orbit*sqrt(1-cos_angle_sat_ksun^2); % distance from the sat to the central axis of shadow
if (dist_sat_ksun > rad_earth)
% satellite is not in the shadow


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



Это очень скучно, поэтому убрано под спойлер
    % STEP 2: if the pass was above horizon,
% select the regions when sat is above horizon and not in shadow
pass_path.reflection_intervals = [];
if isempty(pass_path.time) == 0

flag_in_shadow = 1;
num_regions_not_in_shadow = 0;
for counter_step = 1:size(pass_path.time,1)

if pass_path.in_shadow(counter_step) == 0
% satellite is not in shadow
if flag_in_shadow == 0
% sat is still not in shadow
else
% sat was in shadow and just went out of it
flag_in_shadow = 0;
num_regions_not_in_shadow = num_regions_not_in_shadow + 1;
pass_path.reflection_intervals(num_regions_not_in_shadow,1) = pass_path.time(counter_step);
end
else
% satellite is in shadow
if flag_in_shadow == 0
% sat was not in shadow and just went into it
flag_in_shadow = 1;
pass_path.reflection_intervals(num_regions_not_in_shadow,2) = pass_path.time(counter_step);
else
% sat is still in shadow
end
end

end

% if at the end sat was still not in shadow,
% time of the end of the last interval will be time of the last
% point above the horizon
if flag_in_shadow == 0
pass_path.reflection_intervals(num_regions_not_in_shadow,2) = pass_path.time(counter_step);
end

end


Интервалы сохраняются в pass_path.reflection_intervals в массиве Nx2, где N – количество интервалов видимости, 2 числа – время начала и конца интервала видимости. Вот теперь структура pass_path возвращается в главный скрипт, откуда вызывается вторая функция.



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



            % CHECK 3: if angle (d_rot, n_opt) is close to angle (d_rot, n_mir0), proceed
% this means that geometry can in principle give a reflection
if (angle_drot_nopt > (angle_drot_nmirror-angle_div))...
&& (angle_drot_nopt < (angle_drot_nmirror+angle_div))

% current orientation of the mirror
n_mir = n_mir_parallel...
+ n_mir_perp0*cosd(360*omega_rot*current_time)...
+ n_mir_perp1*sind(360*omega_rot*current_time);

% angle between current orientation of the mirror and n_opt
angle_nmir_nopt = acosd(fn_dot(n_mir, n_opt));

% CHECK 4: if n_mir is close to n_opt, proceed and calculate the reflection
if angle_nmir_nopt < angle_div


Маленькая военная хитрость про положение зеркала

Как видно из кода, расположение зеркала проверяется по двум параметрам – азимутальному (CHECK 3) и радиальному (CHECK 4) углам по отношению к оси вращения спутника. Фишка в том, что радиальный угол меняется при вращении спутника. А вот азимутальный – нет!



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



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



                    % calculate the angle between the direction on the observer
% and direction of the main reflection
cos_angle_nmir_ksun = fn_dot(n_mir, -k_sun);
k_refl = 2*cos_angle_nmir_ksun*n_mir + k_sun;
cos_angle_refl_obs = abs(fn_dot(k_refl, dir_sat_obs));
angle_refl_obs = acosd(cos_angle_refl_obs);

% flux at the observer's place
refl_fraction = 2/(pi*degtorad(gauss_w)^2) * exp (-2*angle_refl_obs^2/gauss_w^2); % now it is Gaussian
int_reflection = flux_sun * sq_mirror * cos_angle_nmir_ksun;
flux = max (int_reflection * refl_fraction ...
* atmosphere_transmission / distance^2, flux_min); % W/m^2


Видно, что если яркость слишком мала (тусклее +6), то в результат запишется +6. Это значит, что мы ничего не увидим. Все результаты записываются в структуру pass, которую функция и возвращает.



Третья функция fn_flares самая простая: она считает вспышки, их максимальную яркость, длительность и расположение на небе. Она реально примитивная, но длинная, поэтому я убрал ее под спойлер.



fn_flares
function pass = fn_flares(pass, param)

% this functions counts the flares into the 'pass' structure
% summarizes data about their magnitudes, durations, etc.
% and saves them in 'pass.flares' field

% initialization before the main loop
counter_flare = 0;
flare_is_now = 0;
clear flares
flares = struct('time',[],'dist',[],'alt',[],'azimuth',[],'inst_magn',[],'vis_magn',[],'dur',[]);

% main loop
for step = 1:length(pass.time)

if pass.magnitude(step) < param.min_m
% flare is present now
if flare_is_now == 0

% flare just began
flare_is_now = 1;
counter_flare = counter_flare + 1;
clear current_flare
current_flare.time_beg = pass.time(step);
% integrate the energy of the flare, [J/m^2]
if step < length(pass.time)
current_flare.energy = pass.flux(step)*...
(pass.time(step+1) - pass.time(step));
else
% this was the last data point, nothing to add
end

% create field 'current_flare.brightest'
% with the data about the brightest point of the flare
current_flare.brightest.magn = pass.magnitude(step);
current_flare.brightest.time = pass.time(step);
current_flare.brightest.azimuth = pass.azimuth(step);
current_flare.brightest.alt = pass.alt(step);
current_flare.brightest.distance = pass.distance(step);
else
% flare still goes on
% increment the energy of the flare
if step < length(pass.time)
current_flare.energy = current_flare.energy + pass.flux(step)*...
(pass.time(step+1) - pass.time(step));
else
% this was the last data point, nothing to add
end
% if here flare is brighter, renew the data about the brightest point
if pass.magnitude(step) < current_flare.brightest.magn
current_flare.brightest.magn = pass.magnitude(step);
current_flare.brightest.time = pass.time(step);
current_flare.brightest.azimuth = pass.azimuth(step);
current_flare.brightest.alt = pass.alt(step);
current_flare.brightest.distance = pass.distance(step);
end
end
else
% no flare now
if flare_is_now == 0

% still no flare
else
% flare just ended, sum up the results
flare_is_now = 0;
current_flare.time_end = pass.time(step-1);
current_flare.duration = current_flare.time_end - current_flare.time_beg;
if current_flare.duration > param.frame_duration
% if the flare is long, we can see it directly
current_flare.vis_magn = current_flare.brightest.magn;
else
% if the flare is short, one should divide its energy
% over the exposure time of an eye/camera to get the
% apparent magnitude
current_flare.vis_magn = ...
min(-2.5*log10(current_flare.energy/param.frame_duration/param.flux_m0), param.min_m);
end

% save the data in the 'flares' structure
flares.time(counter_flare,1) = current_flare.brightest.time;
flares.dist(counter_flare,1) = current_flare.brightest.distance;
flares.azimuth(counter_flare,1) = current_flare.brightest.azimuth;
flares.alt(counter_flare,1) = current_flare.brightest.alt;
flares.inst_magn(counter_flare,1) = current_flare.brightest.magn;
flares.vis_magn(counter_flare,1) = current_flare.vis_magn;
flares.dur(counter_flare,1) = current_flare.duration;
end
end

end

% return the result
flares.num_of_flares = counter_flare;
pass.flares = flares;

end


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



Что значит «короткая вспышка»? Очевидно, такая, которая короче времени экспозиции камеры. А если мы наблюдаем ее глазом? Что считать экспозицией глаза, я не имел ни малейшего понятия, поэтому спросил у Meklon. Он рассказал много всего интересного – но сводилось это к тому, что глаз работает ооочень сложно, чувствительность у него сильно нелинейная, и поэтому четкого ответа на вопрос дать нельзя. Даже примерно. В итоге было решено выбрать какое-нибудь более-менее осмысленное число — я выбрал 30 мс, что соответствует примерно 30 кадрам в секунду.



Итак, если вспышка оказывается короче времени экспозиции, то ее длительность округляется вверх до времени экспозиции, суммарное количество света остается неизменным, ну а яркость, соответственно, падает. Например, если в реальности вспышка длится 10 мс и за это время прилетело 3000 фотонов, то ее моментальная «яркость» составляет 300 фотонов/мс. Мы же увидим, что она продолжалась 30 мс с «яркостью» 100 фотонов/мс – в 3 раза меньше моментальной.



Наконец, главный скрипт calc_one сначала вызывает первую функцию fn_pass, и если спутник виден над горизонтом и не находится в тени, то считает яркость отражений функцией fn_reflection и параметры вспышек функцией fn_flares:



% calculate trajectory of the pass
pass_path = fn_pass(param);

% proceed if the satellite is above the horizon and not always in shadow
if isempty(pass_path.reflection_intervals) == 0
% calculate the reflections and magnitudes of flares
pass = fn_reflection(param, pass_path.reflection_intervals);
pass = fn_flares(pass, param);
end


Оптимизация



Первый прием оптимизации я уже упомянул: это переменный шаг по времени. Функция fn_pass работает с большим шагом (1 с) – для проверки того, виден ли спутник над горизонтом, этого более, чем достаточно. А вот у функции fn_reflection целых три возможных шага по времени: самый маленький (до 1 мс) нужен для расчета профиля вспышки, два других – когда вспышка не видна, но вскоре может появиться. Они определяются в структуре с параметрами:



param.dt1 = 1;                          % time step if the satellite is below horizon or in the shadow, s
param.dt2 = 1; % same if satellite is above horizon, but reflection is far away
param.dt3 = 1 * 1E-2; % same if satellite is above horizon, reflection is close
param.dt4 = 5 * 1E-3; % same if reflection is seen


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



Встроенная функция для скалярного произведения
function c = dot(a,b,dim)

if isinteger(a) || isinteger(b)
error(message('MATLAB:dot:integerClass'));
end

% Special case: A and B are vectors and dim not supplied
if ismatrix(a) && ismatrix(b) && nargin<3
if min(size(a))==1, a = a(:); end
if min(size(b))==1, b = b(:); end
end;

% Check dimensions
if any(size(a)~=size(b)),
error(message('MATLAB:dot:InputSizeMismatch'));
end

if nargin==2,
c = sum(conj(a).*b);
else
c = sum(conj(a).*b,dim);
end


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



function c = fn_dot(a, b)
c = a(1)*b(1) + a(2)*b(2) + a(3)*b(3);
end


После этого скорость вычислений выросла примерно вдвое.



Лирическое отступление: анонимные функции в MATLAB.

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



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



Код можно найти на гитхабе. Расчет вспышек проверялся на Иридиуме и давал осмысленные результаты – яркость вспышек составляла около -9 звездной величины.



Маяк



Спутник «Маяк» планируется запустить в этом октябре на круговую орбиту высотой 600 км. В средней полосе России его будет лучше всего видно примерно с 22:00 до 02:00 по местному времени. Пленочные зеркала Маяка отражают свет гораздо хуже антенн Иридиума, зато заметно превосходят их по площади.



Параметры Маяка

























Параметр Значение
Высота орбиты 400 – 600 км
Наклонение орбиты 97.7°
Направление на восходящий узел -7.5°
Площадь зеркала 3.8 м2
Расходимость отражения (1

https://habrahabr.ru/post/307212/?utm_source=habrahabr&utm_medium=rss&utm_campaign=best

Метки:   Комментарии (0)КомментироватьВ цитатник или сообщество
rss_rss_hh_new

[Из песочницы] Обучение с подкреплением для самых маленьких

Пятница, 19 Августа 2016 г. 10:09 (ссылка)

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



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



Задача



С помощью метода «обучение с подкреплением» необходимо научить тележку отъезжать от стены на максимальное расстояние. Награда представлена в виде значения изменения расстояния от стены до тележки при движении. Измерение расстояния D от стены производится дальномером. Движение в данном примере возможно только при определенном смещении «привода», состоящего из двух стрел S1 и S2. Стрелы представляют собой два сервопривода с направляющими, соединенными в виде «колена». Каждый сервопривод в данном примере может поворачиваться на 6 одинаковых углов. Модель имеет возможность совершить 4 действия, которые представляют собой управление двумя сервоприводами, действие 0 и 1 поворачивают первый сервопривод на определенный угол по часовой и против часовой стрелке, действие 2 и 3 поворачивают второй сервопривод на определенный угол по часовой и против часовой стрелке. На рисунке 1 показан рабочий прототип тележки.





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



На рисунке 2 красным цветом выделена стрела S2, синим цветом – стрела S1, черным цветом – 2 сервопривода.





Рис. 2. Двигатель системы



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





Рис. 3. Схема системы



Диапазон возможных положений для S1 и S2 показан на рисунке 4:





Рис. 4.а. Диапазон положений стрелы S1





Рис. 4.б. Диапазон положений стрелы S2



Пограничные положения привода показаны на рисунке 5:



При S1 = S2 = 5 максимальная дальность от земли.

При S1 = S2 = 0 минимальная дальность до земли.





Рис. 5. Пограничные положения стрел S1 и S2



У «привода» 4 степени свободы. Действие (action) изменяет положение стрел S1 и S2 в пространстве по определённому принципу. Виды действий показаны на рисунке 6.





Рис. 6. Виды действий (Action) в системе



Действие 0 увеличивает значение S1. Действие 1 уменьшает значение S1.

Действие 2 увеличивает значение S2. Действие 3 уменьшает значение S2.



Движение



В нашей задаче тележка приводится в движение всего в 2х случаях:

В положении S1 =0, S2 = 1 действие 3 приводит в движение тележку от стены, система получает положительное вознаграждение, равное изменению расстояния до стены. В нашем примере вознаграждение равно 1.





Рис. 7. Движение системы с положительным вознаграждением



В положении S1 = 0, S2 = 0 действие 2 приводит в движение тележку к стене, система получает отрицательное вознаграждение, равное изменению расстояния до стены. В нашем примере вознаграждение равно -1.





Рис. 8. Движение системы с отрицательным вознаграждением



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

Хочется отметить, что стабильным динамическим состоянием системы будет последовательность действий 0-2-1-3 из состояния S1=S2=0, в котором тележка будет двигаться в положительном направлении при минимальном количестве затраченных действий. Подняли колено – разогнули колено – опустили колено – согнули колено = тележка сдвинулась вперед, повтор. Таким образом, с помощью метода машинного обучения необходимо найти такое состояние системы, такую определенную последовательность действий, награда за которые будет получена не сразу (действия 0-2-1 – награда за которые равна 0, но которые необходимы для получения 1 за последующее действие 3).



Метод Q-Learning



Основой метода Q-Learning является матрица весов состояния системы. Матрица Q представляет собой совокупность всевозможных состояний системы и весов реакции системы на различные действия.

В данной задаче возможных комбинаций параметров системы 36 = 6^2. В каждом из 36 состояний системы возможно произвести 4 различных действия (Action = 0,1,2,3).

На рисунке 9 показано первоначальное состояние матрицы Q. Нулевая колонка содержит индекс строки, первая строка – значение S1, вторая – значение S2, последние 4 колонки равны весам при действиях равных 0, 1, 2 и 3. Каждая строка представляет собой уникальное состояние системы.

При инициализации таблицы все значения весов приравняем 10.





Рис. 9. Инициализация матрицы Q



После обучения модели (~15000 итераций) матрица Q имеет вид, показанный на рисунке 10.





Рис. 10. Матрица Q после 15000 итераций обучения



Обратите внимание, действия с весами, равными 10, невозможны в системе, поэтому значение весов не изменилось. Например, в крайнем положении при S1=S2=0 нельзя выполнить действие 1 и 3, так как это ограничение физической среды. Эти пограничные действия запрещены в нашей модели, поэтому 10тки алгоритм не использует.



Рассмотрим результат работы алгоритма:



Iteration: 14991, was: S1=0 S2=0, action= 0, now: S1=1 S2=0, prize: 0

Iteration: 14992, was: S1=1 S2=0, action= 2, now: S1=1 S2=1, prize: 0

Iteration: 14993, was: S1=1 S2=1, action= 1, now: S1=0 S2=1, prize: 0

Iteration: 14994, was: S1=0 S2=1, action= 3, now: S1=0 S2=0, prize: 1

Iteration: 14995, was: S1=0 S2=0, action= 0, now: S1=1 S2=0, prize: 0

Iteration: 14996, was: S1=1 S2=0, action= 2, now: S1=1 S2=1, prize: 0

Iteration: 14997, was: S1=1 S2=1, action= 1, now: S1=0 S2=1, prize: 0

Iteration: 14998, was: S1=0 S2=1, action= 3, now: S1=0 S2=0, prize: 1

Iteration: 14999, was: S1=0 S2=0, action= 0, now: S1=1 S2=0, prize: 0



Рассмотрим подробнее:

Возьмем итерацию 14991 в качестве текущего состояния.

1. Текущее состояние системы S1=S2=0, этому состоянию соответствует строка с индексом 0. Наибольшим значением является 0.617 (значения равные 10 игнорируем, описано выше), оно соответствует Action = 0. Значит, согласно матрице Q при состоянии системы S1=S2=0 мы производим действие 0. Действие 0 увеличивает значение угла поворота сервопривода S1 (S1 = 1).

2. Следующему состоянию S1=1, S2=0 соответствует строка с индексом 6. Максимальное значение веса соответствует Action = 2. Производим действие 2 – увеличение S2 (S2 = 1).

3. Следующему состоянию S1=1, S2=1 соответствует строка с индексом 7. Максимальное значение веса соответствует Action = 1. Производим действие 1 – уменьшение S1 (S1 = 0).

4. Следующему состоянию S1=0, S2=1 соответствует строка с индексом 1. Максимальное значение веса соответствует Action = 3. Производим действие 3 – уменьшение S2 (S2 = 0).

5. В итоге вернулись в состояние S1=S2=0 и заработали 1 очко вознаграждения.



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





Рис. 11.а. Матрица Q





Рис. 11.б. Матрица Q



Рассмотрим подробнее процесс обучения.



Алгоритм Q-learning
minus = 0;
plus = 0;
initializeQ();
for t in range(1,15000):
epsilon = math.exp(-float(t)/explorationConst);
s01 = s1; s02 = s2
current_action = getAction();
setSPrime(current_action);
setPhysicalState(current_action);
r = getDeltaDistanceRolled();
lookAheadValue = getLookAhead();
sample = r + gamma*lookAheadValue;
if t > 14900:
print 'Time: %(0)d, was: %(1)d %(2)d, action: %(3)d, now: %(4)d %(5)d, prize: %(6)d ' % \
{"0": t, "1": s01, "2": s02, "3": current_action, "4": s1, "5": s2, "6": r}
Q.iloc[s, current_action] = Q.iloc[s, current_action] + alpha*(sample - Q.iloc[s, current_action] ) ;
s = sPrime;
if deltaDistance == 1:
plus += 1;
if deltaDistance == -1:
minus += 1;
print( minus, plus )






Полный код на .



Установим начальное положение колена в крайнее верхнее положение:



s1=s2=5.


Инициализируем матрицу Q, заполнив начальным значением:



initializeQ();


Вычислим параметр epsilon. Это вес «случайности» действия алгоритма в нашем расчёте. Чем больше итераций обучения прошло, тем меньше случайных значений действий будут выбраны:



epsilon = math.exp(-float(t)/explorationConst)


Для первой итерации:



epsilon = 0.996672


Сохраним текущее состояние:



s01 = s1; s02 = s2


Получим «лучшее» значение действия:



current_action = getAction();


Рассмотрим функцию поподробнее.



Функция getAction() выдает значение действия, которому соответствует максимальный вес при текущем состоянии системы. Берется текущее состояние системы в матрице Q и выбирается действие, которому соответствует максимальный вес. Обратим внимание, что в этой функции реализован механизм выбора случайного действия. С увеличением числа итераций случайный выбор действия уменьшается. Это сделано, для того, чтобы алгоритм не зацикливался на первых найденных вариантах и мог пойти по другому пути, который может оказаться лучше.



В исходном начальном положении стрел возможны только два действия 1 и 3. Алгоритм выбрал действие 1.

Далее определим номер строки в матрице Q для следующего состояние системы, в которое система перейдет после выполнения действия, которое мы получили в предыдущем шаге.



setSPrime(current_action);


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

Выполним функции:



setPhysicalState(current_action);
— моделируем реакцию среды на выбранное нами действие. Изменяем положение сервоприводов, смещаем тележку.



r = getDeltaDistanceRolled();
— Вычисляем вознаграждение – расстояние, пройденное тележкой.



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

Теперь самое интересное – для расчета веса текущего шага заглянем в будущее.

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



lookAheadValue = getLookAhead();


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



sample = r + gamma*lookAheadValue;
sample = 7.5
Q.iloc[s, current_action] = Q.iloc[s, current_action] + alpha*(sample - Q.iloc[s, current_action] ) ;
Q.iloc[s, current_action] = 9.75


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

Этот простой трюк дает хорошие результаты сходимости алгоритма.



Масштабирование алгоритма



Данный алгоритм можно расширить на большее количество степеней свободы системы (s_features), и большее количество значений, которые принимает степень свободы (s_states), но в небольших пределах. Достаточно быстро матрица Q займет всю оперативную память. Ниже пример кода построения сводной матрицы состояний и весов модели. При количестве «стрел» s_features = 5 и количестве различных положений стрелы s_states = 10 матрица Q имеет размеры (100000, 9).



Увеличение степеней свободы системы
import numpy as np

s_features = 5
s_states = 10
numActions = 4

data = np.empty((s_states**s_features, s_features + numActions), dtype='int')
for h in range(0, s_features):
k = 0
N = s_states**(s_features-1-1*h)
for q in range(0, s_states**h):
for i in range(0, s_states):
for j in range(0, N):
data[k, h] = i
k += 1

for i in range(s_states**s_features):
for j in range(numActions):
data[i,j+s_features] = 10.0;

data.shape

# (100000L, 9L)




Вывод



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



Спасибо за внимание!
Original source: habrahabr.ru.

https://habrahabr.ru/post/308094/?utm_source=habrahabr&utm_medium=rss&utm_campaign=best

Метки:   Комментарии (0)КомментироватьВ цитатник или сообщество
rss_rss_hh_new

Рекомендательные системы в онлайн-образовании. Продолжение

Вторник, 16 Августа 2016 г. 09:40 (ссылка)

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



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



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







Общие слова



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



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



Рекомендательная система может применяться в разных ситуациях. Два основных случая:




  • обычная рекомендация для конкретного пользователя, которая показывается ему на главной странице и во вкладке ”Рекомендованные уроки”, а также время от времени присылается на почту;

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



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



Пользовательские данные



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



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



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



Графы переходов



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



У каждого пользовательского действия в логе есть временна

https://habrahabr.ru/post/307670/?utm_source=habrahabr&utm_medium=rss&utm_campaign=best

Комментарии (0)КомментироватьВ цитатник или сообщество

Следующие 30  »

<алгоритмы - Самое интересное в блогах

Страницы: [1] 2 3 ..
.. 10

LiveInternet.Ru Ссылки: на главную|почта|знакомства|одноклассники|фото|открытки|тесты|чат
О проекте: помощь|контакты|разместить рекламу|версия для pda