TPCCLIB
Loading...
Searching...
No Matches
Macros | Functions
tgo.c File Reference

Topographical global optimization algorithm. More...

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include "tpcextensions.h"
#include "tpcrand.h"
#include "tpcstatist.h"
#include "tpcnlopt.h"

Go to the source code of this file.

Macros

#define MAX_PARAMETERS   50

Functions

int nloptITGO1 (NLOPT *nlo, const int doLocal, unsigned int tgoNr, unsigned int sampleNr, unsigned int neighNr, TPCSTATUS *status)
int nloptITGO2 (NLOPT *nlo, const int doLocal, unsigned int tgoNr, unsigned int sampleNr, unsigned int neighNr, TPCSTATUS *status)
int nloptIATGO (NLOPT *nlo, const int doLocal, unsigned int maxIterNr, double neighFract, TPCSTATUS *status)

Detailed Description

Topographical global optimization algorithm.

Author
Vesa Oikonen
Copyright
(c) Turku PET Centre

Based on an algorithm by Aimo Törn and Sami Viitanen. Reference: Törn A, Viitanen S. Iterative topographical global optimization. In: C.A. Floudas and P.M. Pardalos (eds.) State of the Art in Global Optimization. Kluwer Academic Publishers, 1996, pp 353-363. https://doi.org/10.1007/978-1-4613-3437-8_22

Remarks
Not well tested!

Definition in file tgo.c.

Macro Definition Documentation

◆ MAX_PARAMETERS

#define MAX_PARAMETERS   50

Max nr of parameters

Definition at line 32 of file tgo.c.

Referenced by nloptITGO1(), and nloptITGO2().

Function Documentation

◆ nloptIATGO()

int nloptIATGO ( NLOPT * nlo,
const int doLocal,
unsigned int maxIterNr,
double neighFract,
TPCSTATUS * status )

Iterative accumulative topographical global optimization algorithm (IATGO).

Precondition
Uses drand(), therefore set seed for a new series of pseudo-random numbers with drandSeed(1) before calling this function.
Postcondition
The last objective function call is usually not done with the best parameter estimates; if objective function simulates data that you need, you must call the function with the final parameters.
Returns
enum tpcerror (TPCERROR_OK when successful).
Author
Vesa Oikonen
See also
nloptMPSO, nloptSimplexARRS
Parameters
nloPointer to NLOPT data. Counter funCalls is initially set to zero, and then increased here.
Precondition
Constraints xlower[] and xupper[] are required. Parameter maxFunCalls is used as one of the stopping criteria. Parameters xtol[] is used for scaling and as one of the stopping criteria. Initial guess can be given in xfull[], but it is not obligatory. Initial step sizes (xdelta[]) are not used.
Parameters
doLocalLocal optimization method: 0=no local optimization, 1=Nelder-Mead.
maxIterNrNumber of TGO iterations; enter 0 to use the default; enter 1 to run TGO just once, ie to use TGO instead of iTGO.
neighFractFraction of samples to define as neighbours (cluster size); enter a large fraction (>0.5) if only few distant local minima are expected.
statusPointer to status data; enter NULL if not needed.

Definition at line 681 of file tgo.c.

701 {
702 FILE *fp=stdout;
703 int verbose=0; if(status!=NULL) {verbose=status->verbose; fp=status->fp;}
704 if(verbose>1) {
705 fprintf(fp, "%s(NLOPT, %d, %u, %g, status)\n", __func__, doLocal, maxIterNr, neighFract);
706 fflush(fp);
707 }
708 if(nlo==NULL ) {
709 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
710 return TPCERROR_FAIL;
711 }
712 if(nlo->totalNr<1 || nlo->xfull==NULL || nlo->_fun==NULL) {
713 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
714 return TPCERROR_NO_DATA;
715 }
716
717 unsigned int dim=nlo->totalNr;
718
719 if(verbose>4) {
720 fprintf(fp, "\nInitial limits and tolerances:\n");
721 for(unsigned int i=0; i<dim; i++)
722 fprintf(fp, "\t%g\t%g\t%e\n", nlo->xlower[i], nlo->xupper[i], nlo->xtol[i]);
723 }
724
725 /* Check if any of the parameters are fixed */
726 unsigned int fixedParNr=nloptLimitFixedNr(nlo);
727 if(verbose>2 && fixedParNr>0) fprintf(fp, "fixedParNr := %d\n", fixedParNr);
728 unsigned int fittedParNr=nlo->totalNr-fixedParNr;
729 if(verbose>2) {fprintf(fp, "fittedParNr := %d\n", fittedParNr); fflush(fp);}
730 if(fittedParNr<1) {
731 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
732 return TPCERROR_NO_DATA;
733 }
734
735 /* Check the tolerations */
736 for(unsigned int i=0; i<dim; i++) {
737 if(nlo->xlower[i]>=nlo->xupper[i]) {nlo->xtol[i]=0.0; continue;}
738 if(!(nlo->xtol[i]>0.0)) {
739 if(verbose>0) {fprintf(stderr, "Error: invalid xtol[].\n"); fflush(stderr);}
740 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
741 return TPCERROR_NO_DATA;
742 }
743 }
744
745 /* Continue input checking */
746 if(neighFract<0.04) neighFract=0.04; else if(neighFract>0.97) neighFract=0.97;
747 unsigned int sampleNr=50*fittedParNr/neighFract;
748 if(maxIterNr==0) maxIterNr=5+2*fittedParNr;
749 unsigned int neighNr=neighFract*sampleNr;
750 if(verbose>2) {
751 fprintf(fp, "sampleNr := %u\n", sampleNr);
752 fprintf(fp, "neighFract := %g\n", neighFract);
753 fprintf(fp, "neighNr := %u\n", neighNr);
754 fprintf(fp, "maxIterNr := %u\n", maxIterNr);
755 fflush(stdout);
756 }
757
758
759 /*
760 * Initialize the sample list with random points
761 */
762 if(verbose>3) {fprintf(fp, "Making initial random samples\n"); fflush(fp);}
763
764 /* Allocate memory */
765 nlo->usePList=1; // All points are stored
766 nlo->funCalls=0;
767 nlo->plist=(double*)malloc(sampleNr*(dim+1)*sizeof(double));
768 if(nlo->plist==NULL) { // will be freed with nloptFree()
769 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
770 return TPCERROR_OUT_OF_MEMORY;
771 }
772
773 /* Fill the list with random points and compute the object function values */
774 {
775 unsigned int badNr=0;
776 /* Add the random points */
777 for(unsigned int si=0; si<sampleNr; si++) {
778 if(nloptRandomPoint(&nlo->plist[si*(dim+1)], nlo->xlower, nlo->xupper, dim, NULL)) {
779 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
780 return TPCERROR_NO_DATA;
781 }
782 nlo->plist[si*(dim+1)+dim]=(*nlo->_fun)(nlo->totalNr, &nlo->plist[si*(dim+1)], nlo->fundata);
783 nlo->funCalls++;
784 if(!isfinite(nlo->plist[si*(dim+1)+dim])) badNr++;
785 }
786 if(verbose>3 && badNr>0) fprintf(fp, "Number of bad points: %d\n", badNr);
787 /* Not all object function values can be bad */
788 if((sampleNr-badNr)<10) {
789 if(verbose>0) {
790 fprintf(stderr, "Error: invalid values inside parameter range.\n"); fflush(stderr);}
791 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_BAD_FIT);
792 return TPCERROR_BAD_FIT;
793 }
794 /* Try new points to replace the failed ones */
795 if(badNr>0) {
796 badNr=0;
797 for(unsigned int si=0; si<sampleNr; si++) if(!isfinite(nlo->plist[si*(dim+1)+dim])) {
798 nloptRandomPoint(&nlo->plist[si*(dim+1)], nlo->xlower, nlo->xupper, dim, NULL);
799 nlo->plist[si*(dim+1)+dim]=(*nlo->_fun)(nlo->totalNr, &nlo->plist[si*(dim+1)], nlo->fundata);
800 //nlo->funCalls++;
801 if(!isfinite(nlo->plist[si*(dim+1)+dim])) badNr++;
802 }
803 if(verbose>4 && badNr>0) fprintf(fp, "Number of bad points after retry: %d\n", badNr);
804 }
805 }
806
807 /* If initial guess is valid, add it to the list */
808 int initGuessAvailable=1;
809 double initGuessCost=nan("");
810 for(unsigned int i=0; i<dim; i++) {
811 if(isnan(nlo->xfull[i])) {initGuessAvailable=0; break;}
812 if(nlo->xfull[i]<nlo->xlower[i]) {initGuessAvailable=0; break;}
813 if(nlo->xfull[i]>nlo->xupper[i]) {initGuessAvailable=0; break;}
814 }
815 if(initGuessAvailable) {
816 initGuessCost=nlo->funval=(*nlo->_fun)(nlo->totalNr, nlo->xfull, nlo->fundata);
817 if(isfinite(initGuessCost)) {
818 nlo->funCalls++;
819 nloptAddP(nlo, nlo->xfull, initGuessCost);
820 } else {
821 initGuessAvailable=0;
822 }
823 }
824 if(verbose>2) {
825 if(initGuessAvailable) fprintf(fp, "valid initial guess with cost=%g\n", initGuessCost);
826 else fprintf(fp, "no valid initial guess provided.\n");
827 }
828
829
830 /*
831 * Start iterations
832 */
833 if(verbose>2) {fprintf(fp, "\nStarting TGO iterations\n"); fflush(fp);}
834 unsigned int iterNr=0, stopCounter=0;
835 double prevBest=nan("");
836 while(iterNr<maxIterNr && nlo->funCalls<nlo->maxFunCalls) {
837 iterNr++;
838 unsigned int sampleNr=nlo->funCalls; // LOCAL sampleNr, increasing in each iteration
839 if(verbose>3) {fprintf(fp, "\nIteration %d with %d samples\n", iterNr, sampleNr); fflush(fp);}
840
841 /* Sort samples based on the evaluated function value */
842 if(nloptSortP(nlo)!=TPCERROR_OK) {
843 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL); return TPCERROR_FAIL;}
844
845 /* Print sampled points */
846 if(verbose>15) {fprintf(fp, "Points so far:\n"); nloptPrintP(nlo, 0, fp);}
847 else if(verbose>6) {fprintf(fp, "Best points so far:\n"); nloptPrintP(nlo, 10, fp);}
848 else if(verbose>4) {
849 fprintf(fp, "best point so far:");
850 for(unsigned int i=0; i<dim; i++) fprintf(fp, " %e", nlo->plist[i]);
851 fprintf(fp, " => %e\n", nlo->plist[dim]); fflush(fp);
852 }
853
854 /* If SD of the best points is minimal, and fit is not improving, then stop */
855 if(verbose>4 && !isnan(prevBest))
856 fprintf(fp, "dif_to_prev_point := %e\n", prevBest-nlo->plist[dim]);
857 if(!isnan(prevBest) && (prevBest-nlo->plist[dim])<1.0E-20) {
858 if(verbose>4) fprintf(fp, "Best point did not improve.\n");
859 if(nloptMeanP(nlo, neighNr, nlo->xfull, nlo->xdelta)==TPCERROR_OK) {
860 unsigned int i=0;
861 for(i=0; i<dim; i++) if(fabs(nlo->xdelta[i])>nlo->xtol[i]) break;
862 if(i==dim) {
863 if(verbose>2) fprintf(fp, "Small SD of best points.\n");
864 stopCounter++;
865 } else stopCounter=0;
866 }
867 } else stopCounter=0;
868 prevBest=nlo->plist[dim];
869
870 /* Stop, if stopping rules bumped into several times */
871 if(verbose>4) fprintf(fp, "stopCounter := %d\n", stopCounter);
872 if(stopCounter>2) {
873 if(verbose>2) fprintf(fp, "\n Required tolerance reached.\n");
874 break;
875 }
876
877 /* Find topographic minima (TM) */
878
879 /* Allocate a list specifying whether point is local minimum (1) or not (0) */
880 /* Allocate a list specifying whether point is part of topographic minimum (>0) or not (0) */
881 unsigned int tmlist[sampleNr];
882 for(unsigned int i=0; i<sampleNr; i++) tmlist[i]=0;
883
884 unsigned int topoNr=0;
885 if(verbose>3) {fprintf(fp, "neighNr := %d\n", neighNr); fflush(fp);}
886
887 /* For each point i, find out if it is a "topographic minimum", that is,
888 better than k neighbour points. */
889 for(unsigned int si=0; si<sampleNr-neighNr; si++) if(tmlist[si]==0) {
890
891 /* If function value is invalid, then this cannot be accepted as TM */
892 if(!isfinite(nlo->plist[si*(dim+1)+dim])) continue;
893
894 /* Compute the distances to other points, scaled using xtol[] */
895 double sdist[sampleNr];
896 for(unsigned int sj=0; sj<sampleNr; sj++) {
897 if(si==sj || !isfinite(nlo->plist[sj*(dim+1)+dim])) { // point itself not included and
898 sdist[sj]=nan(""); continue;} // valid function value required
899 sdist[sj]=0.0;
900 for(unsigned int k=0; k<dim; k++) {
901 if(!(nlo->xtol[k]>0.0)) continue;
902 double d=(nlo->plist[si*(dim+1)+k]-nlo->plist[sj*(dim+1)+k])/nlo->xtol[k];
903 /*if(isfinite(d))*/ sdist[sj]+=d*d;
904 }
905 /* If scaled distance is less than 1, the points are practically equal: leave out */
906 //if(sdist[sj]<1.0) sdist[sj]=nan(""); // would lead to multiple local minima at same point
907 }
908 double sdist2[sampleNr]; // make copy for printing
909 for(unsigned int sj=0; sj<sampleNr; sj++) sdist2[sj]=sdist[sj];
910
911 /* Find its closest neighbours that are not any better */
912 unsigned int nn=0; double prev_mind=nan("");
913 while(1) {
914 /* Find the closest neighbour that was not already found */
915 unsigned int mini=0;
916 double mind=nan("");
917 for(unsigned int sj=0; sj<sampleNr; sj++) {
918 if(!isfinite(mind) || mind>sdist[sj]) {mind=sdist[sj]; mini=sj;}
919 }
920 if(!isfinite(mind)) break;
921 sdist[mini]=nan(""); // this point not used again in search of neighbours
922 if(verbose>100) fprintf(fp, " min_distance[%u] := %g , with fval := %g\n",
923 1+nn, mind, nlo->plist[mini*(dim+1)+dim]);
924 /* If this neighbour is better than the tentative TM, then stop */
925 if(nlo->plist[mini*(dim+1)+dim] < nlo->plist[si*(dim+1)+dim]) break;
926 /* Otherwise, count this as a worse neighbour, belonging to this TM */
927 nn++; tmlist[mini]=1+si;
928 // we do not want to include more than neighNr points into this TM,
929 // but if distance is less than one or not larger than previously then
930 // this is essentially the same point
931 if(isnan(prev_mind)) prev_mind=mind;
932 if(nn>=neighNr && mind>1.0 && mind>prev_mind) break;
933 prev_mind=mind;
934 }
935 /* If the number of worse neighbours is less than required, then stop considering this
936 as a local optimum, and continue with the next sample */
937 if(nn<neighNr) {
938 /* its neighbours then do not belong to TM either */
939 for(unsigned int ni=0; ni<sampleNr; ni++) if(tmlist[ni]==1+si) tmlist[ni]=0;
940 continue;
941 }
942 /* otherwise mark this as topographic minimum (TM) */
943 tmlist[si]=1+si;
944 topoNr++;
945 /* Print TM point */
946 if(verbose>5) {
947 fprintf(fp, "TM point %d\t", 1+si);
948 for(unsigned int i=0; i<dim; i++) fprintf(fp, "%e ", nlo->plist[si*(dim+1)+i]);
949 fprintf(fp, "=> %e\n", nlo->plist[si*(dim+1)+dim]);
950 if(verbose>14) {
951 fprintf(fp, "and its %u neighbours:\n", nn);
952 for(unsigned int sj=0; sj<sampleNr; sj++) if(sj!=si && tmlist[sj]==1+si) {
953 for(unsigned int i=0; i<dim; i++) fprintf(fp, "%e ", nlo->plist[sj*(dim+1)+i]);
954 fprintf(fp, "=> %e\t(%g)\n", nlo->plist[sj*(dim+1)+dim], sdist2[sj]);
955 }
956 }
957 fflush(fp);
958 }
959 }
960 if(verbose>3) {fprintf(fp, "topographic minima: %d\n", topoNr); fflush(fp);}
961 if(topoNr==0) {
962 if(verbose>0) {fprintf(stderr, "Error: no topographical minima found\n"); fflush(stderr);}
963 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_BAD_FIT);
964 return TPCERROR_BAD_FIT;
965 }
966
967 /*
968 * Local optimization from each TM
969 */
970 for(unsigned int si=0; si<sampleNr; si++) if(tmlist[si]==1+si) {
971 if(verbose>5) {fprintf(fp, "LO for TM point %d\n", 1+si); fflush(fp);}
972 if(verbose>6) {
973 fprintf(fp, "TM point %d before LO:\n", 1+si);
974 for(unsigned int i=0; i<dim; i++) fprintf(fp, " %e", nlo->plist[si*(dim+1)+i]);
975 fprintf(fp, " => %e\n", nlo->plist[si*(dim+1)+dim]); fflush(fp);
976 }
977 /* Compute the mean and SD of points belonging to this TM */
978 double meanp[dim], sdp[dim]; unsigned int nn=0;
979 for(unsigned int i=0; i<dim; i++) {
980 if(nlo->xupper[i]>nlo->xlower[i]) {
981 double x[sampleNr];
982 nn=0;
983 for(unsigned int sj=0; sj<sampleNr; sj++)
984 if(tmlist[sj]==1+si) x[nn++]=nlo->plist[sj*(dim+1)+i];
985 statMeanSD(x, nn, meanp+i, sdp+i, NULL);
986 } else {
987 meanp[i]=nlo->xlower[i]; sdp[i]=0.0;
988 }
989 }
990 if(verbose>7) {
991 fprintf(fp, "TM point (n=%u) mean and SD:\n", nn);
992 for(unsigned int i=0; i<dim; i++) fprintf(fp, "\t%g\t%g\n", meanp[i], sdp[i]);
993 fflush(fp);
994 }
995 /* Make sure that SD is not zero, unless parameter is fixed */
996 for(unsigned int i=0; i<dim; i++) if(nlo->xupper[i]>nlo->xlower[i]) {
997 if(sdp[i]<0.1*nlo->xtol[i]) sdp[i]=drandExponential(0.05*nlo->xtol[i]);
998 }
999 if(doLocal==1) { // downhill simplex from the best point so far
1000 doubleCopy(nlo->xfull, &nlo->plist[si*(dim+1)], dim);
1001 doubleCopy(nlo->xdelta, sdp, dim);
1002 int ret=nloptSimplex(nlo, 100*dim, NULL);
1003 if(ret!=TPCERROR_OK) {
1004 if(verbose>1) {fprintf(fp, " LO failed: %s\n", errorMsg(ret)); fflush(fp);}
1005 } else {
1006 if(verbose>6) {
1007 fprintf(fp, "Point after LO:");
1008 for(unsigned int i=0; i<dim; i++) fprintf(fp, " %e ", nlo->xfull[i]);
1009 fprintf(fp, " => %e\n", nlo->funval); fflush(fp);
1010 }
1011 }
1012 } else { // Add random points with Gaussian distribution
1013 unsigned int bi=0; double bval=nan("");
1014 for(unsigned int ni=0; ni<neighNr; ni++) {
1015 double p[dim];
1016 nloptGaussianPoint(p, &nlo->plist[si*(dim+1)], sdp, nlo->xlower, nlo->xupper, dim, NULL);
1017 double fval=(*nlo->_fun)(dim, p, nlo->fundata);
1018 if(isnan(bval) || bval>fval) {bval=fval; bi=nlo->funCalls;}
1019 if(isfinite(fval)) {nlo->funCalls++; nloptAddP(nlo, p, fval);}
1020 }
1021 if(verbose>6) {
1022 fprintf(fp, "Best TM point %d after LO:\n", 1+si);
1023 for(unsigned int i=0; i<dim; i++) fprintf(fp, " %e", nlo->plist[bi*(dim+1)+i]);
1024 fprintf(fp, " => %e\n", nlo->plist[bi*(dim+1)+dim]); fflush(fp);
1025 }
1026 }
1027 } // next TM
1028
1029 } // end of iterations
1030
1031 /* Check the reason for loop exit */
1032 if(verbose>1) {
1033 if(iterNr>=maxIterNr)
1034 fprintf(fp, "\n Exceeded the maximum number for loops.\n");
1035 if(nlo->funCalls>=nlo->maxFunCalls)
1036 fprintf(fp, "\n Exceeded the maximum number for function calls.\n");
1037 fflush(fp);
1038 }
1039
1040 /* Sort samples based on the evaluated function value */
1041 if(nloptSortP(nlo)!=TPCERROR_OK) {
1042 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL); return TPCERROR_FAIL;}
1043 /* Get the best point so far */
1044 for(unsigned int i=0; i<dim; i++) nlo->xfull[i]=nlo->plist[i];
1045 nlo->funval=nlo->plist[dim];
1046
1047
1048 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
1049 return TPCERROR_OK;
1050}
doubleCopy
void doubleCopy(double *t, double *s, const unsigned int n)
Definition doubleutil.c:117
drandExponential
double drandExponential(double mean)
Get pseudo-random number with exponential distribution.
Definition gaussdev.c:178
statMeanSD
int statMeanSD(double *data, unsigned int n, double *mean, double *sd, unsigned int *vn)
Definition mean.c:25
nloptLimitFixedNr
unsigned int nloptLimitFixedNr(NLOPT *d)
Definition nlopt.c:333
nloptMeanP
int nloptMeanP(NLOPT *nlo, unsigned int nr, double *meanp, double *sdp)
Definition nlopt.c:219
nloptPrintP
void nloptPrintP(NLOPT *nlo, unsigned int nr, FILE *fp)
Definition nlopt.c:274
nloptAddP
int nloptAddP(NLOPT *nlo, double *p, double funval)
Definition nlopt.c:149
nloptSortP
int nloptSortP(NLOPT *nlo)
Definition nlopt.c:182
nloptRandomPoint
int nloptRandomPoint(double *p, double *low, double *up, unsigned int n, MERTWI *mt)
Definition rndpoint.c:28
nloptGaussianPoint
int nloptGaussianPoint(double *p, double *mean, double *sd, double *low, double *up, unsigned int n, MERTWI *mt)
Definition rndpoint.c:70
nloptSimplex
int nloptSimplex(NLOPT *nlo, unsigned int maxIter, TPCSTATUS *status)
Definition simplex.c:32
errorMsg
char * errorMsg(tpcerror e)
Definition statusmsg.c:68
statusSet
void statusSet(TPCSTATUS *s, const char *func, const char *srcfile, int srcline, tpcerror error)
Definition statusmsg.c:142
NLOPT::_fun
double(* _fun)(int, double *, void *)
Definition tpcnlopt.h:42
NLOPT::funval
double funval
Definition tpcnlopt.h:50
NLOPT::maxFunCalls
unsigned int maxFunCalls
Definition tpcnlopt.h:46
NLOPT::xupper
double * xupper
Definition tpcnlopt.h:33
NLOPT::xlower
double * xlower
Definition tpcnlopt.h:31
NLOPT::fundata
void * fundata
Definition tpcnlopt.h:44
NLOPT::xfull
double * xfull
Definition tpcnlopt.h:29
NLOPT::funCalls
unsigned int funCalls
Definition tpcnlopt.h:48
NLOPT::xdelta
double * xdelta
Definition tpcnlopt.h:36
NLOPT::usePList
unsigned int usePList
Definition tpcnlopt.h:52
NLOPT::xtol
double * xtol
Definition tpcnlopt.h:39
NLOPT::plist
double * plist
Definition tpcnlopt.h:56
NLOPT::totalNr
unsigned int totalNr
Definition tpcnlopt.h:27
TPCSTATUS::verbose
int verbose
Verbose level, used by statusPrint() etc.
Definition tpcextensions.h:242
TPCSTATUS::fp
FILE * fp
File pointer for writing log information during development and testing.
Definition tpcextensions.h:243
TPCERROR_BAD_FIT
@ TPCERROR_BAD_FIT
Fitting not successful.
Definition tpcextensions.h:228
TPCERROR_FAIL
@ TPCERROR_FAIL
General error.
Definition tpcextensions.h:192
TPCERROR_OUT_OF_MEMORY
@ TPCERROR_OUT_OF_MEMORY
Cannot allocate memory.
Definition tpcextensions.h:193
TPCERROR_OK
@ TPCERROR_OK
No error.
Definition tpcextensions.h:191
TPCERROR_NO_DATA
@ TPCERROR_NO_DATA
File contains no data.
Definition tpcextensions.h:203

◆ nloptITGO1()

int nloptITGO1 ( NLOPT * nlo,
const int doLocal,
unsigned int tgoNr,
unsigned int sampleNr,
unsigned int neighNr,
TPCSTATUS * status )

Iterative Topographical global optimization algorithm 1 (iTGO1).

Remarks
Not well tested!
Precondition
Uses drand(), therefore set seed for a new series of pseudo-random numbers with drandSeed(1) before calling this function.
Postcondition
The last objective function call is usually not done with the best parameter estimates; if objective function simulates data that you need, you must call the function with the final parameters.
Returns
enum tpcerror (TPCERROR_OK when successful).
Author
Vesa Oikonen
Todo
Finalize the code and tests.
See also
nloptIATGO, nloptMPSO, nloptSimplexARRS, nloptITGO2
Parameters
nloPointer to NLOPT data. Counter funCalls is initially set to zero, and then increased here. Parameter maxFunCalls is used as one of the stopping criteria. Parameters xtol[] is used for scaling and as one of the stopping criteria.
doLocalLocal optimization method: 0=no local optimization, 1=Nelder-Mead.
tgoNrNumber of TGO iterations; enter 0 to use the default; enter 1 to run TGO just once, ie to use TGO instead of iTGO. Iterations are needed if sampleNr (below) would otherwise require too much memory.
sampleNrNumber of points to sample in one iteration; may need to be large if the number of iterations (above) is small; enter 0 to let function decide it.
neighNrNumber of neighbours to investigate (cluster size); enter a large number relative to sampleNr (above) if only few local minima are expected; enter 0 to let function decide it.
statusPointer to status data; enter NULL if not needed.

Definition at line 59 of file tgo.c.

81 {
82 int verbose=0; if(status!=NULL) verbose=status->verbose;
83 if(verbose>0) {
84 printf("%s(NLOPT, %d, %u, %u, %u, status)\n", __func__, doLocal, tgoNr, sampleNr, neighNr);
85 fflush(stdout);
86 }
87 if(nlo==NULL ) {
88 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
89 return TPCERROR_FAIL;
90 }
91 if(nlo->totalNr<1 || nlo->xfull==NULL || nlo->_fun==NULL) {
92 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
93 return TPCERROR_NO_DATA;
94 }
95 if(nlo->totalNr>MAX_PARAMETERS) {
96 if(verbose>0) {fprintf(stderr, "Error: too many dimensions to optimize.\n"); fflush(stderr);}
97 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
98 return TPCERROR_FAIL;
99 }
100
101 /* Check if any of the parameters are fixed */
102 unsigned int fixedParNr=nloptLimitFixedNr(nlo);
103 if(verbose>1 && fixedParNr>0) printf("fixedParNr := %d\n", fixedParNr);
104 unsigned int fittedParNr=nlo->totalNr-fixedParNr;
105 if(verbose>1) printf("fittedParNr := %d\n", fittedParNr);
106 if(fittedParNr<1) {
107 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
108 return TPCERROR_NO_DATA;
109 }
110
111 /* Check the tolerations */
112 for(unsigned int i=0; i<nlo->totalNr; i++) {
113 if(nlo->xlower[i]>=nlo->xupper[i]) {nlo->xtol[i]=0.0; continue;}
114 if(!(nlo->xtol[i]>0.0)) {
115 if(verbose>0) {fprintf(stderr, "Error: invalid xtol[].\n"); fflush(stderr);}
116 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
117 return TPCERROR_NO_DATA;
118 }
119 }
120
121 /* Continue input checking */
122 if(sampleNr<10) sampleNr=50*fittedParNr;
123 if(sampleNr&1) sampleNr++; // If number is odd, then add 1
124 if(neighNr<2) neighNr=sampleNr/2;
125 else if(neighNr>sampleNr-1) neighNr=sampleNr-1; // Make sure "neighNr" isn't too big
126 if(tgoNr==0) tgoNr=1+fittedParNr;
127 if(verbose>2) {
128 printf("sampleNr := %d\n", sampleNr);
129 printf("neighNr := %d\n", neighNr);
130 printf("tgoNr := %d\n", tgoNr);
131 fflush(stdout);
132 }
133
134 /* Check if initial guess is valid */
135 int initGuessAvailable=1;
136 double initGuessCost=nan("");
137 for(unsigned int i=0; i<nlo->totalNr; i++) {
138 if(isnan(nlo->xfull[i])) {initGuessAvailable=0; break;}
139 if(nlo->xfull[i]<nlo->xlower[i]) {initGuessAvailable=0; break;}
140 if(nlo->xfull[i]>nlo->xupper[i]) {initGuessAvailable=0; break;}
141 }
142 if(initGuessAvailable) {
143 initGuessCost=nlo->funval=(*nlo->_fun)(nlo->totalNr, nlo->xfull, nlo->fundata);
144 nlo->funCalls++;
145 if(!isfinite(initGuessCost)) initGuessAvailable=0;
146 else if(verbose>2) {
147 printf("valid initial guess with cost=%g\n", initGuessCost);
148 }
149 }
150
151 /* Allocate memory */
152 TGO_POINT *points=(TGO_POINT*)malloc(sampleNr*sizeof(TGO_POINT));
153 if(points==NULL) {
154 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
155 return TPCERROR_OUT_OF_MEMORY;
156 }
157 for(unsigned int i=0; i<sampleNr; i++) {
158 points[i].topomin=0;
159 points[i].fvalue=0.0;
160 }
161
162 /*
163 * Iterative TGO, or non-iterative if tgoNr==1
164 */
165 unsigned int topoNr=0;
166 for(unsigned int tgoi=0; tgoi<tgoNr; tgoi++) {
167
168 if(verbose>3 && tgoNr>1) {printf("TGO iteration # %d: \n", 1+tgoi); fflush(stdout);}
169
170 /* Sample N points in the feasible region and compute the object function values for
171 those points which do not already have it. */
172 unsigned int badNr=0;
173 for(unsigned int si=0; si<sampleNr; si++) if(points[si].topomin==0) {
174 if(tgoi==0 && si==0 && initGuessAvailable) {
175 /* If initial guess was given, then use it as the first point */
176 for(unsigned int i=0; i<nlo->totalNr; i++) points[si].par[i]=nlo->xfull[i];
177 points[si].fvalue=initGuessCost;
178 continue;
179 }
180 nloptRandomPoint(points[si].par, nlo->xlower, nlo->xupper, nlo->totalNr, NULL);
181 points[si].fvalue=(*nlo->_fun)(nlo->totalNr, points[si].par, nlo->fundata);
182 nlo->funCalls++;
183 /* If function return value was not valid, then we'll try twice with new guesses */
184 int tries=0;
185 while(!isfinite(points[si].fvalue) && tries<2) {
186 if(verbose>8) {
187 printf("this point did not give normal return value:\n");
188 for(unsigned int i=0; i<nlo->totalNr; i++) printf(" %10.2e", points[si].par[i]);
189 printf("\n");
190 }
191 nloptRandomPoint(points[si].par, nlo->xlower, nlo->xupper, nlo->totalNr, NULL);
192 points[si].fvalue=(*nlo->_fun)(nlo->totalNr, points[si].par, nlo->fundata);
193 nlo->funCalls++;
194 tries++;
195 }
196 if(!isfinite(points[si].fvalue)) badNr++;
197 }
198 if(verbose>4 && badNr>0) printf("Number of bad points: %d\n", badNr);
199 /* Object functions values must be good for at least NeighNr points */
200 if(tgoi==0 && (sampleNr-badNr)<=neighNr) {
201 if(verbose>0) {
202 fprintf(stderr, "Error: invalid values inside parameter range.\n"); fflush(stderr);}
203 free(points);
204 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_BAD_FIT);
205 return TPCERROR_BAD_FIT;
206 }
207 /* Print sampled points */
208 if(verbose>10) {
209 printf("Sampled points:\n");
210 for(unsigned int si=0; si<sampleNr; si++) {
211 printf("%d\t", 1+si);
212 for(unsigned int i=0; i<nlo->totalNr; i++) printf("%e ", points[si].par[i]);
213 printf("=> %e\n", points[si].fvalue);
214 }
215 fflush(stdout);
216 }
217
218 /* For each point i, find out if it is a "topographic minimum", that is,
219 better than k neighbour points. */
220 topoNr=0;
221 for(unsigned int si=0; si<sampleNr; si++) {
222 points[si].topomin=0;
223
224 /* Compute the distances, scaled using xtol[] */
225 double sdist[sampleNr];
226 for(unsigned int sj=0; sj<sampleNr; sj++) {
227 sdist[sj]=0.0;
228 if(si==sj) {sdist[sj]=1.0E+99; continue;}
229 for(unsigned int k=0; k<nlo->totalNr; k++) {
230 if(!(nlo->xtol[k]>0.0)) continue;
231 double d=(points[si].par[k]-points[sj].par[k])/nlo->xtol[k];
232 if(isfinite(d)) sdist[sj]+=d*d;
233 }
234 /* Distance should be computed as square root, but it does not affect
235 the order of distances; since sqrt() is slow, it is not done. */
236 // sdist[sj]=sqrt(sdist[sj]);
237 }
238
239 /* Find the closest neighbours for sample i */
240 unsigned int neighs[neighNr]; // collect here the neighbours
241 unsigned int ni=0;
242 for(ni=0; ni<neighNr; ni++) {
243 unsigned int mini=0;
244 double minv=sdist[mini];
245 for(unsigned int sj=0; sj<sampleNr; sj++) {
246 if(sdist[sj]<minv) {minv=sdist[sj]; mini=sj;}
247 }
248 sdist[mini]=1.0E+99; // not used again
249 /* If point i is worse than any of the closest neighbours, then point i
250 is certainly not a topographic minimum; then stop this loop and go to
251 the next point i+1 */
252 if(!(points[si].fvalue<points[mini].fvalue)) break;
253 neighs[ni]=mini;
254 }
255 if(ni<neighNr) continue; // continue with the next i
256 /* otherwise mark this as topographic minimum (TM) */
257 points[si].topomin=1; topoNr++;
258 /* Print TM point */
259 if(verbose>6) {
260 printf("TM point %d:\n", topoNr);
261 printf("%d\t", 1+si);
262 for(unsigned int i=0; i<nlo->totalNr; i++) printf("%e ", points[si].par[i]);
263 printf("=> %e\n", points[si].fvalue);
264 if(verbose>7) {
265 printf(" points neighbours:");
266 for(ni=0; ni<neighNr; ni++) printf(" %d", 1+neighs[ni]);
267 printf("\n");
268 }
269 fflush(stdout);
270 }
271
272 }
273 if(verbose>4) {printf(" %d topographical minima\n", topoNr); fflush(stdout);}
274 if(topoNr==0) {
275 if(verbose>0) {fprintf(stderr, "Warning: no topographical minima found\n"); fflush(stderr);}
276 continue;
277 }
278
279 } // end of iTGO
280
281
282 /* If requested, do local optimization for the TMs */
283 if(topoNr==0) {
284 if(verbose>0) {fprintf(stderr, "Error: no topographical minima found\n"); fflush(stderr);}
285 free(points);
286 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_BAD_FIT);
287 return TPCERROR_BAD_FIT;
288 }
289 for(unsigned int si=0; si<sampleNr; si++) if(points[si].topomin>0) {
290 if(verbose>2) printf("LO for TM point %d\n", 1+si);
291 NLOPT lnlo; nloptInit(&lnlo);
292 nloptDuplicate(nlo, &lnlo);
293 for(unsigned int i=0; i<nlo->totalNr; i++) lnlo.xfull[i]=points[si].par[i];
294 lnlo.funval=points[si].fvalue;
295 for(unsigned int i=0; i<nlo->totalNr; i++) lnlo.xdelta[i]=100.0*lnlo.xtol[i];
296 lnlo.maxFunCalls=2000*nlo->totalNr; lnlo.funCalls=0;
297 int ret=nloptSimplex(&lnlo, 100*nlo->totalNr, status);
298 if(ret==TPCERROR_OK) {
299 for(unsigned int i=0; i<nlo->totalNr; i++) points[si].par[i]=lnlo.xfull[i];
300 points[si].fvalue=lnlo.funval;
301 if(verbose>6) {
302 printf("TM point %d after LO:\n", 1+si);
303 printf("%d\t", 1+si);
304 for(unsigned int i=0; i<nlo->totalNr; i++) printf("%e ", points[si].par[i]);
305 printf("=> %e\n", points[si].fvalue); fflush(stdout);
306 }
307 } else if(verbose>7) {printf(" LO failed\n"); fflush(stdout);}
308 nlo->funCalls+=lnlo.funCalls;
309 nloptFree(&lnlo);
310 }
311
312 /* Find the best minimum and return it */
313 {
314 unsigned int mini=0;
315 double minv=points[mini].fvalue;
316 for(unsigned int si=1; si<sampleNr; si++)
317 if(!(points[si].fvalue>minv)) {minv=points[si].fvalue; mini=si;}
318 for(unsigned int i=0; i<nlo->totalNr; i++) nlo->xfull[i]=points[mini].par[i];
319 nlo->funval=points[mini].fvalue;
320 }
321
322 free(points);
323 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
324 return TPCERROR_OK;
325}
nloptInit
void nloptInit(NLOPT *nlo)
Definition nlopt.c:25
nloptFree
void nloptFree(NLOPT *nlo)
Definition nlopt.c:52
nloptDuplicate
int nloptDuplicate(NLOPT *nlo1, NLOPT *nlo2)
Definition nlopt.c:112
MAX_PARAMETERS
#define MAX_PARAMETERS
Definition tgo.c:32

◆ nloptITGO2()

int nloptITGO2 ( NLOPT * nlo,
const int doLocal,
unsigned int tgoNr,
unsigned int sampleNr,
unsigned int neighNr,
TPCSTATUS * status )

Iterative Topographical global optimization algorithm 2 (iTGO2).

Remarks
Not well tested!
Precondition
Uses drand(), therefore set seed for a new series of pseudo-random numbers with drandSeed(1) before calling this function.
Postcondition
The last objective function call is usually not done with the best parameter estimates; if objective function simulates data that you need, you must call the function with the final parameters.
Returns
enum tpcerror (TPCERROR_OK when successful).
Author
Vesa Oikonen
Todo
Finalize the code and tests.
See also
nloptIATGO, nloptMPSO, nloptSimplexARRS, nloptITGO1
Parameters
nloPointer to NLOPT data. Counter funCalls is initially set to zero, and then increased here. Parameter maxFunCalls is used as one of the stopping criteria. Parameters xtol[] is used for scaling and as one of the stopping criteria.
doLocalLocal optimization method: 0=no local optimization, 1=Nelder-Mead.
tgoNrNumber of TGO iterations; enter 0 to use the default; enter 1 to run TGO just once, ie to use TGO instead of iTGO. Iterations are needed if sampleNr (below) would otherwise require too much memory.
sampleNrNumber of points to sample in one iteration; may need to be large if the number of iterations (above) is small; enter 0 to let function decide it.
neighNrNumber of neighbours to investigate (cluster size); enter a large number relative to sampleNr (above) if only few local minima are expected; enter 0 to let function decide it.
statusPointer to status data; enter NULL if not needed.

Definition at line 342 of file tgo.c.

364 {
365 int verbose=0; if(status!=NULL) verbose=status->verbose;
366 if(verbose>0) {
367 printf("%s(NLOPT, %d, %u, %u, %u, status)\n", __func__, doLocal, tgoNr, sampleNr, neighNr);
368 fflush(stdout);
369 }
370 if(nlo==NULL ) {
371 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
372 return TPCERROR_FAIL;
373 }
374 if(nlo->totalNr<1 || nlo->xfull==NULL || nlo->_fun==NULL) {
375 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
376 return TPCERROR_NO_DATA;
377 }
378 if(nlo->totalNr>MAX_PARAMETERS) {
379 if(verbose>0) {fprintf(stderr, "Error: too many dimensions to optimize.\n"); fflush(stderr);}
380 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
381 return TPCERROR_FAIL;
382 }
383
384 /* Check if any of the parameters are fixed */
385 unsigned int fixedParNr=nloptLimitFixedNr(nlo);
386 if(verbose>1 && fixedParNr>0) printf("fixedParNr := %d\n", fixedParNr);
387 unsigned int fittedParNr=nlo->totalNr-fixedParNr;
388 if(verbose>1) printf("fittedParNr := %d\n", fittedParNr);
389 if(fittedParNr<1) {
390 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
391 return TPCERROR_NO_DATA;
392 }
393
394 /* Check the tolerations */
395 for(unsigned int i=0; i<nlo->totalNr; i++) {
396 if(nlo->xlower[i]>=nlo->xupper[i]) {nlo->xtol[i]=0.0; continue;}
397 if(!(nlo->xtol[i]>0.0)) {
398 if(verbose>0) {fprintf(stderr, "Error: invalid xtol[].\n"); fflush(stderr);}
399 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
400 return TPCERROR_NO_DATA;
401 }
402 }
403
404 /* Continue input checking */
405 if(sampleNr<10) sampleNr=50*fittedParNr;
406 if(sampleNr&1) sampleNr++; // If number is odd, then add 1
407 if(neighNr<2) neighNr=sampleNr/2;
408 else if(neighNr>sampleNr-1) neighNr=sampleNr-1; // Make sure "neighNr" isn't too big
409 if(tgoNr==0) tgoNr=1+fittedParNr;
410 if(verbose>2) {
411 printf("sampleNr := %d\n", sampleNr);
412 printf("neighNr := %d\n", neighNr);
413 printf("tgoNr := %d\n", tgoNr);
414 fflush(stdout);
415 }
416
417 /* Check if initial guess is valid */
418 int initGuessAvailable=1;
419 double initGuessCost=nan("");
420 for(unsigned int i=0; i<nlo->totalNr; i++) {
421 if(isnan(nlo->xfull[i])) {initGuessAvailable=0; break;}
422 if(nlo->xfull[i]<nlo->xlower[i]) {initGuessAvailable=0; break;}
423 if(nlo->xfull[i]>nlo->xupper[i]) {initGuessAvailable=0; break;}
424 }
425 if(initGuessAvailable) {
426 initGuessCost=nlo->funval=(*nlo->_fun)(nlo->totalNr, nlo->xfull, nlo->fundata);
427 nlo->funCalls++;
428 if(!isfinite(initGuessCost)) initGuessAvailable=0;
429 else if(verbose>2) {
430 printf("valid initial guess with cost=%g\n", initGuessCost);
431 }
432 }
433
434 /* Allocate memory */
435 TGO_POINT *points=(TGO_POINT*)malloc(sampleNr*sizeof(TGO_POINT));
436 if(points==NULL) {
437 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
438 return TPCERROR_OUT_OF_MEMORY;
439 }
440 for(unsigned int i=0; i<sampleNr; i++) {
441 points[i].topomin=0;
442 points[i].fvalue=0.0;
443 }
444 TGO_POINT *gpoints=(TGO_POINT*)malloc(sampleNr*sizeof(TGO_POINT));
445 if(gpoints==NULL) {
446 free(points);
447 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
448 return TPCERROR_OUT_OF_MEMORY;
449 }
450 for(unsigned int i=0; i<sampleNr; i++) {
451 gpoints[i].topomin=0;
452 gpoints[i].fvalue=0.0;
453 }
454
455 /*
456 * Iterative TGO, or non-iterative if tgoNr==1
457 */
458 unsigned int topoNr=0;
459 if(initGuessAvailable) {
460 /* If initial guess was given, then use it as the first G point */
461 for(unsigned int i=0; i<nlo->totalNr; i++) gpoints[0].par[i]=nlo->xfull[i];
462 gpoints[0].fvalue=initGuessCost;
463 topoNr=1;
464 }
465 for(unsigned int tgoi=0; tgoi<tgoNr; tgoi++) {
466
467 if(verbose>2 && tgoNr>1) {printf("TGO iteration # %d: \n", 1+tgoi); fflush(stdout);}
468
469 while(topoNr<sampleNr) {
470 if(verbose>3) {printf(" topoNr := %d\n", topoNr); fflush(stdout);}
471
472 /* Sample N points and compute the object function values for all points. */
473 unsigned int badNr=0;
474 for(unsigned int si=0; si<sampleNr; si++) {
475 points[si].topomin=0;
476 nloptRandomPoint(points[si].par, nlo->xlower, nlo->xupper, nlo->totalNr, NULL);
477 points[si].fvalue=(*nlo->_fun)(nlo->totalNr, points[si].par, nlo->fundata);
478 nlo->funCalls++;
479 /* If function return value was not valid, then we'll try twice with new guesses */
480 int tries=0;
481 while(!isfinite(points[si].fvalue) && tries<2) {
482 if(verbose>8) {
483 printf("this point did not give normal return value:\n");
484 for(unsigned int i=0; i<nlo->totalNr; i++) printf(" %10.2e", points[si].par[i]);
485 printf("\n");
486 }
487 nloptRandomPoint(points[si].par, nlo->xlower, nlo->xupper, nlo->totalNr, NULL);
488 points[si].fvalue=(*nlo->_fun)(nlo->totalNr, points[si].par, nlo->fundata);
489 nlo->funCalls++;
490 tries++;
491 }
492 if(!isfinite(points[si].fvalue)) badNr++;
493 }
494 if(verbose>4 && badNr>0) printf("Number of bad points: %d\n", badNr);
495 /* Object functions values must be good for at least NeighNr points */
496 if(tgoi==0 && (sampleNr-badNr)<=neighNr) {
497 if(verbose>0) {
498 fprintf(stderr, "Error: invalid values inside parameter range.\n"); fflush(stderr);}
499 free(points); free(gpoints);
500 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_BAD_FIT);
501 return TPCERROR_BAD_FIT;
502 }
503
504 /* For each point i, find out if it is a "topographic minimum", that is,
505 better than k neighbour points. */
506 for(unsigned int si=0; si<sampleNr; si++) {
507 points[si].topomin=0;
508
509 /* Compute the distances, scaled using xtol[] */
510 double sdist[sampleNr];
511 for(unsigned int sj=0; sj<sampleNr; sj++) {
512 sdist[sj]=0.0;
513 if(si==sj) {sdist[sj]=1.0E+99; continue;}
514 for(unsigned int k=0; k<nlo->totalNr; k++) {
515 if(!(nlo->xtol[k]>0.0)) continue;
516 double d=(points[si].par[k]-points[sj].par[k])/nlo->xtol[k];
517 if(isfinite(d)) sdist[sj]+=d*d;
518 }
519 /* Distance should be computed as square root, but it does not affect
520 the order of distances; since sqrt() is slow, it is not done. */
521 // sdist[sj]=sqrt(sdist[sj]);
522 }
523
524 /* Find the closest neighbours for sample i */
525 unsigned int ni=0;
526 for(ni=0; ni<neighNr; ni++) {
527 unsigned int mini=0;
528 double minv=sdist[mini];
529 for(unsigned int sj=0; sj<sampleNr; sj++) {
530 if(sdist[sj]<minv) {minv=sdist[sj]; mini=sj;}
531 }
532 sdist[mini]=1.0E+99; // not used again
533 /* If point i is worse than any of the closest neighbours, then point i is certainly not
534 a topographic minimum; then stop this loop and go to the next point i+1 */
535 if(!(points[si].fvalue<points[mini].fvalue)) break;
536 }
537 if(ni<neighNr) continue; // continue with the next sample i
538 /* otherwise mark this as topographic minimum (TM) */
539 points[si].topomin=1;
540 /* and copy it to G point list (if there is space left) */
541 if(topoNr<sampleNr) {
542 for(unsigned int i=0; i<nlo->totalNr; i++) gpoints[topoNr].par[i]=points[si].par[i];
543 gpoints[topoNr++].fvalue=points[si].fvalue;
544 }
545 /* Print TM point */
546 if(verbose>6) {
547 printf("TM point %d\t", 1+si);
548 for(unsigned int i=0; i<nlo->totalNr; i++) printf("%e ", points[si].par[i]);
549 printf("=> %e\n", points[si].fvalue);
550 fflush(stdout);
551 }
552 }
553 } // inner loop
554
555
556 if(verbose>3) {printf(" process gpoints.\n"); fflush(stdout);}
557 /* For each point i, find out if it is a "topographic minimum", that is,
558 better than k/2 neighbour points. */
559 for(unsigned int si=0; si<sampleNr; si++) {
560 gpoints[si].topomin=0;
561
562 /* Compute the distances, scaled using xtol[] */
563 double sdist[sampleNr];
564 for(unsigned int sj=0; sj<sampleNr; sj++) {
565 sdist[sj]=0.0;
566 if(si==sj) {sdist[sj]=1.0E+99; continue;}
567 for(unsigned int k=0; k<nlo->totalNr; k++) {
568 if(!(nlo->xtol[k]>0.0)) continue;
569 double d=(gpoints[si].par[k]-gpoints[sj].par[k])/nlo->xtol[k];
570 if(isfinite(d)) sdist[sj]+=d*d;
571 }
572 }
573
574 /* Find the closest neighbours for sample i */
575 unsigned int ni=0;
576 for(ni=0; ni<neighNr; ni++) {
577 unsigned int mini=0;
578 double minv=sdist[mini];
579 for(unsigned int sj=0; sj<sampleNr; sj++) {
580 if(sdist[sj]<minv) {minv=sdist[sj]; mini=sj;}
581 }
582 sdist[mini]=1.0E+99; // not used again
583 /* If point i is worse than any of the closest neighbours, then point i is certainly not
584 a topographic minimum; then stop this loop and go to the next point i+1 */
585 if(!(gpoints[si].fvalue<gpoints[mini].fvalue)) break;
586 }
587 if(ni<neighNr) continue; // continue with the next sample i
588 /* otherwise mark this as topographic minimum (TM) */
589 gpoints[si].topomin=1;
590 /* Print TM point */
591 if(verbose>6) {
592 printf("GTM point %d\t", 1+si);
593 for(unsigned int i=0; i<nlo->totalNr; i++) printf("%e ", gpoints[si].par[i]);
594 printf("=> %e\n", gpoints[si].fvalue);
595 fflush(stdout);
596 }
597 }
598 /* Copy TMs into the start of gpoint list and set topoNr accordingly */
599 topoNr=0;
600 for(unsigned int si=0; si<sampleNr; si++) if(gpoints[si].topomin>0) {
601 if(si>topoNr) {
602 for(unsigned int i=0; i<nlo->totalNr; i++) gpoints[topoNr].par[i]=gpoints[si].par[i];
603 gpoints[topoNr].fvalue=gpoints[si].fvalue;
604 }
605 topoNr++;
606 }
607 if(verbose>3) {printf(" %d topographical minima\n", topoNr); fflush(stdout);}
608 if(topoNr==0) {
609 if(verbose>0) {fprintf(stderr, "Warning: no topographical minima found\n"); fflush(stderr);}
610 continue;
611 }
612
613
614 } // end of iTGO
615
616
617 /* If requested, do local optimization for the TMs */
618 if(topoNr==0) {
619 if(verbose>0) {fprintf(stderr, "Error: no topographical minima found\n"); fflush(stderr);}
620 free(points); free(gpoints);
621 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_BAD_FIT);
622 return TPCERROR_BAD_FIT;
623 }
624 for(unsigned int si=0; si<sampleNr; si++) if(gpoints[si].topomin>0) {
625 if(verbose>2) printf("LO for TM point %d\n", 1+si);
626 if(verbose>3) {
627 printf("TM point %d before LO:\n", 1+si);
628 printf("%d\t", 1+si);
629 for(unsigned int i=0; i<nlo->totalNr; i++) printf("%e ", gpoints[si].par[i]);
630 printf("=> %e\n", gpoints[si].fvalue); fflush(stdout);
631 }
632 NLOPT lnlo; nloptInit(&lnlo);
633 nloptDuplicate(nlo, &lnlo);
634 for(unsigned int i=0; i<nlo->totalNr; i++) lnlo.xfull[i]=gpoints[si].par[i];
635 lnlo.funval=points[si].fvalue;
636 for(unsigned int i=0; i<nlo->totalNr; i++) lnlo.xdelta[i]=100.0*lnlo.xtol[i];
637 lnlo.maxFunCalls=2000*nlo->totalNr; lnlo.funCalls=0;
638 int ret=nloptSimplex(&lnlo, 100*nlo->totalNr, status);
639 if(ret==TPCERROR_OK) {
640 for(unsigned int i=0; i<nlo->totalNr; i++) gpoints[si].par[i]=lnlo.xfull[i];
641 gpoints[si].fvalue=lnlo.funval;
642 if(verbose>3) {
643 printf("TM point %d after LO:\n", 1+si);
644 printf("%d\t", 1+si);
645 for(unsigned int i=0; i<nlo->totalNr; i++) printf("%e ", gpoints[si].par[i]);
646 printf("=> %e\n", gpoints[si].fvalue); fflush(stdout);
647 }
648 } else if(verbose>7) {printf(" LO failed\n"); fflush(stdout);}
649 nlo->funCalls+=lnlo.funCalls;
650 nloptFree(&lnlo);
651 }
652
653 /* Find the best minimum and return it */
654 {
655 unsigned int mini=0;
656 double minv=gpoints[mini].fvalue;
657 for(unsigned int si=1; si<sampleNr; si++)
658 if(!(gpoints[si].fvalue>minv)) {minv=gpoints[si].fvalue; mini=si;}
659 for(unsigned int i=0; i<nlo->totalNr; i++) nlo->xfull[i]=gpoints[mini].par[i];
660 nlo->funval=gpoints[mini].fvalue;
661 }
662
663 free(points); free(gpoints);
664 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
665 return TPCERROR_OK;
666}
Generated on for TPCCLIB by doxygen 1.15.0