Java tutorial
/* * Copyright 2013 University of Glasgow. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ package broadwick.stochastic; import com.google.common.base.Throwables; import com.google.common.collect.Table; import com.google.common.collect.TreeBasedTable; import java.util.Collection; import java.util.Collections; import java.util.HashMap; import java.util.HashSet; import java.util.Map; import java.util.Set; import lombok.Getter; import lombok.Setter; import lombok.Synchronized; import lombok.extern.slf4j.Slf4j; import org.apache.commons.lang3.time.StopWatch; /** * Base class for stochastic simulators. Extending classes just have to implement * <code>performStep</code> which gets invoked as long as the simulation is running. They also may override * <code>init</code> but should in that case call * <code>super.init()</code> to avoid unwanted effects. <p> Simulators also should handle the special time * <code>theta</code>: is is a moment in time when * <code>doThetaEvent</code> is supposed to be invoked (e.g. to measure species populations at this moment). Consider * this, if you want to implement a simulator. * <p/> * The <TT>doEvent</TT> are supposed to be invoked when a simulator causes a reaction to fire. * <p/> * If an extending class sticks to these conventions, it can take full benefit of the observers system: One or more * {@link Observer} can be registered at a simulator and observe certain aspects of the simulation (see the * {@link Observer}s javadoc for more information). * <p/> * Note that this class is NOT thread-safe. Consequently, if a simulation program uses multiple threads, it should * acquire a lock on a simulator (using a <TT>synchronized</TT> block) before accessing its state. Note however, that * one can launch many simulations in parallel with as many threads, as long as <SPAN CLASS="textit">each thread has its * own</SPAN> <TT>Simulator</TT>. * <p/> * This class is mostly copied from the FERN projects stochastic classes and modified to fit into this framework. * @see http://www.bio.ifi.lmu.de/FERN/ */ @Slf4j public abstract class StochasticSimulator { /** * Creates a new simulator using a given amount manager. * @param amountManager the amount manager used in the simulator. * @param transitionKernel the transition kernel to be used with the stochastic solver. */ public StochasticSimulator(final AmountManager amountManager, final TransitionKernel transitionKernel) { this(amountManager, transitionKernel, false); } /** * Creates a new simulator. * @param amountManager the amount manager used in the simulator. * @param transitionKernel the transition kernel to be used with the stochastic solver. * @param reverseTime true if we wish to go backwards in time. */ public StochasticSimulator(final AmountManager amountManager, final TransitionKernel transitionKernel, final boolean reverseTime) { this.amountManager = amountManager; this.transitionKernel = transitionKernel; this.reverseTime = reverseTime; this.thetaQueue = new ThetaQueue(); } /** * Starts the simulation up to a given time. It just uses a {@link DefaultController} and calls * {@link StochasticSimulator#start(SimulationController)}. * @param time simulation time */ public final void run(final double time) { if (controller == null) { controller = new DefaultController(time); } run(); } /** * Start the simulation with the given {@link SimulationController}. */ public final void run() { try { // If we haven't set a controller for the process then set the default one with a max time of 5 // (unknonw units). if (controller == null) { controller = new DefaultController(5); } init(); // tell our observers that the simulation has started for (Observer observer : observers) { observer.started(); } final StopWatch sw = new StopWatch(); sw.start(); while (controller.goOn(this)) { // The performStep() method is implemented by the specifc stochastic algorithm (e.g. Gillespie's) performStep(); //amountManager.doHouseKeeping(); // tell our observers that the step has been completed. for (Observer observer : observers) { observer.step(); } sw.split(); log.trace("Performed step in {}", sw.toSplitString()); } // tell our observers that the simulation has finished. for (Observer observer : observers) { observer.finished(); } } catch (Exception e) { log.error("Error running stochastic simulation. {}", Throwables.getStackTraceAsString(e)); throw e; } } /** * Initializes the algorithm. <ul><li>set currentTime=0</li><li>reset the {@link AmountManager}</li><li>recalculate * the propensities</li></ul> Gets called at the very beginning of * <code>start</code> */ public final void init() { currentTime = startTime; } /** * Fires a reaction. It calls the observers and the {@link AmountManager}. * @param mu reaction to be fired. * @param t time at which the firing occurs. */ protected final void doEvent(final SimulationEvent mu, final double t) { doEvent(mu, t, 1); } /** * Perform an event by informing the subscribed observers and the amountManager. * @param event event to be performed. * @param t time at which the firing occurs. * @param n the number of times mu is to be fired. */ protected final void doEvent(final SimulationEvent event, final double t, final int n) { log.trace("Performing event {}", event); for (Observer observer : observers) { observer.observeEvent(event, t, n); } // change the amount of the reactants if (!Double.isInfinite(t)) { amountManager.performEvent(event, n); } } /** * Gets called when the simulator reaches the predetermined time of a theta event. All the observers for the events * that are configured for this time are notified and given a list of events that are triggered. */ protected final void doThetaEvent() { final double nextThetaEventTime = thetaQueue.getNextThetaEventTime(); log.trace("Performing theta events at t = {}", nextThetaEventTime); final Map<Observer, Collection<Object>> nextEvents = thetaQueue.getNextEventDataAndRemove(); for (Map.Entry<Observer, Collection<Object>> entry : nextEvents.entrySet()) { final Observer observer = entry.getKey(); final Collection<Object> events = entry.getValue(); if (events != null) { log.trace("Informing observer of {} theta events", events.size()); observer.theta(nextThetaEventTime, events); } } log.trace("Finished theta events next = {}", thetaQueue.getNextThetaEventTime()); } /** * Add an observer to the engines list of observers. * @param observer the observer. */ public final void addObserver(final Observer observer) { if (observer.getProcess() != this) { throw new IllegalArgumentException("Observer doesn't belong to this simulator!"); } this.getObservers().add(observer); } /** * Theta defines a moment, where the simulator has to invoke <TT>theta</TT> of a observers. It is used e.g. to * determine the amounts of species at one moments. Extending class just have to call * {@link Simulator#doThetaEvent()} which basically calls the observers. * @return the theta */ public final double getNextThetaEventTime() { return thetaQueue.getNextThetaEventTime(); } /** * Register a new theta event that is triggered at a given time. * @param obs the observers which is registering. * @param thetaTime the time the theta event occurs. * @param event the theta event. */ public final void registerNewTheta(final Observer obs, final double thetaTime, final Object event) { thetaQueue.pushTheta(thetaTime, obs, event); } /** * Get the default observer for the stochastic process. The default observer is first one defined. * @return the default observer. */ public final Observer getDefaultObserver() { return observers.toArray(new Observer[observers.size()])[0]; } /** * Manages the registered theta events. Each registered theta event is stored in a table containing the time of the * event the list of observers for that event and the list of events for that time. */ private class ThetaQueue { /** * Construct an empty theta queue. */ public ThetaQueue() { thetas = TreeBasedTable.create(); if (reverseTime) { nextThetaEventTime = Double.NEGATIVE_INFINITY; } else { nextThetaEventTime = Double.POSITIVE_INFINITY; } } /** * Add a new theta event to the registry, including the time, collection of observers and collection of events. * Each event is stored as an object where it is assumed that the observer * @param time the time the theta event occurs. * @param obs the observers. * @param theta the theta event. */ @Synchronized public void pushTheta(final double time, final Observer obs, final Object theta) { log.trace("Adding new {} theta event at t={}", theta.getClass(), time); Collection<Object> events = thetas.get(time, obs); if (events == null) { try { log.trace("No theta events registered at t={}; Creating new list", time); events = new HashSet<>(); events.add(theta); thetas.put(time, obs, events); } catch (UnsupportedOperationException | ClassCastException | IllegalArgumentException | IllegalStateException e) { log.error("Could not register theta. {}", e.getLocalizedMessage()); } } else { log.trace("Found {} theta events for t={}; Adding new event", events.size(), time); events.add(theta); } if (reverseTime) { nextThetaEventTime = Math.max(nextThetaEventTime, time); } else { nextThetaEventTime = Math.min(nextThetaEventTime, time); } } /** * Get the next observer and the collection of events they are subscribed to and remove it from the theta queue. * @return a map of the observers and their subscribed data. */ public Map<Observer, Collection<Object>> getNextEventDataAndRemove() { final Map<Observer, Collection<Object>> nextEventData = thetas.row(nextThetaEventTime); // we have a view of the underlying data that we want to return, copy it first then delete the // underlying data. final Map<Observer, Collection<Object>> eventData = new HashMap<>(nextEventData); log.trace("Found {} configured events and observers at t={}", eventData.size(), nextThetaEventTime); for (Observer obs : eventData.keySet()) { final Collection<Object> removed = thetas.remove(nextThetaEventTime, obs); log.trace("Removed {} items from registered theta list", removed.size()); } // Now we need to update the nextThetaEventTime if (reverseTime) { if (thetas.rowKeySet().isEmpty()) { nextThetaEventTime = Double.NEGATIVE_INFINITY; } else { nextThetaEventTime = Collections.max(thetas.rowKeySet()); } } else { if (thetas.rowKeySet().isEmpty()) { nextThetaEventTime = Double.POSITIVE_INFINITY; } else { nextThetaEventTime = Collections.min(thetas.rowKeySet()); } } return eventData; } @Getter private final Table<Double, Observer, Collection<Object>> thetas; @Getter private double nextThetaEventTime; } /** * Reset propensities when a event has been executed. */ public abstract void reinitialize(); /** * Performs one simulation step. Between steps the terminating condition is checked. The simulators controller is * passed, if an extending class wants to check it within one step. Implementing algorithms cannot be sure that the * propensities are correct at the beginning of a step (e.g. if the volume changed). They should override * {@link Simulator#setAmount(int, long)} and {@link Simulator#setVolume(double)} if they need correct values! */ public abstract void performStep(); /** * Gets the name of the algorithm. * @return name of the algorithm */ public abstract String getName(); /** * Seed the RNG if required. * @param seed the RNG seed. */ public abstract void setRngSeed(final int seed); @Setter private int startTime = 0; @Setter @Getter @SuppressWarnings("PMD.UnusedPrivateField") private TransitionKernel transitionKernel; @Getter private AmountManager amountManager; @Setter @Getter @SuppressWarnings("PMD.UnusedPrivateField") private double currentTime = 0; @Getter @Setter private SimulationController controller = null; private ThetaQueue thetaQueue; @Getter private final Set<Observer> observers = new HashSet<>(1); protected boolean reverseTime = false; }