Content and opinions

The first chapter introduces the reader to the scientific ecosystem. The main modules are covered, but there are some small mistakes (ipython –pylab does not import Matplotlib, Numpy ans Scipy, just a subset of the first that gives a Matlab-like interface aith some renamed Numpy functions) and what I think is a bad habit (importing everything from Pylab and using it as is). Nothing is lost because in the second chapter, Numpy is properly introduced and imported explicitly. There is a link missing between the two chapters because I didn’t understand why Numpy was used this way and not with pylab.

So the second chapter is about the core Numpy data structure: the multidimensional array. The author browses through different ways of creating them, by stacking them, flattening them… The next chapter deals with universal functions, or put it in an other way, functions that run on array element by elements. There are many different way to do so, with specifics that are tackled properly.

Before the chapter of matrices, there is a useful chapter on correlation, convolution and polynomials. Although one may want to go up to Scipy for signal processing, more often than not, Numpy is enough. There enough examples to understand how everything work. Then, the much hated matrix class is introduced. There are many discussions online on whether this class is actually useful, and I won’t delve on this. Suffice to say that the power of this class can be seen in the examples.

The following chapter is about the different submodules inside Numpy: linear algebra, fft, random number. The proper pointers to Scipy are provided, as well as explicit samples. I only can regret the time shifting example is not perfect, as a filter is applied on the amplitude and the inverse FFT is applied on the amplitudes only. So the transformed signal loses all the phase information, which may be why the image is similar but not that similar to the original one.

The seventh chapter tackles different extra functions, mainly finance-related, and I have to say that I don’t know their use enough to comment. Of course, this is why they are in the middle of the book and not in the first pages as the other Numpy functions. Still, there are some useful functions here (some I didn’t know about), as sorting, searching…

When one code scientific code, one often forget about testing. And Numpy has a nice module for scientific testing. It is nice to know that this aspect of science is not forgotten here and has a proper introduction. (also don’t forget about version control!)

The last three chapters introduces additional modules. The first one was partial addressed before, Matplotlib, and if you need something more advanced, there is always the Matplotlib book also published by Packt. Then there are some examples with Scipy, Scikits (soon a new book on Machine Learning with scikit-learn will be available, also by Packt, for which I was a technical reviewer, and it is really great) and other tools. The final chapter is about Pygame, but I don’t code games 100% in Python, so I didn’t really read it!

Conclusion

It’s hard to be mad at the author for the import issue. But I find it also difficult not to, as the philosophy changes depending on the chapter without saying why. Still, for an introductory book to Numpy, it is great if not excellent. A lot of simple examples, a lot of checks, the good pointers to write efficient code…

Changelog

• Fixes on the Simplex state implementation
• Added several Quasi-Newton steps (BFGS, rank 1 update…)

The scikit can be installed with pip/easy_install or downloaded from PyPI

Old announces:

• 0.2
• 0.1

Optimizers

What is great with the Rosenbrock cost function can be summed up in a few points:

1. It is hard to optimize
2. Gradient can be easily computed

I’ve decided to compare the number of function and gradient calls as well as the cost behavior for several usual optimization algorithms. So the contestants will be:

• A Nelder-Mead/Polytope/simplex optimizer
• SSA, a simplex with simulated annealing (think of amotsa from Numerical Recipes)
• GA, a genetic algorithm
• a gradient-based optimizer
• CG, a conjugate-gradient optimizer
• BFGS, a quasi-Newton optimizer

The first 4 are from scikits.optimization, the last 2 are based on proprietary code that cannot be published, but it’s interesting to compare with other tools that are used to compare gradient-free complex cost functions.

Results

I’ve made a small slideshow with the derivative-free algorithms. First you have for each of the three algorithms the number of function calls versus iteration, the cost versus iteration and finally the location of testing parameters.
Click to view slideshow.

This slideshow is for the derivative-based algorithms.
Click to view slideshow.

Conclusion

I was quite surprised by some algorithm behaviors. Clearly, the conjugate gradient algorithm behaves far better than the simple gradient, but the BFGS followed the Rosenbrock valley far better. A good quasi-Newton can be really efficient (not a brent because it needs to solve a linear equation), although a conjugate gradient can be enough in some cases.

For the gradient-free algorithms, SSA really behaved badly. This is mainly due because the hyperparameters that must be adequately tuned. This function is quite simple, but my first trial at setting these parameters was far more efficient for GA or the simplex than for SSA. So I would go for GA for gradient-free optimization: few and easy hyper parameters and a good browse of the search space.

The code for the 4 first tests and display plots

Content and opinions

Numpy Cookbook is a list of several dozens recipes to use Numpy. Or so the title and subtitle say, but in fact it’s nothing like that. More or less than 75% of the book are not about Numpy, but about the Numpy ecosystem, which is far better. This also means that the audience for the book is not Python beginners, but people who know about how the language works and are looking for advice on the scientific side.

The first chapter tackles IPython and more precisely a new feature called the notebook. The specific scientific side of it is the support of Matplotlib, the notebook being a good way of sharing scripts with plots between people. It is indeed a common problem that one may encounter. The next two chapters cover Numpy aspects: how arrays, indexing work and how functions can be applied on them. It is worth mentioning that the last recipe (the Sieve of Erathosthenes) from these chapters doesn’t work. I actually don’t know how the author expected it to work, the algorithm being wrong. There are also several installation recipes, but they tend to be slightly different each time. Why? No real explanation. This happens a lot with all the different installations. Sometimes there are even packages that don’t have such recipe. So no coherency here.

The next chapter handles communication between other Python modules, other languages or frameworks. A small overview of what is called the buffer protocol and the array interface is provided. They are both sides of the same coin, the inner Numpy way of working. Then other languages are covered: Matlab and Octave (but it seems that the author forgot going back from Matlab to Numpy and the fact that Matlab changed its I/O format and went to HDF5), R, Java and the Google App Engine. Mainly these recipes are about installing the interfaces, not about using them and then, there is only a small example of using them. I would have expecting something deeper. One still has to read the manual of these interfaces to use them even for something very basic.

Next two chapters are actual scientific recipes. Based on 1D and 2D signals, the first one tackles some modules in Scipy in addition to memmap, a way of using data from the disk without string it in memory. Unfortunately, this is not dealt in-depth, memmap may have needed more that just a few lines. The second chapter presents more Numpy classes that can be useful, especially masked arrays, as well as universal functions. It is advanced Numpy use and regularly helpful in complex algorithms.

Before the last chapter on Scikits, the author spends two chapters on something scientists often forget: profiling, debugging and quality insurance. Almost all the main modules in these domains are introduced, safe for nosetests, currently the main test framework. It’s a shame that it isn’t one of the recipes, as it supersedes all the modules of the tests recipes. Still, I really appreciate the presence of all these tips that all scientists should apply as they are too often overlooked on scientific books.

As I’ve said, the last chapter presents some scikits, application specific packages (machine learning, image), and pandas, a statistics module. I think of them as threads that readers can follow if they want to dig in a specific direction.

Conclusion

Let’s start with saying that the book is neither a good nor a bad book. It could have been a good book if the author had gone further into the content, and more precisely Numpy. My point of view is that the title, subtitle and everything on the cover cannot represent the book. It doesn’t solve common problems (finding prime number or a palindrome is not a common problem, it’s a pet problem but not an actual real-life problem) and too many recipes are installation recipes. One would have sufficed.

But the book is not a bad one. If the reader doesn’t know much of Numpy, it can help starting using it and as I’ve said, threads can be followed if a specific direction is of interest. It can also a sort of exercise book for a Numpy crash course, provided that some errors are fixed.

Final word: even if I would have liked more structure is the code (class and function), at least the infamous from pylab import * is not advocated in the book.

When the only thing you have is the cost function and when you don’t have dozens of parameters, the first thing that can be tried is a simplex algorithm.

Usage

The only specific thing you need to do before creating the optimizer is create the first simplex. With the Rosenbrock cost function in 2D, this means setting up 3 vectors of 2 parameters:

  startPoint = np.empty((3, 2), np.float)
  startPoint[:,0] = 1.
  startPoint[:,1] = 0.
  startPoint[1,0] -= .1
  startPoint[2,1] += .1

Then, the optimizer instance can be built and called:

  optimi = optimizer.PolytopeOptimizer(function = Rosenbrock(), criterion = criterion.criterion(ftol = .0001, iterations_max=50), x0 = startPoint)
  print optimi.optimize()

As usual, you can use a recorder to display how the optimizer behaved, which is the purpose of next section.

Interactive results check

In the tutorial folder on github, I’ve put an interactive application with Traits and Chaco. I would have prefered that the optimization took place in a background thread and that the main GUI updates itself with new iteration states, but I couldn’t get it to work (I need to practice Traits a bit more). Still, you can explore the results with this app.

The other way to check the results is to save each state, display it with matplotlib and make it a movie:

As usual with the Rosenbrock function, the optimizer finds quickly the minimum valley but then takes a lot of iterations to find the actual solution. As we need to move the simplex in this bent valley, the number of iterations increases rapidly. Still, the optimal value is found in the end.

Conclusion

The Simplex algorithm is a very crude one but it can be effective when the cost function is too complex to derive or if the gradient is too costly to compute. It is still very sensitive to the stopping criterion as well as the starting polytope, but it can be a first way to check your result before trying something different like a gradient-based optimization, genetic algorithms…

Next tutorial will be about Quasi-Newton methods, and then the following one may be on optimization comparisons.

Changelog:

• Fixed Polytope/Simplex/Nelder-Mead
• Fixed the Quadratic Hessian helper class

Additional tutorials will be available in the next weeks.

Old announces:

• 0.1

Implementation

The implementation is not optimized, the only noticeable thing is that I don’t create subnodes for all points, I only split nodes when they have more than n points, and then according to the best dimension (best is defined by the biggest distance on one axis in the node, the root distance being set by hand). The search is of course really close to the cover tree implementation.

Also, when doing the search, I can override the original distance function with a new one. This is because when doing the search, some dimensions may be of greater interest than others. At least, this is the case I have. When creating the tree, I don’t have access to the actual distance, and it could even change during the process.

Results

OK, some comparisons now. I’ve changed the number of nodes to 1 million. I’ve relaunched the cover tree test, and this is the result:

Build time 12:20:49.484375
Out time (linear) 00:00:04.500000
Out time (cover_tree) 00:00:00.062500

The first figure is true, it is more than 12 hours. You have to search a lot to be profitable to use a cover tree.

The kd-tree:

Build time 00:00:01.890625
Out time (linear) 00:00:02.921875
Out time (kd-tree) 00:00:00.875000

Strangely the timing for the linear search is a bit different, but I guess this is due to the long time the program had to run. The timing is not as good, but it is only one test. I would need to test far more points, but due to the high cover tree build time, I didn’t try yet.

Now, if we use another distance (infinite instead of L2):

Out time (linear) 00:00:04.718750
Out time (kd-tree) 00:00:00.890625

The linear time is a bit higher, like the one from the cover tree. What is interesting is that the kdtree time is comparable to the one from the L2 search time. And this is what I wanted to check.

Conclusion

Kd-trees are not something earth-shaking, but sometimes it is still very relevant. Of course, I didn’t optimize any of the cover tree or kd tree implementation, so there is place for improvement. What I really like for the kdtree is that if you change the distance function, it can still work, whereas a cover tree, being implemented on a specific distance, cannot easily.

Kd-tree code on Github is here.

I could have used the Numerical Recipes version, but the licence cost is not worth it, and the code is not generic enough, not explained enough. And also there are some design decisions that are questionable (one method = one responsibility).

Design

The main choice was to use Eigen for the inner computation in the code. This library is perfectly fitted for numerical computation in C++. The parameter set for the function could be any matrix class, but Eigen’s ones are ensured to be compatible inside the code. The simplex is stored as an Eigen dynamic array to cover any dimension possible cases.

There are several tools that are pretty much generic in a simplex algorithm: the state (current parameters, iteration…) and the stopping criterion (a conjunction of a maximum iteration and a relative value difference). So these concepts will be implemented as separate generic classes.

The main class, Simplex, will encompass the actual simplex algorithm. There are three inner central methods: the one that creates the first state (initialize_polytope()), the one that will detect the best parameter set, the worst one and the worst second (find_best_worst_near_worst()) and the one that create a new candidate based on the current simplex and the selected worst parameter set (create_new_parameters()).

The main method is optimize(), building on the class attributes that were set before and the main methods that were exposed before. The algorithm follows strictly the implementation on the Wikipedia page. The sizes of expansion, contraction could have been made template, or even modified at compile time, it is an exercise left for those that want a fully generic class.

A static builder was written to ease the creation of a new instance. It doesn’t do anything else than instantiating the class with the constructor, its only goal is to remove the need for knowing the exact template parameter and to fully specify them. This can be done because of the auto keyword in C++11.

Using it

Let’s start with a function to optimize. I’ve chosen the Rosenbrock function as it is not complex to compute but still difficult to optimize. For the parameter kind, I’ve also used Eigen.

typedef Eigen::Vector2f ParameterTypeRosenbrock;
 
// Rosenbrock cost function
struct Rosenbrock
{
  typedef float DataType;
  typedef ParameterTypeRosenbrock ParameterType;
  float operator()(const ParameterTypeRosenbrock& parameters) const
  {
    return (parameters(1, 0) - parameters(0, 0) * parameters(0, 0)) * (parameters(1, 0) - parameters(0, 0) * parameters(0, 0))
        + (1 - parameters(0, 0)) * (1 - parameters(0, 0));
  }
};

Now, I can create the optimizer with all the different parameters.

int main(int argc, char** argv)
{
  ParameterTypeRosenbrock start_point;
  start_point << 10, 10;
 
  Rosenbrock fun;
 
  long max_iterations = 1000;
  float ftol = 0.00001;
 
  auto optimizer = Optimization::Local::build_simplex( // Builder to generate the correct optimizer
      fun,                                             // Used to infer function type
      Optimization::Local::make_and_criteria(Optimization::Local::IterationCriterion(max_iterations),
          Optimization::Local::RelativeValueCriterion<float>(ftol))); // Composite stoping criterion
 
  optimizer.set_start_point(start_point);  // Starting parameters
  optimizer.set_delta(1);                  // Simplex size
  optimizer.optimize(fun);                 // Optimization start
 
  std::cout << "Starting point: " << start_point << std::endl;
  std::cout << "Starting value: " << fun(start_point) << std::endl;
  std::cout << "Best point: " << optimizer.get_best_parameters() << std::endl; // Retrieve the best parameters
  std::cout << "Best value: " << optimizer.get_best_value() << std::endl;      // Retrieve the best value
}

With C++11, it is easy to find the type of the optimizer, thanks to the builder. It is still cumbersome as the builder function needs the function instance to infer the inner data type as well as the parameter set type.

Starting point: 10
10
Starting value: 8181
Best point: 1
1
Best value: 5.68434e-014

The simplex does not need such a complex generix framework. The implementation usually uses the criterion, so I didn’t have to implement it the way I did. Still, it is interesting for a more complex optimizer like a quasi-Newton or conjugate-gradient.

Conclusion

This code is not world shaking but it explains how I consider optimization. The other optimizers from scikits.optimization could be written in the same way and they could be better understood by someone versed in object oriented languages than usual procedural code.

As usual, it can be found on github: https://github.com/mbrucher/CppOptimization

Content and opinions

The first three chapters help you set up your environment with your own sequencer. The book tries to be sequencer agnostic, although the author mainly resorts to Sonar or Reaper. So the first chapter defines some terms and install some free plugins that will be used through the book. The second one is about what sound is while the third explains how effects can interact with sound to create something.

Afterwards, the author deals with actual effects. The next three chapters are dedicated to frequency and dynamics effects. The different types of equalization processes and each instrument has specific frequency bands that need to be shaped. As the book comes with several EQ plugins, there is room for tests with the different types of equalization and different instruments. Then the first dynamic effect is noise gate. I think it is one of the effects you have to know how to use fast as it is mandatory for real recorded tracks. There are different levels in the noise gate (it can have side chains, more filters), but the simplest one is easy to use. The tutorial does not handle all the different combinations, only the most usual ones. The compressor chapter is a bit more complex, as compression can be used in different ways. Not only ways has a complete tutorial, but the reader is exhorted to test all of them. There are several “checkpoints” where another audio project is used than the ones in the tutorials. In this case, you don’t have a reference where you will apply the tutorial without thinking, you have to test all parameters and different ways of adding a plugin.

The next three chapters are about time effects. I think time effects are what makes the most innovation in tracks. The first plugins are always on all audio tracks we listen to (with the addition of reverberation, which is the object of the second chapter in this part), but the other are used with scarcity and it makes the plus that surprises the listener (in my opinion). So the first chapter is about delay and chorus. A delay can be very complex (from a simple tap to multi-tap with a spatial aspect) whereas the chorus is something more simple, but nonetheless complex to use. They can add a spatial dimension to the audio track, and as such they must be tuned with care. Reverberation is the subject of the second chapter, and it is one of the mandatory effects (as I’ve said before). The author spends several pages on synthetic reverberation, but I prefer convolution. The reason is that you can exactly specify reverberation you want and need, and the CPU cost is not a problem anymore. Convolution is the future, why spend so many pages on synthetic? The last chapter is about other time effects, like phaser, flanger, tremolo… They will be applied to different tracks depending of what you want to achieve. It’s perhaps the most complicated chapter because you need to remember what all these effects sound like and imagine what they sound like on your work.

Almost all the remaining chapters are dedicated to making all the pieces of the puzzle fit together. The tenth chapter is about effects that used for something different from their usual goal (like deessers) and plugin effects that chain several effects (possibly in the order you want). They are incredibly valuable, although they are sequencer specific or not free (so no tutorials, you have to try by yourself). One of the great benefit of a sequencer over traditional hardware is automation, the subject of the next chapter. The author warns of the different aspects of automation and what you can and cannot do as well as why. A small chapter is dedicated to stereo effects. It is quite short and as the other dimensions are tackled before (delay and reverb), it is strange to have the main dimension’s discussion now. The next two chapters help setting the threads between all the plugins. There is nothing actually new, but the tutorials help testing different serial and parallel combinations. It’s what makes a final mix, and it is properly explained.

The final chapter before the conclusion is about mastering effects. The author provides a guide to check your mix and prepare it for the mastering phase (on your own or by someone else), and then he reviews the effects you will need. There is no discussion on the loudness war, it could have been a good talk.

The conclusion gives you some hints to sort your plugins and choose the proper one, it is something one can easily forget.

Conclusion

With all the tracks and plugins on the CD, the explanations of each effect and then the tutorials (for each effect, even if I didn’t say it), this book is what I expected from a book on audio effects. If you want to know what audio effects are, go and buy the book, you won’t regret it.

Implementation

Contrary to kdtrees, the cover tree nodes always represent one point. So you have as many nodes as you have data points. Each of them can have several children, sorted by “level”, or distance hierarchy. Each level halves the distance between the node and its children. When constructing the tree instance, the distance you use is given as parameters (and its type is a template parameter, as well as the point type and the data type used for comparisons).

Construction of the tree is done by searching in the tree for the closest node and the lowest hierarchy possible. In all time, the min and max levels are tracked in the tree. On insertion, the max level is tried, and if it fails, the max level is increased by one. The failure occurs if the new point is furthest of the root point that the maximum distance on the max level. At each level, the set of the nodes that are closer to the new point than the current maximum distance are kept (function populate_set_from_node). If the least distance is greater than the next level maximum level (function find_min_dist, the recursion stops and the new node is inserted in the first node in the set.

When doing a search, a list of distances and points (of type NearestNodesStructure) is created and its size is kept at the number of neighbors k. On entrance of a new level (level_traversal), all nodes that have children that could be close enough (current distance less than maximum level distance + distance to the current k-th neighbors found with find_k_dist) are added to the list. In order to speed search, a map is not used, but a partial sort is done at each iteration (which makes find_k_dist work without sorting the container again). Without this, the performance is worse than a linear knn search.

I didn’t implement other functions like removal or traversal. I don’t need them, but feel free to add them if you need. Sidenote: I will try to find time to refactor the code, because it is not very clean at the moment (especially the knn method).

Results

So some figures. I’ve populated a tree with 100k elements, and then made a 10-neighbors search with the cover tree, and with a standard linear search. Without the construction time, the search is about ten times faster than a linear one. Of course, it all depends on the cover tree balance. You could design a cover tree that has too many levels (but also the data would be unusual) and that would result in an almost linear time.

I’ve done a profile with callgrind, but the majority of the time is in the tree construction. No point in checking the search and displaying it here as the following timings show:

Build time 00:09:38.530397
Out time (linear) 00:00:00.238989
Out time (cover_tree) 00:00:00.012518

What must be optimized is clearly the build time.

The test code is in the main.cpp file on github.

Cover Tree code on Github is here and the publication with detailed algorithms here.