TI/Programowanie dla Fizyków Medycznych/Optymalizacja

import numpy as np
import pylab as py
import time
import scipy.optimize as so

class stoper(object):
    def __init__(self,func):
        self.func=func
    def __call__(self,*args,**kwargs):
        print 'start funkcji '+str(self.func.__name__)
        t0=time.time()
        result=self.func(*args,**kwargs)
        tk=time.time()
        print 'stop funkcji '+str(self.func.__name__)
        print 'czas dzialania: '+str(tk-t0)+' sekund'
        return result

licznikTestowej=0

def testowa(x):
    global licznikTestowej
    licznikTestowej+=1
    return 1/x+np.exp(x)

#xtest=np.arange(0.2,2,0.01)
#ytest=[testowa(x) for x in xtest]
#py.plot(xtest,ytest)
#py.show()

@stoper
def bruteForce(func,xmin,xmax,args=(),xtol=0.01):
    xlist=np.arange(xmin,xmax,xtol)
    ylist=[func(x,*args) for x in xlist]
    #py.plot(xlist,ylist)
    #py.show()
    return xlist[ylist.index(max(ylist))]

#tutaj lepiej bez stopera!
def twoMidPointsR(func,xmin,xmax,args=(),xtol=0.01):
    if xmax-xmin<xtol: return 0.5*(xmax+xmin)
    xL=xmin+(xmax-xmin)/3.0
    xR=xmax-(xmax-xmin)/3.0
    fxL=func(xL,*args)
    fxR=func(xR,*args)
    if fxL>fxR:
        return twoMidPointsR(func,xmin,xR,args,xtol)
    else:
        return twoMidPointsR(func,xL,xmax,args,xtol)

#szuka maximum funkcji
@stoper
def GoldenRatioRearch(func,xmin,xmax,args=(),xtol=0.01):
    golden=0.5*(np.sqrt(5.0)-1.0)
    xL=xmax-golden*(xmax-xmin)
    xR=xmin+golden*(xmax-xmin)
    fxL=func(xL,*args)
    fxR=func(xR,*args)
    while xmax-xmin>xtol: 
        if fxL<fxR:
            xmin=xmin
            xmax=xR
            xR=xL
            fxR=fxL
            xL=xmax-golden*(xmax-xmin)
            fxL=func(xL,*args)
        else:
            xmin=xL
            xmax=xmax
            xL=xR
            fxL=fxR
            xR=xmin+golden*(xmax-xmin)
            fxR=func(xR,*args)
    return 0.5*(xmax+xmin)


#test metody zlotego podzialu
#print GoldenRatioRearch(testowa,0,5,xtol=0.000001)
#print so.fmin(testowa,np.array([1]))[0]

#przypadek z jednym parametrem
def squares(a,func,xlist,ylist):
    return sum([(func(xlist[i],*a)-ylist[i])**2 for i in range(len(xlist))])

#funkcja liniowa
def liniowa(x,a): return x*a

#funkcja kwadratowa
def kwadratowa(x,a,b=0): return a*x*x+b*x


#generujemy przykladowe xlist
xlist=np.arange(0,1,0.001)
#generujemy wartosci funkcji z szumem
ylist=[kwadratowa(x,1.23,-0.73)+0.000001*np.random.randn() for x in xlist]
#py.plot(xlist,ylist)
#py.show()

#lista przeszukiwanych wartosci parametru a
#alist=np.arange(0,10,0.01)
#wartosci squares dla tych parametrow
#squareslist=[squares(a,liniowa,xlist,ylist)for a in alist]
#zobaczmy wykres
#py.plot(alist,squareslist)
#py.show()

#najlepsze dopasowanie metoda brute force
#print bruteForce(squares,0,10,args=(liniowa,xlist,ylist),xtol=0.01)

#najlepsze dopasowanie metoda golden ration funkcji z 1 parametrem
#print GoldenRatioRearch(squares,0,10,args=(liniowa,xlist,ylist),xtol=0.01)
print so.fmin(squares,(1),args=(liniowa,xlist,ylist))

#najlepsze dopasowanie metoda symplexowa z dwoma parametrami
print so.fmin(squares,(1,0),args=(kwadratowa,xlist,ylist))