import numpy as np
import matplotlib.pyplot as plt
import math
import scipy.special
import time
import os
import random
import copy



def placeCities(N,xmax,ymax):
	c = np.zeros([N,2])
	c[:,0],c[:,1] = map(lambda x: random.randint(0,xmax),c[:,0]),map(lambda x: random.randint(0,xmax),c[:,1])
	return c

#plt.figure(1)
#c = placeCities(5,100,100)
#plt.scatter(c[:,0],c[:,1])
#plt.show()

def distances(c):
	A = np.zeros([len(c),len(c)])
	for i in range(0,len(c)):
		for j in range(i+1,len(c)):
			A[i,j] = np.linalg.norm(c[i,:]-c[j,:])
			A[j,i] = A[i,j]
	return A

def genRoute(N):
	r = np.array(range(0,N),dtype=int)
	random.shuffle(r)
	return r

def rlength(r,R):
	return sum(R[r[1:],r[:-1]])

def genNRoutes(NR,NS):
	R = np.zeros([NR,NS],dtype=int)
	R[:,:] = map (lambda x: genRoute(NS),R)
	return R

def shortestRoute(rn,c):
	m,mi = 0,0
	R = distances(c)
	for i in range(0,rn.shape[0]):
		t = rlength(rn[i],R)
		if t < m:
			m,mi = t,i
	return rn[i]

def duell(rn,c):
	R = distances(c)
	d = range(0,rn.shape[0])
	random.shuffle(d)
	s = np.zeros([rn.shape[0]/2,rn.shape[1]],dtype=int)
	j=0
	for i in range(0,len(d),2):
		s[j] = rn[d[i],:] if rlength(rn[d[i],:],R)<rlength(rn[d[i+1],:],R) else rn[d[i+1],:]
		j+=1
	return s

def plotroute(r,c):
	plt.figure(1)
	plt.scatter(c[:,0],c[:,1])
	plt.plot(c[r][:,0],c[r][:,1])
	plt.show()

def nextGen(rn,M):
	parents = copy.deepcopy(rn)
	children = copy.deepcopy(rn)
	for i in range(0,M):
		c,p1,p2 = random.randint(0,len(rn)-1),random.randint(0,rn.shape[1]-1),random.randint(0,rn.shape[1]-1)
		if p1!=p2:
			rn[c,p2],rn[c,p1] = rn[c,p1],rn[c,p2]

	return parents,children

def tsp():
	cpop = 200
	ccities = 15
	c = placeCities(ccities,100,100)
	R = distances(c) #TODO!!! duell
	rn = duell(duell(duell(duell(genNRoutes(16*cpop,ccities),c),c),c),c)
	#plotroute(shortestRoute(rn,c),c)
	ltest = np.zeros(200);
	for k in range(0,200):
		parents,children = nextGen(rn,1)
		g = np.concatenate((parents,children),axis=0)
		rn = duell(g,c)
		ltest[k] = rlength(shortestRoute(rn,c),R)
	plt.figure(3)
	plt.plot(range(0,200),ltest)
	plt.show()
	plotroute(shortestRoute(rn,c),c)

tsp()





























# def V(d):
# 	return (np.pi**(d/2))/scipy.special.gamma((d/2)+1)



# plt.figure(1)
# plt.plot(range(0,11),[V(i) for i in range(0,11)])
# plt.show()


# def U(u,x):
#     return (1/2.0)*u*x**2

# def p(u,u0,kT,x):
#     return np.exp(-u(u0, x)/kT)/np.sqrt(2*np.pi*kT/u0)

# def mA(U,u0,kT,N,r):
#     ra = random(7**5,0,(2**31)-1)
#     x = np.zeros(N+1)
#     counter = 0.0
#     x[0] = float(ra.next())/(2**31-1)
#     for i in range(1,N):
#         a = -r+(2*r)*float(ra.next())/(2**31-1)
#         U1 = U(u0, x[i-1]+a) - U(u0, x[i-1])
#         if U(u0, x[i-1]+a) <= U(u0, x[i-1]):
#             x[i] = x[i-1] + a
#         else:
#             if float(ra.next())/(2**31-1) <= p(U,U1,kT,x[N-1]):
#                 x[i] = x[i-1]+a
#             else:
#                 x[i] = x[i-1]
#                 counter += 1
#     return x, 1-(counter/N)

# print 'Akzeptanzrate: '+str(mA(U,1.0,2.0,10000,1.8)[1])

# x = np.linspace(-100,100,10001)

# plt.figure()
# plt.hist(mA(U,1.0,2.0,10000,1.9)[0],100,normed=1)
# plt.plot(x,p(U,1,2.,x))
# plt.show()


def random(a,c,m):
	n = int(((np.abs(64979*(round(time.time()*1000))*(os.getpid()-83))%104729)/(104729))*m)
	while 1:
		yield n
		n = (a*n+c)%m

r = random(7**5,0,(2**31)-1)
z = [float(r.next())/((2**31)-1) for i in range(0,10000)]
# print np.average(z)
# print np.std(z)
# plt.hist(z)
#plt.show()


def crudemc(f,a,b,N):
	return ((b-a)/N)*sum(np.vectorize(f)(a+(b-a)*[random(7**5,0,(2**31)-1)for i in range(0,N)]))


#Aufgabe 11.2.1:
def crudemc(f,a,b,N):
	r = random(7**5,0,(2**31)-1)
	print (b-a)/N
	return ((b-a)/N)*sum(f([a+float((b-a)*float(r.next()/(2**31))) for i in range(0,N)]))


def hommc(f,a,b,N,fmin,fmax):
	r = random(7**5,0,(2**31)-1)
	hits = 0
	for i in range(0,N):
		hits+=1 if(f(a+((b-a)*float(r.next())/(2**31-1)))-fmin>(fmax-fmin)*float(r.next())/(2**31-1)) else 0
		print (b-a)
		print ((hits/N)*(fmax-fmin)+fmin)
		return (b-a)*((hits/N)*(fmax-fmin)+fmin)






# def f(x):
# 	return 1/(1+np.power(np.cos(x),2))
# def g(x):
# 	return np.power(x,-12.0)-np.power(x,-6.0)

# def dft_forward(a):
# 	tmp = np.zeros(len(a),dtype="complex")
# 	for k in range(0,len(a)):
# 		for i in range(0,len(a)):
# 			tmp[k] += np.exp(-2j*np.pi*i*k/len(a))*a[i]
# 	return tmp

# def dft_backward(a):
# 	tmp = np.zeros(len(a),dtype="complex")
# 	for k in range(0,len(a)):
# 		for i in range(0,len(a)):
# 			tmp[k] += np.exp(2j*np.pi*i*k/len(a))*a[i]
# 		tmp[k] = tmp[k]/len(a)
# 	return tmp



# signal = np.genfromtxt('/Users/blank/signal.txt')
# signal_fft = np.fft.fft((signal[:,1]))
# plt.figure(1)
# # plt.plot(signal)
# plt.plot(signal_fft)
# # print signal[:,1]
# plt.show()


# guitar = scipy.io.wavfile.read('/Users/blank/uni/comphy/guitar.wav')
# plt.figure(1)
# plt.plot(np.fft.rfft(guitar[1]))
# d2 = np.fft.rfft(guitar[1])
# for x in range(0,len(d2)):
# 	if np.abs(d2[x])<70000:
# 		d2[x]=0
# plt.figure(2)
# plt.plot(d2)
# plt.show()
# rev = np.array(np.fft.irfft(d2),dtype="int16")
# scipy.io.wavfile.write("testfile.wav",guitar[0],rev)

# def V(x):
# 	y = np.zeros_like(x)
# 	for i in range(0,length(x)):
# 		y[i] = x[i]*np.e**((-x[i]**2)/4) if x < 5 and x > -5 else np.inf
# 	return y

# def y_deriv(y,x,E,V):
# 	return y[1],2*(V(x)-E)*y[0]

# x=np.linspace(-5,5.0001,0.01)
# E=-0.99
# y=np.array([0,1])
# y0=np.array([0,1])

# y = scipy.integrate.odeint(y_deriv,y0,x,args=(E,V))



# def skonf(dx,x0):
# 	temp = np.zeros(1.0/dx)
# 	temp[x0] = 1
# 	return temp

# def solve(dx,dt,x0,D,T_max):
# 	steps = int(T_max/dt)
# 	c = np.zeros([1.0/dx,steps])
# 	c[:,0] = skonf(dx,x0)
# 	for i in range(0,steps-1):
# 		for j in range(0,int(1/dx)-1):
# 			c[j,i+1] = c[j,i]+D*(dt/(dx**2))*(c[j+1,i]-2.0*c[j,i]+c[j-1.0,i])
# 	return c



# D = ((1.38068*10**-23)*298)/(6*np.pi*(10**-3)*(7.3*10**-10))
# plt.figure(1)
# plt.plot(solve(1.0/100,16,20,D,3600))
# plt.show()



















def y_(y):
	return alpha*(y**2)*(1-y)-0.01*y

def phi(y1, y2):
	return y2,-np.sqrt((y1+9.81)*(80/0.2))
# y_k+1 = y_k + h*f(t_k,y_k)  k = 0,1,2,...
def euler(f,steps,h,y01,y02):
	y1 = np.zeros(steps)
	y2 = np.zeros(steps)
	y1[0] = y01
	y2[0] = y02
	for j in range(1,steps):
		y11,y21 = f(y1[j-1],y2[j-1])
		y1[j] = y1[j-1]+h*y11
		y2[j] = y2[j-1]+h*y21
	return y1,y2

#Aufgabe 8.2.3:
def rk4(f,steps,h,y01,y02):
	y1 = np.zeros(steps)
	y1[0] = y01
	y2 = np.zeros(steps)
	y2[0] = y02
	for i in range(1,steps):
		k11,k12 = f(y1[i-1],y2[i-1])
		k21,k22 = f(y1[i-1]+(h/2)*k11,y2[i-1]+(h/2)*k12)
		k31,k32 = f(y1[i-1]+(h/2)*k21,y2[i-1]+(h/2)*k22)
		k41,k42 = f(y1[i-1]+h*k31,y1[i-1]+h*k32)
		y1[i]=y1[i-1]+(h/6)*(k11+2*k21+2*k31+k41)
		y2[i]=y2[i-1]+(h/6)*(k12+2*k22+2*k32+k42)
	return y1,y2

def test():
	steps = 20
	i,k = rk4(phi,steps,0.05,3000,3000)
	plt.plot(range(0,steps),k)
	plt.show()
#test()

def rk4(f,steps,h,y0):
	y = np.zeros(steps)
	y[0] = y0
	for i in range(1,steps):
		k1 = f(y[i-1])
		k2 = f(y[i-1]+(h/2)*k1)
		k3 = f(y[i-1]+(h/2)*k2)
		k4 = f(y[i-1]+h*k3)
		y[i]=y[i-1]+(h/6)*(k1+2*k2+2*k3+k4)
	return y

def aufg812():
	global alpha
	steps = 100000
	alpha = 1
	i = euler(y_,steps,0.05,1)
	plt.plot(range(0,steps),i,label="alpha = 1")
	alpha = 0.05
	i = euler(y_,steps,0.05,1)
	plt.plot(range(0,steps),i,label="alpha = 0.05")
	alpha = 0.04
	i = euler(y_,steps,0.05,1)
	plt.plot(range(0,steps),i,label="alpha = 0.04")
	alpha = 0.03
	i = euler(y_,steps,0.05,1)
	plt.plot(range(0,steps),i,label="alpha = 0.03")
	plt.legend()
	plt.show()


def aufg813():
	global alpha
	steps = 10000
	alpha = 1
	i = euler(y_,steps,0.05,1)
	plt.plot(range(0,steps),i,label="alpha = 1")
	alpha = 0.05
	i = euler(y_,steps,0.05,1)
	plt.plot(range(0,steps),i,label="alpha = 0.05")
	alpha = 0.04
	i = euler(y_,steps,0.05,1)
	plt.plot(range(0,steps),i,label="alpha = 0.04")
	alpha = 0.03
	i = euler(y_,steps,0.05,1)
	plt.plot(range(0,steps),i,label="alpha = 0.03")
	alpha = 0.001
	i = euler(y_,steps,0.05,1)
	plt.plot(range(0,steps),i,label="alpha = 0.001")
	alpha = -0.001
	i = euler(y_,steps,0.05,1)
	plt.plot(range(0,steps),i,label="alpha = -0.001")
	plt.legend()
	plt.show()
#aufg813()