psifluo: Fits Kolber et al (1998) model to FL3500 data

Documented in read_psi_file_campbell

#' @title Function to import data from PSI FRRF using Doug Campbell double light/dark FRRf protocol
#' @description Imports a PSI FRRF data file generated by applying an FRRf induction/relaxation on top of a level of actinic light (which may  = 0)  and then again after a user defined dark period (ex. 1 s) to allow re-opening of PSII. Reads in data file, parses metadata, subtracts electronic background from fluorescence signal during light flashlets, splits columns into "light" and "dark" vectors,  splits data by sti and str, calculates the actinic light level at each step.
#' @param filename Path to the file to be imported
#' @param calib_file Path to the calibration file to convert lamp voltage settings to photons. A^-2^ . s^-1^. If not provided the default is a calibration function from a Doug Campbell PSI FL3500 instrument. If you are not using his machine the parameters sigmaPSII, cumulative excitation and actinic PAR values associated to an individual machine will be wrongly calibrated.
#' @param dec Decimal separator, default= "."
#' @return Returns a list of five data frames and one numeric vector. 1 - all.data  = raw data input; 2 - meta.data = acquisiton settings; flashlet_energy = the calculated flashlet energy from the calibration file; 3 - data.sti  = induction only; 4 - data.str  = re-opening only; 5 - measuring_steps = voltage and equivalent actinic PAR at each step.
#' @keywords external
#' @export

read_psi_file_campbell<-function(filename,calib_file=NA,dec="."){
  ####  Read in data ####
  # Note: not all data files are of the same format, I have tried
  # to make this as generic as possible by reading in with scan()
  # and then parsing the data. AB 2014-11-02
  data.raw<-scan(filename,character(0),sep='\n',quiet=TRUE);


#sometimes the computer operating system uses commas insteat of points as a
#decimal separator, the next line replaces commas by points to fix that
  data.raw<-gsub(",",".",data.raw)


  # There are 3 dashed lines:
  # 1) The start of the initial header
  # 2) The end of the inital header, start of acquisiton settings
  # 3) The end of acquisiton settings, start of the data matrix
  # Data start 2 lines after the last dashed separator
  data.sep <- which(data.raw=="-----------------------------------");
  data.start <- grep(pattern="^Selected Data Sets:",x=data.raw);
  data.acq <- grep(pattern="^Data",x=data.raw);
  n.data.sets <- data.sep[2]-data.sep[1]-8; # number of data sets

  # Convert scans into data frame
  # Make text data into a list of numeric values
  tst2<-lapply(strsplit(data.raw[(data.sep[length(data.sep)]+2):length(data.raw)],split='\t'),as.numeric);
  # Take list of numeric values and make a data frame, should have n.data.sets+1 b/c one col is time
  all.data<-data.frame(do.call(rbind,lapply(tst2,matrix,ncol=length(tst2[[1]]),byrow=TRUE)));

  #print(all.data)

  # Get sample times
  tst2<-lapply(strsplit(data.raw[(data.start+1):(data.sep[length(data.sep)-1]-1)],split='\t'),as.character);
  data.samples<-data.frame(do.call(rbind,tst2));


  #extracting gain and offset from the data.samples dataframe, this will be attached
  #to the metadata dataframe


  Gain<-as.character(data.samples$X5[1])
  Gain<-substr(Gain,6,nchar(Gain))
  Offset<-as.character(data.samples$X6[1])
  Offset<-substr(Offset,6,nchar(Offset))



   # Extract Meta-Data
  meta.data.names<-unlist(strsplit(data.raw[data.sep[2]+2],split='\t'),use.names=FALSE);
  tst<-unlist(strsplit(data.raw[data.sep[2]+3],split='\t'));
  meta.data<-as.data.frame(as.list(as.numeric(tst[2:length(tst)])));
  names(meta.data)<-meta.data.names[2:length(meta.data.names)];

  #adding the Gain and Offset to the metadata dataframe
  meta.data$Gain<-as.numeric(Gain)
  meta.data$Offset<-as.numeric(Offset)


  #re-ordering the metadata information to make it a bit more organized
meta.data<-cbind.data.frame(meta.data$ST1_num,meta.data$ST1_light,meta.data$ST1_duration,meta.data$ST1_dark,meta.data$Light_curves_repeats,meta.data$Red_light_curves_enable,meta.data$Blue_light_curves_enable,meta.data$F_Voltage,meta.data$F2_Voltage,meta.data$F1_RED_Voltage,meta.data$F2_BLUE_Voltage,meta.data$Red_light_curves_0,meta.data$Red_light_curves_1,meta.data$Red_light_curves_2,meta.data$Red_light_curves_3,meta.data$Red_light_curves_4,meta.data$Red_light_curves_5,meta.data$Red_light_curves_6,meta.data$Red_light_curves_7,meta.data$Red_light_curves_8,meta.data$Red_light_curves_9,meta.data$Red_light_curves_10,meta.data$Blue_light_curves_0,meta.data$Blue_light_curves_1,meta.data$Blue_light_curves_2,meta.data$Blue_light_curves_3,meta.data$Blue_light_curves_4,meta.data$Blue_light_curves_5,meta.data$Blue_light_curves_6,meta.data$Blue_light_curves_7,meta.data$Blue_light_curves_8,meta.data$Blue_light_curves_9,meta.data$Blue_light_curves_10,meta.data$Light_period,meta.data$Dark_period,meta.data$A_Voltage,meta.data$A2_Voltage,meta.data$M_Voltage,meta.data$A1_preillumination,meta.data$A2_preillumination,meta.data$ActinicFlash,meta.data$MeasuringFlash,meta.data$sample_duration,meta.data$AuxDuration,meta.data$FAR_duration,meta.data$first_AUX,meta.data$second_AUX,meta.data$Logarithmic_Relax1_enable,meta.data$Relax1,meta.data$Relax1_time,meta.data$Relax1_dark,meta.data$Fast_sampling_CH1,meta.data$Fast_sampling_CH2,meta.data$PreFlash,meta.data$init,meta.data$Signal_Averages,meta.data$Stirrer_disable,meta.data$delayed_run,meta.data$lightstep_length,meta.data$SkipDefaultDA,meta.data$MaxDAValue,meta.data$MeasurDelay,meta.data$Gain,meta.data$Offset)

names(meta.data)<-c("ST1_num","ST1_light","ST1_duration","ST1_dark","Light_curves_repeats","Red_light_curves_enable","Blue_light_curves_enable","F_Voltage","F2_Voltage","F1_RED_Voltage","F2_BLUE_Voltage","Red_light_curves_0","Red_light_curves_1","Red_light_curves_2","Red_light_curves_3","Red_light_curves_4","Red_light_curves_5","Red_light_curves_6","Red_light_curves_7","Red_light_curves_8","Red_light_curves_9","Red_light_curves_10","Blue_light_curves_0","Blue_light_curves_1","Blue_light_curves_2","Blue_light_curves_3","Blue_light_curves_4","Blue_light_curves_5","Blue_light_curves_6","Blue_light_curves_7","Blue_light_curves_8","Blue_light_curves_9","Blue_light_curves_10","Light_period","Dark_period","A_Voltage","A2_Voltage","M_Voltage","A1_preillumination","A2_preillumination","ActinicFlash","MeasuringFlash","sample_duration","AuxDuration","FAR_duration","first_AUX","second_AUX","Logarithmic_Relax1_enable","Relax1","Relax1_time","Relax1_dark","Fast_sampling_CH1","Fast_sampling_CH2","PreFlash","init","Signal_Averages","Stirrer_disable","delayed_run","lightstep_length","SkipDefaultDA","MaxDAValue","MeasurDelay","Gain","Offset")



  #Add column to separate sample data ("light") and background ("background")
  all.data$Cond<-rep(c("background","Light"));



  # Subtract background from flashes
  cor.data<-all.data[which(all.data$Cond=="Light"),!names(all.data)=="X1" & !names(all.data)=="Cond"]-all.data[which(all.data$Cond=="background"),!names(all.data)=="X1" & !names(all.data)=="Cond"];

  cor.data<-as.data.frame(cor.data)
  #print(str(cor.data))

 #re-adding sample names
  names(cor.data)<-paste("V",ncol(cor.data):1,sep="");



  # Create time column
  cor.data$time<-all.data$X1[which(all.data$Cond=="Light")] * 1000 * 1000;



  # reorder columns so they are in the proper order, first flash is in last column
  cor.data<-cor.data[ncol(cor.data):1];



  # find duration of dark period between the two induction curves
  dark.period <- meta.data[1,"Dark_period"] *1000 * 1000;



  if(is.na(dark.period)){
    stop('Read function can only samples run with D. Campbell Dark/Light protocol.');
  } else {
    # split data curves
    # For now, assuming only two curves
    halfway<-nrow(cor.data)/2;

    data.light<-cor.data[1:halfway,];  #the first set of rows (with time < 2) describe the first 'read' in each column. put this in a dataframe of their own
    names(data.light)<-gsub("V","light_",names(data.light));
    data.dark<-cor.data[(halfway+1):nrow(cor.data),];  #the second set (with time > 2) describes the second 'read'. Put these in their own dataframe, too.  The second read is usually better (according to Doug, Mitchell), so work with this one
    names(data.dark)<-gsub("V","dark_",names(data.dark));
    data.dark<-data.dark[,!(names(data.dark)=="time")];


    # Combine data frames
    merge.data<-cbind(data.light,data.dark);



    # split data set into STI and STR
    sti.length<-meta.data[1,"ST1_duration"]*1000*1000;
    data.sti<-merge.data[which(merge.data$time<=sti.length),];
    data.str<-merge.data[which(merge.data$time>sti.length),];
    data.str$time<-data.str$time-sti.length;
  }


##############################################################################
#adding the cumulative energy to the sti data.frame
##############################################################################

#importing the calibration file (voltages used vs light values)
#checking if the user has defined a calibration file and if not use Doug Campbell
#file
#ADD disclaimer to inform the user that it's only valid for Doug's machine
if(is.na(calib_file)){
    cali<-psifluo::generic_cali
} else{
  cali<-utils::read.csv(calib_file)
}

#Defining some constants
photons_in_umol <- 6.02214129e17
meter_sq_in_angstrom_sq <- 1.0e-20

#flashlet_energy this needs to came from the calibration file (cali) and from the voltage used
#depends on the color used
if (meta.data$Red_light_curves_enable==1){
  power.level <- meta.data[1,"F_Voltage"]
  flashlet_int <- cali[which(cali[,"Power.Level"]==power.level),"Red.Flash"]
}
if (meta.data$Blue_light_curves_enable==1){
  power.level <- meta.data[1,"F2_Voltage"]
  flashlet_int <- cali[which(cali[,"Power.Level"]==power.level),"Blue.Flash"]
}

#multiply by duration of the flashlet (ST1_light), photons per umol and convert to square angstrom
flashlet_energy <- flashlet_int*meta.data[1,"ST1_light"]*photons_in_umol*meter_sq_in_angstrom_sq

#print(flashlet_int)
#print(flashlet_energy)

#adding the cumulative dose

for (i in 1:nrow(data.sti)){
if (i==1)
{
  data.sti$cumul_energy[i]<-flashlet_energy
}else{
  data.sti$cumul_energy[i]<-data.sti$cumul_energy[i-1]+flashlet_energy
}
}

################################################################################
################################################################################



##############################################################################
#creating a data.frame that will store the light intensities used and associated
#voltages
##############################################################################

#format the time data, the time format is determined by the computer OS
#the first thing is to check if there's an AM or PM string the time.data[,3] object
#it will fail if the time is in any other third type format
#turns out it failed imediatly when we tried it with a greek format
#the temporary solution was to "turn it off" while we think about a better solution
#the time will still be stored in the same objects but it will be stored as a text variable

####NEW CODE
#f1<-scan(filename,character(0),sep='\n',quiet=TRUE)
#data.sep <- which(f1=="-----------------------------------")
#l<-strsplit(f1[(data.sep[1]+9):data.sep[2]-1],split='\t')
#time.data <- data.frame(matrix(unlist(l), nrow=n.data.sets, byrow=T),stringsAsFactors=FALSE)

#the code bellow is redundant and it will be corrected once I have time. BJ (2018-02-25)
a<-grepl("AM|PM", data.samples[,3])
#print(a)
if (a[1]==FALSE){
  #formatted.time <- strptime(time.data[,3], format="%d/%m/%Y %H:%M:%S")
  formatted.time <- data.samples[,3]
}else{
  #formatted.time <- strptime(time.data[,3], format="%m/%d/%Y %I:%M:%S %p")
  formatted.time <- data.samples[,3]
}



#Getting the voltage and respective light intensity, it's dependent on the
#calibration  file

#getting the number of light steps
no_light_steps<-meta.data$Light_curves_repeats
#print(no_light_steps)

#flag to see which colour was used
blue<-meta.data[1,"Blue_light_curves_enable"]
red<-meta.data[1,"Red_light_curves_enable"]

#print(blue)
#print(red)

#if (blue==1){
#  lc_voltage<-names(meta.data)[which(grepl('Blue',names(meta.data)))]
#  lc_voltage<-lc_voltage[which(lc_voltage!='Blue_light_curves_enable')]
#  lc_voltage<-lc_voltage[which(lc_voltage!='Blue_light_curves')]

#  voltages<-data.frame(step=rep(NA,no_light_steps), voltage=NA, par=NA)
#  for(i in 1:(no_light_steps)){
#    voltages$step[i]<-lc_voltage[which(lc_voltage==paste('Blue_light_curves',i-1,sep='_'))]
#    voltages$voltage[i]<-meta.data[1, voltages$step[i]]
    #look up PAR using voltage in calibration file - must already be loaded under 'cali'
#    voltages$par[i]<-cali[cali$Power.Level==voltages$voltage[i], paste('Blue','.Actinic',sep='')]
    #print(voltages)
#  }
#}
#if (red==1){
#  lc_voltage<-names(meta.data)[which(grepl('Red',names(meta.data)))]
#  lc_voltage<-lc_voltage[which(lc_voltage!='Red_light_curves_enable')]
#  lc_voltage<-lc_voltage[which(lc_voltage!='Red_light_curves')]

#  voltages<-data.frame(step=rep(NA,no_light_steps), voltage=NA, par=NA)
#  for(i in 1:(no_light_steps)){
#    voltages$step[i]<-lc_voltage[which(lc_voltage==paste('Red_light_curves',i-1,sep='_'))]
#    voltages$voltage[i]<-meta.data[1, voltages$step[i]]
    #look up PAR using voltage in calibration file - must already be loaded under 'cali'
#    voltages$par[i]<-cali[cali$Power.Level==voltages$voltage[i], paste('Red','.Actinic',sep='')]
    #print(voltages)
#  }
#}


##############################################################################

if (blue==1){
  lc_voltage<-names(meta.data)[which(grepl('Blue',names(meta.data)))]
  lc_voltage<-lc_voltage[which(lc_voltage!='Blue_light_curves_enable')]
  lc_voltage<-lc_voltage[which(lc_voltage!='Blue_light_curves')]

  voltages<-data.frame(step=rep(NA,n.data.sets), voltage=NA, par=NA)

  #extracting the PAR values from the voltages used
  for(i in 1:(no_light_steps)){
    voltages$step[i]<-lc_voltage[which(lc_voltage==paste('Blue_light_curves',i-1,sep='_'))]
    voltages$voltage[i]<-meta.data[1, voltages$step[i]]
    #look up PAR using voltage in calibration file - must already be loaded under 'cali'
    voltages$par[i]<-cali[cali$Power.Level==voltages$voltage[i], paste('Blue','.Actinic',sep='')]
    #print(voltages)
  }
  #filling the rest of the dataframe with the voltage and PAR values. They are
  #always the same
  voltages$step<-rep(voltages$step[1:no_light_steps],length(voltages$step)/no_light_steps)
  voltages$voltage<-rep(voltages$voltage[1:no_light_steps],length(voltages$step)/no_light_steps)
  voltages$par<-rep(voltages$par[1:no_light_steps],length(voltages$step)/no_light_steps)
}
if (red==1){
  lc_voltage<-names(meta.data)[which(grepl('Red',names(meta.data)))]
  lc_voltage<-lc_voltage[which(lc_voltage!='Red_light_curves_enable')]
  lc_voltage<-lc_voltage[which(lc_voltage!='Red_light_curves')]

  voltages<-data.frame(step=rep(NA,n.data.sets), voltage=NA, par=NA)
  for(i in 1:(no_light_steps)){
    voltages$step[i]<-lc_voltage[which(lc_voltage==paste('Red_light_curves',i-1,sep='_'))]
    voltages$voltage[i]<-meta.data[1, voltages$step[i]]
    #look up PAR using voltage in calibration file - must already be loaded under 'cali'
    voltages$par[i]<-cali[cali$Power.Level==voltages$voltage[i], paste('Red','.Actinic',sep='')]
    #print(voltages)
  }
  #filling the rest of the dataframe with the voltage and PAR values. They are
  #always the same
  voltages$step<-rep(voltages$step[1:no_light_steps],length(voltages$step)/no_light_steps)
  voltages$voltage<-rep(voltages$voltage[1:no_light_steps],length(voltages$step)/no_light_steps)
  voltages$par<-rep(voltages$par[1:no_light_steps],length(voltages$step)/no_light_steps)
}




#associating the time stamps to the measuring_steps data.frame
voltages$time<-formatted.time

#print(voltages)
#correcting the name of the 1st column of the all.data dataframe to time
names(all.data)[1]<-'time'

################################################################################
################################################################################



return(list(all.data=all.data,meta.data=meta.data,flashlet_energy=flashlet_energy,data.sti=data.sti,
            data.str=data.str,data.samples=data.samples,measuring_steps=voltages));




}