From 470734964c82dd1e35a3dc3d22b164ca7ab57172 Mon Sep 17 00:00:00 2001 From: Katie Harrington Date: Mon, 15 Jul 2024 20:19:55 +0000 Subject: [PATCH 1/2] add new source calculator --- src/pages/2_SATP1_v2.py | 330 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 330 insertions(+) create mode 100644 src/pages/2_SATP1_v2.py diff --git a/src/pages/2_SATP1_v2.py b/src/pages/2_SATP1_v2.py new file mode 100644 index 0000000..3b56f5a --- /dev/null +++ b/src/pages/2_SATP1_v2.py @@ -0,0 +1,330 @@ +import datetime as dt +from functools import partial +import yaml +import os +import numpy as np +from matplotlib import pyplot as plt +import matplotlib.dates as mdates + +import streamlit as st +from streamlit_timeline import st_timeline +from streamlit_ace import st_ace +from streamlit_sortables import sort_items + +from schedlib import policies, core, utils +from schedlib import rules as ru, source as src, instrument as inst +from scheduler_server.configs import get_config + +from schedlib.policies.satp1 import make_geometry +from schedlib.thirdparty import SunAvoidance + +from so3g.proj import quat, CelestialSightLine +from sotodlib import coords, core as todlib_core + +import jax.tree_util as tu + +""" +streamlit run src/Home.py --server.address=localhost --browser.gatherUsageStats=false --server.fileWatcherType=none --server.port 8075 +""" + +geometry = make_geometry() + + +basedir = '/so/home/kmharrin/software/scheduler-scripts/satp1' +schedule_files = { + 50 : os.path.join(basedir, 'master_files/cmb_2024_el50_20240423.txt'), + 60 : os.path.join(basedir, 'master_files/cmb_2024_el60_20240423.txt'), +} + +array_focus = { + 0 : { + 'left' : 'ws3,ws2', + 'middle' : 'ws0,ws1,ws4', + 'right' : 'ws5,ws6', + 'bottom': 'ws1,ws2,ws6', + 'all' : 'ws0,ws1,ws2,ws3,ws4,ws5,ws6', + }, + 45 : { + 'left' : 'ws3,ws4', + 'middle' : 'ws2,ws0,ws5', + 'right' : 'ws1,ws6', + 'bottom': 'ws1,ws2,ws3', + 'all' : 'ws0,ws1,ws2,ws3,ws4,ws5,ws6', + }, + -45 : { + 'left' : 'ws1,ws2', + 'middle' : 'ws6,ws0,ws3', + 'right' : 'ws4,ws5', + 'bottom': 'ws1,ws6,ws5', + 'all' : 'ws0,ws1,ws2,ws3,ws4,ws5,ws6', + }, +} + +SOURCES = ['Moon', 'Jupiter', 'Saturn'] + +def tod_from_block( block, ndet=100 ): + # pretty sure these are in degrees + t, az, alt = block.get_az_alt() + + tod = todlib_core.AxisManager( + todlib_core.LabelAxis('dets', ['a%02i' % i for i in range(ndet)]), + todlib_core.OffsetAxis('samps', len(t)) + ) + tod.wrap_new('timestamps', ('samps', ))[:] = t + + bs = todlib_core.AxisManager(tod.samps) + bs.wrap_new('az', ('samps', ))[:] = np.mod(az,360)* coords.DEG + bs.wrap_new('el', ('samps', ))[:] = az*0 + alt * coords.DEG + bs.wrap_new('roll', ('samps', ))[:] = az*0+block.boresight_angle*coords.DEG + tod.wrap('boresight', bs) + + return tod + +def get_focal_plane(tod): + xid_list = [] + etad_list = [] + + roll = np.mean(tod.boresight.roll) + + for waf in geometry: + xi0, eta0 = geometry[waf]['center'] + R = geometry[waf]['radius'] + phi = np.arange(tod.dets.count) * 2*np.pi / tod.dets.count + qwafer = quat.rotation_xieta(xi0 * coords.DEG, eta0 * coords.DEG) + qdets = quat.rotation_xieta(R * coords.DEG * np.cos(phi), + R * coords.DEG * np.sin(phi)) + if roll != 0: + q_bore_rot = quat.euler(2, -roll * coords.DEG) + qwafer = q_bore_rot * qwafer + qdets = q_bore_rot * qdets + + xid, etad, _ = quat.decompose_xieta(qwafer * qdets) + xid_list.append(xid) + etad_list.append(etad) + + xid = np.concatenate(xid_list) + etad = np.concatenate(etad_list) + return xid, etad + +now = dt.datetime.utcnow() +start_date = now.date() +start_time = now.time() +end_date = start_date + dt.timedelta(days=1) +end_time = start_time + +if 'timing' not in st.session_state: + st.session_state['timing'] = {} + st.session_state['timing']['t0'] = dt.datetime.combine( + start_date, start_time, tzinfo=dt.timezone.utc + ) + st.session_state['timing']['t1'] = dt.datetime.combine( + end_date, end_time, tzinfo=dt.timezone.utc + ) + st.session_state['timing']['sun_avoid'] = 45 + + +left_column, right_column = st.columns(2) + +with left_column: + start_date = st.date_input("Start date", value=start_date) + end_date = st.date_input("End date", value=end_date) + + sun_avoid_angle = st.number_input( + "Sun Avoidance Angle (deg):", + min_value= 0, + max_value= 90, + value=41, + step=1, + ) + + sun_avoid_time = st.number_input( + "Sun Avoidance Time (min)", + min_value= 0, + max_value= 60, + value=33, + step=1, + ) + +with right_column: + start_time = st.time_input("Start time (UTC)", value=start_time) + end_time = st.time_input("End time (UTC)", value=end_time) + +sources = st.multiselect("Sources", SOURCES, SOURCES) + +if st.button('Plot Sources'): + st.session_state['timing']['t0'] = dt.datetime.combine( + start_date, start_time, tzinfo=dt.timezone.utc + ) + st.session_state['timing']['t1'] = dt.datetime.combine( + end_date, end_time, tzinfo=dt.timezone.utc + ) + st.session_state['timing']['sun_avoid_angle'] = sun_avoid_angle + st.session_state['timing']['sun_avoid_time'] = sun_avoid_time + + t0=st.session_state['timing']['t0'] + t1=st.session_state['timing']['t1'] + sun_avoid_angle = st.session_state['timing']['sun_avoid_angle'] + sun_avoid_time = st.session_state['timing']['sun_avoid_time'] + + st.header("Source Availability") + st.write("Lighter line indicates source is cut by sun avoidance") + + fig = plt.figure(figsize=(8,3.75)) + ax = fig.add_subplot(111) + sun = SunAvoidance( + min_angle=sun_avoid_angle, + min_sun_time=sun_avoid_time*60 + ) + + for c, source in enumerate(sources): + src_blocks = src.source_gen_seq(source.lower(), t0, t1) + for block in src_blocks: + t, az, alt = block.get_az_alt(time_step=30) + plt.plot([dt.datetime.utcfromtimestamp(x) for x in t], alt, f'C{c}-', alpha=0.3) + + src_blocks = core.seq_flatten(sun.apply(src_blocks)) + + for b,block in enumerate(src_blocks): + if b == 0: + lab=source + else: + lab=None + t, az, alt = block.get_az_alt(time_step=30) + plt.plot([dt.datetime.utcfromtimestamp(x) for x in t], + alt, f'C{c}-', label=lab) + + locator = mdates.AutoDateLocator() + formatter = mdates.ConciseDateFormatter(locator) + + ax.xaxis.set_major_locator(locator) + ax.xaxis.set_major_formatter(formatter) + plt.xticks() + plt.ylim(0,90) + plt.legend() + + plt.xlabel("Azimuth (deg)") + plt.ylabel("Elevation (deg)") + ax.set_title( + f"{t0.strftime('%Y-%m-%d %H:%M')} to " + f"{t1.strftime('%Y-%m-%d %H:%M')}" + ) + st.pyplot(fig) + +with st.form("my data",clear_on_submit=False): + + st.title("Calibration Targets") + + col1, col2, col3 = st.columns(3) + with col1: + target = st.radio( + "Array Target", + ["all", "left", "middle", "right", "bottom", "custom"], + index=0, + ) + source = st.radio( + "Source to Scan", + SOURCES, + index=0, + ) + with col2: + ws0 = st.checkbox("ws0", value=False) + ws1 = st.checkbox("ws1", value=False) + ws2 = st.checkbox("ws2", value=False) + ws3 = st.checkbox("ws3", value=False) + + elevation = st.number_input( + "Elevation (deg)", + min_value = 48, + max_value = 80, + value=50, + step=1, + ) + boresight = st.number_input( + "boresight (deg)", + min_value = -60, + max_value = 60, + value=0, + step=1, + ) + + with col3: + ws4 = st.checkbox("ws4", value=False) + ws5 = st.checkbox("ws5", value=False) + ws6 = st.checkbox("ws6", value=False) + + min_scan_duration=st.number_input( + "Minimum Scan Duration (min)", + min_value = 0, + max_value = 60, + value=10, + step=1, + ) + + run_calculation = st.form_submit_button("Calculate") + if run_calculation: + if target == "custom": + arr = ['ws0', 'ws1', 'ws2', 'ws3', 'ws4', 'ws5', 'ws6'] + x = [ws0, ws1, ws2, ws3, ws4, ws5, ws6] + target_str = ','.join( [arr[i] for i in range(len(arr)) if x[i]]) + elif target == "all": + target_str = 'ws0,ws1,ws2,ws3,ws4,ws5,ws6' + else: + if int(boresight) not in [-45,0,45]: + st.write( + f":red[Boresight must be -45,0, or 45 to use array focus.]" + ) + target_str = array_focus[boresight][target] + + st.write(f"Target String is {target_str}") + + t0=st.session_state['timing']['t0'] + t1=st.session_state['timing']['t1'] + sun_avoid_angle = st.session_state['timing']['sun_avoid_angle'] + sun_avoid_time = st.session_state['timing']['sun_avoid_time'] + + sun = SunAvoidance( + min_angle=sun_avoid_angle, + min_sun_time=sun_avoid_time*60 + ) + min_dur_rule = ru.make_rule( + 'min-duration', **{'min_duration': min_scan_duration*60}, + ) + src_blocks = sun(src.source_gen_seq(source.lower(), t0, t1)) + #st.write(f"Source Blocks {src_blocks}") + + array_info = inst.array_info_from_query(geometry, target_str) + ces_rule = ru.MakeCESourceScan( + array_info=array_info, + el_bore=elevation, + drift=True, + boresight_rot=boresight, + allow_partial=True, + ) + scan_blocks = ces_rule(src_blocks) + st.write(f"Scan Blocks {scan_blocks}") + + scan_blocks = core.seq_flatten(min_dur_rule(sun(scan_blocks))) + st.write(f"Scan Blocks {scan_blocks}") + + for block in scan_blocks: + fig = plt.figure(figsize=(5,3.75)) + ax = fig.add_subplot(111) + + tod = tod_from_block(block) + xi_fp, eta_fp = get_focal_plane(tod) + ax.scatter(xi_fp, eta_fp, c='k', alpha=0.5) + + csl = CelestialSightLine.az_el( + tod.timestamps, tod.boresight.az, tod.boresight.el, weather='vacuum') + ra, dec, _ = quat.decompose_lonlat(csl.Q) + src_path = coords.planets.SlowSource.for_named_source( + source, tod.timestamps.mean() + ) + ra0, dec0 = src_path.ra, src_path.dec + ra0 = (ra0 - ra[0]) % (2 * np.pi) + ra[0] + # Un-rotate the planet into boresight coords. + xip, etap, _ = quat.decompose_xieta( + ~csl.Q * quat.rotation_lonlat(ra0, dec0) + ) + plt.plot(xip, etap, alpha=0.5) + st.pyplot(fig) \ No newline at end of file From b984a37c9e7aa2d0c43ce3995ec7841dfb7791fa Mon Sep 17 00:00:00 2001 From: Katie Harrington Date: Fri, 16 Aug 2024 18:45:37 +0000 Subject: [PATCH 2/2] add source scan planner --- src/pages/2_SATP1.py | 187 ------------------ ...{1_Sun_Avoidance.py => 2_Sun_Avoidance.py} | 0 ...{2_SATP1_v2.py => 3_SAT_Source_Planner.py} | 37 +++- 3 files changed, 30 insertions(+), 194 deletions(-) delete mode 100644 src/pages/2_SATP1.py rename src/pages/{1_Sun_Avoidance.py => 2_Sun_Avoidance.py} (100%) rename src/pages/{2_SATP1_v2.py => 3_SAT_Source_Planner.py} (89%) diff --git a/src/pages/2_SATP1.py b/src/pages/2_SATP1.py deleted file mode 100644 index 506b323..0000000 --- a/src/pages/2_SATP1.py +++ /dev/null @@ -1,187 +0,0 @@ -import datetime as dt -from functools import partial -import yaml -import pandas as pd - -import streamlit as st -from streamlit_timeline import st_timeline -from streamlit_ace import st_ace -from streamlit_sortables import sort_items - -from schedlib import policies, core, utils -from scheduler_server.configs import get_config - -import jax.tree_util as tu - -SOURCES = ['moon', 'mercury', 'venus', 'mars', 'jupiter', 'saturn', 'uranus', 'neptune'] - -# ==================== -# utility functions -# ==================== - -def seq2visdata(seqs): - # make group - def block2group(path, block): - key = utils.path2key(path) - if key == '': key = 'root' - return { - 'id': key, - 'content': key - } - groups = tu.tree_leaves( - tu.tree_map_with_path( - block2group, - seqs, - is_leaf=lambda x: isinstance(x, list) - ), - is_leaf=lambda x: 'id' in x, - ) - - # make items - def block2item(block, group=""): - res = { - 'id': hash(block), - 'content': block.name, - 'start': block.t0.isoformat(), - 'end': block.t1.isoformat(), - } - if group != "": res['group'] = group - return res - items = tu.tree_leaves( - tu.tree_map_with_path( - lambda path, x: core.seq_map( - partial(block2item, group=utils.path2key(path)), - core.seq_sort(x, flatten=True) - ), - seqs, - is_leaf=lambda x: isinstance(x, list) - ), - is_leaf=lambda x: 'id' in x - ) - return items, groups - -# ==================== -# initialize session state -# ==================== - -if 'user_config_str' not in st.session_state: - st.session_state.user_config_str = "{}" - -if 'commands' not in st.session_state: - st.session_state.commands = "" - -if 'checkpoints' not in st.session_state: - st.session_state.checkpoints = {} - -# ==================== -# sidebar UI -# ==================== - -with st.sidebar: - st.subheader("Schedule") - now = dt.datetime.utcnow() - - start_date = now.date() - start_time = now.time() - end_date = start_date + dt.timedelta(days=1) - end_time = start_time - start_date = st.date_input("Start date", value=start_date) - start_time = st.time_input("Start time (UTC)", value=start_time) - end_date = st.date_input("End date", value=end_date) - end_time = st.time_input("End time (UTC)", value=end_time) - - options = [] - for _src in SOURCES: - options += [ - [_src, 'left_boresight_0', 50, 0, 'left_focal_plane'], - [_src, 'middle_boresight_0', 50, 0, 'middle_focal_plane'], - [_src, 'right_boresight_0', 50, 0, 'right_focal_plane'], - [_src, 'bottom_boresight_0', 50, 0, 'bottom_focal_plane'], - [_src, 'left_boresight_p45', 50, 45, 'left_focal_plane'], - [_src, 'middle_boresight_p45', 50, 45, 'middle_focal_plane'], - [_src, 'right_boresight_p45', 50, 45, 'right_focal_plane'], - [_src, 'bottom_boresight_p45', 50, 45, 'bottom_focal_plane'], - [_src, 'left_boresight_n45', 50, -45, 'left_focal_plane'], - [_src, 'middle_boresight_n45', 50, -45, 'middle_focal_plane'], - [_src, 'right_boresight_n45', 50, -45, 'right_focal_plane'], - [_src, 'bottom_boresight_n45', 50, -45, 'bottom_focal_plane'] - ] - cal_targets_candidate = yaml.safe_load(st.session_state.user_config_str).get('cal_targets', []) - cal_targets = st.multiselect("Calibration Sources", options=options) - user_config = yaml.safe_load(st.session_state.user_config_str) - user_config['cal_targets'] = cal_targets - merge_order = list(set([tar[0] for tar in cal_targets])) + ['baseline'] - st.text("Priority: (descending) ") - merge_order_sorted = sort_items(merge_order) - - user_config = yaml.safe_load(st.session_state.user_config_str) - user_config['cal_targets'] = cal_targets - user_config['merge_order'] = merge_order_sorted - st.session_state.user_config_str = yaml.dump(user_config) - - with st.expander("Customize Source", expanded=False): - source_name = st.selectbox("Name", options=SOURCES) - elevation = st.number_input("Elevation (deg)", value=50.0) - boresight_angle = st.selectbox("Boresight angle", options=[0, 45, -45]) - query = st.multiselect("Array query", options=[ - 'left_boresight_0', 'middle_boresight_0', 'right_boresight_0', 'bottom_boresight_0', - 'left_boresight_p45', 'middle_boresight_p45', 'right_boresight_p45', 'bottom_boresight_p45', - 'left_boresight_n45', 'middle_boresight_n45', 'right_boresight_n45', 'bottom_boresight_n45', - 'ws0', 'ws1', 'ws2', 'ws3', 'ws4', 'ws5', 'ws6', - ]) - tag = st.text_input("Tag", value="") - def on_add(): - user_config = yaml.safe_load(st.session_state.user_config_str) - new_entry = [source_name, ",".join(query), elevation, boresight_angle, tag] - if 'cal_targets' not in user_config: user_config['sources'] = [] - user_config['cal_targets'].append(new_entry) - st.session_state.user_config_str = yaml.dump(user_config) - def on_reset(): - user_config = yaml.safe_load(st.session_state.user_config_str) - user_config['cal_targets'] = [] - st.session_state.user_config_str = yaml.dump(user_config) - st.button("Add source", on_click=on_add) - st.button("Reset sources", on_click=on_reset) - _sources = yaml.safe_load(st.session_state.user_config_str).get('cal_targets', []) - st.table(pd.DataFrame( - _sources, - columns=['source', 'query', 'elevation', 'boresight', 'tag'])) - - with st.expander("Advanced"): - # user_config = st.text_area("Config overwrite:", value=json.dumps(st.session_state.user_config, indent=2), height=300) - user_config_str = st_ace(value=st.session_state.user_config_str, language='yaml') - try: - user_config = yaml.safe_load(user_config_str) - # save a good config on parsing success - st.session_state.user_config_str = user_config_str - except Exception as e: - st.error('Unable to parse config', icon="🚨") - user_config = yaml.safe_load(st.session_state.user_config_str) - - def on_load_schedule(): - t0 = dt.datetime.combine(start_date, start_time).astimezone(dt.timezone.utc) - t1 = dt.datetime.combine(end_date, end_time).astimezone(dt.timezone.utc) - - config = get_config('satp1') - config = utils.nested_update(config, user_config) - policy = policies.SATPolicy.from_config(config) - - seqs = policy.apply(policy.init_seqs(t0, t1)) - commands = policy.seq2cmd(seqs, t0, t1) - st.session_state.checkpoints = policy.checkpoints - st.session_state.commands = commands - - st.button("Generate Schedule", on_click=on_load_schedule) - - -# ==================== -# main page -# ==================== - -for ckpt_name, ckpt_seqs in st.session_state.checkpoints.items(): - with st.expander(f"Checkpoint: {ckpt_name}", expanded=True): - data, groups = seq2visdata(ckpt_seqs) - timeline = st_timeline(data, groups, key=ckpt_name) - -with st.expander("Commands", expanded=False): - st.code(str(st.session_state.commands), language='python') \ No newline at end of file diff --git a/src/pages/1_Sun_Avoidance.py b/src/pages/2_Sun_Avoidance.py similarity index 100% rename from src/pages/1_Sun_Avoidance.py rename to src/pages/2_Sun_Avoidance.py diff --git a/src/pages/2_SATP1_v2.py b/src/pages/3_SAT_Source_Planner.py similarity index 89% rename from src/pages/2_SATP1_v2.py rename to src/pages/3_SAT_Source_Planner.py index 3b56f5a..55a6962 100644 --- a/src/pages/2_SATP1_v2.py +++ b/src/pages/3_SAT_Source_Planner.py @@ -23,9 +23,9 @@ import jax.tree_util as tu -""" +""" How to run this in your own directory streamlit run src/Home.py --server.address=localhost --browser.gatherUsageStats=false --server.fileWatcherType=none --server.port 8075 -""" +"""; geometry = make_geometry() @@ -80,6 +80,28 @@ def tod_from_block( block, ndet=100 ): return tod +def plot_focal_plane(ax, tod): + roll = np.mean(tod.boresight.roll) + + for waf in geometry: + xi0, eta0 = geometry[waf]['center'] + R = geometry[waf]['radius'] + phi = np.arange(tod.dets.count) * 2*np.pi / tod.dets.count + qwafer = quat.rotation_xieta(xi0 * coords.DEG, eta0 * coords.DEG) + qdets = quat.rotation_xieta(R * coords.DEG * np.cos(phi), + R * coords.DEG * np.sin(phi)) + if roll != 0: + q_bore_rot = quat.euler(2, -roll )#* coords.DEG) + qwafer = q_bore_rot * qwafer + qdets = q_bore_rot * qdets + + xi_c, eta_c, _ = quat.decompose_xieta(qwafer ) + xid, etad, _ = quat.decompose_xieta(qwafer * qdets) + + ax.scatter( xid, etad, marker='.') + ax.text( xi_c, eta_c, waf) + + def get_focal_plane(tod): xid_list = [] etad_list = [] @@ -120,7 +142,6 @@ def get_focal_plane(tod): st.session_state['timing']['t1'] = dt.datetime.combine( end_date, end_time, tzinfo=dt.timezone.utc ) - st.session_state['timing']['sun_avoid'] = 45 left_column, right_column = st.columns(2) @@ -202,7 +223,7 @@ def get_focal_plane(tod): plt.ylim(0,90) plt.legend() - plt.xlabel("Azimuth (deg)") + #plt.xlabel("Time") plt.ylabel("Elevation (deg)") ax.set_title( f"{t0.strftime('%Y-%m-%d %H:%M')} to " @@ -311,8 +332,9 @@ def get_focal_plane(tod): ax = fig.add_subplot(111) tod = tod_from_block(block) - xi_fp, eta_fp = get_focal_plane(tod) - ax.scatter(xi_fp, eta_fp, c='k', alpha=0.5) + #xi_fp, eta_fp = get_focal_plane(tod) + #ax.scatter(xi_fp, eta_fp, c='k', alpha=0.5) + plot_focal_plane(ax, tod) csl = CelestialSightLine.az_el( tod.timestamps, tod.boresight.az, tod.boresight.el, weather='vacuum') @@ -326,5 +348,6 @@ def get_focal_plane(tod): xip, etap, _ = quat.decompose_xieta( ~csl.Q * quat.rotation_lonlat(ra0, dec0) ) - plt.plot(xip, etap, alpha=0.5) + ax.plot(xip, etap, alpha=0.5) + ax.set_title( block.t0.isoformat() + f'\n{block.az} throw:{block.throw}' ) st.pyplot(fig) \ No newline at end of file