  # single ABM run type 2 with split&merge starting with synchronic production data in str.df

  # this script reads synchronic production data from str.df, initializes agents' memories, and
  # performs a full ABM run.
  # A vast number of variables control the ABM behavior and the logging/plotting during the ABM
  # (see the first part of the master.R script for an explanation of these variables).
  # The script does not perform any before/after analysis, except that it calculates the agent-averaged
  # cluster variability before and after the ABM to evaluate whether the ABM run was plausible 
  # with regard to stability: if the variability increases or is reduced more than 20% a warning 
  # is printed (even if printMessages = F).
  

  # get number of agents (number of speakers) from the Speaker column
  nrOfAgents = length(unique(str.df$Speaker))
  # initialise memory
  # returns an array of lists (one list per Speaker); each list contains a 
  # number of arrays, each as long as there are tokens for that Speaker in str.df.
  # The number of (parallel) arrays can be changed by changing the function cmem.*();
  # cmem.timedelay() reads from its argument (dataframe with initial speaker data):
  # columns named "P..." (e.g. P1 P2 P3) : features exchanged in this ABM, 
  #   here: DCT-0...2 coefficients of M1 tracks (this could be different number of dims!);
  #   the number of these columns define the dimensionality of the ABM feature space;
  #   note that all plot functions are dim-dependent; if you change the dim, check all plot functions!
  # Word : word label
  # V : phonological class label that is defined by the word and does not change (initially only /s/ or /S/)
  # Age : age of token (a number, higher value = younger), this is overwritten by random values in cmem.timedelay()
  # Vpn : the speaker label
  # Initial : the canonical label, here: sibilant type depedent of Word: 's', 'str' or 'S'
  initMemory = cmem.timedelay(str.df)
  # create population (basically a copy, so that the initMemory stays intact during the ABM)
  myPop = createPopulation(nrOfAgents, initMemory, list(ratiosOfTotal=list()))

  # if doSplitAndMergeBeforeABM is set, do an exhaustive split&merge on all agents' data to see how the 
  # synchronic data distribution influences individual phonological cluster without any 
  # iteractions between agents; note that this does not alter the 'beforeABM' plots, since these
  # show the (un-altered) data points/tracks, but in the log file (and the plot over cluster
  # distributions) we will see the initial phonological clustering at sG = 0.
  if(doSplitAndMergeBeforeABM == T) {
    for(j in 1:nrOfAgents) {
      myPop[[j]] = splitandmergefull(myPop[[j]])
    }
  }
  # Flo 20170303 : Is a single splitandmerge enough?
  # Assume that the initial /s/ cluster consists of 20 different words/contexts, then 
  # the single split allows only one separation of this cluster, although the optimal 
  # clustering might anything between 2 and 10 (2 words per cluster are minimum).
  # These clusterings would then occur during the ABM, but *not* due to interactions but rather 
  # due to the initial data configuration, thus masquerading as change through interactions.
  # We can easily test this by repeating the split until nothing changes, and then the 
  # merge until nothing changes (see function 
  # splitandmergefull() in file splitAndMerge.R)

 
  # make a single data-frame o out of the synchroneous data.
  # this function simply converts the array of lists into a dataframe for easier plotting
  # Later the script will create a log dataframe
  # res2.df which at the beginning contains only 'o' marked with o$Cond == 'orig' to indicate 
  # 'original memory data'; later the script will 
  # add new sections to 'res2.df' marked with res2.df$Cond = <number of past interactions>. 
  # E.g. all lines of res.df with o$Cond="40" represent the agents's memories after 40 simGroups. 
  o = myPoptodf(myPop, "orig")
 
  sG = 0  # SimGroup is zero (no interaction yet)

  # write header for cluster config table in abmParameters; this is just a summary of 
  # equivalence phone clusters after each simGroup; a more detailed protocol (which agent 
  # does what) is in file phonClusterTable_<agentName>.df (see below)
  cat("Number of agents having equivalence labels after each simGroup\n\nABM simGroup ",phonClusterLog.names,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  # store the initial cluster configuration (before interactions)
  cat("\"ABM\" \"SimGroup\" \"Vpn\" \"cGroup\"\n",file=paste(LogDirDate,"/phonClusterTable_",abmRun,".df",sep=""))
  if(doSplitAndMergeBeforeABM == T) {
    # log phone clusters of each agent into phonClusterTable_<agentName>.df
    cgroups = get.cgroups(o)
    if(abmRun!="") {
      ABM = rep(as.character(abmRun),n=nrow(cgroups))
    } else { 
      ABM = rep("ALL",times=nrow(cgroups))
    }
    SimGroup = rep(sG,times=nrow(cgroups))
    CGroups.df = data.frame(ABM=ABM,SimGroup=SimGroup,cgroups)
    write.table(CGroups.df,file=paste(LogDirDate,"/phonClusterTable_",abmRun,".df",sep=""),append=T,col.names=F,row.names=F)
    # make statistics output of equivalence phone clusters
    if(printMessages) cat("ABM: ",abmRun,", Sim Group ",sG," finished (",sG*simGroupSize,"interactions)\n")
    cLine = calc.clusterStat(o)
    if(printMessages) print(cLine)
    write.table(cLine,file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T,col.names=F,row.names=F)
    # add equivalence cluster occurance counters to log file for summarized plot later
    phonClusterLog = rbind(phonClusterLog,cLine)
  }  

  # calculate the 'articulation space' for sibilants based on the population
  # this is a cube in the 3D DCT space which borders are at the 99%-quantile 
  # of the observed data. In this ABM speakers cannot produce sibilants outside this space
  # streched by the factor articulStretch (default: 99.0). The default is to set this 
  # parameter articulStretch to such a high value that there is effectively no 
  # restriction of production; if articulStretch is set to say 2.5 agents cannot produce 
  # tokens which are very far outliers.
  # In this experiment this feature did not turn out to be useful.
  P1q = quantile(o$P1,c(0.01,0.99))
  P2q = quantile(o$P2,c(0.01,0.99))
  P3q = quantile(o$P3,c(0.01,0.99))
  articulSpace = as.matrix(cbind(P1q,P2q,P3q))
  # stretch lower/upper quantiles
  articulSpace[1,] = articulSpace[1,] - articulStretch*abs(articulSpace[2,]-articulSpace[1,])
  articulSpace[2,] = articulSpace[2,] + articulStretch*abs(articulSpace[2,]-articulSpace[1,])
  # if you want to limit plots during the ABM to this 'articulation space', remove the following comments
  #xlim = articulSpace[,1]
  #ylim = articulSpace[,2]
  #zlim = articulSpace[,3]
 
  #########################
  # make some background plot layers beforehand because they take a lot of processing time 
  # and are the same for all iteration groups, and make a snapshot of the data before any interaction 
  # with frame number '0000'
  # for subsequent plotting filter initial tokens to one speaker
  # this just a feature to allow speaker-specific plots during the ABM, which is useful
  # to observe the memory contents of a single agent; if parameter plotSpeaker="", 
  # memory content of all agents are plotted together
  if(plotSpeaker != "") {
    o.plot = o[o$Vpn==plotSpeaker,]
  } else {
    o.plot = o
  }
  # for plotting filter to a single group (here: the sibilant type Initial)
  # this just another plot feature, if you want to observe the distribution of 
  # just one group during ABM; if plotOnlySibilant="", all groups are plotted in different colors
  if(plotOnlySibilant != "") {
    o.plot = o.plot[o.plot$Initial==plotOnlySibilant,]
  }
  if(dispOrig||printPNG) {
    # make plot layer of initial data for reconstructed M1 tracks (takes time)
    if(plotTracks) {
      # aggregate: first reconstruct all tracks in o from DCTs, then 
      # aggregate for mean and sd (to plot confidence interval)
      o.dctTracks.df = build_dct_df.sibilant(o.plot)
      # rename F2Pop as M1smooth
      names(o.dctTracks.df)[1] = "M1smooth"
      # calculate the mean at each time point
      # calculate the SD at each time point, separately per speaker and per Initial class
      if(plotSpeakerAverage) {
        m1m = aggregate(M1smooth ~ times * Initial + Cond, mean, data = o.dctTracks.df)
        m1sd = aggregate(M1smooth ~ times * Initial + Cond, sd, data = o.dctTracks.df)
      } else {
        m1m = aggregate(M1smooth ~ times * Vpn * Initial + Cond, mean, data = o.dctTracks.df)
        m1sd = aggregate(M1smooth ~ times * Vpn * Initial + Cond, sd, data = o.dctTracks.df)
      }     
      # this is for ±1.96 standard deviations, analogous to what they do in SS-ANOVA
      k = 1.96
      # the lower limit
      lower = m1m$M1smooth - k * m1sd$M1smooth
      # the upper limit
      upper = m1m$M1smooth + k * m1sd$M1smooth
      # plot orig data tracks as solid lines
      # group = Initial works here because we have only one track per Vpn and per Initial!
      # /S/ : blue, /s/ : red, /str/ : green
      ylim=c(-2.5,2.5)
      if(plotOnlySibilant=="") {
        col=strongColors[c(1,5,7)]
      } else {
        if(plotOnlySibilant=="s") { col=strongColors[5] }
        if(plotOnlySibilant=="S") { col=strongColors[1] }
        if(plotOnlySibilant=="str") { col=strongColors[7] }
      }
      # plot the mean with a thick solid line
      # plot the upper and lower confidence in this solid lines
      if(!plotSpeakerAverage) {
        plot.orig.m = xyplot(M1smooth ~ times | Vpn, group = Initial, type = c("l", "g"),
          data=m1m,ylim=ylim,col=col,lty=1,lwd=3)
        plot.orig.u = xyplot(upper ~ times | Vpn, group = Initial, type = c("l", "g"),
          data=m1sd,ylim=ylim,col=col,lty=1,lwd=1)
        plot.orig.l = xyplot(lower ~ times | Vpn, group = Initial, type = c("l", "g"),
          data=m1sd,ylim=ylim,col=col,lty=1,lwd=1)
      } else {
        plot.orig.m = xyplot(M1smooth ~ times | Initial, group = Initial, type = c("l", "g"),
          data=m1m,ylim=ylim,col=col,lty=1,lwd=3)
        plot.orig.u = xyplot(upper ~ times | Initial, group = Initial, type = c("l", "g"),
          data=m1sd,ylim=ylim,col=col,lty=1,lwd=1)
        plot.orig.l = xyplot(lower ~ times | Initial, group = Initial, type = c("l", "g"),
        data=m1sd,ylim=ylim,col=col,lty=1,lwd=1)
      }
      plot.orig = plot.orig.m + plot.orig.u + plot.orig.l
    }
    # make the first animation frame labelled '0000' with initial data
    # decide color grouping...
    sGnum = "0000"
    if(colorGrouping=="equivalence") {
      # ... according equivalence labels
      eclassLabels = as.character(o.plot$Initial)  # make vector as long as necessary
      for(spkr in unique(o.plot$Vpn)) {
        temp = o.plot$Vpn == spkr
        eclassLabels[temp] = equal_class(as.character(o.plot$Initial[temp]), as.character(o.plot$V[temp]))
      }
      colorGroupVectorOrig = factor(eclassLabels,levels=phonClusterLog.names)
    } else {
      # ... canonical labels
      colorGroupVectorOrig = factor(o.plot$Initial,levels=phonClusterLog.names) 
    }
    mainString = paste("interactions: 0000, speaker: ",abmRun,", colors: ",colorGrouping," labels")
    if(plotCloud) {
      # iteration 3D plot 
      Claspect=c(1,1)
      cloudOverlaySim = cloud(P3 ~ P1 * P2, data=o.plot,group=colorGroupVector,col=strongColors,scales = list(arrows = FALSE),xlim=xlim,ylim=ylim,zlim=zlim,main=paste("DCT-0...2 space, ",mainString),xlab=xname,ylab=yname,zlab=zname,key=generalLegend)
      if(dispSim) {
        print(cloudOverlaySim)
      }
      if(printPNG) {
        PNGname = paste(dirPNG,"Sim0000.png",sep="/")
        png(PNGname,1024,768)
        print(cloudOverlaySim)
        dev.off()
      }
    } else if(plotTracks) {
      if(dispSim) { 
        print(plot.orig)
      }
      if(printPNG) {
        PNGname = paste(dirPNG,"Sim0000.png",sep="/")
        png(PNGname,1024,768)
        print(plot.orig)
        dev.off()
      }
    } else {
     # default sim plot
     # plot DCT-0 vs DCT-2 only since DCT-1 is not distinguishable
     if(plotSpeakerAverage) {
       # layer with actual data colored 
       xyOverlay02 = xyplot(P3 ~ P1,group=colorGroupVectorOrig,type='p',col=strongColors,data=o.plot,xlim=xlim,ylim=zlim,main=paste("DCT-0 - DCT-2 , ",mainString),xlab=xname,ylab=zname,key=generalLegend)
     } else {
       # layer with actual data colored 
       xyOverlay02 = xyplot(P3 ~ P1|Vpn,group=colorGroupVectorOrig,type='p',col=strongColors,data=o.plot,xlim=xlim,ylim=zlim,main=paste("DCT-0 - DCT-2 , ",mainString),xlab=xname,ylab=zname,key=generalLegend)
     }
     if(dispSim) {
       print(xyOverlay02)
     }
     if(printPNG) {
       PNGname = paste(dirPNG,"Sim0000.png",sep="/")
       png(PNGname,1024,768)
       print(xyOverlay02)
       dev.off()
     }
    }
  }
  ########################
  # end making initial data plots
  ########################
 

  ########################
  # Start ABM simulation loop
  # the loop goes over simGroups where each simGroup models simGroupSize interactions
  ######################## 
  # collect orig distribution and the following sim group distributions in 
  # the log dataframe res2.df; column ABM marks the ABM run (in single run mode = plotSpeaker or empty)
  res2.df = cbind(o,ABM=rep(abmRun,times=nrow(o)))
  sGnum = 1000 # just for naming the animated png pictures (starting with ...)
  # Loop
  for(sG in 1:simGroups) {
    ##################
    # Do the interaction over one sim group with simGroupSize interactions
    ##################
    # This is the core call that actually simulates simGroupSize random interactions between agents
    # and changes their memory accordingly (MyPop is overwritten by the new memory!); 
    # as you can see, the function performMultipleInteraction() receives two 
    # function names (character strings) that define which function is to be used for
    # receiving (and memorizing) a token : percFuncName
    # producing a token : prodFuncName
    # Both functions must be defined in a file in sub dir 'functions' that has been 
    # sourced in Rcmd/pathsAndLibraries.R (at the very beginning of this script).
    # In this example the function 'prod.sibilant.2()' and 'perc.sibilant.2()' are 
    # defined in the file functions/sibilantagentfunctions.R
    # Note that these functions *must match* the memory structure ofthe agents; e.g.
    # if you decide you need another item in the agent's memory (also see cmem....()),
    # then these basic functions have to be adapted to deal with these new data.
    myPop = performMultipleInteraction(myPop, prodFuncName=prodFuncName, percFuncName=percFuncName, nrOfInteractions=simGroupSize)
    # convert new memories of all agents into dataframe and log it into res2.df; marked with the simulation group number in column Cond 
    o1 = myPoptodf(myPop, as.character(sG))
    res2.df = rbind(res2.df,cbind(o1,ABM=rep(abmRun,times=nrow(o1))))
    # for plotting filter current memory tokens to one speaker
    if(plotSpeaker != "") {
      o1.plot = o1[o1$Vpn==plotSpeaker,]
    } else {
      o1.plot = o1
    }
    # for plotting filter current memory tokens to a single sibilant class (Initial)
    if(plotOnlySibilant != "") {
      o1.plot = o1.plot[o1.plot$Initial==plotOnlySibilant,]
    }

    ################
    # make simulation plots of the current memory tokens etc. after the simGroup
    # Aside from the plots to screen (dispSim) we additionally print animation frames to 
    # dirPNG (printPNG) that are later aggregated into an animated GIF;
    # data points of the initial ABM data are printed in faded colors;
    # grouping by colors can be canonical labels (default) or equivalence labels 
    # (colorGrouping="equivalence")
    ################
    # decide color grouping...
    if(colorGrouping=="equivalence") {
      # ... according equivalence labels
      eclassLabels = as.character(o1.plot$Initial)  # make vector as long as necessary
      for(spkr in unique(o1.plot$Vpn)) {
        temp = o1.plot$Vpn == spkr
        eclassLabels[temp] = equal_class(as.character(o1.plot$Initial[temp]), as.character(o1.plot$V[temp]))
      }
      colorGroupVector = factor(eclassLabels,levels=phonClusterLog.names)
    } else {
      # ... canonical labels
      colorGroupVectorOrig = factor(o.plot$Initial,levels=phonClusterLog.names) 
      colorGroupVector = factor(o1.plot$Initial,levels=phonClusterLog.names) 
    }
    mainString = paste("interactions: ",as.character(sG*simGroupSize),", speaker: ",plotSpeaker,", colors: ",colorGrouping," labels")
    if(plotCloud) {
      # iteration 3D plot 
      Claspect=c(1,1)
      if(dispOrig) cloudOverlayOrig = cloud(P3 ~ P1 * P2, data=o.plot,group=colorGroupVectorOrig,col=fadedColors,scales = list(arrows = FALSE),xlim=xlim,ylim=ylim,zlim=zlim,xlab=xname,ylab=yname,zlab=zname)
      cloudOverlaySim = cloud(P3 ~ P1 * P2, data=o1.plot,group=colorGroupVector,col=strongColors,scales = list(arrows = FALSE),xlim=xlim,ylim=ylim,zlim=zlim,main=paste("DCT-0...2 space, ",mainString),xlab=xname,ylab=yname,zlab=zname,key=generalLegend)
      if(dispSim) {
       if(dispOrig) {
        print(cloudOverlaySim + cloudOverlayOrig)
       } else {
        print(cloudOverlaySim)
       }
      }
      if(printPNG) {
        PNGname = paste(dirPNG,paste("Sim",as.character(sGnum),".png",sep=""),sep="/")
        png(PNGname,1024,768)
        if(dispOrig) {
          print(cloudOverlaySim + cloudOverlayOrig)
        } else {
          print(cloudOverlaySim)
        }
        dev.off()
      }
    } else if(plotTracks) {
      # aggregate: first reconstruct all tracks in o1 from DCTs, then 
      # aggregate for mean and sd (to plot confidence interval)
      o1.dctTracks.df = build_dct_df.sibilant(o1.plot)
      # rename F2Pop as M1smooth
      names(o1.dctTracks.df)[1] = "M1smooth"
      # calculate the mean at each time point
      # calculate the SD at each time point, separately per speaker and per Initial class
      if(!plotSpeakerAverage) {
        m1m = aggregate(M1smooth ~ times * Vpn * Initial + Cond, mean, data = o1.dctTracks.df)
        m1sd = aggregate(M1smooth ~ times * Vpn * Initial + Cond, sd, data = o1.dctTracks.df)
      } else {
        m1m = aggregate(M1smooth ~ times * Initial + Cond, mean, data = o1.dctTracks.df)
        m1sd = aggregate(M1smooth ~ times * Initial + Cond, sd, data = o1.dctTracks.df)
      }
      # this is for ±1.96 standard deviations, analogous to what they do in SS-ANOVA
      k = 1.96
      # the lower limit
      lower = m1m$M1smooth - k * m1sd$M1smooth
      # the upper limit
      upper = m1m$M1smooth + k * m1sd$M1smooth
      # plot orig data tracks as dashed lines
      # group = Initial works here because we have only one track per Vpn and per Initial!
      # /S/ : blue, /s/ : red, /str/ : green
      if(plotOnlySibilant=="") {
        #col=c(4,2,3)
        col=strongColors[c(1,5,7)]
      } else {
        if(plotOnlySibilant=="s") { col=strongColors[5] }
        if(plotOnlySibilant=="S") { col=strongColors[1] }
        if(plotOnlySibilant=="str") { col=strongColors[7] }
      }
      ylim=c(-2.5,2.5)
      # plot the mean with a thick dashed line
      # plot the upper and lower confidence in thin dashed lines
      if(!plotSpeakerAverage) {
        plot.end.m = xyplot(M1smooth ~ times | Vpn, group = Initial, type = c("l", "g"),
          data=m1m,ylim=ylim,col=col,lty=2,lwd=3,main=paste("M1 tracks, ",mainString))
        plot.end.u = xyplot(upper ~ times | Vpn, group = Initial, type = c("l", "g"),
          data=m1sd,ylim=ylim,col=col,lty=2,lwd=1)
        plot.end.l = xyplot(lower ~ times | Vpn, group = Initial, type = c("l", "g"),
          data=m1sd,ylim=ylim,col=col,lty=2,lwd=1)
      } else {
        plot.end.m = xyplot(M1smooth ~ times | Initial, group = Initial, type = c("l", "g"),
          data=m1m,ylim=ylim,col=col,lty=2,lwd=3,main=paste("M1 tracks, ",mainString))
        plot.end.u = xyplot(upper ~ times | Initial, group = Initial, type = c("l", "g"),
          data=m1sd,ylim=ylim,col=col,lty=2,lwd=1)
        plot.end.l = xyplot(lower ~ times | Initial, group = Initial, type = c("l", "g"),
          data=m1sd,ylim=ylim,col=col,lty=2,lwd=1)
      }
      plot.end = plot.end.m + plot.end.u + plot.end.l # plot with main title must come first!
      # plot data based on current agent's memories as dashed lines
      # group = Initial works here because we have only one track per Vpn and per Initial!
      # /S/ : blue, /s/ : red, /str/ : green
      if(dispSim) { 
       if(dispOrig) {
         print(plot.end +  plot.orig) # plot with main title must come first!
       } else {
         print(plot.end)
       }
      }
      if(printPNG) {
        PNGname = paste(dirPNG,paste("Sim",as.character(sGnum),".png",sep=""),sep="/")
        png(PNGname,1024,768)
        if(dispOrig) {
          print(plot.end + plot.orig) # plot with main title must come first!
        } else {
          print(plot.end)
        }
        dev.off()
      }
    } else {
     # default sim plot
     # iteration plot 2D, blue = /S/, red = /s/, green = /str/
     # Phonological classes are expressed by superimposed symbols 'x', '+', ...
     # plot DCT-0 vs DCT-2 only since DCT-1 is not distinguishable
     if(plotSpeakerAverage) {
       # the following additional layers can be used to show a different grouping of the 
       # data (here: the derived labels) by superimposed symbols together with the color grouping;
       # note that you must add this layer then in the later prints (see examples there)
       #if(dispOrig) xyOverlayOrig01 = xyplot(P3 ~ P1,type='p',col=fadedColors,data=o.plot,xlim=xlim,ylim=zlim,pch=as.numeric(factor(o.plot$V))+1)
       #xyOverlay01 = xyplot(P3 ~ P1,type='p',col="black",data=o1.plot,xlim=xlim,ylim=zlim,pch=as.numeric(factor(o1.plot$V))+1)

       # initial data layer with faded colors
       if(dispOrig) xyOverlayOrig02 = xyplot(P3 ~ P1,group=colorGroupVectorOrig,type='p',col=fadedColors,data=o.plot,xlim=xlim,ylim=zlim)
       # layer with actual data colored 
       xyOverlay02 = xyplot(P3 ~ P1,group=colorGroupVector,type='p',col=strongColors,data=o1.plot,xlim=xlim,ylim=zlim,main=paste("DCT-0 - DCT-2 , ",mainString),xlab=xname,ylab=zname,key=generalLegend)
     } else {
       # the following additional layer can be used to show a different grouping of the 
       # data (here: the derived labels) by superimposed symbols together with the color grouping;
       # note that you must add this layer then in the later prints (see examples there)
       #if(dispOrig) xyOverlayOrig01 = xyplot(P3 ~ P1|Vpn,type='p',col=rgb(0,0,0,50,maxColorValue=256),data=o.plot,xlim=xlim,ylim=zlim,pch=as.numeric(factor(o.plot$V))+1)
       #xyOverlay01 = xyplot(P3 ~ P1|Vpn,type='p',col="black",data=o1.plot,xlim=xlim,ylim=zlim,pch=as.numeric(factor(o2.plot$V))+1)

       # initial data layer with faded colors 
       if(dispOrig) xyOverlayOrig02 = xyplot(P3 ~ P1|Vpn,group=colorGroupVectorOrig,type='p',col=fadedColors,data=o.plot,xlim=xlim,ylim=zlim)
       # layer with actual data colored 
       xyOverlay02 = xyplot(P3 ~ P1|Vpn,group=colorGroupVector,type='p',col=strongColors,data=o1.plot,xlim=xlim,ylim=zlim,main=paste("DCT-0 - DCT-2 , ",mainString),xlab=xname,ylab=zname,key=generalLegend)
     }

     # 20161104 : print version with larger symbols and labels instead of super-imposed colors
     # (therefore only one layer!)
     #xyOverlay01 = xyplot(P3 ~ P1,type='p',col="black",data=o1.plot,ylab=list("DCT-2",cex=2),xlab=list("DCT-0",cex=2),main=list(paste("M1 tracks : ",as.character(sG*simGroupSize)," interactions"),cex=2.5),scales=list(x=list(cex=2.5),y=list(cex=2.5)),xlim=xlim,ylim=zlim,cex=2,pch=as.numeric(factor(o1.plot$V))+1,panel=function(x,y,...){
	#panel.xyplot(x,y,...);
      #ltext(x=x,y=y,labels=o1.plot$Initial,pos=1,cex=1.8)
      #})

     if(dispSim) {
      if(dispOrig) {
       print(xyOverlay02 + xyOverlayOrig02)
       #print(xyOverlay02 + xyOverlayOrig02 + xyOverlay01)
      } else {
       print(xyOverlay02)
       #print(xyOverlay02 + xyOverlay01)
      }
     }
     if(printPNG) {
       PNGname = paste(dirPNG,paste("Sim",as.character(sGnum),".png",sep=""),sep="/")
       png(PNGname,1024,768)
       if(dispOrig) {
         print(xyOverlay02 + xyOverlayOrig02)
         #print(xyOverlay02 + xyOverlayOrig02 + xyOverlay01)
       } else {
         print(xyOverlay02)
         #print(xyOverlay02 + xyOverlay01)
       }
       dev.off()
     }
    }
    ###################
    # end simulation plots
    ###################
     
    ###################
    # Make some helpful progress print-outs to console
    # This can be anything you thing of being interesting; in this example we
    # print out the current phonological clusters of each agent to see how these
    # are changing from one simGroup to the next
    ###################
    # log equivalence phone clusters of each agent into phonClusterTable_<agentName>.df
    cgroups = get.cgroups(o1)
    if(abmRun!="") {
      ABM = rep(as.character(abmRun),n=nrow(cgroups))
    } else { 
      ABM = rep("ALL",times=nrow(cgroups))
    }
    SimGroup = rep(sG,times=nrow(cgroups))
    CGroups.df = data.frame(ABM=ABM,SimGroup=SimGroup,cgroups)
    write.table(CGroups.df,file=paste(LogDirDate,"/phonClusterTable_",abmRun,".df",sep=""),append=T,col.names=F,row.names=F)
    # make statistics output of categorized phone clusters
    if(printMessages) cat("ABM: ",abmRun,", Sim Group ",sG," finished (",sG*simGroupSize,"interactions)\n")
    cLine = calc.clusterStat(o1)
    if(printMessages) print(cLine)
    write.table(cLine,file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T,col.names=F,row.names=F)
    # add cluster occurance counters to log file for summarized plot later
    phonClusterLog = rbind(phonClusterLog,cLine)

    percentRejections = 100 * get.ratioOfRejections(myPop)
    averageDerivedLabels = sum(table(o1$Vpn,o1$V)!=0)/length(unique(o1$Vpn))
    if(printMessages) cat("Average rejections: ",percentRejections,"% Average derived clusters per agent: ",averageDerivedLabels,"\n")
    sGnum = sGnum + 1
  }
  # end of loop - single run ABM finished
  
  ##################
  # save ABM protocol: the memory content of all agents at every simGroup
  # in file LogDirDate/Dataframe.df
  # (including the initial memory data marked with '$Cond='orig'))
  if(saveDataframe) {
    if(printMessages) cat("Store Memory Log (dataframe res2.df) in ",LogDirDate,"/Dataframe_",abmRun,".df\n",sep="")
    write.table(res2.df,file=paste(LogDirDate,"/Dataframe_",abmRun,".df",sep=""))
  }
  ##################
  # make animated sim GIF and store it in LogDirDate/Animation_abmRun.gif
  if(printPNG && printPaperFigs) {
    if(printMessages) cat("Store animation in ",LogDirDate,"/Animation_",abmRun,".gif\n",sep="") 
    system(paste("convert -delay 150 -loop 0 ",dirPNG,"/Sim*.png ",LogDirDate,"/Animation_",abmRun,".gif",sep=""))
  }

  ##################
  # Check if the ABM run was 'stable':
  # Compare variation in the sibilants of 'sh' words before and after the ABM:
  # If the ABM works in realistic boundaries, there should be no increase of variability
  # in each agent (the general form of data point cloud stays the same, regardless where
  # the cloud is located). Therefore we calculate the variability as the sum of diagonal
  # of the covariance matrix for each agent and them take the mean accross agents
  # calculate agents' gaussian based on the initial data stemming from
  # 'sh' words
  agentSigma = NULL
  for(j in unique(str.df$Speaker)) {
    temp = res2.df$Initial=="S"&res2.df$Cond=="orig"&res2.df$Vpn==j
    Params = as.matrix(res2.df[temp,grep("P.", names(res2.df))])
    S.gaussian = train(Params)
    agentSigma = c(agentSigma,sum(diag(S.gaussian$cov)))
  }
  agentSigma.orig = mean(agentSigma)
  # do the same on the 'sh' word data after the ABM (regardless of their derived labels)
  agentSigma = NULL
  for(j in unique(str.df$Speaker)) {
    temp = res2.df$Initial=="S"&res2.df$Cond==sG&res2.df$Vpn==j
    Params = as.matrix(res2.df[temp,grep("P.", names(res2.df))])
    S.gaussian = train(Params)
    agentSigma = c(agentSigma,sum(diag(S.gaussian$cov)))
  }
  agentSigma.abm = mean(agentSigma)
  if( agentSigma.abm > agentSigma.orig) {
    cat("WARNING: coreABM reports average increase of 'sh' words variability in agents\n        averageSigmaBefore: ",agentSigma.orig,", averageSigmaAfter: ",agentSigma.abm,"\n        ABM might not be stable, check feature distributions after ABM\n")
    cat("WARNING: coreABM reports average increase of 'sh' words variability in agents\n        averageSigmaBefore: ",agentSigma.orig,", averageSigmaAfter: ",agentSigma.abm,"\n        ABM might not be stable, check feature distributions after ABM\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  }
  if( agentSigma.abm < 0.8*agentSigma.orig) {
    cat("WARNING: coreABM reports average decrease of 'sh' words variability of more than 20% in agents\n        averageSigmaBefore: ",agentSigma.orig,", averageSigmaAfter: ",agentSigma.abm,"\n        ABM might not be stable, check feature distributions after ABM\n")
    cat("WARNING: coreABM reports average decrease of 'sh' words variability of more than 20% in agents\n        averageSigmaBefore: ",agentSigma.orig,", averageSigmaAfter: ",agentSigma.abm,"\n        ABM might not be stable, check feature distributions after ABM\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  }

