Jump to content
  • Revista PROGRAMAR: Já está disponível a edição #60 da revista programar. Faz já o download aqui!

JoaoRodrigues

Steepest Descent Algorithm - Maneiras de optimizar?

Recommended Posts

JoaoRodrigues

Caros,

Já lá vai um tempo que posto aqui algo. Andei ocupado a escrever um pequeno script que dado um número de pontos, optimiza a distribuição dos mesmo numa esfera de forma a que a distância entre pontos adjacentes seja sempre a mesma.

A ideia é primeiro gerar os pontos na esfera (uso o Golden Spiral Algorithm para isso) e depois minimizar uma função (E=1/r, onde r é a distância entre 2 pontos) até à convergência. A minimização fica a cargo de um algoritmo chamado Steepest Descent. Basicamente, a cada ciclo da minimização, são calculadas forças a actuar em cada ponto (força é a derivada da energia numa determinada direcção, x y ou z) e depois move-se o ponto de acordo com estas forças.

Isto funciona bastante bem. Para 4 pontos temos um belo tetraedro por exemplo, o que é a distribuição óptima. Para os restantes também funciona porreiro. E é rápido, para poucos pontos (digamos, N < 20 ). Não preciso de usar isto para N maior que 20 sinceramente, ou os casos são (muito) raros. Mas estou curioso em ver se consigo optimizar o código, tirando partido de mais funções de Numpy ou Scipy (sem ser a scipy.optimizer obviamente..).

Ideias? Se for preciso dou explicações mais detalhadas sobre cada passo.

O código está aqui: https://gist.github.com/4698472

Abraços

João

Share this post


Link to post
Share on other sites
ffunenga

Deixo aqui duas possiveis optimizações que eu tentaria fazer à função calc_forces:

def calc_forces(system):

forces = [ [0,0,0] for particle in system ]
# Em vez de [0,0,0] utilizava arrays de numpy de forma a poder melhorar a
# operação que ocorre em cada iteração dos ciclos aninhados abaixo:
#	 forces = [ np.zeros(3) for particle in system ]
nparticles = len(system)

# Cautela com estes dois ciclos aninhados. Para *nparticles* elevado
# (>1000) isto é lento.
for i in xrange(nparticles):
	#(...)
	for j in xrange(i+1, nparticles):
		#(...)
		# As próximas linhas resultam da alteração feita na lista forces
		r = particle_i-particle_j
		f = numpy.linalg.norm(r)**-3
		__t = r*f
		forces[i] += __t
		forces[j] -= __t

def calc_forces(system):
n = len(system)
forces = np.array([np.zeros(3) for particle in range(n)])

idxs = np.arange(n)
idxi,idxj = np.meshgrid(idxs, idxs.T)
__triu = np.triu_indices(n, 1)
i = idxj[__triu]
j = (idxi-idxj)[__triu]

r = system[i]-system[j]
__t = np.array([r[i]*np.linalg.norm(ri)**-3 for ri in r])

#(...) aplicar na array forces as somas/subtrações nos indices i e j

return forces

Share this post


Link to post
Share on other sites
motherFFH

Falta-me o scipy para correr o código e também não conheço o algoritmo. Tendo em conta isto, a única optimização que posso sugerir é remover algumas invariantes dos loops: o r*f como dito acima, e outra que parece importante, optimizar o steepest_descent() para remover cálculos duplicados.

Nessa função os cálculos de #3.1 e # 3.2 só variam se o valor de coordinates tiver mudado, o que só acontece se se entrar em #3.3 ou mais tarde em coordinates = old_coordinates. A estrutura do código pode ser revista para estes valores serem apenas recalculados se o valor de coordinates tiver mudado e também aproveitar os valores previamente calculados de old_coordinates. Também o "new_energy = calc_coulomb(coordinates)" logo em #3.4 devia ser feito à saída do if do #3.3 e antes ser inicializado a cur_energy.

Concretizando um pouco:

def steepest_descent(coordinates, max_steps, step_size, min_step_size, min_energy_delta, verbose=True):
#(...)
from collections import deque
## queue com duas posições apenas: old, cur. item da queue: tuple(energy, forces, force_normal)
results_cache = deque([], 2)
#(...)
# dentro do for, usar os resultados da cache

	if len(results_cache) > 0:
		t = results_cache[-1]
		cur_energy = t[0]
		forces = t[1]
		force_normal = t[2]
		results_cache.append(t) # simula old_collections = collections a nível de resultados guardados
	else:
		# 3.1 Calculate energy
		cur_energy = calc_coulomb(old_coordinates)

		# 3.2 Calculate forces
		forces = calc_forces(coordinates)
		force_normal = numpy.linalg.norm(forces)
		results_cache.append((cur_energy, forces, force_normal))
#(...)
	# 3.3 Move particles
	new_energy = cur_energy ## iniciar aqui new_energy
	if force_normal > 0:
		for p in xrange(nparticles):
			directions = (forces[p]/force_normal)*step_size
			coordinates[p] += directions
			coordinates[p] /= numpy.linalg.norm(coordinates[p])
		new_energy = calc_coulomb(coordinates) ## só calcula nova se coordinates tiverem mudado
		results_cache.clear() ## e faz reset da cache nesse caso
#(...)
	# 3.4 Check new energy and accept move if favourable + accelerate in that direction
	if new_energy >= cur_energy:
		step_size /= 2
		if step_size < min_step_size:
			print "-- Step Size Converged (%1.5f < %1.5f)" %(step_size, min_step_size)
			break
		coordinates = old_coordinates
		## se coordinates reverterem para o valor original, anula a cache referente a coordinates
		## em results_cache ficará a cache de old_coordinates, caso exista...
		if len(results_cache) > 0:
			results_cache.pop()

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

×

Important Information

By using this site you accept our Terms of Use and Privacy Policy. We have placed cookies on your device to help make this website better. You can adjust your cookie settings, otherwise we'll assume you're okay to continue.