<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">""" written by Snehal Chokhandre
This script extracts the preconditioning, load unload before and after preconditioning and a combined stress relaxation
data (ramps and waits) and exports to three text files with time, displacement and load for further calculations. No disp
normalization is done. Extracted force data is raw and filtered. peak and equilibrium time, load and displacement values
are also extracted along with ramp times and expected and actual velocities (rates) are  and saved in a csv file.

input needed = data file name """

import sys
import numpy as np
import matplotlib.pyplot as plt
import os
import csv
from scipy.signal import butter, filtfilt
import BasicUtilities


def USAGE(argv):
    print
    print ('USAGE: ') + argv[0] + ' &lt;Insufficient arguments&gt;'
    print

    
def main(argv):

    if len(argv)!= 2:
        USAGE(argv)
        sys.exit(1)

    input_filename = argv[1]
    infile = open(input_filename, 'r')
    
# --------------------------
# pick all data for move relative and wait commands 

    

    ramps = [[] for x in range(15)]
    waits = [[] for x in range(14)]
    sins = []
    vel = []
   
    i = -1 
    f = -1
    sn = -1
    
    for line in infile: # read each line          
         if line[1:5] == 'Move':
             i += 1
             for _ in xrange(6):
               line = infile.next()
               a = line[0]
               while a.isdigit():
                   ramps[i].append(line)
                   line = infile.next()
                   a = line[0]  
         if line[1:5] == 'Wait':
             f += 1
             for _ in xrange(4):
                 line = infile.next()
                 b = line[0]             
                 while b.isdigit():
                     waits[f].append(line)
                     line = infile.next()
                     b = line[0]     
         if line[1:4] == 'Sin':
             sn += 1
             for _ in xrange(8):
                 line = infile.next()
                 b= line[0]             
                 while b.isdigit():
                     sins.append(line)
                     line = infile.next()
                     b = line[0]         
    infile.close()
    
    infile = open(input_filename, 'r') 
    for line in infile: # read each line  
      if line[0:8] == 'Velocity':   # expected rate
             vel.append(line[16:22]) 
             line = infile.next()
    infile.close()
    
    print '----------'
    print 'Expected rate, mm/s =', vel[13]
   
# ------------------------------
# separate time, disp and load from all move relative and wait command data 
    
    r_time = [[] for x in range(15)]
    r_disp = [[] for x in range(15)]
    r_load = [[] for x in range(15)]    
    w_time = [[] for x in range(14)]
    w_disp = [[] for x in range(14)]
    w_load = [[] for x in range(14)]
    sin_disp = []
    sin_load = []
    sin_time = []
    

    for j in range(0,15): 
     for line in ramps[j]:
         line = line.split()
         line = map(float, line)
         r_time[j].append(line[0])
         r_disp[j].append(abs(line[1]))
         r_load[j].append(abs(line[2]))
         line = infile.next
         

    for p in range(0,14): 
     for line in waits[p]:
         line = line.split()
         line = map(float, line)
         w_time[p].append(line[0])
         w_disp[p].append(abs(line[1]))
         w_load[p].append(abs(line[2]))
         line = infile.next 
         
    for line in sins:
         line = line.split()
         line = map(float, line)
         sin_time.append(line[0])
         sin_disp.append(abs(line[1]))
         sin_load.append(abs(line[2]))
         line = infile.next
         
    
#---------------------------------
# filter ramp load data and sinusoid load data. low pass butterworth filter, 3rd order, 20 hz cutoff freq 

         
    # Filter inputs
    order = 3
    fs = 2.5*1000       # sample rate, Hz
    cutoff = 20 # desired cutoff frequency of the filter, Hz
    
   
#    b, a = butter_lowpass(cutoff, fs, order)  # filter coefficients to check its frequency response
    
#    
    filtered_r_load = [[] for x in range(15)]
    filtered_sin_load = []

#    r_load[12] = r_load[12] + w_load[6] # for patella cc, data file from mach 1 has error!! 
#    r_load[12] = r_load[12].tolist()
    
   # Apply the filter
    for r in range(0,15):
          filtered_r_load[r] = BasicUtilities.butter_lowpass_filter(r_load[r], cutoff, fs, order)
          
          
    filtered_sin_load = BasicUtilities.butter_lowpass_filter(sin_load, cutoff, fs, order)
    

#    w_load[6] = BasicUtilities.butter_lowpass_filter(w_load[6], cutoff, fs, order) # for patella cc, data file from mach1 had error!!
#    w_load[6] = w_load[6].tolist()
#---------------------------------
# find peak/instantaneous and equilibrium load, disp and time (from filtered data, not normalized by sample thickness/length) easier to pick from separated ramps than combined stress-relaxation data 
    
  
    inst_disp = []
    inst_load = [] 
    inst_time_stamp = []
    ramp_low_limit = []
    ramp_up_limit = []
    
    for r in range(11,14): # range depends on the number of move relative commands
           y = 0
           while r_disp[r][y+1]-r_disp[r][y] == 0:
                 y += 1
          
           ramp_low_limit.append(y)  # will be used to find ramp data for rate and duration calculation
           
#           peak = max(r_disp[r])
#           kx = [iv for iv, jg in enumerate(r_disp[r]) if jg == peak] # find all indices at which peak value is present in list
#           yj = kx[0] # pick 1st index at which max disp appears 

           min_ramp = r_disp[r+1][0]
           yj = BasicUtilities.find_nearest(r_disp[r], min_ramp)
           
           
           ramp_up_limit.append(yj)   # will be used to find ramp data for rate and duration calculation
           inst_time_stamp.append(yj) 
           inst_disp.append(r_disp[r][yj])
           inst_load.append(filtered_r_load[r][yj])

    
    eqbm_disp = []
    eqbm_load = [] 
    eqbm_time_stamp = []
    for t in range(12,15):  # range depends on the number of wait commands
         g = 0
         while r_disp[t][g+1]-r_disp[t][g] == 0:
               g += 1
         eqbm_time_stamp.append(g)
         eqbm_disp.append(r_disp[t][g])
         eqbm_load.append(filtered_r_load[t][g])
         

    print 'inst_disp =', inst_disp
    print 'inst_load =', inst_load
    print 'eqbm_disp =', eqbm_disp
    print 'eqbm_load =', eqbm_load
    
# -------------------------------   
# strain rate and ramp duration calculation 
    
    actual_ramp = [[] for x in range(3)]
    actual_time = [[] for x in range(3)]
    for r in range(11,14):
        actual_ramp[r-11].append(r_disp[r][ramp_low_limit[r-11]])
        actual_time[r-11].append(r_time[r][ramp_low_limit[r-11]])
        xy = ramp_low_limit[r-11]+1
        while xy != ramp_up_limit[r-11]:
              actual_ramp[r-11].append(r_disp[r][xy])
              actual_time[r-11].append(r_time[r][xy])
              xy += 1

   
    applied_rates = []
    for r in range(0,3):
        (applied_rate,applied_inter) = np.polyfit(actual_time[r],actual_ramp[r],1)
        applied_rates.append(applied_rate)

    
    print 'Applied rates at each ramp, mm/s =', applied_rates
    
    time_taken = []
    for r in range(0,3):
        end = r_time[r+2][ramp_up_limit[r]]
        start = r_time[r+2][ramp_low_limit[r]]
        time_taken.append(end-start)
            
    print 'Time taken for each ramp, s =', time_taken   
    print '----------'

#-------------------------------
# combine all stress - relaxation load, time and disp  values (unfiltered and filtered load)  for plotting and saving 

    load_stress_relax = []  # unfiltered load combined for stress relaxation 
    filtered_load_stress_relax = [] # filtered load combined for stress relaxation 
    disp_stress_relax = []  # displacement combined for stress relaxation 
    time_stress_relax = []  # time combined for stress relaxation
    
    z = 0
    for h in range(11,14):
          filtered_r_load[h] = filtered_r_load[h].tolist()  # convert array to list 
          filtered_load_stress_relax = filtered_load_stress_relax + filtered_r_load[h]
          disp_stress_relax = disp_stress_relax + r_disp[h]
          load_stress_relax = load_stress_relax + r_load[h]
          z += 1
          k = z + 4
          for z in range(z,k):
                filtered_load_stress_relax = filtered_load_stress_relax + w_load[z+1]
                disp_stress_relax = disp_stress_relax + w_disp[z+1]  
                load_stress_relax = load_stress_relax + w_load[z+1]
   
   
# time resets to zero at the beginning of every command on mach 1. to combine it in one list, time from every following command has to be added to that from the last 
    z = 0
    u = 0

    for h in range(11,14):  
           r_time[h] = [x + u for x in r_time[h]]
           time_stress_relax = time_stress_relax + r_time[h]
           z += 1
           k = z + 4
           for z in range(z,k):
               u = time_stress_relax[len(time_stress_relax) - 1]
               w_time[z+1] = [x+u for x in w_time[z+1]]
               time_stress_relax = time_stress_relax + w_time[z+1]
           u = time_stress_relax[len(time_stress_relax) - 1]
    

# -------------------------------    
# save combined ramp and hold to text file, (both filtered and unfiltered), text file name = folder name + _ramp_hold 

    file_path = os.path.abspath(input_filename)  
    folder_name = os.path.basename(os.path.normpath(file_path[:-8])) # grab folder name 
    
    ramp_hold_file = folder_name + '_ramp_hold'
    fd_data = np.array([time_stress_relax, disp_stress_relax, load_stress_relax, filtered_load_stress_relax])
    fd_data = fd_data.T
    datafile_id = open(ramp_hold_file+'.txt', 'w+')
    np.savetxt(datafile_id, fd_data)
    datafile_id.close()
    
# -------------------------------
# save preconditioning data to text file,(both filtered and unfiltered), text file name = folder name + _pc 
    
    pc_file=folder_name + '_pc'
    pc_fd_data = np.array([sin_time, filtered_sin_load, sin_disp])
    pc_fd_data = pc_fd_data.T
    pc_datafile_id = open(pc_file+'.txt', 'w+')
    np.savetxt(pc_datafile_id, pc_fd_data)
    pc_datafile_id.close()
    
# -------------------------------
# save ramp load unload before and after pc to text file (both filtered and unfiltered), text file name = folder name + _pc_ramps

    pc_ramps_file = folder_name + '_pc_ramps'
    pc_ramps_fd_data = np.array([r_disp[2], r_load[2], filtered_r_load[2], r_disp[3], r_load[3], filtered_r_load[3], r_disp[6],r_load[6],filtered_r_load[6], r_disp[7], r_load[7], filtered_r_load[7]])
    pc_ramps_fd_data = pc_ramps_fd_data.T
    pc_ramps_datafile_id = open(pc_ramps_file+'.txt', 'w+')
    np.savetxt(pc_ramps_datafile_id, pc_ramps_fd_data)
    pc_ramps_datafile_id.close()
    
# -------------------------------   
# write outputs to csv file peak and eqbm load and disp, desired and actual rates, time taken for each ramp. 
         
    output_file = folder_name + '_output'
    out_datafile_id = open(output_file+'.csv', 'w')
    out_writer = csv.writer(out_datafile_id)
    out_writer.writerow(['Expected rate, mm/s', vel[13] ])
    out_writer.writerow([ 'Applied rates at each ramp, mm/s', applied_rates ])
    out_writer.writerow([ 'Time taken for each ramp, s', time_taken ])   
    out_writer.writerow(['inst_disp', inst_disp ])
    out_writer.writerow([ 'inst_load',inst_load ])
    out_writer.writerow([ 'eqbm_disp',eqbm_disp ])
    out_writer.writerow([ 'eqbm_load',eqbm_load ])   
    out_writer.writerow([ 'Experimental reference location sr, 10g',r_disp[9][0] ])
    out_writer.writerow([ 'Experimental reference location pc, 10g',r_disp[0][0] ])
# ------------------------------
# PLOT DATA (to check the extracted data quality, not saved) 
    
    plt.figure(1)
    ax1 = plt.subplot2grid((2,1), (0,0))
    ax2 = plt.subplot2grid((2,1), (1,0))    
    ax1.plot(time_stress_relax,disp_stress_relax,'r')
    ax1.set_ylabel('disp,mm')
    ax1.set_xlabel('time,s')    
    ax2.plot(time_stress_relax,load_stress_relax,'-ro')
    ax2.plot(time_stress_relax,filtered_load_stress_relax,'-go')   
    ax2.set_ylabel('load,gf')
    ax2.set_xlabel('time,s')
    
    plt.figure(2)
    plt.plot(disp_stress_relax, filtered_load_stress_relax)
    plt.plot(inst_disp,inst_load,'go')
    plt.plot(eqbm_disp,eqbm_load,'ro')
    plt.xlabel('disp,mm')
    plt.ylabel('load,gf')
    
    plt.figure(3)
    plt.plot(sin_disp,filtered_sin_load)
    plt.xlabel('disp,mm')
    plt.ylabel('load,gf')
    
    plt.figure(4)
    plt.plot(sin_time,filtered_sin_load)
    plt.xlabel('time,s')
    plt.ylabel('load,gf')
    
    plt.figure(5)
    plt.plot(r_disp[2],filtered_r_load[2])
    plt.plot(r_disp[3],filtered_r_load[3])
    plt.plot(r_disp[6],filtered_r_load[6])
    plt.plot(r_disp[7],filtered_r_load[7])
    plt.xlabel('disp,mm')
    plt.ylabel('load,gf')
    
    
#    plt.figure(6)
#    plt.plot(w_disp[6], w_load[6],'.')
    
    
    plt.show()
    
# ------------------------------


if __name__=="__main__":
    main(sys.argv)




</pre></body></html>