[Freeciv-Dev] (PR#11156) prune CM branches by fitness

["Jason Short" <jdorje@xxxxxxxxxxxxxxxxxxxxx>]

<URL: http://rt.freeciv.org/Ticket/Display.html?id=11156 >

I found a bug causing maximum fitness to be too large. Fixing that makes
things much faster (about 0.15s for "Boston").  However I think there
may be another bug that I'm forgetting to account for the city center in
init_min_production.  This would make things slower again.

jason

? diff
Index: common/aicore/cm.c
===================================================================
RCS file: /home/freeciv/CVS/freeciv/common/aicore/cm.c,v
retrieving revision 1.44
diff -u -r1.44 cm.c
--- common/aicore/cm.c  19 Nov 2004 02:31:35 -0000      1.44
+++ common/aicore/cm.c  23 Nov 2004 08:05:03 -0000
@@ -173,6 +173,7 @@
 struct cm_tile_type {
   int production[NUM_STATS];
   double estimated_fitness; /* weighted sum of production */
+  int max_fitness; /* maximum weighted sum of production */
   bool is_specialist;
   enum specialist_type spec; /* valid only if is_specialist */
   struct tile_vector tiles;  /* valid only if !is_specialist */
@@ -210,6 +211,7 @@
   /* the tile lattice */
   struct tile_type_vector lattice;
   struct tile_type_vector lattice_by_prod[NUM_STATS];
+  struct tile_type_vector lattice_by_max_fitness;
 
   /* the best known solution, and its fitness */
   struct partial_solution best;
@@ -221,6 +223,12 @@
    * fail (for being unhappy, for instance). */
   int min_production[NUM_STATS];
 
+  /* Minimum possible usage in each category (this should be an attempt
+   * at a greatest lower bound).  The fitness gives the minimum possible
+   * fitness "consumed" by the minimum production usage. */
+  int min_usage[NUM_STATS];
+  int min_usage_fitness;
+
   /* the current solution we're examining. */
   struct partial_solution current;
 
@@ -818,6 +826,35 @@
   return compare_tile_type_by_lattice_order(*a, *b);
 }
 
+/****************************************************************************
+  Sort by max fitness.  Since fitness is monotone in the production,
+  if a has higher fitness than b, then a cannot be a child of b, so
+  this respects the partial order -- unless a and b have equal fitness.
+  In that case, use compare_tile_type_by_lattice_order.
+****************************************************************************/
+static int compare_tile_type_by_max_fitness(const void *va, const void *vb)
+{
+  struct cm_tile_type * const *a = va;
+  struct cm_tile_type * const *b = vb;
+  double diff;
+
+  if (*a == *b) {
+    return 0;
+  }
+
+  /* To avoid double->int roundoff problems, we call a result non-zero only
+   * if it's larger than 0.5. */
+  diff = (*b)->max_fitness - (*a)->max_fitness;
+  if (diff > 0.5) {
+    return 1; /* return value is int; don't round down! */
+  }
+  if (diff < -0.5) {
+    return -1;
+  }
+
+  return compare_tile_type_by_lattice_order(*a, *b);
+}
+
 static enum cm_stat compare_key;
 
 /****************************************************************************
@@ -1102,6 +1139,10 @@
 ****************************************************************************/
 static double estimate_fitness(const struct cm_state *state,
                               const int production[NUM_STATS]);
+static int maximum_fitness(const struct cm_state *state,
+                          const int production[]);
+static int minimum_fitness(const struct cm_state *state,
+                          const int production[]);
 
 static void sort_lattice_by_fitness(const struct cm_state *state,
                                    struct tile_type_vector *lattice)
@@ -1111,6 +1152,9 @@
   /* compute fitness */
   tile_type_vector_iterate(lattice, ptype) {
     ptype->estimated_fitness = estimate_fitness(state, ptype->production);
+
+    /* We fill in maximum fitness here too. */
+    ptype->max_fitness = maximum_fitness(state, ptype->production);
   } tile_type_vector_iterate_end;
 
   /* sort by it */
@@ -1457,15 +1501,18 @@
   If we take the max along each production, and it's not better than the
   best in at least one stat, the partial solution isn't worth anything.
 
-  This function computes the max-stats produced by a partial solution.
+  This function computes the max-stats produced by a partial solution.  It
+  also calculates the maximum *fitness* a partial solution may provide in
+  a more sophisticated way.  This value (the maximum fitness) is returned.
 ****************************************************************************/
-static void compute_max_stats_heuristic(const struct cm_state *state,
-                                       const struct partial_solution *soln,
-                                       int production[NUM_STATS],
-                                       int check_choice)
+static int compute_max_stats_heuristic(const struct cm_state *state,
+                                      const struct partial_solution *soln,
+                                      int production[NUM_STATS],
+                                      int check_choice)
 {
   enum cm_stat stat;
   struct partial_solution solnplus; /* will be soln, plus some tiles */
+  int max_fitness, idle;
 
   /* Production is whatever the solution produces, plus the
      most possible of each kind of production the idle workers could
@@ -1481,7 +1528,7 @@
     for (stat = 0; stat < NUM_STATS; stat++) {
       production[stat] += ptype->production[stat];
     }
-    return;
+    return FC_INFINITY;
   }
 
   /* initialize solnplus here, after the shortcut check */
@@ -1498,6 +1545,29 @@
   }
 
   destroy_partial_solution(&solnplus);
+
+  /* Now compute the maximum fitness that this partial solution could ever
+   * provide.  This is an attempt at a least-upper-bound. */
+  max_fitness = -state->min_usage_fitness;
+  idle = soln->idle - 1;
+  max_fitness += tile_type_get(state, check_choice)->max_fitness;
+  tile_type_vector_iterate(&state->lattice_by_max_fitness, ptype) {
+    int index = ptype->lattice_index;
+    int used = soln->worker_counts[index];
+    int touse = MIN(soln->idle, tile_type_num_tiles(ptype) - used);
+
+    idle -= touse;
+    used += touse;
+
+    max_fitness += ptype->max_fitness * used;
+
+    if (idle == 0) {
+      break;
+    }
+  } tile_type_vector_iterate_end;
+  max_fitness += state->parameter.happy_factor;
+
+  return max_fitness;
 }
 
 /****************************************************************************
@@ -1508,11 +1578,20 @@
 ****************************************************************************/
 static bool choice_is_promising(struct cm_state *state, int newchoice)
 {
-  int production[NUM_STATS];
+  int production[NUM_STATS], fitness;
   enum cm_stat stat;
   bool beats_best = FALSE;
 
-  compute_max_stats_heuristic(state, &state->current, production, newchoice);
+  fitness = compute_max_stats_heuristic(state, &state->current, production,
+                                       newchoice);
+
+  if (fitness <= state->best_value.weighted) {
+    /* There's no possibility for this to beat the "best" solution. */
+    freelog(LOG_PRUNE_BRANCH,
+           "--- pruning: insufficient fitness (%d <= %d) [%d idle]",
+           fitness, state->best_value.weighted, state->current.idle);
+    return FALSE;
+  }
 
   for (stat = 0; stat < NUM_STATS; stat++) {
     if (production[stat] < state->min_production[stat]) {
@@ -1541,10 +1620,10 @@
 static void init_min_production(struct cm_state *state)
 {
   int x = CITY_MAP_RADIUS, y = CITY_MAP_RADIUS;
-  int usage[NUM_STATS];
   struct city *pcity = state->pcity;
   bool is_celebrating = base_city_celebrating(pcity);
   struct city backup;
+  int *usage = state->min_usage;
 
   /* make sure the city's numbers make sense (sometimes they don't,
    * somehow) */
@@ -1552,6 +1631,7 @@
   generic_city_refresh(pcity, FALSE, NULL);
 
   memset(state->min_production, 0, sizeof(state->min_production));
+  memset(state->min_usage, 0, sizeof(state->min_usage));
 
   /* If the city is content, then we know the food usage is just
    * prod-surplus; otherwise, we know it's at least 2*size but we
@@ -1595,6 +1675,16 @@
     state->min_production[SHIELD] = 0;
   }
 
+  /* HACK: there's no way to find gold usage other than this. */
+  usage[GOLD] = pcity->tax_total - city_gold_surplus(pcity, pcity->tax_total);
+
+  /* FIXME: are we forgetting about the city center here? */
+  state->min_usage_fitness = minimum_fitness(state, state->min_usage);
+
+  /* HACK: we can't prune from a minimum gold amount because gold isn't
+   * calculated until the end (most tiles produce trade). */
+  usage[GOLD] = 0;
+
   /* we should be able to get a min_production on gold and trade, too;
      also, lux, if require_happy, but not otherwise */
 
@@ -1642,7 +1732,118 @@
   return sum;
 }
 
+/****************************************************************************
+  Return the maximum fitness of the production fields of a tile.
+
+  It is very important that under no circumstances will this tile ever
+  provide more fitness than is returned by this function!  If that were to
+  happen then the CM algorithm will fail.
+****************************************************************************/
+static int maximum_fitness(const struct cm_state *state,
+                          const int production[])
+{
+  const struct city *pcity = state->pcity;
+  const struct player *pplayer = get_player(pcity->owner);
+  enum cm_stat stat;
+  double estimates[NUM_STATS];
+  int iestimates[NUM_STATS], sum;
+  const double f = 0.9999999999; /* 1.0 - epsilon */
+
+  for (stat = 0; stat < NUM_STATS; stat++) {
+    estimates[stat] = production[stat];
+  }
+
+  /* sci/lux/gold get benefit from the tax rates (in percentage).  We
+   * increase the value by 1 here to account for all the possible ways
+   * things could be rounded later (and just to be safe). */
+  estimates[SCIENCE] += (int)(pplayer->economic.science
+                             * estimates[TRADE] / 100.0 + f);
+  estimates[LUXURY] += (int)(pplayer->economic.luxury
+                            * estimates[TRADE] / 100.0 + f);
+  estimates[GOLD] += (int)(pplayer->economic.tax
+                          * estimates[TRADE] / 100.0 + f);
+
+  /* now add in the bonuses (none for food or trade) (in percentage) */
+  estimates[SHIELD] *= pcity->shield_bonus / 100.0;
+  estimates[LUXURY] *= pcity->luxury_bonus / 100.0;
+  estimates[GOLD] *= pcity->tax_bonus / 100.0;
+  estimates[SCIENCE] *= pcity->science_bonus / 100.0;
+
+  /* finally, sum it all up, rounding up to (again) account for possible
+   * later methods of rounding. */
+  sum = 0;
+  for (stat = 0; stat < NUM_STATS; stat++) {
+    iestimates[stat] = estimates[stat] + f;
+    sum += iestimates[stat] * state->parameter.factor[stat];
+  }
+
+#if 0
+  freelog(LOG_NORMAL, "Maximum fitness for %d,%d,%d,%d,%d,%d: %d",
+         production[0], production[1], production[2],
+         production[3], production[4], production[5], sum);
+  freelog(LOG_NORMAL, "  Estimates %f,%f,%f,%f,%f,%f",
+         estimates[0], estimates[1], estimates[2],
+         estimates[3], estimates[4], estimates[5]);
+  freelog(LOG_NORMAL, "  Iestimates %d,%d,%d,%d,%d,%d",
+         iestimates[0], iestimates[1], iestimates[2],
+         iestimates[3], iestimates[4], iestimates[5]);
+  freelog(LOG_NORMAL, "  Science bonus: %d", pcity->science_bonus);
+#endif
+  return sum;
+}
+
+/****************************************************************************
+  Return the maximum fitness of the production fields of a tile.
 
+  It is very important that under no circumstances will this tile ever
+  provide more fitness than is returned by this function!  If that were to
+  happen then the CM algorithm will fail.
+****************************************************************************/
+static int minimum_fitness(const struct cm_state *state,
+                          const int production[])
+{
+  const struct city *pcity = state->pcity;
+  const struct player *pplayer = get_player(pcity->owner);
+  enum cm_stat stat;
+  double estimates[NUM_STATS];
+  int iestimates[NUM_STATS], sum;
+
+  for (stat = 0; stat < NUM_STATS; stat++) {
+    estimates[stat] = production[stat];
+  }
+
+  /* sci/lux/gold get benefit from the tax rates (in percentage).  We
+   * decrease the value by 1 here to account for all the possible ways
+   * things could be rounded later (and just to be safe). */
+  estimates[SCIENCE] += (pplayer->economic.science
+                        * estimates[TRADE] / 100.0 - 1.0);
+  estimates[LUXURY] += (pplayer->economic.luxury
+                       * estimates[TRADE] / 100.0 - 1.0);
+  estimates[GOLD] += (pplayer->economic.tax
+                     * estimates[TRADE] / 100.0 - 1.0);
+
+  /* now add in the bonuses (none for food or trade) (in percentage) */
+  estimates[SHIELD] *= pcity->shield_bonus / 100.0;
+  estimates[LUXURY] *= pcity->luxury_bonus / 100.0;
+  estimates[GOLD] *= pcity->tax_bonus / 100.0;
+  estimates[SCIENCE] *= pcity->science_bonus / 100.0;
+
+  /* Now round up. */
+  for (stat = 0; stat < NUM_STATS; stat++) {
+    iestimates[stat] = estimates[stat] + 0.99999999;
+  }
+
+  /* finally, sum it all up, rounding down to (again) account for possible
+   * later methods of rounding. */
+  sum = 0;
+  for (stat = 0; stat < NUM_STATS; stat++) {
+    int iestimate;
+
+    iestimate = estimates[stat];
+    sum += iestimate * state->parameter.factor[stat];
+  }
+  return sum;
+}
 
 /****************************************************************************
   The high-level algorithm is:
@@ -1719,6 +1920,13 @@
          compare_tile_type_by_stat);
   }
 
+  tile_type_vector_init(&state->lattice_by_max_fitness);
+  tile_type_vector_copy(&state->lattice_by_max_fitness, &state->lattice);
+  qsort(state->lattice_by_max_fitness.p, state->lattice_by_max_fitness.size,
+       sizeof(*state->lattice_by_max_fitness.p),
+       compare_tile_type_by_max_fitness);
+
+
   /* We have no best solution yet, so its value is the worst possible. */
   init_partial_solution(&state->best, numtypes, pcity->size);
   state->best_value = worst_fitness();
@@ -1825,6 +2033,10 @@
 {
   struct cm_state *state = cm_init_state(pcity);
 
+  /* FIXME: For some reason the client may call this function with
+   * an unitialized city! */
+  generic_city_refresh(pcity, FALSE, 0);
+
   cm_find_best_solution(state, param, result);
   cm_free_state(state);
 }