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.sampling;
019
020import org.apache.commons.rng.UniformRandomProvider;
021import org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler;
022import org.apache.commons.rng.sampling.distribution.ZigguratSampler;
023
024/**
025 * Generate vectors <a href="http://mathworld.wolfram.com/SpherePointPicking.html">
026 * isotropically located on the surface of a sphere</a>.
027 *
028 * <p>Sampling in 2 or more dimensions uses:</p>
029 *
030 * <ul>
031 *   <li>{@link UniformRandomProvider#nextLong()}
032 *   <li>{@link UniformRandomProvider#nextDouble()}
033 * </ul>
034 *
035 * <p>Sampling in 1D uses the sign bit from {@link UniformRandomProvider#nextInt()} to set
036 * the direction of the vector.
037 *
038 * @since 1.1
039 */
040public class UnitSphereSampler implements SharedStateObjectSampler<double[]> {
041    /** The dimension for 1D sampling. */
042    private static final int ONE_D = 1;
043    /** The dimension for 2D sampling. */
044    private static final int TWO_D = 2;
045    /** The dimension for 3D sampling. */
046    private static final int THREE_D = 3;
047    /**
048     * The mask to extract the second bit from an integer
049     * (naming starts at bit 1 for the least significant bit).
050     * The masked integer will have a value 0 or 2.
051     */
052    private static final int MASK_BIT_2 = 0x2;
053
054    /** The internal sampler optimised for the dimension. */
055    private final UnitSphereSampler delegate;
056
057    /**
058     * Sample uniformly from the ends of a 1D unit line.
059     */
060    private static final class UnitSphereSampler1D extends UnitSphereSampler {
061        /** The source of randomness. */
062        private final UniformRandomProvider rng;
063
064        /**
065         * @param rng Source of randomness.
066         */
067        UnitSphereSampler1D(UniformRandomProvider rng) {
068            this.rng = rng;
069        }
070
071        @Override
072        public double[] sample() {
073            // Either:
074            // 1 - 0 = 1
075            // 1 - 2 = -1
076            // Use the sign bit
077            return new double[] {1.0 - ((rng.nextInt() >>> 30) & MASK_BIT_2)};
078        }
079
080        @Override
081        public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) {
082            return new UnitSphereSampler1D(rng);
083        }
084    }
085
086    /**
087     * Sample uniformly from a 2D unit circle.
088     * This is a 2D specialisation of the UnitSphereSamplerND.
089     */
090    private static final class UnitSphereSampler2D extends UnitSphereSampler {
091        /** Sampler used for generating the individual components of the vectors. */
092        private final NormalizedGaussianSampler sampler;
093
094        /**
095         * @param rng Source of randomness.
096         */
097        UnitSphereSampler2D(UniformRandomProvider rng) {
098            sampler = ZigguratSampler.NormalizedGaussian.of(rng);
099        }
100
101        @Override
102        public double[] sample() {
103            final double x = sampler.sample();
104            final double y = sampler.sample();
105            final double sum = x * x + y * y;
106
107            if (sum == 0) {
108                // Zero-norm vector is discarded.
109                return sample();
110            }
111
112            final double f = 1.0 / Math.sqrt(sum);
113            return new double[] {x * f, y * f};
114        }
115
116        @Override
117        public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) {
118            return new UnitSphereSampler2D(rng);
119        }
120    }
121
122    /**
123     * Sample uniformly from a 3D unit sphere.
124     * This is a 3D specialisation of the UnitSphereSamplerND.
125     */
126    private static final class UnitSphereSampler3D extends UnitSphereSampler {
127        /** Sampler used for generating the individual components of the vectors. */
128        private final NormalizedGaussianSampler sampler;
129
130        /**
131         * @param rng Source of randomness.
132         */
133        UnitSphereSampler3D(UniformRandomProvider rng) {
134            sampler = ZigguratSampler.NormalizedGaussian.of(rng);
135        }
136
137        @Override
138        public double[] sample() {
139            final double x = sampler.sample();
140            final double y = sampler.sample();
141            final double z = sampler.sample();
142            final double sum = x * x + y * y + z * z;
143
144            if (sum == 0) {
145                // Zero-norm vector is discarded.
146                return sample();
147            }
148
149            final double f = 1.0 / Math.sqrt(sum);
150            return new double[] {x * f, y * f, z * f};
151        }
152
153        @Override
154        public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) {
155            return new UnitSphereSampler3D(rng);
156        }
157    }
158
159    /**
160     * Sample uniformly from a ND unit sphere.
161     */
162    private static final class UnitSphereSamplerND extends UnitSphereSampler {
163        /** Space dimension. */
164        private final int dimension;
165        /** Sampler used for generating the individual components of the vectors. */
166        private final NormalizedGaussianSampler sampler;
167
168        /**
169         * @param rng Source of randomness.
170         * @param dimension Space dimension.
171         */
172        UnitSphereSamplerND(UniformRandomProvider rng, int dimension) {
173            this.dimension = dimension;
174            sampler = ZigguratSampler.NormalizedGaussian.of(rng);
175        }
176
177        @Override
178        public double[] sample() {
179            final double[] v = new double[dimension];
180
181            // Pick a point by choosing a standard Gaussian for each element,
182            // and then normalize to unit length.
183            double sum = 0;
184            for (int i = 0; i < dimension; i++) {
185                final double x = sampler.sample();
186                v[i] = x;
187                sum += x * x;
188            }
189
190            if (sum == 0) {
191                // Zero-norm vector is discarded.
192                // Using recursion as it is highly unlikely to generate more
193                // than a few such vectors. It also protects against infinite
194                // loop (in case a buggy generator is used), by eventually
195                // raising a "StackOverflowError".
196                return sample();
197            }
198
199            final double f = 1 / Math.sqrt(sum);
200            for (int i = 0; i < dimension; i++) {
201                v[i] *= f;
202            }
203
204            return v;
205        }
206
207        @Override
208        public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) {
209            return new UnitSphereSamplerND(rng, dimension);
210        }
211    }
212
213    /**
214     * This instance delegates sampling. Use the factory method
215     * {@link #of(UniformRandomProvider, int)} to create an optimal sampler.
216     *
217     * @param dimension Space dimension.
218     * @param rng Generator for the individual components of the vectors.
219     * A shallow copy will be stored in this instance.
220     * @throws IllegalArgumentException If {@code dimension <= 0}
221     * @deprecated Use {@link #of(UniformRandomProvider, int)}.
222     */
223    @Deprecated
224    public UnitSphereSampler(int dimension,
225                             UniformRandomProvider rng) {
226        this(of(rng, dimension));
227    }
228
229    /**
230     * Private constructor used by deprecated constructor used to prevent partially
231     * initialized object if the construction of the delegate throws.
232     * In future versions the public constructor should be removed and the class made abstract.
233     *
234     * @param delegate Delegate.
235     */
236    private UnitSphereSampler(UnitSphereSampler delegate) {
237        this.delegate = delegate;
238    }
239
240    /**
241     * Private constructor used by sub-class specialisations.
242     * In future versions the public constructor should be removed and the class made abstract.
243     */
244    private UnitSphereSampler() {
245        delegate = null;
246    }
247
248    /**
249     * @return a random normalized Cartesian vector.
250     * @since 1.4
251     */
252    @Override
253    public double[] sample() {
254        return delegate.sample();
255    }
256
257    /**
258     * @return a random normalized Cartesian vector.
259     * @deprecated Use {@link #sample()}.
260     */
261    @Deprecated
262    public double[] nextVector() {
263        return sample();
264    }
265
266    /**
267     * {@inheritDoc}
268     *
269     * @since 1.3
270     */
271    @Override
272    public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) {
273        return delegate.withUniformRandomProvider(rng);
274    }
275
276    /**
277     * Create a unit sphere sampler for the given dimension.
278     *
279     * @param rng Generator for the individual components of the vectors. A shallow
280     * copy will be stored in the sampler.
281     * @param dimension Space dimension.
282     * @return the sampler
283     * @throws IllegalArgumentException If {@code dimension <= 0}
284     *
285     * @since 1.4
286     */
287    public static UnitSphereSampler of(UniformRandomProvider rng,
288                                       int dimension) {
289        if (dimension <= 0) {
290            throw new IllegalArgumentException("Dimension must be strictly positive");
291        } else if (dimension == ONE_D) {
292            return new UnitSphereSampler1D(rng);
293        } else if (dimension == TWO_D) {
294            return new UnitSphereSampler2D(rng);
295        } else if (dimension == THREE_D) {
296            return new UnitSphereSampler3D(rng);
297        }
298        return new UnitSphereSamplerND(rng, dimension);
299    }
300}