#!/usr/bin/python
# -*- coding: utf-8 -*-
"""Invert John Costella’s Magic Kernel by brute force.

I was curious what kind of “pre-emphasis” you need to give to a
sequence of values in order to undo the effects of the convolution
kernel [1 3 3 1], which is [1 1] ★ [1 1] ★ [1 1], a discrete version
of the quadratic uniform cardinal B-spline.  It wasn’t obvious to me
what the solution was, so I thought I’d apply brute-force search.

While this was running, I figured out how to formulate the problem
(for some finite sequence length) as a linear algebra problem and
solved it with Numpy.  This showed me a bug in this code!

After running all night to try to find an input that produced an
impulse [0, 0, 0, 0, 0, 0, 24, 0, 0, 0, 0], an earlier version of this
code came up with

    [24.0, -143.998, 359.994, -239.996, 143.997, -71.998, 23.999, 0.001, -0.0, 0.0, -0.0]

with output

    [0.0, 0.0, -0.0, 0.0, -0.0, 0.0, 24.0, 0.0, -0.0, 0.0, -0.0] 2.33432139019e-08

but it’s looking for the wrong thing.  The first and last few samples
in the output buffer aren’t correct: the first two samples have only
been through the last one and two stages of the pipeline, and the last
three have only been through the first two, one, and zero stages of
the pipeline.  So the values 24 and -144 are totally bogus.  The 360,
-240, 144, -72, 24, 0 sequence is fine though.

So, I thought, I should have tried trimming off those invalid padding
samples from the output.  So I tried that, function
magic_kernel_valid.  Imagine my surprise when, asked to produce

    [0, 0, 0, 1, 0, 0, 0]

in about eight seconds, it converged on

    [0.185, 0.913, -0.22, -0.613, 1.005, -0.956, 0.466, 0.466, -0.838, 0.651, 0.094, -1.399]

which produces the result

    [0.0, -0.0, 0.0, 1.0, 0.0, -0.0, 0.0] 3.85185988877e-31

or, more precisely, including the invalid samples on the ends:

    [0.2652742816945357,
    -0.36099065250966755,
    0.0,
    -1.6653345369377348e-16,
    2.220446049250313e-16,
    0.9999999999999996,
    2.220446049250313e-16,
    -2.220446049250313e-16,
    1.1102230246251565e-16,
    -0.5589317457557798,
    -1.3044639795044046,
    -1.3986356237403164]

So, that’s pretty close to being a valid solution.  Except that
somehow those nonzero samples out at the end might be important!  So
really I only want to search over the input samples corresponding to
valid output samples.  Or, just to be sure, input samples that don't
affect any invalid output samples.

And, sure enough, once that restriction is added back in, it gets slow
again.

"""
from __future__ import print_function
import random


def l2_diff(a, b):
    assert len(a) == len(b), (a, b)
    return sum((ai-bi)**2 for ai, bi in zip(a, b))

def search(f, x, y, norm, mutate):
    last = None

    while True:
        xp = mutate(x, 1 if last is None else last)
        yp = f(xp)
        err = norm(yp, y)
        if last is None or err < last:
            x = xp
            last = err
            yield x, yp, err

def magic_kernel(x):
    x = x[:]
    for i in range(len(x)-3):
        x[i+2] += x[i+3]
        x[i+1] += x[i+2]
        x[i] += x[i+1]
    return x

def magic_kernel_valid(x):
    return magic_kernel(x)[2:-3]

def magic_kernel_pad(x):
    return magic_kernel_valid([0] * 4 + x + [0] * 3)

def mutall(xs, s):
    return [xi + random.random() - 0.5 for xi in xs]

def mutallsqrt(xs, s):
    return [xi + random.normalvariate(mu=0, sigma=0.1*s**0.5) for xi in xs]

def mutint(xs, s):
    return [xi + random.randrange(-1, 2) for xi in xs]

def mutintmore(xs, s):
    return [xi + random.randrange(-3, 4) for xi in xs]

def main():
    y = [0, 0, 0, 1, 0, 0, 0]
    for x, yp, err in search(magic_kernel_pad, [0 for i in range(len(y)-2)], y,
                             l2_diff, mutallsqrt):
        print([round(xi, 3) for xi in x])
        print([round(yi, 3) for yi in yp], err)
        print('#', magic_kernel([0] * 4 + x + [0] * 3))

if __name__ == '__main__':
    main()