Python: Generating a list and then matching question

[MathewsMD]

This seems like a fairly common technique but I'm fairly new to programming and don't quite know what the proper terms to find such an algorithm would be. In my case, it seems like a list of all possible values would approach almost a billion, and I feel like there's a quicker approach to generate the list and match it to a desired value. Essentially, I have a 2D function that I am integrating and I have constants in the function that I want to find all possible values for. Here's the basic integration code:

import numpy as np

from scipy import integrate

i = 0

while i < 1: #run once

i = i + 1

sigma0 = 4

r0 = 0

theta0 = 0

def G(r,theta):

    return (np.e**(-((r**2 + r0**2 - 2*r*r0*(np.cos(theta)*np.cos(theta0) + \
    np.sin(theta)*np.sin(theta0)))/(2*sigma0**2))))*r

R1I = integrate.nquad(G, [[0,4],[0,2*np.pi]])
R2I = integrate.nquad(G, [[4,7],[0,0.5*np.pi]])
R3I = integrate.nquad(G, [[4,7],[0.5*np.pi,np.pi]])
R4I = integrate.nquad(G, [[4,7],[np.pi,1.5*np.pi]])
R5I = integrate.nquad(G, [[4,7],[1.5*np.pi,2*np.pi]])

In the case above, I specified values for sigma0, r0, and theta0, but I don't want to keep these constant (i.e. I don't want to specify them). What I want to do is set constraints and then find all possible values these can 3 variables can take on. Essentially, I want to specify a range for R1I, R2I, R3I, R4I, and R5I, and then find all possible combination of sigma0, r0, and theta0 in the function that would allow this. Making the range as small as possible and on arbitrary values (e.g. 0.055 < R1I < 0.060, 0.008 < R2I < 0.013, 0.016 < R3I < 0.025, 0.006 < R4I < 0.007, 0.020 < R5I < 0.988) would be optimal.

Any advice or recommendations would be greatly appreciated! Thank you!

 

[Carno Raar]

This seems like a fairly common technique but I'm fairly new to programming and don't quite know what the proper terms to find such an algorithm would be.

I'm not a mathmatican and there are many who could give you a better reply than I could. However in the absence of other replies, I'll give it a go.

The field is "numerical analysis" (http://en.wikipedia.org/wiki/Numerical_analysis) and you're probably looking for a way of solving it by successively more accurate approximations using methods such as Newton Rapheson (http://en.wikipedia.org/wiki/Newton's_method).

Take this with a pinch of salt as this is not something I am proficient in.

 

[MathewsMD]

I'm not a mathmatican and there are many who could give you a better reply than I could. However in the absence of other replies, I'll give it a go.

The field is "numerical analysis" (http://en.wikipedia.org/wiki/Numerical_analysis) and you're probably looking for a way of solving it by successively more accurate approximations using methods such as Newton Rapheson (http://en.wikipedia.org/wiki/Newton's_method).

Take this with a pinch of salt as this is not something I am proficient in.

Thank you for the links! Yes, that seems exactly like what I want to do. I guess implementing the exact code to do this is where I'm stuck and trying to work at right now.

 

[wle]

In the case above, I specified values for sigma0, r0, and theta0, but I don't want to keep these constant (i.e. I don't want to specify them).

This isn't a solution to your problem in itself (and I'm not sure if it's still useful to post a reply), but one thing Python let's you do is define nested functions (i.e., functions defined inside other functions) and pass/return functions around like any other variable, which might help you structure your code a bit. First, personally I'd rewrite your function to use intermediate variables for readability:

def g(r, theta):
  num = r**2 + r0**2 - 2 * r * r0 * np.cos(theta - theta0)
  den = 2 * sigma0**2
  return r * np.exp(-num/den)

Then you can wrap the whole definition of g inside another function that takes r0, theta0, and sigma0 as parameters, e.g.:

def make_g(r0, theta0, sigma0):
  def g(r, theta):
    num = r**2 + r0**2 - 2 * r * r0 * np.cos(theta - theta0)
    den = 2 * sigma0**2
    return r * np.exp(-num/den)
  return g

If you do this, then you can type e.g. g = make_g(0, 0, 4) to get a version of the function with r0, theta0, and sigma0 set to 0, 0, and 4.*

I'm not sure what to suggest regarding the five integrals since I'm not sure what you want to do with them, but, for example, you can write a function that takes r0, theta0, and sigma0 as parameters and returns a tuple or list containing the five results if that's what you want:

def f(r0, theta0, sigma0):
  g = make_g(r0, theta0, sigma0)
  r1i = integrate.nquad(g, [[0, 4], [0, 2*np.pi]])
  r2i = integrate.nquad(g, [[4, 7], [0, 0.5*np.pi]])
  r3i = integrate.nquad(g, [[4, 7], [0.5*np.pi, np.pi]])
  r4i = integrate.nquad(g, [[4, 7], [np.pi, 1.5*np.pi]])
  r5i = integrate.nquad(g, [[4, 7], [1.5*np.pi, 2*np.pi]])
  return r1i, r2i, r3i, r4i, r5i

(this is assuming you don't also want to make the constants 0, 4, 2*np.pi etc. variables). If you load or type the contents of the last two code boxes into the Python interactive shell, you can try this function out to see what it does:

>>> f(0, 0, 4)
((39.555852443507625, 4.3915818141457424e-13), (9.808441643120279, 1.0889557747884337e-13), (9.808441643120279, 1.0889557747884337e-13), (9.808441643120279, 1.0889557747884337e-13), (9.808441643120279, 1.0889557747884337e-13))

This returns five pairs of numbers because integrate.nquad returns a pair of numbers (I'm not very familiar with SciPy, but these are presumably an estimate of both the integral and an error margin). You can modify the code above accordingly if, e.g., you just want the integrals.*If you want a term to look up, this sort of object is called a closure.

What I want to do is set constraints and then find all possible values these can 3 variables can take on. Essentially, I want to specify a range for R1I, R2I, R3I, R4I, and R5I, and then find all possible combination of sigma0, r0, and theta0 in the function that would allow this. Making the range as small as possible and on arbitrary values (e.g. 0.055 < R1I < 0.060, 0.008 < R2I < 0.013, 0.016 < R3I < 0.025, 0.006 < R4I < 0.007, 0.020 < R5I < 0.988) would be optimal.

Well I doubt you can find all of them because there is probably an infinite range of them. One thing you can try to do is pick a set of values for two of the variables (e.g. r0 and theta0) and try to find which values (if any) of the third variable (e.g., sigma0) satisfy the constraints.

It's probably a good idea to try to manually guess some solutions (try different values of r0, theta0, and sigma0 yourself) and/or plot some of the integral values as functions of some of your variables to get a feel for what the solutions look like. When you eventually try to determine a whole set of solutions, this could both help you guess good initial values to use as well as give you some assurance that the results you're getting aren't complete nonsense.

If by "making the range as small as possible" you really mean "approximately zero", SciPy already has some root-finding functions that you should be able to use.

Potentially useful things to look up (if you don't know them already):

• SciPy root-finding functions: http://docs.scipy.org/doc/scipy/reference/optimize.html#root-finding.
• matplotlib: http://matplotlib.org/.
• Python language features: tuples, using tuples to return multiple values from a function, lists (may work better for you than tuples), anonymous (lambda) functions, and closures.