flow = c(rep(NA,10),NAins(as.data.frame(runif(5000)), .1)$runif)
    rain = runif (length(flow))
    event = with(rle(!is.na(flow )),cbind(length=lengths[values],position=cumsum(c(1,lengths))[values]));

test_function = function(ndays, event, flow, rain,upboundary){
      flowvolume = rainvolume = raininweek = raininmonth =NULL;
      for (i in 1:(length(event)/2)){
        if (upboundary < event[,'position'][i]){
          flowvolume[i] = sum(flow[(event[,'position'][i]):(event[,'position'][i]+event[,'length'][i]-1)], na.rm = TRUE) # total flow during the non NA period
          rainvolume[i] = sum(rain[(event[,'position'][i]):(event[,'position'][i]+event[,'length'][i]-1)], na.rm = TRUE) # total rain during the non NA period
          raininweek[i] = sum(rain[(event[,'position'][i]-ndays[1]):(event[,'position'][i]-1)], na.rm = TRUE) #total rain imediate before NA with a constrained period of nday[1]
          raininmonth[i] = sum(rain[(event[,'position'][i]-ndays[2]-ndays[1]):(event[,'position'][i]-ndays[1]-1)], na.rm = TRUE) #total rain iprior to nday[1]
        } else {next}
      }
      -summary(lm(flowvolume ~ rainvolume + raininweek + raininmonth))$r.squared # to minimise R2
    }

results= optim(par=c(2,20), lower=c(1,10), upper=c(10,30),method="L-BFGS-B",test_function, event=event, rain=rain, flow=flow,upboundary=30)

################################################################
    # RANDOMLY INSERT A CERTAIN PROPORTION OF NAs INTO A DATAFRAME #
    ################################################################
    NAins <-  NAinsert <- function(df, prop){
      n <- nrow(df)
      m <- ncol(df)
      num.to.na <- ceiling(prop*n*m)
      id <- sample(0:(m*n-1), num.to.na, replace = FALSE)
      rows <- id %/% m + 1
      cols <- id %% m + 1
      sapply(seq(num.to.na), function(x){
        df[rows[x], cols[x]] <<- NA
      }
      )
      return(df)
    }