Adding a planet to a params object does not propagate through basis transformations correctly

[sealauren]

Is it possible to add a planet to an existing params object? I made the following script to add an additional planet, but I run into a basis problem further along in which the tp of the new planet does not get calculated.

Code:

def add_planet(params,NT_per=0.0, NT_phase = 0.0, NT_per_vary=True, ecc_vary = True, dvdt_vary=True, curv_vary=False): params.num_planets+=1 npl = params.num_planets # new planet number print("%d planets in this model..." % npl) params['per'+str(npl)] = radvel.Parameter(value=NT_per,vary=NT_per_vary) params['tc'+str(npl)] = radvel.Parameter(value=t.min() + np.sin(NT_phase - np.pi/2.)*NT_per) # use smart phase params['sesinw'+str(npl)] = radvel.Parameter(value=0.,vary=ecc_vary) params['secosw'+str(npl)] = radvel.Parameter(value=0.,vary=ecc_vary) params['k'+str(npl)] = radvel.Parameter(50.) params['dvdt'] = radvel.Parameter(value=-1e-3,vary=dvdt_vary) params['curv'] = radvel.Parameter(value=0.,vary=curv_vary) time_base = np.median(t) mod = radvel.RVModel(params, time_base=time_base) return mod

Error: Fitting Keplerian for new planet Period: 1914.95263667 Phase: 1.92436073184 3 planets in this model... Initial guess for N_giants=2: parameter value vary per1 0.66931 False tc1 134.132 False secosw1 0 False sesinw1 0 False k1 -0.890949 True per2 2112.42 True tc2 -932.821 True secosw2 0.408161 True sesinw2 0.521555 True k2 170.161 True dvdt -0.001 True curv 0 False gamma 0.1 True jit 1 True per3 1914.95 True tc3 1527.06 True sesinw3 0 True secosw3 0 True k3 50 True

Traceback (most recent call last): File "KGPS_generalized.py", line 322, in method='Nelder-Mead', # Nelder-Mead also works File "/Users/lweiss/anaconda/lib/python2.7/site-packages/scipy/optimize/_minimize.py", line 475, in minimize return _minimize_neldermead(fun, x0, args, callback, *options) File "/Users/lweiss/anaconda/lib/python2.7/site-packages/scipy/optimize/optimize.py", line 532, in _minimize_neldermead fsim[k] = func(sim[k]) File "/Users/lweiss/anaconda/lib/python2.7/site-packages/scipy/optimize/optimize.py", line 292, in function_wrapper return function((wrapper_args + args)) File "/Users/lweiss/anaconda/lib/python2.7/site-packages/radvel-1.1.5-py2.7-macosx-10.6-x86_64.egg/radvel/likelihood.py", line 96, in neglogprob_array return -self.logprob_array(params_array) File "/Users/lweiss/anaconda/lib/python2.7/site-packages/radvel-1.1.5-py2.7-macosx-10.6-x86_64.egg/radvel/posterior.py", line 79, in logprob_array _logprob = self.logprob() File "/Users/lweiss/anaconda/lib/python2.7/site-packages/radvel-1.1.5-py2.7-macosx-10.6-x86_64.egg/radvel/posterior.py", line 47, in logprob _logprob = self.likelihood.logprob() File "/Users/lweiss/anaconda/lib/python2.7/site-packages/radvel-1.1.5-py2.7-macosx-10.6-x86_64.egg/radvel/likelihood.py", line 264, in logprob residuals = self.residuals() File "/Users/lweiss/anaconda/lib/python2.7/site-packages/radvel-1.1.5-py2.7-macosx-10.6-x86_64.egg/radvel/likelihood.py", line 227, in residuals res = self.y - self.params[self.gamma_param].value - self.model(self.x) File "/Users/lweiss/anaconda/lib/python2.7/site-packages/radvel-1.1.5-py2.7-macosx-10.6-x86_64.egg/radvel/model.py", line 214, in call tp = params_synth['tp{}'.format(num_planet)].value KeyError: 'tp3'

It looks like tp3 never gets calculated for the new planet, even though it is in the same basis as the other planets. Is this a radvel bug, or is there just a better way to do this?

[sblunt]

I think you just need to convert from the input basis (you input the new planet parameters in 'tc') to the fitting basis ('tp'). Does that fix it?

[petigura]

I believe the basis object also has a nplanets attribute (which tells the kepler synthesizer how many planets to loop over). Therefore, I don't recommend changing the number of planets on the fly. The way we do this in the BIC comparison is by fixing setting K to zero for the "unused" planets and fixing all the other parameters.

[bjfultn]

It looks like @sealauren was trying to update the params.num_planet parameter so I'm not exactly sure why that doesn't work. I think its probably a Python copying vs. pointer issue. I could potentially look into addressing this, but its going to be much safer to initialize a new Parameters object and push the old parameter values in to the new object. That way you know you are dealing with a new parameter set and you aren't going to be fooled by Python pointers looking like new objects when they really aren't. Is that a viable solution @sealauren?

[bjfultn]

Closing since a solution was found