Browse files

more refactoring

  • Loading branch information...
1 parent c2ae8f5 commit 7b473fa8b3c546396d48f300bf10fe488a0c5183 @thisrod committed Nov 23, 2011
Showing with 30 additions and 189 deletions.
  1. +3 −34 README.markdown
  2. +25 −91 integration.py
  3. +2 −6 physics/kubo.py
  4. +0 −58 weightings.py
View
37 README.markdown
@@ -3,9 +3,6 @@ The point of this is explained in [the paper](http://arxiv.org/abs/quant-ph/0404
Review requests
---------------
-* There must be a better way to write Simulator.__init__
-* Please suggest extra Greens functions to test
-* TestDelta and TestDerivative aren't nearly thorough enough
* Please confirm the assumption 〈H〉 = -t∑<sub>〈ij〉σ</sub> n<sub>ijσ</sub> + U∑<sub>j</sub> n<sub>jj↑</sub> n<sub>jj↓</sub> - μ ∑<sub>jσ</sub>n<sub>jjσ</sub>
Issues
@@ -30,37 +27,9 @@ The n matrices for spin up and spin down are stored in a single array N. The sp
Design notes
-----
-Integration is carried out by four classes. Three of them, <tt>Integrator</tt>, <tt>Noise</tt>, and <tt>Moments</tt>, do the integration. A forth class represents the system of interest, and is different for every integration. Let's call it <tt>Physics</tt>.
+The class <tt>physicalSystem</tt> has subclasses such as <tt>KuboOscillator</tt> and <tt>FermiHubbardGrid</tt>. These remember the physical parameters of the system, and calculate exact solutions for moments, if any are known.
-<tt>Integrator</tt> is a method object, that integrates the SDEs. It
+The class<tt>stateEnsemble</tt> holds a numerical representation of system states, and context such as time and random processes. Subclasses compute derivatives and moments.
-* stores numerical parameters, most notably the timestep,
-* adds up the increments of the system state at each timestep,
-* passes moments to the <tt>Record</tt> at the times it requests.
+<tt>record</tt> stores the results of simulations, and tells the integration routines when to record them.
-Subclasses implement particular algorithms, by the method <tt>increment</tt> which
-
-* computes the increment of the state over a specified timestep.
-
-<tt>Record</tt> averages moments over trajectories. It
-
-* tells what times to record moments
-* records the moments of a labelled trajectory at a given time
-* interpolates moments at times between those recorded
-* computes a weighted average of moments over the recorded trajectories
-
-<tt>Noise</tt> generates Wiener processes. It
-
-* generates an independent process for every trajectory label and integer index
-* computes the value of each process at a given time
-* computes the derivative of each process over any timestep
-
-<tt>Physics</tt> is the object that knows what's going on. It
-
-* remembers physical parameters
-* constructs initial states
-* tells which noise processes are required to compute the derivative of a state
-* computes the derivative of a state, given the noise
-* computes the moments to be recorded for a state.
-
-Currently, <tt>Integrator</tt> associates <tt>Physics</tt> and <tt>Noise</tt>, and passes arrays of numbers between them. Things would be simpler if this association was direct, and <tt>Integrator</tt> was ignorant of <tt>Noise</tt>. This would require the <tt>derivative</tt> method to take labels.
View
116 integration.py
@@ -1,100 +1,34 @@
-"""
-Semi implicit integration, and storage of results.
-"""
-from namespace import *
-from collections import MutableMapping, Callable
-from InterpoList import InterpoList as Interpolation
-from weightings import Weighting
+def integrate_exactly(state, duration, record):
+ "record.after rounds up, record.next assumes time is close to recording time"
+ state = state.advanced_exactly(record.after(state))
+ while record.next(state) <= end:
+ record.add(state)
+ state = state.advanced_exactly(record.next(state))
+ if state.time <= end:
+ record.add(state)
+
-# StepwiseIntegrator computes steps with the increment method, which must be overridden by subclasses.
+class stepwise_integrator(object):
-class Integrator:
+ def __init__(self, timestep):
+ self.step = timestep
- def __init__(self, a_simulation, a_noise_source, timestep):
- self.system = a_simulation
- self.timestep = float(timestep)
- self.noise = a_noise_source
-
- def integrate(self, initial_state, duration, record, *run_labels):
- """stochastically integrate a single sample, starting at time 0"""
- self.noise.set_labels(run_labels)
- t = 0
- state = initial_state
- next_sample_time = 0
- while t <= duration:
- if t > (next_sample_time - 0.5*self.timestep):
- record[(t,) + run_labels] = (state.weight, self.system.moments(state))
- next_sample_time = record.after(next_sample_time)
- t, state = t + self.timestep, state + self.increment(t, state)
+ def __call__(self, state, duration, record)
+ final_time = state.time + duration
+ sample_time = record.after(state.time)
+ while state.time < final_time:
+ while state.time < min(final_time, sample_time):
+ state = self.advance(state)
+ if state.time >= sample_time:
+ record.add(state)
+ sample_time = record.next(state.time)
-class SemiImplicitIntegrator(Integrator):
+class semi_implicit_integrator(stepwise_integrator):
- def increment(self, t, state):
+ def advance(self, state):
halfstep = state
- xis = array([self.noise[i](t, self.timestep) for i in self.system.noise_required(state)])
for i in range(4):
- halfstep = state + 0.5*self.timestep*self.system.derivative(t, halfstep, xis)
- return self.timestep*self.system.derivative(t, halfstep, xis)
-
-
-
-class Record(MutableMapping, Callable):
-
- """ I average the results of trial runs """
-
- # results is a dictionary of run_labels
- # results[run_label] is an Interpolation from times to Weightings
- # moments need to interpret + and * as elemental operations. I.e., ndarrays are in, but lists and tuples are out.
-
- # self[t] is a tuple of weight and moments. self[t] = moments implies the weight is 1.0
-
- # warning: Python ad-hocity: x[2] means x.__getitem__(2), but x[2,3] means x.__getitem__((2,3))
-
- def __init__(self, timestep):
- self.timestep = timestep
- self.results = {}
-
- def demux_key(self, key):
- # Split key into time and run_label
- if type(key) is tuple:
- return float(key[0]), key[1:]
- else:
- return float(key), ()
-
- def demux_value(self, value):
- # Split value into weight and moments
- # Storing moments in tuples ist verboten: see class comment.
- if type(value) is tuple:
- return float(value[0]), value[1]
- else:
- return 1.0, value
-
- def __call__(self, key):
- time, prefix = self.demux_key(key)
- total = Weighting(0, 0)
- for x in [self.results[label](time) for label in self.results if label[0:len(prefix)]==prefix]:
- total = total.combine(x)
- return total.mean
-
- def __delitem__(self, key): pass
-
- def __iter__(self): pass
-
- def __len__(self): pass # Presumably answer the number of items stored?
-
- def __setitem__(self, key, value):
- time, run_label = self.demux_key(key)
- weight, moments = self.demux_value(value)
- if run_label not in self.results:
- self.results[run_label] = Interpolation()
- self.results[run_label][time] = Weighting(moments, weight)
-
- def __getitem__(self, key):
- time, run_label = self.demux_key(key)
- w = self.results[run_label][time]
- return w.weight, w.mean
-
- def after(self, time):
- return time + self.timestep
+ halfstep = state.advanced(0.5*self.step, halfstep)
+ return state.advanced(self.step, halfstep)
View
8 physics/kubo.py
@@ -1,9 +1,8 @@
"""System and representation for the Kubo oscillator. To be made parallel, as a CUDA exercise."""
-from namespace import *
-from physics import *
+from dynamics import *
-class KuboOscillator(physicalSystem):
+class KuboOscillator(object):
def __init__(self, resonance_frequency):
self.w = resonance_frequency
@@ -22,9 +21,6 @@ def set_amplitude(self, x0):
def derivative(self, noise):
return 1j*self.representations*(self.system.w + noise[0:self.size])
- def weight_log_derivative(self, noise):
- return zeros_like(self.weights)
-
def advanced_exactly(self, duration):
final = copy(self)
final.time += duration
View
58 weightings.py
@@ -1,58 +0,0 @@
-from namespace import *
-
-class Weightings(object):
-
- """Stores an ndarray of values, and one of weights. The weights correspond to slices of the values along the leading dimension."""
-
- def __init__(self, values, weights = None):
- assert values.shape[0:1] == weights.shape
- self.means = values
- self.weights = weights
-
- def __repr__(self):
- return "Weightings(%s, weights = %s)" % (repr(self.means), repr(self.weights))
-
- def __eq__(self, other):
- # We care about permutations because of elemental arithmetic
- return self.means == other.means and self.weights == other.weights
-
- def __add__(self, other):
- return Weighting(self.means + other.means, self.weights + other.weights)
-
- def __mul__(self, scalar):
- return Weighting(self.means*scalar, self.weights*scalar)
-
- def __rmul__(self, scalar):
- return Weighting(scalar*self.means, scalar*self.weights)
-
- def __div__(self, scalar):
- return Weighting(self.means/scalar, self.weights/scalar)
-
- def combine(self, other):
- # Fails on zero weights. Note the argument at http://alvyray.com/Memos/MemosCG.htm#ImageCompositing, re alpha channels. Zero weight implies the value doesn't matter, so the round-trip error is spurious.
- broadcast = weights.shape + (len(self.means.shape)-1) * (1,)
- ws = self.weights.copy().reshape(broadcast)
- wo = other.weights.copy().reshape(broadcast)
- weights = ws + wo
- means = (self.means*ws + other.means*wo) / weights
- return Weighting(means, weights.reshape(self.weights))
-
-
-class WeightingIncrement(object):
-
- """An increment to a Weighting, stored as an absolute change to the mean, and a relative change to the weight. When called with a Weighting as argument, it returns an Weighting to be added."""
-
- def __init__(self, values, weights):
- assert values.shape[0:1] == weights.shape
- self.means = values
- self.weights = weights
-
- def __eq__(self, other):
- return self.means == other.means and self.weights == other.weights
-
- def __call__(self, ws):
- return Weightings(self.means + ws.means, (1+self.weights)*ws.weights)
-
-
-def weightclose(self, other):
- return allclose(self.means, other.means) and allclose(self.weights, other.weights)

0 comments on commit 7b473fa

Please sign in to comment.