Experimenting with Genetic Algorithms

September 13, 2016

Genetic algorithm trying to find a series of four mathematical operations (e.g. -3*4/7+9) that would result in the number 42.

Genetic algorithm trying to find a series of four mathematical operations (e.g. -3*4/7+9) that would result in the number 42.

I’m teaching a numerical methods class that’s partly an introduction to programming, and partly a survey of numerical solutions to different types of problems students might encounter in the wild. I thought I’d look into doing a session on genetic algorithms, which are an important precursor to things like networks that have been found to be useful in a wide variety of fields including image and character recognition, stock market prediction and medical diagnostics.

The ai-junkie, bare-essentials page on genetic algorithms seemed a reasonable place to start. The site is definitely readable and I was able to put together a code to try to solve its example problem: to figure out what series of four mathematical operations using only single digits (e.g. +5*3/2-7) would give target number (42 in this example).

The procedure is as follows:

  • Initialize: Generate several random sets of four operations,
  • Test for fitness: Check which ones come closest to the target number,
  • Select: Select the two best options (which is not quite what the ai-junkie says to do, but it worked better for me),
  • Mate: Combine the two best options semi-randomly (i.e. exchange some percentage of the operations) to produce a new set of operations
  • Mutate: swap out some small percentage of the operations randomly,
  • Repeat: Go back to the second step (and repeat until you hit the target).

And this is the code I came up with:

genetic_algorithm2.py

''' Write a program to combine the sequence of numbers 0123456789 and
    the operators */+- to get the target value (42 (as an integer))
'''

'''
Procedure:
    1. Randomly generate a few sequences (ns=10) where each sequence is 8
       charaters long (ng=8).
    2. Check how close the sequence's value is to the target value.
        The closer the sequence the higher the weight it will get so use:
            w = 1/(value - target)
    3. Chose two of the sequences in a way that gives preference to higher
       weights.
    4. Randomly combine the successful sequences to create new sequences (ns=10)
    5. Repeat until target is achieved.

'''
from visual import *
from visual.graph import *
from random import *
import operator

# MODEL PARAMETERS
ns = 100
target_val = 42 #the value the program is trying to achieve
sequence_length = 4  # the number of operators in the sequence
crossover_rate = 0.3
mutation_rate = 0.1
max_itterations = 400


class operation:
    def __init__(self, operator = None, number = None, nmin = 0, nmax = 9, type="int"):
        if operator == None:
            n = randrange(1,5)
            if n == 1:
                self.operator = "+"
            elif n == 2:
                self.operator = "-"
            elif n == 3:
                self.operator = "/"
            else:
                self.operator = "*"
        else:
            self.operator = operator
            
        if number == None:
            #generate random number from 0-9
            self.number = 0
            if self.operator == "/":
                while self.number == 0:
                    self.number = randrange(nmin, nmax)
            else:
                self.number = randrange(nmin, nmax)
        else:
            self.number = number
        self.number = float(self.number)

    def calc(self, val=0):
        # perform operation given the input value
        if self.operator == "+":
            val += self.number
        elif self.operator == "-":
            val -= self.number
        elif self.operator == "*":
            val *= self.number
        elif self.operator == "/":
            val /= self.number
        return val


class gene:

    def __init__(self, n_operations = 5, seq = None):
        #seq is a sequence of operations (see class above)
        #initalize
        self.n_operations = n_operations
        
        #generate sequence
        if seq == None:
            #print "Generating sequence"
            self.seq = []
            self.seq.append(operation(operator="+"))  # the default operation is + some number
            for i in range(n_operations-1):
                #generate random number
                self.seq.append(operation())

        else:
            self.seq = seq

        self.calc_seq()

        #print "Sequence: ", self.seq
    def stringify(self):
        seq = ""
        for i in self.seq:
            seq = seq + i.operator + str(i.number)
        return seq

    def calc_seq(self):
        self.val = 0
        for i in self.seq:
            #print i.calc(self.val)
            self.val = i.calc(self.val)
        return self.val

    def crossover(self, ingene, rate):
        # combine this gene with the ingene at the given rate (between 0 and 1)
        #  of mixing to create two new genes

        #print "In 1: ", self.stringify()
        #print "In 2: ", ingene.stringify()
        new_seq_a = []
        new_seq_b = []
        for i in range(len(self.seq)):
            if (random() < rate): # swap
                new_seq_a.append(ingene.seq[i])
                new_seq_b.append(self.seq[i])
            else:
                new_seq_b.append(ingene.seq[i])
                new_seq_a.append(self.seq[i])

        new_gene_a = gene(seq = new_seq_a)
        new_gene_b = gene(seq = new_seq_b)
                       
        #print "Out 1:", new_gene_a.stringify()
        #print "Out 2:", new_gene_b.stringify()

        return (new_gene_a, new_gene_b)

    def mutate(self, mutation_rate):
        for i in range(1, len(self.seq)):
            if random() < mutation_rate:
                self.seq[i] = operation()
                
            
            

def weight(target, val):
    if val <> None:
        #print abs(target - val)
        if abs(target - val) <> 0:
            w = (1. / abs(target - val))
        else:
            w = "Bingo"
            print "Bingo: target, val = ", target, val
    else:
        w = 0.
    return w

def pick_value(weights):
    #given a series of weights randomly pick one of the sequence accounting for
    # the values of the weights

    # sum all the weights (for normalization)
    total = 0
    for i in weights:
        total += i

    # make an array of the normalized cumulative totals of the weights.
    cum_wts = []
    ctot = 0.0
    cum_wts.append(ctot)
    for i in range(len(weights)):
        ctot += weights[i]/total
        cum_wts.append(ctot)
    #print cum_wts

    # get random number and find where it occurs in array
    n = random()
    index = randrange(0, len(weights)-1)
    for i in range(len(cum_wts)-1):
        #print i, cum_wts[i], n, cum_wts[i+1]
        if n >= cum_wts[i] and n < cum_wts[i+1]:
            
            index = i
            #print "Picked", i
            break
    return index

def pick_best(weights):
    # pick the top two values from the sequences
    i1 = -1
    i2 = -1
    max1 = 0.
    max2 = 0.
    for i in range(len(weights)):
        if weights[i] > max1:
            max2 = max1
            max1 = weights[i]
            i2 = i1
            i1 = i
        elif weights[i] > max2:
            max2 = weights[i]
            i2 = i

    return (i1, i2)
        
    
    
            


# Main loop
l_loop = True
loop_num = 0
best_gene = None

##test = gene()
##test.print_seq()
##print test.calc_seq()

# initialize
genes = []
for i in range(ns):
    genes.append(gene(n_operations=sequence_length))
    #print genes[-1].stringify(), genes[-1].val
    

f1 = gcurve(color=color.cyan)

while (l_loop and loop_num < max_itterations):
    loop_num += 1
    if (loop_num%10 == 0):
        print "Loop: ", loop_num

    # Calculate weights
    weights = []
    for i in range(ns):
        weights.append(weight(target_val, genes[i].val))
        # check for hit on target
        if weights[-1] == "Bingo":
            print "Bingo", genes[i].stringify(), genes[i].val
            l_loop = False
            best_gene = genes[i]
            break
    #print weights

    if l_loop:

        # indicate which was the best fit option (highest weight)
        max_w = 0.0
        max_i = -1
        for i in range(len(weights)):
            #print max_w, weights[i]
            if weights[i] > max_w:
                max_w = weights[i]
                max_i = i
        best_gene = genes[max_i]
##        print "Best operation:", max_i, genes[max_i].stringify(), \
##              genes[max_i].val, max_w
        f1.plot(pos=(loop_num, best_gene.val))


        # Pick parent gene sequences for next generation
        # pick first of the genes using weigths for preference    
##        index = pick_value(weights)
##        print "Picked operation:  ", index, genes[index].stringify(), \
##              genes[index].val, weights[index]
##
##        # pick second gene
##        index2 = index
##        while index2 == index:
##            index2 = pick_value(weights)
##        print "Picked operation 2:", index2, genes[index2].stringify(), \
##              genes[index2].val, weights[index2]
##

        (index, index2) = pick_best(weights)
        
        # Crossover: combine genes to get the new population
        new_genes = []
        for i in range(ns/2):
            (a,b) = genes[index].crossover(genes[index2], crossover_rate)
            new_genes.append(a)
            new_genes.append(b)

        # Mutate
        for i in new_genes:
            i.mutate(mutation_rate)
                

        # update genes array
        genes = []
        for i in new_genes:
            genes.append(i)


print
print "Best Gene:", best_gene.stringify(), best_gene.val
print "Number of iterations:", loop_num
##

When run, the code usually gets a valid answer, but does not always converge: The figure at the top of this post shows it finding a solution after 142 iterations (the solution it found was: +8.0 +8.0 *3.0 -6.0). The code is rough, but is all I have time for at the moment. However, it should be a reasonable starting point if I should decide to discuss these in class.

Citing this post: Urbano, L., 2016. Experimenting with Genetic Algorithms, Retrieved July 27th, 2017, from Montessori Muddle: http://MontessoriMuddle.org/ .
Attribution (Curator's Code ): Via: Montessori Muddle; Hat tip: Montessori Muddle.

Limiting Chemical Reactions

January 1, 2015

Figuring out the limiting reactant in a chemical reaction integrates many of the basic chemistry concepts including: unit conversions from moles to mass and vice versa; the meaning of chemical formulas; and understanding the stoichiometry of chemical reactions. So, since we’ll need a number of these, I wrote a python program to help me design the questions (and figure out the answers).

Program examples come zipped because they require the program file and the elements_database.py library:

Baking Powder and Vinegar (Common Molecules)

Limiting_component-Common.py: This has the baking powder and vinegar reaction limited by 5 g of baking soda. It’s nice because it uses a few pre-defined “common molecules” (which are defined in the elements_database.py library.

You enter the reactants and products and the program checks if the reaction is balanced, then calculates the moles and masses based on the limiting component, and finally double checks to make sure the reaction is mass balanced.

Limiting_component-Common.py

from elements_database import *
import sys

print "LIMITING REACTANT PROGRAM"
print 
print "  Determines the needed mass and moles of reactants and products if reaction is limited by one of the components"

c = common_molecules()

'''Create Reaction'''
rxn = reaction()
# Add reactants (and products)
#   Use rxn.add_reactant(molecule[, stoichiometry])
#       molecule: from molecule class in elements_database
#       stoichiometry: integer number
rxn.add_reactant(c.baking_soda,1)
rxn.add_reactant(c.hydrochloric_acid,1)
rxn.add_product(c.carbon_dioxide)
rxn.add_product(c.water, 1)
rxn.add_product(c.salt)

'''Print out the reaction'''
print 
print "Chemical Formula"
print "  " + rxn.print_reaction()
print 

'''Check if reaction is balanced'''
balanced = rxn.check_for_balance(printout=True)

'''Calculate limits of reaction'''
if balanced:
    rxn.limited_by_mass(c.baking_soda, 5, printout=True)

Outputs results in the Results table (using scientific notation):

LIMITING REACTANT PROGRAM

  Determines the needed mass and moles of reactants and products if reaction is limited by one of the components

Chemical Formula
  NaHCO3 + HCl  --> CO2 + H2O + NaCl 

Check for balance
 --------------------------------------------- 
| Element | Reactants | Products | Difference |
 --------------------------------------------- 
|   Na    |    1      |    -1    |     0      |
|   H     |    2      |    -2    |     0      |
|   C     |    1      |    -1    |     0      |
|   O     |    3      |    -3    |     0      |
|   Cl    |    1      |    -1    |     0      |
 --------------------------------------------- 
Balance is:  True

Given: Limiting component is 5 g of NaHCO3.
  Molar mass = 84.00676928
  Moles of NaHCO3 = 0.0595190130849

 Results
   ------------------------------------------------------------------ 
  | Molecule | Stoich.*| Molar Mass (g) | Moles req. |    Mass (g)   |
   ------------------------------------------------------------------ 
  |NaHCO3    |    1    |    84.0068     |  5.952e-02 |   5.000e+00   |
  |HCl       |    1    |    36.4602     |  5.952e-02 |   2.170e+00   |
  |CO2       |    -1   |    44.0096     |  5.952e-02 |   2.619e+00   |
  |H2O       |    -1   |    18.0156     |  5.952e-02 |   1.072e+00   |
  |NaCl      |    -1   |    58.4418     |  5.952e-02 |   3.478e+00   |
   ------------------------------------------------------------------ 
 * negative stoichiometry means the component is a product

  Final Check: Confirm Mass balance: 
     Reactants:    7.1701 g 
     Products:    -7.1701 g 
     --------------------------
          =      0.0000 g 
     --------------------------

General Example

If not using the common molecules database, you need to define the components in the reaction as molecules yourself. This example reacts magnesium sulfate and sodium hydroxide, and limits the reaction with 20 g of magnesium sulfate.

Limiting_component-General.zip

The main file is:

Chem_Exam-limiting.py

from elements_database import *
import sys

print "LIMITING REACTANT PROGRAM"
print 
print "  Determines the needed mass and moles of reactants and products if reaction is limited by one of the components"

# create reaction

rxn2 = reaction()
# Add reactants (and products)
#   Use rxn.add_reactant(molecule[, stoichiometry])
#       molecule: from molecule class in elements_database
#       stoichiometry: integer number
rxn2.add_reactant(molecule("Mg:1,S:1,O:4"))
rxn2.add_reactant(molecule("Na:1,O:1,H:1"), 2)
rxn2.add_product(molecule("Mg:1,O:2,H:2"))
rxn2.add_product(molecule("Na:2,S:1,O:4"))

'''Print out the reaction'''
print 
print "Chemical Formula"
print "  " + rxn2.print_reaction()
print 

'''Check if reaction is balanced'''
balanced = rxn2.check_for_balance(printout=True)

'''Calculate limits of reaction'''
if balanced:
    rxn2.limited_by_mass(molecule("Mg:1,S:1,O:4"), 20, True)

# print out masses
print "Masses Involved in Reaction"
print "  Reactants:"
for i in rxn2.reactants:
    #print "rxn", i
    print "    {m}: {g} g".format(m=i.molecule.print_formula().ljust(10), g=i.mass)
print "  Products:"
for i in rxn2.products:
    #print "rxn", i
    print "    {m}: {g} g".format(m=i.molecule.print_formula().ljust(10), g=-i.mass)

This program is slightly different from the common molecules example in that, at the end, it prints out masses calculated in a more easily readable format in addition to the other data.

Masses Involved in Reaction
  Reactants:
    MgSO4     : 20.0 g
    NaOH      : 13.2919931774 g
  Products:
    MgO2H2    : 9.69066512026 g
    Na2SO4    : 23.6013280571 g

When I have some time I’ll convert this to JavaScript like the molecular mass example.

Citing this post: Urbano, L., 2015. Limiting Chemical Reactions, Retrieved July 27th, 2017, from Montessori Muddle: http://MontessoriMuddle.org/ .
Attribution (Curator's Code ): Via: Montessori Muddle; Hat tip: Montessori Muddle.

How Long does it Take to Make a Vpython Program?

January 25, 2012

Moving the magnet through a wire coil creates an electric current in the wire. Animation captured from the VPython program: Magnetic Induction - Coils.

My students asked me this question the other day, and while slapping together an animation of electromagnetic induction I gave it some thought.

This program itself is really simple. It took about 15 minutes.

But that’s not counting the half hour I spent searching the web for an image I could use to illustrate magnetic induction and not finding one I could use.

Nor does it count the four hours I spent after I got the animation working to get the program to take screen captures automatically. Of course, I must admit that figuring out the screen captures would have gone a lot quicker if I’d not had to rebuild all my permissions on my hard drive (I’d recently reformatted it), and reinstall ImageMagick and gifsicle to take the screen captures and make animations.

Citing this post: Urbano, L., 2012. How Long does it Take to Make a Vpython Program?, Retrieved July 27th, 2017, from Montessori Muddle: http://MontessoriMuddle.org/ .
Attribution (Curator's Code ): Via: Montessori Muddle; Hat tip: Montessori Muddle.

Algebra in Programming, and a First Box

September 6, 2011

One of the first things you learn in algebra is to use variables to represent numbers. Variables are at the heart of computer programming, and in Python you can use them for more than just numbers. So open the IDLE text editor and get to work.

But first we’ll start with numbers. To assign a number to a variable you just write an equation. Here are a couple (you can copy and paste the code into the IDLE window, but usually typing it in yourself tends to be more meaningful and help you remember):

a = 2
b = 3

Now we can add these two variables together and create a new variable c.

a = 2
b = 3
c = a + b

Which is all very nice, but now we want our program to actually tell us what the result is so we print it to the screen.

a = 2
b = 3
c = a + b
print c

Run this program (F5 or select “Run Module” in the “Run” menu) to get:

Output from the first program: 2+3=5.

Basic Operations and Types of Numbers

Now try some other basic operations:

+ and – are easy as you’ve seen above.

× : for multiplication use *, as in:

a = 2
b = 3
c = a * b
print c

÷ : to divide use a /, as in:

a = 2
b = 3
c = a / b
print c

Now as you know, 2 divided by 3 should give you two thirds, but the running this program outputs 0:

2 divided by 3 to gives 0 because Python thinks we're working with integers.

This is because the Python thinks you’re using integers (a whole number), so it gives you an integer result. If you want a fraction in your result, you need to indicate that you’re using real numbers, or more specifically, rational numbers, which can be integers or fractions, but usually show up as a decimal (these are usually referred to as floating point numbers in programming).

The easiest way to indicate that you don’t just want integers is to make one of your original numbers a decimal:

a = 2.0
b = 3
c = a / b
print c

which produces:

Use a = 2.0 to indicate that we're using rational numbers.

Other Things Can Be Variables

In object oriented programming languages like Python you can assign all sorts of things to variables, not just numbers.

To create a box and give it a variable name you can use the program:

from visual import *
c = box()

which produces:

A box created with VPython. It looks much more interesting if you rotate it to see it in perspective.

A rotated, zoomed-out view of a box.

To rotate the view, hold down and drag the right mouse button. To zoom in or out, hold down the right and left buttons together and drag in and out. Mac users will probably have to use the “option” button to zoom, and the “command” button to rotate.

The line from visual import * tells the computer that it needs to use all the stuff in the module called “visual”, which has all the commands to make 3d objects (the “*” indicates all).

The c = box() line creates the box and assigns it a variable name of c. You don’t just have to use letters as variable names. In programming you want to use variable names that will remind you of what it’s supposed to represent. So you could just as well have named your variable “mybox” and gotten the same result:

from visual import *
mybox = box()

Now objects like this box have properties, like color. To make the box red you set the color property in one of two ways. The first method is to set the color as you create the object:

from visual import *
mybox = box(color = color.red)

The second is to set the property using the variable you’ve created and “dot” (.) notation.

from visual import *
mybox = box()
mybox.color = color.red

In both of these color.red is a variable name that the computer already knows because it was imported when you imported the “visual” module. There are a few other named colors like color.green and color.blue that you can find out more about in the VPython documentation (specifically here).

You can also find out about the other properties boxes have, like length, width and position (pos), as well as all the other objects you can create, such as springs, arrows and spheres.

At this point, its probably worth spending a little time creating new objects, and varying their properties.

Overview

We’ve just covered:

  • Assigning values to variables
  • Basic operations (+,-,×, and ÷)
  • Types of Numbers: Integers versus floating point
  • Assigning 3d objects to variables
  • Setting properties of 3d objects

Next we’ll try to make things move.

References

For how to install and run VPython, check here.

Citing this post: Urbano, L., 2011. Algebra in Programming, and a First Box, Retrieved July 27th, 2017, from Montessori Muddle: http://MontessoriMuddle.org/ .
Attribution (Curator's Code ): Via: Montessori Muddle; Hat tip: Montessori Muddle.

Algebra and Programming with VPython

September 4, 2011

Computer programming is the place where algebra comes to life. Students seem to get really excited when they write even the simplest instructions and see the output on the screen. I’m not sure exactly why this is the case, but I suspect it has something to do with being able to see the transition from abstract programming instructions to “concrete” results.

So I’ve decided to supplement my Algebra classes with an introduction to programming. I’m using the Python programming language, or, more specifically, the VPython variant of the language.

Why VPython? Because it’s free, it’s an easy-to-use high-level language, and it’s designed for 3d output, which seems to be somewhat popular these days. The oohs and aahs of seeing the computer print the result of a+b turn into wows when they create their first box. I’ve used the language quite a bit myself, and there are a lot of other interesting applications available if you search the web.

VPython was created to help with undergraduate physics classes, but since it was made to be usable by non-science majors, it’s really easy for middle and high school students to pick up. In fact, NCSU also has a distance education course for high school physics teachers. They also have some instructional videos available on YouTube that provide a basic introduction.

Image from a game created by middle school student Ryan W.

I use VPython models for demonstrations in my science classes, I’ve had middle school students use it for science projects, and I’ve just started my middle school algebra/pre-algebra students learning it as a programming language and they’re doing very well so far.

What I hope to document here is the series of lessons I’m putting together to tie, primarily, into my middle school algebra class, but should be useful as a general introduction to programming using VPython.

Getting VPython

You’ll need to install Python and VPython on your system. They can be directly downloaded from the VPython website’s download page for Windows, Macintosh or LINUX.

Running a Python program.

Once they’re installed, you’ll have the IDLE (or VIDLE) program somewhere on your system; a short-cut is usually put on the desktop of your Windows system. Run (double-click) this program and the VPython programming editor will pop up any you’re ready to go. You can test it by typing in something simple like:

a = 1
b = 2
c = a + b
print c

Then you run the program by going through the Run–>Run Module in the menu bar.

Which should cause a new window to pop up with:

Python 2.7.1 (r271:86882M, Nov 30 2010, 09:39:13) 
[GCC 4.0.1 (Apple Inc. build 5494)] on darwin
Type "copyright", "credits" or "license()" for more information.
>>> ================================ RESTART ================================
>>> 
3
>>> 

Even better might be to test the 3d rendering, which you can do with the following program:

from visual import *

box()

which creates the following exciting image:

A box created with VPython. It looks much more interesting if you rotate it to see it in perspective.

To rotate the view, hold down and drag the right mouse button. To zoom in or out, hold down the right and left buttons together and drag in and out.

A rotated, zoomed-out view of a box.

Lessons

Lesson 1: Variables, Basic Operations, Real and Integer Numbers and the First Box.

Lesson 2: Creating a graphical calculator: Coordinates, lists, loops and arrays: A Study in Linear Equations

Citing this post: Urbano, L., 2011. Algebra and Programming with VPython, Retrieved July 27th, 2017, from Montessori Muddle: http://MontessoriMuddle.org/ .
Attribution (Curator's Code ): Via: Montessori Muddle; Hat tip: Montessori Muddle.

Creative Commons License
Montessori Muddle by Montessori Muddle is licensed under a Creative Commons Attribution-Noncommercial-Share Alike 3.0 United States License.