diff --git a/orbitize/example_data/H024332.d b/orbitize/example_data/H024332.d new file mode 100644 index 00000000..d7d1682f --- /dev/null +++ b/orbitize/example_data/H024332.d @@ -0,0 +1,94 @@ +# This file contains residual records, extracted from the Hipparcos 2 +# Interactive Data Access Tool (2014). For more information, see: +# https://www.cosmos.esa.int/web/hipparcos/interactive-data-access +# https://www.cosmos.esa.int/web/hipparcos/catalogues +# +# HIP MCE NRES NC isol_n SCE F2 F1 +# 24332 24273 81 1 5 0 1.17 0 +# Hp B-V VarAnn NOB NR +# 6.5960 0.542 0 81 0 +# RAdeg DEdeg Plx pm_RA pm_DE e_RA e_DE e_Plx e_pmRA e_pmDE dpmRA dpmDE e_dpmRA e_dpmDE ddpmRA ddpmDE e_ddpmRA e_ddpmDE upsRA upsDE e_upsRA e_upsDE var +# 78.32313339 37.33763724 38.29 -145.00 -135.14 0.51 0.26 0.55 0.50 0.27 --- --- --- --- --- --- --- --- --- --- --- --- --- +# +# IORB EPOCH PARF CPSI SPSI RES SRES + 191 -1.1740 0.2376 -0.1869 0.9824 0.63 1.81 + 191 -1.1740 0.2375 -0.1869 0.9824 0.70 2.22 + 192 -1.1727 0.2280 -0.1750 0.9846 1.29 1.62 + 192 -1.1727 0.2289 -0.1753 0.9845 0.24 1.94 + 192 -1.1727 0.2252 -0.1706 0.9853 -1.62 1.60 + 192 -1.1727 0.2241 -0.1690 0.9856 3.60 2.00 + 192 -1.1727 0.2220 -0.1662 0.9861 2.60 1.48 + 216 -1.1435 -0.4636 0.5895 0.8078 -1.74 1.73 + 216 -1.1435 -0.4638 0.5896 0.8077 0.24 1.27 + 216 -1.1435 -0.4680 0.5943 0.8043 1.42 1.65 + 290 -1.0536 0.6767 -0.7471 0.6647 -0.89 1.73 + 290 -1.0536 0.6758 -0.7463 0.6656 -4.63 1.81 + 290 -1.0536 0.6776 -0.7479 0.6638 -3.15 1.70 + 290 -1.0536 0.6771 -0.7474 0.6644 0.82 1.96 + 290 -1.0536 0.6782 -0.7485 0.6632 -0.64 1.25 + 371 -0.9551 -0.3388 0.1595 0.9872 -1.68 1.61 + 371 -0.9551 -0.3387 0.1594 0.9872 -1.85 1.40 + 371 -0.9551 -0.3344 0.1540 0.9881 -0.18 1.43 + 372 -0.9539 -0.3252 0.1405 0.9901 -0.38 1.71 + 372 -0.9539 -0.3242 0.1393 0.9902 0.60 1.49 + 372 -0.9539 -0.3207 0.1343 0.9909 -1.39 1.82 + 372 -0.9539 -0.3187 0.1318 0.9913 -0.15 1.74 + 372 -0.9539 -0.3158 0.1277 0.9918 -1.45 1.96 + 372 -0.9539 -0.3145 0.1262 0.9920 -1.43 1.90 + 372 -0.9539 -0.3092 0.1194 0.9928 -0.59 1.91 + 373 -0.9527 -0.3025 0.1093 0.9940 -0.35 1.97 + 373 -0.9527 -0.3013 0.1078 0.9942 -1.47 1.92 + 373 -0.9527 -0.2972 0.1018 0.9948 0.46 1.80 + 385 -0.9381 -0.0032 -0.3088 0.9511 2.46 2.20 + 385 -0.9381 -0.0017 -0.3108 0.9505 0.89 2.34 + 386 -0.9371 0.0058 -0.3229 0.9464 -1.15 2.06 + 386 -0.9371 0.0110 -0.3300 0.9440 0.99 1.85 + 387 -0.9356 0.0269 -0.3539 0.9353 0.37 1.58 + 387 -0.9356 0.0281 -0.3558 0.9346 0.31 1.46 + 387 -0.9356 0.0368 -0.3662 0.9305 1.17 1.65 + 387 -0.9356 0.0371 -0.3676 0.9300 -1.05 1.42 + 677 -0.5832 -0.6617 -0.6426 -0.7662 -0.41 1.91 + 677 -0.5832 -0.6609 -0.6419 -0.7668 -3.56 2.10 + 745 -0.5005 0.6159 0.7602 -0.6497 -1.31 1.25 + 745 -0.5005 0.6152 0.7594 -0.6506 -1.75 1.26 + 745 -0.5005 0.6140 0.7584 -0.6518 -3.36 1.59 + 745 -0.5005 0.6112 0.7559 -0.6547 3.84 1.80 + 783 -0.4543 -0.4032 -0.2511 -0.9680 -0.80 2.98 + 783 -0.4543 -0.4077 -0.2567 -0.9665 -3.22 2.88 + 783 -0.4543 -0.4071 -0.2562 -0.9666 2.74 2.52 + 1071 -0.1044 0.6198 -0.6600 0.7513 0.71 1.83 + 1071 -0.1044 0.6177 -0.6578 0.7532 0.61 1.35 + 1071 -0.1044 0.6183 -0.6584 0.7526 0.60 1.12 + 1128 -0.0351 -0.6614 0.5774 0.8165 -0.30 1.31 + 1128 -0.0351 -0.6621 0.5780 0.8160 0.85 1.28 + 1128 -0.0351 -0.6598 0.5755 0.8178 -3.85 1.85 + 1128 -0.0351 -0.6600 0.5758 0.8176 2.53 1.76 + 1176 0.0232 0.4723 -0.7162 0.6979 0.17 3.08 + 1176 0.0232 0.4724 -0.7163 0.6977 -4.10 3.56 + 1176 0.0232 0.4754 -0.7194 0.6946 1.57 3.41 + 1176 0.0232 0.4761 -0.7195 0.6944 -4.89 3.12 + 1457 0.3647 -0.4952 -0.6027 -0.7980 1.73 3.32 + 1457 0.3647 -0.4940 -0.6016 -0.7988 -1.52 3.16 + 1457 0.3647 -0.4921 -0.5994 -0.8004 -2.06 3.26 + 1457 0.3647 -0.4910 -0.5985 -0.8011 -0.85 3.08 + 1458 0.3658 -0.4904 -0.5940 -0.8045 0.06 3.15 + 1508 0.4266 0.6755 0.7105 -0.7037 0.57 3.49 + 1854 0.8469 0.4172 -0.4309 0.9024 1.16 2.49 + 1854 0.8469 0.4162 -0.4296 0.9030 8.28 2.66 + 1854 0.8469 0.4114 -0.4238 0.9058 1.80 2.28 + 1855 0.8481 0.4098 -0.4218 0.9067 -2.22 1.61 + 1855 0.8481 0.4096 -0.4220 0.9066 -1.14 1.64 + 1892 0.8930 -0.6099 0.6390 0.7692 -2.99 1.25 + 1892 0.8930 -0.6121 0.6410 0.7675 0.48 1.96 + 1892 0.8930 -0.6108 0.6397 0.7686 -1.69 2.03 + 1957 0.9719 0.6484 -0.7590 0.6511 0.47 1.56 + 1958 0.9732 0.6487 -0.7614 0.6482 -1.34 1.27 + 2245 1.3217 -0.0755 -0.2542 -0.9671 3.35 2.64 + 2246 1.3228 -0.0634 -0.2349 -0.9720 2.06 2.23 + 2246 1.3228 -0.0631 -0.2348 -0.9720 0.27 2.16 + 2246 1.3228 -0.0588 -0.2292 -0.9734 1.09 1.51 + 2246 1.3228 -0.0586 -0.2288 -0.9735 -2.37 1.40 + 2246 1.3228 -0.0540 -0.2228 -0.9749 -1.92 1.32 + 2246 1.3228 -0.0535 -0.2220 -0.9751 -1.69 1.34 + 2717 1.8945 0.6215 -0.6608 0.7505 -1.56 2.36 + 2717 1.8945 0.6208 -0.6601 0.7512 -1.05 1.36 diff --git a/orbitize/example_data/HD33632Ab.csv b/orbitize/example_data/HD33632Ab.csv new file mode 100644 index 00000000..0252e0b3 --- /dev/null +++ b/orbitize/example_data/HD33632Ab.csv @@ -0,0 +1,4 @@ +epoch,object,raoff,raoff_err,decoff,decoff_err,inst +58409.0,1,-761,5,-176,4,CHARIS +58423.0,1,-753,5,-178,5,NIRC2 +59092.0,1,-740,5,-95,3,CHARIS \ No newline at end of file diff --git a/orbitize/example_data/gost_22.4.3_806543_2026-07-28-17-46-26_HD_33632.csv b/orbitize/example_data/gost_22.4.3_806543_2026-07-28-17-46-26_HD_33632.csv new file mode 100644 index 00000000..37502d39 --- /dev/null +++ b/orbitize/example_data/gost_22.4.3_806543_2026-07-28-17-46-26_HD_33632.csv @@ -0,0 +1,76 @@ +Target, ra[rad], dec[rad], ra[h:m:s], dec[d:m:s], ObservationTimeAtGaia[UTC], CcdRow[1-7], zetaFieldAngle[rad], scanAngle[rad], Fov[FovP=preceding/FovF=following], parallaxFactorAlongScan, parallaxFactorAcrossScan, ObservationTimeAtBarycentre[BarycentricJulianDateInTCB] +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2014-09-07T20:02:55.945,1,0.004261441385324073,-2.419338622854937,FoVP,-0.7139807546396543,0.721507425468232,2456908.3357978356 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2014-09-07T21:49:30.161,5,-9.544894992248997E-4,-2.4202204741216615,FoVF,-0.7136305588529797,0.7219838864565272,2456908.409812096 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2014-10-11T19:00:12.343,4,0.001844990104153869,-2.6989244787331725,FoVF,-0.5481129956895733,0.715562776608969,2456942.295425042 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-02-26T06:30:37.544,4,-0.0010163849352416324,-0.7842990960205769,FoVP,0.6685663841239335,0.7038981947222167,2457079.7735618614 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-02-26T08:17:11.748,2,0.0043439215473310136,-0.7831533431631229,FoVF,0.6677295026077388,0.7051080241324862,2457079.8475617887 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-03-28T07:55:49.189,1,0.005457876467591252,0.588823806199045,FoVF,-0.6552996516448154,0.7077481533077975,2457109.8298096643 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-03-28T12:09:29.205,7,-0.005608569703553211,0.5823156230017611,FoVP,-0.6507356749560859,0.7108871281275355,2457110.005950688 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-04-18T06:12:04.293,7,-0.006533791344582265,-0.7680639806381448,FoVP,0.3898713902485748,0.7177846905785434,2457130.7559392946 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-04-18T07:58:38.435,6,-0.003103768767807712,-0.7721436609584973,FoVF,0.3921103718707833,0.7157601708853536,2457130.829939783 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-08-31T06:43:00.914,1,0.004865564971320537,-2.4058116280271005,FoVP,-0.6878085837994734,0.7220322118950759,2457265.7795853484 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-08-31T08:29:35.136,5,-1.3801332074555073E-4,-2.407470299670282,FoVF,-0.6870095371493553,0.723135527289165,2457265.8535995437 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-10-21T06:18:54.018,1,0.004471299965674476,-2.8968391351998934,FoVP,-0.40082547260616475,0.7152650870733186,2457316.7675270303 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2015-10-21T08:05:28.195,4,5.314581390129491E-4,-2.893116607845149,FoVF,-0.40346347720117237,0.7129338131237428,2457316.84153943 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-24T02:22:49.398,7,-0.006579356942279041,0.236145041021683,FoVP,-0.03834907693935221,0.694428649438883,2457411.6042917464 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-24T08:23:03.958,6,-0.0038725920968205717,0.2546928790991775,FoVP,-0.052999549896909234,0.6964123003569312,2457411.85444359 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-24T10:09:38.234,6,-0.0031575038856307617,0.2601916339504578,FoVF,-0.057358402889244674,0.6969374272071714,2457411.9284464642 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-24T14:23:18.360,4,-0.0016147811680190876,0.27329939400961767,FoVP,-0.06777602079582802,0.6980729604950042,2457412.1045935275 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-24T16:09:52.595,5,-0.0010341575185062362,0.27881361290059464,FoVF,-0.07216954405264611,0.6985016850502285,2457412.178595899 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-24T20:23:32.633,3,1.8556436156794566E-4,0.2919545529356563,FoVP,-0.08266471202294118,0.699406370240556,2457412.3547418835 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-24T22:10:06.835,4,6.292818786538373E-4,0.29748114711269064,FoVF,-0.08708867642283567,0.6997375811664404,2457412.4287438514 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-25T04:10:20.983,4,0.0018248394530124305,0.3161833991809132,FoVF,-0.10210103063515528,0.7006418920781364,2457412.6788906604 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-25T08:24:00.906,2,0.0023818943010621763,0.32936757429373015,FoVP,-0.1127198448698961,0.7010784128387068,2457412.8550352077 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-25T10:10:35.068,3,0.0025450006018311066,0.33490894028439583,FoVF,-0.1171911823082803,0.7012122611944958,2457412.9290366755 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-25T14:24:14.967,2,0.00276305700835884,0.3481023527273083,FoVP,-0.12785522427154672,0.7014127131565167,2457413.1051808777 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-25T16:10:49.123,3,0.0027827504038084585,0.3536457808391573,FoVF,-0.13234316126613765,0.7014472524242198,2457413.179182254 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-25T20:24:29.019,2,0.002657189307499812,0.3668396702564774,FoVP,-0.1430409223748212,0.7014112029115603,2457413.355326366 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-25T22:11:03.178,3,0.0025316119796161395,0.37238140551992394,FoVF,-0.14754057165499984,0.7013463982839263,2457413.429327759 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-26T02:24:43.093,2,0.002058215448561157,0.3855666646156358,FoVP,-0.15826032082384694,0.701074122304188,2457413.605472041 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-26T04:11:17.266,4,0.0017856834223902193,0.39110281243277567,FoVF,-0.16276672159993388,0.7009102418277547,2457413.67947356 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-26T08:24:57.222,3,9.606596044898878E-4,0.40427001689142494,FoVP,-0.17349660196053956,0.7004027467082757,2457413.855618272 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-26T10:11:31.418,4,5.396729360711204E-4,0.40979655614814264,FoVF,-0.1780047614673837,0.7001403707809352,2457413.9296200266 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-26T14:25:11.439,4,-6.403204596819836E-4,0.42293599705504864,FoVP,-0.18873289084544134,0.6993994168962655,2457414.1057654307 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-26T16:11:45.665,5,-0.0012110681643103644,0.42844879422358995,FoVF,-0.19323782756796007,0.6990394454062332,2457414.179767531 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-26T20:25:25.773,5,-0.002748902097942923,0.4415505126895518,FoVP,-0.20395240322629904,0.6980675596501031,2457414.3559138877 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-26T22:12:00.041,7,-0.003470514929325013,0.4470453367865934,FoVF,-0.20844919124041658,0.6976112166242967,2457414.4299164424 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-01-27T02:25:40.259,6,-0.005368569191757951,0.46009916042917925,FoVP,-0.2191385960607972,0.6964117005468431,2457414.6060640095 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-03-07T19:39:51.913,3,0.0012121922309560186,-0.8542754870348369,FoVP,0.7041869581442349,0.7060162146074431,2457455.320674432 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-03-07T21:26:26.122,1,0.006665552464438079,-0.8542565304598202,FoVF,0.7040186443066007,0.7063161809805514,2457455.3946742327 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-04-12T14:55:55.710,1,0.005942206012869388,0.3628258985486866,FoVF,-0.48780566568372496,0.7085666572351743,2457491.1201509284 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-04-12T19:09:35.826,4,-7.294599410350257E-4,0.35199461456885456,FoVP,-0.4799952817800107,0.7120836520516343,2457491.296295099 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-04-12T20:56:10.002,7,-0.0034748396872719325,0.34740972259971487,FoVF,-0.47668869637983835,0.7135501447550359,2457491.3702956517 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-04-25T13:17:49.888,7,-0.00639600613720418,-0.5340341406546021,FoVP,0.1506039943362212,0.7184719412369589,2457504.051040776 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-04-25T15:04:24.008,7,-0.004297615005516056,-0.5389084680464333,FoVF,0.15350936115572406,0.7170256827047164,2457504.1250416706 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-04-25T19:18:03.957,3,7.591508228285377E-4,-0.5504829060947667,FoVP,0.16036430913827335,0.7135379064159008,2457504.3011863227 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-04-25T21:04:38.197,3,0.0029092713056413024,-0.5553339189506192,FoVF,0.1632188447704601,0.7120539222954,2457504.375188631 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-26T16:07:23.284,1,0.006574733665764367,-3.118795104425183,FoVF,0.10682475654855952,0.7218811268715385,2457596.168649828 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-26T22:07:37.867,2,0.005180252835275456,-3.1363485488714855,FoVF,0.11848620507860806,0.7229014193780243,2457596.4188352157 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-27T02:21:18.090,1,0.004362678166136526,3.134463758074025,FoVP,0.12674718221109432,0.7234981166164549,2457596.5950069344 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-27T04:07:52.375,2,0.004060258859464345,3.1292629484774115,FoVF,0.13022930699956975,0.7237185581011203,2457596.6690198025 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-27T08:21:32.554,1,0.0034392270401076255,3.116877979780593,FoVP,0.13854387784111935,0.724170859952826,2457596.8451910713 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-27T10:08:06.823,3,0.003220174398940672,3.111672752562402,FoVF,0.14204752721561414,0.7243303340855934,2457596.919203779 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-27T14:21:46.971,2,0.0027994552690540557,3.0992787745764105,FoVP,0.1504109179466158,0.7246368681123053,2457597.095374736 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-27T16:08:21.229,3,0.002665340738750663,3.0940704121321847,FoVF,0.1539340462749031,0.7247348503782834,2457597.169387342 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-27T20:22:01.358,2,0.0024486425767270693,3.081670569837051,FoVP,0.16234129836396347,0.7248944717745328,2457597.3455581274 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-27T22:08:35.611,3,0.0024010093550536505,3.0764604318182083,FoVF,0.16588178662740793,0.7249305255718157,2457597.4195706914 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-28T02:22:15.733,2,0.002391968504679075,3.0640580527793175,FoVP,0.17432777634488222,0.7249423130042384,2457597.5957414494 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-28T04:08:49.985,3,0.0024323280144230033,3.058847575269989,FoVF,0.17788344036041448,0.724916102015573,2457597.669754032 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-28T08:22:30.113,2,0.0026345009197534094,3.046446167753184,FoVP,0.18636289987487514,0.7247793795003656,2457597.8459249088 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-28T10:09:04.371,3,0.002764329210817298,3.0412368621777817,FoVF,0.1899314997313748,0.7246906732262365,2457597.9199375724 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-28T14:22:44.517,2,0.003181183779839112,3.0288401122270505,FoVP,0.19843904058849557,0.7244050230103958,2457598.0961087164 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-28T16:09:18.786,3,0.0034019178732009876,3.0236335639869374,FoVF,0.20201829092939788,0.7242537020485526,2457598.170121524 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-28T20:22:58.964,1,0.004036824717549616,3.0112453318580044,FoVP,0.2105484294938558,0.7238189768308214,2457598.346293086 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-28T22:09:33.248,2,0.004349858903315388,3.006043198698712,FoVF,0.21413601060055207,0.7236050378553299,2457598.4203061014 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-07-29T04:09:47.778,1,0.005612764589441238,2.9884715247663682,FoVF,0.22627676495516516,0.7227449326322577,2457598.6704915203 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-09-10T07:53:21.077,3,2.5697815122738374E-4,-2.4290032247557622,FoVP,-0.7192512483154985,0.7205751339300248,2457641.8294886546 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-10-17T19:13:04.480,2,0.0019670414371021298,2.3267191540506578,FoVP,0.4528101981918288,0.7141797978884794,2457679.304956016 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-10-17T20:59:38.636,2,0.0049775136690584825,2.3312166991054286,FoVF,0.44890664506950456,0.7158274316724391,2457679.378968362 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-10-28T09:06:29.801,3,0.003469437222110181,3.065978956067008,FoVF,-0.180726009925592,0.712925193808823,2457689.884508299 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2016-10-28T13:20:09.807,5,-0.002756221733663789,3.0773433937782557,FoVP,-0.18931149002963266,0.7086098902236475,2457690.060677803 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2017-01-28T21:19:07.703,5,-0.0022362656409328804,-0.26838445629601826,FoVP,0.29847350987287763,0.6975925921145726,2457782.3930357434 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2017-01-28T23:05:41.920,4,0.0011818264198154163,-0.2640267638820214,FoVF,0.2953490008385252,0.6997892657822943,2457782.4670374626 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2017-02-12T04:47:43.319,4,-2.97168219826661E-4,0.7062620273648472,FoVP,-0.5481689981488462,0.7034066986201087,2457796.7034035507 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2017-02-12T06:34:17.557,7,-0.0044351025601507554,0.709904993482966,FoVF,-0.5513672776360866,0.7016248404356326,2457796.7774044797 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2017-03-18T10:12:07.737,6,-0.004262440417520913,-0.891777649505845,FoVP,0.697256372367284,0.7091677863936382,2457830.925425397 +HD 33632,1.3669889772966024,0.6516590165506833,05:13:17.448,+37:20:14.321,2017-03-18T11:58:41.940,4,0.0010442329234670706,-0.8928944307050657,FoVF,0.6977166434167678,0.7085186421825921,2457830.9994251765 diff --git a/orbitize/gaia.py b/orbitize/gaia.py index af957b7c..d232610d 100644 --- a/orbitize/gaia.py +++ b/orbitize/gaia.py @@ -7,7 +7,7 @@ from astroquery.gaia import Gaia from astropy import units as u import astropy.io.fits as fits -import astropy.time as time +from astropy.time import Time from astropy.io.ascii import read from astropy.coordinates import get_body_barycentric_posvel import numpy.linalg @@ -61,12 +61,12 @@ def __init__(self, gaia_num, hiplogprob, dr="dr2", query=True, gaia_data=None): self.dr = dr if self.dr == "edr3": - self.gaia_epoch = 2016.0 + self.gaia_epoch = Time("J2016", format='jyear_str') elif self.dr == "dr2": - self.gaia_epoch = 2015.5 + self.gaia_epoch = Time("J2015.5", format='jyear_str') else: raise ValueError("`dr` must be either `dr2` or `edr3`") - self.hipparcos_epoch = 1991.25 + self.hipparcos_epoch = Time("J1991.25", format='jyear_str') if query: query = """SELECT @@ -134,12 +134,12 @@ def compute_lnlike(self, raoff_model, deoff_model, samples, param_idx): """ alpha_H0 = samples[param_idx["alpha0"]] # [deg] - pm_ra = samples[param_idx["pm_ra"]] # [mas/yr] - delta_alpha_from_pm = pm_ra * (self.gaia_epoch - self.hipparcos_epoch) # [mas] + pm_ra = samples[param_idx["pm_ra"]] # [mas/jyr] + delta_alpha_from_pm = pm_ra * (self.gaia_epoch.jyear - self.hipparcos_epoch.jyear) # [mas] delta_H0 = samples[param_idx["delta0"]] # [deg] - pm_dec = samples[param_idx["pm_dec"]] # [mas/yr] - delta_delta_from_pm = pm_dec * (self.gaia_epoch - self.hipparcos_epoch) # [mas] + pm_dec = samples[param_idx["pm_dec"]] # [mas/jyr] + delta_delta_from_pm = pm_dec * (self.gaia_epoch.jyear - self.hipparcos_epoch.jyear) # [mas] # difference in position due to orbital motion between Hipparcos & Gaia epochs alpha_diff_orbit = raoff_model[1, :] - raoff_model[0, :] # [mas] @@ -162,6 +162,7 @@ def compute_lnlike(self, raoff_model, deoff_model, samples, param_idx): # technically this is an angle so we should wrap it, but the precision # of Hipparcos and Gaia is so good that we'll never have to. alpha_resid = alpha_model - alpha_data + alpha_chi2 = (alpha_resid / alpha_unc) ** 2 delta_model = self.hiplogprob.delta0 + self.mas2deg * ( # [deg] @@ -296,9 +297,9 @@ def __init__(self, hip_id, hiplogprob, gost_filepath, hgca_filepath=None): self.gaia_epoch_dec = entry["epoch_dec_gaia"][0] # read in the GOST file to get the estimated Gaia epochs and scan angles gost_dat = read(gost_filepath, converters={"*": [int, float, bytes]}) - self.gaia_epoch = time.Time( + self.gaia_epoch = Time( gost_dat["ObservationTimeAtGaia[UTC]"] - ).decimalyear # in decimal year + ) # in julian year gaia_scan_theta = np.array(gost_dat["scanAngle[rad]"]) gaia_scan_phi = gaia_scan_theta + np.pi / 2 self.gaia_cos_phi = np.cos(gaia_scan_phi) @@ -443,9 +444,9 @@ def _linear_pm_fit( """ # Sovle y = A * x # construct A matrix - A_pmra = cos_phi * (epochs - epoch_ref_ra) / errs + A_pmra = cos_phi * (epochs.jyear - epoch_ref_ra) / errs A_raoff = cos_phi / errs - A_pmdec = sin_phi * (epochs - epoch_ref_dec) / errs + A_pmdec = sin_phi * (epochs.jyear - epoch_ref_dec) / errs A_decoff = sin_phi / errs A_matrix = np.vstack((A_raoff, A_decoff, A_pmra, A_pmdec)).T diff --git a/orbitize/hipparcos.py b/orbitize/hipparcos.py index 02b23c45..94908548 100644 --- a/orbitize/hipparcos.py +++ b/orbitize/hipparcos.py @@ -1,5 +1,4 @@ import numpy as np -from astropy.io import ascii import pandas as pd import emcee from scipy.stats import norm @@ -17,17 +16,17 @@ class PMPlx_Motion(object): parallax and proper motion model (NO orbital motion is added in this class). Args: - times_mjd (np.array of float): times (in mjd) at which we have absolute astrometric + epochs_mjd (np.array of float): times (in mjd) at which we have absolute astrometric measurements alpha0 (float): measured RA position (in degrees) of the object at alphadec0_epoch (see below). delta0 (float): measured Dec position (in degrees) of the object at alphadec0_epoch (see below). - alphadec0_epoch (float): a (fixed) reference time. For stars with Hipparcos data, this - should generally be 1991.25, but you can define it however you want. Absolute + alphadec0_epoch (float): a (fixed) reference time (in jyear). For stars with Hipparcos data, this + should generally be J1991.25, but you can define it however you want. Absolute astrometric data (passed in via an orbitize! data table) should be defined as offsets from the reported position of the object at this epoch (with propagated uncertainties). For example, if you have two absolute astrometric measurements - taken with GRAVITY, as well as a Hipparcos-derived position (at epoch 1991.25), - alphadec0_epoch should be 1991.25, and you should pass in absolute astrometry + taken with GRAVITY, as well as a Hipparcos-derived position (at epoch J1991.25), + alphadec0_epoch should be J1991.25, and you should pass in absolute astrometry in terms of mas *offset* from the Hipparcos catalog position, with propagated errors of your measurement and the Hipparcos measurement. """ @@ -39,7 +38,6 @@ def __init__(self, epochs_mjd, alpha0, delta0, alphadec0_epoch=1991.25): self.delta0 = delta0 epochs = Time(epochs_mjd, format="mjd") - self.epochs = epochs.decimalyear # compute Earth XYZ position in barycentric coordinates bary_pos, _ = get_body_barycentric_posvel("earth", epochs) @@ -61,7 +59,7 @@ def compute_astrometric_model(self, samples, param_idx, epochs=None): indices in an array of fitting parameters (generally set to System.basis.param_idx). epochs: if None, use self.epochs for astrometric predictions. Otherwise, - use this array passed in [in decimalyear]. + use this array passed in [in mjd]. Returns: tuple of: @@ -78,14 +76,14 @@ def compute_astrometric_model(self, samples, param_idx, epochs=None): delta_H0 = samples[param_idx["delta0"]] if epochs is None: - epochs = self.epochs + epochs = Time(self.epochs_mjd, format='mjd') X = self.X Y = self.Y Z = self.Z else: # compute Earth XYZ position in barycentric coordinates bary_pos, _ = get_body_barycentric_posvel( - "earth", Time(epochs, format="decimalyear") + "earth", Time(epochs, format="mjd") ) X = bary_pos.x.value # [au] Y = bary_pos.y.value # [au] @@ -105,7 +103,7 @@ def compute_astrometric_model(self, samples, param_idx, epochs=None): X[i] * np.sin(np.radians(self.alpha0)) - Y[i] * np.cos(np.radians(self.alpha0)) ) - + (epochs[i] - self.alphadec0_epoch) * pm_ra + + (epochs[i].jyear - self.alphadec0_epoch) * pm_ra ) delta_C_array[i] = ( delta_H0 @@ -119,7 +117,7 @@ def compute_astrometric_model(self, samples, param_idx, epochs=None): * np.sin(np.radians(self.delta0)) - Z[i] * np.cos(np.radians(self.delta0)) ) - + (epochs[i] - self.alphadec0_epoch) * pm_dec + + (epochs[i].jyear - self.alphadec0_epoch) * pm_dec ) return alpha_C_st_array, delta_C_array @@ -137,9 +135,9 @@ class HipparcosLogProb(object): are described here for completeness. See Nielsen+ 2020 for more detail. - alpha0: RA offset from the reported Hipparcos position at a particular - epoch (usually 1991.25) [mas] + epoch (usually J1991.25) [mas] - delta0: Dec offset from the reported Hipparcos position at a particular - epoch (usually 1991.25) [mas] + epoch (usually J1991.25) [mas] - pm_ra: RA proper motion [mas/yr] - pm_dec: Dec proper motion [mas/yr] - plx: parallax [mas] @@ -157,7 +155,7 @@ class HipparcosLogProb(object): zeros in the prefix if number is <100,000. (i.e. 27321 should be passed in as '027321'). num_secondary_bodies (int): number of companions in the system - alphadec0_epoch (float): epoch (in decimal year) that the fitting + alphadec0_epoch (float): epoch (in Julian decimal year) that the fitting parameters alpha0 and delta0 are defined relative to (see above). renormalize_errors (bool): if True, normalize the scan errors to get chisq_red = 1, following Nielsen+ 2020 (eq 10). In general, this @@ -173,7 +171,7 @@ def __init__( path_to_iad_file, hip_num, num_secondary_bodies, - alphadec0_epoch=1991.25, + alphadec0_epoch=Time("J1991.25", format="jyear_str").jyear, renormalize_errors=False, ): self.path_to_iad_file = path_to_iad_file @@ -284,7 +282,11 @@ def __init__( n_lines = len(iad) - times = iad[1] + 1991.25 + + self.epochs = Time( + iad[1] * 365.25 + Time("J1991.25",format="jyear_str").jd, format='jd' + ) + self.cos_phi = iad[3] # scan direction self.sin_phi = iad[4] self.R = iad[5] # abscissa residual [mas] @@ -295,18 +297,18 @@ def __init__( if n_lines - len(good_scans) > 0: print("{} Hipparcos scans rejected.".format(n_lines - len(good_scans))) - times = times[good_scans] + self.epochs = self.epochs[good_scans] self.cos_phi = self.cos_phi[good_scans] self.sin_phi = self.sin_phi[good_scans] self.R = self.R[good_scans] self.eps = self.eps[good_scans] + self.epochs_mjd = self.epochs.mjd + # if the star has a type 1 (stochastic) solution, we need to undo the addition of a jitter term in quadrature self.eps = np.sqrt(self.eps**2 - self.var**2) - epochs = Time(times, format="decimalyear") - self.epochs = epochs.decimalyear - self.epochs_mjd = epochs.mjd + self.hipparcos_plxpm_predictor = PMPlx_Motion( self.epochs_mjd, @@ -316,7 +318,7 @@ def __init__( ) if self.renormalize_errors: - D = len(epochs) - 6 + D = len(self.epochs) - 6 G = f2 f = (G * np.sqrt(2 / (9 * D)) + 1 - (2 / (9 * D))) ** (3 / 2) @@ -333,7 +335,7 @@ def __init__( self.hipparcos_plxpm_predictor.X * np.sin(np.radians(self.alpha0)) - self.hipparcos_plxpm_predictor.Y * np.cos(np.radians(self.alpha0)) ) - + (self.epochs - 1991.25) * self.pm_ra0 + + (self.epochs.jyear - self.alphadec0_epoch) * self.pm_ra0 ) changein_delta = ( @@ -347,7 +349,7 @@ def __init__( * np.sin(np.radians(self.delta0)) - self.hipparcos_plxpm_predictor.Z * np.cos(np.radians(self.delta0)) ) - + (self.epochs - 1991.25) * self.pm_dec0 + + (self.epochs.jyear - self.alphadec0_epoch) * self.pm_dec0 ) # compute abcissa point (Nielsen+ Eq 3) diff --git a/orbitize/plot.py b/orbitize/plot.py index afaf85c1..b72de050 100644 --- a/orbitize/plot.py +++ b/orbitize/plot.py @@ -1675,7 +1675,7 @@ def plot_propermotion( axs[0].set_ylabel(r"$\mu_\alpha^*$ [mas/yr]") axs[0].errorbar( - np.nanmedian(system.gaia.hipparcos_epoch), + np.nanmedian(system.gaia.hipparcos_epoch.jyear), system.gaia.hip_pm[0], yerr=system.gaia.hip_pm_err[0], zorder=30, @@ -1685,10 +1685,10 @@ def plot_propermotion( ) hgca_epoch = ( - system.gaia.gaia_epoch_ra + np.nanmedian(system.gaia.hipparcos_epoch) + system.gaia.gaia_epoch_ra + np.nanmedian(system.gaia.hipparcos_epoch.jyear) ) / 2 hgca_epoch_err = ( - system.gaia.gaia_epoch_ra - np.nanmedian(system.gaia.hipparcos_epoch) + system.gaia.gaia_epoch_ra - np.nanmedian(system.gaia.hipparcos_epoch.jyear) ) / 2 axs[0].errorbar( @@ -1716,7 +1716,7 @@ def plot_propermotion( axs[1].yaxis.set_major_formatter(FormatStrFormatter("%.1f")) axs[1].errorbar( - np.nanmedian(system.gaia.hipparcos_epoch), + np.nanmedian(system.gaia.hipparcos_epoch.jyear), system.gaia.hip_pm[1], yerr=system.gaia.hip_pm_err[1], zorder=30, diff --git a/orbitize/results.py b/orbitize/results.py index fea92e13..b3f5c117 100644 --- a/orbitize/results.py +++ b/orbitize/results.py @@ -48,6 +48,19 @@ def __init__( self.ln_evidence = None self.ln_evidence_err = None + # bookkeeping for save_results(): lets repeated calls for the same + # output file append only the newly added rows instead of rewriting + # the whole file from scratch every time + self._saved_filename = None + self._n_rows_saved = 0 + + # backing buffers for add_samples(): allow post/lnlike to grow with + # amortized O(1) cost per row instead of reallocating+copying the + # full array on every call (see add_samples() docstring) + self._post_buf = None + self._lnlike_buf = None + self._n_used = 0 + if self.system is not None: self.tau_ref_epoch = self.system.tau_ref_epoch self.labels = self.system.labels @@ -58,40 +71,88 @@ def __init__( self.param_idx = self.system.param_idx self.standard_param_idx = self.system.basis.standard_basis_idx - def add_samples(self, orbital_params, lnlikes, curr_pos=None): + def add_samples(self, orbital_params, lnlikes, curr_pos=None): """ - Add accepted orbits, their likelihoods, and the orbitize version number + Add accepted orbits, their likelihoods, and the orbitize version number to the results Args: - orbital_params (np.array): add sets of orbital params (could be multiple) + orbital_params (np.array): add sets of orbital params (could be multiple) to results lnlike (np.array): add corresponding lnlike values to results - curr_pos (np.array of float): for MCMC only. A multi-D array of the + curr_pos (np.array of float): for MCMC only. A multi-D array of the current walker positions + .. Note:: ``post``/``lnlike`` are backed by an internal buffer that's + over-allocated and grown by doubling, so repeated calls (e.g. from + ``periodic_save_freq`` during MCMC) append in amortized O(1) time + per row instead of reallocating and copying the full accumulated + array on every call. + Written: Henry Ngo, 2018 API Update: Sarah Blunt, 2021 """ - + # Adding the orbitize version number to the results if self.version_number is None: self.version_number = orbitize.__version__ - # If no exisiting results then it is easy - if self.post is None: - self.post = orbital_params - self.lnlike = lnlikes + n_new = len(orbital_params) - # Otherwise, need to append properly - else: - self.post = np.vstack((self.post, orbital_params)) - self.lnlike = np.append(self.lnlike, lnlikes) + # lazily adopt any existing post/lnlike (e.g. set directly via the + # constructor) as the starting buffer + if self._post_buf is None: + if self.post is not None: + self._post_buf = self.post + self._lnlike_buf = self.lnlike + self._n_used = len(self.post) + else: + self._post_buf = np.empty((0, orbital_params.shape[1]), dtype=orbital_params.dtype) + self._lnlike_buf = np.empty(0, dtype=lnlikes.dtype) + self._n_used = 0 + + needed = self._n_used + n_new + capacity = self._post_buf.shape[0] + if needed > capacity: + new_capacity = max(needed, capacity * 2) + + new_post_buf = np.empty((new_capacity,) + self._post_buf.shape[1:], dtype=self._post_buf.dtype) + new_post_buf[: self._n_used] = self._post_buf[: self._n_used] + self._post_buf = new_post_buf + + new_lnlike_buf = np.empty(new_capacity, dtype=self._lnlike_buf.dtype) + new_lnlike_buf[: self._n_used] = self._lnlike_buf[: self._n_used] + self._lnlike_buf = new_lnlike_buf + + self._post_buf[self._n_used : needed] = orbital_params + self._lnlike_buf[self._n_used : needed] = lnlikes + self._n_used = needed + + self.post = self._post_buf[: self._n_used] + self.lnlike = self._lnlike_buf[: self._n_used] if curr_pos is not None: self.curr_pos = curr_pos + def _write_growable_dataset(self, hf, name, data, n_already_saved): + """ + Writes ``data`` to the dataset ``name`` in the open hdf5 file ``hf``. + + If the dataset doesn't exist yet, it's created as resizable so that + future calls can extend it. If it already exists (i.e. a previous + call to ``save_results`` already wrote the first ``n_already_saved`` + rows of ``data`` to this same file), only the new rows beyond that + point are written, instead of rewriting the whole dataset. + """ + if name in hf: + dset = hf[name] + dset.resize(len(data), axis=0) + dset[n_already_saved:] = data[n_already_saved:] + else: + maxshape = (None,) + data.shape[1:] + hf.create_dataset(name, data=data, maxshape=maxshape, chunks=True) + def save_results(self, filename): """ Save results.Results object to an hdf5 file @@ -105,12 +166,26 @@ def save_results(self, filename): ``post``, ``lnlike``, and ``parameter_labels`` are datasets that are members of the root group. + If called repeatedly with the same ``filename`` (e.g. from + ``periodic_save_freq`` during MCMC), only the rows of ``post``/``lnlike`` + that were added since the last call are written to disk, rather than + rewriting the entire accumulated chain every time. + Written: Henry Ngo, 2018 API Update: Sarah Blunt, 2021 """ - hf = h5py.File(filename, 'w') # Creates h5py file object + # if this is a new target file (or the first save), start fresh; + # otherwise reopen the file we've already been writing to and append + if filename != self._saved_filename: + mode = 'w' + self._n_rows_saved = 0 + self._saved_filename = filename + else: + mode = 'a' + + hf = h5py.File(filename, mode) # Creates/opens h5py file object # Add sampler_name as attribute of the root group hf.attrs['sampler_name'] = self.sampler_name @@ -124,15 +199,22 @@ def save_results(self, filename): # Now add post and lnlike from the results object as datasets if self.post is not None: - hf.create_dataset('post', data=self.post) + self._write_growable_dataset(hf, 'post', self.post, self._n_rows_saved) # hf.create_dataset('data', data=self.data) if self.lnlike is not None: - hf.create_dataset('lnlike', data=self.lnlike) + self._write_growable_dataset(hf, 'lnlike', self.lnlike, self._n_rows_saved) + + self._n_rows_saved = len(self.post) if self.post is not None else 0 if self.curr_pos is not None: + if 'curr_pos' in hf: + del hf['curr_pos'] hf.create_dataset("curr_pos", data=self.curr_pos) - self.system.save(hf) + # the system configuration doesn't change between saves, so it only + # needs to be (re)written the first time we save to this file + if mode == 'w': + self.system.save(hf) hf.close() # Closes file object, which writes file to disk diff --git a/orbitize/sampler.py b/orbitize/sampler.py index 09c23cc3..3e9e3c49 100644 --- a/orbitize/sampler.py +++ b/orbitize/sampler.py @@ -123,12 +123,9 @@ def _logl(self, params): ) if self.system.gaia is not None: - gaiahip_epochs = Time( - np.append( - self.system.gaia.hipparcos_epoch, self.system.gaia.gaia_epoch - ), - format="decimalyear", - ).mjd + gaiahip_epochs = np.append( + self.system.gaia.hipparcos_epoch.mjd, self.system.gaia.gaia_epoch.mjd + ) # compute Ra/Dec predictions at the Gaia epoch raoff_model, deoff_model, _, _ = self.system.compute_all_orbits( @@ -849,13 +846,19 @@ def _logl(self, params, include_logp=False): return super(MCMC, self)._logl(full_params) + logp - def _update_chains_from_sampler(self, sampler, num_steps=None): + def _update_chains_from_sampler(self, sampler, num_steps=None, start_step=0): """ Updates self.post, self.chain, and self.lnlike from the MCMC sampler Args: sampler (emcee.EnsembleSampler or ptemcee.Sampler): sampler object. num_steps (int): if not None, only stores the first num_steps number of steps + start_step (int): if not 0, only (re)computes steps from this index + onwards, so self.post/self.lnlikes end up covering just the + ``[start_step, num_steps)`` chunk of the chain instead of the + full history. Used during periodic saving so that the (costly, + unvectorized) prior recomputation below doesn't get redone over + steps that were already processed in a previous call. """ if num_steps is None: # use all the steps, grab total number of steps from dimension of chains @@ -866,16 +869,16 @@ def _update_chains_from_sampler(self, sampler, num_steps=None): if self.use_pt: # chain is shape: Ntemp x Nwalkers x Nsteps x Nparams - self.post = sampler.chain[0, :, :num_steps].reshape( + self.post = sampler.chain[0, :, start_step:num_steps].reshape( -1, num_params ) # the reshaping flattens the chain # should also be picking out the lowest temperature logps - self.lnlikes = sampler.loglikelihood[0, :, :num_steps].flatten() - self.lnlikes_alltemps = sampler.loglikelihood[:, :, :num_steps] + self.lnlikes = sampler.loglikelihood[0, :, start_step:num_steps].flatten() + self.lnlikes_alltemps = sampler.loglikelihood[:, :, start_step:num_steps] else: # chain is shape: Nwalkers x Nsteps x Nparams - self.post = sampler.chain[:, :num_steps].reshape(-1, num_params) - self.lnlikes = sampler.lnprobability[:, :num_steps].flatten() + self.post = sampler.chain[:, start_step:num_steps].reshape(-1, num_params) + self.lnlikes = sampler.lnprobability[:, start_step:num_steps].flatten() # convert posterior probability (returned by sampler objects) to likelihood (required by orbitize.results.Results) for i, orb in enumerate(self.post): @@ -1001,10 +1004,10 @@ def run_sampler( self._logl, orbitize.priors.all_lnpriors, ntemps=self.num_temps, - threads=self.num_threads, logpargs=[ self.priors, ], + pool=pool ) else: sampler = emcee.EnsembleSampler( @@ -1057,68 +1060,45 @@ def run_sampler( if periodic_save_freq is not None: if (i + 1) % periodic_save_freq == 0: # we've completed i+1 steps - self._update_chains_from_sampler(sampler, num_steps=i + 1) - - # figure out what is the new chunk of the chain and corresponding lnlikes that have been computed before last save - # grab the current posterior and lnlikes and reshape them to have the Nwalkers x Nsteps dimension again - post_shape = self.post.shape - curr_chain_shape = ( - self.num_walkers, - post_shape[0] // self.num_walkers, - post_shape[-1], + # only (re)compute the chunk of the chain since the last + # save, not the full history -- avoids redoing the prior + # recomputation over steps that were already saved + self._update_chains_from_sampler( + sampler, num_steps=i + 1, start_step=saved_upto ) - curr_chain = self.post.reshape(curr_chain_shape) - curr_lnlike_chain = self.lnlikes.reshape(curr_chain_shape[:2]) - # use the reshaped arrays and find the new steps we computed - curr_chunk = curr_chain[:, saved_upto : i + 1] - curr_chunk = curr_chunk.reshape( - -1, curr_chunk.shape[-1] - ) # flatten nwalkers x nsteps dim - curr_lnlike_chunk = curr_lnlike_chain[ - :, saved_upto : i + 1 - ].flatten() # add this current chunk to the results object (which already has all the previous chunks saved) self.results.add_samples( - curr_chunk, curr_lnlike_chunk, curr_pos=self.curr_pos + self.post, self.lnlikes, curr_pos=self.curr_pos ) self.results.save_results(output_filename) saved_upto = i + 1 print("") - self._update_chains_from_sampler(sampler) if periodic_save_freq is None: - # need to save everything + # nothing has been processed/saved yet; need to do the whole chain + self._update_chains_from_sampler(sampler) self.results.add_samples( self.post, self.lnlikes, curr_pos=self.curr_pos ) elif saved_upto < nsteps: - # just need to save the last few - # same code as above except we just need to grab the last few - post_shape = self.post.shape - curr_chain_shape = ( - self.num_walkers, - post_shape[0] // self.num_walkers, - post_shape[-1], - ) - curr_chain = self.post.reshape(curr_chain_shape) - curr_lnlike_chain = self.lnlikes.reshape(curr_chain_shape[:2]) - curr_chunk = curr_chain[:, saved_upto:] - curr_chunk = curr_chunk.reshape( - -1, curr_chunk.shape[-1] - ) # flatten nwalkers x nsteps dim - curr_lnlike_chunk = curr_lnlike_chain[:, saved_upto:].flatten() - + # just need to process and save the leftover chunk since the last save + self._update_chains_from_sampler(sampler, start_step=saved_upto) self.results.add_samples( - curr_chunk, curr_lnlike_chunk, curr_pos=self.curr_pos + self.post, self.lnlikes, curr_pos=self.curr_pos ) + else: + # everything has already been saved; still refresh self.chain, + # self.post, and self.lnlikes to reflect the full run in case + # the caller inspects them after run_sampler() returns + self._update_chains_from_sampler(sampler) if output_filename is not None: self.results.save_results(output_filename) print("Run complete") - # Close pool + if examine_chains: self.examine_chains() diff --git a/tests/end-to-end-tests/hd_33632Ab.py b/tests/end-to-end-tests/hd_33632Ab.py new file mode 100644 index 00000000..fa78e72e --- /dev/null +++ b/tests/end-to-end-tests/hd_33632Ab.py @@ -0,0 +1,82 @@ +""" +Compare to Hsu+ 2024 (case with HGCA, no RVs, and relative astrometry) + +https://iopscience.iop.org/article/10.3847/1538-4357/ad58d3#apjad58d3t5 +""" + +import os +from orbitize import DATADIR, hipparcos, gaia, read_input, system, priors, sampler, results +import matplotlib.pyplot as plt + +# the necessary input data for beta Pic is part of the orbitize! example data! +iad_filepath = os.path.join(DATADIR, "H024332.d") +gost_filepath = os.path.join(DATADIR, "gost_22.4.3_806543_2026-07-28-17-46-26_HD_33632.csv") + +# Create the HGCA and helper Hipparcos object +hipparcos_lnprob = hipparcos.HipparcosLogProb(iad_filepath, 24332, 1) +hgca_lnprob = gaia.HGCALogProb(24332, hipparcos_lnprob, gost_filepath) + +# read in relative astrometry +astrometry_filepath = os.path.join(DATADIR, "HD33632Ab.csv") +data_table = read_input.read_file(astrometry_filepath) + +# set up the system, passing in hgca_lnprob and setting it fit dynamical mass +stellar_mass = 1.11 +stellar_mass_err = 0.09 +plx = 37.8953 +plx_err = 0.0263 + +this_system = system.System( + 1, + data_table, + stellar_mass, + plx, + mass_err=stellar_mass_err, + plx_err=plx_err, + fit_secondary_mass=True, + gaia=hgca_lnprob, +) + +# adjust the prior on mass to be uniform between 0 and 0.1 Msol +# this_system.sys_priors[this_system.param_idx["m1"]] = priors.LogUniformPrior( +# 0, 0.1 +# ) + +# MCMC parameters +n_temps=20 +n_walkers=1000 +n_threads=20 +total_orbits= n_walkers * 50_000 +burn_steps=10_000 +thin=100 + +run_fit = True + +if __name__ == '__main__': + + # create the sampler, run it, and save posteriors + output_filename = "HD_33632_Ab.hdf5" + if os.path.exists(output_filename): + prev_result_filename = output_filename + else: + prev_result_filename = None + + this_sampler = sampler.MCMC(this_system, n_temps, n_walkers, n_threads, prev_result_filename=prev_result_filename) + + periodic_save_freq = 5_000 + + if run_fit: + + this_sampler.run_sampler( + total_orbits, burn_steps=burn_steps, thin=thin, periodic_save_freq=periodic_save_freq, + output_filename=output_filename + ) + + this_sampler.results.save_results(output_filename) + + myResults = results.Results() + myResults.load_results(output_filename) + + # make corner plot + fig = myResults.plot_corner() + plt.savefig("HD_33632_Ab.png", dpi=250) \ No newline at end of file diff --git a/tests/test_abs_astrometry.py b/tests/test_abs_astrometry.py index 84f9dc0c..effdef8d 100644 --- a/tests/test_abs_astrometry.py +++ b/tests/test_abs_astrometry.py @@ -88,10 +88,9 @@ def test_arbitrary_abs_astrom(): path_to_iad_file, hip_num, num_secondary_bodies ) - epochs_astropy = Time( - np.array([0, 0.5, 1.0]) + testHiPIAD.alphadec0_epoch, format="decimalyear" - ) - epochs = epochs_astropy.mjd + epochs = Time( + testHiPIAD.alphadec0_epoch, format="jyear" + ).mjd + np.array([0, 365.25/2, 365.25]) ra_model = np.zeros(epochs.shape) dec_model = np.zeros(epochs.shape) @@ -190,5 +189,5 @@ def test_arbitrary_abs_astrom(): if __name__ == "__main__": - test_1planet() + # test_1planet() test_arbitrary_abs_astrom() diff --git a/tests/test_gaia.py b/tests/test_gaia.py index 9f99a583..fc41e34a 100644 --- a/tests/test_gaia.py +++ b/tests/test_gaia.py @@ -156,11 +156,11 @@ def test_orbit_calculation(): myGaia.ra = myHip.alpha0 + ( myGaia.mas2deg * pm_a - * (myGaia.gaia_epoch - myGaia.hipparcos_epoch) + * (myGaia.gaia_epoch.jyear - myGaia.hipparcos_epoch.jyear) / np.cos(np.radians(myHip.delta0)) ) myGaia.dec = myHip.delta0 + ( - myGaia.mas2deg * pm_d * (myGaia.gaia_epoch - myGaia.hipparcos_epoch) + myGaia.mas2deg * pm_d * (myGaia.gaia_epoch.jyear - myGaia.hipparcos_epoch.jyear) ) test_samples = [sma, ecc, inc, aop, pan, tau, plx, m0, m1, a0, d0, pm_a, pm_d] @@ -194,7 +194,7 @@ def test_orbit_calculation(): myGaia.dec = myHip.delta0 + 1 sma = 2 * (myGaia.dec - myHip.delta0) * deg2arcsec * (plx * mas2arcsec) # [au] - per = 2 * (myGaia.gaia_epoch - myGaia.hipparcos_epoch) # [yr] + per = 2 * (myGaia.gaia_epoch.jyear - myGaia.hipparcos_epoch.jyear) # [yr] mtot = sma**3 / per**2 test_samples[param_idx["sma1"]] = sma @@ -203,7 +203,7 @@ def test_orbit_calculation(): # passes through peri (+sma decl for e=0 orbits) at Hipparcos epoch # -> @ Gaia epoch, primary should be at +sma decl - tau = basis.tp_to_tau(myGaia.hipparcos_epoch, 58849, per) + tau = basis.tp_to_tau(myGaia.hipparcos_epoch.jyear, 58849, per) test_samples[param_idx["tau1"]] = tau # choose sma and mass so that Hipparcos/Gaia difference is only due to orbit. @@ -308,8 +308,8 @@ def test_nointernet(): if __name__ == "__main__": test_nointernet() - # test_dr2_edr3() - # test_system_setup() - # test_valueerror() - # test_orbit_calculation() + test_dr2_edr3() + test_system_setup() + test_valueerror() + test_orbit_calculation() test_hgca() diff --git a/tests/test_hipparcos.py b/tests/test_hipparcos.py index 08299835..5c401a47 100644 --- a/tests/test_hipparcos.py +++ b/tests/test_hipparcos.py @@ -148,7 +148,7 @@ def test_dvd_vs_2021catalog(): np.isclose(new_iadHipLogProb.sin_phi, old_iadHipLogProb.sin_phi, atol=1e-2) ) assert np.all( - np.isclose(new_iadHipLogProb.epochs, old_iadHipLogProb.epochs, atol=1e-2) + np.isclose(new_iadHipLogProb.epochs.value, old_iadHipLogProb.epochs.value, atol=1e-2) ) @@ -268,5 +268,5 @@ def test_save_load_2021(): # test_save_load_dvd() # test_save_load_2021() # test_hipparcos_api() - test_iad_refitting() - # test_dvd_vs_2021catalog() + # test_iad_refitting() + test_dvd_vs_2021catalog()