Golden Redux ------------ 1. Changes in the gridpoints of Golden when the endpoints of the search interval are changed. 2. Remaking the Golden program to more efficiently use memory, with pointers. ------------------------------------------------------------------------ 1. Changes in the gridpoints of Golden when the endpoints of the search interval are changed. This shouldn't have surprised me when it happened in class (but it did!). The explanation is simply that the grid points are calculated as x_i = a + i*(b-a)/(n-1.0) [ i = 0..19; n=20 in the program ] When a = -1.0, b=+3.0, we get x_i = -1.0 + 2.0*i/(19.0) [ i = 0..19; n=20 in the program ] When a = -2.0, b=+2.0, we get x_i = -2.0 + 2.0*i/(19.0) [ i = 0..19; n=20 in the program ] These two sequences could only have values in common if for some value of i, 2.0*i/(19.0) = 1. But this cannot happen, for any integer i<19. ------------------------------------------------------------------------ 2. Remaking the Golden program to more efficiently use memory, with pointers. I mentioned in class a while back that golden.c made very inefficient use of the memory it allocated for Fibonacci numbers and for function evaluation. Here is a slightly modified version "golden.c" that allows us at the end to see how much storage was allocated by the program, and how much was needed, to store the Fibonacci numbers and the function evaluations. /******************************************************************** golden.c Compile: gcc -o golden golden.c -lm Golden Section Search for maximum of a unimodal function f(x) on the interval [a,b]. Specify: f(x), a, b, and n=desired number of pts of grid. Routine rounds n up to the nearest Fibonacci number, and uses that many grid points. Limits: From the Fibonacci relationship F_k = F_{k-1} + F_{k-2} (and starting with F_1 = F_2 = 1) we can conclude (start with n even, and extend to case of n s.t. n+1 even): F_k > 2^(ceiling(k/2)-1) Setting n = rhs, get log(n)/log(2) = k/2-1, approx, or k = 2*(log_2(n) + 1). For specificity we take NMAX = 1024; FMAX = 22; Cleverer C programming, using malloc(), would let the limit be implicitly defined by the amount of available memory. ********************************************************************/ #include #include #define NMAX 1024 #define FMAX 22 typedef double real; int itmp = 0; /********************************************************************/ /* WORK ARRAYS */ int eval[NMAX]; real y[NMAX]; int fib[FMAX]; /********************************************************************/ /* DEFINE NUMBER OF GRID POINTS AND INTERVAL ENDPOINTS HERE */ int n = 18; real a = -1.0, b = +3.0; /*********************************************************************/ /* DEFINE FUNCTION TO BE MAXIMIZED HERE */ real f(real x) { real y; y = exp(-pow(x,2.0)/2.0); itmp++; return(y); } /*******************************************************************/ int init_fib(void) { int i=1; fib[0] = fib[1] = 1; while (fib[i]<=n) { i = i + 1; if (i>=FMAX) { return(0); } fib[i] = fib[i-1] + fib[i-2]; } return(i); } /********************************************************************/ void init_arrays(void) { int i; for (i=0;i 1.5*(b-a)/n ) { diff = evaluate(fib[k-2],n,offset) - evaluate(fib[k-1],n,offset); printf("Checked y[%d] = %f vs. y[%d] = %f; \n", fib[k-2]+offset, evaluate(fib[k-2],n,offset), fib[k-1]+offset, evaluate(fib[k-1],n,offset)); if (diff<0.0) { aprime = X(fib[k-2],n,offset); offset = offset + fib[k-2]; k = k - 1; } else if (diff>0.0) { bprime = X(fib[k-1],n,offset); /* offset doesn't move */ k = k - 1; } else { aprime = X(fib[k-2],n,offset); bprime = X(fib[k-1],n,offset); offset = offset + fib[k-2]; k = k - 3; } printf("Resetting aprime = %f, bprime = %f\n\n", aprime,bprime); /* for (i=0; i Thus at some point, for example, the list has the form NULL <- (5, 0.998616) <-> (8, 0.791305) <-> (13, 0.221284) -> NULL A pointer variable "eval" points to the most recently examined entry in this list. For example, * if "eval" is pointing to the entry (5, 0.998616), and we want to see what the value is at i=13, we have to step "eval" two entries to the right. * if "eval" is currently pointing to (5, 0.998616), and we want to see what the value is at i=6, we need to evaluate the function, obtaining the new entry (6, 0.965967), which we need to insert in the list (between the i=5 and i=8 entries), for future use. etc. The predefined C constant NULL indicates "the end of the list" and keeps us from falling off one end or the other as we are looking up entries. Here is golden-list.c /******************************************************************** golden-list.c Compile: with a makefile that has LDFLAGS = -lm can just say "make golden-list" Golden Section Search for maximum of a unimodal function f(x) on the interval [a,b]. Specify: f(x), a, b, and n=desired number of pts of grid. Routine rounds n up to the nearest Fibonacci number, and uses that many grid points. This implementation uses a "doubly-linked list" to keep track of values of f(x) that have already been computed. ********************************************************************/ #include #include typedef double real; /********************************************************************/ /* DEFINE NUMBER OF GRID POINTS AND INTERVAL ENDPOINTS HERE */ int n = 18; real a = -1.0, b = +3.0; /********************************************************************/ int itmp = 0; /* total number of function evaluations */ int nmax = 0; /* number of nodes of the "eval" linked list */ /********************************************************************/ /* WORK SPACE */ int *fib; /* will be a dynamically allocated array */ typedef struct f_eval_self { int i; real y; struct f_eval_self *prev; struct f_eval_self *next; } f_eval; f_eval *eval; /* function values will be kept here */ /*********************************************************************/ /* DEFINE FUNCTION TO BE MAXIMIZED HERE */ real f(real x) { real y; y = exp(-pow(x,2.0)/2.0); itmp += 1; /* or itmp++ or itmp=itmp+1 */ return(y); } /*******************************************************************/ int init_fib(void) { /* From the Fibonacci relationship F_k = F_{k-1} + F_{k-2} (and starting with F_1 = F_2 = 1) we can conclude (start with n even, and extend to case of n s.t. n+1 even): F_k > 2^(ceiling(k/2)-1) Setting n = rhs, get log(n)/log(2) = k/2-1, approx, or k = 2*(log_2(n) + 1). */ int i, imax; imax = ceil(2.0*log(n)/log(2.0) + 1.0); fib = (int *) malloc ( imax * (sizeof(int)) ); if (fib==NULL) { printf("Cannot allocate enough space for Fibonacci numbers!\n"); exit(1); } i = 1; fib[0] = fib[1] = 1; while (fib[i]<=n) { i = i + 1; if (i>=imax) { return(0); } fib[i] = fib[i-1] + fib[i-2]; } return(i); } /*********************************************************************/ real X(int i, int n, int offset) { real x; x = a + (b-a) * ((i+offset)/(n-1.0)) ; return(x); } /*********************************************************************/ f_eval *new_eval (void) { f_eval *tmp; tmp = (f_eval *) malloc (sizeof(f_eval)); if (tmp==NULL) { printf("Cannot allocate space to grow the evaluation list!\n"); exit(2); } nmax++; return(tmp); } /*********************************************************************/ f_eval *insert_before(f_eval *eval) { f_eval *tmp; tmp = new_eval(); /* ALTERNATE NOTATION: */ (*tmp).next = eval; /* tmp->next = eval */ (*tmp).prev = (*eval).prev; /* tmp->prev = (*eval).prev; */ if ((*eval).prev!=NULL) { (*(*eval).prev).next = tmp; /* eval->prev->next = tmp; */ } (*eval).prev = tmp; /* eval->prev = tmp; */ return(tmp); } /*********************************************************************/ f_eval *insert_after(f_eval *eval) { f_eval *tmp; tmp = new_eval(); (*tmp).prev = eval; (*tmp).next = (*eval).next; if ((*eval).next!=NULL) { (*(*eval).next).prev = tmp; } (*eval).next = tmp; return(tmp); } /*********************************************************************/ real evaluate(int i, int n, int offset) { int ii; f_eval *tmp; /* just for some debug code at the end */ ii = i + offset; /* THIS FUNCTION CHANGED FROM THE FOLLOWING 6-LINER if(eval[ii]!=1) { y[ii] = f(X(i,n,offset)); eval[ii] = 1; } return(y[ii]); TO THE FOLLOWING SEARCH-AND-INSERT ROUTINE */ if (ii>(*eval).i) { while((*eval).next!=NULL) { eval = (*eval).next; if (ii<=(*eval).i) break; } } else if (ii<(*eval).i) { while((*eval).prev!=NULL) { eval = (*eval).prev; if (ii>=(*eval).i) break; } } if (ii>(*eval).i) { eval = insert_after(eval); (*eval).y = f(X(i,n,offset)); (*eval).i = ii; } else if (ii<(*eval).i) { eval = insert_before(eval); (*eval).y = f(X(i,n,offset)); (*eval).i = ii; } /* else if (ii==(*eval).i) { do nothing } */ /* FOLLOWING CAN BE USED TO PRINT CURRENT LINKED LIST */ /* tmp=eval; while((*tmp).prev!=NULL) { tmp = (*tmp).prev; } printf("NULL--"); do { printf("(%d, %f)--", (*tmp).i, (*tmp).y); tmp = (*tmp).next; } while (tmp!=NULL); printf("NULL\n"); */ /* PREVIOUS CAN BE USED TO PRINT CURRENT LINKED LIST */ return((*eval).y); } /*********************************************************************/ int main(void) { int i; int k,kmax; int offset = 0; /* essentially, the index for aprime in the grid */ real aprime, bprime; real diff; int oldn; /* Fill Fibonacci #'s up to n, and round n upward to nearest Fib # */ if ( (kmax = init_fib()) == 0 ) { printf("\n\nFibonacci array not big enough!!\n\n"); return(1); } oldn = n; k = kmax; n = fib[k]-1; /* Initialize endpoints and function evaluation list */ aprime = a; bprime = b; eval = new_eval(); (*eval).i = fib[k-2]; (*eval).y = f(X(fib[k-2],n,offset)); (*eval).prev = NULL; (*eval).next = NULL; /* NOTHING BELOW THIS COMMENT CHANGES IN THIS VERSION */ /* aprime and bprime "bracket" the max. We see if we can shrink aprime toward bprime (and/or the opposite) by evaluating two "middle" values located at Fibonacci #'s between aprime and bprime */ while ( (bprime - aprime) > 1.5*(b-a)/n ) { diff = evaluate(fib[k-2],n,offset) - evaluate(fib[k-1],n,offset); printf("Checked y[%d] = %f vs. y[%d] = %f; \n", fib[k-2]+offset, evaluate(fib[k-2],n,offset), fib[k-1]+offset, evaluate(fib[k-1],n,offset)); if (diff<0.0) { aprime = X(fib[k-2],n,offset); offset = offset + fib[k-2]; k = k - 1; } else if (diff>0.0) { bprime = X(fib[k-1],n,offset); /* offset doesn't move */ k = k - 1; } else { aprime = X(fib[k-2],n,offset); bprime = X(fib[k-1],n,offset); offset = offset + fib[k-2]; k = k - 3; } printf("Resetting aprime = %f, bprime = %f\n\n", aprime,bprime); /* for (i=0; i