dtamayo / spock

🖖 A package to determine whether planetary orbital configurations will live long and prosper
GNU General Public License v3.0
63 stars 15 forks source link

Add regression model to spock #7

Closed MilesCranmer closed 4 years ago

MilesCranmer commented 4 years ago

This can be used the same as StabilityClassifier, even with more than 3 planets, e.g.,:

import sys
sys.path.append('/mnt/home/mcranmer/spock')
from spock import StabilityRegression
import rebound
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
%matplotlib inline

sim = rebound.SimulationArchive(
'/mnt/ceph/users/mcranmer/general_data/orbital_integrations/full_resonant_p1_000003000000/simulation_archives/runs/satest000003121011.bin')[0]

outer = sim.particles[3]

model = StabilityRegression()
out = model.predict(sim, samples=100000)
plt.hist(np.log10(out).ravel(), alpha=0.3, bins=50, label='normal', density=True, range=[4, 12])
simt = sim.copy()
simt.add(m=outer.m, a=outer.a*1.01)
out = model.predict(simt, samples=100000)
plt.hist(np.log10(out).ravel(), alpha=0.3, bins=50, label='extra particle', density=True, range=[4, 12])

plt.hist(np.log10(np.array([2.39455e+06, 2.03797e+06])), bins=50, label='truth', density=True, range=[4, 12])

plt.legend()
plt.show()

Will output this:

Screen Shot 2020-02-24 at 12 26 51 AM
MilesCranmer commented 4 years ago

Pausing merge until later.