Skip to content

Commit

Permalink
adjust yaml vector format
Browse files Browse the repository at this point in the history
  • Loading branch information
aremazeilles committed Jan 20, 2021
1 parent 1130d59 commit 80209b5
Show file tree
Hide file tree
Showing 14 changed files with 1,279 additions and 1,282 deletions.
45 changes: 22 additions & 23 deletions src/pi_FPE/pi_FPE/footplacementestimator.py
Original file line number Diff line number Diff line change
Expand Up @@ -30,16 +30,16 @@ def calcFPE(joint,com,angMom,comIR,events,mass,trlinfo):
omega = np.ones([t.shape[0],3])*np.nan
# get angular momentum matrix
H = np.vstack([angMom.x,angMom.y,angMom.z]).transpose()
# rearrange IR and calculate omega
# rearrange IR and calculate omega
for i in range(comIR.shape[0]):
# get inertial tensor and reshape
# get inertial tensor and reshape
IR = np.vstack([np.hstack([comIR.xx[i],comIR.xy[i],comIR.xz[i]]),
np.hstack([comIR.yx[i],comIR.yy[i],comIR.yz[i]]),
np.hstack([comIR.zx[i],comIR.zy[i],comIR.zz[i]])])
# multiply the angular momentum with the inverse of the inertial tensor
omega[i,:] = np.matmul(np.transpose(H[i,:]),np.linalg.inv(IR))

l_foot = pd.concat([(joint.l_TTII_x+joint.l_heel_x)/2,(joint.l_TTII_y+joint.l_heel_y)/2,(joint.l_TTII_z+joint.l_heel_z)/2],axis=1)
l_foot = pd.concat([(joint.l_TTII_x+joint.l_heel_x)/2,(joint.l_TTII_y+joint.l_heel_y)/2,(joint.l_TTII_z+joint.l_heel_z)/2],axis=1)
walkdir = np.argmax(np.ptp(np.array(l_foot),axis=0))

# Bring the data into the direction of the angular momentum
Expand All @@ -56,7 +56,7 @@ def calcFPE(joint,com,angMom,comIR,events,mass,trlinfo):
yaxis = np.cross(zaxis,xaxis)
# normalize axes
yaxis = yaxis/norm_col(yaxis)
zaxis = zaxis/norm_col(zaxis)
zaxis = zaxis/norm_col(zaxis)
xaxis = xaxis/norm_col(xaxis)
R_proj = np.hstack([xaxis,yaxis,zaxis])
omega_proj = prod_col(transpose_col(R_proj),omega)
Expand All @@ -73,7 +73,7 @@ def calcFPE(joint,com,angMom,comIR,events,mass,trlinfo):
anglearr[nonwalkdir[walkdir]]
angTM=np.tile(np.radians(anglearr),[vcom.shape[0],1])
theta_proj = prod_col(transpose_col(R_proj),angTM)

# create empty matrices
lFPE = np.zeros([vcom.shape[0],1])*np.nan
phi = np.zeros([vcom.shape[0],1])*np.nan
Expand All @@ -99,7 +99,7 @@ def calcFPE(joint,com,angMom,comIR,events,mass,trlinfo):
Jcom[i] = np.dot(np.dot(yaxis[i,:],IR),yaxis[i,:])
if vcom_proj[i,0] == vcom_proj[i,0]:
phi[i] = optimize.fmin(fFPE,phi_in,args = (Jcom[i,0],com_proj[i,2],g,mass,omega_proj[i,1],theta_proj[i,1],vcom_proj[i,0],vcom_proj[i,1]),disp=0)
# phi[i] = optimize.brentq(fFPE,0,3,args = (Jcom[i,0],com_proj[i,2],g,mass,omega_proj[i,1],theta_proj[i,1],vcom_proj[i,0],vcom_proj[i,1]),xtol=1e-6,rtol=1e-7)
# phi[i] = optimize.brentq(fFPE,0,3,args = (Jcom[i,0],com_proj[i,2],g,mass,omega_proj[i,1],theta_proj[i,1],vcom_proj[i,0],vcom_proj[i,1]),xtol=1e-6,rtol=1e-7)
phi_in = phi[i]
lFPE[i] = ((np.cos(theta_proj[i,1])*(com_arr[i,1]))/(np.cos(phi[i])))*np.sin(theta_proj[i,1]+phi[i])
sys.stdout.write("]\n") # end progress bar
Expand All @@ -113,7 +113,7 @@ def calcFPE(joint,com,angMom,comIR,events,mass,trlinfo):
rhsFPE = np.abs(rankle-com_arr)-np.abs(gFPE-com_arr)
dfpe = np.vstack((lhsFPE[lhs.astype(int),:],rhsFPE[rhs.astype(int),:]))
return dfpe

def fFPE(phi,Jcom,hcom,g,m,omega,theta,vx,vy):
# =============================================================================
# Foot placement estimator formula to optimze
Expand All @@ -127,10 +127,10 @@ def fFPE(phi,Jcom,hcom,g,m,omega,theta,vx,vy):
t9 = hcom**2
fpefun = (Jcom/2+(m*t5*t8*t9)/2)*((Jcom*omega+m*t3*t7*hcom*(vx*np.cos(t4)+vy*np.sin(t4)))**2)*1/((Jcom+m*t5*t8*t9)**2)-g*m*t3*t7*hcom+g*m*t3*t7*hcom*np.cos(t4)
return fpefun**2

def store_result(file_out, value):
file = open(file_out, 'w')
file.write('---\ntype: \'vector\'\nvalues:')
file.write('---\ntype: \'vector\'\nvalue:')
for line in value:
file.write('\n '+format(line[0], '.5f'))
file.write('\n')
Expand All @@ -139,9 +139,9 @@ def store_result(file_out, value):

def store_result2(file_out, value):
file = open(file_out, 'w')
file.write('---\ntype: \'vector\'\nvalues:')
file.write('---\ntype: \'vector\'\nvalue:')
for line in value:
file.write('\n '+format(line, '.5f'))
file.write('\n - '+format(line, '.5f'))
file.write('\n')
file.close()
return True
Expand Down Expand Up @@ -169,9 +169,9 @@ def main():
file_in_comIR = sys.argv[4]
file_in_events = sys.argv[5]
file_in_trialinfo = sys.argv[6]
folder_out = sys.argv[7]
folder_out = sys.argv[7]


# check input parameters are good
if not os.path.exists(file_in_joint):
print(colored("Input file {} does not exist".format(file_in_joint), "red"))
Expand Down Expand Up @@ -209,7 +209,7 @@ def main():
if not os.path.isfile(file_in_trialinfo):
print(colored("Input path {} is not a file".format(file_in_trialinfo), "red"))
return -1


if not os.path.exists(folder_out):
print(colored(
Expand All @@ -235,7 +235,7 @@ def main():
if not os.path.isfile(file_in_trialinfo):
print(colored("Output path {} is not a folder".format(file_in_trialinfo), "red"))
return -1

# load events structure
events = read_events(file_in_events)
# load joint data
Expand Down Expand Up @@ -266,7 +266,7 @@ class strct():
mass = 81
# calculate margins of stability
l_mos_ap,l_mos_ml,r_mos_ap,r_mos_ml,apdir,mldir= calcMOS(joint,com,events,trlinfo)
# Estimate the foot placement using the foot placement estimator
# Estimate the foot placement using the foot placement estimator
dfpe = calcFPE(joint,com,angmom,comIR,events,mass,trlinfo)

file_out0 = folder_out + "/pi_FPE_ap.yaml"
Expand All @@ -275,40 +275,39 @@ class strct():
print(colored(
"Foot Placement Estimator: vector with size {}x1, stored in {}".format(int(dfpe.shape[0]),file_out0),
"green"))

file_out1 = folder_out + "/pi_FPE_ml.yaml"
if not store_result2(file_out1,dfpe[:,mldir]):
return -1
print(colored(
"Foot Placement Estimator: vector with size {}x1, stored in {}".format(int(dfpe.shape[0]),file_out1),
"green"))

file_out2 = folder_out + "/pi_l_MOS_ap.yaml"
if not store_result2(file_out2, l_mos_ap):
return -1
print(colored(
"Margins of stability [ap]: vector with size {}x1, stored in {}".format(l_mos_ap.shape[0],file_out2),
"green"))

file_out3 = folder_out + "/pi_r_MOS_ap.yaml"
if not store_result2(file_out3, r_mos_ap):
return -1
print(colored(
"Margins of stability [ap]: vector with size {}x1, stored in {}".format(r_mos_ml.shape[0],file_out3),
"green"))

file_out4 = folder_out + "/pi_l_MOS_ml.yaml"
if not store_result2(file_out4, l_mos_ml):
return -1
print(colored(
"Margins of stability [ap]: vector with size {}x1, stored in {}".format(l_mos_ml.shape[0],file_out4),
"green"))

file_out5 = folder_out + "/pi_r_MOS_ml.yaml"
if not store_result2(file_out5, r_mos_ml):
return -1
print(colored(
"Margins of stability [ap]: vector with size {}x1, stored in {}".format(r_mos_ml.shape[0],file_out5),
"green"))
return 0

17 changes: 8 additions & 9 deletions src/pi_FPM/pi_FPM/footplacementmodel.py
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ def timenormalize_step(signal,start,stop):
del stop[0]
print('Delete first sample of stop')
if start.shape[0]>stop.shape[0]:
start = start[0:stop.shape[0]]
start = start[0:stop.shape[0]]
print('remove trailing part start')
if stop.shape[0]>start.shape[0]:
stop = start[0:start.shape[0]]
Expand All @@ -49,7 +49,7 @@ def timenormalize_step(signal,start,stop):
cycle[:,i] = fsig_int(np.linspace(0,50,51))
sigoi = None
return cycle


def calcFPM(joint,com,events,fs):
# =============================================================================
Expand All @@ -60,7 +60,7 @@ def calcFPM(joint,com,events,fs):
rto = np.round(np.array(events.r_toe_off)*fs)
rhs = np.round(np.array(events.r_heel_strike)*fs)
lto = np.round(np.array(events.l_toe_off)*fs)

meve = np.array([])
arr_ieve = np.hstack((lhs,rto,rhs,lto))
arr_neve = np.hstack((np.ones(lhs.shape[0])*0,np.ones(rto.shape[0]),np.ones(rhs.shape[0])*2,np.ones(lto.shape[0])*3))
Expand All @@ -74,7 +74,7 @@ def calcFPM(joint,com,events,fs):
rto = meve[0,np.where(meve[1,:]==1)]
rhs = meve[0,np.where(meve[1,:]==2)]
lto = meve[0,np.where(meve[1,:]==3)]

# get foot position from joint
lfoot = np.array(pd.concat([(joint.l_TTII_x+joint.l_heel_x)/2,(joint.l_TTII_y+joint.l_heel_y)/2,(joint.l_TTII_z+joint.l_heel_z)/2],axis=1)/1000)
rfoot = np.array(pd.concat([(joint.r_TTII_x+joint.r_heel_x)/2,(joint.r_TTII_y+joint.r_heel_y)/2,(joint.r_TTII_z+joint.r_heel_z)/2],axis=1)/1000)
Expand Down Expand Up @@ -211,7 +211,7 @@ def main():
if not os.path.isfile(file_in_events):
print(colored("Output path {} is not a folder".format(file_in_joint), "red"))
return -1

# load events structure
events = read_events(file_in_events)
# load joint data
Expand All @@ -222,22 +222,21 @@ def main():
fs = 1/np.mean(np.diff(np.array(com.time)))
# calculate the spatiotemporal parameters
FPMl_est,FPMr_est,rsq_l,rsq_r,com_l = calcFPM(joint,com,events,fs)

file_out0 = folder_out + "/pi_l_FPM.yaml"
if not store_result(file_out0, rsq_l[24]):
return -1
print(colored(
"r-squared left leg: {} stored in {}".format(round(rsq_l[24],2), file_out0),
"green"))

file_out1 = folder_out + "/pi_r_FPM.yaml"
if not store_result(file_out1, rsq_r[24]):
return -1
print(colored(
"r-squared right leg: {} stored in {}".format(round(rsq_r[24],2), file_out1),
"green"))
return 0




Expand All @@ -248,4 +247,4 @@ def main():




25 changes: 12 additions & 13 deletions src/pi_LDE/pi_LDE/localdivergence.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,15 +18,15 @@ def calcStateSpace(sig,events,fs,n_dim,delay):
# =============================================================================
from scipy.interpolate import interp1d
print("creating state space...")
# make sorted heel strike list
# make sorted heel strike list
lhs = np.array(events.l_heel_strike)*fs
rhs = np.array(events.r_heel_strike)*fs
meve = np.array([])
arr_ieve = np.hstack((rhs,lhs))
arr_neve = np.hstack((np.zeros(rhs.shape),np.ones(lhs.shape)))
meve = np.vstack((arr_ieve,arr_neve))
meve = meve[:,meve[0,:].argsort()]
# get number of strides and samples
# get number of strides and samples
n_samples = (meve.shape[1]-1)*100
# new signal without trailinabs
idx_truesamp = int(meve[0,0]-1)
Expand All @@ -39,7 +39,7 @@ def calcStateSpace(sig,events,fs,n_dim,delay):
# sligthly different than Sjoerd's method (does not extrapolate data)
fsig_int = interp1d(t_new,sig_new,kind = 'cubic')
sig_int = fsig_int(t_inter)
# create empty statespace array
# create empty statespace array
state = np.zeros([n_samples-(n_dim-1)*delay,n_dim])
#create the state space now:
for i_dim in range(n_dim):
Expand Down Expand Up @@ -80,16 +80,16 @@ def calcLDE(state,ws,fs,period,nnbs):
for i_nn in range(nnbs):
div_tmp = np.subtract(state[i_t:i_t+int(win),:],state[idx_sort[i_nn]:(idx_sort[i_nn]+int(win)),:])
divmat[i_t*nnbs+i_nn,:] = np.sum(div_tmp**2,1)**.5 # track divergence, and store
sys.stdout.write("]\n") # end progress bar
sys.stdout.write("]\n") # end progress bar
divergence = np.nanmean(np.log(divmat),0)
xdiv = np.linspace(1,divergence.shape[0],divergence.shape[0])/fs
xs = xdiv[1:int(np.floor(L1))]
Ps = np.polynomial.polynomial.polyfit(xs,divergence[1:int(np.floor(L1))],1)
Ps = np.polynomial.polynomial.polyfit(xs,divergence[1:int(np.floor(L1))],1)

lde = np.ones([1,2])*np.nan
lde[0] = Ps[1]
L2 = np.round(4*period*fs)

L2 = np.round(4*period*fs)
if L2 < win:
idxl = (np.linspace(0,win-L2-1,int(win-L2))+int(L2))
Pl = np.polynomial.Polynomial.fit(idxl/fs,divergence[int(idxl[0]):int(idxl[-1])],1)
Expand All @@ -116,9 +116,9 @@ def main():

file_in_com = sys.argv[1]
file_in_events = sys.argv[2]
folder_out = sys.argv[3]
folder_out = sys.argv[3]


# check input parameters are good
if not os.path.exists(file_in_com):
print(colored("Input file {} does not exist".format(file_in_com), "red"))
Expand All @@ -144,7 +144,7 @@ def main():
if not os.path.isfile(file_in_events):
print(colored("Output path {} is not a folder".format(file_in_com), "red"))
return -1

# load joint data
com = pd.read_csv(file_in_com)
# load events structure
Expand All @@ -169,7 +169,6 @@ def main():
return 0






4 changes: 2 additions & 2 deletions src/pi_SpatTemp/pi_SpatTemp/SpatTemp.py
Original file line number Diff line number Diff line change
Expand Up @@ -65,9 +65,9 @@ class strct():

def store_result(file_out, value):
file = open(file_out, 'w')
file.write('---\ntype: \'vector\'\nvalues:')
file.write('---\ntype: \'vector\'\nvalue:')
for line in value:
file.write('\n '+format(line, '.5f'))
file.write('\n - '+format(line, '.5f'))
file.write('\n')
file.close()
return True
Expand Down
Loading

0 comments on commit 80209b5

Please sign in to comment.