Como se discute en los comentarios, he aquí una idea para un algoritmo voraz. Básicamente, el objetivo es encontrar de forma iterativa el rectángulo que se minimice $|w-w'|$ en esa situación en particular. Una combinación de un máximo de subarray y mínimo subarray algoritmo encuentra un rectángulo que viene bastante cerca, si se utiliza el promedio de la resultante subarray como $a_i$. El único problema es que los valores por debajo de la media (o por encima, por la negativa rectángulos) ya debería contar como negativo (positivo); no sé si hay una manera simple de resolver esto (y es un poco difícil de explicar, también). Una forma de reducir este problema es recortar el rectángulo si que mejora el resultado. Otra es limitar el rectángulo de longitud, como el tiempo de rectángulos tienden a tener un mayor porcentaje de los valores por debajo de la media. Puede haber mejores opciones, posiblemente dependiendo de la forma de la señal.
De todos modos, necesitaba experimentar con esto de averiguar si realmente funciona, así que aquí está el código (en C, pero debe ser fácil de puerto):
/* Type of index into signal. */
typedef unsigned int Index;
/* This should be a floating-point type. If you decide to use integers, change
the division to round away from zero instead of towards zero. */
typedef double Value;
/* Represents a rectangle at [start,end). */
typedef struct {
Index start, end;
Value height;
} Rectangle;
/* Trims the rectangle by removing samples whose absolute value will increase
when preliminaryHeight is subtracted. */
static Index optimizeRectangle(const Value signal[], Rectangle *rectangle, Value preliminaryHeight)
{
Index result = 0;
/* Divided into two cases for better efficiency. */
if (preliminaryHeight >= 0) {
while (rectangle->start + 1 < rectangle->end) {
Value value = signal[rectangle->start];
if (preliminaryHeight - value < value) {
break;
}
rectangle->start++;
rectangle->height -= value;
result++;
}
while (rectangle->start + 1 < rectangle->end) {
Value value = signal[rectangle->end - 1];
if (preliminaryHeight - value < value) {
break;
}
rectangle->end--;
rectangle->height -= value;
result++;
}
} else {
while (rectangle->start + 1 < rectangle->end) {
Value value = signal[rectangle->start];
if (preliminaryHeight - value > value) {
break;
}
rectangle->start++;
rectangle->height -= value;
result++;
}
while (rectangle->start + 1 < rectangle->end) {
Value value = signal[rectangle->end - 1];
if (preliminaryHeight - value > value) {
break;
}
rectangle->end--;
rectangle->height -= value;
result++;
}
}
return result;
}
/* Finds the rectangle that seems most appropriate at the moment, and whose
length is at most maxResultLength.
Returns the original length of the rectangle, before optimization. This
value should be used for maxResultLength in the next iteration. If it is
zero, the entire signal was zero. */
static Index findRectangle(const Value signal[], Index signalLength, Rectangle *result, Index maxResultLength)
{
result->start = result->end = 0;
result->height = 0;
Value resultHeightAbs = 0;
/* Start index and height of current maximum and minimum subarrays. */
Index posStart = 0, negStart = 0;
Value posHeight = 0, negHeight = 0;
for (Index index = 0; index < signalLength; index++) {
/* If the length of a subarray exceeds maxResultLength, backtrack
to find the next best start index. */
Index nextIndex = index + 1;
if (nextIndex - posStart > maxResultLength && index > 0) {
Index oldStart = posStart;
posStart = index;
posHeight = 0;
Value newHeight = 0;
for (Index newStart = index - 1; newStart > oldStart; newStart--) {
newHeight += signal[newStart];
if (newHeight > posHeight) {
posStart = newStart;
posHeight = newHeight;
}
}
}
if (nextIndex - negStart > maxResultLength && index > 0) {
Index oldStart = negStart;
negStart = index;
negHeight = 0;
Value newHeight = 0;
for (Index newStart = index - 1; newStart > oldStart; newStart--) {
newHeight += signal[newStart];
if (newHeight < negHeight) {
negStart = newStart;
negHeight = newHeight;
}
}
}
/* Extend the subarrays. */
Value value = signal[index];
posHeight += value;
negHeight += value;
/* Throw away the entire maximum subarray if it is negative or zero. */
if (posHeight <= 0) {
posStart = nextIndex;
posHeight = 0;
} else {
/* Update result if current maximum subarray is better. */
if (posHeight > resultHeightAbs) {
result->start = posStart;
result->end = nextIndex;
result->height = posHeight;
resultHeightAbs = posHeight;
}
}
/* Throw away the entire minimum subarray if it is positive or zero. */
if (negHeight >= 0) {
negStart = nextIndex;
negHeight = 0;
} else {
/* Update result if current minimum subarray is better. */
Value negHeightAbs = -negHeight;
if (negHeightAbs > resultHeightAbs) {
result->start = negStart;
result->end = nextIndex;
result->height = negHeight;
resultHeightAbs = negHeightAbs;
}
}
}
Index resultLength = result->end - result->start;
if (!resultLength) {
/* Return now to avoid division by zero. */
return 0;
}
/* Trim rectangle. */
Value normalizedHeight;
Index trimLength;
do {
normalizedHeight = result->height / (result->end - result->start);
trimLength = optimizeRectangle(signal, result, normalizedHeight);
} while (trimLength);
/* Normalize height. */
result->height = normalizedHeight;
return resultLength;
}
/* Subtracts a rectangle from a signal. */
static void subtractRectangle(Value signal[], const Rectangle *rectangle)
{
for (Index index = rectangle->start; index < rectangle->end; index++) {
signal[index] -= rectangle->height;
}
}
/* Decomposes a signal into rectangles. Stores at most resultLength rectangles
in result, and returns the actual number of rectangles written. All
rectangles are subtracted from the signal. */
unsigned int extractRectangles(Value signal[], Index signalLength, Rectangle result[], unsigned int resultLength)
{
Index rectangleLength = signalLength;
for (unsigned int resultIndex = 0; resultIndex < resultLength; resultIndex++) {
Rectangle *rectangle = &(result[resultIndex]);
rectangleLength = findRectangle(signal, signalLength, rectangle, rectangleLength);
if (!rectangleLength) {
return resultIndex;
}
subtractRectangle(signal, rectangle);
}
return resultLength;
}
La función principal es la extractRectangles. Devuelve el rectángulos como una matriz, pero esto se puede adaptar fácilmente. La diversión.
