Permalink
Browse files

Refactoring of derivative calculation

  • Loading branch information...
1 parent e3f1d96 commit 908d43e3aaa00340cf56c82c9fe4a6ee5ffbd389 @thisrod committed Oct 18, 2011
Showing with 42 additions and 48 deletions.
  1. +32 −31 fermi_hubbard.py
  2. +2 −2 integration.py
  3. +5 −0 kubo_tests.py
  4. +3 −15 physics_tests.py
View
@@ -1,8 +1,9 @@
from numpy.random import randn as normal_deviates
-from numpy import array, mat, diagonal, diagflat as make_diagonal, zeros, identity as unit, logical_or
+from numpy import array, zeros, identity as unit, logical_or
from weightings import Weighting
+from ndmat import ndmat
from adjacencies import *
-from math import sqrt
+from math import sqrt, copysign
from operator import mul
@@ -12,50 +13,50 @@ class FermiHubbardSystem:
"""Parameters: sites, repulsion, hopping, chemical_potential
- The state is represented as a Weighting, whose value is a 2x[sites]x[sites] array of the up and down greens functions. Here [sites] is the dimensions of the grid."""
+ The state is represented as a Weighting, whose value is a 2x[sites]x[sites] array of the up and down greens functions. Here [sites] is the dimensions of the grid, stored as a tuple."""
def __init__(self, model = None, **parameters):
if model is not None:
self.__dict__ = model.__dict__.copy()
self.__dict__.update(parameters)
-
- def delta(self, normal_greens, spin, noise):
- def repulsion(self, site):
- return self.repulsion * greens[(1-spin,)+2*site] \
- - abs(self.repulsion) * (0.5 - greens[(spin,)+2*site])
-
- sites = normal_greens.shape[1]
- if spin == 1:
- f = 1
- else:
- f = -self.repulsion/abs(self.repulsion)
- return \
- make_diagonal([self.hopping]*(sites-1), 1) \
- + make_diagonal([self.hopping]*(sites-1), -1) \
- - make_diagonal( \
- self.repulsion_terms(normal_greens, spin) \
- - self.chemical_potential \
- + f * sqrt(2*abs(self.repulsion)) * noise)
+ self.sites = tuple(self.sites)
def derivative(self, time, state, noise):
greens = state.mean
+ noise = array(noise).reshape((2,)+self.sites)
def greens_dot(spin):
- n = normal_greens[spin,:,:]
- noise1 = noise[:len(noise)//2]
- noise2 = noise[-(len(noise)//2):]
- particles = mat(n)
- holes = mat(unit_like(n) - n)
- delta1 = mat(self.delta(normal_greens, spin, noise1))
- delta2 = mat(self.delta(normal_greens, spin, noise2))
- return 0.5*(holes*delta1*particles + particles*delta2*holes)
+ n = greens[spin,::]
+
+ particles = ndmat(n)
+ holes = ndmat(unit_like(n) - n)
+
+ return 0.5*(holes*delta(spin, 0)*particles + particles*delta(spin, 1)*holes)
+ def delta(spin, r):
+
+ def repulsion(site):
+ return self.repulsion * greens[(1-spin,)+2*site] \
+ - abs(self.repulsion) * (0.5 - greens[(spin,)+2*site])
+
+ f = 1 if spin == 1 else copysign(1, -self.repulsion)
+
+ result = ndmat(zeros(2*self.sites))
+ for i in sites(self.sites):
+ result[2*i] = repulsion(i) \
+ + f * sqrt(2*abs(self.repulsion)) * noise[(r,)+i] \
+ - self.chemical_potential
+ for i, j in adjacencies(self.sites):
+ assert(i != j)
+ result[i+j] = self.hopping
+ return result
+
weight_log_dot = \
self.hopping * sum([greens[(0,)+i+j] + greens[(1,)+i+j] for i, j in adjacencies(self.sites)]) \
+ self.repulsion * (greens[0,::]*greens[1,::]).sum() \
- - self.chemical_potential * sum([greens[(0,)+i+i] + [greens[(1,)+i+i] for i in sites(self.sites)])
+ - self.chemical_potential * sum([ greens[(0,)+i+i] + greens[(1,)+i+i] for i in sites(self.sites)])
return Weighting( \
array([greens_dot(0), greens_dot(1)]), \
@@ -72,7 +73,7 @@ def initial(self):
"Answer the state at infinite temperature"
id = zeros(self.sites*2)
for i in sites(self.sites):
- id[i*2] = 1
+ id[i*2] = 0.5
return Weighting(array([id, id]))
View
@@ -12,14 +12,14 @@ def __init__(self, a_simulation, a_noise_source, timestep):
self.timestep = float(timestep)
self.noise = a_noise_source
- def integrate(self, initial_state, duration, record, **run_labels):
+ def integrate(self, initial_state, duration, record, *run_labels):
"""stochastically integrate a single sample, starting at time 0"""
t = 0
state = initial_state
next_sample_time = 0
while t <= duration:
if t > (next_sample_time - 0.5*self.timestep):
- record[t] = self.system.moments(state)
+ record[(t,) + run_labels] = self.system.moments(state)
next_sample_time = record.after(next_sample_time)
t, state = t + self.timestep, state + self.increment(t, state)
View
@@ -39,11 +39,16 @@ def setUp(self):
self.integrator.integrate(self.system.initial(), 3.1, self.moments)
def testSolution(self):
+ self.integrator.integrate(self.system.initial(), 3.1, self.moments)
for t in range(3):
exact = self.system.exact(t, self.noise)
computed = self.moments(t)
self.assertTrue(abs(computed - exact) < 0.05*abs(exact))
+ def testLabels(self):
+ self.integrator.integrate(self.system.initial(), 3.1, self.moments, 1)
+ self.assertTrue((1,) in self.moments.results)
+
if __name__ == '__main__':
run_tests()
View
@@ -5,6 +5,7 @@
from numpy import array, mat, diagonal, diagflat as make_diagonal, zeros, identity as unit, logical_or, isnan
from unittest import TestCase, main as run_tests
from pairs import Pair
+from weightings import Weighting
# Test data
@@ -83,19 +84,6 @@ def testOverlap(self):
for sltr, greens, spin in specimens():
if logical_or(greens[0,:,:] == 0, greens[1,:,:] == 0).all():
self.assertTrue((sltr.repulsion_terms(greens, spin) == 0.5).all())
-
-
-class TestDelta(TestCase):
-
- def testShape(self):
- for sltr, greens, spin in specimens():
- for noise in noise_specimens(greens.shape[1]):
- self.assertTrue(sltr.delta(greens, spin, noise[:noise.size//2]).shape == greens[spin,:,:].shape)
-
- def testTridiagonal(self):
- for sltr, greens, spin in specimens():
- for noise in noise_specimens(greens.shape[1]):
- self.assertTrue(is_tridiagonal(sltr.delta(greens, spin, noise[:noise.size//2])))
class TestGreensDerivative(TestCase):
@@ -126,8 +114,8 @@ def testDerivation(self):
class TestFermiHubbardSystemInterface(TestCase):
def setUp(self):
- self.sltr = FermiHubbardSystem(sites = [2], repulsion = 1, chemical_potential = 0.5, hopping = 2)
- self.state = Pair(1, array([unit(5), unit(5)]))
+ self.sltr = FermiHubbardSystem(sites = [5], repulsion = 1, chemical_potential = 0.5, hopping = 2)
+ self.state = Weighting(array([unit(5), unit(5)]))
def testNoise(self):
self.assertTrue(len(self.sltr.noise_required(self.state)) == 10)

0 comments on commit 908d43e

Please sign in to comment.