Однородный итеративный семплинг в многомерном пространстве?

Question

[DaylightIsBurning]

Есть функция от Nd параметров f(x[1], x[2], x[3],…, x[Nd]). Даны граничные значения для каждого из параметров. Я хочу итеративно вычислять значения этой функции так, что бы чем дальше, тем больше была плотность точек перебора, при чем плотность росла равномерно для всех аргументов, во всём пространстве.

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

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

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

Comments

[DaylightIsBurning]

Итак, в одномерном случае получилось всё красиво:

//return integer with only the highest previosly set bit 1, rest - 0, which results in 2^((int)log2(n))
uint32_t hibit(uint32_t n) {  //13 bitwise instructions -> 13 time units on GPUs
	n |= (n >> 1);
	n |= (n >> 2);
	n |= (n >> 4);
	n |= (n >> 8);
	n |= (n >> 16);
	return n - (n >> 1); // return (n & ~(n >> 1));
}
/* Convert iterator value to the value of the argument in the range (0.0..1.0)
 * Consecutive iterations will produce more and more fine-grained sampling of argument.
*/
float argValue(uint32_t i) {
	uint32_t p, q;
	uint32_t h = hibit(i);
	p = ((i & (h - 1)) << 1) + 1;
	q = h << 1;
	return (float) p / q;
}
int main(int argc, char *argv[]) {
	unsigned N1 = 0xFF;
	for (unsigned i = 1; i < N1; i++) {
		std::cout << argValue(i) << ", ";
	}
	return false;
}
//0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 0.0625, 0.1875, 0.3125, 0.4375, 0.5625, 0.6875, 0.8125, 0.9375, 0.03125, 0.09375...

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

[DaylightIsBurning]

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

//return integer with only the highest previosly set bit 1, rest - 0, which results in 2^((int)log2(n))
uint32_t hibit(uint32_t n) {  //13 bitwise instructions -> 13 time units on GPUs
	n |= (n >> 1);
	n |= (n >> 2);
	n |= (n >> 4);
	n |= (n >> 8);
	n |= (n >> 16);
	return n - (n >> 1); // return (n & ~(n >> 1));
}
/* Convert iterator value to the value of the argument in the range (0.0..1.0)
 * Consecutive iterations will produce more and more fine-grained sampling of argument.
 */
float argValue(uint32_t i) {
	uint32_t p, q;
	uint32_t h = hibit(i);
	p = ((i & (h - 1)) << 1) + 1;
	q = h << 1;
	return (float) p / q;
}

void pq(unsigned i, uint32_t &p, uint32_t &q) {
	uint32_t h = hibit(i);
	p = ((i & (h - 1)) << 1) + 1;
	q = h << 1;
}

int main(int argc, char *argv[]) {
	unsigned N1 = 0xFF + 1;
	int n = 1;
	std::cout << "n\tx\ty\tpx\tqx\tpy\tqy\n";
	uint32_t px, qx, py, qy;
	for (unsigned i = 1; i < N1; i++) {
		for (unsigned j = 1; j < i; j++) {

			pq(j, px, qx);
			pq(i, py, qy);
			std::cout << n++ << '\t' << argValue(j) << '\t' << argValue(i)
							<< '\t' << px << '\t' << qx << '\t' << py << '\t' << qy
							<< '\n';
		}
		for (unsigned j = 1; j <= i; j++) {
			pq(i, px, qx);
			pq(j, py, qy);
			std::cout << n++ << '\t' << argValue(i) << '\t' << argValue(j)
							<< '\t' << px << '\t' << qx << '\t' << py << '\t' << qy
							<< '\n';
		}

	}
	return false;
}

Над приведением его к единственному итератору и обобщению на N-мерный случай еще думаю.

Answer 1

[DaylightIsBurning]

Вроде придумал, для одномерного случая:
{0/1, 1/1}; {1/2}; {1/4, 3/4}; {1/8, 3/8, 5/8, 7/8}; {1/16,3/16,5/16,7/16,9/16,11/16,13/16,15/16};…
То есть в числителе последовательно перебираем все нечетные, в знаменателе чем глубже, тем больше степень двойки. Не надо никаких if и хранения «пройденных» значений.

Answer 2

[Sayonji]

Да точно так же можно, как и при делении на два.
Предположим, аргументы нормированны, то есть от 0 до 1.
Тогда по шагам:
1) точка 0.5; 0.5; 0.5,…
2) точки (0.25a), (0.25b), (0.25c),… | a,b,c = 1...3
3) точки (0.125a), (0.125b), (0.125c),… | a,b,c = 1...7
…
Ну, само собой, каждый раз берем те, которые еще не брали.

Как получать точки для итерации?

iters = new int[Nd] ;
for (i = 1; ; i++) { // номер итерации
    step = 1 / 2^i; // полдлины интервала семплирования
    for ( /* iters пробегает все возможные наборы из чисел от 1 до (2^i - 1) */ ) {
        if ( /* все числа в iters четные */ )
            continue; // уже брали раньше
        sample F(iters[0] * step, iters[1] * step, iters[2] * step, ...)
    }
}

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

Как добиться однородной плотности по всем аргументам. Если у них разные, но не слишком, диапазоны, нормировать их нужно самым большим из них. Т.е. сдвинуть все в [0; x], где x будет равно 1 только для самого длинного диапазона. А потом во втором форе (который будет работать как куча вложенных друг в друга) отсекать лишние хвосты. Т. е. если, например, 0.125 * iters[2] уже больше, чем величина х для третьего аргумента, то по этой переменной дальше не итерируем. Перед входом во второй for эти границы можно просчитать, это совсем несложно.
И последнее. Если диапазоны отличаются значительно, то может иметь смысл как-то посчитать, с какой итерации следует начинать цикл, т.к. часть из них будет пропущена. Тем больше, чем больше различаются диапазоны.
Как-то так, я думаю.

Comments

[Sayonji]

Про последний случай. Мне кажется, что пропускать можно, пока самый короткий делить на самый длинный меньше, чем step. Ну или что-то в этом духе.
Вообще, описанный подход не даст супер однородной плотности. Скажем, если размер области опеределения ~ 10 x 5, то первая точка вообще окажется где-то на одной из границ. Но слишком неприятные уплотнения получаться не должны. Косяки будут возле границ и будут становится незначительными по мере увеличения номера итерации.
А если даже такое неприемлемо, то тут надо думать насчет общего кратного длин диапазонов, и думать, и думать…

[DaylightIsBurning]

небольшие неоднородности на границах — не критично. Я как раз и думаю о том, как реализовать эффективный алгоритм пропуска ненужных итераций. В идеале нужна функция 1, 2, 3, 4, 5, 6, 7, 8… -> 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 0.0625… Всё это в Nd измерениях.

[Sayonji]

Ненужные итерации пропускаются лишь в самом начале. (У меня итерация — это когда уменьшается интервал семплирования)
Я предложил, кажется, вариант в первом комментарии, как сразу вычислить, с какой итерации начинать.
А если вы про пропуск точек, которые уже были использованы, то чем вам не нравится делать continue; при всех числах в iters четных? Вроде бы не очень сильные потери выйдут. Только ~ 1/2^Nd доля наборов является такой.

[DaylightIsBurning]

if — это медленнее, а на некоторых платформах (openCL/GPU) — убийство вообще :(.
В одномерном случае я пока пришел вот к такому:

uint32_t hibit(uint32_t n) {
    n |= (n >>  1);
    n |= (n >>  2);
    n |= (n >>  4);
    n |= (n >>  8);
    n |= (n >> 16);
    return(n & ~(n >> 1));//return n - (n >> 1);
}
int main(int argc, char *argv[])
{
	const int Nd=1;
	const int n[Nd]={66};
	const int d=0;
	for(uint32_t i=1; i<n[d]; i++)
	{
		uint32_t h=hibit(i);
		uint32_t p= (i&(h-1))*2 + 1;
		std::cout<<"h:"<<h<<'\t'<<"p:"<<p<<'\n';
	}

}

[DaylightIsBurning]

точнее,

	for (uint32_t i = 1; i < n[d]; i++) {
		uint32_t h = hibit(i);
		uint32_t p = (i & (h - 1)) * 2 + 1;
		uint32_t q = 2 * h;
		std::cout << "h:" << h << "\tp:" << p << "\tq:" << q << "\tx:"
				<< (float) p / q << '\n';
	}

Answer 3

[DaylightIsBurning]

-