6

New to python coming from MATLAB.

I am using a hyperbolic tangent truncation of a magnitude-scale function. I encounter my problem when applying the 0.5 * math.tanh(r/rE-r0) + 0.5 function onto an array of range values r = np.arange(0.1,100.01,0.01). I get several 0.0 values for the function on the side approaching zero, which cause domain issues when I perform the logarithm:

P1 = [ (0.5*m.tanh(x / rE + r0 ) + 0.5) for x in r] # truncation function

I use this work-around:

P1 = [ -m.log10(x) if x!=0.0 else np.inf for x in P1 ]

which is sufficient for what I am doing but is a bit of a band-aid solution.

As requested for mathematical explicitness:

In astronomy, the magnitude scale works roughly as such:

mu = -2.5log(flux) + mzp # apparent magnitude

where mzp is the magnitude at which one would see 1 photon per second. Therefore, greater fluxes equate to smaller (or more negative) apparent magnitude. I am making models for sources which use multiple component functions. Ex. two sersic functions with different sersic indices with a P1 outer truncation on the inner component and a 1-P1 inner truncation on the outer component. This way, when adding the truncation function to each component, the magnitude as defined by radius, will become very large because of how small mu1-2.5*log(P1) gets as P1 asymptotically approaches zero.

TLDR: What I would like to know is if there is a way of preserving floating points whose accuracy is insufficient to be distinguishable from zero (in particular in the results of functions that asymptotically approach zero). This important because when taking the logarithm of such numbers a domain error is the result.

The last number before the output in the non-logarithmic P1 starts reading zero is 5.551115123125783e-17, which is a common floating point arithmetic rounding error result where the desired result should be zero.

Any input would be greatly appreciated.

@user:Dan without putting my whole script:

xc1,yc1 = 103.5150,102.5461;
Ee1 = 23.6781;
re1 = 10.0728*0.187;
n1 = 4.0234;
# radial brightness profile (magnitudes -- really surface brightness but fine in ex.)
mu1 = [ Ee1 + 2.5/m.log(10)*bn(n1)*((x/re1)**(1.0/n1) - 1) for x in r];

# outer truncation
rb1 = 8.0121
drs1 = 11.4792

P1 = [ (0.5*m.tanh( (2.0 - B(rb1,drs1) ) * x / rb1 + B(rb1,drs1) ) + 0.5) for x in r]

P1 = [ -2.5*m.log10(x) if x!=0.0 else np.inf for x in P1 ] # band-aid for problem

mu1t = [x+y for x,y in zip(P1,mu1)] # m1 truncated by P1

where bn(n1)=7.72 and B(rb1,drs1) = 2.65 - 4.98 * ( r_b1 / (-drs1) );

mu1 is the magnitude profile of the component to be truncated. P1 is the truncation function. Many of the final entries for P1 are zero, which is due to the floating points being undistinguished from zero due to the floating point accuracy.

An easy way to see the problem:

>>> r = np.arange(0,101,1)
>>> P1 = [0.5*m.tanh(-x)+0.5 for x in r]
>>> P1
[0.5, 0.11920292202211757, 0.01798620996209155, 0.002472623156634768, 0.000335350130466483, 4.539786870244589e-05, 6.144174602207286e-06, 8.315280276560699e-07, 1.1253516207787584e-07, 1.5229979499764568e-08, 2.0611536366565986e-09, 2.789468100949932e-10, 3.775135759553905e-11, 5.109079825871277e-12, 6.914468997365475e-13, 9.35918009759007e-14, 1.2656542480726785e-14, 1.7208456881689926e-15, 2.220446049250313e-16, 5.551115123125783e-17, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

Note also the floats before zeros.

Connor
  • 397
  • 2
  • 10
  • @GreenAsJade's suggestion won't fix *this* problem though. It might even compound it. NumPy is not a computer algebra system and loves fixed-precision numbers. –  Jan 17 '14 at 22:52
  • @GreenAsJade: He appears to be using numpy already, with `np.arange()`. – Dan Jan 17 '14 at 22:57
  • Yes, I have looked at numpy. Are you suggesting that I use numpy as an alternative to math.tanh? ex. np.tanh? I have. It gives the same result. – Connor Jan 17 '14 at 23:00
  • 1
    What exactly are you trying to accomplish? A cleaner way of dealing with cases where there isn't enough floating point accuracy to distinguish something from 0.0? Also, what exactly do you mean by "faint" numbers and "magnitude scale"? Could you be more mathematically explicit about what exactly you're doing with these numbers after you calculate them? There may be a more numerically stable mathematical workaround that we can't see because we don't know what these numbers mean or what they're being used for. – Dan Jan 17 '14 at 23:05
  • Precisely. I am looking for a way to distinguish numbers whose floating point accuracy is insufficient to distinguish from zero. I will edit the OP for clarity. Thanks. – Connor Jan 17 '14 at 23:10
  • Maybe it would be easier if you gave the MATLAB version of what you wanted as an example? – DSM Jan 17 '14 at 23:12
  • What are typical values for `rE` and `r0`? It's difficult to figure out a more numerically stable way to do this without knowing what's big and what's small. – Dan Jan 17 '14 at 23:14
  • Also, could you explain what you're using `P1` for, and what precisely the magnitude scale you mentioned is mathematically? It's hard to make a workaround without this additional context. – Dan Jan 17 '14 at 23:16
  • @Dan, Yes. P1, as listed before logarithmic application, would be used to damp the flux of a certain component of the model of a light source. At `rb1`, the P1 is 0.99 and at `rb1+drs1` P1 is 0.01 . It is an asymptotic function between 1 and 0. It follows that to apply this function to apply it to the component (which is distributed by radius) -2.5*log(P1(r)) is added at each corresponding radius interval. – Connor Jan 17 '14 at 23:48
  • Just an FYI, but in the future you might want to ask questions like this on the [computational science stack exchange](http://scicomp.stackexchange.com). It's the same company as StackOverflow, but they have MathJax enabled so it's easier to write about math. – Dan Jan 18 '14 at 02:01

2 Answers2

5

Recall that the hyperbolic tangent can be expressed as (1-e^{-2x})/(1+e^{-2x}). With a bit of algebra, we can get that 0.5*tanh(x)-0.5 (the negative of your function) is equal to e^{-2x}/(1+e^{-2x}). The logarithm of this would be -2*x-log(1+exp(-2*x)), which would work and be stable everywhere.

That is, I recommend you replace:

P1 = [ (0.5*m.tanh( (2.0 - B(rb1,drs1) ) * x / rb1 + B(rb1,drs1) ) + 0.5) for x in r]

P1 = [ -2.5*m.log10(x) if x!=0.0 else np.inf for x in P1 ] # band-aid for problem

With this simpler and more stable way of doing it:

r = np.arange(0.1,100.01,0.01)
#r and xvals are numpy arrays, so numpy functions can be applied in one step
xvals=(2.0 - B(rb1,drs1) ) * r / rb1 + B(rb1,drs1)
P1=2*xvals+np.log1p(np.exp(-2*xvals))
Dan
  • 12,157
  • 12
  • 50
  • 84
  • Sorry, I might have been unclear with all the edits now. You may have noticed that the function I am looking at better resembles 0.5*m.tanh(-x)+0.5 (as I gave in my last edit). The problem is that as the function tends to zero as x becomes large, the numbers start actually reading zero. This is no good when I want to take the logarithm to add it to something in magnitude scales. After ~1e-17, there are just zero-floats. How can I get non-zero numbers here (however small in significant figures, I really just need the exponent to be defined for the logarithm). – Connor Jan 18 '14 at 00:22
  • 1
    @user3208266: Ah! Sorry, I wrote most of this before your last edit. Give me a bit and I'll fix it for asymptotically large instead of asymptotically small values. – Dan Jan 18 '14 at 00:26
  • Before you do, what do you mean. The large values are not the concern since the function is defined confined between y=1 and y=0 asymptotically. The problem is with handling very small floats as x gets big and y gets closer and closer to zero, because python starts not to distinguish between these floats and zero, and so just outputs zero. – Connor Jan 18 '14 at 00:31
  • 1
    I mean large values of `x` in your function above. – Dan Jan 18 '14 at 00:37
  • This post is great too. You can bet you're bottom dollar I'm going to use this someday. – Connor Jan 18 '14 at 00:48
  • 1
    @user3208266: Ok, it's designed for large x values now. – Dan Jan 18 '14 at 00:55
  • If I had the will to do the math to check that this is correct, I'd definitely +1. Ah heck I'll do it anyway. – Cody Piersall Jan 18 '14 at 01:06
  • Excellent. This sounds like a very reasonable solution. I'll check the math and let let you know how it works on Monday with an update on the algebra. ps. don't know who down-voted this or why.. +1 – Connor Jan 18 '14 at 01:07
  • 2
    I'd suggest using `log1p(...)` instead of `log(1 + ...)`. I don't think it makes that much difference in this particular case, but in general it's safer when the `...` is small. – Mark Dickinson Jan 18 '14 at 19:21
  • @Mark Dickinson: I edited my answer to do that instead. Thanks for the tip! – Dan Jan 19 '14 at 00:46
  • @user3208266: If my post is wrong in some way, would you mind editing it to correct it? – Dan Jan 20 '14 at 21:39
  • I gave it a touch-up. Turns out there was just a negative sign missing in one spot. – Connor Jan 20 '14 at 23:09
1

Two things you can try.

(1) brute force approach: find a variable-precision float arithmetic package and use that instead of built-in fixed precision. I am playing with your problem in Maxima [1] and I find that I have to increase the float precision quite a lot in order to avoid underflow, but it is possible. I can post the Maxima code if you want. I would imagine that there is some suitable variable-precision float library for Python.

(2) approximate log((1/2)(1 + tanh(-x)) with a Taylor series or some other kind of approximation in order to avoid the log(tanh(...)) altogether.

[1] http://maxima.sourceforge.net

Robert Dodier
  • 16,905
  • 2
  • 31
  • 48
  • I had never heard of Maxima and thanks for the link. I will continue to look more rigorously for a variable precision float library in Python. – Connor Jan 18 '14 at 00:59
  • I'd recommend gmpy2 for Python: it wraps GNU's MPFR library for variable precision (binary) floating-point. – Mark Dickinson Jan 18 '14 at 19:15