Tutorial #8 ---------------------------------------------------------------------- Solving the eight queens problem in GENESIS using a network of spiking neurons. ---------------------------------------------------------------------- If you've completed tutorial 7, you've had the wonderfully illuminating experience of setting up a network model in GENESIS which gives rise to a lot of colorful blinking lights. While this may be of limited use in impressing supervisors, lab visitors and funding agencies, what we'd really like to do is to use networks within GENESIS to do some interesting computations, and hopefully expand our knowledge of what more-or-less biologically-based neural networks are capable of. In this tutorial we will adapt the network in tutorial 7 to solve a famous optimization problem, the eight queens problem. In so doing we will have to modify the basic network from tutorial 7 in a number of small and not-so-small ways to get it to do what we want it to do. Of course, we're not claiming that this has anything to do with what the brain is doing (try solving the eight queens problem yourself if you disagree!). This tutorial is laid out as follows. First we'll describe the problem. Next we'll look at how to modify the network in tutorial 7 to solve it. Then we'll look at how to set up a special genesis elements to allow us to visualize the solution. Finally we'll discuss how to set up the solution to run in batch mode so the machine can find solutions for us automatically. ---------------------------------------------------------------------- The eight queens problem ---------------------------------------------------------------------- The problem is: suppose you have a chessboard (64 squares arranged as 8 squares by 8 squares) and eight mutually hostile queens. As you recall, a queen can move any number of squares in a horizontal, vertical or diagonal line. In a desperate attempt to keep the peace, you decide that you will arrange the eight queens on the chessboard so that no queen can capture another queen on the next move. This is not as trivial as it sounds (try it with pencil and paper). It is known that this problem has multiple independent solutions (by independent I mean that the variant solutions are not just rotations and reflections of one another). In addition, this problem has been successfully attacked using traditional connectionist-type neural networks. The trick is to see if we can modify our little GENESIS network of spiking cells to get solutions to this problem. We will describe one way of solving the problem below. We don't claim that there are no other ways to arrange the network to solve this problem, and we will point out a few of the variations that the enterprising GENESIS hacker may wish to try out. ---------------------------------------------------------------------- Getting started ---------------------------------------------------------------------- First, copy all the files from the Scripts/tutorials/tutorial7 directory to one of your own directories, called, for instance, ~jrhacker/tutorial8. Rename tutorial7.g to tutorial8.g. We will make a number of modifications to these script files as well as creating some new ones. ---------------------------------------------------------------------- Changing the cells ---------------------------------------------------------------------- The connection pattern we will use to solve the eight queens problem will depend critically on inhibitory synaptic activity. In order to get the maximum inhibitory effect, change the location of the inhibitory synapse from the dendrite of the cell to the soma by modifying "net_cell.p". This is also more biological, in that inhibitory synapses onto cortical pyramidal cells tend to be clustered near the soma. Also, it will prove useful later to increase the decay time constants for the excitatory and inhibitory synaptic channels (synchans). This will make the inhibitory influence of the cells on their neighbors more long-lasting and hence more effective. Increasing the decay time constant of the excitatory channels on the dendrite will have the effect of smoothing out the random spiking pattern of the inputs, leading to a relatively constant excitatory drive upon which a little noise is superimposed. In the file "synchans.g", increase tau2 for Ex_channel from 3 msec to 10 msec. Then increase tau2 for Inh_channel from 20 msec to 50 msec. This is not very realistic, but see the exercises at the end of the tutorial. ---------------------------------------------------------------------- Setting the connection strengths ---------------------------------------------------------------------- One way of solving the eight queens problem is to set the network up so that it has *attractors* (which are patterns of firing which are stable over time, at least in some statistical sense) and then to arrange the connection strengths and polarities (excitatory or inhibitory) so that some of these attractors solve the eight queens problem (ideally, they all will, but we'll get to that). The easiest way to do this is to set up the connection pattern so that a given cell will inhibit other cells which are located on a horizontal, vertical or diagonal line relative to the cell. Edit "net_conn.g" from tutorial 7 as follows: a) Remove all excitatory connections between network cells (but not from the inputs to the network cells!). Remove all the inhibitory connections between cells in the network too. b) Use the "planarconnect" command to set up inhibitory connections between cells in the network. Make the source equal to the spike generator (spikegen) objects on all the cells. Make the destination equal to the inhibitory channels on the somas of the cells (soma/Inh_channel). Make the source mask cover all the cells. Make sure that you add a "-desthole" option to prevent cells from inhibiting themselves. The easiest way to do this is to use the -relative option and specify the limits as a very small box, i.e. planarconnect etc. ... -desthole box -0.01 -0.01 0.01 0.01 Also set the -destmask option to cover a horizontal line of cells relative to each cell. c) Copy the above command. Change the -destmask arguments to cover a vertical line of cells relative to each cell. d) Now you have to set up the inhibitory connections along the diagonals. Make no mistake: this is a pain. The strategy we used is to use a series of planarconnect commands, each one of which connects each cell to neighboring cells on the four diagonal corners. The planarconnect commands are placed in a "for" loop, and each time through the loop you set up connections with cells further from the starting cell. You need four -destmask options, one for each corner. The arguments to the -destmask option should be variables which are set each time GENESIS goes through the loop. Now that all the connections are set up, the weights and delays need to be specified. We use the "planarweight" and "planardelay" commands, described in chapter 17 of the Book of Genesis and in the GENESIS reference manual, to set the weights and delays for all cells in the network from single commands. Since we don't care about delays, we can set them all to zero with the command planardelay /network/cell[]/soma/spike -fixed 0.0 Similarly, the weights can all be a single fixed value, so we can use the -fixed option of planarweights as well. Try a fixed weight of 15.0 for the intrinsic network connections (spikegen to inhibitory synapse) and change the input weights to 5.0. ---------------------------------------------------------------------- Visualizing the solution ---------------------------------------------------------------------- When the simulation is run now, after a time you may see some of the cells of the network firing frequently while others are silent. However, it would be nice if you could display the average firing frequency of the cells, since that is what will encode the solution to the problem. In order to do this, we will create an array of elements called "RC" elements. These are merely leaky integrators which take their inputs (in this case the membrane potentials (Vm) of the cells) and pass them through an RC circuit, giving rise to outputs which reflect the average membrane potentials of the cells (which in turn reflects the average firing rates). The equivalent circuit of the RC object is: Iinject --> ________________________________ <-- state | | | | \ | / ___|___ R \ _______ C / | | | V0 --- | ------- | | | | | ------------------- | ------- ----- --- - The RC object has four useful fields (from our standpoint): V0, R, C, and "state". The "state" field is the output of the object. We can get a reasonable averaging of Vm by using a time constant of 30 msec. The time constant is equal to R*C, and since it doesn't matter what the values of either R or C alone are, we can set R to 1.0 and C to 0.030 to get an RC time constant of 30 msec. V0 is the level that the state field is set to on a reset, and is also what the state field will approach in the absence of any inputs. We set it to the reversal potential of the soma, which is -0.070 (-70 mV; this is EREST_ACT for the soma). Open a new file called "chess_avg.g". Create a prototype RC element called /library/RC with the above parameters. Then create a two-dimensional array of RC elements by using the "createmap" command: create neutral /avg createmap /library/RC /avg 8 8 -delta 0.1 0.1 -origin 0.15 0.15 This will create an array of 64 RC elements arranged physically in a two-d grid. The origin is set to (0.15, 0.15) in order to improve the appearance of the array when displayed in an xview object (which we'll get to below). We then set Iinject by adding INJECT messages from the somas of the cells in the network to the RC elements. By the way, Iinject is not a field of the RC object; it's just included in the diagram for clarity. There is an "inject" field on RC objects, but we won't use it; it's for steady current injection only. Since the input is a voltage, not a current, we would want to inject (Vm - V0)/R (which is a current) into the RC element in order to give an asymptotic state value of Vm for a constant Vm input. Since R = 1.0, this simplifies to injecting (Vm - V0). Ideally we would like to be able to write: addmsg /network/cell[{i}]/soma /avg/RC[{i}] INJECT (Vm - Erev) for each cell i, where Erev equals -0.070. However, this isn't possible in GENESIS, so we have to divide this into two commands: addmsg /network/cell[{i}]/soma /avg/RC[{i}] INJECT Vm addmsg /network/cell[{i}]/soma /avg/RC[{i}] INJECT {-Erev} This is OK since you can always use a constant value (like -Erev) in a message. Now put the above two lines in a for loop and all the messages will be set up correctly (don't forget to declare i as an integer and declare and define Erev). Note that each soma sends messages to only one RC element. Now that we have the RC elements in place, we have to have a way of looking at their outputs. Open a new file called "chess_out.g". Since we defined the "make_xview" function in tutorial 7 to create xview widgets, it's easy to make a new xview widget to display the outputs (states) of the RC elements. Just type make_xview /avg/RC[] state \ /xavg "Average Vm of network cells:" \ 700 450 \ -0.070 -0.060 This creates an xview element called /xavg/drawsquare/cells (look at the definition of make_xview in net_out.g to see why). The last two numbers show that the range of the output is between -70 mV and -60 mV. Anything more negative than -70 mV will look black and anything more positive than -60 mV will look white. If the active cells are spiking frequently the average Vm will probably be higher than -60 mV and thus they will show up as white. Once the network has converged to a solution the active cells should look white while all the rest of the cells will look black. At this point you can add the lines: include chess_avg.g include chess_out.g to tutorial8.g and run the simulation. You will see a random pattern of firing on the input units as in tutorial 7. Initially the network units will show a burst of excitation due to the random inputs, followed by a long quiescent period as the inhibition kicks in. After about 150 msec the inhibition will start to wear off and some cells will become active again. They will inhibit the cells that are along the same vertical, horizontal or diagonal line and a solution will start to develop. The solution is most easily viewed on the xview object in the lower right-hand corner, which displays the output of the RC elements. In many cases, a true eight-queens solution will not be found, since there are a lot of local minima for this network with these connections which only involve 6 or 7 units active. If this is the case, just restart the simulation and try again. About one out of every ten attempts should give a solution with eight cells active i.e. a true eight queens solution. It is rather tedious to keep restarting the simulation and waiting to see if it gives you the desired results. It would be nice if the whole process could be automated. This leads us to the next topic. ---------------------------------------------------------------------- Running the simulation in batch mode, or: >>> Real Modelers Don't Use Graphics <<< , or: "Graphics? We don't need no stinkin' graphics!" ---------------------------------------------------------------------- Now we're going to discuss two of the most fundamental laws of nature for neural modelers. They can be stated as: 1) Don't use graphics. It slows you down. 2) Automate everything. Make the computer do the work. There is a considerable computational overhead associated with using graphical objects in a simulation. Each time step, the output of one or more elements has to be sent to the graphics objects and displayed. If we could eliminate all such objects and do our simulation without graphics, the speed of the simulation would be much greater. Just for fun, go into "tutorial8.g" and change the line setclock 1 {dt * 2} // for outputs to setclock 1 {dt * 20} // for outputs The output displays will be much jerkier and rather unpleasant to look at, but the simulation will run much faster. Note that we haven't changed the simulation time step (which is clock 0), just the rate at which outputs are displayed. Since the graphics objects do no computation of their own, we might think we could just delete the graphical objects (or never load the script files containing their definitions in the first place). While this is true, we would then need a way of evaluating the results that doesn't depend on us looking at a graphical object and noting whether the network has converged to a solution to the problem. Thus, applying law # 1 will require us to devise a solution of law # 2. We will put all the commands relating to running the simulation in batch mode in a new file called "batch.g". In the file tutorial8.g, put the following line at the beginning: int batch = 0 // flag for whether to run in batch mode If batch is 0, the simulation will be run using graphics, while if batch is 1, it will be run in batch mode. The variable "batch" is thus a "flag" or boolean variable, although we are defining it as an integer since there is no boolean type in GENESIS. Now, after the line include chess_avg.g put the lines: include batch.g if (!batch) // send output to Xodus objects, not files include net_out.g include chess_out.g end Remove the other line including net_out.g, of course. Now, the files "net_out.g" and "chess_out.g", which define the XODUS interface for the simulation, will only be included if batch = 0. Finally, at the end of the file, add the lines: if (batch) run_batch quit end If batch = 1, this will cause the simulation to execute the function "run_batch" (see below) and then quit. That's it for tutorial8.g. The rest of this section will discuss the functions in batch.g. Assuming that we've run the simulation without graphics for 0.5 seconds (which is long enough to get convergence to a stable attractor in this network) how can we test the solution to see whether or not it is a solution to the eight queens problem? One way is to loop through all the RC elements after the simulation is over and count which ones have a potential over some threshold. By trial and error, we determined that a threshold of -65 mV is sufficient to distinguish between an active and an inactive cell. Thus, we can write a function called "count_active_cells" as follows: // Count the number of active cells at the end of the simulation. float thresh = -0.065 // threshold for considering a cell "active" function count_active_cells int i, count = 0 for (i = 0; i < 64; i = i + 1) if ({getfield /avg/RC[{i}] state} > {thresh}) count = count + 1 end end return {count} end Now, if this function returns 8 after running the simulation for 0.5 sec of simulated time, we most probably have a solution to the eight-queens problem. If it returns anything else, we definitely do not have such a solution. Since we are running in batch mode, it's important to print the solutions out once they are obtained (either to the screen or to a file or both). The arithmetic statement {getfield /avg/RC[{i}] state} > {thresh} will return either a 0 or a 1, since that's what comparison operators return. If we set i to be equal to all the RC elements from 0 to 63 in a loop, we can print out the entire solution. It's also important to print out a newline after every eight 0's or 1's, so that the appearance of the output is like the chessboard. The following function, "print_solution", does this by using the modulus (remainder) operator, %: function print_solution int i for (i = 0; i < 64; i = i + 1) echo {{getfield /avg/RC[{i}] state} > {thresh}} -n if (((i + 1) % 8) == 0) echo // this prints a blank line end end end If we run this function after the network has settled down into an attractor we will see a display like this: 00001000 00100000 10000000 00000100 00000001 01000000 00010000 00000010 In this case the network has found a solution to the eight queens problem (it doesn't always, as mentioned above). Make a copy of this function, rename it "save_solution" and change it so that the solution is saved to a file rather than printed to the screen. Use the "openfile", "writefile" and "closefile" commands to do this (type "help openfile" etc. for the syntax). Use the file name as an argument to the function. Remember to close the file you open: GENESIS can have no more than 20 files open at any given time. Now we're ready for the main event: defining a function called "run_batch" which will run the simulation in batch mode, test whether the simulation has reached a solution, print a solution to the screen if it finds one, and save all solutions to separate files. We will use the function "count_active_cells", described above, to distinguish solutions from non-solutions. The simplest form of such a function would be something like: function run_batch int solutioncount = 0 int count str filename while (1) // repeat forever reset step 0.5 -t // run for 0.5 seconds count = {count_active_cells} if (count == 8) // we have a solution print_solution filename = "solution." @ {solutioncount} save_solution {filename} solutioncount = solutioncount + 1 end end end This function will run forever; a more reasonable thing to do is to set it up so that it stops after it has found 100 solutions (which should be enough :-). Make this change to "run_batch". Now when you run the simulation in batch mode, it should keep going until 100 solutions have been found, and after leaving GENESIS, the solutions will be in the files "solution.0" to "solution.99". The solution files should look something like this: 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 You can inspect them to ascertain that they are in fact solutions to the eight-queens problem. ---------------------------------------------------------------------- Exercises ---------------------------------------------------------------------- These can be done singly or in combination: 1) Try putting the inhibitory synapses back onto the dendrites (and removing them from the soma) and see if the network can still reach solutions to the eight queens problem, and if so, whether it can reach them as quickly. 2) Lower the decay time constants of the synchans (tau2) to the values in tutorial 7 (3 msec for the excitatory channels, 20 msec for the inhibitory channels). Can you still get convergence to solutions to the eight queens problem? Why or why not? 3) Set up excitatory synapses between the cells in the network that aren't connected by inhibitory synapses. Does the system still reach the desired state(s)? If so, does it reach them faster? Do you have to increase the inhibitory conductances to balance off the increased excitation?