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

["Jason Short" <jdorje@xxxxxxxxxxxxxxxxxxxxx>]

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

This patch adds an additional pruning method for the CM algorithm.  It 
is pruned by fitness.

We can calculate the maximum fitness that can be provided by the 
production fields of a particular tile type.  If science is weighted 1 
and there is a +350% science bonus, then each scientist that produces 3 
base science will produce a maximum of 3 * 3.5 * 1 = 14 fitness. 
Meanwhile if luxury is weighted 0 each elvis will produce 0 fitness.

With a little more fiddling we can use this to calculate a maximum 
fitness that could be provided by any partial solution, and prune this 
branch if this maximum fitness is not at least as large as the best 
fitness.  Calculating the maximum fitness is a tricky and somewhat 
expensive process, involving rounding up at nearly every step.  It's 
surely bug-prone as well as being rather inextensible.  However with 
that change my Boston city that took 0.8 seconds under CM (include in 
PR#10203) now only takes 0.25 seconds.

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 07:35:52 -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,31 @@
   }
 
   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;
+  tile_type_vector_iterate(&state->lattice_by_max_fitness, ptype) {
+    int index = ptype->lattice_index;
+    int used = soln->worker_counts[index];
+
+    if (index == check_choice) {
+      used++;
+    }
+
+    if (idle > 0) {
+      int touse = MIN(soln->idle, tile_type_num_tiles(ptype) - used);
+
+      idle -= touse;
+      used += idle;
+    }
+
+    max_fitness += ptype->max_fitness * used;
+  } tile_type_vector_iterate_end;
+  max_fitness += state->parameter.happy_factor;
+
+  return max_fitness;
 }
 
 /****************************************************************************
@@ -1508,11 +1580,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 +1622,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 +1633,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 +1677,15 @@
     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);
+
+  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 +1733,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 +1921,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 +2034,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);
 }