Skip to content

nonparametric

bin_map(hdu, rbins, x0=None, y0=None, cunit=None)

Radially bin a map into rbins. Code adapted from CLASS

Parameters:

Name Type Description Default
hdu HDUList

hdu containing map to bin

 required
rbins NDArray[floating]

Bin edges in radians

 required
cunit Union[None, np.floating], Default: None

Pixel units. If None, will atempt to infer from imap

 None

Returns:

Name Type Description
bin1d NDArray[floating]

Bin center values
var1d NDArray[floating]

Bin variance estimate
Source code in witch/nonparametric.py 
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
def bin_map(hdu, rbins, x0=None, y0=None, cunit=None):
    """
    Radially bin a map into rbins. Code adapted from CLASS

    Parameters
    ----------
    hdu : fits.HDUList
        hdu containing map to bin
    rbins : NDArray[np.floating]
        Bin edges in radians
    cunit : Union[None, np.floating], Default: None
        Pixel units. If None, will atempt to infer from imap

    Returns
    -------
    bin1d : NDArray[np.floating]
        Bin center values
    var1d : NDArray[np.floating]
        Bin variance estimate
    """

    if cunit is None:
        try:
            cunit = hdu[0].header["CUNIT1"].lower()
        except KeyError as e:
            raise e

    if (
        cunit.lower() == "rad"
        or cunit.lower() == "radian"
        or cunit.lower() == "radians"
    ):
        pixunits = 1
    elif (
        cunit.lower() == "deg"
        or cunit.lower() == "degree"
        or cunit.lower() == "degrees"
    ):
        pixunits = wu.rad_to_deg
    elif (
        cunit.lower() == "arcmin"
        or cunit.lower() == "arcminute"
        or cunit.lower() == "arcminutes"
    ):
        pixunits = wu.rad_to_arcmin
    elif (
        cunit.lower() == "arcsec"
        or cunit.lower() == "arcsecond"
        or cunit.lower() == "arcseconds"
    ):
        pixunits = wu.rad_to_arcsec
    else:
        raise ValueError("Error: cunit {} is not a valid pixel unit".format(cunit))

    pixsize = np.abs(hdu[0].header["CDELT1"]) / pixunits
    x0 = hdu[0].header["CRVAL1"] / pixunits
    y0 = hdu[0].header["CRVAL2"] / pixunits

    if np.abs(hdu[0].header["CDELT1"]) != np.abs(hdu[0].header["CDELT2"]):
        warnings.warn(
            "Warning: non-square pixels: RA: {} Dec{}".format(
                np.abs(hdu[0].header["CDELT1"]), np.abs(hdu[0].header["CDELT2"])
            )
        )

    # The offset is redundent if the binning center is taken to be the map center but frequently it is not
    x = np.linspace(
        -hdu[0].data.shape[1] / 2 * pixsize + hdu[0].header["CRVAL1"] / pixunits,
        hdu[0].data.shape[1] / 2 * pixsize + hdu[0].header["CRVAL1"] / pixunits,
        hdu[0].data.shape[1],
    )
    y = np.linspace(
        -hdu[0].data.shape[0] / 2 * pixsize + hdu[0].header["CRVAL2"] / pixunits,
        hdu[0].data.shape[0] / 2 * pixsize + hdu[0].header["CRVAL2"] / pixunits,
        hdu[0].data.shape[0],
    )

    X, Y = np.meshgrid(x, y)
    R = np.sqrt((X - x0) ** 2 + (Y - y0) ** 2)
    rbins = rbins.append(999999)
    bin1d = np.zeros(len(rbins) - 1)
    var1d = np.zeros(len(rbins) - 1)

    for k in range(len(rbins) - 1):
        pixels = [
            hdu[0].data[i, j]
            for i in range(len(y))
            for j in range(len(x))
            if rbins[k] < R[i, j] <= rbins[k + 1]
        ]
        bin1d[k] = np.mean(pixels)
        var1d[k] = np.var(pixels)

    return bin1d, var1d

broken_power(rs, condlist, rbins, amps, pows, c)

Function which returns a broken powerlaw evaluated at rs.

Parameters:

rs : jax.Array Array of rs at which to compute pl. condlist : tuple tuple which enocdes which rs are evaluated by which parametric function rbins : jax.Array Array of bin edges for power laws amps : jax.Array Amplitudes of power laws pows : jax.Array Exponents of power laws c : float Constant offset for powerlaws

Source code in witch/nonparametric.py 
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
@jax.jit
def broken_power(
    rs: jax.Array,
    condlist: tuple,
    rbins: jax.Array,
    amps: jax.Array,
    pows: jax.Array,
    c: float,
) -> jax.Array:
    """
    Function which returns a broken powerlaw evaluated at rs.

    Parameters:
    -----------
    rs : jax.Array
        Array of rs at which to compute pl.
    condlist : tuple
        tuple which enocdes which rs are evaluated by which parametric function
    rbins : jax.Array
        Array of bin edges for power laws
    amps : jax.Array
        Amplitudes of power laws
    pows : jax.Array                                                                                                                                                                                                                                                                            Exponents of power laws
    c : float
        Constant offset for powerlaws
    """
    cur_c = c  # TODO: necessary?
    funclist = []
    for i in range(len(condlist) - 1, -1, -1):
        funclist.append(
            partial(power, rbin=rbins[i + 1], cur_amp=amps[i], cur_pow=pows[i], c=cur_c)
        )
        cur_c += amps[i] * (rbins[i] ** pows[i] - rbins[i + 1] ** pows[i])
    return jnp.piecewise(rs, condlist, funclist)

get_rbins(model, rmax=3.0 * 60.0, struct_num=0, sig_params=['amp', 'P0'], default=(0, 10, 20, 30, 50, 80, 120, 180))

Function which returns a good set of rbins for a non-parametric fit given the significance of the underlying parametric model.

Parameters:

Name Type Description Default
model Model

Parametric model to calculate rbins on

 required
rmax float

Maximum radius of the rbins

 180
struct_num int, defualt: 0

Structure within model to calculate rbins on

 0
sig_params list[str]

Parameters to consider for computing significance. Only first match will be used.

 ['amp', 'P0']
default tuple[int]

Default rbins to be returned if generation fails.

 (0, 10, 20, 30, 50, 80, 120, 180)

Returns:

Name Type Description
rbins tuple[int]

rbins for nonparametric fit
Source code in witch/nonparametric.py 
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
def get_rbins(
    model,
    rmax: float = 3.0 * 60.0,
    struct_num: int = 0,
    sig_params: list[str] = ["amp", "P0"],
    default: tuple[int] = (0, 10, 20, 30, 50, 80, 120, 180),
) -> tuple[int]:
    """
    Function which returns a good set of rbins for a non-parametric fit given the significance of the underlying parametric model.

    Parameters
    ----------
    model : container.Model
        Parametric model to calculate rbins on
    rmax : float, default: 180
        Maximum radius of the rbins
    struct_num : int, defualt: 0
        Structure within model to calculate rbins on
    sig_params: list[str], default: ["amp", "P0"]
        Parameters to consider for computing significance.
        Only first match will be used.
    default: tuple[int], default: (0, 10, 20, 30, 50, 80, 120, 180)
        Default rbins to be returned if generation fails.

    Returns
    -------
    rbins: tuple[int]
        rbins for nonparametric fit
    """
    sig = 0
    for par in model.structures[struct_num].parameters:
        if par.name in sig_params:
            sig = par.val / par.err
            break
    if sig == 0:
        warnings.warn(
            "Warning: model does not contain any valid significance parameters {}. Returning default bins.".format(
                sig_params
            )
        )
        return default

    if sig < 10:
        warnings.warn(
            "Warning, significance {} too low to calculate bins. Returning default bins.".format(
                sig
            )
        )
        return default

    rbins = [0, 10, 20]
    rmin = 30
    nrbins = int(np.floor(sig / 5)[0] - 3)

    if nrbins == 1:
        rbins = np.array(rbins)
        rbins = np.append(rbins, rmax)

        return rbins

    logrange = np.logspace(np.log10(rmin), np.log10(rmax), nrbins)
    step = logrange[1] - logrange[0]

    while step < 10:
        rbins.append(rmin)
        rmin += 10
        nrbins -= 1
        logrange = np.logspace(np.log10(rmin), np.log10(rmax), nrbins)
        step = logrange[1] - logrange[0]
        if rmin > rmax or nrbins < 1:
            return np.array(rbins)
    rbins = np.array(rbins)
    rbins = np.append(rbins, logrange)

    return tuple(rbins)

power(x, rbin, cur_amp, cur_pow, c)

Function which returns the powerlaw, given the bin-edge constraints. Exists to be partialed.

Parameters:

x : float Dummy variable to be partialed over rbin : float Edge of bin for powerlaw cur_amp : float Amplitude of power law cur_pow : float Power of power law c : float Constant offset

Returns:

Name Type Description
tmp float

Powerlaw evaluated at x
Source code in witch/nonparametric.py 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
@jax.jit
def power(x: float, rbin: float, cur_amp: float, cur_pow: float, c: float):
    """
    Function which returns the powerlaw, given the bin-edge constraints. Exists to be partialed.

    Parameters:
    -----------
    x : float
        Dummy variable to be partialed over
    rbin : float
        Edge of bin for powerlaw
    cur_amp : float
        Amplitude of power law
    cur_pow : float
        Power of power law
    c : float
        Constant offset

    Returns
    -------
    tmp : float
        Powerlaw evaluated at x
    """
    tmp = cur_amp * (x**cur_pow - rbin**cur_pow) + c
    return tmp

profile_to_broken_power(rs, ys, condlist, rbins)

Estimates a non-parametric broken power profile from a generic profile. Note this is an estimation only; in partciular since we fit piece-wise the c's get messed up. This broken powerlaw should then be fit to the data.

Parameters:

Name Type Description Default
rs ArrayLike

Array of radius values for the profile

 required
ys ArrayLike

Profile y values

 required
condlist list[ArrayLike]

List which defines which powerlaws map to which radii. See broken_power

 required
rbins ArrayLike

Array of bin edges defining the broken powerlaws

 required

Returns:

Name Type Description
amps array

Best fit amps for the powerlaws
pows array

Best fit powers for the powerlaws
c float

Best fit c for only the outermost powerlaw
Source code in witch/nonparametric.py 
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
def profile_to_broken_power(
    rs: ArrayLike, ys: ArrayLike, condlist: list[ArrayLike], rbins: ArrayLike
) -> tuple[jnp.array, jnp.array, float]:
    """
    Estimates a non-parametric broken power profile from a generic profile.
    Note this is an estimation only; in partciular since we fit piece-wise
    the c's get messed up. This broken powerlaw should then be fit to the
    data.

    Parameters
    ----------
    rs : ArrayLike
        Array of radius values for the profile
    ys : ArrayLike
        Profile y values
    condlist : list[ArrayLike]
        List which defines which powerlaws map to which radii. See broken_power
    rbins : ArrayLike
        Array of bin edges defining the broken powerlaws

    Returns
    -------
    amps : jnp.array
        Best fit amps for the powerlaws
    pows : jnp.array
        Best fit powers for the powerlaws
    c : float
        Best fit c for only the outermost powerlaw
    """
    rs = jnp.array([x if x != 0 else 1e-1 for x in rs])  # Dont blow up

    rbins = jnp.array(
        [x if x != 0 else jnp.amin(rs) for x in rbins]
    )  # Dont blow up 2.0

    amps = jnp.zeros(len(condlist))
    pows = jnp.zeros(len(condlist))

    for i in range(len(condlist)):
        xdata = rs[condlist[i]]
        ydata = ys[condlist[i]]
        if i == len(condlist) - 1:
            popt, pcov = curve_fit(power, xdata, ydata, method="trf")
        else:
            popt, pcov = curve_fit(
                power,
                xdata,
                ydata,
                method="trf",
                p0=[rbins[::-1][i], np.amax(ydata) * 1e5, -2, 0.0],
            )
        if i == 0:
            c = popt[3]
        amps = amps.at[i].set(popt[1])
        pows = pows.at[i].set(popt[2])

    return amps[::-1], pows[::-1], c