source: branches/HeuristicLab.Problems.GrammaticalOptimization-gkr/HeuristicLab.Problems.GrammaticalOptimization/Problems/PrimePolynomialProblem.cs @ 13042

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Text;
using System.Text.RegularExpressions;
using HeuristicLab.Common;
using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;

namespace HeuristicLab.Problems.GrammaticalOptimization {

  // Predicting Prime Numbers Using Cartesian Genetic Programming, James Alfred Walker and Julian Francis Miller
  // "Euler’s polynomial was the inspiration behind one of the GECCO competitions
  // in 2006. The aim of the GECCO Prime Prediction competition (and some
  // of the work in this paper) was to produce a polynomial f(i) with integer coeffi-
  // cients, such that given an integer input, i, it produces the i
  // th prime number, p(i),
  // for the largest possible value of i. For example, f(1) = 2, f(2) = 3, f(3) = 5,
  // f(4) = 7. Therefore, the function f(i) must produce consecutive prime numbers
  // for consecutive values of i. The requirement that the polynomial must not only
  // produce primes for consecutive input values, but also that the primes themselves
  // must be consecutive, makes the problem considerably more challenging
  // than mathematicians have previously considered. The two approaches described
  // in Section 4.2 were entered in the GECCO Prime Prediction competition and
  // were ranked second overall. The winning entry evolved floating point co-efficients
  // of a polynomial using a Genetic Algorithm (GA), where the output of the polynomial
  // was rounded to produce the prime numbers for consecutive values of i.
  // However, the winning entry was only able to predict correctly a few consecutive
  // prime numbers (9 in total). Unfortunately, the details regarding this have not
  // been published."

  public class PrimePolynomialProblem : ISymbolicExpressionTreeProblem {
    // C represents Koza-style ERC (the symbol is the index of the ERC), the values are initialized below
    private const string grammarString = @"
        G(E):
        E -> a | C | a+E | a-E | a*E | (E) | C+E | C-E | C*E  
        C -> 0..9
        ";

    // S .. Sum (+), N .. Neg. sum (-), P .. Product (*), D .. Division (%)
    private const string treeGrammarString = @"
        G(E):
        E -> a | C | S | N | P | D 
        S -> EE | EEE
        N -> EE | EEE
        P -> EE | EEE
        D -> EE
        C -> 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
        ";

    private readonly IGrammar grammar;
    private double[] erc;
    private double[] vars;
    // source http://primes.utm.edu/lists/small/1000.txt
    private double[] primes = new double[]
    {
      2,      3,      5,      7,     11,     13,     17,     19,     23,     29, 
     31,     37,     41,     43,     47,     53,     59,     61,     67,     71, 
     73,     79,     83,     89,     97,    101,    103,    107,    109,    113, 
    127,    131,    137,    139,    149,    151,    157,    163,    167,    173, 
    179,    181,    191,    193,    197,    199,    211,    223,    227,    229, 
    233,    239,    241,    251,    257,    263,    269,    271,    277,    281, 
    283,    293,    307,    311,    313,    317,    331,    337,    347,    349, 
    353,    359,    367,    373,    379,    383,    389,    397,    401,    409, 
    419,    421,    431,    433,    439,    443,    449,    457,    461,    463, 
    467,    479,    487,    491,    499,    503,    509,    521,    523,    541, 
    547,    557,    563,    569,    571,    577,    587,    593,    599,    601, 
    607,    613,    617,    619,    631,    641,    643,    647,    653,    659, 
    661,    673,    677,    683,    691,    701,    709,    719,    727,    733, 
    739,    743,    751,    757,    761,    769,    773,    787,    797,    809, 
    811,    821,    823,    827,    829,    839,    853,    857,    859,    863, 
    877,    881,    883,    887,    907,    911,    919,    929,    937,    941, 
    947,    953,    967,    971,    977,    983,    991,    997,   1009,   1013, 
   1019,   1021,   1031,   1033,   1039,   1049,   1051,   1061,   1063,   1069, 
   1087,   1091,   1093,   1097,   1103,   1109,   1117,   1123,   1129,   1151, 
   1153,   1163,   1171,   1181,   1187,   1193,   1201,   1213,   1217,   1223, 
   1229,   1231,   1237,   1249,   1259,   1277,   1279,   1283,   1289,   1291, 
   1297,   1301,   1303,   1307,   1319,   1321,   1327,   1361,   1367,   1373, 
   1381,   1399,   1409,   1423,   1427,   1429,   1433,   1439,   1447,   1451, 
   1453,   1459,   1471,   1481,   1483,   1487,   1489,   1493,   1499,   1511, 
   1523,   1531,   1543,   1549,   1553,   1559,   1567,   1571,   1579,   1583, 
   1597,   1601,   1607,   1609,   1613,   1619,   1621,   1627,   1637,   1657, 
   1663,   1667,   1669,   1693,   1697,   1699,   1709,   1721,   1723,   1733, 
   1741,   1747,   1753,   1759,   1777,   1783,   1787,   1789,   1801,   1811, 
   1823,   1831,   1847,   1861,   1867,   1871,   1873,   1877,   1879,   1889, 
   1901,   1907,   1913,   1931,   1933,   1949,   1951,   1973,   1979,   1987, 
   1993,   1997,   1999,   2003,   2011,   2017,   2027,   2029,   2039,   2053, 
   2063,   2069,   2081,   2083,   2087,   2089,   2099,   2111,   2113,   2129, 
   2131,   2137,   2141,   2143,   2153,   2161,   2179,   2203,   2207,   2213, 
   2221,   2237,   2239,   2243,   2251,   2267,   2269,   2273,   2281,   2287, 
   2293,   2297,   2309,   2311,   2333,   2339,   2341,   2347,   2351,   2357, 
   2371,   2377,   2381,   2383,   2389,   2393,   2399,   2411,   2417,   2423, 
   2437,   2441,   2447,   2459,   2467,   2473,   2477,   2503,   2521,   2531, 
   2539,   2543,   2549,   2551,   2557,   2579,   2591,   2593,   2609,   2617, 
   2621,   2633,   2647,   2657,   2659,   2663,   2671,   2677,   2683,   2687, 
   2689,   2693,   2699,   2707,   2711,   2713,   2719,   2729,   2731,   2741, 
   2749,   2753,   2767,   2777,   2789,   2791,   2797,   2801,   2803,   2819, 
   2833,   2837,   2843,   2851,   2857,   2861,   2879,   2887,   2897,   2903, 
   2909,   2917,   2927,   2939,   2953,   2957,   2963,   2969,   2971,   2999, 
   3001,   3011,   3019,   3023,   3037,   3041,   3049,   3061,   3067,   3079, 
   3083,   3089,   3109,   3119,   3121,   3137,   3163,   3167,   3169,   3181, 
   3187,   3191,   3203,   3209,   3217,   3221,   3229,   3251,   3253,   3257, 
   3259,   3271,   3299,   3301,   3307,   3313,   3319,   3323,   3329,   3331, 
   3343,   3347,   3359,   3361,   3371,   3373,   3389,   3391,   3407,   3413, 
   3433,   3449,   3457,   3461,   3463,   3467,   3469,   3491,   3499,   3511, 
   3517,   3527,   3529,   3533,   3539,   3541,   3547,   3557,   3559,   3571, 
   3581,   3583,   3593,   3607,   3613,   3617,   3623,   3631,   3637,   3643, 
   3659,   3671,   3673,   3677,   3691,   3697,   3701,   3709,   3719,   3727, 
   3733,   3739,   3761,   3767,   3769,   3779,   3793,   3797,   3803,   3821, 
   3823,   3833,   3847,   3851,   3853,   3863,   3877,   3881,   3889,   3907, 
   3911,   3917,   3919,   3923,   3929,   3931,   3943,   3947,   3967,   3989, 
   4001,   4003,   4007,   4013,   4019,   4021,   4027,   4049,   4051,   4057, 
   4073,   4079,   4091,   4093,   4099,   4111,   4127,   4129,   4133,   4139, 
   4153,   4157,   4159,   4177,   4201,   4211,   4217,   4219,   4229,   4231, 
   4241,   4243,   4253,   4259,   4261,   4271,   4273,   4283,   4289,   4297, 
   4327,   4337,   4339,   4349,   4357,   4363,   4373,   4391,   4397,   4409, 
   4421,   4423,   4441,   4447,   4451,   4457,   4463,   4481,   4483,   4493, 
   4507,   4513,   4517,   4519,   4523,   4547,   4549,   4561,   4567,   4583, 
   4591,   4597,   4603,   4621,   4637,   4639,   4643,   4649,   4651,   4657, 
   4663,   4673,   4679,   4691,   4703,   4721,   4723,   4729,   4733,   4751, 
   4759,   4783,   4787,   4789,   4793,   4799,   4801,   4813,   4817,   4831, 
   4861,   4871,   4877,   4889,   4903,   4909,   4919,   4931,   4933,   4937, 
   4943,   4951,   4957,   4967,   4969,   4973,   4987,   4993,   4999,   5003, 
   5009,   5011,   5021,   5023,   5039,   5051,   5059,   5077,   5081,   5087, 
   5099,   5101,   5107,   5113,   5119,   5147,   5153,   5167,   5171,   5179, 
   5189,   5197,   5209,   5227,   5231,   5233,   5237,   5261,   5273,   5279, 
   5281,   5297,   5303,   5309,   5323,   5333,   5347,   5351,   5381,   5387, 
   5393,   5399,   5407,   5413,   5417,   5419,   5431,   5437,   5441,   5443, 
   5449,   5471,   5477,   5479,   5483,   5501,   5503,   5507,   5519,   5521, 
   5527,   5531,   5557,   5563,   5569,   5573,   5581,   5591,   5623,   5639, 
   5641,   5647,   5651,   5653,   5657,   5659,   5669,   5683,   5689,   5693, 
   5701,   5711,   5717,   5737,   5741,   5743,   5749,   5779,   5783,   5791, 
   5801,   5807,   5813,   5821,   5827,   5839,   5843,   5849,   5851,   5857, 
   5861,   5867,   5869,   5879,   5881,   5897,   5903,   5923,   5927,   5939, 
   5953,   5981,   5987,   6007,   6011,   6029,   6037,   6043,   6047,   6053, 
   6067,   6073,   6079,   6089,   6091,   6101,   6113,   6121,   6131,   6133, 
   6143,   6151,   6163,   6173,   6197,   6199,   6203,   6211,   6217,   6221, 
   6229,   6247,   6257,   6263,   6269,   6271,   6277,   6287,   6299,   6301, 
   6311,   6317,   6323,   6329,   6337,   6343,   6353,   6359,   6361,   6367, 
   6373,   6379,   6389,   6397,   6421,   6427,   6449,   6451,   6469,   6473, 
   6481,   6491,   6521,   6529,   6547,   6551,   6553,   6563,   6569,   6571, 
   6577,   6581,   6599,   6607,   6619,   6637,   6653,   6659,   6661,   6673, 
   6679,   6689,   6691,   6701,   6703,   6709,   6719,   6733,   6737,   6761, 
   6763,   6779,   6781,   6791,   6793,   6803,   6823,   6827,   6829,   6833, 
   6841,   6857,   6863,   6869,   6871,   6883,   6899,   6907,   6911,   6917, 
   6947,   6949,   6959,   6961,   6967,   6971,   6977,   6983,   6991,   6997, 
   7001,   7013,   7019,   7027,   7039,   7043,   7057,   7069,   7079,   7103, 
   7109,   7121,   7127,   7129,   7151,   7159,   7177,   7187,   7193,   7207,
   7211,   7213,   7219,   7229,   7237,   7243,   7247,   7253,   7283,   7297, 
   7307,   7309,   7321,   7331,   7333,   7349,   7351,   7369,   7393,   7411, 
   7417,   7433,   7451,   7457,   7459,   7477,   7481,   7487,   7489,   7499, 
   7507,   7517,   7523,   7529,   7537,   7541,   7547,   7549,   7559,   7561, 
   7573,   7577,   7583,   7589,   7591,   7603,   7607,   7621,   7639,   7643, 
   7649,   7669,   7673,   7681,   7687,   7691,   7699,   7703,   7717,   7723, 
   7727,   7741,   7753,   7757,   7759,   7789,   7793,   7817,   7823,   7829, 
   7841,   7853,   7867,   7873,   7877,   7879,   7883,   7901,   7907,   7919, 
  };
    private bool[] isPrime = new bool[7920];
    public string Name { get { return "PrimePolynomial"; } }
    public PrimePolynomialProblem() {
      this.grammar = new Grammar(grammarString);
      this.TreeBasedGPGrammar = new Grammar(treeGrammarString);

      // values according to the paper: Monte-Carlo Expression Discovery, Cazenave, 2013
      erc = Enumerable.Range(1, 10).Select(x => (double)x).ToArray(); // the numbers 1..10;
      vars = new double[] { 0 }.Concat(primes.Take(25)).ToArray(); // the idx and all primes < 100;

      // for faster access;
      for (int i = 0; i < isPrime.Length; i++) {
        isPrime[i] = primes.Contains(i);
      }

    }

    public double BestKnownQuality(int maxLen) {
      // estimate
      return 100;
    }

    public IGrammar Grammar {
      get { return grammar; }
    }

    // according to above mentioned paper: ".. number of different primes it generations in a row for integer values of x starting at zero and increasing by one at each step"
    public double Evaluate(string sentence) {
      var interpreter = new ExpressionInterpreter();
      double[] vars = new[] { 0.0 };
      int i = 0;
      var primesFound = new HashSet<int>();
      for (; i <= primes.Length; i++) {
        vars[0] = i;
        var result = (int)interpreter.Interpret(sentence, vars, erc);
        if (result < 0 || result >= isPrime.Length || !isPrime[result]) {
          break;
        }
        primesFound.Add(result);
      }
      return primesFound.Count;
    }

    public string CanonicalRepresentation(string phrase) {
      return phrase;
    }

    public IEnumerable<Feature> GetFeatures(string phrase)
    {
      return new Feature[] {new Feature(phrase, 1.0),};
    }
    public IGrammar TreeBasedGPGrammar { get; private set; }
    public string ConvertTreeToSentence(ISymbolicExpressionTree tree) {
      var sb = new StringBuilder();

      TreeToSentence(tree.Root.GetSubtree(0).GetSubtree(0), sb);

      return sb.ToString();
    }
    public bool IsOptimalPhrase(string phrase) {
      // throw new NotImplementedException();
      return false; // TODO
    }


    private void TreeToSentence(ISymbolicExpressionTreeNode treeNode, StringBuilder sb) {
      if (treeNode.SubtreeCount == 0) {
        // terminal
        sb.Append(treeNode.Symbol.Name);
      } else {
        string op = string.Empty;
        switch (treeNode.Symbol.Name) {
          case "S": op = "+"; break;
          case "N": op = "-"; break;
          case "P": op = "*"; break;
          case "D": op = "%"; break;
          default: {
              Debug.Assert(treeNode.SubtreeCount == 1);
              break;
            }
        }
        // nonterminal
        if (treeNode.SubtreeCount > 1) sb.Append("(");
        TreeToSentence(treeNode.Subtrees.First(), sb);
        foreach (var subTree in treeNode.Subtrees.Skip(1)) {
          sb.Append(op);
          TreeToSentence(subTree, sb);
        }
        if (treeNode.SubtreeCount > 1) sb.Append(")");
      }
    }
  }
}