# ABM simulation script - Australian English s-retraction by Mary Stevens

# F. Schiel, J. Harrington, R. Winkelmann, M. Stevens
# version 2.0 - 2017-02-17 - changed fixed dim to variable dim (max 20!),
#               moved df production to Rcmd/emuR2ABMdataframe.R,
#               moved pre-liminary analysis examples to Rcmd/analysis.R,
#               moved before/after plots to Rdmd/beforeAfterPlots.R,
#               added option doSplitAndMergeBeforeABM and better logging
#               of cluster statistics into LogDirDate/abmParameters.txt,
#               made phonsplit() deterministic (was stochastic before),
#               added 'ABM stability test' at the end of coreABM.
# version 3.3   several improvemnets as in PublicABM.3;
#               fixed perc.sibilant.1() 
# version 3.4   added independent script Rcmd/postAnalysis (specific but 
#               could be used a sa template for other experiments:
#               post-experimental analysis of the log files of multiple ABMs
# version 3.5   bug fix: the default simulation plots calculated the equivalance
#               labels not per agent but over the accumulated memory contents, 
#               resulting in more merged equiv.labels than were actually there. 
# version 3.6   minor bug fix in simulation plot of tracks
#                


############################################################
# CHANGE THIS IF YOU MOVE THE EXPERIMENT TO ANOTHER COMPUTER
# to the installation dir of the package (where you see 'Rcmd') 
pfad2 = "/homes/schiel/ABM/PublicABM.3"
############################################################


# This script exemplifies a Agent Based Modelling after the fashion 
# of Harrington & Schiel 2017, Language
# applied to within-group s-retraction in Australian English and extended
# by the possibility of split&merge of agent-specific phonological categories
# (called 'derived labels' from now on).
# It can be used to adapt the basic ABM technique to other 
# use cases.

# Refer to the following paper for background and details: 
# Harrington, J. & Schiel, F. (in press) /u/-fronting and agent-based modeling: 
# The relationship between the origin and spread of sound change. Language, in press.
# download: http://www.phonetik.uni-muenchen.de/~jmh/papers/harringtonschiel.pdf

# This ABM assumes that the exchanged features between agents are triplets
# of three DCT-0 ... DCT-2 derived from a feature track within a phoneme 
# (e.g. in this example the first spectral moment over the time of a 
# sibilant /s/ or /ʃ/). In the remainder script we use /S/ instead of /ʃ/. You can 
# change the script to another number of dimensions or parameterisation by 
# changing the ... in ...
# But this will make all the track recovery plots where the feature track 
# is re-covered from the DCT triplets stored in memory brake! You'll have either to 
# adapt these plots to the new dimension of DCT coefficients or remove them.
 
# This ABM models mainly within-group interactions; there is no group contact
# in the main hypothesis, but rather a grouping according word types 
# (word initial 's', 'str' and 'S'). 
# However, to give examples for group-related sound change
# modelling we add a group category 'gender' to the main dataframe str.df 
# and show some example plots accross genders and give some comments at 
# possible locations in the script where agent grouping might be applicable.
# In general you can think of the 'sibililant class' or 'word initial'
# expressed in the column str.df$Initial as the grouping factor for this ABM;
# If you intend to model 'real' agent groups you'll have to replace this factor
# by a grouping factor such as 'gender'. 

# In this example script two basic types of ABM experiments are implemented
# which the user can select interactively (or by setting if running as a batch):
# 1. Singe-run ABM runs one simulation of simGroups * simGroupSize interactions,
# every simGroup plots are created and the agent memories are stored; single-runs
# can be performed with ploting group-level distributions/tracks or for a single 
# agent only.
# 2. Multiple-run ABM runs the same ABM multipleABMRuns times and collects plots and 
# statistics in the LogDirDate directory; graphic outputs during the ABM are suppressed,
# but can be stored as annimations (ImageMagick must be installed on your system!);
# this is usually run in the background and over night ;-)

# The term 'Harrington Rule' in the following refers to the basic idea that 
# sound change from one categorical feature distribution in the direction of 
# another category is driven by a simple rule, that only word tokens 
# are memorized by an agent, if they are close to the class distribution 
# within the memory of the agent. In this script the first (extended) method
# based on probabilities and the second method based on MH distances as 
# described in Harrington & Schiel 2017, Language, Section 4, are implemented,
# i.e. memorization is retricted either by:
# 1. the received token must be 'recognized' by the correct phonological class
# (the posteriori probability must be higher than other, concurring phonological
# classes), this works fine if there are concurrying classes (such as s and sh), or
# 2. the log Mahalanobis distance of the received token to the phonological distribution
# must be less or equal to a threshold to prevent outliers to memorized, this is often necessary
# if the phonological classes inthe experiment are acoustically to far apart and do not really
# compete in perception (such as /u/ and /i/).
# The latter constraint is independent of the other phonological classes,
# and just dependent on the closeness *and form* of the Gaussian of the exchanged class.

# Furthermore, memory loss (= the deletion of an exemplar from agent memory) is
# controlled either by age, or by outlier-removal (the largest Mahalanobis distance). 

# The term 'split&merge' in the following refers to a basic data-driven
# feature clustering based on k-means developped by J. Harrington in 2016.

# The software structure is as follows:
# This is the master script residing in sub-dir Rcmd; running the script 
# by issueing  'source("<path>/master.R")' on the R commandline
# should lead the experimenter through the simulation. At numerous times
# the user is asked to confirm plots or printouts before continueing and/or
# is asked to select different ways to simulate (basically a single run or
# multiple runs). In the same location as the master script there must be 
# the script pathsAndLibraries.R which defines file paths (e.g. the 
# installation dir on the local computer), libraries/functions to be loaded,
# etc. Be sure to use an up-to-date R version, since you will probably have 
# to install some packages before the master.R script runs.

# In the same location as the master.R script there are several other scripts
# that are sourced in this example, but are basically independent parts of the 
# experiment. In this example these are:
# analysis.R  : plots of the synchroneous input data before ABM
# coreABM.R   : a single ABM run
# beforeAfterABMPlots.R : plots to compare memory contents before/after a single ABM
# There are of course dependencies between scripts: for instance, it does not make 
# any sense to run beforeAfterABMPlots.R without running coreABM at least one time.
# Also scripts exchange a lot of global variables that control the behavior, 
# most of them are set at the beginning of this master script and stay constant
# throughout the experiment. 
# An independent script is postAnalysis.R that can be used to analyse the results of 
# multiple ABM runs; the time of 'after' analysis can be changed, e.g. if the original 
# multipleABM run had 100 simGroups, but you discover that after 40 simGroups nothing 
# changes, you can perform the postAnalysis at 40 simGroups instead of 100.

# For each run a log dir <date&time> is created in the 
# sub dir defined by LogDir. All results and protocols e.g. animation of the 
# feature space, log files, feature distributions/tracks before/after the ABM etc.
# are store here. Most of the names are selfexplanatory. Beware that these lof dirs
# can be get very large especially for multiple ABM runs.

# The sub dir 'data' stores initial feature data to be loaded by the ABM
# (in most cases dataframes and/or trackobjects) 

# The subdir 'functions' stores low level functions partly developed by 
# R. Winkelmann, J. Harrington and F. Schiel. Some of these exist in several different 
# forms; the master
# script can select which functions are being used via parameters passed to
# the ABM simulation (see parameters 'percFuncName' 'prodFuncName' below)
# In Rcmd/pathsAndLibraries.R these functions are loaded.



# Input data

# One of the main innovative features of Harrington ABM is the usage of real
# synchroneous production data observed from speech recordings. Therefore, the 
# first part of this script is dedicated to create these input data in form of a dataframe
# that can then be used to initialize agent memories. Basically, this means that 
# we make a very simple assumption, namely that the exemplars in the memory of an 
# agent are the same as his observed productions at a certain point in time.
# The input dataframe is called str.df in the following and must be conform with 
# the code (here and in many functions residing in the sub dir functions). Here
# is a brief description that might help to adapt the script to other types of
# ABM experiments:
# Each line of str.df represents one exemplar token observed in the acoustics of 
# a real speaker. The obligatory columns in str.df and their meaning for the ABM 
# are (all other columns are ignored):
# labels    : the phonological label of the spoken phoneme (here: 's' or 'S');
#             this implicitely defines the initial phonological inventory of each speaker
#             (here: every agent start with 2 phoneme classes 's' and 'S'; and to make 
#             things even more complicated, 'labels' is called 'V' in the later part 
#             of the script, sorry)
# Word      : the word class of the word in which the phoneme is contains
#             (here: 'sane', 'stream', 'shame' etc.)
# Speaker   : the speaker ID; if the grouping according gender should work, the 
#             ID must be of the form 'F##' (female) or 'M##' (male)
# Initial   : the phonological labels (here: 's', 'str' or 'S'
#             this never changes because it is defined by the word label, insofar 
#             this is a redundant information but handy to have.
# Gender    : the speaker grouping factor; in other types of ABM experiments (e.g. group 
#             contacts) this can be used to mark a speaker group instead (or something else)
# P1,P2,P3  : the actual DCT coefficients of the track data: DCT-0 ... 2, also refered
#             to as 'Height', 'Slope' and 'Curvature'; this is the 'signal' exchanged
#             the agents; you can chose the number of features (= dimensionality of the ABM 
#             feature space) just by adding more columns named 'P...'.

# Plots

# Pictures are created as PDF/PNG in the LogDir.
# Several switches control the behaviour of the script, i.e. if plots are 
# generated during the simulations and whether and which time-slice frames 
# (*.pgn) are generated in the PNGs subdir that later 
# can be aggregated into animated GIFs (see boolean definitions below). For this
# the software ImageMagick (command line command 'convert') must be installed.
# Ubuntu: sudo apt-get install imagemagick
# Mac: http://imagemagick.org/script/binary-releases.php#macosx


############################################################
# Begin main script
# call with 'source("<pfad2>/Rcmd/master.R")
############################################################
 
############################################################
# read packages and define file locations
source(paste(pfad2,"/Rcmd/pathsAndLibraries.R",sep=""))
############################################################

# Parameters

LogDir = paste(pfad2,"LogDir",sep="/")
                     # Sub dir where the logs are stored from each ABM run
printMessages = T    # Print messages on console during coreABm run (the parameter list is always
                     # printed).
runInteractive = T   # It set to F, the script runs without asking any questions and without plotting graphics 
                     # to the X device; good for multiple runs in the background
runMode = "multiple" # select run mode of ABM when runInteractive=F; can be superceeded in interative mode
                     # 'single' : one ABM is performed and before analysis (Rcmd/analyse.R) and after analysis
                     # (Rcmd/beforeAfterABMPlots.R) are carried out; this mode is useful for trying out parameters
                     # and watching the ABM while running (runInteractive=T,dispSim=T...)
                     # 'multiple' : multipleABMRuns are performed and some basic statistic plots generated;
                     # before analysis (Rcmd/analyse.R) is done automatically, if runInteractive=F; after analysis
                     # (Rcmd/beforeAfterABMPlots.R) is not carried out; instead source the script Rcmd/postAnalysis.R after 
                     # the multiple ABM run.
multipleABMRuns = 100  # number of multiple ABM runs to perform in multiple mode

# The following booleans define mainly the plots/fig creations during the simulations.

printPaperFigs = T   # create all plotted analysis figures as PDF in LogDirDate
                     # including Animation derived from PNGs
dispSim = T          # display plots during ABM simulation on X device; default (both
                     # plotCloud and plotTracks are F) is plotting DCT-0...DCT-2 space
                     # (will be switched off in multiple run mode)
#colorGrouping = "canonical"
colorGrouping = "equivalence"
                     # label coloring in animation plots:
                     # 'canonical' (default) : colors in 2D and 3D plots encode canonical
                     #                         labels 's' 'S' 'str'
                     # 'equivalence' : colors in 2D and 3D plots encode equivalence labels (7)
plotCloud = F        # display 3D plots of DCT space rather than 2-D plots of DCTs
plotTracks = F       # display re-constructed tracks rather than 2-D plots of DCTs
                     # only one of plotCloud or plotTracks should be set to TRUE
dispOrig = F         # display initial values in weak colors as overlay in all plots
                     # (this can make plots a bit fuzzy!)
plotSpeakerAverage = F # where applicable, plot the average across speakers rather than 
                       # a panel plot of speakers separately
plotSpeaker = ""     # constrain plots to this speaker; "" = all speakers
plotOnlySibilant = ""
                     # constrain all plots to a single phonological class e.g. "str". If empty string,
                     # all sibilants are plotted
printPNG = T         # create 'Sim1###.png' files in dirPNG containing the steps of simulation
                     # where ### is 0000 ... simGroups-1 (e.g. Sim1000 - Sim1019)
                     
# Plot parameters

xlim = c(-2.5,2)     # plot dimensions in the DCT0...2 space for all plots
ylim = c(-2,1)
zlim = c(-1.5,0.5)
xname = "DCT-0"      # labelling of the three dimensions of the exchanged data triplets
yname = "DCT-1"
zname = "DCT-2"
MainLab= "M1 DCT space"

# the following are name definitions of equivalence labels and their plot colors (very specific for this example!)
phonClusterLog.names = c("S","S+s","S+s+str","S+str","s","s+str","str")
                    # this vector is used to sort agent-specific equivalence phone cluster statistics
                    # (must be char sorted!), their colors in all plots in the same order
                    # blue,pink,darkgrey,cyan,red,orange,green
                    # useful: col2rgb("orange")
strongColors = c(rgb(0,0,255,250,maxColorValue=255),rgb(255,0,255,250,maxColorValue=255),rgb(60,60,60,250,maxColorValue=255),rgb(0,255,255,250,maxColorValue=255),rgb(255,0,0,250,maxColorValue=255),rgb(255,165,0,250,maxColorValue=255),rgb(0,255,0,250,maxColorValue=255))
fadedStrength = 75  # faded colors: fading factor 1...255 (not faded)
fadedColors = c(rgb(0,0,255,fadedStrength,maxColorValue=255),rgb(255,0,255,fadedStrength,maxColorValue=255),rgb(60,60,60,fadedStrength,maxColorValue=255),rgb(0,255,255,fadedStrength,maxColorValue=255),rgb(255,0,0,fadedStrength,maxColorValue=255),rgb(255,255,0,fadedStrength,maxColorValue=255),rgb(0,255,0,fadedStrength,maxColorValue=255))
generalLegend = list(points=list(pch = 1, col = strongColors), text = list(phonClusterLog.names, col = strongColors),columns = 4, col=strongColors)


# General parameters of ABM (probably used in every ABM)
# (typical values for Mary's s-retraction experiment that leads to a stable
#  situation with phonological re-structuring in agents are in given in [...]) 

# Data source 
loadEmuData = F      # set to T, if you want to process the original data anew and overwrite 
                     # the pre-processed dataframe in sub dir 'data'; otherwise (default) the 
                     # the pre-processed dataframe is loaded.
                     # make production data from scratch (T) requires an intact emuDB.
saveDataframe = T    # save ABM protocol in form of a 'res.df' dataframe to file <date&time>.df in LogDir
                     # (single ABM run), or for a multiple-ABM-run as '###Dataframe.df' in LogDir
                     # where '###' is an agent



speakerNormalize = T # speaker normalize production data before analysis and ABM using Lobanov 
                     # this means that the entire ABM runs in a 'normalized feature space', that is 
                     # we assume that each agent has a perfect way to compensate for speaker 
                     # individual shifts and SD differences; only relative changes within 
                     # an agents data set are modelled. This also changes all dimensions of 
                     # plots to a dimensionless number [T]
simMahalThreshold = 1.5
                   # the log Mahalanobis threshold for acception an incoming token in ABM: 
                   # s-retraction (2017 data) : 1.5-2.0
                   # u-fronting: 2.5 [2.7]
                   # values below 1.5 lead to very rare memorization in agents, since the log Mahal.distance 
                   # is usually around 1.5 - 3.5, and therefore the distributions collapse 
                   # after 20.000 interactions.
simGroups = 20      # number of simulation groups [60]
simGroupSize = 1000 # number of iterations per simGroup [1000]
maxTokens = 410     # max number of word tokens in agent's memory; if exceeded, the agent starts
                    # to 'forget' the oldest tokens of a word class when memorizing a new token
                    # if the chosen percFuncName allows this [140]
splitAndMerge = T   # if set, the perc.sibiblant.2 function will check the agent's memory for 
                    # possible splits or mergers every time a split&merge is triggered (see below). [T]
splitAndMergeInterval = 100
                    # the interval, counted in personal updates of memory, after which an agent 
                    # performs a split&merge check on his memory. If set to 1, every time
                    # a token is received, the check is performed; if set to value > 1, 
                    # the number of received tokens %% splitAndMergeInterval is checked for 0
                    # to trigger split&merge. [20...100]
doSplitAndMergeBeforeABM = T
                    # before the first interaction, do a split&merge on all agents; that way the initial 
                    # plots do not show the /s_str/+/S/ cluster situation for all agents

# Parameters specific to the ABM case (these most likely change in a new experiment)

                    # optional select subset of speakers (see also comment 'filter to agent subset' below)
subsetSpeakers = NULL   # default: all agents found in str.df

                    # speaker/listener selection probabilities: if these vectors exist, the ABM choses speakers 
                    # according to the weight proportions in these vectors, e.g. agent myPop[[1]] is selected
                    # with prob speakerProb[1]/sum(speakerProb); zeros are allowed (= never speaking/listening!)
                    # these two vectors must have the same length as number of agents and do not contain negatives;
                    # if the vector does not exist, equal probabilities are used (backwards compatible) 
# example:
#noOfSpeakers = length(subsetSpeakers)
#noOfNonnativeSpeakers = length(subsetSpeakers)-2
#noOfNativeSpeakers = 2
#speakerProb = rep(1/(2*noOfNonnativeSpeakers),times=noOfNonnativeSpeakers)
#speakerProb = c(speakerProb,0.25,0.25)   # native speakers S100 S101 (teachers) have 50% speaking prob
#listenerProb = rep(1/noOfNonnativeSpeakers,times=noOfNonnativeSpeakers)
#listenerProb = c(listenerProb,0,0)       # native speaker S100 S101 (teachers) never listen


percFuncName = "perc.sibilant.1"
                    # the name of the basic function that handles the perception process of a single token
prodFuncName = "prod.sibilant.3"
                    # the name of the basic function that handles the production process of a single token
                    # These functions must be in some sourced module, most likely in sub dir 'functions'
                    # and they must match the memory structure of the agents in this ABM experiment.
                    # Some implemented functions that work in this script:
                    # [perc.sibilant.3]
                    # prod.sibilant.3 : The same as ...sibilant.2 but extended to variable dim feature sets
                    # [perc.sibilant.2] 
                    # prod.sibilant.2 : The strategy used in M. Steven's s-retraction experiments on Australian English
                    #                   The main difference to the standard ABM is the treatment of /s/ in the context
                    #                   of 'str...' words: agents store never reject them,
                    #                   i.e. they are *always* memorized, and the Harrington Rule (in form of a MH threshold) 
                    #                   is only applied to the 'other' /s/ and the /S/. With other 
                    #                   words: agents try to keep nice sharp and stable phonological models for /s/ and /S/ 
                    #                   perception (and reject tokens that do not fit into their model), but they don't care 
                    #                   much about the /str/ exemplars. This trick together with split&merge prevents that /s/ in 
                    #                   'str...' context and 'other' /s/ collapse (as they would in a standard model), but stay 
                    #                   rather stable, often split from the /s/ class and sometimes even merge with the /S/ class.
                    #                   These functions use a fixed set of three dynamic DCT features; see 
                    #                   ...sibilant.3 functions for a version that can use arbitrary feature dimensions
                    # perc.sibilant.1 : the original strategy of Harrington in the very first ABM on u-fronting in 
                    #                   younger/older speaker of English extended by: perceived tokens are only 
                    #                   accepted, if their posterior phone probability is the highest of all classes regarding the
                    #                   quadratic discriminant phone model of the receiving agent (maxPosteriorMatch=T) or is 
                    #                   higher that 1/nr_of_classes (maxPosteriorMatch=F, not recommended); removed token is the token that 
                    #                   is farest away from phone probability Gaussian (timeDelayStrategy=F) or the oldest token 
                    #                   (timeDelayStrategy=T); split&merge is optional (splitAndMerge=T).
                    #                   perc.sibilant.1 can be used together with prod.sibilant.2|3 (the production strategy
                    #                   is the same!); recommended: timeDelayStrategy=T because otherwise
                    #                   distributions collapse rather quickly; splitAndMerge=T and maxPosteriorMatch=T.
                    # IMPORTANT NOTE: perc.sibilant.1 works only inan ABM with more than one phonological class! Use
                    #                   prod.sibilant.3 otherwise!

timeDelayStrategy = T
                    # if set to FALSE, the script switches to the original token removal strategy of Harrington's
                    # first experiment (= remove the outlier token instead of the oldest token); 
                    # currently this parameter is only effective for percFuncName = "perc.sibilant.1"; 
                    # if you use perc.sibilant.2, timeDelayStrategy is always true!
maxPosteriorMatch = T
                    # Variant of perc.sibilant.1(): if maxPosteriorMatch=T the original 'soft' decision rule (u-fronting)
                    # is replaced by a rule that requires the incoming token to have the maximum posterior of all classes. 
                    # Multiple ABM runs show that this outperforms the older 'soft' rule; therefore this should always be
                    # set to T, if the perc.sibilant.1() function is used!
compressedSD = 1.0  # factor to compress production SD in str words in ABM 
                    # currently this parameter is only effective for percFuncName = "perc.sibilant.2", not recommended!
threshold.str = simMahalThreshold  # = same Harrington Rule as for the other phones
#threshold.str = 99.0  # = no rule: all /str/ tokens are memorized
                    # log MH threshold for sibilant perception of 'str' words in ABM
                    # extrem high number means no Harrington Rule is applied for <str> types 
                    # (all tokens accepted). [simMahalThreshold]
                    # currently this parameter is only effective for percFuncName = "perc.sibilant.2"; we found that this is 
                    # often leads to an 'explosion' of the 'str' token distributions, not recommended!
#articulStretch = 0.5  # factor that stretches the articulatory space for production from a cube defined 
                    # by the 1% - 99% quantiles of start data; productions outside this area are suppressed;
                    # a very high value basically disables the restriction to articulatory space
articulStretch = 99.9 
                    # no restriction (we found that this is not really necessary, and 
                    # makes things harder to explain, only do this in cases where the space explodes)
                    # currently this parameter is only effective for prodFuncName = "prod.sibilant.2" [99.9]
phoneticCoarticulate1 = F
                    # if set to TRUE the production function prod.sibilant.2 will only produce 
                    # /str/ tokens that are on the side of the /str/ Gaussian where the /S/ 
                    # centroid is (this of course causes a down-shift of /str/ M1 tracks, but this is not 
                    # what we want to test here; it would be a 'built-in sound change', not recommended!)
                    # currently this parameter is only effective for prodFuncName = "prod.sibilant.2" [F]
phoneticCoarticulate2 = F
                    # if set to TRUE the production function prod.sibilant.2 will shift the 
                    # /str/ Gaussian by shiftStrGaussionFactor*100 % in the direction of the 
                    # /S/ centroid (this of course causes a down-shift of /str/ M1 tracks, but this is not
                    # what we want to test here; it would be a 'built-in sound change', not recommended!)
                    # currently this parameter is only effective for prodFuncName = "prod.sibilant.2" [F]
shiftStrGaussionFactor = 0.1


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

# date and time
DateTime = format(Sys.time(), "%Y%m%d%H%M%S")

# print out ABM parameters to check

cat("Script master.R invoked with the following parameters:\n\n")
cat("DATA parameters:\n")
cat("root of experiment (pfad2): ",pfad2,"\n")
# create LogDirDate and dirPNG
LogDirDate = paste(LogDir,DateTime,sep="/")
cat("LogDir: ",LogDirDate,"\n")
if(dir.exists(paths=LogDirDate)) {
    stop("Result log dir already exists; rename or delete existing dir and run this script again!")
  }
system(paste("mkdir ",LogDirDate,sep=""))
dirPNG = paste(LogDirDate,"PNGs",sep="/")    
                     # directory to create PNG time-series files of simulation 
                     # this numbered time-slice series of *.png is then converted 
                     # into an animated GIF using e.g. 'convert' and stored in LogDir;
                     # note that the files in dirPNG are over-written by the next ABM run!
system(paste("mkdir ",dirPNG,sep=""))
cat("loadEmuData: ",loadEmuData,"\n")
cat("speakerNormalize: ",speakerNormalize,"\n")
cat("PLOT parameters:\n")
cat("printPaperFigs: ",printPaperFigs,"\n")
cat("dispSim: ",dispSim,"\n")
cat("dispOrig: ",dispOrig,"\n")
cat("plotCloud: ",plotCloud,"\n")
cat("plotTracks: ",plotTracks,"\n")
cat("plotSpeakerAverage: ",plotSpeakerAverage,"\n")
cat("printPNG: ",printPNG,"\n")
cat("dirPNG: ",dirPNG,"\n")
cat("plotSpeaker (EMPTY = ALL): ",plotSpeaker,"\n")
cat("plotOnlySibilant (EMPTY = ALL): ",plotOnlySibilant,"\n")
cat("ABM parameters:\n")
cat("prodFuncName: ",prodFuncName,"\n")
cat("percFuncName: ",percFuncName,"\n")
cat("timeDelayStrategy: ",timeDelayStrategy,"\n")
cat("maxPosteriorMatch: ",maxPosteriorMatch,"\n")
cat("simMahalThreshold: ",simMahalThreshold,"\n")
cat("threshold.str: ",threshold.str,"\n")
cat("phoneticCoarticulate1: ",phoneticCoarticulate1,"\n")
cat("phoneticCoarticulate2: ",phoneticCoarticulate2,"\n")
cat("splitAndMerge: ",splitAndMerge,"\n")
cat("splitAndMergeInterval: ",splitAndMergeInterval,"\n")
cat("doSplitAndMergeBeforeABM: ",doSplitAndMergeBeforeABM,"\n")
cat("simGroupSize: ",simGroupSize,"\n")
cat("articulStretch: ",articulStretch,"\n")
cat("maxTokens: ",maxTokens,"\n")
cat("simGroups: ",simGroups,"\n")
cat("multipleABMRuns: ",multipleABMRuns,"\n")
if(exists("speakerProb")) cat("speakerProb:\n",speakerProb,"\n")
if(exists("listenerProb")) cat("listenerProb:\n",listenerProb,"\n")
cat("saveDataframe: ",saveDataframe,"\n\n")

# store ABM parameters in LogDirDate / abmParameters.txt

  cat("DATA parameters:\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"))
  cat("root of experiment (pfad2): ",pfad2,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("LogDir: ",LogDirDate,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("loadEmuData: ",loadEmuData,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("speakerNormalize: ",speakerNormalize,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("PLOT parameters:\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("printPaperFigs: ",printPaperFigs,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("dispSim: ",dispSim,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("dispOrig: ",dispOrig,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("plotCloud: ",plotCloud,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("plotTracks: ",plotTracks,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("plotSpeakerAverage: ",plotSpeakerAverage,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("printPNG: ",printPNG,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("dirPNG: ",dirPNG,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("plotSpeaker (EMPTY = ALL): ",plotSpeaker,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("plotOnlySibilant (EMPTY = ALL): ",plotOnlySibilant,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("ABM parameters:\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("prodFuncName: ",prodFuncName,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("percFuncName: ",percFuncName,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("timeDelayStrategy: ",timeDelayStrategy,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("maxPosteriorMatch: ",maxPosteriorMatch,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("simMahalThreshold: ",simMahalThreshold,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("threshold.str: ",threshold.str,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("phoneticCoarticulate1: ",phoneticCoarticulate1,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("phoneticCoarticulate2: ",phoneticCoarticulate2,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("splitAndMerge: ",splitAndMerge,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("splitAndMergeInterval: ",splitAndMergeInterval,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("doSplitAndMergeBeforeABM: ",doSplitAndMergeBeforeABM,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("simGroupSize: ",simGroupSize,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("articulStretch: ",articulStretch,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("maxTokens: ",maxTokens,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("simGroups: ",simGroups,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("multipleABMRuns: ",multipleABMRuns,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  if(exists("speakerProb")) cat("speakerProb:\n",speakerProb,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  if(exists("listenerProb")) cat("listenerProb:\n",listenerProb,"\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
  cat("saveDataframe: ",saveDataframe,"\n\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)


# Prepare production data to initialze ABM
if(loadEmuData) {
  
  #############
  # this script basically prepares the objects data/str.df 
  # Usually we just load these below; this part is only necessary after changes 
  # the emuDBs 
  source(paste(pfad2,"/Rcmd/emuR2ABMdataframe.R",sep=""))
  
} else {
  
  #############
  # Here we re-load previously calculated data sets that represent 
  # the initial production data of our agents; in this example
  # we need a dataframe with metadata of the observed tokens and 
  # the DCT coefficients of the feature track per token
  # In the example here str.df consists of a emuR seglist containing 
  # the sibilants extended by 3 columns 'Word' (word label), 'Speaker'
  # (the anonymized speaker label) and 'Initial' (the grouping factor, 
  # here the sibilant type, i.e. 's', 'str' or 'S'.
  #############
  str.df=read.table(file.path(pfad2,"data/str.df"),stringsAsFactors=F)
  
} # end loadEmuData==F

# change 'sh' labels from ʃ to S,because we have troubles with factor ordering
In = as.character(str.df$Initial)
In[In=="ʃ"] = "S"
str.df$Initial=factor(In)
In = as.character(str.df$labels)
In[In=="ʃ"] = "S"
str.df$labels=factor(In)
str.df$Speaker = factor(str.df$Speaker)

#############
# Check some parameters
nrOfAgents = length(unique(str.df$Speaker))
if(length(str.df$Word)/nrOfAgents > maxTokens) {
  cat("WARNING: the average number of words in an agent's memory is larger than the\n          parameter maxTokens (",maxTokens,"); maxTokens has probably no effect\n") 
  cat("WARNING: the average number of words in an agent's memory is larger than the\n          parameter maxTokens (",maxTokens,"); maxTokens has probably no effect\n",file=paste(LogDirDate,"abmParameters.txt",sep="/"),append=T)
} 
Speakers = unique(as.character(str.df$Speaker))
if(plotSpeaker!="" & sum(Speakers==plotSpeaker)!=1) {
  cat("ERROR: plotSpeaker = ",plotSpeaker," is not a speaker - exiting\n")
  stop()
}
if(plotTracks & plotCloud) {
  cat("ERROR: you cannot select both, plotCloud and plotTracks - exiting\n")
  stop()
}
#############
# filter to subset
if(!is.null(subsetSpeakers)) {
  temp = str.df$Speaker %in% subsetSpeakers
  str.df = str.df[temp,]
}

################
# Preliminary data plots
# In this part of the script we look at the initial production data 
# *before* we feed them into the agent's memories. This is important
# to check for possible asymmetries in distribution accross the grouping 
# factors; only if we see assymmetries a direct effect of the Harrington 
# Rule is to be expected.
# We do this in two blocks: looking at DCT space (without saving to PDF)
# and looking at (smoothed) M1 tracks (with the possibility to save these
# plots in PDFs in the sub dir LogDirDate, if parameter printPaperFigs=T). 
#
# All the following plots can be seen as examples on how to look at the 
# synchronic data before ABM; for your special case of ABM you probably 
# have to select some of these, delete others that don't make sense and 
# possibly add new plot routines.
# If your ABM involves group contact, you'll probably look at the plots 
# were data are separated by gender or competence.
# The plots will not run with other dimensions than 3 (p1,P2,P3)
################

if(runInteractive) {
  inpu = "N"
  inpu = readline("Make preliminary plots on data before ABM? (y/N) ")
  if(inpu == "Y" || inpu == "y") {
    source(paste(pfad2,"/Rcmd/analysis.R",sep="")) 
  }
} else {
  # if running in batch mode, always make these plots
  source(paste(pfad2,"/Rcmd/analysis.R",sep=""))
}


##################
# Start of ABM
# The next section contains the two types of ABM simulations:
# a single run or a multiple run (one run for each agent).
# In this example script only one type of ABM is performed:
# Each agent is initialized with the synchronic production data
# as prepared and analysed above. Then each agent gets two phonological 
# class labels assigned to his stored word exemplars: either /s/ or /S/.
# Agents start exchanging tokens randomly and memorizing tokens according 
# the Harrington Rule. After a total of splitAndMergeInterval memorized tokens, the 
# agent checks his complete memory for possible split&merge of derived labels,
# if parameter splitAndMerge=T.
# After a totel of maxTokens memorized tokens, the agent starts to delete the oldest 
# token of a word class after each new memorization into this class.
#
# Example: splitAndMerge=T, splitAndMergeInterval=100, maxTokens=180
# Each agent starts with 10 tokens of 5 word classes (total: 50); after 100 interactions
# that lead to a memorization of the received token (total:150), the agent tests 
# his memory for split&merge for the first time; after another 30 memorized tokens 
# (total: 180), the agent starts deleteting the oldest word class token, when memorizing 
# a new token of that word class. 
##################

cat("Run ABM simulation starting with two acoustic models /s/, /S/\n")

if(runInteractive) {
  inpu = readline("Run single/multiple ABM? (S/m) ")
  if(inpu != "" && inpu != "S" && inpu != "m" && inpu != "s" && inpu != "M") {
    cat("Stopped\n")
    stop()
  }
} else {
  if(runMode=="single") inpu = "S"
  else inpu = "M"
}
##################

# initialize dataframe to log phonological cluster statistics during the ABM;
# see comments at the end of the ABM loop for details 
phonClusterLog = NULL

if(inpu == "S" || inpu == "s" || inpu == "" ) {

  ######################
  # Run the ABM core loop
  # this script will print some basic progress infos on console, and depending on 
  # parameter settings will plot progress plots on screen/file, save the accumulated 
  # simulation data (res2.df) to file (Dataframe.df), and create an annimated GIF
  # (Annimation.gif) in the log dir (LogDirDate)

  abmRun = plotSpeaker # use plotSpeaker in file naming if set

  source(paste(pfad2,"Rcmd/coreABM.R",sep="/"))

  ######################
  # Do some development plots
  # e.g. how do thing evolve over the time course of the ABM; for this we exploit 
  # the Memory Log in res2.df and phonClusterLog

  ######################
  # plot of phonological cluster development over ABM
  # Note that multiple (phonetic?) clusters within a phone class are summerized, e.g.
  # 's_s_s' becomes 's' etc.
  phonClusterPlot = NULL
  for(j in 3:(length(phonClusterLog.names)+2)) {
    phonClusterPlot = rbind(phonClusterPlot,data.frame(simGroup=phonClusterLog$simGroup,numPhonCluster=phonClusterLog[,j],
        group=rep(phonClusterLog.names[j-2],n=simGroups)))                            
  }
  # the following line: colors stable for the 7 equiv.classes even if a class is not shown here
  phonClusterPlot$group = factor(phonClusterPlot$group,levels=phonClusterLog.names)
  plotPhonCluster = xyplot(numPhonCluster ~ simGroup,group=group,col=strongColors,type = c("l", "g"),ylab="Number of agents",data=phonClusterPlot,xlab="Interactions (x 1000)",key=generalLegend)
  if(printPaperFigs) {
    pdf(paste(LogDirDate,"phonClusterPlot.pdf",sep="/"))
    print(plotPhonCluster)
    dev.off()
  }
  if(runInteractive) {
    print(plotPhonCluster)
    inpu = readline("Plot number of agents with equivalence clusters over ABM  ok? Continue? (Y/n) ")
    if(inpu != "" && inpu != "Y" && inpu != "y") {
      cat("Script stopped by user\n")
      stop()
    }
  }
  ######
  # STOP
  ######
  
  ######################
  # plot of mean DCT-0 (= the height of acoustic tracks) in memories development over ABM
  temp = res2.df$Initial=="str"
  P1_over_ABM = aggregate(P1 ~ Cond, data=res2.df[temp,],mean)
  plotHeightOverABM = xyplot(P1 ~ Cond,data=P1_over_ABM,type = c("l", "g"),ylab="Height of M1 tracks (DCT-0)",xlab=paste("Simulation groups (",simGroupSize," each)"),main="Mean absolute height of M1 tracks of 'str' words over ABM")
  if(printPaperFigs) {
    pdf(paste(LogDirDate,"plotHeightOverABM.pdf",sep="/"))
    print(plotHeightOverABM)
    dev.off()
  }
  if(runInteractive) {
    print(plotHeightOverABM)
    inpu = readline("Plot average M1 track height of 'str' words over ABM  ok? Continue? (Y/n) ")
    if(inpu != "" && inpu != "Y" && inpu != "y") {
      cat("Script stopped by user\n")
      stop()
    }
  }
  ######
  # STOP
  ######






  ##################
  # do some optional before/after plots
  # optionally store them into LogDirDate as PDFs
  # these will not run with other dimensions than 3 (P1,P2,P3)
  ##################
  if(runInteractive) {
    inpu = readline("Plot before/after ABM data? Continue? (Y/n) ")
    if(inpu != "" && inpu != "Y" && inpu != "y") {
      cat("Script stopped by user\n")
      stop()
    }
  }
  ######
  # STOP
  ######

  source(paste(pfad2,"/Rcmd/beforeAfterABMPlots.R",sep=""))


  # end single ABM run type 2
  ########################
  
} else {

  ########################
  # multiple ABM 
  # instead of a single run the user can select multiple ABM runs in the following section
  # Multiple runs are intended to run without user inspection over a log time 
  # (possible over-night) and store all results for later inspection/analysis.
  # In this example we simply run one complete ABM for each agent (thus if you 
  # model 20 agents there will be 20 repetitions of the same ABM). In each run
  # we create an animation showing the data of one agent, save the Memory log,
  # collect cluster configurations for statistics, and create some before/after plots.
  ########################
  cat("Run multiple ABMs (one for each speaker), create sim pictures and collect statistics about split/mergers per speaker\n")

  # messages, no plots during coreABM run
  printMessages = T
  dispSim = F
  # get vector of agents
  agentNames = unique(str.df$Speaker)
  nrOfAgents = length(unique(str.df$Speaker))
  abmRun = 1
  # initialize dataframe to log phonological cluster statistics during the ABMs;
  # see comments at the end of the ABM loop for details 
  phonClusterLog = NULL
  # initialize datafram to log mean DCT-values during the ABMs;
  # see comments at the end of the ABM loop for details
  HeightOverABMLog = NULL

  # loop over all speakers: for each speaker run a complete ABM,
  # monitor of phonological classes after each simGroup, res2.df (= memory contents
  # of agents after each SimGroup) is stored in LogDirDate/Dataframe_<agentName>.df
  # for later analysis
  #for(plotSpeaker in c("F01","F05","M08")) {
  #for(plotSpeaker in agentNames) {
  # OR
  # simply run multipleABMRuns times
  plotSpeaker = ""  # we want to see animation over all speakers
  for(i in 1:multipleABMRuns) {

    #cat("Multiple ABM: running ",abmRun," of ",nrOfAgents,", results dir: ",LogDirDate,"\n")
    cat("Multiple ABM: running ",abmRun," of ",multipleABMRuns,", results dir: ",LogDirDate,"\n")

    source(paste(pfad2,"Rcmd/coreABM.R",sep="/"))


    ######################
    # Do some development plots for each run (in case that we later want to see individual results!)
    # e.g. how do thing evolve over the time course of the ABM; for this we exploit
    # the Memory Log in res2.df and phonClusterLog
  
    ######################
    # plot of phonological cluster development in this ABM
    # Note that multiple (phonetic?) clusters within a phone class are summerized, e.g.
    # 's_s_s' becomes 's' etc.
    phonClusterPlot = NULL
    temp = phonClusterLog$ABM==as.character(abmRun)
    for(j in 3:(length(phonClusterLog.names)+2)) {
      phonClusterPlot = rbind(phonClusterPlot,data.frame(simGroup=phonClusterLog$simGroup[temp],numPhonCluster=phonClusterLog[temp,j],
          group=rep(phonClusterLog.names[j-2],n=simGroups)))
    }
    # the following line: colors stable for the 7 equiv.classes even if a class is not shown here
    phonClusterPlot$group = factor(phonClusterPlot$group,levels=phonClusterLog.names)
    plotPhonCluster = xyplot(numPhonCluster ~ simGroup,group=group,col=strongColors,type = c("l", "g"),main="Number of agents with equiv.labels over ABM",ylab="Number of agents",data=phonClusterPlot,xlab="Interactions (x 1000)",key=generalLegend)
    if(printPaperFigs) {
      pdf(paste(LogDirDate,"/phonClusterPlot_",abmRun,".pdf",sep=""))
      print(plotPhonCluster)
      dev.off()
    }

    #####################
    # collect average DCT values in HeightOverABMLog
    #temp = res2.df$Initial=="str"
    HeightOverABMLog = rbind(HeightOverABMLog,data.frame(ABM=abmRun,aggregate(cbind(P1,P2,P3) ~ Cond * Initial, data=res2.df,mean)))
        
    abmRun = abmRun + 1

  } # end of multiple loop

  ######################
  # plot of averaged phonological cluster development over ABM
  # Note that multiple (phonetic?) clusters within a phone class are mapped to one equiv.label, e.g.
  # 's_s_s' becomes 's'

  # save and plot development of phonological cluster configurations during the ABMs
  # average cluster counts over ABM simulations
  write.table(phonClusterLog,file=paste(LogDirDate,"PhonClusterLog.txt",sep="/"))

  # ensemble plot number of agents with S+str cluster over ABMsa (one line per ABM)
  #(the following lines replace '+' by 'M' temporarily, because aggregate/xyplot can't work with '+' 
  phonClusterLog.names.M = gsub("\\+","M",phonClusterLog.names)
  names(phonClusterLog)[3:9] = phonClusterLog.names.M
  plotPhonCluster = xyplot(SMstr ~ simGroup,group=ABM,col=strongColors,type = c("l", "g"),ylim=c(-1,nrOfAgents+1),main=paste("Agents with S+str cluster over ",multipleABMRuns," ABMs",sep=""),ylab="Number of agents",data=phonClusterLog,xlab="Interactions (x 1000)")
  if(runInteractive) {
    print(plotPhonCluster)
    inpu = readline("Plot ensemble number of agents using equivalence cluster S+str over ABMs ok? Continue? (Y/n) ")
    if(inpu != "" && inpu != "Y" && inpu != "y") {
      cat("Script stopped by user\n")
      stop()
    }
  }
  if(printPaperFigs) {
    pdf(paste(LogDirDate,"phonClusterPlotMultipleS+str.pdf",sep="/"))
    print(plotPhonCluster)
    dev.off()
  }
  ######
  # STOP
  ######

  # plot average number of agents over ABMs with equiv.labels in their memories
  phonClusterAggregate = aggregate(cbind(S,SMs,SMsMstr,SMstr,s,sMstr,str) ~ simGroup,FUN=mean,data=phonClusterLog)
  names(phonClusterAggregate)[2:8] = phonClusterLog.names # restitute names!
  names(phonClusterLog)[3:9] = phonClusterLog.names # restitute names!
  # make nice plot
  phonClusterPlot = NULL
  for(j in 2:(length(phonClusterLog.names)+1)) {
    phonClusterPlot = rbind(phonClusterPlot,data.frame(simGroup=phonClusterAggregate[,1],numPhonCluster=phonClusterAggregate[,j],
      group=rep(phonClusterLog.names[j-1],n=simGroups))) 
  }
  # the following line: colors stable for the 7 equiv.classes even if a class is not shown here  
  phonClusterPlot$group = factor(phonClusterPlot$group,levels=phonClusterLog.names)
  plotPhonCluster = xyplot(numPhonCluster ~ simGroup,group=group,col=strongColors,type = c("l", "g"),ylim=c(-1,nrOfAgents+1),ylab=paste("Number of agents averaged over ",multipleABMRuns," simulations",sep=""),main=paste("Average number of agents with equiv.labels over ",multipleABMRuns," ABMs",sep=""),data=phonClusterPlot,xlab="Interactions (x 1000)",key=generalLegend)
  if(runInteractive) {
    print(plotPhonCluster)
    inpu = readline("Plot averaged number of agents using equivalence cluster over ABMs ok? Continue? (Y/n) ")
    if(inpu != "" && inpu != "Y" && inpu != "y") {
      cat("Script stopped by user\n")
      stop()
    }
  }
  if(printPaperFigs) {
    pdf(paste(LogDirDate,"phonClusterPlotMultiple.pdf",sep="/"))
    print(plotPhonCluster)
    dev.off()
  }
  ######
  # STOP
  ######

  ######################
  # save DCT log
  write.table(HeightOverABMLog,paste(LogDirDate,"DistanceValuesOverABMLog.txt",sep="/"))

  ######################
  # ensemble plot of DCT-0 (= the height of acoustic tracks) of 'str' tokens in memories development over all ABMs (colored)
  plotHeightOverABM = xyplot(P1 ~ Cond,group=ABM,data=HeightOverABMLog[HeightOverABMLog$Initial=="str",],type = c("l", "g"),ylab="Height of M1 tracks (DCT-0)",xlab=paste("Simulation groups (",simGroupSize," each)"),main=paste("Absolute height of M1 tracks of 'str' words over ",multipleABMRuns," ABMs"))
  if(printPaperFigs) {
    pdf(paste(LogDirDate,"plotHeightOverABMMultiple.pdf",sep="/"))
    print(plotHeightOverABM)
    dev.off()
  }
  if(runInteractive) {
    print(plotHeightOverABM)
    inpu = readline("Plot average M1 track height of 'str' words over ABMs  ok? Continue? (Y/n) ")
    if(inpu != "" && inpu != "Y" && inpu != "y") {
      cat("Script stopped by user\n")
      stop()
    }
  }
  ######
  # STOP
  ######

  # make before/after pie chart of equivalence phone label distribution of all ABM runs
  # pie before
  
  # get memory contents before ABM; since these are the same in all ABM runs, we
  # simply take the ones from the last ABM run stored in res2.df
  temp = res2.df$Cond=="orig"
  res2.orig = res2.df[temp,]
  eclassLabels = as.character(res2.orig$V)
  for(spkr in unique(res2.orig$Vpn)){
    temp = res2.orig$Vpn == spkr
    eclassLabels[temp] = equal_class(as.character(res2.orig$Initial[temp]), as.character(res2.orig$V[temp]))
  }
  tab.orig = table(eclassLabels)
  col.orig = strongColors[phonClusterLog.names %in% names(tab.orig)]
  # now get the memory contents of all agents after the last SimGroup (sG) of all ABM protocols
  # first collect and filter them to the sG simGroup in a new *.df file (faster);
  # the simGroup is in the perultimate column.
  system(paste("head -n 1 ",LogDirDate,"/Dataframe_1.df > ",LogDirDate,"/DataframeAfter",simGroups,".df",sep=""))
  system(paste("cat ",LogDirDate,"/Dataframe_*.df | grep \'\"",sG,"\" [0-9]*$\' | cut -d ' ' -f 2- >> ",LogDirDate,"/DataframeAfter",simGroups,".df",sep=""))
  read.table(paste(LogDirDate,"/DataframeAfter",simGroups,".df",sep=""),header=T,stringsAsFactors=F) -> res2.abm
  # calculate equiv.class labels: one agent, one time separately!
  eclassLabels = as.character(res2.abm$V)
  for(spkr in unique(res2.abm$Vpn)){
    for(abm in unique(res2.abm$ABM)){
      temp = res2.abm$Vpn == spkr & res2.abm$ABM == abm
      eclassLabels[temp] = equal_class(as.character(res2.abm$Initial[temp]), as.character(res2.abm$V[temp]))
    }
  }
  tab.abm = table(eclassLabels)
  col.abm = strongColors[phonClusterLog.names %in% names(tab.abm)]
  # 2 pie charts on one sheet
  if(runInteractive) {
    par(mfrow=c(1,2))
    pie(tab.orig, names(tab.orig),main="Equivalence Phone Labels before ABM",col=col.orig)
    pie(tab.abm, names(tab.abm),main="Equivalence Phone Labels after ABM",col=col.abm)
    par(mfrow=c(1,1))
    inpu = readline("Pie charts Equivalence Phone Labels before/after ABMs ok? Continue? (Y/n) ")
    if(inpu != "" && inpu != "Y" && inpu != "y") {
      cat("Script stopped by user\n")
      stop()
    }
  }
  if(printPaperFigs) {
    pdf(paste(LogDirDate,"catPhonClusterPieBeforeAfter.pdf",sep="/"))
    par(mfrow=c(1,2))
    pie(tab.orig, names(tab.orig),main="Equivalence Phone Labels before ABM",col=col.orig)
    pie(tab.abm, names(tab.abm),main="Equivalence Phone Labels after ABM",col=col.abm)
    par(mfrow=c(1,1))
    dev.off()
  }

  ######
  # STOP
  ######

}  # end multiple ABM



