Model Development

A new model is created by inheriting from the markov::Model class. This base class defines the virtual pure methods which the new model must implement:

method	purpose
get_identifier()	string identifying the model -- this will be the type='...' value in the input.
get_version()	string identifying the version -- this will be the version='...' value in the input.
configure(...)	method telling qiotoolkit how to read model configuration from input.
calculate_cost(state)	Implements the (parametrized) cost function
calculate_cost_difference(state, transition)	Calculates how much the cost will change if transition is applied.
get_random_state(rng)	Creates a new random starting state
get_random_transition(state, rng)	Proposes a random change to the state
apply_transition(transition, state)	Modifies the state as dictated by the transition

These methods entail the essence of the model: They define the state space and how one can move within it, as well as the association of a cost with each sate. The dynamics of whether to accept or reject a transition are left to the user.

Self consistency

The implementation of calculate_cost_difference() and calculate_cost() must be consistent. That is,

\text{costdiff}(\text{before}, \text{transition}) \equiv \text{cost}(\text{transition}(\text{before})) - \text{cost}(\text{before})

You can use the SelfConsistency Build to verify this holds.

Note

calculate_cost_difference() could be based on a call to calculate_cost() with a modified state, but this is less efficient for models where individual changes affect only a small number of terms.

State and Transition

In some situations it is sufficient to use a base type to represent the transition (or even the state); in others you need to define your own class to hold this information. You may extend markov::State and markov::Transition, respectively, to get the right interfaces. Whichever route you opt for, you need to template the base class of your new model with this two pieces of information:

class MyModel : public ::markov::Model<MyState, MyTransition> {
 public:
  using State_T = MyState;
  using Transition_T = MyTransition;
  ...
};

Note

I find it convenient to typedef these two to State_T and Transition_T to homogenize the interfaces which need to be overloaded.

Graph Model

The special base class ::model::GraphModel makes use of the Graph-Cost correspondency and provides the appropriate interfaces to access the graph structure. If your cost function has the appropriate shape, this can simplify writing the cost function logic substantially.

Example

The following is an example implementation of the above interfaces for a soft-spin Model. It can be found in the cpp/examples/ directoy of the qiotoolkit codebase.

#pragma once

#include "utils/config.h"
#include "markov/state.h"
#include "markov/transition.h"
#include "model/graph_model.h"

namespace examples
{
class SoftSpinState : public ::markov::State
{
 public:
  std::vector<double> spin;
  utils::Structure render() const override { return spin; }
  utils::Structure get_status() const override { return spin; }
  std::string get_class_name() const override { return "SoftSpinState"; }
  static size_t memory_estimate(size_t N)
  {
    return sizeof(SoftSpinState) +
           utils::vector_values_memory_estimate<double>(N);
  }

  static size_t state_only_memory_estimate(size_t N)
  {
    return memory_estimate(N);
  }
};

class SoftSpinTransition : public ::markov::Transition
{
 public:
  SoftSpinTransition() : spin_id(0), new_value(0) {}
  int spin_id;
  double new_value;
  bool operator<(const SoftSpinTransition& trans) const
  {
    if (spin_id == trans.spin_id)
    {
      return new_value < trans.new_value;
    }
    else
    {
      return spin_id < trans.spin_id;
    }
  }
};

class SoftSpin : public ::model::GraphModel<SoftSpinState, SoftSpinTransition>
{
 public:
  using State_T = SoftSpinState;
  using Transition_T = SoftSpinTransition;
  using Graph = ::model::GraphModel<State_T, Transition_T>;

  std::string get_identifier() const override { return "softspin"; }
  std::string get_version() const override { return "0.1"; }

  void configure(const utils::Json& json) override { Graph::configure(json); }

  void configure(Configuration_T& configuration)
  {
    Graph::configure(configuration);
  }

  double calculate_cost(const State_T& state) const override
  {
    double cost = 0;
    for (auto e : edges())
    {
      double term = e.cost();
      for (auto spin_id : e.node_ids())
      {
        term *= state.spin[spin_id];
      }
      cost += term;
    }
    return cost;
  }

  double calculate_cost_difference(
      const State_T& state, const Transition_T& transition) const override
  {
    double diff = 0;
    for (auto edge_id : node(transition.spin_id).edge_ids())
    {
      const auto& e = edge(edge_id);
      double term_before = e.cost();
      double term_after = e.cost();
      for (auto spin_id : e.node_ids())
      {
        term_before *= state.spin[spin_id];
        if (spin_id == transition.spin_id)
        {
          term_after *= transition.new_value;
        }
        else
        {
          term_after *= state.spin[spin_id];
        }
      }
      diff += term_after - term_before;
    }
    return diff;
  }

  State_T get_random_state(utils::RandomGenerator& rng) const override
  {
    State_T state;
    state.spin.resize(nodes().size());
    for (size_t i = 0; i < nodes().size(); i++)
    {
      state.spin[i] = rng.uniform() * 2.0 - 1;
    }
    return state;
  }

  Transition_T get_random_transition(const State_T&,
                                     utils::RandomGenerator& rng) const override
  {
    Transition_T transition;
    transition.spin_id =
        (size_t)floor(rng.uniform() * static_cast<double>(nodes().size()));
    transition.new_value = rng.uniform() * 2.0 - 1;
    return transition;
  }

  void apply_transition(const Transition_T& transition,
                        State_T& state) const override
  {
    state.spin[transition.spin_id] = transition.new_value;
  }

  size_t state_memory_estimate() const override
  {
    return State_T::memory_estimate(nodes().size());
  }

  size_t state_only_memory_estimate() const override
  {
    return State_T::state_only_memory_estimate(nodes().size());
  }
};

}  // namespace examples

template <>
struct std::hash<examples::SoftSpinTransition>
{
  std::size_t operator()(
      const examples::SoftSpinTransition& trans) const noexcept
  {
    return utils::get_combined_hash(trans.spin_id, trans.new_value);
  }
};