Skip to content

Commit 6c55e8e

Browse files
Cameron-Van-EckjackieykmaAlecThomson
authored
New qufitting models (#114)
* New QU-fitting models, and minor bug fix (#72) * Improved model files to be more descriptive, and added 3 new models * Fixed error message not showing correct path to model files * Added triple Faraday thin model * Fixed indentation error * Changed fracpol minimum to 0.0 * Changed fracpol minimum to 0.0 * Changed fracpol minimum to 0.0 * Changed fracpol minimum to 0.0 * Changed fracpol minimum to 0.0 * Changed fracpol minimum to 0.0 * Changed fracpol minimum to 0.0 * Changed fracpol minimum to 0.0 * Changed fracpol minimum to 0.0 * Test models * Standardise test names * Ignore * Testing tests * Verbose output * Updated new models to better LaTeX strings. Made math of models immune to black reformatting, for readability. * Re-introduced Alec's new QUfitting tests. Disabled the values tests, because they were too sensitive to variations between runs. Only tests if QU-fitting runs, but not the output values. * Fixing tests. --------- Co-authored-by: Yik Ki (Jackie) Ma <yikki.ma@anu.edu.au> Co-authored-by: Alec Thomson (S&A, Kensington WA) <alec.thomson@csiro.au>
1 parent 3860e5f commit 6c55e8e

11 files changed

Lines changed: 583 additions & 54 deletions

File tree

‎RMtools_1D/do_QUfit_1D_mnest.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -223,7 +223,7 @@ def load_model(modelNum, verbose=False):
223223
spec.loader.exec_module(mod)
224224
except:
225225
print(
226-
"Model could not be found! Please make sure model is present either in {}/models_ns/, or in {}/RMtools_1D/models_ns/".format(
226+
"Model could not be found! Please make sure model is present either in {}/models_ns/, or in {}/models_ns/".format(
227227
os.getcwd(), RMtools_dir
228228
)
229229
)

‎RMtools_1D/models_ns/m1.py‎

Lines changed: 13 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -13,13 +13,24 @@
1313
# quArr = Complex array containing the Re and Im spectra. #
1414
# -----------------------------------------------------------------------------#
1515
def model(pDict, lamSqArr_m2):
16-
"""Simple Faraday thin source."""
16+
"""
17+
18+
Simple Faraday thin source
19+
20+
Ref:
21+
Sokoloff et al. (1998) Eq 2
22+
O'Sullivan et al. (2012) Eq 8
23+
Ma et al. (2019a) Eq 10
24+
25+
"""
1726

1827
# Calculate the complex fractional q and u spectra
28+
# fmt: off
1929
pArr = pDict["fracPol"] * np.ones_like(lamSqArr_m2)
2030
quArr = pArr * np.exp(
2131
2j * (np.radians(pDict["psi0_deg"]) + pDict["RM_radm2"] * lamSqArr_m2)
2232
)
33+
# fmt: on
2334

2435
return quArr
2536

@@ -31,7 +42,7 @@ def model(pDict, lamSqArr_m2):
3142
# -----------------------------------------------------------------------------#
3243
priors = {
3344
"fracPol": bilby.prior.Uniform(
34-
minimum=0.001, maximum=1.0, name="fracPol", latex_label=r"$p$"
45+
minimum=0.0, maximum=1.0, name="fracPol", latex_label=r"$p$"
3546
),
3647
"psi0_deg": bilby.prior.Uniform(
3748
minimum=0,

‎RMtools_1D/models_ns/m11.py‎

Lines changed: 17 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -14,19 +14,30 @@
1414
# quArr = Complex array containing the Re and Im spectra. #
1515
# -----------------------------------------------------------------------------#
1616
def model(pDict, lamSqArr_m2):
17-
"""Two separate Faraday components, averaged within same telescope beam
18-
(i.e., unresolved), with a common Burn depolarisation term."""
17+
"""
18+
19+
Two separate Faraday thin sources
20+
Averaged within the same telescope beam (i.e., unresolved)
21+
22+
Ref (for individual source component):
23+
Sokoloff et al. (1998) Eq 2
24+
O'Sullivan et al. (2012) Eq 8
25+
Ma et al. (2019a) Eq 10
26+
27+
"""
1928

2029
# Calculate the complex fractional q and u spectra
2130
pArr1 = pDict["fracPol1"] * np.ones_like(lamSqArr_m2)
2231
pArr2 = pDict["fracPol2"] * np.ones_like(lamSqArr_m2)
32+
# fmt: off
2333
quArr1 = pArr1 * np.exp(
2434
2j * (np.radians(pDict["psi01_deg"]) + pDict["RM1_radm2"] * lamSqArr_m2)
2535
)
2636
quArr2 = pArr2 * np.exp(
2737
2j * (np.radians(pDict["psi02_deg"]) + pDict["RM2_radm2"] * lamSqArr_m2)
2838
)
2939
quArr = quArr1 + quArr2
40+
# fmt: on
3041

3142
return quArr
3243

@@ -61,13 +72,14 @@ def converter(parameters):
6172
priors = PriorDict(conversion_function=converter)
6273

6374
priors["fracPol1"] = bilby.prior.Uniform(
64-
minimum=0.001,
75+
minimum=0.0,
6576
maximum=1.0,
6677
name="fracPol1",
6778
latex_label=r"$p_1$",
6879
)
80+
6981
priors["fracPol2"] = bilby.prior.Uniform(
70-
minimum=0.001,
82+
minimum=0.0,
7183
maximum=1.0,
7284
name="fracPol2",
7385
latex_label=r"$p_2$",
@@ -107,7 +119,7 @@ def converter(parameters):
107119
latex_label=r"$\Delta\phi_{1,2}$ (rad m$^{-2}$)",
108120
)
109121
priors["sum_p1_p2"] = Constraint(
110-
minimum=0.001,
122+
minimum=0.0,
111123
maximum=1,
112124
name="sum_p1_p2",
113125
latex_label=r"$p_1+p_2$)",

‎RMtools_1D/models_ns/m111.py‎

Lines changed: 156 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,156 @@
1+
# =============================================================================#
2+
# MODEL DEFINITION FILE #
3+
# =============================================================================#
4+
import bilby
5+
import numpy as np
6+
from bilby.core.prior import Constraint, PriorDict
7+
8+
9+
# -----------------------------------------------------------------------------#
10+
# Function defining the model. #
11+
# #
12+
# pDict = Dictionary of parameters, created by parsing inParms, below. #
13+
# lamSqArr_m2 = Array of lambda-squared values #
14+
# quArr = Complex array containing the Re and Im spectra. #
15+
# -----------------------------------------------------------------------------#
16+
def model(pDict, lamSqArr_m2):
17+
"""
18+
19+
Three separate Faraday thin sources
20+
Averaged within the same telescope beam (i.e., unresolved)
21+
22+
Ref (for individual source component):
23+
Sokoloff et al. (1998) Eq 2
24+
O'Sullivan et al. (2012) Eq 8
25+
Ma et al. (2019a) Eq 10
26+
27+
"""
28+
29+
# Calculate the complex fractional q and u spectra
30+
# fmt: off
31+
pArr1 = pDict["fracPol1"] * np.ones_like(lamSqArr_m2)
32+
pArr2 = pDict["fracPol2"] * np.ones_like(lamSqArr_m2)
33+
pArr3 = pDict["fracPol3"] * np.ones_like(lamSqArr_m2)
34+
quArr1 = pArr1 * np.exp( 2j * (np.radians(pDict["psi01_deg"]) +
35+
pDict["RM1_radm2"] * lamSqArr_m2))
36+
quArr2 = pArr2 * np.exp( 2j * (np.radians(pDict["psi02_deg"]) +
37+
pDict["RM2_radm2"] * lamSqArr_m2))
38+
quArr3 = pArr3 * np.exp( 2j * (np.radians(pDict["psi03_deg"]) +
39+
pDict["RM3_radm2"] * lamSqArr_m2))
40+
quArr = (quArr1 + quArr2 + quArr3)
41+
# fmt: on
42+
43+
return quArr
44+
45+
46+
# -----------------------------------------------------------------------------#
47+
# Priors for the above model. #
48+
# See https://lscsoft.docs.ligo.org/bilby/prior.html for details. #
49+
# #
50+
# -----------------------------------------------------------------------------#
51+
def converter(parameters):
52+
"""
53+
Function to convert between sampled parameters and constraint parameter.
54+
55+
Parameters
56+
----------
57+
parameters: dict
58+
Dictionary containing sampled parameter values, 'RM1_radm2', 'RM2_radm2',
59+
'RM3_radm2', 'fracPol1', 'fracPol2', 'fracPol3'
60+
61+
Returns
62+
-------
63+
dict: Dictionary with constraint parameter 'delta_RM1_RM2_radm2' and 'sum_p1_p2_p3' added.
64+
"""
65+
converted_parameters = parameters.copy()
66+
converted_parameters["delta_RM1_RM2_radm2"] = (
67+
parameters["RM1_radm2"] - parameters["RM2_radm2"]
68+
)
69+
converted_parameters["delta_RM2_RM3_radm2"] = (
70+
parameters["RM2_radm2"] - parameters["RM3_radm2"]
71+
)
72+
converted_parameters["sum_p1_p2_p3"] = (
73+
parameters["fracPol1"] + parameters["fracPol2"] + parameters["fracPol3"]
74+
)
75+
return converted_parameters
76+
77+
78+
priors = PriorDict(conversion_function=converter)
79+
80+
priors["fracPol1"] = bilby.prior.Uniform(
81+
minimum=0.0,
82+
maximum=1.0,
83+
name="fracPol1",
84+
latex_label="$p_1$",
85+
)
86+
priors["fracPol2"] = bilby.prior.Uniform(
87+
minimum=0.0,
88+
maximum=1.0,
89+
name="fracPol2",
90+
latex_label="$p_2$",
91+
)
92+
priors["fracPol3"] = bilby.prior.Uniform(
93+
minimum=0.0,
94+
maximum=1.0,
95+
name="fracPol3",
96+
latex_label="$p_3$",
97+
)
98+
99+
priors["psi01_deg"] = bilby.prior.Uniform(
100+
minimum=0,
101+
maximum=180.0,
102+
name="psi01_deg",
103+
latex_label="$\psi_{0,1}$ (deg)",
104+
boundary="periodic",
105+
)
106+
priors["psi02_deg"] = bilby.prior.Uniform(
107+
minimum=0,
108+
maximum=180.0,
109+
name="psi02_deg",
110+
latex_label="$\psi_{0,2}$ (deg)",
111+
boundary="periodic",
112+
)
113+
priors["psi03_deg"] = bilby.prior.Uniform(
114+
minimum=0,
115+
maximum=180.0,
116+
name="psi03_deg",
117+
latex_label="$\psi_{0,3}$ (deg)",
118+
boundary="periodic",
119+
)
120+
121+
priors["RM1_radm2"] = bilby.prior.Uniform(
122+
minimum=-1100.0,
123+
maximum=1100.0,
124+
name="RM1_radm2",
125+
latex_label="$\phi_1$ (rad m$^{-2}$)",
126+
)
127+
priors["RM2_radm2"] = bilby.prior.Uniform(
128+
minimum=-1100.0,
129+
maximum=1100.0,
130+
name="RM2_radm2",
131+
latex_label="$\phi_2$ (rad m$^{-2}$)",
132+
)
133+
priors["RM3_radm2"] = bilby.prior.Uniform(
134+
minimum=-1100.0,
135+
maximum=1100.0,
136+
name="RM3_radm2",
137+
latex_label="$\phi_3$ (rad m$^{-2}$)",
138+
)
139+
priors["delta_RM1_RM2_radm2"] = Constraint(
140+
minimum=0,
141+
maximum=2200.0,
142+
name="delta_RM1_RM2_radm2",
143+
latex_label="$\Delta\phi_{1,2}$ (rad m$^{-2}$)",
144+
)
145+
priors["delta_RM2_RM3_radm2"] = Constraint(
146+
minimum=0,
147+
maximum=2200.0,
148+
name="delta_RM2_RM3_radm2",
149+
latex_label="$\Delta\phi_{1,2}$ (rad m$^{-2}$)",
150+
)
151+
priors["sum_p1_p2_p3"] = Constraint(
152+
minimum=0.0,
153+
maximum=1,
154+
name="sum_p1_p2_p3",
155+
latex_label="$p_1+p_2+p_3$)",
156+
)

‎RMtools_1D/models_ns/m2.py‎

Lines changed: 15 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -13,15 +13,27 @@
1313
# quArr = Complex array containing the Re and Im spectra. #
1414
# -----------------------------------------------------------------------------#
1515
def model(pDict, lamSqArr_m2):
16-
"""Single Faraday component with Burn depolarisation"""
16+
"""
17+
18+
Single Faraday component with external Faraday dispersion
19+
20+
Ref:
21+
Burn (1966) Eq 21
22+
Sokoloff et al. (1998) Eq B3
23+
O'Sullivan et al. (2012) Eq 11
24+
Ma et al. (2019a) Eq 13
25+
26+
"""
1727

1828
# Calculate the complex fractional q and u spectra
29+
# fmt: off
1930
pArr = pDict["fracPol"] * np.ones_like(lamSqArr_m2)
2031
quArr = (
2132
pArr
2233
* np.exp(2j * (np.radians(pDict["psi0_deg"]) + pDict["RM_radm2"] * lamSqArr_m2))
2334
* np.exp(-2.0 * pDict["sigmaRM_radm2"] ** 2.0 * lamSqArr_m2**2.0)
2435
)
36+
# fmt: on
2537

2638
return quArr
2739

@@ -33,7 +45,7 @@ def model(pDict, lamSqArr_m2):
3345
# -----------------------------------------------------------------------------#
3446
priors = {
3547
"fracPol": bilby.prior.Uniform(
36-
minimum=0.001, maximum=1.0, name="fracPol", latex_label=r"$p$"
48+
minimum=0.0, maximum=1.0, name="fracPol", latex_label=r"$p$"
3749
),
3850
"psi0_deg": bilby.prior.Uniform(
3951
minimum=0,
@@ -50,7 +62,7 @@ def model(pDict, lamSqArr_m2):
5062
),
5163
"sigmaRM_radm2": bilby.prior.Uniform(
5264
minimum=0.0,
53-
maximum=1000.0,
65+
maximum=100.0,
5466
name="sigmaRM_radm2",
5567
latex_label=r"$\sigma_{RM}$ (rad m$^{-2}$)",
5668
),

‎RMtools_1D/models_ns/m3.py‎

Lines changed: 21 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -14,12 +14,24 @@
1414
# quArr = Complex array containing the Re and Im spectra. #
1515
# -----------------------------------------------------------------------------#
1616
def model(pDict, lamSqArr_m2):
17-
"""Two separate Faraday components, averaged within same telescope beam
18-
(i.e., unresolved), with a common Burn depolarisation term."""
17+
"""
18+
19+
Two separate Faraday components with external Faraday dispersion
20+
With a common depolarisation term
21+
Averaged within the same telescope beam (i.e., unresolved)
22+
23+
Ref (for individual source component):
24+
Burn (1966) Eq 21
25+
Sokoloff et al. (1998) Eq B3
26+
O'Sullivan et al. (2012) Eq 11
27+
Ma et al. (2019a) Eq 13
28+
29+
"""
1930

2031
# Calculate the complex fractional q and u spectra
2132
pArr1 = pDict["fracPol1"] * np.ones_like(lamSqArr_m2)
2233
pArr2 = pDict["fracPol2"] * np.ones_like(lamSqArr_m2)
34+
# fmt: off
2335
quArr1 = pArr1 * np.exp(
2436
2j * (np.radians(pDict["psi01_deg"]) + pDict["RM1_radm2"] * lamSqArr_m2)
2537
)
@@ -29,7 +41,7 @@ def model(pDict, lamSqArr_m2):
2941
quArr = (quArr1 + quArr2) * np.exp(
3042
-2.0 * pDict["sigmaRM_radm2"] ** 2.0 * lamSqArr_m2**2.0
3143
)
32-
44+
# fmt: on
3345
return quArr
3446

3547

@@ -62,13 +74,14 @@ def converter(parameters):
6274
priors = PriorDict(conversion_function=converter)
6375

6476
priors["fracPol1"] = bilby.prior.Uniform(
65-
minimum=0.001,
77+
minimum=0.0,
6678
maximum=1.0,
6779
name="fracPol1",
6880
latex_label=r"$p_1$",
6981
)
82+
7083
priors["fracPol2"] = bilby.prior.Uniform(
71-
minimum=0.001,
84+
minimum=0.0,
7285
maximum=1.0,
7386
name="fracPol2",
7487
latex_label=r"$p_2$",
@@ -114,9 +127,10 @@ def converter(parameters):
114127
name="sigmaRM_radm2",
115128
latex_label=r"$\sigma_{RM}$ (rad m$^{-2}$)",
116129
)
130+
117131
priors["sum_p1_p2"] = Constraint(
118-
minimum=0.001,
132+
minimum=0.0,
119133
maximum=1,
120134
name="sum_p1_p2",
121-
latex_label=r"$p_1+p_2$",
135+
latex_label=r"$p_1+p_2$)",
122136
)

0 commit comments

Comments
 (0)