From 91484012deade99502a8641426820c084a1aacd9 Mon Sep 17 00:00:00 2001
From: Raymond Chia <rqchia@janus0.ihpc.uts.edu.au>
Date: Sat, 18 Nov 2023 15:05:58 +1100
Subject: [PATCH] hernandez dsp

---
 modules/.datapipeline.py.swp | Bin 20480 -> 0 bytes
 regress_rr.py                | 137 +++++++++++++++++++++++++++--------
 2 files changed, 108 insertions(+), 29 deletions(-)
 delete mode 100644 modules/.datapipeline.py.swp

diff --git a/modules/.datapipeline.py.swp b/modules/.datapipeline.py.swp
deleted file mode 100644
index 4213d68c350303ec5029eb85195c3fa423f0fba4..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001

literal 20480
zcmeHPYm6jS6)upcvMQkb5zS@Rrn_f*YPxsXuxZnTMOJqL!($&H+dHY$RNb1Xnd+*d
zs(NN;Hy|Vi|EVM(YJ&P3ksrR6XiP-JM~Fe9fq)W;gkU5D@P~w#_?`QxuI}mC9W*9B
zs*-Q&)_t6F?z!hy-9G1J=9@Rn%(7i()4=sX!}#4J2l=mW+;h(h&+Igko1I19riK$=
zOLo{#o2I`QJ7zyiP44zhyT5+St&-q$oJ339$<k>%@$IREr?H>dnIDByjvY+*BDWv#
zblSBuJN9E9_#rprwe{r6no<L$1~#jKe&$W>+H!Em_LkV{TPH4J=bv|Qv!ly?r3Ok3
zlo}{CP->vmK&gRJ1EmH^4g3#iAdAm7K8tyrs^)UDx}P#~KTSP{>Yj|W|NR|##j5_<
zk@mM$d#dWMQw%hpw^aMQYJYg7{Y}+=ShZJ_z<T@}s{N2^zkQ_rb=B^v_TP`Rzem+;
zdb?GHepb!pt<*rNfl>pd21*T-8YneTYM|6Wsew`hr3Ok3{8wtgb`9fPNd7AQ&_(=T
z|NejW1BP)N_!e*va13yOi-9e`bMH5d9|Ly*i@*YKJ+K3~0QfL)D)6^0hVd8R8Q{0T
zFM#`i1>h4v8<+-6zyRJk2mOHKzyrV?z>UD=zy$F6`wZhT;LE@~Z~^e>d(jv8Ja7QG
z6xaeheKuqSz6LA+9&i;<0p2{zFdhf)1#Sf*AOsEoQ$Pc_0C@3C!}t#H72q?#Wxy%G
zYiAh7gTTGODliLd1ztYgFn$9(1pEZ}25=a-9B2Ve;9TH2I8^=w+y(Rj7nlJUumv~=
zco{y3p97Bo4**{UmVr577jP=@3j7#90=@`*0oV(i2he)HGR}Fya5+)vvJP9|nMFmP
zaxc=F7ur2uH>kt9YmwV7Pdd#83zKLy?TEh(HaWSxY9|Y+9OZf-8oOqK&X$v|NQw=X
zF8i?`f{PubvF-5Afuzr6wID|ARA^xOg-ch1Xkv#8yw)0w+16SwXraF!W;G$Al||Vk
zHCQ6UWOxt{Qr;e7N!)WX+d6D>bbN0x$WK`qW$c<LERbH$NC=0Coty!P`>B;BwjWAa
z6*oB{b{T872Wrdb4AO)0R4$&F?1tis*21t3Yl$i-r9iF2(tfYjVqTOm%FV$Xm?ssx
z>Sv2Aig~DPiP<Uhq||!Cq?nv25!P33N^QgCP-q<Zu+DVF>hWwbawWxfnkCW*<n%jD
zb5|qR&!}AyFJ)5+&-w}EBoj9U?OMv?dUvWdryeJ)iF&s=Czy}h7{t1+RIR85^MOI)
z;|_iI)-2a|vM%i?3!Udgo2o_ovFdEK&9H;2`*4RXSFdiOu$>Am6n0QyI|`RjVFwC3
zZxPZ6f>~q*_B;=2n*11fQtk+zTNz5=n7Jm*-wrH1La4MRS7|HBT)cX|lSHad$A#!-
zp=N_IL*3;=!yuw#sE7TGr<yQ5=$;w9g)+a#T^?k1Z87R6svi|)FVIFOAr<vZGrJe7
zi4Z$#pmsas)NJB0rb3PWK#`T+bj6*n+{7w-+m+dNWxrYRs&Y=lyK>Wna!NTFY^bVP
z75g50ffx_oD(`~qj1-Vn;v!`tM`p@89Yz*p(5&g&G-WY!MUwVDc5k=Zu68=kqS(Z^
zDY573#HA&7V9(C*_9WX8x*W0nT;X04^+-m!Fz8*E`<FX3T6e>4b4Dnj#-MrpKUk}n
zC{32G6E2+PhYjZSGHIf60HX>R`9$TWO0VJ`sO+uGR`ws7Fmy9=Bw>Pws(KwkWMC)M
z>0e-p&4Vp8w0|bp2<Pd43+J1xO3$jeQ<c3|W!9?fU&q$^XW7E^=*ZI`jZRg1Qx%t0
zX3Wa0S=o=>H!-Yx(jOT(?F1U7URls_7dD?{D7bj-pj4~3Io|z+=WQ+d6(Ab=Qe8kM
z!Rw0KiL+Qss<T$LA=&uhdj2`1kx_@;W0%~r(Ll0^19rMxW%yE8XZLNatqrEPvpsdm
zKz>PXnXf0qk2pNe*bR2j=X;VQN(RRFE)7rE`T9k{JDb);nq8KxN&2I-v&gmA@M)mq
zgMRN6-(Xl6aaQotsPictMLE5Y==DYCVb5P+elL!aj9t0sii4kSumdx*d#q1hd-eVX
zyK?t|-PV;e`wT;MaaOjgq3LpfB@9g`=%*P^G>uuDcdxKBXYs08ZilI4qx&MmnbSu(
zPzI%liZh;tPi1C_6H0^aU+jBc!1vQxNRe>-c+K=9%_oR#7cXip*HRjVx}0!bLWfIt
z=m$S`=XH4@i1KpS@5O6$U<(n#6<V>=j49E;pC>p-@D)3djb1-=vM35t-6h5bzyV*j
zW;eS%-9TkUvDXhWI_f(-P5m%e_;sMuU{2h}uaSPJW*VhtY~whoB=A$$r=uUKf*;9H
z;<4m4(s3!m#|+b;!4^G??ZuDyX_n^w<fH7Ft_|iC*cpzwJxxLj!^0Inq4!5`SxD8h
zmpKY)0|(+XlY(Y8gj}_QWkWA757Q{|f@qc4W}dEhu|-YTtL-H?oKl<X>bBSp(JaIq
z)U9Plr9j%qwOIODzDUHGN)v(!#6gq={yb9YsU&CynL)q*e~91nPXqM(U;kvJ=clNy
zyh;s}8YneTYM|6Wsew`hr3Ok3lo}{CP->vmK&gTMU=3)0Bi+d_WV-P1zz-jIru28g
z&s-P3{}JXk{-EOjp_<Y2`0f8B@C0x_a0Iv(FoE-d(}B}~A0WnmCy)Ra0~Z2+MjZbI
z;4$EvKn5%Udw`Du7XmLKmj6@Wd%$;r+kivBG;k(x2JjMM`#%Hj2R;Yf2z&%M3pf*a
z3NikNfI9(-{qF;=2U@_Ji1EJ){0jIba3`<|aKHsxz{i2-G5LpqF9BU(H*gVf4)7Y{
z{7(Rn0}lenfv*F11M@%~cmpy2hk)+`$ACS+7T|5f_}>Da1daoD01mJdr~(%QFCxZ&
zAFu`-0j>c)1pF0o{l@`{`F|U@9k>m+9H5y0DZrzM@qZb(8=#o~#{gQxKWYaK){$1;
zPE!t_X0*bSz>Y01uooIEPAKj?p(+uZUtr{RWbigI;e{mssjd9r@cU9kyGWtU$cqT)
zH6la_0{SW*SL^~ubr_pI+<Z0Z*CdLT>tT=s{}een3pI7q4&l1(NWZL<e7FxBvlg6>
zeV*3l*E&^Df!o^;<3Q&VNT>M-uT6{+S}O8_<j7&b7IEyXR%`&jg)mAuIr#Yz5o{jA
zba*(0?e#>*j+g9*yA_#&+)+QSO=i8gF7-mIiYMY@_UyPup4bL+J&SygQXmnBC$DNP
zcdFLvb>Z?GxaNtDx>D(9R$g&~YKbbbkX~K2hC(pqo0et&DDUJY8L%B%g5jlRpxB*j
z0CFlyXw;G)+k*L7X^7bh%b@9uW|(YcGa=(PB9@WIvPI53AO9+(O_?m>F(9Ay_GX2q
zO3fr1Bd0#%HhE-U3Miy@y&|uCB1kJhSiX>vAXK8F{;nsUNJ%>3MH$P;dI&|;M^^-X
zc1yXyd0Jc(C_{{ZHn4N-B(c|IMa}h)b~7|gNHA1|QH8Y@bH&=(L4ercQ7%kS_?HDS
zTHKLnEU8>jrGgYe&WeW!t=DR{;|L-Rwy>6nQoX@!c+oqBOdQ=L$0KOvEZ2@;yL78)
ziBbotA=NokZeytzD3NK>>qnslK@vo0IKmGjBW2KOtT_`=af$s=j9t``5J7f@-C(Fb
z`FtdMq4)JdAI`h$Va5BV_zX|8$Y*4ll9vc`Q7DU&QT$%t5&;?^*~rpqYpF$jwf}Qr
zQDfRM5^`K&=iOC$ccpWK$<r*b&Y2(f2WdHY<Ga4+S<d3}K)<?1sR^&^9s(BwO7w&_
z>UJ3tb^<pSf~wFW5*efb-Bx>Uj!g|!HQPvX7(^R0B+6v1;bC*g&yleWJrPQb$oW8Q
zMP8*#G}~+*xoFFx2c;|dh{y;Hh(I~wQZ){(BC-VwUrASC(IQ`_gR+?}+A#-1D4XGe
zCLs)7TzW}Xp~+~KyJ9b>Y!li$4TjMmL@MaCn$3~LlCNn&ax30UF*^$TW-0cle3yw9
z*$Mm*Fv_@Fqd1*Oe`=mZw+2GT4svDu6`4y1YZl`mCzCxeHV=U1gQ}BxVo5Ra?F1vd
z>M;>|PA&SR+@i-r(}il}UlnrRP9%tqqarWoIv+u(NqobrlbY&qKF%1a2E)_&rQ%4E
zXMSl|4!utFg%F>{T0vx$M>fCIp&YT1B~o9dd?KV(RoC^Ahl2+b(mP$9DcG**(l(Li
zRv*|Tjn%1dbmpqHrYQeLBqOHv`eabwHpsD~&YBX~P=DDGW{92!wd;#A&H^ObxhnOJ
z>XnJ(tAr5fz_OL>%8;Qe8mA<ne4p)w1VYu9%tyS`Wbz$;DU3Z`<lnJQ;oDRErpe(%
zO`RijYN*QUlfnAokH2KqmY;tqW~bXUksjAHTeKYwn(Oe#n9h}(IP1~Ck#Dw|@Wy8|
zIYvw$^XOtYsafh6gBN%>SZ(=?Zgegir@CPiA11xA!l~cjNOw}DV};$wD#qFi3w*K_
zDI9wY;yWzyy3yjqQb^MrO3ox(mU)_^RB6NJjl=tvn2x^+=P}yc!fM3%Zzzdbtdz2z
z0S~dyaY`Q6i%b|DZ&2MTzKb$7x=Nd#o(vL<3V0ruBRafA^H~g%r}HTYO4=0lrZ`Ey
UwD|PszGKp2bK)Cpk+C!S2T0GB<NyEw

diff --git a/regress_rr.py b/regress_rr.py
index 38d987f..a5376e4 100644
--- a/regress_rr.py
+++ b/regress_rr.py
@@ -22,6 +22,7 @@ from itertools import repeat, chain, combinations
 from multiprocessing import Pool, cpu_count
 import tensorflow as tf
 
+from sklearn.decomposition import PCA
 from sklearn.preprocessing import StandardScaler, MinMaxScaler, OneHotEncoder
 from sklearn.preprocessing import PolynomialFeatures, LabelEncoder
 from sklearn.model_selection import KFold, train_test_split
@@ -35,6 +36,7 @@ from modules.datapipeline import get_file_list, load_and_snip, load_data, \
         load_split_data, load_harness_data
 from modules.digitalsignalprocessing import vectorized_slide_win as vsw
 from modules.digitalsignalprocessing import imu_signal_processing
+from modules.digitalsignalprocessing import hernandez_sp, reject_artefact
 from modules.digitalsignalprocessing import do_pad_fft,\
         pressure_signal_processing, infer_frequency
 from modules.utils import *
@@ -64,7 +66,8 @@ from sktime.transformations.panel.rocket import (
         MiniRocketMultivariateVariable,
 )
 
-from config import WINDOW_SIZE, WINDOW_SHIFT, IMU_FS, DATA_DIR, BR_FS
+from config import WINDOW_SIZE, WINDOW_SHIFT, IMU_FS, DATA_DIR, BR_FS\
+        , FS_RESAMPLE
 
 IMU_COLS =  ['acc_x', 'acc_y', 'acc_z', 'gyr_x', 'gyr_y', 'gyr_z']
 
@@ -453,7 +456,6 @@ def load_and_sync_xsens(subject):
 
     return xsens_df
 
-
 def load_tsfresh(subject, project_dir,
                  window_size=12, window_shift=0.2, fs=IMU_FS,
                  overwrite=False):
@@ -556,6 +558,55 @@ def get_test_data(cal_df, activity_df, xsens_df, test_standing):
 
     return pd.DataFrame(activity_list)
 
+def dsp_win_func(w_inds, df, i, cols):
+    time = df['sec'].values
+    if w_inds[-1] == 0: return
+    w_df = df.iloc[w_inds]
+    t0, t1 = time[w_inds][0], time[w_inds][-1]
+    diff = time[w_inds[1:]] - time[w_inds[0:-1]]
+    mask = np.abs(diff)>20
+    diff_chk = np.any(mask)
+    if diff_chk:
+        return
+
+    # cols = ['acc_x', 'acc_y', 'acc_z',
+    #         'gyr_x', 'gyr_y', 'gyr_z']
+
+    data = w_df[cols].values
+
+    if reject_artefact((data-np.mean(data,axis=0))/np.std(data,axis=0)):
+        return
+
+    # DSP
+    pca = PCA(n_components=1, random_state=3)
+
+    # do hernandez sp on datacols for df
+    filt = hernandez_sp(data=data, fs=IMU_FS)[1]
+
+    # pca
+    pca_out = pca.fit_transform(filt)
+
+    std = StandardScaler().fit_transform(pca_out)
+
+    pred = get_max_frequency(std, fs=FS_RESAMPLE)
+
+    # get pss / br estimates
+    # x_time median, pss max_freq, br median
+    sm_out = w_df['BR'].values
+    ps_out = w_df['PSS'].values
+
+    x_vec_time = np.median(time[w_inds])
+
+    fs = 1/np.mean(diff)
+    ps_freq = int(get_max_frequency(ps_out, fs=IMU_FS))
+
+    y_tmp = np.array([x_vec_time, np.nanmedian(sm_out), ps_freq])
+
+    y_hat = pd.DataFrame([ {'sec': x_vec_time, 'pred': pred} ])
+    y_out = pd.DataFrame([y_tmp], columns=['sec', 'br', 'pss'])
+
+    return y_hat, y_out
+
 # save evaluation metrics in single file that handles the models for the
 # subject and config
 class EvalHandler():
@@ -616,18 +667,21 @@ def imu_rr_model(subject,
                  feature_method='tsfresh',
                  train_len:int=3,
                  test_standing=False,
+                 data_input:str='imu+bvp',
                 ):
     # window_size, window_shift, intra, inter
     cal_str = 'cpm'
     fs = IMU_FS
     tmp = []
     imu_cols = ['acc_x', 'acc_y', 'acc_z', 'gyro_x', 'gyro_y', 'gyro_z']
-    bvp_cols= ['bvp']
+    bvp_cols = ['bvp']
 
-    # TODO: 
-        # implement and input args config by data cols
-        # implement and input args config with test_standing
-    data_cols = imu_cols + bvp_cols
+    if 'imu' in data_input and 'bvp' in data_input:
+        data_cols = imu_cols + bvp_cols
+    elif 'imu' in data_input and not 'bvp' in data_input:
+        data_cols = imu_cols
+    elif not 'imu' in data_input and 'bvp' in data_input:
+        data_cols = bvp_cols
 
     do_minirocket = False
     use_tsfresh   = False
@@ -680,7 +734,37 @@ def imu_rr_model(subject,
     # include standing or not
     test_df_tmp = get_test_data(cal_df, activity_df, xsens_df, test_standing)
     test_df = pd.concat([df for df in test_df_tmp['data']], axis=0)
-    
+    x_test_df, y_test_df = get_df_windows(test_df, 
+                                          imu_df_win_task,
+                                          window_size=window_size,
+                                          window_shift=window_shift,
+                                          fs=fs,
+                                         )
+
+    acc_dsp_df, acc_y_dsp_df =  get_df_windows(test_df, dsp_win_func,
+                                               window_size=window_size, 
+                                               window_shift=window_shift,
+                                               fs=fs,
+                                               cols=['acc_x', 'acc_y', 'acc_z'])
+    gyr_dsp_df, gyr_y_dsp_df =  get_df_windows(test_df, dsp_win_func,
+                                               window_size=window_size, 
+                                               window_shift=window_shift,
+                                               fs=fs,
+                                               cols=['gyr_x', 'gyr_y', 'gyr_z'])
+
+    acc_evals = Evaluation(acc_y_dsp_df[lbl_str], acc_dsp_df['pred'])
+    gyr_evals = Evaluation(gyr_y_dsp_df[lbl_str], gyr_dsp_df['pred'])
+    print("acc evals: \n", acc_evals.get_evals())
+    print("gyr evals: \n", gyr_evals.get_evals())
+    plt.subplot(211)
+    plt.plot(acc_y_dsp_df[lbl_str]); plt.plot(acc_dsp_df['pred'])
+    plt.subplot(212)
+    plt.plot(acc_y_dsp_df[lbl_str]); plt.plot(acc_dsp_df['pred'])
+    plt.show()
+
+    # TODO implement evals from dsp to results
+    ipdb.set_trace()
+
     for combi in combinations(cal_df[cal_str].values, train_len):
         combi_str = "-".join([str(x) for x in combi])
         pfh.config[cal_str] = combi_str
@@ -707,17 +791,12 @@ def imu_rr_model(subject,
                                                     window_shift=window_shift,
                                                     fs=fs,
                                                    )
-            x_test_df, y_test_df = get_df_windows(test_df, 
-                                                  imu_df_win_task,
-                                                  window_size=window_size,
-                                                  window_shift=window_shift,
-                                                  fs=fs,
-                                                 )
 
             x_train = make_windows_from_id(x_train_df, imu_cols)
-            x_test  = make_windows_from_id(x_test_df, imu_cols)
             y_train = y_train_df[lbl_str].values.reshape(-1, 1)
+            x_test  = make_windows_from_id(x_test_df, imu_cols)
             y_test  = y_test_df[lbl_str].values.reshape(-1, 1)
+    
 
             print("minirocket transforming...")
             x_train = np.swapaxes(x_train, 1, 2)
@@ -731,19 +810,6 @@ def imu_rr_model(subject,
             x_test  = test_df.iloc[:, 3:].values
             y_test  = test_df[lbl_str].values.reshape(-1, 1)
         else:
-            x_train_df, y_train_df = get_df_windows(train_df,
-                                                    imu_df_win_task,
-                                                    window_size=window_size,
-                                                    window_shift=window_shift,
-                                                    fs=fs,
-                                                   )
-            x_test_df, y_test_df = get_df_windows(test_df, 
-                                                  imu_df_win_task,
-                                                  window_size=window_size,
-                                                  window_shift=window_shift,
-                                                  fs=fs,
-                                                 )
-
             x_train = make_windows_from_id(x_train_df, imu_cols)
             x_test  = make_windows_from_id(x_test_df, imu_cols)
             y_train = y_train_df[lbl_str].values.reshape(-1, 1)
@@ -814,6 +880,15 @@ def arg_parser():
                         default=3,
                         help='minutes of data to use for calibration'
                        )
+    parser.add_argument('-d', '--data_input', type=str,
+                        default='imu',
+                        help='imu, bvp, imu+bvp: select data cols for input'
+                       )
+    parser.add_argument('-ts', '--test_standing', type=int,
+                        default=0,
+                        help='1 or 0 input, choose if standing data will be '\
+                        'recorded or not'
+                       )
     args = parser.parse_args()
     return args
 
@@ -832,6 +907,8 @@ if __name__ == '__main__':
     lbl_str        = args.lbl_str
     train_len      = args.train_len
     overwrite      = args.overwrite
+    data_input     = args.data_input
+    test_standing  = args.test_standing
 
     print(args)
     assert train_len>0,"--train_len must be an integer greater than 0"
@@ -855,7 +932,9 @@ if __name__ == '__main__':
                               mdl_str=mdl_str,
                               overwrite=overwrite,
                               feature_method=feature_method,
-                              train_len=train_len
+                              train_len=train_len,
+                              data_input=data_input,
+                              test_standing=test_standing,
                              )
 
         if mdl_str in ['fnn', 'lstm', 'cnn1d', 'elastic', 'ard', 'xgboost']:
-- 
GitLab