-
Notifications
You must be signed in to change notification settings - Fork 1
/
demo_netflow_regions.py
178 lines (148 loc) · 6.57 KB
/
demo_netflow_regions.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
import ets_fiber_assigner.netflow as nf
import numpy as np
from ics.cobraOps.Bench import Bench
from ics.cobraOps.TargetGroup import TargetGroup
from ics.cobraOps.CobrasCalibrationProduct import CobrasCalibrationProduct
from ics.cobraOps.CollisionSimulator import CollisionSimulator
from ics.cobraOps.cobraConstants import NULL_TARGET_POSITION, NULL_TARGET_ID
from ics.cobraOps import plotUtils
from collections import defaultdict
# make runs reproducible
np.random.seed(20)
# define locations of the input files
catalog_path = "data/"
fsky_pos = catalog_path+"pfs_preliminary_target_cosmology_sky.dat"
tgt = nf.readCalibrationFromFile(fsky_pos, "sky")
# get a complete, idealized focal plane configuration
bench = Bench(layout="full")
# if you have the XML file, you can also generate a more realistic focal plane
# bench = Bench(calibrationProduct=CobrasCalibrationProduct(
# "../ics_cobraOps/python/ics/demos/updatedMaps6.xml"))
# point the telescope at the center of all science targets
raTel, decTel = nf.telescopeRaDecFromFile(fsky_pos)
posang = 0.
otime = "2016-04-03T08:00:00Z"
telescopes = []
# number of distinct observations
nvisit = 6
# generate randomly jittered telescope pointings for every observation
for _ in range(nvisit):
telescopes.append(nf.Telescope(raTel+np.random.normal()*1e-2,
decTel+np.random.normal()*1e-2, posang, otime))
# get focal plane positions for all targets and all visits
tpos = [tele.get_fp_positions(tgt) for tele in telescopes]
# create the dictionary containing the costs and constraints for all classes
# of targets
classdict = {}
classdict["sky"] = {"numRequired": 240,
"nonObservationCost": 1e6, "calib": True}
# optional: slightly increase the cost for later observations,
# to observe as early as possible
vis_cost = [i*10. for i in range(nvisit)]
# optional: penalize assignments where the cobra has to move far out
def cobraMoveCost(dist):
return 0.1*dist
# duration of one observation in seconds
t_obs = 300.
gurobiOptions = dict(seed=0, presolve=1, method=4, degenmoves=0,
heuristics=0.8, mipfocus=0, mipgap=1.0e-04)
# let's pretend that most targets have already been completely observed,
# and that the rest has been partially observed
alreadyObserved={}
for t in tgt:
alreadyObserved[t.ID] = 3
for t in tgt[::10]:
alreadyObserved[t.ID] = 1
forbiddenPairs = []
for i in range(nvisit):
forbiddenPairs.append([])
# define two location regions arbitrarily
ncobras = bench.cobras.nCobras
cobraRegions=np.zeros(ncobras,dtype=np.int32)
cobraRegions[ncobras//2:] = 1
done = False
while not done:
# compute observation strategy
prob = nf.buildProblem(bench, tgt, tpos, classdict, t_obs,
vis_cost, cobraMoveCost=cobraMoveCost,
collision_distance=2., elbow_collisions=True,
gurobi=False, gurobiOptions=gurobiOptions,
alreadyObserved=alreadyObserved,
forbiddenPairs=forbiddenPairs,
cobraLocationGroup=cobraRegions,
minSkyTargetsPerLocation=40,
locationGroupPenalty=1e9)
# print("writing problem to file ", mpsName)
# prob.dump(mpsName)
print("solving the problem")
prob.solve()
# extract solution
res = [{} for _ in range(nvisit)]
for k1, v1 in prob._vardict.items():
if k1.startswith("Tv_Cv_"):
visited = prob.value(v1) > 0
if visited:
_, _, tidx, cidx, ivis = k1.split("_")
res[int(ivis)][int(tidx)] = int(cidx)
print("Checking for trajectory collisions")
ncoll = 0
for ivis, (vis, tp) in enumerate(zip(res, tpos)):
selectedTargets = np.full(len(bench.cobras.centers), NULL_TARGET_POSITION)
ids = np.full(len(bench.cobras.centers), NULL_TARGET_ID)
for tidx, cidx in vis.items():
selectedTargets[cidx] = tp[tidx]
ids[cidx] = ""
for i in range(selectedTargets.size):
if selectedTargets[i] != NULL_TARGET_POSITION:
dist = np.abs(selectedTargets[i]-bench.cobras.centers[i])
simulator = CollisionSimulator(bench, TargetGroup(selectedTargets, ids))
simulator.run()
if np.any(simulator.endPointCollisions):
print("ERROR: detected end point collision, which should be impossible")
coll_tidx = []
for tidx, cidx in vis.items():
if simulator.collisions[cidx]:
coll_tidx.append(tidx)
ncoll += len(coll_tidx)
for i1 in range(0,len(coll_tidx)):
for i2 in range(i1+1,len(coll_tidx)):
if np.abs(tp[coll_tidx[i1]]-tp[coll_tidx[i2]])<10:
forbiddenPairs[ivis].append((coll_tidx[i1],coll_tidx[i2]))
print("trajectory collisions found:", ncoll)
done = ncoll == 0
# write output file
with open("output.txt", "w") as f:
for i, (vis, tp, tel) in enumerate(zip(res, tpos, telescopes)):
print("exposure {}:".format(i))
print(" assigned Cobras: {}".format(len(vis)))
tdict = defaultdict(int)
f.write("# Exposure {}: duration {}s, RA: {}, Dec: {}, PA: {}\n".
format(i+1, t_obs, tel._ra, tel._dec, tel._posang))
f.write("# Target Fiber X Y "
"RA DEC\n")
for tidx, cidx in vis.items():
tdict[tgt[tidx].targetclass] += 1
f.write("{:} {:6d} {:10.5f} {:10.5f} {:10.5f} {:10.5f}\n"
.format(tgt[tidx].ID, cidx+1, tp[tidx].real, tp[tidx].imag,
tgt[tidx].ra, tgt[tidx].dec))
for cls, num in tdict.items():
print(" {}: {}".format(cls, num))
for vis, tp in zip(res, tpos):
selectedTargets = np.full(len(bench.cobras.centers), NULL_TARGET_POSITION)
ids = np.full(len(bench.cobras.centers), NULL_TARGET_ID)
for tidx, cidx in vis.items():
selectedTargets[cidx] = tp[tidx]
ids[cidx] = ""
for i in range(selectedTargets.size):
if selectedTargets[i] != NULL_TARGET_POSITION:
dist = np.abs(selectedTargets[i]-bench.cobras.centers[i])
simulator = CollisionSimulator(bench, TargetGroup(selectedTargets, ids))
simulator.run()
simulator.plotResults(paintFootprints=False)
plotUtils.pauseExecution()
# Animate the trajectory collisions
(problematicCobras,) = np.where(np.logical_and(
simulator.collisions, ~simulator.endPointCollisions))
for cbr in problematicCobras:
simulator.animateCobraTrajectory(cbr)
plotUtils.pauseExecution()