{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. Introduction of Rule-based Modeling\n",
    "\n",
    "E-Cell4 provides the rule-based modeling environment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from ecell4.prelude import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.1. count_species_matches\n",
    "\n",
    "First, `count_species_matches` counts the number of matches between `Species`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "2\n"
     ]
    }
   ],
   "source": [
    "sp1 = Species(\"A(b^1).B(b^1)\")\n",
    "sp2 = Species(\"A(b^1).A(b^1)\")\n",
    "pttrn1 = Species(\"A\")\n",
    "print(count_species_matches(pttrn1, sp1))  # => 1\n",
    "print(count_species_matches(pttrn1, sp2))  # => 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above case, `Species.count` just returns the number of `UnitSpecies` named `A` in `Species` regardless of its sites. To specify the occupancy of a bond:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "pttrn1 = Species(\"A(b)\")\n",
    "pttrn2 = Species(\"A(b^_)\")\n",
    "print(count_species_matches(pttrn1, sp1))  # => 0\n",
    "print(count_species_matches(pttrn2, sp1))  # => 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where `A(b)` suggests that bond `b` is empty, and `A(b^_)` does that bond `b` is occupied. Underscore `_` means wildcard here. Similarly, you can also specify the state of sites."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "sp1 = Species(\"A(b=u)\")\n",
    "pttrn1 = Species(\"A(b)\")\n",
    "pttrn2 = Species(\"A(b=u)\")\n",
    "print(count_species_matches(pttrn1, sp1))  # => 1\n",
    "print(count_species_matches(pttrn2, sp1))  # => 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`A(b)` says nothing about the state, but `A(b=u)` specifies both state and bond. `A(b=u)` means that `UnitSpecies` named `A` has a site named `b` which state is `u` and the bond is empty. Wildcard `_` is acceptable even in a state and name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "2\n"
     ]
    }
   ],
   "source": [
    "sp1 = Species(\"A(b=u^1).B(b=p^1)\")\n",
    "pttrn1 = Species(\"A(b=_^_)\")  # This is equivalent to `A(b^_)` here\n",
    "pttrn2 = Species(\"_(b^_)\")\n",
    "print(count_species_matches(pttrn1, sp1))  # => 1\n",
    "print(count_species_matches(pttrn2, sp1))  # => 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wildcard `_` matches anything, and the pattern matched is not consistent between wildcards even in the `Species`. On the other hand, named wildcards, `_1`, `_2` and so on, confer the consistency within the match."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "sp1 = Species(\"A(b^1).B(b^1)\")\n",
    "pttrn1 = Species(\"_._\")\n",
    "pttrn2 = Species(\"_1._1\")\n",
    "print(count_species_matches(pttrn1, sp1))  # => 1\n",
    "print(count_species_matches(pttrn2, sp1))  # => 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where the first pattern matches in two ways (`A.B` and `B.A`), but the second matches nothing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n"
     ]
    }
   ],
   "source": [
    "sp1 = Species(\"A(b^1).A(b^1)\")\n",
    "pttrn1 = Species(\"_1._1\")\n",
    "print(count_species_matches(pttrn1, sp1))  # => 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.2. ReactionRule.count and generate\n",
    "\n",
    "`ReactionRule` also has a function to count matches agaist the given list of reactants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n"
     ]
    }
   ],
   "source": [
    "rr1 = create_unimolecular_reaction_rule(Species(\"A(p=u)\"), Species(\"A(p=p)\"), 1.0)\n",
    "sp1 = Species(\"A(b^1,p=u).B(b^1)\")\n",
    "print(rr1.count([sp1]))  # => 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`ReactionRule.generate` returns a list of `ReactionRule`s generated based on the matches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['A(b^1,p=u).B(b^1)>A(b^1,p=p).B(b^1)|1']\n"
     ]
    }
   ],
   "source": [
    "print([rr.as_string() for rr in rr1.generate([sp1])])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "`ReactionRule.generate` matters the order of `Species` in the given list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['A(b)+B(b)>A(b^1).B(b^1)|1']\n",
      "[]\n"
     ]
    }
   ],
   "source": [
    "rr1 = create_binding_reaction_rule(Species(\"A(b)\"), Species(\"B(b)\"), Species(\"A(b^1).B(b^1)\"), 1.0)\n",
    "sp1 = Species(\"A(b)\")\n",
    "sp2 = Species(\"B(b)\")\n",
    "print([rr.as_string() for rr in rr1.generate([sp1, sp2])])\n",
    "print([rr.as_string() for rr in rr1.generate([sp2, sp1])])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the current implementation, `ReactionRule.generate` does **not** always return a list of unique `ReactionRule`s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4\n",
      "['A(b,c^1).A(b,c^1)+B(b,c^1).B(b,c^1)>A(b^1,c^2).A(b,c^2).B(b^1,c^3).B(b,c^3)|1', 'A(b,c^1).A(b,c^1)+B(b,c^1).B(b,c^1)>A(b^1,c^2).A(b,c^2).B(b,c^3).B(b^1,c^3)|1', 'A(b,c^1).A(b,c^1)+B(b,c^1).B(b,c^1)>A(b,c^1).A(b^2,c^1).B(b^2,c^3).B(b,c^3)|1', 'A(b,c^1).A(b,c^1)+B(b,c^1).B(b,c^1)>A(b,c^1).A(b^2,c^1).B(b,c^3).B(b^2,c^3)|1']\n"
     ]
    }
   ],
   "source": [
    "sp1 = Species(\"A(b,c^1).A(b,c^1)\")\n",
    "sp2 = Species(\"B(b,c^1).B(b,c^1)\")\n",
    "print(rr1.count([sp1, sp2]))  # => 4\n",
    "print([rr.as_string() for rr in rr1.generate([sp1, sp2])])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`ReactionRules` listed above look different, but all the products suggest the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'A(b,c^1).A(b^2,c^1).B(b^2,c^3).B(b,c^3)'}\n"
     ]
    }
   ],
   "source": [
    "print(set([format_species(rr.products()[0]).serial() for rr in rr1.generate([sp1, sp2])]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is because these `ReactionRule`s are generated based on the diffent matches though they produces the same `Species`. For details, See the section below.\n",
    "\n",
    "Wildcard is also available in `ReactionRule`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['A(p=u^1).B(p^1)>A(p=p^1).B(p^1)|1']\n"
     ]
    }
   ],
   "source": [
    "rr1 = create_unimolecular_reaction_rule(Species(\"A(p=u^_)\"), Species(\"A(p=p^_)\"), 1.0)\n",
    "print([rr.as_string() for rr in rr1.generate([Species(\"A(p=u^1).B(p^1)\")])])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, a wildcard can work as a name of `UnitSpecies`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['A(p=u)>A(p=p)|1']\n",
      "['B(b^1,p=u).C(b^1,p=u)>B(b^1,p=p).C(b^1,p=u)|1', 'B(b^1,p=u).C(b^1,p=u)>B(b^1,p=u).C(b^1,p=p)|1']\n"
     ]
    }
   ],
   "source": [
    "rr1 = create_unimolecular_reaction_rule(Species(\"_(p=u)\"), Species(\"_(p=p)\"), 1.0)\n",
    "print([rr.as_string() for rr in rr1.generate([Species(\"A(p=u)\")])])\n",
    "print([rr.as_string() for rr in rr1.generate([Species(\"B(b^1,p=u).C(b^1,p=u)\")])])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Named wildcards in a state is useful to specify the correspondence between sites."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['AB(a=x, b=y)>B(b=y)+A(a=x)|1']\n",
      "['AB(a=y, b=x)>B(b=x)+A(a=y)|1']\n"
     ]
    }
   ],
   "source": [
    "rr1 = create_unbinding_reaction_rule(Species(\"AB(a=_1, b=_2)\"), Species(\"B(b=_2)\"), Species(\"A(a=_1)\"), 1.0)\n",
    "print([rr.as_string() for rr in rr1.generate([Species(\"AB(a=x, b=y)\")])])\n",
    "print([rr.as_string() for rr in rr1.generate([Species(\"AB(a=y, b=x)\")])])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nameless wildcard `_` does not care about equality between matches. Products are generated in order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_(b)+_(b)>_(b^1)._(b^1)|1\n",
      "['A(b)+A(b)>A(b^1).A(b^1)|1']\n",
      "['A(b)+B(b)>A(b^1).B(b^1)|1']\n"
     ]
    }
   ],
   "source": [
    "rr1 = create_binding_reaction_rule(Species(\"_(b)\"), Species(\"_(b)\"), Species(\"_(b^1)._(b^1)\"), 1.0)\n",
    "print(rr1.as_string())\n",
    "print([rr.as_string() for rr in rr1.generate([Species(\"A(b)\"), Species(\"A(b)\")])])\n",
    "print([rr.as_string() for rr in rr1.generate([Species(\"A(b)\"), Species(\"B(b)\")])])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For its symmetry, the former case above is sometimes preffered to have a half of the original kinetic rate. This is because the number of combinations of molecules in the former is $n(n-1)/2$ even though that in the later is $n^2$, where both numbers of A and B molecules are $n$. This is true for `gillespie` and `ode`. However, in `egfrd` and `spatiocyte`, a kinetic rate is intrinsic one, and not affected by the number of combinations. Thus, in E-Cell4, no modification in the rate is done even for the case. See [Homodimerization and Annihilation](../../Tests/Homodimerization_and_Annihilation.ipynb) for the difference between algorithms.\n",
    "\n",
    "In constrast to nameless wildcard, named wildcard keeps its consistency, and always suggests the same value in the `ReactionRule`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_1(b)+_1(b)>_1(b^1)._1(b^1)|1\n",
      "['A(b)+A(b)>A(b^1).A(b^1)|1']\n",
      "[]\n"
     ]
    }
   ],
   "source": [
    "rr1 = create_binding_reaction_rule(Species(\"_1(b)\"), Species(\"_1(b)\"), Species(\"_1(b^1)._1(b^1)\"), 1.0)\n",
    "print(rr1.as_string())\n",
    "print([rr.as_string() for rr in rr1.generate([Species(\"A(b)\"), Species(\"A(b)\")])])\n",
    "print([rr.as_string() for rr in rr1.generate([Species(\"A(b)\"), Species(\"B(b)\")])])  # => []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Named wildcard is consistent even between `UnitSpecies`' and `site`'s names, technically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A(b=_1)+_1(b)>A(b=_1^1)._1(b^1)|1\n",
      "[]\n",
      "['A(b=B)+B(b)>A(b=B^1).B(b^1)|1']\n"
     ]
    }
   ],
   "source": [
    "rr1 = create_binding_reaction_rule(Species(\"A(b=_1)\"), Species(\"_1(b)\"), Species(\"A(b=_1^1)._1(b^1)\"), 1.0)\n",
    "print(rr1.as_string())\n",
    "print([rr.as_string() for rr in rr1.generate([Species(\"A(b=B)\"), Species(\"A(b)\")])])  # => []\n",
    "print([rr.as_string() for rr in rr1.generate([Species(\"A(b=B)\"), Species(\"B(b)\")])])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.3. NetfreeModel\n",
    "\n",
    "`NetfreeModel` is a `Model` class for the rule-based modeling. The interface of `NetfreeModel` is almost same with `NetworkModel`, but takes into account rules and matches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "rr1 = create_binding_reaction_rule(Species(\"A(r)\"), Species(\"A(l)\"), Species(\"A(r^1).A(l^1)\"), 1.0)\n",
    "\n",
    "m1 = NetfreeModel()\n",
    "m1.add_reaction_rule(rr1)\n",
    "print(m1.num_reaction_rules())\n",
    "\n",
    "m2 = NetworkModel()\n",
    "m2.add_reaction_rule(rr1)\n",
    "print(m2.num_reaction_rules())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python notation explained in [2. How to Build a Model](tutorial2.ipynb) is available too. To get a model as `NetfreeModel`, set `is_netfree` `True` in `get_model`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<ecell4_base.core.NetfreeModel object at 0x145bb6e32998>\n"
     ]
    }
   ],
   "source": [
    "with reaction_rules():\n",
    "    A(r) + A(l) > A(r^1).A(l^1) | 1.0\n",
    "\n",
    "m1 = get_model(is_netfree=True)\n",
    "print(repr(m1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Model.query_reaction_rules` returns a list of `ReactionRule`s agaist the given reactants. `NetworkModel` just returns `ReactionRule`s based on the equality of `Species`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "print(len(m2.query_reaction_rules(Species(\"A(r)\"), Species(\"A(l)\"))))  # => 1\n",
    "print(len(m2.query_reaction_rules(Species(\"A(l,r)\"), Species(\"A(l,r)\"))))  # => 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the other hand, `NetfreeModel` genarates the list by applying the stored `ReactionRule`s in the rule-based way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "print(len(m1.query_reaction_rules(Species(\"A(r)\"), Species(\"A(l)\"))))  # => 1\n",
    "print(len(m1.query_reaction_rules(Species(\"A(l,r)\"), Species(\"A(l,r)\"))))  # => 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`NetfreeModel` does not cache generated objects. Thus, `NetfreeModel.query_reaction_rules` is slow, but needs less memory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A(l,r)+A(l,r)>A(l,r^1).A(l^1,r)|2\n",
      "A(l,r^1).A(l^1,r)+A(l,r)>A(l,r^1).A(l^1,r^2).A(l^2,r)|1\n",
      "A(l,r^1).A(l^1,r)+A(l,r^1).A(l^1,r)>A(l,r^1).A(l^1,r^2).A(l^2,r^3).A(l^3,r)|2\n"
     ]
    }
   ],
   "source": [
    "print(m1.query_reaction_rules(Species(\"A(l,r)\"), Species(\"A(l,r)\"))[0].as_string())\n",
    "print(m1.query_reaction_rules(Species(\"A(l,r^1).A(l^1,r)\"), Species(\"A(l,r)\"))[0].as_string())\n",
    "print(m1.query_reaction_rules(Species(\"A(l,r^1).A(l^1,r)\"), Species(\"A(l,r^1).A(l^1,r)\"))[0].as_string())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`NetfreeModel.expand` expands `NetfreeModel` to `NetworkModel` by iteratively applying `ReactionRule`s agaist the given set of `Species` (called seed)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n",
      "6\n",
      "['A(b)+A(b)>A(b^1).A(b^1)|1', 'A(b)+B(b)>A(b^1).B(b^1)|1', 'B(b)+B(b)>B(b^1).B(b^1)|1', 'A(b^1).A(b^1)>A(b)+A(b)|1', 'A(b^1).B(b^1)>A(b)+B(b)|1', 'B(b^1).B(b^1)>B(b)+B(b)|1']\n"
     ]
    }
   ],
   "source": [
    "with reaction_rules():\n",
    "    _(b) + _(b) == _(b^1)._(b^1) | (1.0, 1.0)\n",
    "\n",
    "m3 = get_model(True)\n",
    "print(m3.num_reaction_rules())\n",
    "\n",
    "m4 = m3.expand([Species(\"A(b)\"), Species(\"B(b)\")])\n",
    "print(m4.num_reaction_rules())\n",
    "print([rr.as_string() for rr in m4.reaction_rules()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To avoid the infinite iteration for expansion, you can limit the maximum number of iterations and of `UnitSpecies` in a `Species`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4\n",
      "['A(l,r)+A(l,r)>A(l,r^1).A(l^1,r)|2', 'A(l,r^1).A(l^1,r)+A(l,r^1).A(l^1,r)>A(l,r^1).A(l^1,r^2).A(l^2,r^3).A(l^3,r)|2', 'A(l,r)+A(l,r^1).A(l^1,r)>A(l,r^1).A(l^1,r^2).A(l^2,r)|2', 'A(l,r)+A(l,r^1).A(l^1,r^2).A(l^2,r)>A(l,r^1).A(l^1,r^2).A(l^2,r^3).A(l^3,r)|2']\n"
     ]
    }
   ],
   "source": [
    "m2 = m1.expand([Species(\"A(l, r)\")], 100, {Species(\"A\"): 4})\n",
    "print(m2.num_reaction_rules())  # => 4\n",
    "print([rr.as_string() for rr in m2.reaction_rules()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.4. Differences between Species, ReactionRule and NetfreeModel\n",
    "\n",
    "The functions explained above is similar, but there are some differences in the criteria."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "['A(b^1).A(b^1)>A(b)+A(b)|1']\n"
     ]
    }
   ],
   "source": [
    "sp1 = Species(\"A(b^1).A(b^1)\")\n",
    "sp2 = Species(\"A(b)\")\n",
    "rr1 = create_unbinding_reaction_rule(sp1, sp2, sp2, 1.0)\n",
    "print(count_species_matches(sp1, sp1))\n",
    "print([rr.as_string() for rr in rr1.generate([sp1])])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Though `count_species_matches` suggests two different ways for matching (left-right and right-left), `ReactionRule.generate` returns only one `ReactionRule` because the order doesn't affect the product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n",
      "['A(b^1,c^2).A(b^3,c^2).B(b^1).B(b^3)>A(b,c^1).A(b^2,c^1).B(b^2)+B(b)|1', 'A(b^1,c^2).A(b^3,c^2).B(b^1).B(b^3)>A(b^1,c^2).A(b,c^2).B(b^1)+B(b)|1']\n"
     ]
    }
   ],
   "source": [
    "sp1 = Species(\"A(b^1).B(b^1)\")\n",
    "rr1 = create_unbinding_reaction_rule(\n",
    "    sp1, Species(\"A(b)\"), Species(\"B(b)\"), 1.0)\n",
    "sp2 = Species(\"A(b^1,c^2).A(b^3,c^2).B(b^1).B(b^3)\")\n",
    "print(count_species_matches(sp1, sp2))\n",
    "print([rr.as_string() for rr in rr1.generate([sp2])])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, `ReactionRule.generate` works similarly with `count_species_matches`. However, `NetfreeModel.query_reaction_rules` returns only one `ReationRule` with higher kinetic rate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['A(b^1,c^2).B(b^1).A(b^3,c^2).B(b^3)>A(b,c^1).A(b^2,c^1).B(b^2)+B(b)|2']\n"
     ]
    }
   ],
   "source": [
    "m1 = NetfreeModel()\n",
    "m1.add_reaction_rule(rr1)\n",
    "print([rr.as_string() for rr in m1.query_reaction_rules(sp2)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`NetfreeModel.query_reaction_rules` checks if each `ReactionRule` generated is the same with others, and summalizes it if possible.\n",
    "\n",
    "As explaned above, `ReactionRule.generate` matters the order of `Species`, but `Netfree.query_reaction_rules` does not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['A(b)+B(b)>A(b^1).B(b^1)|1']\n",
      "['A(b)+B(b)>A(b^1).B(b^1)|1']\n",
      "[]\n",
      "['A(b)+B(b)>A(b^1).B(b^1)|1']\n"
     ]
    }
   ],
   "source": [
    "sp1 = Species(\"A(b)\")\n",
    "sp2 = Species(\"B(b)\")\n",
    "rr1 = create_binding_reaction_rule(sp1, sp2, Species(\"A(b^1).B(b^1)\"), 1.0)\n",
    "m1 = NetfreeModel()\n",
    "m1.add_reaction_rule(rr1)\n",
    "print([rr.as_string() for rr in rr1.generate([sp1, sp2])])\n",
    "print([rr.as_string() for rr in m1.query_reaction_rules(sp1, sp2)])\n",
    "print([rr.as_string() for rr in rr1.generate([sp2, sp1])])  # => []\n",
    "print([rr.as_string() for rr in m1.query_reaction_rules(sp2, sp1)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Named wildcards must be consistent in the context while nameless wildcards are not necessarily consistent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "1\n",
      "1\n",
      "1\n",
      "['A(b)+A(b)>A(b^1).A(b^1)|1']\n",
      "['A(b)+B(b)>A(b^1).B(b^1)|1']\n",
      "['A(b)+A(b)>A(b^1).A(b^1)|1']\n",
      "[]\n"
     ]
    }
   ],
   "source": [
    "sp1 = Species(\"_(b)\")\n",
    "sp2 = Species(\"_1(b)\")\n",
    "sp3 = Species(\"A(b)\")\n",
    "sp4 = Species(\"B(b)\")\n",
    "rr1 = create_binding_reaction_rule(sp1, sp1, Species(\"_(b^1)._(b^1)\"), 1)\n",
    "rr2 = create_binding_reaction_rule(sp2, sp2, Species(\"_1(b^1)._1(b^1)\"), 1)\n",
    "print(count_species_matches(sp1, sp2))  # => 1\n",
    "print(count_species_matches(sp1, sp3))  # => 1\n",
    "print(count_species_matches(sp2, sp2))  # => 1\n",
    "print(count_species_matches(sp2, sp3))  # => 1\n",
    "print([rr.as_string() for rr in rr1.generate([sp3, sp3])])\n",
    "print([rr.as_string() for rr in rr1.generate([sp3, sp4])])\n",
    "print([rr.as_string() for rr in rr2.generate([sp3, sp3])])\n",
    "print([rr.as_string() for rr in rr2.generate([sp3, sp4])])  # => []"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}