source: branches/MCTS-SymbReg-2796/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/ExprHash.cs @ 15414

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
 * 
 * This file is part of HeuristicLab.
 *
 * HeuristicLab 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 3 of the License, or
 * (at your option) any later version.
 *
 * HeuristicLab 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 HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
 */
#endregion
using System;
using System.Collections.Generic;
using System.Diagnostics.Contracts;
using HeuristicLab.Random;
using System.Linq;

namespace HeuristicLab.Algorithms.DataAnalysis.MctsSymbolicRegression {

  // calculates a hash-code for expressions.
  // The hash codes of equivalent expressions x * (x + y) = (y + x)*x = (yx * xx) must be the same.
  // We use a random number for each symbol in the expressions and use symmetric operations to produce
  // results of sub-expressions.
  // Surrogate functions for the functions log, exp, 1/x are necessary.
  // Ideally, the surrogate functions also match identities of the original expressions.
  // Simply using floating point operations and random floats for symbols would be easiest. 
  // However, when we use floating point operations it is necessary to handle the inaccuracies 
  // in the last bits (e.g. for division, log and exp)
  // 
  // Collisions are OK as long as they occur seldom. 
  // It is probably be very unlikely to miss the optimal solution because of a collision. 
  //
  // We only need to identify equivalent structures. The numeric constants are irrelevant.
  // Therefore, equivalent structures with different numeric constants map to the same hash code.

  public static class ExprHash {
    const int MaxStackSize = 100;
    const int MaxVariables = 1000;
    private static double[] varSymbValues;
    static ExprHash() {
      var rand = new MersenneTwister();
      varSymbValues = new double[MaxVariables];
      for (int i = 0; i < MaxVariables; i++) {
        varSymbValues[i] = rand.NextDouble();
      }
    }

    public static ulong GetHash(byte[] code, int nParams) {
      var bits = (ulong)BitConverter.DoubleToInt64Bits(Eval(code, nParams));
      // clear last five bits (insignificant?)
      bits = bits & 0xFFFFFFFFFFFFFFE0;
      return bits;
    }

    private static double Eval(byte[] code, int nParams) {
      // The hash code calculation already preserves commutativity, associativity and distributivity of operations.
      // However, we also need to hash c1*x1 + c2*x1 to the same value as c3*x1.
      // Similarly for x1*x2 + x1*x2 or log(x1) + log(x1)!

      // Calculate sums lazily. Keep all terms and only when the actual sum is necessary remove duplicate terms and calculate sum
      // think about speed later (TODO)

      var stack = new ISet<double>[MaxStackSize];
      var terms = new HashSet<double>(new ApproximateDoubleEqualityComparer()); // the set of arguments for the current operator (+, *)
      int topOfStack = -1;
      int pc = 0;
      int nextParamIdx = -1;
      OpCodes op;
      short arg;
      while (true) {
        ReadNext(code, ref pc, out op, out arg);
        switch (op) {
          case OpCodes.Nop: break;
          case OpCodes.LoadConst0: {
              ++topOfStack;
              stack[topOfStack] = new HashSet<double>( new[] { 0.0 });

              // terms.Add(0.0); // ignore numeric constants in expr-hash

              break;
            }
          case OpCodes.LoadConst1: {
              ++topOfStack;
              stack[topOfStack] = new HashSet<double>(new[] { 1.0 });

              // args.Add(1.0); ignore numeric constants in expr-hash

              break;
            }
          case OpCodes.LoadParamN: {
              ++topOfStack;
              stack[topOfStack] =  new HashSet<double>(new[] { 1.0 });
              break;
            }
          case OpCodes.LoadVar: {
              ++topOfStack;
              stack[topOfStack] = new HashSet<double>(new[] { varSymbValues[arg] });

              // args.Add(varSymbValues[arg]);

              break;
            }
          case OpCodes.Add: {
              // take arguments from stack and put both terms into the set of terms
              // for every other operation we need to evaluate the sum of terms first and put it onto the stack (lazy eval of sums)

              stack[topOfStack - 1].UnionWith(stack[topOfStack]);
              topOfStack--;

              // stack[topOfStack] = t1 + t2; (later)
              break;
            }
          case OpCodes.Mul: {
              var t1 = stack[topOfStack - 1];
              var t2 = stack[topOfStack];
              topOfStack--;
              stack[topOfStack] = new HashSet<double>(new double[] { t1.Sum() * t2.Sum() });
              break;
            }
          case OpCodes.Log: {
              var v1 = stack[topOfStack];
              stack[topOfStack] = new HashSet<double>(new double[] { Math.Log( v1.Sum()) });
              break;
            }
          case OpCodes.Exp: {
              var v1 = stack[topOfStack];
              stack[topOfStack] = new HashSet<double>(new double[] { Math.Exp(v1.Sum()) });
              break;
            }
          case OpCodes.Inv: {
              var v1 = stack[topOfStack];
              stack[topOfStack] = new HashSet<double>(new double[] { 1.0 / v1.Sum() });
              break;
            }
          case OpCodes.Exit:
            Contract.Assert(topOfStack == 0);
            return stack[topOfStack].Sum();
        }
      }
    }

    private static void EvalTerms(HashSet<double> terms, double[] stack, ref int topOfStack) {
      ++topOfStack;
      stack[topOfStack] = terms.Sum();
      terms.Clear();
    }

    private static void ReadNext(byte[] code, ref int pc, out OpCodes op, out short s) {
      op = (OpCodes)Enum.ToObject(typeof(OpCodes), code[pc++]);
      s = 0;
      if (op == OpCodes.LoadVar) {
        s = (short)((code[pc] << 8) | code[pc + 1]);
        pc += 2;
      }
    }
  }
}