diff --git a/RMtools_1D/do_QUfit_1D_mnest.py b/RMtools_1D/do_QUfit_1D_mnest.py index 0987ab2..d7d4e5a 100755 --- a/RMtools_1D/do_QUfit_1D_mnest.py +++ b/RMtools_1D/do_QUfit_1D_mnest.py @@ -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: @@ -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"]), diff --git a/RMtools_1D/do_RMsynth_1D.py b/RMtools_1D/do_RMsynth_1D.py index be198f9..4c8b3d2 100755 --- a/RMtools_1D/do_RMsynth_1D.py +++ b/RMtools_1D/do_RMsynth_1D.py @@ -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) @@ -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!") diff --git a/RMutils/util_misc.py b/RMutils/util_misc.py index 186ec4a..29b10d7 100644 --- a/RMutils/util_misc.py +++ b/RMutils/util_misc.py @@ -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 @@ -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'] diff --git a/VERSION_HISTORY.txt b/VERSION_HISTORY.txt index 176f9f1..935b7f6 100644 --- a/VERSION_HISTORY.txt +++ b/VERSION_HISTORY.txt @@ -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 diff --git a/tests/RMsynth1D_referencevalues.json b/tests/RMsynth1D_referencevalues.json index 354eb55..d9775d0 100644 --- a/tests/RMsynth1D_referencevalues.json +++ b/tests/RMsynth1D_referencevalues.json @@ -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} \ No newline at end of file +{"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}