Retrieval of a few parameters with Nested Sampling
This example uses nested sampling. In archNEMESIS, nested sampling is done through the pyMultiNest package, which itself is a wrapper around the MultiNest code. Detailed instructions on installation can be found here: http://johannesbuchner.github.io/PyMultiNest/install.html
[1]:
from archnemesis import *
Let’s load in an example Neptune setup.
[2]:
runname = 'neptune'
#Reading the input files
Atmosphere,Measurement,Spectroscopy,Scatter,Stellar,Surface,CIA,Layer,Variables,Retrieval = read_input_files(runname)
[3]:
Atmosphere.plot_Atm()
Atmosphere.plot_Dust()
Atmosphere.DUST_UNITS_FLAG = [-1]*Atmosphere.DUST.shape[1]
We’ll generate a spectrum, add some noise to it, and save it to our spx file. We’re only using a small region of the spectrum for speed.
[4]:
%%time
ForwardModel = ForwardModel_0(runname=runname, Atmosphere=Atmosphere,Surface=Surface,Measurement=Measurement,Spectroscopy=Spectroscopy,Stellar=Stellar,Scatter=Scatter,CIA=CIA,Layer=Layer,Variables=Variables)
SPECONV = ForwardModel.nemesisfm()
Normalisation of phase function should be 1.0
Minimum integral of phase function is 1.0024492587068203
Maximum integral of phase function is 1.0029856636186842
CIRSrad :: Aerosol optical depths at 1.8490000175515888 :: [0.37875526 0. ]
CIRSrad :: Performing multiple scattering calculation
CIRSrad :: NF = 6 ; NMU = 5 ; NPHI = 100
CPU times: user 16.4 s, sys: 431 ms, total: 16.8 s
Wall time: 16.7 s
[5]:
new_spectrum = SPECONV * np.random.normal(1,0.02,SPECONV.shape)
new_spectrum_error = new_spectrum*0.02
Measurement.MEAS[:,0] = new_spectrum
Measurement.ERRMEAS[:,0] = new_spectrum_error
#Not actually saving, so saved nested sampling results are still valid
# Measurement.write_spx()
Now, let’s reload, with our new slightly noisy spectrum. We’re keeping the apr file the same. This means that if we were running optimal estimation, we would start at the optimal point and the retrieval would end quickly. However, because nested sampling draws randomly from the prior distribution, it will still fully explore the space defined in the apr file.
[6]:
Atmosphere,Measurement,Spectroscopy,Scatter,Stellar,Surface,CIA,Layer,Variables,Retrieval = read_input_files(runname)
[19]:
Measurement.plot_nadir()
Now let’s run our nested sampling. This would usually take a long time (~40 minutes for me), but this example folder already contains the results from a previous run, so we can just investigate the output.
[8]:
import subprocess
legacy_files=True
retrieval_method=1 #Nested Sampling
NCores=30
subprocess.run(["mpiexec", "-np", str(NCores), "python", "-c",
f"import archnemesis as ans; ans.Retrievals.retrieval_nemesis('{runname}', legacy_files={legacy_files}, retrieval_method={retrieval_method}, NCores={NCores})"])
*****************************************************
MultiNest v3.10
Copyright Farhan Feroz & Mike Hobson
Release Jul 2015
no. of live points = 400
dimensionality = 3
resuming from previous job
*****************************************************
Starting MultiNest
Acceptance Rate: 0.505920
Replacements: 3931
Total Samples: 7770
Nested Sampling ln(Z): -119.706506
Importance Nested Sampling ln(Z): -119.794794 +/- 0.073413
ln(ev)= -119.41738389222705 +/- 0.12783258525578428
Total Likelihood Evaluations: 7770
Sampling finished. Exiting MultiNest
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.5007145404815674
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.490624189376831
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.477888345718384
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.5092155933380127
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.503382921218872
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.5163533687591553
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.5170300006866455
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.53873872756958
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.5503487586975098
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.5413360595703125
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.5618622303009033
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.5623037815093994
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.5251457691192627
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.5424697399139404
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.5078117847442627
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.528520107269287
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.544609785079956
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.5556983947753906
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.517728805541992
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.542154312133789
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.522583484649658
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.5330421924591064
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.50384783744812
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.5602753162384033
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.540872573852539
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.550699472427368
analysing data from chains/.txt
Evidence: -119.4 +- 0.1
Parameter values:
0 : -0.904 +- 0.014
1 : -1.886 +- 0.007
2 : 0.177 +- 0.031
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.6063857078552246
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.5783612728118896
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.6308035850524902
analysing data from chains/.txt
analysing data from chains/.txt
Creating marginal plot ...
Model run OK
Elapsed time (s) = 2.6241729259490967
[8]:
CompletedProcess(args=['mpiexec', '-np', '30', 'python', '-c', "import archnemesis as ans; ans.Retrievals.retrieval_nemesis('neptune', legacy_files=True, retrieval_method=1, NCores=30)"], returncode=0)
Let’s also do an optimal estimation run and compare. This should be quick.
[9]:
legacy_files=True
retrieval_method=0
NCores=4
ans.Retrievals.retrieval_nemesis(runname,legacy_files=legacy_files,retrieval_method=retrieval_method,NCores=NCores)
nemesis :: Calculating Jacobian matrix KK
Calculating numerical part of the Jacobian :: running 4 forward models
Calculating forward model 1/4
Calculating forward model 3/4
Calculating forward model 4/4
Calculating forward model 2/4
Calculated forward model 1/4
Calculated forward model 4/4
Calculated forward model 3/4
Calculated forward model 2/4
nemesis :: Calculating gain matrix
nemesis :: Calculating cost function
calc_phiret: phi1, phi2 = 328.08172199802243, 0.0)
chisq/ny = 1.2967656995969266
Assess:
Average of diagonal elements of Kk*Sx*Kt : 9.594884787809242e-18
Average of diagonal elements of Se : 4.686677258343464e-19
Ratio = 20.47268087583359
Average of Kk*Sx*Kt/Se element ratio : 25.624506543596524
******************* ASSESS WARNING *****************
Insufficient constraint. Solution likely to be exact
****************************************************
nemesis :: Iteration 0/10
nemesis :: Calculating next iterated state vector
Normalisation of phase function should be 1.0
Minimum integral of phase function is 1.0024492587068203
Maximum integral of phase function is 1.0029856636186842
nemesis :: Calculating Jacobian matrix KK
Calculating numerical part of the Jacobian :: running 4 forward models
Calculating forward model 1/4
Calculating forward model 2/4
Calculating forward model 3/4
Calculating forward model 4/4
Calculated forward model 1/4
Calculated forward model 2/4
Calculated forward model 3/4
Calculated forward model 4/4
calc_phiret: phi1, phi2 = 250.63293345692614, 0.04432732195668598)
chisq/ny = 0.9906440057586013
Successful iteration. Updating xn,yn and kk
calc_phiret: phi1, phi2 = 250.63293345692614, 0.04432732195668598)
nemesis :: Iteration 1/10
nemesis :: Calculating next iterated state vector
Normalisation of phase function should be 1.0
Minimum integral of phase function is 1.0024492587068203
Maximum integral of phase function is 1.0029856636186842
nemesis :: Calculating Jacobian matrix KK
Calculating numerical part of the Jacobian :: running 4 forward models
Calculating forward model 1/4
Calculating forward model 2/4
Calculating forward model 3/4
Calculating forward model 4/4
Calculated forward model 2/4
Calculated forward model 1/4
Calculated forward model 3/4
Calculated forward model 4/4
calc_phiret: phi1, phi2 = 224.05712639068108, 0.13206217479772497)
chisq/ny = 0.8856012900817434
Successful iteration. Updating xn,yn and kk
calc_phiret: phi1, phi2 = 224.05712639068108, 0.13206217479772497)
nemesis :: Iteration 2/10
nemesis :: Calculating next iterated state vector
Normalisation of phase function should be 1.0
Minimum integral of phase function is 1.0024492587068203
Maximum integral of phase function is 1.0029856636186842
nemesis :: Calculating Jacobian matrix KK
Calculating numerical part of the Jacobian :: running 4 forward models
Calculating forward model 1/4
Calculating forward model 2/4
Calculating forward model 3/4
Calculating forward model 4/4
Calculated forward model 1/4
Calculated forward model 2/4
Calculated forward model 3/4
Calculated forward model 4/4
calc_phiret: phi1, phi2 = 223.2224322890125, 0.1649909942286262)
chisq/ny = 0.8823021039091403
Successful iteration. Updating xn,yn and kk
calc_phiret: phi1, phi2 = 223.2224322890125, 0.1649909942286262)
nemesis :: Iteration 3/10
nemesis :: Calculating next iterated state vector
Normalisation of phase function should be 1.0
Minimum integral of phase function is 1.0024492587068203
Maximum integral of phase function is 1.0029856636186842
nemesis :: Calculating Jacobian matrix KK
Calculating numerical part of the Jacobian :: running 4 forward models
Calculating forward model 1/4
Calculating forward model 2/4
Calculating forward model 3/4
Calculating forward model 4/4
Calculated forward model 4/4
Calculated forward model 3/4
Calculated forward model 2/4
Calculated forward model 1/4
calc_phiret: phi1, phi2 = 223.1923017069037, 0.17103934096698375)
chisq/ny = 0.8821830106992241
Successful iteration. Updating xn,yn and kk
calc_phiret: phi1, phi2 = 223.1923017069037, 0.17103934096698375)
nemesis :: Iteration 4/10
nemesis :: Calculating next iterated state vector
Normalisation of phase function should be 1.0
Minimum integral of phase function is 1.0024492587068203
Maximum integral of phase function is 1.0029856636186842
nemesis :: Calculating Jacobian matrix KK
Calculating numerical part of the Jacobian :: running 4 forward models
Calculating forward model 1/4
Calculating forward model 2/4
Calculating forward model 3/4
Calculating forward model 4/4
Calculated forward model 4/4
Calculated forward model 1/4
Calculated forward model 2/4
Calculated forward model 3/4
calc_phiret: phi1, phi2 = 223.18688597840355, 0.17222229570918313)
chisq/ny = 0.8821616046577215
Successful iteration. Updating xn,yn and kk
calc_phiret: phi1, phi2 = 223.18688597840355, 0.17222229570918313)
phi, phlimit : 0.0018950172119009561,0.01
Phi has converged
Terminating retrieval
Model run OK
Elapsed time (s) = 32.4851336479187
Let’s compare these two results.
[10]:
lat,lon,ngeom,ny,wave,specret,specmeas,specerrmeas,nx,Var,aprprof,aprerr,retprof,reterr = ans.Files.read_mre('neptune')
[11]:
from archnemesis.NestedSampling_0 import coreretNS
NestedSampling = coreretNS(runname,Variables,Measurement,Atmosphere,Spectroscopy,Scatter,Stellar,Surface,CIA,Layer)
*****************************************************
MultiNest v3.10
Copyright Farhan Feroz & Mike Hobson
Release Jul 2015
no. of live points = 400
dimensionality = 3
resuming from previous job
*****************************************************
Starting MultiNest
Acceptance Rate: 0.505920
Replacements: 3931
Total Samples: 7770
Nested Sampling ln(Z): -119.706506
Importance Nested Sampling ln(Z): -119.794794 +/- 0.073413
analysing data from chains/.txt
ln(ev)= -119.41738389222705 +/- 0.12783258525578428
Total Likelihood Evaluations: 7770
Sampling finished. Exiting MultiNest
Evidence: -119.4 +- 0.1
Parameter values:
0 : -0.904 +- 0.014
1 : -1.886 +- 0.007
2 : 0.176 +- 0.031
[12]:
NestedSampling.compare()
analysing data from chains/.txt
We can see that both methods converge to the same value (within error), and the posterior distributions look fairly similar. However, this is a very simplified example - often the two methods do not agree so much!