[SciPy-User] Scaling up root-finding

[jkhilmer at chemistry.montana.edu]

Chris,

Just as an academic exercise, don't forget:

1. This is a highly parallel problem.
2. You need to dramatically modify the parameters of the root-finding
function: by default they will iterate far too many times.  You probably
only have ~10 bits of precision, correct?
3. Don't store 512x512 polynomials: break it down to maybe 10-20.  Remember
that excessive precision isn't necessary.

I'd be willing to bet that #2 above makes a huge improvement.

Jonathan


On Wed, Sep 4, 2013 at 3:27 PM, Chris Weisiger <cweisiger at msg.ucsf.edu>wrote:

> That...makes a lot of sense, actually. Why bother manually inverting the
> function when I can just switch the order of parameters I pass to
> numpy.polyfit? And then (when correcting actual data) just plug my measured
> values directly into the resulting polynomials?
>
> Intuitively this sounds great. In practice my brain's a bit fried right
> now so I'm going to put off really working on it until tomorrow. We'll see
> how it goes.
>
> Thanks for the advice!
>
> -Chris
>
>
> On Wed, Sep 4, 2013 at 1:56 PM, Charles R Harris <
> charlesr.harris at gmail.com> wrote:
>
>>
>>
>>
>> On Wed, Sep 4, 2013 at 2:35 PM, Chris Weisiger <cweisiger at msg.ucsf.edu>wrote:
>>
>>> We have a somewhat elderly microscopy camera whose response curve is
>>> best-characterized as a degree-2 polynomial (starts out steep, then
>>> flattens out; we hit the camera's maximum bit depth before the function
>>> stops being monotonic though). That is:
>>>
>>> reported photons = a + b * (incident photons) + c * (incident photons)^2
>>>
>>> Of course, each pixel is slightly different, giving 512*512 different
>>> polynomials. I want to linearize the camera's response in post-production,
>>> by inverting the (monotonic) function so that it maps reported photons to
>>> incident photons.
>>>
>>> The obvious method would be to use something from scipy.optimize's
>>> root-finding functions. I'm not remotely a numerical analysis expert, so I
>>> have no feel for which function would be best for this problem.
>>> Performance-wise, a quick test indicates that it would take about 3.5s on
>>> my laptop to process a single frame by manually iterating over every pixel
>>> in the frame and calling scipy.optimize.newton(that pixel's response
>>> function, that pixel's starting value). This isn't great, but it's well
>>> within the realm of feasibility.
>>>
>>> Is there some more efficient method to do this? In the ideal world we'd
>>> be able to handle arbitrary-degree polynomials, but I'd be willing to
>>> accept being restricted to degree-2. Of course, degree-1 is trivial to
>>> implement.
>>>
>>> Thanks for any advice you care to share.
>>>
>>>
>> Depends on the data you use to fit the curve. I always fit the inverse
>> response rather than the response.
>>
>> Chuck
>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20130904/c307fda8/attachment.html>