Optimisation and Search in Chapter 11 (Machine Learning written by Stephen Marsland)

editing now.....

I'm gonna introduce some methods for optimization problems.....

• Newton method
• Lenvenberg_marquard method
• Conjugate Gradient method

#-*- coding:utf-8 -*-

import numpy as np

def Jacobian(x):
	return np.array([.4*x[0],2*x[1]])
	#return np.array([x[0], 0.4*x[1], 1.2*x[2]])

def Hessian(x):
	return np.array([[.2,0],[0,1]])
	#return np.array([[1., 0 ,0], [0, 0.4, 0], [0, 0, 1.2]])

#conjugate gradient
#initial set parameter x0
def CG(x0):
	#initial parameter
	x = x0
	
	#stopping parameter
	epsilon = 1e-3
	iteration = 3
	
	J = -Jacobian(x)
	p_old = J
	p_new = J
	
	#set numerator of beta
	beta_n = np.dot(J.T, J)
	
	counter = 0
	while counter < iteration and all( abs(p_new) > p_old * epsilon ** 2):
		
		#compute initial alpha_k
		alpha = - np.dot(Jacobian(x).T, p_new) / np.dot(np.dot(p_new.T, Hessian(x)), p_new)
		print alpha
		
		#compute initial 
		dp = np.dot(p_new.T, p_new)
		
		#counter
		counter2 = 0
		
		#Newton Raphson iteration
		while counter2 < iteration and alpha**2 * dp > epsilon**2:
			#update alpha
			alpha = -np.dot(Jacobian(x).T, p_new) / np.dot(p_new.T, np.dot(Hessian(x), p_new))
			
			#update x
			x = x + alpha * p_new
			
			#update dp
			dp = np.dot(p_new.T, p_new)
			
			counter2 +=1
		
		#evaluate delta_f of x_new
		J = - Jacobian(x)
		
		#compute beta(k+1) using Fletcher-Reeves formula
		beta_d = beta_n
		beta_n = np.dot(J.T, J)
		beta = beta_n / beta_d
		
		#update p_new = J + beta * P
		p_old = p_new
		p_new = J + beta * p_new
		counter +=1
		
		print p_new
	print x

#x0 = np.array([-2,2,-2])
x0 = np.array([-2, 2])
CG(x0)

The example of CVXOPT

Hi there!

Today, I would like to introduce cvxopt to you.
CVXOPT is a library for resolving convex optimization like linear programs,
and it is written by Python.

please be careful for other packages when you want to install this in Mac.
because you also need to download "Command Line Tools" through Xcode.

I think this is helpful for you….
satomacoto: Pythonで線形計画法を解いてみる

Anyway, I'll show a usage of this.

at this time, I'd like to resolve this problem.

Problem 1:

Problem 2

from cvxopt import matrix, solvers

A = matrix([[1., -1.], [1., 2.]])
c = matrix([-1., -1.])
b = matrix([2., 8.])
sol=solvers.lp(c,A.T,b)
print(sol['x'])

A = matrix( [[1., 1., 1.], [1., 2., 3.], [-1., 0., 0.], [0., -1., 0.], [0., 0., -1.]])
c = matrix([-2., -3., -1.])
b = matrix([10., 12., 0., 0., 0.])
sol=solvers.lp(c,A.T,b)
print(sol['x'])


>>> from cvxopt import matrix, solvers
>>> 
>>> A = matrix([[1., -1.], [1., 2.]])
>>> c = matrix([-1., -1.])
>>> b = matrix([2., 8.])
>>> sol=solvers.lp(c,A.T,b)
     pcost       dcost       gap    pres   dres   k/t
 0: -6.0000e+00 -6.0000e+00  1e+00  2e-01  2e-16  1e+00
 1: -6.0000e+00 -6.0000e+00  1e-02  2e-03  1e-16  1e-02
 2: -6.0000e+00 -6.0000e+00  1e-04  2e-05  1e-16  1e-04
 3: -6.0000e+00 -6.0000e+00  1e-06  2e-07  2e-16  1e-06
 4: -6.0000e+00 -6.0000e+00  1e-08  2e-09  2e-16  1e-08
Optimal solution found.
>>> print(sol['x'])
[ 4.00e+00]
[ 2.00e+00]

>>> 
>>> A = matrix( [[1., 1., 1.], [1., 2., 3.], [-1., 0., 0.], [0., -1., 0.], [0., 0., -1.]])
>>> c = matrix([-2., -3., -1.])
>>> b = matrix([10., 12., 0., 0., 0.])
>>> sol=solvers.lp(c,A.T,b)
     pcost       dcost       gap    pres   dres   k/t
 0: -1.4000e+01 -6.5333e+01  4e+01  2e-01  2e+00  1e+00
 1: -2.0550e+01 -2.7849e+01  5e+00  3e-02  3e-01  4e-01
 2: -2.1989e+01 -2.2885e+01  5e-01  4e-03  4e-02  7e-02
 3: -2.1999e+01 -2.2012e+01  7e-03  5e-05  6e-04  1e-03
 4: -2.2000e+01 -2.2000e+01  7e-05  5e-07  6e-06  1e-05
 5: -2.2000e+01 -2.2000e+01  7e-07  5e-09  6e-08  1e-07
Optimal solution found.
>>> print(sol['x'])
[ 8.00e+00]
[ 2.00e+00]
[-1.49e-08]

Hodgkin-Huxley model

This is the famous model for simulating the action potential of a neuron.

In terms of neuroscience, the mechanism of an action potential is explained as like this.
(1)When a membrane receives a stimulus, stimulus-dependent Na+ channels open, and the intracellular concentration of Na+ is boosted.
(2) A voltage of the membrane goes up. This is called depolarization.
(3)Instead of closing Na+ channels, K+ channels open, and these decrease the concentration of intracellular ion level while the efflux K+ happen. This is called repolarization Opening K+ channels depend on the intracellular Na+ ion level. This is why relatively a reaction to Na+ is slow.
(4) However, the efflux of K+ induce hypolarization which means the voltage is less than a resting potential. This is because of efflux of K+ too much. In this period, a action potential cannot happen until K+ channels close. This is absolute refractory period.
(5) Finally, the membrane potential goes back to the resting potential. From a view of these process, it is possible that a collision causes the cancellation is explained by refractory period because each action potential has each refractory period. When it meet with another absolute refractory period, an action potential does not react any more. This happens each other.

In the Hodgkin Huxley model, this is shown as well. Activation and inactivation of Na+ channel or K+ is described in n(x, t), m(x, t), h(x, t). These are depends on intensity of the voltage.
These terms in (1) modulate intensity of voltage. Increasing the voltage facilitates K+ and Na+‘s roles through n(x, t) and h(x, t) while boosting the voltage inhibits a Na+‘s role for boosting through m(x, t).
These terms function as the dynamics between Na+ and K+ which induce an action potential. Furthermore, delayed reaction of inhibition of K+ channels is also implemented(Figure 2).
Although chronic stimulus is given, the gap period still exists.
This result suggests that n, m, h terms regulate initiation of action potential, and this was expected to be. Finally, re-opening Na+ channels synchronizes with another action potential, and this is driven by closing K+ channels.

#-*- coding:utf- 8 -*-
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt

class HH(object):
	"""
	simple Hodgkin-Huxley model
	[equation]
	dV/dt = 1/Cm * (Im - g_Na * m^3 * h(V-E_Na)- g_k*n^4(V-E_k)-g_l*(V-E_l))
	g: conductance
	E: potentials
	"""
	#Set HH Parameters
	#HH(1, 120, 36, 0.3, 115, -12, 10.613, 0)
	def __init__(self, Cm, g_Na, g_K, g_l, E_Na, E_k, E_l, Vres):
		self.Cm = Cm #uF/cm2
		self.g_Na = g_Na #mS/cm2
		self.g_K = g_K #mS/cm2
		self.g_l = g_l #mS/cm2
		self.E_Na = E_Na #mV
		self.E_k = E_k #mV
		self.E_l = E_l #mV
		self.Vres = Vres #mV
		
		#function for K channel
		self.alpha_n = np.vectorize(lambda v: 0.01*(-v + 10)/(np.exp((-v + 10)/10) - 1) if v != 10 else 0.1)
		self.beta_n = lambda v: 0.125*np.exp(-v/80)
		self.n_inf = lambda v:self.alpha_n(v)/(self.alpha_n(v) + self.beta_n(v))
		
		#Na channel(activating)
		self.alpha_m = np.vectorize(lambda v: 0.1*(-v + 25)/(np.exp((-v + 25)/10) - 1) if v != 25 else 1)
		self.beta_m = lambda v: 4 * np.exp(-v/18)
		self.m_inf = lambda v: self.alpha_m(v)/(self.alpha_m(v) + self.beta_m(v))
		
		#Na channel(inactivating)
		self.alpha_h = lambda v: 0.07*np.exp(-v/20)
		self.beta_h  = lambda v: 1/(np.exp((-v + 30)/10) + 1)
		self.h_inf   = lambda v: self.alpha_h(v)/(self.alpha_h(v) + self.beta_h(v))
		
	def Ch_act(self, v):
		return self.m_inf(v), self.h_inf(v), self.n_inf(v) 
	
	def stimulus(self, time, st):
		self.I = np.zeros(len(self.time))
		for i, t in enumerate(time):
			#if 5 <= t <= st : self.I[i] = 30#uA/cm2
			if 5 <= t <= st : self.I[i] = 40#uA/cm2
			#if 12 <= t <= 13: self.I[i] = 60 #uA/cm2
		return self.I
		
	def dVdt(self, time, dt, stimulus=30):
		self.time = time
		self.dt = dt
		self.st = stimulus
		Vm = np.zeros(len(time)) #mV
		Vm[0] = self.Vres
		m = self.m_inf(self.Vres)
		h = self.h_inf(self.Vres)
		n = self.n_inf(self.Vres)
		
		M = np.zeros(len(time))
		H = np.zeros(len(time))
		N = np.zeros(len(time))
		
		Na_M = np.zeros(len(time))
		Na_H = np.zeros(len(time))
		K_N = np.zeros(len(time))
		
		I = self.stimulus(time, self.st)
		
		#simulate model
		for i in xrange(1, len(time)):
			Na = self.g_Na * (m**3)*h
			K = self.g_K * (n ** 4)
			L = self.g_l
			
			M[i] = m
			H[i] = h
			N[i] = n
		
			m += (self.alpha_m(Vm[i-1])*(1 - m) - self.beta_m(Vm[i-1])*m) * self.dt
  			h += (self.alpha_h(Vm[i-1])*(1 - h) - self.beta_h(Vm[i-1])*h) * self.dt
  			#h = 0.1
  			n += (self.alpha_n(Vm[i-1])*(1 - n) - self.beta_n(Vm[i-1])*n) * self.dt
  			
  			Na_M[i] = Na
  			K_N[i] = K
  			
  			Vm[i] = Vm[i-1] + (I[i-1] - Na*(Vm[i-1] - self.E_Na) - K * (Vm[i-1] - self.E_k) - L *(Vm[i-1] - self.E_l)) / self.Cm * dt	
  		
  		return Vm, I, Na_M, K_N, M, H, N
  		

	
def main():
	T = 55
	#T = 1000
	dt = 0.025
	time= np.arange(0, T+dt, dt)
	#Cm, g_Na, g_K, g_l, E_Na, E_k, E_l, Vres
	#cell = HH(1, 120, 36, 0.3, 115, -12, 10.613, 0)
	cell = HH(1, 120, 36, 0.3, 115, -12, 10.613, 0)
	v = np.arange(-50, 151)
	a, b, c = cell.Ch_act(v)
	
	plt.figure()
	plt.title('Steady state values of ion channel gating variables');
	plt.ylabel('Magnitude')
	plt.xlabel('Voltage (mV)')
	plt.plot(v, a, v, b, v, c);
	plt.legend( (r'm',r'h',r'n') );
	plt.show();
	
	#Vm, I = cell.dVdt(time, dt, 1000)
	Vm, I, Na_M, K_N, M, H, N= cell.dVdt(time, dt, 55)
	## plot membrane potential trace
	plt.figure()
	plt.plot(time, Vm, time, -30+I, time, Na_M*5, time, K_N*5, time, M*100, time, H*100, time, N*100)
	
	plt.title('Hodgkin-Huxley Example')
	plt.legend( (r'Threshold', r'Stimulus', r'Na_M', r'K_N', r'm', r'h', r'n') )
	plt.ylabel('Membrane Potential (mV)')
	plt.xlabel('Time (msec)')
	plt.show()
	
if __name__ == "__main__":
	main()

Leaky Integrate and Fire model

(When I go back to Tokyo, I would like to comment on this code.)

This code is mostly copied from the end of this site....
Neurdon

#-*- coding:utf-8 -*-
from __future__ import division
import numpy as np
import  matplotlib.pyplot  as plt

class SIF(object):
	"""simple integrate and fire model:
	[equation]
	Cm * dV/dt + gl * (V - EI) = lapp
	Cm=(is the membrane capacitante)
	gl=(is the leak conductance)
	EI=(is the resting membrane potential)
	Iapp=(is the applied current)
	"""
	
	def __init__(self, Cm, gl, EI, Vres=-70):
		self._Cm = Cm #ms
		self._gl = gl #uS
		self._EI = EI #mV
		self._Vres = Vres #mV
		
		#fixed threshold
		self._Vth = - 50 #mV
		
	def dVdt(self, Vmb, Iapp):
		#function for the immediate voltage fluctuation  with Cm, gl, EI and Iapp
		Cm = self._Cm
		gl = self._gl
		EI = self._EI
		return (Iapp - gl * (Vmb - EI)) / Cm
		
	def timecourse(self, current, dt=0.02):
		#function f(x) = f(x - 1) + dt * f'(x)
		voltage = np.empty(len(current))
		
		#initial condition
		voltage[0] = self._Vres
		dt = dt #step siez integration
		
		for i in xrange(1, len(current)):
			dVdt = self.dVdt(Vmb = voltage[i - 1], Iapp = current[i - 1])
			voltage[i] = voltage[i - 1] + dt * dVdt
			#action potential threshold
			if voltage[i] > self._Vth:
				voltage[i - 1] = 0.2 -  np.random.rand() * 5 #overshooting
				voltage[i] = self._EI #back to resting potential
		return voltage

def main():
	dt = .02
	current = int(250/dt)*[0] + int(3000/dt)*[2.5] + int(250/dt)*[0]
	cell = SIF(Cm=4.9, gl=.16, EI=-65, Vres=-75)
	voltage = cell.timecourse(current)
	
	time = np.linspace(0, len(current)*(dt/1000), len(current)) # transform in sec
	plt.plot(time, voltage)
	plt.xlabel("Time (ms)")
	plt.ylabel("Voltage (mV)")
	plt.show()

if __name__ == '__main__':
	main()

Neural Modeling with Python(IF, HH, Izhikevich,Diffusional equation, cable theory etc....)

In this log, I am summarizing Neural Models' articles.
Later, I would like to make gists of them with python codes.

Spiking Neuron Models Single Neurons, Populations, Plasticity
comments for models

STEPS: Modeling and Simulating Complex Reaction-Diffusion Systems with Python
comments for modeling diffusional systems in Neuroscience

neurdon
神経活動を微分方程式で表してみよう
this site is summarizing major models in Neuroscience with python codes.

Modeling Spiking Neurons with Python
the presentation of neurons

NEURON and Python
how to use NEURON with python

Python in Computational Neuroscience & Modular toolkit for Data Processing (MDP)
Explanations for computational models, python, MDP

Python: solving 1D diffusion equation
for resolving how to diffuse molecules

Performance Python: Solving The 2D Diffusion Equation With numpy
this is for 2D diffusion equation.

Finite Difference Heat Equation using NumPy
the method for resolving diffusion equation

Project Euler 89 (python, 0.04s)

#-*- coding:utf-8 -*-
from datetime import datetime

def smallest_num(n):
	counter = 0
	num = n
	dict_num = {"0":0, "1":1, "2":2, "3":3, "4":2, "5":1, "6":2, "7":3, "8":4, "9":2}
	if n >= 1000:
		counter += num / 1000
		num = num % 1000
	for i, d in enumerate(str(num)):
		counter += dict_num[d]
	return counter
	
def count_roma(n,last=0):
	dict_roma = {"I":1, "V":5, "X":10, "L":50, "C":100, "D":500, "M":1000}
	I_list = ["V", "X"]
	X_list = ["L", "C"]
	C_list = ["D", "M"]	
	counter = 0
	length = len(n)
	for i, x in enumerate(n):
		if ((length - 1) == i) and last != 1:
			break
		if ((length - 1) == i) and last == 1:
			counter += dict_roma[x]
			break
		if x ==  "C" and n[i +1] in C_list:
			counter -= 100
		elif x == "X" and n[i +1] in X_list:
			counter -= 10
		elif x == "I" and n[i + 1] in I_list:
			counter -= 1
		else:
			counter += dict_roma[x]
	return counter
		
	
def Euler89(file1):
	before, after = 0, 0
	Number = 0
	for d, i in enumerate(file1):
		if d != 999: Number = count_roma(i)
		else:	Number = count_roma(i, 1)
		length = len(str(i))
		if d != 999: before += length - 1
		else: before += length
		after += smallest_num(Number)
	return before - after

def main():
	f = open("roman.txt", "r")
	file = f.readlines()
	f.close()
	
	start = datetime.now()
	answer = Euler89(file)
	end = datetime.now()
	print end - start, answer
	
if __name__ == "__main__":
	main()

replace() を使うと綺麗だし速い。
敗北感。。。

#こっちだと0.000957s
#-*- coding:utf-8 -*-
from datetime import datetime
def Euler89():
	dict_roma = [["DCCCC", "CM"], ["LXXXX", "XC"], ["VIIII", "IX"], ["IIII", "IV"], ["XXXX", "XL"], ["CCCC" , "CD"]]
	f = open("roman.txt", "r")
	file = f.readlines()
	f.close()
	before = ''
	for x in file:
		before += x
	after = before
	for  x in dict_roma:
		after = after.replace(x[0], x[1])
	return len(before) - len(after)

def main():
	start = datetime.now()
	answer = Euler89()
	end = datetime.now()
	print end - start, answer
	
if __name__ == "__main__":
	main()

Project Euler 70(python, 1.776935s)

#-*- coding:utf-8 -*-
import random
from datetime import datetime

#素数判定
def is_prime3(q,k=50):
	q = abs(q)
	if q == 2: return True
	if q < 2 or q&1 == 0: return False
	d = (q-1)>>1
	while d&1 == 0:
		d >>= 1
	for i in xrange(k):
		a = random.randint(1,q-1)
		t = d
		y = pow(a,t,q)
		while t != q-1 and y != 1 and y != q-1: 
			y = pow(y,2,q)
			t <<= 1
		if y != q-1 and t&1 == 0:
			return False
	return True

	
def permutation(n):
	return sorted(list(str(n)))

def Euler70():
	limit = 10 ** 7
	N = 5000
	dict_list = [x for x in xrange(N) if is_prime3(x)]
	length = len(dict_list)
	dict_phi = []
	for i in reversed(xrange(length)):
		for j in reversed(xrange(i+1)):
			A = dict_list[ i ]
			B = dict_list[ j ]
			n = A * B
			if (n < limit):
				phi = (A - 1) * (B -1)
				if permutation( n ) == permutation( phi ):
					dict_phi.append( [float(n)/float(phi), n])
	return min(dict_phi)[1]

def main():
	start = datetime.now()
	answer = Euler70()
	end = datetime.now()
	print end - start, answer

if __name__ == "__main__":
	main()