Design Document for DirectSearch and its Derived Classes

P.L. (Anne) Shepherd

Introduction

Direct searches [13] are used to minimize a function where . The objective function is assumed to be continuously differentiable, but we either do not have or cannot rely upon the gradient.

Two categories of direct searches are pattern searches [4,8,32] and simplex searches [7,18,25]. A pattern search looks for a minimum among the points in a design, each of which is on a node of a rational lattice or grid [32]. The rational lattice, where the distance between nodes is a rational factor of the initial step length (conventionally called ), is a sine qua non of a pattern search, necessary to the proof of pattern searches' global convergence [32]. A simplex search uses a nondegenerate simplex, at whose vertices the function is sampled [7]. There is no theoretical guarantee of convergence to a stationery point for a simplex search that is not also a pattern search; neither of the two most widely used simplex searches, those of Nelder and Mead [18] and of Spendley, Hext, and Himsworth [25], are pattern search methods.

I will assume in writing this document that the reader has some basic familiarity with the mechanical details of these searches; interested readers may wish to consult [13], which gives a historical overview of direct search optimization methods, and for further detail, [2,4,7,8,12,14,18,25,32].

The DirectSearch class was devised as an umbrella class for two separate sets of searches--five pattern searches and three simplex searches--that had been written by two Computer Science Honors students [4,7]. My work here began as an attempt to bring these eight searches into harmony with each other in a manner consistent with good programming practices.

Of primary importance in this project was reworking the interface to make it easier (more intuitive and therefore less error-prone) to use, more uniform among the various searches, less reliant on outside files, and easier to maintain. I also wanted to preserve, as far as possible, the spirit of the software as it was originally written, and of course to maintain the algorithmic integrity of the searches themselves. At the same time, I wished to put the user in a position where testing the searches, comparing them against each other, and switching back and forth among them, could be easily done with meaningful results.

It should be noted that this document is concerned mainly with the decisions and compromises that were made in modifying and augmenting the code I was given. More traditional class documentation, generated by Doc++ [34], is available at [28]. A brief user's manual for the search is available as well, at [27].

Background

In 1999-2000, Dolan [4] and Gurson [7] implemented a number of direct searches in order to study their behavior, test their effectiveness, and, if possible, improve the algorithms. Dolan concentrated on pattern searches and wrote C++ code for five of them:

1. CoordinateSearch, a very simple search along the coordinate directions [21];
2. CompassSearch, a form of coordinate search using a somewhat greedier approach [4];
3. NLessSearch, a minimal-positive-basis search [14];
4. HJSearch, a search based on the method described by Hooke and Jeeves [8]; and
5. EdHJSearch, Dolan's edited version of the Hook and Jeeves search [4].
Gurson studied simplex searches, only one of which (SMDSearch) can also be classified as a pattern search. He implemented three of these, also in C++:
1. SHHSearch, based on the simplex method of Spendley, Hext, and Himsworth [25];
2. NMSearch, based on the Nelder-Mead [18] simplex method, as specified by Lagarias et al. [11]; and
3. SMDSearch, a sequential version of Multidirectional Search [31].

In implementing their searches, both Dolan and Gurson used container classes defined in the files cmat.h and vec.h in implementing their searches; cmat.h holds a Matrix template class, and vec.h contains a Vector template class. These two container classes were written for numerical linear algebra, and suited our purposes well; These files originally were distributed as part of the Template Numerical Toolkit [29]. The versions we use include modifications by Siefert [24].

Goals

Each of Dolan's five searches was implemented as a concrete class derived from the PatternSearch abstract base class. Gurson's three searches, on the other hand, were in three separate classes, with no abstract base class. Although Gurson used Dolan's code as a starting point, in the end there were numerous differences, both in the implementation and in the interface, between the two.

My goals in designing the DirectSearch class were:

1. Most importantly, to bring Dolan's PatternSearch and derived classes, and Gurson's three simplex search classes, under a single umbrella class, with:

• as uniform a public interface as possible;
• a more self-contained form, with less reliance on C-style macros and #define statements;
• increased reliance on code reuse, with a view toward easier maintenance and improved readability; and
• preservation, whenever reasonable, of the original spirit of the software, with no material changes at this time in the conduct of the actual searches.

2. To provide for the searches:
• more extensive internal documentation,
• an external design document, and
• a user's manual.

3. To increase portability, taking into account the dynamic nature of the C++ language at the time of this writing.
4. to get the public interface into a form that would be simple to use and easy to maintain, so that any later implementation changes would be transparent to users.
5. To work with the two major users of the PatternSearch code, Chris Siefert and Tony Padula, who had changed Dolan's software to work with their MAPS [22] and krigifier [19] software. Here, the goal was to standardize DirectSearch so there would be only one version in use by everyone, in order to minimize future maintenance.

Scope of DirectSearch and its derived classes

Scope of the DirectSearch class

The DirectSearch class was conceived as an abstract base class, a convenient umbrella class to bring the pattern and simplex searches together. It allows for more code reuse and thus eases maintenance; however, it does no searching. In formulating DirectSearch, we decided that a DirectSearch object should know, at a minimum, the following things, which would be represented as fields (or data members) of the class:
1. the dimension of the problem,
• implemented as long dimension;
2. the point with the minimum function value found thus far,
• implemented as Vector<double> * minPoint;
3. the value of the function at that point,
• implemented as double minValue;
4. the maximum number of function calls allowed,
• implemented as long maxCalls;
5. the number of function calls/evaluations made thus far,
• implemented as long functionCalls;
6. the stopping step length, i.e. the smallest refinement of the underlying search design that will be allowed before the search is stopped,
• implemented as double stoppingStepLength;
7. where to find its design structure,
• implemented as Matrix<double> * design, a pointer to a Matrix object.
8. how to count the number of function values for the purpose of deciding when to terminate the search,
• implemented using the field exact_count.

Note that by refinement of the search design'' I refer to the fact that many of the searches in the DirectSearch class have the equivalent of a shrink step somewhere in the algorithm. In a pattern search, if no new minimum is found after the pattern (or some reasonable subset of the pattern) has been searched, the pattern itself shrinks; i.e. the underlying rational lattice is made finer. In the simplex searches, it is the simplex that does the shrinking. Simply stated, when the lattice becomes finer or the simplex edge lengths shorter than a given tolerance [4,5,7,12], the search is halted.

In terms of the methods (or member functions) for the class, we thought that a DirectSearch object should be able to do, at a minimum, the following:

1. construct and destroy itself,
• implemented using three constructors and a destructor;
2. initialize its design structure in various ways,
• implemented using
void InitializeDesign(const Matrix<double>*designPtr),
virtual void CleanSlate(long dim, Vector<double>& startPoint),
and the pure virtual method
virtual void ReadInFile(istream & fp) = 0,
as well as with constructors and mutators.
3. allow the user access to the above fields as appropriate,
• implemented in various accessors;
4. print out useful information about the search,
• implemented as
virtual void PrintDesign() const = 0
--note that this method is pure virtual in DirectSearch;
5. execute the search itself,
• implemented as the pure virtual methods
virtual void BeginSearch() = 0 (public), and
virtual void ExploratoryMoves() = 0 (protected); and
6. decide when to stop the search, whether based on design refinement or on the maximum number of function calls sanctioned by the user.
• implemented using a number of methods (see §5.3.

Scope of the PatternSearch class

Dolan formulated PatternSearch as an abstract base class for the five pattern searches she implemented. I modified the PatternSearch class and its child classes from Dolan's originals in order to fit them into the new DirectSearch hierarchy.

What a PatternSearch object should know, in addition to what a DirectSearch object already knows, is:

1. the density/refinement of the underlying lattice (i.e. ),
• implemented as double delta;
2. the number of trial points in the search pattern,
• implemented as long patternLength; and
3. the initial step length of the pattern (i.e. the initial setting for delta),
• implemented as double initialStepLength;

What a PatternSearch object should be able to do, in addition to the methods inherited from DirectSearch, is:

1. find the next point in the pattern,
• implemented by the protected method
virtual void NextPoint(long index,
const Vector<double> &currentPoint,
Vector<double> & nextPoint);

2. scale the pattern when a reduction step is called for,
• implemented by the protected method
virtual void ScalePattern(double scalar); and
3. allow the user access to the above fields as appropriate,
• implemented by various accessors and mutators.

Scope of the SimplexSearch class

I devised the SimplexSearch class as an abstract base class for the three simplex searches originally implemented by Gurson [7] as separate, stand-alone classes. SimplexSearch inherits from DirectSearch but not from PatternSearch.

We decided that at a minimum a SimplexSearch object should know:

1. the centroid through which to reflect the simplex (where applicable),
• implemented as Vectordouble *centroid;
2. , the scalar for the shrink step,
• implemented as double sigma;
3. the coordinates of the vertices of the simplex at any given time in the search,
• implemented as double *simplexValues;
4. the index (in the design/simplex matrix) of the vertex that generated the current minimum function value,
• implemented as long minIndex;
5. the index of the vertex to be replaced,
• implemented as long replacementIndex;
6. the type of simplex to be used to start the search,
• implemented using bool SimplexSpecified and two public methods (see §5.1.2);
7. the lengths of the edges in the simplex used to begin the search,
• implemented as Vectordouble *starting_edgeLengths;
8. when to stop the search, based on the size of the simplex,
• implemented using int toleranceHit and bool Stop_on_std.

What a SimplexSearch object should be able to do, in addition to the methods inherited from DirectSearch, is:

1. allow the user to specify a type of starting simplex (although we provide a default),
• implemented using the bool Stop_on_std field and the pure virtual methods:

• virtual void ChooseRightSimplex() = 0 and
• virtual void ChooseRegularSimplex()= 0;

2. allow the user to access the above fields as appropriate,
• implemented in various accessors;
3. initialize different types of simplices,
• implemented in the abstract base class SimplexSearch with three protected methods:

• virtual void Initialize_Right(Matrixdouble *plex)
• virtual void Initialize_Regular(Matrixdouble *plex)
• virtual void InitGeneralSimplex(const Matrixdouble *plex);

4. shrink a simplex as necessary,
• implemented with the protected method virtual void ShrinkSimplex();
5. replace points in a simplex as necessary,
• implemented with the protected method
virtual void ReplaceSimplexPoint(long index,
const Vector<double> &newPoint);

6. find a centroid of a simplex as needed.
• implemented by the protected method virtual void FindCentroid();

The process of choosing and initializing a simplex will be discussed in §5.

The Interface

As stated in the introduction, my primary concerns in the refactoring of the interface were ease of use, uniformity across the classes, ease of maintenance, less reliance on outside files, and preservation of the spirit and algorithmic integrity of the original software.

With these goals in mind, I made a number of decisions that changed both the public and non-public interfaces. The most obvious change was the implementation of two new abstract base classes, DirectSearch and SimplexSearch. The integration of the two groups of searches naturally led to name changes for a number of fields and methods; for example, the variables field in the pattern searches [4] and the dimensions field in the simplex searches [7] became dimension in DirectSearch. Other considerations, which I will discuss in this section, were:

1. construction and initialization;
2. connecting the user's objective function to our searches;
3. when and how to halt each of the searches;
4. protection of and access to the fields; and
5. simplifying the interface.

Throughout the process I worked with our two known users to ensure that they would be able to make these changes with minimal disruption to their code.

Construction and initialization

In the original software [4,7], the simplest constructor in each class took only one argument, the dimension of the optimization problem. The starting point for the search was introduced in the objective file in the pattern searches, and in the simplex initialization process in the simplex searches. I made the starting point a protected field of the search, and made it a parameter to all constructors except the copy constructor. The user is going to have to provide a starting point one way or another; it is easiest to do so right at the start, in a constructor.

I decided to make the starting point a member of the Vector class (Dolan had it as an array, Gurson as a Vector) because we had to turn it into a Vector anyway. We were already too committed to using the Vector class to change it without making substantial modifications to the public interface, so there was nothing to be gained by hiding it.

Constructor considerations in the simplex searches

With the simplex searches, there are a number of fields that a user might want to set at the beginning of a search; and where the pattern searches have three constructors each, the simplex searches have five. Here, I have tried to provide flexibility for the user without excessive constructor proliferation.

In order to conduct a simplex search it is necessary to know the edge lengths of the simplex that is used to start the search. In Gurson's original searches, these were sent in to the simplex initialization methods (see below). I made starting_edgelengths a protected field of the SimplexSearch class and made it an argument to every constructor except the first constructor (which requires only the dimension and the starting point), and the copy constructor. I left starting_edgelengths out of the first constructor's argument list in order to allow at least one means of construction that was exactly the same for all searches. This constructor sets all edge lengths to a default value (currently 2.0, which is the same as the default step length in the pattern searches). The user can also set the edge lengths by calling the public method SetStartingEdgeLengths(const Vectordouble&).

For more details on how starting_edgelengths is used by SimplexSearch, see the next subsection.

Simplex initialization

With Gurson's original code [7], users had to call initialization methods explicitly before they could begin the actual search. Although this worked perfectly well, I wanted to make the user's steps in executing a basic simplex search as much as possible the same as would be used to execute a basic pattern search.

The user can choose either a regular or a right simplex; the right simplex may be designated to have either fixed-length or variable-length edges. (This refers to the leg'' edges. The lengths of the hypoteneuse'' edges of a right simplex will of course depend on these, and are not specified by the user.) There is currently no provision for the user to specify an arbitrary simplex that is neither right nor regular.

The original methods, in all three searches, were:

   void InitRegularTriangularSimplex
(const Vector<double> *basePoint,
const double edgeLength)
void InitFixedLengthRightSimplex
(const Vector<double> *basePoint,
const double edgeLength)
void InitVariableLengthRightSimplex
(const Vector<double> *basePoint,
const double edgeLength)
void InitGeneralSimplex(const Matrix<double> *plex)


The first three were used to choose and begin initialization of a specific kind of simplex: a regular simplex, a right-angled simplex with edges of fixed (user-specified or default) length; or a right-angled simplex with edges of variable, user-specified length. The utility method InitGeneralSimplex was called as part of the initialization of every new simplex. All four were public methods.

In the new SimplexSearch class and its derived classes, there are two public mask methods:

   void ChooseRegularSimplex()
void ChooseRightSimplex()


and three protected methods. In the abstract base class SimplexSearch, these are:

  void Initialize_Right(Matrix<double> *plex)
void Initialize_Regular(Matrix<double> *plex)
void InitGeneralSimplex(const Matrix<double> *plex)


and in the three concrete search classes they are:

  void InitRegSimplex()
void InitRightSimplex()
void InitGeneralSimplex(const Matrix<double> *plex)


Because starting_edgelengths is now a data member of the class, and is a pointer to a Vector object, it is no longer necessary to differentiate between a fixed-length and a variable-length right simplex. If the entries in *(starting_edgelengths) are all equal, then the simplex will be fixed-length. If a regular simplex is chosen, only the first entry in the Vector is used. I have also pushed the memory allocation and destruction for the Matrix pointer plex down into the protected part of the implementation.

In the original simplex searches the user was required to initialize a simplex up front, before the search could begin. Now if there is none specified, a default right-angled simplex will be initialized, using either a user-defined starting_edgelengths Vector or a default Vector whose entries are all equal, set to a default value (currently 2.0). The bool field SimplexSpecified has been added to manage this. It is set to false by all constructors except the copy constructor. If the user calls ChooseRegularSimplex() a regular simplex is initialized using InitRegSimplex(), and SimplexSpecified is set to true; otherwise, when the user calls BeginSearch() to start the actual search, BeginSearch() calls InitRightSimplex() using whatever values (user-defined or default) are stored in *(starting_edgelengths). The call to ChooseRightSimplex(), similarly, initializes a right simplex. Because this is currently the default setting, it is not necessary for the user to make the explicit call if a right simplex is preferred. In cases where a copy constructor has been used to replicate a search and the user wants a right simplex for one version and a regular simplex for the other, however, this is how it can be done.

Connecting the user's objective function to our searches

There is unfortunately no way of getting around the fact that the user will have to take steps to allow our software to work with the function to be minimized. We have tried to make this as easy as possible. There are now three different ways of connecting a user's objective function to the search programs:

1. Using Dolan's original approach (but with new versions of objective.h and objective.cc), where the user writes a wrapper for the objective function and places it in the objective file in place of our toy function.

2. Using the Evalf() function written by Siefert, and used by Gurson, to call the objective function using system calls.

3. Sending the function name to the search as a parameter to the constructor. That way the whole thing can be done in the main program and the user doesn't even need to #include objective.h'' or compile with objective.cc (but will of course need to compile and link with the file containing the user's objective function). This is discussed below.

More specific information on integrating the user's function into our searches can be found the user's manual for DirectSearch [27].

Using the constructor to specify the objective function

To allow DirectSearch to work properly with Siefert's MAPS program [22], I added to all searches a special constructor using void pointers; for example:

  CompassSearch (long dim,
Vector<double> &startPoint,
double startStep,
double stopStep,
void (*objective)
(long vars, Vector<double> &x,
double &func, bool &flag, void* an_obj),
void * input_obj);


This allows the user's main program to specify the name of the objective function by making it a parameter to the constructor. The user's objective function may be:

• inserted into our objective.h and objective.cc files by name,
• in its own file, which must then be linked with an #include statement and compiled with the search files, or
• a part of a separate class, as in MAPS [22]. More details can be found in the user's manual [27].

Typical usage of this constructor (here, in a CompassSearch) is as follows:


...
long n = 5;
Vector<double> minVec(n, 2); // this will initialize a Vector of 2's
double start = .25           // starting step length
long tol = 10e-8;            // default stopping step length
...

// Note that we typically don't need the constructor's last argument, so
// we send in a NULL for the parameter void* input_obj.
// MyFunction is the user's function to be minimized.

CompassSearch CMPsearcher(n, minVec, start,
tol, MyFunction, NULL);
...


The last argument to this constructor (input_obj, which sets the field some_object, but which is set to NULL in the above example) may instead be used to specify a class object, where the function to be minimized is a method of that class.

In that case, the call to the constructor might look more like this:

  ...

// construct an instance of the class that holds the objective
// function.
MyClass* MyClassObject  = new MyClass();

CompassSearch CMPSearcher(n, minVec, start, tol,
MyFriendFunction,
(void *)MyClassObject);


In this example, the field fcn_name points to MyFriendFunction, a friend'' function of the class MyClass. MyFriendFunction wraps the actual objective function, a member function of MyClass.

I have provided a toy class trigclass and an example program main_obj.cc to illustrate this process. Listings of these two programs are included as Appendix B to this document. The class documentation [28] and user's manual [27] give further details.

This constructor was based on one that Siefert had added to the version of CompassSearch he had modified for use with MAPS [22,23]. The parameters fcn_name and some_object were made protected fields of DirectSearch class. In addition to solving a MAPS compatibility problem, this constructor allows the objective function to be specified by the user's main program without using a separate objective.cc file as a wrapper.''

The objective function still must be either formulated or wrapped so the signature matches that of our fcn_name field. The last field (which sets the field some_object) can be used for additional information as required by the user. Although void pointers are not considered particularly safe programming devices, I decided to add this constructor to allow compatibility with MAPS and to give other users flexibility that could not be otherwise provided.

When and how to halt the searches

The traditional termination criterion of a direct search is the step size: the search is halted when the step falls below some appropriately small tolerance. [5]. By default, all the searches in the DirectSearch class will terminate when their step size (or simplex edge length) falls below approximately the square root of machine epsilon.

The nature of a pattern search is such that a reduction in occurs only when a search of the points in the pattern yields no improvement [4,5,12]. Thus, as Dolan et al. point out [5], decrease in provides a natural measure of progress which can reliably be used to test for convergence.

Function Call budget in pattern searches

Given that we have a good measure for assessing the asymptotic behavior of pattern search, there is, nevertheless, a more pragmatic measure to consider--one that applies to the simplex searches as well. Because the function evaluations in the optimization process are typically expensive, it is often desirable to decide on a maximum number of such evaluations, or function call budget,'' before undertaking the search. A function call budget is also useful when testing and comparing the various searches. Therefore, all searches have a maxCalls field, along with accessor and mutator functions so the user can set the call budget. By default, maxCalls is set to the flag value NO_MAX (currently, NO_MAX = -1), which means that the search will terminate based only on refinement of the design, no matter how many function calls have been made. A user who wishes to set a call budget may use SetMaxCalls(long calls) to do so.

In the PatternSearch class the default behavior, if maxCalls is not set to a user-specified value, is that the search will terminate only after a reduction in delta (the step length of the search) takes place, even if the function call count has somewhat exceeded maxCalls. In this case, the search will terminate following the first shrink step made after maxCalls has been reached.

On occasion, it may be either necessary or desirable to terminate the search immediately when the limit on the number of function calls has been reached. For those cases, I have provided the protected field exact_count, and the method SetMaxCallsExact(long calls). If the function call budget has already been set, the method SetExact() may be used instead. This will cause the search to terminate based on precisely the number of calls specified. The accessor bool IsExact() returns the value of exact_count. The protected utility method bool BreakOnExact() is also used in the implementations.

I have implemented this feature throughout the DirectSearch hierarchy. The default setting for exact_count for all members of PatternSearch class, as well as for the hybrid pattern/simplex search SMDSearch, is false; for the simplex searches NMSearch and SHHSearch the default is true. I have included the exact/inexact option for all searches primarily so that researchers can more meaningfully compare them. In actual optimization, however, the default settings are recommended. The simplex searches do not have the same theoretical basis for stopping only after a shrink as the pattern searches do and, in fact, the Nelder-Mead algorithm can converge without ever actually shrinking [17].

Stopping criteria in the simplex searches

The classic measure for stopping simplex searches is the standard deviation test suggested by Nelder and Mead [18], and used by Gurson in his original code [7].

In deference to the historical precedents of the Nelder-Mead and Spendley, Hext, and Himsworth simplex searches, I have tried as much as possible not to change the search algorithms themselves in my reworking of the classes. But, given a problem I ran into in testing SHHSearch, where use of the standard-deviation criterion resulted in premature convergence of the search, I thought it necessary to provide the user with a safer stopping criterion. The alternative test terminates based on the refinement of the simplex rather than on the standard deviation of the function values at the vertices in the simplex. To allow this option I added the boolean fields delta and Stop_on_std, and the methods Is_Stop_on_std(), Set_Stop_on_std(), and Set_Stop_on_delta() to the SimplexSearch class; and to NMSearch, NMdelta, an overloaded version of delta, along with the method CalculateNMdelta().

Details of the premature-convergence example, as well as a fuller discussion of stopping criteria in the simplex searches, can be found in Appendix A, with source code listings in Appendix C. Appendix A also discusses the Hybrid_NMSearch class, which was developed in response to the concerns about stopping criteria.

Protection and access control

In order to facilitate inheritance, but at the same time protect data from direct manipulation by the user, I made fields in the DirectSearch hierarchy protected rather than public (as they were in the original PatternSearch class) or private (as they were in the original three simplex searches). Although mindful of Stroustrup's [26] comment that

...declaring data members protected is usually a design error. Placing significant amounts of data in a common class for all derived classes to use leaves that data open to corruption.

I felt that, since I was refactoring old code rather than writing new code from scratch, the compromise was justified. I had at one point actually rewritten all the code for the pattern searches using accessor and mutator methods. The result was less-readable--and less efficient--code without much of a payoff in terms of protecting the data. The programs as I received them were designed so that the concrete classes could manipulate the fields directly. Rather than snarl up the code by calling accessor after accessor, I decided to make the fields protected.

I moved the ExploratoryMoves() method into the protected section in all cases. The actual algorithmic work of the searches is done here, and I saw no reason to keep those implementation details in the public part of the classes. The public interface for all the classes now has a method called BeginSearch(), which calls ExploratoryMoves() and does some housekeeping tasks as necessary. Hiding the details with a masking function should make it easier to modify the algorithms or their implementations later on, with little or no change in the public interface.

Simplifying the interface

Originally there were several variables that were assigned in other files (e.g objective.h) or simply by a #define statement. Because these were necessary for the proper execution of the searches, I brought them into the classes as fields. These new fields are:

• maxCalls
• stoppingStepLength
• initialStepLength

Wherever I could, I tried to reduce parameter lists in the public methods (except, of course, for constructors and some accessors and mutators). For example, as discussed above, the simplex initialization methods no longer take any parameters.

I made a few changes in auxiliary files as well. For example, I modified the vec.h file to include the one line we needed from subscript.h, in order to reduce the number of files. I took a version of Chris Siefert's Bound struct (used for constrained optimization) from his MAPS header file and put it into our objective.h files.

Some Implementation Considerations

Portability

At the time of this writing, C++ is still in a dynamic state. The compilers on which we tested this software were in varying stages of implementing the new ANSI C++ Standard[1]. I made a number of decisions based on that reality.

Allocation of memory with operator new

The code as I received it checked for the successful allocation of memory using assert statements (when it checked at all). Though this is not a problem with experimental code, I wanted to put more error-checking into the release version. Here, I ran into an unexpected snag.

Under the ANSI C++ standard [1,26] if operator new fails to allocate the requested amount of memory, it throws an exception rather than returning a null pointer. This means that, unless the programmer explicitly uses the (nothrow) version of new, using assert to check for failed allocation is no longer effective. What happens with a program compiled with standard-compliant compiler is that the system will, upon failure, look for an exception handler, then keep looking until it can look no further, at which point it will call terminate() and--by default--abort(). The appropriate way to deal with a failure of allocation with standard compilers is using try/catch syntax, e.g.:

    try{
myptr = new mytype;
} //try
catch ( bad_alloc exception ) {
...
// handle exception
...
} //catch


Unfortunately, many older C++ compilers do not recognize the try/catch construction. I could find no safe way of testing allocation that worked on compilers written both before and after the acceptance of the standard..

My solution to this awkward (though presumably temporary) situation is to use two sets of templates for memory allocation. By default, the standard try/catch templates will be used. If the user has an older compiler, however, compiling with the flag -DDOLD_ALLOC (or simply doing a #define DOLD_ALLOC) will make the switch to old-style allocators.

I put the templates (one for arrays, one for Vector objects, one for Matrix objects) into one file called Dyn_alloc.h. (I originally had two separate files, but Tony Padula, who was using DirectSearch for his krigifier software, put them into one to keep the file count down. It was a good idea, so I kept it.)

Templates are a better idea here than simple functions because they do the same work with less code: a single template to allocate arrays, for example, works for arrays of doubles, ints, or any simple type. Keeping all allocators in a single file keeps maintenance to a minimum. And keeping all those try-blocks and conditional compilation statements out of the class code makes it easier to read. When the C++ standard is more universally implemented by the compiler-writers, it will be easy to delete the old-style allocators and retire the -DDOLD_ALLOC flag.

The cmat.h and vec.h files

As mentioned earlier, we use Siefert's modified version [24] of Pozo's cmat.h and vec.h from the Template Numerical Toolkit [29]. I have further modified the two files by adding error-handling, as discussed above, for memory allocation. There were only a few instances where new was used in these files, so I just added the appropriate syntax inline, along with the necessary statements for conditional compilation.

We chose these Vector and Matrix templates over those in the Standard Template Library (STL) because, first, they were designed for linear algebra and work very well for our purposes, and second, because the implementation of the STL is not yet complete for all compilers, including the ones with which we have been working.

Namespaces

When these searches were first being written by Dolan and Gurson, namespaces were not supported on the version of gcc installed as our default compiler. Chris Siefert remedies the problem by modifying vec.h and cmat.h to remove the namespace designations. The situation is somewhat better now: namespaces are supported by our current version of gcc (we are using version 2.91.66 as I write this in April 2001), but that is not necessarily the case everywhere. Accordingly, this version of DirectSearch is being released without the use of namespaces. When namespaces are universally supported by compiler-writers, it will be easy enough to adapt future releases of the software to include it.

Efficiency

The changes I have made to bring all these searches together have concentrated on readability and maintainability, admittedly at the expense of some efficiency: inheritance is not cheap in terms of CPU time. I have added one additional level of inheritance for the pattern searches, and two additional levels for the three simplex searches.

A fundamental assumption underlying this entire line of software development is that the majority of the CPU time in a direct search will be taken up not by the search process per se but by the function evaluations themselves. Because of this, I felt the gain in programmer-friendliness provided by a more object-oriented approach more than outweighed the small relative reduction in efficiency that would result from it.

That said, I did make the decision, as discussed above, to make the fields protected rather than making them private and imposing upon the concrete classes the additional cost of using accessor and mutator functions every time they needed to get at the data. I would probably have made different design decisions if I had been starting from scratch rather than reworking given code.

Some further comments on SMDSearch

Gurson [7] devised the sequential multidirectional search (SMDSearch) as a single-processor version of multidirectional search [30,31]. A look at the code will show the reader that the SMDSearch class (derived from SimplexSearch), unlike NMSearch and SHHSearch, uses very little of its parent's code. SMDSearch is an oddity--a simplex search that is also a pattern search--and its implementation is quite different from the implementations of its sister classes. For example, SMDSearch reflects not through a centroid but through a point of the simplex; and it must maintain information on a reflection simplex as well as its primary simplex.

Additionally, the SMDSearch algorithm itself was designed by Gurson to be economical in its use of function evaluations. Whereas the Nelder-Mead [18] and Spendley, Hext, and Himsworth [25] searches evaluate every point of the simplex for each iteration, the greedier SMDSearch makes evaluations only until a new minimizer is found. Because of this, the SMDSearch class overwrites a number of SimplexSearch's methods. For example, SMDSearch does all its own simplex initialization in order to avoid the additional function evaluations that are performed in initializing simplices for the other two classes. Because it is a pattern search, SMDSearch has a delta field and ignores the centroid field that it inherits from SimplexSearch.

I considered having SMDSearch inherit from both the SimplexSearch and PatternSearch classes, but abandoned the idea, not only because multiple inheritance can be tricky to implement (and to maintain), but because Gurson originally wrote the search as a simplex search, and one of my goals in this project was to maintain whenever possible the original spirit of the software.

A few other additions

I added overloaded assignment operators for all classes, and additional accessors where needed. I also added PrintDesign() in the PatternSearch and SimplexSearch classes, first as a testing tool, but later as a way for the user to print out useful information about the search. I have incorporated a slightly modified version of Siefert's MultiCompassSearch class into the package, (as well as the Hybrid_NMSearch discussed above). I hope to add additional searches later.

At Siefert's request, I added an IDnumber field to all searches, as a way of identifying what type of search a pointer is pointing to, and the accessor methods PrintSearchType() and GetID() so the user could recover this information.

Acknowledgements

I would like to thank the many people and groups who helped me with this project. The National Science Foundation supported this project and the research underlying it (NSF Grant CCR-9734044). Liz Dolan and Adam Gurson wrote and tested the original code. Their work is the heart and soul of this software, and I thank them for giving me such excellent, well-though-out material to work with. Chris Siefert and Tony Padula, who were the principal users of the original searches, gave me vast amounts of time, help, and advice on how to improve and streamline the interface; and rewrote their own programs to accommodate the changes. Dr. Robert Michael Lewis and Dr. Paul Stockmeyer graciously agreed to serve on my 710 committee and were most helpful with details and explanations. Professors Bill Bynum, Robert Noonan, and Michael Trosset also gave advice. Dr. Bynum loaded gcc on William and Mary's RISC server, and Dr. Phil Kearns reloaded an old release of gcc on the CS department system, so I could test the new code on different compilers using different operating systems. Finally, my advisor, Dr. Virginia Torczon, gave more help and support than I can put into words.

Appendix A

Stopping Criteria in the Simplex Searches

During the testing of the DirectSearch version of the simplex searches, I ran into a problem with the SHH search. I had set it to to minimize the elliptical paraboloid , whose minimum is, of course, at the origin. The parameters were the default ones, with dimension = 2, a starting step length of 2.0, and a stoppingStepLength of 10e-8. The starting simplex was also the default, a regular right simplex. I made multiple runs, each using a different starting point, and found that on many of these runs the search would converge prematurely at the point (1,1). Eventually, I traced the problem to the Stop() method.

Gurson wrote the SHHSearch and NMSearch classes to terminate based on the standard deviation test originally recommended by Nelder and Mead [18]. These searches would terminate based on this classic test unless first stopped by exhaustion of a function call budget. The Nelder-Mead stopping criterion says to take the mean, , of the function values for all simplex vertices except the current best point, then halt the search if

 (1)

where () is the number of simplex points and is a small tolerance constant [2,7,18].

In the case of my paraboloid, for some of the starting points I used, e.g. (5,5), when the search terminated the simplex looked like this:

 Point: 1 -1 Value: 4 Point: -1 1 Value: 4 Point: 1 1 Value: 4

Because the vertices of the search simplex were all lying on the same level set of the objective function, all the terms being summed in the stopping criterion formula became zero, and the search halted there. The listings for the file fail.cc, along with the versions of objective.h and objective.cc that contain the appropriate function, can be found in Appendix B. The search set up by fail.cc will stall as described above.

This phenomenon was mentioned in Woods [33], and can happen with either the Nelder-Mead or the Spendley, Hext and Himsworth search, with any type of simplex. All that is required is this unfortunate alignment of the simplex vertices on, or very near, a level set of the objective function. The SimplexSearch field delta is analogous to the delta of the pattern searches, but in the case of the simplex searches contains information about the simplex edge lengths rather than the rational lattice'' of the pattern searches. It is for this reason that I did not raise the field up into the DirectSearch class. Granted, in the hybrid SMDSearch, delta means precisely what it does in any pattern search, but the meaning is overloaded to one extent or another with the other two searches currently in the SimplexSearch class.

As discussed in §5.3.2, I have added the boolean fields delta and Stop_on_std, and the methods Is_Stop_on_std(), Set_Stop_on_std(), and Set_Stop_on_delta() to the SimplexSearch class. To NMSearch I have added NMdelta, an overloaded version of delta, along with the method CalculateNMdelta().

In SHHSearch, delta is the length of the longest side of the design simplex. This is intuitive and easy to track and update. Whenever we shrink the simplex, we reduce delta by the same coefficient. This is an obvious extension to the use of in the traditional pattern searches since the shape of the initial simplex is not deformed.

With a Nelder-Mead search, however, this would not be an effective extension. The distinguishing feature of this search is the way the simplex deforms during the progress of the search. One given iteration's longest side may well be the next iteration's shortest; and a Nelder-Mead search may converge without ever actually shrinking the way an SHH search would. It was necessary to come up with a different measure of delta'' for the Nelder-Mead simplex algorithm.

My solution to this problem was to overload the concept of delta'' with the field NMdelta, which represents the mean of the lengths of all edges of the simplex. Unfortunately, finding this value is relatively expensive ( including the underlying vector operations). To minimize the algorithmic cost, I made the stopping criterion for NMSearch a hybrid of the standard-deviation'' and delta'' approaches. At the beginning of the search, the standard-deviation test is used whether the user chooses to set Stop_on_std to true'' or to use the default of stopping on delta. When, later in the search process, this test indicates that the search should be halted, a boolean variable called stopBool is set to true.'' The program enters the following if-clause:

    if (!Stop_on_std && stopBool) {
CalculateNMDelta();
stopBool =  NMdelta < stoppingStepLength;
toleranceHit = int(stopBool);
}
return stopBool;


This way, the faster () standard-deviation test does most of the checking, and the more expensive test of NMdelta is used only when the function values are very close to each other. (This work-around is unnecessary for SHHSearch, where delta is trivial to calculate.) As another cost-cutting measure, although SHHSearch's delta is adjusted automatically whenever the simplex shrinks, we calculate NMSearch's NMdelta only when truly necessary; i.e. inside the Stop() method as above, and when the user calls GetDelta(). (The NMSearch version of GetDelta() has been overloaded to calculate and return the value of NMdelta).

Given these modifications, stopping based on delta is now the default behavior for all the searches, with the exception of NMSearch, where there is no obvious analog to the delta that plays such a critical role in pattern searches. If a user of SHHSearch prefers to use the standard Nelder-Mead stopping criterion, a call to Set_Stop_on_std() is all that is necessary. Setting Set_Stop_on_std() for an SMDSearch has no effect; it will stop based on delta in any case.

Hybrid_NMSearch class

As another approach to this problem, I implemented the Hybrid_NMSearch class. In this search, both an NMSearch object and an EdHJSearch object are constructed. The search begins with an NMSearch and proceeds until the suggested Nelder-Mead stopping criterion is met. Then, if there is room in the function call budget and delta has not been reached, the EdHJSearch search is constructed using the end state of the NMSearch as a starting point. The EdHJSearch finishes the search, terminating based on maxCalls (if applicable) and EdHJSearch::delta. Note that exact_count is set to true'' here by default. If the recommended procedure of stopping only after a shrinking step is preferred, the user must call SetMaxCalls().

The purpose of this arrangement is to avoid the expense of calculating NMdelta and at the same time give the user the option of using the safer delta'' test. It has the added advantage of giving the search the convergence guarantee of a Pattern Search, which the Nelder-Mead search does not enjoy.

Appendix B

Source code listing for main_obj.cc, trigclass.h, and trigclass.cc.

/*
main_obj.cc
example main program using an HJSearch and long constructor

Use with trigclass.h

Anne Shepherd, 8/2000 at
The College of William and Mary, Williamsburg, Virginia,
under advisor Dr. Virginia Torczon

to compile use:
g++ -g -Wall main_obj.cc objective.cc trigclass.cc DirectSearch.cc \
PatternSearch.cc  HJSearch.cc  -lm -o main_obj

*/

#include "objective.h"
#include "HJSearch.h"
#include "trigclass.h"
#include <iostream.h>              // for cout

int main(void)
{
long n = 1;                      // number of variables (dimension of
//                      the search)
double startVal = 1.0;           // starting point for x
double startstep = .25;          // starting step length
double endstep = 10e-8;          // ending step length

/* here, we declare and initialize a pointer to  an object of
* type trigclass.
*/

trigclass* t_ptr = new trigclass(0);

Vector<double> Hminimum_1(n);      // to store the first minimum point later

/* we'll initialize a one-entry Vector whose value is startVal,
* and use it as our starting point.
*/
Vector<double>minVec(n, startVal);

/* now construct an HJSearch object with the long
* constructor.
*
* In this example we minimize the cosine. find_Cos is the name of
* a "friend" function of trigclass.
*/
HJSearch HJ_1(n, minVec, startstep, endstep, find_Cos, (void *)t_ptr);

double HMinVal_1;
long Hcalls_1;
double step_1;

/* start searching */
HJ_1.BeginSearch();

/* retrieve information about HJ_1 */
HJ_1.GetMinPoint(Hminimum_1);
HJ_1.GetMinVal(HMinVal_1);
Hcalls_1 = HJ_1.GetFunctionCalls();
step_1 = HJ_1.GetDelta();

cout << "\nMinimum point found: \n" << Hminimum_1;
cout << "Value: \n" << HMinVal_1 << "\n\n";
cout << "Last used step length:\n" << step_1 * 2 << " in "
<< Hcalls_1 << " function calls.\n" << endl;

delete t_ptr;
return 0;

}//main


/*
trigclass.h

a toy class to illustrate the long DirectSearch class constructor.

Anne Shepherd
College of William and Mary

*/

#ifndef _toy_trig_
#define _toy_trig_

#include <iostream.h>
#include <math.h>
#include "vec.h"

class trigclass
{
public:

trigclass(double x);

void Sin_fcn(long vars, Vector<double> &v, double & f);

void Cos_fcn(long vars, Vector<double> &v, double & f);

void Set_x_val(double new_x);

double Get_y_val();

/* the following two methods are friends rather then members because
* we want our main program to be able to pass pointers to them.
*
* Note that in the calls masked by these methods, the void pointer
* should be cast to a pointer of the class type.
*/

friend void find_Sin(long vars, Vector<double> &v, double & f,
bool & success, void* t_ptr)
{
(*(trigclass*)t_ptr).Sin_fcn(vars, v, f);
success = true;
}

friend void find_Cos(long vars, Vector<double> &v, double & f,
bool & success, void* t_ptr)
{
(*(trigclass*)t_ptr).Cos_fcn(vars, v, f);
success = true;
}

private:

double mySin(Vector<double> & vec1);

double myCos(Vector<double> & vec1);

/* We don't really use this field in this toy class. */
double _x;

\begin{verbatim}
};
#endif


/* trigclass.cc

definitions for toy trig class

Anne Shepherd
W&M 8/00

*/

#include "trigclass.h"

//public:

trigclass::trigclass(double x) : _x(x)
{
}

// Note that these two friend functions have the necessary signature for
// our searches.

void trigclass::Sin_fcn(long vars, Vector<double> &v, double & f)
{
if (vars != 1) {
cerr << "\nError: Dimension must be 1!!\n";
exit(1);
}
f = mySin(v);
return;
}

void trigclass::Cos_fcn(long vars, Vector<double> &v, double & f)
{
if (vars != 1) {
cerr << "\nError: Dimension must be 1!!\n";
exit(1);
}
f = myCos(v);
return;
}

void trigclass::Set_x_val(double new_x)
{
_x = new_x;
}

//private:

double trigclass::mySin(Vector<double> & vec1)
{
double scalar = vec1.begin()[0];
return sin(scalar);
}

double trigclass::myCos(Vector<double> & vec2)
{
double scalar = vec2.begin()[0];
return cos(scalar);
}


Appendix C

Source code listing for fail.cc and associated objective files.

/*
fail.cc
example main program using an SHHSearch. This search will terminate
prematurely.

P.L. (Anne) Shepherd, 8/2000 at
The College of William and Mary, Williamsburg, Virginia,
under advisor Dr. Virginia Torczon

to compile use:

g++ -g -Wall fail.cc objective.cc DirectSearch.cc SimplexSearch.cc \
SHHSearch.cc  -lm -o fail

*/

#include "objective.h"
#include "SHHSearch.h"
#include <iostream.h>              // for cout
#include "vec.h"

int main()
{
long n = 2;                      // number of variables (dimension of
//                      the search)

/* terminates prematurely with startval = 1, 3, 5, 7, 9 ...
finds true min with startval = 2, 4, 6, 8, 10 ...
*/
double startVal = 3.0;          // starting point for x

/* if we change this to 2.01 it doesn't hang up. */
double startstep = 2.0;          // starting step length
double endstep = 10e-8;          // ending step length

Vector<double> Sminimum(n);      // to store the minimum point later

/* we'll initialize an n-entry Vector whose value is startVal,
* and use it as our starting point.
*/
Vector<double>minVec(n, startVal);
cout << "Starting point is: " << minVec << endl;

/* now construct a Search object.  Note that parab_fcn is a
* wrapper for the function  (x^2 + 3y^2). The last argument
* is set to NULL because we won't be using it.
*/
SHHSearch SHH(n, minVec, startstep, endstep, parab_fcn, NULL);

double SMinVal;
long Scalls;
SHH.SetMaxCalls(250);

/* If this flag is set, the search will terminate prematurely
* on (1,1).
*/
SHH.Set_Stop_on_std();

cout << "delta = " << SHH.GetDelta() << endl;
cout << "\nStop_on_std status  = " << SHH.Is_Stop_on_std() << endl;

/* start searching */
SHH.BeginSearch();
SHH.PrintDesign();
SHH.GetMinPoint(Sminimum);
SHH.GetMinVal(SMinVal);
Scalls = SHH.GetFunctionCalls();

cout << "\nMinimum point found: " << Sminimum;
cout << "\nValue: \n" << SMinVal <<  " in ";
cout << Scalls << " function calls.\n" << endl<< "\n\n";
cout << "delta = " << SHH.GetDelta() << endl;

return 0;
}//main


//objective.h
// this example file declares the functions as required
// Liz Dolan, P.L. (Anne) Shepherd, College of William and Mary

#if !defined _userfile_
#define _userfile_

#include "vec.h"
// This can be used for fcn_eval() or EvalF();
// extern char Fname[120];

//Chris Siefert's Bound struct--used for constrained problems
struct Bound {
Vector<double> lower;
Vector<double> upper;
};

void fcn(long, Vector<double> &, double &, bool &, void *);

void fcn_2(long, Vector<double> &, double &, bool &, void *);

void fcn_eval(long vars, Vector<double> &x, double & f, bool & flag,
void * nothing);

void parab_fcn(long, Vector<double> &, double &, bool &, void *);

double mySin(Vector<double> &);

double myCos(Vector<double> &);

double eval_parab(Vector<double> &);

void mcKinnon(long vars, Vector<double> & x, double & f,
bool & flag, void * nothing);

void powell(long vars, Vector<double> &x,
double & f, bool & flag, void* nothing);
#endif

/* objective.cc
*
* Simple examples of how to wrap--or otherwise connect--an
*   objective function so it can work with the DirectSearch
*   class.
*
* Anne Shepherd, 8/2000 (Thanks to Chris Siefert for the EvalF
*   function).  Revised, pls, 1/2001
*/

#include "objective.h"
#include "vec.h"
#include <math.h>
#include <iostream.h>
#include <unistd.h>
#include <sys/wait.h>
#include <sys/types.h>
#include <stdio.h>

double EvalF(const Vector<double> &point, bool &success, void * VFname);

/* Here, we use fcn to wrap the user's
* objective function.
*/
void fcn(long vars, Vector<double> &x, double & f, bool & flag, void* nothing)
{
if (vars < 1) {
cerr << "\nError: Dimension cannot be less than 1!!\n";
exit(1);
}
f = mySin(x);
flag = true;
return;
}

void fcn_2(long vars, Vector<double> &x, double & f,
bool & flag, void* nothing)
{
if (vars < 1) {
cerr << "\nError: Dimension cannot be less than 1!!\n";
exit(1);
}
f = myCos(x);
flag = true;
return;
}

/**
* This uses the EvalF() function to call an executable.  It gets its
* "f" value from the stdout of the objective function.  See user's manual.
* void * VFname is used to send in the name of the executable.
*/

void fcn_eval(long vars, Vector<double> &x, double & f, bool & flag,
void * VFname)
{
Vector<double> vec(x);
f = EvalF(vec, flag, VFname);
}

void parab_fcn(long vars, Vector<double> & x, double & f,
bool & flag, void * nothing)
{
if (vars != 2) {
cerr << "\nError: Dimension must be 2 !!\n";
exit(1);
}
f = eval_parab(x);
flag = true;
return;
}

/* see K.I.M. McKinnon, Convergence of the Nelder-Mead Simplex
* Method to a Nonstationary Point,'' Siam J. Optim., Vol 9, No. 1,
* pp. 148-158. (1998)
*
*  vars must equal 2
*
*  Note that this function is not wrapped--I just wrote it to
*  take the parameters required by the searches.
*/
void mcKinnon(long vars, Vector<double> & x, double & f,
bool & flag, void * nothing)
{
double my_x = x[0];
double my_y = x[1];
if (my_x <= 0) {
f = (360 * (my_x * my_x)) + my_y + (my_y * my_y);
}
else {
f = (6 * (my_x * my_x)) + my_y + (my_y * my_y);
}
flag = true;
}

/** Powell's quartic function, as discussed in Nelder and Mead's 1965
* paper.
*
*    y= (x_1 + 10*x_2)^2 + 5*(x_3 - x_4)^2
*             + (x_2 - 2*x_3)^4 + 10*(x_1 - x_4)^4
*
*  Minimum value is 0. vars must equal 4;
*
*  Nelder and Mead tested this with starting point (3, -1, 0, 1).
*/
void powell(long vars, Vector<double> &x,
double & f, bool & flag, void* nothing)
{
if (vars != 4) {
cerr << "\nThis function required four variables. ";
cerr << "Exiting with value 1.\n";
exit(1);
}
double my_x1 = x[0];
double my_x2 = x[1];
double my_x3 = x[2];
double my_x4 = x[3];

double first, second, third, fourth;

first = (my_x1 + 10 * my_x2) * (my_x1 + 10 * my_x2);
second = 5 * (my_x3 - my_x4) * (my_x3 - my_x4);
third = (my_x2 - (2 * my_x3)) * (my_x2 - (2 * my_x3));
third *= third;
fourth = (my_x1 - my_x4) * (my_x1 - my_x4);
fourth *= fourth;
fourth *= 10;

f = first + second + third + fourth;

flag = true;
}

/*  We'll minimize the sine of a number--here vec1
*  is a Vector of one member. (Or even if it isn't, we just
*  take the sine of the first element.)
*/
double mySin(Vector<double> & vec1)
{
double scalar = vec1.begin()[0];
return sin(scalar);
}

/*
* and here we find the cosine
*/
double myCos(Vector<double> & vec2)
{
double scalar = vec2.begin()[0];
return cos(scalar);
}

/*
* here, it's an elliptical paraboloid--min at
* the origin. n must equal 2.
*/
double eval_parab(Vector<double> & parab_vec)
{
double my_x = parab_vec[0];
double my_y = parab_vec[1];

return ( (my_x * my_x) + ( 3 * (my_y * my_y) ) );
}

#define ERROR HUGE_VAL

/****
* Evalf()
*   written by Chris Siefert, slightly modified by Anne Shepherd
*
* INPUT: The point at which to evaluate the function; a boolean flag;
*     the file name.
* OUTPUT: The function value.
* EFFECT: This function will do the fork/exec dance, sending the
*     program the parameters on its stdin, then after waiting, it will
*     get a  double from the  program's stdout.  If the program exited
*     with an error, it will return ERROR and set success to false.
*     Otherwise success will be true.
****/

double EvalF(const Vector<double> &point, bool &success, void * VFname) {
int childpid, pipe1[2], pipe2[2];
const char* Fname;
long i;
long P = point.dim();
double ret_val;
if (VFname != NULL) {
Fname = (const char *) VFname;
}
FILE *writepipe;

long state;

if ((pipe(pipe1) < 0) || (pipe(pipe2) < 0) ) {
perror("pipe");
exit(-1);
}/*end if*/

if ((childpid = fork()) < 0) {
perror("fork");
exit(-1);
} else if (childpid > 0) {   /*Parent*/
close(pipe1[0]); close(pipe2[1]);

/* Write to child on pipe1[1], read from child on pipe2[0]. */
writepipe=fdopen(pipe1[1],"w");
{perror("fdopen");cerr<<"MPEF   : Cannot open pipes... exiting\n";abort();}

/*Output the vector x*/
for(i=0;i<P;i++) fprintf(writepipe,"%1.8f ", point[i]);
fprintf(writepipe,"\n");
fflush(writepipe);

/*Read in the returned value, and WAIT*/
wait(&state);
/*Close outstanding pipes*/

if(WIFEXITED(state)&&(WEXITSTATUS(state)==0)) {
/*normal exit - no abort, return code 0*/
success=true;
return ret_val;
}/*end if*/
else {
success=false;
success = true;
return ERROR;
}/*end else*/

} else { /*Child*/
close(pipe1[1]); close(pipe2[0]);
/* Read from parent on pipe1[0], write to parent on pipe2[1]. */
dup2(pipe1[0],0); dup2(pipe2[1],1);
close(pipe1[0]); close(pipe2[1]);

if (execlp(Fname, Fname, NULL) < 0) {
perror("execlp");
abort();
}
return ERROR;
/* Never returns */
}/*end else*/

return ERROR;  /*never returns*/
}/*end EvalF*/


Bibliography

1
ANSI.
ANSI C++ standard, November 1997.
ISO/IEC FDIS 14882.

2
Mordecai Avriel.
Nonlinear Programming: Analysis and Methods.
Prentice-Hall, Englewood Cliffs, NJ, 1976.

3
H.M Deitel and P.J. Deitel.
C++: How to Program.
Prentice Hall, 1994.

4
Elizabeth D. Dolan.
Pattern search behavior in nonlinear optimization.
Honors Thesis, Department of Computer Science, College of William & Mary, Williamsburg, Virginia 23187-8795, May 1999. Accepted with Highest Honors.

5
Elizabeth D. Dolan, Robert Michael Lewis, and Virginia J. Torczon.
On the local convergence properties of pattern search.
In review for release as an ICASE technical report. Submitted to SIAM Journal on Optimization.

6
Martin Fowler.
Refactoring: Improving the Design of Existing Code.

7
Simplex search behavior in nonlinear optimization.
Honors Thesis, Department of Computer Science, College of William & Mary, Williamsburg, Virginia 23187-8795, May 2000. Accepted with Highest Honors.

8
Robert Hooke and T. A. Jeeves.
Direct search solution of numerical and statistical problems.
Journal of the Association for Computing Machinery, 8(2):212-229, April 1961.

9
Danny Kalev.
The ANSI/ISO C++ Professional Programmer's Handbook.
Que Corporation, second edition, 1999.

10
Brian W. Kernighan and Rob Pike.
The Practice of Programming.

11
Jeffrey C. Lagarias, James A. Reeds, Margaret H. Wright, and Paul E. Wright.
Convergence properties of the Nelder-Mead simplex method in low dimensions.
SIAM Journal on Optimization, 9(1):112-147, 1998.

12
Robert Michael Lewis, Virginia Torczon, and Michael W. Trosset.
Why pattern search works.
Optima, (59):1-7, October 1998.

13
Robert Michael Lewis, Virginia Torczon, and Michael W. Trosset.
Direct search methods, then and now.
Journal of Computational and Applied Mathematics, 2000.

14
Robert Michael Lewis and Virginia J. Torczon.
Rank ordering and positive bases in pattern search algorithms.
Technical Report 96-71, Institute for Computer Applications in Science and Engineering, Mail Stop 132C, NASA Langley Research Center, Hampton, Virginia 23681-2199, 1996.

15
Steven A. Maguire.
Writing Solid Code.
Microsoft Press, 1993.

16
Steven C. McConnell.
Code Complete.
Microsoft Press, 1993.

17
K.I.M. McKinnon.
Convergence of the Nelder-Mead simplex method to a nonstationary point.
SIAM Journal on Optimization, 9(1):148-158, December 1998.

18
J. A. Nelder and R. Mead.
A simplex method for function minimization.
The Computer Journal, 7(4):308-313, January 1965.

19
Interpolation and pseudorandom function generation.
Honors Thesis, Department of Mathematics, College of William & Mary, Williamsburg, Virginia 23187-8795, May 2000. Accepted with High Honors.

20
Jo Ellen Perry and Harold D. Levin.
An Introduction to Object-Oriented Design in C++.

21
Elijah Polak.
Computational Methods in Optimization: A Unified Approach.

22
Christopher M. Siefert.
Model-assisted pattern search.
Honors Thesis, Department of Computer Science, College of William & Mary, Williamsburg, Virginia 23187-8795, May 2000. Accepted with Highest Honors.

23
Christopher M. Siefert, Virginia Torczon, and http://www.cs.wm.edu/ va/software/maps Michael Trosset.
Model-Assistyed Pattern Search (maps) source code.

24
Christopher M. Siefert and http://www.cs.wm.edu/ va/software/Vectors/ Virginia Torczon.
Vector and matrix classes.

25
W. Spendley, G. R. Hext, and F. R. Himsworth.
Sequential application of simplex designs in optimisation and evolutionary operation.
Technometrics, 4(4):441-461, November 1962.

26
Bjarne Stroustrup.
The C++ Programming Language.
Addison-Wesley, special edition, 2000.

27
P.L. Shepherd. http://www.cs.wm.edu/plshep.
Class Documentation for DirectSearch and its derived classes.

28
P.L. Shepherd. http://www.cs.wm.edu/plshep.
User's Manual for DirectSearch and its derived classes.

29
Roldan Pozo. http://www.math.nist.gov/tnt.
Template numerical toolkit: A numeric library for scientific computing in C++.

30
Virginia Torczon.
Multi-Directional Search: A Direct Search Algorithm for Parallel Machines.
PhD thesis, Department of Mathematical Sciences, Rice University, Houston, Texas, 1989; available as Tech. Rep. 90-07, Department of Computational and Applied Mathematics, Rice University, Houston, Texas 77005-1892.

31
Virginia Torczon.
On the convergence of the multidirectional search algorithm.
SIAM Journal on Optimization, 1(1):123-145, February 1991.

32
Virginia Torczon.
On the convergence of pattern search algorithms.
SIAM Journal on Optimization, 7(1):1-25, February 1997.

33
Daniel J. Woods.
An Interactive Approach for Solving Multi-Objective Optimization Problems.
PhD thesis, Rice University, 1985.

34
Roland Wunderling and http://www.zib.de/Visual/software/doc++/index.html Malte Zockler.
Doc++: A documentation system for C/C++ and Java.

Design Document for DirectSearch and its Derived Classes

This document was generated using the LaTeX2HTML translator Version 99.2beta6 (1.42)

Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.

The command line arguments were:
latex2html -antialias -no_navigation -white -split 0 -t DirectSearch -image_type gif design.tex

The translation was initiated by Anne Shepherd on 2001-05-17

Anne Shepherd 2001-05-17