001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018package org.apache.commons.rng.examples.quadrature;
019
020import org.apache.commons.rng.UniformRandomProvider;
021
022/**
023 * <a href="https://en.wikipedia.org/wiki/Monte_Carlo_integration">Monte-Carlo method</a>
024 * for approximating an integral on a n-dimensional unit cube.
025 */
026public abstract class MonteCarloIntegration {
027    /** RNG. */
028    private final UniformRandomProvider rng;
029    /** Integration domain dimension. */
030    private final int dimension;
031
032    /**
033     * Simulation constructor.
034     *
035     * @param rng RNG.
036     * @param dimension Integration domain dimension.
037     */
038    public MonteCarloIntegration(UniformRandomProvider rng,
039                                 int dimension) {
040        this.rng = rng;
041        this.dimension = dimension;
042    }
043
044    /**
045     * Run the Monte-Carlo integration.
046     *
047     * @param n Number of random points to generate.
048     * @return the integral.
049     */
050    public double integrate(long n) {
051        long inside = 0;
052        for (long i = 0; i < n; i++) {
053            if (isInside(generateU01())) {
054                ++inside;
055            }
056        }
057
058        return inside / (double) n;
059    }
060
061    /**
062     * Indicates whether the given points is inside the region whose
063     * integral is computed.
064     *
065     * @param point Point whose coordinates are random numbers uniformly
066     * distributed in the unit interval.
067     * @return {@code true} if the {@code point} is inside.
068     */
069    protected abstract boolean isInside(double... point);
070
071    /**
072     * @return a value from a random sequence uniformly distributed
073     * in the {@code [0, 1)} interval.
074     */
075    private double[] generateU01() {
076        final double[] rand = new double[dimension];
077
078        for (int i = 0; i < dimension; i++) {
079            rand[i] = rng.nextDouble();
080        }
081
082        return rand;
083    }
084}