Be aware of algorithms that can be useful to you. There are many textbooks on algorithms and data structure. One of the more encyclopedic is the book by Cormen, Leiserson, Rivest, and Stein [23].

Scientific and engineering software almost always needs to be efficient in both time and memory. You should first consider optimizing at a high level, choosing data structures and algorithms that are inherently efficient. At a lower level, you should understand how computers work, and how to efficiently use modern computer architectures.

Choosing good algorithms and data structures is the foundation of designing efficient software. Without this, no amount of low-level tweaking of the code will make it run fast. This is especially true for large problems. By tuning your code you could gain an improvement of anything from, say, a factor of two to a factor of ten, in the speed of an algorithm. But moving from an algorithm that is O(n3) in time and O(n2) in memory to one that is O(n2) in time and O(n) in memory can give you a much greater benefit, especially if n ≍ 10 000 or larger. Sometimes you can get approximate algorithms that are fast, but the speed will depend on how accurately you want the answer. Here you need to know the problem well in order to see how large an error can be tolerated.

Numerical algorithms

Since the development of electronic digital computers (and even well before then) there has been a great deal of development of numerical algorithms. Many are highly efficient, accurate, and robust.

Dynamic memory allocation

In the middle of a routine you realize that you need some scratch space n floating point numbers long, and n is a complicated function of the inputs to the routine. What do you do? Do you add another scratch space argument to your routine along with its length as an argument and check that it is big enough (stopping the program immediately if it is not)? Often a better idea is to allocate the memory needed. In Fortran 90 this is done using the allocate command; in C you use the malloc function; in C++ and Java the new operator will do this task; in Pascal the allocation is done by the new operator, but the syntax is different from C++ or Java. These can be used to dynamically allocate memory. All of these commands return, or set, a pointer to the allocated block of memory.

The allocated block of memory is taken from a global list of available memory which is ultimately controlled by the operating system. This block of memory remains allocated until it is no longer accessible (if garbage collection is used), explicitly de-allocated, or the program terminates. So dynamically allocated memory can be used to hold return values or returned data structures. Dynamically allocated memory can be passed to other routines, and treated like memory that has been statically allocated, or allocated on a stack in most respects.

The data structure that controls the allocation and de-allocation of memory is called a memory heap. A memory heap may contain a pair of linked lists of pointers to blocks of memory.

Picking algorithms vs. keyhole optimization

Optimizing the performance of software requires working with at least two different points of view. One is the global view where the questions are about how to structure the overall architecture of the system. Another view is how the individual routines are written, and even how each line of code is written. All of these views are important. Selecting good algorithms is as important as selecting good hardware and implementing algorithms efficiently.

When the term “optimization” is used in computing, it is often taken to mean something like picking compiler options or “writing tight code” that uses the least number of clock cycles or operations. This is clearly important in writing fast, effective software, but it is only a part of the process, which begins early in the design stage.

Usually the first part of the design stage is the design of the central data structures and databases that will be used. These should be chosen so that there are well-known efficient algorithms for handling these data structures, preferably with readily available implementations that perform correctly, reliably and efficiently. Then the algorithms to carry out the main tasks need to be selected. An important guide to selecting them is their asymptotic complexity or estimate of the time needed. However, this is not the only guide; see the last section of this chapter for more information about how to refine this information and to make sure that the asymptotic estimates are relevant in practice. The last part is the detailed design process where the algorithms are implemented, and this should be done with an eye on efficiency to ensure that the whole system works well.

The development of the Internet has contributed to the development of public libraries of scientific software. Some of this development has occurred through the efforts of individuals, through Internet collaborations (as with the development of Linux), and through government supported software development by academics and others (as with the development of LAPACK). There is a wide range of other software packages for scientific computing, many now written in C/C++, although Fortran and other languages are used for parts of many of these systems: PETSc (which supports the use of MPI for parallel computation), IML++ (an iterative methods library in C++), SparseLib++ (for handling sparse matrices in C++), and PLTMG (for solving partial differential equations).

In parallel with this, there has also been a tremendous development of commercial numerical software. Beginning in 1970 the Numerical Algorithms Group (NAG), based in the UK, developed libraries which have been sold commercially as the NAG libraries since 1976; the Harwell library was also developed in the UK; the IMSL libraries were developed commercially in the US. Another set of numerical libraries, called SLATEC, was developed by the Sandia and Los Alamos US National Laboratories and the US Air Force. These are available through netlib (see the next section). Perhaps the most spectacular example of commercial numerical software is the development of MATLAB. Initially a collection of Fortran 77 routines based on the early LINPACK and EISPACK libraries for dense matrix computations with a text interface, MATLAB has evolved into a full-featured interactive programming language with special support for numerical computation and scientific visualization.

Introduction

Any picture may be conceived as a mathematical object, lying on part of the euclidean plane, each point having its own colour. Then it is a strange and wonderful entity. It is mysterious, for you probably cannot see it. And worse, you cannot even describe it in the type of language with which you normally talk about objects you can see; at least, not without making a lot of assumptions. But we want to be able to see, to describe and to make pictures on paper of fractals and other mathematical objects that we feel ought to be capable of representation as pictures. We want to make mathematical models for real-world images, biological entities such as leaves and many other types of data. To be able to do this we need certain parts of the language of mathematics, related to set theory, metric spaces and topology.

Code space There is a remarkable set, called a code space, which consists of an uncountable infinity of points and which can be embedded in the tiniest real interval. A code space can be reorganized in an endless variety of amazing geometrical, topological, ways, to form sets that look like leaves, ferns, cells, flowers and so on. For this reason we think of a code space as being somehow protoplasmic, plastic, impressionable and capable of diverse re-expressions, like the meristem of a plant; see Figure 1.1. This idea is a theme of this chapter and of the whole book.

In the BSA Chapter 3 we learned that a DP algorithm for pairwise sequence alignment allows a probabilistic interpretation. Indeed, the equivalent equations appear in the logarithmic form of the Viterbi algorithm for the hidden Markov model of a gapped sequence alignment. The hidden states of such a model, called a pair HMM, correspond to the alignment match, the x-gap, and the y-gap positions. The pair HMM state diagram is topologically similar to the diagram of the finite state machine (Durbin et al. (1998), Fig. 4.1), although the pair HMM parameters have clear probabilistic meanings. The optimal finite state machine alignment found by standard DP is equivalent to the most probable path through the pair HMM determined by the Viterbi algorithm. Both global and local optimal DP alignment algorithms have Viterbi counterparts for suitably defined HMMs. Interestingly, the HMM has an advantage over the finite state machine because the HMM can compute the full probability that sequences X and Y could be generated by a given pair HMM; thus, a probabilistic measure can be introduced to help establish evolutionary relationships. This full probabilistic model also defines (i) the posterior distribution over all possible alignments given sequences X and Y and (ii) the posterior probability that a particular symbol x of sequence X is aligned to a given symbol y of sequence Y. However, real biological sequences cannot be considered to be exact realizations of probabilistic models. This explains the difficulties met by the HMM based alignment methods for the similarity search (Durbin et al. (1998), Sect. 4.5), while more simplistic finite state machine methods perform sufficiently well.