Coverage for src/ufig/plugins/gamma_interpolation_table.py: 90%
58 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-17 08:48 +0000
1# Copyright (c) 2014 ETH Zurich, Institute of Astronomy, Claudio Bruderer
2# <claudio.bruderer@phys.ethz.ch>
3"""
4Created on Nov 24, 2013
5@author: Claudio Bruderer
7"""
9import os
11import numpy as np
12from astropy.io import fits
13from cosmic_toolbox import file_utils
14from ivy.plugin.base_plugin import BasePlugin
15from pkg_resources import resource_filename
16from scipy.special import gammaincinv
18import ufig
20NAME = "gamma interpolation table"
23def load_intrinsicTable(
24 sersicprecision, gammaprecision, gammaprecisionhigh, copy_to_cwd=False
25):
26 """
27 Load the table containing values of the inverse gamma function.
29 :param sersicprecision: 2**sersicpresicion points to interpolate the Sersic range
30 [0,10]
31 :param gammaprecision: 2*(2**gammaprecision) to interpolate the range [0,1[ where
32 the inverse cdf is evaluated
33 :param gammaprecisionhigh: Change of point density at 1-2**(-gammaprecisionhigh)
34 (increased precision close to 1)
35 :param copy_to_cwd: Copy the table to the current working directory
36 :return intrinsicTable: Table containing values
37 """
39 filename = (
40 f"intrinsic_table_{sersicprecision}_{gammaprecision}_{gammaprecisionhigh}"
41 )
43 if np.__version__ >= "2.0": 43 ↛ 45line 43 didn't jump to line 45 because the condition on line 43 was always true
44 filename += "_np2"
45 filename = filename + ".fits"
46 resource_directory = resource_filename(ufig.__name__, "res/intrinsictables/")
47 load_filename = os.path.join(resource_directory, filename)
48 local_filename = os.path.join(os.getcwd(), filename)
50 # copy to local directory (local scratch) if not already there
51 if copy_to_cwd and (os.path.exists(local_filename) is False):
52 src = os.path.join(resource_directory, filename)
53 file_utils.robust_copy(src, local_filename)
54 load_filename = local_filename
55 elif copy_to_cwd:
56 load_filename = local_filename
58 try:
59 intrinsicTable = fits.getdata(load_filename, ext=0).astype(np.float32)
60 except OSError:
61 raise OSError(
62 "Include gamma_interpolation_table-module or provide matching parameter"
63 " values to pre-computed tables"
64 ) from None
66 return intrinsicTable
69def compute_intrinsictable(sersicprecision, gammaprecision, gammaprecisionhigh):
70 """
71 Compute the table containing values of the inverse gamma function.
73 :param sersicprecision: 2**sersicpresicion points to interpolate the Sersic range
74 [0,10]
75 :param gammaprecision: 2*(2**gammaprecision) to interpolate the range [0,1[ where
76 the inverse cdf is evaluated
77 :param gammaprecisionhigh: Change of point density at 1-2**(-gammaprecisionhigh)
78 (increased precision close to 1)
79 :return intrinsicTable: Table containing values
80 """
82 # gamma lookup has 1<<sersicprecision, and a more precise fit on 1-1/(1<<3) also
83 # with 1<<sersicprecision elements
84 intrinsicTable = np.empty(
85 (
86 (1 << sersicprecision) + 1,
87 (1 << (gammaprecision + 1)) + 1,
88 ),
89 dtype=np.float32,
90 )
91 intrinsicTable[0, 0 : (1 << gammaprecision)] = (
92 np.power(
93 gammaincinv(
94 2e-15,
95 np.float32(range(0, 1 << gammaprecision))
96 / np.float32(1 << gammaprecision),
97 ),
98 1e-15,
99 )
100 / 1e-15
101 ).astype(np.float32)
102 intrinsicTable[0, (1 << gammaprecision) : (1 << (gammaprecision + 1))] = (
103 np.power(
104 gammaincinv(
105 2e-15,
106 1.0
107 - 1.0 / np.float32(1 << gammaprecisionhigh)
108 + np.float32(range(0, 1 << gammaprecision))
109 / np.float32(1 << (gammaprecision + gammaprecisionhigh)),
110 ),
111 1e-15,
112 )
113 / 1e-15
114 ).astype(np.float32)
115 intrinsicTable[0, 1 << (gammaprecision + 1)] = (
116 np.power(gammaincinv(2e-15, 1.0 - 1e-6), 1e-15) / 1e-15
117 ).astype(np.float32)
119 for i in range(1, (1 << sersicprecision) + 1):
120 n = 10.0 * np.float32(i << (32 - sersicprecision)) / (np.int64(1) << 32)
121 k = gammaincinv(2.0 * n, 0.5)
122 intrinsicTable[i, 0 : (1 << gammaprecision)] = (
123 np.power(
124 gammaincinv(
125 2.0 * n,
126 np.float32(range(0, 1 << gammaprecision))
127 / np.float32(1 << gammaprecision),
128 ),
129 n,
130 )
131 / np.power(k, n)
132 ).astype(np.float32)
133 intrinsicTable[i, (1 << gammaprecision) : (1 << (gammaprecision + 1))] = (
134 np.power(
135 gammaincinv(
136 2.0 * n,
137 1.0
138 - 1.0 / np.float32(1 << gammaprecisionhigh)
139 + np.float32(range(0, 1 << gammaprecision))
140 / np.float32(1 << (gammaprecision + gammaprecisionhigh)),
141 ),
142 n,
143 )
144 / np.power(k, n)
145 ).astype(np.float32)
146 intrinsicTable[i, 1 << (gammaprecision + 1)] = (
147 np.power(gammaincinv(2.0 * n, 1.0 - 1e-6), n) / np.power(k, n)
148 ).astype(np.float32)
150 return intrinsicTable
153class Plugin(BasePlugin):
154 """
155 Loads, or in case explicitly set or non-existent computes on-the-fly, the table with
156 pre-computed values of the inverse gamma function in the range [0,1[ for different
157 Sersic indexes in [0,10[ used in the render_galaxies.py-module to sample the radial
158 cdf.
160 Optionally a new intrinsic_table can furthermore be saved. Desirably not overwriting
161 the ones in res/intrinsictables/ unless specifically stated.
163 :param sersicprecision: 2**sersicpresicion points to interpolate the Sersic range
164 [0,10]
165 :param gammaprecision: 2*(2**gammaprecision) to interpolate the range [0,1[ where
166 the inverse cdf is evaluated
167 :param gammaprecisionhigh: Change of point density at 1-2**(-gammaprecisionhigh)
168 (increased precision close to 1)
169 :param compute_gamma_table_onthefly: whether or not the interpolation table is
170 computed on the fly
172 :return: Table containing values of the inv gamma function for different Sersic
173 values in [0,10] and values in [0,1[
175 """
177 def __call__(self):
178 if self.ctx.parameters.compute_gamma_table_onthefly:
179 self.ctx.intrinsicTable = compute_intrinsictable(
180 self.ctx.parameters.sersicprecision,
181 self.ctx.parameters.gammaprecision,
182 self.ctx.parameters.gammaprecisionhigh,
183 )
184 else:
185 self.ctx.intrinsicTable = load_intrinsicTable(
186 self.ctx.parameters.sersicprecision,
187 self.ctx.parameters.gammaprecision,
188 self.ctx.parameters.gammaprecisionhigh,
189 self.ctx.parameters.copy_gamma_table_to_cwd,
190 )
192 def __str__(self):
193 return NAME
195 def save_new_interpolation_table(self):
196 filename = (
197 f"intrinsic_table_{self.ctx.parameters.sersicprecision}"
198 f"_{self.ctx.parameters.gammaprecision}"
199 f"_{self.ctx.parameters.gammaprecisionhigh}"
200 )
201 if np.__version__ >= "2.0": 201 ↛ 203line 201 didn't jump to line 203 because the condition on line 201 was always true
202 filename += "_np2"
203 resource_directory = resource_filename(ufig.__name__, "res/intrinsictables/")
205 hdu = fits.PrimaryHDU(data=self.ctx.intrinsicTable)
206 try:
207 hdu.writeto(resource_directory + filename + ".fits")
208 except OSError:
209 os.remove(resource_directory + filename + ".fits")
210 hdu.writeto(resource_directory + filename + ".fits")