Coverage for src/ufig/plugins/gamma_interpolation_table.py: 92%
54 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-12 19:08 +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 )
42 filename = filename + ".fits"
43 resource_directory = resource_filename(ufig.__name__, "res/intrinsictables/")
44 load_filename = os.path.join(resource_directory, filename)
45 local_filename = os.path.join(os.getcwd(), filename)
47 # copy to local directory (local scratch) if not already there
48 if copy_to_cwd and (os.path.exists(local_filename) is False):
49 src = os.path.join(resource_directory, filename)
50 file_utils.robust_copy(src, local_filename)
51 load_filename = local_filename
52 elif copy_to_cwd:
53 load_filename = local_filename
55 try:
56 intrinsicTable = fits.getdata(load_filename, ext=0).astype(np.float32)
57 except OSError:
58 raise OSError(
59 "Include gamma_interpolation_table-module or provide matching parameter"
60 " values to pre-computed tables"
61 ) from None
63 return intrinsicTable
66def compute_intrinsictable(sersicprecision, gammaprecision, gammaprecisionhigh):
67 """
68 Compute the table containing values of the inverse gamma function.
70 :param sersicprecision: 2**sersicpresicion points to interpolate the Sersic range
71 [0,10]
72 :param gammaprecision: 2*(2**gammaprecision) to interpolate the range [0,1[ where
73 the inverse cdf is evaluated
74 :param gammaprecisionhigh: Change of point density at 1-2**(-gammaprecisionhigh)
75 (increased precision close to 1)
76 :return intrinsicTable: Table containing values
77 """
79 # gamma lookup has 1<<sersicprecision, and a more precise fit on 1-1/(1<<3) also
80 # with 1<<sersicprecision elements
81 intrinsicTable = np.empty(
82 (
83 (1 << sersicprecision) + 1,
84 (1 << (gammaprecision + 1)) + 1,
85 ),
86 dtype=np.float32,
87 )
88 intrinsicTable[0, 0 : (1 << gammaprecision)] = (
89 np.power(
90 gammaincinv(
91 2e-15,
92 np.float32(range(0, 1 << gammaprecision))
93 / np.float32(1 << gammaprecision),
94 ),
95 1e-15,
96 )
97 / 1e-15
98 ).astype(np.float32)
99 intrinsicTable[0, (1 << gammaprecision) : (1 << (gammaprecision + 1))] = (
100 np.power(
101 gammaincinv(
102 2e-15,
103 1.0
104 - 1.0 / np.float32(1 << gammaprecisionhigh)
105 + np.float32(range(0, 1 << gammaprecision))
106 / np.float32(1 << (gammaprecision + gammaprecisionhigh)),
107 ),
108 1e-15,
109 )
110 / 1e-15
111 ).astype(np.float32)
112 intrinsicTable[0, 1 << (gammaprecision + 1)] = (
113 np.power(gammaincinv(2e-15, 1.0 - 1e-6), 1e-15) / 1e-15
114 ).astype(np.float32)
116 for i in range(1, (1 << sersicprecision) + 1):
117 n = 10.0 * np.float32(i << (32 - sersicprecision)) / (np.int64(1) << 32)
118 k = gammaincinv(2.0 * n, 0.5)
119 intrinsicTable[i, 0 : (1 << gammaprecision)] = (
120 np.power(
121 gammaincinv(
122 2.0 * n,
123 np.float32(range(0, 1 << gammaprecision))
124 / np.float32(1 << gammaprecision),
125 ),
126 n,
127 )
128 / np.power(k, n)
129 ).astype(np.float32)
130 intrinsicTable[i, (1 << gammaprecision) : (1 << (gammaprecision + 1))] = (
131 np.power(
132 gammaincinv(
133 2.0 * n,
134 1.0
135 - 1.0 / np.float32(1 << gammaprecisionhigh)
136 + np.float32(range(0, 1 << gammaprecision))
137 / np.float32(1 << (gammaprecision + gammaprecisionhigh)),
138 ),
139 n,
140 )
141 / np.power(k, n)
142 ).astype(np.float32)
143 intrinsicTable[i, 1 << (gammaprecision + 1)] = (
144 np.power(gammaincinv(2.0 * n, 1.0 - 1e-6), n) / np.power(k, n)
145 ).astype(np.float32)
147 return intrinsicTable
150class Plugin(BasePlugin):
151 """
152 Loads, or in case explicitly set or non-existent computes on-the-fly, the table with
153 pre-computed values of the inverse gamma function in the range [0,1[ for different
154 Sersic indexes in [0,10[ used in the render_galaxies.py-module to sample the radial
155 cdf.
157 Optionally a new intrinsic_table can furthermore be saved. Desirably not overwriting
158 the ones in res/intrinsictables/ unless specifically stated.
160 :param sersicprecision: 2**sersicpresicion points to interpolate the Sersic range
161 [0,10]
162 :param gammaprecision: 2*(2**gammaprecision) to interpolate the range [0,1[ where
163 the inverse cdf is evaluated
164 :param gammaprecisionhigh: Change of point density at 1-2**(-gammaprecisionhigh)
165 (increased precision close to 1)
166 :param compute_gamma_table_onthefly: whether or not the interpolation table is
167 computed on the fly
169 :return: Table containing values of the inv gamma function for different Sersic
170 values in [0,10] and values in [0,1[
172 """
174 def __call__(self):
175 if self.ctx.parameters.compute_gamma_table_onthefly:
176 self.ctx.intrinsicTable = compute_intrinsictable(
177 self.ctx.parameters.sersicprecision,
178 self.ctx.parameters.gammaprecision,
179 self.ctx.parameters.gammaprecisionhigh,
180 )
181 else:
182 self.ctx.intrinsicTable = load_intrinsicTable(
183 self.ctx.parameters.sersicprecision,
184 self.ctx.parameters.gammaprecision,
185 self.ctx.parameters.gammaprecisionhigh,
186 self.ctx.parameters.copy_gamma_table_to_cwd,
187 )
189 def __str__(self):
190 return NAME
192 def save_new_interpolation_table(self):
193 filename = (
194 f"intrinsic_table_{self.ctx.parameters.sersicprecision}"
195 f"_{self.ctx.parameters.gammaprecision}"
196 f"_{self.ctx.parameters.gammaprecisionhigh}"
197 )
198 resource_directory = resource_filename(ufig.__name__, "res/intrinsictables/")
200 hdu = fits.PrimaryHDU(data=self.ctx.intrinsicTable)
201 try:
202 hdu.writeto(resource_directory + filename + ".fits")
203 except OSError:
204 os.remove(resource_directory + filename + ".fits")
205 hdu.writeto(resource_directory + filename + ".fits")