Creating a genome and a cost function

So my goal is to minimize crosscorrelation and autocorrelation, safe for the middle value of autocorrelation. I won’t use a fixed-size sequence, I will encode in the genome the distance between two impulses. This will allow for a really small genome with a lot of possibilities.

For this purpose, I create a class named OrthogonalSequences:

class OrthogonalSequences(object):
  def __init__(self, size, number):
    self.size = size
    self.number = number
    self.genome_length = size * number
 
  @staticmethod
  def convert_slice_to_time(chromosome):
    """"
    Converts a chromosome in a time sequence
    """"
    total = numpy.cumsum(chromosome)
    array = numpy.zeros(total[-1]+1)
    array[total] = 1
    array[0] = 1
    return array
 
  def evaluate(self, chromosomes):
    """"
    Evaluate the cost (square of error) of the genome by evaluating cross- and autocorrelation.
    """"
    value = 0
    for i in range(self.number):
      slice_i = self.convert_slice_to_time(chromosomes[i * self.size:(i+1) * self.size])
      value += .1 * len(slice_i)
      for j in range(i, self.number):
        slice_j = self.convert_slice_to_time(chromosomes[j * self.size:(j+1) * self.size])
        correlation = numpy.correlate(slice_i, slice_j, mode = "full")
        if (i == j):
          correlation[len(slice_i)-1] -= self.size
 
        value += numpy.sum(correlation ** 2)
 
    return value

In fact this class will only be used to evaluate a genome. As I have two combs (and more if you’d like), I need to pass this information to the evaluation function, and thus nothing is more easy than Python for this. I’ve also added a small factor that adds a cost for lengthy combs.

Setting up the optimizer

Now, the only thing remaining is to set all objects and launch the optimization:

if __name__ == "__main__":
  from pyevolve import G1DList, GSimpleGA, Consts
  import sys
 
  size = int(sys.argv[1]) if len(sys.argv) > 1 else 10
  size -= 1
  number = int(sys.argv[2]) if len(sys.argv) > 2 else 2
  rangemax = int(sys.argv[3]) if len(sys.argv) > 3 else 100
  generations = int(sys.argv[4]) if len(sys.argv) > 4 else 50
  population = int(sys.argv[5]) if len(sys.argv) > 5 else 100
 
  genome = G1DList.G1DList(number*size)
  sequences = OrthogonalSequences(size, number)
  ga = GSimpleGA.GSimpleGA(genome)
 
  genome.setParams(rangemin=1, rangemax=rangemax)
  genome.evaluator.set(sequences.evaluate)
 
  ga.setMinimax(Consts.minimaxType["minimize"])
 
  ga.setPopulationSize(population)
  ga.setElitism(True)
  ga.setGenerations(generations)
  ga.evolve(freq_stats = generations/10)
  print ga.bestIndividual()

Results

Now the result is the following for 2 sequences of 10 impulses:

Gen. 1 (2.00%): Max/Min/Avg Fitness(Raw) [436.73(387.20)/317.81(349.20)/363.94(363.94)]
Gen. 5 (10.00%): Max/Min/Avg Fitness(Raw) [417.89(375.20)/340.39(345.20)/348.24(348.24)]
Gen. 10 (20.00%): Max/Min/Avg Fitness(Raw) [419.21(379.20)/339.65(345.20)/349.34(349.34)]
Gen. 15 (30.00%): Max/Min/Avg Fitness(Raw) [417.72(377.20)/341.16(345.20)/348.10(348.10)]
Gen. 20 (40.00%): Max/Min/Avg Fitness(Raw) [419.40(375.20)/337.80(345.20)/349.50(349.50)]
Gen. 25 (50.00%): Max/Min/Avg Fitness(Raw) [419.86(389.20)/341.55(345.20)/349.88(349.88)]
Gen. 30 (60.00%): Max/Min/Avg Fitness(Raw) [418.25(381.20)/341.41(345.20)/348.54(348.54)]
Gen. 35 (70.00%): Max/Min/Avg Fitness(Raw) [418.08(375.20)/340.08(345.20)/348.40(348.40)]
Gen. 40 (80.00%): Max/Min/Avg Fitness(Raw) [418.42(375.20)/339.53(345.20)/348.68(348.68)]
Gen. 45 (90.00%): Max/Min/Avg Fitness(Raw) [418.80(381.20)/340.76(345.20)/349.00(349.00)]
Gen. 50 (100.00%): Max/Min/Avg Fitness(Raw) [417.94(371.20)/338.92(345.20)/348.28(348.28)]
Total time elapsed: 2.478 seconds.
- GenomeBase
        Score:                   345.200000
        Fitness:                 338.919595

        Slot [Evaluator] (Count: 1)
                Name: evaluate
        Slot [Initializator] (Count: 1)
                Name: G1DListInitializatorInteger
                Doc:  Integer initialization function of G1DList

   This initializator accepts the *rangemin* and *rangemax* genome parameters.

        Slot [Mutator] (Count: 1)
                Name: G1DListMutatorSwap
                Doc:  The mutator of G1DList, Swap Mutator
        Slot [Crossover] (Count: 1)
                Name: G1DListCrossoverSinglePoint
                Doc:  The crossover of G1DList, Single Point

   .. warning:: You can't use this crossover method for lists with just one element.

- G1DList
        List size:       18
        List:            [36, 48, 1, 11, 95, 37, 34, 38, 25, 10, 29, 64, 19, 11, 62, 68, 7, 15]

Here are the visual results for this optimization:

The first comb	The second comb
Correlations between the combs

They are not perfect, far from it, but remember that they are only combs with positive impulses. This means that you can’t get a perfect result. An almost-perfect result would be that the crosscorrelation is below 1, but with only one value of 2, it’s still pretty good.

Before the end

The process is very fast for 2 sequences of 10 impulses. If you add more sequences and more impulses, the computation time will rise fastly, but PyEvolve has now the ability to use the multiprocessing module, which means it can do parallel optimizations now.

PyEvolve can do much more, I only barely scratch the surface of its capablities, but I hope this example will make you try this package.

Content and opinions

The book is split in three parts with progress hacking solutions. It is to be mentioned that it is dedicated to hacking the third generation, but the interface should be the same, although some instructions may have changed. Perhaps a new edition of the book may be needed.

iRobot had the good idea of opening their pets to hackers. It is done through a serial interface called ROI. The first encompasses the creation (or the buying) of a connector between a PC and the Roomba. There are several solutions, from a simple serial cable to a Bluetooth interface. Once you have a connector, you need to learn and use the interface to make the robot move and feel. Each captor and motor is explained as well as how it can be interacted with. The main examples are given through a Java library (I would have prefered a Python library, of course).

The second part is perhaps the least interesting. It starts with making an application to steer the robot. It is done with Processing, an obscur language that runs on the JVM. Why not one of the top five script language on it? You may then use the inner beeper to make some music. Indeed, it can be fun, but it’s also not that easy, even with the help of the book. More fun is making art drawings with addition of pens sticked on the Roomba. I was a big fan of these kind of drawngs when I was young. The last chapter exposes the Roomba as a big input device such as a mouse. Of course, it can help doing some sports…

The last part is about using additional platforms to make a more complicated robot. First, you may connect it to the Internet to steer it from the other end of the world, the next step being using WiFi. Then, you may upgrade its brain with a classic Linux board, adding a camera, … OK, this is fun, but it’s not a cleaner anymore. It begins to be even very ugly and not that autonomous unless you stick a big battery pack. But I have to be agree on some aspects of the mod: it’s a good basis for building one’s custom robot. But you may find a skeletton for less (or not).

Conclusion

I intend to keep my Roomba as a vacuum cleaner, not a complete robot with a wifi router on top of it. Still, it is interesting to know that you can do a lot with it, and first of all with a simple serial cable or a Bluetooth connection.

Hacking Roomba
Price: $20.22
By Kurt – John Wiley & Sons, Inc. (00) – Paperback – ISBN 10 0470072717 – ISBN 13 9780470072714

Hacking Roomba: ExtremeTech (Paperback)
by Tod E. Kurt
ISBN:

Price: —
27 used & new available from USD 22.62

| 4.5 | 7

I’ve started a branch available on github.com, and I will some examples in the scikit as well. I may explain them here, but I won’t rewrite what is already published. A future post will explain the changes, and I hope that interested people will understand the modifications and apply them to my former posts. It’s just that I don’t have much time to change everything…

A templated Chamberlin filter

A Chamberlin filter is a IIR (Infinite Impulse Response) filter, and it is possible to make it output high-, band- or low-pass signals. It has only two parameters that are independent.

The equations are very simple:

Each of the series is either the low-pass, the band-pass or the high-pass output.

There are some factors that are pre-computed:

The first is a function of the Cut frequency and the Sampling frequency, a “numerical frequency”. The second one is the numerical attenuation, between 0 and 2.

This is an excerpt of the actual code, I’ve only displayed the relevant code.

template<class Data_Type>
class VariableFilter
{
public:
  typedef Data_Type DataType;
 
  void process(const DataType* in, DataType* out, unsigned long nb_samples)
  {
    for(unsigned long i = 0; i < nb_samples; ++i)
    {
      yh = in[i] - yl - numerical_attenuation * yb;
      yb = numerical_frequency * yh + yb;
      yl = numerical_frequency * yb + yl;
      if(selected == 0)
      {
        out[i] = yl;
      }
      else if(selected == 1)
      {
        out[i] = yb;
      }
      else
      {
        out[i] = yh;
      }
    }
  } 
 
protected:
  void compute_factors()
  {
    numerical_frequency = 2 * std::sin(M_PI * cutoff_frequency / sampling_frequency);
    numerical_attenuation = 2 * attenuation;
  }
};

Now, I’ve written a simple GUI for this plugin.

GUI with Qt

I’ve used the QtAgain skeletton to wrap the filter and to make it available to PyVST. It’s really easy to add other parameters to the skeletton, so I will not write everything here again. There are three parameters, 2 numerical ones, and one three-state variable.

Results

I’ve recorded the filter’s outputs for 6kHz as central frequency, and an attenuation factor of .3. The output signal was generated from a random signal, and then an FFT was performed on the inputs and the outputs and displayed.

Note: I didn’t convert the scale to dB.

As you may have noted, the filter amplifies the signal near the central frequency for each output for the attenuation I’ve selected. The issue is that selecting a higher one has also an impact on the stability of the filter at higher frequencies. The sum of the numerical attenuation and the numerical frequency should always be inferior to 2.

Conclusion

This numerical is really simple to implement, and it has few parameters. The dark side of this filter is that is not always stable and the central frequency seems to be always amplified. Other filters can be better balanced.

The code is available on Launchpad.

scikits.optimization has a special class for handling numerical differentiation, and several tools for conjugate gradients.

Numerical Differentiation

Differentiation is accomplished through helper class. The new function has to be derived from of two classes,
ForwardFiniteDifferences or CenteredFiniteDifferences. In these classes, gradient and Hessian are computed based on the actual cost function. You may also override the gradient so that only the Hessian is numerically computed.

A single parameter, difference, is the delta that will be applied on each parameter to compute the derivatives. This number is arbitrally fixed, and may be changed in numerical differentiation fails.

Here is what a weighted sinc function looks like (the weights are given by self.a):

from scikits.optimization.helpers import ForwardFiniteDifferences
 
class Function(ForwardFiniteDifferences):
  def __init__(self):
    ForwardFiniteDifferences.__init__(self)
    self.a = numpy.array((1, 2))
  def __call__(self, x):
    t = numpy.sqrt(numpy.sum((self.a * x)**2))
    return -numpy.sinc(numpy.pi * t)

Conjugate Gradient

Conjugate gradient is an optimization technique that uses the new gradient as well as the last descent direction or gradient to create a new descent direction. Besides this, the line search must have some properties. You don’t need an exact one, but it must obey some rules, known as the Wolfe-Powell rules. There are two formulations, both available in the scikit. I suggest you read some of the reference book I give at the end of the topic to choose the right line search, as well as the appropriate conjugate gradient.

Besides the well-known Fletcher Reeves conjugate gradient, there are several other factors that can be used, as the Polak-Ribière-Polyak that can also be combined with Fletcher-Reeves one to form what is perhaps the best formulation. The different available conjugate gradients are available in the documentation.

Application

OK, let’s see what it can do. I’ve assembled a Fletcher-Reeves conjugate gradient with Wolfe-Powell rules line search.

mystep = step.FRConjugateGradientStep()
mylinesearch = line_search.WolfePowellRule()
mycriterion = criterion.criterion(ftol = 0.0001, iterations_max = 100)
myoptimizer = optimizer.StandardOptimizer(function = fun,
                                          step = mystep,
                                          line_search = mylinesearch,
                                          criterion = mycriterion,
                                          x0 = numpy.array((.35, .45)))

I didn’t display the result, because in fact it didn’t finished at 100 iterations.

Here are animations on the evolution of the value of the cost function and the position of the parameters.

Knowing that Fletcher-Reeves adds to the new gradient a fraction of the last direction, this circular evolution can be understood. In fact, with a simple gradient, it works better, as well as if a Polak-Ribière-Polyak is used. It all depends on how the different ingredients are mixed, as well as if the Wolfe-Powell rules are close to an exact line search.

Conclusion

I think that trying different solutions is what is fun about optimizations. With conjugate gradient, a lot of different optimizers can be built. With numerical differentitation, you don’t have to rely on analytical gradient. Know though that the more the parameters, the longest it takes to compute the gradient (O(n)) and I don’t even speak about the Hessian (O(n²)).

Here are some useful books that explain the theory and the ideas behind these topics:

Optimization Theory and Methods: Nonlinear Programming

Numerical Optimization

Engineering Optimization: Theory and Practice

Link to the tutorial code: Launchpad.

I’ve been thinking about this for several months now, for work that has thousands of tasks, each task being massively data parallel. Traditionnally, one launches a job through one’s favorite batch scheduler (favorite or mandatory…) with fixed resources and during an estimated amount of time. This may work well in research, but in the industrial world, there often a new job that arises and that needs part of your scarce resources. You may have to stop your work, loose your current advances and/or restart the job with less resources. And then the cycle goes on.

Static resource allocation

How can resource allocation work? Let’s start with a simple case where you have 2 applications with different priorities. One of them has a priority of 70 (it’s supposed to finish in three days) whereas the other one has a priority of 50 (four days left). They share the cluster so that 66% is allocated to the first application and 33% to the second one.

What happens if a third application must be launched with a higher priority, because it has to ne finished by tomorrow? You may stop the other two programs, you may loose a lot of work if you didn’t implement checkpoints (besides, one of them may be an of-the-shelf program you bought yesterday) or suspend it. Either way, this is what you will get:

In fact, even if you use dynamic resource allocation, this is what you must get to have your results by the time you need them, but obviously, you have lost your two other applications. Some batch schedulers allow applications to be suspended, but this is a double-edge sword:

• your cluster must support job suspension, and thus have access to drives to save the job state (which is not possible for medium to large-scaled clusters)
• if your application does not scale to your entire cluster (it happens), although one of the other two applications could go on, it is not possible, all processes are put to sleep

So all things considered, you have to implement dynamic resource allocation.

Dynamic resource allocation

How does this work? Each application must be aware that it can be allocated more resources or deallocated some at all time. To be portable on all clusters, you cannot suspend part of your program, it must really go away. The batch scheduler must also notice that your application has freed some of its resources. You thus have to allocate small jobs that will communicate together (this can be done with MPI-2).

This means that you will have hundreds or thousands of small works. All of them will not have to be connected to the scheduler, only one master must be. Of course, this can easilly be done by using a specific queue. Each application on this queue will thus receive orders from the batch scheduler and act upon it. Another advantage is that also the application gets no resource at one point, it still has a saved state that enable the continuation of a run.

Of course, this is not easy to do. How can this be applied to an of-the-shelf application? Well, in this case, you may create a bogus application on the master queue that will at least allow other applications to be allocated resources beside it.

You do not have to implement this on top of MPI. It can be really hard to do (handling data moves between processors, change the decomposition, …), and you may implement another solution. In my case, I have thousands different tasks that can be run on very few cores, so this is my elementary unit. I don’t need all tasks to communicate between them, so I create each time brand new independent jobs and I also can tell the scheduler it can kill jobs that are not responding before the next allocation phase.

Conclusion

To finish, I’ll say that I know that LSF allows plugins that help dispatch jobs on specific hosts of your cluster (to have the best communication location). There seems to be a way of implementing the needs gathering and the resource assignment, but the documentation is not clear (at all). A specific daemon may be needed. I don’t know if other batch scheduler allow plugins to modify their behavior, if you know of them and their API, please do tell

Content and opinions

There are three book’s authors, three German people, as Germany is one of the countries where building automation is well developped.

The first chapter describes what building automation is. Its sub-categories are also explained (control, measurement and management), and I have to say that there are some differences that I didn’t apprehend before. Automation has several clear benefits: comfort with automatic heating, lights, … but also economies as a good heating control will lead to less consumption. This can really help in commercial or industrial buildings. It’s also astonishing to see that even new bulding do not come with a complete automation system. Heating is almost handled, but lights and blinds are not. Heating is not completely handled in my opinion because settings cannot be given outside your office.

The second chapter is perhaps the less interesting of the five. If you have basis in telecommunications, it will be a rehearsal of your first class. If you have not, you will at least know how it can work.

Now, the book has three chapters on each of the main automation technologies, Konnex, LON and BACnet. Each time, the chapter is independent, and sometimes you even learn something on a technology when reading the next one. When I say independent, I mean that you don’t need to read the other chapters if you want to read only one. It also means that you will see redundancy. This could be explained by each author writting his own chapter without much interaction on their content.

Each time, you will learn every detail of the different busses, with perhaps sometimes too many. For BACnet, you will go through the Ethernet interface in detail, although it is something really common in the ISO model. But to select the right bus, you need to know everything it can provide.

Conclusion

Contrary to my last book, this one was really interesting. I didn’t know that KNX was mainly based on one BCU (Bus Coupling Unit) with different actuators. Although you can also do management with KNX, I concur that it is a control bus at heart. I also appreciate the fact that each bus is open in the sense that there are several manufacturers that can provide appliances, and that their operation is clearly explained.

If you don’t know what building automation is and if you want to build your home, a commercial building, or if you renovate something, this book can help you should an adequate automation bus. Don’t forget that in these times of peak energy, economy associated with better comfort is a horse that can win the race.

Building Automation: Communication Systems with Eib/Knx, Lon Und Bacnet
Price: $87.2

Building Automation: Communication systems with EIB/KNX, LON and BACnet (Signals and Communication Technology) (Hardcover)
by Hermann Merz, Thomas Hansemann, Christof Hübner
ISBN: 3540888284

Price: USD 87.20
27 used & new available from USD 81.05

| | 0

Content and opinions

The book starts with two protocols that can be used, X10 and C-Bus. X10 uses the power lines to communicate between elements, and C-Bus uses a proprietary protocol on a specific bus that must be installed by a specific company. It can also use a wireless protocol. Both seem to be widely available in America, but in Europe, it’s different. The chapter presents a lot of appliances, but I don’t know f they are available on Europe power lines for X10 or even if they are sold in Europe (for both of them). It’s too bad the author didn’t choose a better support (like a real ISO and open standard like KNX). X10 is somewhat limited like C-Bus is although in a different way.

The second chapter is about different hardwares that can be used. The first is the well-known NSLU2, then are other gaming consoles. It seems that safe for the high end consoles, every other hardware is either discontinued or difficult to hack (there are a lot of hacking to do according to the book). For instance Arduino needs some soldering. What about hardware like the Sheeva plug ?

The next chapter explains how media data can be stored and accessed. It’s more or less a sum up of what can be found on the Internet. In fact, the rest of the book is material like this. A real big sum up of what can be found on media use on the Internet, and a lot of shell code to make it work.

Conclusion

It’s not that the book is of no use. If you have X10 appliances or C-Bus ones, if you need a big tutorial on media control, this book might be for you. If you want to do home automation, I’m afraid only the first chapter deals with what home automation is mainly about: light and heating control.

Smart Home Automation with Linux
Price: $25.19

Smart Home Automation with Linux (Expert’s Voice in Linux) (Paperback)
by Steven Goodwin
ISBN: 1430227788

Price: USD 23.09
33 used & new available from USD 20.51

| 5 | 4

Content and opinions

Each “thing” is explained by a two-pages article. The tone is different as each chapter is written by different authors. They are sorted alphabetically, and there is a table of contents sorted by topics.

The topics browse the modern/current approach to programming (in almost every aspect), emphasing the need for tests, automatic tools, good designs or personal enhancement. Of course, they cannot be deepened in so few words, so you have to get another book for more detailed explanations. This also means that for some pieces of advice, you may found them too obvious, as you need details to implement the approach.

Conclusion

This book is a great first approach to good practices. If you are used to reading good practices books, you may not need this one, but it tries efficiently to sum up the different faces of modern programming.

97 Things Every Programmer Should Know: Collective Wisdom from the Experts
Price: $21.59

97 Things Every Programmer Should Know: Collective Wisdom from the Experts (Paperback)
by
ISBN: 0596809484

Price: USD 19.79
44 used & new available from USD 15.04

| 3.5 | 10

Setting up the optimizer and the cost function

I’ve refactored the cost function and I’ve added a gradient method:

class Function(object):
    def __call__(self, x):
        t = numpy.sqrt(numpy.sum(x**2, axis=0))
        return -numpy.sinc(t)
    def gradient(self, x):
        t = numpy.sqrt(numpy.sum(x**2, axis=0))
        return x / t**2 * (- numpy.cos(numpy.pi*t) + numpy.sinc(t))

This is how the scikit works. All cost functions are instance of this kind of class. The __call__ method returns the actual cost, and additional methods provide gradient or hessian. Here, only the gradient is required.

Now, I can build the optimizer:

from scikits.optimization import *
 
mystep = step.GradientStep()
mylinesearch = line_search.GoldenSectionSearch(min_alpha_step = 0.0001, alpha_step = .5)
mycriterion = criterion.criterion(ftol = 0.0000001, iterations_max = 10)

As I’ve said before, the scikit is built around the separation principle, so what I’ve created does not contain any state. Now, I can create the glue that will steer the optimization:

myoptimizer = optimizer.StandardOptimizer(function = fun,
                                          step = mystep,
                                          line_search = mylinesearch,
                                          criterion = mycriterion,
                                          x0 = numpy.array((.8, 1.)))

As for the polytope optimizer, I can begin the optimization by calling the optimize() method.

Conclusion

This is the evolution of the cost function during the optimization:

Cost function values for each iteration

This was a simple usage of a gradient optimizer. The function is of course easy to optimize, as it has a circular symmetry. Another problem is that you need to have access to an analytical gradient. In further tutorials, I’ll introduce numerical approximation classes that will help with this task.

As usual, the code of this tutorial is available on Launchpad.