objective.cc

在linux下面实现的单纯性算法的源代码

/* objective.cc
 * simple quadratic equation in four dimensions with rigged starting point
 */

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

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

/* I defined my stoppingStepLength to be tolerance, which I must define here.
 * Note that the value chosen is machine/compiler dependent as approximately
 *   the square root of machine epsilon
 */
//double tolerance = sqrt(fabs(3.0 * (4.0/3.0 -1.0) -1.0));

/* Initial trial step length, reflecting the refinement of the search lattice
 */
double initialStep = 2.0;

/* Returns initial step length for the search.
 * This function is not technically necessary, as long as initialStepLength 
 * is set to a real number in the header file.
 */
double getInitialStep()
{
  return initialStep;
}

/* The user's function to be minimized.
 * vars is the number of variables in the problem.
 * double * x is a pointer to an array of doubles,
 *   which represent the point to be evaluated
 * double & f will be the value of the function
 *   at that point upon return
 * flag will be 1 if the function can be successfully
 *   evaluated, 0 otherwise
 * Note that pattern searches have enough flexibility
 *   that they need not necessarily be able to get the
 *   value of all possible points
 */
void fcn(int vars, double * x, double & f, int & flag)
{
  //  if(vars!=n) {flag=0; cout << "PROBLEM IN fcn()\n"; return;}
  
  // double tone = 0.0;
  /*for(int i = 0; i < vars; i++)
    {
      tone += (x[i]-4.0) * (x[i]-4.0) + 10;
      }*/
  //  f = (x[0]-2.54)*(x[0]-2.54) + (x[1]+7.5)*(x[1]+7.5) + (x[2]+1.22)*(x[2]+1.22) + (x[3]-6.34)*(x[3]-6.34)+ (x[4]+8.78)*(x[4]+8.78) + 239.78984346;
  //  flag = 1;
  //  return;
      Vector<double> vec(vars);
      for (int i = 0; i < vars; i++ )
         vec[i] = x[i];
      f = EvalF(vec, (bool&)flag);
}

/* initpt returns the user's initial guess
 * vars is the number of variables of the problem
 * double * x is a pointer (NULL),
 *   which will point to a newly allocated array
 *   of doubles that represent the initial point
 */
void initpt(int vars, double *&x)
{
  double grace = 1;
  x = new double[vars];
  for(int i = 0; i < vars; i++)
    {
      x[i] = grace;
    }

}

#define ERROR HUGE_VAL

double EvalF(const Vector<double> &point, bool &success) {
/****
INPUT: The point at which to evaluate the function, a boolean flag.
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.
****/
  int childpid, pipe1[2], pipe2[2];
  long i;
  int P = point.dim();
  double ret_val;
  FILE *readpipe;
  FILE *writepipe;
  
  int 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]. */
    readpipe=fdopen(pipe2[0],"r");
    writepipe=fdopen(pipe1[1],"w");
    if(writepipe==NULL || readpipe==NULL)
      {perror("fdopen");cerr<<"MPEF   : Cannot open pipes... exiting\n";abort();}
    setlinebuf(writepipe);/*important!  This way it won't just hang*/

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

    /*Read in the returned value, and WAIT*/
    fscanf(readpipe,"%lf", &ret_val);    
    wait(&state);
    /*Close outstanding pipes*/
    fclose(writepipe);fclose(readpipe);
    //    close(pipe2[0]);close(pipe1[1]);
    
    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*/