Finding the self-consistent solution to an equation
Question: At the bottom of this question are a set of functions transcribed from a
published neural-network model. When I call `R`, I get the following error:

> RuntimeError: maximum recursion depth exceeded while calling a Python object

Note that within each call to `R`, a recursive call to `R` is made for every
other neuron in the network. This is what causes the recursion depth to be
exceeded. Each return value for `R` depends on all the others (with the
network involving `N = 512` total values.) Does anyone have any idea what
method should be used to compute the self-consistent solution for `R`? Note
that `R` itself is a smooth function. I've tried treating this as a vector
root-solving problem -- but in this case the 512 dimensions are not
independent. With so many degrees of freedom, the roots are never found (using
the `scipy.optimize` functions). Does Python have any tools that can help with
this? Maybe it would be more natural to solve `R` using something like
Mathematica? I don't know how this is normally done.

    
    
    """Recurrent model with strong excitatory recurrence."""
    
    
    import numpy as np
    
    
    l = 3.14
    
    
    def R(x_i):
    
        """Steady-state firing rate of neuron at location x_i.
    
        Parameters
        ----------
        x_i : number
          Location of this neuron.
    
        Returns
        -------
        rate : float
          Firing rate.
    
        """
    
        N = 512
        T = 1
    
        x = np.linspace(-2, 2, N)
        sum_term = 0
        for x_j in x:
            sum_term += J(x_i - x_j) * R(x_j)
    
        rate = I_S(x_i) + I_A(x_i) + 1.0 / N * sum_term - T
    
        if rate < 0:
            return 0
    
        return rate
    
    
    def I_S(x):
        """Sensory input.
    
        Parameters
        ----------
        x : number
          Location of this neuron.
    
        Returns
        -------
        float
          Sensory input to neuron at x.
    
        """
        S_0 = 0.46
        S_1 = 0.66
        x_S = 0
        sigma_S = 1.31
        return S_0 + S_1 * np.exp(-0.5 * (x - x_S) ** 2 / sigma_S ** 2)
    
    
    def I_A(x):
        """Attentional additive bias.
    
        Parameters
        ----------
        x : number
          Location of this neuron.
    
        Returns
        -------
        number
          Additive bias for neuron at x.
    
        """
        x_A = 0
        A_1 = 0.089
        sigma_A = 0.35
        A_0 = 0
        sigma_A_prime = 0.87
        if np.abs(x - x_A) < l:
            return (A_1 * np.exp(-0.5 * (x - x_A) ** 2 / sigma_A ** 2) +
                    A_0 * np.exp(-0.5 * (x - x_A) ** 2 / sigma_A_prime ** 2))
        return 0
    
    
    def J(dx):
        """Connection strength.
    
        Parameters
        ----------
        dx : number
          Neuron i's distance from neuron j.
    
        Returns
        -------
        number
          Connection strength.
    
        """
        J_0 = -2.5
        J_1 = 8.5
        sigma_J = 1.31
        if np.abs(dx) < l:
            return J_0 + J_1 * np.exp(-0.5 * dx ** 2 / sigma_J ** 2)
        return 0
    
    
    if __name__ == '__main__':
    
        pass
    


Answer: This recursion never ends since there is no termination condition before
recursive call, adjusting maximum recursion depth does not help

    
    
    def R(x_i): 
       ...
       for x_j in x:
           sum_term += J(x_i - x_j) * R(x_j)
    

Perhaps you should be doing something like

    
    
    # some suitable initial guess
    state = guess
    
    while True: # or a fixed number of iterations
       next_state = compute_next_state(state)
    
       if some_condition_check(state, next_state):
           # return answer
           return state 
    
       if some_other_check(state, next_state):
           # something wrong, terminate
          raise ...
    

