R/favint.R

In nicolemunoz99/multisensory: Predicting the conditions under which multisensory integration is favored

#' aplot: favorability of integration
#'
#' Plots showing how "favorability of integration" (sensu Munoz &Blumstein)
#'   is calculated. Each plot is of a user-defined parameter area (having dimensions DIM1
#'   and DIM2). Favorability of integration is the proportion of this area in which the
#'   animal integrates.
#'
#' @param zlist A list having length, structure and named components equal to "optint.empty". See
#'   vignette for further requirements of zlist.opt .
#'   COL is Parameter varied between plots.
#'   DIM1 is parameter of dimension-1 of parameter area for calculating favorability of integration.
#'   DIM2 is parameter of dimension-2 of parameter area for calculating favorability of integration.
#'   CONST1 is parameter held constant among plots.
#'   CONST2 is parameter held constant among plots.
#'   CONST3 is parameter held constant among plots.
#'   CONST4 is parameter held constant among plots.
#'   CONST5 is parameter held constant among plots.
#'
#' @return Plot(s) (actually grids produced using function image in the package graphics) of user-specified parameter area (dimensions DIM1 and DIM2) color coded
#'   in order to show which stimuli, if any, the animal uses. The number of plots returned
#'   is equal to length(zlist$COL$value). The "value of integration" of a plot is given by A,
#'   in the title of each plot.
#'
#' @examples
#' aplot(favint.ex1)
#'
#' @export


aplot <- function(zlist){
  # plots "favorability of integration" as a function of 8 parameters (3 of which, are held constant)
  # as specified by user. Requires a list having the structure and named components like "optint.empty"

  if (!is.list(zlist)) stop('zlist must be a list with structure and named components like optint.empty')

  ############################################################
  # associate parameter with plot feature from user input
  ############################################################
  COL <- zlist$COL$param 		# name of parameter varied along x-axis
  DIM1 <- zlist$DIM1$param 	# name of parameter of one dimension of A matrix
  DIM2 <- zlist$DIM2$param 	# name of parameter of the other dimension of A matrix
  CONST1 <- zlist$CONST1$param	# name of parameter held constant
  CONST2 <- zlist$CONST2$param	# name of parameter held constant
  CONST3 <- zlist$CONST3$param	# name of parameter held constant
  CONST4 <- zlist$CONST4$param 		# name of parameter held constant
  CONST5 <- zlist$CONST5$param 		# name of parameter held constant


  if (length(COL) > 5) stop ('Maximum length for COL is 5')
  if (length(CONST4) != 1) stop ('Maximum length for CONST4 is 1')
  if (length(CONST5) != 1) stop ('Maximum length for CONST4 is 1')
  if (length(DIM1) > 100 || length(DIM2) > 100) ("Maximum length for DIM1, DIM2 is 100")
  if (length(CONST1) != 1 || length(CONST2) != 1 || length(CONST3) != 1) stop('CONST1, CONST2,
	  CONST3 must have lengths = 1')

  # Are all parameters included in zlist?
  param.nam = c("bb.b", "u1", "u2", "prior.b", "bb.a", "k1", "k2", "s1")
  tt = unlist(lapply(zlist, function(x) x[[1]]))
  if (any(is.na(match(tt, param.nam)))) stop('zlist must contain name of parameters as in
	  optint.empty (i.e., u1, u2, prior.b, s1, k1, k2, bb.a, bb.b)')

  # Are all parameter values in zlist numeric?
  ss = lapply(zlist, function(x) x[[2]])
  nonnum = sapply(ss, function(x) any(is.numeric(x)))
  if (any(nonnum==F)) {
    exceptions = tt[which(nonnum==F)]
    stop(paste(exceptions, ' has been specified with non-numeric values'))
  }

  ############################################################
  # define parameters based on list from user input
  ############################################################
  param.name <- lapply(zlist, function(x){
    x$param
  })

  # prior probability of world in state b
  prior.b <- unlist(zlist[[which(param.name=="prior.b")]]["value"])


  # fraction of overlap of b and a distributions; assume 10% < overlap < %90
  u1 <- unlist(zlist[[which(param.name=="u1")]]["value"])
  u2 <- unlist(zlist[[which(param.name=="u2")]]["value"])


  # benefit of correct behavior in state "b" (bb.b) and correct behavior in state "a" (bb.a)
  bb.b <- unlist(zlist[[which(param.name=="bb.b")]]["value"])
  bb.a <- unlist(zlist[[which(param.name=="bb.a")]]["value"])

  # cost of using first and second stimulus
  k1 <- unlist(zlist[[which(param.name=="k1")]]["value"])
  k2 <- unlist(zlist[[which(param.name=="k2")]]["value"])

  # magnitude of first stimulus
  s1 <- unlist(zlist[[which(param.name=="s1")]]["value"])

  if(any(prior.b <= 0) || any (prior.b >= 1)) stop ('prior.b must be between 0 and 1')
  if(any(u1 <= 0) || any (u1 >= 1)) stop ('u1 must be between 0 and 1')
  if(any(u2 <= 0) || any (u2 >= 1)) stop ('u2 must be between 0 and 1')
  if(any(bb.b < 0)) stop ('bb.b must be >= 0')
  if(any(bb.a < 0)) stop ('bb.a must be >= 0')
  if(k1 < 0) stop ('k1 must be >= 0')
  if(k2 < 0) stop ('k2 must be >= 0')
  if(any(s1 < -10) || any(s1 > 10) || any(bb.b > 10) || any(bb.a > 10) || any(k1 > 10) || any(k2 > 10))
    warning ('Sensitivity of Favorability of Integration is best visualized for: -10 < s1 < 10;
		0 < bb.a < 10; 0 < bb.b < 10; 0 < k1 < 10; 0 < k2 < 10')


  ############################################################
  # create "param" a table of parameter sets for calculating Value of Information
  ############################################################


  col <- eval(parse(text = COL))
  row <- eval(parse(text = CONST4))
  ser <- eval(parse(text = CONST5))
  xx <- eval(parse(text = DIM1))
  yy <- eval(parse(text = DIM2))
  const1 <- eval(parse(text = CONST1))
  const2 <- eval(parse(text = CONST2))
  const3 <- eval(parse(text = CONST3))

  # x-axis values
  x1 <- sapply(col, function(pp){
    rep(pp, length(row)*length(ser)*length(xx)*length(yy))
  })
  x1.a <- array(x1)

  # values of row-varying parameter
  x1.sub1 <- sapply(row, function(rrr){
    rep(rrr, length(ser)*length(xx)*length(yy))
  })
  x1.sub1.a <- rep(array(x1.sub1), length(col))

  # values of series-varying parameter
  x1.sub2 <- sapply(ser, function(ss){
    rep(ss, length(xx)*length(yy))
  })
  x1.sub2.a <- rep(unlist(x1.sub2), length(col)*length(row))


  # values of dimension 1 of matrix-varying parameter
  x1.sub3.a <- rep(xx, length(col)*length(row)*length(ser)*length(yy))

  # values of dimension 2 of matrix-varying parameter
  x1.sub4 <- sapply(yy, function(yyy){rep(yyy, length(xx))})
  x1.sub4.a <- rep(array(x1.sub4), length(row)*length(col)*length(ser))


  param1 <- cbind(x1.a, x1.sub1.a, x1.sub2.a, x1.sub3.a, x1.sub4.a)
  colnames(param1) <- cn <- sapply(c(COL, CONST4, CONST5, DIM1, DIM2), function(x) x)

  # add values of constants to "param1"
  param <- cbind(param1, rep(const1, nrow(param1)), rep(const2, nrow(param1)), rep(const3, nrow(param1)))
  colnames(param)[6:8] <- c(CONST1, CONST2, CONST3)


  # rownames(param) <- sapply(1:nrow(param), function(x){
  # 	str_c(sapply(cn, function(y){ paste(y, "=", round(param[x,y],2), sep="")}), collapse=" ")
  # })


  ############################################################
  # apply function "voi" to calculate integration results
  # results stored in arrays: vi1.u, vi2.u, int1.u, int2.u, st.use.u
  ############################################################

  vi2.u <- vi1.u <- int2.u <- int1.u <- st.use.u <- array(NA, nrow(param)) # initiate unlisted arrays

  p.list <- vector("list", ncol(param)) # list to hold parameters that will be passed to function "voi"
  names(p.list) <- paste("t.", colnames(param), sep="")


  for(gg in 1:nrow(param)){

    p.list[1:ncol(param)] = param[gg,] # pass parameter values to p.list

    temp = do.call("voi", p.list) # evaluate voi with parameters in p.list

    int1.u[gg] <- temp$int1[1] # is the 1st stimulus used?
    vi1.u[gg] <- temp$vi.1[1] # value of information of 1st stimulus
    int2.u[gg] <- temp$int2[1] # is the 2nd stimulus used?
    vi2.u[gg] <- temp$vi.2[1] # value of information of 2nd stimulus
    st.use.u[gg] <- temp$st.use[1] # 1--> only first stimulus used; 2 --> only 2nd stimulus used; 3 --> integration

    temp = NA # reset
  }



  ############################################################
  # rearrange integration output into lists for plotting
  ############################################################



  skel <- lapply(col, function(a){ # structure for relisting integration arrays back into matrices suitable for graphing
    x = lapply(row, function(b){
      y = lapply(ser, function(c){
        z = matrix(NA, length(xx), length(yy))
        rownames(z) = paste(DIM1, round(xx,2), sep = "=")
        colnames(z) = paste(DIM2, round(yy,2), sep = "=")
        z
      }); names(y) = paste(CONST5, ser, sep = "="); y
    }); names(x) = paste(CONST4, round(row,2), sep = "="); x
  })

  names(skel) <- paste(COL, round(col,2), sep = "=")

  int2 <- relist(int2.u, skel) # execute relisting
  int1 <- relist(int1.u, skel)
  vi2 <- relist(vi2.u, skel)
  vi1 <- relist(vi1.u, skel)
  st.use <- relist(st.use.u, skel)




  ############################################################
  # make plots
  ############################################################

  X.VAL = xx
  Y.VAL = yy

  colors = gray(c(0.2, 0.4, 0.7, 1))
  par(mfrow = c(1,1), cex.axis = 1, mar = c(6, 5, 1, 8), las = 1, lwd = 1, xpd = T)
  for(AA in 1:length(col)){
    for(BB in 1:length(row)){

      image(x = X.VAL, y = Y.VAL, z = st.use[[AA]][[BB]][[1]], xlab = DIM1, ylab = DIM2,
            col = colors, xpd = T, zlim = c(0,3))
      favorab = length(which(st.use[[AA]][[BB]][[1]]==3))/length(st.use[[AA]][[BB]][[1]])
      title(paste("A = ", favorab), cex=1.3)

      #add legend & values of constants
      loc = legend(x = 1.06*max(xx), y = max(yy), rev(c("None", "stim1 only", "stim2 only", "Integration")),
        fill = rev(colors), xpd = NA, cex = 1, title = "stimulus used")
      text(loc$rect$left, 0.5*max(Y.VAL), paste(COL,col[AA], sep=" = "), xpd=NA, cex=1.2, adj=0)
      text(loc$rect$left, 0.4*max(Y.VAL), paste(CONST4, row, sep=" = "), xpd=NA, cex=1.2, adj=0)
      text(loc$rect$left, 0.3*max(Y.VAL), paste(CONST5, ser, sep=" = "), xpd=NA, cex=1.2, adj=0)
      text(loc$rect$left, 0.2*max(Y.VAL), paste(CONST3, const3, sep=" = "), xpd=NA, cex=1.2, adj=0)
      text(loc$rect$left, 0.1*max(Y.VAL), paste(CONST1, const1, sep=" = "), xpd=NA, cex=1.2, adj=0)
      text(loc$rect$left, 0*max(Y.VAL), paste(CONST2, const2, sep=" = "), xpd=NA, cex=1.2, adj=0)

    }
  }

}