buscar_maximo.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""

Ref: http://listas.python.org.ar/pipermail/pyar/2014-January/027667.html

Sea f una función continua en [a, b], tal que f presenta un único
máximo en [a, b]. Por ejemplo f puede ser una parábola con los cuernos
hacia abajo, y f no puede ser de la forma de los cuernos de Sauron
(varios picos mas pequeños).

Escribir un programa que encuentre el x entre [a, b] tal que f(x) sea
el máximo de la función.

  def buscarmax(f, a, b):
    # ...

"""

import math
import timeit

EPS=0.001 # 1e-9

def busqueda_ternaria(f, a, b):
    assert a < b
    a = float(a)
    b = float(b)

    while (b - a) > EPS:
        p1 = a + (b - a) / 3.0
        p2 = a + 2 * (b - a) / 3.0

        if f(p1) < f(p2):
            a = p1
        else:
            b = p2
    return a + (b - a) / 2

# Idea de @ajlopez
def busqueda_segmentos(f, a, b):
    assert a < b

    a = float(a)
    b = float(b)
    N = 32 # 32 es el que mejor resulta con el f de ejemplo, YMMV

    while b - a > EPS:
        #L = [((a + x * (b-a) / (N))) for x in xrange(N+1)]
        #d1 = f(L[0]) < f(L[1])
        d1 = f(a) < f(a + (b-a) / N)
        for x in xrange(N):
            #d2 = f(L[x]) < f(L[x + 1])
            d2 = f(a + (x * (b-a)) / N) < f(a + ((x+1) * (b-a)) / N)
            if d1 != d2:
                a = a + ((x-2) * (b-a)) / N # L[x-2]
                b = a + ((x+1) * (b-a)) / N # L[x+1]
                break
    return a + (b - a) / 2

# Idea de @fisadev
def busqueda_hill_climbing(f, a, b):
    actual = (a + b) / 2  # medio entre a y b
    E = EPS # Para que sea justo ;)
    while True:
        # mejor_vecino = vecino con f más alto en vecinos(actual)
        if f(actual + E) > f(actual - E):
            mejor_vecino = actual + E
        else:
            mejor_vecino = actual - E

        if f(actual) > f(mejor_vecino):
            return actual
        else:
            actual = mejor_vecino

def d(f, x, e):
    return (f(x + e) - f(x)) / e

def f_prima(f, x):
    d1 = d(f, x, EPS)
    d2 = d(f, x, -EPS)

    return d1 + (d2 - d1) / 2.0

def busqueda_derivada(f, a, b):
    assert a < b
    a_orig = a = float(a) + EPS
    b_orig = b = float(b) - EPS

    while (b - a) > EPS:
        m = a + (b - a) / 2.0
        s1 = int(math.copysign(1, f_prima(f, a)))
        s2 = int(math.copysign(1, f_prima(f, m)))
        s3 = int(math.copysign(1, f_prima(f, b)))

        if s1 == s2:
            a = m
        else:
            b = m

    return max(f(a_orig), f(b_orig), a + (b - a) / 2)

def f1(x):
    return -((x-1)*(x-1))

START = -50
END = 50

def test1():
    busqueda_derivada(f1, START, END)

def test2():
    busqueda_ternaria(f1, START, END)

def test3():
    busqueda_segmentos(f1, START, END)

def test4():
    busqueda_hill_climbing(f1, START, END)

if __name__ == '__main__':
    print "por derivada: ", busqueda_derivada(f1, START, END)
    print "por ternaria: ", busqueda_ternaria(f1, START, END)
    print "por segmentos: ", busqueda_segmentos(f1, START, END)
    print "por hill: ", busqueda_hill_climbing(f1, START, END)

    T = 10000

    print "timeit derivada: ", timeit.timeit("test1()", setup="from __main__ import test1", number=T)
    print "timeit ternaria: ", timeit.timeit("test2()", setup="from __main__ import test2", number=T)
    print "timeit segmentos: ", timeit.timeit("test3()", setup="from __main__ import test3", number=T)
    print "timeit hill: ", timeit.timeit("test4()", setup="from __main__ import test4", number=T)