СОбственно вопрос в следующем

нужно вычислять кусочно линейную функцию и обратную к ней

Дан отрезок [A,B] и набор вложенных не пересекающихся отрезков
[ai,bi] i = 1..n

кусочно — линейная функция выглядит так:

                       /
                      /
                     /
              ______/
             /
            /
    _______/
   /
  /
 /
 A  a1     b1  a2  b2  B

формула

f: [A,B] -> [A`,B`]
f(x) = ai — delta(x) <=> x in [ai,bi]
f(x) = x — delta(x) <=> x not in [ai,bi] , foreach i 1..n

, где
delta(x) = d1 + ... + dk , где k = max {i | x>ai && x>bi}
di = bi-ai

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

по сути задача сводится к следующей
дан отрезок [A, B] и набор точек на нем A <= r1 < r2 < ... < rn <= B
дан x in [A, B] как найти k такое что
rk < x < r(k+1)

пока в голове только бинарное дерево из r1 < r2 < ... < rn

Здравствуйте, AlexDoberman, Вы писали:

AD>СОбственно вопрос в следующем

AD>нужно вычислять кусочно линейную функцию и обратную к ней

AD>Дан отрезок [A,B] и набор вложенных не пересекающихся отрезков
AD>[ai,bi] i = 1..n

AD>кусочно — линейная функция выглядит так:


AD>

AD>                       /
AD>                      /
AD>                     /
AD>              ______/
AD>             /
AD>            /
AD>    _______/
AD>   /
AD>  /
AD> /
AD> A  a1     b1  a2  b2  B
AD>

AD>формула

AD>f: [A,B] -> [A`,B`]
AD> f(x) = ai — delta(x) <=> x in [ai,bi]
AD> f(x) = x — delta(x) <=> x not in [ai,bi] , foreach i 1..n

AD>, где
AD> delta(x) = d1 + ... + dk , где k = max {i | x>ai && x>bi}
AD> di = bi-ai

AD>необходимо быстро вычислять значение этой функции и обратной для нее.

AD>по сути задача сводится к следующей
AD>дан отрезок [A, B] и набор точек на нем A <= r1 < r2 < ... < rn <= B
AD>дан x in [A, B] как найти k такое что
AD>rk < x < r(k+1)

AD>пока в голове только бинарное дерево из r1 < r2 < ... < rn
std::sort, std::lower_bound, std::upper_bound

Здравствуйте, AlexDoberman, Вы писали:

>...
AD>по сути задача сводится к следующей
AD>дан отрезок [A, B] и набор точек на нем A <= r1 < r2 < ... < rn <= B
AD>дан x in [A, B] как найти k такое что
AD>rk < x < r(k+1)

AD>пока в голове только бинарное дерево из r1 < r2 < ... < rn

Вместо дерева лучше воспользоваться бинарным поиском по отсортированному массиву.
см. std::lower_bound

Здравствуйте, iriska2, Вы писали:

I>Здравствуйте, AlexDoberman, Вы писали:

AD>>СОбственно вопрос в следующем

AD>>нужно вычислять кусочно линейную функцию и обратную к ней

AD>>Дан отрезок [A,B] и набор вложенных не пересекающихся отрезков
AD>>[ai,bi] i = 1..n

AD>>кусочно — линейная функция выглядит так:


AD>>

AD>>                       /
AD>>                      /
AD>>                     /
AD>>              ______/
AD>>             /
AD>>            /
AD>>    _______/
AD>>   /
AD>>  /
AD>> /
AD>> A  a1     b1  a2  b2  B
AD>>

AD>>формула

AD>>f: [A,B] -> [A`,B`]
AD>> f(x) = ai — delta(x) <=> x in [ai,bi]
AD>> f(x) = x — delta(x) <=> x not in [ai,bi] , foreach i 1..n

AD>>, где
AD>> delta(x) = d1 + ... + dk , где k = max {i | x>ai && x>bi}
AD>> di = bi-ai

AD>>необходимо быстро вычислять значение этой функции и обратной для нее.

AD>>по сути задача сводится к следующей
AD>>дан отрезок [A, B] и набор точек на нем A <= r1 < r2 < ... < rn <= B
AD>>дан x in [A, B] как найти k такое что
AD>>rk < x < r(k+1)

AD>>пока в голове только бинарное дерево из r1 < r2 < ... < rn
I>std::sort, std::lower_bound, std::upper_bound

спасибо на ссылку, но время работы все равно логарифмическое

Здравствуйте, AlexDoberman, Вы писали:
AD>спасибо на ссылку, но время работы все равно логарифмическое

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

Например, если их распределение достаточно равномерно, можно разбить исходный отрезок на N равных отрезков, так чтобы в каждый отрезок попало не более M точек перегиба.
Сперва найти отрезок делением, (+ округление до целого) , а потом снова бинарным поиском O(log(M)).

Здравствуйте, Chorkov, Вы писали:

C>Здравствуйте, AlexDoberman, Вы писали:
AD>>спасибо на ссылку, но время работы все равно логарифмическое

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

C>Например, если их распределение достаточно равномерно, можно разбить исходный отрезок на N равных отрезков, так чтобы в каждый отрезок попало не более M точек перегиба.
C>Сперва найти отрезок делением, (+ округление до целого) , а потом снова бинарным поиском O(log(M)).

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

Здравствуйте, AlexDoberman, Вы писали:

AD>>>пока в голове только бинарное дерево из r1 < r2 < ... < rn
I>>std::sort, std::lower_bound, std::upper_bound

AD>спасибо на ссылку, но время работы все равно логарифмическое

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

You will always get what you always got
  If you always do  what you always did

Здравствуйте, AlexDoberman, Вы писали:

AD>нужно вычислять кусочно линейную функцию и обратную к ней

Какой диапазон? Можно ведь и табличным способом реализовать, в конце концов.

Я бы предложил вот такой общий подход.

1. Отказаться от единственного представления кусочно-линейной кусочно-константной функции и обратной функции.
Пусть будут две кусочно-линейные функции одинакового вида.

y(x) = {
  [a0..b0) -> [y0..y1),
  [b0..a1) -> [y1..y1), --- это константа
  [a1..b1) -> [y1..y2),
  [b1..a2) -> [y2..y2), --- это константа
  ...
}

x(y) = {
  [y0..y1) -> [a0..b0),
  [y1..y2) -> [a1..b1), --- между y1-e и y1 был разрыв
  ...
}

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

2. Любым способом решаем задачу поиска диапазона по точке.
Можно двоичный поиск, можно двоичное дерево, можно Б-дерево большой арности.

3. Каждому диапазону соответствуют 4 числа: [x1..x2)->[y1..y2), но для вычисления нужны только два: f(x) = (x-x1)*(y2-y1)/(x2-x1) = (x-x1)*C

Таким образом, наиболее компактное представление будет таким

// табличное представление
std::vector<float> Xs; // ключи: края диапазонов
std::vector<float> Cs; // коэффициенты наклона
// или такое
std::map<float,float> XtoCs;

float x; // аргумент функции

int i = std::distance(Xs.begin(), std::lower_bound(Xs.begin(), Xs.end(), x));
float y = (x-Xs[i])*Cs[i];

std::map<float,float>::const_iterator it = XtoCs.lower_bound(x);
float y = (x-it.first)*it.second;

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

if(x < Xmin) return Ymin;
if(x > Xmax) return Ymax;
..... // далее как показано выше

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

Если количество очень большое, то я бы посмотрел в сторону Б-деревьев.

Здравствуйте, Кодт, Вы писали:

К>Здравствуйте, AlexDoberman, Вы писали:

AD>>нужно вычислять кусочно линейную функцию и обратную к ней

К>Какой диапазон? Можно ведь и табличным способом реализовать, в конце концов.

количество кусков не большой , не более 100 — 1000
много вызовов функции вычисления

Здравствуйте, AlexDoberman, Вы писали:

AD>количество кусков не большой , не более 100 — 1000
AD>много вызовов функции вычисления

Тогда нужно выбрать между наивным двоичным поиском (std::lower_bound) и каким-нибудь рукодельным оптимизированным (двоичный по степеням 2, затем линейный по небольшому массиву). И профилировать.
Какие-то хитровывернутые структуры данных здесь будут во вред, потому что 1000 флотов хорошо ложится в кэш.
Да, кстати, если есть возможность, — массивы Xs и Cs разместить в памяти с выравниванием на линейку кэша.

Если вызовы функции как-то связаны между собой (например, идёт вычисление значений в некоторой окрестности, или монотонная последовательность значений), то стоит подумать о пакетном вычислении — это сэкономит процедуру поиска.

Если абсциссы кусков расположены регулярно (например, кратны (B-A)/1024), то табличное решение будет отчаянно рулить. Вместо 100-1000 интервалов взять ровно 1024 интервала, — и вперёд.

Я что-то с арифметикой раздружился.

y(x) = y1 + (x-x1)*(y2-y1)/(x2-x1) = x*K + C
То есть, каждому диапазону соответствует не одно, а два числа

struct F
{
  std::vector<float> Xs, Ks, Cs;

  float operator()(float x) const
  {
    i = std::lower_bound(Xs.begin(), Xs.end(), x) - Xs.begin();
    return x*Ks[i]+Cs[i];
  }
};

Кстати, такая схема позволит обобщённо работать и с левым-правым краем.
Дополним Xs[0] = -INF, Xs[n+1] = B. Всего у нас будет не n, а n+2 интервала.

Для x<A мы попадаем в интервал №0, (-oo,A), и вычисляем x*Ks[0]+Cs[0] — какую-то линейную (или константную) левую ветвь.
Для x>B — попадаем в интервал №(n+1), [B,+oo)

В обоих случаях бесконечность в вычислениях не участвует. В отличие от схемы (x-x1)*K+C.

Здравствуйте, Кодт, Вы писали:

К>Здравствуйте, AlexDoberman, Вы писали:

AD>>количество кусков не большой , не более 100 — 1000
AD>>много вызовов функции вычисления

К>Тогда нужно выбрать между наивным двоичным поиском (std::lower_bound) и каким-нибудь рукодельным оптимизированным (двоичный по степеням 2, затем линейный по небольшому массиву). И профилировать.
К>Какие-то хитровывернутые структуры данных здесь будут во вред, потому что 1000 флотов хорошо ложится в кэш.
К>Да, кстати, если есть возможность, — массивы Xs и Cs разместить в памяти с выравниванием на линейку кэша.
спасибо за совет

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

К>Если абсциссы кусков расположены регулярно (например, кратны (B-A)/1024), то табличное решение будет отчаянно рулить. Вместо 100-1000 интервалов взять ровно 1024 интервала, — и вперёд.
не совсем понял. куски у меня произвольные , а если бы они были периодические тогда я бы замутил прямую формулу.

Здравствуйте, AlexDoberman, Вы писали:

К>>Да, кстати, если есть возможность, — массивы Xs и Cs разместить в памяти с выравниванием на линейку кэша.
AD>спасибо за совет

Только обрати внимание на мою поправку про арифметику.


К>>Если вызовы функции как-то связаны между собой (например, идёт вычисление значений в некоторой окрестности, или монотонная последовательность значений), то стоит подумать о пакетном вычислении — это сэкономит процедуру поиска.
AD>я думаю можно кешировать найденный интервал

Это зависит от того, какая автокорреляция у входных данных. Кстати, что за вычислительная задача?


К>>Если абсциссы кусков расположены регулярно (например, кратны (B-A)/1024), то табличное решение будет отчаянно рулить. Вместо 100-1000 интервалов взять ровно 1024 интервала, — и вперёд.
AD>не совсем понял. куски у меня произвольные , а если бы они были периодические тогда я бы замутил прямую формулу.

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


Кстати, ещё один способ ускорения поиска.

vector<float> Xs;    // абсциссы - как обычно
vector<float> Ks,Cs; // линейные функции - как обычно

float X0, DX;
vector<int> Ls, Rs; // индексы абсцисс для регулярной сетки:

int Xk(float x) { return floor((x-X0)/DX); } // обратный закон регулярной сетки
// эти функции нужны только для пояснения
float Xk(int k) { return X0 + k*DX; } // регулярная сетка
float Xl(int k) { return Xs[Ls[k]]; } // абсцисса, ближайшая слева к отрезку регулярной сетки
float Xr(int k) { return Xs[Rs[k]]; } // абсцисса, ближайшая справа к этому отрезку
// Xl(k) <= Xk(k) < Xk(k+1) < Xr(k)

float k;

int k = Kx(x);
//          Xk(k) <= x < Xk(k+1)
// Xl(k) <= Xk(k)        Xk(k+1) < Xr(k)
int l = Ls[k], r = Rs[k]

int i = lower_bound(Xs.begin()+l, Xs.begin()+r, x) - Xs.begin(); // двоичный поиск только в уточнённом множестве абсцисс

float y = x*Ks[i]+Cs[i];

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

int K0; // индекс "нулевого" узла сетки
float X0; // "нулевой" узел - мода
float D; // шаг сетки около "нуля"
float E = 2; // показатель степенного закона

// степенная (линейная, квадратичная...) сетка
float Xk(int k)
{
  return k < K0 ? X0-pow(K0-k,E)*D
                : X0+pow(k-K0,E)*D;
}
int Kx(float x)
{
  return x < X0 ? K0-(int)ceil (pow((X0-x)/D,1/E))
                : K0+(int)floor(pow((x-X0)/D,1/E));
}

// сетка с квадратичным законом выглядит так
//                                K0
//  0         1       2     3   4 5 6   7     8       9        10
// -|---------|-------|-----|---|-|-|---|-----|-------|---------|-
//                                X0

Здравствуйте, AlexDoberman, Вы писали:

AD>по сути задача сводится к следующей
AD>дан отрезок [A, B] и набор точек на нем A <= r1 < r2 < ... < rn <= B
AD>дан x in [A, B] как найти k такое что
AD>rk < x < r(k+1)

AD>пока в голове только бинарное дерево из r1 < r2 < ... < rn


1. Бинарный поиск.
2. Если длины отрезков соизмеримы — то хэширование. То есть наложить на отрезок массив, каждый элемент m[i] которого хранит список отрезков, которые задевают промежуток от i до i+1.