source: branches/QAPGenerators/HeuristicLab.Problems.Instances.QAPGenerator/3.3/PRNG/BetaDistributedRandom.cs @ 14704

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2013 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 HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.Random;

namespace HeuristicLab.Problems.Instances.QAPGenerator {
  [Item("Beta-distributed random number generator", "Creates random numbers that are distributed according to the beta distribution.")]
  [StorableClass]
  public sealed class BetaDistributedRandom : Item, IRandom {
    private double alpha, beta, minimum, maximum;

    [Storable]
    public double Alpha {
      get { return alpha; }
      set { alpha = value; }
    }

    [Storable]
    public double Beta {
      get { return beta; }
      set { beta = value; }
    }

    [Storable]
    public double Minimum {
      get { return minimum; }
      set { minimum = value; }
    }

    [Storable]
    public double Maximum {
      get { return maximum; }
      set { maximum = value; }
    }

    [Storable]
    private IRandom uniform;

    [Storable]
    private IRandom normal;

    [StorableConstructor]
    private BetaDistributedRandom(bool deserializing) : base(deserializing) { }
    private BetaDistributedRandom(BetaDistributedRandom original, Cloner cloner)
      : base(original, cloner) {
      this.alpha = original.alpha;
      this.beta = original.beta;
      this.uniform = cloner.Clone(uniform);
    }
    public BetaDistributedRandom() : this(new MersenneTwister(), 5, 5) { }
    public BetaDistributedRandom(IRandom uniform, double alpha, double beta) : this(uniform, alpha, beta, 0, 1) { }
    public BetaDistributedRandom(IRandom uniform, double alpha, double beta, double minimum, double maximum) {
      if (alpha.IsAlmost(0) || beta.IsAlmost(0)) throw new ArgumentException("Alpha or Beta must be greater than 0.");
      this.uniform = uniform;
      this.normal = new NormalDistributedRandom(uniform, 0, 1);
      this.alpha = alpha;
      this.beta = beta;
      this.minimum = minimum;
      this.maximum = maximum;
    }

    public override IDeepCloneable Clone(Cloner cloner) {
      return new BetaDistributedRandom(this, cloner);
    }

    public void Reset() {
      uniform.Reset();
    }

    public void Reset(int seed) {
      uniform.Reset(seed);
    }

    public int Next() {
      throw new NotImplementedException();
    }

    public int Next(int maxVal) {
      throw new NotImplementedException();
    }

    public int Next(int minVal, int maxVal) {
      throw new NotImplementedException();
    }

    public double NextDouble() {
      if ((alpha <= 1.0) && (beta <= 1.0)) {
        // Use Jonk's algorithm
        while (true) {
          var x = Math.Pow(uniform.NextDouble(), 1.0 / alpha);
          var y = Math.Pow(uniform.NextDouble(), 1.0 / beta);
          if ((x + y) <= 1.0) return minimum + (maximum - minimum) * (x / (x + y));
        }
      }
      var Ga = NextDoubleGammaDistributed(alpha);
      var Gb = NextDoubleGammaDistributed(beta);
      return minimum + (maximum - minimum) * (Ga / (Ga + Gb));
    }

    private double NextDoubleGammaDistributed(double shape) {
      double u, v, x;

      if (shape.IsAlmost(1.0))
        return -Math.Log(1.0 - uniform.NextDouble());

      if (shape < 1.0) {
        while (true) {
          u = uniform.NextDouble();
          v = -Math.Log(1.0 - uniform.NextDouble());
          if (u <= 1.0 - shape) {
            x = Math.Pow(u, 1.0 / shape);
            if (x <= v) return x;
          } else {
            var y = -Math.Log((1 - u) / shape);
            x = Math.Pow(1.0 - shape + shape * y, 1.0 / shape);
            if (x <= (v + y)) return x;
          }
        }
      }
      // shape > 1.0
      var b = shape - 1.0 / 3.0;
      var c = 1.0 / Math.Sqrt(9 * b);
      for (; ; ) {
        do {
          x = normal.NextDouble();
          v = 1.0 + c * x;
        } while (v <= 0.0);

        v = v * v * v;
        u = uniform.NextDouble();
        if (u < 1.0 - 0.0331 * (x * x) * (x * x)) return (b * v);
        if (Math.Log(u) < 0.5 * x * x + b * (1.0 - v + Math.Log(v))) return (b * v);
      }
    }
  }
}