Colonization models: a programming tutorial (Part 1/2)

16th May 2011

6MrMind

2snarles

6JenniferRM

1Anubhav

1Anubhav

1PhilGoetz

New Comment

6 comments, sorted by Click to highlight new comments since: Today at 11:38 AM

Two observations to improve the model.

The first is that in this code, you generate only one civilization per model, and this penalizes new civilizations if the probability for the old ones to be filtered out of existence is low. A more realistic model should create a number N of civilization home planets at each step.

The second and IMO more important one is that you assume a three dimensional sphere as a model of distance, so that for example two points at the exact opposite of the sphere are maximally distant. But if we are to accept the Big Bang model, we live in the surface of an expanding 4d sphere, and this means that two points at the polar opposite are **the same** point. This drastically increases the connectivity of the space.

@MrMind. Indeed, in a later edit I pointed out that a more realistic simulation would allow for the creation of multiple new civilizations at a time step.

However, your suggestion of modelling the geometry of the universe as a 4d sphere is not only relevant, but also easy to implement. I am not so sure how to handle the expansion. I admit to being rusty on my cosmology.

If it was not 10 minutes to my bedtime I'd work through your tutorial, but it is so I won't. Still I wanted there to be a comment and it seems worthwhile to point out that the conjunction of the fermi paradox and R in the same post is *awesome*.

Trying to implement this in Python. With the given probabilities, around 170-180 civs spawn in 2000 epochs, but none of them ever transition to type II.

Is this the expected result or is there a bug in my code?

**Edit**: Was a bug.

I notice I get almost no type 0s/2bs at the end of this simulation.

The 0s quickly die out or (more rarely) ascend to 2a/2b, and the 2bs quickly ascend to 3.

(And this is *before* implementing the '3s kill everything in their path' thing.)

Thank you for being concrete!

In the model we propose in the post, the above parameters are held to be constant throughout the entire history of the universe.

You might also try a model in which these aren't constant, because the presence of neighboring civilizations can change the probability of going from one type to the next.

IntroductionAre we alone in the universe? How likely is our species to survive the transition from a Type 0 to a Type II civilization? The answers to these questions would be of immense interest to our race; however, we have few tools to reason about these questions. This does not stop us from wanting to find answers to these questions, often by employing controversial principles of inference such as 'anthropic reasoning.' The reader can find a wealth of stimulating discussion about anthropic reasoning at Katja Grace's blog, the site from which this post takes its inspiration. The purpose of this post is to give a quantitatively oriented approach to anthropic reasoning, demonstrating how computer simulations and Bayesian inference can be used as tools for exploration.

The central mystery we want to examine is the

Fermi paradox: the fact thatOne explanation for the Fermi paradox is that we are the only intelligent civilization in the universe. A far more chilling explanation is that intelligent civilizations emerge quite frequently, but that all other intelligent civilizations that have come before us ended up destroying themselves before they could manage to make their mark on their universe.

We can reason about which of the above two explanations are more likely if we have the audacity to assume

a modelfor the emergence and development of civilizations in universe 'similar to ours.' In such a model, it is usually useful to distinguish different 'types' of civilizations. Type 0 civilizations are civilizations with similar levels of technology as ourselves. If a Type 0 civilization survives long enough and accumulates enough scientific knowledge, it can make a transition to a Type I civilization--a civilization which has attained mastery of their home planet. A Type I civilization, over time, can transition to a Type II civilization if it colonizes its solar system. We would suppose that a nearby civilization would have to have reached Type II in order for their activities to be prominent enough for us to be able to detect them. In the original terminology, a Type III civilization is one which has mastery of its galaxy, but in this post we take it to mean something else.The simplest model for the emergence and development of civilizations would have to specify the following:

In the model we propose in the post, the above parameters are held to be constant throughout the entire history of the universe. The importance of the model is that after given a particular specification of the parameters, we can apply Bayesian inference to see how well the model explains the Fermi paradox. The idea is to simulate many different histories of universes for a given set of parameters, so as to find the

expected number of observers who observe the Fermi paradoxgiven a particular specification of the parameters. More details about Bayesian inference given in Part 2 of this tutorial.This post is targeted at readers who are interested in simulating the emergence and expansion of intelligent civilizations in 'universes similar to ours' but who lack the programming knowledge to code these simulations. In this post we will guide the reader through the design and production of a relatively simple universe model and the methodology for doing 'anthropic' Bayesian inference using the model.

I. The modelWe are going to simulate the entire history, or time-line, of a universe. We assume that the universe has a finite lifespan; readers are invited to attempt to extend the model to allow for indefinite lifespans. For computational simplicity, we will break up the history of the universe into, say, 1000 discrete epochs, where each epoch might represent twenty million years. The epoch number is represented by a variable

tranging from 1 to 1000. The epocht = 1does not represent the very beginning of the universe, but rather, the earliest period where intelligent life has a significant chance of appearing; likewise, the last epocht=1000 would be the last period in which the cosmology of the universe allows for the emergence of intelligent life.We assume the cosmological model of space as a (the three-dimensional surface) of a four-dimensional unit sphere centered at the origin. That is, the universe is represented as the set of all points

(x,y,z,w)wherex^2 + y^2 + z^2 + w^2 = 1. Note that since w can be calculated fromx, y, z, only the pointsx, y, zneed be stored in the simulation.Now it is important to correctly specify thespeed of light as

speed of lightin terms of the units we have chosen for time and space. Assuming the maximum geodesic of the universe is 100 billion light years, and given that each time step is twenty million years, this would put thes = 4 x 10^-4in our units.One could go further to populate this universe with randomly located habitable planets. However, we would need a lot of planets for a realistic model, so this approach would not be particularly computationally efficient. Instead, we will

assume that the number of habitable planets in a region of space is proportional to its volume.This will allow us to do the simulation without specifying the location of the planets.Now we need to specify the mechanism by which new Type 0 civilizations emerge in our universe. For the sake of simplicity,

at most one new Type 0 civilization appears in each time step.(A more realistic model might employ a Poisson distribution to determine the number of new civilizations at each time step.) The location of the new civilization is located uniformly at random within the uninhabited volume of the universe. We will suppose that no new Type 0 civilization can emerge on planets already inhabited by others civilizations. Thus, the probability that a new Type 0 civilization emerges must depend on how much of the universe is uninhabited. We will posit that this probability is a constant factor,a, times the proportion of universe which is uninhabited. A more advanced model might have the propensity for life factoravary with time, to account for planetary and stellar evolution.Next we need to specify how existing civilizations transition to the next stage, self-destruct, become visible or colonize other planets. We will use four types of civilizations:

all Type III civilizations are expansionist.That is, once a civilization becomes Type III, it starts colonizing the surrounding space as quickly as possible. Space colonized by Type III civilizations issterilized, it can not longer produce new Type 0 civilizations.We have not included the "Type I" civilization in our model, since for our purposes Type I civilizations are virtually identical to Type 0 civilizations. In the same way, Type IIa civilizations would be virtually identical to 'peaceful' Type III civilizations for the purposes of out model, which is why hold all Type III civilizations to be expansionist.

Here will be our model for this process:

k * s.Here we have the speed multiplierkas constant, but a more realistic model might have the speed of colonization increase as the Type III civilization gains more technologyeof transitioning to a Type III civilization.bof self-destructing and probabilitycof transitioning to Type IIa and probabilitydof transitioning to Type IIb (and remains Type 0 with probability 1-b-c-d).a, locate a random pointxwithin the unit sphere which is a candidate site for a new Type 0 civilization. If that point is within the colonization sphere of any existing civilization, no new Type 0 civilization is generated. Otherwise, create a new Type 0 civilization at that point.To sum up, each civilization in our universe has four variables associated with it.

In each step of the simulation, we need to store these quantities in a data structure. The data structure for the simulation should have a slot for each time step of the simulation; in each slot, we will store a list of civilizations with their associated quantities.

After we have run the simulation, we can extract various quantities of interest from the data structure, for example, we might want to know:

II. ImplementationThe program we will be using is the R programming environment, an open-source computing environment commonly used for statistical data analysis. No programming experience is assumed, but readers are expected to be able to learn by example. Conversely, readers with programming skills can skip to the quoted code blocks, reading the explanations when needed.

After you install and run R, you should see a blank terminal. To run the code, simply copy and paste the provided code into the buffer and press 'enter.' We give the code in blocks so as to make copying and pasting more convenient. Here is the first block of code--don't be alarmed! We will walk through the code line-by-line immediately afterwards.

Now we'll explain the code for "Example 1", starting from the top.

ato 1/100. Recall this is the base rate for a Type 0 civilization to emerge in a time step.b, c, d. The semicolon ';' separates multiple commands that are on the same linee, (type IIb -> type III transition)kto 0.9. That is, the maximum speed of colonization is taken to be 90% of the speed of light.civsis a data structure for storing the information of all the civilizations which exist at a particular time step. It will be a list, with an element for each civilization in existence at that timeearth_coordinatesis a vector with three elements: (0.1, 0.2, 0.3). Thec()function forms a vector out of the enclosed numbers. As you might guess, this particular vector is the spatial coordinates of planet Earth in our example.human_civis a data structure containing all the information needed to describe (our) civilization:locationwill be the vector given by earth_coordinates: (0.1,0.2,0.3);type=0indicates a Type 0 civilization,r_vis the radius of our sphere of visibility, which is zero because we haven't reached Type II yet.r_cis the radius of our sphere of colonization, which is zero because we haven't reached Type III yet.zerg_civholds the information for another civilization, which is older than ours by 400 million years. It has a home planet with coordinates (0.2,0.1,0.2), which has reached Type III, and which can be seen by all observers within distance 0.2, and which has colonized all planets within distance 0.1 of its home planet!civs[[1]] = human_civsets the first element of the list civs to be the data structure encapsulation information about our civilization. Yes, you can have the elements of lists be other lists--we will be doing that a lot.civs[[2]] = zerg_civsets the second element of the list civs to be the information about the Zerg civilization(In this code, we will have the parameters a, b, c, d, e, k,s be global variables. This is not the best coding practice, but it simplifies the code for this tutorial.)

It is certainly cause for alarm that there is a hostile Type III civilization in our vicinity. Can we see them? Have we already been colonized by them? To answer these questions, we need to program a distance function. We'll be lazy here and use the Euclidean distance, instead of calculating the length of the geodesic on a 3-sphere.

distanceis the name of the function.= function(x,y)tells R that this is to be a function with two arguments, x and y. The{brace, or 'begin-block symbol', means that the subsequent code will be part of the functionthereturn()command is for returning the 'answer' computed by the function. In this case, the answer is the Euclidean distance of the vectors x and y;sqrt()evaluates the square root of the argument.sum()evaluates the sum of the components of a vector.^2applied to a vector squares each component of the vector.}or 'end-block symbol' indicates that we are done writing the code for the functionAfter you copy and paste the above code, there will be a new function in the environment, distance(). We will use it now

Line 2: look back at how we defined the list human_civ. The

$symbol allows you to access the element of that list labeled as the 'location.'Line 4: Recall x is the coordinates of our planet, and y is the coordinates of the zerg home planet. We call the function we just defined to find the distance between our home planet and the zerg home planet, and then we store the answer in a new variable, d

Line 5: A boolean expression: is d less than the zerg's radius of visibility?

Line 6 & 7: The if() statement executes the code in brackets only if the statement in the parentheses is true

The '<' symbol is obviously the comparison 'less than'. Less obvious: '>=' for greater-or-equal-to, '==' for equal to. A common newbie mistake is to confuse the assignment '=' symbol for the equality '==' symbol--they have very different meanings!

EXERCISE 1: Add some additional civilizations to the list

civsmanually. We code the civilization types as followsNow we start on coding the simulation step. Recall the simulation step consist of two parts. In the first step, we roll to see if any of the existing civilizations transition to another stage, and we update the radii of visibility and colonization for Type II+ civilizations. In the second step, we roll to see if a new Type 0 civilization emerges.

We will code a separate transition function for each Type of civilization. We'll start with the Type IIb and Type III civilizations, they are the easiest.

Notice that there is no return() statement in either function. That is because a function automatically returns the value in the last line (

civin this case) whenever it makes it all the way to the last line.Now, test to see if the transition function works properly on the Zerg civ. Also test the transition for Type 2b planets on any Type 2b planet you created in Exercise 1.

Notice how the transition function works. The transition function returns an updated version of the civ it was given. But it does not update the original civ automatically: you have to make sure that the updated state gets stored in the appropriate place.

Now, we want to make a function which automatically chooses the right transition for a given civ. Also, we can make it keep track of the civilization's age.

EXERCISE: What happens if you call the transition() function on a dead (type=-1) civilization?

In order to handle the transition rolls for type 0 and type IIa civilizations, we need a way to generate random numbers. R has a built-in function to generate random numbers from 0 to 1, the

runif()function, which we demonstrate.The command

runif(1)returns one random number between 0 and 1. See that the boolean statement(runif(1) < e)is true with probability e.Now the transition for a Type 0 civilization is more complicated: it has a probability of dying, of going to Type IIa, of going to Type IIb, or staying type 0, and all these events are mutually exclusive. There is a simple trick for handling all that with a single roll.

EXERCISE: Does the above code indeed transition the civ to: dead with probability b, Type IIa with probability c, Type IIb with probability d?

Finally, we make a function for automatically applying the transition function to every civ in a list. How we will go through the list is to use a

for()loop to access each of the elements civs[[1]], civs[[2]], all the way up to civs[[no_civs]].civs = transition_all(civs)to update the listcivslength()returns the number of elements contained in a list.iinitialized at 1. At the end of every iteration, it incrementsiby one. It stops iterating wheniexceedsno_civs.ith elementcivsEXERCISE: Check to see that this properly updates the civilizations in your list. And by the way, what civilization does Humanity become in the end?

Now the last step of the simulation is to randomly decide if a new civilization emerges.

PLANET-GENERATION PROCEDURE

a, no new civilization emerges, and we exit the procedure.xin the universe. If that point is in an uninhabited region of space we create a new Type 0 civilization with home planet located atx. But if pointxis located within the sphere of colonization of any Type III civilization, no new civilization is generated.For this, we need a way to draw random points on the unit 3-sphere, and a method for checking whether a point lies within the sphere of colonization of an existing civilization. There is a cool trick from probability theory for generating points from the surface of an N-sphere, which takes advantage of the fact that the multivariate Gaussian has spherical symmetry. The

rnorm()function samples from the standard Gaussian with mean 0 and variance 1. After we generate a vector which has 4 independently and identically distributed mean 0 Gaussian components, all we have to do is normalize it so that it has radius 1.EXERCISE: Try out the new function by entering the command

random_in_sphere(), see that you get random points with radius 1.Next we need to code the function

is_uncolonized(x), which returns FALSE if the point x is within the sphere of colonization of an existing civilization, and TRUE otherwise. First, we check to see if there are any existing civilizations in the universe at all. If not, there is nothing to check, and we return TRUE instantly. But if there are any existing civilizations, we will have to see if our point x is within the colonization sphere of any of the existing civilizations. We will go through the list of existing civilizations one by one. using a for() loop. If at any point we find that x is within a colonization sphere, we return FALSE. However, if we make it through the for() loop, that means the point x was in none of the colonization spheres, so we can return TRUE.EXERCISE: Test these functions.

We are almost done with the simulation part of the code. We write a function to implement the planet generation step, and another to handle the simulation step in its entirety.

EXERCISE: Run a few simulation steps. Increase the parameters a, b, etc. if things are too slow.

Now, we'll write a function for generating the entire simulation history. The function will return a list,

history, whosetth element will be the list of civilizationscivsat time steptof the simulation.EXERCISE: Run an entire simulation. How many civilizations are around at time 1000?

Fantastic! Now we can generate entire simulated histories of universes with adjustable parameters! Now what?

The reader quickly finds that it's a cumbersome task to manually inspect all of these

historyobjects generated by the simulation.We leave it as a challenge to any reader experienced in R to come up with a function for visually displaying

historyobjects in a pleasing manner.In the second part of this tutorial, we will generate lots and lots of

historyobjects, and write functions for extracting relevant information from these objects, in order to do Bayesian inference.EXERCISE: Modify the simulation so that all civilizations of type II or lower which come within the colonization sphere of a Type III civilization are destroyed. (Remember that destroyed Type II civs are still 'visible'!)

EXERCISE: Implement the correct distance function for the length of the geodesic between two points on a 3-sphere. Does this affect results in Part 2?