/* NIGHTFALL Light Curve Synthesis Program */
/* Copyright (C) 1998 Rainer Wichmann */
/* */
/* This program is free software; you can redistribute it */
/* and/or modify */
/* it under the terms of the GNU General Public License as */
/* published by */
/* the Free Software Foundation; either version 2 of the License, or */
/* (at your option) any later version. */
/* */
/* This program is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
/* GNU General Public License for more details. */
/* */
/* You should have received a copy of the GNU General Public License */
/* along with this program; if not, write to the Free Software */
/* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <float.h>
#include "Light.h"
#include "LightSimplex.h"
ParStr Param[SDIM] = {
{ LIM_MQ_L, LIM_MQ_H, 0.0, NULL, "Mass Ratio", 1.0 }, /* 0 */
{ LIM_IN_L, LIM_IN_H, 0.0, NULL, "Inclination", 1.0 },
{ LIM_RF_L, LIM_RF_H, 0.0, NULL, "RocheFill (Primary)", 1.0 },
{ LIM_RF_L, LIM_RF_H, 0.0, NULL, "RocheFill (Secondary)", 1.0 },
{ LIM_TB_L, LIM_TB_H, 0.0, NULL, "Temperature (Primary)", 1.0 },
{ LIM_TB_L, LIM_TB_H, 0.0, NULL, "Temperature (Secondary)", 1.0 }, /* 5 */
{ LIM_EX_L, LIM_EX_H, 0.0, NULL, "Excentricity", 1.0 },
{ LIM_PA_L, LIM_PA_H, 0.0, NULL, "Periastron Length", 1.0 },
{ LIM_FR_L, LIM_FR_H, 0.0, NULL, "Asynchronicity (Primary)", 1.0 },
{ LIM_FR_L, LIM_FR_H, 0.0, NULL, "Asynchronicity (Primary)", 1.0 },
{ LIM_SLO_L, LIM_SLO_H, 0.0, NULL, "Longitude (Sp1 P)", 1.0 }, /* 10 */
{ LIM_SLA_L, LIM_SLA_H, 0.0, NULL, "Latitude (Sp1 P)", 1.0 },
{ LIM_SRA_L, LIM_SRA_H, 0.0, NULL, "Radius (Sp1 P)", 1.0 },
{ LIM_SDF_L, LIM_SDF_H, 0.0, NULL, "Dimfactor (Sp1 P)", 1.0 },
{ LIM_SLO_L, LIM_SLO_H, 0.0, NULL, "Longitude (Sp2 P)", 1.0 },
{ LIM_SLA_L, LIM_SLA_H, 0.0, NULL, "Latitude (Sp2 P)", 1.0 }, /* 15 */
{ LIM_SRA_L, LIM_SRA_H, 0.0, NULL, "Radius (Sp2 P)", 1.0 },
{ LIM_SDF_L, LIM_SDF_H, 0.0, NULL, "Dimfactor (Sp2 P)", 1.0 },
{ LIM_SLO_L, LIM_SLO_H, 0.0, NULL, "Longitude (Sp1 S)", 1.0 },
{ LIM_SLA_L, LIM_SLA_H, 0.0, NULL, "Latitude (Sp1 S)", 1.0 },
{ LIM_SRA_L, LIM_SRA_H, 0.0, NULL, "Radius (Sp1 S)", 1.0 }, /* 20 */
{ LIM_SDF_L, LIM_SDF_H, 0.0, NULL, "Dimfactor (Sp1 S)", 1.0 },
{ LIM_SLO_L, LIM_SLO_H, 0.0, NULL, "Longitude (Sp2 S)", 1.0 },
{ LIM_SLA_L, LIM_SLA_H, 0.0, NULL, "Latitude (Sp2 S)", 1.0 },
{ LIM_SRA_L, LIM_SRA_H, 0.0, NULL, "Radius (Sp2 S)", 1.0 },
{ LIM_SDF_L, LIM_SDF_H, 0.0, NULL, "Dimfactor (Sp2 S)", 1.0 }, /* 25 */
{ LIM_MASS_L, LIM_MASS_H, 0.0, NULL, "Total Mass", 1.0 },
{ LIM_DIST_L, LIM_DIST_H, 0.0, NULL, "Distance", 1.0 },
{ LIM_3L_L, LIM_3L_H, 0.0, NULL, "Third Light (U)", 1.0 },
{ LIM_3L_L, LIM_3L_H, 0.0, NULL, "Third Light (B)", 1.0 },
{ LIM_3L_L, LIM_3L_H, 0.0, NULL, "Third Light (V)", 1.0 }, /* 30 */
{ LIM_3L_L, LIM_3L_H, 0.0, NULL, "Third Light (R)", 1.0 },
{ LIM_3L_L, LIM_3L_H, 0.0, NULL, "Third Light (I)", 1.0 },
{ LIM_3L_L, LIM_3L_H, 0.0, NULL, "Third Light (J)", 1.0 },
{ LIM_3L_L, LIM_3L_H, 0.0, NULL, "Third Light (H)", 1.0 },
{ LIM_3L_L, LIM_3L_H, 0.0, NULL, "Third Light (K)", 1.0 }, /* 35 */
{ LIM_3L_L, LIM_3L_H, 0.0, NULL, "Third Light (u)", 1.0 },
{ LIM_3L_L, LIM_3L_H, 0.0, NULL, "Third Light (v)", 1.0 },
{ LIM_3L_L, LIM_3L_H, 0.0, NULL, "Third Light (b)", 1.0 },
{ LIM_3L_L, LIM_3L_H, 0.0, NULL, "Third Light (y)", 1.0 },
{ LIM_OR_L, LIM_OR_H, 0.0, NULL, "Outer Radius", 1.0 }, /* 40 */
{ LIM_IR_L, LIM_IR_H, 0.0, NULL, "Inner Radius", 1.0 },
{ LIM_TH_L, LIM_TH_H, 0.0, NULL, "Thickness", 1.0 },
{ LIM_HR_L, LIM_HR_H, 0.0, NULL, "H/R", 1.0 },
{ LIM_DT_L, LIM_DT_H, 0.0, NULL, "Temperature (Disk)", 1.0 },
{ LIM_HST_L, LIM_HST_H, 0.0, NULL, "T (hot spot)", 1.0 }, /* 45 */
{ LIM_HSL_L, LIM_HSL_H, 0.0, NULL, "Longitude (hot spot)", 1.0 },
{ LIM_HSE_L, LIM_HSE_H, 0.0, NULL, "Extent (hot spot)", 1.0 },
{ LIM_HSD_L, LIM_HSD_H, 0.0, NULL, "Depth (hot spot)", 1.0 }
};
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Downhill siplex code
@long Kallrath J., Linnel A.P. 1987, Astrophysical Journal 313, 346
Nelder J.A., Mead R. 1965, Comput. J. 7, 308
@param (void)
@return (int) status The exit status
@heading Data Fitting
****************************************************************************/
int Simplex()
{
int ErrCode = 0; /* exit status subroutines */
double Y[SDIM+1] = {0.0}; /* merit (result) vector */
double X[SDIM] = {0.0}; /* some Vertex */
double BestEver[SDIM] = {0.0}; /* best vertex found */
double BestChi = 1E38; /* best chi found */
double S[SDIM+1][SDIM]; /* The Simplex */
double Centre[SDIM] = {0.0}; /* Simplex Centre */
int FitDim; /* Number of Fitted Parameters */
int Worst = 0; /* Index of worst Vertex */
int NextWorst = 0; /* Index of second worst Vertex*/
int Best = 0; /* Index of best Vertex */
double Sum; /* temporary sum */
long NumIter = 0; /* Iteration count */
long NumCycle = 0; /* Cycle count */
int i, j; /* Loop Variables */
double Y_Test; /* new merit */
double Y_Store; /* temporary storage */
int VertexTest; /* check whether V. in bounds */
double MeanChiSquare; /* mean of merits */
double VarChiSquare; /* standard deviation of merits*/
float ReflectFactor = -1.0; /* a reflection */
float ExpandFactor = 2.0; /* an expansion */
float ContractFactor = 0.35; /* a contraction */
float ShrinkFactor = 0.5; /* shrinkage */
double Yb, Yw, Ynw; /* temp variables */
long dof; /* degrees of freedom */
long MaxIter = MAX_ITER_SIMPLEX; /* maximum iterations */
float TempZero = TEMP_ANNEAL_SCALE;/* start temperature SA */
float TempScale; /* cooling SA */
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON) {
if (Flags.anneal == OFF) {
printf("\n");
printf(_(" Starting Downhill Simplex Fit\n"));
}
else {
printf("\n");
printf(_(" Starting Simulated Annealing\n"));
}
}
/* ---------------- initialize ---------------------------------- */
for (j = 0; j < SDIM+1; ++j) {
for(i = 0; i < SDIM; ++i) {
S[j][i] = 0.0;
}
}
FitDim = 0;
for (j = 0; j < SDIM; ++j) {
if (Flags.simplex[j] == ON) ++FitDim;
}
dof = 0;
for (j = 0; j < NUM_MAG+2; ++j) {
dof = dof + Flags.Passbands[j];
}
if (Flags.debug[OPTIMIZE] == ON)
fprintf(stderr, _("Data = %5ld, Fit Parameters = %5d\n"), dof, FitDim);
dof = dof - FitDim;
if (FitDim == 0 || dof <= 0) {
#ifdef _WITH_GTK
if ( Flags.interactive == ON) {
make_dialog(_(errmsg[18]));
} else {
WARNING(_(" Dear User,"));
WARNING(_(" We have a minor problem: "));
WARNING(_(" There is nothing to fit, thus we abort"));
WARNING(_(" the fitting procedure"));
}
return(8);
#else
WARNING(_(" Dear User,"));
WARNING(_(" We have a minor problem: "));
WARNING(_(" There is nothing to fit, thus we abort"));
WARNING(_(" the fitting procedure"));
return(8);
#endif
}
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON) {
if (Flags.anneal == OFF) {
printf("\n");
printf(_(" Initialize Simplex \n"));
}
else {
printf("\n");
printf(_(" Initialize Simplex for Simulated Annealing \n"));
}
}
#ifdef _WITH_GTK
if (Flags.interactive == ON && Flags.anneal == OFF)
my_appbar_push(_("Initialize Simplex"));
else if (Flags.interactive == ON && Flags.anneal == ON)
my_appbar_push(_("Initialize Simplex for Simulated Annealing"));;
#endif
/* ----- Init and find the Best Vertex (for annealing) ----------- */
SimplexInitRanges();
ErrCode = SimplexInit(Y, S);
#ifdef _WITH_GTK
if ( Flags.interactive == ON && ErrCode > 0) {
make_dialog(_(errmsg[19]));
}
#endif
for (j = 0; j < (SDIM+1); ++j) { /* init to some valid index */
if (fabs(Y[j]) >= FLT_EPSILON) { BestChi = Y[j]; Best = j; }
}
for (j = 0; j < (SDIM+1); ++j) {
if ( fabs(Y[j]) >= FLT_EPSILON && Y[j] <= BestChi ) {
BestChi = Y[j]; Best = j;
}
}
for (j = 0; j < SDIM; ++j) {
if (Flags.simplex[j] == ON) BestEver[j] = S[Best][j];
}
if (Flags.debug[OPTIMIZE] == ON)
fprintf(stderr, _("FitDim after Init = %5d\n"), FitDim);
MaxIter = MAX_ITER_SIMPLEX;
if ((Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON) )
printf(_(" Simplex Initialized\n"));
/* >>>>>>>>>>>>>> do simulated annealing <<<<<<<<<<<<<<<<<<<<< */
for (j = 0; j < SDIM; ++j) {
if (fabs(Y[j]) >= FLT_EPSILON && Y[j] <= (PENALTY/2.) )
TempZero = TempZero + Y[j];
}
if (TempZero >= FLT_EPSILON) TempZero = TempZero / FitDim;
else TempZero = TEMP_ANNEAL_SCALE;
TempScale = (float)(exp( -log((double)TEMP_ANNEAL_SCALE) /(double)FitDim));
TempScale = (float)((-log(TEMP_RATIO_SCALE)) * TempScale);
if (Flags.anneal == ON) {
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON)
printf(_(" Start Simulated Annealing, T0 = %6.2f, Tscale = %6.2f\n"),
TempZero,
TempScale);
for (j = 0; j < SDIM; ++j) X[j] = BestEver[j];
(void) SaAnneal(TempScale, TempZero, X, &Y[Best],
BestEver, &BestChi);
if (Flags.debug[BUSY] == ON) {
for (j = 0; j < SDIM; ++j) printf(" %7.2f", BestEver[j]);
printf("\n");
}
for (j = 0; j < SDIM; ++j) S[Best][j] = BestEver[j];
Y[Best] = BestChi;
/* ----------- compute L.C. with best vertex ------------------ */
SimplexSetVariables(BestEver);
ErrCode = MainLoop(&Y[Best]);
/* ----------------- print out --------------------------------- */
SimplexPrintVariables(BestChi, 0.0);
return(ErrCode);
}
/* >>>>>>>>>>>>>> do simplex fitting <<<<<<<<<<<<<<<<<<<<<<<<< */
/* ------ find centre of simplex ---------------------------- */
for (j = 0; j < SDIM; ++j) {
Sum = 0.0;
for (i = 0; i < (SDIM+1); ++i) Sum = Sum + S[i][j];
if (Flags.simplex[j] == ON) Centre[j] = Sum / (FitDim+1);
else Centre[j] = 0.;
}
/* --------- start iterating ----------------------------------- */
while (1) {
++NumCycle;
/* ----------- best, worst, next_worst vertices ----------- */
for (j = 0; j < (SDIM+1); ++j) { /* init to some valid index */
if (fabs(Y[j]) >= FLT_EPSILON) {
Best = j;
Worst = j;
if (j != Worst) NextWorst = j;
}
}
Yb = Y[Best];
Yw = Y[Worst];
for (j = 0; j < (SDIM+1); ++j) { /* find best and worst */
if (fabs(Y[j]) >= FLT_EPSILON) {
if (Y[j] <= (Yb-FLT_EPSILON)) { Best = j; Yb = Y[j]; }
if (Y[j] >= (Yw+FLT_EPSILON)) { Worst = j; Yw = Y[j]; }
}
}
Ynw = Y[NextWorst];
for (j = 0; j < (SDIM+1); ++j) { /* find next worst */
if (fabs(Y[j]) >= FLT_EPSILON && j != Worst) {
if (Y[j] >= (Ynw+FLT_EPSILON)) { NextWorst = j; Ynw = Y[j]; }
}
}
if (Flags.anneal == OFF) {
for (j = 0; j < SDIM; ++j) BestEver[j] = S[Best][j];
BestChi = Y[Best];
}
/* >>>>>>>>> compute standard deviation, finish <<<<<<<<<<< */
MeanChiSquare = 0.0; VarChiSquare = 0.0;
for (j = 0; j < SDIM+1; ++j) MeanChiSquare = MeanChiSquare + Y[j];
MeanChiSquare = MeanChiSquare / (FitDim+1);
for (j = 0; j < SDIM+1; ++j) {
if (fabs(Y[j]) >= FLT_EPSILON)
VarChiSquare = VarChiSquare+SQR(Y[j]-MeanChiSquare);
}
VarChiSquare = sqrt(VarChiSquare / FitDim);
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON) {
printf("\n");
printf(_("-- Cycle %3ld -- Best: %10.5g SDV (ChiSquare): %10.5g\n"),
NumCycle, Y[Best], VarChiSquare);
}
if (VarChiSquare <= Flags.SimplexTol || NumCycle > MAX_ITER_SIMPLEX) {
/* ----------- compute L.C. with best vertex -------------- */
for (j = 0; j < SDIM; ++j) X[j] = BestEver[j];
SimplexSetVariables(X);
ErrCode = MainLoop(&Y[Best]);
/* ----------------- print out ----------------------------- */
SimplexPrintVariables(BestChi, VarChiSquare);
if (NumCycle > MaxIter
&& VarChiSquare >= (Flags.SimplexTol-FLT_EPSILON)) {
WARNING(_("No Convergence in Fit, Maximum Evaluation Count exceed"))
return(8);
} else {
return(0);
}
}
/* >>>>>>>>> the next simplex iteration <<<<<<<<<<<< */
++NumIter;
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON)
printf(_("FuncEval %3ld: I reflect my worst point\n"), NumIter );
Y_Test = SimplexFlow(ReflectFactor, S, Y, Centre, Worst,
FitDim, &ErrCode);
if (ErrCode == 8) return(8);
if (Y_Test <= (Y[Best]-FLT_EPSILON)) {
++NumIter;
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON)
printf(_("FuncEval %3ld: Reflection was successfull, I expand further\n"),
NumIter);
Y_Test = SimplexFlow(ExpandFactor, S, Y, Centre, Worst,
FitDim, &ErrCode);
if (ErrCode == 8) return(8);
}
else if (Y_Test >= Y[NextWorst]) {
++NumIter;
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON)
printf(_("FuncEval %3ld: No improvement, I try a contraction\n"),
NumIter);
Y_Store = Y[Worst];
Y_Test = SimplexFlow(ContractFactor, S, Y, Centre, Worst,
FitDim, &ErrCode);
if (ErrCode == 8) return(8);
if (Y_Test >= Y_Store) {
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON)
printf(_("Still no improvement, I shrink towards my best point\n"));
/* shrink the simplex */
/* loop over vertices */
for (i = 0; i < SDIM+1; ++i) {
if (i != Best && fabs(Y[i]) >= FLT_EPSILON) {
/* this is not a best vertex */
/* contract and calculate anew */
++NumIter;
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON)
printf(_(" FuncEval %3ld: Computing Vertex %3d\n"),
NumIter,i);
for (j = 0; j < SDIM; ++j) {
if (Flags.simplex[j] == ON) {
X[j] = S[Best][j] + ShrinkFactor * (S[i][j] - S[Best][j]);
}
}
if (Flags.debug[OPTIMIZE] == ON) {
fprintf(stderr, _("Best:") );
for (j = 0; j < SDIM; ++j) fprintf(stderr," %7.2f", S[Best][j]);
fprintf(stderr,"\n");
fprintf(stderr, _("Old:") );
for (j = 0; j < SDIM; ++j) fprintf(stderr," %7.2f", S[i][j]);
fprintf(stderr,"\n");
fprintf(stderr, _("New:") );
for (j = 0; j < SDIM; ++j) fprintf(stderr," %7.2f", X[j]);
fprintf(stderr,"\n");
}
VertexTest = SimplexCheckVertex(X);
for (j = 0; j < SDIM; ++j) S[i][j] = X[j];
if (Flags.debug[OPTIMIZE] == ON) {
for (j = 0; j < SDIM; ++j) fprintf(stderr,"%7.2f",X[j]);
}
if (Flags.debug[OPTIMIZE] == ON) printf("\n");
if (VertexTest > 0) {
Y[i] = PENALTY;
} else {
SimplexSetVariables(X);
ErrCode = MainLoop(&Y[i]);
if (ErrCode == 8) return(8);
}
}
} /* shrink best value finished */
for (j = 0; j < SDIM; ++j) {
Sum = 0.0;
for (i = 0; i < (SDIM+1); ++i) Sum = Sum + S[i][j];
if (Flags.simplex[j] == ON) Centre[j] = Sum / (FitDim+1);
else Centre[j] = 0.;
}
}
}
} /* check convergence, start new loop */
return(0);
} /* END OF PROGRAM */
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Morph the simplex
@param (double) Factor The morphing factor
@param (double) S[SDIM+1][SDIM] The simplex
@param (double) Y[] A vector
@param (double) Centre[] The centre
@param (int) Worst The index of the worst vector
@param (int) FitDim Dimension of simplex
@param (int*) ErrCode The exit status
@return (double) The merit
@heading Data Fitting
****************************************************************************/
double SimplexFlow(double Factor, double S[][SDIM],
double Y[], double Centre[], int Worst,
int FitDim, int *ErrCode)
{
int j; /* loop counter */
double Y_Test; /* the merit */
double X_Test[SDIM]; /* the parameter vector */
int VertexTest; /* exit status subroutine */
*ErrCode = 0;
/* -------- compute parameter vector --------------------------------- */
for (j = 0; j < SDIM; ++j) {
if (Flags.simplex[j] == ON) {
X_Test[j] = Centre[j]
+ Factor * (S[Worst][j] - Centre[j]);
} else {
X_Test[j] = 0.0;
}
}
if (Flags.debug[OPTIMIZE] == ON) {
fprintf(stderr, _("Centre:") );
for (j = 0; j < SDIM; ++j) fprintf(stderr, " %7.2f", Centre[j]);
fprintf(stderr, "\n");
fprintf(stderr, _("Worst:") );
for (j = 0; j < SDIM; ++j) fprintf(stderr, " %7.2f", S[Worst][j]);
fprintf(stderr, "\n");
fprintf(stderr, _("New:") );
for (j = 0; j < SDIM; ++j) fprintf(stderr, " %7.2f", X_Test[j]);
fprintf(stderr, "\n");
}
/* -------- validate and compute --------------------------------- */
VertexTest = SimplexCheckVertex(X_Test);
if (VertexTest > 0) {
Y_Test = PENALTY;
return(Y_Test);
} else {
SimplexSetVariables(X_Test);
if (Flags.debug[OPTIMIZE] == ON) {
for (j = 0; j < SDIM; ++j) printf(" %7.2f", X_Test[j]);
}
if (Flags.debug[OPTIMIZE] == ON) printf("\n");
*ErrCode = MainLoop(&Y_Test);
if (*ErrCode == 8) return(PENALTY);
}
/* -------- update simplex --------------------------------- */
if (Y_Test <= (Y[Worst]-FLT_EPSILON)) {
Y[Worst] = Y_Test;
for (j = 0; j < SDIM; ++j) {
if (Flags.simplex[j] == ON) {
Centre[j] =
((FitDim+1) * Centre[j] - S[Worst][j] + X_Test[j]) / (FitDim+1);
S[Worst][j] = X_Test[j];
}
}
}
return(Y_Test);
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Initialize ranges
@param (void)
@return (void)
@heading Data Fitting
****************************************************************************/
void SimplexInitRanges()
{
int i;
if (Flags.fill == OFF) {
Param[0].Lower = LIM_MQ_L; Param[0].Upper = LIM_MQ_H;
Param[2].Lower = LIM_RF_L; Param[2].Upper = LIM_RF_H;
Param[3].Lower = LIM_RF_L; Param[3].Upper = LIM_RF_H;
} else {
Param[0].Lower = LIM_MQO_L; Param[0].Upper = LIM_MQO_H;
Param[2].Lower = LIM_RO_L; Param[2].Upper = LIM_RO_H;
Param[3].Lower = LIM_RO_L; Param[3].Upper = LIM_RO_H;
}
if (Flags.blackbody == OFF) {
Param[4].Lower = LIM_TM_L;
Param[4].Upper = LIM_TM_H(Binary[Primary].log_g);
Param[5].Lower = LIM_TM_L;
Param[5].Upper = LIM_TM_H(Binary[Secondary].log_g);
} else {
Param[4].Lower = LIM_TB_L; Param[4].Upper = LIM_TB_H;
Param[5].Lower = LIM_TB_L; Param[5].Upper = LIM_TB_H;
}
for (i=0; i<SDIM; ++i) {
Param[i].Range = Param[i].Upper - Param[i].Lower;
Param[i].t = 1.0;
}
Param[0].par = &Binary[Primary].Mq;
Param[1].par = &Orbit.Inclination;
Param[2].par = &Binary[Primary].RocheFill;
Param[3].par = &Binary[Secondary].RocheFill;
Param[4].par = &Binary[Primary].Temperature;
Param[5].par = &Binary[Secondary].Temperature;
Param[6].par = &Orbit.Excentricity;
Param[7].par = &Orbit.Omega;
Param[8].par = &Binary[Primary].Fratio;
Param[9].par = &Binary[Secondary].Fratio;
Param[10].par = &Spot[Primary][0].longitude;
Param[11].par = &Spot[Primary][0].latitude;
Param[12].par = &Spot[Primary][0].radius;
Param[13].par = &Spot[Primary][0].dimfactor;
Param[14].par = &Spot[Primary][1].longitude;
Param[15].par = &Spot[Primary][1].latitude;
Param[16].par = &Spot[Primary][1].radius;
Param[17].par = &Spot[Primary][1].dimfactor;
Param[18].par = &Spot[Secondary][0].longitude;
Param[19].par = &Spot[Secondary][0].latitude;
Param[20].par = &Spot[Secondary][0].radius;
Param[21].par = &Spot[Secondary][0].dimfactor;
Param[22].par = &Spot[Secondary][1].longitude;
Param[23].par = &Spot[Secondary][1].latitude;
Param[24].par = &Spot[Secondary][1].radius;
Param[25].par = &Spot[Secondary][1].dimfactor;
Param[26].par = &Orbit.TrueMass;
Param[27].par = &Orbit.TrueDistance;
Param[28].par = &Orbit.Third[Umag];
Param[29].par = &Orbit.Third[Bmag];
Param[30].par = &Orbit.Third[Vmag];
Param[31].par = &Orbit.Third[Rmag];
Param[32].par = &Orbit.Third[Imag];
Param[33].par = &Orbit.Third[Jmag];
Param[34].par = &Orbit.Third[Hmag];
Param[35].par = &Orbit.Third[Kmag];
Param[36].par = &Orbit.Third[umag];
Param[37].par = &Orbit.Third[vmag];
Param[38].par = &Orbit.Third[bmag];
Param[39].par = &Orbit.Third[ymag];
#ifdef HAVE_DISK
Param[40].par = &Binary[Disk].Rout;
Param[41].par = &Binary[Disk].Rin;
Param[42].par = &Binary[Disk].Thick;
Param[43].par = &Binary[Disk].HR;
Param[44].par = &Binary[Disk].Temperature;
Param[45].par = &Binary[Disk].tempHS;
Param[46].par = &Binary[Disk].longitudeHS;
Param[47].par = &Binary[Disk].extentHS;
Param[48].par = &Binary[Disk].depthHS;
#endif
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Initialize the simplex
@param (double) Y[] A vector
@param (double) S[SDIM+1][SDIM] The simplex
@return (int) The exit status
@heading Data Fitting
****************************************************************************/
int SimplexInit(double Y[], double S[][SDIM])
{
int i, j; /* loop counter */
int vertexnum; /* vertex counter */
int VertexTest; /* exit status */
int ErrCode = 0; /* exit status */
double p, q; /* for vertex definition */
double FitDim; /* # params to fit */
double X[SDIM] = { 0.0 }; /* a vertex */
double Var[SDIM+1] = { 0.0 }; /* init deltas */
double tmp; /* temporary */
/* ------------------- initialize ------------------------------------- */
FitDim = 0;
for (j = 0; j < SDIM; ++j) {
S[0][j] = 0.0;
if (Flags.simplex[j] == ON) {
++FitDim;
S[0][j] = (*Param[j].par);
}
Var[j+1] = Param[j].Range / 10.;
}
tmp = sqrt(2.0*FitDim);
p = (sqrt(FitDim+1.0) + FitDim -1.0) / tmp;
q = (sqrt(FitDim+1.0) - FitDim) / tmp;
Binary[Secondary].Mq = 1.0 / Binary[Primary].Mq;
/* -------------- generate init vertex ---------------------------------- */
for (j = 0; j < SDIM; ++j)
X[j] = S[0][j];
VertexTest = SimplexCheckVertex(X);
if (VertexTest > 0)
{
Y[0] = PENALTY;
}
else
{
SimplexSetVariables(X);
ErrCode = MainLoop(&Y[0]);
if (ErrCode == 8)
/* return(8); */ Y[0] = PENALTY;
}
/* ---------------- define vertices of simplex ------------------------- */
vertexnum = 0;
for (i = 1; i < SDIM+1; ++i) {
if (Flags.simplex[i-1] == ON) {
for (j = 0; j < SDIM; ++j) S[i][j] = S[0][j];
if(Flags.simplex[j] == ON) S[i][j] = S[i][j] + q * Var[i];
S[i][i-1] = S[0][i-1] + p * Var[i];
for (j = 0; j < SDIM; ++j) X[j] = S[i][j];
VertexTest = SimplexCheckVertex(X);
if (VertexTest > 0) {
Y[i] = PENALTY;
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON)
printf(_(" -- Vertex %3d initialized (invalid state)\n"),
vertexnum);
} else {
SimplexSetVariables(X);
ErrCode = MainLoop(&Y[i]);
if (ErrCode == 8) /* return(8); */ Y[i] = PENALTY;
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON)
printf(_(" -- Vertex %3d initialized \n"),
vertexnum);
}
++vertexnum;
} else {
for (j = 0; j < SDIM; ++j) S[i][j] = 0.0;
Y[i] = 0.0;
}
}
return(0);
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Validate a simplex vertex
@param (double) Y[] The vertex
@return (int) The exit status
@heading Data Fitting
****************************************************************************/
int SimplexCheckVertex(double X[])
{
int ErrVal = 0; /* the error to return */
int i; /* loop counter */
for (i = 0; i < SDIM; ++i) {
if (Flags.simplex[i] == ON) {
if (X[i] <= (Param[i].Lower-FLT_EPSILON)
|| X[i] >= (Param[i].Upper+FLT_EPSILON)) {
if (Flags.debug[OPTIMIZE] == ON)
fprintf(stderr,
_("*** Out of range [%8.4g,%8.4g]: %s, Value = %8.4g\n"),
Param[i].Lower,
Param[i].Upper,
Param[i].Name,
X[i]);
ErrVal = i + 1;
return(ErrVal);
}
}
}
/* Now set the Variables and do the last Test(s) */
SimplexSetVariables(X);
/* check whether geometry can be defined properly */
if (Flags.simplex[2] == ON && Flags.fill == ON) {
if (DefineParamOver() == 1) {
ErrVal = 3;
return(ErrVal);
}
}
/* check whether stars intersect */
if ( ( (Binary[Primary].Fratio-1.0) <= (-FLT_EPSILON)
&& Flags.asynchron1 == ON )
|| ( (Binary[Secondary].Fratio-1.0) <= (-FLT_EPSILON)
&& Flags.asynchron2 == ON ) ) {
if (LightSetupTests() == 1) {
WARNING(_("Stars intersect (Large Fill Factor, Rotation Rate < One)"));
ErrVal = 9;
return(ErrVal);
}
}
/* ---------------------- no problems --------------------------------- */
return (0);
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Set the global variables according to vertex
@param (double) Y[] The vertex
@return (void)
@heading Data Fitting
****************************************************************************/
void SimplexSetVariables(double x[])
{
int j;
if (Flags.debug[OPTIMIZE] == ON) {
fprintf(stderr, _("SimplexSetVariables: Vector\n") );
for (j = 0; j < SDIM; ++j) fprintf(stderr, " %9.4f", x[j]);
fprintf(stderr, "\n");
fprintf(stderr, _("Flags\n") );
for (j = 0; j < 26; ++j) fprintf(stderr, " %3d", Flags.simplex[j]);
fprintf(stderr, "\n");
}
for (j=0; j < SDIM; ++j)
if (Flags.simplex[j] == ON) (*Param[j].par) = x[j];
Binary[Secondary].Mq = 1.0 / Binary[Primary].Mq;
ComputeGravDark();
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Print fitted values
@param (FILE*) file_out The output file
@return (void)
@heading Data Fitting
****************************************************************************/
void printFitted(FILE *file_out)
{
char fitted[2][16];
strcpy(fitted[0], _("CONSTANT") );
strcpy(fitted[1], _("FITTED") );
fprintf(file_out, " \n");
fprintf(file_out, " ----------------------------------------------------------\n");
fprintf(file_out,
_(" | %10s | Mass Ratio | %10.3f |\n"),
fitted[Flags.simplex[0]], Binary[Primary].Mq);
fprintf(file_out,
_(" | %10s | Inclination | %10.3f |\n"),
fitted[Flags.simplex[1]], RTOD*Orbit.Inclination);
fprintf(file_out,
_(" | %10s | RocheFill (Primary) | %10.3f |\n"),
fitted[Flags.simplex[2]], Binary[Primary].RocheFill);
if (Flags.fill == OFF) {
fprintf(file_out,
_(" | %10s | RocheFill (Secondary) | %10.3f |\n"),
fitted[Flags.simplex[3]], Binary[Secondary].RocheFill);
}
fprintf(file_out,
_(" | %10s | Temperature (Primary) | %10.3f |\n"),
fitted[Flags.simplex[4]], Binary[Primary].Temperature);
fprintf(file_out,
_(" | %10s | Temperature (Secondary) | %10.3f |\n"),
fitted[Flags.simplex[5]], Binary[Secondary].Temperature);
fprintf(file_out,
_(" | %10s | Eccentricity | %10.3f |\n"),
fitted[Flags.simplex[6]], Orbit.Excentricity);
fprintf(file_out,
_(" | %10s | Periastron Length | %10.3f |\n"),
fitted[Flags.simplex[7]],Orbit.Omega);
if (Flags.asynchron1 == ON)
fprintf(file_out,
_(" | %10s | Asynchronicity (Primary) | %10.3f |\n"),
fitted[Flags.simplex[8]], Binary[Primary].Fratio);
if (Flags.asynchron2 == ON)
fprintf(file_out,
_(" | %10s | Asynchronicity (Secondary) | %10.3f |\n"),
fitted[Flags.simplex[9]], Binary[Secondary].Fratio);
#ifdef HAVE_DISK
if (Flags.disk == ON)
{
fprintf(file_out, " ----------------------------------------------------------\n");
fprintf(file_out,
_(" | %10s | Outer Radius (Disk) | %10.3f |\n"),
fitted[Flags.simplex[40]], Binary[Disk].Rout);
fprintf(file_out,
_(" | %10s | Inner Radius (Disk) | %10.3f |\n"),
fitted[Flags.simplex[41]], Binary[Disk].Rin);
fprintf(file_out,
_(" | %10s | Thickness (Disk) | %10.3f |\n"),
fitted[Flags.simplex[42]], Binary[Disk].Thick);
fprintf(file_out,
_(" | %10s | H/R (Disk) | %10.3f |\n"),
fitted[Flags.simplex[43]], Binary[Disk].HR);
fprintf(file_out,
_(" | %10s | Temperature (Disk) | %10.3f |\n"),
fitted[Flags.simplex[44]], Binary[Disk].Temperature);
fprintf(file_out,
_(" | %10s | Temperature (hot spot) | %10.3f |\n"),
fitted[Flags.simplex[45]], Binary[Disk].tempHS);
fprintf(file_out,
_(" | %10s | Longitude (hot spot) | %10.3f |\n"),
fitted[Flags.simplex[46]], Binary[Disk].longitudeHS);
fprintf(file_out,
_(" | %10s | Extent (hot spot) | %10.3f |\n"),
fitted[Flags.simplex[47]], Binary[Disk].extentHS);
fprintf(file_out,
_(" | %10s | Depth (hot spot) | %10.3f |\n"),
fitted[Flags.simplex[48]], Binary[Disk].depthHS);
}
#endif
fprintf(file_out, " ----------------------------------------------------------\n");
fprintf(file_out,
_(" | %10s | Total Mass (Kg) | %10.3e |\n"),
fitted[Flags.simplex[26]], Orbit.TrueMass);
fprintf(file_out,
_(" | %10s | Distance (m) | %10.3e |\n"),
fitted[Flags.simplex[27]], Orbit.TrueDistance);
fprintf(file_out, " ----------------------------------------------------------\n");
fprintf(file_out,
_(" | %10s | ThirdLight (U) | %10.3f |\n"),
fitted[Flags.simplex[28]], Orbit.Third[Umag]);
fprintf(file_out,
_(" | %10s | ThirdLight (B) | %10.3f |\n"),
fitted[Flags.simplex[29]], Orbit.Third[Bmag]);
fprintf(file_out,
_(" | %10s | ThirdLight (V) | %10.3f |\n"),
fitted[Flags.simplex[30]], Orbit.Third[Vmag]);
fprintf(file_out,
_(" | %10s | ThirdLight (R) | %10.3f |\n"),
fitted[Flags.simplex[31]], Orbit.Third[Rmag]);
fprintf(file_out,
_(" | %10s | ThirdLight (I) | %10.3f |\n"),
fitted[Flags.simplex[32]], Orbit.Third[Imag]);
fprintf(file_out,
_(" | %10s | ThirdLight (J) | %10.3f |\n"),
fitted[Flags.simplex[33]], Orbit.Third[Jmag]);
fprintf(file_out,
_(" | %10s | ThirdLight (H) | %10.3f |\n"),
fitted[Flags.simplex[34]], Orbit.Third[Hmag]);
fprintf(file_out,
_(" | %10s | ThirdLight (K) | %10.3f |\n"),
fitted[Flags.simplex[35]], Orbit.Third[Kmag]);
fprintf(file_out,
_(" | %10s | ThirdLight (u) | %10.3f |\n"),
fitted[Flags.simplex[36]], Orbit.Third[umag]);
fprintf(file_out,
_(" | %10s | ThirdLight (v) | %10.3f |\n"),
fitted[Flags.simplex[37]], Orbit.Third[vmag]);
fprintf(file_out,
_(" | %10s | ThirdLight (b) | %10.3f |\n"),
fitted[Flags.simplex[38]], Orbit.Third[bmag]);
fprintf(file_out,
_(" | %10s | ThirdLight (y) | %10.3f |\n"),
fitted[Flags.simplex[39]], Orbit.Third[ymag]);
fprintf(file_out, " ----------------------------------------------------------\n");
if (Flags.Spots1 > 0) {
fprintf(file_out,
_(" | Spots on Primary |\n") );
fprintf(file_out,
_(" | Spot 1 |\n") );
fprintf(file_out,
_(" | %10s | Longitude | %10.3f |\n"),
fitted[Flags.simplex[10]], Spot[Primary][0].longitude );
fprintf(file_out,
_(" | %10s | Latitude | %10.3f |\n"),
fitted[Flags.simplex[11]], Spot[Primary][0].latitude );
fprintf(file_out,
_(" | %10s | Radius | %10.3f |\n"),
fitted[Flags.simplex[12]], Spot[Primary][0].radius );
fprintf(file_out,
_(" | %10s | Dimfactor | %10.3f |\n"),
fitted[Flags.simplex[13]], Spot[Primary][0].dimfactor );
if (Flags.Spots1 > 1) {
fprintf(file_out,
_(" | Spot 2 |\n") );
fprintf(file_out,
_(" | %10s | Longitude | %10.3f |\n"),
fitted[Flags.simplex[14]], Spot[Primary][1].longitude );
fprintf(file_out,
_(" | %10s | Latitude | %10.3f |\n"),
fitted[Flags.simplex[15]], Spot[Primary][1].latitude );
fprintf(file_out,
_(" | %10s | Radius | %10.3f |\n"),
fitted[Flags.simplex[16]], Spot[Primary][1].radius );
fprintf(file_out,
_(" | %10s | Dimfactor | %10.3f |\n"),
fitted[Flags.simplex[17]], Spot[Primary][1].dimfactor );
}
fprintf(file_out, " ----------------------------------------------------------\n");
}
if (Flags.Spots2 > 0) {
fprintf(file_out,
_(" | Spots on Secondary |\n") );
fprintf(file_out,
_(" | Spot 1 |\n") );
fprintf(file_out,
_(" | %10s | Longitude | %10.3f |\n"),
fitted[Flags.simplex[18]], Spot[Secondary][0].longitude );
fprintf(file_out,
_(" | %10s | Latitude | %10.3f |\n"),
fitted[Flags.simplex[19]], Spot[Secondary][0].latitude );
fprintf(file_out,
_(" | %10s | Radius | %10.3f |\n"),
fitted[Flags.simplex[20]], Spot[Secondary][0].radius );
fprintf(file_out,
_(" | %10s | Dimfactor | %10.3f |\n"),
fitted[Flags.simplex[21]], Spot[Secondary][0].dimfactor );
if (Flags.Spots2 > 1) {
fprintf(file_out,
_(" | Spot 2 |\n") );
fprintf(file_out,
_(" | %10s | Longitude | %10.3f |\n"),
fitted[Flags.simplex[22]], Spot[Secondary][1].longitude );
fprintf(file_out,
_(" | %10s | Latitude | %10.3f |\n"),
fitted[Flags.simplex[23]], Spot[Secondary][1].latitude );
fprintf(file_out,
_(" | %10s | Radius | %10.3f |\n"),
fitted[Flags.simplex[24]], Spot[Secondary][1].radius );
fprintf(file_out,
_(" | %10s | Dimfactor | %10.3f |\n"),
fitted[Flags.simplex[25]], Spot[Secondary][1].dimfactor );
}
fprintf(file_out, " ----------------------------------------------------------\n");
}
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Driver for print of fitted values
@param (double) Chi The merit function
@param (double) VarChiSquare Variance of merit function
@return (void)
@heading Data Fitting
****************************************************************************/
void SimplexPrintVariables(double Chi, double VarChiSquare)
{
FILE *file_out; /* the output file */
long i; /* loop counter */
int band; /* passband */
int param = 0; /* # of fitted parameters */
long dof = 0; /* degrees of freedom */
/* double Fillout = 0.0; fillout factor */
double ResMean[NUM_MAG+2]; /* mean residuals */
double ResDev[NUM_MAG+2]; /* sdv of residuals */
/* --------------- compute fillout --------------------------------- */
/*
if (Flags.fill == ON)
Fillout = (Binary[Primary].RochePot - Binary[Primary].RCLag1)
/ (Binary[Primary].RCLag2 - Binary[Primary].RCLag1);
*/
/* --------------- compute param --------------------------------- */
param = 0;
for (i = 0; i < SDIM; ++i) {
if (Flags.simplex[i] == ON) ++param;
}
/* --------------- compute mean residuals and sdv------------------- */
for (band=0; band < (NUM_MAG+2); ++band) {
if (Flags.Passbands[band] > 0) {
ResMean[band] = 0.0; ResDev[band] = 0.0;
for (i = 0; i < Flags.Passbands[band]; ++i) {
ResMean[band] =
ResMean[band] + Observed[band][i].residual;
}
ResMean[band] = ResMean[band]/Flags.Passbands[band];
for (i = 0; i < Flags.Passbands[band]; ++i) {
ResDev[band] = ResDev[band] +
SQR(ResMean[band]-Observed[band][i].residual);
}
ResDev[band] = sqrt(ResDev[band])/(Flags.Passbands[band]-1);
}
}
/* >>>>>>>>>>>>>>>>>>> terminal output <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON) {
OutfileHeader(stdout);
printf(_(" Fitted Data\n")); printf("\n");
dof = 0;
for (i = 0; i < NUM_MAG+2; ++i) {
if(Flags.Passbands[i] > 0)
printf(_(" %12s : %10ld data\n"),
Filters[i],
Flags.Passbands[i]);
dof = dof + Flags.Passbands[i];
}
dof = dof - param;
printf("\n");
printf(_(" ChiSquare = %10.5g\n"), Chi);
printf(_(" SDV (ChiSquare) = %10.5g\n"), VarChiSquare);
printf(_(" Degrees of Freedom = %10ld \n"), dof);
printf("\n");
printFitted(stdout);
for (band=0; band < (NUM_MAG+2); ++band) {
if (Flags.Passbands[band] > 0) {
fprintf(stdout,
_(" %12s: Mean Residual %7.4f SDV Residuals %7.4f\n"),
Filters[band], ResMean[band], ResDev[band]);
}
}
}
/* >>>>>>>>>>>>>>>>>>>>>>> file output <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
file_out = fopen(OUT_SIMPLEX, "w");
if (file_out == NULL)
#ifdef _WITH_GTK
{
if ( Flags.interactive == ON) {
make_dialog(_(errmsg[3]));
return;
} else nf_error(_(errmsg[3]));
}
#else
nf_error(_(errmsg[3]));
#endif
/* ----------------- write header -------------------------------------- */
OutfileHeader(file_out);
/* --------------- write fit results ----------------------------------- */
fprintf(file_out, "# \n");
fprintf(file_out, "# -------------------------------------------------\n");
fprintf(file_out, "# \n");
fprintf(file_out, "\n\n");
fprintf(file_out, _(" Fitted Data\n"));
fprintf(file_out, "\n");
dof = 0;
for (i = 0; i < NUM_MAG+2; ++i) {
if(Flags.Passbands[i] > 0)
fprintf(file_out, _(" %12s : %10ld data\n"),
Filters[i], Flags.Passbands[i]);
dof = dof + Flags.Passbands[i];
}
dof = dof - param;
fprintf(file_out, "\n");
fprintf(file_out, _(" ChiSquare = %10.5g\n"), Chi);
fprintf(file_out, _(" SDV (ChiSquare) = %10.5g\n"), VarChiSquare);
fprintf(file_out, _(" Degrees of Freedom = %10ld \n"), dof);
fprintf(file_out, "\n");
for (band=0; band < (NUM_MAG+2); ++band) {
if (Flags.Passbands[band] > 0) {
fprintf(file_out,
_(" %12s: Mean Residual %7.4f SDV Residuals %7.4f\n"),
Filters[band], ResMean[band], ResDev[band]);
}
}
printFitted(file_out);
/* ------------------------ close output file -------------------------- */
(void) fclose(file_out);
return;
}
/* ---------------------------------------------------------------------- */
/* */
/* */
/* random number generator */
/* */
/* */
/* ---------------------------------------------------------------------- */
/* Period parameters */
#define PER_N 624
#define PER_M 397
#define MATRIX_A 0x9908b0df /* constant vector a */
#define UPPER_MASK 0x80000000 /* most significant w-r bits */
#define LOWER_MASK 0x7fffffff /* least significant r bits */
/* for tempering */
#define TEMPERING_MASK_B 0x9d2c5680
#define TEMPERING_MASK_C 0xefc60000
#define TEMPERING_SHIFT_U(y) (y >> 11)
#define TEMPERING_SHIFT_S(y) (y << 7)
#define TEMPERING_SHIFT_T(y) (y << 15)
#define TEMPERING_SHIFT_L(y) (y >> 18)
static unsigned long ptgfsr[PER_N]; /* set initial seeds: PER_N = 624 words */
/****************************************************************************
@package nightfall
@author Makoto Matsumoto and Takuji Nishimura
@version MT19937B
@short Initialize 'Mersenne Twister' Random Number Generator
@long genrand() generate one pseudorandom number with double precision
which is uniformly distributed on [0,1]-interval for each call.
sgenrand(seed) set initial values to the working area of 624 words.
sgenrand(seed) must be called once before calling genrand()
(seed is any integer except 0).
Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
When you use this, send an email to: matumoto@math.keio.ac.jp
with an appropriate reference to your work.
Covered by GNU Public License
@param (unsigned long) seed The seed
@return (void)
@heading Data fitting
****************************************************************************/
void sgenrand(unsigned long seed)
{
/* seed should not be 0 */
int k;
/* setting initial seeds to ptgfsr[N] using */
/* the generator Line 25 of Table 1 in */
/* [KNUTH 1981, The Art of Computer Programming */
/* Vol. 2 (2nd Ed.), pp102] */
ptgfsr[0]= seed & 0xffffffff;
for (k=1; k<PER_N; k++)
ptgfsr[k] = (69069 * ptgfsr[k-1]) & 0xffffffff;
}
/****************************************************************************
@package nightfall
@author Makoto Matsumoto and Takuji Nishimura
@version MT19937B
@short 'Mersenne Twister' Random Number Generator (period 2^19937 - 1)
@long genrand() generate one pseudorandom number with double precision
which is uniformly distributed on [0,1]-interval for each call.
sgenrand(seed) set initial values to the working area of 624 words.
sgenrand(seed) must be called once before calling genrand()
(seed is any integer except 0).
Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
When you use this, send an email to: matumoto@math.keio.ac.jp
with an appropriate reference to your work.
Covered by GNU Public License
@param (void)
@return (double) A random number in [0,1]
@heading Data fitting
****************************************************************************/
double genrand()
{
unsigned long y;
static int k = 1;
static unsigned long mag01[2]={0x0, MATRIX_A};
/* mag01[x] = x * MATRIX_A for x=0,1 */
if(k == PER_N){ /* generate PER_N words at one time */
int kk;
for (kk=0;kk<PER_N-PER_M;kk++) {
y = (ptgfsr[kk]&UPPER_MASK)|(ptgfsr[kk+1]&LOWER_MASK);
ptgfsr[kk] = ptgfsr[kk+PER_M] ^ (y >> 1) ^ mag01[y & 0x1];
}
for (;kk<PER_N-1;kk++) {
y = (ptgfsr[kk]&UPPER_MASK)|(ptgfsr[kk+1]&LOWER_MASK);
ptgfsr[kk] = ptgfsr[kk+(PER_M-PER_N)] ^ (y >> 1) ^ mag01[y & 0x1];
}
y = (ptgfsr[PER_N-1]&UPPER_MASK)|(ptgfsr[0]&LOWER_MASK);
ptgfsr[PER_N-1] = ptgfsr[PER_M-1] ^ (y >> 1) ^ mag01[y & 0x1];
k = 0;
}
y = ptgfsr[k++];
y ^= TEMPERING_SHIFT_U(y);
y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
y &= 0xffffffff; /* you may delete this line if word size = 32 */
y ^= TEMPERING_SHIFT_L(y);
return ( (double)y / (unsigned long)0xffffffff );
}
/* ---------------------------------------------------------------------- */
/* */
/* */
/* simulated annealing routines */
/* */
/* */
/* ---------------------------------------------------------------------- */
static long states_uphill = 0;
static long states_accepted = 0;
static long states_valid = 0;
static long states_invalid = 0;
static long states_generated = 0;
static long states_repeat_cost = 0;
/**************************************************************************
@package nightfall
@author Rainer Wichmann
@version 1.0
@short Reset state of simulated annealing
@heading Data fitting
**************************************************************************/
void SaReset()
{
states_uphill = 0;
states_accepted = 0;
states_valid = 0;
states_invalid = 0;
states_generated = 0;
states_repeat_cost = 0;
}
/**************************************************************************
@package nightfall
@author Rainer Wichmann
@version 1.0
@short Equilibrate
@param (double) t Current temperature
@param (int) n Maximum iterations
@param (double) X[] Starting vertex
@param (double) *Y_Best Best merit so far
@param (double) BestEver[] Best ever vertex
@param (double) *BestChi Best ever merit
@return (int) Iterations done
@heading Data fitting
**************************************************************************/
int SaSweep(double t, int n, double X[], double *Y_Best,
double BestEver[], double *BestChi)
{
double delta; /* fractional change in energy */
double c; /* change in energy */
double p; /* acceptance probability */
double xc; /* delta(param) */
double yj; /* random variable */
double urand; /* random number */
double usign; /* signed random number */
double Y_New = 0.0; /* new merit */
double X_New[SDIM] = { 0.0 }; /* new vertex */
int VertexTest; /* exit status of vertex test */
int ErrCode; /* exit status of main loop */
int i, j, k; /* loop counters */
int equil = 0; /* equilibrium reached ? */
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON)
fprintf(stderr, _(" Enter Sweep, Temp = %8.4g\n"), t);
delta = 1.0;
for (i = 0; i < n; i++) {
do {
/* ------------ new vertex --------------------------------------- */
for (j = 0; j < SDIM; j++) {
if (Flags.simplex[j] == ON) {
urand = genrand();
if (urand >= 0.5) usign = (-1.0);
else usign = 1.0;
yj = usign * Param[j].t *
(pow((1.0 + 1.0/Param[j].t),fabs(2.0*urand - 1.0))-1.0);
xc = Param[j].Range * yj;
X_New[j] = X[j] + xc;
}
}
++states_generated;
/* ------------ validate new vertex ------------------------------ */
VertexTest = SimplexCheckVertex(X_New);
if (VertexTest > 0) {
Y_New = PENALTY;
++states_invalid;
} else {
SimplexSetVariables(X_New);
if (Flags.debug[OPTIMIZE] == ON) {
fprintf(stderr, _("State:") );
for (k = 0; k < SDIM; ++k) fprintf(stderr, " %7.4g", X_New[k]);
fprintf(stderr, "\n");
}
ErrCode = MainLoop(&Y_New);
if (ErrCode == 8) {
Y_New = PENALTY;
++states_invalid;
}
}
} while (Y_New >= PENALTY);
++states_valid;
/* ------------ accept/discard ------------------------------------- */
c = Y_New - *Y_Best;
if ( fabs(c) <= COST_PRECISION) ++states_repeat_cost;
else states_repeat_cost = 0;
if (c <= DBL_EPSILON) { /* energy reduced -> keep */
++states_accepted;
for (j = 0; j < SDIM; j++) {
X[j] = X_New[j];
}
*Y_Best = Y_New;
if (Y_New <= *BestChi) { /* save the best ever */
for (j = 0; j < SDIM; j++) {
BestEver[j] = X_New[j];
}
*BestChi = Y_New;
}
delta = fabs( c );
if ( fabs(Y_New) >= FLT_EPSILON ) delta = delta / Y_New;
else if ( fabs(*Y_Best ) >= FLT_EPSILON) delta = delta / *Y_Best;
else delta = delta / 1.0;
/* equilibrium is defined as a <= 10% change for 10 iterations */
if ( delta <= 0.10 ) ++equil;
else equil = 0;
} else {
p = exp( - (Y_New - *Y_Best) / (t) );
if ( p >= genrand() ) {
++states_accepted;
++states_uphill;
for (j = 0; j < SDIM; j++) {
X[j] = X_New[j];
}
*Y_Best = Y_New;
}
++equil;
}
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON)
printf(_(" Iter %d ChiSqr %10.3f Best %10.3f Equil %d\n"),
i+1, Y_New, *BestChi, equil);
if (equil > 9) break;
}
return (i + 1);
}
/**************************************************************************
@package nightfall
@author Rainer Wichmann
@version 1.0
@short Simulated annealing
@long Uses the 'Very Fast Annealing' algorithm
by L. Ingber, J. Math. Comput. Modelling 12, 967 (1989),
also available on http://www.ingber.com/
(no re-annealing, though)
@param (double) tscale Temperature scale
@param (double) t0 Starting temperature
@param (double) X[] Starting vertex
@param (double) *Y_Best Best merit so far
@param (double) BestEver[] Best ever vertex
@param (double) *BestChi Best ever merit
@return (float) Final temperature
@heading Data fitting
**************************************************************************/
/* simulated annealing */
float SaAnneal(double tscale, double t0, double X[], double *Y_Best,
double BestEver[], double *BestChi)
{
int i, j; /* loop counters */
double t; /* current temperature */
int dwell = 5; /* sweepsize */
int FitDim; /* Number of Fitted Parameters */
double tfactor; /* update factor for temperature */
sgenrand(4357);
FitDim = 0;
for (j = 0; j < SDIM; ++j) {
if (Flags.simplex[j] == ON) ++FitDim;
}
t = t0;
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON)
printf(_(" Start Annealing, Temp = %8.4g\n"), t);
(void) SaSweep( t0, 20 * dwell, X, Y_Best, BestEver, BestChi );
i = 0;
/* termination conditions */
while ( (states_generated < LIMIT_GENERATED )
&& (states_accepted < LIMIT_ACCEPTANCES )
&& (t >= TEMP_MINIMUM )
&& (states_repeat_cost < MAXIMUM_COST_REPEAT)) {
++i;
/* ------------ update temperature ---------------------------------- */
tfactor = exp( -tscale * pow ( (double) i, 1.0/(double)FitDim ) );
t = t0 * tfactor;
for (j = 0; j < SDIM; ++j)
Param[j].t = 1.0 * tfactor;
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON)
printf(_(" Annealing -- Temp = %8.4g Generated = %10ld Accepted = %10ld\n"),
t, states_generated, states_accepted);
(void) SaSweep(t, dwell, X, Y_Best, BestEver, BestChi);
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON) {
printf(_(" End Sweep, Best %10.4f --- Generated = %10ld Accepted = %10ld\n"),
*BestChi, states_generated, states_accepted );
printf("\n"); printf(_(" Best Vertex:\n") );
for (j = 0; j < SDIM; ++j) printf(" %7.4g", BestEver[j]);
printf("\n\n");
}
}
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON) {
printf("\n\n");
printf("**************************************************************\n");
printf("\n\n");
printf(_(" End Annealing -- Final Cost Temperature = %8.4g\n"), t);
printf("\n");
printf(_(" -- Generated %10ld states\n"),states_generated);
printf(_(" -- Rejected %10ld states\n"), states_invalid);
printf(_(" -- Validated %10ld states\n"), states_valid);
printf(_(" -- Accepted %10ld states\n"), states_accepted);
printf(_(" -- Uphill %10ld states\n"), states_uphill);
printf("\n");
if (states_generated > 0)
printf(_(" %7.3f per cent of generated states were valid\n"),
((float)states_valid/(float)states_generated) * 100.0 );
if (states_valid > 0)
printf(_(" %7.3f per cent of valid states were accepted\n"),
((float)states_accepted/(float)states_valid) * 100.0 );
if (states_accepted > 0)
printf(_(" %7.3f per cent of accepted states were uphill steps\n"),
((float)states_uphill/(float)states_accepted) * 100.0 );
printf("\n\n");
printf("**************************************************************\n");
printf("\n\n");
}
return ((float)t);
}
syntax highlighted by Code2HTML, v. 0.9.1