importrandomfromfunctoolsimportpartialimportpdbfromastropyimportunitsasufromastropy.timeimportTimefromssapy.constantsimportEARTH_MUfromssapy.utilsimportnorm,normSq,normedfromssapyimportrvsamplerimportnumpyasnpimportssapyfrom.utilsimportlazy_propertyfrom.constantsimportRGEOkeppropagator=ssapy.propagator.KeplerianPropagator()priorvolume=dict(rv=(2*40*10**6)**3*(2*7000)**3*1.0,# rough guess of total volume of prior for satellite# parameter space (2*RGEO)**3*(2*v_LEO)**3 / m^3 / (m/s)^3equinoctial=2*40*10**6*(4*np.pi)*np.pi*2*np.pi,# equivalent rough guess in equinoctial space# out to 2*RGEO in semimajor axis, 4pi for the# plane of the orbit, pi for ex and ey, 2*pi for lambda.angle=4*np.pi*np.pi*0.01**2*2*40*10**6*7000,# very rough guess for angular space# 4 pi for the position, 0.01 rad/s for the proper motions# 2*RGEO for the distance, 7000 km/s for the range in# allowed los velocities. Vaguely the right ballpark?)# dead tracks: want to prune hypotheses that differ only by dead tracks.# every time a track dies, do the N^2 search for equality among all of the# hypotheses it includes?
[docs]classCircVelocityPrior:"""Gaussian prior that log(GM/r/v^2) = 0. Parameters ---------- sigma : float Prior standard deviation on log(GM/r/v^2), dimensionless. Attributes ---------- sigma Methods ------- __call__(orbit) Returns log prior probability of given orbit. """def__init__(self,sigma=np.log(2)/3.5):self.sigma=sigma
[docs]def__call__(self,orbit,distance,chi=False):"""Return log prior probability of given orbit. Parameters ---------- orbit : Orbit The given orbit. chi : bool If True, return chi corresponding to log probability, rather than log probability. Returns ------- logprior : float The log of the prior probability for given orbit. """loggmorv2=np.log(EARTH_MU/norm(orbit.r)/normSq(orbit.v))chi0=loggmorv2/self.sigmaifchi:res=chi0else:res=-0.5*chi0**2return[res]
[docs]classZeroRadialVelocityPrior:"""Gaussian prior that v_R = 0. Parameters ---------- sigma : float Prior standard deviation on np.dot(v, r)/|v|/|r|, dimensionless. Attributes ---------- sigma Methods ------- __call__(orbit) Returns log prior probability of given orbit. """def__init__(self,sigma=0.25):self.sigma=sigma
[docs]def__call__(self,orbit,distance,chi=False):"""Return log prior probability of given orbit. Parameters ---------- orbit : Orbit The given orbit. chi : bool If True, return chi corresponding to log probability, rather than log probability. Returns ------- logprior : float The log of the prior probability for given orbit. """vdotr=np.dot(normed(orbit.r),normed(orbit.v))chi0=vdotr/self.sigmaifchi:res=chi0else:res=-0.5*chi0**2return[res]
[docs]classGaussPrior:"""Gaussian prior on parameters. Parameters ---------- mu : array_like (6) Mean of prior cinvcholfac : array_like (6, 6) Inverse Cholesky factor for covariance matrix, chi = cinvcholfac*(param-mu) translator : function(orbit) -> param function that translates from orbits to parameters Attributes ---------- mu, cinvcholfac, tranlator Methods ------- __call__(orbit) Returns log prior probability of given orbit. """def__init__(self,mu,cinvcholfac,translator):self.mu=muself.cinvcholfac=cinvcholfacself.translator=translator
[docs]def__call__(self,orbit,distance=None,chi=False):"""Return log prior probability of given orbit. Parameters ---------- orbit : Orbit The given orbit. chi : bool If True, return chi corresponding to log probability, rather than log probability. Returns ------- logprior : float The log of the prior probability for given orbit. """param=self.translator(orbit)chi0=self.cinvcholfac.dot(param-self.mu)ifchi:res=chi0else:res=-0.5*chi0**2returnres
[docs]classVolumeDistancePrior:"""Prior on range like r^2*exp(-r/scale), preferring larger distances where there is more volume. Methods ------- __call__(orbit) Returns log prior probability of given orbit. """def__init__(self,scale=RGEO):self.scale=RGEO
[docs]def__call__(self,orbit,distance,chi=False):"""Return log prior probability of given orbit. Parameters ---------- orbit : Orbit The given orbit. distance : float Distance between object and observer (m) chi : bool If True, return chi corresponding to log probability, rather than log probability. Returns ------- logprior : float The log of the prior probability for given orbit. """distance=distance/self.scalelnp=2*np.log(distance)-distance-(2*np.log(2)-2+1e-9)# - np.log(2*self.scale) --- true normalization factor.# 2*np.log(2) - 2 --- fake normalization factor that takes# max(lnp) = 0, so more chi^2 like.ifchi:res=np.sqrt(-2*lnp)else:res=lnpreturnres
[docs]deffit_arc_blind(arc,verbose=False,mode='rv',priors=None,propagator=None,damp=-1,orbitattr=None,optimizerkw={},lsq=True,factor=2,**kw):"""Fit an orbit to an arc. The arc can be arbitrarily long. Orbits are initially fit to the initial two observations and extended using a geomtrically increasing time baseline. If no new observations are inside of a given baseline, the next obs is added. The fit is least-squares, Levenberg-Marquardt, given specified priors and propagator. If the factor is set to 1, the orbit is built out iteratively by satID. Parameters ---------- arc : array_like (n) numpy ndarray containing several fields necessary for describing observations. rStation_GCRF (m), vStation_GCRF (m), time (astropy.time.Time or gps seconds), ra, dec, pmra, pmdec are all required fields. verbose : bool True for more verbose output mode : str Mode for fitting, must be one of 'rv', 'equinoctial', or 'angle'. In the first case, the parameters defining the orbit are taken to be the position and velocity of the object. In the second case, the parameters are taken to be the equinoctial elements. In the third case, the parameters are taken to be the angular position and proper motion of the object, together with its line-of-sight range and velocity. priors : list of objects implementing ssa Prior interface any priors to be applied in fitting propagator : ssapy.propagator.Propagator instance propagator to use in fitting damp : float damp chi residuals according to pseudo-Huber loss function. Default of -1 means do not damp orbitattr : list(str) names of additional orbit propagation keywords to guess optimizerkw : dict any extra keywords to pass to the optimizer factor : float factor for geometrically increasing time baseline (default 2) **kw : dict any extra keywords to pass to optimizer.optimize() """assertfactor>=1,"Geometric factor must be greater than or equal to 1"iflen(arc)==0:raiseValueError('len(arc) must be > 0')ifpriorsisNone:priors=[rvsampler.EquinoctialExEyPrior(0.3)]satids_ordered,times_ordered=time_ordered_satIDs(arc,with_time=True)paramInit,epoch=rvsampler.circular_guess(arc[np.isin(arc['satID'],satids_ordered[:2])])meanpm='pmra'inarc.dtype.namesparamInit=make_param_guess(paramInit,arc,mode=mode,orbitattr=orbitattr)nparam=6iforbitattrisNoneelse6+len(orbitattr)ifepoch.gps!=paramInit[nparam]:raiseValueError('inconsistent epoch and guess?')optimizerclass=make_optimizer(mode,paramInit,lsq=lsq)ifpropagatorisNone:propagator=keppropagatorgroupind=[0]whilegroupind[-1]<len(times_ordered)-1:dt=times_ordered[groupind[-1]]-times_ordered[0]newgroupind=np.max(np.flatnonzero(times_ordered<=times_ordered[0]+factor*dt))ifnewgroupind<=groupind[-1]:newgroupind=groupind[-1]+1groupind.append(newgroupind)forindingroupind:arc0=arc[np.isin(arc['satID'],satids_ordered[:ind+1])]prob=rvsampler.RVProbability(arc0,epoch,priors=priors,propagator=propagator,meanpm=meanpm,damp=damp)optimizer=optimizerclass(prob,paramInit[:nparam],orbitattr=orbitattr,**optimizerkw)state=optimizer.optimize(**kw)paramInit=[xforxinstate]+[xforxinparamInit[nparam:]]ifgetattr(optimizer.result,'hithyperbolicorbit',False):print('hyperbolic orbit in blind')chi2=np.sum(optimizer.result.residual**2)ifgetattr(optimizer.result,'hithyperbolicorbit',False):print('hyperbolic orbit')chi2+=1e9return(chi2,np.concatenate([state,paramInit[nparam:]]),optimizer.result)
[docs]deffit_arc(arc,guess,verbose=False,propagator=keppropagator,mode='rv',priors=None,damp=-1,optimizerkw={},lsq=True,orbitattr=None,**kw):"""Fit an orbit to an arc. See documentation for fit_arc_blind. fit_arc implements the same interface, but takes a additional guess & epoch arguments. These specify the initial guess in appropriate units given the chosen mode. epoch specifies the time at which this initial guess applies. Parameters ---------- arc : array_like (n) numpy ndarray containing several fields necessary for describing observations. rStation_GCRF (m), vStation_GCRF (m), time (astropy.time.Time or gps seconds), ra, dec, pmra, pmdec are all required fields. guess : array_like (n_par) parameters for initial guess guess[nparam] should be the epoch in GPS guess[nparam:] should be extra parameters not optimized over verbose : bool True for more verbose output propagator : ssapy.propagator.Propagator instance propagator to use in fitting mode : str Mode for fitting, must be one of 'rv', 'equinoctial', or 'angle'. In the first case, the parameters defining the orbit are taken to be the position and velocity of the object. In the second case, the parameters are taken to be the equinoctial elements. In the third case, the parameters are taken to be the angular position and proper motion of the object, together with its line-of-sight range and velocity. priors : list of objects implementing ssa Prior interface any priors to be applied in fitting damp : float Damping used in pseudo-Huber loss function. Default of -1 means no damping. orbitattr : list(str) Names of orbit attributes used for propagation (mass, area, ...) optimizerkw : dict any extra keywords to pass to the optimizer **kw : dict any extra keywords to pass to optimizer.optimize() """ifpriorsisNone:priors=[rvsampler.EquinoctialExEyPrior(0.3)]paramInit=guessmeanpm='pmra'inarc.dtype.namesnparam=6+(0iforbitattrisNoneelselen(orbitattr))epoch=Time(paramInit[nparam],format='gps')prob=rvsampler.RVProbability(arc,epoch,priors=priors,propagator=propagator,meanpm=meanpm,damp=damp)optimizerclass=make_optimizer(mode,paramInit,lsq=lsq)optimizer=optimizerclass(prob,paramInit[:nparam],orbitattr=orbitattr,**optimizerkw)state=optimizer.optimize(**kw)chi2=np.sum(optimizer.result.residual**2)ifgetattr(optimizer.result,'hithyperbolicorbit',False):chi2+=1e9return(chi2,np.concatenate([state,paramInit[nparam:]]),optimizer.result)
[docs]deffit_arc_with_gaussian_prior(arc,mu,cinvcholfac,verbose=False,propagator=keppropagator,mode='rv',lsq=True,optimizerkw={},orbitattr=None,**kw):"""Fit an orbit to an arc. See documentation for fit_arc_blind. fit_arc_with_gaussian_prior implements the same interface, but takes an additional mu, cinvcholfac, and epoch arguments. These specify a gaussian prior on the parameters. The mean and covariance of the Gaussian are given by mu and cinvcholfac, the inverse Cholesky matrix corresponding to the covariance matrix. The initial guess is taken to be the mean of the Gaussian prior. Presently, any other priors are assumed to have been folded into the Gaussian. epoch specifies the time at which Gaussian prior parameters are applicable. For greatest speed, the arc should contain a single observation and the epoch should be equal to the time of that observation; then no propagations are needed in this function. Parameters ---------- arc : array_like (n) numpy ndarray containing several fields necessary for describing observations. rStation_GCRF (m), vStation_GCRF (m), time (astropy.time.Time or gps seconds), ra, dec, pmra, pmdec are all required fields. mu : array_like (n_par) parameters for initial guess cinvcholfac : array_like (n_par, n_par) inverse cholesky matrix for covariance matrix verbose : bool True for more verbose output propagator : ssapy.propagator.Propagator instance propagator to use in fitting mode : str Mode for fitting, must be one of 'rv', 'equinoctial', or 'angle'. In the first case, the parameters defining the orbit are taken to be the position and velocity of the object. In the second case, the parameters are taken to be the equinoctial elements. In the third case, the parameters are taken to be the angular position and proper motion of the object, together with its line-of-sight range and velocity. orbitattr : list(str) list of strings of names of additional propagation attributes to fit (mass, area, cr, cd, ...) optimizerkw : dict any extra keywords to pass to the optimizer **kw : dict any extra keywords to pass to optimizer.optimize() """nparam=6+(0iforbitattrisNoneelselen(orbitattr))translator=partial(orbit_to_param,mode=mode,orbitattr=orbitattr,fitonly=True)ifmode=='angle':translator=partial(translator,rStation=mu[-6:-3],vStation=mu[-3:])priors=[GaussPrior(mu[:nparam],cinvcholfac,translator=translator)]paramInit=mumeanpm='pmra'inarc.dtype.namesepoch=Time(paramInit[nparam],format='gps')prob=rvsampler.RVProbability(arc,epoch,priors=priors,propagator=propagator,meanpm=meanpm)optimizerclass=make_optimizer(mode,paramInit,lsq=lsq)optimizer=optimizerclass(prob,paramInit[:nparam],orbitattr=orbitattr,**optimizerkw)state=optimizer.optimize(**kw)chi2=np.sum(optimizer.result.residual**2)ifgetattr(optimizer.result,'hithyperbolicorbit',False):chi2+=1e9return(chi2,np.concatenate([state,paramInit[nparam:]]),optimizer.result)
[docs]defdata_for_satellite(data,satIDList):"""Extract data for specific satID. This helper routine takes a set of observations containing a column 'satID', and returns the subset of those observations for which 'satID' is in the list satIDList. Parameters ---------- data : array_like (n), including column 'satID' data from which to select matching rows satIDList : list list of satIDs to select from data """m=np.zeros(len(data),dtype='bool')forsatIDinsatIDList:ifsatID==-1:continuem0=data['satID']==satIDifnp.sum(m0)==0:raiseValueError('satID %d not found!'%satID)m|=m0returndata[m]
[docs]defparam_to_orbit(param,mode='rv',orbitattr=None):"""Convert parameters to Orbits. Note: Orbits are converted so that they are all bound orbits to avoid challenges with hyperbolic orbit propagation. Parameters ---------- param : array_like (n, n_param) parameters for orbits mode : str parameter type, one of 'rv', 'angle', 'orbit' """param=np.atleast_1d(param)iflen(param.shape)==1:param=param[None,:]squeeze=Trueelse:squeeze=Falsenparam=6iforbitattrisnotNone:nparam+=len(orbitattr)else:orbitattr=[]tt=Time(param[:,nparam],format='gps')propkw=dict()forname,valinzip(orbitattr,param[:,6:nparam].T):ifname.startswith('log10'):name=name[5:]val=10.**valpropkw[name]=valifmode=='rv':orbit=ssapy.orbit.Orbit(r=param[:,:3],v=param[:,3:6],t=tt,propkw=propkw)elifmode=='equinoctial':aa=np.clip(param[:,0],1,np.inf)orbit=ssapy.orbit.Orbit.fromEquinoctialElements(aa,param[:,1],param[:,2],param[:,3],param[:,4],param[:,5],t=tt,propkw=propkw)elifmode=='angle':rv=ssapy.compute.radecRateObsToRV(param[:,0],param[:,1],param[:,2],param[:,3],param[:,4],param[:,5],param[:,-6:-3],param[:,-3:])orbit=ssapy.orbit.Orbit(r=rv[0],v=rv[1],t=tt,propkw=propkw)else:raiseValueError('unknown mode')vmax=np.sqrt(2*EARTH_MU/norm(orbit.r))m=norm(orbit.v)>vmaxifnp.any(m):vnew=orbit.v.copy()vnew[m,:]=vnew[m,:]*vmax[m,None]/norm(vnew[m,:])[:,None]*0.999orbit=ssapy.Orbit(r=orbit.r,v=vnew,t=orbit.t,propkw=propkw)ifsqueeze:orbit=orbit[0]returnorbit
[docs]deforbit_to_param(orbit,mode='rv',orbitattr=None,rStation=None,vStation=None,fitonly=False):"""Convert Orbits to parameters. Note: Orbits are converted so that they are all bound orbits to avoid challenges with hyperbolic orbit propagation. Parameters ---------- param : array_like (n, n_param) parameters for orbits mode : str parameter type, one of 'rv', 'angle', 'orbit' """iforbit.r.ndim==1:squeeze=Trueorbit=ssapy.Orbit(orbit.r[None,:],orbit.v[None,:],np.atleast_1d(orbit.t),propkw={k:np.atleast_1d(v)fork,vinorbit.propkw.items()})else:squeeze=Falsevmax=np.sqrt(2*EARTH_MU/norm(orbit.r))m=norm(orbit.v)>vmaxifnp.any(m):vnew=orbit.v.copy()vnew[m,:]=vnew[m,:]*vmax[m,None]/norm(vnew[m,:])[:,None]*0.999orbit=ssapy.Orbit(r=orbit.r,v=vnew,t=orbit.t,propkw=orbit.propkw)tt=orbit.t[:,None]iforbitattrisNone:orbitattr=list()orbitparam=[]forninorbitattr:ifn.startswith('log10'):val=np.log10(orbit.propkw[n[5:]])else:val=orbit.propkw[n]orbitparam.append(val[:,None])ifmode=='rv':param=np.hstack([orbit.r,orbit.v]+orbitparam+[tt])elifmode=='equinoctial':eqel=[e[:,None]foreinorbit.equinoctialElements]param=np.hstack([*eqel]+orbitparam+[tt])elifmode=='angle':rdrate=ssapy.compute.rvObsToRaDecRate(orbit.r,orbit.v,rStation,vStation)rdrate=[r[:,None]forrinrdrate]rStation=np.tile(rStation,len(orbit)).reshape(len(orbit),3)vStation=np.tile(vStation,len(orbit)).reshape(len(orbit),3)param=np.hstack([*rdrate]+orbitparam+[tt]+[rStation,vStation])else:raiseValueError('unknown mode')iffitonly:nparam=6+len(orbitattr)param=param[:,:nparam]ifsqueeze:param=param[0]returnparam
[docs]classTrackBase:""" Set of observations fit together as an orbital track. This class represents a set of observations that are being modeled as the motion of a single object through space. Parameters ---------- satIDs : array_like list of detection IDs considered part of this track data : array_like observations that could be part of this track. Observations that are actually part of this track are specified in satIDs. data is must contain a number of fields used in fitting, like ra, dec, obsPos, obsVel, and time. volume : float the volume in 6-D orbital element space in which objects could be found, needed to normalize false-positive vs. new satellite odds default to a very rough number only applicable for RV parameterization propagator : ssa Propagator instance the propagator to use when propagating the orbit. Default to None (Keplerian) Attributes ---------- satIDs : array_like, int the detection IDs that are part of this Track data : array_like observations that could be part of this track; see constructor volume : float the prior volume in 6-D parameter where satellites may be observed times : array_like, float the GPS times at which this object was observed mode : str the type of parameterization to use to fit the track one of 'rv', 'angle', 'equinoctial' lnprob : float the log probability of this association of observations to an orbit contains contributions both from priors and likelihood. """# new_track_lnpenalty = -13 # is this valuable?!def__init__(self,satIDs,data,volume=None,mode='rv',propagator=None,orbitattr=None):self.satIDs=list(satIDs)self.data=dataself.propagator=(propagatorifpropagatorisnotNoneelsessapy.propagator.KeplerianPropagator())ifvolumeisnotNone:self.volume=volumeelse:self.volume=priorvolume[mode]self.times=data_for_satellite(data,satIDs)['time'].gpsself.mode=modeself.orbitattr=orbitattr@lazy_propertydeflnprob(self):lnprob=-np.log(self.volume)# prior termifnp.any(~np.isfinite(self.covar)):lnprob+=0# i.e., we localized to 1 m^3 x 1 (m/s)^3...# for ~70 obs, localization is usually ~(100 m)^3*(0.01 m/s)^3# (at least formally!)# we still need to think more about treating the fact that the# errors should be more like ~40" in position# maybe translating from 3 position measurements to a# position measurement (with 40" uncertainties) and a velocity# measurement (with formal uncertainties only?)else:# SVD seems more robust than direct np.linalg.det_,ss,_=np.linalg.svd(2*np.pi*self.covar)determinant=np.product(ss)ifdeterminant<=0:determinant=1ifself.chi2<1e8:print('Weird determinant!',self.satIDs)# pdb.set_trace()lnprob+=0.5*np.log(determinant)# rarely, the fit screws up, and gets something with a huge covariance# in all cases I've encountered, due to hitting the hyperbolic orbit# boundary. We adjust the covariance volume here to be the same as# the prior volume, and add a huge number to the chi2 in# fit_arc.iflnprob>0:lnprob=-np.log(self.volume)ifself.chi2<1e8:print('warning, huge covariance',self.satIDs)# normalizing factor# if len(self.satIDs) == 1:# self.initial_lnprob = lnprob# lnprob -= self.initial_lnprob# lnprob += TrackBase.new_track_lnpenaltylnprob+=-self.chi2/2returnlnprob
[docs]defpredict(self,arc0,return_sigma=False,return_nwrap=True):"""Predict the location and uncertainty of this object at future times. Parameters ---------- arc0 : array_like, structure array must contain fields: time: future times (GPS) (float) rStation_GCRF: GCRF position of observer at future time (3, float) vStation_GCRF: GCRF velocity of observer at future time (3, float) return_sigma : bool if True, also return the sigma points used to compute the future covariance return_nwrap : bool if True, also return the number of orbital wrappings of the sigma points. Returns ------- mean, covar, fsigma mean : array_like (4, n) ra, dec, dra/dt*cos(dec), ddec/dt at future times (rad, rad/s) covar : array_like (4, 4, n) covariance of mean parameters at future times fsigma : array_like (4, 13, n) 13 sigma points used to compute covar fsigma[:, 0, :] is identical to mean. """nparam=6ifself.orbitattrisnotNone:nparam+=len(self.orbitattr)fixed_dimensions=np.arange(len(self.param))>=nparamfsigma=ssapy.utils.sigma_points(partial(self.propagaterdz,arc0=arc0,return_nwrap=return_nwrap),self.param,self.covar,fixed_dimensions=fixed_dimensions)mean=fsigma[:,0,...]dmean=fsigma[:,1:,...]-mean[:,None,...]dmean[0,...]=((dmean[0,...]+np.pi)%(2*np.pi))-np.piiflen(dmean.shape)==2:covar=np.cov(dmean,ddof=0)else:covar=np.array([np.cov(dmean[:,:,i],ddof=0)foriinrange(dmean.shape[-1])])# next step here would be to do some tangent# plane projections (stereoscopic?) so that we don't have# issues at the pole or with wraps in ra.res=mean,covarifreturn_sigma:res=res+(fsigma,)returnres
[docs]defgate(self,arc,return_nwrap=False):"""Compute chi^2 that this track could be associated with arc. Parameters ---------- arc : array_like, structure array must contain fields: time: future time (GPS) (float) rStation_GCRF: GCRF position of observer at future time (3, float) vStation_GCRF: GCRF velocity of observer at future time (3, float) ra: observed ra (float) pmra: observed proper motion in ra (float) dec: observed dec (float) pmdec: observed proper motion in dec (float) dra: uncertainty in ra (float) ddec: uncertainty in dec (float) dpmra: uncertainty in proper motion in ra (float) dpmdec uncertainty in proper motion in dec (float) return_nwrap : bool if True, also return the standard deviation of the number of orbits the object could have made at the time of observation. Returns ------- chi2, nwrapsig chi2 : float The chi-squared implied by associating this observation with this track. nwrapsig : float if return_nwrap, the standard deviation of the number of orbits the object could have made at the time of observation """arc0=arc[0]# just take the first tracklet at the moment...mean,covar=self.predict(arc0,return_nwrap=True)mean=mean[:-1]# drop nwrapcovar,nwrapcovar=covar[:-1,:-1],covar[-1,-1]cd=np.cos(arc0['dec'].to(u.rad).value)covar[:,0]*=cdcovar[0,:]*=cdcovarobs=np.diag([arc0['dra'].to(u.rad).value**2,arc0['ddec'].to(u.rad).value**2,arc0['dpmra'].to(u.rad/u.s).value**2,arc0['dpmdec'].to(u.rad/u.s).value**2])covar=covar+covarobsarea=np.linalg.det(2*np.pi*covar[0:2,0:2])**0.5if(area<=0)or(area>np.pi/8):covar=np.diag([np.inf,np.inf,np.inf,np.inf])try:cinv=np.linalg.inv(covar)except:cinv=1e-30measvec=np.array([arc0['ra'].to(u.rad).value,arc0['dec'].to(u.rad).value,arc0['pmra'].to(u.rad/u.s).value,arc0['pmdec'].to(u.rad/u.s).value])dvec=measvec-meandvec[0]=(((dvec[0]+np.pi)%(2*np.pi))-np.pi)*cdchi2=np.dot(dvec,cinv).dot(dvec)# we should choose a space other than ra/dec, where things# are coming closer to following lines, and so are more Gaussian.# note however: satellites are not following great circles,# so even in gnomonic projection about central point, our# points won't quite follow a straight line.# but it should at least be much better than whatever we're doing now!# or stereoscopic, with NP centered at the observation? Great circles# through projection point are lines.# simpler: line up track orbital plane with projection equator.res=chi2ifreturn_nwrap:res=(res,np.sqrt(nwrapcovar))returnres
[docs]defpropagaterdz(self,param,arc0=None,return_nwrap=False):"""Propagate param to times and locations specified in arc0. Parameters ---------- param : array_like parameters corresponding to self.mode for orbit to propagate arc0 : arary_like, must contain observation fields array containing fields specifying times and observer locations to which track fit should be propagated. return_nwrap : bool if True, additionally return nwrap Returns ------- array([rr, pmrr, dd, pmdd, nwrap]) rr : ra at future times, radians dd : dec at future times, radians pmrr : proper motion in ra at future times, rad / s pmdd : proper motion in dec at future times, rad / s nwrap : number of orbits completed. Only provided if return_nwrap is True. """orb=param_to_orbit(param,mode=self.mode,orbitattr=self.orbitattr)rr,dd,zz,pmrr,pmdd,dzzdt,nwrap=radeczn(orb,arc0,propagator=self.propagator)res=[rr,dd,pmrr,pmdd]ifreturn_nwrap:res=res+[nwrap]returnnp.array(res)
[docs]defupdate(self,time,rStation=None,vStation=None):"""Shift this track to a new time. Only implemented for numerical propagators where time-shifting can provide a big speed up. """return
[docs]classTrack(TrackBase):"""Set of observations to be fit as an object moving through space. Subclasses TrackBase, which could have other implementations (for example, a subclass that doesn't fit the whole track simultaneously, but instead approximates past information with a single Gaussian prior on the parameters). Parameters ---------- satIDs : see TrackBase data : see TrackBase guess : array_like, float guess of initial parameters fitting orbit initial_lnprob : float adjust probability of this track by initial_lnprob (default: 0) priors : rvsampler.Prior instance list of priors to use when fitting this track mode : str fitting mode to use. One of 'rv', 'angle', or 'equinoctial' parameters have different meanings in these modes propagator : Propagator instance ssa Propagator instance to use for propagation. Default None (Keplerian) """def__init__(self,satIDs,data,guess=None,initial_lnprob=0,priors=None,mode='rv',propagator=None,orbitattr=None):super().__init__(satIDs,data,mode=mode,propagator=propagator,orbitattr=orbitattr)self.initial_lnprob=initial_lnprobifguessisnotNone:chi2,param,res=fit_arc(data_for_satellite(data,self.satIDs),guess,maxfev=100,mode=mode,priors=priors,propagator=self.propagator,orbitattr=orbitattr)else:chi2,param,res=fit_arc_blind(data_for_satellite(data,self.satIDs),maxfev=100,mode=mode,priors=priors,propagator=self.propagator,orbitattr=orbitattr)self.chi2=chi2self.param=paramself.covar=getattr(res,'covar',np.array([[np.inf]]))self.priors=priorsself.success=getattr(res,'success',False)
[docs]defgaussian_approximation(self,propagator=None):"""Get Gaussian approximation of this track. Returns ------- TrackGauss, giving Gaussian approximation of this track. """# determinant = np.sqrt(np.linalg.det(2*np.pi*self.covar))# this is supposed to be a volume# when it's small enough, or when it's "close enough" to actually being# Gaussian, we want to make a Gaussian approximation.# but we really want, e.g., a fractional uncertainty in a.# I don't really know how to do this.# let's just say that we need at least four observations.# can do some tests with real tracks.iflen(self.satIDs)>=4:prop=propagatorifpropagatorisnotNoneelseself.propagatorifprop!=self.propagator:# refitting with a better propagatort=Track(self.satIDs,self.data,guess=self.param,initial_lnprob=self.initial_lnprob,priors=self.priors,mode=self.mode,propagator=prop,orbitattr=self.orbitattr)self.param=t.paramself.covar=t.covarself.chi2=t.chi2returnTrackGauss(self.satIDs,self.data,self.param,self.covar,self.chi2,mode=self.mode,propagator=prop,orbitattr=self.orbitattr)else:# don't bother approximating if the track is too short.returnself
# maybe we could just always approximate the track?# since the initial track is "on epoch", i.e., not propagated,# it's going to be some ~cone in xyz, vxvyvz.# that's not _great_ to approximate as a Gaussian, but maybe not# horrible?
[docs]defaddto(self,satid):"""Add a detection to this track. Parameters ---------- satID : int or list(int) the detection ID of the observation to add to this track Returns ------- A new Track, including the additional observation. """newsatidlist=satidifisinstance(satid,list)else[satid]returnTrack(self.satIDs+newsatidlist,self.data,guess=self.param,initial_lnprob=self.initial_lnprob,priors=self.priors,mode=self.mode,propagator=self.propagator,orbitattr=self.orbitattr)
[docs]classTrackGauss(TrackBase):"""Set of observations to be fit as an object moving through space. Subclasses TrackBase. This implementation does not fit the entire track simultaneously, but instead approximates past tracklets with a single Gaussian prior on the parameters. This information is then updated with each new observation in a Kalman-filter approach. Parameters ---------- satIDs : see TrackBase data : see TrackBase param : array_like, float mean of Gaussian covar : array_like, float covariance of Gaussian chi2 : float chi2 of fit at center of Gaussian mode : str fitting mode to use. One of 'rv', 'angle', or 'equinoctial' parameters have different meanings in these modes propagator : Propagator instance ssa Propagator instance to use for propagation. Default None (Keplerian) priors : list of Priors instances priors used in initially fitting this orbit. Their influence should now be completely recorded in param + covar, but they are recorded here for reference. """def__init__(self,satIDs,data,param,covar,chi2,mode='rv',propagator=None,sigma_points=None,priors=None,orbitattr=None):super().__init__(satIDs,data,mode=mode,propagator=propagator,orbitattr=orbitattr)self.param=paramself.covar=covarself.chi2=chi2ifnp.any(~np.isfinite(covar)):finitecovar=np.diag(np.ones(covar.shape[0],dtype='f4'))self.chi2+=1e9else:finitecovar=self.covarself.cinvcholfac=np.linalg.inv(np.linalg.cholesky(finitecovar))nparam=6iforbitattrisNoneelse6+len(orbitattr)ifsigma_pointsisNone:fixed=np.arange(len(param))>=nparamsigma_points=ssapy.utils.sigma_points(None,param,finitecovar,fixed_dimensions=fixed)self.sigma_points=sigma_pointsself.priors=priors
[docs]defgaussian_approximation(self,propagator):"""Get Gaussian approximation of this track. Returns ------- self. No op; this track already is a Gaussian approximation. """ifpropagator!=self.propagator:raiseNotImplementedErrorreturnself
[docs]defupdate(self,t,rStation=None,vStation=None):"""Move this track to a different epoch. Does not bother updating if t is different from the current track time by less than one microsecond. Parameters ---------- t : astropy.time.Time epoch for new TrackGauss rStation : array_like (3) position of station at new time (m, GCRF); required if mode == 'angle' vStation : array_like (3) velocity of station at new time (m, GCRF); required if mode == 'angle' """nparam=6ifself.orbitattrisNoneelse6+len(self.orbitattr)oldepoch=Time(self.param[nparam],format='gps')ifnp.abs((t-oldepoch).to(u.s).value)<1e-6:# don't bother propagating to a new time for times less than one# microsecond.returnorb=param_to_orbit(self.sigma_points,mode=self.mode,orbitattr=self.orbitattr)neworb=orb.at(t,propagator=self.propagator)ifisinstance(rStation,u.Quantity):rStation=rStation.to(u.m).valuevStation=vStation.to(u.m/u.s).valuenew_sigma=orbit_to_param(neworb,mode=self.mode,rStation=rStation,vStation=vStation,orbitattr=self.orbitattr)new_param=new_sigma[0]new_covar=np.cov((new_sigma[1:,:nparam]-new_param[:nparam]).T,ddof=0)self.param=new_paramself.covar=new_covarself.cinvcholfac=np.linalg.inv(np.linalg.cholesky(self.covar))self.sigma_points=new_sigma
[docs]defat(self,t,rStation=None,vStation=None):"""Get Gaussian approximation of this track at different epoch. This just copies the track and then updates it to be at the new time. Parameters ---------- t : astropy.time.Time epoch for new TrackGauss rStation : array_like (3) position of station at new time (m, GCRF); required if mode == 'angle' vStation : array_like (3) velocity of station at new time (m, GCRF); required if mode == 'angle' Returns ------- Gaussian-approximated track (TrackGauss) at new time. """newtrack=TrackGauss(self.satIDs,self.data,self.param,self.covar,self.chi2,mode=self.mode,propagator=self.propagator,sigma_points=self.sigma_points,orbitattr=self.orbitattr)newtrack.update(t,rStation=rStation,vStation=vStation)returnnewtrack
[docs]defaddto(self,satid):"""Add a detection to this track. This also shifts the track's epoch to the epoch of the added observation. Parameters ---------- satID : int the detection ID of the observation to add to this track Returns ------- A new TrackGauss, including the additional observation. """newsatidlist=satidifisinstance(satid,list)else[satid]arc=data_for_satellite(self.data,newsatidlist)nparam=6ifself.orbitattrisNoneelse6+len(self.orbitattr)mu=self.param[:nparam]oldepoch=Time(self.param[nparam],format='gps')ind=np.flatnonzero(arc['satID']==newsatidlist[0])lastind=ind[np.argmax(arc['time'][ind].gps)]newepoch=arc['time'][lastind]ifnewepoch.gps!=oldepoch.gps:newtrack=self.at(newepoch,rStation=arc['rStation_GCRF'][lastind],vStation=arc['vStation_GCRF'][lastind])param=newtrack.paramcinvcholfac=newtrack.cinvcholfacelse:param=self.paramcinvcholfac=self.cinvcholfacchi2new,paramnew,resnew=fit_arc_with_gaussian_prior(arc,param,cinvcholfac,mode=self.mode,propagator=self.propagator,orbitattr=self.orbitattr)covarnew=getattr(resnew,'covar',np.array([[np.inf]]))ifnp.any(~np.isfinite(covarnew)):covarnew=np.array([[np.inf]])returnTrackGauss(self.satIDs+newsatidlist,self.data,paramnew,covarnew,self.chi2+chi2new,mode=self.mode,propagator=self.propagator,orbitattr=self.orbitattr)
[docs]deffit_arc_blind_via_track(arc,propagator=None,verbose=False,reset_if_too_uncertain=False,mode='rv',priors=None,order='forward',approximate=False,orbitattr=None,factor=1):"""Fit a series of detections using the Track interface. Parameters ---------- arc : ndarry with many fields The observations to fit. Must contain many fields (ra, dec, ...) propagator : Propagator instance Propagator to use for orbit propagation. Default None (Keplerian). verbose : bool print additional information during fitting if True reset_if_too_uncertain : bool if the uncertainty in the number of wraps grows to be greater than pi, start track fitting over, ignoring old tracklets. factor : float factor for geometrically increasing time baseline (default 1) order : str order to process observations; default 'forward'. 'backward' processes observations in reverse chronological order. Returns ------- tracklist A list of Tracks, giving the fits to the track as each observation in added sequentially in time to the fit. """assertfactor>=1,"Geometric factor must be greater than or equal to 1"satids_ordered,times_ordered=time_ordered_satIDs(arc,with_time=True,order=order)tracklist=[]prop=Noneifapproximateelsepropagatortrack0=Track([satids_ordered[0]],arc,propagator=prop,mode=mode,priors=priors,orbitattr=orbitattr)tracklist.append(track0)resetnwrap=np.pi/3groupind=[0]whilegroupind[-1]<len(times_ordered)-1:dt=times_ordered[groupind[-1]]-times_ordered[0]newgroupind=np.max(np.flatnonzero(times_ordered<=times_ordered[0]+factor*dt))ifnewgroupind<=groupind[-1]:newgroupind=groupind[-1]+1groupind.append(newgroupind)foriinrange(1,len(groupind)):arc0=arc[np.isin(arc['satID'],satids_ordered[groupind[i]])]chi2,nwrapsig=tracklist[-1].gate(arc0[0:1],return_nwrap=True)reset=((reset_if_too_uncertainand(nwrapsig>resetnwrap))or(tracklist[-1].chi2>=1e9))ifverbose:ifreset:resetstr='resetting'else:resetstr=''dt=(times_ordered[groupind[i]]-times_ordered[groupind[i-1]])print(f'{chi2:5.1f}{nwrapsig:8.4f} dt:{dt/60/60/24:5.2f}{resetstr}')ifreset:track=Track(arc0['satID'],arc,propagator=prop,mode=mode,priors=priors,orbitattr=orbitattr)else:newsatidlist=satids_ordered[groupind[i-1]+1:groupind[i]+1]track=tracklist[-1].addto(newsatidlist)ifapproximate:track=track.gaussian_approximation(propagator)tracklist.append(track)returntracklist
[docs]defcombinatoric_lnprior(nsat,ntrack,ndet):"""Prior probability of assigning detections to tracks of satellites. This function intends to give the prior probability that ndet detections of nsat total satellites correspond to ntrack satellites. Parameters ---------- nsat : int number of satellites that could possibly be observed ntrack : int number of tracks in current solution ndet : int number of detections of these satellites """# I believe the desired factor is:# 1/N_sat^(N_detections - 1) * (N_sat - 1)! / (N_sat - N_objects)!fromscipyimportspeciallnprior=-(ndet-1)*np.log(nsat)lnprior+=special.gammaln(nsat-1)-special.gammaln(nsat-ntrack)# may have to think about an approximation to N_sat!/(N_sat - N_objects)!# that doesn't rely on as much cancellation as the above...returnlnprior
[docs]deftime_ordered_satIDs(data,with_time=False,order='forward'):"""Give satIDs in data in order of observation. Optionally includes first times for each ID. Parameters ---------- data : array_like the data to consider with_time : bool Option to also output list of corresponding gps times order : 'forward' or 'backward' Returns ------- list of satIDs in data, ordered by observation time list of corresponding gps times of first obs per ID """s=np.argsort(data['time'].gps)iforder=='backward':s=s[::-1]usedtracklet={}out=[]times=[]forsatidindata['satID'][s]:ifusedtracklet.get(satid,False):continueusedtracklet[satid]=Trueout.append(satid)ifwith_time:times.append(data[data['satID']==satid]['time'].gps[0])res=outifwith_time:res=(res,)+(times,)returnres
[docs]classHypothesis:"""Assignment of detections to multiple tracks. This class represents an assignment of detections to multiple tracks. The goal of Multiple Hypothesis Tracking to to find the Hypothesis (track assignment) that has the highest likelihood. Parameters ---------- tracks : list of Tracks The tracks in the hypothesis. Each track must correspond to different observations, and all observations under consideration must be accounted for Attributes ---------- tracks : list the tracks in this hypothesis lnprob : float the log probability of this hypothesis nsat : int the total number of satellites in the sky that could be observed """def__init__(self,tracks,nsat=1000):self.tracks=tracksself.lnprob=np.sum([t.lnprobfortintracks])ndet=self.ntracklet()self.lnprob+=combinatoric_lnprior(nsat,len(tracks),ndet)# not clear this has value.self.nsat=nsatdef__repr__(self):returnself.summarize()
[docs]defsummarize(self,verbose=False):"""Summarize the Hypothesis as a string. Parameters ---------- verbose : bool Include information about each track in the hypothesis. Returns ------- A string summarizing the hypothesis. """out='Hypothesis with %d tracks:'%len(self.tracks)out+=', lnprob: %7.4e'%self.lnprobifverbose:out+='\n'fortrackinself.tracks:out+=str(track)+'\n'out+='\n'returnout
[docs]defdifference(self,hypothesis):"""Print the difference between two hypotheses. Often two hypotheses are very similar. This makes it easier to inspect only the differences between them. Parameters ---------- hypothesis : Hypothesis hypothesis to which to compare this hypothesis """matchinother=np.zeros(len(hypothesis.tracks),dtype='bool')matchinsame=np.zeros(len(self.tracks),dtype='bool')fori,tinenumerate(self.tracks):try:ind=hypothesis.tracks.index(t)matchinsame[i]=Truematchinother[ind]=TrueexceptValueError:continueprint('dlnprob: %5.3f, %d tracks in common'%(self.lnprob-hypothesis.lnprob,sum(matchinsame)))print('Tracks only in 1:')print(self.summarize(verbose=False))foriinnp.flatnonzero(~matchinsame):print(self.tracks[i])print('Tracks only in 2:')print(hypothesis.summarize(verbose=False))foriinnp.flatnonzero(~matchinother):print(hypothesis.tracks[i])
[docs]defntracklet(self):"""Number of observations in tracks in the hypothesis."""n=0fortinself.tracks:n+=len(t.satIDs)returnn
[docs]@staticmethoddefaddto(hypothesis,track,oldtrack=None,**kw):"""Add a new track to a Hypothesis. Parameters ---------- hypothesis : Hypothesis hypothesis to which to add track : Track track to add to hypothesis oldtrack : Track track to remove from hypothesis, if the new track updates the old track. Returns ------- new Hypothesis, with track added, replacing oldtrack if oldtrack is not None. """newtracks=[tfortinhypothesis.tracks]ifoldtrackisnotNoneandoldtrack!=[]:oldtrackind=hypothesis.tracks.index(oldtrack)newtracks[oldtrackind]=trackelse:newtracks.append(track)returnHypothesis(newtracks,nsat=hypothesis.nsat,**kw)
[docs]classMHT:"""Multiple Hypothesis Tracking of many hypotheses explaining data. This class manages the assessment of which hypothesis are most likely to explain data, building a set of tracks observation by observation. Initialize the MHT with the data, do mht.run(), and inspect the most likely hypotheses using [h.lnprob for h in mht.hypotheses]. Parameters ---------- data : array_like[N] data to model. Must contain angles of observations, observer locations and velocities, times, detection IDs, etc. nsat : int total number of satellites in the sky. Affects the overall prior and the preference for new tracks vs. additional assignments to existing tracks. truth : dict when true assignments are known, this argument can provide debug information about when the tracker has lost the true assignment hypotheses : list of Hypothesis initialize the MHT with this existing list of Hypotheses. Default None. propagator : instance of ssa Propagator propagator to use. Default to None (Keplerian) fitonly : ndarray[N] bool entries in data to fit """def__init__(self,data,nsat=1000,truth=None,hypotheses=None,propagator=None,mode='rv',fitonly=None,approximate=False,orbitattr=None,priors=None):self.data=data.copy()iffitonlyisNone:fitonly=np.ones(len(data),dtype='bool')self.satids=time_ordered_satIDs(data[fitonly])self.nsat=nsatself.truth=truthself.track2hyp={}self.mode=modeself.nfit=0self.orbitattr=orbitattrself.priors=priorsself.propagator=(propagatorifpropagatorisnotNoneelsessapy.propagator.KeplerianPropagator())self._newly_dead_tracks=[]self.approximate=approximateifhypothesesisNone:self.hypotheses=[Hypothesis([],nsat=self.nsat)]else:forhypinhypotheses:fortinhyp.tracks:iftnotinself.track2hyp:self.track2hyp[t]=[]self.track2hyp[t].append(hyp)self.hypotheses=hypotheses
[docs]defrun(self,first=None,last=None,verbose=False,order='forward',**kw):"""Run the MHT analysis. Parameters ---------- first : int First tracklet to consider last : int Last tracklet to consider verbose : bool print extra progress information while running **kw : dict extra keyword arguments to pass to MHT.prune. """tsatid=self.satidsiffirstisNone:first=0iflastisNone:last=len(tsatid)tsatid=tsatid[first:last]iforder=='backward':tsatid=tsatid[::-1]fori,satidinenumerate(tsatid):ifverbose:ndead=sum([getattr(t,'dead',False)fortinself.track2hyp])print(('Tracklet %5d of %5d, %5d live tracks ''(%5d dead), %4d fit, %5d hypotheses')%(first+i,len(self.satids),len(self.track2hyp)-ndead,ndead,self.nfit,len(self.hypotheses)))self.add_tracklet(satid)self.prune(satid,**kw)
[docs]defadd_tracklet(self,satid,**kw):"""Add an observation from data to the MHT analysis. Parameters ---------- satid : int the ID of the observation to add. """newhypotheses=[]tracklist=list(self.track2hyp.keys())newtrack2hyp={t:[]fortintracklist}self.nfit=0tooshorttogate=0self._newly_dead_tracks=[]arc=data_for_satellite(self.data,[satid])if(len(arc)<1)and((('dpmra'notinarc.dtype.names)or(np.isfinite(arc['dpmra'][0])))):print('Just one observation, skipping tracklet %d'%satid)# used to match these, but not generate new tracks for these.# this is problematic since then hypotheses aren't directly# comparable to one another; they don't explain the same# numbers of tracksreturnfortrackintracklist:ifarc['time'].gps[0]intrack.times:continueifgetattr(track,'dead',False):continuetrack.update(arc['time'][0],rStation=arc['rStation_GCRF'][0],vStation=arc['vStation_GCRF'][0])gatechi2,nwrapsig=track.gate(arc,return_nwrap=True)outofgatechi2=8**2ifself.truthisnotNone:goodtrack=np.all([self.truth[s]==self.truth[satid]forsintrack.satIDs])ifgoodtrackandgatechi2>outofgatechi2:print('warning, excluding real track by gate',gatechi2)# pdb.set_trace()elifgoodtrack:# print('real gate chi2: %5.2f' % gatechi2)passifgatechi2>outofgatechi2:# if it's hard to add this tracklet to this track, don't# bother.continueifnwrapsig>np.pi/3:# at this point, there's real danger of fitting the wrong# number of orbits to this point. Don't do that; kill the# track instead, and hope we get better sampling of it later.# we could do better explicitly trying multiple sets of orbits# there's currently an assumption that a track is uniquely# defined by a set of satID; this would require opening another# dimension, or, better, playing some games with chi2 and covar# in the track, letting those try to track the multimodal# posteriors.track.dead=Trueself._newly_dead_tracks.append(track)continueifgatechi2<1e-9:tooshorttogate+=1# pdb.set_trace()ifgatechi2==0.:# pdb.set_trace()passnewtrack=track.addto(satid)self.nfit+=1ifnewtrack.chi2-track.chi2>outofgatechi2:continueifself.approximate:newtrack=newtrack.gaussian_approximation(self.propagator)if(gatechi2>1e-9)&(len(track.satIDs)>=2):# print(track.satIDs, '%5.1f %5.1f %5.1f %5.1f %5.1f' % (gatechi2, track.chi2, track.lnprob(),# newtrack.chi2, newtrack.lnprob()))# if (newtrack.chi2 > 1e9) & (len(track.satIDs) >= 3):# pdb.set_trace()passnewtrack2hyp[newtrack]=[]hypaffected=self.track2hyp[track]forhypothesisinhypaffected:newhypothesis=Hypothesis.addto(hypothesis,newtrack,track,**kw)newhypotheses.append(newhypothesis)fort0innewhypothesis.tracks:newtrack2hyp[t0].append(newhypothesis)prop=Noneifself.approximateelseself.propagatornewtrack=Track([satid],self.data,mode=self.mode,propagator=prop,orbitattr=self.orbitattr,priors=self.priors)forhinself.hypotheses:newhypothesis=Hypothesis.addto(h,newtrack,**kw)ifnewhypothesis.ntracklet()!=h.ntracklet()+1:pdb.set_trace()newhypotheses.append(newhypothesis)fortinnewhypothesis.tracks:iftnotinnewtrack2hyp:newtrack2hyp[t]=[]newtrack2hyp[t].append(newhypothesis)# if self.flag_inconsistency(newtrack2hyp, newhypotheses):# pdb.set_trace()self.hypotheses=newhypothesesself.track2hyp=newtrack2hyp
[docs]defprune_tracks(self,satid,nconfirm=6):"""Prune tracks that are now nearly identical from the MHT analysis. As the MHT adds observations to tracks, eventually some tracks have a number of confirming observations. For these long tracks, we really only need to keep different Tracks and Hypotheses for cases where there are recent disagreements with how those tracks continue; if hypotheses agree that the track continues in the same way, we don't need to continue tracking the past differences. So we prune out cases of former disagreement. This tries to identify such overlapping tracks and keep only one of them. Parameters ---------- satid : int the satID of the most-recently added detection to the MHT nconfirm : int if two tracks are identical in their last nconfirm observations delete the tracks not in the most likely hypothesis. Returns ------- boolean array indicating hypotheses to keep """lnprob=np.array([hyp.lnprobforhypinself.hypotheses])s=np.argsort(-lnprob)# best firstbesthyp=self.hypotheses[s[0]]lengthenedtrack=np.flatnonzero([satidintrack.satIDsfortrackinbesthyp.tracks])iflen(lengthenedtrack)!=1:pdb.set_trace()lengthenedtrack=besthyp.tracks[lengthenedtrack[0]]nsatid=len(lengthenedtrack.satIDs)keephyp=np.ones(len(self.hypotheses),dtype='bool')ifnsatid<=nconfirm+1:# no op; no sufficiently long tracks for this to be meaningful.returnkeephypkeeptogether=lengthenedtrack.satIDs[:nsatid-nconfirm]# relies on ordering of self.track2hyp being fixed...deltrack=[]ndel=0fori,trackinenumerate(self.track2hyp):ntogether=0forsidinkeeptogether:ntogether+=(sidintrack.satIDs)if(ntogether!=0)and(ntogether!=len(keeptogether)):deltrack.append(track)fortrackindeltrack:forhypinself.track2hyp[track]:hyp.deleteme=Truendel+=1keephyp=~np.array([getattr(hyp,'deleteme',False)forhypinself.hypotheses])ifndel>0:print('%d hypotheses pruned with differences more than %d frames ago.'%(ndel,nconfirm))# import pdb# pdb.set_trace()returnkeephyp
[docs]defprune_stale_hypotheses(self,newdeadtracks):"""Prune hypotheses that are different from better hypotheses only in dead tracks. Note: implentation relies on identical tracks always having the same id. The MHT code tries to guarantee this---otherwise it must do extra work matching identical tracks to new detections, etc. But using id(Track) feels horrible to me. Parameters ---------- newdeadtracks : list of Track list of tracks that died since the last pruning. Returns ------- boolean mask indicating hypotheses to keep """iflen(newdeadtracks)==0:returnnp.ones(len(self.hypotheses),dtype='bool')hyps=list(set(sum([self.track2hyp[t]fortinnewdeadtracks],[])))ntrack=[len(h.tracks)forhinhyps]maxntrack=max(ntrack)idarr=np.zeros((maxntrack,len(hyps)),dtype='i4')-1fori,hypinenumerate(hyps):idarr[:ntrack[i],i]=[id(t)ifnotgetattr(t,'dead',False)else-1fortinhyp.tracks]idarr=np.sort(idarr,axis=0)lnprob=[h.lnprobforhinhyps]s=np.lexsort([lnprob]+[idarr[i]foriinrange(len(idarr))])todelete=0foriinrange(1,len(hyps)):ifnp.all(idarr[:,s[i]]==idarr[:,s[i-1]]):hyps[s[i-1]].deleteme=Truetodelete+=1keephyp=~np.array([getattr(hyp,'deleteme',False)forhypinself.hypotheses])iftodelete>0:print('%d hypotheses pruned, differing only by dead tracks.'%todelete)returnkeephyp
[docs]defprune(self,satid,nkeepmax=10000,pkeep=1e-9,keeponlytrue=False,nconfirm=6):"""Prune unlikely hypotheses from the MHT. Parameters ---------- satid : int the detection ID most recently added to the MHT nkeepmax : int keep no more than nkeepmax hypotheses pkeep : float keep no hypotheses more than pkeep times less likely than the most likely hypothesis keeponlytrue : bool keep only the single hypothesis known to be true for speed comparison purposes. nconfirm : int only keep one variant of tracks that agree in the last nconfirm detections """lnprob=np.array([hyp.lnprobforhypinself.hypotheses])s=np.argsort(-lnprob)# best firstkeep=np.zeros(len(self.hypotheses),dtype='bool')bestprob=lnprob[s[0]]# keep ones with high enough probability relative to bestkeep[lnprob>bestprob+np.log(pkeep)]=1# don't keep more than nkeepmaxkeep[s[nkeepmax:]]=0ifnconfirm>0:keep=keep&self.prune_tracks(satid,nconfirm=nconfirm)keep=keep&self.prune_stale_hypotheses(self._newly_dead_tracks)ifnp.sum(keep)==0:raiseValueError('should not be possible!')# keep at least nkeepmin# keep[s[:nkeepmin]] = 1ifself.truthisnotNone:tracksat=[[[self.truth[s]forsint.satIDs]fortinh.tracks]forhinself.hypotheses]nomixup=[[all(s==t[0]forsint)fortinh]forhintracksat]nomixup=np.array([all(h)forhinnomixup])split=np.array([len(set(t[0]fortinh))!=len(h)forhintracksat])truthind=np.argmax(~split&nomixup)truthisalive=nomixup&~split&keepifnp.any(nomixup&~split):print('truth: dlnprob: %5.3f, nhyp: %d'%(bestprob-lnprob[truthind],np.sum(lnprob>lnprob[truthind])))ifnp.sum(truthisalive)!=1:print('warning: true solution is no longer included.')ifkeeponlytrue:keep=keep&nomixup&~split# print(len(self.hypotheses), np.sum(keep))self.hypotheses=[self.hypotheses[i]foriinnp.flatnonzero(keep)]live={h:Trueforhinself.hypotheses}fortrackinself.track2hyp:self.track2hyp[track]=([hyp0forhyp0inself.track2hyp[track]iflive.get(hyp0,False)])fortrackinlist(self.track2hyp):iflen(self.track2hyp[track])==0:delself.track2hyp[track]
# if self.flag_inconsistency(self.track2hyp, self.hypotheses):# pdb.set_trace()
[docs]@staticmethoddefflag_inconsistency(track2hyp,hyp):"""Debug method checking to see if there are any inconsistencies between the tracks tracked by the MHT and the hypotheses tracked by the MHT. Should only be needed for debugging. Every hypothesis in track2hyp[track] should include the Track track. Every hypothesis in hyp should be in the list of track2hyp[track] for each of its tracks. Every track should be in a hypothesis. Parameters ---------- track2hyp : dict[Track -> list(Hypothesis)] hyp: list(Hypothesis) Returns ------- True if everything is consistent. """consistent=0fortintrack2hyp:forhintrack2hyp[t]:iftnotinh.tracks:consistent|=1iflen(track2hyp[t])==0:consistent|=8forhinhyp:fortinh.tracks:ifhnotintrack2hyp[t]:consistent|=2ntracklet=np.array([sum([len(t.satIDs)fortinh.tracks],0)forhinhyp])iflen(ntracklet)>0andnp.any(ntracklet[0]!=ntracklet):consistent|=4ifconsistent!=0:pdb.set_trace()returnconsistent
[docs]defsummarize_tracklets(data,posuncfloor=None,pmuncfloor=None):"""Add fields to default data set to accommodate MHT fitting. Parameters ---------- data : array_like data to add fields to posuncfloor : float or None add a positional uncertainty floor to the ra/dec uncertainties (deg) pmuncfloor : float or None add a proper motion uncertainty floor to the ra/dec proper motion uncertainties (deg/s) Returns ------- array parallel to data with additional fields dra, ddec, pmra, pmdec, dpmra, dpmdec """datao=data.copy()data=data.copy()s=np.lexsort([data['time'],data['satID']])usat,first=np.unique(data['satID'][s],return_index=True)last=np.concatenate([first[1:],[len(data)]])newfields=[('dra',u.deg,None),('ddec',u.deg,None),('pmra',u.deg/u.s,None),('pmdec',u.deg/u.s,None),('dpmra',u.deg/u.s,None),('dpmdec',u.deg/u.s,None),('t_baseline',u.s,None)]forfield,unit,diminnewfields:dtype='f8'ifunitisnotNoneelse'object'ifdimisNone:data[field]=np.zeros(len(data),dtype=dtype)*unitelse:data[field]=np.zeros((len(data),dim),dtype=dtype)*unitforf,linzip(first,last):meanpos,dmeanpos,pm,dpm=summarize_tracklet(data[s[f:l]])r0=np.mean(data['rStation_GCRF'][s[f:l],:],axis=0)dr0=data['rStation_GCRF'][s[l-1]]-data['rStation_GCRF'][s[f]]dt0=data['time'][s[l-1]]-data['time'][s[f]]v0=(dr0/dt0).to(u.m/u.s)data['rStation_GCRF'][s[f:l]]=r0ifnp.any(~np.isfinite(data['vStation_GCRF'][s[f:l]])):data['vStation_GCRF'][s[f:l]]=v0# this isn't quite the same as taking the mean velocity over# the observations. But on Aero-B it makes a maximum difference# of less than one part in a million, which seems negligible# for _velocity_ measurements. Note that for _positions_ this# would be significant!else:data['vStation_GCRF'][s[f:l]]=(np.mean(data['vStation_GCRF'][s[f:l],:],axis=0)[None,:])data['time'][s[f:l]]=Time(np.mean(data['time'][s[f:l]].gps),format='gps')data['ra'][s[f:l]]=meanpos[0]data['dec'][s[f:l]]=meanpos[1]data['dra'][s[f:l]]=dmeanpos[0]data['ddec'][s[f:l]]=dmeanpos[1]data['pmra'][s[f:l]]=pm[0]data['pmdec'][s[f:l]]=pm[1]data['dpmra'][s[f:l]]=dpm[0]data['dpmdec'][s[f:l]]=dpm[1]data['t_baseline'][s[f:]]=dt0.to(u.s)ifposuncfloorisnotNone:data['dra']=np.sqrt(posuncfloor**2+data['dra']**2)data['ddec']=np.sqrt(posuncfloor**2+data['ddec']**2)ifpmuncfloorisnotNone:data['dpmra']=np.sqrt(pmuncfloor**2+data['dpmra']**2)data['dpmdec']=np.sqrt(pmuncfloor**2+data['dpmdec']**2)data=data[s[first]]s=np.argsort(data['time'].gps)returndata[s]
[docs]defsummarize_tracklet(arc):"""Compute mean ra/dec, uncertainty in mean, proper motion, uncertainty for a short arc. Parameters ---------- arc : the short arc of detections Returns ------- [(ra, dec), (dra, ddec), (pmra, pmdec), (dpmra, dpmdec)] ra, dec: the mean position in right ascension and declination dra, ddec: the uncertainty in the mean RA and declination pmra, pmdec: the proper motion in right ascension and declination dpmra, dpmdec: the uncertainty in the proper motion in RA and declination """unit=ssapy.utils.lb_to_unit(arc['ra'].to(u.rad).value,arc['dec'].to(u.rad).value)meanunit=np.mean(unit,axis=0)# defines tangent planepole=np.array([0,0,1])raunit=normed(np.cross(pole,meanunit))decunit=normed(np.cross(meanunit,raunit))ratp=np.einsum('ij,j',unit,raunit)dectp=np.einsum('ij,j',unit,decunit)iflen(unit)<=1:return((arc['ra'],arc['dec']),(arc['sigma'],arc['sigma']),(0.,0.),(np.inf,np.inf))meantime=np.mean(arc['time'].gps)dt=arc['time'].gps-meantimedpos=arc['sigma'].to(u.rad).valuerap,racovar=np.polyfit(dt,ratp,1,w=1./dpos,cov='unscaled')rasig=np.sqrt(np.diag(racovar))decp,deccovar=np.polyfit(dt,dectp,1,w=1./dpos,cov='unscaled')decsig=np.sqrt(np.diag(deccovar))meanra,meandec=ssapy.utils.unit_to_lb(meanunit[None,:])rap[1]+=meanra[0]decp[1]+=meandec[0]out=[(rap[1]*u.rad,decp[1]*u.rad),(rasig[1]*u.rad,decsig[1]*u.rad),(rap[0]*u.rad/u.s,decp[0]*u.rad/u.s),(rasig[0]*u.rad/u.s,decsig[0]*u.rad/u.s)]returnout
# need to write code that ~iterates MHT to try to bring in detections that weren't# matched.defiterate_mht(data,oldmht,nminlength=20,trimends=2,**kw):bestind=np.argmax([h.lnprobforhinoldmht.hypotheses])besthyp=oldmht.hypotheses[bestind]tracks=[tfortinbesthyp.tracksifnotgetattr(t,'dead',False)and(len(t.satIDs)>nminlength)]# want to cut off the ends of the tracksiftrimends>0:newtracks=[]fori,tinenumerate(tracks):newt=Track(t.satIDs[trimends:-trimends],t.data,guess=t.param,priors=t.priors,mode=t.mode,propagator=t.propagator,orbitattr=t.orbitattr)print(i,len(tracks),newt.chi2/4/len(newt.satIDs))ifoldmht.approximate:newt=newt.gaussian_approximation()newtracks.append(newt)tracks=newtracksinitialhyp=Hypothesis(tracks,nsat=besthyp.nsat)usedsatid=set(sum([t.satIDsfortintracks],[]))keepdata=np.array([snotinusedsatidforsindata['satID']])newmht=MHT(data,hypotheses=[initialhyp],nsat=oldmht.nsat,propagator=oldmht.propagator,mode=oldmht.mode,fitonly=keepdata,approximate=oldmht.approximate,priors=oldmht.priors,truth=oldmht.truth,orbitattr=oldmht.orbitattr)newmht.run(**kw)returnnewmht