yaflpy.pyx 56 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643
  1. # -*- coding: utf-8 -*-
  2. """
  3. Copyright 2021 anonimous <shkolnick-kun@gmail.com> and contributors.
  4. Licensed under the Apache License, Version 2.0 (the "License");
  5. you may not use this file except in compliance with the License.
  6. You may obtain a copy of the License at
  7. http://www.apache.org/licenses/LICENSE-2.0
  8. Unless required by applicable law or agreed to in writing,
  9. software distributed under the License is distributed on an "AS IS" BASIS,
  10. WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  11. See the License for the specific language governing permissions
  12. and limitations under the License.
  13. """
  14. #==============================================================================
  15. # YAFL C API
  16. #==============================================================================
  17. #cython: language_level=3
  18. #distutils: language=c
  19. from libc cimport stdint
  20. #------------------------------------------------------------------------------
  21. IF YAFLPY_USE_64_BIT:
  22. cdef extern from "yafl_config.h":
  23. ctypedef double yaflFloat
  24. ctypedef stdint.int32_t yaflInt
  25. ELSE:
  26. cdef extern from "yafl_config.h":
  27. ctypedef float yaflFloat
  28. ctypedef stdint.int32_t yaflInt
  29. #------------------------------------------------------------------------------
  30. cdef extern from "yafl_math.c":
  31. ctypedef enum yaflStatusEn:
  32. #Warning flag masks
  33. YAFL_ST_MSK_REGULARIZED = 0x01 #YAFL_ST_R
  34. YAFL_ST_MSK_GLITCH_SMALL = 0x02 #YAFL_ST_S
  35. YAFL_ST_MSK_GLITCH_LARGE = 0x04 #YAFL_ST_L
  36. YAFL_ST_MSK_ANOMALY = 0x08 #YAFL_ST_A
  37. #Everthing is OK
  38. YAFL_ST_OK = 0x00
  39. #Wagnings
  40. YAFL_ST_R = 0x01
  41. YAFL_ST_S = 0x02
  42. YAFL_ST_SR = 0x03
  43. YAFL_ST_L = 0x04
  44. YAFL_ST_LR = 0x05
  45. YAFL_ST_SL = 0x06
  46. YAFL_ST_SLR = 0x07
  47. YAFL_ST_A = 0x08
  48. YAFL_ST_AR = 0x09
  49. YAFL_ST_SA = 0x0a
  50. YAFL_ST_SAR = 0x0b
  51. YAFL_ST_LA = 0x0c
  52. YAFL_ST_LAR = 0x0d
  53. YAFL_ST_SLA = 0x0e
  54. YAFL_ST_SLAR = 0x0f
  55. # Error threshold value (greater values are errors)
  56. YAFL_ST_ERR_THR = 0x010
  57. # Invalid argument numer
  58. YAFL_ST_INV_ARG_1 = 0x100
  59. YAFL_ST_INV_ARG_2 = 0x110
  60. YAFL_ST_INV_ARG_3 = 0x120
  61. YAFL_ST_INV_ARG_4 = 0x130
  62. YAFL_ST_INV_ARG_5 = 0x140
  63. YAFL_ST_INV_ARG_6 = 0x150
  64. YAFL_ST_INV_ARG_7 = 0x160
  65. YAFL_ST_INV_ARG_8 = 0x170
  66. YAFL_ST_INV_ARG_9 = 0x180
  67. YAFL_ST_INV_ARG_10 = 0x190
  68. YAFL_ST_INV_ARG_11 = 0x1a0
  69. #--------------------------------------------------------------------------
  70. #cdef yaflStatusEn yafl_math_set_u(yaflInt sz, yaflFloat *res, yaflFloat *u)
  71. #------------------------------------------------------------------------------
  72. cdef extern from "yafl.c":
  73. ctypedef _yaflKalmanBaseSt yaflKalmanBaseSt
  74. ctypedef yaflStatusEn (* yaflKalmanFuncP)(yaflKalmanBaseSt *, \
  75. yaflFloat *, yaflFloat *)
  76. ctypedef yaflStatusEn (* yaflKalmanResFuncP)(yaflKalmanBaseSt *, \
  77. yaflFloat *, yaflFloat *, \
  78. yaflFloat *)
  79. ctypedef yaflStatusEn (* yaflKalmanScalarUpdateP)(yaflKalmanBaseSt *, \
  80. yaflInt)
  81. ctypedef yaflFloat (* yaflKalmanRobFuncP)(yaflKalmanBaseSt *, yaflFloat)
  82. ctypedef struct _yaflKalmanBaseSt:
  83. yaflKalmanFuncP f #
  84. yaflKalmanFuncP h #
  85. yaflKalmanResFuncP zrf #
  86. yaflFloat * x #
  87. yaflFloat * y #
  88. yaflFloat * Up #
  89. yaflFloat * Dp #
  90. yaflFloat * Uq #
  91. yaflFloat * Dq #
  92. yaflFloat * Ur #
  93. yaflFloat * Dr #
  94. yaflInt Nx #
  95. yaflInt Nz #
  96. #==========================================================================
  97. # UD-factorized EKF definitions
  98. #==========================================================================
  99. ctypedef struct yaflEKFBaseSt:
  100. yaflKalmanBaseSt base
  101. yaflKalmanFuncP jf #
  102. yaflKalmanFuncP jh #
  103. yaflFloat * H #
  104. yaflFloat * W #
  105. yaflFloat * D #
  106. cdef yaflStatusEn yafl_ekf_base_predict(yaflKalmanBaseSt * self)
  107. cdef yaflStatusEn \
  108. yafl_ekf_base_update(yaflKalmanBaseSt * self, yaflFloat * z, \
  109. yaflKalmanScalarUpdateP scalar_update)
  110. #--------------------------------------------------------------------------
  111. cdef yaflStatusEn \
  112. yafl_ekf_bierman_update_scalar(yaflKalmanBaseSt * self, yaflInt i)
  113. #--------------------------------------------------------------------------
  114. cdef yaflStatusEn \
  115. yafl_ekf_joseph_update_scalar(yaflKalmanBaseSt * self, yaflInt i)
  116. #==========================================================================
  117. ctypedef struct yaflEKFAdaptiveSt:
  118. yaflEKFBaseSt base
  119. yaflFloat chi2
  120. #--------------------------------------------------------------------------
  121. cdef yaflStatusEn \
  122. yafl_ekf_adaptive_bierman_update_scalar(yaflKalmanBaseSt * self, \
  123. yaflInt i)
  124. #--------------------------------------------------------------------------
  125. cdef yaflStatusEn \
  126. yafl_ekf_adaptive_joseph_update_scalar(yaflKalmanBaseSt * self, \
  127. yaflInt i)
  128. #--------------------------------------------------------------------------
  129. # For demonstration purposes only
  130. cdef yaflStatusEn \
  131. yafl_ekf_do_not_use_this_update_scalar(yaflKalmanBaseSt * self, yaflInt i)
  132. #==========================================================================
  133. ctypedef struct yaflEKFRobustSt:
  134. yaflEKFBaseSt base
  135. yaflKalmanRobFuncP g
  136. yaflKalmanRobFuncP gdot
  137. #--------------------------------------------------------------------------
  138. cdef yaflStatusEn \
  139. yafl_ekf_robust_bierman_update_scalar(yaflKalmanBaseSt * self, \
  140. yaflInt i)
  141. #--------------------------------------------------------------------------
  142. cdef yaflStatusEn \
  143. yafl_ekf_robust_joseph_update_scalar(yaflKalmanBaseSt * self, \
  144. yaflInt i)
  145. #==========================================================================
  146. ctypedef struct yaflEKFAdaptiveRobustSt:
  147. yaflEKFRobustSt base
  148. yaflFloat chi2
  149. #--------------------------------------------------------------------------
  150. cdef yaflStatusEn \
  151. yafl_ekf_adaptive_robust_bierman_update_scalar(yaflKalmanBaseSt * self, \
  152. yaflInt i)
  153. #--------------------------------------------------------------------------
  154. cdef yaflStatusEn \
  155. yafl_ekf_adaptive_robust_joseph_update_scalar(yaflKalmanBaseSt * self, \
  156. yaflInt i)
  157. #==========================================================================
  158. # UD-factorized UKF definitions
  159. #==========================================================================
  160. ctypedef _yaflUKFBaseSt yaflUKFBaseSt
  161. ctypedef yaflStatusEn (* yaflUKFSigmaAddP)(yaflUKFBaseSt *, yaflFloat *, \
  162. yaflFloat *, yaflFloat)
  163. ctypedef struct yaflUKFSigmaSt:
  164. yaflInt np
  165. yaflUKFSigmaAddP addf
  166. #--------------------------------------------------------------------------
  167. ctypedef yaflStatusEn (* yaflUKFSigmaGenWeigthsP)(yaflUKFBaseSt *)
  168. ctypedef yaflStatusEn (* yaflUKFSigmaGenSigmasP)(yaflUKFBaseSt *)
  169. ctypedef struct yaflUKFSigmaMethodsSt:
  170. yaflUKFSigmaGenWeigthsP wf
  171. yaflUKFSigmaGenSigmasP spgf
  172. #--------------------------------------------------------------------------
  173. ctypedef struct _yaflUKFBaseSt:
  174. yaflKalmanBaseSt base
  175. yaflUKFSigmaSt * sp_info
  176. const yaflUKFSigmaMethodsSt * sp_meth
  177. yaflKalmanFuncP xmf
  178. yaflKalmanResFuncP xrf
  179. yaflKalmanFuncP zmf
  180. yaflFloat * zp
  181. yaflFloat * Sx
  182. yaflFloat * Pzx
  183. yaflFloat * sigmas_x
  184. yaflFloat * sigmas_z
  185. yaflFloat * wm
  186. yaflFloat * wc
  187. #--------------------------------------------------------------------------
  188. cdef yaflStatusEn yafl_ukf_post_init(yaflUKFBaseSt * self) #static inline
  189. cdef yaflStatusEn yafl_ukf_gen_sigmas(yaflUKFBaseSt * self) #static inline
  190. cdef yaflStatusEn yafl_ukf_base_predict(yaflUKFBaseSt * self)
  191. cdef yaflStatusEn \
  192. yafl_ukf_base_update(yaflUKFBaseSt * self, yaflFloat * z, \
  193. yaflKalmanScalarUpdateP scalar_update)
  194. #==========================================================================
  195. cdef yaflStatusEn \
  196. yafl_ukf_bierman_update_scalar(yaflKalmanBaseSt * self, yaflInt i)
  197. #==========================================================================
  198. ctypedef struct yaflUKFAdaptivedSt:
  199. yaflUKFBaseSt base
  200. yaflFloat chi2
  201. cdef yaflStatusEn \
  202. yafl_ukf_adaptive_bierman_update_scalar(yaflKalmanBaseSt * self, \
  203. yaflInt i)
  204. #==========================================================================
  205. ctypedef struct yaflUKFRobustSt:
  206. yaflUKFBaseSt base
  207. yaflKalmanRobFuncP g
  208. yaflKalmanRobFuncP gdot
  209. #--------------------------------------------------------------------------
  210. cdef yaflStatusEn \
  211. yafl_ukf_robust_bierman_update_scalar(yaflKalmanBaseSt * self, \
  212. yaflInt i)
  213. #==========================================================================
  214. ctypedef struct yaflUKFAdaptiveRobustSt:
  215. yaflUKFRobustSt base
  216. yaflFloat chi2
  217. #--------------------------------------------------------------------------
  218. cdef yaflStatusEn \
  219. yafl_ukf_adaptive_robust_bierman_update_scalar(yaflKalmanBaseSt * self, \
  220. yaflInt i)
  221. #==========================================================================
  222. ctypedef struct yaflUKFSt:
  223. yaflUKFBaseSt base
  224. yaflFloat * Us
  225. yaflFloat * Ds
  226. #--------------------------------------------------------------------------
  227. cdef yaflStatusEn yafl_ukf_update(yaflUKFBaseSt * self, yaflFloat * z)
  228. #==========================================================================
  229. ctypedef struct yaflUKFFullAdapiveSt:
  230. yaflUKFSt base
  231. yaflFloat chi2
  232. yaflStatusEn yafl_ukf_adaptive_update(yaflUKFBaseSt * self, yaflFloat * z)
  233. #==========================================================================
  234. # Van der Merwe sigma point generator
  235. #==========================================================================
  236. ctypedef struct yaflUKFMerweSt:
  237. yaflUKFSigmaSt base
  238. yaflFloat alpha
  239. yaflFloat beta
  240. yaflFloat kappa
  241. cdef const yaflUKFSigmaMethodsSt yafl_ukf_merwe_spm
  242. #==============================================================================
  243. #Extension API
  244. #==============================================================================
  245. # We need numpy for Pythonic interfaces
  246. cimport numpy as np
  247. import numpy as np#WTF??
  248. import scipy.stats as st
  249. # We need traceback to print pythonic callback exceptions
  250. import traceback as tb
  251. #==============================================================================
  252. IF YAFLPY_USE_64_BIT:
  253. NP_DTYPE = np.float64
  254. ELSE:
  255. NP_DTYPE = np.float32
  256. #Status masks
  257. ST_MSK_REGULARIZED = YAFL_ST_MSK_REGULARIZED
  258. ST_MSK_GLITCH_SMALL = YAFL_ST_MSK_GLITCH_SMALL
  259. ST_MSK_GLITCH_LARGE = YAFL_ST_MSK_GLITCH_LARGE
  260. ST_MSK_ANOMALY = YAFL_ST_MSK_ANOMALY
  261. ST_MSK_ERROR = 0xFF0
  262. ST_MSK_WARNING = 0xF
  263. #Everthing is OK
  264. ST_OK = YAFL_ST_OK
  265. #Error threshold
  266. ST_ERR_THR = YAFL_ST_ERR_THR
  267. #Errors:
  268. # Invalid argument numer
  269. ST_INV_ARG_1 = YAFL_ST_INV_ARG_1
  270. ST_INV_ARG_2 = YAFL_ST_INV_ARG_2
  271. ST_INV_ARG_3 = YAFL_ST_INV_ARG_3
  272. ST_INV_ARG_4 = YAFL_ST_INV_ARG_4
  273. ST_INV_ARG_5 = YAFL_ST_INV_ARG_5
  274. ST_INV_ARG_6 = YAFL_ST_INV_ARG_6
  275. ST_INV_ARG_7 = YAFL_ST_INV_ARG_7
  276. ST_INV_ARG_8 = YAFL_ST_INV_ARG_8
  277. ST_INV_ARG_9 = YAFL_ST_INV_ARG_9
  278. ST_INV_ARG_10 = YAFL_ST_INV_ARG_10
  279. ST_INV_ARG_11 = YAFL_ST_INV_ARG_11
  280. #==============================================================================
  281. # UD-factorized EKF API
  282. #==============================================================================
  283. #------------------------------------------------------------------------------
  284. # Kalman filter basic union
  285. #------------------------------------------------------------------------------
  286. ctypedef union yaflPyKalmanBaseUn:
  287. yaflKalmanBaseSt base
  288. yaflEKFBaseSt ekf
  289. yaflEKFAdaptiveSt ekf_adaptive
  290. yaflEKFRobustSt ekf_robust
  291. yaflEKFAdaptiveRobustSt ekf_ada_rob
  292. yaflUKFBaseSt ukf
  293. yaflUKFAdaptivedSt ukf_adaptive
  294. yaflUKFRobustSt ukf_robust
  295. yaflUKFAdaptiveRobustSt ukf_ada_rob
  296. yaflUKFSt ukf_full
  297. yaflUKFFullAdapiveSt ukf_full_adaptive
  298. #------------------------------------------------------------------------------
  299. # Kalman filter C-structure with Python callback
  300. #------------------------------------------------------------------------------
  301. ctypedef struct yaflPyKalmanBaseSt:
  302. # Kalman filter base union
  303. yaflPyKalmanBaseUn base
  304. # Python/Cython self
  305. void * py_self
  306. #------------------------------------------------------------------------------
  307. cdef int _U_sz(int dim_u):
  308. return max(1, (dim_u * (dim_u - 1))//2)
  309. #------------------------------------------------------------------------------
  310. # Basic Filter class
  311. #------------------------------------------------------------------------------
  312. cdef class yaflKalmanBase:
  313. # Kalman filter C-self
  314. cdef yaflPyKalmanBaseSt c_self
  315. # Kalman filter memory views
  316. cdef yaflFloat [::1] v_x
  317. cdef yaflFloat [::1] v_y
  318. cdef yaflFloat [::1] v_z
  319. cdef yaflFloat [::1] v_Up
  320. cdef yaflFloat [::1] v_Dp
  321. cdef yaflFloat [::1] v_Uq
  322. cdef yaflFloat [::1] v_Dq
  323. cdef yaflFloat [::1] v_Ur
  324. cdef yaflFloat [::1] v_Dr
  325. # Kalman filter numpy arrays
  326. cdef np.ndarray _x
  327. cdef np.ndarray _y
  328. cdef np.ndarray _z
  329. cdef np.ndarray _Up
  330. cdef np.ndarray _Dp
  331. cdef np.ndarray _Uq
  332. cdef np.ndarray _Dq
  333. cdef np.ndarray _Ur
  334. cdef np.ndarray _Dr
  335. # Callback info
  336. cdef yaflFloat _dt
  337. cdef dict _fx_args
  338. cdef object _fx
  339. cdef dict _hx_args
  340. cdef object _hx
  341. cdef object _residual_z
  342. def __init__(self, int dim_x, int dim_z, yaflFloat dt, \
  343. fx, hx, residual_z = None):
  344. #Store dimensions
  345. self.c_self.base.base.Nx = dim_x
  346. self.c_self.base.base.Nz = dim_z
  347. #Setup callbacks
  348. self.c_self.py_self = <void *>self
  349. self._dt = dt
  350. self._fx_args = {}
  351. self._hx_args = {}
  352. if not callable(fx):
  353. raise ValueError('fx must be callable!')
  354. self.c_self.base.base.f = <yaflKalmanFuncP>yafl_py_kalman_fx
  355. self._fx = fx
  356. if not callable(hx):
  357. raise ValueError('hx must be callable!')
  358. self.c_self.base.base.h = <yaflKalmanFuncP>yafl_py_kalman_hx
  359. self._hx = hx
  360. if residual_z:
  361. if not callable(residual_z):
  362. raise ValueError('residual_z must be callable!')
  363. self.c_self.base.base.zrf = <yaflKalmanResFuncP>yafl_py_kalman_zrf
  364. self._residual_z = residual_z
  365. else:
  366. self.c_self.base.base.zrf = <yaflKalmanResFuncP>0
  367. self._residual_z = None
  368. # Allocate memories and setup the rest of c_self
  369. self._z = np.zeros((dim_z,), dtype=NP_DTYPE)
  370. self.v_z = self._z
  371. self._x = np.zeros((dim_x,), dtype=NP_DTYPE)
  372. self.v_x = self._x
  373. self.c_self.base.base.x = &self.v_x[0]
  374. self._y = np.zeros((dim_z,), dtype=NP_DTYPE)
  375. self.v_y = self._y
  376. self.c_self.base.base.y = &self.v_y[0]
  377. self._Up = np.zeros((_U_sz(dim_x),), dtype=NP_DTYPE)
  378. self.v_Up = self._Up
  379. self.c_self.base.base.Up = &self.v_Up[0]
  380. self._Dp = np.ones((dim_x,), dtype=NP_DTYPE)
  381. self.v_Dp = self._Dp
  382. self.c_self.base.base.Dp = &self.v_Dp[0]
  383. self._Uq = np.zeros((_U_sz(dim_x),), dtype=NP_DTYPE)
  384. self.v_Uq = self._Uq
  385. self.c_self.base.base.Uq = &self.v_Uq[0]
  386. self._Dq = np.ones((dim_x,), dtype=NP_DTYPE)
  387. self.v_Dq = self._Dq
  388. self.c_self.base.base.Dq = &self.v_Dq[0]
  389. self._Ur = np.zeros((_U_sz(dim_z),), dtype=NP_DTYPE)
  390. self.v_Ur = self._Ur
  391. self.c_self.base.base.Ur = &self.v_Ur[0]
  392. self._Dr = np.ones((dim_z,), dtype=NP_DTYPE)
  393. self.v_Dr = self._Dr
  394. self.c_self.base.base.Dr = &self.v_Dr[0]
  395. #==========================================================================
  396. #Decorators
  397. @property
  398. def x(self):
  399. return self._x
  400. @x.setter
  401. def x(self, value):
  402. self._x[:] = value
  403. #--------------------------------------------------------------------------
  404. @property
  405. def y(self):
  406. return self._y
  407. @y.setter
  408. def y(self, value):
  409. raise AttributeError('yaflKalmanBase does not support this!')
  410. #--------------------------------------------------------------------------
  411. @property
  412. def Up(self):
  413. return self._Up
  414. @Up.setter
  415. def Up(self, value):
  416. self._Up[:] = value
  417. #--------------------------------------------------------------------------
  418. @property
  419. def Dp(self):
  420. return self._Dp
  421. @Dp.setter
  422. def Dp(self, value):
  423. self._Dp[:] = value
  424. #--------------------------------------------------------------------------
  425. @property
  426. def Uq(self):
  427. return self._Uq
  428. @Uq.setter
  429. def Uq(self, value):
  430. self._Uq[:] = value
  431. #--------------------------------------------------------------------------
  432. @property
  433. def Dq(self):
  434. return self._Dq
  435. @Dq.setter
  436. def Dq(self, value):
  437. self._Dq[:] = value
  438. #--------------------------------------------------------------------------
  439. @property
  440. def Ur(self):
  441. return self._Ur
  442. @Ur.setter
  443. def Ur(self, value):
  444. self._Ur[:] = value
  445. #--------------------------------------------------------------------------
  446. @property
  447. def Dr(self):
  448. return self._Dr
  449. @Dr.setter
  450. def Dr(self, value):
  451. self._Dr[:] = value
  452. #==========================================================================
  453. def _predict(self):
  454. raise NotImplementedError('yaflKalmanBase is the base class!')
  455. def predict(self, dt=None, **fx_args):
  456. old_dt = self._dt
  457. if dt:
  458. if np.isnan(dt):
  459. raise ValueError('Invalid dt value (nan)!')
  460. self._dt = <yaflFloat>dt
  461. self._fx_args = fx_args
  462. res = self._predict()
  463. if res > YAFL_ST_ERR_THR:
  464. raise ValueError('Bad return value on yaflKalmanBase.predict!')
  465. self._dt = old_dt
  466. return res
  467. #==========================================================================
  468. def _update(self):
  469. raise NotImplementedError('yaflKalmanBase is the base class!')
  470. def update(self, z, **hx_args):
  471. self._z[:] = z
  472. self._hx_args = hx_args
  473. res = self._update()
  474. if res > YAFL_ST_ERR_THR:
  475. raise ValueError('Bad return value on yaflKalmanBase.update!')
  476. return res
  477. #------------------------------------------------------------------------------
  478. cdef yaflStatusEn yafl_py_kalman_fx(yaflPyKalmanBaseSt * self, \
  479. yaflFloat * new_x, yaflFloat * old_x):
  480. try:
  481. if not isinstance(<object>(self.py_self), yaflKalmanBase):
  482. raise ValueError('Invalid py_self type (must be subclass of yaflKalmanBase)!')
  483. py_self = <yaflKalmanBase>(self.py_self)
  484. fx = py_self._fx
  485. if not callable(fx):
  486. raise ValueError('fx must be callable!')
  487. dt = py_self._dt
  488. if np.isnan(dt):
  489. raise ValueError('Invalid dt value (nan)!')
  490. fx_args = py_self._fx_args
  491. if not isinstance(fx_args, dict):
  492. raise ValueError('Invalid fx_args type (must be dict)!')
  493. nx = self.base.base.Nx
  494. if nx <= 0:
  495. raise ValueError('nx must be > 0!')
  496. _new_x = np.asarray(<yaflFloat[:nx]> new_x) #np.float64_t
  497. _old_x = np.asarray(<yaflFloat[:nx]> old_x) #np.float64_t
  498. _new_x[:] = fx(_old_x, dt, **fx_args)
  499. return YAFL_ST_OK
  500. except Exception as e:
  501. print(tb.format_exc())
  502. return YAFL_ST_INV_ARG_1
  503. #------------------------------------------------------------------------------
  504. cdef yaflStatusEn yafl_py_kalman_hx(yaflPyKalmanBaseSt * self, \
  505. yaflFloat * z, yaflFloat * x):
  506. try:
  507. if not isinstance(<object>(self.py_self), yaflKalmanBase):
  508. raise ValueError('Invalid py_self type (must be subclass of yaflKalmanBase)!')
  509. py_self = <yaflKalmanBase>(self.py_self)
  510. hx = py_self._hx
  511. if not callable(hx):
  512. raise ValueError('hx must be callable!')
  513. hx_args = py_self._hx_args
  514. if not isinstance(hx_args, dict):
  515. raise ValueError('Invalid hx_args type (must be dict)!')
  516. nx = self.base.base.Nx
  517. if nx <= 0:
  518. raise ValueError('nx must be > 0!')
  519. nz = self.base.base.Nz
  520. if nx <= 0:
  521. raise ValueError('nz must be > 0!')
  522. _x = np.asarray(<yaflFloat[:nx]> x) #np.float64_t
  523. _z = np.asarray(<yaflFloat[:nz]> z) #np.float64_t
  524. _z[:] = hx(_x, **hx_args)
  525. return YAFL_ST_OK
  526. except Exception as e:
  527. print(tb.format_exc())
  528. return YAFL_ST_INV_ARG_1
  529. #------------------------------------------------------------------------------
  530. cdef yaflStatusEn yafl_py_kalman_zrf(yaflPyKalmanBaseSt * self, yaflFloat * res, \
  531. yaflFloat * sigma, yaflFloat * pivot):
  532. try:
  533. if not isinstance(<object>(self.py_self), yaflKalmanBase):
  534. raise ValueError('Invalid py_self type (must be subclass of yaflKalmanBase)!')
  535. py_self = <yaflKalmanBase>(self.py_self)
  536. residual_z = py_self._residual_z
  537. if not callable(residual_z):
  538. raise ValueError('residual_z must be callable!')
  539. nz = self.base.base.Nz
  540. if nz <= 0:
  541. raise ValueError('nx must be > 0!')
  542. _res = np.asarray(<yaflFloat[:nz]> res) #np.float64_t
  543. _sigma = np.asarray(<yaflFloat[:nz]> sigma) #np.float64_t
  544. _pivot = np.asarray(<yaflFloat[:nz]> pivot) #np.float64_t
  545. _res[:] = residual_z(_sigma, _pivot)
  546. return YAFL_ST_OK
  547. except Exception as e:
  548. print(tb.format_exc())
  549. return YAFL_ST_INV_ARG_1
  550. #==============================================================================
  551. # Basic Filter class
  552. #==============================================================================
  553. cdef class yaflExtendedBase(yaflKalmanBase):
  554. # Kalman filter memory views
  555. cdef yaflFloat [:, ::1] v_H
  556. cdef yaflFloat [:, ::1] v_W
  557. cdef yaflFloat [::1] v_D
  558. # Kalman filter numpy arrays
  559. cdef np.ndarray _H
  560. cdef np.ndarray _W
  561. cdef np.ndarray _D
  562. # Callback info
  563. cdef object _jfx
  564. cdef object _jhx
  565. #The object will be Extensible
  566. cdef dict __dict__
  567. def __init__(self, int dim_x, int dim_z, yaflFloat dt, \
  568. fx, jfx, hx, jhx, residual_z = None):
  569. super().__init__(dim_x, dim_z, dt, fx, hx, residual_z)
  570. if not callable(jfx):
  571. raise ValueError('jfx must be callable!')
  572. self.c_self.base.ekf.jf = <yaflKalmanFuncP>yafl_py_ekf_jfx
  573. self._jfx = jfx
  574. if not callable(jhx):
  575. raise ValueError('jhx must be callable!')
  576. self.c_self.base.ekf.jh = <yaflKalmanFuncP>yafl_py_ekf_jhx
  577. self._jhx = jhx
  578. # Allocate memories and setup the rest of c_self
  579. self._H = np.zeros((dim_z, dim_x), dtype=NP_DTYPE)
  580. self.v_H = self._H
  581. self.c_self.base.ekf.H = &self.v_H[0, 0]
  582. self._W = np.zeros((dim_x, 2 * dim_x), dtype=NP_DTYPE)
  583. self.v_W = self._W
  584. self.c_self.base.ekf.W = &self.v_W[0,0]
  585. self._D = np.ones((2 * dim_x,), dtype=NP_DTYPE)
  586. self.v_D = self._D
  587. self.c_self.base.ekf.D = &self.v_D[0]
  588. #==========================================================================
  589. def _predict(self):
  590. return yafl_ekf_base_predict(&(self.c_self.base.base))
  591. #==========================================================================
  592. def _update(self):
  593. raise NotImplementedError('yaflExtendedBase is the base class!')
  594. #------------------------------------------------------------------------------
  595. # State transition function Jacobian
  596. cdef yaflStatusEn yafl_py_ekf_jfx(yaflPyKalmanBaseSt * self, yaflFloat * w, \
  597. yaflFloat * x):
  598. try:
  599. if not isinstance(<object>(self.py_self), yaflExtendedBase):
  600. raise ValueError('Invalid py_self type (must be subclass of yaflExtendedBase)!')
  601. py_self = <yaflExtendedBase>(self.py_self)
  602. jfx = py_self._jfx
  603. if not callable(jfx):
  604. raise ValueError('jfx must be callable!')
  605. dt = py_self._dt
  606. if np.isnan(dt):
  607. raise ValueError('Invalid dt value (nan)!')
  608. fx_args = py_self._fx_args
  609. if not isinstance(fx_args, dict):
  610. raise ValueError('Invalid fx_args type (must be dict)!')
  611. #How about handling exceptions here???
  612. py_self._W[:, :self.base.base.Nx] = jfx(py_self._x, dt, **fx_args)
  613. return YAFL_ST_OK
  614. except Exception as e:
  615. print(tb.format_exc())
  616. return YAFL_ST_INV_ARG_1
  617. #------------------------------------------------------------------------------
  618. # State transition function Jacobian
  619. cdef yaflStatusEn yafl_py_ekf_jhx(yaflPyKalmanBaseSt * self, yaflFloat * h, \
  620. yaflFloat * x):
  621. try:
  622. if not isinstance(<object>(self.py_self), yaflExtendedBase):
  623. raise ValueError('Invalid py_self type (must be subclass of yaflExtendedBase)!')
  624. py_self = <yaflExtendedBase>(self.py_self)
  625. jhx = py_self._jhx
  626. if not callable(jhx):
  627. raise ValueError('jhx must be callable!')
  628. hx_args = py_self._hx_args
  629. if not isinstance(hx_args, dict):
  630. raise ValueError('Invalid hx_args type (must be dict)!')
  631. #How about handling exceptions here???
  632. py_self._H[:,:] = jhx(py_self._x, **hx_args)
  633. return YAFL_ST_OK
  634. except Exception as e:
  635. print(tb.format_exc())
  636. return YAFL_ST_INV_ARG_1
  637. #==============================================================================
  638. cdef class Bierman(yaflExtendedBase):
  639. def _update(self):
  640. return yafl_ekf_base_update(&self.c_self.base.base, &self.v_z[0], \
  641. yafl_ekf_bierman_update_scalar)
  642. #------------------------------------------------------------------------------
  643. cdef class Joseph(yaflExtendedBase):
  644. def _update(self):
  645. return yafl_ekf_base_update(&self.c_self.base.base, &self.v_z[0], \
  646. yafl_ekf_joseph_update_scalar)
  647. #==============================================================================
  648. # Adaptive filter basic class
  649. #==============================================================================
  650. cdef class yaflAdaptiveBase(yaflExtendedBase):
  651. def __init__(self, int dim_x, int dim_z, yaflFloat dt, \
  652. fx, jfx, hx, jhx, **kwargs):
  653. super().__init__(dim_x, dim_z, dt, fx, jfx, hx, jhx, **kwargs)
  654. #Init chi2 with scipy.stats.chi2.ppf(0.999, 1)
  655. self.c_self.base.ekf_adaptive.chi2 = 10.8275662
  656. #==========================================================================
  657. #Decorators
  658. @property
  659. def chi2(self):
  660. return self.c_self.base.ekf_adaptive.chi2
  661. @chi2.setter
  662. def chi2(self, value):
  663. self.c_self.base.ekf_adaptive.chi2 = <yaflFloat>value
  664. #==============================================================================
  665. cdef class AdaptiveBierman(yaflAdaptiveBase):
  666. def _update(self):
  667. return yafl_ekf_base_update(&self.c_self.base.base, &self.v_z[0], \
  668. yafl_ekf_adaptive_bierman_update_scalar)
  669. #------------------------------------------------------------------------------
  670. cdef class AdaptiveJoseph(yaflAdaptiveBase):
  671. def _update(self):
  672. return yafl_ekf_base_update(&self.c_self.base.base, &self.v_z[0], \
  673. yafl_ekf_adaptive_joseph_update_scalar)
  674. #------------------------------------------------------------------------------
  675. cdef class DoNotUseThisFilter(yaflAdaptiveBase):
  676. """
  677. WARNING!!!
  678. DO NOT USE THIS variant of Adaptive Joseph filter !!!
  679. It was implemented to show some flaws of the corresponding algorithm!
  680. """
  681. def _update(self):
  682. return yafl_ekf_base_update(&self.c_self.base.base, &self.v_z[0], \
  683. yafl_ekf_do_not_use_this_update_scalar)
  684. #==============================================================================
  685. # Robust filter basic class
  686. #==============================================================================
  687. cdef class yaflRobustBase(yaflExtendedBase):
  688. cdef object _gz
  689. cdef object _gdotz
  690. def __init__(self, int dim_x, int dim_z, yaflFloat dt, \
  691. fx, jfx, hx, jhx, gz=None, gdotz=None, **kwargs):
  692. super().__init__(dim_x, dim_z, dt, fx, jfx, hx, jhx, **kwargs)
  693. if gz:
  694. if not callable(gz):
  695. raise ValueError('gz must be callable!')
  696. if not gdotz:
  697. raise ValueError('gdotz must be passed!')
  698. if not callable(gdotz):
  699. raise ValueError('gdotz must be callable!')
  700. self._gz = gz
  701. self._gdotz = gdotz
  702. self.c_self.base.ekf_robust.g = <yaflKalmanRobFuncP>yafl_py_ekf_rob_gz
  703. self.c_self.base.ekf_robust.gdot = <yaflKalmanRobFuncP>yafl_py_ekf_rob_gdotz
  704. else:
  705. self.c_self.base.ekf_robust.g = <yaflKalmanRobFuncP>0
  706. self.c_self.base.ekf_robust.gdot = <yaflKalmanRobFuncP>0
  707. #------------------------------------------------------------------------------
  708. # Influence limiting function
  709. cdef yaflFloat yafl_py_ekf_rob_gz(yaflPyKalmanBaseSt * self, yaflFloat nu):
  710. try:
  711. if not isinstance(<object>(self.py_self), yaflRobustBase):
  712. raise ValueError('Invalid py_self type (must be subclass of yaflRobustBase)!')
  713. py_self = <yaflRobustBase>(self.py_self)
  714. gz = py_self._gz
  715. if not callable(gz):
  716. raise ValueError('gz must be callable!')
  717. hx_args = py_self._hx_args
  718. if not isinstance(hx_args, dict):
  719. raise ValueError('Invalid hx_args type (must be dict)!')
  720. #How about handling exceptions here???
  721. ret = gz(nu, **hx_args)
  722. if type(ret) != float:
  723. raise ValueError('gz must return float!')
  724. return <yaflFloat>ret
  725. except Exception as e:
  726. print(tb.format_exc())
  727. return <yaflFloat>0.0
  728. #------------------------------------------------------------------------------
  729. # Influence limiting function derivative
  730. cdef yaflFloat yafl_py_ekf_rob_gdotz(yaflPyKalmanBaseSt * self, yaflFloat nu):
  731. try:
  732. if not isinstance(<object>(self.py_self), yaflRobustBase):
  733. raise ValueError('Invalid py_self type (must be subclass of yaflRobustBase)!')
  734. py_self = <yaflRobustBase>(self.py_self)
  735. gdotz = py_self._gdotz
  736. if not callable(gdotz):
  737. raise ValueError('gdotz must be callable!')
  738. hx_args = py_self._hx_args
  739. if not isinstance(hx_args, dict):
  740. raise ValueError('Invalid hx_args type (must be dict)!')
  741. #How about handling exceptions here???
  742. ret = gdotz(nu, **hx_args)
  743. if type(ret) != float:
  744. raise ValueError('gdotz must return float!')
  745. return <yaflFloat>ret
  746. except Exception as e:
  747. print(tb.format_exc())
  748. return <yaflFloat>0.0
  749. #==============================================================================
  750. cdef class RobustBierman(yaflRobustBase):
  751. def _update(self):
  752. return yafl_ekf_base_update(&self.c_self.base.base, &self.v_z[0], \
  753. yafl_ekf_robust_bierman_update_scalar)
  754. #------------------------------------------------------------------------------
  755. cdef class RobustJoseph(yaflRobustBase):
  756. def _update(self):
  757. return yafl_ekf_base_update(&self.c_self.base.base, &self.v_z[0], \
  758. yafl_ekf_robust_joseph_update_scalar)
  759. #==============================================================================
  760. # Adaptive robust filter basic class
  761. #==============================================================================
  762. cdef class yaflAdaptiveRobustBase(yaflRobustBase):
  763. def __init__(self, int dim_x, int dim_z, yaflFloat dt, \
  764. fx, jfx, hx, jhx, **kwargs):
  765. super().__init__(dim_x, dim_z, dt, fx, jfx, hx, jhx, **kwargs)
  766. #Init chi2 with scipy.stats.chi2.ppf(0.997, 1)
  767. self.c_self.base.ekf_ada_rob.chi2 = 8.807468393511947
  768. #==========================================================================
  769. #Decorators
  770. @property
  771. def chi2(self):
  772. return self.c_self.base.ekf_ada_rob.chi2
  773. @chi2.setter
  774. def chi2(self, value):
  775. self.c_self.base.ekf_ada_rob.chi2 = <yaflFloat>value
  776. #==============================================================================
  777. cdef class AdaptiveRobustBierman(yaflAdaptiveRobustBase):
  778. def _update(self):
  779. return yafl_ekf_base_update(&self.c_self.base.base, &self.v_z[0], \
  780. yafl_ekf_adaptive_robust_bierman_update_scalar)
  781. #------------------------------------------------------------------------------
  782. cdef class AdaptiveRobustJoseph(yaflAdaptiveRobustBase):
  783. def _update(self):
  784. return yafl_ekf_base_update(&self.c_self.base.base, &self.v_z[0], \
  785. yafl_ekf_adaptive_robust_joseph_update_scalar)
  786. #==============================================================================
  787. # UD-factorized UKF API
  788. #==============================================================================
  789. # Sigma points generator basic definitions
  790. #------------------------------------------------------------------------------
  791. ctypedef union yaflPySigmaBaseUn:
  792. yaflUKFSigmaSt base
  793. yaflUKFMerweSt merwe
  794. #------------------------------------------------------------------------------
  795. ctypedef struct yaflPySigmaSt:
  796. # Sigma point generator base structure
  797. yaflPySigmaBaseUn base
  798. # Python/Cython self
  799. void * py_self
  800. #------------------------------------------------------------------------------
  801. cdef class yaflSigmaBase:
  802. # Sigma point generator C-self
  803. cdef yaflPySigmaSt c_self
  804. # Callback info
  805. cdef object _addf
  806. #The object will be Extensible
  807. cdef dict __dict__
  808. def __init__(self, yaflInt dim_x, addf=None):
  809. #Setup callbacks
  810. self.c_self.py_self = <void *>self
  811. if addf:
  812. if not callable(addf):
  813. raise ValueError('addf must be callable!')
  814. self.c_self.base.base.addf = <yaflUKFSigmaAddP>yafl_py_sigma_addf
  815. self._addf = addf
  816. else:
  817. self.c_self.base.base.addf = <yaflUKFSigmaAddP>0
  818. pnum = self.get_np(dim_x)
  819. self.c_self.base.base.np = pnum
  820. cdef yaflInt get_np(self, int dim_x):
  821. raise NotImplementedError('yaflSigmaBase is the base class!')
  822. cdef const yaflUKFSigmaMethodsSt * get_spm(self):
  823. raise NotImplementedError('yaflSigmaBase is the base class!')
  824. @property
  825. def pnum(self):
  826. return self.c_self.base.base.np
  827. @pnum.setter
  828. def pnum(self, value):
  829. raise AttributeError('yaflSigmaBase does not support this!')
  830. #------------------------------------------------------------------------------
  831. # UD-factorized UKF definitions
  832. #------------------------------------------------------------------------------
  833. cdef class yaflUnscentedBase(yaflKalmanBase):
  834. # Kalman filter memory views
  835. cdef yaflFloat [::1] v_zp
  836. cdef yaflFloat [::1] v_Sx
  837. cdef yaflFloat [:, ::1] v_Pzx
  838. cdef yaflFloat [:, ::1] v_sigmas_x
  839. cdef yaflFloat [:, ::1] v_sigmas_z
  840. cdef yaflFloat [::1] v_wm
  841. cdef yaflFloat [::1] v_wc
  842. # Kalman filter numpy arrays
  843. cdef np.ndarray _zp
  844. cdef np.ndarray _Sx
  845. cdef np.ndarray _Pzx
  846. cdef np.ndarray _sigmas_x
  847. cdef np.ndarray _sigmas_z
  848. cdef np.ndarray _wm
  849. cdef np.ndarray _wc
  850. # Callback info
  851. cdef object _points
  852. cdef object _mean_x
  853. cdef object _resudual_x
  854. cdef object _mean_z
  855. #The object will be Extensible
  856. cdef dict __dict__
  857. def __init__(self, int dim_x, int dim_z, yaflFloat dt, hx, fx, points,\
  858. x_mean_fn=None, z_mean_fn=None, \
  859. residual_x=None, residual_z=None):
  860. super().__init__(dim_x, dim_z, dt, fx, hx, residual_z)
  861. if x_mean_fn:
  862. if not callable(x_mean_fn):
  863. raise ValueError('x_mean_fn must be callable!')
  864. self.c_self.base.ukf.xmf = <yaflKalmanFuncP>yafl_py_ukf_xmf
  865. self._mean_x = x_mean_fn
  866. else:
  867. self.c_self.base.ukf.xmf = <yaflKalmanFuncP>0
  868. self._mean_x = None
  869. if residual_x:
  870. if not callable(residual_x):
  871. raise ValueError('residual_x must be callable!')
  872. self.c_self.base.ukf.xrf = <yaflKalmanResFuncP>yafl_py_ukf_xrf
  873. self._residual_x = residual_x
  874. else:
  875. self.c_self.base.ukf.xrf = <yaflKalmanResFuncP>0
  876. self._residual_x = None
  877. if z_mean_fn:
  878. if not callable(z_mean_fn):
  879. raise ValueError('z_mean_fn must be callable!')
  880. self.c_self.base.ukf.zmf = <yaflKalmanFuncP>yafl_py_ukf_zmf
  881. self._mean_z = z_mean_fn
  882. else:
  883. self.c_self.base.ukf.zmf = <yaflKalmanFuncP>0
  884. self._mean_z = None
  885. #Setup sigma points generator
  886. if not isinstance(points, yaflSigmaBase):
  887. raise ValueError('Invalid points type (must be subclass of yaflSigmaBase)!')
  888. _points = <yaflSigmaBase>points
  889. self._points = _points
  890. self.c_self.base.ukf.sp_info = &_points.c_self.base.base
  891. self.c_self.base.ukf.sp_meth = _points.get_spm()
  892. # Allocate memories and setup the rest of c_self
  893. # Sigma points and weights
  894. pnum = _points.pnum
  895. self._wm = np.zeros((pnum,), dtype=NP_DTYPE)
  896. self.v_wm = self._wm
  897. self.c_self.base.ukf.wm = &self.v_wm[0]
  898. self._wc = np.zeros((pnum,), dtype=NP_DTYPE)
  899. self.v_wc = self._wc
  900. self.c_self.base.ukf.wc = &self.v_wc[0]
  901. self._sigmas_x = np.zeros((pnum, dim_x), dtype=NP_DTYPE)
  902. self.v_sigmas_x = self._sigmas_x
  903. self.c_self.base.ukf.sigmas_x = &self.v_sigmas_x[0, 0]
  904. self._sigmas_z = np.zeros((pnum, dim_z), dtype=NP_DTYPE)
  905. self.v_sigmas_z = self._sigmas_z
  906. self.c_self.base.ukf.sigmas_z = &self.v_sigmas_z[0, 0]
  907. #Rest of rhe UKF
  908. self._zp = np.zeros((dim_z,), dtype=NP_DTYPE)
  909. self.v_zp = self._zp
  910. self.c_self.base.ukf.zp = &self.v_zp[0]
  911. self._Pzx = np.zeros((dim_z, dim_x), dtype=NP_DTYPE)
  912. self.v_Pzx = self._Pzx
  913. self.c_self.base.ukf.Pzx = &self.v_Pzx[0, 0]
  914. self._Sx = np.zeros((dim_x,), dtype=NP_DTYPE)
  915. self.v_Sx = self._Sx
  916. self.c_self.base.ukf.Sx = &self.v_Sx[0]
  917. #Call C-post init
  918. yafl_ukf_post_init(&self.c_self.base.ukf)
  919. #==========================================================================
  920. #Decorators
  921. @property
  922. def Pzx(self):
  923. return self._Pzx
  924. #--------------------------------------------------------------------------
  925. @property
  926. def zp(self):
  927. return self._zp
  928. #--------------------------------------------------------------------------
  929. @property
  930. def sigmas_x(self):
  931. return self._sigmas_x
  932. @sigmas_x.setter
  933. def sigmas_x(self, value):
  934. raise AttributeError('yaflUnscentedBase does not support this!')
  935. #--------------------------------------------------------------------------
  936. @property
  937. def sigmas_z(self):
  938. return self._sigmas_z
  939. @sigmas_z.setter
  940. def sigmas_z(self, value):
  941. raise AttributeError('yaflUnscentedBase does not support this!')
  942. #--------------------------------------------------------------------------
  943. @property
  944. def wc(self):
  945. return self._wc
  946. @wc.setter
  947. def wc(self, value):
  948. raise AttributeError('yaflUnscentedBase does not support this!')
  949. #--------------------------------------------------------------------------
  950. @property
  951. def wm(self):
  952. return self._wm
  953. @wm.setter
  954. def wm(self, value):
  955. raise AttributeError('yaflUnscentedBase does not support this!')
  956. #==========================================================================
  957. def _predict(self):
  958. return yafl_ukf_base_predict(&self.c_self.base.ukf)
  959. #==========================================================================
  960. def _update(self):
  961. raise NotImplementedError('yaflUnscentedBase is the base class!')
  962. #==============================================================================
  963. cdef yaflStatusEn yafl_py_sigma_addf(yaflPyKalmanBaseSt * self, yaflFloat * delta, \
  964. yaflFloat * pivot, yaflFloat mult):
  965. try:
  966. if not isinstance(<object>(self.py_self), yaflUnscentedBase):
  967. raise ValueError('Invalid py_self type (must be subclass of yaflUnscentedBase)!')
  968. py_self = <yaflUnscentedBase>(self.py_self)
  969. _addf = py_self._points._addf
  970. if not callable(_addf):
  971. raise ValueError('_addf must be callable!')
  972. nx = self.base.base.Nx
  973. if nx <= 0:
  974. raise ValueError('nx must be > 0!')
  975. _delta = np.asarray(<yaflFloat[:nx]> delta) #np.float64_t
  976. _pivot = np.asarray(<yaflFloat[:nx]> pivot) #np.float64_t
  977. _delta[:] = _addf(_delta, _pivot, mult)
  978. return YAFL_ST_OK
  979. except Exception as e:
  980. print(tb.format_exc())
  981. return YAFL_ST_INV_ARG_1
  982. #------------------------------------------------------------------------------
  983. cdef yaflStatusEn yafl_py_ukf_xmf(yaflPyKalmanBaseSt * self, \
  984. yaflFloat * res, yaflFloat * sigmas):
  985. try:
  986. if not isinstance(<object>(self.py_self), yaflUnscentedBase):
  987. raise ValueError('Invalid py_self type (must be subclass of yaflUnscentedBase)!')
  988. py_self = <yaflUnscentedBase>(self.py_self)
  989. mean_x = py_self._mean_x
  990. if not callable(mean_x):
  991. raise ValueError('mean_x must be callable!')
  992. nx = self.base.base.Nx
  993. if nx <= 0:
  994. raise ValueError('nx must be > 0!')
  995. _points = py_self._points
  996. if not isinstance(py_self, yaflSigmaBase):
  997. raise ValueError('Invalid _points type (must be subclass of yaflSigmaBase)!')
  998. pnum = _points.base.base.np
  999. if pnum <= 0:
  1000. raise ValueError('pnum must be > 0!')
  1001. _sigmas = np.asarray(<yaflFloat[:pnum, :nx]> sigmas) #np.float64_t
  1002. _res = np.asarray(<yaflFloat[:nx]> res) #np.float64_t
  1003. _res[:] = mean_x(_sigmas, _points._wm)
  1004. return YAFL_ST_OK
  1005. except Exception as e:
  1006. print(tb.format_exc())
  1007. return YAFL_ST_INV_ARG_1
  1008. #------------------------------------------------------------------------------
  1009. cdef yaflStatusEn yafl_py_ukf_xrf(yaflPyKalmanBaseSt * self, yaflFloat * res, \
  1010. yaflFloat * sigma, yaflFloat * pivot):
  1011. try:
  1012. if not isinstance(<object>(self.py_self), yaflUnscentedBase):
  1013. raise ValueError('Invalid py_self type (must be subclass of yaflUnscentedBase)!')
  1014. py_self = <yaflUnscentedBase>(self.py_self)
  1015. residual_x = py_self._residual_x
  1016. if not callable(residual_x):
  1017. raise ValueError('residual_x must be callable!')
  1018. nx = self.base.base.Nx
  1019. if nx <= 0:
  1020. raise ValueError('nx must be > 0!')
  1021. _res = np.asarray(<yaflFloat[:nx]> res) #np.float64_t
  1022. _sigma = np.asarray(<yaflFloat[:nx]> sigma) #np.float64_t
  1023. _pivot = np.asarray(<yaflFloat[:nx]> pivot) #np.float64_t
  1024. _res[:] = residual_x(_sigma, _pivot)
  1025. return YAFL_ST_OK
  1026. except Exception as e:
  1027. print(tb.format_exc())
  1028. return YAFL_ST_INV_ARG_1
  1029. #------------------------------------------------------------------------------
  1030. cdef yaflStatusEn yafl_py_ukf_zmf(yaflPyKalmanBaseSt * self, \
  1031. yaflFloat * res, yaflFloat * sigmas):
  1032. try:
  1033. if not isinstance(<object>(self.py_self), yaflUnscentedBase):
  1034. raise ValueError('Invalid py_self type (must be subclass of yaflUnscentedBase)!')
  1035. py_self = <yaflUnscentedBase>(self.py_self)
  1036. mean_z = py_self._mean_z
  1037. if not callable(mean_z):
  1038. raise ValueError('mean_z must be callable!')
  1039. nz = self.base.base.Nz
  1040. if nz <= 0:
  1041. raise ValueError('nz must be > 0!')
  1042. _points = py_self._points
  1043. if not isinstance(py_self, yaflSigmaBase):
  1044. raise ValueError('Invalid _points type (must be subclass of yaflSigmaBase)!')
  1045. pnum = _points.base.base.np
  1046. if pnum <= 0:
  1047. raise ValueError('pnum must be > 0!')
  1048. _sigmas = np.asarray(<yaflFloat[:pnum, :nz]> sigmas) #np.float64_t
  1049. _res = np.asarray(<yaflFloat[:nz]> res) #np.float64_t
  1050. _res[:] = mean_z(_sigmas, _points._wm)
  1051. return YAFL_ST_OK
  1052. except Exception as e:
  1053. print(tb.format_exc())
  1054. return YAFL_ST_INV_ARG_1
  1055. #==============================================================================
  1056. cdef class UnscentedBierman(yaflUnscentedBase):
  1057. """
  1058. UD-factorized UKF implementation
  1059. """
  1060. def _update(self):
  1061. return yafl_ukf_base_update(&self.c_self.base.ukf, &self.v_z[0], \
  1062. yafl_ukf_bierman_update_scalar)
  1063. #==============================================================================
  1064. cdef class UnscentedAdaptiveBierman(yaflUnscentedBase):
  1065. def __init__(self, int dim_x, int dim_z, yaflFloat dt, hx, fx, points, \
  1066. **kwargs):
  1067. super().__init__(dim_x, dim_z, dt, hx, fx, points, **kwargs)
  1068. #Init chi2 with scipy.stats.chi2.ppf(0.99999, 1)
  1069. self.c_self.base.ukf_adaptive.chi2 = 19.511420964666268
  1070. #==========================================================================
  1071. #Decorators
  1072. @property
  1073. def chi2(self):
  1074. return self.c_self.base.ukf_adaptive.chi2
  1075. @chi2.setter
  1076. def chi2(self, value):
  1077. self.c_self.base.ukf_adaptive.chi2 = <yaflFloat>value
  1078. """
  1079. UD-factorized UKF implementation
  1080. """
  1081. def _update(self):
  1082. return yafl_ukf_base_update(&self.c_self.base.ukf, &self.v_z[0], \
  1083. yafl_ukf_adaptive_bierman_update_scalar)
  1084. #==============================================================================
  1085. # Robust filter basic class
  1086. #==============================================================================
  1087. cdef class yaflRobustUKFBase(yaflUnscentedBase):
  1088. cdef object _gz
  1089. cdef object _gdotz
  1090. def __init__(self, int dim_x, int dim_z, yaflFloat dt, \
  1091. hx, fx, points, gz=None, gdotz=None, **kwargs):
  1092. super().__init__(dim_x, dim_z, dt, hx, fx, points, **kwargs)
  1093. if gz:
  1094. if not callable(gz):
  1095. raise ValueError('gz must be callable!')
  1096. if not gdotz:
  1097. raise ValueError('gdotz must be passed!')
  1098. if not callable(gdotz):
  1099. raise ValueError('gdotz must be callable!')
  1100. self._gz = gz
  1101. self._gdotz = gdotz
  1102. self.c_self.base.ukf_robust.g = <yaflKalmanRobFuncP>yafl_py_ukf_rob_gz
  1103. self.c_self.base.ukf_robust.gdot = <yaflKalmanRobFuncP>yafl_py_ukf_rob_gdotz
  1104. else:
  1105. self.c_self.base.ukf_robust.g = <yaflKalmanRobFuncP>0
  1106. self.c_self.base.ukf_robust.gdot = <yaflKalmanRobFuncP>0
  1107. #------------------------------------------------------------------------------
  1108. # Influence limiting function
  1109. cdef yaflFloat yafl_py_ukf_rob_gz(yaflPyKalmanBaseSt * self, yaflFloat nu):
  1110. try:
  1111. if not isinstance(<object>(self.py_self), yaflRobustUKFBase):
  1112. raise ValueError('Invalid py_self type (must be subclass of yaflRobustBase)!')
  1113. py_self = <yaflRobustUKFBase>(self.py_self)
  1114. gz = py_self._gz
  1115. if not callable(gz):
  1116. raise ValueError('gz must be callable!')
  1117. hx_args = py_self._hx_args
  1118. if not isinstance(hx_args, dict):
  1119. raise ValueError('Invalid hx_args type (must be dict)!')
  1120. #How about handling exceptions here???
  1121. ret = gz(nu, **hx_args)
  1122. if type(ret) != float:
  1123. raise ValueError('gz must return float!')
  1124. return <yaflFloat>ret
  1125. except Exception as e:
  1126. print(tb.format_exc())
  1127. return <yaflFloat>0.0
  1128. #------------------------------------------------------------------------------
  1129. # Influence limiting function derivative
  1130. cdef yaflFloat yafl_py_ukf_rob_gdotz(yaflPyKalmanBaseSt * self, yaflFloat nu):
  1131. try:
  1132. if not isinstance(<object>(self.py_self), yaflRobustUKFBase):
  1133. raise ValueError('Invalid py_self type (must be subclass of yaflRobustBase)!')
  1134. py_self = <yaflRobustUKFBase>(self.py_self)
  1135. gdotz = py_self._gdotz
  1136. if not callable(gdotz):
  1137. raise ValueError('gdotz must be callable!')
  1138. hx_args = py_self._hx_args
  1139. if not isinstance(hx_args, dict):
  1140. raise ValueError('Invalid hx_args type (must be dict)!')
  1141. #How about handling exceptions here???
  1142. ret = gdotz(nu, **hx_args)
  1143. if type(ret) != float:
  1144. raise ValueError('gdotz must return float!')
  1145. return <yaflFloat>ret
  1146. except Exception as e:
  1147. print(tb.format_exc())
  1148. return <yaflFloat>0.0
  1149. #==============================================================================
  1150. cdef class UnscentedRobustBierman(yaflRobustUKFBase):
  1151. def _update(self):
  1152. return yafl_ukf_base_update(&self.c_self.base.ukf, &self.v_z[0], \
  1153. yafl_ukf_robust_bierman_update_scalar)
  1154. #==============================================================================
  1155. # Adaptive robust UKF base
  1156. #==============================================================================
  1157. cdef class yaflAdaptiveRobustUKFBase(yaflRobustUKFBase):
  1158. def __init__(self, int dim_x, int dim_z, yaflFloat dt, \
  1159. hx, fx, points, **kwargs):
  1160. super().__init__(dim_x, dim_z, dt, hx, fx, points, **kwargs)
  1161. #Init chi2 with scipy.stats.chi2.ppf(0.997, 1)
  1162. self.c_self.base.ukf_ada_rob.chi2 = 8.807468393511947
  1163. #==========================================================================
  1164. #Decorators
  1165. @property
  1166. def chi2(self):
  1167. return self.c_self.base.ukf_ada_rob.chi2
  1168. @chi2.setter
  1169. def chi2(self, value):
  1170. self.c_self.base.ukf_ada_rob.chi2 = <yaflFloat>value
  1171. #==============================================================================
  1172. cdef class UnscentedAdaptiveRobustBierman(yaflAdaptiveRobustUKFBase):
  1173. def _update(self):
  1174. return yafl_ukf_base_update(&self.c_self.base.ukf, &self.v_z[0], \
  1175. yafl_ukf_adaptive_robust_bierman_update_scalar)
  1176. #==============================================================================
  1177. # Full UKF, not sequential square root version of UKF
  1178. #==============================================================================
  1179. cdef class Unscented(yaflUnscentedBase):
  1180. cdef yaflFloat [::1] v_Us
  1181. cdef yaflFloat [::1] v_Ds
  1182. def __init__(self, int dim_x, int dim_z, yaflFloat dt, hx, fx, points, \
  1183. **kwargs):
  1184. super().__init__(dim_x, dim_z, dt, hx, fx, points, **kwargs)
  1185. self._Us = np.zeros((_U_sz(dim_z),), dtype=NP_DTYPE)
  1186. self.v_Us = self._Us
  1187. self.c_self.base.ukf_full.Us = &self.v_Us[0]
  1188. self._Ds = np.ones((dim_z,), dtype=NP_DTYPE)
  1189. self.v_Ds = self._Ds
  1190. self.c_self.base.ukf_full.Ds = &self.v_Ds[0]
  1191. #==========================================================================
  1192. #Decorators
  1193. @property
  1194. def Us(self):
  1195. return self._Us
  1196. @Us.setter
  1197. def Us(self, value):
  1198. raise AttributeError('Unscented does not support this!')
  1199. #--------------------------------------------------------------------------
  1200. @property
  1201. def Ds(self):
  1202. return self._Ds
  1203. @Ds.setter
  1204. def Ds(self, value):
  1205. raise AttributeError('Unscented does not support this!')
  1206. #==========================================================================
  1207. """
  1208. UD-factorized UKF implementation
  1209. """
  1210. def _update(self):
  1211. return yafl_ukf_update(&self.c_self.base.ukf, &self.v_z[0])
  1212. #==============================================================================
  1213. # Full adaptive UKF, not sequential square root version of UKF
  1214. #==============================================================================
  1215. cdef class UnscentedAdaptive(Unscented):
  1216. def __init__(self, int dim_x, int dim_z, yaflFloat dt, hx, fx, points, \
  1217. aplha = 0.000001, **kwargs):
  1218. super().__init__(dim_x, dim_z, dt, hx, fx, points, **kwargs)
  1219. self.c_self.base.ukf_full_adaptive.chi2 = \
  1220. st.chi2.ppf(1.0 - aplha, dim_z)
  1221. #==========================================================================
  1222. #Decorators
  1223. @property
  1224. def chi2(self):
  1225. return self.c_self.base.ukf_full_adaptive.chi2
  1226. @chi2.setter
  1227. def chi2(self, value):
  1228. self.c_self.base.ukf_full_adaptive.chi2 = <yaflFloat>value
  1229. #==========================================================================
  1230. """
  1231. UD-factorized UKF implementation
  1232. """
  1233. def _update(self):
  1234. return yafl_ukf_adaptive_update(&self.c_self.base.ukf, &self.v_z[0])
  1235. #==============================================================================
  1236. cdef class MerweSigmaPoints(yaflSigmaBase):
  1237. """
  1238. Van der Merwe sigma point generator implementation
  1239. """
  1240. def __init__(self, yaflInt dim_x, yaflFloat alpha, yaflFloat beta, \
  1241. yaflFloat kappa=0.0, **kwargs):
  1242. super().__init__(dim_x, **kwargs)
  1243. self.c_self.base.merwe.alpha = alpha
  1244. self.c_self.base.merwe.beta = beta
  1245. self.c_self.base.merwe.kappa = kappa
  1246. cdef yaflInt get_np(self, int dim_x):
  1247. return (2 * dim_x + 1)
  1248. cdef const yaflUKFSigmaMethodsSt * get_spm(self):
  1249. return &yafl_ukf_merwe_spm
  1250. #==========================================================================
  1251. #Decorators
  1252. @property
  1253. def alpha(self):
  1254. return self.c_self.base.merwe.alpha
  1255. @alpha.setter
  1256. def alpha(self, value):
  1257. raise AttributeError('MerweSigmaPoints does not support this!')
  1258. @property
  1259. def beta(self):
  1260. return self.c_self.base.merwe.beta
  1261. @beta.setter
  1262. def beta(self, value):
  1263. raise AttributeError('MerweSigmaPoints does not support this!')
  1264. @property
  1265. def kappa(self):
  1266. return self.c_self.base.merwe.kappa
  1267. @kappa.setter
  1268. def kappa(self, value):
  1269. raise AttributeError('MerweSigmaPoints does not support this!')