minmax_rejection

make_subject_trialwise.R

library("dplyr")
library(matrixStats)
source("two_force_sensor_functions.R")

pre_process_xpd = function(xpd) {
    # preprocessing xpd-files:
    #
    # adding variable
    #     small
    #     switch:  mini-block switch 
    #     count: row counter
    #     ref_time: reference time for the mini-block (first stim presentation)
    
    #find switches
    xpd$count = c(1:nrow(xpd))
    xpd$small = as.numeric(xpd$number<5)
    
    next_small = c(xpd$small[2:nrow(xpd)], NA)
    switch = abs(xpd$small - next_small)
    xpd$switch = c(1, switch[1:length(switch)-1])
    
    xpd$ref_row = NA
    reference_row = NA
    for (i in c(1:nrow(xpd))) {
        if (xpd[i, "switch"] == 1) {
            reference_row = xpd[i, "count"]
        }
        xpd[i, "ref_row"] = reference_row
    }
    
    return(xpd)
}

minmax_rejection =function(data, minmax_difference, n_samples) {
    # returns the id of first sample of the first period in which a 
    # a specified difference between mininum and maximum occured 
    #
    # nsamples: the period to be test for each sample
    
    # make a matrix: in each row the data copied and shifted by one 
    # sample. nrow = n_samples
    len = length(data)
    mtx = matrix(nrow = n_samples, ncol = len ) 
    for (r in 1:n_samples) {
        mtx[r, 1:(len+1-r)] = data[r:len]
    }    
    # test columwise min max difference
    found = which(colMaxs(mtx, na.rm = T) - colMins(mtx, na.rm = T) > minmax_difference)
    if (length(found) ==0)  return(NA)
    else return(min(found))
}

data_path = "~/raw_data/16-number_force_two_sensors/data"

# loading and preparing data

subject_id = c(c(1:21))
for ( id in subject_id ) {
    filename_in = paste0(data_path, "/number_force_two_sensors_", id, ".Rdata")
    filename_out = paste0("data/subject_trialwise", id, ".Rdata")
    
    print(filename_in)
    load(filename_in)
    
    d = extract.forces(force1.df = data$force1, force2.df = data$force2, 
                       zero_times = data$xpd$time_udp, before = 200, after=1000)
    
    
    xpd = pre_process_xpd(data$xpd)
    
    
    # correcting for baseline --> force_std
    baseline1 = rowMeans( d$force1[, 180:200], na.rm=T)
    baseline2 = rowMeans( d$force2[, 180:200], na.rm=T)
    force1_std = d$force1 - baseline1
    force2_std = d$force2 - baseline2
    
    is_good_function = function(d) {
        is.na(apply(d, 1, minmax_rejection, minmax_difference = 200, n_samples = 100))
    }
    good_samples1 = is_good_function(force1_std)
    good_samples2 = is_good_function(force2_std)
    
    n_good_trials = sum(good_samples1 & good_samples2, na.rm = T)
    print(n_good_trials)
    
    xpd$good_trials = good_samples1 & good_samples2 & (xpd$trial >=0) & (xpd$oddball==0)
    #copy mtx twice for left & right hand
    xpd$hand = "l"
    tmp = xpd
    tmp$hand = "r"
    
    design = rbind(xpd, tmp)
    data = rbind(force1_std, force2_std)
    
    save(design, data, file = filename_out)
}