title: integrate.odeint pausing (perhaps near removable discontinuity?)
body: <p>I have a system of equations that I am trying to integrate with <code>integrate.odeint</code>.  </p>
  
  <p>Within that system are terms like <code>XiXj[i,j]*XiXj[j,k]/X[j]</code>.  It is possible for the denominator to be zero (particularly in the initial condition).  However, both terms in the numerator should go to zero there, and the whole term should be interpreted as zero (and values going into the equations are nonnegative).</p>
  
  <p>What I have done to handle this is to replace  the <code>1/X[j]</code> with <code>Xinv[j]</code> which is set to be <code>1/X[j]</code> when <code>X[j]&gt;eps</code> for some smallish <code>eps</code>, and use <code>1</code> if it is smaller than <code>eps</code>.</p>
  
  <p>However, when I integrate, I find that <code>odeint</code> will sit for a very long time calculating derivatives at what appears to be a single value of <code>t</code>.  I think this is because of it trying to calculate things near a discontinuity in what I've done.</p>
  
  <p>Larger values of <code>eps</code> or assigning a minimum step size reduces the probability of encountering this issue, but it can still happen.  </p>
  
  <p>I'm hoping to figure out how to prevent this from happening.  Unfortunately, I haven't been able to force this to happen in a simple example.  So I've stripped down my example as far as I can.</p>
  
  <p>Some background on what's being solved: I'm thinking about nodes in a networkx network (so running this code will require networkx).  Nodes are susceptible, infected or recovered.  I have equations for the derivatives of the susceptible (X) and infected (Y) states of each node.  As correlations form across edges, I am also tracking the probability each edge is between S and I (XY) or S and S (XX) nodes.</p>
  
  <p>Here's a system of equations.  I've set g_{ij} to be 1 and gamma_i is simply gamma in the code below.
  <a href="https://i.stack.imgur.com/YUbBG.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/YUbBG.png" alt="enter image description here"></a></p>
  
  <p>It should be possible to just copy and paste the code and it will run it.  As it runs, the derivative function outputs the time and some simple measure of the derivatives.  You'll notice that it pauses for a long time at particular times.  Sometimes it starts up again.  Sometimes I give up waiting.</p>
  
  <p>(note, the network it solves for is random, so if by some unlikely chance you don't see the behavior, run it again.  You may want to pickle the graph and reuse a graph to make things consistent between tests.)</p>
  
  <pre><code>from scipy import integrate
  import scipy
  import networkx as nx
  
  
  def _dSIR_pair_based_(V, t, G, nodelist, index_of_node, tau, gamma):
      #V is a vector of states
      #t is time
      #G is a networkx Graph
      #nodelist is a list of nodes in G
      #index_of_node is a dict index_of_node[u] is i such that nodelist[i]=u
      #    {u:i for i, u in enumerate(nodelist)}
      #tau is transmission rate
      #gamma is recovery rate
  
      N=G.order()
      X = V[0:N] #probability node is susceptible
      eps = 0.1**(5)
      Xinv = scipy.array([1/v if v&gt;eps else 1 for v in X]) #there are places where we divide by X[i] which may == 0.
                                                           #In those cases the numerator is (very) 0, so it's easier
                                                           #to set this up as mult by inverse with a dummy value when
                                                           #it is 1/0.
      Y = V[N:2*N] #prob node is infected
      Yinv = scipy.array([1/v if v&gt;eps else 1 for v in Y])
  
      XY = V[2*N: 2*N+N**2]  #prob edge exists and is from susceptible to infected 
      XX = V[2*N+N**2:]      #prob edge exists and is between two susceptibles.
  
      #print X.shape, Y.shape, XY.shape, XX.shape, N
      XY.shape = (N,N)   #get them into the right shape
      XX.shape = (N,N)
  
  
      YX = XY.T #not really needed, but helps keep consistent with equations as written in text
  
      dX = scipy.zeros(N)
      dY = scipy.zeros(N)
      dXY = scipy.zeros((N,N))
      dXX = scipy.zeros((N,N))
  
  
      #I could make the below more efficient, but I think this sequence of for loops is easier to read, or at least understand.
      #I expect it to run quickly regardless.  Will avoid (premature) optimization for now.
      for u in nodelist:
          i = index_of_node[u]
          dY[i] += -gamma*Y[i] 
          for v in G.neighbors(u):
              j = index_of_node[v]
              dX[i] += -tau*XY[i,j]
              dY[i] += tau*XY[i,j]
  
              dXY[i,j] +=  - (tau+gamma)*XY[i,j] 
  
              for w in G.neighbors(u):
                  if w == v: #skip these
                      continue
                  #so w != v.  
                  k= index_of_node[w]
                  dXY[i,j] += tau * XX[i,j] * XY[j,k]*Xinv[j]  -  tau * YX[k,i] * XY[i,j]*Xinv[i]
                  dXX[i,j] += -tau * XX[i,j] * XY[j,k]*Xinv[j] -  tau * YX[k,i] * XX[i,j]*Xinv[i]  
  
      dXY.shape = (N**2,1)
      dXX.shape = (N**2,1)
  
      dV = scipy.concatenate((dX[:, None], dY[:,None], dXY, dXX), axis=0).T[0]
      print t, sum(dV)
      return dV
  
  
  
  
  
  def SIR_pair_based(G, nodelist, Y0, tau, gamma, tmin = 0, tmax = 100, tcount = 1001):
  
      times = scipy.linspace(tmin,tmax,tcount)
      X0 = 1-Y0 #if not initially infected, then initially susceptible
      N = len(Y0)
  
      XY0 = X0[:,None]*Y0[None,:]
      XX0 = X0[:,None]*X0[None,:]
  
      XY0.shape=(N**2,1)
      XX0.shape=(N**2,1)
  
      V0 = scipy.concatenate((X0[:,None], Y0[:,None], XY0, XX0), axis=0).T[0]
      index_of_node = {node:i for i, node in enumerate(nodelist)}
  
      V = integrate.odeint(_dSIR_pair_based_, V0, times, args = (G, nodelist, index_of_node, tau, gamma))#, mxstep=10)#(times[1]-times[0])/1000)
  
  
  
  
  def test_SIR_pair_based():
      G = nx.fast_gnp_random_graph(100,0.02)
      nodelist = G.nodes()
      Y0 = scipy.array([1 if node&lt;10 else 0 for node in nodelist])
  
      SIR_pair_based(G, nodelist, Y0, 1, 0.5, tmax = 10)
  
  
  test_SIR_pair_based()
  </code></pre>
  
  <p>If I plot <code>sum(X)</code>, <code>sum(Y)</code> and <code>1-sum(X+Y)</code> I get the following for a larger value of <code>eps=0.01</code> (when the random graph is such that it does get through the process - sometimes it still gets stuck):</p>
  
  <p><a href="https://i.stack.imgur.com/i3Odm.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/i3Odm.png" alt="enter image description here"></a></p>