Log In | Get Help   
Home My Page Projects Code Snippets Project Openings Mareframe
Summary Activity Forums Tracker Lists Tasks Docs Surveys News SCM Files
[mareframe] Diff of /trunk/gadget/simann.cc
[mareframe] / trunk / gadget / simann.cc Repository:
ViewVC logotype

Diff of /trunk/gadget/simann.cc

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 11, Thu Jul 23 19:00:38 2015 UTC revision 12, Fri Jul 24 18:36:24 2015 UTC
# Line 186  Line 186 
186  #ifndef NO_OPENMP  #ifndef NO_OPENMP
187  extern Ecosystem** EcoSystems;  extern Ecosystem** EcoSystems;
188  #endif  #endif
 //extern StochasticData* data;  
   
189    
190    /*sequential code replaced at seq_optimize_template.h*/
191  //void OptInfoSimann::OptimiseLikelihood() {  //void OptInfoSimann::OptimiseLikelihood() {
192  //  //
193  //  //set initial values  //  //set initial values
# Line 507  Line 506 
506      handle.logMessage(LOGINFO, "Error starting Simulated Annealing optimisation with f(x) = infinity");      handle.logMessage(LOGINFO, "Error starting Simulated Annealing optimisation with f(x) = infinity");
507      converge = -1;      converge = -1;
508      iters = 1;      iters = 1;
 //    EcoSystem->writeOptValuesOMP();  
509      return;      return;
510    }    }
511    
# Line 652  Line 650 
650                x[i] = storage[bestId].trialx[i];                x[i] = storage[bestId].trialx[i];
651              funcval = trialf;              funcval = trialf;
652              nacc++;              nacc++;
 //            vns[param[l+bestId]]++;  
             //nacp[param[l+bestId]]+=numThr;  //JMB - not sure about this ...  
 //            for (i = 0; i < numThr; ++i)  
 //              nacp[param[l+i]]++;  
 //            nacp[param[l+bestId]]++;  
   
653            } else {            } else {
654              //Accept according to metropolis condition              //Accept according to metropolis condition
655              p = expRep((trialf - funcval) / t);              p = expRep((trialf - funcval) / t);
# Line 669  Line 661 
661                funcval = trialf;                funcval = trialf;
662                naccmet++;                naccmet++;
663                nacp[param[(l+bestId)%nvars]]++;                nacp[param[(l+bestId)%nvars]]++;
 //              vns[param[l+bestId]]++;  
 //              nacp[param[l+bestId]]+=numThr;  
 //              for (i = 0; i < numThr; ++i)  
 //                nacp[param[l+i]]++;  
 //                nacp[param[l+bestId]]++;  
664              } else {              } else {
665                //Reject point                //Reject point
666                nrej++;                nrej++;
# Line 683  Line 670 
670            if (isZero(trialf)) {            if (isZero(trialf)) {
671              handle.logMessage(LOGINFO, "Error in Simulated Annealing optimisation after", iters, "function evaluations, f(x) = 0");              handle.logMessage(LOGINFO, "Error in Simulated Annealing optimisation after", iters, "function evaluations, f(x) = 0");
672              converge = -1;              converge = -1;
 //            EcoSystem->writeOptValuesOMP();  
673              delete[] storage;              delete[] storage;
674              return;              return;
675            }            }
# Line 697  Line 683 
683              if (scale) {              if (scale) {
684                for (i = 0; i < nvars; ++i)                for (i = 0; i < nvars; ++i)
685                  scalex[i] = bestx[i] * init[i];                  scalex[i] = bestx[i] * init[i];
               //FIXME EcoSystems[]->  
686                EcoSystem->storeVariables(-fopt, scalex);                EcoSystem->storeVariables(-fopt, scalex);
687              } else              } else
688                  EcoSystem->storeVariables(-fopt, bestx);                  EcoSystem->storeVariables(-fopt, bestx);
# Line 708  Line 693 
693            }            }
694          }          }
695        }        }
   
       //Adjust vm so that approximately half of all evaluations are accepted  
 //      for (i = 0; i < nvars; ++i) {  
 ////              cout << "nacp["<<i<<"]:"<<nacp[i]<<endl;  
 //        ratio = nsdiv * nacp[i];  
 //        nacp[i] = 0;  
 //        if (ratio > uratio) {  
 //          vm[i] = vm[i] * (1.0 + cs * (ratio - uratio));  
 //        } else if (ratio < lratio) {  
 //          vm[i] = vm[i] / (1.0 + cs * (lratio - ratio));  
 //        }  
 //  
 //        if (vm[i] < rathersmall)  
 //          vm[i] = rathersmall;  
 //        if (vm[i] > (upperb[i] - lowerb[i]))  
 //          vm[i] = upperb[i] - lowerb[i];  
 //  
 ////        cout << "vm["<<i<<"]:"<<vm[i]<<endl;  
 //  
 //      }  
696      }      }
697    
698      //Check termination criteria      //Check termination criteria
# Line 758  Line 723 
723        converge = 1;        converge = 1;
724        score = EcoSystem->SimulateAndUpdate(bestx);        score = EcoSystem->SimulateAndUpdate(bestx);
725        handle.logMessage(LOGINFO, "\nSimulated Annealing finished with a likelihood score of", score);        handle.logMessage(LOGINFO, "\nSimulated Annealing finished with a likelihood score of", score);
 //      EcoSystem->writeOptValuesOMP();  
726        delete[] storage;        delete[] storage;
727        return;        return;
728      }      }
# Line 773  Line 737 
737      for (i = 0; i < nvars; ++i)      for (i = 0; i < nvars; ++i)
738        x[i] = bestx[i];        x[i] = bestx[i];
739    }    }
 //  EcoSystem->writeOptValuesOMP();  
740  }  }
741  #endif  #endif
742    
# Line 812  Line 775 
775                  DoubleVector& x, DoubleVector& lowerb, DoubleVector& upperb, DoubleVector& vm, unsigned* seed)                  DoubleVector& x, DoubleVector& lowerb, DoubleVector& upperb, DoubleVector& vm, unsigned* seed)
776  {  {
777          int i, k;          int i, k;
 //      double old = trialx[param[l]];  
 //      cout << "[" << endl;  
778            for (i = 0; i < nvars; ++i) {            for (i = 0; i < nvars; ++i) {
779                  if (i == param[l]) {                  if (i == param[l]) {
780                    trialx[i] = x[i] + ((randomNumber(&*seed) * 2.0) - 1.0) * vm[i];                    trialx[i] = x[i] + ((randomNumber(&*seed) * 2.0) - 1.0) * vm[i];
# Line 831  Line 792 
792                    }                    }
793                  } else                  } else
794                    trialx[i] = x[i];                    trialx[i] = x[i];
 //              cout  << trialx[i]<<" ";  
795            }            }
   
 //        cout << endl <<l<< "-old[" <<param[l] << "]:" << old <<"-" << trialx[param[l]]<< " - " << vm[param[l]]<<endl;  
796  }  }
797    
798  void buildNewParams_f(Siman& seed, DoubleVector& params) {  void buildNewParams_f(Siman& seed, DoubleVector& params) {
# Line 850  Line 808 
808  #ifndef NO_OPENMP  #ifndef NO_OPENMP
809          int id = omp_get_thread_num ();          int id = omp_get_thread_num ();
810          trialf = EcoSystems[id]->SimulateAndUpdate(params);          trialf = EcoSystems[id]->SimulateAndUpdate(params);
 //      cout << id << "-";  
811  #else  #else
812          trialf = EcoSystem->SimulateAndUpdate(params);          trialf = EcoSystem->SimulateAndUpdate(params);
813  #endif  #endif
 //      cout << "trailf:" << trialf<< endl;  
814    return -trialf;    return -trialf;
815  }  }
816    
# Line 872  Line 828 
828                  double uratio = seed.getUratio();                  double uratio = seed.getUratio();
829                  double lratio = seed.getLratio();                  double lratio = seed.getLratio();
830    
         //      cout << "uratio:" << uratio << endl;  
         //      cout << "lratio:" << lratio << endl;  
   
         //cout << "adjust vm" << endl;  
831                  for (i = 0; i < seed.getNvars(); i++) {                  for (i = 0; i < seed.getNvars(); i++) {
 //                      cout << "nacp[" << i << "]:" << seed.getNacp(i) << endl;  
832                          ratio = nsdiv * seed.getNacp(i);                          ratio = nsdiv * seed.getNacp(i);
833                          seed.setNacp(i,0);                          seed.setNacp(i,0);
834                          if (ratio > uratio) {                          if (ratio > uratio) {
# Line 906  Line 857 
857    
858                  DoubleVector* bestx = seed.getBestx();                  DoubleVector* bestx = seed.getBestx();
859    
                 //funcval = fopt;  
860                  for (i = 0; i < seed.getNvars(); i++)                  for (i = 0; i < seed.getNvars(); i++)
861                          x[i] = (*bestx)[i];                          x[i] = (*bestx)[i];
862          }          }
# Line 945  Line 895 
895            //Accept the new point if the new function value better            //Accept the new point if the new function value better
896            bool ret = true;            bool ret = true;
897            int i;            int i;
 //        cout << "mustAccept:" << trialf << "|" << funcval<< endl;  
898            int aux = siman.getParam()[siman.getL()];            int aux = siman.getParam()[siman.getL()];
899          if ((trialf - funcval) > verysmall) {          if ((trialf - funcval) > verysmall) {
900                  siman.incrementNacp(aux);                  siman.incrementNacp(aux);
# Line 955  Line 904 
904                  double p, pp;                  double p, pp;
905                  //Accept according to metropolis condition                  //Accept according to metropolis condition
906    
 //              cout << "t:" << siman.getT() << endl;  
   
907                  p = expRep((trialf - funcval) / siman.getT());                  p = expRep((trialf - funcval) / siman.getT());
908                  pp = randomNumber(siman.getSeedM());                  pp = randomNumber(siman.getSeedM());
909                  if (pp < p) {                  if (pp < p) {
 //                      cout << pp << "<" << p << endl;  
910                          siman.incrementNacp(aux);                          siman.incrementNacp(aux);
911                          siman.incrementNaccmet();                          siman.incrementNaccmet();
912  //                      //Accept point  //                      //Accept point
 //                      for (i = 0; i < nvars; i++)  
 //                              x[i] = trialx[i];  
 //                      funcval = trialf;  
 //                      naccmet++;  
 //                      nacp[param[l]]++;  
913                  } else {                  } else {
914                          //Reject point                          //Reject point
915                          ret = false;                          ret = false;
916                          siman.incrementNrej();                          siman.incrementNrej();
917                  }                  }
918          }          }
 //      cout << "nacp[" << aux << "]:" << siman.getNacp(aux) << endl;  
919          // JMB added check for really silly values          // JMB added check for really silly values
920          if (isZero(trialf)) {          if (isZero(trialf)) {
921                  handle.logMessage(LOGINFO, "Error in Simulated Annealing optimisation after",                  handle.logMessage(LOGINFO, "Error in Simulated Annealing optimisation after",
922                                  iters, "function evaluations, f(x) = 0");                                  iters, "function evaluations, f(x) = 0");
   
923                  siman.setConverge(-1);                  siman.setConverge(-1);
924                  return false;                  return false;
925          }          }
# Line 989  Line 928 
928          if (siman.getL() == 0)          if (siman.getL() == 0)
929                  siman.incrementNS();                  siman.incrementNS();
930          if (siman.getNS() >= siman.getNs()){          if (siman.getNS() >= siman.getNs()){
   
931                  siman.setNS(0);                  siman.setNS(0);
   
932                  siman.incrementNT();                  siman.incrementNT();
933                  adjustVm(siman);                  adjustVm(siman);
934                  siman.randomizeParams();                  siman.randomizeParams();
 //              if (seed.getNT() >= seed.getNt()){  
 //  
 //              }  
   
935          }          }
936    
937          return ret;          return ret;
# Line 1015  Line 948 
948           bool quit = false;           bool quit = false;
949            if (siman.getNT() >= siman.getNt())            if (siman.getNT() >= siman.getNt())
950            {            {
 //                cout << "nt:" << siman.getNT();  
                   //Check termination criteria  
951                    int i;                    int i;
952                    DoubleVector fstar = siman.getFstar();                    DoubleVector fstar = siman.getFstar();
953    
# Line 1067  Line 998 
998          handle.logMessage(LOGINFO, "Number of rejected points", siman.getNrej());          handle.logMessage(LOGINFO, "Number of rejected points", siman.getNrej());
999          *score = EcoSystem->SimulateAndUpdate(*bestX);          *score = EcoSystem->SimulateAndUpdate(*bestX);
1000          handle.logMessage(LOGINFO, "\nSimulated Annealing finished with a likelihood score of", *score);          handle.logMessage(LOGINFO, "\nSimulated Annealing finished with a likelihood score of", *score);
 //      EcoSystem->writeOptValues();  
1001    }    }
1002  };  };
1003    
1004  // Required  // Required
1005  std::ostream &operator<<(std::ostream &os, const DoubleVector &p)  std::ostream &operator<<(std::ostream &os, const DoubleVector &p)
1006  {  {
 //  os <<  "[ AUCH";  
 //  for (int i = 0; i< NPARAMS; ++i) {  
 //    std::cout << p[i] << ' ';  
 //  }  
 //  std::cout << "]";  
1007          os << "";          os << "";
1008    return os;    return os;
1009  }  }
# Line 1089  Line 1014 
1014          //set initial values          //set initial values
1015    
1016          double tmp, p, pp;          double tmp, p, pp;
1017          double fopt, funcval, trialf;          double funcval, trialf;
1018          int a, i, j, k, l, quit;          int a, i, j, k, l, quit;
1019          int rchange, rcheck, rnumber; //Used to randomise the order of the parameters          int rchange, rcheck, rnumber; //Used to randomise the order of the parameters
1020    
1021          handle.logMessage(LOGINFO,          handle.logMessage(LOGINFO,
1022                          "\nStarting Simulated Annealing optimisation algorithm----\n");                          "\nStarting Simulated Annealing optimisation algorithm\n");
1023          int nvars = EcoSystem->numOptVariables();          int nvars = EcoSystem->numOptVariables();
1024          DoubleVector x(nvars);          DoubleVector x(nvars);
1025          DoubleVector init(nvars);          DoubleVector init(nvars);
# Line 1152  Line 1077 
1077    
1078          //the function is to be minimised so switch the sign of funcval (and trialf)          //the function is to be minimised so switch the sign of funcval (and trialf)
1079          funcval = -funcval;          funcval = -funcval;
 //      offset = EcoSystem->getFuncEval(); //number of function evaluations done before loop  
 //      nacc++;  
1080          cs /= lratio;  //JMB save processing time          cs /= lratio;  //JMB save processing time
 //      nsdiv = 1.0 / ns;  
         fopt = funcval;  
1081          for (i = 0; i < tempcheck; i++)          for (i = 0; i < tempcheck; i++)
1082                  fstar[i] = funcval;                  fstar[i] = funcval;
1083    
         //unsigned seed = 1234;  
   
1084          Siman s(seed, seedM, seedP, nvars, nt, ns, param, &x, &lowerb, &upperb, vm, t, rt, (1.0 / ns),          Siman s(seed, seedM, seedP, nvars, nt, ns, param, &x, &lowerb, &upperb, vm, t, rt, (1.0 / ns),
1085                          tempcheck, simanneps, fstar, lratio, uratio, cs, &bestx, scale, &converge, &score);                          tempcheck, simanneps, fstar, lratio, uratio, cs, &bestx, scale, &converge, &score);
1086    
1087          ReproducibleSearch<Siman, DoubleVector, ControlClass, evaluate_f, buildNewParams_f>          ReproducibleSearch<Siman, DoubleVector, ControlClass, evaluate_f, buildNewParams_f>
1088          pa(s, x, simanniter);          pa(s, x, simanniter);
1089    
 //      cout << "Sequential run" << endl;  
 //      pa.seq_opt();  
   
1090  #ifndef NO_OPENMP  #ifndef NO_OPENMP
1091            // OpenMP parallelization
1092          int numThr = omp_get_max_threads ( );          int numThr = omp_get_max_threads ( );
1093          pa.paral_opt_omp(numThr,numThr);          pa.paral_opt_omp(funcval,numThr,numThr);
1094  #else  #else
1095          pa.seq_opt();          // sequential code
1096            pa.seq_opt(funcval);
1097  #endif  #endif
1098    
1099  }  }
Legend:
Removed from v.11  
changed lines
  Added in v.12
root@forge.cesga.es
 ViewVC Help
Powered by ViewVC 1.0.0  

Powered By FusionForge