Skip to content

Added Stokes I model parameter errors as standard output for rmsynth1… #2

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
AlecThomson merged 1 commit into from
Apr 30, 2022
Merged

Added Stokes I model parameter errors as standard output for rmsynth1… #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion RMtools_1D/do_QUfit_1D_mnest.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -253,6 +253,8 @@ def run_qufit(dataFile, modelNum, outDir="", polyOrd=2, nBits=32,
MPI.Finalize()
return



# If no Stokes I present, create a dummy spectrum = unity
if noStokesI:
if mpiRank==0:
Expand Down Expand Up @@ -471,7 +473,8 @@ def run_qufit(dataFile, modelNum, outDir="", polyOrd=2, nBits=32,
"ln(EVIDENCE) ": toscalar(lnEvidence),
"dLn(EVIDENCE)": toscalar(dLnEvidence),
"nFree": toscalar(nFree),
"Imodel": toscalar(IfitDict["p"]),
"Imodel": toscalar(",".join([str(x) for x in IfitDict["p"]])),
"Imodel_errs": toscalar(",".join([str(x) for x in IfitDict["perror"]])),
"IfitChiSq": toscalar(IfitDict["chiSq"]),
"IfitChiSqRed": toscalar(IfitDict["chiSqRed"]),
"IfitPolyOrd": toscalar(IfitDict["polyOrd"]),
Expand Down
11 changes: 7 additions & 4 deletions RMtools_1D/do_RMsynth_1D.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -314,9 +314,9 @@ def run_rmsynth(data, polyOrd=2, phiMax_radm2=None, dPhi_radm2=None,
Ifreq0 = modStokesI_interp(freq0_Hz)
dirtyFDF *= (Ifreq0) # FDF is in fracpol units initially, convert back to flux

if modStokesI is None:
#Need to renormalize the Stokes I parameters here to the actual reference frequency.
fitDict=renormalize_StokesI_model(fitDict,freq0_Hz)
# if modStokesI is None:
# #Need to renormalize the Stokes I parameters here to the actual reference frequency.
# fitDict=renormalize_StokesI_model(fitDict,freq0_Hz)

# Calculate the theoretical noise in the FDF !!Old formula only works for wariance weights!
weightArr = np.where(np.isnan(weightArr), 0.0, weightArr)
Expand All @@ -333,15 +333,18 @@ def run_rmsynth(data, polyOrd=2, phiMax_radm2=None, dPhi_radm2=None,
lam0Sq = lam0Sq_m2)
mDict["Ifreq0"] = toscalar(Ifreq0)
mDict["polyCoeffs"] = ",".join([str(x) for x in fitDict["p"]])
mDict["polyCoefferr"] = ",".join([str(x) for x in fitDict["perror"]])
mDict["poly_reffreq"] = fitDict['reference_frequency_Hz']
mDict['polyOrd'] = fitDict['polyOrd']
mDict["IfitStat"] = fitDict["fitStatus"]
mDict["IfitChiSqRed"] = fitDict["chiSqRed"]
mDict["fit_function"] = fit_function
mDict["lam0Sq_m2"] = toscalar(lam0Sq_m2)
mDict["freq0_Hz"] = toscalar(freq0_Hz)
mDict["fwhmRMSF"] = toscalar(fwhmRMSF)
mDict["dQU"] = toscalar(nanmedian(dQUArr))
mDict["dFDFth"] = toscalar(dFDFth)
mDict["units"] = units
mDict['polyOrd'] = fitDict['polyOrd']

if (fitDict["fitStatus"] >= 128) and verbose:
log("WARNING: Stokes I model contains negative values!")
Expand Down
8 changes: 7 additions & 1 deletion RMutils/util_misc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -385,6 +385,7 @@ def fit_StokesI_model(freqArr,IArr,dIArr,polyOrd,fit_function="log"):
fitDict["chiSq"] = mp.fnorm
fitDict["chiSqRed"] = mp.fnorm/fitDict["dof"]
fitDict["nIter"] = mp.niter
fitDict["perror"] = mp.perror


return fitDict
Expand All @@ -411,7 +412,12 @@ def renormalize_StokesI_model(fitDict,new_reference_frequency):
and fixes the fit parameters such that the the model is the same. This is
important because the initial Stokes I fitted model uses an arbitrary
reference frequency, and it may be desirable for users to know the exact
reference frequency of the model."""
reference frequency of the model.

This function is depreciated because it can't propagate the errors in
the fit parameters."""
print("The renormalize_StokesI_model function is depreciated because it can't propagate errors.\n"
"If this message appears, it's been invoked (perhaps by legacy code?)")
#Renormalization ratio:
x=new_reference_frequency/fitDict['reference_frequency_Hz']
(a,b,c,d,f,g)=fitDict['p']
Expand Down
6 changes: 6 additions & 0 deletions VERSION_HISTORY.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,12 @@ Github main branch:
This is because the fitter was found to be vulnerable to numerical issues
when fitting log-polynomials in 32-bits (thanks to Shannon Vanderwoude for
discovering the problem).
-added parameter errors of Stokes I fit as output to Stokes I fitDict, to
RMsynth1D outputs, and QU-fitting outputs. These are the 1-sigma errors
computed by the mpfit routine.
-the renormalize_StokesI_model function is now depreciated, since it cannot
transform the errors. It is left in for the time being in case it is needed
for any legacy reasons.

1.1.2:
-modified the generation of maxPI and peakRM maps in RMsynth3D so that
Expand Down
2 changes: 1 addition & 1 deletion tests/RMsynth1D_referencevalues.json
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1 +1 @@
{"dFDFcorMAD": 0.012824177742004395, "dFDFrms": 0.03170722723007202, "phiPeakPIchan_rm2": 201.0, "dPhiPeakPIchan_rm2": 0.24814919763068932, "ampPeakPIchan": 0.6994333267211914, "ampPeakPIchanEff": 0.6993762345202659, "dAmpPeakPIchan": 0.005892556709695188, "snrPIchan": 118.69776757012694, "indxPeakPIchan": 267, "peakFDFimagChan": -0.5605906844139099, "peakFDFrealChan": 0.41826435923576355, "polAngleChan_deg": 153.36356353759766, "dPolAngleChan_deg": 0.24135154639379308, "polAngle0Chan_deg": 44.00376296895752, "dPolAngle0Chan_deg": 1.3731897339320969, "phiPeakPIfit_rm2": 200.29003676576465, "dPhiPeakPIfit_rm2": 0.248062934269283, "ampPeakPIfit": 0.6996765531831387, "ampPeakPIfitEff": 0.699619480830623, "dAmpPeakPIfit": 0.005892556709695188, "snrPIfit": 118.73904446807298, "indxPeakPIfit": 266.7633455885882, "peakFDFimagFit": -0.5593282099678656, "peakFDFrealFit": 0.4190415975286484, "polAngleFit_deg": 153.42004263949445, "dPolAngleFit_deg": 0.24126764607950102, "polAngle0Fit_deg": 48.26124570608431, "dPolAngle0Fit_deg": 1.372712376102901, "Ifreq0": 1.0, "polyCoeffs": "0.0,0.0,0.0,0.0,0.0,1.0", "IfitStat": 4, "IfitChiSqRed": 0.0, "lam0Sq_m2": 0.10327484831236765, "freq0_Hz": 932874912.6426204, "fwhmRMSF": 58.90951156616211, "dQU": 0.10000000149011612, "dFDFth": 0.005892556709695188, "units": "Jy/beam", "polyOrd": 2, "min_freq": 800000000.0, "max_freq": 1088000000.0, "N_channels": 288, "median_channel_width": 1003456.0, "fracPol": 0.6996765531831387, "sigmaAddQ": 0.39312059810649663, "dSigmaAddMinusQ": 0.2859949907172898, "dSigmaAddPlusQ": 0.13645904926608626, "sigmaAddU": 0.22680054259742569, "dSigmaAddMinusU": 0.20949437831506873, "dSigmaAddPlusU": 0.20581573041436466}
{"dFDFcorMAD": 0.012824177742004395, "dFDFrms": 0.03170722723007202, "phiPeakPIchan_rm2": 201.0, "dPhiPeakPIchan_rm2": 0.24814919763068932, "ampPeakPIchan": 0.6994333267211914, "ampPeakPIchanEff": 0.6993762345202659, "dAmpPeakPIchan": 0.005892556709695188, "snrPIchan": 118.69776757012694, "indxPeakPIchan": 267, "peakFDFimagChan": -0.5605906844139099, "peakFDFrealChan": 0.41826435923576355, "polAngleChan_deg": 153.36356353759766, "dPolAngleChan_deg": 0.24135154639379308, "polAngle0Chan_deg": 44.00376296895752, "dPolAngle0Chan_deg": 1.3731897339320969, "phiPeakPIfit_rm2": 200.29003676576465, "dPhiPeakPIfit_rm2": 0.248062934269283, "ampPeakPIfit": 0.6996765531831387, "ampPeakPIfitEff": 0.699619480830623, "dAmpPeakPIfit": 0.005892556709695188, "snrPIfit": 118.73904446807298, "indxPeakPIfit": 266.7633455885882, "peakFDFimagFit": -0.5593282099678656, "peakFDFrealFit": 0.4190415975286484, "polAngleFit_deg": 153.42004263949445, "dPolAngleFit_deg": 0.24126764607950102, "polAngle0Fit_deg": 48.26124570608431, "dPolAngle0Fit_deg": 1.372712376102901, "Ifreq0": 1.0, "polyCoeffs": "0.0,0.0,0.0,0.0,0.0,1.0", "polyCoefferr": "0.0,0.0,0.0,19.15748709659049,0.6786377104199252,0.08796823965028172", "poly_reffreq": 944000000.2222222, "IfitStat": 4, "IfitChiSqRed": 0.0, "fit_function": "log", "lam0Sq_m2": 0.10327484831236765, "freq0_Hz": 932874912.6426204, "fwhmRMSF": 58.90951156616211, "dQU": 0.10000000149011612, "dFDFth": 0.005892556709695188, "units": "Jy/beam", "polyOrd": 2, "min_freq": 800000000.0, "max_freq": 1088000000.0, "N_channels": 288, "median_channel_width": 1003456.0, "fracPol": 0.6996765531831387, "sigmaAddQ": 0.39312059810649663, "dSigmaAddMinusQ": 0.2859949907172898, "dSigmaAddPlusQ": 0.13645904926608626, "sigmaAddU": 0.22680054259742569, "dSigmaAddMinusU": 0.20949437831506873, "dSigmaAddPlusU": 0.20581573041436466}