From f4628a1c9f8a72bcd50f8a09d8e2fc0d95a749a4 Mon Sep 17 00:00:00 2001 From: "David R. Miller" Date: Sat, 18 Sep 2021 02:25:02 -0700 Subject: [PATCH] initial upload --- Makefile | 228 ++++++++++ README.md | 322 +++++++++++++ biosim4.cbp | 85 ++++ biosim4.ini | 208 +++++++++ src/analysis.cpp | 375 +++++++++++++++ src/basicTypes.cpp | 173 +++++++ src/basicTypes.h | 127 ++++++ src/biosim4.layout | 10 + src/createBarrier.cpp | 184 ++++++++ src/endOfGeneration.cpp | 36 ++ src/endOfSimStep.cpp | 92 ++++ src/executeActions.cpp | 245 ++++++++++ src/feedForward.cpp | 98 ++++ src/genome-compare.cpp | 175 +++++++ src/genome-neurons.h | 93 ++++ src/genome.cpp | 428 ++++++++++++++++++ src/getSensor.cpp | 357 +++++++++++++++ src/grid.cpp | 50 ++ src/grid.h | 71 +++ src/imageWriter.cpp | 260 +++++++++++ src/imageWriter.h | 55 +++ src/indiv.cpp | 30 ++ src/indiv.h | 40 ++ src/main.cpp | 28 ++ src/params.cpp | 300 ++++++++++++ src/params.h | 90 ++++ src/peeps.cpp | 82 ++++ src/peeps.h | 54 +++ src/random.cpp | 94 ++++ src/random.h | 39 ++ src/sensors-actions.h | 95 ++++ src/signals.cpp | 61 +++ src/signals.h | 51 +++ src/simulator.cpp | 171 +++++++ src/simulator.h | 53 +++ src/spawnNewGeneration.cpp | 195 ++++++++ src/survival-criteria.cpp | 345 ++++++++++++++ src/unitTestBasicTypes.cpp | 312 +++++++++++++ ...itTestConnectNeuralNetWiringFromGenome.cpp | 42 ++ src/unitTestGridVisitNeighborhood.cpp | 35 ++ tools/graph-nnet.py | 73 +++ tools/graphlog-genomeLength.gp | 28 ++ tools/graphlog.gp | 31 ++ 43 files changed, 5921 insertions(+) create mode 100644 Makefile create mode 100644 README.md create mode 100644 biosim4.cbp create mode 100644 biosim4.ini create mode 100644 src/analysis.cpp create mode 100644 src/basicTypes.cpp create mode 100644 src/basicTypes.h create mode 100644 src/biosim4.layout create mode 100644 src/createBarrier.cpp create mode 100644 src/endOfGeneration.cpp create mode 100644 src/endOfSimStep.cpp create mode 100644 src/executeActions.cpp create mode 100644 src/feedForward.cpp create mode 100644 src/genome-compare.cpp create mode 100644 src/genome-neurons.h create mode 100644 src/genome.cpp create mode 100644 src/getSensor.cpp create mode 100644 src/grid.cpp create mode 100644 src/grid.h create mode 100644 src/imageWriter.cpp create mode 100644 src/imageWriter.h create mode 100644 src/indiv.cpp create mode 100644 src/indiv.h create mode 100644 src/main.cpp create mode 100644 src/params.cpp create mode 100644 src/params.h create mode 100644 src/peeps.cpp create mode 100644 src/peeps.h create mode 100644 src/random.cpp create mode 100644 src/random.h create mode 100644 src/sensors-actions.h create mode 100644 src/signals.cpp create mode 100644 src/signals.h create mode 100644 src/simulator.cpp create mode 100644 src/simulator.h create mode 100644 src/spawnNewGeneration.cpp create mode 100644 src/survival-criteria.cpp create mode 100644 src/unitTestBasicTypes.cpp create mode 100644 src/unitTestConnectNeuralNetWiringFromGenome.cpp create mode 100644 src/unitTestGridVisitNeighborhood.cpp create mode 100644 tools/graph-nnet.py create mode 100644 tools/graphlog-genomeLength.gp create mode 100644 tools/graphlog.gp diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..082fadd --- /dev/null +++ b/Makefile @@ -0,0 +1,228 @@ +#------------------------------------------------------------------------------# +# This makefile was generated by 'cbp2make' tool rev.147 # +#------------------------------------------------------------------------------# + + +WORKDIR = `pwd` + +CC = gcc +CXX = g++ +AR = ar +LD = g++ +WINDRES = windres + +INC = -I/usr/include/opencv4 +CFLAGS = -Wall -fexceptions +RESINC = +LIBDIR = +LIB = /usr/lib/x86_64-linux-gnu/libopencv_core.so /usr/lib/x86_64-linux-gnu/libopencv_video.so /usr/lib/x86_64-linux-gnu/libopencv_videoio.so +LDFLAGS = -lX11 + +INC_DEBUG = $(INC) +CFLAGS_DEBUG = $(CFLAGS) -g -fopenmp +RESINC_DEBUG = $(RESINC) +RCFLAGS_DEBUG = $(RCFLAGS) +LIBDIR_DEBUG = $(LIBDIR) +LIB_DEBUG = $(LIB) +LDFLAGS_DEBUG = $(LDFLAGS) -fopenmp +OBJDIR_DEBUG = obj/Debug +DEP_DEBUG = +OUT_DEBUG = bin/Debug/biosim4 + +INC_RELEASE = $(INC) +CFLAGS_RELEASE = $(CFLAGS) -O3 +RESINC_RELEASE = $(RESINC) +RCFLAGS_RELEASE = $(RCFLAGS) +LIBDIR_RELEASE = $(LIBDIR) +LIB_RELEASE = $(LIB) +LDFLAGS_RELEASE = $(LDFLAGS) -O3 -s +OBJDIR_RELEASE = obj/Release +DEP_RELEASE = +OUT_RELEASE = bin/Release/biosim4 + +OBJ_DEBUG = $(OBJDIR_DEBUG)/src/signals.o $(OBJDIR_DEBUG)/src/main.o $(OBJDIR_DEBUG)/src/params.o $(OBJDIR_DEBUG)/src/peeps.o $(OBJDIR_DEBUG)/src/random.o $(OBJDIR_DEBUG)/src/simulator.o $(OBJDIR_DEBUG)/src/spawnNewGeneration.o $(OBJDIR_DEBUG)/src/survival-criteria.o $(OBJDIR_DEBUG)/src/unitTestBasicTypes.o $(OBJDIR_DEBUG)/src/unitTestConnectNeuralNetWiringFromGenome.o $(OBJDIR_DEBUG)/src/unitTestGridVisitNeighborhood.o $(OBJDIR_DEBUG)/src/genome-compare.o $(OBJDIR_DEBUG)/src/analysis.o $(OBJDIR_DEBUG)/src/basicTypes.o $(OBJDIR_DEBUG)/src/createBarrier.o $(OBJDIR_DEBUG)/src/endOfGeneration.o $(OBJDIR_DEBUG)/src/endOfSimStep.o $(OBJDIR_DEBUG)/src/executeActions.o $(OBJDIR_DEBUG)/src/feedForward.o $(OBJDIR_DEBUG)/src/genome.o $(OBJDIR_DEBUG)/src/getSensor.o $(OBJDIR_DEBUG)/src/grid.o $(OBJDIR_DEBUG)/src/imageWriter.o $(OBJDIR_DEBUG)/src/indiv.o + +OBJ_RELEASE = $(OBJDIR_RELEASE)/src/signals.o $(OBJDIR_RELEASE)/src/main.o $(OBJDIR_RELEASE)/src/params.o $(OBJDIR_RELEASE)/src/peeps.o $(OBJDIR_RELEASE)/src/random.o $(OBJDIR_RELEASE)/src/simulator.o $(OBJDIR_RELEASE)/src/spawnNewGeneration.o $(OBJDIR_RELEASE)/src/survival-criteria.o $(OBJDIR_RELEASE)/src/unitTestBasicTypes.o $(OBJDIR_RELEASE)/src/unitTestConnectNeuralNetWiringFromGenome.o $(OBJDIR_RELEASE)/src/unitTestGridVisitNeighborhood.o $(OBJDIR_RELEASE)/src/genome-compare.o $(OBJDIR_RELEASE)/src/analysis.o $(OBJDIR_RELEASE)/src/basicTypes.o $(OBJDIR_RELEASE)/src/createBarrier.o $(OBJDIR_RELEASE)/src/endOfGeneration.o $(OBJDIR_RELEASE)/src/endOfSimStep.o $(OBJDIR_RELEASE)/src/executeActions.o $(OBJDIR_RELEASE)/src/feedForward.o $(OBJDIR_RELEASE)/src/genome.o $(OBJDIR_RELEASE)/src/getSensor.o $(OBJDIR_RELEASE)/src/grid.o $(OBJDIR_RELEASE)/src/imageWriter.o $(OBJDIR_RELEASE)/src/indiv.o + +all: debug release + +clean: clean_debug clean_release + +before_debug: + test -d bin/Debug || mkdir -p bin/Debug + test -d $(OBJDIR_DEBUG)/src || mkdir -p $(OBJDIR_DEBUG)/src + +after_debug: + +debug: before_debug out_debug after_debug + +out_debug: before_debug $(OBJ_DEBUG) $(DEP_DEBUG) + $(LD) $(LIBDIR_DEBUG) -o $(OUT_DEBUG) $(OBJ_DEBUG) $(LDFLAGS_DEBUG) $(LIB_DEBUG) + +$(OBJDIR_DEBUG)/src/signals.o: src/signals.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/signals.cpp -o $(OBJDIR_DEBUG)/src/signals.o + +$(OBJDIR_DEBUG)/src/main.o: src/main.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/main.cpp -o $(OBJDIR_DEBUG)/src/main.o + +$(OBJDIR_DEBUG)/src/params.o: src/params.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/params.cpp -o $(OBJDIR_DEBUG)/src/params.o + +$(OBJDIR_DEBUG)/src/peeps.o: src/peeps.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/peeps.cpp -o $(OBJDIR_DEBUG)/src/peeps.o + +$(OBJDIR_DEBUG)/src/random.o: src/random.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/random.cpp -o $(OBJDIR_DEBUG)/src/random.o + +$(OBJDIR_DEBUG)/src/simulator.o: src/simulator.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/simulator.cpp -o $(OBJDIR_DEBUG)/src/simulator.o + +$(OBJDIR_DEBUG)/src/spawnNewGeneration.o: src/spawnNewGeneration.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/spawnNewGeneration.cpp -o $(OBJDIR_DEBUG)/src/spawnNewGeneration.o + +$(OBJDIR_DEBUG)/src/survival-criteria.o: src/survival-criteria.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/survival-criteria.cpp -o $(OBJDIR_DEBUG)/src/survival-criteria.o + +$(OBJDIR_DEBUG)/src/unitTestBasicTypes.o: src/unitTestBasicTypes.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/unitTestBasicTypes.cpp -o $(OBJDIR_DEBUG)/src/unitTestBasicTypes.o + +$(OBJDIR_DEBUG)/src/unitTestConnectNeuralNetWiringFromGenome.o: src/unitTestConnectNeuralNetWiringFromGenome.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/unitTestConnectNeuralNetWiringFromGenome.cpp -o $(OBJDIR_DEBUG)/src/unitTestConnectNeuralNetWiringFromGenome.o + +$(OBJDIR_DEBUG)/src/unitTestGridVisitNeighborhood.o: src/unitTestGridVisitNeighborhood.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/unitTestGridVisitNeighborhood.cpp -o $(OBJDIR_DEBUG)/src/unitTestGridVisitNeighborhood.o + +$(OBJDIR_DEBUG)/src/genome-compare.o: src/genome-compare.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/genome-compare.cpp -o $(OBJDIR_DEBUG)/src/genome-compare.o + +$(OBJDIR_DEBUG)/src/analysis.o: src/analysis.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/analysis.cpp -o $(OBJDIR_DEBUG)/src/analysis.o + +$(OBJDIR_DEBUG)/src/basicTypes.o: src/basicTypes.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/basicTypes.cpp -o $(OBJDIR_DEBUG)/src/basicTypes.o + +$(OBJDIR_DEBUG)/src/createBarrier.o: src/createBarrier.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/createBarrier.cpp -o $(OBJDIR_DEBUG)/src/createBarrier.o + +$(OBJDIR_DEBUG)/src/endOfGeneration.o: src/endOfGeneration.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/endOfGeneration.cpp -o $(OBJDIR_DEBUG)/src/endOfGeneration.o + +$(OBJDIR_DEBUG)/src/endOfSimStep.o: src/endOfSimStep.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/endOfSimStep.cpp -o $(OBJDIR_DEBUG)/src/endOfSimStep.o + +$(OBJDIR_DEBUG)/src/executeActions.o: src/executeActions.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/executeActions.cpp -o $(OBJDIR_DEBUG)/src/executeActions.o + +$(OBJDIR_DEBUG)/src/feedForward.o: src/feedForward.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/feedForward.cpp -o $(OBJDIR_DEBUG)/src/feedForward.o + +$(OBJDIR_DEBUG)/src/genome.o: src/genome.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/genome.cpp -o $(OBJDIR_DEBUG)/src/genome.o + +$(OBJDIR_DEBUG)/src/getSensor.o: src/getSensor.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/getSensor.cpp -o $(OBJDIR_DEBUG)/src/getSensor.o + +$(OBJDIR_DEBUG)/src/grid.o: src/grid.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/grid.cpp -o $(OBJDIR_DEBUG)/src/grid.o + +$(OBJDIR_DEBUG)/src/imageWriter.o: src/imageWriter.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/imageWriter.cpp -o $(OBJDIR_DEBUG)/src/imageWriter.o + +$(OBJDIR_DEBUG)/src/indiv.o: src/indiv.cpp + $(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c src/indiv.cpp -o $(OBJDIR_DEBUG)/src/indiv.o + +clean_debug: + rm -f $(OBJ_DEBUG) $(OUT_DEBUG) + rm -rf bin/Debug + rm -rf $(OBJDIR_DEBUG)/src + +before_release: + test -d bin/Release || mkdir -p bin/Release + test -d $(OBJDIR_RELEASE)/src || mkdir -p $(OBJDIR_RELEASE)/src + +after_release: + +release: before_release out_release after_release + +out_release: before_release $(OBJ_RELEASE) $(DEP_RELEASE) + $(LD) $(LIBDIR_RELEASE) -o $(OUT_RELEASE) $(OBJ_RELEASE) $(LDFLAGS_RELEASE) $(LIB_RELEASE) + +$(OBJDIR_RELEASE)/src/signals.o: src/signals.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/signals.cpp -o $(OBJDIR_RELEASE)/src/signals.o + +$(OBJDIR_RELEASE)/src/main.o: src/main.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/main.cpp -o $(OBJDIR_RELEASE)/src/main.o + +$(OBJDIR_RELEASE)/src/params.o: src/params.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/params.cpp -o $(OBJDIR_RELEASE)/src/params.o + +$(OBJDIR_RELEASE)/src/peeps.o: src/peeps.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/peeps.cpp -o $(OBJDIR_RELEASE)/src/peeps.o + +$(OBJDIR_RELEASE)/src/random.o: src/random.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/random.cpp -o $(OBJDIR_RELEASE)/src/random.o + +$(OBJDIR_RELEASE)/src/simulator.o: src/simulator.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/simulator.cpp -o $(OBJDIR_RELEASE)/src/simulator.o + +$(OBJDIR_RELEASE)/src/spawnNewGeneration.o: src/spawnNewGeneration.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/spawnNewGeneration.cpp -o $(OBJDIR_RELEASE)/src/spawnNewGeneration.o + +$(OBJDIR_RELEASE)/src/survival-criteria.o: src/survival-criteria.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/survival-criteria.cpp -o $(OBJDIR_RELEASE)/src/survival-criteria.o + +$(OBJDIR_RELEASE)/src/unitTestBasicTypes.o: src/unitTestBasicTypes.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/unitTestBasicTypes.cpp -o $(OBJDIR_RELEASE)/src/unitTestBasicTypes.o + +$(OBJDIR_RELEASE)/src/unitTestConnectNeuralNetWiringFromGenome.o: src/unitTestConnectNeuralNetWiringFromGenome.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/unitTestConnectNeuralNetWiringFromGenome.cpp -o $(OBJDIR_RELEASE)/src/unitTestConnectNeuralNetWiringFromGenome.o + +$(OBJDIR_RELEASE)/src/unitTestGridVisitNeighborhood.o: src/unitTestGridVisitNeighborhood.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/unitTestGridVisitNeighborhood.cpp -o $(OBJDIR_RELEASE)/src/unitTestGridVisitNeighborhood.o + +$(OBJDIR_RELEASE)/src/genome-compare.o: src/genome-compare.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/genome-compare.cpp -o $(OBJDIR_RELEASE)/src/genome-compare.o + +$(OBJDIR_RELEASE)/src/analysis.o: src/analysis.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/analysis.cpp -o $(OBJDIR_RELEASE)/src/analysis.o + +$(OBJDIR_RELEASE)/src/basicTypes.o: src/basicTypes.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/basicTypes.cpp -o $(OBJDIR_RELEASE)/src/basicTypes.o + +$(OBJDIR_RELEASE)/src/createBarrier.o: src/createBarrier.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/createBarrier.cpp -o $(OBJDIR_RELEASE)/src/createBarrier.o + +$(OBJDIR_RELEASE)/src/endOfGeneration.o: src/endOfGeneration.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/endOfGeneration.cpp -o $(OBJDIR_RELEASE)/src/endOfGeneration.o + +$(OBJDIR_RELEASE)/src/endOfSimStep.o: src/endOfSimStep.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/endOfSimStep.cpp -o $(OBJDIR_RELEASE)/src/endOfSimStep.o + +$(OBJDIR_RELEASE)/src/executeActions.o: src/executeActions.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/executeActions.cpp -o $(OBJDIR_RELEASE)/src/executeActions.o + +$(OBJDIR_RELEASE)/src/feedForward.o: src/feedForward.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/feedForward.cpp -o $(OBJDIR_RELEASE)/src/feedForward.o + +$(OBJDIR_RELEASE)/src/genome.o: src/genome.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/genome.cpp -o $(OBJDIR_RELEASE)/src/genome.o + +$(OBJDIR_RELEASE)/src/getSensor.o: src/getSensor.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/getSensor.cpp -o $(OBJDIR_RELEASE)/src/getSensor.o + +$(OBJDIR_RELEASE)/src/grid.o: src/grid.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/grid.cpp -o $(OBJDIR_RELEASE)/src/grid.o + +$(OBJDIR_RELEASE)/src/imageWriter.o: src/imageWriter.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/imageWriter.cpp -o $(OBJDIR_RELEASE)/src/imageWriter.o + +$(OBJDIR_RELEASE)/src/indiv.o: src/indiv.cpp + $(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c src/indiv.cpp -o $(OBJDIR_RELEASE)/src/indiv.o + +clean_release: + rm -f $(OBJ_RELEASE) $(OUT_RELEASE) + rm -rf bin/Release + rm -rf $(OBJDIR_RELEASE)/src + +.PHONY: before_debug after_debug clean_debug before_release after_release clean_release + diff --git a/README.md b/README.md new file mode 100644 index 0000000..b9a5c0d --- /dev/null +++ b/README.md @@ -0,0 +1,322 @@ +# biosim4 + +## What is this? + +This pile of code was used to simulate biological creatures that evolve through natural selection. +The results of the experiments are summarized in this YouTube video: + +      "I programmed some creatures. They evolved." + +      https://www.youtube.com/watch?v=N3tRFayqVtk + +This code lacks a friendly interface, so compiling and executing the program may +require attention to details. If you ask questions in the Issues, +I'll try to help if I can. + +Document Contents +----------------- + +* [Code walkthrough](#CodeWalkthrough) +* [Main data structures](#MainDataStructures) +* [Config file](#ConfigFile) +* [Program output](#ProgramOutput) +* [Main program loop](#MainProgramLoop) +* [Sensory inputs and action outputs](#SensoryInputsAndActionOutputs) +* [Basic value types](#BasicValueTypes) +* [Pheromones](#Pheromones) +* [Useful utility functions](#UsefulUtilityFunctions) +* [Installing the code](#InstallingTheCode) +* [Building the executable](#BuildingTheExecutable) +* [System requirements](#SystemRequirements) +* [Compiling](#Compiling) +* [Bugs](#Bugs) +* [Execution](#Execution) +* [Tools directory](#ToolsDirectory) +* [Build log](#BuildLog) + + +Code walkthrough +-------------------- + + +### Main data structures + +The code in the src directory compiles to a single console program named biosim4. When it is +invoked, it will read parameters from a config file named biosim4.ini by default. A different +config file can be specified on the command line. + +The simulator will then configure a 2D arena where the creatures live. Class Grid (see grid.h and grid.cpp) +contains a 2D array of 16-bit indexes, where each nonzero index refers to a specific individual in class Peeps (see below). +Zero values in Grid indicate empty locations. Class Grid does not know anything else about the world; it only +stores indexes to represent who lives where. + +The population of creatures is stored in class Peeps (see peeps.h and peeps.cpp). Class Peeps contains +all the individuals in the simulation, stored as instances of struct Indiv in a std::vector container. +The indexes in class Grid are indexes into the vector of individuals in class Peeps. Class Peeps keeps a +container of struct Indiv, but otherwise does not know anything about the internal workings of individuals. + +Each individual is represented by an instance of struct Indiv (see indiv.h and indiv.cpp). Struct Indiv +contains an individual's genome, its corresponding neural net brain, and a redundant copy of the individual's +X,Y location in the 2D grid. It also contains a few other parameters for the individual, such as its +"responsiveness" level, oscillator period, age, and other personal parameters. Struct Indiv knows how +to convert an individual's genome into its neural net brain at the beginning of the simulation. +It also knows how to print the genome and neural net brain in text format to stdout during a simulation. +It also has a function Indiv::getSensor() that is called to compute the individual's input neurons for +each simulator step. + +All the simulator code lives in the BS namespace (short for "biosim".) + + +### Config file + +The config file, named biosim4.ini by default, contains all the tunable parameters for a +simulation run. The biosim4 executable reads the config file at startup, then monitors it for +changes during the simulation. Although it's not foolproof, many parameters can be modified during +the simulation run. Class ParamManager (see params.h and params.cpp) manages the configuration +parameters and makes them available to the simulator through a read-only pointer provided by +ParamManager::getParamRef(). + +See the provided biosim4.ini for documentation for each parameter. Most of the parameters +in the config file correspond to members in struct Params (see params.h). A few additional +parameters may be stored in struct Params. See the documentation in params.h for how to +support new parameters. + + + +### Program output + +Depending on the parameters in the config file, the following data can be produced: + +* The simulator will append one line to logs/epoch.txt after the completion of +each generation. Each line records the generation number, number of individuals +who survived the selection criterion, an estimate of the population's genetic +diversity, average genome length, and number of deaths due to the "kill" gene. +This file can be fed to tools/graphlog.gp to produce a graphic plot. + +* The simulator will display a small number of sample genomes at regular +intervals to stdout. Parameters in the config file specify the number and interval. +The genomes are displayed in hex format and also in a mnemonic format that can +be fed to tools/graph-nnet.py to produce a graphic network diagram. + +* Movies of selected generations will be created in the images/ directory. Parameters +in the config file specify the interval at which to make movies. Each movie records +a single generation. + +* At intervals, a summary is printed to stdout showing the total number of neural +connections throughout the population from each possible sensory input neuron and to each +possible action output neuron. + + +### Main program loop + +The simulator starts with a call to simulator() in simulator.cpp. After initializing the +world, the simulator executes three nested loops: the outer loop for each generation, +an inner loop for each simulator step within the generation, and an innermost loop for +each individual in the population. The innermost loop is thread-safe so that it can +be parallelized by OpenMP. + +At the end of each simulator step, a call is made to endOfSimStep() in single-thread +mode (see endOfSimStep.cpp) to create a video frame representing the locations of all +the individuals at the end of the simulator step. The video frame is pushed on to a +stack to be converted to a movie later. Also some housekeeping may be done for certain +selection scenarios. See the comments in endOfSimStep.cpp for more information. + +At the end of each generation, a call is made to endOfGeneration() in single-thread +mode (see endOfGeneration.cpp) to create a video from the saved video frames. +Also a new graph might be generated showing the progress of the simulation. See +endOfGeneraton.cpp for more information. + + +### Sensory inputs and action outputs + +See the YouTube video (link above) for a description of the sensory inputs and action +outputs. Each sensory input and each action output is a neuron in the individual's +neural net brain. + +The header file sensors-actions.h contains enum Sensor which enumerates all the possible sensory +inputs and enum Action which enumerates all the possible action outputs. +In enum Sensor, all the sensory inputs before the enumerant NUM_SENSES will +be compiled into the executable, and all action outputs before NUM_ACTIONS +will be compiled. By rearranging the enumerants in those enums, you can select +a subset of all possible sensory and action neurons to be compiled into the +simulator. + + +### Basic value types + +There are a few basic value types: + +* enum Compass represents eight-way directions with enumerants N=0, NE, E, SW, S, SW, W, NW, CENTER. + +* struct Dir is an abstract representation of the values of enum Compass. + +* struct Coord is a signed 16-bit integer X,Y coordinate pair. It is used to represent a location +in the 2D world, or can represent the difference between two locations. + +* struct Polar holds a signed 32-bit integer magnitude and a direction of type Dir. + +Various conversions and math are possible between these basic types. See unitTestBasicTypes.cpp +for examples. Also see basicTypes.h for more information. + + +### Pheromones + +A simple system is used to simulate pheromones emitted by the individuals. Pheromones +are called "signals" in simulator-speak (see signals.h and signals.cpp). Struct Signals +holds a single layer that overlays the 2D world in class Grid. Each location can contain +a level of pheromone (there's only a single kind of pheromone supported at present). The +pheromone level at any grid location is stored as an unsigned 8-bit integer, where zero means no +pheromone, and 255 is the maximum. Each time an individual emits a pheromone, it increases +the pheromone values in a small neighborhood around the individual up to the maximum +value of 255. Pheromone levels decay over time if they are not replenished +by the individuals in the area. + + +### Useful utility functions + +The utility function visitNeighborhood() in grid.cpp can be used to execute a +user-defined lambda or function over each location +within a circular neighborhood defined by a center point and floating point radius. The function +calls the user-defined function once for each location, passing it a Coord value. Only locations +within the bounds of the grid are visited. The center location is included among the visited +locations. For example, a radius of 1.0 includes only the center location plus four neighboring locations. +A radius of 1.5 includes the center plus all the eight-way neighbors. The radius can be arbitrarily large +but large radii require lots of CPU cycles. + + + + +## Installing the code +-------------------- + +Copy the directory structure to a location of your choice. + + +## Building the executable +-------------------- + + +### System requirements + +This code is known to run in the following environment: + +* Ubuntu 21.04 +* cimg-dev 2.8.4 or later +* libopencv-dev 4.2 or later +* gcc 9.3 or 10.3 +* python-igraph 0.8.3 (used only by tools/graph-nnet.py) +* gnuplot 5.2.8 (used only by tools/graphlog.gp) + +The code also runs in distributions based on Ubuntu 20.04, but only if the default version of +cimg-dev is replaced with version 2.8.4 or later. + + +### Compiling + +Two ways to compile: + +* The file named "biosim4.cbp" is a configuration file for the Code::Blocks IDE version 20.03. + +* A Makefile is provided which +was created from biosim4.cbp with cbp2make, but is not tested. A default "make" will generate a debug and a +release version. + + +## Bugs +-------------------- + +If you try to compile the simulator under a distribution based on Ubuntu 20.04, you will encounter this +bug in the version of CImg.h (package cimg-dev) provided by the package maintainer: + +      https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=951965 + +In biosim4, CImg.h is used only as a convenient interface to OpenCV +to generate movies of the simulated creatures in their 2D world. You have several +choices if you want to proceed with Ubuntu 20.04: + +* You can strip out the code that generates the movies and just run the simulator without the movies. Most of +that graphics code is in imageWriter.cpp and imageWriter.h. + +* You can upgrade your CImg.h to version 2.8.4 or later by getting it from the appropriate Debian repository. +Sorry I don't have the instructions at hand to do this. + +* You could convert the CImg.h function calls to use OpenCV directly. Sorry I don't have a guide for how +to do that. + + +## Execution +-------------------- + +Edit the config file (default "biosim4.ini") for the parameters you want for the simulation run, then execute the Debug +or Release executable in the bin directory. Optionally specify the name of the config file as the first +command line argument, e.g.: + +``` +./bin/Release/biosim4 [biosim4.ini] +``` + + +## Tools directory +-------------------- + +tools/graphlog.gp takes the generated log file logs/epoch-log.txt +and generates a graphic plot of the simulation run in images/log.png. You may need to adjust +the directory paths in graphlog.gp for your environment. graphlog.gp can be invoked manually, +or if the option "updateGraphLog" is set to true +in the simulation config file, the simulator will try to invoke tools/graphlog.gp automatically +during the simulation run. Also see the parameter named updateGraphLogStride in the config file. + +tools/graph-nnet.py takes a text file (hardcoded name "net.txt") and generates a neural net +connection diagram using igraph. The file net.txt contains an encoded form of one genome, and +must be the same format as the files +generated by displaySampleGenomes() in src/analysis.cpp which is called by simulator() in +src/simulator.cpp. The genome output is printed to stdout automatically +if the parameter named "displaySampleGenomes" is set to nonzero in the config file. +An individual genome can be copied from that output stream and renamed "net.txt" in order to run +graph-nnet.py. + + +## Build log +-------------------- + +In case it helps for debugging the build process, here is a build log from Code::Blocks running under Ubuntu 21.04: + + +``` +-------------- Clean: Release in biosim4 (compiler: GNU GCC Compiler)--------------- + +Cleaned "biosim4 - Release" + +-------------- Build: Release in biosim4 (compiler: GNU GCC Compiler)--------------- + +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/analysis.cpp -o obj/Release/src/analysis.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/basicTypes.cpp -o obj/Release/src/basicTypes.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/createBarrier.cpp -o obj/Release/src/createBarrier.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/endOfGeneration.cpp -o obj/Release/src/endOfGeneration.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/endOfSimStep.cpp -o obj/Release/src/endOfSimStep.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/executeActions.cpp -o obj/Release/src/executeActions.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/feedForward.cpp -o obj/Release/src/feedForward.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/genome-compare.cpp -o obj/Release/src/genome-compare.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/genome.cpp -o obj/Release/src/genome.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/getSensor.cpp -o obj/Release/src/getSensor.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/grid.cpp -o obj/Release/src/grid.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/imageWriter.cpp -o obj/Release/src/imageWriter.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/indiv.cpp -o obj/Release/src/indiv.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/main.cpp -o obj/Release/src/main.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/params.cpp -o obj/Release/src/params.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/peeps.cpp -o obj/Release/src/peeps.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/random.cpp -o obj/Release/src/random.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/signals.cpp -o obj/Release/src/signals.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/simulator.cpp -o obj/Release/src/simulator.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/spawnNewGeneration.cpp -o obj/Release/src/spawnNewGeneration.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/survival-criteria.cpp -o obj/Release/src/survival-criteria.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/unitTestBasicTypes.cpp -o obj/Release/src/unitTestBasicTypes.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/unitTestConnectNeuralNetWiringFromGenome.cpp -o obj/Release/src/unitTestConnectNeuralNetWiringFromGenome.o +g++ -Wall -fexceptions -fopenmp -O3 -I/usr/include/opencv4 -c /home/dm/sw/biosim4-git/src/unitTestGridVisitNeighborhood.cpp -o obj/Release/src/unitTestGridVisitNeighborhood.o +g++ -o bin/Release/biosim4 obj/Release/src/analysis.o obj/Release/src/basicTypes.o obj/Release/src/createBarrier.o obj/Release/src/endOfGeneration.o obj/Release/src/endOfSimStep.o obj/Release/src/executeActions.o obj/Release/src/feedForward.o obj/Release/src/genome-compare.o obj/Release/src/genome.o obj/Release/src/getSensor.o obj/Release/src/grid.o obj/Release/src/imageWriter.o obj/Release/src/indiv.o obj/Release/src/main.o obj/Release/src/params.o obj/Release/src/peeps.o obj/Release/src/random.o obj/Release/src/signals.o obj/Release/src/simulator.o obj/Release/src/spawnNewGeneration.o obj/Release/src/survival-criteria.o obj/Release/src/unitTestBasicTypes.o obj/Release/src/unitTestConnectNeuralNetWiringFromGenome.o obj/Release/src/unitTestGridVisitNeighborhood.o -lX11 -lgomp -pthread -O3 -s /usr/lib/x86_64-linux-gnu/libopencv_core.so /usr/lib/x86_64-linux-gnu/libopencv_video.so /usr/lib/x86_64-linux-gnu/libopencv_videoio.so +Output file is bin/Release/biosim4 with size 778.42 KB +Process terminated with status 0 (0 minute(s), 11 second(s)) +0 error(s), 0 warning(s) (0 minute(s), 11 second(s)) +``` + + diff --git a/biosim4.cbp b/biosim4.cbp new file mode 100644 index 0000000..20a38bf --- /dev/null +++ b/biosim4.cbp @@ -0,0 +1,85 @@ + + + + + + diff --git a/biosim4.ini b/biosim4.ini new file mode 100644 index 0000000..2ae06c5 --- /dev/null +++ b/biosim4.ini @@ -0,0 +1,208 @@ +# biosim4.ini +# biosim4.ini is the default config file for the simluator. +# The config filename is determined in simulator() in simulator.cpp. +# The config file is parsed by class ParamManager, see params.cpp and params.h. +# Although not foolproof, the config file can be modifed during a simulation +# run and the param manager will make any new params available to the simulator +# after the end of the current simulator step or after the end of the current +# generation. + +# numThreads must be 1 or greater. Best value is less than or equal to +# the number of CPU cores. +numThreads = 10 + +# sizeX, sizeY define the size of the 2D world. Minimum size is 16,16. +# Maximum size is 32767, 32767. +sizeX = 128 +sizeY = 128 + +# Population at the start of each generation. Maximum value = 32766. +population = 3000 + +# Number of simulation steps per generation. Range 1..INT_MAX. +stepsPerGeneration = 300 + +# The simulator will stop when the generation number == maxGenerations. +# Range 1..INT_MAX +maxGenerations = 2000000 + +# genomeInitialLengthMin and genomeInitialLengthMax should be set to +# the same value. (For future use, the max length might be larger to +# allow mutations that lengthen the genome.) Range 1..INT_MAX and +# must be no larger than genomeMaxLength. The range of genomeMaxLength +# is genomeInitialLengthMax..INT_MAX. +genomeInitialLengthMin = 48 +genomeInitialLengthMax = 48 +genomeMaxLength = 300 + +# maxNumberNeurons is the maximum number of internal neurons that may +# be addressed by genes in the genome. Range 1..INT_MAX. +maxNumberNeurons = 5 + +# If killEnable is true and the "kill" action neuron is enabled in +# sensors-actions.h and compiled in, then agents are permitted to +# kill their neighbor in the adjacent location in the direction of +# forward movement. If false, the neighbors are safe. +killEnable = false + +# If sexualReproduction is false, newborns inherit the genes from a +# single parent. If true, newborns inherit a mixture of genes from +# two parents. +sexualReproduction = true + +# If chooseParentByFitness is false, then every agent that survives the +# selection criterion has equal chance of reproducing. If true, then +# preference is given to those parents who passed the selection criterion +# with a greater score. Fitness scores are determined in survival-criteria.cpp. +chooseParentsByFitness = true + +# pointMutationRate is the probability per gene of having a single-bit +# mutation during spawning. Range 0.0 .. 1.0. A reasonable range is +# 0.0001 to 0.01. +pointMutationRate = 0.001 + +# geneInsertionDeletionRate and deletionRatio are for future use to +# allow mutations that lengthen or shorten the genome. Ignored for now. +geneInsertionDeletionRate = 0.0 +deletionRatio = 0.5 + +# responsivenessCurveKFactor is a small positive integer that determines +# the shape of the curve that determines how reactive an agent is to its +# sensory inputs. Typical values are # 1, 2, 3, or 4, but greater values +# are allowed experimentally. +responsivenessCurveKFactor = 2 + +# populationSensorRadius is the radius in which the population sensor +# looks for neighbors. Floating point value. A value of 1.5 includes +# all the immediate eight-neighborhood. Larger values incur exponentially +# increasing processor overhead. Range 0.5 up to (float)max(sizeX, sizeY). +populationSensorRadius = 2.5 + +# longProbeDistance is the default distance that the long-probe sensors +# are able to see. Applies to long-probe population sensor and long-probe +# signal (pheromone) sensor. Range 1..INT_MAX. +longProbeDistance = 24 + +# shortProbeBarrierDistance is the distance that the short-probe sensor +# can see. Range 1..INT_MAX. +shortProbeBarrierDistance = 4 + +# signalSensorRadius is the radius in which the signal (pheromone) sensor +# looks for pheromones. Floating point value. A value of 1.5 includes +# all the immediate eight-neighborhood. Larger values incur exponentially +# increasing processor overhead. Range 0.5 up to (float)max(sizeX, sizeY). +signalSensorRadius = 2.0 + +# signalLayers defines the number of pheromone layers. Must be 1 for now. +# Values > 1 are for future use. +signalLayers = 1 + +# imageDir is the relative or absolute directory path where generation +# movies are created. +imageDir = images + +# logDir is the relative or absolute directory path where text log files +# are created. +logDir = logs + +# displayScale scales the generation movie. Typical values are +# 1 for actual size, or 2, 4, 8, 16, or 32 to scale up the movie. +displayScale = 8 + +# agentSize controls the size of the dot used to represent an agent +# in the generation movie. Typical value is displayScale / 2. +agentSize = 4 + +# If videoSaveFirstFrames is 0, then only the parameter videoStride controls +# how often generation movies are made. If videoSaveFirstFrames is nonzero, +# then generation movies will also be generated for every generation from 0 +# through videoSaveFirstFrames (because the first few generations are often +# the most interesting). Range 1..INT_MAX. +videoSaveFirstFrames = 2 + +# updateGraphLog can be set to true to cause the simulator program to +# invoke graphlog.gp to update the simulation progress graph. If true, +# then updateGraphLogStride controls how often it is invoked. If +# updateGraphLog is false, then the simulator program will not invoke +# graphlog.gp. +updateGraphLog = true + +# If saveVideo is true, the simulator program will create generation +# movies in the directory named by imageDir at the intervals set by +# videoSaveFirstFrames and videoStride. +saveVideo = true + +# videoStride determines how often generation movies will be created. +# Also see saveVideo and videoSaveFirstFrames. Range 1..INT_MAX. +videoStride = 25 + +# updateGraphLogStride determines how often the simulation progress graph +# is updated by direct invocation of graphlog.gp. Ignored if updateGraphLog +# is false. updateGraphLogStride may be a positive integer from 1 to INT_MAX, +# or may be set to the string videoStride to use the value of videoStride. +updateGraphLogStride = videoStride + +# genomeAnalysisStride determines how often the simulator will print genomic +# statistics. The stats are printed to stdout when the generation number +# modulo genomeAnalysisStride == 0. The value may be a positive integer from +# 1 to INT_MAX, or may be set to the string videoStride to use the value of +# videoStride. +genomeAnalysisStride = videoStride + +# When the genomic statistics are printed (see genomeAnalysisStride), the +# method used to measure genome diversity in the population is determined +# by genomeComparisonMethod. May be set to 0 for Jaro-Winkler method (useful +# for future use if genomes are allowed to grow or shrink in size); or 1 +# for a Hamming measure bit-by-bit, or 2 for a Hamming measure byte-by-byte. +# Typically set to 1. +genomeComparisonMethod = 1 + +# When genomic statistics are printed (see genomeAnalysisStride), the number +# of genomes sampled from the population and printed to stdout is determined +# by displaySampleGenomes. Range 0 to population size. +displaySampleGenomes = 10 + +# challenge determines the selection criterion for reproduction. This is +# typically always under active development. See survival-criteria.cpp for +# more information. +# 0 = circle +# 1 = right half +# 2 = right quarter +# 3 = neighbor count +# 4 = center weighted +# 40 = center unweighted +# 5 = corner 1/4 radius +# 6 = corner 1/4 radius weighted +# 7 = migrate distance +# 8 = center sparse +# 9 = left eighth +# 10 = radioactive walls +# 11 = against any wall +# 12 = touch any wall any time +# 13 = east-west eighths +# 14 = near barrier +# 15 = pairs +# 16 = contact location sequence +# 17 = altruism, circle + NE corner +challenge = 17 + +# The simulator supports a feature called "barriers." Barriers are locations +# in the simulated 2D world where agents may not occupy. The value of +# barrierType is typically under active development. See createBarrier.cpp +# for more information. Also see replaceBarrierType. +# 0 = none +# 1 = vertical bar constant location +# 2 = vertical bar random locations +# 3 = five staggered vertical bars +# 4 = horiz bar constant location north center +# 5 = floating islands +# 6 = sequence of spots +barrierType = 0 + +# The barrier type will automatically change from barrierType to +# replaceBarrierType at the generation specified by +# replaceBarrierTypeGenerationNumber. Range 0..INT_MAX, or set to -1 to +# disable. +replaceBarrierType = 0 +replaceBarrierTypeGenerationNumber = -1 + diff --git a/src/analysis.cpp b/src/analysis.cpp new file mode 100644 index 0000000..cc8f7e6 --- /dev/null +++ b/src/analysis.cpp @@ -0,0 +1,375 @@ +// analysis.cpp -- various reports + +#include +#include +#include +#include +#include +#include +#include "simulator.h" + +namespace BS { + +// This converts sensor numbers to descriptive strings. +std::string sensorName(Sensor sensor) +{ + switch(sensor) { + case AGE: return "age"; break; + case BOUNDARY_DIST: return "boundary dist"; break; + case BOUNDARY_DIST_X: return "boundary dist X"; break; + case BOUNDARY_DIST_Y: return "boundary dist Y"; break; + case LAST_MOVE_DIR_X: return "last move dir X"; break; + case LAST_MOVE_DIR_Y: return "last move dir Y"; break; + case LOC_X: return "loc X"; break; + case LOC_Y: return "loc Y"; break; + case LONGPROBE_POP_FWD: return "long probe population fwd"; break; + case LONGPROBE_BAR_FWD: return "long probe barrier fwd"; break; + case BARRIER_FWD: return "short probe barrier fwd-rev"; break; + case BARRIER_LR: return "short probe barrier left-right"; break; + case OSC1: return "osc1"; break; + case POPULATION: return "population"; break; + case POPULATION_FWD: return "population fwd"; break; + case POPULATION_LR: return "population LR"; break; + case RANDOM: return "random"; break; + case SIGNAL0: return "signal 0"; break; + case SIGNAL0_FWD: return "signal 0 fwd"; break; + case SIGNAL0_LR: return "signal 0 LR"; break; + case GENETIC_SIM_FWD: return "genetic similarity fwd"; break; + default: assert(false); break; + } +} + + +// Converts action numbers to descriptive strings. +std::string actionName(Action action) +{ + switch(action) { + case MOVE_EAST: return "move east"; break; + case MOVE_WEST: return "move west"; break; + case MOVE_NORTH: return "move north"; break; + case MOVE_SOUTH: return "move south"; break; + case MOVE_FORWARD: return "move fwd"; break; + case MOVE_X: return "move X"; break; + case MOVE_Y: return "move Y"; break; + case SET_RESPONSIVENESS: return "set inv-responsiveness"; break; + case SET_OSCILLATOR_PERIOD: return "set osc1"; break; + case EMIT_SIGNAL0: return "emit signal 0"; break; + case KILL_FORWARD: return "kill fwd"; break; + case MOVE_REVERSE: return "move reverse"; break; + case MOVE_LEFT: return "move left"; break; + case MOVE_RIGHT: return "move right"; break; + case MOVE_RL: return "move R-L"; break; + case MOVE_RANDOM: return "move random"; break; + case SET_LONGPROBE_DIST: return "set longprobe dist"; break; + default: assert(false); break; + } +} + + +// This converts sensor numbers to mnemonic strings. +// Useful for later processing by graph-nnet.py. +std::string sensorShortName(Sensor sensor) +{ + switch(sensor) { + case AGE: return "Age"; break; + case BOUNDARY_DIST: return "ED"; break; + case BOUNDARY_DIST_X: return "EDx"; break; + case BOUNDARY_DIST_Y: return "EDy"; break; + case LAST_MOVE_DIR_X: return "LMx"; break; + case LAST_MOVE_DIR_Y: return "LMy"; break; + case LOC_X: return "Lx"; break; + case LOC_Y: return "Ly"; break; + case LONGPROBE_POP_FWD: return "LPf"; break; + case LONGPROBE_BAR_FWD: return "LPb"; break; + case BARRIER_FWD: return "Bfd"; break; + case BARRIER_LR: return "Blr"; break; + case OSC1: return "Osc"; break; + case POPULATION: return "Pop"; break; + case POPULATION_FWD: return "Pfd"; break; + case POPULATION_LR: return "Plr"; break; + case RANDOM: return "Rnd"; break; + case SIGNAL0: return "Sg"; break; + case SIGNAL0_FWD: return "Sfd"; break; + case SIGNAL0_LR: return "Slr"; break; + case GENETIC_SIM_FWD: return "Gen"; break; + default: assert(false); break; + } +} + + +// Converts action numbers to mnemonic strings. +// Useful for later processing by graph-nnet.py. +std::string actionShortName(Action action) +{ + switch(action) { + case MOVE_EAST: return "MvE"; break; + case MOVE_WEST: return "MvW"; break; + case MOVE_NORTH: return "MvN"; break; + case MOVE_SOUTH: return "MvS"; break; + case MOVE_X: return "MvX"; break; + case MOVE_Y: return "MvY"; break; + case MOVE_FORWARD: return "Mfd"; break; + case SET_RESPONSIVENESS: return "Res"; break; + case SET_OSCILLATOR_PERIOD: return "OSC"; break; + case EMIT_SIGNAL0: return "SG"; break; + case KILL_FORWARD: return "Klf"; break; + case MOVE_REVERSE: return "Mrv"; break; + case MOVE_LEFT: return "MvL"; break; + case MOVE_RIGHT: return "MvR"; break; + case MOVE_RL: return "MRL"; break; + case MOVE_RANDOM: return "Mrn"; break; + case SET_LONGPROBE_DIST: return "LPD"; break; + default: assert(false); break; + } +} + + +// List the names of the active sensors and actions to stdout. +// "Active" means those sensors and actions that are compiled into +// the code. See sensors-actions.h for how to define the enums. +void printSensorsActions() +{ + unsigned i; + std::cout << "Sensors:" << std::endl; + for (i = 0; i < Sensor::NUM_SENSES; ++i) { + std::cout << " " << sensorName((Sensor)i) << std::endl; + } + std::cout << "Actions:" << std::endl; + for (i = 0; i < Action::NUM_ACTIONS; ++i) { + std::cout << " " << actionName((Action)i) << std::endl; + } + std::cout << std::endl; +} + + +// Format: 32-bit hex strings, one per gene +void Indiv::printGenome() const +{ + constexpr unsigned genesPerLine = 8; + unsigned count = 0; + for (Gene gene : genome) { + if (count == genesPerLine) { + std::cout << std::endl; + count = 0; + } else if (count != 0) { + std::cout << " "; + } + + assert(sizeof(Gene) == 4); + uint32_t n; + std::memcpy(&n, &gene, sizeof(n)); + std::cout << std::hex << std::setfill('0') << std::setw(8) << n; + ++count; + } + std::cout << std::dec << std::endl; +} + + +///* +//Example format: +// +// ACTION_NAMEa from: +// ACTION_NAMEb from: +// SENSOR i +// SENSOR j +// NEURON n +// NEURON m +// Neuron x from: +// SENSOR i +// SENSOR j +// NEURON n +// NEURON m +// Neuron y ... +//*/ +//void Indiv::printNeuralNet() const +//{ +// for (unsigned action = 0; action < Action::NUM_ACTIONS; ++action) { +// bool actionDisplayed = false; +// for (auto & conn : nnet.connections) { +// +// assert((conn.sourceType == NEURON && conn.sourceNum < p.maxNumberNeurons) +// || (conn.sourceType == SENSOR && conn.sourceNum < Sensor::NUM_SENSES)); +// assert((conn.sinkType == ACTION && conn.sinkNum < Action::NUM_ACTIONS) +// || (conn.sinkType == NEURON && conn.sinkNum < p.maxNumberNeurons)); +// +// if (conn.sinkType == ACTION && (conn.sinkNum) == action) { +// if (!actionDisplayed) { +// std::cout << "Action " << actionName((Action)action) << " from:" << std::endl; +// actionDisplayed = true; +// } +// if (conn.sourceType == SENSOR) { +// std::cout << " " << sensorName((Sensor)(conn.sourceNum)); +// } else { +// std::cout << " Neuron " << (conn.sourceNum % nnet.neurons.size()); +// } +// std::cout << " " << conn.weightAsFloat() << std::endl; +// } +// } +// } +// +// for (size_t neuronNum = 0; neuronNum < nnet.neurons.size(); ++neuronNum) { +// bool neuronDisplayed = false; +// for (auto & conn : nnet.connections) { +// if (conn.sinkType == NEURON && (conn.sinkNum) == neuronNum) { +// if (!neuronDisplayed) { +// std::cout << "Neuron " << neuronNum << " from:" << std::endl; +// neuronDisplayed = true; +// } +// if (conn.sourceType == SENSOR) { +// std::cout << " " << sensorName((Sensor)(conn.sourceNum)); +// } else { +// std::cout << " Neuron " << (conn.sourceNum); +// } +// std::cout << " " << conn.weightAsFloat() << std::endl; +// } +// } +// } +//} +// + + +// This prints a neural net in a form that can be processed with +// graph-nnet.py to produce a graphic illustration of the net. +void Indiv::printIGraphEdgeList() const +{ + for (auto & conn : nnet.connections) { + if (conn.sourceType == SENSOR) { + std::cout << sensorShortName((Sensor)(conn.sourceNum)); + } else { + std::cout << "N" << std::to_string(conn.sourceNum); + } + + std::cout << " "; + + if (conn.sinkType == ACTION) { + std::cout << actionShortName((Action)(conn.sinkNum)); + } else { + std::cout << "N" << std::to_string(conn.sinkNum); + } + + std::cout << " " << std::to_string(conn.weight) << std::endl; + } +} + + +float averageGenomeLength() +{ + unsigned count = 100; + unsigned numberSamples = 0; + unsigned long sum = 0; + + while (count-- > 0) { + sum += peeps[randomUint(1, p.population)].genome.size(); + ++numberSamples; + } + return sum / numberSamples; +} + + +// The epoch log contains one line per generation in a format that can be +// fed to graphlog.gp to produce a chart of the simulation progress. +// ToDo: remove hardcoded filename. +void appendEpochLog(unsigned generation, unsigned numberSurvivors, unsigned murderCount) +{ + std::ofstream foutput; + + if (generation == 0) { + foutput.open(p.logDir + "/epoch-log.txt"); + foutput.close(); + } + + foutput.open(p.logDir + "/epoch-log.txt", std::ios::app); + + if (foutput.is_open()) { + foutput << generation << " " << numberSurvivors << " " << geneticDiversity() + << " " << averageGenomeLength() << " " << murderCount << std::endl; + } else { + assert(false); + } +} + + +// Print stats about pheromone usage. +void displaySignalUse() +{ + if (Sensor::SIGNAL0 > Sensor::NUM_SENSES && Sensor::SIGNAL0_FWD > Sensor::NUM_SENSES && Sensor::SIGNAL0_LR > Sensor::NUM_SENSES) { + return; + } + + unsigned long long sum = 0; + unsigned count = 0; + + for (int16_t x = 0; x < p.sizeX; ++x) { + for (int16_t y = 0; y < p.sizeY; ++y) { + unsigned magnitude = signals.getMagnitude(0, { x, y }); + if (magnitude != 0) { + ++count; + sum += magnitude; + } + } + } + std::cout << "Signal spread " << ((double)count / (p.sizeX * p.sizeY)) << "%, average " << ((double)sum / (p.sizeX * p.sizeY)) << std::endl; +} + + +// Print how many connections occur from each kind of sensor neuron and to +// each kind of action neuron over the entire population. This helps us to +// see which sensors and actions are most useful for survival. +void displaySensorActionReferenceCounts() +{ + std::vector sensorCounts(Sensor::NUM_SENSES, 0); + std::vector actionCounts(Action::NUM_ACTIONS, 0); + + for (unsigned index = 1; index <= p.population; ++index) { + if (peeps[index].alive) { + const Indiv &indiv = peeps[index]; + for (const Gene &gene : indiv.nnet.connections) { + if (gene.sourceType == SENSOR) { + assert(gene.sourceNum < Sensor::NUM_SENSES); + ++sensorCounts[(Sensor)gene.sourceNum]; + } + if (gene.sinkType == ACTION) { + assert(gene.sinkNum < Action::NUM_ACTIONS); + ++actionCounts[(Action)gene.sinkNum]; + } + } + } + } + + std::cout << "Sensors in use:" << std::endl; + for (unsigned i = 0; i < sensorCounts.size(); ++i) { + if (sensorCounts[i] > 0) { + std::cout << " " << sensorCounts[i] << " - " << sensorName((Sensor)i) << std::endl; + } + } + std::cout << "Actions in use:" << std::endl; + for (unsigned i = 0; i < actionCounts.size(); ++i) { + if (actionCounts[i] > 0) { + std::cout << " " << actionCounts[i] << " - " << actionName((Action)i) << std::endl; + } + } +} + + +void displaySampleGenomes(unsigned count) +{ + unsigned index = 1; // indexes start at 1 + for (index = 1; count > 0 && index <= p.population; ++index) { + if (peeps[index].alive) { + std::cout << "---------------------------\nIndividual ID " << index << std::endl; + peeps[index].printGenome(); + std::cout << std::endl; + + //peeps[index].printNeuralNet(); + peeps[index].printIGraphEdgeList(); + + + std::cout << "---------------------------" << std::endl; + --count; + } + } + + displaySensorActionReferenceCounts(); +} + +} // end namespace BS diff --git a/src/basicTypes.cpp b/src/basicTypes.cpp new file mode 100644 index 0000000..6ac9fc1 --- /dev/null +++ b/src/basicTypes.cpp @@ -0,0 +1,173 @@ +// basicTypes.cpp + +#include +#include "basicTypes.h" + +namespace BS { + +// This rotates a Dir value by the specified number of steps. There are +// eight steps per full rotation. Positive values are clockwise; negative +// values are counterclockwise. E.g., rotate(4) returns a direction 90 +// degrees to the right. + +Dir Dir::rotate(int n) const +{ + constexpr uint8_t rotateRight[9] { 3, 0, 1, 6, 4, 2, 7, 8, 5 }; + constexpr uint8_t rotateLeft[9] { 1, 2, 5, 0, 4, 8, 3, 6, 7 }; + uint8_t n9 = asInt(); + if (n < 0) { + while (n++ < 0) { + n9 = rotateLeft[n9]; + } + } else if (n > 0) { + while (n-- > 0) { + n9 = rotateRight[n9]; + } + } + return Dir{(Compass)n9}; +} + + +/* + A normalized Coord is a Coord with x and y == -1, 0, or 1. + A normalized Coord may be used as an offset to one of the + 8-neighbors. + + A Dir value maps to a normalized Coord using + + Coord { (d%3) - 1, (trunc)(d/3) - 1 } + + 0 => -1, -1 + 1 => 0, -1 + 2 => 1, -1, + 3 => -1, 0 + 4 => 0, 0 + 5 => 1 0 + 6 => -1, 1 + 7 => 0, 1 + 8 => 1, 1 +*/ +Coord Dir::asNormalizedCoord() const +{ + const int d = asInt(); + return Coord{(int16_t)((d % 3) - 1), (int16_t)((d / 3) - 1)}; +} + + +Polar Dir::asNormalizedPolar() const +{ + return Polar{1, dir9}; +} + + +/* + A normalized Coord has x and y == -1, 0, or 1. + A normalized Coord may be used as an offset to one of the + 8-neighbors. + We'll convert the Coord into a Dir, then convert Dir to normalized Coord. +*/ +Coord Coord::normalize() const +{ + return asDir().asNormalizedCoord(); +} + + +// Calculate the angle, round to the nearest 360/8 degree slice, then +// convert the slice to a Dir8Compass value. +Dir Coord::asDir() const +{ + if (x == 0 && y == 0) { + return Dir{Compass::CENTER}; + } + + const float TWO_PI = 3.1415927f * 2.0f; + float angle = std::atan2((float)y, (float)x); + + if (angle < 0.0f) { + angle = 3.1415927f + (3.1415927f + angle); + } + //assert(angle >= 0.0 && angle <= TWO_PI); + + angle += (TWO_PI / 16.0f); // offset by half a slice + if (angle > TWO_PI) { + angle -= TWO_PI; + } + unsigned slice = (unsigned)(angle / (TWO_PI/8.0f)); // find which division it's in + /* + We have to convert slice values: + + 3 2 1 + 4 0 + 5 6 7 + + into Dir8Compass value: + + 6 7 8 + 3 4 5 + 0 1 2 + */ + const Compass dirconversion[8] { + Compass::E, Compass::NE, Compass::N, Compass::NW, + Compass::W, Compass::SW, Compass::S, Compass::SE }; + return Dir{dirconversion[slice]}; +} + + +Polar Coord::asPolar() const +{ + return Polar{(int)length(), asDir()}; +} + + +/* + Compass values: + + 6 7 8 + 3 4 5 + 0 1 2 +*/ +Coord Polar::asCoord() const +{ + if (dir == Compass::CENTER) { + return Coord{0, 0}; + } + + constexpr double S = (3.14159265359 * 2) / 8; // radians per slice + double compassToRadians[] { 5*S, 6*S, 7*S, 4*S, 0, 0*S, 3*S, 2*S, 1*S }; +// uint8_t asint = dir.asInt(); +// double sliceradians = compassToRadians[asint]; +// double angle = std::cos(sliceradians); +// double angleXmag = mag * angle; +// double adj = angleXmag + 0.5; +// int16_t trunc = (int16_t)adj; + int16_t x = (int16_t)std::round(mag * std::cos(compassToRadians[dir.asInt()])); + int16_t y = (int16_t)std::round(mag * std::sin(compassToRadians[dir.asInt()])); + return Coord{x, y}; +} + + +// returns -1.0 (opposite directions) .. +1.0 (same direction) +// returns 1.0 if either vector is (0,0) +float Coord::raySameness(Coord other) const +{ + float mag1 = std::sqrt(x * x + y * y); + float mag2 = std::sqrt(other.x * other.x + other.y * other.y); + if (mag1 == 0.0 || mag2 == 0.0) { + return 1.0; // anything is "same" as zero vector + } + float dot = x * other.x + y * other.y; + float cos = dot / (mag1 * mag2); + assert(cos >= -1.0001 && cos <= 1.0001); + cos = std::min(std::max(cos, -1.0), 1.0); // clip + return cos; +} + + +// returns -1.0 (opposite directions) .. +1.0 (same direction) +// returns 1.0 if self is (0,0) or d is CENTER +float Coord::raySameness(Dir d) const +{ + return raySameness(d.asNormalizedCoord()); +} + +} // end namespace BS diff --git a/src/basicTypes.h b/src/basicTypes.h new file mode 100644 index 0000000..451099e --- /dev/null +++ b/src/basicTypes.h @@ -0,0 +1,127 @@ +#ifndef BASICTYPES_H_INCLUDED +#define BASICTYPES_H_INCLUDED + +/* +Basic types used throughout the project: + +Compass - an enum with enumerants N=0, NE, E, SW, S, SW, W, NW, CENTER + +Dir, Coord, Polar, and their constructors: + + Dir - abstract type for 8 directions plus center + ctor Dir(Compass = CENTER) + + Coord - signed int16_t pair, absolute location or difference of locations + ctor Coord() = 0,0 + + Polar - signed magnitude and direction + ctor Polar(Coord = 0,0) + +Conversions + + uint8_t = Dir.asInt() + + Dir = Coord.asDir() + Dir = Polar.asDir() + + Coord = Dir.asNormalizedCoord() + Coord = Polar.asCoord() + + Polar = Dir.asNormalizedPolar() + Polar = Coord.asPolar() + +Arithmetic + + Dir.rotate(int n = 0) + + Coord = Coord + Dir + Coord = Coord + Coord + Coord = Coord + Polar + + Polar = Polar + Coord (additive) + Polar = Polar + Polar (additive) + Polar = Polar * Polar (dot product) +*/ + +#include +#include +#include +#include "random.h" + +namespace BS { + +extern bool unitTestBasicTypes(); + +enum class Compass :uint8_t { SW = 0, S, SE, W, CENTER, E, NW, N, NE }; + +struct Dir; +struct Coord; +struct Polar; + +// Supports the eight directions in enum class Compass plus CENTER. +struct __attribute__((packed)) Dir { + static Dir random8() { return Dir(Compass::N).rotate(randomUint(0, 7)); } + + Dir(Compass dir = Compass::CENTER) : dir9{dir} {}; + Dir operator=(Compass d) { dir9 = d; return *this; } + uint8_t asInt() const { return (uint8_t)dir9; } + Coord asNormalizedCoord() const; // (-1, -0, 1, -1, 0, 1) + Polar asNormalizedPolar() const; + + Dir rotate(int n = 0) const; + Dir rotate90DegCW() const { return rotate(2); } + Dir rotate90DegCCW() const { return rotate(-2); } + Dir rotate180Deg() const { return rotate(4); } + + bool operator==(Compass d) const { return asInt() == (uint8_t)d; } + bool operator!=(Compass d) const { return asInt() != (uint8_t)d; } + bool operator==(Dir d) const { return asInt() == d.asInt(); } + bool operator!=(Dir d) const { return asInt() != d.asInt(); } +private: + Compass dir9; +}; + + +// Coordinates range anywhere in the range of int16_t. Coordinate arithmetic +// wraps like int16_t. Can be used, e.g., for a location in the simulator grid, or +// for the difference between two locations. +struct __attribute__((packed)) Coord { + Coord(int16_t x0 = 0, int16_t y0 = 0) : x{x0}, y{y0} { } + bool isNormalized() const { return x >= -1 && x <= 1 && y >= -1 && y <= 1; } + Coord normalize() const; + unsigned length() const { return (int)(std::sqrt(x * x + y * y)); } // round down + Dir asDir() const; + Polar asPolar() const; + + bool operator==(Coord c) const { return x == c.x && y == c.y; } + bool operator!=(Coord c) const { return x != c.x || y != c.y; } + Coord operator+(Coord c) const { return Coord{(int16_t)(x + c.x), (int16_t)(y + c.y)}; } + Coord operator-(Coord c) const { return Coord{(int16_t)(x - c.x), (int16_t)(y - c.y)}; } + Coord operator*(int a) const { return Coord{(int16_t)(x * a), (int16_t)(y * a)}; } + Coord operator+(Dir d) const { return *this + d.asNormalizedCoord(); } + Coord operator-(Dir d) const { return *this - d.asNormalizedCoord(); } + + float raySameness(Coord other) const; // returns -1.0 (opposite) .. 1.0 (same) + float raySameness(Dir d) const; // returns -1.0 (opposite) .. 1.0 (same) +public: + int16_t x; + int16_t y; +}; + + +// Polar magnitudes are signed 32-bit integers so that they can extend across any 2D +// area defined by the Coord class. +struct __attribute__((packed)) Polar { + explicit Polar(int mag0 = 0, Compass dir0 = Compass::CENTER) + : mag{mag0}, dir{Dir{dir0}} { } + explicit Polar(int mag0, Dir dir0) + : mag{mag0}, dir{dir0} { } + Coord asCoord() const; +public: + int mag; + Dir dir; +}; + +} // end namespace BS + +#endif // BASICTYPES_H_INCLUDED diff --git a/src/biosim4.layout b/src/biosim4.layout new file mode 100644 index 0000000..061489b --- /dev/null +++ b/src/biosim4.layout @@ -0,0 +1,10 @@ + + + + + + + + + + diff --git a/src/createBarrier.cpp b/src/createBarrier.cpp new file mode 100644 index 0000000..10aeac6 --- /dev/null +++ b/src/createBarrier.cpp @@ -0,0 +1,184 @@ +// createBarrier.cpp + +#include +#include "simulator.h" + +namespace BS { + +// This generates barrier points, which are grid locations with value +// BARRIER. A list of barrier locations is saved in private member +// Grid::barrierLocations and, for some scenarios, Grid::barrierCenters. +// Those members are available read-only with Grid::getBarrierLocations(). +// This function assumes an empty grid. This is typically called by +// the main simulator thread after Grid::init() or Grid::zeroFill(). + +// This file typically is under constant development and change for +// specific scenarios. + +void Grid::createBarrier(unsigned barrierType) +{ + barrierLocations.clear(); + barrierCenters.clear(); // used only for some barrier types + + auto drawBox = [&](int16_t minX, int16_t minY, int16_t maxX, int16_t maxY) { + for (int16_t x = minX; x <= maxX; ++x) { + for (int16_t y = minY; y <= maxY; ++y) { + grid.set(x, y, BARRIER); + barrierLocations.push_back( {x, y} ); + } + } + }; + + switch(barrierType) { + case 0: + return; + + // Vertical bar in constant location + case 1: + { + int16_t minX = p.sizeX / 2; + int16_t maxX = minX + 1; + int16_t minY = p.sizeY / 4; + int16_t maxY = minY + p.sizeY / 2; + + for (int16_t x = minX; x <= maxX; ++x) { + for (int16_t y = minY; y <= maxY; ++y) { + grid.set(x, y, BARRIER); + barrierLocations.push_back( {x, y} ); + } + } + } + break; + + // Vertical bar in random location + case 2: + { + int16_t minX = randomUint(20, p.sizeX - 20); + int16_t maxX = minX + 1; + int16_t minY = randomUint(20, p.sizeY / 2 - 20); + int16_t maxY = minY + p.sizeY / 2; + + for (int16_t x = minX; x <= maxX; ++x) { + for (int16_t y = minY; y <= maxY; ++y) { + grid.set(x, y, BARRIER); + barrierLocations.push_back( {x, y} ); + } + } + } + break; + + // five blocks staggered + case 3: + { + int16_t blockSizeX = 2; + int16_t blockSizeY = p.sizeX / 3; + + int16_t x0 = p.sizeX / 4 - blockSizeX / 2; + int16_t y0 = p.sizeY / 4 - blockSizeY / 2; + int16_t x1 = x0 + blockSizeX; + int16_t y1 = y0 + blockSizeY; + + drawBox(x0, y0, x1, y1); + x0 += p.sizeX / 2; + x1 = x0 + blockSizeX; + drawBox(x0, y0, x1, y1); + y0 += p.sizeY / 2; + y1 = y0 + blockSizeY; + drawBox(x0, y0, x1, y1); + x0 -= p.sizeX / 2; + x1 = x0 + blockSizeX; + drawBox(x0, y0, x1, y1); + x0 = p.sizeX / 2 - blockSizeX / 2; + x1 = x0 + blockSizeX; + y0 = p.sizeY / 2 - blockSizeY / 2; + y1 = y0 + blockSizeY; + drawBox(x0, y0, x1, y1); + return; + } + break; + + // Horizontal bar in constant location + case 4: + { + int16_t minX = p.sizeX / 4; + int16_t maxX = minX + p.sizeX / 2; + int16_t minY = p.sizeY / 2 + p.sizeY / 4; + int16_t maxY = minY + 2; + + for (int16_t x = minX; x <= maxX; ++x) { + for (int16_t y = minY; y <= maxY; ++y) { + grid.set(x, y, BARRIER); + barrierLocations.push_back( {x, y} ); + } + } + } + break; + + // Three floating islands -- different locations every generation + case 5: + { + float radius = 3.0; + unsigned margin = 2 * (int)radius; + + auto randomLoc = [&]() { +// return Coord( (int16_t)randomUint((int)radius + margin, p.sizeX - ((float)radius + margin)), +// (int16_t)randomUint((int)radius + margin, p.sizeY - ((float)radius + margin)) ); + return Coord( (int16_t)randomUint(margin, p.sizeX - margin), + (int16_t)randomUint(margin, p.sizeY - margin) ); + }; + + Coord center0 = randomLoc(); + Coord center1; + Coord center2; + + do { + center1 = randomLoc(); + } while ((center0 - center1).length() < margin); + + do { + center2 = randomLoc(); + } while ((center0 - center2).length() < margin || (center1 - center2).length() < margin); + + barrierCenters.push_back(center0); + //barrierCenters.push_back(center1); + //barrierCenters.push_back(center2); + + auto f = [&](Coord loc) { + grid.set(loc, BARRIER); + barrierLocations.push_back(loc); + }; + + visitNeighborhood(center0, radius, f); + //visitNeighborhood(center1, radius, f); + //visitNeighborhood(center2, radius, f); + } + break; + + // Spots, specified number, radius, locations + case 6: + { + unsigned numberOfLocations = 5; + float radius = 5.0; + + auto f = [&](Coord loc) { + grid.set(loc, BARRIER); + barrierLocations.push_back(loc); + }; + + unsigned verticalSliceSize = p.sizeY / (numberOfLocations + 1); + + for (unsigned n = 1; n <= numberOfLocations; ++n) { + Coord loc = { (int16_t)(p.sizeX / 2), + (int16_t)(n * verticalSliceSize) }; + visitNeighborhood(loc, radius, f); + barrierCenters.push_back(loc); + } + } + break; + + default: + assert(false); + } +} + +} // end namespace BS diff --git a/src/endOfGeneration.cpp b/src/endOfGeneration.cpp new file mode 100644 index 0000000..6e3afa2 --- /dev/null +++ b/src/endOfGeneration.cpp @@ -0,0 +1,36 @@ +// endOfGeneration.cpp + +#include +#include +#include +#include +#include +#include "simulator.h" +#include "imageWriter.h" + +namespace BS { + +// At the end of each generation, we save a video file (if p.saveVideo is true) and +// print some genomic statistics to stdout (if p.updateGraphLog is true). + +void endOfGeneration(unsigned generation) +{ + { + if (p.saveVideo && + ((generation % p.videoStride) == 0 + || generation <= p.videoSaveFirstFrames + || (generation >= p.replaceBarrierTypeGenerationNumber + && generation <= p.replaceBarrierTypeGenerationNumber + p.videoSaveFirstFrames))) { + imageWriter.saveGenerationVideo(generation); + } + } + + { + if (p.updateGraphLog && (generation == 1 || ((generation % p.updateGraphLogStride) == 0))) { +#pragma GCC diagnostic ignored "-Wunused-result" + std::system(p.graphLogUpdateCommand.c_str()); + } + } +} + +} // end namespace BS diff --git a/src/endOfSimStep.cpp b/src/endOfSimStep.cpp new file mode 100644 index 0000000..e888c60 --- /dev/null +++ b/src/endOfSimStep.cpp @@ -0,0 +1,92 @@ +// endOfSimStep.cpp + +#include +#include +#include "simulator.h" +#include "imageWriter.h" + +namespace BS { + +/* +At the end of each sim step, this function is called in single-thread +mode to take care of several things: + +1. We may kill off some agents if a "radioactive" scenario is in progress. +2. We may flag some agents as meeting some challenge criteria, if such + a scenario is in progress. +3. We then drain the deferred death queue. +4. We then drain the deferred movement queue. +5. We fade the signal layer(s) (pheromones). +6. We save the resulting world condition as a single image frame (if + p.saveVideo is true). +*/ + +void endOfSimStep(unsigned simStep, unsigned generation) +{ + if (p.challenge == CHALLENGE_RADIOACTIVE_WALLS) { + // During the first half of the generation, the west wall is radioactive, + // where X == 0. In the last half of the generation, the east wall is + // radioactive, where X = the area width - 1. There's an exponential + // falloff of the danger, falling off to zero at the arena half line. + int16_t radioactiveX = (simStep < p.stepsPerGeneration / 2) ? 0 : p.sizeX - 1; + + for (uint16_t index = 1; index <= p.population; ++index) { // index 0 is reserved + Indiv &indiv = peeps[index]; + int16_t distanceFromRadioactiveWall = std::abs(indiv.loc.x - radioactiveX); + if (distanceFromRadioactiveWall < p.sizeX / 2) { + float chanceOfDeath = 1.0 / distanceFromRadioactiveWall; + if (randomUint() / (float)RANDOM_UINT_MAX < chanceOfDeath) { + peeps.queueForDeath(indiv); + } + } + } + } + + // If the individual is touching any wall, we set its challengeFlag to true. + // At the end of the generation, all those with the flag true will reproduce. + if (p.challenge == CHALLENGE_TOUCH_ANY_WALL) { + for (uint16_t index = 1; index <= p.population; ++index) { // index 0 is reserved + Indiv &indiv = peeps[index]; + if (indiv.loc.x == 0 || indiv.loc.x == p.sizeX - 1 + || indiv.loc.y == 0 || indiv.loc.y == p.sizeY - 1) { + indiv.challengeBits = true; + } + } + } + + // If this challenge is enabled, the individual gets a bit set in their challengeBits + // member if they are within a specified radius of a barrier center. They have to + // visit the barriers in sequential order. + if (p.challenge == CHALLENGE_LOCATION_SEQUENCE) { + float radius = 9.0; + for (uint16_t index = 1; index <= p.population; ++index) { // index 0 is reserved + Indiv &indiv = peeps[index]; + for (unsigned n = 0; n < grid.getBarrierCenters().size(); ++n) { + unsigned bit = 1 << n; + if ((indiv.challengeBits & bit) == 0) { + if ((indiv.loc - grid.getBarrierCenters()[n]).length() <= radius) { + indiv.challengeBits |= bit; + } + break; + } + } + } + } + + peeps.drainDeathQueue(); + peeps.drainMoveQueue(); + signals.fade(0); // takes layerNum todo!!! + + // saveVideoFrameSync() is the synchronous version of saveVideFrame() + if (p.saveVideo && + ((generation % p.videoStride) == 0 + || generation <= p.videoSaveFirstFrames + || (generation >= p.replaceBarrierTypeGenerationNumber + && generation <= p.replaceBarrierTypeGenerationNumber + p.videoSaveFirstFrames))) { + if (!imageWriter.saveVideoFrameSync(simStep, generation)) { + std::cout << "imageWriter busy" << std::endl; + } + } +} + +} // end namespace BS diff --git a/src/executeActions.cpp b/src/executeActions.cpp new file mode 100644 index 0000000..c021b9b --- /dev/null +++ b/src/executeActions.cpp @@ -0,0 +1,245 @@ +// executeActions.cpp -- Executes the actions computed for a single simStep for a specified individual + +#include +#include +#include +#include +#include +#include "simulator.h" +#include "omp.h" + +namespace BS { + +// Given a factor in the range 0.0..1.0, return a bool with the +// probability of it being true proportional to factor. For example, if +// factor == 0.2, then there is a 20% chance this function will +// return true. +bool prob2bool(float factor) +{ + assert(factor >= 0.0 && factor <= 1.0); + return (randomUint() / (float)RANDOM_UINT_MAX) < factor; +} + + +// This takes a probability from 0.0..1.0 and adjusts it according to an +// exponential curve. The steepness of the curve is determined by the K factor +// which is a small positive integer. This tends to reduce the activity level +// a bit (makes the peeps less reactive and jittery). +float responseCurve(float r) +{ + const float k = p.responsivenessCurveKFactor; + return std::pow((r - 2.0), -2.0 * k) - std::pow(2.0, -2.0 * k) * (1.0 - r); +} + + +/********************************************************************************** +Action levels are driven by sensors or internal neurons as connected by an agent's +neural net brain. Each agent's neural net is reevaluated once each simulator +step (simStep). After evaluating the action neuron outputs, this function is +called to execute the actions according to their output levels. This function is +called in multi-threaded mode and operates on a single individual while other +threads are doing to the same to other individuals. + +Action (their output) values arrive here as floating point values of arbitrary +range (because they are the raw sums of zero or more weighted inputs) and will +eventually be converted in this function to a probability 0.0..1.0 of actually +getting executed. + +For the various possible action neurons, if they are driven by a sufficiently +strong level, we do this: + + MOVE_* actions- queue our agent for deferred movement with peeps.queueForMove(); the + queue will be executed at the end of the multithreaded loop in a single thread. + SET_RESPONSIVENESS action - immediately change indiv.responsiveness to the action + level scaled to 0.0..1.0 (because we have exclusive access to this member in + our own individual during this function) + SET_OSCILLATOR_PERIOD action - immediately change our individual's indiv.oscPeriod + to the action level exponentially scaled to 2..2048 (TBD) + EMIT_SIGNALn action(s) - immediately increment the signal level at our agent's + location using signals.increment() (using a thread-safe call) + KILL_FORWARD action - queue the other agent for deferred death with + peeps.queueForDeath() + +The deferred movement and death queues will be emptied by the caller at the end of the +simulator step by endOfSimStep() in a single thread after all individuals have been +evaluated multithreadedly. +**********************************************************************************/ + +void executeActions(Indiv &indiv, std::array &actionLevels) +{ + // Only a subset of all possible actions might be enabled (i.e., compiled in). + // This returns true if the specified action is enabled. See sensors-actions.h + // for how to enable sensors and actions during compilation. + auto isEnabled = [](enum Action action){ return (int)action < (int)Action::NUM_ACTIONS; }; + + // Responsiveness action - convert neuron action level from arbitrary float range + // to the range 0.0..1.0. If this action neuron is enabled but not driven, will + // default to mid-level 0.5. + if (isEnabled(Action::SET_RESPONSIVENESS)) { + float level = actionLevels[Action::SET_RESPONSIVENESS]; // default 0.0 + level = (std::tanh(level) + 1.0) / 2.0; // convert to 0.0..1.0 + indiv.responsiveness = level; + } + + // For the rest of the action outputs, we'll apply an adjusted responsiveness + // factor (see responseCurve() for more info). Range 0.0..1.0. + float responsivenessAdjusted = responseCurve(indiv.responsiveness); + + // Oscillator period action - convert action level nonlinearly to + // 2..4*p.stepsPerGeneration. If this action neuron is enabled but not driven, + // will default to 1.5 + e^(3.5) = a period of 34 simSteps. + if (isEnabled(Action::SET_OSCILLATOR_PERIOD)) { + auto periodf = actionLevels[Action::SET_OSCILLATOR_PERIOD]; + float newPeriodf01 = (std::tanh(periodf) + 1.0) / 2.0; // convert to 0.0..1.0 + unsigned newPeriod = 1 + (int)(1.5 + std::exp(7.0 * newPeriodf01)); + assert(newPeriod >= 2 && newPeriod <= 2048); + indiv.oscPeriod = newPeriod; + } + + // Set longProbeDistance - convert action level to 1..maxLongProbeDistance. + // If this action neuron is enabled but not driven, will default to + // mid-level period of 17 simSteps. + if (isEnabled(Action::SET_LONGPROBE_DIST)) { + constexpr unsigned maxLongProbeDistance = 32; + float level = actionLevels[SET_LONGPROBE_DIST]; + level = (std::tanh(level) + 1.0) / 2.0; // convert to 0.0..1.0 + level = 1 + level * maxLongProbeDistance; + indiv.longProbeDist = (unsigned)level; + } + + // Emit signal0 - if this action value is below a threshold, nothing emitted. + // Otherwise convert the action value to a probability of emitting one unit of + // signal (pheromone). + // Pheromones may be emitted immediately (see signals.cpp). If this action neuron + // is enabled but not driven, nothing will be emitted. + if (isEnabled(Action::EMIT_SIGNAL0)) { + constexpr float emitThreshold = 0.5; // 0.0..1.0; 0.5 is midlevel + float level = actionLevels[Action::EMIT_SIGNAL0]; + level = (std::tanh(level) + 1.0) / 2.0; // convert to 0.0..1.0 + level *= responsivenessAdjusted; + if (level > emitThreshold && prob2bool(level)) { + signals.increment(0, indiv.loc); + } + } + + // Kill forward -- if this action value is > threshold, value is converted to probability + // of an attempted murder. Probabilities under the threshold are considered 0.0. + // If this action neuron is enabled but not driven, the neighbors are safe. + if (isEnabled(Action::KILL_FORWARD) && p.killEnable) { + constexpr float killThreshold = 0.5; // 0.0..1.0; 0.5 is midlevel + float level = actionLevels[Action::KILL_FORWARD]; + level = (std::tanh(level) + 1.0) / 2.0; // convert to 0.0..1.0 + level *= responsivenessAdjusted; + if (level > killThreshold && prob2bool((level - ACTION_MIN) / ACTION_RANGE)) { + Coord otherLoc = indiv.loc + indiv.lastMoveDir; + if (grid.isInBounds(otherLoc) && grid.isOccupiedAt(otherLoc)) { + Indiv &indiv2 = peeps.getIndiv(otherLoc); + if (indiv2.alive) { + assert((indiv.loc - indiv2.loc).length() == 1); + peeps.queueForDeath(indiv2); + } + } + } + } + + // ------------- Movement action neurons --------------- + + // There are multiple action neurons for movement. Each type of movement neuron + // urges the individual to move in some specific direction. We sum up all the + // X and Y components of all the movement urges, then pass the X and Y sums through + // a transfer function (tanh()) to get a range -1.0..1.0. The absolute values of the + // X and Y values are passed through prob2bool() to convert to -1, 0, or 1, then + // multiplied by the component's signum. This results in the x and y components of + // a normalized movement offset. I.e., the probability of movement in either + // dimension is the absolute value of tanh of the action level X,Y components and + // the direction is the sign of the X, Y components. For example, for a particular + // action neuron: + // X, Y == -5.9, +0.3 as raw action levels received here + // X, Y == -0.999, +0.29 after passing raw values through tanh() + // Xprob, Yprob == 99.9%, 29% probability of X and Y becoming 1 (or -1) + // X, Y == -1, 0 after applying the sign and probability + // The agent will then be moved West (an offset of -1, 0) if it's a legal move. + + float level; + Coord offset; + Coord lastMoveOffset = indiv.lastMoveDir.asNormalizedCoord(); + + // moveX,moveY will be the accumulators that will hold the sum of all the + // urges to move along each axis. (+- floating values of arbitrary range) + float moveX = isEnabled(Action::MOVE_X) ? actionLevels[Action::MOVE_X] : 0.0; + float moveY = isEnabled(Action::MOVE_Y) ? actionLevels[Action::MOVE_Y] : 0.0; + + if (isEnabled(Action::MOVE_EAST)) moveX += actionLevels[Action::MOVE_EAST]; + if (isEnabled(Action::MOVE_WEST)) moveX -= actionLevels[Action::MOVE_WEST]; + if (isEnabled(Action::MOVE_NORTH)) moveY += actionLevels[Action::MOVE_NORTH]; + if (isEnabled(Action::MOVE_SOUTH)) moveY -= actionLevels[Action::MOVE_SOUTH]; + + if (isEnabled(Action::MOVE_FORWARD)) { + level = actionLevels[Action::MOVE_FORWARD]; + moveX += lastMoveOffset.x * level; + moveY += lastMoveOffset.y * level; + } + if (isEnabled(Action::MOVE_REVERSE)) { + level = actionLevels[Action::MOVE_REVERSE]; + moveX -= lastMoveOffset.x * level; + moveY -= lastMoveOffset.y * level; + } + if (isEnabled(Action::MOVE_LEFT)) { + level = actionLevels[Action::MOVE_LEFT]; + offset = indiv.lastMoveDir.rotate90DegCCW().asNormalizedCoord(); + moveX += offset.x * level; + moveY += offset.y * level; + } + if (isEnabled(Action::MOVE_RIGHT)) { + level = actionLevels[Action::MOVE_RIGHT]; + offset = indiv.lastMoveDir.rotate90DegCW().asNormalizedCoord(); + moveX += offset.x * level; + moveY += offset.y * level; + } + if (isEnabled(Action::MOVE_RL)) { + level = actionLevels[Action::MOVE_RL]; + if (level < 0.0) { + offset = indiv.lastMoveDir.rotate90DegCCW().asNormalizedCoord(); + } else if (level > 0.0) { + offset = indiv.lastMoveDir.rotate90DegCW().asNormalizedCoord(); + } else { + offset = { 0, 0 }; + } + moveX += offset.x * level; + moveY += offset.y * level; + } + + if (isEnabled(Action::MOVE_RANDOM)) { + level = actionLevels[Action::MOVE_RANDOM]; + offset = Dir::random8().asNormalizedCoord(); + moveX += offset.x * level; + moveY += offset.y * level; + } + + // Convert the accumulated X, Y sums to the range -1.0..1.0 and scale by the + // individual's responsiveness (0.0..1.0) (adjusted by a curve) + moveX = std::tanh(moveX); + moveY = std::tanh(moveY); + moveX *= responsivenessAdjusted; + moveY *= responsivenessAdjusted; + + // The probability of movement along each axis is the absolute value + int16_t probX = (int16_t)prob2bool(std::abs(moveX)); // convert abs(level) to 0 or 1 + int16_t probY = (int16_t)prob2bool(std::abs(moveY)); // convert abs(level) to 0 or 1 + + // The direction of movement (if any) along each axis is the sign + int16_t signumX = moveX < 0.0 ? -1 : 1; + int16_t signumY = moveY < 0.0 ? -1 : 1; + + // Generate a normalized movement offset, where each component is -1, 0, or 1 + Coord movementOffset = { (int16_t)(probX * signumX), (int16_t)(probY * signumY) }; + + // Move there if it's a valid location + Coord newLoc = indiv.loc + movementOffset; + if (grid.isInBounds(newLoc) && grid.isEmptyAt(newLoc)) { + peeps.queueForMove(indiv, newLoc); + } +} + +} // end namespace BS + diff --git a/src/feedForward.cpp b/src/feedForward.cpp new file mode 100644 index 0000000..2ad53dc --- /dev/null +++ b/src/feedForward.cpp @@ -0,0 +1,98 @@ +// feedForward.cpp - reads sensors, returns actions, for an individual + +#include +#include +#include +#include "simulator.h" + +namespace BS { + +/******************************************************************************** +This function does a neural net feed-forward operation, from sensor (input) neurons +through internal neurons to action (output) neurons. The feed-forward +calculations are evaluated once each simulator step (simStep). + +There is no back-propagation in this simulator. Once an individual's neural net +brain is wired at birth, the weights and topology do not change during the +individual's lifetime. + +The data structure Indiv::neurons contains internal neurons, and Indiv::connections +holds the connections between the neurons. + +We have three types of neurons: + + input sensors - each gives a value in the range SENSOR_MIN.. SENSOR_MAX (0.0..1.0). + Values are obtained from getSensor(). + + internal neurons - each takes inputs from sensors or other internal neurons; + each has output value in the range NEURON_MIN..NEURON_MAX (-1.0..1.0). The + output value for each neuron is stored in Indiv::neurons[] and survives from + one simStep to the next. (For example, a neuron that feeds itself will use + its output value that was latched from the previous simStep.) Inputs to the + neurons are summed each simStep in a temporary container and then discarded + after the neurons' outputs are computed. + + action (output) neurons - each takes inputs from sensors or other internal + neurons; In this function, each has an output value in an arbitrary range + (because they are the raw sums of zero or more weighted inputs). + The values of the action neurons are saved in local container + actionLevels[] which is returned to the caller by value (thanks RVO). +********************************************************************************/ + +std::array Indiv::feedForward(unsigned simStep) +{ + // This container is used to return values for all the action outputs. This array + // contains one value per action neuron, which is the sum of all its weighted + // input connections. The sum has an arbitrary range. Return by value assumes compiler + // return value optimization. + std::array actionLevels; + actionLevels.fill(0.0); // undriven actions default to value 0.0 + + // Weighted inputs to each neuron are summed in neuronAccumulators[] + std::vector neuronAccumulators(nnet.neurons.size(), 0.0); + + // Connections were ordered at birth so that all connections to neurons get + // processed here before any connections to actions. As soon as we encounter the + // first connection to an action, we'll pass all the neuron input accumulators + // through a transfer function and update the neuron outputs in the indiv, + // except for undriven neurons which act as bias feeds and don't change. The + // transfer function will leave each neuron's output in the range -1.0..1.0. + + bool neuronOutputsComputed = false; + for (Gene & conn : nnet.connections) { + if (conn.sinkType == ACTION && !neuronOutputsComputed) { + // We've handled all the connections from sensors and now we are about to + // start on the connections to the action outputs, so now it's time to + // update and latch all the neuron outputs to their proper range (-1.0..1.0) + for (unsigned neuronIndex = 0; neuronIndex < nnet.neurons.size(); ++neuronIndex) { + if (nnet.neurons[neuronIndex].driven) { + nnet.neurons[neuronIndex].output = std::tanh(neuronAccumulators[neuronIndex]); + } + } + neuronOutputsComputed = true; + } + + // Obtain the connection's input value from a sensor neuron or other neuron + // The values are summed for now, later passed through a transfer function + float inputVal; + if (conn.sourceType == SENSOR) { + inputVal = getSensor((Sensor)conn.sourceNum, simStep); + } else { + inputVal = nnet.neurons[conn.sourceNum].output; + } + + // Weight the connection's value and add to neuron accumulator or action accumulator. + // The action and neuron accumulators will therefore contain +- float values in + // an arbitrary range. + if (conn.sinkType == ACTION) { + actionLevels[conn.sinkNum] += inputVal * conn.weightAsFloat(); + } else { + neuronAccumulators[conn.sinkNum] += inputVal * conn.weightAsFloat(); + } + } + + return actionLevels; +} + +} // end namespace BS + diff --git a/src/genome-compare.cpp b/src/genome-compare.cpp new file mode 100644 index 0000000..da2b93c --- /dev/null +++ b/src/genome-compare.cpp @@ -0,0 +1,175 @@ +// genome-compare.cpp -- compute similarity of two genomes + +#include +#include "simulator.h" + +namespace BS { + +// Approximate gene match: Has to match same source, sink, with similar weight +// +bool genesMatch(const Gene &g1, const Gene &g2) +{ + return g1.sinkNum == g2.sinkNum + && g1.sourceNum == g2.sourceNum + && g1.sinkType == g2.sinkType + && g1.sourceType == g2.sourceType + && g1.weight == g2.weight; +} + + +// The jaro_winkler_distance() function is adapted from the C version at +// https://github.com/miguelvps/c/blob/master/jarowinkler.c +// under a GNU license, ver. 3. This comparison function is useful if +// the simulator allows genomes to change length, or if genes are allowed +// to relocate to different offsets in the genome. I.e., this function is +// tolerant of gaps, relocations, and genomes of unequal lengths. +// +float jaro_winkler_distance(const Genome &genome1, const Genome &genome2) { + float dw; + auto max = [](int a, int b) { return a > b ? a : b; }; + auto min = [](int a, int b) { return a < b ? a : b; }; + + const auto &s = genome1; + const auto &a = genome2; + + int i, j, l; + int m = 0, t = 0; + int sl = s.size(); // strlen(s); + int al = a.size(); // strlen(a); + + constexpr unsigned maxNumGenesToCompare = 20; + sl = min(maxNumGenesToCompare, sl); // optimization: approximate for long genomes + al = min(maxNumGenesToCompare, al); + + std::vector sflags(sl, 0); + std::vector aflags(al, 0); + int range = max(0, max(sl, al) / 2 - 1); + + if (!sl || !al) + return 0.0; + + /* calculate matching characters */ + for (i = 0; i < al; i++) { + for (j = max(i - range, 0), l = min(i + range + 1, sl); j < l; j++) { + if (genesMatch(a[i], s[j]) && !sflags[j]) { + sflags[j] = 1; + aflags[i] = 1; + m++; + break; + } + } + } + + if (!m) + return 0.0; + + /* calculate character transpositions */ + l = 0; + for (i = 0; i < al; i++) { + if (aflags[i] == 1) { + for (j = l; j < sl; j++) { + if (sflags[j] == 1) { + l = j + 1; + break; + } + } + if (!genesMatch(a[i], s[j])) + t++; + } + } + t /= 2; + + /* Jaro distance */ + dw = (((float)m / sl) + ((float)m / al) + ((float)(m - t) / m)) / 3.0f; + return dw; +} + + +// Works only for genomes of equal length +float hammingDistanceBits(const Genome &genome1, const Genome &genome2) +{ + assert(genome1.size() == genome2.size()); + + const unsigned int *p1 = (const unsigned int *)genome1.data(); + const unsigned int *p2 = (const unsigned int *)genome2.data(); + const unsigned numElements = genome1.size(); + const unsigned bytesPerElement = sizeof(genome1[0]); + const unsigned lengthBytes = numElements * bytesPerElement; + const unsigned lengthBits = lengthBytes * 8; + unsigned bitCount = 0; + + for (unsigned index = 0; index < genome1.size(); ++p1, ++p2, ++index) { + bitCount += __builtin_popcount(*p1 ^ *p2); + } + + // For two completely random bit patterns, about half the bits will differ, + // resulting in c. 50% match. We will scale that by 2X to make the range + // from 0 to 1.0. We clip the value to 1.0 in case the two patterns are + // negatively correlated for some reason. + return 1.0 - std::min(1.0, (2.0 * bitCount) / (float)lengthBits); +} + + +// Works only for genomes of equal length +float hammingDistanceBytes(const Genome &genome1, const Genome &genome2) +{ + assert(genome1.size() == genome2.size()); + + const unsigned int *p1 = (const unsigned int *)genome1.data(); + const unsigned int *p2 = (const unsigned int *)genome2.data(); + const unsigned numElements = genome1.size(); + const unsigned bytesPerElement = sizeof(genome1[0]); + const unsigned lengthBytes = numElements * bytesPerElement; + unsigned byteCount = 0; + + for (unsigned index = 0; index < genome1.size(); ++p1, ++p2, ++index) { + byteCount += (unsigned)(*p1 == *p2); + } + + return byteCount / (float)lengthBytes; +} + + +// Returns 0.0..1.0 +// +// ToDo: optimize by approximation for long genomes +float genomeSimilarity(const Genome &g1, const Genome &g2) +{ + switch (p.genomeComparisonMethod) { + case 0: + return jaro_winkler_distance(g1, g2); + case 1: + return hammingDistanceBits(g1, g2); + case 2: + return hammingDistanceBytes(g1, g2); + default: + assert(false); + } +} + + +// returns 0.0..1.0 +// Samples random pairs of individuals regardless if they are alive or not +float geneticDiversity() +{ + if (p.population < 2) { + return 0.0; + } + + // count limits the number of genomes sampled for performance reasons. + unsigned count = std::min(1000U, p.population); // todo: !!! p.analysisSampleSize; + int numSamples = 0; + float similaritySum = 0.0f; + + while (count > 0) { + unsigned index0 = randomUint(1, p.population - 1); // skip first and last elements + unsigned index1 = index0 + 1; + similaritySum += genomeSimilarity(peeps[index0].genome, peeps[index1].genome); + --count; + ++numSamples; + } + float diversity = 1.0f - (similaritySum / numSamples); + return diversity; +} + +} // end namespace BS diff --git a/src/genome-neurons.h b/src/genome-neurons.h new file mode 100644 index 0000000..cdbf5c1 --- /dev/null +++ b/src/genome-neurons.h @@ -0,0 +1,93 @@ +#ifndef GENOME_H_INCLUDED +#define GENOME_H_INCLUDED + +#include +#include +#include +#include "sensors-actions.h" +#include "random.h" + +namespace BS { + +// Each gene specifies one synaptic connection in a neural net. Each +// connection has an input (source) which is either a sensor or another neuron. +// Each connection has an output, which is either an action or another neuron. +// Each connection has a floating point weight derived from a signed 16-bit +// value. The signed integer weight is scaled to a small range, then cubed +// to provide fine resolution near zero. + +constexpr uint8_t SENSOR = 1; // always a source +constexpr uint8_t ACTION = 1; // always a sink +constexpr uint8_t NEURON = 0; // can be either a source or sink + +struct Gene { //__attribute__((packed)) Gene { + uint16_t sourceType:1; // SENSOR or NEURON + uint16_t sourceNum:7; + uint16_t sinkType:1; // NEURON or ACTION + uint16_t sinkNum:7; + int16_t weight; + + static constexpr float f1 = 8.0; + static constexpr float f2 = 64.0; + //float weightAsFloat() { return std::pow(weight / f1, 3.0) / f2; } + float weightAsFloat() const { return weight / 8192.0; } + static int16_t makeRandomWeight() { return randomUint(0, 0xefff) - 0x8000; } +}; + + +// An individual's genome is a set of Genes (see Gene comments above). Each +// gene is equivalent to one connection in a neural net. An individual's +// neural net is derived from its set of genes. +typedef std::vector Genome; + + +// An individual's "brain" is a neural net specified by a set +// of Genes where each Gene specifies one connection in the neural net (see +// Genome comments above). Each neuron has a single output which is +// connected to a set of sinks where each sink is either an action output +// or another neuron. Each neuron has a set of input sources where each +// source is either a sensor or another neuron. There is no concept of +// layers in the net: it's a free-for-all topology with forward, backwards, +// and sideways connection allowed. Weighted connections are allowed +// directly from any source to any action. + +// Currently the genome does not specify the activation function used in +// the neurons. (May be hardcoded to std::tanh() !!!) + +// When the input is a sensor, the input value to the sink is the raw +// sensor value of type float and depends on the sensor. If the output +// is an action, the source's output value is interpreted by the action +// node and whether the action occurs or not depends on the action's +// implementation. + +// In the genome, neurons are identified by 15-bit unsigned indices, +// which are reinterpreted as values in the range 0..p.genomeMaxLength-1 +// by taking the 15-bit index modulo the max number of allowed neurons. +// In the neural net, the neurons that end up connected get new indices +// assigned sequentially starting at 0. + + +struct NeuralNet { + std::vector connections; // connections are equivalent to genes + + struct Neuron { + float output; + bool driven; // undriven neurons have fixed output values + }; + std::vector neurons; +}; + +// When a new population is generated and every individual is given a +// neural net, the neuron outputs must be initialized to something: +//constexpr float initialNeuronOutput() { return (NEURON_RANGE / 2.0) + NEURON_MIN; } +constexpr float initialNeuronOutput() { return 0.5; } + +extern Gene makeRandomGene(); +extern Genome makeRandomGenome(); +extern void unitTestConnectNeuralNetWiringFromGenome(); +extern float genomeSimilarity(const Genome &g1, const Genome &g2); // 0.0..1.0 +extern float geneticDiversity(); // 0.0..1.0 + +} // end namespace BS + +#endif // GENOME_H_INCLUDED diff --git a/src/genome.cpp b/src/genome.cpp new file mode 100644 index 0000000..2911186 --- /dev/null +++ b/src/genome.cpp @@ -0,0 +1,428 @@ +// genome.cpp + +#include +#include +#include +#include +#include +#include +#include "simulator.h" +#include "random.h" + +namespace BS { + +// This structure is used while converting the connection list to a +// neural net. This helps us to find neurons that don't feed anything +// so that they can be removed along with all the connections that +// feed the useless neurons. We'll cull neurons with .numOutputs == 0 +// or those that only feed themselves, i.e., .numSelfInputs == .numOutputs. +// Finally, we'll renumber the remaining neurons sequentially starting +// at zero using the .remappedNumber member. +struct Node { + uint16_t remappedNumber; + uint16_t numOutputs; + uint16_t numSelfInputs; + uint16_t numInputsFromSensorsOrOtherNeurons; +}; + + +// Two neuron renumberings occur: The original genome uses a uint16_t for +// neuron numbers. The first renumbering maps 16-bit unsigned neuron numbers +// to the range 0..p.maxNumberNeurons - 1. After culling useless neurons +// (see comments above), we'll renumber the remaining neurons sequentially +// starting at 0. +typedef std::map NodeMap; // key is neuron number 0..p.maxNumberNeurons - 1 + +typedef std::list ConnectionList; + + +// Returns by value a single gene with random members. +// See genome.h for the width of the members. +// ToDo: don't assume the width of the members in gene. +Gene makeRandomGene() +{ + Gene gene; + + gene.sourceType = randomUint() & 1; + gene.sourceNum = (uint16_t)randomUint(0, 0x7fff); + gene.sinkType = randomUint() & 1; + gene.sinkNum = (uint16_t)randomUint(0, 0x7fff); + gene.weight = Gene::makeRandomWeight(); + + return gene; +} + + +// Returns by value a single genome with random genes. +Genome makeRandomGenome() +{ + Genome genome; + + unsigned length = randomUint(p.genomeInitialLengthMin, p.genomeInitialLengthMax); + for (unsigned n = 0; n < length; ++n) { + genome.push_back(makeRandomGene()); + } + + return genome; +} + + +// Convert the indiv's genome to a renumbered connection list. +// This renumbers the neurons from their uint16_t values in the genome +// to the range 0..p.maxNumberNeurons - 1 by using a modulo operator. +// Sensors are renumbered 0..Sensor::NUM_SENSES - 1 +// Actions are renumbered 0..Action::NUM_ACTIONS - 1 +void makeRenumberedConnectionList(ConnectionList &connectionList, const Genome &genome) +{ + connectionList.clear(); + for (auto const &gene : genome) { + connectionList.push_back(gene); + auto &conn = connectionList.back(); + + if (conn.sourceType == NEURON) { + conn.sourceNum %= p.maxNumberNeurons; + } else { + conn.sourceNum %= Sensor::NUM_SENSES; + } + + if (conn.sinkType == NEURON) { + conn.sinkNum %= p.maxNumberNeurons; + } else { + conn.sinkNum %= Action::NUM_ACTIONS; + } + } +} + + +// Scan the connections and make a list of all the neuron numbers +// mentioned in the connections. Also keep track of how many inputs and +// outputs each neuron has. +void makeNodeList(NodeMap &nodeMap, const ConnectionList &connectionList) +{ + nodeMap.clear(); + + for (const Gene &conn : connectionList) { + if (conn.sinkType == NEURON) { + auto it = nodeMap.find(conn.sinkNum); + if (it == nodeMap.end()) { + assert(conn.sinkNum < p.maxNumberNeurons); + nodeMap.insert(std::pair(conn.sinkNum, {} )); + it = nodeMap.find(conn.sinkNum); + assert(it->first < p.maxNumberNeurons); + it->second.numOutputs = 0; + it->second.numSelfInputs = 0; + it->second.numInputsFromSensorsOrOtherNeurons = 0; + } + + if (conn.sourceType == NEURON && (conn.sourceNum == conn.sinkNum)) { + ++(it->second.numSelfInputs); + } else { + ++(it->second.numInputsFromSensorsOrOtherNeurons); + } + assert(nodeMap.count(conn.sinkNum) == 1); + } + if (conn.sourceType == NEURON) { + auto it = nodeMap.find(conn.sourceNum); + if (it == nodeMap.end()) { + assert(conn.sourceNum < p.maxNumberNeurons); + nodeMap.insert(std::pair(conn.sourceNum, {} )); + it = nodeMap.find(conn.sourceNum); + assert(it->first < p.maxNumberNeurons); + it->second.numOutputs = 0; + it->second.numSelfInputs = 0; + it->second.numInputsFromSensorsOrOtherNeurons = 0; + } + ++(it->second.numOutputs); + assert(nodeMap.count(conn.sourceNum) == 1); + } + } +} + + +// During the culling process, we will remove any neuron that has no outputs, +// and all the connections that feed the useless neuron. +void removeConnectionsToNeuron(ConnectionList &connections, NodeMap &nodeMap, uint16_t neuronNumber) +{ + for (auto itConn = connections.begin(); itConn != connections.end(); ) { + if (itConn->sinkType == NEURON && itConn->sinkNum == neuronNumber) { + // Remove the connection. If the connection source is from another + // neuron, also decrement the other neuron's numOutputs: + if (itConn->sourceType == NEURON) { + --(nodeMap[itConn->sourceNum].numOutputs); + } + itConn = connections.erase(itConn); + } else { + ++itConn; + } + } +} + + +// If a neuron has no outputs or only outputs that feed itself, then we +// remove it along with all connections that feed it. Reiterative, because +// after we remove a connection to a useless neuron, it may result in a +// different neuron having no outputs. +void cullUselessNeurons(ConnectionList &connections, NodeMap &nodeMap) +{ + bool allDone = false; + while (!allDone) { + allDone = true; + for (auto itNeuron = nodeMap.begin(); itNeuron != nodeMap.end(); ) { + assert(itNeuron->first < p.maxNumberNeurons); + // We're looking for neurons with zero outputs, or neurons that feed itself + // and nobody else: + if (itNeuron->second.numOutputs == itNeuron->second.numSelfInputs) { // could be 0 + allDone = false; + // Find and remove connections from sensors or other neurons + removeConnectionsToNeuron(connections, nodeMap, itNeuron->first); + itNeuron = nodeMap.erase(itNeuron); + } else { + ++itNeuron; + } + } + } +} + + +// This function is used when an agent is spawned. This function converts the +// agent's inherited genome into the agent's neural net brain. There is a close +// correspondence between the genome and the neural net, but a connection +// specified in the genome will not be represented in the neural net if the +// connection feeds a neuron that does not itself feed anything else. +// Neurons get renumbered in the process: +// 1. Create a set of referenced neuron numbers where each index is in the +// range 0..p.genomeMaxLength-1, keeping a count of outputs for each neuron. +// 2. Delete any referenced neuron index that has no outputs or only feeds itself. +// 3. Renumber the remaining neurons sequentially starting at 0. +void Indiv::createWiringFromGenome() +{ + NodeMap nodeMap; // list of neurons and their number of inputs and outputs + ConnectionList connectionList; // synaptic connections + + // Convert the indiv's genome to a renumbered connection list + makeRenumberedConnectionList(connectionList, genome); + + // Make a node (neuron) list from the renumbered connection list + makeNodeList(nodeMap, connectionList); + + // Find and remove neurons that don't feed anything or only feed themself. + // This reiteratively removes all connections to the useless neurons. + cullUselessNeurons(connectionList, nodeMap); + + // The neurons map now has all the referenced neurons, their neuron numbers, and + // the number of outputs for each neuron. Now we'll renumber the neurons + // starting at zero. + + assert(nodeMap.size() <= p.maxNumberNeurons); + uint16_t newNumber = 0; + for (auto & node : nodeMap) { + assert(node.second.numOutputs != 0); + node.second.remappedNumber = newNumber++; + } + + // Create the indiv's connection list in two passes: + // First the connections to neurons, then the connections to actions. + // This ordering optimizes the feed-forward function in feedForward.cpp. + + nnet.connections.clear(); + + // First, the connections from sensor or neuron to a neuron + for (auto const &conn : connectionList) { + if (conn.sinkType == NEURON) { + nnet.connections.push_back(conn); + auto &newConn = nnet.connections.back(); + // fix the destination neuron number + newConn.sinkNum = nodeMap[newConn.sinkNum].remappedNumber; + // if the source is a neuron, fix its number too + if (newConn.sourceType == NEURON) { + newConn.sourceNum = nodeMap[newConn.sourceNum].remappedNumber; + } + } + } + + // Last, the connections from sensor or neuron to an action + for (auto const &conn : connectionList) { + if (conn.sinkType == ACTION) { + nnet.connections.push_back(conn); + auto &newConn = nnet.connections.back(); + // if the source is a neuron, fix its number + if (newConn.sourceType == NEURON) { + newConn.sourceNum = nodeMap[newConn.sourceNum].remappedNumber; + } + } + } + + // Create the indiv's neural node list + nnet.neurons.clear(); + for (unsigned neuronNum = 0; neuronNum < nodeMap.size(); ++neuronNum) { + nnet.neurons.push_back( {} ); + nnet.neurons.back().output = initialNeuronOutput(); + nnet.neurons.back().driven = (nodeMap[neuronNum].numInputsFromSensorsOrOtherNeurons != 0); + } +} + + +// --------------------------------------------------------------------------- + + +// This applies a point mutation at a random bit in a genome. +void randomBitFlip(Genome &genome) +{ + int method = 1; + + unsigned byteIndex = randomUint(0, genome.size() - 1) * sizeof(Gene); + unsigned elementIndex = randomUint(0, genome.size() - 1); + uint8_t bitIndex8 = 1 << randomUint(0, 7); + + if (method == 0) { + ((uint8_t *)&genome[0])[byteIndex] ^= bitIndex8; + } else if (method == 1) { + float chance = randomUint() / (float)RANDOM_UINT_MAX; // 0..1 + if (chance < 0.2) { // sourceType + genome[elementIndex].sourceType ^= 1; + } else if (chance < 0.4) { // sinkType + genome[elementIndex].sinkType ^= 1; + } else if (chance < 0.6) { // sourceNum + genome[elementIndex].sourceNum ^= bitIndex8; + } else if (chance < 0.8) { // sinkNum + genome[elementIndex].sinkNum ^= bitIndex8; + } else { // weight + genome[elementIndex].weight ^= (1 << randomUint(1, 15)); + } + } else { + assert(false); + } +} + + +// If the genome is longer than the prescribed length, and if it's longer +// than one gene, then we remove genes from the front or back. This is +// used only when the simulator is configured to allow genomes of +// unequal lengths during a simulation. +void cropLength(Genome &genome, unsigned length) +{ + if (genome.size() > length && length > 0) { + if (randomUint() / (float)RANDOM_UINT_MAX < 0.5) { + // trim front + unsigned numberElementsToTrim = genome.size() - length; + genome.erase(genome.begin(), genome.begin() + numberElementsToTrim); + } else { + // trim back + genome.erase(genome.end() - (genome.size() - length), genome.end()); + } + } +} + + +// Inserts or removes a single gene from the genome. This is +// used only when the simulator is configured to allow genomes of +// unequal lengths during a simulation. +void randomInsertDeletion(Genome &genome) +{ + float probability = p.geneInsertionDeletionRate; + if (randomUint() / (float)RANDOM_UINT_MAX < probability) { + if (randomUint() / (float)RANDOM_UINT_MAX < p.deletionRatio) { + // deletion + if (genome.size() > 1) { + genome.erase(genome.begin() + randomUint(0, genome.size() - 1)); + } + } else if (genome.size() < p.genomeMaxLength) { + // insertion + //genome.insert(genome.begin() + randomUint(0, genome.size() - 1), makeRandomGene()); + genome.push_back(makeRandomGene()); + } + } +} + + +// This function causes point mutations in a genome with a probability defined +// by the parameter p.pointMutationRate. +void applyPointMutations(Genome &genome) +{ + unsigned numberOfGenes = genome.size(); + while (numberOfGenes-- > 0) { + if ((randomUint() / (float)RANDOM_UINT_MAX) < p.pointMutationRate) { + randomBitFlip(genome); + } + } +} + + +// This generates a child genome from one or two parent genomes. +// If the parameter p.sexualReproduction is true, two parents contribute +// genes to the offspring. The new genome may undergo mutation. +// Must be called in single-thread mode between generations +Genome generateChildGenome(const std::vector &parentGenomes) +{ + // random parent (or parents if sexual reproduction) with random + // mutations + Genome genome; + + uint16_t parent1Idx; + uint16_t parent2Idx; + + // Choose two parents randomly from the candidates. If the parameter + // p.chooseParentsByFitness is false, then we choose at random from + // all the candidate parents with equal preference. If the parameter is + // true, then we give preference to candidate parents according to their + // score. Their score was computed by the survival/selection algorithm + // in survival-criteria.cpp. + if (p.chooseParentsByFitness && parentGenomes.size() > 1) { + parent1Idx = randomUint(1, parentGenomes.size() - 1); + parent2Idx = randomUint(0, parent1Idx - 1); + } else { + parent1Idx = randomUint(0, parentGenomes.size() - 1); + parent2Idx = randomUint(0, parentGenomes.size() - 1); + } + + const Genome &g1 = parentGenomes[parent1Idx]; + const Genome &g2 = parentGenomes[parent2Idx]; + + if (g1.size() == 0 || g2.size() == 0) { + std::cout << "invalid genome" << std::endl; + assert(false); + } + + auto overlayWithSliceOf = [&](const Genome &gShorter) { + uint16_t index0 = randomUint(0, gShorter.size() - 1); + uint16_t index1 = randomUint(0, gShorter.size()); + if (index0 > index1) { + std::swap(index0, index1); + } + std::copy(&gShorter[index0], &gShorter[index1], &genome[index0]); + }; + + if (p.sexualReproduction) { + if (g1.size() > g2.size()) { + genome = g1; + overlayWithSliceOf(g2); + assert(genome.size() > 0); + } else { + genome = g2; + overlayWithSliceOf(g1); + assert(genome.size() > 0); + } + + // Trim to length = average length of parents + unsigned sum = g1.size() + g2.size(); + // If average length is not an integral number, add one half the time + if ((sum & 1) && (randomUint() & 1)) { + ++sum; + } + cropLength(genome, sum / 2); + assert(genome.size() > 0); + } else { + genome = g2; + assert(genome.size() > 0); + } + + randomInsertDeletion(genome); + assert(genome.size() > 0); + applyPointMutations(genome); + assert(genome.size() > 0); + assert(genome.size() <= p.genomeMaxLength); + + return genome; +} + +} // end namespace BS diff --git a/src/getSensor.cpp b/src/getSensor.cpp new file mode 100644 index 0000000..dbb2dac --- /dev/null +++ b/src/getSensor.cpp @@ -0,0 +1,357 @@ +// getSensor.cpp + +#include +#include +#include +#include "simulator.h" + +namespace BS { + +float getPopulationDensityAlongAxis(Coord loc, Dir dir) +{ + // Converts the population along the specified axis to the sensor range. The + // locations of neighbors are scaled by the inverse of their distance times + // the positive absolute cosine of the difference of their angle and the + // specified axis. The maximum positive or negative magnitude of the sum is + // about 2*radius. We don't adjust for being close to a border, so populations + // along borders and in corners are commonly sparser than away from borders. + // An empty neighborhood results in a sensor value exactly midrange; below + // midrange if the population density is greatest in the reverse direction, + // above midrange if density is greatest in forward direction. + + double sum = 0.0; + auto f = [&](Coord tloc) { + if (tloc != loc && grid.isOccupiedAt(tloc)) { + Coord offset = tloc - loc; + double anglePosCos = offset.raySameness(dir); + double dist = std::sqrt((double)offset.x * offset.x + (double)offset.y * offset.y); + double contrib = (1.0 / dist) * anglePosCos; + sum += contrib; + } + }; + + visitNeighborhood(loc, p.populationSensorRadius, f); + double maxSumMag = 6.0 * p.populationSensorRadius; + assert(sum >= -maxSumMag && sum <= maxSumMag); + + double sensorVal; + sensorVal = sum / maxSumMag; // convert to -1.0..1.0 + sensorVal = (sensorVal + 1.0) / 2.0; // convert to 0.0..1.0 + + return sensorVal; +} + + +// Converts the number of locations (not including loc) to the next barrier location +// along opposite directions of the specified axis to the sensor range. If no barriers +// are found, the result is sensor mid-range. Ignores agents in the path. +float getShortProbeBarrierDistance(Coord loc0, Dir dir, unsigned probeDistance) +{ + unsigned countFwd = 0; + unsigned countRev = 0; + Coord loc = loc0 + dir; + unsigned numLocsToTest = probeDistance; + // Scan positive direction + while (numLocsToTest > 0 && grid.isInBounds(loc) && !grid.isBarrierAt(loc)) { + ++countFwd; + loc = loc + dir; + --numLocsToTest; + } + if (numLocsToTest > 0 && !grid.isInBounds(loc)) { + countFwd = probeDistance; + } + // Scan negative direction + numLocsToTest = probeDistance; + loc = loc0 - dir; + while (numLocsToTest > 0 && grid.isInBounds(loc) && !grid.isBarrierAt(loc)) { + ++countRev; + loc = loc - dir; + --numLocsToTest; + } + if (numLocsToTest > 0 && !grid.isInBounds(loc)) { + countRev = probeDistance; + } + + float sensorVal = ((countFwd - countRev) + probeDistance); // convert to 0..2*probeDistance + sensorVal = (sensorVal / 2.0) / probeDistance; // convert to 0.0..1.0 + return sensorVal; +} + + +float getSignalDensity(unsigned layerNum, Coord loc) +{ + // returns magnitude of the specified signal layer in a neighborhood, with + // 0.0..maxSignalSum converted to the sensor range. + + unsigned countLocs = 0; + unsigned long sum = 0; + Coord center = loc; + + auto f = [&](Coord tloc) { + ++countLocs; + sum += signals.getMagnitude(layerNum, tloc); + }; + + visitNeighborhood(center, p.signalSensorRadius, f); + double maxSum = (float)countLocs * SIGNAL_MAX; + double sensorVal = sum / maxSum; // convert to 0.0..1.0 + + return sensorVal; +} + + +float getSignalDensityAlongAxis(unsigned layerNum, Coord loc, Dir dir) +{ + // Converts the signal density along the specified axis to sensor range. The + // values of cell signal levels are scaled by the inverse of their distance times + // the positive absolute cosine of the difference of their angle and the + // specified axis. The maximum positive or negative magnitude of the sum is + // about 2*radius*SIGNAL_MAX (?). We don't adjust for being close to a border, + // so signal densities along borders and in corners are commonly sparser than + // away from borders. + + double sum = 0.0; + auto f = [&](Coord tloc) { + if (tloc != loc) { + Coord offset = tloc - loc; + double anglePosCos = offset.raySameness(dir); + double dist = std::sqrt((double)offset.x * offset.x + (double)offset.y * offset.y); + double contrib = (1.0 / dist) * anglePosCos * signals.getMagnitude(layerNum, loc); + sum += contrib; + } + }; + + visitNeighborhood(loc, p.signalSensorRadius, f); + double maxSumMag = 6.0 * p.signalSensorRadius * SIGNAL_MAX; + assert(sum >= -maxSumMag && sum <= maxSumMag); + double sensorVal = sum / maxSumMag; // convert to -1.0..1.0 + sensorVal = (sensorVal + 1.0) / 2.0; // convert to 0.0..1.0 + + return sensorVal; +} + + +// Returns the number of locations to the next agent in the specified +// direction, not including loc. If the probe encounters a boundary or a +// barrier before reaching the longProbeDist distance, returns longProbeDist. +// Returns 0..longProbeDist. +unsigned longProbePopulationFwd(Coord loc, Dir dir, unsigned longProbeDist) +{ + assert(longProbeDist > 0); + unsigned count = 0; + loc = loc + dir; + unsigned numLocsToTest = longProbeDist; + while (numLocsToTest > 0 && grid.isInBounds(loc) && grid.isEmptyAt(loc)) { + ++count; + loc = loc + dir; + --numLocsToTest; + } + if (numLocsToTest > 0 && (!grid.isInBounds(loc) || grid.isBarrierAt(loc))) { + return longProbeDist; + } else { + return count; + } +} + + +// Returns the number of locations to the next barrier in the +// specified direction, not including loc. Ignores agents in the way. +// If the distance to the border is less than the longProbeDist distance +// and no barriers are found, returns longProbeDist. +// Returns 0..longProbeDist. +unsigned longProbeBarrierFwd(Coord loc, Dir dir, unsigned longProbeDist) +{ + assert(longProbeDist > 0); + unsigned count = 0; + loc = loc + dir; + unsigned numLocsToTest = longProbeDist; + while (numLocsToTest > 0 && grid.isInBounds(loc) && !grid.isBarrierAt(loc)) { + ++count; + loc = loc + dir; + --numLocsToTest; + } + if (numLocsToTest > 0 && !grid.isInBounds(loc)) { + return longProbeDist; + } else { + return count; + } +} + + +// Returned sensor values range SENSOR_MIN..SENSOR_MAX +float Indiv::getSensor(Sensor sensorNum, unsigned simStep) const +{ + float sensorVal; + + switch (sensorNum) { + case Sensor::AGE: + // Converts age (units of simSteps compared to life expectancy) + // linearly to normalized sensor range 0.0..1.0 + sensorVal = (float)age / p.stepsPerGeneration; + break; + case Sensor::BOUNDARY_DIST: + { + // Finds closest boundary, compares that to the max possible dist + // to a boundary from the center, and converts that linearly to the + // sensor range 0.0..1.0 + int distX = std::min(loc.x, (p.sizeX - loc.x) - 1); + int distY = std::min(loc.y, (p.sizeY - loc.y) - 1); + int closest = std::min(distX, distY); + int maxPossible = std::max(p.sizeX / 2 - 1, p.sizeY / 2 - 1); + sensorVal = (float)closest / maxPossible; + break; + } + case Sensor::BOUNDARY_DIST_X: + { + // Measures the distance to nearest boundary in the east-west axis, + // max distance is half the grid width; scaled to sensor range 0.0..1.0. + int minDistX = std::min(loc.x, (p.sizeX - loc.x) - 1); + sensorVal = minDistX / (p.sizeX / 2.0); + break; + } + case Sensor::BOUNDARY_DIST_Y: + { + // Measures the distance to nearest boundary in the south-north axis, + // max distance is half the grid height; scaled to sensor range 0.0..1.0. + int minDistY = std::min(loc.y, (p.sizeY - loc.y) - 1); + sensorVal = minDistY / (p.sizeY / 2.0); + break; + } + case Sensor::LAST_MOVE_DIR_X: + { + // X component -1,0,1 maps to sensor values 0.0, 0.5, 1.0 + auto lastX = lastMoveDir.asNormalizedCoord().x; + sensorVal = lastX == 0 ? 0.5 : + (lastX == -1 ? 0.0 : 1.0); + break; + } + case Sensor::LAST_MOVE_DIR_Y: + { + // Y component -1,0,1 maps to sensor values 0.0, 0.5, 1.0 + auto lastY = lastMoveDir.asNormalizedCoord().y; + sensorVal = lastY == 0 ? 0.5 : + (lastY == -1 ? 0.0 : 1.0); + break; + } + case Sensor::LOC_X: + // Maps current X location 0..p.sizeX-1 to sensor range 0.0..1.0 + sensorVal = (float)loc.x / (p.sizeX - 1); + break; + case Sensor::LOC_Y: + // Maps current Y location 0..p.sizeY-1 to sensor range 0.0..1.0 + sensorVal = (float)loc.y / (p.sizeY - 1); + break; + case Sensor::OSC1: + { + // Maps the oscillator sine wave to sensor range 0.0..1.0; + // cycles starts at simStep 0 for everbody. + float phase = (simStep % oscPeriod) / (float)oscPeriod; // 0.0..1.0 + float factor = -std::cos(phase * 2.0f * 3.1415927f); + assert(factor >= -1.0f && factor <= 1.0f); + factor += 1.0f; // convert to 0.0..2.0 + factor /= 2.0; // convert to 0.0..1.0 + sensorVal = factor; + // Clip any round-off error + sensorVal = std::min(1.0, std::max(0.0, sensorVal)); + break; + } + case Sensor::LONGPROBE_POP_FWD: + { + // Measures the distance to the nearest other individual in the + // forward direction. If non found, returns the maximum sensor value. + // Maps the result to the sensor range 0.0..1.0. + sensorVal = longProbePopulationFwd(loc, lastMoveDir, longProbeDist) / (float)longProbeDist; // 0..1 + break; + } + case Sensor::LONGPROBE_BAR_FWD: + { + // Measures the distance to the nearest barrier in the forward + // direction. If non found, returns the maximum sensor value. + // Maps the result to the sensor range 0.0..1.0. + sensorVal = longProbeBarrierFwd(loc, lastMoveDir, longProbeDist) / (float)longProbeDist; // 0..1 + break; + } + case Sensor::POPULATION: + { + // Returns population density in neighborhood converted linearly from + // 0..100% to sensor range + unsigned countLocs = 0; + unsigned countOccupied = 0; + Coord center = loc; + + auto f = [&](Coord tloc) { + ++countLocs; + if (grid.isOccupiedAt(tloc)) { + ++countOccupied; + } + }; + + visitNeighborhood(center, p.populationSensorRadius, f); + sensorVal = (float)countOccupied / countLocs; + break; + } + case Sensor::POPULATION_FWD: + // Sense population density along axis of last movement direction, mapped + // to sensor range 0.0..1.0 + sensorVal = getPopulationDensityAlongAxis(loc, lastMoveDir); + break; + case Sensor::POPULATION_LR: + // Sense population density along an axis 90 degrees from last movement direction + sensorVal = getPopulationDensityAlongAxis(loc, lastMoveDir.rotate90DegCW()); + break; + case Sensor::BARRIER_FWD: + // Sense the nearest barrier along axis of last movement direction, mapped + // to sensor range 0.0..1.0 + sensorVal = getShortProbeBarrierDistance(loc, lastMoveDir, p.shortProbeBarrierDistance); + break; + case Sensor::BARRIER_LR: + // Sense the nearest barrier along axis perpendicular to last movement direction, mapped + // to sensor range 0.0..1.0 + sensorVal = getShortProbeBarrierDistance(loc, lastMoveDir.rotate90DegCW(), p.shortProbeBarrierDistance); + break; + case Sensor::RANDOM: + // Returns a random sensor value in the range 0.0..1.0. + sensorVal = randomUint() / (float)UINT_MAX; + break; + case Sensor::SIGNAL0: + // Returns magnitude of signal0 in the local neighborhood, with + // 0.0..maxSignalSum converted to sensorRange 0.0..1.0 + sensorVal = getSignalDensity(0, loc); + break; + case Sensor::SIGNAL0_FWD: + // Sense signal0 density along axis of last movement direction + sensorVal = getSignalDensityAlongAxis(0, loc, lastMoveDir); + break; + case Sensor::SIGNAL0_LR: + // Sense signal0 density along an axis perpendicular to last movement direction + sensorVal = getSignalDensityAlongAxis(0, loc, lastMoveDir.rotate90DegCW()); + break; + case Sensor::GENETIC_SIM_FWD: + { + // Return minimum sensor value if nobody is alive in the forward adjacent location, + // else returns a similarity match in the sensor range 0.0..1.0 + Coord loc2 = loc + lastMoveDir; + if (grid.isInBounds(loc2) && grid.isOccupiedAt(loc2)) { + const Indiv &indiv2 = peeps.getIndiv(loc2); + if (indiv2.alive) { + sensorVal = genomeSimilarity(genome, indiv2.genome); // 0.0..1.0 + } + } + break; + } + default: + assert(false); + break; + } + +if (std::isnan(sensorVal) || sensorVal < -0.01 || sensorVal > 1.01) { + std::cout << "sensorVal=" << (int)sensorVal << " for " << sensorName((Sensor)sensorNum) << std::endl; + sensorVal = std::max(0.0f, std::min(sensorVal, 1.0f)); // clip +} + + assert(!std::isnan(sensorVal) && sensorVal >= -0.01 && sensorVal <= 1.01); + + return sensorVal; +} + +} // end namespace BS diff --git a/src/grid.cpp b/src/grid.cpp new file mode 100644 index 0000000..e094dd4 --- /dev/null +++ b/src/grid.cpp @@ -0,0 +1,50 @@ +// grid.cpp + +#include +#include +#include "simulator.h" + +namespace BS { + +// Allocates space for the 2D grid +void Grid::init(uint16_t sizeX, uint16_t sizeY) +{ + auto col = Column(sizeY); + data = std::vector(sizeX, col); +} + + +// Finds a random unoccupied location in the grid +Coord Grid::findEmptyLocation() const { + Coord loc; + bool found = false; + + while (!found) { + loc.x = randomUint(0, p.sizeX - 1); + loc.y = randomUint(0, p.sizeY - 1); + if (grid.isEmptyAt(loc)) { + break; + } + } + return loc; +} + + +// This is a utility function used when inspecting a local neighborhood around +// some location. This function feeds each valid (in-bounds) location in the specified +// neighborhood to the specified function. Locations include self (center of the neighborhood). +void visitNeighborhood(Coord loc, float radius, std::function f) +{ + for (int dx = -std::min(radius, loc.x); dx <= std::min(radius, (p.sizeX - loc.x) - 1); ++dx) { + int16_t x = loc.x + dx; + assert(x >= 0 && x < p.sizeX); + int extentY = (int)sqrt(radius * radius - dx * dx); + for (int dy = -std::min(extentY, loc.y); dy <= std::min(extentY, (p.sizeY - loc.y) - 1); ++dy) { + int16_t y = loc.y + dy; + assert(y >= 0 && y < p.sizeY); + f( Coord { x, y} ); + } + } +} + +} // end namespace BS diff --git a/src/grid.h b/src/grid.h new file mode 100644 index 0000000..b0b3ad2 --- /dev/null +++ b/src/grid.h @@ -0,0 +1,71 @@ +#ifndef GRID_H_INCLUDED +#define GRID_H_INCLUDED + +// The grid is the 2D arena where the agents live. + +#include +#include +#include +#include "basicTypes.h" + +namespace BS { + +// Grid is a somewhat dumb 2D container of unsigned 16-bit values. +// Grid understands that the elements are either EMPTY, BARRIER, or +// otherwise an index value into the peeps container. +// The elements are allocated and cleared to EMPTY in the ctor. +// Prefer .at() and .set() for random element access. Or use Grid[x][y] +// for direct access where the y index is the inner loop. +// Element values are not otherwise interpreted by class Grid. + +const uint16_t EMPTY = 0; // Index value 0 is reserved +const uint16_t BARRIER = 0xffff; + +class Grid { +public: + // Column order here allows us to access grid elements as data[x][y] + // while thinking of x as column and y as row + struct Column { + Column(uint16_t numRows) : data { std::vector(numRows, 0) } { } + void zeroFill() { std::fill(data.begin(), data.end(), 0); } + uint16_t& operator[](uint16_t rowNum) { return data[rowNum]; } + uint16_t operator[](uint16_t rowNum) const { return data[rowNum]; } + size_t size() const { return data.size(); } + private: + std::vector data; + }; + + void init(uint16_t sizeX, uint16_t sizeY); + void zeroFill() { for (Column &column : data) column.zeroFill(); } + uint16_t sizeX() const { return data.size(); } + uint16_t sizeY() const { return data[0].size(); } + bool isInBounds(Coord loc) const { return loc.x >= 0 && loc.x < sizeX() && loc.y >= 0 && loc.y < sizeY(); } + bool isEmptyAt(Coord loc) const { return at(loc) == EMPTY; } + bool isBarrierAt(Coord loc) const { return at(loc) == BARRIER; } + // Occupied means an agent is living there. + bool isOccupiedAt(Coord loc) const { return at(loc) != EMPTY && at(loc) != BARRIER; } + bool isBorder(Coord loc) const { return loc.x == 0 || loc.x == sizeX() - 1 || loc.y == 0 || loc.y == sizeY() - 1; } + uint16_t at(Coord loc) const { return data[loc.x][loc.y]; }; + uint16_t at(uint16_t x, uint16_t y) const { return data[x][y]; } + + void set(Coord loc, uint16_t val) { data[loc.x][loc.y] = val; } + void set(uint16_t x, uint16_t y, uint16_t val) { data[x][y] = val; } + Coord findEmptyLocation() const; + void createBarrier(unsigned barrierType); + const std::vector &getBarrierLocations() const { return barrierLocations; } + const std::vector &getBarrierCenters() const { return barrierCenters; } + // Direct access: + Column & operator[](uint16_t columnXNum) { return data[columnXNum]; } + const Column & operator[](uint16_t columnXNum) const { return data[columnXNum]; } +private: + std::vector data; + std::vector barrierLocations; + std::vector barrierCenters; +}; + +extern void visitNeighborhood(Coord loc, float radius, std::function f); +extern void unitTestGridVisitNeighborhood(); + +} // end namespace BS + +#endif // GRID_H_INCLUDED diff --git a/src/imageWriter.cpp b/src/imageWriter.cpp new file mode 100644 index 0000000..e324e4b --- /dev/null +++ b/src/imageWriter.cpp @@ -0,0 +1,260 @@ +// imageWriter.cpp + +#include +#include +#include +#include +#include +#include +#include +#include "simulator.h" +#include "imageWriter.h" +#define cimg_use_opencv 1 +#include "CImg.h" + +namespace BS { + +cimg_library::CImgList imageList; + +// Pushes a new image frame onto .imageList. +// +void saveOneFrameImmed(const ImageFrameData &data) +{ + using namespace cimg_library; + + CImg image(p.sizeX * p.displayScale, p.sizeY * p.displayScale, + 1, // Z depth + 3, // color channels + 255); // initial value + uint8_t color[3]; + std::stringstream imageFilename; + imageFilename << p.imageDir << "frame-" + << std::setfill('0') << std::setw(6) << data.generation + << '-' << std::setfill('0') << std::setw(6) << data.simStep + << ".png"; + + // Draw barrier locations + + color[0] = color[1] = color[2] = 0x88; + for (Coord loc : data.barrierLocs) { + image.draw_rectangle( + loc.x * p.displayScale - (p.displayScale / 2), ((p.sizeY - loc.y) - 1) * p.displayScale - (p.displayScale / 2), + (loc.x + 1) * p.displayScale, ((p.sizeY - (loc.y - 0))) * p.displayScale, + color, // rgb + 1.0); // alpha + } + + // Draw agents + + constexpr uint8_t maxColorVal = 0xb0; + constexpr uint8_t maxLumaVal = 0xb0; + + auto rgbToLuma = [](uint8_t r, uint8_t g, uint8_t b) { return (r+r+r+b+g+g+g+g) / 8; }; + + for (size_t i = 0; i < data.indivLocs.size(); ++i) { + int c = data.indivColors[i]; + color[0] = (c); // R: 0..255 + color[1] = ((c & 0x1f) << 3); // G: 0..255 + color[2] = ((c & 7) << 5); // B: 0..255 + + // Prevent color mappings to very bright colors (hard to see): + if (rgbToLuma(color[0], color[1], color[2]) > maxLumaVal) { + if (color[0] > maxColorVal) color[0] %= maxColorVal; + if (color[1] > maxColorVal) color[1] %= maxColorVal; + if (color[2] > maxColorVal) color[2] %= maxColorVal; + } + + image.draw_circle( + data.indivLocs[i].x * p.displayScale, + ((p.sizeY - data.indivLocs[i].y) - 1) * p.displayScale, + p.agentSize, + color, // rgb + 1.0); // alpha + } + + //image.save_png(imageFilename.str().c_str(), 3); + imageList.push_back(image); + + //CImgDisplay local(image, "biosim3"); +} + + +// Starts the image writer asynchronous thread. +ImageWriter::ImageWriter() + : droppedFrameCount{0}, busy{true}, dataReady{false}, + abortRequested{false} +{ + startNewGeneration(); +} + + +void ImageWriter::startNewGeneration() +{ + imageList.clear(); + skippedFrames = 0; +} + + +uint8_t makeGeneticColor(const Genome &genome) +{ + return ((genome.size() & 1) + | ((genome.front().sourceType) << 1) + | ((genome.back().sourceType) << 2) + | ((genome.front().sinkType) << 3) + | ((genome.back().sinkType) << 4) + | ((genome.front().sourceNum & 1) << 5) + | ((genome.front().sinkNum & 1) << 6) + | ((genome.back().sourceNum & 1) << 7)); +} + + +// This is a synchronous gate for giving a job to saveFrameThread(). +// Called from the same thread as the main simulator loop thread during +// single-thread mode. +// Returns true if the image writer accepts the job; returns false +// if the image writer is busy. Always called from a single thread +// and communicates with a single saveFrameThread(), so no need to make +// a critical section to safeguard the busy flag. When this function +// sets the busy flag, the caller will immediate see it, so the caller +// won't call again until busy is clear. When the thread clears the busy +// flag, it doesn't matter if it's not immediately visible to this +// function: there's no consequence other than a harmless frame-drop. +// The condition variable allows the saveFrameThread() to wait until +// there's a job to do. +bool ImageWriter::saveVideoFrame(unsigned simStep, unsigned generation) +{ + if (!busy) { + busy = true; + // queue job for saveFrameThread() + // We cache a local copy of data from params, grid, and peeps because + // those objects will change by the main thread at the same time our + // saveFrameThread() is using it to output a video frame. + data.simStep = simStep; + data.generation = generation; + data.indivLocs.clear(); + data.indivColors.clear(); + data.barrierLocs.clear(); + data.signalLayers.clear(); + //todo!!! + for (uint16_t index = 1; index <= p.population; ++index) { + const Indiv &indiv = peeps[index]; + if (indiv.alive) { + data.indivLocs.push_back(indiv.loc); + data.indivColors.push_back(makeGeneticColor(indiv.genome)); + } + } + + auto const &barrierLocs = grid.getBarrierLocations(); + for (Coord loc : barrierLocs) { + data.barrierLocs.push_back(loc); + } + + // tell thread there's a job to do + { + std::lock_guard lck(mutex_); + dataReady = true; + } + condVar.notify_one(); + return true; + } else { + // image saver thread is busy, drop a frame + ++droppedFrameCount; + return false; + } +} + + +// Synchronous version, always returns true +bool ImageWriter::saveVideoFrameSync(unsigned simStep, unsigned generation) +{ + // We cache a local copy of data from params, grid, and peeps because + // those objects will change by the main thread at the same time our + // saveFrameThread() is using it to output a video frame. + data.simStep = simStep; + data.generation = generation; + data.indivLocs.clear(); + data.indivColors.clear(); + data.barrierLocs.clear(); + data.signalLayers.clear(); + //todo!!! + for (uint16_t index = 1; index <= p.population; ++index) { + const Indiv &indiv = peeps[index]; + if (indiv.alive) { + data.indivLocs.push_back(indiv.loc); + data.indivColors.push_back(makeGeneticColor(indiv.genome)); + } + } + + auto const &barrierLocs = grid.getBarrierLocations(); + for (Coord loc : barrierLocs) { + data.barrierLocs.push_back(loc); + } + + saveOneFrameImmed(data); + return true; +} + + +// ToDo: put save_video() in its own thread +void ImageWriter::saveGenerationVideo(unsigned generation) +{ + if (imageList.size() > 0) { + std::stringstream videoFilename; + videoFilename << p.imageDir.c_str() << "/gen-" + << std::setfill('0') << std::setw(6) << generation + << ".avi"; + cv::setNumThreads(2); + imageList.save_video(videoFilename.str().c_str(), + 25, + "H264"); + if (skippedFrames > 0) { + std::cout << "Video skipped " << skippedFrames << " frames" << std::endl; + } + } + startNewGeneration(); +} + + +void ImageWriter::abort() +{ + busy =true; + abortRequested = true; + { + std::lock_guard lck(mutex_); + dataReady = true; + } + condVar.notify_one(); +} + + +// Runs in a thread; wakes up when there's a video frame to generate. +// When this wakes up, local copies of Params and Peeps will have been +// cached for us to use. +void ImageWriter::saveFrameThread() +{ + busy = false; // we're ready for business + std::cout << "Imagewriter thread started." << std::endl; + + while (true) { + // wait for job on queue + std::unique_lock lck(mutex_); + condVar.wait(lck, [&]{ return dataReady && busy; }); + // save frame + dataReady = false; + busy = false; + + if (abortRequested) { + break; + } + + // save image frame + saveOneFrameImmed(imageWriter.data); + + //std::cout << "Image writer thread waiting..." << std::endl; + //std::this_thread::sleep_for(std::chrono::seconds(2)); + + } + std::cout << "Image writer thread exiting." << std::endl; +} + +} // end namespace BS diff --git a/src/imageWriter.h b/src/imageWriter.h new file mode 100644 index 0000000..5ef45f8 --- /dev/null +++ b/src/imageWriter.h @@ -0,0 +1,55 @@ +#ifndef IMAGEWRITER_H_INCLUDED +#define IMAGEWRITER_H_INCLUDED + +// Creates a graphic frame for each simStep, then +// assembles them into a video at the end of a generation. + +#include +#include +#include +#include +#include "indiv.h" +#include "params.h" +#include "peeps.h" + +namespace BS { + +// This holds all data needed to construct one image frame. The data is +// cached in this structure so that the image writer can work on it in +// a separate thread while the main thread starts a new simstep. +struct ImageFrameData { + unsigned simStep; + unsigned generation; + std::vector indivLocs; + std::vector indivColors; + std::vector barrierLocs; + typedef std::vector> SignalLayer; // [x][y] + std::vector signalLayers; // [layer][x][y] +}; + + +struct ImageWriter { + ImageWriter(); + void startNewGeneration(); + bool saveVideoFrame(unsigned simStep, unsigned generation); + bool saveVideoFrameSync(unsigned simStep, unsigned generation); + void saveGenerationVideo(unsigned generation); + void abort(); + void saveFrameThread(); // runs in a thread + std::atomic droppedFrameCount; +private: + std::atomic busy; + std::mutex mutex_; + std::condition_variable condVar; + bool dataReady; + + ImageFrameData data; + bool abortRequested; + unsigned skippedFrames; +}; + +extern ImageWriter imageWriter; + +} // end namespace BS + +#endif // IMAGEWRITER_H_INCLUDED diff --git a/src/indiv.cpp b/src/indiv.cpp new file mode 100644 index 0000000..55962e0 --- /dev/null +++ b/src/indiv.cpp @@ -0,0 +1,30 @@ +// indiv.cpp + +#include +#include +#include "simulator.h" + +namespace BS { + +// This is called when any individual is spawned. +// The responsiveness parameter will be initialized here to maximum value +// of 1.0, then depending on which action activation function is used, +// the default undriven value may be changed to 1.0 or action midrange. +void Indiv::initialize(uint16_t index_, Coord loc_, const Genome &&genome_) +{ + index = index_; + loc = loc_; + //birthLoc = loc_; + grid.set(loc_, index_); + age = 0; + oscPeriod = 34; // ToDo !!! define a constant + alive = true; + lastMoveDir = Dir::random8(); + responsiveness = 0.5; // range 0.0..1.0 + longProbeDist = p.longProbeDistance; + challengeBits = (unsigned)false; // will be set true when some task gets accomplished + genome = std::move(genome_); + createWiringFromGenome(); +} + +} // end namespace BS diff --git a/src/indiv.h b/src/indiv.h new file mode 100644 index 0000000..65a6652 --- /dev/null +++ b/src/indiv.h @@ -0,0 +1,40 @@ +#ifndef INDIV_H_INCLUDED +#define INDIV_H_INCLUDED + +// Indiv is the structure that represents one individual agent. + +#include +#include +#include +#include "basicTypes.h" +#include "genome-neurons.h" + +namespace BS { + +// Also see class Peeps. + +struct Indiv { + bool alive; + uint16_t index; // index into peeps[] container + Coord loc; // refers to a location in grid[][] + Coord birthLoc; + unsigned age; + Genome genome; + NeuralNet nnet; // derived from .genome + float responsiveness; // 0.0..1.0 (0 is like asleep) + unsigned oscPeriod; // 2..4*p.stepsPerGeneration (TBD, see executeActions()) + unsigned longProbeDist; // distance for long forward probe for obstructions + Dir lastMoveDir; // direction of last movement + unsigned challengeBits; // modified when the indiv accomplishes some task + std::array feedForward(unsigned simStep); // reads sensors, returns actions + float getSensor(Sensor, unsigned simStep) const; + void initialize(uint16_t index, Coord loc, const Genome &&genome); + void createWiringFromGenome(); // creates .nnet member from .genome member + void printNeuralNet() const; + void printIGraphEdgeList() const; + void printGenome() const; +}; + +} // end namespace BS + +#endif // INDIV_H_INCLUDED diff --git a/src/main.cpp b/src/main.cpp new file mode 100644 index 0000000..b5bdd2e --- /dev/null +++ b/src/main.cpp @@ -0,0 +1,28 @@ +#include + +// This is included here only for the purpose of unit testing of basic types +#include "basicTypes.h" + +// This is the only call needed to start the simulator. No header file +// because it's just one declaration. Pass it argc and argv as if from main(). +// If there is no command line argument, the simulator will read the default +// config file ("biosim4.ini" in the current directory) to get the simulation +// parameters for this run. If there are one or more command line args, then +// argv[1] must contain the name of the config file which will be read instead +// of biosim4.ini. Any args after that are ignored. The simulator code is +// in namespace BS (for "biosim"). +namespace BS { + void simulator(int argc, char **argv); +} + + +int main(int argc, char **argv) +{ + BS::unitTestBasicTypes(); // called only for unit testing of basic types + + // Start the simulator with optional config filename (default "biosim4.ini"). + // See simulator.cpp and simulator.h. + BS::simulator(argc, argv); + + return 0; +} diff --git a/src/params.cpp b/src/params.cpp new file mode 100644 index 0000000..603894f --- /dev/null +++ b/src/params.cpp @@ -0,0 +1,300 @@ +// params.cpp +// See params.h for notes. + +#include +#include +#include +#include +#include +#include +#include +#include +#include "params.h" + +// To add a new parameter: +// 1. Add a member to struct Params in params.h. +// 2. Add a member and its default value to privParams in ParamManager::setDefault() +// in params.cpp. +// 3. Add an else clause to ParamManager::ingestParameter() in params.cpp. +// 4. Add a line to the user's parameter file (default name biosim4.ini) + +namespace BS { + +void ParamManager::setDefaults() +{ + privParams.sizeX = 128; + privParams.sizeY = 128; + privParams.challenge = 0; + + privParams.genomeInitialLengthMin = 16; + privParams.genomeInitialLengthMax = 16; + privParams.genomeMaxLength = 20; + privParams.logDir = "./logs/"; + privParams.imageDir = "./images/"; + privParams.population = 100; + privParams.stepsPerGeneration = 100; + privParams.maxGenerations = 100; + privParams.barrierType = 0; + privParams.replaceBarrierType = 0; + privParams.replaceBarrierTypeGenerationNumber = (uint32_t)-1; + privParams.numThreads = 1; + privParams.signalLayers = 1; + privParams.maxNumberNeurons = privParams.genomeMaxLength / 2; + privParams.pointMutationRate = 0.0001; + privParams.geneInsertionDeletionRate = 0.0001; + privParams.deletionRatio = 0.7; + privParams.killEnable = false; + privParams.sexualReproduction = true; + privParams.chooseParentsByFitness = true; + privParams.populationSensorRadius = 2.0; + privParams.signalSensorRadius = 1; + privParams.responsiveness = 0.5; + privParams.responsivenessCurveKFactor = 2; + privParams.longProbeDistance = 16; + privParams.shortProbeBarrierDistance = 3; + privParams.valenceSaturationMag = 0.5; + privParams.saveVideo = true; + privParams.videoStride = 1; + privParams.videoSaveFirstFrames = 0; + privParams.displayScale = 1; + privParams.agentSize = 2; + privParams.genomeAnalysisStride = 1; + privParams.displaySampleGenomes = 0; + privParams.genomeComparisonMethod = 1; + privParams.updateGraphLog = false; + privParams.updateGraphLogStride = 16; + privParams.graphLogUpdateCommand = "/usr/bin/gnuplot --persist ./tools/graphlog.gp"; +} + + +void ParamManager::registerConfigFile(const char *filename) +{ + configFilename = std::string(filename); +} + + +bool checkIfUint(const std::string &s) +{ + return s.find_first_not_of("0123456789") == std::string::npos; +} + + +bool checkIfInt(const std::string &s) +{ + //return s.find_first_not_of("-0123456789") == std::string::npos; + std::istringstream iss(s); + int i; + iss >> std::noskipws >> i; // noskipws considers leading whitespace invalid + // Check the entire string was consumed and if either failbit or badbit is set + return iss.eof() && !iss.fail(); +} + + +bool checkIfFloat(const std::string &s) +{ + std::istringstream iss(s); + double d; + iss >> std::noskipws >> d; // noskipws considers leading whitespace invalid + // Check the entire string was consumed and if either failbit or badbit is set + return iss.eof() && !iss.fail(); +} + + +bool checkIfBool(const std::string &s) +{ + return s == "0" || s == "1" || s == "true" || s == "false"; +} + + +bool getBoolVal(const std::string &s) +{ + if (s == "true" || s == "1") + return true; + else if (s == "false" || s == "0") + return false; + else + return false; +} + + +void ParamManager::ingestParameter(std::string name, std::string val) +{ + std::transform(name.begin(), name.end(), name.begin(), + [](unsigned char c){ return std::tolower(c); }); + //std::cout << name << " " << val << '\n' << std::endl; + + bool isUint = checkIfUint(val); + unsigned uVal = isUint ? (unsigned)std::stol(val.c_str()) : 0; + bool isInt = checkIfInt(val); + int iVal = isInt ? std::stoi(val.c_str()) : 0; + bool isFloat = checkIfFloat(val); + double dVal = isFloat ? std::stod(val.c_str()) : 0.0; + bool isBool = checkIfBool(val); + bool bVal = getBoolVal(val); + + do { + if (name == "sizex" && isUint && uVal >= 2 && uVal <= (uint16_t)-1) { + privParams.sizeX = uVal; break; + } + else if (name == "sizey" && isUint && uVal >= 2 && uVal <= (uint16_t)-1) { + privParams.sizeY = uVal; break; + } + else if (name == "challenge" && isUint && uVal < (uint16_t)-1) { + privParams.challenge = uVal; break; + } + else if (name == "genomeinitiallengthmin" && isUint && uVal > 0 && uVal < (uint16_t)-1) { + privParams.genomeInitialLengthMin = uVal; break; + } + else if (name == "genomeinitiallengthmax" && isUint && uVal > 0 && uVal < (uint16_t)-1) { + privParams.genomeInitialLengthMax = uVal; break; + } + else if (name == "logdir") { + privParams.logDir = val; break; + } + else if (name == "imagedir") { + privParams.imageDir = val; break; + } + else if (name == "population" && isUint && uVal > 0 && uVal < (uint32_t)-1) { + privParams.population = uVal; break; + } + else if (name == "stepspergeneration" && isUint && uVal > 0 && uVal < (uint16_t)-1) { + privParams.stepsPerGeneration = uVal; break; + } + else if (name == "maxgenerations" && isUint && uVal > 0 && uVal < 0x7fffffff) { + privParams.maxGenerations = uVal; break; + } + else if (name == "barriertype" && isUint && uVal < (uint32_t)-1) { + privParams.barrierType = uVal; break; + } + else if (name == "replacebarriertype" && isUint && uVal < (uint32_t)-1) { + privParams.replaceBarrierType = uVal; break; + } + else if (name == "replacebarriertypegenerationnumber" && isInt && iVal >= -1) { + privParams.replaceBarrierTypeGenerationNumber = (iVal == -1 ? (uint32_t)-1 : iVal); break; + } + else if (name == "numthreads" && isUint && uVal > 0 && uVal < (uint16_t)-1) { + privParams.numThreads = uVal; break; + } + else if (name == "signallayers" && isUint && uVal < (uint16_t)-1) { + privParams.signalLayers = uVal; break; + } + else if (name == "genomemaxlength" && isUint && uVal > 0 && uVal < (uint16_t)-1) { + privParams.genomeMaxLength = uVal; break; + } + else if (name == "maxnumberneurons" && isUint && uVal > 0 && uVal < (uint16_t)-1) { + privParams.maxNumberNeurons = uVal; break; + } + else if (name == "pointmutationrate" && isFloat && dVal >= 0.0 && dVal <= 1.0) { + privParams.pointMutationRate = dVal; break; + } + else if (name == "geneinsertiondeletionrate" && isFloat && dVal >= 0.0 && dVal <= 1.0) { + privParams.geneInsertionDeletionRate = dVal; break; + } + else if (name == "deletionratio" && isFloat && dVal >= 0.0 && dVal <= 1.0) { + privParams.deletionRatio = dVal; break; + } + else if (name == "killenable" && isBool) { + privParams.killEnable = bVal; break; + } + else if (name == "sexualreproduction" && isBool) { + privParams.sexualReproduction = bVal; break; + } + else if (name == "chooseparentsbyfitness" && isBool) { + privParams.chooseParentsByFitness = bVal; break; + } + else if (name == "populationsensorradius" && isFloat && dVal > 0.0) { + privParams.populationSensorRadius = dVal; break; + } + else if (name == "signalsensorradius" && isFloat && dVal > 0.0) { + privParams.signalSensorRadius = dVal; break; + } + else if (name == "responsiveness" && isFloat && dVal >= 0.0) { + privParams.responsiveness = dVal; break; + } + else if (name == "responsivenesscurvekfactor" && isUint && uVal >= 1 && uVal <= 20) { + privParams.responsivenessCurveKFactor = uVal; break; + } + else if (name == "longprobedistance" && isUint && uVal > 0) { + privParams.longProbeDistance = uVal; break; + } + else if (name == "shortprobebarrierdistance" && isUint && uVal > 0) { + privParams.shortProbeBarrierDistance = uVal; break; + } + else if (name == "valencesaturationmag" && isFloat && dVal >= 0.0) { + privParams.valenceSaturationMag = dVal; break; + } + else if (name == "savevideo" && isBool) { + privParams.saveVideo = bVal; break; + } + else if (name == "videostride" && isUint && uVal > 0) { + privParams.videoStride = uVal; break; + } + else if (name == "videosavefirstframes" && isUint) { + privParams.videoSaveFirstFrames = uVal; break; + } + else if (name == "displayscale" && isUint && uVal > 0) { + privParams.displayScale = uVal; break; + } + else if (name == "agentsize" && isFloat && dVal > 0.0) { + privParams.agentSize = dVal; break; + } + else if (name == "genomeanalysisstride" && isUint && uVal > 0) { + privParams.genomeAnalysisStride = uVal; break; + } + else if (name == "genomeanalysisstride" && val == "videoStride") { + privParams.genomeAnalysisStride = privParams.videoStride; break; + } + else if (name == "displaysamplegenomes" && isUint) { + privParams.displaySampleGenomes = uVal; break; + } + else if (name == "genomecomparisonmethod" && isUint) { + privParams.genomeComparisonMethod = uVal; break; + } + else if (name == "updategraphlog" && isBool) { + privParams.updateGraphLog = bVal; break; + } + else if (name == "updategraphlogstride" && isUint && uVal > 0) { + privParams.updateGraphLogStride = uVal; break; + } + else if (name == "updategraphlogstride" && val == "videoStride") { + privParams.updateGraphLogStride = privParams.videoStride; break; + } + else { + std::cout << "Invalid param: " << name << " = " << val << std::endl; + } + } while (0); +} + + +void ParamManager::updateFromConfigFile() +{ + // std::ifstream is RAII, i.e. no need to call close + std::ifstream cFile(configFilename.c_str()); + if (cFile.is_open()) { + std::string line; + while(getline(cFile, line)){ + line.erase(std::remove_if(line.begin(), line.end(), isspace), + line.end()); + if(line[0] == '#' || line.empty()) { + continue; + } + auto delimiterPos = line.find("="); + auto name = line.substr(0, delimiterPos); + std::transform(name.begin(), name.end(), name.begin(), + [](unsigned char c){ return std::tolower(c); }); + auto value0 = line.substr(delimiterPos + 1); + auto delimiterComment = value0.find("#"); + auto value = value0.substr(0, delimiterComment); + auto rawValue = value; + value.erase(std::remove_if(value.begin(), value.end(), isspace), + value.end()); + //std::cout << name << " " << value << '\n' << std::endl; + ingestParameter(name, value); + } + } + else { + std::cerr << "Couldn't open config file " << configFilename << ".\n" << std::endl; + } +} + +} // end namespace BS diff --git a/src/params.h b/src/params.h new file mode 100644 index 0000000..e589527 --- /dev/null +++ b/src/params.h @@ -0,0 +1,90 @@ +#ifndef PARAMS_H_INCLUDED +#define PARAMS_H_INCLUDED + +// Global simulator parameters + +#include + +// To add a new parameter: +// 1. Add a member to struct Params in params.h. +// 2. Add a member and its default value to privParams in ParamManager::setDefault() +// in params.cpp. +// 3. Add an else clause to ParamManager::ingestParameter() in params.cpp. +// 4. Add a line to the user's parameter file (default name biosim4.ini) + +namespace BS { + +enum class RunMode { STOP, RUN, PAUSE, ABORT }; +extern RunMode runMode; + +// A private copy of Params is initialized by ParamManager::init(), then modified by +// UI events by ParamManager::uiMonitor(). The main simulator thread can get an +// updated, read-only, stable snapshot of Params with ParamManager::paramsSnapshot. + +struct Params { + unsigned population; // >= 0 + unsigned stepsPerGeneration; // > 0 + unsigned maxGenerations; // >= 0 + unsigned numThreads; // > 0 + unsigned signalLayers; // >= 0 + unsigned genomeMaxLength; // > 0 + unsigned maxNumberNeurons; // > 0 + double pointMutationRate; // 0.0..1.0 + double geneInsertionDeletionRate; // 0.0..1.0 + double deletionRatio; // 0.0..1.0 + bool killEnable; + bool sexualReproduction; + bool chooseParentsByFitness; + float populationSensorRadius; // > 0.0 + unsigned signalSensorRadius; // > 0 + float responsiveness; // >= 0.0 + unsigned responsivenessCurveKFactor; // 1, 2, 3, or 4 + unsigned longProbeDistance; // > 0 + unsigned shortProbeBarrierDistance; // > 0 + float valenceSaturationMag; + bool saveVideo; + unsigned videoStride; // > 0 + unsigned videoSaveFirstFrames; // >= 0, overrides videoStride + unsigned displayScale; + unsigned agentSize; + unsigned genomeAnalysisStride; // > 0 + unsigned displaySampleGenomes; // >= 0 + unsigned genomeComparisonMethod; // 0 = Jaro-Winkler; 1 = Hamming + bool updateGraphLog; + unsigned updateGraphLogStride; // > 0 + unsigned challenge; + unsigned barrierType; // >= 0 + unsigned replaceBarrierType; // >= 0 + unsigned replaceBarrierTypeGenerationNumber; // >= 0 + + // These must not change after initialization + uint16_t sizeX; // 2..0x10000 + uint16_t sizeY; // 2..0x10000 + unsigned genomeInitialLengthMin; // > 0 and < genomeInitialLengthMax + unsigned genomeInitialLengthMax; // > 0 and < genomeInitialLengthMin + std::string logDir; + std::string imageDir; + std::string graphLogUpdateCommand; +}; + +class ParamManager { +public: + const Params &getParamRef() const { return privParams; } // for public read-only access + void setDefaults(); + void registerConfigFile(const char *filename); + void updateFromConfigFile(); +private: + Params privParams; + std::string configFilename; + time_t lastModTime; // when config file was last read + void ingestParameter(std::string name, std::string val); +}; + +// Returns a copy of params with default values overridden by the values +// in the specified config file. The filename of the config file is saved +// inside the params for future reference. +Params paramsInit(int argc, char **argv); + +} // end namespace BS + +#endif // PARAMS_H_INCLUDED diff --git a/src/peeps.cpp b/src/peeps.cpp new file mode 100644 index 0000000..6c2b210 --- /dev/null +++ b/src/peeps.cpp @@ -0,0 +1,82 @@ +// peeps.cpp +// Manages a container of individual agents of type Indiv and their +// locations in the grid container + +#include +#include +#include +#include +#include "simulator.h" + +namespace BS { + + +Peeps::Peeps() +{ +} + + +void Peeps::init(unsigned population) +{ + // Index 0 is reserved, so add one: + individuals.resize(population + 1); +} + + +// Safe to call during multithread mode. +// Indiv will remain alive and in-world until end of sim step when +// drainDeathQueue() is called. +void Peeps::queueForDeath(const Indiv &indiv) +{ + #pragma omp critical + { + deathQueue.push_back(indiv.index); + } +} + + +// Called in single-thread mode at end of sim step. This executes all the +// queued deaths, removing the dead agents from the grid. +void Peeps::drainDeathQueue() +{ + for (uint16_t index : deathQueue) { + auto & indiv = peeps[index]; + grid.set(indiv.loc, 0); + indiv.alive = false; + } + deathQueue.clear(); +} + + +// Safe to call during multithread mode. Indiv won't move until end +// of sim step when drainMoveQueue() is called. +void Peeps::queueForMove(const Indiv &indiv, Coord newLoc) +{ + #pragma omp critical + { + auto record = std::make_pair(uint16_t(indiv.index), Coord(newLoc)); + moveQueue.push_back(record); + } +} + + +// Called in single-thread mode at end of sim step. This executes all the +// queued movements. Each movement is typically one 8-neighbor cell distance +// but this function can move an individual any arbitrary distance. +void Peeps::drainMoveQueue() +{ + for (auto moveRecord : moveQueue) { + auto & indiv = peeps[moveRecord.first]; + Coord newLoc = moveRecord.second; + Dir moveDir = (newLoc - indiv.loc).asDir(); + if (grid.isEmptyAt(newLoc)) { + grid.set(indiv.loc, 0); + grid.set(newLoc, indiv.index); + indiv.loc = newLoc; + indiv.lastMoveDir = moveDir; + } + } + moveQueue.clear(); +} + +} // end namespace BS diff --git a/src/peeps.h b/src/peeps.h new file mode 100644 index 0000000..b1b859e --- /dev/null +++ b/src/peeps.h @@ -0,0 +1,54 @@ +#ifndef PEEPS_H_INCLUDED +#define PEEPS_H_INCLUDED + +#include +#include +#include "basicTypes.h" +#include "grid.h" +#include "params.h" +#include "indiv.h" + +namespace BS { + +class Indiv; +extern class Grid grid; + +// This class keeps track of alive and dead Indiv's and where they +// are in the Grid. +// Peeps allows spawning a live Indiv at a random or specific location +// in the grid, moving Indiv's from one grid location to another, and +// killing any Indiv. +// All the Indiv instances, living and dead, are stored in the private +// .individuals member. The .cull() function will remove dead members and +// replace their slots in the .individuals container with living members +// from the end of the container for compacting the container. +// Each Indiv has an identifying index in the range 1..0xfffe that is +// stored in the Grid at the location where the Indiv resides, such that +// a Grid element value n refers to .individuals[n]. Index value 0 is +// reserved, i.e., .individuals[0] is not a valid individual. +// This class does not manage properties inside Indiv except for the +// Indiv's location in the grid and its aliveness. +class Peeps { +public: + Peeps(); // makes zero individuals + void init(unsigned population); + void queueForDeath(const Indiv &); + void drainDeathQueue(); + void queueForMove(const Indiv &, Coord newLoc); + void drainMoveQueue(); + unsigned deathQueueSize() const { return deathQueue.size(); } + // getIndiv() does no error checking -- check first that loc is occupied + Indiv & getIndiv(Coord loc) { return individuals[grid.at(loc)]; } + const Indiv & getIndiv(Coord loc) const { return individuals[grid.at(loc)]; } + // Direct access: + Indiv & operator[](uint16_t index) { return individuals[index]; } + Indiv const & operator[](uint16_t index) const { return individuals[index]; } +private: + std::vector individuals; // Index value 0 is reserved + std::vector deathQueue; + std::vector> moveQueue; +}; + +} // end namespace BS + +#endif // PEEPS_H_INCLUDED diff --git a/src/random.cpp b/src/random.cpp new file mode 100644 index 0000000..3b31fab --- /dev/null +++ b/src/random.cpp @@ -0,0 +1,94 @@ +// random.cpp +// Provides a random number generator for the main thread +// and child threads + +#include +#include +#include +#include +#include +#include "random.h" +#include "omp.h" + +namespace BS { + + +// Default is determinstic +RandomUintGenerator::RandomUintGenerator(bool deterministic) +{ + if (deterministic) { + // for Marsaglia + rngx = 123456789; + rngy = 362436000; + rngz = 521288629; + rngc = 7654321; + + // for Jenkins: + a = 0xf1ea5eed, b = c = d = 123456789; + } else { + randomize(); + } +} + + +void RandomUintGenerator::randomize() +{ + std::mt19937 generator(time(0)); // mt19937 is a standard mersenne_twister_engine + + // for Marsaglia + do { rngx = generator(); } while (rngx == 0); + do { rngy = generator(); } while (rngy == 0); + do { rngz = generator(); } while (rngz == 0); + do { rngc = generator(); } while (rngc == 0); + + // for Jenkins: + a = 0xf1ea5eed, b = c = d = generator(); +} + + +// This algorithm is from http://www0.cs.ucl.ac.uk/staff/d.jones/GoodPracticeRNG.pdf +// where it is attributed to G. Marsaglia. +// +uint32_t RandomUintGenerator::operator()() +{ + if (false) { + // Marsaglia + uint64_t t, a = 698769069ULL; + rngx = 69069 * rngx + 12345; + rngy ^= (rngy << 13); + rngy ^= (rngy >> 17); + rngy ^= (rngy << 5); /* y must never be set to zero! */ + t = a * rngz + rngc; + rngc = (t >> 32);/* Also avoid setting z=c=0! */ + return rngx + rngy + (rngz = t); + } else { + // Jenkins + #define rot32(x,k) (((x)<<(k))|((x)>>(32-(k)))) + uint32_t e = a - rot32(b, 27); + a = b ^ rot32(c, 17); + b = c + d; + c = d + e; + d = e + a; + return d; + } +} + + +// Sure, there's a bias when using modulus operator where (max - min) is not +// a power of two, but we don't care if we generate one value a little more +// often than another. Our randomness does not have to be any better quality +// than the randomness of a shotgun. We do care about speed, because this will +// get called inside deeply nested inner loops. +// +unsigned RandomUintGenerator::operator()(unsigned min, unsigned max) +{ + assert(max >= min); + return ((*this)() % (max - min + 1)) + min; +} + + +// The globally accessible random number generator. Threads can be +// given their own private copies of this. +RandomUintGenerator randomUint; + +} // end namespace BS diff --git a/src/random.h b/src/random.h new file mode 100644 index 0000000..b328c49 --- /dev/null +++ b/src/random.h @@ -0,0 +1,39 @@ +#ifndef RANDOM_H_INCLUDED +#define RANDOM_H_INCLUDED + +// See random.cpp for notes. + +#include +#include + +namespace BS { + +extern int randomInt(int min = 0, int max = INT_MAX); +extern uint16_t randomU16(uint16_t min = 0, uint16_t max = (uint16_t)-1); +extern float randomFloat01(); // 0.0..1.0 +extern float randomFloat11(); // -1.0..1.0 + +class RandomUintGenerator{ +private: + // for Marsaglia + uint32_t rngx; + uint32_t rngy; + uint32_t rngz; + uint32_t rngc; + // for Jenkins + uint32_t a, b, c, d; +public: + RandomUintGenerator(bool deterministic = false); + RandomUintGenerator& operator=(const RandomUintGenerator &rhs) = default; + void randomize(); + uint32_t operator()(); + unsigned operator()(unsigned min, unsigned max); +}; + +// The globally accessible random number generator (not thread safe) +extern RandomUintGenerator randomUint; +constexpr uint32_t RANDOM_UINT_MAX = 0xffffffff; + +} // end namespace BS + +#endif // RANDOM_H_INCLUDED diff --git a/src/sensors-actions.h b/src/sensors-actions.h new file mode 100644 index 0000000..30b768a --- /dev/null +++ b/src/sensors-actions.h @@ -0,0 +1,95 @@ +#ifndef SENSORS_ACTIONS_H_INCLUDED +#define SENSORS_ACTIONS_H_INCLUDED + +// This file defines which sensor input neurons and which action output neurons +// are compiled into the simulator. This file can be modified to create a simulator +// executable that supports only a subset of all possible sensor or action neurons. + +#include + +namespace BS { + +// Neuron Sources (Sensors) and Sinks (Actions) + +// These sensor, neuron, and action value ranges are here for documentation +// purposes. Most functions now assume these ranges. We no longer support changes +// to these ranges. +constexpr float SENSOR_MIN = 0.0; +constexpr float SENSOR_MAX = 1.0; +constexpr float SENSOR_RANGE = SENSOR_MAX - SENSOR_MIN; + +constexpr float NEURON_MIN = -1.0; +constexpr float NEURON_MAX = 1.0; +constexpr float NEURON_RANGE = NEURON_MAX - NEURON_MIN; + +constexpr float ACTION_MIN = 0.0; +constexpr float ACTION_MAX = 1.0; +constexpr float ACTION_RANGE = ACTION_MAX - ACTION_MIN; + + +// Place the sensor neuron you want enabled prior to NUM_SENSES. Any +// that are after NUM_SENSES will be disabled in the simulator. +// If new items are added to this enum, also update the name functions +// in analysis.cpp. +// I means data about the individual, mainly stored in Indiv +// W means data about the environment, mainly stored in Peeps or Grid +enum Sensor { + LOC_X, // I distance from left edge + LOC_Y, // I distance from bottom + BOUNDARY_DIST_X, // I X distance to nearest edge of world + BOUNDARY_DIST, // I distance to nearest edge of world + BOUNDARY_DIST_Y, // I Y distance to nearest edge of world + GENETIC_SIM_FWD, // I genetic similarity forward + LAST_MOVE_DIR_X, // I +- amount of X movement in last movement + LAST_MOVE_DIR_Y, // I +- amount of Y movement in last movement + LONGPROBE_POP_FWD, // W long look for population forward + LONGPROBE_BAR_FWD, // W long look for barriers forward + POPULATION, // W population density in neighborhood + POPULATION_FWD, // W population density in the forward-reverse axis + POPULATION_LR, // W population density in the left-right axis + OSC1, // I oscillator +-value + AGE, // I + BARRIER_FWD, // W neighborhood barrier distance forward-reverse axis + BARRIER_LR, // W neighborhood barrier distance left-right axis + RANDOM, // random sensor value, uniform distribution + SIGNAL0, // W strength of signal0 in neighborhood + SIGNAL0_FWD, // W strength of signal0 in the forward-reverse axis + SIGNAL0_LR, // W strength of signal0 in the left-right axis + NUM_SENSES, // <<------------------ END OF ACTIVE SENSES MARKER +}; + + +// Place the action neuron you want enabled prior to NUM_ACTIONS. Any +// that are after NUM_ACTIONS will be disabled in the simulator. +// If new items are added to this enum, also update the name functions +// in analysis.cpp. +// I means the action affects the individual internally (Indiv) +// W means the action also affects the environment (Peeps or Grid) +enum Action { + MOVE_X, // W +- X component of movement + MOVE_Y, // W +- Y component of movement + MOVE_FORWARD, // W continue last direction + MOVE_RL, // W +- component of movement + MOVE_RANDOM, // W + SET_OSCILLATOR_PERIOD, // I + SET_LONGPROBE_DIST, // I + SET_RESPONSIVENESS, // I + EMIT_SIGNAL0, // W + MOVE_EAST, // W + MOVE_WEST, // W + MOVE_NORTH, // W + MOVE_SOUTH, // W + MOVE_LEFT, // W + MOVE_RIGHT, // W + MOVE_REVERSE, // W + NUM_ACTIONS, // <<----------------- END OF ACTIVE ACTIONS MARKER + KILL_FORWARD, // W +}; + +extern std::string sensorName(Sensor sensor); +extern std::string actionName(Action action); +extern void printSensorsActions(); // list the names to stdout + +} // end namespace BS + +#endif // SENSORS_ACTIONS_H_INCLUDED diff --git a/src/signals.cpp b/src/signals.cpp new file mode 100644 index 0000000..9de92bc --- /dev/null +++ b/src/signals.cpp @@ -0,0 +1,61 @@ +// signals.cpp +// Manages layers of pheremones + +#include +#include "simulator.h" + +namespace BS { + +void Signals::init(uint16_t numLayers, uint16_t sizeX, uint16_t sizeY) +{ + data = std::vector(numLayers, Layer(sizeX, sizeY)); +} + + +// Increases the specified location by centerIncreaseAmount, +// and increases the neighboring cells by neighborIncreaseAmount + +// Is it ok that multiple readers are reading this container while +// this single thread is writing to it? todo!!! +void Signals::increment(uint16_t layerNum, Coord loc) +{ + constexpr float radius = 1.5; + constexpr uint8_t centerIncreaseAmount = 2; + constexpr uint8_t neighborIncreaseAmount = 1; + +#pragma omp critical + { + visitNeighborhood(loc, radius, [layerNum](Coord loc) { + if (signals[layerNum][loc.x][loc.y] < SIGNAL_MAX) { + signals[layerNum][loc.x][loc.y] = + std::min(SIGNAL_MAX, + signals[layerNum][loc.x][loc.y] + neighborIncreaseAmount); + } + }); + + if (signals[layerNum][loc.x][loc.y] < SIGNAL_MAX) { + signals[layerNum][loc.x][loc.y] = + std::min(SIGNAL_MAX, + signals[layerNum][loc.x][loc.y] + centerIncreaseAmount); + } + } +} + + +// Fades the signals +void Signals::fade(unsigned layerNum) +{ + constexpr unsigned fadeAmount = 1; + + for (int16_t x = 0; x < p.sizeX; ++x) { + for (int16_t y = 0; y < p.sizeY; ++y) { + if (signals[layerNum][x][y] >= fadeAmount) { + signals[layerNum][x][y] -= fadeAmount; // fade center cell + } else { + signals[layerNum][x][y] = 0; + } + } + } +} + +} // end namespace BS diff --git a/src/signals.h b/src/signals.h new file mode 100644 index 0000000..a096f50 --- /dev/null +++ b/src/signals.h @@ -0,0 +1,51 @@ +#ifndef SIGNALS_H_INCLUDED +#define SIGNALS_H_INCLUDED + +// Container for pheromones. + +#include +#include +#include "basicTypes.h" + +namespace BS { + +// Usage: uint8_t magnitude = signals[layer][x][y]; +// or magnitude = signals.getMagnitude(layer, Coord); + + +constexpr unsigned SIGNAL_MIN = 0; +constexpr unsigned SIGNAL_MAX = UINT8_MAX; + +struct Signals { + struct Column { + Column(uint16_t numRows) : data { std::vector(numRows, 0) } { } + uint8_t& operator[](uint16_t rowNum) { return data[rowNum]; } + uint8_t operator[](uint16_t rowNum) const { return data[rowNum]; } + void zeroFill() { std::fill(data.begin(), data.end(), 0); } + private: + std::vector data; + }; + + struct Layer { + Layer(uint16_t numCols, uint16_t numRows) : data { std::vector(numCols, Column(numRows)) } { } + Column& operator[](uint16_t colNum) { return data[colNum]; } + const Column& operator[](uint16_t colNum) const { return data[colNum]; } + void zeroFill() { for (Column &col : data) { col.zeroFill(); } } + private: + std::vector data; + }; + + void init(uint16_t layers, uint16_t sizeX, uint16_t sizeY); + Layer& operator[](uint16_t layerNum) { return data[layerNum]; } + const Layer& operator[](uint16_t layerNum) const { return data[layerNum]; } + uint8_t getMagnitude(uint16_t layerNum, Coord loc) const { return (*this)[layerNum][loc.x][loc.y]; } + void increment(uint16_t layerNum, Coord loc); + void zeroFill() { for (Layer &layer : data) { layer.zeroFill(); } } + void fade(unsigned layerNum); +private: + std::vector data; +}; + +} // end namespace BS + +#endif // SIGNALS_H_INCLUDED diff --git a/src/simulator.cpp b/src/simulator.cpp new file mode 100644 index 0000000..dfe8bf7 --- /dev/null +++ b/src/simulator.cpp @@ -0,0 +1,171 @@ +// simulator.cpp - Main thread +// This file contains simulator(), the top-level entry point of the simulator. +// simulator() is called from main.cpp with a copy of argc and argv. +// If there is no command line argument, the simulator will read the default +// config file ("biosim4.ini" in the current directory) to get the simulation +// parameters for this run. If there are one or more command line args, then +// argv[1] must contain the name of the config file which will be read instead +// of biosim4.ini. Any args after that are ignored. The simulator code is +// in namespace BS (for "biosim"). + +#include +#include +#include +#include +#include +#include "simulator.h" // the simulator data structures +#include "imageWriter.h" // this is for generating the movies + +namespace BS { + +extern void initializeGeneration0(); +extern unsigned spawnNewGeneration(unsigned generation, unsigned murderCount); +extern void displaySampleGenomes(unsigned count); +extern void executeActions(Indiv &indiv, std::array &actionLevels); +extern void endOfSimStep(unsigned simStep, unsigned generation); +extern void endOfGeneration(unsigned generation); + +RunMode runMode = RunMode::STOP; +Grid grid; // The 2D world where the creatures live +Signals signals; // A 2D array of pheromones that overlay the world grid +Peeps peeps; // The container of all the individuals in the population +ImageWriter imageWriter; // This is for generating the movies + +// The paramManger maintains a private copy of the parameter values, and a copy +// is available read-only through global variable p. Although this is not +// foolproof, you should be able to modify the config file during a simulation +// run and modify many of the parameters. See params.cpp and params.h for more info. +ParamManager paramManager; +const Params &p { paramManager.getParamRef() }; // read-only params + + +/********************************************************************************************** +Execute one simStep for one individual. + +This executes in its own thread, invoked from the main simulator thread. First we execute +indiv.feedForward() which computes action values to be executed here. Some actions such as +signal emission(s) (pheromones), agent movement, or deaths will have been queued for +later execution at the end of the generation in single-threaded mode (the deferred queues +allow the main data structures (e.g., grid, signals) to be freely accessed read-only in all threads). + +In order to be thread-safe, the main simulator-wide data structures and their +accessibility are: + + grid - read-only + signals - (pheromones) read-write for the location where our agent lives + using signals.increment(), read-only for other locations + peeps - for other individuals, we can only read their index and genome. + We have read-write access to our individual through the indiv argument. + +The other important variables are: + + simStep - the current age of our agent, reset to 0 at the start of each generation. + For many simulation scenarios, this matches our indiv.age member. + randomUint - global random number generator, a private instance is given to each thread +**********************************************************************************************/ +void simStepOneIndiv(Indiv &indiv, unsigned simStep) +{ + ++indiv.age; // for this implementation, tracks simStep + auto actionLevels = indiv.feedForward(simStep); + executeActions(indiv, actionLevels); +} + + +/******************************************************************************** +Start of simulator + +All the agents are randomly placed with random genomes at the start. The outer +loop is generation, the inner loop is simStep. There is a fixed number of +simSteps in each generation. Agents can die at any simStep and their corpses +remain until the end of the generation. At the end of the generation, the +dead corpses are removed, the survivors reproduce and then die. The newborns +are placed at random locations, signals (pheromones) are updated, simStep is +reset to 0, and a new generation proceeds. + +The paramManager manages all the simulator parameters. It starts with defaults, +then keeps them updated as the config file (biosim4.ini) changes. + +The main simulator-wide data structures are: + grid - where the agents live (identified by their non-zero index). 0 means empty. + signals - multiple layers overlay the grid, hold pheromones + peeps - an indexed set of agents of type Indiv; indexes start at 1 + +The important simulator-wide variables are: + generation - starts at 0, then increments every time the agents die and reproduce. + simStep - reset to 0 at the start of each generation; fixed number per generation. + randomUint - global random number generator + +The threads are: + main thread - simulator + simStepOneIndiv() - child threads created by the main simulator thread + imageWriter - saves image frames used to make a movie (possibly not threaded + due to unresolved bugs when threaded) +********************************************************************************/ +void simulator(int argc, char **argv) +{ + printSensorsActions(); // show the agents' capabilities + + // Simulator parameters are available read-only through the global + // variable p after paramManager is initialized. + // Todo: remove the hardcoded parameter filename. + paramManager.setDefaults(); + paramManager.registerConfigFile(argc > 1 ? argv[1] : "biosim4.ini"); + paramManager.updateFromConfigFile(); + + // Allocate container space. Once allocated, these container elements + // will be reused in each new generation. + grid.init(p.sizeX, p.sizeY); // the land on which the peeps live + signals.init(p.signalLayers, p.sizeX, p.sizeY); // where the pheromones waft + peeps.init(p.population); // the peeps themselves + + // If imageWriter is to be run in its own thread, start it here: + //std::thread t(&ImageWriter::saveFrameThread, &imageWriter); + + // Unit tests: + //unitTestConnectNeuralNetWiringFromGenome(); + //unitTestGridVisitNeighborhood(); + + unsigned generation = 0; + initializeGeneration0(); // starting population + + runMode = RunMode::RUN; + while (runMode == RunMode::RUN && generation < p.maxGenerations) { // generation loop + unsigned murderCount = 0; // for reporting purposes + for (unsigned simStep = 0; simStep < p.stepsPerGeneration; ++simStep) { + + // multithreaded loop: index 0 is reserved, start at 1 +#pragma omp parallel for num_threads(p.numThreads) default(shared) firstprivate(randomUint) lastprivate(randomUint) schedule(auto) + for (unsigned indivIndex = 1; indivIndex <= p.population; ++indivIndex) { + if (peeps[indivIndex].alive) { + simStepOneIndiv(peeps[indivIndex], simStep); + } + } + // In single-thread mode: this executes deferred, queued deaths and movements, + // updates signal layers (pheromone), etc. + murderCount += peeps.deathQueueSize(); + endOfSimStep(simStep, generation); + } + + endOfGeneration(generation); + paramManager.updateFromConfigFile(); + unsigned numberSurvivors = spawnNewGeneration(generation, murderCount); + if (numberSurvivors > 0 && (generation % p.genomeAnalysisStride == 0)) { + displaySampleGenomes(p.displaySampleGenomes); + } + if (numberSurvivors == 0) { + generation = 0; // start over + } else { + ++generation; + } + } + displaySampleGenomes(3); // final report, for debugging + + std::cout << "Simulator waiting..." << std::endl; + std::cout << "Simulator exit." << std::endl; + + // If imageWriter is in its own thread, stop it and wait for it here: + //imageWriter.abort(); + //t.join(); +} + +} // end namespace BS diff --git a/src/simulator.h b/src/simulator.h new file mode 100644 index 0000000..8da51e4 --- /dev/null +++ b/src/simulator.h @@ -0,0 +1,53 @@ +#ifndef SIMULATOR_H_INCLUDED +#define SIMULATOR_H_INCLUDED + +// Main header for the simulator. Also see simulator.cpp. + +#include "basicTypes.h" // types Dir, Coord, Polar and their values +#include "params.h" // params from the config file plus some extra stuff +#include "indiv.h" // data structure for an individual +#include "grid.h" // the 2D world where the peeps live +#include "signals.h" // a 2D array of pheromones that overlay the world grid +#include "peeps.h" // the 2D world where the peeps live +#include "random.h" + +namespace BS { + +// Some of the survival challenges to try. Some are interesting, some +// not so much. Fine-tune the challenges by tweaking the corresponding code +// in survival-criteria.cpp. +constexpr unsigned CHALLENGE_CIRCLE = 0; +constexpr unsigned CHALLENGE_RIGHT_HALF = 1; +constexpr unsigned CHALLENGE_RIGHT_QUARTER = 2; +constexpr unsigned CHALLENGE_STRING = 3; +constexpr unsigned CHALLENGE_CENTER_WEIGHTED = 4; +constexpr unsigned CHALLENGE_CENTER_UNWEIGHTED = 40; +constexpr unsigned CHALLENGE_CORNER = 5; +constexpr unsigned CHALLENGE_CORNER_WEIGHTED = 6; +constexpr unsigned CHALLENGE_MIGRATE_DISTANCE = 7; +constexpr unsigned CHALLENGE_CENTER_SPARSE = 8; +constexpr unsigned CHALLENGE_LEFT_EIGHTH = 9; +constexpr unsigned CHALLENGE_RADIOACTIVE_WALLS = 10; +constexpr unsigned CHALLENGE_AGAINST_ANY_WALL = 11; +constexpr unsigned CHALLENGE_TOUCH_ANY_WALL = 12; +constexpr unsigned CHALLENGE_EAST_WEST_EIGHTHS = 13; +constexpr unsigned CHALLENGE_NEAR_BARRIER = 14; +constexpr unsigned CHALLENGE_PAIRS = 15; +constexpr unsigned CHALLENGE_LOCATION_SEQUENCE = 16; +constexpr unsigned CHALLENGE_ALTRUISM = 17; +constexpr unsigned CHALLENGE_ALTRUISM_SACRIFICE = 18; + +extern ParamManager paramManager; // manages simulator params from the config file plus more +extern const Params &p; // read-only simulator config params +extern Grid grid; // 2D arena where the individuals live +extern Signals signals; // pheromone layers +extern Peeps peeps; // container of all the individuals +extern void simulator(int argc, char **argv); + +// Feeds in-bounds Coords to a function: given a center location and a radius, this +// function will call f(Coord) once for each location inside the specified area. +extern void visitNeighborhood(Coord loc, float radius, std::function f); + +} // end namespace BS + +#endif // SIMULATOR_H_INCLUDED diff --git a/src/spawnNewGeneration.cpp b/src/spawnNewGeneration.cpp new file mode 100644 index 0000000..647395c --- /dev/null +++ b/src/spawnNewGeneration.cpp @@ -0,0 +1,195 @@ +// spawnNewGeneration.cpp + +#include +#include +#include +#include +#include +#include "simulator.h" + +namespace BS { + +extern std::pair passedSurvivalCriterion(const Indiv &indiv, unsigned challenge); + + +// Requires that the grid, signals, and peeps containers have been allocated. +// This will erase the grid and signal layers, then create a new population in +// the peeps container at random locations with random genomes. +void initializeGeneration0() +{ + // The grid has already been allocated, just clear and reuse it + grid.zeroFill(); + grid.createBarrier(p.replaceBarrierTypeGenerationNumber == 0 + ? p.replaceBarrierType : p.barrierType); + + // The signal layers have already been allocated, so just reuse them + signals.zeroFill(); + + // Spawn the population. The peeps container has already been allocated, + // just clear and reuse it + for (uint16_t index = 1; index <= p.population; ++index) { + peeps[index].initialize(index, grid.findEmptyLocation(), makeRandomGenome()); + } +} + + +// Requires a container with one or more parent genomes to choose from. +// Called from spawnNewGeneration(). This requires that the grid, signals, and +// peeps containers have been allocated. This will erase the grid and signal +// layers, then create a new population in the peeps container with random +// locations and genomes derived from the container of parent genomes. +void initializeNewGeneration(const std::vector &parentGenomes, unsigned generation) +{ + extern Genome generateChildGenome(const std::vector &parentGenomes); + + // The grid, signals, and peeps containers have already been allocated, just + // clear them if needed and reuse the elements + grid.zeroFill(); + grid.createBarrier(generation >= p.replaceBarrierTypeGenerationNumber + ? p.replaceBarrierType : p.barrierType); + signals.zeroFill(); + + // Spawn the population. This overwrites all the elements of peeps[] + for (uint16_t index = 1; index <= p.population; ++index) { + peeps[index].initialize(index, grid.findEmptyLocation(), generateChildGenome(parentGenomes)); + } +} + + +// At this point, the deferred death queue and move queue have been processed +// and we are left with zero or more individuals who will repopulate the +// world grid. +// In order to redistribute the new population randomly, we will save all the +// surviving genomes in a container, then clear the grid of indexes and generate +// new individuals. This is inefficient when there are lots of survivors because +// we could have reused (with mutations) the survivors' genomes and neural +// nets instead of rebuilding them. +// Returns number of survivor-reproducers. +// Must be called in single-thread mode between generations. +unsigned spawnNewGeneration(unsigned generation, unsigned murderCount) +{ + unsigned sacrificedCount = 0; // for the altruism challenge + + extern void appendEpochLog(unsigned generation, unsigned numberSurvivors, unsigned murderCount); + extern std::pair passedSurvivalCriterion(const Indiv &indiv, unsigned challenge); + extern void displaySignalUse(); + + // This container will hold the indexes and survival scores (0.0..1.0) + // of all the survivors who will provide genomes for repopulation. + std::vector> parents; // + + // This container will hold the genomes of the survivors + std::vector parentGenomes; + + if (p.challenge != CHALLENGE_ALTRUISM) { + // First, make a list of all the individuals who will become parents; save + // their scores for later sorting. Indexes start at 1. + for (uint16_t index = 1; index <= p.population; ++index) { + std::pair passed = passedSurvivalCriterion(peeps[index], p.challenge); + // Save the parent genome if it results in valid neural connections + // ToDo: if the parents no longer need their genome record, we could + // possibly do a move here instead of copy, although it's doubtful that + // the optimization would be noticeable. + if (passed.first && peeps[index].nnet.connections.size() > 0) { + parents.push_back( { index, passed.second } ); + } + } + } else { + // For the altruism challenge, test if the agent is inside either the sacrificial + // or the spawning area. We'll count the number in the sacrificial area and + // save the genomes of the ones in the spawning area, saving their scores + // for later sorting. Indexes start at 1. + + bool considerKinship = true; + std::vector sacrificesIndexes; // those who gave their lives for the greater good + + for (uint16_t index = 1; index <= p.population; ++index) { + // This the test for the spawning area: + std::pair passed = passedSurvivalCriterion(peeps[index], CHALLENGE_ALTRUISM); + if (passed.first && peeps[index].nnet.connections.size() > 0) { + parents.push_back( { index, passed.second } ); + } else { + // This is the test for the sacrificial area: + passed = passedSurvivalCriterion(peeps[index], CHALLENGE_ALTRUISM_SACRIFICE); + if (passed.first && peeps[index].nnet.connections.size() > 0) { + if (considerKinship) { + sacrificesIndexes.push_back(index); + } else { + ++sacrificedCount; + } + } + } + } + + unsigned generationToApplyKinship = 10; + constexpr unsigned altruismFactor = 10; // the saved:sacrificed ratio + + if (considerKinship) { + if (generation > generationToApplyKinship) { + // Todo: optimize!!! + float threshold = 0.7; + + std::vector> survivingKin; + for (unsigned passes = 0; passes < altruismFactor; ++passes) { + for (uint16_t sacrificedIndex : sacrificesIndexes) { + // randomize the next loop so we don't keep using the first one repeatedly + unsigned startIndex = randomUint(0, parents.size() - 1); + for (unsigned count = 0; count < parents.size(); ++count) { + const std::pair &possibleParent = parents[(startIndex + count) % parents.size()]; + const Genome &g1 = peeps[sacrificedIndex].genome; + const Genome &g2 = peeps[possibleParent.first].genome; + float similarity = genomeSimilarity(g1, g2); + if (similarity >= threshold) { + survivingKin.push_back(possibleParent); + // mark this one so we don't use it again? + break; + } + } + } + } + std::cout << parents.size() << " passed, " + << sacrificesIndexes.size() << " sacrificed, " + << survivingKin.size() << " saved" << std::endl; // !!! + parents = std::move(survivingKin); + } + } else { + // Limit the parent list + unsigned numberSaved = sacrificedCount * altruismFactor; + std::cout << parents.size() << " passed, " << sacrificedCount << " sacrificed, " << numberSaved << " saved" << std::endl; // !!! + if (parents.size() > 0 && numberSaved < parents.size()) { + parents.erase(parents.begin() + numberSaved, parents.end()); + } + } + } + + // Sort the indexes of the parents by their fitness scores + std::sort(parents.begin(), parents.end(), + [](const std::pair &parent1, const std::pair &parent2) { + return parent1.second > parent2.second; + }); + + // Assemble a list of all the parent genomes. These will be ordered by their + // scores if the parents[] container was sorted by score + for (const std::pair &parent : parents) { + parentGenomes.push_back(peeps[parent.first].genome); + } + + std::cout << "Gen " << generation << ", " << parentGenomes.size() << " survivors" << std::endl; + appendEpochLog(generation, parentGenomes.size(), murderCount); + //displaySignalUse(); // for debugging only + + // Now we have a container of zero or more parents' genomes + + if (parentGenomes.size() != 0) { + // Spawn a new generation + initializeNewGeneration(parentGenomes, generation + 1); + } else { + // Special case: there are no surviving parents: start the simulation over + // from scratch with randomly-generated genomes + initializeGeneration0(); + } + + return parentGenomes.size(); +} + +} // end namespace BS diff --git a/src/survival-criteria.cpp b/src/survival-criteria.cpp new file mode 100644 index 0000000..c6c7001 --- /dev/null +++ b/src/survival-criteria.cpp @@ -0,0 +1,345 @@ +// survival-criteria.cpp + +#include +#include +#include "simulator.h" + +namespace BS { + +// Returns true and a score 0.0..1.0 if passed, false if failed +std::pair passedSurvivalCriterion(const Indiv &indiv, unsigned challenge) +{ + if (!indiv.alive) { + return { false, 0.0 }; + } + + switch(challenge) { + + // Survivors are those inside the circular area defined by + // safeCenter and radius + case CHALLENGE_CIRCLE: + { + Coord safeCenter { (int16_t)(p.sizeX / 4.0), (int16_t)(p.sizeY / 4.0) }; + float radius = p.sizeX / 4.0; + + Coord offset = safeCenter - indiv.loc; + float distance = offset.length(); + return distance <= radius ? + std::pair { true, (radius - distance) / radius } + : std::pair { false, 0.0 }; + } + + // Survivors are all those on the right side of the arena + case CHALLENGE_RIGHT_HALF: + return indiv.loc.x > p.sizeX / 2 ? + std::pair { true, 1.0 } + : std::pair { false, 0.0 }; + + // Survivors are all those on the right quarter of the arena + case CHALLENGE_RIGHT_QUARTER: + return indiv.loc.x > p.sizeX / 2 + p.sizeX / 4 ? + std::pair { true, 1.0 } + : std::pair { false, 0.0 }; + + // Survivors are all those on the left eighth of the arena + case CHALLENGE_LEFT_EIGHTH: + return indiv.loc.x < p.sizeX / 8 ? + std::pair { true, 1.0 } + : std::pair { false, 0.0 }; + + // Survivors are those not touching the border and with exactly the number + // of neighbors defined by neighbors and radius, where neighbors includes self + case CHALLENGE_STRING: + { + unsigned minNeighbors = 22; + unsigned maxNeighbors = 2; + float radius = 1.5; + + if (grid.isBorder(indiv.loc)) { + return { false, 0.0 }; + } + + unsigned count = 0; + auto f = [&](Coord loc2){ + if (grid.isOccupiedAt(loc2)) ++count; + }; + + visitNeighborhood(indiv.loc, radius, f); + if (count >= minNeighbors && count <= maxNeighbors) { + return { true, 1.0 }; + } else { + return { false, 0.0 }; + } + } + + // Survivors are those within the specified radius of the center. The score + // is linearly weighted by distance from the center. + case CHALLENGE_CENTER_WEIGHTED: + { + Coord safeCenter { (int16_t)(p.sizeX / 2.0), (int16_t)(p.sizeY / 2.0) }; + float radius = p.sizeX / 3.0; + + Coord offset = safeCenter - indiv.loc; + float distance = offset.length(); + return distance <= radius ? + std::pair { true, (radius - distance) / radius } + : std::pair { false, 0.0 }; + } + + // Survivors are those within the specified radius of the center + case CHALLENGE_CENTER_UNWEIGHTED: + { + Coord safeCenter { (int16_t)(p.sizeX / 2.0), (int16_t)(p.sizeY / 2.0) }; + float radius = p.sizeX / 3.0; + + Coord offset = safeCenter - indiv.loc; + float distance = offset.length(); + return distance <= radius ? + std::pair { true, 1.0 } + : std::pair { false, 0.0 }; + } + + // Survivors are those within the specified outer radius of the center and with + // the specified number of neighbors in the specified inner radius. + // The score is not weighted by distance from the center. + case CHALLENGE_CENTER_SPARSE: + { + Coord safeCenter { (int16_t)(p.sizeX / 2.0), (int16_t)(p.sizeY / 2.0) }; + float outerRadius = p.sizeX / 4.0; + float innerRadius = 1.5; + unsigned minNeighbors = 5; // includes self + unsigned maxNeighbors = 8; + + Coord offset = safeCenter - indiv.loc; + float distance = offset.length(); + if (distance <= outerRadius) { + unsigned count = 0; + auto f = [&](Coord loc2){ + if (grid.isOccupiedAt(loc2)) ++count; + }; + + visitNeighborhood(indiv.loc, innerRadius, f); + if (count >= minNeighbors && count <= maxNeighbors) { + return { true, 1.0 }; + } + } + return { false, 0.0 }; + } + + // Survivors are those within the specified radius of any corner. + // Assumes square arena. + case CHALLENGE_CORNER: + { + assert(p.sizeX == p.sizeY); + float radius = p.sizeX / 8.0; + + float distance = (Coord(0, 0) - indiv.loc).length(); + if (distance <= radius) { + return { true, 1.0 }; + } + distance = (Coord(0, p.sizeY - 1) - indiv.loc).length(); + if (distance <= radius) { + return { true, 1.0 }; + } + distance = (Coord(p.sizeX - 1, 0) - indiv.loc).length(); + if (distance <= radius) { + return { true, 1.0 }; + } + distance = (Coord(p.sizeX - 1, p.sizeY - 1) - indiv.loc).length(); + if (distance <= radius) { + return { true, 1.0 }; + } + return { false, 0.0 }; + } + + // Survivors are those within the specified radius of any corner. The score + // is linearly weighted by distance from the corner point. + case CHALLENGE_CORNER_WEIGHTED: + { + assert(p.sizeX == p.sizeY); + float radius = p.sizeX / 4.0; + + float distance = (Coord(0, 0) - indiv.loc).length(); + if (distance <= radius) { + return { true, (radius - distance) / radius }; + } + distance = (Coord(0, p.sizeY - 1) - indiv.loc).length(); + if (distance <= radius) { + return { true, (radius - distance) / radius }; + } + distance = (Coord(p.sizeX - 1, 0) - indiv.loc).length(); + if (distance <= radius) { + return { true, (radius - distance) / radius }; + } + distance = (Coord(p.sizeX - 1, p.sizeY - 1) - indiv.loc).length(); + if (distance <= radius) { + return { true, (radius - distance) / radius }; + } + return { false, 0.0 }; + } + + // This challenge is handled in endOfSimStep(), where individuals may die + // at the end of any sim step. There is nothing else to do here at the + // end of a generation. All remaining alive become parents. + case CHALLENGE_RADIOACTIVE_WALLS: + return { true, 1.0 }; + + // Survivors are those touching any wall at the end of the generation + case CHALLENGE_AGAINST_ANY_WALL: + { + bool onEdge = indiv.loc.x == 0 || indiv.loc.x == p.sizeX - 1 + || indiv.loc.y == 0 || indiv.loc.y == p.sizeY - 1; + + if (onEdge) { + return { true, 1.0 }; + } else { + return { false, 0.0 }; + } + } + + // This challenge is partially handled in endOfSimStep(), where individuals + // that are touching a wall are flagged in their Indiv record. They are + // allowed to continue living. Here at the end of the generation, any that + // never touch a wall will die. All that touched a wall at any time during + // their life will become parents. + case CHALLENGE_TOUCH_ANY_WALL: + if (indiv.challengeBits != 0) { + return { true, 1.0 }; + } else { + return { false, 0.0 }; + } + + // Everybody survives and are candidate parents, but scored by how far + // they migrated from their birth location. + case CHALLENGE_MIGRATE_DISTANCE: + { + //unsigned requiredDistance = p.sizeX / 2.0; + float distance = (indiv.loc - indiv.birthLoc).length(); + distance = distance / (float)(std::max(p.sizeX, p.sizeY)); + return { true, distance }; + } + + // Survivors are all those on the left or right eighths of the arena + case CHALLENGE_EAST_WEST_EIGHTHS: + return indiv.loc.x < p.sizeX / 8 || indiv.loc.x >= (p.sizeX - p.sizeX / 8)? + std::pair { true, 1.0 } + : std::pair { false, 0.0 }; + + // Survivors are those within radius of any barrier center. Weighted by distance. + case CHALLENGE_NEAR_BARRIER: + { + float radius; + //radius = 20.0; + radius = p.sizeX / 2; + //radius = p.sizeX / 4; + + const std::vector barrierCenters = grid.getBarrierCenters(); + float minDistance = 1e8; + for (Coord center : barrierCenters) { + float distance = (indiv.loc - center).length(); + if (distance < minDistance) { + minDistance = distance; + } + } + if (minDistance <= radius) { + return { true, 1.0 - (minDistance / radius) }; + } else { + return { false, 0.0 }; + } + } + + // Survivors are those not touching a border and with exactly one neighbor which has no other neighbor + case CHALLENGE_PAIRS: + { + bool onEdge = indiv.loc.x == 0 || indiv.loc.x == p.sizeX - 1 + || indiv.loc.y == 0 || indiv.loc.y == p.sizeY - 1; + + if (onEdge) { + return { false, 0.0 }; + } + + unsigned count = 0; + for (int16_t x = indiv.loc.x - 1; x < indiv.loc.x + 1; ++x) { + for (int16_t y = indiv.loc.y - 1; y < indiv.loc.y + 1; ++y) { + Coord tloc = { x, y }; + if (tloc != indiv.loc && grid.isInBounds(tloc) && grid.isOccupiedAt(tloc)) { + ++count; + if (count == 1) { + for (int16_t x1 = tloc.x - 1; x1 < tloc.x + 1; ++x1) { + for (int16_t y1 = tloc.y - 1; y1 < tloc.y + 1; ++y1) { + Coord tloc1 = { x1, y1 }; + if (tloc1 != tloc && tloc1 != indiv.loc && grid.isInBounds(tloc1) && grid.isOccupiedAt(tloc1)) { + return { false, 0.0 }; + } + } + } + } else { + return { false, 0.0 }; + } + } + } + } + if (count == 1) { + return { true, 1.0 }; + } else { + return { false, 0.0 }; + } + } + + // Survivors are those that contacted one or more specified locations in a sequence, + // ranked by the number of locations contacted. There will be a bit set in their + // challengeBits member for each location contacted. + case CHALLENGE_LOCATION_SEQUENCE: + { + unsigned count = 0; + unsigned bits = indiv.challengeBits; + unsigned maxNumberOfBits = sizeof(bits) * 8; + + for (unsigned n = 0; n < maxNumberOfBits; ++n) { + if ((bits & (1 << n)) != 0) { + ++count; + } + } + if (count > 0) { + return { true, count / (float)maxNumberOfBits }; + } else { + return { false, 0.0 }; + } + } + break; + + // Survivors are all those within the specified radius of the NE corner + case CHALLENGE_ALTRUISM_SACRIFICE: + { + //float radius = p.sizeX / 3.0; // in 128^2 world, holds 1429 agents + float radius = p.sizeX / 4.0; // in 128^2 world, holds 804 agents + //float radius = p.sizeX / 5.0; // in 128^2 world, holds 514 agents + + float distance = (Coord(p.sizeX - p.sizeX / 4, p.sizeY - p.sizeY / 4) - indiv.loc).length(); + if (distance <= radius) { + return { true, (radius - distance) / radius }; + } else { + return { false, 0.0 }; + } + } + + // Survivors are those inside the circular area defined by + // safeCenter and radius + case CHALLENGE_ALTRUISM: + { + Coord safeCenter { (int16_t)(p.sizeX / 4.0), (int16_t)(p.sizeY / 4.0) }; + float radius = p.sizeX / 4.0; // in a 128^2 world, holds 3216 + + Coord offset = safeCenter - indiv.loc; + float distance = offset.length(); + return distance <= radius ? + std::pair { true, (radius - distance) / radius } + : std::pair { false, 0.0 }; + } + + default: + assert(false); + } +} + +} // end namespace BS diff --git a/src/unitTestBasicTypes.cpp b/src/unitTestBasicTypes.cpp new file mode 100644 index 0000000..a419ed4 --- /dev/null +++ b/src/unitTestBasicTypes.cpp @@ -0,0 +1,312 @@ +// unitTestBasicTypes.cpp +// This tests the types Dir, Coord, and Polar, and enum Compass. +// See basicTypes.h for more info about the basic types. + +#include +#include +#include "basicTypes.h" + +namespace BS { + +// returns true only if a and b are the same to within +// 4 digits accuracy +bool areClosef(float a, float b) +{ + return std::abs(a - b) < 0.0001; +} + + +bool unitTestBasicTypes() +{ + // Dir + + // ctor from Compass + // .asInt() + // copy assignment + Dir d1 = Dir(Compass::N); + Dir d2 = Dir(Compass::CENTER); + d1 = d2; + assert(d1.asInt() == (int)Compass::CENTER); + d1 = Dir(Compass::SW); + assert(d1.asInt() == 0); + d1 = Dir(Compass::S); + assert(d1.asInt() == 1); + d1 = Dir(Compass::SE); + assert(d1.asInt() == 2); + d1 = Dir(Compass::W); + assert(d1.asInt() == 3); + d1 = Dir(Compass::CENTER); + assert(d1.asInt() == 4); + d1 = Dir(Compass::E); + assert(d1.asInt() == 5); + d1 = Dir(Compass::NW); + assert(d1.asInt() == 6); + d1 = Dir(Compass::N); + assert(d1.asInt() == 7); + d1 = Dir(Compass::NE); + assert(d1.asInt() == 8); + + assert(Dir(Compass::SW).asInt() == 0); + assert(Dir(Compass::S).asInt() == 1); + assert(Dir(Compass::SE).asInt() == 2); + assert(Dir(Compass::W).asInt() == 3); + assert(Dir(Compass::CENTER).asInt() == 4); + assert(Dir(Compass::E).asInt() == 5); + assert(Dir(Compass::NW).asInt() == 6); + assert(Dir(Compass::N).asInt() == 7); + assert(Dir(Compass::NE).asInt() == 8); + assert(Dir((Compass)8).asInt() == 8); + assert(Dir((Compass)((Dir((Compass)8)).asInt())).asInt() == 8); + assert(Dir((Compass)(Dir(Compass::NE).asInt())).asInt() == 8); + d2 = Compass::E; + d1 = d2; + assert(d1.asInt() == 5); + d2 = d1; + assert(d1.asInt() == 5); + + // .operator=() from Compass + d1 = Compass::SW; + assert(d1.asInt() == 0); + d1 = Compass::SE; + assert(d1.asInt() == 2); + + // [in]equality with Compass + d1 = Compass::CENTER; + assert(d1 == Compass::CENTER); + d1 = Compass::SE; + assert(d1 == Compass::SE); + assert(Dir(Compass::W) == Compass::W); + assert(Dir(Compass::W) != Compass::NW); + + // [in]equality with Dir + d1 = Compass::N; + d2 = Compass::N; + assert(d1 == d2); + assert(d2 == d1); + d1 = Compass::NE; + assert(d1 != d2); + assert(d2 != d1); + + // .rotate() + assert(d1.rotate(1) == Compass::E); + assert(d1.rotate(2) == Compass::SE); + assert(d1.rotate(-1) == Compass::N); + assert(d1.rotate(-2) == Compass::NW); + assert(Dir(Compass::N).rotate(1) == d1); + assert(Dir(Compass::SW).rotate(-2) == Compass::SE); + + // .asNormalizedCoord() + Coord c1 = Dir(Compass::CENTER).asNormalizedCoord(); + assert(c1.x == 0 && c1.y == 0); + d1 = Compass::SW; + c1 = d1.asNormalizedCoord(); + assert(c1.x == -1 && c1.y == -1); + c1 = Dir(Compass::S).asNormalizedCoord(); + assert(c1.x == 0 && c1.y == -1); + c1 = Dir(Compass::SE).asNormalizedCoord(); + assert(c1.x == 1 && c1.y == -1); + c1 = Dir(Compass::W).asNormalizedCoord(); + assert(c1.x == -1 && c1.y == 0); + c1 = Dir(Compass::E).asNormalizedCoord(); + assert(c1.x == 1 && c1.y == 0); + c1 = Dir(Compass::NW).asNormalizedCoord(); + assert(c1.x == -1 && c1.y == 1); + c1 = Dir(Compass::N).asNormalizedCoord(); + assert(c1.x == 0 && c1.y == 1); + c1 = Dir(Compass::NE).asNormalizedCoord(); + assert(c1.x == 1 && c1.y == 1); + + // .asNormalizedPolar() + d1 = Compass::SW; + Polar p1 = d1.asNormalizedPolar(); + assert(p1.mag == 1 && p1.dir == Compass::SW); + p1 = Dir(Compass::S).asNormalizedPolar(); + assert(p1.mag == 1 && p1.dir == Compass::S); + p1 = Dir(Compass::SE).asNormalizedPolar(); + assert(p1.mag == 1 && p1.dir == Compass::SE); + p1 = Dir(Compass::W).asNormalizedPolar(); + assert(p1.mag == 1 && p1.dir == Compass::W); + p1 = Dir(Compass::E).asNormalizedPolar(); + assert(p1.mag == 1 && p1.dir == Compass::E); + p1 = Dir(Compass::NW).asNormalizedPolar(); + assert(p1.mag == 1 && p1.dir == Compass::NW); + p1 = Dir(Compass::N).asNormalizedPolar(); + assert(p1.mag == 1 && p1.dir == Compass::N); + p1 = Dir(Compass::NE).asNormalizedPolar(); + assert(p1.mag == 1 && p1.dir == Compass::NE); + + // Coord + // ctor from int16_t,int16_t + c1 = Coord(); + assert(c1.x == 0 && c1.y == 0); + c1 = Coord(1, 1); + assert(c1.x == 1 && c1.y == 1); + c1 = Coord(-6, 12); + assert(c1.x == -6 && c1.y == 12); + + // copy assignment + Coord c2 = Coord(9, 101); + assert(c2.x == 9 && c2.y == 101); + c1 = c2; + assert(c1.x == 9 && c2.y == 101); + + // .isNormalized() + assert(!c1.isNormalized()); + assert(Coord(0, 0).isNormalized()); + assert(Coord(0, 1).isNormalized()); + assert(Coord(1, 1).isNormalized()); + assert(Coord(-1, 0).isNormalized()); + assert(Coord(-1, -1).isNormalized()); + assert(!Coord(0, 2).isNormalized()); + assert(!Coord(1, 2).isNormalized()); + assert(!Coord(-1, 2).isNormalized()); + assert(!Coord(-2, 0).isNormalized()); + + // .normalize() + // .asDir() + c1 = Coord(0, 0); + c2 = c1.normalize(); + assert(c2.x == 0 && c2.y == 0); + assert(c2.asDir() == Compass::CENTER); + c1 = Coord(0, 1).normalize(); + assert(c1.x == 0 && c1.y == 1); + assert(c1.asDir() == Compass::N); + c1 = Coord(-1, 1).normalize(); + assert(c1.x == -1 && c1.y == 1); + assert(c1.asDir() == Compass::NW); + c1 = Coord(100, 5).normalize(); + assert(c1.x == 1 && c1.y == 0); + assert(c1.asDir() == Compass::E); + c1 = Coord(100, 105).normalize(); + assert(c1.x == 1 && c1.y == 1); + assert(c1.asDir() == Compass::NE); + c1 = Coord(-5, 101).normalize(); + assert(c1.x == 0 && c1.y == 1); + assert(c1.asDir() == Compass::N); + c1 = Coord(-500, 10).normalize(); + assert(c1.x == -1 && c1.y == 0); + assert(c1.asDir() == Compass::W); + c1 = Coord(-500, -490).normalize(); + assert(c1.x == -1 && c1.y == -1); + assert(c1.asDir() == Compass::SW); + c1 = Coord(-1, -490).normalize(); + assert(c1.x == 0 && c1.y == -1); + assert(c1.asDir() == Compass::S); + c1 = Coord(1101, -1090).normalize(); + assert(c1.x == 1 && c1.y == -1); + assert(c1.asDir() == Compass::SE); + c1 = Coord(1101, -3).normalize(); + assert(c1.x == 1 && c1.y == 0); + assert(c1.asDir() == Compass::E); + + // .length() + assert(Coord(0, 0).length() == 0); + assert(Coord(0, 1).length() == 1); + assert(Coord(-1, 0).length() == 1); + assert(Coord(-1, -1).length() == 1); // round down + assert(Coord(22, 0).length() == 22); + assert(Coord(22, 22).length() == 31); // round down + assert(Coord(10, -10).length() == 14); // round down + assert(Coord(-310, 0).length() == 310); + + // .asPolar() + p1 = Coord(0, 0).asPolar(); + assert(p1.mag == 0 && p1.dir == Compass::CENTER); + p1 = Coord(0, 1).asPolar(); + assert(p1.mag == 1 && p1.dir == Compass::N); + p1 = Coord(-10, -10).asPolar(); + assert(p1.mag == 14 && p1.dir == Compass::SW); // round down mag + p1 = Coord(100, 1).asPolar(); + assert(p1.mag == 100 && p1.dir == Compass::E); // round down mag + + // operator+(Coord), operator-(Coord) + c1 = Coord(0, 0) + Coord(6, 8); + assert(c1.x == 6 && c1.y == 8); + c1 = Coord(-70, 20) + Coord(10, -10); + assert(c1.x == -60 && c1.y == 10); + c1 = Coord(-70, 20) - Coord(10, -10); + assert(c1.x == -80 && c1.y == 30); + + // operator*(int) + c1 = Coord(0, 0) * 1; + assert(c1.x == 0 && c1.y == 0); + c1 = Coord(1, 1) * -5; + assert(c1.x == -5 && c1.y == -5); + c1 = Coord(11, 5) * -5; + assert(c1.x == -55 && c1.y == -25); + + // operator+(Dir), operator-(Dir) + c1 = Coord(0, 0); + c2 = c1 + Dir(Compass::CENTER); + assert(c2.x == 0 && c2.y == 0); + c2 = c1 + Dir(Compass::E); + assert(c2.x == 1 && c2.y == 0); + c2 = c1 + Dir(Compass::W); + assert(c2.x == -1 && c2.y == 0); + c2 = c1 + Dir(Compass::SW); + assert(c2.x == -1 && c2.y == -1); + + c2 = c1 - Dir(Compass::CENTER); + assert(c2.x == 0 && c2.y == 0); + c2 = c1 - Dir(Compass::E); + assert(c2.x == -1 && c2.y == 0); + c2 = c1 - Dir(Compass::W); + assert(c2.x == 1 && c2.y == 0); + c2 = c1 - Dir(Compass::SW); + assert(c2.x == 1 && c2.y == 1); + + // raySameness() + c1 = Coord { 0, 0 }; + c2 = Coord { 10, 11 }; + d1 = Compass::CENTER; + assert(c1.raySameness(c2) == 1.0); // special case - zero vector + assert(c2.raySameness(c1) == 1.0); // special case - zero vector + assert(c2.raySameness(d1) == 1.0); // special case - zero vector + c1 = c2; + assert(c1.raySameness(c2) == 1.0); + assert(areClosef(Coord(-10,-10).raySameness(Coord(10,10)), -1.0)); + c1 = Coord{0,11}; + c2 = Coord{20,0}; + assert(areClosef(c1.raySameness(c2), 0.0)); + assert(areClosef(c2.raySameness(c1), 0.0)); + c1 = Coord{0,444}; + c2 = Coord{113,113}; + assert(areClosef(c1.raySameness(c2), 0.707106781)); + c2 = Coord{113,-113}; + assert(areClosef(c1.raySameness(c2), -0.707106781)); + + // Polar + // ctor from mag, dir + p1 = Polar(); + assert(p1.mag == 0 && p1.dir == Compass::CENTER); + p1 = Polar(0, Compass::S); + assert(p1.mag == 0 && p1.dir == Compass::S); + p1 = Polar(10, Compass::SE); + assert(p1.mag == 10 && p1.dir == Compass::SE); + p1 = Polar(-10, Compass::NW); + assert(p1.mag == -10 && p1.dir == Compass::NW); + + // .asCoord() + c1 = Polar(0, Compass::CENTER).asCoord(); + assert(c1.x == 0 && c1.y == 0); + c1 = Polar(10, Compass::CENTER).asCoord(); + assert(c1.x == 0 && c1.y == 0); + c1 = Polar(20, Compass::N).asCoord(); + assert(c1.x == 0 && c1.y == 20); + //c1 = Polar(12, Compass::W).asCoord(); + p1 = Polar(12, Compass::W); + c1 = p1.asCoord(); + assert(c1.x == -12 && c1.y == 0); + c1 = Polar(14, Compass::NE).asCoord(); + assert(c1.x == 10 && c1.y == 10); + c1 = Polar(-14, Compass::NE).asCoord(); + assert(c1.x == -10 && c1.y == -10); + c1 = Polar(14, Compass::E).asCoord(); + assert(c1.x == 14 && c1.y == 0); + c1 = Polar(-14, Compass::E).asCoord(); + assert(c1.x == -14 && c1.y == 0); + + return true; +} + +} // end namespace BS diff --git a/src/unitTestConnectNeuralNetWiringFromGenome.cpp b/src/unitTestConnectNeuralNetWiringFromGenome.cpp new file mode 100644 index 0000000..dc5b2a3 --- /dev/null +++ b/src/unitTestConnectNeuralNetWiringFromGenome.cpp @@ -0,0 +1,42 @@ +// unitTestConnectNeuralNetWiringFromGenome.cpp + +#include +#include "simulator.h" + +namespace BS { + +void unitTestConnectNeuralNetWiringFromGenome() +{ + Indiv indiv; + Genome genome1 { +// { SENSOR, 0, NEURON, 0, 0.0 }, +// { SENSOR, 1, NEURON, 2, 2.2 }, +// { SENSOR, 13, NEURON, 9, 3.3 }, +// { NEURON, 4, NEURON, 5, 4.4 }, +// { NEURON, 4, NEURON, 4, 5.5 }, +// { NEURON, 5, NEURON, 9, 6.6 }, +// { NEURON, 0, NEURON, 0, 7.7 }, +// { NEURON, 5, NEURON, 9, 8.8 }, +// { SENSOR, 0, ACTION, 1, 9.9 }, +// { SENSOR, 2, ACTION, 12, 10.1 }, +// { NEURON, 0, ACTION, 1, 11.0 }, +// { NEURON, 4, ACTION, 2, 12.0 } + }; + + indiv.genome = { genome1 }; + + indiv.createWiringFromGenome(); + + for (auto & conn : indiv.nnet.connections) { + std::cout << (conn.sourceType == SENSOR ? "SENSOR" : "NEURON") << " " + << conn.sourceNum << " -> " + << (conn.sinkType == ACTION ? "ACTION" : "NEURON") << " " + << conn.sinkNum << " at " << conn.weight << std::endl; + } + +// for (unsigned n = 0; n < indiv.nnet.neurons.size(); ++n) { +// std::cout << n << " " << indiv.nnet.neurons[n].numInputs << " inputs" << std::endl; +// } +} + +} // end namespace BS diff --git a/src/unitTestGridVisitNeighborhood.cpp b/src/unitTestGridVisitNeighborhood.cpp new file mode 100644 index 0000000..efc2a21 --- /dev/null +++ b/src/unitTestGridVisitNeighborhood.cpp @@ -0,0 +1,35 @@ +// unitTestGridVisitNeighborhood + +#include +#include "simulator.h" + +namespace BS { + +void unitTestGridVisitNeighborhood() +{ + // prints each coord: + auto printLoc = [&](Coord loc){ std::cout << loc.x << ", " << loc.y << std::endl; }; + + std::cout << "Test loc 10,10 radius 1" << std::endl; + visitNeighborhood(Coord { 10, 10 }, 1.0, printLoc); + + std::cout << "\nTest loc 0,0 radius 1" << std::endl; + visitNeighborhood(Coord { 0, 0 }, 1.0, printLoc); + + std::cout << "\nTest loc 10,10 radius 1.4" << std::endl; + visitNeighborhood(Coord { 10, 10 }, 1.4, printLoc); + + std::cout << "\nTest loc 10,10 radius 1.5" << std::endl; + visitNeighborhood(Coord { 10, 10 }, 1.5, printLoc); + + std::cout << "\nTest loc 1,1 radius 1.4" << std::endl; + visitNeighborhood(Coord { 1, 1 }, 1.4, printLoc); + + std::cout << "\nTest loc 10,10 radius 2.0" << std::endl; + visitNeighborhood(Coord { 10, 10 }, 2.0, printLoc); + + std::cout << "\nTest loc p.sizeX-1, p.sizeY-1 radius 2.0" << std::endl; + visitNeighborhood(Coord { (int16_t)(p.sizeX-1), (int16_t)(p.sizeY-1) }, 2.0, printLoc); +} + +} // end namespace BS diff --git a/tools/graph-nnet.py b/tools/graph-nnet.py new file mode 100644 index 0000000..50df1fe --- /dev/null +++ b/tools/graph-nnet.py @@ -0,0 +1,73 @@ +#!/usr/bin/python3 + +import igraph + +# load data into a graph +g = igraph.Graph.Read_Ncol('net.txt', names=True, weights=True) + +for v in g.vs: + v['size'] = 35 + v['label'] = v['name'] + if v['name'] in ['Lx', 'Ly', 'EDx', 'EDy', 'ED', 'Bfd', 'Blr', 'Gen', 'LMx', 'LMy', 'LPf', 'LPb', 'Pop', 'Pfd', 'Plr', 'Osc', 'Age', 'Rnd', 'Sg', 'Sfd', 'Slr']: + v['color'] = 'lightblue' + elif v['name'] in ['MvX', 'MvY', 'MvE', 'MvW', 'MvN', 'MvS', 'Mfd', 'MvL', 'MvR', 'MRL', 'Mrv', 'Mrn', 'OSC', 'LPD', 'Res', 'SG', 'Klf' ]: + v['color'] = 'lightpink' + else: + v['color'] = 'lightgrey' + +# convert edge weights to color and size +for e in g.es: + #print(e['weight']) + if e['weight'] < 0: + e['color'] = 'lightcoral' + elif e['weight'] == 0: + e['color'] = 'grey' + else: + e['color'] = 'green' + + width = abs(e['weight']) + e['width'] = 1 + 1.25 * (width / 8192.0) + + +# plot graph + +print(len(g.vs)) + +if len(g.vs) < 6: + bbox = (300,300) + layout = 'fruchterman_reingold' +elif len(g.vs) < 12: + bbox = (400,400) + layout = 'fruchterman_reingold' +elif len(g.vs) < 18: + bbox = (500,500) + layout = 'fruchterman_reingold' +elif len(g.vs) < 24: + bbox = (520,520) + layout = 'fruchterman_reingold' +elif len(g.vs) < 26: + bbox = (800,800) + layout = 'fruchterman_reingold' +elif len(g.vs) < 50: + bbox = (1000,1000) + layout = 'fruchterman_reingold' +elif len(g.vs) < 130: + bbox = (1200,1000) + layout = 'fruchterman_reingold' +elif len(g.vs) < 150: + bbox = (4000,4000) + layout = 'fruchterman_reingold' + for v in g.vs: + v['size'] = v['size'] * 1.5 +elif len(g.vs) < 200: + bbox = (4000,4000) + layout = 'kamada_kawai' + for v in g.vs: + v['size'] = v['size'] * 2 +else: + bbox = (8000,8000) + layout = 'fruchterman_reingold' + +igraph.plot(g, "net.svg", edge_curved=True, bbox=bbox, margin=64, layout=layout) + + diff --git a/tools/graphlog-genomeLength.gp b/tools/graphlog-genomeLength.gp new file mode 100644 index 0000000..2c89bb9 --- /dev/null +++ b/tools/graphlog-genomeLength.gp @@ -0,0 +1,28 @@ +#!/usr/bin/gnuplot --persist + +# Requires a text file named "epoch-log.txt" in the log directory + +set term png size 2000, 400 +set output "/home/dm/sw/biosim4/images/log.png" + +# Left Y axis gets scaled to the max survivors. +# Right Y axis gets scaled to 0..255. +# 1:2 Survivors 0..N => 0..N +# 1:3 Genome length 0..50 => 0..255 +# 1:4 Diversity 0..0.7 => 0..255 +# 1:5 Anxiety 0..255 => 0..255 + +set mxtics +set ytics autofreq nomirror tc lt 2 +set y2range [ 0:256 ] +set y2tics autofreq nomirror tc lt 1 +set grid +set key lmargin + +ScaleGenomeLength(y)= y*2 +ScaleDiversity(d)= 350*d + +plot "/home/dm/sw/biosim4/logs/epoch-log.txt" using 1:2 with lines lw 1 linecolor 2 title "Survivors", \ + "" using 1:(ScaleDiversity($3)) with lines lw 1 linecolor 1 title "Diversity" axes x1y2, \ + "" using 1:(ScaleGenomeLength($4)) with lines lw 1 linecolor 6 title "Genome Len" axes x1y2 + diff --git a/tools/graphlog.gp b/tools/graphlog.gp new file mode 100644 index 0000000..91304ce --- /dev/null +++ b/tools/graphlog.gp @@ -0,0 +1,31 @@ +#!/usr/bin/gnuplot --persist + +# Requires a text file named "epoch-log.txt" in the log directory + +set term png size 2000, 400 +set output "../images/log.png" + +# Left Y axis gets scaled to the max survivors. +# Right Y axis gets scaled to 0..255. +# 1:2 Survivors 0..N => 0..N +# 1:3 Genome length 0..50 => 0..255 +# 1:4 Diversity 0..0.7 => 0..255 +# 1:5 Anxiety 0..255 => 0..255 + +set mxtics +set ytics autofreq nomirror tc lt 2 +set yrange [ 0:8000 ] +set y2range [ 0:1 ] +set y2tics autofreq nomirror tc lt 1 +set grid +set key lmargin + +ScaleSurvivors(s) = s +ScaleGenomeLength(y)= y*2 +ScaleDiversity(d)= d +#ScaleMurders(m) = m + +plot "../logs/epoch-log.txt" \ + using 1:(ScaleSurvivors($2)) with lines lw 2 linecolor 2 title "Survivors", \ + "" using 1:(ScaleDiversity($3)) with lines lw 2 linecolor 1 title "Diversity" axes x1y2 +