From 7da879f10ce0d3cc2051deb343cbf341a7ae5d87 Mon Sep 17 00:00:00 2001 From: Steph Merritt <97111051+astronomerritt@users.noreply.github.com> Date: Tue, 5 Nov 2024 12:12:17 +0000 Subject: [PATCH] Adding adler-day draft script (#180) * adler day draft script * linting * adler cli now runs adler_run.py * Moving database, fixing bug. --------- Co-authored-by: jrob93 --- pyproject.toml | 2 +- src/adler/adler_day.py | 117 +++++++++++++++++++++++++++ src/adler/{adler.py => adler_run.py} | 0 tests/data/adler-day-testing.db | Bin 0 -> 16384 bytes 4 files changed, 118 insertions(+), 1 deletion(-) create mode 100644 src/adler/adler_day.py rename src/adler/{adler.py => adler_run.py} (100%) create mode 100644 tests/data/adler-day-testing.db diff --git a/pyproject.toml b/pyproject.toml index 7874b80..b5742c4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -46,7 +46,7 @@ dev = [ ] [project.scripts] -adler = "adler.adler:main" +adler = "adler.adler_run:main" [build-system] requires = [ diff --git a/src/adler/adler_day.py b/src/adler/adler_day.py new file mode 100644 index 0000000..2a5ae3b --- /dev/null +++ b/src/adler/adler_day.py @@ -0,0 +1,117 @@ +from adler.objectdata.AdlerPlanetoid import AdlerPlanetoid +from adler.science.PhaseCurve import PhaseCurve +import os +import pandas as pd +import astropy.units as u +import numpy as np + + +def run_adler_day(ssobject_list): + """Runs the Adler day process for a list of SSObjectIds. + + Parameters + ----------- + ssobject_list : list of str + List of strings of the SSObjectIDs we want to investigate. + + """ + + phase_model = "HG12_Pen16" + + for ssoid in ssobject_list: + + # date range is start of survey to "today's" date, these are dummy dates + # this should also pull in the most recent alert packet, which should be in Cassandra by now? + planetoid = AdlerPlanetoid.construct_from_cassandra(ssoid, date_range=[60290.0, 61902]) + + try: + dirname = os.path.dirname(__file__) + database_path = os.path.join(dirname, "../../tests/data/adler-day-testing.db") + planetoid.attach_previous_adler_data(database_path) + previousAdlerData = True # this could be a flag on the AdlerPlanetoid object + except ValueError: + print("No data found for this object in AdlerDatabase.") + previousAdlerData = False + + for filt in planetoid.filter_list: + obs = planetoid.observations_in_filter(filt) + df_obs = pd.DataFrame(obs.__dict__) + + ######## DO WE FIT? ######## + + # here we assume that there have been enough new points to justify refitting the phase curve + # ie: that triage has already been done in adler-night + # whether we refit depends on some complex considerations (n_obs, phase angle range, minimum phase angle) + # however: could do a basic check: five new data points? + + ######## FIT NEW PHASE CURVE ######## + + if previousAdlerData: + params = planetoid.PreviousAdlerData.get_phase_parameters_in_filter(filt, phase_model) + H = params.H + G = params.phase_parameter_1 + else: + sso = planetoid.SSObject_in_filter(filt) + H = sso.H + G = sso.G12 + + if np.isnan(H) or np.isnan(G): + print( + "No H and/or phase parameter value found for this object and filter in SSObject or Adler databases." + ) + print("Using guess values. Fit may take longer.") + H = 15.0 + G = 0.5 + + pc = PhaseCurve( + H=H * u.mag, + phase_parameter_1=G, + model_name=phase_model, + ) + + pc_fit = pc.FitModel( + np.array(df_obs["phaseAngle"]) * u.deg, + np.array(df_obs["reduced_mag"]) * u.mag, + np.array(df_obs["magErr"]) * u.mag, + ) + + pc_fit = pc.InitModelSbpy(pc_fit) + pc_fit.__dict__["H"] = pc_fit.__dict__["H"].value + + planetoid.AdlerData.populate_phase_parameters(filt, **pc_fit.__dict__) + + ######## OUTLIER DETECTION: ONE METHOD ######## + + # can do phase curve fit including all but last point for outlier detection + # use H and G from alert if first time, otherwise pull from AdlerData + + # check for that phase curve model what the predicted values are, check alert point against those within threshold + # can later compare to LSST-predicted mag + + # of course this only checks the most recent point - what if the outlier was ten alerts ago but we didn't know? + # what if it was an alert we triaged out? + + # write new row to Adler database + # planetoid.AdlerData.write_row_to_database("adler-day-testing.db") + + +if __name__ == "__main__": + + # an assumption that these are the ssObjectIds of interest + # in the future, there will be a function which reads from the Kafka topic of triaged alerts from adler-night + # and returns this list. + # why are some of them negative? Ken doesn't know + ssobject_list = [ + "685835657114313227", + "-6411717652251119994", + "3881603647193383121", + "6711448832195578320", + "2959207197656523004", + "-8724852777827460529", + "-7394109909111663666", + "2154070653759973672", + "-7805034407650925178", + "-4490542680786866656", + ] + + run_adler_day(ssobject_list) diff --git a/src/adler/adler.py b/src/adler/adler_run.py similarity index 100% rename from src/adler/adler.py rename to src/adler/adler_run.py diff --git a/tests/data/adler-day-testing.db b/tests/data/adler-day-testing.db new file mode 100644 index 0000000000000000000000000000000000000000..91e515e2c83b15582ab0824a852af1e658c8be5d GIT binary patch literal 16384 zcmeI12~ZSQ8pme@0dWK|Dq=ha-6(Pl{rYY(eoYiK;^;mTk6^NR3Eksc?f(s~k-&bNZ*@l^M!%+-VWz}r9x{ImqesAWZfB$!U z@Be!IpGStpXyh{^wR6-lazfHkB9%&bxm+TV_z53t;bS}o2p7MCRgheze(W)atG<9VWk<+8=&T4IN-D&tRgsaACA`q<_-uPan(P|r&)&CDw zlv=BvqlwXIRmfC*D`Ma#V&M81+9rSvN(_Kn1E5w9sBHw;@ZK2scOoRS;* z0BkXa0WhoyV1p6^V3+|gOb-~=2(aP3{^1&wH39IHw_M+f7y@96F${q5O#mB|7y#o9 zfbn|3_(p&Y@AVJYpsX3dmg`#)LjY_sMm^www!+*ab^N7n2C~+l{wR}?_?)5wq5`4< zq5`4|WS) zwk2)Gwf)38)he&eS?OKjEC1(*PU+-9jFL#~B(}m)2iwvWjwljiFq$G5isKl9W^f$C zNiT$XVXPcs11U5R5ekyQ5yN_5*gy>9(GMkCPb5r)ROzVN@0WLj+Wz)CzYN?C$5;Ms za%Zej zP?BEj9|=dE77ac1R|m972eu37V#xBP3{yWJv=dR`Srz4l*IuePJG>=>_HedcEw+&4wL!um5#b z_UAn2nY_8Tb`;D=2&npW)nll3EKI+9-WS&WeacZA|4`ohS-Fe%4Qu%1M5ilA_j2` zP0MmPAVN6Ja>h=2?sM6noQJ@w**(H8eLoCpCv_Qt=0AXgTOIZuyZ)F*W2Sdc zEX##?y)sY9eR!zcwaq)_vkQ<#t>4pj_P0F3&h;&7^AbKwyD3XozJaRKR8;Aa`>;K8 zd}Q{T2Ru?-?wPbh4$H8NkUsObKrQj=;b0Y92e9+alwjp8+Sl)$#;P>UI#$Q@QUP8 zK7dDeF1l(wD9F1v@X6`SGPsz(aq76F8Ys=(8CzM9qEpJu03=uz(HKFKLL(9cOCTII z2gp(^Nl?Z?*zW7aGlOGbd$$2zPiJPsb+@boN&kq33$Cfl-o)N@HrIjJv^<8 ztIgZ^Pmhi+Ut3hEn`ZBUT0N1{LgNs^@T=>p1R`*pX|4q+gb5R@fse3Smp^JP4qrWX zmGu~Z32xq37Z=^jhW-9$9WGW_@pxjtwe5mVKuDV{H!D_kh3achhL-JJ4__yCYwk?DJ%fb$Emy}9XCNv0$139BOQ^gy;l7Jw z7uY?1xyNpAADyb5IW5Sr7=;m}kPc}IS#C~|C`611FLmPJ2M8M$y~1oWH>@t4;TO~p6hMRYC==FUQbBV9AXIrp~Z8(jX8uS zbG@FBq*;P_n<*@v>)V>sNloQ?Js?5i!VF~Wq?XS0)@A^j%=LOe9H(()T#s5j*IStZ zXfoI90TG6?jBy^tES>Aym@^1X=K4B9Ou-2YZo-%gEuQP8?=aU3d>E%-af(7lJ`3mi zI{)9^CQ9PiOYr|wW%JvwulN7Y1^d9wx=+LrCGrtQ;jzv$$=E95H?3J!nkZ#I4u6QpW8IZ!dGmVv28z`9ku9D@Gl?2 zscUh|5B*xjdlg(ebTnih?3!EMt+-29II}t~c?7Wy&OIEw_2;ZK-s{TA^wRE2;8==N z_sh?}gX*hyt7F@%;nwkG9e-UDuT#Bm&YH)hJErOhFcx8?k-*~axVst4R+H|y4voWh zj1;CDV<)w=JML-*ph5lb; z!u|%18Qr7B-SKidfgiOy-#ul`&JJzW>nbU%%bAA0wVTxuM;_VAX%kG5#2SI`; APXGV_ literal 0 HcmV?d00001