diff --git a/Simulation.sln b/Simulation.sln index 045a3b2c1dc..265863d4ce3 100644 --- a/Simulation.sln +++ b/Simulation.sln @@ -65,6 +65,16 @@ Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "QCIExe", "src\Simulation\Si EndProject Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "TargetedExe", "src\Simulation\Simulators.Tests\TestProjects\TargetedExe\TargetedExe.csproj", "{D292BF18-3956-4827-820E-254C3F81EF09}" EndProject +Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "src", "src", "{EBDE31D8-BB73-4E7E-B035-DE92657F3700}" +EndProject +Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Simulation", "Simulation", "{5497B844-8266-4B66-B51B-AE26148E9F78}" +EndProject +Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Simulators.Tests", "Simulators.Tests", "{C64D5562-CF69-4FF0-8A44-19FE8BEB8CE4}" +EndProject +Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "TestProjects", "TestProjects", "{877C3E74-5533-4517-8EB1-CA24EBAB4A25}" +EndProject +Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "IntrinsicTests", "src\Simulation\Simulators.Tests\TestProjects\IntrinsicTests\IntrinsicTests.csproj", "{D5D41201-101F-4C0D-B6A0-201D8FC3AB91}" +EndProject Global GlobalSection(SolutionConfigurationPlatforms) = preSolution Debug|Any CPU = Debug|Any CPU @@ -477,6 +487,22 @@ Global {D292BF18-3956-4827-820E-254C3F81EF09}.RelWithDebInfo|Any CPU.Build.0 = Release|Any CPU {D292BF18-3956-4827-820E-254C3F81EF09}.RelWithDebInfo|x64.ActiveCfg = Release|Any CPU {D292BF18-3956-4827-820E-254C3F81EF09}.RelWithDebInfo|x64.Build.0 = Release|Any CPU + {D5D41201-101F-4C0D-B6A0-201D8FC3AB91}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {D5D41201-101F-4C0D-B6A0-201D8FC3AB91}.Debug|Any CPU.Build.0 = Debug|Any CPU + {D5D41201-101F-4C0D-B6A0-201D8FC3AB91}.Debug|x64.ActiveCfg = Debug|Any CPU + {D5D41201-101F-4C0D-B6A0-201D8FC3AB91}.Debug|x64.Build.0 = Debug|Any CPU + {D5D41201-101F-4C0D-B6A0-201D8FC3AB91}.MinSizeRel|Any CPU.ActiveCfg = Debug|Any CPU + {D5D41201-101F-4C0D-B6A0-201D8FC3AB91}.MinSizeRel|Any CPU.Build.0 = Debug|Any CPU + {D5D41201-101F-4C0D-B6A0-201D8FC3AB91}.MinSizeRel|x64.ActiveCfg = Debug|Any CPU + {D5D41201-101F-4C0D-B6A0-201D8FC3AB91}.MinSizeRel|x64.Build.0 = Debug|Any CPU + {D5D41201-101F-4C0D-B6A0-201D8FC3AB91}.Release|Any CPU.ActiveCfg = Release|Any CPU + {D5D41201-101F-4C0D-B6A0-201D8FC3AB91}.Release|Any CPU.Build.0 = Release|Any CPU + {D5D41201-101F-4C0D-B6A0-201D8FC3AB91}.Release|x64.ActiveCfg = Release|Any CPU + {D5D41201-101F-4C0D-B6A0-201D8FC3AB91}.Release|x64.Build.0 = Release|Any CPU + {D5D41201-101F-4C0D-B6A0-201D8FC3AB91}.RelWithDebInfo|Any CPU.ActiveCfg = Debug|Any CPU + {D5D41201-101F-4C0D-B6A0-201D8FC3AB91}.RelWithDebInfo|Any CPU.Build.0 = Debug|Any CPU + {D5D41201-101F-4C0D-B6A0-201D8FC3AB91}.RelWithDebInfo|x64.ActiveCfg = Debug|Any CPU + {D5D41201-101F-4C0D-B6A0-201D8FC3AB91}.RelWithDebInfo|x64.Build.0 = Debug|Any CPU EndGlobalSection GlobalSection(SolutionProperties) = preSolution HideSolutionNode = FALSE @@ -508,6 +534,10 @@ Global {55833C6C-6E91-4413-9F77-96B3A09666B8} = {09C842CB-930C-4C7D-AD5F-E30DE4A55820} {C015FF41-9A51-4AF0-AEFC-2547D596B10A} = {09C842CB-930C-4C7D-AD5F-E30DE4A55820} {D292BF18-3956-4827-820E-254C3F81EF09} = {09C842CB-930C-4C7D-AD5F-E30DE4A55820} + {5497B844-8266-4B66-B51B-AE26148E9F78} = {EBDE31D8-BB73-4E7E-B035-DE92657F3700} + {C64D5562-CF69-4FF0-8A44-19FE8BEB8CE4} = {5497B844-8266-4B66-B51B-AE26148E9F78} + {877C3E74-5533-4517-8EB1-CA24EBAB4A25} = {C64D5562-CF69-4FF0-8A44-19FE8BEB8CE4} + {D5D41201-101F-4C0D-B6A0-201D8FC3AB91} = {877C3E74-5533-4517-8EB1-CA24EBAB4A25} EndGlobalSection GlobalSection(ExtensibilityGlobals) = postSolution SolutionGuid = {929C0464-86D8-4F70-8835-0A5EAF930821} diff --git a/src/Simulation/Common/Extensions.cs b/src/Simulation/Common/Extensions.cs index c3e3e00cb87..b36867a7388 100644 --- a/src/Simulation/Common/Extensions.cs +++ b/src/Simulation/Common/Extensions.cs @@ -49,5 +49,57 @@ where op.IsSubclassOf(typeof(T)) factory.Register(op.BaseType, op); } } + + internal static long NextLongBelow(this System.Random random, long upperExclusive) + { + // Don't allow sampling non-positive numbers, so that we don't break + // Math.Log2. + if (upperExclusive <= 0) + { + throw new ArgumentException( + $"Must be positive, got {upperExclusive}.", + nameof(upperExclusive) + ); + } + long SampleNBits(int nBits) + { + // Note that we can assume that nBytes is never more than 8, + // since we got there by looking at the bit length of + // upperExclusive, which is itself a 64-bit integer. + var nBytes = (nBits + 7) / 8; + var nExcessBits = nBytes * 8 - nBits; + var bytes = new byte[nBytes]; + random.NextBytes(bytes); + + // ToInt64 requires an array of exactly eight bytes. + // We can use IsLittleEndian to check which side we + // need to pad on to get there. + var padded = new byte[8]; + bytes.CopyTo(padded, + // ToInt64 requires an array of exactly eight bytes. + // We can use IsLittleEndian to check which side we + // need to pad on to get there. + System.BitConverter.IsLittleEndian + ? 0 + : 8 - nBytes + ); + return System.BitConverter.ToInt64(padded) >> nExcessBits; + }; + + var nBits = (int) (System.Math.Log(upperExclusive, 2) + 1); + var sample = SampleNBits(nBits); + while (sample >= upperExclusive) + { + sample = SampleNBits(nBits); + } + return sample; + } + + internal static long NextLong(this System.Random random, long lower, long upper) + { + var delta = upper - lower; + return lower + random.NextLongBelow(delta + 1); + } + } } diff --git a/src/Simulation/Common/SimulatorBase.cs b/src/Simulation/Common/SimulatorBase.cs index 7d216c00edf..bdb125468bf 100644 --- a/src/Simulation/Common/SimulatorBase.cs +++ b/src/Simulation/Common/SimulatorBase.cs @@ -31,6 +31,13 @@ public abstract class SimulatorBase : AbstractFactory, IOperat public event Action? OnLog = null; public event Action>? OnException = null; + + protected readonly int randomSeed; + protected readonly Lazy randomGenerator; + public int Seed => randomSeed; + protected System.Random RandomGenerator => randomGenerator.Value; + + /// /// An event fired whenever a simulator has additional diagnostic data @@ -55,8 +62,12 @@ public abstract class SimulatorBase : AbstractFactory, IOperat /// public StackFrame[]? CallStack { get; private set; } - public SimulatorBase(IQubitManager? qubitManager = null) + public SimulatorBase(IQubitManager? qubitManager = null, int? seed = null) { + this.randomSeed = seed ?? Guid.NewGuid().GetHashCode(); + this.randomGenerator = new Lazy( + () => new System.Random(Seed) + ); this.QubitManager = qubitManager; this.InitBuiltinOperations(this.GetType()); @@ -420,6 +431,55 @@ public GetQubitsAvailableToBorrow(SimulatorBase m) : base(m) sim.QubitManager.GetFreeQubitsCount(); } + /// + /// Implements the DrawRandomInt operation from the + /// Microsoft.Quantum.Random namespace. + /// + public class DrawRandomInt : Random.DrawRandomInt + { + private SimulatorBase sim; + + public DrawRandomInt(SimulatorBase m) : base(m) + { + sim = m; + } + + public override Func<(long, long), long> Body => arg => + { + var (min, max) = arg; + if (max <= min) + { + throw new ExecutionFailException($"Max must be greater than min, but {max} <= {min}."); + } + return sim.RandomGenerator.NextLong(min, max); + }; + } + + /// + /// Implements the DrawRandomInt operation from the + /// Microsoft.Quantum.Random namespace. + /// + public class DrawRandomDouble : Random.DrawRandomDouble + { + private SimulatorBase sim; + + public DrawRandomDouble(SimulatorBase m) : base(m) + { + sim = m; + } + + public override Func<(double, double), double> Body => arg => + { + var (min, max) = arg; + if (max <= min) + { + throw new ExecutionFailException($"Max must be greater than min, but {max} <= {min}."); + } + var delta = max - min; + return min + delta * sim.RandomGenerator.NextDouble(); + }; + } + public virtual void StartOperation(ICallable operation, IApplyData inputValue) { OnOperationStart?.Invoke(operation, inputValue); diff --git a/src/Simulation/QsharpCore/Diagnostics/Facts.qs b/src/Simulation/QsharpCore/Diagnostics/Facts.qs new file mode 100644 index 00000000000..846a65d22cf --- /dev/null +++ b/src/Simulation/QsharpCore/Diagnostics/Facts.qs @@ -0,0 +1,51 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +namespace Microsoft.Quantum.Diagnostics { + + /// # Summary + /// Declares that a classical condition is true. + /// + /// # Input + /// ## actual + /// The condition to be declared. + /// ## message + /// Failure message string to be printed in the case that the classical + /// condition is false. + /// + /// # See Also + /// - Microsoft.Quantum.Diagnostics.Contradiction + /// + /// # Example + /// The following Q# snippet will fail: + /// ```Q# + /// Fact(false, "Expected true."); + /// ``` + function Fact(actual : Bool, message : String) : Unit { + if (not actual) { fail message; } + } + + /// # Summary + /// Declares that a classical condition is false. + /// + /// # Input + /// ## actual + /// The condition to be declared. + /// ## message + /// Failure message string to be printed in the case that the classical + /// condition is true. + /// + /// # See Also + /// - Microsoft.Quantum.Diagnostics.Fact + /// + /// # Example + /// The following Q# code will print "Hello, world": + /// ```Q# + /// Contradiction(2 == 3, "2 is not equal to 3."); + /// Message("Hello, world."); + /// ``` + function Contradiction(actual : Bool, message : String) : Unit { + if (actual) { fail message; } + } + +} diff --git a/src/Simulation/QsharpCore/Diagnostics/Properties/NamespaceInfo.qs b/src/Simulation/QsharpCore/Diagnostics/Properties/NamespaceInfo.qs new file mode 100644 index 00000000000..60ad51909db --- /dev/null +++ b/src/Simulation/QsharpCore/Diagnostics/Properties/NamespaceInfo.qs @@ -0,0 +1,7 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +/// # Summary +/// This namespace contains functions and operations useful for diagnostic +/// purposes, including assert operations and claim functions. +namespace Microsoft.Quantum.Diagnostics {} diff --git a/src/Simulation/QsharpCore/Intrinsic.qs b/src/Simulation/QsharpCore/Intrinsic.qs index 92980459031..0a22cb388f1 100644 --- a/src/Simulation/QsharpCore/Intrinsic.qs +++ b/src/Simulation/QsharpCore/Intrinsic.qs @@ -5,27 +5,10 @@ namespace Microsoft.Quantum.Intrinsic { open Microsoft.Quantum.Math; open Microsoft.Quantum.Convert; - /// # Summary - /// The random operation takes an array of doubles as input, and returns - /// a randomly-selected index into the array as an `Int`. - /// The probability of selecting a specific index is proportional to the value - /// of the array element at that index. - /// Array elements that are equal to zero are ignored and their indices are never - /// returned. If any array element is less than zero, - /// or if no array element is greater than zero, then the operation fails. - /// - /// # Input - /// ## probs - /// An array of floating-point numbers proportional to the probability of - /// selecting each index. - /// - /// # Output - /// An integer $i$ with probability $\Pr(i) = p_i / \sum_i p_i$, where $p_i$ - /// is the $i$th element of `probs`. + @Deprecated("Microsoft.Quantum.Random.DrawCategorical") operation Random (probs : Double[]) : Int { body intrinsic; - } - + } @Deprecated("Microsoft.Quantum.Diagnostics.AssertMeasurement") operation Assert (bases : Pauli[], qubits : Qubit[], result : Result, msg : String) : Unit diff --git a/src/Simulation/QsharpCore/Random/Convienence.qs b/src/Simulation/QsharpCore/Random/Convienence.qs new file mode 100644 index 00000000000..aa35dc082c5 --- /dev/null +++ b/src/Simulation/QsharpCore/Random/Convienence.qs @@ -0,0 +1,115 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +namespace Microsoft.Quantum.Random { + open Microsoft.Quantum.Intrinsic; + open Microsoft.Quantum.Math; + open Microsoft.Quantum.Diagnostics; + + /// # Summary + /// Draws a random sample from a categorical distribution specified by a + /// list of probablities. + /// + /// # Description + /// The probability of selecting a specific index is proportional to the value + /// of the array element at that index. + /// Array elements that are equal to zero are ignored and their indices are never + /// returned. If any array element is less than zero, + /// or if no array element is greater than zero, then the operation fails. + /// + /// # Input + /// ## probs + /// An array of floating-point numbers proportional to the probability of + /// selecting each index. + /// + /// # Output + /// An integer $i$ with probability $\Pr(i) = p_i / \sum_i p_i$, where $p_i$ + /// is the $i$th element of `probs`. + /// + /// # See Also + /// - Microsoft.Quantum.Random.CategoricalDistribution + operation DrawCategorical(probs : Double[]) : Int { + // There are nicer ways of doing this, but they require the full + // standard library to be available. + mutable sum = 0.0; + for (prob in probs) { + Fact(prob >= 0.0, "Probabilities must be positive."); + set sum += prob; + } + + let variate = DrawRandomDouble(0.0, sum); + mutable acc = 0.0; + for (idx in 0..Length(probs) - 1) { + set acc += probs[idx]; + if (variate <= acc) { + return idx; + } + } + + return Length(probs) - 1; + } + + /// # Summary + /// Draws a random Pauli value. + /// + /// # Output + /// Either `PauliI`, `PauliX`, `PauliY`, or `PauliZ` with equal + /// probability. + operation DrawRandomPauli() : Pauli { + return [PauliI, PauliX, PauliY, PauliZ][DrawRandomInt(0, 3)]; + } + + /// # Summary + /// Given an array of data and an a distribution over its indices, + /// attempts to choose an element at random. + /// + /// # Input + /// ## data + /// The array from which an element should be chosen. + /// ## indexDistribution + /// A distribution over the indices of `data`. + /// + /// # Ouput + /// A tuple `(succeeded, element)` where `succeeded` is `true` if and only + /// if the sample chosen from `indexDistribution` was a valid index into + /// `data`, and where `element` is an element of `data` chosen at random. + /// + /// # Example + /// The following Q# snippet chooses an element at random from an array: + /// ```Q# + /// let (succeeded, element) = MaybeChooseElement( + /// data, + /// DiscreteUniformDistribution(0, Length(data) - 1) + /// ); + /// Fact(succeeded, "Index chosen by MaybeChooseElement was not valid."); + /// ``` + operation MaybeChooseElement<'T>(data : 'T[], indexDistribution : DiscreteDistribution) : (Bool, 'T) { + let index = indexDistribution::Sample(); + if (index >= 0 and index < Length(data)) { + return (true, data[index]); + } else { + return (false, Default<'T>()); + } + } + + /// # Summary + /// Given a success probability, returns a single Bernoulli trial that + /// is true with the given probability. + /// + /// # Input + /// ## successProbability + /// The probability with which `true` should be returned. + /// + /// # Output + /// `true` with probability `successProbability` and `false` with + /// probability `1.0 - successProbability`. + /// + /// # Example + /// The following Q# snippet samples flips from a biased coin: + /// ```Q# + /// let flips = DrawMany(DrawRandomBool, 10, 0.6); + /// ``` + operation DrawRandomBool(successProbability : Double) : Bool { + return DrawRandomDouble(0.0, 1.0) <= successProbability; + } +} diff --git a/src/Simulation/QsharpCore/Random/Internal.qs b/src/Simulation/QsharpCore/Random/Internal.qs new file mode 100644 index 00000000000..2dad23ac5e2 --- /dev/null +++ b/src/Simulation/QsharpCore/Random/Internal.qs @@ -0,0 +1,12 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +namespace Microsoft.Quantum.Random { + + // This operation duplicates an operation from the Q# standard libraries + // and is used to support partially applying all inputs to a given + // operation without actually executing it. + internal operation Delay<'TInput, 'TOutput>(op : ('TInput => 'TOutput), input : 'TInput, delay : Unit) : 'TOutput { + return op(input); + } +} diff --git a/src/Simulation/QsharpCore/Random/Intrinsic.qs b/src/Simulation/QsharpCore/Random/Intrinsic.qs new file mode 100644 index 00000000000..4ac4fe640ed --- /dev/null +++ b/src/Simulation/QsharpCore/Random/Intrinsic.qs @@ -0,0 +1,62 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +namespace Microsoft.Quantum.Random { + + /// # Summary + /// Draws a random integer in a given inclusive range. + /// + /// # Input + /// ## min + /// The smallest integer to be drawn. + /// ## max + /// The largest integer to be drawn. + /// + /// # Output + /// An integer in the inclusive range from `min` to `max` with uniform + /// probability. + /// + /// # Remarks + /// Fails if `max <= min`. + /// + /// # Example + /// The following Q# snippet randomly rolls a six-sided die: + /// ```Q# + /// let roll = DrawRandomInt(1, 6); + /// ``` + /// + /// # See Also + /// - Microsoft.Quantum.DiscreteUniformDistribution + operation DrawRandomInt(min : Int, max : Int) : Int { + body intrinsic; + } + + /// # Summary + /// Draws a random real number in a given inclusive interval. + /// + /// # Input + /// ## min + /// The smallest real number to be drawn. + /// ## max + /// The largest real number to be drawn. + /// + /// # Output + /// A random real number in the inclusive interval from `min` to `max` with + /// uniform probability. + /// + /// # Remarks + /// Fails if `max <= min`. + /// + /// # Example + /// The following Q# snippet randomly draws an angle between $0$ and $2 \pi$: + /// ```Q# + /// let angle = DrawRandomDouble(0.0, 2.0 * PI()); + /// ``` + /// + /// # See Also + /// - Microsoft.Quantum.ContinuousUniformDistribution + operation DrawRandomDouble(min : Double, max : Double) : Double { + body intrinsic; + } + +} diff --git a/src/Simulation/QsharpCore/Random/Normal.qs b/src/Simulation/QsharpCore/Random/Normal.qs new file mode 100644 index 00000000000..01f11b710f7 --- /dev/null +++ b/src/Simulation/QsharpCore/Random/Normal.qs @@ -0,0 +1,54 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +namespace Microsoft.Quantum.Random { + open Microsoft.Quantum.Intrinsic; + open Microsoft.Quantum.Math; + + internal operation DrawStandardNormalVariate() : Double { + let (u1, u2) = (DrawRandomDouble(0.0, 1.0), DrawRandomDouble(0.0, 1.0)); + return Sqrt(-2.0 * Log(u1)) * Cos(2.0 * PI() * u2); + } + + /// # Summary + /// Returns a normal distribution with mean 0 and variance 1. + /// + /// # Example + /// The following draws 10 samples from the standard normal distribution: + /// ```Q# + /// let samples = DrawMany((StandardNormalDistribution())::Sample, 10, ()); + /// ``` + /// + /// # See Also + /// - Microsoft.Quantum.Random.NormalDistribution + function StandardNormalDistribution() : ContinuousDistribution { + return ContinuousDistribution(DrawStandardNormalVariate); + } + + internal function StandardTransformation(mean : Double, variance : Double, variate : Double) : Double { + return mean + Sqrt(variance) * variate; + } + + /// # Summary + /// Returns a normal distribution with a given mean and variance. + /// + /// # Example + /// The following draws 10 samples from the normal distribution with mean + /// 2 and standard deviation 0.1: + /// ```Q# + /// let samples = DrawMany( + /// (NormalDistribution(2.0, PowD(0.1, 2.0)))::Sample, + /// 10, () + /// ); + /// ``` + /// + /// # See Also + /// - Microsoft.Quantum.Random.StandardNormalDistribution + function NormalDistribution(mean : Double, variance : Double) : ContinuousDistribution { + return TransformedContinuousDistribution( + StandardTransformation(mean, variance, _), + StandardNormalDistribution() + ); + } + +} diff --git a/src/Simulation/QsharpCore/Random/Types.qs b/src/Simulation/QsharpCore/Random/Types.qs new file mode 100644 index 00000000000..703701814d7 --- /dev/null +++ b/src/Simulation/QsharpCore/Random/Types.qs @@ -0,0 +1,123 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +namespace Microsoft.Quantum.Random { + open Microsoft.Quantum.Intrinsic; + open Microsoft.Quantum.Math; + + /// # Summary + /// Represents a univariate probability distribution over real numbers. + /// + /// # See Also + /// - Microsoft.Quantum.Random.ComplexDistribution + /// - Microsoft.Quantum.Random.DiscreteDistribution + /// - Microsoft.Quantum.Random.BigDiscreteDistribution + newtype ContinuousDistribution = ( + Sample : (Unit => Double) + ); + + + /// # Summary + /// Represents a univariate probability distribution over complex numbers. + /// + /// # See Also + /// - Microsoft.Quantum.Random.ContinuousDistribution + /// - Microsoft.Quantum.Random.DiscreteDistribution + /// - Microsoft.Quantum.Random.BigDiscreteDistribution + newtype ComplexDistribution = ( + Sample : (Unit => Complex) + ); + + /// # Summary + /// Represents a univariate probability distribution over integers. + /// + /// # See Also + /// - Microsoft.Quantum.Random.ContinuousDistribution + /// - Microsoft.Quantum.Random.ComplexDistribution + /// - Microsoft.Quantum.Random.BigDiscreteDistribution + newtype DiscreteDistribution = ( + Sample : (Unit => Int) + ); + + /// # Summary + /// Represents a univariate probability distribution over integers of + /// arbitrary size. + /// + /// # See Also + /// - Microsoft.Quantum.Random.ContinuousDistribution + /// - Microsoft.Quantum.Random.ComplexDistribution + /// - Microsoft.Quantum.Random.DiscreteDistribution + newtype BigDiscreteDistribution = ( + Sample : (Unit => BigInt) + ); + + /// # Summary + /// Returns a discrete categorical distribution, in which the probability + /// for each of a finite list of given outcomes is explicitly specified. + /// + /// # Input + /// ## probs + /// The probabilities for each outcome from the categorical distribution. + /// These probabilities may not be normalized, but must all be non-negative. + /// + /// # Output + /// The index `i` with probability `probs[i] / sum`, where `sum` is the sum + /// of `probs` given by `Fold(PlusD, 0.0, probs)`. + /// + /// # Example + /// The following Q# code will display 0 with probability 30% and 1 with + /// probability 70%: + /// ```Q# + /// let dist = CategoricalDistribution([0.3, 0.7]); + /// Message($"Got sample: {dist::Sample()}"); + /// ``` + function CategoricalDistribution(probs : Double[]) : DiscreteDistribution { + return DiscreteDistribution(Delay(DrawCategorical, probs, _)); + } + + /// # Summary + /// Internal-only operation for sampling from transformed distributions. + /// Should only be used via partial application. + internal operation SampleTransformedContinuousDistribution( + transform : (Double -> Double), + distribution : ContinuousDistribution + ) : Double { + return transform(distribution::Sample()); + } + + /// # Summary + /// Given a continuous distribution, returns a new distribution that + /// transforms the original by a given function. + /// + /// # Input + /// ## transform + /// A function that transforms variates of the original distribution to the + /// transformed distribution. + /// ## distribution + /// The original distribution to be transformed. + /// + /// # Output + /// A new distribution related to `distribution` by the transformation given + /// by `transform`. + /// + /// # Example + /// The following two distributions are identical: + /// ```Q# + /// let dist1 = ContinuousUniformDistribution(1.0, 2.0); + /// let dist2 = TransformedContinuousDistribution( + /// PlusD(1.0, _), + /// ContinuousUniformDistribution(0.0, 1.0) + /// ); + /// ``` + function TransformedContinuousDistribution( + transform : (Double -> Double), + distribution : ContinuousDistribution + ) : ContinuousDistribution { + return ContinuousDistribution(Delay( + SampleTransformedContinuousDistribution, + (transform, distribution), + _ + )); + } + +} diff --git a/src/Simulation/QsharpCore/Random/Uniform.qs b/src/Simulation/QsharpCore/Random/Uniform.qs new file mode 100644 index 00000000000..e022d9f6d46 --- /dev/null +++ b/src/Simulation/QsharpCore/Random/Uniform.qs @@ -0,0 +1,71 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +namespace Microsoft.Quantum.Random { + open Microsoft.Quantum.Intrinsic; + open Microsoft.Quantum.Math; + open Microsoft.Quantum.Diagnostics; + + /// # Summary + /// Returns a uniform distribution over a given inclusive interval. + /// + /// # Input + /// ## min + /// The smallest real number to be drawn. + /// ## max + /// The largest real number to be drawn. + /// + /// # Output + /// A distribution whose random variates are real numbers in the inclusive + /// interval from `min` to `max` with uniform probability. + /// + /// # Remarks + /// Fails if `max <= min`. + /// + /// # Example + /// The following Q# snippet randomly draws an angle between $0$ and $2 \pi$: + /// ```Q# + /// let angleDistribution = ContinuousUniformDistribution(0.0, 2.0 * PI()); + /// let angle = angleDistribution::Sample(); + /// ``` + /// + /// # See Also + /// - Microsoft.Quantum.DrawRandomDouble + function ContinuousUniformDistribution( + min : Double, max : Double + ) : ContinuousDistribution { + Fact(max > min, $"Max must be larger than min, but {max} <= {min}."); + return ContinuousDistribution(Delay(DrawRandomDouble, (min, max), _)); + } + + /// # Summary + /// Returns a uniform distribution over a given inclusive range. + /// + /// # Input + /// ## min + /// The smallest integer to be drawn. + /// ## max + /// The largest integer to be drawn. + /// + /// # Output + /// A distribution whose random variates are integers in the inclusive + /// range from `min` to `max` with uniform probability. + /// + /// # Remarks + /// Fails if `max <= min`. + /// + /// # Example + /// The following Q# snippet randomly rolls a six-sided die: + /// ```Q# + /// let dieDistribution = DiscreteUniformDistribution(1, 6); + /// let dieRoll = dieDistribution::Sample(); + /// ``` + /// + /// # See Also + /// - Microsoft.Quantum.DrawRandomDouble + function DiscreteUniformDistribution(min : Int, max : Int) : DiscreteDistribution { + Fact(max > min, $"Max must be larger than min, but {max} <= {min}."); + return DiscreteDistribution(Delay(DrawRandomInt, (min, max), _)); + } + +} diff --git a/src/Simulation/Simulators.Tests/OperationsTestHelper.cs b/src/Simulation/Simulators.Tests/OperationsTestHelper.cs index 37f0122a25c..7fa82ca1077 100644 --- a/src/Simulation/Simulators.Tests/OperationsTestHelper.cs +++ b/src/Simulation/Simulators.Tests/OperationsTestHelper.cs @@ -259,7 +259,7 @@ internal static void ctrlOnReleasedQubitTest(SimulatorBase sim, Action<(IQArray< /// internal static void ctrlOnReleasedCtrlQubitTest(SimulatorBase sim, Action<(IQArray, Qubit)> operationControlled) { - var random = new Random(); + var random = new System.Random(); ctrlTestShell(sim, operationControlled, (enabled, ctrlQs, q) => { diff --git a/src/Simulation/Simulators.Tests/QuantumSimulatorTests/VerifyGates.cs b/src/Simulation/Simulators.Tests/QuantumSimulatorTests/VerifyGates.cs index 59c7ae0e0d7..0cc32328fdf 100644 --- a/src/Simulation/Simulators.Tests/QuantumSimulatorTests/VerifyGates.cs +++ b/src/Simulation/Simulators.Tests/QuantumSimulatorTests/VerifyGates.cs @@ -31,7 +31,7 @@ public State((double, double) alpha, (double, double) beta) public partial class QuantumSimulatorTests { public const int seed = 19740212; - public static Random r = new Random(seed); + public static System.Random r = new System.Random(seed); public static double sqrt1_2 = Sqrt(1.0 / 2.0); diff --git a/src/Simulation/Simulators.Tests/TestProjects/IntrinsicTests/IntrinsicTests.csproj b/src/Simulation/Simulators.Tests/TestProjects/IntrinsicTests/IntrinsicTests.csproj new file mode 100644 index 00000000000..4fbf6b19d13 --- /dev/null +++ b/src/Simulation/Simulators.Tests/TestProjects/IntrinsicTests/IntrinsicTests.csproj @@ -0,0 +1,31 @@ + + + + netcoreapp3.1 + false + true + + false + false + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/Simulation/Simulators.Tests/TestProjects/IntrinsicTests/Random/Tests.qs b/src/Simulation/Simulators.Tests/TestProjects/IntrinsicTests/Random/Tests.qs new file mode 100644 index 00000000000..ab8416de9e6 --- /dev/null +++ b/src/Simulation/Simulators.Tests/TestProjects/IntrinsicTests/Random/Tests.qs @@ -0,0 +1,250 @@ +namespace Microsoft.Quantum.Tests { + open Microsoft.Quantum.Intrinsic; + open Microsoft.Quantum.Canon; + open Microsoft.Quantum.Random; + open Microsoft.Quantum.Diagnostics; + open Microsoft.Quantum.Convert; + open Microsoft.Quantum.Math; + + // Uses Welford's method to compute the mean and variance of an array + // of samples. + internal function SampleMeanAndVariance(samples : Double[]) : (Double, Double) { + mutable meanAcc = 0.0; + mutable varAcc = 0.0; + for (idx in 0..Length(samples) - 1) { + let sample = samples[idx]; + let oldMeanAcc = meanAcc; + let delta = (sample - meanAcc); + set meanAcc += delta / IntAsDouble(idx + 1); + set varAcc += delta * (sample - oldMeanAcc); + } + + return (meanAcc, varAcc / IntAsDouble(Length(samples) - 1)); + } + + internal operation EstimateMeanAndVariance(dist : ContinuousDistribution, nSamples : Int) : (Double, Double) { + mutable samples = new Double[nSamples]; + for (idx in 0..nSamples - 1) { + set samples w/= idx <- dist::Sample(); + } + return SampleMeanAndVariance(samples); + } + + internal operation CheckMeanAndVariance( + name : String, + distribution : ContinuousDistribution, + nSamples : Int, + (expectedMean : Double, expectedVariance : Double), + tolerance : Double + ) : Unit { + let (mean, variance) = EstimateMeanAndVariance( + distribution, + nSamples + ); + Fact( + expectedMean - tolerance <= mean and + mean <= expectedMean + tolerance, + $"Mean of {name} distribution should be {expectedMean}, was {mean}." + ); + Fact( + expectedVariance - tolerance <= variance and + variance <= expectedVariance + tolerance, + $"Variance of {name} distribution should be {expectedVariance}, was {variance}." + ); + } + + /// # Summary + /// Checks that @"microsoft.quantum.random.drawrandomdouble" obeys ranges. + @Test("QuantumSimulator") + @Test("ToffoliSimulator") + operation CheckDrawRandomDoubleObeysRanges() : Unit { + for (j in 0..10000) { + let random = DrawRandomDouble(0.0, 1.0); + if (random < 0.0 or random > 1.0) { + fail $"DrawRandomDouble(0.0, 1.0) returned {random}, outside the allowed interval."; + } + } + } + + /// # Summary + /// Checks that @"microsoft.quantum.random.drawrandomdint" obeys ranges. + @Test("QuantumSimulator") + @Test("ToffoliSimulator") + operation CheckDrawRandomIntObeysRanges() : Unit { + let randomInt = DrawRandomInt(0, 45); + if (randomInt > 45 or randomInt < 0) { + fail $"DrawRandomInt(0, 45) returned {randomInt}, outside the allowed range."; + } + } + + /// # Summary + /// Checks that @"microsoft.quantum.random.continuousuniformdistribution" has the + /// expected moments. + @Test("QuantumSimulator") + operation CheckContinuousUniformDistributionHasRightMoments() : Unit { + CheckMeanAndVariance( + "uniform", + ContinuousUniformDistribution(0.0, 1.0), + 1000000, + (0.5, 1.0 / 12.0), + 0.02 + ); + } + + /// # Summary + /// Checks that @"microsoft.quantum.random.standardnormaldistribution" has the + /// expected moments. + @Test("QuantumSimulator") + operation CheckStandardNormalDistributionHasRightMoments() : Unit { + CheckMeanAndVariance( + "standard normal", + StandardNormalDistribution(), + 1000000, + (0.0, 1.0), + 0.02 + ); + } + + /// # Summary + /// Checks that @"microsoft.quantum.random.normaldistribution" has the + /// expected moments. + @Test("QuantumSimulator") + operation CheckNormalDistributionHasRightMoments() : Unit { + CheckMeanAndVariance( + "normal(-2.0, 5.0)", + NormalDistribution(-2.0, 5.0), + 1000000, + (-2.0, 5.0), + 0.02 + ); + } + + /// # Summary + /// Checks that @"microsoft.quantum.random.drawrandombool" has the right + /// first moment. Note that since DrawRandomBool represents a Bernoulli + /// trial, it is entirely characterized by its first moment; we don't need + /// to check variance here. + @Test("QuantumSimulator") + operation CheckDrawRandomBoolHasRightExpectation() : Unit { + // NB: DrawMany isn't available yet, since it's in the + // Microsoft.Quantum.Standard package, not QSharpCore. + let prHeads = 0.65; + let nFlips = 1000000; + let stdDev = Sqrt(IntAsDouble(nFlips) * prHeads * (1.0 - prHeads)); + let expected = IntAsDouble(nFlips) * prHeads; + let nAllowedStdDev = 4.0; + mutable nHeads = 0; + for (idx in 0..nFlips - 1) { + if (DrawRandomBool(prHeads)) { + set nHeads += 1; + } + } + + let delta = IntAsDouble(nHeads) - expected; + + Fact( + -nAllowedStdDev * stdDev <= delta and + delta <= nAllowedStdDev * stdDev, + "First moment of Bernoulli distribution was incorrect." + ); + } + + /// # Summary + /// Checks that DrawCategorical never draws elements with probability zero. + @Test("QuantumSimulator") + operation CheckImpossibleEventsAreNotDrawn() : Unit { + let distribution = CategoricalDistribution([0.5, 0.0, 0.5]); + let nTrials = 100000; + for (idxTrial in 0..nTrials - 1) { + let variate = distribution::Sample(); + Fact( + variate != 1, + "A variate of 1 was drawn from a categorical distribution, despite having a probability of 0." + ); + } + } + + // We define a couple callables to help us run continuous tests on discrete + // distributions as well. + + internal operation DrawDiscreteAsContinuous(discrete : DiscreteDistribution, delay : Unit) : Double { + return IntAsDouble(discrete::Sample()); + } + + internal function DiscreteAsContinuous(discrete : DiscreteDistribution) : ContinuousDistribution { + return ContinuousDistribution(DrawDiscreteAsContinuous(discrete, _)); + } + + @Test("QuantumSimulator") + operation CheckCategoricalMomentsAreCorrect() : Unit { + let categorical = DiscreteAsContinuous( + CategoricalDistribution([0.2, 0.5, 0.3]) + ); + let expected = 0.0 * 0.2 + 1.0 * 0.5 + 2.0 * 0.3; + let variance = PowD(0.0 - expected, 2.0) * 0.2 + + PowD(1.0 - expected, 2.0) * 0.5 + + PowD(2.0 - expected, 2.0) * 0.3; + + CheckMeanAndVariance( + "categorical([0.2, 0.5, 0.3])", + categorical, + 1000000, + (expected, variance), + 0.04 + ); + } + + @Test("QuantumSimulator") + operation CheckRescaledCategoricalMomentsAreCorrect() : Unit { + let categorical = DiscreteAsContinuous( + CategoricalDistribution([2.0, 5.0, 3.0]) + ); + let expected = 0.0 * 0.2 + 1.0 * 0.5 + 2.0 * 0.3; + let variance = PowD(0.0 - expected, 2.0) * 0.2 + + PowD(1.0 - expected, 2.0) * 0.5 + + PowD(2.0 - expected, 2.0) * 0.3; + + CheckMeanAndVariance( + "categorical([0.2, 0.5, 0.3])", + categorical, + 1000000, + (expected, variance), + 0.04 + ); + } + + @Test("QuantumSimulator") + operation CheckCategoricalHistogramIsCorrect() : Unit { + let categorical = CategoricalDistribution([0.2, 0.5, 0.3]); + mutable counts = new Int[3]; + let nSamples = 1000000; + + for (idx in 0..nSamples - 1) { + let sample = categorical::Sample(); + set counts w/= sample <- counts[sample] + 1; + } + + Fact(190000 <= counts[0] and counts[0] <= 210000, $"counts[0] was {counts[0]}, expected about 200000."); + Fact(490000 <= counts[1] and counts[1] <= 510000, $"counts[1] was {counts[1]}, expected about 500000."); + Fact(290000 <= counts[2] and counts[2] <= 310000, $"counts[2] was {counts[2]}, expected about 300000."); + } + + @Test("QuantumSimulator") + operation CheckDiscreteUniformMomentsAreCorrect() : Unit { + let (min, max) = (-3, 7); + let expected = 0.5 * (IntAsDouble(min + max)); + let variance = (1.0 / 12.0) * ( + PowD(IntAsDouble(max - min + 1), 2.0) - 1.0 + ); + CheckMeanAndVariance( + $"discrete uniform ({min}, {max})", + DiscreteAsContinuous( + DiscreteUniformDistribution(min, max) + ), + 1000000, + (expected, variance), + 0.1 + ); + } + +} diff --git a/src/Simulation/Simulators/QCTraceSimulator/QCTraceSimulator.Primitive.random.cs b/src/Simulation/Simulators/QCTraceSimulator/QCTraceSimulator.Primitive.random.cs index 2a5578ed962..14523f44531 100644 --- a/src/Simulation/Simulators/QCTraceSimulator/QCTraceSimulator.Primitive.random.cs +++ b/src/Simulation/Simulators/QCTraceSimulator/QCTraceSimulator.Primitive.random.cs @@ -19,7 +19,7 @@ public TracerRandom(QCTraceSimulatorImpl m) : base(m) public override Func, Int64> Body => (p) => { - return CommonUtils.SampleDistribution(p, core.random.NextDouble()); + return CommonUtils.SampleDistribution(p, core.RandomGenerator.NextDouble()); }; } } diff --git a/src/Simulation/Simulators/QCTraceSimulator/QCTraceSimulatorImpl.cs b/src/Simulation/Simulators/QCTraceSimulator/QCTraceSimulatorImpl.cs index 53088993426..c3bd7f2de57 100644 --- a/src/Simulation/Simulators/QCTraceSimulator/QCTraceSimulatorImpl.cs +++ b/src/Simulation/Simulators/QCTraceSimulator/QCTraceSimulatorImpl.cs @@ -19,8 +19,6 @@ public partial class QCTraceSimulatorImpl : SimulatorBase protected readonly QCTraceSimulatorConfiguration configuration; private readonly QCTraceSimulatorCore tracingCore; private readonly double[] gateTimes; - - private readonly Random random = new Random(DateTime.Now.Millisecond); protected readonly QCTraceSimulatorCoreConfiguration tCoreConfig; private Dictionary primitiveOperationsIdToNames = diff --git a/src/Simulation/Simulators/QuantumProcessor/QuantumProcessorDispatcher.cs b/src/Simulation/Simulators/QuantumProcessor/QuantumProcessorDispatcher.cs index fa530a28316..c941ec7e133 100644 --- a/src/Simulation/Simulators/QuantumProcessor/QuantumProcessorDispatcher.cs +++ b/src/Simulation/Simulators/QuantumProcessor/QuantumProcessorDispatcher.cs @@ -12,10 +12,6 @@ namespace Microsoft.Quantum.Simulation.QuantumProcessor public partial class QuantumProcessorDispatcher : SimulatorBase { private const int PreallocatedQubitCount = 256; - /// - /// Random number generator used for Microsoft.Quantum.Intrinsic.Random - /// - public readonly System.Random random; public override string Name => "QuantumProcessorDispatcher"; @@ -35,9 +31,14 @@ public IQuantumProcessor QuantumProcessor /// An instance of a class implementing interface. If the parameter is null is used. /// A seed to be used by Q# Microsoft.Quantum.Intrinsic.Random operation. public QuantumProcessorDispatcher(IQuantumProcessor? quantumProcessor = null, IQubitManager? qubitManager = null, int? randomSeed = null) - : base(qubitManager ?? new QubitManagerTrackingScope(PreallocatedQubitCount, mayExtendCapacity:true, disableBorrowing:false)) + : base( + qubitManager ?? new QubitManagerTrackingScope( + PreallocatedQubitCount, + mayExtendCapacity: true, disableBorrowing: false + ), + randomSeed + ) { - random = new System.Random(randomSeed == null ? DateTime.Now.Millisecond : randomSeed.Value); QuantumProcessor = quantumProcessor ?? new QuantumProcessorBase(); OnOperationStart += QuantumProcessor.OnOperationStart; OnOperationEnd += QuantumProcessor.OnOperationEnd; diff --git a/src/Simulation/Simulators/QuantumProcessor/random.cs b/src/Simulation/Simulators/QuantumProcessor/random.cs index 0b0452c23c8..4badfd98802 100644 --- a/src/Simulation/Simulators/QuantumProcessor/random.cs +++ b/src/Simulation/Simulators/QuantumProcessor/random.cs @@ -20,7 +20,7 @@ public QuantumProcessorDispatcherRandom(QuantumProcessorDispatcher m) : base(m) public override Func, Int64> Body => (p) => { - return CommonUtils.SampleDistribution(p, Simulator.random.NextDouble()); + return CommonUtils.SampleDistribution(p, Simulator.RandomGenerator.NextDouble()); }; } } diff --git a/src/Simulation/Simulators/QuantumSimulator/QuantumSimulator.cs b/src/Simulation/Simulators/QuantumSimulator/QuantumSimulator.cs index e073a4ac39b..a4f7bf96ece 100644 --- a/src/Simulation/Simulators/QuantumSimulator/QuantumSimulator.cs +++ b/src/Simulation/Simulators/QuantumSimulator/QuantumSimulator.cs @@ -31,8 +31,6 @@ public partial class QuantumSimulator : SimulatorBase, IDisposable [DllImport(QSIM_DLL_NAME, ExactSpelling = true, CallingConvention = CallingConvention.Cdecl, EntryPoint = "seed")] private static extern void SetSeed(uint id, UInt32 seedValue); - public uint Seed{ get; private set; } - /// /// Creates a an instance of a quantum simulator. /// @@ -43,14 +41,15 @@ public QuantumSimulator( bool throwOnReleasingQubitsNotInZeroState = true, UInt32? randomNumberGeneratorSeed = null, bool disableBorrowing = false) - : base(new QSimQubitManager(throwOnReleasingQubitsNotInZeroState, disableBorrowing : disableBorrowing)) + : base( + new QSimQubitManager(throwOnReleasingQubitsNotInZeroState, disableBorrowing : disableBorrowing), + (int?)randomNumberGeneratorSeed + ) { - Seed = (randomNumberGeneratorSeed.HasValue) - ? randomNumberGeneratorSeed.Value - : (uint)Guid.NewGuid().GetHashCode(); - Id = Init(); - SetSeed(this.Id, this.Seed); + // Make sure that the same seed used by the built-in System.Random + // instance is also used by the native simulator itself. + SetSeed(this.Id, (uint)this.Seed); ((QSimQubitManager)QubitManager).Init(Id); } diff --git a/src/Simulation/Simulators/ToffoliSimulator/Random.cs b/src/Simulation/Simulators/ToffoliSimulator/Random.cs index 1d9d4a3520d..aa64984ec4e 100644 --- a/src/Simulation/Simulators/ToffoliSimulator/Random.cs +++ b/src/Simulation/Simulators/ToffoliSimulator/Random.cs @@ -3,6 +3,7 @@ using System; using System.Linq; +using Microsoft.Quantum.Simulation.Common; using Microsoft.Quantum.Simulation.Core; namespace Microsoft.Quantum.Simulation.Simulators @@ -12,9 +13,9 @@ public partial class ToffoliSimulator /// /// Implementation of the Random operation for the Toffoli simulator. /// - public class ToffSimRandom : Quantum.Intrinsic.Random + public class ToffSimRandom : Microsoft.Quantum.Intrinsic.Random { - Random random = new Random(); + private ToffoliSimulator Simulator; /// /// Constructs a new operation instance. @@ -22,36 +23,14 @@ public class ToffSimRandom : Quantum.Intrinsic.Random /// The simulator that this operation affects. public ToffSimRandom(ToffoliSimulator m) : base(m) { + Simulator = m; } /// /// The implementation of the operation. /// public override Func, Int64> Body => (probs) => - { - if (probs.Any(d => d < 0.0)) - { - throw new ArgumentOutOfRangeException("probs", "All probabilities must be greater than or equal to zero"); - } - var sum = probs.Sum(); - if (sum <= 0.0) - { - throw new ArgumentOutOfRangeException("probs", "At least one probability must be greater than zero"); - } - - var threshhold = random.NextDouble() * sum; - for (int i = 0; i < probs.Length; i++) - { - threshhold -= probs[i]; - if (threshhold <= 0.0) - { - return i; - } - } - - // This line should never be reached. - return probs.Length - 1; - }; + CommonUtils.SampleDistribution(probs, Simulator.RandomGenerator.NextDouble()); } } }