Немного о тетраэдальных сетках. Разбиение границы. Минусы тетраэдальной сетки.

Фактически, этот пост - продолжение этого.
Остановились мы на том, что построили тетраэдальную сетку из кубической. Следующая стадия - построить разбиение границы. Главная проблема в том, что разбиение границы должно быть согласованно с разбиением самой области - то есть, если у нас есть треугольник на границе, то этот треугольник полностью принадлежит одному из тетраэдров. Естественно есть множество способов сделать это. Простейший - берём прямоугольник разбиваем его на 4 треугольника (по 3 вершины разными способами), далее пробегаем по всем тетраэдрам и проверяем каждый тетраэдр на наличие этого треугольника. Проблема очевидна - на больших сетках это будет занимать гигантское количество времени.
Теперь если подумать о более рациональных способах, то можно сделать так: хранить отдельно параллелепипеды, находящиеся у границы сетки. Тогда сначала можно будет сделать разбиение внутренних, а затем внешних вместе с границей - тогда пробегать по тетраэдрам не придётся вообще(ну даже если придётся, то для каждого всего по 6 тетраэдрам, а не по всем, что есть - уже существенный плюс). Ну либо наоборот, порядок не важен.

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

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

Немного о тетраэдальный сетках. Алгоритмы построения.

С начала кратко, о чём будет пост: о построение тетраэдальных сеток (сеток, элементы которых - тетраэдры), получение этих сеток из сеток на треугольных призмах и параллелепипедальных сеток (далее буду называть их кубическими -  это короче и проще выговаривать и читать).
Теперь то, что привело меня к этому: для курсового проекта требовалось строить тетраэдальные сетки, причём большие (в итоге получилось 132651 узлов и 750000 тетраэдров), что руками делать не очень быстро. Есть различные методы построения, но я выбрал достаточно простой (форма области была не важна) - строится кубическая сетка, затем она разбивается на тетраэдры. И столкнулся с некоторыми проблемами - так, как я сначала  разбивал сетку - она получала не комфорной, ну и первом делом я полез в гугл и, что самое печально я не нашёл того, чего искал. Возможно я конечно плохо искал, но всё равно. Поэтому путь это будет хотя бы тут.
Ну чтож, а теперь начнём с краткого обзора алгоритмов построения тетраэдальных сеток(и соответсвено треугольных, если смотреть двухмерное разбиение).
1. Алгоритм исчерпывания
Выходные данные: граничный фронт (множество треугольников - разбиение границы области)
Принцип работы алгоритма достаточно простой и понятный, берём треугольник из фронта, откладываем по нормали к нему точку на определённую длину, если нет пересечений с тетраэдрами то добавляем тетраэдр в разбиение, если есть пытаемся модифицировать его(переносим точку) и добавляем в разбиение. Исходный треугольник удаляем из фронта, полученные(из тетраэдра) добавляем в него. Заканчиваем когда больше не можем строить тетраэдры.
Плюсы  и минусы: явным минусом является необходимость начального фронта - его тоже ещё надо построить, а так же нахождение оптимальной для задачи стратегии выбора шага по нормали и модификации  тетраэдров. Плюсам является то, что не важна геометрия области.

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

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

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

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

 Из неё получаются следующие три тетраэдра: (2, 4, 5, 6), (1, 2, 3, 6) и (1, 2, 4, 6).
Теперь перейдём к параллелепипедам: по суть мы разбиваем один параллелепипед на две призмы, а эти призмы на тетраэдры. Главная проблема - это правильная локальная нумерация призм друг по отношению к другу. Возьмём эту призму и добавим к ней ещё два узла, тем самым получим парллелепипед:

На кривости рисунка из-за копипаста обращать внимание не будем. Если взять симметричную нумерацию(7->2, 8->5), как с делал изначально сетка получится не комфорная, что в некоторых(даже думаю в большинстве) случаев будет проблемой. Правильно же локальную нумерацию во второй(левой) призме ввести следующим образом(глоб. -> лок.): 3 -> 1, 1 -> 2, 7 -> 3, 6 -> 4, 4 -> 5, 8 -> 6. И тогда, используя предыдущее разбиение мы получим наши заветные 6 тетраэдров: (2, 4, 5, 6), (1, 2, 3, 6), (1, 2, 4, 6), (1, 6, 4, 8), (3, 1, 7 8) и (3, 1, 6, 8).

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

Visual Studio: небольшие странности с dll

Без большого объявления подробностей, надо было сделать две вещи:

1. Использовать dll
2. Компилировать с помощью bat-файла

 Экспортировалась из dll функция с простым прототипом:

extern "C" __declspec(dllexport) unsigned int doit(string& dir)

Возможно конечно я не правильно что-то делал(не extern "C", например), но если я компилировал и dll и использующую его программу в IDE любой из 2008, 2010 и 2012 Visual Studio всё работало прекрасно. Если же компиляция производилась, с помощью cl, то приложение вылетало. Два дня я запускал приложение и компиляцию на различных компьютерах и в итоге выяснил интересную вещь: если компилировать в IDE всё работает, если компилировать из-под консоли на Windows 7 работает только в случае, если используется VS 2012, в случае 2008 и 2010 - приложение вылетает. На Windows XP с точностью наоборот - на 2008 и 2010 (а даже 2005) работает, на 2012 - вылетает. Как только компилировать не пробовал. В итоге отчаялся и решил поэкспериментировать с самой прогой и изменил портотип на:

extern "C" __declspec(dllexport) unsigned int doit(char* windir)

И всё заработало! И в IDE и из-под консоли. Я не знаю из-за чего возникли эти странности, но вот так вот оно получилось.
Будьте аккуратнее с тем, что экспортируете.

Случайный поиск или на удивление быстрый поиск

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

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

Собственно алгоритм простой:

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

Алгоритм очень простой, теперь поговорим про две реализации, первая - которая приходит первой на ум, вторая - немного оптимизированный вариант. На счёт генерации - стандартный генератор C/С++ примерно воспроизводит равномерное распределение, поэтому использовать можно его.

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

Вторая реализация
Поскольку точки могут генерироваться независимости друг от друга, а так же вычисляться значение функции, то первое что пришло на ум - распараллелить. Осталось это сделать так, чтобы прога не занимала слишком много памяти(3 миллиарда double - это всё-таки уже много, если я не ошибся 22.4 Гб, если размер типа 8 байт). Поэтому реализовывалось всё так:
У нас есть n потоков, каждый поток должен посчитать по m точек, за раз он считает N точек(тоже параллельно), выбирает из них минимум и сравнивает с текущим приближением, далее, если надо меняет текущее приближение и генерирует ещё N точек и так пока не сгенериует их m штук. В конце у нас остаётся n чисел - минимум из каждого потока, их них выбирается уже глобальный минимум.
Число N можно высчитывать как-нибудь из сходя из общего числа точек на поток, но я просто бал фиксированно N = 10000 (были и другие варианты). Потому что меньших N (при N = 100) затраты на создание потоков были больше чем выгода.
Так что если нам надо хранить три числа(координаты точки и значение функции), то уйдёт всего 3*(n + n*N) ячеек памяти, что намного меньше чем 3*n*m (учитываем, что число потоков n на персональном компьютере невелико - от 2 до 16, примерно, а число m может быть несколько сотен миллионов, то выигрыш в памяти, по сравнению в всеобщим распараллеливанием виден).
Теперь зачем нужно делать это порциями по N штук, а не посчитать по m на потоке - по тем же самым причинам время и память, в рассматриваемом случае минимум из N ищется очень быстро(если ещё упорядочивать при добавлении то ещё быстрее) и памяти надо дополнительной, под 3*N, а не 3*m.
Плюсы: скорость - вот тут и было моё удивление, меньше за минуту на 4-х потоках обрабатывает до миллиарда точек.
Минусы: в реализации конечно не сильно труден, но значительно труднее чем первый метод.

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

Оценка времени работы алгоритмов. Виды роста затрат.

В последнее время было лень что-то писать. Очень лень. Но я собрался и решил написать о довольно лёгкой и в тоже время нужной теме - временные оценки для алгоритмов. Их цель - примерно оценить алгоритм по времени работы, в зависимости от входных данных и дать возможность сравнить алгоритмы, не имея их реализации. В прицепе, у всех алгоритмов есть особенности и если они[алгоритмы] легко и быстро реализуемы, то в конкретном случае сравнение на "практике" будет намного информативнее.
Об обозначениях: обычно используется O-натация, то есть обозначается как О(f(x)) или o(f(x)). Первое можно интерпретировать, как c*f(x), c>1, во втором случае 0<c<=1. Сама функция f(x) может быть, как непрерывной. так и дискретной, под x здесь понимается размер входных данных. Например x=n - длина массива, то функция будет дискретной, x=l - длина отрезка, функция будет непрерывной.
И теперь к самой теме: характеры роста, то есть вид функции f(x). Перечислять их буду в порядке возрастания временных затрат: сначала самые быстрые, затем все более медленные.

1. Полиномиальный рост
• f(x) = log(x) - логарифмический рост, считается одной из лучших оценок. Поскольку log(x) ~ x^a, 0<a<1 при больших x, то этот рост можно назвать минимальным из полиномиальных. Поскольку оценка определяется с точностью умножения на константу, то основание логарифма не важно.
• f(x) = x^a, a >1 - общий вид полиномиального роста, является приемлемой оценкой. Чем меньше показатель степени a, тем оценка лучше - времени меньше. Как было указанно выше  рост вида x*log(x) можно привести к этому типу роста
2. Экспоненциальный рост
Здесь небольшая путаница в терминологии,  при оценки алгоритмов часто любой не полиномиальный рост называют экспоненциальным.

• f(x) = e^x - экспоненциальный рост, рост в геометрической прогрессии. Является не очень хорошей оценкой для работы, однако обычно не критической. Если есть полиномиальный алгоритм, решающий эту задачу - лучше выбирать его.
• f(n) = n!, f(x) = Г(x) - факториальный рост, его спокойно можно включать в следующую группу, однако такие оценки тоже часто встречаются. Поясняю, почему здесь две функции - факториал достаточно известная оценка, но она определенна только для дискретных величин, для непрерывных общение факториала - гамма-функция.
3. Алгоритмы, от которых лучше бежать подальше
Врятле такие оценки вы когда-нибудь встретите, они здесь больше для сравнения - может быть и хуже.

• f(n) = sf(n) - суперфакториальный рост. sf(n) = 1! * 2! * .. * n! - произведение факториалов.  Растёт очень быстро.
• f(n) = (n!)^(n!) - я не думаю, что у такой зависимости есть название вообще. При увеличении n на единицу порядок порядка n увеличивается на единицу. Убунтовский калькулятор (9!)^(9!) на нетбуке уже считает оооочень долго(за 10 минут не посчитал), а виндосовский уже (7!)^(7!) не считает - переполнение.

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

Добавил (n!)^(n!) на график, только при n=1,2,3,4. Для сравнения: с (4!)^(4!) смогло сравняться только 32!. Ну а что дальше - понятно.

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

Стандартный random C++

Что-то вдруг захотелось посмотреть, как работает функция rand() в C++. Не алгоритм генерации, а распределение результатов. Ну и забавы ради посмотрел что получилось.
Результаты получились довольной скучные - распределение примерно равномерное. Как проводились тесты - создаётся массив, зануляется, затем 10000000 прибавляется единица к случайному элементу массива, ну и затем выводится относительная частота. Объяснение немного странное, поэтому в конце поста есть код программы.

1. Без перемешивания, 32768 элементов.
 
 Под перемешиванием я имею ввиду функцию srand





 Первые 100 элементов из этой последовательности:














Перемешивание в начале генерации послеовательности, 32768 элементов
 
 Аналогично первые 100:











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

Код программы(для последнего случая):

#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <sstream>
#include <time.h>

using namespace std;

string convertInt(int number){
   stringstream ss;
   ss << number;
   return ss.str();
}

int main(){
 string f_n;
 FILE* outp;
 unsigned long int *mas = new unsigned long int [32769];
 int rand_N = 10000000;
 
 for(int n = 2; n < 32769; n *= 2){

  for(int i = 0; i < n; i++)
   mas[i] = 0;

  f_n = convertInt(n);
  f_n = f_n + "_tsr.txt";

  outp = fopen(f_n.c_str(), "w");
  for(int i = 0; i < rand_N; i++){
   mas[rand()%n]++;
   srand(time(0));
  }

  for(int i = 0; i < n; i++)
   fprintf(outp, "%d\t%f\n", i, 1.0*mas[i]/rand_N);
  fclose(outp);
 }
 return 0;
}

Для тестирования использовалась Visual Studio 2008 с её компилятором. Не думаю, что на gcc отличия будут большие(в среднем конечно) и получится примерно та же картинка.

Немного об абстракции

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

 Разнообразие объектов в зависимости от уровня абстракции   

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

При отсутствии абстракции сложность уходит в бесконечность.

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