|
| 1 | +{ |
| 2 | + "cells": [ |
| 3 | + { |
| 4 | + "cell_type": "markdown", |
| 5 | + "metadata": {}, |
| 6 | + "source": [ |
| 7 | + "# On dynamic range: the basics\n", |
| 8 | + "## Intro to (single) floating point precision\n", |
| 9 | + "\n", |
| 10 | + "### Background: \n", |
| 11 | + "See https://en.wikibooks.org/wiki/Floating_Point/Floating-Point_Numbers for an introduction to floating point precision, range, significands and exponents.\n", |
| 12 | + "In particular, https://en.wikibooks.org/wiki/Floating_Point/Floating_Point_Formats is a useful reference page on the precision of single vs. double floats (IEEE754 standard)." |
| 13 | + ] |
| 14 | + }, |
| 15 | + { |
| 16 | + "cell_type": "markdown", |
| 17 | + "metadata": {}, |
| 18 | + "source": [ |
| 19 | + "### A few warm-up examples\n", |
| 20 | + "Just to check that we're understanding both the theory and the Python implementation correctly..." |
| 21 | + ] |
| 22 | + }, |
| 23 | + { |
| 24 | + "cell_type": "code", |
| 25 | + "execution_count": null, |
| 26 | + "metadata": {}, |
| 27 | + "outputs": [], |
| 28 | + "source": [ |
| 29 | + "import numpy as np" |
| 30 | + ] |
| 31 | + }, |
| 32 | + { |
| 33 | + "cell_type": "markdown", |
| 34 | + "metadata": {}, |
| 35 | + "source": [ |
| 36 | + "A single-precision float has a 24-bit significand. This gives it a relative precision of 1 in $2^{24}$, or about 1 in 16 million:" |
| 37 | + ] |
| 38 | + }, |
| 39 | + { |
| 40 | + "cell_type": "code", |
| 41 | + "execution_count": null, |
| 42 | + "metadata": {}, |
| 43 | + "outputs": [], |
| 44 | + "source": [ |
| 45 | + "print(\"1 in {0} ({0:g})\".format(2**(24)))\n", |
| 46 | + "print(\"{:.24f}\".format(1./2**(24)))" |
| 47 | + ] |
| 48 | + }, |
| 49 | + { |
| 50 | + "cell_type": "markdown", |
| 51 | + "metadata": {}, |
| 52 | + "source": [ |
| 53 | + "Note that double-floats have a relative precision of 1 in $2^{53}$, i.e." |
| 54 | + ] |
| 55 | + }, |
| 56 | + { |
| 57 | + "cell_type": "code", |
| 58 | + "execution_count": null, |
| 59 | + "metadata": {}, |
| 60 | + "outputs": [], |
| 61 | + "source": [ |
| 62 | + "print(\"1 in {0} ({0:g})\".format(2**(53)))" |
| 63 | + ] |
| 64 | + }, |
| 65 | + { |
| 66 | + "cell_type": "markdown", |
| 67 | + "metadata": {}, |
| 68 | + "source": [ |
| 69 | + "This has some basic effects on the accuracy of addition (see also [Kahan summation](https://en.wikipedia.org/wiki/Kahan_summation_algorithm)):" |
| 70 | + ] |
| 71 | + }, |
| 72 | + { |
| 73 | + "cell_type": "code", |
| 74 | + "execution_count": null, |
| 75 | + "metadata": {}, |
| 76 | + "outputs": [], |
| 77 | + "source": [ |
| 78 | + "2**24, 2**24+1, 2**24+0.5 #Double precision, gives exact results for these examples" |
| 79 | + ] |
| 80 | + }, |
| 81 | + { |
| 82 | + "cell_type": "code", |
| 83 | + "execution_count": null, |
| 84 | + "metadata": {}, |
| 85 | + "outputs": [], |
| 86 | + "source": [ |
| 87 | + "np.float32(2**24) + np.float32(1) #Single precision, gets rounded down" |
| 88 | + ] |
| 89 | + }, |
| 90 | + { |
| 91 | + "cell_type": "code", |
| 92 | + "execution_count": null, |
| 93 | + "metadata": {}, |
| 94 | + "outputs": [], |
| 95 | + "source": [ |
| 96 | + "np.float32(2**24) + np.float32(2) #Beyond 2**24 the smallest rounding unit is 2" |
| 97 | + ] |
| 98 | + }, |
| 99 | + { |
| 100 | + "cell_type": "markdown", |
| 101 | + "metadata": {}, |
| 102 | + "source": [ |
| 103 | + "Of course, this scales up and down as the exponent varies:" |
| 104 | + ] |
| 105 | + }, |
| 106 | + { |
| 107 | + "cell_type": "code", |
| 108 | + "execution_count": null, |
| 109 | + "metadata": {}, |
| 110 | + "outputs": [], |
| 111 | + "source": [ |
| 112 | + "2**21 + 1./2**3 # double precision, gives exact results for these examples" |
| 113 | + ] |
| 114 | + }, |
| 115 | + { |
| 116 | + "cell_type": "code", |
| 117 | + "execution_count": null, |
| 118 | + "metadata": {}, |
| 119 | + "outputs": [], |
| 120 | + "source": [ |
| 121 | + "np.float32(2**21) + np.float32(1./2**3)" |
| 122 | + ] |
| 123 | + }, |
| 124 | + { |
| 125 | + "cell_type": "code", |
| 126 | + "execution_count": null, |
| 127 | + "metadata": {}, |
| 128 | + "outputs": [], |
| 129 | + "source": [ |
| 130 | + "1./2**21" |
| 131 | + ] |
| 132 | + }, |
| 133 | + { |
| 134 | + "cell_type": "code", |
| 135 | + "execution_count": null, |
| 136 | + "metadata": {}, |
| 137 | + "outputs": [], |
| 138 | + "source": [ |
| 139 | + "2**2" |
| 140 | + ] |
| 141 | + }, |
| 142 | + { |
| 143 | + "cell_type": "code", |
| 144 | + "execution_count": null, |
| 145 | + "metadata": {}, |
| 146 | + "outputs": [], |
| 147 | + "source": [ |
| 148 | + "print(\"{:.24g}\".format(np.float32(2**2) + np.float32(1./2**21)))" |
| 149 | + ] |
| 150 | + }, |
| 151 | + { |
| 152 | + "cell_type": "markdown", |
| 153 | + "metadata": {}, |
| 154 | + "source": [ |
| 155 | + "### Beware of aggregates\n", |
| 156 | + "One obvious example for transient detection is averaging of visibilities. Assume for now that earlier parts of the data-acquisition / reduction process are sufficiently high-precision that we can treat them as real-number operations, but we drop down to single-float precision for visibility storage. \n", |
| 157 | + "\n", |
| 158 | + "We'll simply treat the combined signal as a real number, with the assumption that the important bits will generalise to complex components.\n", |
| 159 | + "\n", |
| 160 | + "Suppose we take $2^4=32$ integrations and then attempt to average them. Now add a transient signal to one of those integrations. What is the flux of a transient that's on the borderline of being lost due to numerical precision considerations? For single-precision with a ratio of $2^{24}$, we get a borderline flux of $2^{24-4}=2^{20}$, for a dynamic range of around a million:" |
| 161 | + ] |
| 162 | + }, |
| 163 | + { |
| 164 | + "cell_type": "code", |
| 165 | + "execution_count": null, |
| 166 | + "metadata": {}, |
| 167 | + "outputs": [], |
| 168 | + "source": [ |
| 169 | + "2**20" |
| 170 | + ] |
| 171 | + }, |
| 172 | + { |
| 173 | + "cell_type": "markdown", |
| 174 | + "metadata": {}, |
| 175 | + "source": [ |
| 176 | + "So for a steady source flux of 1Jy, the borderline transient source fluence (in the careless / naive case) is now $\\approx$1 $\\mu$Jy:" |
| 177 | + ] |
| 178 | + }, |
| 179 | + { |
| 180 | + "cell_type": "code", |
| 181 | + "execution_count": null, |
| 182 | + "metadata": {}, |
| 183 | + "outputs": [], |
| 184 | + "source": [ |
| 185 | + "1./2**20" |
| 186 | + ] |
| 187 | + }, |
| 188 | + { |
| 189 | + "cell_type": "code", |
| 190 | + "execution_count": null, |
| 191 | + "metadata": {}, |
| 192 | + "outputs": [], |
| 193 | + "source": [ |
| 194 | + "steady_src_flux = 1.\n", |
| 195 | + "n_integrations = 2**4\n", |
| 196 | + "steady_data = steady_src_flux*np.ones(n_integrations, dtype=np.float64)\n", |
| 197 | + "transient_signal = np.zeros_like(steady_data)\n", |
| 198 | + "transient_fluence = 1./2**20\n", |
| 199 | + "transient_signal[0] = transient_fluence\n", |
| 200 | + "print(transient_signal, transient_signal.dtype) " |
| 201 | + ] |
| 202 | + }, |
| 203 | + { |
| 204 | + "cell_type": "code", |
| 205 | + "execution_count": null, |
| 206 | + "metadata": {}, |
| 207 | + "outputs": [], |
| 208 | + "source": [ |
| 209 | + "data = steady_data + transient_signal\n", |
| 210 | + "data" |
| 211 | + ] |
| 212 | + }, |
| 213 | + { |
| 214 | + "cell_type": "markdown", |
| 215 | + "metadata": {}, |
| 216 | + "source": [ |
| 217 | + "Let's recover the transient flux:" |
| 218 | + ] |
| 219 | + }, |
| 220 | + { |
| 221 | + "cell_type": "code", |
| 222 | + "execution_count": null, |
| 223 | + "metadata": {}, |
| 224 | + "outputs": [], |
| 225 | + "source": [ |
| 226 | + "np.mean(data) - steady_src_flux" |
| 227 | + ] |
| 228 | + }, |
| 229 | + { |
| 230 | + "cell_type": "markdown", |
| 231 | + "metadata": {}, |
| 232 | + "source": [ |
| 233 | + "And now working with single-precision:" |
| 234 | + ] |
| 235 | + }, |
| 236 | + { |
| 237 | + "cell_type": "code", |
| 238 | + "execution_count": null, |
| 239 | + "metadata": {}, |
| 240 | + "outputs": [], |
| 241 | + "source": [ |
| 242 | + "single_prec_data = np.asarray(data, dtype=np.float32)\n", |
| 243 | + "single_prec_data " |
| 244 | + ] |
| 245 | + }, |
| 246 | + { |
| 247 | + "cell_type": "code", |
| 248 | + "execution_count": null, |
| 249 | + "metadata": {}, |
| 250 | + "outputs": [], |
| 251 | + "source": [ |
| 252 | + "np.mean(single_prec_data)# Lost the transient to numerical precision" |
| 253 | + ] |
| 254 | + }, |
| 255 | + { |
| 256 | + "cell_type": "code", |
| 257 | + "execution_count": null, |
| 258 | + "metadata": {}, |
| 259 | + "outputs": [], |
| 260 | + "source": [ |
| 261 | + "np.mean(single_prec_data) - steady_src_flux" |
| 262 | + ] |
| 263 | + }, |
| 264 | + { |
| 265 | + "cell_type": "markdown", |
| 266 | + "metadata": {}, |
| 267 | + "source": [ |
| 268 | + "So, if we were searching for the faint transient in the time-averaged (single-precision) image, we'd have no hope - the averaging process has diluted the flux to below the $\\frac{1}{2^{24}}$ precision-ratio limit. The takeaway being that we need to think in terms of **fluence ratios** (or equivalently, time-averaged flux-ratios) for the search-averaging period." |
| 269 | + ] |
| 270 | + }, |
| 271 | + { |
| 272 | + "cell_type": "markdown", |
| 273 | + "metadata": {}, |
| 274 | + "source": [ |
| 275 | + "If we search for transient on the timescale of a single integration, we'll still detect it, so simply need to beware of aggregate statistics:" |
| 276 | + ] |
| 277 | + }, |
| 278 | + { |
| 279 | + "cell_type": "code", |
| 280 | + "execution_count": null, |
| 281 | + "metadata": {}, |
| 282 | + "outputs": [], |
| 283 | + "source": [ |
| 284 | + "residuals = single_prec_data - np.mean(single_prec_data)\n", |
| 285 | + "residuals" |
| 286 | + ] |
| 287 | + }, |
| 288 | + { |
| 289 | + "cell_type": "markdown", |
| 290 | + "metadata": {}, |
| 291 | + "source": [ |
| 292 | + "Note, we could alternatively use higher-precision for aggregate computation:" |
| 293 | + ] |
| 294 | + }, |
| 295 | + { |
| 296 | + "cell_type": "code", |
| 297 | + "execution_count": null, |
| 298 | + "metadata": {}, |
| 299 | + "outputs": [], |
| 300 | + "source": [ |
| 301 | + "np.mean(single_prec_data, dtype=np.float64)" |
| 302 | + ] |
| 303 | + }, |
| 304 | + { |
| 305 | + "cell_type": "code", |
| 306 | + "execution_count": null, |
| 307 | + "metadata": {}, |
| 308 | + "outputs": [], |
| 309 | + "source": [ |
| 310 | + "# For comparison:\n", |
| 311 | + "transient_fluence / n_integrations" |
| 312 | + ] |
| 313 | + }, |
| 314 | + { |
| 315 | + "cell_type": "markdown", |
| 316 | + "metadata": {}, |
| 317 | + "source": [ |
| 318 | + "## Summary:\n", |
| 319 | + "For single-precision floating point, the 'instantaneous' dynamic range is 1 in 16 million. This is preserved so long as care is taken when dealing with aggregate quantities...\n", |
| 320 | + "\n", |
| 321 | + "... Another example of which is the fourier transform (see next notebook)." |
| 322 | + ] |
| 323 | + } |
| 324 | + ], |
| 325 | + "metadata": { |
| 326 | + "language_info": { |
| 327 | + "name": "python", |
| 328 | + "pygments_lexer": "ipython3" |
| 329 | + } |
| 330 | + }, |
| 331 | + "nbformat": 4, |
| 332 | + "nbformat_minor": 1 |
| 333 | +} |
0 commit comments